7.3. Convergence of iterative methods
We have seen that both the Jacobi and Gauss-Seidel methods are iterated until the estimates of the solution converge to a given tolerance. In Example 7.1, the Jacobi method required 49 iterations to converge to the solution of a system of linear equations, whereas in Example 7.2 the Gauss-Seidel method required only 20 iterations for the same system. Not all iterative methods converge for every linear system, so it is important to establish conditions under which convergence occurs.
Let \(\mathbf{e}^{(k)} = \mathbf{x} - \mathbf{x}^{(k)}\) be the error between the exact solution \(\mathbf{x}\) and the estimate \(\mathbf{x}^{(k)}\). The error from one estimate to the next is updated using the iteration matrix for the method
\[ \mathbf{e}^{(k+1)} = T\mathbf{e}^{(k)}. \]
Assuming that \(T\) has \(n\) linearly dependent eigenvectors \(\mathbf{v}_1, \ldots, \mathbf{v}_n\) then
\[ \mathbf{e}^{(0)} =\alpha_1 \mathbf{v}_1 +\alpha_2 \mathbf{v}_2 +\cdots +\alpha_n \mathbf{v}_n =\sum_{i=1}^n \alpha_i \mathbf{v}_i, \]
where \(\alpha_i\) are scalars. Since \(\mathbf{v}_i\) are eigenvectors of \(T\) then \(T\mathbf{v}_i = \lambda_i \mathbf{v}_i\) where \(\lambda_i\) are the corresponding eigenvalues, iterating \(\mathbf{e}^{(k)}\) gives
\[\begin{split} \begin{align*}
\mathbf{e}^{(1)} &= T\mathbf{e}^{(0)} =T\left(\sum_{i=1}^n \alpha_i \mathbf{v}_i \right)=\sum_{i=1}^n \alpha_i T\mathbf{v}_i = \sum_{i=1}^n \alpha_i \lambda_i \mathbf{v}_i , \\
\mathbf{e}^{(2)} &=T\mathbf{e}^{(1)} =T\left(\sum_{i=1}^n \alpha_i \lambda_i \mathbf{v}_i \right)=\sum_{i=1}^n \alpha_i \lambda_i T\mathbf{v}_i =\sum_{i=1}^n \alpha_i \lambda_i^2 \mathbf{v}_i , \\
&\vdots \\
\mathbf{e}^{(k+1)} &=\sum_{i=1}^n \alpha_i \lambda_i^{k+1} \mathbf{v}_i .
\end{align*} \end{split}\]
Extracting the first summation term and factorising gives
\[\begin{split} \begin{align*}
\mathbf{e}^{(k+1)} &= \alpha_1 \lambda_1^{k+1} \mathbf{v}_1 +\sum_{i=2}^n \alpha_i \lambda_i^{k+1} \mathbf{v}_i \\
&= \lambda_1^{k+1} \left(\alpha_1 \mathbf{v}_1 +\sum_{i=2}^n \alpha_i \mathbf{v}_i {\left(\frac{\lambda_i }{\lambda_1 }\right)}^{k+1} \right)
\end{align*} \end{split}\]
Let \(\lambda_1\) be the eigenvalue of largest magnitude then \(\lambda_i / \lambda_1 < 1\) and
\[ \begin{align*}
\lim_{k\to \infty } \mathbf{e}^{(k+1)} =\alpha_1 \lambda_1^{k+1} \mathbf{v}_1 .
\end{align*} \]
This means that as the number of iterations increases, the error varies by a factor of \(\lambda_1^{k+1}\) where \(\lambda_1\) is the largest eigenvalue of \(T\) which is known as the spectral radius. For a method to converge to the solution, we require the this limit to tend to zero which only happens when
\[ |\lambda_1| < 1. \]
Definition 7.4 (Spectral radius)
The spectral radius of \(A\) denoted by \(\rho(A)\) is
(7.5)\[ \rho(A) = \max_i(| \lambda_i |). \]
where \(\lambda_i\) are the eigenvalues of \(A\).
Theorem 7.1 (Convergence criterion for a stationary iterative method)
An iterative method of the form
\[ \mathbf{x}^{(k+1)} = T\mathbf{x}^{(k)} + \mathbf{c} \]
converges for every initial approximation (\mathbf{x}^{(0)}) if and only if
\[ \rho(T)<1. \]
Furthermore, smaller values of (\rho(T)) generally lead to faster convergence.
7.3.1. Iteration matrices for the Jacobi and Gauss-Seidel methods
Given a linear system \(A \mathbf{x} = \mathbf{b}\), we can split up the coefficient matrix \(A\) into three matrices such that
\[ A = L + D + U, \]
where \(L\) and \(U\) arethe lower triangular and upper triangular parts of \(A\) and \(D\) is the diagonal matrix, i.e.,
\[\begin{split} \begin{pmatrix}
a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{pmatrix} =
\begin{pmatrix}
0 & 0 & 0 \\
a_{21} & 0 & 0 \\
a_{31} & a_{32} & 0
\end{pmatrix} +
\begin{pmatrix}
a_{11} & 0 & 0 \\
0 & a_{22} & 0 \\
0 & 0 & a_{33}
\end{pmatrix} +
\begin{pmatrix}
0 & a_{12} & a_{13} \\
0 & 0 & a_{23} \\
0 & 0 & 0
\end{pmatrix}. \end{split}\]
Rewriting \(A \mathbf{x} = \mathbf{b}\) using \(L\), \(D\) and \(U\)
(7.6)\[ (L + D + U) \mathbf{x} = \mathbf{b}, \]
which can be rearranged to
\[ \mathbf{x} = -D^{-1}(L + U)\mathbf{x} + D^{-1} \mathbf{b}. \]
The Jacobi method assumes the values on the left-hand side are the new approximations and the values on the right-hand side are the current approximations so
\[ \mathbf{x}^{(k+1)} =-D^{-1}(L + U)\mathbf{x}^{(k)} + D^{-1} \mathbf{b}. \]
The iteraiton matrix \(T_J\) satisfies
\[ \mathbf{x}^{(k+1)} = T \mathbf{x}^{(k)} + \mathbf{c}, \]
so
(7.7)\[ T_J = -D^{-1}(L + U). \]
The iteration matrix for the Gauss-Seidel method can be derived in a similar way, rearranging equation (7.6)
\[ \mathbf{x} = -(L + D)^{-1} U \mathbf{x} + (L + D)^{-1} \mathbf{b}, \]
and assuming the left-hand term are the new approximations and the right-hand terms are the current approximations then
\[ \mathbf{x}^{(k+1)} = -(L + D)^{-1} U \mathbf{x}^{(k)} + (L + D)^{-1} \mathbf{b}, \]
and the iteration matrix for the Gauss-Seidel method is
(7.8)\[ T_{GS} = -(L + D)^{-1}U. \]
Example 7.3
Show that the Jacobi and Gauss-Seidel methods are convergent of the system of linear equations from Example 7.1 (shown below)
\[\begin{split} \begin{align*}
4x_1 +3x_2 &=-2, \\
3x_1 +4x_2 -x_3 &=-8, \\
-x_2 +4x_3 &=14.
\end{align*} \end{split}\]
Solution
The coefficient matrix for this linear system is
\[\begin{split} A = \begin{pmatrix} 4 & 3 & 0 \\ 3 & 4 & -1 \\ 0 & -1 & 4 \end{pmatrix}. \end{split}\]
so
\[\begin{split} \begin{align*}
L &= \begin{pmatrix} 0 & 0 & 0 \\ 3 & 0 & 0 \\ 0 & -1 & 0 \end{pmatrix}, &
D &= \begin{pmatrix} 4 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 4 \end{pmatrix}, &
U &= \begin{pmatrix} 0 & 3 & 0 \\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{pmatrix}.
\end{align*} \end{split}\]
The iteration matrices for the Jacobi and Gauss-Seidel methods are given in equations (7.7) and (7.8) which for this system are
\[\begin{split} \begin{align*}
T_J &= -D^{-1} ( L + U) \\
&= -\begin{pmatrix} 4 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 4 \end{pmatrix} ^{-1} \left(
\begin{pmatrix} 0 & 0 & 0 \\ 3 & 0 & 0 \\ 0 & -1 & 0 \end{pmatrix} +
\begin{pmatrix} 0 & 3 & 0 \\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{pmatrix} \right) \\
&= - \begin{pmatrix} \frac{1}{4} & 0 & 0 \\ 0 & \frac{1}{4} & 0 \\ 0 & 0 & \frac{1}{4} \end{pmatrix}
\begin{pmatrix} 0 & 3 & 0 \\ 3 & 0 & -1 \\ 0 & -1 & 0 \end{pmatrix} \\
&= \begin{pmatrix} 0 & -\frac{3}{4} & 0 \\ -\frac{3}{4} & 0 & \frac{1}{4} \\ 0 & \frac{1}{4} & 0 \end{pmatrix}, \\
T_{GS} &= - (L + D)^{-1} U \\
&= -\left(
\begin{pmatrix} 0 & 0 & 0 \\ 3 & 0 & 0 \\ 0 & -1 & 0 \end{pmatrix} +
\begin{pmatrix} 4 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 4 \end{pmatrix}
\right)^{-1}
\begin{pmatrix} 0 & 3 & 0 \\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{pmatrix}
\\
&= - \begin{pmatrix} 3 & 0 & 0 \\ 3 & 4 & 0 \\ 0 & -1 & 4 \end{pmatrix}^{-1}
\begin{pmatrix} 0 & 3 & 0 \\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{pmatrix} \\
&= \begin{pmatrix} -\frac{1}{4} & 0 & 0 \\ \frac{3}{16} & \frac{1}{4} & 0 \\ \frac{3}{64} & -\frac{1}{16} & \frac{1}{4} \end{pmatrix}
\begin{pmatrix} 0 & 3 & 0 \\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{pmatrix} \\
& =
\begin{pmatrix} 0 & -\frac{3}{4} & 0 \\ 0 & \frac{9}{16} & \frac{1}{4} \\ 0 & \frac{9}{64} & \frac{1}{16} \end{pmatrix}
\end{align*} \end{split}\]
Calculating the spectral radius for these iteration matrices gives \(\rho(T_J )=0.7906\) and \(\rho (T_{GS})=0.625\) which are both less than 1 so both of these methods are convergent for this system. Furthermore, the Gauss-Seidel method will converge faster than the Jacobi method since it has a smaller spectral radius.
7.3.2. Diagonal dominance test
From the previous section we know that an iterative method converges if
\[ \rho(T) < 1, \]
where \(T\) is the iteration matrix and \(\rho(T)\) is the spectral radius of \(T_J\). Since
\[ \rho(T) \leq \| T \|, \]
for any matrix norm, it is sufficient to show that
\[ \| T \|_\infty < 1. \]
where \(\|T\|_\infty\) is the matrix infinity norm
\[ \| T \|_\infty = \max_i \sum_{j=1}^n |(T)_{ij}|. \]
The iteration matrix for the Jacobi method is
\[\begin{split} \begin{align*}
T_J &= -D^{-1}(L + U) \\
&= -
\begin{pmatrix}
\frac{1}{a_{11}} & 0 & \cdots & 0 \\
0 & \frac{1}{a_{22}} & \ddots & \vdots \\
\vdots & \ddots & \ddots & 0 \\
0 & \cdots & 0 & \frac{1}{a_{nn}}
\end{pmatrix}
\begin{pmatrix}
0 & a_{12} & \cdots & a_{1n} \\
a_{21} & 0 & \ddots & \vdots \\
\vdots & \ddots & \ddots & a_{n-1,n} \\
a_{n1} & \cdots & a_{n,n-1} & 0
\end{pmatrix} \\
&=
\begin{pmatrix}
0 & -\frac{a_{12}}{a_{22}} & \cdots & -\frac{a_{1n}}{a_{11}} \\
-\frac{a_{21}}{a_{22}} & 0 & \ddots & \vdots \\
\vdots & \ddots & \ddots & -\frac{a_{n-1,n}}{a_{n-1,n-1}} \\
-\frac{a_{n1}}{a_{nn}} & \cdots & -\frac{a_{n,n-1}}{a_{nn}} & 0
\end{pmatrix}
\end{align*} \end{split}\]
so
\[\begin{split} \begin{align*}
\| T_J \| &= \max_i \sum_{j=1}^n |(T_J)_{ij}| \\
&= \max_i \sum_{j=1,j\ne i}^n \left| \frac{a_{ij}}{a_{ii}} \right| \\
&= \max_i \frac{1}{|a_{ii}|} \sum_{j=1,j\ne i}^n |a_{ij}|.
\end{align*} \end{split}\]
If \(A\) is strictly diagonally dominant so that
\[ |a_{ii}| > \sum_{j=1,j \ne i}^n |a_{ij}|, \qquad i = 1, \ldots, n,\]
then \(\| T_J \|_\infty < 1\) and the Jacobi method converges. It can also be shown that the diaongal dominance test applies to the Gauss-Seidel method.
Theorem 7.2 (Diagonal dominance test for convergence)
The Jacobi and Gauss-Seidel methods are guaranteed to converge for any initial approximation \(\mathbf{x}^{(0)}\) if the coefficient matrix \(A=[a_{ij}]\) is strictly diagonally dominant, that is,
\[ |a_{ii}| > \sum_{j=1,j \neq i}^n |a_{ij}|, \qquad i = 1,\ldots,n. \]
In other words, the magnitude of each diagonal element must be greater
than the sum of the magnitudes of all the other elements in the same row.
Note
Strict diagonal dominance is a sufficient condition for convergence, but it is not necessary. Some matrices that are not diagonally dominant still yield convergent Jacobi iterations.
Example 7.4
Use the diagonal dominance test to determine whether the Jacobi and Gauss-Seidel methods will converge for the following linear system
\[\begin{split} \begin{align*}
x_1 + 4x_2 + x_3 &= 6, \\
x_1 + 2x_2 + 6x_3 &= 9, \\
8x_1 - x_2 + x_3 &= 7.
\end{align*} \end{split}\]
Solution
The coefficient matrix is
\[\begin{split} A =
\begin{pmatrix}
1 & 4 & 1 \\
1 & 2 & 6 \\
8 & -1 & 1
\end{pmatrix} \quad
\begin{matrix} 1 < 4 + 1 \\ 2 < 1 + 6 \\ 1 < 8 + 1 \end{matrix}. \end{split}\]
\(A\) is clearly not diagonally dominant. However, if we swap rows 1 and 3 we get the coefficient matrix
\[\begin{split} A =
\begin{pmatrix}
8 & -1 & 1 \\
1 & 4 & 1 \\
1 & 2 & 6 \\
\end{pmatrix} \quad
\begin{matrix} 8 > 1 + 1 \\ 4 > 1 + 1 \\ 6 > 1 + 2 \end{matrix}. \end{split}\]
which is diagonally dominant.
