Neural Ordinary Differential Equations

神经常微分方程（Neural Ordinary Differential Equations）

Ricky T. Q. Chen*, Yulia Rubanova*, Jesse Bettencourt*, David Duvenaud

University of Toronto, Vector Institute

{rtqichen, rubanova, jessebett, duvenaud}@cs.toronto.edu

Ricky T. Q. Chen*、Yulia Rubanova*、Jesse Bettencourt*、David Duvenaud

多伦多大学，Vector 研究所

{rtqichen, rubanova, jessebett, duvenaud}@cs.toronto.edu

Abstract

摘要

We introduce a new family of deep neural network models. Instead of specifying a discrete sequence of hidden layers, we parameterize the derivative of the hidden state using a neural network. The output of the network is computed using a black-box differential equation solver. These continuous-depth models have constant memory cost, adapt their evaluation strategy to each input, and can explicitly trade numerical precision for speed. We demonstrate these properties in continuous-depth residual networks and continuous-time latent variable models. We also construct continuous normalizing flows, a generative model that can train by maximum likelihood, without partitioning or ordering the data dimensions. For training, we show how to scalably backpropagate through any ODE solver, without access to its internal operations. This allows end-to-end training of ODEs within larger models.

我们提出一类新的深度神经网络模型。它不再指定一个离散的隐藏层序列，而是用神经网络参数化隐藏状态的导数。网络的输出由一个黑箱微分方程求解器计算。这些连续深度（continuous-depth）模型的内存开销恒定，能针对每个输入调整求值策略，并且可以显式地用数值精度换取速度。我们在连续深度残差网络和连续时间隐变量模型中展示了这些性质。我们还构建了连续归一化流（continuous normalizing flows），这是一种可以用最大似然训练的生成模型，无需对数据维度做划分或排序。在训练方面，我们展示了如何在不访问求解器内部运算的情况下，以可扩展的方式对任意常微分方程（ODE）求解器做反向传播。这使得 ODE 可以在更大的模型中端到端训练。

Residual Network	ODE Network

Figure 1: Left: A Residual network defines a discrete sequence of finite transformations. Right: A ODE network defines a vector field, which continuously transforms the state. Both: Circles represent evaluation locations.

残差网络	ODE 网络

图 1：*左：*残差网络定义了一个由有限变换组成的离散序列。*右：*ODE 网络定义了一个向量场，它连续地变换状态。*两者：*圆圈表示求值位置。

1 Introduction

1 引言

Models such as residual networks, recurrent neural network decoders, and normalizing flows build complicated transformations by composing a sequence of transformations to a hidden state:

残差网络、循环神经网络解码器、归一化流（normalizing flows）等模型，都是通过对隐藏状态依次施加一串变换来构建复杂变换的：

$$\mathbf{h}_{t+1}=\mathbf{h}_{t}+f(\mathbf{h}_{t},\theta_{t}) \tag{1}$$

where $t\in\{0\dots T\}$ and $\mathbf{h}_{t}\in\mathbb{R}^{D}$. These iterative updates can be seen as an Euler discretization of a continuous transformation (Lu et al. 2017; Haber and Ruthotto 2017; Ruthotto and Haber 2018).

其中 $t\in\{0\dots T\}$，$\mathbf{h}_{t}\in\mathbb{R}^{D}$。这些迭代更新可以看作对某个连续变换的欧拉离散化（Lu et al. 2017; Haber and Ruthotto 2017; Ruthotto and Haber 2018）。

What happens as we add more layers and take smaller steps? In the limit, we parameterize the continuous dynamics of hidden units using an ordinary differential equation (ODE) specified by a neural network:

如果我们加入更多层、步长取得更小，会发生什么？在极限情况下，我们用一个由神经网络指定的常微分方程（ordinary differential equation，ODE）来参数化隐藏单元的连续动力学：

$$\frac{d\mathbf{h}(t)}{dt}=f(\mathbf{h}(t),t,\theta) \tag{2}$$

Starting from the input layer $\mathbf{h}(0)$, we can define the output layer $\mathbf{h}(T)$ to be the solution to this ODE initial value problem at some time $T$. This value can be computed by a black-box differential equation solver, which evaluates the hidden unit dynamics $f$ wherever necessary to determine the solution with the desired accuracy. Figure 1 contrasts these two approaches.

从输入层 $\mathbf{h}(0)$ 出发，可以把输出层 $\mathbf{h}(T)$ 定义为这个 ODE 初值问题在某个时刻 $T$ 的解。这个值可以由黑箱微分方程求解器计算，求解器会在需要的地方对隐藏单元动力学 $f$ 求值，以达到所要求的精度。图 1 对比了这两种方法。

Defining and evaluating models using ODE solvers has several benefits:

用 ODE 求解器定义和计算模型有几项好处：

Memory efficiency

In Section 2, we show how to compute gradients of a scalar-valued loss with respect to all inputs of any ODE solver, without backpropagating through the operations of the solver. Not storing any intermediate quantities of the forward pass allows us to train our models with constant memory cost as a function of depth, a major bottleneck of training deep models.

内存效率

在第 2 节中，我们展示如何计算标量损失关于任意 ODE 求解器所有输入的梯度，而不需要对求解器内部的运算做反向传播。由于不存储前向传播的任何中间量，我们训练模型的内存开销相对深度是恒定的；而内存正是训练深层模型的一个主要瓶颈。

Adaptive computation

Euler’s method is perhaps the simplest method for solving ODEs. There have since been more than 120 years of development of efficient and accurate ODE solvers (Runge 1895; Kutta 1901; Hairer et al. 1987). Modern ODE solvers provide guarantees about the growth of approximation error, monitor the level of error, and adapt their evaluation strategy on the fly to achieve the requested level of accuracy. This allows the cost of evaluating a model to scale with problem complexity. After training, accuracy can be reduced for real-time or low-power applications.

自适应计算

欧拉法大概是求解 ODE 最简单的方法。此后，高效而精确的 ODE 求解器已经发展了 120 多年（Runge 1895; Kutta 1901; Hairer et al. 1987）。现代 ODE 求解器能对近似误差的增长给出保证，监控误差水平，并即时调整求值策略以达到所要求的精度。这让模型的求值成本可以随问题的复杂度伸缩。训练完成后，还可以为实时或低功耗应用降低精度。

Scalable and invertible normalizing flows

An unexpected side-benefit of continuous transformations is that the change of variables formula becomes easier to compute. In Section 4, we derive this result and use it to construct a new class of invertible density models that avoids the single-unit bottleneck of normalizing flows, and can be trained directly by maximum likelihood.

可扩展、可逆的归一化流

连续变换带来一个意外的附带好处：变量替换公式（change of variables formula）变得更容易计算。在第 4 节中，我们推导这一结果，并用它构建一类新的可逆密度模型。这类模型避免了归一化流的单单元瓶颈，并可以直接用最大似然训练。

Continuous time-series models

Unlike recurrent neural networks, which require discretizing observation and emission intervals, continuously-defined dynamics can naturally incorporate data which arrives at arbitrary times. In Section 5, we construct and demonstrate such a model.

连续时间序列模型

循环神经网络需要把观测和输出的时间间隔离散化；而连续定义的动力学可以自然地纳入在任意时刻到达的数据。在第 5 节中，我们构建并演示了这样一个模型。

2 Reverse-mode automatic differentiation of ODE solutions

2 ODE 解的反向模式自动微分

The main technical difficulty in training continuous-depth networks is performing reverse-mode differentiation (also known as backpropagation) through the ODE solver. Differentiating through the operations of the forward pass is straightforward, but incurs a high memory cost and introduces additional numerical error.

训练连续深度网络的主要技术难点，是穿过 ODE 求解器做反向模式微分（也就是反向传播）。直接对前向传播中的运算求导很简单，但内存开销高，还会引入额外的数值误差。

We treat the ODE solver as a black box, and compute gradients using the adjoint sensitivity method (Pontryagin et al. 1962). This approach computes gradients by solving a second, augmented ODE backwards in time, and is applicable to all ODE solvers. This approach scales linearly with problem size, has low memory cost, and explicitly controls numerical error.

我们把 ODE 求解器当作黑箱，用伴随灵敏度方法（adjoint sensitivity method）（Pontryagin et al. 1962）计算梯度。这种方法通过沿时间反向求解第二个增广 ODE 来计算梯度，适用于所有 ODE 求解器。它的计算量随问题规模线性增长，内存开销低，并能显式控制数值误差。

Consider optimizing a scalar-valued loss function $L()$, whose input is the result of an ODE solver:

考虑优化一个标量损失函数 $L()$，它的输入是 ODE 求解器的结果：

$$L(\mathbf{z}({t_{\textnormal{1}}}))=L\left(\mathbf{z}({t_{\textnormal{0}}})+\int_{{t_{\textnormal{0}}}}^{{t_{\textnormal{1}}}}f(\mathbf{z}(t),t,\theta)dt\right)=L\left(\textnormal{ODESolve}(\mathbf{z}({t_{\textnormal{0}}}),f,{t_{\textnormal{0}}},{t_{\textnormal{1}}},\theta)\right) \tag{3}$$

Figure 2: Reverse-mode differentiation of an ODE solution. The adjoint sensitivity method solves an augmented ODE backwards in time. The augmented system contains both the original state and the sensitivity of the loss with respect to the state. If the loss depends directly on the state at multiple observation times, the adjoint state must be updated in the direction of the partial derivative of the loss with respect to each observation.

图 2：ODE 解的反向模式微分。伴随灵敏度方法沿时间反向求解一个增广 ODE。增广系统同时包含原始状态，以及损失对状态的灵敏度。如果损失直接依赖于多个观测时刻的状态，伴随状态就必须沿损失对每个观测的偏导数方向更新。

To optimize $L$, we require gradients with respect to $\theta$. The first step is to determining how the gradient of the loss depends on the hidden state $\mathbf{z}(t)$ at each instant. This quantity is called the adjoint $\mathbf{a}(t)={{\partial{L}}}/{{\partial\mathbf{z}(t)}}$. Its dynamics are given by another ODE, which can be thought of as the instantaneous analog of the chain rule:

为了优化 $L$，我们需要关于 $\theta$ 的梯度。第一步是确定损失的梯度在每一时刻如何依赖于隐藏状态 $\mathbf{z}(t)$。这个量称为伴随（adjoint）$\mathbf{a}(t)=\partial L/\partial\mathbf{z}(t)$。它的动力学由另一个 ODE 给出，可以看作链式法则的瞬时版本：

$$\frac{d\mathbf{a}(t)}{dt}=-\mathbf{a}(t)^{\mkern-1.5mu\mathsf{T}}{}\frac{\partial f(\mathbf{z}(t),t,\theta)}{\partial\mathbf{z}} \tag{4}$$

We can compute ${{\partial L}}/{{\partial\mathbf{z}({t_{\textnormal{0}}})}}$ by another call to an ODE solver. This solver must run backwards, starting from the initial value of ${{\partial L}}/{{\partial\mathbf{z}({t_{\textnormal{1}}})}}$. One complication is that solving this ODE requires the knowing value of $\mathbf{z}(t)$ along its entire trajectory. However, we can simply recompute $\mathbf{z}(t)$ backwards in time together with the adjoint, starting from its final value $\mathbf{z}({t_{\textnormal{1}}})$.

我们可以再调用一次 ODE 求解器来计算 $\partial L/\partial\mathbf{z}(t_{0})$。这个求解器必须反向运行，从初值 $\partial L/\partial\mathbf{z}(t_{1})$ 出发。一个麻烦之处在于，求解这个 ODE 需要知道 $\mathbf{z}(t)$ 在整条轨迹上的值。不过，我们可以从 $\mathbf{z}(t)$ 的终值 $\mathbf{z}(t_{1})$ 出发，把它与伴随一起沿时间反向重新算出来。

