Equation discovery began as search over model structures. Automated reverse engineering evolved both model forms and parameters to reproduce measured time series [2], and symbolic regression by genetic programming famously distilled conserved quantities and free-form laws from physical data [3]. Sparse identification of nonlinear dynamics inverted the paradigm: with a fixed library of candidate terms, the governing equation becomes a sparse coefficient vector recoverable by thresholded least squares [4]. The software ecosystem around this idea is mature [14], extensions reach partial differential equations and noisy derivative estimation [9], and relaxed regularizers such as SR3 unify robust and structural regression with improved convergence behaviour [12]. Coordinate discovery adds a further layer, learning observability transformations under which the sparse representation becomes compact [13], and low-data-limit variants target control applications [15]. AI Feynman demonstrated that physics-inspired constraints, such as symmetries and separability, can make the symbolic search space tractable [10], while modern evolutionary engines such as PySR scale the search with typed operators and parsimony pressure [11].

Library discovery depends on a coordinate system and derivative estimates in which the target law can be represented. Coordinate-learning and noise-aware variants address some of these issues [9], [13]. A decomposition stage is a design option for slow-fast data, but full-state multiscale ODEs need not be nonsparse by construction. The chain used here already has a sparse full-state polynomial law.

Neural models of dynamics learn flows rather than equations. Neural ordinary differential equations parameterize the vector field and integrate it with differentiable solvers [5]. Universal differential equations embed neural components inside known mechanistic structure, enabling hybrid models whose unknown parts are learned [16]. Physics-informed networks constrain the learner with residual losses from known partial differential equations [17], a programme whose scope and open problems are charted in [18]. Energy-based structure offers another constraint: Hamiltonian networks conserve a learned energy by construction [19], and Lagrangian networks extend this to general coordinates [20]. Linear-operator alternatives such as dynamic mode decomposition trade expressivity for spectral interpretability [21].

Neural models can absorb fast transients into latent states and predict observations without identifying a unique symbolic slow law. Hybrid models recover some structure [16], but symbolic postprocessing is still a separate proposed step. The HERMES architecture assigns representation learning, decomposition, and refinement to such possible stages; the experiment here tests only a smoothing–screening–reweighting regression core on observed full-state trajectories.

Attracting slow manifolds require more than scale separation: normally hyperbolic transverse directions must be stable, with suitable compactness and regularity conditions [6], [7]. Regularized noisy-data differentiation [22], robust losses [23], stability selection [24], and information criteria [25], [26] address different statistical issues. Their guarantees depend on assumptions and tuning; citing them does not establish that the proposed composition inherits all of them.

Multiscale model reduction has its own data-driven tradition, from equation-free computation to manifold learning, but the outputs are numerically defined objects rather than closed-form equations. The present design accepts that limitation for its neural stages and refuses it at the output: whatever the intermediate representations, the deliverable is a symbolic equation whose complexity can be counted and compared. This stance forces the model-selection question to the front. A discovery method that cannot say how complex its answer is, or how that complexity trades against predictive fidelity, is not yet a scientific instrument, and the trade-off is therefore formalized as a first-class object in Section 4.

Table 2 summarizes selected architectures. HERMES’s row describes proposed stages, not all executed components. Sparse regression can use noise-aware derivatives, coordinate discovery can learn representations, and symbolic models can incorporate physical constraints. The present question is whether the specific regression implementation improves the displayed outcomes under the stated preprocessing and budget.

HERMES (this losses, stability selection slow-fast Pareto front work) decomposition

Note: The HERMES row describes the proposed full design; only a regression-side realization is evaluated.

The design is organized around three research questions. They concern, respectively, whether a neural decomposition can prepare slow coordinates well enough for identification, whether hybridizing library regression with seeded evolutionary search improves the discovery front, and how robustness degrades under noise and weak separation. Each question is paired with a falsifiable hypothesis.

RQ1 asks whether a learned slow-fast decomposition supports identification of reduced terms under noisy partial observations.

RQ2 asks whether seeded symbolic refinement improves the complexity–accuracy front relative to its parent methods.

RQ3 asks how recovery changes with noise, separation, and record length.

These are design-level questions. The current campaign sweeps observation noise for three fixed full-state systems; it does not execute neural decomposition, evolutionary refinement, or systematic separation/record-length sweeps.

