14.7 独立成分分析和探索投影寻踪

多元数据通常可被看成是对一个隐含数据来源的多个间接的测量，而这个数据来源一般无法被直接测量。以下为一些例子：

• 教育性和心理学的测试使用对问卷的回答来测量被试主体隐含的智力水平和其他心理能力。
• 脑电图（EEG）扫描通过放置在头部不同位置的传感器记录到的电磁信号来间接地测量大脑不同部位的神经活动。
• 股票的交易价格不停地波动，反映了许多无法测量因子的影响，例如市场信心、外部影响、以及其他难以识别或测量的驱动因素。

因子分析（factor analysis）是统计学领域中发展出的一个用来识别这些隐含数据源的经典方法。因子分析模型通常会结合高斯分布的假设，这在某种程度上来说限制了模型的适用性。最近出现的独立成分分析（independent component analysis, ICA）是对因子分析的一个强劲的竞争方法，在下文中将会说明它的成功建立在隐含数据源的非高斯性质上。

14.7.1 潜变量和因子分析

奇异值分解 $\mathbf{X}=\mathbf{U}\mathbf{D}\mathbf{V}^T$（式 14.54）有潜变量（latent variable）的表达形式。记 $\mathbf{S}=\sqrt{N}\mathbf{U}$ 和 $\mathbf{A}^T=\mathbf{D}\mathbf{V}^T/\sqrt{N}$，则有 $\mathbf{X}=\mathbf{S}\mathbf{A}^T$，所以 $\mathbf{X}$ 的每一列都是 $\mathbf{S}$ 列向量的线性组合。由于 $\mathbf{U}$ 是正交矩阵，并且和之前一样假设 $\mathbf{X}$ 的每列均值为零（因此 $\mathbf{U}$ 的每列均值也为零），那么可推导出 $\mathbf{S}$ 的列向量均值为零、互不相关、并且方差为一。从随机变量的角度来看，我们可以将奇异值分解（SVD），或对应的主成分分析（PCA）表述为对一个潜变量模型的估计。

或简洁地写为 $X=\mathbf{A}S$。彼此相关的每个 $X_j$ 可被表达成一组不相关的、方差为一的随机变量 $S_\ell$ 的线性扩展。不过这并不是一个令人满意的结果，因为给定任意一个正交 $p\times p$ 矩阵 $\mathbf{R}$，都有：

而其中的 $\operatorname{Cov}(S^*)=\mathbf{R}\operatorname{Cov}(S)\mathbf{R}^T=\mathbf{I}$。因此这样的分解有很多，也就无法将某组特定的潜变量识别为唯一的隐含数据源。不过奇异值分解（SVD）的确有如下的性质：对于每个秩 $q<p$ 截断的分解都是对 $\mathbf{X}$ 的最优估计。

经典的因子分析模型主要由心理测量学（psychometrics）的研究者提出，它在一定程度上缓解了这些问题，可参考 Mardia et al. (1979)。当 $q<p$，一个因子分析模型的形式为：

或写为 $X=\mathbf{A}S+\varepsilon$。其中的 $S$ 是 $q<p$ 个隐含潜变量或因子的向量，$\mathbf{A}$ 是 $p\times q$ 的因子载荷（factor loading）矩阵，$\varepsilon_j$ 是互补相关均质为零的干扰项。可以将这些潜变量 $S_\ell$ 理解为所有 $X_j$ 中的变动的共同来源，而且决定了它们的相关性结构；同时互不相关的 $\varepsilon_j$ 是每个 $X_j$ 所独有的，包含了其他模型未解释的变动来源。一般会假设 $S_\ell$ 和 $\varepsilon_j$ 都为高斯分布的随即变量，然后用最大似然估计方法拟合模型。模型所有的参数都会出现在协方差矩阵中：

