Sparse inference and active learning of stochastic differential equations from data

Automatic machine learning of empirical models from experimental data has recently become possible as a result of increased availability of computational power and dedicated algorithms. Despite the successes of non-parametric inference and neural-network-based inference for empirical modelling, a physical interpretation of the results often remains challenging. Here, we focus on direct inference of governing differential equations from data, which can be formulated as a linear inverse problem. A Bayesian framework with a Laplacian prior distribution is employed for finding sparse solutions efficiently. The superior accuracy and robustness of the method is demonstrated for various cases, including ordinary, partial, and stochastic differential equations. Furthermore, we develop an active learning procedure for the automated discovery of stochastic differential equations. In this procedure, learning of the unknown dynamical equations is coupled to the application of perturbations to the measured system in a feedback loop. We show that active learning can significantly improve the inference of global models for systems with multiple energetic minima.


I. INTRODUCTION
Throughout the natural sciences, mathematical models are frequently formulated as differential equations.For example, with stochastic, ordinary, and partial differential equations (SDEs, ODEs, and PDEs).In physics, governing differential equations are often derived from first principles, for instance, from conservation of energy, momentum, and thermodynamic considerations.However, for complex systems studied, e.g., in biophysics, climate science, and neuroscience, first principles determining the system properties are typically not fully known.For example, because such systems are in a driven nonequilibrium state, highly nonlinear, and because dynamics may occur on multiple scales that are not well separated.In these cases, one can resort to phenomenological, effective descriptions that may result from some level of coarse graining and are based on experimental data.Recently, the increased availability of computational power has made it possible to construct such models in an automated fashion, which is known as data-driven discovery of governing equations.
Various approaches have been developed for inferring the differential equations that govern a non-linear dynamical system directly from measured data [1][2][3][4][5][6].In a popular approach called "symbolic regression", function libraries are employed to automatically extract the terms in a governing equation that best represents the measured data according to some optimization criterion [1,2].Recently, the use of sparse regression techniques for symbolic regression has received considerable scientific attention [3,4].In symbolic regression, the physical quantity , which is for illustration taken to be a scalar here, is assumed to obey an equation of the general form where ˇ can be, e.g., a time derivative ˇ = for ODEs and PDEs.F ( , x, , ) is an unknown function whose arguments x represent space coordinates while represents time and is a constant parameter.The aim of symbolic regression is to estimate the function F (. ..) from a data set z, which could be a measured sequence of values of at different time-space coordinates.The vector ž is either measured or estimated from z, e.g., with a discrete difference scheme.For inference of F (• • • ), a socalled "library" Θ(z) is constructed from a suitable set of functions of z, e.g., various powers of z, combinations of partial derivatives, or trigonometric functions.Assuming that the governing equation (1) can be expressed as a linear superposition of library terms, we write where ξ is a weight vector.The inference of the governing equation is thus reduced to a regression problem for the optimal ξ, given ž and Θ(z).In general, solving the inverse problem in Eq. ( 2) is not straight-forward since the matrix Θ should represent many equation terms and can have a large condition number (Θ).In Ref. [3], a method called sparse identification of nonlinear dynamics (SINDy) has been proposed.The method works iteratively.At each iteration, ξ is first obtained from a least-squares optimization involving Eq. ( 2) and ξ is subsequently thresholded such that values smaller than a cutoff are set to zero.The iteration is continued until convergence conditions are satisfied.SINDy has been shown to be a powerful and versatile method that is applicable for inference of various types of ODEs [3].However, the method requires the user to manually select the thresholds .For the identification of PDEs, an alternative algorithm called train sequential threshold ridge regression (TrainSTRidge) has been described in Ref. [4].This method is a variant of a least-squares optimization procedure for ridge regression called Sequential Threshold Ridge regression (STRidge).In STRidge, the vector ξ is first calculated by using ridge regression with a fixed regularization parameter.Then, all elements in ξ that have a smaller absolute value than a threshold are set to zero.Both, the regularization parameter and the threshold need to be provided by the user in STRidge.TrainSTRidge [4] employs L0 regularization and a training step to automatically determine the threshold while the regularization parameter remains to be set by the expert user.Conversely, a method called threshold sparse Bayesian regression, which also was employed for identification of PDEs [7], requires no input of a regularization parameters but some thresholds remain to be provided by the user.
The first aim of this work is to provide a method to solve the inverse problem associated with data-driven discovery of governing physical equations, Eq. ( 2), by combining a Bayesian approach with a automatic thresholding procedure.We call this method automatic threshold sparse Bayesian learning (ATSBL).Our algorithm does not require any manual fine-tuning of parameters to correctly infer governing differential equations from measured data.The method can be employed to identify ODEs, PDEs, and SDEs.
The case of SDEs requires particular attention, since the above-mentioned methods of equation inference are mainly designed for deterministic processes and some moderate amount of additive noise.The question of how to reconstruct the force fields for stochastic processes has been investigated in numerous studies, e.g., for application in soft matter physics and biophysics [5,[8][9][10][11][12][13]. Recently, sophisticated methods have been proposed for dealing with discretization and the inference problem in the context of SDEs for second-order dynamics [14][15][16][17].
Here, we focus on the use of symbolic regression for the inference of analytical expressions of SDEs of the overdamped Langevin-type.One approach to symbolic regression in this context is based on dividing the phase space into small hypercubes which are also called bins in the one-dimensional case.Average values of the state variables and of their derivatives are estimated in each hypercube and the regression is defined with respect to these averages [5].This kind of averaging generally depends on the chosen discretization and the averaging may lead to a substantial loss of information.Furthermore, application of this method to non-stationary processes requires a large ensemble of trajectories and considerable numerical effort to sample the time-dependent probability distribution in phase space.The difficulties related to the averaging in phase space motivate the investigation of the question to what extent the above-mentioned inference tools can be used in the context of noisy data without the need to perform ad hoc averaging, and, eventually, how the robustness of the inference methods may be improved in this context.We show that imposing Laplacian or Gaussian prior distributions on the inferred models is generally sufficient to identify the correct SDEs directly from trajectories without phase-space binning and we provide a comparison of the accuracy of results obtained with the two types of prior distributions.A remarkable performance of the Laplacian prior is demonstrated with several examples, including Brownian motion in time-dependent potentials.
A major challenge for the inference of SDEs is that the phase space is often sampled very inhomogeneously in available data.This problem is encountered, e.g., for systems where the long-term dynamics is dominated by transitions between different, locally stable states, while the short-term dynamics are dominated by fluctuations around individual stable states.In such cases, the inferred equation may be meaningful only locally, i.e, within the region covered by the measurement trajectory, and it may be a priori impossible to infer the global dynamics from a given data set.To enable an automatic inference of a global model under these conditions, we consider the question of how to design an external perturbation to the system, also called "control force", such that the state variables are forced to explore the full phase space in a shortened sampling time.Established Umbrella sampling routines used for this purpose rely on quadratic control forces and involve nontrivial design steps for the control force [18][19][20][21].See, e.g., Refs.[22,23] for alternative approaches.This kind of methodology has proved useful, e.g., in the context of computational studies of nucleation [24] and growth [25] processes.We develop an alternative adaptive control technique that recursively infers the governing equation and adapts the external control solely based on inferred equations.The adaptation loop consists of inference of the governing equation and a subsequent update of the control force such that it is directly opposite to the inferred force.No parameters need to be tuned for designing the control with this adaptive scheme.Using the adaptive control scheme, we demonstrate a substantial improvement of the inference of SDEs for several different simulations of Brownian motion.This work is organized as follows.The Methods section provides details on the the construction of function libaries and the casting of the inference problem into a system of linear equations.The inference algorithm is summarized and it is explained how Laplacian prior distributions can be used to impose the sparsity condition on the inferred models.In the Results section, the performance of the described method is illustrated by means of numerical examples and a comparison with previously described methods is presented.An adaptive sampling technique for improving the inference of SDEs is proposed and the usefulness of this approach is demonstrated.
Measurement data from a system of interest is presumed to be recorded as a time series of states, for example, a time-dependent position vector.In a data-driven approach to model a system, the data is used to automatically infer the a priori unknown dynamical equations that govern the observed process.In this work, inference is based on libraries of candidate functions for the governing equations.The data used for inference of differential equations is assumed to contain additive noise but no systematic errors.
For inference of ODEs, we generalize the introductory example for a scalar variable , Eq. ( 2), to a system with components that are assumed to be sampled with the same regular time interval for all ℓ ∈ {1 . . .} components.To distinguish discrete measurements from continuous variables, a subscript notation is employed in the following.The ℓ-th component measured in an ordered time series [ 1 , . . ., ] is written as Vectors or arrays containing multiple variable measurements, e.g., at different time points, are denoted with bold letters.The whole data can then be written in matrix form as Our approach also requires derivatives of the measured data.For simplicity, finite-difference approximations are used throughout this work.Approximate derivatives are denoted by the operator D ••• , which represents here a fourth-order finite central difference scheme.For example, a time derivative of the ℓ-th state component, z ℓ , at the -th timepoint is written as ℓ ( ) For the entire dataset, we write the time derivative as A governing ODE for the vector containing the trajectory of the ℓ-th state component may be written as a linear combination of elementary functions of all {z ℓ ′ }, e.g., as where the indices ℓ ′ , ℓ ′′ , and ℓ ′′′ cover the system dimensions, ⊙ denotes an element-wise product, and represents a constant.F ℓ can also depend on time, but we focus mostly on autonomous differential equations in the following.Since F ℓ (•) represents a linear combination of functions that can be calculated from the data, F ℓ (•) can be expressed with the help of a library matrix Θ(Z) multiplied with a sparse vector ξ ℓ .Thus, we obtain for Eq.(3) in discretized form where the terms of the library matrix Θ(Z) are calculated from the measurement data by evaluating the functions of {z ℓ ′ } and the non-zero elements of ξ ℓ characterize the dynamics of the system.Since Eq. ( 4) refers to ODEs, no derivative terms are contained in the library on the righthand side of the equation.Given D z ℓ and Θ(z), the aim is to calculate a sparse vector ξ ℓ with a minimal number of non-zero coefficients corresponding to a minimal number of terms necessary to describe the dynamics.For inference of PDEs, the library matrix Θ has to contain partial-derivative terms.Thus, data is required that allows the numerical estimation of derivative expressions with respect to two or more variables, for example, with respect to time and space.Usually, measurements therefore consist of discrete space-time series recordings of system variables.For example, an array Z representing the -dimensional state vector that is measured at time points in positions of one space coordinate is written as . . . . . . . . . . . .
with a sparse coefficient vector ξ ℓ to be determined.Note that a robust estimation of derivatives from noisy data is an important prerequisite for data-driven inference of ODEs and PDEs in this framework.The fourth order finite-difference approximations employed here may be supplemented or replaced with other methods, including denoising procedures and Gaussian process regression models.

B. Stochastic differential equations
We focus on Langevin-type SDEs to describe the time evolution of continuous, real state variables X( ), rep-resenting, e.g., the position of a Brownian particle in space [26].Trajectories, denoted by X( ), are timeordered sequences of values of space coordinates x.The general form of the considered SDEs is where we employ the Einstein sum convention and ℓ ( ) denotes the ℓ-th component of the system state at time .The trajectories X are calculated by making use of Ito's interpretation of stochastic integrals [26].The ℓ (X( ), ) represent the deterministic parts of the differential equations.For example, for a Brownian particle undergoing overdamped motion in the presence of conservative forces with a potential (x, ), we have The stochastic perturbations are assumed to result from a Wiener process with a noise source Γ ℓ ( ) and d ℓ = Γ ℓ ( ) d .The noise is assumed to obey a Gaussian distribution with a vanishing mean and a -correlated variance as respectively.The coefficient matrix ℎ ℓ,ℓ ′ in Eq. ( 6) scales the magnitude of the stochastic perturbations and is assumed to be diagonal, for simplicity.Further noise sources, e.g., resulting from an experimental measurement of a trajectory, are not explicitly considered throughout this work.The Fokker-Planck equation that corresponds to Eq. ( 6) and describes the evolution of a probability density function (x, ) is given by where the Fokker-Planck operator ˆ acting on (x, ) has the form The functions ℓ,ℓ ′ (x, ) are called Kramers-Moyal (KM) coefficients or drift and diffusion coefficients.Under the assumption of perfect knowledge of the trajectories X( ), the KM coefficients can be calculated from the incremental changes Δ ℓ ( ) ≡ ℓ ( + ) − ℓ ( ) in an infinitesimal time interval as (2) where . . .X( )=x denotes averages over the stochastic trajectories.The KM coefficients are related to the functions ℓ and ℎ ℓ,ℓ ′ in the Langevin equation as We consider only diagonal diffusion matrices, but the KM coefficients can depend explicitly on space and time.To estimate the coefficients, -dimensional trajectories ℓ, , ℓ ∈ {1, . . ., } are sampled with a small, regular time step at time points ∈ {1, . . ., }. Trajectory samples ℓ, are distinguished from the original stochastic variable ℓ ( ) by the index , representing the -th time point.Therewith, two new sequences are constructed as where is a small time step [5].The F (1)   ℓ and F (2)   ℓ are constructed with sample trajectories from random processes that are not differentiable.Use of these quantities for estimation of the KM coefficients in the spirit of Eq. ( 10) makes it necessary to first sample the stochastic process extensively to then approximate the average . . .X( )=x over different realizations of the process.Note that in basing the estimation on Eq. ( 10), we are neglecting two problems that occur for time series measured in the "real world".Firstly, measurement noise may render the assumption of a Markov process invalid on small scales [27].Secondly, the finite sampling interval cannot be made arbitrarily small in practice and therefore the estimated KM coefficients deviate systematically from the true coefficients [28,29].Procedures for correcting finite-sampling-time errors are available for various stochastic processes [30][31][32].While the focus of this work is on the inference problem for governing equations, finite sampling-time corrections should be employed in practical applications.
In the following, we employ two different methods for estimating the drift and diffusion coefficients.Firstly, a method is described in the next subsection that is based on binning of the data in phase space to produce histograms.Secondly, we compare the results obtained from data binning with results from direct estimation of the KM coefficients.

Estimation of KM coefficients from binned data
A classical method for the characterization of stationary, Markovian time series resulting from Langevin dynamics is based on binning of the trajectory data in space intervals [5,33,34].For this approach, we focus on problems with only one space dimension ( = 1).To estimate probability distributions, the data from multiple sample trajectories of the stochastic process is grouped into bins and the values in each bin are averaged as where ¯ , ¯ (2) , and ¯ (2) are bin-wise averages.The estimated probability for finding trajectory parts in the -th bin, , is normalized as =1 = 1 with 0 ≤ ≤ 1. Histograms resulting from data binning directly yield the curves for the drift and diffusion coefficients, see Refs.[5,33].The equations for the KM coefficients, (1) ( ) and (2) ( ), are inferred by finding analytical expressions for F(1,2) as functions of X.For this purpose, a library Θ ∈ R × is constructed from the binned data, where is the number of bins and is the number of terms in the library.For example, Θ( X) = [1, X, X ⊙ X, , X ⊙ X ⊙ X, sin( X), . . .] where ⊙ again denotes an element-wise product.If the library contains all the function expressions necessary to describe the KM coefficients analytically, the governing equations can be written as where W (1) and W (2) are two sparse vectors whose nonzero entries correspond to the library terms to be included in the sought-for analytical expressions for the KM coefficients.Equation (14a) yields (1) ( ) and Eq.(14b) yields (2) ( ).The inverse problems of finding optimal W (1,2) in Eq. ( 14) have the same form as the problem in Eq. ( 2).The binning of trajectories can produce significant errors in sparsely sampled regions, both in the interior and at the boundaries of the sampled phase space.We propose that the identification of SDEs can be improved by removal or filtering of the bins with high uncertainty.To substantiate this suggestion, we implement the inference procedure for unfiltered histograms and, additionally, implement a straight-forward extension that essentially consists of fixing a small probability threshold, below which all the data is discarded.While the probability threshold can can be determined in different ways, we employ here an automatic heuristic that was originally designed for edge detection in images [35].The procedure that is described in Ref. [35] consists of dividing the data according to probability thresholds to maximize the Shannon and Tsallis entropy, respectively.Maximization of the Shannon entropy produces thresholds that divide the data into "foreground" and "background", corresponding to signal-dominated and noise-dominated phase-space regions, respectively.The threshold value determining the "background" is then improved in a second step by maximizing the Tsallis entropy, whose pseudo additivity reportedly improves the analysis of data containing longrange correlations, see also Ref. [36].While we found that this method for determining a probability threshold is useful in practice, its theoretical underpinnings are to our knowledge not entirely clear.Thus, a manual selection of the probability threshold based on the results may be preferable in some cases.

Estimation of KM coefficients without data binning
A more direct approach for estimating the KM coefficients is based on the use of the trajectories F (1)   ℓ and F (2)  ℓ without binning or filtering.Since we do not intend to study transient initial dynamics, we mostly employ as input data a single, long trajectory generated from the stochastic process.For inference of the KM coefficients from the space-time trajectories, we construct a library Θ ∈ R × , where is the length of the trajectory and is the number of terms in the library.For example, where ℓ ′ covers all components of the stochastic process.Note that the library is constructed such that 1 ℓ, and 2 ℓ, at the -th time point depend only on functions involving coordinates { ℓ ′ , } ℓ ′ at the same time point.Thus, a velocity dependence or a history dependence of the estimators for the drift and diffusion coefficients is excluded.Under the assumption that the library contains all necessary terms describing the drift and diffusion coefficients, the coefficients for the ℓ-th component of the stochastic process can be inferred with where ℓ ∈ {1 . . .} and the vectors W (1,2) ℓ are non-zero in those entries that correspond to the terms in the libary that are required for the analytical description of the KM coefficient.The determination of the W (1,2)   ℓ is again an inverse optimization problem.

C. Solution of the inference problems with automatic threshold sparse Bayesian learning
For identification of the relevant library terms as, e.g., for Eq. ( 15), we propose an algorithm that we call Automatic threshold sparse Bayesian learning (ATSBL).The method consists of two main steps.First, the inverse problem is solved with an efficient algorithm called Bayesian compressive sensing using Laplace priors (BCSL) [37].Since the library is large, the solution vector generated by the BCSL algorithm typically still contains quite a few non-vanishing but small entries.Therefore, in a second step, the negligible contributions to the resulting governing equations are removed by an automatic thresholding procedure [3,5,7].These two steps of the method are detailed below.
1. Bayesian compressive sensing using Laplace priors (BCSL) We consider a generic linear equation system involving a given vector g and matrix Φ and an unknown, sparse vector w as where the vector s represents noise or measurement errors.Here, w can be thought of as a solution vector appearing in an iterative solution procedure for Eq. ( 5), Eq. ( 14), or Eq.(15).Various methods can be used to calculate sparse solution vectors w from Eq. ( 16).In particular research on compressive sensing, which deals with the reconstruction of sparse signals from underdetermined systems, has yielded broadly applicable, efficient methods for finding sparse solution vectors w.Among these are Bayesian methods based on the relevance vector machine (RVM) [38,39].Very sparse result vectors are obtained if a Laplace distribution is used as a prior probability distribution for w.Here, we employ a method called Bayesian compressive sensing using Laplace priors (BCSL) [37].Specifically, we employ a variant of BCSL that interatively calculates approximate solutions, which is very computationally efficient and yields accurate results for our type of applications.Briefly, the mathematical basis of BCSL is as follows, see Ref. [37].The method is based on a three-stage hierarchical model.It is assumed that the errors s are drawn from a zero-mean Gaussian distribution with unknown variance 1/ .Therefore, the likelihood function for finding a vector g is given by The unknown vector w is assigned a prior distribution, which represents our knowledge on the nature of this quantity.To encode sparsity, one would like to employ a Laplace prior (w| ) = /2 exp(− | |/2) with a hyperparameter .However, the evaluation of integrals using this choice of a Laplace prior is not readily achieved since the Laplace prior is not conjugate to the Gaussian likelihood, Eq. (17).Therefore, an auxiliary vector of non-negative hyperparameters γ with the same dimension as w is employed to express the prior as the convolution of the two different distributions (w|γ) = Π exp (− 2 /(2 ))/ √ 2 and (γ| ) = Π [ exp (− /2)/2].These two distributions together result in a Laplace prior after marginalizing out γ as see Ref. [40].Overall, the joint probability density results as where the parameters and are both assumed to obey Gamma distributions.To infer values for the most probable solution vector w as well as the hyperparameters, an evidence procedure is employed wherein the posterior probability (w, γ, , |g) is maximized with respect to w, γ, , and , given the data.By making use of the expression ) together with Eq. ( 19), w is determined by simply maximizing (g|w, ) (w|γ).This calculation yields for the result vector the expression w * = ΣΦ T g with Σ = ( Φ T Φ + Λ) −1 with being a square matrix that contains the (1/ ) on the diagonal and is zero otherwise.This step corresponds to a Ridge regression that depends on the unknown values of γ, , and .Determination of these hyperparameters proceeds by maximizing with respect to γ, , and .Here, (γ, , , g) is calculated from the right hand side of Eq. ( 19) by integrating out w.With the fast, approximate version of BCSL, the equations determining the optimal values of the hyperparameters are solved iteratively, where only one entry of the vector γ is adjusted in every step.

