Function Repository Resource:

MatrixGeometricMean

Source Notebook

Compute the geometric mean of two matrices

Contributed by: Jan Mangaldan

ResourceFunction["MatrixGeometricMean"][a,b]

gives the geometric mean of the matrices a and b.

Details

The geometric mean of two matrices a and b is the unique positive definite matrix solution x of the equation Image, and is the generalization of the usual geometric mean for scalars.
The matrices a and b can be numerical or symbolic, but must be Hermitian and positive definite.
The geometric mean of two matrices can have complex-valued elements.

Examples

Basic Examples (2) 

The geometric mean of two exact 2×2 symmetric positive definite matrices:

In[1]:=
ResourceFunction["MatrixGeometricMean"][( {
    {2, 1},
    {1, 2}
   } ), ( {
    {3, 2},
    {2, 3}
   } )] // Simplify
Out[1]=
Image

The geometric mean is also symmetric and positive definite:

In[2]:=
SymmetricMatrixQ[%] && PositiveDefiniteMatrixQ[%]
Out[2]=
Image

Scope (2) 

Two symmetric positive definite matrices:

In[3]:=
ma = HilbertMatrix[4];
mb = Array[Min, {4, 4}];

Compute the geometric mean with machine arithmetic:

In[4]:=
ResourceFunction["MatrixGeometricMean"][N[ma], N[mb]]
Out[4]=
Image

Compute the geometric mean with 24-digit precision arithmetic:

In[5]:=
ResourceFunction["MatrixGeometricMean"][N[ma, 24], N[mb, 24]]
Out[5]=
Image

Compute the geometric mean of two random Hermitian positive definite matrices:

In[6]:=
ma = ConjugateTranspose[#] . # &[RandomComplex[1 + I, {4, 4}]];
mb = ConjugateTranspose[#] . # &[RandomComplex[1 + I, {4, 4}]];
ResourceFunction["MatrixGeometricMean"][ma, mb]
Out[7]=
Image

Properties and Relations (6) 

MatrixGeometricMean of two 1×1 matrices is equivalent to GeometricMean:

In[8]:=
a = RandomReal[];
b = RandomReal[];
{ResourceFunction["MatrixGeometricMean"][{{a}}, {{b}}], GeometricMean[{a, b}]}
Out[9]=
Image

The geometric mean is symmetric in its arguments:

In[10]:=
n = 3;
a = Transpose[#] . # &[RandomReal[1, {n, n}]];
b = Transpose[#] . # &[RandomReal[1, {n, n}]];
ResourceFunction["MatrixGeometricMean"][a, b] == ResourceFunction["MatrixGeometricMean"][b, a]
Out[11]=
Image

The geometric mean of a matrix and the identity is equivalent to the square root of the matrix:

In[12]:=
n = 3;
a = Transpose[#] . # &[RandomReal[1, {n, n}]];
ResourceFunction["MatrixGeometricMean"][a, IdentityMatrix[n]] == MatrixPower[a, 1/2]
Out[13]=
Image

If two matrices commute, their geometric mean is equivalent to the square root of their product:

In[14]:=
a = {{34, 12}, {12, 41}};
b = {{2100, 300}, {300, 2275}};
a . b == b . a
Out[15]=
Image
In[16]:=
ResourceFunction["MatrixGeometricMean"][a, b] == MatrixPower[a . b, 1/2] // Simplify
Out[16]=
Image

The geometric mean can be expressed in terms of MatrixPower:

In[17]:=
n = 3;
a = Transpose[#] . # &[RandomReal[1, {n, n}]];
b = Transpose[#] . # &[RandomReal[1, {n, n}]];
eps = 10^(-8);
mgm = ResourceFunction["MatrixGeometricMean"][a, b];
diff1 = Norm[mgm - MatrixPower[a . Inverse[b], 1/2] . b, Infinity];
diff2 = Norm[mgm - b . MatrixPower[Inverse[b] . a, 1/2], Infinity];
{diff1 < eps, diff2 < eps}
Out[18]=
Image

The geometric mean can be expressed as an integral:

In[19]:=
n = 3;
a = Transpose[#] . # &[RandomReal[1, {n, n}]];
b = Transpose[#] . # &[RandomReal[1, {n, n}]];
Norm[ResourceFunction["MatrixGeometricMean"][a, b] - 2/\[Pi] NIntegrate[
     Inverse[(1 + t) Inverse[a] + (1 - t) Inverse[b]]/Sqrt[
     1 - t^2], {t, -1, 1}], \[Infinity]] // Chop
Out[20]=
Image

Version History

  • 1.0.0 – 06 May 2022

Source Metadata

Related Resources

Author Notes

The matrix geometric mean for a pair of matrices is uniquely defined. The matrix geometric mean of three or more HPD matrices can be defined in many inequivalent ways, and all of them are nontrivial to implement robustly. Thus, that generalization is left for a future version.

License Information