论文笔记 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$$