Automatic thresholding
Solution of the inverse problem (16) with BCSL typically yields vectors w that contain only a few large entries, but also a number of very small, non-zero entries.Removal of these negligible entries is desirable and we improve the solution sparsity with an iterative thresholding procedure [4].The pseudocode 1 illustrates how the thresholding procedure proposed in in Ref. [4] is combined with BCSL proposed in Ref. [37].Briefly, the thresholding algorithm works as follows.The input is given by g, the library matrix Θ, an initial increment tol for the threshold , and the maximum number of iterations iters .The data g and Θ is spilt into two parts for training and test, respectively.Usually, 80% of the data is used for training and 20% for testing.Thresholds are calculated iteratively from the training data and the validity of the thresholds is evaluated based on the error resulting from their application to the test data.The core part of the algorithm is a loop for iterative calculation of the sparse vector w and the threshold .In each iteration step, the approximate, fast BCSL routine is first employed to obtain an estimate of w from the training data.The quality of this solution estimate is evaluated by calculating the resulting error with the test data where the penalty factor of the solution norm is chosen = 10 −3 (Θ) as suggested for the original algorithm [4].
If the error of the current solution is smaller than the error of previous iterations, the new solution is accepted and the threshold is increased.In the opposite case, the threshold is decreased and the increment tol is refined.The final solution w best is the sparse vector that determines the terms in the governing differential equations, SDEs, ODEs, and PDEs.

