Python解决矩阵的PLU分解及求矩阵的逆

Python解决矩阵的PLU分解及求矩阵的逆

关于PLU的分解基础知识就不叙述了,可以自己去看矩阵分析的书,大体上和高斯消去法差不多。

PLU分解被经常用在 A x = b Ax=b Ax=b的求解上
在这里 x x x是一个 n n n维的列向量,因为此时L,U分别是下三角和上三角矩阵,所以对求解带来了极大的便利性

细心的朋友可能会发散思维, A x = b Ax=b Ax=b是不是可以用于求矩阵A的逆,确实如此

A x = b P A x = P b P A = L U L U x = P b Ax=b\\ PAx = Pb\\ PA=LU \quad LUx=Pb Ax=bPAx=PbPA=LULUx=Pb
此时另 U x = y Ux=y Ux=y
L y = P b Ly=Pb Ly=Pb
到这里,我们就可以先求 y y y,再求 x x x了;若是要求逆的话, b b b就是单位矩阵的每一列,用这个方法 n n n次求解 x x x,此时的 x x x就是 A − 1 A^{-1} A1的每一列

程序如下

# @Time    : 2021/9/24 10:49
# @Author  : AlanLiang
# @FileName: PLU_Factorization.py
# @Software: VS Code
# @Github  :https://github.com/AlanLiangC
# @Mail    : liangao21@mail.ucas.ac.cn


import numpy as np
import copy

# 定义类 功能均在类内实现
class PLU():
    def __init__(self,Mat):
        self.Mat = copy.deepcopy(Mat) #深层拷贝,防止初始矩阵内存被重复占用,下同
        self.Mat_stable = copy.deepcopy(Mat)
        self.Mat_R = Mat.shape[0]
        self.Mat_C = Mat.shape[1]
        self.Mat_P = np.array([i for i in range(self.Mat_R)])
        self.Mat_P_result = np.zeros([Mat.shape[0],Mat.shape[1]],dtype = np.float32)
        self.Mat_L = np.eye(Mat.shape[0],dtype = np.float32)
        self.Mat_U = np.zeros([Mat.shape[0],Mat.shape[1]],dtype = np.float32)
        self.Mat_inversion = np.zeros([Mat.shape[0],Mat.shape[1]],dtype = np.float32)
    
    # 更新矩阵,始终将该次计算的最大主元移到最上方
    def update_mat(self,row):
        Mat_this = abs(self.Mat)
        max_num = np.max(Mat_this[row:,row])
        max_index = np.where(Mat_this[row:,row] == max_num)[0] + row
        if row == max_index :
            return self.Mat
        else:
            change = self.Mat[max_index,:]
            self.Mat[max_index,:] = self.Mat[row,:]
            self.Mat[row,:] = change
            change_P = self.Mat_P[max_index]
            self.Mat_P[max_index] = self.Mat_P[row]
            self.Mat_P[row] = change_P
        return self.Mat


    # 分解矩阵,当detail=Trur时输出分解过程
    def Factorization(self,detail = True):
        for i in range(self.Mat_R - 1):
            self.Mat = self.update_mat(i)
            if detail:
                print("第{}次更新".format(i),"\n",self.Mat,"\n",self.Mat_P)
            for j in range(self.Mat_R - i - 1):
                if self.Mat[i+j+1,i] == 0:
                    self.Mat = self.Mat
                else:
                    ratio = self.Mat[i+j+1,i] / self.Mat[i,i]
                    self.Mat[i+j+1,i+1:] = self.Mat[i+j+1,i+1:] - ratio*self.Mat[i,i+1:]
                    self.Mat[i+j+1,i] = ratio
            if detail:
                print("第{}次计算".format(i),"\n",self.Mat)
            
    # 输出分解后的P,L,U矩阵
    def get_result(self):
        for i in range(self.Mat_R):
            self.Mat_P_result[i,self.Mat_P[i]] = 1
            for j in range(self.Mat_C):
                if i > j :
                    self.Mat_L[i,j] = self.Mat[i,j]
                else:
                    self.Mat_U[i,j] = self.Mat[i,j]
        print("P:","\n",self.Mat_P_result)
        print("L:","\n",self.Mat_L)
        print("U:","\n",self.Mat_U)
        return self.Mat_P_result,self.Mat_L,self.Mat_U



    # 进行矩阵求逆,当detail=True时输出运行过程,以便进行bug定位
    # 求逆过程为逐行求逆 类比于Ax=b,x为A逆的每一列,b为单位阵I的每一列
    def get_inversion(self,detail = True):
        B = np.eye(self.Mat_R)
        for i in range(self.Mat_R):
            sub_B = np.dot(self.Mat_P_result,B[:,i])
            y = np.zeros([self.Mat_R,1],dtype=np.float32)
            y[0] = sub_B[0]
            for m in range(self.Mat_R-1):
                sum_y = 0
                for n in range(m+1):
                    sum_y += self.Mat_L[m+1,n] * y[n]
                y[m+1] = sub_B[m+1] - sum_y
            
            if detail:
                print("此时y为:",y,"\n")
                print("此时的乘积Ly为","\n",np.dot(self.Mat_L,y))

            x = np.zeros([self.Mat_R,1],dtype=np.float32)
            x[self.Mat_R-1] = y[self.Mat_R-1] / self.Mat_U[self.Mat_R-1,self.Mat_R-1]
            for p in range(self.Mat_R-1):
                sum_x = 0
                for q in range(p+1):
                    sum_x += self.Mat_U[-p-2,-q-1] * x[-q-1]
                x[-p-2] = (y[-p-2] - sum_x) / self.Mat_U[-p-2,-p-2]

            if detail:    
                print("此时x为:",x,"\n")
                print("此时的乘积Ux为","\n",np.dot(self.Mat_U,x))
                print("此时的乘积Ax为","\n",np.dot(self.Mat_stable,x))



            self.Mat_inversion[:,i] = x.reshape(self.Mat_R)
        print("逆矩阵为:","\n",self.Mat_inversion)

        return self.Mat_inversion

测试程序

mat = np.array([[1,2,4,17],
            [3,6,-12,3],
            [2,3,-3,2],
            [0,2,-2,6]],dtype=np.float32)

# 创建对象
test = PLU(mat)
# 分解矩阵,当detail=Trur时输出分解过程
test.Factorization(detail=False)
# 输出分解后的P,L,U矩阵
P,L,U = test.get_result()
# 进行矩阵求逆,当detail=True时输出运行过程,以便进行bug定位
Inversion = test.get_inversion(detail=False)

输出

P: 
 [[0. 1. 0. 0.]
 [0. 0. 0. 1.]
 [1. 0. 0. 0.]
 [0. 0. 1. 0.]]
L: 
 [[ 1.          0.          0.          0.        ]
 [ 0.          1.          0.          0.        ]
 [ 0.33333334  0.          1.          0.        ]
 [ 0.6666667  -0.5         0.5         1.        ]]
U: 
 [[  3.   6. -12.   3.]
 [  0.   2.  -2.   6.]
 [  0.   0.   8.  16.]
 [  0.   0.   0.  -5.]]
逆矩阵为: 
 [[ 0.35        0.35       -0.19999997 -1.1       ]
 [-0.375      -0.5416667   1.          1.        ]
 [-0.075      -0.24166667  0.4         0.2       ]
 [ 0.1         0.1        -0.2        -0.1       ]]
  • 4
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值