Discovery of Nonlinear Dynamics with Automated Basis Function Generation

Abstract: 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.

• The paper introduces AutoSINDy, a three-stage pipeline that mines candidate functions with symbolic regression, removes redundant terms, and estimates sparse coefficients with bootstrap-validated SINDy.
• AutoSINDy recovered the ground-truth equations in 92.8% of 180 noisy benchmark trials, including 90% recovery at 5% noise, outperforming standard SINDy and standalone PySR.
• The method achieved 68.3% excellent long-horizon simulations with zero hard integration crashes, but its scalability, performance on real data, and evaluation for chaotic systems remain open challenges.

Motivation and problem statement

Sparse identification of nonlinear dynamics (SINDy) recovers governing ODEs by sparse regression over a user-specified library of candidate functions, but its success is bounded by the quality of that library: if a critical basis function is omitted, the true model cannot be recovered regardless of the downstream optimizer, while enriching the library indiscriminately induces multicollinearity that destabilizes coefficient estimation. Symbolic regression (SR), conversely, searches the full expression space without a library but is noise-sensitive, computationally intensive, and prone to overfitted, unstable expressions, and because it fits each state variable independently it cannot enforce coupled structure across a multi-dimensional system. The paper identifies the construction of a minimal yet expressive function library directly from data, without domain knowledge, as the central open problem in equation discovery, and proposes AutoSINDy to address it.

The AutoSINDy framework

AutoSINDy is a three-stage "Discovery-then-Solve" pipeline. The key structural insight is that symbolic regression need not solve the complete discovery problem; it need only surface plausible candidate functional forms, which are then curated and handed to SINDy for precise sparse coefficient identification.

Stage 1 — Discover. PySR is applied to $K$ short, randomly sampled data chunks rather than the full trajectory, which improves tractability, exposes different dynamical regimes, and prevents over-representation of any attractor region. For each chunk and each target derivative $\dot{x}_j$, the entire Pareto-optimal front of expressions is harvested, not just the best model, since structurally correct atomic sub-terms often appear in sub-optimal Pareto expressions. Nested trigonometric compositions are forbidden to keep the search comparable with traditional methods.

Stage 2 — Curate. Raw expressions are parsed with SymPy and decomposed into additive atoms with numeric prefactors stripped. Configurable expansion strategies (severe, gentle, hybrid) determine whether grouped polynomial and trigonometric sub-expressions are expanded into monomials. Candidates are then sorted in ascending SymPy operator count and admitted greedily via pairwise Pearson correlation pruning (or, alternatively, variance inflation factor filtering), so a complex term is retained only if no simpler term already explains the same variance. This complexity-ordered traversal enforces a simplicity bias directly at the library level, independent of the downstream optimizer.

Stage 3 — Identify. The curated library is fitted with STLSQ or SR3 inside a bootstrap ensemble following the E-SINDy paradigm. A hard-cutoff mask retains a term only if its empirical inclusion probability across $B$ bootstrap replicates exceeds $\kappa = 0.80$, with the final coefficient being the ensemble mean over replicates where the term was active. Both separate per-variable libraries (the default, producing lower complexity) and a unified pooled library are supported.

The design provides two complementary layers of noise rejection: structural redundancy is removed before regression, and statistically unreliable coefficients are removed after it.

Experimental design

The framework is evaluated on six canonical systems spanning linear, limit-cycle, and chaotic regimes, from the harmonic oscillator (true canonical complexity 3) to a Complex Lorenz system with a mixed polynomial–trigonometric coupling term $x_1\sin(x_0+x_2)$ (complexity 15). Each of the $6 \times 6 \times 5 = 180$ trials per method (540 total) uses a 5000-snapshot trajectory with additive Gaussian noise at relative amplitudes $\sigma \in \{0, 0.01, \dots, 0.05\}$. Two evaluation choices deserve emphasis: all generalization metrics are computed on a completely unseen, noise-free validation trajectory from a different initial condition, and both instantaneous derivative prediction and long-horizon simulation are reported, since a model can fit local derivatives well yet diverge catastrophically under integration. Equation recovery is a strict binary criterion ($R^2 > 0.99$ on the validation trajectory), and complexity is measured as the SymPy operator count of the fully expanded equation — a metric the authors argue is more meaningful than the standard $\ell_0$ term count, which treats $x_1$ and a fourth-degree monomial as equivalent.

Results

The headline result is a 92.8% ground-truth recovery rate across all 180 trials, versus 58.9% Excellent-tier derivative prediction for Standard SINDy with an enriched polynomial–Fourier library and 28.3% for standalone PySR. In long-horizon simulation the gap widens: AutoSINDy achieves Excellent simulation in 68.3% of trials against 18.9% (SINDy) and 8.3% (PySR), and — a strong claim supported across all systems and noise levels — produces zero hard integration crashes in all 180 trials, while Standard SINDy diverges in 51.7% and Standard PySR in 75.0% of trials.

