14.5 主成分、主曲线和主曲面

在第 3.4.1 节曾介绍过主成分，是为了阐释岭回归的收缩原理。主成分是数据的相互不相关、按方差排序的一系列投影。在下一节中我们将主成分表述为对 $N$ 个样本点 $x_i\in {\mathbb{R}}^p$ 的线性流形近似。然后在第 14.5.2 节介绍一些分线性的推广。第 14.9 节将介绍其他一些近期提出的非线性流形近似。

14.5.1 主成分

${\mathbb{R}}^p$ 上一组数据的 主成分（principal component） 为这组数据所有的秩 $q\le p$ 提供了一系列最优的线性近似。

记观测样本为 $x_1$、$x_2$、……、$x_N$，考虑表述这些样本的一个秩为 $q$ 的线性模型

其中 $\mu$ 是 ${\mathbb{R}}^p$ 上的一个位置向量，${\mathbf{V}}_q$ 是由 $q$ 个正交单位向量作为列向量的 $p\times q$ 矩阵，$\lambda$ 是一个长度为 $q$ 的参数向量。这是秩为 $q$ 的一个仿射（线性）超平面的参数化的表述。

图 14.20 和图 14.21 分别演示了 $q=1$ 和 $q=2$ 的情况。在数据上拟合这个模型，也就是对 重构误差（reconstruction error） 进行最小化：

我们可对参数 $\mu$ 和 ${\lambda }_i$ 进行部分最优化（练习 14.7），得出结果：

这样我们就只需要再对正交矩阵 ${\mathbf{V}}_q$ 求解：

为了计算简便，我们可以假设 $\overline{x}=0$（或者也可以用中心化的版本 ${\tilde{x}}_i=x_i-\overline{x}$ 代替原来观测的特征变量）。这个 $p\times p$ 的矩阵 ${\mathbf{H}}_q={\mathbf{V}}_q{\mathbf{V}}_q^T$ 是一个 投影矩阵（projection matrix），它把每个样本点 $x_i$ 映射到对应的秩为 $q$ 的重构点 ${\mathbf{H}}_q{x}_i$ 上，即为 $x_i$ 在由 ${\mathbf{V}}_q$ 列向量张成的子空间上的正交投影。这个解可以下的方式表达。将（中心化的）样本输入变量按行叠加成一个 $N\times p$ 的特征变量矩阵 $X$。我们可以对 $\mathbf{X}$ 进行 奇异值分解（singular value decomposition）：

这是数值分析中一个常见的的分解，已有一些计算这个分解的算法（例如 Golub and Van Loan, 1983）。其中的 $\mathbf{U}$ 是一个 $N\times p$ 的正交矩阵（即 ${\mathbf{U}}^T\mathbf{U}={\mathbf{I}}_p$），它的列向量 ${\mathbf{u}}_j$ 被称为 左奇异向量（left singular vector）；\mathbf{V}$是一个$p\times p$的正交矩阵（{\mathbf{V}}^T\mathbf{V}={\mathbf{I}}_p），它的列向量$\mathbf{v}_j$ 被称为 右奇异向量（right singular vector）；\mathbf{D} 是一个 $p\times p$ 的对角矩阵，对角线上的元素 $d_1\ge d_2\ge \cdots \ge d_p\ge 0$ 被称为 奇异值（singular value）。给定一个秩 $q$，式 14.53 的解 ${\mathbf{V}}_q$ 由 $\mathbf{V}$ 的前 $q$ 列组成。矩阵 $\mathbf{U}\mathbf{D}$ 的列被称为 $\mathbf{X}$ 的主成分（见第 3.5.1 节）。前 $q$ 个主成分即给出了式 14.52 中的 $N$ 个 ${\hat{\lambda }}_i$ 最优解（${\hat{\lambda }}_i$ 对应着 $N\times q$ 的矩阵 ${\mathbf{U}}_q{\mathbf{D}}_q$ 的 $N$ 个行向量）。

图 14.20 在 ${\mathbb{R}}^2$ 上展示了一维的主成分直线。每个数据点 $x_i$ 都在这条直线上有一个对应的最近的点，即为 $u_{i1}{d}_1{v}_1$。其中的 $v_1$ 是这条直线的方向向量，${\hat{\lambda }}_i=u_{i1}{d}_1$ 是这个点沿这条直线到原点的距离。

类似地，图 14.21 展示了半球数据在二维主成分平面上的拟合（左图）。右图展示了数据样本在前两个主成分空间上的投影。这个投影是上一节中的自组织映射方法中初始化配置的基础。这个方法比较成功地分隔了聚类簇。由于半球是非线性的，所以非线性的投影可能效果会更好，这会在下一节展开介绍。

