Successive Over-Relaxation (SOR) Method#
The Successive Over-Relaxation (SOR) method, denoted as SOR(\(\omega\)), is a highly effective stationary iterative method that improves upon the Gauss-Seidel iteration by introducing a relaxation parameter, \(\omega\), designed to accelerate convergence.
The motivation behind SOR is to overcome the limitation of Jacobi and Gauss-Seidel methods, which lack a parameter that can be adjusted to accelerate convergence in challenging cases. The SOR iteration is equivalent to the Gauss-Seidel method when \(\omega = 1\).
Defining the SOR Iteration#
The SOR method is a splitting method that modifies the successive update generated by the Gauss-Seidel iteration (\(\mathbf{x}^{(k+1)}_{\text{GS}}\)) by taking a weighted average of the current iterate (\(\mathbf{x}^{(k)}_{\text{SOR}}\)) and the Gauss-Seidel increment (\(\Delta \mathbf{x}_{\text{GS}}\)).
If we define the Gauss-Seidel update as \(\mathbf{x}^{(k+1)}_{\text{GS}} = \mathbf{x}^{(k)} + \Delta \mathbf{x}_{\text{GS}}\), then the SOR update is:
In terms of the matrix components \(A = D - L - U\) (where \(D\) is diagonal, \(-L\) is strictly lower triangular, and \(-U\) is strictly upper triangular), the Gauss-Seidel update is:
Substituting this into the relaxed form gives the defining equation for the SOR update:
SOR as a Splitting Method#
The rearranged matrix form shows that SOR fits the general splitting framework \(\mathbf{x}^{(k+1)} = M^{-1} N \mathbf{x}^{(k)} + M^{-1}\mathbf{b}\). By separating the \(\mathbf{x}^{(k+1)}\) terms to the left-hand side, we get:
Thus, the splitting \(A = M - N\) is defined by:
The SOR iteration matrix, \(G_{\text{SOR}} = L_{\omega}\), is given by:
The Relaxation Parameter \(\omega\)#
The parameter \(\omega\) distinguishes SOR into three regimes:
Gauss-Seidel: When \(\omega = 1\), the SOR splitting recovers the Gauss-Seidel method.
Underrelaxation (\(\omega < 1\)): Used to damp the increment (\(\Delta \mathbf{x}\)). If the standard iteration is diverging, applying a smaller increment may stabilize the process.
Overrelaxation (\(\omega > 1\)): This is the typical use case, where the goal is to accelerate convergence by moving further in a potentially good direction.
The choice of \(\omega\) depends heavily on the properties of the matrix \(A\).
Python Code Example#
Here is a simple implementation of the SOR method in Python:
import numpy as np
def sor(A, b, omega, x0=None, tol=1e-10, max_iter=1000):
"""
Solve the linear system Ax = b using the SOR iteration method.
Parameters:
-----------
A : ndarray
Coefficient matrix (n x n)
b : ndarray
Right-hand side vector (n,)
omega : float
Relaxation parameter (0 < omega < 2)
x0 : ndarray, optional
Initial guess (n,). If None, uses zero vector.
tol : float, optional
Tolerance for convergence (default: 1e-10)
max_iter : int, optional
Maximum number of iterations (default: 1000)
Returns:
--------
x : ndarray
Approximate solution vector
iterations : int
Number of iterations performed
"""
n = len(b)
x = np.zeros(n) if x0 is None else x0.copy()
for k in range(max_iter):
x_old = x.copy()
# Update each component using SOR iteration
for i in range(n):
sigma = np.dot(A[i, :i], x[:i]) + np.dot(A[i, i+1:], x[i+1:])
x[i] = (1 - omega) * x[i] + omega * (b[i] - sigma) / A[i, i]
# Check convergence
if np.linalg.norm(x - x_old, ord=np.inf) < tol:
return x, k + 1
print(f"Warning: Maximum iterations ({max_iter}) reached without convergence")
return x, max_iter
Convergence Conditions#
Theorem 33 (SOR Convergence Conditions)
The SOR method converges if and only if the spectral radius of the iteration matrix \(\rho(G_{\text{SOR}})\) is less than one.
Necessary Condition#
Theorem 34 (Necessary Condition for SOR Convergence)
A necessary condition for the convergence of the SOR method is \(0 < \omega < 2\).
Proof. This can be proven by analyzing the determinant of the iteration matrix. Since \(\det(G_{\text{SOR}}) = (1 - \omega)^n\), and the spectral radius must satisfy
convergence requires \(|\omega - 1| < 1\), which yields \(0 < \omega < 2\).
Sufficient Condition#
The convergence of SOR is guaranteed under certain conditions on the matrix \(A\).
Theorem 35 (Sufficient Condition for SOR Convergence)
If \(A\) is a Symmetric Positive Definite matrix, SOR converges for any relaxation parameter \(\omega\) such that \(0 < \omega < 2\). This also implies that Gauss-Seidel (\(\omega=1\)) converges if \(A\) is SPD.
Orderings#
The selection of the ordering, or permutation, of variables is a fundamental concern for the Successive Over-Relaxation (SOR) method, as this iterative approach, like Gauss-Seidel, relies on computing components successively based on a fixed prescribed order. The sequence in which components are updated influences both the rate of convergence and the potential for parallelism.
We can identify several common types of orderings used with SOR:
1. Sequential (Natural) Ordering: 1 to \(n\)#
This is the standard, fixed, prescribed order in which the entries of the solution vector \(\mathbf{x}^{(k+1)}\) are computed. This ordering corresponds to iterating through the indices from \(j=1\) up to \(j=n\) in sequence.
In the context of matrices derived from partial differential equations (PDEs), common sequential orderings (such as row-wise or column-wise enumeration of grid points) are often referred to as natural orderings. These typically represent consistent orderings.
For matrices derived from convection-diffusion problems, optimal convergence speed using SOR often depends critically on ordering the variables in the direction of the convection flow.
2. Symmetric SOR (SSOR): 1 to \(n\) followed by \(n\) to 1#
The Symmetric Successive Over-Relaxation (SSOR) method uses a pair of sweeps: a forward sweep (1 to \(n\)) followed by a backward sweep (\(n\) to 1).
An SSOR step consists of the standard SOR step, followed by a backward SOR step. This backward sweep is defined analogously to the backward Gauss-Seidel sweep, where coordinate corrections are performed in the reverse order, \(n, n-1, \dots, 1\).
This symmetrization is crucial because, if the matrix \(A\) is symmetric and positive definite (SPD), the iteration matrix resulting from the forward/backward sweeps is similar to an SPD matrix. This property makes SSOR highly effective as a preconditioner in iterative acceleration methods, particularly the Conjugate Gradient (CG) method.
3. Randomized Permutation#
The SOR iteration can also be implemented using a different randomized permutation at every iteration.
4. Alternative Structural Orderings#
Beyond simple sequential or randomized orderings, specific structural permutations can be vital for accelerating SOR, especially when solving systems derived from PDEs. These include red-black and multicolor orderings, which are techniques that partition the variables (or nodes in the associated graph) into independent sets such that no two variables in the same set are directly coupled in the matrix.
For example, in a 2D grid corresponding to a discretized PDE, a red-black ordering colors the grid points alternately (like a chessboard), so that all “red” points can be updated independently of each other, followed by all “black” points. This allows for greater parallelism and can improve convergence rates. Multicolor orderings generalize this idea to more than two colors, further dividing the variables into multiple independent sets, which is especially useful for matrices with more complex sparsity patterns. By updating all variables of the same color simultaneously, these orderings exploit the matrix structure to enhance the efficiency and scalability of the SOR method.
Optimality and Performance#
The primary difficulty in implementing SOR is determining the optimal relaxation parameter, \(\omega_{\text{opt}}\), which minimizes \(\rho(L_{\omega})\). When \(\omega_{\text{opt}}\) is used, SOR can significantly outperform Jacobi and Gauss-Seidel. The SOR method is a classic example of achieving substantial acceleration by leveraging knowledge about the matrix structure to tune a parameter, transforming a linearly convergent method into one with potentially a faster convergence rate.