数值分析——解线性方程组方法补充

数值分析——解线性方程组补充

• 前言
• LU直接分解法 — Doolittle_Decomposition
• • 代码实现
  • test example
• 三对角矩阵的特殊解法
• • 代码实现
  • test example
• 迭代法解线性方程组
• • Jacobi迭代法
  • • 代码实现
    • test example
  • Gauss-Seidel迭代法
  • • 代码实现
    • test example
  • 超松弛迭代法SOR
  • 未完待续

前言

​  首先很多内容都跟矩阵论的部分重合了，相关的关于LU分解、LDV分解、LU解线性方程组、求行列式、求逆等都可以在矩阵论专栏中对应查看，这里仅补充一些其他方法。

LU直接分解法 — Doolittle_Decomposition

• 道立特分解属于一种先设定好U矩阵第一行和L矩阵第一列的值，便可以开始迭代求解L和U的其他值的算法，分解结果如下（以三阶矩阵为例）：
  [ a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 ] = [ 1 0 0 l 21 1 0 l 31 l 32 1 ] [ u 11 u 12 u 13 0 u 22 u 23 0 0 u 33 ] begin{bmatrix} a_{11} & a_{12} & a_{13} \ a_{21} & a_{22} & a_{23} \ a_{31} & a_{32} & a_{33} end{bmatrix}= begin{bmatrix} 1 & 0 & 0 \ l_{21} & 1& 0 \ l_{31} & l_{32} & 1 end{bmatrix} begin{bmatrix} u_{11} & u_{12} & u_{13} \ 0 & u_{22} & u_{23} \ 0 & 0 & u_{33} end{bmatrix} ⎣⎡​a11​a21​a31​​a12​a22​a32​​a13​a23​a33​​⎦⎤​=⎣⎡​1l21​l31​​01l32​​001​⎦⎤​⎣⎡​u11​00​u12​u22​0​u13​u23​u33​​⎦⎤​
• 算法流程（先求解U的对应行再求解L的对应列）：
  1. 求解U的第i行第j个元素： u i j = a i j − ∑ k = 1 i − 1 l i k ∗ u k j u_{ij}=a_{ij}-sum_{k=1}^{i-1}l_{ik}*u_{kj} uij​=aij​−∑k=1i−1​lik​∗ukj​
  2. 求解L的第j列第i个元素： l j i = ( a j i − ∑ k = 1 i − 1 l j k ∗ u k i ) / u i i l_{ji}=(a_{ji}-sum_{k=1}^{i-1}l_{jk}*u_{ki})/u_{ii} lji​=(aji​−∑k=1i−1​ljk​∗uki​)/uii​

代码实现

@staticmethod
def domain_row_transform(arr, eps=1e-6, Test=False):
    """domain_row_transform 选取主元并进行行变换

    Args:
        arr ([np.darray]): [input matrix]
        eps ([float]): numerical threshold.
        test (bool, optional): [show checking information]. Defaults to False.

    Returns:
        [A]: [resort one of arr]
    
    Author: Junno
	Date: 2022-04-27
    """
    assert len(arr.shape) == 2
    A = np.copy(arr)
    M, N = A.shape
    for i in range(M):
        for j in range(i, N):
            if Test:
                print('During process in row:{}, col:{}'.format(i, j))
            if sum(abs(A[i:, j])) > eps:
                zero_index = []
                non_zero_index = []
                for k in range(i, M):
                    if abs(A[k, j]) > eps:
                        non_zero_index.append(k)
                    else:
                        zero_index.append(k)

                non_zero_index = sorted(
                    non_zero_index, key=lambda x: abs(A[x, j]))
                sort_idnex = non_zero_index+zero_index

                if Test:
                    print('Sorting index for {} cols:n'.format(i), sort_idnex)

                if sort_idnex[0] != i:
                    A[[sort_idnex[0], i], :] = A[[i, sort_idnex[0]], :]
                if Test:
                    print('After resortingn', A)
    return A

@classmethod
def Doolittle_Decomposition(self, target, domain_selection=True, test=False):
    """Doolittle_Decomposition equal to LU_Decomposition

    Args:
        target ([np.darray]): [input matrix]
        domain_selection (bool, optional): [operate domain selection on target]. Defaults to True.
        test (bool, optional): [show testing information]. Defaults to False.

    Returns:
        [resort input, L, U] if domain_selection=True else [L,U]

Author: Junno
	Date: 2022-04-27
    """
    assert len(target.shape) == 2
    arr = np.copy(target)
    M, N = arr.shape
    L = np.zeros_like(arr, dtype=arr.dtype)
    U = np.zeros_like(arr, dtype=arr.dtype)
    if domain_selection:
        arr = self.domain_row_transform(arr)
    for i in range(M):
        if test:
            print(i)
            print(L, U, sep='n')

        for j in range(i, N):
            if test:
                print('U', i, j, arr[i, j], L[i, :i], U[:i, j])
            U[i, j] = arr[i, j]-L[i, :i]@U[:i, j]
            if i == j:
                L[i, j] = 1
            else:
                if test:
                    print('L', j, i, arr[j, i],
                          L[j, :i], U[:i, i], U[i, i])
                L[j, i] = (arr[j, i]-(L[j, :i]@U[:i, i]))/U[i, i]
    if domain_selection:
        print('return resort input matrix, L, U')
        return arr, L, U
    else:
        return L, U

test example

• Tips：为避免大数吃小数的问题和舍入误差问题，加入了选主元的选项，可以对原输入矩阵进行行变换，不改变结果的正确性，可以自行选择

A=np.random.randn(5,5)
A
>>>
array([[-1.50948856, -1.61033161, -0.20671825,  0.4286645 , -0.2572608 ],
       [-1.28615783,  0.87991375,  1.71080037, -1.61832726,  1.52242031],
       [-0.98092031, -0.32273876,  2.20229521,  0.853808  , -2.15168542],
       [-2.02003405,  0.05052931, -1.01116616,  0.26437957, -1.65870163],
       [-0.83327297,  0.32495716, -1.02878403, -0.96607898, -0.59284746]])
       
L,U=Doolittle_Decomposition(A,domain_selection=False,test=False)
L
>>>
[[    1.        0.        0.        0.        0.   ]
 [    4.688     1.        0.        0.        0.   ]
 [    5.484 -3863.106     1.        0.        0.   ]
 [    8.449 -8279.836     2.143     1.        0.   ]
 [   -8.763 10132.862    -2.623   -11.842     1.   ]]

U
>>>
[[   -0.165     0.618    -0.885    -0.515    -1.011]
 [    0.        0.001     3.434     2.532     6.176]
 [    0.        0.    13269.089  9782.795 23865.142]
 [    0.        0.        0.        0.216    -4.26 ]
 [    0.        0.        0.        0.      -43.798]]

print(L@U)
>>>
array([[-1.50948856, -1.61033161, -0.20671825,  0.4286645 , -0.2572608 ],
       [-1.28615783,  0.87991375,  1.71080037, -1.61832726,  1.52242031],
       [-0.98092031, -0.32273876,  2.20229521,  0.853808  , -2.15168542],
       [-2.02003405,  0.05052931, -1.01116616,  0.26437957, -1.65870163],
       [-0.83327297,  0.32495716, -1.02878403, -0.96607898, -0.59284746]])

三对角矩阵的特殊解法

• 对于一般的线性方程组，做高斯消去法的乘法复杂度大概为 O ( N 3 ) O(N^3) O(N3)， 当N超过100时就比较吃力了。在很多物理仿真的求解计算中，最常见的问题便是求解线性方程组，但是它们的维度一般都成千上万。但同时也存在一个比较有用的特性，即稀疏特性，现在有很多先进的算法来求解这些问题。这里我们介绍一种特殊情况下，即对角占优下三对角线性方程组的快速解法——Thomas分解或者追赶法，它的复杂度可以与N同阶，为 O ( 3 N − 2 ) O(3N-2) O(3N−2)。

• 对角占优的三对角方程组形如下：

[ b 1 c 1 0 ⋯ 0 a 2 b 2 c 2 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 a n − 1 b n − 1 c n − 1 0 0 0 a n b n ] [ x 1 x 2 ⋮ x n ] = [ f 1 f 2 ⋮ f n ] begin{bmatrix} b_{1} & c_{1} & 0 & cdots & 0 \ a_{2} & b_{2} & c_{2} & cdots & 0\ 0 & ddots & ddots &ddots&0\ 0 & 0 & a_{n-1} & b_{n-1} & c_{n-1}\ 0& 0& 0& a_{n} & b_{n} end{bmatrix} begin{bmatrix} x_1\ x_2\ vdots\ x_{n} end{bmatrix} = begin{bmatrix} f_1\ f_2\ vdots\ f_{n} end{bmatrix} ⎣⎢⎢⎢⎢⎡​b1​a2​000​c1​b2​⋱00​0c2​⋱an−1​0​⋯⋯⋱bn−1​an​​000cn−1​bn​​⎦⎥⎥⎥⎥⎤​⎣⎢⎢⎢⎡​x1​x2​⋮xn​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​f1​f2​⋮fn​​⎦⎥⎥⎥⎤​

