3.1. Diagonal scaling

The earliest preconditioning ideas came from consideration of the first of these criteria, and are largely in the ‘given a matrix’ category. Taking $P$ to be a diagonal matrix is very simple; if in fact $P=\text{diag}(A)$ , then this is widely called Jacobi preconditioning. For a well-conditioned matrix that is just badly scaled, this is generally all that is needed, but in general diagonal preconditioning is likely to achieve very little in terms of reducing iterations or computation time except in specific situations as below. There are two related types of problem where Jacobi preconditioning (or just diagonal scaling) is important.

Firstly, the ‘identity’ matrix for finite element computation is the mass matrix, $Q$ , which is the Gram matrix of the finite element basis functions $\{{\it\phi}_{i},i=1,\ldots ,n\}$ , that is, it has entries

(3.4) $$\begin{eqnarray}Q=\{q_{i,j},i,j=1,\ldots ,n\},\quad q_{i,j}=\int _{{\rm\Omega}}{\it\rho}{\it\phi}_{i}(\mathbf{x}){\it\phi}_{j}(\mathbf{x})\end{eqnarray}$$

for a domain ${\rm\Omega}$ . Here the potentially variable coefficient ${\it\rho}$ might, for example, represent material density in a structure, and if this is highly variable or if the discretization has a large variation in the size (length, area, volume) of the elements – as is inevitable on an irregular or adapted mesh – then this variability is also represented in the diagonal entries of $Q$ . Thus $P=\text{diag}(Q)$ is provably a good preconditioner in this case (Wathen Reference Wathen1986). For example the gallery(’wathen’,n,n) matrix in matlab is a finite element mass matrix with randomly variable density for a so-called serendipity finite element. The precise eigenvalue bounds of the form (3.2) or (3.3), together with tabulated values for $\ell$ and ${\rm\Upsilon}$ for many types of finite elements, are derived in Wathen (Reference Wathen1986).

We briefly comment that in problems as above with analytic a priori eigenvalue bounds $\ell$ and ${\rm\Upsilon}$ which are (reasonably) tight, the possibility also arises to use a fixed but small number of Chebyshev (semi-)iterations as a (linear) preconditioner; see Wathen and Rees (Reference Wathen and Rees2009).

The second important situation arises also through highly variable coefficients; the classic example arises from the operator

$$\begin{eqnarray}{\rm\nabla}\cdot (K{\rm\nabla}u),\end{eqnarray}$$

where the permeability $K$ has large jumps. For continuous finite element approximation of such terms, Graham and Hagger (Reference Graham and Hagger1999) were the first to identify the utility of diagonal scaling to render the preconditioned problem more homogeneous. For such an elliptic PDE, Jacobi preconditioning is highly unlikely to be an adequate preconditioner on its own, but when used together with an appropriate preconditioner for a homogeneous elliptic PDE operator (see below), the combination can be much better than just the elliptic preconditioner by itself. In the literature there are several developments of the idea of Graham and Hagger, including for non-isotropic permeabilities.

3.2. Incomplete Cholesky factorization

It is well known that Cholesky factorization ( $LL^{T}$ or $LDL^{T}$ factorization where $L$ is lower triangular and $D$ is diagonal) of a symmetric and positive definite matrix leads immediately via forward and backward substitution to direct solution of a matrix system. For sparse matrices, unfortunately, the factor $L$ is usually much less sparse (i.e., it has many more non-zeros) than $A$ . Such fill-in is usually unavoidable, though widely used sparse direct methods use sophisticated algorithms in attempting to minimize it (HSL 2013, Davis Reference Davis2004, Amestoy, Duff, L’Excellent and Koster Reference Amestoy, Duff, L’Excellent and Koster2000, Li Reference Li2005). In the context of PDE problems, the essential rule is that a sparse direct method can be a good solver for a problem on a two-dimensional spatial domain, but is generally computationally infeasible for a three-dimensional problem. If one preselects (or dynamically selects) the sparsity pattern of $L$ and drops terms which would fill outside this pattern, then one can form an incomplete Cholesky factorization

$$\begin{eqnarray}A=LL^{T}+E,\end{eqnarray}$$

where $E$ represents the error due to the dropped terms. One can additionally threshold so that small entries are not kept in the incomplete factors. With $P=LL^{T}$ , it is clear that forward and backward substitution still enable easy solution of systems with the preconditioner, and the main question remaining is whether the preconditioner leads to fast convergence of cg.

In general, many different variants of this basic idea exist, and software for various incomplete factorization-based algorithms is widely available. This is certainly a plus-point, but unless $A$ is diagonally dominant, incomplete Cholesky algorithms do not always work well, indeed they do not always exist without modification; see Meijerink and van der Vorst (Reference Meijerink and van der Vorst1977), and see Benzi (Reference Benzi2002) for discussion of this as well as many other aspects of incomplete factorization preconditioning. With regard to robustness, there are certain advantages to computing an incomplete triangular factorization not using Cholesky but via orthogonalization in the inner product generated by $A$ ; see Benzi and Tuma (Reference Benzi and Tuma2003). This approach is called RIF (robust incomplete factorization) and shares some similarity with the AINV sparse approximate inverse discussed below.

Let us make some of these statements precise in the context of solving a simple discretized elliptic boundary value problem, namely the common 5-point finite difference approximation of the Dirichlet problem for the (negative) Laplacian

$$\begin{eqnarray}-{\rm\nabla}^{2}u=f\quad \text{in}\;{\rm\Omega},\quad u~\text{given on}\;\partial {\rm\Omega},\end{eqnarray}$$

which arises from second-order central differencing on a regular grid with spacing $h$ in each of the two spatial dimensions. The choice of sign ensures that the operator and thus the arising matrix are positive rather than negative definite. If the domain, ${\rm\Omega}$ , is the unit square (so the boundary $\partial {\rm\Omega}$ is four coordinate-aligned line segments), then the arising matrix $A$ is of dimension $N\approx h^{-2}$ for small $h$ and also ${\it\kappa}^{A}={\mathcal{O}}(h^{-2})$ ; see Elman et al. (Reference Elman, Silvester and Wathen2014b , Section 1.6). Thus (3.1) predicts that ${\mathcal{O}}(h^{-1})$ cg iterations would be required for convergence of unpreconditioned cg. The simplest incomplete Cholesky preconditioner for which $L$ has the same sparsity pattern as the lower triangular part of $A$ also leads to ${\it\kappa}^{P^{-1}A}={\mathcal{O}}(h^{-2})$ , though there is significant clustering of a large number of the eigenvalues around 1; the bound (3.1) based only on the extreme eigenvalues is pessimistic in this situation. A variant, modified incomplete factorization (due to Gustafsson Reference Gustafsson1978 but related to rather earlier work by Dupont, Kendall and Rachford Reference Dupont, Kendall and Rachford1968) reduces this to ${\it\kappa}^{P^{-1}A}={\mathcal{O}}(h^{-1})$ , and generally leads to slightly faster cg convergence. Applying (3.1) in this case predicts that ${\mathcal{O}}(h^{-1/2})$ preconditioned cg iterations would be required for convergence, and this is close to what is seen in practical computation. This implies that for a sequence of such problems with smaller and smaller $h$ , it will take an increasing number of iterations to achieve acceptable convergence as $h$ gets smaller (and so the matrix dimension gets larger) – a rather undesirable feature.

