Nonparametric inference of stochastic differential equations based on the relative entropy rate

The information detection of complex systems from data is currently undergoing a revolution, driven by the emergence of big data and machine learning methodology. Discovering governing equations and quantifying dynamical properties of complex systems are among central challenges. In this work, we devise a nonparametric approach to learn the relative entropy rate from observations of stochastic differential equations with different drift functions.The estimator corresponding to the relative entropy rate then is presented via the Gaussian process kernel theory. Meanwhile, this approach enables to extract the governing equations. We illustrate our approach in several examples. Numerical experiments show the proposed approach performs well for rational drift functions, not only polynomial drift functions.


Introduction
Dynamical properties of complex systems such as complexity and sensitivity are ubiquitous, which are significant contents in the research of stochastic dynamical systems [1][2][3]. Less well understood, however, is how to better quantify dynamical properties of stochastic dynamics. The relative entropy rate is an effective tool to characterise the complexity and sensitivity of dynamical properties [4,5]. Traditional methods are affected by factors such as unknown governing equations. Owing to the decline cost of data storage and computation, as well as the development of machine learning, data-driven discovery methodologies have made great progress. Combining with data-driven discovery methodologies, it is possible to compute the relative entropy rate from time series data.
There exist many different forms in regard to data-driven methods, such as parametric and nonparametric approaches [6][7][8][9][10]. The recent sparse identification of nonlinear dynamics method, which proposed by Brunton and co-workers [11], is a scriptures of parametric approaches. Opper [12] employs the parametric techniques to compute the relative entropy rate. Whereas some complex systems are too complicated to model precisely via parametric representation. The nonparametric techniques, for example neural networks [13,14], Gaussian processes [15,16] thus are applied in data-driven modeling. They outperform the parametric methods for obtaining the data-driven models with accuracy and extensiveness, although the parametric methods provide concrete manifestation of models.
The relative entropy or Kullback-Leibler divergence has recently reemerged of machine learning as a cost function describing the difference between two probability distributions [17,18]. It provides a precise characterization for approximating underlying probability distribution by probability distribution from data. The Kullback-Leibler divergence has been applied to infinite dimensional probability measures of stochastic processes, such as probability inference [19] and optimal control problems [20][21][22]. The relative entropy rate [23], which is the long-term average limit of the Kullback-Leibler divergence, is a distance measure of two stochastic processes. It plays an important role in quantifying dynamical properties. Dupuis et.al. [4] discussed the uncertainty and sensitivity bounds of stochastic dynamics from observation data using the relative entropy rate. Pantazis and Katsoulakis [5] applied the relative entropy rate as a suitable information-theoretic object to analyze the sensitivity of the probability distribution of stochastic processes in regard to perturbations in the parameters of the complex dynamics. Techakesari et.al. [24] use the relative entropy rate to design hybrid system filters in the presence of (parameterised) model uncertainty.
In this work, we propose a nonparametric approach to learn the relative entropy rate from observations of stochastic differential equations with different drift functions. Current efforts are mainly focused on deriving a simple estimator for relative entropy rate via the Gaussian process kernel theory. Meanwhile, this approach enables to extract the governing equations. Our method requires only one sample data of the stochastic differential equation. It performs well for rational drift functions, not only polynomial drift functions.
The remainder of this paper is structured as follows. In section 2, we introduce the relative entropy rate of stochastic differential equations and the variational formulation. In section 3, we propose a nonparametric technique for the calculation of the relative entropy rate based on one path data of stochastic differential equations. Some numerical examples are presented in section 4, followed in section 5 by conclusion.

The relative entropy rate
In this work, we consider a stochastic differential equation for the dynamics of a d-dimensional diffusion process X t ∈ R d given by with initial data X 0 = x 0 ∈ R d , where g(·) ∈ R d is the drift function, the diffusion function σ(·) is the d × k dimensional matrix and B t is the standard Brownian motion in R k .
The generator A of this stochastic differential equation (1) is [25] for all ϕ is in Sobolev space where H is the Hessian operator with H = ∇∇ T . Thus, the probability density p(x, t) of the solution process X t satisfies the Fokker-Planck equation with initial condition p(x, 0) = δ(x − x 0 ), where A * is the adjoint operator of the generator A in Hilbert space L 2 (R d ), given by Here, H(σσ T ϕ) is interpreted as matrix multiplication of H = ∇∇ T and σσ T ϕ (note that ϕ is a scalar function). We call p the stationary probability density of the solution process X t if it satisfies the stationary Fokker-Planck equation A * p = 0.

