Chapter 90. Gauss-Seidel Method

The Gauss-Seidel method is a stationary iterative method for solving a linear system

It is closely related to the Jacobi method. Both methods solve each equation for one unknown and repeat the process. The difference is that Gauss-Seidel uses the newest available values immediately. When the method computes $x_i^{(k+1)}$, it uses already updated values $x_1^{(k+1)}, \ldots, x_{i-1}^{(k+1)}$ and old values $x_{i+1}^{(k)}, \ldots, x_n^{(k)}$.

This small change often improves convergence. It also makes the method less naturally parallel than Jacobi.

Gauss-Seidel is guaranteed to converge for important matrix classes, including strictly diagonally dominant matrices and symmetric positive definite matrices. More generally, convergence is controlled by the spectral radius of its iteration matrix.

90.1 The Basic System

Let

The system

means

Assume each diagonal entry is nonzero:

Solving the $i$-th equation for $x_i$ gives

The Gauss-Seidel method uses this identity as an update rule.

90.2 Gauss-Seidel Update Formula

Starting with an initial guess

the Gauss-Seidel method defines

The first sum uses new values from the current iteration.

The second sum uses old values from the previous iteration.

This is the defining feature of Gauss-Seidel.

90.3 Comparison with Jacobi

Jacobi computes all new components from old components:

Gauss-Seidel computes components sequentially and immediately reuses updated values:

Gauss-Seidel may be viewed as a Jacobi-like method with immediate feedback.

90.4 Example

Consider

Solving each equation for its diagonal variable gives

The Gauss-Seidel iteration is

Start with

First update:

Then use this new value immediately:

Thus

Second update:

The exact solution is

The iterates move toward the exact solution.

90.5 Matrix Splitting

Write

where:

Then

becomes

Move the strictly upper triangular part to the right:

The Gauss-Seidel iteration is

Since $D+L$ is lower triangular, each iteration requires one triangular solve.

90.6 Fixed-Point Form

The iteration may be written as

Thus the iteration matrix is

The fixed-point form is

where

If the sequence converges to a vector $x$, then

Multiplying by $D+L$, we recover

Hence the limit solves the original system.

90.7 Error Recurrence

Let $x$ be the exact solution. Since

and

subtracting gives

With

we obtain

Therefore,

Convergence depends on whether powers of $G_{GS}$ approach zero.

90.8 Convergence Criterion

Gauss-Seidel converges for every initial guess if and only if

where $\rho(G_{GS})$ is the spectral radius of the iteration matrix. For a general splitting $A=M-N$, the corresponding stationary iteration converges when the spectral radius of $M^{-1}N$ is less than one.

For Gauss-Seidel,

so

Thus the convergence condition becomes

90.9 Strict Diagonal Dominance

A matrix is strictly row diagonally dominant if

$$ |a_{ii}|

\sum_{j\ne i}|a_{ij}| $$

for every row $i$.

If $A$ is strictly diagonally dominant, then the Gauss-Seidel method converges. This is a standard sufficient condition.

Strict diagonal dominance means that each equation is controlled mainly by its own variable. The off-diagonal coupling is not strong enough to prevent contraction of the iteration error.

This condition is sufficient, but not necessary. Gauss-Seidel may converge even when strict diagonal dominance fails.

90.10 Symmetric Positive Definite Matrices

Another important convergence result applies to symmetric positive definite matrices.

If $A$ is symmetric positive definite, then Gauss-Seidel converges.

This class appears frequently in applications, including:

The SPD case is especially important because it also supports conjugate gradient methods.

90.11 Geometric Interpretation

Gauss-Seidel may be interpreted as solving one coordinate direction at a time.

At each step, the method adjusts $x_i$ so that the $i$-th equation is satisfied using the newest known values.

After updating $x_1$, the first equation is satisfied relative to the current state. After updating $x_2$, the second equation is satisfied relative to the updated $x_1$, and so on.

The method sweeps through the equations. Each sweep reduces error when the coupling structure is favorable.

90.12 Residual Form

The residual at step $k$ is

Gauss-Seidel can be interpreted as applying a lower-triangular correction.

From

subtract

from both sides:

Thus

Therefore,

Gauss-Seidel updates the current approximation by applying a lower-triangular solve to the residual.

90.13 Algorithm

Input:

For

update each component in order:

Here the same vector $x$ is overwritten in place.

After a full sweep, compute the residual

Stop if

Return the current $x$.

90.14 Pseudocode