• 对角占优条件：

  1. ∣ b 1 ∣ > ∣ c 1 ∣ > 0 |b_1| > |c_1| >0 ∣b1​∣>∣c1​∣>0
  2. ∣ b i ∣ ≥ ∣ a i ∣ + ∣ c i ∣ , a i , c i ≠ 0 , ( i = 2 , 3 , . . . , n − 1 ) |b_i| geq |a_i| +|c_i|, a_i,c_i neq0,(i=2,3,...,n-1) ∣bi​∣≥∣ai​∣+∣ci​∣,ai​,ci​​=0,(i=2,3,...,n−1)
  3. ∣ b n ∣ > ∣ a n ∣ > 0 |b_n| > |a_n|>0 ∣bn​∣>∣an​∣>0

• 类似于LU分解解线性方程组的做法，尝试将三对角矩阵分解为上下三角矩阵，形如：
  A = [ b 1 c 1 0 ⋯ 0 a 2 b 2 c 2 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 a n − 1 b n − 1 c n − 1 0 0 0 a n b n ] = L U = [ α 1 0 0 ⋯ 0 γ 2 α 2 0 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 γ n − 1 α n − 1 0 0 0 0 γ n α n ] [ 1 β 1 0 ⋯ 0 0 1 β 2 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 0 1 β n − 1 0 0 0 0 1 ] A=begin{bmatrix} b_{1} & c_{1} & 0 & cdots & 0 \ a_{2} & b_{2} & c_{2} & cdots & 0\ 0 & ddots & ddots &ddots&0\ 0 & 0 & a_{n-1} & b_{n-1} & c_{n-1}\ 0& 0& 0& a_{n} & b_{n} end{bmatrix}\ =LU=begin{bmatrix} alpha_{1} & 0 & 0 & cdots & 0 \ gamma_{2} & alpha_{2} & 0 & cdots & 0\ 0 & ddots & ddots &ddots&0\ 0 & 0 & gamma_{n-1} & alpha_{n-1} & 0\ 0& 0& 0& gamma_{n} & alpha_{n} end{bmatrix} begin{bmatrix} 1 & beta_{1} & 0 & cdots & 0 \ 0 & 1 & beta_{2} & cdots & 0\ 0 & ddots & ddots &ddots&0\ 0 & 0 & 0 &1 & beta_{n-1}\ 0& 0& 0& 0 & 1 end{bmatrix} A=⎣⎢⎢⎢⎢⎡​b1​a2​000​c1​b2​⋱00​0c2​⋱an−1​0​⋯⋯⋱bn−1​an​​000cn−1​bn​​⎦⎥⎥⎥⎥⎤​=LU=⎣⎢⎢⎢⎢⎡​α1​γ2​000​0α2​⋱00​00⋱γn−1​0​⋯⋯⋱αn−1​γn​​0000αn​​⎦⎥⎥⎥⎥⎤​⎣⎢⎢⎢⎢⎡​10000​β1​1⋱00​0β2​⋱00​⋯⋯⋱10​000βn−1​1​⎦⎥⎥⎥⎥⎤​

• 比较待定系数，可以逐步求解出各系数，（推导过程详见教材第七章P184）此时 A x = f Ax=f Ax=f等价于两个三角方程组 L y = f Ly=f Ly=f和 U x = y Ux=y Ux=y，得到追赶法求解过程：

1. 计算 β i beta_i βi​： β 1 = c 1 / b 1 , β i = c i / ( b i − a i β i − 1 ) , i = ( 2 , 3 , . . . , n − 1 ) beta_1=c_1/b_1,beta_i=c_i/(b_i-a_i beta_{i-1}), i=(2,3,...,n-1) β1​=c1​/b1​,βi​=ci​/(bi​−ai​βi−1​),i=(2,3,...,n−1)
2. 解 L y = f Ly=f Ly=f： y 1 = f 1 / b 1 , y i = ( f i − a i y i − 1 ) / ( b i − a i β i − 1 ) , i = ( 2 , 3 , . . . , n ) y_1=f_1/b_1,y_i=(f_i-a_i y_{i-1})/(b_i-a_i beta_{i-1}), i=(2,3,...,n) y1​=f1​/b1​,yi​=(fi​−ai​yi−1​)/(bi​−ai​βi−1​),i=(2,3,...,n)
3. 解 U x = y Ux=y Ux=y： x n = y n , x i = y i − β i x i + 1 , ( i = n − 1 , n − 2 , . . . , 2 , 1 ) x_n=y_n,x_i=y_i-beta_i x_{i+1}, (i=n-1,n-2,...,2,1) xn​=yn​,xi​=yi​−βi​xi+1​,(i=n−1,n−2,...,2,1)