Fortunately, there are much better preconditioners for this and related elliptic PDE problems which guarantee ${\it\kappa}^{P^{-1}A}={\mathcal{O}}(1)$ , so that the number of required cg iterations is essentially constant for larger-dimensional discrete problems coming from approximation on finer grids. Such preconditioning methods are said to have $P$ spectrally equivalent to $A$ . The most important class of such methods arises from multigrid ideas, which we consider next.

Of course, even in the more general context of problems not coming from PDEs, incomplete Cholesky and its variants can be of considerable utility; in particular, algebraic preconditioners of this form can very often be applied.

Before moving on, we comment that preconditioning techniques have been introduced that turn out to be related to incomplete factorization. A particular example is the ‘support graph’ preconditioner introduced originally in a conference talk by Vaidya, who wrote commercial software but did not publish on the subject. The support graph idea comes out of graph-theoretic notions and has a somewhat discrete mathematical flavour. Later analysis both clarified the scope of the ideas – they are now essentially extended to cover most symmetric and diagonally dominant matrices – and made connections with incomplete factorization; see Bern et al. (Reference Bern, Gilbert, Hendrickson, Nguyen and Toledo2006), Spielman and Teng (Reference Spielman and Teng2014) and references therein.

3.3. Multigrid

Since the groundbreaking paper by Brandt (Reference Brandt1977), multigrid methods have been at the core of scientific computing. Earlier contributions, notably by Fedorenko (Reference Fedorenko1964) and Bakhvalov (Reference Bakhvalov1966), are now considered to have introduced the key ideas of combining a simple smoothing iteration (such as Jacobi or Gauss–Seidel iteration) with coarse grid correction on successively coarser grids for the finite difference Laplacian, but Brandt’s vision was much broader: he showed just how fast multigrid solvers could be and introduced many ideas to keep researchers busy for decades! For a positive definite elliptic PDE operator, it has been clear for a long time that an appropriate multigrid method will provide an optimal solver, that is, a solver such that the work to compute the solution scales linearly with the dimension of the discretized problem.

Classically, multigrid is a stationary iteration (albeit a non-trivial one) for solving $Ax=b$ when $A$ is a discrete Laplacian matrix (arising from finite difference or finite elements). For V- or W-cycles, iterates satisfy single-cycle contraction of the form

(3.5) $$\begin{eqnarray}\Vert x-x_{k}\Vert _{A}\leq {\it\eta}\;\Vert x-x_{k-1}\Vert _{A},\end{eqnarray}$$

where the contraction factor, ${\it\eta}$ , is typically around 0.1 (see, e.g., Elman et al. (Reference Elman, Silvester and Wathen2014b , Section 2.5)). It follows that the linear multigrid operator (represented by a complicated and fairly dense matrix which one would never want to use in computations) is an excellent candidate for a preconditioner, since from (3.5) it readily follows that

$$\begin{eqnarray}1-{\it\eta}\leq \frac{x^{T}Ax}{x^{T}Px}\leq 1+{\it\eta}\end{eqnarray}$$

(Elman et al. Reference Elman, Silvester and Wathen2014b , Lemma 4.2). That is, applying one multigrid V-cycle, say, is equivalent to solution of a preconditioner system with preconditioner $P$ which guarantees rapid cg convergence because ${\it\kappa}^{P^{-1}A}\leq (1+{\it\eta})/(1-{\it\eta})$ . There is little point in using multigrid cycles as a preconditioner for cg in this way for simple problems: one might as well simply use multigrid iteration alone. For more complicated problems, however, the observation that a single multigrid cycle is an efficient spectrally equivalent approximation of an elliptic PDE operator such as the Laplacian is very useful; see Section 5. Note also the important point that arises here for the first time in this review: a preconditioner can be a linear operation rather than specifically a known matrix. Thus $\{\text{3 V-cycles}\}$ , for example, is a suitable preconditioner. On the other hand, iterating multigrid (or any other iteration) to some residual tolerance will probably ensure that different numbers of V-cycles are used at different Krylov subspace iterations and so the preconditioner is not a fixed linear operator, and the above convergence analysis cannot apply.

Multigrid – and multilevel methodology more generally – has now become a huge subject (see for example Trottenberg, Oosterlee and Schüler Reference Trottenberg, Oosterlee and Schüler2000), with specialized techniques for the treatment of several specific features which might arise in a PDE (or other) problem. In the past couple of decades, algebraic multigrid methods (AMG), which try to mimic a geometric multigrid approach but use only information from the entries of a matrix rather than geometric information regarding grids, boundary conditions, etc., have attracted much attention, and there is now widely available software, for example HSL_MI20 (Boyle, Mihajlovic and Scott Reference Boyle, Mihajlovic and Scott2010), ML (Gee et al. Reference Gee, Siefert, Hu, Tuminaro and Sala2006) and BoomerAMG (Henson and Yang Reference Henson and Yang2000). In particular, in situations where a problem is fundamentally of the character of an elliptic differential operator but there are complicating features, an AMG cycle can be an excellent preconditioner. This can be true even if AMG is non-convergent, since it may contract on all but a few modes; the few corresponding outlier eigenvalues are easily dealt with by cg. It can simply be much more convenient to employ an AMG cycle as (part of) a preconditioner for a Krylov subspace method than to have to ‘hard-code’ bespoke techniques for each different variant problem that comes up. So common is this combined use of AMG and a Krylov subspace method that it is rarely even mentioned that AMG is used as a preconditioner, but this is almost always the case in practice.

We should mention that multigrid ideas have been tried in contexts other than just for PDE problems; see for example Giles (Reference Giles2008) and Livne and Brandt (Reference Livne and Brandt2012).

