1. 统计推断-基于神经网络与Langevin扩散的自适应潜变量建模与优化

• 在先验分布估计方面，借鉴了机器学习文献中各种深度生成模型的思想，通过深度神经网络将低维的已知分布变换为潜变量分布 $\pi (u)$ 。这种方式使潜变量分布具有高度灵活性，从而用户无需在受限的参数形式下显式定义 $\pi (u)$ 。
• 在模型参数估计方面，已有的两种方法是变分自编码器和交替反向传播 ，旨在优化某个目标函数，该目标函数可视为对真实似然函数的近似：前者利用变分下界，后者依赖蒙特卡洛采样。相比之下，我们的贡献在于：

1. 提出了一种基于真实似然函数的估计方法并具备坚实的收敛性保证。
2. 设计了一种高效的计算算法，使该框架适用于复杂模型和大规模数据集。更为重要的是，与近似方法相比，该方法并未牺牲统计精度。凭借这些优势，我们的目标是帮助数据分析人员更专注于潜变量建模，而不必再为大规模计算而设计繁琐的定制算法。

该框架之所以能具备上述特性，源于对 Langevin 采样算法与随机逼近方法的系统性分析。这两类技术在潜变量模型中的应用已在文献中得到广泛探讨，但这些工作大多忽略了 Langevin 算法的偏差。本文的创新在于将 Langevin 算法的有限样本偏差纳入优化过程的一部分（研究了 Langevin 算法在有限迭代步和有限样本条件下的偏差以及其与随机逼近方法之间的相互作用）。这一分析为现代统计与机器学习模型中常用的采样与优化技术提供了理论支撑。值得注意的是，这些结果是普遍性的，而不仅限于潜变量模型。

本文其余部分的组织如下：第二节对所提出的推断框架进行概述；第三节介绍相应的计算算法；第四节给出算法的理论保证；第五节和第六节分别通过模拟实验和真实数据分析展示方法的实际表现；第七节进行总结与讨论。

2. 本文提出的框架

统计建模的目标之一是使模型分布 $f_{\theta }(x)$ 与真实数据分布 $f^{\ast }(x)$ 相匹配，其中随机向量 $X\in \mathcal{X}\subseteq {\mathbb{R}}^p$ 。对于潜变量模型，$f_{\theta }(x)$ 通过一个不可观测的潜变量 $U\in \mathcal{U}\subseteq {\mathbb{R}}^d$ 定义为
$$(1)f_{\theta }(x)={\int }_{\mathcal{U}}{f}_{\theta }(x\mid u)\pi (u)\mathrm{d}u$$其中我们假设数据点 $X$ 是由潜变量 $U$ 生成的，生成机制遵循已知的条件分布 $f_{\theta }(x\mid u)$ ，其形式中可能包含未知参数向量 $\theta$ ，并且 $U$具有边际分布 $\pi (u)$ 。尽管形式简单，公式（1）实际上包含了许多著名且功能强大的统计模型作为特殊情形。

2.1. 隐式模型表示 $\pi (u)$

我们假设潜变量分布 $\pi$ 完全未知（如果已有 $\pi$ 的先验知识也可以很容易将其纳入模型）。然而直接寻找一个密度估计器 $\hat{\pi }$ 非常困难，因为这需要在一个给定的分布族中优化 $\pi$ 的某个特定函数；由于多维密度函数受约束较多，因此考虑对 $\pi$ 进行另一种表示方式。

我们通过一个灵活的隐式模型来表示 $\pi (u)$ ，采用两步法：

1. 首先，从一个固定的分布 ${\pi }_0(z)$ 中采样一个变量 $Z\in \mathcal{Z}\subseteq {\mathbb{R}}^r$ 。这里，${\pi }_0$ 可以是任何常见分布，例如标准正态 $N\left(0,{\mathbf{I}}_r\right)$ 。
2. 然后，通过一个确定性变换 $h:\mathcal{Z}\to \mathcal{U}$ 将 $Z$ 映射到潜变量 $U=h(Z)$ 。这种表示方法在深度学习中很常见，比如生成对抗网络和变分自编码器（VAE）。当 $U$ 和 $Z$ 是一维变量时，这种映射就称为概率积分变换。

