Iterative Methods for Sparse Linear Systems

The front matter includes the title page, copyright page, TOC, preface to the second edition, and preface to the first edition

This chapter gives an overview of the relevant concepts in linear algebra that are useful in later chapters. It begins with a review of basic matrix theory and introduces the elementary notation used throughout the book. The convergence analysis of iterative methods requires a good level of knowledge in mathematical analysis and in linear algebra. Traditionally, many of the concepts presented specifically for these analyses have been geared toward matrices arising from the discretization of partial differential equations (PDEs) and basic relaxation-type methods. These concepts are now becoming less important because of the trend toward projection-type methods, which have more robust convergence properties and require different analysis tools. The material covered in this chapter will be helpful in establishing some theory for the algorithms and defining the notation used throughout the book.

1.1 Matrices

For the sake of generality, all vector spaces considered in this chapter are complex, unless otherwise stated. A complex n × m matrix A is an n × m array of complex numbers aij ,i=1,…,n,j=1,…,m. The set of all n × m matrices is a complex vector space denoted by ℂ^{n × m}. The main operations with matrices are the following:

• Addition: C = A + B, where A, B, and C are matrices of size n × m and cij = aij + bij ,i=1,2,…,n,j=1,2,…,m.

• Multiplication by a scalar: C = α A, where cij =α aij ,i=1,2,…,n,j=1,2,…,m.

Partial differential equations (PDEs) constitute by far the biggest source of sparse matrix problems. The typical way to solve such equations is to discretize them, i.e., to approximate them by equations that involve a finite number of unknowns. The matrix problems that arise from these discretizations are generally large and sparse; i.e., they have very few nonzero entries. There are several different ways to discretize a PDE. The simplest method uses finite difference approximations for the partial differential operators. The finite element method replaces the original function with a function that has some degree of smoothness over the global domain but is piecewise polynomial on simple cells, such as small triangles or rectangles. This method is probably the most general and well understood discretization technique available. In between these two methods, there are a few conservative schemes called finite volume methods, which attempt to emulate continuous conservation laws of physics. This chapter introduces these three different discretization methods.

2.1 Partial Differential Equations

Physical phenomena are often modeled by equations that relate several partial derivatives of physical quantities, such as forces, momentums, velocities, energy, temperature, etc. These equations rarely have a closed-form (explicit) solution. In this chapter, a few types of PDEs are introduced, which will serve as models throughout the book. Only one- and two-dimensional problems are considered, and the space variables are denoted by x in the case of one-dimensional problems and x₁ and x₂ for two-dimensional problems. In two dimensions, x denotes the “vector” of components (x₁, x₂).

As described in the previous chapter, standard discretizations of partial differential equations (PDEs) typically lead to large and sparse matrices. A sparse matrix is defined, somewhat vaguely, as a matrix with very few nonzero elements. But, in fact, a matrix can be termed sparse whenever special techniques can be utilized to take advantage of the large number of zero elements and their locations. These sparse matrix techniques begin with the idea that the zero elements need not be stored. One of the key issues is to define data structures for these matrices that are well suited for efficient implementation of standard solution methods, whether direct or iterative. This chapter gives an overview of sparse matrices, their properties, their representations, and the data structures used to store them.

3.1 Introduction

The natural idea to take advantage of the zeros of a matrix and their location was initiated by engineers in various disciplines. In the simplest case involving banded matrices, special techniques are straightforward to develop. Electrical engineers dealing with electrical networks in the 1960s were the first to exploit sparsity to solve general sparse linear systems for matrices with irregular structure. The main issue, and the first addressed by sparse matrix technology, was to devise direct solution methods for linear systems. These had to be economical in terms of both storage and computational effort. Sparse direct solvers can handle very large problems that cannot be tackled by the usual “dense” solvers.

Essentially, there are two broad types of sparse matrices: structured and unstructured. A structured matrix is one whose nonzero entries form a regular pattern, often along a small number of diagonals. Alternatively, the nonzero elements may lie in blocks (dense submatrices) of the same size, which form a regular pattern, typically along a small number of (block) diagonals. A matrix with irregularly located entries is said to be irregularly structured.

The first iterative methods used for solving large linear systems were based on relaxation of the coordinates. Beginning with a given approximate solution, these methods modify the components of the approximation, one or a few at a time and in a certain order, until convergence is reached. Each of these modifications, called relaxation steps, is aimed at annihilating one or a few components of the residual vector. Now these techniques are rarely used separately. However, when combined with the more efficient methods described in later chapters, they can be quite successful. Moreover, there are a few application areas where variations of these methods are still quite popular.

4.1 Jacobi, Gauss—Seidel, and Successive Overrelaxation

