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

MINRES-QLP: A KR YLO V SUBSP A CE METHOD F OR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS ∗ SOU-CHENG T. CHOI † , CHRISTOPHER C. P AIGE ‡ , AND MICHAEL A. SAUNDERS § Abstract. CG, SYMMLQ, and MINRES are Krylov subspace methods for solving symmetric systems of linear equations. When these metho ds 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 w ould giv e a least-squares solution but not necessarily the minimum-length (pseudoinv erse) solution. This understanding motiv ates 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 pro cess (where R is upper- tridiagonal). MINRES-QLP uses a QLP decomposition (where rotations on the right reduce R to low er-tridiagonal form). On ill-conditioned systems (singular or not), MINRES-QLP can give more accurate solutions than MINRES. W e derive preconditioned MINRES-QLP , new stopping rules, and better estimates of the solution and residual norms, the matrix norm, and the condition n umber. Key words. MINRES, Krylo v subspace method, Lanczos pro cess, conjugate-gradient method, minimum-residual metho d, singular least-squares problem, sparse matrix AMS sub ject classiﬁcations. 15A06, 65F10, 65F20, 65F22, 65F25, 65F35, 65F50, 93E24 DOI. xxx/xxxxxxxxx 1. In tro duction. W e are concerned with iterative methods for solving a sym- metric linear system Ax = b or the related least-squares (LS) problem min k x k 2 s.t. x ∈ arg min x k Ax − b k 2 , (1.1) where A ∈ R n × n is symmetric and p ossibly singular, b ∈ R n , A 6 = 0, and b 6 = 0. Most of the results in our discussion are directly extendable to problems with complex Hermitian matrices A and complex vectors b . The solution of ( 1.1 ), called the minimum-length or pseudoinverse solution [ 18 ], is formally giv en b y x † = ( A T A ) † A T b = ( A 2 ) † Ab = ( A † ) 2 Ab , where A † denotes the pseudoin verse of A . The pseudoinv erse is con tinuous under perturbations E for which rank ( A + E ) = rank ( A ) [ 49 ], and x † is contin uous under the same condition. Problem ( 1.1 ) is then w ell-p osed [ 19 ]. Let A = U Λ U T b e an eigenv alue decomp osition of A , with U orthogonal and Λ ≡ diag( λ 1 , . . . , λ n ). W e deﬁne the condition num b er of A to b e κ ( A ) = max | λ i | min λ i 6 =0 | λ i | , and w e say that A is ill-conditioned if κ ( A )  1. Hence a singular matrix could b e w ell-conditioned or ill-conditioned. ∗ Received by the editors March 7, 2010. Draft MINRESQLP56 of Marc h 30, 2015. Revised March 31, 2011. http://www.siam.org/journals/sisc/ † Institute for Computational and Mathematical Engineering, Stanford University , Stanford, CA 94305-4121 (scc hoi@stanford.edu). This author’s research was partially supp orted b y National Sci- ence F oundation grant CCR-0306662. ‡ Computer Science, McGill Universit y , Mon treal, Quebec, Canada, H3A 2A7 (paige@cs.mcgill.ca). This author’s research was partially supp orted by NSERC of Canada grant OGP0009236. § Department of Management Science and Engineering, Stanford Universit y , Stanford, CA 94305- 4026 (saunders@stanford.edu). This author’s researc h was partially supp orted b y National Science F oundation gran t CCR-0306662, Oﬃce of Nav al Researc h grants N00014-02-1-0076 and N00014-08- 1-0191, and AHPCRC. 1 2 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS SYMMLQ and MINRES [ 39 ] are Krylov subspace metho ds for solving symmetric indeﬁnite 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 div erge to a multiple of a nullv ector of A [ 10 , Prop osition 2.15] and [ 10 , Lemma 2.17]. MINRES seems more desirable to users b ecause its residual norms are monotonically decreasing. On singular compatible systems, MINRES returns x † (see Theorem 3.1 ). On singular incompatible systems, MINRES is reliable if terminated with a suitable stopping rule inv olving k Ar k k (see Lemma 3.3 ), but the solution will not b e x † . Here w e dev elop a new solver of this type named MINRES-QLP [ 10 ]. The aim is to deal reliably with compatible or incompatible systems and to return the unique solution of ( 1.1 ). W e give theoretical reasons wh y MINRES-QLP impro ves the accuracy of MINRES on ill-conditioned systems, and illustrate with numerical examples. Incompatible symmetric systems could arise from discretized semideﬁnite Neu- mann b oundary v alue problems [ 27 , section 4], and from any other singular systems in volving measurement errors in b . Another potential application is large symmetric indeﬁnite low-rank T o eplitz LS problems as described in [ 16 , section 4.1]. 1.1. Notation. The letters i , j , k denote integer indices, c and s cosine and sine of some angle θ , e k the k th unit vector, e a v ector of all ones, and other low er-case letters such as b , u , and x (p ossibly with in teger subscripts) denote c olumn vectors. Upp er-case letters A , T k , V k , . . . denote matrices, and I k is the iden tity matrix of order k . Low er-case Greek letters denote scalars; in particular, ε ≈ 10 − 16 denotes the ﬂoating-p oin t precision. If a quan tit y δ k is mo diﬁed one or more times, we denote its v alues b y δ k , δ (2) k , δ (3) k , . . . . The sym b ol k · k denotes the 2-norm of a vector or matrix. F or an incompatible system, Ax ≈ b is shorthand for the LS problem min x k Ax − b k . 1.2. Ov erview. In sections 2 – 4 we brieﬂy review the Lanczos process, MINRES , and QLP decomp osition b efore introducing MINRES-QLP in section 5 . W e deriv e norm estimates in section 6 and preconditioned MINRES-QLP in section 7 . Numerical exp erimen ts are describ ed in section 8 . 2. The Lanczos pro cess. Given A and b , the Lanczos process [ 30 ] computes v ectors v k and tridiagonal matrices T k according to v 0 ≡ 0, β 1 v 1 = b , and then 1 p k = Av k , α k = v T k p k , β k +1 v k +1 = p k − α k v k − β k v k − 1 for k = 1 , 2 , . . . , ` , where w e choose β k > 0 to give k v k k = 1. In matrix form, AV k = V k +1 T k , T k ≡          α 1 β 2 β 2 α 2 . . . . . . . . . β k β k α k β k +1          ≡  T k β k +1 e T k  , V k ≡  v 1 · · · v k  . (2.1) In exact arithmetic, the columns of V k are orthonormal and the pro cess stops with k = ` and β ` +1 = 0 for some ` ≤ n , and then AV ` = V ` T ` . F or deriv ation purposes w e assume that this happ ens, though in practice it is unlikely unless V k is reorthog- onalized for each k . In any case, ( 2.1 ) holds to machine precision and the computed v ectors satisfy k V k k 1 ≈ 1 (even if k  n ). 1 Numerically , p k = Av k − β k v k − 1 , α k = v T k p k , β k +1 v k +1 = p k − α k v k is slightly b etter [ 38 ]. MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 3 2.1. Prop erties of the Lanczos process. The k th Krylo v subspace generated b y A and b is deﬁned to b e K k ( A, b ) = span { b, Ab, A 2 b, . . . , A k − 1 b } = span( V k ). The follo wing prop erties should b e kept in mind: 1. If A is c hanged to A − σ I for some scalar shift σ , T k b ecomes T k − σ I and V k is unaltered, showing that singular systems are commonplace. Shifted problems app ear in inv erse iteration or Ra yleigh quotient iteration. 2. T k has full column rank k for all k < ` . 3. If A is indeﬁnite, some T k migh t b e singular for k < ` , but then b y the Sturm sequence prop ert y (see [ 18 ]), T k has exactly one zero eigen v alue and the strict in terlacing prop ert y implies that T k ± 1 are nonsingular. Hence T k cannot b e singular twice in a row (whether A is singular or not). 4. T ` is nonsingular if and only if b ∈ range( A ). (See app endix A .) 3. MINRES . Algorithm MINRES [ 39 ] is a natural w ay of using the Lanczos pro cess to solv e Ax = b or min x k Ax − b k . F or k < ` , if x k = V k y k for some vector y k , the asso ciated residual is r k ≡ b − Ax k = b − AV k y k = β 1 v 1 − V k +1 T k y k = V k +1 ( β 1 e 1 − T k y k ) . (3.1) T o make r k small, it is clear that β 1 e 1 − T k y k should be small. A t this iteration k , MINRES minimizes the residual sub ject to x k ∈ K k ( A, b ) by choosing y k = arg min y ∈ R k k T k y − β 1 e 1 k . (3.2) This subproblem is pro cessed b y the expanding QR factorization: Q 0 ≡ 1 and Q k,k +1 ≡  I k − 1 c k s k s k − c k  , Q k ≡ Q k,k +1  Q k − 1 1  , Q k  T k β 1 e 1  =  R k t k 0 φ k  , (3.3) where c k and s k form the Householder reﬂector Q k,k +1 that annihilates β k +1 in T k to giv e upp er-tridiagonal R k , with R k and t k b eing unaltered in later iterations. When k < ` , the unique solution of ( 3.2 ) satisﬁes R k y k = t k . Instead of solving for y k , MINRES solves R T k D T k = V T k b y forward substitution, obtaining the last column d k of D k at iteration k . At the same time, it up dates x k via x 0 ≡ 0 and x k = V k y k = D k R k y k = D k t k = x k − 1 + τ k d k , τ k ≡ e T k t k . (3.4) When k = ` , we can form T ` but nothing else expands. In place of ( 3.1 ) and ( 3.3 ) we hav e r ` = V ` ( β 1 e 1 − T ` y ` ) and Q ` − 1  T ` β 1 e 1  =  R ` t `  and it is natural to choose y ` from the subproblem min k T ` y ` − β 1 e 1 k ≡ min k R ` y ` − t ` k . (3.5) There are tw o cases to consider: 1. If T ` is nonsingular, R ` y ` = t ` has a unique solution. Since AV ` y ` = V ` T ` y ` = b , the problem is solv ed b y x ` = V ` y ` with residual r ` = 0 (the system is compatible, even if A is singular). Theorem 3.1 pro v es that x ` = x † . 2. If T ` is singular, A and R ` are singular ( R `` = 0) and b oth Ax = b and R ` y ` = t ` are incompatible. This case was not handled b y MINRES in [ 39 ]. Theorem 3.2 pro ves that the MINRES p oin t x ` − 1 is a least-squares solution (but not necessarily x † ). Theorem 5.1 prov es that the MINRES-QLP p oin t x ` = V ` y † ` = x † , where y † ` is the min-length solution of ( 3.5 ). 4 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS 3.1. F urther details of MINRES . T o describ e MINRES-QLP thoroughly , w e need further details of the MINRES QR factorization ( 3.3 ). F or 1 ≤ k < ` ,  R k 0  =              γ 1 δ 2  3 γ (2) 2 δ (2) 3 . . . . . . . . .  k . . . δ (2) k γ (2) k 0              ,  t k φ k  ≡            τ 1 τ 2 . . . . . . τ k φ k            = β 1            c 1 s 1 c 2 . . . . . . s 1 · · · s k − 1 c k s 1 · · · s k − 1 s k            (3.6) (where the sup erscripts are deﬁned in section 1.1 ). With φ 0 ≡ β 1 > 0, the full action of Q k,k +1 in ( 3.3 ), including its eﬀect on later columns of T j , k < j ≤ ` , is described b y  c k s k s k − c k  γ k δ k +1 0 β k +1 α k +1 β k +2     φ k − 1 0  = " γ (2) k δ (2) k +1  k +2 0 γ k +1 δ k +2     τ k φ k # , (3.7) where s k = β k +1 / k [ γ k β k +1 ] k > 0, giving γ 1 , γ (2) k > 0 with R j nonsingular for eac h j ≤ k < ` . Th us the d j in ( 3.4 ) can b e found from R T k D T k = V T k : ( d 1 = v 1 /γ 1 , d 2 = ( v 2 − δ 2 d 1 ) /γ (2) 2 , d j = ( v j − δ (2) j d j − 1 −  j d j − 2 ) /γ (2) j , j = 3 , . . . , k . (3.8) Also, τ k = φ k − 1 c k and φ k = φ k − 1 s k > 0. Hence from ( 3.1 )–( 3.3 ), k r k k = k T k y k − β 1 e 1 k = φ k ⇒ k r k k = k r k − 1 k s k , (3.9) whic h is nonincreasing and tending to zero if Ax = b is compatible. Remark 3.1. If k < ` and T k is singular, we have γ k = 0 , s k = 1 , and k r k k = k r k − 1 k (not a strict de cr e ase), but this c annot happ en twic e in a r ow (cf. se ction 2.1 ). Remark 3.2. If T ` is singular, MINRES sets th e last element of y ` to b e zer o. The ﬁnal p oint and r esidual stay as x ` − 1 and r ` − 1 with k r ` − 1 k = φ ` − 1 = β 1 s 1 · · · s ` − 1 > 0 . 3.2. Compatible systems. The following theorem assures us that MINRES is a useful solver for compatible linear systems ev en if A is singular. Theorem 3.1 ( [ 10 , Theorem 2.25] ). If b ∈ range( A ) , the ﬁnal MINRES p oint x ` is the minimum-length solution of Ax = b (and r ` = b − Ax ` = 0 ). Pr o of . If b ∈ range( A ), the Lanczos pro cess gives AV ` = V ` T ` with nonsingular T ` , and MINRES terminates with Ax ` = b and x ` = V ` y ` = Aq , where q = V ` T − 1 ` y ` . If some other p oin t b x satisﬁes A b x = b , let p = b x − x ` . W e ha ve Ap = 0 and x T ` p = q T Ap = 0. Hence k b x k 2 = k x ` + p k 2 = k x ` k 2 + 2 x T ` p + k p k 2 ≥ k x ` k 2 . 3.3. Incompatible systems. F or a singular LS problem Ax ≈ b , the optimal residual v ector b r is unique, but inﬁnitely many solutions x give that residual. In the follo wing example, MINRES ﬁnds a least-squares solution (with optimal residual) but not the minimum-length solution. Example 3.1. L et A = diag(1 , 1 , 0) and b = e . The minimum-length solution to Ax ≈ b is x † = [1 1 0] T with r esidual b r = b − Ax † = e 3 and A b r = 0 . MINRES r eturns the solution x ] = e (with r esidual r ] = b − Ax ] = e 3 = b r and Ar ] = 0 ). Theorem 3.2 ( [ 10 , Theorem 2.27] ). If b / ∈ range( A ) , then k Ar ` − 1 k = 0 and the MINRES x ` − 1 is an LS solution (but not ne c essarily x † ). Pr o of . Since b / ∈ range( A ), T ` is singular and R `` = γ ` = 0. By Lemma 3.3 below, A ( Ax ` − 1 − b ) = − Ar ` − 1 = −k r ` − 1 k γ ` v ` = 0. Thus x ` − 1 is an LS solution. MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 5 3.4. Norm estimates in MINRES . F or incompatible systems, r k ( 3.1 ) will nev er b e zero. Ho wev er, all LS solutions satisfy A 2 x = Ab , so that Ar = 0. W e therefore need a new stopping condition based on the size of k Ar k k . In applications requiring nullv ectors, k Ax k k is also useful. W e present eﬃcient recurrence relations for k Ar k k and k Ax k k in the following Lemma, whic h w as not considered in the framework of MINRES when it was originally designed for nonsingular systems [ 39 ]. Lemma 3.3 ( Ar k , k Ar k k , k Ax k k for MINRES ). If k < ` , Ar k = k r k k ( γ k +1 v k +1 + δ k +2 v k +2 ) ( wher e δ k +2 v k +2 = 0 if k = ` − 1) , ψ 2 k ≡ k Ar k k 2 = k r k k 2  [ γ k +1 ] 2 + [ δ k +2 ] 2  ( wher e δ k +2 = 0 if k = ` − 1) , ω 2 k ≡ k Ax k k 2 = ω 2 k − 1 + τ 2 k , ω 0 ≡ 0 . Pr o of . F or k < ` , R k is nonsingular. F rom ( 3.1 )–( 3.4 ) with R k y k = t k w e hav e r k = V k +1 Q T k  t k φ k  −  R k 0  y k  = φ k V k +1 Q T k e k +1 , (3.10) Ar k = φ k V k +2 T k +1 Q T k e k +1 , Q k T k +1 T = Q k  T k +1 β k +2 e k +1  = Q k  T k β k +1 e k 0 β k +1 e T k α k +1 β k +2  , e T k +1 Q k T k +1 T =  0 γ k +1 δ k +2  , see ( 3.7 ), where AV k +1 = V k +1 T k +1 and we take δ k +2 = 0 if k = ` − 1, so Ar k = φ k V k +2  0 γ k +1 δ k +2  T = φ k ( γ k +1 v k +1 + δ k +2 v k +2 ) , ψ 2 k ≡ k Ar k k 2 = k r k k 2  [ γ k +1 ] 2 + [ δ k +2 ] 2  . F or the recurrence relations of Ax k and its norm, w e ha v e Ax k = AV k y k = V k +1 T k y k = V k +1 Q T k  R k 0  y k = V k +1 Q T k  t k 0  , ω 2 k ≡ k Ax k k 2 = k t k k 2 = k t k − 1 k 2 + τ 2 k = ω 2 k − 1 + τ 2 k . Note that even using ﬁnite precision the expression for ψ 2 k is extremely accurate for the versions of the Lanczos algorithm giv en in section 2 , since (taking k v j k = 1 with negligible error), k Ar k k 2 = φ 2 k ([ γ k +1 ] 2 + 2 γ k +1 δ k +2 v T k +1 v k +2 + [ δ k +2 ] 2 ), where from ( 3.7 ) | δ k +2 | ≤ β k +2 , while from [ 38 , (18)] β k +2 | v T k +1 v k +2 | ≤ O ( ε ) k A k , and with | γ k +1 | ≤ k A k , see [ 38 , (19)], w e see that | γ k +1 δ k +2 v T k +1 v k +2 | ≤ O ( ε ) k A k 2 . T ypically k Ar k k is not monotonic, while clearly k r k k and k Ax k k ar e monotonic. In the eigensystem A = U Λ U T , let U =  U 1 U 2  , where the eigenv ectors U 1 cor- resp ond to nonzero eigenv alues. Then P A ≡ U 1 U T 1 and P ⊥ A ≡ U 2 U T 2 are orthogonal pro jectors [ 53 ] onto the range and n ullspace of A . F or general linear LS problems, Chang et al. [ 7 ] c haracterize the dynamics of k r k k and k A T r k k in three phases deﬁned in terms of the ratios among k r k k , k P A r k k , and k P ⊥ A r k k , and propose tw o new stop- ping criteria for iterative solv ers. The exp ositions [ 1 , 26 ] show that these estimates are cheaply computable in CGLS and LSQR [ 40 , 41 ]. 6 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS 3.5. Eﬀects of rounding errors in MINRES . MINRES should stop if R k is singular (which theoretically implies k = ` and A is singular). Singularity was not discussed b y P aige and Saunders [ 39 ], but they did raise the question: Is MINRES stable when R k is ill-conditioned? Their concern was that k D k k could b e large in ( 3.8 ), and there could then b e cancellation in forming x k − 1 + τ k d k in ( 3.4 ). Sleijp en, V an der V orst, and Mo dersitzki [ 47 ] analyzed the eﬀects of rounding errors in MINRES and rep orted examples of apparen t failure with a matrix of the form A = QD Q T , where D is an ill-conditioned diagonal matrix and Q in volv es a single plane rotation. W e were unable to repro duce MINRES ’s p erformance on the t wo examples deﬁned in Figure 4 of their pap er, but w e mo diﬁed the examples by using an n × n Householder transformation for Q , and then observ ed similar diﬃculties with MINRES —see Figure 8.2 . The recurred residual norm φ M k is a go o d appro ximation to the directly computed k r M k k until the last few iterations. The recurred norms φ M k then k eep decreasing but the directly computed norms k r M k k b ecome stagnan t or ev en increase (see the lo wer subplots in Figure 8.2 ). Remark 3.3. Note that we do want φ k to ke ep de cr e asing on c omp atible systems, so that the test φ k ≤ tol ( k A kk x k k + k b k ) with tol ≥ ε wil l eventual ly b e satisﬁe d even if the c ompute d k r k k is no longer as smal l as φ k . The analysis in [ 47 ] fo cuses on the rounding errors in volv ed in the n low er tri- angular solv es R T k D T k = V T k (one solv e for eac h ro w of D k ), compared to the single upp er triangular solv e R k y k = t k (follo wed by x k = V k y k ) that w ould be possible at the ﬁnal k if all of V k w ere stored as in GMRES [ 44 ]. W e shall see that a key feature of MINRES-QLP is that a single low er triangular solv e suﬃces with no need to store V k , muc h the same as in SYMMLQ . 4. Orthogonal decomp ositions for singular matrices. In 1999 Stew art pro- p osed the pivote d QLP de c omp osition [ 51 ], which is equiv alent to tw o consecutive QR factorizations with column in terchanges, ﬁrst on A , then on R T : Q R A Π R =  R S 0 0  , Q L  R T 0 S T 0  Π L =  ˆ R 0 0 0  , (4.1) giving nonnegative diagonal elements, where Π R and Π L are p erm utations c hosen to maximize the next diagonal element of R and ˆ R at each stage. This giv es A = QLP , where Q = Q T R Π L , L =  ˆ R T 0 0 0  , P = Q L Π T R , with Q and P orthogonal. Stew art demonstrated that the diagonals of L (the L - values ) give b etter singular-v alue estimates than the diagonals of R (the R -values ), and the accuracy is particularly go od for the extreme singular v alues σ 1 and σ n : R ii ≈ σ i , L ii ≈ σ i , σ 1 ≥ max i L ii ≥ max i R ii , min i R ii ≥ min i L ii ≥ σ n . (4.2) The ﬁrst p erm utation Π R in piv oted QLP is imp ortan t. The main purp ose of the second p erm utation Π L is to ensure that the L -v alues presen t themselv es in decreasing order, which is not alw ays necessary . If Π R = Π L = I , it is simply called the QLP de c omp osition . MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 7 5. MINRES-QLP . W e now develop MINRES-QLP for solving ill-conditioned or singular symmetric systems Ax ≈ b . The Lanczos framew ork is the same as in MIN- RES , but w e handle T ` in ( 3.5 ) with extra care when it is rank-deﬁcien t. In this case, the normal approach to solving ( 3.5 ) is via a QLP decomposition of T ` to obtain the (unique) minim um-length solution y ` [ 51 , 18 ]. Thus in MINRES-QLP we use a QLP decomp osition of T k in subproblem ( 3.2 ) for all k ≤ ` . This is the MINRES QR ( 3.3 ) follo wed b y an LQ factorization of R k : Q k T k =  R k 0  , R k P k = L k , so that Q k T k P k =  L k 0  , (5.1) where Q k and P k are orthogonal, R k is upper tridiagonal and L k is lo wer tridiagonal. When k < ` , R k and L k are nonsingular. MINRES-QLP obtains the same solution as MINRES , but by a diﬀerent pro cess (and with diﬀeren t rounding errors). Deﬁning u b y y = P k u , we see from ( 3.3 ) that Q k ( T k y − β 1 e 1 ) =  L k 0  u −  t k φ k  , and ( 3.2 ) is solved b y 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 ≡ V k P k . W e will see that only the last three columns of V k are needed to up date x k . 5.1. The QLP factorization of T k . The QLP decomp osition of each T k m ust b e without permutations in order to ensure inexp ensiv e up dating of the factors as k increases. Our experience is that the desired rank-revealing prop erties ( 4.2 ) tend to b e retained, p erhaps b ecause of the tridiagonal structure of T k and the conv ergence prop erties of the underlying Lanczos pro cess. F or k < ` , the QLP decomp osition of T k ( 5.1 ) giv es nonsingular tridiagonal R k and L k . As in MINRES , Q k is a pro duct of Householder reﬂectors, see ( 3.3 ) and ( 3.7 ), while P k in volv es a pro duct of p airs of essentially 2 × 2 reﬂectors: Q k = Q k,k +1 · · · Q 3 , 4 Q 2 , 3 Q 1 , 2 , P k = P 1 , 2 P 1 , 3 P 2 , 3 · · · P k − 2 ,k P k − 1 ,k . F or MINRES-QLP to be eﬃcient, in the k th iteration ( k ≥ 3) the application of the left reﬂector Q k,k +1 is follo w ed immediately by the righ t reﬂectors P k − 2 ,k , P k − 1 ,k , so that only the last 2 × 2 principal submatrix of the transformed T k will b e changed in future iterations. These ideas can b e understo od more easily from Figure 5.1 and the follo wing compact form, whic h represents the actions of right reﬂectors on T k (additional to Q k,k +1 ( 3.7 )):    γ (5) k − 2  k ϑ k − 1 γ (4) k − 1 δ (2) k γ (2) k      c k 2 s k 2 1 s k 2 − c k 2     1 c k 3 s k 3 s k 3 − c k 3   =    γ (6) k − 2 ϑ (2) k − 1 γ (4) k − 1 δ (3) k η k γ (3) k      1 c k 3 s k 3 s k 3 − c k 3   =    γ (6) k − 2 ϑ (2) k − 1 γ (5) k − 1 η k ϑ k γ (4) k    . (5.2) 5.2. The diagonals of L k . Figure 5.2 sho ws the relation b et w een the singular v alues of A and the diagonal elemen ts of R k and L k with k = 19. This illustrates ( 4.2 ) for matrix ID 1177 from [ 54 ] with n = 25. 8 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS T 5 R 5 L 5 T 5 L 5 • P 3 , 5 • P 1 , 2 Q 4 , 5 • • P 2 , 3 • P 4 , 5 • P 1 , 3 Q 1 , 2 • Q 3 , 4 • Q 2 , 3 • • P 2 , 4 • P 3 , 4 Q 5 , 6 • • P 3 , 5 • P 2 , 3 Q 3 . 4 • • P 2 , 4 • P 4 , 5 Q 4 , 5 • Q 1 , 2 • • P 1 , 2 Q 2 , 3 • • P 3 , 4 Q 5 , 6 • • P 1 , 3 Fig. 5.1 . QLP with left and right reﬂe ctors interle aved on T 5 . This ﬁgur e c an b e r epr oduc e d with the help of QLPfig5.m . 0 10 20 0 2 4 6 8 nonz ero σ ( A ) = | λ ( A ) | , n = 25 0 10 20 0 2 4 6 8 γ Q k , k = 19 0 10 20 10 −8 10 −4 10 0 | γ Q k − σ ( A ) | 0 10 20 0 2 4 6 8 γ M k , k = 19 0 10 20 10 −4 10 −2 10 0 | γ M k − σ ( A ) | Fig. 5.2 . Upp er left: Nonzer o singular values of A sorted in decr e asing or der. Upp er midd le and right: The diagonals γ M k of R k (r e d circles) fr om MINRES ar e plotte d as r ed cir cles above or b e- low the ne ar est singular value of A . They appr oximate the extr eme nonzero singular values of A wel l. L ower: The diagonals γ Q k of L k (r e d cir cles) fr om MINRES-QLP appr oximate the extreme nonzer o singular values of A even b etter. A n implication is that the r atio of the lar gest and smal lest diagonals of L k pr ovides a go od estimate of κ ( A ) . T o r epro duc e this ﬁgur e, run test minresqlp fig3(2) . MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 9 5.3. Solving the subproblem. With y k = P k u k , subproblem ( 3.2 ) b ecomes u k = arg min u ∈ R k      L k 0  u −  t k φ k      , (5.3) where t k and φ k are as in ( 3.3 ) and ( 3.6 ). At the start of iteration k , the ﬁrst k − 3 elemen ts of u k , denoted by µ j for j ≤ k − 3, are known from previous iterations; see the 10th matrix in Figure 5.1 . The remainder depend on the rank of L k . 1. If rank( L k ) = k (so k < ` , or k = ` and b ∈ range( A )), w e need to solv e the last three equations of L k u k = t k :    γ (6) k − 2 ϑ (2) k − 1 γ (5) k − 1 η k ϑ k γ (4) k       µ (3) k − 2 µ (2) k − 1 µ k    =   ¯ τ k − 2 ¯ τ k − 1 τ k   ≡    τ k − 2 − η k − 2 µ (4) k − 4 − ϑ k − 2 µ (3) k − 3 τ k − 1 − η k − 1 µ (3) k − 3 τ k    . (5.4) 2. If k = ` and b 6∈ range( A ), the last row and column of L k are zero, and w e only need to solv e the last tw o equations of L k − 1 u k − 1 = t k − 1 , where L k =  L k − 1 0 0  , u k =  u k − 1 0  , " γ (6) k − 2 ϑ (2) k − 1 γ (5) k − 1 #" µ (3) k − 2 µ (2) k − 1 # =  ¯ τ k − 2 ¯ τ k − 1  . (5.5) The corresp onding solution estimate is x k = V k y k = V k P k u k = W k u k , where W k ≡ V k P k =  V k − 1 P k − 1 v k  P k − 2 ,k P k − 1 ,k (5.6) = h W (4) k − 3 w (3) k − 2 w (2) k − 1 v k i P k − 2 ,k P k − 1 ,k = h W (4) k − 3 w (4) k − 2 w (3) k − 1 w (2) k i , W T k W k = I k , range( W k ) = K k ( A, b ) , (5.7) and we up date x k − 2 and compute x k b y short-recurrence orthogonal steps: x (2) k − 2 = x (2) k − 3 + w (4) k − 2 µ (3) k − 2 , where x (2) k − 3 ≡ W (4) k − 3 u (3) k − 3 , (5.8) x k = x (2) k − 2 + w (3) k − 1 µ (2) k − 1 + w (2) k µ k . (5.9) 5.4. T ermination. When k = ` , Q k,k +1 is not formed or applied, see ( 3.3 ) and ( 3.7 ), and the QR factorization stops. In MINRES-QLP , we still need to apply P k − 2 ,k P k − 1 ,k on the right to obtain the minimum-length solution; see Figure 5.1 . Theorem 5.1 ( [ 10 , Theorem 3.1] ). In MINRES-QLP , x ` = x † . Pr o of . When b ∈ range( A ), the pro of is the same as that for Theorem 3.1 . When b / ∈ range( A ), for all u = [ u ` − 1 µ k ] T ∈ R ` that solves ( 5.3 ), MINRES-QLP returns the min-length LS solution u ` = [ u ` − 1 0] T b y the construction in ( 5.5 ). F or an y x ∈ range( W ` ) = K ` ( A, b ) by ( 5.7 ), k Ax − b k = k AW ` u − b k = k AV ` P ` u − b k = k V ` T ` P ` u − β 1 V ` e 1 k = k T ` P ` u − β 1 e 1 k =     Q ` − 1 T ` P ` u −  t ` − 1 φ ` − 1      =      L ` − 1 0 0 0  u −  t ` − 1 φ ` − 1      . Since L ` − 1 is nonsingular, φ ` − 1 = min k Ax − b k can b e ac hiev ed b y x ` = W ` u ` = W ` − 1 u ` − 1 and k x ` k = k W ` − 1 u ` − 1 k = k u ` − 1 k b y ( 5.7 ). Thus x ` is the min-length LS 10 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS solution of k Ax − b k in K ` ( A, b ), i.e., x ` = arg min {k x k | A 2 x = Ab, x ∈ K ` ( A, b ) } . Lik ewise y ` = P ` u ` is the min-length LS solution of k T ` y − β 1 e 1 k and so y ` ∈ range( T ` ), i.e. y ` = T ` z for some z . Thus x ` = V ` y ` = V ` T ` z = AV ` z ∈ range( A ). W e know that x † = arg min {k x k | A 2 x = Ab, x ∈ R n } is unique and x † ∈ range( A ). Since x ` ∈ range( A ), we must hav e x ` = x † . 5.5. T ransfer from MINRES to MINRES-QLP . On well-conditioned sys- tems, MINRES and MINRES-QLP b eha ve very similarly . How ev er, MINRES-QLP re- quires one more v ector of storage, and eac h iteration needs 4 more axp y’s ( y ← α x + y ) and 3 more vector scalings ( x ← αx ). Thus it would b e a desirable feature to inv ok e 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 T k is not yet to o ill-conditioned. The MINRES and MINRES-QLP solution estimates are the same, so from ( 3.4 ), ( 5.9 ), and ( 5.3 ): x M k = x k ⇐ ⇒ D k t k = W k u k = W k L − 1 k t k . Now from ( 3.8 ), ( 5.1 ), and ( 5.6 ), D k L k = ( V k R − 1 k )( R k P k ) = V k P k = W k , (5.10) and the last three columns of W k can b e obtained from the last three columns of D k and L k . (Th us, 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, w e need to generate x (2) k − 2 using ( 5.8 ): x (2) k − 2 = x M k − w (3) k − 1 µ (2) k − 1 − w (2) k µ k . It is clear from ( 5.10 ) that w e still need to do the righ t transformation R k P k = L k in the MINRES phase and k eep the last 3 × 3 principal submatrix of L k for each k so that we are ready to transfer to MINRES-QLP when necessary . W e then obtain a short recurrence for k x k k (see section 6.5 ) and for this computation we sav e ﬂops relativ e to the original MINRES algorithm, where k x k k is computed directly . In the implementation, the MINRES iterates transfer to MINRES-QLP iterates when an estimate of the condition n um b er of T k (see ( 6.3 )) exceeds an input param- eter tr anc ond . Thus, tr anc ond > 1 /ε leads to MINRES iterates throughout, while tr anc ond = 1 generates MINRES-QLP iterates from the start. 5.6. Comparison of Lanczos-based solv ers. W e compare MINRES-QLP with CG , SYMMLQ , and MINRES in T ables 5.1 – 5.2 in terms of subproblem deﬁnitions, basis, solution estimates, ﬂops and memory . A careful implementation of SYMMLQ pro vides a point in K k +1 ( A, b ) as sho wn. All solv ers need storage for v k , v k +1 , x k , and a pro duct p k = Av k eac h iteration. Some additional work-v ectors are needed for eac h metho d (e.g., d k − 1 and d k for MINRES , giving 7 work-v ectors in total). 6. Stopping conditions and norm estimates. This section deriv es sev eral norm estimates that are computed in MINRES-QLP . As before, we assume exact arithmetic throughout, so that V k and Q k are orthonormal. T able 6.1 summarizes ho w the norm estimates are used to formulate three groups of stopping conditions. The second NRBE test k Ar k k ≤ k A kk r k k tol is from Stewart [ 50 ] with symmetric A . 6.1. Residual and residual norm. First we derive recurrence relations for r k and its norm φ k ≡ k r k k . Lemma 6.1 ( r k and k r k k for MINRES-QLP and monotonicity of k r k k ). • If k < ` , then rank( L k ) = k , r k = s 2 k r k − 1 − φ k c k v k +1 , and φ k = φ k − 1 s k > 0 . • If rank( L ` ) = ` , then r ` = 0 . • If rank( L ` ) = ` − 1 , then r ` = r ` − 1 6 = 0 , and k r ` k = φ ` − 1 > 0 . MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 11 T able 5.1 Subpr oblems deﬁning x k for CG, SYMMLQ, MINRES, and MINRES-QLP. Metho d Subproblem F actorization Estimate of x k cgLanczos T k y k = β 1 e 1 Cholesky: x C k = V k y k [ 24 , 39 , 48 ] T k = L k D k L T k ∈ K k ( A, b ) SYMMLQ y k +1 = arg min y ∈ R k +1 k y k LQ: x L k = V k +1 y k +1 [ 39 , 45 ] s.t. T T k y = β 1 e 1 T k T Q T k =  L k 0  ∈ K k +1 ( A, b ) MINRES y k = arg min y ∈ R k k T k y − β 1 e 1 k QR: x M k = V k y k [ 39 ] Q k T k =  R k 0  ∈ K k ( A, b ) MINRES-QLP y k = arg min y ∈ R k k y k QLP: x Q k = V k y k [ 10 ] s.t. y ∈ arg min k T k y − β 1 e 1 k Q k T k P k =  L k 0  ∈ K k ( A, b ) T able 5.2 Bases, subpr oblem solutions, stor age, and work for e ach metho d. Metho d New basis z k , t k , u k x k estimate v ecs ﬂops cgLanczos W k ≡ V k L − T k L k D k z k = β 1 e 1 x C k = W k z k 5 8 n SYMMLQ W k ≡ V k +1 Q T k  I k 0  L k z k = β 1 e 1 x L k = W k z k 6 9 n MINRES D k ≡ V k R − 1 k t k = β 1  I k 0  Q k e 1 x M k = D k t k 7 9 n MINRES-QLP W k ≡ V k P k L k u k = β 1  I k 0  Q k e 1 x Q k = W k u k 8 14 n Pr o of . If k < ` , the residual is the same as for MINRES . W e hav e k r k k = φ k = φ k − 1 s k > 0; see ( 3.6 )–( 3.9 ). Also from r k = φ k V k +1 Q T k e k +1 ( 3.10 ) we hav e r k = φ k  V k v k +1   Q T k − 1 1    I k − 1 c k s k s k − c k     0 0 1   b y ( 3.3 ) , = φ k  V k v k +1   Q T k − 1 1  s k e k − c k  = φ k  V k v k +1   s k Q T k − 1 e k − c k  = φ k s k V k Q T k − 1 e k − φ k c k v k +1 = φ k − 1 s 2 k V k Q T k − 1 e k − φ k c k v k +1 = s 2 k r k − 1 − φ k c k v k +1 b y ( 3.10 ) . If T ` is nonsingular, r ` = 0. Otherwise Q ` − 1 ,` has made the last row of R ` zero, so the last row and column of L ` are zero; see ( 5.5 ). Th us r ` = r ` − 1 6 = 0; see Remark 3.2 . 6.2. Norm of Ar k . Next we derive recurrence relations for Ar k and its norm ψ k ≡ k Ar k k , and we sho w that Ar k is orthogonal to K k ( A, b ). Lemma 6.2 ( Ar k and ψ k ≡ k Ar k k for MINRES-QLP ). • If k < ` , then rank( L k ) = k , Ar k = k r k k ( γ k +1 v k +1 + δ k +2 v k +2 ) and ψ k = k r k kk [ γ k +1 δ k +2 ] k , wher e δ k +2 = 0 if k = ` − 1 . • If rank( L ` ) = ` , then Ar ` = 0 and ψ ` = 0 . • If rank( L ` ) = ` − 1 , then Ar ` = Ar ` − 1 = 0 , and k ψ ` k = ψ ` − 1 = 0 . Pr o of . F or the ﬁrst case, the pro of is essentially the same as the pro of of Lemma 3.3 . F or the other tw o cases, the results follow directly from Lemma 6.1 . 12 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS T able 6.1 Stopping c onditions in MINRES-QLP. NRBE me ans normwise r elative b ackwar d err or, and tol , maxit , maxc ond , and maxxnorm ar e input par ameters. Al l norms and κ ( A ) ar e estimate d by MINRES-QLP. Lanczos NRBE Regularization attempts β k +1 ≤ n k A k ε k r k k / ( k A kk x k k + k b k ) ≤ tol κ ( A ) ≥ maxc ond k = maxit k Ar k k / ( k A kk r k k ) ≤ tol k x k k ≥ maxxnorm 6.3. Matrix norms. F or Lanczos-based algorithms, k A k ≥ k V T k +1 AV k k = k T k k . Deﬁne A (0) ≡ 0 , A ( k ) ≡ max j =1 ,...,k n k T j e j k o = max n A ( k − 1) , k T k e k k o for k ≥ 1 . (6.1) Then k A k ≥ k T k k ≥ A ( k ) . Clearly , A ( k ) is monotonically increasing and is thus an impro ving estimate for k A k as k increases. By the prop ert y of QLP decomposition in ( 4.2 ) and ( 5.2 ), w e could easily extend ( 6.1 ) to include the largest diagonal of L k : A (0) ≡ 0 , A ( k ) ≡ max {A ( k − 1) , k T k e k k , γ (6) k − 2 , γ (5) k − 1 , | γ (4) k |} for k ≥ 1 . (6.2) Some other schemes inspired b y Larsen [ 31 , section A.6.1], Higham [ 25 ], and Chen and Demmel [ 8 ] follo w. F or the latter scheme, w e use an implemen tation b y Kaustuv [ 28 ] for estimating the norms of the rows of A . 1. [ 31 ] k T k k 1 ≥ k T k k 2. [ 31 ] q k T k T T k k 1 ≥ k T k k 3. [ 31 ] k T j k ≤ k T k k for small j = 5 or 20 4. [ 25 ] Ma tlab function NORMEST ( A ), which is based on the p o wer metho d 5. [ 8 ] max i k h i k / √ m , where h T i is the i th row of AZ , each column of Z ∈ R n × m is a random v ector of ± 1’s, and m is a small integer (e.g., m = 10). Figure 6.1 plots estimates of k A k for 12 matrices from the Florida sparse matrix collection [ 54 ] whose sizes n v ary from 25 to 3002. In particular, scheme 3 abov e with j = 20 giv es signiﬁcan tly more accurate estimates than other schemes for the 12 matri- ces w e tried. Ho wev er, the choice of j is not alwa ys clear and the scheme adds a little to the cost of MINRES-QLP . Hence we propose incorp orating it into MINRES-QLP (or other Lanzcos-based iterativ e metho ds) if very accurate k A k is needed. Otherwise ( 6.2 ) uses quantities readily a v ailable from MINRES-QLP and giv es us satisfactory estimates for the order of k A k . 6.4. Matrix condition n um b ers. W e again apply the prop erty of the QLP decomp osition in ( 4.2 ) and ( 5.2 ) to estimate κ ( T k ), which is a low er b ound for κ ( A ): γ min ← min { γ 1 , γ (2) 2 } , γ min ← min { γ min , γ (6) k − 2 , γ (5) k − 1 , | γ (4) k |} for k ≥ 3 , κ (0) ≡ 1 , κ ( k ) ≡ max  κ ( k − 1) , A ( k ) γ min  for k ≥ 1 . (6.3) 6.5. Solution norms. W e deriv e a recurrence relation for k x k k whose cost is as lo w as computing the norm of a 3- or 4- v ector. Since k x k k = k V k P k u k k = k u k k , w e can estimate k x k k b y computing χ k ≡ k u k k . Ho wev er, the last tw o elemen ts of u k c hange in u k +1 (and a new element µ k +1 is MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 13 1 2 3 4 5 6 7 8 9 10 11 12 10 −15 10 −10 10 −5 10 0 | m ax i {k T i e i k , | γ i |} − k A k | / k A k | k T k k 1 − k A k | / k A k ¯ ¯ ¯ ¯ q k T k T T k k 1 − k A k ¯ ¯ ¯ ¯ / k A k | σ 1 ( T 5 ) − k A k | / k A k | σ 1 ( T 2 0 ) − k A k | / k A k | N ORME S T ( A ) − k A k | / k A k Chen a nd Demme l 1 2 3 4 5 6 7 8 9 10 11 12 10 −2 10 −1 10 0 10 1 | k T k | k F − k A k F | / k A k F Fig. 6.1 . R elative err ors in diﬀer ent estimates of k A k . This ﬁgur e can b e repr o duc e d by testminresQLPNormA8 . added). W e therefore maintain χ k − 2 b y up dating it and then using it according to χ (2) k − 2 = k [ χ (2) k − 3 µ (3) k − 2 ] k , χ k = k [ χ (2) k − 2 µ (2) k − 1 µ k ] k cf. ( 5.8 ) and ( 5.9 ) . Th us χ (2) k − 2 increases monotonically but we cannot guarantee that k x k k and its recurred estimate χ k are increasing, and indeed they are not in som e examples (see Figure 8.1 ). 6.6. Pro jection norms. Sometimes the pro jection of the righ t-hand side vector b on to K k ( A, b ) is required (for example, see [ 46 ]). A simple recurrence relation is ω 2 k ≡ k Ax k k 2 = ω 2 k − 1 + τ 2 k and w e can deriv e it in the same wa y as sho wn in Lemma 3.3 . With ω 0 ≡ 0 we hav e ω k ≡ k Ax k k = k [ ω k − 1 τ k ] k . 7. Preconditioned MINRES and MINRES-QLP . It is often ask ed: How can w e construct a preconditioner for a linear system solver so that the same problem is solv ed with fewer iterations? Previous work on preconditioning the symmetric solvers CG , SYMMLQ , or MINRES includes [ 43 , 37 , 17 , 12 , 14 , 35 , 42 , 34 , 20 , 2 , 52 ]. W e hav e the same question for singular symmetric systems Ax ≈ b . Tw o-sided preconditioning is needed to preserv e symmetry . W e can still solv e compatible sys- tems, but w e will no longer obtain the minimum-length solution. F or incompatible systems, preconditioning alters the “least squares” norm. T o av oid this diﬃculty w e m ust w ork with larger equiv alent systems that are compatible. W e consider eac h case in turn, using a positive-deﬁnite preconditioner M = C C T with MINRES and MINRES-QLP to solve symmetric compatible systems Ax = b . Implicitly , w e are solv- ing equiv alent symmetric systems C − 1 AC − T y = C − 1 b , where C T x = y . As usual, it is p ossible to work with M itself, so without loss of generality we can assume C = M 1 2 . 7.1. Deriv ation. W e derive preconditioned MINRES for compatible Ax = b b y applying MINRES to the equiv alent problem ¯ A ¯ x = ¯ b , where ¯ A ≡ M − 1 2 AM − 1 2 , ¯ b ≡ M − 1 2 b , and x = M − 1 2 ¯ x . 7.1.1. Preconditioned Lanczos pro cess. Let v k denote the Lanczos v ectors of K ( ¯ A, ¯ b ). With v 0 = 0 and β 1 v 1 = ¯ b , for k = 1 , 2 , . . . w e deﬁne z k = β k M 1 2 v k , q k = β k M − 1 2 v k , so that M q k = z k . (7.1) 14 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS Then β k = k β k v k k = k M − 1 2 z k k = k z k k M − 1 = k q k k M = q q T k z k , where the square ro ot is well deﬁned b ecause M is p ositiv e deﬁnite, and the Lanczos iteration is p k = ¯ Av k = M − 1 2 AM − 1 2 v k = M − 1 2 Aq k /β k , α k = v T k p k = q T k Aq k /β 2 k , β k +1 v k +1 = M − 1 2 AM − 1 2 v k − α k v k − β k v k − 1 . Multiplying the last equation by M 1 2 w e get z k +1 = β k +1 M 1 2 v k +1 = AM − 1 2 v k − α k M 1 2 v k − β k M 1 2 v k − 1 = 1 β k Aq k − α k β k z k − β k β k − 1 z k − 1 . The last expression inv olving consecutiv e z j ’s replaces the three-term recurrence in v j ’s. In addition, w e need to solve a linear system M q k = z k ( 7.1 ) each iteration. 7.1.2. Preconditioned MINRES . F rom ( 3.4 ) and ( 3.8 ) we hav e the follo wing recurrence for the k th column of D k = V k R − 1 k and ¯ x k : d k =  v k − δ (2) k d k − 1 −  k d k − 2  /γ (2) k , ¯ x k = ¯ x k − 1 + τ k d k . Multiplying the ab o ve t wo equations by M − 1 2 on the left and deﬁning ¯ d k = M − 1 2 d k , w e can up date the solution of our original problem by ¯ d k =  1 β k q k − δ (2) k ¯ d k − 1 −  k ¯ d k − 2   γ (2) k , x k = M − 1 2 ¯ x k = x k − 1 + τ k ¯ d k . W e list the algorithm in [ 10 , T able 3.4]. 7.1.3. Preconditioned MINRES-QLP . A preconditioned MINRES-QLP can b e deriv ed very similarly . The additional work is to apply right reﬂectors P k to R k , and the new subproblem bases are W k ≡ V k P k , with ¯ x k = W k u k . Multiplying the new basis and solution estimate by M − 1 2 on the left, w e obtain W k ≡ M − 1 2 W k = M − 1 2 V k P k , x k = M − 1 2 ¯ x k = M − 1 2 W k u k = W k u k = x (2) k − 2 + µ (2) k − 1 ¯ w (3) k − 1 + µ k ¯ w (2) k . Algorithm 1 lists all steps. Note that ¯ w k is written as w k for all relev ant k . Also, the output x solves Ax ≈ b but the other outputs are asso ciated with ¯ A ¯ x ≈ ¯ b . R emark. The requiremen t of p ositiv e-deﬁnite preconditioners M in MINRES and MINRES-QLP may seem unnatural for a problem with indeﬁnite A b ecause w e cannot ac hieve M − 1 2 AM − 1 2 ≈ I . How ever, as sho wn in [ 17 ], we can achiev e M − 1 2 AM − 1 2 ≈  I − I  using an appro ximate blo c k-LDL T factorization A ≈ LD L T to get M = L | D | L T , where D is indeﬁnite with blo c ks of order 1 and 2, and | D | has the same eigensystem as D except negative eigen v alues are changed in sign. SQMR [ 15 ] without preconditioning is analytically equiv alent to MINRES . Unlik e MINRES , SQMR can w ork directly with an indeﬁnite preconditioner (suc h as blo c k- LDL T ). How ever, in ﬁnite precision, SQMR needs “lo ok-ahead” to prev ent numerical breakdo wn. MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 15 Algorithm 1: Preconditioned MINRES-QLP to solve ( A − σ I ) x ≈ b . input : A, b, σ, M 1 z 0 = 0, z 1 = b , Solv e M q 1 = z 1 , β 1 = p b T q 1 [Initialize] 2 w 0 = w − 1 = 0, x − 2 = x − 1 = x 0 = 0 3 c 0 , 1 = c 0 , 2 = c 0 , 3 = − 1, s 0 , 1 = s 0 , 2 = s 0 , 3 = 0, φ 0 = β 1 , τ 0 = ω 0 = χ − 2 = χ − 1 = χ 0 = 0 4 δ 1 = γ − 1 = γ 0 = η − 1 = η 0 = η 1 = ϑ − 1 = ϑ 0 = ϑ 1 = µ − 1 = µ 0 = 0, A = 0 , κ = 1 5 k = 0 6 while no stopping condition is satisﬁed do 7 k ← k + 1 8 p k = Aq k − σ q k , α k = 1 β 2 k q T k p k [Preconditioned Lanczos] 9 z k +1 = 1 β k p k − α k β k z k − β k β k − 1 z k − 1 10 Solv e M q k +1 = z k +1 , β k +1 = q q T k +1 z k +1 11 if k = 1 then ρ k = k [ α k β k +1 ] k else ρ k = k [ β k α k β k +1 ] k 12 δ (2) k = c k − 1 , 1 δ k + s k − 1 , 1 α k [Previous left reflection...] 13 γ k = s k − 1 , 1 δ k − c k − 1 , 1 α k [on middle two entries of T k e k ...] 14  k +1 = s k − 1 , 1 β k +1 [produces first two entries in T k +1 e k +1 ] 15 δ k +1 = − c k − 1 , 1 β k +1 16 c k 1 , s k 1 , γ (2) k ← SymOrtho( γ k , β k +1 ) [Current left reflection] 17 c k 2 , s k 2 , γ (6) k − 2 ← SymOrtho( γ (5) k − 2 ,  k ) [First right reflection] 18 δ (3) k = s k 2 ϑ k − 1 − c k 2 δ (2) k , γ (3) k = − c k 2 γ (2) k , η k = s k 2 γ (2) k 19 ϑ (2) k − 1 = c k 2 ϑ k − 1 + s k 2 δ (2) k 20 c k 3 , s k 3 , γ (5) k − 1 ← SymOrtho( γ (4) k − 1 , δ (3) k ) [Second right reflection...] 21 ϑ k = s k 3 γ (3) k , γ (4) k = − c k 3 γ (3) k [to zero out δ (3) k ] 22 τ k = c k 1 φ k − 1 [Last element of t k ] 23 φ k = s k 1 φ k − 1 , ψ k − 1 = φ k − 1 k [ γ k δ k +1 ] k [Update k r k k , k Ar k − 1 k ] 24 if k = 1 then γ min = γ 1 else γ min ← min { γ min , γ (6) k − 2 , γ (5) k − 1 , | γ (4) k |} 25 A ( k ) = max {A ( k − 1) , ρ k , γ (6) k − 2 , γ (5) k − 1 , | γ (4) k |} [Update k A k ] 26 ω k = k [ ω k − 1 τ k ] k , κ ← A ( k ) /γ min [Update k Ax k k , κ ( A ) ] 27 w k = − ( c k 2 /β k ) q k + s k 2 w (3) k − 2 [Update w k − 2 , w k − 1 , w k ] 28 w (4) k − 2 = ( s k 2 /β k ) q k + c k 2 w (3) k − 2 29 if k > 2 then w (2) k = s k 3 w (2) k − 1 − c k 3 w k , w (3) k − 1 = c k 3 w (2) k − 1 + s k 3 w k 30 if k > 2 then µ (3) k − 2 = ( τ k − 2 − η k − 2 µ (4) k − 4 − ϑ k − 2 µ (3) k − 3 ) /γ (6) k − 2 [Update µ k − 2 ] 31 if k > 1 then µ (2) k − 1 = ( τ k − 1 − η k − 1 µ (3) k − 3 − ϑ (2) k − 1 µ (3) k − 2 ) /γ (5) k − 1 [Update µ k − 1 ] 32 if γ (4) k 6 = 0 then µ k = ( τ k − η k µ (3) k − 2 − ϑ k µ (2) k − 1 ) /γ (4) k else µ k = 0 [Compute µ k ] 33 x (2) k − 2 = x (2) k − 3 + µ (3) k − 2 w (3) k − 2 [Update x k − 2 ] 34 x k = x (2) k − 2 + µ (2) k − 1 w (3) k − 1 + µ k w (2) k [Compute x k ] 35 χ (2) k − 2 = k [ χ (2) k − 3 µ (3) k − 2 ] k [Update k x k − 2 k ] 36 χ k = k [ χ (2) k − 2 µ (2) k − 1 µ k ] k [Compute k x k k ] 37 x = x k , φ = φ k , ψ = φ k k [ γ k +1 δ k +2 ] k , χ = χ k , A = A ( k ) , ω = ω k output : x, φ, ψ , χ, A , κ, ω 38 [ c, s ← SymOrtho ( a, b ) is a stable form for computing r = √ a 2 + b 2 , c = a r , s = b r ] 16 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS 7.2. Preconditioning singular Ax = b . F or singular compatible systems, MIN- RES and MINRES-QLP ﬁnd the minimum-length solution (see Theorems 3.1 and 5.1 ). If M is nonsingular, the preconditioned system is also compatible and the solvers return its minim um-length solution. The unpreconditioned solution solv es Ax ≈ b , but is not necessarily a minimum-length solution. Example 7.1. L et A =  1 1 0 0 1 1 1 0 0 1 0 1 0 0 1 0  and b =  6 9 6 3  . Then rank( A ) = 3 and Ax = b is a singular c omp atible system. The minimum-length solution is x † = [ 2 4 3 2 ] T . By binormalization [ 33 ] we c onstruct the matrix D = diag ([ 0 . 84201 0 . 81228 0 . 30957 3 . 2303 ]) . The minimum-length solution of the diagonal ly pr e c onditione d pr oblem DAD y = D b is y † = [ 3 . 5739 3 . 6819 9 . 6909 0 . 93156 ] T . Then x = D y † = [ 3 . 0092 2 . 9908 3 . 0000 3 . 0092 ] T is a solution of Ax = b , but x 6 = x † . 7.3. Preconditioning singular Ax ≈ b . W e prop ose the following tec hniques for obtaining minimum-residual solutions of singular incompatible problems. In each case we use an equiv alen t but larger c omp atible system to whic h MINRES ma y be applied. Even if the larger system is singular, Theorem 3.1 shows that the minimum- length solution of the larger system will be obtained. The required x will b e part of this solution. Preconditioning still gives a minimum-residual solution of Ax ≈ b , and in some cases x will b e x † . If the systems are ill-conditioned, it will be safer and more eﬃcien t to apply MINRES-QLP to the original incompatible system. Ho wev er, preconditioning will give an x that is “minimum length” in a diﬀerent norm. 7.3.1. Augmen ted system. When A is singular, so is the augmented system  I A A  r x  =  b 0  , (7.2) but it is alwa ys compatible. Preconditioning with symmetric p ositiv e-deﬁnite M giv es us a solution [ r x ] in which r is unique, but x may not b e x † . 7.3.2. A giant KKT system. Problem ( 1.1 ) is equiv alen t to min r, x x T x sub ject to ( 7.2 ), which is an equality-con trained conv ex quadratic program. The corresp ond- ing KKT system [ 36 , section 16.1] is b oth symmetric and compatible:     I A − I A I A A         r x y z     =     0 0 b 0     . (7.3) Although this is still a singular system, the upper-left 3 × 3 blo c k-submatrix is nonsin- gular 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. 7.3.3. Regularization. If the rank of a given matrix A is ill-determined, w e ma y wan t to solve the r e gularize d problem [ 13 , 22 ] with parameter δ > 0: min x      A δ I  x −  b 0      2 . (7.4) The matrix  A δ I  has full rank and is alw ays b etter conditioned than A . LSQR [ 40 , 41 ] ma y be applied, and its iterates x k will reduce k r k k 2 + δ 2 k x k k 2 monotonically . Alter- nativ ely , w e could transform ( 7.4 ) in to the follo wing symmetric compatible systems and apply MINRES or MINRES-QLP . They tend to reduce k Ar k − δ 2 x k k monotonically . MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 17 Normal equation: ( A 2 + δ 2 I ) x = Ab. (7.5) Augmen ted system:  I A A − δ 2 I  r x  =  b 0  . A t w o-la y ered problem: If we eliminate v from the system  I A 2 A 2 − δ 2 A 2  x v  =  0 Ab  . (7.6) w e obtain ( 7.5 ). Thus x is also a solution of our regularized problem ( 7.4 ). This is equiv alent to the tw o-lay ered formulation (4.3) in Bobro vniko v a and V av asis [ 5 ] (with A 1 = A , A 2 = D 1 = D 2 = I , b 1 = b , b 2 = 0, δ 1 = 1, δ 2 = δ 2 ). A key prop ert y is that x → x † as δ → 0. A KKT-lik e system: If w e deﬁne y = − Av and r = b − Ax − δ 2 y , then w e can sho w (b y eliminating r and y from the following system) that x in     I A − I A I A δ 2 I A         r x y v     =     0 0 b 0     (7.7) is also a solution of ( 7.6 ) and thus of ( 7.4 ). The upper-left 3 × 3 blo c k- submatrix of ( 7.7 ) is nonsingular, and the correct limiting b eha vior o ccurs: x → x † as δ → 0. In fact, ( 7.7 ) reduces to ( 7.3 ). 7.4. General preconditioners. The construction of preconditioners is usually problem-dep enden t. If not muc h is kno wn about the structure of A , w e can only consider general metho ds such as diagonal preconditioning and incomplete Cholesky factorization. These metho ds require access to the nonzero elements of A . (They are not applicable if A exists only as an op erator for returning the pro duct Ax .) F or a comprehensiv e survey of preconditioning techniques, see Benzi [ 3 ]. W e discuss a few metho ds for symmetric A that also require access to the nonzero A ij . 7.4.1. Diagonal preconditioning. If A has entries that are v ery diﬀeren t in magnitude, diagonal scaling migh t impro ve its condition. When A is diagonally dom- inan t and nonsingular, w e can deﬁne D = diag( d 1 , . . . , d n ) with d j = 1 / | A j j | 1 / 2 . Instead of solving Ax = b , we solve D ADy = Db , where D AD is still diagonally dominan t and nonsingular with all entries ≤ 1 in magnitude, and x = D y . More generally , if A is not diagonally dominant and p ossibly singular, w e can safeguard division-by-zero errors b y c ho osing a parameter δ > 0 and deﬁning d j ( δ ) = 1 / max { δ , √ | A j j | , max i 6 = j | A ij |} , j = 1 , . . . , n. (7.8) Example 7.2. 1. If A = " − 1 10 − 8 10 − 8 1 10 4 10 4 0 0 # , then κ ( A ) ≈ 10 4 . L et δ = 1 , D =  1 10 − 2 10 − 2 1  in ( 7.8 ) . Then D AD = " − 1 10 − 10 10 − 10 10 − 4 1 1 0 0 # and κ ( D AD ) ≈ 1 . 18 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS 2. A = " 10 − 4 10 − 8 10 − 8 10 − 4 10 − 8 10 − 8 0 0 # c ontains mostly very smal l entries, and κ ( A ) ≈ 10 10 . L et δ = 10 − 8 and D = " 10 2 10 2 10 8 10 8 # . Then DAD = " 1 10 − 4 10 − 4 1 10 2 10 2 0 0 # and κ ( D AD ) ≈ 10 2 . (The choic e of δ makes a critic al diﬀer enc e in this c ase: with δ = 1 , we have D = I .) 7.4.2. Binormalization (BIN). Livne and Golub [ 33 ] scale a symmetric ma- trix b y a series of k diagonal matrices on b oth sides until all rows and columns of the scaled matrix hav e unit 2-norm: D AD = D k · · · D 1 AD 1 · · · D k . See also Bradley [ 6 ]. Example 7.3. If A =  10 − 8 1 1 10 − 8 10 4 10 4 0  , then κ ( A ) ≈ 10 12 . With just one swe ep of BIN, we obtain D = diag (8 . 1 e − 3 , 6 . 6 e − 5 , 1 . 5) , D AD ≈ h 6 . 5 e − 1 5 . 3 e − 1 0 5 . 3 e − 1 0 1 0 1 0 i and κ ( DAD ) ≈ 2 . 6 even though the r ows and c olumns have not c onver ge d to one in the two-norm. In c ontr ast, diagonal sc aling ( 7.8 ) deﬁne d by δ = 1 and D = diag(1 , 10 − 4 , 10 − 4 ) r e duc es the c ondition numb er to appr oximately 10 4 . 7.4.3. Incomplete Cholesky factorization. F or a sparse symmetric p ositiv e deﬁnite matrix A , w e could compute a preconditioner b y the incomplete Cholesky factorization that preserves the sparsity pattern of A . This is known as IC0 in the literature. Often there exists a p erm utation P suc h that the IC0 factor of P AP T is more sparse than that of A . When A is semideﬁnite or indeﬁnite, IC0 may not exist, but a simple v arian t that ma y work is the incomplete Cholesky-inﬁnity factorization [ 55 , section 5]. 8. Numerical exp erimen ts. W e compare the computed results of MINRES- QLP and v arious other Krylov subspace methods to solutions computed directly b y the eigen v alue decomposition (EVD) and the truncated eigen v alue decompositions (TEVD) of A . F or TEVD we ha ve x t ≡ X | λ i | >t k A k ε 1 λ i u i u T i b, k A k = max | λ i | , κ t ( A ) = max | λ i | min | λ i | >t k A k ε | λ i | with parameter t > 0. Often t is set to 1, and sometimes to a moderate n umber suc h as 10 or 100; it deﬁnes a cut-oﬀ point relative to the largest eigen v alue of A . F or example, if most eigen v alues are of order 1 in magnitude and the rest are of order k A k ε ≈ 10 − 16 , we exp ect TEVD to work b etter when the small eigenv alues are excluded, while EVD (with t = 0) could return an explo ding solution. In the tables of results, Ma tlab MINRES and Ma tlab SYMMLQ are Ma tlab ’s implemen tation of MINRES and SYMMLQ respectively . They incorp orate lo c al r e- ortho gonalization of the Lanczos vector v 2 , which could enhance the accuracy of the computations if b is close to an eigenv ector of A [ 32 ]: Second Lanczos iteration: β 1 v 1 = b, and q 2 ≡ β 2 v 2 = Av 1 − α 1 v 1 Lo cal reorthogonalization: q 2 ← q 2 − ( v T 1 q 2 ) v 1 . Lac king the correct stopping condition for singular problems, Ma tlab SYMMLQ w orks more than necessary and then selects the smallest residual norm from all com- puted iterates; it would sometimes report that the metho d did not con v erge although the selected estimate app eared to b e reasonably accurate. MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 19 MINRES SOL and SYMMLQ SOL are implemen tations based on [ 39 ]. MINRES + and SYMMLQ + are the same but with additional stopping conditions for singular incompatible systems (see Lemma 3.3 and [ 10 , Prop osition 2.12]). The computations in this section w ere p erformed on a Windo ws XP mac hine with a 3.2GHz In tel P en tium D Processor 940 and 3GB RAM ( ε ≈ 10 − 16 ) . T ests were p erformed with each solver on ﬁve types of problem: 1. symmetric, nonsingular linear systems 2. symmetric, singular linear systems 3. mildly incompatible symmetric (and singular) systems (meaning k r k is rather small with resp ect to k b k ) 4. symmetric (and singular) LS problems 5. Hermitian systems. W e present a few examples that illustrate the k ey features of MINRES-QLP . F or a larger set of tests and results, such as applying MINRES-QLP and other Krylov metho ds to Hermitian systems with preconditioning, we refer to [ 10 , Chapter 4]. F or a compatible system, we generate a random vector b that is in the range of the test matrix ( b ≡ Ay , y i ∼ i.i.d. U (0 , 1), i.e., y 1 , . . . , y n are independent and iden tically distributed random v ariables, whose v alues are dra wn from the standard uniform distribution with supp ort [0 , 1]). F or an LS problem, we generate a random b that is not in the range of the test matrix ( b i ∼ i.i.d. U (0 , 1) often suﬃces). If A is Hermitian, then v H Av is real for all complex vectors v . Numerically (in double precision), α k = v H k Av k is lik ely to hav e small imaginary parts in the ﬁrst few Lanczos iterations and sno wball to ha ve large imaginary parts in later iterations. This w ould result in a p oor estimation of k T k k F or k A k F , and unnecessary errors in the Lanczos vectors. Thus we made sure to typ e c ast α k = real( v H k Av k ) in MINRES-QLP and MINRES SOL . W e could sa y from the results that the Lanczos-based metho ds hav e built-in regularization features [ 29 ], often matching the TEVD solutions very well. 8.1. A Laplacian system Ax ≈ b (almost compatible). Our ﬁrst example in volv es a singular indeﬁnite Laplacian matrix A of order n = 400. It is block- tridiagonal with each blo c k b eing a tridiagonal matrix T of order N = 20 with all nonzeros equal to 1: A =       T T T T . . . . . . . . . T T T       n × n , T =       1 1 1 1 . . . . . . . . . 1 1 1       N × N . (8.1) Ma tlab ’s eig( A ) rep orts the following data: 205 positive eigen v alues in the inter- v al [6 . 1e − 2 , 8 . 87], 39 almost-zero eigen v alues in [ − 2 . 18e − 15 , 3 . 71e − 15], 156 negative eigen v alues in [ − 2 . 91 , − 6 . 65e − 2], numerical rank = 361. W e used a right-hand side with a small incompatible comp onen t: b = Ay + 10 − 8 z with y i and z i ∼ i.i.d. U (0 , 1). Results are summarized in T able 8.1 . In the column lab eled “C?”, the v alue “Y” denotes that the asso ciated algorithm in the row has con verged to the desired NRBE tolerances within maxit iterations (cf. T able 6.1 ); otherwise, we hav e v alues “N” and “N?”, where “N?” indicates that the algorithm could hav e conv erged if more relaxed stopping conditions w ere used. The column “ Av ” sho ws the total num b er of matrix-vector pro ducts, and column “ x (1)” lists the 20 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS T able 8.1 Finite element pr oblem Ax ≈ b with b almost c omp atible. L aplacian on a 20 × 20 grid, n = 400 , maxit = 1200 , shift = 0 , tol = 1 . 0 e − 15 , maxnorm = 100 , maxc ond = 1 e 15 , k b k = 87 . T o r epr o duce this example, run test minresqlp eg7 1(24) . Method C? Av x (1) k x k k e k k r k k Ar k k A k κ ( A ) EVD – – − 7 . 39e5 4 . 12e7 4 . 1e7 1 . 7e − 7 7 . 8e − 7 8 . 9e0 1 . 1e17 TEVD – – 3 . 89e − 1 1 . 15e1 0 . 0e0 1 . 7e − 8 1 . 4e − 12 8 . 9e0 1 . 5e2 Ma tlab SYMMLQ N? 371 3 . 89e − 1 1 . 15e1 1 . 4e − 7 1 . 8e − 7 5 . 8e − 7 – – SYMMLQ SOL N 447 − 3 . 08e0 9 . 63e1 9 . 5e1 1 . 4e2 4 . 4e2 9 . 6e1 1 . 3e1 SYMMLQ + N 447 2 . 94e6 4 . 27e8 4 . 3e8 1 . 8e2 6 . 5e2 8 . 6e0 1 . 3e1 Ma tlab MINRES N 1200 − 7 . 50e5 2 . 10e7 2 . 1e7 1 . 5e7 9 . 1e7 – – MINRES SOL N 1200 9 . 89e5 6 . 10e7 6 . 1e7 2 . 3e7 1 . 5e8 1 . 8e2 1 . 5e1 MINRES + N 611 1 . 02e0 9 . 28e1 9 . 2e1 1 . 7e − 8 2 . 5e − 11 8 . 6e0 6 . 9e13 MINRES-QLP Y 612 3 . 89e − 1 1 . 15e1 3 . 7e − 11 1 . 7e − 8 9 . 3e − 11 8 . 7e0 4 . 3e13 Ma tlab LSQR Y 1462 3 . 89e − 1 1 . 15e1 2 . 3e − 13 1 . 7e − 8 3 . 3e − 13 – – LSQR SOL Y 1464 3 . 89e − 1 1 . 15e1 2 . 4e − 13 1 . 7e − 8 3 . 9e − 13 1 . 5e2 6 . 4e3 Ma tlab GMRES(30) N? 1200 3 . 90e − 1 1 . 15e1 5 . 2e − 2 3 . 4e − 3 9 . 4e − 4 – – SQMR N 1200 − 2 . 58e8 3 . 74e10 3 . 7e10 4 . 6e3 2 . 3e4 – – Ma tlab QMR N? 798 3 . 89e − 1 1 . 15e1 5 . 2e − 7 1 . 9e − 8 2 . 6e − 8 – – Ma tlab BICG N? 790 3 . 89e − 1 1 . 15e1 4 . 7e − 7 3 . 9e − 8 1 . 9e − 7 – – Ma tlab BICGST AB N? 2035 3 . 89e − 1 1 . 15e1 4 . 2e − 7 1 . 7e − 8 4 . 3e − 13 – – ﬁrst element of the ﬁnal solution estimate x for each algorithm. F or GMRES , the in teger in parentheses is the v alue of the restart parameter. MINRES SOL gives a larger solution than MINRES-QLP . This example has a residual norm of ab out 1 . 7 × 10 − 8 , so it is not clear whether to classify it as a linear system or an LS problem. T o the credit of Ma tlab SYMMLQ , it thinks the system is linear and returns a go od solution. F or MINRES-QLP , the ﬁrst 410 iterations are in standard “ MINRES mo de”, with a transfer to “ MINRES-QLP mo de” for the last 202 iterations. LSQR [ 40 , 41 ] con verges to the minim um-length solution but with more than t wice the n um b er of iterations of MINRES-QLP . The other solvers fall short in some wa y . 8.2. A Laplacian LS problem min k Ax − b k . This example uses the same Laplacian matrix A ( 8.1 ) but with a clearly incompatible b = 10 × rand( n, 1), i.e., b i ∼ i.i.d. U (0 , 10). The residual norm is about 17. Results are summarized in T able 8.2 . MINRES gives an LS solution, while MINRES-QLP is the only solv er that matc hes the solution of TEVD. The other solv ers do not p erform satisfactorily . 8.3. Regularizing eﬀect of MINRES-QLP . This example illustrates the reg- ularizing eﬀect of MINRES-QLP with the stopping condition χ k ≤ maxxnorm . F or k ≥ 18 in Figure 8.1 , we observe the following v alues: χ 18 = k  2 . 51 3 . 87e − 11 1 . 38 × 10 2  k = 1 . 38 × 10 2 , χ 19 = k  2 . 51 − 8 . 00e − 10 − 1 . 52 × 10 2  k = 1 . 52 × 10 2 , χ 20 = k  2 . 51 1 . 62e − 10 − 1 . 62 × 10 6  k = 1 . 62 × 10 6 > maxxnorm ≡ 10 4 . Because the last v alue exceeds maxxnorm , MINRES-QLP regards the last diagonal elemen t of L k as a singular v alue to b e ignored (in the spirit of truncated SVD solutions). It discards the last element of u 20 and up dates χ 20 ← k  2 . 51 1 . 62e − 10 0  k = 2 . 51 . MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 21 T able 8.2 Finite element pr oblem min k Ax − b k . L aplacian on a 20 × 20 grid, n = 400 , maxit = 500 , shift = 0 , tol = 1 . 0 e − 14 , maxnorm = 1 e 4 , maxc ond = 1 e 14 , k b k = 120 . T o r epr o duc e this example, run test minresqlp eg7 1(25) . Method C? Av x (1) k x k k e k k r k k Ar k k A k κ ( A ) EVD – – − 7 . 39e14 4 . 12e16 4 . 1e16 1 . 8e2 7 . 9e2 8 . 9e0 1 . 1e17 TEVD – – − 8 . 75e0 1 . 43e2 0 . 0e0 1 . 7e1 4 . 1e − 12 8 . 9e0 1 . 5e2 Ma tlab SYMMLQ N 1 2 . 74e − 1 1 . 52e1 1 . 4e2 6 . 0e1 2 . 9e2 – – SYMMLQ SOL N 228 − 7 . 70e2 9 . 93e3 9 . 9e3 7 . 0e3 3 . 4e4 6 . 8e1 9 . 7e0 SYMMLQ + N 228 − 7 . 70e2 9 . 93e3 9 . 9e3 7 . 0e3 3 . 4e4 7 . 6e0 9 . 7e0 Ma tlab MINRES N 500 2 . 80e14 4 . 07e16 4 . 1e16 2 . 3e2 1 . 4e3 – – MINRES SOL N 500 − 1 . 46e14 2 . 11e16 2 . 1e16 1 . 1e2 6 . 6e2 1 . 5e2 1 . 4e1 MINRES + N 381 3 . 88e1 6 . 90e3 6 . 9e3 1 . 7e1 1 . 2e − 5 7 . 9e0 1 . 6e10 MINRES-QLP Y 382 − 8 . 75e0 1 . 43e2 1 . 7e − 6 1 . 7e1 1 . 7e − 5 8 . 6e0 3 . 5e10 Ma tlab LSQR Y 1000 − 8 . 75e0 1 . 43e2 2 . 0e − 5 1 . 7e1 1 . 4e − 5 – – LSQR SOL Y 1000 − 8 . 75e0 1 . 43e2 2 . 3e − 5 1 . 7e1 1 . 1e − 5 1 . 2e2 4 . 4e3 Ma tlab GMRES(30) N 500 − 8 . 84e0 1 . 25e2 4 . 8e1 1 . 7e1 8 . 2e − 1 – – SQMR N 500 9 . 58e15 1 . 39e18 1 . 4e18 1 . 2e11 6 . 7e11 – – Ma tlab QMR N 556 − 7 . 30e0 1 . 92e2 1 . 4e2 1 . 7e1 1 . 2e1 – – Ma tlab BICG N 2 1 . 40e0 1 . 71e1 1 . 4e2 6 . 0e1 2 . 6e2 – – Ma tlab BICGST AB N 104 − 1 . 12e1 1 . 40e2 9 . 6e1 2 . 6e1 1 . 8e1 – – The full truncation strategy used in the implemen tation is justiﬁed b y the fact that x k = W k u k with W k orthogonal. When k x k k b ecomes large, the last element of u k is treated as zero. If k x k k is still large, the second-to-last element of u k is treated as zero. If k x k k is stil l large, the third-to-last element of u k is treated as zero. 8.4. Eﬀects of rounding errors in MINRES-QLP . The recurred residual norms φ M k in MINRES usually approximate the directly computed ones k r M k k very w ell un til k r M k k becomes small. W e observ e that φ M k con tinues to decrease in the last few iterations, ev en though k r M k k has become stagnant. This is desirable in the sense that the stopping rule will cause termination, although the ﬁnal solution is not as accurate as predicted. W e presen t similar plots of MINRES-QLP in the following examples, with the cor- resp onding quan tities as φ Q k and k r Q k k . W e observ e that except in v ery ill-conditioned LS problems, φ Q k appro ximates k r Q k k very closely . Figure 8.2 illustrates four singular compatible linear systems. Figure 8.3 illustrates four singular LS problems. 9. Conclusion. MINRES constructs its k th solution estimate from the recursion x k = D k t k = x k − 1 + τ k d k ( 3.4 ), where n separate triangular systems R T k D T k = V T k are solved to obtain the n elements of eac h direction d 1 , . . . , d k . (Only d k is obtained during iteration k , but it has n elements.) In con trast, MINRES-QLP constructs its k th solution estimate using orthogonal steps: x Q k = ( V k P k ) u k ; see ( 5.3 )–( 5.9 ). Only one triangular system L k u k = Q k ( β 1 e 1 ) is inv olv ed for each k . Th us MINRES-QLP o v ercomes the p oten tial instability predicted by the MINRES authors [ 39 ] and analyzed b y Sleijp en et al. [ 47 ]. The additional work and storage are mo derate, and maximum eﬃciency is retained by transferring from MINRES to the MINRES-QLP iterates only when the estimated condition num b er of A exceeds a sp eciﬁed v alue. 22 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS 5 10 15 20 10 −14 10 −7 10 0 k φ M k φ Q k N RB E φ M k N RB E φ Q k 5 10 15 20 10 −14 10 −7 10 0 k ψ M k ψ Q k NR B E ψ M k NR B E ψ Q k 5 10 15 20 10 1 10 2 k k x M k k χ Q k Fig. 8.1 . R e curr e d φ k ≈ k r k k , ψ k ≈ k Ar k k , and k x k k for MINRES and MINRES-QLP. The matrix A (ID 1177 fr om [ 54 ]) is p ositive semideﬁnite, n = 25 , and b is random with k b k ' 1 . 7 . Both solvers c ould have achieve d essential ly the TEVD solution of Ax ' b at iter ation 11 . However, the stringent tol = 10 − 14 on the r e curr ed normwise r elative b ackwar d err ors (NRBE in T able 6.1 ) pr events them fr om stopping “in time”. MINRES ends with an explo ding solution, yet MINRES- QLP brings it b ack to the TEVD solution at iteration 20 . L eft: φ M k and φ Q k (r e curr e d k r k k of MINRES and MINRES-QLP) and their NRBE. Midd le: ψ M k and ψ Q k (r e curr e d k Ar k k ) and their NRBE. R ight: k x M k k (norms of solution estimates fr om MINRES) and χ Q k (r e curr e d k x k k fr om MINRES-QLP) with maxxnorm = 10 4 . This ﬁgur e can be r epr o duce d by test minresqlp fig7 1(2) . MINRES and MINRES-QLP are readily applicable to Hermitian matrices, once α k is typecast as a real scalar in ﬁnite-precision arithmetic. F or both algorithms, w e deriv ed recurrence relations for k Ar k k and k Ax k 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 ﬁnd appro ximate solutions to min k Ax − b k that serve as a form of r e gularization . Krylo v subspace metho ds also ha ve regularization prop erties [ 23 , 21 , 29 ]. Since MINRES- QLP monitors more carefully the rank of T k , whic h could b e k or k − 1, w e may say that regularization is a stronger feature in MINRES-QLP , as w e ha v e sho wn in our n umerical examples. It is imp ortan t to develop robust techniques for estimating an a priori b ound for the solution norm since the MINRES-QLP appro ximations are not monotonic as is the case in CG and LSQR . Ideally , w e would also lik e to determine a practical threshold asso ciated with the stopping condition γ (4) k = 0 in order to handle cases when γ (4) k is n umerically small but not exactly zero. These are topics for future research. 10. Soft w are and reproducible research. Ma tlab 7.6, F ortran 77, and F or- tran 90 implementations of MINRES with new stopping conditions k Ar k k ≤ tol k A kk r k k and k Ax k k ≤ tol k A kk x k k , and a Ma tlab 7.6 implementation of MINRES-QLP are a v ailable from SOL [ 48 ]. F ollowing the philosophy of reproducible computational research as adv o cated in [ 11 , 9 ], for each ﬁgure and example in this pap er we mention either the source or the sp eciﬁc Ma tlab command. Our Ma tlab scripts are av ailable at SOL [ 48 ]. MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 23 5 10 15 20 25 10 −15 10 −10 10 −5 10 0 k κ ≈ 10 8 ,  A  =3 ,  x ≈ 10 8 ,  b ≈ 10 1 ,t o l =1 0 − 14  r M k  φ M k  r Q k  φ Q k 5 10 15 20 25 30 10 −15 10 −10 10 −5 10 0 κ ≈ 10 8 ,  A  =3 ,  x ≈ 10 1 ,  b ≈ 10 2 ,t o l =1 0 − 14 k 5 10 15 20 25 10 −15 10 −10 10 −5 10 0 κ ≈ 10 10 ,  A  =3 ,  x ≈ 10 10 ,  b ≈ 10 1 ,t o l =1 0 − 14 k 5 10 15 20 25 30 35 10 −15 10 −10 10 −5 10 0 κ ≈ 10 10 ,  A  =3 ,  x ≈ 10 1 ,  b ≈ 10 2 ,t o l =1 0 − 14 k Four symmetric positive semidefinite ill−conditioned systems. Fig. 8.2 . Solving Ax = b with semideﬁnite A similar to an example of Sleijp en et al. [ 47 ]. A = Q diag([0 5 , η , 2 η , 2 : 1 789 : 3]) Q of dimension n = 797 , nul lity 5, and norm k A k = 3 , wher e Q = I − (2 /n ) ww T is a Householder matrix gener ate d by v = [0 5 , 1 , . . . , 1] T , w = v / k v k . These plots il lustr ate and c ompar e the eﬀect of r ounding err ors in MINRES and MINRES-QLP. The upp er part of e ach plot shows the c omputed and r e curre d residual norms, and the lower p art shows the com pute d and r e curre d normwise r elative b ackwar d err ors (NRBE, deﬁne d in T able 6.1 ). MINRES and MINRES-QLP terminate when the r e curr e d NRBE is less than the given tol = 10 − 14 . Upp er left: η = 10 − 8 and thus κ ( A ) ≈ 10 8 . Also b = e and ther efor e k x k  k b k . The gr aphs of MINRES’s dire ctly c ompute d r esidual norms k r M k k and re curr ently c ompute d r esidual norms φ M k start to diﬀer at the level of 10 − 1 starting at iter ation 21, while the values φ Q k ≈ k r Q k k fr om MINRES-QLP de cr e ase monotonic al ly and stop ne ar 10 − 6 at iter ation 26. Upp er right: A gain η = 10 − 8 but b = Ae . Thus k x k = k e k = O ( k b k ) . The MINRES gr aphs of k r M k k and φ M k start to diﬀer when they r e ach a much smal ler level of 10 − 10 at iter ation 30. The MINRES-QLP φ Q k ’s ar e exc el lent appr oximations of φ Q k , with b oth r e aching 10 − 13 at iter ation 33. L ower left: η = 10 − 10 and thus A is even mor e il l-c onditione d than the matrix in the upp er plots. Here b = e and k x k is again explo ding. MINRES ends with k r M k k ≈ 10 2 , which me ans no c onver genc e, while MINRES-QLP r e aches a r esidual norm of 10 − 4 . L ower right: η = 10 − 10 and b = Ae . The ﬁnal MINRES r esidual norm k r M k k ≈ 10 − 8 , which is satisfactory but not as ac cur ate as φ M k claims at 10 − 13 . MINRES-QLP again has φ Q k ≈ k r Q k k ≈ 10 − 13 at iter ation 37. This ﬁgur e c an b e r epr o duc e d by DPtestSing7.m . 24 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS 5 10 15 20 10 −10 10 −5 10 0 k κ ≈ 10 2 ,  A  =3 ,  x ≈ 10 11 ,  b ≈ 10 1 ,t o l =1 0 − 9  Ar M k  ψ M k  Ar Q k  ψ Q k 5 10 15 20 25 10 −10 10 −5 10 0 κ ≈ 10 4 ,  A  =3 ,  x ≈ 10 10 ,  b ≈ 10 1 ,t o l =1 0 − 9 k 5 10 15 20 25 30 10 −10 10 −5 10 0 κ ≈ 10 6 ,  A  =3 ,  x ≈ 10 11 ,  b ≈ 10 1 ,t o l =1 0 − 9 k 5 10 15 20 25 30 10 −10 10 −5 10 0 κ ≈ 10 8 ,  A  =3 ,  x ≈ 10 11 ,  b ≈ 10 1 ,t o l =1 0 − 9 k Four (singular) symmetric least squares problems. Fig. 8.3 . Solving Ax = b with semideﬁnite A similar to an example of Sleijp en et al. [ 47 ]. A = Q diag([0 5 , η , 2 η , 2 : 1 789 : 3]) Q of dimension n = 797 with k A k = 3 , wher e Q = I − (2 /n ) ee T is a Householder matrix gener ate d by e = [1 , . . . , 1] T . (We are not plotting the NRBE quantities b e c ause k A kk r k k ≈ 6 thr oughout the iter ations in this example.) Upp er left: η = 10 − 2 and thus cond( A ) ≈ 10 2 . Also b = e and ther efor e k x k  k b k . The gr aphs of MIN RES’s dir e ctly compute d k Ar M k k and r e curr ently compute d ψ M k , and also ψ Q k ≈ k Ar Q k k fr om MINRES-QLP, match very wel l throughout the iter ations. Upp er right: Here, η = 10 − 4 and A is mor e il l-c onditione d than the last example (upp er left). The ﬁnal MINRES residual norm ψ M k ≈ k Ar M k k is slightly larger than the ﬁnal MINRES-QLP r esidual norm ψ Q k ≈ k Ar Q k k . The MINRES-QLP ψ Q k ar e exc el lent appr oximations of k Ar Q k k . L ower left: η = 10 − 6 and cond( A ) ≈ 10 6 . MINRES’s ψ M k and k Ar M k k diﬀer starting at iter ation 21. Eventual ly, k Ar M k k ≈ 3 , which me ans no c onverg enc e. MINRES-QLP r eaches a r esidual norm of ψ Q k = k Ar Q k k = 10 − 2 . L ower right: η = 10 − 8 . MINRES p erforms even worse than in the last example (lower left). MINRES-QLP re aches a minimum k Ar Q k k ≈ 10 − 7 but tol = 10 − 8 do es not shut it down so on enough. The ﬁnal ψ Q k = k Ar Q k k = 10 − 2 . The values of ψ Q k and k Ar Q k k diﬀer only at iter ations 27–28. This ﬁgur e c an b e r epr o duc e d by DPtestLSSing5.m . Ac kno wledgemen ts. W e thank Jan Modersitzki, Gerard Sleijp en, Henk V an der V orst, and also Kaustuv for providing us with their Ma tlab scripts, which ha v e aided us in pro ducing parts of Figures 6.1 , 8.2 , and 8.3 . W e also thank Michael F riedlander, Rasmus Larsen, and Lek-Heng Lim for their ﬁnest examples of work, discussion, and support. W e are most grateful to our anonymous reviewers for their insigh tful suggestions for improving this manuscript. Last but not least, w e dedicate this pap er to the memory of our colleague and friend, Gene Golub. MINRES-QLP FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS 25 App endix A. Proof that T ` is nonsingular iﬀ b ∈ range( A ) (section 2.1 ) . If T ` is nonsingular, we ha ve AV ` T − 1 ` e 1 = V ` e 1 = β − 1 1 b . Conv ersely , if b ∈ range( A ), then range( V ` ) ⊆ range( A ) and null( A ) ∩ range( V ` ) = { 0 } . W e also know that rank( V ` ) = ` and rank( T ` ) = rank( V ` T ` ) = rank( AV ` ) = rank( V ` ) − dim[null( A ) ∩ range( V ` )]; see [ 4 , F act 2.10.4 ii]. Th us rank( T ` ) = ` and so T ` is nonsingular.) REFERENCES [1] M. Arioli and S. Gratton. Least-squares problems, normal equations, and stopping criteria for the conjugate gradien t metho d. T echnical Report RAL-TR-2008-008, Rutherford Appleton Laboratory , Oxfordshire, UK, 2008. [2] Z.-Z. Bai and G.-Q. Li. Restrictiv ely preconditioned conjugate gradient methods for systems of linear equations. IMA J. Numer. A nal. , 23(4):561–580, 2003. [3] M. Benzi. Preconditioning techniques for large linear systems: a survey . J. Comput. Phys. , 182(2):418–477, 2002. [4] D. S. Bernstein. Matrix mathematics: theory, facts, and formulas . Princeton University Press, Princeton, NJ, 2nd edition, 2009. [5] E. Y. Bobrovnik ov a and S. A. V av asis. Accurate solution of weigh ted least squares by iterative methods. SIAM J. Matrix Anal. Appl. , 22(4):1153–1174, 2001. [6] A. M. Bradley . A lgorithms for the Equilibr ation of Matric es and Their Applic ation to Limite d- Memory Quasi-Newton Metho ds . PhD thesis, ICME, Stanford Universit y , 2010. [7] X.-W. Chang, C. C. Paige, and D. Titley-P eloquin. Stopping criteria for the iterative solution of linear least squares problems. Submitted to SIAM J. Matrix Anal. Appl. , 2008. [8] T.-Y. Chen and J. W. Demmel. Balancing sparse matrices for computing eigenv alues. Line ar Algebr a Appl. , 309(Issues 1-3):261–287, 2000. [9] S.-C. Choi, D. L. Donoho, A. G. Flesia, X. Huo, O. Levi, and D. Shi. Ab out Beamlab—a to olbox for new multiscale metho dologies. http://www- stat.stanford.edu/ ~ beamlab/ , 2002. [10] S.-C. T. Choi. Iter ative Methods for Singular Line ar Equations and L east-Squar es Pr oblems . PhD thesis, ICME, Stanford University , 2006. [11] J. Claerb out. Hyp ertext documents about reproducible research. http://sepwww.stanford. edu/doku.php?id=sep:research:reproducible . [12] F. A. Dul. MINRES and MINERR are better than SYMMLQ in eigenpair computations. SIAM J. Sci. Comput. , 19(6):1767–1782, 1998. [13] L. Eld´ en. Algorithms for the regularization of ill-conditioned least squares problems. Nor disk Tidskr. Informationsb ehand ling (BIT) , 17(2):134–145, 1977. [14] B. Fischer, A. Ramage, D. J. Silvester, and A. J. W athen. Minimum residual metho ds for augmented systems. BIT , 38(3):527–543, 1998. [15] R. W. F reund and N. M. Nach tigal. A new Krylov-subspace metho d for symmetric indeﬁnite linear systems. In W. F. Ames, editor, Pr o ce e dings of the 14th IMACS World Congr ess on Computational and Applied Mathematics , pages 1253–1256. IMACS, 1994. [16] K. A. Galliv an, S. Thirumalai, P . V an Do oren, and V. V ermaut. High p erformance algorithms for To eplitz and blo c k To eplitz matrices. In Pro c. F ourth Confer enc e of the International Line ar Algebr a So ciety (R otter dam, 1994) , volume 241/243, pages 343–388, 1996. [17] P . E. Gill, W. Murray , D. B. Poncele´ on, and M. A. Saunders. Preconditioners for indeﬁnite systems arising in optimization. SIAM J. Matrix Anal. Appl. , 13(1):292–311, 1992. [18] G. H. Golub and C. F. V an Loan. Matrix Computations . Johns Hopkins Univ ersity Press, Baltimore, MD, 3rd edition, 1996. [19] J. Hadamard. Sur les probl` emes aux d ´ eriv´ ees partielles et leur signiﬁcation physique. Princ eton University Bul letin , XI II(4):49–52, 1902. [20] W. W. Hager. Iterative metho ds for nearly singular linear systems. SIAM J. Sci. Comput. , 22(2):747–766, 2000. [21] M. Hanke and J. G. Nagy . Restoration of atmospherically blurred images by symmetric indef- inite conjugate gradient techniques. Inverse Pr oblems , 12(2):157–173, 1996. [22] P . C. Hansen. T runcated singular v alue decomp osition solutions to discrete ill-p osed problems with ill-determined numerical rank. SIAM J. Sci. Statist. Comput. , 11(3):503–518, 1990. [23] P . C. Hansen and D. P . O’Leary . The use of the L-curve in the regularization of discrete ill-p osed problems. SIAM J. Sci. Comput. , 14(6):1487–1503, 1993. [24] M. R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear systems. J. R ese ar ch Nat. Bur. Standar ds , 49:409–436, 1952. [25] N. J. Higham. A ccur acy and Stability of Nu meric al Algorithms . SIAM, Philadelphia, P A, 2nd 26 S.-C. T. CHOI, C. C. P AIGE, AND M. A. SAUNDERS edition, 2002. [26] P . Jir´ anek and D. Titley-Peloquin. Estimating the bac kward error in LSQR. SIAM J. Matrix Anal. Appl. , 31(4):2055–2074, 2010. [27] E. F. Kaasschieter. Preconditioned conjugate gradients for solving singular systems. J. Comput. Appl. Math. , 24(1-2):265–275, 1988. [28] Kaustuv. IPSOL: An Interior Point Solver for Nonc onvex Optimization Problems . PhD thesis, SCCM Program, Stanford Univ ersit y , 2008. [29] M. Kilmer and G. W. Stewart. Iterativ e regularization and MINRES. SIAM J. Matrix Anal. Appl. , 21(2):613–628, 1999. [30] C. Lanczos. An iteration metho d for the solution of the eigenv alue problem of linear diﬀeren tial and integral op erators. J. Rese ar ch Nat. Bur. Standar ds , 45:255–282, 1950. [31] R. M. Larsen. Eﬃcient Algorithms for Helioseismic Inversion . PhD thesis, Dept of Computer Science, Universit y of Aarhus, 1998. [32] J. G. Lewis. Algorithms for Sparse Matrix Eigenvalue Problems . PhD thesis, Dept of Computer Science, Stanford Universit y , 1976. [33] O. E. Livne and G. H. Golub. Scaling by binormalization. Numer. Algorithms , 35(1):97–120, 2004. [34] M. F. Murphy , G. H. Golub, and A. J. W athen. A note on preconditioning for indeﬁnite linear systems. SIAM J. Sci. Comput. , 21(6):1969–1972, 2000. [35] M. G. Neytchev a and P . S. V assilevski. Preconditioning of indeﬁnite and almost singular ﬁnite element elliptic equations. SIAM J. Sci. Comput. , 19(5):1471–1485, 1998. [36] J. No cedal and S. J. W right. Numeric al Optimization . Springer, New Y ork, 2nd edition, 2006. [37] Y. Nota y . Solving p ositiv e (semi)deﬁnite linear systems by preconditioned iterative methods. In Pre c onditione d Conjugate Gr adient Metho ds (Nijme gen, 1989) , volume 1457 of L e ctur e Notes in Math. , pages 105–125. Springer, Berlin, 1990. [38] C. C. P aige. Error analysis of the Lanczos algorithm for tridiagonalizing a symmetric matrix. J. Inst. Math. Appl. , 18(3):341–349, 1976. [39] C. C. Paige and M. A. Saunders. Solution of sparse indeﬁnite systems of linear equations. SIAM J. Numer. Anal. , 12(4):617–629, 1975. [40] C. C. P aige and M. A. Saunders. LSQR: an algorithm for sparse linear equations and sparse least squares. ACM T r ans. Math. Software , 8(1):43–71, 1982. [41] C. C. Paige and M. A. Saunders. Algorithm 583; LSQR: Sparse linear equations and least- squares problems. ACM T r ans. Math. Software , 8(2):195–209, 1982. [42] W.-Q. Ren and J.-X. Zhao. Iterative methods with preconditioners for indeﬁnite systems. J. Comput. Math. , 17(1):89–96, 1999. [43] M. Rozlo ˇ zn ´ ık and V. Simoncini. Krylo v subspace metho ds for saddle p oin t problems with indeﬁnite preconditioning. SIAM J. Matrix Anal. Appl. , 24(2):368–391 (electronic), 2002. [44] Y. Saad and M. H. Sc h ultz. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Statist. Comput. , 7(3):856–869, 1986. [45] M. A. Saunders. Solution of sparse rectangular systems using LSQR and Craig. BIT , 35(4):588– 604, 1995. [46] M. A. Saunders. Computing pro jections with LSQR. BIT , 37(1):96–104, 1997. [47] G. L. G. Sleijp en, H. A. V an der V orst, and J. Mo dersitzki. Diﬀerences in the eﬀects of rounding errors in Krylov solvers for symmetric indeﬁnite linear systems. SIAM J. Matrix Anal. Appl. , 22(3):726–751, 2000. [48] Systems Optimization Lab oratory (SOL), Stanford Universit y , downloadable soft ware. http: //www.stanford.edu/group/SOL/software.html . [49] G. W. Stewart. On the contin uity of the generalized inv erse. SIAM J. Appl. Math. , 17:33–45, 1969. [50] G. W. Stew art. Researc h, developmen t and LINP ACK. In J. R. Rice, editor, Mathematical Softwar e III , pages 1–14. Academic Press, New Y ork, 1977. [51] G. W. Stewart. The QLP approximation to the singular value decomposition. SIAM J. Sci. Comput. , 20(4):1336–1348, 1999. [52] K.-C. T oh, K.-K. Pho on, and S.-H. Chan. Block preconditioners for symmetric indeﬁnite linear systems. Internat. J. Numer. Metho ds Engr g. , 60(8):1361–1381, 2004. [53] L. N. T refethen and D. Bau, II I. Numeric al Line ar Algebr a . SIAM, Philadelphia, P A, 1997. [54] Univ ersit y of Florida sparse matrix collection. http://www.cise.ufl.edu/research/sparse/ matrices/ . [55] Y. Zhang. Solving large-scale linear programs by in terior-point metho ds under the MA TLAB environmen t. Optim. Metho ds Softw. , 10(1):1–31, 1998.

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

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment