The direct methods (Gauss--Jordan) for solving linear systems, using elementary row operations, lead to exact solutions. However, given systems of equations are often subject to errors due to roundoff and other
factors. The Gaussian elimination is sensitive to roundoff errors, and this sensitivity can lead to inaccurate or even wildly wrong answers. Also, even if Gaussian elimination does not go astray, we cannot improve on a solution once we have found it. For
example, if we use Gaussian elimination to calculate a solution to four decimal places,
there is no way to obtain the solution to eight decimal places except to start over again
and work with increased accuracy (we do not discuss this topic and closely related rate of convergence at this moment, which is related to numerical methods).
In contrast, we can achieve additional accuracy with iterative methods simply by
doing more iterations.
The methods, starting from an initial guess x(0), iteratively apply some formulas to detect the solution of the system. For this reason, these methods are often indicated as iterative methods.
In this section, we explore one-step iterative methods that may be beneficial. When applied, they proceed iteratively by successively generating sequences of vectors that approach a solution to a linear system. In many instances (such as when the coefficient matrix is sparse, that is, contains many zero entries), iterative methods can be faster and more accurate than direct methods. Furthermore, the amount of memory required for
a sparse coefficient n×n matrix A is proportional to n. Also, iterative methods can be stopped whenever the approximate solution they generate is sufficiently accurate.
Iterative methods
To demonstrate applications of iterative methods in linear algebra,
we consider the system of n linear equations with n unknowns x1, x2, … , xn (here n is a positive integer):
where coefficients 𝑎i,j of the system and entries bi are given (known) scalars, which we assume to be real numbers. It is convenient to rewrite the system (1) in compact (vector) form:
for some (vector) function f. Any solution of (3) is called a fixed-point for the iteration function f.
Starting with initial guess x(0), an one-step iterative method produces an improved approximation
x(k+1) from the previous approximation x(k), say
It is our hope that the iterative formula above provides a sequence of numbers { x(k) } that converges to the true solution x = c. However, we can stop iteration when we want to. There are known some sufficient conditions on function f that guarantee convergence of fixed point iteration procedure.
If we consider a system of linear equations A x = b, the starting solution x(0),
and the approximated solution at the k-th step x(k), the solution x = c of the system A x = b, then an approximated method is said to converge to the solution of the system when
\[
\lim_{k\to\infty} \left\| {\bf x}^{(k)} - {\bf c} \right\| = {\bf 0} .
\]
On the contrary, if the above limit does not exist, then the method is said to diverge.
When a sequence of errors ek = x(k) − c convergences fast enough to zero,
iterative methods are competitive with direct methods, which require much
more memory and have a higher computational complexity.
There are several references on iterative methods to solve a system of linear
equations; for example, P.G. Ciarlet, J.J. Dongarra et al., C.T. Kelley, G. Meurant, and Y. Saad.
There are infinite many ways to introduce an equivalent fixed point
problem for the given equation (1) or (2).
In order to apply an iterative method for solving system of equations (1), or its equivalent compact form (2), we need to rearrange this equation to the fixed point form (3). This can be achieved by splitting the matrix A = S + M, where S is a nonsingular matrix. S is also called a preconditioner, and its choice is crucial in numerical analysis. Then
the equation A x = b is the same as S x = b −M x. Application of the inverse matrix
S−1 yields the fixed point equation (3), which becomes
There is no guarantee that this method will work. A successful splitting A = S + M satisfies two different requirements:
The new vector x(k+1) should be easy to compute. Hence, S should be a simple
(and invertible!) matrix; it may be diagonal or triangular.
The sequence { x(k) } should converge to the true solution x = c. If we subtract the iteration
in equation (5) from the true equation x = S−1b − S−1M x, the result is a formula involving only the errors ek = c − x(k):
Equation (6) is a difference equation of first order. It starts with the initial error e₀ and after k steps it produces the new error ek. The question of convergence is exactly
the same as the question of stability: x(k) → c exactly when ek → 0.
The foregoing result holds for any standard norm on ℝn because they all equivalent.
Theorem 1:
Let x(0) be an arbitrary vector in ℝn and define x(k+1) = S−1b − S−1Mx(k) for k = 0, 1, 2,, …. If x = c is the solution to Eq.(4), then a necessary and sufficient condition for x(k) → c is that all eigenvalues of matrix S−1M be less than 1, i.e., | λ | < 1.
In other words,
the iterative method in equation (5) is convergent if and only if every
eigenvalue of S−1M satisfies | λ | < 1. Its rate of convergence depends on the maximum size of eigenvalues of matrix S−1M. So the largest eigenvalue will eventually be dominant.
.
We will prove the theorem only in the case where B = −S−1M has n linearly independent eigenvectors. The case where B is not diagonalizable is later. If u1, u2, … , un are n linearly independent eigenvectors of B, we can write
\[
{\bf x}^{(k)} - {\bf c} = c_1 {\bf u}_1 + c_2 {\bf u}_2 + \cdots + c_n {\bf u}_n
\]
and it follows from x(k) − c = Bk(x(0) − c) that
\begin{align*}
{\bf x}^{(k)} - {\bf c} &= {\bf B}^k \left( c_1 {\bf u}_1 + c_2 {\bf u}_2 + \cdots + c_n {\bf u}_n \right)
\\
&= c_1 \lambda_1^k {\bf u}_1 + c_2 \lambda_2^k {\bf u}_2 + \cdots + c_n \lambda_n^k {\bf u}_n .
\end{align*}
Hence, x(k) − c → 0 if and only if
| λi | < 1 for i = 1, 2, … , n.
In this section, we discuss three versions of stationary iterative method: Jacobi's scheme, the Gauss--Seidel method, and the SOR method.
Since they converge slowly, stationary iterative methods are no longer
used alone, but rather used within other algorithms.
The Jacobi method
Carl Jacobi
The simplest iterative method that was first applied for solving systems of linear equations is known as the Jacobi’s method,
Carl Gustav Jacob Jacobi (1804--1851) was a German mathematician who made fundamental contributions to elliptic functions, dynamics, differential equations, determinants, and number theory.
Although much of his
work was in applied mathematics,
Jacobi believed in the importance of
doing mathematics for its own sake.
A fine teacher, he held positions
at the Universities of Berlin and
Königsberg (now it is city of Kaliningrad) and was one of the most
famous mathematicians in Europe.
In numerical linear algebra, the Jacobi method (a.k.a. the
Jacobi iteration method) is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The method was discovered by Carl Gustav Jacob Jacobi in 1845.
The idea of this iteration method is to extract the main diagonal of the matrix of the given vector equation
\[
{\bf A}\,{\bf x} = {\bf b} ,
\]
and consider the rest of matrix A as perturbation:
\[
{\bf A} = \Lambda + {\bf M} .
\]
Here A is the given square n × n matrix, Λ is the diagonal matrix built from A (it is assumed that the diagonal entries of A do not vanish), and M is the matrix obtained from A by annihilating its diagonal entries. So
Since the new equation x = d + B x is of the form x = f(x), it is suitable for iteration procedure: x(k+1) = d + B x(k). Jacobi’s method uses diagonal part of matrix A for converting the given system of equation into fixed point form.
Here is the algorithm for Jacobi iteration method:
Input: initial guess x(0) to the solution, (diagonal dominant) matrix A, right-hand side
vector b, convergence criterion
Output: solution when convergence is reached
Comments: pseudocode based on the element-based formula above
k = 0
while convergence not reached do
for i := 1 step until n do
σ = 0
for j := 1 step until n do
if j ≠ i then
σ = σ + aij xj(k)
end
end
xi(k+1) = (bi − σ) / aii
end
increment k
end
We now need an initial approximation to the solution. It turns out that it does not
matter what this initial approximation is, so we might as well take x₁ = 0, x₂ = 0. We
use these values in Equations (1.1) to get new values of x₁ and x₂:
We repeat this process (using the old values of x₁
and x₂ to get the new values of x₁ and x₂), producing the sequence of approximations
given in the following Table
The successive vectors \( \displaystyle \quad \begin{bmatrix} x_1^{(k)} \\ x_2^{(k)} \end{bmatrix} \quad \) are called iterates. For instance, when k = 5, the fifth iterate is
\( \displaystyle \quad \begin{bmatrix} 1.8110 \\ 0.740849 \end{bmatrix} \quad \)
Here is Python code to produce the result you see above, first using a module, which is used in other examples also.
the 46th iteration: [1.99999996 1. ]
the 47th iteration: [2. 0.99999998]
the 48th iteration: [1.99999998 1. ]
the 49th iteration: [2. 0.99999999]
the 50th iteration: [1.99999999 1. ]
In code above, parameter n_ is for the number of iterations you want - the display is always the last 5 iterations; fp_ is for floating point, the default being "True" - if you input "False" as the last argument, you will get explicit solution, which produces fractions when input coefficients are from ℚ.
Clear[A, M, B, b, diag, d, \[CapitalLambda], x, l, u, x1, x2, x3];
Although exact calculations involve rational numbers, their usage in actual computations is entirely impractical because the number of bits required
to maintain exactness increases exponentially with the number of iterations. Therefore, next we demonstrate Jacobi's iterative method with floating point arithmetic.
Now we repeat calculations by making 50 iterative steps.
It took less than a second for Mathematica to execute the code with 50 steps. We do not present all output because it is expressed with huge numbers not suitable for human's eyes. For instance,
Expressing the rational numbers above as floating point reals, we get the same result, which is more suitable for human eyes.
So, we repeat calculations using floating point representation.
Here are the component parts of the code above, producing the same output in steps.
A = ({
{6, -5},
{4, -7}
});
d = DiagonalMatrix[Diagonal[A]];
l = LowerTriangularize[A, -1];
u = UpperTriangularize[A, 1];
b = ({
{7},
{1}
});
x[1] = ({
{0},
{0}Sdemo3
});
Grid[{{"A", "d", "l", "u", "b", "x0"},
MatrixForm[#] & /@ {A, d, l, u, b, x[1]}}, Frame -> All]
Starting point does not matter for final destination of Jacobi's ietrating scheme. Let's test that with calculations using different starting points, random floating point numbers and integers.
Below, a function is provided that works with a list of starting points.
The eigenvalues of matrix −Λ−1M are λ₁ ≅ 1.44914 and λ₂ ≅ −1.44914. Since the magnitudes of these eigenvalues exceed 1, we conclude that the Jacobi's iteration scheme diveges for this matrix.
We repeat this process (using the old values of x₁
and x₂ to get the new values of x₁ and x₂), producing the sequence of approximations
given in Table
Table of first five Jacobi iterates
n
0
1
2
3
4
5
x₁
0
\( \displaystyle \frac{4}{5} \)
\( \displaystyle \frac{11}{5} = 2.2 \)
\( \displaystyle \frac{13}{25} = 0.52 \)
\( \displaystyle -\frac{341}{50} = -6.82 \)
\( \displaystyle -\frac{823}{250} = -3.292 \)
x₂
0
\( \displaystyle -\frac{5}{2} \)
\( \displaystyle \frac{11}{10} =1.1 \)
\( \displaystyle -\frac{127}{20} = -6.35 \)
\( \displaystyle -\frac{341}{100} = -3.41 \)
\( \displaystyle \frac{2887}{200} = -14.435 \)
It is hard to believe that this iteration scheme provides the correct answer:
x₁ = 2
and x₂ = 1.
the 48th iteration: [-1.08216395e+08 -5.41081974e+07]
the 49th iteration: [-6.49298361e+07 -1.89378693e+08]
the 50th iteration: [-2.27254431e+08 -1.13627216e+08]
the 51th iteration: [-1.36352658e+08 -3.97695257e+08]
the 52th iteration: [-4.77234308e+08 -2.38617154e+08]
the 53th iteration: [-2.86340584e+08 -8.35160041e+08]
A = {{5, -6}, {7, -4}};
b = {4, 10};
LinearSolve[A, b]
Unlike direct methods that converge to the theoretical solution in a finite time, iterative methods provide approximations since they converge,
under some conditions, to the exact solution of the system of linear equations in infinite steps.
The Jacobi method to be applied, matrix A should
not display zeros on its main diagonal.
It is also convenient to organize the unknown values x1, x2, … , xn into a column vector, denoted by x.
To start an iteration method, you need an initial guess vector
Since convergence of Jacobi's scheme does not depend on your initial choice x(0), it could be any vector, for instance, zero.
At the (k + 1)st iteration x(k+1) = d + Bx(k), the vector x(k+1) is calculated by
The vector x(k) is used in calculating x(k+1). Consequently, these two vectors must be
stored separately.
If the diagonal elements of matrix A are much larger than the off-diagonal elements, the
entries of B should all be small and the Jacobi iteration should converge. We introduce a very important definition that reflects diagonally dominant property of a matrix.
A square matrix is said to be diagonally dominant if, for every row of the matrix, the magnitude of the diagonal entry in a row is larger than or equal to the sum of the magnitudes of all the other (non-diagonal) entries in that row. More precisely, the n×n matrix A = [𝑎ij] is diagonally dominant if
\[
\left\vert a_{ii} \right\vert \ge \sum_{j\ne i} \left\vert a_{ij} \right\vert \qquad \forall i.
\]
If A is diagonally dominant, the matrix B = −Λ−1M of the Jacobi iteration will have the property
It follows, then, that if A is diagonally dominant, the Jacobi iteration will converge to
the solution of A x = b.
Theorem 2:
If a system of n linear equations in n variables has a strictly diagonally dominant coefficient matrix, then it has a unique solution and both the Jacobi and the
Gauss--Seidel method converge to it.
Although we don't need to prove this theorem because it is a special (finite-dimensional) case of more general fixed point theorem, we present its proof..
For j = 1, 2, … , n, let
\[
\alpha_j = \sum_{i=1}^{j-1} \left\vert a_{j,i} \right\vert , \qquad \beta_j = \sum_{i=j+1}^{n} \left\vert a_{j,i} \right\vert \qquad\mbox{and} \qquad M_j = \frac{\beta_j}{| a_{j,i} | - \alpha_j} .
\]
Since A is diagonally dominant, it follows that
\[
\left\vert a_{j,j} \right\vert > \alpha_j + \beta_j , \qquad j=1,2,\ldots , n.
\]
Consequently, Mj < 1 for j = 1, 2, … , n. Thus,
\[
M = \max_{1 \le j \le n} M_j < 1.
\]
We will show that
\[
\| {\bf B} \|_{\infty} = \max_{{\bf x} \ne 0} \frac{\| {\bf B}\,{\bf x} \|_{\infty}}{\| {\bf x} \|_{\infty}} \le M < 1 .
\]
Let x be a nonzero vector in ℝn and let y =
B x. Choose k so that
\[
\| {\bf y} \|_{\infty} = \max_{1 \le i \le n} \, | y_i | = | y_k | .
\]
It follows from the definition of B that
\[
{\bf y} = {\bf B} \,{\bf x} = {\bf S}^{-1} {\bf U}\, {\bf x} ,
\]
and hence
\[
{\bf y} = {\bf B}\,{\bf x} = \Lambda^{-1} \left( {\bf L}\,{\bf y} + {\bf U}\, {\bf x} \right) ,
\]
where Λ is the diagonal matrix extracted from A and
L = S − Λ is strictly lower triangular matrix. Comparing the kth coordinates of each side, we see that
\[
y_k = -\frac{1}{a_{k,k}} \left( \sum_{i=1}^{k-1} a_{k,i} y_i + \sum_{i=k+1}^n a_{k,i} x_i \right) ,
\]
and hence
\[
\| {\bf y} \|_{\infty} = \left\vert y_k \right\vert \le \frac{1}{| a_{k,k} |} \left( \alpha_k \| {\bf y} \|_{\infty} + \beta_k \| {\bf x} \|_{\infty} \right) .
\]
From this inequality, it follows
\[
\frac{\| {\bf B}\,{\bf x} \|_{\infty}}{\| {\bf x} \|_{\infty}} =
\frac{\| {\bf y} \|_{\infty}}{\| {\bf x} \|_{\infty}} \le M_k \le M .
\]
So
\[
\| {\bf B} \|_{\infty} = \max_{{\bf x} \ne 0} \frac{\| {\bf B}\,{\bf x} \|_{\infty}}{\| {\bf x} \|_{\infty}} \le M < 1 ,
\]
and hence the iteration will converge to the solution of A x = b.
Example 3:
We consider system of linear equations
\[
\begin{split}
5\,x_1 - 4\, x_2 + x_3 &= 5, \\
5\,x_1 + 7\, x_2 -4\, x_3 &= 21, \\
3\, x_1 - 7\, x_2 + 8\, x_3 &= -9 ,
\end{split}
\tag{3.1}
\]
with matrix not being diagonally dominant:
\[
{\bf A} = \begin{bmatrix} 5 & -4 & \phantom{-}1 \\
5 & \phantom{-}7 & -4 \\ 3 & -7 & \phantom{-}8 \end{bmatrix} .
\]
Indeed,
\[
5 = |-4| +1, \qquad 7 < 5 + |-4| = 9, \qquad 8 < 3 + |-7| = 10 .
\]
We can check that the given system has the unique solution x₁ = 2, x₂ = 1, x₃ = −1.
A = {{5, -4, 1}, {5, 7, -4}, {3, -7, 8}};
b = {5, 21, -9};
LinearSolve[A, b]
{2, 1, -1}
However, the code above will not work if you define vector b as column vector. To overcome this, use the following modification.
b = ({
{5},
{21},
{-9}
});
LinearSolve[A, b][[All,1]]
Upon extracting the diagonal matrix from A, we split the given system (3.1) into the following form:
\[
\left( \Lambda + {\bf M} \right) {\bf x} = {\bf b} \qquad \iff \qquad {\bf x} =
\Lambda^{-1} {\bf b} - \Lambda^{-1} {\bf M} \,{\bf x} ,
\]
where
\[
{\bf d} = \Lambda^{-1} {\bf b} = \begin{bmatrix} \frac{1}{5} & 0 & 0 \\ 0 & \frac{1}{7} & 0 \\ 0 & 0 & \frac{1}{8} \end{bmatrix} \begin{pmatrix} \phantom{-}5 \\ 21 \\ -9 \end{pmatrix} = \begin{bmatrix} 1 \\ 3 \\ - \frac{9}{8} \end{bmatrix} , \qquad {\bf B} = - \Lambda^{-1} {\bf M} = \begin{bmatrix} 0 & \frac{4}{5} & - \frac{1}{5} \\ - \frac{5}{7} & 0 & \frac{4}{7} \\ - \frac{3}{8} & \frac{7}{8} & 0 \end{bmatrix} .
\]
Matrix B has one real eigenvalue λ₁ ≈ -0.362726,
and two complex conjugate eigenvalues λ₂ ≈ 0.181363 +
0.308393j, λ₃ ≈ 0.181363 - 0.308393j, where j is the imaginary unit of the complex plane ℂ, so j² = −1.
Since all eigenvalues have absolute value less than 1, we expect Jacobi's method to be convergent,
Sqrt[0.181363^2 + 0.308393^2]
0.357769
The code above shows that the absolute value of complex number λ = x + jy, with x = 0.181363 and y = 0.308393, is
\( \displaystyle |\lambda | = | x + {\bf j}\,y | = \sqrt{x^2 + y^2} \approx 0.357769 < 1 . \)
Although we know exact solution, x₁ = 2, x₂ = 1,
x₃ = −1, of the given system of linear equations (3.1), we know that convergence of Jacobi's iteration scheme does not depend on the initial guess. So we start with zero:
\[
{\bf x}^{(0)} = \left[ 0, 0, 0 \right] .
\]
d = {1,3,-9/8};
B= {{0, 4/5, -1/5},{-5/7,0,4/7},{-3/8,7/8,0}};
x3 = d + B.x2
{117/56, 59/56, -(67/64)}
N[%]
{2.08929, 1.05357, -1.04688}
\[
{\bf x}^{(3)} = {\bf d} + {\bf B}\, {\bf x}^{(2)} = \begin{bmatrix} \frac{117}{56} \\ \frac{59}{56} \\ -\frac{67}{64} \end{bmatrix} =
\begin{bmatrix} 2.08929 \\ 1.05357 \\ -1.04688 \end{bmatrix} .
\]
Now it is right time to organize our calculations in a loop.
Theorem 3:
If the Jacobi or the Gauss-Seidel method converges for a system of n linear
equations in n unknowns, then it must converge to the solution of the system.
The idea of the proof is quite simple: to apply an iterative method, we first transferred the linear algebra problem A x = b to a fixed point problem x = f(x), and only then applied the iterative method to the latter. To prove Theorem 3, we go backward: from the fixed point equation we go to the linear algebra one.
Convergence means that "as iterations increase, the values of the iterates get
closer and closer to a limiting value:' This means that x(k) converge to r, As usual, we denote by x = c the true solution of the system A x = b.
We must show that r = c. Recall that the Jacobi's (GS is similar) iterates satisfy the equation
\[
\Lambda\,{\bf x}^{(k+1)} = {\bf b} + \left( \Lambda - {\bf A} \right) {\bf x}^{(k)} , \qquad k=0,,2,\ldots .
\tag{P3.1}
\]
The Gauss-Seidel method leads to a similar equation where A is further decomposed into upper and lower triangular parts. Application of limit to both sides of Eq.(P3.1), we obtain
\[
\Lambda\,{\bf r} = {\bf b} + \left( \Lambda - {\bf A} \right) {\bf r} \qquad \iff \qquad {\bf 0} = {\bf b} - {\bf A}\,{\bf r} .
\]
Therefore, vector \( \displaystyle {\bf r} = \lim_{k\to \infty} {\bf x}^{(k)} \) satisfies the equation A x = b, so it must be equal to c.
Example 4:
Let us reconsider the system four equations with four unknowns:
Although we know one of solutions, x₁ = 2, x₂ = 1, x₃ = 3, x₄ = 4, we demonstrate the iterative method---Jacobi's one. Since the nullspace of A is spanned on the vector [12, 31, 124, 120], the general solution of Eq.(4.1) depends on real parameter t ∈ ℝ:
\[
{\bf x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} 2 \\ 1 \\ 3 \\ 4 \end{bmatrix} + t \begin{bmatrix} 12 \\ 31 \\ 124 \\ 120 \end{bmatrix} \in \mathbb{R}^{4\times 1} .
\]
The Jacobi scheme generates approximants according to the formula
\[
{\bf x}^{(k+1)} = {\bf d} + {\bf B}\,{\bf x}^{(k)} , \qquad k=0, 1, 2, \ldots .
\tag{4.3}
\]
Eigenvalues of matrix B can be determined with Mathematica
N[Eigenvalues[B]]
{1., -0.576902 + 0.664232 I, -0.576902 - 0.664232 I, 0.153804}
Since one of the eigenvalues λ₁ = 1, we do not expect convergence of Jacobi's iteration procedure. Other eigenvalues λ₂ = 0.153804, λ₃ = -0.576902 + 0.664232j, λ₄ = -0.576902 − 0.664232j, have magnitudes less than 1.
Sqrt[0.5769020858621225^2 + 0.6642324328522395^2]
0.879784
Starting with zero vector x(0) = 0, we obtain the first approximant to be x(1) = d. The second and third approximations become
\[
{\bf x}^{(2)} = {\bf d} + {\bf B}\,{\bf x}^{(1)} = {\bf d} + {\bf B}\,{\bf d} = \begin{pmatrix} \frac{209}{63} \\ -\frac{89}{112} \\ \frac{1}{2} \\ \frac{37}{84} \end{pmatrix}
\]
and
\[
{\bf x}^{(3)} = {\bf d} + {\bf B}\,{\bf x}^{(2)} = {\bf d} + {\bf B} \left( {\bf d} + {\bf B}\,{\bf b} \right) = {\bf b} + {\bf B}\,{\bf b} + {\bf B}^2 {\bf b} = \begin{pmatrix} \frac{40}{189} \\ \frac{175}{289} \\ -\frac{1235}{1764} \\ \frac{425}{1512} \end{pmatrix} .
\]
diag = {{9, 0, 0, 0}, {0, 8, 0, 0}, {0, 0, 7, 0}, {0, 0, 0, 6}};
M = A - diag ;
d = {1, 7/8, -(11/7), 3/2};
B = -Inverse[diag] . M ;
x1 = d + B . d
{209/63, -(89/112), 1/2, 37/84}
x2 = d + B . d + B . B . d
{40/189, 175/288, -(1235/1764), 425/1512}
Now it is right time to organize our calculations in a loop.
The numerical values in the table above confirm that Jacobi's iteration diverges for this system of equations.
We are going to make a numerical experiment and check whether Jacobi's iteration depends on the initial guess, which we take closer to the true values:
x(0) = [1.5, 0.5, 2.5, 3.5]. Using iteration (4.3), we calculate Jacobi's approximants.
This table of iterative values shows that the first coordinate plays around true value x₁ = 2, but we cannot claim convergence of other coordinates, except probably x₂ = 1.
Input: square matrix A and vector bn is the size of A and vector b
whle precision condition do
for i=1 : n do
s=0
for j=1 : n do
if j ≠ i then
s = s + 𝑎i,jxj
endif
endfor
yi = (bi−s)/𝑎i,i
endfor
x=y
endwhile
Gauss--Seidel method
Philipp Seidel
Carl Gauss
In numerical linear algebra, the Gauss–Seidel method, also known as the Liebmann method or the method of successive displacement, is an iterative method used to solve a system of linear equations. It is named after the German mathematicians Carl Friedrich Gauss (1777--1855) and Philipp Ludwig von Seidel (1821--1896), and is similar to the Jacobi method. Though it can be applied to any matrix with non-zero elements on the diagonals, convergence is only guaranteed if the matrix is either strictly diagonally dominant, or symmetric and positive definite. It was only mentioned in a private letter from Gauss to his student Gerling in 1823. A publication was not delivered before 1874 by Seidel
Philipp Ludwig
von Seidel
worked in analysis, probability
theory, astronomy, and optics.
Unfortunately, he suffered from
eye problems and retired at a young
age. The paper in which he described
the method now known as Gauss
Seidel was published in 1874.
The Gauss-Seidel method, which is a modification of the Jacobi algorithm,
uses each new value as soon as it is available. This method is based on splitting the matrix A of Eq.(2) into three parts:
Therefore, in general formula (4), M = U is the strictly upper triangular matrix extracted from A, and S is a lower triangular part of A. Then we have an equivalent to Eq.(2) system of equations
Although structure of the inverse to the lower triangular matrix S = Λ + L can be written explicitly, its practical implementation becomes problematic with size of matrix growth. The n × n matrix S−1 contains n divisions, and its entries contain alternated series, which could be a nightmare for numerical evaluations.
Instead of standard approach (5), the Gauss--Seidel method suggests to
extract the strictly lower triangular matrix from S,
Now we demonstrate the iterative procedure above in more detail by showing step by step its application for solving the system of linear equations (2). The first coordinate of x(k+1) is given by Jacobi's formula
It is interesting to compare the Jacobi and Gauss–Seidel iterations noting that the latter is using the coordinates of x(k+1) as soon as
they are calculated rather than in the next iteration. The program for the Gauss–Seidel
iteration is actually simpler than the program for the Jacobi iteration. The vectors x(k)
and x(k+1) are both stored in the same vector, x. As a coordinate of x(k+1) is calculated,
it replaces the corresponding coordinate of x(k).
Here is algorithm for Gauss-Seidel iteration scheme:
algorithm Gauss–Seidel method is
inputs: A, b
output: φ
Choose an initial guess φ to the solution
repeat until convergence
for i from 1 until n do
σ ← 0
for j from 1 until n do
if j ≠ i then
σ ← σ + 𝑎i,jφj
end if
end (j-loop)
φi ← (bi − σ) / 𝑎i,i
end (i-loop)
check if convergence is reached
end (repeat
Another algorithm using while loop:
input A and bn is the size of A
while precision conditions do
for i = 1 : n do
s = 0
for j = 1 : n do
if j ≠ i then
s = s + 𝑎i,jxj
end if
end for
xi = (1/𝑎i,i) (bi − s)
end for
end while
Example 5:
We consider a linear system with three unknowns:
\[
\begin{bmatrix} \phantom{-}4 & 2 & -2 \\ -2 & 5 & \phantom{-}2 \\ \phantom{-}4 & 3 & -2 \end{bmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} 12 \\ 4 \\ 14 \end{pmatrix}
\tag{5.1}
\]
that has the true solution x = 1, y = 2, z = −2.
A = {{4, 2, -2}, {-2, 5, 2}, {4, 3, -2}};
b = {12, 4, 14};
LinearSolve[A, b]
{1, 2, -2}
If you use column vector b instead of row vector, LinearSolve command should be applied as follows.
The errors of the Gauss--Seidel iterations at 50-th run are
\[
\left\vert x_1^{(50)} -1 \right\vert = 0.0602767 , \qquad \left\vert x_2^{(50)} -2 \right\vert = 0.0181273 , \qquad \left\vert x_3^{(50)} +2 \right\vert = 0.095263 .
\]
As a result of this numerical experiment we conclude that both methods, Jacobi and GS, are very fast and provide approximately the same accuracy, but Jacobi's iteration scheme gives slightly better approximation. Also, a usage of floating foint instead of exact operations with rational numbers does not effect the time of execution.
End of Example 5
Example 6:
Let us reconsider the system from Example 3:
\[
\begin{split}
5\,x_1 - 4\, x_2 + x_3 &= 5, \\
5\,x_1 + 7\, x_2 -4\, x_3 &= 21, \\
3\, x_1 - 7\, x_2 + 8\, x_3 &= -9 ,
\end{split}
\tag{3.1}
\]
with matrix not being diagonally dominant:
\[
{\bf A} = \begin{bmatrix} 5 & -4 & \phantom{-}1 \\
5 & \phantom{-}7 & -4 \\ 3 & -7 & \phantom{-}8 \end{bmatrix} .
\]
Now it is right time to organize our calculations in a loop.
This example shows that the Gauss--Seidel method converges very slowly: after 50 runs it gives the answer with about 25% accuracy. So Jacobi's method performs much better.
End of Example 6
Example 7:
Let us reconsider the system four equations with four unknowns:
Matrix B has rank 3 and all its eigenvalues are real numbers: λ₁ = 1, λ₂ = 0, λ₃ ≈ -0.874609, λ₄ ≈ 0.181487. Therefore, we do not expect convergence of Gauss-Seidel iterations.
Upon choosing initial guess to be zero, we perform iterations according to the difference equation (6.4).
Now it is right time to organize our calculations in a loop.
Therefore, we conclude that except the first coordinates, all other entries do not converge.
The following table presents Gauss--Seidel iterations with different initial guess.
x0 = {1.5, 0.5, 2.5, 3.5} {1.5, 0.5, 2.5, 3.5}
{1.5, 0.5, 2.5, 3.5}
x1 = d + B . x0
{1.61111, 0.923611, 2.41667, 3.59259}
Table of Jacobi iterates
n
0
1
2
3
4
5
10
x₁
1.5
1.61111
2.0751
1.81082
2.06765
1.84769
2.00357
x₂
0.5
0.923611
0.778807
0.946414
0.807257
0.930314
0.843734
x₃
2.5
2.41667
2.54586
2.44308
2.53482
2.45492
2.51133
x₄
3.5
3.59259
3.44905
3.56324
3.46131
3.550083.55008
3.48741
We repeat calculations with floating point arithmetic for 50 runs.
End of Example 7
Successive overrelaxation or SOR method
In numerical linear algebra, the method of successive over-relaxation (SOR) is a variant of the Gauss–Seidel method for solving a linear system of equations, resulting in faster convergence. It was devised simultaneously by David M. Young Jr. (1923--2008) and by Stanley P. Frankel (1919--1978) in 1950 for the purpose of automatically solving linear systems on digital computers.
When SOR method is applied to the traditional fixed point equation (3), it generates a sequence of approximations according to the formula
where values of overrelaxation factor ω > 1 are used to speed up convergence of a slow-converging process, while values of ω < 1 are often used to help establish convergence of a diverging iterative process or speed up the convergence of an overshooting process. With ω = 1, we recover Gauss-Seidel method. The optimal choice of ω never exceeds 2. It is often in the neighborhood of 1.9.
To describe application of the overrelaxation for solving system of linear equations (1) or (2), we follow Gauss-Seidel approach and decompose matrix A of coefficients into sum of three matrices:
\[
{\bf A} = \Lambda + {\bf L} + {\bf U} ,
\]
where Λ is the diagonal matrix extracted from A, L is the strictly lower triangular part of A, and U is the strictly upper triangular part of A:
More specifically, if x(k) is the solution at the step
k and xGS(k+1) is the update
to the step k + 1 according to the Gauss--Siedel’s method, the update formula of the
SOR method is:
where ω is a real parameter to be set. Obviously if ω = 1, the SOR algorithm degenerates
into the Gauss-Siedel’s method. The explicit update formula of the SOR method can
be simply obtained from that of Gauss-Seidel’s method by adding the contribution
due to x(k):
Inputs: A, b, ω
Output: φ
Choose an initial guess φ to the solution
repeat until convergence
for i from 1 until n do
set σ to 0
for j from 1 until n do
if j ≠ i then
set σ to σ + 𝑎i,j φj
end if
end (j-loop)
set φi to (1 − ω)φi + ω(bi − σ) / 𝑎i,i
end (i-loop)
check if convergence is reached
end (repeat)
Another algorithm using while loop:
input A and b
input ω
n is the size of A
while precision conditions do
for i = 1 : n do
s = 0
for j = 1 : n do
if j ne; i then
s = s + 𝑎i,j xj
end if
end for
xi = ω (bi − s)/ 𝑎i,i + (1 − ω ) xi
end for
end while
Matlab algorithm:
function [x,r,k]=RM(A,b,Omega,kmax,eps1)
% kmax maximum number of iterations;
% eps1 given tolerance;
% Omega: relaxation parameter
% stopping criteria: norm(r)/norm(b)<=eps1 or k>kmax
% Output: x the kth iterate, r the kth residual
% k the number of iterations
% Initial guess is x=0; r=b
n=length(b);r=b;x=zeros(n,1);nb=norm(b);
E=tril(A,-1);D=diag(diag(A));M=(1/Omega)*D-E;
k=1;mr=1;mx=1;e=1;
while e>eps1 & k<=kmax
dx=M\r;
x=x+dx;% saxpy operation
r=r-A*dx;% GAXPY operation
e=norm(r)/nb;
k=k+1;
end
Theorem 4 (Kahan, 1958):
A necessary condition for the SOR method to converge is ω ∈ (0, 2).
Example 8:
Let us reconsider the system from Example 3:
\[
\begin{split}
5\,x_1 - 4\, x_2 + x_3 &= 5, \\
5\,x_1 + 7\, x_2 -4\, x_3 &= 21, \\
3\, x_1 - 7\, x_2 + 8\, x_3 &= -9 ,
\end{split}
\tag{3.1}
\]
with matrix not being diagonally dominant:
\[
{\bf A} = \begin{bmatrix} 5 & -4 & \phantom{-}1 \\
5 & \phantom{-}7 & -4 \\ 3 & -7 & \phantom{-}8 \end{bmatrix} .
\]
The SOR scheme yields
\begin{align*}
x_1^{(k+1)} &= \frac{\omega}{5} \left( 5 + 4\,x_2^{(k)} - x_3^{(k)} \right) + \left( 1 - \omega \right) x_1^{(k)} ,
\\
x_2^{(k+1)} &= \frac{\omega}{7} \left( 21 + 4\, x_3^{(k)} - 5 \, x_1^{(k+1)} \right) + \left( 1 - \omega \right) x_2^{(k)} ,
\\
x_3^{(k+1)} &= \frac{\omega}{8} \left( -9 - 3\,x_1^{(k+1)} +7 \, x_2{(k+1)} \right) + \left( 1 - \omega \right) x_3^{(k)} ,
\end{align*}
Now it is right time to organize our calculations in a loop.
Table of SOR iterates with ω = 0.7
n
0
1
2
3
4
5
x₁
0
\( \displaystyle \frac{}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{1}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₂
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₃
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
Table of SOR iterates with ω = 1.7
n
0
1
2
3
4
5
x₁
0
\( \displaystyle \frac{}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{1}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₂
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₃
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
End of Example 8
Example 9:
Let us reconsider the system four equations with four unknowns:
The matrix of this system
\[
{\bf A} = \begin{bmatrix} \phantom{-}9 & -8 & \phantom{-}5 & -4 \\ \phantom{-}1 & \phantom{-}8 & -5 & \phantom{-}3 \\ -2 & -4 & \phantom{-}7 & -6 \\ \phantom{-}2 & -4 & -5 & \phantom{-}6 \end{bmatrix}
\tag{4.2}
\]
is singular of rank 3:
The SOR scheme yields
\begin{align*}
x_1^{(k+1)} &= \frac{\omega}{9} \left( 9 + 8\,x_2^{(k)} -5\, x_3^{(k)} + 4\, x_4^{(k)} \right) + \left( 1 - \omega \right) x_1^{(k)} ,
\\
x_2^{(k+1)} &= \frac{\omega}{8} \left( 7 + 5\, x_3^{(k)} - 3\, x_4^{(k)} - x_1^{(k+1)} \right) + \left( 1 - \omega \right) x_2^{(k)} ,
\\
x_3^{(k+1)} &= \frac{\omega}{7} \left( -11 + 6\,x_4^{(k)} +2\, x_1^{(k+1)} +4 \, x_2{(k+1)} \right) + \left( 1 - \omega \right) x_3^{(k)} ,
\\
x_4^{(k+1)} &= \frac{\omega}{6} \left( 9 - 2\,x_1^{(k+1)} +4 \, x_2{(k+1)} + 5\,x_3^{(k+1)} \right) + \left( 1 - \omega \right) x_4^{(k)} ,
\end{align*}
Now it is right time to organize our calculations in a loop.
Table of SOR iterates with ω = 0.8
n
0
1
2
3
4
5
x₁
0
\( \displaystyle \frac{}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{1}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₂
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₃
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₄
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
Table of SOR iterates with ω = 1.8
n
0
1
2
3
4
5
x₁
0
\( \displaystyle \frac{}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{1}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₂
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₃
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
x₄
0
\( \displaystyle -\frac{1}{} \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
\( \displaystyle \frac{}{} \approx \)
End of Example 9
Example 10:
We demonstrate performance of three iterating methods for solving the following system of equations:
\[
{\bf A}\,{\bf x} = {\bf b} ,
\tag{10.1}
\]
where
\[
{\bf A} = \begin{bmatrix} \phantom{-}4 & -1 & 2 \\ -2 & \phantom{-}4 & 5 \\ \phantom{-}1 & \phantom{-}2 & 5 \end{bmatrix}, \qquad {\bf b} = \begin{pmatrix} -2 \\ -4 \\ -5 \end{pmatrix} .
\]
This system has the true solution x = [1, 2, −2]. We start with Jacobi's scheme by extracting the diagonal matrix from A:
\[
\Lambda = \begin{bmatrix} 4 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 5 \end{bmatrix} , \qquad {\bf M} = \begin{bmatrix} \phantom{-}0 & -1 & 2 \\ -2 & \phantom{-}0 & 5 \\ \phantom{-}1 & \phantom{-}2 & 0 \end{bmatrix} .
\]
Although the given system has matrix A that is not diagonally dominant, Jacobi's scheme converges because eigenvalues of matrix M do not exceed 1.
A = {{4, -1, 2}, {-2, 4, 5}, {1, 2, 5}};
diag = DiagonalMatrix[{4, 4, 5}];
B = Inverse[diag] . (A - diag)
In this exercise, you will see that the order of the equations in a system can affect the performance of an iterative solution method. Upon swapping the order of equations in the previous exercise, we get
\[
\begin{split}
2\,x_1 + 6\, x_2 &= -16 ,
\\
5\, x_1 + x_2 &= 2 .
\end{split}
\]
Compute five iterations with both the Jacobi and Gauss--Seidel methods using zero as an initial guess.
Show that if a 2 x 2 coefficient matrix
\[
{\bf A} = \begin{bmatrix} a_{1,1} & a_{1,2} \\ a_{2,1} & a_{2,2} \end{bmatrix}
\]
is strictly diagonally dominant, then the eigenvalues of the iteration matrices B corresponding to the Jacobi and Gauss-Seidel methods are of magnitude < 1 (which guarantees convergence of iteration methods).
Apply Jacobi's iteration scheme to the system
\[
\begin{split}
x_1 + x_2 & = 0 ,
\\
x_1 - x_2 &= 2 .
\end{split}
\]
Do you think that the iteration method converges to the true solution x₁ = 1, x₂ = −1?
Apply Jacobi's iteration scheme for solving
\[
\begin{bmatrix} 2 & 1 \\ 2 & 3 \end{bmatrix} \begin{pmatrix} x_1 \\ x_2
\end{pmatrix} = \begin{pmatrix} \phantom{-}1 \\ -1 \end{pmatrix} .
\]
Recall that norm of the error is ∥e(k)∥ =
∥Be(k-1)∥ = |λmax| ∥e(k-1)∥, where λmax is the largest (in magnitude) eigenvalue of iteration matrix B. In general (for larger systems), this approximate relationship often is not the case for the first few iterations). Show that the error for this system tends to zero.
Use the SOR method to solve the system
\[
\begin{bmatrix} 2 & 1 \\ 3 & 4 \end{bmatrix} \begin{pmatrix} x_1 \\ x_2
\end{pmatrix} = \begin{pmatrix} \phantom{-}1 \\ -1 \end{pmatrix} .
\]
In the previous exercise, find the characteristic equation for the corresponding matrix B depending on parameter ω. What are its eigenvalues?
Solve the system of equations with a direct method:
\[
\begin{cases}
5 x - 2 y + 3 z &= 8,
\\
3 x + 8 y + 2 z &= 23,
\\
2 x + y + 6 z &= 2.
\end{cases}
\]
In order to understand differences and relative advantages among st the three iterative methods, perform 8 iterations using Jacobi's scheme, the Gauss--Seidel method, and SOR iteration procedure.
Let us consider the following system of linear equations:
\[
\begin{cases}
5 x - 2 y + z &= 2,
\\
9 x -3 y + 2 z &= 7,
\\
4 x + y + z &= 9.
\end{cases}
\]
Show that this system cannot
be tackled by means of Jacobi’s nor Gauss-Seidel’s methods. On the other hand, a
tuning of ω can lead to the detection of a good approximation of the solution.
Consider the system of equations
\[
\begin{bmatrix} 1 & -5 & 2 \\ 7 & 1 & 3 \\ 2 & 1 & 8 \end{bmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} -17 \\ 11 \\ -9 \end{pmatrix} .
\]
Show that the given system is not strictly diagonally dominant. Nonetheless, if we swap the first and second
equations, we obtain an equivalent system whose associated matrix is suitable for Jacobi's iteration scheme. Show this by performing 10 iterations.
Ciarlet, P.G., Introduction to Numerical Linear Algebra and Optimization.
Cambridge University Press, Cambridge, 1989.
Dongarra, J.J., Duff, I.S., Sorensen, D.C., and H. van der Vorst, H., Numerical
Linear Algebra for High-Performance Computers. SIAM, 1998.
Gauss, Carl, (1903), Werke (in German), vol. 9, Göttingen: Köninglichen Gesellschaft der Wissenschaften. direct link
Golub, Gene H.; Van Loan, Charles F. (1996), Matrix Computations (3rd ed.), Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9.
Kahan, W., Gauss–Seidel methods of solving large systems of linear equations, Doctoral Thesis, University of Toronto,
Toronto, Canada, 1958.
Seidel, Ludwig (1874). "Über ein Verfahren, die Gleichungen, auf welche die Methode der kleinsten Quadrate führt, sowie lineäre Gleichungen überhaupt, durch successive Annäherung aufzulösen" [On a process for solving by successive approximation the equations to which the method of least squares leads as well as linear equations generally]. Abhandlungen der Mathematisch-Physikalischen Klasse der Königlich Bayerischen Akademie der Wissenschaften (in German). 11 (3): 81–108.
Strang, G., Introduction to Linear Algebra, Wellesley-Cambridge Press, 5-th edition, 2016.
