python实现矩阵的四种分解

课程作业随便写写，不一定对

1.满秩分解

#满秩分解
import numpy as np
import sympy
from sympy import Matrix
from sympy.matrices import dense
from sympy import *
import math
# a = [[1, 2, 1, 1], [2, 1, -2, -2], [1, -1, -4, 3]]
# a = [[1,3,2,1,4],[2,6,1,0,7],[3,9,3,1,11]]
a = [[1,4,-1,5,6],[2,0,0,0,-14],[-1,2,-4,0,1],[2,6,-5,5,-7]]
# Matrix convert to array
A_mat = Matrix(a)
A_arr1 = np.array(A_mat.tolist())
A_rref = np.array(A_mat.rref()[0].tolist())#最简行阶梯矩阵r = np.linalg.matrix_rank(a)#求秩count = 0
k = []#被选中的列向量
for i in range(A_rref.shape[0]):for j in range(A_rref.shape[1]):#遇到1就说明找到了A矩阵的列向量的下标，这些下标的列向量组成B矩阵，然后再找下一行if A_rref[i][j] == 1:count += 1k.append(j)breakB = A_arr1[:,k]
C = A_rref[:r]#秩就是C矩阵的行数B = Matrix(B)
C = Matrix(C)

2.QR分解

#正交三角分解（QR）
a = [[-3, 1, -2],[1, 1, 1],[1, -1, 0],[1, -1, 1]]# a = [[1,1,-1],
#                   [1,0,0],
#                   [0,1,0],
#                   [0,0,1]]
A_mat = Matrix(a)#α向量组成的矩阵A
# A_gs= GramSchmidt(A_mat)
A_arr = np.array(A_mat)
L = []
for i in range(A_mat.shape[1]):L.append(A_mat[:,i])
#求Q
A_gs = GramSchmidt(L)#α的施密特正交化得到β
A_gs_norm = GramSchmidt(L,True)#β的单位化得到vA = []for i in range(A_mat.shape[1]):for j in range(A_mat.shape[0]):A.append(A_gs_norm[i][j])#把数排成一行A_arr = np.array(A)
A_arr = A_arr.reshape((A_mat.shape[0],A_mat.shape[1]),order = 'F')#用reshape重新排列（‘F’为竖着写）
#得到最后的Q
U = Matrix(A_arr)#求RC = []
for i in range(A_mat.shape[1]):for j in range(A_mat.shape[1]):if i > j:C.append(0)elif i == j:t = np.array(A_gs[i])m = np.dot(t.T,t)C.append(sympy.sqrt(m[0][0]))else:t = np.array(A_mat[:,j])x = np.array(A_gs_norm[i])n = np.dot(t.T,x)
#             print(n)C.append(n[0][0])
# C_final为R          
C_arr = np.array(C)
# print(C_arr)
C_arr = C_arr.reshape((A_mat.shape[1],A_mat.shape[1]))
R = Matrix(C_arr)

3.奇异值分解

# 奇异值分解
# a = np.array([[1,1],
#      [0,0],
#      [1,1]])
a = np.array([[2,0],[0,-1j],[0,0]])
a = np.mat(a)#AA^H
W1 = np.dot(a,a.H)
#A^HA
W2 = np.dot(a.H,a)#对任意一个矩阵都有W1和W2是半正定Hermite矩阵，所以特征值大于等于0
m, U = np.linalg.eig(W1)#m为W1的特征值，u为W1特征向量
n, v = np.linalg.eig(W2)#n为W2的特征值，v为特征向量k = []#k用来储存W1中不为0的特征值的特征向量count = 0#计算不等于0的特征值的个数
for i in range(len(m)):if m[i] < 0.000001 and m[i] > 0:m[i] = 0k.append(sqrt(m[i]))if m[i] != 0:count += 1 
k = np.array(k,dtype='float')d = np.mat(np.diag(k[k>0]),dtype='float')#V1必须用公式，W2的特征向量v不一定是V，V由v1和v2组成。v2是W2为0的特征值的特征向量
v1 = np.dot(np.dot(a.H,U[:,0:count]),np.linalg.inv(d.H))V = np.concatenate((v1,v[:,count:len(m) - count]),axis = 1)#拼接v1,v2D = np.diag(k)[:,:a.shape[1]]#D的形状和A矩阵一致print("U = \n{}\n\n D = \n{}\n\n V = \n{}".format(U,D,V))

4.谱分解 

特征值可能会出现2.000000000001和2.0的情况，这样会导致两个数不相等，从而引起后面G有问题，但其实应该处理成是一样的数。

#谱分解
a = [[-2j,4,-2],[-4,-2j,-2j],[2,-2j,-5j]]
# a = [[1/2,0,3j/2],[0,2,0],[-3j/2,0,1/2]]
a = np.mat(a)
#求特征值和特征向量
m,g = np.linalg.eig(a)
g = np.array(g)
for i in range(m.shape[0]):if abs(m[i].real) < 0.000001:m[i] = complex(0,m[i].imag)if abs(m[i].imag) < 0.000001:m[i] = complex(m[i].real,0j)for x in range(g.shape[0]):for y in range(g.shape[1]):if abs(g[x][y].real) < 0.000001:g[x][y] = complex(0,g[x][y].imag)if abs(g[x][y].imag) < 0.000001:g[x][y] = complex(g[x][y].real,0j)# print(g)
# print(m)
k = m.argsort()#将m中的元素从小到大排列，提取其对应的index(索引)，然后输出到k
g = g[:,k]
m = np.array(m[k])
# print(m)
# print(g)
stop_index = [0]#把相等的特征值分片在固定区间内（[0,2,5]表示index0-1元素相等，index2-4相等）j = 1
i = 0
while i < m.shape[0]:
#     print('i = {}'.format(i))while j < m.shape[0]:
#         print('j = {}'.format(j))
#         print(np.round(m[i]) != np.round(m[j]))if abs(m[i] - m[j]) >= 0.000001:#判断2数是否相等，<0.000001表示相等stop_index.append(j)
#             print("stop_index  = {}".format(stop_index))index = jbreakj += 1i = indexindex += 1stop_index.append(m.shape[0])  
G = []
# k = 
#stop_index有多少元素，就有几个相等特征值的区间（区间里相等特征值的个数就是特征值的重根数）
for i in range(len(stop_index) - 1):#     print(i)q = np.zeros([g.shape[1],g.shape[1]])p = g[:,stop_index[i]:stop_index[ i + 1]].Tp = np.matrix(p)
#     print(p)
#     print('+++++++++++++++++++++++++++++++++++++++')for j in p:
#         print(j)q = q + np.dot(j.T,j.T.H)
#         print('np.dot = {}\n'.format(np.dot(j.H,j)))q = np.array(q)
#     print('q = {}'.format(q))for x in range(q.shape[0]):for y in range(q.shape[1]):if abs(q[x][y].real) < 0.000001:q[x][y] = complex(0,q[x][y].imag)G.append(q)
G = np.array(G)

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！