PCA 与 Eigenfaces：从方差最大化到人脸识别

上一篇我们建立了矩阵微积分——梯度、Jacobian、Hessian——工具链的第三件利器。至此，Part 1 “拆”的四件核心工具（分解、度量、微分、迭代）的前三件已经齐备。现在是时候把它们用起来了。

真实世界的数据往往是高维的——一张 $112 \times 92$ 的灰度人脸有 10304 个像素，一段基因组测序有上百万维。在这样的高维空间中，存储、计算和可视化都变得不可行，而且大量维度之间高度相关，携带的独立信息远少于维度数暗示的那么多。如何找到数据中真正重要的少数几个方向，用它们忠实地表示整个数据集？

PCA（Principal Component Analysis，主成分分析） 给出了一个优雅的答案：找到方差最大的方向作为新坐标轴。核心洞察是——数据沿哪个方向”散得最开”，哪个方向就携带最多信息。如果数据点沿某个方向几乎不变，那个方向就是冗余的，可以丢弃。

数学上，PCA 完美地将前面学到的三件工具编织在一起：对数据的协方差矩阵做特征分解（Art. 2 特征分解），等价地对中心化数据矩阵做 SVD（Art. 3 SVD），而优化目标的推导用到了 Lagrange 乘子法（Art. 5 矩阵微积分 中矩阵微积分的直接应用）。

本文将从三个等价视角推导 PCA，然后展示它最著名的应用——Eigenfaces（特征脸）人脸识别，最后讨论 PCA 的局限。

三种等价视角

PCA 可以从三个不同的角度理解，它们在数学上严格等价。

视角一：最大方差

给定 $n$ 个 $d$ 维数据点 $\mathbf{x}_1, \ldots, \mathbf{x}_n \in \mathbb{R}^d$。中心化后的数据点为 $\tilde{\mathbf{x}}_i = \mathbf{x}_i - \bar{\mathbf{x}}$，其中 $\bar{\mathbf{x}} = \frac{1}{n}\sum_{i=1}^n \mathbf{x}_i$ 是均值。

问题：找一个单位方向 $\mathbf{w} \in \mathbb{R}^d$（$\|\mathbf{w}\| = 1$），使数据在该方向上的投影方差最大。

数据点 $\tilde{\mathbf{x}}_i$ 在方向 $\mathbf{w}$ 上的投影是标量 $\tilde{\mathbf{x}}_i^T \mathbf{w}$。投影的方差为：

$\text{Var} = \frac{1}{n-1}\sum_{i=1}^n (\tilde{\mathbf{x}}_i^T \mathbf{w})^2 = \mathbf{w}^T \underbrace{\left(\frac{1}{n-1}\sum_{i=1}^n \tilde{\mathbf{x}}_i \tilde{\mathbf{x}}_i^T\right)}_{C} \mathbf{w} = \mathbf{w}^T C \mathbf{w}$

其中 $C \in \mathbb{R}^{d \times d}$ 是协方差矩阵（covariance matrix）：

$C = \frac{1}{n-1}\tilde{X}^T \tilde{X}$

这里 $\tilde{X} \in \mathbb{R}^{n \times d}$ 是中心化数据矩阵（每行是一个 $\tilde{\mathbf{x}}_i^T$）。$C$ 是对称半正定矩阵。

于是 PCA 的优化问题是：

$\max_{\mathbf{w}} \; \mathbf{w}^T C \mathbf{w} \quad \text{s.t.} \quad \mathbf{w}^T \mathbf{w} = 1$

用 Lagrange 乘子法：构造 $\mathcal{L}(\mathbf{w}, \lambda) = \mathbf{w}^T C \mathbf{w} - \lambda(\mathbf{w}^T \mathbf{w} - 1)$，对 $\mathbf{w}$ 求导令其为零：

$\nabla_\mathbf{w} \mathcal{L} = 2C\mathbf{w} - 2\lambda\mathbf{w} = \mathbf{0}$

$\boxed{C\mathbf{w} = \lambda\mathbf{w}}$

