PaperlayerView exact revision on arXiv

Source-linked AI summary

MINRES-QLP: a Krylov subspace method for indefinite or singular symmetric systems

Sou-Cheng T. Choi, Christopher C. Paige, Michael A. Saunders

arXiv:1003.4042v3math.NAcs.CEcs.NAstat.CO

TL;DR

For singular incompatible symmetric systems, MINRES may find a least-squares solution that is not the minimum-length solution. The paper introduces MINRES-QLP, using a QLP decomposition within the Lanczos framework, and reports improved accuracy over MINRES on ill-conditioned systems. It also develops norm estimates and stopping rules, while identifying solution-norm bounds and practical thresholds as future work.

  • Problem

    For incompatible singular symmetric systems, existing methods may fail or return a least-squares solution that is not the minimum-length pseudoinverse solution.

    Source: §1 · p. 2

  • Method

    MINRES-QLP applies a QLP decomposition to the reduced Lanczos subproblem and derives norm estimates and stopping conditions for singular problems.

    Source: §5 · p. 7, §9 · p. 22

  • Results

    MINRES-QLP can improve accuracy over MINRES on ill-conditioned systems, and the numerical examples show stronger regularization as it monitors the rank of T_k.

    Source: §1 · p. 2, §9 · p. 22

  • Takeaways & Limitations

    MINRES-QLP returns the minimum-length solution for compatible or incompatible systems, and transfers from MINRES only when the estimated condition number exceeds a specified value.

    Source: §1 · p. 2, §9 · p. 21

    Estimating an a priori solution-norm bound and choosing a practical threshold for the stopping condition remain future research topics.

    Source: §9 · p. 22

Abstract

from arXiv · show

CG, SYMMLQ, and MINRES are Krylov subspace methods for solving symmetric systems of linear equations. When these methods are applied to an incompatible system (that is, a singular symmetric least-squares problem), CG could break down and SYMMLQ's solution could explode, while MINRES would give a least-squares solution but not necessarily the minimum-length (pseudoinverse) solution. This understanding motivates us to design a MINRES-like algorithm to compute minimum-length solutions to singular symmetric systems. MINRES uses QR factors of the tridiagonal matrix from the Lanczos process (where R is upper-tridiagonal). MINRES-QLP uses a QLP decomposition (where rotations on the right reduce R to lower-tridiagonal form). On ill-conditioned systems (singular or not), MINRES-QLP can give more accurate solutions than MINRES. We derive preconditioned MINRES-QLP, new stopping rules, and better estimates of the solution and residual norms, the matrix norm, and the condition number.

1. Introduction.

The paper develops MINRES-QLP to obtain minimum-length solutions for compatible and incompatible symmetric systems, addressing cases where MINRES may return a least-squares solution that is not the pseudoinverse solution. It also sets out the paper’s notation and roadmap.

1. Introduction.

  • The paper uses x† for the minimum-length or pseudoinverse solution and defines ill-conditioning through the ratio of the largest to smallest nonzero eigenvalue magnitudes.A singular matrix can be either well-conditioned or ill-conditioned; the pseudoinverse solution is continuous under perturbations that preserve rank.Source: §1 · p. 1
  • For singular incompatible systems, SYMMLQ iterates can diverge toward a nullvector, while MINRES can return a least-squares solution that is not the minimum-length solution.MINRES residual norms decrease monotonically, and an appropriate stopping rule involving ||Ar_k|| can make it reliable on these systems.Source: §1 · p. 2
  • MINRES-QLP is designed to solve compatible or incompatible symmetric systems reliably and return the unique minimum-length solution, with theoretical and numerical support for improved accuracy over MINRES on ill-conditioned systems.The paper also notes applications of incompatible systems, including discretized semidefinite Neumann problems and singular systems with measurement errors in b.Source: §1 · p. 2

1.2. Overview.

  • The paper reviews Lanczos, MINRES, and QLP decomposition before introducing MINRES-QLP, deriving norm estimates and preconditioning, and presenting numerical experiments.For incompatible systems, Ax ≈ b denotes minimizing ||Ax − b||.Source: §1.2 · p. 2, §1.1 · p. 2

2. The Lanczos process.

The Lanczos process builds an orthonormal Krylov basis and a tridiagonal matrix that reduce the symmetric-system problem to a smaller subproblem. The rank and singularity properties of these tridiagonal matrices govern when the reduced problems are solvable.

2. The Lanczos process.

  • Lanczos generates vectors v_k and tridiagonal matrices T_k from A and b, with ||v_k|| = 1 and, in exact arithmetic, AV_ℓ = V_ℓT_ℓ when the process terminates.In practice, termination at β_{ℓ+1} = 0 is unlikely without reorthogonalization, though the Lanczos relation holds to machine precision.Source: §2 · p. 2

