# Iterative solvers and preconditioning

A successful family of methods, usually referred to as Conjugate
Gradient-like algorithms, or Krylov subspace methods, can be viewed as
Galerkin or least-squares methods applied to a linear system $Ax =b$.
This view is different from the standard approaches to deriving the
classical Conjugate Gradient method in the literature. Nevertheless,
the fundamental ideas of least squares and Galerkin approximations
can be used to derive the most
popular and successful methods for linear systems, and this is the
topic of the present chapter.  Such a view may increase the general
understanding of variational methods and their applicability.


Given a linear system

<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>

$$
\begin{equation}
Ax = b,\quad x,b\in\mathbb{R}^n,\ A\in\mathbb{R}^{n,n}
\label{_auto1} \tag{1}
\end{equation}
$$

and a start vector $x^0$,
we want to construct an iterative solution method that produces approximations
$x^1,x^2,\ldots$, which hopefully converge to the exact solution $x$.
In iteration no. $k$ we seek an approximation

<!-- Equation labels as ordinary links -->
<div id="linalg:xk:update"></div>

$$
\begin{equation}
x^{k+1} = x^{k} + u,\quad u = \sum_{j=0}^k c_j{q}_j
\label{linalg:xk:update} \tag{2}
\end{equation}
$$

where ${q}_j\in\mathbb{R}^n$ are known vectors
<!-- spanning the space $V_{k+1}$, -->
and $c_j$ are constants to be determined.
To be specific, let $q_0,\ldots,q_k$ be basis vectors for $V_{k+1}$:

$$
V_{k+1} = \hbox{span}\{q_0,\ldots,q_k\}{\thinspace .}
$$

The associated inner product $(\cdot,\cdot )$ is here
the standard Euclidean inner product on $\mathbb{R}^n$.

The corresponding error in the equation $Ax =b$, the residual, becomes

$$
r^{k+1} = b -Ax^{k+1} =
r^{k} -\sum_{j=0}^kc_jA{q}_j {\thinspace .}
$$

# Conjugate gradient-like iterative methods
<div id="ch:linalg:CGmethods"></div>


## The Galerkin method

Galerkin's method states that the error in the equation, the residual,
is orthogonal to the space $V_{k+1}$ where we seek the approximation to
the problem.

The Galerkin (or projection) method aims at finding $u=\sum_jc_jq_j\in V_{k+1}$
such that

$$
(r^{k+1},v)=0,\quad\forall v\in V_{k+1}{\thinspace .}
$$

This statement is
equivalent to the residual being orthogonal to each basis vector:

<!-- Equation labels as ordinary links -->
<div id="linalg2:cg:eq1"></div>

$$
\begin{equation}
(r^{k+1},{q}_i )=0,\quad i=0,\ldots,k{\thinspace .}
\label{linalg2:cg:eq1} \tag{3}
\end{equation}
$$

