Inversion Free Algorithms for Computing the Principal Square Root of a Matrix

New algorithms are presented about the principal square root of an n×nmatrixA. In particular, all the classical iterative algorithms require matrix inversion at every iteration. The proposed inversion free iterative algorithms are based on the Schulz iteration or the Bernoulli substitution as a special case of the continuous time Riccati equation. It is certified that the proposed algorithms are equivalent to the classical Newton method. An inversion free algebraic method, which is based on applying the Bernoulli substitution to a special case of the continuous time Riccati equation, is also proposed.


Introduction
Let  be a real or complex × matrix with no eigenvalues on R − (the closed negative real axis).Then there exists a unique  ×  matrix  such that  2 =  and the eigenvalues of  lie in the segment { : −/2 < arg() < /2}.We refer to  as the principle square root of .
The computation of the principal square root of a matrix is a problem of great interest and practical importance with numerous applications in mathematical and engineering problems.Due to the importance of the problem, many iterative algorithms have been proposed and successfully employed for calculating the principal square root of a matrix [1][2][3][4][5][6][7] without seeking the eigenvalues and eigenvectors of the matrix; these algorithms require matrix inversion at every iteration.Blocked Schur Algorithms for Computing the Matrix Square Root are proposed in [8] where the matrix is reduced to upper triangular form and a recurrence relation enables the square root of the triangular matrix to be computed a column or superdiagonal at a time.In this paper new inversion free iterative algorithms for computing the principal square root of a matrix are proposed.
The paper is organized as follows: in Section 2 a survey of classical iterative algorithms is presented; all the algorithms use matrix inversion in every iteration.In Section 3, the inversion free iterative algorithms are developed; the algorithms are derived by the Schulz iteration for computing the inversion of a matrix or by eliminating matrix inversion from the convergence criteria of the algorithms.In Section 4, an algebraic method for computing the principal square root of a matrix is presented; the method is associated with the solution of a special case of the continuous time Riccati equation, which arises in linear estimation [9] and requires only one matrix inversion (this is not an iterative method).Simulation results are presented in Section 5. Section 6 summarizes the conclusions.

Inversion Use Iterative Algorithms
In this section, a survey of classical iterative algorithms, which use matrix inversion in every iteration is presented.Among the algorithms, which have been proposed for computing the principle square root of a matrix, the Newton methods have been studied for many years.Since the convergence is proved but numerical instability has been observed [4], several Newton variants have been derived.Newton methods have also been used for computing the th root a matrix [1,2].The conclusion that a nonsingular and diagonalizable matrix for any positive integer  has an th root and is stated in [10].
Algorithm 1(a).This algorithm is the classical Newton method [6] Note that we are also able to use the initial condition  0 = .Note that all the algorithms are applicable for  = 0, 1, . . .and that the convergence is achieved, when ‖ +1 −   ‖ < , where   is the sequence that converges to the principle square root of  and  is a small positive number and ‖‖ denotes the spectral norm of the matrix  (i.e., the largest singular value of ).

Algorithm 2(a)
. This algorithm is a variant of the classical Newton method [7] Note that we are also able to use the initial condition  0 = .

