算法杂记-SVD，PCA和KPCA

Posted 2020-08-05 simplex

SVD

定义

　　假设\(A\)为\(M\times N\)矩阵，则存在\(M\times M\)维正交矩阵\(U=[u_1,u_2,\cdots,u_m]\)，\(N\times N\)维正交矩阵\(V=[v_1,v_2,\cdots,v_n]\)和\(M\times N\)对角矩阵\(\Sigma=diag{\sigma_1,\sigma_2,\cdots,\sigma_p}\)，使得\(A=U^T\Sigma V\)，这种矩阵分解形式成为\(A\)的奇异值分解（SVD），而\(\sigma_i\)称作\(A\)的奇异值。

解释

　　\(M\times N\)的矩阵\(A\)可以看做是\(M\times N\)维线性空间的一个向量。这样我们总能选择\(M\times N\)个基向量将该矩阵展开。最简单的一组基可以选成\(e(i,j)\)，即只在第\(i\)行第\(j\)列元素值为1，其余元素均为0的矩阵。这样任意矩阵在这组基下的展开系数刚好就是对应矩阵元素。

　　理论上来说任选一组完备的基向量都可以用来展开该线性空间中所有的向量，不过对于不同的情况，展开系数的复杂度可能差别非常大。所以我们一般希望根据不同的情况来选择不同的基向量，好比三维空间中的直角坐标，球坐标，柱坐标这些。而SVD分解就给我们提供了一种根据给定矩阵\(A\)来选取比较合适的基向量的做法，同时也给出了矩阵\(A\)在这组基向量下的展开系数。

　　另外，我们知道\(M\times N\)维空间可以由\(M\)维空间和\(N\)维空间做张量积得到，而对应的基向量也是由\(M\)维空间和\(N\)维空间的基向量做张量积得到。由于矩阵\(U\)和\(V\)都是正交的，所以他们的列向量构成了\(M\)维和\(N\)维空间的基向量。故\( u_i v_j^T\)（注意下标不一样）构成\(M\times N\)维空间的一组基向量，而\(\Sigma\)则是用该组基向量展开\(A\)得到的系数，可以认为是基向量的权重。

　　可见此时\(\Sigma\)最多只有\(p\)个非零值，对应的基向量为\( u_i v_i^T\)。相比用\(e(i,j)\)展开要简单很多。根据矩阵乘法易知，\(A=\sum_{i=1}^{r}\sigma_i u_i v_i^T\)，给定了\(U\)和\(V\)，\(A\)就完全由奇异值\(\sigma_i\)来决定，这样我们可以选择丢掉某些很小的奇异值以取得压缩数据的效果。

应用

　　SVD最主要的应用就是对数据进行压缩，只保留最主要的数据。由\(A=\sum_{i=1}^{r}\sigma_i u_i v_i^T\)可知，最多只需要存储矩阵\(U\)和\(V\)的各\(r\)个列向量以及\(r\)个奇异值便可以完全恢复出\(A\)。如果还嫌\(r\)太大，则可以只保留最大的前\(p\)个奇异值及其对应的\(U\)和\(V\)的列向量，因为奇异值\(\sigma_i\)的大小代表了对应基向量\(u_iv_i^{T}\)的权重。

　　另外， 从等式\(A=U^T\Sigma V\)可以导出两个矩阵\(A_1=AV^T=U^T\Sigma=\sum_{i=1}^{r}\sigma_i u_i\)和\(A_2=UA=\Sigma V=\sum_{i=1}^{r}\sigma_i v_i^T\)。这两个矩阵可以认为是对原始矩阵\(A\)进行了列/行压缩得到的\(M\times r\)和\(r\times N\)矩阵。该方法可以用来替代下文提到的PCA对数据进行处理，比如scikit-learn中的pca算法就是基于SVD实现的。

PCA

　　已知\(M\times N\)维数据矩阵\(X=[x_1,x_2,\cdots,x_N]\)，其中\(x_i\)是长度为\(M\)的列向量。我们想要把数据在保留主要信息的情况下进行压缩，方式是寻找一个线性变换 \(V\)，令\(Y=VX\)，变换后的数据\(Y\)应该可以比较清晰的表征出各个数据包含的信息量，使删除某些数据成为可能。我们可以用PCA方法来寻找\(V\)，和SVD不同，PCA的思路是计算\(X\)的协方差矩阵，然后对其进行对角化，寻找协方差矩阵的主轴方向（特征向量），认为主轴长度（即特征值）大的方向上方差也大，从而对应数据信息量也很大。对于较小的主轴长度，可以将其忽略而不至于对数据造成太大影响。具体做法如下：

