论文笔记 Stochastic Gradient Hamiltonian Monte Carlo （SGHMC）

随机梯度哈密顿量蒙特卡罗

0.摘要

Hamiltonian Monte Carlo (HMC) 采样方法提供了一种机制，用于在 Metropolis-Hastings 框架中定义具有高接受概率的远距离建议，从而能够比标准随机游走提案更有效地探索状态空间。近年来，此类方法的普及显着增长。然而，HMC 方法的一个限制是模拟哈密顿动力系统所需的梯度计算——这种计算在涉及大样本量或流数据的问题中是不可行的。相反，我们必须依赖从数据子集计算的噪声梯度估计。在本文中，我们探讨了这种随机梯度 HMC 方法的特性。令人惊讶的是，随机近似的自然实现可能非常糟糕。为了解决这个问题，我们**引入了一种变体，该变体使用二阶朗之万动力学，其摩擦项抵消了噪声梯度的影响，将所需的目标分布保持为不变分布。**模拟数据的结果验证了我们的理论。我们还将我们的方法应用于使用神经网络的分类任务和在线贝叶斯矩阵分解。

1.介绍

哈密​​顿蒙特卡罗 (HMC) (Duane et al., 1987; Neal, 2010) 采样方法提供了强大的马尔可夫链蒙特卡罗 (MCMC) 采样算法。这些方法根据我们期望样本的目标分布（势能） 和由一组“动量”辅助变量参数化的动能项来定义哈密顿函数。基于对动量变量的简单更新，可以模拟哈密顿动力学系统，该系统能够提出远距离状态的建议。目标分布在这些动态下是不变的；在实践中，需要对连续时间系统进行离散化，从而需要进行 Metropolis-Hastings (MH) 校正，但仍具有较高的接受概率。基于 HMC 在快速探索状态空间方面的吸引力，HMC 方法最近越来越受欢迎（Neal，2010；Hoffman & Gelman，2011；Wang 等人，2013）。

然而，HMC 的一个限制是必须计算势能函数的梯度以模拟哈密顿动力系统。我们越来越多地面临具有数百万到数十亿观察值的数据集，或者数据以流的形式进入的数据集，我们需要在线进行推断，例如在线广告或推荐系统。在这些越来越常见的海量批处理或流数据场景中，这种梯度计算是不可行的，因为它们利用了整个数据集，因此不适用于“大数据”问题。最近，在各种机器学习算法中，我们见证了利用基于小批量数据的梯度噪声估计来扩展算法的许多成功（Robbins & Monro, 1951; Hoffman et al., 2013; Welling &德，2011）。大多数这些发展都是基于优化的算法（Robbins & Monro, 1951; Nemirovski et al., 2009），一个问题是是否可以通过基于采样的算法来获得类似的效率，这些算法保持了贝叶斯的许多理想的理论特性推理。在采样环境中应用此类方法的一种尝试是最近提出的随机梯度朗之万动力学 (SGLD) (Welling & Teh, 2011; Ahn et al., 2012; Patterson & Teh, 2013)。该方法建立在不包括 HMC 的关键动量项的一阶朗之万动力学之上。

在本文中，我们探讨了将 HMC 的状态空间探索效率与随机梯度的大数据计算效率相结合的可能性。这样的算法将启用大规模的在线贝叶斯采样算法，并有可能快速探索后验。作为第一步，我们考虑简单地对 HMC 应用随机梯度修改并评估噪声梯度的影响。我们证明了随机梯度注入系统的噪声不再导致哈密顿动力学，期望的目标分布作为平稳分布。因此，即使在离散化动力系统之前，我们也需要纠正这种影响。可以通过 MH 步骤校正注入的梯度噪声，尽管这本身需要对整个数据集进行昂贵的计算。在实践中，人们可能会建议在 MH 校正之前进行长时间的模拟运行，但这会导致低接受率，因为哈密顿量与注入噪声的偏差很大。使用 (Korattikara et al., 2014; Bardenet et al., 2014) 的最新结果可能会提高该 MH 步骤的效率。在本文中，我们改为介绍一种随机梯度 HMC 方法，其中在动量更新中添加了摩擦。我们假设注入的噪声是高斯的，根据中心极限定理，并分析相应的动力学。我们表明，使用这种二阶朗之万动力学使我们能够将所需的目标分布保持为平稳分布。也就是说，摩擦抵消了注入噪声的影响。对于离散系统，我们考虑让步长趋向于零，这样就不需要 MH 步长，这给我们带来了显着的计算优势。根据经验，我们证明即使设置为较小的固定值，我们也有良好的性能。补充材料中提供了这种小方法的理论计算与精度权衡。

