逐次超松弛迭代法

数值线性代数中，逐次超松弛（successive over-relaxation，SOR）迭代法是高斯-赛德尔迭代的一种变体，用于求解线性方程组。类似方法也可用于任何缓慢收敛的迭代过程。

SOR迭代法由David M. Young Jr.和Stanley P. Frankel在1950年同时独立提出，目的是在计算机上自动求解线性方程组。之前，人们已经为计算员的计算开发过超松弛法，如路易斯·弗莱·理查德森的方法以及R. V. Southwell开发的方法。但这些方法需要一定专业知识确保求解的收敛，不适用于计算机编程。David M. Young Jr.的论文对这些方面进行了探讨。^[1]

形式化

给定n个线性方程组成的方系统：

${\displaystyle A\mathbf {x} =\mathbf {b} }$

其中

${\displaystyle A={\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}},\qquad \mathbf {x} ={\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}},\qquad \mathbf {b} ={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{bmatrix}}.}$

则A可分解为对角矩阵D、严格上下三角矩阵U、L：

${\displaystyle A=D+L+U,}$

其中

${\displaystyle D={\begin{bmatrix}a_{11}&0&\cdots &0\\0&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &a_{nn}\end{bmatrix}},\quad L={\begin{bmatrix}0&0&\cdots &0\\a_{21}&0&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &0\end{bmatrix}},\quad U={\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\0&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &0\end{bmatrix}}.}$

线性方程组可重写为

${\displaystyle (D+\omega L)\mathbf {x} =\omega \mathbf {b} -[\omega U+(\omega -1)D]\mathbf {x} }$

其中${\displaystyle \omega >1}$是常数，称作松弛因子（relaxation factor）。

逐次超松弛迭代法可以通过迭代逼近x的精确解，可分析地写作

${\displaystyle \mathbf {x} ^{(k+1)}=(D+\omega L)^{-1}{\big (}\omega \mathbf {b} -[\omega U+(\omega -1)D]\mathbf {x} ^{(k)}{\big )}=L_{\omega }\mathbf {x} ^{(k)}+\mathbf {c} ,}$

其中${\displaystyle \mathbf {x} ^{(k)}}$是${\displaystyle \mathbf {x} }$的第k次迭代值，${\displaystyle \mathbf {x} ^{(k+1)}}$是${\displaystyle \mathbf {x} }$下一次迭代所得的值。 利用${\displaystyle (D+\omega L)}$的三角形，可用向前替换法依次计算${\displaystyle \mathbf {x} ^{(k+1)}}$的元素：