其中的 $\mathbf{D}_\varepsilon=\operatorname{diag}[\operatorname{Var}(\varepsilon_1),\dots,\operatorname{Var}(\varepsilon_p)]$。因为假设了 $S_j$ 为高斯分布并且互不相关，统计学意义上它们是独立的随机变量。所以可以将一组教育性的测试得分视为一些独立的隐含因子的体现，例如智力、自驱力、等等。矩阵 $\mathbf{A}$ 的列被称为 因子载荷（factor loading），会被用于对因子的命名和解释。

不过不幸的是识别问题（identifiability，式 14.79）仍未解决，其原因是给定任意一个 $q\times q$ 的正交矩阵 $\mathbf{R}$，式 14.81 中的 $\mathbf{A}$ 和 $\mathbf{A}\mathbf{R}^T$ 都是等价的。研究者可以去寻找更容易解释的不同旋转版本的因子，所以这使因子分析在应用中存在一定的主观性。这方面的问题导致一些研究者对因子分析持怀疑态度，也可能是它在当代的统计学中不常见的原因。我们不会在这里深入讨论具体细节，简略来说奇异值分解（SVD）在对式 14.81 的估计中起到了关键作用。例如，如果假设所有的 $\operatorname{Var}(\varepsilon_j)$ 都相等，那么 SVD 的前 $q$ 个成分就等同于由 $\mathbf{A}$ 所决定的子空间。

由于每个 $X_j$ 都有独自的扰动项 $\varepsilon_j$，所以因子分析可被看作是 $X_j$ 相关性结构的模型，而不是协方差结构的模型。这一点从式 14.81 的标准化协方差结构中可明显地看得出来（练习 14.14。虽然这并不是本节讨论的重点，但这一点是因子分析和主成分分析（PCA）的一个重要的区别。练习 14.15 介绍了一个简单的例子，其中因子分析和 PCA 的结果会由于上述两者的区别而产生巨大的差异。

14.7.2 独立成分分析

独立成分分析（independent component analysis，ICA） 模型的形式与式 14.78 完全一样，不过其中 $S_\ell$ 的假设为统计意义的独立（independent），而不是互不相关（uncorrelated）。直观来说，互不相关约束了一个多元概率分布的二阶交叉矩（协方差），而一般来说统计意义的独立约束了所有的交叉矩。这些额外的矩约束条件可以让我们能够唯一地识别出 $\mathbf{A}$ 中的元素。由于多元高斯分布可以被其二阶矩（协方差矩阵）所定义，所以它是一个例外情况，与之前（因子分析）一样，任意高斯分布的成分都在不考虑旋转变换下可以被确定。所以如果假设 $S_\ell$ 是独立的非高斯分布，式 14.78 和 14.80 中的识别问题也可以被解决。

本节将介绍式 14.78 中的完整 $p$ 成分模型，其中的 $S_\ell$ 是方差为一的独立变量；同时也有因子分析模型式 14.80 的 ICA 版本。我们的描述是基于 Hyvärinen and Oja (2000) 的综述文章。

我们希望可以找出 $X=\mathbf{A}S$ 中的混合矩阵 $\mathbf{A}$。不失一般性，可以假设 $X$ 已经被标准化为 $\operatorname{Cov}(X)=\mathbf{I}$；这一般是通过上述的奇异值分解（SVD）实现的。由于 $S$ 的协方差矩阵也是 $\mathbf{I}$，这就意味着 $\mathbf{A}$ 是正交矩阵。所以对 ICA 问题的求解就可以转化为寻找一个使得向量 $S=\mathbf{A}^TX$ 中的成分随机变量彼此独立（非高斯分布）的正交矩阵 $\mathbf{A}$。

图 14.37 展示了 ICA 在分离两个混合信号中的作用。这是经典的 鸡尾酒会问题（cocktail party problem） 的一个例子，不同的麦克风 $X_j$ 采集到了不同独立的声音源 $S_\ell$（音乐、不同人的说话、等等）的混合。通过利用原始信号源独立和非高斯的性质，ICA 可以实现 盲源分离（blind source separation）。

很多 ICA 的常见方法是基于熵的概念。一个密度函数为 $g(y)$ 的随即变量 $Y$ 的 微分熵（differential entropy） $H$ 的定义为：