主成分也有很多其他良好的性质，例如在所有的特征线型组合中 $\mathbf{X}{v}_1$ 的方差最高；而在满足参数向量（$v_2$）正交于 $v_1$ 的所有线型组合中 $\mathbf{X}{v}_2$ 的方差最高；以此类推。

例：手写数字

主成分在降维和压缩中是一个很有用的工具。我们在第一章中介绍过的手写数据数别问题中演示下它的作用。

图 14.22 展示了共 658 个样本中的 130 个手写数字“3”，每个图片是一个数字化的 $16\times 16$ 的灰度图。可看出这些图在书写风格、笔画粗细、和方向上都有明显的不同。我们将这些图片视为在 ${\mathbb{R}}^{256}$ 空间上的点 $x_i$，然后通过 SVD（式 14.54）计算它们的主成分。

图 14.23 展示了这些数据的前两个主成分。为这前两个主成分 $u_{i1}{d}_1$ 和 $u_{i2}{d}_2$，计算出 $5\mathrm{\%}$、$25\mathrm{\%}$、$50\mathrm{\%}$、$75\mathrm{\%}$、和 $95\mathrm{\%}$ 分位点，并用它们定义出叠加在图上的方形网格。圆圈点代表了那些离网格顶点最近的那些图片，其中距离测度主要考虑的是这两个投射坐标，但也给了正交空间中的其他成分一些权重。右图展示了对应着这些原型点的图片。这样我们可以直观地看出前两个主成分的性质。可见 $v_1$（水平方向坐标）主要解释了数字“3”最后一笔尾巴的长度，而 $v_2$（垂直方向坐标）主要解释了数字的笔画粗细。以式 14.49 中的参数化模型的角度，这个两个成分的模型的形式为：

这里将前两个主成分方向 $v_1$ 和 $v_2$ 用图片来表示。尽管可能一共有 256 个主成分，但大约前 50 个主成分占据了“3”中 $90\mathrm{\%}$ 的变化，前 12 个主成分占据了 $63\mathrm{\%}$。

将 $\mathbf{X}$ 的每列随机打乱可得到一个没有相关性的数据集，图 14.24 展示了与这个数据集对比的奇异值曲线。数字化图片的像素点本身就是相关的，而且由于这些图片都是同一个数字这个相关性就更强了。主成分的一个相对小的子集可被用作高维数据在低维空间上一组优秀的特征变量。

例：普洛克路斯忒斯变换和形状平均

图 14.25 中同时展示了橙色和绿色两组数据点。在这个例子中，这些点代表了两个手写字母“S”的数字化版本，它们是从一位“Suresh”的签名中抽取出的。图 14.26 展示了提取这两个字母“S”的整体签名图片（第三和第四图）。这些签名是通过触摸屏设备动态地捕捉的，这在现代的超市中很常见¹。每个“S”由 $N=96$ 个点来表示，我们将其各记为 $N\times 2$ 的矩阵 ${\mathbf{X}}_1$ 和 ${\mathbf{X}}_2$。在这些点之间存在着对应关系：${\mathbf{X}}_1$ 和 ${\mathbf{X}}_2$ 的第 $i$ 行都应该代表着两个“S”的同一个位置。以形态测量学（morphometrics）的语言来说，这些点代表了两个对象的 地标（landmark）。一般来说，如何寻找这些相应地标点的问题是困难而且比较主观的。在这个特殊的例子中，我们对每个签名的速度信号使用了 动态时间规整（dynamic time warping），不过在这里不会展开说明。

右图（图 14.25）对绿色的点使用了一个形变和旋转，使其尽可能地与橙色点匹配，即所谓的 普洛克路斯忒斯变换（Procrustes transformation）²（可参考 Mardia et al., 1979）。

考虑下面这个问题：

其中 ${\mathbf{X}}_1$ 和 ${\mathbf{X}}_2$ 都是对应点的 $N\times p$ 矩阵，$\mathbf{R}$ 是一个 $p\times p$ 的正交矩阵³，$\mu$ 是一个长度 $p$ 的位置坐标向量。这里的 $\Vert \mathbf{X}{\Vert }_F^2=\mathrm{trace}({\mathbf{X}}^T\mathbf{X})$ 是平方 弗罗比尼乌斯矩阵范数（Frobenius matrix norm）。

令 ${\overline{x}}_1$ 和 ${\overline{x}}_2$ 分别为矩阵的列均值向量，并且 ${\tilde{\mathbf{X}}}_1$ 和 ${\tilde{\mathbf{X}}}_2$ 分别是这两个矩阵被移除了均值的版本。进行一个 SVD 分解：${\tilde{\mathbf{X}}}_1^T{\tilde{\mathbf{X}}}_1=\mathbf{U}\mathbf{D}{\mathbf{V}}^T$。则式 14.56 的解为（练习 14.8）：

