# An Algorithm for Computing Geometric Mean of Two Hermitian Positive Definite Matrices via Matrix Sign

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

---

## Abstract

Using the relation between a principal matrix square root and its inverse with the geometric mean, we
present a fast algorithm for computing the geometric mean of two Hermitian positive definite matrices. The algorithm
is stable and possesses a high convergence order. Some experiments are included to support the proposed computational
algorithm.


---

## Body

## 1. Introduction

This paper tries to provide a fast algorithm for finding the geometric mean of two Hermitian positive definite matrices via the application of matrix sign function. It is known that the scalar arithmetic-geometric meanagm
(
a
,
b
) of two (nonnegative) numbers a and b is defined by starting with a
0
=
a and b
0
=
b and then iterating
(1)
a
k
+
1
=
1
2
[
a
k
+
b
k
]
,
b
k
+
1
=
a
k
b
k
,
k
=
0,1
,
2
,
…
,
until a
k
=
b
k to the desired precision. The scalar sequences of {
a
k
}, {
b
k
} converge to each other. Note that in (1), a
k
b
k is the geometric mean of two positive numbers a
k and b
k per computing cycle. Although this iterative formula is quite easy and reliable for two scalars, its extension for general square and nonsingular matrices is not an easy task. In this work, we are concerned with the matrix case of two square Hermitian positive definite matrices.A right definition of the matrix geometric meanGM
(
A
,
B
) of two positive definite matrices A and B can be expressed by
(2)
GM
(
A
,
B
)
∶
=
A
(
A
-
1
B
)
1
/
2
,
where given a square matrix M having no nonpositive real eigenvalues, M
1
/
2 denotes the unique solution of the quadratic matrix equation
(3)
X
2
=
M
,
whose eigenvalues lie in the right half plane. The definition (2) was given in the seventies by Pusz and Woronowicz [1]. There are some other important definitions for computing the matrix geometric mean of two matrices; see for example, the work of Lawson-Lim [2]. Note that when just two matrices are involved the theory is well developed, but in case of finding the matrix geometric mean of more than two matrices, the formulation is kind of hard; see for more [3].A variant formulation for the scalar case of geometric mean could be expressed by(4)
gm
(
a
,
b
)
=
a
b
=
(
a
b
a
a
)
1
/
2
=
a
1
/
2
(
b
a
)
1
/
2
a
1
/
2
=
a
1
/
2
(
a
-
1
/
2
b
a
-
1
/
2
)
1
/
2
a
1
/
2
.As could be seen in (4), the computation of scalar geometric mean is fully related to the square root of the positive scalar a and its inverse. A similar and significant definition for a matrix geometric mean has been developed and suggested by Bhatia in [4] as follows:
(5)
A
#
B
=
A
1
/
2
(
A
-
1
/
2
B
A
-
1
/
2
)
1
/
2
A
1
/
2
,
for two Hermitian positive definite (HPD) matrices A and B. We use the notation # to show the geometric mean of two HPD matrices.Note that the formulas (2) and (5) coincide and are equal. For more, refer to [5].The rest of this paper is organized as follows. In Section2 we combine the root-finding Newton’s method and a Chebyshev-Halley-type method to devise a fast iterative formula for solving nonlinear equations. In Section 3 we provide the link between nonlinear equation solvers and the calculation of some special matrix functions. This will illustrate how the new algorithms could be constructed and implemented. An implementation of the proposed iterative formula in symbolic software Mathematica [6] along a discussion about the stability of the scheme will be provided therein. Note that the idea of computing the geometric mean using the sign function can also be found in [7]. In Section 4 we show the numerical results and highlight the benefit of the proposed technique. Conclusions are drawn in Section 5.


## 2. Construction of a New Method

