古典迭代法可以用于求解线性方程组，特点是格式简单，收敛性与迭代初值无关，与迭代矩阵的谱半径有关。

古典迭代法实质上是基于矩阵的近似逆设计的不动点迭代法，与此不同的，共轭梯度法等则是最优化在矩阵求解问题上的具体体现。

1. 迭代格式与收敛性

求解非奇异的线性方程组

\[ A x = b \]

将方阵\(A\)分解为三个部分

\[ A = D-L-U \]

  • \(D\)为对角阵
  • \(L\)为下三角阵，对角线全零
  • \(U\)为上三角阵，对角线全零

以下的几类迭代法都有如下形式（实质即一种不动点迭代法）

\[ x^{(k+1)} = B x^{(k)} + g \]

最终希望收敛到 \(Ax=b\) 的精确解 \(x\)（在向量范数意义下）

\[ \lim_{k\to \infty} \Vert x^{(k)} - x \Vert = 0 \]

此时\(x\)满足如下等价形式

\[ x = B x + g \]

收敛的充要条件与 \(B\) 的谱半径有关：

\[ \rho(B) < 1 \]

收敛的充分条件有很多，例如 \(B\) 的某个矩阵范数小于 1 即可。收敛与否与迭代初值无关。

2. Jacobi 迭代法

\[ \begin{aligned} A x &= (D-L-U)x = b\\ D x &= (L+U)x + b\\ x &= D^{-1}(L+U) x + D^{-1} b \end{aligned} \]

即取

\[ \begin{aligned} B &= D^{-1}(L+U)\\ g &= D^{-1}b \end{aligned} \]

得到 Jacobi 迭代

