Stein variational gradient descent（SVGD）
最近读了Liu Qiang组关于SVGD的几篇文章，和老师讨论了很久、查阅了很多文献后大概了解了它的思路。这篇文章整合一下我所理解的SVGD。由于太笨，不知道Katex怎么像Latex一样写公式标号····所以本文所有公式都没有标号Orz
近似推断被广泛用于概率机器学习与统计中，Stein variational gradient descent （SVGD）是由Liu Qiang等提出的一种近似推断算法。不同于MCMC，它是一种确定性的算法。不同于变分推断（VI），它采用粒子方法直接对目标概率分布进行逼近。
1.Stein’s Identity
首先Liu Qiang组关于SVGD的idea出发点都是Stein’s Identity。Stein’s Identity讲了这么一件事情：
E p [ ∇ x log q ( x ) f ( x ) + ∇ x f ( x ) ] = 0 i f p ( x ) = q ( x ) \mathbb{E}_{p}\left[\nabla_{x} \log q(x) f(x)+\nabla_{x} f(x)\right]=0\ \ \ if\ \ \ p(x)=q(x) Ep[∇xlogq(x)f(x)+∇xf(x)]=0 if p(x)=q(x)
我们考虑 R d \mathbb{R}^{d} Rd上两个光滑分布 p p p和 q q q，上式成立的充分条件是 p ( x ) = q ( x ) p(x)=q(x) p(x)=q(x). 这一等式的成立非常好证，只需要将积分展开用分部积分法就可以证明( ∫ u ( x ) v ′ ( x ) d x = u ( x ) v ( x ) − ∫ u ′ ( x ) v ( x ) d x \int u(x) v^{\prime}(x) d x=u(x) v(x)-\int u^{\prime}(x) v(x) d x ∫u(x)v′(x)dx=u(x)v(x)−∫u′(x)v(x)dx).
那么，当 p ( x ) ≠ q ( x ) p(x) \neq q(x) p(x)=q(x)时，上式是不等于0的，这就有些像一个距离度量，我们进一步想到上式是不是也可以像KL divergence一样衡量两个分布之间的距离呢？并且通过优化这个距离做一些事情呢？这就是Liu Qiang组的研究工作。
2.向量情况下的距离度量
根据Stein’s identity我们可以定义一种度量来衡量分布p和q之间的距离，首先考虑向量情况下的推导，即 f ( x ) f(x) f(x)属于 R 1 × 1 \mathbb{R}^{1×1} R1×1
我们定义算子
A
p
f
(
x
)
=
s
p
(
x
)
f
(
x
)
+
∇
x
f
(
x
)
\mathcal{A}_{p} f(x)=s_{p}(x) f(x)+\nabla_{x} f(x)
Apf(x)=sp(x)f(x)+∇xf(x)，以及
s
q
=
∇
x
log
q
(
x
)
\boldsymbol{s}_{q}=\nabla_{x} \log q(x)
sq=∇xlogq(x)，之前的Stein’s identity就可以写成：
E
x
∼
p
[
A
q
f
(
x
)
]
=
0
w
h
e
r
e
(
q
=
p
)
\mathbb{E}_{x \sim p}\left[\mathcal{A}_{q} f(x)\right] =0\ \ where\ \ (q=p)
Ex∼p[Aqf(x)]=0 where (q=p)
我们定义距离 S ( p , q ) = max f ∈ F ( E p [ A q f ( x ) ] ) 2 {S}(p, q)=\max _{f \in \mathcal{F}}\left(\mathbb{E}_{p}\left[\mathcal{A}_{q} f(x)\right]\right)^{2} S(p,q)=maxf∈F(Ep[Aqf(x)])2
这里 F \mathcal{F} F是光滑函数的集合，但是从 F \mathcal{F} F中找到 f f f不是一件容易的事情，所以我们还需要别的trick（也就是kernel）。关于这里为什么是max，因为我们希望找到的度量，一定是在 p ( x ) ≠ q ( x ) p(x) \neq q(x) p(x)=q(x)时让这个距离值尽可能的大，这样才能在同一种度量下分的更开。
3.矩阵情况下的距离度量
下面我们尝试将 A p f ( x ) \mathcal{A}_{p} f(x) Apf(x)推广到矩阵的形式，也就是 f ( x ) f(x) f(x)推广到 R d × 1 \mathbb{R}^{d×1} Rd×1的形式。
这里为了区分我们把 f ( x ) f(x) f(x)属于 R d × 1 \mathbb{R}^{d×1} Rd×1记作 ϕ ( x ) \phi (x) ϕ(x)，令 ϕ ( x ) = [ ϕ 1 ( x ) , ⋯ , ϕ d ( x ) ] ⊤ \phi(x)=\left[\phi_{1}(x), \cdots, \phi_{d}(x)\right]^{\top} ϕ(x)=[ϕ1(x),⋯,ϕd(x)]⊤. 也就是说 ϕ ( x ) \phi (x) ϕ(x)代表着一组光滑函数，是一个向量值函数。
此时， A p ϕ ( x ) = s p ( x ) ϕ ⊤ ( x ) + ∇ x ϕ ( x ) \mathcal{A}_{p} \phi(x)=s_{p}(x) \phi^{\top}(x)+\nabla_{x} \phi(x) Apϕ(x)=sp(x)ϕ⊤(x)+∇xϕ(x)
Stein’s identity可以写成：
E
x
∼
p
[
A
q
ϕ
(
x
)
]
=
E
x
∼
p
[
A
q
ϕ
(
x
)
−
A
p
ϕ
(
x
)
]
=
E
x
∼
p
[
(
s
q
(
x
)
−
s
p
(
x
)
)
ϕ
⊤
(
x
)
]
\mathbb{E}_{x \sim p}\left[\mathcal{A}_{q} \phi(x)\right]=\mathbb{E}_{x \sim p}\left[\mathcal{A}_{q} \phi(x)-\mathcal{A}_{p} \phi(x)\right]=\mathbb{E}_{x \sim p}\left[(s_q(x)-s_p(x)) \phi^{\top}(x)\right]
Ex∼p[Aqϕ(x)]=Ex∼p[Aqϕ(x)−Apϕ(x)]=Ex∼p[(sq(x)−sp(x))ϕ⊤(x)]
由于 s q ( x ) s_q(x) sq(x)属于 R d × 1 \mathbb{R}^{d×1} Rd×1， ϕ ⊤ ( x ) \phi^{\top}(x) ϕ⊤(x)属于 R 1 × d \mathbb{R}^{1×d} R1×d，上式应该是一个 d × d d×d d×d的矩阵。但是我们想要的是距离度量，矩阵的期望是没有任何意义的，所以我们可以取标量为：
E p [ trace ( A q ϕ ( x ) ) ] = E p [ ( s q ( x ) − s p ( x ) ) ⊤ ϕ ( x ) ] \mathbb{E}_{p}\left[\operatorname{trace}\left(\mathcal{A}_{q}\phi(x)\right)\right]=\mathbb{E}_{p}\left[\left(\boldsymbol{s}_{q}(x)-\boldsymbol{s}_{p}(x)\right)^{\top} \phi(x)\right] Ep[trace(Aqϕ(x))]=Ep[(sq(x)−sp(x))⊤ϕ(x)]
此时， s q ⊤ ( x ) s_q^{\top}(x) sq⊤(x)属于 R 1 × d \mathbb{R}^{1×d} R1×d， ϕ ( x ) \phi(x) ϕ(x)属于 R d × 1 \mathbb{R}^{d×1} Rd×1，我们得到了一个标量，这是我们所希望的距离值。同时，也可以很直接地看到我们取标量后的值其实就是原矩阵形式的迹。
总结一下，到现在为止，我们将Stein’s identity从向量形式拓展到了矩阵形式，在矩阵形式下，我们依旧可以将 S ( p , q ) {S}(p, q) S(p,q)写成一个标量，为之后kernel的引入做好了铺垫。
4.Kernel的引入
4.1 RKHS
在引入kernel到Stein’s identity之前，我们先回顾一下RKHS(reproducing kernel hilbert space)。
每一个函数都可以看做一个无限维的向量，那么二元函数 K ( x , y ) K(x,y) K(x,y)就可以看做是一个无限维的矩阵。当 K ( x , y ) K(x,y) K(x,y)满足正定性和对称性的时候我们把它叫做核函数。与矩阵特征值和特征向量类似，核函数存在特征值和特征函数 (将函数看做无限维向量)。也就是
∫ K ( x , y ) ψ ( x ) d x = λ ψ ( y ) \int K(\mathrm{x}, \mathrm{y}) \psi(\mathrm{x}) d \mathrm{x}=\lambda \psi(\mathrm{y}) ∫K(x,y)ψ(x)dx=λψ(y)
那么根据Mercer定理，我们有
K
(
x
,
y
)
=
∑
i
=
0
∞
λ
i
ψ
i
(
x
)
ψ
i
(
y
)
K(\mathrm{x}, \mathrm{y})=\sum_{i=0}^{\infty} \lambda_{i} \psi_{i}(\mathrm{x}) \psi_{i}(\mathrm{y})
K(x,y)=i=0∑∞λiψi(x)ψi(y)
这里一个核函数对应无穷个特征值
{
λ
i
}
i
=
1
∞
\left\{\lambda_{i}\right\}_{i=1}^{\infty}
{λi}i=1∞和无穷个特征方程
{
ψ
i
}
i
=
1
∞
\left\{\psi_{i}\right\}_{i=1}^{\infty}
{ψi}i=1∞
将
{
λ
i
ψ
i
}
i
=
1
∞
\left\{\sqrt{\lambda_{i}} \psi_{i}\right\}_{i=1}^{\infty}
{λi
ψi}i=1∞作为一组正交基构建一个希尔伯特空间
H
\mathcal{H}
H 。这个空间中的任何一个函数（向量）都可以表示为这组基的线性组合。如
f
=
(
f
1
,
f
2
,
…
)
H
T
f=\left(f_{1}, f_{2}, \ldots\right)_{\mathcal{H}}^{T}
f=(f1,f2,…)HT
用
K
(
x
,
⋅
)
K(\mathrm{x}, \cdot)
K(x,⋅)表示固定核函数的一个参数为
x
x
x，即矩阵第
x
x
x行的一元函数或无限维向量。那么，我们有
K
(
x
,
⋅
)
=
(
λ
1
ψ
1
(
x
)
,
λ
2
ψ
2
(
x
)
,
…
)
H
T
K
(
y
,
⋅
)
=
(
λ
1
ψ
1
(
y
)
,
λ
2
ψ
2
(
y
)
,
…
)
H
T
<
K
(
x
,
⋅
)
,
K
(
y
,
⋅
)
>
H
=
∑
i
=
0
∞
λ
i
ψ
i
(
x
)
ψ
i
(
y
)
=
K
(
x
,
y
)
\begin{array}{c} K(\mathrm{x}, \cdot)=\left(\sqrt{\lambda_{1}} \psi_{1}(\mathrm{x}), \sqrt{\lambda_{2}} \psi_{2}(\mathrm{x}), \ldots\right)_{\mathcal{H}}^{T} \\ K(\mathrm{y}, \cdot)=\left(\sqrt{\lambda_{1}} \psi_{1}(\mathrm{y}), \sqrt{\lambda_{2}} \psi_{2}(\mathrm{y}), \ldots\right)_{\mathcal{H}}^{T} \\ <K(\mathbf{x}, \cdot), K(\mathbf{y}, \cdot)>_{\mathcal{H}}=\sum_{i=0}^{\infty} \lambda_{i} \psi_{i}(\mathbf{x}) \psi_{i}(\mathbf{y})=K(\mathbf{x}, \mathbf{y}) \end{array}
K(x,⋅)=(λ1
ψ1(x),λ2
ψ2(x),…)HTK(y,⋅)=(λ1
ψ1(y),λ2
ψ2(y),…)HT<K(x,⋅),K(y,⋅)>H=∑i=0∞λiψi(x)ψi(y)=K(x,y)
这就是RKHS，核心思想就是核函数的分解与重构。
4.2 Kernelized Stein Discrepancy（KSD）
现在，在RKHS的基础上，我们将kernel引入Stein’s identity
定义距离kernelized Stein discrepancy（KSD）:
S
(
p
,
q
)
=
E
x
,
y
∼
p
[
(
s
q
(
x
)
−
s
p
(
x
)
)
T
k
(
x
,
y
)
(
s
q
(
x
)
−
s
p
(
x
)
)
]
S(p, q)=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{x})\right)^{T} k(\mathbf{x}, \mathbf{y})\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{x})\right)\right]
S(p,q)=Ex,y∼p[(sq(x)−sp(x))Tk(x,y)(sq(x)−sp(x))]
此时的 S ( p , q ) S(p, q) S(p,q)本质上就是上一节矩阵情况中最后式子的变换， k ( x , y ) ( s q ( x ) − s p ( x ) ) k(\mathbf{x}, \mathbf{y})\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{x})\right) k(x,y)(sq(x)−sp(x))可以看成 ϕ ( x ) \phi(x) ϕ(x)，这一点从维度分析上易知。同时，我们还可以换一种理解方式， ϕ ( x ) \phi(x) ϕ(x)原来是一组函数，我们要找这组函数，那么现在变成了 k ( x , y ) ( s q ( x ) − s p ( x ) ) k(\mathbf{x}, \mathbf{y})\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{x})\right) k(x,y)(sq(x)−sp(x))，我们同样要找的是 x x x固定时的一组列向量函数。
可行性证明
但是，这个时候出现了一个问题，我们原来的 S ( p , q ) = max f ∈ F ( E p [ A q f ( x ) ] ) 2 {S}(p, q)=\max _{f \in \mathcal{F}}\left(\mathbb{E}_{p}\left[\mathcal{A}_{q} f(x)\right]\right)^{2} S(p,q)=maxf∈F(Ep[Aqf(x)])2中 E p \mathbb{E}_{p} Ep内的式子只与 q q q有关，而引入了kernel后， S ( p , q ) S(p, q) S(p,q)中 E p \mathbb{E}_{p} Ep内的式子与 q , p q,p q,p都有关了， p p p是我们不希望出现的，这会导致不可解。根据变换证明，这个问题是不存在的，引入kernel后的 E p \mathbb{E}_{p} Ep中同样可以只与 q q q有关，证明如下：
S ( p , q ) = E x , y ∼ p [ ( s q − s p ) T k ( x , y ) ( s q − s p ) ] = E x , y ∼ p [ ( s q − s p ) T ( k ( x , y ) s q + ∇ y k ( x , y ) − k ( x , y ) s p − ∇ y k ( x , y ) ) ] = E x , y ∼ p [ ( s q − s p ) T v ( x , y ) ] \begin{aligned} S(p, q) &=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[\left(s_{q}-s_{p}\right)^{T} k(\mathbf{x}, \mathbf{y})\left(s_{q}-s_{p}\right)\right] \\ &=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[\left(s_{q}-s_{p}\right)^{T}\left(k(\mathbf{x}, \mathbf{y}) s_{q}+\nabla_{y} k(\mathbf{x}, \mathbf{y})-k(\mathbf{x}, \mathbf{y}) s_{p}-\nabla_{\mathbf{y}} k(\mathbf{x}, \mathbf{y})\right)\right] \\ &=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[\left(s_{q}-s_{p}\right)^{T} v(\mathbf{x}, \mathbf{y})\right] \end{aligned} S(p,q)=Ex,y∼p[(sq−sp)Tk(x,y)(sq−sp)]=Ex,y∼p[(sq−sp)T(k(x,y)sq+∇yk(x,y)−k(x,y)sp−∇yk(x,y))]=Ex,y∼p[(sq−sp)Tv(x,y)]
这里 v ( x , y ) = k ( x , y ) s q ( y ) + ∇ y k ( x , y ) = A q k x ( y ) v(\mathbf{x}, \mathbf{y})=k(\mathbf{x}, \mathbf{y}) s_{q}(\mathbf{y})+\nabla_{\mathbf{y}} k(\mathbf{x}, \mathbf{y})=\mathcal{A}_{q} k_{\mathbf{x}}(\mathbf{y}) v(x,y)=k(x,y)sq(y)+∇yk(x,y)=Aqkx(y)， k x ( ⋅ ) = k ( x , ⋅ ) k_{\mathbf{x}} (\cdot)=k(\mathbf{x}, \cdot) kx(⋅)=k(x,⋅)
为了让
p
p
p消失，将
s
p
s_p
sp带入
S
(
p
,
q
)
=
E
x
,
y
∼
p
[
s
q
T
v
(
x
,
y
)
−
(
∇
x
ln
p
(
x
)
)
T
v
(
x
,
y
)
]
=
E
x
,
y
∼
p
[
s
q
T
v
(
x
,
y
)
]
−
∫
d
x
d
y
p
(
x
)
p
(
y
)
(
∇
x
ln
p
(
x
)
)
T
v
(
x
,
y
)
=
E
x
,
y
∼
p
[
s
q
T
v
(
x
,
y
)
]
+
E
x
,
y
∼
p
[
tr
∇
x
v
(
x
,
y
)
]
=
E
x
,
y
∼
p
[
u
q
(
x
,
y
)
]
\begin{aligned} S(p, q) &=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[s_{q}^{T} v(\mathbf{x}, \mathbf{y})-\left(\nabla_{\mathbf{x}} \ln p(\mathbf{x})\right)^{T} v(\mathbf{x}, \mathbf{y})\right] \\ &=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[s_{q}^{T} v(\mathbf{x}, \mathbf{y})\right]-\int d \mathbf{x} d \mathbf{y} p(\mathbf{x}) p(\mathbf{y})\left(\nabla_{\mathbf{x}} \ln p(\mathbf{x})\right)^{T} v(\mathbf{x}, \mathbf{y}) \\ &=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[s_{q}^{T} v(\mathbf{x}, \mathbf{y})\right]+\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[\operatorname{tr} \nabla_{\mathbf{x}} v(\mathbf{x}, \mathbf{y})\right] \\ &=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[u_{q}(\mathbf{x}, \mathbf{y})\right] \end{aligned}
S(p,q)=Ex,y∼p[sqTv(x,y)−(∇xlnp(x))Tv(x,y)]=Ex,y∼p[sqTv(x,y)]−∫dxdyp(x)p(y)(∇xlnp(x))Tv(x,y)=Ex,y∼p[sqTv(x,y)]+Ex,y∼p[tr∇xv(x,y)]=Ex,y∼p[uq(x,y)]
这里 u q ( x , y ) = s q ( x ) T k ( x , y ) s q ( y ) + s q ( x ) T ∇ y k ( x , y ) + ( ∇ x k ( x , y ) ) T s q ( y ) + tr ( ∇ x ∇ y k ( x , y ) ) u_{q}(\mathbf{x}, \mathbf{y})=s_{q}(\mathbf{x})^{T} k(\mathbf{x}, \mathbf{y}) s_{q}(\mathbf{y})+s_{q}(\mathbf{x})^{T} \nabla_{\mathbf{y}} k(\mathbf{x}, \mathbf{y})+\left(\nabla_{\mathbf{x}} k(\mathbf{x}, \mathbf{y})\right)^{T} s_{q}(\mathbf{y})+\operatorname{tr}\left(\nabla_{\mathbf{x}} \nabla_{\mathbf{y}} k(\mathbf{x}, \mathbf{y})\right) uq(x,y)=sq(x)Tk(x,y)sq(y)+sq(x)T∇yk(x,y)+(∇xk(x,y))Tsq(y)+tr(∇x∇yk(x,y))
虽然很长Orz但是只和q有关了ao。到现在为止，我们证明了引入kernel不影响求解，引入kernel后与引入前的原式形式等价。
易于求解证明
引入kernel自然是有引入的好处，好处在刚开始就已经说过了因为易与求解。光滑函数函数集 F \mathcal{F} F中直接找使距离最大的 f f f是很难做到的，但是引入kernel后我们就很容易找到了。我们下面来看为什么容易找到。
令 β ( y ) = E x ∼ p [ A q k y ( x ) ] = E x ∼ p [ ( s q ( x ) − s p ( x ) ) k y ( x ) ] \boldsymbol{\beta}(\mathbf{y})=\mathbb{E}_{\mathbf{x} \sim p}\left[\mathcal{A}_{q} k_{\mathbf{y}}(\mathbf{x})\right]=\mathbb{E}_{\mathbf{x} \sim p}\left[\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{x})\right) k_{\mathbf{y}}(\mathbf{x})\right] β(y)=Ex∼p[Aqky(x)]=Ex∼p[(sq(x)−sp(x))ky(x)]
根据RKHS的分析，我们可以将KSD写成
S
(
p
,
q
)
=
E
x
,
y
∼
p
[
(
s
q
(
x
)
−
s
p
(
y
)
)
T
k
(
x
,
y
)
(
s
q
(
x
)
−
s
p
(
y
)
)
]
=
E
x
,
y
∼
p
[
(
s
q
(
x
)
−
s
p
(
y
)
)
T
<
k
(
x
,
⋅
)
,
k
(
⋅
,
y
)
>
H
(
s
q
(
x
)
−
s
p
(
y
)
)
]
=
∑
i
=
1
d
<
E
x
∼
p
[
(
s
q
i
(
x
)
−
s
p
i
(
x
)
)
k
(
x
,
⋅
)
]
,
E
y
∼
p
[
(
s
q
i
(
y
)
−
s
p
i
(
y
)
)
k
(
⋅
y
)
]
>
H
=
∑
i
=
1
d
<
β
i
,
β
i
>
H
=
∥
β
∥
H
d
2
\begin{aligned} S(p, q) &=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{y})\right)^{T} k(\mathbf{x}, \mathbf{y})\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{y})\right)\right] \\ &=\mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{y})\right)^{T}<k(\mathbf{x}, \cdot), k(\cdot, \mathbf{y})>_{\mathcal{H}}\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{y})\right)\right] \\ &=\sum_{i=1}^{d}<\mathbb{E}_{\mathbf{x} \sim p}\left[\left(s_{q}^{i}(\mathbf{x})-s_{p}^{i}(\mathbf{x})\right) k(\mathbf{x}, \cdot)\right], \mathbb{E}_{\mathbf{y} \sim p}\left[\left(s_{q}^{i}(\mathbf{y})-s_{p}^{i}(\mathbf{y})\right) k(\cdot \mathbf{y})\right]>_{\mathcal{H}} \\ &=\sum_{i=1}^{d}<\beta_{i}, \beta_{i}>_{\mathcal{H}}\\ &=\|\boldsymbol{\beta}\|^2_{\mathcal{H}^{d}} \end{aligned}
S(p,q)=Ex,y∼p[(sq(x)−sp(y))Tk(x,y)(sq(x)−sp(y))]=Ex,y∼p[(sq(x)−sp(y))T<k(x,⋅),k(⋅,y)>H(sq(x)−sp(y))]=i=1∑d<Ex∼p[(sqi(x)−spi(x))k(x,⋅)],Ey∼p[(sqi(y)−spi(y))k(⋅y)]>H=i=1∑d<βi,βi>H=∥β∥Hd2
也就是说，我们可以将 S ( p , q ) S(p,q) S(p,q)看成希尔伯特空间中 β \beta β函数的内积了。此时易知，任意的希尔伯特空间中的归一化（归一化目的：防止max导致无穷大）函数 ϕ \phi ϕ与 β \beta β的内积都小于 β \beta β与自己的内积。那么，我们想找的max，就可以写成
∥ β ∥ H d = S ( p , q ) = max ϕ ∈ H d { E x ∼ p [ tr ( A q ϕ ( x ) ) ] , s.t. ∥ ϕ ∥ H d ≤ 1 } \|\boldsymbol{\beta}\|_{\mathcal{H}^{d}}=S(p, q)=\max _{\mathbf{\phi} \in \mathcal{H}^{d}}\left\{\mathbb{E}_{\mathbf{x} \sim p}\left[\operatorname{tr}\left(\mathcal{A}_{q} \mathbf{\phi}(\mathbf{x})\right)\right], \quad \text { s.t. } \quad\|\mathbf{\phi}\|_{\mathcal{H}^{d}} \leq 1\right\} ∥β∥Hd=S(p,q)=ϕ∈Hdmax{Ex∼p[tr(Aqϕ(x))], s.t. ∥ϕ∥Hd≤1}
此时最大值 ϕ ∗ = β / ∥ β ∥ H d \mathbf{\phi}^{*}=\boldsymbol{\beta} /\|\boldsymbol{\beta}\|_{\mathcal{H}^{d}} ϕ∗=β/∥β∥Hd
这也就是说我们用kernel可以直接导出想找的 ϕ \phi ϕ
5.SVGD算法
那么有了KSD之后，我们能做什么呢？一个简单的想法是用它来优化KL散度，导师说这个式子 ∇ ϵ K L ( q [ T ] ∥ p ) ∣ ϵ = 0 = − E x ∼ q [ trace ( A p ϕ ( x ) ) ] \left.\nabla_{\epsilon} \mathrm{KL}\left(q_{[T]} \| p\right)\right|_{\epsilon=0}=-\mathbb{E}_{x \sim q}\left[\operatorname{trace}\left(\mathcal{A}_{p} \phi(x)\right)\right] ∇ϵKL(q[T]∥p)∣∣ϵ=0=−Ex∼q[trace(Apϕ(x))]在learning领域是很常见的，也就是说KL散度变分求导是等于KSD的。所以我其实在想是不是这些论文都是根据KL散度反推回去的…其实它们都是拿结果往回推起源…然后再填充证明…
5.1 KL divergence的联系
下面我们来证明 ∇ ϵ K L ( q [ T ] ∥ p ) ∣ ϵ = 0 = − E x ∼ q [ trace ( A p ϕ ( x ) ) ] \left.\nabla_{\epsilon} \mathrm{KL}\left(q_{[T]} \| p\right)\right|_{\epsilon=0}=-\mathbb{E}_{x \sim q}\left[\operatorname{trace}\left(\mathcal{A}_{p} \phi(x)\right)\right] ∇ϵKL(q[T]∥p)∣∣ϵ=0=−Ex∼q[trace(Apϕ(x))]
我们定义 z = T ϵ ( x ) z={T}_{\epsilon}(x) z=Tϵ(x), T ϵ ( x ) = T ( x ) {T}_{\epsilon}(x)=T(x) Tϵ(x)=T(x), T ( x ) = x + ϵ ϕ ( x ) {T}(x)=x+\epsilon \phi(x) T(x)=x+ϵϕ(x)，当 x ∼ q x \sim q x∼q时， q [ T ] q_{[T]} q[T]表示z的概率密度。
同样地，取反函数，当 x ∼ p x \sim p x∼p时， p [ T − 1 ] p_{[T^{-1}]} p[T−1]表示 z = T − 1 ( x ) z={T}^{-1}(x) z=T−1(x)的概率密度。
也可以理解成用微小扰动 ϵ ϕ ( x ) \epsilon \phi(x) ϵϕ(x)来让 q q q逼近 p p p
我们有
K
L
(
q
[
T
]
∥
p
)
=
K
L
(
q
∥
p
[
T
−
1
]
)
∇
ϵ
K
L
(
q
[
T
]
∥
p
)
=
−
E
x
∼
q
[
∇
ϵ
log
p
[
T
−
1
]
(
x
)
]
\begin{array}{c} \mathrm{KL}\left(q_{[T]} \| p\right)=\mathrm{KL}\left(q \| p_{\left[T^{-1}\right]}\right) \\ \nabla_{\epsilon} \mathrm{KL}\left(q_{[T]} \| p\right)=-\mathbb{E}_{x \sim q}\left[\nabla_{\epsilon} \log p_{\left[T^{-1}\right]}(x)\right] \end{array}
KL(q[T]∥p)=KL(q∥p[T−1])∇ϵKL(q[T]∥p)=−Ex∼q[∇ϵlogp[T−1](x)]
将
∇
ϵ
log
p
[
T
−
1
]
(
x
)
\nabla_{\epsilon} \log p_{\left[T^{-1}\right]}(x)
∇ϵlogp[T−1](x)展开得到
∇
ϵ
log
p
[
T
−
1
]
(
x
)
=
s
p
(
T
(
x
)
)
⊤
∇
ϵ
T
(
x
)
+
trace
(
(
∇
x
T
(
x
)
)
−
1
⋅
∇
ϵ
∇
x
T
(
x
)
)
\nabla_{\epsilon} \log p_{\left[T^{-1}\right]}(x)=s_{p}(\boldsymbol{T}(x))^{\top} \nabla_{\epsilon} \boldsymbol{T}(x)+\operatorname{trace}\left(\left(\nabla_{x} \boldsymbol{T}(x)\right)^{-1} \cdot \nabla_{\epsilon} \nabla_{x} \boldsymbol{T}(x)\right)
∇ϵlogp[T−1](x)=sp(T(x))⊤∇ϵT(x)+trace((∇xT(x))−1⋅∇ϵ∇xT(x))
当扰动
ϵ
\epsilon
ϵ很小，
ϵ
=
0
\epsilon=0
ϵ=0的时候
T
(
x
)
=
x
,
∇
ϵ
T
(
x
)
=
ϕ
(
x
)
,
∇
x
T
(
x
)
=
I
,
∇
ϵ
∇
x
T
(
x
)
=
∇
x
ϕ
(
x
)
\boldsymbol{T}(x)=x, \quad \nabla_{\epsilon} \boldsymbol{T}(x)=\boldsymbol{\phi}(x), \quad \nabla_{x} \boldsymbol{T}(x)=I, \quad \nabla_{\epsilon} \nabla_{x} \boldsymbol{T}(x)=\nabla_{x} \boldsymbol{\phi}(x)
T(x)=x,∇ϵT(x)=ϕ(x),∇xT(x)=I,∇ϵ∇xT(x)=∇xϕ(x).
代入上式可得 ∇ ϵ log p [ T − 1 ] ( x ) = trace ( A p ϕ ( x ) ) \nabla_{\epsilon} \log p_{\left[T^{-1}\right]}(x)=\operatorname{trace}\left(\mathcal{A}_{p} \phi(x)\right) ∇ϵlogp[T−1](x)=trace(Apϕ(x))
所以有 ∇ ϵ K L ( q [ T ] ∥ p ) ∣ ϵ = 0 = − E x ∼ q [ trace ( A p ϕ ( x ) ) ] \left.\nabla_{\epsilon} \mathrm{KL}\left(q_{[T]} \| p\right)\right|_{\epsilon=0}=-\mathbb{E}_{x \sim q}\left[\operatorname{trace}\left(\mathcal{A}_{p} \phi(x)\right)\right] ∇ϵKL(q[T]∥p)∣∣ϵ=0=−Ex∼q[trace(Apϕ(x))]
即证明了KL散度变化最快的方向就是KSD所对应的向量函数 ϕ ∗ = β / ∥ β ∥ H d \mathbf{\phi}^{*}=\boldsymbol{\beta} /\|\boldsymbol{\beta}\|_{\mathcal{H}^{d}} ϕ∗=β/∥β∥Hd
5.2 Algorithm
有了梯度方向，有了kernel，我们就可以设计算法了，就是SVGD的实现。
算法流程如下所示：
算法的实现比较清晰明了，相比之前的kernel引入推导更容易快速理解一些。 p ( x ) p(\mathbf{x}) p(x) 是我们想要逼近的 R d \mathbb{R}^{d} Rd 分布, 我们想要用若干粒子来做 p ( x ) p(\mathbf{x}) p(x)上的采样。选取粒子群 { x i } i = 1 n ⊂ R d \left\{\mathbf{x}_{i}\right\}_{i=1}^{n} \subset \mathbb{R}^{d} {xi}i=1n⊂Rd . 使用梯度下降法，设置学习率为 ϵ \epsilon ϵ, 每一步的梯度为 ϕ ( x i ) \phi\left(\mathbf{x}_{i}\right) ϕ(xi).不断迭代直至 l l l达到指定次数或梯度变化达到阈值。
总的来说，这个idea感觉很赞，我觉得从 E x ∼ p [ ( s q ( x ) − s p ( x ) ) ϕ ⊤ ( x ) ] \mathbb{E}_{x \sim p}\left[(s_q(x)-s_p(x)) \phi^{\top}(x)\right] Ex∼p[(sq(x)−sp(x))ϕ⊤(x)]看到 E x , y ∼ p [ ( s q ( x ) − s p ( x ) ) T k ( x , y ) ( s q ( x ) − s p ( x ) ) ] \mathbb{E}_{\mathbf{x}, \mathbf{y} \sim p}\left[\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{x})\right)^{T} k(\mathbf{x}, \mathbf{y})\left(s_{q}(\mathbf{x})-s_{p}(\mathbf{x})\right)\right] Ex,y∼p[(sq(x)−sp(x))Tk(x,y)(sq(x)−sp(x))]这一步真的很强，毕竟是ICML、NeurIPS连着中的组（
References
[1]Liu Q , Wang D . Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm[C]// 2016.
[2]Liu Q , Jason D.Lee . A Kernelized Stein Discrepancy for Goodness-of-fit Tests[C]// 2016.
[3]Liu Q , Stein Variational Gradient Descent as Gradient Flow[C]// 2016.
[4]Stein variational gradient descent https://zhuanlan.zhihu.com/p/114489837