It should also be mentioned in this section that a multigrid V- or W-cycle can usually be quite easily made into a symmetric linear operator (see Elman et al. Reference Elman, Silvester and Wathen2014b , page 98), though the so-called full multigrid cycle, for example, cannot. Positive definiteness is not usually an issue for an underlying positive definite problem.

3.4. Domain decomposition

For a very large-dimensional problem, it might be a fairly natural idea to split into separate subproblems and solve each separately (and usually in parallel) to give an approximation to the solution of the overall problem. In the original context of solving an elliptic boundary value problem on a domain ${\rm\Omega}$ , this is the basic idea of domain decomposition: the subproblems being problems expressed on subdomains ${\rm\Omega}_{k},k=1,\ldots ,N$ with $\cup {\rm\Omega}_{k}={\rm\Omega}$ . In its simplest form, it corresponds to a block diagonal preconditioner with each separate subproblem represented by a single diagonal block which may or may not be overlapping, that is, which may or may not share some variables with other subproblems. An important issue is how to treat the variables either in the overlap regions or on the introduced boundaries between subdomains. Domain decomposition was first employed abstractly by Hermann Schwarz in the nineteenth century to prove an existence result for an elliptic PDE, hence methods of this type are often called Schwarz methods.

Like multigrid, domain decomposition is a huge subject which we have no intention of covering with any depth here; see for example the books by Quarteroni and Valli (Reference Quarteroni and Valli1999), Smith, Bjørstat and Gropp (Reference Smith, Bjørstat and Gropp1996), Toselli and Widlund (Reference Toselli and Widlund2004), Mathew (Reference Mathew2008) or the earlier survey by Chan and Mathew (Reference Chan and Mathew1994). What is clear, however, even from this brief description is that domain decomposition gives rise to preconditioners that are well suited to parallel computation.

To give a brief idea of what can be achieved on a simple problem, we consider the Poisson equation on a domain ${\rm\Omega}\subset \mathbb{R}^{d}$ approximated by standard (conforming) finite elements or finite differences, with a mesh spacing $h$ . If ${\rm\Omega}$ is split into non-overlapping subdomains of width $H$ , then there are basic block diagonal preconditioners without a coarse grid solve for which one can establish that

$$\begin{eqnarray}{\it\kappa}^{P^{-1}A}\leq CH^{-2}(1+\log ^{2}(H/h)).\end{eqnarray}$$

The more sophisticated domain decomposition schemes which also include a coarse grid (or ‘wire basket’ solve) can achieve

$$\begin{eqnarray}{\it\kappa}^{P^{-1}A}\leq C(1+\log ^{2}(H/h)).\end{eqnarray}$$

In either case, $C$ is a generic constant not dependent on $H,h$ , and the bounds translate into convergence estimates via (3.1). The size $H$ can be considered to be the coarse grid size in essence; the bounds above assume exact solution of the problems on subdomains.

With the importance of parallel and distributed computing, the domain decomposition concept remains very important and much is known for more advanced methods and for different problems; see the regular Domain Decomposition conference proceedings available at www.ddm.org (2014).

3.5. Sparse approximate inverses

For any matrix norm, minimization of

$$\begin{eqnarray}\Vert I-AR\Vert\end{eqnarray}$$

over all matrices $R\in \mathbb{R}^{n\times n}$ will clearly yield $R=A^{-1}$ . Minimization with restrictions on the entries where $R$ may have non-zero values leads to sparse approximate inverses. That is,

$$\begin{eqnarray}R=\text{argmin}_{R\in S_{R}}\Vert I-AR\Vert ,\end{eqnarray}$$

where $S_{R}$ describes the subset of entries of $R$ which are allowed to be non-zero, is a sparse approximate inverse for $A$ . To be consistent with the nomenclature in this article, $P=R^{-1}$ is the preconditioner here; solution of systems with $P$ is, of course, simply multiplication by $R$ in this situation. $S_{R}$ is strictly a subspace of matrices with restricted sparsity pattern.

The Frobenius norm

$$\begin{eqnarray}\Vert B\Vert _{F}^{2}=\mathop{\sum }_{i=1,\ldots ,n}~\mathop{\sum }_{j=1,\ldots ,n}b_{i,j}^{2}\end{eqnarray}$$

is by far the most widely considered, since the minimization problem then reduces to least-squares problems for the columns of $R$ separately (thus computable in parallel). Further, these least-squares problems are all of small dimension when the specification of $S_{R}$ ensures that $R$ is a sparse matrix. We again refer to Benzi (Reference Benzi2002) for a thorough description.

As posed in this simple form, the sparse approximate inverse, $R$ , is not generally symmetric when $A$ is symmetric. However, we introduce this type of preconditioner in this section because symmetry can be enforced by computing $R$ in factored form $R=LL^{T}$ , where $L$ is sparse and lower triangular; see Kolotilina and Yeremin (Reference Kolotilina and Yeremin1993). Symmetry of a sparse approximate inverse can also be achieved by the AINV algorithm of Benzi, Meyer and Tuma (Reference Benzi, Meyer and Tuma1996) and the more robust (stabilized) SAINV variant (Benzi, Cullum and Tuma Reference Benzi, Cullum and Tuma2000), which is guaranteed not to break down if $A$ is positive definite. See again Benzi (Reference Benzi2002).

The main issue in the construction of any sparse approximate inverse is the choice of the sparsity pattern $S_{R}$ . The simplest idea is to use the sparsity of $A$ , whilst there have been several attempts to dynamically define $S_{R}$ , most notably in the SPAI algorithm of Grote and Huckle (Reference Grote and Huckle1997). An extension which can be used in a number of different scenarios is to use a target matrix, $T$ , and hence minimize $\Vert T-AR\Vert _{F}$ ; see Holland, Wathen and Shaw (Reference Holland, Wathen and Shaw2005). A further comment is that sparse approximate inverse algorithms are clearly quite widely applicable, but generally do not give spectrally equivalent preconditioners for elliptic PDEs. They do, however, make excellent smoothers for use with multigrid (Tang and Wan Reference Tang and Wan2000, Bröker, Grote, Mayer and Reusken Reference Bröker, Grote, Mayer and Reusken2001).

3.6. Hierarchical representations

For symmetric and positive definite matrices arising via low-order finite element approximation from elliptic PDEs, Yserentant (Reference Yserentant1986) introduced the hierarchical basis preconditioner nearly three decades ago. This uses the idea that the coefficient matrix would have a much smaller condition number if a different (but less computationally convenient) basis for the approximation space were employed. By using a hierarchical basis – which would give rise to a fairly dense matrix – both this good conditioning and a rapid change of basis matrix based on a simple tree structure allow the sparse matrix arising from the usual locally defined finite element basis functions to be effectively preconditioned.

