论文地址：点这里。作者还提供了Stein变分梯度下降法的源码。
Note：
源码不涉及深度学习，所以PyTorch用户或者TF用户都可以使用。

Stein变分梯度下降(SVGD)可以理解是一种和随机梯度下降(SGD)一样的优化算法。在强化学习算法中，Soft-Q-Learning使用了SVGD去优化，而Soft-AC选择了SGD去做优化。

Stein Variational Gradient Descent:A General Purpose Bayesian Inference Algorithm

Abstract
1 Introduction
2 Background
3 Variational Inference Using Smooth Transforms
- 3.1 Stein Operator as the Derivative of KL Divergence
- 3.2 Stein Variational Gradient Descent
4 Related Works
5 Experiments
6 Conclusion
源码解析

Abstract

本文旨在提供一种新的变分推断(VI)算法用来做梯度下降(优化)。通过优化两个分布的KL散度，粒子不断迭代，使之从初始的分布一步步达到目标分布。算法涉及到Stein特征的平滑转换以及核差异(KSD)。
Note：

关于KSD的介绍，可以参考我的另一篇论文笔记。

1 Introduction

贝叶斯推断是一种模型估计的方法，核心就是那个贝叶斯公式。该公式中有一个归一化因子p(x)=∫p(x)dxp(x)=\int p(x)\mathrm{d}xp(x)=∫p(x)dx，当高维问题时候，计算将变得困难。因此MCMC和VI这两种方法就应运而生，用于解决这个计算复杂度问题。
马尔科夫链蒙特卡洛(MCMC)的缺陷在于收敛慢且难收敛，优点是偏差较小，毕竟是个需要采样的方法(强化学习最初的MC方法也有这2个特点)。
变分推断(VI)是一种优化算法，用随机梯度下降最小化目标分布与估计分布之间的KL散度，其优点是优化速度快适合于大数据，缺陷是存在较大的偏差，优化结果的好坏与否取决于最先定义的数簇(如经典的平均场变分数簇)偏差大小。
最大后验概率(MAP)是一种依靠先验知识的计算后验概率的方法，这是一种MAP贝叶斯，并不是全贝叶斯。
作者提出一种新型通用的变分推断算法——Stein变分梯度下降。可以看成是一种用于贝叶斯推断的梯度下降方法。它有以下特点：

使用粒子去做分布的估计。
是一种梯度下降算法，可以用在任何使用SGD做优化的地方。
是一种最小化目标分布和估计分布间KL散度的算法，不断地驱动粒子去拟合(fit)目标分布。
算法引入了平滑转换、KSD(Stein方法+RKHS)。
是一种全贝叶斯算法。
当使用单个粒子的时候，该算法就相当于对MAP做梯度下降；当使用所有的粒子的时候，就是全贝叶斯算法。
是一种超越经典VI的新VI算法。
Stein梯度下降最大的贡献在于：它在再生核希尔伯特空间下给出了使得KL散度下降最快的确定性方向。

2 Background

这里是对贝叶斯中的后验分布p(x)=pˉ(x)/Zp(x)=\bar{p}(x)/Zp(x)=pˉ(x)/Z以及KSD的介绍。KSD有个优点在于通过使用Stein得分函数∇xlog⁡p(x)\nabla_x\log p(x)∇xlogp(x)来避免计算归一化因子ZZZ。

3 Variational Inference Using Smooth Transforms

VI做模型估计：
目标分布p(x)p(x)p(x)，估计分布q(x)q(x)q(x)，数簇Q={q(x)}\mathcal{Q}=\{q(x)\}Q={q(x)}，通过优化KL散度来得到估计模型q∗q^*q∗：
q∗=arg min⁡q∈Q{KL(q∣∣p)=Eq[log⁡q(x)]−Eq[log⁡pˉ(x)]+log⁡Z}(4)q^*=\argmin_{q\in\mathcal{Q}}\{KL(q||p)=\mathbb{E}_q[\log q(x)]-\mathbb{E}_q[\log\bar{p}(x)]+\log Z\}\tag{4} q∗=q∈Qargmin{KL(q∣∣p)=Eq[logq(x)]−Eq[logpˉ(x)]+logZ}(4)
VI就是通过这种方式来免除归一化因子的计算。问题就是这个数簇Q\mathcal{Q}Q很难去选择一个合适通用的。
作者指出，合适的Q\mathcal{Q}Q应该具备以下原则：

