专题：线性方程组求解

清雨影的Blog
搜索文档…
README
机器学习
专题：线性方程组求解
简介
这里完整介绍一些求 Ax⃗=b⃗\bold{A}\vec{x}=\vec{b}Ax=b
 的方法，特点是A矩阵的规模比较大。 正常来说，求解线性方程组可以通过求逆矩阵，高斯消元，克拉默法则（大学线代第一节课就学这个没用的东西，不会真的有人用吧，不会吧不会吧）等等方式去求解。 但是A很大的时候这个方法就失效了，这个时候一般用迭代法，本文以介绍迭代法为主。
这里解的方程组一定要正定，负定的矩阵解起来注定是失稳的。如果是正半定，或者说奇异矩阵（存在为0的特征值），那么解将会在某个子空间里而不唯一。
其实从二次型的视角就很容易看出：
（请注意，本文中的推导仅仅是验证实验用，性能并非最优）
Jacobi 迭代与 Gauss-Seidel
Jacobi 迭代
推导过程：
我们把A矩阵分尸三块：A=D−L−UA = D - L - UA=D−L−U
​
其中，D是对角线，L是上三角取负，U是下三角取负。
这样，Ax=bAx=bAx=b
 就可以写成 (D−L−U)x=b(D-L-U)x=b(D−L−U)x=b
 
移动一下：Dx=(L+U)x+bDx=(L+U)x + bDx=(L+U)x+b
 
由于D是对角矩阵，可无成本求逆：x=D−1(L+U)x+D−1bx=D^{-1}(L+U)x+D^{-1}bx=D−1(L+U)x+D−1b
 
于是我们利用下列公式迭代更新：
x(k+1)=D−1(L+U)x(k)+D−1bx^{(k+1)} = D^{-1}(L+U)x^{(k)} + D^{-1}bx(k+1)=D−1(L+U)x(k)+D−1b
即可求得解，至于为什么是这样，可看收敛性分析。
%accordion% Numpy实现 %accordion%
1
def jacobi_iteration(A,b,init_x,iters=10):
2
    dim = init_x.shape[0]
3
    x = init_x
4
    result = np.zeros((iters+1,dim))
5
    result[0,:] = x.reshape(-1)
6
​
7
    diag_A = np.diag(A)
8
    D_inv = np.diag(1.0/diag_A)
9
    LU = np.diag(diag_A) - A
10
    B = D_inv @ LU
11
    D_inv_b = D_inv @ b
12
​
13
    for i in range(iters):
14
        x = B @ x + D_inv_b
15
        result[i+1,:] = x.reshape(-1)
16
    return result
Copied!
%/accordion%
Gauss-Seidel
Jacobi迭代中，x(k+1)=D−1(L+U)x(k)+D−1bx^{(k+1)} = D^{-1}(L+U)x^{(k)} + D^{-1}bx(k+1)=D−1(L+U)x(k)+D−1b
， 这里的 x(k+1)x^{(k+1)}x(k+1)
是向量，Jacobi迭代更新的时候是并行更新的。 但是Gauss-Seidel则是一个个元素更新的。
更新x1(k+1)x^{(k+1)}_{1}x1(k+1)​
的时候会用到x0(k+1)x^{(k+1)}_{0}x0(k+1)​
，同理，更新x2(k+1)x^{(k+1)}_{2}x2(k+1)​
的时候会用到x0(k+1)x^{(k+1)}_{0}x0(k+1)​
和x1(k+1)x^{(k+1)}_{1}x1(k+1)​
。
（然而并没什么卵用，不必Jacobi快多少，还不能并行计算了。。。）
%accordion% Numpy实现 %accordion%
1
def gauss_seidel_iteration(A,b,init_x,iters=10):
2
    dim = init_x.shape[0]
3
    x = init_x
4
    result = np.zeros((iters+1,dim))
5
    result[0,:] = x.reshape(-1)
6
​
7
    diag_A = np.diag(A)
8
    D_inv = np.diag(1.0/diag_A)
9
    LU = np.diag(diag_A) - A
10
    B = D_inv @ LU
11
    D_inv_b = D_inv @ b
12
​
13
    for i in range(iters):
14
        for d in range(dim):
15
            x[d] = B[d,:] @ x + D_inv_b[d]
16
        result[i+1,:] = x.reshape(-1)
17
    return result
Copied!
%/accordion%
收敛性分析
如果矩阵A每一行的“除了对角线元素之外，其它元素的和”都要小于对角线元素的话，则这个矩阵一定可以被Jacobi迭代求解（我也不知为何如此，不过应该可以从这个条件推出谱半径小于1）。
还有一个就是B矩阵 （B=D−1(L+U)B=D^{-1}(L+U)B=D−1(L+U)
） 的谱半径，谱半径等于最大的特征值。每个向量都可以被分解为特征向量的和（也就是表达在某个以特征向量为基的空间），如果谱半径大于1，那么求解的时候就会发散。谱半径不仅要小于1，而且要越小越好，这样收敛才够迅速。
而为什么我们关心B的谱半径呢？
我们站在神的视角，已经知道了真实的解x，也就是说x绝对满足： x=Bx+D−1bx=Bx + D^{-1}bx=Bx+D−1b
，那么不妨设置误差 e(k)=x(k)−xe^{(k)} = x^{(k)} - xe(k)=x(k)−x
​
每一次迭代的时候，我们看作：e(k+1)+x=B(e(k)+x)+D−1be^{(k+1)} + x = B(e^{(k)}+x) + D^{-1}be(k+1)+x=B(e(k)+x)+D−1b
​
所以有：e(k+1)+x=Bx+D−1b+Be(k)e^{(k+1)} + x = Bx + D^{-1}b + Be^{(k)}e(k+1)+x=Bx+D−1b+Be(k)
​
我们知道 x=Bx+D−1bx = Bx + D^{-1}bx=Bx+D−1b
 所以消去等号两侧的两项，得到：
