# Sort of problems where SOR is faster than Gauss-Seidel?

Is there any simple rule of thumb to say if it is worth to do SOR instead of Gauss-Seidel? ( and possible way how to estimate realxation parameter $\omega$)

I mean just by looking on the matrix, or knowledge of particular problem the matrix represents?

I was reading answer on this questions: Are there any heuristics for optimizing the successive over-relaxation (SOR) method? but it's a bit too sophisticated. I don't see simple heuristics how to estimate spectral radius just looking on the matrix ( or problem which it represents ).

I would like something much simpler - just few examples of matrices (problems) for which SOR converge faster.

I was experimenting with SOR for matrix of this king: $A = I + C + R$ where $I$ is identity matrix, $C_{ij}=c$ $\forall i,j$ and $R_{ij}$s are random numbers from unifrom distribution such that $|R_{ij}|<r$. I was thinking that there would be some dependence of optimal $\omega$ on parameters $c,r$.

EDIT: I used very small $c,r$ to make sure tha $A$ is strongly diagonally dominant. ( $|c|<0.1$, $r<2|c|$ for matrix of dimension 5-10 ). I should also say that these $A$ was real and symmetric.

However, I found that Gauss-Seidel ($\omega=1$) is almost always the best (?). Does this means that there has to be some more correlation between $A_{ij}$s to get advantage of SOR ? Or I did something wrong?

I know, that SOR is not the most efficient solver (in comparison to CG,GMRES ... ) but it is simple to implement and paraelize, and modify for particular problem. Sure good for prototyping.

The convergence of classical iterative solvers for linear systems is determined by the spectral radius of the iteration matrix, $\rho(\mathbf{G})$. For a general linear system, it is difficult to determine an optimal (or even good) SOR parameter due to the difficulty in determining the spectral radius of the iteration matrix. Below I have included many additional details, including an example of a real problem where the optimal SOR weight is known.

The spectral radius is defined as the absolute value of the largest magnitude eigenvalue. A method will converge if $\rho<1$ and a smaller spectral radius means faster convergence. SOR works by altering the matrix splitting used to derive the iteration matrix based on the choice of a weighting parameter $\omega$, hopefully decreasing the spectral radius of the resulting iteration matrix.

# Matrix splitting

For the discussion below, I will assume that the system to be solved is given by

$$\mathbf{A}x=b,$$

with an iteration of the form

$$x^{(k+1)}=v+\mathbf{G}x^{(k)},$$

where $v$ is a vector, and iteration number $k$ is denoted $x^{(k)}$.

SOR takes a weighted average of the old iteration and a Gauss-Seidel iteration. The Gauss-Seidel method relies on a matrix splitting of the form

$$\mathbf{A} = \mathbf{D} + \mathbf{L} + \mathbf{U}$$

where $\mathbf{D}$ is the diagonal of $\mathbf{A}$, $\mathbf{L}$ is a lower triangular matrix containing all elements of $\mathbf{A}$ strictly below the diagonal and $\mathbf{R}$ is an upper triangular matrix containing all elements of $\mathbf{A}$ strictly above the diagonal. The Gauss-Seidel iteration is then given by

$$x^{(k+1)}=(\mathbf{D}+\mathbf{L})^{-1}b+\mathbf{G}_{\rm G-S}x^{(k)}$$

and the iteration matrix is

$$\mathbf{G}_{\rm G-S}=-(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U}.$$

SOR can then be written as

$$x^{(k+1)}=\omega(\mathbf{D}+\omega\mathbf{L})^{-1}b+\mathbf{G}_{\rm SOR}x^{(k)}$$

where

$$\mathbf{G}_{\rm SOR} = (\mathbf{D}+\omega\mathbf{L})^{-1}((1-\omega)\mathbf{D}-\omega\mathbf{U}).$$

Determining the convergence rate of the iterative scheme really boils down to determining the spectral radius of these iteration matrices. In general, this is a hard problem unless you know something specific about the structure of the matrix. There are very few examples that I'm aware of where the optimal weighting coefficient is computable. In practice, $\omega$ must be determined on the fly based on the observed (presumed) convergence of the running algorithm. This works in some cases, but fails in others.

# Optimal SOR

One realistic example where the optimal weighting coefficient is known arises in the context of solving a Poisson equation:

$$\nabla^2 u = f~{\rm in}~\Omega\\ u = g~{\rm on}~\partial\Omega$$

Discretizing this system on a square domain in 2D using second order finite differences with uniform grid spacing results in a symmetric banded matrix with 4 on the diagonal, -1 immediately above and below the diagonal, and two more bands of -1 some distance from the diagonal. There are some differences due to boundary conditions, but that is the basic structure. Given this matrix, the provably optimal choice for SOR coefficient is given by

$$\omega = \dfrac{2}{1+\sin(\pi \Delta x/L)}$$

where $\Delta x$ is the grid spacing and $L$ is the domain size. Doing so for a simple case with a known solution gives the following error versus iteration number for these two methods:

As you can see, SOR reaches machine precision in about 100 iterations at which point Gauss-Seidel is about 25 orders of magnitude worse. If you want to play around with this example, I've included the MATLAB code I used below.

clear all
close all

%number of iterations:
niter = 150;

%number of grid points in each direction
N = 16;
% [x y] = ndgrid(linspace(0,1,N),linspace(0,1,N));
[x y] = ndgrid(linspace(-pi,pi,N),linspace(-pi,pi,N));
dx = x(2,1)-x(1,1);
L = x(N,1)-x(1,1);

