阅读笔记八：Efficient Learning of the Parameters of Non-Linear Models using Differentiable Resampling

文中多次提到对数似然方程，所以我先介绍一下对数似然方程。
当总体X为连续型随机变量时，设其分布密度为$f(x;{\theta }_1,{\theta }_1,...,{\theta }_m)$，其中${\theta }_1,{\theta }_1,...,{\theta }_m$为未知参数。又设$x_1,x_2,...,x_n$为总体的一个样本，称:$$L({\theta }_1,{\theta }_1,...,{\theta }_m)=\prod _{i=1}^{n}f(x;{\theta }_1,{\theta }_1,...,{\theta }_m)$$为样本的似然函数，简记为$L_n$
当总体X为离散型随机变量时，设其分布律为$PX=x=p(x;{\theta }_1,{\theta }_1,...,{\theta }_m)$，则称$$L(x_1,x_2,...,x_n;{\theta }_1,{\theta }_1,...,{\theta }_m)=\prod _{i=1}^{n}p(x;{\theta }_1,{\theta }_1,...,{\theta }_m)$$为样本的似然函数
若似然函数$L(\theta ;x)$为$\theta$的连续函数，且关于$\theta$的各分量的偏导数存在。设$\theta$是m维变量，且$\Theta \in R^m$为开区域，则由极值的一阶必要条件，得到$$\frac{\partial L(\theta ;x)}{\partial {\theta }_i}=0,i=1,2,...,m$$通常称为似然方程，由于独立同分布的样本的似然函数上$L(\theta ;x)$具有连乘积的形式，故对$L(\theta ;x)$取对数后再求偏导数是方便的，因此实用上常采用与似然方程等价的形式 ：$$\frac{\partial lnL(\theta ;x)}{\partial {\theta }_i}=0,i=1,2,...,m$$称为对数似然方程。

简介：

状态空间模型(SSMs)已被用于在广泛的研究领域建模动力系统。状态空间模型(SSMs)由两个随机过程表示: { X t } t ≥ 0 \{X_t\}_{t≥0} {Xt​}t≥0​和 { Y t } t ≥ 0 \{Y_t\}_{t≥0} {Yt​}t≥0​，其中$X_t$表示根据马尔可夫过程$p(x_t|x_{t-1})$演化的隐藏状态，$Y_t$表示观测值。用式子表示为：$$X_t|X_{t-1}\sim p(x_t|x_{t-1},\theta )（公式1）$$ $$Y_t|X_t\sim p(y_t|x_t,\theta )（公式2）$$初始$X_0$的初始密度表示$\mu \theta (x0)$，SSM由参数空间$\Theta$中包含的未知静态参数$\theta$参数化。

在这篇文章中，我们专注于对状态空间模型(SSMs)中的贝叶斯参数估计使用马尔可夫链蒙特卡洛 (p-MCMC)，这种方法结合了两种蒙特卡罗方法（一种是Markov Chain Monte Carlo (MCMC)，另一种是Sequential MonteCarlo (SMC) ），它们使用重复采样技术来获得目标分布$\pi (\theta )$的数值估计，但是这种方法的精确推断是不好处理的。
Markov Chain Monte Carlo (MCMC)算法，例如如Metropolis-Hastings (M-H)，在该算法中经常使用随机游走抽样。但是当估计大量参数时，这种随机游走抽样的做法很难使MCMC达到平稳分布。
Hamiltonian Monte Carlo (HMC)算法是一种从特定问题哈密顿系统模拟生成MCMC中使用的样本的方法。HMC在目标分布复杂、超参数敏感由用户决定的情况下被认为是最有效的。
粒子滤波器最初是用来计算$\theta$的无偏差估计，马尔可夫链蒙特卡洛（p-MCMC）中使用Metropolis-Hastings (M-H)随机游走算法进行采样将出现上面Markov Chain Monte Carlo (MCMC)中描述的相同问题（也就是很难达到平稳分布）。 于是又提出了重参数化技巧，通过提前对噪声向量进行采样，并将$\theta$的似然性定义为该采样噪声向量的确定性函数，将采样操作重新制定为可微函数。然而，重采样仍然是有问题的，因为重采样后所有的权重都是相等的。最后，人们又引出了一种新的网络架构——Soft resampling。它训练粒子转换器，然后取代传统的重采样。改架构涉及到分布$q(n)=\alpha w_{1:t}^{(\theta ，i)}+(1-\alpha )1/N，其中\alpha \in [0,1]$，表示一个权衡参数。如果$\alpha =1$，则使用常规重采样;如果$\alpha =0$，则执行子采样。我们就得到了一个新的权重：$$w_{1:t}^{\prime (\theta ,i)}=\frac{w_{1:t}^{(\theta ,i)}}{\alpha w_{1:t}^{(\theta ,i)}+(1-\alpha )1/N}（公式3）$$
通过改变α值，该方法用有偏差的梯度估计来交换重采样质量。
**我们相比之前的最优输运思想进行重采样、停止梯度重采样的方法等不同的地方在于我们关注的是确保重采样是可微的，而不必改变重采样的操作方式。**我们的核心是解决在重采样步骤中使用的随机数。 在上面条件的对重采样的输入会使对后续的粒子导数计算是父粒子的函数。然后，可以在p-MCMC框架内有效地估计和利用梯度。

