主成分分析（PCA）的原理

前言

主成分分析（Principal components analysis，以下简称PCA）是最重要的降维方法之一。在数据压缩消除冗余和数据噪音消除等领域都有广泛的应用。

一、PCA的思想

PCA 顾名思义，就是找出数据里最主要的方面，用数据里最主要的方面来代替原始数据。具体的，假如我们的数据集是 $n$ 维的，共有 $m$ 个数据 $(x^{(1)},x^{(2)},...,x^{(m)})$ 。我们希望将这 $m$ 个数据的维度从 $n$ 维降到 $n^{\prime }$ 维，希望这 $m$ 个 $n^{\prime }$ 维的数据集尽可能的代表原始数据集。我们知道数据从 $n$ 维降到 $n^{\prime }$ 维肯定会有损失，但是我们希望损失尽可能的小。那么如何让这 $n^{\prime }$ 维的数据尽可能表示原来的数据呢？

我们先看看最简单的情况，也就是 $n=2，n^{\prime }=1$ ，即将数据从二维降维到一维。如下图所示。我们希望找到某一个维度方向，它可以代表这两个维度的数据。图中列了两个向量方向，$u_1$ 和 $u_2$ ，那么哪个向量可以更好的代表原始数据集呢？从直观上也可以看出，$u_1$ 比 $u_2$ 好。

为什么 $u_1$ 比 $u_2$ 好呢？可以有两种解释，第一种解释是样本点到这个直线（向量）的距离足够近，第二种解释是样本点在这个直线上的投影能尽可能的分开。

假如我们把 $n^{\prime }$ 从 $1$ 维推广到任意维，则我们希望的降维的标准为：样本点到这个超平面的距离足够近，或者说样本点在这个超平面上的投影能尽可能的分开。

基于上面的两种标准，我们可以得到PCA的两种等价推导。

二、PCA的推导：基于最小投影距离

我们首先看第一种解释的推导，即样本点到这个超平面的距离足够近。

假设 $m$ 个 $n$ 维数据 $(x^{(1)},x^{(2)},...,x^{(m)})$ 都已经进行了中心化（使用原始数据减去原始数据的均值，得到的数据即为中心化后的数据，几何意义是将原始数据做了一个平移，平移之后数据的中心是(0, 0)），即 ${\Sigma }_{i=1}^m{x}^{(i)}=0$ 。经过投影变换后得到的新坐标系为 { w 1 , w 2 , . . . , w n } \{w_1, w_2, ..., w_n\} {w1​,w2​,...,wn​} ,其中 $w$ 是标准正交基，即 $||w|{|}_2=1,w_i^T{w}_j=0$ 。

如果我们将数据从 $n$ 维降到 $n^{\prime }$ 维，即丢弃新坐标系中的部分坐标，则新的坐标系为 { w 1 , w 2 , . . . , w n ′ } \{w_1, w_2, ..., w_{n'}\} {w1​,w2​,...,wn′​} ，样本点 $x^{(i)}$ 在 $n^{\prime }$ 维坐标系中的投影为：$z^{(i)}=(z_1^{(i)},z_2^{(i)},...,z_{n^{\prime }}^{(i)}{)}^T$ 。其中， $z_j^{(i)}=w_j^T{x}^{(i)}$ 是 $x^{(i)}$ 在低维坐标系里第 $j$ 维的坐标。

如果我们用 $z^{(i)}$ 来恢复原始数据 $x^{(i)}$ ,则得到的恢复数据 ${\overline{x}}^{(i)}={\sum }_{j=1}^{n^{\prime }}{z}_j^{(i)}{w}_j=W{z}^{(i)}$，其中，$W$ 为标准正交基组成的矩阵。

现在我们考虑整个样本集，我们希望所有的样本到这个超平面的距离足够近，即最小化下式：
$$\sum _{i=1}^m||{\overline{x}}^{(i)}-x^{(i)}|{|}_2^2$$
将这个式子进行整理，可以得到：
$$\sum _{i=1}^m||{\overline{x}}^{(i)}-x^{(i)}|{|}_2^2=\sum _{i=1}^m||W{z}^{(i)}-x^{(i)}|{|}_2^2\text{⋯⋯①}$$
$$=\sum _{i=1}^m(W{z}^{(i)}{)}^T(W{z}^{(i)})-2\sum _{i=1}^m(W{z}^{(i)}{)}^T{x}^{(i)}+\sum _{i=1}^m(x^{(i)}{)}^T{x}^{(i)}\text{⋯⋯②}$$
$$=\sum _{i=1}^m(z^{(i)}{)}^T{z}^{(i)}-2\sum _{i=1}^m(z^{(i)}{)}^T{W}^T{x}^{(i)}+\sum _{i=1}^m(x^{(i)}{)}^T{x}^{(i)}\text{⋯⋯③}$$
$$=\sum _{i=1}^m(z^{(i)}{)}^T{z}^{(i)}-2\sum _{i=1}^m(z^{(i)}{)}^T{z}^{(i)}+\sum _{i=1}^m(x^{(i)}{)}^T{x}^{(i)}\text{⋯⋯④}$$
$$=-\sum _{i=1}^m(z^{(i)}{)}^T(z^{(i)})+\sum _{i=1}^m(x^{(i)}{)}^T{x}^{(i)}\text{⋯⋯⑤}$$
$$=-\sum _{i=1}^{m}tr((z^{(i)})(z^{(i)}{)}^T)+\sum _{i=1}^m(x^{(i)}{)}^T{x}^{(i)}\text{⋯⋯⑥}$$
$$=-tr\left(W^T\right(\sum _{i=1}^m{x}^{(i)}(x^{(i)}{)}^T)W)+\sum _{i=1}^m(x^{(i)}{)}^T{x}^{(i)}\text{⋯⋯⑦}$$
$$=-tr(W^{T}X{X}^{T}W)+\sum _{i=1}^m(x^{(i)}{)}^T{x}^{(i)}\text{⋯⋯⑧}$$
其中：

第①步用到了 ${\overline{x}}^{(i)}=W{z}^{(i)}$ ；

第②步用到了平方和展开式；

第③步用到了矩阵转置公式 $(AB{)}^T=B^T{A}^T$ 和 $W^{T}W=I$ ；

第④步用到了 $z^{(i)}=W^T{x}^{(i)}$ ；

第⑤步合并同类项；

第⑥步用到了向量 $z^{(i)}$ 的内积等于矩阵 $(z^{(i)})(z^{(i)}{)}^T$ 的迹；

第⑦步用到了 $z^{(i)}=W^T{x}^{(i)}$ ；

第⑧步用到了 ${\sum }_{i=1}^m{x}^{(i)}(x^{(i)}{)}^T=X{X}^T$ ，其中 $X$ 为特征值矩阵。

注意到第⑦步中 ${\sum }_{i=1}^m{x}^{(i)}(x^{(i)}{)}^T$ 是数据集的特征间的协方差矩阵，$W$ 的每一个向量 $w_j$ 是标准正交基。而 ${\sum }_{i=1}^m(x^{(i)}{)}^T{x}^{(i)}$ 是一个常量。最小化上式等价于：
$$\underset{W}{\underbrace{\arg \min }}-tr(W^{T}X{X}^{T}W)\text{subject to WTW=I}$$
这个最小化不难，直接观察也可以发现最小值对应的 $W$ 由协方差矩阵 $X{X}^T$ 最大的 $n^{\prime }$ 个特征值对应的特征向量组成。当然用数学推导也很容易。利用拉格朗日函数可以得到：
$$J(W)=-tr(W^{T}X{X}^{T}W)+\lambda (W^{T}W-I)$$
对 $W$ 求导有 $-X{X}^{T}W+\lambda W=0$ ，整理后得：
$$X{X}^{T}W=\lambda W$$
这样可以更清楚的看出，$W$ 为 $X{X}^T$ 的 $n^{\prime }$ 个特征向量组成的矩阵，而 $\lambda$ 为 $X{X}^T$ 的若干特征值组成的矩阵，特征值在主对角线上，其余位置为 $0$ 。当我们将数据集从 $n$ 维降到 $n^{\prime }$ 维时，需要找到最大的 $n^{\prime }$ 个特征值对应的特征向量。这 $n^{\prime }$ 个特征向量组成的矩阵 $W$ 即为我们需要的矩阵。对于原始数据集，我们只需要用 $z^{(i)}=W^T{x}^{(i)}$ ，就可以把原始数据集降维到最小投影距离的 $n^{\prime }$ 维数据集。

从这里我们就可以看出来PCA是求前k个最大的特征值对应的特征向量。

三、PCA的推导：基于最大投影方差

假设 $m$ 个 $n$ 维数据 $(x^{(1)},x^{(2)},...,x^{(m)})$ 都已经进行了中心化，即 ${\sum }_{i=1}^m{x}^{(i)}=0$ 。经过投影变换后得到的新坐标系为 { w 1 , w 2 , . . . , w n } \{w_1, w_2, ..., w_n\} {w1​,w2​,...,wn​} ，其中 $w$ 是标准正交基，即 $||w|{|}_2=1,w_i^T{w}_j=0$ 。

如果我们将数据从 $n$ 维降到 $n^{\prime }$ 维，即丢弃新坐标系中的部分坐标，则新的坐标系为 $w_1,w_2,...,w_{n^{\prime }}$ ，样本点 $x^{(i)}$ 在 $n^{\prime }$ 维坐标系中的投影为：$z(i)=(z_1^{(i)},z_2^{(i)},...,z_{n^{\prime }}^{(i)}{)}^T$ 。其中，$z_j^{(i)}=w_j^T{x}^{(i)}$ 是 $x^{(i)}$ 在低维坐标系里第 $j$ 维的坐标。

对于任意一个样本 $x^{(i)}$，在新的坐标系中的投影为 $W^T{x}^{(i)}$ ,在新坐标系中的投影方差为 $W^T{x}^{(i)}(x^{(i)}{)}^{T}W$ ，要使所有的样本的投影方差和最大，也就是最大化 ${\sum }_{i=1}^m{W}^T{x}^{(i)}(x^{(i)}{)}^{T}W$ 的迹，即：
$$\underset{W}{\underbrace{\arg \max }}tr(W^{T}X{X}^{T}W)\text{subject to WTW=I}$$
观察第二小节的基于最小投影距离的优化目标，可以发现完全一样，只是一个是加负号的最小化，一个是最大化。

我们利用拉格朗日函数可以得到：
$$J(W)=tr(W^{T}X{X}^{T}W)+\lambda (W^{T}W-I)$$
对 $W$ 求导有 $X{X}^{T}W+\lambda W=0$ ，整理后得：
$$X{X}^{T}W=(-\lambda )W$$
和上面一样可以看出，$W$ 为 $X{X}^T$ 的 $n^{\prime }$ 个特征向量组成的矩阵，而 $-\lambda$ 为 $X{X}^T$ 的若干特征值组成的矩阵，特征值在主对角线上，其余位置为 $0$ 。当我们将数据集从 $n$ 维降到 $n^{\prime }$ 维时，需要找到最大的 $n^{\prime }$ 个特征值对应的特征向量。这 $n^{\prime }$ 个特征向量组成的矩阵 $W$ 即为我们需要的矩阵。对于原始数据集，我们只需要用 $z^{(i)}=W^T{x}^{(i)}$，就可以把原始数据集降维到最小投影距离的 $n^{\prime }$ 维数据集。

四、PCA算法流程

从上面两小节我们可以看出，求样本 $x^{(i)}$ 的 $n^{\prime }$ 维的主成分其实就是求样本集的协方差矩阵 $X{X}^T$ 的前 $n^{\prime }$ 个特征值对应的特征向量矩阵 $W$ ，然后对于每个样本 $x^{(i)}$，做如下变换 $z^{(i)}=W^T{x}^{(i)}$，即达到PCA的降维目的。

下面我们看看具体的算法实现流程。

输入：$n$ 维样本集 $D=(x^{(1)},x^{(2)},...,x^{(m)})$ ，要降维到的维数 $n^{\prime }$ 。

输出：降维后的样本集 $D^{\prime }$ 。

• ① 对所有的样本进行中心化： $x^{(i)}=x^{(i)}-\frac{1}{m}{\sum }_{j=1}^m{x}^{(j)}$ ；

• ② 计算样本的协方差矩阵 $X{X}^T$ ；

• ③ 对矩阵 $X{X}^T$ 进行特征值分解；

• ④ 取出最大的 $n^{\prime }$ 个特征值对应的特征向量 $(w_1,w_2,...,w_{n^{\prime }})$， 将所有的特征向量标准化后，组成特征向量矩阵 $W$ ；

• ⑤ 对样本集中的每一个样本 $x^{(i)}$，转化为新的样本 $z^{(i)}=W^T{x}^{(i)}$ ；

• ⑥ 得到输出样本集 $D^{\prime }=(z^{(1)},z^{(2)},...,z^{(m)})$ 。

有些时候，我们不指定降维后的 $n^{\prime }$ 的值，而是换种方式，指定一个降维到的主成分比重阈值 t 。这个阈值 $t$ 在 (0,1] 之间。假如我们的 $n$ 个特征值为 ${\lambda }_1\ge {\lambda }_2\ge ...\ge {\lambda }_n$ ，则 $t$ 可以通过下式得到：
$$t\le \frac{{\sum }_{i=1}^{n^{\prime }}{\lambda }_i}{{\sum }_{i=1}^n{\lambda }_i}$$

五、PCA实例

下面举一个简单的例子，说明PCA的实现过程。

假设我们的数据集有 10 个二维数据(2.5, 2.4), (0.5, 0.7), (2.2, 2.9), (1.9, 2.2), (3.1, 3.0), (2.3, 2.7), (2, 1.6), (1, 1.1), (1.5, 1.6), (1.1, 0.9)，需要将它们用PCA降到1维特征。

首先我们对样本中心化，这里样本的均值为 (1.81, 1.91)，所有的样本减去这个均值向量后，即中心化后的数据集为 (0.69, 0.49), (-1.31, -1.21), (0.39, 0.99), (0.09, 0.29), (1.29, 1.09), (0.49, 0.79), (0.19, -0.31), (-0.81, -0.81), (-0.31, -0.31), (-0.71, -1.01) 。

现在我们开始求样本的协方差矩阵，由于我们是二维的，则协方差矩阵为：
$$X{X}^T=\left(\begin{array}{cc}cov(x_1,x_1)& cov(x_1,x_2)\\ cov(x_2,x_1)& cov(x_2,x_2)\end{array}\right)$$
对于我们的数据源，求出协方差矩阵为：
$$X{X}^T=\left(\begin{array}{cc}0.616555556& 0.615444444\\ 0.615444444& 0.716555556\end{array}\right)$$
求出特征值为 (0.0490833989, 1.28402771) ，对应的特征向量分别为： $(0.735178656,0.677873399{)}^T,(-0.677873399,-0.735178656{)}^T$ ，由于最大的 $k=1$ 个特征值为 1.28402771 ，对应的前 $k=1$ 个特征向量为 $(-0.677873399,-0.735178656{)}^T$ ，则我们的 $w=(-0.677873399,-0.735178656{)}^T$ 。

我们对所有的数据集进行投影 $z^{(i)}=W^T{x}^{(i)}$ ，得到PCA降维后的10个一维数据集为：(-0.827970186， 1.77758033， -0.992197494， -0.274210416， -1.67580142， -0.912949103， 0.0991094375， 1.14457216, 0.438046137， 1.22382056) 。

六、核主成分分析KPCA介绍

在上面的PCA算法中，我们假设存在一个线性的超平面，可以让我们对数据进行投影。但是有些时候，数据不是线性的，不能直接进行PCA降维。这里就需要用到和支持向量机一样的核函数的思想，先把数据集从 n 维映射到线性可分的高维 $N>n$ ，然后再从 N 维降维到一个低维度 n' , 这里的维度之间满足 $n^{\prime }<n<N$ 。

使用了核函数的主成分分析一般称之为核主成分分析(Kernelized PCA，以下简称KPCA)。假设高维空间的数据是由 n 维空间的数据通过映射 $\varphi$ 产生。

则对于 n 维空间的特征分解：
$$\sum _{i=1}^m{x}^{(i)}(x^{(i)}{)}^{T}W=\lambda W$$
映射为：
$$\sum _{i=1}^m\varphi (x^{(i)})\varphi (x^{(i)}{)}^{T}W=\lambda W$$
通过在高维空间进行协方差矩阵的特征值分解，然后用和PCA一样的方法进行降维。一般来说，映射 $\varphi$ 不用显式的计算，而是在需要计算的时候通过核函数完成。由于KPCA需要核函数的运算，因此它的计算量要比PCA大很多。

七、总结

作为一个非监督学习的降维方法，PCA只需要将特征值分解，就可以对数据进行压缩，去噪。因此在实际场景应用很广泛。为了克服PCA的一些缺点，出现了很多PCA的变种方法，比如第六小节中的为解决非线性降维的KPCA，还有解决内存限制的增量PCA方法Incremental PCA，以及解决稀疏数据降维的Sparse PCA等。

PCA算法的主要优点有：

• ① 仅仅需要用方差来衡量数据包含的信息量，不受数据集以外的因素影响。

• ② 各主成分之间正交，可消除原始数据成分间的相互影响的因素。

• ③ 计算方法简单，主要运算是特征值分解，易于实现。

主要缺点有：

• ① 主成分各个特征维度的含义具有一定的模糊性，不如原始样本特征的解释性强。

• ② 方差小的非主成分也可能含有对样本差异的重要信息，因降维丢弃后可能对后续数据处理有影响。

八、参考文献

① 主成分分析（PCA）原理总结 - 刘建平
② 机器学习中的数学(2) – 矩阵中心化、标准化的意义和作用