准确性：广度足够以至于可以近似一系列目标分布。
计算可行性：容易去计算，不要整个难以实现的。
可解决性：可以和KL散度优化结合起来。

在此基础上，作者引入平滑变换T:X→X,X⊂RdT:\mathcal{X}\to\mathcal{X},\mathcal{X}\subset\mathbb{R}^dT:X→X,X⊂Rd，也就是说将Q\mathcal{Q}Q看成一个一对一平滑变换z=T(x),x∈Xz=T(x),x\in\mathcal{X}z=T(x),x∈X得到的数簇，且x∼q0(x)x\sim q_0(x)x∼q0(x)。
Note：

一对一的原因在于：①必须保证是函数。②如果是多对一的话，就不能保证是增量转换(单调)了，比如如果T是个y=x2y=x^2y=x2，那么当你粒子xxx经过好几轮平滑转化之后，qqq已经发生改变了，这时候你从qqq采样的xxx经过平滑转换可能就会回到之前发生过的某个qqq分布了，这是不允许的，因为这样就有点回环曲折的意思，不能保证最快的优化了。

那么经过TTT之后的数簇是怎么样的呢？利用坐标变换公式可得：
q[T](z)=q(x)⋅∣det⁡(∂x∂z)∣=q(T−1(z))⋅∣det⁡(∇zT−1(z))∣q_{[T]}(z)=q(x)\cdot|\det(\frac{\partial x}{\partial z})|=q(T^{-1}(z))\cdot|\det(\nabla_zT^{-1}(z))| q[T](z)=q(x)⋅∣det(∂z∂x)∣=q(T−1(z))⋅∣det(∇zT−1(z))∣其中det⁡(∇zT−1(z))\det(\nabla_zT^{-1}(z))det(∇zT−1(z))是矩阵T−1(z)T^{-1}(z)T−1(z)的雅克比行列式。
或者反过来：
q[T−1](x)=q(T(x))⋅∣det⁡(∇xT(x))∣q_{[T^{-1}]}(x) = q(T(x))\cdot|\det(\nabla_xT(x))| q[T−1](x)=q(T(x))⋅∣det(∇xT(x))∣作者指出平滑变化的引入使得准确性和计算可行性满足了，可是可解决性还不能满足，因此一种方法是将TTT参数化，但参数化将面临3个问题：

需要选择合适的参数表示来平衡三个原则。
参数化表示不一定能满足T的一对一性。
参数化不一定使得Jacobian可以得到有效计算。

参数化既然这么麻烦，作者索性放弃选用了一种不需要参数化表示TTT的方法——一种迭代式构建增量变化的方法，在RKHS下，可表现出最快的梯度下降方向：

该方法不需要计算雅克比矩阵。
形式上简单，类似于经典的梯度下降算法，可计算性强。

3.1 Stein Operator as the Derivative of KL Divergence

作者提出了转换的表达式：T(x)=x+ϵ⋅ϕ(x)T(x)=x+\epsilon\cdot\phi(x)T(x)=x+ϵ⋅ϕ(x)。
Note：

这个表达式长得很像梯度下降吧，意味着ϕ(x)\phi(x)ϕ(x)代表着方向。
ϵ\epsilonϵ是一个很小很小的数。
ϕ(x)\phi(x)ϕ(x)可微，是一个向量。
当∣ϵ∣|\epsilon|∣ϵ∣很小的时候。T(x)≈xT(x)\approx xT(x)≈x，所以其雅可比矩阵∇xT−1\nabla_xT^{-1}∇xT−1就是个单位阵，显然雅克比行列式不等于0，根据反函数理论，保证了TTT的单调性，也就保证了一对一。