这种方法的优点包括：（a）$h(Z)$ 可以表示非常广泛的分布；（b）在许多情况下，潜变量 $U$ 实际上只位于一个低维流形上，因此 $Z$ 的维度 $r$ 可以远小于 $U$ 的维度 $d$ ；（c）映射 $h$ 没有额外约束，这大大简化了优化问题。

通过 $U=h(Z)$ ，可以定义 $(X,Z)$ 的联合分布：$$Z\sim {\pi }_0,X\mid Z\sim p_{X\mid Z}(x\mid z)=f_{\theta }(x\mid h(z))$$其中 $p_{X\mid Z}(x\mid z)$ 是给定 $Z=z$ 时 $X$ 的条件密度，$p_{Z\mid X}(z\mid x)$ 则表示给定 $X=x$ 时 $Z$ 的条件密度。在这种表述下，潜变量模型的推断问题就转化为三个任务：1．估计映射 $h$；2．估计参数向量 $\theta$；3．估计条件分布 $p_{Z\mid X}(z\mid x)$。

2.2. 神经网络拟合

过去十年中，深度神经网络（DNN）在复杂人工智能任务中表现出巨大成功，例如图像识别和语音理解。在统计领域，也有越来越多的工作使用 DNN 来解决问题，例如变量选择。为了利用 DNN 的优势，我们用一个 DNN $h_{\eta }$ 来近似映射 $h$ ，其定义如下：
$$\begin{array}{rl}& (2)h_{\eta }=T_D\circ \sigma \circ T_{D-1}\circ \cdots \circ \sigma \circ T_1\\ & (3)T_1\in {\mathcal{A}}_d^N,T_2,\dots ,T_{D-1}\in {\mathcal{A}}_N^N,T_D\in {\mathcal{A}}_N^r\\ & \sigma ((x_1,\dots ,x_N{)}^{\mathrm{T}})={\left(\sigma \left(x_1\right),\dots ,\sigma \left(x_N\right)\right)}^{\mathrm{T}}\end{array}$$其中，$N$ 和 $D$ 分别表示 DNN 的宽度和深度， ${\mathcal{A}}_{n_1}^{n_2}$ 表示从 ${\mathbb{R}}^{n_1}$ 到 ${\mathbb{R}}^{n_2}$ 的仿射变换，其形式为$x_{n_1\times 1}\mapsto W_{n_2\times n_1}x+b_{n_2\times 1}$, $\sigma$ 是非线性激活函数，$\eta$ 表示 DNN 中的全部参数。在公式（3）中，每一层使用固定宽度 $N$ 仅是为了简化符号；在实际应用中，每一层 $T_i$ 的维度可以不同。DNN 模型在选择激活函数上具有很大的灵活性。常用的激活函数包括整流线性单元（ReLU，Glorot，Bordes，and Bengio，2011）， σ ( x ) = max ⁡ { x , 0 } \sigma(x)=\max \{x, 0\} σ(x)=max{x,0} ，以及其平滑版本 σ ( x ) = log ⁡ ( 1 + exp ⁡ ( x ) ) ≈ max ⁡ { x , 0 } \sigma(x)=\log (1+\exp (x)) \approx \max \{x, 0\} σ(x)=log(1+exp(x))≈max{x,0} 。

