第八篇：矩阵分解

正定矩阵（后面会详细介绍）

奇异矩阵

正交矩阵

QR分解（正交三角分解）

施密特正交化过程

把一组线性无关向量组化为规范正交向量组，继而得到正交阵。

QR分解步骤

1.写出矩阵的列向量。

2.把列向量按照施密特正交化方法进行正交。

3.把正交化后的向量组再进行单位化。

4.得出矩阵的QR分解。

from sympy import pprint,Symbol,linsolve,solve,symarray,Eq,Expr,roots
from sympy.matrices import Matrix,zeros,diag,eye,GramSchmidt
from sympy.abc import lamda,x,y
A = Matrix([[1,2,2],[2,1,2],[1,2,1]])
print(\'A矩阵是否是可逆的非奇异矩阵：{}\'.format(A.rank() == A.shape[0]))
Q,R = A.QRdecomposition()
print(Q) #正交矩阵
print(R) #非奇异上三角矩阵
print(A == Q*R)
print()

\'\'\'矩阵正交化\'\'\'
B = [Matrix([1,2,2]),Matrix([2,1,2]),Matrix([1,2,1])]
smit = GramSchmidt(B,True) #标准正交，就是正交化后再单位化。
smit_len = len(smit)
for i in range(smit_len-1):
    print(\'向量smit[{}]：\'.format(i))
    pprint(smit[i].norm())
    print(\'simt[{}]的模长是：{}\'.format(i, smit[i].norm()))
    for j in range(i+1,smit_len):
        print(\'向量smit[{}]与向量smit[{}]的点乘等于：{}\'.format(i,j,smit[i].dot(smit[j])))

import numpy as np
import scipy.linalg  as lg
A = np.mat([[1,2,2],[2,1,2],[1,2,1]])
Q,R = np.linalg.qr(A,\'complete\')
q,r = lg.qr(A)

View Code

import numpy as np

def gram_schmidt(A):
    """Gram-schmidt正交化"""
    Q=np.zeros_like(A)
    cnt = 0
    for a in A.T:
        u = np.copy(a)
        for i in range(0, cnt):
            u -= np.dot(np.dot(Q[:, i].T, a), Q[:, i]) # 减去待求向量在以求向量上的投影
        e = u / np.linalg.norm(u)  # 归一化
        Q[:, cnt] = e
        cnt += 1
    R = np.dot(Q.T, A)
    return (Q, R)

def givens_rotation(A):
    """Givens变换"""
    (r, c) = np.shape(A)
    Q = np.identity(r)
    R = np.copy(A)
    (rows, cols) = np.tril_indices(r, -1, c)
    for (row, col) in zip(rows, cols):
        if R[row, col] != 0:  # R[row, col]=0则c=1,s=0,R、Q不变
            r_ = np.hypot(R[col, col], R[row, col])  # d
            c = R[col, col]/r_
            s = -R[row, col]/r_
            G = np.identity(r)
            G[[col, row], [col, row]] = c
            G[row, col] = s
            G[col, row] = -s
            R = np.dot(G, R)  # R=G(n-1,n)*...*G(2n)*...*G(23,1n)*...*G(12)*A
            Q = np.dot(Q, G.T)  # Q=G(n-1,n).T*...*G(2n).T*...*G(23,1n).T*...*G(12).T
    return (Q, R)

def householder_reflection(A):
    """Householder变换"""
    (r, c) = np.shape(A)
    Q = np.identity(r)
    R = np.copy(A)
    for cnt in range(r - 1):
        x = R[cnt:, cnt]
        e = np.zeros_like(x)
        e[0] = np.linalg.norm(x)
        u = x - e
        v = u / np.linalg.norm(u)
        Q_cnt = np.identity(r)
        Q_cnt[cnt:, cnt:] -= 2.0 * np.outer(v, v)
        R = np.dot(Q_cnt, R)  # R=H(n-1)*...*H(2)*H(1)*A
        Q = np.dot(Q, Q_cnt)  # Q=H(n-1)*...*H(2)*H(1)  H为自逆矩阵
    return (Q, R)

np.set_printoptions(precision=4, suppress=True)
A = np.array([[1,1,-1],[1,0,0],[0,1,0],[0,0,1]],dtype=float)

(Q, R) = gram_schmidt(A)
print(Q)
print(R)
print(np.dot(Q,R))

(Q, R) = givens_rotation(A)
print(Q)
print(R)
print(np.dot(Q,R))

(Q, R) = householder_reflection(A)
print(Q)
print(R)
def to_int(x):
    ls = []
    for i in x:
        ls.append(round(i,2))
    return ls
qr = np.dot(Q,R)
QR = np.array(list(map(lambda x:to_int(x),qr)))
print(qr)
print(QR)
print(A == QR)
print(isinstance(A[1,0],float))
print(Q)


from sympy import pprint,Symbol,linsolve,solve,symarray,Eq,Expr,roots
from sympy.matrices import Matrix,zeros,diag,eye
from sympy.abc import lamda,x,y
A = Matrix([[1,1,-1],[1,0,0],[0,1,0],[0,0,1]])
print(\'A矩阵是否是可逆的非奇异矩阵：{}\'.format(A.rank() == A.shape[0]))
Q,R = A.QRdecomposition()
# print(Q) #正交矩阵
# print(R) #非奇异上三角矩阵
pprint(Q.applyfunc(lambda x:x.evalf(n=10)))
print(A == Q*R)

求QR分解的三个方法

from sympy import pprint,Symbol,linsolve,solve,symarray,Eq,Expr,roots
from sympy.matrices import Matrix,zeros,diag,eye,GramSchmidt
from sympy.abc import lamda,x,y
A = Matrix([[1,1,-1],[1,0,0],[0,1,0],[0,0,1]])
print(\'A矩阵是否是可逆的非奇异矩阵：{}\'.format(A.rank() == A.shape[0]))
Q,R = A.QRdecomposition()
print(Q) #正交矩阵
print(R) #非奇异上三角矩阵
pprint(Q.applyfunc(lambda x:x.evalf(n=10)))
print(A == Q*R)

import numpy as np
B = np.mat([[1,1,-1],[1,0,0],[0,1,0],[0,0,1]])
q,r = np.linalg.qr(B)
print(q)
print()
print(r)

SVD分解

矩阵的奇异值分解在矩阵特征值问题，最小二乘法问题及广义逆矩阵问题等方面有重要应用。

矩阵的等价标准型

奇异值的定义

奇异值分解定理

奇异值分解方法—利用矩阵求解

from sympy import pprint,Symbol,linsolve,solve,symarray,Eq,Expr,roots
from sympy.matrices import Matrix,zeros,diag,eye
from sympy.abc import lamda,x,y,z,r
def mySvd(A:Matrix):
    from sympy.matrices import Matrix,diag,eye,zeros
    def sort_eig_vet(sym_eig_vet):
        sort_eig_vets = sorted(sym_eig_vet, key=lambda x: x[0], reverse=True)  # 按照特征值的大小排序特征向量
        V2M = Matrix()  # 生成空矩阵
        for i in sort_eig_vets:
            for v in i[2]:
                V2M = V2M.row_join(v.normalized())
        return V2M
    eig_vets = (A.T * A).eigenvects()  # 求特征向量
    V = sort_eig_vet(eig_vets)
    norm_M = diag(*sorted(A.singular_values(), reverse=True))
    norm_r = norm_M.rank()
    dt = norm_M[0:norm_r, 0:norm_r]
    U = eye(A.shape[0])
    U[:,0:norm_r] = A * V[:, 0:norm_r] * dt.inv()
    VH = V.T
    D = zeros(U.shape[1],VH.shape[0])
    D[:norm_r,:norm_r] = dt
    # ev = (A*A.T).eigenvects()
    # U = sort_eig_vet(ev)
    return U,D,VH

A = Matrix([[1,0,1],[0,1,1],[0,0,0]])
U,D,VH = mySvd(A)
# print(diag(*A.singular_values())) #求A的奇异值
B = U*D*VH
print(A == B)

\'\'\'使用numpy\'\'\'
import numpy as np
B = np.mat([[1,0,1],[0,1,1],[0,0,0]])
u,d,vh = np.linalg.svd(B)
# print(-0.11266613006591797<0.0)
# D = np.zeros((3,3))
# s = np.diag(d)
# D[:3,:3] = s
# # print(vh.shape)
# print(np.allclose(B,np.dot(u,np.dot(D,vh))))


import pandas as pd
from scipy.io import loadmat
# 读取数据，使用自己数据集的路径。
train_data_mat = loadmat("../data/train_data2.mat")
train_data = train_data_mat["Data"]
print(train_data.shape)
# 数据必需先转为浮点型，否则在计算的过程中会溢出，导致结果不准确
train_dataFloat = train_data / 255.0
# 计算特征值和特征向量
eval_sigma1,evec_u = np.linalg.eigh(train_dataFloat.dot(train_dataFloat.T))

#降序排列后，逆序输出
eval1_sort_idx = np.argsort(eval_sigma1)[::-1]
# 将特征值对应的特征向量也对应排好序
eval_sigma1 = np.sort(eval_sigma1)[::-1]
evec_u = evec_u[:,eval1_sort_idx]
# 计算奇异值矩阵的逆
eval_sigma1 = np.sqrt(eval_sigma1)
eval_sigma1_inv = np.linalg.inv(np.diag(eval_sigma1))
# 计算右奇异矩阵
evec_part_v = eval_sigma1_inv.dot((evec_u.T).dot(train_dataFloat))