Arbitrary Precision, Nearly-Singular Matrices, and LinearSolve

Question

I have been trying to solve a nonlinear eigenvalue problem, $\mathbf{M}(\lambda) \mathbf{v} = 0$, in Mathematica using Newton's method. The core of the algorithm relies upon an inverse iteration, roughly solving $\mathbf{M}(\lambda_i) \mathbf{v}_{i+1} = \mathbf{M}'(\lambda_i)\mathbf{v}_i$ many times. Unfortunately, there are certain aspects of the problem (precise cancellations in the construction of $\mathbf{M}(\lambda_i)$) that force me to go beyond MachinePrecision. Since the matrices I'm working with become nearly singular as I approach the eigenvalue, I see dramatic performance hits in using LinearSolve to find $\mathbf{v}_{i+1}$. While I appreciate that inverting a nearly-singular matrix leads to large errors in the result (fortunately, inverse iteration is robust to this), I do not understand the performance impact or how to tell Mathematica to chill out about the errors. In fact, I'm having trouble getting any information about the inner workings of LinearSolve (other than it apparently uses arbitrary-precision BLAS routines).

A simple working example of this is:

prec = 60;
dim = 1000;
pertSize=1/10^9;

vec = RandomReal[{-1000, 1000}, dim, WorkingPrecision -> prec];
mat = Append[#[[;; -2]], #[[-2]] + RandomReal[{0, pertSize}, dim, WorkingPrecision -> prec]] &[RandomReal[{-1000, 1000}, {dim, dim}, WorkingPrecision -> prec]];
Block[{$MinPrecision = prec}, res = LinearSolve[mat, vec]; // Timing]

I construct a nearly-singular random matrix mat by replacing its last row by its second-to-last row and perturbing it by a small, random vector. The 1/10^9 thus controls how close the matrix is to being singular. On my laptop, the LinearSolve generally takes ~2-4 seconds with pertSize=1/10^9. However, with pertSize=1/10^10 the runtime jumps to ~120 seconds. The condition numbers for these matrices are not dramatically different, so I am having trouble understanding what is causing this ~60x slowdown.

I would love to understand what is causing this and how to force Mathematica not to worry too much about precisely finding the inverse.

My suspicion was that, because of the large number of "precision changing" operations inverting a nearly-singular matrix involves, the slowdown is mainly due to Mathematica having to repeatedly recast numbers back up to $MinPrecision. This explanation doesn't really make sense, though, since setting prec = MachinePrecision and playing around with pertSize, I see a similar performance impact once the matrix becomes singular.

An alternative is that LinearSolve switches to a different algorithm once the matrix becomes singular "enough." However, I do not know how to check this since I have not been able to find almost any information on LinearSolve, other than its anemic documentation. Also, it doesn't seem like fiddling with LinearSolve's Method helps.

I'd appreciate any insight or wisdom!

Answer 1

I have only a partial answer.

First, it's rare that a method will worry about finding the inverse. Something else is going on.

Some idea of what's happening is discussed here:

• Mark Sofroniou and Giulia Spaletta, Runtime Complexity Estimation for Accurate Solution of Linear System, Series in Applied Sciences. Volume 1, 2018 (also available here)

Your system leads to Method -> "IterativeRefinement" being selected. To quote Sofroniou and Spaletta: "The iterative refinement method becomes less efficient as the condition number of the system increases." I believe that's what's happening in your case: Convergence slows down as pertSize decreases. With random inputs, I observed that it sometimes happens with larger pertSize, presumably because the condition number happened to be larger than average.

Sofroniou and Spaletta suggests that the precision is managed internally by LinearSolve, so attempts to control precision with $MinPrecision may be futile. A message below suggests that some steps are done in machine precision. [Side note: Setting $MaxPrecision to the same number as $MinPrecision turns off precision tracking and speeds up arbitrary-precision arithmetic. Doing this in this case doesn't seem to matter much, since LinearSolve is the boss.]

The iterative refinement method has options. Unfortunately little is known about them:

Method -> {"IterativeRefinement",
  "ConditionNumber" -> Automatic, (* positive real; unknown effect *)
  "ConvergenceTest" -> True, (* True or False; unknown effect *)
  "MatrixDecomposition" -> Automatic, (* ? *)
  "MaxIterations" -> Automatic, (* positive integer, Infinity or Automatic *)
  "Norm" -> Infinity} (* "Norm" -> p: Norm[#, p], p = 1 or Infinity *)

The setting "ConditionNumber" -> 10.^16 yields a warning, and the call to LinearSolve returns unevaluated. Since the system with Automatic instead of 10.^16 is solved by LinearSolve, the warning seems a fatal error:

LinearAlgebra`LinearSolve`IterativeRefinement::lcondn: The matrix condition number 1.`*^16 is too large for starting precision MachinePrecision. The method may not converge.

Setting "MaxIterations" -> 10 has a predictable effect, but if the maximum is reached, you get this warning that also seems a fatal error:

LinearAlgebra`LinearSolve`IterativeRefinement::maxit: Maximum number of iterations 10 reached without convergence.

Setting ConvergenceTest has no effect on this or the "ConditionNumber" option.

I could not guess a valid "MatrixDecomposition" setting.

Finally, if you're still interested in this sort of problem and you have Premier Service, I'd suggest contacting WRI support and get their help. They can probably point out how to use these options. (They've done that for me in the past.)