gauss_seidel(A, b, x, tol, max_iter):
    n = length(b)

    for k = 0 to max_iter - 1:
        for i = 1 to n:
            s1 = 0
            s2 = 0

            for j = 1 to i - 1:
                s1 = s1 + A[i,j] * x[j]

            for j = i + 1 to n:
                s2 = s2 + A[i,j] * x[j]

            x[i] = (b[i] - s1 - s2) / A[i,i]

        r = b - A * x

        if norm(r) / norm(b) <= tol:
            return x

    return x

Unlike Jacobi, the method updates $x$ in place.

This is why Gauss-Seidel needs less vector storage but has less natural parallelism.

90.15 Storage Requirements

Gauss-Seidel needs:

Jacobi usually stores both $x^{(k)}$ and $x^{(k+1)}$.

Gauss-Seidel can overwrite $x$, so it often needs one fewer vector.

For very large sparse systems, this memory difference may matter.

90.16 Computational Cost

For a dense $n\times n$ matrix, one sweep costs

For a sparse matrix, one sweep costs approximately

The cost per iteration is similar to Jacobi.

The difference lies in convergence and parallelism.

Gauss-Seidel usually reduces error faster per sweep, but Jacobi is easier to parallelize.

90.17 Sparse Implementation

For a sparse matrix stored by rows, one Gauss-Seidel update has the form:

for i = 1 to n:
    s = 0
    diag = 0

    for each nonzero (j, aij) in row i:
        if j == i:
            diag = aij
        else:
            s = s + aij * x[j]

    x[i] = (b[i] - s) / diag

This works because values $x[j]$ for $j<i$ have already been updated during the current sweep, while values for $j>i$ still come from the previous sweep.

The ordering of unknowns therefore affects the method.

90.18 Dependence on Ordering

Gauss-Seidel depends on the order in which variables are updated.

Changing the order changes the splitting of $A$ into lower and upper parts.

For some matrices, one ordering may converge faster than another.

Ordering matters especially for sparse matrices from grids or graphs.

Common orderings include:

This dependence is both a weakness and a tool.

90.19 Parallel Gauss-Seidel

Basic Gauss-Seidel is sequential because each update depends on previous updates in the same sweep.

Parallel versions are possible when variables can be grouped so that updates within a group do not depend on each other.

For example, red-black Gauss-Seidel splits grid points into two sets. All red points are updated first, then all black points.

Graph coloring generalizes this idea.

If two variables are not coupled by a nonzero matrix entry, they may be updated simultaneously.

90.20 Symmetric Gauss-Seidel

Symmetric Gauss-Seidel performs two sweeps:

1. a forward sweep,
2. a backward sweep.

The forward sweep uses the usual order:

The backward sweep uses the reverse order:

This produces a more symmetric operation.

Symmetric Gauss-Seidel is often used as a preconditioner, especially for symmetric systems.

90.21 Gauss-Seidel as Preconditioner

From the residual form,

This suggests the preconditioner

Applying the preconditioner means solving a lower triangular system.

This is more expensive than Jacobi preconditioning but may be stronger because it includes lower-triangular coupling.

For symmetric methods, symmetric Gauss-Seidel or SSOR is often preferred.

90.22 Relation to Coordinate Descent

For symmetric positive definite $A$, solving

is equivalent to minimizing the quadratic function

Gauss-Seidel may be interpreted as exact coordinate minimization of this quadratic.

Each update chooses one coordinate $x_i$ while holding the others fixed, minimizing $\phi$ along that coordinate direction.

This interpretation explains why Gauss-Seidel converges for symmetric positive definite systems.

90.23 Failure Modes

Gauss-Seidel may fail or perform poorly for several reasons.

A robust implementation should check diagonal entries and enforce a maximum iteration count.

90.24 Residual Versus Error

As with all iterative solvers, the residual does not directly equal the error.

For

and

we have

If $A$ is invertible,

Thus

If $A$ is ill-conditioned, a small residual may still correspond to a significant solution error.

90.25 Practical Use

Gauss-Seidel is simple, but it is rarely the best standalone method for large difficult systems.

It remains useful because it is:

In modern large-scale linear algebra, Gauss-Seidel is often used inside larger algorithms rather than as the final solver.

90.26 Summary

The Gauss-Seidel method solves

by the iteration

In matrix form,

The iteration matrix is

The method converges for every initial guess exactly when

Strict diagonal dominance and symmetric positive definiteness are important sufficient conditions for convergence. Gauss-Seidel often converges faster than Jacobi, but it is more sequential and depends on variable ordering.