The relative entropy rate for stochastic differential equations
The relative entropy or Kullback-Leibler divergence between the probability measures P g and P r of two solution processes for stochastic differential equation (1) with different drifts g and r is defined as Here, the diffusion matrix of two solution processes is both D(x) . = σ(x)σ(x) T , p g (x, t) is the probability density of solution process with drift function g and ||u(x)|| 2 A . = u(x) · A(x)u(x) for some positive definite matrix A. Suppose that the stationary probability exists. The probability density p g (x, t) then converges to the stationary probability density p g as t → ∞. Hence, we consider the relative entropy rate as follows In our work, we suppose that the drift r(x) and the diffusion D(x) are known, but the expression form of the drift g(x) is unknown. We, however, want to estimate the relative entropy rate via the observation data of the process X t with drift g(x) on a large time T .
Based on equation (6), the estimator of the drift g(x) and the stationary probability density p g (x) are crucial to the calculation of the relative entropy rate. In order to simplify the estimation problem, we assume that the diffusion D(x) = σ 2 I and drift g(x) satisfies a potential condition such as g(x) = −∇ψ(x).
Then the stationary probability density fulfills p g (x) ∝ e − 2ψ(x) σ 2 and one can estimate the density from observation data.
We will introduce a different way to calculate the relative entropy rate based on the variational formulation and generalised potential condition. To be specific, suppose that the expression of drift g(x) is Then the relative entropy rate (6) becomes Moreover, the stationary Fokker-Planck equation for drift g is defined as where p g (x) is the stationary probability density and the operator A * r , corresponding to known drift r(x), can express as

Variational formulation
In this section, we will give a brief introduction about the variational formulation. Here, assume that the specific form of drift g is unknown. We first estimate the drift g from the observation data using the variational formulation for the stationary Fokker-Planck equation (9) through relative entropy rate. Suppose that the stationary probability density p g (x) is given and we search for a estimator of the drift g(x) by minimizing the relative entropy rate (6). Introducing a Lagrange multiplier function ψ(x), we may derive the drift g(x) from the following Lagrange functional The Fokker-Planck operator A * g is defined in (4) for drift g(x) and the operator A * r is in (10). For more details see appendix.
Furthermore, making a variation of Lagrange functional (11) with respect to g(x) − r(x), one can obtain g(x) − r(x) = D(x)∇ψ(x). Inserting this result back into (11), the variational representation of relative entropy rate for a unknown potential ψ is where the generator A r is adjoint operator of A * r , (10) which satisfies ψ(x)A * p(x)dx = p(x)Aψ(x)dx. In addition, the explicit expression of generator A r is We next introduce the variational bound for the Lagrange functional (12) where equality is achieved while ψ = ψ * . We are surprised to discover that, from the perspective of variational bound, the minimisation of variational representation (12) can give us a estimator for the potential ψ and also compute the relative entropy rate According to the form of Lagrange functional ε g [ψ], we all know that applying the observation data which are the ergodic samples of the process with the drift g, one can estimate every potential ψ under the stationary probability density condition. The estimator of the relative entropy rate, then, can be obtained.
Subsequently, we give a simple proof for variational bound (13) as follows ≥0.
The first equation comes from the variational representation (12) and the second is due to the algebraic formula (a − b) 2 = a 2 − 2ab + b 2 . Owing to the relationship of generator A r and its adjoint operator A * r , as well as the usage of integration by parts, the third equation is established. Inserting (9) into it later, we obtain the second to last equation.
Before moving on to next section, let us discuss another stochastic presentation (14). Consider The dP r+Dψ /dP r represents the Radon-Nykodim derivative of two probability measures. According to (5) and (6), we see Since d(P g , P r+Dψ ) ≥ 0, we thus have The equality holds in (18) when P g = P r+Dψ . In next section, we will present a nonparametric estimator for the potential ψ and the relative entropy rate.