H1 predicts useful term recovery after a specified learned decomposition. It remains untested without that stage and independently verified coordinates.

H2 predicts a benefit of seeded refinement at matched search vocabulary and budget. It remains untested. H3 predicts a recovery boundary under noise and separation sweeps. The current chain failure is a limitation of one implementation, not an estimated general boundary.

We consider slow-fast systems in the standard form of Eq. (1), where x collects the slow variables, z the fast variables, eps the ratio of timescales, and mu a parameter. Only a subset of the state is observed, through the measurement model of Eq. (2), with eta_k denoting additive noise whose level is summarized by a signal-to-noise ratio. The observed record therefore mixes slow trend, fast transients, and measurement noise in unknown proportions. The unknown objects of interest are the slow manifold Z and the reduced drift F, and Table 1 collects the notation used throughout.

ẋ = f(x, z; μ), εż = g(x, z; μ), 0 < ε ≪ 1. (1)

y = h(x(t ), z(t )) + η , k = 1, … , K. (2)

k k k k

x, z Slow and fast state variables Latent decomposition eps Timescale-separation ratio System parameter mu System parameter Sweep axis y_k Observation at time t_k Data h Observation operator Partial observability SNR Signal-to-noise ratio of the record Stress factor U Slow coordinates recovered in Stage A Regression input Z, Z_theta True and learned slow manifold Invariance object F Reduced drift on the manifold Discovery target Theta(U) Library of candidate terms Stage B basis xi Sparse coefficient vector Stage B output c(F) Symbolic node count of F Complexity measure H, L_roll Rollout horizon and loss Fidelity measure

For a smooth slow-fast system, the critical manifold at ε=0 satisfies the first relation of Eq. (3). Under appropriate compact normal hyperbolicity, a nearby locally invariant manifold can persist for small positive ε [6], [7]; attraction requires stable transverse fast dynamics. The finite-ε invariance relation contains DZεf, so setting g=0 alone is a quasi-steady approximation, not exact invariance. Eq. (4) defines the leading-order reduced drift. Recovering this object from noisy partial observations is the proposed discovery problem; the current experiments instead identify full-state ODEs.

g(x, Z0(x; μ); μ) = 0, g(x, Zε(x; μ); μ) = εDxZε(x; μ)f(x, Zε(x; μ); μ).

ẋ = F (x; μ) = f(x, Z (x; μ); μ). (4) A future implementation could report candidate laws and their validation-error/complexity front. Structural agreement with a known law belongs to evaluation, not training selection. Support error and coefficient error are separate: an error metric restricted to true support can ignore additional false terms, which must be assessed through support precision and exact matching.

Stage A is proposed to learn slow and fast coordinates from delayed observations, reconstruct measured channels, and predict trajectories. A learned manifold head Zθ can be regularized by a finite-ε invariance residual, rather than requiring fast drift to vanish exactly. Eq. (5) states that residual for compatible learned drift models. These coordinates, networks, and regularization weights must be specified and validated before any manifold-recovery claim. Robust loss options [23] do not resolve identifiability by themselves.

P = ∑∥g (u , Z (u )) − εDZ (u )f (u , Z (u ))∥ . (5) man ∥ θ i θ i θ i θ i θ i ∥

Proposed slow/fast heads, gating, and ensemble disagreement are modeling choices. Delay embedding [8] can supply an alternative state representation under its assumptions, but not uniquely named physical coordinates or sparse reduced laws. The neural output and its uncertainty remain unexecuted design components in this study.

The proposed sparse stage uses a declared library and derivative estimates. For a Huber loss ρδ, Eq. (6) sums losses over residual components, with optional nonnegative weights; it does not square the Huber loss again. Sequential thresholding is an algorithmic approximation to sparse selection and is not guaranteed to solve the stated nonconvex objective globally. Stability selection [24] and BIC-style scoring [26] require their own assumptions, which may fail for autocorrelated trajectories. The implemented reweighted STLSQ core is a different realization; Huber optimization and stability-selection guarantees are not assigned to it.

min ∑ w ρ (U̇ − [Θ(U)ξ] ) + λ ∥ ξ ∥ , w ≥ 0. (6)

ij δ ij ij 0 ij

Noise may affect support and coefficient estimates through derivative error and library conditioning. The proposed uncertainty-weighted regression requires validated weights; it is not tested here. The executed core uses correlation screening and reweighted thresholding on smoothed full-state data.