It is known that a common way for improving the convergence speed of iterative methods for solving nonlinear scalar smooth equations is to combine the already developed schemes. Hence, let us first combine the well-known method of Newton into a special method of Chebyshev-Halley-type scheme [8] to derive a new iterative scheme as follows:
(6)
y
k
=
x
k
-
f
(
x
k
)
f
′
(
x
k
)
,
x
k
+
1
=
y
k
-
(
1
+
1
2
(
L
k
1
-
(
(
4
/
5
)
L
k
)
)
)
f
(
y
k
)
f
′
(
y
k
)
,
wherein L
k
=
f
′′
(
y
k
)
f
(
y
k
)
/
f
′
(
y
k
)
2. For obtaining a background about nonlinear equations, one may consult [9, 10].Theorem 1.
Letα
∈
D be a simple zero of a sufficiently differentiable function f
:
D
⊆
R
→
R for an open interval D, which contains x
0 as an initial approximation of α. Then, the iterative expression (6) without memory satisfies the error equation below
(7)
e
k
+
1
=
(
2
c
2
5
5
-
c
2
3
c
3
)
e
k
6
+
O
(
e
k
7
)
,
wherein c
j
=
f
(
j
)
(
α
)
/
(
j
!
f
′
(
α
)
), e
k
=
x
k
-
α.Proof.
The proof of this theorem is based on Taylor’s series expansion of the iterative method (6) around the solution in the kth iterate. To save the space and not to be distracted from the main topic, we here exclude the proof. The steps of the proof are similar to those taken in [11].The iterative method (6) reaches sixth-order convergence using five functional evaluations and thus achieves the efficiency index 6
1
/
5
≈
1.430, which is higher than that of Newton; that is, 2
1
/
2
≈
1.414. Furthermore, applying (6) for solving the polynomial equation g
(
x
)
≡
x
2
-
1
=
0 has global convergence (except the points lying on the imaginary axis). This global behavior which could easily be shown by drawing the basins of attraction (6) for g
(
x
) is useful in practical matrix problems so as to allow us deal with all kinds of HPD matrices with any spectra.Remark 2.
Note that since the global behavior can be observed from the basins of attraction, then we do not pursue this fact by theoretical analysis for (6).A recent discussion about the relation between matrix sign and nonlinear equation solvers is given in [12, 13]. We state that all of such extensions in the scalar case and studying their orders are of symbolic computational nature. Such an extension to the matrix environment will also be done in the next section using symbolic computations.


## 3. Construction of an Algorithm for Geometric Mean