二、粒子滤波的背景

在第2节中，我们描述了一个通用的粒子滤波器，然后描述了与区分采样和重采样步骤相关的困难。
力度变化
假设我们考虑了$t$个时间步，在$y_{1:t}$的每一个增量处获得数据。状态序列$x_{1:t}$随时间增长，其中$x_t$有$n_x$维。力度变化和可能性用$\theta$(它有$n_{\theta }$维)来参数化，如下所示：$$p\left(y_{1:t},x_{1:t}\mid \theta \right)=p\left(y_1\mid x_1,\theta \right)p\left(x_1\mid \theta \right)\prod _{\tau =2}^{t}p\left(y_{\tau }\mid x_{\tau },\theta \right)p\left(x_{\tau }\mid x_{\tau -1},\theta \right).（公式4）$$
2.1粒子滤波
在每个时间间隔$t$，粒子滤波器从分布$q(x_{1:t}|y_{1:t},\theta )$中提取N个样本(粒子)，该分布由状态和测量的序列参数化。这些样本在统计上是独立的，每个样本都代表了系统状态序列的不同假设。$x_{1:t}$的样本集合表示动态系统的概率密度函数。第$i$个样本有一个相关的权重$w_t^{(\theta ，i)}$，它表示第i个样本$x_t^{(\theta ，i)}$是系统真实状态的相对概率。t = 0时的权值设为$1/N$。分布递归构造为:$$q(x_{1:t}|y_{1:t},\theta )=q(x_1|y_1,\theta )\prod _{\tau =2}^{t}q(x_{\tau }|x_{\tau -1},y_{\tau },\theta ),（公式5）$$我们就可以对联合分布$p\left(y_{1:t},x_{1:t}\mid \theta \right)$进行估计,如下所示：$$\int p\left(y_{1:t},x_{1:t}\mid \theta \right)f\left(x_{1:t}\right)d{x}_{1:t}\approx \frac{1}{N}\sum _{i=1}^N{w}_{1:t}^{(\theta ,i)}f\left(x_{1:t}^{(i)}\right)（公式6）$$当t>1时，这是一个不存在偏差的估计：$$\begin{array}{rl}{w}_{1:t}^{(\theta ,i)}& =\frac{p\left(y_1\mid x_1^{(\theta ,i)},\theta \right)p\left(x_1^{(\theta ,i)}\mid \theta \right){\prod }_{\tau =2}^{t}p\left(y_{\tau }\mid x_{\tau }^{(\theta ,i)},\theta \right)p\left(x_{\tau }^{(\theta ,i)}\mid x_{\tau -1}^{(\theta ,i)},\theta \right)}{q\left(x_1^{(\theta ,i)}\mid y_1,\theta \right){\prod }_{\tau =2}^{t}q\left(x_{\tau }^{(\theta ,i)}\mid x_{\tau -1}^{(\theta ,i)},y_{\tau },\theta \right)}\\ & =w_{1:t-1}^{(\theta ,i)}\frac{p\left(y_t\mid x_t^{(\theta ,i)},\theta \right)p\left(x_t^{(\theta ,i)}\mid x_{t-1}^{(\theta ,i)},\theta \right)}{q\left(x_t^{(\theta ,i)}\mid x_{t-1}^{(\theta ,i)},y_t\right)},\end{array}（公式7）$$当t=1时，非归一化权值$w_{1:t}^{(\theta ,i)}$的递归公式与增加的权重：$$\begin{gathered} \sigma \left(x_k^{(\theta ,i)},x_{k-1}^{(\theta ,i)},\theta \right)=\frac{p\left(y_t\mid x_t^{(\theta ,i)},\theta \right)p\left(x_t^{(\theta ,i)}\mid x_{t-1}^{(\theta ,i)},\theta \right)}{q\left(x_t^{(\theta ,i)}\mid x_{t-1}^{(\theta ,i)},y_t\right)}. \\ \sigma \left(x_{1:1}^{(\theta ,i)}\right)=\frac{p\left(y_1\mid x_1^{(\theta ,i)},\theta \right)p\left(x_1^{(\theta ,i)}\mid \theta \right)}{q\left(x_1^{(\theta ,i)}\mid y_1\right)}（公式8） \end{gathered}$$

