1 回顾特征值分解的几何意义

在上一篇 chat 中，我们讲了通过特征值分解(EVD)的方法对样本的特征提取主成分，从而实现数据的降维。在介绍奇异值分解(SVD)之前，我们再着重挖掘一下特征值分解的几何意义。

1.1 分解过程回顾

我们最开始获得的是一组原始的 m×nm×n 数据样本矩阵 AA ，其中，mm 表示特征的个数， nn 表示样本的个数。通过与自身转置相乘：AATAAT 得到了样本特征的 mm 阶协方差矩阵 CC ，通过求取协方差矩阵 CC 的一组标准正交特征向量 q1,q2,...qmq1,q2,...qm 以及对应的特征值 λ1,λ2,...,λmλ1,λ2,...,λm。

我们这里处理的就是协方差矩阵 CC，对 CC 进行特征值分解，将矩阵分解成了 C=[q1q2...qm] λ1λ2...λm qT1qT2...qTm C=[q1q2...qm][λ1λ2...λm][q1Tq2T...qmT]

最终，我们选取前 kk 个特征向量构成数据压缩矩阵 PP 的各行，通过 PAPA 达到数据压缩的目的。

1.2 几何意义剖析

以上是回顾前文的内容，不难发现，为了完成矩阵的特征值分解，最最关键还是要回归到这个基本性质上来：Cqi=λiqiCqi=λiqi。

我们为什么又提这个呢？结合主成分分析的推导过程我们知道，协方差矩阵 CC 之所以能够分解，是因为在原始空间 RmRm 中，我们原本默认是用 e1,e2,...,eme1,e2,...,em 这组默认基向量来表示我们空间中的任意一个向量 aa，如果我们采用基变换，将 aa 用 q1,q2,...,qmq1,q2,...,qm 这组标准正交基来表示后，CaCa 的乘法运算就变得很简单了，只需要在各个基向量的方向上对应伸长 λiλi 倍即可，如图 1 所示：


实际上，我们之前也重点分析过，因为协方差矩阵具备对称性、正定性，保证了他可以被对角化，并且特征值一定为正，从而使得特征值分解的过程一定能够顺利完成。

因此利用特征值分解进行主成分分析，核心就是获取协方差矩阵，然后对其进行矩阵分解，获得一组特征值和其对应的方向。

2 从 Av=σuAv=σu 入手奇异值分解

但是，如果我们不进行协方差矩阵 CC 的求取，绕开它直接对原始的数据采样矩阵 AA 进行矩阵分解，从而进行降维操作，行不行？

如果继续沿用上面的办法，显然是不行的，特征值分解对矩阵的要求很严，首先得是一个方阵，其次在方阵的基础上，还得满足可对角化的要求。但是原始的 m×nm×n 数据采样矩阵 AA 连方阵这个最基本的条件都不满足，是根本无法进行特征值分解的。

找不到类似 Ap=λpAp=λp 的核心等式了，岂不是无能为力了？怎料，天无绝人之路，这里，我首先给大家介绍一个对于任意 m×nm×n 矩阵的更具普遍意义的一般性质：

对于一个 m×nm×n，秩为 rr 的矩阵 AA，这里我们暂且假设 m>nm>n，于是就有 r≤n

Avi=σiuiAvi=σiui，这个等式非常神奇，我们仔细的揭开里面的迷雾，展现他的精彩之处：矩阵 AA是一个 m×nm×n 的矩阵，他所表示的线性变换是将 nn 维原空间中的向量映射到更高维的 mm 维目标空间中，而 Avi=σiuiAvi=σiui 这个等式意味着，在原空间中找到一组新的标准正交向量 [v1v2...vn][v1v2...vn]， 在目标空间中存在着对应的一组标准正交向量 [u1u2...un][u1u2...un]，此时 vivi 与 uiui 线性无关。当矩阵 AA 作用在原空间上的某个基向量 vivi 上时，其线性变换的结果就是：对应在目标空间中的 uiui 向量沿着自身方向伸长 σiσi 倍，并且任意一对 (vi,ui)(vi,ui) 向量都满足这种关系(显然特征值分解是这里的一种特殊情况，即两组标准正交基向量相等)。如图2所示：


