通过Python实现矩阵分解
更新于2022年2月22日 浏览：1553
矩阵分解可以一定程度上避免进行求逆带来的庞大计算，是在矩阵规模不大的时候所采用的方法，但矩阵规模较大的时候还是比较推荐采用迭代的方法进行求解的，这里主要讲一讲常用的矩阵分解方法，然后把代码给出来，同样的代码基于python并且封装成了一个类，可以直接实例化并且调用，方法主要按照网上来写的，稳定性不作保证。
这里主要有Choleski分解，LU分解，QR分解，高斯消去和SVD分解，同时对于上三角和下三角的矩阵求解我单独列了出来，因为在上面的分解方法很多都是分解为这两种东西然后求解减少计算量的。
# -*- coding: utf-8 -*-"""Created on 2022/1/25 10:45@Author: 123@Email: 985040320@qq.com@Describe: **加入文件描述, 要实现的功能, 注意事项等**"""import timeimport numpy as npclass DirectSolver:
def __init__(self, A, b):
self.A = A.copy()
self.b = b.copy()
self.n = len(self.A)
self.upper = np.zeros((self.n, self.n))
self.lower = np.zeros((self.n, self.n))
def CholeskiSolver(self):
# 分解 [A] = [L] [L^T]
self.b = self.b.flatten()
for k in range(self.n):
self.A[k, k] = np.sqrt(self.A[k, k] - np.dot(self.A[k, 0:k], self.A[k, 0:k]))
for i in range(k + 1, self.n):
self.A[i, k] = (self.A[i, k] - np.dot(self.A[i, 0:k], self.A[k, 0:k])) / self.A[k, k]
for k in range(1, self.n):
self.A[0:k, k] = 0.0
# 求解 [L]{y} = {b}
for k in range(self.n):
self.b[k] = (self.b[k] - np.dot(self.A[k, 0:k], self.b[0:k])) / self.A[k, k]
# 求解 [L^T]{x} = {y}
for k in range(self.n - 1, -1, -1):
self.b[k] = (self.b[k] - np.dot(self.A[k + 1:self.n, k], self.b[k + 1:self.n])) / self.A[k, k]
return self.b
def LUSolver(self):
self.upper[0, :] = self.A[0, :]
for i in range(1, self.n + 1):
self.lower[i - 1, i - 1] = 1.0
for k in range(1, self.n):
if abs(self.upper[k - 1, k - 1]) > 1.0e-10:
for i in range(k + 1, self.n + 1):
# 下三角分解
for j in range(1, i):
total = 0
for l in range(1, j):
total = total - self.lower[i - 1, l - 1] * self.upper[l - 1, j - 1]
self.lower[i - 1, j - 1] = (self.A[i - 1, j - 1] + total) / self.upper[j - 1, j - 1]
# 上三角分解
for j in range(1, self.n + 1):
total = 0
for l in range(1, i):
total = total - self.lower[i - 1, l - 1] * self.upper[l - 1, j - 1]
self.upper[i - 1, j - 1] = self.A[i - 1, j - 1] + total
else:
print('有0向量在第', k, '行')
break
# 前向计算得到y
self.b = self.LowerSolve(self.lower, self.b)
# 后向计算得到x
self.b = self.UpperSolve(self.upper, self.b)
return self.b.flatten()
def QRSolver(self):
# householder转化为正交矩阵 Q 与非奇异上三角矩阵 R 的乘积
(r, c) = np.shape(self.A)
Q = np.identity(r)
R = np.copy(self.A)
for i in range(r - 1):
x = R[i:, i]
e = np.zeros_like(x)
e[0] = np.linalg.norm(x)
u = x - e
v = u / np.linalg.norm(u)
Q_i = np.identity(r)
Q_i[i:, i:] -= 2.0 * np.outer(v, v) # 求个外集组成矩阵
R = np.dot(Q_i, R)
Q = np.dot(Q, Q_i)
x = np.arange(c)
B = np.dot(Q.T, self.b)
# x = np.linalg.solve(R, B) #采用np的求解器作为测试
B = self.UpperSolve(R, B)
return B.flatten()
def GaussSolver(self):
# 采用高斯消去法求解
B = np.hstack((self.A, self.b))
for i in range(0, (B.shape[0]) - 1):
if B[i, i] == 0:
print("终断运算：")
break
else:
for j in range(i + 1, B.shape[0]):
B[j:j + 1, :] = B[j:j + 1, :] - (B[j, i] / B[i, i]) * B[i, :]
# 创建矩阵存放答案 初始化为0
x = np.array(np.zeros(B.shape[0], dtype=float))
size = x.shape[0] - 1
b = size + 1
x[size] = B[size, b] / B[size, size]
for i in range(size - 1, -1, -1):
try:
x[i] = (B[i, b] - np.sum(np.multiply(B[i, i + 1:b], x[i + 1:b]))) / (B[i, i])
except:
print("错误")
return x
def SvdSolver(self):
[k, f, c] = np.linalg.svd(self.A, full_matrices=1, compute_uv=1)
v = np.transpose(c)
d = np.eye(len(k)) * 1 / f
u = np.transpose(k)
x_1 = np.matmul(v, d)
x_2 = np.matmul(x_1, u)
x = np.matmul(x_2, b)
return x.flatten()
# 上三角矩阵求解程序
@staticmethod
def UpperSolve( A, b):
A = A.copy()
b = b.copy()
for i in range(len(A[0, :]), 0, -1):
total = b[i - 1]
if i < len(A):
for j in range(i + 1, len(A) + 1):
total = total - A[i - 1, j - 1] * b[j - 1]
b[i - 1] = total / A[i - 1, i - 1]
return b
# 下三角矩阵求解程序
@staticmethod
def LowerSolve(A, b):
A = A.copy()
b = b.copy()
for i in range(1, len(A) + 1):
total = b[i - 1]
if i > 1:
for j in range(1, i):
total = total - A[i - 1, j - 1] * b[j - 1]
b[i - 1] = total / A[i - 1, i - 1]
return bif __name__ == "__main__":
# cholesky矩阵分解算法
A = np.array([[4., -1., 1.], [-1., 4.25, 2.75], [1., 2.75, 5.]])
b = np.array([[4.], [6.], [7.25]])
t1 = time.time()
cls = DirectSolver(A, b) # 创建一个求解器的实例cls
x = cls.CholeskiSolver() # 调用Choleski法求解
t2 = time.time()
print("CholeskiSolver", x, f"用时{t1-t2}s")
# LU矩阵分解算法
cls = DirectSolver(A, b) # 创建一个求解器的实例cls
t1 = time.time()
x = cls.LUSolver() # 调用LU法求解
t2 = time.time()
print("LUSolver", x, f"用时{t1-t2}s")
# QR矩阵分解算法
cls = DirectSolver(A, b)
t1 = time.time()
x = cls.QRSolver()
t2 = time.time()
print("QRSolver", x, f"用时{t1-t2}s")
# 高斯消去法
cls = DirectSolver(A, b)
t1 = time.time()
x = cls.GaussSolver()
t2 = time.time()
print("GaussSolver", x, f"用时{t1-t2}s")
# SVD法
cls = DirectSolver(A, b)
t1 = time.time()
x = cls.SvdSolver()
t2 = time.time()
print("SvdSolver", x,f"用时{t1-t2}s")
t1 = time.time()
x=np.linalg.solve(A, b).flatten()
t2 = time.time()
print("numpy自带求解器", x, f"用时{t1-t2}s")
喜欢的朋友可以给个关注或者联系我
默认 最新
当前暂无评论，小编等你评论哦！点赞 评论 收藏
关注
