博客笔记思路：

1. 关于扩散模型的直觉理解

见《Diffusion 扩散模型知识体系梳理（一）》

2. 从 DDPM 的角度看离散 Diffusion

见《Diffusion 扩散模型知识体系梳理（一）》

3. 从 Score-based 的角度再看 Diffusion

见《Diffusion 扩散模型知识体系梳理（二）》

4. 从 SDE 的角度统一

4.1. 理论简介

很好，非常好，相信目前为止，你已经对 DDPM 和 Score-based 这两种处理扩散模型的思路都各自有了一定的认识和见解了，下面我将要分享的是，有一种理论，可以将这两者在数学本质上进行神奇的统一，也就是随机微分方程 Stochastic Differential Equations, SDE 角度的演化 [8]。什么意思呢? 也就是说，无论是离散的 DDPM，还是通过添加不同量级噪声的 Score-Based Model，都可以抽象出一个过程，遵循随机微分方程的演化原则，从样本过渡到纯噪声。用来训练的变量也是 Score Function，损失函数做了类似的调整，采样推理的方法则变成了求解 Reverse SDE。下面我将会详细展开介绍。

4.2. 训练框架的统一

• 首先，让我们用一种较为统一的方式来描述一下随机微分方程的通用形式，下式表示 Ito SDE (伊藤随机过程)：$$d\text{x}=\mathbf{f}(\text{x},t)dt+g(t)d\text{w}$$ 其中 $\mathbf{f}(\text{x},t)$ 是一个向量值函数，即输出的是一个向量场，表示的是系统演化的漂移系数 Drift Coeffcient。$g(t)$ 表示的则是扩散系数 Diffusion Coeffcient，反映掺杂的随机布朗运动的系数函数。

• 需要注意一点，上述方程的解，是一个预测的时间序列 { x ( t ) } t ∈ [ 0 , T ] \{\textbf{x}(t)\}_{t\in [0,T]} {x(t)}t∈[0,T]​，其中关联 NCSN 的思路，可以做出如下类比（注意两者的方向是一致的，都是噪声递增，但起点编号不同：SDE 从 $t=0$ 开始，NCSN 从 $i=1$ 开始）：

  • 时间步 $t=0,1,...,T$ $\stackrel{类比}{\to }$ 加噪的量级索引 $i=1,2,...,L$；
  • 第 $t$ 时间步的状态概率分布 $p_t(\text{x})$ $\stackrel{类比}{\to }$ 加噪量级为 ${\sigma }_i$ 时刻的概率分布（$p_0(\text{x})=p(\text{x})$）；
  • 时间序列长度 $T\to \infty$ $\stackrel{类比}{\to }$ 噪声的波动性 ${\sigma }_L^2\to \infty$；
  • Time-Dependent Score-Based $s_{\theta }(\text{x},t)={\nabla }_{\text{x}}\log p_t(\text{x})$ $\stackrel{类比}{\to }$ Noise-Conditional Score-Based $s_{\theta }(\text{x},i)={\nabla }_{\text{x}}\log p_{{\sigma }_i}(\text{x})$；

• 有趣的一点是，SDE 的逆过程也一定是 SDE（这一结论最早由 Anderson [9] 给出），且表达式为，关于这一点的数学推导，见 Appendix C：$$d\mathbf{x}=\underset{\text{漂移项}}{\underbrace{\left[\mathbf{f}(\mathbf{x},t)-g^2(t){\nabla }_{\mathbf{x}}\log p_t(\mathbf{x})\right]}}d\tau +\underset{\text{扩散项}}{\underbrace{g(t)d\bar{\mathbf{w}}}}$$ 其中 $\tau$ 表示逆向时间（从 $T$ 回退到 $0$），$d\bar{\mathbf{w}}$ 是逆向布朗运动。此时我们的 Score-SDE 模型在采样的时候就不需要 Langevin Dynamics 了，而可以直接求解这个 Reverse SDE 即可。

