Systems Optimization Laboratory

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).

• 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.