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

Authors: Nicholas Assimakis; Maria Adam
Journal: International Journal of Mathematics and Mathematical Sciences
 (2014)

---

## Abstract

New algorithms are presented about the principal square root of ann×n matrix A. 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.


---

## Body

## 1. Introduction

LetA be a real or complex n×n matrix with no eigenvalues on R- (the closed negative real axis). Then there exists a unique n×n matrix X such that X2=A and the eigenvalues of X lie in the segment {z:-π/2<arg(z)<π/2}. We refer to X as theprinciple square root of A.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–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 Section2 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.


## 2. 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 pth root a matrix [1, 2]. The conclusion that a nonsingular and diagonalizable matrix for any positive integer n has an nth root and is stated in [10].Algorithm  1(a). This algorithm is the classical Newton method [6]
(1)Yn+1=12(Yn+Yn-1A),Y0=A,Yn⟶A1/2.
Note that we are also able to use the initial condition Y0=I.Note that all the algorithms are applicable forn=0,1,… and that the convergence is achieved, when ∥Qn+1-Qn∥<ε, where Qn is the sequence that converges to the principle square root of A and ε is a small positive number and ∥M∥ denotes the spectral norm of the matrix M (i.e., the largest singular value of M).Algorithm  2(a). This algorithm is a variant of the classical Newton method [7]
(2)Zn+1=12(Zn+AZn-1),Z0=A,Zn⟶A1/2.
Note that we are also able to use the initial condition Z0=I.Algorithm  3(a). This algorithm is a variant of the classical Newton method, where the principle square root of A-1 is computed as well [6]
(3)Mn+1=12(Mn+Nn-1),M0=A,Mn⟶A1/2,Nn+1=12(Nn+Mn-1),N0=I,Nn⟶A-1/2.Algorithm  4(a). This algorithm is proposed in [7]
(4)Gn+1=4Gn(I+Gn)-2,G0=A,Gn⟶I,Rn+1=12Rn(I+Gn),R0=I,Rn⟶A1/2.Algorithm  5(a). This algorithm is proposed in [5]
(5)Xn+1=12I+Xn-1-12Xn-2,X0=12(I+A),Xn⟶I,Wn+1=WnXn,W0=I,Wn⟶A1/2.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
(6)dP(t)dt=F(t)P(t)+P(t)FT(t)+Q(t)-P(t)HT(t)R-1(t)H(t)P(t),P(0)=P0,
where P(t) is the filtering error-covariance matrix, F(t) and H(t) are the system dynamic and output matrices, respectively. Q(t) and R(t) are the plant and measurement noise covariance matrices, respectively.In the special case of the time invariant model whereF(t)=F=0, H(t)=H, Q(t)=Q=A, R(t)=R and HT(t)R-1(t)H=I, and using the initial condition P0=0, the Riccati equation of interest takes the following form
(7)dP(t)dt=A-P(t)P(t),P(0)=0.
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 matrixA 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 A”.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]
(8)P(t)=V(t)U-1(t).
Substituting (8) in (7), we have
(9)dU(t)dt=V(t),U(0)=I,dV(t)dt=AU(t),V(0)=0.
Then, by discretising the above equations, we have:
(10)Uk+1=Uk+Vk,U0=I,(11)Vk+1=AUk+Vk,V0=0,(12)Pk+1=Vk+1+Uk+1-1,P0=0.
Thus, by applying the Bernoulli substitution [11] in a special case of the continuous time Riccati equation, the following algorithm is derived
(13)Un+1=Un+Vn,U0=I,Vn+1=AUn+Vn,V0=0,Pn=VnUn-1,Pn⟶A1/2.
Note that the following variant uses iterations only for the first quantity involved
(14)Un+1=2Un+(A-I)Un-1,U1=I,U0=I,Vn+1=Un+1+(A-I)Un,Pn=VnUn-1,Pn⟶A1/2.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
(15)[Uk+1Vk+1]=Φ[UkVk],Φ=[IIAI].
Let us denote:
(16)S(0)2k=S(k)=[akbkAbkak],S(0)=[IIAI].
Then, using the doubling principle [12], that is, calculating the power
(17)S(k+1)=S2(k)=S(k)S(k),
it is easy to prove that the following equations hold
(18)ak+1=akak+bkAbk,bk+1=akbk+bkak,Abk+1=Abkak+akAbk,ak+1=Abkbk+akak.
Due to the form of matrix S(0) given in (16), it is obvious that the matrices ak and bk are polynomials in A of the form
(19)ak=∑i=0l(k)xiAi,bk=∑j=0m(k)yjAj,
where xi and yj are scalar coefficients and l(k) and m(k) are functions of the number of iterations, which need not be specified. Then, directly from the above equation we derive
(20)akA=Aak,(21)bkA=Abk,(22)akbk=bkak,
since
(23)akbk=∑i=0l(k)xiAi∑j=0m(k)yjAj=∑j=0m(k)yjAj∑i=0l(k)xiAi=bkak.
Furthermore, it is easily seen from (15) and (16) that the following state space equation holds
(24)[U2kV2k]=S(0)2k[U0V0]=[akbkAbkak][I0].
Then, it is obvious that:
(25)U2k=ak,V2k=Abk,P2k=Pk=V2kU2k-1.
Then, using (18), (20)–(22) and (24)-(25), we have
(26)ak+1=ak2+Abk2,a0=I,(27)bk+1=2akbk,b0=I,(28)Pk+1=Abk+1ak+1-1,P0=0.
Thus, by using the doubling principle [12], the following algorithm is derived
(29)an+1=an2+Abn2,a0=I,bn+1=2anbn,b0=I,Pn=Abnan-1,Pn⟶A1/2.Algorithm  8(a). In the following, using (22) and (26)–(28), we derive
(30)Pk+1=2(Pk-1+PkA-1)-1,
with initial condition
(31)P1=2A(I+A)-1=2(A-1+I)-1.
Note that, where we are obliged to use the initial condition P1 by (31), which is derived from (26)–(28) instead of using the initial condition P0 in (28), due to the fact that (30) requires inversion of matrix Pk and P0=0 is a singular matrix.At this point, we are able to use (30) with initial condition P0=A, resulting to the same sequence of Pk as in (26)–(28).Then, it is obvious that the algorithm in (26)–(28) can be written in the following equivalent form:
(32)Pn+1=2(Pn-1+PnA-1)-1,P0=A,Pn⟶A1/2.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:
(33)Yn=Zn=Mn=ANn=Gn=Wn=anbn-1=Pn-1A.
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.


