> Sparse Matrices and Iterative Solvers > Jacobi and Relaxation Methods

Jacobi and Relaxation Methods

Jacobi Method

The Jacobi method is the simplest of the iterative methods, and relies on the fact that the matrix is diagonally dominant. Starting from the problem definition:

$A\mathbf{x} = \mathbf{b}$

we decompose $A$ in to $A = L + D + U$ , where $L$ is lower triangular, $D$ is diagonal, $U$ is upper triangular.

$A\mathbf{x} = L\mathbf{x} + D\mathbf{x} + U\mathbf{x} = \mathbf{b}$

We then assume that we have an initial guess at the solution $\mathbf{x}^0$ , and try to find a new estimate $\mathbf{x}^1$ . Assuming that the diagonal $D$ dominates over $L$ and $U$ , a sensible choice would be to insert $x^0$ and the unknown $x^1$ into the equation like so:

$L\mathbf{x}^0 + D\mathbf{x}^1 + U\mathbf{x}^0 = \mathbf{b}$

we can rearrange to get an equation for $x^1$ . This is easily solved as we can take the inverse of the diagonal matrix by simply inverting each diagonal element individually:

$D\mathbf{x}_1 = \mathbf{b} - (L+U)\mathbf{x}_0$

Thus we end up with the general Jacobi iteration:

$\mathbf{x}_{k+1} = D^{-1}(\mathbf{b} - (L+U)\mathbf{x}_k)$

Relaxation methods

The Jacobi method is an example of a relaxation method, where the matrix $A$ is split into a dominant part $M$ (which is easy to solve), and the remainder $N$ . That is, $A = M - N$

$M\mathbf{x}_{k+1} = N\mathbf{x}_k + \mathbf{b}$

$\mathbf{x}_{k+1} = M^{-1}N\mathbf{x}_k + M^{-1}\mathbf{b}$

This can be rearranged in terms of the residual $\mathbf{r}_k = \mathbf{b} - A \mathbf{x}_k$ to the update equation

$\mathbf{x}_{k+1} = \mathbf{x}_{k} + M^{-1}\mathbf{r}_k$

For the Jacobi method $M = D$ and $N = -(L + U)$ . Other relaxation methods include Gauss-Seidel, where $M = (D + L)$ and $N = -U$ , and successive over-relaxation (SOR), where $M = \frac{1}{\omega} D + L$ and $N = -(\frac{\omega - 1}{\omega} D + U)$ , where $\omega$ is the relaxation parameter that is within the range $0 \le \omega \le 2$ .

For any relaxation method to converge we need $\rho(M^{-1}N) < 1$ , where $\rho()$ is the spectral radius of $M^{-1} N$ , which is defined as the largest eigenvalue $\lambda$ of a a given matrix $G$ :

$\rho(G) = \max{|\lambda|: \lambda \in \lambda(G)}$

For the SOR method, the relaxation parameter $\omega$ is generally chosen to minimise $\rho(M^{-1}N)$ , so that the speed of convergence is maximised. In some cases this optimal $\omega$ is known, for example for finite difference discretisation of the Poisson equation. However, in many cases sophisticated eigenvalue analysis is required to determine the optimal $\omega$ .

Problems

Question

This exercise involves the manipulation and solution of the linear system resulting from the finite difference solution to Poisson's equation in two dimensions. Let $A$ be a sparse symmetric positive definite matrix of dimension $(N-1)^2 \times (N-1)^2$ created using scipy.sparse (for a given $N$ ) by the function buildA as follows:

import numpy as np
import scipy.sparse as sp

def buildA(N):
    dx = 1 / N
    nvar = (N - 1)**2;
    e1 = np.ones((nvar), dtype=float);
    e2 = np.copy(e1)
    e2[::N-1] = 0
    e3 = np.copy(e1)
    e3[N-2::N-1] = 0
    A = sp.spdiags(
        (-e1, -e3, 4*e1, -e2, -e1),
        (-(N-1), -1, 0, 1, N-1), nvar, nvar
    )
    A = A / dx**2;
    return A

and let $\mathbf{f}_1$ and $\mathbf{f}_2$ be the vectors defined in buildf1 and buildf2

def buildf1(N):
    x = np.arange(0, 1, 1/N).reshape(N, 1)
    y = x.T
    f = np.dot(np.sin(np.pi*x), np.sin(np.pi*y))
    return f[1:,1:].reshape(-1,1)

def buildf2(N):
    x = np.arange(0, 1, 1/N).reshape(N, 1)
    y = x.T
    f = np.dot(np.maximum(x,1-x), np.maximum(y,1-y))
    return f[1:,1:].reshape(-1, 1)

We will consider manipulation of the matrix $A$ and solution of the linear systems $A\mathbf{U}_i=\mathbf{f}_i$ . The solution to this linear system corresponds to a finite difference solution to Poisson's equation $-\nabla^2 u = f$ on the unit square with zero Dirichlet boundary conditions where $f$ is either $\sin(\pi x) \sin (\pi y)$ or $\max(x,1-x) \max(y,1-y)$ . PDEs of this type occur (usually with some additional reaction and or convection terms) very frequently in mathematical modelling of physiological processes, and even in image analysis.

Write a function to solve a linear system using the Jacobi method. In terms of $N$ , how many iterations does it take to converge? (Try $N=4,8,16,32,64$ .)

Expand for solution

Solution

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg
import scipy.optimize
import matplotlib.pylab as plt

def jacobi(A, b, x0=None, tol=1e-5, max_iter=1000):
    if x0 is None:
        x0 = np.zeros_like(b)
    x = np.copy(x0)
    b_norm = np.linalg.norm(b)

    # jacobi method: M = D
    M = A.diagonal().reshape(-1, 1)
    invM = 1/M

    # main relaxation iteration
    for i in range(max_iter):
        r = b - A @ x
        error = np.linalg.norm(r) / b_norm
        if error < tol:
            break
        x += invM * r
    return x, i

num = 20
iterations = np.empty((num, 2), dtype=int)
iterations[:] = np.nan
Ns = np.logspace(0.5, 1.5, num=num, dtype=int)
for j, buildf in enumerate((buildf1, buildf2)):
    for i, N in enumerate(Ns):
        A = buildA(N)
        f = buildf(N)
        max_iter = 10*N
        x, iters = jacobi(A, f, max_iter=max_iter)
        if iters < max_iter:
            iterations[i, j] = iters

plt.plot(Ns, iterations)
plt.xlabel('N')
plt.ylabel('iterations')
plt.show()