Principal Component Analysis
主成分分析
主成分分析 (PCA) 是一种降维方法，可以从两种角度推导——最大化投影方差与最小化重构误差。
首先引入记号。设数据集为 \(\{\mathbf x_n\}_{n=1}^N\)，其中 \(\mathbf x_n\in\mathbb R^d\)，也可以记作矩阵形式： \[ \mathbf X=\begin{bmatrix}\mathbf x_1^T\\\vdots\\\mathbf x_N^T\end{bmatrix}\in\mathbb R^{N\times d} \] 不失一般性地，我们假设数据的均值为 \(0\)，否则可以事先进行中心化。于是有样本协方差矩阵： \[ \mathbf S=\frac{1}{N}\sum_{n=1}^N\mathbf x_n\mathbf x_n^T=\frac{1}{N}\mathbf X^T\mathbf X\in\mathbb R^{d\times d} \] 显然 \(\mathbf S\) 是一个对称阵，即 \(\mathbf S^T=\mathbf S\). 我们将看到，PCA 的本质其实就是求解 \(\mathbf S\) 的特征值与特征向量。
角度一·最大化投影方差
为了降维，我们考虑将数据点投影到一个 \(r\) 维子空间中，其中 \(r<d\). 为了尽可能少地损失数据的信息，我们希望最大化投影后数据在各个维度上的方差。
首先考虑一维情形。设投影空间的基向量为 \(\mathbf u_1\in\mathbb R^d\)，满足 \(\mathbf u_1^T\mathbf u_1=1\)，则每个 \(\mathbf x_n\) 将被投影为 \(\mathbf u_1^T\mathbf x_n\). 投影后数据的样本均值仍为 \(0\)，而样本方差为： \[ \frac{1}{N}\sum_{n=1}^N\left(\mathbf u_1^T\mathbf x_n\right)^2=\mathbf u_1^T\mathbf S\mathbf u_1 \] 于是优化问题为： \[ \begin{align} \max_{\mathbf u_1}&~~\mathbf u_1^T\mathbf S\mathbf u_1\\ \text{s.t.}&~~\mathbf u_1^T\mathbf u_1=1 \end{align} \] 引入拉格朗日乘子 \(\lambda_1\)，构建拉格朗日函数： \[ L(\mathbf u_1,\lambda_1)=\mathbf u_1^T\mathbf S\mathbf u_1+\lambda_1(1-\mathbf u_1^T\mathbf u_1) \] 对 \(\mathbf u_1\) 求导并令为零，可得： \[ \frac{\partial L}{\partial \mathbf u_1}=2\mathbf S\mathbf u_1-2\lambda_1\mathbf u_1=\mathbf 0\implies \mathbf S\mathbf u_1=\lambda_1\mathbf u_1 \] 这说明 \(\lambda_1,\mathbf u_1\) 正是 \(\mathbf S\) 的特征值和特征向量。进一步地，为了使得优化目标 \(\mathbf u_1^T\mathbf S\mathbf u_1=\lambda_1\mathbf u_1^T\mathbf u_1=\lambda_1\) 最大，\(\lambda_1,\mathbf u_1\) 应该取 \(\mathbf S\) 的最大特征值及其对应的特征向量。我们称 \(\mathbf u_1\) 为第一个主成分。
对于更多维的情形，我们以增量的方式定义下一个主成分，并要求新的主成分与已求得的主成分正交。例如，对于第二个主成分 \(\mathbf u_2\)，优化问题写作： \[ \begin{align} \max_{\mathbf u_2}&~~\mathbf u_2^T\mathbf S\mathbf u_2\\ \text{s.t.}&~~\mathbf u_1^T\mathbf u_2=0,\;\mathbf u_2^T\mathbf u_2=1 \end{align} \] 引入拉格朗日乘子 \(\lambda_2,\lambda_{12}\)，构建拉格朗日函数： \[ L(\mathbf u_2,\lambda_2,\lambda_{12})=\mathbf u_2^T\mathbf S\mathbf u_2+\lambda_2(1-\mathbf u_2^T\mathbf u_2)-\lambda_{12}\mathbf u_1^T\mathbf u_2 \] 对 \(\mathbf u_2\) 求导并令为零，可得： \[ \frac{\partial L}{\partial\mathbf u_2}=2\mathbf S\mathbf u_2-2\lambda_2\mathbf u_2-\lambda_{12}\mathbf u_1=0 \] 上式左乘 \(\mathbf u_1^T\) 得： \[ 2\mathbf u_1^T\mathbf S\mathbf u_2-2\lambda_2\mathbf u_1^T\mathbf u_2-\lambda_{12}\mathbf u_1^T\mathbf u_1=0-0-\lambda_{12}=0 \] 其中第一项 \(\mathbf u_1^T\mathbf S\mathbf u_2=(\mathbf S\mathbf u_1)^T\mathbf u_2=\lambda_1\mathbf u_1^T\mathbf u_2=0\). 将 \(\lambda_{12}=0\) 代回，我们得到： \[ \mathbf S\mathbf u_2=\lambda_2\mathbf u_2 \] 这说明 \(\lambda_2,\mathbf u_2\) 也是 \(\mathbf S\) 的特征值和特征向量。于是，为了让优化目标最大，\(\lambda_2,\mathbf u_2\) 应该是 \(\mathbf S\) 的次大特征值及其对应特征向量。
以此类推，我们可以得到结论：使得投影方差最大化的 \(r\) 维子空间的基就是 \(\mathbf S\) 最大的 \(r\) 个特征值对应的特征向量。
角度二·最小化重构误差
第二个推导角度的出发点是在一个 \(r\;(r<d)\) 维的线性子空间中近似表示数据，并希望近似的误差最小。具体而言，设 \(r\) 维子空间的一组正交基为 \((\mathbf u_1,\ldots,\mathbf u_r)\)，其中 \(\mathbf u_i\in\mathbb R^d\)，则该子空间中的任一点可唯一表示为： \[ \tilde{\mathbf x}_n=\sum_{i=1}^r\gamma_{ni}\mathbf u_i \] 其中 \(\gamma_{n1},\ldots,\gamma_{nr}\) 是表示系数。将这组基扩充为 \(\mathbb R^d\) 中的完备正交基 \((\mathbf u_1,\ldots,\mathbf u_d)\)，那么每个 \(\mathbf x_n\) 在这组基下都有唯一的线性表示： \[ \mathbf x_n=\sum_{i=1}^d\alpha_{ni}\mathbf u_i=\sum_{i=1}^d(\mathbf x_n^T\mathbf u_i)\mathbf u_i \] 其中 \(\alpha_{ni}=\mathbf x_n^T\mathbf u_i\) 是表示系数。于是，重构误差写作： \[ J=\frac{1}{N}\sum_{n=1}^N\Vert\mathbf x_n-\tilde{\mathbf x}_n\Vert^2=\frac{1}{N}\sum_{n=1}^N\left\Vert\sum_{i=1}^r(\mathbf x_n^T\mathbf u_i-\gamma_{ni})\mathbf u_i+\sum_{i=r+1}^d(\mathbf x_n^T\mathbf u_i)\mathbf u_i\right\Vert^2 \] 对 \(\gamma_{ni}\) 求导并令为零，解得： \[ \gamma_{ni}=\mathbf x_n^T\mathbf u_i \] 于是优化目标简化为： \[ \begin{align} J&=\frac{1}{N}\sum_{n=1}^N\left\Vert\sum_{i=r+1}^d\left(\mathbf x_n^T\mathbf u_i\right)\mathbf u_i\right\Vert^2\\ &=\frac{1}{N}\sum_{n=1}^N\left[\sum_{i=r+1}^d\left(\mathbf x_n^T\mathbf u_i\right)\mathbf u_i\right]^T\left[\sum_{j=r+1}^d\left(\mathbf x_n^T\mathbf u_j\right)\mathbf u_j\right]\\ &=\frac{1}{N}\sum_{n=1}^N\sum_{i=r+1}^d\left(\mathbf x_n^T\mathbf u_i\right)^2\\ &=\sum_{i=r+1}^d\mathbf u_i^T\mathbf S\mathbf u_i \end{align} \] 至此，我们的优化问题变成了： \[ \begin{align} \min_{\{\mathbf u_i\}}&~~\sum_{i=r+1}^d\mathbf u_i^T\mathbf S\mathbf u_i\\ \text{s.t.}&~~\mathbf u_i^T\mathbf u_i=1\\ &~~\mathbf u_i^T\mathbf u_j=0,\;i\neq j\\ \end{align} \] 这与上一节得到的优化问题是类似的，只不过这里是最小化，因此解就是 \(\mathbf S\) 最小的 \(d-r\) 个特征值对应的特征向量。相应地，\(\mathbf u_1,\ldots,\mathbf u_r\) 张成的空间就是 \(\mathbf S\) 最大的 \(r\) 个特征值对应的特征向量张成的空间，因此直接取 \(\mathbf u_1,\ldots,\mathbf u_r\) 为 \(\mathbf S\) 最大的 \(r\) 个特征值对应的特征向量即可。
计算上的考量
通过上面两小节的推导，我们发现 PCA 就是要求解样本协方差矩阵 \(\mathbf S\) 的前 \(r\) 大的特征值与对应的特征向量。然而众所周知，求解特征值与特征向量并不是一件容易的事，所以我们会有一些计算上的考量。
首先，我们最好保证数据在各个维度上的取值范围差不多，否则 PCA 的解将偏向于取值范围大的维度，因为它们具有更大的方差。因此实践中我们最好对数据进行标准化，而非仅仅中心化。等价地，这相当于我们用相关矩阵、而非协方差矩阵去求解。
其次，在高维度场景下，即 \(d>N\) 时，为了节省计算量，我们可以求解 Gram 矩阵 \(\mathbf X\mathbf X^T\in\mathbb R^{N\times N}\) 的特征值与特征向量而非协方差矩阵 \(\mathbf X^T\mathbf X\in\mathbb R^{d\times d}\) 的特征值与特征向量。这是因为矩阵论的知识告诉我们，\(\mathbf X\mathbf X^T\) 与 \(\mathbf X^T\mathbf X\) 有着相同的非零特征值，并且特征向量也有简单的变换关系。具体而言，假设 \(\mathbf X\mathbf X^T\) 的特征值分解为： \[ (\mathbf X\mathbf X^T)\mathbf U=\mathbf U\boldsymbol\Lambda \] 其中 \(\mathbf U\) 是特征向量构成的矩阵，\(\boldsymbol\Lambda\) 是特征值构成的对角阵。那么左乘 \(\mathbf X^T\) 得： \[ (\mathbf X^T\mathbf X)(\mathbf X^T\mathbf U)=(\mathbf X^T\mathbf U)\boldsymbol\Lambda \] 这说明 \(\mathbf X^T\mathbf X\) 的特征值也是 \(\boldsymbol\Lambda\)，而特征向量构成的矩阵为 \(\mathbf X^T\mathbf U\). 我们将会看到，这种计算方式也让我们能够在 PCA 中用上核技巧。
最后，根据特征值与奇异值的关系，求解 \(\mathbf X^T\mathbf X\) 的特征值与特征向量，其实等价于求解 \(\mathbf X\) 的奇异值与右奇异向量。具体而言，设有奇异值分解 \(\mathbf X=\mathbf U\boldsymbol\Sigma\mathbf V^T\)，那么： \[ \mathbf X^T\mathbf X=(\mathbf U\boldsymbol\Sigma\mathbf V^T)^T(\mathbf U\boldsymbol\Sigma\mathbf V^T)=\mathbf V\boldsymbol\Sigma^2\mathbf V^T \] 所以我们直接求解 \(\mathbf X\) 最大的 \(r\) 个奇异值对应的右奇异向量即可。
核主成分分析
在 SVM 和 RKHS 中，我们介绍了核技巧及其在 SVM 中的应用。自然地可以想到，核技巧也可以用在 PCA 中。具体而言，设有非线性映射 \(\phi\) 将数据点 \(\mathbf x_n\in\mathbb R^d\) 映射为特征 \(\phi(\mathbf x_n)\in\mathbb R^D\)，其中 \(D\) 可以是无穷。设核函数为 \(k(\cdot,\cdot)\)，满足： \[ k(\mathbf x_n,\mathbf x_m)=\phi(\mathbf x_n)^T\phi(\mathbf x_m) \] 考虑在特征空间中实施 PCA，那么这就相当于在数据空间中做了非线性的主成分分析。
我们暂时假设映射后特征的均值依旧是 \(0\)，则特征空间中的样本协方差矩阵为： \[ \mathbf C=\frac{1}{N}\sum_{n=1}^N\phi(\mathbf x_n)\phi(\mathbf x_n)^T=\frac{1}{N}\begin{bmatrix}\phi(\mathbf x_1)&\cdots&\phi(\mathbf x_N)\end{bmatrix}\begin{bmatrix}\phi(\mathbf x_1)^T\\\vdots\\\phi(\mathbf x_N)^T\end{bmatrix}\in\mathbb R^{D\times D} \] 不过，为了使用上核技巧，我们需要的是 \(\phi(\mathbf x_n)^T\phi(\mathbf x_n)\) 而非 \(\phi(\mathbf x_n)\phi(\mathbf x_n)^T\). 好在根据前文的讨论，求解 \(\mathbf C\) 的特征值与特征向量，可以转化为求解 Gram 矩阵的特征值与特征向量： \[ \begin{align} \mathbf K&=\begin{bmatrix}\phi(\mathbf x_1)^T\\\vdots\\\phi(\mathbf x_N)^T\end{bmatrix}\begin{bmatrix}\phi(\mathbf x_1)&\cdots&\phi(\mathbf x_N)\end{bmatrix}\\\\ &=\begin{bmatrix} \phi(\mathbf x_1)^T\phi(\mathbf x_1)&\cdots&\phi(\mathbf x_1)^T\phi(\mathbf x_N)\\ \vdots&\ddots&\vdots\\ \phi(\mathbf x_N)^T\phi(\mathbf x_1)&\cdots&\phi(\mathbf x_N)^T\phi(\mathbf x_N) \end{bmatrix}\\\\ &=\begin{bmatrix} k(\mathbf x_1,\mathbf x_1)&\cdots&k(\mathbf x_1,\mathbf x_N)\\ \vdots&\ddots&\vdots\\ k(\mathbf x_1,\mathbf x_1)&\cdots&k(\mathbf x_1,\mathbf x_N) \end{bmatrix}\in\mathbb R^{N\times N} \end{align} \] 假设求解出 \(\mathbf K\) 的特征向量为 \(\mathbf u_1,\ldots,\mathbf u_N\)，则 \(\mathbf C\) 的第 \(i\) 个特征向量为： \[ \mathbf v_i=\begin{bmatrix}\phi(\mathbf x_1)&\cdots&\phi(\mathbf x_N)\end{bmatrix}\mathbf u_i=\sum_{n=1}^Nu_{in}\phi(\mathbf x_n) \] 其中 \(u_{in}\) 表示 \(\mathbf u_i\) 的第 \(n\) 个维度。\(\mathbf v_i\) 模长的平方为： \[ \Vert\mathbf v_i\Vert^2=\mathbf v_i^T\mathbf v_i=\sum_{n=1}^N\sum_{m=1}^Nu_{in}u_{im}k(\mathbf x_n,\mathbf x_m)=\mathbf u_i^T\mathbf K\mathbf u_i \] 于是数据点 \(\mathbf x\) 在 \(\mathbf v_i\) 上的投影为： \[ y_i(\mathbf x)=\frac{\mathbf v_i^T\phi(\mathbf x)}{\mathbf v_i^T\mathbf v_i}=\frac{1}{\mathbf u_i^T\mathbf K\mathbf u_i}\sum_{n=1}^Nu_{in}k(\mathbf x_n,\mathbf x) \] 现在还有一个遗留问题——我们假设映射后特征的均值依旧是 \(0\)，但这并不总是正确（即便数据的均值是 \(0\)）。假设中心化后的特征为： \[ \tilde\phi(\mathbf x_n)=\phi(\mathbf x_n)-\frac{1}{N}\sum_{l=1}^N\phi(\mathbf x_l) \] 那么相应的 Gram 矩阵 \(\tilde{\mathbf K}\) 的第 \(n\) 行第 \(m\) 列为： \[ \begin{align} \tilde K_{nm}=&~\tilde\phi(\mathbf x_n)^T\tilde\phi(\mathbf x_m)\\ =&~\phi(\mathbf x_n)^T\phi(\mathbf x_m)-\frac{1}{N}\sum_{l=1}^N\phi(\mathbf x_n)^T\phi(\mathbf x_l)-\frac{1}{N}\sum_{l=1}^N\phi(\mathbf x_l)^T\phi(\mathbf x_m)+\frac{1}{N^2}\sum_{j=1}^N\sum_{l=1}^N\phi(\mathbf x_j)^T\phi(\mathbf x_l)\\ =&~k(\mathbf x_n,\mathbf x_m)-\frac{1}{N}\sum_{l=1}^Nk(\mathbf x_n,\mathbf x_l)-\frac{1}{N}\sum_{l=1}^Nk(\mathbf x_m,\mathbf x_l)+\frac{1}{N^2}\sum_{j=1}^N\sum_{l=1}^Nk(\mathbf x_j,\mathbf x_l) \end{align} \] 写作矩阵形式，即： \[ \tilde{\mathbf K}=\mathbf K-\mathbf 1_N\mathbf K-\mathbf K\mathbf 1_N+\mathbf 1_N\mathbf K\mathbf 1_N \] 其中 \(\mathbf 1_N\) 表示所有元素均为 \(1/N\) 的 \(N\times N\) 的矩阵。于是我们求解 \(\tilde{\mathbf K}\) 的特征值与特征向量即可。