关于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}
A−1的每一列
# @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 ]]
