Scalable Gradients for Stochastic Differential Equations

Introduction

Deterministic dynamical systems can often be modeled by ordinary differential equations (ODEs). The adjoint sensitivity method can efficiently compute gradients of ODE solutions with constant memory cost. This method was well-known in the physics, numerical analysis, and control communities for decades . Recently, it was combined with modern reverse-mode automatic differentiation packages, enabling ODEs with millions of parameters to be fit to data and allowing more flexible density estimation and time series models .

Stochastic differential equations (SDEs) generalize ODEs, adding instantaneous noise to their dynamics . They are a natural model for phenomena governed by many small and unobserved interactions, such as motion of molecules in a liquid , allele frequencies in a gene pool , or prices in a market . Previous attempts on fitting SDEs mostly relied on methods with poor scaling properties. The pathwise approach , a form of forward-mode automatic differentiation, scales poorly in time with the number of parameters and states in the model. On the other hand, simply differentiating through the operations of an SDE solver scales poorly in memory.

In this work, we generalize the adjoint method to stochastic dynamics defined by SDEs. We give a simple and practical algorithm for fitting SDEs with tens of thousands of parameters, while allowing the use of high-order adaptive time-stepping SDE solvers. We call this approach the stochastic adjoint sensitivity method.

There are two main difficulties in generalizing the adjoint formulation for ODEs to SDEs. The first is mathematical: SDEs are defined using nonstandard integrals that usually rely on Itô calculus. The adjoint method requires solving the dynamics backwards in time from the end state. However, it is not clear exactly what “running the SDE backwards” means in the context of stochastic calculus, and when it correctly reconstructs the forward trajectory. We address this problem in Section 3, deriving a backward Stratonovich SDE whose dynamics compute the necessary gradient.

The second difficulty is computational: To retrace the steps, one needs to reconstruct the noise sampled on the forward pass, ideally without storing it. In Section 4, we give an algorithm that allows querying a Brownian motion sample at any time point arbitrarily-precisely, while only storing a single random seed.

We combine our adjoint approach with a gradient-based stochastic variational inference scheme for efficiently marginalizing over latent SDE models with arbitrary differentiable likelihoods. This model family generalizes several existing families such as latent ODEs , Gaussian state-space models , and deep Kalman filters , and can naturally handle irregularly-sampled times series and missing observations. We train latent SDEs on toy and real datasets, demonstrating competitive performance compared to existing approaches for dynamics modeling.

Background: Stochastic Flows

The adjoint sensitivity method is an efficient approach to solve control problems relying on the adjoint (co-state) system . Chen et al. used this method to compute the gradient with respect to parameters of a neural ODE, which is a particular model among many others inspired by the theory of dynamical systems . The method, shown in Algorithm 1, is scalable, since the most costly computation is a vector-Jacobian product defining its backwards dynamics. In addition, since the gradient is obtained by solving another ODE, no intermediate computation is stored as in the case of regular backpropagation .

2 Stochastic Differential Equations

3 Neural Stochastic Differential Equations

Similar to neural ODEs, one can consider drift and diffusion functions defined by neural networks, a model known as the neural SDE .

Among works on neural SDEs, none has enabled an efficient training framework. In particular, Tzen and Raginsky and Liu et al. considered computing the gradient by simulating the forward dynamics of an explicit Jacobian matrix. This Jacobian has size of either the square of the number of parameters, or the number of parameters times the number of states, building on the pathwise approach . In contrast, our approach only requires a small number of cheap vector-Jacobian products, independent of the dimension of the parameter and state vectors. These vector-Jacobian products have the same asymptotic time cost as evaluating the drift and diffusion functions, and can be easily computed by modern automatic differentiation libraries .

4 Backward Stratonovich Integral

Our stochastic adjoint sensitivity method involves stochastic processes running both forward and backward in time. The Stratonovich stochastic integral, due to its symmetry, gives nice expressions for the backward dynamics and is more convenient for our purpose. Our results can also be straightforwardly applied to Itô SDEs, relying on a simple conversion rule (see e.g. [64, Sec. 2]).

