随机梯度 MCMC 推断

〖摘要〗本文提出了一个采用蒙特卡洛方法的新框架，用于从小批量迭代学习的大规模数据集学习。当我们对步长进行退火时，将适量噪声添加到标准随机梯度优化算法中，其结果表明迭代将收敛到来自真实后验分布的样本。优化和贝叶斯后验采样之间的这种无缝过渡提供了针对过拟合的内在保护。我们还提出了一种后验统计量的实用蒙特卡洛估计方法，该方法监视 “采样阈值” 并在超过阈值后收集样本。我们基于自然梯度将该方法应用于高斯、逻辑斯谛回归和独立组份分析的混合模型。

〖原文〗 Welling, M. and Teh, Y.W. (2011) ‘Bayesian learning via stochastic gradient Langevin dynamics’, in Proceedings of the 28th international conference on machine learning (ICML-11), pp. 681–688.

1 引言

近年来，越来越多的超大规模机器学习数据集，范围从互联网流量和网络数据、计算机视觉、自然语言处理到生物信息学。现在，这些大规模数据推动了机器学习的越来越多的进步，这提供了学习大型复杂模型以解决许多有用的应用问题的机会。最近在大规模机器学习中取得的成功主要是基于优化的方法。虽然有专门为某些类型的模型设计的复杂算法，但最成功的一类算法是随机优化或 Robbins-Monro 算法。这些算法在每次迭代中处理小批数据，通过在代价函数中采用小梯度步骤来更新模型参数。这些算法通常在在线设置中运行，其中数据批次在处理后被丢弃，并且只执行一次数据传递，从而大大减少了内存需求。

大规模机器学习进展中相对 “落后” 的一类方法是贝叶斯方法。这部分与贝叶斯在线参数估计的负面结果有关 (Andrieu 等, 1999 ^[2])，但也与典型马尔可夫链蒙特卡罗 (MCMC) 算法的每次迭代都需要对整个数据集进行计算这一事实有关。尽管如此，贝叶斯方法在捕获学习参数的不确定性并避免过拟合方面的能力很有吸引力。可以说，对于大型数据集，几乎不会出现过拟合。或者，随着我们可以访问更大的数据集和更多的计算资源，我们开始对构建更复杂的模型感兴趣，因此始终需要量化参数不确定性的数量。

在本文中，我们提出了一种从大规模数据集进行贝叶斯学习的方法。我们的方法结合了随机优化似然的 Robbins-Monro 类型算法，以及将噪声注入参数更新的 Langevin 动力学，使得参数轨迹最终收敛到完整的后验分布，而不仅仅是最大后验峰值。生成的算法一开始类似于随机优化，然后自动过渡到使用 Langevin 动力学模拟后验样本的算法。

在 第 2 节 中，我们介绍了我们方法的两个组成部分：随机优化和 Langevin 动力学。 第 3 节 描述了我们的算法以及它如何收敛到后验分布。 第 4 节 描述了一种实用的方法来估计我们的算法何时从随机优化过渡到 Langevin 动力学。 第 5 节 在几个模型上演示了我们的算法， 第 6 节 总结。

2. 预备知识

令 $θ$ 表示参数向量，$p(θ)$ 是先验分布，$p(x|θ)$ 是给定 $θ$ 参数化时模型产生数据 $x$ 的似然（虽表示为 $p(x|\theta)$ 形式，但并不是概率分布）。一组 $N$ 个数据项 $X = \{x_i\}^N_{i=1}$ 的后验分布是：$p(θ|X) \propto p(θ) \prod^N_{i=1} p(x_i|θ)$。在优化文献中，先验对参数进行正则化，而似然项通常构成了要优化的代价函数，常见任务是找到最大后验参数 $\theta^*$。

一类流行的方法称为随机优化 (Robbins & Monro, 1951 ^[8])，操作如下。在每次迭代 $t$ 时，如果给定 $n$ 个项构成的数据子集 $X_t = \{x_{t1},\ldots , x_{tn}\}$，则参数更新如下：

$$\Delta \theta_t = \frac{\epsilon_t}{2} \left(\nabla \log p(\theta_t) + \frac{N}{n} \sum^n_{i=1} \nabla \log p(x_{ti}|\theta_t) \right) \tag{1}$$

其中 $\epsilon_t$ 是一系列步长。

随机优化的朴素想法是： 用在子集上计算的梯度近似整个数据集的真实梯度。在经过多次迭代后，多次随机抽取的子集接近于使用了整个数据集，并且由于使用子集而产生的梯度噪声，也会在迭代后被平均掉。对于子集梯度的近似足够准确的大型数据集来说，这会比直接使用整个数据集显著节省计算量。

