lineax icon indicating copy to clipboard operation
lineax copied to clipboard

Implementation of LSMR for iterative least squares.

Open f0uriest opened this issue 11 months ago • 1 comments

Tests seem to be passing locally. I added LSMR in the places it seemed to make sense (mainly places where I saw SVD which solves similar problems. Let me know if there are other tests that should be added.

A few things that could use some input:

  • There are a large number of state variables to keep track of. Right now I'm just using nested tuples but there's probably a more elegant way.

Resolves #85

f0uriest avatar Mar 16 '24 05:03 f0uriest

The original version has 7 possible exit conditions:

  1. exceeds maxiter
  2. x solve Ax=b to user specified tolerance (ie, no least squares needed, problem appears well posed)
  3. x is the least squares solution to user specified tolerance
  4. cond(A) larger than user specified conlim
  5. x solve Ax=b to tolerance=eps (same as 2, but guards against users setting tolerances to zero?)
  6. x is the least squares solution to tolerance=eps (same as 3?)
  7. cond(A) > 1/eps (same as 4?)

I think we could get rid of 5-7 provided that maxiter is always finite. 2 and 3 could probably be combined into a standard "successful" exit assuming the user doesn't really care whether its the "true" solution or least squares. The main new one would be the condition number one. We could maybe shoehorn that into lx.RESULTS.singular or lx.RESULTS.breakdown? or could add a new one just for that.

f0uriest avatar Apr 27 '24 19:04 f0uriest

I know it's bad form to comment just when you want to 👍🏻 , but having LSMR in lineax basically made my problem go from "untractable and 10x slower with jax than scipy" to tractable.

The problem here is one where I have a large (complex) linear system and I never want to construct the design matrix, so iterative methods are all I can really use. QR or SVD had huge memory overheads and were very slow. The Jax implementation for GMRES (in lineax or jax.scipy.sparse.linalg.gmres) doesn't seem to accept non-square matrices, although the scipy.sparse.linalg.gmres implementation does.

So this is a big 👍🏻 to this feature.

andycasey avatar Oct 15 '24 17:10 andycasey