# 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 '13 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 '13 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 '13 at 9:31
Can you write it down, please? Does it reduce the whole computational time? –  Fazlollah Soleymani Aug 2 '13 at 10:08

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.

Max@Abs[A.final - IdentityMatrix@Length[A]]
(* -> 049.48507777383387 *)
(* -> 048.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.