超松弛迭代法

超松弛迭代法（successive over-relaxation，SOR）在 Gauss–Seidel 的逐分量修正上乘一个松弛参数。参数为一时回到 GS；大于一时沿当前坐标的修正方向多走一些，小于一时少走一些。能否加速取决于矩阵及参数，过大的超调会使迭代不收敛。

一次超调的具体数值

对系统 $2x-y=1$、$-x+2y=1$，精确解为 $(1,1)$。取松弛参数 $\omega =1.2$，按先 $x$ 后 $y$ 的顺序更新： $$x^{k+1}=(1-\omega )x^k+\frac{\omega }{2}(1+y^k),y^{k+1}=(1-\omega )y^k+\frac{\omega }{2}(1+x^{k+1}).$$ 从零开始，第一轮为 $(0.6,0.96)$，第二轮为 $(1.056,1.0416)$。第二轮超过精确值并不自动意味着发散；收敛应看全局误差或残差，而不是每个分量是否一直从同一侧接近。

一般逐分量公式

对 $Ax=b$，设对角元非零。先利用新下三角、旧上三角形成当前坐标的 GS 值，再与旧值组合： $$x_i^{new}=(1-\omega )x_i^{old}+\frac{\omega }{A_{ii}}(b_i-\sum _{j<i}{A}_{ij}{x}_j^{new}-\sum _{j>i}{A}_{ij}{x}_j^{old}).$$ 若 $A=D+L+U$，迭代矩阵是 $$B_{\omega }=(D+\omega L{)}^{-1}((1-\omega )D-\omega U).$$ 它由三角求解应用，不需显式形成逆。对所有初值收敛等价于 $\rho (B_{\omega })<1$。Netlib 的 SOR 章节说明了松弛参数、谱与收敛条件。

对称正定情形的完整收敛论证

若 $A$ 实对称正定，二次能量 $\varphi (x)=x^{T}Ax/2-b^{T}x$ 严格凸且水平集有界。某个坐标更新前的残量记为 $g_i=(Ax-b{)}_i$，其中已经包含当前轮此前的更新。SOR 位移为 ${\delta }_i=-\omega g_i/A_{ii}$。精确展开能量变化： $$\varphi (x+{\delta }_i{e}_i)-\varphi (x)=g_i{\delta }_i+\frac{1}{2}{A}_{ii}{\delta }_i^2=-\frac{\omega (2-\omega )}{2A_{ii}}{g}_i^2.$$ 当 $0<\omega <2$，每次坐标更新都下降能量。能量有下界，所以各坐标所遇残量的平方总和有限，残量趋零，坐标位移也趋零。迭代有聚点；因一整轮位移趋零，聚点同时满足每个坐标的残量为零，即 $Ax=b$。正定性使该解唯一，所有聚点相同，因此整个序列收敛。

这个证明同时表明范围的作用。若 $\omega =2$，能量下降量为零。主例从零开始依次得到 $(1,2),(2,1),(0,0)$，三轮循环不收敛。更一般地，$det{B}_{\omega }=(1-\omega {)}^n$，若所有特征值模都小于一，必须有 $|1-\omega |<1$，因此实参数 $0<\omega <2$ 也是收敛的必要范围；对非正定矩阵，它不再是充分条件。

本图的最优参数只针对给定二阶矩阵。参数一对应GS；轻度超松弛可改善谱半径，继续增大反而变慢。

本例的最优参数与适用范围

主例的迭代矩阵为 $$B_{\omega }=\left(\begin{array}{cc}1-\omega & \omega /2\\ \omega (1-\omega )/2& 1-\omega +{\omega }^2/4\end{array}\right).$$ 其迹为 $T=2(1-\omega )+{\omega }^2/4$，行列式为 $(1-\omega {)}^2$。特征根从实数合并为共轭对的临界点由 $T^2=4(1-\omega {)}^2$ 给出，在 $1<\omega <2$ 内得到 $${\omega }_{\ast }=8-4\sqrt{3}=\frac{4}{2+\sqrt{3}}\approx \mathrm{1.07179677.}$$ 临界点后两根互为共轭，模为 $\omega -1$ 并随参数增加；临界点前的较大实根下降，故此处达到最小谱半径 ${\omega }_{\ast }-1\approx 0.07179677$。参数一的谱半径为 $1/4$，所以本例有明显改善。一般最优参数依赖矩阵结构，某些 Poisson 网格公式有额外排序与结构条件，不能照搬为所有 SPD 系统的公式。

有解变式取对角系统，误差每轮乘 $1-\omega$。此时最优值是一，一步到解；任何超松弛都没有帮助。它说明“越大越快”不成立，最优参数来自具体耦合结构。

可运行的正定矩阵实现

代码以 NumPy 检查实对称正定性，并使用 $0<\omega <2$；从零出发，以真实残差验收，返回近似与轮数，不修改输入。一次 Cholesky 检查成本约 $O(n^3)$，稠密每轮更新为 $O(n^2)$，矩阵存储为 $O(n^2)$。若大稀疏问题的正定性已由模型保证，通常不为检查而额外做完整稠密分解。

Python 3

import numpy as np


def sor(A, b, omega=1.2, tol=1e-12, maxiter=1000):
    A, b = np.array(A, dtype=float, copy=True), np.array(b, dtype=float, copy=True)
    if b.ndim != 1 or not b.size or A.shape != (b.size, b.size):
        raise ValueError("square matrix matching b required")
    if not np.isfinite(A).all() or not np.isfinite(b).all() or not np.array_equal(A, A.T):
        raise ValueError("finite symmetric data required")
    if not np.isfinite(omega) or not 0 < omega < 2 or not np.isfinite(tol) or tol <= 0:
        raise ValueError("0 < omega < 2 and positive tolerance required")
    if not isinstance(maxiter, int) or maxiter < 0:
        raise ValueError("nonnegative integer iteration limit required")
    np.linalg.cholesky(A)
    x = np.zeros_like(b); target = tol*max(1.0, np.linalg.norm(b, np.inf))
    for k in range(maxiter+1):
        if np.linalg.norm(b-A@x, np.inf) <= target:
            return x, k
        if k == maxiter:
            break
        old = x.copy()
        for i in range(b.size):
            gs = (b[i]-A[i, :i]@x[:i]-A[i, i+1:]@old[i+1:])/A[i, i]
            x[i] = (1-omega)*old[i]+omega*gs
        if not np.isfinite(x).all():
            raise ArithmeticError("nonfinite SOR iterate")
    raise RuntimeError("true residual did not reach tolerance")


print(sor([[2, -1], [-1, 2]], [1, 1]))  # near [1, 1]
print(sor([[2, -1], [-1, 2]], [1, 1], 8-4*np.sqrt(3)))