# How to use adaptive precision in matrix computations?

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.

-
This is a nice idea. However, I do think it's more of an idea than a concrete question. What is called for here is a detailed analysis of how precision is carried through every operation used in the calculation, which will be specific to each algorithm for which one might wish to implement this. I'm sure that this specific problem has been well studied, so the existing literature may offer a good starting point rather than having to conduct this analysis from scratch. –  Oleksandr R. Aug 2 at 9:16
Thanks for your comment. I do agree. However, I think some built-in functions of MMA already do such an action. An example is the function FindRoot[], which applies the precision of the input data (function and the initial guess) and then improve it per cycle. Also, applying the adaptive precision in matrix calculations (such as the above concrete question/idea) is very good in high precision computing environment. It should reduce the computational time dramatically for large scale problems. –  Fazlollah Soleymani Aug 2 at 9:21
Mathematica's precision tracking is either done on a case-by-case basis for e.g. special functions and linear algebra, or in an approximate and heuristic way for e.g. FindRoot. These heuristics are not really safe to use when you need accuracy up to 0.5-1 ULPs as in most numerical linear algebra contexts; a careful analysis is called for instead. Anyway, a possible non-adaptive approach would be to get your initial result in machine precision, then SetPrecision` the machine-precision estimate to the final precision (plus a wide margin of safety) and continue from there. –  Oleksandr R. Aug 2 at 9:31
Can you write it down, please? Does it reduce the whole computational time? –  Fazlollah Soleymani Aug 2 at 10:08