1. Problem formulation

Successive over-relaxation (SOR) is designed to solve linear systems of the form \[ Ax=b, \] where \(A\) is typically sparse and structured, and \(b\) is known. Such systems often appear after discretizing partial differential equations (PDEs) on a grid.

1.1 Linear systems suited to SOR

SOR is most useful when \(A\) exhibits properties that support efficient iterative refinement. Common favorable cases include:

  • Sparsity, especially with a banded or stencil-like sparsity pattern.
  • Diagonal dominance or similar structural conditions that make diagonal entries representative of the system’s scale.
  • Matrices that can be split into components corresponding to “known/updated” and “unknown/not-yet-updated” portions during an iteration.

In PDE discretizations (e.g., finite difference stencils), the resulting linear system often naturally aligns with these requirements.

1.2 Iterative solution viewpoint

Rather than solving \(Ax=b\) directly, SOR constructs a sequence of approximations \(\{x^{(k)}\}\) that converges to the true solution \(x\). Each iteration uses the most recently computed components as soon as they are available, so the method progressively corrects earlier estimates within the same iteration.

This “incremental update” perspective is central: the quality of the iteration depends not only on the relaxation factor but also on how updates are ordered and how the matrix is partitioned.

1.3 Ordering and splitting of matrices

To define SOR, the matrix is implicitly partitioned according to an update order. The standard approach rewrites \(A\) into a form \[ A = D + L + U, \] where:

  • \(D\) is the diagonal part,
  • \(L\) is the strictly lower-triangular part,
  • \(U\) is the strictly upper-triangular part.

The iteration then uses values from \(L\) and \(D\) that have already been updated (or can be treated as updated) while treating \(U\) as not yet updated for the current iteration step. This ordering choice has a significant impact on convergence speed.

2. Gauss–Seidel as the baseline

SOR is commonly introduced as an acceleration of the Gauss–Seidel method. Understanding Gauss–Seidel provides both the baseline update rule and the intuition behind relaxation.

2.1 Gauss–Seidel update rule

Given \(x^{(k)}\), Gauss–Seidel computes components sequentially. For the \(i\)-th component, \[ x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum_{j<i} a_{ij}x_j^{(k+1)}-\sum_{j>i} a_{ij}x_j^{(k)}\right), \] where newly computed values \(x_j^{(k+1)}\) for \(j<i\) are used immediately, while the remaining components use the previous iterate.

This scheme can be interpreted as repeatedly applying a linear transformation derived from the decomposition \(A=D+L+U\).

2.2 Convergence intuition

Gauss–Seidel tends to converge when the coupling in the matrix does not overwhelm the diagonal contribution. Informally, each update “balances” the residual for one row using the latest information available. If the error propagation associated with the iteration shrinks over time, convergence follows.

However, plain Gauss–Seidel may converge slowly, especially for large, grid-based PDE systems where error modes decay at very different rates.

2.3 Relationship to matrix splittings

In splitting form, Gauss–Seidel can be expressed through an iteration operator constructed from \((D+L)\) and \(U\). The same viewpoint extends to SOR: SOR changes how the iteration operator acts by inserting a relaxation factor, thereby altering the error-reduction characteristics.

3. Successive over-relaxation method

SOR modifies the Gauss–Seidel update by blending the newly computed value with the previous estimate using a parameter \(\omega\), called the relaxation factor.

3.1 Relaxation factor and weighted updates

The central idea is to compute the Gauss–Seidel candidate \( \tilde{x}_i \) and then update \[ x_i^{(k+1)} = (1-\omega)x_i^{(k)} + \omega \tilde{x}_i, \] where \(0<\omega\le 2\) is typical in classical settings.

  • When \(\omega=1\), this reduces to Gauss–Seidel.
  • If \(\omega>1\), the method “over-relaxes” the update, potentially accelerating convergence.
  • If \(\omega<1\), it “under-relaxes,” often improving robustness at the cost of slower progress.

3.2 Standard SOR iteration scheme