=
+ tol end } return w best Algorithm 1: Pseudocode for automatic threshold sparse Bayesian learning (ATSBL), which combines BCSL [37] and TrainSTRidge [4] to achieve parameter-free inference of highly sparse solutions to inverse problems.

D. Quality score for identified governing equations
The error of the inference procedure can be directly quantified by comparison of the results with a known set of original differential equations in test cases.For this purpose, we define the deviation of identified coefficient (DIC) as where every is a coefficient of one term in the identified equation and ′ is the related coefficient in the original equation that was used to generate the test data.Here, at least one of the coefficients in each pair { , ′ } is required to be non-zero and the sum only runs over these coefficients.represents the number of these coefficients.The DIC lies in the range [0, ∞] where 0 indicates a perfectly identified equation.

A. Inference of SDEs from noisy trajectories
We now illustrate the data-driven identification of SDEs by the example of overdamped Brownian motion of a particle inside a one-dimensional double-well potential with coordinate .The drift and diffusion coefficients of this system are given by (1) The trajectory data that is to be used for inferring the governing equation is generated by integrating the Langevin equation with the Euler-Maruyama method.A trajectory X is shown in Fig. (1-a-i) (10 6 time steps).The trajectories F (1) and F (2) are shown in Fig.
(1-a-ii, iii).To visualize the -dependence of the estimator for the drift coefficient, we plot F (1) against X, see Fig. (1-ci).Similarly, F (2) is plotted against X as estimator of the diffusion coefficient (2) ( ) in Fig. (1-c-ii).Both plots exhibit large fluctuations around the true drift and diffusion coefficients and the resulting averages are clearly prone to errors, particularly at the boundaries of the sampled domain.
Using the trajectory data, we next construct a library consisting of 11 terms for the drift coefficient and 6 terms for the diffusion coefficient as illustrated in Fig. (1-b).Then, we employ ATSBL to identify W (1) and W (1) directly from the trajectory without binning.The identified -dependent functions for the drift and diffusion coefficients are plotted in Fig. (1-c-i, ii).They agree well with the original functions used for creating the data.The identified equations with estimated uncertainties are shown in Fig. (1-c).
The main distinction of ATSBL as compared to established inference techniques is the assumption of a Laplacian distribution for the prior of the library coefficients.The more direct, albeit theoretically less sparsitypromoting procedure is to employ a Gaussian prior, corresponding to a Ridge regression with fixed regularization parameter, for inference of the solution vector w in Eq. ( 16) prior to automatic thresholding, as done, e.g., in Ref. [4].To compare the performance of these two approaches for inference of SDEs, we evaluate the deviation of the identified coefficients, DIC, as a function of the number of data points used for inference.The results shown in Fig. (1-d) indicate that the Laplacian prior is preferable over the Gaussian prior since it requires less data and results in a smaller DIC.To further establish the robustness of ATSBL, we consider the convergence of the iterative thresholding procedure for each of the two prior distributions.The result shown in Fig. (1-e-i, ii) demonstrate a better convergence achieved in the case of the Laplacian prior.For both, Gaussian and Laplacian prior distributions, the threshold and the error oscillate during the iteration process, which is due to the adaptive step size during the thresholding.

