要求完成课堂上讲的关于矩阵分解的LU、QR（Gram-Schmidt）、Orthogonal Reduction (Householder reduction 和Givens reduction)和 URV程序实现，要求如下：
1、一个综合程序，根据选择参数的不同，实现不同的矩阵分解；在此基础上，实现Ax=b方程组的求解，以及计算A的行列式；
2、可以用matlab、Python等编写程序，需附上简单的程序说明，比如参数代表什么意思，输入什么，输出什么等等，附上相应的例子；
3、一定是可执行文件，例如 .m文件等,不能是word或者txt文档。附上源代码，不能为直接调用matlab等函数库;
import numpy as np
import scipy
import copy
def det_P(P: np.ndarray):
"""
计算LU分解中的P矩阵的行列式
:param P: LU分解中的P矩阵
:return: P矩阵的行列式
"""
h, w = P.shape
t = 0
for i in range(h):
if P[i][i] != 1:
t = t + 1
for j in range(h):
if P[i][j] == 1:
s = j
for j in range(i + 1, h):
if P[j][i] == 1:
P[i][i] = 1
P[i][j] = 0
P[j][i] = 0
P[j]▼显示 = 1
if t % 2 == 0:
return 1
else:
return -1
def det_triangle(matrix: np.ndarray):
"""
计算三角矩阵的行列式
:param matrix: 上三角矩阵或下三角矩阵
:return: 上三角矩阵或下三角矩阵的行列式
"""
h, w = matrix.shape
sum = 1
for i in range(h):
sum = sum * matrix[i][i]
return sum
def det(matrix: np.ndarray):
"""
利用LU分解计算行列式
:param matrix:
:return: 行列式
"""
P, L, U = matrix_decompose(matrix, "LU")
print("*" * 30 + "计算行列式" + "*" * 30)
det_p = det_P(P)
det_l = det_triangle(L)
det_u = det_triangle(U)
return np.rint(det_p * det_u * det_l)
def matrix_decompose(matrix: np.ndarray, mode):
"""
实现矩阵的LU、QR（Gram-Schmidt）、Orthogonal Reduction (Householder reduction 和Givens reduction)和 URV
:param matrix: 需要分解的矩阵
:param mode:
1、LU 实现矩阵的LU分解
2、Gram-Schmidt 采用Gram-Schmidt的方式实现矩阵的QR分解
3、Householder reduction 采用Householder reduction的方式实现矩阵的QR分解
4、Givens reduction 采用Givens reduction的方式实现矩阵的QR分解
5、URV 实现矩阵的URV分解
:return:
若为LU分解，则一次返回P、L、U矩阵，使得PA=LU
若为QR分解，则返回Q、R矩阵，使得A=QR
若为URV分解，则返回R，R，U矩阵，使得A=URV^T
"""
h, w = matrix.shape
if h < w:
print("*" * 30 + mode + "(h<w 类型)" + "*" * 30)
elif h == w:
print("*" * 30 + mode + "(h=w 类型)" + "*" * 30)
else:
print("*" * 30 + mode + "(h>w 类型)" + "*" * 30)
if mode == 'LU':
h, w = matrix.shape
A = copy.deepcopy(matrix)
P = np.identity(h)
L = np.zeros((h, w))
U = np.zeros((w, w))
min_len = min(A.shape)
# 高斯消元法
for i in range(min_len):
index = i
for j in range(i + 1, h):
if np.abs(A[j, i]) > np.abs(A[index, i]):
index = j
temp = copy.deepcopy(P[i, :])
P[i, :] = P[index, :]
P[index, :] = temp
temp = copy.deepcopy(A[i, :])
A[i, :] = A[index, :]
A[index, :] = temp
print("行交换({},{}):
".format(i, index), A)
for j in range(i + 1, h):
c = A[j][i] / A[i][i]
for k in range(i, w):
A[j][k] = A[j][k] - c * A[i][k]
A[j][i] = c
print("行消去({},{}以下):
{}".format(i, i, A))
print("*A*
", A)
for i in range(0, h):
if i < w:
L[i, i] = 1
for j in range(0, min(i, w)):
L[i, j] = A[i, j]
for i in range(0, h):
for j in range(i, w):
U[i, j] = A[i, j]
print("*" * 30 + "分解计算完成" + "*" * 30)
print("P:
", P)
print("L:
", L)
print("U:
", U)
print("*" * 30 + "验证PA=LU" + "*" * 30)
restored_A = L @ U
print("P@A:
", P @ matrix)
print("restored_A:
", np.rint(restored_A))
return P, L, U
elif mode == 'Gram-Schmidt':
u = []
h, w = matrix.shape
A = copy.deepcopy(matrix)
R = np.zeros((w, w))
Q = np.zeros((h, w))
q = []
# 施密特正交化
for i in range(w):
u.append(A[:, i].reshape((-1, 1)))
for i in range(0, w):
sum = np.zeros((h, 1))
for j in range(0, i):
rji = (q[j].T @ u[i])[0][0]
R[j][i] = rji
sum += rji * q[j]
qi = u[i] - sum
rii = np.linalg.norm(qi)
R[i][i] = rii
qi = qi / rii
q.append(qi)
print("h:{},w:{}".format(h, w))
print(Q.shape)
for i in range(w):
Q[:, i] = q[i].reshape(h)
print("*" * 30 + "分解计算完成" + "*" * 30)
print("Q:
", Q)
print("R:
", R)
print("*" * 30 + "验证A=QR" + "*" * 30)
A_restored = Q @ R
A_restored[(A_restored < 1e-12) & (A_restored > -1e-12)] = 0
print("A_restored:
", np.rint(A_restored))
print("matrix:
", matrix)
return Q, R
elif mode == 'Householder reduction':
# 计算PA=R 其中R为上三角矩阵
A_h, A_w = matrix.shape
Plist = []
R = np.zeros((A_h, A_w))
l = min(A_h, A_w)
A = copy.deepcopy(matrix)
if A_h > A_w:
l = l + 1
for i in range(l - 1):
print("epoch:{}".format(i + 1))
min_len = min(A.shape)
print("A的大小为{},A的最短边大小为{}".format(A.shape, min_len))
A_1 = A[:, 0].reshape(-1, 1) # A的第一列
A_1_len = A_1.shape[0]
print("A1的长度:
", A_1_len)
print("A的第一列:
", A_1)
A_1_norm = np.linalg.norm(A_1) # A的第一列的模
print("A的第一列的模:
", A_1_norm)
e1 = np.identity(A_1_len)[:, 0].reshape((-1, 1))
print("e1:
", e1)
u = A_1 - A_1_norm * e1 # u=A1-||A1||e1
print("u=A1-||A1||e1:
", u)
P_1 = np.identity(A_1_len) - 2 * (u @ u.T / (u.T @ u)) # P1=I-2*(UUT/UTU)
print("P1=I-2*(UUT/UTU):
", P_1)
P_ = np.identity(A_h)
P_[i:, i:] = P_1 # 扩充大小
print("P_为P左上角添加1:
", P_)
Plist.append(P_)
A = P_1 @ A # 将A的第一列主元下方的元素消去
A[(A < 1e-12) & (A > -1e-12)] = 0
print("A:
", A)
# 计算R矩阵
R[i, i:] = A[0, :]
R[i:, i] = A[:, 0]
print("R:
", R)
A_ = A[1:, 1:]
print("A_:
", A_)
A = A_
if A_.shape[0] == 1:
R[i + 1][i + 1:] = A[0][0:]
R[(R < 1e-12) & (R > -1e-12)] = 0
print("上三角矩阵R:
", R)
P = np.identity(A_h)
for p in Plist:
P = p @ P
print("P:
", P)
Q = np.linalg.inv(P)
Q[(Q < 1e-12) & (Q > -1e-12)] = 0
print("*" * 30 + "分解计算完成" + "*" * 30)
print("Q:
", Q)
print("R:
", R)
print("*" * 30 + "验证A=QR" + "*" * 30)
A_restored = Q @ R
A_restored[(A_restored < 1e-12) & (A_restored > -1e-12)] = 0
print("A_restored:
", np.rint(A_restored))
print("matrix:
", matrix)
return Q, R
elif mode == 'Givens reduction':
h, w = matrix.shape
Plist = []
A = copy.deepcopy(matrix)
# 构造旋转矩阵
for i in range(0, h):
for j in range(0, min(i, w)):
P_ = np.identity(h)
a = A[j][j]
b = A[i][j]
print("a:{},b:{}".format(a, b))
if a != 0 or b != 0:
d = np.sqrt(a ** 2 + b ** 2)
c = a / d
s = b / d
print("c:{},s:{}".format(c, s))
P_[i][i] = c
P_[j][j] = c
P_[i][j] = -s
P_[j][i] = s
Plist.append(P_)
print("消去({},{})".format(i, j))
print("旋转矩阵P_:
", P_)
A = P_ @ A
A[(A < 1e-12) & (A > -1e-12)] = 0
print("A:
", A)
P = np.identity(h)
for p in Plist:
P = p @ P
print("P:
", P)
Q = np.linalg.inv(P)
R = A
Q[(Q < 1e-12) & (Q > -1e-12)] = 0
print("*" * 30 + "分解计算完成" + "*" * 30)
print("Q:
", Q)
print("R:
", R)
print("*" * 30 + "验证A=QR" + "*" * 30)
A_restored = Q @ R
A_restored[(A_restored < 1e-12) & (A_restored > -1e-12)] = 0
print("A_restored:
", np.rint(A_restored))
print("matrix:
", matrix)
return Q, R
elif mode == 'URV':
A = matrix
m, n = A.shape
r = np.linalg.matrix_rank(A)
# 奇异值分解
U, sigma, VT = np.linalg.svd(A, full_matrices=True)
# 构造C和零矩阵
C = np.diag(sigma[:r])
R = np.zeros((m, n))
R[:r, :r] = C
# 截断U和V以获得正确的大小
U = U[:, :m]
V = VT.T[:, :n]
print("*" * 30 + "验证A=URV^T" + "*" * 30)
print((U @ R @ VT))
return U, R, V.T
else:
raise ValueError("mode error")
if mode in ['Givens reduction', 'Householder reduction', 'Gram-Schmidt']:
print("*" * 30 + "验证QR分解" + "*" * 30)
Q_, R_ = np.linalg.qr(matrix)
print("Q_:
", Q_)
print("R_:
", R_)
elif mode == 'LU':
print("*" * 30 + "验证LU分解" + "*" * 30)
P_, L_, U_ = scipy.linalg.lu(matrix)
print("P_:
", P_)
print("L_:
", L_)
print("U_:
", U_)
def AXb(A: np.ndarray, b: np.ndarray):
"""
解方程AX=b
:param A: 系数矩阵A
:param b: b
:return: 待求未知量X
"""
n = A.shape[0]
Ab = np.hstack((A, b))
Rank_A = np.linalg.matrix_rank(A)
Rank_Ab = np.linalg.matrix_rank(Ab)
if Rank_A == Rank_Ab and Rank_A == n:
P, L, U = matrix_decompose(A, "LU")
print("*" * 30 + "解方程AX=b" + "*" * 30)
# LUX=Pb
Pb = P @ b
print("P@b
", Pb)
# LY=Pb 计算Y
y_list = []
for i in range(n):
sum = 0
for j in range(0, i):
sum += L[i][j] * y_list[j]
y_list.append(Pb[i][0] - sum)
print("y_list:
", y_list)
# UX=Y 计算X
x_list = []
for i in range(n - 1, -1, -1):
sum = 0
for j in range(n - 1, i, -1):
sum += U[i][j] * x_list[n - j - 1]
x_list.append((y_list[i] - sum) / U[i][i])
x_list.reverse()
X = np.array(x_list).reshape((-1, 1))
print("X:
", X)
print("*" * 30 + "验证AX=b" + "*" * 30)
print(np.linalg.solve(A, b))
else:
print("不存在唯一解")
if __name__ == '__main__':
mode = "Householder reduction"
# # h<w 类型
# matrix1 = np.array(
# [[0, -20, 23, 34, 11, 20, 2], [3, 27, 124, 52, 1, 49, 123], [4, 11, 8, 23, 9, 64, 76],
# [12, 45, 12, 86, 95, 19, 99]],
# dtype=np.float32)
# print("matrix1:
{}
matrix1.shape:
{}".format(matrix1, matrix1.shape))
# matrix_decompose(matrix1, mode)
# # h=w 类型
# matrix2 = np.array([[0, -20, 23], [3, 27, -4], [4, 11, -2]], dtype=np.float32)
# print("matrix2:
{}
matrix2.shape:
{}".format(matrix2, matrix2.shape))
# matrix_decompose(matrix2, mode)
# # h>w 类型
# matrix3 = np.array(
# [[0, -20, 23], [3, 27, 124], [4, 11, 8], [12, 45, 12], [123, 1, 80], [12, 99, 662]],
# dtype=np.float32)
# b = np.array([[0], [0], [0], [5], [5], [5]])
# print("matrix3:
{}
matrix3.shape:
{}".format(matrix3, matrix3.shape))
# matrix_decompose(matrix3, b, mode)
# LU分解
matrix4 = np.array(
[[0, -20, 23, 319, -89], [3, 27, -4, 987, 34], [4, 11, -2, -237, 10], [12, 67, 234, 66, -1099],
[32, 562, 64, 88, 24]],
dtype=np.float32)
P, L, U = matrix_decompose(matrix4, "LU")
print(det(matrix4))
print(np.linalg.det(matrix4))
# matrix6 = np.array([[1, 2, -3, 4], [4, 8, 12, -8], [2, 3, 2, 1], [-3, -1, 1, -4]], dtype=np.float32)
# b = np.array(([1], [2], [3], [4]), dtype=np.float32)
# AXb(matrix6, b)