基于随机梯度哈密顿蒙特卡洛(SGHMC)的非凸随机优化问题
字数 4619 2025-12-05 13:19:48

基于随机梯度哈密顿蒙特卡洛(SGHMC)的非凸随机优化问题

问题描述

考虑如下随机优化问题：

\[\min_{x \in \mathbb{R}^n} F(x) = \mathbb{E}_{\xi \sim \mathcal{D}} [f(x; \xi)] \]

其中，目标函数 \(F(x)\) 是非凸的，并且我们无法直接计算其精确梯度 \(\nabla F(x)\)，只能通过从分布 \(\mathcal{D}\) 中采样小批量（mini-batch）数据 \(\{\xi_i\}\) 来获得随机梯度估计 \(g(x) = \frac{1}{m} \sum_{i=1}^m \nabla f(x; \xi_i)\)。这里的随机性可能来源于数据噪声、随机子采样等。我们的目标是找到一个近似稳定点（例如，梯度范数较小的点），而非全局最优解。由于目标函数非凸且存在随机噪声，传统的随机梯度下降（SGD）可能收敛缓慢或陷入较差的局部极小点。本题将介绍随机梯度哈密顿蒙特卡洛 方法，它通过引入动量 和可控的随机噪声，帮助优化过程更好地探索参数空间，并可能逃离较差的局部极小点。

解题过程

我们将循序渐进地讲解SGHMC的核心思想和算法步骤，重点是其如何结合物理学中的哈密顿动力学 和随机优化 的需求。

步骤1：理解问题与SGD的局限性

  1. 问题本质：我们需要最小化一个期望风险 \(F(x)\)，但只能通过随机梯度 \(g(x)\) 来获取信息。\(g(x)\) 是真实梯度 \(\nabla F(x)\) 的无偏估计，但带有方差 \(\mathbb{E}[||g(x) - \nabla F(x)||^2] = \sigma^2 > 0\)。
  2. SGD更新：\(x_{t+1} = x_t - \eta g(x_t)\)，其中 \(\eta\) 是学习率。
  3. 局限性：
    • 方差导致震荡：梯度估计的噪声会导致优化路径在最小值附近剧烈震荡，需要小心翼翼地降低学习率。
    • 难以逃离局部极小点：在非凸问题上，SGD容易陷入“尖锐”的局部极小点或鞍点，缺乏有效的探索机制。

步骤2：从物理学中汲取灵感——哈密顿动力学

  1. 核心类比：将参数 \(x\) 视为一个粒子 的位置。在物理学中，要描述一个粒子的运动，不仅需要位置 \(x\)，还需要动量 \(v\)（通常 \(v = \dot{x}\)，即速度）。
  2. 哈密顿量：一个封闭物理系统的总能量 \(H(x, v)\) 通常是守恒的。一个常见的简单形式是动能加势能：\(H(x, v) = U(x) + \frac{1}{2} v^T M^{-1} v\)。
    • 在本优化问题中，我们令势能 \(U(x) = F(x)\)，即我们的目标函数。
    • \(M\) 是质量矩阵（常设为单位矩阵 \(I\)），因此哈密顿量简化为 \(H(x, v) = F(x) + \frac{1}{2} v^T v\)。
  3. 哈密顿方程：描述系统随时间的演化。

\[ \begin{aligned} dx &= v \, dt \\ dv &= -\nabla F(x) \, dt \end{aligned} \]

这组方程定义了一个**确定性**的轨迹。如果沿着这个轨迹运动，系统的总能量 $H(x, v)$ 是守恒的。这意味着粒子可以在势能面 $F(x)$ 上高效地“滑动”，而不会像SGD那样因摩擦力而停止在最近的谷底。

步骤3：引入摩擦与噪声——朗之万动力学

纯哈密顿动力学是保守的，会永远运动下去。为了使其最终稳定在低势能区域（即函数最小值附近），我们需要引入：

  1. 摩擦项：\(-\gamma v \, dt\)，其中 \(\gamma > 0\) 是摩擦系数。这会使动量 \(v\) 衰减，模拟能量耗散。
  2. 随机噪声项：\(\sqrt{2\gamma} \, dW\)，其中 \(W\) 是标准的维纳过程（布朗运动），\(\sqrt{2\gamma}\) 是噪声强度。这是涨落耗散定理的要求，确保系统最终能收敛到正确的平稳分布（玻尔兹曼分布）。

