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;

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!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.