降维--深入理解PCA

1 引言

主成分分析（Principal Component Analysis, PCA）是一种经典又常用的数据降维算法（注意这里的降维是指特征提取 ，有时也称子空间学习，还有一支叫特征选择，有兴趣可参这篇博客），它的主要思想是寻找数据分布方差最大的投影方向，可结合下图理解。

图1

如图1标注出了信号和噪音的方差方向。我们找数据变化大的方向，变化大则含信息量大，所以认为它代表信号方向，而认为小的数据波动是由噪音引起。事实上，有个叫信噪比的指标：$SNR={\sigma }_{signal}^2/{\sigma }_{noise}^2$，$SNR$越大($\gg 1$)，说明信息越准确噪音越少。图中用x,y坐标轴来描述这些数据点的位置，如果我们忽略次要因素，抓主要矛盾，就可以用一个沿着信号方向的轴（投影方向），来刻画这些点，虽然可能会损失一些有用的信息，但是这个过程中我们简化了问题，实现了降维：从二维到一维。

降维的好处：降维可以去除数据中的冗余特征，相应也可减少数据存储空间，更重要的是，方便后续各种算法对数据的进一步处理、分析！

图2

如图2A，把一张“彩色的纸”卷起来。如图2B，在纸上随机选取一些数据。如图2C，把纸摊开。不同的颜色代表不同的类别，在现实生活中，数据的呈现形式往往是这个B样子，同类数据间的欧式距离甚至比非同类的距离都大，所以最好就是找到简单又好分的数据本质的分布空间如C，这种B到C的降维技术就是传说中的流形学习(manifold learning)，典型代表有LLE，LPP，NPE等，在此按住不提。总之流行学习是关注局部结构，而PCA关注的是全局。

2 目标函数及求解

上一节简要介绍了PCA的思想，提到了PCA就是找方差大的方向，这一节就是解决怎样来找这些方向。在机器学习模型中，常针对一个要求解的问题，提出若干准则，以此提出一个目标函数，再在数学上对目标函数进行优化，得出问题的解。PCA的准则是，让投影后的数据分布方差尽可能大，所以PCA的目标函数可以这样定义：
$$\max \frac{1}{n}\sum _{i=1}^n(y_i-\overline{y}{)}^2=\underset{w^{T}w=1}{\max }\frac{1}{n}\sum _{i=1}^n(w^T{x}_i{)}^2=\underset{w^{T}w=1}{\max }{w}^T(\frac{1}{n}X{X}^T)w=\underset{w^{T}w=1}{\max }{w}^T{C}_{x}w\text{\hspace{1em}(2.1)}$$ 其中，$w\in {\mathbb{R}}^d$表示投影方向，$y_i=w^T{x}_i$是投影后的数据点，$X=[x_1,\dots ,x_n]\in {\mathbb{R}}^{d\times n}$是数据矩阵，且我们假设它已经被零均值化处理过，即$\overline{x}=\frac{1}{n}{\sum }_i{x}_i=0$，也自然有$\overline{y}=0$。$C_x$就是协方差矩阵，有时也用修正后的，即$C_x=\frac{1}{n-1}X{X}^T$。关于正交约束，可以这样简单理解，一个作用是使目标式子有解，限制$w$的变化范围，另一个作用是选择标准化的坐标轴，接下来会再提到。

上面是学了一个投影方向（也称子空间），即把数据降成一维，如果要学多个投影方向，很自然地，我们可以让它们的和最大，即：
$$\underset{\forall j,w_j^T{w}_j=1}{\max }\sum _{j=1}^c{w}_j^T{C}_x{w}_j=\underset{W^{T}W=I_c}{\max }Tr(W^T{C}_{x}W)\text{\hspace{1em}(2.2)}$$ 这里，$W=[w_1,\dots ,w_c]\in {\mathbb{R}}^{d\times c}$是投影矩阵，$I_c\in {\mathbb{R}}^{c\times c}$表示单位矩阵，正交约束使得各投影方向自身是标准化的，且它们之间是两两正交的，为什么要正交，可简单类比我们常用x,y轴表示平面上的点，用xyz轴表示三维空间中的点，都是为了使问题简化而已。