Nonparametric inference for relative entropy rate
In this framework, our goal is to construct an nonparametric estimator to the the relative entropy rate. We however know from the above theory of the potential ψ and stationary probability density p g is essential to calculating the relative entropy rate. Hence, we first need to a smooth estimator for stationary probability density p g , and then compute the potential ψ from data using a nonparametric approach via the variational representation (12).
There are many density estimation methods, among which the kernel density estimation is commonly used. The drawback of the kernel density estimation is nevertheless that it neglects the temporal ordering of observation data, as well as the time scale between them. Here, we employ the empirical distribution to replace the exact probability p g ,p The data x 1 , x 2 , · · · , x n are random, ergodic samples of the process with the exact probability p g . On the other hand, one usually build a parametric estimator of the potential ψ, which is represented by a finite set of basis functions. However, this representation of parametric method is not sufficient to express many functions, such as rational functions. We therefore work on a nonparametric estimate, which is a more general approach to these functions.
We now demonstrate the nonparametric estimate method. Introduce a penalty term firstly, which is selected as a quadratic form 1 2 k ω 2 k /λ k , where ω k are the weights corresponding to the basis functions of the parametric representation and λ k are hyper-parameters, to regularize the potential estimator due to a limited amount of observational data. We can also consider the penalty term from the perspective of pseudo-Bayesian. The exp{− 1 2 k ω 2 k /λ k } is viewed as a Gaussian prior distribution over weights ω k . Assume that the parametric of the variational representation (12) is ε emp [ψ ω ]. At the same time, the exp{−Cε emp [ψ ω ]} can be interpreted as a likelihood, where C represents the proportional weight between the penalty term and the observation data. The potential function ψ can thus be treated as a Gaussian process according to the above interpretation. Inspired by the viewpoint of Gaussian process, we will convert the parameter representation form into the kernel function form using the kernel trick, which makes the information for ψ expressed fully. To this end, we define where λ k and ϕ k are the orthonormal eigenvalues and eigenfunctions, respectively. This is also the covariance kernel of the prior of a Gaussian process about the potentials ψ.
The regularized functional of the potential function ψ can be defined as via the kernel approach. The K −1 (x, x ) represents the inverse of the kernel operator. In addition, the penalty term in (21) can be proved to be equivalent to the reproducing kernel Hilbert space norm of the drifts ψ defined by the kernel K. We will further derive a specific representation of the drift function estimator.
The variation of (21) with respect to ψ yields Multiplying both sides by the operator K, one obtains where the generator acts on the kernel function K as From above equations (22) and (23), we see that if we know ∇ψ(x) at all observation data x = x i , the potential function ψ(x) can be calculated for all x. The key to this work, however, is to evaluate the gradient of the potential ψ(x). We next compute ∇ψ(x) at all observation data points via performing the gradient of equation (22) and setting x = x i . Thus, there are a series of linear equations We so far obtain the drift function value ∇ψ * (x i ) at every point. Further, substituting into these drift function value, we can calculate the relative entropy rate