而这个最小距离被称为 普洛克路斯忒斯距离（Procrustes distance）。从解的形式可见，我们可以将每个矩阵对它的列中心向量进行去中心化，然后可完全地忽略掉位置向量。以下部分将遵照这样的假设。

普洛克路斯忒斯尺度距离（Procrustes distance with scaling） 是求解一个稍微更一般性的问题：

其中 $\beta >0$ 是一个正数。$\mathbf{R}$ 的解与之前一样，而 $\hat{\beta }=\mathrm{trace}(D)/\Vert {\mathbf{X}}_1{\Vert }_F^2$。

与普洛克路斯忒斯距离相关联的是 $L$ 个形状的 普洛克路斯忒斯平均（Procrustes average），它是对下式求解：

也就是寻找一个与所有的形状以平方普洛克路斯忒斯距离最接近的形状 $\mathbf{M}$。它可通过一个简单的交替算法求解：

0. 进行初始化（例如）：$\mathbf{M}={\mathbf{X}}_1$。
1. 固定 $\mathbf{M}$，求解 $L$ 个普洛克路斯忒斯旋转问题，得到 ${\mathbf{X}}_{\ell }^{\prime }\leftarrow \mathbf{X}{\hat{\mathbf{R}}}_{\ell }$。
2. 令 $\mathbf{M}\to \frac{1}{L}\sum _{\ell =1}^L{\mathbf{X}}_{\ell }^{\prime }$
   重复进行步骤 1 和 2，直到准则函数（式 14.59）收敛。

图 14.26 展示了三个形状的简单示例。注意我们仅能够找到对旋转变化下等价的解；或者可以通过一个约束来达到唯一性，例如要求 $\mathbf{M}$ 是上三角矩阵。在定义式 14.59 中可很容易地加入尺度缩放（scaling）；见练习 14.9。

更一般地，可以通过下式定义一组形状的 仿射不变平均（affine-invariant average）：

其中 ${\mathbf{A}}_{\ell }$ 为任意的 $p\times p$ 非奇异矩阵。这里需要通过一个标准化条件来规避平凡解（trivial solution），例如 ${\mathbf{M}}^T\mathbf{M}=\mathbf{I}$。这个解不需要进行迭代就可以计算出来（练习 14.10）。

1. 令 ${\mathbf{H}}_{\ell }={\mathbf{X}}_{\ell }({\mathbf{X}}_{\ell }^T{\mathbf{X}}_{\ell }{)}^{-1}{\mathbf{X}}_{\ell }^T$ 是被 ${\mathbf{X}}_{\ell }$ 所定义的秩为 $p$ 的投影矩阵。
2. $\mathbf{M}$ 则是由 $\overline{\mathbf{H}}=\frac{1}{L}\sum _{\ell =1}^L{\mathbf{H}}_{\ell }$ 的 $p$ 个最大的特征向量组成的 $N\times p$ 矩阵。

14.5.2 主曲线和主曲面

主曲线（principal curve） 是主成分（直）线的一般化，为 ${\mathbb{R}}^p$ 上的一组数据点提供了一个平滑的一维曲线近似。主曲面（principal surface） 是更进一步的一般化，提供了一个二维或更高维的曲面流形（manifold）近似。

我们先定义一些随机变量 $X\in {\mathbb{R}}^p$ 的主曲线，然后再考虑有限数据样本的情况。令 $f(\lambda )$ 为一个 ${\mathbb{R}}^p$ 上参数化的平滑曲线。因此 $f(\lambda )$ 是取值为 $p$ 个坐标的向量函数，每个（坐标）都是一个有单个（输入）参数 $\lambda$ 的平滑函数。举个例子，参数 $\lambda$ 可以被选择为与某个固定原点沿着曲线的弧长（arc-length）。对每个数据点取值 $x$，记 ${\lambda }_f(x)为曲线上距离$x$最近的点（的弧长）。那么满足下式的$f(\lambda)$则被称为随机向量$X$ 分布的主成分曲线：

也就是说 $f(\lambda )$ 是所有投射到 $\lambda$ 处的数据点的平均，可理解为这个投射点 $f(\lambda )$ “负责”了这些数据点。这个性质也被称为 自洽性（self-consistency）。尽管在实际问题中，连续的多元随机分布有无穷多个主曲线（Duchamp and Stuetzle, 1996），我们主要关注的是平滑的主曲线。图 14.27 展示了一个主曲线例子。