${\displaystyle x_{i}^{(k+1)}=(1-\omega )x_{i}^{(k)}+{\frac {\omega }{a_{ii}}}\left(b_{i}-\sum _{j<i}a_{ij}x_{j}^{(k+1)}-\sum _{j>i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n.}$

收敛性

SOR迭代法迭代矩阵${\displaystyle C_{\omega }}$的谱半径${\displaystyle \rho (C_{\omega })}$。 图表显示的是雅可比迭代矩阵的谱半径${\displaystyle \mu :=\rho (C_{\text{Jac}})}$。

松弛因子${\displaystyle \omega }$的选择并不容易，取决于系数矩阵的性质。1947年，亚历山大·马雅科维奇·奥斯特洛夫斯基证明，若A是对称正定矩阵，则${\displaystyle \forall 0<\omega <2,\ \rho (L_{\omega })<1}$。因此，迭代过程将收敛，但更高的收敛速度更有价值。

收敛速度

SOR法的收敛速度可通过分析得出。假设^[2][3]

• 松弛因子适当：${\displaystyle \omega \in (0,2)}$
• 雅可比法迭代矩阵${\displaystyle C_{\text{Jac}}:=I-D^{-1}A}$只有实特征值
• 雅可比法收敛：${\displaystyle \mu :=\rho (C_{\text{Jac}})<1}$
• 矩阵分解${\displaystyle A=D+L+U}$满足${\displaystyle \forall z\in \mathbb {C} \setminus \{0\},\ \lambda \in \mathbb {C} ,\ \operatorname {det} (\lambda D+zL+{\tfrac {1}{z}}U)=\operatorname {det} (\lambda D+L+U)}$

则收敛速度可表为

${\displaystyle \rho (C_{\omega })={\begin{cases}{\frac {1}{4}}\left(\omega \mu +{\sqrt {\omega ^{2}\mu ^{2}-4(\omega -1)}}\right)^{2}\,,&0<\omega \leq \omega _{\text{opt}}\\\omega -1\,,&\omega _{\text{opt}}<\omega <2\end{cases}}}$

最优松弛因子是

${\displaystyle \omega _{\text{opt}}:=1+\left({\frac {\mu }{1+{\sqrt {1-\mu ^{2}}}}}\right)^{2}=1+{\frac {\mu ^{2}}{4}}+O(\mu ^{3})\,.}$

特别地，${\displaystyle \omega =1}$时SOR法即退化为高斯-赛德尔迭代，有${\displaystyle \rho (C_{\omega })=\mu ^{2}=\rho (C_{\text{Jac}})^{2}}$。 对最优的${\displaystyle \omega }$，有${\displaystyle \rho (C_{\omega })={\frac {1-{\sqrt {1-\mu ^{2}}}}{1+{\sqrt {1-\mu ^{2}}}}}={\frac {\mu ^{2}}{4}}+O(\mu ^{3})}$，表明SOR法的效率约是高斯-赛德尔迭代的4倍。

最后一条假设对三对角矩阵也满足，因为${\displaystyle Z(\lambda D+L+U)Z^{-1}=\lambda D+zL+{\tfrac {1}{z}}U}$对对角阵Z，其元素${\displaystyle Z_{ii}=z^{i-1}}$、${\displaystyle \operatorname {det} (\lambda D+L+U)=\operatorname {det} (Z(\lambda D+L+U)Z^{-1})}$。

算法

由于此算法中，元素可在迭代过程中被覆盖，所以只需一个存储向量，不需要向量索引。

输入：A, b, ω
输出：φ

选择初始解φ
repeat until convergence
    for i from 1 until n do
        set σ to 0
        for j from 1 until n do
            if j ≠ i then
                set σ to σ + aij φj
            end if
        end (j-loop)
        set φi to (1 − ω)φi + ω(bi − σ) / aii
    end (i-loop)
    check if convergence is reached
end (repeat)

注意：${\displaystyle (1-\omega )\phi _{i}+{\frac {\omega }{a_{ii}}}(b_{i}-\sigma )}$也可写作${\displaystyle \phi _{i}+\omega \left({\frac {b_{i}-\sigma }{a_{ii}}}-\phi _{i}\right)}$，这样每次外层for循环可以省去一次乘法。

例子

解线性方程组

${\displaystyle {\begin{aligned}4x_{1}-x_{2}-6x_{3}+0x_{4}&=2,\\-5x_{1}-4x_{2}+10x_{3}+8x_{4}&=21,\\0x_{1}+9x_{2}+4x_{3}-2x_{4}&=-12,\\1x_{1}+0x_{2}-7x_{3}+5x_{4}&=-6.\end{aligned}}}$

择松弛因子${\displaystyle \omega =0.5}$与初始解${\displaystyle \phi =(0,0,0,0)}$。由SOR算法可得下表，在38步取得精确解(3, −2, 2, 1)。

迭代	${\displaystyle x_{1}}$	${\displaystyle x_{2}}$	${\displaystyle x_{3}}$	${\displaystyle x_{4}}$
01	0.25	−2.78125	1.6289062	0.5152344
02	1.2490234	−2.2448974	1.9687712	0.9108547
03	2.070478	−1.6696789	1.5904881	0.76172125
...	...	...	...	...
37	2.9999998	−2.0	2.0	1.0
38	3.0	−2.0	2.0	1.0

用Common Lisp的简单实现：

;; 默认浮点格式设为long-float，以确保在更大范围数字上正确运行
(setf *read-default-float-format* 'long-float)

(defparameter +MAXIMUM-NUMBER-OF-ITERATIONS+ 100
  "The number of iterations beyond which the algorithm should cease its
   operation, regardless of its current solution. A higher number of
   iterations might provide a more accurate result, but imposes higher
   performance requirements.")

(declaim (type (integer 0 *) +MAXIMUM-NUMBER-OF-ITERATIONS+))

(defun get-errors (computed-solution exact-solution)
  "For each component of the COMPUTED-SOLUTION vector, retrieves its
   error with respect to the expected EXACT-SOLUTION vector, returning a
   vector of error values.
   ---
   While both input vectors should be equal in size, this condition is
   not checked and the shortest of the twain determines the output
   vector's number of elements.
   ---
   The established formula is the following:
     Let resultVectorSize = min(computedSolution.length, exactSolution.length)
     Let resultVector     = new vector of resultVectorSize
     For i from 0 to (resultVectorSize - 1)
       resultVector[i] = exactSolution[i] - computedSolution[i]
     Return resultVector"
  (declare (type (vector number *) computed-solution))
  (declare (type (vector number *) exact-solution))
  (map '(vector number *) #'- exact-solution computed-solution))

(defun is-convergent (errors &key (error-tolerance 0.001))
  "Checks whether the convergence is reached with respect to the
   ERRORS vector which registers the discrepancy betwixt the computed
   and the exact solution vector.
   ---
   The convergence is fulfilled if and only if each absolute error
   component is less than or equal to the ERROR-TOLERANCE, that is:
   For all e in ERRORS, it holds: abs(e) <= errorTolerance."
  (declare (type (vector number *) errors))
  (declare (type number            error-tolerance))
  (flet ((error-is-acceptable (error)
          (declare (type number error))
          (<= (abs error) error-tolerance)))
    (every #'error-is-acceptable errors)))

(defun make-zero-vector (size)
  "Creates and returns a vector of the SIZE with all elements set to 0."
  (declare (type (integer 0 *) size))
  (make-array size :initial-element 0.0 :element-type 'number))

(defun successive-over-relaxation (A b omega
                                   &key (phi (make-zero-vector (length b)))
                                        (convergence-check
                                          #'(lambda (iteration phi)
                                              (declare (ignore phi))
                                              (>= iteration +MAXIMUM-NUMBER-OF-ITERATIONS+))))
  "Implements the successive over-relaxation (SOR) method, applied upon
   the linear equations defined by the matrix A and the right-hand side
   vector B, employing the relaxation factor OMEGA, returning the
   calculated solution vector.
   ---
   The first algorithm step, the choice of an initial guess PHI, is
   represented by the optional keyword parameter PHI, which defaults
   to a zero-vector of the same structure as B. If supplied, this
   vector will be destructively modified. In any case, the PHI vector
   constitutes the function's result value.
   ---
   The terminating condition is implemented by the CONVERGENCE-CHECK,
   an optional predicate
     lambda(iteration phi) => generalized-boolean
   which returns T, signifying the immediate termination, upon achieving
   convergence, or NIL, signaling continuant operation, otherwise. In
   its default configuration, the CONVERGENCE-CHECK simply abides the
   iteration's ascension to the ``+MAXIMUM-NUMBER-OF-ITERATIONS+'',
   ignoring the achieved accuracy of the vector PHI."
  (declare (type (array  number (* *)) A))
  (declare (type (vector number *)     b))
  (declare (type number                omega))
  (declare (type (vector number *)     phi))
  (declare (type (function ((integer 1 *)
                            (vector number *))
                           *)
                 convergence-check))
  (let ((n (array-dimension A 0)))
    (declare (type (integer 0 *) n))
    (loop for iteration from 1 by 1 do
      (loop for i from 0 below n by 1 do
        (let ((rho 0))
          (declare (type number rho))
          (loop for j from 0 below n by 1 do
            (when (/= j i)
              (let ((a[ij]  (aref A i j))
                    (phi[j] (aref phi j)))
                (incf rho (* a[ij] phi[j])))))
          (setf (aref phi i)
                (+ (* (- 1 omega)
                      (aref phi i))
                   (* (/ omega (aref A i i))
                      (- (aref b i) rho))))))
      (format T "~&~d. solution = ~a" iteration phi)
      ;; Check if convergence is reached.
      (when (funcall convergence-check iteration phi)
        (return))))
  (the (vector number *) phi))

;; Summon the function with the exemplary parameters.
(let ((A              (make-array (list 4 4)
                        :initial-contents
                        '((  4  -1  -6   0 )
                          ( -5  -4  10   8 )
                          (  0   9   4  -2 )
                          (  1   0  -7   5 ))))
      (b              (vector 2 21 -12 -6))
      (omega          0.5)
      (exact-solution (vector 3 -2 2 1)))
  (successive-over-relaxation
    A b omega
    :convergence-check
    #'(lambda (iteration phi)
        (declare (type (integer 0 *)     iteration))
        (declare (type (vector number *) phi))
        (let ((errors (get-errors phi exact-solution)))
          (declare (type (vector number *) errors))
          (format T "~&~d. errors   = ~a" iteration errors)
          (or (is-convergent errors :error-tolerance 0.0)
              (>= iteration +MAXIMUM-NUMBER-OF-ITERATIONS+))))))

上述伪代码的简单Python实现。

import numpy as np
from scipy import linalg

def sor_solver(A, b, omega, initial_guess, convergence_criteria):
    """
    This is an implementation of the pseudo-code provided in the Wikipedia article.
    Arguments:
        A: nxn numpy matrix.
        b: n dimensional numpy vector.
        omega: relaxation factor.
        initial_guess: An initial solution guess for the solver to start with.
        convergence_criteria: The maximum discrepancy acceptable to regard the current solution as fitting.
    Returns:
        phi: solution vector of dimension n.
    """
    step = 0
    phi = initial_guess[:]
    residual = linalg.norm(A @ phi - b)  # Initial residual
    while residual > convergence_criteria:
        for i in range(A.shape[0]):
            sigma = 0
            for j in range(A.shape[1]):
                if j != i:
                    sigma += A[i, j] * phi[j]
            phi[i] = (1 - omega) * phi[i] + (omega / A[i, i]) * (b[i] - sigma)
        residual = linalg.norm(A @ phi - b)
        step += 1
        print("Step {} Residual: {:10.6g}".format(step, residual))
    return phi

# An example case that mirrors the one in the Wikipedia article
residual_convergence = 1e-8
omega = 0.5  # Relaxation factor

A = np.array([[4, -1, -6, 0],
              [-5, -4, 10, 8],
              [0, 9, 4, -2],
              [1, 0, -7, 5]])

b = np.array([2, 21, -12, -6])

initial_guess = np.zeros(4)

phi = sor_solver(A, b, omega, initial_guess, residual_convergence)
print(phi)

对称逐次超松弛

对对称矩阵A，其中

${\displaystyle U=L^{T},\,}$

有对称逐次超松弛迭代法（SSOR）：

${\displaystyle P=\left({\frac {D}{\omega }}+L\right){\frac {\omega }{2-\omega }}D^{-1}\left({\frac {D}{\omega }}+U\right),}$

迭代法为

${\displaystyle \mathbf {x} ^{k+1}=\mathbf {x} ^{k}-\gamma ^{k}P^{-1}(A\mathbf {x} ^{k}-\mathbf {b} ),\ k\geq 0.}$

SOR与SSOR法都来自David M. Young Jr.

其他应用

主条目：理查德森外推法

任何迭代法都可应用相似技巧。若原迭代的形式为

${\displaystyle x_{n+1}=f(x_{n})}$

则可将其改为

${\displaystyle x_{n+1}^{\mathrm {SOR} }=(1-\omega )x_{n}^{\mathrm {SOR} }+\omega f(x_{n}^{\mathrm {SOR} }).}$

但若将x视作完整向量，则上述解线性方程组的公式不是这种公式的特例。此公式基础上，计算下一个向量的公式是

${\displaystyle \mathbf {x} ^{(k+1)}=(1-\omega )\mathbf {x} ^{(k)}+\omega L_{*}^{-1}(\mathbf {b} -U\mathbf {x} ^{(k)}),}$

其中${\displaystyle L_{*}=L+D}$。${\displaystyle \omega >1}$用于加快收敛速度，${\displaystyle \omega <1}$可使发散的迭代收敛或加快过调（overshoot）过程的收敛。有多种方法可根据观察到的收敛过程行为，自适应地调整松弛因子${\displaystyle \omega }$。这些方法通常只对一部分问题有效。

另见

• 雅可比法
• 置信传播
• 矩阵分裂

注释

1. Young, David M., Iterative methods for solving partial difference equations of elliptical type (PDF), PhD thesis, Harvard University, 1950-05-01 [2009-06-15]
2. Hackbusch, Wolfgang. 4.6.2. Iterative Solution of Large Sparse Systems of Equations | SpringerLink. Applied Mathematical Sciences 95. 2016. ISBN 978-3-319-28481-1. doi:10.1007/978-3-319-28483-5 （英国英语）.
3. Greenbaum, Anne. 10.1. Iterative Methods for Solving Linear Systems. Frontiers in Applied Mathematics 17. 1997. ISBN 978-0-89871-396-1. doi:10.1137/1.9781611970937 （英国英语）.

参考文献

• Abraham Berman, Robert J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, 1994, SIAM. ISBN 0-89871-321-8.
• Black, Noel. Successive Overrelaxation Method. MathWorld.
• A. Hadjidimos, Successive overrelaxation (SOR) and related methods, Journal of Computational and Applied Mathematics 123 (2000), 177–199.
• Yousef Saad, Iterative Methods for Sparse Linear Systems, 1st edition, PWS, 1996.
• Netlib's copy of "Templates for the Solution of Linear Systems", by Barrett et al.
• Richard S. Varga 2002 Matrix Iterative Analysis, Second ed. (of 1962 Prentice Hall edition), Springer-Verlag.
• David M. Young Jr. Iterative Solution of Large Linear Systems, Academic Press, 1971. (reprinted by Dover, 2003)

外部链接

• Module for the SOR Method
• Tridiagonal linear system solver based on SOR, in C++