Gauss-Seidel method

With $A=D+L+U$, the Gauss-Seidel method solves $(D+L)x^{(k+1)}=b-U{x}^{(k)}$. It uses each newly computed coordinate immediately in the remaining triangular solve. It is successive over-relaxation with parameter one, and converges from every start when $A$ is a symmetric positive-definite matrix.

To solve a linear system $Ax=b$, split its matrix into the diagonal matrix $D$ and the strictly lower and upper parts $L,U$, retaining their actual signs. Every diagonal entry must be nonzero. A sweep updates coordinates in order:The Jacobi method uses the old vector throughout; the Gauss-Seidel method reuses each fresh coordinate. Thus the update order is part of the method.^[ref]

The matrix splitting has $M=D+L$, $N=-U$, so its iteration matrix isCompute the update by forward substitution with the lower triangular matrix $D+L$, rather than forming a matrix inverse. For a sparse matrix, a sweep only visits stored entries and can overwrite the working vector. Dependencies on fresh coordinates constrain parallel execution. This makes the Gauss-Seidel method an inexpensive solver component, including in a multigrid method.^[ref]

For a real symmetric positive-definite matrix $A=D+L+L^T$, let $e_k$ be the Gauss-Seidel method error and $d_k=e_k-e_{k+1}$. Its recurrence gives $A{e}_k=(D+L)d_k$. Expanding the two squared energy norms yields $2d_k^{T}A{e}_k-d_k^{T}A{d}_k=d_k^{T}D{d}_k$. Every diagonal entry of $D$ is positive, so the error energy strictly decreases unless $d_k=0$. In that case $A{e}_k=0$ and the exact solution has already been reached. This identity exhibits the method as a sequence of energy-reducing coordinate corrections.

Consider the original linear systemDirect substitution verifies the solution. The matrix is a symmetric positive-definite matrix: completing the square in its quadratic form gives $4x^2+2xy+3y^2=4(x+y/4{)}^2+11y^2/4>0$ for every nonzero vector. It is also a strictly diagonally dominant matrix.

Starting from $(x^{(0)},y^{(0)})=(0,0)$, update the first row and then the second:The first sweep gives $x^{(1)}=1/4$ and $y^{(1)}=(2-1/4)/3=7/12$. The next gives $x^{(2)}=(1-7/12)/4=5/48$ and $y^{(2)}=(2-5/48)/3=91/144$. The second row uses $1/4$ during the first sweep, even though the initial first coordinate was zero.

With a stale first coordinate, the Jacobi method would instead give $(1/4,2/3)$ after the first sweep. Both methods have the same fixed point, but changing when coordinates become available changes the entire sequence of vectors.

For the Gauss-Seidel iteration for a two-variable positive-definite system, define the forward error by $e_x^{(k)}=x^{(k)}-x_{\ast }$ and $e_y^{(k)}=y^{(k)}-y_{\ast }$. Subtract the two exact solution equations from the two update equations:The fresh first coordinate passes the old second-coordinate error through two successive factors. The iteration matrix is thereforeIts eigenvalues are $0$ and $1/12$, as seen from the determinant. The zero eigenvalue explains why the old first-coordinate error disappears completely: the first row overwrites it. After one sweep, every remaining forward error lies on the eigenvector line spanned by $(-3,1{)}^T$.

For this zero start, $e_y^{(0)}=-7/11$, so the exact vectors for $k\ge 1$ areThe Euclidean norm of the forward error decreases by exactly $1/12$ per sweep from $k=1$ onward. The first decrease is different because the initial error has not yet reached that eigenvector line. Each row update moves horizontally or vertically in the coordinate plane; sweep endpoints then approach the solution along one straight line.

The linear-system residual is an independent way to check the calculation. Completing a sweep satisfies the second row exactly, but changing $y$ disturbs the first row:Thus a zero residual in the last updated row alone does not certify a solution; inspect the complete linear-system residual.

For the same matrix, the stale schedule of the Jacobi method gives the iteration matrixHere its eigenvalues are $\pm 1/\sqrt{12}$, whereas the fresh schedule gives $\rho (H)=\rho (B{)}^2$. This comparison belongs to this particular linear system; it is not a universal relation between the two methods' spectral radii.

A strictly diagonally dominant matrix or a real symmetric positive-definite matrix guarantees convergence of the Gauss-Seidel method from every initial vector. These are sufficient conditions.^[ref]

For strict diagonal dominance, putInduct through a sweep. Old coordinates have absolute error at most $E$; previously updated coordinates have absolute error at most $qE\le E$. The component update and the triangle inequality therefore bound the new error by $qE$. This proves $\Vert e^{(k+1)}{\Vert }_{\infty }\le q\Vert e^{(k)}{\Vert }_{\infty }$, a maximum norm contraction.

For a symmetric positive-definite matrix, write $A=D+L+L^T$. The Householder-John theorem applies to $M=D+L$ and $N=-L^T$, because $M^T+N=D$ is a positive-definite matrix. Hence the iteration matrix has spectral radius less than one, which gives the convergence criterion for a stationary iteration.

There is also a geometric reason fresh updates help: the Gauss-Seidel method is exact coordinate descent on the convex quadratic functionHolding the other coordinates fixed, minimizing in coordinate $i$ gives the corresponding row update. A change ${\delta }_i$ lowers $F$ by exactly $a_{ii}{\delta }_i^2/2$. Each fresh coordinate is therefore the minimizer for the current vector, and a complete sweep lowers the quadratic form error unless every coordinate already solves its row. This explains the method's beauty without implying monotone decrease in every vector norm.

An invertible matrix and nonzero diagonal entries make the update well defined, but do not guarantee convergence. For the original linear systemthe Gauss-Seidel method updates $x^{(k+1)}=1-2y^{(k)}$ and $y^{(k+1)}=1-2x^{(k+1)}=-1+4y^{(k)}$. Starting at zero gives $(1,-1)$, then $(3,-5)$. The second-coordinate forward error is $e_y^{(k)}=-4^k/3$, and the iteration matrix $H=\left(\begin{array}{cc}0& -2\\ 0& 4\end{array}\right)$ has spectral radius $4$. The linear system has a unique solution while this particular stationary iterative method for a linear system diverges.