接下来就是本文最重要的地方之一了。

Lemma1
如果X∈X,X⊂RX\in\mathcal{X},\mathcal{X}\subset\mathbb{R}X∈X,X⊂R，FFF是平滑变换，则：
∂det⁡(F(X))∂F(X)T=det⁡(F(X))⋅F−1(X)\frac{\partial{\det(F(X))}}{\partial{F(X)}^T}=\det(F(X))\cdot F^{-1}(X) ∂F(X)T∂det(F(X))=det(F(X))⋅F−1(X)证明如下：

Note:

求导对象xxx必须是方阵。

Lemma2
如果p(x)、q(x)p(x)、q(x)p(x)、q(x)都是光滑的函数，x∈Xx\in\mathcal{X}x∈X，T=Tϵ(x)T=T_\epsilon(x)T=Tϵ(x)是个一对一的转换函数，且对x,ϵx,\epsilonx,ϵ都可微，sp(x)s_p(x)sp(x)为Stein得分函数，q[T]q_{[T]}q[T]是x∼qx\sim qx∼q时z=Tϵ(x)z=T_\epsilon(x)z=Tϵ(x)的概率密度函数，则：
∇ϵKL(q[T]∣∣p)=Eq[spT(T(x))∇ϵT(x)+tr((∇xT(x))−1⋅(∇ϵ∇xT(x)))]\nabla_\epsilon KL(q[T]||p)=\mathbb{E}_q[s_p^T(T(x))\nabla_\epsilon T(x)+ tr((\nabla_xT(x))^{-1}\cdot(\nabla_\epsilon\nabla_xT(x)))] ∇ϵKL(q[T]∣∣p)=Eq[spT(T(x))∇ϵT(x)+tr((∇xT(x))−1⋅(∇ϵ∇xT(x)))]证明如下：

Note：

画红色的转置TTT和迹trtrtr是为了输出格式的需要。
原论文补充材料中的推导对最后的结果漏了一个转置符TTT。

Theorem 3.1.
T(x)=x+ϵ⋅ϕ(x)T(x)=x+\epsilon\cdot\phi(x)T(x)=x+ϵ⋅ϕ(x)，平滑变换z=T(x)z=T(x)z=T(x)以及未知分布q[T](z)q_{[T]}(z)q[T](z)，当x∼p(x)x\sim p(x)x∼p(x)时，有
∇ϵKL(q[T]∣∣p)∣ϵ=0=−Ex∼q[tr(Apϕ(x))](5)\nabla_\epsilon KL(q_{[T]}||p)|_{\epsilon=0}=-\mathbb{E}_{x\sim q}[tr(\mathcal{A}_p\phi(x))]\tag{5} ∇ϵKL(q[T]∣∣p)∣ϵ=0=−Ex∼q[tr(Apϕ(x))](5)证明如下：

Note:

Apϕ(x)=spϕ(x)T+∇xϕ(x)\mathcal{A}_p\phi(x)=s_p\phi(x)^T+\nabla_x\phi(x)Apϕ(x)=spϕ(x)T+∇xϕ(x)
这里和KSD的原论文相反，这里是qqq未知，为估计分布，ppp已知，为目标分布。
等式左边涉及到泛函导数、变分法的知识，变分法是专门用于处理泛函问题的工具，相当于微积分法是专门用于处理函数问题的工具，我们要找到使得KL散度下降最快的函数TTT，这就是典型的泛函问题，我们求的是泛函极值问题，所以该问题就是变分问题。
等式右边是个和KSD相关的东西，等式左边代表着KL散度的方向导数。这个怎么理解呢？根据泛函分析，左边的表达式如果等于0的话，意味着此时KL散度取得极值，但由于ϵ=0\epsilon=0ϵ=0之后仍然会有变量xxx的存在，所以这一次你右边不会等于一个0，那ϵ=0\epsilon=0ϵ=0出来的是什么呢？