大量模拟实验验证了我们的理论结果并证明了

• (i) 精确 HMC
• (ii) 随机梯度 HMC 的简单实现（简单地用随机梯度替换梯度）
• (iii) 我们提出的方法之间的差异 加入摩擦。

我们还比较了 SGLD 的一阶朗之万动力学。 最后，我们将我们提出的方法应用于使用贝叶斯神经网络的分类任务和标准电影数据集的在线贝叶斯矩阵分解。 我们的实验结果证明了所提出算法的有效性。

2.HMC

假设我们想从给定一组独立观察 $x\in D$ 的后验分布中进行采样：
$$p(\theta |D)\propto \exp (-U(\theta ))\text{\hspace{1em}(1)}$$
其中势能函数 $U$ 由下式给出
$$U=-\sum _{x\in D}\log p(x|\theta )-\log p(\theta )\text{\hspace{1em}(2)}$$
Hamiltonian (Hybrid) Monte Carlo (HMC) (Duane et al., 1987; Neal, 2010) 提供了一种在 Metropolis-Hastings (MH) 框架中提出样本 $\theta$ 的方法，与标准随机游走建议相比， 这些建议是从基于引入一组辅助动量变量 $r$ 的哈密顿系统生成的。 也就是说，为了从 $p(\theta |D)$ 中采样，HMC 考虑从由 $(\theta ,r)$ 定义的联合分布中生成样本
$$\pi (\theta ,r)\propto \exp (-U(\theta )-\frac{1}{2}{r}^T{M}^{-1}r)\text{\hspace{1em}(3)}$$
如果我们简单地丢弃得到的 $r$ 样本，则样本具有边际分布 $p(\theta |D)$。 这里，$M$ 是质量矩阵，与 $r$ 一起定义了动能项。 $M$ 通常设置为单位矩阵 $I$，但是当我们有更多关于目标分布的信息时，可以用来对采样器进行预处理。 哈密顿函数定义为 $H(\theta ,r)=U(\theta )+\frac{1}{2}{r}^T{M}^{-1}r$。 直观地说，$H$ 用位置变量 $\theta$ 和动量变量 $r$ 测量物理系统的总能量。

