Jacobi method

Writing $A=D+L+U$, where $D$ is diagonal and $L,U$ are strictly triangular, the Jacobi iteration is

The Jacobi method solves a linear system by updating every coordinate from the same old vector. It requires every diagonal entry $a_{ii}$ to be nonzero. Isolating coordinate $i$ in row $i$ gives the simultaneous ruleHere $B$ is the iteration matrix. A sweep computes a complete new vector before replacing the old one, so its coordinate updates can run independently.^[ref]

For a sparse matrix, a sweep needs work proportional to the stored nonzero entries, while Gaussian elimination may create fill-in. The Jacobi method is therefore a simple starting point for large linear systems, and the weighted Jacobi method supplies a useful building block of a multigrid method. Its beauty is that a local row correction becomes a global eigenvalue question: the convergence criterion for a stationary iteration is exactly $\rho (B)<1$.

Consider the original linear system with a symmetric matrix that is also a tridiagonal matrix:Here $D=4I$ and the Jacobi method readsStarting from $x^{(0)}=0$, the first three vectors areEvery right-hand side uses superscript $k$, including the middle row, even when a new first coordinate has already been computed.

For comparison, a first Gauss-Seidel method sweep in coordinate order $1,2,3$ givesThe new first coordinate enters the middle update, and the new second coordinate enters the final update. Accidentally overwriting old coordinates during a Jacobi method sweep therefore changes the mathematical method.

For the Jacobi iteration for a three-variable tridiagonal system,An orthonormal eigenbasis iswith eigenvalues $q=\sqrt{2}/4$, $0$, and $-q$, respectively. These are error modes: their shapes persist under the iteration matrix, while their coefficients change. Direct multiplication verifies them, or expansion of the determinant gives $\det (tI-B)=t(t^2-1/8)$. Thus $\rho (B)=q<1$ and convergence holds for every starting vector.

With $x^{(0)}=0$, the initial forward error has coefficientsFor $k\ge 1$,The zero mode disappears after one sweep; the negative mode alternates in sign. The Euclidean norm drops by exactly $q$ per sweep from $k=1$ onward. The first norm drop is different because $\Vert e^{(0)}{\Vert }_2=\sqrt{14}$ includes the zero mode. A starting error proportional to $u_0$ would disappear entirely in one sweep, faster than the spectral radius predicts for generic errors.

The curves are computed iterates and exact formulas for this linear system. Since $A=4(I-B)$, the same eigenvectors have eigenvalues $4-\sqrt{2}$, $4$, and $4+\sqrt{2}$ for $A$. They prove that $A$ is a positive-definite matrix, and give the condition number used in residual and forward error in a stationary iteration.

A strictly diagonally dominant matrix guarantees convergence of the Jacobi method for every initial vector.^[ref]

Indeed, using the maximum norm $\Vert v{\Vert }_{\infty }={\max }_i|v_i|$, the induced operator norm is the largest absolute row sum. The Jacobi method has $b_{ii}=0$ and $b_{ij}=-a_{ij}/a_{ii}$ for $j\ne i$, soThis gives a cheap sufficient check without computing any eigenvalues. For the Jacobi iteration for a three-variable tridiagonal system, the absolute row sums of $B$ are $1/4,1/2,1/4$, giving $\Vert B{\Vert }_{\infty }=1/2$. The true spectral radius $\sqrt{2}/4$ gives a sharper eventual rate.

Strict diagonal dominance is sufficient, not necessary: $A=\left(\begin{array}{cc}1& 2\\ 1/10& 1\end{array}\right)$ violates it in the first row, but its Jacobi method has eigenvalues $\pm 1/\sqrt{5}$ and converges. Conversely, being a positive-definite matrix alone is insufficient: the Jacobi convergence for a three-by-three equicorrelation matrix with $\alpha =3/4$ has positive eigenvalues $5/2,1/4,1/4$ for $A$, but spectral radius $3/2$ for its iteration matrix. The extra tridiagonal matrix structure in Jacobi convergence for a symmetric positive-definite tridiagonal matrix is therefore substantive.

For $A={\mu }^{-1}I+U-U^T$, with $U$ the four-dimensional unit superdiagonal, the Jacobi method has eigenvalues $2i\mu \cos (j\pi /5)$, $j=1,\dots ,4$. The Gauss-Seidel method has two zero eigenvalues and $-{\mu }^2(3\pm \sqrt{5})/2$. Therefore both methods converge for every starting vector exactly when $0<|\mu |<(\sqrt{5}-1)/2$. This equality of convergence ranges follows from their spectra, not from positive-definiteness of the nonsymmetric coefficient matrix.

With relaxation parameter $\omega$, weighted Jacobi updates byIts iteration matrix is $H_{\omega }=I-\omega D^{-1}A$.

Equivalently, the weighted Jacobi method takes a weighted average of the old vector and its full Jacobi method update: $x^{(k+1)}=(1-\omega )x^{(k)}+\omega x_{\mathrm{Jacobi}}^{(k+1)}$. If $Bv=\lambda v$, its iteration matrix multiplies that eigenvector by $1-\omega +\omega \lambda$. For the diagonal-4 Jacobi iteration for a three-variable tridiagonal system, the factors are $1-\omega$ and $1-\omega \pm \omega \sqrt{2}/4$; a Poisson equation smoothing parameter must not be inferred from this different matrix.

For the five-point discretization of $-\Delta +c$ on the unit square, the weighted iteration eigenvalues are $1-\omega [4{\sin }^2(\pi kh/2)+4{\sin }^2(\pi lh/2)+c{h}^2]/(4+c{h}^2)$. A fixed parameter converges for every mesh when $0<\omega \le 1$, but $\omega =1$ poorly damps the highest frequencies.

For the five-point Dirichlet Poisson matrix, weighted Jacobi iteration has sine-mode multiplier $1-\omega +(\omega /2)(\cos {\theta }_r+\cos {\theta }_s)$. Call a mode high-frequency if at least one cosine is nonpositive. Then the sum lies in $[-2,1]$, so for $\omega =4/5$ the multiplier lies in $[-3/5,3/5]$. This uniform damping makes the iteration a useful smoother for a multigrid method, even though low-frequency multipliers approach one on fine grids. Balancing $|1-2\omega |$ and $|1-\omega /2|$ shows why $4/5$ is the minimax parameter for this high-frequency interval. The two-dimensional high-frequency set gives a different optimum from the familiar one-dimensional value $2/3$.

For the centered second-difference matrix with diagonal $-2$ and off-diagonals $1$, the weighted-Jacobi eigenvalues areThe method converges for every grid size exactly when $0<\omega \le 1$.

Multiplying both $A$ and $b$ by $-1$ leaves $D^{-1}A$ and the weighted Jacobi method unchanged. Thus the same formula applies to the positive second-difference matrix $T_n$ with diagonal $2$ and neighboring entries $-1$, the finite difference form of the one-dimensional Poisson equation with Dirichlet boundary conditions. Its eigenvectors have components $\sin (jk\pi /(n+1))$, and $D^{-1}{T}_n$ has eigenvalues $1-\cos (k\pi /(n+1))$. At $\omega =2/3$, upper-half modes have factors between $-1/3$ and $1/3$, while low-frequency factors remain close to one on fine grids. This explains why the weighted Jacobi method helps a multigrid method by reducing oscillatory errors, even when it is a slow solver on its own.^[ref]

For modes with $k/(n+1)\in [1/2,1]$ and large $n$, the worst attenuation factor isIt is minimized by $\omega =2/3$, for which ${\mu }_{\omega }=1/3$.

For the one-dimensional Dirichlet second-difference grid, high sine modes have eigenvalues $1-\omega t_k$, with $t_k=1-\cos (k\pi /(n+1))$ and $k\ge \lceil (n+1)/2\rceil$. Let $a$ and $b$ be the smallest and largest of these $t_k$. Convexity gives attenuation $\max (|1-\omega a|,|1-\omega b|)$. Balancing its endpoint values with opposite signs gives the displayed exact minimax choice. For odd $n$, $a=1$; for even $n$, $a=1+\sin (\pi /[2(n+1)])$, and always $b=1+\cos (\pi /(n+1))$. Only in the large-grid or uniform-in-grid limit do these become ${\omega }_{\ast }=2/3$, ${\mu }_{\ast }=1/3$.

For a three-by-three matrix with diagonal entries one and every off-diagonal entry $\alpha$, the Jacobi method iteration matrix has eigenvalue $-2\alpha$ on the all-ones vector and eigenvalue $\alpha$ with multiplicity two on its orthogonal complement. Its spectral radius is $2|\alpha |$, so Jacobi iteration converges exactly when $|\alpha |<1/2$.

If $A=D+L+L^T$ is symmetric positive definite and tridiagonal, set $S=\mathrm{diag}(1,-1,1,-1,\dots )$. Thenis positive definite. Applying the Householder-John theorem to $M=D$ and $N=-(L+L^T)$ proves that the Jacobi iteration converges.