如上图所示，可以把KL散度理解为一个二元的函数F(x,ϵ)F(x,\epsilon)F(x,ϵ)，实际上的图应该是个立体的，有许多个小山坡，这只是简易表示。分析：式(5)左边这是一个对ϵ\epsilonϵ的偏导数，即只关于ϵ\epsilonϵ的斜率。式(5)表达的就是点P的情况，也就是说KL散度关于ϵ\epsilonϵ的斜率是个和变量xxx有关的值，但可以是这条曲线上的P1、P2、P3点的斜率值，具体哪个取决于xxx的值，但如上面所说，实际上是立体图，也就是说确定了xxx，即确定了哪个P点之后，比如确定了P1，那么立体图上这个斜率是有无数多种的，沿着不同方向有不同的变化率，这其实就是方向导数！所以式(5)就是方向导数，一看到方向导数，就想到了梯度，它是方向导数的最大值，也是变化率最大的方向，巧的是式(5)右边就是个能进行优化求出最大值的表达式，当ϕ(x)\phi(x)ϕ(x)取特定的值，右边就能达到最大。至于这个“负号”，并不碍事，我们在利用式(5)的时候，只是表示从上升最快转成下降最快而已，这样便于KL梯度下降，梯度的反方向就是下降最快的方向。

Lemma 3.2.
如果所有的方向ϕ(x)\phi(x)ϕ(x)都在球内：B={ϕ∈Hd:∣∣ϕ∣∣Hd2≤S(q,p)}\mathcal{B}=\{\phi\in\mathcal{H}^d:||\phi||^2_{\mathcal{H}^d}\leq\mathbb{S}(q,p)\}B={ϕ∈Hd:∣∣ϕ∣∣Hd2≤S(q,p)}，则在这些方向中能最大化式(5)中那个负梯度的方向为：
ϕ∗(⋅)=Ex∼q[sp(x)k(x,⋅)+∇xk(x,⋅)]=Ex∼q[Apk(x,⋅)](6)\phi^*(\cdot)=\mathbb{E}_{x\sim q}[s_p(x)k(x,\cdot)+\nabla_xk(x,\cdot)]=\mathbb{E}_{x\sim q}[\mathcal{A}_pk(x,\cdot)]\tag{6} ϕ∗(⋅)=Ex∼q[sp(x)k(x,⋅)+∇xk(x,⋅)]=Ex∼q[Apk(x,⋅)](6)Note:

ϕ=ϕ(⋅),k(x,⋅)\phi=\phi(\cdot),k(x,\cdot)ϕ=ϕ(⋅),k(x,⋅)都表示向量。
sps_psp为Stein得分函数。
这个引理来自于KSD那篇论文的Theorem3.8。
式(6)和论文中的不太一样，从我的理解角度看，论文中表达的意思应该就是我这里式(6)写的那样。
当ϕ=ϕ∗\phi=\phi^*ϕ=ϕ∗的时候，式(5)就相当于：∇ϵKL(q[T]∣∣p)∣ϵ=0=−S(p,q)\nabla_\epsilon KL(q_{[T]}||p)|_{\epsilon=0}=-\mathbb{S}(p,q)∇ϵKL(q[T]∣∣p)∣ϵ=0=−S(p,q)。换句话说当ϕ=ϕ∗\phi=\phi^*ϕ=ϕ∗的时候，是KL散度下降最快的方向。