• 此时，我的模型训练的目标，还是 $s_{\theta }(\text{x},t)={\nabla }_{\text{x}}\log p_t(\text{x})$，训练的损失函数为： $$\mathcal{J}(\theta )={\mathbb{E}}_{t\sim \mathcal{U}(0,T)}{\mathbb{E}}_{p_t(\mathbf{x})}\left[\lambda (t){\Vert {\nabla }_{\mathbf{x}}\log p_t(\mathbf{x})-{\mathbf{s}}_{\theta }(\mathbf{x},t)\Vert }_2^2\right],\text{\hspace{1em}(11)}$$ 需要注意的是，此处我们采样的时间是在 $(0,T)$ 之间的均匀采样获得，系数 $\lambda (t)$ 满足下面的关系 $$\lambda (t)\propto \frac{1}{\mathbb{E}[||{\nabla }_{\text{x}(t)}\log p(\text{x}(t)|\text{x}(0))|{|}_2^2]}$$ 其中有一组概念需要澄清：

  • $p_t(\text{x})$ 是 Marginal Distribution of all real data after adding noise to $t$；
  • $p(\text{x}(t)|\text{x}(0))$ 是 Conditional Distribution of $\text{x}(t)$ obtained by gradually adding noise to the real data $\text{x}(0)$；

4.3. 如何求解 Reverse SDE

• 好的，现在我们已经知道了训练的目标和损失函数，接下来最关键的问题就是：训练完成之后，如何利用学到的 $s_{\theta }(\text{x},t)$ 来进行采样生成？答案就是数值求解 Reverse SDE。回顾一下我们的 Reverse SDE：$$d\mathbf{x}=\left[\mathbf{f}(\mathbf{x},t)-g^2(t){\nabla }_{\mathbf{x}}\log p_t(\mathbf{x})\right]dt+g(t)d\bar{\mathbf{w}}$$ 其中 ${\nabla }_{\mathbf{x}}\log p_t(\mathbf{x})$ 就用我们训练好的 $s_{\theta }(\text{x},t)$ 来替代。那么问题就变成了：如何数值求解一个 SDE？

• 最直接的方法是 Euler-Maruyama 方法，这是 SDE 版本的欧拉法。对于一般的 SDE $d\mathbf{x}=\mathbf{a}(\mathbf{x},t)dt+b(t)d\mathbf{w}$，其离散化形式为：$${\mathbf{x}}_{t-\Delta t}={\mathbf{x}}_t+\mathbf{a}({\mathbf{x}}_t,t)\Delta t+b(t)\sqrt{|\Delta t|}\mathbf{z},\mathbf{z}\sim \mathcal{N}(0,\mathbf{I})$$ 代入 Reverse SDE 的具体形式，我们得到每一步的更新规则：$${\mathbf{x}}_{t-\Delta t}={\mathbf{x}}_t+\left[\mathbf{f}({\mathbf{x}}_t,t)-g^2(t)s_{\theta }({\mathbf{x}}_t,t)\right]\Delta t+g(t)\sqrt{|\Delta t|}\mathbf{z}$$ 从 $t=T$（纯噪声）出发，逐步迭代到 $t=0$（生成样本），这就是最基本的采样流程。其中需要注意的一点是，由于随机项模拟的是维纳过程（Wiener Process/布朗运动），其方差与时间步长 $\Delta t$ 成正比，因此其标准差（即步长缩放）必须是 $\sqrt{|\Delta t|}$。

• 当然，Euler-Maruyama 方法的精度是一阶的，步长 $\Delta t$ 不能太大，否则误差会累积。实践中有几种改进策略：

  1. Predictor-Corrector 方法：Yang Song et al. [8] 在 Score SDE 论文中提出的经典方案。每一步先用数值 SDE 求解器（如 Euler-Maruyama）做一次"预测"（Predictor），然后再用若干步朗之万动力学做"校正"（Corrector），利用当前噪声尺度下的 $s_{\theta }(\mathbf{x},t)$ 对样本进行微调。这样可以在不增加太多计算量的前提下显著提升生成质量。
  2. 高阶 SDE 求解器：类比 ODE 中的 Runge-Kutta 方法，也可以构造 SDE 的高阶数值格式（如 Milstein 方法），但由于需要计算 score function 的导数，实际应用中并不常见。
  3. 自适应步长：根据当前 score 的变化幅度动态调整 $\Delta t$，在 score 变化剧烈的区域用小步长，平坦区域用大步长。