现在就要来求解上述提出的目标式子，更一般地，我们直接对式(2.2)进行操作。设$\Lambda \in {\mathbb{R}}^{c\times c}$为拉格朗日乘子，式(2.2)的拉格朗日函数为：
$$J(W)=Tr(W^T{C}_{x}W)-Tr\left(\Lambda (W^{T}W-I_c)\right)\text{\hspace{1em}(2.3)}$$ 求导并置零得：
$$\frac{\partial J(W)}{\partial W}=2C_{x}W-W(\Lambda +{\Lambda }^T)=0\text{\hspace{1em}(2.4)}$$
记$\tilde{\Lambda }=\frac{1}{2}(\Lambda +{\Lambda }^T)$，有$C_{x}W=W\tilde{\Lambda }$。所以，$W$就是$C_X$的特征向量矩阵，$\tilde{\Lambda }$是对应的特征值组成的对角矩阵。$J(W)=Tr(W^{T}W\tilde{\Lambda })=Tr(\tilde{\Lambda })$，要最大化$J(W)$，自然取$C_x$的前$c$个最大的特征值，相应地，$W$就是由$C_x$的前$c$大特征值对应的特征向量按列组成的矩阵。通过投影矩阵$W$，数据从$d$维降低到$c$维。

关于主成分个数($c$)，通常根据需要保留的方差百分比，即需要保留的信息量来确定。比如我们要保留99%的信息量，将$\lambda$从大到小排列，即${\lambda }_1>{\lambda }_2>\dots >{\lambda }_n$，找到最小的$k$，使满足：
$$\frac{{\sum }_{j=1}^k{\lambda }_j}{{\sum }_{j=1}^n{\lambda }_j}\ge 0.99$$

3 再理解

上一节，我们走完了机器学习中传统的目标引出和公式推导过程，这一节中，我们尝试从另外几个角度来剖析PCA，最后殊途同归。

3.1 特征值分解

前面已经引出了协方差的概念，在概率论和统计学中，协方差用于衡量变量之间的总体误差，方差是协方差的特例，即两个变量相同情况下的协方差。在协方差矩阵中，对角线元素对应了各变量自身的方差大小，非对角线元素就是不同变量之间的协方差。如果两个变量是统计独立的，那么两者之间的协方差就是0。

所以我们可以直接拿协方差矩阵做文章，简单分析一下，我们想要的协方差矩阵应该是介个样纸的：对角线元素尽可能大，非对角线元素尽量为0，即协方差矩阵最好是一个对角阵。
原始数据$X$，对应原始协方差矩阵$C_x$，$C_x$通常含有冗余特征，即非对角，现在我们想要通过一个投影变换$Y=W^{T}X$，使得降维后新数据的协方差矩阵$C_y=W^T{C}_{x}W$ (近似)为对角阵。
注意$C_x$实对称，是一个正规矩阵(满足$A^{T}A=A{A}^T,A\in {\mathbb{R}}^{n\times n}$)，根据矩阵论的知识，一定可以对角化：$D=P^T{C}_{x}P$，其中$P\in {\mathbb{R}}^{d\times d}$由$C_x$的特征向量按列组成，且$P$是个正交阵，$D\in {\mathbb{R}}^{d\times d}$是对应的特征值组成的对角矩阵。
嘿！这不就是我们想要的对角矩阵嘛，对比$C_y=W^T{C}_{x}W$，我们可令$W$由$P$的某$c$列按列组成，因为我们想要对角元素尽可能大，这$c$列自然是对应$C_x$的前$c$大特征值的列(特征向量)，跟上一节中推导的结果一致！

3.2 奇异值分解(SVD)

SVD(singular value decomposition)，即奇异值分解，跟上一小节的特征值分解有很大关联，我们不妨先来对比一下两者：

方法	表示	说明
特征值分解	$A=P\Lambda P^{-1}$	$A$是方阵，$P$可逆但不需要是正交阵，$\Lambda$不需要是半正定阵（因为可以有$\lambda <0$）
SVD	$A=U\Sigma V^T$	$A$不需要为方阵， $U,V$是正交阵，$\Sigma$为半正定阵（因为奇异值$\ge 0$）

【注】方阵$A$半正定：$\forall x\in {\mathbb{R}}^n$,有$x^{T}Ax\ge 0$

令$\tilde{Y}=\frac{1}{\sqrt{n}}{X}^T$，则有${\tilde{Y}}^T\tilde{Y}=\frac{1}{n}X{X}^T=C_x$。对$\tilde{Y}$进行SVD：$\tilde{Y}=U\Sigma V^T$，那么$V$就包含了$Y^{T}Y$也即$C_x$的特征向量。SVD这个角度更多地是为了求解的方便（即已知我们要求解的投影方向就是协方差矩阵的特征向量），在数值计算上，SVD比特征值分解要稳定一些，而且对一个$n\times n$的矩阵进行特征分解的计算复杂度高达$O(n^3)$。

补充：不难发现，对于一个对称半正定矩阵$A$（保证了可对角化和特征值$\ge 0$），特征分解和奇异值分解是等价的，因为$A=P\Lambda P^T$既可以看做是特征分解，也可以看做是$U=V$时的奇异值分解，此时$A$的特征值=奇异值。无巧不成书，$C_x$正是这样一个对称半正定阵！所以我们对$C_x$进行特征分解或者SVD都是可以得到投影方向滴。

3.3 数据重构

也可从数据重构的角度来看待PCA，主要思想是将投影后的数据重新投影回原空间，使得数据的重构误差最小，公式化表述就是：
$$\underset{W\in {\mathbb{R}}^{d\times c}}{\min }\Vert X-W{W}^{T}X{\Vert }_F^{2}s.t.w_i^T{w}_j=0(i\ne j),w_i^T{w}_i=1\text{\hspace{1em}(3.1)}$$ 首先，如果我们保留100%的信息量，即选择的主成分个数$c$就等于原始维度$d$，相当于不降维，$W$是$d\times d$大小的正交矩阵，有$W{W}^T=I_d$，此时重构误差为0。下面，我们来证明式(3.1)等价于式(2.2)。

$$\begin{gathered} (3.1)\iff \underset{W^{T}W=I_c}{\min }Tr(X-W{W}^{T}X{)}^T(X-W{W}^{T}X) \\ \iff \underset{W^{T}W=I_c}{\max }Tr(X^{T}W{W}^{T}X) \\ \iff \underset{W^{T}W=I_c}{\max }\frac{1}{n}Tr(W^{T}X{X}^{T}W)\iff (2.2) \end{gathered}$$
【注】上面第一个等价就是按定义展开，第二个是去掉了一些无关的常数项，不影响最大或最小问题的解，第三个用到了迹的性质Tr(AB)=Tr(BA)。

3.4 低秩近似

还有一种观点，认为PCA是找一个低秩矩阵来近似数据矩阵，用欧式距离来衡量的话，就是：
$$\underset{rank(Z)=c}{\min }\Vert X-Z{\Vert }_F^2.\text{\hspace{1em}(3.2)}$$ 根据满秩分解，$Z\in {\mathbb{R}}^{d\times n}$可以分解为：$Z=U{V}^T$，其中$U\in {\mathbb{R}}^{d\times c},V\in {\mathbb{R}}^{n\times c}$，且有$U^{T}U=I_c$。式(2.2)可以改写为：
$$\underset{U^{T}U=I_c,V}{\min }\Vert X-U{V}^T{\Vert }_F^2.\text{\hspace{1em}(3.2)}$$ 这里，$U$可以看做是$d$维空间的$c$个基，$V$可以看做是数据点在该空间中的(低维)表达。下面我们还是来把式(3.2)跟式(2.2)联系起来。式(3.2)对$V$求偏导并置零，可得$V=X^{T}U$，代回(3.2)中可得：
$$\underset{U^{T}U=I_c}{\max }Tr(U^{T}X{X}^{T}U).\text{\hspace{1em}(3.3)}$$ 又回到最初的起点，PCA现在跟低秩近似也联系起来了。

4 拓展

4.1 PCA与Kmeans

PCA是经典的降维算法，Kmeans是经典的聚类算法，看似八竿子打不着，私下会不会有何联系呢。
Kmeans大意是先随机给定一些聚类中心，然后依据数据点距聚类中心的远近给数据分配簇标签，然后计算新的聚类中心，然后再重新分配，然后再计算，再分配…最终聚类中心会移动到一个比较好的位置，使得数据点距聚类中心的距离的平方和比较小。对，就是根据距离的平方和作为准则的，所以目标函数可以这样写：
$$\underset{C_k}{\min }\sum _{k=1}^K\sum _{i\in C_k}\Vert x_i-m_k{\Vert }^2\text{\hspace{1em}(4.1)}$$ 这就是我们常见的Kmeans目标式，其中$K$表示簇的个数，$C_k$表示第$k$个簇的数据点集合，$m_k={\sum }_{i\in C_k}\frac{1}{n_k}{x}_i$表示第$k$个簇的聚类中心，$n_k$表示第$k$个簇中的数据点个数。
仅通过(4.1)好像看不出什么名堂，那我们就来变一变，看能不能写出更简便的表达式，因为是要找跟PCA之间的关系，所以联想PCA的目标式，看能不能跟(2.2)套一下近乎。下面是公式推导：

