主成分分析 -- Principal Components Analysis

如今，数据挖掘和机器学习常面对高维数据（$X\in {\mathbb{R}}^d$）。为了应对高维数据，数据预处理通常需要对数据做降维处理（$Y=f(X)\in {\mathbb{R}}^p,p\ll d$）。在降低数据维度的同时，要尽可能的保留原有的数据信息。主成分分析便是一种常用的降维方法。它在降维的同时尽可能的保留原数据的方差信息。

假设我们有下图中的二维数据（$X=[x_1,x_2{]}^T$），我们想将其降维到一维，即用一个坐标轴$v$来描述。如图所示，有无数个$v$轴可以选择（图中仅表现3个），选择哪条轴可以最大程度地保留数据$X$的信息呢？PCA选择数据在该方向上投影的范围最广的轴$v_3$，即数据的投影在$v_3$方向上的点最分散（方差最大）。数学表达即：
$$arg\underset{v}{\max }var(v^{T}X)$$

此处，需要引入协方差的概念和一些推导公式。协方差$Cov(X,Y)=E[(X-E[X])(Y-E[Y])]$是计算两个随机变量$X$和$Y$之间的相关性。两个随机变量越线性相关，其协方差越大；线性无关，其协方差为零。

假定我们有原始数据$X_{d\times n}$，其中$n$是样本数量，$d$是数据维度。那么$X$的协方差矩阵$S_X$就是计算数据不同维度间的协方差，如下所示：
$$S_X=\left[\begin{array}{ccc}Cov(X_{1,:},X_{1,:})& \dots & Cov(X_{1,:},X_{d,:})\\ ⋮& & ⋮\\ Cov(X_{d,:},X_{1,:})& \dots & Cov(X_{d,:},X_{d,:})\end{array}\right]$$
计算公式：
$$S_X=\frac{1}{n}\sum _{i=1}^n(X_{:,i}-\mu )(X_{:,i}-\mu {)}^T,where\mu =\frac{1}{n}\sum _{i=1}^n{X}_{:,i}$$
下面我们证明$S_X=\frac{1}{n}\bar{X}{\bar{X}}^T$，其中：$\bar{X}=\left[\begin{array}{ccc}{\bar{X}}_{:,1}& \dots & {\bar{X}}_{:,n}\end{array}\right]$，${\bar{X}}_{:,i}=X_{:,i}-\mu$。

$$\begin{array}{rl}\frac{1}{n}\bar{X}{\bar{X}}^T& =\frac{1}{n}\left[\begin{array}{ccc}{X}_{:,1}-\mu & \dots & X_{:,1}-\mu \end{array}\right]{\left[\begin{array}{ccc}{X}_{:,1}-\mu & \dots & X_{:,1}-\mu \end{array}\right]}^T\\ & =\frac{1}{n}\left[\begin{array}{ccc}{X}_{:,1}-\mu & \dots & X_{:,1}-\mu \end{array}\right]\left[\begin{array}{c}(X_{:,1}-\mu {)}^T\\ ⋮\\ (X_{:,1}-\mu {)}^T\end{array}\right]\\ & =\frac{1}{n}\left[\begin{array}{ccc}{X}_{1,1}-{\mu }_1& \dots & X_{1,n}-{\mu }_1\\ ⋮& & ⋮\\ X_{d,1}-{\mu }_d& \dots & X_{d,n}-{\mu }_d\end{array}\right]\left[\begin{array}{ccc}{X}_{1,1}-{\mu }_1& \dots & X_{d,1}-{\mu }_d\\ ⋮& & ⋮\\ X_{1,n}-{\mu }_1& \dots & X_{d,n}-{\mu }_d\end{array}\right]\\ & =\frac{1}{n}\left[\begin{array}{ccc}{\sum }_{i=1}^n(X_{1,i}-{\mu }_1)(X_{1,i}-{\mu }_1)& \dots & {\sum }_{i=1}^n(X_{1,i}-{\mu }_1)(X_{d,i}-{\mu }_d)\\ ⋮& & ⋮\\ {\sum }_{i=1}^n(X_{d,i}-{\mu }_d)(X_{1,i}-{\mu }_1)& \dots & {\sum }_{i=1}^n(X_{d,i}-{\mu }_d)(X_{d,i}-{\mu }_d)\end{array}\right]\end{array}$$

$$\begin{array}{rl}{S}_X& =\frac{1}{n}\sum _{i=1}^n(X_{:,i}-\mu )(X_{:,i}-\mu {)}^T\\ & =\frac{1}{n}\sum _{i=1}^n\left[\begin{array}{c}{X}_{1;i}-{\mu }_1\\ ⋮\\ X_{d;i}-{\mu }_d\end{array}\right]\left[\begin{array}{ccc}{X}_{1;i}-{\mu }_1& \dots & X_{d;i}-{\mu }_d\end{array}\right]\\ & =\frac{1}{n}\sum _{i=1}^n\left[\begin{array}{ccc}(X_{1,i}-{\mu }_1)(X_{1,i}-{\mu }_1)& \dots & (X_{1,i}-{\mu }_1)(X_{d,i}-{\mu }_d)\\ ⋮& & ⋮\\ (X_{d,i}-{\mu }_d)(X_{1,i}-{\mu }_1)& \dots & (X_{d,i}-{\mu }_d)(X_{d,i}-{\mu }_d)\end{array}\right]\\ & =\frac{1}{n}\left[\begin{array}{ccc}{\sum }_{i=1}^n(X_{1,i}-{\mu }_1)(X_{1,i}-{\mu }_1)& \dots & {\sum }_{i=1}^n(X_{1,i}-{\mu }_1)(X_{d,i}-{\mu }_d)\\ ⋮& & ⋮\\ {\sum }_{i=1}^n(X_{d,i}-{\mu }_d)(X_{1,i}-{\mu }_1)& \dots & {\sum }_{i=1}^n(X_{d,i}-{\mu }_d)(X_{d,i}-{\mu }_d)\end{array}\right]\\ & =\frac{1}{n}\bar{X}{\bar{X}}^T\end{array}$$

我们再回到目标函数$arg\underset{v}{\max }var(v^{T}X)$：
$$\begin{array}{rl}var(v^{T}X)& =\frac{{\sum }_{i=1}^n(v^T{X}_{:,i}-\frac{1}{n}{\sum }_{i=1}^n{v}^T{X}_{:,i}{)}^2}{n}\\ & =\frac{{\sum }_{i=1}^n(v^T{X}_{:,i}-\frac{1}{n}{\sum }_{i=1}^n{v}^T{X}_{:,i})(v^T{X}_{:,i}-\frac{1}{n}{\sum }_{i=1}^n{v}^T{X}_{:,i}{)}^T}{n}\\ & =\frac{{\sum }_{i=1}^n{v}^T(X_{:,i}-\frac{1}{n}{\sum }_{i=1}^n{X}_{:,i})(v^T(X_{:,i}-\frac{1}{n}{\sum }_{i=1}^n{X}_{:,i}){)}^T}{n}\\ & =\frac{{\sum }_{i=1}^n{v}^T(X_{:,i}-\frac{1}{n}{\sum }_{i=1}^n{X}_{:,i})(X_{:,i}-\frac{1}{n}{\sum }_{i=1}^n{X}_{:,i}{)}^{T}v}{n}\\ & =v^T\frac{{\sum }_{i=1}^n(X_{:,i}-\frac{1}{n}{\sum }_{i=1}^n{X}_{:,i})(X_{:,i}-\frac{1}{n}{\sum }_{i=1}^n{X}_{:,i}{)}^T}{n}v\\ & =v^T{S}_{X}v\end{array}$$
即目标函数可转化为：
$$arg\underset{v}{\max }{v}^T{S}_{X}v$$
细心的你可能发现了，上式是一个开口朝上的函数，最大值无穷大。又如何能得到目标函数最大化时$v$的取值呢？仔细观察上面的二维数据图示，我们发现，我们只关心数据在$v$方向的投影，即$v$的指向。$v$向量的大小并不重要。因此，为目标函数增加约束条件（令$v$为单位向量）：
$$arg\underset{v}{\max }{v}^T{S}_{X}vs.t.v^{T}v=1$$
我们可以利用拉格朗日乘数来将上式有约束的优化问题转化为无约束的优化问题：
$$L(v)=-v^T{S}_{X}v+\lambda (v^{T}v-1)$$
然后，再求$L$对$v$的导数（导数为零，求得$L$的最小化时，$v$的取值）：
$$\frac{\partial L}{\partial v}=-2S_{X}v+2\lambda v=0\to S_{X}v=\lambda v$$
注意观察上式右边的式子，是不是和特征值/特征向量的定义很像？至此，我们发现对原始数据$X$降维就是求$X$的协方差矩阵的特征向量。至于$X$的投影最大化时$v$的取值就是特征值$\lambda$取最大值时对应的特征向量（$arg\underset{v}{\max }{v}^T{S}_{X}v=arg\underset{v}{\max }{v}^T\lambda v=arg\underset{v}{\max }\lambda v^{T}v=arg\underset{v}{\max }\lambda$）。

$S_X$是$d\times d$的矩阵，即最多有$d$对特征值和特征向量。PCA中将$S_X$的特征向量叫做主成分（principal component）。特征值第一大的特征向量叫第一主成分，特征值第二大的特征向量叫第二主成分，以此类推。PCA便是希望选择头$p$个最大特征值对应的主成分，在$p$维中描述$d$维中的原始数据$X$（通常$p\ll d$）。可以采用以下计算公式来选择$p$值的大小：
$$期待保留原始数据方差信息的比例=\frac{{\sum }_{i=1}^p{\lambda }_i}{{\sum }_{i=1}^d{\lambda }_i}$$
例如：我们期待保留原始数据$90\%$的方差信息，$\frac{{\sum }_{i=1}^p{\lambda }_i}{{\sum }_{i=1}^d{\lambda }_i}=0.9$。$d$和${\lambda }_i$已知，可将$p$值求解出。请再次观察二维数据$X$的图示，会发现有两个数据点被标注为红色了。原因是，如果只选择一个主成分，即在一维中描述二维数据$X$，这两个点会被投影到一维上的同一个点上。原始数据的部分信息丢失了。所以，数据降维不可避免地会导致数据信息的丢失，但是为了降低数据的复杂度，又不得不对数据进行降维处理。通常，小部分的信息丢失是可以忍受的。

Dual PCA

以上内容是针对维度小于样本数的数据。但是，如今越来越多数据的维度特别高，但是数据量却很少。例如：基因数据、高精度图像等。对于数据$X_{d\times n}$，用PCA对其降维就需要先计算其协方差矩阵$S_X=\frac{1}{n}\bar{X}{\bar{X}}^T$，然后再求解$S_X$的特征值和特征向量。$X$的维度是$d\times n$，$S_X$的维度是$d\times d$。对于$d=100,000$，$n=100$的数据而言，对$S_X$进行特征分解计算量非常大。有没有什么办法能减少计算量呢？

答案是肯定的，奇异值分解（SVD: Singular Value Decomposition）。对于任意矩阵$A_{m\times n}$，我们可以对其做奇异值分解：
$$A=U\Sigma V^T$$
其中，$U$是$A{A}^T$的特征向量，其维度是$m\times m$；$\Sigma$是对角矩阵，对角线上的值是$A{A}^T$和$A^{T}A$的特征值的平方（$A{A}^T$的特征值和$A^{T}A$的特征值相同），呈降序排列，其维度是$m\times n$；$V$是$A^{T}A$的特征向量，其维度是$n\times n$。

我们再整理一下我们的问题。对于维数远大于样本数的数据$X_{d\times n}$，我们想要通过其协方差矩阵的特征向量来对其降维。即想要得到$\frac{1}{n}\bar{X}{\bar{X}}^T$的特征向量，对应于SVD中的$U$。我们对$\bar{X}$奇异值分解公式左右均乘以$V$，得到：
$$\bar{X}V=U\Sigma V^{T}V$$

$V$是$A^{T}A$的特征向量，即$V$的每一列都是单位向量。那么，$V^{T}V=I_{n\times n}$。所以，上式可以转化为：
$$\bar{X}V=U\Sigma$$
再在公式左右均乘以${\Sigma }^{-1}$。得到：
$$\bar{X}V{\Sigma }^{-1}=U$$

上式将对$U$的求解可以转化为对$V$和${\Sigma }^{-1}$的求解。其中，$V$是${\bar{X}}^T\bar{X}$的特征向量，$\Sigma$对角线上是${\bar{X}}^T\bar{X}$的特征值的平方（降序）。

至此，对于$X_{d\times n}$（$d\gg n$），对庞大的协方差矩阵$\frac{1}{n}\bar{X}{\bar{X}}^T$做特征分解，可以转化为对小矩阵${\bar{X}}^T\bar{X}$的特征分解，以此降低计算量。

可以参考案例Start Face Recognition by PCA from scratch来理解PCA。