# 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

## 1 Answer

### 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$.

### Comments

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.

-
I agree it is not 1950 anymore :o), however, I disagree that it does not make sense to use stationery iterative solvers anymore. We can achieve multigrid textbook efficiency using a stationary iterative solver for an engineering application solver based on high-order nonlinear free surface solvers (both potential flow and euler equations). The efficiency was just as good as a preconditioned GMRES krylov subspace method within attainable accuracy (our recent pub is found here onlinelibrary.wiley.com/doi/10.1002/fld.2675/abstract serving as proof-of-concept). –  Allan P. Engsig-Karup Jan 17 '12 at 18:58
You are using Gauss-Seidel as a smoother for multigrid (which is where methods like SOR belong). If multigrid is performing well, an outer Krylov method is also not necessary (though your paper doesn't show those comparisons). As soon as multigrid starts to lose efficiency (e.g. more than 5 iterations to reach discretization error), it is usually worthwhile to wrap a Krylov method around the multigrid cycle. –  Jed Brown Jan 17 '12 at 19:21
The whole method is a p-multigrid with GS type smoothing, however, the complete method can be written as a stationary iterative method since all operators are constant. You can view it as a preconditioned Richardson method with M a preconditioner constructed from the multigrid method. Analysis have been done but is not yet published. Actually, this work went in the other direction that you propose. The krylov method in this work (a GMRES) was discarded and then it was turned into a high-order multigrid method as we found that this was just as efficient (and With reduced memory requirements). –  Allan P. Engsig-Karup Jan 17 '12 at 21:38
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
Note that in multigrid smoothers, it is sometimes preferable (architecture permitting), to make the high-order/low-order coupling multiplicative. This also extends the "preconditioned Richardson" formulation. (I had a discussion at a conference last week with a guy who wanted to view essentially all methods as preconditioned Richardson with nested iteration, which I don't think is particular benefit over other statements of solver composition. I don't know if it's relevant to you, but your points reminded me of the discussion.) –  Jed Brown Jan 17 '12 at 22:00