Computing the gradients with respect to the parameters $\theta$ requires evaluating a third integral, which depends on both $\mathbf{z}(t)$ and $\mathbf{a}(t)$:

计算关于参数 $\theta$ 的梯度，需要再求第三个积分，它同时依赖于 $\mathbf{z}(t)$ 和 $\mathbf{a}(t)$：

$$\frac{dL}{d\theta}=-\int^{t_{\textnormal{0}}}_{t_{\textnormal{1}}}\mathbf{a}(t)^{\mkern-1.5mu\mathsf{T}}{}\frac{\partial f(\mathbf{z}(t),t,\theta)}{\partial\theta}dt \tag{5}$$

The vector-Jacobian products $\mathbf{a}(t)^{T}\frac{\partial f}{\partial\mathbf{z}}$ and $\mathbf{a}(t)^{T}\frac{\partial f}{\partial\theta}$ in (4) and (5) can be efficiently evaluated by automatic differentiation, at a time cost similar to that of evaluating $f$. All integrals for solving $\mathbf{z}$, $\mathbf{a}$ and $\frac{\partial L}{\partial\theta}$ can be computed in a single call to an ODE solver, which concatenates the original state, the adjoint, and the other partial derivatives into a single vector. Algorithm 1 shows how to construct the necessary dynamics, and call an ODE solver to compute all gradients at once.

式 (4) 和式 (5) 中的向量-雅可比积（vector-Jacobian product）$\mathbf{a}(t)^{T}\frac{\partial f}{\partial\mathbf{z}}$ 和 $\mathbf{a}(t)^{T}\frac{\partial f}{\partial\theta}$ 可以用自动微分高效计算，时间开销与计算 $f$ 相当。求解 $\mathbf{z}$、$\mathbf{a}$ 和 $\frac{\partial L}{\partial\theta}$ 所需的所有积分，都可以通过一次 ODE 求解器调用完成：把原始状态、伴随和其他偏导数拼接成一个向量即可。算法 1 展示了如何构造所需的动力学，并调用 ODE 求解器一次算出所有梯度。

Algorithm 1 Reverse-mode derivative of an ODE initial value problem

Input: dynamics parameters $\theta$, start time $t_{0}$, stop time $t_{1}$, final state $\mathbf{z}(t_{1})$, loss gradient $\partial L/\partial\mathbf{z}(t_{1})$

1. $s_{0}=[\mathbf{z}(t_{1}),\frac{\partial L}{\partial\mathbf{z}(t_{1})},\mathbf{0}_{|\theta|}]$ $\triangleright$ Define initial augmented state
2. def aug_dynamics($[\mathbf{z}(t),\mathbf{a}(t),\cdot],t,\theta$): $\triangleright$ Define dynamics on augmented state
3. 　return $[f(\mathbf{z}(t),t,\theta),-\mathbf{a}(t)^{\mathsf{T}}\frac{\partial f}{\partial\mathbf{z}},-\mathbf{a}(t)^{\mathsf{T}}\frac{\partial f}{\partial\theta}]$ $\triangleright$ Compute vector-Jacobian products
4. $[\mathbf{z}(t_{0}),\frac{\partial L}{\partial\mathbf{z}(t_{0})},\frac{\partial L}{\partial\theta}]=\textnormal{ODESolve}(s_{0},\textnormal{aug\_dynamics},t_{1},t_{0},\theta)$ $\triangleright$ Solve reverse-time ODE
5. return $\frac{\partial L}{\partial\mathbf{z}(t_{0})},\frac{\partial L}{\partial\theta}$ $\triangleright$ Return gradients

算法 1 ODE 初值问题的反向模式导数

输入：动力学参数 $\theta$，起始时间 $t_{0}$，终止时间 $t_{1}$，终态 $\mathbf{z}(t_{1})$，损失梯度 $\partial L/\partial\mathbf{z}(t_{1})$

1. $s_{0}=[\mathbf{z}(t_{1}),\frac{\partial L}{\partial\mathbf{z}(t_{1})},\mathbf{0}_{|\theta|}]$ $\triangleright$ 定义初始增广状态
2. def aug_dynamics($[\mathbf{z}(t),\mathbf{a}(t),\cdot],t,\theta$): $\triangleright$ 定义增广状态上的动力学
3. 　return $[f(\mathbf{z}(t),t,\theta),-\mathbf{a}(t)^{\mathsf{T}}\frac{\partial f}{\partial\mathbf{z}},-\mathbf{a}(t)^{\mathsf{T}}\frac{\partial f}{\partial\theta}]$ $\triangleright$ 计算向量-雅可比积
4. $[\mathbf{z}(t_{0}),\frac{\partial L}{\partial\mathbf{z}(t_{0})},\frac{\partial L}{\partial\theta}]=\textnormal{ODESolve}(s_{0},\textnormal{aug\_dynamics},t_{1},t_{0},\theta)$ $\triangleright$ 求解反向时间 ODE
5. return $\frac{\partial L}{\partial\mathbf{z}(t_{0})},\frac{\partial L}{\partial\theta}$ $\triangleright$ 返回梯度

Most ODE solvers have the option to output the state $\mathbf{z}(t)$ at multiple times. When the loss depends on these intermediate states, the reverse-mode derivative must be broken into a sequence of separate solves, one between each consecutive pair of output times (Figure 2). At each observation, the adjoint must be adjusted in the direction of the corresponding partial derivative ${{\partial L}}/{{\partial\mathbf{z}(t_{i})}}$.

大多数 ODE 求解器都可以在多个时刻输出状态 $\mathbf{z}(t)$。当损失依赖于这些中间状态时，反向模式导数必须拆成一连串单独的求解，每对相邻输出时刻之间各求解一次（图 2）。在每个观测处，伴随都要沿相应偏导数 $\partial L/\partial\mathbf{z}(t_{i})$ 的方向调整。

The results above extend those of Stapor et al. 2018. An extended version of Algorithm 1 including derivatives w.r.t. ${t_{\textnormal{0}}}$ and ${t_{\textnormal{1}}}$ can be found in Appendix C. Detailed derivations are provided in Appendix B. Appendix D provides Python code which computes all derivatives for scipy.integrate.odeint by extending the autograd automatic differentiation package. This code also supports all higher-order derivatives. We have since released a PyTorch (Paszke et al. 2017) implementation, including GPU-based implementations of several standard ODE solvers at github.com/rtqichen/torchdiffeq.

上述结果扩展了 Stapor et al. 2018 的结果。附录 C 给出了算法 1 的扩展版本，其中包含关于 $t_{0}$ 和 $t_{1}$ 的导数。附录 B 给出了详细推导。附录 D 提供了 Python 代码，它通过扩展 autograd 自动微分包，为 scipy.integrate.odeint 计算所有导数。这段代码也支持所有高阶导数。此后我们还发布了 PyTorch（Paszke et al. 2017）实现，其中包括几种标准 ODE 求解器的 GPU 实现，地址为 github.com/rtqichen/torchdiffeq。

3 Replacing residual networks with ODEs for supervised learning

3 在监督学习中用 ODE 替代残差网络

In this section, we experimentally investigate the training of neural ODEs for supervised learning.

本节通过实验考察神经 ODE 在监督学习中的训练。

Software

To solve ODE initial value problems numerically, we use the implicit Adams method implemented in LSODE and VODE and interfaced through the scipy.integrate package. Being an implicit method, it has better guarantees than explicit methods such as Runge-Kutta but requires solving a nonlinear optimization problem at every step. This setup makes direct backpropagation through the integrator difficult. We implement the adjoint sensitivity method in Python’s autograd framework (Maclaurin et al. 2015). For the experiments in this section, we evaluated the hidden state dynamics and their derivatives on the GPU using Tensorflow, which were then called from the Fortran ODE solvers, which were called from Python autograd code.

软件

为了数值求解 ODE 初值问题，我们使用 LSODE 和 VODE 中实现、通过 scipy.integrate 包调用的隐式 Adams 方法。作为隐式方法，它比龙格-库塔（Runge-Kutta）等显式方法有更好的保证，但每一步都需要求解一个非线性优化问题。这种设置使得直接穿过积分器做反向传播很困难。我们在 Python 的 autograd 框架（Maclaurin et al. 2015）中实现了伴随灵敏度方法。本节的实验中，隐藏状态动力学及其导数在 GPU 上用 Tensorflow 计算，由 Fortran 编写的 ODE 求解器调用，而这些求解器又由 Python autograd 代码调用。

Model Architectures

	Test Error	# Params	Memory	Time
1-Layer MLP†	1.60%	0.24 M	-	-
ResNet	0.41%	0.60 M	$\mathcal{O}(L)$	$\mathcal{O}(L)$
RK-Net	0.47%	0.22 M	$\mathcal{O}(\tilde{L})$	$\mathcal{O}(\tilde{L})$
ODE-Net	0.42%	0.22 M	$\boldsymbol{\mathcal{O}(1)}$	$\mathcal{O}(\tilde{L})$

Table 1: Performance on MNIST. †From LeCun et al. 1998.

模型架构

	测试误差	参数量	内存	时间
单层 MLP†	1.60%	0.24 M	-	-
ResNet	0.41%	0.60 M	$\mathcal{O}(L)$	$\mathcal{O}(L)$
RK-Net	0.47%	0.22 M	$\mathcal{O}(\tilde{L})$	$\mathcal{O}(\tilde{L})$
ODE-Net	0.42%	0.22 M	$\boldsymbol{\mathcal{O}(1)}$	$\mathcal{O}(\tilde{L})$

表 1：MNIST 上的性能。† 取自 LeCun et al. 1998。

We experiment with a small residual network which downsamples the input twice then applies 6 standard residual blocks He et al. 2016b, which are replaced by an ODESolve module in the ODE-Net variant. We also test a network with the same architecture but where gradients are backpropagated directly through a Runge-Kutta integrator, referred to as RK-Net. Table 1 shows test error, number of parameters, and memory cost. $L$ denotes the number of layers in the ResNet, and $\tilde{L}$ is the number of function evaluations that the ODE solver requests in a single forward pass, which can be interpreted as an implicit number of layers. We find that ODE-Nets and RK-Nets can achieve around the same performance as the ResNet.

我们实验了一个小型残差网络：它先对输入做两次下采样，再接 6 个标准残差块（He et al. 2016b）；在 ODE-Net 变体中，这些残差块被一个 ODESolve 模块替代。我们还测试了一个架构相同、但梯度直接穿过龙格-库塔积分器反向传播的网络，称为 RK-Net。表 1 给出了测试误差、参数量和内存开销。$L$ 表示 ResNet 的层数，$\tilde{L}$ 是 ODE 求解器在一次前向传播中请求的函数求值次数，可以理解为隐式的层数。我们发现，ODE-Net 和 RK-Net 的性能与 ResNet 大致相当。

Error Control in ODE-Nets

ODE solvers can approximately ensure that the output is within a given tolerance of the true solution. Changing this tolerance changes the behavior of the network. We first verify that error can indeed be controlled in Figure 3a. The time spent by the forward call is proportional to the number of function evaluations (Figure 3b), so tuning the tolerance gives us a trade-off between accuracy and computational cost. One could train with high accuracy, but switch to a lower accuracy at test time.

ODE-Net 中的误差控制

ODE 求解器可以近似保证输出与真实解的误差在给定容差（tolerance）以内。改变这个容差会改变网络的行为。我们先在图 3a 中验证误差确实可以控制。前向调用花费的时间与函数求值次数成正比（图 3b），因此调节容差就能在精度和计算成本之间权衡。可以用高精度训练，测试时再切换到较低精度。

Figure 3: Statistics of a trained ODE-Net. (NFE = number of function evaluations.)

图 3：一个训练好的 ODE-Net 的统计量。（NFE = 函数求值次数。）

