Optimum $\omega$

For some classes of matrix we are able to choose $\omega$ to minimise the spectral radius of M_GS. One such class is block-tridiagonal matrices, whose diagonal blocks are themselves diagonal, matrices which arise naturally in the discretization of PDE's. However, we restrict ourselves here to the case when A is just tridiagonal; it turns out that the optimum $\omega$ has precisely the same form.

The trick is to scale the rows and columns of the characteristic polynomial for SOR so that the off-diagonal elements are the same as those for the Jacobi case. We row scale the determinant by premultipling its matrix by a matrix $E=\, {\rm diag}\, [e_i]$ say, and column scale it by postmultiplying by $F=\, {\rm diag}\, [f_i]$ say to obtain

$$\det \{ E[(\lambda -1+\omega)D-\lambda \omega L-\omega U ]F\} =\det [(\lambda -1+\omega)EDF-\lambda \omega ELF-\omega EUF ].$$ (4.9)

The corresponding Jacobi polynomial is

$$\det [\lambda D-L-U].$$ (4.10)

We must choose E and F so that

$$\lambda \omega ELF=L\mbox{ and }\omega EUF=U.$$ (4.11)

If we can make this choice in such a way that

$$(\lambda -1+\omega)EDF=kD,$$ (4.12)

where k is independent of i, any value $\lambda _J$ of $\lambda$which makes (7.10) vanish must be a value of k which makes (7.9) vanish.

Thus, in component form, we have

$$\begin{eqnarray*}\lambda \omega e_ia_{i,i-1}f_{i-1} & = & a_{i,i-1},\quad i=2,3,... ...omega)e_ia_{ii}f_i & = & \lambda _Ja_{ii},\quad i=1,2,\ldots ,n, \end{eqnarray*}$$

where k has been replaced by $\lambda _J$ in the last of these. These equations simplify to

$$\begin{eqnarray*}e_if_{i-1} & = & \frac{1}{\lambda \omega }\\ e_if_{i+1} & = & ... ...{\omega }\\ e_if_i & = & \frac{\lambda _J}{\lambda -1+\omega }. \end{eqnarray*}$$

Dividing the last by the first and the second by the last, we obtain

$$\frac{f_i}{f_{i-1}}=\frac{\lambda \omega \lambda _J}{\lambda ... ...ac{f_{i+1}}{f_i}=\frac{\lambda -1+\omega }{\omega \lambda _J}.$$

But these must hold for $i=2,3,\ldots ,n-1$ so

$$\frac{\lambda \omega \lambda _J}{\lambda -1+\omega }=\frac{\lambda -1+\omega }{\omega \lambda _J},$$

giving

$$\lambda \omega ^2\lambda _J^2=(\lambda -1+\omega )^2.$$ (4.13)

This is known as Young's formula after its inventor. It gives the SOR eigenvalues in terms of the Jacobi eigenvalues for the corresponding coefficient matrix.

Note 1. $\omega =1$ gives back Gauss-Seidel, and Young's formula becomes $\lambda _{GS}=\lambda _J^2$, so that $\rho (M_{GS})=\rho ^2(M_J)$, and Gauss-Seidel converges twice as quickly as Jacobi.

Note 2. Writing Young's formula as a quadratic in $\sqrt{\lambda }$we obtain

$$(\sqrt{\lambda })^2-\omega \lambda _J\sqrt{\lambda } +\omega -1=0.$$ (4.14)

For convergence of SOR we require $\vert\sqrt{\lambda }\vert<1$. But the product of the two roots for $\sqrt{\lambda }$ in this quadratic is $\omega -1$, so we must have $\vert\omega -1\vert<1$, or $0<\omega <2$ for convergence.

In order to find the optimum value of $\omega$ we solve (7.14) with fixed $\lambda _J$ satisfying $0<\lambda _J^2<1$, that is, for the case when $\lambda _J$ is real and the Jacobi method is convergent. After some manipulation(!), we find

$$\frac{2\sqrt{\lambda }}{\lambda _J}=\omega \pm \sqrt{(\omega -\omega _1)(\omega -\omega _2)},$$ (4.15)

where

$$\omega _1=\frac{2}{1+\sqrt{1-\lambda _J^2}},\quad \omega _2=\frac{2}{1-\sqrt{1-\lambda _J^2}},$$

for which $0<\omega _1<\omega _2$.

If $0\leq \omega \leq \omega _1, ~ \lambda$ is real, and since we are interested in $\rho (M_{SOR})$ we need the larger $\vert\lambda \vert$, so we consider

$$\begin{eqnarray*}\frac{2\vert\sqrt{\lambda }\vert}{\vert\lambda _J\vert} & = & \... ...ega _1)^2}{4(\omega -\omega _1)(\omega -\omega _2)}}\\ & < & 0, \end{eqnarray*}$$

so $\vert\sqrt{\lambda }\vert$ is a decreasing function of $\omega$.

If $\omega _1<\omega \leq \omega _2$, $\lambda$ is complex, and $\vert\lambda \vert=\omega -1>0$, since $\omega _1>1$. Thus we can sketch the variation of $\vert\lambda \vert$ with the parameter $\omega$: Clearly, $\vert\lambda \vert$ is minimal for each eigenvalue $\lambda _J$ when

$$\omega =\omega _1=\frac{2}{1+\sqrt{1-\lambda _J^2}},$$

with the optimum $\vert\lambda \vert _{\mbox{opt}}=\omega _1-1$. As $\lambda _J$varies (i.e., as we run through the eigenvalues of M_J) the largest value of $\vert\lambda \vert _{\mbox{opt}}$ occurs for the largest value of $\omega _1$, which itself occurs for the largest value of $\lambda _J$, so that

$$\omega _{\mbox{opt}}=\frac{2}{1+\sqrt{1-\rho ^2(M_J)}},\quad \rho _{\mbox{opt}}=\omega _{\mbox{opt}}-1.$$

For example, if $\rho (M_J)=0.8,~ \omega _{\mbox{opt}}=1.25,~ \rho _{\mbox{opt}}=0.25,~ \rho _{GS}=0.64$. [filename: d306/nla.tex]