B. Inference of SDEs with time-dependent drift coefficient
In the previous section, an example is provided of how the KM coefficients can be obtained by performing a regression directly with the trajectory data.The direct use of the trajectory data becomes particularly important for the treatment of the more complex situation of a time-varying force.In such a situation, the probability distributions change over time and a histogram-based approximation of the dynamic distributions can be technically challenging and requires the availability of many sample trajectories for the same conditions.In order to explore the validity of our approach in this situation, we consider the example of a particle diffusing within a timedependent one-dimensional potential.The drift and diffusion coefficients of this system are given by (1) The potential has a double-well shape, where the positions of the two minima vary in time.The two minima start at separate positions and merge periodically into one minimum before separating again to reach the initial positions.Each time the potential wells come close to each other, the transition probability becomes large and the particle is likely to change from one well to another.This gives rise to stochastic oscillations between the two potential wells.To test whether the underlying equations can be inferred with a library constructed from a single trajectory, we assume that the frequency at which the potential changes is known.Inference of this frequency from the data is in principle also possible, but requires excessive computational power since high-order terms with explicit time-dependence must be accounted for in the library.We construct a library consisting of a time-dependent and a time-independent part.The first half of the library is simply a polynomial basis, the second half corresponds to the polynomial basis multiplied with a cos( ) factor.The results of the inference procedure are shown in Fig. (2-a).The inferred equation is in agreement with the correct equation.For illustration, we plot snapshots of the drift and diffusion for the original equations and the inferred equations in Fig. (2-a-ii, iii).Note that the inference procedure for first-order SDEs shows a remarkable performance even though only one sample trajectory is used for inference.