2.1. Properties of the Lanczos process.

  • The Krylov subspace K_k(A,b) is the span of b, Ab, through A^{k−1}b, and equals the span of V_k.A shift from A to A − σI shifts T_k by the same amount while leaving V_k unchanged.Source: §2.1 · p. 3
  • For k < ℓ, T_k has full column rank; if an intermediate T_k is singular, it has exactly one zero eigenvalue and adjacent T matrices are nonsingular.T_ℓ is nonsingular exactly when b ∈ range(A).Source: §2.1 · p. 3

3. MINRES.

MINRES uses the Lanczos reduction and an expanding QR factorization to minimize residuals over Krylov subspaces. Its behavior differs between compatible and incompatible systems, and residual-based stopping and finite-precision effects require care.

3. MINRES.

  • At iteration k, MINRES minimizes ||Ax_k − b|| over x_k ∈ K_k(A,b), using the Lanczos relation to express the residual through the reduced tridiagonal problem.An expanding QR factorization produces upper-tridiagonal R_k, and the iteration updates x_k from the reduced solve.Source: §3 · p. 3
  • On compatible systems, MINRES returns the minimum-length solution with zero residual; on incompatible systems, its final least-squares solution need not be minimum-length.MINRES-QLP instead returns the minimum-length solution in the incompatible case.Source: §3 · p. 3

3.1. Further details of MINRES.

  • For k < ℓ, a singular T_k can prevent a strict residual decrease for that step, but this cannot occur on consecutive iterations.The residual norm is nonincreasing and tends to zero when Ax = b is compatible.Source: §3.1 · p. 4

3.2. Compatible systems.

  • When b ∈ range(A), the final MINRES iterate solves Ax = b and is the minimum-length solution, even if A is singular.In this case T_ℓ is nonsingular.Source: §3.2 · p. 4

3.3. Incompatible systems.

  • For incompatible systems, MINRES can attain the optimal residual without producing the minimum-length solution, as illustrated by A = diag(1, 1, 0) and b = e.The example’s minimum-length solution is [1 1 0]^T, whereas MINRES returns e; both have residual e_3.Source: §3.3 · p. 4

3.4. Norm estimates in MINRES.

  • Because incompatible systems never have zero residual, the paper proposes stopping based on ||Ar_k|| and gives recurrence estimates for ||Ar_k|| and ||Ax_k||.The norms ||r_k|| and ||Ax_k|| are monotonic, while ||Ar_k|| typically is not.Source: §3.4 · p. 5

3.5. Effects of rounding errors in MINRES.

  • Ill-conditioned R_k can make MINRES vulnerable to cancellation, and reported rounding-error examples show its recurred residual norm may stagnate or increase.MINRES-QLP’s key feature is that one lower-triangular solve suffices without storing all of V_k.Source: §3.5 · p. 6

4. Orthogonal decompositions for singular matrices.

The section describes Stewart’s pivoted QLP decomposition and the unpivoted QLP variant. Its lower-triangular factor provides singular-value estimates that can be more accurate than those from the upper-triangular factor.

4. Orthogonal decompositions for singular matrices.

  • Stewart’s pivoted QLP decomposition applies consecutive QR factorizations with column interchanges, first to A and then to R^T.The permutations are selected to maximize the next diagonal element at each stage.Source: §4 · p. 6
  • The diagonal L-values can estimate singular values more accurately than the R-values, particularly for the largest and smallest singular values.Source: §4 · p. 6
  • Without either permutation, the decomposition is called QLP; the second permutation in the pivoted version orders the L-values decreasingly, which is not always necessary.Source: §4 · p. 6

5. MINRES-QLP.

MINRES-QLP applies a QLP decomposition to the Lanczos subproblem to handle ill-conditioned or singular symmetric systems. At termination it returns the minimum-length least-squares solution, and its L_k diagonals can estimate the extreme nonzero singular values of A.

5. MINRES-QLP.

  • MINRES-QLP uses a QLP decomposition of each Lanczos subproblem T_k, applying an LQ factorization after MINRES’s QR factorization.The QLP handles rank-deficient subproblems to obtain the unique minimum-length solution.Source: §5 · p. 7

5.1. The QLP factorization of Tk.

  • The QLP factorization uses unpermuted factors for inexpensive updates, with rank-revealing properties tending to be retained.Left and right reflectors are interleaved; at each iteration, the right reflectors follow the left reflector so future updates affect only the last 2 × 2 principal submatrix.Source: §5.1 · p. 7, Figure 5.1 · p. 8

