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

Algorithm 1

ClearAll["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]