%desired solution:
U = sin(x/2).*cos(y);

% Right hand side for the Poisson equation (computed from U to produce the
% desired known solution)
Ix = 2:N-1;
Iy = 2:N-1;
f = zeros(size(U));
f(Ix,Iy) = (-4*U(Ix,Iy)+U(Ix-1,Iy)+U(Ix+1,Iy)+U(Ix,Iy-1)+U(Ix,Iy+1));

figure(1)
clf
contourf(x,y,U,50,'linestyle','none')
title('True solution')

%initial guess (must match boundary conditions)
U0 = U;
U0(Ix,Iy) = rand(N-2);

%Gauss-Seidel iteration:
UGS = U0; EGS = zeros(1,niter);
for iter=1:niter
for iy=2:N-1
for ix=2:N-1
UGS(ix,iy) = -1/4*(f(ix,iy)-UGS(ix-1,iy)-UGS(ix+1,iy)-UGS(ix,iy-1)-UGS(ix,iy+1));
end
end

%error:
EGS(iter) = sum(sum((U-UGS).^2))/sum(sum(U.^2));
end

figure(2)
clf
contourf(x,y,UGS,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow

%SOR iteration:
USOR = U0; ESOR = zeros(1,niter);
w = 2/(1+sin(pi*dx/L));
for iter=1:niter
for iy=2:N-1
for ix=2:N-1
USOR(ix,iy) = (1-w)*USOR(ix,iy)-w/4*(f(ix,iy)-USOR(ix-1,iy)-USOR(ix+1,iy)-USOR(ix,iy-1)-USOR(ix,iy+1));
end
end

%error:
ESOR(iter) = sum(sum((U-USOR).^2))/sum(sum(U.^2));
end

figure(4)
clf
contourf(x,y,USOR,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow

figure(5)
clf
semilogy(EGS,'b')
hold on
semilogy(ESOR,'r')
title('L2 relative error')
xlabel('Iteration number')
legend('Gauss-Seidel','SOR','location','southwest')

• Do you know of any good/well-known techniques that are used to compute the SOR parameter on the fly? I have heard before that these techniques use estimates of the spectral radius -- could you explain how they use the spectral radius, or provide a good reference? – nukeguy Jul 21 '17 at 16:54
• Oh, I see that this is addressed in the linked question scicomp.stackexchange.com/questions/851/… . Never mind my questions, but if you have more to add, please feel free to do so. – nukeguy Jul 21 '17 at 16:56
• @Doug Lipinski I thought that f should be multiplied by dx*dy. This factor coming from the discrete second derivative (see here for example). Btw, when I do it the algorithm does not work properly. Do you know why? – shamalaia Mar 5 '18 at 1:05

This side of things isn't really my speciality, but I don't think this is a super-fair test for many realistic applications.

I'm not sure what values you were using for c and r, but I suspect you were working with extremely ill-conditioned matrices. (Below is some Python code showing that these might not be the most invertible matrices.)

>>> import numpy
>>> for n in [100, 1000]:
...     for c in [.1, 1, 10]:
...         for r in [.1, 1, 10]:
...             print numpy.linalg.cond(
...                 numpy.eye(n) +
...                 c * numpy.ones((n, n)) +
...                 r * numpy.random.random((n, n))
...             )
...
25.491634739
2034.47889101
2016.33059429
168.220149133
27340.0090644
5532.81258852
1617.33518781
42490.4410689
5326.3865534
6212.01580004
91910.8386417
50543.9269739
24737.6648458
271579.469212
208913.592289
275153.967337
17021788.5576
117365.924601


If you actually needed to invert matrices this ill-conditioned, you'd a) use a specialized method, and b) should probably just go find a new field 😉

For well-conditioned matrices of any size, SOR is likely to be faster. For real problems where speed matters, it would be rare to use SOR -- on the sophisticated side, there's much better these days; on the slow but reliable side, SOR isn't the best you can do.

• hi, I don't say that my "test" is fair. I would not even say that it is a test, it's just my naieve attempt to get some idea how SOR and Gauss-Seidel behaves experimentaly. Assume that I'm a complete noob in this field. My parameters was in range $0.01<|c|<0.1$ and $r<2|c|$. To make sure that the matrix is strongly diagonally domiannt (I used smaller matrices of dimension ~10 ) – Prokop Hapala Feb 15 '14 at 22:05
• I was going to say strongly diagonally dominant. – meawoppl Apr 17 '14 at 15:57

OK, so for symetric matrices of this king:

1 t 0 0 0 0 t t 0 0
t 1 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 t 0 t
0 0 0 1 0 0 0 0 0 t
0 0 0 0 1 t 0 0 0 0
0 0 0 0 t 1 0 t 0 0
t 0 0 0 0 0 1 0 0 0
t 0 t 0 0 t 0 1 0 0
0 0 0 0 0 0 0 0 1 t
0 0 t t 0 0 0 0 t 1


SOR does converge faster than Gauss-Seidel if number of $t$s in each row is small ( much smaller than dimension of A ) and if all $t$s are similar. I was using $t$s generated like this:

$t_i = c + random(-r,r)$

If $t$s vary much and are centred arround 0 ( $c=0,r=0.1$ ) than Gauss-Seidel is faster. Gauss-Seidel is also faster if each row is more than half filled by $t$s. This also means that SOR is better for very big and very sparse matrices.

( This is just emperical observation, nothing rigorous )