Gauss–Seidel method

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 and Philipp Ludwig von Seidel, 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.

Description

The Gauss–Seidel method is an iterative technique for solving a square system of linear equations with unknown :
$A\mathbf x = \mathbf b .$

It is defined by the iteration
$L_* \mathbf^ = \mathbf - U \mathbf^,$
where $\mathbf^$ is the -th approximation or iteration of $\mathbf,\,\mathbf^$ is the next or -th iteration of $\mathbf$, and the matrix is decomposed into a lower triangular component $L_*$, and a strictly upper triangular component $U$; that is, $A = L_* + U$.

In more detail, write out , and in their components:

$A=\begin a_ & a_ & \cdots & a_ \\ a_ & a_ & \cdots & a_ \\ \vdots & \vdots & \ddots & \vdots \\a_ & a_ & \cdots & a_ \end, \qquad \mathbf = \begin x_ \\ x_2 \\ \vdots \\ x_n \end , \qquad \mathbf = \begin b_ \\ b_2 \\ \vdots \\ b_n \end.$

Then the decomposition of into its lower triangular component and its strictly upper triangular component is given by:

$A = \underbrace_ + \underbrace_ .$

The system of linear equations may be rewritten as:
:$L_* \mathbf = \mathbf - U \mathbf$

The Gauss–Seidel method now solves the left hand side of this expression for x, using previous value for x on the right hand side. Analytically, this may be written as:
$\mathbf^ = L_*^ \left(\mathbf - U \mathbf^\right).$

However, by taking advantage of the triangular form of $L_*$, the elements of can be computed sequentially using forward substitution:
$x^_i = \frac \left(b_i - \sum_^a_x^_j - \sum_^a_x^_j \right),\quad i=1,2,\dots,n.$

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 extremely similar to that of the Jacobi method.

The computation of uses the elements of that have already been computed, and only the elements of that have not been computed in the -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 SOR (successive over-relaxation) with $\omega=1$.

Convergence

The convergence properties of the Gauss–Seidel method are dependent on the matrix ''A''. Namely, the procedure is known to converge if either:
* is symmetric positive-definite, or
* is strictly or irreducibly diagonally dominant.

The Gauss–Seidel method sometimes converges even if these conditions are not satisfied.

Algorithm

Since elements can be overwritten as they are computed in this algorithm, only one storage vector is needed, and vector indexing is omitted. The algorithm goes as follows:

Inputs: ,

repeat until convergence
for from 1 until do

for from 1 until do
if ≠ then

end if
end (-loop)

end (-loop)
check if convergence is reached
end (repeat)

Examples

An example for the matrix version

A linear system shown as $A \mathbf = \mathbf$ is given by:
$A= \begin 16 & 3 \\ 7 & -11 \\ \end \quad \text \quad b= \begin 11 \\ 13 \end.$

We want to use the equation
$\mathbf^ = L_*^ (\mathbf - U \mathbf^)$
in the form
$\mathbf^ = T \mathbf^ + C$
where:
:$T = - L_*^ U \quad \text \quad C = L_*^ \mathbf.$

We must decompose $A$ into the sum of a lower triangular component $L_*$ and a strict upper triangular component $U$:
$L_*= \begin 16 & 0 \\ 7 & -11 \\ \end \quad \text \quad U = \begin 0 & 3 \\ 0 & 0 \end.$

The inverse of $L_*$ is:
$L_*^ = \begin 16 & 0 \\ 7 & -11 \end^ = \begin 0.0625 & 0.0000 \\ 0.0398 & -0.0909 \\ \end.$

Now we can find:
\begin
T &= -
\begin
0.0625 & 0.0000 \\
0.0398 & -0.0909
\end
\begin
0 & 3 \\
0 & 0
\end
=
\begin
0.000 & -0.1875 \\
0.000 & -0.1194
\end,
\\ ex C &=
\begin
0.0625 & 0.0000 \\
0.0398 & -0.0909
\end
\begin
11 \\
13
\end
=
\begin
0.6875 \\
-0.7439
\end.
\end

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

First of all, we have to choose $\mathbf^$: we can only guess. The better the guess, the quicker the algorithm will perform.

We choose a starting point:
$x^ = \begin 1.0 \\ 1.0 \end.$

We can then calculate:
\begin
x^ &=
\begin
0.000 & -0.1875 \\
0.000 & -0.1193
\end
\begin
1.0 \\
1.0
\end
+
\begin
0.6875 \\
-0.7443
\end
=
\begin
0.5000 \\
-0.8636
\end.
\\ ex x^ &=
\begin
0.000 & -0.1875 \\
0.000 & -0.1193
\end
\begin
0.5000 \\
-0.8636
\end
+
\begin
0.6875 \\
-0.7443
\end
=
\begin
0.8494 \\
-0.6413
\end.
\\ ex x^ &=
\begin
0.000 & -0.1875 \\
0.000 & -0.1193
\end
\begin
0.8494 \\
-0.6413 \\
\end
+
\begin
0.6875 \\
-0.7443
\end
=
\begin
0.8077 \\
-0.6678
\end.
\\ ex x^ &=
\begin
0.000 & -0.1875 \\
0.000 & -0.1193
\end
\begin
0.8077 \\
-0.6678
\end
+
\begin
0.6875 \\
-0.7443
\end
=
\begin
0.8127 \\
-0.6646
\end.
\\ ex x^ &=
\begin
0.000 & -0.1875 \\
0.000 & -0.1193
\end
\begin
0.8127 \\
-0.6646
\end
+
\begin
0.6875 \\
-0.7443
\end
=
\begin
0.8121 \\
-0.6650
\end.
\\ ex x^ &=
\begin
0.000 & -0.1875 \\
0.000 & -0.1193
\end
\begin
0.8121 \\
-0.6650
\end
+
\begin
0.6875 \\
-0.7443
\end
=
\begin
0.8122 \\
-0.6650
\end.
\\ ex x^ &=
\begin
0.000 & -0.1875 \\
0.000 & -0.1193
\end
\begin
0.8122 \\
-0.6650
\end
+
\begin
0.6875 \\
-0.7443
\end
=
\begin
0.8122 \\
-0.6650
\end.
\end