Stage C refines the shortlist by evolutionary symbolic search in the style of PySR [11]. Each candidate support seeds the initial population, so the search begins from a structurally informed region rather than from random expressions, while standard mutation and crossover operators explore replacements, simplifications, and eliminations. Complexity is counted as the number of symbolic nodes, which makes parsimony explicit and comparable across candidates. The search optimizes the objective of Eq. (7), a rollout loss plus a complexity penalty, evaluated by integrating each candidate equation and comparing it with withheld slow-trajectory segments. Multi-objective bookkeeping retains the full Pareto set rather than a single winner.

min ℒroll(F) + λcc(F). (7)

Seeding symbolic search may improve efficiency, but neither a better Pareto front nor domination of unseeded search follows from initialization alone. A matched-budget ablation must test it. The current component ablations concern regression only.

A future Pareto output should retain candidates nondominated in validation rollout error and complexity. Scalarization by Eq. (7) samples only part of a potentially nonconvex front. A knee selection needs an explicit normalization and distance/tie rule rather than an informal corner description. No such front is produced in the executed study.

This framing has two consequences. It forbids the common practice of reporting the single best-recovered equation while the method actually produces many, and it renders the method honest about degeneracy: multiple simple laws may explain finite data, and their coexistence on the front is a finding about the data rather than a defect of the algorithm. The executed campaign accordingly scores the single delivered model per realization against the known true support and coefficients. In this respect the pipeline is closer to a measurement instrument for structure than to a point estimator.

Algorithm 1 describes the proposed full architecture. Exact networks, losses, parameter ranges, supports, budgets, integrators, and selection rules remain necessary for reimplementation. Section 5 specifies a separate regression-only realization.

Algorithm 1. HERMES - hybrid equation recovery for noisy multiscale data.

Input: observations {(t_k, y_k)}, library Theta, operator set O, thresholds schedule T, generations G, complexity weight lam_c 1: Stage A: train slow-fast networks with manifold penalty of Eq. (5); ensemble across seeds to obtain U, Z_theta, uncertainty weights 2: estimate dU/dt on smoothed slow trajectories (regularized method) 3: Stage B: for each threshold in schedule T do 4: solve thresholded robust regression of Eq. (6) on subsamples 5: retain terms by stability selection; score supports by BIC 6: end for; output ranked candidate supports S_1, S_2, ...

7: Stage C: seed evolutionary search with each support 8: for generation = 1 to G do 9: mutate / crossover symbolic candidates over operator set O 10: evaluate L_roll on withheld segments; update Pareto front 11: end for 12: apply an explicitly specified normalized selection rule; report the full front Output: Pareto front of symbolic laws with the knee selection

Computational cost includes neural training and integration, library construction, repeated least-squares solves, screening, bagging or stability subsamples, and candidate rollouts. For K rows and p library terms, a dense least-squares solve commonly requires O(Kp²+p³) work; thresholds and iterations multiply this cost. Evolutionary cost depends on population size, expression complexity, integration horizon, and solver steps. Memory must include trajectories, the K×p library, models, and candidate populations. The small regression runtime does not measure full-architecture cost.

The following are assumptions of the proposed architecture, not conditions all checked by the full-state campaign:

I. A1 requires a suitable compact normally hyperbolic critical manifold; attraction additionally requires stable fast transverse directions [6].

II. A2 requires adequate observations for a valid embedding or identifiable coordinate model; an embedding alone does not recover unique slow physical coordinates [8].

III. A3 requires the target law to be representable in the library or reachable operator vocabulary.

IV. A4 specifies an approximately additive stationary noise model for the intended record.

V. A5 requires data covering the time constants and excitation needed to identify the target terms.

Complementarily, the pipeline is expected to fail in identifiable situations, declared here as scope rather than observed. Anticipating these failures disciplines the claims made on the pipeline's behalf. The stress grid treats each condition as a cell to be probed deliberately.

Without useful separation, slow-fast reduction may be inappropriate; that does not imply the full dynamics lack an identifiable model. An omitted functional form cannot be recovered by the fixed library. Derivative error and ill-conditioning may prevent support recovery under some noise levels. Short or insufficiently excited records can leave multiple candidate laws indistinguishable.

