3.2 Iterative methods for solving linear systems

3.2 Iterative methods for solving linear systems#

Gaussian elimination is considered a direct method, which produces the exact solution if we ignore finite precision. There are other methods for solving linear systems that follow the mechanism of iterative methods such as Newton’s method and secant method. Starting from an initial guess of the solution vector \(\boldsymbol{x}^0\), such methods generate a sequence of estimates that converge to the true solution under certain conditions.

Jacobi method#

Consider the following linear system:

(19)#\[\begin{split}\begin{align} 4x_1 - 3x_2 &= -1 \\ 2x_1 + 5x_2 &= 19. \end{align}\end{split}\]

Solving for \(x_1\) from the first equation, and \(x_2\) from the second equation gives:

\[\begin{align*} x_1 &= \frac{-1+3x_2}{4} \\ x_2 &= \frac{19-2x_1}{5}. \end{align*}\]

This suggests an iterative scheme to update the estimated solution if we start from an initial guess \(\boldsymbol{x}^0=(x_1^0, x_2^0)\). If we add a superscript \((k)\) to the \(x\) variables on the right hand side, and \((k+1)\) to those on the left, where \((k)\) and \((k+1)\) denote the iteration numbers, then the iterative scheme, known as Jacobi method, is:

\[\begin{align*} x_1^{(k+1)} &= \frac{-1+3x_2^{(k)}}{4} \\ x_2^{(k+1)} &= \frac{19-2x_1^{(k)}}{5}. \end{align*}\]

The following Python code runs \(10\) steps of Jacobi method starting from the initial guess \([0,0]^T\).

import numpy as np

x = np.array([0,0])
nsteps = 10
for i in range(nsteps):
    x = np.array([(-1+3*x[1])/4, (19-2*x[0])/5])
    print('x1 = {0:10.6f}, x2 = {1:10.6f}'.format(x[0], x[1]))

x1 =  -0.250000, x2 =   3.800000
x1 =   2.600000, x2 =   3.900000
x1 =   2.675000, x2 =   2.760000
x1 =   1.820000, x2 =   2.730000
x1 =   1.797500, x2 =   3.072000
x1 =   2.054000, x2 =   3.081000
x1 =   2.060750, x2 =   2.978400
x1 =   1.983800, x2 =   2.975700
x1 =   1.981775, x2 =   3.006480
x1 =   2.004860, x2 =   3.007290

After \(6\) steps, the solution gets close to the true solution \([2,3]^T\), and stays there, indicating the method works.

Now consider the equivalent linear system, with the two equations swapped:

\[\begin{align*} 2x_1 + 5x_2 &= 19 \\ 4x_1 - 3x_2 &= -1. \end{align*}\]

Derive the Jacobi iterative scheme

\[\begin{align*} x_1^{(k+1)} &= \frac{19-5x_2^{(k)}}{2} \\ x_2^{(k+1)} &= \frac{1+4x_1^{(k)}}{3} \\ \end{align*}\]

and run the scheme for \(10\) iterations.

x = np.array([0,0])
nsteps = 10
for i in range(nsteps):
    x = np.array([(19-5*x[1])/2, (1+4*x[0])/3])
    print('x1 = {0:10.6f}, x2 = {1:10.6f}'.format(x[0], x[1]))

x1 =   9.500000, x2 =   0.333333
x1 =   8.666667, x2 =  13.000000
x1 = -23.000000, x2 =  11.888889
x1 = -20.222222, x2 = -30.333333
x1 =  85.333333, x2 = -26.629630
x1 =  76.074074, x2 = 114.111111
x1 = -275.777778, x2 = 101.765432
x1 = -244.913580, x2 = -367.370370
x1 = 927.925926, x2 = -326.218107
x1 = 825.045267, x2 = 1237.567901

This time Jacobi method diverges! Jacobi method does not work for any linear system. The coefficient matrix in linear system (19) satisfies a property called strictly diagonally dominance, which is defined below:

Definition 9

An \(n\times n\) matrix \(A\) is said to be strictly diagonally dominant if for each \(i\), \(1\le i\le n\), \(|a_{ii}|>\sum_{j\ne i}|a_{ij}|\), i.e., the absolute value of each diagonal elment of \(A\) is strictly greater than the absolute sum of the rest of the elemtents in the same row.

