16. Iterative Solutions of Linear Systems
Dated: 21-04-2025
A linear system can be solved either by matrix factorization or by iteration.
Benefit of the iteration is that, we can stop the solution when we reach a good enough approximation and is a lot faster for coefficient matrices which are large and sparse (containing a lot of \(0\) entries).
General Framework for the Iterative Solution of \(A \vec X = \vec b\)
If \(A^{-1}\) exists and \(\vec x^*\) is the solution to \(A \vec x = \vec b\) then the goal for iterative solution is to produce a sequence of vectors \(\vec x^{(1)}, \vec x^{(2)}, \ldots, \vec x^{(k)}\) where \(k\) is large enough such that \(\vec x^{(k)} \approx \vec x^*\).
For a recursive algorithm, we write \(A = M - N\).
\[A \vec x = \vec b\]
\[\implies (M - N) \vec x = \vec b\]
\[\implies M \vec x = N \vec x + \vec b\]
If a sequence \(\{x^{(k)}\}\) satisfies
\[M \vec x^{(k + 1)} = N \vec x^{(k)} + \vec b \quad k \in \mathbb W\]
and converges to some \(\vec x^*\) then
\[A \vec x^* = \vec b\]
approaches
\[M \vec x^* = N \vec x^* + \vec b\]
There are two iterative methods
Jacobi’s Method
This method assumes that \([a_{ij}] \in A \neq 0\) where \(i = j\).
Choose \(M\) to be a diagonal matrix formed from the entries of \(A\).
Also, as we know
\[N = M - A\]
\[\implies M \vec x^{(k+1)} = (M - A)\vec x^{(k)} + \vec b \quad k \in \mathbb W\]
For simplicity, we take \(\vec x^{(0)}\) (the null vector) as the initial approximation.
Example
\[
\begin{aligned}
10x_1 + \quad x_2 - \quad x_3 &= \quad 18 \\
x_1 + 15x_2 + \quad x_3 &= - \,12 \\
-x_1 + \quad x_2 + 20x_3 &= \quad17
\end{aligned}
\]
Solution
For some \(k\), let
\[\vec x^{(k)} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \langle x_1, x_2, x_3\rangle\]
\[\vec x^{(k + 1)} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \langle y_1, y_2, y_3\rangle\]
\[A = \begin{bmatrix} 10 & 1 & -1 \\ 1 & 15 & 1 \\ -1 & 1 & 20 \end{bmatrix}\]
\[M = \begin{bmatrix} 10 & 0 & 0 \\ 0 & 15 & 0 \\ 0 & 0 & 20 \end{bmatrix}\]
\[N = M - A = \begin{bmatrix} 10 & 0 & 0 \\ 0 & 15 & 0 \\ 0 & 0 & 20 \end{bmatrix} - \begin{bmatrix} 10 & 1 & -1 \\ 1 & 15 & 1 \\ -1 & 1 & 20 \end{bmatrix} = \begin{bmatrix} 0 & -1 & 1 \\ -1 & 0 & -1 \\ 1 & -1 & 0 \end{bmatrix}\]
Apply the recursion
\[M \vec x^{(k+1)} = (M - A)\vec x^{(k)} + \vec b \quad k \in [0, 6]\]
\[
\begin{bmatrix} 10 & 0 & 0 \\ 0 & 15 & 0 \\ 0 & 0 & 20 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} 0 & -1 & 1 \\ -1 & 0 & -1 \\ 1 & -1 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} + \begin{bmatrix} 18 \\ -12 \\ 17 \end{bmatrix}
\]
\[
\Rightarrow
\begin{bmatrix}
10y_1 \\
15y_2 \\
20y_3
\end{bmatrix}
=
\begin{bmatrix}
0x_1 - 1x_2 + 1x_3 \\
-1x_1 + 0x_2 - 1x_3 \\
1x_1 - 1x_2 + 0x_3
\end{bmatrix}
+
\begin{bmatrix}
18 \\
-12 \\
17
\end{bmatrix}
\]
\[
\Rightarrow
\begin{bmatrix}
10y_1 \\
15y_2 \\
20y_3
\end{bmatrix}
=
\begin{bmatrix}
0x_1 - 1x_2 + 1x_3 + 18 \\
-1x_1 + 0x_2 - 1x_3 - 12 \\
1x_1 - 1x_2 + 0x_3 + 17
\end{bmatrix}
\]
\[
\implies
\begin{aligned}
10 y_1 = -x_2 + x_3 + 18\\
15 y_2 = -x_1 - x_3 - 12\\
20 y_3 = x_1 - x_2 + 17
\end{aligned}
\]
\[
\implies
\begin{aligned}
y_1 = \frac{-x_2 + x_3 + 18}{10}\\
y_2 = \frac{-x_1 - x_3 - 12}{15}\\
y_3 = \frac{x_1 - x_2 + 17}{20}
\end{aligned}
\]
1st Iteration
\[k = 0\]
\[\vec {x}^{(0)} = \langle x_1, x_2, x_3\rangle = \langle 0, 0, 0 \rangle\]
\[\vec{x}^{(1)} = \langle y_1, y_2, y_3 \rangle = \langle 18/10, -12/15, 17/20 \rangle = \langle 1.8, -0.8, 0.85 \rangle\]
2nd Iteration
\[k = 1\]
\[\vec {x}^{(1)} = \langle x_1, x_2, x_3\rangle = \langle 1.8, -0.8, 0.85 \rangle\]
\[
\begin{aligned}
y_1 &= [ -(-0.8) + (0.85) + 18 ]/10 = 1.965 \\
y_2 &= [ -(1.8) - (0.85) - 12 ]/15 = -0.9767 \\
y_3 &= [ (1.8) - (-0.8) + 17 ]/20 = 0.98
\end{aligned}
\]
\[\vec {x}^{(2)} = \langle 1.965, -0.9767, 0.98 \rangle\]
Here are results after computing with MATLAB.
\[
\begin{array}{ccccccc}
\vec {x}^{(0)} &
\vec {x}^{(1)} &
\vec{x}^{(2)} &
\vec {x}^{(3)} &
\vec {x}^{(4)} &
\vec {x}^{(5)} &
\vec {x}^{(6)} \\
\begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} &
\begin{bmatrix} 1.8 \\ -0.8 \\ 0.85 \end{bmatrix} &
\begin{bmatrix} 1.965 \\ -0.9767 \\ 0.98 \end{bmatrix} &
\begin{bmatrix} 1.9957 \\ -0.9963 \\ 0.9971 \end{bmatrix} &
\begin{bmatrix} 1.9993 \\ -0.9995 \\ 0.9996 \end{bmatrix} &
\begin{bmatrix} 1.9999 \\ -0.9999 \\ 0.9999 \end{bmatrix} &
\begin{bmatrix} 2.0000 \\ -1.0000 \\ 1.0000 \end{bmatrix}
\end{array}
\]
The Gauss-seidel Method
Same as Jacobi's method except that \(M\) is in lower triangular form.
Example
We are going to use the previous example again.
\[M = \begin{bmatrix} 10 & 0 & 0 \\ 1 & 15 & 0 \\ -1 & 1 & 20 \end{bmatrix}\]
\[N = M - A = \begin{bmatrix} 10 & 0 & 0 \\ 1 & 15 & 0 \\ -1 & 1 & 20 \end{bmatrix} - \begin{bmatrix} 10 & 1 & -1 \\ 1 & 15 & 1 \\ -1 & 1 & 20 \end{bmatrix} = \begin{bmatrix} 0 & -1 & 1 \\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{bmatrix}\]
Applying the recursion, we get
\[
\begin{bmatrix} 10 & 0 & 0 \\ 1 & 15 & 0 \\ -1 & 1 & 20 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} 0 & -1 & 1 \\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} + \begin{bmatrix} 18 \\ -12 \\ 17 \end{bmatrix}
\]
\[
\implies
\begin{bmatrix}
10y_1 & 0 & 0 \\
y_1 & 15y_2 & 0 \\
-y_1 & y_2 & 20y_3
\end{bmatrix}
=
\begin{bmatrix}
-x_2 + x_3 \\
-x_3 \\
0
\end{bmatrix}
+
\begin{bmatrix}
18 \\
-12 \\
17
\end{bmatrix}
\]
\[
\implies
\begin{bmatrix}
10y_1 & 0 & 0 \\
y_1 & 15y_2 & 0 \\
-y_1 & y_2 & 20y_3
\end{bmatrix}
=
\begin{bmatrix}
-x_2 + x_3 + 18 \\
-x_3 - 12 \\
0 + 17
\end{bmatrix}
\]
\[
\implies
\begin{aligned}
10y_1 &= -x_2 + x_3 + 18 & \implies y_1 &= \frac{-x_2 + x_3 + 18}{10} \\
y_1 + 15y_2 &= -x_3 - 12 & \implies y_2 &= \frac{-y_1 - x_3 - 12}{15} \\
-y_1 + y_2 + 20y_3 &= 17 & \implies y_3 &= \frac{y_1 - y_2 + 17}{20}
\end{aligned}
\]
\[
\begin{array}{ccccccc}
\vec {x}^{(0)} &
\vec {x}^{(1)} &
\vec{x}^{(2)} &
\vec {x}^{(3)} &
\vec {x}^{(4)} &
\vec {x}^{(5)} &
\vec {x}^{(6)} \\
\begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} &
\begin{bmatrix} 1.8 \\ -0.92 \\ 0.986 \end{bmatrix} &
\begin{bmatrix} 1.9906 \\ -0.9984 \\ 0.9995 \end{bmatrix} &
\begin{bmatrix} 1.9998 \\ -0.9999 \\ 1.0000 \end{bmatrix} &
\begin{bmatrix} 2.0000 \\ -1.0000 \\ 1.0000 \end{bmatrix} &
\begin{bmatrix} 2.0000 \\ -1.0000 \\ 1.0000 \end{bmatrix} &
\begin{bmatrix} 2.0000 \\ -1.0000 \\ 1.0000 \end{bmatrix}
\end{array}
\]
Comparison of Jacobi’s and Gauss-seidel Method
There is no better solution.
It mostly depends on the situation.
However, Jacobi's method is faster for parallel computing.
Condition for the Convergence of both Iterative Methods
An matrix \(A\) is said to be strictly diagonally dominant if the absolute value of each diagonal entry exceeds the sum of the absolute values of the other entries in the same row.
This guarantees convergence.
The following matrix is not strictly diagonally dominant.
\[
\begin{bmatrix}
-6 & 2 & -3 \\
1 & 4 & -2 \\
3 & -5 & 8
\end{bmatrix}
\qquad
\begin{aligned}
|-6| &> |2| + |-3| \\
|4| &> |1| + |-2| \\
|8| &= |3| + |-5|
\end{aligned}
\]
In such cases, it might be possible to convert such matrices by moving the rows up and down.
References
Read more about notations and symbols.