LSQR: Sparse Equations and Least Squares
- AUTHORS: Chris Paige, Michael Saunders.
- CONTRIBUTORS: James Howse, Michael Friedlander, John Tomlin, Miha Grcar, Jeffery Kline, Dominique Orban, Austin Benson, Victor Minden, Matthieu Gomez, Tim Holy.
- CONTENTS: Implementation of a conjugate-gradient type method
for solving sparse linear equations and
sparse least-squares problems:
\begin{align*}
\text{Solve } & Ax=b
\\ \text{or minimize } & \|Ax-b\|^2
\\ \text{or minimize } & \|Ax-b\|^2 + \lambda^2 \|x\|^2,
\end{align*}
where the matrix \(A\) may be square or rectangular
(over-determined or under-determined), and may have any rank.
It is represented by a routine for computing \(Av\) and \(A^T u\)
for given vectors \(v\) and \(u\).
The scalar \(\lambda\) is a damping parameter. If \(\lambda > 0\), the solution is "regularized" in the sense that a unique solution always exists, and \(\|x\|\) is bounded.
The method is based on the Golub-Kahan bidiagonalization process. It is algebraically equivalent to applying CG to the normal equation \( (A^T A + \lambda^2 I) x = A^T b, \) but has better numerical properties, especially if \(A\) is ill-conditioned.
NOTE: LSQR reduces \(\|r\|\) monotonically (where \(r = b - Ax\) if \(\lambda=0\)). This is desirable on compatible systems \(Ax=b\). On least-squares problems, if an approximate solution is acceptable (stopping tolerances quite large), LSMR may be a preferable solver because it reduces both \(\|r\|\) and \(\|A^T r\|\) monotonically and may be able to terminate significantly earlier.
If \(A\) is symmetric, use SYMMLQ, MINRES, or MINRES-QLP.
If \(A\) has low column rank and \(\lambda=0\), the solution is not unique. LSQR returns the solution of minimum length. Thus for under-determined systems, it solves the problem \(\min \|x\| \text{ subject to } Ax=b\). More generally, it solves the problem \(\min \|x\| \text{ subject to } A^T Ax=A^T b\), where \(A\) may have any shape or rank.
For \(\min \|x\| \text{ subject to } Ax=b\), consider using CRAIG.
- REFERENCES:
C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear equations and sparse least squares, TOMS 8(1), 43-71 (1982).
C. C. Paige and M. A. Saunders, Algorithm 583; LSQR: Sparse linear equations and least-squares problems, TOMS 8(2), 195-209 (1982). - RELEASE:
f77 and Matlab files are well tested.
C, C++ versions are Beta.
Windows DLL and .NET (C#) versions are Beta.
26 Oct 2012: f90 test program updated.
05 Aug 2013: Complex f90 version added.
30 Sep 2015: Link to Julia version added (Matthieu Gomez and Tim Holy).
05 Jul 2016: Link to second C++ implementation added (Tom Vercauteren).
DOWNLOADS
- Fortran 77 files
- Matlab files
24 Dec 2010: Function handles supported.
04 Sep 2011: Documentation slightly updated.
08 Mar 2019: colnorms.m included to give a diagonal preconditioner for LSQR. - C files 1:
Contributed Sep 1999 by James Howse (jhowse@lanl.gov), LANL.
An elaborate implementation with memory management. - C files 2:
Contributed Aug 2007 by
Michael Friedlander, UBC.
A more direct line-by-line translation from f77. - Fortran 90 files (for real \(A\)):
Second-generation f90 implementation.
Separate modules are used for LSQR, Check routines, and
example Test problems.
Contributed Sep 2007 by Michael Saunders (saunders@stanford.edu).
Revised 24 Oct 2007, 19 Dec 2008, 26 Oct 2012. - Fortran 90 files (for real or complex \(A\)): Complex implementation contributed by Austin Benson (arbenson@stanford.edu) and Victor Minden (victorminden@gmail.com), ICME, Stanford University. Includes real implementation.
- C++ files 1:
Contributed Oct 2007 by John Tomlin (tomlin@yahoo-inc.com),
Yahoo! Research.
C++ translation of C files in lsqr/c/clsqr1.zip. - C++ files 2: Contributed Jul 2016 by Tom Vercautereen, University College London http://cmictig.cs.ucl.ac.uk/.
- Windows DLL and .NET files:
Contributed Oct 2007 by Miha Grcar (miha.grcar@ijs.si),
Jozef Stefan Institute.
Windows DLL and .NET (C#) wrappers for LSQR based on the C++ version of LSQR. - Python files 1: Contributed Nov 2009 by Jeffery Kline.
- Python files 2: Contributed by Dominique Orban (dominique.orban@gerad.ca).
- Julia version: Contributed 2015 by Matthieu Gomez, Princeton University, and Tim Holy, Washington University in St Louis.
- GPU version included in MAGMA. See p2 of MAGMA handout.