• 不过，SDE 求解器的一个本质局限在于：每一步都需要注入随机噪声 $\mathbf{z}$，这意味着同一个初始噪声 ${\mathbf{x}}_T$ 每次生成的结果都不同，而且随机性的引入也会增加方差，需要更多的步数才能收敛到高质量的样本。这就自然引出了一个问题——能不能去掉随机性，用确定性的方式来采样？ 答案是肯定的，这就是下一节要讨论的 Probability Flow ODE。

4.4. 从 SDE 到 Probability Flow ODE

• 这是 Score SDE 框架中一个极其优雅的结论：对于任意一个 Forward SDE，都存在一个对应的常微分方程（ODE），使得该 ODE 的解在每个时刻 $t$ 的边缘分布 $p_t(\mathbf{x})$ 与原 SDE 完全一致。这个 ODE 被称为 Probability Flow ODE，其形式为：$$\frac{d\mathbf{x}}{dt}=\mathbf{f}(\mathbf{x},t)-\frac{1}{2}{g}^2(t){\nabla }_{\mathbf{x}}\log p_t(\mathbf{x})$$ 对比 Reverse SDE：$$d\mathbf{x}=\left[\mathbf{f}(\mathbf{x},t)-g^2(t){\nabla }_{\mathbf{x}}\log p_t(\mathbf{x})\right]d\tau +g(t)d\bar{\mathbf{w}}$$ 可以发现，Probability Flow ODE 的漂移项恰好是 Reverse SDE 漂移项中 score 部分的"一半"，同时完全去掉了扩散项（随机噪声）。这个 $\frac{1}{2}$ 系数的来源可以这样理解：回顾 FPE 方程 $\frac{\partial p_t}{\partial t}=-\nabla \cdot (\mathbf{f}{p}_t)+\frac{1}{2}{g}^2(t){\nabla }^2{p}_t$，其中扩散项 $\frac{1}{2}{g}^2{\nabla }^2{p}_t$ 可以改写为 $-\nabla \cdot (-\frac{1}{2}{g}^2{p}_t\nabla \log p_t)$（利用 ${\nabla }^{2}p=\nabla \cdot (p\nabla \log p)$）。所以 FPE 等价于连续性方程 $\frac{\partial p_t}{\partial t}=-\nabla \cdot (\tilde{\mathbf{v}}{p}_t)$，其中 $\tilde{\mathbf{v}}=\mathbf{f}-\frac{1}{2}{g}^2\nabla \log p_t$，这正是 Probability Flow ODE 的向量场。换言之，SDE 中随机噪声项对概率质量的"扩散"效应，恰好等价于一个确定性的漂移 $-\frac{1}{2}{g}^2\nabla \log p_t$，所以去掉随机项后只需保留 score 的一半就能维持相同的边缘分布演化。

• 这个 ODE 的意义非常深远：

  1. 确定性采样：给定同一个初始噪声 ${\mathbf{x}}_T$，ODE 的解是唯一确定的，这意味着生成过程是完全可复现的。
  2. 精确的对数似然计算：利用连续正则化流（Continuous Normalizing Flow, CNF [10]）的理论，可以通过瞬时变量变换公式（Instantaneous Change of Variables）精确计算 $\log p_0(\mathbf{x})$：$$\log p_0(\mathbf{x}(0))=\log p_T(\mathbf{x}(T))+{\int }_0^T\nabla \cdot {\tilde{\mathbf{f}}}_{\theta }(\mathbf{x}(t),t)dt$$ 其中 ${\tilde{\mathbf{f}}}_{\theta }(\mathbf{x},t)=\mathbf{f}(\mathbf{x},t)-\frac{1}{2}{g}^2(t)s_{\theta }(\mathbf{x},t)$ 就是 Probability Flow ODE 的向量场。这在 SDE 框架下是做不到的。
  3. 可逆编码：由于 ODE 是可逆的，我们可以将任意数据点 ${\mathbf{x}}_0$ 沿正向 ODE 编码到隐空间 ${\mathbf{x}}_T$，再从 ${\mathbf{x}}_T$ 解码回来，实现无损的编码-解码。这为图像编辑、插值等下游任务提供了天然的支持。
  4. 高效求解器：ODE 的数值求解比 SDE 成熟得多，可以直接使用 Runge-Kutta、Adams-Bashforth 等经典方法，也可以利用自适应步长控制器（如 dopri5）来自动选择步长，在保证精度的同时最小化函数评估次数（NFE）。

