# A Numerical Method for Computing the Principal Square Root of a Matrix

Authors: F. Soleymani; S. Shateyi; F. Khaksar Haghani
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

It is shown how the mid-point iterative method with cubical rate of convergence can be applied for finding the principal matrix square root. Using an identity between matrix sign function and matrix square root, we construct a variant of mid-point method which is asymptotically stable in the neighborhood of the solution. Finally, application of the presented approach is illustrated in solving a matrix differential equation.


---

## Body

## 1. Introductory Notes

Let us consider a scalar functionf and a matrix A
∈
C
n
×
n and specify f
(
A
) to be a matrix function of the same dimensions as A. Higham in [1] indicates that if f is defined on the spectrum of A, then one has the following properties for the matrix function f
(
A
):(i)
f
(
A
) commutes with A;
(ii)
f
(
A
T
)
=
f
(
A
)
T;
(iii)
f
(
X
A
X
-
1
)
=
X
f
(
A
)
X
-
1;
(iv)
the eigenvalues off
(
A
) are f
(
λ
i
), where λ
i are the eigenvalues of A;
(v)
ifX commutes with A, then X commutes with f
(
A
).In this paper, we are interested in numerical computation of matrix square root, which is one of the most fundamental matrix functions with potential applications. An example involving the computation of matrix square roots appears in a semidefinite programming algorithm [2]. A step length is calculated using a line search algorithm that involves computing the unique symmetric positive definite (SPD) square roots of a sequence of SPD matrices, where each matrix differs little from the previous one. That is, given a SPD A
∈
R
n
×
n and the SPD square root X of A, one wishes to find the SPD square root X
¯ of A
¯, where ∥
A
¯
-
A
∥ is small. This has to be done repeatedly with A changing slightly each time.It is known that [3] A
∈
C
n
×
n has a square root if and only if in the ascent sequence of integers d
1
,
d
2
,
…
, defined by
(1)
d
i
=
dim
⁡
(
null
(
A
i
)
)
-
dim
⁡
(
null
(
A
i
-
1
)
)
,
no two terms are the same odd integer.Among the square roots of a matrix, the principal square root is distinguished by its usefulness in theory and applications. LetA
∈
C
n
×
n have no eigenvalues on R
- (the closed negative real axis). There is a unique square root X of A that all of its eigenvalues lie in the open right half-plane, and it is a primary matrix function of A [4]. We refer to X as the principal square root of A and write X
=
A
1
/
2. Note that if A is real, then A
1
/
2 is real.Björck and Hammarling [5] developed a numerical method for the principal square root via Schur form. Among all the available numerical procedures, the Newton’s matrix method and its variants have received many attractions in order to tackle this problem. To be more precise, applying the famous Newton’s method to the matrix equation
(2)
F
(
X
)
≡
X
2
-
A
=
0
would yield the following (simplified) matrix iteration method [4]:
(3)
X
k
+
1
=
1
2
[
X
k
+
A
X
k
-
1
]
,
for suitable initial values X
0, wherein the commutativity of matrices X
k and A has been considered to simplify the procedure as much as possible. Unfortunately, this iteration becomes unstable in most cases as pointed out by Higham [4].To remedy this, Denman and Beavers (DB) in [6] proposed a quadratically convergent variant of Newton’s matrix method (3) which is numerically stable in what follows:
(4)
Y
0
=
A
,
Z
0
=
I
,
k
=
0,1
,
…
,
Y
k
+
1
=
1
2
[
Y
k
+
Z
k
-
1
]
,
Z
k
+
1
=
1
2
[
Z
k
+
Y
k
-
1
]
.
This method generates the sequence {
Y
k
}, which converges to A
1
/
2. For further readings, one may refer to [7, 8].In the next section, we first apply the following mid-point cubically convergent scheme of Frontini and Sormani (given for scalar nonlinear equations [9])
(5)
y
k
=
x
k
-
2
-
1
f
′
(
x
k
)
-
1
f
(
x
k
)
,
x
k
+
1
=
x
k
-
f
′
(
y
k
)
-
1
f
(
x
k
)
,
to derive a new iterative scheme for the numerical computation of A
1
/
2. In Section 2, we theoretically prove that the scheme is convergent under some proper conditions. A variant of the mid-point rule (5) is brought forward to provide a stable iteration for finding both the principal matrix square root and the principal inverse matrix square root. We note that our results in this work are some extensions to [10]. In order to derive a higher order new version of mid-point method for computing square root, we use some of the methods in the literature and obtain a new scheme for A
1
/
2 as well. Section 3 is devoted to the application of the discussed matrix methods in solving some examples. Finally, Section 4 draws a conclusion of this paper.


