Principal Component Analysis
主成分分析
主成分分析 (Principal Component Analysis, PCA) 是一种经典的数据降维方法。
推导：最小化近似误差
设有数据集 \(\{\mathbf x_n\}_{n=1}^N\)，其中每个数据点 \(\mathbf x_n\in\mathbb R^d\). 不失一般性地，假设数据的均值为 \(0\)，否则可以事先进行中心化。主成分分析希望找到 \(r\;(r<d)\) 个标准正交向量 \(\mathbf w_1,\mathbf w_2,\ldots,\mathbf w_r\in\mathbb R^d\) 作为基底，用它们的线性组合近似表示原数据： \[ \mathbf x_n\xleftarrow{\text{approximate}} \hat{\mathbf x}_n=\sum_{i=1}^rz_n^i\mathbf w_i=\mathbf W\mathbf z_n,\quad n=1,2,\ldots,N \] 其中 \(\mathbf W=[\mathbf w_1,\mathbf w_2,\cdots,\mathbf w_r]\in\mathbb R^{d\times r}\)，满足 \(\mathbf W^\mathsf T\mathbf W=\mathbf I_r\)；\(\mathbf z_n=(z_n^1,z_n^2,\ldots,z_n^r)^\mathsf T\in\mathbb R^r\) 是表示坐标，也就是数据 \(\mathbf x_n\) 降维后的结果。换句话说，我们将数据点 \(\mathbf x_n\) 从 \(\mathbb R^d\) 映射到了由 \(\mathbf w_1,\mathbf w_2,\ldots,\mathbf w_r\) 张成的 \(r\) 维子空间中。
为了找到最优的基底和表示坐标，我们最小化近似误差： \[ \mathcal L(\mathbf W,\{\mathbf z_n\})=\frac{1}{N}\sum_{n=1}^N\left\Vert\mathbf x_n-\mathbf W\mathbf z_n\right\Vert^2 \] 首先求解各个 \(\mathbf z_n\). 上式对 \(\mathbf z_n\) 求导并令为零，有： \[ \nabla_{\mathbf z_n}\mathcal L(\mathbf W,\{\mathbf z_n\})=-\frac{2}{N}\sum_{n=1}^N\mathbf W^\mathsf T(\mathbf x_n-\mathbf W\mathbf z_n)=-\frac{2}{N}\sum_{n=1}^N\left(\mathbf W^\mathsf T\mathbf x_n-\mathbf z_n\right)=\mathbf 0 \] 解得： \[ \mathbf z_n=\mathbf W^\mathsf T\mathbf x_n \] 代回近似误差并化简，得： \[ \begin{align} \mathcal L(\mathbf W)&=\frac{1}{N}\sum_{n=1}^N\left\Vert\mathbf x_n-\mathbf W\mathbf W^\mathsf T\mathbf x_n\right\Vert^2\\ &=\frac{1}{N}\sum_{n=1}^N\mathbf x_n^\mathsf T(\mathbf I_d-\mathbf W\mathbf W^\mathsf T)(\mathbf I_d-\mathbf W\mathbf W^\mathsf T)\mathbf x_n\\ &=\frac{1}{N}\sum_{n=1}^N\mathbf x_n^\mathsf T(\mathbf I_d-\mathbf W\mathbf W^\mathsf T)\mathbf x_n\\ &=\frac{1}{N}\sum_{n=1}^N\mathbf x_n^\mathsf T\mathbf x_n-\frac{1}{N}\sum_{n=1}^N\mathbf x_n^\mathsf T\mathbf W\mathbf W^\mathsf T\mathbf x_n\\ &=\mathrm{tr}(\mathbf S)-\mathrm{tr}(\mathbf W^\mathsf T\mathbf S\mathbf W) \end{align} \] 这里我们引入了记号： \[ \mathbf S=\frac{1}{N}\sum_{n=1}^N\mathbf x_n\mathbf x_n^\mathsf T \] 表示样本协方差矩阵。由于第一项与 \(\mathbf W\) 无关，因此现在优化问题化作： \[ \begin{align} \max_{\mathbf W}&~~\mathrm{tr}(\mathbf W^\mathsf T\mathbf S\mathbf W)\\ \text{s.t.}&~~\mathbf W^\mathsf T\mathbf W=\mathbf I_r \end{align} \] 这是经典的 Rayleigh 商问题的矩阵版本。引入拉格朗日乘子 \(\mathbf H\in\mathbb R^{r\times r}\)，构造拉格朗日函数： \[ L(\mathbf W,\mathbf H)=\mathrm{tr}(\mathbf W^\mathsf T\mathbf S\mathbf W)-\mathrm{tr}(\mathbf H^\mathsf T(\mathbf W^\mathsf T\mathbf W-\mathbf I_r)) \] 对 \(\mathbf W\) 求导并令为零，得： \[ \begin{gather} \nabla_{\small\mathbf W}L(\mathbf W,\mathbf H)=2\mathbf S\mathbf W-\mathbf W\mathbf H^\mathsf T\mathbf H=\mathbf 0_{d\times r}\\ \mathbf S\mathbf W=\frac{1}{2}\mathbf W\mathbf H^\mathsf T\mathbf H \end{gather} \] 这看着很像实对称矩阵 \(\mathbf S\) 的谱分解。为了看得更清楚，注意到 \(\frac{1}{2}\mathbf H^\mathsf T\mathbf H\) 是一个实对称矩阵，因此可分解为 \(\frac{1}{2}\mathbf H^\mathsf T\mathbf H=\mathbf P\mathbf\Lambda\mathbf P^\mathsf T\)，其中 \(\mathbf P\) 是正交矩阵，\(\mathbf\Lambda=\mathrm{diag}(\lambda_1,\lambda_2,\ldots,\lambda_r)\) 是对角阵。代入得： \[ \begin{gather} \mathbf S\mathbf W=\mathbf W\mathbf P\mathbf\Lambda\mathbf P^\mathsf T\\ \mathbf S(\mathbf W\mathbf P)=(\mathbf W\mathbf P)\mathbf\Lambda \end{gather} \] 这说明，\(\lambda_1,\lambda_2,\ldots,\lambda_r\) 是 \(\mathbf S\) 的 \(r\) 个特征值，\(\mathbf W\mathbf P\) 各列是对应的 \(r\) 个特征向量。但是 \(\mathbf S\) 一共有 \(d\) 个特征值，为了确定究竟选取哪 \(r\) 个特征值，考察此时的优化目标： \[ \mathrm{tr}(\mathbf W^\mathsf T\mathbf S\mathbf W)=\mathrm{tr}\left((\mathbf W\mathbf P)^\mathsf T\mathbf S(\mathbf W\mathbf P)\right)=\mathrm{tr}\left((\mathbf W\mathbf P)^\mathsf T(\mathbf W\mathbf P)\mathbf\Lambda\right)=\mathrm{tr}(\mathbf\Lambda)=\sum_{i=1}^r\lambda_i \] 可见为了最大化优化目标，\(\lambda_1,\lambda_2,\ldots,\lambda_r\) 应取为 \(\mathbf S\) 最大的 \(r\) 个特征值。又鉴于正交矩阵 \(\mathbf P\) 的任意性，可知 \(\mathbf w_1,\mathbf w_2,\ldots,\mathbf w_r\) 张成的空间就是 \(\lambda_1,\lambda_2,\ldots,\lambda_r\) 对应特征向量所张成的空间，该空间中的任意一组标准正交基都可以是我们要找的解。为了唯一确定一个解，我们可以引入额外的要求。注意到对 \(r=1,2,\ldots\)，解空间是从小到大依次嵌套的，因此我们可以要求当前解应该保留降到更低维度时的解。具体而言，当 \(r=1\) 时，有 \(\mathbf w_1\) 是 \(\lambda_1\) 对应特征向量；当 \(r=2\) 时，保留 \(\mathbf w_1\)，则 \(\mathbf w_2\) 必须与 \(\mathbf w_1\) 正交，由于实对称矩阵的特征向量都是正交的，所以 \(\mathbf w_2\) 正是 \(\lambda_2\) 对应特征向量……以此类推，最终必然有 \(\mathbf w_1,\mathbf w_2,\ldots,\mathbf w_r\) 分别是 \(\lambda_1,\lambda_2,\ldots,\lambda_r\) 对应的特征向量。
总结一下，PCA 的执行过程为：
- 中心化数据 \(\{\mathbf x_n\}_{n=1}^N\)，计算样本协方差矩阵 \(\mathbf S\)；
- 求解 \(\mathbf S\) 最大的 \(r\) 个特征值所对应的特征向量，排列为投影矩阵 \(\mathbf W=[\mathbf w_1,\mathbf w_2,\ldots,\mathbf w_r]\)；
- 执行降维 \(\mathbf z_n=\mathbf W^\mathsf T\mathbf x_n\)；
- 若要从降维结果重构原数据，执行 \(\hat{\mathbf x}_n=\mathbf W\mathbf z_n\).
解释：最大化投影方差
PCA 寻找到的子空间除了使得近似误差最小以外，还可以解释为使得数据投影后的方差最大。为了看清这一点，考察降维投影后的数据 \(\mathbf z_n=\mathbf W^\mathsf T\mathbf x_n\)，由于我们假设了原数据的均值为零，因此易知投影后数据的均值也是零，故其协方差矩阵为： \[ \mathbf\Sigma=\frac{1}{N}\sum_{n=1}^N\mathbf z_n\mathbf z_n^\mathsf T=\frac{1}{N}\sum_{n=1}^N\mathbf W^\mathsf T\mathbf x_n\mathbf x_n^\mathsf T\mathbf W=\mathbf W^\mathsf T\mathbf S\mathbf W \] 值得注意的是 \(\mathbf\Sigma\) 其实是一个对角阵，因为： \[ \mathbf\Sigma_{ij}=\mathbf w_i^\mathsf T\mathbf S\mathbf w_j=\lambda_j\mathbf w_i^\mathsf T\mathbf w_j=\begin{cases}\lambda_i,&\text{if}~~i=j\\0,&\text{otherwise}\end{cases} \] 也就是说投影后数据的各个维度是不相关（没有线性关系）的。
回顾上文推导的近似误差表达式： \[ \mathcal L(\mathbf W)=\mathrm{tr}(\mathbf S)-\mathrm{tr}(\mathbf W^\mathsf T\mathbf S\mathbf W) \] 也即： \[ \mathrm{tr}(\mathbf S)=\mathcal L(\mathbf W)+\mathrm{tr}(\mathbf\Sigma) \] 可见原数据各维方差之和等于近似误差加上投影数据各维方差之和。由于原数据方差是定值，因此最小化近似误差与最大化投影数据方差是等价的。几何直观上，考察单个数据点，则原数据方差体现为数据点到原点的距离，近似误差就是投影垂线的长度，而投影数据方差可看作投影垂足到原点的距离，因此上式描述的其实就是勾股定理。显然，当斜边（原数据方差）固定时，两条直角边（近似误差、投影数据方差）呈负相关关系。
实践：三点计算考量
通过前面的推导，我们知道 PCA 就是要求解 \(\mathbf X^\mathsf T\mathbf X\) 前 \(r\) 大的特征值对应的特征向量，其中 \(\mathbf X=(\mathbf x_1^\mathsf T,\mathbf x_2^\mathsf T,\ldots,\mathbf x_N^\mathsf T)^\mathsf T\in\mathbb R^{N\times d}\) 表示样本构成的矩阵。不过在实践中，我们还有几点计算上的考量。
首先，我们最好保证数据在各个维度上的取值范围差不多，否则 PCA 的解将偏向于取值范围大的维度，因为它们具有更大的方差。因此实践中我们最好对数据进行标准化，而非仅仅中心化。
其次，在高维场景下，即 \(d>N\) 时，为了节省计算量，我们可以求解 Gram 矩阵 \(\mathbf X\mathbf X^\mathsf T\in\mathbb R^{N\times N}\) 而非 \(\mathbf X^\mathsf T\mathbf X\in\mathbb R^{d\times d}\) 的特征值与特征向量。这是因为二者有相同的非零特征值，并且特征向量也有简单的变换关系。具体而言，设有： \[ \mathbf X\mathbf X^\mathsf T\mathbf U=\mathbf U\boldsymbol\Lambda \] 其中 \(\mathbf U\) 是特征向量构成的矩阵，\(\boldsymbol\Lambda\) 是特征值构成的对角阵。那么左乘 \(\mathbf X^\mathsf T\) 得： \[ \mathbf X^\mathsf T\mathbf X(\mathbf X^\mathsf T\mathbf U)=(\mathbf X^\mathsf T\mathbf U)\boldsymbol\Lambda \] 这说明 \(\mathbf X^\mathsf T\mathbf X\) 的特征值也是 \(\boldsymbol\Lambda\)，而特征向量构成的矩阵为 \(\mathbf X^\mathsf T\mathbf U\). 不过这些特征向量尚未归一化，因为： \[ (\mathbf X^\mathsf T\mathbf U)^\mathsf T(\mathbf X^\mathsf T\mathbf U)=\mathbf U^\mathsf T\mathbf X\mathbf X^\mathsf T\mathbf U=\mathbf U^\mathsf T\mathbf U\mathbf\Lambda=\mathbf\Lambda \] 因此归一化后的特征向量矩阵应为： \[ \mathbf V=\mathbf X^\mathsf T\mathbf U\mathbf\Lambda^{-1/2} \] 我们将会看到，这种计算方式也让我们能够在 PCA 中用上核技巧。
最后，根据特征值与奇异值的关系，求解 \(\mathbf X^\mathsf T\mathbf X\) 的特征值与特征向量，等价于求解 \(\mathbf X\) 的奇异值与右奇异向量。这是因为，设有奇异值分解 \(\mathbf X=\mathbf U\boldsymbol\Sigma\mathbf V^\mathsf T\)，那么： \[ \mathbf X^\mathsf T\mathbf X=(\mathbf U\boldsymbol\Sigma\mathbf V^\mathsf T)^\mathsf T(\mathbf U\boldsymbol\Sigma\mathbf V^\mathsf T)=\mathbf V\boldsymbol\Sigma^2\mathbf V^\mathsf T \] 所以我们也可直接求解 \(\mathbf X\) 最大的 \(r\) 个奇异值对应的右奇异向量。
补充：线性自编码器
如果我们将降维过程视作编码器，重构过程视作解码器，那么先降维、再重构的过程就构成了一个自编码器 (autoencoder)： \[ \min_{\text{encode},\,\text{decode}}~\mathcal L=\frac{1}{N}\sum_{n=1}^N\left\Vert\mathbf x_n-\text{decode}(\text{encode}(\mathbf x_n))\right\Vert^2 \] 对比 PCA 的优化目标： \[ \min_{\mathbf W}~\mathcal L=\frac{1}{N}\sum_{n=1}^N\left\Vert\mathbf x_n-\mathbf W\mathbf W^\mathsf T\mathbf x_n\right\Vert^2 \] 可见 PCA 就是在「编解码器均为正交变换且互为转置」的限制下的自编码器。
事实上，即便抛弃正交和转置的限制，我们仍然能够证明 PCA 与线性自编码器是等价的。具体而言，取： \[ \begin{gather} \text{encode}(\mathbf x;\mathbf W_\text{enc})=\mathbf W_\text{enc}\mathbf x\in\mathbb R^r\\ \text{decode}(\mathbf z;\mathbf W_\text{dec})=\mathbf W_\text{dec}\mathbf z\in\mathbb R^d \end{gather} \] 其中 \(\mathbf W_\text{enc}\in\mathbb R^{r\times d},\;\mathbf W_\text{dec}\in\mathbb R^{d\times r}\). 那么重构损失写作： \[ \mathcal L(\mathbf W_\text{enc},\mathbf W_\text{dec})=\frac{1}{N}\sum_{n=1}^N\left\Vert\mathbf x_n-\mathbf W_\text{dec}\mathbf W_\text{enc}\mathbf x_n\right\Vert^2=\frac{1}{N}\left\Vert\mathbf X^\mathsf T-\mathbf W_\text{dec}\mathbf W_\text{enc}\mathbf X^\mathsf T\right\Vert_F^2 \] 这是一个经典的最佳低秩逼近问题。根据 Eckart–Young–Mirsky 定理，设 \(\mathbf X\) 的奇异值分解为 \(\mathbf X=\mathbf U\mathbf\Sigma\mathbf V^\mathsf T\)，其中 \(\mathbf U\in\mathbb R^{N\times k},\mathbf\Sigma\in\mathbb R^{k\times k},\mathbf V\in\mathbb R^{d\times k}\)，则其最优解应满足： \[ \mathbf W_\text{dec}\mathbf W_\text{enc}\mathbf X^\mathsf T=\mathbf V_r\mathbf\Sigma_r\mathbf U_r^\mathsf T \] 其中 \(\mathbf U_r\in\mathbb R^{N\times r},\mathbf\Sigma_r\in\mathbb R^{r\times r},\mathbf V_r\in\mathbb R^{d\times r}\). 容易验证符合条件的一族解为： \[ \mathbf W_\text{dec}=\mathbf V_r\mathbf R,\;\mathbf W_\text{enc}=\mathbf R^{-1}\mathbf V_r^\mathsf T \] 其中 \(\mathbf R\) 是任意可逆矩阵。这说明求解线性自编码器得到的基与 PCA 解出的基仅差一个可逆变换，二者张成的子空间是相同的。
核主成分分析
在 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)^\mathsf T\phi(\mathbf x_m)\). 考虑在特征空间中实施 PCA，那么这就相当于在数据空间中做了非线性的主成分分析。
我们暂且假设映射后特征的均值依旧是 \(0\)，则特征空间中的样本协方差矩阵为： \[ \mathbf C=\frac{1}{N}\sum_{n=1}^N\phi(\mathbf x_n)\phi(\mathbf x_n)^\mathsf T=\frac{1}{N}\begin{bmatrix}\phi(\mathbf x_1)&\cdots&\phi(\mathbf x_N)\end{bmatrix}\begin{bmatrix}\phi(\mathbf x_1)^\mathsf T\\\vdots\\\phi(\mathbf x_N)^\mathsf T\end{bmatrix}\in\mathbb R^{D\times D} \] 不过，为了使用上核技巧，我们需要的是 \(\phi(\mathbf x_n)^\mathsf T\phi(\mathbf x_n)\) 而非 \(\phi(\mathbf x_n)\phi(\mathbf x_n)^\mathsf T\). 好在根据前文的讨论，求解 \(\mathbf C\) 的特征值与特征向量，可以转化为求解 Gram 矩阵的特征值与特征向量： \[ \begin{align} \mathbf K&=\begin{bmatrix}\phi(\mathbf x_1)^\mathsf T\\\vdots\\\phi(\mathbf x_N)^\mathsf T\end{bmatrix}\begin{bmatrix}\phi(\mathbf x_1)&\cdots&\phi(\mathbf x_N)\end{bmatrix}\\\\ &=\begin{bmatrix} \phi(\mathbf x_1)^\mathsf T\phi(\mathbf x_1)&\cdots&\phi(\mathbf x_1)^\mathsf T\phi(\mathbf x_N)\\ \vdots&\ddots&\vdots\\ \phi(\mathbf x_N)^\mathsf T\phi(\mathbf x_1)&\cdots&\phi(\mathbf x_N)^\mathsf 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,\mathbf u_2,\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^\mathsf 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^\mathsf T\mathbf K\mathbf u_i \] 因此数据点 \(\mathbf x\) 在 \(\mathbf v_i\) 上的投影为： \[ y_i(\mathbf x)=\frac{\mathbf v_i^\mathsf T\phi(\mathbf x)}{\mathbf v_i^\mathsf T\mathbf v_i}=\frac{1}{\mathbf u_i^\mathsf T\mathbf K\mathbf u_i}\sum_{n=1}^Nu_{in}k(\mathbf x_n,\mathbf x) \] 于是 \(\mathbf y=(y_1(\mathbf x),y_2(\mathbf x),\ldots,y_r(\mathbf x))^\mathsf T\) 就是 \(\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)^\mathsf T\tilde\phi(\mathbf x_m)\\ =&~\phi(\mathbf x_n)^\mathsf T\phi(\mathbf x_m)-\frac{1}{N}\sum_{l=1}^N\phi(\mathbf x_n)^\mathsf T\phi(\mathbf x_l)-\frac{1}{N}\sum_{l=1}^N\phi(\mathbf x_l)^\mathsf T\phi(\mathbf x_m)+\frac{1}{N^2}\sum_{j=1}^N\sum_{l=1}^N\phi(\mathbf x_j)^\mathsf 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\in\mathbb R^{N\times N}\) 表示所有元素均为 \(1/N\) 的矩阵。现在只需将上文的 \(\mathbf K\) 替换为 \(\tilde{\mathbf K}\) 进行计算即可。
参考资料
- Murphy, Kevin P. Probabilistic Machine Learning: An Introduction. ↩︎
- Bishop, Christopher M., and Nasser M. Nasrabadi. Pattern recognition and machine learning. ↩︎
- Bourlard, Hervé, and Yves Kamp. Auto-association by multilayer perceptrons and singular value decomposition. Biological cybernetics 59.4 (1988): 291-294. ↩︎
- LynnHo. Matrix-Calculus-Tutorial. https://github.com/LynnHo/Matrix-Calculus-Tutorial/blob/master/examples.md ↩︎