5  Iterative techniques for solving systems of linear equations

Consider the system of linear equations \[ A{\mathbf{x}}={\mathbf{b}} \tag{5.1}\] where \[ A=\left( \begin{array}{cccc} a_{11} &a_{12} &\dots &a_{1n} \\ a_{21} &a_{22} &\dots &a_{2n} \\ \vdots &\vdots & &\vdots \\ a_{n1} &a_{n2} &\dots &a_{nn} \end{array}\right), \qquad \mathbf{x}=\left( \begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right), \qquad \mathbf{b}=\left( \begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{array}\right) \tag{5.2}\] for a vector \(\mathbf{x}\) of unknowns. An iterative method for a linear system starts with an initial approximation \(\mathbf{x}^{(0)}\) and generates a sequence of vectors \(\{\mathbf{x}^{(k)}\}\) that converges to the solution \(\mathbf{x}\). In constructing an iterative method, we first convert the system \(A\mathbf{x}=\mathbf{b}\) to the form \[ \mathbf{x}=T\mathbf{x}+ \mathbf{c} \tag{5.3}\] where \(T\) is a fixed matrix and \(\mathbf{c}\) is a fixed vector. Then, we compute a sequence of approximations using the formula \[ \mathbf{x}^{(k+1)}=T\mathbf{x}^{(k)} + \mathbf{c}   \text{ for }   k=0, 1, 2, ... \tag{5.4}\] Note that this can be viewed as the fixed point iteration for function \(\mathbf{g}(\mathbf{x})=T\mathbf{x}^{(k)} + \mathbf{c}\).

Let us say a few words about the convergence of sequences of vectors generated by formula Eq. 5.4 (a more detailed analysis for concrete methods will be given later). Subtracting Eq. 5.3 from Eq. 5.4, we obtain \[ \mathbf{x}^{(k+1)}-\mathbf{x}=T(\mathbf{x}^{(k)}-\mathbf{x}) \tag{5.5}\] from which it follows that \[ \Vert\mathbf{x}^{(k+1)}-\mathbf{x}\Vert=\Vert T(\mathbf{x}^{(k)}-\mathbf{x})\Vert \leq \Vert T\Vert , \Vert\mathbf{x}^{(k)}-\mathbf{x}\Vert , \tag{5.6}\] where we have used the property of the natural matrix norm, Theorem 4.2. It follows from Eq. 5.6 and the contraction mapping theorem that the sequence will converge linearly provided that \(\Vert T\Vert \, \Vert < 1\).

5.1 The Jacobi method

Example 5.1 Consider the system \[ \begin{split} 10x_{1} + 2x_{2}-x_{3} &= 7, \\ x_{1} +8 x_{2} +3x_{3}&= -4, \\ -2x_{1}-x_{2}+10x_{3} &= 9, \end{split} \tag{5.7}\] Its unique solution is \(\mathbf{x}=(1, -1,1)^{t}\). First, we convert the system to the form \(\mathbf{x}=T\mathbf{x}+ \mathbf{c}\) by solving the first equation for \(x_{1}\), the second for \(x_{2}\) and the third for \(x_{3}\): \[ \begin{split} x_{1} &= -\frac{2}{10}x_{2}+\frac{1}{10}x_{3} + \frac{7}{10}, \\ x_{2} &=-\frac{1}{8}x_{1} -\frac{3}{8}x_{3}-\frac{1}{2}, \\ x_{3} &=\frac{2}{10} x_{1}+\frac{1}{10}x_{2} +\frac{9}{10} , \end{split} \tag{5.8}\] Thus, \(\mathbf{x}=T\mathbf{x}+ \mathbf{c}\) with \[ T=\left( \begin{array}{ccc} 0 &-\frac{2}{10} &\frac{1}{10} \\ -\frac{1}{8} &0 &-\frac{3}{8} \\ \frac{2}{10} &\frac{1}{10} &0 \end{array} \right), \qquad \mathbf{c}=\left( \begin{array}{c} \frac{7}{10} \\ -\frac{1}{2} \\ \frac{9}{10} \end{array} \right). \tag{5.9}\] Let the initial approximation be \(\mathbf{x}^{(0)}=(0, 0,0)^{t}\). Then, \[ \begin{split} x_{1}^{(1)} &= -\frac{2}{10}x_{2}^{(0)}+\frac{1}{10}x_{3}^{(0)} + \frac{7}{10}=\frac{7}{10}, \\ x_{2}^{(1)} &=-\frac{1}{8}x_{1}^{(0)} -\frac{3}{8}x_{3}^{(0)}- \frac{1}{2}=- \frac{1}{2}, \\ x_{3}^{(1)} &=\frac{2}{10} x_{1}^{(0)}+\frac{1}{10}x_{2}^{(0)} + \frac{9}{10}=\frac{9}{10} . \end{split} \tag{5.10}\] and \[ \begin{split} x_{1}^{(2)} &= -\frac{2}{10}x_{2}^{(1)}+\frac{1}{10}x_{3}^{(1)} + \frac{7}{10}=0.89, \\ x_{2}^{(2)} &=-\frac{1}{8}x_{1}^{(1)} -\frac{3}{8}x_{3}^{(1)}- \frac{1}{2}=-0.925, \\ x_{3}^{(2)} &=\frac{2}{10} x_{1}^{(1)}+\frac{1}{10}x_{2}^{(1)} + \frac{9}{10}=0.99. \end{split} \tag{5.11}\] One can see that just two iterations give a fairly good approximation to the solution \(\mathbf{x}=(1, -1, 1)^{t}\). Further calculations yield:

\(k\)	\(x^{(k)}_{1}\)	\(x^{(k)}_{2}\)	\(x^{(k)}_{3}\)
0	0.0	0.0	0.0
1	0.7	-0.5	0.9
2	0.89	-0.925	0.99
3	0.984	-0.9825	0.9855
4	0.99505	-0.9925625	0.99855
5	0.9983675	-0.9988375	0.99975375
6	0.999742875	-0.9997035938	0.9997897500

The method we have used is called the Jacobi iterative method. In general, it consists of solving the \(i\)th equation of the system \(A\mathbf{x}=\mathbf{b}\) for \(x_{i}\) to obtain \[ x_{i}=\sum_{j=1, j\neq i}^{n} \left(-\frac{a_{ij}x_{j}}{a_{ii}}\right)+ \frac{b_{i}}{a_{ii}} \ \ \ \text{ for } \ \ i=1,2,...,n \tag{5.12}\] and calculating \(\mathbf{x}^{(k+1)}\) from \(\mathbf{x}^{(k)}\) for \(k\geq 0\) by \[ x_{i}^{(k+1)}=\frac{1}{a_{ii}}\left( \sum_{j=1, j\neq i}^{n} \left(-a{ij}x_{j}^{(k)}\right)+ b_{i}\right)    \text{ for }   i=1,2,...,n. \tag{5.13}\] Let us write the Jacobi method in matrix form \(\mathbf{x}^{(k+1)}=T\mathbf{x}^{(k)} + \mathbf{c}\). To do this, we decompose the matrix \(A\) in the form \(A=D+L+U\), where \(D\) is the diagonal part of \(A\), \(L\) is the strictly lower-triangular part of \(A\), and \(U\) is the strictly upper-triangular part of \(A\): \[ A =\left( \begin{array}{cccc} a_{11} &a_{12} &\dots &a_{1n} \\ a_{21} &a_{22} &\dots &a_{2n} \\ \vdots &\vdots & &\vdots \\ a_{n1} &a_{n2} &\dots &a_{nn} \end{array} \right) =\left( \begin{array}{cccc} a_{11} &0 &\dots &0 \\ 0 &a_{22} &\ddots &\vdots \\ \vdots &\ddots &\ddots &0 \\ 0 &\dots &0 &a_{nn} \end{array} \right) + \left( \begin{array}{cccc} 0 &\dots &\dots &0 \\ a_{12} &\ddots & &\vdots \\ \vdots &\ddots &\ddots &\vdots \\ a_{n2} &\dots &a_{n,n-1} &0 \end{array} \right)+ \left( \begin{array}{cccc} 0 &a_{12}&\dots &a_{1n} \\ \vdots &\ddots &\ddots &\vdots \\ \vdots & &\ddots &a_{n-1,n} \\ 0 &\dots &\dots &0 \end{array} \right). \tag{5.14}\] Now we have \[ (D+L+U)\mathbf{x}=\mathbf{b}   \Leftrightarrow   D\mathbf{x}=-(L+U)\mathbf{x}+\mathbf{b}   \Leftrightarrow   \mathbf{x}=-D^{-1}(L+U)^\mathbf{x}+D{-1}\mathbf{b}. \tag{5.15}\] Thus, we have converted the initial system \(A\mathbf{x}=\mathbf{b}\) to the form \(\mathbf{x}=T_{J}\mathbf{x}+ \mathbf{c}_{J}\) with \(T_{J}=-D^{-1}(L+U)\) and \(\mathbf{c}_{J}=D^{-1}\mathbf{b}\), so that the Jacobi iterative method has the form \[ \mathbf{x}^{(k+1)}=T_{J}^\mathbf{x}{(k)} + \mathbf{c}_{J} . \tag{5.16}\] The Jacobi method requires that \(a_{ii}\neq 0\) for each \(i=1, 2, ..., n\). Therefore, if one of the elements \(a_{ii}\) is zero, we need to reorder the equation so that \(a_{ii}\neq 0\) for each \(i=1, 2, ..., n\).

5.2 The Gauss-Seidel method

How to improve the Jacobi method? In the Example 5.1, we used the iteration procedure \[ \begin{split} x_{1}^{(k+1)} &= -\frac{2}{10}x_{2}^{(k)}+\frac{1}{10}x_{3}^{(k)} + \frac{7}{10}, \\ x_{2}^{(k+1)} &=-\frac{1}{8}x_{1}^{(k)} -\frac{3}{8}x_{3}^{(k)}- \frac{1}{2}, \\ x_{3}^{(k+1)} &=\frac{2}{10} x_{1}^{(k)}+\frac{1}{10}x_{2}^{(k)} + \frac{9}{10} . \end{split} \tag{5.17}\] When we calculate \(x_{2}^{(k+1)}\), we use \(x_{1}^{(k)}\) and \(x_{3}^{(k)}\), which have been calculated at the previous step. We may notice however that this stage we already know \(x_{1}^{(k+1)}\) which is assumed to be a better approximation to \(x_{1}\). It is natural therefore to replace \(x_{1}^{(k)}\) in the second equation by \(x_{1}^{(k+1)}\). Similarly, we may replace \(x_{1}^{(k)}\) and \(x_{2}^{(k)}\) in the third equation by \(x_{1}^{(k+1)}\) and \(x_{2}^{(k+1)}\). This yields the formula \[ \begin{split} x_{1}^{(k+1)} &= -\frac{2}{10}x_{2}^{(k)}+\frac{1}{10}x_{3}^{(k)} + \frac{7}{10}, \\ x_{2}^{(k+1)} &=-\frac{1}{8}x_{1}^{(k+1)} -\frac{3}{8}x_{3}^{(k)} -\frac{1}{2}, \\ x_{3}^{(k+1)} &=\frac{2}{10} x_{1}^{(k+1)}+\frac{1}{10}x_{2}^{(k+1)} +\frac{9}{10} . \end{split} \tag{5.18}\] Eq. 5.18 with the initial approximation \(\mathbf{x}^{(0)}=(0, 0, 0)^{t}\) generate the sequence

\(k\)	\(x^{(k)}_{1}\)	\(x^{(k)}_{2}\)	\(x^{(k)}_{3}\)
0	0.0	0.0	0.0
1	0.7	-0.5875	0.9812500000
2	0.915625	-0.9824218750	0.9848828125
3	0.9949726562	-0.9937026368	0.9996242675
4	0.9987029542	-0.9996969695	0.9997708938
5	0.9999164833	-0.9999036455	0.9999929322
6	0.9999800223	-0.9999948524	0.9999965193

This modification of the Jacobi technique is called the Gauss-Seidel iterative method. Comparing \(\mathbf{x}^{(6)}\) obtained using the Jacobi and Gauss-Seidel methods, we see that the Gauss-Seidel method produces a better approximation to the exact solution \(\mathbf{x}=(1,-1,1)^{t}\).

The general formula for the Gauss-Seidel method is (cf. Eq. 5.13) \[ x_{i}^{(k+1)}=\frac{1}{a_{ii}}\left( -\sum_{j=1}^{i-1} a_{ij}x_{j}^{(k+1)}- \sum_{j=i+1}^{n} a_{ij}x_{j}^{(k)} + b_{i}\right) \tag{5.19}\] for \(i=1, \dots, n\). Let us rewrite the Gauss-Seidel method in matrix form. To do this, we multiply both sides of Eq. 5.19 by \(a_{ii}\) and collect all terms of the \(k\)th iteration. This yields \[ a_{i1}x_{1}^{(k+1)}+a_{i2}x_{2}^{(k+1)}+...+a_{ii}x_{i}^{(k+1)}= -a_{i, i+1}x_{i+1}^{(k)}-...- a_{in}x_{n}^{(k)}+b_{i} \tag{5.20}\] for \(i=1, \dots, n\), or, equivalently, \[ \begin{aligned} &&a_{11}x_{1}^{(k+1)} \qquad\qquad\qquad\qquad\qquad\quad\quad\quad = -a_{12}x_{2}^{(k)}-a_{13}x_{3}^{(k)}-...- a_{1n}x_{n}^{(k)}+b_{1}\\ &&a_{21}x_{1}^{(k+1)} +a_{22}x_{2}^{(k+1)} \qquad\qquad\qquad \quad\quad = \qquad\qquad \ -a_{23}x_{3}^{(k)}-...- a_{2n}x_{n}^{(k)}+b_{2} \\ &&\vdots \\ &&a_{n1}x_{1}^{(k+1)} +a_{n2}x_{2}^{(k+1)}+...+a_{n,n}x_{n}^{(k+1)} \ = \qquad\qquad\qquad\qquad\qquad\qquad \qquad\quad \ \ \ b_{n} . \end{aligned} \tag{5.21}\] Therefore, \[ (D+L)\mathbf{x}^{(k+1)}=-U\mathbf{x}^{(k)}+\mathbf{b}   \Leftrightarrow   \ \mathbf{x}^{(k+1)}=-(D+L)^{-1}U\mathbf{x}^{(k)}+(D+L)^{-1}\mathbf{b}   \ \tag{5.22}\] for \(k=1, 2, \dots.\) Introducing the notation \(T_{G}=-(D+L)^{-1}U\) and \(\mathbf{c}_{G}=(D+L)^{-1}\mathbf{b}\), we rewrite the Gauss-Seidel method in the form \[ \mathbf{x}^{(k+1)}=T_{G}\mathbf{x}^{(k)} + \mathbf{c}_{G} . \tag{5.23}\] Note that a lower-triangular matrix is nonsingular iff its diagonal elements are nonzero. In particular, \(D+L\) is nonsingular (i.e. its inverse exists) iff \(a_{ii}\neq 0\) for each \(i=1, 2, ..., n\).

5.3 The convergence of iterative techniques.

We will use the following notation \[ E_{i}^{(k)} = \vert x_{i}^{(k)} - x_{i}\vert \tag{5.24}\] and \[ E^{(k)} = \Vert \mathbf{x}^{(k)}-\mathbf{x}\Vert_{\infty}=\max\{E_{1}^{(k)}, \dots, E_{n}^{(k)}\}. \tag{5.25}\] Then \(\mathbf{x}^{(k)}\to\mathbf{x}\) is equivalent to \(E^{(k)}\to 0\) as \(k\to\infty\). We assume that \(a_{ii}\neq 0\) for all \(i\) (otherwise both methods will not work).

For the Jacobi method we have \[ x_{i}=\frac{1}{a_{ii}}\left( \sum_{j=1, j\neq i}^{n} \left(-a_{ij}x_{j}\right)+ b_{i}\right) \tag{5.26}\] and \[ x_{i}^{(k+1)}=\frac{1}{a_{ii}}\left( \sum_{j=1, j\neq i}^{n} \left(-a_{ij}x_{j}^{(k)}\right)+ b_{i}\right) \text{ for } i=1,2,...,n \tag{5.27}\] Subtracting Eq. 5.26 from Eq. 5.27 and using the triangle inequality, we obtain \[ \begin{split} E_{i}^{(k+1)}&=\frac{1}{\vert a_{ii}\vert}\left\vert \sum_{j=1, j\neq i}^{n} a_{ij}\left(x_{j}^{(k)}-x_{j}\right)\right\vert \leq \frac{1}{\vert a_{ii}\vert} \sum_{j=1, j\neq i}^{n} \left\vert a_{ij}\right\vert \left\vert x_{j}^{(k)}-x_{j}\right\vert, \\ &=\sum_{j=1, j\neq i}^{n}\frac{\left\vert a_{ij}\right\vert}{\vert a_{ii}\vert}E_{j}^{(k)} \leq \sum_{j=1, j\neq i}^{n}\frac{\left\vert a_{ij}\right\vert}{\vert a_{ii}\vert}E^{(k)} \leq \mu E^{(k)}, \end{split} \tag{5.28}\] where \[ \mu=\max_{i}\sum_{j=1, j\neq i}^{n}\frac{\left\vert a_{ij}\right\vert}{\vert a_{ii}\vert} . \tag{5.29}\] Therefore, \[ E^{(k+1)}\leq \mu , E^{(k)}, \tag{5.30}\] and if \(\mu < 1\), then \(\mathbf{x}^{(k)}\) converges to \(\mathbf{x}\) strongly linearly, with rate of convergence \(\mu\).

Note that \(\mu\) depends only on the coefficient matrix \(A\). Note also that the terms in Eq. 8.34 are the absolute values of the entries of \(T_{J}\), summed in each row, so that \(\mu=\Vert T_{J}\Vert_{\infty}\). In the example, we considered in the previous chapter, \[ A=\left( \begin{array}{ccc} 10 &2 &-1 \\ 1 &8 &3 \\ -2 &-1 &10 \end{array}\right), \qquad \mathbf{b}=\left( \begin{array}{c} 7 \\ -4 \\ 9 \end{array}\right), \qquad T_{J}=\left( \begin{array}{ccc} 0 &-\frac{2}{10} &\frac{1}{10} \\ -\frac{1}{8} &0 &-\frac{3}{8} \\ \frac{2}{10} &\frac{1}{10} &0 \end{array} \right). \tag{5.31}\] So, it follows from Eq. 8.34 that in our example \(\mu = 1/2\). In general, the condition \(\mu < 1\) is equivalent to \[ \sum_{j=1, j\neq i}^{n}\vert a_{ij}\vert < \vert a_{ii}\vert \quad \text{ for \ all } \ \ i. \tag{5.32}\] Matrices \(A\) with this property are said to be strictly diagonally dominant. It follows that every strictly diagonally dominant square matrix is invertible. In fact it is not difficult to prove this directly, using elementary linear algebra.

For the Gauss-Seidel scheme the situation is slightly more complicated. Here we have \[ x_{i}^{(k+1)}=\frac{1}{a_{ii}}\left( -\sum_{j=1}^{i-1} a_{ij}x_{j}^{(k+1)}- \sum_{j=i+1}^{n} a_{ij}x_{j}^{(k)} + b_{i}\right) \tag{5.33}\] for \(i=1,2, \dots, n\). Subtracting Eq. 5.26 from Eq. 5.33, we obtain \[ \begin{split} E_{i}^{(k+1)}&=\frac{1}{\vert a_{ii}\vert}\left\vert \sum_{j=1}^{i-1} a_{ij}\left(x_{j}^{(k+1)}-x_{j}\right)+ \sum_{j=i+1}^{n} a_{ij}\left(x_{j}^{(k)}-x_{j}\right)\right\vert \\ &\leq \frac{1}{\vert a_{ii}\vert}\left( \sum_{j=1}^{i-1} \left\vert a_{ij}\right\vert \left\vert x_{j}^{(k+1)}-x_{j}\right\vert +\sum_{j=i+1}^{n} \left\vert a_{ij}\right\vert \left\vert x_{j}^{(k)}-x_{j}\right\vert\right) \\ &=\sum_{j=1}^{i-1}\frac{\left\vert a_{ij}\right\vert}{\vert a_{ii}\vert}E_{j}^{(k+1)} +\sum_{j=i+1}^{n}\frac{\left\vert a_{ij}\right\vert}{\vert a_{ii}\vert}E_{j}^{(k)} \leq \alpha_{i} E^{(k+1)}+\beta_{i} E^{(k)}, \end{split} \tag{5.34}\] where \[ \begin{split} \alpha_{i}&=\sum{j=1}^{i-1}\frac{\left\vert a_{ij}\right\vert}{\vert a_{ii}\vert},\\ \beta_{i}&=\sum{j=i+1}^{n}\frac{\left\vert a_{ij}\right\vert}{\vert a_{ii}\vert},\\ \alpha_{1}&=\beta_{n}=0. \end{split} \tag{5.35}\] Evidently, \(\mu = \max_{i}(\alpha_{i}+\beta_{i})\). From our analysis of the Jacobi method it is appropriate to assume \(\mu < 1\), so that \(\alpha_{i}+\beta_{i} < 1\) for all \(i\) and also \(1-\alpha_{i} > 0\) for all \(i\) . Suppose now that \(E^{(k+1)}=E_{m}^{(k+1)}\) for some \(m\) (recalling that \(E^{(k+1)}\) is the maximum of the \(E_{i}^{(k+1)}\)). Then we have \[ E^{(k+1)}=E_{m}^{(k+1)}\leq \alpha_{m} E^{(k+1)}+\beta_{m} E^{(k)}. \tag{5.36}\] Hence, \[ E^{(k+1)} \leq \frac{\beta_{m}}{1-\alpha_{m}}E^{(k)} \leq \eta \, E^{(k)}, \tag{5.37}\] where \[ \eta=\max_{i} \frac{\beta_{i}}{1-\alpha_{i}}. \tag{5.38}\] Finally, we note that for each \(i\) we have \[ \alpha_{i}+\beta_{i}-\frac{\beta_{i}}{1-\alpha_{i}}= \frac{\alpha_{i}[1-(\alpha_{i}+\beta_{i})]}{1-\alpha_{i}}\geq \frac{\alpha_{i}[1-\mu]}{1-\alpha_{i}}\geq 0. \tag{5.39}\] This implies that \(\mu\geq \eta\). Indeed, \[ \alpha_{i}+\beta_{i} \geq \frac{\beta_{i}}{1-\alpha_{i}} \quad \Rightarrow \quad \mu = \max{i}(\alpha_{i}+\beta_{i})\geq \frac{\beta_{i}}{1-\alpha_{i}} \quad \Rightarrow \quad \mu \geq \max_{i}\frac{\beta_{i}}{1-\alpha_{i}}=\eta . \tag{5.40}\]

So, the Gauss-Seidel imethod converges strongly linearly with rate \(\eta\), which is at least as fast as that of the Jacobi method. Although straightforward, computation of \(\eta\) is more complicated than \(\mu\). In our example we have

\(i\)	\(\alpha_{i}\)	\(\beta_{i}\)	\(1-\alpha_{i}\)	\(\beta_{i}/(1-\alpha_{i})\)
1	0	3/10	1	3/10
2	1/8	3/8	7/8	3/7
3	3/10	0	7/10	0

Therefore \(\eta = 3/7\approx 0.43\) (compared with \(\mu = 0.5\)).

Thus, we have proved the following:

Theorem 5.1 If \(A\) is strictly diagonally dominant, then for any \(\mathbf{x}^{(0)}\in\mathbb{R}^{n}\), both the Jacobi and Gauss-Seidel methods generate sequences \(\{\mathbf{x}^{(k)}\}\) that converge to the unique solution of \(A\mathbf{x}=\mathbf{b}\).

Note that matrix \(A\) in our example is strictly diagonally dominant, so it is not unexpected that both Jacobi and Gauss-Seidel iterations converge to the solution of the system.

Remark. The above theorem gives us a sufficient condition for convergence of sequences generated by both methods. However, this condition is not necessary. If it is not satisfied, the Jacobi and/or Gauss-Seidel methods may still produce converging sequences. So, if \(A\) is not strictly diagonally dominant, we cannot predict whether the Jacobi method or Gauss-Seidel method will converge.