• 至此，我们可以看到 SDE 框架的全貌：训练时，通过 Denoising Score Matching 学习 $s_{\theta }(\mathbf{x},t)$；采样时，既可以求解 Reverse SDE（随机采样，多样性更好），也可以求解 Probability Flow ODE（确定性采样，效率更高、支持似然计算）。两条路径殊途同归，共享同一个训练好的 score 模型。

---

Appendix

Appendix C：Reverse SDE 的局部高斯推导

1. Fokker-Planck 方程初始条件
   好的，下面我们来由标准的前向 Ito SDE 来推导逆向的随机微分方程，也就是 Reverse SDE：$$\begin{array}{rl}\text{Forward}:d\mathbf{x}& =\mathbf{f}(\mathbf{x},t)dt+g(t)d{\mathbf{w}}_t\\ \text{Reverse}:d\mathbf{x}& =\left[\mathbf{f}(\mathbf{x},t)-g^2(t){\nabla }_{\mathbf{x}}\log p_t(\mathbf{x})\right]d\tau +g(t)d\bar{\mathbf{w}}\end{array}$$ 我们想要得到上述结论，需要重要的一个结论就是 Fokker-Planck 方程：$$\frac{\partial p_t}{\partial t}=-\nabla \cdot (\mathbf{f}{p}_t)+\frac{1}{2}{g}^2(t){\nabla }^2{p}_t$$ 它反映了边缘概率分布的变化趋势，是构造流场的关键信息。那么我们首先需要构造出边缘概率分布的表达式，考虑一个光滑测试函数 $\varphi (\text{x})$（通常为紧支集，无穷可微），其中我们把 $\text{x}(t)$ 当成一个随着时间 $t$ 变化的随机变量，将其在 $p_t(\text{x})$ 的分布下求期望：$${\mathbb{E}}_{p_t(\text{x})}[\varphi (\text{x}(t))]={\int }_{{\mathbb{R}}^d}\varphi (\text{x})p_t(\text{x})d\text{x}$$ 将其对时间求导，可以得到：$$\frac{d}{dt}{\mathbb{E}}_{p_t(\text{x})}[\varphi (\text{x}(t))]={\int }_{{\mathbb{R}}^d}\varphi (\text{x})\frac{\partial p_t(\text{x})}{\partial t}d\text{x}$$ 下面我们主要考虑如何将 $\varphi (\text{x})$ 进行展开，由于 Ito 过程中均为微分项，故考虑 $d\varphi$ 展开如下（先考虑标量 Ito Equation）：$$\begin{array}{rl}d\varphi (x(t))& =\frac{d\varphi }{dx}dx+\frac{1}{2}\frac{d^2\varphi }{d{x}^2}(dx{)}^2+o(dx)\\ & =\frac{d\varphi }{dx}(f(x,t)dt+g(t)dw)+\frac{1}{2}\frac{d^2\varphi }{d{x}^2}(f(x,t)dt+g(t)dw{)}^2\end{array}$$ 其中需要注意的是，由于布朗运动时间增长的独立性：$$dw=w(t+\Delta t)-w(t)\sim \mathcal{N}(0,dt)$$ 也即在期望层面上有： $$(dx{)}^2=(g(t)dw{)}^2=g^2(t)dt,其中\mathbb{E}[(dw{)}^2]=dt$$ 此时，将标量 Ito 方程拓展至多维 Ito，满足方程：$$d(\varphi (\text{x}))={\nabla }_{\text{x}}\varphi \cdot d\text{x}+\frac{1}{2}\sum _{i,j=1}^d(\frac{{\partial }^2\varphi }{\partial x_i\partial y_i}{)}_{d\times d}(d\text{x}{)}_i(d\text{x}){)}_j$$

