Chapter 89. Jacobi Method

The Jacobi method is a stationary iterative method for solving a linear system

Ax=b.

It is one of the simplest iterative solvers in numerical linear algebra. The method updates each unknown by using the values from the previous iteration. This makes it easy to understand, easy to implement, and naturally parallel.

The method is most useful as a basic model for iterative solution, as a smoother in multigrid methods, and as a simple preconditioner. By itself, it may converge slowly, but it introduces the main ideas behind stationary methods: matrix splitting, iteration matrices, residuals, convergence criteria, and stopping rules.

89.1 The Basic Problem

Let

A= \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \qquad x= \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix}, \qquad b= \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{bmatrix}.

The system

Ax=b

is equivalent to the equations

a_{i1}x_1+a_{i2}x_2+\cdots+a_{in}x_n=b_i, \qquad i=1,\ldots,n.

Assume that

a_{ii}\ne 0

for every $i$ . Then the $i$ -th equation can be solved for $x_i$ :

x_i = \frac{1}{a_{ii}} \left( b_i-\sum_{j\ne i}a_{ij}x_j \right).

The Jacobi method turns this identity into an iteration.

89.2 Jacobi Update Formula

Starting from an initial guess

x^{(0)},

the Jacobi method defines

x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i-\sum_{j\ne i}a_{ij}x_j^{(k)} \right), \qquad i=1,\ldots,n.

All entries of $x^{(k+1)}$ are computed from entries of $x^{(k)}$ , not from entries already updated during the same iteration.

This is the defining feature of the Jacobi method.

The update is simultaneous:

x^{(k)} \longmapsto x^{(k+1)}.

Each component uses old values from the previous iterate.

89.3 Example

Consider the system

\begin{aligned} 4x_1+x_2 &= 9,\\ x_1+3x_2 &= 7. \end{aligned}

Solving each equation for its diagonal variable gives

x_1=\frac{9-x_2}{4}, \qquad x_2=\frac{7-x_1}{3}.

Thus the Jacobi iteration is

x_1^{(k+1)}=\frac{9-x_2^{(k)}}{4}, \qquad x_2^{(k+1)}=\frac{7-x_1^{(k)}}{3}.

Start with

x^{(0)}= \begin{bmatrix} 0\\ 0 \end{bmatrix}.

Then

x_1^{(1)}=\frac{9}{4}=2.25, \qquad x_2^{(1)}=\frac{7}{3}\approx 2.3333.

x^{(1)} = \begin{bmatrix} 2.25\\ 2.3333 \end{bmatrix}.

Next,

x_1^{(2)} = \frac{9-2.3333}{4} \approx 1.6667,

x_2^{(2)} = \frac{7-2.25}{3} \approx 1.5833.

The exact solution is

x= \begin{bmatrix} 2\\ 5/3 \end{bmatrix}.

The iterates move toward this solution.

89.4 Matrix Splitting

Write the matrix $A$ as

A = D + L + U,

where:

Symbol	Meaning
$D$	diagonal part of $A$
$L$	strictly lower triangular part of $A$
$U$	strictly upper triangular part of $A$

Then

Ax=b

becomes

(D+L+U)x=b.

Move the off-diagonal terms to the right side:

Dx=-(L+U)x+b.

The Jacobi iteration is

Dx^{(k+1)} = -(L+U)x^{(k)}+b.

Since $D$ is diagonal, solving with $D$ is cheap:

x^{(k+1)} = -D^{-1}(L+U)x^{(k)}+D^{-1}b.

Thus the Jacobi iteration matrix is

G_J=-D^{-1}(L+U).

89.5 Fixed-Point Form

The Jacobi method has the fixed-point form

x^{(k+1)}=G_Jx^{(k)}+c,

where

G_J=-D^{-1}(L+U), \qquad c=D^{-1}b.

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

x=G_Jx+c.

Substituting the definitions gives

x=-D^{-1}(L+U)x+D^{-1}b.

Multiplying by $D$ ,

Dx=-(L+U)x+b.

Therefore,

(D+L+U)x=b.

So the limit is a solution of the original system.

89.6 Error Recurrence

Let $x$ be the exact solution. The exact solution satisfies

x=G_Jx+c.

The computed iterate satisfies

x^{(k+1)}=G_Jx^{(k)}+c.

Subtracting gives

x-x^{(k+1)} = G_J(x-x^{(k)}).

Thus the error

e^{(k)}=x-x^{(k)}

satisfies

e^{(k+1)}=G_Je^{(k)}.

By repeated substitution,

e^{(k)}=G_J^k e^{(0)}.

The convergence of Jacobi is therefore controlled entirely by powers of the iteration matrix.

