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

Introduction

We are concerned with iterative methods for solving a symmetric linear system $Ax=b$ or the related least-squares (LS) problem

The solution of (1), called the minimum-length or pseudoinverse solution , is formally given by $x^{\dagger}=(A^{T}\!A)^{\dagger}A^{T}\!b=(A^{2})^{\dagger}Ab=(A^{\dagger})^{2}Ab$, where $A^{\dagger}$ denotes the pseudoinverse of $A$. The pseudoinverse is continuous under perturbations $E$ for which $\mathop{\operator@font rank}\nolimits{(A+E)}=\mathop{\operator@font rank}\nolimits{(A)}$ , and $x^{\dagger}$ is continuous under the same condition. Problem (1) is then well-posed .

Let $A=U\Lambda U^{T}\!$ be an eigenvalue decomposition of $A$, with $U$ orthogonal and $\Lambda\equiv\mathop{\operator@font diag}\nolimits(\lambda_{1},\ldots,\lambda_{n})$. We define the condition number of $A$ to be $\kappa(A)=\smash[b]{\frac{\max|\lambda_{i}|}{\min_{\lambda_{i}\neq 0}|\lambda_{i}|}}$, and we say that $A$ is ill-conditioned if $\kappa(A)\gg 1$. Hence a singular matrix could be well-conditioned or ill-conditioned.

SYMMLQ and MINRES are Krylov subspace methods for solving symmetric indefinite systems $Ax=b$. SYMMLQ is reliable on compatible systems even if $A$ is ill-conditioned or singular, while on (singular) incompatible problems its iterates $x_{k}$ diverge to a multiple of a nullvector of $A$ [10, Proposition 2.15] and [10, Lemma 2.17]. MINRES seems more desirable to users because its residual norms are monotonically decreasing. On singular compatible systems, MINRES returns $x^{\dagger}$ (see Theorem 1). On singular incompatible systems, MINRES is reliable if terminated with a suitable stopping rule involving $\|Ar_{k}\|$ (see Lemma 3), but the solution will not be $x^{\dagger}$.

Here we develop a new solver of this type named MINRES-QLP . The aim is to deal reliably with compatible or incompatible systems and to return the unique solution of (1). We give theoretical reasons why MINRES-QLP improves the accuracy of MINRES on ill-conditioned systems, and illustrate with numerical examples.

Incompatible symmetric systems could arise from discretized semidefinite Neumann boundary value problems [27, section 4], and from any other singular systems involving measurement errors in $b$. Another potential application is large symmetric indefinite low-rank Toeplitz LS problems as described in [16, section 4.1].

2 Overview

In sections 2–4 we briefly review the Lanczos process, MINRES, and QLP decomposition before introducing MINRES-QLP in section 5. We derive norm estimates in section 6 and preconditioned MINRES-QLP in section 7. Numerical experiments are described in section 8.

The Lanczos process

Given $A$ and $b$, the Lanczos process computes vectors $v_{k}$ and tridiagonal matrices $\underline{T_{k}}$ according to $v_{0}\equiv 0$, $\beta_{1}v_{1}=b$, and thenNumerically, $p_{k}=Av_{k}-\beta_{k}v_{k-1}$, $\alpha_{k}=v_{k}^{T}\!p_{k}$, $\beta_{k+1}v_{k+1}=p_{k}-\alpha_{k}v_{k}$ is slightly better .

The $k$th Krylov subspace generated by $A$ and $b$ is defined to be $\mathcal{K}_{k}(A,b)=\mbox{\rm span}\{b,Ab,A^{2}b,\dots,A^{k-1}b\}=\mbox{\rm span}(V_{k})$. The following properties should be kept in mind:

If $A$ is changed to $A-\sigma I$ for some scalar shift $\sigma$, $T_{k}$ becomes $T_{k}-\sigma I$ and $V_{k}$ is unaltered, showing that singular systems are commonplace. Shifted problems appear in inverse iteration or Rayleigh quotient iteration.

MINRES

To make $r_{k}$ small, it is clear that $\beta_{1}e_{1}-\underline{T_{k}}y_{k}$ should be small. At this iteration $k$, MINRES minimizes the residual subject to $x_{k}\in\mathcal{K}_{k}(A,b)$ by choosing

This subproblem is processed by the expanding QR factorization: $Q_{0}\equiv 1$ and

where $c_{k}$ and $s_{k}$ form the Householder reflector $Q_{k,k+1}$ that annihilates $\beta_{k+1}$ in $\underline{T_{k}}$ to give upper-tridiagonal $R_{k}$, with $R_{k}$ and $t_{k}$ being unaltered in later iterations.

Also, $\tau_{k}=\phi_{k-1}c_{k}$ and $\phi_{k}=\phi_{k-1}s_{k}>0$. Hence from (3)–(5),

which is nonincreasing and tending to zero if $Ax=b$ is compatible.

2 Compatible systems

The following theorem assures us that MINRES is a useful solver for compatible linear systems even if $A$ is singular.

3 Incompatible systems

For a singular LS problem $Ax\approx b$, the optimal residual vector ${\widehat{r\mkern 3.0mu}\mkern-3.0mu}{}$ is unique, but infinitely many solutions $x$ give that residual. In the following example, MINRES finds a least-squares solution (with optimal residual) but not the minimum-length solution.