DNN 的有效性在很大程度上源于其通用逼近性质。为了简化表示，我们将映射 $h$ 视为一个函数向量，$h=\left(h_1,\dots ,h_d\right)$ ，其中每个 $h_i$ 是定义在 ${\mathbb{R}}^r$ 上的 $r$ 维连续函数；同样地，根据输出维度，$h_{\eta }$ 可分解为 $h_{\mathbf{\eta }}=\left(h_{1,\eta },\dots ,h_{d,\eta }\right)$ 。Leshno 等人（1993）的定理 1 表明：只要隐含单元数足够，每个 $h_{\eta }$ 的分量 $h_{i,\eta }$ 在连续函数类中是稠密的，前提是激活函数 $\sigma (x)$ 连续、局部有界，且不是多项式函数。大多数常用的激活函数，如 σ ( x ) = max ⁡ { x , 0 } \sigma(x)=\max \{x, 0\} σ(x)=max{x,0} 和 $\sigma (x)=\log (1+\exp (x))$ ，都满足这一条件。因此，DNN 模型 $h_{\mathbf{\eta }}$（公式（2））可以以任意预设精度逼近任意连续函数 $h$ 。当然，在有限样本情况下，需要对 $h_{\eta }$ 的宽度和深度进行调节，以在偏差与方差之间取得平衡。本文不讨论这种超参数的调节，而是关注在有限 $N$ 和 $D$ 下，如何求解最优参数 $\eta$ ，因为本文的结论适用于任意固定架构的 DNN 模型。

如果 DNN 参数 $\eta$ 和条件分布参数 $\mathbf{\theta }$ 已知，则观测数据的边际密度函数可以表示为
$$(4)f_{\beta }(x)={\int }_{\mathcal{U}}{f}_{\theta }(x\mid u){\pi }_{\eta }(u)\mathrm{d}u={\int }_{\mathcal{Z}}{f}_{\theta }\left(x\mid h_{\eta }(z)\right){\pi }_0(z)\mathrm{d}z$$其中 $\beta =(\eta ,\theta )$ 包含模型中的所有参数。在本文的后续部分，如果不会引起混淆，我们将省略积分区域的符号。

给定观测数据 $X_1,\dots ,X_n$，感兴趣的参数 $\beta$ 可以根据最大似然原则进行估计：
$$(5)\hat{\beta }=\underset{\beta }{\mathrm{argmax}}\frac{1}{n}\sum _{i=1}^n\ell \left(\beta ;X_i\right)=\underset{\beta }{\mathrm{argmax}}\frac{1}{n}\sum _{i=1}^n\log \left[f_{\beta }\left(X_i\right)\right],$$其中 $\ell (\beta ;x)=\log \left[f_{\beta }(x)\right]$ 是对数似然函数。

3．算法

优化问题（5）的主要困难在于计算 $\ell (\beta ;x)$ 及其关于 $\beta$ 的梯度。对于固定的 $x$ ，精确计算 $f_{\beta }(x)$ 是不可行的，因为一般来说积分（4）是高维的且没有闭式解。因此，如何高效地计算估计量 $\hat{\beta }$ 成为本文剩余部分的主要关注点。

3.1 随机梯度方法

为简化讨论，我们先考虑对 $\ell (\beta ;x)$ 的单个固定 $x$ 的处理，因为方法可以很容易推广到求和形式 $n^{-1}{\sum }_{i=1}^n\ell \left(\beta ;X_i\right)$ 。

当无法精确计算对数似然函数 $$(\ast )\ell (\beta ;x)=\log \left[f_{\beta }(x)\right]=\log {\int }_{\mathcal{Z}}{f}_{\theta }\left(x\mid h_{\eta }(z)\right){\pi }_0(z)\mathrm{d}z$$ 时，一种常见的做法是找到一个＂低成本＂的随机量来近似它，然后使用随机逼近方法进行优化。这种方法通常称为随机梯度方法。

例如：很容易看出下式是 $f_{\beta }(x)$ 的无偏估计 $$(6){\tilde{f}}_{\beta }(x,V)=\frac{1}{M}\sum _{i=1}^M{f}_{\theta }\left(x\mid h_{\eta }\left(V_i\right)\right),V_1,\dots ,V_M\stackrel{\mathrm{iid}}{\sim }{\pi }_0$$其中 $V=\left(V_1,\dots ,V_M\right)$ 为蒙特卡洛样本，$M$ 为样本量。

评价：原则上，$\tilde{\ell }(\beta ;x,V)=\log \left[{\tilde{f}}_{\beta }(x,V)\right]$ 可以通过足够大的 $M$ 任意精确地逼近 $\ell (\beta ;x)$ ，且随机梯度方法收敛的条件已有充分研究。然而，尽管理论上简单，但公式（6）通常方差较高，因为大部分求和项接近于零。这会导致收敛速度极慢，因此直接使用 $\tilde{\ell }(\beta ;x,V)$ 并不实用。

