逐次超松弛迭代法

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

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

形式化

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

$A𝐱=𝐛$

其中

$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right],𝐱=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right],𝐛=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right].$

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

$A=D+L+U,$

其中

$D=\left[\begin{array}{cccc}{a}_{11}& 0& \cdots & 0\\ 0& a_{22}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & a_{nn}\end{array}\right],L=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ a_{21}& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & 0\end{array}\right],U=\left[\begin{array}{cccc}0& a_{12}& \cdots & a_{1n}\\ 0& 0& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0\end{array}\right].$

线性方程组可重写为

$(D+\omega L)𝐱=\omega 𝐛-[\omega U+(\omega -1)D]𝐱$

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

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

${𝐱}^{(k+1)}=(D+\omega L{)}^{-1}(\omega 𝐛-[\omega U+(\omega -1)D]{𝐱}^{(k)})=L_{\omega }{𝐱}^{(k)}+𝐜,$

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

$x_i^{(k+1)}=(1-\omega )x_i^{(k)}+\frac{\omega }{a_{ii}}(b_i-\sum _{j<i}{a}_{ij}{x}_j^{(k+1)}-\sum _{j>i}{a}_{ij}{x}_j^{(k)}),i=1,2,\dots ,n.$

收敛性

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

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

收敛速度

SOR法的收敛速度可通过分析得出。假设

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

则收敛速度可表为

$\rho (C_{\omega })=\{\begin{array}{cc}\frac{1}{4}{(\omega \mu +\sqrt{{\omega }^2{\mu }^2-4(\omega -1)})}^2,& 0<\omega \le {\omega }_{\text{opt}}\\ \omega -1,& {\omega }_{\text{opt}}<\omega <2\end{array}$

最优松弛因子是

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

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

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

算法

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

输入：Template:Mvar, Template:Mvar, Template:Mvar
输出：Template:Mvar

选择初始解Template:Mvar
repeat until convergence
    for Template:Mvar from 1 until Template:Mvar do
        set Template:Mvar to 0
        for Template:Mvar from 1 until Template:Mvar do
            if Template:Mvar ≠ Template:Mvar then
                set Template:Mvar to Template:Math
            end if
        end (Template:Mvar-loop)
        set Template:Math to Template:Math
    end (Template:Mvar-loop)
    check if convergence is reached
end (repeat)

注意：$(1-\omega ){\varphi }_i+\frac{\omega }{a_{ii}}(b_i-\sigma )$也可写作${\varphi }_i+\omega (\frac{b_i-\sigma }{a_{ii}}-{\varphi }_i)$，这样每次外层for循环可以省去一次乘法。

例子

解线性方程组

$\begin{array}{cc}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{array}$

择松弛因子$\omega =0.5$与初始解$\varphi =(0,0,0,0)$。由SOR算法可得下表，在38步取得精确解Template:Nowrap。

迭代	$x_1$	$x_2$	$x_3$	$x_4$
Template:01	0.25	−2.78125	1.6289062	0.5152344
Template:02	1.2490234	−2.2448974	1.9687712	0.9108547
Template: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，其中

$U=L^T,$

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

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

迭代法为

${𝐱}^{k+1}={𝐱}^k-{\gamma }^k{P}^{-1}(A{𝐱}^k-𝐛),k\ge 0.$

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

其他应用

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

$x_{n+1}=f(x_n)$

则可将其改为

$x_{n+1}^{\mathrm{S}\mathrm{O}\mathrm{R}}=(1-\omega )x_n^{\mathrm{S}\mathrm{O}\mathrm{R}}+\omega f(x_n^{\mathrm{S}\mathrm{O}\mathrm{R}}).$

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

${𝐱}^{(k+1)}=(1-\omega ){𝐱}^{(k)}+\omega L_{*}^{-1}(𝐛-U{𝐱}^{(k)}),$

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

相關條目

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

参考文献

• Abraham Berman, Robert J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, 1994, SIAM. Template:Isbn.
• Template:MathWorld
• A. Hadjidimos, Successive overrelaxation (SOR) and related methods Template:Wayback, Journal of Computational and Applied Mathematics 123 (2000), 177–199.
• Yousef Saad, Iterative Methods for Sparse Linear Systems Template:Wayback, 1st edition, PWS, 1996.
• Netlib Template:Wayback'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++