C. Data binning for inference of SDEs from short trajectories
An inference method that relies on direct use of sample trajectories for a regression can become unreliable when confronted with short trajectories in an inhomogeneous force field.In such a situation, we find that it is more appropriate to employ data binning.We illustrate this procedure with Brownian motion of a particle in a one-dimensional double-well potential where the diffusion coefficient depends on space.The drift and diffusion coefficients of the model are given by (1) We first consider short trajectories that have 2 • The distributions approximated by the binned data clearly deviate from the known functions (1) ( ) and (2) ( ) in undersampled regions.Therefore, the binned data is filtered to remove data points with high uncertainty.This filtering is done as described in the Methods section by discarding bins below a probability threshold * that is determined by entropy maximization [35], see Fig. (2-b-iv).To assess if binning and filtering is also beneficial for inference from long trajectories, we also use data from trajectories with 2 • 10 7 time steps.Figure (2-b-vii) shows the errors of the identified coefficients.
For short trajectories, binning is advantageous in combination with a filtering procedure to suppress data with high uncertainty.The reason for this result can be understood from inspection of Fig. (2b-v, vi), where the inferred functions match the correct functions only in the most populated regions of phase space.Thus, the exclusion of data points with high uncertainty prevents overfitting and improves the inference of the underlying dynamical equations if the trajectory is not long enough to allow a sufficient sampling of the whole phase space.Conversely, data binning with or without filtering is disadvantageous for the analysis of long trajectories that sample the whole phase space, see Fig.

