线性代数Python计算：消元法与矩阵初等变换

对线性方程组
$$\{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=b_2\\ \cdots \cdots \cdots \\ a_{m1}{x}_1+a_{m2}{x}_2+\cdots +a_{mn}{x}_n=b_n\end{array}$$
交换其中的两个方程（或交换方程组的两个未知量位置）、用非零常数乘一个方程、用非零常数乘一个方程加到另一个方程上，得到的方程组与原方程组是同解的。用这三种操作，消掉方程组的各方程中的一些未知量，得到方程组的解的过程称为消元法。令系数矩阵$\mathbf{A}=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \cdots & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right)$，未知数矩阵$\mathbf{x}=\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right)$，常数矩阵$\mathbf{b}=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right)$则线性方程组可表示为
$$\mathbf{Ax}=\mathbf{b}.$$
线性方程的消元解法形式化地可通过对增广矩阵
$$(\mathbf{A},\mathbf{b})=\left(\begin{array}{ccccc}{a}_{11}& a_{12}& \cdots & a_{1n}& b_1\\ a_{21}& a_{22}& \cdots & a_{2n}& b_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}& b_n\end{array}\right)$$
作有限次的第一、二、三种初等变换（交换两行或两列、用非零数乘一行、用一个数乘一行加到另一行）得到一个行最简阶梯形矩阵。
消元法解线性方程组通常分两个过程：消元过程和回代过程。先看消元过程：
从$i=1$开始，只要系数矩阵对应部分（竖线左边部分）第$i$行至最后一行内有非零行就作如下操作：若$a_{ii}=0$，先在$\left(\begin{array}{c}{a}_{ii}\\ a_{i+1i}\\ ⋮\\ a_{mi}\end{array}\right)$中找第一个非零元素，譬如第$a_{ji}\ne 0$，交换第$i$行与第$j$行。否则再在$\left(\begin{array}{cccc}{a}_{ii+1}& a_{ii+2}& \cdots & a_{in}\end{array}\right)$中找第一个非零元素。找到了，譬如$a_{ij}\ne 0$，交换第$i$列与第$j$列。否则，意味着第$i$行元素全为零，将此行与以下的某一非零行（其中的首非零元素列标$j$必不小于$i$）交换第$i$列与第$j$列，使得$a_{ii}\ne 0$。用$1/a_{ii}$乘第$i$行，使得$a_{ii}=1$。对每一个$j>i$，用$-a_{ij}$乘第$i$行加到第$j$行，使得$\left(\begin{array}{c}{a}_{i+1i}\\ a_{i+2i}\\ ⋮\\ a_{mi}\end{array}\right)$全为零。$i$自增1，进入下一轮消元。消元过程结束时，得到如下的行阶梯矩阵:
$$\left(\begin{array}{ccccccc}1& a_{12}& \cdots & a_{1r}& \cdots & a_{1n}& b_1\\ 0& 1& \cdots & a_{2r}& \cdots & a_{2n}& b_2\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 1& \cdots & a_{rn}& b_r\\ 0& 0& \cdots & 0& \cdots & 0& b_{r+1}\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 0& \cdots & 0& b_m\end{array}\right)$$
表示成Python代码如下：

import numpy as np                              #导入numpy
def rowLadder(A, m, n):
    rank=0                                      #非零行数初始化为0
    zero=m                                      #全零行首行下标
    i=0                                         #当前行号
    order=np.array(range(n))                    #未知量顺序
    while i<min(m,n) and i<zero:                #自上向下处理每一行
        flag=False                              #A[i,i]非零标志初始化为False
        index=np.where(abs(A[i:,i])>1e-10)      #当前列A[i,i]及其以下的非零元素下标
        if len(index[0])>0:                     #存在非零元素
            rank+=1                             #非零行数累加1
            flag=True                           #A[i,i]非零标记
            k=index[0][0]                       #非零元素最小下标
            if k>0:                             #若非第i行
                P1(A,i,i+k)                     #交换第i，k+i行
        else:                                   #A[i:,i]内全为0
            index=np.where(abs(A[i,i:n])>1e-10) #当前行A[i,i:n]的非零元素下标           
            if len(index[0])>0:                 #存在非零元素，交换第i，k+i列
                rank+=1
                flag=True
                k=index[0][0]
                P1(A,i,i+k,row=False)           #列交换
                order[[i, k+i]]=order[[k+i, i]] #跟踪未知量顺序
        if flag:                                #A[i,i]不为零，A[i+1:m,i]消零
            P2(A,i,1/A[i,i])
            for t in range(i+1, zero):
                P3(A,i,t,-A[t,i])
            i+=1                                #下一行
        else:                                   #将全零元素行交换到矩阵底部
            P1(A,i,zero-1)
            zero-=1                             #更新全零行首行下标
    return rank, order