The relation between Section2 and our main aim for computing matrix geometric mean is not clear at the first sight. To illustrate this, we recall that many of the important matrix functions such as matrix square root, that is, the solution to the matrix equation (3), and matrix sign function which is the solution to the following matrix equation:
(8)
X
2
=
I
,
can be calculated approximately by matrix iteration functions. Such fixed-point type methods are convergent under special conditions to the aimed matrix function and are basically originated from root-finding methods [7].On the other hand, the important definition of the matrix geometric mean (5) requires the computation of the matrix square root of A and its inverse. Hence, in order to design an efficient algorithm for (5), we wish to construct an iterative expression in which we compute both A
1
/
2 and A
-
1
/
2 at the same time.Toward this aim, we apply the new nonlinear equation solver (6) for solving the matrix equation (8). This application would yield the following matrix iteration in its reciprocal form:
(9)
X
k
+
1
=
(
48
X
k
+
272
X
k
3
+
272
X
k
5
+
48
X
k
7
)
×
[
7
I
+
148
X
k
2
+
330
X
k
4
+
148
X
k
6
+
7
X
k
8
]
-
1
,
with a proper X
0 for finding the sign matrix.Remark 3.
The iteration (9) is not a member of the Padé family of iterations introduced in [14] for matrix sign.In order to connect (9) with our aim, we remind an important identity as follows [7]:
(10)
sign
⁡
(
[
0
A
I
0
]
)
=
[
0
A
1
/
2
A
-
1
/
2
0
]
,
which indicates an important relationship between principal matrix square root A
1
/
2 and the matrix sign function.The identity (10) has an advantage which is the computation of the principal inverse matrix square root along with the principal matrix square root at the same time. Let us in what follows first study the stability behavior of (9).Lemma 4.
LetA
=
[
a
i
,
j
]
n
×
n have no eigenvalues on R
-. Then, the sequence {
X
k
}
k
=
0
k
=
∞ generated by (9) using X
0
=
[
0
A
I
0
] is asymptotically stable.Proof.
LetΔ
k be a numerical perturbation introduced at the kth iterate of (9). We obtain
(11)
X
~
k
=
X
k
+
Δ
k
.
All terms that are quadratic in their errors are removed from our analysis. This formal manipulation is meaningful if Δ
k is sufficiently small. We have
(12)
X
~
k
+
1
=
(
48
X
~
k
+
272
X
~
k
3
+
272
X
~
k
5
+
48
X
~
k
7
)
×
[
7
I
+
148
X
~
k
2
+
330
X
~
k
4
+
148
X
~
k
6
+
7
X
~
k
8
]
-
1
=
(
48
(
X
k
+
Δ
k
)
+
272
(
X
k
+
Δ
k
)
3
+
272
(
X
k
+
Δ
k
)
5
+
48
(
X
k
+
Δ
k
)
7
)
×
[
7
I
+
148
(
X
k
+
Δ
k
)
2
+
330
(
X
k
+
Δ
k
)
4
+
148
(
X
k
+
Δ
k
)
6
+
7
(
X
k
+
Δ
k
)
8
]
-
1
.
Note that the commutativity between X
k and Δ
k is not used throughout this lemma. To simplify this, the following identity will be used (for any nonsingular matrix B and the matrix C):
(13)
(
B
+
C
)
-
1
≈
B
-
1
-
B
-
1
C
B
-
1
.
Simplifying yields
(14)
X
~
k
+
1
≈
(
S
+
5
2
Δ
k
+
3
2
S
Δ
k
S
)
(
I
+
2
S
Δ
k
+
2
Δ
k
S
)
-
1
≈
(
S
+
5
2
Δ
k
+
3
2
S
Δ
k
S
)
(
I
-
2
S
Δ
k
-
2
Δ
k
S
)
,
where S
2
=
I, S
-
1
=
S, and for enough large k, we assumed X
k
≈
S. After some algebraic manipulations and using Δ
k
+
1
=
X
~
k
+
1
-
X
k
+
1
≈
X
~
k
+
1
-
S, we conclude that
(15)
Δ
k
+
1
≈
1
2
Δ
k
-
1
2
S
Δ
k
S
.
Applying (15) recursively, we have
(16)
∥
Δ
k
+
1
∥
≤
1
2
k
+
1
∥
Δ
0
-
S
Δ
0
S
∥
.
From (16), we can conclude that the perturbation at the iterate k
+
1 is bounded. This allows to conclude that a perturbation at a given step is bounded at the subsequent steps. Therefore, the sequence {
X
k
} generated by (9) is asymptotically stable.Note that Lemma4 just shows the asymptotical stability. In general, the fixed-point iteration X
k
+
1
=
g
(
X
k
) is stable in a neighborhood of a fixed point X if Fréchet derivative d
g
X has bounded powers [7]. Furthermore, if X
k
+
1
=
g
(
X
k
) is any super-linearly convergent iteration for S
=
sign
⁡
(
X
0
), then d
g
S
(
E
)
=
L
S
(
E
)
=
(
1
/
2
)
(
E
-
S
E
S
), where L
S is the Fréchet derivative of the matrix sign function at S. Hence d
g
S is idempotent (d
g
S
∘
d
g
S
=
d
g
S) and the iteration is stable, thus all sign iterations, such as (9), are automatically stable.It is now easy to deduce the following convergence theorem for (9).Theorem 5.
LetA
=
[
a
i
,
j
]
n
×
n have no eigenvalues on R
-. If X
0
=
[
0
A
I
0
], then the iterative method (9) converges to S
=
sign
⁡
(
[
0
A
I
0
]
).Proof.
The convergence of rational iterations can be analyzed in terms of the convergence of the eigenvalues of the matricesX
k. The reason for this is that if X has a Jordan decomposition X
=
Z
J
Z
-
1, then R
(
X
)
=
Z
R
(
J
)
Z
-
1. Let A have a Jordan canonical form arranged as Z
-
1
A
Z
=
Λ, where Z is a nonsingular matrix. It is also known that [7]
(17)
sign
⁡
(
Λ
)
=
sign
⁡
(
Z
-
1
A
Z
)
=
Z
-
1
sign
⁡
(
A
)
Z
.
If we define D
k
=
Z
-
1
X
k
Z, then from the method (9), we obtain
(18)
D
k
+
1
=
(
48
D
k
+
272
D
k
3
+
272
D
k
5
+
48
D
k
7
)
×
[
7
I
+
148
D
k
2
+
330
D
k
4
+
148
D
k
6
+
7
D
k
8
]
-
1
.
Notice that if D
0 is a diagonal matrix then, based on an inductive proof, all successive D
k are diagonal too. From (18), it is enough to prove that {
D
k
} converges to sign
⁡
(
Λ
), in order to ensure the convergence of the sequence generated by (9).
We can write (18) as n uncoupled scalar iterations to solve g
(
x
)
=
x
2
-
1
=
0, given by
(19)
d
k
+
1
i
=
48
d
k
i
+
272
d
k
i
3
+
272
d
k
i
5
+
48
d
k
i
7
7
+
148
d
k
i
2
+
330
d
k
i
4
+
148
d
k
i
6
+
7
d
k
i
8
,
where d
k
i
=
(
D
k
)
i
,
i and 1
≤
i
≤
n. From (18) and (19), it is enough to study the convergence of {
d
k
i
} to sign
⁡
(
λ
i
), for all 1
≤
i
≤
n.
From (19) and since the eigenvalues of A are not pure imaginary, we have that sign
⁡
(
λ
i
)
=
s
i
=
±
1. Thus, we attain
(20)
d
k
+
1
i
-
1
d
k
+
1
i
+
1
=
-
(
-
1
+
d
k
i
)
6
(
7
+
d
k
i
(
-
6
+
7
d
k
i
)
)
(
1
+
d
k
i
)
6
(
7
+
d
k
i
(
6
+
7
d
k
i
)
)
.
Since |
d
0
i
|
=
|
λ
i
|
>
0, we have
(21)
lim
⁡
k
→
∞
|
d
k
+
1
i
-
1
d
k
+
1
i
+
1
|
=
0
,
and lim
⁡
k
→
∞
|
d
k
i
|
=
1
=
|
sign
⁡
(
λ
i
)
|. This shows that {
d
k
i
} is convergent. Now, it could be easy to conclude that lim
⁡
k
→
∞
D
k
=
sign
⁡
(
Λ
). Finally, we have
(22)
lim
⁡
k
→
∞
X
k
=
Z
(
lim
⁡
k
→
∞
D
k
)
Z
-
1
=
Z
sign
⁡
(
Λ
)
Z
-
1
=
sign
⁡
(
A
)
.
The proof is complete.The iteration (9) requires one matrix inversion before computing step and obtains both A
1
/
2 and A
-
1
/
2 which are of interest in (5). The implementation of (9) 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.An implementation of (9) to compute A
#
B for two HPD matrices in the programming package Mathematica is brought forward as in (Algorithm 1).Algorithm 1
C
learAll["Global’*"]
FunM[fun_, X_]:= Module[{faux, dim, mataux, JordanD, sim, JordanF, eps,
fdiag, diagQ, fauxD}, (dim = Length@X;
faux[xx_, i_, j_]:=
Which[i <= j, 1/Abs[i - j]! (D[fun,{x, Abs[i - j]}]) /. x -> xx, True, 0];
mataux[Y_]:= Table[faux[Y[[i, j]], i, j],{i, 1, dim},{j, 1, dim}];
JordanD = JordanDecomposition[X]
/
/ N; sim = JordanD[[
1
]];
JordanF = JordanD[[
2
]]; eps = 1*10
⋀
-
10;
diagQ = Norm[JordanF
-
DiagonalMatrix[Diagonal[JordanF]]];
fauxD[xx_]:= (fun)
/
. x -> xx;
fdiag:= DiagonalMatrix[Map[fauxD, Diagonal[JordanF]]];
Which[diagQ < eps, sim.fdiag.Inverse[sim], True,
sim.mataux[JordanF].Inverse[sim]])]
MGM[A_, B_, maxIterations_, tolerance_]:= Module[{k = 0},
{n, n}
= Dimensions[A]; Id = SparseArray[
{
{
i_, i_}
-> 1.},{n, n}];
A1 = SparseArray@ArrayFlatten[
{
{
0, N@A},{Id, 0
}
}
]; Y[
0
] = A1;
R[
0
] = 1; Id2 = SparseArray[
{
{
i_, i_}
-> 1.},{2 n, 2 n}];
{Quiet@While[k < maxIterations
&&
R[k] >= tolerance,
Y2 = Y[k].Y[k]; Y3 = Y2.Y[k]; Y4 = Y3.Y[k];
Y5 = Y4.Y[k]; Y6 = Y5.Y[k]; Y7 = Y6.Y[k]; Y8 = Y7.Y[k];
l1 = SparseArray[7 Id2 + 148 Y2 + 330 Y4 + 148 Y6 + 7 Y8];
l2 = SparseArray@ArrayFlatten[
{
{
Inverse@l1[[1;; n, 1;; n]], 0},
{0, Inverse@l1[[n + 1;; 2 n, n + 1;; 2 n]]
}
}
];
Y[k + 1] = SparseArray[(48 Y[k] + 272 Y3 + 272 Y5 + 48 Y7).l2];
R[k + 1] = Norm[Y[k + 1] - Y[k], Infinity]/ Norm[Y[k + 1], Infinity]; k++];
AS = Y[k][[1;; n, n + 1;; 2 n]]; IAS = Y[k][[n + 1;; 2 n, 1;; n]];
Mat = (IAS.B.IAS)};
AS.FunM[Sqrt[x], Mat].AS]In Algorithm1, the four-argument function MGM computes A
#
B by entering the four arguments “matrix A,” “matrix B,” “the maximum number of iterations that (9) is allowed to take,” and the “tolerance” of the stopping termination in the infinity norm ∥
X
k
+
1
-
X
k
∥
∞
/
∥
X
k
+
1
∥
∞
≤
tolerance.Note that for computing the principal matrix square root of(
A
-
1
/
2
B
A
-
1
/
2
)
1
/
2, we have used the Jordan Canonical approach, which has been provided in the general form for computing matrix functions in the code FunM.