2. 多维 Ito Equation 展开
   其中先关注第一项，代入前向 SDE 后，由于 $\mathbb{E}(d\text{w})=0$，即有 $${\nabla }_{\text{x}}\varphi \cdot d\text{x}={\nabla }_{\text{x}}\varphi \cdot \text{f}dt+{\nabla }_{\text{x}}\varphi \cdot g(t)d\text{w}={\nabla }_{\text{x}}\varphi \cdot \text{f}dt$$ 至于后一项，则同样由布朗运动的性质，由于 $g^2(t)d{\text{w}}_{i}d{\text{w}}_j=g^2(t){\delta }_{ij}dt$，将第二项用拉普拉斯算子的形式展开，忽略所有的二阶无穷小量（如 $dtd\text{w}$ 这类）：$$\frac{1}{2}\sum _{i,j=1}^d(\frac{{\partial }^2\varphi }{\partial x_i\partial y_i}{)}_{d\times d}(d\text{x}{)}_i(d\text{x}){)}_j=\frac{1}{2}{g}^2(t){\nabla }_{\text{x}}^2\varphi dt$$ 所以综上所述，将 $d(\varphi (\text{x}))$ 对时间积分，得到 $$\varphi (\text{x}(t))-\varphi (\text{x}(0))={\int }_0^t({\nabla }_{\text{x}}\varphi \cdot \text{f}+\frac{1}{2}{g}^2(s){\nabla }_{\text{x}}^2\varphi )ds$$ 等式两边同时求期望后不难发现，由于我们的漂移项与扩散项光滑有界，由富比尼定理可以将积分与期望进行交换，即：$${\mathbb{E}}_{p_t}[\varphi (\text{x}(t))]={\int }_0^t({\mathbb{E}}_{p_t}[{\nabla }_{\text{x}}\varphi \cdot \text{f}+\frac{1}{2}{g}^2(s){\nabla }_{\text{x}}^2\varphi ])ds$$ 最后，我们再次求导并代入最开始的关于边缘概率分布的等式 $$\frac{d}{dt}{\mathbb{E}}_{p_t(\text{x})}[\varphi (\text{x}(t))]={\int }_{{\mathbb{R}}^d}\varphi (\mathbf{x})\frac{\partial p_t(\mathbf{x})}{\partial t}d\mathbf{x}={\int }_{{\mathbb{R}}^d}\left(\nabla \varphi (\mathbf{x})\cdot \mathbf{f}(\mathbf{x},t)+\frac{1}{2}{g}^2(t){\nabla }^2\varphi (\mathbf{x})\right)p_t(\mathbf{x})d\mathbf{x}$$

3. 分部积分与散度定理
   观察上面的等式，我们不难发现，既然测试函数是人为给定的，那么我们就需要想办法把 $\varphi$ 从散度算子中提出来。很自然地，我们就会想到使用散度定理（严谨定义不在此详细介绍）：$${\int }_{\Omega }{\nabla }_{\text{x}}\cdot \text{F}(\text{x})d\text{x}={\int }_{\partial \Omega }\text{F}(\text{x})\cdot \text{n}(\text{x})dS$$ 在此处，由于 $\varphi$ 是紧支集，${\int }_{\partial {\mathbb{R}}^d}\cdots dS=0$，故等式右边恒等于零，于是我们尝试用这个等式进行化简。

   ① 构造向量场 $\text{F}(\text{x})=\varphi (\text{x})\text{f}(\text{x},t)p_t(\text{x})$，那么分部积分得到 $${\nabla }_{\text{x}}\cdot \text{F}={\nabla }_{\text{x}}\varphi \cdot (\text{f}{p}_t)+\varphi {\nabla }_{\text{x}}\cdot (\text{f}{p}_t)$$ 于是有： $$0={\int }_{{\mathbb{R}}^d}{\nabla }_{\text{x}}\cdot \text{F}d\text{x}=\underset{\text{Part 1}}{\underbrace{{\int }_{{\mathbb{R}}^d}{\nabla }_{\text{x}}\varphi \cdot (\text{f}{p}_t)d\text{x}}}+{\int }_{{\mathbb{R}}^d}\varphi {\nabla }_{\text{x}}\cdot (\text{f}{p}_t)d\text{x}$$
   ② 同样，对于拉普拉斯算子，需要做两次分部积分，令 $\text{G}(\text{x})=p_t(\text{x}){\nabla }_{\text{x}}\varphi (\text{x}),\text{H}(\text{x})=\varphi (\text{x}){\nabla }_{\text{x}}{p}_t(\text{x})$，则分部积分得到：$$\begin{array}{r}{\nabla }_{\text{x}}\cdot \text{G}(\text{x})={\nabla }_{\text{x}}{p}_t(\text{x}){\nabla }_{\text{x}}\varphi (\text{x})+p_t(\text{x}){\nabla }_{\text{x}}^2\varphi (\text{x})\\ {\nabla }_{\text{x}}\cdot \text{H}(\text{x})={\nabla }_{\text{x}}\varphi (\text{x}){\nabla }_{\text{x}}{p}_t(\text{x})+\varphi (\text{x}){\nabla }_{\text{x}}^2{p}_t(\text{x})\end{array}$$ 所以代入散度定理的公式：$$\begin{array}{r}0=\underset{\text{Part 2}}{\underbrace{{\int }_{{\mathbb{R}}^d}{p}_t{\nabla }_{\text{x}}^2\varphi d\text{x}}}+{\int }_{{\mathbb{R}}^d}{\nabla }_{\text{x}}{p}_t\cdot {\nabla }_{\text{x}}\varphi d\text{x}\\ 0={\int }_{{\mathbb{R}}^d}{\nabla }_{\text{x}}{p}_t\cdot {\nabla }_{\text{x}}\varphi d\text{x}+{\int }_{{\mathbb{R}}^d}\varphi {\nabla }_{\text{x}}^2{p}_{t}d\text{x}\end{array}$$ 最后，将上述结果，代入期望演化方程，由于测试函数的任意性，被积函数也就相应的必须相等，即可得到 FPE 方程：$$\frac{\partial p_t}{\partial t}=-\nabla \cdot (\mathbf{f}{p}_t)+\frac{1}{2}{g}^2(t){\nabla }^2{p}_t$$

