逐次超松弛迭代法

数值线性代数中，逐次超松弛（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.}$

松弛因子${\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 }$。这些方法通常只对一部分问题有效。

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

• 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++