当ϕ\phiϕ取最优值ϕ∗\phi^*ϕ∗的时候，只要变量沿着方向ϕ∗\phi^*ϕ∗移动ϵ⋅S(q,p)\epsilon\cdot\mathbb{S(q,p)}ϵ⋅S(q,p)，就可以使得KL散度以最快的速度下降(从这里也看出了ϕ(x)\phi(x)ϕ(x)代表着方向)。同时，平滑转换T=x+ϵ⋅ϕ∗(x)T=x+\epsilon\cdot\phi^*(x)T=x+ϵ⋅ϕ∗(x)。总的流程用一个动态过程表示为：
起始分布q0q_0q0；
第一次梯度下降：T0∗(x)=x+ϵ0⋅ϕq0,p∗(x);q1(x)=q0[T0∗](x)T^*_0(x)=x+\epsilon_0\cdot\phi^*_{q_0,p}(x);q_1(x)=q_{0[T_0^*]}(x)T0∗(x)=x+ϵ0⋅ϕq0,p∗(x);q1(x)=q0[T0∗](x)；
第二次梯度下降：T1∗(x)=x+ϵ1⋅ϕq1,p∗(x);q2(x)=q1[T1∗](x)T^*_1(x)=x+\epsilon_1\cdot\phi^*_{q_1,p}(x);q_2(x)=q_{1[T_1^*]}(x)T1∗(x)=x+ϵ1⋅ϕq1,p∗(x);q2(x)=q1[T1∗](x)；
… … …直至收敛
总结起来就是：
ql+1=ql[Tl∗]，Tl∗(x)=x+ϵl⋅ϕql,p∗(x)(7)q_{l+1}=q_{l[T^*_l]}，T^*_l(x)=x+\epsilon_l\cdot\phi^*_{q_l,p}(x)\tag{7} ql+1=ql[Tl∗]，Tl∗(x)=x+ϵl⋅ϕql,p∗(x)(7)Note:

T,qT,qT,q的不断迭代也显示了平滑转换一对一性的重要，因此ϵ\epsilonϵ要足够小。
最终KL散度会因为不断减小而为0，qqq会收敛到目标分布ppp，当且仅当ϕp,q∞∗(x)=0\phi^*_{p,q^\infty}(x)=0ϕp,q∞∗(x)=0

Functional Gradient(泛函梯度)
设f、g∈Hd\mathbf{f}、\mathbf{g}\in\mathcal{H}^df、g∈Hd，则F[f]F[f]F[f]就是泛函，∇fF[f]\nabla_\mathbf{f}F[\mathbf{f}]∇fF[f]就是泛函梯度，其有如下结论：F[f+ϵg(x)]=F[f]+ϵ<∇fF[f],g>Hd+O(ϵ2)F[\mathbf{f}+\epsilon\mathbf{g}(x)]=F[f]+\epsilon<\nabla_\mathbf{f}F[\mathbf{f}],\mathbf{g}>_{\mathcal{H}^d}+O(\epsilon^2)F[f+ϵg(x)]=F[f]+ϵ<∇fF[f],g>Hd+O(ϵ2)。
Theorem 3.3.
设T(x)=x+f(x),x∈Rd,f∈HdT(x)=x+\mathbf{f}(x),x\in\mathbb{R}^d,\mathbf{f}\in\mathcal{H}^dT(x)=x+f(x),x∈Rd,f∈Hd，则：
∇fKL(q[T]∣∣p)∣f=0=−ϕq,p∗(x)\nabla_{\mathbf{f}}KL(q_{[T]}||p)|_{\mathbf{f}=0}=-\phi^*_{q,p}(x) ∇fKL(q[T]∣∣p)∣f=0=−ϕq,p∗(x)Note:

这个式子我最终没有证明出来，看过作者在补充材料的证明是有问题的，要么就是矩阵向量前后大小对不上，要么就是最后部分漏东西，最后f=0\mathbf{f}=0f=0带进去也是得不出ϕ∗=Ex∼q(x)[Apk(x,⋅)]\phi^*=\mathbb{E}_{x\sim q(x)}[\mathcal{A}_pk(x,\cdot)]ϕ∗=Ex∼q(x)[Apk(x,⋅)]这个结果的，只能证明到Ex∼q(x)[tr(Apk(x,⋅))]\mathbb{E}_{x\sim q(x)}[tr(\mathcal{A}_pk(x,\cdot))]Ex∼q(x)[tr(Apk(x,⋅))]。
以下是我的证明：接下来需要用到向量矩阵的Taylor展开(和高数中的略有不同)：
总之感觉作者在凑这个ϕ∗\phi^*ϕ∗，这个Theorem 3.3应该算错的，因此下面关于这个定理的就不分析了。不过好在这个定理的好坏并不影响Algorithm1，即Stein变分梯度下降。
另外需要注意的是，我们一般向量内积或者矩阵内积的结果都是实值，默认都是大小一样的且向量必须同为行或者列向量，矩阵内积必须均为方阵。如果<A,B><A,B><A,B>运算的A和B(向量和向量、向量和矩阵)大小不一样，则输出的结果就不一定是实值了，可能是向量或者矩阵。

3.2 Stein Variational Gradient Descent

算法的伪代码如下：

伪码简洁明了，左边是公式(7)，不断更新粒子(可以理解成样本)来改变初始分布q0q_0q0。
粒子{xil}i=1n\{x_i^l\}^n_{i=1}{xil}i=1n表示第lll次迭代的第iii个粒子，一个nnn个。粒子最开始是从分布q0q_0q0中采样的，最初的分布qqq是不准确的并且可以是任意的。也就是说，该算法不依赖于初始的分布。
准确的说，粒子的统计分布就是q(x)q(x)q(x)。比如第lll次迭代，统计粒子的均值、方差构成高斯分布p(x)p(x)p(x)。伪代码所做的就是用这些粒子去拟合目标分布p(x)p(x)p(x)。
右边是公式(6)的近似表示，用样本均值作为期望的估计值，这是常用的做法了。
算法模仿了我们所熟知的梯度上升，之前说过，ϕ\phiϕ代表了梯度方向，为KL散度下降最快的方向。该算法的梯度上升模式是粒子级别的，我们之前接触比较多的神经网络，梯度上升都是参数级别的。
公式(8)的Part1和Part2类似于与探索与利用这种相互排斥的操作。首先Part1的∇\nabla∇部分，粒子会朝着ppp分布概率高的地方移动，但这样容易陷入局部最优，不要以为SVGD就不会陷入局部最优了，SVGD确实保证了每一步都向着KL减小最快的方向移动，这当然包括局部最小和全局最小，所以是否陷入局部最优和算法优化速度快慢是两回事。至于Part2的∇\nabla∇部分，作者举了个RBF核的例子：k(x,x′)=exp⁡(−1h∣∣x−x′∣∣2)k(x,x')=\exp(-\frac{1}{h}||x-x'||^2)k(x,x′)=exp(−h1∣∣x−x′∣∣2)，经过Part2的处理后：∑j2h(x−xj)k(xj,x)\sum_j\frac{2}{h}(x-x_j)k(x_j, x)∑jh2(x−xj)k(xj,x)，这一项代表着粒子将会朝着远离当前迭代轮数lll的粒子，从而减轻局部最优的风险。同样是这个高斯核，当h→0h\to0h→0，那么这个排斥作用就没了，只剩下了Part1，那么粒子都会朝着最大化log⁡p(x)\log p(x)logp(x)移动，这也就是MAP，从而陷入局部最优。
当n=1的时候，算法就演变成了对MAP做梯度上升。
算法中给出了确定性的下降方法，远远的优于传统的MCMC算法MCMC算法中使用随机性的方法来增加粒子的多样性，在Stein 变分中根本就不需要。