4. 化简得到 Reverse SDE
   首先考虑拉普拉斯算子的另一种等效形式：$$\nabla \cdot (p\nabla \log p)=\nabla \cdot (p\frac{\nabla p}{p})=\nabla \cdot (\nabla p)={\nabla }^{2}p$$ 将其代入 FPE 后合并散度算子，即可得到：$$\frac{\partial p_t}{\partial t}=-\nabla \cdot \left[\mathbf{f}{p}_t-\frac{1}{2}{g}^2(t)p_t\nabla \log p_t\right]$$ 假设逆向 SDE 的形式为（注意：由于时间反向，符号与正向 FPE 相反）：$$d\mathbf{x}={\mathbf{f}}_{\text{rev}}(\mathbf{x},t)dt+g(t)d\bar{\mathbf{w}}$$ 这里的时间 $t$ 是从 $T$ 回退到 $0$。需要注意的是，逆向 SDE 中的扩散系数 $g(t)$ 与正向相同且符号不变，这是因为扩散项对应 FPE 中的二阶导数项 $\frac{1}{2}{g}^2{\nabla }^2{p}_t$，而二阶空间导数在时间反转下不改变符号（它只描述概率质量的局部扩散速率，与时间方向无关），因此逆向过程的随机扰动强度与正向一致。在物理上，逆向演化的概率密度演化满足：$$\frac{\partial p_t}{\partial t}=+\nabla \cdot {\mathbf{J}}_{\text{rev}}$$ 为了使正向和逆向过程描述同一个密度演化过程，必须满足：$${\mathbf{J}}_{\text{rev}}=-{\mathbf{J}}_{\text{fwd}},\text{a.k.a.}{\mathbf{f}}_{\text{rev}}{p}_t-\frac{1}{2}{g}^2(t)p_t\nabla \log p_t=-\left(\mathbf{f}{p}_t-\frac{1}{2}{g}^2(t)p_t\nabla \log p_t\right)$$ 消去 $p_t$ 并解出 ${\mathbf{f}}_{\text{rev}}$：$${\mathbf{f}}_{\text{rev}}=-\mathbf{f}+g^2\nabla \log p_t$$
   于是此时再代入，即可得到最终版本的 Reverse SDE：$$d\mathbf{x}=\left[\mathbf{f}(\mathbf{x},t)-g^2(t){\nabla }_{\mathbf{x}}\log p_t(\mathbf{x})\right]d\tau +g(t)d\bar{\mathbf{w}}$$

---

Reference

[8] Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., & Poole, B. (2021). Score-Based Generative Modeling through Stochastic Differential Equations. ICLR 2021.

[9] Anderson, B. D. O. (1982). Reverse-Time Diffusion Equation Models. Stochastic Processes and their Applications, 12(3), 313–326.

[10] Chen, R. T. Q., Rubanova, Y., Bettencourt, J., & Duvenaud, D. (2018). Neural Ordinary Differential Equations. NeurIPS 2018.