# Are there any heuristics for optimizing the successive over-relaxation (SOR) method?

As I understand it, successive over relaxation works by choosing a parameter $0\leq\omega\leq2$ and using a linear combination of a (quasi) Gauss-Seidel iteration and the value at the previous timestep... that is

${u}^{k+1} = (\omega){u_{gs}}^{k+1} + (1-\omega)u^{k}$

I state 'quasi' because ${u_{gs}}^{k+1}$ includes the latest information updated according to this rule, at any timestep. (note that at $\omega=1$, this is exactly gauss-seidel).

In any case, I have read that on optimal choice for $\omega$ (such that the iteration converges faster than any other) approaches 2 for the poisson problem as the spatial resolution approaches zero. Does a similar trend exist for other symmetric, diagonally dominant problems? That is, is there a way to choose omega optimally without embedding it into an adaptive optimization scheme? Are there other heuristics for other types of problems? What kinds of problems would under-relaxation ($\omega<1$) be optimal?

-
Not quite your question, but see Salakhutdinov and Roweis, Adaptive Overrelaxed Bound Optimization Methods 2003, 8p. (Adaptive speedups have high bang per buck, but are afaik impossible to analyze, so off-topic here.) –  denis Dec 16 '14 at 17:23

### Damped Jacobi

Suppose the matrix $A$ has diagonal $D$. If the spectrum of $D^{-1}A$ lies in the interval $[a,b]$ of the positive real axis, then the iteration matrix of Jacobi with damping factor $\omega$ $$B_\text{Jacobi} = I - \omega D^{-1} A$$ has spectrum in the range $[1 - \omega b,1 - \omega a]$, so minimizing the spectral radius with $$\omega_{\text{opt}} = \frac 2 {a + b}$$ gives a convergence factor of $$\rho_\text{opt} = 1 - \frac{2a}{a+b} = \frac{b-a}{a+b}.$$ If $a \ll b$, then this convergence factor is very poor, as expected. Note that it is relatively easy to estimate $b$ using a Krylov method, but quite expensive to estimate $a$.

### Successive over-relaxation (SOR)

Young (1950) proved an optimal result for SOR applied to matrices with so-called Property A, consistent ordering, and positive real eigenvalues of $D^{-1}A$. Given a maximal eigenvalue $\mu_\max$ of the undamped Jacobi iteration matrix $I - D^{-1} A$ ($\mu_\max < 1$ is guaranteed by the assumptions in this case), the optimal damping factor for SOR is $$\omega_\text{opt} = 1 + \left( \frac{\mu_\max}{1 + \sqrt{1 - \mu_\max^2}} \right)^2$$ which results in a convergence rate of $$\rho_\text{opt} = \omega_\text{opt} - 1.$$ Note that $\omega_\text{opt}$ approaches 2 when $\mu_\max \to 1$.

It's not 1950 any more and it really doesn't make sense to use stationary iterative methods as solvers. Instead, we use them as smoothers for multigrid. In this context, we only care to target the upper end of the spectrum. Optimizing the relaxation factor in SOR causes SOR to produce very little damping of high frequencies (in exchange for better convergence on lower frequencies), so it is usually better to use standard Gauss-Seidel, corresponding to $\omega = 1$ in SOR. For nonsymmetric problems and problems with highly variable coefficients, under-relaxed SOR ($\omega <1$) may have better damping properties.
Estimating both eigenvalues of $D^{-1}A$ is expensive, but the largest eigenvalue can be estimated quickly using a few Krylov iterations. Polynomial smoothers (preconditioned with Jacobi) are more effective than multiple iterations of damped Jacobi and are easier to configure, so they should be preferred. See this answer for more on polynomial smoothers.
It is sometimes claimed that SOR should not be used as a preconditioner for Krylov methods such as GMRES. This comes from the observation that the optimal relaxation parameter should place all eigenvalues of the iteration matrix $$B_\text{SOR} = 1 - \left(\frac 1 \omega D + L\right)^{-1} A$$ on a circle centered at the origin. The spectrum of the preconditioned operator $(\frac 1 \omega D + L)^{-1} A$ has eigenvalues on a circle of the same radius, but centered on 1. For poorly conditioned operators, the radius of the circle is quite close to 1, so GMRES sees eigenvalues close to the origin at a range of angles, which is usually not good for convergence. In practice, GMRES may converge reasonably when preconditioned with SOR, especially for problems that are already fairly well conditioned, but other preconditioners are often more effective.
The use of $p$- and $hp$-multigrid is of course independent of whether a Krylov method is used on the outside. The relative costs of various operations are of course different for GPUs compared to CPUs, and there is variability between implementations. Preconditioned Richardson is just a defect correction method. So are the Newton and Picard (if written as such) nonlinear methods. Other nonlinear methods (NGMRES, BFGS, etc) also use history, and can be better depending on the relative strength of the nonlinearity. –  Jed Brown Jan 17 '12 at 21:54