《机器学习的数学基础》第10章"主成分分析的降维"

第10章 主成分分析的降维

Dimensionality Reduction with Principal Component Analysis

直接处理高维数据（例如图像）会带来一些困难：它很难分析、难以解释、几乎无法可视化，并且（从实际角度来看）存储这些数据向量的代价可能很高。然而，高维数据往往具有一些可以利用的性质。例如，高维数据通常是 过完备的（overcomplete），即许多维度是冗余的，可以由其他维度的组合来解释。此外，高维数据中的各个维度往往是相关的，因此数据实际上存在一个 内在的低维结构。降维 就是利用这种结构和相关性，使我们能够以更紧凑的方式表示数据，理想情况下还能避免信息丢失。我们可以将降维看作是一种压缩技术，类似于 jpeg 或 mp3，它们分别是图像和音乐的压缩算法。

在本章中，我们将讨论 主成分分析（PCA），这是一种线性的降维算法。PCA 由 Pearson (1901) 和 Hotelling (1933) 提出，至今已有一百多年历史，但仍然是数据压缩和数据可视化中最常用的技术之一。它还被广泛用于识别高维数据中的简单模式、潜在因子以及数据结构。在信号处理领域，PCA 也被称为 Karhunen-Loève 变换。在本章中，我们将从最基本的原理推导 PCA，依赖于我们对 基和基变换（第 2.6.1 和 2.7.2 节）、投影（第 3.8 节）、特征值（第 4.2 节）、高斯分布（第 6.5 节）以及 约束优化（第 7.2 节）的理解。

降维通常利用高维数据（例如图像）的一个性质：它们往往位于低维子空间上。图 10.1 给出了一个二维的示例说明。虽然图 10.1(a) 中的数据并不完全落在一条直线上，但它在 \(x_2\)-方向上的变化很小，因此我们几乎可以把它看作是落在一条直线上——几乎没有信息损失；见图 10.1(b)。为了描述图 10.1(b) 中的数据，只需要 \(x_1\)-坐标即可，此时数据位于 \(\mathbb{R}^2\) 的一个一维子空间中。

10.1 问题设定

在主成分分析（PCA）中，我们希望找到数据点 \(x_n\) 的投影 \(\tilde{x}_n\)，使得这些投影尽可能地接近原始数据点，但其内在维度却显著降低。图 10.1 给出了这种情形的一个示意。更具体地说，我们考虑一个独立同分布的数据集

\[ X = \{x_1, \dots, x_N\}, \quad x_n \in \mathbb{R}^D, \]

其均值为 0，并具有数据协方差矩阵 (式 6.42)

\[ S = \frac{1}{N} \sum_{n=1}^N x_n x_n^\top. \tag{10.1} \]

此外，我们假设存在一个低维压缩表示（编码）

\[ z_n = B^\top x_n \in \mathbb{R}^M, \tag{10.2} \]

其中投影矩阵（projection matrix）定义为

\[ B := [b_1, \dots, b_M] \in \mathbb{R}^{D \times M}. \tag{10.3} \]

我们假设矩阵 \(B\) 的列向量是正交归一的（定义 3.7），即当且仅当 \(i \neq j\) 时，\(b_i^\top b_j = 0\)，并且 \(b_i^\top b_i = 1\)。我们希望找到一个 \(M\) 维子空间

\[ U \subseteq \mathbb{R}^D, \quad \dim(U) = M < D, \]

并将数据投影到该子空间上。我们记投影后的数据为 \(\tilde{x}_n \in U\)，而它们相对于子空间 \(U\) 的基向量 \(b_1, \dots, b_M\) 的坐标（即编码）为 \(z_n\)。我们的目标是找到投影\(\tilde{x}_n \in \mathbb{R}^D\)（或等价地，找到编码 \(z_n\) 以及基向量 \(b_1, \dots, b_M\)），使得这些投影尽可能接近原始数据 \(x_n\)，并使因压缩导致的损失最小化。

个人注：关于投影详见章节 “3.8 正交投影”

投影是一种特殊的线性变换。得出的最终投影的坐标还是原空间的表达！

例 10.1（坐标表示/编码）

考虑 \(\mathbb{R}^2\)，其标准基为

\[ e_1 = [1, 0]^\top, \quad e_2 = [0, 1]^\top. \]

在第 2 章中，我们知道任意 \(x \in \mathbb{R}^2\) 都可以表示为这些基向量的线性组合，例如：

\[ \begin{bmatrix} 5 \\ 3 \end{bmatrix} = 5 e_1 + 3 e_2. \tag{10.4} \]

然而，当我们只考虑如下形式的向量：

\[ \tilde{x} = \begin{bmatrix} 0 \\ z \end{bmatrix} \in \mathbb{R}^2, \quad z \in \mathbb{R}, \tag{10.5} \]

它们总可以写作 \(0 e_1 + z e_2\)。因此，要表示这些向量，仅需记住/存储 \(\tilde{x}\) 在基向量 \(e_2\) 上的坐标（即编码）\(z\)。更准确地说，这些 \(\tilde{x}\) 向量的集合（配合标准的向量加法和数乘运算）构成一个向量子空间 \(U\)（参见第 2.4 节），其维度为

\[ \dim(U) = 1, \quad 因为 \; U = \text{span}[e_2]. \]

在 第 10.2 节 中，我们将寻找能尽可能保留信息、并最小化压缩损失的低维表示。另一种推导 PCA 的方法将在 第 10.3 节 中给出，在那里我们将通过最小化原始数据 \(x_n\) 与其投影 \(\tilde{x}_n\) 之间的平方重建误差

\[ \|x_n - \tilde{x}_n\|^2 \]

来得到结果。

图 10.2 展示了 PCA 中所考虑的设定，其中 \(z\) 表示压缩数据 \(\tilde{x}\) 的低维表示，并扮演“瓶颈”的角色，控制着信息在 \(x\) 与 \(\tilde{x}\) 之间的流动量。在 PCA 中，我们考虑原始数据 \(x\) 与其低维编码 \(z\) 之间的线性关系，即 \[ z = B^\top x, \quad \tilde{x} = Bz, \]

其中 \(B\) 是一个合适的矩阵。基于将 PCA 理解为数据压缩技术的动机，我们可以将图 10.2 中的箭头解释为表示“编码器”和“解码器”的一对操作。矩阵 \(B\) 表示的线性映射可以看作 解码器，它将低维编码 \(z \in \mathbb{R}^M\) 映射回原始数据空间 \(\mathbb{R}^D\)。类似地，\(B^\top\) 可以看作 编码器，它将原始数据 \(x\) 编码为低维（压缩）表示 \(z\)。

在本章中，我们将多次使用 MNIST 手写数字数据集 作为示例。该数据集包含 60,000 个手写数字（0 到 9）的样本。每个数字是一个大小为 \(28 \times 28\) 的灰度图像，即包含 784 个像素，因此我们可以将数据集中的每张图片看作一个向量

\[ x \in \mathbb{R}^{784}. \]

这些数字的示例如图 10.3 所示。

10.2 最大方差视角

图 10.1 给出了一个示例，展示了如何用单一坐标表示一个二维数据集。在图 10.1(b) 中，我们选择忽略了数据的 \(x_2\) 坐标，因为它并没有带来太多信息，从而使压缩后的数据仍与原始数据（图 10.1(a)）非常接近。我们本可以选择忽略 \(x_1\) 坐标，但那样得到的压缩数据将与原始数据差别很大，许多信息都会丢失。

如果我们将数据中的信息含量理解为数据“填充空间”的程度，那么可以通过观察数据的分布范围来描述其中包含的信息量。根据 第 6.4.1 节，我们知道方差是数据分布范围的一个指标。由此可以将 PCA 推导为一种维度约简算法，它通过在低维表示中最大化方差来保留尽可能多的信息。图 10.4 说明了这一点。

个人注：所以为了获得最大的信息量，从而保留方差最大的方向。

结合 第 10.1 节 的设定，我们的目标是找到一个矩阵 \(B\)（见 (10.3)），当数据被投影到由 \(B\) 的列向量 \(b_1, \dots, b_M\) 所张成的子空间时，能够尽可能保留信息。在数据压缩后保留最多信息等价于在低维编码中捕获最大的方差（Hotelling, 1933）。

备注（居中数据）： 在 (10.1) 中的数据协方差矩阵里，我们假设了数据是中心化的。这个假设不会带来普遍性损失。假设 \(\mu\) 是数据的均值。利用 第 6.4.4 节 中讨论过的方差性质，我们有： \[ V_z[z] = V_x[B^\top (x - \mu)] = V_x[B^\top x - B^\top \mu] = V_x[B^\top x], \tag{10.6} \]

即低维编码的方差与数据的均值无关。因此，我们可以不失一般性地假设数据均值为 0。在此假设下，低维编码的均值也为 0，因为

\[ E_z[z] = E_x[B^\top x] = B^\top E_x[x] = 0. \]

10.2.1 最大方差方向

我们采用一种 逐步方法（a sequential approach） 来最大化低维编码的方差。首先，我们寻找一个向量

\[ b_1 \in \mathbb{R}^D \]

