课程作业随便写写,不一定对
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 += 1
k.append(j)
break
B = 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)#β的单位化得到v
A = []
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)
#求R
C = []
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] = 0
k.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,v2
D = 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 = j
break
j += 1
i = index
index += 1
stop_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]].T
p = 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)