我们引入一个在优化方法中起核心作用的函数。考虑对对数似然函数 $(\ast ):\ell (\beta ;x)=\log \left[f_{\beta }(x)\right]$的一个代理函数，定义为
$$(7)\mathcal{L}(\beta ,\tilde{\beta };x)=\int \log \left[\frac{f_{\theta }\left(x\mid h_{\eta }(z)\right){\pi }_0(z)}{p_{\tilde{\beta }}(z\mid x)}\right]p_{\tilde{\beta }}(z\mid x)\mathrm{d}z$$其中 $p_{\tilde{\beta }}$ 表示当 $\beta$ 固定为 $\tilde{\beta }$ 时的条件密度 $p_{Z\mid X}(z\mid x)=f_{\theta }(x\mid h_{\mathbf{\eta }}(z))$ 。

评价：函数 $\mathcal{L}(\beta ,\tilde{\beta };x)$ 具有以下简单而重要的性质：引理1． $\mathcal{L}(\beta ,\tilde{\beta };x)$ 是 $\ell (\beta ;x)$ 的一致下界，即对任意 $\tilde{\beta }$ ，有$\mathcal{L}(\beta ,\tilde{\beta };x)\le \ell (\beta ;x)$。此外，当 $\tilde{\beta }=\beta$ 时，$\mathcal{L}(\beta ,\beta ;x)=\ell (\beta ;x)$。

需要注意的是，公式（7）是变分推断和深度学习文献中证据下界的一个特例。然而，一般的 ELBO 只满足引理1的第一个性质，而额外的恒等式 $\mathcal{L}(\beta ,\beta ;x)=\ell (\beta ;x)$ 是本文算法的关键。

迭代更新方法：定义 $\mathcal{L}(\beta ,\tilde{\beta };x)$ 对 $\beta$ 的偏导为$g(\beta ,\tilde{\beta };x)=\frac{\partial \mathcal{L}(\beta ,\tilde{\beta };x)}{\partial \beta }$，则我们的迭代更新方法为$${\beta }_{t+1}={\beta }_t+{\alpha }_t\cdot g_t\left({\beta }_t;x\right)$$其中 ${\beta }_t=\left({\eta }_t,{\theta }_t\right)$ 是第 $t$ 次迭代的估计值，$t\in \mathbb{N},{\alpha }_t$ 为步长，且$$(8)g_t(\beta ;x)=g\left(\beta ,{\beta }_t;x\right)=\int \frac{\partial \log \left[f_{\theta }\left(x\mid h_{\eta }(z)\right)\right]}{\partial \beta }{p}_{{\beta }_t}(z\mid x)\mathrm{d}z$$为梯度。

随机梯度：然而，精确计算 $(8):g_t(\beta ;x)$ 是不可行的，因此我们使用随机梯度 $\tilde{g}\left(\beta ;x,W_t\right)$ 进行替代：
$$(9){\beta }_{t+1}={\beta }_t+{\alpha }_t\cdot \tilde{g}\left({\beta }_t;x,W_t\right)$$其中 $W_t=\left(W_t^{(1)},\dots ,W_t^{\left(M_t\right)}\right)$ 为一系列随机向量，$M_t$ 为序列长度。一个常见选择是从 $p_{{\beta }_t}(z\mid x)$ 中采样足够大的马尔可夫链蒙特卡洛样本，使得 $\tilde{g}\left(\beta ;x,W_t\right)$ 近似无偏地估计 $g_t(\beta ;x)$ 。然而，通常很难事先确定 MCMC 样本量，因此在下一节中我们将介绍一种对超参数更透明的替代采样策略。

3.2. Langevin 算法计算 $\tilde{g}\left(\beta ;x,W_t\right)$

在所提出的 ALMOND 框架中，随机梯度 $\tilde{g}\left(\beta ;x,W_t\right)$ 的计算基于 Langevin 算法。