在 Avi=σiuiAvi=σiui 的基础上，我们明白该如何往下走了：A[v1v2...vn]=[σ1u1σ2u2...σnun]A[v1v2...vn]=[σ1u1σ2u2...σnun] ，转换成：

A[v1v2...vn]=[u1u2...un] σ1σ2...σn A[v1v2...vn]=[u1u2...un][σ1σ2...σn]

3 着手尝试分解

此时感觉还差一点，因为我们发现 un+1,un+2,...,umun+1,un+2,...,um 并没有包含在式子里，我们把他们加进去，把 un+1,un+2,...,umun+1,un+2,...,um 加到矩阵右侧，形成完整的 mm 阶方阵 U=[u1u2...unun+1...um]U=[u1u2...unun+1...um]，在对角矩阵下方加上 mnmn 个全零行，形成 m×nm×n 的矩阵 Σ= σ1σ2...σn Σ=[σ1σ2...σn] ，很明显由于 ΣΣ 矩阵最下面的 mnmn 行全零，因此右侧的计算结果不变，等式依然成立。

此时就有：AV=UΣAV=UΣ ，由于 VV 的各列是标准正交向量，因此 V1=VTV1=VT，移到等式右侧，得到了一个矩阵分解的式子：

A=UΣVTA=UΣVT，其中 UU 和 VV 是由标准正交向量构成的 mm 阶和 nn 阶方阵，而 ΣΣ 是一个 m×nm×n 的对角矩阵(注意，不是方阵)。

4 分析分解过程中的细节

从大的框架宏观来看，这个结论非常漂亮，在原空间和目标空间中各找一组标准正交基，就轻轻松松的把对角化的一系列要求轻松化解。直接得到了数据采样矩阵 AA 的矩阵分解形式。

但是此时还有一个最关键的地方似乎没有明确：就是方阵 UU 和方阵 VV 该如何取得，以及 ΣΣ 矩阵中的各个值应该为多少，我们借助在对称矩阵那一节中储备的基础知识来一一化解。

我们还是从 A=UΣVTA=UΣVT 式子入手。首先，获取转置矩阵 AT=(UΣVT)T=VΣUTAT=(UΣVT)T=VΣUT，我们在此基础上可以获取两个对称矩阵：

第一个对称矩阵是： ATA=VΣUTUΣVTATA=VΣUTUΣVT，由于 UU 的各列是标准正交向量，因此 UTU=IUTU=I ，式子最终化简为了：ATA=VΣ2VTATA=VΣ2VT。

同理，第二个对称矩阵：AAT=UΣVTVΣUT=UΣ2UTAAT=UΣVTVΣUT=UΣ2UT。

这里我们结合对称矩阵那一节中的一个重要结论，揭示一下这里面的所有细节： ATAATA 是 nn 阶对称方阵， AATAAT 是 mm 阶对称方阵。他们的秩相等，为 r=Rank(A)r=Rank(A)。因此他们拥有完全相同的 rr 个非零特征值，从大到小排列为：λ1,λ2,...λrλ1,λ2,...λr，两个对称矩阵的剩余 nrnr 个和 mrmr 个特征值为 00。这进一步也印证了 ATA=VΣ2VTATA=VΣ2VT 和 AAT=UΣ2UTAAT=UΣ2UT 对角线上的非零特征值是完全一样的。

同时，由对称矩阵的性质可知： ATAATA 一定含有 nn 个标准正交特征向量，对应特征值从大到小的顺序排列为： [v1v2...vn][v1v2...vn]，而 AATAAT 也一定含有 mm 个标准正交特征向量，对应特征值从大到小依次排列为 [u1u2...um][u1u2...um]。这里的 vivi 和 uiui 一一对应。

对应的 ΣΣ 矩阵也很好求，求出 AATAAT 或 ATAATA 的非零特征值，从大到小排列为： λ1,λ2,...,λrλ1,λ2,...,λr，ΣΣ 矩阵中对角线上的非零值 σiσi 则依次为： √λ1，√λ2,...,√λrλ1，λ2,...,λr。因此，ΣΣ 矩阵对角线上 σrσr 以后的元素均为 00 了。

整个推导分析过程结束，我们隐去零特征值，最终得到了最完美的SVD分解结果： A=[u1u2...um] σ1σ2...σr vT1vT2...vTn A=[u1u2...um][σ1σ2...σr][v1Tv2T...vnT]，这里 r≤n