## 2. Extension of the Mid-Point Method

Multipoint nonlinear solvers such as (5) belong to the class of the most efficient methods for nonlinear scalar equations. Recent interests in the research and development of this type of methods have arisen from their capability to overcome theoretical limits of one-point methods concerning the convergence order and computational efficiency.Let us apply (5) for solving the nonlinear matrix equation (2). We attain the following scheme:
(6)
X
k
+
1
=
[
3
A
X
k
+
X
k
3
]
[
A
+
3
X
k
2
]
-
1
,
considering the fact that the initial approximation is of the following forms:
(7)
X
0
=
α
I
,
or
(8)
X
0
=
β
A
,
wherein α and β are suitable positive parameters to ensure convergence. The derivation of (6) may not be clear yet. Hence, we give a Mathematica code written for obtaining it as fast as possible in a symbolical language as follows:ClearAll["Global‘*"]f[X_]:= X^2 - A;Y=X - f[X]/(2f’[X]) // FullSimplify;FullSimplify[X - f[X]/f’[Y]]Lemma 1.
LetA
=
[
a
i
j
]
n
×
n be a nonsingular complex matrix with no eigenvalues on R
-. If A
X
0
=
X
0
A is valid, then for the sequence of {
X
k
}
k
=
0
k
=
∞ of (6), one has that
(9)
A
X
k
=
X
k
A
holds, for all k
=
1,2
,
….Proof.
The proof is by induction. LetX
0 be one of the initial values (7) or (8); then it is obvious that A
X
0
=
X
0
A. For the inductive hypothesis, we take k
>
0 and A
X
k
=
X
k
A and show that
(10)
A
X
k
+
1
=
A
[
3
A
X
k
+
X
k
3
]
[
A
+
3
X
k
2
]
-
1
=
[
3
A
2
X
k
+
A
X
k
3
]
[
A
+
3
X
k
2
]
-
1
=
[
3
(
X
k
A
)
A
+
X
k
2
(
X
k
A
)
]
[
A
+
3
X
k
2
]
-
1
=
[
3
X
k
A
+
X
k
3
]
[
A
-
1
A
+
3
A
-
1
X
k
2
]
-
1
=
[
3
X
k
A
+
X
k
3
]
[
A
A
-
1
+
3
X
k
2
A
-
1
]
-
1
=
X
k
+
1
A
.
Note that, from A
X
k
=
X
k
A, we obtain that A
-
1
X
k
=
X
k
A
-
1. This completes the proof.Let(
M
,
d
) be a metric space. We say that a sequence {
X
k
} of points of M converges to a point X
* of M if
(11)
lim
⁡
k
→
∞
d
(
X
k
,
X
*
)
=
0
.
The convergence properties will depend on the choice of a distance in M, but for a given distance the speed of convergence of the sequence {
X
k
} is characterized by the speed of convergence of the sequence of nonnegative numbers d
(
X
k
,
X
*
) [11]. In what follows, we investigate the rate of convergence for (6).Theorem 2.
LetA
=
[
a
i
,
j
]
n
×
n have no eigenvalues on R
-. If the initial approximation X
0 is according to the rules (7) or (8), then the iterative method (6) converges to the principal A
1
/
2 with local third order of convergence (without considering the round-off error).Proof.
Using the methodology of [10], we now prove the convergence order of (6). If A is considered to be diagonalizable and nonsingular, then
(12)
V
-
1
A
V
=
D
=
diag
⁡
{
λ
1
,
…
,
λ
n
}
,
where V is a nonsingular matrix. In addition assume that M
k
=
V
-
1
X
k
V. From the iterative method (6), we attain
(13)
M
k
+
1
=
[
3
D
M
k
+
M
k
3
]
[
D
+
3
M
k
2
]
-
1
.
{
M
k
} is the sequence of real block diagonal matrices M
k
=
diag
⁡
{
m
1
(
k
)
,
…
,
m
n
(
k
)
}. Note that if M
0 is a diagonal matrix, then all successive M
k are diagonal too. The recursive relation (13) is equivalent to n scalar equations as follows:
(14)
m
i
(
k
+
1
)
=
[
3
λ
i
m
i
(
k
)
+
(
m
i
(
k
)
)
3
]
[
λ
i
+
3
(
m
i
(
k
)
)
2
]
-
1
,
1
≤
i
≤
n
.
The sequence {
m
i
(
k
)
} is real and it is obvious that
(15)
lim
⁡
k
→
∞
M
k
=
D
1
/
2
=
diag
⁡
{
λ
1
,
…
,
λ
n
}
.
Then the fact that M
k
=
V
-
1
X
k
V yields lim
⁡
k
→
∞
X
k
=
V
D
1
/
2
V
-
1
=
A
1
/
2. Therefore,
(16)
M
k
+
1
-
D
1
/
2
=
1
8
(
3
M
k
+
D
1
/
2
)
M
k
-
3
(
M
k
-
D
1
/
2
)
3
.
Using a transformation by V, we have
(17)
X
k
+
1
-
A
1
/
2
=
1
8
(
3
X
k
+
A
1
/
2
)
X
k
-
3
(
X
k
-
A
1
/
2
)
3
.
Taking norm of (17), we obtain that
(18)
∥
X
k
+
1
-
A
1
/
2
∥
≤
1
8
∥
3
X
k
+
A
1
/
2
∥
∥
X
k
-
1
∥
3
∥
X
k
-
A
1
/
2
∥
3
.
Therefore, the method (6) possesses the local third order of convergence. Equation (11) reveals, in a metric space, how convergence could be understood from (17). Furthermore since A has no nonpositive real eigenvalues, so the mid-point approach computes the principal matrix square root.
We remark that all iterates{
X
k
} generated by (6) are well defined when the function f
(
x
)
=
a
+
3
x
2 is not zero for all nonnegative numbers x. In fact, the initial values actually determine a unique real number x
k for each k
≥
1. In general, well defined means that some object was given a description up to some arbitrary choices but that the arbitrary choices have no material effect. Accordingly, using the initial choices A and I, the sequence is well defined.Under the assumptions of Theorem2. (without considering the round-off error), and when A is nonsingular, X
k tends to A
1
/
2. This means that X
k
-
1 tends to A
-
1
/
2, that is, the principal inverse matrix square root.Unfortunately, scheme (6) has turned out to be unstable when the size of the input matrix A is high or is ill-conditioned. A similar reasoning as in [4] could reveal this point using X
0
=
A. Therefore, (6) can only be employed for symmetric positive definite matrices with extremely well-conditioned structure.Iannazzo in [12] showed that in the case that X
0 is the identity matrix, the Newton’s method converges for any matrix A
having eigenvalues with modulus less than 1 and with positive real parts. Similar results can be deduced for (6).In order to reach third order of convergence by the mid-point approach in solving (2) as well as to obtain a stable sequence, we use a fact (identity) in what follows:
(19)
sign
⁡
(
[
0
A
I
0
]
)
=
[
0
A
1
/
2
A
-
1
/
2
0
]
,
which indicates an important relationship between matrix square root A
1
/
2 and the matrix sign function sign
⁡
(
A
).We remark that (19) has an extra important advantage which is the computation of the inverse matrix square root along with the matrix square root at the same time.Applying (5) to solve the matrix equation G
(
X
)
≡
X
2
-
I
=
0 yields in its reciprocal form
(20)
X
k
+
1
=
[
I
+
3
X
k
2
]
[
X
k
(
3
I
+
X
k
2
)
]
-
1
.
This method which is similar to the application of Halley’s method for solving X
2
=
I (see, e.g., [13, 14]) computes the matrix sign of the nonsingular matrix
(21)
[
0
A
I
0
]
and in the convergence phase gives
(22)
[
0
A
1
/
2
A
-
1
/
2
0
]
.
It is easy to deduce the following convergence theorem for (20).Theorem 3.
LetA
=
[
a
i
,
j
]
n
×
n have no eigenvalues on R
-. If X
0 is (21), then the iterative method (20) converges to (22).Proof.
The proof of this theorem is based on the diagonalization property of the sign matrixS. It is hence omitted.Lemma 4.
The sequence{
X
k
}
k
=
0
k
=
∞ generated by (20) using X
0
= (21) is asymptotically stable.Proof.
LetΔ
k be the numerical perturbation introduced at the kth iterate of (20). Next, one has
(23)
X
~
k
=
X
k
+
Δ
k
.
Here, we perform a first-order error analysis; that is, we formally neglect quadratic terms such as (
Δ
k
)
2. We have
(24)
X
~
k
+
1
=
[
I
+
3
X
~
k
2
]
[
X
~
k
(
3
I
+
X
~
k
2
)
]
-
1
=
[
I
+
3
(
X
k
+
Δ
k
)
2
]
[
(
X
k
+
Δ
k
)
(
3
I
+
(
X
k
+
Δ
k
)
2
)
]
-
1
,
where the following identity will be used (for any nonsingular matrix B and the matrix C):
(25)
(
B
+
C
)
-
1
≈
B
-
1
-
B
-
1
C
B
-
1
.
Further simplifying yields to
(26)
X
~
k
+
1
≈
(
I
+
3
(
S
+
Δ
k
)
2
)
(
(
S
+
Δ
k
)
(
3
I
+
(
S
+
Δ
k
)
2
)
)
-
1
≈
(
4
I
+
3
S
Δ
k
+
3
Δ
k
S
)
(
(
S
+
Δ
k
)
(
4
I
+
S
Δ
k
+
Δ
k
S
)
)
-
1
≈
(
4
I
+
3
S
Δ
k
+
3
Δ
k
S
)
(
1
4
S
-
5
16
S
Δ
k
S
-
1
16
Δ
k
)
,
where S
2
=
I, S
-
1
=
S, and, for enough large k, we have considered X
k
≈
S. After some algebraic manipulations and using Δ
k
+
1
=
X
~
k
+
1
-
X
k
+
1
=
X
~
k
+
1
-
S, we conclude that
(27)
Δ
k
+
1
≈
1
2
Δ
k
-
1
2
S
Δ
k
S
.
Applying (27) recursively, and after some algebraic manipulations, we have
(28)
∥
Δ
k
+
1
∥
≤
1
2
k
+
1
∥
Δ
0
-
S
Δ
0
S
∥
.
From (28), we can conclude that the perturbation at the iterate k
+
1 is bounded. Therefore, the sequence {
X
k
} generated by (20) is asymptotically stable.The iteration (20) requires one matrix inversion per computing step and obtains both A
1
/
2 and A
-
1
/
2 which are of interest in most practical problems. The implementation of (20) for computing principal square roots requires a sharp attention so as to save much effort. Since the intermediate matrices are all sparse (at least half of the entries are zero), then one could simply use sparse approximation techniques to save up the memory and time.


