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