注：在《SVD（奇异值分解）小结 》中分享了 SVD 原理，但其中只是利用了numpy.linalg.svd函数应用了它，并没有提到如何自己编写代码实现它，在这里，我再分享一下如何自已写一个SVD 函数。但是这里会利用到 SVD 的原理，如果大家还不明白它的原理，可以去看看《SVD（奇异值分解）小结 》，或者自行百度 /google。数据集：https://pan.baidu.com/s/1ZmpUSIscy4VltcimwwIWew。

1、SVD 算法实现#

1.1 SVD 原理简单回顾#

有一个$m\times n$的实数矩阵$A$，我们可以将它分解成如下的形式

其中$U$和$V$均为单位正交阵，即有$U{U}^T=I$和$V{V}^T=I$，$U$称为左奇异矩阵，$V$称为右奇异矩阵，$\Sigma$仅在主对角线上有值，我们称它为奇异值，其它元素均为 0。上面矩阵的维度分别为$U\in R^{m\times m},\Sigma \in R^{m\times n}$,$V\in R^{n\times n}$。

正常求上面的$U,V,\Sigma$不便于求，我们可以利用如下性质

1.2 SVD 算法#

据1.1 小节，对式（1-3）和式（1-4）做特征值分解，即可得到奇异值分解的结果。但是样分开求存在一定的问题，由于做特征值分解的时候，特征向量的正负号并不影响结果，比如，我们利用式（1-3）和（1-4）做特征值分解

如果在计算过程取，取上面的$u_i$组成左奇异矩阵$U$，取$-v_i$组成右奇异矩阵$V$，此时$A\ne U\Sigma V^T$。因此求$v_i$时，要根据$u_i$来求，这样才能保证$A=U\Sigma V^T$。因此，我们可以得出如下 1.1 计算 SVD 的算法。它主要是先做特性值分解，再根据特征值分解得到的左奇异矩阵$U$间接地求出部分的右奇异矩阵$V^{\prime }\in R^{m\times n}$。

1. 计算特征值： 特征值分解$A{A}^T$，其中$A\in R^{m\times n}$为原始样本数据

   $$A{A}^T=U\Sigma {\Sigma }^T{U}^T$$

   得到左奇异矩阵$U\in R^{m\times m}$和奇异值矩阵${\Sigma }^{\prime }\in R^{m\times m}$

2. 间接求部分右奇异矩阵： 求$V^{\prime }\in R^{m\times n}$

   利用$A=U{\Sigma }^{\prime }{V}^{\prime }$可得

   $$\begin{array}{}{V}^{\prime }=(U{\Sigma }^{\prime }{)}^{-1}A=({\Sigma }^{\prime }{)}^{-1}{U}^{T}A& \text{(1-4)}\end{array}$$

3. 返回$U,{\Sigma }^{\prime },V^{\prime }$，分别为左奇异矩阵，奇异值矩阵，右奇异矩阵。

注：这里得到的${\Sigma }^{\prime }$和$V^{\prime }$与式（1-2）所得到的$\Sigma ,V$有区别，它们的维度不一样。${\Sigma }^{\prime }$是只取了前$m$个奇异值形成的对角方阵，即${\Sigma }^{\prime }\in R^{m\times m}$；$V^{\prime }$不是一个方阵，它只取了$V\in R^{m\times n}$的前$m$行（假设$m<n$），即有$V^{\prime }=V(:m,\cdot )$。这样一来，我们同样有类似式（1-1）的数学关系成立，即

$$\begin{array}{}A=U{\Sigma }^{\prime }(V^{\prime }{)}^T& \text{(1-5)}\end{array}$$

我们可以利用此关系重建原始数据。

2、SVD 的 Python 实现#

以下代码的运行环境为python3.6+jupyter5.4。

2.1 SVD 实现过程#

读取数据#

这里面的数据集大家随便找一个数据就好，如果有需要我的数据集，可以下在面留言。

Copy Highlighter-hljs
import numpy as np
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)

特征值分解#

Copy Highlighter-hljs
# 数据必需先转为浮点型，否则在计算的过程中会溢出，导致结果不准确
train_dataFloat = train_data / 255.0
# 计算特征值和特征向量
eval_sigma1,evec_u = np.linalg.eigh(train_dataFloat.dot(train_dataFloat.T))

计算右奇异矩阵#

Copy Highlighter-hljs
#降序排列后，逆序输出
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))

上面的计算出的evec_u, eval_sigma1, evec_part_v分别为左奇异矩阵，所有奇异值，右奇异矩阵。

2.2 SVD 降维后重建数据#

取不同个数的奇异值，重建图片，计算出均方误差，如图 2-1 所示。从图中可以看出，随着奇异值的增加，均方误差（MSE）在减小，且奇异值和的比率正快速上升，在 100 维时，奇异值占总和的 53%。

图 2-1 奇值分解维度和均方误差变化图

注： 均方误差 MSE 有如下计算公式

$$\text{MSE}=\frac{1}{n}\left(\right(y_1-y_1^{\prime }{)}^2+(y_2-y_2^{\prime }{)}^2+\cdots +(y_n-y_n^{\prime }{)}^2)$$

我们平时听到的$\text{RMSE}=\sqrt{\text{MSE}}$。

将图和 10、50、100 维的图进行比较，如图 2-2 所示。在直观上，100 维时，能保留较多的信息，此时能从图片中看出车辆形状。

图 2-2 原图与降维重建后的图比较

总结#

SVD 与特征值分解（EVD）非常类似，应该说 EVD 只是 SVD 的一种特殊怀况。我们可以通过它们在实际的应用中返过来理解特征值 / 奇异值的含义：特征值 / 奇异值代表着数据的信息量，它的值越大，信息越多。

最近作业是真的多呀，冒着生命危险来分享，希望能给大家带来帮助😄