Let \(\tilde{x}_i\) be the Gauss–Seidel value computed from the current iterate components (using updated values where available). A standard SOR scheme proceeds in the same order as Gauss–Seidel: \[ \tilde{x}_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j&lt;i} a_{ij}x_j^{(k+1)}-\sum_{j&gt;i} a_{ij}x_j^{(k)}\right), \] then \[ x_i^{(k+1)}=(1-\omega)x_i^{(k)}+\omega \tilde{x}_i. \]

In compact matrix form, SOR also corresponds to an iteration based on a modified splitting that incorporates \(\omega\).

3.3 Typical stopping criteria

Practical implementations stop when the solution appears sufficiently accurate. Common criteria include:

- Residual norm test: \(\|Ax^{(k)}-b\|\le \text{tol}\).
- Relative residual: \(\|Ax^{(k)}-b\|/\|b\|\le \text{tol}\).
- Update size test: \(\|x^{(k+1)}-x^{(k)}\|\le \text{tol}\cdot \|x^{(k+1)}\|\).

For ill-conditioned problems or large systems, residual-based tests are often preferred since they measure the error in the equation itself.

3.4 Implementation details and indexing

SOR is sensitive to implementation order. Typical considerations:

  • Index order: The \(i\)-loop order defines which components are treated as “already updated.” Changing ordering can change convergence speed.
  • Storage/access patterns: For sparse matrices, efficient traversal of nonzeros in each row matters for cost.
  • Initialization: Convergence can depend on the initial guess, though for many PDE applications the influence is secondary to matrix structure and the choice of \(\omega\).
  • Boundary treatment: For PDE-derived systems, indices associated with boundary nodes are typically fixed or eliminated so that the iteration updates only interior degrees of freedom.

4. Convergence and stability

Convergence behavior is governed by the relaxation parameter, the matrix splitting, and the structure of \(A\).

4.1 Role of the relaxation parameter ω

The parameter \(\omega\) controls how aggressively SOR corrects the current iterate. If \(\omega\) is too small, the method resembles a damped iteration and can converge slowly. If \(\omega\) is too large, the iteration may amplify errors rather than damp them, leading to divergence or oscillation.

For many classical matrix classes, there exists an interval of \(\omega\) for which convergence occurs, with an optimal value near which convergence is fastest.

4.2 Conditions for convergence

A full characterization depends on the matrix class, but typical sufficient conditions include:

  • Diagonal dominance, especially in a way that supports the underlying Gauss–Seidel convergence.
  • Symmetric positive definite structure combined with a consistent splitting.
  • Appropriate ordering that aligns the lower/upper parts with the matrix’s coupling.

When these conditions hold, SOR’s convergence can be expected and its performance can be tuned via \(\omega\).

4.3 Spectral radius interpretation

Iterative methods often admit a formulation where the error at step \(k\) is transformed by an iteration matrix \(T\). Convergence generally requires the spectral radius \(\rho(T)\) to be less than one: \[

\|e^{(k)}\|\to 0 \quad \text{if} \quad \rho(T)<1,

\] where \(e^{(k)}=x^{(k)}-x\).

SOR changes the iteration matrix through \(\omega\), so tuning \(\omega\) effectively reshapes the error propagation and can reduce \(\rho(T)\).

4.4 Effect of matrix structure (e.g., diagonally dominant cases)

Matrix structure affects both the stability of the iteration and the relative decay rates of different error modes. For diagonally dominant systems, each row update is strongly anchored by the diagonal entry \(a_{ii}\), which helps prevent the weighted update from drifting away.

Conversely, when diagonal dominance is weak or the matrix contains strong long-range coupling, the method may require more careful \(\omega\) selection or may converge very slowly.

5. Choosing the relaxation factor

Selecting \(\omega\) is crucial for SOR’s efficiency. Several approaches range from theory-driven bounds to practical tuning workflows.

5.1 Heuristic tuning approaches