这正是 $C$ 的特征值问题。最优的 $\mathbf{w}$ 是 $C$ 的特征向量，而投影方差 $\mathbf{w}^T C \mathbf{w} = \lambda \mathbf{w}^T \mathbf{w} = \lambda$。要最大化方差，取最大特征值对应的特征向量——这就是第一主成分（1st principal component）。

类推：第 $k$ 主成分是在与前 $k-1$ 个主成分正交的约束下，方差最大的方向。由于 $C$ 对称，它的特征向量自然正交（谱定理），所以 $d$ 个主成分就是 $C$ 的 $d$ 个特征向量，按特征值从大到小排列：

$\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_d \geq 0$

视角二：最小重建误差

换一个角度：我们想用 $k$ 维子空间近似原始 $d$ 维数据，使重建误差最小。

设 $W_k = [\mathbf{w}_1, \ldots, \mathbf{w}_k] \in \mathbb{R}^{d \times k}$ 是 $k$ 个正交基向量构成的矩阵。数据点 $\tilde{\mathbf{x}}_i$ 投影到这个子空间后的重建为 $\hat{\mathbf{x}}_i = W_k W_k^T \tilde{\mathbf{x}}_i$。重建误差是：

$\mathcal{E} = \sum_{i=1}^n \|\tilde{\mathbf{x}}_i - W_k W_k^T \tilde{\mathbf{x}}_i\|^2 = \|\tilde{X} - \tilde{X} W_k W_k^T\|_F^2$

由于 $W_k W_k^T$ 是投影矩阵，可以证明（利用 $\tilde{X}$ 的 SVD）：

$\mathcal{E} = \sum_{i=k+1}^d \lambda_i$

即重建误差等于被丢弃的特征值之和。最小化重建误差 $\iff$ 最大化保留的特征值之和 $\iff$ 选择前 $k$ 个最大的特征值对应的特征向量——与视角一的结论完全一致。

视角三：SVD

对中心化数据矩阵 $\tilde{X} \in \mathbb{R}^{n \times d}$ 做 SVD：

$\tilde{X} = U \Sigma V^T$

其中 $U \in \mathbb{R}^{n \times n}$, $\Sigma \in \mathbb{R}^{n \times d}$, $V \in \mathbb{R}^{d \times d}$。

协方差矩阵变为：

$C = \frac{1}{n-1}\tilde{X}^T\tilde{X} = \frac{1}{n-1}V\Sigma^T U^T U \Sigma V^T = \frac{1}{n-1}V\Sigma^T\Sigma V^T = V \underbrace{\frac{\Sigma^T\Sigma}{n-1}}_{\Lambda} V^T$

这正是 $C = V\Lambda V^T$ 的特征分解。因此：

• $V$ 的列（SVD 的右奇异向量）= $C$ 的特征向量 = 主成分方向
• 特征值 $\lambda_i = \sigma_i^2 / (n-1)$，其中 $\sigma_i$ 是 $\tilde{X}$ 的奇异值
• 数据在主成分上的坐标（即”得分”scores）为 $\tilde{X}V_k = U_k \Sigma_k$

SVD 视角的实践优势：直接对 $\tilde{X}$ 做 SVD，不需要显式计算 $d \times d$ 的协方差矩阵——当 $n \ll d$ 时（如人脸图像，像素维度远大于样本数）这尤为重要。

信息保留率

保留前 $k$ 个主成分的信息保留率（也叫方差解释率，explained variance ratio）：

$\rho_k = \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^d \lambda_i} = \frac{\sum_{i=1}^k \sigma_i^2}{\sum_{i=1}^d \sigma_i^2}$

实践中常选 $k$ 使 $\rho_k \geq 0.95$（保留 95% 以上的方差）。如果前几个特征值远大于其余，说明数据有很强的低维结构。

可视化：2D 数据的 PCA

下面的交互组件展示 PCA 在 2D 点云上的效果。数据呈椭圆分布——PC1（红色）沿方差最大的方向，PC2（紫色）沿正交方向。滑动 $k$ 从 2 到 1，观察投影到 PC1 后数据如何压缩到一维，以及方差解释率如何变化。

PCA 可视化：主成分与投影