This chapter begins by reviewing the basic iterative methods for solving linear systems. Given an n × n real matrix A and a real n-vector b, the problem considered is as follows: Find x belonging to ℝⁿ such that Ax=b. 4.1 Equation (4.1) is a linear system, A is the coefficient matrix, b is the right-hand side vector, and x is the vector of unknowns. Most of the methods covered in this chapter involve passing from one iterate to the next by modifying one or a few components of an approximate vector solution at a time. This is natural since there are simple criteria when modifying a component in order to improve an iterate. One example is to annihilate some component(s) of the residual vector b − Ax. The convergence of these methods is rarely guaranteed for all matrices, but a large body of theory exists for the case where the coefficient matrix arises from the finite difference discretization of elliptic partial differential equations (PDEs).

Most of the existing practical iterative techniques for solving large linear systems of equations utilize a projection process in one way or another. A projection process represents a canonical way of extracting an approximation to the solution of a linear system from a subspace. This chapter describes these techniques in a very general framework and presents some theory. The one-dimensional case is covered in detail at the end of the chapter, as it provides a good preview of the more complex projection processes to be seen in later chapters.

5.1 Basic Definitions and Algorithms

Consider the linear system Ax=b, 5.1 where A is an n × n real matrix. In this chapter, the same symbol A is often used to denote the matrix and the linear mapping in ℝⁿ that it represents. The idea of projection techniques is to extract an approximate solution to the above problem from a subspace of ℝⁿ. If K is this subspace of candidate approximants, also called the search subspace, and if m is its dimension, then, in general, m constraints must be imposed to be able to extract such an approximation. A typical way of describing these constraints is to impose m (independent) orthogonality conditions. Specifically, the residual vector b − Ax is constrained to be orthogonal to m linearly independent vectors. This defines another subspace ℒ of dimension m, which will be called the subspace of constraints or left subspace for reasons that will be explained below. This simple framework is common to many different mathematical methods and is known as the Petrov-Galerkin conditions.

There are two broad classes of projection methods: orthogonal and oblique. In an orthogonal projection technique, the subspace ℒ is the same as K .

The next two chapters explore a few methods that are considered currently to be among the most important iterative techniques available for solving large linear systems. These techniques are based on projection processes, both orthogonal and oblique, onto Krylov sub-spaces, which are subspaces spanned by vectors of the form p(A)υ, where p is apolynomial. In short, these techniques approximate A⁻¹b by p(A)b, where p is a “good” polynomial. This chapter covers methods derived from, or related to, the Arnoldi orthogonalization. The next chapter covers methods based on Lanczos biorthogonalization.

6.1 Introduction

Recall from the previous chapter that a general projection method for solving the linear system Ax=b 6.1 extracts an approximate solution x_m from an affine subspace x0 + Km of dimension m by imposing the Petrov-Galerkin condition b−A xm ⊥ ℒm , where ℒ_m is another subspace of dimension m. Here, x₀ represents an arbitrary initial guess to the solution. A Krylov subspace method is a method for which the subspace Km is the Krylov subspace Km (A, r0 ) =span { r0 ,A r0 , A2 r0 ,…, Am−1 r0 } , where r₀ = b − Ax₀. When there is no ambiguity, Km (A, r0 ) will be denoted by Km . The different versions of Krylov subspace methods arise from different choices of the subspace ℒ_m and from the ways in which the system is preconditioned, a topic that will be covered in detail in later chapters.

The previous chapter considered a number of Krylov subspace methods that relied on some form of orthogonalization of the Krylov vectors in order to compute an approximate solution. This chapter will describe a class of Krylov subspace methods that are instead based on a biorthogonalization algorithm due to Lanczos. These are projection methods that are intrinsically nonorthogonal. They have some appealing properties, but are harder to analyze theoretically.

7.1 Lanczos Biorthogonalization

The Lanczos biorthogonalization algorithm is an extension to nonsymmetric matrices of the symmetric Lanczos algorithm seen in the previous chapter. One such extension, the Arnoldi procedure, has already been seen. However, the nonsymmetric Lanczos algorithm is quite different in concept from Arnoldi's method because it relies on biorthogonal sequences instead of orthogonal sequences.

7.1.1 The Algorithm

The algorithm proposed by Lanczos for nonsymmetric matrices builds a pair of biorthogonal bases for the two subspaces Km (A, υ1 ) =span { υ1 ,A υ1 ,…, Am−1 υ1 } and Km ( AT , w1 ) =span { w1 , AT w1 ,…, ( AT ) m−1 w1 } .