主点（principal point） 是一个重要的相关概念。假设有一组 $k$ 个原型点，并且随机分布的支撑集（support）上可确定每个点距离最近的原型点，即对这个点“负责”的原型点。这就会将特征空间划分成了所谓的 沃罗诺伊分区图（Voronoi region）。使随机变量 $X$ 到其对应原型点的预期距离最小的那组 $k$ 个原型点被称为这个随机分布的主点。每个主点都是自洽的，即主点是它对应的沃罗诺伊分区中 $X$ 的均值。例如当 $k=1$ 时，一个圆形的（多维）正态分布的主点是它的均值向量；当 $k=2$ 时，主点是在穿过均值向量的射线上位置对称的两个点。随机分布的主点可类比于 K 均值聚类得出的中心点。而主曲线可看成是 $k=\mathrm{\infty }$ 个主点，不过需要约束这些点形成一个平滑的曲线，类似于自组织映射模型是将 K 均值聚类中心点约束在一个平滑流形上。

为得出一个随机分布的主曲线 $f(\lambda )$，记它的坐标函数为 $f(\lambda )=[f_1(\lambda ),f_2(\lambda ),\dots ,f_p(\lambda )]$，记 $X^T=(X_1,X_2,\dots ,X_p)$。则考虑以下的交替步骤：

第一步固定 $\lambda$（函数），然后按自洽性的约束（式 14.61）更新主曲线函数。第二步固定主曲线，然后为每个数据点找出主曲线上距离它最近的点。在有限的数据集上，主曲线算法以线型主成分作为初始化，交替进行式 14.62 中的两个步骤，直到结果收敛。在步骤（a）中会使用散点图平滑器模型来估计条件期望，在每个 $X_j$ 处的平滑结果是（邻域上）弧长 $\hat{\lambda }(X)$ 的函数；在步骤（b）中，要为每个观测数据点分别做投影运算。证明普遍的收敛性比较困难，但如果在散点图平滑器中使用的是最小二乘拟合，则可证明这个计算过程收敛到第一个线型主成分，并且等价于寻找一个矩阵的最大特征量的 幂迭代方法（power method）。

主曲面与主曲线的形式完全一样，不过只是有更高的维度。最常用的是二维主曲面，坐标函数为：

则上述步骤（a）中的估计需要使用二维曲面平滑器得出。很少会使用高于二维的主曲面，因为在可视化上没有多少帮助，并且高维上的平滑器也比较难处理。

图 14.28 展示了半球数据拟合的主曲面的结果。右图中为数据点对估计出的非线性坐标 ${\hat{\lambda }}_1(x_i)$ 和 ${\hat{\lambda }}_2(x_i)$ 的函数。图中类别的分隔很清楚。

主曲面与自组织映射非常相似。如果在估计每个坐标函数 ${\text{f}}_j({\lambda }_1,{\lambda }_2)$ 时使用一个核曲面平滑器，那么这将与批处理版本的自组织映射（式 14.48）有一样的表达式。自组织映射中的权重 $w_k$ 即为核的权重。不过它们也有一处差异：主曲面为每个数据点 $x_i$ 分别地估计原型点 $f({\lambda }_1(x_i),{\lambda }_2(x_i))$，而自组织映射中所有的数据点共同使用少数几个原型点。因此，只有当自组织映射的原型点增加到非常大的数量时，自组织映射和主曲面才会趋于一致。

两者也存在一个概念上的差异。主曲面以坐标函数的方式给出了对整个流形的平滑的参数化，而自组织映射是离散的，只给出了用于近似数据的原型点估计。 主曲面的平滑参数化保留了局部距离：在图 14.28 中可见，红色的簇要比绿色的和蓝色的更紧凑。在简单的示例中可看出估计的坐标函数本身包含了这些信息，参考练习 14.13。

14.5.3 谱聚类

传统的聚类方法，例如 K 均值，使用了圆形的或椭圆的度量来为数据点分组。因此当聚类簇是非凸集合时，例如 图 14.29 左上图中的同心环，这些方法不再有效。谱聚类（spectral cluster） 是这些标准聚类方法的一个推广，专为这些场景而设计。它与推广了多维尺度分析的局部多维尺度方法（第 14.9 节）有紧密的联系。

首先，有一个所有观测样本点之间两两相似度 $s_{i{i}^{\prime }}\ge 0$ 的 $N\times N$ 矩阵。将观测样本点表述为一个无向的 相似度图（similarity graph） $G=⟨V,E⟩$。其中的 $N$ 个顶点（vertex）$v_i$ 代表了观测样本点；当两个顶点之间的相似度为正（或高于某个阈值）时，则这顶点对之间被一个边（edge）连接。这些边的权重为 $s_{i{i}^{\prime }}$。于时聚类问题就被表述成了一个图分区（graph partition）问题，将一些连接的成分定义为聚类簇。我们想要对图进行分区，使得不同组之间的边的权重比较低，而组内部的边权重比较高。谱聚类的思想就是构建一个表现了观测样本之间局部邻域关系的相似度图。