\[ x^{(k+1)} = D^{-1}(L+U) x^{(k)} + D^{-1}b \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def JacobiSolver(A, b, ep):
    m, n = A.shape
    n2, = b.shape
    if(m != n or n != n2):
        print("size error")
        return b

    LplusU = -1 * np.array(np.tril(A, -1)) - 1* np.array(np.triu(A, 1))
    D_inv = np.diag(1/np.diag(A))

    x = np.zeros(m)
    iter = 0
    iter_max = 100000
    error = 0
    while iter < iter_max:
        xold = x.copy()
        for i in range(m):
            x[i] = D_inv[i][i]*(b[i] + LplusU[i].dot(xold))

        # print(x)
        iter += 1
        error = np.linalg.norm(x-xold, ord=np.inf)
        if(error < ep):
            return (x,iter,error)

    print("iter >= iter_max")
    return (x,iter,error)

3. G-S 迭代法

作为 Jacobi 迭代的改进，假设每一次更新从上到下进行，那么更新后面的元素时可以直接利用新算出来的值，得到 G-S 迭代法。

此时 \(x\) 每个分量的更新顺序必须是从上往下的，不可以改变。

\[ \begin{aligned} Jacobi: x^{(k+1)} &= D^{-1}L x^{(k)} + D^{-1}U x^{(k)} + D^{-1}b\\ G-S: x^{(k+1)} &= D^{-1}L x^{(k+1)} + D^{-1}U x^{(k)} + D^{-1}b \end{aligned} \]

或者整理成统一形式

\[ x^{(k+1)} = (D-L)^{-1} U x^{(k)} + (D-L)^{-1}b \]

（编程不要采用，因为无法实际算出\((D-L)^{-1}\)，只是用于理论推导）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def GaussSeidelSolver(A, b, ep):
    m,n = A.shape
    n2, = b.shape
    if(m != n or n != n2):
        print("size error")
        return b

    LplusU = -1 * np.array(np.tril(A, -1)) - 1* np.array(np.triu(A, 1))
    D_inv = np.diag(1/np.diag(A))

    x = np.zeros(m)
    iter = 0
    iter_max = 100000
    error = 0
    while iter < iter_max:
        xold = x.copy()
        for i in range(m):
            x[i] = D_inv[i][i]*(b[i] + LplusU[i].dot(x))

        # print(x)
        iter += 1
        error = np.linalg.norm(x-xold, ord=np.inf)
        if(error < ep):
            return (x,iter,error)

    print("iter >= iter_max")
    return (x,iter,error)

4. SOR 迭代法

超松弛迭代法，可以视作 G-S 迭代法的推广或加速。

引入松弛因子\(w > 1\) ，极限情况\(w=1\) 退化为 G-S 迭代法。

\[ \begin{aligned} G-S: x^{(k+1)} &= D^{-1}L x^{(k+1)} + D^{-1}U x^{(k)} + D^{-1}b\\ &= x^{(k)} + \left[D^{-1}L x^{(k+1)} + D^{-1}U x^{(k)} + D^{-1}b - x^{(k)}\right]\\ SOR: x^{(k+1)} &= x^{(k)} + w \left[D^{-1}L x^{(k+1)} + D^{-1}U x^{(k)} + D^{-1}b - x^{(k)}\right]\\ &= (1-w)x^{(k)} + w \left[D^{-1}L x^{(k+1)} + D^{-1}U x^{(k)} + D^{-1}b\right] \end{aligned} \]

推导一下迭代格式

\[ \begin{aligned} (D-wL)x^{(k+1)} &= (1-w)D x^{(k)} + w U x^{(k)} + w b\\ x^{(k+1)} &= (D-wL)^{-1}((1-w)D+wU)x^{(k)} + w(D-wL)^{-1}b\\ &= L_w x^{(k)} + w(D-wL)^{-1} b \end{aligned} \]

其中迭代矩阵

\[ L_w = (D-wL)^{-1}((1-w)D+wU) \]

可以证明 SOR 迭代法收敛的必要条件

\[ 0 < w < 2 \]

因此实际取参数\(w\in(1,2)\)，对于具体问题，存在最优的\(w\)使得收敛速度最快，称为最佳松弛因子。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def SORSolver(A, b,ep,w):
    m, n = A.shape
    n2, = b.shape
    if(m != n or n != n2):
        print("size error")
        return b

    if(w >= 2 or w <= 0):
        print("keep w in (0,2)")
        return b

    LplusU = -1 * np.array(np.tril(A, -1)) - 1* np.array(np.triu(A, 1))
    D_inv = np.diag(1/np.diag(A))

    x = np.zeros(m)
    iter = 0
    iter_max = 100000
    error = 0
    while iter < iter_max:
        xold = x.copy()
        for i in range(m):
            dx = D_inv[i][i]*(b[i] + LplusU[i].dot(x))
            x[i] = (1-w)*x[i]+w*dx

        # print(x)
        iter += 1
        error = np.linalg.norm(x-xold, ord=np.inf)
        if(error < ep):
            return (x,iter,error)

    print("iter >= iter_max")
    return (x,iter,error)

5. 近似逆

古典迭代法实质是取了\(A\)的一个可计算的近似逆\(M\)，这里称\(M\)为\(A\)的近似(左)逆，含义为

\[ \Vert I-MA \Vert < \Delta \le 1 \]

即\(MA\)在矩阵范数意义下(或谱半径意义下)近似为单位阵\(I\)。那么

\[ \begin{aligned} A x &= b \\ (I-MA)x &= x - M b \\ x &= (I-MA) x + M b \\ \end{aligned} \]

由此诱导出迭代格式

\[ x^{(k+1)} = (I-MA) x^{(k)} + M b \]

这里的\(I-MA\)即前文中的迭代矩阵\(B\)，因此收敛性充要条件就是\(I-MA\)的谱半径小于 1，含义为： 当可计算的\(M\)可以作为\(A\)的一个近似逆时，对应的迭代格式收敛；并且\(M\)越接近\(A^{-1}\)，收敛速度越快。

对应于前文中的迭代法，可以发现

  • Jacobi 迭代法：取对角阵的逆\(D^{-1}\)来近似\(A^{-1}\)
  • G-S 迭代法：取对角线和下三角部分之和的逆\((D-L)^{-1}\)来近似\(A^{-1}\)
  • SOR 迭代法: 取\(w(D-wL)^{-1}\)来近似\(A^{-1}\)