信息论中有一个著名的结论：在方差相同的所有随机变量中，高斯分布变量的熵最大。最后，随机向量 $Y$ 成分之间的 互信息（mutual information） $I(Y)$ 是对其相关性的一个自然的测量：

$I(Y)$ 被称为 $Y$ 的密度函数 $g(y)$ 与其独立版本 $\prod_{j=1}^pg_j(y_j)$ 之间的 KL（Kullback–Leibler） 距离，其中的 $g_j(y_j)$ 为 $Y_j$ 的边际密度函数。 如果 $X$的协方差矩阵为 $\mathbf{I}$，并且 $Y=\mathbf{A}^TX$，$\mathbf{A}$ 为正交矩阵，则易证明：

寻找使 $I(Y)=I(\mathbf{A}^TX)$ 最小化的 $\mathbf{A}$，就是寻找使得其成分之间“最独立”的正交变换。从式 14.84 可看出，这等价于对 $Y$ 的各个成分的熵之和进行最小化，这也可理解成是对它们与高斯分布的偏离的最大化。

方便起见，Hyvärinen and Oja (2000) 没有使用熵 $H(Y_j$，而使用 负熵（negentropy） $J(Y_j)$，其定义为：

其中的 $Z_j$ 是一个与 $Y_j$ 方差相等的高斯分布随机变量。负熵取值非负，衡量了 $Y_j$ 与高斯分布的差别程度。作者提出了对负熵的简化近似，可被应用在数据上的计算和最优化。 图 14.37、图 14.38 和图 14.39 中的 ICA 结果都使用了这个近似：

其中 $G(u)=\frac{1}{a}\log\cosh(au)$，$1\leq a\leq2$。在应用与一个 $x_i$ 的样本上时，用样本的平均值来代替期望。这也是这几位作者提供的 FastICA 软件中的选项之一。更经典（但更不稳健）的测度是基于四阶矩，所以就是通过峰度（kurtosis）衡量与高斯分布的偏离程度。更多细节可参考 Hyvärinen and Oja (2000)。第 14.7.4 节会介绍他们寻找最优方向的近似牛顿算法。

综上所述，多元数据上的 ICA 是在寻找一系列正交投影，使得投影的数据尽可能地偏离高斯分布。对已经白化处理（whitened）的数据来说，这等同于寻找尽可能独立的组成成分。

ICA 本质上从一个因子分析结果出发，寻找可得到彼此独立的成分的那个旋转版本。从这个视角看，ICA 与心理测量学中传统的方差最大（varimax）和四次方最大（quartimax）旋转方法一样，只是另一个因子旋转方法。

示例：手写数字

我们再次回到第 14.5.1 节中通过 PCA 分析的手写数字“3”的识别问题。图 14.39 对比了（标准化的）主成分分析和独立成分分析的各自前五个成分，并都按同样的标准化单位展示。每个图都是从 256 维空间到二维空间上的投影。PCA 成分都看起来服从了联合高斯分布，而 ICA 成分都有一个长尾分布。这也并不出乎意料，因为 PCA 着重在方差，而 ICA 特异地寻找非高斯分布。由于所有的成分都已被标准化，所以图中并没有体现主成分的方差递减。

在每个 ICA 成分中都高亮标记了两组极值点以及一对中心点，图 14.40 展示了中心点对应的手写数字。这个图演示了每个成分的本性。例如第五个 ICA 成分捕捉的是数字 “3” 的长条收尾笔划。

示例：脑电图时序数据

ICA 已成为脑动力学（brain dynamics）研究中的一个重要工具。这里介绍的示例使用了 ICA 来分离多渠道脑电图数据中信号的不同成分（Onton and Makeig, 2006）。

被试者戴着一个嵌入了 100 个脑电图电极组成的网格，它们记录着头部不同位置的大脑活动。图 14.41¹ 的上图展示的是一个被试者在进行一个 30 分钟的普通“2-back”学习任务时这些电极中的九个的 15 秒钟的输出结果。被试者会以大概 1500 毫秒的时间间隔看到一个字母（B、H、J、C、F、或 K），然后需要通过一或两个按钮来反馈当前这个字母与两个之前的（上上个）字母相同还是不同。基于这些回答，被试者会得到或扣除分数，并且偶尔会得到奖励分数或扣除惩罚分数。时序数据展示出了脑电图信号中的空间上的相关性：邻近传感器的信号看起来非常相似。

这里有一个关键的假设：每个头部电极记录的信号是来自不同大脑皮层活动和非大脑皮层的其他来源的独立电势的混合信号。更多关于 ICA 在此领域的综述，可参考本章相关文献。

图 14.41 的下图展示了一部分 ICA 成分。不同颜色的头部图形用热力图表达了估计出的分离系数向量 $\hat{\alpha}_j$，标识了大脑活动的位置。对应的时序图展示了训练学习得出的 ICA 成分的活动形状。

例如，被试者在每次回答反馈信号（彩色垂直线）后会眨眼，造成了成分 IC1 和 IC2 的位置和伪信号。IC12 是与心率相关的伪信号。IC4 和 IC7 代表了前额 $\theta$ 频段的活动，出现在正确回答反馈后的一小段时间内。关于这个例子更细节的讨论以及 ICA 在脑电图模型中的应用，可参考 Onton and Makeig (2006)。

14.7.3 探索投影寻踪

Friedman and Tukey (1974) 提出了 探索投影寻踪（exploratory projection pursuit，EPP），一个用于高维数据可视化的图形化的探索方法。他们的看法是从高维数据到低维（一维或二维）的投影大多看起来是高斯分布的。而大家感兴趣的结构，例如聚类或长尾，都可能会由非高斯分布的投影展示出来。于是他们提出了对一些投影指标（project idnex）进行最优化，每个指标衡量了与高斯分布不同层面的背离。自从他们最早提出这个方法后，又出现了很多不同的改进方法（Huber, 1985；Friedman, 1987），以及很多不同的指标，包括了熵；一个交互性的图像软件 Xgobi 实现了这些方法（Swayne et al., 1991，现称为 GGobi）²。这些投影指标与之前的 $J(Y_j)$ 的形式完全相同，其中的 $Y_j=a_j^TX$ 是 $X$ 成分的一个标准化的线性组合。实际上，对交叉熵的一些近似和置换与投影寻踪（projection pursuit）中提出的指标完全一致。一般来说，投影寻踪中的方向 $a_j$ 没有正交性的约束。Friedman (1987) 先将数据转换使其在选定的投影下看起来是高斯分布的，然后在寻找后续的方向。尽管它们的出发点不同，但 ICA 和 EPP 至少在本节所描述的形式下是非常相似的。

14.7.4 ICA 的直接方法 😱

根据定义，独立成分的联合（乘积）密度函数为：

本节将介绍一个使用广义加性模型（第 9.1 节）直接对这个密度函数进行估计的方法。更完整的细节可以参考 Hastie and Tibshirani (2003)，R 中的 ProDenICA 包实现了这个方法，可从 CRAN 获取。

为了表现与高斯分布的背离程度，我们将每个 $f_j$ 表达为一个倾斜的（tilted）高斯密度函数：

其中的 $\phi$ 是标准高斯密度函数，$g_j$ 可令 $f_j$ 满足密度函数所要求的标准化条件。和之前一样假设 $X$ 已进行了白化处理（whitened），则观测数据 $X=\mathbf{A}S$ 的对数似然函数为：

我们希望对上式进行最大化，约束为 $\mathbf{A}$ 是正交矩阵，并且 $g_j$ 满足式 14.89 中密度函数的约束条件。在不对 $g_j$ 加以任何其他约束的条件下，模型 14.90 是过参数化的，所以我们对一个正则化版本进行最大化：

我们在式 14.91 中（为每个 $j$）及减去了两个惩罚项，这个想法来自于 Siverman (1986, Section 5.4.4)：

• 第一个惩罚项使得任意的解 $\hat{g}_j$ 都满足概率密度函数的约束条件 $\int\phi(t)e^{\hat{g}_j(t)}dt=1$。
• 第二项是对粗糙程度的惩罚，保证了解 $\hat{g}_j$ 是一个四次样条函数，每个观测值 $s_{ij}=a_j^Tx_i$ 都是它的节点。

进一步，可证明解的密度函数 $\hat{f}_j=\phi e^{\hat{g}_j}$ 都是均值为零、方差为一（练习 14.18）。随着我们逐渐提高 $\lambda_j$，这些解将逐步逼近于标准高斯分布 $\phi$。

我们用一种交替形式对式 14.91 进行最优化来拟合出函数 $g_j$ 和方向 $a_j$，如在算法 14.3 所述。

算法 14.3：乘积密度函数 ICA 算法（ProDenICA）

1. 初始化矩阵 $\mathbf{A}$（随机的高斯分布矩阵，然后进行正交化）。
2. 交替进行以下两步，直到 $\mathbf{A}$ 收敛。
   1. 给定 $\mathbf{A}$，（为每个 $j$ 分别地）针对 $g_j$ 对式 14.91 最优化。
   2. 给定 $g_j$，$j=1,\dots,p$，通过单步骤的不动点算法（fix point algorithm）寻找最优的 $\mathbf{A}$。

步骤 2.1 可看作为一个半参数化的密度估计，它可以通过对广义加性模型的一个新颖的应用来求解。方便其间，我们从 $p$ 个分别的最优化问题中选取一个：

虽然式 14.92 中的第二个积分可得出一个平滑的样条函数，第一个积分会带来一些问题，需要进行近似。我们创建一个密集的网格，共有 $L$ 个取值 $s_\ell^*$，递增幅度为 $\Delta$，网格覆盖着所有的观测样本点 $s_i$，然后计算每个网格区间中 $s_i$ 的个数：

我们通常选择 $L$ 为 1000，一般都是充足的。然后我们可以将式 14.92 近似为：

这个最后的表达式可视为与一个带惩罚项的泊松对数似然函数成比例，其输出变量为 $y_\ell^*/\Delta$，惩罚参数为 $\lambda/\Delta$，均值为 $\mu(s)=\phi(s)e^{g(s)}$。这是一个广义加性样条模型（generalized additive spline model，Hastie and Tibshirani, 1990; Efron and Tibshirani, 1996），抵消（offset）项为 $\log\phi(s)$。模型的拟合可以使用牛顿算法，计算量为 $O(L)$。虽然理论上使用的是四次样条，不过我们在实践中发现三次样条就已经可以了。有 $p$ 个调整参数 $\lambda_j$ 需要被确定，在实践中我们令它们取值相等，然后通过有效自由度 $df(\lambda)$ 来指定函数的平滑程度。我们的软件用 $df=5$ 作为默认值。

算法 14.3 的步骤 2.2 需要固定 $\hat{g}_j$ 后针对 $\mathbf{A}$ 对式 14.91 进行最优化。 式子中只有求和中的第一项与 $\mathbf{A}$ 有关，而且由于 $\mathbf{A}$ 是正交的，所有涉及 $\phi$ 的项的求和也不依赖于 $\mathbf{A}$（练习 14.19）。 因此我们只需要对下式最大化：

$C(\mathbf{A})$ 是拟合密度函数和高斯密度函数的对数似然函数的比值，并可被看作是对负熵（14.86）的一个估计，与式 14.87 一样每个 $\hat{g}_j$ 是一个对比函数。步骤 2.1 中的不动点更新是一个修改后的牛顿更新步骤（练习 14.20）。

1. 对每个 $j$，更新： $$a_j \leftarrow \operatorname{E}\{ X\hat{g}'_j(a_j^T X) - \operatorname{E}[\hat{g}''_j(a_j^T X)]a_j\} \tag{14.96}$$ 其中的 $\operatorname{E}$ 代表着对样本 $x_i$ 的期望。由于 $\hat{g}_j$ 是一个拟合的四次（或三次）样条函数，一阶和二阶导数已经可得。
2. 使用对陈平方根转换 $(\mathbf{A}\mathbf{A}^T)^{-\frac{1}{2}}\mathbf{A}$ 对 $\mathbf{A}$ 进行正交化。如果 $\mathbf{A}=\mathbf{U}\mathbf{D}\mathbf{V}^T$ 是 $\mathbf{A}$ 的一个奇异值分解（SVD），则容易推导出更新步骤为 $\mathbf{A}\leftarrow\mathbf{U}\mathbf{V}^T$。

例：模拟数据

图 14.42 展示了在一个模拟数据上 ProDenICA、FastICA、和另一个半参数化的方法 KernelICA（Bach and Jordan, 2002）的结果。左图展示了 18 个用作对比的基础概率分布。我们为每个概率分布生成了一对独立的成分变量（$N=1024$），以及一个随机的 $\mathbb{R}^2$ 上的取值在 1 到 2 之间³的混合矩阵。我们使用了 FastICA 在 R 中的实现，负熵的准则函数 14.87，以及 ProDenICA。我们使用了 KernelICA 作者的 MATLAB 代码⁴。由于最优化准则函数是非凸的，我们在每个方法中使用了五个随机的起始点。每个算法都给出了一个正交的混合矩阵 $\mathbf{A}$（样本数据已被白化），可以将这些矩阵与实际真实的正交化混合矩阵 $\mathbf{A}_0$ 进行比较。我们使用了 Amari 指标（Bach and Jordan, 2002）作为两个框架之间相似性的测量：

其中 $r_{ij}=(\mathbf{A}_0\mathbf{A}^{-1})_{ij}$。图 14.42 的右图对比了真实值和估计混合矩阵之间的平均 Amari 指标（对数尺度下）。 ProDenICA 在所有的场景中都优于 FastICA 和 KernelICA，并且在大多数的混合模拟中有明显的优势。

本节练习

练习 14.14

Pre- and post-multiply equation (14.81) by a diagonal matrix containing the inverse variances of the Xj . Hence obtain an equivalent decomposition for the correlation matrix, in the sense that a simple scaling is applied to the matrix A.

练习 14.15

Generate 200 observations of three variates X1 , X2 , X3 according $$\begin{align} X_1 & \sim Z_1 \\ X_2 &= X_1 + 0.001 \cdot Z_2 \\ X_3 &= 10 \cdot Z_3 \end{align}\tag{14.117}$$ where Z1 , Z2 , Z3 are independent standard normal variates. Compute the leading principal component and factor analysis directions. Hence show that the leading principal component aligns itself in the maximal variance direction X3 , while the leading factor essentially ignores the uncorrelated component X3 , and picks up the correlated component X2 + X1 (Geoffrey Hinton, personal communication).

练习 14.18

Consider the regularized log-likelihood for the density estimation problem arising in ICA,

The solution ĝ is a quartic smoothing spline, and can be written as ĝ(s) = q̂(s) + q̂⊥ (s), where q is a quadratic function (in the null space of the penalty). Let q(s) = θ0 + θ1 s + θ2 s2 . By examining the stationarity condi- tions for θ̂k , k = 1, 2, 3, show that the solution fˆ = φeĝ is a density, and has R ′′mean2 zero and variance one. If we used a second-derivative penalty {g (t)} (t)dt instead, what simple modification could we make to the problem to maintain the three moment conditions?

练习 14.19

If A is p × p orthogonal, show that the first term in (14.91) on page 567

with aj the jth column of A, does not depend on A.

练习 14.20

Fixed point algorithm for ICA (Hyvärinen et al., 2001). Consider maximizing C(a) = E{g(aT X)} with respect to a, with ||a|| = 1 and Cov(X) = I. Use a Lagrange multiplier to enforce the norm constraint, and write down the first two derivatives of the modified criterion. Use the approximation

to show that the Newton update can be written as the fixed-point update (14.96).