5.2. The diagonals of Lk.

  • The diagonals of L_k approximate A’s extreme nonzero singular values better than those of R_k, and their ratio estimates κ(A).For matrix ID 1177, the comparison is illustrated at k = 19.Source: §5.2 · pp. 7–8, Figure 5.2 · p. 8

5.3. Solving the subproblem.

  • When L_k is rank deficient at termination, MINRES-QLP solves the last two equations of the preceding subproblem to construct its minimum-length solution.For full rank, it instead solves the last three equations of L_k u_k = t_k.Source: §5.3 · p. 9

5.4. Termination.

  • At termination, MINRES-QLP returns x†, the minimum-length least-squares solution for both compatible and incompatible systems.For incompatible systems, its construction selects the minimum-length solution among the least-squares solutions in the Krylov subspace.Source: §5.4 · pp. 9–10

5.5. Transfer from MINRES to MINRES-QLP.

  • MINRES can transfer to MINRES-QLP when the estimated condition number of T_k exceeds trancond, balancing the methods’ behavior with MINRES-QLP’s added work and storage.A trancond above 1/ε keeps MINRES throughout, while trancond = 1 uses MINRES-QLP from the start.Source: §5.5 · p. 10

5.6. Comparison of Lanczos-based solvers.

  • The solver comparison covers subproblem definitions, bases, solution estimates, flop counts, and memory requirements.The methods all require storage for v_k, v_k+1, x_k, and A v_k, with additional work vectors varying by method.Source: §5.6 · p. 10

6. Stopping conditions and norm estimates.

This section derives recurrence-based residual and norm estimates for stopping tests, along with estimates of matrix, condition, solution, and projection norms.

6. Stopping conditions and norm estimates.

  • MINRES-QLP uses norm estimates and input parameters to formulate three groups of stopping conditions, including normwise relative backward error tests.The stopping-condition table identifies tol, maxit, maxcond, and maxxnorm as input parameters; its norms and κ(A) are estimated by MINRES-QLP.Source: §6 · p. 10, Table 6.1 · p. 12

6.1. Residual and residual norm.

  • The residual recurrence gives φ_k = ∥r_k∥, with zero residual at termination for full rank and a nonzero residual when the terminal rank is deficient.The section also derives a recurrence for A r_k and its norm, which is zero at termination in both terminal-rank cases.Source: §6.1 · p. 10, §6.2 · p. 11

6.3. Matrix norms.

  • For 12 test matrices, the scheme using ∥T_j∥ with j = 20 gave significantly more accurate ∥A∥ estimates than the other tested schemes.The authors note that choosing j is not always clear and that the scheme adds some cost; their readily available estimate is satisfactory for estimating the order of ∥A∥.Source: §6.3 · p. 12

6.4. Matrix condition numbers.

  • The estimate κ(T_k) is a lower bound for κ(A).MINRES-QLP uses the QLP decomposition property to form this condition-number estimate.Source: §6.4 · p. 12

6.5. Solution norms.

  • MINRES-QLP computes a recurrence for ∥x_k∥ at a cost comparable to the norm of a three- or four-element vector.Because its basis is orthonormal, ∥x_k∥ can be estimated as ∥u_k∥; neither ∥x_k∥ nor its recurred estimate is guaranteed to increase monotonically.Source: §6.5 · pp. 12–13

6.6. Projection norms.

  • The projection norm ∥A x_k∥ can be updated by the recurrence ω_k = ∥[ω_k−1 τ_k]∥, starting with ω_0 = 0.Source: §6.6 · p. 13

7. Preconditioned MINRES and MINRES-QLP.

The section extends MINRES and MINRES-QLP to preconditioned symmetric systems while preserving symmetry through equivalent systems. Its derivations use a positive-definite preconditioner and adapt the Lanczos and QLP bases.

7. Preconditioned MINRES and MINRES-QLP.

  • Two-sided preconditioning preserves symmetry, but for incompatible systems it changes the least-squares norm and may change the minimum-length solution.The section therefore uses larger equivalent compatible systems for singular least-squares problems.Source: §7 · p. 13

7.1. Derivation.

  • Preconditioned MINRES is derived by applying MINRES to an equivalent system with matrix M^-1A.The passage introducing this derivation is incomplete, so it does not support further details of the transformed system.Source: §7.1 · p. 13

7.1.1. Preconditioned Lanczos process.

  • The preconditioned Lanczos iteration replaces the recurrence in Lanczos vectors with one involving consecutive z vectors and requires solving Mq_k = z_k each iteration.M is positive definite, so the square root used in the iteration is well-defined.Source: §7.1.1 · p. 14

