查看“︁逐次超松弛迭代法”︁的源代码

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

SOR迭代法由David M. Young Jr.和Stanley P. Frankel在1950年同时独立提出，目的是在计算机上自动求解线性方程组。之前，人们已经为[[计算员]]的计算开发过超松弛法，如[[路易斯·弗莱·理查德森]]的方法以及R. V. Southwell开发的方法。但这些方法需要一定专业知识确保求解的收敛，不适用于计算机编程。David M. Young Jr.的论文对这些方面进行了探讨。<ref>{{citation
 |last=Young
 |first=David M.
 |title=Iterative methods for solving partial difference equations of elliptical type
 |url=http://www.ma.utexas.edu/CNA/DMY/david_young_thesis.pdf
 |date=1950-05-01
 |series=PhD thesis, Harvard University
 |access-date=2009-06-15
 |archive-date=2011-06-07
 |archive-url=https://web.archive.org/web/20110607135504/http://www.ma.utexas.edu/CNA/DMY/david_young_thesis.pdf
 |dead-url=no
 }}</ref>

==形式化==
给定''n''个线性方程组成的方系统：

:<math>A\mathbf x = \mathbf b</math>

其中

:<math>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}.</math>

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

:<math>A=D+L+U, </math>
其中
:<math>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}. </math>

线性方程组可重写为

:<math>(D+\omega L) \mathbf{x} = \omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x} </math>

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

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

:<math> \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}, </math>

其中<math>\mathbf{x}^{(k)}</math>是<math>\mathbf{x}</math>的第''k''次迭代值，<math>\mathbf{x}^{(k+1)}</math>是<math>\mathbf{x}</math>下一次迭代所得的值。
利用<math>(D+\omega L)</math>的三角形，可用[[三角矩阵#向前与向后替换|向前替换]]法依次计算<math>\mathbf{x}^{(k+1)}</math>的元素：

:<math> x^{(k+1)}_i  = (1-\omega)x^{(k)}_i + \frac{\omega}{a_{ii}} \left(b_i - \sum_{j<i} a_{ij}x^{(k+1)}_j - \sum_{j>i} a_{ij}x^{(k)}_j \right),\quad i=1,2,\ldots,n. </math>

== 收敛性 ==
[[File:Spectral Radius.svg|thumb|SOR迭代法迭代矩阵<math> C_\omega </math>的谱半径<math> \rho(C_\omega) </math>。
图表显示的是雅可比迭代矩阵的谱半径<math> \mu := \rho(C_\text{Jac}) </math>。]]
松弛因子<math>\omega</math>的选择并不容易，取决于系数矩阵的性质。1947年，[[亚历山大·马雅科维奇·奥斯特洛夫斯基]]证明，若''A''是[[对称矩阵|对称]][[正定矩阵]]，则<math>\forall 0<\omega<2,\ \rho(L_\omega)<1</math>。因此，迭代过程将收敛，但更高的收敛速度更有价值。

=== 收敛速度 ===
SOR法的收敛速度可通过分析得出。假设<ref>{{Cite book|title=Iterative Solution of Large Sparse Systems of Equations {{!}} SpringerLink|volume = 95|last=Hackbusch|first=Wolfgang|year=2016|isbn=978-3-319-28481-1|location=|pages=|language=en-gb|chapter=4.6.2|doi=10.1007/978-3-319-28483-5|series = Applied Mathematical Sciences}}</ref><ref>{{Cite book|title=Iterative Methods for Solving Linear Systems|series = Frontiers in Applied Mathematics|volume = 17|last=Greenbaum|first=Anne|year=1997|isbn=978-0-89871-396-1|language=en-gb|chapter=10.1|doi=10.1137/1.9781611970937}}</ref>
* 松弛因子适当：<math> \omega \in (0,2) </math>
* [[雅可比法]]迭代矩阵<math> C_\text{Jac}:= I-D^{-1}A </math>只有实特征值
* [[雅可比法]]收敛：<math> \mu := \rho(C_\text{Jac}) < 1 </math>
* 矩阵分解<math> A=D+L+U </math>满足<math>\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) </math>
则收敛速度可表为
:<math>
\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}
</math> 
最优松弛因子是
:<math>
\omega_\text{opt} := 1+ \left( \frac{\mu}{1+\sqrt{1-\mu^2}} \right)^2 = 1 + \frac{\mu^2}{4} + O(\mu^3)\,.
</math>
特别地，<math>\omega = 1</math>时SOR法即退化为[[高斯-赛德尔迭代]]，有<math>\rho(C_\omega)=\mu^2=\rho(C_\text{Jac})^2</math>。
对最优的<math>\omega</math>，有<math>\rho(C_\omega)=\frac{1-\sqrt{1-\mu^2}}{1+\sqrt{1-\mu^2}} = \frac{\mu^2}{4} + O(\mu^3)</math>，表明SOR法的效率约是高斯-赛德尔迭代的4倍。

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

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

 输入：{{mvar|A}}, {{mvar|b}}, {{mvar|ω}}
 输出：{{mvar|φ}}
 <!-- the leading whitespace in the next line is significant -->
 
 选择初始解{{mvar|φ}}
 '''repeat''' until convergence
     '''for''' {{mvar|i}} '''from''' 1 '''until''' {{mvar|n}} '''do'''
         set {{mvar|σ}} to 0
         '''for''' {{mvar|j}} '''from''' 1 '''until''' {{mvar|n}} '''do'''
             '''if''' {{mvar|j}} &ne; {{mvar|i}} '''then'''
                 set {{mvar|σ}} to {{math|''σ'' + ''a<sub>ij</sub>'' ''φ<sub>j</sub>''}}
             '''end if'''
         '''end''' ({{mvar|j}}-loop)
         set {{math|''φ<sub>i</sub>''}} to {{math|(1 − ''ω'')''φ<sub>i</sub>'' + ''ω''(''b<sub>i</sub>'' − ''σ'') / ''a<sub>ii</sub>''}}
     '''end''' ({{mvar|i}}-loop)
     check if convergence is reached
 '''end''' (repeat)

