Gauss--Seidel method

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

input A and b
n 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,j xj
         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.

If you use column vector b instead of row vector, LinearSolve command should be applied as follows.

Jacobi's algorithm requires extraction of the diagonal matrix A = λ + M, where \[ {\bf A} = \begin{bmatrix} \phantom{-}4 & 2 & -2 \\ -2 & 5 & \phantom{-}2 \\ \phantom{-}4 & 3 & -2 \end{bmatrix} , \qquad \Lambda = \begin{bmatrix} 4 & 0 & \phantom{-}0 \\ 0 & 5 & \phantom{-}0 \\ 0 & 0 & -2 \end{bmatrix} , \quad {\bf M} = \begin{bmatrix} \phantom{-}0 & 2 & -2 \\ -2 & 0 & \phantom{-}2 \\ \phantom{-}4 & 3 & \phantom{-}0 \end{bmatrix} . \]

We rewrite the given system of equations in equivalent form \[ \Lambda\,{\bf x} = {\bf b} - {\bf M}\,{\bf x} \qquad \iff \qquad {\bf x} = {\bf d} + {\bf B}\,{\bf x} , \tag{5.2} \] where \[ {\bf d} = \Lambda^{-1} {\bf b} = \begin{pmatrix} 3 \\ \frac{4}{5} \\ -7 \end{pmatrix} , \qquad {\bf B} = -\Lambda^{-1} {\bf M} = \begin{bmatrix} 0 & - \frac{1}{2} & \frac{1}{2} \\ \frac{2}{5} & 0 & - \frac{2}{5} \\ 2 & \frac{3}{2} & 0 \end{bmatrix} . \] Although matrix B is not diagonally dominant, its eigenvalues have magnitude less than 1. So we expect that Jacobi's iteration scheme converges.

Upon choosing the initial guess to be zero, we obtain the first approximation to be \[ {\bf x}^{(1)} = {\bf d} = \begin{pmatrix} 3 \\ \frac{4}{5} \\ -7 \end{pmatrix} . \] The second iteration step yields \[ {\bf x}^{(2)} = {\bf d} + {\bf B}\,{\bf x}^{(1)} = \begin{pmatrix} - \frac{9}{10} \\ \frac{24}{5} \\ \frac{1}{5} \end{pmatrix} = \begin{pmatrix} -0.9 \\ 4.8 \\ 0.2 \end{pmatrix} . \]

With third iteration, we have \[ {\bf x}^{(3)} = {\bf d} + {\bf B}\,{\bf x}^{(2)} = {\bf d} + {\bf B}\,{\bf d} + {\bf B}^2 {\bf d} = \begin{pmatrix} \frac{7}{10} \\ \frac{9}{25} \\ -\frac{8}{5} \end{pmatrix} = \begin{pmatrix} 0.7 \\ 0.36 \\ -1.6 \end{pmatrix} . \] Application of jacMeth code from Example 1 yields

The errors of Jacobi's iterations at 50-th run are \[ \left\vert x_1^{(50)} -1 \right\vert = 0.039291 , \qquad \left\vert x_2^{(50)} -2 \right\vert = 0.0101124 , \qquad \left\vert x_3^{(50)} +2 \right\vert = 0.069383 . \]

We repeat calculations with 50 runs. It took less than a second to execute the code.

Now we apply the Gauss--Seidel method with the same initial choice---the origin. At the first step, we have \begin{align*} x_1^{(1)} &= 3 - \frac{1}{2}\, x_2^{(0)} + \frac{1}{2}\, x_3^{(0)} = 3, \\ x_2^{(1)} &= \frac{4}{5} + \frac{2}{5}\, x_1^{(1)} - \frac{2}{5}\, x_3^{(0)} = 2 \\ x_3^{(1)} &= -7 + 2\,x_1^{(1)} + \frac{3}{2}\,x_2^{(1)} = 2 . \end{align*} Next iteration gives \begin{align*} x_1^{(2)} &= 3 - \frac{1}{2}\, x_2^{(1)} + \frac{1}{2}\, x_3^{(1)} = 3, \\ x_2^{(2)} &= \frac{4}{5} + \frac{2}{5}\, x_1^{(2)} - \frac{2}{5}\, x_3^{(1)} = \frac{6}{5} , \\ x_3^{(2)} &= -7 + 2\,x_1^{(2)} + \frac{3}{2}\,x_2^{(2)} = \frac{4}{5} \end{align*} We manually repeat a few iterations before application of looping procedure. \begin{align*} x_1^{(3)} &= 3 - \frac{1}{2}\, x_2^{(2)} + \frac{1}{2}\, x_3^{(2)} = \frac{14}{5} , \\ x_2^{(3)} &= \frac{4}{5} + \frac{2}{5}\, x_1^{(3)} - \frac{2}{5}\, x_3^{(2)} = \frac{8}{5} , \\ x_3^{(3)} &= -7 + 2\,x_1^{(3)} + \frac{3}{2}\,x_2^{(3)} = 1 , \end{align*} and \begin{align*} x_1^{(4)} &= 3 - \frac{1}{2}\, x_2^{(3)} + \frac{1}{2}\, x_3^{(3)} = \frac{27}{10} , \\ x_2^{(4)} &= \frac{4}{5} + \frac{2}{5}\, x_1^{(3)} - \frac{2}{5}\, x_3^{(3)} = \frac{37}{25} , \\ x_3^{(4)} &= -7 + 2\,x_1^{(4)} + \frac{3}{2}\,x_2^{(4)} = \frac{31}{50} ; \end{align*} \begin{align*} x_1^{(5)} &= 3 - \frac{1}{2}\, x_2^{(4)} + \frac{1}{2}\, x_3^{(4)} = \frac{257}{100} = 2.57 , \\ x_2^{(5)} &= \frac{4}{5} + \frac{2}{5}\, x_1^{(5)} - \frac{2}{5}\, x_3^{(4)} = \frac{79}{50} = 1.58 , \\ x_3^{(5)} &= -7 + 2\,x_1^{(5)} + \frac{3}{2}\,x_2^{(5)} = \frac{51}{100} = 0.51 . \end{align*} Sequential iterates are generated with the following Mathematica code:

We repeat GS-iterative procedure with 50 runs. It tAkes less than a second to execute 50 runs.

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 point instead of exact operations with rational numbers does not effect the time of execution.

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.

We repeat calculations with floating point arithmetic for 50 runs.

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.

Example 7: 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{6.2} \] is singular of rank 3. We split the matrix into lower triangular and upper triangular parts: \[ {\bf A} = {\bf L} + {\bf U} , \] where \[ {\bf L} = \begin{bmatrix} \phantom{-}9 & 0 & \phantom{-}0 & 0 \\ \phantom{-}1 & \phantom{-}8 & \phantom{-}0 & 0 \\ -2 & -4 & \phantom{-}7 & 0 \\ \phantom{-}2 & -4 & -5 & 6 \end{bmatrix} , \qquad {\bf U} = \begin{bmatrix} 0 & -8 & \phantom{-}5 & -4 \\ 0 & \phantom{-}0 & -5 & \phantom{-}3 \\ 0 & \phantom{-}0 & \phantom{-}0 & -6 \\ 0 & \phantom{-}0 & \phantom{-}0 & \phantom{-}0 \end{bmatrix} . \] This allows us to rewrite equation (6.2) in fixed point form: \[ {\bf L}\,{\bf x} = {\bf b} - {\bf U}\,{\bf x} \qquad \iff \qquad {\bf x} = {\bf L}^{-1} {\bf b} - {\bf L}^{-1} {\bf U}\,{\bf x} . \tag{6.3} \] Upon introducing auxiliary matrices, \[ {\bf d} = {\bf L}^{-1} {\bf b} = \begin{pmatrix} 1 \\ \frac{3}{4} \\ -\frac{6}{7} \\ \frac{20}{21} \end{pmatrix} , \qquad {\bf B} = - {\bf L}^{-1} {\bf U} = \begin{bmatrix} 0 & \frac{8}{9} & - \frac{5}{9} & \frac{4}{9} \\ 0 & - \frac{1}{9} &\frac{25}{36} & -\frac{31}{72} \\ 0 & \frac{4}{21} & \frac{5}{21} & \frac{31}{42} \\ 0 & - \frac{40}{189} & \frac{160}{189} & \frac{34}{189} \end{bmatrix} , \] we rewrite Eq.(6.3) in compact form: \[ {\bf x} = {\bf d} + {\bf B}\, {\bf x} . \]

Eq.(6.3) is used to generate Gauss-Seidel iteration: \[ {\bf x}^{(k+1)} = {\bf L}^{-1} {\bf b} - {\bf L}^{-1} {\bf U}\,{\bf x}^{(k)} \qquad \iff \qquad {\bf x}^{(k+1)} = {\bf d} + {\bf B}\,{\bf x}^{(k)} , \qquad k=0,1,2,\ldots . \tag{6.4} \]

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.

Table of Jacobi iterates

n	0	1	2	3	4	5	10
x₁	0	1	\( \displaystyle \frac{485}{189} \approx 2.56614 \)	0.925058	2.31112	1.08992	1.94983
x₂	0	\( \displaystyle -\frac{3}{4} \)	\( \displaystyle -\frac{64}{189} \approx -0.338624 \)	0.534972	-0.243334	0.434793	-0.0439163
x₃	0	\( \displaystyle -\frac{6}{7} \)	\( \displaystyle -\frac{95}{441} \approx -0.21542 \)	-0.291808	-0.291808	-0.733657	-0.422114
x₄	0	\( \displaystyle \frac{20}{21} \)	\( \displaystyle \frac{950}{3969} \approx .239355 \)	0.88474	0.324231	0.815174	0.469016

Running GS-iterations 50 times yields

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.

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.

1. Anton, Howard (2005), Elementary Linear Algebra (Applications Version) (9th ed.), Wiley International
2. Beezer, R., A First Course in Linear Algebra, 2015.
3. Beezer, R., A Second Course in Linear Algebra, 2013.