$$\begin{gathered} (4.1)\iff \min \sum _{k=1}^K\sum _{i\in C_k}(x_i^T{x}_i-2x_i^T{m}_k+m_k^T{m}_k) \\ \iff \min \sum _{i=1}^n{x}_i^T{x}_i-2\sum _{k=1}^K\sum _{i\in C_k}{x}_i^T{m}_k+\sum _{k=1}^K{n}_k{m}_k^T{m}_k \\ \iff \min \sum _{i=1}^n{x}_i^T{x}_i-\sum _{k=1}^K\sum _{i\in C_k}{x}_i^T{m}_k(\because n_k{m}_k^T{m}_k=\sum _{i\in C_k}{x}_i^T{m}_k) \\ \iff \min \sum _{i=1}^n{x}_i^T{x}_i-\sum _{k=1}^K\frac{1}{n_k}\sum _{i,j\in C_k}{x}_i^T{x}_j \\ \iff \min Tr(X^{T}X)-Tr(H^T{X}^{T}XH) \\ \iff \underset{H^{T}H=I,H\ge 0}{\max }Tr(H^{T}KH)(4.2) \end{gathered}$$

其中引入了一个指示矩阵（或称标签矩阵）$H\in {\mathbb{R}}^{n\times k}$，如果$x_i$是第$j$个簇的，则$H_{ij}=1/\sqrt{n_k}$，否则为0。比如现在有5个样例分别属于两个簇，前两个属于第一个簇，后三个属于第二个簇，则$H$长这样：
$${\left[\begin{array}{ccccc}1/\sqrt{2}& 1/\sqrt{2}& 0& 0& 0\\ 0& 0& 1/\sqrt{3}& 1/\sqrt{3}& 1/\sqrt{3}\end{array}\right]}^T$$
$K=X^{T}X$是相似度矩阵，这里是采用的标准内积线性核，也可以扩展到其它核。式(4.2)相比式(4.1)不仅看起来简明多了，而且看上去跟我们目标式(2.2)长得很像，那式(4.2)的解跟式(2.2)的解怎么联系起来呢？还记得上一节中说的特征值分解和SVD分解吧，我们再把它们拿出来瞧瞧（此处忽略$C_x$中的$\frac{1}{n}$，且不考虑$C_x$不可逆的情况）：

$$\begin{gathered} X=U\Sigma V^T \\ {C}_x=X{X}^T=U\Lambda U^T \\ K=X^{T}X=V\Lambda V^T \end{gathered}$$

上面$U\in {\mathbb{R}}^{d\times c}$就是PCA中要求的投影矩阵$W$，$V$是式(4.2)松弛非负约束后的指示矩阵$H$。从第一行中的式子可以知道，$V=X^{T}U{\Sigma }^{-1}$，这说明PCA实际上是在进行一种松弛化的Kmeans聚类，簇的指示矩阵就是投影后的数据矩阵（再缩放一下）。

拓展的拓展：1）kernel Kmeans跟谱聚类（spectral clustering）也有关系；2）$X{X}^T$也叫总体散度矩阵，$XH{H}^T{X}^T$是类间散度矩阵，可参见线性判别分析LDA，它是另一种经典的降维算法。PCA是无监督的，LDA是有监督。所以Kmeas实际上在${\min }_{H}Tr(S_w)$，这里$S_w=S_t-S_b$是类内散度矩阵。3）LDA是已知标签，调整投影矩阵，Kmeans是调整标签。

4.2 Robust PCA

关于鲁棒PCA这个方向，也有很多研究成果，这里只简单地提一种思路，是关于用${\ell }_1$-norm替代${\ell }_2$-norm的，这是一种在机器学习中内十分常用的手法，用于实现鲁棒性或稀疏性。涉及到范数，我们先把式(2.2)改写一下：
$$\underset{W^{T}W=I_c}{\min }\frac{1}{n}\sum _{i=1}^n\Vert W^T{x}_i{\Vert }_2^2\text{\hspace{1em}(4.1)}$$ 那么鲁棒PCA就是：
$$\underset{W^{T}W=I_c}{\min }\frac{1}{n}\sum _{i=1}^n\Vert W^T{x}_i{\Vert }_1\text{\hspace{1em}(4.1)}$$ 关于求解的问题超出本文范畴，感兴趣地可以看相关论文。

主要参考文献

【1】Shlens J. A Tutorial on Principal Component Analysis[J]. 2014, 51(3):219-226.
【2】http://ufldl.stanford.edu/tutorial/unsupervised/PCAWhitening/
【3】Ding C H Q, He X. Principal Component Analysis and Effective K-Means Clustering[C]. SIAM ICDM. DBLP, 2004:126–143.
【4】Masaeli M, Yan Y, Cui Y, et al. Convex Principal Feature Selection[C]. SIAM ICDM. DBLP, 2010:619-628.
【5】Nie F, Yuan J, Huang H. Optimal mean robust principal component analysis[C]. ICML. 2014:1062-1070.
【6】Nie F, Huang H, Ding C, et al. Robust Principal Component Analysis with Non-Greedy L1-Norm Maximization[C]. IJCAI. 2011:1433-1438.