SVD(奇异值分解)Python实现

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

1、SVD算法实现

1.1 SVD原理简单回顾

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

\[A = U\Sigma V^T \tag{1-1} \]

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

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

\[AA^T=U\Sigma V^TV\Sigma^TU^T=U\Sigma \Sigma^TU^T \tag{1-2} \]

\[A^TA=V\Sigma^TU^TU\Sigma V^T=V\Sigma^T\Sigma V^T \tag{1-3} \]

1.2 SVD算法

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

\[AA^T\mathbf{u}_i = \sigma_i \mathbf{u}_i\quad \text{or} \quad AA^T(-\mathbf{u}_i) = \sigma_i (-\mathbf{u}_i)\\ A^TA\mathbf{v}_i = \sigma_i \mathbf{v}_i\quad \text{or} \quad A^TA(-\mathbf{v}_i) = \sigma_i (-\mathbf{v}_i) \]

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

算法1.1:SVD

输入:样本数据
输出:左奇异矩阵,奇异值矩阵,右奇异矩阵
  1. 计算特征值: 特征值分解\(AA^T\),其中\(A \in \mathbf{R}^{m\times n}\)为原始样本数据

    \[AA^T=U\Sigma \Sigma^TU^T \]

    得到左奇异矩阵\(U \in \mathbf{R}^{m \times m}\)和奇异值矩阵\(\Sigma' \in \mathbf{R}^{m \times m}\)

  2. 间接求部分右奇异矩阵:\(V' \in \mathbf{R}^{m \times n}\)

    利用\(A=U\Sigma'V'\)可得

    \[V' = (U\Sigma')^{-1}A = (\Sigma')^{-1}U^TA \tag{1-4} \]

  3. 返回\(U,\ \Sigma',\ V'\),分别为左奇异矩阵,奇异值矩阵,右奇异矩阵。

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

\[ A = U\Sigma' (V')^T\tag{1-5} \]

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

2、SVD的Python实现

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

2.1 SVD实现过程

读取数据

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

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)

特征值分解

# 数据必需先转为浮点型,否则在计算的过程中会溢出,导致结果不准确
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))

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

2.2 SVD降维后重建数据

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

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

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

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

我们平时听到的\(\text{RMSE}=\sqrt{\text{MSE}}\)

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

fig2-2 exp2_svd2.svg
图2-2 原图与降维重建后的图比较

总结

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

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

posted @ 2018-12-03 15:02  EndlessCoding  阅读(31070)  评论(25编辑  收藏  举报