As expected, the algorithm converges to the exact solution:
$\mathbf = A^ \mathbf \approx \begin 0.8122\\ -0.6650 \end.$

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

Another example for the matrix version

Another linear system shown as $A \mathbf = \mathbf$ is given by:

$A= \begin 2 & 3 \\ 5 & 7 \\ \end \quad \text \quad b = \begin 11 \\ 13 \\ \end.$

We want to use the equation
$\mathbf^ = L_*^ (\mathbf - U \mathbf^)$
in the form
$\mathbf^ = T \mathbf^ + C$
where:
:$T = - L_*^ U \quad \text \quad C = L_*^ \mathbf.$

We must decompose $A$ into the sum of a lower triangular component $L_*$ and a strict upper triangular component $U$:
$L_*= \begin 2 & 0 \\ 5 & 7 \\ \end \quad \text \quad U = \begin 0 & 3 \\ 0 & 0 \\ \end.$

The inverse of $L_*$ is:
$L_*^ = \begin 2 & 0 \\ 5 & 7 \\ \end^ = \begin 0.500 & 0.000 \\ -0.357 & 0.143 \\ \end .$

Now we can find:
\begin
T &= -
\begin
0.500 & 0.000 \\
-0.357 & 0.143 \\
\end
\begin
0 & 3 \\
0 & 0 \\
\end
=
\begin
0.000 & -1.500 \\
0.000 & 1.071 \\
\end,
\\ ex C &=
\begin
0.500 & 0.000 \\
-0.357 & 0.143 \\
\end
\begin
11 \\
13 \\
\end
=
\begin
5.500 \\
-2.071 \\
\end.
\end

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

First of all, we have to choose $\mathbf^$: we can only guess. The better the guess, the quicker will perform the algorithm.

We suppose:
$x^ = \begin 1.1 \\ 2.3 \end.$

We can then calculate:
\begin
x^ &=
\begin
0 & -1.500 \\
0 & 1.071 \\
\end
\begin
1.1 \\
2.3 \\
\end
+
\begin
5.500 \\
-2.071 \\
\end
=
\begin
2.050 \\
0.393 \\
\end.
\\ ex x^ &=
\begin
0 & -1.500 \\
0 & 1.071 \\
\end
\begin
2.050 \\
0.393 \\
\end
+
\begin
5.500 \\
-2.071 \\
\end
=
\begin
4.911 \\
-1.651
\end.
\\ ex x^ &= \cdots.
\end

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
$\mathbf = A^ \mathbf = \begin -38\\ 29 \end$
is not guaranteed and, in this case, will not occur.

An example for the equation version

Suppose given equations where ''x''_''n'' are vectors of these equations and starting point ''x''₀.
From the first equation solve for ''x''₁ in terms of $x_, x_, \dots, x_n.$ For the next equations substitute the previous values of ''x''s.

To make it clear consider an example.
$\begin 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$

Solving for $x_1, x_2, x_3$ and $x_4$ gives:
$\begin 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$

Suppose we choose as the initial approximation, then the first approximate solution is given by
$\begin 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$

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.

The exact solution of the system is .

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., -1., 2., 0.
1., 11., -1., 3.
., -1., 10., -1.
., 3., -1., 8.)
# initialize the RHS vector
b = np.array( .0, 25.0, -11.0, 15.0

print("System of equations:")
for i in range(A.shape :
row = "*x"_for_j_in_range(A.shape[1.html" ;"title=".html" ;"title=""*x" for j in range(A.shape[1">"*x" for j in range(A.shape[1">.html" ;"title=""*x" for j in range(A.shape[1">"*x" for j in range(A.shape[1 print("[] = []".format(" + ".join(row), b[i]))

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

print(f"Solution: ")
error = np.dot(A, x) - b
print(f"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= 15Iteration 1: 0. 0. 0. 0.Iteration 2: 0.6 2.32727273 -0.98727273 0.87886364Iteration 3: 1.03018182 2.03693802 -1.0144562 0.98434122Iteration 4: 1.00658504 2.00355502 -1.00252738 0.99835095Iteration 5: 1.00086098 2.00029825 -1.00030728 0.99984975Iteration 6: 1.00009128 2.00002134 -1.00003115 0.9999881 Iteration 7: 1.00000836 2.00000117 -1.00000275 0.99999922Iteration 8: 1.00000067 2.00000002 -1.00000021 0.99999996Iteration 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
$x^_i = \frac \left(b_i - \sum_a_x^_j - \sum_a_x^_j \right),\quad \begin i=1,2,\ldots,n \\ k=0,1,2,\ldots \end$

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

See also

* Gaussian belief propagation
* Iterative method. Linear systems
*Kaczmarz method (a "row-oriented" method, whereas Gauss-Seidel is "column-oriented." See, for example
this paper
)
*Matrix splitting
*Richardson iteration

Notes

References

* .
* .
*

External links

*
Gauss–Seidel from www.math-linux.com
From Holistic Numerical Methods Institute
BicksonMatlab code

{{DEFAULTSORT:Gauss-Seidel Method
Numerical linear algebra
Articles with example pseudocode
Relaxation (iterative methods)
Articles with example Python (programming language) code
Articles with example MATLAB/Octave code