There are a number of techniques for converting a nonsymmetric linear system into a symmetric one. One such technique solves the equivalent linear system A^TAx = A^Tb, called the normal equations. Often, this approach is avoided in practice because the coefficient matrix A^TA is much worse conditioned than A. However, the normal equations approach may be adequate in some situations. Indeed, there are even applications in which it is preferred to the usual Krylov subspace techniques. This chapter covers iterative methods that are either directly or implicitly related to the normal equations.

8.1 The Normal Equations

In order to solve the linear system Ax = b when A is nonsymmetric, we can solve the equivalent system AT A x= AT b, 8.1 which is symmetric positive definite (SPD). This system is known as the system of the normal equations associated with the least-squares problem minimize ‖b−Ax‖ 2 . 8.2 Note that (8.1) is typically used to solve the least-squares problem (8.2) for overdetermined systems, i.e., when A is a rectangular matrix of size n × m, m < n.

A similar well-known alternative sets x = A^Tu and solves the following equation for u: A AT u=b. 8.3 Once the solution u is computed, the original unknown x can be obtained by multiplying u by A^T.

Although the methods seen in previous chapters are well founded theoretically, they are all likely to suffer from slow convergence for problems that arise from typical applications such as fluid dynamics and electronic device simulation. Preconditioning is a key ingredient for the success of Krylov subspace methods in these applications. This chapter discusses the preconditioned versions of the iterative methods already seen, but without being specific about the particular preconditioners used. The standard preconditioning techniques will be covered in the next chapter.

9.1 Introduction

Lack of robustness is a widely recognized weakness of iterative solvers relative to direct solvers. This drawback hampers the acceptance of iterative methods in industrial applications despite their intrinsic appeal for very large linear systems. Both the efficiency and robustness of iterative techniques can be improved by using preconditioning. A term introduced in Chapter 4, preconditioning is simply a means of transforming the original linear system into one with the same solution, but that is likely to be easier to solve with an iterative solver. In general, the reliability of iterative techniques, when dealing with various applications, depends much more on the quality of the preconditioner than on the particular Krylov subspace accelerators used. We will cover some of these preconditioners in detail in the next chapter. This chapter discusses the preconditioned versions of the Krylov subspace algorithms already seen, using a generic preconditioner.

To begin with, it is worthwhile to consider the options available for preconditioning a system. The first step in preconditioning is to find a preconditioning matrix M. The matrix M can be defined in many different ways but it must satisfy a few minimal requirements. From a practical point of view, the most important requirement for M is that it be inexpensive to solve linear systems Mx = b.

Finding a good preconditioner to solve a given sparse linear system is often viewed as a combination of art and science. Theoretical results are rare and some methods work surprisingly well, often despite expectations. A preconditioner can be defined as any subsidiary approximate solver that is combined with an outer iteration technique, typically one of the Krylov subspace iterations seen in previous chapters. This chapter covers some of the most successful techniques used to precondition a general sparse linear system. Note at the outset that there are virtually no limits to available options for obtaining good preconditioners. For example, preconditioners can be derived from knowledge of the original physical problems from which the linear system arises. However, a common feature of the preconditioners discussed in this chapter is that they are built from the original coefficient matrix.

10.1 Introduction

Roughly speaking, a preconditioner is any form of implicit or explicit modification of an original linear system that makes it easier to solve by a given iterative method. For example, scaling all rows of a linear system to make the diagonal elements equal to one is an explicit form of preconditioning. The resulting system can be solved by a Krylov subspace method and may require fewer steps to converge than the original system (although this is not guaranteed). As another example, solving the linear system M−1 Ax= M−1 b, where M⁻¹ is some complicated mapping that may involve fast Fourier transforms (FFT), integral calculations, and subsidiary linear system solutions, may be another form of preconditioning. Here, it is unlikely that the matrices M and M⁻¹ A can be computed explicitly. Instead, the iterative processes operate with A and with M⁻¹ whenever needed. In practice, the preconditioning operation M⁻¹ should be inexpensive to apply to an arbitrary vector.

Parallel computing has recently gained widespread acceptance as a means of handling very large computational tasks. Since iterative methods are appealing for large linear systems of equations, it is no surprise that they are the prime candidates for implementations on parallel architectures. There have been two traditional approaches for developing parallel iterative techniques thus far. The first extracts parallelism whenever possible from standard algorithms. The advantage of this viewpoint is that it is easier to understand in general since the underlying method has not changed from its sequential equivalent. The second approach is to develop alternative algorithms that have enhanced parallelism. This chapter will give an overview of implementations and will emphasize methods in the first category. The later chapters will consider alternative algorithms that have been developed specifically for parallel computing environments.

11.1 Introduction

Because of the increased importance of three-dimensional models and the high cost associated with sparse direct methods for solving these problems, iterative techniques play a major role in application areas. The main appeal of iterative methods is their low storage requirement. Another advantage is that they are far easier to implement on parallel computers than sparse direct methods because they only require a rather small set of computational kernels. Increasingly, direct solvers are being used in conjunction with iterative solvers to develop robust preconditioners.