7.1.2. Preconditioned MINRES.

  • Preconditioned MINRES updates the original-system solution using the transformed basis and preconditioned Lanczos quantities.The supplied passages give only fragments of the update recurrence.Source: §7.1.2 · p. 14

7.1.3. Preconditioned MINRES-QLP.

  • Preconditioned MINRES-QLP applies right reflectors to R_k and uses W_k = V_kP_k as its solution basis.For indefinite A, an approximate block-LDLT factorization can yield a positive-definite preconditioner M = L|D|L^T.Source: §7.1.3 · p. 14

7.2. Preconditioning singular Ax = b.

For singular compatible systems, preconditioned MINRES and MINRES-QLP return minimum-length solutions of the preconditioned system, but the corresponding original-system solution need not be minimum-length.

7.2. Preconditioning singular Ax = b.

  • If M is nonsingular, the preconditioned system remains compatible and the solvers return its minimum-length solution.This does not ensure the corresponding solution of Ax ≈ b is minimum-length.Source: §7.2 · p. 16
  • The unpreconditioned solution solves Ax ≈ b but is not necessarily the minimum-length solution.Source: §7.2 · p. 16
  • In the example, diagonal preconditioning produces a solution of Ax = b that differs from the stated minimum-length solution.The preconditioned solution is y† = [3.5739, 3.6819, 9.6909, 0.93156]^T, and x = Dy† = [3.0092, 2.9908, 3.0000, 3.0092]^T, whereas x† = [2, 4, 3, 2]^T.Source: §7.2 · p. 16

7.3. Preconditioning singular Ax ≈b.

For singular incompatible systems, the section proposes solving larger equivalent compatible systems to obtain minimum-residual solutions, while warning that preconditioning can change the norm in which a solution is minimum-length. It also surveys general preconditioner options and their scope.

7.3. Preconditioning singular Ax ≈b.

  • Larger equivalent compatible systems let MINRES produce minimum-residual solutions for singular incompatible problems, though the resulting solution may be minimum-length in a different norm.For ill-conditioned systems, the authors say applying MINRES-QLP to the original incompatible system is safer and more efficient.Source: §7.3 · p. 16

7.3.1. Augmented system.

  • The augmented system is compatible even when A is singular; its residual component is unique, but its x component may not be x†.Source: §7.3.1 · p. 16

7.3.2. A giant KKT system.

  • The KKT formulation is symmetric and compatible, and its nonsingular upper-left 3×3 block makes r, x, and y unique.The passage states that preconditioning this KKT system would give x as the minimum-length solution of the original problem.Source: §7.3.2 · p. 16

7.3.3. Regularization.

  • Regularization is proposed when A’s rank is ill-determined; the two-layered and KKT-like formulations have x → x† as δ → 0.LSQR reduces ∥r_k∥² + δ²∥x_k∥² monotonically, while MINRES or MINRES-QLP tend to reduce ∥Ar_k − δ²x_k∥ monotonically.Source: §7.3.3 · pp. 16–17

7.4. General preconditioners.

  • General diagonal and incomplete-Cholesky methods require access to A’s nonzero elements, so they do not apply when A is available only through products Ax.Source: §7.4 · p. 17

7.4.1. Diagonal preconditioning.

  • For diagonally dominant nonsingular A, diagonal scaling forms D with d_j = 1/|A_jj|^1/2 and solves DADy = Db, with x = Dy.The scaled matrix remains diagonally dominant and nonsingular, with entries no greater than 1 in magnitude.Source: §7.4.1 · p. 17

7.4.2. Binormalization (BIN).

  • In one example, binormalization gives κ(DAD) ≈ 2.6, while diagonal scaling with δ = 1 gives a condition number of approximately 10^4.Source: §7.4.2 · p. 18

7.4.3. Incomplete Cholesky factorization.

  • IC0 preserves the sparsity pattern for sparse symmetric positive-definite A, but may not exist when A is semidefinite or indefinite.An incomplete Cholesky-infinity variant may work for semidefinite or indefinite matrices.Source: §7.4.3 · p. 18

8. Numerical experiments.

The experiments compare Lanczos-based solvers with direct eigenvalue solutions across several symmetric and Hermitian problem types. Results include MINRES-QLP matching the TEVD solution in a clearly incompatible Laplacian least-squares case, alongside reported regularization and rounding-error behavior.

8. Numerical experiments.

  • The tests cover nonsingular and singular systems, mildly incompatible systems, symmetric least-squares problems, and Hermitian systems.The authors describe their examples as illustrating key features and refer to [10, Chapter 4] for a larger set of tests.Source: §8 · p. 19