使得投影数据的方差最大化。换句话说，我们希望最大化低维表示 \(z \in \mathbb{R}^M\) 的第一个坐标 \(z_1\) 的方差：

\[ V_1 := V[z_1] = \frac{1}{N} \sum_{n=1}^N z_{1n}^2, \tag{10.7} \]

其中我们利用了数据的 i.i.d. 假设，并定义 \(z_{1n}\) 为样本 \(x_n \in \mathbb{R}^D\) 的低维表示 \(z_n \in \mathbb{R}^M\) 的第一个坐标。注意到，\(z_n\) 的第一个分量由下式给出：

\[ z_{1n} = b_1^\top x_n, \tag{10.8} \]

即，它是向量 \(x_n\) 在由 \(b_1\) 张成的一维子空间上的正交投影的坐标（见第 3.8 节）。将 (10.8) 代入 (10.7)，得到：

\[ \begin{align} V_1 &= \frac{1}{N} \sum_{n=1}^N (b_1^\top x_n)^2 \\ &= b_1^\top \left(\frac{1}{N} \sum_{n=1}^N x_n x_n^\top \right) b_1 \tag{10.9a} \\ &= b_1^\top S b_1, \tag{10.9b} \end{align} \]

其中 \(S\) 是在 (10.1) 中定义的数据协方差矩阵。在 (10.9a) 中，我们利用了向量内积关于参数对称的性质，即

\[ b_1^\top x_n = x_n^\top b_1. \]

注意，如果任意放大向量 \(b_1\) 的长度，就会使 \(V_1\) 增大。例如，将 \(b_1\) 的长度变为原来的两倍，可能使 \(V_1\) 增大到原来的四倍。因此，我们必须施加约束条件 \(\|b_1\|^2 = 1\)，从而将问题转化为一个约束优化问题： \[ \max_{b_1} \; b_1^\top S b_1 \\ \quad \text{subject to } \|b_1\|^2 = 1. \tag{10.10} \]

根据 第 7.2 节，该约束优化问题的拉格朗日函数为：

\[ L(b_1, \lambda_1) = b_1^\top S b_1 + \lambda_1 (1 - b_1^\top b_1). \tag{10.11} \]

分别对 \(b_1\) 和 \(\lambda_1\) 求偏导数：

\[ \frac{\partial L}{\partial b_1} = 2 b_1^\top S - 2 \lambda_1 b_1^\top, \quad \frac{\partial L}{\partial \lambda_1} = 1 - b_1^\top b_1. \tag{10.12} \]

令偏导为 0，可得：

\[ S b_1 = \lambda_1 b_1, \tag{10.13} \]

\[ b_1^\top b_1 = 1. \tag{10.14} \]

对比 第 4.4 节 的特征值分解定义，可以发现：

• \(b_1\) 是数据协方差矩阵 \(S\) 的特征向量；
• 拉格朗日乘子 \(\lambda_1\) 即对应的特征值。

因此，方差目标 (10.10) 可以改写为：

\[ V_1 = b_1^\top S b_1 = \lambda_1 b_1^\top b_1 = \lambda_1. \tag{10.15} \]

也就是说，将数据投影到某一维子空间后的方差，等于该子空间所对应的基向量 \(b_1\) 的特征值。为了最大化低维编码的方差，我们应选择数据协方差矩阵最大特征值所对应的特征向量作为基向量。这个特征向量被称为 第一主成分（first principal component）。

个人注：在标准 PCA 中，主成分方向就是协方差矩阵的特征向量（必须是单位长度），\(\lambda_1\) 是沿该方向的方差。

主成分本身是“投影出来的新坐标”

由于 \(S\) 是实对称矩阵，不同特征值对应的特征向量必定正交。再加上我们通常会把每个特征向量单位化，得到一组正交单位向量： \[ b_i^\top b_j = \delta_{ij} \]

我们可以通过将坐标 \(z_{1n}\) 映射回数据空间，来观察主成分 \(b_1\) 在原始数据空间中的作用/贡献，从而得到投影后的数据点： \[ \tilde{x}_n = b_1 z_{1n} = b_1 b_1^\top x_n \in \mathbb{R}^D. \tag{10.16} \]

备注： 尽管 \(\tilde{x}_n\) 是一个 \(D\)-维向量，但它只需要一个坐标 \(z_{1n}\) 就能在基向量 \(b_1 \in \mathbb{R}^D\) 下表示出来。

10.2.2 最大方差的 M 维子空间

假设我们已经找到了前 \(m-1\) 个主成分，即协方差矩阵 \(S\) 的前 \(m-1\) 个最大特征值所对应的特征向量。由于 \(S\) 是对称矩阵，根据 谱定理（定理 4.15），我们可以用这些特征向量构造出一个 \(\mathbb{R}^D\) 中维度为 \(m-1\) 的正交规范特征子空间。

一般而言，第 \(m\) 个主成分可以通过从数据中去除前 \(m-1\) 个主成分 \(b_1, \ldots, b_{m-1}\) 的影响而得到，从而在剩余信息中继续寻找主成分。这样，我们得到新的数据矩阵：

\[ \hat{X} := X - \sum_{i=1}^{m-1} b_i b_i^\top X \;=\; X - B_{m-1} X, \tag{10.17} \]

其中 \(X = [x_1, \ldots, x_N] \in \mathbb{R}^{D \times N}\) 的列向量为数据点， 而

\[ B_{m-1} := \sum_{i=1}^{m-1} b_i b_i^\top \]

是一个投影矩阵，将数据投影到由 \(b_1, \ldots, b_{m-1}\) 张成的子空间上。

个人注：注意定义，\(b_i b_i^\top X\)是\(X\)投影到\(b_i\)后的坐标（原空间的）。

备注（符号约定）：在本章中，我们不采用将数据点 \(x_1, \ldots, x_N\) 作为矩阵行向量的惯例，而是将它们作为矩阵 \(X\) 的列向量。这意味着我们的数据矩阵 \(X\) 是一个 \(D \times N\) 矩阵，而非常见的 \(N \times D\) 矩阵。我们这样选择的原因是，在后续代数运算中可以避免频繁转置或将列向量重定义为行向量来左乘矩阵。

个人注：在线性回归分析中，拟合的模型是\(X \theta\)，\(X\)作为数据的时候，放在左边，这时候可以把数据看作行。

为了找到第 \(m\) 个主成分，我们希望最大化：

\[ V_m = V[z_m] = \frac{1}{N} \sum_{n=1}^N z_{mn}^2 = \frac{1}{N} \sum_{n=1}^N (b_m^\top \hat{x}_n)^2 = b_m^\top \hat{S} b_m, \tag{10.18} \]

其中约束为 \(\|b_m\|^2 = 1\)。这里我们遵循 (10.9b) 的相同步骤，并定义 \(\hat{S}\) 为变换后数据集\(\hat{X} := \{\hat{x}_1, \ldots, \hat{x}_N\}\) 的协方差矩阵。与第一主成分类似，我们得到一个约束优化问题，其最优解 \(b_m\) 是 \(\hat{S}\) 最大特征值所对应的特征向量。

结果表明，\(b_m\) 也是 \(S\) 的特征向量。更一般地说，\(S\) 和 \(\hat{S}\) 的特征向量集合是相同的。由于 \(S\) 和 \(\hat{S}\) 都是对称矩阵，根据谱定理 (4.15)，它们都有 \(D\) 个两两正交的特征向量。下面我们来证明：每一个 \(S\) 的特征向量同时也是 \(\hat{S}\) 的特征向量。

假设我们已经得到了 \(\hat{S}\) 的前 \(m-1\) 个特征向量 \(b_1, \ldots, b_{m-1}\)。考虑 \(S\) 的某个特征向量 \(b_i\)，即：

\[ S b_i = \lambda_i b_i. \]

一般情况下：

\[ \hat{S} b_i = \frac{1}{N} \hat{X} \hat{X}^\top b_i = \frac{1}{N} (X - B_{m-1} X)(X - B_{m-1} X)^\top b_i \tag{10.19a} \]

\[ = (S - S B_{m-1} - B_{m-1} S + B_{m-1} S B_{m-1}) b_i. \tag{10.19b} \]

我们区分两种情况：

• 情况 1： \(i \geq m\)。 此时，\(b_i\) 不是前 \(m-1\) 个主成分中的一个，因此与它们正交，即

\[ B_{m-1} b_i = 0. \]

• 情况 2： \(i < m\)。 此时，\(b_i\) 是前 \(m-1\) 个主成分之一，因此属于由 \(b_1, \ldots, b_{m-1}\) 张成的子空间。由于它们构成该子空间的 ONB，有：

\[ B_{m-1} b_i = b_i. \]

总结如下：

\[ B_{m-1} b_i = \begin{cases} b_i, & i < m \\ 0, & i \geq m \end{cases} \tag{10.20} \]

在 情况 1 (\(i \geq m\)) 下，将 (10.20) 代入 (10.19b)，得：

\[ \hat{S} b_i = (S - B_{m-1} S) b_i = S b_i = \lambda_i b_i, \]

即 \(b_i\) 也是 \(\hat{S}\) 的特征向量，特征值为 \(\lambda_i\)。特别地：