于是，我们得到朗之万动力学 方程：

\[\begin{aligned} dx &= v \, dt \\ dv &= -\nabla F(x) \, dt - \gamma v \, dt + \sqrt{2\gamma} \, dW \end{aligned} \]

可以证明，这个随机微分方程（SDE）的平稳分布是 \(p(x, v) \propto \exp(-H(x, v)) = \exp(-F(x) - \frac{1}{2} v^T v)\)。这意味着，在长期运行下，系统在 \((x, v)\) 空间中的分布正比于 \(e^{-F(x)}\)。因此，采样得到的 \(x\) 更有可能出现在 \(F(x)\) 值较小的区域。

步骤4：处理随机梯度——随机梯度朗之万动力学（SGLD）的启发

在实际优化中，我们无法得到精确的 \(\nabla F(x)\)，只有带噪声的估计 \(g(x)\)。如果我们简单地把 \(g(x)\) 代入朗之万方程：

\[dv = -g(x) \, dt - \gamma v \, dt + \sqrt{2\gamma} \, dW \]

那么此时更新中的随机性来源有两个：

  • 算法固有的、我们有意引入的噪声 \(\sqrt{2\gamma} dW\)。
  • 梯度估计本身引入的噪声 \(g(x) - \nabla F(x)\)。

如果这两部分噪声不匹配，系统就不会收敛到我们期望的平稳分布 \(e^{-H(x, v)}\)。梯度噪声的协方差记为 \(B(x)\)。为了使系统仍能收敛到正确的分布，我们需要在更新方程中引入一个额外的“修正”项。

步骤5：随机梯度哈密顿蒙特卡洛（SGHMC）算法

SGHMC巧妙地将梯度估计的噪声 \(B(x)\) 融入到动力学方程中，并重新定义噪声项，使得最终的离散时间算法形式简单。

  1. 修改的动力学方程：将朗之万方程改写为：

\[ \begin{aligned} dx &= v \, dt \\ dv &= -\color{blue}{g(x)} \, dt - \gamma v \, dt + \sqrt{2\gamma} \, dW \\ &= -\color{blue}{\nabla F(x)} \, dt - \gamma v \, dt + \sqrt{2\gamma} \, dW + \color{red}{(\nabla F(x) - g(x)) \, dt} \end{aligned} \]

注意，红色部分 $\nabla F(x) - g(x)$ 是均值为0、协方差为 $B(x)$ 的随机噪声。
  1. 关键简化：我们不去显式估计复杂的 \(B(x)\)。相反，我们注意到，如果我们能精确知道摩擦系数 \(\gamma\) 和梯度噪声的协方差 \(B\) 满足 \(B = 2\gamma I\)，那么红色部分的噪声和固有的 \(dW\) 噪声就可以合并。实际上，我们不知道 \(B\)。SGHMC的策略是设置一个足够大的摩擦系数 \(\hat{\gamma}\)，并定义 \(\beta = \hat{\gamma} I\)。然后，我们将动力学方程重写为：

\[ \begin{aligned} dx &= v \, dt \\ M dv &= -\nabla F(x) \, dt - \beta v \, dt + \sqrt{2(\beta - \hat{B})} \, dW + \color{green}{\sqrt{2\hat{B}} \, dW} \end{aligned} \]

其中，$\hat{B}$ 是**我们对梯度噪声协方差 $B$ 的一个估计**（实践中常设为 $\hat{B} = \frac{1}{2} \eta V$，$V$ 是梯度方差估计，$\eta$ 是学习率），绿色部分是**新增的可控噪声项**。
  1. 最终算法（离散化形式）：对上述SDE应用欧拉-丸山离散化，并定义学习率 \(\epsilon = \eta\)，质量矩阵 \(M = I\)，得到SGHMC的迭代更新公式：

\[ \begin{aligned} v_{t+1} &= (1 - \gamma \eta) v_t - \eta g(x_t) + \sqrt{2\gamma \eta} \, \mathcal{N}(0, I) \\ x_{t+1} &= x_t + v_{t+1} \end{aligned} \]

其中：
*   $v_t$ 是动量（速度）。
*   $\gamma$ 是摩擦系数（控制动量衰减）。
*   $\eta$ 是学习率（离散化的步长）。
*   $g(x_t)$ 是在 $x_t$ 处计算的随机梯度。
*   $\mathcal{N}(0, I)$ 是标准正态分布噪声。