Figure 3c) shows a surprising result: the number of evaluations in the backward pass is roughly half of the forward pass. This suggests that the adjoint sensitivity method is not only more memory efficient, but also more computationally efficient than directly backpropagating through the integrator, because the latter approach will need to backprop through each function evaluation in the forward pass.

图 3c 给出了一个出人意料的结果：反向传播中的求值次数大约是前向传播的一半。这说明伴随灵敏度方法不仅比直接穿过积分器反向传播更省内存，计算上也更高效，因为后者需要对前向传播中的每一次函数求值做反向传播。

Network Depth

It’s not clear how to define the ‘depth‘ of an ODE solution. A related quantity is the number of evaluations of the hidden state dynamics required, a detail delegated to the ODE solver and dependent on the initial state or input. Figure 3d shows that he number of function evaluations increases throughout training, presumably adapting to increasing complexity of the model.

网络深度

如何定义 ODE 解的“深度”并不清楚。一个相关的量是所需的隐藏状态动力学求值次数，这个细节交给 ODE 求解器决定，并取决于初始状态或输入。图 3d 显示，函数求值次数在整个训练过程中不断增加，这大概是在适应模型不断增加的复杂度。

4 Continuous Normalizing Flows

4 连续归一化流

The discretized equation (1) also appears in normalizing flows (Rezende and Mohamed 2015) and the NICE framework (Dinh et al. 2014). These methods use the change of variables theorem to compute exact changes in probability if samples are transformed through a bijective function $f$:

离散化的方程 (1) 也出现在归一化流（Rezende and Mohamed 2015）和 NICE 框架（Dinh et al. 2014）中。这些方法利用变量替换定理，在样本经过双射函数 $f$ 变换时精确计算概率的变化：

$${\mathbf{z}}_{1}=f({\mathbf{z}}_{0})\;\implies\log p({\mathbf{z}}_{1})=\log p({\mathbf{z}}_{0})-\log\left|\det\frac{\partial f}{\partial{\mathbf{z}}_{0}}\right| \tag{6}$$

An example is the planar normalizing flow (Rezende and Mohamed 2015):

平面归一化流（planar normalizing flow）（Rezende and Mohamed 2015）就是一个例子：

$${\mathbf{z}}(t+1)={\mathbf{z}}(t)+uh(w^{\mkern-1.5mu\mathsf{T}}{}{\mathbf{z}}(t)+b),\quad\log p({\mathbf{z}}(t+1))=\log p({\mathbf{z}}(t))-\log\left|1+u^{\mkern-1.5mu\mathsf{T}}{}\frac{\partial h}{\partial{\mathbf{z}}}\right| \tag{7}$$

Generally, the main bottleneck to using the change of variables formula is computing of the determinant of the Jacobian ${{\partial f}}/{{\partial{\mathbf{z}}}}$, which has a cubic cost in either the dimension of $\mathbf{z}$, or the number of hidden units. Recent work explores the tradeoff between the expressiveness of normalizing flow layers and computational cost (Kingma et al. 2016; Tomczak and Welling 2016; Berg et al. 2018).

一般来说，使用变量替换公式的主要瓶颈是计算雅可比矩阵 $\partial f/\partial\mathbf{z}$ 的行列式，其开销是 $\mathbf{z}$ 的维度或隐藏单元数的三次方。近来的工作探索了归一化流层的表达能力与计算成本之间的权衡（Kingma et al. 2016; Tomczak and Welling 2016; Berg et al. 2018）。

Surprisingly, moving from a discrete set of layers to a continuous transformation simplifies the computation of the change in normalizing constant:

令人意外的是，从一组离散的层转向连续变换，反而简化了归一化常数变化的计算：

Theorem 1 (Instantaneous Change of Variables).

Let ${\mathbf{z}}(t)$ be a finite continuous random variable with probability $p({\mathbf{z}}(t))$ dependent on time. Let $\frac{d{\mathbf{z}}}{dt}=f({\mathbf{z}}(t),t)$ be a differential equation describing a continuous-in-time transformation of ${\mathbf{z}}(t)$. Assuming that $f$ is uniformly Lipschitz continuous in ${\mathbf{z}}$ and continuous in $t$, then the change in log probability also follows a differential equation,

定理 1（瞬时变量替换，Instantaneous Change of Variables）

设 $\mathbf{z}(t)$ 是一个有限的连续随机变量，其概率 $p(\mathbf{z}(t))$ 随时间变化。设 $\frac{d\mathbf{z}}{dt}=f(\mathbf{z}(t),t)$ 是描述 $\mathbf{z}(t)$ 随时间连续变换的微分方程。假设 $f$ 关于 $\mathbf{z}$ 一致利普希茨连续（uniformly Lipschitz continuous），关于 $t$ 连续，那么对数概率的变化也服从一个微分方程：

$$\frac{\partial\log p({\mathbf{z}}(t))}{\partial t}=-\textnormal{tr}\left(\frac{df}{d{\mathbf{z}}(t)}\right) \tag{8}$$

Proof in Appendix A. Instead of the log determinant in (6), we now only require a trace operation. Also unlike standard finite flows, the differential equation $f$ does not need to be bijective, since if uniqueness is satisfied, then the entire transformation is automatically bijective.

证明见附录 A。我们现在只需要一次求迹运算，而不是式 (6) 中的对数行列式。此外，与标准的有限流不同，微分方程 $f$ 不必是双射，因为只要满足解的唯一性，整个变换自动就是双射。

As an example application of the instantaneous change of variables, we can examine the continuous analog of the planar flow, and its change in normalization constant:

作为瞬时变量替换的一个应用示例，我们可以考察平面流的连续版本，以及它的归一化常数变化：

$$\frac{d{\mathbf{z}}(t)}{dt}=uh(w^{\mkern-1.5mu\mathsf{T}}{}{\mathbf{z}}(t)+b),\quad\frac{\partial\log p({\mathbf{z}}(t))}{\partial t}=-u^{\mkern-1.5mu\mathsf{T}}{}\frac{\partial h}{\partial{\mathbf{z}}(t)} \tag{9}$$

Given an initial distribution $p({\mathbf{z}}(0))$, we can sample from $p({\mathbf{z}}(t))$ and evaluate its density by solving this combined ODE.

给定初始分布 $p(\mathbf{z}(0))$，求解这个联合 ODE，就可以从 $p(\mathbf{z}(t))$ 中采样并计算其密度。

Using multiple hidden units with linear cost

While $\det$ is not a linear function, the trace function is, which implies $\textnormal{tr}(\sum_{n}J_{n})=\sum_{n}\textnormal{tr}(J_{n})$. Thus if our dynamics is given by a sum of functions then the differential equation for the log density is also a sum:

以线性开销使用多个隐藏单元

$\det$ 不是线性函数，但迹是线性的，即 $\textnormal{tr}(\sum_{n}J_{n})=\sum_{n}\textnormal{tr}(J_{n})$。因此，如果动力学由若干函数之和给出，那么对数密度的微分方程也是一个求和：

$$\frac{d{\mathbf{z}}(t)}{dt}=\sum_{n=1}^{M}f_{n}({\mathbf{z}}(t)),\quad\frac{d\log p({\mathbf{z}}(t))}{dt}=\sum_{n=1}^{M}\textnormal{tr}\left(\frac{\partial f_{n}}{\partial{\mathbf{z}}}\right) \tag{10}$$

This means we can cheaply evaluate flow models having many hidden units, with a cost only linear in the number of hidden units $M$. Evaluating such ‘wide’ flow layers using standard normalizing flows costs $\mathcal{O}(M^{3})$, meaning that standard NF architectures use many layers of only a single hidden unit.

这意味着我们可以低成本地计算具有很多隐藏单元的流模型，开销只与隐藏单元数 $M$ 呈线性关系。用标准归一化流计算这种“宽”流层的开销是 $\mathcal{O}(M^{3})$，所以标准归一化流架构使用很多层、每层只有一个隐藏单元。

Time-dependent dynamics

We can specify the parameters of a flow as a function of $t$, making the differential equation $f({\mathbf{z}}(t),t)$ change with $t$. This is parameterization is a kind of hypernetwork (Ha et al. 2016). We also introduce a gating mechanism for each hidden unit, $\frac{d{\mathbf{z}}}{dt}=\sum_{n}\sigma_{n}(t)f_{n}({\mathbf{z}})$ where $\sigma_{n}(t)\in(0,1)$ is a neural network that learns when the dynamic $f_{n}({\mathbf{z}})$ should be applied. We call these models continuous normalizing flows (CNF).

随时间变化的动力学

我们可以把流的参数设为 $t$ 的函数，让微分方程 $f(\mathbf{z}(t),t)$ 随 $t$ 变化。这种参数化是一种超网络（hypernetwork）（Ha et al. 2016）。我们还为每个隐藏单元引入了门控机制，$\frac{d\mathbf{z}}{dt}=\sum_{n}\sigma_{n}(t)f_{n}(\mathbf{z})$，其中 $\sigma_{n}(t)\in(0,1)$ 是一个神经网络，学习何时应当施加动力学 $f_{n}(\mathbf{z})$。我们把这些模型称为连续归一化流（continuous normalizing flows，CNF）。

4.1 Experiments with Continuous Normalizing Flows

4.1 连续归一化流实验

(a) Target　(b) NF: K=2, K=8, K=32　(c) CNF: M=2, M=8, M=32　(d) Loss vs. K/M (rows 1, 2, 3)

Figure 4: Comparison of normalizing flows versus continuous normalizing flows. The model capacity of normalizing flows is determined by their depth (K), while continuous normalizing flows can also increase capacity by increasing width (M), making them easier to train.

（a）目标　（b）NF：K=2、K=8、K=32　（c）CNF：M=2、M=8、M=32　（d）损失与 K/M 的关系（第 1、2、3 行）

图 4：归一化流与连续归一化流的比较。归一化流的模型容量由深度（K）决定，而连续归一化流还可以通过增加宽度（M）来提升容量，因而更容易训练。

We first compare continuous and discrete planar flows at learning to sample from a known distribution. We show that a planar CNF with $M$ hidden units can be at least as expressive as a planar NF with $K=M$ layers, and sometimes much more expressive.

我们先比较连续和离散的平面流学习从已知分布中采样的能力。我们表明，有 $M$ 个隐藏单元的平面 CNF，表达能力至少与有 $K=M$ 层的平面 NF 相当，有时还强得多。

Density matching

We configure the CNF as described above, and train for 10,000 iterations using Adam (Kingma and Ba 2014). In contrast, the NF is trained for 500,000 iterations using RMSprop (Hinton et al. 2012), as suggested by Rezende and Mohamed 2015. For this task, we minimize $\mathrm{KL}\left(q(\mathbf{x})\middle\|p(\mathbf{x})\right)$ as the loss function where $q$ is the flow model and the target density $p(\cdot)$ can be evaluated. Figure 4 shows that CNF generally achieves lower loss.

密度匹配

我们按上文所述配置 CNF，用 Adam（Kingma and Ba 2014）训练 10,000 次迭代。相比之下，按照 Rezende and Mohamed 2015 的建议，NF 用 RMSprop（Hinton et al. 2012）训练 500,000 次迭代。这个任务的损失函数是最小化 $\mathrm{KL}\left(q(\mathbf{x})\middle\|p(\mathbf{x})\right)$，其中 $q$ 是流模型，目标密度 $p(\cdot)$ 可以直接计算。图 4 显示，CNF 通常能达到更低的损失。

Maximum Likelihood Training

(a) Two Circles　(b) Two Moons (each panel: Density and Samples at 5%, 20%, 40%, 60%, 80%, 100%; NF; Target)

Figure 6: Visualizing the transformation from noise to data. Continuous-time normalizing flows are reversible, so we can train on a density estimation task and still be able to sample from the learned density efficiently.

最大似然训练

