An Inversion-Free Method for Finding Positive Definite Solution of a Rational Matrix Equation

A new iterative scheme has been constructed for finding minimal solution of a rational matrix equation of the form X + A*X −1 A = I. The new method is inversion-free per computing step. The convergence of the method has been studied and tested via numerical experiments.


Introduction
In this paper, we will discuss the following nonlinear matrix equation: where is an × nonsingular complex matrix, is the unit matrix of the appropriate size, and ∈ C × is an unknown Hermitian positive definite (HPD) matrix that should be found. It was proved in [1] that if (1) has an HPD solution, then all its Hermitian solutions are positive definite and, moreover, it has the maximal solution and the minimal solution in the sense that ≤ ≤ for any HPD solution .
A lot of papers have been published regarding the iterative HPD solutions of such nonlinear rational matrix equations in the literature due to their importance in some practical problems arising in control theory, dynamical problems, and so forth (see [2,3]).
The most common iterative method for finding the maximal solution of (1) is the following fixed-point iteration [4]: The maximal solution of (1) can be obtained through = − , where is the minimal solution of the dual equation + −1 * = . In 2010, Monsalve and Raydan in [5] proposed the following iteration method (also known as Newton's method) for finding the minimal solution: which is an inversion-free scheme. Since, −1 should be computed only once in contrast to the matrix iteration (2). Note that − * = −1 * , and similar notations are used throughout.
Remark 1. We remark that there are several other well-known iterative methods for solving (1) rather than Newton's method (3). To the best of our knowledge, the procedure of extending higher-order iterative methods for finding the solution of (1) has not been exploited up to now. Hence, we hope that this interlink among the fields of root-finding and solving (1) may lead to discovering novel and innovative techniques. The Scientific World Journal The rest of this paper is organized as follows. In Section 2, we develop and analyze a new inversion-free method for finding roots of a special map . In Section 3, we provide some numerical comparisons by employing some experiments in machine precision. Some concluding remarks will be drawn in Section 4.

A New Iterative Method
An equivalent formulation of (1) is to find an HPD matrix such that ( ) = 0. Toward this goal, we write ( ) := + * −1 − = 0. Furthermore, we have Now, using a change of variable as = −1 − * , we could simplify (4) as follows: In order to obtain an iterative method for finding the minimal HPD solution of (1), it is now enough to solve the well-known matrix equation −1 − = 0, at which = ( −1 − * − ) * . One of such ways to challenge this matrix inversion problem is via applying the Schulz-type iteration methods (see, e.g., [6][7][8][9][10]). Applying Chebyshev's method [11] yields We remark that there is a tight relationship between iterative methods for nonlinear systems and the construction of higher-order methods for matrix equations ( [12,13]).
The matrix iteration (6) requires −1 to be computed only once at the beginning of the iteration and this makes the iterative method fall in the category of inversion-free algorithms for solving (1).
In the meantime, it is easy to show that the zeros of the map ( ) = −1 − are equal to the zeros of the map To be more precise, we are finding the inverse of the matrix which matches the minimal HPD solution of (1).

Remark 2.
Following Remark 1, we applied Chebyshev's method for (1) in this work and will study its theoretical behavior. The extension of the other well-known root-finding schemes for finding the minimal solution of (1) will remain for future studies. Proof. The initial matrix * is Hermitian, and = − * ( − ) −1 . Thus, 0 = − * −1 − − * * −1 is also Hermitian; that is, * 0 = 0 . Now using inductive argument, we have By considering ( ) * = , ( ≥ ) we now show that Note that = ( ) * has been used in (8). Now the conclusion holds for any +1. Thus, the proof is complete.

Theorem 4. By considering that and
are nonsingular matrices, the sequence { } generated by (6) is convergent to the minimal solution using the initial matrix 0 = * .
Proof. Let us consider = − * ( − ) −1 . We therefore have Taking a generic matrix operator norm from both sides of (9), we obtain On the other hand, Chebyshev's method for matrix inversion problem is convergent if the initial approximation reads ‖ − 0 ‖ < 1. That is to say, ‖ −[ − * ( − ) −1 ] 0 ‖ < 1. This together with the initial matrix 0 = * gives − * * < 1, which is true when is the minimal HPD solution of (1). Note that since The Scientific World Journal 3 we obtain that −1 0 > −1 ; thus 0 < . And subsequently using mathematical induction, it would be observed that { } tends to .
The only problem that happens in this process is the fact that the convergence order is -linear. In fact, although Chebyshev's method for matrix inversion has third local order of convergence, this rate will not be preserved for finding the minimal HPD solution of (1).
The reason is that the matrix , which we must compute its inverse by Chebyshev's method, is dependent on the itself. That is to say, the unknown is located in the essence of the matrix = − * ( − ) −1 .

Theorem 5.
The sequence of matrices produced by (6) satisfies the following error inequality: Proof. First since lim → ∞ = and by using (6) we have Note that we have used the fact that ( − ) = * −1 . Relation (14) yields Consequently, one has the error inequality (13). This shows the -linear order of convergence for finding the minimal HPD solution of (1). We thus have which is guaranteed since

Numerical Comparisons
In this section, we mainly investigate the performance of the new method (6) for matrix equation (1). All experiments were run on a Pentium IV computer, using Mathematica 8 [14]. We report the number of required iterations (Iter) for converging. In our implementations, we stop all considered methods when the infinity norm of two successive iterates is less than given tolerance. Note that recently Zhang in [15] studied a way to accelerate the beginning of such iterative methods for finding the minimal solution of (1) via applying multiple Newton's method for matrix inversion. This technique could be given by for any 1 ≤ ≤ 2. Subsequently, we could improve the behavior of the new method (6) using (18) as provided in Algorithm 1.
We compare Algorithm 1, denoted by PM, with (2) denoted by M1, (3) denoted by M2, and the method proposed by El-Sayed and Al-Dbiban [16] denoted by M3, which is a modification of the method presented by Zhan in [17], as follows: 4 The Scientific World Journal Example 1 (see [18]). In this experiment, we compare the results of different methods for finding the minimal solution of (1) when the matrix is defined by The results are given in Figure 1(a) in terms of the number of iterations when the stopping criterion is ‖ +1 − ‖ ∞ ≤ 10 −8 .
Note that PM and M2 converge to , whereas all the other schemes converge to ; thus we use other schemes to find the maximal solution of the dual equation ( ) = + −1 * − in our written codes so as to have fair comparisons. Here = 2 has been chosen for PM (with = 19). This = 19 for the number of iterations in the inner finite loop of Algorithm 1 has been considered in the numerical report.
Furthermore, we have chosen this number empirically. In fact, varying shows us that we even can obtain better or worse results than the reported ones in different examples.
In Example 1, we have used = 2. In fact we have chosen this value for since we are solving an operator equation in essence. To be more precise, we wish to consider the solution of the operator equation to be of multiplicity 2. This consideration makes the algorithm converge faster at the initial phase of the process and when we are enough close to the solution, then we flash back to the ordinary methods, that is, treat the solution as a simple zero (solution) of the operator equation.
Example 2 (see [19]). Applying the stopping criterion ‖ +1 − ‖ ∞ ≤ 10 −12 , we compare the behavior of various methods for the following test matrix: The results are illustrated in Figure 1(b), wherein = 1.2 has been chosen for PM (with = 1).
The Scientific World Journal 5

Conclusions
We have studied the fact that the minimal HPD solution of (1) is equivalent to the roots of a nonlinear map. This special map has been solved by the well-known Chebyshev method as a matrix inversion problem. The developed method requires the computation of one matrix inverse at the beginning of the process and it is hence an inversion-free method. The convergence and the rate of convergence have been studied for this scheme. Furthermore, using a proper acceleration technique from the literature, we have further speeded up the process of finding the HPD solution of (1).