注意，在标准实现中，为了简化，常常会忽略对 $\hat{B}$ 的显式估计，而将 $\gamma$ 当作一个可以调节的超参数。**算法中人为注入的噪声 $\sqrt{2\gamma \eta} \, \mathcal{N}(0, I)$ 正是SGHMC帮助逃离局部极小点、进行探索的关键**。

步骤6：算法步骤总结

  1. 初始化：选择初始参数 \(x_0\)，初始动量 \(v_0 = 0\)。设定超参数：学习率 \(\eta\)，摩擦系数 \(\gamma\)，总迭代步数 \(T\)，批量大小 \(m\)。
  2. 迭代循环（对于 \(t = 0, 1, ..., T-1\)）：
    a. 采样：从数据分布 \(\mathcal{D}\) 中随机抽取一个小批量样本 \(\{\xi_1, ..., \xi_m\}\)。
    b. 计算随机梯度：\(g_t = \frac{1}{m} \sum_{i=1}^m \nabla f(x_t; \xi_i)\)。
    c. 采样随机噪声：\(\zeta_t \sim \mathcal{N}(0, I)\)。
    d. 更新动量：\(v_{t+1} = (1 - \gamma \eta) v_t - \eta g_t + \sqrt{2 \gamma \eta} \, \zeta_t\)。
    e. 更新参数：\(x_{t+1} = x_t + v_{t+1}\)。
  3. 输出：最终的参数 \(x_T\)，或者最后若干次迭代的平均（Polyak-Ruppert平均）以提高稳定性。

核心思想总结：
SGHMC通过为优化过程引入物理动量 和受控的随机噪声，将SGD转变为一个在参数空间中进行持续探索的过程。动量项帮助平滑梯度噪声的影响，加速收敛；而人为注入的噪声（与摩擦系数平衡）使算法不满足于停留在某个临界点附近，而是以一定的概率“穿越”能量壁垒，从而有可能找到更优的解区域。这种方法特别适合非凸、有噪声的随机优化场景。

