Learning Stochastic Dynamics from Data

Authors: ziheng  guo
igor  cialenco
ming  zhong
Venue: ICLR_cc_2024_Workshop_AI4DiffEqtnsInSci
Published: 2024-02-11 06:43:55

Abstract: We present a noise guided trajectory based system identification method for inferring the dynamical structure from observation generated by stochastic differential equations.  Our method can handle various kinds of noise, including the case when the components of the noise are correlated.  Our method can also learn both the noise level and drift term together from trajectory.  We present various numerical tests for showcasing the superior performance of our learning algorithm.

Main paper: 
 introduction: 
 Stochastic Differential Equation (SDE) is a fundamental modeling tool in various science and engineering fields Evans (2013) ; Särkkä" & Solin (2019) . Compared with traditional deterministic models which often fall short in capturing the stochasticity in the nature, the significance of SDE lies in the ability to model complex systems influenced by random perturbations. Hence, SDE provides insights into the behavior of such systems under uncertainty. By incorporating a random component, typically through a Brownian motion, SDE provides a more realistic and flexible framework for simulating and predicting the behavior of these complex and dynamic systems.We consider the models following the SDE of the following form, dx t = f (x t )dt + dw t , where the state vector x t ∈ R d , the drift term f : R d → R d , and the stochastic noise w t is a Brownian motion with covariance matrix D(x) which also depends on the state. In the field of mathematical finance, there are several SDE models that are widely used, for example, Black-Scholes Model Black & Scholes (1973) in option pricing, Vasicek Model Vasicek (1977) for analyzing interest rate dynamic and Heston Model Heston (2015) for volatility studies. Apart from the field of finance, SDE has its application in physics Coffey et al. (2004) ; Ebeling et al. (2008) and biology Székely & Burrage (2014) ; Dingli & Pacheco (2011) where SDE is used to study evolution of particles subjected to random forces and modeling physical and biological systems.The calibration of these models is crucial for their effective application in areas like derivative pricing and asset allocation in finance, as well as in the analysis of particle systems in quantum mechanics and the study of biological system behaviors in biology. This requires the use of diverse statistical and mathematical techniques to ensure the models' outputs align with empirical data, thereby enhancing their predictive and explanatory power. Since these SDE models have explicit function form of both drift and diffusion terms, the calibration or estimation of parameters is mostly done by minimizing the least square error between observation and model prediction Mrázek & Pospíšil (2017) ; Abu-Mostafa (2001) . SDE has also been studied for more general cases where the drift term does not have an explicit form where it satisfies f = f (x t , θ) with θ being a vector a unknown parameters. These models can be estimated with the maximum likelihood estimator approach presented in Phillips (1972) . The other approach is to obtain the likelihood function by Radon Nikodym derivative over the whole observation {x t } T t=0 . Inspired by theorem 7.4 in Särkkä" & Solin (2019) and the Girsanov theorem, we derive a likelihood function that incorporates state-dependent correlated noise. By utilizing such likelihood function, we are able to capture the essential structure of f from data with complex noise structure.The remainder of the paper is structured as follows, section 2 outlines the framework we use to learn the drift term and the noise level (if not state dependent), we demonstrates the effectiveness of our learning by testing it on various cases summarized in section 4 and additional examples presented in Appendix B, we conclude our paper in section 5 with a few pointers for ongoing and future developments. 
learning framework: 
 We consider the following SDEEQUATIONwhere f : R d → R d is a drift term, and w represents the Brownian noise with a symmetric positive definite covariance matrix D = D(x) : R d → R d×d .We consider the scenario when we are given continuous observation data in the form of {x t , dx t } t∈[0,T ] for x 0 ∼ µ 0 . We will find the minimizer to the following loss functionE H ( f ) = E x0∈µ0 1 2T T t=0 < f (x t ), D -1 (x t ) f (x t ) > dt -2 < f (x t ), D -1 (x t )dx t > ,(2) for f ∈ H; the function space H is designed to be convex and compact and it is also data-driven. This loss functional is derived from the Girsanov theorem as well as inspiration from Theorem 7.4 from Särkkä" & Solin (2019). In the case of uncorrelated noise, i.e. D(x) = σ 2 (x)I, where I is the d × d identity matrix and σ : R d → R + is a scalar function depends on the state, representing the noise level, 2 can be simplified toEQUATIONMoreover, we provide three different performance measures of our estimated drift. First, if we have access to original drift function f , then we will use the following error to compute the difference between f (our estimator) to f with the following normEQUATIONHere the weighted measure ρ is defined on Ω, where it defines the region of x explored due to the dynamics defined by equation 1; therefore ρ is given as followsρ(x) = E x0∼µ0 1 T T t=0 δ xt (x), where x t evolves from x 0 by equation 1.(5)The norm given by equation 4 is useful only from the theoretical perspective, i.e. showing convergence. In real life situation, f is most likely non-accessible. Hence we will look at a performance measure that compare the difference between X(f , x 0 , T ) = {x t } t∈[0,T ] (the observed trajectory evolves from x 0 ∼ µ 0 with the unknown f ) and X( f , x 0 , T ) = {x t } t∈[0,T ] (the estimated trajectory evolves from the same x 0 with the learned f and the same random noise as used by the original dynamics). Then, the difference between the two trajectories is measured as followsEQUATIONHowever, comparing two sets of trajectories (even with the same initial condition) on the same random noise is not realistic. We compare the distribution of the trajectories over different initial conditions and all possible noise at some chosen time snapshots using the Wasserstein distance. 
