Source-linked AI summary
LSMR: An iterative algorithm for sparse least-squares problems
David Fong, Michael Saunders
TL;DR
LSMR addresses iterative solution of sparse or operator-based linear systems and least-squares problems, including regularized variants. It combines Golub–Kahan bidiagonalization with MINRES applied to the normal equation and provides cheaply computable stopping quantities. Experiments report backward-error estimates close to optimal, often allowing earlier termination than LSQR, while reorthogonalization remains a scope boundary for future work.
Problem
LSMR addresses solving sparse or operator-based linear systems and least-squares problems, including regularized least squares.
Method
LSMR uses Golub–Kahan bidiagonalization and is analytically equivalent to MINRES applied to the normal equation.
Results
LSMR’s Stewart backward-error estimate is experimentally close to the optimal backward error at each iterate and often permits significantly earlier termination than LSQR.
Takeaways & Limitations
LSMR offers safer early termination because ∥A^T r_k∥ decreases monotonically, while ∥r_k∥ is also observed to decrease monotonically in practice.
Takeaways & Limitations
The paper leaves partial reorthogonalization techniques within LSMR for future research.
Abstract
from arXiv · showhide
An iterative method LSMR is presented for solving linear systems $Ax=b$ and least-squares problem $\min \norm{Ax-b}_2$, with $A$ being sparse or a fast linear operator. LSMR is based on the Golub-Kahan bidiagonalization process. It is analytically equivalent to the MINRES method applied to the normal equation $A\T Ax = A\T b$, so that the quantities $\norm{A\T r_k}$ are monotonically decreasing (where $r_k = b - Ax_k$ is the residual for the current iterate $x_k$). In practice we observe that $\norm{r_k}$ also decreases monotonically. Compared to LSQR, for which only $\norm{r_k}$ is monotonic, it is safer to terminate LSMR early. Improvements for the new iterative method in the presence of extra available memory are also explored.
1. Introduction.
LSMR is introduced for sparse or operator-based linear systems and least-squares problems. It combines Golub–Kahan bidiagonalization with MINRES-like behavior, supporting safer early termination and practical backward-error stopping.
- Problem setting: LSMR targets unsymmetric equations, linear least squares, and regularized least squares for matrices used through products Av and A^T u.The framework allows either m ≤ n or m ≥ n and includes λ ≥ 0 regularization.
- Method and comparison: LSMR is based on Golub–Kahan bidiagonalization and is analytically equivalent to MINRES applied to the normal equation.This differs from LSQR, which is equivalent to conjugate gradients on the normal equation.
- Method and comparison: The quantity ∥A^T r_k∥ decreases monotonically, while experiments observe that ∥r_k∥ also decreases monotonically and remains close to LSQR’s value.Because LSQR guarantees monotonicity only for ∥r_k∥, early termination is considered safer with LSMR.
- Stopping criteria: A cheaply computable LSMR backward error is close to the optimal backward error across many sparse least-squares test problems.The paper characterizes this closeness as an unexpected and desirable advantage for stopping decisions.
- Paper scope: The paper surveys the algorithm, norms and stopping criteria, singular systems, complexity, regularization, backward errors, and numerical experiments.It also states notation for norms, condition numbers, matrices, vectors, and scalars.
2. Derivation of LSMR.
The derivation uses Golub–Kahan bidiagonalization to express the normal-equation residual in a small least-squares subproblem. LSMR solves that subproblem through QR factorizations and forward substitutions.
- Golub–Kahan representation: Golub–Kahan generates expanding Krylov bases V_k and U_k, with V_k also linked to the symmetric Lanczos process for A^T A.The iterates are represented as x_k = V_k y_k.
- Reduced problem: LSMR chooses y_k to minimize ∥A^T r_k∥ rather than the residual norm minimized in LSQR’s subproblem.Using A^T r_k = A^T b − A^T A V_k y_k leads directly to the LSMR subproblem.
- Reduced problem: Efficient solution of this reduced least-squares subproblem is the heart of LSMR.The algorithm applies successive QR factorizations to obtain triangular systems.
- Iterative solution: Forward substitution produces auxiliary vectors that connect the reduced triangular solves to the iterate x_k.The construction is arranged so the resulting quantities can be updated iteratively.
2.4. Recurrence for xk.
The recurrence section turns the triangular-factorization representation into short updates for x_k and auxiliary vectors. Plane rotations expose the scalar relations needed to extend the recurrence one iteration at a time.
- Iterate recurrence: The iterate is represented through successive triangular solves, yielding x_k = W_k t_k = W̄_k z_k.This establishes x_0 = 0 and expresses each iterate through forward-substitution vectors.
- Iterate recurrence: The update x_k = x_{k−1} + ζ_k w̄_k adds one scaled auxiliary direction to the previous iterate.When k increases, existing quantities remain unchanged apart from one additional term.
- Rotations: The first QR factorization uses plane rotations on adjacent rows to transform the evolving bidiagonal structure.The rotations define the scalar quantities used in the recurrence.
- Rotations: A second sequence of rotations produces the barred triangular factors and defines the next auxiliary vectors through their last rows.These vectors support the continuing forward-substitution updates.
- Rotations: The two rotations are summarized by their effects on the computation, organizing the transformations that maintain the recurrence.
2.6. The two rotations.
The LSMR algorithm initializes the Golub–Kahan quantities and repeatedly applies two rotations before updating its recurrence vectors and iterate. Norms and condition estimates are maintained with constant work per iteration.
- Iteration: The scaled definitions h_k = ρ_k w_k and h̄_k = ρ_k ρ̄_k w̄_k make the forward substitutions more efficient.The resulting updates are used in the algorithm’s sixth pseudocode step.
- Iteration: Each iteration continues bidiagonalization, constructs and applies rotation P_k, constructs and applies rotation P̄_k, then updates h, h̄, and x.These are the core repeated stages of the LSMR pseudocode.
- Initialization: Initialization sets β_1u_1 = b, α_1v_1 = A^T u_1, ρ_0 = 1, and ρ̄_0 = 1.The starting state also uses x_0 = 0 and h_1 = v_1.
- Monitoring and stopping: LSMR derives ∥r_k∥, ∥A^T r_k∥, ∥x_k∥, and estimates of ∥A∥ and cond(A) for stopping rules using O(1) computation per iteration.
3. Norms and stopping rules.
LSMR derives cheaply updated norm estimates and stopping criteria from Golub-Kahan and QR-factorization quantities. Its key residual-gradient quantity decreases monotonically, while additional rules address consistency, conditioning, and singular termination.
- Norm estimates: The residual norm can be estimated from only the last two elements of ˜bk and the last element of ˜tk.The required vector ˜tk is computed by forward substitution, enabling a cheap update at each iteration.
- Norm estimates: ∥ATrk∥ = |¯ζk+1| is monotonically decreasing for LSMR.This follows from LSMR’s MINRES equivalence on the normal equations.
- Norm estimates: The estimate ∥xk∥ can be updated cheaply because only the last diagonal of eRk and the bottom 2 × 2 part of ˆRk change each iteration.Experiments indicate that ∥xk∥ > ∥xk−1∥ is true or nearly true at every iteration.
- Norm estimates: LSMR estimates ∥A∥ and cond(A) from the bidiagonal matrix Bk and updates these estimates cheaply.The estimates use singular-value interlacing and ratios of quantities derived from the factorization.
- Termination: When the Golub-Kahan process terminates because αk+1 = 0, a least-squares solution is obtained and, in the corresponding exact case, ∥rk∥ = 0.The latter condition implies Axk = b.
- Stopping rules: The stopping rules S1, S2, and S3 address consistent systems, inconsistent systems, and conditioning, respectively.S1 uses residual norm, S2 uses ∥ATrk∥, and S3 stops when cond(A) reaches CONLIM.
4. Characteristics of the solution on singular systems.
For rank-deficient A, the normal equations are singular but consistent, and LSMR returns the minimum-norm least-squares solution. Its storage and work requirements are compared with LSQR and MINRES for sparse least-squares problems.
- Singular systems: The minimum-norm result follows because differences between solutions lie in N(A) and are orthogonal to the Golub-Kahan basis vectors.The proof uses ATAp = 0 to establish Ap = 0.
- Singular systems: At convergence, LSMR reduces its projected problem to the same least-squares subproblem solved by LSQR.Because Bk has full rank, the converged yk equals the LSQR solution for the projected problem.
- Complexity: The comparison with LSQR and MINRES considers vector storage and per-iteration floating-point multiplications, excluding storage for A and b.The requirements are summarized for A ∈ R^{m×n}, where m may be considerably larger than n.
5. Regularized least squares.
LSMR is extended to regularized least-squares problems with λ ≥ 0 while retaining the main unregularized algorithmic structure. Regularization changes the residual norm and requires modified rotation updates and stopping procedures.
- Regularized formulation: The regularized LSMR algorithm solves the regularized least-squares problem for a given λ.Its main algorithm follows the unregularized case after the regularization-specific transformations.
- Regularized formulation: Regularization changes the residual norm and the rotations used during each iteration.The effects of regularization on the rotations are summarized separately before the algorithmic steps.
- Algorithm: The rotation updates transform αk, βk+1, and related quantities needed for the regularized recurrence.The displayed recurrences define the updated scalars used by subsequent steps.
- Termination: The regularized algorithm terminates when one of the stopping criteria from section 3.6 is satisfied.The stopping criteria are applied after updating the iteration quantities.
- Backward-error connection: For inconsistent problems with uncertainty in A, normwise backward error measures the perturbation needed to make an approximate solution exact.The perturbation E is defined through the perturbed normal equations.
6. Backward errors.
The paper studies backward errors for approximate least-squares solutions and develops a cheaply computable estimate used by LSMR. Experiments show that this estimate is close to optimal and often enables earlier termination than LSQR.
- Backward-error definitions: The normwise backward error µ(x) measures the smallest perturbation to A that makes x solve the perturbed least-squares problem exactly.It is characterized through a perturbation satisfying the perturbed normal equations.
- Backward-error definitions: Evaluating µ(x) directly is generally too expensive, so the method uses approximations.Several approximate backward-error expressions are discussed, including Stewart’s E1 and E2 estimates.
- E2 estimate: LSQR and LSMR compute ∥E2∥ for each iterate from accurately estimated ∥rk∥ and ∥ATrk∥ at almost no cost.Because E2xk = 0 under orthogonality of Vk, xk and rk are exact for the perturbed problem.
- Stopping rule: Stopping rule S2 requires ∥E2∥ ≤ ATOL∥A∥, giving LSMR an advantage for early termination.The rule is tied to uncertainty in A for inconsistent problems.
- Optimal backward error: The computed eµ(xk) is experimentally very close to the true optimal backward error µ(xk) for each LSQR and LSMR iterate.The estimate is obtained using an economy-size sparse QR factorization.
- Empirical comparison: LSMR’s rule S2 terminates significantly sooner than LSQR on most inconsistent test cases.The comparison is reported for the experiments referenced in Theorem 6.1 and Figures 7.2 and 7.3.
7. Numerical results.
Numerical experiments on sparse least-squares and square systems find that LSMR generally provides monotonic residual-related measures and near-optimal backward-error estimates, while LSQR can have smaller solution errors. Reorthogonalization and partial reorthogonalization can improve behavior when memory is available.
- Least-squares problems: LSMR residual norms appear monotonic and nearly as small as LSQR residual norms on almost all tested problems.A rare case shows substantially larger LSMR residuals during early iterations.
- Backward error: LSMR’s computable backward-error estimate eµ(xk) is usually monotonic, no larger than LSQR’s, and nearly equal to its optimal backward error.The experiments also report ∥ELSQR 1 ∥ and ∥ELSMR 2 ∥ as often close to the optimal backward error.
- Solution accuracy: LSQR typically has smaller solution errors than LSMR, although both methods converge to the minimum-norm least-squares solution for singular systems.For selected square systems, LSMR residual norms can exhibit prolonged plateaus even while ∥A^T r_k∥ decreases appropriately.
- Reorthogonalization: Reorthogonalization can make a large difference in difficult cases, while partial reorthogonalization with l = 10 or 50 provides partial speedups and exposes a storage-time tradeoff.The benefit depends strongly on the computational cost of Av and A^T u; related techniques for both vector sets remain future work.
8. Summary.
The summary presents LSMR as offering practical stopping advantages over LSQR while supporting accurate Golub–Kahan computations with selective reorthogonalization.
- The summary positions LSMR as having advantages over the widely adopted LSQR algorithm.
- LSMR often terminates significantly sooner than LSQR because its Stewart backward-error estimate is experimentally close to the optimal backward error.The stopping framework estimates residuals, normal-equation residuals, solution norms, matrix norms, and condition numbers cheaply.
- The method provides theoretical and practical stopping criteria for linear systems and least-squares problems, including optional Tikhonov regularization.
- Full reorthogonalization experiments show that reorthogonalizing either V_k or U_k preserves high accuracy, so reorthogonalizing both is unnecessary.The paper identifies this as potentially useful beyond LSMR for other applications of the Golub–Kahan process.
- Matlab, Python, and Fortran 90 implementations of LSMR are available, and each permits local reorthogonalization of V_k.
Appendix A. Proof of Lemma 3.1.
The appendix proves Lemma 3.1 by induction, establishing equality between transformed scalar quantities and agreement among leading components of transformed vectors.
- The appendix begins by describing the effects of the rotations P_k and eP_{k−1} on the transformed quantities.
- The appendix defines cumulative products s^(k) and ¯s^(k), deriving their relation from equations (3.3) and (2.4).
- An induction argument is used to show that ˜τ_i = ˜β_i for every index i.The base case i = 1 is introduced before subsequent equalities are justified from earlier relations and the induction hypothesis.
- By induction, ˜τ_i = ˜β_i holds for i = 1, 2, …, and the first k−1 elements of ˜b_k and ˜t_k are equal at iteration k.