Common heuristics include:

  • Start near \(\omega=1\) (Gauss–Seidel) and gradually increase \(\omega\) while monitoring residual decrease per iteration.
  • Use a small pilot run: compare residual norms after a fixed number of iterations for a set of candidate \(\omega\) values.
  • Avoid extreme values: if \(\omega\) causes residual growth or sign of oscillatory divergence, it is likely too large.

This approach is often effective when the matrix resembles well-studied PDE discretizations.

5.2 Bounds and guideline estimates

For some structured problems, theoretical estimates relate the optimal \(\omega\) to eigenvalues of relevant operators or to grid parameters. While these formulas do not always transfer directly to arbitrary sparse matrices, they provide useful starting points:

  • Optimal \(\omega\) typically increases with the problem size in many grid-based PDE settings.
  • There are often known intervals (e.g., depending on the matrix class) where convergence is expected.

When explicit spectral information is unavailable, such bounds can still guide initial choices.

5.3 Practical selection workflows

A practical workflow for engineering use is:

  1. Choose an initial guess and an ordering for the unknowns.
  2. Run a short sequence of trials for several \(\omega\) values.
  3. Select the \(\omega\) that yields the fastest residual reduction (or lowest residual after a fixed iteration budget).
  4. Reconfirm with a longer run using the chosen stopping criterion.
  5. If convergence is slow, consider reordering variables or switching to a related method (e.g., preconditioned Krylov methods or multigrid) rather than only adjusting \(\omega\).

6. Performance considerations

SOR’s runtime depends on both the number of iterations and the per-iteration cost.

6.1 Iteration count versus per-iteration cost

Each SOR iteration typically costs roughly the cost of a Gauss–Seidel sweep: one row update per unknown, requiring access to the nonzeros in that row. Increasing \(\omega\) may reduce the number of sweeps, but if it pushes the method near instability, it can increase iterations or fail to converge. The goal is to minimize total time, not just iteration count.

6.2 Impact of sparse versus dense representations

For sparse matrices, SOR can exploit efficient iteration over nonzeros in each row. For dense matrices, the per-iteration cost becomes large, and SOR often loses appeal compared with direct solvers or more advanced iterative strategies. Therefore, SOR’s typical advantage emerges in large sparse systems from discretized PDEs.

6.3 Scaling with grid size (common PDE use cases)

In grid-based problems, discretization refinement increases the number of unknowns. For simple stationary methods like SOR, iteration counts often grow as the grid becomes finer because long-wavelength error components decay slowly. This scaling behavior helps explain why SOR is frequently used as a smoother inside multigrid methods or combined with acceleration strategies for large problems.

SOR belongs to a family of iterative approaches that differ mainly in how they update components and how they damp or amplify corrections.

7.1 Under-relaxation versus over-relaxation

  • Under-relaxation (\(\omega<1\)) dampens updates, which can reduce overshoot and improve stability for certain difficult systems.
  • Over-relaxation (\(\omega>1\)) can accelerate convergence by applying a stronger correction, but it increases sensitivity to matrix structure and ordering.

When SOR exhibits oscillatory residual behavior, under-relaxing can sometimes restore stability.

7.2 Comparison to Jacobi and Gauss–Seidel

  • Jacobi uses only values from the previous iterate to update all components simultaneously. It is easy to parallelize but often converges more slowly.
  • Gauss–Seidel uses immediate updates within the iteration, typically improving convergence speed relative to Jacobi.
  • SOR generalizes Gauss–Seidel by introducing the weighted correction, aiming to further accelerate convergence.

These comparisons highlight why SOR is often viewed as a refinement of Gauss–Seidel rather than an unrelated method.

7.3 SOR as a special case within broader iterative frameworks

SOR can be interpreted within broader concepts such as:

  • stationary iterative methods derived from matrix splittings,
  • iterative refinement schemes that modify iteration operators,
  • relaxation-based smoothers in multilevel solvers.

In these contexts, SOR’s update structure remains consistent, but its role may be limited to a smoothing step rather than serving as the full solver.

8. Practical examples

Practical SOR usage is best understood through representative PDE discretizations, where the algebraic structure is predictable.

