# The arithmetic-geometric mean for symmetric positive definite matrices

A while back, I wanted to see if the notion of the arithmetic-geometric mean could be extended to a pair of symmetric positive definite matrices. (I considered positive definite matrices only since the notion of the matrix square root is a bit intricate for other kinds of matrices.)

I expected that some complications would arise since, unlike scalar multiplication, matrix multiplication is noncommutative. Another complication would be that the product of two symmetric matrices need not be symmetric (though the positive definiteness is retained, so one can still speak of a principal matrix square root).

By analogy with the scalar AGM, I considered the iteration

$$\mathbf A_0=\mathbf A \; ; \; \mathbf B_0=\mathbf B$$ $$\mathbf A_{i+1}=\frac12(\mathbf A_i+\mathbf B_i) \; ; \; \mathbf B_{i+1}=\sqrt{\mathbf A_i \mathbf B_i}$$

I cranked up a short Mathematica routine:

matAGM[u_, v_] := First[FixedPoint[
{Apply[Plus, #]/2, MatrixPower[Apply[Dot, #], 1/2]} &, {u, v}]] /;
MatrixQ[u, InexactNumberQ] && MatrixQ[v, InexactNumberQ]


and decided to try it out on randomly generated SPD matrices.

(A numerical note: Mathematica uses the numerically stable Schur decomposition in computing matrix functions like the matrix square root.)

I found that for all of the randomly generated pairs of SPD matrices I tried, the process was convergent (though the rate of convergence is apparently not as fast as the scalar AGM). As expected, the order matters: matAGM[A, B] and matAGM[B, A] are usually not equal (and both results are unsymmetric) unless A and B commute (for the special case of diagonal A and B, the result is the diagonal matrix whose entries are the arithmetic-geometric means of the corresponding entries of the pair.)

I now have three questions:

1. How do I prove or disprove that this process converges for any pair of SPD matrices? If it is convergent, what is the rate of convergence?

2. Is there any relationship between matAGM[A, B] and matAGM[B, A] if the two matrices A and B do not commute?

3. Is there any relationship between this matrix arithmetic-geometric mean and the usual scalar arithmetic-geometric mean? Would, say, arithmetic-geometric means of the eigenvalues of the two matrices have anything to do with this?

More digging around has me convinced that I should indeed be considering the formulation of the geometric mean by Pusz and Woronowicz:

$$\mathbf A\#\mathbf B=\mathbf A^{1/2}(\mathbf A^{-1/2}\mathbf B\mathbf A^{-1/2})^{1/2}\mathbf A^{1/2}$$

as more natural; the proof of convergence is then simplified, as shown in the article Willie linked to. However, I'm still wondering why the original "unnatural" formulation seems to be convergent (or else, I'd like to see a pair of SPD matrices that cause trouble for the unnatural iteration). I am also interested in how elliptic integrals might crop up in here, just as they did for the scalar version of the AGM.

-

Nice question!

What follows below is not really an answer, but too long to fit as a comment. It would be helpful if you could write your Mathematica algorithm is traditional mathematical notation so it is easier to see what is going on.

First, I am sure you are aware of the standard definitions of matrix means. I recall them below for the benefit of others who consider your question.

1. The arithmetic mean is: $(A+B)/2$ as expected. It satisfies several of the standard desirable properties of a mean of spd matrices (nonnegative, if $A \le B$, then $A \le AM(A,B) \le B$, etc)

2. The standard geometric mean is: $$GM(A,B) = A^{1/2}(A^{-1/2}BA^{-1/2})^{1/2}A^{1/2}$$ This particular choice of GM satisfies: $GM(A,B)=GM(B,A)$.

Now if you iterate the AGM, using the above two definitions of the means, then I think convergence, etc. can be proven---though you might have to resort to some kind of Fixed-point theorem. If I get time, I will think more about this.

-
Hmm, I didn't use that definition of the geometric mean, but that would make for an interesting variation... :) as I said, it's been a while since I tried this experiment, but I'll look into this when I get to my computer with Mathematica. – J. M. Dec 23 '10 at 0:13
@J.M. I think Slowsolver's version of the GM is perhaps more natural. That GM can be defined as the "geodesic half way point" between $A$ and $B$ by putting a Riemannian metric on the space of SPD matrices $$g(X,Y) = \mathop{tr}(A^{-1}XA^{-1}Y)$$ for $X,Y$ in the tangent space at $A$ (so they are arbitrary symmetric matrices). This metric is formally similar to the $dx\otimes dx / x^2$ metric on $\mathbb{R}_+$ with regards to which the usual geometric mean on $\mathbb{R}_+$ is a "geodesic midpoint". – Willie Wong Dec 23 '10 at 0:50
Thanks @Willie for convincing me further to try out that version of the GM; I'll look into it when I'm on my Mathematica machine again. – J. M. Dec 23 '10 at 0:57

I'll use the characterisation of the geometric mean as given by Slowsolver. As I mentioned in a comment to his answer, the geometric mean $GM(A,B)$ is the geodesic midpoint corresponding to the Riemannian metric on the space of symmetric positive definite matrices with

$$g(X,Y) = \mathop{tr}(A^{-1}XA^{-1}Y)~,\quad X,Y\in T_A(SPD)$$

where the tangent space $T_A(SPD)$ is the space of symmetric matrices. This interpretation has the following advantage: both the AM and GM are "translation invariant". Let $P$ be a non-singular matrix, then $AM(PAP^T,PBP^T) = P(AM(A,B))P^T$ trivially, and similarly for $GM$ using that $$g_A(X,Y) = g_{PAP^T}(PXP^T, PYP^T)$$ so that a geodesic $\gamma$ joining $A,B$ becomes a geodesic $P\gamma P^T$ joining $PAP^T$ and $PBP^T$. By first conjugating with $P = A^{-1/2}$, we can WLOG assume $A = I$. By then conjugating with a suitable orthogonal matrix $O$ to diagonalize $B$, we can WLOG assume $A = I$ and $B$ is diagonal.

The result on convergence now follows from the convergence of the AGM for the single, real variable case.

1. The rate of convergence is the same as the scalar case, up to a constant factor depending on the initial values of $A,B$ (which determines $P$).
2. In this definition $GM(A,B) = GM(B,A)$, so the question is moot.
3. This is in some sense the natural generalisation of the scalar case. Notice that the scalar geometric mean corresponds to the geodesic midpoint for the multiplicative Lie group $(\mathbb{R}_+,\times)$ with the invariant metric. A first generalisation is to the abelian Lie group $\mathbb{R}^k_+$ with componentwise multiplication. (This corresponds to the diagonal case.) What we showed above is that the entire AGM sequence for arbitrary starting SPD matrices $A$ and $B$ is contained, up to conjugation by $P\cdot P^T$, in one such abelian case. Note that this is slightly less related to eigenvalues per se, since the transformation we use above does not preserve eigenvalues.

EDIT: It just occurred to me to do a literature search. In Lawson, J. D. & Lim, Y. "The geometric mean, matrices, metrics, and more". Amer. Math. Monthly, 2001, 108, 797-812, many properties of this geometric mean on SPD matrices are derived in Sections 2 and 3. In particular, the crucial property (as slowsolver remarks below) of being able to set $A = I$ and $B$ diagonal, is the content of Lemma 3.1 of that paper. That is followed by the harmonic-geometric-arithmetic-mean inequality for SPD matrices. By starting with that and playing with monotonicity properties, you presumably also get another proof of the convergence result.

-
Interesting, thanks for the analysis Willie! (I still have to do numerical experiments though.) I still can't help but wonder if this is related to the usual scalar AGM. (I have a feeling there is, but the connection is not transparent to me.) – J. M. Dec 23 '10 at 2:18
@J.M. I'm not quite sure what you mean by related. But at the very least you can consider the usual scalar AGM on positive numbers as a special case of the above, by consider the scalar numbers as $1\times 1$ matrices. – Willie Wong Dec 23 '10 at 2:28
@Willie: Admittedly, I could just be babbling... but it would be interesting if this construction turns out to be equivalent to evaluating elliptic integrals for SPD matrix arguments (i.e., the result is similar to a diagonal matrix whose elements are elliptic integrals). – J. M. Dec 23 '10 at 2:31
@J.M. oh.. that is a bit beyond my familiarity. My mental picture of the above is more like $GL(n)$ bundle over $\mathbb{R}_+^n$ (but not quite). Sorry, can't help there. – Willie Wong Dec 23 '10 at 3:50
@J.M. The Lawson and Lim paper also refers to the following two which you may find interesting. (1) R.D. Nussbaum and J.E. Cohen, "The arithmetic-geometric mean and its generalizations for noncommuting linear operators", Ann. Scuola. Norm. Sup. Pisa Cl. Sci. 15 (1988) 239--308 (2) H. Maass, "Siegel's modular forms and Dirichlet series", Lecture Notes in Math. 216. The latter contains a few words on the geometry of the space of SPD matrices when endowed with the above Riemannian metric. – Willie Wong Dec 23 '10 at 20:30

You may want to take a look at exercise 6 on page 223 of Borwein and Borwein. (I am guessing anyone who asks about the AGM is familiar with this book.) There they discus a matrix version of the AGM and asks one to prove a connection between it and an elliptic integral.

They also give a paper by Stickel from 1985: "Fast compuation of Matrix Exponentials and Logarithms" as a reference. The journal name is Analysis.

-

I'm sorry that I'm putting this in an answer rather than a comment, but I'm too new to comment on your question. My thoughts for proving this process converges (question 1.) were along the lines of:

1. $A\,\&\,B$ both SPD $\implies AB$ and $A+B$ are also SPD
2. If $A$ and $B$ are symmetric $n\times n$ matrices, $A$'s eigenvalues $\{\lambda_n^A\} \subset [0, a]$, and $B$'s $\{\lambda_n^B\} \subset [0, b]$, then the eigenvalues $\{\lambda_N^{AB}\}$ of $AB$ are $\subset [0, ab]$. See, for instance, this discussion. Moreover, $\{\lambda_n^{A+B}\} \subset [0, a + b]$; see Terry Tao's notes on this.
3. I think convergence might follow by examining the AGM of the eigenvalues; that is, look at $AGM(\lambda_k^A, \lambda_k^B)$ for $k = 1,\ldots,n$. The difference between the arithmetic and geometric means gets arbitrarily small, and the eigenvalues are all positive; this iteration leads to a sequence of SPD matrices whose largest eigenvalue $\lambda_1 \to 0$, so the AGM in question does converge (perhaps).

Not sure, but thought I'd give my $\$0.02$. As @Slowsolver said, nice question! I also second the request of writing your algorithm in more traditional notation, since I am not as well-versed in Mathematica as others may be. Finally, I might add the tag "numerical-linear-algebra" to your question. EDIT: turns out I was wrong in part 1, as others have pointed out to me. My apologies! In light of their corrections, it may make sense to use the geometric mean given by SlowSolver. - If$A,B$are symmetric, that doesn't guarantee$AB$is symmetric. Consider$A = \mathop{diag}(1,2)$and$B = (2,1; 1,2)$.$BA = (2,2;1,4)$is not symmetric. The "right" definition of the product would be something more like$A^{1/2}BA^{1/2}\$, which while still non-commutative, at leasts maps SPD matrices to SPD matrices. – Willie Wong Dec 22 '10 at 18:41
thus, if A and B are spd, the product need not be spd; though it will have +ve eigenvalues. – user1709 Dec 22 '10 at 21:30
Thanks for pointing out my mistakes folks! I appreciate it. – user4689 Dec 23 '10 at 15:22