# Systems Optimization Laboratory

## LSMR: Sparse Equations and Least Squares

• AUTHORS: David Fong, Michael Saunders.
• CONTRIBUTORS: Dominique Orban, Austin Benson, Victor Minden.
• 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 MINRES 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.

Special feature: Both $$\|r\|$$ and $$\|A^T r\|$$ decrease monotonically, where $$r = b - Ax$$ is the current residual. For LSQR, only $$\|r\|$$ is monotonic. With loose stopping tolerances, LSMR may be able to terminate significantly sooner than LSQR.

Special application: To find a null vector of a singular (square or rectangular) matrix $$A$$, apply LSMR to the system $$\min \|A^T x - b\|$$ with any nonzero vector $$b$$ (e.g. a random $$b$$). At a minimizer, the residual vector $$r = b - A^T x$$ will satisfy $$Ar=0$$. See [1] for examples.

• REFERENCES:
[1] S.-C. T. Choi (2006). Iterative Methods for Singular Linear Equations and Least-Squares Problems, PhD thesis, ICME, Stanford University.

[2] D. C.-L. Fong and M. A. Saunders, LSMR: An iterative algorithm for sparse least-squares problems, SIAM J. Sci. Comput. 33:5, 2950-2971, published electronically Oct 27, 2011.

• RELEASE:
14 Apr 2010: Matlab implementation.
03 Jun 2010: Python implementation.
17 Jul 2010: f90 implementation (Beta).
07 Sep 2010: f90 local reorthogonalization fixed.
26 Oct 2012: f90 test program updated.
05 Aug 2013: Complex f90 version added.