D. Active sampling improves the identification of SDEs
We have so far restricted our attention to the extraction of estimates from data that was generated prior to the analysis, e.g., in experiments.Thereby, we have assumed that the size of the data set is large enough to allow some form of inference of the governing equations.In a different scenario one might have the ability to perturb the studied system, either in a computer simulation or in an experimental setup, while simultaneously recording the data.Then, one may enhance the sampling efficiency by means of an appropriately designed perturbation that is applied to the system.Generally, this methodology is expected to be useful whenever the system exhibits an energy landscape with multiple local minima that can trap the trajectory for long times.We describe an adaptive control method where the inference of the dynamical equations together with a simultaneous perturbation of the system recursively results in a global exploration of the phase space to provide sufficient sampling everywhere.
Since the probability distribution tends to be peaked around local energy minima, the dynamical equations can be estimated locally near these minima.To take advantage of this local estimation while iteratively extending the sampled region, we re-sample repeatedly while applying in each sampling round a control force that is opposite to a force from the system that is estimated locally from previous rounds.The difficulty with a straight-forward application of this method is that the control force can admit large deviations away from the initial estimation region.This effect produces large errors, slows down convergence, and may even lead to divergence problems.We overcome this problem by weighting the control force with a Gaussian distribution such that the control force vanishes away from the current estimation region.Such a control force expels the trajectory from the local minimum where the estimation has been performed and the trajectory eventually reaches another local minimum.
The method, which we call automatic iterative sampling optimization (AISO), is illustrated in Fig. (3-a, iiii) and the pseudocode is provided in Algorithm 2. At each iteration, the underlying dynamical equations are estimated from the data accumulated during all previous iterations.The negative of the inferred drift term is employed locally as control force.The center and width of the Gaussian weight of the control force is calculated only from the mean and standard deviation of the trajectory of the previous step.Thus, we define our control force acting on the component ℓ as where the index indicates that values are to be taken at the iteration number ; ℓ and ℓ stand for the mean and variance of the trajectory extracted from the step in each iteration.After a sufficient number of iterations, the data points accumulated from all iterations are combined and the equation of motion is extracted from the accumulated data.This procedure is repeated for a predefined number of iteration steps.For the examples presented in the following, the iteration step number has been fixed to = 10, since convergence has been achieved within less than 10 steps in these cases.
For a first demonstration of our method, we employ a three-well potential ( ) = 6 − 6 4 + 0.5 3 + 8 2 with a constant diffusion coefficient for simulating the trajectory of a particle in one dimension.The drift and diffusion coefficients are Next, we also consider a two-dimensional drift field, consisting of a radially symmetric component and a shear component in the , -plane.The drift and diffusion coefficients are given by Using these driving forces, we simulate trajectories with As the algorithm proceeds through more iterations, the coefficients of the control potential approach the coefficients of the correct drift field, and the expulsion from each local minimum becomes more efficient, Fig. (3b-iii, vii).This results in an enhancement of rare events where the particle crosses the saddle points, as illustrated in Fig. (3-b-i, iv) by the controlled and uncontrolled trajectories.The error is quantified by calculating the coefficients ˜ (1) ( ) and ˜ (2) ( ) in each iteration.The DIC reduces from 1 to nearly 0.01 during the iterations, see Fig. (3-b-ii, v).Thus, the terms of the identified equations approach those given in the original equations.