　　假设\(X\)已经是去中心化的，那么\(X\)的协方差矩阵\(C=XX^T\)或\(C=X^TX\)。后面可以看到\(C\)的不同表示方式实际上实现了对\(X\)不同维度的数据压缩（即对数据点数和特征数做压缩）。我们以前者为例，则\(C\)是\(M\times M\)的矩阵。

　　接下来，对\(C\)进行相似对角化，即求出其\(M\)个特征值\(\lambda_1 \sim \lambda_M\)及其对应特征向量\(v_1 \sim v_M\)，这里假设特征值和特征向量已经按从大到小的顺序进行排列。我们令\(V=[v_1,v_2,\cdots,v_M]\)，计算\(Y=V^TX=[y_1,y_2,\cdots,y_n]\)，其中\(y_i=[v_1^Tx_i,v_2^Tx_i,\cdots, v_m^Tx_i]\)。\(Y\)就是我们想要寻找的数据的新的表示形式。易知\(Y\)满足如下性质：

•     \(y_i\)互不相关，即\(YY^T\)为对角阵（\(YY^T=(V^TX)(V^TX)^T=V^T(XX^T)V\)，刚好是\(C\)的相似对角化矩阵）
•     若\(i>j\)则\(var(y_i)>var(y_j)\)（\(var(y_i)\)是特征值，已经从大到小排列）

　　\(Y\)的计算过程可以认为是把原始数据往各个“主轴”方向上进行投影，而主轴长度则用来衡量各个主轴的权重。这样，权重较小的那些指标就可以被忽略。具体操作起来就是在\(Y\)中删去对应较小特征值的那些特征向量，比如只保留\(r\)个特征向量，那么最后得到的\(y_i\)的长度也只有\(r\)，矩阵\(Y\)是\(r*N\)的，相比之前的\(X\)，实现了在第一维上面的压缩。此外，注意到协方差矩阵是对称阵，所以\(V\)是正交矩阵，要想恢复原始数据，可以直接计算\(X=V^TY=V^TVX\)。当\(V\)已经删除若干特征向量后，恢复的数据也会相应变得简单。

　　如前所述，如果需要对\(X\)的另一维进行压缩，只需要计算协方差矩阵\(X^TX\)即可。


SVD和PCA的区别

　　SVD直接对数据进行分解，不需要数据是方阵。PCA方法需要先求出数据的协方差矩阵（是个方阵），然后对其做特征值分解。此外，用SVD可以直接实现对行/列的压缩，分别对应对数据特征指标和数据量的压缩。而PCA要想实现对行/列的压缩，需要分别求协方差矩阵\(XX^T\)和\(X^TX\)。

KPCA

　　KPCA是为了解决原始数据结构较为复杂，从而线性不可分的情况而发展出来的一种方法。和SVM一样，它也是利用核函数将已知低维数据映射到高维，然后在高维做PCA。整个高维空间的计算过程涉及内积运算，可以利用核函数在原空间进行计算。

　　待压缩数据\(X\)同上文，我们选择一个核函数\(Z: R^M\times R^M \mapsto R\)使得\(Z(x_i,x_j)=<\Phi(x_i),\Phi(x_j)>\)对所有\(x_i,j \in R^M\)均成立。其中\(\Phi\)定义为\(\Phi:R^M \mapsto F\)，将\(R^M\)中的元素\(x\)映射为定义了内积\(<,>\)的希尔伯特空间\(F\)(可能是无穷维的)中的元素\(\Phi(x)\)。\(F\)又叫做特征空间。接下来证明可以无论\(F\)的维度是多少，都可以在其中对数据做PCA。

假设还是打算压缩数据的第一维，即把数据缩减成\(r\times N\)的。根据PCA方法，数据\(X\)在\(F\)空间内的协方差矩阵可以表示为：

\(\bar{C}=\Phi(X)\Phi(X)^T=\sum_{j=1}^{N}\Phi(x_j)\Phi(x_j)^T\)

注意这时矩阵\(\bar{C}\)不再是\(M\times M\)的，而有可能是无穷维的。所以直接对角化行不通，需要寻求变通的方法。

\(\bar{C}\)的特征向量\(v\)满足方程：

\(\bar{C}v=\lambda v\)

将\(\bar{C}=\sum_{j=1}^{N}\Phi(x_j)\Phi(x_j)^T\)带入上式可得：

\(v=\frac{1}{\lambda}\sum_{j=1}^{N}\Phi(x_j)<\Phi(x_j),v>\)