## 4. Experiments

We test the contributed method (9) denoted by PM using Mathematica 8 in machine precision. Apart from this scheme, another iterative method which is known as DB method [15] and given by
(23)
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
]
is considered. This method generates the sequences {
Y
k
} and {
Z
k
} which converge to A
1
/
2 and A
-
1
/
2, respectively.We remark that the first and the second substeps of (6) result in the quadratic Newton’s method (NM) and cubic Chebyshev-Halley’s method (CHM) for matrix sign [7] in what follows:
(24)
X
k
+
1
=
1
2
(
X
k
+
X
k
-
1
)
,
X
k
+
1
=
(
16
X
k
+
24
X
k
3
)
[
3
I
+
30
X
k
2
+
7
X
k
4
]
-
1
.Example 6.
As the first example, we consider the matrixA
=
(
2
1
1
2
), B
=
(
38
1
1
2
), whereas its exact matrix geometric mean is given by [16]: A
#
B
=
(
8
1
1
2
).The proposed approach converges to the solution matrix in 2 iterations, which shows a completely fast convergence.Example 7.
We now consider the two HPD matrices as follows:(25)
A
=
(
2
1
1
⋱
⋱
⋱
⋱
1
1
2
)
n
×
n
,
B
=
(
3
2
2
5
2
5
⋱
⋱
⋱
⋱
2
5
2
5
3
2
)
n
×
n
,
when n
=
250 and compute A
#
B.We compare the behavior of different methods and report the numerical results usingl
∞ for all norms involved with the stopping criterion ∥
X
k
+
1
-
X
k
∥
∞
/
∥
X
k
+
1
∥
∞
≤
ϵ
=
1
0
-
5 in Figure 1. As could be seen, the numerical results are in harmony with the theoretical aspects of Section 3 and show a fast convergence for the proposed method (9) instead of DB, NM, and CHM.Figure 1
The log plot for the comparison of number of iterations in Example7.


## 5. Conclusions

Based on the quadratical Newton’s scheme and the cubical method of Chebyshev-Halley, we have developed an iterative method with sixth order of convergence for solving nonlinear scalar smooth equations. The computational efficiency showed its superiority in contrast to NM. Then, the method with global behavior for finding matrix sign function has been extended for computing the principal matrix square root and its inverse. This procedure was followed so as to build a fast algorithm for finding the matrix geometric mean of two HPD matrices. We also have studied the asymptotical stability for the proposed technique.To illustrate the new technique some numerical examples were presented. Computational results have justified robust and efficient convergence behavior of the proposed method. Similar numerical experimentations, carried out for a number of problems of different types, confirmed the above conclusions to a great extent.We conclude the paper with the remark that in many numerical applications high precision in computations is required. The results of numerical experiments show that the high order efficient methods such as (9) associated with a multiple-precision arithmetic are very useful, because they yield a clear reduction in the number of iterates.



---

*Source: 978629-2014-08-06.xml*