The following theorem provides a sufficient condition for Jacobi method to converge.

Theorem 15

If the coefficient matrix of a linear system is strictly diagonally dominant, then Jacobi method applied to the linear system converges to the unique solution for any initial guess.

It can be easily verified that once we swap the two equations of linear system (19), the coefficient matrix is no longer strictly diagonally dominant.

Jacobi method can be derived from another perspective. Rewrite the linear system \(A\boldsymbol{x}=\boldsymbol{b}\) as

(20)#\[(L+D+U)\boldsymbol{x}=\boldsymbol{b}\]

where \(L\) is the strictly lower triangular part of \(A\), \(D\) is the diagonal of \(A\), and \(U\) is the strictly upper triangular part of \(A\).

Further rewriting Eqn. (20) as

\[\begin{equation*} \boldsymbol{x}=D^{-1}(\boldsymbol{b}-(L+U)\boldsymbol{x}) \end{equation*}\]

and adding superscripts to denote iteration steps:

\[\begin{equation*} \boldsymbol{x}^{(k+1)}=D^{-1}(\boldsymbol{b}-(L+U)\boldsymbol{x}^{(k)}) \end{equation*}\]

we derived Jacobi scheme in matrix form.

The following Python function jacobi implements Jacobi iteration.

# Jacobi method
def jacobi(A, b, x_0, k):
    """
    Perform k steps of Jacobi method
    A: the matrix
    b: the right-hand-side
    x_0: initial guess x0
    k: number of steps
    """
    d = np.diag(A)
    r  = A-np.diag(d)
    # Initialize the solution vector
    x = x_0.copy()
    for j in range(k):
        x = (b-np.dot(r,x))/d
        
    return x

Convergence of Jacobi iteration for linear system (19) is shown in the following animation. The initial guess is \([-1, 1]^T\), and only after a few steps, the solution approaches the true solution \([2,3]\).

import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

Gauss-Seidel method#

In the Jacobi scheme, if we slightly tweak the algorithm by using newly updated solutions immediately in the current iteration, instead of waiting until the next, then the algoirithm is called Gauss-Seidel method. Consider linear system (19) again. If we use updated \(x_1\) instantly to update \(x_2\), then the algorithm can be written as

\[\begin{align*} x_1^{(k+1)} &= \frac{-1+3x_2^{(k)}}{4} \\ x_2^{(k+1)} &= \frac{19-2x_1^{(k+1)}}{5}. \end{align*}\]

Equivalently, Gauss-Seidel iteration can be derived by matrix notation. Eqn. (20) can be rewritten as

\[\begin{equation*} (L+D)\boldsymbol{x} = \boldsymbol{b}-U\boldsymbol{x}. \end{equation*}\]

Add superscripts to distinguish the solutions at the current and next iterations:

\[\begin{equation*} (L+D)\boldsymbol{x}^{(k+1)} = \boldsymbol{b}-U\boldsymbol{x}^{(k)} \end{equation*}\]

which is

\[\begin{equation*} \boldsymbol{x}^{(k+1)} = D^{-1}(\boldsymbol{b}-U\boldsymbol{x}^{(k)}-L\boldsymbol{x}^{(k+1)}). \end{equation*}\]

The sufficient condition for Jacobi method to converge also applies to Gauss-Seidel.

We implement Gauss-Seidel by the following Python function Gauss_Seidel.

# Gauss-Seidel method
def Gauss_Seidel(A, b, x_0, k):
    """
    Run k steps of Gauss_Seidel method
    A: the matrix
    b: the right-hand-side
    x_0: initial guess x0
    k: the number of steps
    """ 
    # Get the size of the system
    n = A.shape[0]
    # Initialize the solution vector
    x = x_0.copy()
    
    for l in range(k):
        for i in range(n):
            sum = 0.
            for j in range(n):            
                if i != j:
                    sum += A[i,j]*x[j]
            x[i] = (b[i]-sum)/A[i,i]   
    return x

Animation of Gauss-Seidel applied to linear system (19) is shown below. Gauss-Seidel converges within only a few steps.

The following animation compares Jacobi and Gauss-Seidel methods applied to linear system (19). It is found that Gauss-Seidel converges faster than Jacobi, which makes sense since Gauss-Seidel utilizes new information more quickly compared to Jacobi.

Successive over-relaxation (SOR)#