为了确保收敛到局部最大值，除了其他技术假设外，一个主要要求是步长满足性质

$$\sum^\infty_{t=1} \epsilon_t = \infty \qquad \sum^\infty_{t=1} \epsilon^2_t < \infty \tag{2}$$

直观上，第一个约束确保参数将到达高概率区域，无论它被初始化到多远，而第二个约束确保参数将收敛到峰值而不是仅仅围绕它跳动。通常，步长 $\epsilon_t = a(b + t)^{-\gamma}$ 是多项式衰减的，其中 $\gamma \in (0.5, 1]$。

随机优化解决的问题之一是模型参数的不确定性量化问题，最大似然估计或最大后验估计都无法捕获参数不确定性，并且可能会过拟合。贝叶斯方法捕获参数不确定性的典型方法是马尔可夫链蒙特卡罗 (MCMC) 技术（Robert & Casella，2004 ^[9]）。而在本文中，我们将考虑一种被称为 Langevin 动力学的 MCMC 技术（Neal，2010 ^[7]）。和之前一样，它也采用梯度步骤，但会将高斯噪声注入参数更新，使其不会坍缩为最大后验解j

$$\begin{align*} \Delta \theta_t &= \frac{\epsilon}{2}(\nabla \log p(\theta_t) + \sum^N_{i=1} \nabla \log p(x_i | \theta_t) ) + \eta_t\\ \eta_t &\sim \mathcal{N}(0, \epsilon) \tag{3} \end{align*}$$

梯度的步长和所注入的噪声方差之间是平衡的，这样后验样本的方差才能与后验方差能够保持匹配。 Langevin 动力学的动机和出发点是作为随机微分方程（其平衡分布等于后验分布）的离散化形式。为了纠正离散化产生的误差，可以将 式（3） 视为提议分布，并使用 Metropolis-Hastings 进行校正。有趣的是，当我们减小 $\epsilon$ 值时，离散化误差随之减小了，因此拒绝率接近于零。不过，典型 MCMC 实践允许一个调整步长的初始适应阶段，然后固定步长以确保此后的平稳马尔可夫链。

更复杂的技术使用具有动量变量的哈密顿动力学来允许参数移动更大的距离，而不会出现 Langevin 动力学的低效随机游走行为（Neal，2010 ^[7]）。但据我们所知，迄今为止提出的所有 MCMC 方法，都需要在每次迭代时对整个数据集进行计算，从而造成大型数据集的计算代价非常高。

3. 随机梯度朗之万动力学

鉴于随机梯度算法 式（1） 和 Langevin 动力学 式（3） 之间的相似性，考虑结合两种方法的思想是很自然的。这允许高效使用大型数据集，同时允许以贝叶斯方式捕获参数不确定性。方法很简单：使用 Robbins-Monro 随机梯度，添加一定量的与所用步长平衡的高斯噪声，并允许步长变为零。提议的更新很简单：

$$\Delta \theta_t = \frac{\epsilon_t}{2} \left(\nabla \log p(\theta_t) + \frac{N}{n} \sum^n_{i=1} \nabla \log p(x_{ti}|\theta_t) \right) + \eta_t \\ \eta_t \sim \mathcal{N}(0, \epsilon_t) \tag{4}$$

其中步长以满足 式（2） 的速率向零减小。这允许对梯度中的随机性进行平均，以及渐近变为零的 MH 拒绝率，这样我们就可以简单地忽略 MH 接受步骤，这需要对整个数据集的概率进行评估。

在本节的其余部分，我们将给出一个直观的论证，说明为什么当 $t \rightarrow \infty$ 时，$\theta_t$ 会接近后验分布中的样本。特别是，我们将证明对于大 $t$，更新 式（4） 将接近 Langevin 动力学 式（3） ，它收敛于后验分布。让

$$g(θ) = \nabla \log p(θ) + \sum^N_{i=1} \nabla \log p(x_i|θ) \tag{5}$$

是 $θ$ 处对数概率的真实梯度，并且

$$h_t(θ) = \nabla \log p(θ) + \frac{N}{n} \sum^n_{i=1} \nabla \log p(x_{ti}|θ) − g(θ) \tag{6}$$

随机梯度则为 $g(θ)+h_t(θ)$，其中 $h_t(θ)$ 是一个零均值随机变量（由于在步骤 $t$ 中选择的数据项的随机性），方差有限 $V(θ)$，并且 式（4） 是，

