上节(https://www.jianshu.com/p/671a94ce586b)说到高斯消元法，今天从高斯列主元消元法开始，拓展到线性方程组的两种迭代方法：雅可比迭代法和高斯-赛德尔迭代法。同样，能力有限，希望读者指出我的问题，用代码和公式和我深入交流。

2.列主元消去法

参考的教材是《数值分析》（李庆扬等）的第148页到151页。
列主元消去法是使用于主元为0或很小的情况，思路是选取这一列中绝对值最大的值所在的行与当前行进行替换，具体可以参考教材。再写代码方面，有了第一节的东西，写列主元消去法就很简单了，只要在高斯消去法之前加上换行就行。代码如下：

def pivotTriangularMatrix(matA):
    '''
    先找主元，找到后跟高斯顺序消元操作一样
    '''
    numRows = matA.shape[0]
    for row in np.arange(0,numRows-1):
        Index = np.argmax(abs(matA[:,row]))
        matA[[row,Index],:] = matA[[Index,row],:] # 交换。也可以用np.copy来做交换
        
        # 交换行结束，开始消元
        for k in range(row+1,numRows):  # 在第row行的情况下，遍历剩下的行
            if matA[k,row] != 0:
            #第k(k小于row)行的新值等于该行减去第row行的值乘于一个系数。这个系数存在目的就是消元
                matA[k,:]=matA[k,:]-(matA[k,row]/matA[row,row])*matA[row,:]
    return matA # 此时输入的矩阵也被改变  

用教材148页的例题4来测试一下这个代码。

test_mat3 = np.array([[0.001,2.0,3.0,1.0],[-1.0,3.712,4.623,2.0],[-2.0,1.072,5.643,3.0]],dtype=float)
result = Gauss_solve(pivotTriangularMatrix(test_mat3)) # Gauss_solve函数是上一节的函数
# 得出结果
result = array([[-0.49039646],
                [-0.05103518],
                [ 0.36752025]])

得出的结果和书中给出的答案是吻合的。暂时没有找其它算例测试。

3.今日加餐

今天舍友问我下面这个代码。刚好今天没去旁听。主要问这里面提到的独特的索引方式是什么。其实就是一组布尔值。这里举个例子，具体解释在代码里面。

实现取斜对角线.png

i,n = 1,3
idx1 = np.arange(i,n+i)
'''
假设i=1,n=3
idx1 = array([1, 2, 3])
'''
idx1>n-1
'''
输出为 array([False, False,  True])，这样就拿到了True所在的索引2
那么idx1[idx1>n-1]在这里就相当于idx1[2]，这个值等于3
'''
# 举例
test_mat1 = np.array([[2,0,6],[1,4,3],[3,2,1]],dtype=float)
test_mat[idx0,idx1]
'''
测试矩阵为：
test_mat1 = array([[2., 0., 6.],
                   [1., 4., 3.],
                   [3., 2., 1.]])
输出值为 array([0., 3., 3.]) 。实现了上图中取斜对角的操作
'''

二、线性方程组的迭代求解法

1.雅可比迭代法

该处参考的是《数值分析》（李庆扬等）第6章的6.1.1和6.2.1。书上关于雅可比迭代的内容有更多细节。这里我按照我对Gradient Descent的理解来编写该处代码。
这是一个迭代的计算方法，从我的理解来看，这类方法与线性方程组直接求解法有2个较大区别：1.对于解，会有一个初始值，会从这个初始值开始按照迭代公式运算下去。2.迭代法都要有迭代结束条件，要么是迭代解与真实解之间的差小于人为设定的一个阈值，要么是达到人为设定的最大迭代次数。
以下是我的代码，参考的是《数值分析》（李庆扬等）188页2.3这个公式，这个公式是一个展开的公式，我代码里写成矩阵运算的形式，这样可以减少使用for循环。

def Jocobi(A,b,initial,delta,maxTimes):
    '''
    输入：A是系数矩阵，N阶方阵
          b是N*1列向量
          initial是解的初始值，N*1大小
          delta是人为设置的一个阈值，指到后期两次迭代之间解之间的差距已经很小了，用作迭代结束条件
          maxTimes指人为设定的最大迭代次数
          
          教材6.2.1将矩阵A分解成D,L,U，这里与之对应
    输出：迭代后的解析解
    '''
    D = np.diag(np.diag(A)) # 获得D矩阵
    L = np.tril(A,-1)       # 获得L矩阵
    U = np.triu(A,1)        # 获得U矩阵
    d = np.linalg.inv(D)    # 对D矩阵求逆
    G = -np.dot(d,L+U)      # D的逆矩阵乘于（L+U）矩阵
    f = np.dot(d,b)         # D的逆矩阵b向量
    X = np.dot(G,initial)+f # 初次的解
    for i in range(maxTimes):
        if np.linalg.norm(X - initial) > delta:
            initial = X
            X = np.dot(G,initial) + f
    return X

利用第180页的例题1来对这个函数进行验证。解的初始值设为0解。

import numpy as np
A = np.array([[8,-3,2],[4,11,-1],[6,3,12]],dtype=float)
b = np.array([[20],[33],[36]],dtype=float)
initial = np.zeros((3,1))
'''
A = array([[ 8., -3.,  2.],
           [ 4., 11., -1.],
           [ 6.,  3., 12.]])
b = array([[20.],
           [33.],
           [36.]])
initial = array([[0.],
                [0.],
                [0.]])
'''
# 迭代求解
X = Jocobi(A,b,initial,0.001,10)
'''
X = array([[3.00028157],
           [1.99991182],
           [0.99974048]])
'''

进过10次迭代后解出的结果与书上181页的结果一致。

2.高斯-赛德尔迭代法

有了雅可比迭代的代码，写高斯赛德尔迭代代码就非常简单，你可能会惊讶，这两个代码怎么是一样的。看下图可以看到这两的区别。参考教材的6.2.2节。

image.png

以下是我的代码：

def Seidel(A,b,initial,delta,maxTimes):
    '''
    输入：A是系数矩阵，N阶方阵
          b是N*1列向量
          initial是解的初始值，N*1大小
          delta是人为设置的一个阈值，指到后期两次迭代之间解之间的差距已经很小了，用作迭代结束条件
          maxTimes指人为设定的最大迭代次数
          
          教材6.2.1将矩阵A分解成D,L,U，这里与之对应
    输出：迭代后的解析解
    '''
    # 先分解
    D = np.diag(np.diag(A)) 
    L = np.tril(A,-1)
    U = np.triu(A,1)
    
    d = np.linalg.inv( D + L ) # 这里是和雅克比不同的地方
    
    G = -np.dot(d,U)
    f = np.dot(d,b)
    X = np.dot( G ,initial ) + f
    for i in range(maxTimes):
        if np.linalg.norm( X - initial ) > delta:
            initial = X
            X = np.dot( G,initial )+f       
    return X

同样用前面雅可比的算例来测试一下。

X = Seidel(A,b,initial,0.0001,10)
'''
X = array([[3.00000201],
           [1.9999987 ],
           [0.99999932]])
'''

同样的例题，在同样迭代10次的情况下，与前面雅可比迭代的结果有一点点偏差，但这都在正常范围之内。与书上的结果基本吻合。

今日小结

今天可太忙了，从早到晚。时间全靠课间挤，关于迭代法也没有考虑收敛条件的问题。这个系列对我来说还挺好玩，打算有始有终做下去。原来，写东西（还不用负责）这件事情能这么快乐。
今晚毕，明天继续。。。。。