为了更明确的表述，令 $x_i\in {\mathbb{R}}^p$ 为一组 $N$ 个样本点，$d_{i{i}^{\prime }}$ 为 $x_i$ 和 $x_{i^{\prime }}$ 之间的欧式距离。我们会将径向基格拉姆矩阵（radial-kernel gram matrix）作为相似度矩阵，即 $s_{i{i}^{\prime }}=\exp (-d_{i{i}^{\prime }}^2/c)$，其中 $c$ 为一个尺度参数。

有很多可以定义相似度矩阵以及体现了局部性质的相似度图的方法。最普遍使用的是 互 K 近邻图（mutual K-nearest-neighbor graph）。定义 ${\mathcal{N}}_K$ 为邻近的点对（pair）的对称集合，即当点 $i$ 是点 $i^{\prime }$ 的 $K$ 最近邻之一时，点对 $(i,i^{\prime })$ 属于 ${\mathcal{N}}_K$，反之亦然。这样就可以将所有的（对称的）近邻点对连接了起来，并给它们的边赋予权重 $w_{i{i}^{\prime }}=s_{i{i}^{\prime }}$；否则，其他的边的权重则为 0。这样就将所有不在 ${\mathcal{N}}_K$ 中的点对之间的相似度设为了零，可得到这个修改后的相似度矩阵的图。

或者也可以使用包含了所有点对之间的边 $w_{i{i}^{\prime }}=s_{i{i}^{\prime }}$ 的全连接图，然后通过尺度参数 $c$ 来控制局部的性质。

相似度图中的边的权重的矩阵 $\mathbf{W}=\{w_{i{i}^{\prime }}\}$ 被称为 邻接矩阵（adjacency matrix）。顶点 $i$ 的 度（degree） 为 $g_i=\sum _{i^{\prime }}{w}_{i{i}^{\prime }}$，即与之连接的边的权重之和。令 $\mathbf{G}$ 表示对角元素为 $g_i$ 的对角线矩阵。

最后，图的 拉普拉斯矩阵（graph Laplacian） 的定义为：

这被称为 非标准化拉普拉斯矩阵（unnormalized graph Laplacian）；已有一些标准化的拉普拉斯矩阵，它们是针对节点度 $g_i$ 而进行的标准化，例如，$\tilde{\mathbf{L}}=\mathbf{I}-{\mathbf{G}}^{-1}\mathbf{W}$。

谱聚类是要找出 $\mathbf{L}$ 的 $m$ 个最小的特征值对应的 $m$ 个特征向量 ${\mathbf{Z}}_{N\times m}$（不包括常数特征向量的平凡解）。使用例如 K 均值的标准聚类方法，就可以对 $\mathbf{Z}$ 的行进行聚类，从而得出一个对原始数据点的聚类。

图 14.29 展示了一个例子。左上图是用不同颜色标记三个环形簇的共 450 个模拟数据点。很明显，K 均值聚类会难以识别外层的簇。我们使用了 10 个最近邻相似度图，并且在左下图中展示了图的拉普拉斯矩阵的第二和第三最小特征值对应的特征向量。右上图展示了最小的 15 个特征值。（左下图）展示出的那两个特征向量就已经可以识别出三个簇了，右下图中的特征向量矩阵 $\mathbf{Y}$ 每行的散点图可明显看出三个簇的分隔。在这些转换后的样本点上再使用例如 K 均值的聚类方法，就可以很轻松地识别出三个组。

谱聚类为什么在这里有效？对每个向量 $\mathbf{f}$：

式 14.64 说明了如果邻近的一对样本点，它们的坐标 $f_i$ 和 $f_{i^{\prime }}$ 比较接近，那么 ${\mathbf{f}}^T\mathbf{L}\mathbf{f}$ 的取值就可能比较小。

由于任意图都满足 ${\mathbf{1}}^T\mathbf{L}\mathbf{1}=0$，所以常数向量就是特征值为零的特征向量平凡解。而这里隐含了一个结论：如果这个图是连通的⁴，则这是它唯一的零特征向量（练习 14.21）。对这个结论进行推广，则可见对于一个有 $m$ 个元件（connected component）的图，可重新排列节点的顺序使得 $\mathbf{L}$ 称为一个由每个元件组成的分块对角矩阵。那么 $\mathbf{L}$ 有 $m$ 个零特征值的特征向量，而元件的指示向量张成了零特征值的特征向量空间。在实际操作中会有强连接和弱连接，所以小特征值可近似为零特征值。