### 2.1. A Novel Variant of Mid-Point Scheme

The analysis of the method from applying mid-point method toX
2
=
A was of limited interest, because it yields an iteration method that is not in general stable. The method with cubic convergence derived via the matrix sign function (20) is more useful. However, this method is somewhat not challenging, since it could be a special case of a family of iteration methods due to Kenney and Laub [15] and derived from the (explicitly known) rational Padé approximations for h
(
ζ
)
=
(
1
-
ζ
)
-
1
/
2. These methods are surveyed in Higham [1, Section 5.4] and the application to the matrix square root is discussed in Section 6.7. General results on the stability and rate of convergence for this family of methods are given in Theorems 6.11 and 6.12 of [1].To move forward, we recall that recently authors in [16] proposed an efficient and interesting variant of mid-point method for nonlinear equations and extended it for polar decomposition, while authors of [17] developed that variant for matrix sign function as follows:
(29)
X
k
+
1
=
(
I
+
18
Y
k
+
13
Z
k
)
[
X
k
(
7
I
+
Y
k
)
(
I
+
3
Y
k
)
]
-
1
,
where Y
k
=
X
k
X
k
,
Z
k
=
Y
k
Y
k.Theorem 5.
LetA
=
[
a
i
,
j
]
n
×
n have no eigenvalues on R
-. If one uses X
0
= (21), then the iterative method (29) converges to (22).Proof.
The proof of this theorem is similar to the proof of Theorem3.Lemma 6.
The sequence{
X
k
}
k
=
0
k
=
∞ generated by (29) using X
0
= (21) is asymptotically stable.Proof.
The proof of this lemma is similar to the proof of Lemma4.Note that in general the inverse of the nonsingular matrix(30)
[
0
K
H
0
]
is of the form
(31)
[
0
H
-
1
K
-
1
0
]
.
In fact, all matrices X
k
(
7
I
+
Y
k
)
(
I
+
3
Y
k
) in (29) would be of the form
(32)
[
0
†
‡
0
]
.
This itself might be helpful to reduce the computational time of matrix inversion problem.An implementation of (29) to compute A
1
/
2 and A
-
1
/
2 for a square matrix A having no eigenvalues on R
-, in the programming package Mathematica, is brought forward as follows:sqrtMatrix2[A_,maxIterations_,tolerance_]:= Module[{k=0},{n,n}=Dimensions[A];Id=SparseArray[{{i_,i_}  ->1.},{n,n}];A1=SparseArray@ArrayFlatten[{{0,A},{Id,0}}];Y[0]=A1;R[0]=1;Id2=SparseArray[{{i_,i_}  ->1.},{2n,2n}];While[k<maxIterations  &&  R[k]>=tolerance,Y2=SparseArray[Y[k].Y[k]];l1=SparseArray[Y[k].(7Id2+Y2).(Id2+3Y2)];l2=SparseArray@ArrayFlatten[{{0,Inverse@l1[[n+1;;2n,1;;n]]},{Inverse@l1[[1;;n,n+1;;2n]],0}}];Y[k+1]=SparseArray[(Id2+18Y2+13Y2.Y2).l2];R[k+1]=Norm[Y[k+1]-Y[k],Infinity]/Norm[Y[k+1],Infinity];k++];Y[k]];The above three-argument functionsqrtMatrix2 computes the principal matrix square root and its inverse at the same time by entering the three arguments “matrix A with no eigenvalue in R
-,” “the maximum number of iterations that (29) could take,” and the tolerance of the stopping termination in the infinity norm ∥
X
k
+
1
-
X
k
∥
∞
/
∥
X
k
+
1
∥
∞
≤
t
o
l
e
r
a
n
c
e. ClearlysqrtMatrix2[A,maxIterations,tolerance] [[1;;n,n+1;;2n]] and sqrtMatrix2[A,maxIterations, tolerance][[n+1;;2n,1;;n]] give A
1
/
2 and A
-
1
/
2, respectively, wherein n is the size of the input matrix A.