E. Identification of ordinary and partial differential equations
It is next shown that the sparse inference scheme based on Laplace priors that is implemented with ATSBL can Function: AISO() c 0 = 0; %Initialize control force to zero N = 10; %Fix number of iteration steps i = 0; %Initialize iterator to zero while < do i= i+1; %For all degrees of freedom ℓ, ℓ ′ ∈ {1 . . .} %Add control force and generate data ℓ ( ) = ℓ (X( ), ) + −1 ℓ + ℓ ′ ℎ ℓ,ℓ ′ (X( ), ) ℓ ′ ( ) %Subtract control force and collect data ℓ ( ) ← ℓ ( ) − −1 ℓ %Concatenate data Θ ← (Θ, Θ( ℓ ( )) g ← (g, ℓ ( ) ) %Estimate libary coefficients with ATSBL w = ATSBL(Θ, g, tol , iters ) end return w Algorithm 2: Pseudocode for identification of SDEs with automatic iterative sampling optimization (AISO) also be used for data-driven discovery of ordinary and partial differential equations.The identification of ODEs from trajectory data is demonstrated with a Lorenz system, which is a paradigm for chaotic behaviour [41].The Lorenz equations are given by where the subcript denotes a time derivative and the parameters are fixed as = 10, = 28, and = 3/8.We numerically integrate these equations to obtain a trajectory as shown in Fig. (4-a).The chaotic system involves two attractors.For data-driven system identification, we utilize three identical libraries Θ for each of the variables, , and .Θ is constructed from the simulated trajectory and includes 56 terms containing up to fourth powers of all variables.Time-derivatives are calculated using fourth-order central-difference approximation.The general ODEs constructed from the library as in Eq. ( 4) are represented by three linear equation systems.The estimated equations resulting from an application of the inference procedure to noise-free data have small errors that are in magnitude comparable to the time step Δ = 2 • 10 −4 , see Fig. (4-b).The same inference procedure is then repeated for a trajectory with additive Gaussian noise.The standard deviation of the noise in each coordinate is chosen to be 2 % of the standard deviation of the noise-free data in the same coordinate.For this case, the ODEs identified with ATSBL still contain all the correct terms and the errors in the inferred system parameters are in the percent range, see Fig. (4-b).
Finally, we demonstrate data-driven discovery of PDEs with ATSBL.Reaction-diffusion equations have attracted interest as prototypic models for pattern formation in biochemical systems, where constituents are locally transformed into each other through chemical reactions and transported in space by diffusion.Here, we consider the popular − system, given by where is equal to 2. A two-dimensional, planar, rectangular area with periodic boundary conditions is considered.The initial values of and are shown in Figs.
(5a-i, ii).The reaction-diffusion equations are solved numerically by using a spectral method.Snapshots of and are shown in Figs.
(5-a-iii, iv).For inference of the governing PDEs, a library matrix Θ is constructed containing 35 terms each for and .Then, using ATSBL, the reaction-diffusion equations are inferred from the simulated data, as illustrated in Fig. 5(b).For noise-free data, the identified equations deviate from the original equations only at the fourth decimal place and this error is due to discretization.However, if and are corrupted with additive noise, identification of the correct PDEs becomes challenging [4].In Ref. [4], it has therefore been suggested to include a denoising step prior to the inference step.Accordingly, we employ a curvelet denoising method [42], which permits reconstruction of the reaction-diffusion equations from data with 2% noise with ATSBL, as illustrated in Fig. 5(b).

IV. SUMMARY AND OUTLOOK
Data-driven, automatic discovery of governing equations has become a viable tool for studying complex systems if first-principle derivations are intractable, e.g., for biological systems or epidemiological data.The aim is here generally to construct an analytical model that characterizes the observed dynamics and extends to parameter-and phase space regions that are hard to access experimentally.
Our main contribution is an inference method that makes use of Laplacian prior distributions in a Bayesian framework to find a minimal set of governing equations without the need for user input.We establish the validity of this approach and compare it to other methods.Regarding data-driven discovery of Langevin-type SDEs, we show that the proposed sparse method converges faster than other methods based on ridge regression.Maximum likelihood methods for the estimation of parameters in SDEs are not considered here, see Ref. [8] for an introduction to those methods.For the studied Langevin SDEs, we find that a binning of the trajectory data for inference of the drift and diffusion coefficients is only advantageous if the phase space is sampled sparsely.In that case, the error of the inference procedure can vary strongly in phase space since the relative uncertainties of the probabilities vary.A filtering procedure consisting of the exclusion of data with high uncertainty results in a significantly improved inference accuracy.
Next, we investigate how active-learning procedures can be useful in situations where the inference of globally valid equations becomes difficult because most trajectories are trapped in local potential minima.This problem can be solved with well-established umbrellasampling methods where quadratic bias potentials are employed to reduce the energetic barriers in the original potential landscape [18,43].However, an appropriate parameterization of such bias potentials can be challenging.For example, if the additional potentials are intended to smoothen an unknown, rough potential landscape.Instead, we employ data-driven identification of governing equations for calculating time-dependent external perturbations that force the trajectory to explore the full phase space.The main feature of our method is that the parameters that determine the control force correspond to the parameters that define the potential landscape.The combination of iterative inference with system perturbations can significantly improve the speed and accuracy of the overall inference procedure.We therefore hope that the suggested active learning procedure will extend the applicability of data-driven methods, in particular in the context of computer simulations.
A central challenge related to the improvement of the library-based methodology for identification of analytical models is to find automated approaches for tailoring the employed function space to the problem at hand.Recent methodological advances suggest that a possible solution is the integration of physical constraints, such as symmetries, conservation laws, or even thermodynamics, into a generic framework for statistical learning of governing equations [44].Data-driven identification of analytical models thus has the potential to become a popular tool for closing the gap between non-parametric, empirical modeling and first-principles-based modeling in the coming years.  1and F (2) generated with discrete differences from the same trajectory.(b) The library matrix Θ is constructed by evaluating a given set of functions for all values of the trajectory.Thereby, one obtains linear equation systems that relate the known sequences F (1,2) to unknown, sparsely populated coefficient vectors W (1,2) .The determination of the non-zero entries of W (1,2) yields a set of library functions that together describe the drift and diffusion coefficients (1) and (2) .(c) Exemplary results of the inference procedure.Despite the large noise amplitude, accurate predictions can be made directly from the trajectory data.  (1and (2) .For long trajectories, data binning does not reduce the error.

10 5
time steps with one time step being Δ = 5 • 10 −3 .Parts of the trajectories on the potential maps are shown in Fig. (3-b-i, iv).Results for the intermediate iteration steps are shown together with the drift field in Fig. (3b-iii, vi).

FIG. 1 .
FIG.1.Data-driven discovery of a one-dimensional SDE with automatic threshold sparse Bayesian learning (ATSBL).(ai) Trajectory of a particle undergoing overdamped diffusive motion in a double-well potential (10 6 time steps).(a-ii,a-II) Values of the F(1) and F(2) generated with discrete differences from the same trajectory.(b) The library matrix Θ is constructed by evaluating a given set of functions for all values of the trajectory.Thereby, one obtains linear equation systems that relate the known sequences F(1,2) to unknown, sparsely populated coefficient vectors W(1,2) .The determination of the non-zero entries of W(1,2) yields a set of library functions that together describe the drift and diffusion coefficients(1) and(2) .(c) Exemplary results of the inference procedure.Despite the large noise amplitude, accurate predictions can be made directly from the trajectory data.(d) Comparison of the use of a Laplacian and Gaussian prior distribution in the inference procedure.The deviation of the identified coefficient (DIC) for the drift coefficient is plotted against the number of data points used for training.The Laplace prior in ATSBL decreases the error and reduces the required sample size.(e) Convergence rate of the thresholding procedure for Laplacian and Gaussian prior distributions.(e-i) Laplace priors result in fast threshold convergence.(e-ii) The error defined in Eq. (22) decreases during the iterations.Errors achieved with Gaussian-and Laplacian priors are comparable.

FIG. 2 .
FIG.1.Data-driven discovery of a one-dimensional SDE with automatic threshold sparse Bayesian learning (ATSBL).(ai) Trajectory of a particle undergoing overdamped diffusive motion in a double-well potential (10 6 time steps).(a-ii,a-II) Values of the F(1) and F(2) generated with discrete differences from the same trajectory.(b) The library matrix Θ is constructed by evaluating a given set of functions for all values of the trajectory.Thereby, one obtains linear equation systems that relate the known sequences F(1,2) to unknown, sparsely populated coefficient vectors W(1,2) .The determination of the non-zero entries of W(1,2) yields a set of library functions that together describe the drift and diffusion coefficients(1) and(2) .(c) Exemplary results of the inference procedure.Despite the large noise amplitude, accurate predictions can be made directly from the trajectory data.(d) Comparison of the use of a Laplacian and Gaussian prior distribution in the inference procedure.The deviation of the identified coefficient (DIC) for the drift coefficient is plotted against the number of data points used for training.The Laplace prior in ATSBL decreases the error and reduces the required sample size.(e) Convergence rate of the thresholding procedure for Laplacian and Gaussian prior distributions.(e-i) Laplace priors result in fast threshold convergence.(e-ii) The error defined in Eq. (22) decreases during the iterations.Errors achieved with Gaussian-and Laplacian priors are comparable.

FIG. 3 .FIG. 4 .FIG. 5 .
FIG.3.Active learning with automatic iterative sampling optimization (AISO).(a) Schematic presentation of automatic iterative sampling optimization for the case of Brownian diffusion (a-i).Initially, the particle is trapped in a local energetic minimum and the functional form of the potential can therefore only be inferred locally.(a-ii) After the first iteration step, the potential hypersurface near the estimated minimum is flattened and the particle can thus explore other regions of phase space.The same procedure is repeated iteratively and the control is always applied at the minimum estimated during the previous iteration.(a-iii) Schematic representation of the main feedback control loop.(b-i) Trajectory of a particle undergoing Brownian motion in a one-dimensional three-well potential.The green curve shows a trapped trajectory while the blue curve shows a trajectory in the presence of control forces.(b-ii) The deviation of the inferred coefficients (DIC) decreases during the iterations.(b-iii) The identified drift field converges to the correct function during the iteration.(b-iv) Trajectory of a particle undergoing diffusion in a two-dimensional force field.The green curve exemplifies a trapped trajectory for plain sampling.The blue curve shows an example of a trajectory in the presence of control forces.The color of the background only represents part of the force field, namely a Mexican hat potential ( , ) = −( 2 + 2 )/2+ ( 2 + 2 ) 2 /4 that generates radial forces.(b-v) Evolution of the of the DIC during the iterations.(b-vi) Streamlines of the identified drift field (pink) and streamlines of the correct drift field (black) after the first iteration step.(b-vii) Streamlines of identified drift field (pink) and streamlines of the correct drift field (black) after the tenth iteration step.The identified force field at the end of the iteration closely matches the original one.