谱聚类是寻找非凸簇的一个重要的方法。如果使用了标准化的拉普拉斯矩阵，那么可以从另一个角度来理解这个方法。定义矩阵 $\mathbf{P}={\mathbf{G}}^{-1}\mathbf{W}$，则可按转移矩阵（transition probability matrix）$\mathbf{P}$ 在图上进行随机行走。那么随机行走会更倾向于在几个组的节点内部移动，而很少在组之间跨越，这也就是谱聚类所得出的节点分组。

在实践中应用谱聚类时，会有一些待解决的问题。首先是要选择相似度图的类型，例如是全连通的图还是 K 最近邻的图，以及相对应的参数，例如最近邻的数量 $k$ 或核的尺度参数 $c$。其次也需要选择从 $\mathbf{L}$ 中提取特征向量的数量，最后与所有聚类方法一样，还需要选择簇的个数。在图 14.29 的例子中，在 $k\in [5,200]$ 的区间上我们得到了不错的结果，其中的上限值 $200$ 对应了一个全连通的图。在 $k<5$ 时的结果不理想。 从图 14.29 的右上图中可见最小的三个特征值与其他特征值之间没有明显的区分。因此并无法明确应该选择多少个特征向量。

14.5.4 核主成分

谱聚类可与 核主成分（kernel principal component） 相关联，后者是线型主成分的一个非线性的版本。标准的线型主成分（PCA）是从协方差矩阵的特征向量得出的，并且可给出数据样本中方差最大的方向。核主成分（Schölkopf et al., 1999）扩展了 PCA 的范畴，模仿出假设对特征进行非线性转换的扩展后可能得到的特征，然后在这个转换过的特征空间上适用 PCA。

第 18.5.2 节中说明了一个数据样本矩阵 $\mathbf{X}$ 的主成分变量 $\mathbf{Z}$ 可从内积（格拉姆，Gram）矩阵 $\mathbf{K}=\mathbf{X}{\mathbf{X}}^T$ 中计算得出。具体来说，计算下式中双重中心化版本的格拉姆矩阵的特征值分解：

其中 $\mathbf{M}=\mathbf{1}{\mathbf{1}}^T/N$，并且令 $\mathbf{Z}=\mathbf{U}\mathbf{D}$。练习 18.15 说明了如何计算新观测样本点在这个空间上的投影。