关于后验的估计

关于后验$p(x_{1:t}|y_{1:t},\theta )$的估计我们可以计算如下:
$$\int p(x_{1:t}|y_{1:t},\theta )f(x_{1:t})d{x}_{1:t}=\int \frac{p(y_{1:t},x_{1:t}|\theta )}{p(y_{1:t}|\theta )}f(x_{1:t})d{x}_{1:t}.（公式10）$$注意，如果$f(x_{1:t})=1$时：$$p(y_{1:t}|\theta )=\int p(y_{1:t},x_{1:t}|\theta )d{x}_{1:t}\approx \frac{1}{N}\sum _{i=1}^N{w}_{1:t}^{(\theta ,i)}（公式11）$$这样一来，公式6也可以写成：$$\begin{array}{rl}\int p\left(x_{1:t}\mid y_{1:t},\theta \right)f\left(x_{1:t}\right)d{x}_{1:t}& \approx \frac{1}{\frac{1}{N}{\sum }_{i=1}^N{w}_{1:t}^{(\theta ,i)}f\left(x_{1:t}^{(\theta ,i)}\right)}\frac{1}{N}\sum _{i=1}^N{w}_{1:t}^{(\theta ,i)}f\left(x_{1:t}^{(\theta ,i)}\right)\\ & =\sum _{i=1}^N{\tilde{w}}_{1:t}^{(\theta ,i)}f\left(x_{1:t}^{(\theta ,i)}\right)\end{array}（公式12)$$ 其中$${\tilde{w}}_{1:t}^{(\theta ,i)}=\frac{w_{1:t}^{(\theta ,i)}}{{\sum }_{i=1}^N{w}_{1:t}^{(\theta ,i)}}（公式13）$$
2.4重采样：
随着时间的推移，规范化的权重将变得越来越倾斜，这时就有人提出检测有效样本，而$N_{eff}$就可以用来确定是否需要重新采样，$N_{eff}$表达式为：$$N_{eff}=\frac{1}{{\sum }_{i=1}^N({\tilde{w}}_{1:t}^{(\theta ,i)}{)}^2}（公式14）$$许多重采样方法随机复制具有较高权重的粒子，同时消除较低权重的粒子，多项重采样是常用的方法，它涉及到从当前粒子集中抽取与其权重成比例的N倍的粒子，为了保持总非归一化权重不变，我们给每个新重采样的样本分配了一个非归一化权重：$$\frac{1}{N}{w}_{1:t}^{(\theta ,i)}（公式15）$$
注意重采样后归一化权值为：$\frac{1}{N}$

三、梯度和可能性的计算

