solvers.rst

OOPS linear equations solvers

This note describes the linear equation solvers implemented in OOPS. The solvers are designed to be as generic as possible, to allow their use with a variety of different vector and matrix classes.

Generic code design

The solvers are implemented as generic (templated) functions. The functions are templated on a type VECTOR and one or more matrix types (AMATRIX, PMATRIX, etc.). Thus, we require that all vector arguments are of the same class, but the matrices do not have to be derived from a single class.

In addition to the vector and matrix arguments, there are two arguments specifying the maximum number of iterations to be performed, and the required reduction in residual norm. Iteration will stop when either the iteration limit is reached or the residual norm is reduced below the required factor. The return value from all the solvers is the achieved reduction in residual norm.

As a typical example, consider the IPCG solver:

template <typename VECTOR,
          typename AMATRIX,
          typename PMATRIX>

double IPCG (VECTOR & x,
             const VECTOR & b,
             const AMATRIX & A,
             const PMATRIX & precond,
             const int maxiter,
             const double tolerance );

For all the solvers, VECTOR is expected to implement basic linear algebra operations:

double dot_product(VECTOR &, VECTOR &);
operator(=)
operator(+=)
operator(-=)
operator(*=) // double times VECTOR
axpy // u.axpy(a,v) sets u = u + a*v

The matrix classes are expected to implement:

void apply(const VECTOR&, VECTOR&) const

This function represents application (multiplication) of the matrix to the first argument, with the result returned in the second. For all the current algorithms, application of the preconditioner may be approximate, for example as the result of an iterative solution method.

The algorithms

The following algorithms are available:

• IPCG: Inexact-Preconditioned Conjugate Gradients.
• DRIPCG: "Derber-Rosati" Inexact-Preconditioned Conjugate Gradients.
• GMRESR.
• DRGMRESR: A "Derber-Rosati" version of GMRESR.

IPCG

Inexact Preconditioned Conjugate Gradients GolubYe99 is a slight variation on the well-known Preconditioned Conjugate Gradients (PCG) algorithm. Given an initial vector x₀, a symmetric, positive definite matrix A, and preconditioners E_k ≈ A^− 1, IPCG solves Ax = b as follows:

$$\begin{aligned} \begin{eqnarray} \mathbf{s}_k &=& \mathbf{E}_k \mathbf{r}_k \\\ \beta_k &=& \frac{ \mathbf{s}_k^{\rm T} (\mathbf{r}_k -\mathbf{r}_{k-1})} { \mathbf{s}_k^{\rm T} \mathbf{r}_{k-1} } \qquad\mbox{for $k>0$} \\\ \mathbf{d}_k &=& \left\{ \begin{array}{ll} \mathbf{s}_k & \qquad\mbox{if $k=0$} \\\ \mathbf{s}_k + \beta_k \mathbf{d}_{k-1} & \qquad\mbox{if $k>0$} \end{array} \right. \\\ \mathbf{w}_k &=& \mathbf{A} \mathbf{d}_k \\\ \alpha_k &=& (\mathbf{s}_k^{\rm T} \mathbf{r}_k )/(\mathbf{d}_k^{\rm T} \mathbf{w}_k )\\\ \mathbf{x}_{k+1} &=& \mathbf{x}_k + \alpha_k \mathbf{d}_k \\\ \mathbf{r}_{k+1} &=& \mathbf{r}_k - \alpha_k \mathbf{w}_k \end{eqnarray} \end{aligned}$$

The algorithm differs from PCG only in the definition of β_k, which PCG defines as:

$$\begin{aligned} \begin{equation} \beta_k = \frac{\mathbf{s}_k^{\rm T} \mathbf{r}_k } {\mathbf{s}_k^{\rm T} \mathbf{r}_{k-1} } \qquad\mbox{for $k>0$}. \\\ \end{equation} \end{aligned}$$

This slight modification requires additional storage for the vector r_k − 1, but has the advantage of significantly improving the convergence properties of the algorithm in the case that the preconditioner varies from iteration to iteration. In particular, s_k in equation IPCG-precond-eqn can be determined as the result of a truncated iterative solution of the equation Ãs_k = r_k, for some approximation Ã ≈ A.

Convergence results for IPCG are presented by GolubYe99. Further results are given by Knyazev08.

DRIPCG

In many applications in variational data assimilation, the matrix A takes the particular form:

A = B^− 1 + C

Furthermore, the matrix B^− 1 may be ill-conditioned or unavailable. In this case, DerberRosati89 showed that the PCG algorithm could be modified in such a way that application of B^− 1 is not required during the iterations. The algorithm requires an additional initial vector, x̂₀ = B^− 1x₀. However application of B^− 1 can be avoided for the initial vector if the initial guess is taken as x₀ = 0, or if both x₀ and x̂₀ are retained from a previous application of the algorithm.

DerberRosati89 modified PCG. However, their approach can be applied more widely to a range of linear equation solvers. The essence of the approach is to introduce auxiliary vectors:

$$\begin{aligned} \mathbf{\hat x}_k &=& \mathbf{B}^{-1} \mathbf{x}_k \\\ \mathbf{\hat s}_k &=& \mathbf{B}^{-1} \mathbf{s}_k \\\ \mathbf{\hat d}_k &=& \mathbf{B}^{-1} \mathbf{d}_k \end{aligned}$$

Defining F_k = B^− 1E_k, and r₀ = b − x̂₀ − Cx₀ (i.e. r₀ = b − Ax₀), we can write the IPCG algorithm as:

$$\begin{aligned} \begin{eqnarray} \mathbf{\hat s}_k &=& \mathbf{F}_k \mathbf{r}_k \\\ \mathbf{s}_k &=& \mathbf{B} \mathbf{r}_k \\\ \beta_k &=& \frac{ \mathbf{s}_k^{\rm T} (\mathbf{r}_k -\mathbf{r}_{k-1})} { \mathbf{s}_k^{\rm T} \mathbf{r}_{k-1} } \qquad\mbox{for $k>0$} \\\ \mathbf{d}_k &=& \left\{ \begin{array}{ll} \mathbf{s}_k & \qquad\mbox{if $k=0$} \\\ \mathbf{s}_k + \beta_k \mathbf{d}_{k-1} & \qquad\mbox{if $k>0$} \end{array} \right. \\\ \mathbf{\hat d}_k &=& \left\{ \begin{array}{ll} \mathbf{\hat s}_k & \qquad\mbox{if $k=0$} \\\ \mathbf{\hat s}_k + \beta_k \mathbf{\hat d}_{k-1} & \qquad\mbox{if $k>0$} \end{array} \right. \\\ \mathbf{w}_k &=& \mathbf{\hat d}_k + \mathbf{C} \mathbf{d}_k \\\ \alpha_k &=& (\mathbf{s}_k^{\rm T} \mathbf{r}_k )/(\mathbf{d}_k^{\rm T} \mathbf{w}_k )\\\ \mathbf{x}_{k+1} &=& \mathbf{x}_k + \alpha_k \mathbf{d}_k \\\ \mathbf{\hat x}_{k+1} &=& \mathbf{\hat x}_k + \alpha_k \mathbf{\hat d}_k \\\ \mathbf{r}_{k+1} &=& \mathbf{r}_k - \alpha_k \mathbf{w}_k \end{eqnarray} \end{aligned}$$

Note that no applications of B^− 1 are required during the iteration. Note also that x_k is not used during the iteration, so that the equation for x_k + 1 can be removed. After some number N of iterations, we can recover x_N from x̂_N by multiplying the latter by B.

The Derber Rosati algorithm is sometimes called "Double" PCG. We have adopted this nomenclature for algorithms that include similar modifications. thus, we call the algorithm described above Derber-Rosati Inexact-Preconditioned Conjugate Gradients, or DRIPCG. The algorithm is closely related to CGMOD (Gratton, personal communication).

DRIPCG is algebraically equivalent to IPCG provided that the preconditioners are related by F_k = B^− 1E_k. A common preconditioning is to choose E_k = B, in which case F_k = I.

GMRESR

GMRESR VanDerVorst94 is a robust algorithm for square, non-symmetric systems. Like IPCG, it allows the preconditioner to vary from iteration to iteration. The algorithm starts with r₀ = b − Ax₀, and iterates the following steps for k = 0, 1, 2, ….

$$\begin{aligned} \begin{eqnarray} \mathbf{z} &=& \mathbf{E}_k \mathbf{r}_k \\\ \mathbf{c} &=& \mathbf{A}\mathbf{z} \\\ && \mbox{for} \quad j = 0,1,\ldots,k-1 \nonumber \\\ && \qquad \alpha = \mathbf{c}_j^{\rm T} \mathbf{c} \\\ && \qquad \mathbf{c} = \mathbf{c} - \alpha \mathbf{c}_j \\\ && \qquad \mathbf{z} = \mathbf{z} - \alpha \mathbf{u}_j \\\ && \mbox{end for} \nonumber \\\ \mathbf{c}_k &=& \frac{\mathbf{c}}{\Vert \mathbf{c} \Vert_2} \\\ \mathbf{u}_k &=& \frac{\mathbf{z}}{\Vert \mathbf{c} \Vert_2} \\\ \beta_k &=& \mathbf{c}_k^{\rm T} \mathbf{r}_k \\\ \mathbf{x}_{k+1} &=& \mathbf{x}_k + \beta_k \mathbf{u}_k \\\ \mathbf{r}_{k+1} &=& \mathbf{r}_k - \beta_k \mathbf{c}_k \end{eqnarray} \end{aligned}$$

For a symmetric matrix and constant symmetric positive definite (SPD) preconditioner, GMRESR is algebraically equivalent to PCG. In this case, the explicit orthogonalization of c_k against earlier vectors mitigates the effects of rounding error, resulting in somewhat faster convergence and a preservation of the super-linear convergence properties of PCG.

The storage requirements of GMRESR are significant, since the vectors c_k and u_k must be retained for all subsequent iterations. Note that this is twice the storage required for a fully-orthogonalizing PCG algorithm such as CONGRAD Fisher98.

DRGMRESR

A "Derber-Rosati" version of GMRESR is easy to derive. As in the case of DRIPCG, we define F_k = B^− 1E_k, and calculate the starting point as r₀ = b − x̂₀ − Cx₀, where x̂₀ = B^− 1x₀. Defining also the auxilliary vectors:

$$\begin{aligned} \mathbf{\hat z} &=& \mathbf{B}^{-1} \mathbf{z} \\\ \mathbf{\hat u}_k &=& \mathbf{B}^{-1} \mathbf{u}_k \\\ \end{aligned}$$

we have:

$$\begin{aligned} \begin{eqnarray} \mathbf{\hat z} &=& \mathbf{F}_k \mathbf{r}_k \\\ \mathbf{z} &=& \mathbf{B} \mathbf{\hat z}_k \\\ \mathbf{c} &=& \mathbf{\hat z} + \mathbf{C}\mathbf{z} \\\ && \mbox{for} \quad j = 0,1,\ldots,k-1 \nonumber \\\ && \qquad \alpha = \mathbf{c}_j^{\rm T} \mathbf{c} \\\ && \qquad \mathbf{c} = \mathbf{c} - \alpha \mathbf{c}_j \\\ && \qquad \mathbf{z} = \mathbf{z} - \alpha \mathbf{u}_j \\\ && \qquad \mathbf{\hat z} = \mathbf{\hat z} - \alpha \mathbf{\hat u}_j \\\ && \mbox{end for} \nonumber \\\ \mathbf{c}_k &=& \frac{\mathbf{c}}{\Vert \mathbf{c} \Vert_2} \\\ \mathbf{u}_k &=& \frac{\mathbf{z}}{\Vert \mathbf{c} \Vert_2} \\\ \mathbf{\hat u}_k &=& \frac{\mathbf{\hat z}}{\Vert \mathbf{c} \Vert_2} \\\ \beta_k &=& \mathbf{c}_k^{\rm T} \mathbf{r}_k \\\ \mathbf{\hat x}_{k+1} &=& \mathbf{\hat x}_k + \beta_k \mathbf{\hat u}_k \\\ \mathbf{r}_{k+1} &=& \mathbf{r}_k - \beta_k \mathbf{c}_k \end{eqnarray} \end{aligned}$$

As in the case of DRIPCG, after N iterations, we can recover the solution x_N from x̂_N by multiplying the latter by B.