数据点 → 主成分方向 → 投影降维

保留主成分数: k =2

方差解释率

90.0%

10.0%

PC1: λ₁ = 4.564PC2: λ₂ = 0.505

保留: 100.0%

当 $k = 1$ 时，绿色点是原始数据（浅蓝）投影到 PC1（红线）上的结果。虚线表示重建误差——每个点到其投影的距离。注意 PC1 方向的方差解释率已经很高，因为椭圆的长轴方向集中了大部分方差。

数值例子：手算 2D PCA

为了将上述理论落地，我们对一个具体的 2D 数据集手算完整的 PCA 流程。

数据（6 个二维点）：

$\mathbf{x}_1 = (2, 3),\; \mathbf{x}_2 = (3, 5),\; \mathbf{x}_3 = (5, 7),\; \mathbf{x}_4 = (4, 5),\; \mathbf{x}_5 = (6, 8),\; \mathbf{x}_6 = (4, 6)$

第一步：中心化

均值：

$\bar{\mathbf{x}} = \frac{1}{6}(2+3+5+4+6+4,\; 3+5+7+5+8+6) = (4, \; 5.667)$

中心化数据 $\tilde{\mathbf{x}}_i = \mathbf{x}_i - \bar{\mathbf{x}}$：

点	$\tilde{x}_1$	$\tilde{x}_2$
1	$-2$	$-2.667$
2	$-1$	$-0.667$
3	$1$	$1.333$
4	$0$	$-0.667$
5	$2$	$2.333$
6	$0$	$0.333$

第二步：协方差矩阵

$C = \frac{1}{n-1}\tilde{X}^T\tilde{X} = \frac{1}{5}\begin{bmatrix} \sum \tilde{x}_{i1}^2 & \sum \tilde{x}_{i1}\tilde{x}_{i2} \\ \sum \tilde{x}_{i1}\tilde{x}_{i2} & \sum \tilde{x}_{i2}^2 \end{bmatrix}$

计算各项：

• $\sum \tilde{x}_{i1}^2 = 4 + 1 + 1 + 0 + 4 + 0 = 10$
• $\sum \tilde{x}_{i1}\tilde{x}_{i2} = 5.333 + 0.667 + 1.333 + 0 + 4.667 + 0 = 12$
• $\sum \tilde{x}_{i2}^2 = 7.111 + 0.444 + 1.778 + 0.444 + 5.444 + 0.111 = 15.333$

$C = \frac{1}{5}\begin{bmatrix} 10 & 12 \\ 12 & 15.333 \end{bmatrix} = \begin{bmatrix} 2 & 2.4 \\ 2.4 & 3.067 \end{bmatrix}$

第三步：特征分解

特征方程 $\det(C - \lambda I) = 0$：

$(2 - \lambda)(3.067 - \lambda) - 2.4^2 = 0$

$\lambda^2 - 5.067\lambda + (6.133 - 5.76) = 0$

$\lambda^2 - 5.067\lambda + 0.373 = 0$

$\lambda = \frac{5.067 \pm \sqrt{25.67 - 1.493}}{2} = \frac{5.067 \pm \sqrt{24.18}}{2} = \frac{5.067 \pm 4.917}{2}$

$\lambda_1 = 4.992, \quad \lambda_2 = 0.075$

第四步：信息保留率

$\rho_1 = \frac{\lambda_1}{\lambda_1 + \lambda_2} = \frac{4.992}{5.067} = 98.5\%$

仅保留第一主成分就解释了 98.5% 的方差——这组数据几乎在一条直线上。

第五步：主成分方向

对 $\lambda_1 = 4.992$，解 $(C - \lambda_1 I)\mathbf{w} = \mathbf{0}$：

$\begin{bmatrix} 2 - 4.992 & 2.4 \\ 2.4 & 3.067 - 4.992 \end{bmatrix}\mathbf{w} = \begin{bmatrix} -2.992 & 2.4 \\ 2.4 & -1.925 \end{bmatrix}\mathbf{w} = \mathbf{0}$

从第一行：$-2.992 w_1 + 2.4 w_2 = 0$，得 $w_2 / w_1 = 2.992/2.4 = 1.247$。归一化后：