For the discrete Laplacian on a two-dimensional domain, Yserentant’s hierarchical preconditioner reduces ${\it\kappa}={\mathcal{O}}(h^{-2})$ to ${\it\kappa}={\mathcal{O}}(\log h^{-1})$ for the discrete Laplacian (see Yserentant Reference Yserentant1986 or Elman and Zhang Reference Elman and Zhang1995 for a more linear algebraic proof). Unfortunately, for three-dimensional domains only ${\it\kappa}={\mathcal{O}}(h^{-1})$ is possible with the straightforward use of the hierarchical basis (Ong Reference Ong1997), though an improved variant of this idea due to Bramble, Pasciak and Xu (Reference Bramble, Pasciak and Xu1990) can yield ${\it\kappa}={\mathcal{O}}(1)$ in three dimensions as well. These ideas, originally introduced as based on the multilevel splitting of finite element spaces, were the precursor to a whole abstract theory based on such splittings; the review article by Xu (Reference Xu1992) remains a largely comprehensive description of this abstract theory, which covers many approaches including multigrid and domain decomposition discussed above for symmetric and positive definite matrices arising from elliptic PDEs.

A different hierarchical structure can also be of utility. We wish to mention some ideas and the associated computational technology which enables certain dense matrices – such as the inverse matrices of discrete positive definite elliptic PDE operators – to be represented to machine accuracy in a highly efficient format. There are a number of related but differing approaches (Chandrasekaran et al. Reference Chandrasekaran, Dewilde, Gu, Pals, Sun, van der Veen and White2005, Greengard and Rokhlin Reference Greengard and Rokhlin1987, Hackbusch Reference Hackbusch1999, Grasedyck and Hackbusch Reference Grasedyck and Hackbusch2003, Hackbusch, Khoromskij and Sauter Reference Hackbusch, Khoromskij, Sauter, Bungartz, Hoppe and Zenger2000, Bebendorf Reference Bebendorf2008), but the essential observation can be thought of in terms of the Green’s function for an elliptic operator: nearby entries need accurate representation, but variables distant from a given entry can be treated together. This notion is most directly represented in the ${\mathcal{H}}$ -matrix representation of Hackbusch: off-diagonal blocks are represented as low-rank blocks using outer products of vectors or matrices with very few columns. By employing a blocking hierarchy, the computational complexity of storing (and also multiplying with) such matrix approximations can be as low as ${\mathcal{O}}(n\log n)$ .

We mention these ideas here since such approximate representations of dense matrices have been suggested as preconditioners rather in the same vein as sparse approximate inverses. They can also be utilized in other more demanding situations, as we describe later.

5.1. Shifted elliptic operators: Helmholtz problems

It might seem that taking a symmetric and positive definite matrix, $K$ , and simply shifting it by subtracting some multiple of the identity should lead to matrices for which easy solution algorithms exist. In general, however, this is not the case. Certainly the matrix $K-{\it\sigma}I$ can have resonances (and thus be singular) when ${\it\sigma}$ is an eigenvalue of $K$ , but the difficulty lies deeper than just this observation.

The classic (and important) PDE problem here is the Helmholtz equation $-{\rm\nabla}^{2}u-{\it\sigma}u=f$ , which arises when looking for periodic (wave-like) solutions of the linear wave equation; the parameter ${\it\sigma}$ is then related to the frequency. Supposing that $K$ is the discretized form of the positive definite operator $-{\rm\nabla}^{2}$ (i.e., the negative Laplacian), then for ${\it\sigma}>{\it\lambda}_{\min }(K)$ it is clear that

(5.4) $$\begin{eqnarray}A=K-{\it\sigma}I\end{eqnarray}$$

is indefinite and symmetric. Many of the important physical applications involve high frequencies and thus large values of ${\it\sigma}$ . Having a good preconditioner for $K$ – we have identified many in the sections above – is of little help in constructing a preconditioner for $A$ . For finite element discretization, the identity operator becomes the mass matrix $Q$ (3.4), but via diagonal approximation or since the grids employed are generally fairly regular for such wave problems, there is little lost in description by considering (5.4).

In fact, for Helmholtz problems there is usually a radiation boundary condition for outgoing waves or an absorbing boundary condition (or ‘perfectly matched layer’) is applied so that there is little reflection of outgoing waves back into the computational domain. These add some complex entries in the matrix which generally render it complex and symmetric (not Hermitian). Nevertheless, it is the indefiniteness of Helmholtz problems which remains the most important feature, hence why we consider them in this section. This aspect can, however, mean that non-symmetric Krylov subspace iterative methods are required; sqmr is quite popular in this context.

Magolu Monga Made, Beauwens and Warzée (Reference Magolu Monga Made, Beauwens and Warzée2000) and Laird and Giles (Reference Laird and Giles2002) were perhaps the first to recognize the potential of preconditioning (5.4) with a matrix of the form $K+{\it\mu}I$ where ${\it\mu}$ is positive. In the case of real ${\it\mu}$ , such a positive definite matrix is readily approximated by standard multigrid cycles which can be used as a positive definite preconditioner for (5.4). Complex values for ${\it\mu}$ have also been explored (see also van Gijzen, Erlangga and Vuik Reference van Gijzen, Erlangga and Vuik2007, for example), and good choices of the parameter can be found from Fourier analysis of the spectrum of the preconditioned operator/matrix. As well as indicating the speed of Krylov subspace iteration convergence, such analysis can also assist in the choice of multigrid components (Erlangga, Oosterlee and Vuik Reference Erlangga, Oosterlee and Vuik2006). Elman, Ernst and O’Leary (Reference Elman, Ernst and O’Leary2001) also propose a multigrid approach for Helmholtz problems employing gmres as a smoother. Domain decomposition approaches have also been suggested (Gander, Magoulès and Nataf Reference Gander, Magoulès and Nataf2002).