comparison to others: 
 When the noise level becomes a constant, i.e. σ(x) = σ > 0, we end up a much simpler lossE Simpler H ( f ) = E x0∈µ0 1 2T σ 2 T t=0 || f (x t )|| 2 dt -2 < f (x t ), dx t > ,which has been investigated in Lu et al. (2022) . System identification of the drift term has been studied in many different scenarios, e.g. identification by enforcing sparsity such as SINDy Brunton et al. (2016) , neural network based methods such as NeuralODE Chen et al. ( 2018) and PINN Raissi et al. (2019) , regression based Cucker & Smale (2002) , and high-dimensional reduction framework Lu et al. (2019) . The uniqueness of our method is that we incorporate the covariance matrix into the learning and hence improving the estimation especially when the noise is correlated. 
examples: 
 In this section, we demonstrate the application of our trajectory-based method for estimating drift functions, showcasing a variety of examples. We explore drift functions ranging from polynomials in one and two dimensions to combinations of polynomials and trigonometric functions, as well as a deep learning approach in one dimension. The observations, serving as the input dataset for testing our method, are generated by the Euler-Maruyama scheme, utilizing the drift functions as we just mentioned. The basis space H is constructed employing either B-spline or piecewise polynomial methods for maximum degree p-max equals 2. For higher order dimensions where d ≥ 2, each basis function is derived through a tensor grid product, utilizing one-dimensional basis defined by knots that segment the domain in each dimension.The parameters for the following examples are listed in table 1. Due to space constraints, additional examples are provided in appendix B. The estimation results are evaluated using several different metrics. We record the noise terms, dw t , from the trajectory generation process and compare the trajectories produced by the estimated drift functions, f , under identical noise conditions. We examine trajectory-wise errors using equation 6 with relative trajectory error and plot both f and f to calculate the relative L 2 error using 4, where ρ is derived by 5. Additionally, we assess the distribution-wise discrepancies between observed and estimated results, computing the Wasserstein distance at various time steps. We initiate our numerical study with a one-dimensional (d = 1) drift function that incorporates both polynomial and trigonometric components, given by f = 2 + 0.08x -0.05 sin(x) + 0.02 cos 2 (x).Figure 1 illustrates the comparison between the true drift function f and the estimated drift function f , alongside a comparison of trajectories. Notably, Figure 1 (a) on the left includes a background region depicting the histogram of x t , which represents the distribution of observations over the domain of x. This visualization reveals that in regions where x has a higher density of observations-indicated by higher histogram values-the estimation of f tends to be more accurate.Conversely, in less dense regions of the dataset (two ends of the domain), the estimation accuracy of f diminishes. Table 2 presents a detailed quantitative analysis of the estimation results, including the L 2 norm difference between f and f , as well as the trajectory error. Furthermore, the table compares the distributional distances between x t and xt at selected time steps, with the Wasserstein distance results included.We continue our numerical investigation with a one-dimensional (d = 1) drift function which is given by f = 0.08x. especially towards the two end points of the interval. However, the errors happen mostly during the two end points of the data interval, where the distribution of the data appears to be small, i.e. few data present in the learning. We are able to recover most the trajectory.In this example, we incorporate a non-diagonal covariance matrix into a 2 dimensional SDE system. As specified in table 1, all parameters remain unchanged except for D and M . We change total trajectory observation number M to 1000 for faster calculation andD = 0.6 0.2 0.2 0.8This change in D implies that the Brownian motions within the system are correlated. The drift function is defined using the notation,f = [f 1 (x) f 2 (x)] ⊤ and x = [x 1 x 2 ]⊤ , where f i : R 2 → R and x i ∈ R for i ∈ {1, 2}. For polynomial drift function f , we set f 1 = 0.4x 1 -0.1x 1 x 2 and f 2 = -0.8x 2 + 0.2x 2 1 . The estimation result can be found in figure 3a and 3b and table 3.  