$\mathbf{w}_1 = \frac{1}{\sqrt{1 + 1.247^2}}(1, 1.247)^T \approx (0.625, 0.780)^T$

验证：$C\mathbf{w}_1 = (2 \times 0.625 + 2.4 \times 0.780,\; 2.4 \times 0.625 + 3.067 \times 0.780) = (3.122, 3.892)$，而 $\lambda_1 \mathbf{w}_1 = 4.992 \times (0.625, 0.780) = (3.120, 3.894)$。微小差异来自四舍五入。 ✓

第一主成分方向约为 $(0.625, 0.780)$——接近 $45°$ 但稍偏向 $x_2$ 轴，与数据的线性趋势一致。

Eigenfaces：PCA 的经典应用

1991 年，Turk 和 Pentland 在 Eigenfaces for Recognition 中将 PCA 应用于人脸识别，这是 PCA 作为降维工具最著名的成功案例。

核心思想

一张灰度人脸图像可以看作一个高维向量。假设图像大小为 $h \times w$ 像素，将每张图展平（flatten）为一个 $d = h \times w$ 维的向量 $\mathbf{x}_i \in \mathbb{R}^d$。例如，$112 \times 92$ 的图像对应 $d = 10304$ 维向量。

给定 $n$ 张训练图像 $\mathbf{x}_1, \ldots, \mathbf{x}_n$，Eigenfaces 方法的步骤是：

1. 计算均值脸

$\bar{\mathbf{x}} = \frac{1}{n}\sum_{i=1}^n \mathbf{x}_i$

均值脸是所有训练图像的”平均长相”。如果将 $\bar{\mathbf{x}}$ 重新 reshape 为 $h \times w$ 图像，它看起来像一张模糊的、特征平均的脸。

2. 中心化

$\tilde{\mathbf{x}}_i = \mathbf{x}_i - \bar{\mathbf{x}}$

每张脸减去均值脸——剩下的是”与平均脸的偏差”。

3. 对中心化数据做 PCA

理论上应该计算 $d \times d$ 协方差矩阵 $C = \frac{1}{n-1}\tilde{X}^T\tilde{X}$ 并求其特征向量。但 $d$ 通常远大于 $n$（例如 $d = 10304$, $n = 400$），计算 $10304 \times 10304$ 矩阵的特征分解代价巨大。

Turk-Pentland 技巧：改为计算 $n \times n$ 的小矩阵 $L = \tilde{X}\tilde{X}^T$。如果 $L\mathbf{v}_i = \mu_i \mathbf{v}_i$，那么 $\mathbf{u}_i = \tilde{X}^T \mathbf{v}_i$ 是 $C$ 的特征向量（可验证 $C\mathbf{u}_i = \frac{\mu_i}{n-1}\mathbf{u}_i$）。这将计算量从 $O(d^3)$ 降到 $O(n^3)$——当 $n \ll d$ 时大幅加速。

等价地，这就是对 $\tilde{X}$ 做 SVD 的 SVD 视角。

4. 选取前 $k$ 个主成分——“Eigenfaces”

前 $k$ 个特征向量 $\mathbf{w}_1, \ldots, \mathbf{w}_k \in \mathbb{R}^d$，每个都是与原始图像同维度的向量。将它们 reshape 回 $h \times w$ 图像，就得到了 eigenfaces——它们看起来像幽灵般的人脸，每个捕获了面部变化的一个”模式”（如光照方向、五官轮廓、面部对称性等）。

5. 人脸表示与识别

任何人脸图像（中心化后）可以近似表示为 eigenfaces 的线性组合：

$\tilde{\mathbf{x}} \approx \sum_{j=1}^k c_j \mathbf{w}_j, \quad c_j = \mathbf{w}_j^T \tilde{\mathbf{x}}$

系数向量 $\mathbf{c} = (c_1, \ldots, c_k)^T \in \mathbb{R}^k$ 是这张脸在”脸空间”（face space）中的 $k$ 维坐标。原始的 $d = 10304$ 维空间被压缩到了 $k$（通常 $k \approx 100 \text{-} 200$）维。

