In [8]:
import numpy as np
import scipy as sp
from scipy import linalg as la
import matplotlib.pyplot as plt
import scipy.sparse.linalg
%matplotlib inline
#%load_ext memory_profiler
import matplotlib as mpl
mpl.rcParams['font.size'] = 14
mpl.rcParams['axes.labelsize'] = 20
mpl.rcParams['xtick.labelsize'] = 14
mpl.rcParams['ytick.labelsize'] = 14
M=8
Welcome to another edition of our Jupyter Notebooks. A few notebooks back, we saw that the Conjugate Gradient Method, an iterative method, was very useful to solve $A\,\mathbf{x}=\mathbf{b}$ but it only worked when $A$ was positive definite and symmetric. So now we need an iterative method that works with nonsymmetric linear system of equations, and for that we have the Generalized Minimum Residual Method (GMRes). It works really well for finding the solution of large, sparse (and dense as well), nonsymmetric linear systems of equations. Of course, it will also have trouble for ill-conditioned linear system of equations. But it is really easy to add a left or right or both preconditioners!
Least Squares is used to solve overdetemined linear systems of equations $A\,\mathbf{x} = \mathbf{b}$. That is, for example, a linear system of equations where there are more equations than unknowns. It finds the best $\overline{\mathbf{x}}$ so that it minimizes the euclidean length of $\mathbf{r} = \mathbf{b} - A\,\mathbf{x}$.
So, you might be wondering, what does Least Squares have to do with GMRes? WELL, since you're dying to know, I'll tell you: the backward error of the system in GMRes is minimized at each iteration step using a Least Squares formulation.
GMRes is a member of the family of Krylov methods. It finds an approximation of $\mathbf{x}$ restricted to live on the Krylov sub-space $\mathcal{K_k}$, where $\mathcal{K_k}=\{\mathbf{r}_0, A\,\mathbf{r}_0, A^2\,\mathbf{r}_0, \cdots, A^{k-1}\,\mathbf{r}_0\}$ and $\mathbf{r}_0 = \mathbf{b} - A\,\mathbf{x}_0$ is the residual vector of the initial guess.
The idea behind this method is to look for improvements to the initial guess $\mathbf{x}_0$ in the Krylov space. At the $k$-th iteration, we enlarge the Krylov space by adding $A^k\,\mathbf{r}_0$, reorthogonalize the basis, and then use least squares to find the best improvement to add to $\mathbf{x}_0$.
The algorithm is as follows:
Generalized Minimum Residual Method
$\mathbf{x}_0$ = initial guess
$\mathbf{r}$ = $\mathbf{b} - A\,\mathbf{x}_0$
$\mathbf{q}_1$ = $\mathbf{r} / \|\mathbf{r}\|_2$
for $k = 1, ..., m$
$\qquad \ \ \mathbf{y} = A\,\mathbf{q}_k$
$\qquad$ for $j = 1,2,...,k$
$\qquad \qquad$ $h_{jk} = \mathbf{q}_j^*\,\mathbf{y}$
$\qquad \qquad$ $\mathbf{y} = \mathbf{y} - h_{jk}\, \mathbf{q}_j$
$\qquad$ end
$\qquad \ h_{k+1,k} = \|y\|_2 \qquad$ (If $h_{k+1,k} = 0$ , skip next line and terminate at bottom.)
$\qquad \ \mathbf{q}_{k+1} = \mathbf{y}/h_{k+1,k}$
$\qquad$ Minimize $\left\|\widehat{H}_k\, \mathbf{c}_k - [\|\mathbf{r}\|_2 \ 0 \ 0 \ ... \ 0]^T \right\|_2$ for $\mathbf{c}_k$
$\qquad$ $\mathbf{x}_k = Q_k \, \mathbf{c}_k + \mathbf{x}_0$
end
Now we have to implement it.
In [2]:
# This is a very instructive implementation of GMRes.
def GMRes(A, b, x0=np.array([0.0]), m=10, flag_display=True, threshold=1e-12):
n = len(b)
if len(x0)==1:
x0=np.zeros(n)
r0 = b - np.dot(A, x0)
nr0=np.linalg.norm(r0)
out_res=np.array(nr0)
Q = np.zeros((n,n))
H = np.zeros((n,n))
Q[:,0] = r0 / nr0
flag_break=False
for k in np.arange(np.min((m,n))):
y = np.dot(A, Q[:,k])
if flag_display:
print('||y||=',np.linalg.norm(y))
for j in np.arange(k+1):
H[j][k] = np.dot(Q[:,j], y)
if flag_display:
print('H[',j,'][',k,']=',H[j][k])
y = y - np.dot(H[j][k],Q[:,j])
if flag_display:
print('||y||=',np.linalg.norm(y))
# All but the last equation are treated equally. Why?
if k+1<n:
H[k+1][k] = np.linalg.norm(y)
if flag_display:
print('H[',k+1,'][',k,']=',H[k+1][k])
if (np.abs(H[k+1][k]) > 1e-16):
Q[:,k+1] = y/H[k+1][k]
else:
print('flag_break has been activated')
flag_break=True
# Do you remember e_1? The canonical vector.
e1 = np.zeros((k+1)+1)
e1[0]=1
H_tilde=H[0:(k+1)+1,0:k+1]
else:
H_tilde=H[0:k+1,0:k+1]
# Solving the 'SMALL' least square problem.
# This could be improved with Givens rotations!
ck = np.linalg.lstsq(H_tilde, nr0*e1)[0]
if k+1<n:
x = x0 + np.dot(Q[:,0:(k+1)], ck)
else:
x = x0 + np.dot(Q, ck)
# Why is 'norm_small' equal to 'norm_full'?
norm_small=np.linalg.norm(np.dot(H_tilde,ck)-nr0*e1)
out_res = np.append(out_res,norm_small)
if flag_display:
norm_full=np.linalg.norm(b-np.dot(A,x))
print('..........||b-A\,x_k||=',norm_full)
print('..........||H_k\,c_k-nr0*e1||',norm_small);
if flag_break:
if flag_display:
print('EXIT: flag_break=True')
break
if norm_small<threshold:
if flag_display:
print('EXIT: norm_small<threshold')
break
return x,out_res
In [3]:
A = np.array([[1,1,0],[0,1,0],[0,1,1]])
b = np.array([1,2,3])
x0 = np.zeros(3)
In [4]:
# scipy gmres
x_scipy = scipy.sparse.linalg.gmres(A,b,x0)[0]
# our gmres
x_our, _ = GMRes(A, b)
# numpy solve
x_np= np.linalg.solve(A,b)
# Showing the solutions
print('--------------------------------')
print('x_scipy',x_scipy)
print('x_our',x_our)
print('x_np',x_np)
In [5]:
A = np.array([[0,0,0,1],[1,0,0,0],[0,1,0,0],[0,0,1,0]])
b = np.array([1,0,1,0])
x_our, _ = GMRes(A, b, m=10)
norm_full=np.linalg.norm(b-np.dot(A,x_our))
print(norm_full)
In [6]:
A = np.random.rand(10,10)+10*np.eye(10)
b = np.random.rand(10)
x_our, out_res = GMRes(A, b, m=10,flag_display=True)
norm_full=np.linalg.norm(b-np.dot(A,x_our))
print(norm_full)
In [7]:
plt.figure(figsize=(M,M))
plt.semilogy(out_res,'.k',markersize=20,label='residual')
plt.grid(True)
plt.xlabel(r'$k$')
plt.ylabel(r'$\|\mathbf{b}-A\,\mathbf{x}_k\|_2$')
plt.grid(True)
plt.show()
%timeit%timeit and %memit and verify the results.
In [ ]: