Gauss-Seidel iteration

You have just finished Gauss and Gauss–Jordan elimination — direct methods that march through a matrix once and hand you the exact answer. The Gauss–Seidel method is a different animal entirely: instead of eliminating unknowns, it starts from a guess, however crude, and keeps correcting it, equation by equation, until the numbers stop moving. Each correction uses whatever values are freshest — even ones updated moments earlier in the very same pass. For a large, sparse, diagonally dominant system, this cheap repeated correction can beat elimination outright, since it never has to form or store the full matrix reduction. This article gives you the update formula, the diagonal-dominance test that guarantees convergence, and enough hand-worked iterations to make both routine.

The core insight

Gauss–Seidel solves $A\mathbf x=\mathbf b$ by repeatedly recomputing each unknown from the others using the most recently updated values, converging to the true solution whenever $A$ is diagonally dominant.

The core idea

Take the system $A\mathbf x=\mathbf b$ with $n$ equations and solve each equation for its own diagonal unknown. This is the same rearrangement you would do by hand for a $3\times3$ system:

x_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j\neq i}a_{ij}x_j\right),\qquad i=1,\dots,n.

where: $a_{ii}$ is the diagonal coefficient of row $i$ , $a_{ij}$ ( $j\neq i$ ) are the off-diagonal coefficients, and $b_i$ is the constant on the right-hand side of equation $i$ . Starting from an initial guess $\mathbf x^{(0)}$ — often just $\mathbf x^{(0)}=\mathbf 0$ — Gauss–Seidel sweeps through the unknowns in order $1,2,\dots,n$ and updates each one using the newest available values on both sides:

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).

where: the first sum uses values already updated this iteration ( $k+1$ ), and the second sum uses values not yet touched this pass, still at the old iteration $k$ . That single detail — reusing fresh values immediately — is what separates Gauss–Seidel from the Jacobi method, which would use only $k$ -th-iteration values on the entire right-hand side. Because Gauss–Seidel folds in the newest information as soon as it exists, it typically converges in fewer iterations than Jacobi for the same system. The figure below shows the difference in information flow within one sweep.

Figure

Diagram comparing Jacobi and Gauss Seidel update flow within one iteration, showing Gauss Seidel feeding each new value forward into the next unknown while Jacobi holds all values back until the sweep ends

Convergence. The iteration is not guaranteed for an arbitrary matrix $A$ ; you need a sufficient condition first. The syllabus-level test is diagonal dominance:

|a_{ii}|>\sum_{j\neq i}|a_{ij}|\quad\text{for every row } i,

i.e., each diagonal entry's magnitude exceeds the sum of the other entries' magnitudes in its own row. If $A$ is strictly diagonally dominant, Gauss–Seidel converges to the exact solution from any starting guess. When the given ordering fails this test, try permuting the rows — reordering equations so the largest coefficient in each equation sits on the diagonal sometimes restores dominance without changing the solution set. Stop iterating once successive approximations are close enough: $\max_i|x_i^{(k+1)}-x_i^{(k)}|<\varepsilon$ for a chosen tolerance $\varepsilon$ . Compared with the direct methods of the previous article, Gauss–Seidel earns its keep mainly on large, sparse systems — the kind that come from discretizing a physical problem, where most entries of $A$ are already zero. Each sweep costs roughly $O(n^2)$ arithmetic operations for a dense matrix (fewer still if most entries are zero and are skipped), against roughly $O(n^3)$ for a full Gauss or Gauss–Jordan elimination on the same size system. Elimination also risks fill-in: zero entries can turn nonzero partway through the row reduction, which both costs memory and destroys the sparsity that made the problem cheap to store in the first place. Gauss–Seidel never eliminates anything, so the zero pattern of $A$ is never disturbed — only the guess $\mathbf x^{(k)}$ changes from sweep to sweep. For a small, dense, well-behaved system, however, a single elimination pass is often less work overall than iterating to convergence, so the choice between the two is a genuine engineering trade-off, not a blanket rule that one method is always better.

Derivation AWhy diagonal dominance forces the error to shrink

The goal is to show, from the update formula itself and nothing more advanced than the triangle inequality, that diagonal dominance guarantees the error in each unknown gets smaller with every iteration, so that after enough sweeps $\mathbf x^{(k)}$ must converge to the exact solution $\mathbf x$ regardless of how bad the starting guess was.

Step 1. Define the error at unknown $i$ after iteration $k$ as $e_i^{(k)}=x_i^{(k)}-x_i$ , the gap between the current iterate and the true value. Subtracting the exact relation $x_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j\neq i}a_{ij}x_j\right)$ from the Gauss–Seidel update $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)$ cancels $b_i$ entirely and leaves the errors talking directly to each other: $e_i^{(k+1)}=-\frac{1}{a_{ii}}\left(\sum_{j<i}a_{ij}e_j^{(k+1)}+\sum_{j>i}a_{ij}e_j^{(k)}\right)$ .

Step 2. Apply the triangle inequality to bound the size of the new error. Since $|e_j^{(k+1)}|$ and $|e_j^{(k)}|$ are both at most the largest error present at that stage, $|e_i^{(k+1)}|\le\frac{1}{|a_{ii}|}\sum_{j\neq i}|a_{ij}|\cdot\max_j\big(\text{current error size}\big)$ . Write $M=\max_j|e_j^{(k)}|$ for the worst error before the sweep and $M'=\max_j|e_j^{(k+1)}|$ for the worst error after it; every term on the right is bounded by $M$ (a slight overestimate, but it makes the bound safe), so $|e_i^{(k+1)}|\le \left(\dfrac{\sum_{j\neq i}|a_{ij}|}{|a_{ii}|}\right)M$ for each row $i$ .

Step 3. Define $r_i=\dfrac{\sum_{j\neq i}|a_{ij}|}{|a_{ii}|}$ , the ratio of off-diagonal weight to diagonal weight in row $i$ , and let $r=\max_i r_i$ . Diagonal dominance means $|a_{ii}|>\sum_{j\neq i}|a_{ij}|$ for every row, which is exactly the statement $r_i<1$ for every $i$ , hence $r<1$ . Taking the worst row on both sides of the bound in Step 2 gives $M'\le rM$ — the worst-case error after one full sweep is at most $r$ times the worst-case error before it.

Step 4. Chain the bound across iterations: after $k$ sweeps, $M^{(k)}\le r^k M^{(0)}$ . Since $r<1$ is guaranteed by diagonal dominance, $r^k\to0$ as $k\to\infty$ , so every component of the error is squeezed to zero no matter what $\mathbf x^{(0)}$ you started from. This is why diagonal dominance is sufficient for convergence: it is exactly the condition that forces the contraction factor $r$ below $1$ . It is not necessary — some non-dominant systems still converge — but it is the safe test to apply by hand.

Worked example

Solve, using three Gauss–Seidel iterations from $\mathbf x^{(0)}=(0,0,0)$ :

10x+2y+z=9,\qquad 2x+20y-2z=-44,\qquad -2x+3y+10z=22.

Step 1. Check diagonal dominance before doing any arithmetic. Row 1: $|10|>|2|+|1|=3$ ✓. Row 2: $|20|>|2|+|-2|=4$ ✓. Row 3: $|10|>|-2|+|3|=5$ ✓. All three rows pass, so convergence is guaranteed and the update formulas are

x=\frac{9-2y-z}{10},\qquad y=\frac{-44-2x+2z}{20},\qquad z=\frac{22+2x-3y}{10}.

Step 2. Iteration 1, starting from $y^{(0)}=z^{(0)}=0$ : $x^{(1)}=\dfrac{9-0-0}{10}=0.9$ ; $y^{(1)}=\dfrac{-44-2(0.9)+2(0)}{20}=\dfrac{-45.8}{20}=-2.29$ ; $z^{(1)}=\dfrac{22+2(0.9)-3(-2.29)}{10}=\dfrac{30.67}{10}=3.067$ .

Step 3. Iteration 2, using $x^{(1)}=0.9,\,y^{(1)}=-2.29,\,z^{(1)}=3.067$ : $x^{(2)}=\dfrac{9-2(-2.29)-3.067}{10}=\dfrac{10.513}{10}=1.0513$ ; $y^{(2)}=\dfrac{-44-2(1.0513)+2(3.067)}{20}=\dfrac{-39.969}{20}=-1.998$ ; $z^{(2)}=\dfrac{22+2(1.0513)-3(-1.998)}{10}=\dfrac{30.097}{10}=3.010$ .

Step 4. Iteration 3, using $x^{(2)}=1.0513,\,y^{(2)}=-1.998,\,z^{(2)}=3.010$ : $x^{(3)}=\dfrac{9-2(-1.998)-3.010}{10}=\dfrac{9+3.996-3.010}{10}=\dfrac{9.986}{10}=0.9986$ ; $y^{(3)}=\dfrac{-44-2(0.9986)+2(3.010)}{20}=\dfrac{-44-1.9972+6.020}{20}=\dfrac{-39.977}{20}=-1.999$ ; $z^{(3)}=\dfrac{22+2(0.9986)-3(-1.999)}{10}=\dfrac{22+1.9972+5.997}{10}=\dfrac{29.994}{10}=2.999$ . After three iterations, $\mathbf x^{(3)}\approx(0.999,-1.999,2.999)$ , essentially the exact solution $(1,-2,3)$ (check: $10(1)+2(-2)+1(3)=10-4+3=9$ ✓; $2(1)+20(-2)-2(3)=2-40-6=-44$ ✓; $-2(1)+3(-2)+10(3)=-2-6+30=22$ ✓). Note how the third iteration barely moved the estimate compared with the first — the contraction factor $r<1$ from the derivation above is already shrinking the error fast. Final answer (three iterations): $(x,y,z)\approx(1.00,-2.00,3.00)$ .

Competitive stretch

If the same system had been handed to you with rows reordered so that row 1 read $-2x+3y+10z=22$ first, dominance would fail on that row ( $|-2|<|3|+|10|$ ), and you would need to permute the equations back to the ordering above — largest coefficient of each variable on its own diagonal — before iterating. Always check dominance on the equations as given; never assume the labelling in the problem is already convergence-friendly.

Common traps

Using the old value of $x_j$ for an index $j<i$ that was already updated this pass — that is Jacobi, not Gauss–Seidel; the entire point is reusing fresh values immediately.
Declaring a system "cannot be solved by Gauss–Seidel" just because it fails diagonal dominance as first written — try reordering the rows before giving up.
Treating diagonal dominance as necessary rather than sufficient; some non-dominant systems still converge, but dominance is the guarantee you can check by hand and cite in an exam.
Stopping at a fixed iteration count without checking $|x_i^{(k+1)}-x_i^{(k)}|$ is actually small — a slowly converging system (small margin of dominance, $r$ close to $1$ ) needs more sweeps than a hasty two-step calculation.

Ask the tutor: Want another system checked for diagonal dominance and worked through iteration by iteration?

Why bother with Gauss Seidel when Gauss elimination always gives the exact answer directly?

Direct elimination costs roughly $O(n^3)$ arithmetic operations overall and, for a large sparse matrix, can introduce nonzero "fill-in" entries where zeros used to be, which eats memory. Gauss–Seidel costs about $O(n^2)$ per sweep and never touches the zero entries, so for large, sparse, diagonally dominant systems — the kind that arise from discretizing physical problems — a handful of cheap iterations can reach an acceptable approximation faster than one expensive elimination.

What happens if I run Gauss Seidel on a system that is not diagonally dominant?

The method might still converge — diagonal dominance is a sufficient condition, not a necessary one — but there is no guarantee, and for some non-dominant systems the iterates will visibly grow larger and diverge instead of settling down. The safe, syllabus-level approach is to check dominance first; if it fails, try reordering the equations to restore it before trusting the iteration.

How do I know when to stop iterating?

Compare each unknown's new value with its previous value and stop once the largest change, $\max_i|x_i^{(k+1)}-x_i^{(k)}|$ , drops below a tolerance you have chosen, such as $0.01$ or $0.001$ . In the worked example above, the change from iteration 2 to iteration 3 was already under $0.01$ in every component, which is why three iterations were enough to quote the answer to three significant figures.

Gauss-Seidel iteration

The core idea

Worked example

Competitive stretch

Common traps

Before this

Where this leads