A Higher Order Iterative Method for Computing the Drazin Inverse

Algorithm 3

PM[X_]  :=

With[{Id = SparseArray[{{i_, i_}  -> 1.}, {n, n}]],

         X1 = A.X; X2 = −7 Id + X1.(9 Id + X1.(−5 Id + X1)); X3 = X1.X2;

         (−1/8) X.X2.(12 Id + X3.(6 Id + X3))

          ];

InitialMatrix[A_]  := 1/SingularValueList[A, 1] [[1]]^∧2  ConjugateTranspose[A];

InitialDrazin[A_]  := 2/Tr[MatrixPower[A, k + 1]]  MatrixPower[A, k];

DrazinInverse[A_, tolerance_]  := If[k == 0,

      Module[{X0 = InitialMatrix[A]}, FixedPoint[(PM[#]  &), X0,

        SameTest -> (Norm[#1 - #2, Infinity]  <= tolerance  &)]],

      Module[{X0 = InitialDrazin[A], FixedPoint[(PM[#]  &), X0,

        SameTest -> (Norm[#1 −  #2, Infinity]  <= tolerance  &)]]

      ];