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 += 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)