代码实现

• 实际计算时只需要传入两个次对角线和主对角线的三个向量值便可以求解，内存消耗更小！

def check_SDD_vec(a, b, c, eps=1e-4, test=False):
        """check_SDD_vec 向量形式的严格对角占优条件判断

        Args:
            a ([np.darray]): [lower-secondary diagonal]
            b ([np.darray]): [diagonal]
            c ([np.darray]): [upper-secondary diagonal]
            eps ([type], optional): [threshold]. Defaults to 1e-4.
            test (bool, optional): [show testing information]. Defaults to False.

        Returns:
            [bool]: [whether or not]

		Author: Junno
    	Date: 2022-04-28
        """
        assert len(a) == len(c) and len(a)+1 == len(b)
        N = len(b)
        a_ = np.concatenate(([eps], a), axis=0)
        c_ = np.concatenate((c, [eps]), axis=0)
        if test:
            print(a_, c_)
        sum_a_c = np.abs(a_)+np.abs(c_)
        b_ = np.abs(b)
        for i in range(N):
            if b_[i] < sum_a_c[i]:
                return False
        return True

def check_SDD_M(arr,eps=1e-4):
    """check_SDDM 判断输入矩阵是否满足严格对角占优条件
    Args:
        arr ([np.darray]): [input matrix]
        eps ([type], optional): [threshold]. Defaults to 1e-4.

    Returns:
        [bool]: [whether or not]

	Author: Junno
    Date: 2022-04-28
    """
    assert arr.shape[0]==arr.shape[1]
    return check_SDD_vec(np.diagonal(arr,offset=-1),np.diagonal(arr,offset=0),np.diagonal(arr,offset=1),eps=eps)

def Thomas_Linear_Equation_Solve(a, b, c, d, eps=1e-4, test=False):
    """Thomas_Linear_Equation_Solve 

    Args:
        a ([list or np.darray]): [lower-secondary diagonal]
        b ([list or np.darray]): [diagonal]
        c ([list or np.darray]): [upper-secondary diagonal]
        d ([list or np.darray]): [b of Ax=b]
        eps ([type], optional): [threshold]. Defaults to 1e-4.
        test (bool, optional): [show testing information]. Defaults to False.

    Returns:
        x [list]: [solution of Ax=b]

    Author: Junno

    Date: 2022-04-28
    """
    assert len(a) == len(c) and len(a)+1 == len(b)
    if not check_SDD_vec(a, b, c, eps=eps):
        print('The answer may be wrong because the input are not meet SDD condition completely')
    N = len(b)
    if test:
        print(a, b, c, d)
    beta = np.zeros_like(b)
    x = np.zeros_like(b)
    y = np.zeros_like(b)
    beta[0] = c[0]/b[0]
    y[0] = d[0]/b[0]
    for i in range(1, N):
        if i < N-1:
            beta[i] = c[i]/(b[i]-a[i-1]*beta[i-1])
        y[i] = (d[i]-a[i-1]*y[i-1])/(b[i]-a[i-1]*beta[i-1])
    for i in range(N-1, -1, -1):
        if i == N-1:
            x[i] = y[i]
        else:
            x[i] = y[i]-beta[i]*x[i+1]

    return x

test example

A=np.array([[3,1,0,0],[1,2,1,0],[0,1,3,2],[0,0,3,4]],dtype=np.float32).reshape((4,4))
A
>>>
array([[3., 1., 0., 0.],
       [1., 2., 1., 0.],
       [0., 1., 3., 2.],
       [0., 0., 3., 4.]], dtype=float32)


# create testing example
x=np.random.randn(4,1)
print(x)
>>>
[[-0.35812746]
 [ 0.44621921]
 [-0.45669147]
 [ 1.85002579]]
 
f=A@x
print(f)
>>>
[[-0.62816316]
 [ 0.07761949]
 [ 2.77619638]
 [ 6.03002876]]

Thomas_Linear_Equation_Solve(np.diagonal(A,offset=-1),np.diagonal(A,offset=0),np.diagonal(A,offset=1),f,test=False)
>>>
array([-0.35812744,  0.44621915, -0.45669138,  1.8500258 ], dtype=float32)