Let $A=\mathop{\operator@font diag}\nolimits(1,1,0)$ and $b=e$. The minimum-length solution to $Ax\approx b$ is $x^{\dagger}=[1\ 1\ 0]^{T}\!$ with residual ${\widehat{r\mkern 3.0mu}\mkern-3.0mu}{}=b-Ax^{\dagger}=e_{3}$ and $A{\widehat{r\mkern 3.0mu}\mkern-3.0mu}{}=0$. MINRES returns the solution $x^{\sharp}=e$ (with residual $r^{\sharp}=b-Ax^{\sharp}=e_{3}={\widehat{r\mkern 3.0mu}\mkern-3.0mu}{}$ and $Ar^{\sharp}=0$).

4 Norm estimates in MINRES

For incompatible systems, $r_{k}$ (3) will never be zero. However, all LS solutions satisfy $A^{2}x=Ab$, so that $Ar=0$. We therefore need a new stopping condition based on the size of $\|Ar_{k}\|$. In applications requiring nullvectors, $\|Ax_{k}\|$ is also useful. We present efficient recurrence relations for $\|Ar_{k}\|$ and $\|Ax_{k}\|$ in the following Lemma, which was not considered in the framework of MINRES when it was originally designed for nonsingular systems .

For the recurrence relations of $Ax_{k}$ and its norm, we have

Note that even using finite precision the expression for $\psi_{k}^{2}$ is extremely accurate for the versions of the Lanczos algorithm given in section 2, since (taking $\|v_{j}\|=1$ with negligible error), $\|Ar_{k}\|^{2}=\phi_{k}^{2}([\gamma_{k+1}]^{2}+2\gamma_{k+1}\delta_{k+2}v_{k+1}^{T}v_{k+2}+[\delta_{k+2}]^{2})$, where from (9) $|\delta_{k+2}|\leq\beta_{k+2}$, while from [38, (18)] $\beta_{k+2}|v_{k+1}^{T}v_{k+2}|\leq O(\varepsilon)\|A\|$, and with $|\gamma_{k+1}|\leq\|A\|$, see [38, (19)], we see that $|\gamma_{k+1}\delta_{k+2}v_{k+1}^{T}v_{k+2}|\leq O(\varepsilon)\|A\|^{2}$.

Typically $\|Ar_{k}\|$ is not monotonic, while clearly $\|r_{k}\|$ and $\|Ax_{k}\|$ are monotonic. In the eigensystem $A=U\Lambda U^{T}\!$, let $U={\begin{bmatrix}U_{1}\!&\!U_{2}\end{bmatrix}}$, where the eigenvectors $U_{1}$ correspond to nonzero eigenvalues. Then $P_{A}\equiv U_{1}U_{1}^{T}\!$ and $P^{\perp}_{A}\equiv U_{2}U_{2}^{T}\!$ are orthogonal projectors onto the range and nullspace of $A$. For general linear LS problems, Chang et al. characterize the dynamics of $\|r_{k}\|$ and $\|A^{T}r_{k}\|$ in three phases defined in terms of the ratios among $\|r_{k}\|$, $\|P_{A}r_{k}\|$, and $\|P^{\perp}_{A}r_{k}\|$, and propose two new stopping criteria for iterative solvers. The expositions show that these estimates are cheaply computable in CGLS and LSQR .

5 Effects of rounding errors in MINRES

Sleijpen, Van der Vorst, and Modersitzki analyzed the effects of rounding errors in MINRES and reported examples of apparent failure with a matrix of the form $A=QDQ^{T}\!$, where $D$ is an ill-conditioned diagonal matrix and $Q$ involves a single plane rotation. We were unable to reproduce MINRES’s performance on the two examples defined in Figure 4 of their paper, but we modified the examples by using an $n\times n$ Householder transformation for $Q$, and then observed similar difficulties with MINRES—see Figure 5. The recurred residual norm $\phi^{M}_{k}$ is a good approximation to the directly computed $\|r^{M}_{k}\|$ until the last few iterations. The recurred norms $\phi^{M}_{k}$ then keep decreasing but the directly computed norms $\|r^{M}_{k}\|$ become stagnant or even increase (see the lower subplots in Figure 5).

Note that we do want $\phi_{k}$ to keep decreasing on compatible systems, so that the test $\phi_{k}\leq\mathit{tol}(\|A\|\|x_{k}\|+\|b\|)$ with $\mathit{tol}\geq\varepsilon$ will eventually be satisfied even if the computed $\|r_{k}\|$ is no longer as small as $\phi_{k}$.

The analysis in focuses on the rounding errors involved in the $n$ lower triangular solves $R_{k}^{T}\!D_{k}^{T}\!=V_{k}^{T}\!$ (one solve for each row of $D_{k}$), compared to the single upper triangular solve $R_{k}y_{k}=t_{k}$ (followed by $x_{k}=V_{k}y_{k}$) that would be possible at the final $k$ if all of $V_{k}$ were stored as in GMRES . We shall see that a key feature of MINRES-QLP is that a single lower triangular solve suffices with no need to store $V_{k}$, much the same as in SYMMLQ.

Orthogonal decompositions for singular matrices

In 1999 Stewart proposed the pivoted QLP decomposition , which is equivalent to two consecutive QR factorizations with column interchanges, first on $A$, then on $R^{T}$:

giving nonnegative diagonal elements, where $\Pi_{R}$ and $\Pi_{L}$ are permutations chosen to maximize the next diagonal element of $R$ and $\hat{R}$ at each stage. This gives $A=QLP$, where

with $Q$ and $P$ orthogonal. Stewart demonstrated that the diagonals of $L$ (the $L$-values) give better singular-value estimates than the diagonals of $R$ (the $R$-values), and the accuracy is particularly good for the extreme singular values $\sigma_{1}$ and $\sigma_{n}$:

