Today we will talk about iterative methods. List of topics (we will possibly not cover all) includes
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).
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:
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]: