The Jacobi and Gauss-Seidel Methods

In Jacobi's method we take N=D, the diagonal part of A and P=L+U, where L and U are the strictly lower and upper parts of -Arespectively. The iteration is

$$\mbox{\boldmath$ x $ } _{m+1}=D^{-1}(L+U)\mbox{\boldmath$ x $ } _m+D^{-1}\mbox{\boldmath$ f $ },$$

provided D is nonsingular. The iteration matrix in this case is M_J=D^-1(L+U), and this is convergent if $\left\Vert \kern.05em M_J \kern.05em \right\Vert<1$ for some norm. The obvious ones to try are $\left\Vert \kern.05em . \kern.05em \right\Vert _\infty$ and $\left\Vert \kern.05em . \kern.05em \right\Vert _1$.

$$\left\Vert \kern.05em M_J \kern.05em \right\Vert _{\infty }=\... ...1}{j\not= i}}^n\left\vert\frac{a_{ij}}{a_{ii}}\right\vert < 1,$$

or, for $i=1,2,\ldots ,m$

$$\vert a_{ii}\vert>\sum _{\stackrel{j=1}{j\not= i}}^n\vert a_{ij}\vert,$$

which means that A is strictly diagonally dominant. The column sum norm gives a similar result, viz.,

$$\vert a_{jj}\vert>\sum _{\stackrel{i=1}{i\not= j}}^n\vert a_{ij}\vert,\quad j=1,2,\ldots ,m.$$

Note that these conditions are sufficient but not necessary; the method may converge even if both of these norms are greater than one. In such a case, we should eventually have recourse to the spectral radius of M_J.

For small systems, we can evaluate $\rho (M_J)$ from the characteristic equation, which takes the distinctive form $\det (\lambda I-M_J)=\det \{ \lambda I-D^{-1}(L+U)\} ~ \Rightarrow~ \det [\lambda D-L-U]=0$, or displayed

$$\left\vert \begin{array}{cccc}a_{11}\lambda & a_{12} & \cdot ... ...& \cdot & a_{m,m-1} & a_{nn}\lambda \end{array} \right\vert=0.$$

The solution is computed using the component form

$$x_i^{(m+1)}=-\frac{1}{a_{ii}}\sum _{\stackrel{j=1}{j\not= i}}^na_{ij}x_j^{(m)} +\frac{f_i}{a_{ii}},$$

for $i=1,2,\ldots ,m$ in any order.

The Gauss-Seidel method is derived from the component form of Jacobi's method

$$x_i^{(m+1)}=-\frac{1}{a_{ii}}\sum _{j=1}^{i-1}a_{ij}x_j^{(m)}... ...ac{1}{a_{ii}}\sum _{j=i+1}^na_{ij}x_j^{(m)}+\frac{f_i}{a_{ii}}$$

simply by replacing the superscript of x_j in the first member of the right hand side by m+1. This corresponds to using the most up-to-date values available, provided we proceed in the natural order, $i=1,2,\ldots ,n$. Thus the component form of Gauss-Seidel is

$$x_i^{(m+1)}=-\frac{1}{a_{ii}}\sum _{j=1}^{i-1}a_{ij}x_j^{(m+1... ...c{1}{a_{ii}}\sum _{j=i+1}^na_{ij}x_j^{(m)}+\frac{f_i}{a_{ii}}.$$

As in the Jacobi case, we analyse the convergence with the help of the matrix representation. From the component form this becomes

$$\mbox{\boldmath$ x $ } _{m+1}=D^{-1}L\mbox{\boldmath$ x $ } _... ...D^{-1}U\mbox{\boldmath$ x $ } _m+D^{-1}\mbox{\boldmath$ f $ },$$ (4.5)

giving $(D-L)\mbox{\boldmath$\space x $ } _{m+1}=U\mbox{\boldmath$\space x $ } _m+\mbox{\boldmath$\space f $ }$, which we can write as

$$\mbox{\boldmath$ x $ } _{m+1}=(D-L)^{-1}U\mbox{\boldmath$ x $ } _m+(D-L)^{-1}\mbox{\boldmath$ f $ } ,$$

provided that D-L is nonsingular, in which case the iteration matrix is M_GS=(D-L)^-1U. Taking row or column sum norms is not so easy as in the Jacobi case, but we obtain a comparable result in the following theorem.

Theorem 4.4.1   The Gauss-Seidel method converges to the solution of $A\mbox{\boldmath$\space x $ } =\mbox{\boldmath$\space f $ }$ if

$$r=\stackrel{\max }{_i}\sum _{\stackrel{j=1}{_{j\not= i}}}^n\left\vert\frac{a_{ij}}{a_{ii}}\right\vert<1,$$

(that is, A is strictly diagonally dominant).

Proof.  Assuming that if the method converges then $\mbox{\boldmath$\space x $ } _m\to \mbox{\boldmath$\space x $ }$, we have for the component form and its limit

$$\begin{eqnarray*}x_i^{(m+1)} & = & -\sum _{j=1}^{i-1}\frac{a_{ij}}{a_{ii}}x_j^{(... ...}x_j-\sum _{j=i+1}^n\frac{a_{ij}}{a_{ii}}x_j+\frac{f_i}{a_{ii}}, \end{eqnarray*}$$

which yields after subtraction and setting $\mbox{\boldmath$\space e $ } _m=\mbox{\boldmath$\space x $ } _m-\mbox{\boldmath$\space x $ }$

$$e_i^{(m+1)}=-\sum _{j=1}^{i-1}\frac{a_{ij}}{a_{ii}}e_j^{(m+1)}-\sum _{j=i+1}^n\frac{a_{ij}}{a_{ii}}e_j^{(m)}.$$ (4.6)

We now prove by induction on i, the component index, that

$$\left\Vert \kern.05em \mbox{\boldmath$ e $ } _{m+1} \kern.05e... ... $ } _m \kern.05em \right\Vert _{\infty },\quad m=0,1,\ldots .$$

From (7.6), $e_1^{(m+1)}=-\sum _{j=2}^n\frac{a_{1j}}{a_{11}}e_j^{(m)}$, so

$$\begin{eqnarray*}\vert e_1^{(m+1)}\vert & \leq & \sum _{j=2}^n\left\vert \frac{a... ...5em \mbox{\boldmath$ e $ } _m \kern.05em \right\Vert _{\infty }. \end{eqnarray*}$$

Assume that $\vert e_k^{(m+1)}\vert \leq r\left\Vert \kern.05em \mbox{\boldmath$\space e $ } _m \kern.05em \right\Vert _{\infty }$ for $k=1,2,\ldots ,i-1$. From (7.6)

$$\begin{eqnarray*}\vert e_i^{(m+1)}\vert & \leq & \sum _{j=1}^{i-1}\left\vert \fr... ...5em \mbox{\boldmath$ e $ } _m \kern.05em \right\Vert _{\infty }. \end{eqnarray*}$$

Hence $\left\Vert \kern.05em \mbox{\boldmath$\space e $ } _{m+1} \kern.05em \right\Ver... ...rt \kern.05em \mbox{\boldmath$\space e $ } _m \kern.05em \right\Vert _{\infty }$, since all the components of $\mbox{\boldmath$\space e $ } _{m+1}$ are bounded by this, and $\left\Vert \kern.05em \mbox{\boldmath$\space e $ } _m \kern.05em \right\Vert _{... ...ern.05em \mbox{\boldmath$\space e $ } _0 \kern.05em \right\Vert _{\infty }\to 0$ as $m\to \infty$.The characteristic equation for Gauss-Seidel is

$$\begin{eqnarray*}\det [\lambda I-(D-L)^{-1}U] & = & 0\\ \Rightarrow~ \det [\lambda (D-L)-U] & = & 0, \end{eqnarray*}$$

or in displayed form

$$\left\vert \begin{array}{cccc}a_{11}\lambda & a_{12} & \cdot ... ... & a_{n,n-1}\lambda & a_{nn}\lambda \end{array} \right\vert=0.$$

We return to the case when A is positive definite, but now consider complex A, for which the definition is slightly different.

Definition 4.4.2   An Hermitian matrix A is positive definite if $\mbox{\boldmath$\space x $ } ^*A\mbox{\boldmath$\space x $ } >0$ for every $\mbox{\boldmath$\space x $ } \not= \mbox{\bf0}$.

Theorem 4.4.3 (Keller)   Let A be Hermitian positive definite and let N be any nonsingular matrix such that N+N^*-A is positive definite. Then the iteration

$$N\mbox{\boldmath$ x $ } _{m+1}=(N-A)\mbox{\boldmath$ x $ } _m+\mbox{\boldmath$ f $ }$$

is convergent.

Proof.  Let the iteration matrix, M=N^-1(N-A), have an eigenvalue $\lambda$, with corresponding eigenvector $\mbox{\boldmath$\space u $ }$, then

$$\begin{eqnarray*}N^{-1}(N-A)\mbox{\boldmath$ u $ } & = & \lambda \mbox{\boldmath... ...& (1-\lambda )\mbox{\boldmath$ u $ } ^*N\mbox{\boldmath$ u $ } . \end{eqnarray*}$$

Since A is positive definite, $\lambda \not= 1$. Hence

$$\begin{eqnarray*}\frac{1}{1-\lambda } & = & \frac{\mbox{\boldmath$ u $ } ^*N\mbo... ...}{\mbox{\boldmath$ u $ } ^*A\mbox{\boldmath$ u $ } }\\ & > & 0, \end{eqnarray*}$$

so $\vert\lambda \vert<1$, and convergence is established.

Corollary 4.4.4   The Gauss-Seidel method converges to the solution of $A\mbox{\boldmath$\space x $ } =\mbox{\boldmath$\space f $ }$ if A is Hermitian positive definite.

Proof.  For the Gauss-Seidel method, N=D-L and N-A=U, so

$$\begin{eqnarray*}N+N^*-A & = & N^*+(N-A)\\ & = & D^*-L^*+U\\ & = & D^*, \end{eqnarray*}$$

since $A^*=A~ \Rightarrow~ L^*=U$. Now $0<\mbox{\boldmath$\space e $ } _i^*A\mbox{\boldmath$\space e $ } _i=a_{ii}$since A is positive definite, so $N+N^*-A=D^*=\, {\rm diag}\, [\bar{a}_{ii}]=\, {\rm diag}\, [a_{ii}]$ and

$$\mbox{\boldmath$ x $ } ^*(N+N^*-A)\mbox{\boldmath$ x $ }=\mbo... ...*\mbox{\boldmath$ x $ }=\sum _{i=1}^na_{ii}\vert x_i\vert^2>0,$$

and so N+N^*-A is positive definite.