The first permutation $\Pi_{R}$ in pivoted QLP is important. The main purpose of the second permutation $\Pi_{L}$ is to ensure that the $L$-values present themselves in decreasing order, which is not always necessary. If $\Pi_{R}=\Pi_{L}=I$, it is simply called the QLP decomposition.

MINRES-QLP

and (4) is solved by $L_{k}u_{k}=t_{k}$ and $y_{k}=P_{k}u_{k}$. The MINRES-QLP estimate of $x$ is therefore $x_{k}=V_{k}y_{k}=V_{k}P_{k}u_{k}=W_{k}u_{k},$ with theoretically orthonormal $W_{k}\equiv V_{k}P_{k}$.

We will see that only the last three columns of $V_{k}$ are needed to update $x_{k}$.

The QLP decomposition of each $\underline{T_{k}}$ must be without permutations in order to ensure inexpensive updating of the factors as $k$ increases. Our experience is that the desired rank-revealing properties (16) tend to be retained, perhaps because of the tridiagonal structure of $\underline{T_{k}}$ and the convergence properties of the underlying Lanczos process.

For MINRES-QLP to be efficient, in the $k$th iteration ($k\geq 3$) the application of the left reflector $Q_{k,k+1}$ is followed immediately by the right reflectors $P_{k-2,k},P_{k-1,k}$, so that only the last $2\times 2$ principal submatrix of the transformed $\underline{T_{k}}$ will be changed in future iterations. These ideas can be understood more easily from Figure 1 and the following compact form, which represents the actions of right reflectors on $\underline{T_{k}}$ (additional to $Q_{k,k+1}$ (9)):

Figure 2 shows the relation between the singular values of $A$ and the diagonal elements of $R_{k}$ and $L_{k}$ with $k=19$. This illustrates (16) for matrix ID 1177 from with $n=25$.

3 Solving the subproblem

With $y_{k}=P_{k}u_{k}$, subproblem (4) becomes

where $t_{k}$ and $\phi_{k}$ are as in (5) and (8). At the start of iteration $k$, the first $k\!-\!3$ elements of $u_{k}$, denoted by $\mu_{j}$ for $j\leq k\!-\!3$, are known from previous iterations; see the 10th matrix in Figure 1. The remainder depend on the rank of $L_{k}$.

The corresponding solution estimate is $x_{k}=V_{k}y_{k}=V_{k}P_{k}u_{k}=W_{k}u_{k}$, where

and we update $x_{k-2}$ and compute $x_{k}$ by short-recurrence orthogonal steps:

4 Termination

5 Transfer from MINRES to MINRES-QLP

On well-conditioned systems, MINRES and MINRES-QLP behave very similarly. However, MINRES-QLP requires one more vector of storage, and each iteration needs 4 more axpy’s ($y\leftarrow\alpha x+y$) and 3 more vector scalings ($x\leftarrow\alpha x$). Thus it would be a desirable feature to invoke MINRES-QLP from MINRES only if $A$ is ill-conditioned or singular. The key idea is to transfer to MINRES-QLP at an iteration where $\underline{T_{k}}$ is not yet too ill-conditioned. The MINRES and MINRES-QLP solution estimates are the same, so from (6), (25), and (19): $x_{k}^{M}=x_{k}\Longleftrightarrow D_{k}t_{k}=W_{k}u_{k}=W_{k}L_{k}^{-1}t_{k}$. Now from (12), (17), and (22),

and the last three columns of $W_{k}$ can be obtained from the last three columns of $D_{k}$ and $L_{k}$. (Thus, we transfer the three MINRES basis vectors $d_{k-2},d_{k-1},d_{k}$ to $w_{k-2},w_{k-1},w_{k}$.) In addition, we need to generate $\smash{x_{k-2}^{(2)}}$ using (24):

It is clear from (26) that we still need to do the right transformation $R_{k}P_{k}=L_{k}$ in the MINRES phase and keep the last $3\times 3$ principal submatrix of $L_{k}$ for each $k$ so that we are ready to transfer to MINRES-QLP when necessary. We then obtain a short recurrence for $\|x_{k}\|$ (see section 6.5) and for this computation we save flops relative to the original MINRES algorithm, where $\|x_{k}\|$ is computed directly.

In the implementation, the MINRES iterates transfer to MINRES-QLP iterates when an estimate of the condition number of $T_{k}$ (see (29)) exceeds an input parameter $\mathit{trancond}$. Thus, $\mathit{trancond}>1/\varepsilon$ leads to MINRES iterates throughout, while $\mathit{trancond}=1$ generates MINRES-QLP iterates from the start.

6 Comparison of Lanczos-based solvers

We compare MINRES-QLP with CG, SYMMLQ, and MINRES in Tables 1–2 in terms of subproblem definitions, basis, solution estimates, flops and memory. A careful implementation of SYMMLQ provides a point in $\mathcal{K}_{k+1}(A,b)$ as shown. All solvers need storage for $v_{k}$, $v_{k+1}$, $x_{k}$, and a product $p_{k}=Av_{k}$ each iteration. Some additional work-vectors are needed for each method (e.g., $d_{k-1}$ and $d_{k}$ for MINRES, giving 7 work-vectors in total).

Stopping conditions and norm estimates

This section derives several norm estimates that are computed in MINRES-QLP. As before, we assume exact arithmetic throughout, so that $V_{k}$ and $Q_{k}$ are orthonormal. Table 3 summarizes how the norm estimates are used to formulate three groups of stopping conditions. The second NRBE test $\|Ar_{k}\|\leq\|A\|\|r_{k}\|\mathit{tol}$ is from Stewart with symmetric $A$.