8.1 2D Poisson-type stencil problem (overview)

A common test case is the 2D Poisson equation discretized on a uniform grid using finite differences. The resulting linear system is sparse with a stencil pattern that couples each interior grid point to its immediate neighbors. For such systems:

  • the diagonal entries are typically dominant relative to off-diagonals,
  • the matrix splitting aligns well with the natural ordering of grid points,
  • SOR can significantly reduce the number of iterations compared with Gauss–Seidel.

The effectiveness depends on how interior nodes are indexed and how boundary nodes are incorporated.

8.2 Boundary-condition handling (high level)

Boundary conditions are handled by modifying the linear system before iteration. Typical strategies include:

  • Dirichlet conditions: boundary values are inserted into the right-hand side, leaving unknowns only at interior points.
  • Neumann or mixed conditions: discretization alters row coefficients near boundaries, but the unknowns remain organized into a sparse system.

The key is that the SOR iteration updates only the degrees of freedom that remain unknown in the system matrix used for iteration.

8.3 Example workflow: from setup to convergence check

A typical workflow:

  1. Assemble \(A\) and \(b\) from a discretization, ensuring consistent indexing.
  2. Choose \(\omega\) (often starting with \(\omega=1\) and then tuning).
  3. Initialize \(x^{(0)}\) (e.g., zeros or a coarse-grid estimate).
  4. For each iteration \(k\):
  • sweep through unknown indices in the chosen order,
  • compute Gauss–Seidel candidates \(\tilde{x}_i\),
  • apply the relaxation blend to obtain \(x^{(k+1)}\).
  1. Compute the residual \(r^{(k+1)}=Ax^{(k+1)}-b\).
  2. Stop when the residual norm satisfies the tolerance or when the update is sufficiently small.

This workflow is straightforward to implement and provides clear diagnostics when convergence is slow or unstable.

9. Pitfalls and troubleshooting

Even when SOR is theoretically applicable, practical issues can degrade performance or prevent convergence.

9.1 Divergence and oscillatory behavior

Symptoms include residual norms that grow over iterations or alternate between improvements and worsening. Common causes:

  • \(\omega\) chosen too large for the matrix class,
  • an ordering that worsens the effective splitting,
  • insufficient diagonal dominance or strong coupling.

Troubleshooting often starts by reducing \(\omega\) toward 1 and checking whether convergence resumes.

9.2 Sensitivity to ω and ordering

For some systems, small changes in \(\omega\) can make a large difference in convergence rate. Ordering can also affect the lower/upper triangular roles during a sweep. If results vary strongly with reindexing, it may indicate that the matrix has anisotropic coupling or that the chosen update order is poorly aligned with its structure.

A practical remedy is to try alternative orderings or use problem-informed traversal (e.g., natural grid ordering) when possible.

9.3 Numerical precision and stopping tolerance issues

Floating-point effects can matter:

  • If the tolerance is set below the achievable accuracy, iterations may stagnate.
  • Large systems can accumulate rounding errors in residual computations.
  • Using an update-based stopping criterion without checking residuals can lead to premature termination in some cases.

A robust approach is to use residual-based stopping and to ensure tolerance is appropriate for the problem’s conditioning and the arithmetic precision.

10. Summary and references

10.1 Key takeaways

  • SOR accelerates Gauss–Seidel by applying a weighted update controlled by \(\omega\).
  • Convergence and speed depend on matrix structure, update ordering, and the chosen relaxation factor.
  • Practical success requires careful stopping criteria and monitoring residual behavior, especially when tuning \(\omega\).
  • For large PDE-derived systems, SOR can be effective but may still scale poorly on its own, motivating use as a component in broader solvers.

10.2 Further reading suggestions

For deeper study, readers often consult texts on numerical linear algebra and PDE discretization, focusing on stationary iterative methods, matrix splitting theory, and practical tuning of relaxation parameters. Topics of interest include spectral analysis of iteration operators, convergence conditions for sparse structured matrices, and integration of relaxation methods into multigrid frameworks.