程序的第2~32行定义函数rowLadder，其三个参数分别为A表示某线性代数的增广矩阵，m和n分别表示方程组拥有的方程个数（增广矩阵的行数）和方程组的未知量个数。第3~6行分别设置4个变量：行阶梯矩阵非零行数rank（初始化为0），零行起始行号（自该行起后全为零行）zero（初始化为m，即末行后），当前行号i（初始化为0，即第1行）和用来跟踪未知量顺序的order（初始化为range(n)，即自然顺序 { 0 , 1 , ⋯   , n − 1 } \{0,1,\cdots,n-1\} {0,1,⋯,n−1}）。
第7~31行的while循环自上而下逐行处理，直至当前行遇到全零行，即i$\ge$zero（i<min(m,n)是用来限制“高”矩阵，即m>n的情形，A[i,i]的列标超出允许范围的“哨兵”）。其中，第8~23行确定A[i,i]$\ne 0$。为此，第8行设置标志flag，初始化为Fales，该部分操作完成，若A[i,i]不等于零，flag为True，否则为Fales。按上节的形式化阐述知，该部分操作分两步：先查看第i列从A[i,i]起至A[m-1,i]为止的各元素中是否有非零元。若有，则将其所在行与第i行交换，确定A[i,i]非零。方法是先在第9行调用numpy的where函数寻求A[i:,i]中不等于零（abs(A[i:,i])>1e-10，即A[i:,i]中元素绝对值大于$\frac{1}{10^{10}}$）的元素所在位置下标。
$$\text{where(bool)}$$
本质上是一个在数组中按条件查找特定元素的函数。参数bool是一个与待操作数组相关的布尔表达式，如此处的abs(A[i:,i])>1e-10。该函数的返回值是一个数组，其中的元素是指示数组A中满足条件bool的元素下标。注意A是一个2维数组，故得到满足条件的元素下标也存储为一个2维数组，赋予index。若A[i:,i]中有非零元素，则index非空，此由第10行的if语句检测之。若是，则意味着可将确定多一非零行，第11行将rank自增1。且第12行将flag设为True。A[i:,i]中非零元素的最小下标应存储于index[0,0]处，第13行将其置为k。若k为0意味着A[i,i]本身非零。否则第14行的if语句测得此情形，即A[i+k]非零。第13行调用我们在博文《生成初等矩阵》中定义的第一种初等变换函数P1，交换A的第i行和第i+k行。
若第10行中测得A[i:,i]中无非零元，则需要考察第i行中自A[i,i]至A[i,n-1]为止的元素中是否有非零元，若有，则找出列标最小者，将其所在列于第i列交换。第17行调用numpy的where函数计算A[i,i:n]中诸非零元下标，赋予index。第18行测得index中确有非零元下标，在19、20行更新rank和flag后，第21行取得其中最小下标k。第22行调用P1函数交换A的第i列和第i+k列。这样，无论是执行了第11~15行的交换两行的操作还是第19~23行的交换两列的操作，都使得A[i,i]非零（此时flag为True，rank增加了1）。否则（即第10行的检测和第18行的检测均不成功），意味着A[i]是一个零行（flag保持为False，rank保持不变）。
第14~31行的if-else语句通过检测flag的值决定是将A[i,i]以下元素清零（第25~28行）还是将零行A[i]换到矩阵底部。若flag为True，第25行调用第二种初等变换函数P2（见博文《生成初等矩阵》），将A[i,i]置为1。第26~27行的for循环调用第三种初等变换P3函数（见博文《生成初等矩阵》）将A[i:,i]中A[i,i]以下的元素清零，然后第28行i自增1准备考察下一行。若flag为False，第30行调用第一种初等变换函数P1交换A的第i行和第zero-1行，第31行zero自减1。
while循环结束，A成为行阶梯矩阵。rank记录下该阶梯矩阵的非零行数，order跟踪未知量顺序的变化（执行过第22行交换两列的操作），第32行将两者作为返回值返回。
例1 运用上述定义的rowLadder函数，将线性方程组$$\{\begin{array}{l}{x}_1-x_2-x_3=2\\ 2x_1-x_2-3x_3=1\\ 3x_1+2x_2-5x_3=0\end{array}$$
的增广矩阵变换为行阶梯矩阵。

import numpy as np                  #导入numpy
A=np.array([[1,-1,-1],              #系数矩阵
            [2,-1,-3],
            [3,2,-5]],dtype='float')
b=np.array([2,1,0])                 #常数矩阵
B=np.hstack((A,b.reshape(3,1)))     #增广矩阵
rank, order=rowLadder(B,3,3)        #计算行阶梯阵
print(B)
print(rank)
print(order)

程序的第2~5行分别设置系数矩阵$\mathbf{A}$和常数矩阵$\mathbf{b}$为A和b。第6行调用numpy的hstack函数将A和b横向地连接成增广矩阵B。第7行调用上述程序定义的rowLadder函数，传递B，3（3个方程），3（3个未知量），将B变换成行阶梯矩阵，返回其中的非零行数rank和变换后未知量的顺序order。运行程序，输出

[[ 1. -1. -1.  2.]
 [ 0.  1. -1. -3.]
 [ 0.  0.  1.  3.]]
3
[0 1 2]

前三行为B在变换后得到的行阶梯阵$\left(\begin{array}{cccc}1& -1& -1& 2\\ 0& 1& -1& -3\\ 0& 0& 1& 3\end{array}\right)$，它对应所得到的同解方程组$\{\begin{array}{l}{x}_1-x_2-x_3=2\\ x_2-x_3=-3\\ x_3=3\end{array}$。第4行输出3，意味着阶梯阵中非零行数为3，第5行输出自然序列0，1，2意味着变换过程中没有作列交换。
线性方程组的消元解法第二个过程为回代。回代过程的操作和消元过程中的几乎一致，从$i=r$（=rank）开始依次用$-a_{i-1,i}$乘以第$i$行加到前$i-1$行，$i$自减1，直至$i=0$。结束时，得到最简行阶梯矩阵：

import numpy as np                              #导入numpy
def simplestLadder(A,rank):
    for i in range(rank-1,0,-1):                #自下而上逐行处理
        for j in range(i-1, -1,-1):             #自下而上将A[i,i]上方元素消零
            P3(A,i,j,-A[j,i])

操作应从行阶梯阵的最后非零行（rank）开始，逐行地将A[i,i]（=1）上方的元素通过第三种变换，即用-A[j,i]乘以第i行加到第j行上，将A[i,j]消为零。最后得到一个最简行阶梯矩阵。程序的第2~5行定义的simplestLadder函数，完成这样的操作。参数A中存储着一个行阶梯矩阵，rank表示A的非零行数。函数体中仅含一个嵌套的for循环：外层循环（第3~5行）是对A的前rank个非零行自下向上（i从rank减少到1）扫描，内层循环（第4~5行）调用第三种初等变换函数P3（见博文《生成初等矩阵》），完成将第i列中A[i,i]（=1）上方的元素A[j,i]化为0（j从i-1减少到0）。
例2 下列代码对例1中的方程组完成消元后继续回代。

import numpy as np                              #导入numpy
np.set_printoptions(suppress=True)              #设置数组输出精度
A=np.array([[1,-1,-1],                          #系数矩阵
            [2,-1,-3],
            [3,2,-5]],dtype='float')
b=np.array([2,1,0])                             #常数矩阵
B=np.hstack((A,b.reshape(3,1)))                 #增广矩阵
rank, order=rowLadder(B,3,3)                    #转换为行阶梯阵
simplestLadder(B,rank)                          #转换为最简行阶梯阵
print(B)

程序的第1~8行几乎和前一段程序的代码一样，此时B变换为行阶梯阵$\left(\begin{array}{cccc}1& -1& -1& 2\\ 0& 1& -1& -3\\ 0& 0& 1& 3\end{array}\right)$。第9行调用上述程序定义的simplestLadder函数，对B作变换。运行程序输出

[[1. 0. 0. 5.]
 [0. 1. 0. 0.]
 [0. 0. 1. 3.]]

即回代变换后的矩阵B为$\left(\begin{array}{cccc}1& 0& 0& 5\\ 0& 1& 0& 0\\ 0& 0& 1& 3\end{array}\right)$，对应同解方程组$\{\begin{array}{l}{x}_1=5\\ x_2=0\\ x_3=3\end{array}$。这就是原方程组的解。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
代码诚可贵，原理价更高。若为AI学，读正版书好。