First we derive recurrence relations for $r_{k}$ and its norm $\phi_{k}\equiv\|r_{k}\|$.

Next we derive recurrence relations for $Ar_{k}$ and its norm $\psi_{k}\equiv\|Ar_{k}\|$, and we show that $Ar_{k}$ is orthogonal to $\mathcal{K}_{k}(A,b)$.

For the first case, the proof is essentially the same as the proof of Lemma 3. For the other two cases, the results follow directly from Lemma 5. ∎

3 Matrix norms

For Lanczos-based algorithms, $\|A\|\geq\|V_{k+1}^{T}AV_{k}\|=\|\underline{T_{k}}\|$. Define

Then $\|A\|\geq\|\underline{T_{k}}\|\geq\mathcal{A}^{(k)}$. Clearly, $\mathcal{A}^{(k)}$ is monotonically increasing and is thus an improving estimate for $\|A\|$ as $k$ increases. By the property of QLP decomposition in (16) and (5.1), we could easily extend (27) to include the largest diagonal of $L_{k}$:

Some other schemes inspired by Larsen [31, section A.6.1], Higham , and Chen and Demmel follow. For the latter scheme, we use an implementation by Kaustuv for estimating the norms of the rows of $A$.

$\sqrt{\|\underline{T_{k}}^{T}\underline{T_{k}}\|_{1}}\geq\|T_{k}\|$

$\|T_{j}\|\leq\|T_{k}\|$ for small $j=5$ or $20$

Matlab function NORMEST$(A)$, which is based on the power method

Figure 3 plots estimates of $\|A\|$ for 12 matrices from the Florida sparse matrix collection whose sizes $n$ vary from 25 to 3002. In particular, scheme 3 above with $j=20$ gives significantly more accurate estimates than other schemes for the 12 matrices we tried. However, the choice of $j$ is not always clear and the scheme adds a little to the cost of MINRES-QLP. Hence we propose incorporating it into MINRES-QLP (or other Lanzcos-based iterative methods) if very accurate $\|A\|$ is needed. Otherwise (28) uses quantities readily available from MINRES-QLP and gives us satisfactory estimates for the order of $\|A\|$.

4 Matrix condition numbers

We again apply the property of the QLP decomposition in (16) and (5.1) to estimate $\kappa(\underline{T_{k}})$, which is a lower bound for $\kappa(A)$:

5 Solution norms

We derive a recurrence relation for $\|x_{k}\|$ whose cost is as low as computing the norm of a $3$- or $4$- vector.

Since $\|x_{k}\|=\|V_{k}P_{k}u_{k}\|=\|u_{k}\|$, we can estimate $\|x_{k}\|$ by computing $\chi_{k}\equiv\|u_{k}\|$. However, the last two elements of $u_{k}$ change in $u_{k+1}$ (and a new element $\mu_{k+1}$ is added). We therefore maintain $\chi_{k-2}$ by updating it and then using it according to

Thus $\chi_{k-2}^{(2)}$ increases monotonically but we cannot guarantee that $\|x_{k}\|$ and its recurred estimate $\chi_{k}$ are increasing, and indeed they are not in some examples (see Figure 4).

6 Projection norms

Sometimes the projection of the right-hand side vector $b$ onto $\mathcal{K}_{k}(A,b)$ is required (for example, see ). A simple recurrence relation is $\omega_{k}^{2}\equiv\|Ax_{k}\|^{2}=\omega_{k-1}^{2}+\tau_{k}^{2}$ and we can derive it in the same way as shown in Lemma 3. With $\omega_{0}\equiv 0$ we have $\omega_{k}\equiv\|Ax_{k}\|=\|[\omega_{k-1}\ \;\tau_{k}]\|$.

Preconditioned MINRES and MINRES-QLP

It is often asked: How can we construct a preconditioner for a linear system solver so that the same problem is solved with fewer iterations? Previous work on preconditioning the symmetric solvers CG, SYMMLQ, or MINRES includes .

We have the same question for singular symmetric systems $Ax\approx b$. Two-sided preconditioning is needed to preserve symmetry. We can still solve compatible systems, but we will no longer obtain the minimum-length solution. For incompatible systems, preconditioning alters the “least squares” norm. To avoid this difficulty we must work with larger equivalent systems that are compatible. We consider each case in turn, using a positive-definite preconditioner $M=CC^{T}$ with MINRES and MINRES-QLP to solve symmetric compatible systems $Ax=b$. Implicitly, we are solving equivalent symmetric systems $C^{-1}AC^{-T}y=C^{-1}b$, where $C^{T}\!x=y$. As usual, it is possible to work with $M$ itself, so without loss of generality we can assume $C=M^{\frac{1}{2}}$.

We derive preconditioned MINRES for compatible $Ax=b$ by applying MINRES to the equivalent problem $\bar{A}\bar{x}=\bar{b}$, where $\bar{A}\equiv M^{-\frac{1}{2}}AM^{-\frac{1}{2}}$, $\bar{b}\equiv M^{-\frac{1}{2}}b$, and $x=M^{-\frac{1}{2}}\bar{x}$.

Let $v_{k}$ denote the Lanczos vectors of $\mathcal{K}(\bar{A},\bar{b})$. With $v_{0}=0$ and $\beta_{1}v_{1}=\bar{b}$, for $k=1,2,\ldots$ we define

Then $\beta_{k}=\|\beta_{k}v_{k}\|=\|M^{-\frac{1}{2}}\!z_{k}\|=\|z_{k}\|_{M^{-1}}=\|q_{k}\|_{M}=\sqrt{q_{k}^{T}\!z_{k}},$ where the square root is well defined because $M$ is positive definite, and the Lanczos iteration is