对权重进行微分可以得到近似于可能性梯度的结果：$$\frac{d}{d\theta }p(y_{1:t}\mid \theta )=\frac{1}{N}\sum _{i=1}^N\frac{d}{d\theta }{w}_{1:t}^{(\theta ,i)}（公式16）$$为了数值的稳定性，通常最好在对数中传播值。将链式法则应用于(公式11)和(公式16)得到:$$\begin{array}{rl}\frac{d}{d\theta }\log p\left(y_{1:t}\mid \theta \right)& =\frac{1}{p\left(y_{1:t}\mid \theta \right)}\sum _{i=1}^N{w}_{1:t}^{(\theta ,i)}\frac{d}{d\theta }\log w_{1:t}^{(\theta ,i)}\\ & \approx \sum _{i=1}^N{\tilde{w}}_{1:t}^{(\theta ,i)}\frac{d}{d\theta }\log w_{1:t}^{(\theta ,i)}\end{array}（公式17）$$对数权值可以递归计算为：$$log{w}_{1:t}^{(\theta ,i)}=log{w}_{1:t-1}^{(\theta ,i)}+log\text{ σ }(x_k^{(\theta ,i)},x_{t-1}^{(\theta ,i)},\theta )（公式18）$$对公式18两边同时求微分可以得到$$\frac{d}{d\theta }log{w}_{1:t}^{(\theta ,i)}=\frac{d}{d\theta }log{w}_{1:t-1}^{(\theta ,i)}+\frac{d}{d\theta }log\text{ σ }(x_k^{(\theta ,i)},x_{t-1}^{(\theta ,i)},\theta )（公式18）$$
其中$$\frac{d}{d\theta }log\text{ σ }(x_k^{(\theta ,i)},x_{t-1}^{(\theta ,i)},\theta )=\frac{d}{d\theta }log\text{ p }(x_k^{(\theta ,i)},x_{t-1}^{(\theta ,i)},\theta )+\frac{d}{d\theta }log\text{ p }(y_k\mid x_{t-1}^{(\theta ,i)})-\frac{d}{d\theta }log\text{ q }(x_k^{(\theta ,i)}\mid x_{t-1}^{(\theta ,i)},\theta ,y_t)（公式19）$$如果我们能区分单一测量的可能性，转换模型和研究计划，我们可以近似计算下一个时间间隔的对数似然导数，从而递归逼近每个时间间隔的对数似然导数。
如果粒子滤波使用跃迁模型作为动力，权重更新中的可能性并不明确地依赖于$\theta$,我们最初可以假设$\frac{d}{d\theta }log\sigma =0$,这样的话，使用(公式18)的归纳论证将表明权导数总是零，因此θ的似然梯度近似为零。但是这种推论是有问题的，因为有些可能性是依赖$\theta$的，我们应用链式法则：
$$\frac{d}{d\theta }logp(y_t|x_t^{(\theta ,i)})=\frac{d}{dx}logp(y_t\mid x){\mid }_{x=x_t^{(\theta ,i)}d}\frac{d}{d\theta }{x}_t^{(\theta ,i)}（公式20）$$由于$x_t^{(\theta ，i)}$是从建议中采样的随机变量,我们使用重参数化技巧——设${ϵ}_t^{(i)}$为从建议中采样时使用的标准N(0,1)随机变量的向量,如果${ϵ}_t^{(i)}$是已知的，那么$x_t^{(\theta ，i)}$是$x_{t-1}^{(\theta ，i)}$的一个可微分的确定性函数，于是我们认为$$\begin{array}{rl}\frac{d}{d\theta }p\left(y_{1:t}\mid \theta \right)& =\frac{d}{d\theta }\int p\left(y_{1:t},{ϵ}_{1:t}\mid \theta \right)d{ϵ}_{1:t}\\ & =\int \frac{d}{d\theta }p\left(y_{1:t},{ϵ}_{1:t}\mid \theta \right)d{ϵ}_{1:t}\\ & \approx \frac{1}{N}\sum _{i=1}^N\frac{d}{d\theta }p\left(y_{1:t}\mid {ϵ}_{1:t}^{(i)},\theta \right)（公式21）\end{array}$$其中${ϵ}_t^{(i)}\sim p({ϵ}_{1:t}\mid \theta )$是固定的,并且公式21是可以计算出来的微分。

四、计算导数