\[ \hat{S} b_m = S b_m = \lambda_m b_m. \tag{10.21} \]

这表明，\(b_m\) 既是 \(S\) 的特征向量，也是 \(\hat{S}\) 的特征向量。更具体地说，\(\lambda_m\) 是 \(\hat{S}\) 的最大特征值，同时也是 \(S\) 的第 \(m\) 大特征值，它们的特征向量相同。

在 情况 2 (\(i < m\)) 下，将 (10.20) 代入 (10.19b)，得：

\[ \hat{S} b_i = (S - S B_{m-1} - B_{m-1} S + B_{m-1} S B_{m-1}) b_i = 0 = 0 b_i. \tag{10.22} \]

这意味着 \(b_1, \ldots, b_{m-1}\) 也是 \(\hat{S}\) 的特征向量，但它们对应的特征值为 0，因此张成了 \(\hat{S}\) 的零空间。

综上所述，\(S\) 的每个特征向量都是 \(\hat{S}\) 的特征向量。但如果这些特征向量属于前 \(m-1\) 个主成分子空间，那么它们在 \(\hat{S}\) 中对应的特征值为 0。由 (10.21) 以及约束 \(b_m^\top b_m = 1\)，数据投影到第 \(m\) 个主成分上的方差为：

\[ V_m = b_m^\top S b_m = \lambda_m b_m^\top b_m = \lambda_m. \tag{10.23} \]

这意味着，当数据投影到一个 \(M\) 维子空间时，其方差等于数据协方差矩阵中与该子空间对应的特征向量的特征值之和

例 10.2（MNIST 中数字 “8” 的特征值） 取出 MNIST 训练集中所有的数字 “8”，我们计算这些数据的协方差矩阵的特征值。图 10.5(a) 展示了数据协方差矩阵最大的 200 个特征值。我们看到，只有很少的一部分特征值与 0 有显著差异。因此，当把数据投影到由相应特征向量所张成的子空间时，大部分的方差仅由少数几个主成分捕捉，如图 10.5(b) 所示。

总体而言，为了在 \(\mathbb{R}^D\) 中找到一个尽可能保留信息的 \(M\) 维子空间，PCA 告诉我们，应选择矩阵 \(B_{in}\)（公式 (10.3) 中）对应的 \(M\) 个特征向量——这些特征向量对应于数据协方差矩阵 \(S\) 的最大的 \(M\) 个特征值。PCA 在前 \(M\) 个主成分中能捕捉到的最大方差为

\[ V_M = \sum_{m=1}^M \lambda_m, \tag{10.24} \]

其中 \(\lambda_m\) 是数据协方差矩阵 \(S\) 的最大的 \(M\) 个特征值。相应地，通过 PCA 压缩数据所丢失的方差为

\[ J_M := \sum_{j=M+1}^D \lambda_j = V_D - V_M. \tag{10.25} \]

与其使用这些绝对量，我们也可以定义相对方差：

• 被捕捉的相对方差：\(\frac{V_M}{V_D}\)，
• 压缩导致的相对方差损失：\(1 - \frac{V_M}{V_D}\)。

10.3 投影视角

在本节中，我们将把 PCA 推导为一种直接最小化平均重构误差的算法。这种视角使我们能够把 PCA 理解为一种最优的线性自编码器。这里我们将大量借鉴第 2 章和第 3 章的内容。

we will derive PCA as an algorithm that directly mini-mizes the average reconstruction error. This perspective allows us to interpret PCA as implementing an optimal linear auto-encoder.

在上一节中，我们是通过最大化投影空间中的方差来推导 PCA 的，以尽可能保留信息。而在这里，我们将考察原始数据点 \(x_n\) 与其重构 \(\tilde{x}_n\) 之间的差异向量，并最小化这一距离，使得 \(x_n\) 与 \(\tilde{x}_n\) 尽可能接近。图 10.6 展示了这一设置。

10.3.1 设置与目标

假设在 \(\mathbb{R}^D\) 中有一个（有序的）正交归一基（ONB）\(B = (b_1, \dots, b_D),\)即当且仅当 \(i=j\) 时 \(b_i^\top b_j = 1\)，否则为 0。根据第 2.5 节，我们知道，对于 \(\mathbb{R}^D\) 的一组基 \((b_1,\dots,b_D)\)，任意 \(x \in \mathbb{R}^D\) 都可以表示为基向量的线性组合，即 \[ x = \sum_{d=1}^D \zeta_d b_d = \sum_{m=1}^M \zeta_m b_m + \sum_{j=M+1}^D \zeta_j b_j, \tag{10.26} \] 其中 \(\zeta_d \in \mathbb{R}\) 为相应的坐标。我们感兴趣的是寻找向量 \(\tilde{x} \in \mathbb{R}^D\)，它位于某个低维子空间\(U \subseteq \mathbb{R}^D, \ \dim(U)=M\)，使得 \[ \tilde{x} = \sum_{m=1}^M z_m b_m \in U \subseteq \mathbb{R}^D \tag{10.27} \] 与原始数据 \(x\) 尽可能相似。注意，此时我们必须假设 \(\tilde{x}\) 的坐标 \(z_m\) 与 \(x\) 的坐标 \(\zeta_m\) 不完全相同。

接下来，我们正是利用这种形式的 \(\tilde{x}\)，来寻找最优的坐标 \(z\) 和基向量 \(b_1,\dots,b_M\)，使得 \(\tilde{x}\) 与原始数据点 \(x\) 尽可能接近，即最小化两者之间的欧几里得距离 \(\|x - \tilde{x}\|\)。图 10.7 展示了这一情况。

为了简化推导，不失一般性地，我们假设数据集\(X = \{x_1,\dots,x_N\}, \quad x_n \in \mathbb{R}^D,\)是零均值的，即 \(E[X]=0\)。如果不做零均值的假设，最终我们也会得到完全相同的解，但符号表达会变得非常繁琐。

我们的目标是：找到把 \(X\) 投影到低维子空间 \(U \subseteq \mathbb{R}^D\)（\(\dim(U)=M\)）的最佳线性投影，其中子空间由正交归一基向量 \(b_1,\dots,b_M\) 张成。我们称该子空间 \(U\) 为主子空间(principal subspace.)。数据点的投影表示为

\[ \tilde{x}_n := \sum_{m=1}^M z_{mn} b_m = B z_n \in \mathbb{R}^D, \tag{10.28} \] 其中\(z_n := [z_{1n}, \dots, z_{Mn}]^\top \in \mathbb{R}^M\)是 \(\tilde{x}_n\) 在基 \((b_1,\dots,b_M)\) 下的坐标向量。换句话说，我们希望 \(\tilde{x}_n\) 尽可能接近原始点 \(x_n\)。

在这里，我们采用的相似度度量是 \(x\) 与 \(\tilde{x}\) 之间的平方欧几里得距离 \(\|x - \tilde{x}\|^2\)。因此，我们定义目标函数为最小化平均平方欧几里得距离（重构误差）（Pearson, 1901）： \[ J_M := \frac{1}{N} \sum_{n=1}^N \|x_n - \tilde{x}_n\|^2, \tag{10.29} \] 其中我们明确指出，数据投影的子空间维度是 \(M\)。为了找到最优的线性投影，我们需要确定主子空间的正交归一基，以及投影点在该基下的坐标向量 \(z_n \in \mathbb{R}^M\)。

为此，我们采用两步策略：

1. 在给定 ONB \((b_1,\dots,b_M)\) 的前提下，先优化坐标 \(z_n\)；
2. 再寻找最优的 ONB。

10.3.2 寻找最优坐标

我们先从寻找投影 \(\tilde{x}_n\) 的最优坐标 \(z_{1n}, \dots, z_{Mn}\) 开始（\(n=1,\dots,N\)）。考虑图 10.7(b)，其中主子空间由单个向量 \(b\) 张成。几何上，寻找最优坐标 \(z\) 对应于：在基向量 \(b\) 上找到 \(\tilde{x}\) 的表示，使得 \(\tilde{x}-x\) 的距离最小。从图 10.7(b) 可以直观看出，这将是正交投影，接下来我们将严格证明这一点。

假设 \(U \subseteq \mathbb{R}^D\) 的 ONB（正交归一基）为 \((b_1,\dots,b_M)\)。要在该基下找到最优坐标 \(z_m\)，我们需要计算偏导数： \[ \begin{align} \frac{\partial J_{M}}{\partial z_{in}} &= \frac{\partial J_{M}}{\partial \tilde{x}_{n}} \frac{\partial \tilde{x}_{n}}{\partial z_{in}}, \tag{10.30a}\\[0.5em] \frac{\partial J_{M}}{\partial \tilde{x}_{n}} &= -\frac{2}{N}\,(x_{n}-\tilde{x}_{n})^{\top} \in \mathbb{R}^{1 \times D}, \tag{10.30b} \end{align} \] 另一方面， \[ \frac{\partial \tilde{x}_n}{\partial z_{in}} = b_i, \quad i=1,\dots,M. \tag{10.30c} \] 结合 (10.28)，可得： \[ \frac{\partial J_M}{\partial z_{in}} = -\frac{2}{N} (x_n - \sum_{m=1}^M z_{mn} b_m)^\top b_i. \tag{10.31a} \] 由于 \(b_i^\top b_i = 1\)，化简得： \[ \frac{\partial J_M}{\partial z_{in}} = -\frac{2}{N}(x_n^\top b_i - z_{in}). \tag{10.31b} \] 令该偏导为 0，立即得到最优坐标：\(U^\perp = \text{span}[b_{M+1}, \dots, b_D]\)是\(U = \text{span}[b_1,\dots,b_M]\)的正交补（参见第 3.6 节）。