我们用一个抽象图来示意一下分解的结果，有助于大家加深印象，如图3所示：


由此，我们顺利的得到了任意 m×nm×n 矩阵 AA 的 SVD 分解形式。

下面我们在理论基础上继续推进，讨论如何利用python工具进行求解，以及利用SVD进行数据压缩。

5 行压缩数据降维

我们直接从矩阵 AA 的奇异值分解： A=UΣVTA=UΣVT 入手，分析如何进行行压缩数据降维。

等式两侧同时乘以左奇异矩阵的转置： UTUT，得到：UTA=UTUΣVT=ΣVTUTA=UTUΣVT=ΣVT，重点是左侧的 UTAUTA，我们把矩阵 AA 记作 nn 个 mm 维列向量并排放置的形式，我们展开来看：

[u1u2...um]TA= uT1uT2...uTm [colA1colA2...colAn]= uT1colA1uT1colA2...uT1colAnuT2colA1uT2colA2...uT2colAn............uTmcolA1uTmcolA2...uTmcolAn [u1u2...um]TA=[u1Tu2T...umT][colA1colA2...colAn]=[u1TcolA1u1TcolA2...u1TcolAnu2TcolA1u2TcolA2...u2TcolAn............umTcolA1umTcolA2...umTcolAn]

这是我们刚讲过的基变换方法，很熟悉了吧，colAicolAi 原本使用的是默认的一组基向量 [e1e2...em][e1e2...em] ，我们通过上面的基变换，将其用 [u1u2...um][u1u2...um] 这一组标准正交基来表示，由于这一组标准正交基本质上也是由协方差对称矩阵 AATAAT 得到，因此将各列做基变换后，数据分布从行的角度来看就彼此无关了。

此时，我们可以把每一列看作是一个样本，各行是不同的特征，各行彼此无关，我们可以按照熟悉的方法，选择最大的 kk 个奇异值对应的 kk 个标准正交向量，形成行压缩矩阵：UTk×m= uT1uT2...uTk Uk×mT=[u1Tu2T...ukT]。

通过 UTk×mcolAi= uT1colAiuT2colAi...uTkcolAi Uk×mTcolAi=[u1TcolAiu2TcolAi...ukTcolAi]，就实现了列向量从 mm 维到 kk 维的数据降维，完成了主成分的提取。

6 列压缩数据降维

奇异值分解的精彩之处就在于他可以从两个维度进行数据降维，分别提取主成分，前面介绍的是对行进行压缩降维，那么下面我们来说说如何对列进行压缩降维。

还是牢牢抓住 A=UΣVTA=UΣVT 进行启发，我们对式子两边同时乘以右奇异矩阵 VV，得到 AV=UΣAV=UΣ，我们还是聚焦等式的左侧：AVAV。

我们将其整体进行转置的预处理，得到 (AV)T=VTAT(AV)T=VTAT，我们把矩阵 AA 记作 rowAT1rowAT2...rowATm [rowA1TrowA2T...rowAmT]，那么同样的道理：

VTAT= vT1vT2...vTn [rowA1rowA2...rowAm]= vT1rowA1vT1rowA2...vT1rowAmvT2rowA1vT2rowA2...vT2rowAm............vTnrowA1vTnrowA2...uTnrowAm VTAT=[v1Tv2T...vnT][rowA1rowA2...rowAm]=[v1TrowA1v1TrowA2...v1TrowAmv2TrowA1v2TrowA2...v2TrowAm............vnTrowA1vnTrowA2...unTrowAm]

类比上面的行压缩过程，在 VV 中，我们从大到小取前 kk 个特征值对应的标准正交特征向量，就构成了另一个压缩矩阵：VTk×n= vT1vT2...vTk Vk×nT=[v1Tv2T...vkT]。很明显，通过 VTk×nATVk×nTAT 就能实现将 ATAT 矩阵的各列由 nn 维压缩到 kk 维。而大家不要忘了， ATAT 的列不正是 AA 矩阵的各行吗。

由此，各行向量的维数由 nn 压缩到了 kk 维，实现了列压缩的数据降维。

刚才介绍的行压缩和列压缩进行数据降维在推荐系统中有非常强的实际应用价值，我们在最后一章会举例详细说明。