为了传播粒子重量的导数，我们需要计算下面的式子：
粒子导数：$\frac{d{x}_t^{(\theta ,i)}}{d\theta }（公式22）$
概率密度函数的导数：$\frac{d}{d\theta }log\text{ q }(x_t^{(\theta ,i)}\mid x_{t-1}^{(\theta ,i)},\theta ,y_t)（公式23）$
对数概率密度函数的导数：$\frac{d}{d\theta }log\text{ p }(x_t^{(\theta ,i)}\mid x_{t-1}^{(\theta ,i)},\theta )（公式24）$
单测量似然对数概率密度函数的导数：$\frac{d}{d\theta }log\text{ p }(y_t\mid x_{t-1}^{(\theta ,i)})（公式25）$
4.1粒子导数的计算
$$q\left(x_t^{(\theta ,i)}\mid x_{t-1}^{(\theta ,i)},\theta ,y_t\right)=\mathcal{N}\left(x_t^{(\theta ,i)};\mu \left(x_{t-1}^{(\theta ,i)},\theta ,y_t\right),C\left(x_{t-1}^{(\theta ,i)},\theta ,y_t\right)\right)（公式26）$$其中$\mu (\cdot )和C(\cdot )$是旧粒子状态的函数，测量值和参数。如果我们提前对噪声进行抽样${ϵ}_t^{(i)}\sim \mathcal{N}(\cdot ;0,I_{n_X})$,那么新的粒子状态可以写成一个确定性函数：$$x_k^{(\theta ,i)}=f\left(x_{t-1}^{(\theta ,i)},\theta ,y_t,{ϵ}_t^i\right)\triangleq \mu \left(x_{t-1}^{(\theta ,i)},\theta ,y_t\right)+\sqrt{C\left(x_{t-1}^{(\theta ,i)},\theta ,y_t\right)}\times {ϵ}_t^i\text{. }（公式27）$$计算的时候注意一点——$x_{t-1}^{(\theta ，i)}$本身是θ的函数。于是我们公式22计算公式为：$$\begin{array}{rl}\frac{d{x}_t^{(\theta ,i)}}{d\theta }& =\frac{d}{d\theta }f\left(x_{t-1}^{(\theta ,i)},\theta ,y_t,{ϵ}_t^i\right)\\ & =\frac{\partial f}{\partial x_{t-1}^{(\theta ,i)}}\frac{d{x}_{t-1}^{(\theta ,i)}}{d\theta }+\frac{\partial f}{\partial \theta }\frac{d\theta }{d\theta }\\ & =\frac{\partial f}{\partial x_{t-1}^{(\theta ,i)}}\frac{d{x}_{t-1}^{(\theta ,i)}}{d\theta }+\frac{\partial f}{\partial \theta }.（公式28）\end{array}$$
4.2概率密度函数的导数的计算：
为了跟对数概率密度函数的导数区分，可以写成：$$logq(x_t^{(\theta ,i)}\mid x_{t-1}^{(\theta ,i)},\theta ,y_t)=Q(x_{t-1}^{(\theta ,i)},\theta ,y_t,{ϵ}^{(}i{)}_t)（公式29）$$为了书写方便我们去掉固定值$y_t、{ϵ}^{(}i{)}_t$ $$\begin{gathered} Q(x_{t-1}^{(\theta ,i)},\theta )\triangleq logq（f(x_{t-1}^{(\theta ,i)},\theta )\mid x_{t-1}^{(\theta ,i)} \\ =log\mathcal{N}(f(x_{t-1}^{(\theta ,i)},\theta );\mu (x_{t-1}^{(\theta ,i)},\theta ),C(x_{t-1}^{(\theta ,i)},\theta ))（公式30） \end{gathered}$$我们假设它是满足高斯的,于是我们就有：$$\begin{array}{rl}\frac{d}{d\theta }Q\left(x_{t-1}^{(\theta ,i)},\theta \right)=& \frac{\partial }{\partial f}\log \mathcal{N}(f;\mu ,C)\left(\frac{df}{d\theta }+\frac{d\mu }{d\theta }+\frac{dC}{d\theta }\right)\\ =& \frac{\partial }{\partial f}\log \mathcal{N}(f;\mu ,C)\left(\frac{\partial f}{\partial x_{t-1}^{(\theta ,i)}}\frac{d{x}_{t-1}^{(\theta ,i)}}{d\theta }+\frac{\partial f}{\partial \theta }\right)+\\ & \frac{\partial }{\partial \mu }\log \mathcal{N}(f;\mu ,C)\left(\frac{\partial \mu }{\partial x_{t-1}^{(\theta ,i)}}\frac{d{x}_{t-1}^{(\theta ,i)}}{d\theta }+\frac{\partial \mu }{\partial \theta }\right)+\\ & \frac{\partial }{\partial C}\log \mathcal{N}(f;\mu ,C)\left(\frac{\partial C}{\partial x_{t-1}^{(\theta ,i)}}\frac{d{x}_{t-1}^{(\theta ,i)}}{d\theta }+\frac{\partial C}{\partial \theta }\right).（公式31）\end{array}$$
4.3对数概率密度函数的导数的计算
首先我们令:$$\begin{array}{rl}P\left(x_{t-1}^{(\theta ,i)},\theta ,y_t,{ϵ}_t^i\right)& \triangleq \log p\left(f\left(x_{t-1}^{(\theta ,i)},\theta ,y_t,{ϵ}_t^i\right)\mid x_{t-1}^{(\theta ,i)},\theta \right)\\ & =\log \mathcal{N}\left(f\left(x_{t-1}^{(\theta ,i)}\right);a\left(x_{t-1}^{(\theta ,i)},\theta \right),\Sigma (\theta )\right)（公式32）\end{array}$$我们假设过渡模型具有加强高斯噪声，它与x^{(θ，i)}_{t−1}无关。然后我们就得到$$\begin{array}{rl}\frac{d}{d\theta }P\left(x_{t-1}^{(\theta ,i)},\theta \right)=& \frac{\partial }{\partial f}\log \mathcal{N}(f;a,\Sigma )\left(\frac{\partial f}{\partial x_{t-1}^{(\theta ,i)}}\frac{d{x}_{t-1}^{(\theta ,i)}}{d\theta }+\frac{\partial f}{\partial \theta }\right)+\\ & \frac{\partial }{\partial a}\log \mathcal{N}(f;a,\Sigma )\left(\frac{\partial a}{\partial x_{t-1}^{(\theta ,i)}}\frac{d{x}_{t-1}^{(\theta ,i)}}{d\theta }+\frac{\partial a}{\partial \theta }\right)+\\ & \frac{\partial }{\partial \Sigma }\log \mathcal{N}(f;a,\Sigma )\left(\frac{\partial \Sigma }{\partial \theta }\right)（公式33）\end{array}$$其中$a=a(x_{t-1}^{(\theta ,i)},\theta )、\Sigma =\Sigma (\theta )$。
4.4单测量似然对数概率密度函数的导数的计算
首先令$$\begin{array}{rl}L\left(x_t^{(\theta ,i)},\theta ,y_t\right)& \triangleq \log p\left(y_t\mid x_t^{(\theta ,i)},\theta \right)\\ & =\log \mathcal{N}\left(y_t;h\left(x_t^{(\theta ,i)},\theta \right),R(\theta )\right)（公式34）\end{array}$$假设概率是高斯分布，方差独立于$x_t^{(\theta ,i)}$,于是我们有$$\begin{array}{rl}\frac{d}{d\theta }& L\left(x_t^{(\theta ,i)},\theta ,y_t\right)=\frac{\partial }{\partial h}\log \mathcal{N}\left(y_k;h,R\right)\left(\frac{\partial h}{\partial x_t^{(\theta ,i)}}\frac{d{x}_t^{(\theta ,i)}}{d\theta }+\frac{\partial h}{\partial \theta }\right)+\frac{\partial }{\partial R}\log \mathcal{N}\left(y_t;h,R\right)\frac{dR}{d\theta }.（公式35）\end{array}$$其中$h=h(x_t^{(\theta ,i)},\theta )、R=R(\theta )$。

总结

概述了在粒子滤波器中执行重采样步骤时如何扩展重新参数化技巧以使用公共随机数。这限制了计算HMC和NUTS中用于p-MCMC中提出新参数的梯度时遇到的不连续问题，如简介所说，机器学习领域最近引入了不同的可微重采样方法。在使用NUTS进行参数估计的背景下比较这些方法将是未来工作的另一个有趣的方向（总结还没有完善，后续继续写）。