;注意：<math>(1-\omega)\phi_i + \frac{\omega}{a_{ii}} (b_i - \sigma)</math>也可写作<math>\phi_i + \omega \left( \frac{b_i - \sigma}{a_{ii}} - \phi_i\right)</math>，这样每次外层''for''循环可以省去一次乘法。

==例子==
解线性方程组

:<math>
  \begin{align}
     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{align}
</math>

择松弛因子<math>\omega = 0.5</math>与初始解<math>\phi = (0, 0, 0, 0)</math>。由SOR算法可得下表，在38步取得精确解{{nowrap|(3, &minus;2, 2, 1)}}。

{| class="wikitable" border="1"
|-
! 迭代
! <math>x_1</math>
! <math>x_2</math>
! <math>x_3</math>
! <math>x_4</math>
|-
| {{0}}1
| 0.25
| &minus;2.78125
| 1.6289062
| 0.5152344
|-
| {{0}}2
| 1.2490234
| &minus;2.2448974
| 1.9687712
| 0.9108547
|-
| {{0}}3
| 2.070478
| &minus;1.6696789
| 1.5904881
| 0.76172125
|-
| ...
| ...
| ...
| ...
| ...
|-
| 37
| 2.9999998
| &minus;2.0
| 2.0
| 1.0
|-
| 38
| 3.0
| &minus;2.0
| 2.0
| 1.0
|}

用Common Lisp的简单实现：

<syntaxhighlight lang="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+))))))
</syntaxhighlight>

上述伪代码的简单Python实现。
<syntaxhighlight lang="python3">
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)
</syntaxhighlight>

==对称逐次超松弛==
对对称矩阵''A''，其中

:<math> U=L^T,\,</math>

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

:<math>P=\left(\frac{D}{\omega}+L\right)\frac{\omega}{2-\omega}D^{-1}\left(\frac{D}{\omega}+U\right),</math>

迭代法为

:<math>\mathbf{x}^{k+1}=\mathbf{x}^k-\gamma^k P^{-1}(A\mathbf{x}^k-\mathbf{b}),\ k \ge 0.</math>

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

==其他应用==
{{main|理查德森外推法}}
任何迭代法都可应用相似技巧。若原迭代的形式为

:<math>x_{n+1}=f(x_n)</math>

则可将其改为

:<math>x^\mathrm{SOR}_{n+1}=(1-\omega)x^{\mathrm{SOR}}_n+\omega f(x^\mathrm{SOR}_n).</math>

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

:<math> \mathbf{x}^{(k+1)} = (1-\omega)\mathbf{x}^{(k)} + \omega L_*^{-1} (\mathbf{b} - U\mathbf{x}^{(k)}),</math>

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

==相關條目==
*[[雅可比法]]
*[[置信传播]]
*[[矩阵分裂]]

==注释==
{{Reflist}}

==参考文献==
* Abraham Berman, Robert J. Plemmons, ''Nonnegative Matrices in the Mathematical Sciences'', 1994, SIAM. {{isbn|0-89871-321-8}}.
*{{MathWorld|urlname=SuccessiveOverrelaxationMethod|title=Successive Overrelaxation Method|author=Black, Noel|author2=Moore, Shirley|name-list-style=amp}}
* A. Hadjidimos, ''[http://www.sciencedirect.com/science/article/pii/S0377042700004039 Successive overrelaxation (SOR) and related methods] {{Wayback|url=http://www.sciencedirect.com/science/article/pii/S0377042700004039 |date=20240426063654 }}'', Journal of Computational and Applied Mathematics 123 (2000), 177–199.
* Yousef Saad, ''[http://www-users.cs.umn.edu/%7Esaad/books.html Iterative Methods for Sparse Linear Systems] {{Wayback|url=http://www-users.cs.umn.edu/%7Esaad/books.html |date=20141223032757 }}'', 1st edition, PWS, 1996.
* [http://www.netlib.org/utk/papers/templates/node11.html Netlib] {{Wayback|url=http://www.netlib.org/utk/papers/templates/node11.html |date=20160304112841 }}'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)

==外部链接==
* [https://web.archive.org/web/20040317151200/http://math.fullerton.edu/mathews/n2003/SORmethodMod.html Module for the SOR Method]
* [https://web.archive.org/web/20140310122757/http://michal.is/projects/tridiagonal-system-solver-sor-c/ Tridiagonal linear system solver] based on SOR, in C++

[[Category:数值线性代数]]
[[Category:带有伪代码示例的条目]]
[[Category:带有Python代码示例的条目]]