（三）SVD奇异值分解及python实现（待修改）

最新推荐文章于 2023-04-23 17:56:17 发布

alihkj

最新推荐文章于 2023-04-23 17:56:17 发布

阅读量567

点赞数

本文链接：https://blog.csdn.net/alihkj/article/details/104964827

版权

奇异值分解（Singular Value Decomposition）是线性代数中一种重要的矩阵分解，奇异值分解则是特征分解在任意矩阵上的推广。在信号处理、统计学等领域有重要应用。

一原理

对于矩阵 $A_{m \times n}$ ， $\Sigma V^{-1}$
$U_{m \times m},V_{n \times n}$ 是标准正交矩阵，列向量为单位长度，且相互正交， $UU^T=I_m,VV^T=I_n$ 。 $\Sigma_{m \times n}$ 对角矩阵，对角元素非负，且从大到小排列。 $A^TA$ 的特征向量为V的列。 $AA^T$ 的特征向量为U的列。 $A^TA或AA^T$ 的特征值求平方根，构成 $\Sigma$ .
$A^TA V =V\Sigma^TU^TU\Sigma V^T V =V\Sigma ^T \Sigma \\ AA^T U = U\Sigma V^T V\Sigma^TU^T V =U\Sigma \Sigma ^T$
一个m × n的矩阵至多有 p = min(m,n)个不同的奇异值

二算法

输入：样本数据 $A_{m\times n}$
输出：左奇异矩阵，奇异值矩阵，右奇异矩阵
1.特征分解 $(AA^T)u_i=\lambda_iu_i$ ，得到的特征矩阵为左奇异矩阵 $U$ ，特征值对角矩阵为 $\Sigma \Sigma^T$ 。 $\Sigma^{'}=\sqrt{\Sigma \Sigma^T}$ 是 $m\times m$ 的矩阵， $\Sigma_r$ 取前 $r$ 项。
2.由 $A=U\Sigma^{'} V^{T'}$ 得 $V^{T'}=(U\Sigma^{'})^{-1}A= (\Sigma^{'})^{-1}U^TA$ ，求得 $V^{T'}$ 。
3.返回 $U_{m\times m},V^{T'}_{n\times m},\Sigma^{'}_{m\times m}$
3.重建数据 $A=U_{m\times r} \Sigma_{r\times r} V^T_{n\times r}$

三 SVD和PCA对比

数据 $A_{m\times n}$ ，在PCA做降维数据 $A_{m\times n}P_{n\times r}=\tilde{A} _{m\times r} \tag{1}$ 奇异值分解 $A_{m\times n} \approx U_{m\times r} \Sigma_{r\times r} V^T_{n\times r}$ 两边同时右乘以 $V_{n\times r}$ ，得到 $A_{m\times n} V_{n\times r} \approx U_{m\times r} \Sigma_{r\times r} V^T_{n\times r} V_{n\times r}$ $A_{m\times n} V_{n\times r} \approx U_{m\times r} \Sigma_{r\times r} \tag{2}$ $(1) (2)$ 两个式子对比， $P_{n\times r}$ 相当于 $V_{n\times r}$ ，PCA得到新的数据 $\tilde{A} _{m\times r}$ 相当于 $U_{m\times r} \Sigma_{r\times r}$ 。都是对矩阵的列进行了压缩(行压缩使用 $U_{m\times r}^T$ 左乘)。其实PCA几乎可以说是对SVD的一个包装，如果我们实现了SVD，那也就实现了PCA了，而且更好的地方是，有了SVD，我们就可以得到两个方向的PCA，如果我们对 $A^TA$ 进行特征值的分解，只能得到一个方向的PCA。

四举例

奇异值可以被看作成一个矩阵的代表值，或者说，奇异值能够代表这个矩阵的信息。当奇异值越大时，它代表的信息越多。因此，我们取前面若干个最大的奇异值，就可以基本上还原出数据本身。
（1）使用np.linalg.svd

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
 

#读取图片
img_eg = mpimg.imread("2.jpg")
print(img_eg.shape)
#(650, 650, 3)

#奇异值分解
img_temp = img_eg.reshape(650, 650* 3)#先将图片变成600×1200，再做奇异值分解。
U,Sigma,VT = np.linalg.svd(img_temp)#从svd函数中得到的奇异值sigma它是从大到小排列的。

#取前部分奇异值重构图片
# 取前60个奇异值
sval_nums = 60
img_restruct1 = (U[:,0:sval_nums]).dot(np.diag(Sigma[0:sval_nums])).dot(VT[0:sval_nums,:])
img_restruct1 = img_restruct1.reshape((650, 650, 3))
 
# 取前120个奇异值
sval_nums = 120
img_restruct2 = (U[:,0:sval_nums]).dot(np.diag(Sigma[0:sval_nums])).dot(VT[0:sval_nums,:])
img_restruct2 = img_restruct2.reshape((650, 650, 3))

#将照片显示出来，对比效果
fig, ax = plt.subplots(1,3,figsize = (24,32))
 
ax[0].imshow(img_eg)
ax[0].set(title = "src")
ax[1].imshow(img_restruct1.astype(np.uint8))
ax[1].set(title = "nums of sigma = 60")
ax[2].imshow(img_restruct2.astype(np.uint8))
ax[2].set(title = "nums of sigma = 120")

（2）自己构造函数，奇异值矩阵和右奇异矩阵与上述方法不同，但不影响结果

import numpy as np
import pandas as pd
from scipy.io import loadmat
 
# 读取数据，使用自己数据集的路径。
train_data_mat = loadmat("train_data2.mat")
train_data = train_data_mat["Data"]
print(train_data.shape)

# 数据必需先转为浮点型，否则在计算的过程中会溢出，导致结果不准确
train_dataFloat = train_data / 255.0

def SVD(train_dataFloat):
    # 计算特征值和特征向量
    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.diag(np.sqrt(eval_sigma1))
    eval_sigma1_inv = np.linalg.inv(eval_sigma1)
    # 计算右奇异矩阵
    evec_part_vT = eval_sigma1_inv.dot((evec_u.T).dot(train_dataFloat))
    evec_part_v=np.transpose(evec_part_vT)
    return evec_u,eval_sigma1,evec_part_vT

evec_u,eval_sigma1,evec_part_vT=SVD(train_dataFloat)