7 对矩阵整体进行数据压缩

这里我们不谈按行还是按列，我们从矩阵的整体视角，再谈一个数据压缩的方式，我们的思路有点类似级数的概念，将一个 m×nm×n 的原始数据矩阵 AA 分解成若干个同等维度矩阵乘以各自权重后相加的形式。这个如何实现？

依旧从奇异值分解的表达式入手，A=UΣVTA=UΣVT ，我们展开成完整的矩阵形式： A=[u1u2...um] σ1σ2...σr vT1vT2...vTn A=[u1u2...um][σ1σ2...σr][v1Tv2T...vnT]，这里 r≤n

这里展开就得到了：

A=σ1u1vT1+σ2u2vT2+σ3u3vT3+...+σrurvTrA=σ1u1v1T+σ2u2v2T+σ3u3v3T+...+σrurvrT，不难发现，每一个 uivTiuiviT 的结果都是一个等维的 m×nm×n 形矩阵，并且他们彼此之间正交，前面的 σiσi 则是对应矩阵的权重，σ1>σ2>σ3>...>σrσ1>σ2>σ3>...>σr，依序代表了“重要性”的程度，因此我们可以按照主成分贡献率的最低要求，选择前 kk 项进行叠加，用来对原始数据矩阵 AA 进行近似：Aσ1u1vT1+σ2u2vT2+σ3u3vT3+...+σkukvTkAσ1u1v1T+σ2u2v2T+σ3u3v3T+...+σkukvkT

原理示意图如图 4 所示：


这种思想和处理方式在图像压缩的应用中很有用处。

8 利用 python 进行实际操作

8.1 利用 python 进行奇异值分解

举一个 7×57×5 的矩阵 AA 为例，A= 00022000330001111100222005550011100 A=[00022000330001111100222005550011100]，这是一个看上去很有规律的矩阵。

我们这里就不按照推导SVD原理的过程：先求 AATAAT, ATAATA，再接着获得 UU ，VV，ΣΣ 矩阵这样一步步计算。我们通过python提供的工具，直接一次性获得奇异值分解的所有结果。

代码片段：

import numpy as np A=[[0, 0, 0, 2, 2], [0, 0, 0, 3, 3], [0, 0, 0, 1, 1], [1, 1, 1, 0, 0], [2, 2, 2, 0, 0], [5, 5, 5, 0, 0], [1, 1, 1, 0, 0]] U, sigma, VT = np.linalg.svd(A) print(U) print(sigma) print(VT)

运行结果：

[[ 0.00000000e+00 5.34522484e-01 8.41650989e-01 5.59998398e-02 -5.26625636e-02 1.14654380e-17 2.77555756e-17] [ 0.00000000e+00 8.01783726e-01 -4.76944344e-01 -2.09235996e-01 2.93065263e-01 -8.21283146e-17 -2.77555756e-17] [ 0.00000000e+00 2.67261242e-01 -2.52468946e-01 5.15708308e-01 -7.73870662e-01 1.88060304e-16 0.00000000e+00] [ -1.79605302e-01 1.38777878e-17 7.39748546e-03 -3.03901436e-01 -2.04933639e-01 8.94308074e-01 -1.83156768e-01] [ -3.59210604e-01 2.77555756e-17 1.47949709e-02 -6.07802873e-01 -4.09867278e-01 -4.47451355e-01 -3.64856984e-01] [ -8.98026510e-01 5.55111512e-17 -8.87698255e-03 3.64681724e-01 2.45920367e-01 -6.85811202e-17 1.25520829e-18] [ -1.79605302e-01 1.38777878e-17 7.39748546e-03 -3.03901436e-01 -2.04933639e-01 5.94635264e-04 9.12870736e-01]] [ 9.64365076e+00 5.29150262e+00 7.40623935e-16 4.05103551e-16 2.21838243e-32] [[ -5.77350269e-01 -5.77350269e-01 -5.77350269e-01 0.00000000e+00 0.00000000e+00] [ -2.46566547e-16 1.23283273e-16 1.23283273e-16 7.07106781e-01 7.07106781e-01] [ -7.83779232e-01 5.90050124e-01 1.93729108e-01 -2.77555756e-16 -2.22044605e-16] [ -2.28816045e-01 -5.64364703e-01 7.93180748e-01 1.11022302e-16 -1.11022302e-16] [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 -7.07106781e-01 7.07106781e-01]]

