茶年 · 中文维基百科

逐次超松弛迭代法

Wikipedia · 2026-01-28 维基百科 👁 1 阅读

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

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

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

:A\mathbf x = \mathbf b

其中

: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：

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

线性方程组可重写为

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

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

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

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

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

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

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

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

  • 松弛因子适当： \omega \in (0,2)
  • 雅可比法迭代矩阵 C_\text{Jac}:= I-D^{-1}A 只有实特征值
  • 雅可比法收敛： \mu := \rho(C_\text{Jac})
  • 矩阵分解 A=D+L+U 满足\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)

则收敛速度可表为
:
\rho(C_\omega) =
\begin{cases}
\frac{1}{4} \left( \omega \mu + \sqrt{\omega^2 \mu^2-4(\omega-1)} \right)^2\,,
& 0
最优松弛因子是
:
\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 + \tfrac{1}{z}U对对角阵Z，其元素Z_{ii}=z^{i-1}、 \operatorname{det}(\lambda D + L + U) = \operatorname{det}(Z(\lambda D + L + U)Z^{-1}) 。

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

输入：, ,
输出：

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

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

例子
解线性方程组

:
\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}

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

用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) = 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=\left(\frac{D}{\omega}+L\right)\frac{\omega}{2-\omega}D^{-1}\left(\frac{D}{\omega}+U\right),

迭代法为

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

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

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

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

则可将其改为

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

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

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

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

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

注释
参考文献

  • Abraham Berman, Robert J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, 1994, SIAM. .

*

外部链接

#带有伪代码示例的条目 #带有Python代码示例的条目 #数值线性代数

评论 (0)

  • 还没有评论，来抢沙发吧。