Fenrier Lab

理解主成分分析(1)——最大方差投影与数据重建

一方面,高维数据的存储和计算对计算机的性能提出了更高的要求,另一方面,对于人来讲,高维数据提供的信息不太直观,所以如何对数据进行降维一直是一项很有意义的话题。

最大方差投影与主成分

让我们先从一个二维数据集开始,设

\[S = \{x_i\mid x_i \in R^2 , i=1,2,,,n\}\]

把上述点集画到平面上,假设产生了下面的图像

现在我们想要为每个点赋予一个具有代表性的值来替代它的坐标,从而降低该数据集的维度,并且这个被赋予的值应该在最大程度上表现出被它替代的点的特征。那什么是点的特征呢?对于一个点来讲,与其他点的不同就是它所拥有的特征,所以如果它的替代值能尽可能地与其他点的替代值相异,那么这种替换就是有意义的。当然,这种降维方式损失了不少信息,但我们先不管这个。一个想法是,将点投影到坐标轴上,然后用坐标轴上的值作为点的替代值,如下面显示的图像

前一张的图像以横轴作为投影直线,后一张以纵轴为投影直线。不难发现,将点集投影到横轴上具有更强的分散性,即点与点的区分度越大。那么怎样来衡量这种区分度呢?我们很自然地会想到方差,显然方差越大越好,上述数据的横坐标方差大于纵坐标方差,所以如果非要对点集 \(S\) 进行降维,那么把它们投影到横轴上比投影到纵轴上要更合理一点。这样一来我们得到了一个新的数据集合,它是由原数据集的所有横坐标组成的。

但是现在有一种情况,那就是给定另外一数据集,如果坐标轴(无论是横轴还是纵轴)不是数据的最佳投影轴,如下图所示,事实上,此时最佳投影轴是斜率为1的直线

此时就不能单纯的用点的横纵坐标值来替代了,一种可行的方案是,先把点投影到斜率为 1 的直线上,然后计算投影点与原点的带符号距离,这个带符号距离就最佳的替代值。但其实找到这个值不用这么麻烦,考虑一条与上述投影线相垂直的直线,假设它的方程为

\[\omega^T x = 0\]

其中令 \(\mid \omega\mid=0\) ,容易证明

\[f(x) =\omega^T x\]

其实就是点 \(x\) 到直线 \(\omega^T x = 0\) 的带符号距离。所以从另一个角度来看,将点映射成它到某直线的带符号距离也是一种降维方法,并且和投影方法是等价的,前者的优势在于计算更简单。

当拿到样本数据集之后,一开始并不知道具体的目标直线,所以只有任意构造一条直线,然后计算得到一系列带符号距离

\[Y = [y_1 \quad y_2 \quad\cdots \quad y_n]\]

其中 \(y_i=f(x_i)\)。于是我们的目标其实就是寻找 \(\omega\) 使得 \(Var(Y)\) 具有最大值。为了一般性讨论,我们考虑更高维的数据,假设特征维度为 d ,那么

\[y_i =\omega^T x_i = \omega_1 x_{i1} + \omega_2 x_{i2} + ... + \omega_d x_{id}\]

写成矩阵形式为

\[Y = \left[ \begin{aligned} y_1 \\y_2\\\vdots\\y_n \end{aligned} \right] =\omega_1\left[ \begin{aligned} x_{11} \\x_{21}\\\vdots\\x_{n1} \end{aligned} \right] + \omega_2\left[ \begin{aligned} x_{12} \\x_{22}\\\vdots\\x_{n2} \end{aligned} \right] + ... + \omega_d\left[ \begin{aligned} x_{1d} \\x_{2d}\\\vdots\\x_{nd} \end{aligned} \right]\]

如果再定义

\[X_i = \left[ \begin{aligned} x_{1i} \\x_{2i}\\...\\x_{ni} \end{aligned} \right]\]

那么就有

\[Y =\omega \cdot X = \sum_{i=1}^d \omega_i X_i\]

我们在这里将所有特征点作为一个矩阵来考虑,即

\[X = [X_1 \quad X_2 \quad ...\quad X_d]\]

然后假设其协方差矩阵为

\[\Sigma = \left[ \begin{aligned} \sigma_{11}\quad&\sigma_{12}\quad..&\sigma_{1d} \\ \sigma_{21}\quad&\sigma_{22}\quad..&\sigma_{2d}\\ ..\quad&..\,\,\,\quad..&..\\ \sigma_{d1}\quad&\sigma_{d2}\quad..&\sigma_{dd}\\ \end{aligned} \right]\]

其中 \(\sigma_{ij} = Cov(X_1, C_2)\)。下面我们将证明 \(Var(Y)\) 的最大值为 \(\Sigma\) 的最大特征值。根据方差的性质

\[Var(Y) =\sum_{i=1}^d\sum_{j=1}^d \omega_i \omega_jCov(X_i,X_j)= \sum_{i=1}^d\sum_{j=1}^d \omega_i \omega_j \sigma_{ij} = \omega^T \Sigma \omega\]

容易知道协方差矩阵 \(\Sigma\) 为实对称矩阵,于是它拥有一组可以作为单位正交基的特征向量

\[e_1,e_2...e_d\]

满足

\[e_i^T e_j = \left\{ \begin{aligned} 1\quad i = j\\0\quad i \ne j \end{aligned}\right.\]

以及相对应的特征值

\[\lambda_1, \lambda_2 ... \lambda_d\]

并不失一般性的假设

\[\lambda_1 \ge \lambda_2 \ge ... \ge \lambda_d\]

利用这组基向量,我们可以将 \(\omega\) 分解成

\[\omega = \sum_{j=1}^d \alpha_j e_j\]

那么

\[Var(Y) = \sum_{i=1}^d \sum_{j=1}^d \alpha_i \alpha_j e_i^T \Sigma e_j = \sum_{i=1}^d \alpha_i^2 \lambda_i \le \lambda_1\sum_{i=1}^n\alpha_i^2\]

其中 \(\lambda_1\) 是最大的特征值,上述公式取等号的条件是,与 \(\lambda_1\) 相对应的系数 \(\alpha_1\) 为 1,其他的系数为 0,这时有 \(\omega = e_1\)。再考虑到约束

\[\|\omega\|^2=\sum_{i=1}^d \sum_{j=1}^d \alpha_i \alpha_j e_i e_j = \sum_{i=1}^n \alpha_i^2 =1\]

于是可以得出结论

\[Var(Y) \le \lambda_1\]

并且当 \(\omega = e_1\) 时,等号成立,对应的映射值就为

\[Y = e_1 \cdot X\]

这里的证明就完成了,下面我们换一种思路,采用更直接的方式来计算最大方差,即

\[\max_{\omega} \quad Var(Y)\\ s.t.\quad \|\omega\| = 1\]

由于前面的推导,我们已知

\[Var(Y) = \omega^T \Sigma \omega\]

于是利用拉格朗日乘子法,建立拉格朗日函数

\[L(\omega, \lambda) = \omega^T \Sigma \omega - \lambda(\omega^T \omega -1)\]

对 \(\omega\) 求梯度得到

\[\Delta_{\omega} L = \omega^T \Sigma - \lambda \omega^T\]

再令 \(\Delta_{\omega} =0\) ,就有

\[\Sigma \omega = \lambda \omega\]

可以看到,\(\lambda\) 就是协方差矩阵的特征值,将上式代入 \(Var(Y)\) ,便得到

\[Var(Y) = \lambda \omega^T \omega = \lambda\]

于是我们再一次证明了 \(Var(Y)\) 的最大值就等于协方差矩阵的最大特征值,相对应的 \(\omega\) 就是协方差矩阵的最大特征值对应的特征向量。这里的 \(Y\) 就是对 \(X\) 的一种表征,它只有一个维度,并且与最大的特征值相关,所以称之为 \(X\) 的第一主成分,可用 \(Y_1\) 表示。但是显然,虽然第一主成分能够一定程度上代表原始特征,却失去了原始特征的很多信息,为了尽量补全信息,下面我们将考虑更多的次要成分。

次要成分与数据重建

在我们计算第一主成分的时候,考虑的是用一条直线(在高维情况下为超平面)作为基线,然后将特征点投影上去,用距离作为特征点的代表,并且这条直线使得距离有很好的散布性。同样的道理,计算第二个主成分,也是要找到一条直线或者超平面来做投影,但是它必须要和计算第一主成分时的基线正交,否则得到的第二主成分将有一部分含有第一主成分的信息。于是计算第二主成分的系数向量 \(\omega^{(2)}\) 必须满足条件

\[\omega^{(2)} \cdot e_1 = 0\]

以及当然

\[\|\omega^{(2)}\| = 1\]

接下来我们将证明,满足上述条件的 \(\omega^{(2)}\) 必然使 \(Y_2 = \omega^{(2)} \cdot X\) 的最大值为 \(\Sigma\) 第二大的特征值 \(\lambda_2\),并且 \(\omega^{(2)} = e_2\)。首先对 \(\omega^{(2)}\) 进行分解

\[\omega^{(2)} = \sum_{i=1}^d \alpha_i^{(2)} e_i\]

考虑到 \(\omega^{(2)}\) 与 \(e_1\) 正交,所以

\[\omega^{(2)} \cdot e_1 = 0\quad \Rightarrow \quad\sum_{i=1}^d \alpha_i^{(2)} e_i \cdot e_1 = 0 \quad \Rightarrow \quad \alpha_1^{(2)} = 0\]

也就是说,第一个方向上的分量为 0。然后根据方差的性质

\[\begin{aligned} Var(Y_2) &= \sum_{i=1}^d \sum_{j=1}^d \alpha_i^{(2)} \alpha_j^{(2)} e_i^T \Sigma e_j \\&= \sum_{i=1}^d \sum_{j=1}^d \alpha_i^{(2)} \alpha_j^{(2)} e_i^T \lambda_j e_j\\ &= \sum_{i=1}^d \alpha_i^{(2)} \alpha_i^{(2)} \lambda_i\\ &= \sum_{i=2}^d \alpha_i^{(2)} \alpha_i^{(2)} \lambda_i\\ &\le \lambda_2 \sum_{i=2}^d \alpha_i^{(2)} \alpha_i^{(2)} \end{aligned}\]

同样,根据 \(\|\omega^{(2)}\|=1\),可得

\[\sum_{i=1}^d \alpha_i^{(2)} \alpha_i^{(2)} =1\]

于是有

\[Var(Y_2) \le \lambda_2\]

其中取等号的条件是,不等式中 \(\alpha_2^{(2)} = 1\), 其余系数为 0。于是我们就证明了,当 \(\omega_{(2)} = e_2\) 时,\(Var(Y_2)\) 取得最大值 \(\lambda_2\)。

基于同样的方法,我们还能寻找第三主成分,第四主成分等等。将这些主成分 \(Y_1, Y_2, Y_3 ...\) 按列合并,形成新的数据矩阵

\[Y = [Y_1 \quad Y_2 \quad ...\quad Y_m] \quad m \le d\]

即为原特征矩阵 \(X\) 的降维变换。

最小重建误差

我们逐步寻找原始数据的第一、第二、第三、第四等等主成分的过程,实际上是在逐步重建数据,比如原始的数据由矩阵 \(X\) 表示

\[X = [X_1 \quad X_2 \quad ...\quad X_d]\]

它的每一行表示一个数据点,即

\[x_i = [X_{1i} \quad X_{2i} \quad \cdots \quad X_{di}]\]

其中 \(X_{ji}\) 表示 \(X_j\) 的第 i 个分量,而重建之后的数据

\[Y = [Y_1 \quad Y_2 \quad ... Y_m] \quad m \le d\]

它的每一行表示一个重建之后的数据点,我们这里用 \(\hat x_i\) 表示

\[\hat x_i = [Y_{1i}\quad Y_{2i} \quad \cdots \quad Y_{mi}]\]

其中 \(Y_{ji}\) 表示 \(Y_j\) 的第 i 个分量。

然后我们用主成分分析方法跑一个实际的例子,使用著名的 鸢尾花数据集 。这个数据集包含 3 种鸢尾花的子属,记录了花瓣与花萼的长宽尺寸作为特征。

(单独使用一种特征绘制的散点图,不同颜色代表不同种类)

分别计算这四个特征的方差值为

\[\begin{aligned} Var(SL) &= 0.685694\\ Var(SW) &= 0.188004\\ Var(PL) &= 3.11318\\ Var(PW) &= 0.582414 \end{aligned}\]

这里使用 \(SL, SW, PL,PW\) 分别表示 sepal length, sepal width, petal length, petal width。可以发现,虽然从图上看,petal width 的区分度要明显高于 sepal length,但是它们的方差却是反过来的。细想一下,造成这一现象的原因在于我们没有统一两个特征的度量,这样的话,数值上较大的特征计算的方差肯定会大一点,于是将这些特征归一化处理之后(也就是对每类特征,减去平均值,再除以最大值),再来看看方差

\[\begin{aligned} Var(SL) &= 0.162107\\ Var(SW) &= 0.103771\\ Var(PL) &= 0.315483\\ Var(PW) &= 0.343918 \end{aligned}\]

这就要显得合理许多。

(数据归一化之后的特征分布)

然后将归一化的特征代入主成分分析算法,得到的四阶主成分

它们的方差(或者说协方差矩阵的特征值)分别为

\[\begin{aligned} Var(Y_1) &= 0.785774\\ Var(Y_2) &= 0.10246\\ Var(Y_3) &= 0.0307906\\ Var(Y_4) &= 0.00600804 \end{aligned}\]

可以看到其中第一主成分占据了相当大的一部分,其他成分显得更加次要。通过主成分分析,我们找到了训练特征中最具差异化的特征表达,它或许不具备对应的现实意义,例如花瓣的宽度,只是这些特征的加权组合,却能够更好地表达特征。

本文遵守 CC-BY-NC-4.0 许可协议。

Creative Commons License

欢迎转载,转载需注明出处,且禁止用于商业目的。