Today we will talk about **iterative methods**. List of topics (we will possibly not cover all) includes

- Concept of iterative methods
- Richardson iteration, Chebyshev acceleration
- Jacobi/Gauss-Seidel methods
- Idea of Krylov methods
- Conjugate gradient methods for symmetric positive definite matrices
- Generalized minimal residual methods
- The Zoo of iterative methods: BiCGStab, QMR, IDR, ...
- The idea of preconditioners
- Jacobi/Gauss-Seidel as preconditioner, successive overrelaxation
- Incomplete LU

Suppose we **change** $\tau$ every step, i.e.
$$
x_{k+1} = x_k - \tau_k (A x_k - f).
$$

Then, $$e_{k+1} = (I - \tau_k A) e_k = (I - \tau_k A) (I - \tau_{k-1} A) e_{k-1} = \ldots = p(A) e_0, $$
where $p(A)$ is a **matrix polynomial** (simplest matrix function)

$$
p(A) = (I - \tau_k A) \ldots (I - \tau_0 A),
$$
and $p(0) = 1$.

Plot some Chebyshev polynomials

```
In [1]:
```import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
x = np.linspace(-1, 1, 128)
p = np.polynomial.Chebyshev((0, 0, 0, 0, 0, 0, 0, 1), (-1, 1)) #These are Chebyshev series, a proto of "chebfun system" in MATLAB
plt.plot(x, p(x))

```
Out[1]:
```

```
In [5]:
```import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import scipy as sp
import scipy.sparse
import scipy.sparse.linalg as spla
import scipy
from scipy.sparse import csc_matrix
n = 20
ex = np.ones(n);
lp1 = sp.sparse.spdiags(np.vstack((ex, -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
rhs = np.ones(n)
ev1, vec = spla.eigs(lp1, k=2, which='LR')
ev2, vec = spla.eigs(lp1, k=2, which='SR')
lam_max = ev1[0]
lam_min = ev2[0]
tau_opt = 2.0/(lam_max + lam_min)
fig, ax = plt.subplots()
#plt.close(fig)
niters = 1000
x = np.zeros(n)
res_all = []
for i in xrange(niters):
rr = lp1.dot(x) - rhs
x = x - tau_opt * rr
res_all.append(np.linalg.norm(rr))
#Convergence of an ordinary Richardson (with optimal parameter)
plt.semilogy(res_all)

```
Out[5]:
```

Chebyshev method is efficient when you know that your spectrum is on an interval.

But what happens if the spectrum has not that nice structure?

For example, if the spectrum is contained on **two intervals**, the solution can be

written in terms of **Zolotarev polynomials** and **elliptic functions**, and for three intervals

you have to go to the **hyperelliptic functions**.

The Chebyshev method produces the approximation in the **Krylov subspace** of the form

$$
K(A, f) = \mathrm{Span}(f, Af, A^2 f, \ldots, )
$$
Then we can minimize the **residual**:

over all $x \in \mathcal{K}_n(A, f)$, where the norms can be different (and even strange).

- In was proposed in by Hestenes and Stiefel in 1952.
- In exact arithmetics it should give
**exact solution**at $n$ iterations - In floating-point arithmetics it did not - people thought it is
**unstable**compared to Gaussian elimination - Only decades later it was realized that it is a wonderful
**iterative method**

For non-symmetric matrices, there are two methods that are typically used: GMRES and BiCGStab.

GMRES idea is very simple: we construct orthogonal basis in the Krylov subspace.

Given the basis $q_1, \ldots, q_n$ we compute the residual $r_ n = Ax_n - f$ and orthogonalize it to the basis.

A Python realization is available as `scipy.sparse.linalg.gmres`

(although quite buggy).

Now we can do a short demo

```
In [11]:
```import numpy as np
import scipy.sparse.linalg
n = 100
ex = np.ones(n);
lp1 = -sp.sparse.spdiags(np.vstack((ex, -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
rhs = np.ones(n)
ee = sp.sparse.eye(n)
#lp2 = sp.kron(lp1, ee) + sp.kron(ee, lp1)
#rhs = np.ones(n * n)
res_all = []
res_all_bicg = []
def my_print(r):
res_all.append(r)
def my_print2(x): #For BiCGStab they have another callback, please rewrite
res_all_bicg.append(np.linalg.norm(lp1.dot(x) - rhs))
sol = scipy.sparse.linalg.gmres(lp1, rhs, callback=my_print, restart=100)
plt.semilogy(res_all, marker='x',label='GMRES')
sol2 = scipy.sparse.linalg.bicgstab(lp1, rhs, callback=my_print2)
plt.xlabel('Iteration number')
plt.ylabel('Residual')
plt.semilogy(res_all_bicg, label='BiCGStab', marker='o')
plt.legend(loc='best')

```
Out[11]:
```

```
In [8]:
```scipy.sparse.linalg.gmres?

If $\mathrm{cond}(A)$ is large, we need to use **preconditioners**.

Instead of solving $Ax = f$ we solve

$$ AP^{-1} x = f, $$where $P y = z$ is easy to be solved (i.e. $P$ is diagonal, triangular, circulant ...).

Another possibility to used **left preconditioner**:

$$ P A x = P f.$$

Finally, the **incomplete LU** is typically the method of choice if you do not anything about your sparse matrix.

You start your Gaussian elimination and either:

- Throw away small elements (ILUT)
- Throw away elements that are not in the original sparse matrix ILU($0$)
- Throw away elements which are in the
**neighbourhood**on the graph ILU(k)

ILU for Laplace2D: short demo

```
In [17]:
```import scipy
n = 250
ex = np.ones(n);
lp1 = -sp.sparse.spdiags(np.vstack((ex, -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
ee = scipy.sparse.eye(n, n)
res_all = []
res_all_bicg = []
lp2 = scipy.sparse.kron(lp1, ee) + scipy.sparse.kron(ee, lp1)
rhs = np.ones(n * n)
def my_print(r):
res_all.append(r)
def my_print2(x): #For BiCGStab they have another callback, please rewrite
res_all_bicg.append(np.linalg.norm(lp2.dot(x) - rhs))
M = scipy.sparse.linalg.spilu(lp2, drop_tol=1e-3)
#Create prec
lo = scipy.sparse.linalg.LinearOperator(lp2.shape, lambda x: M.solve(x))
sol = scipy.sparse.linalg.gmres(lp2, rhs, callback=my_print, M=lo, restart=50)
plt.semilogy(res_all, marker='x',label='GMRES')
sol2 = scipy.sparse.linalg.bicgstab(lp2, rhs, callback=my_print2, M=lo)
plt.xlabel('Iteration number')
plt.ylabel('Residual')
plt.semilogy(res_all_bicg, label='BiCGStab', marker='o')
plt.legend(loc='best')

```
Out[17]:
```

```
In [123]:
```from IPython.core.display import HTML
def css_styling():
styles = open("./styles/custom.css", "r").read()
return HTML(styles)
css_styling()

```
Out[123]:
```