The campaign contains models labeled Lotka–Volterra [27] and Duffing [28] with two states each, and a three-state linear chain. These historical works identify the model families; they do not validate the specific modern equations or the present simulation. The stated libraries, sample counts, and initial conditions are two-degree/501 samples at Δt=0.02 for Lotka–Volterra, cubic/626 samples at Δt=0.02 for Duffing, and degree-two/851 samples at Δt=0.1 for the chain. Chain rates are 2.0, 0.2, and 0.02, corresponding to time constants 0.5, 5, and 50. Gaussian noise levels are 0,1,5,10% of per-state standard deviation. Twenty realizations per system and level give 240 system-realization records overall and 60 at each level. The underlying state-equation count is (2+2+3)×20=140 per level, not 60. The reported table aggregation must therefore be clarified from the original records before modelwise and equationwise interpretations are assigned.

The methods share observation records but differ in preprocessing. Ridge regression [29] uses the full library; STLSQ [4] uses fixed thresholds 0.2,0.02,0.05 for the respective systems. The ensemble-lite variant bags 25 fits on 80% of time rows with a 60% support-retention rule; it is not the full Ensemble-SINDy method [30]. HERMES uses Savitzky–Golay filtering/differentiation [31], a window bank {5,9,15,25,41,61}, correlation screening, and reweighted STLSQ with adaptive thresholds and κ=0.5. Its described generalized-cross-validation selection adapts a criterion developed for linear smoothing/regularization [32]. Baselines use raw second-order finite differences. Matched-smoothed baselines are needed to isolate benefits of sparse selection from derivative preprocessing.

The reported metrics include support precision/recall, norm-based coefficient relative error on the true support, and a 10% maximum per-term coefficient tolerance for exact recovery. False extra coefficients are excluded from the true-support coefficient metric and must be assessed separately. Whether exactness requires all equations in a model or is averaged equationwise remains unresolved by the source summaries. The 5% level is the designated endpoint. The stated master seed 20261011, 36.3 s runtime, and software configuration, including SciPy [33] and matplotlib [34], require code and logs. Neural decomposition, hidden-state recovery, evolutionary refinement, and Pareto output are not executed; no release identifier or dated preregistration is established.

At 5% noise Table 3 reports pooled exact rate 0.667 for HERMES and 0.000 for all three baselines. The source interpretation is 40 successes among 60 realization records, requiring confirmation of modelwise aggregation. Coefficient errors are 0.0551±0.0647 for HERMES, 0.4615±0.5710 for STLSQ, 0.4601±0.5696 for ensemble-lite, and 0.8926±0.1073 for ridge. These summaries demonstrate a difference between the tested implementations, but their pooling units and raw records are needed for precise inference. Reported p-values are retained descriptively pending verification of test method and assumptions.

At zero noise the plotted STLSQ curve is 0%, whereas the source prose claims a pooled STLSQ rate of 0.667 and success on both single-scale systems. The contradiction and aggregation unit require the original per-system indicators and plotting code; there are sixty system-realization records and 140 underlying state equations per level (Fig. 1).

Note: Pooled summaries and per-system exact rates are preserved; modelwise versus equationwise aggregation requires original records.

The per-system summaries show exact rate 1.000 for HERMES on Lotka–Volterra and Duffing at 5% noise, and 0.000 on the chain. At 10% Duffing’s rate 0.900 means 18 of twenty successes and two failures if the unit is a realization, not one failure. The chain’s reported precision 1.000 and recall 0.833 are consistent with omitting a weak true term. The claimed third-equation omission of coefficient −0.02 and selected threshold 2.0 require the delivered support records; aggregate precision/recall alone do not identify them.

HERMES has the lowest displayed error on Lotka–Volterra and Duffing. Its error rises on the multiscale chain, where the displayed STLSQ and ensemble-lite errors are approximately order one. The plotted comparisons are descriptive; the underlying records are needed to resolve the pooling unit in Table 3.

Table 4 compares the implemented core with three removals. Without smoothing, the 5% exact rate is 0.633 versus 0.667, and the 10% rate is 0.350 versus 0.633. Removing screening gives 0.000 at 5%; removing reweighting gives 0.033, with precision 0.677 and recall 0.944. The lower precision directly indicates false extra terms, so the loss of exactness cannot be attributed solely to coefficient accuracy. These are scoped component contrasts, not proof that each stage is universally necessary or that neural/evolutionary hybridization works.

