Symbolic Discovery of Stochastic Differential Equations with Genetic Programming
Authors: Sigur de Vries, Sander W. Keemink, Marcel A. J. van Gerven
Organizations: Department of Machine Learning and Neural Computing, Donders Institute for Brain, Cognition and Behaviour, Radboud University, Nijmegen, the Netherlands
Automated scientific discovery aims to improve scientific understanding through machine learning. A central approach in this field is symbolic regression, which uses genetic programming or sparse regression to learn interpretable mathematical expressions to explain observed data. Conventionally, the focus of symbolic regression is on identifying ordinary differential equations. The general view is that noise only complicates the recovery of deterministic dynamics. However, explicitly learning a symbolic function of the noise component in stochastic differential equations enhances modelling capacity, increases knowledge gain and enables generative sampling. We introduce a method for symbolic discovery of stochastic differential equations based on genetic programming, jointly optimizing drift and diffusion functions via the maximum likelihood estimate. Our results demonstrate accurate recovery of governing equations, efficient scaling to higher-dimensional systems, robustness to sparsely sampled problems and generalization to stochastic partial differential equations. This work extends symbolic regression toward interpretable discovery of stochastic dynamical systems, contributing to the automation of science in a noisy and dynamic world.
Figures & tables
Figure 1: Overview of genetic programming of stochastic differential equations. (a) Stochastic time series are observed to which stochastic differential equations are fitted with symbolic expressions. The resulting stochastic differential equation offers interpretability and generative modelling. The symbolic equations for the drift f(x) and diffusion g(x) are optimized with genetic programming. (b) In genetic programming, trees are adapted with crossover that swaps subtrees, illustrated with the exchange of the red and blue subtrees. (c) Besides crossover, mutation also changes the structure of individual trees. Various aspects can be mutated, where this example shows the change of an operator in green. (d) Examples of systems on which the proposed method can be applied to. The curves in each subpanel illustrate the behaviour of the state following stochastic system dynamics.
Figure 2: GP-SDE accurately recovers the governing equations of stochastic dynamical systems. The mean squared error (MSE) between the true and learned functions for the drift and diffusion on test data are presented for ten different seeds. The average of the ten runs in indicated with a black cross for every method. Furthermore, for each run, a star indicates that the correct structure of the equations were found for all variables, while a circle indicates that the structure of at least one equation was incorrect. The methods include Kramers-Moyal expansion followed by sparse regression (KM-SR), genetic programming of ordinary differential equations (GP-ODE), which only learns a function for the drift, and genetic programming of stochastic differential equations (GP-SDE). (a, b) Results for the double well problem with an additive, linear and non-linear multiplicative diffusion term. (c, d) Results for the van der Pol oscillator. (e, f) Results for the Rössler attractor. (g, h) Results for the Lorenz96 model with 5, 10 and 20 variables. (i, j) Results for the Lotka-Volterra model with a sampling rate of 0.02, 0.2 and 0.5. Both GP-ODE and GP-SDE are extended with multi-step integration, labelled as GP-ODE-MS and GP-SDE-MS respectively.
Method
fx(x,y,z)
gx(x,y,z)
fy(x,y,z)
gy(x,y,z)
fz(x,y,z)
gz(x,y,z)
True
−x−y
0.100x
x+0.200y
0.100y
xz−5.700z+0.200
0.100z
GP-ODE
−0.999x−0.999y
−
x+0.196y
−
xz−5.790z+0.282
−
KM-SR
−0.950x−1.486y
0.306
1.003x+0.198y
0.146y
0.868xz−4.849z+0.157
0.064z+0.055
GP-SDE
−x−0.981y
0.083x
x+0.191y
0.083y
xz−5.670z+0.188
0.091z−0.008
Table 1: The equations corresponding to the true system and best solutions discovered of the Rössler attractor with genetic programming for ordinary differential equations (GP-ODE), Kramers-Moyal expansion and sparse regression (KM-SR) and genetic programming for stochastic differential equations (GP-SDE). For each variable x , y and z , the equation for the drift and diffusion are shown. GP-ODE only evolves a function for the drift, therefore no diffusion equations are shown.
Figure 3: Simulation of the discovered models for the Rössler attractor. (a) Given a fixed initial condition, three trajectories of the Rössler attractor are simulated with different Wiener processes (black). The best ordinary differential equation evolved by genetic programming (GP-ODE) is integrated from the same initial condition (orange). (b) The same three trajectories of the true system are shown, together with the mean (green) and standard deviation (shaded area) computed over 50 trajectories sampled from the best model found by method based on the Kramers-Moyal expansion and sparse regression (KM-SR), where the initial condition was fixed but different Wiener processes were used in each trajectory. (c) Same approach as in (b), but with the best stochastic differential equation discovered with genetic programming (GP-SDE). For each method, the corresponding equations are presented in Table 1 .
Table 2: The system corresponding to the true stochastic partial differential equations and the solutions evolved with GP-SDE are presented.
Figure 4: GP-SDE can recover stochastic partial differential equations. (a) Evolution of the Fisher-KPP equation with respect to x and time, given the ground truth equation and identified system with GP-SDE, presented in Table 2 . (b) Snapshots of the evolution of the two-dimensional heat transfer with respect to x , y and time. Again, the true and learned system with GP-SDE are presented.
Figure 5: Runtimes of the methods. The computational time of GP-ODE, KM-SR and GP-SDE on the Lorenz96 model with increasing dimensionality. The runtimes are averaged over the ten seeds for every experiment, after running one seed to compile the algorithms. KM-SR is evaluated with four and sixteen bins ( b ).
Appendix figures & tables2 assets
Supplementary material from the paper’s appendix.
Appendix
Environment
b
β
α
λ
d
Double well problem
[5, 10, 25, 50, 100]
[1, 10, 20]
[0.001, 0.01, 0.1]
[0.05, 0.1, 0.2]
3
Van der Pol oscillator
[5, 10, 25, 50, 100]
[1, 10, 20]
[0.001, 0.01, 0.1]
[0.05, 0.1, 0.2]
3
Rössler attractor
[5, 10, 25]
[1, 10, 20]
[0.001, 0.01, 0.1]
[0.05, 0.1, 0.2]
2
5-D Lorenz96 model
[5, 10, 15]
[1, 10, 20]
[0.001, 0.01, 0.1]
[0.05, 0.1, 0.2]
2
10-D Lorenz96 model
[4]
[1, 10, 20]
[0.001, 0.01, 0.1]
[0.05, 0.1, 0.2]
2
20-D Lorenz96 model
[2]
[1, 10, 20]
[0.001, 0.01, 0.1]
[0.05, 0.1, 0.2]
2
Appendix
Table 3: Hyperparameter sets of the Kramers-Moyal expansion and sparse regression algorithm. In every run, hyperparameter optimization is performed to maximize the accuracy of the recovery. The definitions of the hyperparameters are given in Sec. 3.3.3 .
Environment
p
g
m
s
p∗
g∗
η
Double well problem
500
50
5
15
100
15
0.1
Van der Pol oscillator
500
50
5
15
100
15
0.1
Rössler attractor
500
50
5
15
100
15
0.1
5-D Lorenz96 model
1000
100
5
20
200
15
0.1
10-D Lorenz96 model
2000
200
5
20
200
15
0.1
20-D Lorenz96 model
2000
200
5
20
200
15
0.1
Appendix
Table 4: Hyperparameters of the genetic programming algorithm for both identifying ordinary and stochastic differential equations. The hyperparameters are defined in Sec. 3.2 .
Institute of Structural Engineering, ETH Zurich, Zurich, Switzerland · ETH AI Center, ETH Zurich, Zurich, Switzerland · Swiss Data Science Center, Zurich, Switzerland
Data Science and Analytics Institute, University of Oklahoma, Norman, OK, USA · School of Industrial and Systems Engineering, University of Oklahoma, Norman, OK, USA