（a）双圆　（b）双月（每组子图：5%、20%、40%、60%、80%、100% 处的密度与样本；NF；目标）

图 6：从噪声到数据的变换过程可视化。连续时间归一化流是可逆的，因此我们可以在密度估计任务上训练，同时仍能高效地从学到的密度中采样。

A useful property of continuous-time normalizing flows is that we can compute the reverse transformation for about the same cost as the forward pass, which cannot be said for normalizing flows. This lets us train the flow on a density estimation task by performing maximum likelihood estimation, which maximizes $\mathbb{E}_{p(\mathbf{x})}[\log q(\mathbf{x})]$ where $q(\cdot)$ is computed using the appropriate change of variables theorem, then afterwards reverse the CNF to generate random samples from $q(\mathbf{x})$.

连续时间归一化流有一个有用的性质：计算逆变换的开销与前向传播大致相同，而归一化流做不到这一点。这让我们可以在密度估计任务上用最大似然估计训练流，即最大化 $\mathbb{E}_{p(\mathbf{x})}[\log q(\mathbf{x})]$，其中 $q(\cdot)$ 用相应的变量替换定理计算；训练完成后，再把 CNF 反过来，从 $q(\mathbf{x})$ 中生成随机样本。

For this task, we use 64 hidden units for CNF, and 64 stacked one-hidden-unit layers for NF. Figure 6 shows the learned dynamics. Instead of showing the initial Gaussian distribution, we display the transformed distribution after a small amount of time which shows the locations of the initial planar flows. Interestingly, to fit the Two Circles distribution, the CNF rotates the planar flows so that the particles can be evenly spread into circles. While the CNF transformations are smooth and interpretable, we find that NF transformations are very unintuitive and this model has difficulty fitting the two moons dataset in Figure 6(b).

这个任务中，CNF 使用 64 个隐藏单元，NF 使用 64 个堆叠的单隐藏单元层。图 6 展示了学到的动力学。我们没有画出初始的高斯分布，而是展示经过一小段时间后的变换分布，以显示初始平面流的位置。有意思的是，为了拟合双圆分布，CNF 旋转了平面流，使粒子能均匀地分散到圆上。CNF 的变换平滑且可解释，而我们发现 NF 的变换很不直观，在图 6(b) 中也难以拟合双月数据集。

5 A generative latent function time-series model

5 一种生成式隐函数时间序列模型

Applying neural networks to irregularly-sampled data such as medical records, network traffic, or neural spiking data is difficult. Typically, observations are put into bins of fixed duration, and the latent dynamics are discretized in the same way. This leads to difficulties with missing data and ill-defined latent variables. Missing data can be addressed using generative time-series models (Álvarez and Lawrence 2011; Futoma et al. 2017; Mei and Eisner 2017; Soleimani et al. 2017a) or data imputation  (Che et al. 2018). Another approach concatenates time-stamp information to the input of an RNN (Choi et al. 2016; Lipton et al. 2016; Du et al. 2016; Li 2017).

把神经网络应用到病历、网络流量或神经元放电数据这类不规则采样的数据上很困难。通常的做法是把观测放进固定时长的区间，并以同样的方式把潜在动力学离散化。这会带来数据缺失和隐变量定义不清的问题。数据缺失可以用生成式时间序列模型（Álvarez and Lawrence 2011; Futoma et al. 2017; Mei and Eisner 2017; Soleimani et al. 2017a）或数据插补（Che et al. 2018）来处理。另一种方法是把时间戳信息拼接到 RNN 的输入中（Choi et al. 2016; Lipton et al. 2016; Du et al. 2016; Li 2017）。

We present a continuous-time, generative approach to modeling time series. Our model represents each time series by a latent trajectory. Each trajectory is determined from a local initial state, $\mathbf{z}_{t_{0}}$, and a global set of latent dynamics shared across all time series. Given observation times $t_{0},t_{1},\dots,t_{N}$ and an initial state $\mathbf{z}_{t_{0}}$, an ODE solver produces $\mathbf{z}_{t_{1}},\dots,\mathbf{z}_{t_{N}}$, which describe the latent state at each observation.We define this generative model formally through a sampling procedure:

我们提出一种连续时间的生成式时间序列建模方法。模型用一条隐轨迹（latent trajectory）表示每个时间序列。每条轨迹由一个局部的初始状态 $\mathbf{z}_{t_{0}}$，以及一组所有时间序列共享的全局潜在动力学决定。给定观测时刻 $t_{0},t_{1},\dots,t_{N}$ 和初始状态 $\mathbf{z}_{t_{0}}$，ODE 求解器给出描述每个观测时刻隐状态的 $\mathbf{z}_{t_{1}},\dots,\mathbf{z}_{t_{N}}$。我们通过如下采样过程正式定义这个生成模型：

$$\begin{aligned}\mathbf{z}_{t_{0}} &\sim p(\mathbf{z}_{t_{0}}) \\ \mathbf{z}_{t_{1}},\mathbf{z}_{t_{2}},\dots,\mathbf{z}_{t_{N}} &=\textnormal{ODESolve}(\mathbf{z}_{t_{0}},f,\theta_{f},t_{0},\dots,t_{N}) \\ \textnormal{each}\quad\mathbf{x}_{t_{i}} &\sim p(\mathbf{x}|\mathbf{z}_{t_{i}},\theta_{\mathbf{x}})\end{aligned} \tag{11}$$

Figure 7: Computation graph of the latent ODE model.

图 7：隐 ODE 模型的计算图。

Function $f$ is a time-invariant function that takes the value $\mathbf{z}$ at the current time step and outputs the gradient: ${{\partial\mathbf{z}(t)}}/{{\partial t}}=f(\mathbf{z}(t),\theta_{f})$. We parametrize this function using a neural net. Because $f$ is time-invariant, given any latent state $\mathbf{z}(t)$, the entire latent trajectory is uniquely defined. Extrapolating this latent trajectory lets us make predictions arbitrarily far forwards or backwards in time.

函数 $f$ 是一个不随时间变化的函数，它接收当前时间步的值 $\mathbf{z}$，输出梯度：$\partial\mathbf{z}(t)/\partial t=f(\mathbf{z}(t),\theta_{f})$。我们用神经网络参数化这个函数。由于 $f$ 不随时间变化，给定任一隐状态 $\mathbf{z}(t)$，整条隐轨迹就唯一确定了。外推这条隐轨迹，就可以对任意远的未来或过去做预测。

Training and Prediction

We can train this latent-variable model as a variational autoencoder (Kingma and Welling 2014; Rezende et al. 2014), with sequence-valued observations. Our recognition net is an RNN, which consumes the data sequentially backwards in time, and outputs $q_{\phi}(\mathbf{z}_{0}|\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{N})$. A detailed algorithm can be found in Appendix E. Using ODEs as a generative model allows us to make predictions for arbitrary time points $t_{1}$…$t_{M}$ on a continuous timeline.

训练与预测

我们可以把这个隐变量模型当作观测为序列的变分自编码器（Kingma and Welling 2014; Rezende et al. 2014）来训练。我们的识别网络是一个 RNN，它沿时间反向依次读入数据，输出 $q_{\phi}(\mathbf{z}_{0}|\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{N})$。详细算法见附录 E。用 ODE 作为生成模型，让我们能在连续时间轴上对任意时刻 $t_{1}$…$t_{M}$ 做预测。

Figure 8: Fitting a latent ODE dynamics model with a Poisson process likelihood. Dots show event times. The line is the learned intensity $\lambda(t)$ of the Poisson process.

图 8：用泊松过程似然拟合隐 ODE 动力学模型。圆点表示事件发生的时刻，曲线是学到的泊松过程强度 $\lambda(t)$。

Poisson Process likelihoods

The fact that an observation occurred often tells us something about the latent state. For example, a patient may be more likely to take a medical test if they are sick. The rate of events can be parameterized by a function of the latent state: $p(\textnormal{event at time } t\,|\,\mathbf{z}(t))=\lambda(\mathbf{z}(t))$. Given this rate function, the likelihood of a set of independent observation times in the interval $[t_{\textnormal{start}},t_{\textnormal{end}}]$ is given by an inhomogeneous Poisson process (Palm 1943):

泊松过程似然

“发生了一次观测”这件事本身，往往就透露了隐状态的某些信息。例如，病人生病时更可能去做医学检查。事件发生率可以参数化为隐状态的函数：$p(\textnormal{event at time } t\,|\,\mathbf{z}(t))=\lambda(\mathbf{z}(t))$。给定这个率函数，区间 $[t_{\textnormal{start}},t_{\textnormal{end}}]$ 内一组独立观测时刻的似然由非齐次泊松过程（Palm 1943）给出：

$$\log p(t_{1}\dots t_{N}|\,t_{\textnormal{start}},t_{\textnormal{end}})=\sum_{i=1}^{N}\log\lambda(\mathbf{z}(t_{i}))-\int_{t_{\textnormal{start}}}^{t_{\textnormal{end}}}\!\lambda(\mathbf{z}(t))dt$$

(a) Recurrent Neural Network　(b) Latent Neural Ordinary Differential Equation　(c) Latent Trajectories

Figure 9: (9(a)): Reconstruction and extrapolation of spirals with irregular time points by a recurrent neural network. (9(b)): Reconstructions and extrapolations by a latent neural ODE. Blue curve shows model prediction. Red shows extrapolation. (9(c)) A projection of inferred 4-dimensional latent ODE trajectories onto their first two dimensions. Color indicates the direction of the corresponding trajectory.

（a）循环神经网络　（b）隐神经常微分方程　（c）隐轨迹

图 9：（9(a)）循环神经网络对不规则时间点上螺旋线的重建与外推。（9(b)）隐神经 ODE 的重建与外推。蓝色曲线为模型预测，红色为外推部分。（9(c)）推断出的 4 维隐 ODE 轨迹在前两个维度上的投影。颜色表示相应轨迹的方向。

We can parameterize $\lambda(\cdot)$ using another neural network. Conveniently, we can evaluate both the latent trajectory and the Poisson process likelihood together in a single call to an ODE solver. Figure 8 shows the event rate learned by such a model on a toy dataset.

我们可以用另一个神经网络来参数化 $\lambda(\cdot)$。方便的是，隐轨迹和泊松过程似然可以在一次 ODE 求解器调用中一起计算。图 8 展示了这样一个模型在一个玩具数据集上学到的事件率。

A Poisson process likelihood on observation times can be combined with a data likelihood to jointly model all observations and the times at which they were made.

观测时刻上的泊松过程似然可以与数据似然结合，对所有观测及其发生时刻联合建模。

5.1 Time-series Latent ODE Experiments

5.1 时间序列隐 ODE 实验

We investigate the ability of the latent ODE model to fit and extrapolate time series. The recognition network is an RNN with 25 hidden units. We use a 4-dimensional latent space. We parameterize the dynamics function $f$ with a one-hidden-layer network with 20 hidden units. The decoder computing $p(\mathbf{x}_{t_{i}}|\mathbf{z}_{t_{i}})$ is another neural network with one hidden layer with 20 hidden units. Our baseline was a recurrent neural net with 25 hidden units trained to minimize negative Gaussian log-likelihood. We trained a second version of this RNN whose inputs were concatenated with the time difference to the next observation to aid RNN with irregular observations.

我们考察隐 ODE 模型拟合和外推时间序列的能力。识别网络是一个有 25 个隐藏单元的 RNN。我们使用 4 维隐空间。动力学函数 $f$ 用一个有 20 个隐藏单元的单隐藏层网络参数化。计算 $p(\mathbf{x}_{t_{i}}|\mathbf{z}_{t_{i}})$ 的解码器是另一个有 20 个隐藏单元的单隐藏层神经网络。基线是一个有 25 个隐藏单元的循环神经网络，训练目标是最小化高斯负对数似然。我们还训练了这个 RNN 的第二个版本，在其输入中拼接到下一次观测的时间差，以帮助 RNN 处理不规则观测。