A more recent proposal by Engquist and Ying (Reference Engquist and Ying2011a ) for variable and high-frequency regimes involves ‘sweeping’ from the absorbing boundary, eliminating the variables locally in layers. The resulting preconditioner is thus a block incomplete $LDL^{T}$ factorization. A key aspect of the implementation is that hierarchical ( ${\mathcal{H}}$ -)matrix compression techniques as described in Section 3 above can be employed for the intermediate layers; this leads to significant reductions in storage and computational work and leads to excellent overall complexities such as ${\mathcal{O}}(n\log ^{2}n)$ for computation of $P=LDL^{T}$ and application of $P^{-1}$ to a vector for a two-dimensional domain problem. Krylov subspace iterative convergence is observed to occur in few iterations, making this a very competitive technique. Unfortunately, although the idea – and a closely related one employing frontal sparse factorization (Engquist and Ying Reference Engquist and Ying2011b ) – is observed to work well for three-dimensional domains as well, the ${\mathcal{H}}$ -matrix compression does not appear to yield the same small computational complexity: algebraic factors in $n$ replace logarithmic factors in estimates of the complexity for setting up the preconditioner (Poulson, Engquist, Li and Ying Reference Poulson, Engquist, Li and Ying2013). That is, for a three-dimensional domain problem, the computation of $P=LDL^{T}$ takes time which grows like $n^{{\it\gamma}}$ for some constant ${\it\gamma}>1$ whereas, once $P$ has been computed, the action of $P^{-1}$ can be computed in ${\mathcal{O}}(n\log n)$ work at each iteration (similarly to the FFT). For a problem with many different right-hand side vectors there could thus be some amelioration in the work required to set up the preconditioner.

We should alert the reader to the article by Ernst and Gander (Reference Ernst, Gander, Graham, Hou, Lakkis and Scheichl2012), which explains the general difficulty of iterative solution of Helmholtz problems.

The Maxwell equations describing electromagnetics have wave-like solutions and hence some of the character of Helmholtz problems. However, they also (in various formulations) have additional semidefinite features, as discussed in the section above, and additional indefinite saddle point structure, and are therefore considered below.

The idea of sweeping through the discrete variables for a problem in a particular order is also important for preconditioner construction for non-symmetric matrices which derive from PDEs involving advection/convection (or transport), as described below.

5.2. Saddle point systems

The second main class of symmetric and indefinite problems come from minimization of quadratic functionals subject to linear constraints. In fact, linearizations of much more general constrained optimization problems quite generally result in linear systems of the same form,

(5.5) $$\begin{eqnarray}{\mathcal{A}}x=\left[\begin{array}{@{}cc@{}}A & B^{T}\\ B & 0\end{array}\right]\left[\begin{array}{@{}c@{}}u\\ p\end{array}\right]=\left[\begin{array}{@{}c@{}}f\\ g\end{array}\right]\!,\end{eqnarray}$$

which are called saddle point systems; see Benzi, Golub and Liesen (Reference Benzi, Golub and Liesen2005) for a comprehensive treatment of such problems. The block $A\in \mathbb{R}^{n\times n}$ is usually positive definite on the kernel of $B\in \mathbb{R}^{m\times n}$ with $m<n$ . This is a sufficient condition for invertibility of ${\mathcal{A}}$ provided $B$ is of full rank. Often $A$ is positive definite, though in the context of sequential quadratic approximations in optimization (and possibly elsewhere), indefiniteness can arise; see Gould and Simoncini (Reference Gould and Simoncini2009) for spectral bounds in this less usual case. The variables $u$ are the primary variables and $p$ are the Lagrange multipliers. Sometimes the zero block is replaced with a regularization matrix $-C\in \mathbb{R}^{m\times m}$ , where $C$ is symmetric and positive semidefinite; if $C\neq 0$ , then some authors refer to a regularized or generalized saddle point problem.

Many problems give rise to matrices of this type: as well as the ubiquitous appearance of such problems in optimization, several important PDE problems, including the Stokes problem describing slow viscous incompressible flow, are of this form. The incompressible Navier–Stokes equations give matrices of a similar form but in which $A$ is non-symmetric. For these applications, the primary variables are the fluid velocity and the Lagrange multiplier variables are the fluid pressure. For other applications and references, see again Benzi et al. (Reference Benzi, Golub and Liesen2005) or Pestana and Wathen (Reference Pestana and Wathen2015b ).

The appearance of the zero diagonal block in (5.5) causes problems for most of the preconditioning approaches identified in Section 3 above. This is also generally true with regularization, $C\neq 0$ . Multigrid approaches, for example, require more complex smoothers; see for example John and Tobiska (Reference John and Tobiska2000) and Manservisi (Reference Manservisi2006). Domain decomposition can be applied, but then the subdomain problems remain indefinite. Incomplete factorization ideas can sometimes be employed if, for example, $A$ and $C$ are diagonal or nearly so (Bergamaschi, Gondzio and Zilli Reference Bergamaschi, Gondzio and Zilli2004).

A simple observation in the case that $A$ is non-singular is that if

(5.6) $$\begin{eqnarray}P=\left[\begin{array}{@{}cc@{}}A & 0\\ 0 & S\end{array}\right]\!,\quad S=BA^{-1}B^{T},\end{eqnarray}$$

then the preconditioned matrix $P^{-1}{\mathcal{A}}$ has only three distinct eigenvalues independently of $m,n$ ; thus minres (or gmres in the non-symmetric case) will converge in three iterations to the solution (Murphy, Golub and Wathen Reference Murphy, Golub and Wathen2000, Korzak Reference Korzak1999). $P$ in (5.6) is rarely a practical preconditioner for ${\mathcal{A}}$ in its exact form, but this observation does give a paradigm for the construction of preconditioners and, in particular, highlight the utility of approximations to the Schur complement, $S$ . An important point to note is that $P$ is nothing like an inverse of ${\mathcal{A}}$ (nor an approximation thereof). Nevertheless, good approximations (preconditioners) $\widehat{A}$ for the block $A$ and $\widehat{S}$ for the Schur complement, $S$ , lead to block diagonal preconditioners

(5.7) $$\begin{eqnarray}P=\left[\begin{array}{@{}cc@{}}\widehat{A} & 0\\ 0 & \widehat{S}\end{array}\right]\end{eqnarray}$$

for ${\mathcal{A}}$ , for which minres converges in few iterations. Precisely, if

$$\begin{eqnarray}\ell \leq \frac{u^{T}Au}{u^{T}\widehat{A}u}\leq {\rm\Upsilon}\;\;\text{and}\;\;{\it\gamma}\leq \frac{p^{T}Sp}{p^{T}\widehat{S}p}\leq {\rm\Gamma}\;\;\text{for all}\;\,u\in \mathbb{R}^{n}-\{0\}\;\,\text{and}\;\,p\in \mathbb{R}^{m}-\{0\}\end{eqnarray}$$