可见特征向量\(v\)属于向量\(\Phi(x_j)\)张成的子空间，因此至多有\(N\)个线性无关的特征向量。

　　考虑等价的方程：

\(\lambda <\Phi(x_j),v>=<\Phi(x_j),\bar{C}v>\)

然后把\(v\)表示成\(\Phi(x_j)\)的展开形式：

\(v=\sum_{j=1}^{N}\alpha_j\Phi(x_j)\)

将上式代入等价的特征值方程可得：

\(\lambda K\alpha=K^2\alpha\)

其中\(K\)是\(N\times N\)矩阵，满足\(K_{ij}=<\Phi(x_i),\Phi(x_j)>\)，\(\alpha=[\alpha_1,\alpha_2,\cdots,\alpha_N]\)

　　所以，我们通过求解特征值方程\(\lambda\alpha=K\alpha\)，便可以求出\(\lambda\)和\(\alpha\)。显然满足\(\lambda\alpha=K\alpha\)的解必然满足\(\lambda K\alpha=K^2\alpha\)，而且可以证明后者的附加解并不会对展开式\(v=\sum_{j=1}^{N}\alpha_j\Phi(x_j)\)有影响。另外， 展开式中的系数\(\alpha_j\)刚好就是高维空间\(F\)中原始数据往主轴方向的投影，即压缩后的数据。这样，我们将原本的求解\(\bar{C}\)的特征向量\(v\)的问题转化成为等价地求\(K\)的特征向量\(\alpha\)，求出了\(\alpha\)便得到了压缩后的数据。但是需要厘清的是\(K\)的的特征值并不是主轴，主轴仍然是\(v\)，\(K\)的意义只是在于帮助我们计算\(F\)空间上的数据\(\Phi(x)\)在\(v\)上的投影。

　　通过令\(<v,v>=1\)，我们可以把对应的特征向量\(\alpha\)进行归一化：

\(1=<v,v>=\sum_{i,j=1}^{N}\alpha_i\alpha_j<\Phi(x_i),\Phi(x_j)>=<\alpha,K\alpha>=\lambda<\alpha,\alpha>\)

　　对于任意新的数据点\(\Phi(x‘)\)，它在主轴\(v\)上的投影为：

\(<v,\Phi(x‘)>=\sum_{j=1}^N\alpha_j<\Phi(x_j),\Phi(x‘)>\)

可见在\(F\)空间进行PCA分解只涉及到计算内积\(<\Phi(x_i),\Phi(x_j)>\)，而这个可以在原空间用核函数计算，所以KPCA是可行的。

　　最后，考虑到在上述推导中，我们假设数据已经是在\(F\)空间上去中心化的，所以对于任意数据还需要对\(\Phi(X)\)进行去中心化。但是由于\(\Phi(X)\)的具体形式未知，所以实际执行起来有困难。解决方法是转而对\(K\)进行去中心化：

　　去中心化后的数据为\(\Phi^c(x)=\Phi(x)-\frac{1}{N}\sum_{i=1}^N\Phi(x_i)\)，据此计算去中心化的矩阵\(K^c_{ij}=<\Phi^c(x_i),\Phi^c(x_j)>\)

　　将\(\Phi^c(x)\)的表达式代入并利用内积的线性性质，最后可得：

\(k^c=K-1_NK-K1_N+1_NK1_N\)

　　其中\(1_N\)即所有元素都为\(1/N\)的\(N\times N\)矩阵。

　　个人觉得KPCA方法的两个关键点一是把协方差矩阵的特征向量\(v\)表示成数据\(\Phi(x)\)的线性叠加，一是利用PCA中压缩后的数据并不是来自特征向量\(v\)，而是来自数据往\(v\)上的投影，即\(<v,\Phi(x)>\)。和PCA不同，KPCA虽然原理上还是要利用协方差矩阵的主轴对数据进行筛选，但是具体计算时并没有对其做分解，而是利用上述两点转而计算矩阵\(K\)。注意原空间的协方差矩阵和\(K\)的维度是不同的，一个是\(M\times M\)，另一个是\(N\times N\)。在映射到高维空间时，协方差矩阵可能成为无穷维，而数据压缩实际上是在\(K\)对应的那一维上进行的。比如，有100组二维数据，想要压缩到1维，PCA是计算\(2\times 2\)的协方差矩阵的特征值，而KPCA则是计算\(100\times 100\)的\(K\)矩阵的特征值，但是会丢弃99个特征向量，只保留1个最长的主轴。