conclusion: 
 We have introduced a novel methodology for learning the drift and diffusion components through a trajectory-based loss function guided by noise. This loss function is derived from the negativelog ratio of likelihood functions, quantifying the ratio of probabilities of observing two stochastic processes sharing the same initial starting point. Our approach accommodates various noise structures, provided that the covariance information is available. Even in cases where the noise structure is unknown and the covariance is a scalar, our method remains effective through the utilization of quadratic variation. However, ongoing research aims to explore scenarios where the noise exhibits state dependence and matrix structure. Evaluation on high-dimensional drift terms and complex noise structures is currently underway. Our loss function solely relies on observable data (trajectories) and accounts for potential correlated noise. 
 
 Appendix: In this section, we will discuss in details how the algorithm is implemented for our learning framework presented in section 2. In real applications, continuous data is rarely obtainable, we will have to deal with discrete observation data, i.e. {x m l } L,M l,m=1 with x m 0 being i.i.d sample from µ 0 , whereEQUATIONfor f ∈ H. Moreover, we also assume that H is a finite-dimensional function space, i.e. dim(H) = n < ∞. Then for any f ∈ H, f (x) = n i=1 a i ψ i (x), where a i ∈ R d is a constant vector coefficient and ψ i : Ω → R is a basis of H and the domain Ω is constructed by finding out the min / max of the components of x t ∈ R d for t ∈ [0, T ]. We consider two scenarios for constructing ψ i , one is to use pre-determined basis such as piecewise polynomials, Clamped B-spline, Fourier basis, or a mixture of all of the aforementioned ones; the other is to use neural networks, where the basis functions are also trained from data. Next, we can put the basis representation of f back to equation 7, we obtain the following loss based on the coefficientsEQUATIONIn the case of covariance matrix D being a diagonal matrix, i.e.D(x) =     σ 2 1 (x) 0 • • • 0 0 σ 2 2 (x) • • • 0 . . . . . . . . . . . . 0 0 • • • σ 2 d (x)     ∈ R d×d , σ i > 0, i = 1, • • • , d.Then equation 8 can be re-written asE L,M,H ({a η } n i=1 ) = 1 2T M L-1,M l,m=1 n i=1 n j=1 d k=1 (a i ) k (a j ) k σ 2 k (x m l ) ψ i (x m l )ψ j (x m l )(t l+1 -t l ) -2 n i=1 d k=1 (a i ) k (x m l+1 -x m l ) k σ 2 k (x m l ) ψ i (x m l ) , Here (x) k is the k th component of any vector x ∈ R d . We define α k = [(a 1 ) k • • • (a n ) k ] ⊤ ∈ R n , and A k ∈ R n×n as A k (i, j) = 1 2T M L-1,M l,m=1 n i=1 n j=1 (a i ) k (a j ) k σ 2 k (x m l ) ψ i (x m l )ψ j (x m l )(t l+1 -t l ), and b k ∈ R n as b k (i) = n i=1 (a i ) k (x m l+1 -x m l ) k σ 2 k (x m l ) ψ i (x m l ) .Then equation 8 can be re-written asE L,M,H ({a η } n i=1 ) = d k=1 (α ⊤ k A k α k -2α ⊤ k b k ). Since each α ⊤ k A k α k -2α ⊤ k b k is decoupled from each other, we just need to solve simultaneously A k αk -b k = 0, k = 1, • • • , d.Then we can obtain f (x) = n i=1 âi ψ k (x). However when D does not have a diagonal structure, we will have to resolve to gradient descent methods to minimize equation 8 in order to find the coefficients {a i } n i=1 for a total number of nd parameters. However, if a data-driven basis is desired, then we set H to be the space neural networks with the same depth, same number of neurons and same activation functions on the hidden layers. Furthermore, we find f from minimizing equation 7 using any deep learning optimizer such as Stochastic Gradient Descent or Adam from well-known deep learning packages.For d = 1, we also worked on polynomial drift function f = 2 + 0.08x -0.01x 2 . The estimation results are depicted in 4 and detailed in Table 4 . where f i : R 2 → R and x i ∈ R for i ∈ {1, 2}.For polynomial drift function f , we set f 1 = 0.4x 1 -0.1x 1 x 2 f 2 = -0.8x 2 + 0.2x 2 1 . Figure 5 , Figure 6a , 6b and Table 5 shows evaluation of the polynomial drift function estimation result.For trigonometric drift function f , we set f 1 = 2 sin(0.2x 1 ) + 1.5 cos(0.1x 2 ) f 2 = 3 sin(0.3x 1 ) cos(0.1x 2 ).  


Decision: Thanks to reviewers for their feedback. The paper proposes a trajectory based method for modeling dynamical systems described by a SDE. The authors are required to clarify the ambiguities addressed by reviewer oJVv, making sure that paper is easy to read. They can provide more details in appendices if they help better understand the main text. I vote accept (as poster) for this work.

Reviews: review: title: title: This work presents a noise guided trajectory based system identification method for inferring the dynamical structure from observation generated by stochastic differential equations.
main_review: main_review: This work presents a noise guided trajectory based system identification method for inferring the dynamical structure from observation generated by stochastic differential equations.
The paper is well written, and some conclusions are drawn from various numerical tests for highlighting the superior performance of our learning algorithm.
title: title: Review for Learning stochastic dynamics from data
main_review: main_review: The authors have presented a very neat study on learning stochastic dynamical equations from data. They have provided the mathematical framework and implementation details for real world data in a detailed manner. 

In section 4.2 the authors seems to have misplaced the Figure number. Please rectify it. The figures in section 4.1 and section 4.2 should be clearly identified. Currently it is misguided. The drift function used in the section 4.2 seems to be inconsistent with section heading of deep neural networks. I request the authors to look into that and avoid such discrepancies. 

The results section should be rewritten in a clear manner avoiding any misinterpretations.