Deriving the Gauss-Seidel Iteration

From a Linear System to an Iteration

For a system \(\mathbf{A}\mathbf{x}=\mathbf{b}\), a direct method seeks the solution by a finite sequence of operations. An iterative method instead begins with a chosen vector \(\mathbf{x}^{(0)}\) and constructs successive approximations \(\mathbf{x}^{(1)},\mathbf{x}^{(2)},\ldots\). The Gauss-Seidel method is one such method. Its defining feature is the order in which it uses the components: during a sweep, each component is updated in turn, and the newly computed values are used immediately in the remaining updates.

Work over \(\mathbb{R}\), with \(\mathbf{A}\in\mathbb{R}^{n\times n}\), \(\mathbf{x},\mathbf{b}\in\mathbb{R}^n\), and \(a_{ii}\ne0\) for every \(i\). The nonzero diagonal condition matters: each component update divides by \(a_{ii}\). Write the system row by row as \(a_{ii}x_i=b_i-\sum_{j\ne i}a_{ij}x_j\). Gauss-Seidel turns these equations into a sequence of updates by using new values for components with indices below \(i\), and the old values for components with indices above \(i\).

The order is part of the method: the updates are performed from component 1 to component \(n\). In particular, when \(i=1\), there are no new components to use, while when \(i=n\), every other component in that equation has already been updated. The method is well-defined for every initial vector when all diagonal entries are nonzero. This condition does not by itself say that the iterates approach the solution; convergence is a separate question.

The Matrix Form of a Sweep

Let \(\mathbf{D}\) contain the diagonal entries of \(\mathbf{A}\), let \(\mathbf{L}\) contain its entries strictly below the diagonal, and let \(\mathbf{U}\) contain its entries strictly above the diagonal. Then \(\mathbf{A}=\mathbf{D}+\mathbf{L}+\mathbf{U}\). The entries of \(\mathbf{L}\) and \(\mathbf{U}\) retain their signs as they appear in \(\mathbf{A}\). The component update can now be written as one matrix equation.

Proof. In row \(i\), the product \((\mathbf{D}+\mathbf{L})\mathbf{x}^{(k+1)}\) is \(a_{ii}x_i^{(k+1)}+\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}\). The row \(i\) entry of \(\mathbf{U}\mathbf{x}^{(k)}\) is \(\sum_{j=i+1}^{n}a_{ij}x_j^{(k)}\). Consequently, row \(i\) of the matrix equation is precisely the component update after multiplying that update by \(a_{ii}\) and moving the new lower-index terms to the left. This holds for every \(i\), so the matrix equation and the component sweep are equivalent.

To check invertibility, consider \((\mathbf{D}+\mathbf{L})\mathbf{z}=\mathbf{0}\). Its first row gives \(a_{11}z_1=0\), so \(z_1=0\). If \(z_1,\ldots,z_{i-1}=0\), row \(i\) gives \(a_{ii}z_i=0\), so \(z_i=0\). Induction shows that the only solution is \(\mathbf{z}=\mathbf{0}\). Also, for any right-hand side, each row can be solved in order by dividing by the nonzero \(a_{ii}\), giving a unique solution. Thus \(\mathbf{D}+\mathbf{L}\) is invertible. Multiplying the matrix equation by its inverse gives the displayed iteration formula. \(\square\)

The matrix form summarizes the same update rule; it does not change it. In particular, the lower-triangular part \(\mathbf{D}+\mathbf{L}\) acts on the new iterate, while the strictly upper-triangular part acts on the old iterate. Since \(\mathbf{D}+\mathbf{L}\) is lower triangular with nonzero diagonal, each sweep can be carried out by forward substitution, rather than by forming its inverse.

Working Through the Component Updates

This last check illustrates how to interpret an intermediate iterate. It is an approximation, not automatically an exact solution. Also, Gauss-Seidel does not solve every row simultaneously using the same vector: it satisfies each row update in sequence, with a mixture of new and old components during the sweep.

For this example the matrix split makes the update especially visible:

The first and second entries of \(\mathbf{U}\mathbf{x}^{(k)}\) use old components, while solving with \(\mathbf{D}+\mathbf{L}\) in order incorporates the new components as soon as they are available.

Fixed Points and Solutions

The matrix form also gives a precise connection between the iteration and the original system. A fixed point of an iteration is a vector that is unchanged by one application of its update rule. This connection identifies what Gauss-Seidel is trying to reach, without asserting that a particular sequence of iterates will reach it.

Proof. The update is \((\mathbf{D}+\mathbf{L})\mathbf{x}^{(k+1)}=\mathbf{b}-\mathbf{U}\mathbf{x}^{(k)}\). If \(\mathbf{x}\) is a fixed point, set \(\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}=\mathbf{x}\). Then \((\mathbf{D}+\mathbf{L})\mathbf{x}=\mathbf{b}-\mathbf{U}\mathbf{x}\). Adding \(\mathbf{U}\mathbf{x}\) to both sides gives \((\mathbf{D}+\mathbf{L}+\mathbf{U})\mathbf{x}=\mathbf{b}\), which is \(\mathbf{A}\mathbf{x}=\mathbf{b}\). Conversely, if \(\mathbf{A}\mathbf{x}=\mathbf{b}\), use \(\mathbf{A}=\mathbf{D}+\mathbf{L}+\mathbf{U}\) and rearrange to obtain \((\mathbf{D}+\mathbf{L})\mathbf{x}=\mathbf{b}-\mathbf{U}\mathbf{x}\). Since \(\mathbf{D}+\mathbf{L}\) is invertible, this is exactly the update with the same vector on both sides as the old and new iterate. Thus \(\mathbf{x}\) is a fixed point. If \(\mathbf{A}\) is invertible, the system has a unique solution, so this fixed point is unique. \(\square\)

Invertibility of \(\mathbf{A}\) is not needed for the equivalence: every fixed point solves the system, and every solution is a fixed point. Invertibility is needed for the stated uniqueness. Neither the fixed-point equivalence nor uniqueness establishes that the sequence from an arbitrary initial vector approaches a fixed point.

Common Mistake: Reusing an Old Component

In the full three-component example, a common error is to calculate the second component from the old first component, as if every component were being updated from \(\mathbf{x}^{(0)}\). That gives the incorrect Gauss-Seidel calculation \(x_2^{(1)}=(5-0-0)/3=5/3\). It fails to use the just-computed value \(x_1^{(1)}=2\). The correct calculation is \(x_2^{(1)}=(5-2-0)/3=1\). The third update must then use both new values: \(x_3^{(1)}=(6-2-1)/4=3/4\).

Using old values for every component is the defining update pattern of the Jacobi iteration, not Gauss-Seidel. A reliable way to avoid mixing the methods is to write the update equations in index order and mark each component as soon as it has been computed. For \(x_i^{(k+1)}\), indices below \(i\) use the new iterate, and indices above \(i\) use the old one.

Finally, a zero diagonal entry prevents the displayed component rule from being used in that row because it requires division by \(a_{ii}\). This does not by itself show that the system has no solution; it means this particular ordering and direct component update are not available as written. One must not silently divide by zero or infer a conclusion about the system’s solvability from the failure of the update formula.