cs.LGOct 26, 2026

Distilling neural networks into equations to predict planetary instability

Authors: Simon Alford, Edward Gu, Omar Costilla-Reyes, Dan Tamayo, Miles Cranmer, Kevin Ellis

Organizations: Department of Computer Science, Cornell University, Ithaca, NY 14850, USA · Computer Science and Artificial Intelligence Laboratory, MIT, Cambridge, MA 02139, USA · Department of Physics, Harvey Mudd College, Claremont, CA 91711, USA · Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge, UK · 

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

The problem of predicting when a multi-planet system will become dynamically unstable has defied analytical progress over the past three hundred years. Recent machine learning advances achieve high predictive performance of instability dynamics, but as black-box systems cannot be interpreted to connect to existing scientific understanding. Here we show analytic equations that predict when a multi-planet system will undergo instability. We discover these equations by first training a neural network on a dataset of N-body simulations, and then distilling its latent knowledge into closed-form algebraic expressions. Our two-stage pipeline first learns an interpretable variable basis, then searches for equations using symbolic regression. The resulting equations outperform existing analytic estimates of instability and apply across wide ranges of mass ratios, spacing, and orbital eccentricities. They capture structure in the prediction landscape with a small number of variables amenable to dynamical interpretation. These results demonstrate how symbolic learning can be combined with neural networks to extract scientific insight from complex dynamical systems.

Figures & tables

Figure 1: Predicting planetary instability. ( A ) Planetary systems evolve chaotically and can eventually undergo instability, such as a collision or ejection. ( B ) We train a two-stage neural network, consisting of a feature extractor and prediction network, on a dataset of NN -body integrations, and distill it into equations via two steps. ( C ) First, we distill a sparse linear feature space via regularization, pruning, and fine-tuning. ( D ) Second, we run symbolic regression to distill the predictor network into equations.
Figure 2: Discovered equations and their accuracy. ( A ) Discovered equations of higher complexity (horizontal axis) have lower error (vertical axis) on held-out data at/near resonance. The red outline shows the highest-accuracy equation. ( B ) The learned features used in the equations shown, comprising means and standard deviations of sparse linear combinations of the planets’ orbital elements. ( C ) RMSE (measured over unstable systems where Tinst<109T_{\rm inst}<10^{9} ) and stability classification accuracy (predicting whether Tinst>109T_{\rm inst}>10^{9} ) for discovered equations compared to neural network, pure symbolic regression, and existing analytic expressions when tested on systems in and near mean-motion resonances (Resonant), randomly initialized systems (Random), or on a grid of systems generated by systematically sweeping the period ratio of adjacent planets (Grid, see also Fig. 4 ).
Figure 3: Generalization to 5-planet systems. Instability time as a function of the interplanetary separation Δ\Delta for equal-mass, equally spaced, 5-planet systems. Distilled equations were learned from 3-planet systems before testing on five planets. Petit+20 includes a manual correction for 5-planet systems.
Figure 4: Visualizing model predictions. Model predictions, ground truth, and errors plotted as one varies the orbital period ratios of adjacent planets: P1/P2P_{1}/P_{2} (x-axis) and P2P_{2} / P3P_{3} (y-axis). Systems experiencing instability during the initial integration are colored white. ( A ) Instability time ground truth (via simulation) and predictions from distilled equations, neural network, and Petit+2020 ( Petit et al., 2020 ) . ( B ) RMSE of the distilled equations over the grid. ( C ) RMSE difference, neural network minus equation model. Regions where the equations are more accurate appear in blue, while regions where the neural network is more accurate appear in purple. ( D ) Predictions of four equations as one moves along the Pareto frontier of increasing equation complexity, illustrating how added complexity incorporates additional structure into the predictions (see corresponding labeled points on Pareto frontier in Fig. 2 A).
Appendix figures & tables14 assets

Supplementary material from the paper’s appendix.

Appendix