​e(k+1)=Be(k)e^{(k+1)} = Be^{(k)}e(k+1)=Be(k)
​
所以，很直观的看到，只有谱半径小于1，才能保证误差不断减少，即所谓收敛。
梯度下降
梯度下降迭代法的思路是优化 x=argminz∥Az−b∥x = argmin_z \| Az - b \|x=argminz​∥Az−b∥
​
所以Loss Function就是 L=xTAx−bTxL = x^TAx - b^TxL=xTAx−bTx
​
∂L∂x=Ax−b\frac{\partial{L}}{\partial{x}}=Ax-b∂x∂L​=Ax−b
梯度方向有了，下面就是学习率了，和机器学习任务不一样，这里的学习率是可以有最优解的， 原理大概是如果沿着梯度方向前后走，并且记录路径上L的大小，那么我们会得到一个抛物线，抛物线是有极小值的啊！！！！
∂∂α(f(α)TAf(α)−bTf(α))=0\frac{\partial}{\partial{\alpha}} (f(\alpha)^TAf(\alpha)-b^Tf(\alpha)) = 0∂α∂​(f(α)TAf(α)−bTf(α))=0
其中 f(α)=x−α(Ax−b)f(\alpha) = x - \alpha (Ax-b)f(α)=x−α(Ax−b)
​
或者说，走到新的点以后，其梯度应该和当前的梯度正交，这两种表述是等价的：
最优解是：
α=xTxxTAx\alpha = \frac{x^Tx}{x^TAx}α=xTAxxTx​
过程太难打暂略。
根据这个性质，梯度下降的算法将会以互相垂直的折线路径快速逼近最优解。
%accordion% Numpy实现 %accordion%
1
def gradient_desc(A,b,init_x,iters=10):
2
    dim = init_x.shape[0]
3
    x = init_x
4
    result = np.zeros((iters+1,dim))
5
    result[0,:] = x.reshape(-1)
6
    for i in range(iters):
7
        ax = A @ x
8
        grad = ax - b
9
        flat_x = x.reshape(-1)
10
        lr = np.dot(flat_x,flat_x)/np.dot(ax.reshape(-1),flat_x)
11
        x -= grad * lr
12
        result[i+1,:] = x.reshape(-1)
13
    return result
Copied!
%/accordion%
共轭梯度法
这个方法的思路是，不走回头路，如果我每一次新的搜索空间都和之前的所有搜索空间正交，那么可以更加高效。
但是求出一个正交的子空间（比如用施密特正交）复杂度他娘的是 O(n3)O(n^3)O(n3)
 啊。
Warm Start
在解微分方程（组）的时候，由于上一个时间片的状态和这个时间片的状态相差并不多，所以可以用上一个时间片的状态作为这个时间片的初始值输入，这样可以大大增加求解的效率。
MultiGrid
References
测试代码
%accordion% 测试框架 %accordion%
1
import numpy as np
2
from matplotlib import pyplot as plt
3
​
4
def solver_visualizer(A:np.ndarray,b:np.ndarray,x:np.ndarray,steps:int=100):
5
    """
6
    A: 2x2 matrix
7
    b: 2x1 vector
8
    x: nx2 vectors
9
    """
10
    ax_min = x[:,0].min()
11
    ax_max = x[:,0].max()
12
    rx = ax_max - ax_min
13
    ay_min = x[:,1].min()
14
    ay_max = x[:,1].max()
15
    ry = ay_max - ay_min
16
    X,Y = np.meshgrid(
17
        np.linspace(ax_min - 0.2 * rx,ax_max + 0.2 * rx,steps),
18
        np.linspace(ay_min - 0.2 * ry,ay_max + 0.2 * ry,steps),
19
    )
20
    E = np.linalg.norm(
21
        A @ np.array([X.reshape(-1),Y.reshape(-1)]) - b, 
22
        axis=0
23
    ).reshape(
24
        X.shape
25
    )
26
    plt.figure(figsize=(10,6))
27
    plt.contourf(X,Y,E)
28
    plt.contour(X,Y,E)
29
    plt.plot(x[:,0],x[:,1],'-o')
30
    plt.show()
31
​
32
A = np.array([[1,0.2],[-0.3,2]])
33
b = np.array([0.6,0.8]).reshape(2,1)
34
init_x =  b / np.diag(A).reshape(2,1)
35
​
36
# Example
37
​
38
solver_visualizer(A,b,gauss_seidel_iteration(A,b,init_x,10))
Copied!
%/accordion%
参考文献
​An Introduction to the Conjugate Gradient Method Without the Agonizing Pain​
​共轭梯度法通俗讲义​
软件库收集
​MKL Sparse Solver Routines​
​PARDISO​
​Sparse Matrix Solvers on the GPU: Conjugate Gradients and Multigrid​
​cuSOLVER​