$$\Delta \theta_t = \frac{\epsilon_t}{2} (g(\theta_t) + h_t(\theta_t)) + \eta_t,\\ \eta_t \sim \mathcal{N}(0, \epsilon_t) \tag{7}$$

式（7） 中的随机性有两个来源：方差为 $\epsilon_t$ 的注入高斯噪声，以及方差为 $(\frac{\epsilon_t}{2})^2V(\theta_t)$ 的随机梯度中的噪声。第一个观察结果是，对于大的 $t$，$\epsilon_t \rightarrow 0$，注入的噪声将主导随机梯度噪声，因此 式（7） 将有效地成为 Langevin 动力学 式（3） 。第二个观察是当 $\epsilon_t \rightarrow 0$ 时，Langevin 动力学的离散化误差可以忽略不计，因此 MH 拒绝概率将接近 $0$，我们可以简单地忽略这一步。

换句话说， 式（4） 、 式（7） 有效地定义了一个非平稳马尔可夫链，使得对于所有大的 $t$，第 $t$ 步转移算子将具有 $θ$ 上的后验分布作为其平衡分布。我们要解决的下一个问题是参数序列 $\theta_1, \theta_2,\ldots$，会收敛到后验分布。因为马尔可夫链不是固定的并且步长减小到 $0$，所以目前还不清楚是不是这种情况。为了证明这确实是真的，我们将证明子序列 $\theta_{t1} , \theta_{t2} ,\ldots$ 将按预期收敛到后验，因此整个序列也将收敛。

首先确定一个 $\epsilon_0$ 使得 $0 < \epsilon_0 \ll 1$。由于 $\{\epsilon_t\}$ 满足步长性质 式（2） ，我们可以找到一个子序列 $t_1 < t_2 < \ldots$ 使得 $\sum^{t_{s+1}}_{t=t_s+1} \epsilon_t \rightarrow \epsilon_0$ 作为 $s \rightarrow \infty$。由于每个步骤的注入噪声是独立的，对于足够大的 $s$，总注入噪声 $\| \sum^{t_{s+1}}_{t=t_s+1} \eta_t\|_2$，在步骤 $t_s$ 和 $t_{s+1}$ 之间将为 $\mathcal{O}(\sqrt{\epsilon_0})$。我们现在表明，由于这些步骤中梯度的随机性，总噪声将由总注入噪声决定。由于 $\epsilon_0 \ll 1$，对于 $t_s$ 和 $t_{s+1}$ 之间的 $t$，我们可以取 $\|\theta_t − \theta_{t_s} \|_2 \ll 1$。假设梯度 $g(·)$ 平滑变化（例如，它们在 第 5 节 的模型中是 Lipschitz 连续的），总随机梯度为：

$$\sum^{t_{s+1}}_{t_s+1} \frac{\epsilon_t}{2} (g(\theta_t) + h_t(\theta_t)) \\ = \frac{\epsilon_0}{2} g(\theta_{t_s} ) + \mathcal{O}(\epsilon_0) + \sum^{t_{s+1}}_{t_s+1} \frac{\epsilon_t}{2} h_t(\theta_t) \tag{8}$$