Table C1 : Table of evaluation metrics—RMSE, accuracy, false positive rate, false negative rate, bias, ROC-AUC, and log-likelihood (a measure of model quality combining RMSE and classification)—for all methods on the resonant, random, and period ratio grid datasets (“Pure SR 2” is Pure SR with no intermediate features). A positive classification corresponds to classifying a system as unstable.
Figure C1 : Predictions of models and baselines on the 5-planet dataset from Obertas et al. (2017) , also used in Fig. 3 .
Figure C2 : Predictions of Pure SR baselines compared to ground truth and distilled equations on the period ratio grid from Fig. 4 . The distilled equations better capture the dynamic range of instability times observed.
Figure C3 : Pareto frontier showing model RMSE on the resonant test set as a function of complexity when comparing architecture choices: either varying the number of terms kk in the sparse linear combination for features in the distilled ϕ\phi , implementing ψ\psi as a linear model with varying number of features, or using pure symbolic regression (with or without intermediate features) to generate equations. All models here use mean squared error (MSE) loss, rather than the combined classification and regression loss used for the main results.
Figure C4 : ROC curves for the distilled equations and baselines, along with ROC-AUC scores for each method. Stars highlight metrics when the natural classification threshold log⁡10Tinst=9\log_{10}T_{\rm inst}=9 is used, corresponding to metrics in Table C1 .
Figure D1 : Comparison of the predicted standard deviation of the neural network versus the distilled equations on the 5-planet dataset from Obertas et al. (2017) , also used in Fig. 3 . Using a neural network to predict the uncertainty of the equations improves prediction performance on the 5-planet dataset.
Figure D2 : Comparison of the predicted standard deviation of the neural network versus the distilled equations on the period ratio grid from Fig. 4 .
Figure D3 : Comparison of the predicted standard deviation of the neural network versus the distilled equations on the resonant and random datasets.
Table E1 : Full table of learned features distilled from neural ϕ\phi network using L1L_{1} regularization and pruning. Each feature is a linear combination of k=2k=2 input variables. Asterisked features are not used in any downstream equation.
Figure E1 : Heat map panels of feature activations on the period ratio grid for the first six learned means/standard deviations (remaining four are shown in Fig. E2 ). Colors in each panel span the 20th–80th-percentile range of that feature across the panel.
Figure E2 : Heat map panels of feature activations on the period ratio grid for the remaining four learned means/standard deviations (first six are shown in Fig. E1 ). Colors in each panel span the 20th–80th-percentile range of that feature across the panel.
Table E2 : Full table of the Pareto frontier yielded by running symbolic regression to distill neural ψ\psi into equations, including complexity and RMSE on the resonant test set. Frontier is selected based on validation loss, hence the lack of monotonicity of the test RMSE.
Figure E3 : First and second rows: density plots showing predicted versus true instability times for equations of increasing complexity and the neural network on the resonant test set. Orange lines show average predictions, indicating overall bias of each model as a function of true log instability time. A bias is observed for both the neural network and the distilled equations, but the bias is larger for the equations. Third row: Mean residual of the complexity 26 equation on the resonant test set as a function of the value of additive subterms of the equation or predicted instability time (right). The compression of predicted instability time, especially for low predicted instability times, can be traced to the subterm in the leftmost plot.
Figure E4 : Heat map of model predictions over a sweep of resonant angle and middle planet orbital period, demonstrating equation sensitivity to changes in instability time near the separatrices surrounding MMRs.
More in this topic

Dynamical Systems

Dynamical systems studies how states evolve over time, including stability, feedback, chaos, and transitions in physical, biological, and engineered systems.

See all 482 papers

Explore similar work

CardsList
Aug 22, 2026

Symbolic Neural ODEs: Learning interpretable models from time-series data

