PaperNotes(8)-Stein Variational Gradient Descent A General Purpose Bayesian Inference Algorithm

SVGD

Abstract

本文提出了一种通用的变分推理算法，这算法与优化中应用的梯度下降法相对应。我们的方法通过应用函数梯度下降算法将kl散度最小化，从而迭代地传输一组元素用以匹配目标分布。我们对各种现实模型和数据集进行了实证研究（Empirical studies），在这些模型和数据集上，我们的方法与现有最先进的方法能够一比高下。我们方法的推导是基于一个新的理论结果：（1）它将光滑变换下的kl散度导数与Stein的恒等式联系起来，（2）以及最近提出的一个具有独特意义的Kernelized Stein差分。

1 Introduction

贝叶斯推理为复杂数据建模和不确定性推理提供了强有力的工具，但是贝叶斯推断的致命缺点是：难以准确计算后验概率。==马尔可夫链蒙特卡罗（MCMC）==已被广泛应用于后验样本的近似提取，但往往速度慢，收敛困难。变分推理将贝叶斯推理问题框架化为确定性优化问题，通过最小化其KL散度，使目标分布近似为更简单的分布。这使得变分方法可以通过使用现成的优化技术有效地解决贝叶斯推断问题，并且使用随机梯度下降技巧（例如1）很容易适用于大数据集（即“大数据”）。相比之下，将MCMC扩展到大数据设置更具挑战性（参见例如2、3）。

同时，变分推理的准确度和计算成本都取决于定要近似的分布。简单的近似集，如传统的 平均场方法 中使用的近似集，限制性太强，无法与真实的后验分布相似，而更高级的选择给后续的优化任务带来了更多的困难。因此，有效的变分方法往往需要在基于模型的基础上推导出来，这是开发通用的、用户友好的、适用于不同类型模型、和应用领域的非ML专家可以获取的变分工具的一个主要障碍。

这种情况与寻找后验模式的最大后验（MAP）优化任务（有时被称为 poor man’s Bayesian estimator，与 full Bayesian inference相对应–用于估计全后验概率）形成了对比。对于这种情况，（随机）梯度下降的变体是一种简单的、通用的，然而强大的工具。最近人们对创建用户友好的变分推理工具（如4-7）越来越感兴趣，但仍需要更多的努力来开发更有效的通用算法。
在本文工作中，我们提出了一种新的通用变分推理算法，它可以作为应用于全贝叶斯推理中的梯度下降的自然对应（见算法1）。我们的算法使用一组 particles进行近似，在这组 particles上执行（函数）梯度下降以最小化kl散度，并驱动 particles拟合真实的后验分布。我们的算法形式简单，可以在任何时候应用(任何可以应用梯度下降的地方应用)。事实上，当只使用一个 particle时，它会降低梯度下降到MAP的，而更多粒子时会自动变为全贝叶斯采样方法。
在我们的算法的基础上是一个新的理论结果：它连接了kl散度导数w.r.t .平滑变量变换 和最近引入的Kernelized Stein差分[8–10]，这使得我们能够 在复制核希尔伯特空间（RKhs）的单位球内 为 使kl散度下降最快 的 最佳平滑扰动方向 推导出一个闭式解。这一新的结果具有独立的意义，可以在机器学习和统计得到广泛的应用(不仅仅在变分推理的学中)。

2 Background

预备知识：
x是$\mathcal{X}\subset {\mathbb{R}}^d$中的一个随机变量， { D k } \left\{D_{k}\right\} {Dk​}是独立同分布的观测集合。x的贝叶斯推断涉及到：在先验概率$p_0(x)$的条件下，后验分布的推断$p(x):=\overline{p}(x)/Z$，其中$\overline{p}(x):=p_0(x){\prod }_{k=1}^{N}p\left(D_k|x\right)$，$Z=\int \overline{p}(x)dx$是一个十分难计算的归一化常数。为了方便起见，我们已经放弃了p(x)中数据 { D k } \left\{D_{k}\right\} {Dk​}的调节。

