GaussSeidel迭代 Jacobi迭代及其收敛性

迭代法求解线性方程组

目的是对于一个可逆矩阵$A$，利用迭代法求解

$$Ax=b.$$

$GaussSeidel$ 和 $Jacobi$ 都是将 $A$ 矩阵分成 $D+L+U$ 矩阵的形式，构造出$x=Gx+d$的迭代矩阵，来迭代求解出$x$，理想情况下当 $x==Gx+d$ 时算法收敛，且此时$x=A^{-1}b$。其中$DLU$ 分别为 $A$ 矩阵的对角线矩阵，下三角矩阵和上三角矩阵。

在上述中介绍到两种方法都是通过定义迭代矩阵 $x=Gx+d$ 来求解 $x$ 的。因此接下来分别介绍，GaussSeidel 和Jacobi分别时怎样定义迭代矩阵的。

$GaussSeidel$

公式推导

矩阵形式

$$A=D+L+U$$
$$Ax=b$$
$$(D+L+U)x=b$$
$$(D+L)x=-Ux+b$$
$$x=-{(D+L)}^{-1}Ux+{(D+L)}^{-1}b$$
$$⟹x=Gx+d$$
由上述推导可知，在$GaussSeidel$中 $G=-{(D+L)}^{-1}U$， $d={(D+L)}^{-1}b$

非矩阵形式

$$x_i^k=\frac{1}{a_{ii}}(b-\sum _{0<j<i}{a}_{ij}{x}_j^k-\sum _{j>i}{a}_{ij}{x}_j^{k-1}))$$
其中 $x_i^k$ 代表 $x$ 的第$i$个分量在第$k$次迭代中的结果

代码示例

使用非矩阵形式实现在稀疏矩阵csr_matrix下实现的$GaussSeidel$。

# A csr_matrix 
import sys
import copy
import numpy as np
import networkx as nx
from scipy.sparse import spdiags, tril, triu, coo_matrix, csr_matrix

def gauss_seidel(A, b, x0 = None, maxiter=20):
    # x0 = spsolve(A, b)
    if x0 is None:
        x0 =  np.array( [0.0]*A.shape[0]  )
    
    for _ in range(maxiter):
        n, _, indptr, indices, data = A.shape[0], A.nnz, A.indptr, A.indices, A.data
        for row in range(n):
            aii, sumi = 0.0, 0.0
            for j in range(indptr[row], indptr[row+1]):
                if row == indices[j]:
                    aii = 1.0/data[j]
                else:
                    sumi += data[j] * x0[ indices[j] ]
            x0[row] = (aii * (b[row] - sumi) )
    return x0
if __name__ == '__main__':  
    A = csr_matrix(([10,-2,-1, -2,10,-1, -1,-2,5], ([0,0,0,1,1,1,2,2,2], [0,1,2,0,1,2,0,1,2]) ) , shape=(3,3), dtype=np.float32)
    x = np.ones(A.shape[0])
    b = A.dot(x)
    xk = gauss_seidel(A, b)
    print(x)
    print(xk)

输出如下

$Jacobi$

公式推导

矩阵形式

$$A=D+L+U$$
$$Ax=b$$
$$(D+L+U)x=b$$
$$Dx=-（L+U)x+b$$
$$x=-D^{-1}(L+U)x+D^{-1}b$$
$$⟹x=Gx+d$$
由上述推导可知，在Jacobi中 $G=-D^{-1}(L+U)$， $d=D^{-1}b$

非矩阵形式

$$x_i^k=\frac{1}{a_{ii}}(b-\sum _{j!=i}{a}_{ij}{x}_j^{k-1})$$

代码示例

使用矩阵形式实现在稀疏矩阵csr_matrix下实现的$Jacobi$。

# A csr_matrix
import sys
import copy
import numpy as np
from scipy.sparse import spdiags, tril, triu, coo_matrix, csr_matrix
def jacobi(A, b, x0 = None, maxiter=20):
    if x0 is None:
        x0 =  np.array( [0.0]*A.shape[0]  )

    D = spdiags(A.diagonal(), np.array([0]), A.shape[0], A.shape[1]) # default with <class 'scipy.sparse.dia.dia_matrix'> 

    L = tril(A) - D
    U = triu(A) - D 

    dia = generate_dia_inv(D)
    B = -( dia ).dot(L+U)
    f = ( dia ).dot(b)

    for _ in range(maxiter):
        x0 = B.dot(x0) + f
    return x0
def generate_dia_inv(mat):
    res_dia = mat.diagonal()
    res_dia = 1.0/res_dia
    res_mat = csr_matrix((res_dia.shape[0], res_dia.shape[0]))
    res_mat.setdiag(res_dia)
    return res_mat
if __name__ == '__main__':  
    A = csr_matrix(([10,-2,-1, -2,10,-1, -1,-2,5], ([0,0,0,1,1,1,2,2,2], [0,1,2,0,1,2,0,1,2]) ) , shape=(3,3), dtype=np.float32)
    x = np.ones(A.shape[0])
    b = A.dot(x)
    xk = jacobi(A, b)
    print(x)
    print(xk)

输出如下

收敛性

因为两种方法的根本思想是一致的，都是通过迭代公式 $x=Gx+d$ 来求解 $x$。所以直接通过此迭代公式来看两种方法收敛的条件。
令误差向量 $e^i=x^i-x$，其中 $x$ 是 $Ax=b$的精确解， $e^i$是$GaussSeidel$或者$Jacoci$第$i$次迭代所求解。

收敛性证明
由于 $x^i=G{x}^{i-1}+d$
且 $x=Gx+d$
所以有，$e^i=G{e}^{i-1}$
且 $e^{i-1}=G{e}^{i-2}$
$...$
因此，$e^i=G^i{e}^0$
当迭代法收敛时，就意味着 $\underset{i\to \infty }{\lim }{e}^i\to 0$
由特征值和特征向量的定义可知$G{\zeta }_i={\lambda }_i{\zeta }_i$，且特征向量之间线性无关
因此$e^0=\Sigma {\delta }_k{\zeta }_k$, 可以这样表示，且 ${\delta }_k$是定值
则$e^1=G{e}^0=G\sum {\delta }_k{\zeta }_k=\sum {\delta }_{k}G{\zeta }_k=\sum {\delta }_k{\lambda }_k{\zeta }_k$，由 $G{\zeta }_i={\lambda }_i{\zeta }_i$得。
则$e^2=\sum {\delta }_k{\lambda }_k^2{\zeta }_k$ , $e^3=\sum {\delta }_k{\lambda }_k^3{\zeta }_k$ …
那么 $\underset{i\to \infty }{\lim }{e}^i=\underset{i\to \infty }{\lim }\sum {\lambda }_k^i{\delta }_k{\zeta }_k$
其中${\delta }_k{\zeta }_k$是定值
则当 $abs({\lambda }_k)>1⟹\underset{i\to \infty }{\lim }{e}^i\to \infty$
当 $abs({\lambda }_k)<1⟹\underset{i\to \infty }{\lim }{e}^i\to 0$
综上，当迭代矩阵$G$的所有特征值 $abs({\lambda }_k)<1$时， 该迭代方法收敛。