Note: Pooling unit requires confirmation from the sixty system-realization records at each level.

At zero noise, the source prose claims STLSQ and ensemble-lite each recovered Lotka–Volterra and Duffing exactly, corresponding to pooled 0.667 if the unit is a system realization. Fig. 1 instead plots STLSQ at 0% and ensemble-lite near 66.7%. STLSQ noiseless exact success is therefore unestablished until per-system exact indicators, the aggregation rule, and plotting code reconcile text and figure. The source reports no exact chain recovery. HERMES’s Duffing rate of 0.900 at 10% implies two failures among twenty realizations if that is the unit. Reported chain smoothing windows include 15 samples for the fast state and 61 for the slow state; the bank also contains 5 and 9, so 15 is not its shortest window. Recorded windows and validation errors are needed to assess their selection.

The noiseless single-scale coefficient errors favor STLSQ: 0.0005 versus 0.0018 on Lotka–Volterra, and 0.0006 versus 0.0142 on Duffing. This is consistent with a smoothing-related trade-off for these records. Chain failure at zero noise does not prove that denoising, alternative thresholds, longer excitation, or another estimator cannot recover it. Conditioning and support-selection effects should be examined before declaring a structural impossibility.

The source reports Wilcoxon comparisons [35], Holm adjustments [36], and Friedman summaries [37]. The claim that all p-values are exact is inconsistent with their scale: forty nonzero same-sign binary differences give a two-sided sign-randomization value 2/2⁴⁰≈1.819e−12, rather than the reported 2.540e−10. With sixty nonzero same-sign continuous differences, W=0 has exact value 2/2⁶⁰≈1.735e−18, rather than 1.630e−11. The reported values resemble normal approximations, but the actual zero/tie handling and software method require code. Tables retain source values as descriptive, unverified summaries. Shared system records, equations, trajectories, and noise realizations must be handled at the appropriate independent unit before confirmatory inference.

The tested smoothing–screening–reweighting combination has the highest reported pooled exact rate at 5% noise, alongside lower true-support coefficient error. The baselines have different derivative preprocessing, which limits method-level attribution. Removal contrasts show differences among these variants. The multiscale chain remains unrecovered, and reported omitted-term details require raw support traces.

The core has larger noiseless single-scale coefficient errors than STLSQ and lower chain recall at 5% (0.833 versus 0.992), while selecting fewer false terms. H1 about learned slow-coordinate recovery and H2 about seeded refinement remain untested. Three fixed systems and a noise sweep cannot estimate a general separation or record-length boundary for H3.

The empirical scope is full-state polynomial ODE identification for three fixed systems with twenty realizations per cell. The hidden-state, neural-decomposition, symbolic-refinement, and Pareto claims remain proposed architecture components. Chain failure informs limitations of the tested realization without falsifying an unexecuted decomposition.

The reported proof-of-concept outcomes require the original trajectory generation, noise realizations, delivered supports, coefficient vectors, and analysis code for verification.

Slow-manifold approximation errors depend on separation, stability, regularity, domain, and observation horizon [6], [7]. Partial observations need not identify physical coordinates. Several candidate laws can fit finite records, and library or operator restrictions define what can be recovered. The proposed Pareto representation may expose ambiguity but does not eliminate it. None of these full-architecture issues is resolved by the present full-state regression experiment.

The proposed decomposition cannot remove closure error merely by penalizing residuals: unresolved fast modes may influence a slow approximation, and proposed robust losses or uncertainty weights have not been executed or validated here. Partial observability adds a further layer, since delay embedding under its assumptions can recover a representation of attractor dynamics, not unique named physical coordinates or an identified sparse law [8]. Real records also contain drift, unmodeled inputs, and instrument artifacts; the current full-state simulated campaign does not establish transfer to them.

Stage outputs can make failures easier to inspect only if they are actually generated and evaluated. The current ablations do not test neural-stage attribution or symbolic-refinement fronts. Matched preprocessing and independent rollout validation would make the regression comparison more informative.

Longer or cleaner records may help when remeasurement is feasible. For archived records, the proposed hybrid pipeline would still need an explicit approximation horizon and an independently held-out rollout test before it could be treated as a reduced model. Fast variables need not disappear after a transient, and a finite-horizon validation target cannot guarantee global asymptotic fidelity.