提高计算效率的方法
核心思想：下采样。
下采样在信号中是对于一个样值序列隔几个样值采样一次。对图像来说，下采样就是对于一副图像I尺寸为M×NM\times NM×N，对起进行s倍下采样，即得到(M/s)×(N/s)(M/s)\times(N/s)(M/s)×(N/s)尺寸的分辨率图像。通俗来说就是采样少部分的样本来近似达到大样本的效果，并且减小计算量。
公式(8)计算的瓶颈在于∇xlog⁡p(x),x∈{xi}i=1n\nabla_x\log p(x),x\in\{x_i\}^n_{i=1}∇xlogp(x),x∈{xi}i=1n。其中p(x)∝p0(x)∏k=1Np(Dk∣x)p(x)\propto p_0(x)\prod^N_{k=1}p(D_k|x)p(x)∝p0(x)∏k=1Np(Dk∣x)，当数据量NNN很大的时候，计算复杂度就会很高。作者提出用下采样Ω⊂{1,⋯,N}\Omega\subset\{1,\cdots,N\}Ω⊂{1,⋯,N}来近似估计∇xlog⁡p(x)\nabla_x\log p(x)∇xlogp(x)。本文中的下采样是这样的：
log⁡p(x)≈log⁡p0(x)+N∣Ω∣∑k∈Ωlog⁡p(Dk∣x).(9)\log p(x)\approx\log p_0(x)+\frac{N}{|\Omega|}\sum_{k\in\Omega}\log p(D_k|x).\tag{9} logp(x)≈logp0(x)+∣Ω∣Nk∈Ω∑logp(Dk∣x).(9)设∇xlog⁡p(x)\nabla_x\log p(x)∇xlogp(x)在全粒子上的统计结果为RRR，在下采样的表现为GGG，它的原理就是：
NΩ=RG\frac{N}{\Omega}=\frac{R}{G} ΩN=GR用局部代替整体。
Note：

遇到大数据NNN且需要计算∑i=1N\sum_{i=1}^N∑i=1N的时候，可采用下采样减小计算复杂度。

4 Related Works

略

5 Experiments

略

6 Conclusion

总结：

Stein变分梯度下降是这样的：首先我们需要用粒子拟合的估计分布q(x)q(x)q(x)通过不断迭代靠近目标分布p(x)p(x)p(x)，每一次迭代都是一次平滑变化TTT，平滑变化有很多种，我们要找到使得KL散度下降最快的变化方向(梯度)。因此我们的计算过程其实就优化函数中的出一个最优函数，这就是泛函问题中求泛函极值问题，我们通常使用变分法求解。
本文提出的Stein变分梯度下降是结合了Stein方法，并在RKHS空间下进行梯度下降的(优化)，即使用了KSD核差异技术。
Stein变分梯度下降算法使用了粒子层次化的方式不断使2个分布接近一致。从上图2种优化算法的参数轨迹曲线可以看出。
SVGD显然比我们最先接触到的机器学习经典优化算法——SGD具有更快的收敛速度。SVGD是确定性的优化算法，而SGD是随机性的优化算法，其优化方向取决于样本，有随机性，因此稳定性会打折扣，同时收敛速度收到削弱。
SVGD的S指的是Stein；SGD的S指的是Stochastic。
SVGD是一种变分推断VI的新方法，可广泛用于快速贝叶斯推断。

Stein变分梯度下降的优缺点：
优点：

Stein方法通过Stein得分函数有效化解掉归一化因子。
Stein特征可以消除未知分布无法计算的问题。
引入KSD，推导出使得KL散度下降最快的方向就是核差异方向ϕ∗\phi^*ϕ∗，这是一种确定性方向，有效增加优化的速度和效率(分析见上)。不像传统的随机性方法，需要通过随机性来增加粒子多样性。
梯度方向ϕ∗\phi^*ϕ∗的近似包含正则化，避免陷入局部最优。
算法不依赖于初始化粒子是怎么样的，无论如何都会走向复杂的分布，一般都是接近全局最优的局部最优。

缺点：

计算复杂度高。

源码解析

关于裸的SVGD代码地址在本文顶部
main.py如下：

A = np.array([[0.2260,0.1652],[0.1652,0.6779]])
mu = np.array([-0.6871,0.8010])
model = MVN(mu, A)
x0 = np.random.normal(0,1, [10,2]);
theta = SVGD().update(x0, model.dlnprob, n_iter=1000, stepsize=0.01)
print("ground truth: ", mu)
print("svgd: ", np.mean(theta,axis=0))