备注（带正交归一基向量的正交投影）

简要回顾第 3.8 节的正交投影。如果 \((b_1,\dots,b_D)\) 是 \(\mathbb{R}^D\) 的正交归一基，那么 \[ \tilde{x} = b_j(b_j^\top b_j)^{-1} b_j^\top x = b_j b_j^\top x \in \mathbb{R}^D \tag{10.33} \] 就是 \(x\) 在第 \(j\) 个基向量张成的一维子空间上的正交投影。而此时坐标 \[ z_j = b_j^\top x \] 正是该投影在基向量 \(b_j\) 下的坐标，因为 \(z_j b_j = \tilde{x}\)。图 10.8(b) 展示了这种情况。

更一般地，若我们希望把 \(x\) 投影到 \(\mathbb{R}^D\) 的一个 M 维子空间上，且该子空间的 ONB 为 \(b_1,\dots,b_M\)，则正交投影为： \[ \tilde{x} = B(B^\top B)^{-1}B^\top x = BB^\top x, \tag{10.34} \] 其中 \[ B = [b_1,\dots,b_M] \in \mathbb{R}^{D\times M}. \] 在有序基 \((b_1,\dots,bM)\) 下，该投影的坐标为 \[ z = B^\top x, \] 这与第 3.8 节的讨论一致。

我们可以将这些坐标理解为在新的坐标系 \((b_1,\dots,b_M)\) 下对投影向量的表示。需要注意的是，虽然 \(\tilde{x} \in \mathbb{R}^D\)，但仅需用 \(M\) 个坐标 \(z_1,\dots,z_M\) 来表示；其余相对于基 \((b_{M+1},\dots,b_D)\) 的 \(D-M\) 个坐标恒为 0。

到目前为止，我们已经证明：在给定 ONB 的情况下，可以通过正交投影找到 \(\tilde{x}\) 的最优坐标。接下来，我们将确定最优基是什么。

10.3.3 寻找主子空间的基

为了确定主子空间的基向量 \(b_1, \ldots, b_M\)，我们利用目前已有的结果重新表述损失函数 (10.29)。这样会使得寻找基向量更为容易。为了重新写损失函数，我们使用之前的推导，得到 \[ \tilde{x}_{n} = \sum_{m=1}^{M} z_{mn} \mathbf{b}_{m} \overset{(10.32)}{=} \sum_{m=1}^{M} \left(x_{n}^{\top} \mathbf{b}_{m}\right) \mathbf{b}_{m}. \tag{10.35} \] 接下来，我们利用点积的对称性，得到 \[ \tilde{x}_n = \left(\sum_{m=1}^M b_m b_m^\top \right) x_n. \tag{10.36} \] 由于我们通常可以把原始数据点 \(x_n\) 表示为所有基向量的线性组合，因此有 \[ \begin{align} x_{n} &= \sum_{d=1}^{D} z_{dn} \mathbf{b}_{d} \overset{(10.32)}{=} \sum_{d=1}^{D} \left(x_{n}^{\top} \mathbf{b}_{d}\right) \mathbf{b}_{d} = \left(\sum_{d=1}^{D} \mathbf{b}_{d} \mathbf{b}_{d}^{\top}\right) x_{n}, \tag{10.37a} \\[0.75em] &= \left(\sum_{m=1}^{M} \mathbf{b}_{m} \mathbf{b}_{m}^{\top}\right) x_{n} + \left(\sum_{j=M+1}^{D} \mathbf{b}_{j} \mathbf{b}_{j}^{\top}\right) x_{n}, \tag{10.37b} \end{align} \] 其中，我们将含 \(D\) 项的和拆分为含 \(M\) 项的部分与含 \(D-M\) 项的部分。由此可得，位移向量 \(\tilde{x}_n - x_n\)，即原始数据点与其投影的差向量为 \[ \tilde{x}_n - x_n = \sum_{j=M+1}^D b_j b_j^\top x_n \tag{10.38a} \]

\[ = \sum_{j=M+1}^D (x_n^\top b_j)b_j. \tag{10.38b} \] 这意味着，该差向量正好是数据点在主子空间的正交补上的投影：我们将 (10.38a) 中的矩阵 \(\sum_{j=M+1}^D b_j b_j^\top\) 识别为执行该投影的投影矩阵。因此，位移向量 \(x_n - \tilde{x}_n\) 位于与主子空间正交的子空间中，如图 10.9 所示。

备注（低秩近似）在 (10.38a) 中，我们看到将 \(x\) 投影到 \(\tilde{x}\) 的投影矩阵为 \[ \sum_{m=1}^M b_m b_m^\top = BB^\top. \tag{10.39} \] 由构造可知，作为秩一矩阵 \(b_m b_m^\top\) 的和，\(BB^\top\) 是对称的，并且秩为 \(M\)。因此，平均平方重构误差可以写为 \[ \begin{align} \frac{1}{N} \sum_{n=1}^N \|x_n - \tilde{x}_n\|^2 = \frac{1}{N} \sum_{n=1}^N \|x_n - BB^\top x_n\|^2 \tag{10.40a} \\ = \frac{1}{N} \sum_{n=1}^N \|(I - BB^\top)x_n\|^2. \tag{10.40b} \end{align} \] 寻找使原始数据 \(x_n\) 与其投影 \(\tilde{x}_n\) 之间差异最小的正交归一基向量 \(b_1, \ldots, b_M\)，等价于寻找单位矩阵 \(I\) 的最佳秩-\(M\) 近似 \(BB^\top\)（参见第 4.6 节。现在我们已经具备了重新表述损失函数 (10.29) 的所有工具：

\[ J_M = \frac{1}{N} \sum_{n=1}^{N} \|\boldsymbol{x}_n - \tilde{\boldsymbol{x}}_n\|^2 \overset{\text{(10.38b)}}{=} \frac{1}{N} \sum_{n=1}^{N} \left\| \sum_{j=M+1}^{D} (\boldsymbol{b}_j^\top \boldsymbol{x}_n) \boldsymbol{b}_j \right\|^2. \tag{10.41} \] 我们现在显式计算平方范数，并利用基向量 \(b_j\) 构成正交归一基 (ONB) 的事实，得到： \[ \begin{align} J_M & = \frac{1}{N} \sum_{n=1}^{N} \sum_{j=M+1}^{D} (\boldsymbol{b}_j^\top \boldsymbol{x}_n)^2 = \frac{1}{N} \sum_{n=1}^{N} \sum_{j=M+1}^{D} \boldsymbol{b}_j^\top \boldsymbol{x}_n \boldsymbol{b}_j^\top \boldsymbol{x}_n \tag{10.42a} \\ & = \frac{1}{N} \sum_{n=1}^{N} \sum_{j=M+1}^{D} \boldsymbol{b}_j^\top \boldsymbol{x}_n \boldsymbol{x}_n^\top \boldsymbol{b}_j, \tag{10.42b} \end{align} \] 其中，在最后一步中我们利用了点积的对称性，将 \(b_j^\top x_n = x_n^\top b_j\) 写出。 接下来交换求和次序： \[ J_M = \sum_{j=M+1}^D b_j^\top \Big( \frac{1}{N}\sum_{n=1}^N x_n x_n^\top \Big) b_j =: \sum_{j=M+1}^D b_j^\top S b_j, \tag{10.43a} \]

\[ = \sum_{j=M+1}^D \operatorname{tr}(b_j^\top S b_j) = \sum_{j=M+1}^D \operatorname{tr}(S b_j b_j^\top) = \operatorname{tr}\Big(S \sum_{j=M+1}^D b_j b_j^\top \Big), \tag{10.43b} \] 这里我们利用了迹算子 \(\operatorname{tr}(\cdot)\) 的性质（见公式 (4.18)）：其是线性的，并且对循环置换不变。由于我们假设数据集已经中心化，即 \(E[X] = 0\)，因此矩阵 \(S\) 就是数据的协方差矩阵。而 (10.43b) 中的投影矩阵由秩一矩阵 \(b_j b_j^\top\) 的和构成，因此其秩为 \(D-M\)。等式 (10.43a) 表明：我们可以将平均平方重构误差等价地表述为数据协方差矩阵在主子空间正交补上的投影。因此，最小化平均平方重构误差就等价于最小化数据在被舍弃子空间（即主子空间正交补）上的投影方差。换句话说，这也等价于最大化我们保留的主子空间上的投影方差。这就把“投影损失”的视角和第 10.2 节讨论的 PCA 最大方差公式直接联系起来。因此，我们会得到和最大方差视角相同的解。于是，这里不再重复推导，而是结合投影视角总结之前的结果。当投影到 \(M\) 维主子空间时，平均平方重构误差为： \[ J_M = \sum_{j=M+1}^D \lambda_j, \tag{10.44} \] 其中 \(\lambda_j\) 是数据协方差矩阵的特征值。因此，为了最小化 (10.44)，我们需要舍弃最小的 \(D-M\) 个特征值，其对应的特征向量就构成主子空间的正交补的基。于是，主子空间的基就是协方差矩阵对应于最大 \(M\) 个特征值的特征向量 \(b_1, \ldots, b_M\)。