### 2.2.  Extension to Matrixpth Root

Computing thepth roots of matrices can arise in certain computations. The unique matrix X such that X
p
=
A, (p
∈
Z
+) and whose eigenvalues are in the segment
(33)
S
=
{
z
∈
C
∖
{
0
}
:
|
arg
(
z
)
|
<
π
p
}
is called the principal pth root and is denoted by A
1
/
p.Many authors have investigated methods for computingpth root of matrices. The methods are normally based on iteration or Schur normal form; see, for example, [18] and the references therein. In general, number of pth roots for a given matrix may be zero, finite, or infinite.In this subsection, we provide a combination of Algorithm 7.9 in [1] and our method (29) to compute principal matrix pth root in what follows.GivenA
∈
C
n
×
n having no eigenvalues on R
-, p
>
2, the following lines compute X
=
A
1
/
p using (29):(1)
computeA
1
/
2 by (29) and set B
=
A
1
/
2;
(2)
setC
=
B
/
θ, where θ is any upper bound for ρ
(
B
) (e.g., θ
=
∥
B
∥);
(3)
use the coupled Newton iteration (35) with X
0
=
I to compute
(34)
X
=
{
C
2
/
p
,
p
even
,
(
C
1
/
p
)
2
,
p
odd
;
(4)
X
←
∥
B
∥
2
/
p
X,where
(35)
X
k
+
1
=
X
k
(
(
p
-
1
)
I
+
M
k
p
)
,
X
0
=
I
,
M
k
+
1
=
(
(
p
-
1
)
I
+
M
k
p
)
-
p
M
k
,
M
0
=
C
,
when X
k
→
A
1
/
p, M
k
→
I.


## 2.1. A Novel Variant of Mid-Point Scheme

The analysis of the method from applying mid-point method toX
2
=
A was of limited interest, because it yields an iteration method that is not in general stable. The method with cubic convergence derived via the matrix sign function (20) is more useful. However, this method is somewhat not challenging, since it could be a special case of a family of iteration methods due to Kenney and Laub [15] and derived from the (explicitly known) rational Padé approximations for h
(
ζ
)
=
(
1
-
ζ
)
-
1
/
2. These methods are surveyed in Higham [1, Section 5.4] and the application to the matrix square root is discussed in Section 6.7. General results on the stability and rate of convergence for this family of methods are given in Theorems 6.11 and 6.12 of [1].To move forward, we recall that recently authors in [16] proposed an efficient and interesting variant of mid-point method for nonlinear equations and extended it for polar decomposition, while authors of [17] developed that variant for matrix sign function as follows:
(29)
X
k
+
1
=
(
I
+
18
Y
k
+
13
Z
k
)
[
X
k
(
7
I
+
Y
k
)
(
I
+
3
Y
k
)
]
-
1
,
where Y
k
=
X
k
X
k
,
Z
k
=
Y
k
Y
k.Theorem 5.
LetA
=
[
a
i
,
j
]
n
×
n have no eigenvalues on R
-. If one uses X
0
= (21), then the iterative method (29) converges to (22).Proof.
The proof of this theorem is similar to the proof of Theorem3.Lemma 6.
The sequence{
X
k
}
k
=
0
k
=
∞ generated by (29) using X
0
= (21) is asymptotically stable.Proof.
The proof of this lemma is similar to the proof of Lemma4.Note that in general the inverse of the nonsingular matrix(30)
[
0
K
H
0
]
is of the form
(31)
[
0
H
-
1
K
-
1
0
]
.
In fact, all matrices X
k
(
7
I
+
Y
k
)
(
I
+
3
Y
k
) in (29) would be of the form
(32)
[
0
†
‡
0
]
.
This itself might be helpful to reduce the computational time of matrix inversion problem.An implementation of (29) to compute A
1
/
2 and A
-
1
/
2 for a square matrix A having no eigenvalues on R
-, in the programming package Mathematica, is brought forward as follows:sqrtMatrix2[A_,maxIterations_,tolerance_]:= Module[{k=0},{n,n}=Dimensions[A];Id=SparseArray[{{i_,i_}  ->1.},{n,n}];A1=SparseArray@ArrayFlatten[{{0,A},{Id,0}}];Y[0]=A1;R[0]=1;Id2=SparseArray[{{i_,i_}  ->1.},{2n,2n}];While[k<maxIterations  &&  R[k]>=tolerance,Y2=SparseArray[Y[k].Y[k]];l1=SparseArray[Y[k].(7Id2+Y2).(Id2+3Y2)];l2=SparseArray@ArrayFlatten[{{0,Inverse@l1[[n+1;;2n,1;;n]]},{Inverse@l1[[1;;n,n+1;;2n]],0}}];Y[k+1]=SparseArray[(Id2+18Y2+13Y2.Y2).l2];R[k+1]=Norm[Y[k+1]-Y[k],Infinity]/Norm[Y[k+1],Infinity];k++];Y[k]];The above three-argument functionsqrtMatrix2 computes the principal matrix square root and its inverse at the same time by entering the three arguments “matrix A with no eigenvalue in R
-,” “the maximum number of iterations that (29) could take,” and the tolerance of the stopping termination in the infinity norm ∥
X
k
+
1
-
X
k
∥
∞
/
∥
X
k
+
1
∥
∞
≤
t
o
l
e
r
a
n
c
e. ClearlysqrtMatrix2[A,maxIterations,tolerance] [[1;;n,n+1;;2n]] and sqrtMatrix2[A,maxIterations, tolerance][[n+1;;2n,1;;n]] give A
1
/
2 and A
-
1
/
2, respectively, wherein n is the size of the input matrix A.


## 2.2.  Extension to Matrixpth Root

Computing thepth roots of matrices can arise in certain computations. The unique matrix X such that X
p
=
A, (p
∈
Z
+) and whose eigenvalues are in the segment
(33)
S
=
{
z
∈
C
∖
{
0
}
:
|
arg
(
z
)
|
<
π
p
}
is called the principal pth root and is denoted by A
1
/
p.Many authors have investigated methods for computingpth root of matrices. The methods are normally based on iteration or Schur normal form; see, for example, [18] and the references therein. In general, number of pth roots for a given matrix may be zero, finite, or infinite.In this subsection, we provide a combination of Algorithm 7.9 in [1] and our method (29) to compute principal matrix pth root in what follows.GivenA
∈
C
n
×
n having no eigenvalues on R
-, p
>
2, the following lines compute X
=
A
1
/
p using (29):(1)
computeA
1
/
2 by (29) and set B
=
A
1
/
2;
(2)
setC
=
B
/
θ, where θ is any upper bound for ρ
(
B
) (e.g., θ
=
∥
B
∥);
(3)
use the coupled Newton iteration (35) with X
0
=
I to compute
(34)
X
=
{
C
2
/
p
,
p
even
,
(
C
1
/
p
)
2
,
p
odd
;
(4)
X
←
∥
B
∥
2
/
p
X,where
(35)
X
k
+
1
=
X
k
(
(
p
-
1
)
I
+
M
k
p
)
,
X
0
=
I
,
M
k
+
1
=
(
(
p
-
1
)
I
+
M
k
p
)
-
p
M
k
,
M
0
=
C
,
when X
k
→
A
1
/
p, M
k
→
I.


## 3. Applications

We test the contributed method (29) denoted by PM1 using Mathematica 8 in machine precision [19]. Apart from these schemes, several iterative methods (3), (4), (20) denoted by PM and the method of Meini originally developed in [20] based on the cyclic reduction algorithm (CR) have been considered as follows:
(36)
Z
0
=
2
(
I
+
A
)
,
Y
0
=
I
-
A
,
Y
k
+
1
=
-
Y
k
Z
k
-
1
Y
k
,
Z
k
+
1
=
Z
k
-
2
Y
k
Z
k
-
1
Y
k
,
k
=
0,1
,
…
.
This method generates the sequence {
Z
k
}, which converges to 4
A
1
/
2. Note that (36) only computes A
1
/
2 while (4), (20), and (29) produce A
1
/
2 and A
-
1
/
2.Note that PM is equal to Halley’s third order method [14].The computational order of convergence for computing the principal matrix square root by matrix iterations can be estimated by(37)
ρ
=
ln
⁡
(
∥
X
k
2
-
A
∥
-
1
∥
X
k
+
1
2
-
A
∥
)
ln
⁡
(
∥
X
k
-
1
2
-
A
∥
-
1
∥
X
k
2
-
A
∥
)
,
wherein X
k
-
1, X
k, X
k
+
1 are the last three approximations for finding A
1
/
2 in the convergence phase.Example 7.
As the first example, we consider a matrix(38)
A
=
(
1
4
+
π
2
81
π
9
π
9
1
)
,
whereas its eigenvalues are {
1.15549,0.216359
}. The principal matrix square root of A is
(39)
X
=
A
1
/
2
≈
(
0.566105
0.226654
0.226654
0.973975
)
.We compare the behavior of different methods and report the numerical results usingl
∞ for all norms involved with the stopping criterion
(40)
E
k
+
1
=
∥
X
k
+
1
-
X
k
∥
∞
∥
X
k
+
1
∥
∞
≤
ϵ
=
1
0
-
8
,
in Table 1. The numerical results are in harmony with the theoretical aspects of Section 2.Table 1
Results of comparisons for Example7.
Newton
DB
CR
PM
PM1
IT
6
6
5
4
3
E
k
+
1
9.33948
×
1
0
-
15
9.31636
×
1
0
-
15
9.22389
×
1
0
-
15
3.29612
×
1
0
-
12
7.94163
×
1
0
-
11
ρ
1.99877
1.99908
2.00031
2.97414
4.02282A differential equation is a mathematical equation for an unknown function of one or several variables that relate the values of the function itself and its derivatives of various orders. A matrix differential equation contains more than one function stacked into vector form with a matrix relating the functions to their derivatives. We now illustrate the application of the new solver for solving matrix differential equations.Example 8.
In this test, we compare the behavior of different methods for solving the following matrix differential equation:(41)
∂
2
y
(
t
)
∂
t
2
+
A
y
(
t
)
=
0
,
y
(
0
)
=
y
0
y
′
(
0
)
=
y
0
′
(
0
)
,
where its exact solution can be expressed as
(42)
y
(
t
)
=
cos
⁡
(
A
1
/
2
t
)
y
0
(
0
)
+
A
-
1
/
2
sin
(
A
1
/
2
t
)
y
0
′
(
0
)
,
and the sparse SPD 100
×
100 matrix A in Mathematica is defined as follows:
A=SparseArray[{{i_,i_}  ->10.,{i_,j_}  /;
Abs[i-j]==1  ->-5.},{n,n}];The interesting point in this example is that we need bothA
1
/
2 and A
-
1
/
2 at the same time and this paves up the application of (29). We compare the behavior of different methods in Figure 1 for finding the principal matrix square root of A. The results are promising for the proposed method in terms of the number of iterations. Figure 1 reveals that the fastest method is (29), while the Newton’s matrix method is numerically unstable and this makes it diverge after the seventh iterate.Figure 1
The comparison of different methods for findingA
1
/
2 in Example 8.Note again that among the compared methods only the method of DB and our proposed method PM1 are able to findA
-
1
/
2 at the same time which is fruitful in solving the matrix differential equation (41). Due to page limitation, we could not present the final solution of (41) explicitly here. Clearly, using the obtained principal matrix square root and the principal inverse matrix square root produces the final solution.


## 4. Conclusion

In this paper, we have developed and discussed the mid-point iterative root-finding method for solving some special and important matrix problems. We discussed under what conditions the obtained method is convergent to the principal matrix square root, that is, whenA has no nonpositive real eigenvalues.Furthermore, an asymptotically stable variant of the mid-point method using an identity between the matrix square root and the matrix sign function has been derived. The numerical experiments which we have performed, reported at the end of the paper, show that the presented procedure is a reliable tool for computing the principal matrix square root and its inverse.



---

*Source: 525087-2014-08-27.xml*