Iterative Solution of Symmetric Quasi-Definite Linear Systems

The front matter includes the title page, series page, copyright page, dedication, TOC, list of algorithms, list of theorems, and preface.

The solution of linear systems of equations has always been and remains a fundamental subproblem of countless computational methods for numerical science and engineering. The method of choice to solve systems is typically guided by their properties, including symmetry, definiteness, sparsity, conditioning, and block structure. Numerous applications, including constrained optimization and fluid dynamics, involve systems with the specific saddle-point block structure; i.e., their coefficient has the form

where the block F is square but need not be symmetric, the block G is also square and can be zero, and B₁ may or may not be equal to B₂. The special case of symmetric saddle-point systems arises when F is symmetric, B₁ = B₂, and G is either zero or symmetric and negative semidefinite. At the time of this writing, a large number of special cases of saddle-point systems, including the symmetric case, are understood sufficiently well that efficient numerical methods are known to solve them, though it is these authors' opinion that such systems still have much to reveal. The state-of-the-art paper by Benzi, Golub, and Liesen (2005) covers the main properties of general saddle-point systems and families of associated numerical methods. This book does not attempt to give a complete overview of saddle-point systems, of their properties, or of numerical methods, as they are the topic of ongoing research. Rather, we focus our efforts on a special case of symmetric saddle-point systems that arises when numerical methods are some how regularized. Thus, we recommend that the interested reader or student who is not yet familiar with saddle-point systems begin with Sections 1–9 of Benzi, Golub, and Liesen (2005) before going further.

2.1 ▪ Important Note about Notation

For the remainder of the book, we use a notation that slightly differs from that of the introductory chapter. The change of notation states that, rather than viewing an SQD system as following from a linear least-squares problem in which the variables are denoted x and the residuals r, we mean to view any SQD system, regardless of the application from which it arose, as the optimality conditions of a certain linear least-squares problem that can be deduced from the system. The philosophical difference is that the linear least-squares problem in question was not necessarily apparent in the original application. Thus, in order to focus on the SQD system, we denote its unknowns (x,y), and it will become apparent that it is y that plays the role of the variables in the least-squares problem, while x plays the role of the residuals.

This chapter covers some background material and known properties of SQD matrices. It also some what anticipates certain concepts and should be seen as a broad overview.

4.1 ▪ The Lanczos Process

In what follows, we often refer to Lanczos processes in various contexts. To establish the notation, consider a generic symmetric linear system with matrix H and right-hand side d. A Lanczos (1950, 1952) process applied to H and d constructs a sequence of vectors {s_k} according to the recursion

with s₀ := 0 and ω_k chosen at each iteration so that ‖s_k ‖₂ = 1. The notation in (4.1) is intentionally nonstandard to avoid confusion with scalar and vector quantities defined in the rest of the book, the latter adhering more closely to standard notation found in the literature. It is possible to derive (4.1) from the standard Gram-Schmidt orthogonalization process by exploiting the symmetry of H, and therefore the vectors s_k are theoretically orthonormal. After k iterations, the Lanczos process has generated k+1 vectors. We gather the first k Lanczos vectors into the matrix S_k := [s₁ … s_k].

This chapter presents a family of four methods based on the normal and Schur-complement equations. Three methods are generalizations of known methods in appropriate metrics: Lsqr, Craig, and Lsmr. The last one is new and may be viewed as an alternative to Lsmr when m < n. It also serves as an essential tool to explain the behavior of Minres on (1.5) in Section 6.4. For each method, we give implementation details in order to be complete and to provide a self-contained reference. The implementation details of Lsqr and Lsmr are our interpretation of the descriptions in Paige and Saunders (1982) and Fong and Saunders (2011). The implementation details of Craig were pieced together from various hints scattered across the literature, and those of Craig-mr are new, although they essentially mirror Lsmr.

In this chapter, we investigate the relation between the methods of Chapter 5 and two well-known Lanczos-based methods applied directly to the SQD system (1.5) with appropriate right-hand side: the conjugate-gradient method and the minimum residual method.

7.1 ▪ General Principle

One way to implement the methods of Chapter 5 is by modifying existing implementations of the standard Lsqr and Lsmr. Each code has two lines in its initialization that compute the first Golub-Kahan vectors ū₁ and . It suffices to replace those two lines with those of Algorithm 4.2 that compute u₁ and ν₁ and the two lines in the main loop that compute ū_j+1 and with the two lines in the main loop of Algorithm 4.2 using appropriate callbacks to solve systems with M and N.

In this final chapter in our journey through iterative methods for SQD systems, we provide a few final comments and also point the interested reader in promising research directions.

The back matter includes bibliography, index, and back cover.