Matrix-based formula

The solution is obtained iteratively via

$${\displaystyle L_{*}\mathbf {x} ^{(k+1)}=\mathbf {b} -U\mathbf {x} ^{(k)},}$$

where the matrix ${\displaystyle A}$ is decomposed into a lower triangular component ${\displaystyle L_{*}}$, and a strictly upper triangular component ${\displaystyle U}$ such that ${\displaystyle A=L_{*}+U}$.^[4] More specifically, the decomposition of ${\displaystyle A}$ into ${\displaystyle L_{*}}$ and ${\displaystyle U}$ is given by:

$${\displaystyle A=\underbrace {\begin{bmatrix}a_{11}&0&\cdots &0\\a_{21}&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}} _{\textstyle L_{*}}+\underbrace {\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\0&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &0\end{bmatrix}} _{\textstyle U}.}$$

Why the matrix-based formula works

The system of linear equations may be rewritten as:

${\displaystyle {\begin{alignedat}{1}A\mathbf {x} &=\mathbf {b} \\(L_{*}+U)\mathbf {x} &=\mathbf {b} \\L_{*}\mathbf {x} +U\mathbf {x} &=\mathbf {b} \\L_{*}\mathbf {x} &=\mathbf {b} -U\mathbf {x} \end{alignedat}}}$

The Gauss–Seidel method now solves the left hand side of this expression for ${\displaystyle \mathbf {x} }$, using previous value for ${\displaystyle \mathbf {x} }$ on the right hand side. Analytically, this may be written as:

$${\displaystyle \mathbf {x} ^{(k+1)}=L_{*}^{-1}\left(\mathbf {b} -U\mathbf {x} ^{(k)}\right).}$$

Element-based formula

However, by taking advantage of the triangular form of ${\displaystyle L_{*}}$, the elements of ${\displaystyle \mathbf {x} ^{(k+1)}}$ can be computed sequentially for each row ${\displaystyle i}$ using forward substitution:^[5]

$${\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum _{j=i+1}^{n}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\dots ,n.}$$

Notice that the formula uses two summations per iteration which can be expressed as one summation ${\displaystyle \sum _{j\neq i}a_{ij}x_{j}}$ that uses the most recently calculated iteration of ${\displaystyle x_{j}}$. The procedure is generally continued until the changes made by an iteration are below some tolerance, such as a sufficiently small residual.

Discussion

The element-wise formula for the Gauss–Seidel method is similar to that of the Jacobi method.

The computation of ${\displaystyle \mathbf {x} ^{(k+1)}}$ uses the elements of ${\displaystyle \mathbf {x} ^{(k+1)}}$ that have already been computed, and only the elements of ${\displaystyle \mathbf {x} ^{(k)}}$ that have not been computed in the (k+1)-th iteration. This means that, unlike the Jacobi method, only one storage vector is required as elements can be overwritten as they are computed, which can be advantageous for very large problems.

However, unlike the Jacobi method, the computations for each element are generally much harder to implement in parallel, since they can have a very long critical path, and are thus most feasible for sparse matrices. Furthermore, the values at each iteration are dependent on the order of the original equations.

Gauss-Seidel is the same as successive over-relaxation with ${\displaystyle \omega =1}$.

An example for the matrix version

A linear system shown as ${\displaystyle A\mathbf {x} =\mathbf {b} }$ is given by:

$${\displaystyle A={\begin{bmatrix}16&3\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad b={\begin{bmatrix}11\\13\end{bmatrix}}.}$$

We want to use the equation

$${\displaystyle \mathbf {x} ^{(k+1)}=L_{*}^{-1}(\mathbf {b} -U\mathbf {x} ^{(k)})}$$

in the form

$${\displaystyle \mathbf {x} ^{(k+1)}=T\mathbf {x} ^{(k)}+C}$$

where:

${\displaystyle T=-L_{*}^{-1}U\quad {\text{and}}\quad C=L_{*}^{-1}\mathbf {b} .}$

We must decompose ${\displaystyle A}$ into the sum of a lower triangular component ${\displaystyle L_{*}}$ and a strict upper triangular component ${\displaystyle U}$:

$${\displaystyle L_{*}={\begin{bmatrix}16&0\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad U={\begin{bmatrix}0&3\\0&0\end{bmatrix}}.}$$

The inverse of ${\displaystyle L_{*}}$ is:

$${\displaystyle L_{*}^{-1}={\begin{bmatrix}16&0\\7&-11\end{bmatrix}}^{-1}={\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\\\end{bmatrix}}.}$$

Now we can find:

$${\displaystyle {\begin{aligned}T&=-{\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\end{bmatrix}}{\begin{bmatrix}0&3\\0&0\end{bmatrix}}={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1194\end{bmatrix}},\\[1ex]C&={\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\end{bmatrix}}{\begin{bmatrix}11\\13\end{bmatrix}}={\begin{bmatrix}0.6875\\-0.7439\end{bmatrix}}.\end{aligned}}}$$

Now we have ${\displaystyle T}$ and ${\displaystyle C}$ and we can use them to obtain the vectors ${\displaystyle \mathbf {x} }$ iteratively.

First of all, we have to choose ${\displaystyle \mathbf {x} ^{(0)}}$: we can only guess. The better the guess, the quicker the algorithm will perform.

We choose a starting point:

$${\displaystyle x^{(0)}={\begin{bmatrix}1.0\\1.0\end{bmatrix}}.}$$

We can then calculate:

$${\displaystyle {\begin{aligned}x^{(1)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}1.0\\1.0\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.5000\\-0.8636\end{bmatrix}}.\\[1ex]x^{(2)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.5000\\-0.8636\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8494\\-0.6413\end{bmatrix}}.\\[1ex]x^{(3)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8494\\-0.6413\\\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8077\\-0.6678\end{bmatrix}}.\\[1ex]x^{(4)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8077\\-0.6678\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8127\\-0.6646\end{bmatrix}}.\\[1ex]x^{(5)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8127\\-0.6646\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8121\\-0.6650\end{bmatrix}}.\\[1ex]x^{(6)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8121\\-0.6650\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.\\[1ex]x^{(7)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.\end{aligned}}}$$

As expected, the algorithm converges to the exact solution:

$${\displaystyle \mathbf {x} =A^{-1}\mathbf {b} \approx {\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.}$$

In fact, the matrix A is strictly diagonally dominant (but not positive definite).

Another example for the matrix version

Another linear system shown as ${\displaystyle A\mathbf {x} =\mathbf {b} }$ is given by:

$${\displaystyle A={\begin{bmatrix}2&3\\5&7\\\end{bmatrix}}\quad {\text{and}}\quad b={\begin{bmatrix}11\\13\\\end{bmatrix}}.}$$

We want to use the equation

$${\displaystyle \mathbf {x} ^{(k+1)}=L_{*}^{-1}(\mathbf {b} -U\mathbf {x} ^{(k)})}$$

in the form

$${\displaystyle \mathbf {x} ^{(k+1)}=T\mathbf {x} ^{(k)}+C}$$

where:

${\displaystyle T=-L_{*}^{-1}U\quad {\text{and}}\quad C=L_{*}^{-1}\mathbf {b} .}$

We must decompose ${\displaystyle A}$ into the sum of a lower triangular component ${\displaystyle L_{*}}$ and a strict upper triangular component ${\displaystyle U}$:

$${\displaystyle L_{*}={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}\quad {\text{and}}\quad U={\begin{bmatrix}0&3\\0&0\\\end{bmatrix}}.}$$

The inverse of ${\displaystyle L_{*}}$ is:

$${\displaystyle L_{*}^{-1}={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}^{-1}={\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}.}$$

Now we can find:

$${\displaystyle {\begin{aligned}T&=-{\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}{\begin{bmatrix}0&3\\0&0\\\end{bmatrix}}={\begin{bmatrix}0.000&-1.500\\0.000&1.071\\\end{bmatrix}},\\[1ex]C&={\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}{\begin{bmatrix}11\\13\\\end{bmatrix}}={\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}.\end{aligned}}}$$

Now we have ${\displaystyle T}$ and ${\displaystyle C}$ and we can use them to obtain the vectors ${\displaystyle \mathbf {x} }$ iteratively.

First of all, we have to choose ${\displaystyle \mathbf {x} ^{(0)}}$: we can only guess. The better the guess, the quicker will perform the algorithm.

We suppose:

$${\displaystyle x^{(0)}={\begin{bmatrix}1.1\\2.3\end{bmatrix}}.}$$

We can then calculate:

$${\displaystyle {\begin{aligned}x^{(1)}&={\begin{bmatrix}0&-1.500\\0&1.071\\\end{bmatrix}}{\begin{bmatrix}1.1\\2.3\\\end{bmatrix}}+{\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}={\begin{bmatrix}2.050\\0.393\\\end{bmatrix}}.\\[1ex]x^{(2)}&={\begin{bmatrix}0&-1.500\\0&1.071\\\end{bmatrix}}{\begin{bmatrix}2.050\\0.393\\\end{bmatrix}}+{\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}={\begin{bmatrix}4.911\\-1.651\end{bmatrix}}.\\[1ex]x^{(3)}&=\cdots .\end{aligned}}}$$

If we test for convergence we'll find that the algorithm diverges. In fact, the matrix A is neither diagonally dominant nor positive definite. Then, convergence to the exact solution

$${\displaystyle \mathbf {x} =A^{-1}\mathbf {b} ={\begin{bmatrix}-38\\29\end{bmatrix}}}$$

is not guaranteed and, in this case, will not occur.

An example for the equation version

Suppose given k equations where x_n are vectors of these equations and starting point x₀. From the first equation solve for x₁ in terms of ${\displaystyle x_{n+1},x_{n+2},\dots ,x_{n}.}$ For the next equations substitute the previous values of xs.

To make it clear consider an example.

$${\displaystyle {\begin{array}{rrrrl}10x_{1}&-x_{2}&+2x_{3}&&=6,\\-x_{1}&+11x_{2}&-x_{3}&+3x_{4}&=25,\\2x_{1}&-x_{2}&+10x_{3}&-x_{4}&=-11,\\&3x_{2}&-x_{3}&+8x_{4}&=15.\end{array}}}$$

Solving for ${\displaystyle x_{1},x_{2},x_{3}}$ and ${\displaystyle x_{4}}$ gives:

$${\displaystyle {\begin{aligned}x_{1}&=x_{2}/10-x_{3}/5+3/5,\\x_{2}&=x_{1}/11+x_{3}/11-3x_{4}/11+25/11,\\x_{3}&=-x_{1}/5+x_{2}/10+x_{4}/10-11/10,\\x_{4}&=-3x_{2}/8+x_{3}/8+15/8.\end{aligned}}}$$

Suppose we choose (0, 0, 0, 0) as the initial approximation, then the first approximate solution is given by

$${\displaystyle {\begin{aligned}x_{1}&=3/5=0.6,\\x_{2}&=(3/5)/11+25/11=3/55+25/11=2.3272,\\x_{3}&=-(3/5)/5+(2.3272)/10-11/10=-3/25+0.23272-1.1=-0.9873,\\x_{4}&=-3(2.3272)/8+(-0.9873)/8+15/8=0.8789.\end{aligned}}}$$

Using the approximations obtained, the iterative procedure is repeated until the desired accuracy has been reached. The following are the approximated solutions after four iterations.

${\displaystyle x_{1}}$	${\displaystyle x_{2}}$	${\displaystyle x_{3}}$	${\displaystyle x_{4}}$
0.6	2.32727	−0.987273	0.878864
1.03018	2.03694	−1.01446	0.984341
1.00659	2.00356	−1.00253	0.998351
1.00086	2.0003	−1.00031	0.99985

The exact solution of the system is (1, 2, −1, 1).

An example using Python and NumPy

The following numerical procedure simply iterates to produce the solution vector.

import numpy as np

ITERATION_LIMIT = 1000

# initialize the matrix
A = np.array(
    [
        [10.0, -1.0, 2.0, 0.0],
        [-1.0, 11.0, -1.0, 3.0],
        [2.0, -1.0, 10.0, -1.0],
        [0.0, 3.0, -1.0, 8.0],
    ]
)
# initialize the RHS vector
b = np.array([6.0, 25.0, -11.0, 15.0])

print("System of equations:")
for i in range(A.shape[0]):
    row = [f"{A[i,j]:3g}*x{j+1}" for j in range(A.shape[1])]
    print("[{0}] = [{1:3g}]".format(" + ".join(row), b[i]))

x = np.zeros_like(b, np.float_)
for it_count in range(1, ITERATION_LIMIT):
    x_new = np.zeros_like(x, dtype=np.float_)
    print(f"Iteration {it_count}: {x}")
    for i in range(A.shape[0]):
        s1 = np.dot(A[i, :i], x_new[:i])
        s2 = np.dot(A[i, i + 1 :], x[i + 1 :])
        x_new[i] = (b[i] - s1 - s2) / A[i, i]
    if np.allclose(x, x_new, rtol=1e-8):
        break
    x = x_new

print(f"Solution: {x}")
error = np.dot(A, x) - b
print(f"Error: {error}")

Produces the output:

System of equations:
[ 10*x1 +  -1*x2 +   2*x3 +   0*x4] = [  6]
[ -1*x1 +  11*x2 +  -1*x3 +   3*x4] = [ 25]
[  2*x1 +  -1*x2 +  10*x3 +  -1*x4] = [-11]
[  0*x1 +   3*x2 +  -1*x3 +   8*x4] = [ 15]
Iteration 1: [ 0.  0.  0.  0.]
Iteration 2: [ 0.6         2.32727273 -0.98727273  0.87886364]
Iteration 3: [ 1.03018182  2.03693802 -1.0144562   0.98434122]
Iteration 4: [ 1.00658504  2.00355502 -1.00252738  0.99835095]
Iteration 5: [ 1.00086098  2.00029825 -1.00030728  0.99984975]
Iteration 6: [ 1.00009128  2.00002134 -1.00003115  0.9999881 ]
Iteration 7: [ 1.00000836  2.00000117 -1.00000275  0.99999922]
Iteration 8: [ 1.00000067  2.00000002 -1.00000021  0.99999996]
Iteration 9: [ 1.00000004  1.99999999 -1.00000001  1.        ]
Iteration 10: [ 1.  2. -1.  1.]
Solution: [ 1.  2. -1.  1.]
Error: [  2.06480930e-08  -1.25551054e-08   3.61417563e-11   0.00000000e+00]

Program to solve arbitrary no. of equations using Matlab

The following code uses the formula

$${\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j<i}a_{ij}x_{j}^{(k+1)}-\sum _{j>i}a_{ij}x_{j}^{(k)}\right),\quad {\begin{array}{l}i=1,2,\ldots ,n\\k=0,1,2,\ldots \end{array}}}$$

function x = gauss_seidel(A, b, x, iters)
    for i = 1:iters
        for j = 1:size(A,1)
            x(j) = (b(j) - sum(A(j,:)'.*x) + A(j,j)*x(j)) / A(j,j);
        end
    end
end