for positive constants $\ell ,{\rm\Upsilon},{\it\gamma}$ and ${\rm\Gamma}$ , then there are positive constants $a,b,c$ and $d$ satisfying the preconditioned minres convergence bound (5.3). In the context of PDE problems, for example, this means that preconditioned minres is a solver of optimal computational complexity provided that $\widehat{A}$ is spectrally equivalent to $A$ and $\widehat{S}$ is spectrally equivalent to $S$ . A multigrid cycle can be an excellent choice for $\widehat{A}$ . See Pestana and Wathen (Reference Pestana and Wathen2015b ) for a general proof which highlights the role of inf-sup stability associated with the minimal singular value of $B$ (in appropriate norms).

For many problems there are natural Schur complement approximations which arise from problem formulation and do not require separate approximation of the matrices $B$ and $A$ ; see Pestana and Wathen (Reference Pestana and Wathen2015b ). In the case of the Stokes equations, the natural Schur complement approximation is particularly simple: it is the mass matrix (3.4) (Verfürth Reference Verfürth1984, Silvester and Wathen Reference Silvester and Wathen1994). For any such problem with a natural Schur complement approximation, (5.7) is an excellent preconditioning approach since $\widehat{A}$ can usually be any of the approximations in Section 3. It should be noted that once such approximations $\widehat{A}$ and $\widehat{S}$ are identified, then there are more involved ways to make use of them than the simple block diagonal preconditioner (5.7). When matrix symmetry is not an issue, a block triangular preconditioner can be employed; see Section 6.

Sometimes it is actually useful to directly approximate the matrices $B$ and $A$ in constructing $\widehat{S}$ , for example in the context of optimization with constraints which are PDEs (Rees, Dollar and Wathen Reference Rees, Dollar and Wathen2010a ), though care must be taken to satisfy the sufficient conditions of Braess and Peisker (Reference Braess and Peisker1986).

For saddle point problems arising particularly in the more general context of optimization, constraint preconditioning has attracted significant attention (Keller, Gould and Wathen Reference Keller, Gould and Wathen2000). Here

$$\begin{eqnarray}P=\left[\begin{array}{@{}cc@{}}H & B^{T}\\ B & 0\end{array}\right]\end{eqnarray}$$

is used as a preconditioner for ${\mathcal{A}}$ . Note that the blocks $B$ and $B^{T}$ which come from the constraints are kept, whilst $A$ is replaced by another symmetric and (usually) positive definite matrix $H$ . It might appear that solution of systems with $P$ may not be significantly easier than with ${\mathcal{A}}$ , but by leaving $H$ implicitly defined, this can be achieved (Dollar and Wathen Reference Dollar and Wathen2006, Dollar, Gould, Schilders and Wathen Reference Dollar, Gould, Schilders and Wathen2006). One significant advantage of keeping the constraint blocks in the preconditioner is that the preconditioned problem is now equivalent to a positive definite problem on the constraint manifold, hence a projected cg algorithm can be used (Gould, Hribar and Nocedal Reference Gould, Hribar and Nocedal2001); this is rather useful since a constraint preconditioner is symmetric and indefinite, precluding simple use of minres. Keller et al. (Reference Keller, Gould and Wathen2000) prove that when $H$ is positive definite then $P^{-1}{\mathcal{A}}$ has $2m$ eigenvalues of $1$ , each of geometric multiplicity $2$ (thus the minimum polynomial generically has a factor $({\it\lambda}-1)^{2}$ ), and that the remaining $n-m$ eigenvalues are real and interlace those of $H^{-1}A$ . Thus if $H$ is a good preconditioner for $A$ and systems with $P$ can be readily solved, then this can be an attractive strategy.

Different challenges arise and different strategies are needed for saddle point problems (5.5) with semidefinite block $A$ . The constraint preconditioning paradigm can be applied, in which case finding a matrix $H$ for the semidefinite matrix $A$ may require ideas as in Section 4 above. The matrix ${\mathcal{A}}$ is still invertible provided $A$ is positive definite on the kernel of $B$ and $B$ is of full rank, as has been mentioned. One possible approach is therefore via an augmentation or augmented system preconditioner

$$\begin{eqnarray}P=\left[\begin{array}{@{}cc@{}}A+B^{T}W^{-1}B & 0\\ 0 & W\end{array}\right]\end{eqnarray}$$

for some appropriate symmetric and positive definite matrix $W$ . This approach is described and analysed by Greif and Schötzau (Reference Greif and Schötzau2006), who were motivated by the solution of Maxwell equations describing electromagnetic phenomena – a key PDE problem with this structure. For related but earlier ideas associated with the semidefiniteness of the $\text{div}$ operator, see Arnold, Falk and Winther (Reference Arnold, Falk and Winther1997).

For the various Maxwell equation formulations, semidefiniteness arises through $\text{curl}~\text{curl}$ terms, indefiniteness because of a $\text{div}$ constraint, and wave-like solutions are always expected, which means that fine computational grids must be used, especially for high-frequency problems. Fortunately, Maxwell’s equations are linear and self-adjoint in most formulations, which is why we bring them up in this section. Mixed finite element approximation is widely used for these PDEs. Representative of the large and still growing literature on preconditioning for such problems are Greif and Schötzau (Reference Greif and Schötzau2007), who develop block preconditioners, and Tsuji, Engquist and Ying (Reference Tsuji, Engquist and Ying2012), who extend sweeping preconditioner ideas to these problems.

A slightly more detailed review of preconditioning for saddle point problems can be found in Benzi and Wathen (Reference Benzi, Wathen, Schilders, van der Vorst and Rommes2009).

5.3. A word of caution

As has already been mentioned in Section 2, for any non-singular square matrix $A$ , the solution of a linear system can be achieved by using an indefinite symmetric system of twice the dimension of the form (2.5) or even

$$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}0 & A\\ A^{T} & 0\end{array}\right]\left[\begin{array}{@{}c@{}}r\\ x\end{array}\right]=\left[\begin{array}{@{}c@{}}b\\ 0\end{array}\right]\!.\end{eqnarray}$$

Since the eigenvalues of the coefficient matrix here are exactly symmetric about the origin (they are $\pm {\it\sigma}_{i}(A)$ , where ${\it\sigma}_{i}(A)$ are the singular values of $A$ ), generically, the convergence of minres would stagnate at every other iteration, so the computational effort in solution is exactly the same as for the normal equations (2.4) (Fischer, Ramage, Silvester and Wathen Reference Fischer, Ramage, Silvester and Wathen1998). For (2.5) a similar issue arises.

