Robust and Scalable SDE Learning: A Functional Perspective
Authors: Scott Alexander Cameron
Tyron Luke Cameron
Arnu Pretorius
Stephen J. Roberts
Venue: ICLR 2022 Poster (2022)
Keywords: SDE Learning
Parallelization
Importance Sampling

Abstract: Stochastic differential equations provide a rich class of flexible generative
models, capable of describing a wide range of spatio-temporal processes. A host
of recent work looks to learn data-representing SDEs, using neural networks and
other flexible function approximators. Despite these advances, learning remains
computationally expensive due to the sequential nature of SDE integrators. In
this work, we propose an importance-sampling estimator for probabilities of
observations of SDEs for the purposes of learning. Crucially, the approach we
suggest does not rely on such integrators. The proposed method produces
lower-variance gradient estimates compared to algorithms based on SDE
integrators and has the added advantage of being embarrassingly parallelizable.
This facilitates the effective use of large-scale parallel hardware for massive
decreases in computation time.


ROBUST AND SCALABLE SDE LEARNING: A FUNCTIONAL PERSPECTIVE
Scott Cameron
Tyron Cameron Discovery Insure South Africa
Arnu Pretorius InstaDeep Ltd. South Africa
Stephen Roberts Oxford University United Kingdom
Oxford University, InstaDeep Ltd. United Kingdom
S . cameron@ nst adeep. com
ABSTRACT
Stochastic differential equations provide rich class of flexible generative models, capable of describing a wide range of spatio-temporal processes. A host of re- cent work looks learn data-representing SDEs, using neural networks and other flexible function approximators. Despite these advances. learning remains compu- tationally expensive due the sequential nature of SDE integrators. In this work, we propose an iepootance-samplingg estimator for probabilities of observations of SDEs for the purposes of learning. Crucially, the approach we suggest does not rely such integrators. The proposed method produces lower-variance gradient estimates compared to algorithms based on SDE integrators and has the added advantage of being embarrassingly parallelizable. This facilitates the effective use large-scale parallel hardware for massive decreases in computation time.
INTRODUCTION
Stochastic differential equations (SDEs) natural extension to ordinary differential equations which allows modelling of noisy and uncertain driving forces. These models are particularly appealing due their flexibility in expressing highly complex relationships with simple equations, while retaining high degree of interpretability. Much work has been done over the last century focussing on understanding and modelling with SDEs, particularly in dynamical systems and quantitative finance (Pavliotis, 2014; Malliavin & Thalmaier, 2006). Examples of such work are the Langevin model of stochastic dynamics of particles in fluid, stochastic turbulence (Kumar et al., 2020), and the Black- Seccl-es model. More recently. SDE models have gathered interest in machine learning and current work focuses on efficient and scalable methods of inferring SDEs from data.
1.1 THE SDE LEARNING PROBLEM
Consider an Itō SDE (Pavliotis, 2014) of the following form
dX f(X,t) + g(X,t)
(I)
and of observations F1.N ((t,, n=1² In the simplest case, these observations would be directly realisation of the process at discrete points in time. In many physical applications, the observations themselves would have some uncertainty associated with them, described by observation noise model. In more complex models, the SDE would be driving some unobserved parameters describing further complex processes (Li et al., 2020). For example, discrete events may be modeled by an inhomogenious Poisson processes, where the rate parameter is evolving stochastically over time. In any case, the generating process can be summarized by the graphical model in the following diagram and the corresponding learning problem can be described as: "given observations Ni infer f and g".
W
X
N

---

For the most part, the model can be trained using standard techniques, such variational inference, Markov chain Monte Carlo or optimizing an appropriate loss function. Generally, this requires estimating an expectation over comple--pths,, either to estimate the marginal probability of the obser- vations as a function of f and g, or to estimate the mismatch of sample-paths with the observations. To estimate these expectations, most current work makes use of integration algorithms. These algorithms generate trajectories of X for a given f and which can be used to estimate the required expectation. For example, the Euler-Maruyama algorithm (Li et al., 2020; Jia & Benson, 2019; Liu al., 2019) generates trajectories with the following update rule
X1+ar + Xx+ + f(X,,t) At + gXx,,, 8) At).
(2)
The expectation of some functional F estimated using the sample mean of its value over many independent trajectories (A)
Unfortunately, integrators such Equation (2) introduce sequential dependence between each point and the next, which inhibits parallelization. Furthermore, it introduces difficulties in probabilistic modelling of the functions I' and g: order to properly estimate the expectation, the sample paths should be generated under consistent realisations of f and g. Models which implicitly describe f as Gaussian process would need to sample f(X,, t) conditioned on all the previously sampled values of f along the trajectory, which may be computationally prohibitive.
The integrator generates trajectories sampled from the prior distribution over paths. In the initial stages of learning, the path-space prior distribution usually has minimal overlap with the path-space posterior and typical trajectories do not pass anywhere near the observations, resulting in high variance gradients. This grows considerably worse the dimension of the problem and the time between observations increases.
The central idea of this paper is to instead generate trajectories from distribution much closer to the path-space posterior distribution when estimating these expectations. In the simplest case where the observations are directly of the process X, this means sampling trajectories which pass directly through the observations. Alternatively, one can marginalize out the observation noise distribution or efficiently estimate some expected loss function - by importance sampling
p(ZEN ( (I:N) ((FI:) px (
Ex [p(EEN x))
(4)
where q is an importance distribution which may depend on the observations themselves. For large class of observation distributions, such as additive Gaussian noise multinomial observations, natural choice for would be the corresponding conjugate distribution, in which case the importance weight p/q becomes a constant. In complex graphical models, might be parameterized and optimized variationally, or Gibbs sampling may be used with a collection of simpler conditional distributions. Taking advantage of this, shall, without loss of generality, concentrate only estimating the quantity Px (ZI:N) = Ex II, 8(zn ((.)))
In this work, we introduce a novel algorithm for efficient estimation of probabilities in SDE models for the purposes of learning. Our algorithm has a number of favourable charateristics:
Our algorithm exhibits lower gradient variance in estimated expectations than algorithms based on SDE solvers, since it enforces that sampled trajectories pass through the observed data, producing accurate results even in hundreds of dimensions.
Since our algorithm does not rely on an SDE solver, it completely removes all sequential dependence between sampled points. Our algorithm is therefore trivially parallelizable, and can take full advantage of modern parallel hardware such as multi-GPU training.
Since estimator can be calculated in single forward pass of the function f. proba- bilistic approaches such as variational Gaussian processes are viable candidates for the representation of f.
2

---

1.2 RELATED WORK
Li et al. (2020) propose a method for calculating gradients of functions of the integrated process with respect to parameters of f and g via a backwards SDE called the adjoint method. This approach requires reconstructing the Wiener process backwards in time when integrating the backward equation and the authors propose a memory efficient algorithm to do so. Their approach uses a constant amount of memory independent of the number of integration steps similar to the adjoint method for NeuralODEs (Chen et al., 201 Tzen Raginsky (2019) discuss the computation of gradients by simulation of forward SDE in the derivatives with respect to parameters of the model using the chain rule of stochastic calculus. This approach does not require proppppopgggting through the SDE solver. Kidger et al. (2021) propose method of learning SDEs by treating the SDE as the generator in a generative adversarial network setup, using a second SDE for the discriminator, All of these approaches use SDE integration algorithms to estimate path-space expectations.
Batz et al. (2018) consider a nonparametric approach to SDE learning, using Gaussian processes for the drift and diffusion functions. They initially consider gradient matching approximation, in which the drift Gaussian process may be fit directly by conventional methods, and thereafter use an appctation-maxiiizztiio algorithm, and estimate the transition probability for sparse observations using a linear approximation of the process. Yildiz et al. (2018) propose an alternative approach to Gaussian process-based SDE learning using the Marer-Maruyama integrator, approximating / by the predictive mean of a Gaussian process conditioned on a MAP estimate for set of inducing points. Unfortunately, this approach completely ignores any uncertainty in the posterior over f.
It should be noted that each of the above mentioned studies rely on SDE integrators to generate trajectories and therefore face the same difficulties associated with them. Furthermore, almost all of the experiments therein are on system of dimension between and 3, with a small number of experiments of dimension as high as We are unaware of any work which trains neural SDEs in higher dimensions, and we believe this is testament to the inherent difficulty of learning SDEs in high dimensions.
Much other work focuses on learning specific classes of SDE or ODE models, such as symplectic or Hamiltonian systems (Zhong et al., 2020; Greydanus et al., 2019), graph-based models (Poli et al., 2021; Sanchez-Gonzaler et al., 2019), and controlled differential equations (Morrill et al., 202 I). Many of these ideas are independent of the learning algorithm and can potentially make use of our algorithm when applied to SDEs.
2 PATH-SPACE IMPORTANCE SAMPLING
It is well known that solutions of linear SDEs are Gaussian processes (see Pavliotis, 2014, Section 3.7). For these processes, the joint probability density over some finite collection points is available in closed form. Furthermore, the posterior distribution of the process conditioned on some set of observations can be sampled from exactly. Unfortunately, this not the case in general for non-linear SDEs, and one has to resort to approximation methods in order to calculate probabilities. Most commonly, one uses an SDE integrator to estimate expectations; however, other methods, such as Laplace approximations and perturbation theory are sometimes used (Karimi & McAuley, 2016; Dass et al., 2017; Hutzenthaler & Jentzen, 2020).
Alternatively, if we are able to find a second process which follows similar SDE to the of interest, it is possible to express expectations and probabilities as importance sampling estimators. Linear SDEs are an attractive candidate for this due to their Gaussianity and ease of sampling.
2.1 STATE-INDEPENDENT DIFFUSION
For now, consider the simpler scenario in which the diffusion coefficient does not depend on the state variable; i.e. g(x t) = att for some function σ. The process defined by
dY = ol)
(5)
is Gaussian with mean zero and conditional variance
Ey (x((t)) - (o(t))-((to)))) = f a(t)oT(t)dd
(6)
3

---

Sample paths of Y can be generated efficiently by simulating Brownian motion. The process Y conditioned on a set of observations of the process a Brownian bridge. If we are able to express expectations under this process instead of our original SDE, then we can
quantities of interest efficiently estimate such quantities via Monte Carlo sampling. As it turns out, expectations and probabilities for SDE with general drift term can in fact be expressed as such. This relation is given in the following theorem.
Theorem Let X and Y be the stochastic processes generated by the following SDEs dX = f(X,t) t) dt o(t) dW,
(7)
dY o() dW.
(8)
Further assume that (t) Lebesgue-almost-er -aaggu---llnotteeewwwwhrr behaved such that X (t) has unique probubility density for all t.
everywhere invertible and f is sufficiently well-
The probability density of observations {()) N n=1" under the process is given by the condi- tional expectation
Px(11NN) pv (I [MMM] (Y =In)==1
(9)
where
1 (0)) ((Y, S[V] =|== / (( dr,
(10)
2
This result follows from Girsanov's 's theorem (Girsanov, 1960; Malliavin Thalmaier, 2006, Sec- tion 1.5) and the definition of conditional expectation. Intuitively speaking, the first term in Equa- tion 10 encourages f to line up with typical paths which pass through the observations, while the second term regularizes the magnitude of
Theorem allows us to develop an importance sampling algorithm for estimating probabilities under the process X by simulating the process Y, where the importance weights are given by () In this approach, one can generate entire sample paths of before calculating f(Y,t), allowing the forward pass of f to be performed in parallel. This approach is described in Algorithm 1.
Algorithm Path Integral Importance Sampling
1: for L,....A K do
2:
Sample (Yr) s.t. for each n, Yt, 2n
Y is Brownian bridge
3:
ft - f(Y,, t) for each
4:
a Σ, 2 55 (()) (Y++At Y))
5:
Σ1 fT
At
6:
S: + 3/3
7: end for
8: return PY (E 1/2 exp(
An important property of this algorithm is that the sampled paths directly pass through observations. This means that the probability estimates directly relate the drift function to plausible trajectories instead of allowing the particles to wander around randomly. See Figure for an illustration. In the initial stages of training, simply integrating the SDE with drift f will typically not lead to trajectories which pass near the data, especially in high dimensions. This can have significant impact on the variance of the gradients and can inhibit the model's ability learn. Training paths that pass through the data can greatly mitigate this high variance, As with other algorithms that generate SDE sample paths, the discretization of the integrals introduces small bias. However, this bias vanishes in the At + 0 limit.
2.1.1 A NOTE ON PARALLELIZABIIITY
With the exception of the sampling step in Line of Algorithm 1, each of the other operations, including the sums, be performed in parallel over feature and batch dimensions, well
4

---

obsanations
(a) SDE integration
(b) Brownian bridge importance samples
Independently sampled paths. shows paths sampled from the SDE using an integrator, while b shows samples from a Brownian bridge which pass exactly through the observations. In both cases the drift function is given by neural network which is randomly initialized with the same seed. The drift vector field is shown in grey.
Figure 1:
the & and indices. Sampling a Brownian bridge, or any conditional Gauss-Markov process, can be performed by first sampling the unconditional process, and then linearly interpolating with the observations. The linear interpolation can performed independently, and therefore in parallel, over the time dimension. The only operation that is not completely parallelizable is the sampling of the unconditional base process in this case standard Brownian motion. Brownian paths can be sampled by cumulative summation of independent normal random numbers. However, this cumulative sum is extremely fast, and can be performed independently between each consecutive pair of observations. In our tests, sampling whole Brownian trajectories was about orders of magnitude faster than the forward pass of the neural network f, and so does not create any noticeable performance impact.
Once the Brownian paths sampled, the rest of Algorithm 1 may be implemented in parallel without difficulty. Contrast this to standard integrator-bar base methods; even if the Brownian paths are sampled upfront, the integrator must still perform each of the forward passes of in sequence.
See Appendix for more details on sampling Gaussian bridges.
2.2 STATE-DEPENDENT DIFFUSION
The assumption that g independent of the is not in fact required for the validity of Theorem the reason this assumption is required is to ensure that the process Y is Gaussian, so as to enable easy conditional sampling and calculation of py (2 Unfortunately, one cannot simply find an analogue of Theorem to calculate expectations of general process expectations with respect to a constant diffusion process.
To address the question of state-dependent diffusion, it is enlightening to consider simple example of such process. Perhaps the most common example is geometric Brownian motion
dX = ,xX di o dW.
(11)
The typical way in which one would solve this is to introduce transformation = log( ()) and apply Itô's lemma (see Malliavin & Thalmaier, 2006, Section 1.5; Pavliotis, 2014, Section 3.5) to obtain the SDE
dz (μ ) dt ddW.
(12)
Similarly, one might address a problem of a general state-dependent diffusion coefficient by consider- ing a transformation be a map which is invertible for all t. That is, there exists a function 7 such that if T(x, t), then æ (g. t). Transforming process with constant diffusion o,, by the map T-1 gives process with diffusion coefficient
space in which the transformed equation has constant diffusion. Let T(,
OT,1 Oyj ōk.j
(ZT
(x,
Okii
(13)
dr
ik
5

---

where repeated indices imply summation. With such transformation matrix that can be written as a Jacobian. While the restrictions on the diffusion matrix may be seen as an inherent limitation of a path-space importance sampling approach, it may be argued that constant or simple diffusion functions would cover the majority of cases of interest, and almost all real-world problems require
represent any diffusion
form of data pre-processing and, potentially, link functions.
This transformation allows us to infer the SDE transformed space, with constant diffusion, and, if desired, transform the equation back to the original space using Itô's lemma. This approach is described in Algorithm 2, where
have used the symbol f to represent the drift in the transformed space as to avoid ambiguity with the drift in the original space. One important point to note is that one should be careful allow gradients propagate through the sampling step on Line in Algorithm so to allow learning of the function T.
Algorithm 2 Transformed-Staae Path Integral Importance Sampling
1: for = 1....,K K do
2:
Sample (Y:)
for each n, Yr,, ((ras
Y is Brownian bridge
3:
fe - f(Y,) for each
4:
a+ E, PT 0 (Y++ii - Yi)
5:
8 E, St S
6:
S 1/3
7:
end for
8: return
II, det ((((ta))))
py ((JINN,IIN)) K E, exp(S.)
Ozn
One would typically apply this method infer rather than the original f for performance reasons, and then only reconstruct the SDE in the original space if required.
2.2.1 RECONSTRUCTION OF THE SDE IN THE UNTRANSFORMED SPACE
Generation of new sample paths can be performed easily in the transformed space and then simply transforming back by applying 7-1, However, sometimes one needs more information for the purposes of analysing the learned equation. In this case the SDE can be transformed back using Itô's lemma. The drift and diffusion coefficients are
0T,11
Ø.. (v. + 2 2 Ōy,ūy1 01ks
8T.
[((x,
+
(14)
Øt
dyk
ÖT
Ii, (I,
Ok.j
(15)
Qyk
All expression evaluated at y T(x, t). When T' is parameterized by a neural network with ReLU non-linearisies, the second derivative term in Equation (14) is almost-everywhere zero. The remaining terms are Jacobian vector products, and be performed efficiently even in high dimensions."
3 EXPERIMENTS
In this section, we provide some experiments to illustrate and verify the performance and capabilities of our proposed algorithm. Further details on the experiments are given Appendix
3.1 THE LORENZ SYSTEM
To compare the effect of our algorithm on learning, we follow Li al. (2020) and train neural network based SDE model on the Lorenz system. In this model the drift function is given by
However, this may require forward-mode instead of backward-mode differentiation.

---

multi-layer perceptron, and the diffusion term is simple scalar. We generated 16 independent paths, consisting of 200 observations each, uniformly sampled over a time interval of [0, 10]. We trained the model on this data by two different methods. In the first method, we minimize the mean-squared error, estimated by integrating out the SDE, and calculating gradients using the adjoint method (Li et al., 2020).² In the second method, we directly optimize the log-probability calculated using Algorithm 1. In both cases we used precisely the same optimization parameters, the same time-step size At, and the same number of sampled trajectories calculate averages. Both models used the same seed for the neural network initialization, and the same number of gradient steps were performed in both cases. These models were trained on an Nvidia RTX 3070, with 8Gb of memory.
To assess the performance of these approaches. we record two metrics. The first is the mean-square difference between the learned drift vector field, and the Lorenz vector field, normalized by the mean- square magnitude of the Lorenz vector field, measured at the observations. Written mathematically:
En - [Lor()n)) 2
(16)
En (fLo
This metric is shown in Figure 2b. The second metric record is the mean-squared error of sampled trajectories compared to the observations; this is shown in Figure 2a. However, we note that this is not fair performance metric since the adjoint-based method is directly optimizing for it. On the x-axis we show the training time in minutes.
4,0
path an adjoint
path Int adjsint
3.5
400-
3,0
2.5
531 moc
2.0
1.5
100
0.0
00
100
120
20
60
80
100
120
time (minutes)
time (mirutes)
(a)
(b)
Figure 2: Learning curves of neural SDE on the Lorenz system. shows the normalized mean- squared deviation of the learned drift from the ground truth. a shows the mean- sean-squared error of sample paths from the observations.
In Figure 2 we have limited the range of the x-axis in order for the lines to be more clearly visible. The total training time for the -dseed MSE optimization was factor 52 longer than the time taken optimizing the probability estimates. We note that, while our approach was significantly faster, it did use approximately four times as much memory as the adjoint-basee method."
The poor performance of the integrator- based learning algorithm in Figure 2b despite its low mean- squared-error suggests that directly optimizing for the discrepancy of typical SDE sample-paths compared to observations not a robust approach to SDE learning.
3.2 GRADIENT VARIANCE
We applied both our algorithm and an SDE integrator to neural network- SDE model on a real-world dataset and recorded the gradients, sampled over many independent trajectories. We thereafter computed the variance of the gradients individually per parameter, for various time lengths between observations. We used discretization step size of 0.01. These results are shown in Figure 3.
We also attempted training model without the adjoint method and instead directly hackpropagating through the SDE integrator. However, a single forward pass required more memory than available on the GPU.
'The adjoint algorithm was specifically designed to have constant memory complexity in the number of integration steps. Our algorithm does not exhibit this feature.
7

---

Importance sampling
1010
30"
108
20oles a1rano0
00'4
7071
10-11
10:11
10147 10-2
32
20
12F
107
time behween aiservaations
between ossernation
(a)
(b)
Figure 3: Gradient variances function of the time between observations. shows gradient variances of the observation weea-squared-error computed with an SDE integrator. shows the gradient variance of the log-probability of observations computed with Algorithm 1.
In each case used the same neural network, initialized with fixed seed, and the exact same values for all parameters. We used the same loss functions as Section 3.1. The data set we used was the Hungarian chickenpox cases dataset from the UCI machine learning repository:" it has 20 features per observation. We scaled the features by their standard deviations before performing the experiment improve numerical stability.
Our findings are as follows:
For very short times, found that both algorithms produced comparable gradient variances: typically less than 10-7
For longer time spans (3 orders of magnitude larger than the discretization step size), the gradient variance of our algorithm increased by few orders of magnitude; the maximum we observed was on the order of 10-4, However, the SDE integrator approach resulted in huge increases in variances, with some as large as 10
Between the above two extremes, we found that the gradient variance of our approach did increase with the time between observations, but stayed below 10-5, while the gradient variance from the SDE integrator exhibited very poor scaling with variances of up to 10th for times between observations of 1.0.
Considering that both approaches produced comparable variances on very small timescales, we feel that this experiment represents fair comparison.
3.3 UNCERTAINTY QUANTIFICATION WITH VARIATIONAL GAUSSIAN PROCESSES
We attempt to test the viability of using our method for an SDE model with GP prior over the drift function and fitting the posterior variationally. To assess this, we used data trajectories from van der Pol oscillator. For the GP prior, we used zero mean function and radial basis function kernel, with independent priors over each output dimension the drift vector field. These choices are not necessarily optimal; however, the purpose of this experiment is merely to illustrate the feasibility of using variational GP to learn the vector field an SDE model.
We used grid-based inducing point approach and the kernel interpolation method of Wilson & Nickisch (2015). The combination of these techniques allows one to take advantage the Kronecker product structure of the covariance matrices and the Toeplitz structure of each factor for accelerated sampling. The model was trained variationally, sampling the variational posterior jointly along entire trajectories.
*This dataset can be found at https::/arde ve. ics.uct Leduee//nnn /ml /ttp///ataassaaaeeeoeeoaa Hunga ant Chi Cases
8

---

not wholly Deviation predictive ground truth
0.25
0.70
0.35
14.
0.30
0.00
standard dhetelitt from menn
(a)
(b)
(c)
Figure 4: GP posterior compared to ground truth vector field. shows the GP vector field coloured by uncertainty. The arrows show the integral curves of the posterior mean and the coloured are scaled linearly with the log predictive variance. shows the integral curves of the ground truth van der Pol oscillator. shows the deviation of GP posterior from the ground truth vector field in standard deviations. The red lines show the 10%, 25%, 75%, and 90% quantiles.
The approximate GP posterior is shown in Figure compared the ground truth vector field the van der Pol oscillator in Figure 4b. Comparing the GP to the van der Pol vector field, we found that the GP model was fairly accurate in areas with many data points, while reverting to the prior far away from observations. To assess the quality of the uncertainty quantification, we further calculated the discrepancy between the variational posterior and the ground truth vector field. Over a dense grid, covering an enlarged area around the observations, calculated the discrepancy as
Svdp (x) - p(x)
5
(17)
(σx)
histogram is shown in Figure 4c.
Under Gaussianity assumptions, one expects to see the quantiles line up well with some associated number of standard deviations. However, this will only be the case near to the observations. In regions where the GP reverts to the prior, which has stationary kernel, the mismatch between the GP and the polynomial growth of the van der Pol oscillator results in heavy tails. We expect that better choice of kernel, which accounts for this polynomial behaviour, to show residuals with lighter tails. We found that the distribution of these ds was fatter-tailed than a Gaussian distribution; the middle 50% (25% 75% quantiles) slightly passed one standard deviation either side, opposed to the 68% a standard Gaussian distribution. The middle 80% (10% 90% quantiles) were within 3 standard deviations. While the posterior uncertainty was somewhat overconfident, unreasonable.
We further note that, while the complexities involved in variational GP fitting can result in longer training times, we still found that our method applied to a GP was significantly faster than using an SDE integrator with neural network.
4 DISCUSSION
In this work we introduced an algorithm to efficiently estimate probabilities of observations under the assumption of SDE model for the purposes of learning. The proposed approach produces estimates which are accurate, even in high dimensions, and can fully utilize parallel hardware for significant speedups.
The proposed algorithm is capable of providing robust probability estimates which can be used in models both with or without an explicit observation noise component. This is in contrast to SDE integrator- based learning approaches, which must incorporate and estimate explicit error term due to their inability to generate paths which are guaranteed to pass through the observations.
Our algorithm produces more stable estimates than alternative integrator-based approaches, particu- larly when the observations are sparsely distributed over long time intervals and in high dimensions. This in turn results in more stable and robust learning. With the rise in applicability of large- -caee (neural network) SDE modeling in science, simulation, and domains such as finance, we see our approach an important contribution reliable scaling of these methods.

---

ACKNOWLEDGMENTS
The authors would like to thank InstaDeep Ltd. for their financial support of this work.
REPRODUCIBILITY
Efforts have been made to describe all parameters and details of the algorithms used in this paper, either in the content or in the appendix. The code included in the supplementary material includes Dockerfile and instructions for running some of the experiments. Where relevant, fixed seeds for random number generators can be found in the code during initialization.
REFERENCES
Philipp Batz, Andreas Ruttor, and Manfred Opper. Approximate Bayes learning of stochastic differential equations. Phys. Rev. E, 98:022109, Aug 2018. doi: 10. 10.1103/PhysRevE.90 wtt.seE.....02209.
URL https:// http:s://liik... aps owp..oog///////////////// 1100 1103/PhysRevee 98 022109.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary
differential equations. 31, 2018. URL https::/proceedding wtt.:rocceednnnnn..eri.apaaa neurips cc/ sapppppr/
2018/ /ttee/lee//////eeeeee 69186699666666 wwwpldffedeee2244eeeeoeeee
Sarat C. Dass, Jaeyong Lee, Kyoungjae Lee, and Jonghun Park. Laplace based approximate posterior
inference for differential equation models, differential equations. Statistics and Computing. 27: 679-698, 2017. doi: 10.
V. Girsanov. On transforming a certain class of stochastic processes by absolutely continuous substitution of measures. Theory of Probability & Its Applications, 5(3):285-301, 1960. doi: 10.1137/110500 URL https: /tt./or..gg///////////o//oo 10.100.1111 http///..111/111111
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian Neural Networks. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Gar- nett (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Asso- ciates, Inc., 2019. URL https: http:::://ooocpcpeoeeeeeeooe (twp://cceeeeeeeeeeeeeeee888833
wwweeriss.ccppppeeeeeeeeeeepeeee pctppppe///19///////c/p 2019//119/9///////////// Paper.pds
Martin Hutzenthaler and Arnulf Jentzen. On perturbation theory and on strong convergence rates for stochastic ordinary and partial differential equations with nonglobally monotone coefficients. The Annals of Probability, 48(1), Jan 2020. ISSN 0091-1798. doi: 10.12114/9-aop1345. URL
http:/ /ttp//xxx.a doi. org/ htt.g/..1114..... 19-AAP1345.
Junteng Jia and Austin R Benson. Neural jump stochastic differential equations. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Gamett (ods.), Advances in Neural Information Processing Systems, volume 32. Curran Asso- ciates, Inc., 2019. URL https::/proceedeeee
neurips cc/ paper/2019/Eile. pdf.
Hadiseh Karimi and Kimberley B. McAuley. Bayesian Estimation Stochastic Differential Equation Models via Laplace Approximation. IFAC-PapersOnLine, 49(7):119-11144 2016. ISSN 2405-8963. doi: https://doi.org/ org/10. O. (0166...faco..206607.3331.1 URL https: //www scienced i rect. com/ sci ence/ artii irticti / 933055888666669 6305585444 11th IFAC Symposium on Dynamics and Control of Process SystemsIncluding Biosystems DYCOPS-CAB 2016.
Patrick Kidger, James Foster, Xuechen Li, and Terry J Lyons. Neural SDEs as Diiieni--immnniionall GANs. In Marina Meila and Tong Zhang (eds.), Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 5453-5463, URL https: http://pcccceeeeeeeeeeee mlr wwt.:ppwwsssseeeeeeeeeess
PMLR, 18-24 Jul 2021.
html.
Prashant Kumar, Martin Schmelzer, and Richard P. Dwight. Stochastic turbulence modeling in RANS simulations via multilevel Monte Carlo. Computers & Fluids, 201:104420, 2020. ISSN 0045-7930. doi: com/ http://aaaaaneerrtrree//////////
URL https: //www http://wwwww.nc..c sciencedireett.
10

---

Xuechen Li, Ting-Kam Leonard Wong. Ricky T. Q. Chen, and David Duvenaud. Scalable gradients for stochastic differential equations. International Conference on Artificial Intelligence and Statistics, 2020.
Xuanqing Liu, Tesi Xiao, Si Si, Qin Cao, S. Kumar, and Cho-Jui Hsieh. Neural SDE: Stabilizing Neural ODE Networks with Stochastic Noise. ArXiv, abs/1906.02355. 2019.
Paul Malliavin and Anton Thalmaier. Stochastic Calculus of Variations in Mathematical Finance. Springer, 01 2006. ISBN 3-540-43431-3. doi: 10.1007/3-540-307777
James Morrill, Cristopher Salvi, Patrick Kidger, and James Foster. Neural rough differential equations for long time series. In Marina Meila and Tong Zhang (eds.). Proceedings of the 38th International Conference Machine Learning. volume 139 of Proceedings of Machine Learning Research, pp. 7829-7888 7838. PMLR, 18-24 Jul 2021. URL https: httpsss//ppopeeeedeo
www.mprresse..e...mm v139/
morrill211 html.
GA Pavliotis. Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations. Springer, 2014.
Michael Poli, Stefano Massaroli, Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Graph neural ordinary differential equations. 2021.
Alvaro Lanchez-Gonzalez. Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian Graph Networks with ODE Integrators, 2019.
Belinda Tzen and M. Raginsky, Neural stochastic differential equations: Deep latent Gaussian models 2019.
in the diffusion limit. ArXiv,
Andrew Wilson and Hannes Nickisch. Kernel Interpolation for Scalable Structured Gaussian Pro- cesses (KISS-GP). In Francis Bach and David Blei (eds.), Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp. 1775-1784, Lille, France, 07-09 Jul 2015. PMLR. URL pxess/v37/wi1sst
mlr.
html.
C. Yildiz, M. Heinonen, Intosalmi, H. Mannerstrom, and H. Lahdesmaki. Learning stochastic differential equations with gaussian processes without gradient matching. In IEEE International Workshop on Machine Learning for Signal Processing. MLSP, 2018.
Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Symplectic ode-net: Learning hamiltonian dynamics with control. In International Conference on Learning Representations, 2020. URL
net http//oormmmmddmmo rKDS.
A THEOREM 1
We have that X and Y are processes obeying the SDEs
dX f(X) di (X)) dY = ggYY) dw.
(18)
(19)
dropping the time argument for notational simplicity. Expectations with respect to SDEs can be expressed as path integrals. In the Itô convention, the path integral representation of an SDE has the following form
dr n
N
Ax, f (I
2
Ex [[[]]
lim N
Maa] II
exp
Σ
2 At Atng(Tn-))
2
g(En-1) 1) At,
(20)
= DDFlu))J exp {
{-1 ((a(t),(t))
1-2
(21)
11

---

where J is Jacobian factor which depends on 9 but not f and the Lagrangian is given by
④ [(x)
C(x, .)
(22)
g(x)
The derivative of x: here may be thought of distributionally; however, the expression for the measure Dr - f -}} need only be defined the Cameron-Martin subspace", guarantecing its unique extension to the whole Wiener space. In contrast to this, the expression inside the expectation F[æ] should be defined almost-everywhere, and hence care should be taken when à appears in an expectation.
We may expand the quadratic separate the Lagrangian into drift-free term plus terms dependent on the drift
C
(23)
g(z)
2f (x)g (x)}
(24)
+ (
(25)
The first term, on line 23, the Lagrangian for the process Y. which has drift term. Shuffling the rest of the terms out of the measure and into the expectation, we have
Ex[ = / DrF exp { (
2
dt al..
(26)
g(x))
where
s(x] / - (X,)
(27)
1 1/2 ((xi) (xx)) dt
(28)
2
/ 1 / (xx))
(29)
(30)
Note that on the third line we have used the distributional definition of X; i.e. X, dt = dX,. Lastly, using the definition of conditional expectation, have that
EA [[A]pp(B = bb/a] = EA (p b| B
EA [F[A] F(A] B B= b)
(31)
By trivial substitution, we have
(x(xIN) = py (ELN) [M (t = J
(32)
B SAMPLING BRIDGES
Let Y be Gauss Gauss-Markov Markov process with kernel k. Since the process is Markov, need only sample
the process conditioned on two end- points. Let y = (Ho-yr) be the known values the process at times 7 = (fo, T). Let Y be the process defined by
y, = Y. + ((,))
for ISIST.
(33)
where K is the matrix whose entries are k(t;, (z). This distribution of Ý is the same the distribution of conditioned on the values at times т. Note that once Y, has been sampled - which may
'The Hilbert space functions whose derivatives are square integrable. The Wiener space is the space of continuous functions. This is not Hilbert space, instead it is a Banach space with the supremum norm. The measure of the Cameron-Martin space is in fact zero; the sample paths of the SDE are almost-surely nowhere differentiable.
12

---

require cumulative sum this transformation be performed independently, and hence in parallel, for each For each consecutive pair of observations, the bridge process can be sampled independently.
In the simplest case, where Y a Brownian motion starting from yo at time to = 0, this reduces to
Y Y
kIT.) M(TT) (yT Yr)
(34)
Yi (vr Yr).
(35)
And similarly for any other processes which can be sampled with a fixed starting point.
In the more general case we can invert the two by two matrix to get an explicit formula for Y,
y, Y, tt),t) T) k(t, (to. T) (yo - Yo) k(lo, (T,
(36)
k(lo. 7)
(tt, to) - k(t, to (0) (to, 1) k(lo, 40 (T. (T, k(4,, T)2 (ur YT). integral, or derivative of the process etc.
(37)
This may be extended straightforwardly to other cases, such as conditioning on the values of the
C EXPERIMENTS
C.I PROBABILITY ESTIMATES
As sanity check, compare the pog-bbbbbbylity estimates of our algorithm to the actual log- probabilities of an Ornstein-Uhlenhed process. The probability density for this process is available in closed form for any number of observations. We estimate the probability of set of data points evenly spaced over a time interval of [0, 10), taken from sine function In = sin("/10 x t.). For these estimates we used 1000 time steps to estimate the integral S, and averaged over 100 trajectories, independently sampled from a Brownian bridge passing through the observations. We repeated this process 1000 times and report the mean and standard deviation of these estimates. Figure shows these estimates compared to the analytic probability given by the Onsseenn-Uhlenbed process for various numbers of observations as a function of the dimension of the problem.
probability estimates ground truth
-20
bag
-000
-000
Figure 5: Log-probability estimates for Omstein- process. The solid lines indicate the mean value of the estimator, and the shaded regions indicate three standard deviations from the mean. The dashed lines are the analytic probbabbiiyl values.
An interesting observation is that the variance of the estimator seems to decrease with the number of data points. This atypical of probability estimators, since the magnitude of the log-prood probbabiiyy increases linearly with the number of observations (not to be confused with the number of samples
used to estimate the log probability). In almost all cases, with the exception observations in more than 40 dimensions, algorithm gives accurate log-probability estimates with small standard
13

---

deviations. The error of this exception is approximately 30%, but is significantly reduced by adding more observations, or using more samples. We believe this error is acceptable considering the small number of observations and the large scale. For 10 observations in 100 dimensions the standard deviation of our estimator is only 0.4% of the magnitude of the log-probability and the error of the mean estimate is 2% of the analytic hog-probability. This small bias can be further reduced by using more samples or smaller At. A single estimate in the 100 dimensional case took less than half second compute on laptop CPU, and therefore feasible use in the inner loop of learning algorithm.
C.2 LORENZ
C.2.1 DATA SET
We generated data by simulating a Lorenz system with the exact same parameters as Li et al. (2020). We generated 16 independent trajectories, with initial values sampled from standard normal distribution, and discretely sampled the values the process over 200 equally spaced time points ranging from 0 to 10. Unlike Li al. (2020), we did not add any observation noise the data.
We deliberately used a sparser data set than Li et al. (2020), due to personal experience of the difficulty of fitting SDE models on sparse observations.
C.2.2 MODEL
The model we used to fit to the observations was
dX, fo())) dt + a dW,
(38)
where is scalar and fo is fully-connected neural network with ReLU nonlinearities and parame- ters 0. We used layer sizes of 3 + 32 + 256 + 32 + 3.
C.2.3 TRAINING PARAMETERS
We used precisely the same parameters for each model, the only difference in the training procedures being the loss functions used. In the first case, the loss function per trajectory was
Liuse (0; 21:N) KN
Σ( _In)² k.n
(39)
where = 200 is the number observations per trajectory, and K is the number of sampled paths used to estimate the MSE, and (*) is sampled using an SDE integrator starting at ((k) then averaged this loss function over the 16 independent observed trajectories.
t1 21 We
In the second case, the loss function per trajectory was
1
Lpi xI.N)
log ( [0),
(40)
N
where Px (#1.N 0) is calculated using Algorithm 1 with K trajectory samples. Again this loss was averaged over the 16 independent observed trajectories. In both cases we used K = 64 and a time step size of A/ 10-2
For each model we used the Adam optimizer with a learning rate of 10 and ran the optimization algorithm for 10 iterations. The model using our importance sampling approach completed the full number of iterations in 34 minutes, while the integrator-base approach took 29 hours and 46 minutes to complete the same number iterations.
Our approach consistently maintained 100% GPU utilization, except during MSE evaluation which used an integrator, while the integrator based approach was unable to surpass 25% GPU utilization in this experiment.
C.3 VARIATIONAL. GAUSSIAN PROCESSES
Our approach allows more straightforward implementation SDE models using a Gaussian process prior for the drift function. Some previous approaches often resort approximations such as
14

---

using only the predictive mean MAP estimates (Yildiz al., 2018). For dense observations, directly discretizing the SDE and fitting the GP to the observation deltas may be viable. For sparse observations this is much more difficult. Our proposed method makes variational GP methods slightly simpler, since the values of f (() can be sampled jointly along the path. However, it is important to note that various complexities still arise such the difficulty of efficiently sampling high-dimensional Gaussian distributions.
We used the kernel interpolation method of Wilson & Nickisch (2015) which stores a grid of inducing locations per dimension. Storing inducing locations per dimension for a Kronecker product kernel effectively gives a variational distribution over m inducing points with only requirements, where m is the number of inducing points per dimension and d is the dimension. When these inducing locations are evenly spaces, the covariance factors have a Toeplitz structure which allows sampling in linearithmic time complexity.
memory
To fit the model we used gradient ascent on the ELBO which is calculated as follows:
ELBO Ey-g [ogg((>[[)] KT [(()))())())
(41)
The KL divergence is between two Gaussian distributions calculated at the inducing locations in the usual way. For the likelihood term, we generated trajectories (as in Algorithm 1), followed by jointly sampling the variational GP posterior along these trajectories and proceed with the rest of the computation Algorithm using these samples. Our optimization objective is calculated
L
1 Σ K
log login ((k)) KI. ((T()))))))) K
(42)
where each (k) is an independent sample over the whole path of Y. In other words, ( (Y,,) and ((() (Y²,) are correlated, but (Y,, and (Y,,) are not when k + j.
In our experiment we used 32 evenly spaced inducing points per dimension, and jointly sampled the variational posterior over 64 trajectories of 1024 time steps. The model was trained with the Adam optimizer with initial learning rate of 0. and an exponential decay with factor 0.999 for 10a training steps. However, the square deviation of the predictive mean from the ground truth converged well before 10 iterations. Training for 10 iterations took 34 minutes, 10 iterations took about and half hours.
15


Reviews:
[{"review_id": "SVwyiv3Ota9", "reviewer": "Reviewer_p6RL", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "- **Clarity.** The paper is mostly written clearly and succinctly. However, it could be more self-contained, possibly with additional background provided in the supplementary materials.\n- **Reproducibility.** Low-level details necessary to reproduce the results from this paper are not all included. However, source code is provided with the supplementary material, which appears to be well-structured and clean.\n- **Novelty.** The main technical novelty comes from the application of Theorem 1, as summarized in Algorithm 1. In particular, the authors propose to estimate the original SDE of interest with an linear SDE, which is an attractive alternative since it is easier to estimate their probabilities and expectations with importance sampling, and has the added benefit that the sampled paths pass through the observations. This appears to provide significant improvement, as evidenced by the empirical results.\n\n## Concerns\n\nOverall, I find the empirical evaluations feel a bit rushed.\n\n- Experiment described in Section 3.1: are you making use of parallelism in this experiment? If so, the details of this need to be described. If not, can you pinpoint and provide some explanation for the factor of 52 speed-up? From what I can tell, even without any parallelism, the method is able to provide speed-ups of this magnitude? Furthermore, the plots might be clearer visualized in log-y scale. Also, instead of truncating the training time along the horizontal axis, you might consider instead recording the amount of time required to attain a particular MSE value. As it is currently, the differing lengths of the curves make it slightly messy. Finally, the only baseline that is compared against is the adjoint method of Li et al. 2020. I would be interested in comparing the proposed method against some of the related methods highlighted in the Related Works section as baselines.\n- Experiment described in Section 3.3: I find the normality test carried out starting on line 271 to be lacking in precision. Firstly it is entirely qualitative and therefore somewhat subjective. There are simple statistical tools available to carry out normality tests. Line 280 concerning training speeds: the speed-up is mentioned in passing. Though not surprising, it would be good to see this quantified more rigorously with a table or figure, rather than as an offhand remark.\n\n## Miscellaneous Issues\n\n- The statement of contribution at the end of page 2 contains four bullet points. Only the first bullet point is the actual contribution. The remaining three are merely descriptions of the first and only contribution.\n- Line 14: \"which\" → \"that\"\n- Line 46: \"Gaussian process, [...]\" - extraneous comma\n- Line 146: \"which\" → \"that\"\n- Line 183: \"which\" → \"that\"\n- Line 222: \"The second metric we record the [...]\" - missing \"is\"\n- Figure 2: If I understand correctly, the two sub-figures shown here are of the same vector field and observations. Is that right? This would be made much more clear by having the subplots share the same vertical and horizontal axes.", "summary_of_the_paper": "This paper proposes an importance-sampling method for estimating the likelihood of observations in stochastic differential equation (SDE) models. Unlike classical integrator-based methods, the proposed approach imposes a structure that forces sampled trajectories to pass through observations thereby improving accuracy, especially in high dimensions. Furthermore, it removes sequential dependencies between sampled points thus allowing for speed-ups through parallelism.", "summary_of_the_review": "Given the numerous advantages of the proposed method over existing works, I am inclined to recommend acceptance of this paper. However, given the questions I raised concerning the empirical evaluations, I am reluctant to recommend a clear acceptance.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "7Uf2qI1LJjl", "reviewer": "Reviewer_k19Y", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "+ The paper indeed addresses a timely topic. Given that continuous-time dynamical system literature is currently dominated by neural ordinary differentia equations, contributions like this would pave the way for more stochastic differential equation publications.\n+ The results seem pretty impressive. However, the experiments are done only on very few dimensional systems. More experiments on bigger systems would help us identify trade-offs and potential issues with the methodology.\n+ The paper is written clearly and the very obvious motivation is communicated very well. A few notes to further improve the readability:\n    - What does the sentence in line 48 mean? More specifically, what prior/posterior distributions do you refer to?\n    - More details on the GP-SDE model in Sec 3.3 would be nice. Also, the experiment seems a bit repetitive.\n    - Alg 1 should be written more verbosely, e.g., the set of $t$'s (in line 3) and $\\Delta t$ can be explicitly given.\n\nThese being said, I'm not able to verify Th 1. Despite checking the references and more resources, I could not immediately see (9-10). If the authors can provide a proof (even if it is straightforward) or other reviewers can comment on it, I would be very happy (including the sum that appears in the last line of Alg 1). Similarly, (13) and the paragraph above is not clear. As I understand, an invertible mapping $T$ is utilized to switch between processes with constant and state dependent diffusions. Yet, I'm not able to immediately see (13). It would be much better to clearly define all the terms/processes and show where the (state) independence comes from.", "summary_of_the_paper": "The paper presents a new approach to compute the density of a given observation sequence under arbitrary drift and diffusion functions. The method relies on Girsanov's theorem, which requires evaluating the density of the observations under Wiener process. The authors further propose an inversion trick to handle state dependent diffusion functions. The experiments show that the method indeed learns accurate drift and diffusion functions much faster than standard integration-based approaches and has much smaller gradient variance", "summary_of_the_review": "The proposed approach concerns an interesting topic and provides excellent improvements over existing approaches. However, I'm not able to mathematically verify Theorem 1 and the derivation in Section 2.2. This is why I preliminarily vote for a reject but would be very happy to re-consider my score if additional details are provided by the authors and discussions with other reviewers.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "-ZdDDQBLqQ", "reviewer": "Reviewer_68Rx", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "While the proposed method improves the speed and efficiency of SDE learning, the presentation is not strong enough. Below are my specific comments.\n\n* The scope of the work is limited to state-independent or transformable to state-independent diffusion. This is not properly reflected in the title and abstract of the paper and should be inferred by the reader throughout the paper.\n* Methodological novelty is limited. The paper is heavily based on a simple application of theorem 1.\n* Comparisons are done only on the toy examples, if the authors believe that the proposed method suffices for most application then results on higher dimensions, various simulations, and real world data should be presented. Provided that the technical novelty is limited, empirical validation is necessary for the publication.\n* Comparisons against the adjoint method is not completely fair, since it solves a way more general problem (state-dependent diffusion). Are there methods tailored to state-independent diffusion to be used for comparison?\n* Why gradient variance results is shown on a separate dataset than learning curves? Please include more than one example for the learning curves experiment and include standard deviation across multiple experiments. Also, can you show this as a function of dimensions, and when adding or removing observation noise?\n\n\n\n\nHere are some minor comments:\n\n* Introduction is very short and many applications of SDE are not mentioned.\n* Line 25: \"In many physical applications ...\" please provide examples and include citations.\n* Line 27: \"In more complex processes ...\"  please provide examples and include citations.\n* Line 46: \"Models that implicitly describe f as a Gaussian process\"  please provide examples and include citations.\n* Line 56: \"Can marginalize out the observation noise ...\" please provide examples of this, how do you marginalize out the observation noise?\n* Theorem 1: please spell out the conditions of Girsanov's theorem explicitly and what exactly should hold for function f and sigma.\n* Instances of informal language can be found throughout the text such as \"wandering around?\" or \"if need be\".\n* Line 184: \"It may be argued ...\" please provide evidence for this. Can you support the argument \"almost all real world problems ...\"  by providing examples from the literature?\n* Some of the text in the supplementary must be brought to the main. In the current form, it's not clear how exactly the experiments are done. I had a grasp on the experiments section after reading the supplementary.\n* Variational GP experiment is very unclear, what is the generative model? What is data (include equations)? What is the inference procedure?", "summary_of_the_paper": "The paper proposes an importance sampling method for the estimation of probabilities in SDE models, applicable to SDEs with state-independent (or transformable to state-independent) diffusion. By taking the advantage of Girsanov's theorem, the probability density of the observations is rewritten as a weighted probability with respect to a base Brownian motion process. Sampling results are presented qualitatively and compared to integration based sampling. Results also suggest that the proposed method is faster than integration based methods and the gradients have lower variances. In addition, a variational GP examples is demonstrated showing that the method can provide reasonable uncertainty estimation.", "summary_of_the_review": "In summary, the methods presented can be of potential interest to the community, however the presentation is weak, scope and novelty is limited, and experiments are rather minimal.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "1: The contributions are neither significant nor novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes an importance-sampling estimator for probabilities of observations of SDEs. The proposed approach has several advantages over conventional methods: it does not require an SDE solver, it has lower gradient variance, and shows nice results with a Gaussian process representation of the function. Reviewers were somewhat split on this paper, with some concerns that experiments were limited. On balance, however, the paper makes several nice contributions, the experiments are in line with related works, and the authors did a good job of clarifying Theorem 1 in the rebuttal. We note that Reviewer K19Y changed their opinion to accept (although they forgot to update the score). Please carefully account for all reviewer comments in the final version.
Decision: Accept (Poster)