由于参数在 $t_s$ 和 $t_{s+1}$ 之间变化不大，$h_t(\theta_t)$ 中的随机性将由小批量选择的随机性决定。假设这些是随机和独立选择的，每个 $t$ 的 $h_t(\theta_t)$ 基本上是独立同分布的（如果通过对整个数据集进行随机分区选择小批量，$h_t(\theta_t)$ 将是负相关的，并且不会改变结果这里）。因此，\sum^{t_{s+1}_{t=t_s +1} \frac{\epsilon_t}{2} h_t(\theta_t) 的方差为 $\mathcal{O}(\sum_t \frac{\epsilon^2_t}{4})$ 并且

$$= \frac{\epsilon_0}{2} g(\theta_{t_s}) + \mathcal{O}(\epsilon_0) + \mathcal{O} \left( \sqrt{\sum^{t_{s+1}}_{t=t_s +1} \frac{\epsilon^2_t}{4}} \right) \\ = \frac{\epsilon_0}{2} g(\theta_{t_s} ) + \mathcal{O}(\epsilon_0)$$

最后一个方程表明，总随机梯度步长大约是 $\theta_{t_s}$ 处的精确梯度步长，步长为 $\epsilon_0$，偏差由 $\mathcal{O}(\epsilon_0)$ 决定。由于这又受总注入噪声 $\mathcal{O}(\sqrt{\epsilon_0})$ 的支配，这意味着序列 $\theta_{t1}, \theta_{t2},\ldots$ 将以固定步长 $\epsilon_0$ 逼近朗之万动力学生成的序列，因此它将收敛于后验分布。另请注意，它将具有无限的有效样本量

这个论点的含义是我们可以使用随机梯度 Langevin 动力学作为 “随时” 和通用的算法。在初始阶段，随机梯度噪声将占主导地位，算法将模仿有效的随机梯度上升算法。在后期，注入的噪声将占主导地位，因此该算法将模仿 Langevin 动力学 MH 算法，该算法将在两者之间平滑过渡。然而，缺点是为了保证算法的工作，重要的是将步长减小到零，这样算法的混合速率将随着迭代次数的增加而减慢。为了解决这个问题，我们可以在步长减小到 MH 拒绝率被认为可以忽略的临界水平以下时保持步长不变，或者使用此算法进行老化，但切换到不同的 MCMC 算法，以更有效地利用整个数据集稍后。这些替代方案可以表现得更好，但需要进一步手动调整，并且超出了本文的范围。本文的重点是演示一种实用算法，该算法仅使用小批量数据即可实现正确的贝叶斯学习。

4. 后验采样

本节将考虑使用随机梯度 Langevin 动力学算法从后验分布中生成样本。我们首先推导出算法何时从随机优化过渡到 Langevin 动力学的估计。这个想法是我们应该只在它进入后采样阶段后才开始收集样本，这要等到它成为 Langevin 动力学之后才会发生。然后我们讨论算法如何随数据集大小 $N$ 进行缩放，并粗略估计算法遍历整个后验所需的迭代次数。最后，我们讨论如何使用获得的样本来形成后验期望的蒙特卡罗估计。

4.1 过渡到 Langevin 动力学阶段

我们首先推广我们的方法以允许进行预处理，这可以通过更好地使步长适应后验的局部结构来显著加快速度（Roberts & Stramer, 2002 ^[10]; Girolami & Calderhead, 2011 ^[4]）。例如，某些维度可能具有更大的曲率，从而导致更大的梯度。在这种情况下，对称预处理矩阵 $M$ 可以将所有维度转换为相同的尺度。预条件随机梯度 Langevin 动力学很简单

$$\Delta \theta_t = \frac{\epsilon_t}{2} M \left( g(\theta_t) + h_t(\theta_t) \right) + \eta_t,\\ \eta_t \sim \mathcal{N}(0, \epsilon_t M )$$

如前所述，算法是处于随机优化阶段还是 Langevin 动力学阶段取决于注入噪声的方差，即 $\epsilon_t M$ 与随机梯度的方差。由于随机梯度是当前 mini-batch 的总和，如果它的大小 $n$ 足够大，中心极限定理就会生效，真实梯度 $g(\theta_t)$ 周围的变化 $h_t(\theta_t)$ 将变为正态分布。然后可以根据经验协方差估计其协方差矩阵：

$$V(\theta_t) ≡ V[h_t(\theta_t)] \approx \frac{N^2}{n^2} \sum^n_{i=1} (s_{ti} − \overline{s_t})(s_{ti} − \overline{s_t})^⊤ \tag{9}$$

其中 $s_{ti} = \nabla \log p(x_{ti}|\theta_t) + \frac{1}{N} \nabla \log p(\theta_t)$ 是数据项 $i$ 在迭代 $t$ 时的得分，$\overline{s_t} = \frac{1}{n} \sum^n_{i=1} s_{ti}$ 是经验平均值。请注意，$V(\theta_t) = \frac{N^2}{n} V_s$，其中 $V_s$ 是分数 $\{s_{ti}\}$ 的经验协方差，因此缩放为 $\frac{N^2}{n}$。从这里我们看到随机梯度步长的方差是 $\frac{\epsilon^2_t N^2}{4n} MV_sM$，所以为了让注入的噪声在所有方向上占主导地位，我们需要条件

$$\frac{\epsilon_t N^2}{4n} λ_{max} (M^{\frac{1}{2}} V_s M^{\frac{1}{2}}) = α \ll 1 \tag{10}$$

其中 $λ_{max}(A)$ 是 $A$ 的最大特征值。换句话说，如果我们选择一个步长，使得样本阈值 $α \ll 1$，算法将处于其 Langevin 动力学阶段，并将近似从后验采样。

我们现在可以通过 Fisher 信息将采样阈值处的步长与后验方差相关联，Fisher 信息与 $V_s$ 相关，因为 $I_F \approx N V_s$，并且与后验方差 $\Sigma_θ \approx I^{-1} F$ 相关。使用这些关系以及 式(10)，我们看到采样阈值处的步长为 $\epsilon_t \approx \frac{4αn}{N} λ_{min}(\Sigma_θ)$。由于 Langevin 动力学通过随机游走探索后验，使用这个步长意味着我们需要 $N/n$ 步的顺序来遍历后验，即我们处理整个数据集。所以我们看到这个方法不是灵丹妙药。然而，该方法的优势在于它的便利性：随机优化可以平滑地自动过渡到后验采样，而无需更改更新方程。即使不测量采样阈值，也可以享受防止过拟合和执行贝叶斯学习能力的好处。仅当需要用有限的样本集合忠实地表示后验分布时，测量采样阈值才重要。

4.2 估计后验期望

由于 $\theta_1, \theta_2,\ldots$ 收敛到后验分布，我们可以通过简单地取样本平均值 $\frac{1}{T} \sum^T_{t=1} f(\theta_t)$ 来估计某个函数 $f(θ)$ 的后验期望 $E[f (θ)]$（通常在 MCMC 中，我们可以移除初始老化阶段，比如使用采样阈值进行估计）。由于 $f(\theta_t)$ 是 $E[f (θ)]$ 的渐近无偏估计量，因此该样本平均值将是一致的。但观察到，由于步长减小，马尔可夫链的混合率也随之降低，简单样本平均会过分强调样本间相关性较高的序列尾端，导致方差较大估算器。相反，我们建议使用步长来对样本进行加权：

$$E[f (θ)] \approx \frac{ \sum^T_{t=1} \epsilon_t f(\theta_t)}{\sum^T_{t=1} \epsilon_t} \tag{11}$$

由于 $\sum^{\infty}_{t=1} \epsilon_t = \infty$，该估计量也将是一致的。直觉是马尔可夫链混合的速率与步长成正比，因此我们期望 $\{\theta_1,\ldots , \theta_T \}$ 与 $\sum^T_{t=1} \epsilon_t$ 成正比，并且每个 $\theta_t$ 将贡献一个与 $\epsilon_t$ 成正比的有效样本量。

5. 实验

5.1 简单示范

我们首先在一个仅涉及两个参数的简单示例中演示随机梯度 Langevin 算法的工作原理。为了使后验多模态更有趣一些，我们使用了具有绑定均值的高斯混合：

$$\theta_1 \sim \mathcal{N}(0, \sigma^2_1) ; \theta_2 \sim \mathcal{N}(0, \sigma^2_2)\\ x_i \sim \frac{1}{2} \mathcal{N}(\theta_1, \sigma^2_x) + \frac{1}{2} \mathcal{N}(\theta_1 + \theta_2, \sigma^2_x)$$

其中 $\sigma^2_1 = 10$，$4\sigma^2_2 = 14$ 和 $\sigma^2_x = 2$。$100$ 个数据点是从 $\theta_1 = 0$ 和 $\theta_2 = 1$ 的模型中提取的。在这个参数设置下有一个峰值，但在 $\theta_1 = 1$, $\theta_2 = −1$, 参数间负相关性强。我们运行批量大小为 $1$ 的随机梯度 Langevin 算法，并使用 $10000$ 次扫描整个数据集。步长为 $\epsilon_t = a(b + t)^{-\gamma}$，其中 $\gamma = .55$，$a$ 和 $b$ 的设置使得 $\epsilon_t$ 在运行期间从 $.01$ 减小到 $.0001$。我们从 图 1 中看到，估计的后验分布非常准确。在 图 2 中，我们看到随机梯度 Langevin 算法确实有两个阶段：第一个阶段随机梯度噪声主导注入噪声，第二个阶段则相反。为了探索作为步长函数的拒绝率的缩放比例，我们重新进行了实验，步长从 $10^{−2}$ 呈指数下降到 $10^{−8}$。在原始实验中，步长的动态范围不够宽，无法进行目视检查。 图 2（右） 显示随着步长减小，拒绝概率降低到零。

图 1: 真实和估计的后验分布。

图 2: 左：随机梯度噪声和注入噪声的方差。右图：拒绝概率与步长。我们报告每次扫描数据集时每次迭代的平均拒绝概率。

5.2 逻辑斯谛回归

我们将随机梯度 Langevin 算法应用于贝叶斯逻辑回归模型。给定相应输入向量 $x_i$ 的第 $i$ 个输出 $y_i \in \{−1, +1\}$ 的概率被建模为：

$$p(y_i|x_i) = σ(y_iβ^{T}x_i) \tag{12}$$

其中 $β$ 是参数，$σ(z) = \frac{1 }{1+\exp(−z)}$ 。通过将 $1$ 作为 $x_i$ 中的条目，偏置参数被吸收到 $β$ 中。我们对标度为 $1$ 的 $β$ 使用拉普拉斯先验。对数似然的梯度为：

$$\frac{\partial}{\partial β} \log p(y_i|x_i) = σ(−y_i β^{T} x_i) y_i x_i \tag{13}$$

而先验的梯度只是 $−\text{sign}(β)$，它是按元素应用的。

我们将我们的推理算法应用于由 (Lin 等, 2008 ^[6]) 从 UCI 成人数据集派生的 a9a 数据集。它包含 $32561$ 个观察值和 $123$ 个特征，我们使用的批量大小为 $10$。 图 3 显示了 $50$ 次运行的结果，模型在每次运行中随机 $80\%$ 的数据集上进行训练，并在另外 $20\%$ 的数据集上进行测试。我们看到联合概率和精度都迅速增加，联合概率在最多 $10$ 次迭代后收敛，而精度在通过数据集的不到 $1$ 次迭代后收敛，证明了随机梯度 Langevin 动力学的效率。

图 3: 每个数据项的平均对数联合概率（左）和测试集的准确度（右）作为整个数据集扫描次数的函数。红色虚线表示 10 次迭代后的准确性。结果取 50 次运行的平均值；蓝色虚线表示 1 个标准差。

5.3 独立成分分析

在下文中，我们将简要回顾一种流行的基于随机（自然）梯度优化的 ICA 算法（Amari 等人，1996 年 ^[1]）。我们从假设独立的重尾边际分布的概率模型开始

$$p(\mathbf{x}, W) = |\det(W)| \left [ \prod_i pi(\mathbf{w}_{i}^T \mathbf{x}) \right ] \prod_{ij} \mathcal{N}(W_{ij}; 0, λ)$$

我们在权重上使用了高斯先验。已经发现，如果我们使用自然梯度，梯度下降的效率可以显著提高。这是通过梯度与项 $W^TW$ 的后乘法实现的（Amari 等人，1996 ^[1]）。如果我们选择 $p_i(y_i) = \frac{1}{4 \cos h^2( \frac{1}{2} y_i)}$ 且 $y_i = \mathbf{w}^T_i \mathbf{x}$，我们得到

$$\mathcal{D}_W \doteq \nabla_W \log[p(X, W)] W^TW = \left( N \mathbf{I} − \sum^{N}_{n=1} \tanh( \frac{1}{2} \mathbf{y}_n) \mathbf{y}^T_n \right ) W − λ WW^TW \tag{15}$$

$W^TW$ 项就像一个预处理矩阵（参见第 4.1 节），$M_{ij,kl} = δ_{ik}(W^TW )_{jl}$ 在交换下是对称的 $(i \leftrightarrow k, j \leftrightarrow l)$。可以证明，$M$ 的逆由 $M^{−1} = \boldsymbol{\delta}(W^TW)^{−1}$ 给出，矩阵平方根为 $\sqrt{M} = \boldsymbol{\delta} \sqrt{W^TW}$，其中 $\sqrt{W^TW} = U \Gamma^{\frac{1}{2}} U^T$ 如果 $W^TW = U \Gamma U^T$。

因此，朗之万动力学的更新方程变为，

$$W_{t+1} = W_t + \frac{1}{2} \epsilon_t \mathcal{D}_W + \boldsymbol{\eta}_t \sqrt{W^TW} \tag{16}$$

其中 $t$ 的每个元素都服从方差 $\epsilon_t$ 的正态分布：$η_{ij,t} \sim \mathcal{N}[0, \epsilon_t]$。我们的随机版本简单地近似于对数据案例求和的梯度部分与对大小为 $n$ 的小批量求和，并将结果乘以 $N/n$ 以使其恢复到正确的比例。我们还根据 $\epsilon_t \approx a(b + t)^{-\gamma}$ 对步长进行退火

这个采样器提出了一个新的状态 $W^*$，如等式 16 所示。但是请注意，我们对所有数据案例求和并且我们不退火步长。其次，我们需要根据所有数据案例接受或拒绝建议的步骤，以保证细节平衡。建议分布由（抑制对 $t$ 的依赖）给出，

$$q(W \rightarrow W^*) = \mathcal{N} \left[ W^*; W + \frac{1}{2} \epsilon \mathcal{D}_W ； \epsilon M \right] \tag{17}$$

其中指数中的二次函数可以方便地计算为

$$\frac{-1}{2\epsilon} \text{tr} \left[(δW − \frac{1}{2} \epsilon \mathcal{D}_W)(W^TW )^{−1}(δW − \frac{1}{2} \epsilon \mathcal{D}_W)^T \right] \tag{18}$$

其中 $δW = W^* − W$ 和归一化常数需要量 $\det M = \det(W^TW)^D$。然后通过通常的 Metropolis Hastings 规则给出接受/拒绝步骤：

$$p(\text{accept}) = \min \left[ 1, \frac{p(W^*)q(W^* \rightarrow W)}{p(W)q(W \rightarrow W^*)} \right] \tag{19}$$

最后，为了计算 式（10） 的采样阈值，我们可以使用

$$M^{\frac{1}{2}} \Nu(s) M^{\frac{1}{2}} = \text{cov}_n \left[\left( \frac{1}{N} \nabla \log p(W ) + \nabla \log p(x_i|W) \right) (W^TW )^{\frac{1}{2}} \right] \tag{20}$$

使用 $\text{cov}_n$ 样本协方差在 $n$ 个数据案例的小批量上。

为了展示我们的 ICA 贝叶斯变体的实用性，我们为独立组件定义了以下“不稳定性”指标：

$$\mathcal{I}_i = \sum_j \operatorname{var}(W_{ij}) \operatorname{var}(x_j) \tag{21}$$

其中 $ \operatorname{var}(W_{ij})$ 是在后验样本上计算的，而 $ \operatorname{var}(x_j)$ 是在数据案例上计算的。我们用 $x_j$ 的方差缩放权重条目 $W_{ij}$ 的方差的原因是源的方差 $y_i = \sum_j W_{ij}x_j$ 对于所有 $i$ 都近似相等，因为它们适合分布 $p_i(y_i) = \frac{1}{4 \cosh^2 (\frac{1}{2}y_i)}$。

图 4: 左侧两个数字：随时间变化的随机朗之万动力学和修正的朗之万动力学的阿马里距离。粗线代表在线平均值。删除了前几百次迭代以显示波动的规模。右二图：第 5.3.1 节计算的 6 个独立分量的不稳定指数，用于随机 Langevin 动力学和校正的 Langevin 动力学。

图 5. 在 $W_{1,1}-W_{1,2}$ 和 $W_{1,1} - W_{2,1} 轴上测量的随机 Langevin 和校正 Langevin 动力学的人工数据集的后验密度估计。

5.3.1 人工数据

在第一个实验中，我们在六个通道中生成了 $1000$ 个数据案例 IID。三个通道具有高峰态分布，而其他三个通道呈正态分布。我们运行批量大小为 $100$ 的随机 Langevin 动力学，共进行 $500,000$ 次迭代和多项式退火计划 $\epsilon_t = \frac{4}{N} t^{−0.55}$。在大约 $10,000$ 次迭代后，达到了 $α = 0.1$ 的采样阈值。在这一点上，我们将“混合距离”记录为 $D_0 = \epsilon_t$，并且仅当最后一次采样时间的总和 $\sum_t \epsilon_t$ 超过 $D_0$ 时才收集样本（换句话说，随着 $\epsilon_t$ 的减小，我们每单位时间收集的样本更少）。我们注意到，简单地收集所有样本对最终结果没有明显影响。 $W$ 的最后估计用于初始化校正的 Langevin 动力学（这样做是为了迫使采样器进入相同的局部最大值），之后我们还收集了 $500,000$ 个样本。对于校正后的 Langevin，我们使用常数步长 $\epsilon = \frac{0.1}{N}$。

图 4 的左两个图分别显示了随时间变化的随机和校正朗之万动力学的 Amari 距离 (Amari 等, 1996 ^[1])。右边的两个图显示了我们提出的不稳定指数的排序值。 图 5 显示了 $W$ 的后验分布的二维边际密度估计。 ICA 无法确定高斯分量，这一事实通过查看后验分布得到验证。事实上，随机 Langevin 算法混合了多种峰值，这些峰值可能对应于高斯分量的不同线性组合。在较小程度上，经过修正的 Langevin 还探索了两种模式。由于后验分布的复杂结构，稳定性指数在高斯分量的两种采样算法之间变化很大（实际上在不同的运行中也会变化）。我们验证了最后三个分量对应于稳定的高峰度分量。

5.3.2 MEG数据

我们从 http://www.cis.hut.fi/projects/ica/eegmeg/MEG_data.html 下载了 MEG 数据集。 有 $122$ 个通道和 $17730$ 个时间点，我们从中提取了前 $10$ 个通道用于我们的实验。为了初始化采样算法，我们首先运行 fastICA (Hyvarinen, 1999 ^[5]) 以找到去混合矩阵 $W$ 的初始估计。然后我们运行随机 Langevin 并校正 Langevin 动力学以从后部采样。这些设置与之前的实验非常相似，随机 Langevin 的时间表为 $\epsilon_t = \frac{0.1}{N} t^{−0.55}$，校正的 Langevin 的恒定步长为 $1/N$。我们在 $800$ 秒内获得了 $500,000$ 个随机朗之万样本，在 $9000$ 秒内获得了 $100,000$ 个校正朗之万样本。我们在视觉上验证了随机 Langevin 和校正的 Langevin 动力学的二维边缘分布非常相似。不稳定性值如 图 6 所示。由于没有高斯分量，我们看到两种采样算法的稳定性指数非常相似。经验证，最稳定的分量对应于高峰度源（峰度 = 15.4），而最不稳定的分量更接近峰度为 $3.4$（$2$ 对应于高斯）的高斯噪声。这些发现证实随机 Langevin 程序产生准确的后验分布，与完善的 MCMC 程序完全一致。

图 6: 随机 Langevin（左）和校正 Langevin（右）分别为 MEG 数据集的 10 个成分的不稳定性指数。

6.讨论

如果我们衡量 “每单位计算获得的预测精度”，则随机梯度优化是最有效的算法之一（Bottou & Bousquet，2008 ^[3]）。由于二次采样噪声，参数估计值围绕其 MAP 值波动。常识是必须将这些步长退火到零才能达到稳定点。然而，我们认为不应优化超出后验分布的范围。后验代表精度的内在统计尺度，并且尝试以更高的精度确定参数值会以额外的计算代价运行过拟合的风险。

来自后验分布的 MCMC 采样当然可以解决过拟合问题。然而，一般的 MCMC 算法需要在每次迭代中查看所有数据，因此失去了随机逼近方法的优势。本文首次提供了一个出人意料的简单解决方案，代表了两全其美：坚持使用随机梯度，但仍然从后验中采样。

但也许随机梯度 Langevin 动力学的最大优势是随机优化无缝过渡到后验采样这一事实。通过简单地添加具有正确方差的高斯噪声，我们的方法会自动执行 “提前停止”，而无需担心。事实上，我们已经表明，通过多项式退火计划，获得的样本将忠实地渐进地表示后验分布。

我们相信这项工作只是暂时的第一步，可以进一步研究基于随机梯度的高效 MCMC 采样。有趣的研究方向包括更强的理论提供收敛性的可靠证明，基于小批量数据推导 MH 拒绝步骤，将算法扩展到动态系统的在线估计，以及基于更复杂的哈密顿蒙特卡洛方法推导算法不受随机游走行为的影响。

参考文献

• [1] Amari, S., Cichocki, A., and Yang, H.H. A new algorithm for blind signal separation. In Neural Information Processing Systems, volume 8, pp. 757-763, 1996.
• [2] Andrieu, C., de Freitas, N., and Doucet, A. Sequential MCMC for Bayesian model selection. In Proceedings of the IEEE Signal Processing Workshop on Higher-Order Statistics, pp. 130-134, 1999.
• [3] Bottou, L. and Bousquet, O. The tradeoffs of large scale learning. In Advances in Neural Information Processing Systems, volume 20, pp. 161-168, 2008.
• [4] Girolami, M. and Calderhead, B. Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society B, 73:1-37, 2011.
• [5] Hyvarinen, A. Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10(3):626-634, 1999.
• [6] Lin, C.-J., Weng, R. C., and Keerthi, S. S. Trust region Newton method for large-scale logistic regression. Journal of Machine Learning Research, 9:627-650, 2008.
• [7] Neal, R. M. MCMC using Hamiltonian dynamics. In Brooks, S., Gelman, A., Jones, G., and Meng, X.-L. (eds.), Handbook of Markov Chain Monte Carlo. Chapman & Hall / CRC Press, 2010.
• [8] Robbins, H. and Monro, S. A stochastic approximation method. Annals of Mathematical Statistics, 22(3):400- 407, 1951.
• [9] Robert, C. P. and Casella, G. Monte Carlo statistical methods. Springer Verlag, 2004.
• [10] Roberts, G. O. and Stramer, O. Langevin diffsions and metropolis-hastings algorithms. Methodology and Computing in Applied Probability, 4:337-357, 2002.