【机器学习】深入剖析主成分分析（PCA）与协方差矩阵

其他机器学习系列文章见于专题：机器学习进阶之路——学习笔记整理，欢迎大家关注。

1. 线性降维方法

  在高维情形下出现的数据样本稀疏、距离计算困难问题，是所有机器学习方法共同面临的严重障碍，被称为“维数灾难”（curse of dimensionality）。

  缓解维数灾难的一个重要途径是降维（dimension reduction），即通过某种数学变换将原始高维属性空间转变为一个低维“子空间”，在这个子空间中样本密度大幅提高，距离计算也变得更为容易。

  对原始高维空间进行线性变换，来进行降维的方法称为线性降维方法。它们都符合以下基本形式：
$$\mathbf{Z}={\mathbf{W}}^{\mathrm{T}}\mathbf{X}$$

  其中，样本矩阵为$\mathbf{X}=\left({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_m\right)\in {\mathbb{R}}^{d\times m}$，投影后样本矩阵为$Z=W^{T}X$，$\mathbf{W}\in {\mathbb{R}}^{d\times d^{\prime }}$是投影矩阵，将维度从$d$降到$d^{\prime }$，$\mathbf{Z}\in {\mathbb{R}}^{d^{\prime }\times m}$是样本在新坐标系中的表达。

  不同线性降维方法的区别在于对低维子空间的性质有不同的要求，相当于对投影矩阵$W$施加不同的约束。我们下面将详细介绍的PCA就是在$W^{T}W=I$的约束下的线性降维方法。

2. PCA概念

  主成分分析（Principal Component Analysis，PCA）是最常用的一种降维方法，通过一个投影矩阵将可能存在相关性和冗余的特征转换为一组更低维度的线性不相关的特征，转换后的特征就叫做主成分。

3. PCA原理

  在降维的过程中，我们希望损失的信息尽可能少，也就是希望保留的信息尽可能多。PCA用方差来度量信息量，在某个维度上，样本分布越分散，方差越大，信息越多。因此，PCA对投影矩阵的第一个要求是使投影后的样本在各维度上方差尽可能大。

  然而，如果单纯只选择方差最大的方向，会导致选择的基向量方向差不多，彼此相关性大，表示的信息几乎是重复的。所以为了使降维后的维度能尽可能地表达信息，第二个要求是不希望投影后的特征之间存在（线性）相关性。

  综上，PCA的优化目标是：

1. 降维后在新维度上的同一维度方差最大；
2. 不同维度之间相关性为0。

  根据线性代数可知，投影值的协方差矩阵的对角线代表了投影值在各个维度上的方差，其他元素代表各个维度之间的相关性（协方差）。基于优化目标，我们希望协方差矩阵是个对角线上值很大的对角矩阵（对角矩阵意味着非对角元素为0）。

  相同维度的多维随机变量$\mathbf{X}={\left[X_1,X_2,\dots ,X_n\right]}^T$和$\mathbf{Y}={\left[Y_1,Y_2,\dots ,Y_n\right]}^T$的协方差矩阵的第$(i,j)$项（第$(i,j)$项是一个协方差）的定义为：
$$\begin{array}{rl}\mathrm{cov}(X_i,Y_j)& =\mathrm{E}[(X_i-\mathrm{E}[X_i])(Y_j-\mathrm{E}[Y_j])]\\ & =\mathrm{E}[X_i{Y}_j-X_i\mathrm{E}[Y_j]-\mathrm{E}[X_i]Y_j+\mathrm{E}[X_i]\mathrm{E}[Y_j]]\\ & =\mathrm{E}[X_i{Y}_j]-\mathrm{E}[X_i]\mathrm{E}[Y_j]-\mathrm{E}[X_i]\mathrm{E}[Y_j]+\mathrm{E}[X_i]\mathrm{E}[Y_j]\\ & =\mathrm{E}[X_i{Y}_j]-\mathrm{E}[X_i]\mathrm{E}[Y_j]\end{array}$$

  如果能对$\mathrm{X}$和$\mathrm{Y}$分别进行中心化，那么$\mathrm{E}[X]=0$，$\mathrm{E}[Y]=0$，将大大简化协方差的计算。这解释了为什么要对数据样本进行中心化。

  假设数据样本已经进行了中心化，得到样本矩阵为$\mathbf{X}=\left({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_m\right)\in {\mathbb{R}}^{d\times m}$，投影后样本矩阵为$Z=W^{T}X$，$\mathbf{W}\in {\mathbb{R}}^{d\times d^{\prime }}$是投影矩阵，将维度从$d$降到$d^{\prime }$，$\mathbf{Z}\in {\mathbb{R}}^{d^{\prime }\times m}$是样本在新坐标系中的表达。投影后的协方差矩阵为：
$$\begin{array}{l}C=\frac{1}{m-1}{\sum }_{i=1}^m{z}_i{z}_i^T\\ =\frac{1}{m-1}Z{Z}^T\\ =\frac{1}{m_1}\left(W^{T}X\right){\left(W^{T}X\right)}^T\\ =\frac{1}{m-1}{W}^{T}X{X}^{T}W\end{array}$$

  在这里必须要声明的是，如果我们用$m$代表样本数，$d$代表维度数。

  1）如果$X$的维度如果是$m\times d$的话，$W$的维度依然是$d\times d^{\prime }$，投影后样本矩阵为$Z=XW$，那么PCA要计算的协方差矩阵是$X^{T}X$。

  2）如果$X$的维度如果是$d\times m$的话，$W$的维度依然是$d\times d^{\prime }$，投影后样本矩阵为$Z=W^{T}X$，那么PCA要计算的协方差矩阵是$X{X}^T$。

  本文按照《机器学习（西瓜书）》里面的写法，采用第二种表达方式，如果要进行代码实现，这种表达方式下的$X$是我们通常理解的样本矩阵的转置。

  要使投影后的协方差矩阵为对角矩阵，就要找到能使投影前的协方差矩阵$X{X}^T$对角化的矩阵$W$（由于协方差矩阵$X{X}^T$是一个实对称矩阵，那么必然存在一个可逆矩阵可使其对角化，且相似对角阵上的元素即为矩阵本身特征值）。然后对投影前协方差矩阵$X{X}^T$进行特征值分解，将求得的特征值排序：${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_d$，取前$d^{\prime }$个特征值对应的特征向量构成$W=\left(w_1,w_2,\dots ,w_{d^{\prime }}\right)$，这就是主成分分析的解。

  严格来说，协方差矩阵是$\frac{1}{m-1}X{X}^T$，但前面常数项不影响，为方便描述我们指的的协方差矩阵是$X{X}^T$。

4. 严格推导过程

  上一节中提到，PCA的优化目标为：

    目标1：降维后同一维度方差最大；

    目标2：不同维度之间相关性为0。

  我们可以将PCA的优化目标转化为在不同维度之间相关性为0的约束条件下，求解使同一维度方差最大化的投影矩阵的问题。

  从两个角度，可以等价地解释方差最大化：

  1. 最近重构性：样本点到投影的超平面（直线的高维推广）的距离都足够近；

  2. 最大可分性：样本点在这个超平面上的投影都尽可能分开。

  根据最近重构性和最大可分性，能够得到主成分分析的两种等价推导。

  （一）基于最近重构性推导PCA

  假定数据样本已经进行了中心化，即${\sum }_{i=1}^m{x}_i=0$，样本点$x_i={\left(x_{i1},x_{i2},\dots ,x_{id}\right)}^T$是$d\times 1$维的列向量，$X=\left(x_1,\dots ,x_i,\dots ,x_m\right)$由$m$个$d$维列向量$x_i$构成，$X$的维度为$d\times m$。如果丢弃部分坐标维度，将维度从$d$降低到$d^{\prime }$，则标准正交基构成的投影矩阵$W=\left(w_1,\dots ,w_i,\dots ,w_{d^{\prime }}\right)$的维度为$d\times d^{\prime }$，由$d^{\prime }$个$d$维列向量$w_i$构成，$w_i={\left(w_{i1},w_{i2},\dots ,w_{id}\right)}^T$是标准正交基向量，${\Vert w_i\Vert }_2=1,w_i^T{w}_j=0(i̸=j)$。投影变换后的新坐标系为$\left(w_1,\dots ,w_i,\dots ,w_{d^{\prime }}\right)$，样本点$x_i$在低维坐标系中的投影为$z_i=\left(z_{i1},z_{i2},\dots ,z_{i{d}^{\prime }}\right)$，$z_{ij}=w_j^T{x}_i$是$x_i$在低维坐标系下第$j$维的坐标。若基于$z_j$来重构$x_i$，则得到${\hat{x}}_i={\sum }_{j=1}^{d^{\prime }}{z}_{ij}{w}_j=W{z}_i$。

  考虑整个训练集，我们想要使原样本点$x_i$与基于投影重构的样本点${\hat{x}}_i$之间的距离为：
$$\begin{array}{l}\sum _{i=1}^m{\Vert {\hat{x}}_i-x_i\Vert }_2^2=\sum _{i=1}^m{\Vert \sum _{j=1}^{d^{\prime }}{z}_{ij}{w}_j-x_i\Vert }_2^2\\ =\sum _{i=1}^m{\left(W{z}_i-x_i\right)}^2\\ =\sum _{i=1}^m\left[{\left(W{z}_i\right)}^T\left(W{z}_i\right)-2{\left(W{z}_i\right)}^T{x}_i+x_i^T{x}_i\right]\\ =\sum _{i=1}^m\left(z_i^T{W}^{T}W{z}_i-2z_i^T{W}^T{x}_i+x_i^T{x}_i\right)\\ \stackrel{W^{T}W=I,z_i=W^T{x}_i}{=}\sum _{i=1}^m\left(z_i^T{z}_i-2z_i^T{z}_i+x_i^T{x}_i\right)\\ =\sum _{i=1}^m\left(-z_i^T{z}_i+x_i^T{x}_i\right)\\ \stackrel{z_i=W^T{x}_i}{=}-\sum _{i=1}^m{\left(W^T{x}_i\right)}^T\left(W^T{x}_i\right)+\sum _{i=1}^m{x}_i^T{x}_i\\ =-\sum _{i=1}^m{x}_i^{T}W{W}^T{x}_i+\sum _{i=1}^m{x}_i^T{x}_i\end{array}\text{\hspace{1em}(1)}$$
  有人可能会有疑问，该式中的$W{W}^T$不是等于$I$吗，那么式（1）不就等于0了吗？注意，$W$是$d\times d^{\prime }$维度的矩阵，不是正交矩阵，尽管$W^{T}W=I$，但$W{W}^T̸=I$，不能和式（1）的第二项抵消。

  由于$x_i^{T}W{W}^T{x}_i$是一个标量，不是一个向量，可以用迹$tr\left(x_i^{T}W{W}^T{x}_i\right)$来代替，(1)式就变成了：
$$\begin{array}{l}-tr\left(\sum _{i=1}^m{x}_i^{T}W{W}^T{x}_i\right)+\sum _{i=1}^m{x}_i^T{x}_i\\ \stackrel{tr\left(AB\right)=tr\left(BA\right)}{=}-tr\left(W^T\left(\sum _{i=1}^m{x}_i{x}_i^T\right)W\right)+\sum _{i=1}^m{x}_i^T{x}_i\\ =-tr\left(W^{T}X{X}^{T}W\right)+\sum _{i=1}^m{x}_i^T{x}_i\end{array}\text{\hspace{1em}(2)}$$

  由于$x_i$的维度是$d\times 1$，$\sum _{i=1}^m{x}_i^T{x}_i$是一个常量，且$W$的每一个向量$w_j$是标准正交基，因此，最小化上式中的距离等价于：
$$\underset{W}{\underbrace{\arg \min }}-\mathrm{tr}\left(W^{T}X{X}^{T}W\right)\text{ s.t. }{W}^{T}W=I\text{\hspace{1em}(3)}$$

  （二）基于最大可分性推导PCA

  样本点$x_i$在空间中超平面上的投影为$W^T{x}_i$，如果要让所有样本点的投影尽可能分开，应使投影后的方差最大化，如图所示。

使所有样本的投影尽可能分开，则需最大化投影点的方差

  由于投影值的协方差矩阵的对角线代表了投影值在各个维度上的方差，则所有维度上的方差和可写成协方差矩阵的迹：（这个过程与第2节中描述的一致）
$$\begin{array}{l}\sum _{i=1}^m\left(W^T{x}_i\right){\left(W^T{x}_i\right)}^T=\sum _{i=1}^m{W}^T{x}_i{x}_i^{T}W\\ =tr\left(W^{T}X{X}^{T}W\right)\end{array}\text{\hspace{1em}(4)}$$

  因此，最大化方差等价于：
$$\underset{W}{\underbrace{\arg \max }}\mathrm{tr}\left(W^{T}X{X}^{T}W\right)\text{ s.t. }{W}^{T}W=I\text{\hspace{1em}(5)}$$

  (5)式与(3)式等价。

5. PCA求解

  对(3)式和(5)式中的优化目标，利用拉格朗日乘子法可得：
$$J(W)=-tr\left(W^{T}X{X}^{T}W\right)+\lambda \left(W^{T}W-I\right)\text{\hspace{1em}(6)}$$
  对$X$求导，由于$\frac{\partial tr\left(A^{T}B\right)}{\partial A}=B$，可得：
$$\frac{\partial J(W)}{\partial W}=-X{X}^{T}W+\lambda W\text{\hspace{1em}(7)}$$
令导数为0，得：
$$X{X}^{T}W=\lambda W\text{\hspace{1em}(8)}$$

  根据线性代数中的特征值分解$Ax=\lambda x$可知上式是一个类似的问题。于是，只需要对协方差矩阵$X{X}^T$进行特征值分解，将求得的特征值排序：${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_d$，取前$d^{\prime }$个特征值对应的特征向量构成$W=\left(w_1,w_2,\dots ,w_{d^{\prime }}\right)$，这就是主成分分析的解。

   降维后低维空间的维数$d^{\prime }$通常是事先指定的，还可以设置一个重构阈值，例如$t=95$，然后选取使下式成立的最小$d^{\prime }$值：
$$\frac{{\sum }_{i=1}^{d^{\prime }}{\lambda }_i}{{\sum }_{i=1}^d{\lambda }_i}\ge t\text{\hspace{1em}(9)}$$

6. PCA算法描述

   输入：训练数据集 D = { x 1 , x 2 , ⋯   , x m } D= \left\{ x _ { 1 } , x _ { 2 } , \cdots , x _ { m } \right\} D={x1​,x2​,⋯,xm​}，低维空间维数$d^{\prime }$。

   过程：

      （1）对所有样本进行中心化：
$$x_i=x_i-\frac{1}{m}\sum _{j=1}^m{x}_j$$

      （2）计算样本的协方差矩阵$X{X}^T$；

      （3）对协方差矩阵$X{X}^T$做特征值分解，求出其特征值及其对应的特征向量；

      （4）取最大的$d^{\prime }$个特征值对应的特征向量$\left(w_1,w_2,\dots ,w_{d^{\prime }}\right)$；

      （5）对样本集$D$中的每一个样本$x_i$，转化为新样本$z_i=W^T{x}_i$，得到输出样本集 D ′ = { z 1 , z 2 , ⋯   , z m } D'= \left\{ z _ { 1 } , z _ { 2 } , \cdots , z_ { m } \right\} D′={z1​,z2​,⋯,zm​}。

   输出：投影矩阵$W=\left(w_1,w_2,\dots ,w_{d^{\prime }}\right)$和降维后的样本集 D ′ = { z 1 , z 2 , ⋯   , z m } D'= \left\{ z _ { 1 } , z _ { 2 } , \cdots , z_ { m } \right\} D′={z1​,z2​,⋯,zm​}。

7. PCA优缺点

  PCA算法的主要优点：

1. 舍弃一部分信息量少的特征， 使样本的采样密度增大；
2. 最小的特征值所对应的特征向量往往与噪声有关，舍弃它们可以起到去噪的作用（在信号处理中任务信号具有较大方差，噪声拥有较小方差）；
3. 计算方法简单，主要运算是特征值分解，易于实现。

  主要缺点：

1. 主成分各个特征维度的含义具有一定的模糊性，不如原始样本特征的解释性强；
2. 方差小的非主成分也可能含有对样本差异的重要信息，因降维丢弃可能对后续数据处理有影响。

参考文献：

1. 《机器学习》第十章降维与度量学习——周志华
2. 主成分分析（PCA）原理总结
3. 协方差矩阵与PCA深入原理剖析
4. PCA算法原理：为什么用协方差矩阵
5. 主成分分析（PCA）原理详解
6. 图文并茂的PCA教程