Multiplying the last equation by $M^{\frac{1}{2}}$ we get

The last expression involving consecutive $z_{j}$’s replaces the three-term recurrence in $v_{j}$’s. In addition, we need to solve a linear system $Mq_{k}=z_{k}$ (30) each iteration.

1.2 Preconditioned MINRES

From (6) and (12) we have the following recurrence for the $k$th column of $D_{k}=V_{k}R_{k}^{-1}$ and $\bar{x}_{k}$:

Multiplying the above two equations by $M^{-\frac{1}{2}}$ on the left and defining $\bar{d}_{k}=M^{-\frac{1}{2}}d_{k}$, we can update the solution of our original problem by

We list the algorithm in [10, Table 3.4].

1.3 Preconditioned MINRES-QLP

A preconditioned MINRES-QLP can be derived very similarly. The additional work is to apply right reflectors $P_{k}$ to $R_{k}$, and the new subproblem bases are $W_{k}\equiv V_{k}P_{k}$, with $\bar{x}_{k}=W_{k}u_{k}$. Multiplying the new basis and solution estimate by $M^{-\frac{1}{2}}$ on the left, we obtain

Algorithm 1 lists all steps. Note that $\bar{w}_{k}$ is written as $w_{k}$ for all relevant $k$. Also, the output $x$ solves $Ax\approx b$ but the other outputs are associated with $\bar{A}\bar{x}\approx\bar{b}$.

The requirement of positive-definite preconditioners $M$ in MINRES and MINRES-QLP may seem unnatural for a problem with indefinite $A$ because we cannot achieve $M^{-\frac{1}{2}}\!AM^{-\frac{1}{2}}\approx I$. However, as shown in , we can achieve $M^{-\frac{1}{2}}\!AM^{-\frac{1}{2}}\approx\left[\begin{smallmatrix}I\\ &-I\end{smallmatrix}\right]$ using an approximate block-LDL${}^{\text{T}}$ factorization $A\approx LDL^{T}\!$ to get $M=L|D|L^{T}\!$, where $D$ is indefinite with blocks of order 1 and 2, and $|D|$ has the same eigensystem as $D$ except negative eigenvalues are changed in sign.

SQMR without preconditioning is analytically equivalent to MINRES. Unlike MINRES, SQMR can work directly with an indefinite preconditioner (such as block-LDL${}^{\text{T}}$). However, in finite precision, SQMR needs “look-ahead” to prevent numerical breakdown.

2 Preconditioning singular A​x=b𝐴𝑥𝑏Ax=b

For singular compatible systems, MINRES and MINRES-QLP find the minimum-length solution (see Theorems 1 and 4). If $M$ is nonsingular, the preconditioned system is also compatible and the solvers return its minimum-length solution. The unpreconditioned solution solves $Ax\approx b$, but is not necessarily a minimum-length solution.

Let $A=\left[\begin{smallmatrix}1&1&0&0\\ 1&1&1&0\\ 0&1&0&1\\ 0&0&1&0\end{smallmatrix}\right]$ and $b=\left[\begin{smallmatrix}6\\ 9\\ 6\\ 3\end{smallmatrix}\right].$ Then $\mathop{\operator@font rank}\nolimits(A)=3$ and $Ax=b$ is a singular compatible system. The minimum-length solution is $x^{\dagger}=\left[\begin{smallmatrix}2&4&3&2\end{smallmatrix}\right]^{T}\!$. By binormalization we construct the matrix $D=\mathop{\operator@font diag}\nolimits([\begin{smallmatrix}0.84201&0.81228&0.30957&3.2303\end{smallmatrix}])$. The minimum-length solution of the diagonally preconditioned problem $DADy\!=\!Db$ is $y^{\dagger}\!=\!\left[\begin{smallmatrix}3.5739&3.6819&9.6909&0.93156\end{smallmatrix}\right]^{T}\!$​​. Then $x=Dy^{\dagger}=\left[\begin{smallmatrix}3.0092&2.9908&3.0000&3.0092\end{smallmatrix}\right]^{T}\!$ is a solution of $Ax=b$, but $x\neq x^{\dagger}$.

3 Preconditioning singular A​x≈b𝐴𝑥𝑏Ax\approx b

We propose the following techniques for obtaining minimum-residual solutions of singular incompatible problems. In each case we use an equivalent but larger compatible system to which MINRES may be applied. Even if the larger system is singular, Theorem 1 shows that the minimum-length solution of the larger system will be obtained. The required $x$ will be part of this solution. Preconditioning still gives a minimum-residual solution of $Ax\approx b$, and in some cases $x$ will be $x^{\dagger}$. If the systems are ill-conditioned, it will be safer and more efficient to apply MINRES-QLP to the original incompatible system. However, preconditioning will give an $x$ that is “minimum length” in a different norm.

When $A$ is singular, so is the augmented system

but it is always compatible. Preconditioning with symmetric positive-definite $M$ gives us a solution $\left[\begin{smallmatrix}r\\ x\end{smallmatrix}\right]$ in which $r$ is unique, but $x$ may not be $x^{\dagger}$.

3.2 A giant KKT system

Problem (1) is equivalent to $\min_{r,\,x}x^{T}\!x$ subject to (31), which is an equality-contrained convex quadratic program. The corresponding KKT system [36, section 16.1] is both symmetric and compatible:

Although this is still a singular system, the upper-left $3\times 3$ block-submatrix is nonsingular and therefore $r$, $x$, and $y$ are unique and a preconditioner applied to the KKT system would give $x$ as the minimum-length solution of our original problem.

3.3 Regularization

If the rank of a given matrix $A$ is ill-determined, we may want to solve the regularized problem with parameter $\delta>0$:

The matrix $\left[\begin{smallmatrix}A\\ \delta I\end{smallmatrix}\right]$ has full rank and is always better conditioned than $A$. LSQR may be applied, and its iterates $x_{k}$ will reduce $\|r_{k}\|^{2}+\delta^{2}\|x_{k}\|^{2}$ monotonically. Alternatively, we could transform (33) into the following symmetric compatible systems and apply MINRES or MINRES-QLP. They tend to reduce $\|Ar_{k}-\delta^{2}x_{k}\|$ monotonically.

we obtain (34). Thus $x$ is also a solution of our regularized problem (33). This is equivalent to the two-layered formulation (4.3) in Bobrovnikova and Vavasis (with $A_{1}=A$, $A_{2}=D_{1}=D_{2}=I$, $b_{1}=b$, $b_{2}=0$, $\delta_{1}=1$, $\delta_{2}=\delta^{2}$). A key property is that $x\rightarrow x^{\dagger}$ as $\delta\rightarrow 0$.

If we define $y=-Av$ and $r=b-Ax-\delta^{2}y$, then we can show (by eliminating $r$ and $y$ from the following system) that $x$ in

is also a solution of (35) and thus of (33). The upper-left $3\times 3$ block-submatrix of (36) is nonsingular, and the correct limiting behavior occurs: $x\rightarrow x^{\dagger}$ as $\delta\rightarrow 0$. In fact, (36) reduces to (32).

4 General preconditioners

The construction of preconditioners is usually problem-dependent. If not much is known about the structure of $A$, we can only consider general methods such as diagonal preconditioning and incomplete Cholesky factorization. These methods require access to the nonzero elements of $A$. (They are not applicable if $A$ exists only as an operator for returning the product $Ax$.)

For a comprehensive survey of preconditioning techniques, see Benzi . We discuss a few methods for symmetric $A$ that also require access to the nonzero $A_{ij}$.

If $A$ has entries that are very different in magnitude, diagonal scaling might improve its condition. When $A$ is diagonally dominant and nonsingular, we can define $D=\mathop{\operator@font diag}\nolimits(d_{1},\ldots,d_{n})$ with $d_{j}=1/|A_{jj}|^{1/2}$. Instead of solving $Ax=b$, we solve $DADy=Db$, where $DAD$ is still diagonally dominant and nonsingular with all entries $\leq 1$ in magnitude, and $x=Dy$.

More generally, if $A$ is not diagonally dominant and possibly singular, we can safeguard division-by-zero errors by choosing a parameter $\delta>0$ and defining

If $A=\left[\begin{smallmatrix}-1&10^{-8}&\\ 10^{-8}&1&10^{4}\\ &10^{4}&0\\ &&&0\end{smallmatrix}\right]$, then $\kappa(A)\approx 10^{4}$. Let $\delta=1$, $D=\left[\begin{smallmatrix}1&&&\\ &10^{-2}&&\\ &&10^{-2}&\\ &&&1\end{smallmatrix}\right]$ in (37). Then $DAD=\left[\begin{smallmatrix}-1&10^{-10}&&\\ 10^{-10}&10^{-4}&1&\\ &1&0&\\ &&&0\end{smallmatrix}\right]$ and $\kappa(DAD)\approx 1$.

$A=\left[\begin{smallmatrix}10^{-4}&10^{-8}&&\\ 10^{-8}&10^{-4}&10^{-8}&\\ &10^{-8}&0&\\ &&&0\end{smallmatrix}\right]$ contains mostly very small entries, and $\kappa(A)\approx 10^{10}$. Let $\delta=10^{-8}$ and $D=\left[\begin{smallmatrix}10^{2}&&&\\ &10^{2}&&\\ &&10^{8}&\\ &&&10^{8}\end{smallmatrix}\right]$. Then $DAD=\left[\begin{smallmatrix}1&10^{-4}&&\\ 10^{-4}&1&10^{2}&\\ &10^{2}&0&\\ &&&0\end{smallmatrix}\right]$ and $\kappa(DAD)\approx 10^{2}$. (The choice of $\delta$ makes a critical difference in this case: with $\delta=1$, we have $D=I$.)

4.2 Binormalization (BIN)

Livne and Golub scale a symmetric matrix by a series of $k$ diagonal matrices on both sides until all rows and columns of the scaled matrix have unit $2$-norm: $DAD=D_{k}\cdots D_{1}AD_{1}\cdots D_{k}$. See also Bradley .

If $A=\left[\begin{smallmatrix}10^{-8}&1&\\ 1&10^{-8}&10^{4}\\ &10^{4}&0\end{smallmatrix}\right]$, then $\kappa(A)\approx 10^{12}$. With just one sweep of BIN, we obtain $D=\mathop{\operator@font diag}\nolimits(8.1\text{e}{-3},6.6\text{e}{-5},1.5)$, $DAD\approx\left[\begin{smallmatrix}6.5\text{e}{-1}&5.3\text{e}{-1}&0\\ 5.3\text{e}{-1}&0&1\\ 0&1&0\end{smallmatrix}\right]$ and $\kappa(DAD)\approx 2.6$ even though the rows and columns have not converged to one in the two-norm. In contrast, diagonal scaling (37) defined by $\delta=1$ and $D=\mathop{\operator@font diag}\nolimits(1,10^{-4},10^{-4})$ reduces the condition number to approximately $10^{4}$.

