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

对线性方程组
{a11x1+a12x2+⋯+a1nxn=b1a21x1+a22x2+⋯+a2nxn=b2⋯⋯⋯am1x1+am2x2+⋯+amnxn=bn\begin{cases}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\\\quad\quad\quad\cdots\quad\cdots\quad\cdots\quad\\a_{m1}x_1+a_{m2}x_2+\cdots+a_{mn}x_n=b_n\end{cases}⎩⎪⎪⎪⎨⎪⎪⎪⎧a11x1+a12x2+⋯+a1nxn=b1a21x1+a22x2+⋯+a2nxn=b2⋯⋯⋯am1x1+am2x2+⋯+amnxn=bn
交换其中的两个方程（或交换方程组的两个未知量位置）、用非零常数乘一个方程、用非零常数乘一个方程加到另一个方程上，得到的方程组与原方程组是同解的。用这三种操作，消掉方程组的各方程中的一些未知量，得到方程组的解的过程称为消元法。令系数矩阵A=(a11a12⋯a1na21a22⋯a2n⋮⋮⋯⋮am1am2⋯amn)\boldsymbol{A}=\begin{pmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\ a_{21}&a_{22}&\cdots&a_{2n}\\ \vdots&\vdots&\cdots&\vdots\\ a_{m1}&a_{m2}&\cdots&a_{mn}\end{pmatrix}A=⎝⎜⎜⎜⎛a11a21⋮am1a12a22⋮am2⋯⋯⋯⋯a1na2n⋮amn⎠⎟⎟⎟⎞，未知数矩阵x=(x1x2⋮xn)\boldsymbol{x}=\begin{pmatrix}x_1\\x_2\\\vdots\\x_n\end{pmatrix}x=⎝⎜⎜⎜⎛x1x2⋮xn⎠⎟⎟⎟⎞，常数矩阵b=(b1b2⋮bn)\boldsymbol{b}=\begin{pmatrix}b_1\\b_2\\\vdots\\b_n\end{pmatrix}b=⎝⎜⎜⎜⎛b1b2⋮bn⎠⎟⎟⎟⎞则线性方程组可表示为
Ax=b.\boldsymbol{Ax}=\boldsymbol{b}.Ax=b.
线性方程的消元解法形式化地可通过对增广矩阵
(A,b)=(a11a12⋯a1nb1a21a22⋯a2nb2⋮⋮⋱⋮⋮am1am2⋯amnbn)(\boldsymbol{A},\boldsymbol{b})=\left(\begin{array}{cccc:c}a_{11}&a_{12}&\cdots&a_{1n}&b_1\\a_{21}&a_{22}&\cdots&a_{2n}&b_2\\\vdots&\vdots&\ddots&\vdots&\vdots\\a_{m1}&a_{m2}&\cdots&a_{mn}&b_n\end{array}\right)(A,b)=⎝⎜⎜⎜⎛a11a21⋮am1a12a22⋮am2⋯⋯⋱⋯a1na2n⋮amnb1b2⋮bn⎠⎟⎟⎟⎞
作有限次的第一、二、三种初等变换（交换两行或两列、用非零数乘一行、用一个数乘一行加到另一行）得到一个行最简阶梯形矩阵。
消元法解线性方程组通常分两个过程：消元过程和回代过程。先看消元过程：
从i=1i=1i=1开始，只要系数矩阵对应部分（竖线左边部分）第iii行至最后一行内有非零行就作如下操作：若aii=0a_{ii}=0aii=0，先在(aiiai+1i⋮ami)\begin{pmatrix}a_{ii}\\a_{i+1i}\\\vdots\\a_{mi}\end{pmatrix}⎝⎜⎜⎜⎛aiiai+1i⋮ami⎠⎟⎟⎟⎞中找第一个非零元素，譬如第aji≠0a_{ji}\not=0aji=0，交换第iii行与第jjj行。否则再在(aii+1aii+2⋯ain)\begin{pmatrix}a_{ii+1}&a_{ii+2}&\cdots&a_{in}\end{pmatrix}(aii+1aii+2⋯ain)中找第一个非零元素。找到了，譬如aij≠0a_{ij}\not=0aij=0，交换第iii列与第jjj列。否则，意味着第iii行元素全为零，将此行与以下的某一非零行（其中的首非零元素列标jjj必不小于iii）交换第iii列与第jjj列，使得aii≠0a_{ii}\not=0aii=0。用1/aii1/a_{ii}1/aii乘第iii行，使得aii=1a_{ii}=1aii=1。对每一个j>ij>ij>i，用−aij-a_{ij}−aij乘第iii行加到第jjj行，使得(ai+1iai+2i⋮ami)\begin{pmatrix}a_{i+1i}\\a_{i+2i}\\\vdots\\a_{mi}\end{pmatrix}⎝⎜⎜⎜⎛ai+1iai+2i⋮ami⎠⎟⎟⎟⎞全为零。iii自增1，进入下一轮消元。消元过程结束时，得到如下的行阶梯矩阵:
(1a12⋯a1r⋯a1nb101⋯a2r⋯a2nb2⋮⋮⋱⋮⋱⋮⋮00⋯1⋯arnbr00⋯0⋯0br+1⋮⋮⋱⋮⋱⋮⋮00⋯0⋯0bm)\left( \begin{array}{cccccc:c} 1&a_{12}&\cdots&a_{1r}&\cdots&a_{1n}&b_1\\ 0&1&\cdots&a_{2r}&\cdots&a_{2n}&b_2\\ \vdots&\vdots&\ddots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&1&\cdots&a_{rn}&b_r\\ 0&0&\cdots&0&\cdots&0&b_{r+1}\\ \vdots&\vdots&\ddots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&0&\cdots&0&b_{m}\end{array}\right)⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛10⋮00⋮0a121⋮00⋮0⋯⋯⋱⋯⋯⋱⋯a1ra2r⋮10⋮0⋯⋯⋱⋯⋯⋱⋯a1na2n⋮arn0⋮0b1b2⋮brbr+1⋮bm⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
表示成Python代码如下：

import numpy as np                              #导入numpy
def rowLadder(A, m, n):rank=0                                      #非零行数初始化为0zero=m                                      #全零行首行下标i=0                                         #当前行号order=np.array(range(n))                    #未知量顺序while i<min(m,n) and i<zero:                #自上向下处理每一行flag=False                              #A[i,i]非零标志初始化为Falseindex=np.where(abs(A[i:,i])>1e-10)      #当前列A[i,i]及其以下的非零元素下标if len(index[0])>0:                     #存在非零元素rank+=1                             #非零行数累加1flag=True                           #A[i,i]非零标记k=index[0][0]                       #非零元素最小下标if k>0:                             #若非第i行P1(A,i,i+k)                     #交换第i，k+i行else:                                   #A[i:,i]内全为0index=np.where(abs(A[i,i:n])>1e-10) #当前行A[i,i:n]的非零元素下标           if len(index[0])>0:                 #存在非零元素，交换第i，k+i列rank+=1flag=Truek=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≥\geq≥zero（i<min(m,n)是用来限制“高”矩阵，即m>n的情形，A[i,i]的列标超出允许范围的“哨兵”）。其中，第8~{}23行确定A[i,i]≠0\not=0=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]中元素绝对值大于11010\frac{1}{10^{10}}10101）的元素所在位置下标。
where(bool)\text{where(bool)}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函数，将线性方程组{x1−x2−x3=22x1−x2−3x3=13x1+2x2−5x3=0\begin{cases}x_1-x_2-x_3=2\\ 2x_1-x_2-3x_3=1\\ 3x_1+2x_2-5x_3=0 \end{cases}⎩⎪⎨⎪⎧x1−x2−x3=22x1−x2−3x3=13x1+2x2−5x3=0
的增广矩阵变换为行阶梯矩阵。

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行分别设置系数矩阵A\boldsymbol{A}A和常数矩阵b\boldsymbol{b}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在变换后得到的行阶梯阵(1−1−1201−1−30013)\begin{pmatrix}1&-1&-1&2\\0&1&-1&-3\\0&0&1&3\end{pmatrix}⎝⎛100−110−1−112−33⎠⎞，它对应所得到的同解方程组{x1−x2−x3=2x2−x3=−3x3=3\begin{cases}x_1-x_2-x_3=2\\\quad\quad x_2-x_3=-3\\ \quad\quad\quad\quad x_3=3\end{cases}⎩⎪⎨⎪⎧x1−x2−x3=2x2−x3=−3x3=3。第4行输出3，意味着阶梯阵中非零行数为3，第5行输出自然序列0，1，2意味着变换过程中没有作列交换。
线性方程组的消元解法第二个过程为回代。回代过程的操作和消元过程中的几乎一致，从i=ri=ri=r（=rank）开始依次用−ai−1,i-a_{i-1,i}−ai−1,i乘以第iii行加到前i−1i-1i−1行，iii自减1，直至i=0i=0i=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的非零行数。函数体中仅含一个嵌套的{\bf{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变换为行阶梯阵(1−1−1201−1−30013)\begin{pmatrix}1&-1&-1&2\\0&1&-1&-3\\0&0&1&3\end{pmatrix}⎝⎛100−110−1−112−33⎠⎞。第9行调用上述程序定义的simplestLadder函数，对B作变换。运行程序输出

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

即回代变换后的矩阵B为(100501000013)\begin{pmatrix}1&0&0&5\\0&1&0&0\\0&0&1&3\end{pmatrix}⎝⎛100010001503⎠⎞，对应同解方程组{x1=5x2=0x3=3\begin{cases}x_1\quad\quad\quad\quad=5\\\quad\quad x_2\quad\quad=0\\ \quad\quad\quad\quad x_3=3\end{cases}⎩⎪⎨⎪⎧x1=5x2=0x3=3。这就是原方程组的解。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！

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

线性代数Python计算导引
线性代数是深度学习的数学基础之一,理论完备,方法经典.Python是当下AI系统首选开发工具,易学好用.教学中攒下两者结合的多个课题,覆盖大学理工科<线性代数>课程内容,写成博文以飨读者. ...
线性代数Python计算：向量组的最大无关组计算
对给定的mmm-维向量组 α1=(a11a21⋮am1),α2=(a12a22⋮am2),⋯,αn=(a1na2n⋮amn)\boldsymbol{\alpha}_1=\begin{pmatrix}a ...
线性代数Python计算：线性方程组的通解
对齐次线性方程组 {a11x1+a12x2+⋯+a1nxn=0a21x1+a22x2+⋯+a2nxn=0⋯⋯⋯am1x1+am2x2+⋯+amnxn=0(1)\begin{cases}a_{11}x_ ...
麻省理工大学线性代数1806(2)消元法及矩阵消元法矩阵行变换、列变换置换矩阵逆矩阵如沐春风、如饮甘露、醍醐灌顶的线性代数
目前为止发现的最适合人工智能的最简易.最深刻的线性代数课程麻省理工公开课:线性代数http://open.163.com/special/opencourse/daishu.html 只要你会数学中 ...
线性代数Python计算：向量空间坐标变换
向量空间中两组基AAA和BBB之间相互线性表示构成的矩阵为过渡阵.若两组基中之一为自然基,譬如AAA为自然基,则基BBB的各向量在自然基下的坐标即构成基AAA到基BBB的过渡阵P\boldsymbol ...
线性代数Python计算：线性变换的值域与核
设数域PPP上的向量空间PnP^nPn的线性变换TTT,在某个基下其变换矩阵A=(α1,α2,⋯,αn)\boldsymbol{A}=(\boldsymbol{\alpha}_1,\boldsymbo ...
线性代数Python计算：向量的模及向量间的夹角
numpy的dot函数计算两个向量α\boldsymbol{\alpha}α和β\boldsymbol{\beta}β的内积: dot(a,b)\text{dot(a,b)}dot(a,b) 两个参数 ...
利用python计算复合材料ABD矩阵以及压缩载荷的代码
原始理论参考美国HIL-handbook的经典层合板理论.经验证,此方法可靠,希望能帮助更多的同仁们. #This program is to caculate the ABD matrix# &qu ...
线性代数Python计算：解可逆系数矩阵线性方程组
numpy的linalg模块中solve函数解可逆系数矩阵的线性方程组.该函数的调用接口为 solve(A, b) \text{solve(A, b)} solve(A, b) 参数A表示系数矩阵 A ...

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

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

最新文章

热门文章