A是目标分布的协方差，mu是目标分布的均值，model是一个二维高斯分布。
x0是初始分布q0q_0q0为一维正态分布采样的10个粒子。将10个粒子送进SVGD的核心代码(Algorithm1)中，输出经过1000次迭代后的的10个粒子分布情况，ϵ=0.01\epsilon=0.01ϵ=0.01，他们的样本均值就是估计分布的均值。
三次试验结果如下：
①
ground truth: [-0.6871 0.801 ]
svgd: [-0.68622667 0.80359776]
②
ground truth: [-0.6871 0.801 ]
svgd: [-0.68865996 0.7990209 ]
③
ground truth: [-0.6871 0.801 ]
svgd: [-0.6861909 0.80438672]

可见SVGD的有效性。

论文笔记之Stein变分梯度下降相关推荐

吴恩达机器学习MATLAB代码笔记（1）梯度下降
吴恩达机器学习MATLAB代码笔记(1)梯度下降单变量线性回归 1.标记数据点(Plotting the Date) fprintf('Plotting Data') data = load('D: ...
最佳论文！牛津大学揭示梯度下降复杂度理论
点上方计算机视觉联盟获取更多干货仅作学术分享,不代表本公众号立场,侵权联系删除转载于:机器之心 AI博士笔记系列推荐周志华<机器学习>手推笔记正式开源!可打印版本附pdf下载链接梯 ...
机器学习/算法面试笔记1——损失函数、梯度下降、优化算法、过拟合和欠拟合、正则化与稀疏性、归一化、激活函数
正值秋招,参考网络资源整理了一些面试笔记,第一篇包括以下7部分. 1.损失函数 2.梯度下降 3.优化算法 4.过拟合和欠拟合 5.正则化与稀疏性 6.归一化 7.激活函数损失函数损失函数分为经验 ...
机器学习笔记：线性规划，梯度下降
主要内容来自stanford Andrew Ng视频课程的总结. 讲的非常好,还有相应的习题,课程能够在以下站点搜索到. https://www.coursera.org/ 机器学习的目的是在训练数据 ...
论文笔记之Soft Q-learning
论文地址,点这里源码地址,点这里参考: soft Q-learning-v1 soft Q-learning-v2 soft Q-learning-v3 本人阅读目的: 这篇文章是Soft Act ...
LeCun称梯度下降是最优雅的 ML 算法，Marcus：我不同意
作者|李梅编辑|陈彩娴来源|AI科技评论前几天刚跟马斯克吵完架的Gary Marcus,又双叒叕跟人吵起来了,这次的吵架对象是Yann LeCun. 一向喜欢给深度学习泼冷水的Marcus,在今 ...
LeCun称梯度下降是最优雅的 ML 算法，Marcus：我不同意！
点击下方卡片,关注"CVer"公众号 AI/CV重磅干货,第一时间送达点击进入-> CVer 微信技术交流群转载自:AI科技评论作者|李梅编辑|陈彩娴前几天刚跟 ...
《深度学习》李宏毅 -- task3误差和梯度下降
误差和梯度下降 1. 误差的来源 1.1 error (误差) 1.2 bias(偏差) 1.3 variance (方差) 1.4 举例说明 1.5 偏差方差困境 2 估测 2.1 估测X的偏差 2 ...
论文笔记 Inference in Deep Gaussian Processes using Stochastic Gradient Hamiltonia使用随机梯度哈密顿量蒙特卡罗推理深度高斯过程
使用随机梯度哈密顿量蒙特卡罗推理深度高斯过程 0.摘要深度高斯过程 (DGP) 是高斯过程的层次概括,它将经过良好校准的不确定性估计与多层模型的高度灵活性相结合. 这些模型的最大挑战之一是精确推断是 ...

论文笔记之Stein变分梯度下降