The artificial construction of indefinite symmetric systems is therefore unlikely to achieve other than equivalent problem statements, except possibly in the solution of least-squares problems by augmented matrix methods (Björck Reference Björck1996, Chapter 7). To go further, a symmetric indefinite matrix of any form which has eigenvalues symmetrically arranged about the origin will generically lead to staircasing behaviour, where progress is made only on alternate minres iterations; cg for (2.4) generally requires the same computational effort for equation solution in such circumstances. In this sense, it is highly unfortunate that the tractable minres convergence bound (5.2) applies only when the positive and negative eigenvalues extend over intervals of equal length. If in fact $a=d$ and $b=c$ , so that these intervals are centred on the origin, then the undesirable staircasing effect is always implied! In general, it is advantageous in terms of convergence speed if there is some asymmetry in the distribution of positive and negative eigenvalues of preconditioned indefinite symmetric matrices. However, it is then considerably more complicated to take account of this asymmetry in a minres convergence bound (but see Wathen, Fischer and Silvester Reference Wathen, Fischer and Silvester1995).

6.1. Incomplete factorization

Just as a Cholesky factorization for a symmetric and positive definite matrix can be incompletely formed, an LU factorization (equivalent to Gauss elimination) for a non-symmetric matrix can be performed incompletely by prespecifying allowable fill or thresholding or both. Perhaps such ILU preconditioning methods and their variants are the most widely used for non-symmetric problems, and software is widely available.

Because the product of triangular matrix factors is computed for ILU, for problems with an implied directionality such as convection–diffusion, ordering of the variables can lead to one of the factors being less important (and perhaps more nearly diagonal).

Although, in general, there is little we can say about ILU preconditioning for non-symmetric matrices, gmres with ILU is a ‘go-to’ technique for very many practical problems, and therefore much largely empirical research continues to explore every aspect of the various possibilities, for example, multilevel variants, block variants, and parallel implementation issues. We again refer to Benzi (Reference Benzi2002), as the huge literature is fairly evenly cited. Because of its significance, particularly in the important area of petroleum reservoir modelling, we specifically mention the nested ILU methods originally derived by Appleyard and Cheshire (Reference Appleyard and Cheshire1983), which use partial eliminations and are designed to take advantage of parallel computation; related ideas were developed, for example, by Saad (Reference Saad1996). See also Vassilevski (Reference Vassilevski2008) for a comprehensive treatment of these methods in a more general framework.

6.2. Multigrid

There is no inherent obstacle to applying multigrid or multilevel ideas for non-symmetric problems. However, smoothing is generally not so easy and requires more work than for symmetric problems, and care must be taken so that coarse grid correction does not reintroduce high frequencies. Without attention to these aspects, the methods simply do not work at all well. Since such issues vary with differing applications, we briefly discuss only the important problem associated with the convection–diffusion equation.

Smoothing – or perhaps more appropriately in this setting, ‘sweeping’ – should ideally be something like a Gauss–Seidel iteration with the variables ordered in the direction of the convection. For problems with recirculation, one can try multidirectional sweeping, so that all parts of the flow have at least one of the directional sweeps that is approximately in the convection direction. Coarse grid correction then requires care (in grid transfers and/or coarse grid operator). In the convection–diffusion geometric multigrid method due to Ramage (Reference Ramage1999), for example, the coarse grid operator is recomputed on each coarse grid with – crucially – streamline upwinding appropriate to that grid size. The point is that when a grid under-resolves any layer-like phenomenon which arises with convection – this must happen at some level if one has a true multigrid concept – then oscillatory solutions are commonly the result. The reintroduction of such high frequencies would require more significant sweeping or smoothing, and it is preferable to provide stabilization appropriate to each level of grid so that smooth coarse grid corrections result.

Multigrid for hyperbolic problems is thought about rather differently than for elliptic problems, and the idea of directional sweeping, rather than smoothing, is key. The coarse grids allow one to sweep errors out of the boundary of the domain more rapidly than using just a single grid. Again, transport is the main issue, not smoothing.

Algebraic multigrid ideas have been tried and can work well for some problems. The computation of ‘coarse grids’ (or in reality, spaces of coarse variables) by purely algebraic means can identify layer phenomena, for example, and automatically semi-coarsen. In particular, widely used coarsening algorithms seem to detect anisotropy and lead to appropriate coarse spaces/grids. All is not rosy, however, and AMG software can quite simply fail. There is a particular danger that simply does not arise with geometric multigrid ideas, namely that the coarse grid can have a significant proportion of the variables of the fine grid; that is, coarsening does not lead to a sufficiently significant reduction in the number of variables. If this occurs, then many levels are required and efficiency can be severely compromised. Active research continues on AMG for non-symmetric problems; see for example Notay (Reference Notay2012).

Multigrid methodology is generally applied as a preconditioner for non-symmetric problems.

6.3. Sparse approximate inverses

Norm minimization readily leads to sparse approximate inverses similarly to the symmetric case. A biconjugation based method (AINV) and its more robust stabilized form (SAINV) are, however, perhaps the most widely considered of such methods in the non-symmetric case (Benzi and Tuma Reference Benzi and Tuma1998, Benzi et al. Reference Benzi, Cullum and Tuma2000).

6.4. Symmetric/skew splitting

Every non-symmetric real matrix can be uniquely split into its symmetric and skew-symmetric parts:

$$\begin{eqnarray}A=\underbrace{{\textstyle \frac{1}{2}}(A+A^{T})}_{H}+\underbrace{{\textstyle \frac{1}{2}}(A-A^{T})}_{S}.\end{eqnarray}$$

(The $H$ stands for Hermitian, since the same idea obviously applies in the complex case when transpose is replaced by conjugate transpose.) When $H$ is positive definite, preconditioning of $A$ by $H$ gives a shifted skew-symmetric matrix for which a special three-term recurrence Krylov subspace method exists (Concus and Golub Reference Concus and Golub1976, Widlund Reference Widlund1978), and convergence depends on the size of $1+|{\it\lambda}_{\max }(H^{-1}S)|$ (Widlund Reference Widlund1978, Freund and Ruscheweyh Reference Freund and Ruscheweyh1986). Although this approach has been considered for some applications – see for example Elman and Schultz (Reference Elman and Schultz1986) in the context of convection–diffusion problems – the so-called HSS (Hermitian skew-Hermitian splitting) methods originally introduced by Bai, Golub and Ng (Reference Bai, Golub and Ng2003) have been rather more widely investigated. They proved the convergence of an alternate iteration using shifted Hermitian and shifted skew-Hermitian parts of the form

$$\begin{eqnarray}\displaystyle ({\it\alpha}I+H)x^{(k+1/2)} & = & \displaystyle ({\it\alpha}I-S)x^{(k)}+b,\nonumber\\ \displaystyle ({\it\alpha}I+S)x^{(k+1)} & = & \displaystyle ({\it\alpha}I-H)x^{(k+1/2)}+b\nonumber\end{eqnarray}$$