## 3. 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 Yn-1 in Algorithm 1 by the Schulz iteration for Yn-1 as described in [13]
(34)Yn+1=12(Yn+XnA),Y0=A,Yn⟶A1/2,Xk+1=Xk(2I-YnXk),X0=I∥Yn∥∞,Xk⟶Xn,
where ∥Yn∥∞ denotes the infinity norm of the matrix Yn.Algorithm  2(b). This algorithm is an inversion free variant of Algorithm  2(a), derived using the Schulz iteration idea
(35)Zn+1=12(Zn+AXn),Z0=A,Zn⟶A1/2,Xk+1=Xk(2I-ZnXk),X0=I∥Zn∥∞,Xk⟶Xn.Algorithm  3(b). This algorithm is an inversion free variant of Algorithm  3(a), derived using the Schulz iteration idea
(36)Mn+1=12(Mn+Xn),M0=A,Mn⟶A1/2,Xk+1=Xk(2I-NnXk),X0=I∥Nn∥∞,Xk⟶Xn,Nn+1=12(Nn+Yn),N0=I,Nn⟶A-1/2,Yk+1=Xk(2I-MnYk),Y0=I∥Mn∥∞,Yk⟶Yn.Algorithm  4(b). This algorithm is an inversion free variant of Algorithm  4(a), derived using the Schulz iteration idea
(37)Gn+1=4GnXn2,G0=A,Gn⟶I,Xn=(I+Gn),Yk+1=Yk(2I-XnYk),Y0=I∥Xn∥∞,Yk⟶Xn,Rn+1=12Rn(I+Gn),R0=I,Rn⟶A1/2.Algorithm  5(b). This algorithm is an inversion free variant of Algorithm  5(a), derived using the Schulz iteration idea
(38)Xn+1=12I+Yn-12Yn2,X0=12(I+A),Xn⟶I,Yk+1=Yk(2I-XnYk),Y0=I∥Xn∥∞,Yk⟶Yn,Wn+1=WnXn,W0=I,Wn⟶A1/2.Algorithm  6(b). This algorithm is an inversion free variant of Algorithm  6(a), derived using another approach to eliminate the matrix inversion requirement
(39)Un+1=Un+Vn,U0=I,Vn+1=AUn+Vn,V0=0,δn=Vn+1Un-VnUn+1⟶0,VnUn-1⟶A1/2.
The convergence criterion depends on the quantity
(40)δn=Vn+1Un-VnUn+1.
It is easy to prove that
(41)δn+1=(I-A)δn,δ0=A.
Hence, if ∥I-A∥<1, then ∥δn+1∥<∥δn∥<∥A∥; thus there exists l:∥δl∥<ε that can be calculated off-line.Finally, we observe that(42)δn=(I-A)nA.
Hence, if all the eigenvalues of A lie inside the unit circle, then the quantity δn tends to zero as n tends to infinite. Furthermore, if the eigenvalues of A lie outside the unit circle, then we are able use the algorithm to calculate the principal square root of A-1, from where we are able to find the principal square root of A (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
(43)an+1=an2+Abn2,a0=I,bn+1=2anbn,b0=I,bn+1an-bnan+1⟶0,Abnan-1⟶A1/2.
Note that the following variant uses another convergence criterion
(44)an+1=an2+Abn2,a0=I,bn+1=2anbn,b0=I,en+1=en2,e0=I-A,en=an2-Abn2⟶0,Abnan-1⟶A1/2.
The convergence criterion depends on the quantity
(45)en=an2-Abn2.
It is easy to prove that
(46)en+1=en2,e0=I-A.
Hence, if ∥I-A∥<1, then ∥en+1∥<∥en∥<1; thus there exists l:∥el∥<ε that can be calculated off-line.Finally, we observe that(47)en=(I-A)2n.
Hence, if all the eigenvalues of A lie inside the unit circle, then the quantity en tends to zero as n tends to infinite. Furthermore, if the eigenvalues of A lie outside the unit circle, then we are able use the algorithm to calculate the principal square root of A-1, from where we are able to find the principal square root of A (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:
(48)Pn+1=2Rn,P0=A,Pn⟶A1/2,Zn=(Xn+PnA-1),Xk+1=Xk(2I-ZnXk),X0=I∥Zn∥∞,Xk⟶Rn,Yk+1=Xk(2I-PnYk),Y0=I∥Pn∥∞,Yk⟶Xn.
All the inversion use and inversion free iterative algorithms are summarized in Table 1.Table 1
Principal square root iterative algorithms.
Inversion use
Inversion free
1a
Y
n
+
1
=
1
2
(
Y
n
+
Y
n
-
1
A
)
,
Y
0
=
A
,
Y
n
⟶
A
1
/
2
1b
Z
n
+
1
=
1
2
(
Z
n
+
A
X
n
)
,
{
∑
W
W
w
}
Z
0
=
A
,
Z
n
⟶
A
1
/
2
X
k
+
1
=
X
k
(
2
I
-
Z
n
X
k
)
,
X
0
=
I
∥
Z
n
∥
∞
,
X
k
⟶
X
n
2a
Z
n
+
1
=
1
2
(
Z
n
+
A
Z
n
-
1
)
,
Z
0
=
A
,
Z
n
⟶
A
1
/
2
2b
Z
n
+
1
=
1
2
(
Z
n
+
A
X
n
)
,
Z
0
=
A
,
Z
n
⟶
A
1
/
2
X
k
+
1
=
X
k
(
2
I
-
Z
n
X
k
)
,
X
0
=
I
∥
Z
n
∥
∞
,
X
k
⟶
X
n
3a
M
n
+
1
=
1
2
(
M
n
+
N
n
-
1
)
,
M
0
=
A
,
M
n
⟶
A
1
/
2
N
n
+
1
=
1
2
(
N
n
+
M
n
-
1
)
,
{
∑
W
W
w
}
N
0
=
I
,
N
n
⟶
A
-
1
/
2
3b
M
n
+
1
=
1
2
(
M
n
+
X
n
)
,
M
0
=
A
,
M
n
⟶
A
1
/
2
X
k
+
1
=
X
k
(
2
I
-
N
n
X
k
)
,
X
0
=
I
∥
N
n
∥
∞
,
X
k
⟶
X
n
N
n
+
1
=
1
2
(
N
n
+
Y
n
)
,
N
0
=
I
,
N
n
⟶
A
-
1
/
2
Y
k
+
1
=
X
k
(
2
I
-
M
n
Y
k
)
,
Y
0
=
I
∥
M
n
∥
∞
,
Y
k
⟶
Y
n
4a
G
n
+
1
=
4
G
n
(
I
+
G
n
)
-
2
,
G
0
=
A
,
G
n
⟶
I
R
n
+
1
=
1
2
R
n
(
I
+
G
n
)
,
R
0
=
I
,
R
n
⟶
A
1
/
2
4b
G
n
+
1
=
4
G
n
X
n
2
,
G
0
=
A
,
G
n
⟶
I
X
n
=
(
I
+
G
n
)
Y
k
+
1
=
Y
k
(
2
I
-
X
n
Y
k
)
,
Y
0
=
I
∥
X
n
∥
∞
,
Y
k
⟶
X
n
R
n
+
1
=
1
2
R
n
(
I
+
G
n
)
,
R
0
=
I
,
R
n
⟶
A
1
/
2
5a
X
n
+
1
=
1
2
I
+
X
n
-
1
-
1
2
X
n
-
2
,
X
0
=
1
2
(
I
+
A
)
,
X
n
⟶
I
W
n
+
1
=
W
n
X
n
,
W
0
=
I
,
W
n
⟶
A
1
/
2
5b
X
n
+
1
=
1
2
I
+
Y
n
-
1
2
Y
n
2
,
X
0
=
1
2
(
I
+
A
)
,
X
n
⟶
I
Y
k
+
1
=
Y
k
(
2
I
-
X
n
Y
k
)
,
Y
0
=
I
∥
X
n
∥
∞
,
Y
k
⟶
Y
n
W
n
+
1
=
W
n
X
n
,
W
0
=
I
,
W
n
⟶
A
1
/
2
6a
U
n
+
1
=
U
n
+
V
n
,
U
0
=
I
V
n
+
1
=
A
U
n
+
V
n
,
V
0
=
0
P
n
=
V
n
U
n
-
1
P
n
⟶
A
1
/
2
6b
U
n
+
1
=
U
n
+
V
n
,
U
0
=
I
V
n
+
1
=
A
U
n
+
V
n
,
V
0
=
0
δ
n
=
V
n
+
1
U
n
-
V
n
U
n
+
1
⟶
0
V
n
U
n
-
1
⟶
A
1
/
2
7a
a
n
+
1
=
a
n
2
+
A
b
n
2
,
a
0
=
I
b
n
+
1
=
2
a
n
b
n
,
b
0
=
I
P
n
=
A
b
n
a
n
-
1
P
n
⟶
A
1
/
2
7b
a
n
+
1
=
a
n
2
+
A
b
n
2
,
a
0
=
I
b
n
+
1
=
2
a
n
b
n
,
b
0
=
I
e
n
+
1
=
e
n
2
,
e
0
=
I
-
A
e
n
=
a
n
2
-
A
b
n
2
⟶
0
A
b
n
a
n
-
1
⟶
A
1
/
2
8a
P
n
+
1
=
2
(
P
n
-
1
+
P
n
A
-
1
)
-
1
,
P
0
=
A
,
P
n
⟶
A
1
/
2
8b
P
n
+
1
=
2
R
n
,
P
0
=
A
,
P
n
⟶
A
1
/
2
Z
n
=
(
X
n
+
P
n
A
-
1
)
X
k
+
1
=
X
k
(
2
I
-
Z
n
X
k
)
,
X
0
=
I
∥
Z
n
∥
∞
,
X
k
⟶
R
n
Y
k
+
1
=
X
k
(
2
I
-
P
n
Y
k
)
,
Y
0
=
I
∥
P
n
∥
∞
,
Y
k
⟶
X
nRemark 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 (A-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 n×n dimensional matrix with complexity of the order O(n3), 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).


## 4. 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–17]. In fact the following matrix is formed:
(49)Φ~=[0IA0].
Then, Φ~ can be written as
(50)Φ~=WLW-1,
where
(51)L=[Λ00-Λ]
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
(52)W=[W11W12W21W22]
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
(53)X=W21W11-1.
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–16].


## 5. Simulation Results

Simulation results are given to illustrate the efficiency of the proposed methods. The proposed algorithms compute accurate solutions as verified trough the following simulation examples.All the algorithms are considered inversion free iterative algorithms, when Moore-Penrose matrix inversion is applied.Example 3.
Consider the 2 × 2 symmetric and positive definite matrix(54)A=[0.90.50.51.1].
All the inversion use and inversion free iterative algorithms have been applied with convergence criterion ε=10-6 and compute the same principal square root of A:
(55)A1/2=[0.91260.25920.25921.0163].
The proposed algebraic method has also been applied and computes the accurate principal square root of  A.Example 4.
Consider the2×2 matrix
(56)A=[1.00.50.24.0].
All the iterative algorithms have been applied and compute the same principal square root of A:
(57)A1/2=[0.99440.16710.06691.9972].
Algorithm  7(b) does not converge to the principal square root of A. But, the algorithm has been used to calculate the principal square root of A-1, from where we are able to find the principal square root of A (we need only two matrix inversions).
The proposed algebraic method has also been applied and computes the accurate principal square root ofA.Example 5.
This example is taken from [2]. Consider the 3×3 matrix
(58)A=[0.60.30.10.20.70.10.10.10.8].
All the iterative algorithms have been applied and compute the same principal square root of A:
(59)A1/2=[0.75720.18830.05440.12480.82080.05440.05670.05220.8911].
Algorithm 6(b) is not stable: it converges to the right value of the principal square root of A 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 A.
The proposed algebraic method has also been applied and computes the accurate principal square root ofA.Example 6.
Consider the4×4 matrix
(60)A=[10.9900.9810.9470.89010.9800.7650.9810.98010.3950.9420.9610.9451].
All the iterative algorithms have been applied and compute the same principal square root of A:
(61)A1/2=[0.59050.41750.42001.59530.26000.62280.51390.58070.71440.58380.6952-0.21910.40770.40020.38040.7892].
Algorithms 6(a), 6(b), 7(a), and 7(b) are not stable.
The proposed algebraic method has also been applied and computes the accurate principal square root ofA.Example 7.
Consider the6×6 matrix
(62)A=[4891271263191310233152753211246213328].
All the iterative algorithms have been applied and compute the same principal square root of A:(63)A1/2=[1.89752.46611.2370-0.12790.52200.09310.00181.13961.17930.4085-0.19242.27850.13710.57113.01040.27450.55810.17480.13801.25900.09042.4991.2067-0.27510.3989-0.0282-0.02100.31481.78381.31160.3744-0.21780.46650.54870.18012.8800].Algorithm 4(b) does not converge to the principal square root of A. But, it has been used to calculate the principal square root of A-1, from where we are able to find the principal square root of A (we need only two matrix inversions).
The proposed algebraic method has also been applied and computes the accurate principal square root ofA.
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 ofA-1, from where we are able to find the principal square root of A (we need only two matrix inversions). Finally, the proposed algebraic method has also been applied and computes the accurate principal square root of A.


## 6. 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 ann×n matrix achieving O(n3) 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.



---

*Source: 613840-2014-06-18.xml*