Successive over-relaxation is obtained by further modifying the Gauss-Seidel iterative scheme so that a new estimate of solution is obtained by a linear combination of the current solution and the Gauss-Seidel formula, with \(1-\omega\) and \(\omega\) as the coefficients respectively. For linear system (19), the SOR scheme is

\[\begin{align*} x_1^{(k+1)} &= (1-\omega)x_1^{(k)} + \omega\frac{-1+3x_2^{(k)}}{4} \\ x_2^{(k+1)} &= (1-\omega)x_2^{(k)} + \omega\frac{19-2x_1^{(k+1)}}{5}. \end{align*}\]

The additional parameter \(\omega\) is called the relaxation factor. Note that if \(\omega=1\), then SOR is reduced to Gauss-Seidel.

SOR is not guaranteed to converge. A suffcient condition for SOR to converge when applied to a symmetric and positive-definite coefficient matrix \(A\) is \(0<\omega<2\) (refer to Chapter 4 of Kincaid and Cheney [2009] for a proof).

To derive the matrix notation for SOR, first multiply Eqn. (20) by \(\omega\):

\[\begin{equation*} \omega(L+D+U)\boldsymbol{x} = \omega\boldsymbol{b}. \end{equation*}\]

Split \(\omega D\) on the left into two terms:

\[\begin{equation*} (\omega(L+U)+(\omega-1)D+D)\boldsymbol{x} = \omega\boldsymbol{b} \end{equation*}\]

which can be rewritten as

\[\begin{equation*} (\omega L+D)\boldsymbol{x} = \omega\boldsymbol{b}-\omega U\boldsymbol{x} - (\omega-1)D\boldsymbol{x}. \end{equation*}\]

With the superscripts added, the equation can further be written as

\[\begin{equation*} (\omega L+D)\boldsymbol{x}^{(k+1)} = \omega\boldsymbol{b}-\omega U\boldsymbol{x}^{(k)} - (\omega-1)D\boldsymbol{x}^{(k)}. \end{equation*}\]

Rearranging it, we have

\[\begin{equation*} \boldsymbol{x}^{(k+1)} = D^{-1}\omega(\boldsymbol{b}- U \boldsymbol{x}^{(k)}-L \boldsymbol{x}^{(k+1)}) + (1-\omega)\boldsymbol{x}^{(k)}. \end{equation*}\]

The choice of \(\omega\) is important, since it determines how fast SOR converges. In practice, one can run pilot simulations with a wide range of \(\omega\)’s without convergence, use the best \(\omega\), and let it converge.

The following Python code SOR performs SOR for a specified number of steps.

# SOR method
def SOR(A, b, x_0, omega, k):
    """
    Run k steps of the SOR method
    A: the matrix
    b: the right-hand-side
    x_0: initial guess x0
    omega: the relaxation parameter
    k: the number of steps
    """ 
    # Get the size of the system
    n = A.shape[0]
    # Initialize the solution vector
    x = x_0.copy()
    
    for l in range(k):
        for i in range(n):
            sum = 0.
            for j in range(n):            
                if i != j:
                    sum += A[i,j]*x[j]
            x[i] = omega*(b[i]-sum)/A[i,i] + (1.-omega)*x[i]   
    return x

The following animation displays the convergence of SOR applied to linear system (19) with \(\omega=0.5\).

We compare the convergence speeds of Jacobi method, Gauss-Seidel method and SOR with two different \(\omega\)’s in the following animation. It demonstrates that SOR outperforms both Jacobi and Gauss-Seidel when \(\omega=0.8\), while it slows down significantly when \(\omega=0.2\).

Exercise 3.2-1

Apply Jacobi, Gauss-Seidel and SOR (\(\omega=1.2\)) to the following linear system

\[\begin{align*} 6x-3y+z=-20\\ x+5y-2z=12\\ 3x-2y+7z=-24. \end{align*}\]

Let the algorithms stop if the 2-norm of the difference between two successive solutions is less than \(10^{-5}\).

KC09: D. Kincaid and W. Cheney. Numerical Analysis: Mathematics of Scientific Computing. American Mathematical Society, 3rd edition, 2009.

3.2 Iterative methods for solving linear systems

Contents

3.2 Iterative methods for solving linear systems#

Jacobi method#

Gauss-Seidel method#

Successive over-relaxation (SOR)#