例 10.3（MNIST 数字嵌入） 图 10.10 将 MNIST 数字 “0” 和 “1” 的训练数据可视化，这些数据被嵌入到由前两个主成分所张成的向量子空间中。我们可以观察到“0”（蓝点）和“1”（橙点）之间有相对清晰的分离，同时也能看到每个簇内部的变化情况。图中用红色标出了四个数字 “0” 和 “1” 的嵌入点，并展示了它们对应的原始数字。图像揭示了“0”的数据集内部变化明显大于“1”的数据集内部变化。

10.4 特征向量计算与低秩近似

在前面的章节中，我们得到了主子空间的基，它们是数据协方差矩阵中与最大特征值相关联的特征向量： \[ S = \frac{1}{N} \sum_{n=1}^{N} x_n x_n^\top \;=\; \frac{1}{N} X X^\top , \quad X = [x_1,\ldots, x_N] \in \mathbb{R}^{D \times N} \tag{10.45-10.46} \] 需要注意的是，\(X\) 是一个 \(D \times N\) 的矩阵，即它是“常见”数据矩阵的转置（Bishop, 2006; Murphy, 2012）。为了得到 \(S\) 的特征值（及对应的特征向量），我们有两种方法：

1. 对 \(S\) 做特征分解（见第 4.2 节），直接计算其特征值与特征向量。
2. 使用奇异值分解（见第 4.5 节）。由于 \(S\) 是对称矩阵，并且可以分解为 \(X X^\top\)（忽略因子 \(\tfrac{1}{N}\)），因此 \(S\) 的特征值就是 \(X\) 的奇异值的平方。

更具体地，\(X\) 的奇异值分解为： \[ X_{D \times N} = U_{D \times D} \; \Sigma_{D \times N} \; V^\top_{N \times N} \tag{10.47} \] 其中，\(U \in \mathbb{R}^{D \times D}\) 和 \(V^\top \in \mathbb{R}^{N \times N}\) 是正交矩阵，\(\Sigma \in \mathbb{R}^{D \times N}\) 是一个只有奇异值 \(\sigma_{ii} \geq 0\) 为非零的矩阵。于是我们有： \[ S = \frac{1}{N} X X^\top = \frac{1}{N} U \Sigma V^\top V \Sigma^\top U^\top = \frac{1}{N} U \Sigma \Sigma^\top U^\top \tag{10.48} \] 根据第 4.5 节的结果，矩阵 \(U\) 的列向量就是 \(X X^\top\)（因此也是 \(S\)）的特征向量。进一步地，\(S\) 的特征值 \(\lambda_d\) 与 \(X\) 的奇异值的关系为： \[ \lambda_d = \frac{\sigma_d^2}{N} \tag{10.49} \] 这种 \(S\) 的特征值与 \(X\) 的奇异值之间的关系，建立了最大方差视角（第 10.2 节）与奇异值分解之间的联系。

10.4.1 使用低秩矩阵近似的 PCA

为了最大化投影数据的方差（或最小化平均平方重构误差），PCA 选择式 (10.48) 中矩阵 \(U\) 的列向量作为数据协方差矩阵 \(S\) 的前 \(M\) 大特征值对应的特征向量。这样，我们可以将 \(U\) 识别为公式 (10.3) 中的投影矩阵 \(B\)，它把原始数据投影到一个维度为 \(M\) 的低维子空间上。

Eckart–Young 定理（见第 4.6 节的定理 4.25）提供了一种直接的方法来估计低维表示。考虑 \(X\) 的最佳秩-\(M\) 近似： \[ \tilde{X}_M := \arg\min_{\text{rk}(A) \leq M} \|X - A\|_2 \;\;\in \mathbb{R}^{D \times N} \tag{10.50} \] 其中，\(\|\cdot\|_2\) 是在公式 (4.93) 中定义的谱范数。Eckart–Young 定理指出，\(\tilde{X}_M\) 可以通过截断奇异值分解到前 \(M\) 个奇异值来获得。换句话说： \[ \tilde{X}_M = U_M \; \Sigma_M \; V_M^\top \;\;\in \mathbb{R}^{D \times N} \tag{10.51} \] 其中：

• \(U_M := [u_1, \ldots, u_M] \in \mathbb{R}^{D \times M}\)，
• \(V_M := [v_1, \ldots, v_M] \in \mathbb{R}^{N \times M}\)，
• \(\Sigma_M \in \mathbb{R}^{M \times M}\) 是对角矩阵，对角线上包含了 \(X\) 的前 \(M\) 个最大奇异值。

10.4.2 实际问题

求解特征值和特征向量在许多需要矩阵分解的基础机器学习方法中也非常重要。理论上，如第 4.2 节所述，我们可以通过特征多项式的根来解特征值。然而，当矩阵大于 \(4 \times 4\) 时，这种方法就不可行了，因为这需要求解 5 次或更高次多项式的根。而 Abel–Ruffini 定理（Ruffini, 1799；Abel, 1826）表明，对于 5 次或更高次的多项式，不存在代数解。因此，在实际应用中，我们通过迭代方法来求解特征值或奇异值，这些方法已在所有现代线性代数库中实现。

在许多应用中（例如本章介绍的 PCA），我们只需要前几个特征向量。如果先计算完整分解，然后丢弃除前几个之外的所有特征向量，就会造成计算资源的浪费。事实证明，如果我们只对前几个特征向量（对应最大特征值）感兴趣，那么直接优化这些特征向量的迭代方法比完整的特征分解（或 SVD）更加高效。

在极端情况下，如果我们只需要第一个特征向量，那么一种简单的方法——幂迭代（power iteration） 就非常高效。幂迭代方法选择一个不在 \(S\) 零空间中的随机向量 \(x_0\)，并进行以下迭代： \[ x_{k+1} = \frac{S x_k}{\| S x_k \|}, \quad k = 0,1,\ldots \tag{10.52} \] 这意味着在每次迭代中，向量 \(x_k\) 都会被矩阵 \(S\) 乘一次，然后再进行归一化，即始终保持 \(\|x_k\| = 1\)。这列向量最终会收敛到 \(S\) 的最大特征值对应的特征向量。

最初的 Google PageRank 算法（Page 等人, 1999）就使用了这种算法来根据网页的超链接关系对网页进行排序。

10.5 高维中的 PCA

为了执行 PCA，我们需要计算数据协方差矩阵。在 \(D\) 维空间中，数据协方差矩阵是一个 \(D \times D\) 的矩阵。计算该矩阵的特征值和特征向量在计算上代价很高，因为其计算复杂度随 \(D\) 呈立方增长。因此，正如我们之前讨论的那样，在非常高维的情况下，PCA 是不可行的。

例如，如果我们的数据点 \(x_n\) 是具有 10,000 个像素的图像（例如 \(100 \times 100\) 的图像），那么我们就需要对一个 \(10,000 \times 10,000\) 的协方差矩阵进行特征分解。接下来，我们将针对 数据点数量远小于维度数量的情况（即 \(N \ll D\)） 提供一种解决方案。

假设我们有一个中心化后的数据集 \(x_1, \ldots, x_N\)，其中 \(x_n \in \mathbb{R}^D\)。那么数据协方差矩阵为： \[ S = \frac{1}{N} XX^\top \in \mathbb{R}^{D \times D}, \tag{10.53} \] 其中 \(X = [x_1, \ldots, x_N]\) 是一个 \(D \times N\) 的矩阵，每一列都是一个数据点。

现在我们假设 \(N \ll D\)，即数据点的数量小于数据的维度。如果没有重复的数据点，那么协方差矩阵 \(S\) 的秩为 \(N\)，因此它有 \(D-N+1\) 个特征值为 0。直观上，这意味着数据中存在冗余。接下来，我们将利用这一点，把原来的 \(D \times D\) 协方差矩阵转化为一个 \(N \times N\) 的协方差矩阵，其特征值全为正数。