基于随机梯度哈密顿蒙特卡洛(SGHMC)的非凸随机优化问题 问题描述 考虑如下随机优化问题： \[ \min_ {x \in \mathbb{R}^n} F(x) = \mathbb{E}_ {\xi \sim \mathcal{D}} [ f(x; \xi) ] \] 其中，目标函数 \(F(x)\) 是 非凸的 ，并且我们无法直接计算其精确梯度 \(\nabla F(x)\)，只能通过从分布 \(\mathcal{D}\) 中采样小批量（mini-batch）数据 \(\{\xi_ i\}\) 来获得 随机梯度估计 \(g(x) = \frac{1}{m} \sum_ {i=1}^m \nabla f(x; \xi_ i)\)。这里的随机性可能来源于数据噪声、随机子采样等。我们的目标是找到一个 近似稳定点 （例如，梯度范数较小的点），而非全局最优解。由于目标函数非凸且存在随机噪声，传统的随机梯度下降（SGD）可能收敛缓慢或陷入较差的局部极小点。本题将介绍 随机梯度哈密顿蒙特卡洛 方法，它通过引入 动量 和 可控的随机噪声 ，帮助优化过程更好地探索参数空间，并可能逃离较差的局部极小点。 解题过程 我们将循序渐进地讲解SGHMC的核心思想和算法步骤，重点是其如何结合 物理学中的哈密顿动力学 和 随机优化 的需求。 步骤1：理解问题与SGD的局限性 问题本质 ：我们需要最小化一个期望风险 \(F(x)\)，但只能通过 随机梯度 \(g(x)\) 来获取信息。\(g(x)\) 是真实梯度 \(\nabla F(x)\) 的无偏估计，但带有方差 \(\mathbb{E}[ ||g(x) - \nabla F(x)||^2 ] = \sigma^2 > 0\)。 SGD更新 ：\(x_ {t+1} = x_ t - \eta g(x_ t)\)，其中 \(\eta\) 是学习率。 局限性 ： 方差导致震荡 ：梯度估计的噪声会导致优化路径在最小值附近剧烈震荡，需要小心翼翼地降低学习率。 难以逃离局部极小点 ：在非凸问题上，SGD容易陷入“尖锐”的局部极小点或鞍点，缺乏有效的探索机制。 步骤2：从物理学中汲取灵感——哈密顿动力学 核心类比 ：将参数 \(x\) 视为一个 粒子 的 位置 。在物理学中，要描述一个粒子的运动，不仅需要位置 \(x\)，还需要 动量 \(v\)（通常 \(v = \dot{x}\)，即速度）。 哈密顿量 ：一个封闭物理系统的总能量 \(H(x, v)\) 通常是守恒的。一个常见的简单形式是动能加势能：\(H(x, v) = U(x) + \frac{1}{2} v^T M^{-1} v\)。 在本优化问题中，我们令 势能 \(U(x) = F(x)\)，即我们的目标函数。 \(M\) 是质量矩阵（常设为单位矩阵 \(I\)），因此哈密顿量简化为 \(H(x, v) = F(x) + \frac{1}{2} v^T v\)。 哈密顿方程 ：描述系统随时间的演化。 \[ \begin{aligned} dx &= v \, dt \\ dv &= -\nabla F(x) \, dt \end{aligned} \] 这组方程定义了一个 确定性 的轨迹。如果沿着这个轨迹运动，系统的总能量 \(H(x, v)\) 是守恒的。这意味着粒子可以在势能面 \(F(x)\) 上高效地“滑动”，而不会像SGD那样因摩擦力而停止在最近的谷底。 步骤3：引入摩擦与噪声——朗之万动力学 纯哈密顿动力学是保守的，会永远运动下去。为了使其最终稳定在低势能区域（即函数最小值附近），我们需要引入： 摩擦项 ：\(-\gamma v \, dt\)，其中 \(\gamma > 0\) 是摩擦系数。这会使动量 \(v\) 衰减，模拟能量耗散。 随机噪声项 ：\(\sqrt{2\gamma} \, dW\)，其中 \(W\) 是标准的维纳过程（布朗运动），\(\sqrt{2\gamma}\) 是噪声强度。这是 涨落耗散定理 的要求，确保系统最终能收敛到正确的 平稳分布 （玻尔兹曼分布）。 于是，我们得到 朗之万动力学 方程： \[ \begin{aligned} dx &= v \, dt \\ dv &= -\nabla F(x) \, dt - \gamma v \, dt + \sqrt{2\gamma} \, dW \end{aligned} \] 可以证明，这个随机微分方程（SDE）的平稳分布是 \(p(x, v) \propto \exp(-H(x, v)) = \exp(-F(x) - \frac{1}{2} v^T v)\)。这意味着，在长期运行下，系统在 \((x, v)\) 空间中的分布正比于 \(e^{-F(x)}\)。因此， 采样得到的 \(x\) 更有可能出现在 \(F(x)\) 值较小的区域 。 步骤4：处理随机梯度——随机梯度朗之万动力学（SGLD）的启发 在实际优化中，我们无法得到精确的 \(\nabla F(x)\)，只有带噪声的估计 \(g(x)\)。如果我们简单地把 \(g(x)\) 代入朗之万方程： \[ dv = -g(x) \, dt - \gamma v \, dt + \sqrt{2\gamma} \, dW \] 那么此时更新中的随机性来源有两个： 算法固有的、我们有意引入的噪声 \(\sqrt{2\gamma} dW\)。 梯度估计本身引入的噪声 \(g(x) - \nabla F(x)\)。 如果这两部分噪声不匹配，系统就不会收敛到我们期望的平稳分布 \(e^{-H(x, v)}\)。梯度噪声的协方差记为 \(B(x)\)。为了使系统仍能收敛到正确的分布，我们需要在更新方程中引入一个额外的“修正”项。 步骤5：随机梯度哈密顿蒙特卡洛（SGHMC）算法 SGHMC巧妙地将梯度估计的噪声 \(B(x)\) 融入到动力学方程中，并重新定义噪声项，使得最终的离散时间算法形式简单。 修改的动力学方程 ：将朗之万方程改写为： \[ \begin{aligned} dx &= v \, dt \\ dv &= -\color{blue}{g(x)} \, dt - \gamma v \, dt + \sqrt{2\gamma} \, dW \\ &= -\color{blue}{\nabla F(x)} \, dt - \gamma v \, dt + \sqrt{2\gamma} \, dW + \color{red}{(\nabla F(x) - g(x)) \, dt} \end{aligned} \] 注意，红色部分 \(\nabla F(x) - g(x)\) 是均值为0、协方差为 \(B(x)\) 的随机噪声。 关键简化 ：我们 不 去显式估计复杂的 \(B(x)\)。相反，我们注意到， 如果我们能精确知道摩擦系数 \(\gamma\) 和梯度噪声的协方差 \(B\) 满足 \(B = 2\gamma I\)，那么红色部分的噪声和固有的 \(dW\) 噪声就可以合并 。实际上，我们不知道 \(B\)。SGHMC的策略是 设置一个足够大的摩擦系数 \(\hat{\gamma}\) ，并定义 \(\beta = \hat{\gamma} I\)。然后，我们将动力学方程重写为： \[ \begin{aligned} dx &= v \, dt \\ M dv &= -\nabla F(x) \, dt - \beta v \, dt + \sqrt{2(\beta - \hat{B})} \, dW + \color{green}{\sqrt{2\hat{B}} \, dW} \end{aligned} \] 其中，\(\hat{B}\) 是 我们对梯度噪声协方差 \(B\) 的一个估计 （实践中常设为 \(\hat{B} = \frac{1}{2} \eta V\)，\(V\) 是梯度方差估计，\(\eta\) 是学习率），绿色部分是 新增的可控噪声项 。 最终算法（离散化形式） ：对上述SDE应用欧拉-丸山离散化，并定义学习率 \(\epsilon = \eta\)，质量矩阵 \(M = I\)，得到SGHMC的迭代更新公式： \[ \begin{aligned} v_ {t+1} &= (1 - \gamma \eta) v_ t - \eta g(x_ t) + \sqrt{2\gamma \eta} \, \mathcal{N}(0, I) \\ x_ {t+1} &= x_ t + v_ {t+1} \end{aligned} \] 其中： \(v_ t\) 是动量（速度）。 \(\gamma\) 是摩擦系数（控制动量衰减）。 \(\eta\) 是学习率（离散化的步长）。 \(g(x_ t)\) 是在 \(x_ t\) 处计算的随机梯度。 \(\mathcal{N}(0, I)\) 是标准正态分布噪声。 注意，在标准实现中，为了简化，常常会忽略对 \(\hat{B}\) 的显式估计，而将 \(\gamma\) 当作一个可以调节的超参数。 算法中人为注入的噪声 \(\sqrt{2\gamma \eta} \, \mathcal{N}(0, I)\) 正是SGHMC帮助逃离局部极小点、进行探索的关键 。 步骤6：算法步骤总结 初始化 ：选择初始参数 \(x_ 0\)，初始动量 \(v_ 0 = 0\)。设定超参数：学习率 \(\eta\)，摩擦系数 \(\gamma\)，总迭代步数 \(T\)，批量大小 \(m\)。 迭代循环 （对于 \(t = 0, 1, ..., T-1\)）： a. 采样 ：从数据分布 \(\mathcal{D}\) 中随机抽取一个小批量样本 \(\{\xi_ 1, ..., \xi_ m\}\)。 b. 计算随机梯度 ：\(g_ t = \frac{1}{m} \sum_ {i=1}^m \nabla f(x_ t; \xi_ i)\)。 c. 采样随机噪声 ：\(\zeta_ t \sim \mathcal{N}(0, I)\)。 d. 更新动量 ：\(v_ {t+1} = (1 - \gamma \eta) v_ t - \eta g_ t + \sqrt{2 \gamma \eta} \, \zeta_ t\)。 e. 更新参数 ：\(x_ {t+1} = x_ t + v_ {t+1}\)。 输出 ：最终的参数 \(x_ T\)，或者最后若干次迭代的平均（Polyak-Ruppert平均）以提高稳定性。 核心思想总结 ： SGHMC通过为优化过程引入 物理动量 和 受控的随机噪声 ，将SGD转变为一个在参数空间中进行 持续探索 的过程。动量项帮助平滑梯度噪声的影响，加速收敛；而人为注入的噪声（与摩擦系数平衡）使算法不满足于停留在某个临界点附近，而是以一定的概率“穿越”能量壁垒，从而有可能找到更优的解区域。这种方法特别适合 非凸、有噪声 的随机优化场景。