• 从整体上看，Langevin 算法是一种从目标分布中获取近似样本的方法。
• 基于 Langevin 采样，我们得到的 $\tilde{g}\left(\beta ;x,W_t\right)$ 是 $g_t(\beta ;x)$ 的 有偏估计，但其偏差可以通过增加计算时间而变得任意小。
• 直观上，在优化初期我们只需粗略的近似来节省计算开销，同时参数仍然可能沿着上升方向前进；而随着优化的进行，会使用更高精度的采样来保证在最优点附近的梯度更新更加准确。

这些启发性解释将在下一节得到严格化：我们将仔细研究 Langevin 算法中的偏差来源，提出控制偏差的方法，并讨论在何种条件下式（9）的迭代能够保证收敛。下一节给出了 Langevin 算法的必要背景知识，

而这里我们主要描述其计算方法: 目标是$\tilde{g}\left(\beta ;x,W_t\right)$ 近似无偏地估计 $g_t(\beta ;x)=\int \frac{\partial \log \left[f_{\theta }\left(x\mid h_{\eta }(z)\right)\right]}{\partial \beta }{p}_{{\beta }_t}(z\mid x)\mathrm{d}z$

定义
$$v_t(z)=\frac{\partial }{\partial z}\log \left[f_{{\theta }_t}\left(x\mid h_{{\eta }_t}(z)\right){\pi }_0(z)\right],G(\beta ;x,z)=\frac{\partial }{\partial \beta }\log \left[f_{\theta }\left(x\mid h_{\eta }(z)\right)\right],$$则有下式可以近似$g_t(\beta ;x)$,$$\begin{gathered} (10)\tilde{g}\left(\beta ;x,W_t\right)=\frac{1}{M_t}\sum _{i=1}^{M_t}G\left(\beta ;x,W_t^{(i)}\right)=\frac{1}{M_t}\sum _{i=1}^{M_t}\frac{\partial }{\partial \beta }\log \left[f_{\theta }\left(x\mid h_{\eta }(W_t^{(i)})\right)\right]. \\ (11)W_t^{(k)}\leftarrow W_t^{(k-1)}+{\gamma }_t\cdot v_t\left(W_t^{(k-1)}\right)+\sqrt{2{\gamma }_t}\cdot {\xi }^{(k)} \end{gathered}$$(11) 式里面$W_t^{(k)}$迭代公式的目标分布是如下条件分布：$$p_{{\beta }_t}(z\mid x)\propto f_{{\theta }_t}\left(x\mid h_{{\eta }_t}(z)\right){\pi }_0(z)$$于是 Langevin 动力学写成：$$W_t^{(k+1)}=W_t^{(k)}+{\gamma }_t\cdot {\nabla }_z\log \left(f_{{\theta }_t}\left(x\mid h_{{\eta }_t}\left(W_t^{(k)}\right)\right){\pi }_0\left(W_t^{(k)}\right)\right)+\sqrt{2{\gamma }_t}{\xi }^{(k)}$$ 由于我们设计的 drift 使用了 ${\nabla }_z\log p_{{\beta }_t}(z\mid x)=:v_t(z)$ ，那么其平稳分布自然就是$p_{{\beta }_t}(z\mid x)$ 。

接下来，我们描述随机向量序列 $W_t$ 的计算方法，以及整个随机梯度方法的概要（见 算法 1）。

我们的方法可以与 交替反向传播算法（Han et al. 2017, 算法 1）进行比较，但存在一个关键区别：在我们的 Langevin 算法中，参数 ${\gamma }_t$ 和 $M_t$ 随迭代次数 $t$ 而变化，而不是取固定值。这一点事实证明是整个优化过程最终收敛的关键。相关讨论将在 第 4.3 节展开。在附录 A 中，我们介绍了若干有助于实际实现的技巧，包括通过 VAE 计算得到的初始值 $({\beta }_0,q_0)$，以及 Langevin 算法的 warm-start 策略。在实际应用中，这些技巧能够显著加速优化收敛，并减少计算成本。