To define the backward Stratonovich integral, we consider the backward Wiener process { ^

where $\Pi=\{0=t_{N}<\cdots<t_{0}=T\}$ is the partition.

5 Stochastic Flow of Diffeomorphisms

It is well known that an ODE defines a flow of diffeomorphisms . Here we consider the stochastic analog for the Stratonovich SDE

The following theorem shows that these maps are diffeomorphisms (after choosing a suitable modification) and that they satisfy backward SDEs.

The coefficients in (52) and (87) differ by only a negative sign. This symmetry is due to our use of the Stratonovich integral (see Figure 2).

Sensitivity via Stochastic Adjoint

We present our main contribution: a stochastic analog of the adjoint sensitivity method for SDEs. We use (87) to derive another backward Stratonovich SDE, which we call the stochastic adjoint process. The direct implication is a gradient computation algorithm that works by solving a set of dynamics in reverse time, and relies on cheap vector-Jacobian products without storing any intermediate quantities.

We first present the SDE for the Jacobian matrix of the backward flow.

The proof included in Appendix 9.2 relies on Itô’s lemma in the Stratonovich form [41, Theorem 2.4.1]. We stress that this lemma considers only the case where the endpoint $z$ is fixed and deterministic.

Now, we extend to the case where the endpoint is not deterministic, but rather computed from the forward flow. To achieve this, we compose the state process and the loss function. Consider $A_{s,t}(z)=\partial\mathcal{L}(\Phi_{s,t}(z))/\partial z$. The chain rule gives $A_{s,t}(z)=\nabla\mathcal{L}(\Phi_{s,t}(z))\nabla\Phi_{s,t}(z)$. Let

Since the drift and diffusion functions of this augmented system are $C_{b}^{\infty,1}$, the system has a unique strong solution. Let $s=0$ and $t=T$. Since (166) admits a strong solution, we may write

where $W_{\cdot}=\{W_{t}\}_{0\leq t\leq T}$ denotes the path of the Wiener process and

is a deterministic measurable function (the Itô map) [69, Chapter V, Definition 10.9]. Intuitively, $\mathsf{F}$ can be thought as a black box that computes the solution to the backward SDE system (166) given the position $z$ at time $T$ and the realized Wiener process sample. Similarly, we let $\mathsf{G}$ be the solution map for the forward flow (52). The next theorem follows immediately from (156) and the definition of $\mathsf{F}$.

For $P$-almost all $\omega\in\Omega$, we have

where $\mathsf{G}(z,W_{\cdot})=\Phi_{0,T}(z)$.

This is a consequence of composing A0,T​(z)= ^

This shows that one can obtain the gradient by “composing” the backward SDE system (166) with the original forward SDE (52) and ends our continuous-time analysis.

2 Numerical Approximation

In practice, we compute solutions to SDEs with numerical solvers $\mathsf{F}_{h}$ and $\mathsf{G}_{h}$, where $h=T/L$ denotes the mesh size of a fixed grid. The approximate algorithm thus outputs $\mathsf{F}_{h}(\mathsf{G}_{h}(z,W_{\cdot}),W_{\cdot})$. The following theorem provides sufficient conditions for convergence.

Suppose the schemes $\mathsf{F}_{h}$ and $\mathsf{G}_{h}$ satisfy the following conditions: (i) $\mathsf{F}_{h}(z,W_{\cdot})\rightarrow\mathsf{F}(z,W_{\cdot})$ and $\mathsf{G}_{h}(z,W_{\cdot})\rightarrow\mathsf{G}(z,W_{\cdot})$ in probability as $h\rightarrow 0$, and (ii) for any $M>0$, we have $\sup_{|z|\leq M}|\mathsf{F}_{h}(z,W_{\cdot})-\mathsf{F}(z,W_{\cdot})|\rightarrow 0$ in probability as $h\rightarrow 0$. Then, for any starting point $z$ of the forward flow, we have

See Appendix 9.3 for the proof. Usual schemes such as the Euler-Maruyama scheme (more generally Itô-Taylor schemes) converge pathwise (i.e. almost surely) from any fixed starting point and satisfies $(i)$. While $(ii)$ is strong, we note that the SDEs considered here have smooth coefficients, and thus their solutions enjoy nice regularity properties in the starting position. Therefore, it is reasonable to expect that the corresponding numerical schemes to also behave nicely as a function of both the mesh size and the starting position. To the best of our knowledge, this property is not considered at all in the literature on numerical methods for SDEs (where the initial position is fixed), but is crucial in the proof of Theorem 3.3. In Appendix 9.4, we prove that condition $(ii)$ holds for the Euler-Maruyama scheme. Detailed analysis for other schemes is beyond the scope of this paper.

3 The Algorithm

So far we have derived the gradient of the loss with respect to the initial state. We can extend these results to give gradients with respect to parameters of the drift and diffusion functions by treating them as an additional part of the state whose dynamics has zero drift and diffusion. We summarize this in Algorithm 2, assuming access only to a black-box solver sdeint. All terms in the augmented dynamics, such as $a_{t}\partial f/\partial\theta$ and $a_{t}\partial\sigma/\partial\theta$ can be cheaply evaluated by calling $\texttt{vjp}(a_{t},f,\theta)$ and $\texttt{vjp}(a_{t},\sigma,\theta)$, respectively.

Fortunately, if we restrict our SDE to have diagonal noise, then even though the backward SDE for the stochastic adjoint will not in general have diagonal noise, it will satisfy a commutativity property . In that case, we can safely adopt certain numerical schemes of strong order 1.0 (e.g. Milstein and stochastic Runge-Kutta ) without approximating multiple integrals or the Lévy area during simulation. We formally show this in Appendix 9.5.

One may also consider numerical schemes with high weak order . However, analysis of this scenario is beyond the current scope.

4 Software and Implementation

We have implemented several common SDE solvers in PyTorch with adaptive time-stepping using a PI controller . Following torchdiffeq , we have created a user-friendly subclass of torch.autograd.Function that facilitates gradient computation using our stochastic adjoint framework for SDEs that are subclasses of torch.nn.Module. We include a short code snippet covering the main idea of the stochastic adjoint in Appendix 9.13. The complete codebase can be found at https://github.com/google-research/torchsde.

Virtual Brownian Tree

Our formulation of the adjoint can be numerically integrated efficiently, since simulating its dynamics only requires evaluating cheap vector-Jacobian products, as opposed to whole Jacobians. However, the backward-in-time nature introduces a new difficulty: The same Wiener process sample path used in the forward pass must be queried again during the backward pass. Naïvely storing Brownian motion increments implies a large memory consumption and complicates the usage of adaptive time-stepping integrators, where the evaluation times in the backward pass may be different from those in the forward pass.

To overcome this issue, we combine Brownian trees with splittable pseudorandom number generators (PRNGs) to give an algorithm that can query values of a Wiener process sample path at arbitrary times. This algorithm, which we call the virtual Brownian tree, has $\mathcal{O}(1)$ memory cost, and time cost logarithmic with respect to the inverse error tolerance.

Lévy’s Brownian bridge states that given a start time $t_{s}$ and end time $t_{e}$ along with their respective Wiener process values $w_{s}$ and $w_{e}$, the marginal of the process at time $t\in(t_{s},t_{e})$ is a normal distribution:

We can recursively apply this formula to evaluate the process at the midpoint of any two distinct timestamps where the values are already known. Constructing the whole sample path of a Wiener process in this manner results in what is known as the Brownian tree . Storing this tree would be memory-intensive, but we show how to reconstruct any node in this tree as desired.

2 Brownian Trees using Splittable Seeds

We assume access to a splittable PRNG , which has an operation split that deterministically generates two keys from an existing key. Given a key, the function BrownianBridge samples deterministically from (177). To obtain the Wiener process value at a specific time, we must first know or sample the values at the initial and terminal times. Then, the virtual Brownian tree recursively samples from the midpoint of Brownian bridges, each sample using a key split from that of its parent node. The algorithm terminates when the most recently sampled time is close enough to the desired time. We outline the full procedure in Algorithm 3.

This algorithm has constant memory cost. For a fixed-step-size solver taking $L$ steps, the tolerance that the tree will need to be queried at scales as $1/L$. Thus the per-step time complexity scales as $\log L$. Our implementation uses an efficient count-based PRNG which avoids passing large random states, and instead simply passes integers. Table 1 compares the asymptotic time complexity of this approach against existing alternatives.

Latent Stochastic Differential Equations

The algorithms presented in Sections 3 and 4 allow us to efficiently compute gradients of scalar objectives with respect to SDE parameters, letting us fit SDEs to data. This raises the question: Which loss to optimize?

Simply fitting SDE parameters to maximize likelihood will in general cause overfitting, and will result in the diffusion function going to zero. In this section, we show how to do efficient variational inference in SDE models, and optimize the marginal log-likelihood to fit both prior (hyper-)parameters and the parameters of a tractable approximate posterior over functions.

In particular, we can parameterize both a prior over functions and an approximate posterior using SDEs:

If both processes share the same diffusion function $\sigma$, then the KL divergence between them is finite (under additional mild regularity conditions; see Appendix 9.7), and can be estimated by sampling paths from the approximate posterior process. Then, the evidence lower bound (ELBO) can be written as:

and the expectation is taken over the approximate posterior process defined by (approx. post.). The likelihoods of observations $x_{1},\dots,x_{N}$ at times $t_{1},\dots,t_{N}$ depend only on latent states $z_{t}$ at corresponding times.

To compute the gradient with respect to prior parameters $\theta$ and variational parameters $\phi$, we need only augment the forward SDE with an extra scalar variable whose drift is $\frac{1}{2}|u(Z_{t},t)|^{2}$ and diffusion is zero. The backward dynamics can be derived analogously using (166). We include a detailed derivation in Appendix 9.7. Thus, a stochastic estimate of the gradients of the loss w.r.t. all parameters can be computed in a single pair of forward and backward SDE solves.

The variational parameters $\phi$ can either be optimized individually for each sequence, or if multiple time series are sharing parameters, then an encoder network can be trained to input the observations and output $\phi$. This architecture, shown in Figure 4, can be viewed as an infinite-dimensional Variational AutoEncoder (VAE) , whose latent is an SDE-induced stochastic process. We may generalize the above to cases where the diffusion is parameterized, which is then analogous to learning the prior of the latent code in VAEs.

Related Work

Gradient computation is closely related to sensitivity analysis. Computing gradients with respect to parameters of vector fields of an SDE has been extensively studied in the stochastic control literature . In particular, for low dimensional problems, this is done effectively using dynamic programming and finite differences . However, both approaches scale poorly with the dimensionality of the parameter vector.

Closely related to the current approach is the pathwise method , which is also a continuous-time analog of the reparameterization trick . Existing methods in this regime all require simulating a (forward) SDE where each step requires computing entire Jacobian matrices. This computational cost is prohibitive for high-dimensional systems with a large number of parameters.

Based on the Euler discretization, Giles and Glasserman considered simply performing reverse-mode automatic differentiation through all intermediate steps. They named this method the adjoint approach, which, by modern standards, is a form of “backpropagation through the operations of a numerical solver”. This approach, widely adopted in the field of finance for calibrating market models , has high memory cost, and relies on a fixed Euler-Maruyama discretization. Recently, this approach was also used by Hegde et al. to learn parameterized drift and diffusion functions of an SDE. In scientific computing, Innes et al. considered backpropagating through high-order implicit SDE solvers.

Ryder et al. perform variational inference over the state and parameters for Euler-discretized latent SDEs and optimize the model with backpropagation. This approach should not be confused with the formulation of variational inference for non-discretized SDEs presented in previous works and our work, as it is unclear whether the limit of their discretization corresponds to that obtained by operating with continuous-time SDEs using Girsanov’s theorem.

Our stochastic adjoint process relies on the notion of backward SDEs devised by Kunita , which is based on two-sided filtrations. This is different from the more traditional notion of backward SDEs where only a single filtration is defined . Based on the latter notion, forward-backward SDEs (FBSDEs) have been proposed to solve stochastic optimal control problems . However, simulating FBSDEs is costly due to the need to estimate conditional expectations in the backward pass .

Recent works considered the problem of inferring an approximate posterior SDE given observed data under a prior SDE with the same diffusion coefficient . The special case with constant diffusion coefficients was considered more than a decade ago . Notably, computing the KL divergence between two SDEs over a finite time horizon was well-explored in the control literature . We include background on this topic in Appendix 9.6.

Bayesian learning and parameter estimation for SDEs have a long history . Techniques which don’t require positing a variational family such as the extended Kalman filter and Markov chain Monte Carlo have been considered in the literature .

Experiments

The aim of this section is threefold. We first empirically verify our theory by comparing the gradients obtained by our stochastic adjoint framework against analytically derived gradients for problems having closed-form solutions. We then fit latent SDE models with our framework on two synthetic datasets, verifying that the variational inference framework allows learning a generative model of time series. Finally, we learn dynamics parameterized by neural networks with a latent SDE from a motion capture dataset, demonstrating competitive performance compared to existing approaches.

We report results based on an implementation of Brownian motion that stores all intermediate queries. The virtual Brownian tree allowed training with much larger batch sizes on GPUs, but was not necessary for our small-scale experiments. Notably, our adjoint approach, even when combined with the Brownian motion implementation that stores noise, was able to reduce the memory usage by $1/2$-$1/3$ compared to directly backpropagating through solver operations on the tasks we considered.

We consider three test problems (examples 1-3 from ; details in Appendix 9.8), all of which have closed-form solutions. We compare the gradient computed from simulating our stochastic adjoint process using the Milstein scheme against the exact gradient. Figure 5(a) shows that for test example 2, the error between the adjoint gradient and analytical gradient decreases with step size.

For all three test problems, the mean squared error across dimensions tends to be smaller as the absolute tolerance of the adaptive solver is reduced (e.g. see Fig. 5 (b)). However, the Number of Function Evaluations (NFEs) tends to be much larger than that in the ODE case .

Additionally, for two out of three test problems, we found that our adjoint approach with the Milstein scheme and fixed step size can be much more time-efficient than regular backpropagation through operations of the Milstein and Euler schemes (see e.g. Fig. 5(c)). Backpropagating through the Euler scheme gives gradients of higher error compared to the Milstein method. On the other hand, directly backpropagating through the Milstein solve requires evaluating high-order derivatives and can be costly.

Results for examples 1 and 3 are in Appendix 9.9.

2 Synthetic Datasets

We trained latent SDEs with our adjoint framework to recover (1) a 1D Geometric Brownian motion, and (2) a 3D stochastic Lorenz attractor process. The main objective is to verify that the learned posterior can reconstruct the training data, and that the learned priors are not deterministic. We jointly optimize the evidence lower bound (178) with respect to parameters of the prior and posterior distributions at the initial latent state $z_{0}$, the prior and posterior drift, the diffusion function, the encoder, and the decoder. We include the details of datasets and architectures in Appendix 9.10.

For the stochastic Lorenz attractor, not only is the model able to reconstruct the data well, but also the learned prior process can produce bimodal samples in both data and latent space. This is showcased in the last row of Figure 6 where the latent and data space samples cluster around two modes. This is hard to achieve using a latent ODE with a unimodal Gaussian initial approximate posterior. We include additional visualizations in Appendix 9.11.

3 Motion Capture Dataset

To demonstrate that latent SDEs can learn complex dynamics from real-world datasets, we evaluated their predictive performance on a 50-dimensional motion capture dataset. The dataset, from Gan et al. , consists of 23 walking sequences of subject 35 partitioned into 16 training, 3 validation, and 4 test sequences. We follow the preprocessing of Wang et al. .

In designing the recognition network, we follow Yıldız et al. and use a fully connected network to encode the first three observations of each sequence and thereafter predicted the remaining sequence. This encoder is chosen for fair comparison to existing models, and could be extended to a recurrent or attention model . The overall architecture is described in Appendix 9.12 and is similar to that of ODE2VAE , with a similar number of parameters. We also use a fixed step size $1/5$ of smallest interval between any two observations .

We train latent ODE and latent SDE models with the Adam optimizer and its default hyperparameter settings, with an initial learning rate of $0.01$ that is exponentially decayed with rate $0.999$ during each iteration. We perform validation over the number of training iterations, KL penalty , and KL annealing schedule. All models were trained for at most $400$ iterations, where we start to observe severe overfitting for most model instances. We report the test MSE on future observations following Yıldız et al. . We believe that the improved performance is due to the strong regularization in path space, as removing the KL penalty improve training error but caused validation error to deteriorate.

Discussion

We presented a generalization of the adjoint sensitivity method to compute gradients through solutions of SDEs. In contrast to existing approaches, this method has nearly the same time and memory complexity as simply solving the SDE. We showed how our stochastic adjoint framework can be combined with a gradient-based stochastic variational inference scheme for training latent SDEs.

It is worthwhile to mention that SDEs and the commonly used GP models define two distinct classes of stochastic processes, albeit having a nonempty intersection (e.g. Ornstein-Uhlenbeck processes fall under both). Computationally, the cost of fitting GPs lies in the matrix inversion, whereas the computational bottleneck of training SDEs is the sequential numerical solve. Empirically, another avenue of research is to reduce the variance of gradient estimates. In the future, we may adopt techniques such as control variates or antithetic paths.

On the application side, our method opens up a broad set of opportunities for fitting any differentiable SDE model, such as Wright-Fisher models with selection and mutation parameters , derivative pricing models in finance, or infinitely-deep Bayesian neural networks . In addition, the latent SDE model enabled by our framework can be extended to include domain knowledge and structural or stationarity constraints in the prior process for specific applications.

On the theory side, there remain fundamental questions to be answered. Convergence rates of numerical gradients estimated with general schemes are unknown. Additionally, since our analyses are based on strong orders of schemes, it is natural to question whether convergence results still hold when we consider weak errors, and moreover if the method could be reformulated more coherently with rough paths theory .

We thank Yulia Rubanova, Danijar Hafner, Mufan Li, Shengyang Sun, Kenneth R. Jackson, Simo Särkkä, Daniel Lacker, and Philippe Casgrain for helpful discussions. We thank Çağatay Yıldız for helpful discussions regarding evaluation settings of the mocap task. We also thank Guodong Zhang, Kevin Swersky, Chris Rackauckas, and members of the Vector Institute for helpful comments on an early draft of this paper.

References

Appendix

2 Proof of Theorem 3.1

3 Proof of Theorem 3.3

We show that both $I_{h}^{(1)}$ and $I_{h}^{(2)}$ converge to in probability as $h\to 0$. For simplicity, we suppress $z$ and $W_{\cdot}$.

Bounding $I_{h}^{(1)}$. Let $\epsilon>0$ be given. Since $G_{h}\rightarrow G$ in probability, there exist $M_{1}>0$ and $h_{0}>0$ such that

Now, suppose $h\leq\min\{h_{0},h_{1}\}$. Then, by the union bound, with probability at least $1-4\epsilon$, we have

Thus, we have shown that $I_{h}^{(1)}$ converges to in probability as $h\rightarrow 0$.

Bounding $I_{h}^{(2)}$. The idea is similar. By condition (ii), we have

in probability. Using this and condition (i), for given $\epsilon>0$, there exist $M>0$ and $h_{2}>0$ such that for all $h\leq h_{2}$, we have

with probability at least $1-\epsilon$. On this event, we have

Thus $I_{h}^{(2)}$ also converges to in probability as $h\to 0$. ∎

4 Euler-Maruyama Scheme Satisfies Local Uniform Convergence

Here we verify that the Euler-Maruyama scheme satisfies condition $(ii)$ when $d=1$. Our proof can be extended to the case where $d>1$ assuming an $L^{p}$ estimate of the error; see the discussion after the proof of Proposition 9.1.

Usual convergence results in stochastic numerics only control the error for a single fixed starting point. Here, we strengthen the result to local uniform convergence. Our main idea is to apply a Sobolev inequality argument [54, Part II]. To do so, we need some preliminary results about the Euler-Maruyama discretization of the original SDE and its derivative. We first recall a theorem characterizing the expected squared error for general schemes.

for a constant $C$ that does not depend on $h$ or $z$.

We refer the reader to for the precise definitions of orders of accuracy and the proof. Given this theorem, we establish an estimate regarding errors of the discretization and its derivative with respect to the initial position.

where $C_{1}$ is a constant independent of $z$ and $h$.

Since the coefficients of the SDE are of class $C_{b}^{\infty,1}$, we may differentiate the SDE in $z$ to get the SDE for the derivative $\nabla_{z}Z^{z}_{t}$ . Specifically, letting $Y_{t}^{z}=\nabla_{z}Z_{t}^{z}$, we have

Note that the augmented process $(\mathsf{F}(z),\nabla_{z}\mathsf{F}(z))$ satisfies an SDE with $C_{b}^{\infty,1}$ coefficients. By the chain rule, one can easily show that the derivative of the Euler-Maruyama discretization $\mathsf{F}_{h}(z)$ is the discretization of the derivative process $Y_{t}^{z}$. Thus, $(\mathsf{F}_{h}(z),\nabla_{z}\mathsf{F}_{h}(z))$ is simply the discretization of $(\mathsf{F}(z),\nabla_{z}\mathsf{F}(z))$.

Since the Euler-Maruyama scheme has orders of accuracy $(p_{1},p_{2})=(1.5,1.0)$ [51, Section 1.1.5], by Theorem 9.2, we have

for some constant $C_{1}$ that does not depend on $z$ or $h$. ∎

We also recall a variant of the Sobolev inequality which we will apply for $d=1$.

For any $p>d$, there exists a universal constant $c_{p}$ such that

where $\alpha>1/2$ is a fixed constant. Since $\mathsf{H}_{h}^{\alpha}$ is continuously differentiable a.s., by Theorem 9.4,

Without loss of generality, we may let the compact set be $A=\{z:|z|\leq M\}$ where $M>0$. Then,

It remains to estimate $\left\|\mathsf{H}_{h}^{\alpha}\right\|_{1,2}$. Starting from the definition of $\left\|\cdot\right\|_{1,p}$, a standard estimation yields

where $C_{2}$ is a deterministic constant depending only on $\alpha$ (but not $z$ and $h$).

Now we take expectation on both sides. By Lemma 9.3, we have

where the last integral is finite since $\alpha>1/2$.

Generically extending the argument to other numerical schemes, however, is technically non-trivial. We plan to address this question in future research.

5 Stochastic Adjoint has Commutative Noise when Original SDE has Diagonal Noise

By definitions of $f$ and $g_{i}$, the Jacobian matrices $\nabla f(x,s)$ and $\nabla g_{i}(x,s)$ have the forms

Thus, the backward Stratonovich SDEs for the adjoint processes of the state and parameters have the forms

Now assume the original SDE has diagonal noise. Then, $m=d$ and Jacobian matrix $\nabla\sigma_{i}(z)$ has the form

Consider the adjoint process for the augmented state along with the backward flow of the backward Stratonovich SDE (87), whose overall state we denote by ^

for all $j_{1},j_{2}\in[m]$ and $k\in[d]$. When an SDE has commutative noise, the computationally intensive double Itô integrals (and the Lévy areas) need not be simulated by having the numerical scheme take advantage of the following property of iterated integrals :

where the Brownian motion increment $\Delta W^{(i)}=W_{t}^{(i)}-W_{s}^{(i)}$ for $i\in[m]$ can be easily sampled.

We show the diffusion function (376) satisfies the commutativity condition (377) with a proof by exhaustion:

Both LHS and RHS are zero unless $j_{1}=j_{2}=k$, since for $\Sigma_{i,j_{2}}(x)\frac{\partial\Sigma_{k,j_{1}}(x)}{\partial x_{i}}$ to be non-zero, $i=j_{1}=j_{2}=k$.

𝑑1…2𝑑k=d+1\dots,2d. Similar to the case above.

2𝑑1…2𝑑𝑝k=2d+1\dots,2d+p. Write $k=2d+l$, where $l\in[p]$. Both LHS and RHS are zero unless $j_{1}=j_{2}=l$, since for $\Sigma_{i,j_{2}}(x)\frac{\partial\Sigma_{k,j_{1}}(x)}{\partial x_{i}}$ to be non-zero $i=l$ or $i=d+l$ and $j_{1}=j_{2}=l$.

This concludes that the commutativity condition holds. Finally, we comment that the Milstein scheme for the stochastic adjoint of diagonal noise SDEs can be implemented such that during each iteration of the backward solve, vjp is only called a number of times independent of the dimensionality of the original SDE.

6 Background on Latent SDE

Recall the approximate posterior process that we intend to learn is governed by the SDE:

Suppose there exists a measurable function $u(z,t)$ such that

$\sigma(z,t)u(z,t)=h_{\phi}(z,t)-h_{\theta}(z,t)$, and

Novikov’s condition ensures that the process

Moreover, since a simple rewrite shows that

we conclude that the $Q$-law of (382) (or equivalently (379)) is the same as the $P$-law of the prior process.

Let $x_{t_{1}},\dots,x_{t_{N}}$ be observed data at times $t_{1},\dots,t_{N}$, whose conditionals only depend on the respective latent states $z_{t_{1}},\dots,z_{t_{N}}$. Since the $Q$-law of the approximate posterior is the same as the $P$-law of the prior,

where the second line follows from the definition of $Q$ and third line follows from Jensen’s inequality. In the last equality we used the fact that the Itô integral $\int_{0}^{\cdot}u(Z_{t})^{\top}dW_{t}$ is a martingale.

7 Stochastic Adjoint for Latent SDE

To simulate the variational lower bound (178) with Monte Carlo in the forward pass along with the original dynamics, we need only extend the original augmented state with an extra variable $L_{t}$ such that the new drift and diffusion functions for the new augmented state $Y_{t}=(Z_{t},\theta,L_{t})$ are

By (166), the backward SDEs of the adjoint processes become

In this case, neither does one need to simulate the backward SDE of the extra variable nor does one need to simulate its adjoint. Moreover, when considered as a single system for the augmented adjoint state, the diffusion function of the backward SDE (389) satisfies the commutativity property (377).

8 Test Problems

In the following, $\alpha,\beta$, and $p$ are parameters of SDEs, and $x_{0}$ is a fixed initial value.

In each numerical experiment, we duplicate the equation 10 times to obtain a system of SDEs where each dimension had their own parameter values sampled from the standard Gaussian distribution and then passed through a sigmoid to ensure positivity. Moreover, we also sample the initial value for each dimension from a Gaussian distribution.

9 Results for Example 1 and 3

10 Toy Datasets Configuration

Consider a geometric Brownian motion SDE:

We use $\mu=1$, $\sigma=0.5$, and $x_{0}=0.1+\epsilon$ as the ground-truth model, where $\epsilon\sim\mathcal{N}(0,0.03^{2})$. We sample $1024$ time series, each of which is observed at intervals of 0.02 from time 0 to time 1. We corrupt this data using Gaussian noise with mean zero and standard deviation $0.01$.

To recover the dynamics, we use a GRU-based latent SDE model where the GRU has 1 layer and 100 hidden units, the prior and posterior drift functions are MLPs with 1 hidden layer of 100 units, and the diffusion function is an MLP with 1 hidden layer of 100 hidden units and the sigmoid activation applied at the end. The drift function in the posterior is time-inhomogenous in the sense that it takes in a context vector of size 1 at each observation that is output by the GRU from running backwards after processing all future observations. The decoder is a linear mapping from a 4 dimensional latent space to observation space. For all nonlinearities, we use the softplus function. We fix the observation model to be Gaussian with noise standard deviation $0.01$.

We optimize the model jointly with respect to the parameters of a Gaussian distribution for initial latent state distribution, the prior and posterior drift functions, the diffusion function, the GRU encoder, and the decoder. We use a fixed discretization with step size of $0.01$ in both the forward and backward pass. We use the Adam optimizer with an initial learning rate of $0.01$ that is decay by a factor of $0.999$ after each iteration. We use a linear KL annealing schedule over the first 50 iterations.

10.2 Stochastic Lorenz Attractor

Consider a stochastic Lorenz attractor SDE with diagonal noise:

We use $\sigma=10$, $\rho=28$, $\beta=8/3$, $(\alpha_{x},\alpha_{y},\alpha_{z})=(.15,.15.,.15)$, and $(x_{0},y_{0},z_{0})$ sampled from the standard Gaussian distribution as the ground-truth model. We sample $1024$ time series, each of which is observed at intervals of $0.025$ from time 0 to time 1. We normalize these samples by their mean and standard deviation across each dimension and corrupt this data by Gaussian noise with mean zero and standard deviation $0.01$.

We use the same architecture and training procedure for the latent SDE model as in the geometric Brownian motion section, except that the diffusion function consists of four small neural networks, each for a single dimension of the latent SDE.

11 Additional Visualization

See Figure 8 for additional visualization on the synthetic Lorenz attractor dataset. See Figure 9 for visualization on the synthetic geometric Brownian motion dataset. We comment that for the second example, the posterior reconstructs the data well, and the prior process exhibit behavior of the data. However, from the third row, we can observe that the prior process is learned such that most of the uncertainty is account for in the initial latent state. We leave the investigation of more interpretable prior process for future work.

12 Model Architecture for Learning from Motion Capture Dataset

We use a latent SDE model with an MLP encoder which takes in the first three frames and outputs the mean and log-variance of the variational distribution of the initial latent state and a context vector. The decoder has a similar architecture as that for the ODE2VAE model and projects the $6$-dimensional latent state into the $50$-dimensional observation space. The posterior drift function takes in a $3$-dimensional context vector output by the encoder and the current state and time, whereas the prior drift only takes in the current state and time. The diffusion function is composed of multiple small neural nets, each producing a scalar for the corresponding dimension such that the posterior SDE has diagonal noise. We use the same observation likelihood as that of the ODE2VAE model . We comment that the overall parameter count of our model (11605) is smaller than that of ODE2VAE for the same task (12157).

The latent ODE baseline was implemented with a similar architecture, except is does not have the diffusion and prior drift components, and its vector field defining the ODE does not take in a context vector. Therefore, the model has slightly fewer parameters (10573) than the latent SDE model. See Figure 10 for overall details of the architecture.

The main hyperparameter we tuned was the coefficient for reweighting the KL. For both the latent ODE and SDE, we considered training the model with a reweighting coefficient in $\{1,0.1,0.01,0.001\}$, either with or without a linear KL annealing schedule that increased from to the prescribed value over the first $200$ iterations of training.

13 Stochastic Adjoint Implementation

We include the core implementation of the stochastic adjoint, assuming access to a callable Brownian motion bm, an Euler-Maruyama integrator ito_int_diag for diagonal noise SDEs, and several helper functions whose purposes can be inferred from their names.