Inserting the expression for $r^{k+1}$ in ([3](#linalg2:cg:eq1))
gives a linear system
for $c_j$:

<!-- Equation labels as ordinary links -->
<div id="_auto2"></div>

$$
\begin{equation}
\sum_{j=0}^k (A{q}_i,{q}_j )c_j = (r^{k},{q}_i),\quad i=0,\ldots,k
{\thinspace .}
\label{_auto2} \tag{4}
\end{equation}
$$

## The least squares method

The idea of the least-squares method is to
minimize the square of the norm of the
residual with respect to the free parameters $c_0,\ldots,c_k$.
That is, we minimize $(r^{k+1},r^{k+1})$:

$$
{\partial\over\partial c_i} (r^{k+1},r^{k+1}) =
2({\partial r^{k+1}\over
\partial c_i},r^{k+1}) =0,\quad i=0,\ldots,k{\thinspace .}
$$

Since $\partial r^{k+1} /\partial c_i = -A{q}_i$, this approach leads to
the following linear system:

<!-- Equation labels as ordinary links -->
<div id="ch:linalg:ls:eq1"></div>

$$
\begin{equation}
\sum_{j=0}^k (A{q}_i,A{q}_j)c_j = (r^{k+1},A{q}_i),\quad i=0,\ldots,k{\thinspace .}
\label{ch:linalg:ls:eq1} \tag{5}
\end{equation}
$$

## Krylov subspaces


To obtain a complete algorithm, we need to establish a rule to update
the basis $\mathcal{B}=\{q_0,\ldots,{q}_k\}$ for the next iteration.
That is, we need to
compute a new basis vector ${q}_{k+1}\in V_{k+2}$ such that

<!-- Equation labels as ordinary links -->
<div id="linalg2:cg:Basis"></div>

$$
\begin{equation}
\mathcal{B} =\{ {q}_0,\ldots,{q}_{k+1}\}
\label{linalg2:cg:Basis} \tag{6}
\end{equation}
$$

is a basis for the space $V_{k+2}$ that is used in the next iteration.
The present family of methods applies the *Krylov space*, where $V_k$
is defined as

<!-- Equation labels as ordinary links -->
<div id="_auto3"></div>

$$
\begin{equation}
V_{k} = \mbox{span} \left\lbrace r^0,Ar^0,A^2r^0,\ldots
A^{k-1}r^0 \right\rbrace {\thinspace .}
\label{_auto3} \tag{7}
\end{equation}
$$

Some frequent names
of the associated iterative methods are therefore
{Krylov subspace iterations}, Krylov projection methods,
or simply Krylov methods.
It is a fact that $V_{k} \subset V_{k+1}$ and that $r^0
Ar^0,\ldots A^{k-1}r^0$ are linearly independent
vectors.

## Computation of the basis vectors

A potential formula for updating ${q}_{k+1}$, such that
${q}_{k+1}\in V_{k+2}$, is

<!-- Equation labels as ordinary links -->
<div id="linalg:q:update1"></div>

$$
\begin{equation}
{q}_{k+1} =  r^{k+1} + \sum_{j=0}^k\beta_j{q}_j{\thinspace .}
\label{linalg:q:update1} \tag{8}
\end{equation}
$$

(Since $r^{k+1}$ involves $Aq_k$, and $q_k\in V_{k+1}$, multiplying
by $A$ raises the dimension of the Krylov space by 1, so $Aq_k\in V_{k+2}$.)
The free parameters $\beta_j$ can be used to enforce desirable
orthogonality properties of ${q}_0,\ldots,{q}_{k+1}$. For example,
it is convenient to require that the coefficient matrices in the linear
systems for $c_0,\ldots,c_k$ are diagonal.
Otherwise, we must solve a $(k+1)\times (k+1)$ linear system in each iteration.
If $k$ should approach
$n$, the systems for the coefficients $c_i$ are of
the same size as our original system $Ax =b$!
A diagonal matrix, however, ensures an efficient closed form solution for
$c_0,\ldots,c_k$.


To obtain a diagonal coefficient matrix, we require in Galerkin's method
that

$$
(A{q}_i,{q}_j)=0\quad \mbox{when}\ i\neq j,
$$

whereas we in the least-squares method require

$$
(A{q}_i,A{q}_j)=0\quad \mbox{when}\ i\neq j{\thinspace .}
$$

We can define
the inner product

<!-- Equation labels as ordinary links -->
<div id="_auto4"></div>

$$
\begin{equation}
\langle u,v\rangle\equiv (Au,v ) =u^TAv,
\label{_auto4} \tag{9}
\end{equation}
$$

provided $A$ is symmetric and positive definite. Another
useful inner product is

<!-- Equation labels as ordinary links -->
<div id="_auto5"></div>

$$
\begin{equation}
[u,v ]\equiv (Au,Av) = u^TA^TAv {\thinspace .}
\label{_auto5} \tag{10}
\end{equation}
$$

These inner products will be referred to as the $A$ product, with the
associated $A$ norm, and the $A^TA$ product, with the associated
$A^TA$ norm.

The orthogonality condition on the ${q}_i$ vectors are then
$\langle {q}_{k+1},{q}_i\rangle =0$ in the Galerkin method
and  $[{q}_{k+1},{q}_i]=0$ in the least-squares method, where $i$
runs from 0 to $k$.
A standard Gram-Schmidt process can be used for constructing
${q}_{k+1}$ orthogonal to ${q}_0,\ldots,{q}_k$. This leads to the determination
of the $\beta_0,\ldots,\beta_k$ constants as

<!-- Equation labels as ordinary links -->
<div id="linsys:betaiG"></div>

$$
\begin{equation}
\beta_i = { \langle r^{k+1},{q}_i\rangle\over\langle{q}_i,{q}_i\rangle}
\quad\hbox{(Galerkin)} \label{linsys:betaiG} \tag{11}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto6"></div>

$$
\begin{equation} 
\beta_i = { [r^{k+1},{q}_i]\over [{q}_i,{q}_i] }
\quad\hbox{(least squares)}
\label{_auto6} \tag{12}
\end{equation}
$$

for $i=0,\ldots,k$.

## Computation of a new solution vector

The orthogonality condition on the basis vectors ${q}_i$ leads to
the following solution for $c_0,\ldots,c_k$:

<!-- Equation labels as ordinary links -->
<div id="linsys:alpha:G"></div>

$$
\begin{equation}
c_i = { (r^{k},{q}_i)\over \langle {q}_i,{q}_i\rangle}
\quad\hbox{(Galerkin)} \label{linsys:alpha:G} \tag{13}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="linsys:alpha:LS"></div>

$$
\begin{equation} 
c_i = { (r^{k},A{q}_i)\over [ {q}_i,{q}_i]}
\quad\hbox{(least squares)} \label{linsys:alpha:LS} \tag{14}
\end{equation}
$$

In iteration $k$, $(r^{k},{q}_i)=0$ and
$(r^{k},A{q}_i)=0$, for $i=0,\ldots,k-1$, in the Galerkin and
least squares case, respectively.  Hence, $c_i =0$, for $i=0,\ldots,
k-1$. In other words,

$$
x^{k+1} = x^{k} + c_k{q}_k {\thinspace .}
$$

When $A$ is symmetric and positive definite, one can show that also
$\beta_i=0$, for $0=1,\ldots,k-1$, in both the Galerkin and least
squares methods.  This means that $x^k$ and
${q}_{k+1}$ can be updated using only ${q}_k$ and not the previous
${q}_0,\ldots,{q}_{k-1}$ vectors.  This property has of course dramatic
effects on the storage requirements of the algorithms as the number of
iterations increases.

For the suggested algorithms to work, we must require that the
denominators in ([13](#linsys:alpha:G)) and ([14](#linsys:alpha:LS)) do
not vanish.  This is always fulfilled for the least-squares method,
while a (positive or negative) definite matrix $A$ avoids break-down
of the Galerkin-based iteration (provided ${q}_i \neq 0$).

The Galerkin solution method for linear systems was originally devised
as a *direct method* in the 1950s.  After $n$ iterations the exact
solution is found in exact arithmetic, but at a higher cost than when
using Gaussian elimination. Naturally, the method did not receive
significant popularity before researchers discovered (in the beginning
of the 1970s) that the method could produce a good approximation to
$x$ for $k\ll n$ iterations.  The method, called the Conjugate
Gradient method, has from then on caught considerable interest as an
iterative scheme for solving linear systems arising from PDEs
discretized such that the coefficient matrix becomes sparse.

Finally, we mention how to terminate the iteration.
The simplest criterion is $||r^{k+1}||\leq\epsilon_r$, where
$\epsilon_r$ is a small prescribed tolerance.
Sometimes it is more appropriate to use a relative residual,
$||r^{k+1}||/||r^0||\leq\epsilon_r$.
Termination criteria for Conjugate Gradient-like methods is a subject
on its own.

## Summary of the least squares method

In the algorithm below, we have summarized
the computational steps in the least-squares method.
Notice that we update the residual recursively instead of using
$r^k=b - Ax^k$ in each iteration since we then avoid a possibly
expensive matrix-vector product.

1. given a start vector $x^0$,
   compute $r^0=b - Ax^0$ and set ${q}_0 =r^0$.

2. for $k=0,1,2,\ldots$ until termination criteria are fulfilled:

a. $c_k = {(r^{k},A{q}_k)/ [{q}_k,{q}_k]}$

b. $x^{k+1} = x^{k} + c_k{q}_k$

c. $r^{k+1} = r^{k} - c_kA{q}_k$

d. if $A$ is symmetric then

1. $\beta_k = {[r^{k+1},{q}_k]/ [{q}_k,{q}_k]}$

a. ${q}_{k+1} = r^{k+1} - \beta_k{q}_k$



b. else

1. $\beta_j = {[r^{k+1},{q}_j]/ [{q}_j,{q}_j]},\quad j=0,\ldots,k$

2. ${q}_{k+1} = r^{k+1} - \sum_{j=0}^k\beta_j{q}_j$



### Remark

The algorithm above is just a summary of the steps in the derivation
of the least-squares method and should not be directly used for
practical computations without further developments.

## Truncation and restart

When $A$ is nonsymmetric, the storage requirements of ${q}_0,\ldots,{q}_k$
may be prohibitively large. It has become a standard trick to
either *truncate* or *restart* the algorithm.
In the latter case one restarts the algorithm every $K+1$-th step, i.e.,
one aborts the iteration and starts the algorithm again with $x^0=x^K$.
The other alternative is to truncate the sum $\sum_{j=0}^k\beta_j{q}_j$
and use only the last $K$ vectors:

$$
x^{k+1} = x^{k} + \sum_{j=k-K}^k \beta_j{q}_j{\thinspace .}
$$

Both the restarted and truncated version of the algorithm require
storage of only $K+1$ basis vectors ${q}_{k-K},\ldots,{q}_k$.
The basis vectors are also often called *search direction vectors*.
The truncated version of the least-squares algorithm above
is widely known as *Generalized Minimum Residuals*, shortened as GMRES,
or GMRES$(K)$ to explicitly indicate the
number of search direction vectors.
In the literature one encounters the name
*Generalized Conjugate Residual method*, abbreviated CGR, for
the restarted version of Orthomin. When $A$ is symmetric, the method
is known under the name *Conjugate Residuals*.


## Summary of the Galerkin method


In case of Galerkin's method, we assume that $A$ is symmetric and
positive definite. The resulting computational procedure
is the famous Conjugate Gradient
method.
Since $A$ must be symmetric, the recursive update of
${q}_{k+1}$ needs only one previous search direction vector ${q}_k$, that is,
$\beta_j=0$ for $j<k$.

1. Given a start vector $x^0$,
   compute $r^0=b - Ax^0$ and set ${q}_0 =r^0$.

2. for $k=1,2,\ldots$ until termination criteria are fulfilled:

a. $c_k = {(r^{k},{q}_k) / \langle {q}_k,{q}_k\rangle}$

b. $x^{k} = x^{k-1} + c_k{q}_k$

c. $r^{k} = r^{k-1} - c_kA{q}_k$

d. $\beta_k = {\langle r^{k+1},{q}_k\rangle / \langle{q}_k,{q}_k\rangle}$

e. ${q}_{k+1} = r^{k+1} - \beta_k{q}_k$


The previous remark that the listed algorithm is just a summary of the
steps in the solution procedure, and not an efficient algorithm
that should be implemented in its present form, must be repeated here.
In general, we recommend to rely on some high-quality linear algebra library
that offers an implementation of the Conjugate gradient method.

**The computational nature of Conjugate gradient-like methods.**

Looking at the Galerkin and least squares algorithms above,
one notice that the matrix $A$ is only used in matrix-vector
products. This means that it is sufficient
to store only the nonzero entries of $A$.
The rest of the algorithms consists of vector operations of the
type $y \leftarrow ax + y$,
the slightly more general variant $q\leftarrow ax +y$, as well as
inner products.



## A framework based on the error

Let us define the error $e^k = x - x^k$. Multiplying this equation by $A$
leads to the well-known relation between the error and the residual
for linear systems:

<!-- Equation labels as ordinary links -->
<div id="linalg:erroreq"></div>

$$
\begin{equation}
Ae^k = r^k {\thinspace .}
\label{linalg:erroreq} \tag{15}
\end{equation}
$$

Using $r^k = Ae^k$ we can reformulate the Galerkin and
least-squares methods in terms of the error.  The Galerkin method can
then be written

<!-- Equation labels as ordinary links -->
<div id="_auto7"></div>

$$
\begin{equation}
(r^k,{q}_i ) = (Ae^k,{q}_i) = \langle e^k,{q}_i\rangle = 0,\quad
i=0,\ldots,k{\thinspace .}
\label{_auto7} \tag{16}
\end{equation}
$$

For the least-squares method we obtain

<!-- Equation labels as ordinary links -->
<div id="_auto8"></div>

$$
\begin{equation}
(r^k,A{q}_i) = [e^k,{q}_i] =0,\quad i=0,\ldots,k{\thinspace .}
\label{_auto8} \tag{17}
\end{equation}
$$

This means that

$$
\begin{align*}
\langle e^k, v\rangle  &= 0\quad\forall v\in V_{k+1}\hbox{ (Galerkin)}\\
\lbrack e^k, v \rbrack  &= 0\quad\forall v\in V_{k+1}\hbox{ (least-squares)}
\end{align*}
$$

In other words, the error is $A$-orthogonal to the space $V_{k+1}$ in the
Galerkin method, whereas the error is $A^TA$-orthogonal to $V_{k+1}$ in
the least-squares method. When the error is orthogonal to a space, we
find the best approximation in the associated norm to the solution in
that space. Specifically here, it means that for a symmetric and positive
definite $A$, the Conjugate gradient method finds the optimal adjustment
in $V_{k+1}$ of the vector $x^k$ (in the $A$-norm)
in the update for $x^{k+1}$. Similarly, the
least squares formulation finds the optimal adjustment in $V_{k+1}$
measured in the $A^TA$-norm.

A lot of Conjugate gradient-like methods were developed in the 1980s and
1990s, some of the most popular methods do not fit
directly into the framework presented here.  The theoretical
similarities between the methods are covered in the literature, and
we refer to it for algorithms and practical
comments related to widespread methods, such as the SYMMLQ method (for
symmetric indefinite systems), the Generalized Minimal Residual
(GMRES) method, the BiConjugate Gradient (BiCG) method, the
Quasi-Minimal Residual (QMR) method, and the BiConjugate Gradient
Stabilized (BiCGStab) method.  When $A$ is symmetric and positive
definite, the Conjugate gradient method is the optimal choice with
respect to computational efficiency, but when $A$ is nonsymmetric,
the performance of the methods is strongly problem dependent, and one
needs to experiment.


# Preconditioning
<div id="ch:linalg2:preconditioning"></div>


## Motivation and Basic Principles


<div id="ch:linalg2:condno"></div>
The Conjugate Gradient method has been subject to extensive analysis,
and its convergence properties are well understood.
To reduce the initial error $e^0 =x -x^0$ with a factor
$0 <\epsilon\ll 1$ after $k$ iterations, or more precisely,
$||e^k||_{A}\leq\epsilon ||e^0||_{A}$, it can be shown that
$k$ is bounded by

$$
{1\over2}\ln{2\over\epsilon}\sqrt{\kappa},
$$

where $\kappa$ is the ratio of the largest and smallest eigenvalue of
$A$. The quantity $\kappa$  is commonly referred to as
the spectral *condition number*.
Common finite element and finite difference discretizations of
Poisson-like PDEs lead to $\kappa\sim h^{-2}$, where $h$ denotes the
mesh size. This implies that the Conjugate Gradient method converges
slowly in PDE problems with fine meshes, as the number of iterations is
proportional to $h^{-1}$.

To speed up the Conjugate Gradient method, we should
manipulate the eigenvalue distribution. For instance, we could reduce
the condition number $\kappa$. This can be achieved by so-called
*preconditioning*. Instead of applying the iterative method to
the system $Ax =b$, we multiply by a matrix $M^{-1}$ and
apply the iterative method to the mathematically equivalent system

<!-- Equation labels as ordinary links -->
<div id="_auto9"></div>

$$
\begin{equation} M^{-1}Ax = M^{-1}b {\thinspace .} \label{_auto9} \tag{18}
\end{equation}
$$

The aim now is to construct a nonsingular *preconditioning matrix or algorithm*
such that $M^{-1}A$
has a more favorable condition number than $A$.
We remark that we use the notation $M^{-1}$ here to indicate that it should resemble the inverse 
of the matrix $A$. 

For increased flexibility we can write
$M^{-1} = C_LC_R$ and transform the system according to

<!-- Equation labels as ordinary links -->
<div id="_auto10"></div>

$$
\begin{equation}
C_LAC_R y = C_Lb,\quad  y =C_R^{-1}x ,
\label{_auto10} \tag{19}
\end{equation}
$$

where $C_L$ is the *left* and $C_R$ is the *right*
preconditioner. If the original coefficient matrix $A$ is symmetric
and positive definite, $C_L = C_R^T$ leads to preservation of these
properties in the transformed system. This is important when applying
the Conjugate Gradient method to the preconditioned linear system.
Even if $A$ and $M$ are symmetric and positive definite,
$M^{-1}A$ does not necessarily inherit these properties.  It
appears that for practical purposes one can express the iterative
algorithms such that it is sufficient to work with a single
preconditioning matrix $M$ only.  We
shall therefore speak of preconditioning in terms of the left
preconditioner $M$ in the following.

## Use of the preconditioning matrix in the iterative methods

Optimal convergence for the Conjugate Gradient method is
achieved when the coefficient matrix $M^{-1}A$ equals the identity matrix
$I$ and only one iteration is required. In the algorithm we need to perform matrix-vector products
$M^{-1}Au$ for an arbitrary $u\in\mathbb{R}^n$.
This means that we have to solve a linear system with $M$ as coefficient
matrix in each iteration since we implement the product $y =M^{-1}Au$
in a two step fashion: First we compute $v = Au$ and then we
solve the linear system $My =v$ for $y$. The optimal choice $M =A$
therefore involves the solution of $Ay =v$ in each iteration, which
is a problem of the same complexity as our original system $Ax =b$.
The strategy must hence be to compute an $M^{-1} \approx A^{-1} $ such that the
algorithmic operations involved in the inversion of $M$ are cheap.


The preceding discussion motivates the following
demands on the preconditioning matrix $M$:

  * $M^{-1}$ should be a good approximation to $A$,

  * $M^{-1}$ should be inexpensive to compute,

  * $M^{-1}$ should be sparse in order to minimize storage requirements,

  * linear systems with $M$ as coefficient matrix must be efficiently solved.

Regarding the last property, such systems must be solved in
$\mathcal{O}(n)$ operations, that is, a complexity of the same order
as the vector updates in the Conjugate Gradient-like algorithms.
These four properties are contradictory and some sort of compromise
must be sought.

## Classical iterative methods as preconditioners
<div id="ch:linalg:SORprecond"></div>


The simplest possible iterative method for solving $Ax=b$ is

$$
x^{k+1} = x^{k} + r^{k} {\thinspace .}
$$

Applying this method to the preconditioned system
$M^{-1}Ax =M^{-1}b$ results in the scheme

$$
x^{k+1} = x^{k} + M^{-1}r^{k} ,
$$

which is nothing but the formula for classical iterative methods such as
the Jacobi method, the Gauss-Seidel method, SOR (Successive over-relaxation),
and SSOR (Symmetric successive over-relaxation).
This motivates for choosing $M$ as one iteration of these classical methods.
In particular, these methods provide simple formulas for $M$:

  * Jacobi preconditioning: $M =D$.

  * Gauss-Seidel preconditioning: $M = D + L$.

  * SOR preconditioning: $M = \omega^{-1}D +L$.

  * SSOR preconditioning:  $M = (2-\omega )^{-1} \left( \omega^{-1}D + L\right) \left( \omega^{-1}D\right)^{-1} \left( \omega^{-1}D + U\right)$

Turning our attention to the four requirements of the preconditioning
matrix, we realize that the suggested $M$ matrices do not demand
additional storage, linear systems with $M$ as coefficient matrix are
solved effectively in $\mathcal{O}(n)$ operations, and $M$ needs no
initial computation.  The only questionable property is how well $M$
approximates $A$, and that is the weak point of using classical
iterative methods as preconditioners.

The Conjugate Gradient method can only utilize the Jacobi and SSOR
preconditioners among the classical iterative methods, because the $M$
matrix in that case is on the form $M^{-1}=C_LC_L^T$, which is
necessary to ensure that the coefficient matrix of the preconditioned
system is symmetric and positive definite.  For certain PDEs, like the
Poisson equation, it can be shown that the SSOR preconditioner reduces
the condition number with an order of magnitude, i.e., from
$\mathcal{O}(h^{-2})$ to $\mathcal{O}(h^{-1})$, provided we use the
optimal choice of the relaxation parameter $\omega$.  According to
experiments, however, the performance of the SSOR preconditioned
Conjugate Gradient method is not very sensitive to the choice of
$\omega$. 

## Incomplete factorization preconditioners
<div id="linalg:ILU"></div>


Imagine that we choose $M =A$ and solve systems $My =v$ by a
direct method. Such methods typically first compute the LU
factorization $M =\bar L\bar U$ and thereafter perform two
triangular solves. The lower and upper triangular factors $\bar L$ and
$\bar U$ are computed from a Gaussian elimination procedure.
Unfortunately, $\bar L$ and $\bar U$ contain nonzero values, so-called
*fill-in*, in many locations where the original matrix $A$ contains
zeroes. This decreased sparsity of $\bar L$ and $\bar U$ increases
both the storage requirements and the computational efforts related to
solving systems $My =v$.  An idea to improve the situation is to
compute *sparse* versions of the factors $\bar L$ and $\bar U$. This
is achieved by performing Gaussian elimination, but neglecting the
fill-in (!). In this way, we can compute approximate factors $\widehat L$
and $\widehat U$ that become as sparse as $A$.  The storage
requirements are hence only doubled by introducing a preconditioner,
and the triangular solves become an $\mathcal{O}(n)$ operation since
the number of nonzeroes in the $\widehat L$ and $\widehat U$ matrices
(and $A$) is $\mathcal{O}(n)$ when the underlying PDE is discretized
by finite difference or finite element methods.  We call $M
=\widehat L\widehat U$ an *incomplete LU factorization*
preconditioner, often just referred to as the ILU preconditioner.

Instead of throwing away all fill-in entries, we can add them to the
main diagonal. This yields the *modified incomplete LU factorization*
(MILU).  One can also allow a certain amount of fill-in
to improve the quality of the preconditioner, at the expense of
more storage for the factors. Libraries offering ILU preconditioners
will then often have a tolerance for how small a fill-in element in
the Gaussian elimination must be to be dropped, and a parameter that
governs the maximum number of nonzeros per row in the factors.
A popular ILU implementation is the open source C code
[SuperLU](https://fs.hlrs.de/projects/craydoc/docs_merged/books/S-6532-30/S-6532-30.pdf). This preconditioner is available to Python programmers through the
`scipy.sparse.linalg.spilu` function.

The general algorithm for MILU preconditioning follows the steps of
traditional exact Gaussian elimination, except that
we restrict the computations to the nonzero entries in $A$.
The factors $\widehat L$ and $\widehat U$ can be stored directly in the
sparsity structure of $A$, that is, the algorithm overwrites a copy
$M$ of $A$ with
its MILU factorization. The steps in the MILU factorizations are
listed below.


1. Given a sparsity pattern as an index set $\mathcal{I}$,
   copy $M_{i,j}\leftarrow A_{i,j}$, $i,j=1,\ldots,n$

2. for $k=1,2,\ldots, n$

a. for $i=k+1,\ldots,n$

     * if $(i,k)\in\mathcal{I}$ then

a. $M_{i,k}\leftarrow M_{i,k}/M_{k,k}$


     * else

a. $M_{i,k} = 0$


     * $r=M_{i,k}$

     * for $j=k+1,\ldots,n$

       * if $j=i$ then

         * $M_{i,j}\leftarrow M_{i,j} - rM_{k,j} + \sum_{p=k+1}^n
           \left( M_{i,p} - rM_{k,p}\right)$


       * else

         * if $(i,j)\in\mathcal{I}$ then

           * $M_{i,j}\leftarrow M_{i,j} - rM_{k,j}$


         * else

           * $M_{i,j} =0$






We also remark here that the algorithm above needs careful refinement
before it should be implemented in a code. For example, one will not
run through a series of $(i,j)$ indices and test for each of them if
$(i,j)\in\mathcal{I}$. Instead one should run more directly through
the sparsity structure of $A$.


## Preconditioners developed for solving PDE problems

The above mentioned preconditioners are general purpose 
preconditioners that may be applied to solve any linear system
and may speed up the solution process significantly. For linear
systems arising from PDE problems there are however often more
efficient preconditioners that exploit the fact that the matrices
come from discretized PDEs. These preconditioners often enable
the solution of linear systems involving millions of unknowns
in a matter of seconds or less on a regular desktop computer. 
This means that solving these huge linear systems often 
takes less time than printing them file. Moreover, these preconditioners
may be used on parallel computers in a scalable fashion such that
simulations involving billions of degrees of freedom are available
e.g. by cloud services or super-computing centers. 

The most efficient preconditioners are the so-called multilevel preconditioners
(e.g. multigrid  and domain decomposition) which in particular for
Poisson type problems yields error reduction of a factor 10 per iteration 
and hence a typical iteration count of the iterative method of 5-10. 
These preconditioners utilize the fact that different parts of the
solution can be localized or represented on a coarser resolution during
the computations before being glued together to a close match to
the true solution. While a proper description of these methods
are beyond the scope here, we mention that black-box preconditioners
are implemented in open source frameworks such as PETSc, Hypre, and Trilinos
and these frameworks are used by most available finite element software (such as FEniCS). 


While these multilevel preconditioners are available and efficient for PDEs similar
to the Poisson problem there is currently no black-box solver that handles
general systems of PDEs. We mention, however, a promising approach for 
tackling some systems of PDEs, namely the so-called operator preconditioning
framework. This framework utilize tools from functional analysis to deduct
appropriate block preconditioners typically involving blocks composed
of for instance multilevel preconditioners in a systematic construction.   
The framework has extended the use of Poisson type multilevel preconditioners
to many systems of PDEs with success.

ipynb
223130c7e385343da9aa92550588bb0746370ce7
Jupyter Notebook
2021-10-05T10:22:29.000Z
2021-08-05T23:14:15.000Z
47ef0a3592b823ae71a874ee35850114f16b6d8b
MIT
mbarzegary/finite-element-intro
9- cg.ipynb
47ef0a3592b823ae71a874ee35850114f16b6d8b
MIT
mbarzegary/finite-element-intro
9- cg.ipynb
47ef0a3592b823ae71a874ee35850114f16b6d8b
MIT
mbarzegary/finite-element-intro
9- cg.ipynb
2022-02-14T15:20:11.000Z
2021-01-26T13:18:02.000Z