We present a machine learning framework for identifying sparse, interpretable models of dynamical systems directly from time-series data. Our approach parameterizes the underlying vector field using a neural architecture and trains it by minimizing a multi-step prediction loss over a finite horizon. To ensure numerical tractability, we optimize a mean absolute error objective averaged across prediction steps, and progressively increase the horizon during training. A key feature of this formulation is that it enforces consistency under repeated composition of the learned dynamics. As a result, the identified models exhibit significantly improved stability compared with approaches based on one-step regression of the vector field. When combined with sparsity-promoting regularization, this leads to parsimonious models that generalize beyond the training data. We demonstrate accurate recovery of systems exhibiting a wide range of behaviors, including stable and unstable fixed points, periodic orbits, and chaotic attractors. For chaotic systems, while long-term trajectory prediction is inherently limited by sensitivity to initial conditions, we show that multi-step training yields models with accurate short-term dynamics and strong agreement in long-time statistical properties, including mean, variance, and Lyapunov exponents. Moreover, we establish theoretical bounds linking trajectory error to statistical accuracy, providing a step toward a principled explanation for this behavior.

Nibodh Boddupalli, Jeff Moehlis

Dept. of Mechanical Engineering, University of California, Santa Barbara, CA 93106

Nonlinear System IdentificationSystem IdentificationNeural ODEsSparse Identification of Nonlinear DynamicsOrdinary Differential EquationsChaotic Time-Series Prediction
AskPDF
Apr 1, 2026

Predicting Dynamics of Ultra-Large Complex Systems by Inferring Governing Equations

Predicting the behavior of ultra-large complex systems, from climate to biological and technological networks, is a central unsolved challenge. Existing approaches face a fundamental trade-off: equation discovery methods provide interpretability but fail to scale, while neural networks scale but operate as black boxes and often lose reliability over long times. Here, we introduce the Sparse Identification Graph Neural Network, a framework that overcome this divide by allowing to infer the governing equations of large networked systems from data. By defining symbolic discovery as edge-level information, SIGN decouples the scalability of sparse identification from network size, enabling efficient equation discovery even in large systems. SIGN allows to study networks with over 100,000 nodes while remaining robust to noise, sparse sampling, and missing data. Across diverse benchmark systems, including coupled chaotic oscillators, neural dynamics, and epidemic spreading, it recovers governing equations with high precision and sustains accurate long-term predictions. Applied to a data set of time series of temperature measurements in 71,987 sea surface positions, SIGN identifies a compact predictive network model and captures large-scale sea surface temperature conditions up to two years in advance. By enabling equation discovery at previously inaccessible scales, SIGN opens a path toward interpretable and reliable prediction of real-world complex systems.

Qi Shao, Duxin Chen, Jiawen Chen +5

School of Mathematics, Southeast University, Nanjing, 210096, P. R. China. · School of Electronic Engineering and Computer Science, Queen Mary University of London, London, E1 4NS, UK. · School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, UK. +4

Dynamical SystemsNonlinear System IdentificationSparse Identification of Nonlinear DynamicsSpatiotemporal ForecastingSpatiotemporal GNNs
AskPDF
Aug 16, 2025

Universal Learning of Nonlinear Dynamics

We study the fundamental problem of one-step prediction of a marginally stable unknown nonlinear dynamical system. We describe an algorithm for this problem, based on the technique of spectral filtering, which learns a mapping from past observations to the next based on a spectral representation of the system. Using techniques from online convex optimization, we prove vanishing prediction error for any nonexpansive nonlinear dynamical system with finitely many marginally stable modes, with rates governed by a novel quantitative control-theoretic notion of learnability. The main technical component of our method is a new spectral filtering algorithm for linear dynamical systems, which incorporates past observations and applies to general noisy and marginally stable systems. This generalizes the original spectral filtering algorithm to both asymmetric dynamics as well as incorporating noise correction, and is of independent interest.

Evan Dogariu, Anand Brahmbhatt, Elad Hazan

Courant Institute of Mathematical Sciences, New York University · Princeton University

Dynamical SystemsNonlinear System IdentificationSpectral Methods
AskPDF