8.1. A Laplacian system Ax ≈b (almost compatible).

  • For the almost-compatible Laplacian example, MINRES SOL gives a larger solution than MINRES-QLP, while LSQR reaches the minimum-length solution in more than twice MINRES-QLP’s iterations.MINRES-QLP runs 410 iterations in MINRES mode and 202 in MINRES-QLP mode; the residual norm is about 1.7 × 10^-8.Source: §8.1 · p. 20

8.2. A Laplacian LS problem min ∥Ax −b∥.

  • MINRES-QLP is the only solver matching TEVD in the clearly incompatible Laplacian least-squares example; Lanczos-based methods often match TEVD solutions.The example uses the same Laplacian with an incompatible right-hand side and residual norm about 17.Source: §8 · p. 19, §8.2 · p. 20

8.3. Regularizing effect of MINRES-QLP.

  • The stopping condition χ_k ≤ maxxnorm triggers regularization by treating a large solution component as zero, with later components also truncated if needed.The strategy relies on x_k = W_ku_k with W_k orthogonal.Source: §8.3 · pp. 20–21

8.4. Effects of rounding errors in MINRES-QLP.

  • In rounding-error tests, the MINRES residual estimate can keep decreasing after the directly computed residual stagnates, but the final solution may be less accurate than predicted.Except in very ill-conditioned least-squares problems, the MINRES-QLP estimates closely track the corresponding quantities.Source: §8.4 · p. 21

9. Conclusion.

MINRES-QLP combines QLP-based solution estimates with new norm recurrences and stopping conditions for singular systems. The conclusion highlights its recovery from MINRES instability, stronger rank-based regularization, and open needs in norm estimation and stopping thresholds.

9. Conclusion.

  • MINRES-QLP forms iterates through orthogonal steps and one triangular solve per iteration, unlike MINRES’s separate triangular systems for each direction.MINRES uses the recursion x_k = x_(k-1) + τ_k d_k, whereas MINRES-QLP uses x_k = (V_kP_k)u_k and solves L_ku_k = Q_k(β_1e_1).Source: §9 · p. 21
  • In the illustrated semidefinite case, MINRES’s solution explodes while MINRES-QLP returns to the TEVD solution at iteration 20, illustrating its ability to overcome potential instability.The paper describes the additional work and storage as moderate and transfers to MINRES-QLP only when the estimated condition number exceeds a specified value.Source: §9 · p. 21, Figure 8.1 · p. 22
  • Recurrences for ∥Ar_k∥ and ∥Ax_k∥ support new stopping conditions for singular problems and apply to Hermitian matrices.The paper also describes stronger regularization through MINRES-QLP’s more careful monitoring of the rank of T_k.Source: §9 · p. 22
  • MINRES-QLP’s solution approximations are nonmonotone, leaving robust a priori solution-norm bounds and a practical stopping threshold as research needs.The conclusion also identifies cases where γ_k^(4) is numerically small but not exactly zero as future work.Source: §9 · p. 22

10. Software and reproducible research.

The paper reports numerical comparisons of MINRES and MINRES-QLP and provides implementations and scripts for reproducibility. Across the described semidefinite examples, MINRES-QLP retains more reliable behavior in several ill-conditioned cases, although one stopping test does not shut it down promptly.

10. Software and reproducible research.

  • The authors provide MINRES implementations in Matlab 7.6, Fortran 77, and Fortran 90, a Matlab 7.6 MINRES-QLP implementation, and scripts at SOL.They state that each figure and example is accompanied by its source or a specific Matlab command.Source: §10 · p. 22
  • With b=Ae and η=10^-8, both recurred quantities reach 10^-13 at iteration 33; at η=10^-10, MINRES’s final residual is approximately 10^-13 at iteration 37.The η=10^-10 case is described as more ill-conditioned than the upper plots.Source: §10 · p. 23
  • With η=10^-10 and b=e, MINRES ends with residual norm approximately 10^2, while MINRES-QLP reaches 10^-4.The passage identifies the MINRES result as no convergence.Source: §10 · p. 23
  • For η=10^-2, directly computed ∥Ar_k∥ values match MINRES-QLP’s throughout iterations; at η=10^-4 MINRES’s final residual is slightly larger, and at η=10^-6 MINRES does not converge (ψM≈3).Source: §10 · p. 24
  • At η=10^-8, MINRES-QLP’s minimum ∥Ar_k∥ is approximately 10^-7, but tol=10^-8 does not stop it soon enough..Source: §10 · p. 24
Loading 1003.4042v3…