Bi-directional spiral dataset

We generated a dataset of 1000 2-dimensional spirals, each starting at a different point, sampled at 100 equally-spaced timesteps. The dataset contains two types of spirals: half are clockwise while the other half counter-clockwise. To make the task more realistic, we add gaussian noise to the observations.

双向螺旋数据集

我们生成了一个包含 1000 条二维螺旋线的数据集，每条起点不同，在 100 个等间隔时间步上采样。数据集包含两类螺旋：一半顺时针，另一半逆时针。为了让任务更贴近现实，我们在观测中加入了高斯噪声。

Time series with irregular time points

To generate irregular timestamps, we randomly sample points from each trajectory without replacement ($n=\{30,50,100\}$). We report predictive root-mean-squared error (RMSE) on 100 time points extending beyond those that were used for training. Table 2 shows that the latent ODE has substantially lower predictive RMSE.

不规则时间点的时间序列

为了生成不规则的时间戳，我们从每条轨迹中无放回地随机采样若干点（$n=\{30,50,100\}$）。我们报告在训练所用时间点之后延伸的 100 个时间点上的预测均方根误差（RMSE）。表 2 显示，隐 ODE 的预测 RMSE 明显更低。

# Observations	30/100	50/100	100/100
RNN	0.3937	0.3202	0.1813
Latent ODE	0.1642	0.1502	0.1346

Table 2: Predictive RMSE on test set

观测数	30/100	50/100	100/100
RNN	0.3937	0.3202	0.1813
隐 ODE	0.1642	0.1502	0.1346

表 2：测试集上的预测 RMSE

Figure 9 shows examples of spiral reconstructions with 30 sub-sampled points. Reconstructions from the latent ODE were obtained by sampling from the posterior over latent trajectories and decoding it to data-space. Examples with varying number of time points are shown in Appendix F. We observed that reconstructions and extrapolations are consistent with the ground truth regardless of number of observed points and despite the noise.

图 9 给出了用 30 个子采样点重建螺旋线的示例。隐 ODE 的重建结果是从隐轨迹的后验中采样，再解码到数据空间得到的。观测点数不同的示例见附录 F。我们观察到，无论观测点数多少、是否有噪声，重建和外推结果都与真实值一致。

Latent space interpolation

Figure 9(c) shows latent trajectories projected onto the first two dimensions of the latent space. The trajectories form two separate clusters of trajectories, one decoding to clockwise spirals, the other to counter-clockwise. Figure 10 shows that the latent trajectories change smoothly as a function of the initial point $\mathbf{z}(t_{0})$, switching from a clockwise to a counter-clockwise spiral.

隐空间插值

图 9(c) 展示了投影到隐空间前两个维度上的隐轨迹。这些轨迹形成两个分开的簇，一簇解码为顺时针螺旋，另一簇解码为逆时针螺旋。图 10 显示，隐轨迹随初始点 $\mathbf{z}(t_{0})$ 平滑变化，从顺时针螺旋切换为逆时针螺旋。

Figure 10: Data-space trajectories decoded from varying one dimension of $\mathbf{z}_{t_{0}}$. Color indicates progression through time, starting at purple and ending at red. Note that the trajectories on the left are counter-clockwise, while the trajectories on the right are clockwise.

图 10：改变 $\mathbf{z}_{t_{0}}$ 的一个维度后解码得到的数据空间轨迹。颜色表示时间进程，从紫色开始，到红色结束。注意，左边的轨迹是逆时针的，右边的轨迹是顺时针的。

6 Scope and Limitations

6 适用范围与局限

Minibatching

The use of mini-batches is less straightforward than for standard neural networks. One can still batch together evaluations through the ODE solver by concatenating the states of each batch element together, creating a combined ODE with dimension $D\times K$. In some cases, controlling error on all batch elements together might require evaluating the combined system $K$ times more often than if each system was solved individually. However, in practice the number of evaluations did not increase substantially when using minibatches.

小批量

使用小批量不像在标准神经网络中那样直接。仍然可以把每个批次元素的状态拼接起来，构成一个维度为 $D\times K$ 的联合 ODE，从而通过 ODE 求解器批量求值。在某些情况下，同时控制所有批次元素的误差，可能需要对联合系统求值的次数比单独求解每个系统多 $K$ 倍。不过实践中，使用小批量时求值次数并没有明显增加。

Uniqueness

When do continuous dynamics have a unique solution? Picard’s existence theorem (Coddington and Levinson 1955) states that the solution to an initial value problem exists and is unique if the differential equation is uniformly Lipschitz continuous in $\mathbf{z}$ and continuous in $t$. This theorem holds for our model if the neural network has finite weights and uses Lipshitz nonlinearities, such as tanh or relu.

唯一性

连续动力学什么时候有唯一解？皮卡存在性定理（Picard’s existence theorem）（Coddington and Levinson 1955）指出，如果微分方程关于 $\mathbf{z}$ 一致利普希茨连续、关于 $t$ 连续，那么初值问题的解存在且唯一。只要神经网络的权重有限，并使用 tanh 或 relu 这类利普希茨非线性，这个定理对我们的模型就成立。

Setting tolerances

Our framework allows the user to trade off speed for precision, but requires the user to choose an error tolerance on both the forward and reverse passes during training. For sequence modeling, the default value of 1.5e-8 was used. In the classification and density estimation experiments, we were able to reduce the tolerance to 1e-3 and 1e-5, respectively, without degrading performance.

设置容差

我们的框架允许用户在速度和精度之间权衡，但要求用户在训练时为前向和反向传播都选定误差容差。序列建模使用默认值 1.5e-8。在分类和密度估计实验中，我们分别把容差降到 1e-3 和 1e-5，性能没有下降。

Reconstructing forward trajectories

Reconstructing the state trajectory by running the dynamics backwards can introduce extra numerical error if the reconstructed trajectory diverges from the original. This problem can be addressed by checkpointing: storing intermediate values of $\mathbf{z}$ on the forward pass, and reconstructing the exact forward trajectory by re-integrating from those points. We did not find this to be a practical problem, and we informally checked that reversing many layers of continuous normalizing flows with default tolerances recovered the initial states.

重建前向轨迹

通过反向运行动力学来重建状态轨迹，如果重建的轨迹偏离原轨迹，就可能引入额外的数值误差。这个问题可以用检查点（checkpointing）解决：在前向传播时保存 $\mathbf{z}$ 的中间值，再从这些点重新积分，精确重建前向轨迹。我们没有发现这在实践中构成问题；我们还非正式地检查过，用默认容差把多层连续归一化流反过来运行，能够恢复初始状态。

7 Related Work

7 相关工作

The use of the adjoint method for training continuous-time neural networks was previously proposed (LeCun et al. 1988; Pearlmutter 1995), though was not demonstrated practically. The interpretation of residual networks He et al. 2016a as approximate ODE solvers spurred research into exploiting reversibility and approximate computation in ResNets (Chang et al. 2017; Lu et al. 2017). We demonstrate these same properties in more generality by directly using an ODE solver.

用伴随方法训练连续时间神经网络此前已有人提出（LeCun et al. 1988; Pearlmutter 1995），但没有在实践中得到验证。把残差网络（He et al. 2016a）解释为近似 ODE 求解器，推动了在 ResNet 中利用可逆性和近似计算的研究（Chang et al. 2017; Lu et al. 2017）。我们直接使用 ODE 求解器，在更一般的意义上展示了同样的性质。

Adaptive computation

One can adapt computation time by training secondary neural networks to choose the number of evaluations of recurrent or residual networks (Graves 2016; Jernite et al. 2016; Figurnov et al. 2017; Chang et al. 2018). However, this introduces overhead both at training and test time, and extra parameters that need to be fit. In contrast, ODE solvers offer well-studied, computationally cheap, and generalizable rules for adapting the amount of computation.

自适应计算

可以训练辅助神经网络来选择循环网络或残差网络的求值次数，从而自适应地调整计算时间（Graves 2016; Jernite et al. 2016; Figurnov et al. 2017; Chang et al. 2018）。但这在训练和测试时都会引入额外开销，还需要拟合额外的参数。相比之下，ODE 求解器提供了经过充分研究、计算成本低、可推广的规则来调整计算量。

Constant memory backprop through reversibility

Recent work developed reversible versions of residual networks (Gomez et al. 2017; Haber and Ruthotto 2017; Chang et al. 2017), which gives the same constant memory advantage as our approach. However, these methods require restricted architectures, which partition the hidden units. Our approach does not have these restrictions.

借助可逆性实现恒定内存的反向传播

近来的工作开发了残差网络的可逆版本（Gomez et al. 2017; Haber and Ruthotto 2017; Chang et al. 2017），它们与我们的方法一样具有恒定内存的优势。然而，这些方法需要受限的架构，要对隐藏单元做划分。我们的方法没有这些限制。

Learning differential equations

Much recent work has proposed learning differential equations from data. One can train feed-forward or recurrent neural networks to approximate a differential equation (Raissi and Karniadakis 2018; Raissi et al. 2018a; Long et al. 2017), with applications such as fluid simulation (Wiewel et al. 2018). There is also significant work on connecting Gaussian Processes (GPs) and ODE solvers (Schober et al. 2014). GPs have been adapted to fit differential equations (Raissi et al. 2018b) and can naturally model continuous-time effects and interventions (Soleimani et al. 2017b; Schulam and Saria 2017). Ryder et al. 2018 use stochastic variational inference to recover the solution of a given stochastic differential equation.

学习微分方程

近来很多工作提出从数据中学习微分方程。可以训练前馈或循环神经网络来逼近微分方程（Raissi and Karniadakis 2018; Raissi et al. 2018a; Long et al. 2017），应用包括流体模拟（Wiewel et al. 2018）。把高斯过程（GP）与 ODE 求解器联系起来的工作也很多（Schober et al. 2014）。GP 已被改造用于拟合微分方程（Raissi et al. 2018b），并能自然地建模连续时间效应和干预（Soleimani et al. 2017b; Schulam and Saria 2017）。Ryder et al. 2018 用随机变分推断恢复给定随机微分方程的解。

Differentiating through ODE solvers

The dolfin library (Farrell et al. 2013) implements adjoint computation for general ODE and PDE solutions, but only by backpropagating through the individual operations of the forward solver. The Stan library (Carpenter et al. 2015) implements gradient estimation through ODE solutions using forward sensitivity analysis. However, forward sensitivity analysis is quadratic-time in the number of variables, whereas the adjoint sensitivity analysis is linear (Carpenter et al. 2015; Zhang and Sandu 2014). Melicher et al. 2017 used the adjoint method to train bespoke latent dynamic models.

对 ODE 求解器求导

dolfin 库（Farrell et al. 2013）为一般的 ODE 和 PDE 解实现了伴随计算，但只能通过对前向求解器的各个运算做反向传播来实现。Stan 库（Carpenter et al. 2015）用前向灵敏度分析实现了穿过 ODE 解的梯度估计。然而，前向灵敏度分析的时间与变量数呈平方关系，而伴随灵敏度分析是线性的（Carpenter et al. 2015; Zhang and Sandu 2014）。Melicher et al. 2017 用伴随方法训练专门定制的潜在动力学模型。

In contrast, by providing a generic vector-Jacobian product, we allow an ODE solver to be trained end-to-end with any other differentiable model components. While use of vector-Jacobian products for solving the adjoint method has been explored in optimal control (Andersson 2013; Andersson et al. In Press, 2018), we highlight the potential of a general integration of black-box ODE solvers into automatic differentiation (Baydin et al. 2018) for deep learning and generative modeling.

相比之下，我们通过提供通用的向量-雅可比积，让 ODE 求解器可以与任何其他可微模型组件一起端到端训练。用向量-雅可比积求解伴随方法在最优控制中已有探索（Andersson 2013; Andersson et al. In Press, 2018），我们则强调，把黑箱 ODE 求解器全面纳入自动微分（Baydin et al. 2018），对深度学习和生成建模有很大潜力。