迭代法解线性方程组

• 对于大型的稀疏矩阵，有时候迭代求解方法在效率和成本上较直接求解更优。对于方程组 A x = b Ax=b Ax=b，可以通过变形得到 x = B x + f x=B x+f x=Bx+f，再通过迭代方法，逐步逼近真实解， x k + 1 = B x k + f x^{k+1}=B x^k +f xk+1=Bxk+f。其中，B称为迭代传递矩阵。

• 如果 lim ⁡ k → ∞ x ( k ) lim_{k rightarrow infty}x^{(k)} limk→∞​x(k)存在，称迭代过程收敛。引进误差向量， ϵ k + 1 = x k + 1 − x ∗ epsilon^{k+1}=x^{k+1}-x^* ϵk+1=xk+1−x∗，递推可以得到， ϵ k = B x k − 1 = B k ϵ 0 epsilon^{k}=B x^{k-1}=B^k epsilon^0 ϵk=Bxk−1=Bkϵ0。即若想要迭代收敛，需要知道B满足什么情况时有 B k → 0    ( k → ∞ ) B^k rightarrow 0 ,, (k rightarrow infty) Bk→0(k→∞)。

• 迭代法基本定理：当 ρ ( B ) < 1 rho(B) <1 ρ(B)<1时，上述迭代法在任意初始向量 x 0 x_0 x0​及任意f开始迭代后均收敛。其中， ρ ( B ) rho(B) ρ(B)表示B的谱半径。可以由Jordan分解简单证明得到，可以参考教材P206-207。