Numerical Experiments
We begin our verification of nonparametric estimate approach by describing a few experimental examples. Without special emphasis, we apply the radial basis function kernel where the length scale l is the hyper-parameter. In this section, we also show the nonparametric learning results of the drift function in order to further illustrate the advantages of this method.
Example 1. Consider a scalar stochastic dynamical system with the polynomial drift term  where W t is the standard Brownian motion and the diffusion function is σ(x) = 1. The known drift term r(x) = 4x − 4x 3 corresponds to a stochastic double-well system. Another drift function, however, which generates the samples through long-term observation is g β (x) = 4x − 4x 3 + ∇ψ(x), where ∇ψ(x) = −β(x 2 + 2x + 1) depends on β. The drifts g and r are equal when β = 0.
In the numerical simulation, we take the time step δt = 0.001, and use the Euler scheme to sample the stochastic differential equation with drifts g β in regard to the different values of parameter β. Here, we apply N = 10000 data points x k = X t k , uniformly sampled from a trajectory with the observation time length T = 10 5 . Thus the discrete times are t k = kτ for k = 1, · · · , N with τ = 10. Employing this data points, we can estimate the drifts ψ(x) and their gradients, and further compute the relative entropy rate in term of (14). The numerical results will later be presented in Fig.1 and Fig. 2.
In Fig.1, we compare the accurate and the learning results with regards to potential function ψ(x) for some β values and discover that the learning results have a good performance. As shown in Fig.2, the exact relative entropy rate and estimators corresponding to different values of parameter β are plotted. We can see that the exact relative entropy rate agrees well with the estimate result directly obtained from nonparametric method.
Example 2. Consider the transcription factor activator (TF-A) monomer concentration of stochastic differential equation in the gene regulation system [26] where the diffusion function is σ(x) = 1 and W t is the standard Brownian motion. The known drift term here is r(x) = 0, and another drift which generates the samples through long-term observation is g β (x) = ∇ψ(x) = β((6x 2 /(x 2 + 10)) − x + 0.4) depending on β. If parameter β = 0, the drifts g and r then are equal.
To further validate our method, taking the time step δt = 0.001 in the numerical simulation, and using the Euler scheme to sample the stochastic differential equation with drifts g β in regard to the different values of parameter β, we generate a trajectory with the observation time length T = 10 5 . Here, we utilize N = 10000 data points x k = X t k , uniformly sampled from this generated trajectory. Thus the discrete times are t k = kτ, k = 1, · · · , N with τ = 10. Employing this data points, we can obtain the potentials ψ(x) and their gradients, and meanwhile compute the relative entropy rate in term of (14).  As shown in Fig.3 and Fig.4, the comparison regarding the learning and accuracy results of the drift function ψ(x) are given, as well as the exact relative entropy rate and estimators for different values of parameter β are plotted. We can see that the evaluation of the drift function ψ(x) is pretty and the exact relative entropy rate agrees well with the estimate result directly obtained from nonparametric method.
Example 3. In order to show the superiority of our method, we also consider the trigonometric polynomial function Here, W t is the standard Brownian motion and the diffusion function σ(x) = 1. Taking the known drift term r(x) = 0, and another drift which generates the samples through long-term observation g β (x) = ∇ψ(x) = β(sin(X t ) − sin(X t ). 3 ) depending on β. If parameter β = 0, the drifts g and r then are equal.
We take the time step δt = 0.001 in the numerical simulation, and then apply the Euler scheme to sample the stochastic differential equation with drifts g β in regard to the different values of parameter β, and generate a trajectory with the observation time length T = 10 5 . Here, we utilize N = 10000 data points x k = X t k , uniformly sampled from this generated trajectory, so that the discrete times are t k = kτ, k = 1, · · · , N with τ = 10. Employing this data points, we obtain the potentials ψ(x) and their gradients, and compute the relative entropy rate in term of (14). We show the results in Fig.5 and Fig.6.
The comparison regarding the learning and accuracy results of the drift function ψ(x) are given in Fig.5, as well as the exact relative entropy rate and estimators for different values of parameter β are plotted in Fig.6. We see that the evaluation of the drift function ψ(x) is pretty and the exact relative entropy rate agrees well with the estimate result.

Conclusion and Discussion
In conclusion, we have presented a nonparametric technique to learn the relative entropy rate which is applied to describe dynamical properties. The study of dynamical properties of systems is significant for stochastic dynamical systems. Here, we offer a good methodology for their study to learn from time series data. Moreover, the nonparametric method not only provides a precise estimation for our functions and the relative entropy rate, but also adapts to a wider range of functions types, such as rational function, compared with parametric approach. This performance has been well demonstrated in examples.   In addition, this work motivates a number of future extensions. Our framework can be extended in other Markov processes, including Markov chains or Lévy processes. They would be interesting in realistic application. On the other hand, the stationary density is worthy of attention because of the need for a long time observation. We can generalize this issue with finite time window T and evaluate further the marginal densities p g (x, t).
Moreover, under the induced probability measure Q by Girsanov theorem, Y t also satisfies the stochastic differential equation where B Q t = t 0 u(s)ds + B t is a Brownian motion with respect to the induced probability measure Q. Then by weak uniqueness of the solutions for stochastic differential equation (1) and (33), the likelihood function can be represented as To this end, we could use Itô formula to represent the stochastic integral as the a Riemann integral and the relative entropy is where D = σσ T and p g (x, t) is the probability density of the stochastic differential equation (1) with drift vector filed g. Assuming that p g (x, t) converges to the stationary density p g (x) for t → ∞, we shall consider the relative entropy rate D(P g , P r ) = 1 2 Adding the Lagrange multiplier term, equation (36) becomes our optimal objective functional (11).