Algorithm 3(a)
. This algorithm is a variant of the classical Newton method, where the principle square root of  −1 is computed as well [6] ( Algorithm 4(a).This algorithm is proposed in [7] Algorithm 5(a).This algorithm is proposed in [5] Algorithm 6(a).It will be shown that the problem of computing the principal square root of a matrix is equivalent to the problem of solving a related Riccati equation.This algorithm takes advantage of this relation and is derived via the solution of the related Riccati equation.The problem of computing the principal square root of a matrix is associated with a well-known problem of estimation theory, namely the problem of solving a continuous time Riccati equation.The Riccati equation arises in linear estimation [3], namely in the implementation of the Kalman-Bucy filter and it is formulated as  ()  =  ()  () +  ()   () +  () −  ()   ()  −1 ()  ()  () , (0) =  0 , where () is the filtering error-covariance matrix, () and () are the system dynamic and output matrices, respectively.() and () are the plant and measurement noise covariance matrices, respectively.
In the special case of the time invariant model where () =  = 0, () = , () =  = , () =  and   () −1 () = , and using the initial condition  0 = 0, the Riccati equation of interest takes the following form It is well known [11,12] that there exists a steady state solution of the Riccati equation.
The following statement is now obviously true: "The problem of computing the principal square root of matrix  is equivalent to the problem of solving the Riccati equation (7), the steady state solution of which is equal to the principal square root of matrix ".
The idea is to apply the Bernoulli substitution [11] in the special case of the continuous time Riccati equation (7): an integration free solution of (7) can be obtained using the Bernoulli substitution [11]  () =  ()  −1 () .
Substituting ( 8) in (7), we have Then, by discretising the above equations, we have: Thus, by applying the Bernoulli substitution [11] in a special case of the continuous time Riccati equation, the following algorithm is derived Note that the following variant uses iterations only for the first quantity involved Algorithm 7(a).This algorithm is derived via the solution of the related Riccati equation using the doubling principle [12].
It is obvious that ( 10) and ( 11) can be easily written in a state space form as follows Let us denote: Then, using the doubling principle [12], that is, calculating the power it is easy to prove that the following equations hold Due to the form of matrix (0) given in (16), it is obvious that the matrices   and   are polynomials in  of the form where   and   are scalar coefficients and ℓ() and () are functions of the number of iterations, which need not be specified.Then, directly from the above equation we derive since Furthermore, it is easily seen from ( 15) and ( 16) that the following state space equation holds Then, it is obvious that: Then, using ( 18), ( 20)-( 22) and ( 24)-( 25), we have Thus, by using the doubling principle [12], the following algorithm is derived Algorithm 8(a).In the following, using ( 22) and ( 26)-(28), we derive with initial condition Note that, where we are obliged to use the initial condition  1 by (31), which is derived from ( 26)-(28) instead of using the initial condition  0 in (28), due to the fact that (30) requires inversion of matrix   and  0 = 0 is a singular matrix.At this point, we are able to use (30) with initial condition  0 = , resulting to the same sequence of   as in ( 26)-(28).
Then, it is obvious that the algorithm in ( 26)-( 28) can be written in the following equivalent form: Remark 1. (1) All the above algorithms are equivalent to each other.In fact, the relations between the quantities involved in the above algorithms hold: The proof is derived by induction and is trivial.
(2) The relationship between the Riccati equation and the principal square root algorithms is certified: it is shown that the algorithm derived by solving a special case of the continuous time Riccati equation is equivalent to the classical Newton method.

Inversion Free Iterative Algorithms
In this section, the inversion free iterative algorithms are developed by using the Schulz iteration for computing the inversion of a matrix or by eliminating matrix inversion from the convergence criteria of the algorithms.
Algorithm 1(b).This algorithm is an inversion free variant of the Newton method in Algorithm 1(a); the basic idea is to replace the computation of  −1  in Algorithm 1 by the Schulz iteration for  −1  as described in [13] where ‖  ‖ ∞ denotes the infinity norm of the matrix   .
Algorithm 2(b).This algorithm is an inversion free variant of Algorithm 2(a), derived using the Schulz iteration idea Algorithm 3(b).This algorithm is an inversion free variant of Algorithm 3(a), derived using the Schulz iteration idea Algorithm 4(b).This algorithm is an inversion free variant of Algorithm 4(a), derived using the Schulz iteration idea Algorithm 5(b).This algorithm is an inversion free variant of Algorithm 5(a), derived using the Schulz iteration idea Algorithm 6(b).This algorithm is an inversion free variant of Algorithm 6(a), derived using another approach to eliminate the matrix inversion requirement The convergence criterion depends on the quantity It is easy to prove that Hence, if ‖ − ‖ < 1, then ‖ +1 ‖ < ‖  ‖ < ‖‖; thus there exists ℓ : ‖ ℓ ‖ <  that can be calculated off-line.
Finally, we observe that Hence, if all the eigenvalues of  lie inside the unit circle, then the quantity   tends to zero as  tends to infinite.Furthermore, if the eigenvalues of  lie outside the unit circle, then we are able use the algorithm to calculate the principal square root of  −1 , from where we are able to find the principal square root of  (we need only two matrix inversions).

Algorithm 7(b).
This algorithm is an inversion free variant of Algorithm 7(a), derived using another approach to eliminate the matrix inversion requirement International Journal of Mathematics and Mathematical Sciences 5 Note that the following variant uses another convergence criterion The convergence criterion depends on the quantity It is easy to prove that Hence, if ‖−‖ < 1, then ‖ +1 ‖ < ‖  ‖ < 1; thus there exists ℓ : ‖ ℓ ‖ <  that can be calculated off-line.Finally, we observe that Hence, if all the eigenvalues of  lie inside the unit circle, then the quantity   tends to zero as  tends to infinite.Furthermore, if the eigenvalues of  lie outside the unit circle, then we are able use the algorithm to calculate the principal square root of  −1 , from where we are able to find the principal square root of  (we need only two matrix inversions).

Algorithm 8(b)
. This algorithm is an inversion free variant of Algorithm 8(a), derived using the Schulz iteration idea: All the inversion use and inversion free iterative algorithms are summarized in Table 1.

Remark 2. (1)
All the algorithms presented in Section 3 are inversion free algorithms; they do not require matrix inversion at every iteration.The elimination of the matrix inversion requirement can be achieved using the Schulz iteration.
(2) Algorithms 6(b) and 7(b) use another technique to eliminate the matrix inversion requirement.These two algorithms require only one matrix inversion, after the convergence of the algorithms in order to compute the final result.It is obvious that this matrix inversion can be substituted by the Schulz iteration.
(3) Algorithm 8(b) uses only one inversion ( −1 ).It is obvious that this matrix inversion can be substituted by the Schulz iteration.
(4) Note that all the proposed inversion free iterative algorithms compute the principal square root of an  ×  dimensional matrix with complexity of the order ( 3 ), due to the fact that the proposed inversion algorithms require matrix additions and matrix multiplications.The algorithms avoid possible problems in matrix inversion (and even at every iteration).

Algebraic Method
In this section an algebraic method for computing the principal square root of a matrix is presented.The method is associated with the solution of a special case of the continuous time Riccati equation which arises in linear estimation [9] and requires only one matrix inversion (this is not an iterative method).
We are able to solve the Riccati equation ( 7) using the algebraic method proposed in [14][15][16][17].In fact the following matrix is formed: Then, Φ can be written as where is a block-diagonal matrix containing the eigenvalues of Φ, with Λ diagonal matrix with all the eigenvalues of Φ lying the right half-plane eigenvalues of matrix Φ and is the matrix containing the corresponding eigenvectors of Φ.
Then, the solution of the Riccati equation is given in terms of the eigenvalues of matrix Φ and formulated The method requires only one matrix inversion in order to compute the final result (this is not an iterative method).It is obvious that this matrix inversion can be substituted by the Schulz iteration.
We are also able to compute all the square roots (not only the principal square root) of a matrix, using the ideas in [14][15][16].

Simulation Results
Simulation results are given to illustrate the efficiency of the proposed methods.The proposed algorithms compute The proposed algebraic method has also been applied and computes the accurate principal square root of .
All the iterative algorithms have been applied and compute the same principal square root of : Algorithm 7(b) does not converge to the principal square root of .But, the algorithm has been used to calculate the principal square root of  −1 , from where we are able to find the principal square root of  (we need only two matrix inversions).The proposed algebraic method has also been applied and computes the accurate principal square root of .
All the iterative algorithms have been applied and compute the same principal square root of : Algorithm 4(b) does not converge to the principal square root of .But, it has been used to calculate the principal square root of  −1 , from where we are able to find the principal square root of  (we need only two matrix inversions).
The proposed algebraic method has also been applied and computes the accurate principal square root of.
We are able to conclude that Algorithms 6 and 7 are not always stable.Note that if any algorithm is not stable we are able to try to compute the principal square root of  −1 , from where we are able to find the principal square root of  (we need only two matrix inversions).Finally, the proposed algebraic method has also been applied and computes the accurate principal square root of .

Conclusions
In this paper, a survey of classical iterative algorithms for computing the principal square root of a matrix is presented; these algorithms require matrix inversion at every iteration.New algorithms for computing the principal square root of a matrix are proposed.The novelty of this work is the derivation of inversion free algorithms.The elimination of the matrix inversion requirement can be achieved by using the Schulz iteration.The proposed algorithms are derived by using the Newton method or by applying the Bernoulli substitution (in a special case of the continuous time Riccati equation) or by using the doubling principle.An inversion free algebraic algorithm is also derived.
Note that, among the proposed algorithms, the one per step iterative algorithm, the one doubling iterative algorithm and the algebraic algorithm, are associated with the solution of a special case of the continuous time Riccati equation, which arises in linear estimation.The relations between these algorithms and the classical Newton method are established.
All the proposed inversion free iterative algorithms compute the principal square root of an  ×  matrix achieving ( 3 ) complexity, due to the fact that the proposed inversion algorithms require matrix additions and matrix multiplications.The algorithms avoid matrix inversion.This means that they avoid possible problems in matrix inversion (at every iteration) and that they are easily programmable.

Table 1 :
Principal square root iterative algorithms.
is not stable: it converges to the right value of the principal square root of  after 7 iterations, but after 144 iterations it becomes to diverge and then it converges after 161 iterations to a value, which is not the principal square root of .The proposed algebraic method has also been applied and computes the accurate principal square root of .