How to use adaptive precision in matrix computations?

Question

I wish to compute the pseudo inverse of rectangular (or square) matrices by the cubically method of Chebyshev given by $X_{k+1}=X_k(3I-AX_k(3I-AX_k))$ where $X_0=\frac{1}{\|A\|_F^2}A^*$. The procedure is simple as I wrote in the following sample example:

Clear["Global`*"]
SeedRandom[1234];
m = 400; n = 500;
Id = SparseArray[{{i_, i_} -> 1.}, {m, m}, 0.];
A = RandomReal[{-100, 100}, {m, n}];
X[0] = 1/Norm[A, "Frobenius"]^2 ConjugateTranspose[A];
max = 50; k = 1; stop = 1;
While[k <= max && stop > 10^-2, {
    AX = A.X[k - 1];
    X[k] = X[k - 1].(3 Id - AX.(3 Id - AX));
    stop = Norm[X[k] - X[k - 1], "Frobenius"]/Norm[X[k], "Frobenius"];
    }; k++]; // AbsoluteTiming

{(k - 1), Norm[Chop[X[k - 1] - PseudoInverse[A]], 1]}

Running the above piece of code, will produce the pseudo inverse in 12 number of iterationsin 1.8 seconds (in my computer). The only problem of the Chebyshev method is that it is slow at the beginning of the process. I mean the third order of convergence can only be seen in the last four iterations. On the other hand, my aim and question is to reduce the computational time/effort of implementing this method in MMA.

Since the method is asymptotically stable, I think one way for reducing the computational time is to use adaptive precision. So, is there any way to apply adaptive precision in the implementation of the above piece of code? I mean, it would be too nice if we could run the method in Single Precision at the beginning of the process and then when the stopping termination ("stop") reaches a pre-tolerance, then we back to the machine precision (double precision)?

Note that I used $SetAccuracy[exp,digits]$, but I think it is so costly for matrix-matrix multiplications and is not useful. Any idea that could be helpful is fully appreciated.

Answer 1

Yes, it is possible to do this with adaptive precision. Unfortunately, Mathematica doesn't support anything less than double precision, so you can't make many gains at the beginning of the process unless you want to defer to a LibraryLink program and do it all in C instead. But, if you're interested in a high-precision result, the approach can definitely be useful. Here I'll show the result for a 50-digit calculation.

First of all:

SeedRandom[1234];
With[{m = 200, n = 300, prec = 50},
  A = RandomReal[{-100, 100}, {m, n}, WorkingPrecision -> prec];
  Id = SetPrecision[SparseArray[{{i_, i_} -> 1}, {m, m}, 0], prec];
 ];

Now, we get a machine-precision estimate:

estimate = With[{
     A = SetPrecision[A, MachinePrecision],
     Id = SetPrecision[Id, MachinePrecision]
    },
    Module[{
      X, k, tol = 1, maxit = 50,
      pinv = PseudoInverse[A]
     },
     X[0] = 1/Norm[A, "Frobenius"]^2 ConjugateTranspose[A];
     k = 0;
     While[k <= maxit && tol > Sqrt[10^-Precision[A]], k++;
      With[{AX = A.X[k - 1], threeId = 3 Id},
       X[k] = X[k - 1].(threeId - AX.(threeId - AX))
      ];
      tol = Norm[X[k] - X[k - 1], "Frobenius"]/Norm[X[k], "Frobenius"];
      Print@Row[{"k = ", k, ", tol = ", tol, ", norm = ", Norm[X[k] - pinv, 1]}];
     ];
     X[k]
    ]
   ];

k = 1, tol = 0.663897, norm = 0.0294513
k = 2, tol = 0.658436, norm = 0.0292652
k = 3, tol = 0.642703, norm = 0.0287461
k = 4, tol = 0.601449, norm = 0.0274857
k = 5, tol = 0.522307, norm = 0.0248225
k = 6, tol = 0.440229, norm = 0.0201895
k = 7, tol = 0.365339, norm = 0.0135799
k = 8, tol = 0.258934, norm = 0.00563918
k = 9, tol = 0.109346, norm = 0.000564096
k = 10, tol = 0.00969869, norm = 7.2773*10^-7
k = 11, tol = 0.0000108769, norm = 1.67539*10^-15
k = 12, tol = 2.42866*10^-14, norm = 1.36833*10^-16

Not counting the pseudoinverse operation using PseudoInverse (which we of course just use to check that our result is converging to the right value), it took 0.125 seconds, which is almost negligible.

Now we would like to "polish" the result up to 50 places of accuracy. First of all though, how long does the 50-digit PseudoInverse take?

pinv = PseudoInverse[A]; // AbsoluteTiming
(* -> 66.859375 seconds *)

Now ours:

final = Block[{
    $MinPrecision = Precision[A],
        $MaxPrecision = Precision[A]
   },
   Module[{
     X, k, tol = 1, maxit = 50
    },
    X[0] = SetPrecision[estimate, Precision[A]];
    k = 0;
    While[k <= maxit && tol > Sqrt[10^-Precision[A]], k++;
     With[{AX = A.X[k - 1], threeId = 3 Id},
      X[k] = X[k - 1].(threeId - AX.(threeId - AX))
     ];
     tol = Norm[X[k] - X[k - 1], "Frobenius"]/Norm[X[k], "Frobenius"];
     Print@Row[{"k = ", k, ", tol = ", tol, ", norm = ", Norm[X[k] - pinv, 1]}];
    ];
    X[k]
   ]
  ];

k = 1, tol = 6.560807*10^-16, norm = 7.7(...)23527321325725384674*10^-17
k = 2, tol = 4.554664*10^-46, norm = 7.7(...)24243494475169487393*10^-17

It took 8.891 seconds. So, we have done the same as PseudoInverse, but in 8.891 + 0.125 = 9.016 seconds, rather than 66.859.

How about accuracy?

Max@Abs[A.final - IdentityMatrix@Length[A]]
(* -> 0``49.48507777383387 *)

Max@Abs[A.pinv - IdentityMatrix@Length[A]]
(* -> 0``48.37090027415175 *)

Well, we got the right answer by both methods, so it's something of a moot point. But, with the iterative method, we managed to lose only half a digit of precision, whereas using PseudoInverse directly lost more than 1.5 digits.