for $Ax=b$ . Though convergence is guaranteed independently of the choice of the parameter ${\it\alpha}>0$ , it can be slow. Nevertheless, many variants of using this idea as a preconditioner have been suggested, as have related structural splitting ideas, and good results are reported for some; see for example Benzi, Ng, Niu and Wang (Reference Benzi, Ng, Niu and Wang2011a ).

6.5. Non-symmetric saddle point systems

Saddle point problems must essentially be symmetric as in Section 5.2, but matrices of the structure (5.5) arise in which the block $A$ is non-symmetric. A key example is provided by linearizations of the incompressible Navier–Stokes equations; see Elman et al. (Reference Elman, Silvester and Wathen2014b , Chapters 8 and 9). It is these problems that we refer to as non-symmetric saddle point problems.

For a non-symmetric saddle point matrix, the results of Murphy et al. (Reference Murphy, Golub and Wathen2000) still hold, and block diagonal preconditioners can still be applied. However, since there is no symmetry to lose, the block triangular preconditioners

$$\begin{eqnarray}P=\left[\begin{array}{@{}cc@{}}\widehat{A} & 0\\ B & \pm \widehat{S}\end{array}\right]\quad \text{or}\quad P=\left[\begin{array}{@{}cc@{}}\widehat{A} & B^{T}\\ 0 & \pm \widehat{S}\end{array}\right]\end{eqnarray}$$

will give in the ‘ideal’ case $\widehat{A}=A,\widehat{S}=S=BA^{-1}B^{T}$ , a preconditioned system with either just the two distinct eigenvalues $\pm 1$ or even just the single eigenvalue of 1, depending on whether the $+$ or the $-$ sign is selected, respectively (Murphy et al. Reference Murphy, Golub and Wathen2000). Even with the single eigenvalue 1, $P$ is not the inverse, so something has to give! In fact it is diagonalizability that is lost in selecting the $-$ sign; the minimum polynomial is then $({\it\lambda}-1)^{2}$ , so two iterations would be required, exactly as for the $+$ case. However, in the practical case when the approximations $\widehat{A},\widehat{S}$ are not ideal, it is usually the $-$ sign which is chosen, so that all eigenvalues are clustered around 1, rather than the two clusters around $\pm 1$ that arise with the $+$ choice, though it is not clear that this is to be preferred (Ymbert, Embree and Sifuentes Reference Ymbert, Embree and Sifuentes2015).

Approximations of the Schur complement remain a key aspect of any such block preconditioning approach. For the common Oseen linearization of the incompressible Navier–Stokes problem, the matrix block $A$ comes from a vector convection–diffusion operator, so that an appropriate multigrid approach, as outlined above, or any other such good approximation of (preconditioner for) convection–diffusion, is required as $\widehat{A}$ . To identify an approximate Schur complement requires more in-depth consideration, but there are now two excellent candidates: the pressure convection–diffusion (PCD) (Kay, Loghin and Wathen Reference Kay, Loghin and Wathen2002) and the least-squares commutator (LSC) preconditioners. Both are derived, described and analysed in Elman et al. (Reference Elman, Silvester and Wathen2014b , Chapter 9). The PCD method is in some ways simpler, but it requires the construction of a convection–diffusion operator on the pressure space (the space of Lagrange multipliers), which is not required for the problem formulation itself. A purely algebraic way to compute this operator based on sparse approximate inverse technology was, however, described by Elman et al. (Reference Elman, Howle, Shadid, Shuttleworth and Tuminaro2006). The LSC method is derived via a commutator argument and is defined purely in terms of the matrix blocks which arise naturally for the Oseen problem. For a comparison and full consideration of advantages and disadvantages of these two approaches, see Elman et al. (Reference Elman, Silvester and Wathen2014b , Chapter 9).

The augmentation or augmented Lagrangian approach has also been applied to non-symmetric saddle point problems. Although some concerns have been expressed regarding possible distortion of the norm in which convergence occurs, good results have been reported even for large Reynolds number (Benzi, Olshanskii and Wang Reference Benzi, Olshanskii and Wang2011b ).

Earlier approaches for this problem involved segregated treatment of pressure and velocity and predate the preconditioned Krylov subspace iteration technology. It is perhaps worth commenting here that this example is one of the first for which block preconditioners comprising blocks designed for different sets of variables allow treatment of different ‘physics’ within the preconditioner, whilst gmres or some other Krylov subspace method is applied to the complete coupled problem. This is an important paradigm.

Such is its importance that it is in the context of non-symmetric saddle point systems, arising from approximation of the Oseen problem, that perhaps the most refined analysis of preconditioned gmres convergence has been pursued. There are several eigenvalue analyses for different preconditioners (see, e.g., Elman and Silvester Reference Elman and Silvester1996), although these do not lead to rigorous estimates as explained above. It is in the context of this problem, however, that field-of-values convergence analysis has been successfully employed to rigorously bound preconditioned gmres convergence and to elucidate its dependence on or independence of the mesh size, $h$ , and of the Reynolds number; see for example the analysis of Loghin and Wathen (Reference Loghin and Wathen2004) for PCD, and Benzi and Olshanskii (Reference Benzi and Olshanskii2011) for block triangular preconditioning of an augmented Lagrangian formulation.

6.6. Circulant and semi-circulant preconditioning

Our opening example already indicates the potential for circulant preconditioning, at least for symmetric Toeplitz matrices. In fact the applicability of such ideas is rather wider, though the theoretical justification is less complete for non-symmetric problems.

The analysis for symmetric and positive definite Toeplitz matrices is well described in the books by Chan and Jin (Reference Chan and Jin2007) and Ng (Reference Ng2004), who also indicate how similar ideas follow for block Toeplitz matrices with Toeplitz blocks. Indeed, the use of circulant preconditioners for non-symmetric Toeplitz matrices does not really depend on symmetry. However, unless the normal equations are used, convergence guarantees are not readily available, simply because of the general lack of descriptive convergence bounds for non-symmetric matrices. Good computational results have been observed in many situations, and in fact even proved for certain problems via bounding of the eigenvector condition number as well as the eigenvalues and using (6.1) (Otto Reference Otto1996). A recent observation of Pestana and Wathen (Reference Pestana and Wathen2015a ), that guaranteed convergence bounds can be achieved for minres, for circulant preconditioned non-symmetric Toeplitz matrices via a variable reordering, is also relevant.