8 Conclusion

8 结论

We investigated the use of black-box ODE solvers as a model component, developing new models for time-series modeling, supervised learning, and density estimation. These models are evaluated adaptively, and allow explicit control of the tradeoff between computation speed and accuracy. Finally, we derived an instantaneous version of the change of variables formula, and developed continuous-time normalizing flows, which can scale to large layer sizes.

我们研究了把黑箱 ODE 求解器用作模型组件，为时间序列建模、监督学习和密度估计开发了新模型。这些模型以自适应方式求值，并允许显式控制计算速度与精度之间的权衡。最后，我们推导了变量替换公式的瞬时版本，并开发了可以扩展到大层宽的连续时间归一化流。

9 Acknowledgements

9 致谢

We thank Wenyi Wang and Geoff Roeder for help with proofs, and Daniel Duckworth, Ethan Fetaya, Hossein Soleimani, Eldad Haber, Ken Caluwaerts, Daniel Flam-Shepherd, and Harry Braviner for feedback. We thank Chris Rackauckas, Dougal Maclaurin, and Matthew James Johnson for helpful discussions. We also thank Yuval Frommer for pointing out an unsupported claim about parameter efficiency.

我们感谢 Wenyi Wang 和 Geoff Roeder 在证明上的帮助，感谢 Daniel Duckworth、Ethan Fetaya、Hossein Soleimani、Eldad Haber、Ken Caluwaerts、Daniel Flam-Shepherd 和 Harry Braviner 的反馈。感谢 Chris Rackauckas、Dougal Maclaurin 和 Matthew James Johnson 的有益讨论。我们还感谢 Yuval Frommer 指出了一处关于参数效率的缺乏依据的说法。

References

参考文献

Álvarez and Lawrence (2011) Mauricio A Álvarez and Neil D Lawrence. Computationally efficient convolved multiple output Gaussian processes. Journal of Machine Learning Research, 12(May):1459–1500, 2011.

Amos and Kolter (2017) Brandon Amos and J Zico Kolter. OptNet: Differentiable optimization as a layer in neural networks. In International Conference on Machine Learning, pages 136–145, 2017.

Andersson (2013) Joel Andersson. A general-purpose software framework for dynamic optimization. PhD thesis, 2013.

Andersson et al. (In Press, 2018) Joel A E Andersson, Joris Gillis, Greg Horn, James B Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, In Press, 2018.

Baydin et al. (2018) Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: a survey. Journal of machine learning research, 18(153):1–153, 2018.

Berg et al. (2018) Rianne van den Berg, Leonard Hasenclever, Jakub M Tomczak, and Max Welling. Sylvester normalizing flows for variational inference. arXiv preprint arXiv:1803.05649, 2018.

Carpenter et al. (2015) Bob Carpenter, Matthew D Hoffman, Marcus Brubaker, Daniel Lee, Peter Li, and Michael Betancourt. The Stan math library: Reverse-mode automatic differentiation in c++. arXiv preprint arXiv:1509.07164, 2015.

Chang et al. (2017) Bo Chang, Lili Meng, Eldad Haber, Lars Ruthotto, David Begert, and Elliot Holtham. Reversible architectures for arbitrarily deep residual neural networks. arXiv preprint arXiv:1709.03698, 2017.

Chang et al. (2018) Bo Chang, Lili Meng, Eldad Haber, Frederick Tung, and David Begert. Multi-level residual networks from dynamical systems view. In International Conference on Learning Representations, 2018. URL https://openreview.net/forum?id=SyJS-OgR-.

Che et al. (2018) Zhengping Che, Sanjay Purushotham, Kyunghyun Cho, David Sontag, and Yan Liu. Recurrent neural networks for multivariate time series with missing values. Scientific Reports, 8(1):6085, 2018. URL https://doi.org/10.1038/s41598-018-24271-9.

Choi et al. (2016) Edward Choi, Mohammad Taha Bahadori, Andy Schuetz, Walter F. Stewart, and Jimeng Sun. Doctor AI: Predicting clinical events via recurrent neural networks. In Proceedings of the 1st Machine Learning for Healthcare Conference, volume 56 of Proceedings of Machine Learning Research, pages 301–318. PMLR, 18–19 Aug 2016. URL http://proceedings.mlr.press/v56/Choi16.html.

Coddington and Levinson (1955) Earl A Coddington and Norman Levinson. Theory of ordinary differential equations. Tata McGraw-Hill Education, 1955.

Dinh et al. (2014) Laurent Dinh, David Krueger, and Yoshua Bengio. NICE: Non-linear independent components estimation. arXiv preprint arXiv:1410.8516, 2014.

Du et al. (2016) Nan Du, Hanjun Dai, Rakshit Trivedi, Utkarsh Upadhyay, Manuel Gomez-Rodriguez, and Le Song. Recurrent marked temporal point processes: Embedding event history to vector. In International Conference on Knowledge Discovery and Data Mining, pages 1555–1564. ACM, 2016.

Farrell et al. (2013) Patrick Farrell, David Ham, Simon Funke, and Marie Rognes. Automated derivation of the adjoint of high-level transient finite element programs. SIAM Journal on Scientific Computing, 2013.

Figurnov et al. (2017) Michael Figurnov, Maxwell D Collins, Yukun Zhu, Li Zhang, Jonathan Huang, Dmitry Vetrov, and Ruslan Salakhutdinov. Spatially adaptive computation time for residual networks. arXiv preprint, 2017.

Futoma et al. (2017) J. Futoma, S. Hariharan, and K. Heller. Learning to Detect Sepsis with a Multitask Gaussian Process RNN Classifier. ArXiv e-prints, 2017.

Gomez et al. (2017) Aidan N Gomez, Mengye Ren, Raquel Urtasun, and Roger B Grosse. The reversible residual network: Backpropagation without storing activations. In Advances in Neural Information Processing Systems, pages 2211–2221, 2017.

Graves (2016) Alex Graves. Adaptive computation time for recurrent neural networks. arXiv preprint arXiv:1603.08983, 2016.

Ha et al. (2016) David Ha, Andrew Dai, and Quoc V Le. Hypernetworks. arXiv preprint arXiv:1609.09106, 2016.

Haber and Ruthotto (2017) Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse Problems, 34(1):014004, 2017.

Hairer et al. (1987) E. Hairer, S.P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I – Nonstiff Problems. Springer, 1987.

He et al. (2016a) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016a.

He et al. (2016b) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Identity mappings in deep residual networks. In European conference on computer vision, pages 630–645. Springer, 2016b.

Hinton et al. (2012) Geoffrey Hinton, Nitish Srivastava, and Kevin Swersky. Neural networks for machine learning lecture 6a overview of mini-batch gradient descent, 2012.

Jernite et al. (2016) Yacine Jernite, Edouard Grave, Armand Joulin, and Tomas Mikolov. Variable computation in recurrent neural networks. arXiv preprint arXiv:1611.06188, 2016.

Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.

Kingma and Welling (2014) Diederik P. Kingma and Max Welling. Auto-encoding variational Bayes. International Conference on Learning Representations, 2014.

Kingma et al. (2016) Diederik P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling. Improved variational inference with inverse autoregressive flow. In Advances in Neural Information Processing Systems, pages 4743–4751, 2016.

Kutta (1901) W. Kutta. Beitrag zur näherungsweisen Integration totaler Differentialgleichungen. Zeitschrift für Mathematik und Physik, 46:435–453, 1901.

LeCun et al. (1988) Yann LeCun, D Touresky, G Hinton, and T Sejnowski. A theoretical framework for back-propagation. In Proceedings of the 1988 connectionist models summer school, volume 1, pages 21–28. CMU, Pittsburgh, Pa: Morgan Kaufmann, 1988.

LeCun et al. (1998) Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278–2324, 1998.

Li (2017) Yang Li. Time-dependent representation for neural event sequence prediction. arXiv preprint arXiv:1708.00065, 2017.

Lipton et al. (2016) Zachary C Lipton, David Kale, and Randall Wetzel. Directly modeling missing data in sequences with RNNs: Improved classification of clinical time series. In Proceedings of the 1st Machine Learning for Healthcare Conference, volume 56 of Proceedings of Machine Learning Research, pages 253–270. PMLR, 18–19 Aug 2016. URL http://proceedings.mlr.press/v56/Lipton16.html.

Long et al. (2017) Z. Long, Y. Lu, X. Ma, and B. Dong. PDE-Net: Learning PDEs from Data. ArXiv e-prints, 2017.

Lu et al. (2017) Yiping Lu, Aoxiao Zhong, Quanzheng Li, and Bin Dong. Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations. arXiv preprint arXiv:1710.10121, 2017.

Maclaurin et al. (2015) Dougal Maclaurin, David Duvenaud, and Ryan P Adams. Autograd: Reverse-mode differentiation of native Python. In ICML workshop on Automatic Machine Learning, 2015.

Mei and Eisner (2017) Hongyuan Mei and Jason M Eisner. The neural Hawkes process: A neurally self-modulating multivariate point process. In Advances in Neural Information Processing Systems, pages 6757–6767, 2017.

Melicher et al. (2017) Valdemar Melicher, Tom Haber, and Wim Vanroose. Fast derivatives of likelihood functionals for ODE based models using adjoint-state method. Computational Statistics, 32(4):1621–1643, 2017.

Palm (1943) Conny Palm. Intensitätsschwankungen im fernsprechverker. Ericsson Technics, 1943.

Paszke et al. (2017) Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. 2017.

Pearlmutter (1995) Barak A Pearlmutter. Gradient calculations for dynamic recurrent neural networks: A survey. IEEE Transactions on Neural networks, 6(5):1212–1228, 1995.

Pontryagin et al. (1962) Lev Semenovich Pontryagin, EF Mishchenko, VG Boltyanskii, and RV Gamkrelidze. The mathematical theory of optimal processes. 1962.

Raissi and Karniadakis (2018) M. Raissi and G. E. Karniadakis. Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics, pages 125–141, 2018.

Raissi et al. (2018a) Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Multistep neural networks for data-driven discovery of nonlinear dynamical systems. arXiv preprint arXiv:1801.01236, 2018a.

Raissi et al. (2018b) Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Numerical Gaussian processes for time-dependent and nonlinear partial differential equations. SIAM Journal on Scientific Computing, 40(1):A172–A198, 2018b.

Rezende et al. (2014) Danilo J Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic backpropagation and approximate inference in deep generative models. In Proceedings of the 31st International Conference on Machine Learning, pages 1278–1286, 2014.

Rezende and Mohamed (2015) Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows. arXiv preprint arXiv:1505.05770, 2015.

Runge (1895) C. Runge. Über die numerische Auflösung von Differentialgleichungen. Mathematische Annalen, 46:167–178, 1895.

Ruthotto and Haber (2018) Lars Ruthotto and Eldad Haber. Deep neural networks motivated by partial differential equations. arXiv preprint arXiv:1804.04272, 2018.

Ryder et al. (2018) T. Ryder, A. Golightly, A. S. McGough, and D. Prangle. Black-box Variational Inference for Stochastic Differential Equations. ArXiv e-prints, 2018.

Schober et al. (2014) Michael Schober, David Duvenaud, and Philipp Hennig. Probabilistic ODE solvers with Runge-Kutta means. In Advances in Neural Information Processing Systems 25, 2014.

Schulam and Saria (2017) Peter Schulam and Suchi Saria. What-if reasoning with counterfactual Gaussian processes. arXiv preprint arXiv:1703.10651, 2017.

Soleimani et al. (2017a) Hossein Soleimani, James Hensman, and Suchi Saria. Scalable joint models for reliable uncertainty-aware event prediction. IEEE transactions on pattern analysis and machine intelligence, 2017a.

