逐次超松弛迭代法

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

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

形式化

编辑

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

 

其中

 

A可分解为对角矩阵D、严格上下三角矩阵UL

 

其中

 

线性方程组可重写为

 

其中 是常数,称作松弛因子(relaxation factor)。

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

 

其中  的第k次迭代值,  下一次迭代所得的值。 利用 的三角形,可用向前替换法依次计算 的元素:

 

收敛性

编辑
 
SOR迭代法迭代矩阵 的谱半径 。 图表显示的是雅可比迭代矩阵的谱半径 

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

收敛速度

编辑

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

  • 松弛因子适当: 
  • 雅可比法迭代矩阵 只有实特征值
  • 雅可比法收敛: 
  • 矩阵分解 满足 

则收敛速度可表为

 

最优松弛因子是

 

特别地, 时SOR法即退化为高斯-赛德尔迭代,有 。 对最优的 ,有 ,表明SOR法的效率约是高斯-赛德尔迭代的4倍。

最后一条假设对三对角矩阵也满足,因为 对对角阵Z,其元素  

算法

编辑

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

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

选择初始解φ
repeat until convergence
    for i from 1 until n do
        set σ to 0
        for j from 1 until n do
            if ji 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)
注意: 也可写作 ,这样每次外层for循环可以省去一次乘法。

例子

编辑

解线性方程组

 

择松弛因子 与初始解 。由SOR算法可得下表,在38步取得精确解(3, −2, 2, 1)

迭代        
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,其中

 

对称逐次超松弛迭代法SSOR):

 

迭代法为

 

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

其他应用

编辑

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

 

则可将其改为

 

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

 

其中  用于加快收敛速度, 可使发散的迭代收敛或加快过调(overshoot)过程的收敛。有多种方法可根据观察到的收敛过程行为,自适应地调整松弛因子 。这些方法通常只对一部分问题有效。

另见

编辑

注释

编辑
  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)

外部链接

编辑