在 PCA 中，我们最终得到如下的特征向量方程： \[ S b_m = \lambda_m b_m, \quad m = 1, \ldots, M, \tag{10.54} \] 其中 \(b_m\) 是主子空间的一个基向量。利用 (10.53) 中定义的 \(S\)，我们得到： \[ S b_m = \frac{1}{N} XX^\top b_m = \lambda_m b_m. \tag{10.55} \] 我们现在从左边乘上 \(X^\top \in \mathbb{R}^{N \times D}\)，得到： \[ \frac{1}{N} \underbrace{\boldsymbol{X}^\top \boldsymbol{X}}_{N \times N} \underbrace{\boldsymbol{X}^\top \boldsymbol{b}_m}_{=:\boldsymbol{c}_m} = \lambda_m \boldsymbol{X}^\top \boldsymbol{b}_m \iff \frac{1}{N} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{c}_m = \lambda_m \boldsymbol{c}_m, \tag{10.56} \] 令 \[ c_m := X^\top b_m, \] 那么我们得到一个新的特征值/特征向量方程： \[ X^\top X c_m = \lambda_m c_m. \] 这里 \(\lambda_m\) 仍然是特征值，这验证了我们在第 4.5.3 节中的结论：\(XX^\top \;\text{的非零特征值等于}\; X^\top X \;\text{的非零特征值。}\) 因此，对于矩阵 \(\frac{1}{N} X^\top X \in \mathbb{R}^{N \times N}\)，我们得到了与 \(\lambda_m\) 相关的特征向量： \[ c_m := X^\top b_m. \] 假设没有重复的数据点，该矩阵的秩为 \(N\)，因此它是可逆的。这也意味着 \(\frac{1}{N} X^\top X\) 拥有与数据协方差矩阵 \(S\) 相同的非零特征值。但它是一个 \(N \times N\) 的矩阵，因此比原始的 \(D \times D\) 协方差矩阵计算特征值和特征向量要高效得多。

现在我们已经得到了 \(\frac{1}{N} X^\top X\) 的特征向量，接下来我们需要恢复 原始协方差矩阵 \(S\) 的特征向量，因为 PCA 仍然需要它们。此时我们知道 \(\frac{1}{N} X^\top X\) 的特征向量。如果我们对特征值方程左乘 \(X\)，就得到：

\[ \frac{1}{N} \underbrace{\boldsymbol{X} \boldsymbol{X}^\top}_{S} \boldsymbol{X} \boldsymbol{c}_m = \lambda_m \boldsymbol{X} \boldsymbol{c}_m \tag{10.57} \]

这正是数据协方差矩阵 \(S\)，同时我们也得出 \(X c_m\) 是 \(S\) 的一个特征向量。

备注：如果我们要应用第 10.6 节中介绍的 PCA 算法，需要将 \(S\) 的特征向量 \(X c_m\) 归一化，使其范数为 1。

10.6 PCA 在实践中的关键步骤

下面我们将通过一个运行示例逐步讲解 PCA 的各个步骤，该示例在图 10.11 中进行了总结。我们有一个二维数据集（图 10.11(a)），希望使用 PCA 将其投影到一个一维子空间上。

1. 均值去除 首先通过计算数据集的均值 \(\mu\)，并从每个数据点中减去它，从而使数据居中。这保证了数据集的均值为 0（图 10.11(b)）。均值去除并不是严格必须的，但它能减少数值计算问题的风险。

2. 标准化 对每个维度 \(d = 1, \ldots, D\)，用该维度的标准差 \(\sigma_d\) 去除数据点的单位。这样，数据变得无量纲，并且在每个坐标轴方向上的方差均为 1，如图 10.11(c) 中的两个箭头所示。至此，数据的标准化过程完成。

3. 协方差矩阵的特征分解 计算数据的协方差矩阵，以及它的特征值和对应的特征向量。由于协方差矩阵是对称的，谱定理（定理 4.15）保证我们可以找到一组正交规范基（ONB）的特征向量。在图 10.11(d) 中，特征向量根据其对应特征值的大小进行了缩放。较长的向量张成了主子空间，我们记作 \(U\)。数据协方差矩阵被表示为椭圆。

4. 投影 我们可以将任意数据点 \(x_* \in \mathbb{R}^D\) 投影到主子空间上。为了正确完成这一过程，需要用训练数据在第 \(d\) 个维度的均值 \(\mu_d\) 和标准差 \(\sigma_d\) 对 \(x_*\) 进行标准化： \[ x_*^{(d)} \; \leftarrow \; \frac{x_*^{(d)} - \mu_d}{\sigma_d}, \quad d = 1, \ldots, D \tag{10.58} \] 其中 \(x_*^{(d)}\) 是 \(x_*\) 的第 \(d\) 个分量。投影得到： \[ \tilde{x}_* = BB^\top x_* \tag{10.59} \] 其在主子空间基底下的坐标为： \[ z_* = B^\top x_* \tag{10.60} \] 这里，\(B\) 是一个矩阵，列向量为与数据协方差矩阵最大特征值对应的特征向量。PCA 返回的是坐标式 (10.60)，而不是投影点 \(x_*\)。

由于我们对数据集进行了标准化，公式 (10.59) 只在标准化数据集的上下文中产生投影。若要得到原始数据空间（即标准化之前）的投影，就需要“反标准化”，即先乘以标准差，再加回均值： \[ \tilde{x}_*^{(d)} \; \leftarrow \; \tilde{x}_*^{(d)} \sigma_d + \mu_d, \quad d = 1, \ldots, D \tag{10.61} \] 图 10.11(f) 展示了原始数据空间中的投影效果。

例 10.4（MNIST 数字：重建） 接下来，我们将 PCA 应用于 MNIST 手写数字数据集。该数据集包含 60,000 个手写数字（0 到 9）的样本。每个数字是一个大小为 \(28 \times 28\) 的图像，即它包含 784 个像素，因此我们可以把该数据集中的每幅图像看作是 \(\mathbf{x} \in \mathbb{R}^{784}\) 的一个向量。部分数字的示例如图 10.3 所示。

为了说明问题，我们将 PCA 应用于 MNIST 数据集中的一个子集，并且专注于数字 “8”。我们使用了 5,389 张训练集中的 “8” 图像，并根据本章介绍的方法确定了主子空间。然后我们使用所学习到的投影矩阵去重建一组测试图像，结果如图 10.12 所示。图 10.12 的第一行显示了来自测试集的四个原始数字。后续的各行展示了这四个数字在使用维度分别为 1、10、100 和 500 的主子空间时的重建结果。我们可以看到，即使只使用一维的主子空间，也能得到一个勉强算得上不错的原始数字重建，但图像会模糊而且比较通用。随着主成分（PCs）数量的增加，重建结果逐渐变得更清晰，细节也逐渐被捕捉到。当使用 500 个主成分时，我们几乎得到了完美的重建。如果我们选择 784 个主成分，就可以毫无压缩损失地恢复出原始数字。

图 10.13 展示了平均平方重建误差： \[ \frac{1}{N}\sum_{n=1}^N \|x_n - \tilde{x}_n\|^2 = \sum_{i=M+1}^D \lambda_i, \tag{10.62} \] 它是主成分个数 \(M\) 的函数。我们可以看到，主成分的重要性下降得非常快，增加更多的主成分只能带来微小的收益。这与图 10.5 中的观察完全一致：我们发现投影数据的大部分方差仅由少数几个主成分捕获。当使用大约 550 个主成分时，我们基本上就能完整重建训练集中的 “8” 图像（在图像边缘的一些像素上没有变化，因为它们始终是黑色的）。

10.7 潜变量视角

Latent Variable Perspective

在前面的章节中，我们通过最大方差和投影的视角推导了 PCA，而没有涉及任何概率模型的概念。这样做的一个优点是：它让我们避开了概率论中复杂的数学细节；但另一方面，概率模型则能为我们提供更多的灵活性和有用的见解。更具体地说，概率模型可以：

• 提供一个似然函数，使我们能够显式地处理带噪声的观测数据（之前我们甚至没有讨论过噪声问题）；
• 允许我们通过边际似然来进行贝叶斯模型比较（参见第 8.6 节）；
• 将 PCA 看作是一个生成模型，从而能够模拟新的数据；
• 让我们能够直接建立与相关算法的联系；
• 通过应用贝叶斯定理来处理随机缺失的数据维度；
• 给予我们一个衡量新数据点“新颖性”的概念；
• 提供一种有原则的方式来扩展模型，例如推广为 PCA 混合模型；
• 使之前推导得到的 PCA 结果成为一个特殊情形；
• 允许通过边缘化模型参数来实现完整的贝叶斯推断。

通过引入一个连续值潜变量 \(z \in \mathbb{R}^M\)，我们可以将 PCA 表述为一个概率潜变量模型(probabilistic latent-variable model)。Tipping 和 Bishop (1999) 提出了这种潜变量模型，称为概率 PCA (Probabilistic PCA, PPCA)。PPCA 解决了上面大多数问题，而我们之前通过最大化投影空间方差或最小化重构误差得到的 PCA 解，实际上就是在无噪声设定下最大似然估计的一个特殊情况。

个人注：见10.8节：“在我们对 PPCA 的讨论中，我们假设模型的参数，即 \(B, \mu\) 和似然参数 \(\sigma^2\)，是已知的。”

10.7.1 生成过程与概率模型

在 PPCA 中，我们显式写出线性降维的概率模型。为此，我们假设存在一个连续潜变量 \(z \in \mathbb{R}^M\)，其先验分布为标准正态： \[ p(z) = \mathcal{N}(0, I), \] 并假设潜变量与观测数据 \(x\) 存在线性关系： \[ x = Bz + \mu + \epsilon \in \mathbb{R}^D, \tag{10.63} \] 其中 \(\epsilon \sim \mathcal{N}(0, \sigma^2 I)\) 表示高斯观测噪声，\(B \in \mathbb{R}^{D \times M}\) 和 \(\mu \in \mathbb{R}^D\) 描述了潜变量到观测变量的线性/仿射映射。于是，PPCA 将潜变量和观测变量联系起来： \[ p(x|z, B, \mu, \sigma^2) = \mathcal{N}(x|Bz + \mu, \sigma^2 I). \tag{10.64} \] 整体上，PPCA 诱导的生成过程如下(PPCA induces the following generative process:)： \[ z_n \sim \mathcal{N}(z|0, I), \tag{10.65} \]