Soleimani et al. (2017b) Hossein Soleimani, Adarsh Subbaswamy, and Suchi Saria. Treatment-response models for counterfactual reasoning with continuous-time, continuous-valued interventions. arXiv preprint arXiv:1704.02038, 2017b.

Stam (1999) Jos Stam. Stable fluids. In Proceedings of the 26th annual conference on Computer graphics and interactive techniques, pages 121–128. ACM Press/Addison-Wesley Publishing Co., 1999.

Stapor et al. (2018) Paul Stapor, Fabian Froehlich, and Jan Hasenauer. Optimization and uncertainty analysis of ODE models using second order adjoint sensitivity analysis. bioRxiv, page 272005, 2018.

Tomczak and Welling (2016) Jakub M Tomczak and Max Welling. Improving variational auto-encoders using Householder flow. arXiv preprint arXiv:1611.09630, 2016.

Wiewel et al. (2018) Steffen Wiewel, Moritz Becher, and Nils Thuerey. Latent-space physics: Towards learning the temporal evolution of fluid flow. arXiv preprint arXiv:1802.10123, 2018.

Zhang and Sandu (2014) Hong Zhang and Adrian Sandu. Fatode: a library for forward, adjoint, and tangent linear integration of ODEs. SIAM Journal on Scientific Computing, 36(5):C504–C523, 2014.

Appendix A Proof of the Instantaneous Change of Variables Theorem

附录 A 瞬时变量替换定理的证明

Theorem (Instantaneous Change of Variables).

Let $\mathbf{z}(t)$ be a finite continuous random variable with probability $p(\mathbf{z}(t))$ dependent on time. Let $\frac{d\mathbf{z}}{dt}=f(\mathbf{z}(t),t)$ be a differential equation describing a continuous-in-time transformation of $\mathbf{z}(t)$. Assuming that $f$ is uniformly Lipschitz continuous in $\mathbf{z}$ and continuous in $t$, then the change in log probability also follows a differential equation:

定理（瞬时变量替换）

设 $\mathbf{z}(t)$ 是一个有限的连续随机变量，其概率 $p(\mathbf{z}(t))$ 随时间变化。设 $\frac{d\mathbf{z}}{dt}=f(\mathbf{z}(t),t)$ 是描述 $\mathbf{z}(t)$ 随时间连续变换的微分方程。假设 $f$ 关于 $\mathbf{z}$ 一致利普希茨连续，关于 $t$ 连续，那么对数概率的变化也服从一个微分方程：

$$\frac{\partial\log p(\mathbf{z}(t))}{\partial t}=-\textnormal{tr}\left(\frac{df}{d\mathbf{z}}(t)\right)$$

Proof.

To prove this theorem, we take the infinitesimal limit of finite changes of $\log p({\mathbf{z}}(t))$ through time. First we denote the transformation of ${\mathbf{z}}$ over an $\varepsilon$ change in time as

证明

为证明这个定理，我们对 $\log p(\mathbf{z}(t))$ 随时间的有限变化取无穷小极限。首先，把 $\mathbf{z}$ 在时间变化 $\varepsilon$ 后的变换记为

$${\mathbf{z}}(t+\varepsilon)=T_{\varepsilon}({\mathbf{z}}(t)) \tag{14}$$

We assume that $f$ is Lipschitz continuous in ${\mathbf{z}}(t)$ and continuous in $t$, so every initial value problem has a unique solution by Picard’s existence theorem. We also assume ${\mathbf{z}}(t)$ is bounded. These conditions imply that $f$, $T_{\varepsilon}$, and $\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}$ are all bounded. In the following, we use these conditions to exchange limits and products.

我们假设 $f$ 关于 $\mathbf{z}(t)$ 利普希茨连续、关于 $t$ 连续，因此由皮卡存在性定理，每个初值问题都有唯一解。我们还假设 $\mathbf{z}(t)$ 有界。这些条件意味着 $f$、$T_{\varepsilon}$ 和 $\frac{\partial}{\partial\mathbf{z}}T_{\varepsilon}$ 都有界。下文用这些条件来交换极限与乘积的次序。

We can write the differential equation $\frac{\partial\log p({\mathbf{z}}(t))}{\partial t}$ using the discrete change of variables formula, and the definition of the derivative:

利用离散的变量替换公式和导数的定义，可以把微分方程 $\frac{\partial\log p(\mathbf{z}(t))}{\partial t}$ 写成：

$$\begin{aligned}\frac{\partial\log p({\mathbf{z}}(t))}{\partial t} &=\lim_{\varepsilon\rightarrow 0^{+}}\frac{\log p({\mathbf{z}}(t))-\log\left|\det\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right|-\log p({\mathbf{z}}(t))}{\varepsilon} \\ &=-\lim_{\varepsilon\rightarrow 0^{+}}\frac{\log\left|\det\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right|}{\varepsilon} \\ &=-\lim_{\varepsilon\rightarrow 0^{+}}\frac{\frac{\partial}{\partial\varepsilon}\log\left|\det\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right|}{\frac{\partial}{\partial\varepsilon}\varepsilon}\qquad\qquad\textnormal{(by L'Hôpital's rule)} \\ &=-\lim_{\varepsilon\rightarrow 0^{+}}\frac{\frac{\partial}{\partial\varepsilon}\left|\det\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right|}{\left|\det\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right|}\qquad\qquad\quad\left(\left.\frac{\partial\log({\mathbf{z}})}{\partial{\mathbf{z}}}\right|_{\mathbf{z}=1}=1\right) \\ &=-\underbrace{\left(\lim_{\varepsilon\rightarrow 0^{+}}\frac{\partial}{\partial\varepsilon}\left|\det\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right|\right)}_{\text{bounded}}\underbrace{\left(\lim_{\varepsilon\rightarrow 0^{+}}\frac{1}{\left|\det\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right|}\right)}_{=1} \\ &=-\lim_{\varepsilon\rightarrow 0^{+}}\frac{\partial}{\partial\varepsilon}\left|\det\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right|\end{aligned} \tag{15}$$

The derivative of the determinant can be expressed using Jacobi’s formula, which gives

行列式的导数可以用雅可比公式（Jacobi’s formula）表示，于是得到

$$\begin{aligned}\frac{\partial\log p({\mathbf{z}}(t))}{\partial t} &=-\lim_{\varepsilon\rightarrow 0^{+}}\textnormal{tr}\left(\text{adj}\left(\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right)\frac{\partial}{\partial\varepsilon}\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right) \\ &=-\textnormal{tr}\left(\underbrace{\left(\lim_{\varepsilon\rightarrow 0^{+}}\text{adj}\left(\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right)\right)}_{=I}\left(\lim_{\varepsilon\rightarrow 0^{+}}\frac{\partial}{\partial\varepsilon}\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right)\right) \\ &=-\textnormal{tr}\left(\lim_{\varepsilon\rightarrow 0^{+}}\frac{\partial}{\partial\varepsilon}\frac{\partial}{\partial{\mathbf{z}}}T_{\varepsilon}({\mathbf{z}}(t))\right)\end{aligned} \tag{21}$$

Substituting $T_{\varepsilon}$ with its Taylor series expansion and taking the limit, we complete the proof.

把 $T_{\varepsilon}$ 替换为它的泰勒级数展开并取极限，证明即告完成。

$$\begin{aligned}\frac{\partial\log p({\mathbf{z}}(t))}{\partial t} &=-\textnormal{tr}\left(\lim_{\varepsilon\rightarrow 0^{+}}\frac{\partial}{\partial\varepsilon}\frac{\partial}{\partial{\mathbf{z}}}\left({\mathbf{z}}+\varepsilon f({\mathbf{z}}(t),t)+\mathcal{O}(\varepsilon^{2})+\mathcal{O}(\varepsilon^{3})+\dots\right)\right) \\ &=-\textnormal{tr}\left(\lim_{\varepsilon\rightarrow 0^{+}}\frac{\partial}{\partial\varepsilon}\left(I+\frac{\partial}{\partial{\mathbf{z}}}\varepsilon f({\mathbf{z}}(t),t)+\mathcal{O}(\varepsilon^{2})+\mathcal{O}(\varepsilon^{3})+\dots\right)\right) \\ &=-\textnormal{tr}\left(\lim_{\varepsilon\rightarrow 0^{+}}\left(\frac{\partial}{\partial{\mathbf{z}}}f({\mathbf{z}}(t),t)+\mathcal{O}(\varepsilon)+\mathcal{O}(\varepsilon^{2})+\dots\right)\right) \\ &=-\textnormal{tr}\left(\frac{\partial}{\partial{\mathbf{z}}}f({\mathbf{z}}(t),t)\right)\end{aligned} \tag{24}$$

∎

A.1 Special Cases

A.1 特殊情形

Planar CNF.

Let $f({\mathbf{z}})=uh(w^{\mathbf{z}}+b)$, then $\frac{\partial f}{\partial{\mathbf{z}}}=u\frac{\partial h}{\partial{\mathbf{z}}}^{\mkern-1.5mu\mathsf{T}}{}$. Since the trace of an outer product is the inner product, we have

平面 CNF

令 $f(\mathbf{z})=uh(w^{\mathbf{z}}+b)$，则 $\frac{\partial f}{\partial\mathbf{z}}=u\frac{\partial h}{\partial\mathbf{z}}^{\mathsf{T}}$。由于外积的迹等于内积，我们有

$$\frac{\partial\log p({\mathbf{z}})}{\partial t}=-\textnormal{tr}\left(u\frac{\partial h}{\partial{\mathbf{z}}}^{\mkern-1.5mu\mathsf{T}}{}\right)=-u^{\mkern-1.5mu\mathsf{T}}{}\frac{\partial h}{\partial{\mathbf{z}}} \tag{28}$$

This is the parameterization we use in all of our experiments.

我们所有实验都使用这种参数化。

Hamiltonian CNF.

The continuous analog of NICE (Dinh et al. 2014) is a Hamiltonian flow, which splits the data into two equal partitions and is a volume-preserving transformation, implying that $\frac{\partial\log p({\mathbf{z}})}{\partial t}=0$. We can verify this. Let

哈密顿 CNF

NICE（Dinh et al. 2014）的连续版本是哈密顿流（Hamiltonian flow）。它把数据分成相等的两部分，是一种保体积变换，这意味着 $\frac{\partial\log p(\mathbf{z})}{\partial t}=0$。我们可以验证这一点。令

$$\begin{bmatrix}\frac{d{\mathbf{z}}{1:d}}{dt}\
\frac{d{\mathbf{z}}{d+1:D}}{dt}\end{bmatrix}=\begin{bmatrix}f({\mathbf{z}}{d+1:D})\
g({\mathbf{z}}{1:d})\end{bmatrix} \tag{29}$$

Then because the Jacobian is all zeros on its diagonal, the trace is zero. This is a volume-preserving flow.

由于雅可比矩阵的对角线全为零，迹为零。这是一个保体积的流。

A.2 Connection to Fokker-Planck and Liouville PDEs

A.2 与福克-普朗克方程和刘维尔偏微分方程的联系

The Fokker-Planck equation is a well-known partial differential equation (PDE) that describes the probability density function of a stochastic differential equation as it changes with time. We relate the instantaneous change of variables to the special case of Fokker-Planck with zero diffusion, the Liouville equation.

福克-普朗克方程（Fokker-Planck equation）是一个著名的偏微分方程（PDE），描述随机微分方程的概率密度函数如何随时间变化。我们把瞬时变量替换与福克-普朗克方程在零扩散时的特殊情形，即刘维尔方程（Liouville equation），联系起来。

As with the instantaneous change of variables, let $\mathbf{z}(t)\in\mathbb{R}^{D}$ evolve through time following $\frac{d\mathbf{z}(t)}{dt}=f(\mathbf{z}(t),t)$. Then Liouville equation describes the change in density of $\mathbf{z}$–a fixed point in space–as a PDE,