识别时，计算查询图像的系数 $\mathbf{c}_{\text{query}}$，与数据库中每张脸的系数比较欧氏距离，最近的就是识别结果。

为什么有效

Eigenfaces 之所以有效，本质上是因为人脸图像的协方差矩阵具有快速衰减的特征值谱。虽然图像名义上是万维向量，但人脸的结构约束（两只眼睛、一个鼻子、对称性等）使得有效维度远低于像素维度。前 100-200 个 eigenfaces 就能保留 95% 以上的方差——这与 Art. 5 矩阵微积分 中讨论的 intrinsic dimensionality 现象一脉相承。

PCA 的局限

PCA 强大但并非万能。了解其局限有助于在适当场景选择替代方法。

1. 只能捕获线性关系

PCA 找到的是方差最大的线性子空间。如果数据的内在结构是非线性的（如流形），PCA 可能表现不佳。

例如，一个弯曲的”瑞士卷”（Swiss Roll）数据集，PCA 投影会将本该分开的点混在一起，因为它只能找到全局的线性方向，无法跟踪局部的弯曲结构。

替代方法：Kernel PCA（将数据映射到高维特征空间再做 PCA）、t-SNE、UMAP 等非线性降维方法。Kernel 方法将在 Art. 20 Kernel 中深入讨论。

2. 对离群值敏感

协方差矩阵 $C$ 中每个元素都涉及平方项（$\tilde{x}_i \tilde{x}_j$），因此离群值对 $C$ 的影响被放大。少量离群点可能严重扭曲主成分方向。

替代方法：Robust PCA（Art. 12 Robust PCA）将数据矩阵分解为 $X = L + S$，其中 $L$ 是低秩的”正常”部分，$S$ 是稀疏的离群部分。

3. 正交约束导致缺乏可解释性

PCA 的主成分是正交的——这在数学上优雅，但在某些场景下导致主成分难以解释。例如，在文本挖掘中，主成分的元素有正有负，无法直接对应”主题”。

替代方法：NMF（非负矩阵分解）（Art. 9 NMF）要求因子矩阵的元素非负，产生的基向量更容易解释为”部分的叠加”（parts-based representation）。

4. 假设方差 = 信息

PCA 将”方差最大的方向”等同于”最重要的方向”。但在分类任务中，方差最大的方向可能是噪声方向，而区分类别的方向可能方差很小。

替代方法：LDA（Linear Discriminant Analysis）最大化类间方差与类内方差之比，而非总方差。

总结与展望

PCA 是分解工具的第一个经典应用——三件工具在这里完美交汇：

• 特征分解（Art. 2 特征分解）：$C\mathbf{w} = \lambda\mathbf{w}$ 是 PCA 的核心方程
• SVD（Art. 3 SVD）：$\tilde{X} = U\Sigma V^T$ 提供了数值稳定的计算路径和 $n \ll d$ 时的加速
• 矩阵微积分（Art. 5 矩阵微积分）：Lagrange 乘子法推导出 PCA = 协方差矩阵的特征值问题

关键公式回顾：

• 协方差矩阵：$C = \frac{1}{n-1}\tilde{X}^T\tilde{X}$
• PCA = 特征值问题：$C\mathbf{w} = \lambda\mathbf{w}$，最大方差方向 = 最大特征值的特征向量
• SVD 与 PCA 的关系：$\lambda_i = \sigma_i^2 / (n-1)$，主成分方向 = 右奇异向量 $V$，坐标 = $U_k\Sigma_k$
• 信息保留率：$\rho_k = \sum_{i=1}^k \lambda_i \big/ \sum_{i=1}^d \lambda_i$
• Eigenfaces：人脸 $\approx$ 均值脸 + eigenfaces 的线性组合

至此，我们完成了 Part 1 “拆”的核心弧线：从特征分解（结构）到 SVD（推广）到范数（度量）到微积分（优化）到 PCA（应用）。下一篇开始 Phase 3——随机化 SVD，当矩阵大到精确 SVD 不可行时，如何用随机采样以概率保证逼近前 $k$ 个奇异向量。