\[ x_n | z_n \sim \mathcal{N}(x|Bz_n + \mu, \sigma^2 I). \tag{10.66} \] 要生成一个符合模型参数的数据点，可以按照祖先采样（ancestral sampling）的方式：首先从 \(p(z)\) 中采样一个潜变量 \(z_n\)，然后将 \(z_n\) 代入式 (10.64)，再采样得到对应的 \(x_n \sim p(x|z_n, B, \mu, \sigma^2)\)。这个生成过程使我们能够写出完整的概率模型（即所有随机变量的联合分布，参见第 8.4 节）：

\[ p(x, z | B, \mu, \sigma^2) = p(x|z, B, \mu, \sigma^2) p(z), \tag{10.67} \] 并且根据第 8.5 节的结果，可以立刻得到如图 10.14 所示的图模型。

备注：请注意连接潜在变量 \(z\) 和观测数据 \(x\) 的箭头方向：箭头是从 \(z\) 指向 \(x\) 的，这意味着 PPCA 模型假设高维观测 \(x\) 来源于一个低维的潜在原因 \(z\)。最终，我们显然对在给定一些观测值的情况下推断 \(z\) 感兴趣。为了实现这一点，我们将应用贝叶斯推断来“反转”箭头（以一种隐含的方式），从观测推到潜在变量。

例 10.5（使用潜在变量生成新数据）

图 10.15 展示了在使用二维主子空间时，PCA 找到的 MNIST 数字 “8” 的潜在坐标（蓝点）。我们可以在该潜在空间中查询任意一个向量 \(z_*\)，并生成一幅类似数字 “8” 的图像 \(\tilde{x}_* = Bz_*\)。我们展示了八个这样的生成图像及其在潜在空间中的对应表示。根据我们在潜在空间查询的位置不同，生成的图像也有所不同（形状、旋转、大小等）。如果我们在远离训练数据的地方进行查询，就会出现越来越多的伪影，例如左上角和右上角的数字。需要注意的是，这些生成图像的内在维度仅仅是 2。

10.7.2 似然函数与联合分布

利用第 6 章中的结果，我们通过对潜在变量 \(z\) 积分（见 8.4.3 节）来得到该概率模型的似然函数： $$ \[\begin{align} p(x \mid B, \mu, \sigma^2) & = \int p(x \mid z, B, \mu, \sigma^2)\, p(z)\, dz \tag{10.68a} \\ & = \int \mathcal{N}\!\big(x \mid Bz+\mu, \sigma^2 I\big)\, \mathcal{N}(z \mid 0, I)\, dz. \tag{10.68b} \end{align}\] \[ 根据 6.5 节，我们知道该积分的解是一个高斯分布，其均值为 \] _x[x] = z[Bz+] + = \[ 其协方差矩阵为 \] \[\begin{align} \mathbb{V}[x] & = \mathbb{V}_z[Bz+\mu] + \mathbb{V}_\epsilon[\epsilon] = \mathbb{V}_z[Bz] + \sigma^2 I \tag{10.70ab}\\ & = B\, \mathbb{V}_z[z]\, B^\top + \sigma^2 I = BB^\top + \sigma^2 I. \tag{10.70b} \end{align}\] $$ (10.68b) 中的似然函数可用于模型参数的极大似然（MLE）或最大后验（MAP）估计。

备注。 我们不能使用 (10.64) 中的条件分布来进行极大似然估计，因为它仍然依赖于潜在变量。极大似然（或 MAP）估计所需的似然函数应该只依赖于观测数据 \(x\) 和模型参数，而不能依赖于潜在变量。

根据 6.5 节，我们知道高斯随机变量 \(z\) 及其线性/仿射变换 \(x = Bz\) 是联合高斯分布的。我们已经知道其边缘分布为 \[ p(z) = \mathcal{N}(z \mid 0, I), \quad p(x) = \mathcal{N}(x \mid \mu, BB^\top + \sigma^2 I). \] 缺失的交叉协方差为 \[ \mathrm{Cov}[x, z] = \mathrm{Cov}_z[Bz+\mu, z] = B\, \mathrm{Cov}_z[z,z] = B. \tag{10.71} \] 因此，PPCA 的概率模型（即潜在变量与观测变量的联合分布）可明确写为 \[ p(x, z \mid B, \mu, \sigma^2) = \mathcal{N} \!\left( \begin{bmatrix} x \\ z \end{bmatrix} \,\middle|\, \begin{bmatrix} \mu \\ 0 \end{bmatrix}, \begin{bmatrix} BB^\top + \sigma^2 I & B \\ B^\top & I \end{bmatrix} \right), \tag{10.72} \] 其中均值向量的维度为 \(D+M\)，协方差矩阵的大小为 \((D+M) \times (D+M)\)。

10.7.3 后验分布

在 (10.72) 中的联合高斯分布 \(p(x,z \mid B, \mu, \sigma^2)\)，使我们能够直接利用第 6.5.1 节的高斯条件分布规则来确定后验分布 \(p(z \mid x)\)。给定观测值 \(x\)，潜在变量的后验分布为 \[ p(z \mid x) = \mathcal{N}(z \mid m, C), \tag{10.73} \] 其中 \[ m = B^\top (BB^\top + \sigma^2 I)^{-1}(x - \mu), \tag{10.74} \]

\[ C = I - B^\top (BB^\top + \sigma^2 I)^{-1}B. \tag{10.75} \] 需要注意的是，后验协方差并不依赖于观测数据 \(x\)。对于数据空间中的一个新观测 \(x_\ast\)，我们使用 (10.73) 来确定对应潜在变量 \(z_\ast\) 的后验分布。协方差矩阵 \(C\) 允许我们评估嵌入的置信度。若 \(C\) 的行列式很小（行列式度量体积），说明潜在嵌入 \(z_\ast\) 相对确定。若我们得到的后验分布 \(p(z_\ast \mid x_\ast)\) 方差很大，则可能遇到一个异常点。

个人注：虽然 \(C\) 公式上是常量，但当数据 \(x_*\) 远离模型学到的子空间或均值 \(\mu\) 时，后验分布的有效方差（比如沿特定方向的不确定性）会显得更大。

不过，我们可以进一步研究该后验分布，以理解在此后验下哪些其他数据点 \(x\) 是合理的。为此，我们利用 PPCA 的生成过程，该过程允许通过生成新的、在该后验下合理的数据来探索潜在变量的后验分布：

1. 从潜在变量的后验分布 (10.73) 中采样一个潜在变量 \[ z_\ast \sim p(z \mid x_\ast)。 \]
2. 从 (10.64) 中的条件分布采样一个重构向量 \[ \tilde{x}_\ast \sim p(x \mid z_\ast, B, \mu, \sigma^2)。 \] 如果我们重复上述过程多次，就可以探索潜在变量 \(z_\ast\) 的后验分布 (10.73) 及其对观测数据的含义。这个采样过程实际上就是在假设一些在该后验分布下合理的数据。

10.8 延伸阅读

我们从两个角度推导了 PCA：（a）在投影空间中最大化方差；（b）最小化平均重构误差。然而，PCA 还可以从不同的视角进行解释。让我们回顾一下已完成的工作：我们取一个高维数据 \(x \in \mathbb{R}^D\)，并使用一个矩阵 \(B^\top\) 来找到其低维表示 \(z \in \mathbb{R}^M\)。矩阵 \(B\) 的列向量是数据协方差矩阵 \(S\) 的特征向量，这些特征向量对应于最大的特征值。得到低维表示 \(z\) 之后，我们可以在原始数据空间中重构一个近似的高维版本： \[ x \approx \tilde{x} = Bz = BB^\top x \in \mathbb{R}^D, \] 其中 \(BB^\top\) 是一个投影矩阵。我们还可以把 PCA 看作是一种线性自编码器（如图 10.16 所示）。自编码器将数据 \(x_n \in \mathbb{R}^D\) 编码为一个代码 \(z_n \in \mathbb{R}^M\)，再将其解码为与 \(x_n\) 相似的 \(\tilde{x}_n\)。从数据到代码的映射称为编码器，从代码回到原始数据空间的映射称为解码器。如果我们考虑线性映射，即代码由\(z_n = B^\top x_n \in \mathbb{R}^M\)给出，并且我们希望最小化数据 \(x_n\) 与其重构 \(\tilde{x}_n = Bz_n\) 的平均平方误差（\(n=1,\dots,N\)），则得到： \[ \frac{1}{N} \sum_{n=1}^N \|x_n - \tilde{x}_n\|^2 = \frac{1}{N} \sum_{n=1}^N \|x_n - BB^\top x_n\|^2. \tag{10.76} \] 这意味着我们得到的目标函数与 10.3 节中 (10.29) 相同，因此通过最小化平方自编码损失，我们得到的就是 PCA 解。如果我们将 PCA 的线性映射替换为非线性映射，就得到非线性自编码器。一个典型的例子是深度自编码器，其中线性函数被深度神经网络替代。在这种情况下，编码器也被称为识别网络或推理网络，而解码器也被称为生成器。