The first considerations for high performance implementations of iterative methods involved implementations on vector computers. These efforts started in the mid-1970s, when the first vector computers appeared. Currently, there is a larger effort to develop new practical iterative methods that are not only efficient in a parallel environment, but also robust. Often, however, these two requirements seem to be in conflict.

This chapter covers a few alternative methods for preconditioning a linear system. These methods are suitable when the desired goal is to maximize parallelism. The simplest approach is the diagonal (or Jacobi) preconditioning. Often, this preconditioner is not very useful, since the number of iterations of the resulting iteration tends to be much larger than the more standard variants, such as incomplete LU (ILU) or symmetric successive overre-laxation (SSOR). When developing parallel preconditioners, one should be aware that the benefits of increased parallelism are not outweighed by the increased number of computations. The main question to ask is whether or not it is possible to find preconditioning techniques that have a high degree of parallelism, as well as good intrinsic qualities.

12.1 Introduction

As seen in the previous chapter, a limited amount of parallelism can be extracted from the standard preconditioners, such as ILU and SSOR. Fortunately, a number of alternative techniques can be developed that are specifically targeted at parallel environments. These are preconditioning techniques that would normally not be used on a standard machine, but only for parallel computers. There are at least three such types of techniques discussed in this chapter. The simplest approach is to use a Jacobi or, even better, a block Jacobi approach. In the simplest case, a Jacobi preconditioner may consist of the diagonal or a block diagonal of A. To enhance performance, these preconditioners can themselves be accelerated by polynomial iterations, i.e., a second level of preconditioning called polynomial preconditioning.

A different strategy altogether is to enhance parallelism by using graph theory algorithms, such as graph-coloring techniques. These consist of coloring nodes such that two adjacent nodes have different colors.

The convergence of preconditioned Krylov subspace methods for solving systems arising from discretized partial differential equations (PDEs) tends to slow down considerably as these systems become larger. This deterioration in the convergence rate, compounded with the increased operation count per step due to the sheer problem size, results in a severe loss of efficiency. In contrast, the class of methods to be described in this chapter is capable of achieving convergence rates that are, in theory, independent of the mesh size. One significant difference with the preconditioned Krylov subspace approach is that multigrid (MG) methods were initially designed specifically for the solution of discretized elliptic PDEs. The method was later extended in different ways to handle other PDE problems, including nonlinear ones, as well as problems not modeled by PDEs. Because these methods exploit more information on the problem than do standard preconditioned Krylov subspace methods, their performance can be vastly superior. On the other hand, they may require implementations that are specific to the physical problem at hand, in contrast with preconditioned Krylov subspace methods, which attempt to be general purpose.

13.1 Introduction

MG techniques exploit discretizations with different mesh sizes of a given problem to obtain optimal convergence from relaxation techniques. At the foundation of these techniques is the basic and powerful principle of divide and conquer. Though most relaxation-type iterative processes, such as Gauss-Seidel, may converge slowly for typical problems, it can be noticed that the components of the errors (or residuals) in the directions of the eigenvectors of the iteration matrix corresponding to the large eigenvalues are damped very rapidly. These eigenvectors are known as the oscillatory or high frequency modes.

As multiprocessing technology steadily gains ground, new classes of numerical methods that can take better advantage of parallelism are emerging. Among these techniques, domain decomposition methods are undoubtedly the best known and perhaps the most promising for certain types of problems. These methods combine ideas from partial differential equations (PDEs), linear algebra, and mathematical analysis, and techniques from graph theory. This chapter is devoted to decomposition methods, which are based on the general concepts of graph partitionings.

14.1 Introduction

Domain decomposition methods refer to a collection of techniques that revolve around the principle of divide and conquer. Such methods have been primarily developed for solving PDEs over regions in two or three dimensions. However, similar principles have been exploited in other contexts of science and engineering. In fact, one of the earliest practical uses for domain decomposition approaches was in structural engineering, a discipline that is not dominated by PDEs. Although this chapter considers these techniques from a purely linear algebra viewpoint, the basic concepts, as well as the terminology, are introduced from a model PDE.

Consider the problem of solving the Laplace equation on an L-shaped domain Ω partitioned as shown in Figure 14.1. Domain decomposition or substructuring methods attempt to solve the problem on the entire domain Ω= ⋃ i=1 s Ωi from problem solutions on the subdomains Ω_i. There are several reasons why such techniques can be advantageous. In the case of the above picture, one obvious reason is that the subproblems are much simpler because of their rectangular geometry.

The back matter includes bibliography and index