让 $k\left(x,x^{\prime }\right):\mathcal{X}\times \mathcal{X}\to \mathbb{R}$是一个正定核，$k\left(x,x^{\prime }\right)$的再生产核希尔伯特空间$(\mathrm{R}\mathrm{K}\mathrm{H}\mathrm{S})\mathcal{H}$是 { f : f ( x ) = ∑ i = 1 m a i k ( x , x i ) , a i ∈ R , m ∈ \left\{f : f(x)=\sum_{i=1}^{m} a_{i} k\left(x, x_{i}\right), \quad a_{i} \in \mathbb{R}, \quad m \in\right. {f:f(x)=∑i=1m​ai​k(x,xi​),ai​∈R,m∈$\mathbb{N},x_i\in \mathcal{X}\}$线性闭合空间。这个空间具有内积计算：$⟨f,g{⟩}_{\mathcal{H}}={\sum }_{ij}{a}_i{b}_{j}k\left(x_i,x_j\right)$ for $g(x)={\sum }_i{b}_{i}k\left(x,x_i\right)$。表示着在${\mathcal{H}}^d$中，向量函数
$\mathbf{f}={\left[f_1,\dots ,f_d\right]}^{\top }$ with $f_i\in \mathcal{H}$的空间，可以进行内积计算$⟨\mathbf{f},\mathbf{g}{⟩}_{{\mathcal{H}}^d}={\sum }_{i=1}^d{⟨f_i,g_i⟩}_{\mathcal{H}}$。我们假设所有向量都是列向量。让${\nabla }_x\mathbf{f}=\left[{\nabla }_x{f}_1,\dots ,{\nabla }_x{f}_d\right]$。

Stein’s Identity and Kernelized Stein Discrepancy
Stein’s Identity在我们的框架中扮演着一个基本角色，让$p(x)$是$\mathcal{X}\subseteq {\mathbb{R}}^d$上的一个连续可微分（也称为光滑）密度。$\varphi (x)={\left[{\varphi }_1(x),\cdots ,{\varphi }_d(x)\right]}^{\top }$光滑的向量函数。 Stein’s identity表明对于充分规则$\varphi$有如下等式成立：
$${\mathbb{E}}_{x\sim p}\left[{\mathcal{A}}_p\varphi (x)\right]=0,\text{ where }{\mathcal{A}}_p\varphi (x)={\nabla }_x\log p(x)\varphi (x{)}^{\top }+{\nabla }_x\varphi (x)$$
其中：${\mathcal{A}}_p$为Stein operator，其作用于函数φ并产生一个零均值函数。。。。。

3 Variational Inference Using Smooth Transforms

变分推理使用更简单的分布$q^{\ast }(x)$近似目标分布$p(x)$。$q^{\ast }(x)$在预先定义的集合 Q = { q ( x ) } \mathcal{Q}=\{q(x)\} Q={q(x)}中通过最小化如下KL散度找到：
q ∗ = arg ⁡ min ⁡ q ∈ Q { K L ( q ∥ p ) ≡ E q [ log ⁡ q ( x ) ] − E q [ log ⁡ p ‾ ( x ) ] + log ⁡ Z } q^{*}=\underset{q \in \mathcal{Q}}{\arg \min }\left\{\mathrm{KL}(q \| p) \equiv \mathbb{E}_{q}[\log q(x)]-\mathbb{E}_{q}[\log \overline{p}(x)]+\log Z\right\} q∗=q∈Qargmin​{KL(q∥p)≡Eq​[logq(x)]−Eq​[logp​(x)]+logZ}

其中，在解决优化问题时，我们不需要计算常数$logZ$。集合$\mathcal{Q}$的选择是很关键的，不同的选择定义了不同类型的变分推理方法。最好的$\mathcal{Q}$的集合应该是在以下三个因素中寻找平衡：i）精度，足够宽，能够接近大类目标分布；ii）可牵引性，由易于推断的简单分布组成；以及iii）可解性，以便有效地解决后续的kl最小化问题。

在我们的工作中，由可追踪的推断分布中获取的分布组成了集合$\mathcal{Q}$，也就是说，$\mathcal{Q}$是由来自$z=T(x)$的随机变量组成的分布组成，其中$T:\mathcal{X}->\mathcal{X}$是光滑的一到一的映射；x是从可追踪的推断分布$q_0(x)$中抽取出来的。通过改变变量方程，z的变量概率分布为：

××××××××××这样的q原则上可以近似任意的分布：可以证明，在没有atoms（即，没有一个点携带正质量）的任何两个分布之间总是存在一个可测量的变换T；此外，对于Lipschitz连续密度p和q，在它们之间总是存在着转换，而转换最少与p和q一样流畅。对此主题进行深入讨论，详见我Villani[12]。

然而，在实践中，我们需要适当地限制变换集T，使 (4) 中相应的变分优化切实可行。一种方法是考虑具有一定参数形式的T并优化相应的参数[例如，13、14]。然而，这给选择合适的参数族来平衡精度、可伸缩性和可解性带来了一个困难，特别是考虑到T必须是一一映射，并且雅可比矩阵必须高效可计算.

相反，我们提出了一种 迭代 构造 增量 变换 的新算法，它可以有效地在RKHS中的T上执行 最速下降 。我们的算法不需要显式地指定参数形式，也不需要计算雅可比矩阵，并且有一个特别简单的形式，模仿典型的梯度下降算法，使得它即使对于非专家的变分推理也很容易实现。

3.1 Stein Operator as the Derivative of KL Divergence

为了解释如何将（4）中的kl散度最小化，我们考虑了由 标识映射 的 小扰动 形成的 增量变换 ：$T(x)=x+ϵ\varphi (x)$，其中$\varphi (x)$是一个表征扰动方向的光滑函数，标量$ϵ$表示扰动幅度。当$|ϵ|$充分小的时候，T的雅克比矩阵是满秩的（近似于单位阵），因此，利用 反函数定理 ，T是一一映射。

下面的结果形成了我们的方法的基础，在 Stein算子 与 KL散度关于微小扰动 $ϵ$ 的导数 之间形成了一个有意义的联系。

定理3.1

引理3.2

引理3.2 表面一个可迭代的过程，其将一个原始分布$q_0$转换为目标分布p：我们一开始将$T_0^{\ast }(x)=x+{ϵ}_0\ast {\varphi }_{q_0,p}^{\ast }(x)$应用于$q_0$ ，使得$q_0与p$之间的KL散度下降${ϵ}_0\ast {\mathcal{D}}^2(q_0,p)$，这就是形成了一个新的分布。

函数梯度为了进一步了解这个过程，我们现将重新解释(6)式子为RKHDS中的函数梯度。对于$f\in {\mathcal{H}}^d$任意方程$F[f]$，其梯度为${▽}_{f}F[f]$是${\mathcal{H}}^d$中的方程，例如“$F[f+ϵg(x)=F[f]+ϵ<{▽}_{f}F[f],g>+O({ϵ}^2)$，对于任意的$g\in {\mathcal{H}}^d$和$ϵ\in \mathcal{R}$

3.2 Stein Variational Gradient Descent

为了在实践中实现迭代过程（7），需要近似计算（6）中${\varphi }_{q,p}^{\ast }(x)$的期望值。要做到这一点，我们可以首先从初始分布$q_0$中抽样一组样本 { x i 0 } i = 1 n \{x_i^0\}_{i=1}^n {xi0​}i=1n​,用式（7）的经验版本迭代更新样本，经验版本：streian算子在$q_l$下的期望近似为样本在$l$轮迭代中的经验均值。这过程在算法1中总结出来了。这允许我们（确定性地）传输一组点（来自初始分布）以匹配我们的目标分布$p(x)$，有效地提供$p(x)$提供了一种采样方法（目标是采样？）。我们可以看到，这个过程的实现完全不依赖于初始分布$q_0$，实际上，我们可以从一组任意点 { x i } i = 1 n \{x_i\}_{i=1}^n {xi​}i=1n​开始，这组任意的点可能由复杂（随机或确定性）黑盒过程生成。

我们期望 { x i l } i = 1 n \{x_i^l\}_{i=1}^n {xil​}i=1n​随着n的增加形式越来越接近$q_l$.为此，用$\Phi$表示$q_l->q_{l+1}$非线性映射。写做$q_{l+1}={\Phi }_l(q_l)$，其中 $q_l$通过$q_{l,[T_l^{\ast }]},{\varphi }_{q_l,p}$进入映射$\Phi$。那么算法1 的更新 可以看做${\hat{q}}_{l+1}=\Phi ({\hat{q}}_l)$,也就是从经验测量${\hat{q}}_l->{\hat{q}}_{l+1}$的映射。因为${\hat{q}}_0$随着n增加收敛于$q_0$,当$\Phi$是连续映射时，$\widehat{q_l}$也应当收敛与$q_l$。在相互作用粒子系统的平均场理论（如15）中，已经建立了关于这种收敛性的严格理论结果：。。。此外，对于任何固定的$i_0$，每个粒子$x_{i_0}^l$的分布也趋向于$q_l$，并且与任何其他有限的粒子子集（n→∞）无关，这一现象称为混沌传播[16]。我们为今后的工作留下了具体的理论分析。

算法1 在样本的水平上模仿了梯度动力学,(8)式${\hat{\varphi }}^{\ast }(x)$中的两项起不同的作用:第一项通过平滑的梯度方向将粒子推向$p(x)$的高概率区域，这是由核函数加权的所有点的梯度的加权和。第二项作为一个排斥力(repulsive force)，防止所有点一起崩溃(collapse together)成$p(x)$的局部模式；要看到这一点，考虑RBF内核$k(x,x^{\prime })=exp(-\frac{1}{h}||x-x^{\prime }|{|}^2)$，对$x_j^l$求梯度之后为${\sum }_j\frac{1}{h}(x-x_j)k(x_j,x)$.当$k(x_j,x)$很大时能够将x拉离它的邻居点$x_j$.如果我们让带宽$h\to 0$，排斥项消失，式（8）化简为一组独立的典型梯度上升链，以使$logp(x)$最大化（即MAP：最大后验概率），所有粒子都将崩溃为局部模式。

另一个有趣的现象 ，只使用一个样本（n=1），在这种情况下，算法1将任何满足${\nabla }_{x}k(x,x)=0$的核简化为典型的单链 梯度上升（MAP）。这表明我们的算法可以很好地推广到有监督的学习任务中，即使只有极少量的n个粒子，因为MPA的梯度上升（n=1）已经在实践中证明是非常成功的。这一特性将我们的粒子法与典型的蒙特卡罗方法区分开来，后者需要对多个点进行平均。这里的关键区别在于，我们使用确定性排斥力(deterministic repulsive force)，而不是蒙特卡罗随机性，来获得分布近似的不同点。
Complexity and Efficient Implementation 主要花费在梯度计算上面$p(x)\propto p_0(x){\prod }_{k=1}^{N}p(D_k|x)$当N很大时要计算梯度很费时间。可以通过对数据进行采样，得到子样本小批次集合$\Omega$近似计算：
$${\nabla }_{x}logp(x)\approx log{p}_0(x)+\frac{N}{\Omega }\sum _{k\in \Omega }logp(D_k|x)$$

通过将n个粒子的梯度评价并行化，可以获得额外的加速。

式（8）还需要计算内核矩阵 { k ( x i , x j ) } \{k(x_i,x_j)\} {k(xi​,xj​)}，其成本为$O(n^2)$；在实践中，与梯度评估成本相比，该成本相对较小，因为在实践中使用相对较小的n（例如几百个）就足够了。如果需要非常大的n，可以通过对粒子进行次采样或使用核$k(x,x_0)[17]的随机特征展开来近似求式子(8)中的和。

4 Related Works

5 Experiments

我们在toy和现实世界的例子中测试了我们的算法，在这些例子中我们发现我们的方法优于各种基线方法。我们的代码可以在https://github.com/dartml/stein找到。

对于我们所有的实验，我们使用RBF核，把带宽设为$h=me{d}^2/logn$ ,med 是点与点之间距离的中点（？）；基于（一个式子） 启发。所以对于每一个x$x_i$，来自它自身梯度的贡献，以及来自其他点的影响，彼此平衡。注意，通过这种方式，带宽h在迭代过程中实际上是自适应变化的。我们使用adagrad作为步长，并使用先验分布初始化粒子。 除非另有规定。

Toy Example on 1D Gaussian Mixture
将目标分布设置为：$p(x)=1/3\mathcal{N}(x;-1,1)+2/3\mathcal{N}(x;2,1)$，初始分布为：$q_0(x)=\mathcal{N}(x;-10,1)$.这就产生了一个挑战，因为$p(x)和q_0(x)$的概率质量彼此相距很远（几乎没有重叠）。图1显示了我们方法中样本（n=1）的分布是如何在不同的迭代中演变的。我们发现，尽管$q_0(x)$和p（x）之间有很小的重叠，我们的方法可以将样本推向目标分布，甚至可以恢复离初始点较远的模式。我们发现了其他基于样本的算法，比如dai等人。[7]在这个toy例子中，由于$q_0(x)$的错误选择，倾向于经历degeneracy。

图2比较了我们的方法与蒙特卡罗取样法，当使用获得的样本去估计具有不同测试函数h（·）的期望$E_p(h(x))$。我们发现，我们的方法的MSE比蒙特卡罗精确采样的性能更接近或更好。这可能是因为由于repulsive force，我们的粒子比I.I.D.样品更分散，因此具有更高的估计精度。正式确定我们方法的错误率仍然是一个悬而未决的问题。

Bayesian Logistic Regression
我们考虑使用与Gershman等人[5]相同的设置对二分类问题进行贝叶斯逻辑回归。它用高斯先验$p_0(w|\alpha )=\mathcal{N}(w,\alpha -1）$和$p_0(\alpha )=gamma(\alpha ,1,0.01)$指定回归权重w。该推理被应用参数后验概率的计算$p(x|D)$，$x=[w,log\alpha ]$。我们将我们的算法与Gershman等人[5]使用的8个数据集（n>500）上的（nuts）1[29]和非参数变分推理（npv）[5]进行了比较。[5]并发现他们在这些（相对简单的）数据集上给出了非常相似的结果；更多细节见附录。

我们用581012个数据点和54个特性进一步测试二分类covertype数据集。这个数据集太大，加速需要随机梯度下降。由于nuts和npv在其代码中没有小批量选项，因此我们用welling和teh[2]等人的stochastic gradient Langevin
dynamics (SGLD)与dai[7]等人的particle mirror descent (PMD)、doubly stochastic variational inference (DSVI)、a parallel version of SGLD进行了比较。×××××。图3，一些比较结果。

Bayesian Neural Network --具体代码是什么
我们将该算法与Hern_ndez Lobato和Adams[30]在贝叶斯神经网络上提出的概率反向传播（PBP）算法进行了比较。我们的实验设置几乎是一致的，除了我们使用一个gamma（1，0.1）先验作为协方差的逆，并且不使用缩放输出层输入的技巧。我们使用具有一个隐藏层的神经网络，对大多数数据集取50个隐藏单元，除了蛋白质和年份等数据采用取100个单元；所有数据集被随机分为90%用于训练和10%用于测试，结果是20次实验的平均。除了蛋白质和年份数据，是5次和1次的平均。我们使用$relu(x）=max(0,x)$作为激活函数，其weak derivative为I[x>0]（stein特征也适用于weak derivative；见stein等人[31]）。使用作者代码6的默认设置重复PBP实验。对于我们的算法，我们只使用20个样本，并且使用动量为标准的ADAGRAD作为深度学习算法。小批量是100，除了Year数据我们使用批次大小为1000。

我们发现我们的算法在精度和速度方面都比PBP有了持续的改进；这是令人鼓舞的，因为PBP是专门为贝叶斯神经网络设计的。我们还发现，我们的结果与同一数据集（例如，32-34）上报告的最新结果具有可比性，这些数据集利用了一些我们也可以从中受益的先进技术。

6 Conclusion

提出了一种快速可扩展贝叶斯的通用变分推理算法。
推论。未来的方向包括对我们的方法有更多的理论理解，在深度学习模型中有更多的实际应用，以及第3.1节中我们基本定理的其他潜在应用。

变分推断：用q(x)近似目标分布p(x)。

方法：在集合 Q = { q ( x ) } \mathcal{Q}=\{q(x)\} Q={q(x)}搜索使$q(x)和p(x)$之间KL散度最小的$q(x)$。

集合$\mathcal{Q}$的类型 决定了 不同类型的变分推断方法。

本文的$\mathcal{Q}$：其中的q(x)由 随机变量z 组成，$z=T(x)$。T是一个光滑的一一映射，x 是从 一个容易计算的分布$q_0(x)$中采样的。

下面的问题就是围绕 T如何计算的：迭代 构造 增量变换。

---

7.公式整理

7.1 预备知识

$x$是$\mathcal{X}\subset {\mathbb{R}}^d$中的一个l连续型随机变量， { D k } \left\{D_{k}\right\} {Dk​}是 $x$独立同分布的观测集合

$x$贝叶斯推断的后验概率为：
$$p(x):=\overline{p}(x)/Z$$其中，$\overline{p}(x):=p_0(x){\prod }_{k=1}^{N}p\left(D_k|x\right)$（似然函数），$Z=\int \overline{p}(x)dx$ 为归一化常数。

令 $k\left(x,x^{\prime }\right):\mathcal{X}\times \mathcal{X}\to \mathbb{R}$是一个正定核，$k\left(x,x^{\prime }\right)$的$\mathcal{H}$空间是$f$线性张成的闭合空间。 { f : f ( x ) = ∑ i = 1 m a i k ( x , x i ) , a i ∈ R , m ∈ \left\{f : f(x)=\sum_{i=1}^{m} a_{i} k\left(x, x_{i}\right), \quad a_{i} \in \mathbb{R}, \quad m \in\right. {f:f(x)=∑i=1m​ai​k(x,xi​),ai​∈R,m∈$\mathbb{N},x_i\in \mathcal{X}\}$，$\mathcal{H}$空间内定义了内积计算：

$$⟨f,g{⟩}_{\mathcal{H}}=\sum _{ij}{a}_i{b}_{j}k\left(x_i,x_j\right)$$

其中： $g(x)={\sum }_i{b}_{i}k\left(x,x_i\right)$（内积计算只要f的系数，和g的系数和单位组成）

${\mathcal{H}}^d$表示向量函数$\mathbf{f}={\left[f_1,\dots ,f_d\right]}^{\top }$ 的空间，其中， $f_i\in \mathcal{H}$。${\mathcal{H}}^d$中定义的内积计算：
$$⟨\mathbf{f},\mathbf{g}{⟩}_{{\mathcal{H}}^d}=\sum _{i=1}^d{⟨f_i,g_i⟩}_{\mathcal{H}}$$

假定所有的向量都是列向量，令，
$${\nabla }_x\mathbf{f}=\left[{\nabla }_x{f}_1,\dots ,{\nabla }_x{f}_d\right]$$

7.2 Stein’s Identity and Kernelized Stein Discrepancy

（Stein 特征 和 核化Stein差异）
p(x) 是$\mathcal{X}\subseteq {\mathbb{R}}^d$上的一个连续可微分密度，$\varphi (x)={\left[{\varphi }_1(x),\cdots ,{\varphi }_d(x)\right]}^{\top }$光滑的向量函数。 Stein’s identity表明对于充分规则$\varphi$有如下等式成立：
$${\mathbb{E}}_{x\sim p}\left[{\mathcal{A}}_p\varphi (x)\right]=0,$$$${\mathcal{A}}_p\varphi (x)={\nabla }_x\log p(x)\varphi (x{)}^{\top }+{\nabla }_x\varphi (x)........(1)$$
其中：${\mathcal{A}}_p$为Stein operator(Stein算子)，作用于函数 $\varphi$ 使结果在$p(x)$下的期望0。通过$\varphi$上的温和边界条件进行部分积分证明stein特征的性质。（温和边界条件为：还是部分积分为）
(1)$p(x)\varphi (x)=0,\forall x\in \partial \mathcal{X}$ 当 $\mathcal{X}$是紧致的
(2)${\lim }_{r\to \infty }{\oint }_{B_r}p(x)\varphi (x{)}^{\top }\mathbf{n}(x)dS=0$ 当 $\mathcal{X}={\mathbb{R}}^d$，${\oint }_{B_r}$为以原点为中心的，半径为r的球面积分，$n(x)$为$B_r$归一化。
如果以上条件满足的话，$\varphi$作为概率函数p下的Stein类。

q(x) 是$\mathcal{X}\subseteq {\mathbb{R}}^d$上的另一个连续可微分密度，考虑stein 算子作用在$\varphi$上，且在$x\sim q$ 下的期望：
$${\mathbb{E}}_{x\sim q}[{\mathcal{A}}_p\varphi (x)]\ne 0$$

这个期望的大小衡量的$q(x)与p(x)$之间的差异，可以被定义为个差异度量–Strein差异。
合适的函数集合$\mathcal{F}$中寻找$\varphi$以最大化stein特征的违反度。（定义了一个范数用于衡量两个分布之间的差异）
D ( q , p ) = max ⁡ ϕ ∈ F { E x ∼ q [ trace ⁡ ( A p ϕ ( x ) ) ] } \mathbb{D}(q, p)=\max _{\phi \in \mathcal{F}}\left\{\mathbb{E}_{x \sim q}\left[\operatorname{trace}\left(\mathcal{A}_{p} \phi(x)\right)\right]\right\} D(q,p)=ϕ∈Fmax​{Ex∼q​[trace(Ap​ϕ(x))]}

集合$\mathcal{F}$很难选，所以采用核化Strein差异方法。
D ( q , p ) = max ⁡ ϕ ∈ H d { E x ∼ q [ trace ⁡ ( A p ϕ ( x ) ) ] ,  s.t.  ∥ ϕ ∥ H d ≤ 1 } . . . . . . ( 2 ) \mathbb{D}(q, p)=\max _{\phi \in \mathcal{H}^{d}}\left\{\mathbb{E}_{x \sim q}\left[\operatorname{trace}\left(\mathcal{A}_{p} \phi(x)\right)\right], \quad \text { s.t. } \quad\|\phi\|_{\mathcal{H}^{d}} \leq 1\right\}......(2) D(q,p)=ϕ∈Hdmax​{Ex∼q​[trace(Ap​ϕ(x))], s.t. ∥ϕ∥Hd​≤1}......(2)

将$\varphi$的搜索空间由集合$\mathcal{F}$换到了$\mathcal{H}$空间的单位球中，且说，优化问题有闭合解。
(2)最优解为[8-10]没有推导过,最优解是求出了D(q,p)，但是，D本身没有最最小化：
$$\varphi (x)={\varphi }_{q,p}^{\ast }(x)/{\Vert {\varphi }_{q,p}^{\ast }\Vert }_{{\mathcal{H}}^d}$$

$${\varphi }_{q,p}^{\ast }(\cdot )={\mathbb{E}}_{x\sim q}\left[{\mathcal{A}}_{p}k(x,\cdot )\right]->\mathbb{D}(q,p)={\Vert {\varphi }_{q,p}^{\ast }\Vert }_{{\mathcal{H}}^d}......(3)$$

$k(x,x^{\prime })$可以为RBF核：
$$k\left(x,x^{\prime }\right)=\exp \left(-\frac{1}{h}{\Vert x-x^{\prime }\Vert }_2^2\right)$$

7.3 光滑变换实现变分推断

变分推断：用简单分布$q(x)$近似目标分布$p(x)$。从预定义集合$\mathcal{Q}=q(x)$中寻找使q与pKL散度最小的q。
q ∗ = arg ⁡ min ⁡ q ∈ Q { K L ( q ∥ p ) ≡ E q [ log ⁡ q ( x ) ] − E q [ log ⁡ p ‾ ( x ) ] + log ⁡ Z } . . . . . . ( 4 ) q^{*}=\underset{q \in \mathcal{Q}}{\arg \min }\left\{\mathrm{KL}(q \| p) \equiv \mathbb{E}_{q}[\log q(x)]-\mathbb{E}_{q}[\log \overline{p}(x)]+\log Z\right\}......(4) q∗=q∈Qargmin​{KL(q∥p)≡Eq​[logq(x)]−Eq​[logp​(x)]+logZ}......(4)
集合$\mathcal{Q}$的类型决定了 不同类型的变分推断方法。本文的$\mathcal{Q}$：其中的分布由随机变量z 组成，$z=T(x)$。$T:\mathcal{X}->\mathcal{X}$是光滑的一一映射，x 是从简单分布$q_0(x)$中采样的。依据随机变量的函数变换可得：
$$q_{[T]}(z)=q\left({\mathbf{T}}^{-1}(z)\right)\cdot \left|\det \left({\nabla }_z{\mathbf{T}}^{-1}(z)\right)\right|$$

T不同,$\mathcal{Q}$集合可以近似任意的概率分布，但是为了使KL散度优化的式（4）好解，要限制T的选择。本文提出了一个新的算法用于 迭代 构造 增量 变换，能够在RKHS空间中的T上达到最速下降。

7.3.1 Stein算子作为KL散度的导数

考虑单位映射上的微小扰动变换T：
$$T(x)=x+ϵ\varphi (x)$$

其中:$\varphi (x)$是一个表征扰动方向的光滑函数，标量$ϵ$表示扰动幅度。

Stein算子与 KL散度在扰动幅度$ϵ$处的导数有莫名的联系。

定理3.1 --增量扰动T构成的 $q_{[T]}(z)$与目标分布$p(x)$之间的KL散度对$ϵ$求梯度，与，KSD定义差一个负号(证明见附录1)。
$${{\nabla }_{ϵ}\mathrm{K}\mathrm{L}\left(q_{[T]}\Vert p\right)|}_{ϵ=0}=-{\mathbb{E}}_{x\sim q}\left[\mathrm{trace}\left({\mathcal{A}}_p\varphi (x)\right)\right]......(5)$$
因为(2)式子的最优解是：$\varphi (x)=\frac{{\varphi }^{\ast }(x)}{||{\varphi }^{\ast }(x)||}$,其中 ${\varphi }_{q,p}^{\ast }(\cdot )={\mathbb{E}}_{x\sim q}\left[{\mathcal{A}}_{p}k(x,\cdot )\right]$,可以求得差异$\mathbb{D}(q,p)=||{\varphi }^{\ast }(x)||$（右边搜索$\varphi (x)$求最大=>左边的梯度就最大）上式左边求最大，就是右边求最大，就是KL散度的梯度的值最大，对应方向：变化最快的方向，KL散度下降最大的方向。此时最优扰动方向${\varphi }^{\ast }(x)$,它可以使$q_{[T]}快速靠近p$：
$${\varphi }_{q,p}^{\ast }(\cdot )={\mathbb{E}}_{x\sim q}\left[{\mathcal{A}}_{p}k(x,\cdot )\right]......(6)$$

引理3.2 带入最优解时(5)式
$${{\nabla }_{ϵ}\mathrm{K}\mathrm{L}\left(q_{[T]}\Vert p\right)|}_{ϵ=0}=-{\mathbb{D}}^2(q,p)$$
引理3.2展示了将 一个初始分布$q_0$ 通过 迭代得到 目标分布$p$ 的迭代式(${\varphi }^{\ast }(x)$包含了目标分布$p$的信息):

${\mathbf{T}}_0^{\ast }(x)=x+{ϵ}_0\cdot {\mathbf{\varphi }}_{q_0,p}^{\ast }(x)$ where $q_1(x)=q_{0\left[{\mathbf{T}}_0\right]}(x)$,

重复这个过程直至T为单位阵：
$$q_{\ell +1}=q_{\ell |T_{\ell }|},\text{ where }{T}_{\ell }^{\ast }(x)=x+{ϵ}_{\ell }\cdot {\varphi }_{q_{\ell ,p}}^{\ast }(x)......(7)$$
函数梯度 为了进一步了解这个过程，我们现将重新解释(6)式子为RKHS中的函数梯度。对于$f\in {\mathcal{H}}^d$任意方程$F[f]$，其梯度为${▽}_{f}F[f]$，函数F的傅里叶展开为：
$$F[f+ϵg(x)=F[f]+ϵ<{▽}_{f}F[f],g>+O({ϵ}^2)$$

对于任意的$g\in {\mathcal{H}}^d$和$ϵ\in \mathcal{R}$

定理3.3–让$\mathbf{T}(x)=x+\mathbf{f}(x),$ 其中 $\mathbf{f}\in {\mathcal{H}}^d$，$q_{[T]}$是$z=T(x)$的密度，$x\sim q$
$${{\nabla }_f\mathrm{KL}\left(q_{[T]}\Vert p\right)|}_{f=0}=-{\varphi }_{q,p}^{\ast }(x)$$
${\Vert {\varphi }_{q,p}^{\ast }\Vert }_{{\mathcal{H}}^d}=\mathbb{D}(q,p)$是RKHS范数。
上式的好处：只算了$ϵ=0$出的梯度，在$ϵ\ne 0$处，的梯度要计算$T^{\prime }(x{)}^{-1}$,这是十分难以计算的。

7.3.2 Strein 变分梯度下降

附录1：

$$\begin{gathered} KL(q_{[T]}||p)=E_{q_{[T]}}[log(q_{[T]}(z))]-E_{q_{[T]}}[log(p(z))] \\ =\int q_{[T]}(z)log(q_{[T]}(z))dz-\int q_{[T]}(z)log(p(z))dz......(1) \end{gathered}$$
随机变量函数变换$q_{[T]}(z)=q(T^{-1}(z))T^{-1}(z{)}^{\prime }$,$T^{-1}(z{)}^{\prime }=T^{\prime }(x{)}^{-1}$：
$$\begin{gathered} \int q_{[T]}(z)log(q_{[T]}(z))dz=\int q(T^{-1}(z))T^{-1}(z{)}^{\prime }log(q(T^{-1}(z))T^{-1}(z{)}^{\prime })dz \\ =\int q(T^{-1}(z))log(q(T^{-1}(z))T^{-1}(z{)}^{\prime })d{T}^{-1}(z) \\ =\int q(T^{-1}(z))log(q(T^{-1}(z))d{T}^{-1}(z)+\int q(T^{-1}(z))log(T^{-1}(z{)}^{\prime }))d{T}^{-1}(z) \\ =\int q(x)logq(x)dx+\int q(x)log(T^{\prime }(x{)}^{-1})dx......(2) \end{gathered}$$

$$\begin{gathered} -\int q_{[T]}(z)log(p(z))dz=-\int q(T^{-1}(z))T^{-1}(z{)}^{\prime }log(p(z))dz \\ =-\int q(T^{-1}(z))log(p(z))d{T}^{-1}(z)=-\int q(x)log(p(x+ϵ\varphi (x)))dx......(3) \end{gathered}$$
所以：
$$KL(q_{[T]}||p)=(2)+(3)=\int q(x)logq(x)dx+\int q(x)log(T^{\prime }(x{)}^{-1})dx-\int q(x)log(p(x+ϵ\varphi (x)))dx......(4)$$
(4)对$ϵ$求导(只有后两项与$ϵ$有关$T^{\prime }(x)=1+ϵ{\varphi }^{\prime }(x)$)
$${▽}_{ϵ}KL(q_{[T]}||p)=-\int q(x)\frac{1}{T^{\prime }(x)}\varphi (x)dx-\int q(x)\frac{1}{p(x+ϵ\varphi (x))}{p}^{\prime }(x+ϵ\varphi (x)){\varphi }^{\prime }(x)dx....(5)$$
当 $ϵ=0$ (5)可以化为：
$${▽}_{ϵ}KL(q_{[T]}||p){|}_{ϵ=0}=-\int q(x){\varphi }^{\prime }(x)dx-\int q(x)\frac{1}{p(x)}{p}^{\prime }(x)\varphi (x)dx....(5)$$
$x$为一维情况下：
$${▽}_{ϵ}KL(q_{[T]}||p){|}_{ϵ=0}=-E_{x\sim q}[{▽}_{x}logp(x)\varphi (x)+{▽}_x\varphi (x)]$$
$x$不为1维的情况：
$${▽}_{ϵ}KL(q_{[T]}||p){|}_{ϵ=0}=-E_{x\sim q}[trace({▽}_{x}logp(x)\varphi (x)+{▽}_x\varphi (x))]=-E_{x\sim q}[trace(A_p\varphi (x))]$$