李保滨矩阵分析大作业2022：LU、QR、URV分解、Householder、Givens变换的程序实现
2024-04-13 01:40:36
介绍
本文为2022年秋季学期国科大李保滨老师的矩阵分析与应用课程大作业实现，编程语言使用python
具体作业要求：
完成课堂上讲的关于矩阵分解的LU、QR（Gram-Schmidt）、正交规约(Householder reduction 和Givens reduction)和 URV程序实现，要求如下：
1、一个综合程序，根据选择参数的不同，实现不同的矩阵分解；在此基础上，实现Ax=b方程组的求解，以及计算A的行列式；
2、可以用matlab、Python等编写程序，需附上简单的程序说明，比如参数代表什么意思，输入什么，输出什么等等，附上相应的例子；
注意：本文因时间仓促，中间实现难免存在疏漏与错误，还请劳烦大家提出宝贵建议！
文章目录
- 介绍
- 本文实现的功能
- 一、LU分解
- 1.存在LU分解的条件
- 2.LU分解结果
- 3.解方程的方法
- 4.程序代码实现
- 5.使用示例
- 二、QR分解（施密特正交化）
- 1.存在QR分解的条件
- 2.QR分解结果
- 3.解方程的方法
- 4.程序代码实现
- 5.使用示例
- 三、HouseHolder规约（反射算子）
- 1.存在HouseHolder规约的条件
- 2.HouseHolder规约结果
- 3.解方程的方法
- 4.程序代码实现
- 5.使用示例
- 四、Givens规约（旋转算子）
- 1.存在Givens规约的条件
- 2.Givens规约结果
- 3.解方程的方法
- 4.程序代码实现
- 5.使用示例
- 五、URV分解
- 1.存在URV分解的条件
- 2.URV分解结果
- 3.解方程的方法
- 4.程序代码实现
- 5.使用示例
- 六. 主文件（满足作业要求）
- 1.实现功能
- 2.运行步骤
- 3.程序代码实现
- 补充说明
- 1.对非法输入报警
- 2.行列式不存在提示
- 总结
本文实现的功能
- 可以自己选择实现五种中的任意一种分解
- 除LU分解外，其他分解均支持非方阵的分解，以及相应的方程和行列式求解
- 可以自己选择是否要解方程
- 解方程分为LU、QR、URV三种不同的解法，行列式的求解均是使用对应的分解求出的
补充说明：
- 对于行列式不存在的A矩阵得到的行列式为“nonexistence”
- 两种正交规约得到的是转换后的Q和R矩阵而非P和T矩阵
- 如果遇到了不符合要求的矩阵（比如LU中不是满秩），则程序终止，并报出错误原因
- 对于方程求解，如果有解且唯一，则列出解；如果有解且不唯一（无穷解），则随机列出一个解；如果无解且存在最小二乘解，则列出最小二乘解；如果无解且计算不出最小二乘解，则只显示无解（no solution）
- 置换矩阵P的行列式符号使用置换次数求，正交矩阵Q的行列式符号利用其LU分解求
一、LU分解
1.存在LU分解的条件
矩阵是可逆的方阵
2.LU分解结果
将矩阵分解为下三角矩阵和上三角矩阵的乘积
PA=LU，其中L是下三角矩阵，U是上三角矩阵，A是初等行变换矩阵
3.解方程的方法
求解方程Ax=b：
1.首先两边同乘以P得到PA=Pb，也就是LUx=Pb
2.令Ly=Pb, Ux=y
3.首先利用回代法求解Ly=Pb
4.再利用回代法求解Ux=y得到解x
4.程序代码实现
python文件为matrix_LU.py，主要函数包含两个：
1.LU_Factor(matrix)函数：输入矩阵A，输出LU分解的结果P、L、U矩阵
2.solve_LU_eqt(P,L,U,b)函数：输入LU分解后的矩阵P、L、U与等式右边向量b得到方程Ax=b的解x
matrix_LU.py文件代码如下，具体的实现方法与步骤写在注释中：
import numpy as np
from numpy.linalg import matrix_rank#控制小数点精度
np.set_printoptions( precision=3, suppress=True )#交换矩阵两行
def exchage_2_row(matrix,i,j):t=np.copy(matrix)matrix[i]=matrix[j]matrix[j]=t[i]return matrix#一般的LU分解，PA=LU
def LU_Factor(matrix):#矩阵运算都用numpy来处理matrix=np.array(matrix)(len_r,len_c)=matrix.shape#首先保证是方阵才能LU分解if len_r!=len_c:print('-------------- warning -----------------')print('warning: LU factor only accept square matrix!!!')return -1#不可逆的话也不可以分解if matrix_rank(matrix)!=len_r:print('-------------- warning -----------------')print('warning: LU factor only accept nonsingular matrix!!!')return -1 ## 第一步：增加b列向量得到增广矩阵#初始化行交换顺序矩阵ex_r_list=np.array([i+1 for i in range(len_r)])extra_matirx=np.concatenate((matrix,np.expand_dims(ex_r_list,axis=0).T),axis=1).astype(np.float32)## 第二步：利用第1类和第3类行变换将矩阵上三角化L_matrix=np.zeros(matrix.shape)idx_c=0for idx_r in range(len_r-1):#最后一次就不用操作了#部分主元法，找一个绝对值最大元素的一行,并将其交换到当前的最顶端choice_list=np.abs(extra_matirx[idx_r:len_r,idx_c])pivot=idx_r+np.argmax(choice_list)#如果最大主元都是0的话，那就肯定不可逆if extra_matirx[pivot,idx_c]==0:print('-------------- warning -----------------')print('warning: LU factor only accept nonsingular matrix!!!')return -1 #可以进行下去的话，就交换主元使其到最前面extra_matirx=exchage_2_row(extra_matirx,idx_r,pivot)L_matrix=exchage_2_row(L_matrix,idx_r,pivot)#高斯消去for i in range(idx_r+1,len_r):if extra_matirx[i,idx_c]!=0 :#得到第3类变换的系数L_matrix[i,idx_c]=alpha=float(extra_matirx[i,idx_c])/extra_matirx[idx_r,idx_c]extra_matirx[i,:len_r] = extra_matirx[i,:len_r]- extra_matirx[idx_r,:len_r]*alphaidx_c+=1## 第三步：检查是否可逆（U矩阵对角线是否全非0）for idx in range(len_r):if extra_matirx[idx,idx]==0:return -1## 第四步：得到分解结果U_matrix=extra_matirx[:,:len_r]P_matrix=np.zeros(matrix.shape)for idx in range(len_r):#L矩阵对角全化为1L_matrix[idx,idx]=1#算出P矩阵P_matrix[idx,int(extra_matirx[idx,-1])-1]=1return P_matrix,L_matrix,U_matrixdef solve_LU_eqt(P,L,U,b):dim=len(b)#先得到PbPb=np.dot(P,b)## 令Ly=Pb, Ux=y#先求解Ly=Pb，利用回代法求解y=np.copy(Pb)for idx in range(dim):#从前往后y[idx]=Pb[idx]#得到等式和for i in range(idx):y[idx] -= L[idx][i]*y[i]#减去前面的分量#同样是回代法求解Ux=y，不过是上三角形式x=np.copy(y)for idx in range(dim-1,-1,-1):#从后往前x[idx]=y[idx]#得到等式和for i in range(idx+1,dim):x[idx] -= U[idx][i]*x[i]#减去前面的分量x[idx] = float(x[idx])/U[idx][idx]#注意U的对角不是1return x
5.使用示例
测试代码如下：
a=np.array([[1,2,-3,4],[4,8,12,-8],[2,3,2,1],[-3,-1,1,-4]])b=np.array([3,60,1,5])x,y,z=LU_Factor(a)print("P:",x)print("L:",y)print("U:",z)print("x:",solve_LU_eqt(x,y,z,b))
终端输出为：
P: [[0. 1. 0. 0.][0. 0. 0. 1.][1. 0. 0. 0.][0. 0. 1. 0.]]
L: [[ 1. 0. 0. 0. ][-0.75 1. 0. 0. ][ 0.25 0. 1. 0. ][ 0.5 -0.2 0.333 1. ]]
U: [[ 4. 8. 12. -8.][ 0. 5. 10. -10.][ 0. 0. -6. 6.][ 0. 0. 0. 1.]]
x:[ 12. 6. -13. -15.]
二、QR分解（施密特正交化）
1.存在QR分解的条件
矩阵的秩等于列数，可以是非方阵
2.QR分解结果
将矩阵分解为正交矩阵和上三角矩阵的乘积
A=QR，其中Q是正交矩阵（行大于等于列），R是上三角矩阵
3.解方程的方法
求解方程Ax=b：
注意：如果方程本身无解则求出来的可能是最小二乘解
1.首先利用QR分解得到A=QR，即QRx=b
2.将QTQ^TQT当作Q−1Q^-1Q−1乘到等式右边Rx=QTQ^TQTb
3.去掉R和b多余的行（因为A不一定是方阵）
4.利用回代法求解Rx=Q.Tb
5.若4中方程有解，则讨论是唯一解还是无穷解
6.若4中方程无解则返回无解
4.程序代码实现
python文件为matrix_Smit.py，主要函数包含两个：
1.QR_Factor(matrix)函数：输入矩阵A，输出QR分解的结果Q、R矩阵
2.solve_QR_eqt(Q,R,b)函数：输入QR分解后的矩阵Q、R与等式右边向量b得到方程Ax=b的解x
matrix_Smit.py文件代码如下，具体的实现方法与步骤写在注释中：
import numpy as np
import random
from numpy.linalg import matrix_rank#控制小数点精度
np.set_printoptions( precision=3, suppress=True )#classic施密特正交化实现的QR分解
#可以应用于非方阵！！！
#其中Q是正交矩阵（行大于等于列），R是上三角矩阵
def QR_Factor(matrix):#矩阵运算都用numpy来处理matrix=np.array(matrix).astype(np.float32)(len_r,len_c)=matrix.shape#首先保证矩阵的秩等于列数（因为行大于等于列）rk=matrix_rank(matrix)if rk!=len_c:print('-------------- warning -----------------')print('warning: Smit QR factor only accept column full rank matrix!!!')return -1#初始化Q，R矩阵Q_matrix=matrix.copy()R_matrix=np.zeros(matrix.shape)ok_cnt=0#表示已经正交化过的列数#对于matrix每一列做正交化for x in matrix.T:#对于每一列u=np.copy(x)for idx in range(ok_cnt):#计算内积并填入R矩阵中inner_product=np.dot(Q_matrix[:,idx],x)R_matrix[idx,ok_cnt]=inner_product#减去其他维度的分量以实现正交u -= inner_product*Q_matrix[:,idx]norm=np.linalg.norm(u)#模长R_matrix[ok_cnt,ok_cnt]=norm#模长填入对角线Q_matrix[:,ok_cnt]=u/norm#归一化得到标准正交向量ok_cnt+=1 #去掉R多余的行，保证是方阵if len_r!=len_c:R_matrix=np.delete(R_matrix,[len_c,len_r-1],axis=0)return Q_matrix,R_matrix#判断方程Ax=b是否有解
def is_solvable(A,b):extra_matirx=np.concatenate((A,np.expand_dims(b,axis=0).T),axis=1)#如果增广矩阵的秩等于A的秩就有解return matrix_rank(A)==matrix_rank(extra_matirx)#输入的Q可能不是方阵
#结果可能是最小二乘解！！！！！
def solve_QR_eqt(Q,R,b):dim=R.shape[1]#等价于求解 Rx=Q.TbQTb=np.dot(Q.T,b)#去掉R和b多余的行if R.shape[0]!=R.shape[1]:R=np.delete(R,[R.shape[1],R.shape[0]-1],axis=0)QTb=QTb[:R.shape[1]]#回代法求解Rx=Q.Tb,上三角形式#但是这个方程不一定有解，要分情况讨论if is_solvable(R,QTb):#1.方程有解x=np.copy(QTb)for idx in range(dim-1,-1,-1):#从后往前x[idx]=QTb[idx]#得到等式和for i in range(idx+1,dim):x[idx] -= R[idx][i]*x[i]#减去前面的分量if int(R[idx][idx])==0:#2.如果R对角线存在0，那就是有无穷解#由于是无穷解，所以随便给出一个解即可x[idx] = random.randrange(-100,100)/10.0else:#3.唯一解x[idx] = float(x[idx])/R[idx][idx]#注意U的对角不是1else:#方程无解return -1return x
5.使用示例
测试代码如下：
matrix1=np.array([[1,0,-1],[1,2,1],[1,1,-3],[0,1,1]])b=np.array([1,1,1,1])#print(QR_Factor(matrix1))x,y=QR_Factor(matrix1)print("Q:",x)print("R:",y)print("x:",solve_QR_eqt(x,y,b))
终端输出为：
Q: [[ 0.577 -0.577 0.408][ 0.577 0.577 0.408][ 0.577 0. -0.816][ 0. 0.577 0. ]]
R: [[ 1.732 1.732 -1.732][ 0. 1.732 1.732][ 0. 0. 2.449]]
x: [0.667 0.333 0. ]
三、HouseHolder规约（反射算子）
1.存在HouseHolder规约的条件
没有限制，可以是非方阵
2.HouseHolder规约结果
将矩阵A规约为PA=T，其中P是正交矩阵(方阵)，T是伪上三角矩阵（可以不是方阵）
为了整个作业的统一性，将其变换为QR分解的形式：A=QR。其中Q=PTP^TPT，R=T
3.解方程的方法
由于其可以化为QR分解的形式，所以求解方程的方法同前面所述QR分解中的解方程方法
4.程序代码实现
python文件为matrix_Householder.py，主要函数为：
1.Householder_Reduction(matrix)函数：输入矩阵A，输出QR分解的结果Q、R矩阵
matrix_Householder.py文件代码如下，具体的实现方法与步骤写在注释中：
import numpy as np
import math#控制小数点精度
np.set_printoptions( precision=3, suppress=True )#对称矩阵实现的正交规约 PA=T
#可以应用于非方阵！！！
#其中P是正交矩阵(方阵)，T是伪上三角矩阵（可以不是方阵）
def Householder_Reduction(matrix):#矩阵运算都用numpy来处理matrix=np.array(matrix).astype(np.float32)(len_r,len_c)=matrix.shapeP_matrix=np.identity(len_r)#分别处理两种不是方阵的情况if len_r > len_c:iter_num=len_celse:iter_num=len_r-1for idx in range(iter_num):#首先构造ei向量e = np.zeros(len_r-idx)e[0] = 1I = np.identity(len_r-idx)#注意u要是列向量u = matrix[idx:,idx].T - math.sqrt(np.dot(matrix[idx:,idx].T,matrix[idx:,idx]))*e.Texp_u = np.expand_dims(u,axis=0)#得到反射算子Rflt = I - (2.0/np.dot(u.T,u))*np.matmul(exp_u.T,exp_u)#拓展成正常大小exp_Rflt = np.identity(len_r)exp_Rflt[idx:,idx:] = Rflt#更新P_matrix = np.dot(exp_Rflt,P_matrix)matrix = np.dot(exp_Rflt,matrix)Q_matrix=P_matrix.TR_matrix=matrixreturn Q_matrix,R_matrix
5.使用示例
测试代码如下：
matrix=np.array([[4,-3,4],[2,-14,-3],[-2,14,0],[1,-7,15]])Q,R=Householder_Reduction(matrix)print('Q:',Q)print('R:',R)
终端输出为：
Q: [[ 0.8 0.6 -0. -0. ][ 0.4 -0.533 -0.333 -0.667][-0.4 0.533 0.133 -0.733][ 0.2 -0.267 0.933 -0.133]]
R: [[ 5. -15. 5.][ -0. 15. -0.][ 0. -0. 15.][ -0. -0. -0.]]
四、Givens规约（旋转算子）
1.存在Givens规约的条件
没有限制，可以是非方阵
2.Givens规约结果
将矩阵A规约为PA=T，其中P是正交矩阵(方阵)，T是伪上三角矩阵（可以不是方阵）
为了整个作业的统一性，将其变换为QR分解的形式：A=QR。其中Q=PTP^TPT，R=T
3.解方程的方法
由于其可以化为QR分解的形式，所以求解方程的方法同前面所述QR分解中的解方程方法
4.程序代码实现
python文件为matrix_Givens.py，主要函数为：
1.Givens_Reduction(matrix)函数：输入矩阵A，输出QR分解的结果Q、R矩阵
matrix_Givens.py文件代码如下，具体的实现方法与步骤写在注释中：
import numpy as np#控制小数点精度
np.set_printoptions( precision=3, suppress=True )#旋转矩阵实现的正交规约 PA=T
#可以应用于非方阵！！！
#其中P是正交矩阵(方阵)，T是伪上三角矩阵（可以不是方阵）
def Givens_Reduction(matrix):#矩阵运算都用numpy来处理matrix=np.array(matrix)(len_r,len_c)=matrix.shape#首先初始化P，T矩阵P_matrix=np.identity(len_r)T_matrix=np.copy(matrix)#按顺序每一列进行变换for idx_c in range(len_c):for idx_r in range(len_r-1,idx_c,-1):#首先得到要变换的两个数x=T_matrix[idx_c,idx_c]y=T_matrix[idx_r,idx_c]norm=np.sqrt(x**2+y**2)#直接构造旋转矩阵c=x/norms=y/normP_idx=np.eye(len_r)#先生成对角矩阵P_idx[idx_c,idx_c]=P_idx[idx_r,idx_r]=cP_idx[idx_c,idx_r]=sP_idx[idx_r,idx_c]=-s#累乘可以得到最终的P和T矩阵P_matrix=np.matmul(P_idx,P_matrix)T_matrix=np.matmul(P_idx,T_matrix)#由于要求是求矩阵分解，所以转化为QR分解的形式Q_matrix=P_matrix.T#由于P是正交方阵，所以转置=逆return Q_matrix,T_matrix
5.使用示例
测试代码如下：
matrix1=np.array([[1,0,-1],[1,2,1],[1,1,-3],[0,1,1]])Q,R=Givens_Reduction(matrix1)print('Q:',Q)print('R:',R)
终端输出为：
Q: [[ 0.577 -0.577 0.408 -0.408][ 0.577 0.577 0.408 0.408][ 0.577 0. -0.816 -0. ][ 0. 0.577 0. -0.816]]
R: [[ 1.732 1.732 -1.732][ 0. 1.732 1.732][ 0. 0. 2.449][ 0. -0. -0. ]]
五、URV分解
1.存在URV分解的条件
有限制，可以是非方阵
2.URV分解结果
对于m∗nm*nm∗n矩阵A，可以分解为A=URVTV^TVT，其中U是正交矩阵(m∗mm*mm∗m)，V是正交矩阵(n∗nn*nn∗n)，R为m∗nm*nm∗n的矩阵
3.解方程的方法
求解方程Ax=b：
注意：如果方程本身无解则求出来的可能是最小二乘解
1.首先URV分解得到A=URVTV^TVT
2.令y=VTV^TVTx,先求解URy=b,也就是Ry=UTU^TUTb
3.去掉R和UTU^TUTb多余的行
4.若Ry=UTU^TUTb有解，则讨论是唯一解还是无穷解
5.若Ry=UTU^TUTb无解则返回无解
6.如果有解，则利用x=Vy求出最后的x
4.程序代码实现
python文件为matrix_URV.py，主要函数包含两个：
1.URV_Factor(matrix)函数：输入矩阵A，输出URV分解的结果U、R、V矩阵
2.solve_URV_eqt(P,L,U,b)函数：输入URV分解后的矩阵U、R、V与等式右边向量b得到方程Ax=b的解x
matrix_URV.py文件代码如下，具体的实现方法与步骤写在注释中：
注意：由于本代码需要使用householder的方法，所以需要将matrix_Householder.py文件放在和本文件同目录下
import numpy as np
from numpy.linalg import matrix_rank
import random
from matrix_Householder import Householder_Reduction#控制小数点精度
np.set_printoptions( precision=3, suppress=True )#实现URV分解,A=URVT
#其中U是正交矩阵(m*m)，V是正交矩阵(n*n)R为m*n的矩阵
def URV_Factor(matrix):#矩阵运算都用numpy来处理matrix=np.array(matrix).astype(np.float32)## 利用两次householder归约得到URV分解#先做一次Q_mtx_1,R_mtx_1=Householder_Reduction(matrix)rank_Q1=matrix_rank(R_mtx_1)tmp_mtx_1=R_mtx_1[:rank_Q1,:]#换个方向再做一次Q_mtx_2,R_mtx_2=Householder_Reduction(tmp_mtx_1.T)rank_Q2=matrix_rank(R_mtx_2)tmp_mtx_2=R_mtx_2[:rank_Q2,:]#得到URV矩阵U_matrix=Q_mtx_1V_matrix=Q_mtx_2R_matrix=np.zeros(matrix.shape).astype(np.float32)R_matrix[:rank_Q2,:rank_Q2]=tmp_mtx_2.Treturn U_matrix,R_matrix,V_matrix#判断方程Ax=b是否有解
def is_solvable(A,b):extra_matirx=np.concatenate((A,np.expand_dims(b,axis=0).T),axis=1)#如果增广矩阵的秩等于A的秩就有解return matrix_rank(A)==matrix_rank(extra_matirx)#求解Ax=b，URVTx=b
#结果可能是最小二乘解！！！！！
def solve_URV_eqt(U,R,V,b):dim=R.shape[1]#令y=VTx,先求解URy=b,也就是Ry=UTbUTb=np.dot(U.T,b)#去掉R和UTb多余的行if R.shape[0]!=R.shape[1]:R=np.delete(R,[R.shape[1],R.shape[0]-1],axis=0)UTb=UTb[:R.shape[1]]#但是这个方程不一定有解，要分情况讨论if is_solvable(R,UTb):#1.方程有解y=np.copy(UTb)for idx in range(dim):#从前往后，R是下三角y[idx]=UTb[idx]#得到等式和for i in range(idx):y[idx] -= R[idx][i]*y[i]#减去前面的分量if int(R[idx][idx])==0:#2.如果R对角线存在0，那就是有无穷解#由于是无穷解，所以随便给出一个解即可y[idx] = random.randrange(-100,100)/10.0else:#3.唯一解y[idx] = float(y[idx])/R[idx][idx]#注意U的对角不是1else:#方程无解return -1#再利用x=Vy求出最后的xx=np.dot(V,y)return x
5.使用示例
测试代码如下：
matrix=np.array([[4,-3,4],[2,-14,-3],[-2,14,0]])b=np.array([5,-15,0])u,r,v=URV_Factor(matrix)print("U:",u)print("R:",r)print("V:",v)print("x:",solve_URV_eqt(u,r,v,b))
终端输出为：
U: [[ 0.816 0.577 0. ][ 0.408 -0.577 -0.707][-0.408 0.577 -0.707]]
R: [[ 14.86 -0. 0. ][-12.927 7.587 -0. ][ 0.291 1.626 1.33 ]]
V: [[ 0.33 0.562 -0.759][-0.934 0.311 -0.176][ 0.137 0.767 0.627]]
x: [-4.2 -0.6 5. ]
六. 主文件（满足作业要求）
1.实现功能
1.可以自行通过数字来选择矩阵分解的方法
2.可以自行选择是否需要解方程,并且可以给出最小二乘解，也可以提示方程无解
3.可以求解矩阵的行列式
2.运行步骤
主文件为main_equation.py
1.直接运行主函数main_equation.py
2.输入1-5之间的数来选择矩阵分解的类型
Please choose factoration method:
1.LU 2.QR(Smit) 3.Householder 4.Givens 5.URV
my factor choice(1-5):
3.选择1或2来表示是否需要解方程Ax=b
Please choose factoration method:
1.LU 2.QR(Smit) 3.Householder 4.Givens 5.URV
my factor choice(1-5):1
do you want to solve equations? 1.yes 2.no
my solve equation choice(1or2):
4.输入A矩阵
示例1：[[1,2,3],[4,5,6],[7,8,9]]
示例2 : [[4,-3,4],[2,-14,-3],[-2,14,0],[1,-7,15]]
注意：矩阵必须按照上面一样写成一行，并且中间不能有空格
Please choose factoration method:
1.LU 2.QR(Smit) 3.Householder 4.Givens 5.URV
my factor choice(1-5):1
do you want to solve equations? 1.yes 2.no
my solve equation choice(1or2):1
Please type marix A !
matirx A:
5.（如果3中选择了1）输入方程中的b
Please choose factoration method:
1.LU 2.QR(Smit) 3.Householder 4.Givens 5.URV
my factor choice(1-5):1
do you want to solve equations? 1.yes 2.no
my solve equation choice(1or2):1
Please type marix A !
matirx A:[[1,2,-3,4],[4,8,12,-8],[2,3,2,1],[-3,-1,1,-4]]
Please type vector b !
vector b:
6.得到结果（方程的解、矩阵分解、行列式det（A））
Please choose factoration method:
1.LU 2.QR(Smit) 3.Householder 4.Givens 5.URV
my factor choice(1-5):1
do you want to solve equations? 1.yes 2.no
my solve equation choice(1or2):1
Please type marix A !
matirx A:[[1,2,-3,4],[4,8,12,-8],[2,3,2,1],[-3,-1,1,-4]]
Please type vector b !
vector b:[3,60,1,5]
-------------- output -----------------
euqation solution:
[ 12. 6. -13. -15.]
-------------- output -----------------
P:
[[0. 1. 0. 0.][0. 0. 0. 1.][1. 0. 0. 0.][0. 0. 1. 0.]]
L:
[[ 1. 0. 0. 0. ][-0.75 1. 0. 0. ][ 0.25 0. 1. 0. ][ 0.5 -0.2 0.333 1. ]]
U:
[[ 4. 8. 12. -8.][ 0. 5. 10. -10.][ 0. 0. -6. 6.][ 0. 0. 0. 1.]]
det(A): 120.0
3.程序代码实现
python文件为main_equation.py，里面包含了使用矩阵分解求解矩阵行列式的方法det_Q函数，具体实现方法在代码中有注释说明。
在运行本文件前需要将前面五个文件都放在同一目录下，即matrix_LU.py、matrix_Smit.py、matrix_Givens.py、matrix_Householder.py、matrix_URV.py。
main_equation.py文件代码如下，具体的实现方法与步骤写在注释中：
import numpy as np
import sys
from numpy.linalg import matrix_rank
from matrix_LU import LU_Factor,solve_LU_eqt
from matrix_Smit import QR_Factor,solve_QR_eqt
from matrix_Givens import Givens_Reduction
from matrix_Householder import Householder_Reduction
from matrix_URV import URV_Factor,solve_URV_eqtdef load_data(str,type='m'):if type=='m':#返回矩阵matrix=[]v_group=str[2:-2].split('],[')#得到行向量组for grp in v_group:#解码行向量nums_list=grp.split(',')row_vector=[]for item in nums_list:row_vector.append(float(item))matrix.append(row_vector)#行向量一行一行输入矩阵return matrixelif type=='v':#返回向量nums_list=str[1:-1].split(',')vector=[]for item in nums_list:vector.append(float(item))return vectordef factor_func(type_choice,matrix):if type_choice==1:return LU_Factor(matrix)elif type_choice==2:return QR_Factor(matrix)elif type_choice==3:return Householder_Reduction(matrix)elif type_choice==4:return Givens_Reduction(matrix)elif type_choice==5:return URV_Factor(matrix)#求对角线的乘积
def diag_mul(matrix):dim=matrix.shape[1]mul=1for idx in range(dim):mul*=matrix[idx][idx]return mul#判断方程Ax=b是否有解
def is_solvable(A,b):extra_matirx=np.concatenate((A,np.expand_dims(b,axis=0).T),axis=1)#如果增广矩阵的秩等于A的秩就有解return matrix_rank(A)==matrix_rank(extra_matirx)#确定置换矩阵P的行列式（符号）
def det_P(matrix):len_=matrix.shape[0]order_list=[]for i in range(len_):for j in range(len_):if matrix[i][j]==1:order_list.append(j+1)#重建顺序序列#重新升序排列，记录交换次序ex_cnt=0#交换次数for i in range(len(order_list)-1):#冒泡排序for j in range(len(order_list)-i-1):if order_list[j] > order_list[j+1]:t = order_list[j] order_list[j] = order_list[j+1]order_list[j+1] = tex_cnt+=1return (-1)**ex_cnt#求正交矩阵行列式
def det_Q(matrix):#首先对矩阵进行PLU分解P,L,U=LU_Factor(matrix)return det_P(P)*np.sign(diag_mul(U))def print_factor(choice,mtx_l,mtx_r,mtx_m=0):if choice==1:print('P:')print(mtx_l)print('L:')print(mtx_m)print('U:')print(mtx_r)#U行列式就是上三角矩阵U的对角元素相乘,再乘上P行列式det_A=det_P(mtx_l)*diag_mul(mtx_r)print('det(A):',round(det_A,4))elif choice==5:print('U:')print(mtx_l)print('R:')print(mtx_m)print('V:')print(mtx_r)if mtx_m.shape[0]==mtx_m.shape[1]:#如果是方阵,注意要乘上两个正交矩阵的行列式det_A=det_Q(mtx_l)*diag_mul(mtx_m)*det_Q(mtx_r)print('det(A):',round(diag_mul(mtx_m),4))else:print('det(A):nonexistence')print('warning:not square matrix do not have det(A)!!!')else:print('Q:')print(mtx_l)print('R:')print(mtx_r)#行列式就是上三角矩阵R的对角元素相乘，再乘上Q的行列式if mtx_l.shape[0]==mtx_l.shape[1] and mtx_r.shape[0]==mtx_r.shape[1]:det_A=det_Q(mtx_l)*diag_mul(mtx_r)print('det(A):',round(det_A,4))else:print('det(A):nonexistence')print('warning:not square matrix do not have det(A)!!!')#最小二乘解也认为是解
def all_sovle_eqt(choice,mtx_l,mtx_r,b,mtx_m=0):if choice==1:return solve_LU_eqt(mtx_l,mtx_m,mtx_r,b)elif choice==5:return solve_URV_eqt(mtx_l,mtx_m,mtx_r,b)else:return solve_QR_eqt(mtx_l,mtx_r,b)if __name__=="__main__":#首先获取分解类型print('Please choose factoration method:\n''1.LU 2.QR(Smit) 3.Householder 4.Givens 5.URV')choice=int(input('my factor choice(1-5):'))#选择是否求解方程print('do you want to solve equations? 1.yes 2.no')solve_eqt=int(input('my solve equation choice(1or2):'))#然后得到矩阵输入print('Please type marix A !')#输入示例：'[[1,2,3],[4,5,6],[7,8,9]]',不能有空格A_string=input('matirx A:')A=load_data(A_string,type='m')#解码得到A#得到分解结果if choice==1 or choice==5:#PLU URVif factor_func(choice,A)==-1:#出错的话就终止程序sys.exit()#错误信息在函数里已打印mtx_l,mtx_m,mtx_r=factor_func(choice,A)else: #QRif factor_func(choice,A)==-1:#出错的话就终止程序sys.exit()mtx_l,mtx_r=factor_func(choice,A)mtx_m=0if solve_eqt==1:print('Please type vector b !')#输入示例：'[1,2,3]'B_string=input('vector b:')b=load_data(B_string,type='v')#解码得到b#求解方程 x=all_sovle_eqt(choice,mtx_l,mtx_r,b,mtx_m)print('-------------- output -----------------')if is_solvable(A,b):print('euqation solution:')print(x)else:print('no solution')if x.all()!=-1:#如果最小二乘解可以解出来print('but the least square solution is')print(x)print('-------------- output -----------------')#输出分解以及行列式结果if choice==1 or choice==5:print_factor(choice,mtx_l,mtx_r,mtx_m)else:print_factor(choice,mtx_l,mtx_r)补充说明
1.对非法输入报警
下面以LU分解为例，说明输入的A不符合要求的时候程序会报警并终止
输入1: 不可逆矩阵
Please choose factoration method:
1.LU 2.QR(Smit) 3.Householder 4.Givens 5.URV
my factor choice(1-5):1
do you want to solve equations? 1.yes 2.no
my solve equation choice(1or2):2
Please type marix A !
matirx A:[[1,2,3],[2,4,6],[6,3,5]]
-------------- warning -----------------
warning: LU factor only accept nonsingular matrix!!!
输入2: 非方阵
Please choose factoration method:
1.LU 2.QR(Smit) 3.Householder 4.Givens 5.URV
my factor choice(1-5):1
do you want to solve equations? 1.yes 2.no
my solve equation choice(1or2):2
Please type marix A !
matirx A:[[1,2,3],[4,5,6]]
-------------- warning -----------------
warning: LU factor only accept square matrix!!!
2.行列式不存在提示
如果输入矩阵不是方阵则会提示行列式不存在
Please choose factoration method:
1.LU 2.QR(Smit) 3.Householder 4.Givens 5.URV
my factor choice(1-5):5
do you want to solve equations? 1.yes 2.no
my solve equation choice(1or2):1
Please type marix A !
matirx A:[[4,-3,4],[2,-14,-3],[-2,14,0],[1,-7,15]]
Please type vector b !
vector b:[5,-15,0,30]
-------------- output -----------------
no solution
but the least square solution is
[-0.8 0.2 2.2]
-------------- output -----------------
U:
[[ 0.8 0.6 -0. -0. ][ 0.4 -0.533 -0.333 -0.667][-0.4 0.533 0.133 -0.733][ 0.2 -0.267 0.933 -0.133]]
R:
[[ 16.583 -0. 0. ][-13.568 6.396 -0. ][ 4.523 9.594 10.607][ 0. 0. 0. ]]
V:
[[ 0.302 0.64 -0.707][-0.905 0.426 -0. ][ 0.302 0.64 0.707]]
det(A):nonexistence
warning:not square matrix do not have det(A)!!!
总结
1.本大作业属于对课上所学的应用实践，可以较好地锻炼同学们的编程能力
2.本作业中针对有无穷解的情况只随机给出了其中的一个解，如果进一步做的更好可以比较完整地表示出无穷解的具体情况
3.如果本文中出现某些问题或者错误大家可以联系我一起讨论，也可以在文下评论。
李保滨矩阵分析大作业2022：LU、QR、URV分解、Householder、Givens变换的程序实现相关推荐
- 中国科学院 导师推荐 计算机,中国科学院计算机与控制学院硕士生导师：李保滨副教授...
李保滨 男 副教授 电子邮件:libb@ucas.ac.cn 通信地址:北京市中关村东路80号青年公寓6号楼206 邮政编码: 100190 研究领域 滤波器设计及其在图像处理中的应用 向量小波的设计 ...
- 哈工大CSAPP大作业 2022
计算机系统 大作业 题 目 程序人生-Hello's P2P 专 业 计算学部 学 号 120L021920 班 级 2003006 学 生 李启明 指 导 教 师 吴锐 计算机科学与技术学院 202 ...
- 程序人生-哈工大计算机系统大作业2022春
计算机系统 大作业 题 目 程序人生-Hello's P2P 专 业 计算学部 学 号 7203610401 班 级 2036012 学 生 王元辰 指 导 ...
- 哈尔滨工业大学计算机系统大作业2022春
计算机系统 大作业 题 目 程序人生-Hello's P2P 专 业 计算机 学 号 120L02**** 班 级 2003005 学 生 无敌飞龙 指 导 教 ...
- 哈尔滨工业大学CSAPP大作业-2022春
计算机系统 大作业 题 目 程序人生-Hello's P2P 专 业 人工智能(未来技术) 学 号 7203610202 班 级 2036015 学 生 熊 ...
- 哈工大计算机系统大作业2022春
计算机系统 大作业 题 目 程序人生-Hello's P2P 专 业 计算机科学与技术学院 学 号 7203610110 班 级 2036011 学 生 王 ...
- HIT计算机系统大作业2022
计算机系统 大作业 计算机科学与技术学院 2021年5月 摘 要 学习计算机系统后,从底层对一个简单的Hello.c如何一步步编译,链接,加载,执行等进行浅层的描述. 关键词:计算机系统,底层. 目 ...
- csapp程序人生大作业
正在上传-重新上传取消 计算机系统 大作业 题 目 程序人生-Hello's P2P 专 业 计算学部 学 号 120L020311 班 级 2003003 学 ...
- HIT计算机系统CSAPP大作业
HIT计算机系统CSAPP大作业 摘 要 一.第1章 概述 1.1 Hello简介 ·P2P过程 ·020过程 1.2 环境与工具 1.2.1 硬件环境 1.2.2 软件环境 1.2.3 开发工具 1 ...
最新文章
- ICML 2018 | 清华排名国内居首：大会论文接收情况一览
- mysql事务在提交后才发送给数据库执行_从一个线上问题分析binlog与内部XA事务提交过程...
- STM32，CAN总线过滤器的设置详细讲解
- 谷歌40人发表59页长文：为何真实场景中ML模型表现不好？
- 《嵌入式系统可靠性设计技术及案例解析》读书笔记（五）
- LiveVideoStackCon 2022 上海站 专题抢先看（3）
- Spring和JSF集成：MVC螺母和螺栓
- 小程序成长之路（四）-- 深入腾讯云（环境搭建）
- 【例题+习题】【数值计算方法复习】【湘潭大学】（二）
- 考研数学一基础技巧题汇总
- Guava Cache 使用笔记
- 从拉新、促活/留存和营收说起，做运营到底是在做什么？
- 程序设计的最基本的三种结构
- 奉子成婚，永远不可能成为潮流
- 1.1 区块链专业术语（中英对照）
- 在线流程图和思维导图开发技术详解（六）
- [ArcGIS] 空间分析（八） 水文分析
- 腾讯 美团 百度 网易游戏 2015校园招聘南京笔试面试之总结分析
- VMware NAT 模式配置端口映射
- java class dex_class文件与dex文件分析
热门文章
- doubango编译过程中遇到的:Cannot open include file: 'com/sun/star/beans/XPropertySet.hpp
- 舆情监测意思及监测工作流程详介
- 习题:有100个学生种100棵树，其中高中生每人种3棵树，初中生每人种1棵树，小学生每3人种1棵树，问高中生、初中生、小学生各有多少人？
- with respect to是什么意思？
- 不是那个层次的人，也便没机会领略…
- 神经网络初认识——BP神经网络(7月18，再次认识）
- 无法打开物理文件 XXX.mdf,操作系统错误 5:5(拒绝访问。)的解决办法
- Mac上鼠标滚轮方向是和Win相反的，系统中设置后触摸板的方向又跟着变了
- java在spring mvc中的图片接收与发送处理
- G711 编码解码及丢包隐藏处理(PLC)