核主成分就是对这个过程的模拟，将核矩阵 $\mathbf{K}=\{K(x_i,x_{i^{\prime }})\}$ 看作是隐含特征的内积 $⟨\varphi (x_i),\varphi (x_{x^{\prime }}⟩$ 的矩阵，并计算它的特征向量。第 $m$ 个成分 ${\mathbf{z}}_m$（即 $\mathbf{Z}$ 的第 $m$ 列）中的元素（中心化处理后）可写为 $z_{im}=\sum _{j=1}^N{\alpha }_{jm}K(x_i,x_j)$，其中的 ${\alpha }_{jm}=u_{jm}/d_m$（练习 14.16）。

通过将 ${\mathbf{z}}_m$ 视为主成分函数 $g_m\in {\mathcal{H}}_K$ 在样本上的取值，我们可以更进一步理解核主成分的原理，其中 ${\mathcal{H}}_K$ 为核函数 $K$ 生成的再生核希尔伯特空间（第 5.8.1 节）。第一个主成分函数 $g_1$ 是下式的解：

其中的 ${\mathrm{Var}}_{\mathcal{T}}$ 代表着训练集 $\mathcal{T}$ 上的样本方差。范数的约束 $\Vert g_1{\Vert }_{{\mathcal{H}}_K}=1$ 控制了函数 $g_1$ 的大小和粗糙程度，而这些是由核函数 $K$ 决定的。与在回归问题中类似，可证明式 14.66 的解是有限维度的，其表达式为 $g_1(x)=\sum _{j=1}^N{c}_{j}K(x,x_j)$。练习 14.17 证明了这个解就是由上述的 $\alpha$ 所确定的：${\hat{c}}_j={\alpha }_{j1},j=1,\dots ,N$。第二个主成分函数的定义与之类似，只是有一个额外的约束 $⟨g_1,g_2{⟩}_{{\mathcal{H}}_K}=0$，以此类推⁵。

Schölkopf et al. (1999) 在手写数字分类问题中演示了用核主成分作为分类的特征变量，并且说明了当用这些特征变量取代线型主成分后可提升分类模型的表现。

请注意如果使用了下式的径向核函数：

则核矩阵 $\mathbf{K}$ 与谱聚类中的相似度矩阵 $\mathcal{S}$ 有同样的形式。而边权重矩阵 $\mathbf{W}$ 是一个对 $\mathbf{K}$ 局部化的版本，将所有非最近邻点对之间的相似度设置为零。

核主成分要找出 $\hat{\mathcal{K}}$ 最大特征值对应的特征向量；这就等价于找出下式矩阵最小特征值对应的特征向量：

这与拉普拉斯矩阵（式 14.63）几乎相同，区别在于 $\tilde{\mathbf{K}}$ 的中心化和 $\mathbf{G}$ 的对角线上是节点的度（degree）。

图 14.30 演示了在图 14.29 的例子中核主成分方法的表现结果。左上图中使用了 $c=2$ 的径向基函数，也是谱聚类中所使用的参数值。这时并没有将不同组分离出来，但是当 $c=10$ 时（右上图），第一个主成分可很好地分离出不同组。在左下图中，我们将谱聚类中的最近邻径向核矩阵 $\mathbf{W}$ 用于核主成分。在右下图中，我们将核矩阵本身作为相似度矩阵来构建谱聚类中的拉普拉斯矩阵（式 14.63）。在这两个情景中，它们的投影都无法分离出其中的两个组。调整参数 $c$ 也没有作用。

从这个例子中可看出核主成分对尺度和核的性质非常敏感。同时也可见核函数的最近邻截取对谱聚类也非常重要。

14.5.5 稀疏主成分

我们通常通过方向向量 $v_j$，也被称作 载荷（loading），来解释主成分分析，从而查看到起作用的变量。例如被加入到式 14.55 中的图片。如果载荷是稀疏的，对结果的解释通常会更容易一些。本节会简单讨论一些推导稀疏载荷下的主成分的方法。它们都基于套索（$L_1$）惩罚项。

首先我们有一个 $N\times p$ 的数据样本矩阵 $\mathbf{X}$，它的列向量已经被中心化。 所提出的方法或是聚焦于主成分的最大方差的性质，或是利用对重构误差的最小化。 Joliffe et al. (2003) 提出的 SCoTLASS 方法采取的是第一个路线，对下式求解：

绝对值形式的约束条件使得一部分载荷倾向等于零，因此 $v$ 会是稀疏的。以同样的方式来计算后续的稀疏主成分，额外加上主成分 $k$ 需要与前 $k-1$ 个主成分正交的约束。遗憾的是这个问题是非凸的，所以难以计算。

Zou et al. (2006) 从另一方面主成分的回归/重构性质入手，这与第 14.5.1 节 中的方法类似。令 $x_i$ 为 $\mathbf{X}$ 的第 $i$ 行。则对于单个主成分，其 稀疏主成分（sparse principal component） 方法是下式的解：

从这个表达式的细节可见：

• 如果 $\lambda$ 和 ${\lambda }_1$ 都为零，并且 $N>p$，易证明出 $v=\theta$ 而且就是最大的主成分的方向向量。
• 如果 $p\gg N$，除非 $\lambda >0$，否则解不一定唯一。给定任意的 $\lambda >0$ 和 ${\lambda }_1=0$，$v$ 的解与最大的主成分方向向量成比例。
• 对 $v$ 的第二个惩罚项导致了载荷中的稀疏性。

对于多个主成分，稀疏主成分方法对下式求解：

使得 ${\mathbf{\Theta }}^T\mathbf{\Theta }={\mathbf{I}}_K$。其中的 $\mathbf{V}$ 是一个 $p\times K$ 的列为 $v_k$ 的矩阵，$\mathbf{\Theta }$ 的维度也是 $p\times K$。

准则函数 14.71 在 $\mathbf{V}$ 和 $\mathbf{\Theta }$ 上不是联合凸性的，不过在每个参数固定的情况下它对零一个参数是凸性的⁶。在 $\mathbf{\Theta }$ 固定的情况下对 $\mathbf{V}$ 进行最小化，这等价于 $K$ 个弹性网络（elastic net）问题（第 18.4 节），可高效地求解。从另一方面，在 $\mathbf{V}$ 固定的情况下对 $\mathbf{\Theta }$ 进行最小化，是普洛克路斯忒斯（Procrustes）问题（式 14.56）的一个版本，可通过一个简单的奇异值分解（SVD）求解（练习 14.12。交替进行这两个步骤，直到结果收敛。

图 14.31 展示了来自 Sjöstrand et al. (2007) 的一个稀疏主成分（式 14.71）的示例。这项研究的样本是 569 位老年人，其中有胼胝体（corpus callosum，CC）的正中矢状面（mid-sagittal）的剖面图形状，以及一些与之相关的临床指标⁷。在这个例子中，对形状（shape）数据使用了主成分分析，这在形态测量学（morphometrics）中是一个常用的工具。在这种应用中，会沿着形状（shape）的周界（circumference）生成一些地标（landmark）点；图 14.32 展示出了一个例子。再用普洛克路斯忒斯（Procrustes）分析来校准（align）这些地标点，从而允许形状中存在旋转和缩放的形变（参考第 14.5.1 节）。主成分中使用的特征变量向量，就是把每个地标点的一系列坐标对放到同一个向量中。

在这个例子的分析中同时计算了标准主成分和稀疏主成分，并且分辨出了与不同临床指标有显著关联的对应成分。在图中同时叠加了对应着显著的主成分的形状变化（红色曲线），和胼胝体（CC）平均的形状（黑色曲线）。与行走速度缓慢相关联的胼胝体，会在连接着运动控制和大脑认知中心的区域比较狭窄（出现了萎缩）。与语言表达不流畅相关联的胼胝体，会在连接着声音、视觉和认知中心的区域比较狭窄。图中的稀疏主成分形状与样本均值形状的差异更微小，而且可能更突出了关键区分的位置。

本节练习

练习 14.7

推导第 14.5.1 节 中的式 14.51 和 14.52。并说明 $\hat{\mu }$ 不唯一，以及一系列等价的解。

练习 14.8

Derive the solution (14.57) to the Procrustes problem (14.56). Derive also the solution to the Procrustes problem with scaling (14.58).

练习 14.9

Write an algorithm to solve

Apply it to the three S’s, and compare the results to those shown in Fig- ure 14.26.

练习 14.10

Derive the solution to the affine-invariant average problem (14.60). Apply it to the three S’s, and compare the results to those computed in Exercise 14.9.

练习 14.12

Consider the sparse PCA criterion (14.71).

1. Show that with $\mathbf{\Theta }$ fixed, solving for $\mathbf{V}$ amounts to K separate elastic-net regression problems, with responses the K elements of ${\mathbf{\Theta }}^T{x}_i$
2. Show that with V fixed, solving for Θ amounts to a reduced-rank version of the Procrustes problem, which reduces to $$\begin{array}{}\text{(14.116)}& \underset{\mathbf{\Theta }}{max}\mathrm{trace}({\mathbf{\Theta }}^T\mathbf{M})\text{ subject to }{\mathbf{\Theta }}^T\mathbf{\Theta }={\mathbf{I}}_K\end{array}$$ 其中的 $\mathbf{M}$ 和 $\mathbf{\Theta }$ 都是 $p\times K$ 的矩阵，$K\le p$。如果 $\mathbf{M}=\mathbf{U}\mathbf{D}{\mathbf{Q}}^T$ 是矩阵 $\mathbf{M}$ 的奇异值分解（SVD），请证明最优解 $\mathbf{\Theta }=\mathbf{U}{\mathbf{Q}}^T$。

练习 14.13

Generate 200 data points with three features, lying close to a helix. In detail, define X 1 = cos(s) + 0.1 · Z 1 , X 2 = sin(s) + 0.1 · Z 2 , X 3 = s + 0.1 · Z 3 where s takes on 200 equally spaced values between 0 and 2π, and Z 1 , Z 2 , Z 3 are independent and have standard Gaussian distributions.

1. Fit a principal curve to the data and plot the estimated coordinate functions. Compare them to the underlying functions cos(s), sin(s) and s.
2. Fit a self-organizing map to the same data, and see if you can discover the helical shape of the original point cloud.

练习 14.16

Consider the kernel principal component procedure outlined in Section 14.5.4. Argue that the number M of principal components is equal to the rank of K, which is the number of non-zero elements in D. Show that the mth component z m (mth column of Z) can be written (up to centering) as z im = j=1 α jm K(x i , x j ), where α jm = u jm /d m . Show that the mapping of a new observation $x_0$ to the mth component is given by $z_{0m}=\sum _{j=1}^N{\alpha }_{jm}K(x_0,x_j)$。

练习 14.17

Show that with $g_1(x)=\sum _{j=1}^N{c}_{j}K(x,x_j)$, the solution to (14.66) is given by ĉ j = u j1 /d 1 , where u 1 is the first column of U in (14.65), and d 1 the first diagonal element of D. Show that the second and subsequent principal component functions are defined in a similar manner (hint: see Section 5.8.1.)

练习 14.21

Consider an undirected graph with non-negative edge weights w ii ′ and graph Laplacian L. Suppose there are m connected components A 1 , A 2 , . . . , A m in the graph. Show that there are m eigenvectors of L corresponding to eigenvalue zero, and the indicator vectors of these components I A 1 , I A 2 , . . . , I A m span the zero eigenspace.