• 由Jordan的证明过程还可以得到，当 ρ ( B ) rho(B) ρ(B)越小时，迭代收敛的速度越快，若给定求解精度为 1 0 − s 10^{-s} 10−s，则至少需要的迭代次数可以由以下式子确定：
  { [ ρ ( B ) ] k ≤ 1 0 − s k ≥ s ln ⁡ 10 − ln ⁡ ρ ( B ) begin{cases} [rho(B)]^k leq 10^{-s}\ k geq frac{s ln{10}}{-ln{rho(B)}}\ end{cases} {[ρ(B)]k≤10−sk≥−lnρ(B)sln10​​

• 由于任意阶矩阵范数都是谱半径的上届，所以可以有一些容易计算的矩阵范数得到收敛性的大概估计。设某一矩阵范数表示为 ∣ ∣ B ∣ ∣ v = q < 1 ||B||_v=q <1 ∣∣B∣∣v​=q<1，第k步的误差估计可以表示为：
  ∣ ∣ x ∗ − x k ∣ ∣ v ≤ q 1 − q ∣ ∣ x k − x k − 1 ∣ ∣ v ||x^*-x^k||_v leq frac{q}{1-q} ||x^k-x^{k-1}||_v ∣∣x∗−xk∣∣v​≤1−qq​∣∣xk−xk−1∣∣v​

Jacobi迭代法

• 设有方程组 ∑ j = 1 n a i j x j = b i sum_{j=1}^n a_{ij}x_j=b_i ∑j=1n​aij​xj​=bi​，记为 A x = b Ax=b Ax=b，将A分解为 A = D + L + U A=D+L+U A=D+L+U，格式如下：
  D = [ a 11 a 22 ⋱ a n n ] D=begin{bmatrix} a_{11} & & & \ & a_{22} & & \ & &ddots& \ & & &a_{nn} \ end{bmatrix} D=⎣⎢⎢⎡​a11​​a22​​⋱​ann​​⎦⎥⎥⎤​

L = [ 0 a 21 0 ⋮ ⋱ ⋱ a n 1 ⋯ a n , n − 1 0 ] L=begin{bmatrix} 0 & & & \ a_{21} & 0 & & \ vdots & ddots &ddots& \ a_{n1} & cdots & a_{n,n-1} &0 \ end{bmatrix} L=⎣⎢⎢⎢⎡​0a21​⋮an1​​0⋱⋯​⋱an,n−1​​0​⎦⎥⎥⎥⎤​

U = [ 0 a 12 ⋯ a 1 n 0 ⋱ ⋮ ⋱ a n − 1 , n 0 ] U=begin{bmatrix} 0 & a_{12} & cdots & a_{1n} \ & 0& ddots & vdots \ & &ddots&a_{n-1,n} \ & & &0 \ end{bmatrix} U=⎣⎢⎢⎢⎡​0​a12​0​⋯⋱⋱​a1n​⋮an−1,n​0​⎦⎥⎥⎥⎤​

• 对于方程组 ∑ j = 1 n a i j x j = b i sum_{j=1}^n a_{ij}x_j=b_i ∑j=1n​aij​xj​=bi​，可以改写为 x i = 1 a i i ( b i − ∑ j = 1 , j ≠ i n a i j x j ) , ( i = 1 , 2 , . . . , n ) x_i=frac{1}{a_{ii}}(b_i-sum_{j=1,j neq i}^n a_{ij}x_j),(i=1,2,...,n) xi​=aii​1​(bi​−j=1,j​=i∑n​aij​xj​),(i=1,2,...,n)若是以 x = B x + f x=Bx+f x=Bx+f，则可以对照写出矩阵形式：
  { B = − D − 1 ( L + U ) f = D − 1 b begin{cases} B=-D^{-1}(L+U)\ f=D^{-1}b\ end{cases} {B=−D−1(L+U)f=D−1b​

代码实现

def Jacobi_iteration(A, b, N=0, eps=1e-4, x_0=None, norm='infty', own_method=False, test=False):
    """Jacobi_iteration 

    Args:
        A ([np.darray]): [A]
        b ([np.darray]): [b]
        N (int, optional): [iteration time, it will give a suitable one when N is 0]. Defaults to 0.
        eps ([float], optional): [error threshold]. Defaults to 1e-4.
        x_0 ([np.darray], optional): [initial x for iteration]. Defaults to None.
        norm (str, optional): [norm method to calculate fitting error,'infty','L1']. Defaults to 'infty'.
        own_method (bool, optional): [use own method to calculate inv or norm, etc]. Defaults to False.
        test (bool, optional): [show testing information]. Defaults to False.

    Raises:
        ValueError: [The spetrum of Iteration Matrix must be smaller than 1]

    Returns:
        [x]: [solution of Ax=b]
    
    Author: Junno
        
    Date: 2022-04-29
    """ 
    assert A.shape[0] == b.shape[0]
    # use np method or your own
    Inv = np.linalg.inv if not own_method else Matrix_Solutions.Primary_Row_Transformation_Inv
    SVD = np.linalg.svd if not own_method else Matrix_Solutions.svd
    # substract D,L,U from A
    D, L, U = np.diag(np.diag(A)), np.tril(A, k=-1), np.triu(A, k=1)
    # calculate B and f
    B = -Inv(D)@(L+U)
    f = Inv(D)@b
    # initial x for iteration
    if x_0 == None:
        x_0 = np.zeros_like(b)
    
    # get eig of B
    u, s, v = SVD(B)
    # spectrum of B
    eig_B = np.diag(s)
    max_eig_B = np.max(eig_B)
    # check convergence condition
    if max_eig_B > 1:
        print('The spetrum of Iteration Matrix must be smaller than 1')
        raise ValueError
        
    # matrix norm
    if norm == 'infty':
        F = Matrix_Solutions.m_infinite_norm
    elif norm == 'L1':
        F = Matrix_Solutions.m_L1_norm

    q = F(B)
    if test:
        print('q =', q)
    
    # iteration times
    N = int(np.ceil(np.log10(1/eps)*np.log(10) /
            (-np.log(max_eig_B)))) if N == 0 else N
    
    # iteration step
    while N:
        x = B@x_0+f
        if q < 1 and q/(1-q)*F(x-x_0) < eps:
            break
        x_0 = x
        N -= 1

    if test:
        print("The speed of iteration is:{:.2f}".format(-np.log(max_eig_B)))
        
    return x

test example

A=np.array([[-1.35 ,  0.209, -0.218, -0.123],
       [ 0.179,  1.616, -0.104, -0.497],
       [ 0.333,  0.278,  2.417, -0.226],
       [ 1.128,  0.346,  0.104,  1.237]],dtype=np.float32).reshape((4,4))
x=np.random.randn(A.shape[0],1)
print(x)
>>>
[[-1.872]
 [ 0.693]
 [ 0.809]
 [-1.577]]
 
b=A@x
print(b)
>>>
[[ 2.691]
 [ 1.485]
 [ 1.881]
 [-3.738]]

ans=Jacobi_iteration(A, b, N=20, norm='infty', own_method=False, test=False)
print(ans)
>>>
[[-1.872]
 [ 0.693]
 [ 0.809]
 [-1.577]]

Matrix_Solutions.m_L1_norm(ans-x)
>>>
4.621872773391544e-09

Gauss-Seidel迭代法

• Jacobi迭代的第k步可以表示为： x i k + 1 = 1 a i i ( b i − ∑ j = 1 , j ≠ i n a i j x j k ) x_i^{k+1}=frac{1}{a_{ii}}(b_i-sum_{j=1,j neq i}^n a_{ij}x_j^{k}) xik+1​=aii​1​(bi​−j=1,j​=i∑n​aij​xjk​) 在计算k+1步中，要用到所有第k步的结果，但是在计算第i个x的第k+1步之前，新计算出的 x 1 k + 1 , x 1 k + 1 , . . . , x i − 1 k + 1 x_1^{k+1},x_1^{k+1},...,x_{i-1}^{k+1} x1k+1​,x1k+1​,...,xi−1k+1​没有被用到。 Gauss-Seidel迭代法便是用上这些新计算出的结果来参与计算。表示如下：
  x i k + 1 = 1 a i i ( b i − ∑ j = 1 i − 1 a i j x j k + 1 − ∑ j = i + 1 n a i j x j k ) x_i^{k+1}=frac{1}{a_{ii}}(b_i-sum_{j=1}^{i-1} a_{ij}x_j^{k+1}-sum_{j=i+1}^n a_{ij}x_j^{k}) xik+1​=aii​1​(bi​−j=1∑i−1​aij​xjk+1​−j=i+1∑n​aij​xjk​)
  写成矩阵形式，即：
  D x k + 1 = b − L x k + 1 − U x k ( D + L ) x k + 1 = b − U x k x k + 1 = − ( D + L ) − 1 U x k + ( D + L ) − 1 b Dx^{k+1}=b-Lx^{k+1}-Ux^{k}\ (D+L)x^{k+1}=b-Ux^{k}\ x^{k+1}=-(D+L)^{-1}Ux^k+(D+L)^{-1}b Dxk+1=b−Lxk+1−Uxk(D+L)xk+1=b−Uxkxk+1=−(D+L)−1Uxk+(D+L)−1b
  按照 x = B x + f x=Bx+f x=Bx+f的格式，即有：
  { B = − ( D + L ) − 1 U f = ( D + L ) − 1 b begin{cases} B=-(D+L)^{-1}U\ f=(D+L)^{-1}b\ end{cases} {B=−(D+L)−1Uf=(D+L)−1b​

代码实现

# 整合Jacobi和Gauss_Seidel到一起
def Serial_iteration(A, b, N=0, eps=1e-4, x_0=None, iter_way='Gauss_Seidel', norm='infty', own_method=False, test=False):
    """Serial_iteration:Jacobi and Gauss_Seidel

    Args:
        A ([np.darray]): [A]
        b ([np.darray]): [b]
        N (int, optional): [iteration time, it will give a suitable one when N is 0]. Defaults to 0.
        eps ([float], optional): [error threshold]. Defaults to 1e-4.
        x_0 ([np.darray], optional): [initial x for iteration]. Defaults to None.
        norm (str, optional): [iteration method,'Gauss_Seidel','Jacobi']. Defaults to 'Gauss_Seidel'.
        norm (str, optional): [norm method to calculate fitting error,'infty','L1']. Defaults to 'infty'.
        own_method (bool, optional): [use own method to calculate inv or norm, etc]. Defaults to False.
        test (bool, optional): [show testing information]. Defaults to False.

    Raises:
        ValueError: [The spetrum of Iteration Matrix must be smaller than 1]

    Returns:
        [x]: [solution of Ax=b]
    
    Author: Junno
        
    Date: 2022-04-29
    """ 
    assert A.shape[0] == b.shape[0]
    # use np method or your own
    Inv = np.linalg.inv if not own_method else Matrix_Solutions.Primary_Row_Transformation_Inv
    SVD = np.linalg.svd if not own_method else Matrix_Solutions.svd
    # substract D,L,U from A
    D, L, U = np.diag(np.diag(A)), np.tril(A, k=-1), np.triu(A, k=1)
    # calculate B and f
    if iter_way == 'Jacobi':
        B = -Inv(D)@(L+U)
        f = Inv(D)@b
    elif iter_way == 'Gauss_Seidel':
        B = -Inv(D+L)@U
        f = Inv(D+L)@b
    # initial x for iteration
    if x_0 == None:
        x_0 = np.zeros_like(b)
    
    # get eig of B
    u, s, v = SVD(B)
    # spectrum of B
    eig_B = np.diag(s)
    max_eig_B = np.max(eig_B)
    # check convergence condition
    if max_eig_B > 1:
        print('The spetrum of Iteration Matrix must be smaller than 1')
        raise ValueError
        
    # matrix norm
    if norm == 'infty':
        F = Matrix_Solutions.m_infinite_norm
    elif norm == 'L1':
        F = Matrix_Solutions.m_L1_norm

    q = F(B)
    if test:
        print('q =', q)
    
    # iteration times
    N = int(np.ceil(np.log10(1/eps)*np.log(10) /
            (-np.log(max_eig_B)))) if N == 0 else N

    # iteration step
    while N:
        x = B@x_0+f
        if q < 1 and q/(1-q)*F(x-x_0) < eps:
            print('last iteration error: ',q/(1-q)*F(x-x_0))
            break
        x_0 = x
        N -= 1

    if test:
        print("The speed of iteration is: {:.2f}".format(-np.log(max_eig_B)))
        
    return x

test example

# use example above
ans=Serial_iteration(A, b, N=0, eps=1e-4, iter_way='Gauss_Seidel',norm='infty', own_method=False, test=True)
print(ans)
>>>
q = 0.41703337
last iteration error:  4.200151608806934e-05
The speed of iteration is: 0.96
[[-1.873]
 [ 0.694]
 [ 0.809]
 [-1.576]]

超松弛迭代法SOR

• 对于Gauss-Seidel迭代法， x i k + 1 = 1 a i i ( b i − ∑ j = 1 i − 1 a i j x j k + 1 − ∑ j = i + 1 n a i j x j k ) x_i^{k+1}=frac{1}{a_{ii}}(b_i-sum_{j=1}^{i-1} a_{ij}x_j^{k+1}-sum_{j=i+1}^n a_{ij}x_j^{k}) xik+1​=aii​1​(bi​−j=1∑i−1​aij​xjk+1​−j=i+1∑n​aij​xjk​)
  引入 x i k x_i^k xik​，有
  { x i k + 1 = x i k + Δ x i Δ x i = 1 a i i ( b i − ∑ j = 1 i − 1 a i j x j k − ∑ j = i n a i j x j k ) begin{cases} x_i^{k+1}=x_i^k+Delta x_i\ Delta x_i=frac{1}{a_{ii}}(b_i-sum_{j=1}^{i-1} a_{ij}x_j^{k}-sum_{j=i}^n a_{ij}x_j^{k})\ end{cases} {xik+1​=xik​+Δxi​Δxi​=aii​1​(bi​−∑j=1i−1​aij​xjk​−∑j=in​aij​xjk​)​
  引入松弛因子 ω omega ω， x i k + 1 = x i k + ω a i i ( b i − ∑ j = 1 i − 1 a i j x j k − ∑ j = i n a i j x j k ) large x_i^{k+1}=x_i^k+frac{omega}{a_{ii}}(b_i-sum_{j=1}^{i-1} a_{ij}x_j^{k}-sum_{j=i}^n a_{ij}x_j^{k}) xik+1​=xik​+aii​ω​(bi​−∑j=1i−1​aij​xjk​−∑j=in​aij​xjk​)，选取合适的松弛因子就可以加速收敛过程，上式即为SOR的表达公式。

• 推到SOR的矩阵形式：

  • a i i x i k + 1 = ( 1 − ω ) a i i x i k + ω ( b i − ∑ j = 1 i − 1 a i j x j k − ∑ j = i + 1 n a i j x j k ) a_{ii}x_i^{k+1}=(1-omega)a_{ii}x_i^{k}+omega(b_i-sum_{j=1}^{i-1} a_{ij}x_j^{k}-sum_{j=i+1}^n a_{ij}x_j^{k}) aii​xik+1​=(1−ω)aii​xik​+ω(bi​−j=1∑i−1​aij​xjk​−j=i+1∑n​aij​xjk​)
  • D x k + 1 = ( 1 − ω ) D x k + ω ( b − L x k + 1 − U x k ) Dx^{k+1}=(1-omega)Dx^{k}+omega(b-Lx^{k+1}-Ux^{k}) Dxk+1=(1−ω)Dxk+ω(b−Lxk+1−Uxk)
  • x k + 1 = ( D + ω L ) − 1 [ ( 1 − ω ) D − ω U ] x k + w ( D + ω L ) − 1 b x^{k+1}=(D+omega L)^{-1}[(1-omega)D-omega U]x^k+w(D+omega L)^{-1}b xk+1=(D+ωL)−1[(1−ω)D−ωU]xk+w(D+ωL)−1b

• 充分条件证明：

未完待续