Noise robustness is equally pronounced. Even at zero noise, Standard SINDy recovers only 83% and PySR 67%, indicating failure modes that are structural rather than purely noise-driven. At the highest noise level ($\dot{x}_j$0), AutoSINDy maintains 90% recovery versus 40% and 20% for the baselines.

Per-system results highlight the mechanism of failure in each baseline:

System	AutoSINDy recovery	Standard SINDy	Standard PySR	SINDy median complexity (truth)
Harmonic oscillator	100%	40%	16.7%	68.5 (3)
Damped pendulum	96.7%	46.7%	33.3%	58.5 (4)
Van der Pol	100%	80%	0%	176 (6)
Duffing	70%	0%	16.7%	11 (6)
Complex Lorenz	90%	93.3%	3.3%	35.5 (15)

Standard SINDy's Van der Pol median complexity of 176 operators against a ground truth of 6 is the most striking instance of library bloat; the enriched basis admits hundreds of correlated alternatives among which the optimizer distributes coefficient mass under noise. Standard PySR's 0% recovery on Van der Pol despite a median derivative $\dot{x}_j$1 of 0.8 illustrates the complementary failure: it identifies $\dot{x}_j$2 but cannot estimate its coefficient precisely enough to remain stable under integration. On Complex Lorenz, Standard SINDy's Fourier library genuinely cannot represent the compound coupling term, so its nominally high 93% recovery rests on a 39-operator polynomial–Fourier approximation — a case where the recovery proxy and true structural identification diverge, a caveat worth noting. AutoSINDy recovers a median complexity of 16 there, within one operator of the ground truth.

AutoSINDy's discovery overhead is modest (median ~24 s, comparable to PySR since both invoke the same engine, and 10–30× Standard SINDy), but it is consistently the fastest method at simulation time (median ~1.2 s), because compact right-hand sides produce less stiff ODE systems. Standard SINDy's simulation times show outliers exceeding $\dot{x}_j$3 s on two-dimensional systems.

Discussion

The paper's analysis of why the hybrid outperforms both components is well supported. Standalone PySR lacks regularization, so slightly incorrect coefficients accumulate into divergence; standalone SINDy suffers both over-completeness under noise and genuine structural incompleteness for terms outside its fixed basis. AutoSINDy resolves both by mining system-specific functional forms and pruning to a non-redundant basis before any regression occurs. The zero hard-crash rate is interpreted as evidence that discovered equations remain within a physically plausible manifold even when coefficients carry error — a substantive claim, though it rests on six synthetic systems with Gaussian noise.

The Duffing oscillator marks a clear boundary of the method: its bistable potential concentrates trajectory density near equilibria where the $\dot{x}_j$4 derivative is most noise-sensitive, producing a bimodal outcome in which 70% of trials recover the exact dynamics and the remainder fail with incorrect coefficient signs, with no intermediate regime.

Limitations and open questions

The authors are explicit about several limitations. Discovery time scales linearly in the number of state variables and chunks, which may be prohibitive for high-dimensional systems; proposed mitigations (amortized symbolic regression, parallelization, warm-starting) remain untested. On chaotic systems, correct structural recovery does not guarantee accurate long-horizon simulation — Complex Lorenz simulation $\dot{x}_j$5 is negative for all methods — and the paper notes that evaluation should move to attractor-geometry metrics insensitive to phase drift. All experiments use synthetic ODE data with controlled Gaussian noise; extension to PDEs and to real experimental time series with irregular sampling, partial observability, and unknown ground truth is left as an open validation question. It is also worth noting that the recovery-rate proxy can mislead when a fixed library happens to approximate a missing compound term well, as in the Complex Lorenz comparison, and that the reported advantage depends on hyperparameter choices ($\dot{x}_j$6, $\dot{x}_j$7, chunk size) whose sensitivity is not systematically ablated.

Conclusion

AutoSINDy operationalizes the complementarity between symbolic regression and sparse identification: PySR serves as an unsupervised basis-mining engine, a collinearity-aware curation pipeline with a simplicity bias converts its output into a well-conditioned library, and ensemble SINDy performs robust sparse coefficient recovery. Across 180 noisy trials on six benchmark systems, the method achieves a 92.8% recovery rate, near-ground-truth symbolic complexity, zero hard simulation crashes, and the fastest simulation times, without any prior specification of the system's mathematical form. The principal open questions are scalability to high-dimensional systems, applicability to PDEs and real experimental data, and evaluation of chaotic systems with trajectory-independent metrics.