4.3 Incomplete Cholesky factorization

For a sparse symmetric positive definite matrix $A$, we could compute a preconditioner by the incomplete Cholesky factorization that preserves the sparsity pattern of $A$. This is known as IC0 in the literature. Often there exists a permutation $P$ such that the IC0 factor of $PAP^{T}\!$ is more sparse than that of $A$.

When $A$ is semidefinite or indefinite, IC0 may not exist, but a simple variant that may work is the incomplete Cholesky-infinity factorization [55, section 5].

Numerical experiments

We compare the computed results of MINRES-QLP and various other Krylov subspace methods to solutions computed directly by the eigenvalue decomposition (EVD) and the truncated eigenvalue decompositions (TEVD) of $A$. For TEVD we have

with parameter $t>0$. Often $t$ is set to $1$, and sometimes to a moderate number such as $10$ or $100$; it defines a cut-off point relative to the largest eigenvalue of $A$. For example, if most eigenvalues are of order 1 in magnitude and the rest are of order $\|A\|\varepsilon\approx 10^{-16}$, we expect TEVD to work better when the small eigenvalues are excluded, while EVD (with $t=0$) could return an exploding solution.

In the tables of results, Matlab MINRES and Matlab SYMMLQ are Matlab’s implementation of MINRES and SYMMLQ respectively. They incorporate local reorthogonalization of the Lanczos vector $v_{2}$, which could enhance the accuracy of the computations if $b$ is close to an eigenvector of $A$ :

Lacking the correct stopping condition for singular problems, Matlab SYMMLQ works more than necessary and then selects the smallest residual norm from all computed iterates; it would sometimes report that the method did not converge although the selected estimate appeared to be reasonably accurate.

MINRES SOL and SYMMLQ SOL are implementations based on . MINRES+ and SYMMLQ+ are the same but with additional stopping conditions for singular incompatible systems (see Lemma 3 and [10, Proposition 2.12]).

The computations in this section were performed on a Windows XP machine with a 3.2GHz Intel Pentium D Processor 940 and 3GB RAM ($\varepsilon\approx 10^{-16}$) . Tests were performed with each solver on five types of problem:

mildly incompatible symmetric (and singular) systems (meaning $\|r\|$ is rather small with respect to $\|b\|$)

We present a few examples that illustrate the key features of MINRES-QLP. For a larger set of tests and results, such as applying MINRES-QLP and other Krylov methods to Hermitian systems with preconditioning, we refer to [10, Chapter 4].

For a compatible system, we generate a random vector $b$ that is in the range of the test matrix ($b\equiv Ay$, $y_{i}\sim i.i.d.\ U(0,1)$, i.e., $y_{1},\ldots,y_{n}$ are independent and identically distributed random variables, whose values are drawn from the standard uniform distribution with support $$). For an LS problem, we generate a random$b$that is not in the range of the test matrix ($b_{i}\sim i.i.d.\ U(0,1)$ often suffices).

If $A$ is Hermitian, then $v^{H}\!Av$ is real for all complex vectors $v$. Numerically (in double precision), $\alpha_{k}=$ $v_{k}^{H}Av_{k}$ is likely to have small imaginary parts in the first few Lanczos iterations and snowball to have large imaginary parts in later iterations. This would result in a poor estimation of $\|T_{k}\|_{F}$ or $\|A\|_{F}$, and unnecessary errors in the Lanczos vectors. Thus we made sure to typecast $\alpha_{k}=$ $\operatorname{real}(v_{k}^{H}Av_{k})$ in MINRES-QLP and MINRES SOL.

We could say from the results that the Lanczos-based methods have built-in regularization features , often matching the TEVD solutions very well.

Our first example involves a singular indefinite Laplacian matrix $A$ of order $n=400$. It is block-tridiagonal with each block being a tridiagonal matrix $T$ of order $N=20$ with all nonzeros equal to 1:

Matlab’s eig($A$) reports the following data: $205$ positive eigenvalues in the interval $[6.1\text{e}{-2},8.87]$, $39$ almost-zero eigenvalues in $[-2.18\text{e}{-15},3.71\text{e}{-15}]$, $156$ negative eigenvalues in $[-2.91,-6.65\text{e}{-2}]$, numerical rank $=361$.

We used a right-hand side with a small incompatible component: $b=Ay+10^{-8}z$ with $y_{i}$ and $z_{i}\sim i.i.d.\ U(0,1)$. Results are summarized in Table 4. In the column labeled “C?”, the value “Y” denotes that the associated algorithm in the row has converged to the desired NRBE tolerances within $\mathit{maxit}$ iterations (cf. Table 3); otherwise, we have values “N” and “N?”, where “N?” indicates that the algorithm could have converged if more relaxed stopping conditions were used. The column “$Av$” shows the total number of matrix-vector products, and column “$x(1)$” lists the first element of the final solution estimate $x$ for each algorithm. For GMRES, the integer in parentheses is the value of the restart parameter.

MINRES SOL gives a larger solution than MINRES-QLP. This example has a residual norm of about $1.7\times 10^{-8}$, so it is not clear whether to classify it as a linear system or an LS problem. To the credit of Matlab SYMMLQ, it thinks the system is linear and returns a good solution. For MINRES-QLP, the first 410 iterations are in standard “MINRES mode”, with a transfer to “MINRES-QLP mode” for the last 202 iterations. LSQR converges to the minimum-length solution but with more than twice the number of iterations of MINRES-QLP. The other solvers fall short in some way.

