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 .
Understanding natural and engineered systems often relies on symbolic formulations, such as differential equations, which provide interpretability and transferability beyond black-box models. We introduce Latent Grammar Flow (LGF), a neuro-symbolic generative framework for discovering ordinary differential equations from data. LGF embeds equations as grammar-based representations into a discrete latent space and forces semantically similar equations to be positioned closer together with a behavioural loss. Then, a discrete flow model guides the sampling process to recursively generate candidate equations that best fit the observed data. Domain knowledge and constraints, such as stability, can be either embedded into the rules or used as conditional predictors.
Karin Yu, Eleni Chatzi, Georgios Kissas
Institute of Structural Engineering, ETH Zurich, Zurich, Switzerland · ETH AI Center, ETH Zurich, Zurich, Switzerland · Swiss Data Science Center, Zurich, Switzerland
Modeling real-world systems requires accounting for noise - whether it arises from unpredictable fluctuations in financial markets, irregular rhythms in biological systems, or environmental variability in ecosystems. While the behavior of such systems can often be described by stochastic differential equations, a central challenge is understanding how noise influences the inference of system parameters and dynamics from data. Traditional symbolic regression methods can uncover governing equations but typically ignore uncertainty. Conversely, Gaussian processes provide principled uncertainty quantification but offer little insight into the underlying dynamics. In this work, we bridge this gap with a hybrid symbolic regression-probabilistic machine learning framework that recovers the symbolic form of the governing equations while simultaneously inferring uncertainty in the system parameters. The framework combines deep symbolic regression with Gaussian process-based maximum likelihood estimation to separately model the deterministic dynamics and the noise structure, without requiring prior assumptions about their functional forms. We verify the approach on numerical benchmarks, including harmonic, Duffing, and van der Pol oscillators, and validate it on an experimental system of coupled biological oscillators exhibiting synchronization, where the algorithm successfully identifies both the symbolic and stochastic components. The framework is data-efficient, requiring as few as 100-1000 data points, and robust to noise - demonstrating its broad potential in domains where uncertainty is intrinsic and both the structure and variability of dynamical systems must be understood.
Discovering governing equations from observational data remains a fundamental challenge in scientific modeling, particularly when the underlying mathematical structure is unknown. Traditional sparse identification methods like SINDy excel at discovering parsimonious models but require researchers to specify candidate basis functions a priori, a limitation that often leads to model failure when critical terms are omitted or when systems exhibit unconventional dynamics. Purely symbolic regression approaches offer unlimited flexibility but struggle with noise sensitivity and frequently produce overly complex, unstable equations. We present AutoSINDy, a hybrid Discovery-then-Solve framework that combines the exploratory power of symbolic regression with the robust sparsity-promoting capabilities of SINDy. Our method operates in three stages: (1) PySR-based symbolic regression discovers candidate functional forms from bootstrapped data chunks; (2) a curation pipeline decomposes, expands, and filters these expressions using collinearity analysis to construct a minimal yet comprehensive library; and (3) SINDy identifies sparse governing equations from this custom-tailored library. Extensive experiments across canonical nonlinear systems demonstrate that AutoSINDy consistently recovers ground-truth equations even under high observational noise, achieving a ground-truth recovery rate of 92.8% across all trials. Compared with standard SINDy using enriched libraries and standalone symbolic regression, AutoSINDy achieves higher predictive accuracy, superior generalization to unseen trajectories, and substantially lower symbolic complexity.
Mohammad Amin Basiri, Charles Nicholson
Data Science and Analytics Institute, University of Oklahoma, Norman, OK, USA · School of Industrial and Systems Engineering, University of Oklahoma, Norman, OK, USA