这样就非常简单的一次性获得了分解的结果，这里要强调一点，通过程序中获得的 sigma 变量不是一个矩阵，而是由 5 个奇异值按照从大到小顺序组成的列表，而最后一项，打印出来的不是 VV 矩阵，而是转置后的矩阵 VTVT。

8.2 行和列的数据压缩实践

下面我们利用奇异值分解的结果进行行压缩和列压缩：

我们观察这一组奇异值，我们发现前两个奇异值在数量级上占有绝对的优势，因此我们选择 k=2k=2 进行行压缩和列压缩。

依照上面介绍的知识点，我们利用 UT2×7AU2×7TA 将矩阵 AA 的行数由 77 行压缩成了 22 行。利用 VT2×5ATV2×5TAT，将矩阵 ATAT 的行由 55 行压缩成了 22 行，换句话说就是将矩阵 AA 的列由 55 列压缩成了 22 列。

代码片段：

import numpy as np A=[[0, 0, 0, 2, 2], [0, 0, 0, 3, 3], [0, 0, 0, 1, 1], [1, 1, 1, 0, 0], [2, 2, 2, 0, 0], [5, 5, 5, 0, 0], [1, 1, 1, 0, 0]] U, sigma, VT = np.linalg.svd(A) U_T_2x7 = U.T[:2,:] print(np.dot(U_T_2x7,A)) VT_2x5=VT[:2,:] print(np.dot(VT_2x5,np.mat(A).T).T)

运行结果：

[[ -5.56776436e+00 -5.56776436e+00 -5.56776436e+00 0.00000000e+00 0.00000000e+00] [ 3.60822483e-16 3.60822483e-16 3.60822483e-16 3.74165739e+00 3.74165739e+00]] [[ 0.00000000e+00 2.82842712e+00] [ 0.00000000e+00 4.24264069e+00] [ 0.00000000e+00 1.41421356e+00] [ -1.73205081e+00 -7.39557099e-32] [ -3.46410162e+00 -1.47911420e-31] [ -8.66025404e+00 -2.95822839e-31] [ -1.73205081e+00 -7.39557099e-32]]

于是我们成功的分别将矩阵 AA 的行和列进行了压缩。

8.3 利用数据压缩进行矩阵近似

最后，我们来实践一下如何对矩阵 AA 的整体进行数据压缩。同样，我们取前两个奇异值 σ1，σ2σ1，σ2，利用 σ1u1vT1+σ2u2vT2σ1u1v1T+σ2u2v2T 进行矩阵 AA 的近似。

代码片段：

import numpy as np A=[[0, 0, 0, 2, 2], [0, 0, 0, 3, 3], [0, 0, 0, 1, 1], [1, 1, 1, 0, 0], [2, 2, 2, 0, 0], [5, 5, 5, 0, 0], [1, 1, 1, 0, 0]] U, sigma, VT = np.linalg.svd(A) A_1 = sigma[0]*np.dot(np.mat(U[:, 0]).T, np.mat(VT[0, :])) A_2 = sigma[1]*np.dot(np.mat(U[:, 1]).T, np.mat(VT[1, :])) print(A_1+A_2)

运行结果：

[[ -6.97395509e-16 3.48697754e-16 3.48697754e-16 2.00000000e+00 2.00000000e+00] [ -1.04609326e-15 5.23046632e-16 5.23046632e-16 3.00000000e+00 3.00000000e+00] [ -3.48697754e-16 1.74348877e-16 1.74348877e-16 1.00000000e+00 1.00000000e+00] [ 1.00000000e+00 1.00000000e+00 1.00000000e+00 5.19259273e-17 5.19259273e-17] [ 2.00000000e+00 2.00000000e+00 2.00000000e+00 1.03851855e-16 1.03851855e-16] [ 5.00000000e+00 5.00000000e+00 5.00000000e+00 2.07703709e-16 2.07703709e-16] [ 1.00000000e+00 1.00000000e+00 1.00000000e+00 5.19259273e-17 5.19259273e-17]]

从运行结果来看，我们用较少的数据量实现了不错的矩阵近似效果。