89.7 Convergence Criterion

The Jacobi method converges for every initial guess if and only if

\rho(G_J)<1,

where $\rho(G_J)$ is the spectral radius of $G_J$ . The spectral radius is the largest absolute value of the eigenvalues of $G_J$ .

This is the standard convergence condition for stationary linear iterations. The same criterion is commonly used for Jacobi and related stationary methods.

\rho(G_J)<1,

then

G_J^k \to 0

k\to\infty.

Consequently,

e^{(k)}=G_J^ke^{(0)}\to 0.

The iterates converge to the exact solution.

89.8 Diagonal Dominance

A common sufficient condition for convergence is strict diagonal dominance.

A matrix $A$ is strictly row diagonally dominant if

$$ |a_{ii}|

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

for every row $i$ .

This means each diagonal entry is larger in magnitude than the sum of the off-diagonal entries in the same row.

For strictly diagonally dominant systems, the Jacobi method is commonly introduced as a convergent iterative algorithm.

Strict diagonal dominance is sufficient, but not necessary. Some matrices that are not strictly diagonally dominant still allow Jacobi to converge.

89.9 Why Diagonal Dominance Helps

The Jacobi update divides each equation by its diagonal coefficient.

If the diagonal coefficient dominates the off-diagonal coefficients, then the new value of $x_i$ depends only weakly on the other variables.

In that case, the iteration matrix has small enough entries to contract errors.

For example, if

$$ |a_{ii}|

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

then the sum of the absolute values in row $i$ of $G_J$ is less than $1$ :

\sum_{j=1}^n |(G_J)_{ij}| = \sum_{j\ne i} \left| \frac{a_{ij}}{a_{ii}} \right| <1.

Thus

\|G_J\|_\infty < 1.

Since

\rho(G_J)\le \|G_J\|_\infty,

we get

\rho(G_J)<1.

Therefore, Jacobi converges.

89.10 Residual Form

The residual at step $k$ is

r^{(k)}=b-Ax^{(k)}.

The Jacobi correction can be written using the residual.

Since

x^{(k+1)} = x^{(k)}+D^{-1}r^{(k)},

Jacobi updates the current approximation by adding a diagonally scaled residual.

This form is useful in implementation and interpretation.

It says that Jacobi corrects each variable according to the residual in its corresponding equation, scaled by the diagonal entry.

89.11 Algorithm

Input:

A,\quad b,\quad x^{(0)},\quad \tau,\quad k_{\max}.

Repeat for

k=0,1,\ldots,k_{\max}-1:

Compute, for each $i$ ,

x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i-\sum_{j\ne i}a_{ij}x_j^{(k)} \right).

Compute the residual

r^{(k+1)}=b-Ax^{(k+1)}.

Stop if

\frac{\|r^{(k+1)}\|}{\|b\|} \le \tau.

Return

x^{(k+1)}.

89.12 Pseudocode

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

    for k = 0 to max_iter - 1:
        x_new = zero vector of length n

        for i = 1 to n:
            s = 0

            for j = 1 to n:
                if j != i:
                    s = s + A[i,j] * x[j]

            x_new[i] = (b[i] - s) / A[i,i]

        r = b - A * x_new

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

        x = x_new

    return x

The new vector must be stored separately from the old vector. Otherwise the method becomes Gauss-Seidel rather than Jacobi.

89.13 Storage Requirements

Jacobi needs:

Object	Storage
Matrix $A$	problem-dependent
Right-hand side $b$	$n$
Current vector $x^{(k)}$	$n$
New vector $x^{(k+1)}$	$n$
Residual $r^{(k)}$ , optional	$n$

For sparse matrices, the matrix is stored using a sparse format such as compressed sparse row.

The method does not require factorizing $A$ . This makes it memory-light compared with direct methods.

89.14 Computational Cost

For a dense $n\times n$ matrix, each Jacobi iteration costs

O(n^2)

operations.

For a sparse matrix, each iteration costs approximately

O(\operatorname{nnz}(A)),

where $\operatorname{nnz}(A)$ is the number of nonzero entries.

The cost per iteration is low, but many iterations may be required.

This is the main tradeoff.

89.15 Parallelism

Jacobi is naturally parallel.

Each component

x_i^{(k+1)}

depends only on values from the previous iterate $x^{(k)}$ . Therefore all components of $x^{(k+1)}$ can be computed independently after $x^{(k)}$ is known.

This makes Jacobi attractive on parallel hardware.

However, parallel efficiency alone does not guarantee overall efficiency. If convergence is slow, a faster converging method with less parallelism may still be better.

89.16 Weighted Jacobi

Weighted Jacobi modifies the update by mixing the old iterate with the ordinary Jacobi iterate.

Let

\tilde{x}^{(k+1)} = D^{-1}\bigl(b-(L+U)x^{(k)}\bigr)

be the ordinary Jacobi update.

Weighted Jacobi defines

x^{(k+1)} = (1-\omega)x^{(k)} + \omega \tilde{x}^{(k+1)}.

Here $\omega$ is a relaxation parameter.

Equivalently,

x^{(k+1)} = x^{(k)}+\omega D^{-1}r^{(k)}.

When

\omega=1,

weighted Jacobi becomes ordinary Jacobi.

For some problems, choosing

0<\omega<1

improves smoothing behavior.

89.17 Jacobi as a Smoother

In multigrid methods, Jacobi is often used as a smoother.

A smoother reduces high-frequency error components quickly, even if it reduces low-frequency components slowly.

This property is useful because multigrid handles different error scales on different grids.

Weighted Jacobi is especially common in this context.

The method may be poor as a standalone solver but useful as a component inside a stronger algorithm.

89.18 Jacobi Preconditioning

Jacobi preconditioning uses the diagonal of $A$ as a preconditioner.

Set

P=D.

Then the preconditioned system is

D^{-1}Ax=D^{-1}b.

This scales each equation by its diagonal entry.

Jacobi preconditioning is cheap and easy to apply. It may help when diagonal scaling is the main issue.

It is weak for strongly coupled systems, where off-diagonal entries play a major role.

89.19 Comparison with Gauss-Seidel

Jacobi and Gauss-Seidel are closely related.

Feature	Jacobi	Gauss-Seidel
Update style	simultaneous	sequential
Uses newest values immediately	no	yes
Needs separate new vector	yes	no
Parallelism	high	lower
Typical convergence speed	slower	faster
Implementation simplicity	high	high

Gauss-Seidel often converges faster because it uses updated values during the same iteration. Jacobi is easier to parallelize because each update is independent.

89.20 Failure Modes

Jacobi may fail for several reasons.

Failure mode	Cause
Divergence	$\rho(G_J)\ge 1$
Slow convergence	$\rho(G_J)$ close to $1$
Division by zero	some $a_{ii}=0$
Poor scaling	diagonal entries too small
Misleading residual	ill-conditioned $A$

Before applying Jacobi, one should check that diagonal entries are nonzero and that convergence is plausible.

89.21 Residual Versus Error

A small residual does not always mean a small error.

For the residual

r^{(k)}=b-Ax^{(k)},

the error satisfies

e^{(k)}=A^{-1}r^{(k)}.

Thus

\|e^{(k)}\| \le \|A^{-1}\|\|r^{(k)}\|.

If $A$ is ill-conditioned, then $\|A^{-1}\|$ may be large. The residual can be small while the solution error remains significant.

This is not specific to Jacobi. It applies to all linear solvers.

89.22 Practical Use

Jacobi is rarely the best standalone method for difficult linear systems.

It remains important because it is:

Role	Reason
Teaching method	exposes core iteration ideas
Parallel method	updates are independent
Preconditioner	diagonal scaling is cheap
Multigrid smoother	damps high-frequency error
Baseline solver	simple reference implementation

In production solvers, Jacobi is often replaced or supplemented by Krylov methods, multigrid methods, domain decomposition, or stronger preconditioners.

89.23 Implementation Notes

A robust implementation should:

Issue	Practical handling
Zero diagonal	reject matrix or reorder equations
Small diagonal	warn or scale system
Sparse storage	avoid looping over zeros
Residual check	use relative residual
Maximum iterations	always enforce a limit
Parallel update	separate old and new vectors
Weighted update	expose $\omega$ as parameter

For sparse matrices, it is inefficient to compute sums over all $j=1,\ldots,n$ . Instead, loop only over nonzero entries in each row.

89.24 Simple Sparse Row Update

For a sparse row, the update is computed by separating the diagonal entry from the off-diagonal entries:

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_new[i] = (b[i] - s) / diag

This structure avoids unnecessary operations on zero entries.

It is the natural implementation for compressed sparse row storage.

89.25 Summary

The Jacobi method is the iteration

x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i-\sum_{j\ne i}a_{ij}x_j^{(k)} \right).

In matrix form,

x^{(k+1)} = -D^{-1}(L+U)x^{(k)}+D^{-1}b.

Its iteration matrix is

G_J=-D^{-1}(L+U).

The method converges for every initial guess exactly when

\rho(G_J)<1.

Jacobi is simple, parallel, and memory-light. Its main weakness is slow or absent convergence. It is often more useful as a building block than as a final solver.