2 A Laplacian LS problem min⁡‖A​x−b‖norm𝐴𝑥𝑏\min\|Ax-b\|

This example uses the same Laplacian matrix $A$ (38) but with a clearly incompatible $b=10\times\operatorname{rand}(n,1)$, i.e., $b_{i}\sim i.i.d.\ U(0,10)$. The residual norm is about $17$. Results are summarized in Table 5. MINRES gives an LS solution, while MINRES-QLP is the only solver that matches the solution of TEVD. The other solvers do not perform satisfactorily.

3 Regularizing effect of MINRES-QLP

This example illustrates the regularizing effect of MINRES-QLP with the stopping condition $\chi_{k}\leq\mathit{maxxnorm}$. For $k\geq 18$ in Figure 4, we observe the following values:

Because the last value exceeds $\mathit{maxxnorm}$, MINRES-QLP regards the last diagonal element of $L_{k}$ as a singular value to be ignored (in the spirit of truncated SVD solutions). It discards the last element of $u_{20}$ and updates

The full truncation strategy used in the implementation is justified by the fact that $x_{k}=W_{k}u_{k}$ with $W_{k}$ orthogonal. When $\|x_{k}\|$ becomes large, the last element of $u_{k}$ is treated as zero. If $\|x_{k}\|$ is still large, the second-to-last element of $u_{k}$ is treated as zero. If $\|x_{k}\|$ is still large, the third-to-last element of $u_{k}$ is treated as zero.

4 Effects of rounding errors in MINRES-QLP

The recurred residual norms $\phi_{k}^{M}$ in MINRES usually approximate the directly computed ones $\|r_{k}^{M}\|$ very well until $\|r_{k}^{M}\|$ becomes small. We observe that $\phi_{k}^{M}$ continues to decrease in the last few iterations, even though $\|r_{k}^{M}\|$ has become stagnant. This is desirable in the sense that the stopping rule will cause termination, although the final solution is not as accurate as predicted.

We present similar plots of MINRES-QLP in the following examples, with the corresponding quantities as $\phi_{k}^{Q}$ and $\|r_{k}^{Q}\|$. We observe that except in very ill-conditioned LS problems, $\phi_{k}^{Q}$ approximates $\|r_{k}^{Q}\|$ very closely.

Figure 5 illustrates four singular compatible linear systems.

Figure 6 illustrates four singular LS problems.

Conclusion

MINRES constructs its $k$th solution estimate from the recursion $x_{k}=D_{k}t_{k}=x_{k-1}+\tau_{k}d_{k}$ (6), where $n$ separate triangular systems $R_{k}^{T}\!D_{k}^{T}\!=V_{k}^{T}$ are solved to obtain the $n$ elements of each direction $d_{1},\ldots,d_{k}$. (Only $d_{k}$ is obtained during iteration $k$, but it has $n$ elements.)

In contrast, MINRES-QLP constructs its $k$th solution estimate using orthogonal steps: $x_{k}^{Q}=(V_{k}P_{k})u_{k}$; see (19)–(25). Only one triangular system $L_{k}u_{k}=Q_{k}(\beta_{1}e_{1})$ is involved for each $k$.

Thus MINRES-QLP overcomes the potential instability predicted by the MINRES authors and analyzed by Sleijpen et al. . The additional work and storage are moderate, and maximum efficiency is retained by transferring from MINRES to the MINRES-QLP iterates only when the estimated condition number of $A$ exceeds a specified value.

MINRES and MINRES-QLP are readily applicable to Hermitian matrices, once $\alpha_{k}$ is typecast as a real scalar in finite-precision arithmetic. For both algorithms, we derived recurrence relations for $\|Ar_{k}\|$ and $\|Ax_{k}\|$ and used them to formulate new stopping conditions for singular problems.

TEVD or TSVD are commonly known to use rank-$k$ approximations to $A$ to find approximate solutions to $\min\|Ax-b\|$ that serve as a form of regularization. Krylov subspace methods also have regularization properties . Since MINRES-QLP monitors more carefully the rank of $T_{k}$, which could be $k$ or $k\!-\!1$, we may say that regularization is a stronger feature in MINRES-QLP, as we have shown in our numerical examples.

It is important to develop robust techniques for estimating an a priori bound for the solution norm since the MINRES-QLP approximations are not monotonic as is the case in CG and LSQR. Ideally, we would also like to determine a practical threshold associated with the stopping condition $\gamma_{k}^{(4)}=0$ in order to handle cases when $\gamma_{k}^{(4)}$ is numerically small but not exactly zero. These are topics for future research.

Software and reproducible research

Matlab 7.6, Fortran 77, and Fortran 90 implementations of MINRES with new stopping conditions $\|Ar_{k}\|\!\leq\!\mathit{tol}\|A\|\|r_{k}\|$ and $\|Ax_{k}\|\leq\mathit{tol}\|A\|\|x_{k}\|$, and a Matlab 7.6 implementation of MINRES-QLP are available from SOL .

Following the philosophy of reproducible computational research as advocated in , for each figure and example in this paper we mention either the source or the specific Matlab command. Our Matlab scripts are available at SOL .

We thank Jan Modersitzki, Gerard Sleijpen, Henk Van der Vorst, and also Kaustuv for providing us with their Matlab scripts, which have aided us in producing parts of Figures 3, 5, and 6. We also thank Michael Friedlander, Rasmus Larsen, and Lek-Heng Lim for their finest examples of work, discussion, and support. We are most grateful to our anonymous reviewers for their insightful suggestions for improving this manuscript. Last but not least, we dedicate this paper to the memory of our colleague and friend, Gene Golub.

References