与瞬时变量替换中一样，设 $\mathbf{z}(t)\in\mathbb{R}^{D}$ 按 $\frac{d\mathbf{z}(t)}{dt}=f(\mathbf{z}(t),t)$ 随时间演化。刘维尔方程以 PDE 的形式描述 $\mathbf{z}$（空间中的一个固定点）处密度的变化：

$$\frac{\partial p(\mathbf{z},t)}{\partial t}=-\sum_{i=1}^{D}\frac{\partial}{\partial\mathbf{z}_{i}}\left[f_{i}(\mathbf{z},t)p(\mathbf{z},t)\right] \tag{30}$$

However, (30) cannot be easily used as it requires the partial derivatives of $\frac{p(\mathbf{z},t)}{\partial\mathbf{z}}$, which is typically approximated using finite difference. This type of PDE has its own literature on efficient and accurate simulation (Stam 1999).

然而，(30) 不容易直接使用，因为它需要偏导数 $\frac{p(\mathbf{z},t)}{\partial\mathbf{z}}$，而这通常用有限差分近似。这类 PDE 的高效、精确模拟有专门的文献（Stam 1999）。

Instead of evaluating $p(\cdot,t)$ at a fixed point, if we follow the trajectory of a particle $\mathbf{z}(t)$, we obtain

如果不在固定点上计算 $p(\cdot,t)$，而是跟随一个粒子 $\mathbf{z}(t)$ 的轨迹，就得到

$$\begin{split}\frac{\partial p(\mathbf{z}(t),t)}{\partial t}&=\underbrace{\frac{\partial p(\mathbf{z}(t),t)}{\partial\mathbf{z}(t)}\frac{\partial\mathbf{z}(t)}{\partial t}}{\text{partial derivative from first argument, }\mathbf{z}(t)}+\underbrace{\vphantom{\frac{\partial p(\mathbf{z}(t),t)}{\partial\mathbf{z}(t)}}\frac{\partial p(\mathbf{z}(t),t)}{\partial t}}{\text{partial derivative from second argument, }t}\
&=\cancel{\sum_{i=1}^{D}\frac{\partial p(\mathbf{z}(t),t)}{\partial\mathbf{z}{i}(t)}\frac{\partial\mathbf{z}{i}(t)}{\partial t}}-\sum_{i=1}^{D}\frac{\partial f_{i}(\mathbf{z}(t),t)}{\partial\mathbf{z}{i}}p(\mathbf{z}(t),t)-\cancel{\sum{i=1}^{D}f_{i}(\mathbf{z}(t),t)\frac{\partial p(\mathbf{z}(t),t)}{\partial\mathbf{z}{i}(t)}}\
&=-\sum{i=1}^{D}\frac{\partial f_{i}(\mathbf{z}(t),t)}{\partial\mathbf{z}_{i}}p(\mathbf{z}(t),t)\end{split} \tag{31}$$

We arrive at the instantaneous change of variables by taking the log,

取对数，就得到瞬时变量替换：

$$\frac{\partial\log p(\mathbf{z}(t),t)}{\partial t}=\frac{1}{p(\mathbf{z}(t),t)}\frac{\partial p(\mathbf{z}(t),t)}{\partial t}=-\sum_{i=1}^{D}\frac{\partial f_{i}(\mathbf{z}(t),t)}{\partial\mathbf{z}_{i}} \tag{32}$$

While still a PDE, (32) can be combined with $\mathbf{z}(t)$ to form an ODE of size $D+1$,

(32) 虽然仍是 PDE，但可以与 $\mathbf{z}(t)$ 合并，构成一个大小为 $D+1$ 的 ODE：

$$\frac{d}{dt}\begin{bmatrix}\mathbf{z}(t)\
\log p(\mathbf{z}(t),t)\end{bmatrix}=\begin{bmatrix}f(\mathbf{z}(t),t)\
-\sum_{i=1}^{D}\frac{\partial f_{i}(\mathbf{z}(t),t)}{\partial t}\end{bmatrix} \tag{33}$$

Compared to the Fokker-Planck and Liouville equations, the instantaneous change of variables is of more practical impact as it can be numerically solved much more easily, requiring an extra state of $D$ for following the trajectory of $\mathbf{z}(t)$. Whereas an approach based on finite difference approximation of the Liouville equation would require a grid size that is exponential in $D$.

与福克-普朗克方程和刘维尔方程相比，瞬时变量替换的实际意义更大，因为它容易数值求解得多，只需额外的 $D$ 维状态来跟踪 $\mathbf{z}(t)$ 的轨迹。而基于刘维尔方程有限差分近似的方法，所需网格大小随 $D$ 呈指数增长。

Appendix B A Modern Proof of the Adjoint Method

附录 B 伴随方法的一个现代证明

We present an alternative proof to the adjoint method (Pontryagin et al. 1962) that is short and easy to follow.

我们给出伴随方法（Pontryagin et al. 1962）的另一个证明，它简短且易于理解。

B.1 Continuous Backpropagation

B.1 连续反向传播

Let $\mathbf{z}(t)$ follow the differential equation $\frac{d\mathbf{z}(t)}{dt}=f(\mathbf{z}(t),t,\theta)$, where $\theta$ are the parameters. We will prove that if we define an adjoint state

设 $\mathbf{z}(t)$ 服从微分方程 $\frac{d\mathbf{z}(t)}{dt}=f(\mathbf{z}(t),t,\theta)$，其中 $\theta$ 是参数。我们将证明，如果定义伴随状态

$$\mathbf{a}(t)=\frac{dL}{d\mathbf{z}(t)} \tag{34}$$

then it follows the differential equation

那么它服从微分方程

$$\frac{d\mathbf{a}(t)}{dt}=-\mathbf{a}(t)\frac{\partial f(\mathbf{z}(t),t,\theta)}{\partial\mathbf{z}(t)} \tag{35}$$

For ease of notation, we denote vectors as row vectors, whereas the main text uses column vectors.

为了记号方便，这里把向量记作行向量，而正文中使用列向量。

The adjoint state is the gradient with respect to the hidden state at a specified time $t$. In standard neural networks, the gradient of a hidden layer $\mathbf{h}_{t}$ depends on the gradient from the next layer $\mathbf{h}_{t+1}$ by chain rule

伴随状态是关于指定时刻 $t$ 的隐藏状态的梯度。在标准神经网络中，由链式法则，隐藏层 $\mathbf{h}_{t}$ 的梯度依赖于下一层 $\mathbf{h}_{t+1}$ 的梯度：

$$\frac{dL}{d\mathbf{h}_{t}}=\frac{dL}{d\mathbf{h}_{t+1}}\frac{d\mathbf{h}_{t+1}}{d\mathbf{h}_{t}}. \tag{36}$$

With a continuous hidden state, we can write the transformation after an $\varepsilon$ change in time as

对于连续的隐藏状态，可以把时间变化 $\varepsilon$ 后的变换写成

$$\mathbf{z}(t+\varepsilon)=\int_{t}^{t+\varepsilon}f(\mathbf{z}(t),t,\theta)dt+\mathbf{z}(t)=T_{\varepsilon}(\mathbf{z}(t),t) \tag{37}$$

and chain rule can also be applied

链式法则同样适用：

$$\frac{dL}{\partial\mathbf{z}(t)}=\frac{dL}{d\mathbf{z}(t+\varepsilon)}\frac{d\mathbf{z}(t+\varepsilon)}{d\mathbf{z}(t)}\quad\quad\text{or}\quad\quad\mathbf{a}(t)=\mathbf{a}(t+\varepsilon)\frac{\partial T_{\varepsilon}(\mathbf{z}(t),t)}{\partial\mathbf{z}(t)} \tag{38}$$

The proof of (35) follows from the definition of derivative:

(35) 的证明可由导数的定义得到：

$$\begin{aligned}\frac{d\mathbf{a}(t)}{dt} &=\lim\limits_{\varepsilon\rightarrow 0^{+}}\frac{\mathbf{a}(t+\varepsilon)-\mathbf{a}(t)}{\varepsilon} \\ &=\lim\limits_{\varepsilon\rightarrow 0^{+}}\frac{\mathbf{a}(t+\varepsilon)-\mathbf{a}(t+\varepsilon)\frac{\partial}{\partial\mathbf{z}(t)}T_{\varepsilon}(\mathbf{z}(t))}{\varepsilon}\qquad\qquad\qquad\qquad\qquad\;\quad\;\textnormal{(by Eq 38)} \\ &=\lim\limits_{\varepsilon\rightarrow 0^{+}}\frac{\mathbf{a}(t+\varepsilon)-\mathbf{a}(t+\varepsilon)\frac{\partial}{\partial\mathbf{z}(t)}\left(\mathbf{z}(t)+\varepsilon f(\mathbf{z}(t),t,\theta)+\mathcal{O}(\varepsilon^{2})\right)}{\varepsilon}\;\;\;\quad\textnormal{(Taylor series around $\mathbf{z}(t)$)} \\ &=\lim\limits_{\varepsilon\rightarrow 0^{+}}\frac{\mathbf{a}(t+\varepsilon)-\mathbf{a}(t+\varepsilon)\left(I+\varepsilon\frac{\partial f(\mathbf{z}(t),t,\theta)}{\partial\mathbf{z}(t)}+\mathcal{O}(\varepsilon^{2})\right)}{\varepsilon} \\ &=\lim\limits_{\varepsilon\rightarrow 0^{+}}\frac{-\varepsilon\mathbf{a}(t+\varepsilon)\frac{\partial f(\mathbf{z}(t),t,\theta)}{\partial\mathbf{z}(t)}+\mathcal{O}(\varepsilon^{2})}{\varepsilon} \\ &=\lim\limits_{\varepsilon\rightarrow 0^{+}}-\mathbf{a}(t+\varepsilon)\frac{\partial f(\mathbf{z}(t),t,\theta)}{\partial\mathbf{z}(t)}+\mathcal{O}(\varepsilon) \\ &=-\mathbf{a}(t)\frac{\partial f(\mathbf{z}(t),t,\theta)}{\partial\mathbf{z}(t)}\end{aligned} \tag{39}$$

We pointed out the similarity between adjoint method and backpropagation (eq. 38). Similarly to backpropagation, ODE for the adjoint state needs to be solved backwards in time. We specify the constraint on the last time point, which is simply the gradient of the loss wrt the last time point, and can obtain the gradients with respect to the hidden state at any time, including the initial value.

我们指出了伴随方法与反向传播的相似之处（式 38）。与反向传播类似，伴随状态的 ODE 需要沿时间反向求解。我们在最后一个时间点上给定约束条件，即损失关于最后时间点的梯度，就能得到关于任意时刻隐藏状态的梯度，包括初值。

$$\underbrace{\mathbf{a}(t_{N})=\frac{dL}{d\mathbf{z}(t_{N})}}_{\text{initial condition of adjoint diffeq.}}\quad\quad\quad\quad\underbrace{\mathbf{a}(t_{0})=\mathbf{a}(t_{N})+\int_{t_{N}}^{t_{0}}\frac{d\mathbf{a}(t)}{dt}\;dt=\mathbf{a}(t_{N})-\int_{t_{N}}^{t_{0}}\mathbf{a}(t)^{T}\frac{\partial f(\mathbf{z}(t),t,\theta)}{\partial\mathbf{z}(t)}}_{\text{gradient wrt. initial value}} \tag{46}$$

Here we assumed that loss function $L$ depends only on the last time point $t_{N}$. If function $L$ depends also on intermediate time points $t_{1},t_{2},\dots,t_{N-1}$, etc., we can repeat the adjoint step for each of the intervals $[t_{N-1},t_{N}]$, $[t_{N-2},t_{N-1}]$ in the backward order and sum up the obtained gradients.