为了提出样本，HMC 模拟哈密顿动力学
$$\{\begin{array}{l}d\theta =M^{-1}rdt\\ dr=-\nabla U(\theta )dt\end{array}\text{\hspace{1em}(4)}$$
为使等式 (4) 更加具体，2D 中的常见类比如下 (Neal, 2010)。 想象一个冰球在不同高度的无摩擦冰面上滑动。 势能项基于当前圆盘位置处的表面高度，而动能基于圆盘的动量 $r$ 及其质量 $M$。如果表面是平坦的 (任意 $\theta$ 都有$\nabla U(\theta )=0$ )，圆盘匀速运动。 对于正斜率$(\nabla U(\theta )>0)$，动能随着势能的增加而减小，直到动能为 $0(r=0)$。 然后冰球滑下山坡，增加其动能并减少势能。 回想一下，在 HMC 中，位置变量是直接的兴趣变量，而动量变量是人工构造（辅助变量）。

在任何区间 $s$ 上，方程4的哈密顿动力学, 定义了从时间 t 的状态到时间 t + s 的状态的映射。 重要的是，这种映射是可逆的，这对于表明动力学保持不变很重要。 同样，动力学保留了总能量 $H$，因此建议总是被接受。 然而，在实践中，我们通常不能从方程(4)的连续系统精确模拟, 而是考虑一个离散系统。 一种常见的方法是“跨越式”方法，该方法在 Alg.1 中进行了概述。 由于离散化引入的不准确性，必须执行 MH 步骤（即接受率不再是 1）。 然而，即使对于离最后状态很远的提案，接受率仍然很高。

HMC 最近取得了许多进展，使算法更加灵活，适用于各种环境。 “No U-Turn”采样器 (Hoffman & Gelman, 2011) 和 Wang 等人提出的方法。 （2013）允许自动调整步长，和模拟步数，m。 黎曼流形 HMC (Girolami & Calderhead, 2011) 利用黎曼几何来适应质量 M，使算法能够利用曲率信息执行更有效的采样。 我们尝试在专注于计算复杂性的正交方向上改进 HMC，但这些自适应 HMC 技术可能会与我们提出的方法相结合，以获得进一步的好处。

3. 随机梯度 HMC

在本节中，我们研究使用随机梯度实现 HMC 的含义，并提出对随机梯度估计引入的噪声更鲁棒的哈密顿动力学变体。 在所有场景中，不是直接使用方程式 (2)计算代价高昂的梯度 $\nabla U(\theta )$。这需要检查整个数据集 D，我们考虑基于从 $D$ 中随机均匀采样的小批量 $\tilde{D}$ 的噪声估计：
$$\nabla \tilde{U}(\theta )=-\frac{|D|}{|\tilde{D}|}\sum _{x\in \tilde{D}}\nabla \log p(x|\theta )-\nabla \log p(\theta ),\tilde{D}\subset D\text{\hspace{1em}(5)}$$
我们假设我们的观察 $x$ 是独立的，并且根据中心极限定理，将这个噪声梯度近似为
$$\nabla \tilde{U}(\theta )\approx \nabla U(\theta )+N(0,V(\theta ))\text{\hspace{1em}(6)}$$
这里，$V$ 是随机梯度噪声的协方差，它可以取决于当前模型参数和样本大小。 请注意，我们在方程式(6)中使用了符号的滥用, 其中 加入$N(\mu ,\Sigma )$ 表示引入了一个随机变量，该变量服从该多元高斯分布。 随着$\tilde{D}$的大小增加，这种高斯近似变得更加准确。 显然，我们希望 $minibatch$ 很小，以获得我们所追求的计算增益。 根据经验，在广泛的设置中，只需考虑数百个数据点的小批量大小就足以使中心极限定理近似准确（Ahn 等人，2012 年）。 在我们感兴趣的应用中，这种大小的小批量仍然代表着梯度计算成本的显着降低。

3.1 朴素随机梯度 HMC

随机梯度 HMC 最直接的方法是由 $\nabla \tilde{U}(\theta )$ 简单地替换算法中的 $\nabla U(\theta )$ 如算法1所示 。 参考方程(6)，这会在动量更新中引入噪声，变为 $\nabla r=-ϵ\nabla \tilde{U}(\theta )=-ϵ\nabla U(\theta )+N(0,{ϵ}^{2}V)$。 得到的离散时间系统可以看作是以下连续随机微分方程的离散化：
$$\{\begin{array}{l}d\theta =M^{-1}rdt\\ dr=-\nabla U(\theta )dt+N(0,2B(\theta )dt)\end{array}\text{\hspace{1em}(4)}$$
这里，$B(\theta )=\frac{1}{2}ϵV(\theta )$ 是梯度噪声贡献的扩散矩阵。 与最初的 HMC 公式一样，返回到连续时间系统以导出该方法的属性是有用的。 为了获得有关此设置的一些直觉，请考虑与第二节相同的曲棍球类比。 在这里，我们可以想象冰球在同一个冰面上，但也有一些随机的风。 这种风可能会将冰球吹得比预期的更远。 形式上，如定理 3.1 的推论 3.1 所给出，当 B 不为零时，等式（3）的 $\pi (\theta ,r)$在方程式(7)描述的动力学下不再是不变的。

定理 3.1:
由方程式（7）的动力学可知，令 $p_t(\theta ,r)$ 为 $(\theta ,r)$ 在时间 $t$ 的分布，定义 $p_t$ 的熵为 $h(p_t)=-{\int }_{\theta ,r}f(p_t(\theta ,r))d\theta dr$，其中 $f(x)=x\ln x$，假定 $p_t$是密度和梯度在无穷远处消失的分布。 此外，假设梯度消失的速度比 $\frac{1}{\ln p_t}$ 快。 然后， $p_t$ 的熵随着时间的推移而增加，增加的速度为
$${\partial }_{t}h(p_t(\theta ,r))={\int }_{\theta ,r}{f}^{\prime \prime }(p_t)({\nabla }_r{p}_t(\theta ,r){)}^{T}B(\theta ){\nabla }_r{p}_t(\theta ,r)d\theta dr\text{\hspace{1em}(8)}$$
方程 (8) 意味着${\partial }_{t}h(p_t(\theta ,r))\ge 0$ 因为 $B(\theta )$ 是一个半正定矩阵。

直观地说，定理 3.1 是正确的，因为无噪声哈密顿动力学保留了熵，而如果我们假设

• (i) $B(\theta )$ 是正定的，则附加噪声项严格增加熵（由于 Fisher 的正常满秩特性，这是一个合理的假设信息）
• (ii) 对于所有 $t$，${\nabla }_r{p}_t(\theta ,r)\ne 0$。 然后，共同地，熵随着时间严格增加。 这暗示了分布 $p_t$ 趋向于均匀分布的事实，这可能与目标分布 $\pi$ 相去甚远。

定理 3.1 和推论 3.1 的证明在补充材料中。因为在方程（7）的动态下不再是不变的，因此我们必须在考虑动力系统离散化引入的误差之前引入一个校正步骤。对于 MH 步骤的正确性（基于整个数据集），我们诉诸于 Neal (2010) 的 HMC 数据分割技术的相同论点。这种方法同样考虑小批量数据并按顺序模拟每批数据的（连续）哈密顿动力学。重要的是，Neal (2010) 暗示了这样一个事实，即从拆分数据情景中得到的 $H$ 可能与模拟后的全数据情景相去甚远，这导致较低的接受率，从而减少了模拟中明显的计算增益。从经验上看，正如我们在图 2 中所展示的，我们看到即使是来自噪声系统的有限长度模拟也可能与无噪声系统的模拟有很大差异。尽管此处考虑的基于小批量的 HMC 技术与 Neal (2010) 的技术略有不同，但我们在定理 3.1 中开发的理论围绕着所得的等式 (7) 3.1 的不变分布的高熵特性。 为在我们的实验和 Neal (2010) 的实验中观察到的 H 偏差提供了一些直觉。

基于使用噪声梯度的模拟的 $H$ 轨迹的不良行为导致复杂的计算与效率的权衡。一方面，在大型数据集中，在短暂的模拟运行后插入 MH 步骤的计算量非常大（其中 H 的偏差不太明显，接受率应该是合理的）。这些 MH 步骤中的每一个都需要使用所有数据进行昂贵的计算，从而抵消了考虑噪声梯度的计算收益。另一方面，MH 步骤之间的长时间模拟运行会导致非常低的接受率。每个拒绝都对应于使用所提出的 算法1 变体进行的浪费（嘈杂）梯度计算和模拟。 未来研究的一个可能方向是考虑使用 Korattikara 等人的最近结果(2014) 和 Bardenet 等人（2014）表明可以使用数据子集进行 MH。但是，我们改为在 Sec3,2 中考虑，对哈密顿动力学的直接修改，缓解了随机梯度引入的噪声问题。特别是，我们的修改使我们能够再次实现连续哈密顿动力系统的不变分布。

3.2. 具有摩擦力的随机梯度 HMC

在 3.1节，我们展示了具有随机梯度的 HMC 需要频繁且昂贵的 MH 校正步骤，否则，长时间的模拟运行具有低的接受概率。 相反，理想情况下，我们希望尽量减少注入噪声对动力学本身的影响，以缓解这些问题。 为此，我们考虑对等式（7）进行修改。 在动量更新中添加“摩擦”项：
$$\{\begin{array}{l}d\theta =M^{-1}rdt\\ dr=-\nabla U(\theta )dt-B{M}^{-1}rdt+N(0,2Bdt)\end{array}\text{\hspace{1em}(9)}$$
在这里和本文的其余部分，为了符号的简单性，我们省略了$B$对 $\theta$ 的依赖。 让我们再做一个曲棍球类比。 想象一下，我们现在打的是街头曲棍球而不是冰球，曲棍球会引入沥青摩擦。 仍然有随机的风在吹，但是表面的摩擦力阻止了冰球跑得很远。 也就是说，摩擦项 $B{M}^{-1}r$ 有助于降低能量 $H(\theta ,r)$，从而降低噪声的影响。 这种类型的动力系统通常被称为物理学中的二阶朗之万动力学（Wang & Uhlenbeck，1945）。 重要的是，我们注意到 SGLD (Welling & Teh, 2011) 中使用的 Langevin 动力学是一阶的，当摩擦项很大时，可以将其视为二阶动力学的极限情况。 有关此比较的更多详细信息，请参见本节末尾。

定理3.2
$\pi (\theta ,r)\propto \exp (-H(\theta ,r))$是方程描述的动力学的唯一平稳分布。
**证明：**令$G=\left[\begin{array}{cc}0& -I\\ I& 0\end{array}\right]$,$D=\left[\begin{array}{cc}0& 0\\ 0& B\end{array}\right]$,其中 $G$ 是反对称矩阵，$D$ 是对称（扩散）矩阵。 方程(9) 可以写成以下分解形式 (Yin & Ao, 2006; Shi et al., 2012)
$$d\left[\begin{array}{c}\theta \\ r\end{array}\right]=-\left[\begin{array}{cc}0& -I\\ I& B\end{array}\right]\left[\begin{array}{c}-\nabla U(\theta )\\ M^{-1}r\end{array}\right]dt+N(0,2Ddt)=-[D+G]\nabla H(\theta ,r)dt+N(0,2Ddt)$$
该动力系统下的分布演化由 Fokker-Planck 方程控制
∂ t p t ( θ , r ) = ∇ T { [ D + G ] [ p t ( θ , r ) ] ∇ H ( θ , r ) + ∇ p t ( θ , r ) } (10) \partial _tp_t(\theta , r)=\nabla ^T\{[D+G][p_t(\theta , r)]\nabla H(\theta , r)+\nabla p_t(\theta , r)\}\tag{10} ∂t​pt​(θ,r)=∇T{[D+G][pt​(θ,r)]∇H(θ,r)+∇pt​(θ,r)}(10)
有关详细信息，请参阅补充材料。 我们可以验证$\pi (\theta ,r)$ 在等式10下是不变的。 通过计算$[e^{H((\theta ,r)}\nabla H((\theta ,r))+\nabla e^{-H(\theta ,r)}]=0$。此外，由于扩散噪声的存在，(10)是方程的唯一平稳分布 。

总之，我们已经证明了方程式（9）给出的动力学。 具有与等式（4）的原始哈密顿动力学相似的不变性。 即使存在噪音。 关键是使用二阶朗之万动力学引入摩擦项。 我们修改后的动量更新也可以看作类似于部分动量更新（Horowitz，1991；Neal，1993），这也对应于二阶朗之万动力学。 在无噪声梯度的情况下，这种部分动量更新被证明不会大大改善 HMC (Neal, 2010)。 然而，正如我们已经证明的那样，这个想法在我们的随机梯度场景中至关重要，以抵消噪声梯度的影响。 我们将得到的方法称为随机梯度 HMC (SGHMC)。

连接到一阶朗之万动力学

正如我们之前讨论的，方程式（9）中引入的动力学， 与 SGLD 中使用的一阶朗之万动力学有关 (Welling & Teh, 2011)。 特别是，SGLD 的动力学可以看作是具有大摩擦项的二阶朗之万动力学。 为了直观地证明这种联系，在等式(9)中让 $B{M}^{-1}=\frac{1}{dt}$ 。 因为摩擦和动量噪声项非常大，动量变量 $r$ 的变化比 $\theta$ 快得多。 因此，相对于快速变化的动量， $\theta$ 可以认为是固定的。 我们可以简单地研究这个案例：
$$dr=-\nabla U(\theta )dt-B{M}^{-1}rdt+N(0,2Bdt)\text{\hspace{1em}(11)}$$

$r$ 的快速演化导致快速收敛到方程(11)的平稳分布。 由 $N(M{B}^{-1}\nabla U(\theta ),M)$ 给出。 现在让我们使用$r$~$N(M{B}^{-1}\nabla U(\theta ),M)$ 考虑 $\theta$ 的变化。又因为 $B{M}^{-1}=\frac{1}{dt}$ ，所以有
$$d\theta =-M^{-1}\nabla U(\theta )d{t}^2+N(0,2M^{-1}d{t}^2)\text{\hspace{1em}(12)}$$
这与 SGLD 的动态完全一致，其中 $M^{-1}$ 作为预处理矩阵（Welling & Teh，2011）。 直观地说，这意味着当摩擦很大时，动力学不依赖于由 $dr$ 表示的过去梯度的衰减序列，减少为一阶朗之万动力学。

到目前为止，在我们所考虑的所有事情中，我们都假设我们知道噪声模型 $B$。显然，在实践中情况并非如此。 想象一下，我们只是有一个估计 $\tilde{B}$ 。 正如将变得清楚的那样，改为引入用户指定的摩擦项 $C\ge \tilde{B}$ 并考
虑以下动态是有益
$$\{\begin{array}{l}d\theta =M^{-1}rdt\\ dr=-\nabla U(\theta )dt-C{M}^{-1}rdt+N(0,(C-\tilde{B})dt)+N(0,2Bdt)\end{array}\text{\hspace{1em}(13)}$$
生成的 SGHMC 算法显示在算法2 中。 请注意，该算法纯粹是根据用户指定的或可计算的量。 为了理解我们对动力学的选择，我们从完美估计 $B$ 的理想场景开始。

提案 3.1：
如果 $\tilde{B}=B$，则方程(13)的动力学, 产生平稳分布 $\pi (\theta ,r)\propto e^{-H(\theta ,r)}$。
证明: 动量更新简化为 $dr=-\nabla U(\theta )dt-C{M}^{-1}rdt+N(0,2Cdt)$，其中摩擦项 $C{M}^{-1}$ 和噪声项 $N(0,2Cdt)$。 注意到定理 3.2 的证明仅依赖于噪声和摩擦的匹配，因此在定理 3.2 中使用 $C$ 代替 $B$ 直接得出结果。

现在考虑在 $B$ 的不准确估计的更现实场景中引入 $C$ 项和修正动力学的好处。例如，最简单的选择是 $\tilde{B}=0$。虽然真正的随机梯度噪声 $B$ 显然非零，因为 步长$ϵ\to 0$，$B=\frac{1}{2}ϵV$ 变为 0，$C$ 占主导地位。 也就是说，动力学再次受可控注入噪声 $N(0;2Cdt)$ 和摩擦力 $C{M}^{-1}$ 的控制。 也可以设置 $\tilde{B}=\frac{1}{2}ϵ\tilde{V}$ ，其中 $\tilde{V}$ 是使用经验 Fisher 信息估计的，如 (Ahn et al., 2012) 中的 SGLD。

计算复杂度:略

用动量连接到 SGD：
在随机梯度下降 (SGD) 中添加动量项是常见的做法。 在概念上，带动量的SGD 和 SGHMC 之间存在明确的关系，在这里我们将这种联系形式化。 令 $v=M^{-1}r$，我们首先重写 算法 中的更新规则为：
$$\{\begin{array}{l}\Delta \theta =v\\ \Delta v=-{ϵ}^2{M}^{-1}\tilde{U}(\theta )-ϵM^{-1}Cv+N(0,2{ϵ}^3{M}^{-1}(C-\tilde{B})M^{-1})\end{array}\text{\hspace{1em}(14)}$$
定义 $\eta ={ϵ}^2{M}^{-1},\alpha =ϵM^{-1}C,\hat{\beta }=ϵM^{-1}\hat{B}$，那么更新规则重写为：
$$\{\begin{array}{l}\Delta \theta =v\\ \Delta v=-\eta \nabla \tilde{U}(\theta )-\alpha v+N(0,2(\alpha -\tilde{B})\eta )\end{array}\text{\hspace{1em}(15)}$$
与具有动量法的 SGD 相比，从方程式(15)中可以清楚地看出,$\eta$ 对应于学习率和 $1-\alpha$ 动量项。 当噪声被去除（通过 $C=\hat{B}=0$）时，SGHMC 自然地减少为具有动量的随机梯度方法。 我们可以使用等式 (15)的更新规则，运行SGHMC，借鉴SGD参数设置的经验，有动力指导我们选择SGHMC设置。 例如，我们可以设置一个固定的小数（例如，0:01 或 0:1），选择学习率 ，然后固定 $\hat{\beta }=\eta \hat{V}/2$。 更复杂的策略涉及使用动量调度（Sutskever et al., 2013）。 我们详细说明了如何在补充材料中选择这些参数。