PCA 的另一种解释与信息论相关。我们可以将代码看作是原始数据点的一个较小或压缩版本。当我们利用代码重构原始数据时，并不能完全恢复原始数据，而是得到一个稍有失真或带噪声的版本。这意味着我们的压缩是有损的。直观上，我们希望最大化原始数据与低维代码之间的相关性。更正式地，这与互信息相关。通过最大化互信息（信息论的核心概念，MacKay, 2003），我们将得到与 10.3 节中讨论的 PCA 相同的解。

在我们对 PPCA 的讨论中，我们假设模型的参数，即 \(B, \mu\) 和似然参数 \(\sigma^2\)，是已知的。Tipping 和 Bishop (1999) 描述了如何在 PPCA 框架下推导这些参数的极大似然估计（注意本章的符号体系与他们的不同）。当将 \(D\) 维数据投影到 \(M\) 维子空间时，极大似然参数为： \[ \mu_{\text{ML}} = \frac{1}{N} \sum_{n=1}^N x_n, \tag{10.77} \]

\[ B_{\text{ML}} = T(\Lambda - \sigma^2 I)^{\tfrac{1}{2}} R, \tag{10.78} \]

\[ \sigma^2_{\text{ML}} = \frac{1}{D-M} \sum_{j=M+1}^D \lambda_j, \tag{10.79} \] 其中，\(T \in \mathbb{R}^{D \times M}\) 包含数据协方差矩阵的 \(M\) 个特征向量，\(\Lambda = \text{diag}(\lambda_1,\dots,\lambda_M) \in \mathbb{R}^{M \times M}\) 是一个对角矩阵，其对角元为主轴对应的特征值，\(R \in \mathbb{R}^{M \times M}\) 是一个任意正交矩阵。极大似然解 \(B_{\text{ML}}\) 在任意正交变换下是唯一的，例如我们可以将 \(B_{\text{ML}}\) 右乘任意旋转矩阵 \(R\)，因此 (10.78) 本质上是一个奇异值分解（见 4.5 节）。该证明的概要由 Tipping 和 Bishop (1999) 给出。在 (10.77) 中给出的 \(\mu\) 的极大似然估计就是数据的样本均值。根据 (10.79)，观测噪声方差 \(\sigma^2\) 的极大似然估计是主子空间正交补中的平均方差，即无法通过前 \(M\) 个主成分捕获的剩余方差被当作观测噪声来处理。在无噪声极限 \(\sigma \to 0\) 时，PPCA 和 PCA 给出的解是相同的：由于数据协方差矩阵 \(S\) 是对称的，它可以被对角化（见 4.4 节），即存在一个由 \(S\) 的特征向量组成的矩阵 \(T\)，使得

\[ S = T \Lambda T^{-1}. \tag{10.80} \] 在 PPCA 模型中，数据协方差矩阵就是高斯似然 \(p(x \mid B, \mu, \sigma^2)\) 的协方差矩阵，即 \(BB^\top + \sigma^2 I\)（见 (10.70b)）。当 \(\sigma \to 0\) 时，我们得到 \(BB^\top\)，因此这个数据协方差必须等于 PCA 的数据协方差（及其在 (10.80) 中给出的分解），即 \[ \mathrm{Cov}[X] = T \Lambda T^{-1} = BB^\top \;\;\Longleftrightarrow\;\; B = T \Lambda^{\tfrac{1}{2}} R, \tag{10.81} \] 这就是说，在 \(\sigma = 0\) 的情况下，我们得到了 (10.78) 中的极大似然估计。由 (10.78) 和 (10.80) 可以清楚地看出，(P)PCA 实际上执行的是数据协方差矩阵的分解。在流式场景中（数据按顺序到达），推荐使用迭代的期望最大化（EM）算法来进行极大似然估计（Roweis, 1998）。为了确定潜在变量的维度（即代码的长度，或投影到的低维子空间的维度），Gavish 和 Donoho (2014) 提出了一种启发式方法：如果我们能够估计数据的噪声方差 \(\sigma^2\)，就应该舍弃所有小于

\[ \frac{4\sigma \sqrt{D}}{\sqrt{3}} \] 的奇异值。另一种方法是使用（嵌套的）交叉验证（见 8.6.1 节）或贝叶斯模型选择准则（见 8.6.2 节）来确定数据的内在维度的合理估计（Minka, 2001b）。

与第 9 章中关于线性回归的讨论类似，我们可以在模型参数上放置一个先验分布，并将其积分消去。这样做有两个好处：（a）避免参数的点估计以及相关问题（见 8.6 节）；（b）能够自动选择潜在空间的适当维度 \(M\)。在 Bishop (1999) 提出的贝叶斯 PCA 中，我们在模型参数上放置一个先验 \(p(\mu, B, \sigma^2)\)。该生成过程允许我们对模型参数进行积分，而不是对其条件化，从而解决了过拟合问题。由于这种积分在解析上是不可解的，Bishop (1999) 建议使用近似推断方法，例如 MCMC 或变分推断。关于这些近似推断技术的更多细节，可以参考 Gilks 等人 (1996) 和 Blei 等人 (2017) 的工作。

在 PPCA 中，我们考虑了线性模型 \[ p(x_n|z_n) = \mathcal{N}\big(x_n \mid Bz_n + \mu, \ \sigma^2 I \big) \] 其中先验分布为 \[ p(z_n) = \mathcal{N}(0, I), \] 即所有观测维度受到相同强度的噪声影响。如果我们允许每个观测维度 \(d\) 具有不同的方差 \(\sigma_d^2\)，就得到了因子分析（FA）（Spearman, 1904; Bartholomew 等, 2011）。这意味着 FA 比 PPCA 给似然函数提供了更多灵活性，但仍然要求数据由模型参数 \(B, \mu\) 来解释。然而，FA 不再允许存在封闭形式（closed-form）的最大似然解，因此需要使用迭代算法（如期望最大化算法 EM）来估计模型参数。在 PPCA 中，所有的平稳点都是全局最优解，而在 FA 中则不再成立。与 PPCA 相比，FA 在缩放数据时不变，但在旋转数据时会得到不同的解。另一个与 PCA 密切相关的算法是独立成分分析（ICA）（Hyvärinen 等, 2001）。 同样从潜变量的角度出发： \[ p(x_n|z_n) = \mathcal{N}\big(x_n \mid Bz_n + \mu, \ \sigma^2 I \big), \] 但这次我们将 \(z_n\) 的先验分布改为非高斯分布。ICA 可用于盲源分离。想象你在一个嘈杂的火车站，周围很多人都在说话。你的耳朵就像麦克风，它们会线性混合火车站内的不同语音信号。盲源分离的目标就是识别出混合信号中的各个组成部分。如前面在 PPCA 的最大似然估计中讨论过的，原始 PCA 解对任何旋转都是不变的。因此，PCA 可以识别信号所在的最佳低维子空间，但无法识别信号本身（Murphy, 2012）。ICA 通过修改潜在源的先验分布 \(p(z)\)，要求其为非高斯分布，从而解决了这个问题。关于 ICA 的更多细节可参考 Hyvärinen 等 (2001) 和 Murphy (2012) 的著作。PCA、因子分析和 ICA 是三种线性模型下的降维方法示例。Cunningham 和 Ghahramani (2015) 对线性降维方法提供了更广泛的综述。我们这里讨论的 (P)PCA 模型可以扩展出一些重要方向。在第 10.5 节中，我们解释了当输入维度 \(D\) 远大于数据点个数 \(N\) 时如何做 PCA。通过利用“PCA 可以通过计算（大量的）内积来实现”的洞见，这个思路可以进一步扩展到无限维特征的情况。核技巧（kernel trick）正是核 PCA 的基础，它允许我们在无限维特征之间隐式地计算内积（Schölkopf 等, 1998; Schölkopf 和 Smola, 2002）。还有一些非线性降维方法是从 PCA 推导出来的（Burges (2010) 提供了很好的综述）。我们在本节前面讨论的“PCA 的自编码器视角”可以将 PCA 表示为深度自编码器的特例。在深度自编码器中，编码器和解码器由多层前馈神经网络表示，而神经网络本身就是非线性映射。如果我们把这些神经网络的激活函数设为恒等函数，模型就等价于 PCA。

另一种非线性降维方法是 Lawrence (2005) 提出的高斯过程潜变量模型（GP-LVM）。GP-LVM 从我们用来推导 PPCA 的潜变量视角出发，将潜变量 \(z\) 与观测 \(x\) 之间的线性关系替换为高斯过程（GP）。在 GP-LVM 中，我们不是估计映射的参数（像在 PPCA 中那样），而是边缘化掉模型参数，并对潜变量 \(z\) 进行点估计。类似于贝叶斯 PCA，Titsias 和 Lawrence (2010) 提出的贝叶斯 GP-LVM 在潜变量 \(z\) 上维持一个分布，并使用近似推断方法将其积分掉。