OrthoReg: Orthogonal Regularization for Hybrid Symbolic-Neural Dynamical Systems
arXiv:2606.19145v1 [cs.LG] 17 Jun 2026
Till Richter Technical University of Munich Helmholtz Munich [email protected]
Niki Kilbertus Technical University of Munich Helmholtz Munich [email protected]
Abstract Dynamical systems are fundamental to modeling the natural world, yet modeling them involves a persistent trade-off: manually prescribed mechanistic models are interpretable by design but often overly simplistic and misspecified; in contrast, flexible data-driven neural methods lack physical insight. Hybrid modeling aims for the best of both worlds by combining a prescribed or symbolic, physics-based component with a flexible neural network. A critical challenge, however, is that the neural component may relearn mechanistic parts, yielding redundant and uninterpretable models, especially when the symbolic structure itself is discovered from data. Existing methods based on standard L2 regularization rely on a projection argument that breaks when the symbolic component is learned through sparse discovery, allowing the neural augmentation to overlap with symbolic structure. We introduce OrthoReg (Orthogonal Regularization), which directly penalizes overlap between the symbolic and neural components, preventing symbolic structure from being absorbed by the neural residual. This yields a complementary decomposition: the symbolic part captures what the library can express, and the neural part captures what remains. On benchmark dynamical systems with partial library mismatch, OrthoReg improves symbolic recovery and out-of-distribution behavior.
1
Introduction
Dynamical systems modeling has long been a cornerstone across the sciences, especially for the natural and life sciences. Applications range from healthcare data [Choi et al., 2016, Hess et al., 2024, Seedat et al., 2022, Richter et al., 2026], climate modeling [Rolnick et al., 2022, Eyring et al., 2024], and power systems [Toubeau et al., 2018], just to name a few. However, it faces a fundamental trade-off: symbolic models, traditionally specified by hand, provide interpretability by design, but typically cannot capture complex unknown phenomena; flexible neural networks instead excel at fitting data from dynamical systems [Chen et al., 2018] but lack physical insight. mechanistic fully prescribed, interpretable by design, no data required, usually simple, assumed to be the correct mechanistic system description
hybrid
data-driven fully data driven, complex black-box model, approximate system behavior within data support
Hybrid modeling becomes essential when domain knowledge is partial: epidemiological models may capture core transmission dynamics but miss behavioral feedback mechanisms that modulate contact rates; mechanical systems follow known laws of motion but exhibit complex friction and damping effects not easily expressible in closed form; climate models encode fundamental physics but require data-driven corrections for sub-grid processes. In these scenarios, a purely symbolic approach underperforms due to missing phenomena, while purely neural models sacrifice the interpretability and physical consistency that domain experts require for scientific insight and decision-making. Hybrid modeling approaches [Rackauckas et al., 2020, Yin et al., 2021, Zou et al., 2024] combine Preprint.
Figure 1: Symbolic–neural decompositions under library mismatch. The symbolic and augmented function spaces can overlap. Pure symbolic models are restricted to the symbolic span; standard hybrids may redundantly explain the shared region with both fphy and faug ; OrthoReg discourages this redundancy by pushing faug away from Fphy , yielding a complementary decomposition. physical priors (predetermined symbolic or parametrized expressions) with learned neural corrections expected to capture phenomena that are unknown or too complex to model directly. They promise the best of both worlds, but still require substantial prior knowledge in crafting the mechanistic part. This work focuses on sparse library-based discovery of the symbolic component, as in SINDy-style methods [Brunton et al., 2016]. In this regime, the symbolic component is fitted with a continuous sparsity penalty (e.g. L1 regularisation), where coefficient shrinkage leaves in-library residuals – the practically relevant case our analysis targets. Our goal is to make the symbolic and neural components complementary under the empirical inner product. In their landmark paper, Yin et al. [2021] present the APHYNITY framework for hybrid dynamical systems modeling when the symbolic structure (but not exact parameter values) is known a priori. Like a large body of related existing work [Rackauckas et al., 2020, Mouli et al., 2024], APHYNITY decomposes the (autonomous) vector field of an ordinary differential equation (ODE) as f = fphy + faug , where fphy ∈ Fphy = span{ϕj }M j=1 captures dynamics within a predetermined library of “symbolic” functions {ϕj }M (e.g., polynomials, trigonometric functions), while faug is supposed j=1 to capture the residual dynamics via flexible neural networks. When the symbolic structure is fixed and Fphy is a closed linear subspace, the two components are provably separated via simple L2 regularization of faug (Theorem A.1). In this setting, optimizing out faug reduces to an L2 projection onto the fixed subspace Fphy , so the residual (and hence faug ) is orthogonal to Fphy , which underlies the analysis of Yin et al. [2021]. APHYNITY gives a projection-based view of hybrid modeling when the physical model class is fixed: the augmentation has a minimal norm subject to fitting the observed dynamics, and in finitedimensional L2 subspaces the residual is orthogonal to the physical component. This fixed-class picture changes when the symbolic component is learned through gradient-based sparse discovery. Continuous sparsity penalties, such as L1 -regularized library learning, shrink coefficients away from the unregularized projection and can leave residual structure inside the symbolic span. A neural augmentation can then relearn symbolic directions even when small in norm. OrthoReg targets this regime by directly penalizing empirical overlap between the augmentation and the symbolic library. Example 1.1 (Counterexample: sparsity forces overlap). Let F be a function space with inner product ⟨·, ·⟩D and induced norm ∥ · ∥D . Let ϕ1 , ϕ2 ∈ F be orthonormal w.r.t. ⟨·, ·⟩D and f = ϕ1 + ϕ2 . Consider min ∥f − (w1 ϕ1 + w2 ϕ2 + g)∥2D + λ2 ∥g∥2D + µ∥w∥1 . 2 w∈R , g∈F
2λ2 2) If 0 < µ < 1+λ , the unique minimizer satisfies w1 = w2 = 1 − µ(1+λ and g = 2λµ2 (ϕ1 + ϕ2 ), 2λ2 2 hence g ̸⊥ span{ϕ1 , ϕ2 } although span{ϕ1 , ϕ2 } contains f .
Thus, hybrid dynamical models that include symbolic discovery require explicit mechanisms to prevent neural-symbolic overlap. Figure 1 illustrates this failure mode and how OrthoReg restores a clean symbolic–neural split. In this work, we introduce a theoretically grounded and practically effective method to learn hybrid dynamical systems, where the mechanistic component is discovered from data via symbolic regression while ensuring that the residual neural component remains orthogonal to the symbolic part. The key contribution is to bridge the gap between theoretical projection theory arguments and practical optimization to leverage both symbolic and data-driven modeling. Concretely, we provide: 2
• Theoretical analysis of symbolic–neural overlap under sparse library-based discovery: L2 regularization can induce overlap between the two components (Theorem A.2), whereas OrthoReg yields a clean in-library / orthogonal-complement split under the empirical inner product (Theorem B.2). • An algorithmic solution realizing the orthogonality constraint within the joint sparse-discovery objective. • Empirical validation on benchmark dynamical systems with partial library mismatch, showing improved out-of-distribution generalization and sparse symbolic recovery. 1
2
Related Work
We review data-driven methods for learning governing dynamics and the interpretability– expressiveness trade-off motivating OrthoReg. (Dynamic) symbolic regression. Symbolic regression recovers closed-form expressions from data using genetic programming [Koza, 1994, Schmidt and Lipson, 2009], neural-guided search [Petersen et al., 2021, Udrescu and Tegmark, 2020], or sparse library-based regression such as SINDy [Brunton et al., 2016]. Sparse-library methods are highly interpretable when the true dynamics lie in the candidate set, but their expressiveness is limited by library design. Recent work improves symbolic discovery through physical constraints such as unit consistency [Tenachi et al., 2023], transformer-based symbolic generation [Lample and Charton, 2020, Biggio et al., 2021, Valipour et al., 2021, Kamienny et al., 2022, Vastl et al., 2024], extensions to ODEs and trajectory data [Becker et al., 2023, d’Ascoli et al., 2024, Sun et al., 2023], and methods for noisy, sparse, distributional, high-dimensional, or guided discovery settings [Qian et al., 2022, Dakhmouche et al., 2025, Liang and Yang, 2025, Tian et al., 2025, Hu et al., 2025]. However, when relevant effects are not representable by a compact symbolic library, purely symbolic models either fail or lose interpretability, motivating hybrid decompositions with a symbolic core and a flexible residual. Physics-informed neural networks. PINNs [Raissi et al., 2019] embed differential equations as soft constraints and Universal ODEs [Rackauckas et al., 2020] parameterize unknown vector-field components with neural networks; surveys [Cuomo et al., 2022, Hao et al., 2022] position these as central scientific-ML paradigms. Both rely on a fixed symbolic part and do not address overlap when the symbolic component is itself discovered under sparsity. Neural, symbolic, and hybrid methods. Hybrid methods combine symbolic interpretability with neural expressiveness. Gray-box discovery has been studied for PDEs via sparse regression [Rudy et al., 2017] and physics-informed symbolic objectives [Chen et al., 2021, Kiyani et al., 2023]. For ODEs, APHYNITY formalizes f = fphy + faug but assumes fixed symbolic structure [Yin et al., 2021]; related work alternates neural learning with symbolic distillation [Grigorian et al., 2024] or studies PINN/PIKAN optimization for gray-box identification [Daryakenari et al., 2026]. Pure neural ODE models can encode sparsity [Aliee et al., 2022], conservation/manifold constraints [Greydanus et al., 2019, Matsubara and Yaguchi, 2022, White et al., 2023], or meta-learned physics–ML trade-offs [Mouli et al., 2024], but generally lack explicit symbolic recovery. In contrast, OrthoReg targets the overlap problem that arises when the symbolic component is itself discovered jointly with a neural residual.
3
Background and Problem Setting
Let F be a class of measurable vector fields f : Rn → Rn with finite L2 (ν) norm, where ν is the data-generating measure on states. We use the population inner product ⟨f, g⟩ν = Ex∼ν [f (x)⊤ g(x)] PN and the empirical semi-inner product ⟨f, g⟩D = N1 i=1 f (xi )⊤ g(xi ) on the sampled states D = {xi }N i=1 , with induced norm and seminorm ∥ · ∥ν and ∥ · ∥D . For a closed linear subspace Fphy ⊆ F, we denote the corresponding population projection by PFν phy ; empirically, PFDphy denotes the projection under ∥ · ∥D , understood up to equivalence on D. 1 Code: https://github.com/richtertill/OrthoReg.
3
The functions f ∈ F are interpreted as vector fields of autonomous, first-order differential equations dx = f (x), dt
with solution trajectories x : R → Rn .
Following prior work [Yin et al., 2021, Rackauckas et al., 2020], we assume an additive decomposition f = fphy + faug ,
fphy ∈ Fphy , faug ∈ F
of vector fields of interest into a “physical” (or symbolic/mechanistic) component and an “augmented” (or neural/residual) component. The space Fphy ⊆ F of candidate symbolic components is typically restricted to functions that can be represented in closed form using known functions, so as to be amenable to direct interpretation and dissemination by humans. Most existing methods assume fphy to be either known exactly, or to be given as a parametric family in which only a (usually small) set of parameters is unknown. Practically, this is often implemented via a linear combination of non-linear basis functions: M X Fphy = αj ϕj αj ∈ R for j ∈ {1, . . . , M } , (1) j=1
where the dictionary functions ϕj : Rn → Rn are fixed. The dynamics governing most real systems are not perfectly described by such simple closed-form expressions, but contain higher-order effects or complex interactions that are rarely captured by simple interpretable mathematical expressions. To capture such residual effects, the augmentation faug ∈ F is supposed to be flexible and expressive, albeit potentially not easily interpretable. Hence, a natural choice to represent faug is via flexible function approximators such as neural networks, giving rise to the term “neural component.” Crucially, the neural component should only capture effects that cannot be captured by the symbolic component. In the current formulation, one could simply set faug ≡ f and fphy ≡ 0. However, this would undermine the entire idea of hybrid modeling. When the physical model class is fixed, Yin et al. [2021] provide thorough theoretical guarantees showing that a relatively simple norm-based regularization scheme is sufficient to ensure that faug “only captures what is necessary, but not more.” A simplified vector-field version of this norm-regularized principle is min
fphy ∈Fphy ,faug ∈F
∥f − fphy − faug ∥2D + λ2 ∥faug ∥2D .
(2)
For fixed fphy , writing r = f − fphy , the minimizer over faug is fˆaug =
1 1 r= (f − fphy ). 1 + λ2 1 + λ2
Substituting back gives min
faug ∈F
∥f − fphy − faug ∥2D + λ2 ∥faug ∥2D =
λ2 ∥f − fphy ∥2D . 1 + λ2
Thus, for λ2 > 0, optimizing over fphy is equivalent to empirical least-squares projection onto Fphy . If Fphy is a closed linear subspace, the Hilbert-space projection theorem gives PFDphy (f ) up to D-equivalence, and the residual f − P D (f ), hence fˆaug , is orthogonal to Fphy under ⟨·, ·⟩D Fphy
[Lax, 2014]. APHYNITY extends this projection view beyond linear subspaces via existence and uniqueness guarantees under proximinality and Chebyshevness [Yin et al., 2021]; see Section A. A natural extension beyond a fully known fphy , or the structure being known up to a small set of parameters, is to discover the symbolic component itself from a larger dictionary, as in SINDy [Brunton et al., 2016]. After fixing the candidate basis functions {ϕj }M j=1 we select only a small support set S ⊂ {1, . . . , M } of basis functions that enter the expression with non-zero coefficients. The induced function space is Fphy (S) := span{ϕj | j ∈ S}. In practice, the set S is fitted via sparse regression methods (e.g., L1 regularization ∥ · ∥1 or more involved iterated sparse regressions as in SINDy) to encourage small supports S. 4
While at first this appears to be a natural extension to APHYNITY, on closer inspection this changes the projection geometry underlying its argument. When sparse discovery is implemented through continuous shrinkage penalties, the symbolic fit is no longer an unregularised projection onto the selected span; the residual can therefore retain components in the symbolic library, which the neural augmentation may relearn (cf. Figure 1 and Theorem 1.1). This is the fundamental gap our work addresses: expressive data-driven sparse symbolic discovery requires additional mechanisms to keep symbolic and neural components from overlapping. The closest faithful adaptation of APHYNITY to this setting is therefore a sparse symbolic component combined with an L2 -regularized neural augmentation, which we call the “Hybrid L2 ” baseline.
4
Method: OrthoReg for Hybrid Modeling
We observe state trajectories of an unknown dynamical system and aim to learn an autonomous vector field f : Rn → Rn such that ẋ = f (x). Following the hybrid modeling setup in Section 3, we model fˆ(x) = fˆphy (x; w) + fˆaug (x; ϑ), fˆphy (x; w) =
M X
wj ϕj (x),
(3)
j=1
ˆ where {ϕj }M j=1 is a fixed symbolic library, w is sparse, and faug (·; ϑ) is a flexible neural augmentation. Learning from trajectories. The data consists of state observations over time (trajectories). Depending on the setting, derivative targets may or may not be available. We therefore consider two standard data-fit losses for learning fˆ. If (approximate) derivatives yi ≈ ẋ(ti ) are available (e.g., from a simulator or numerical differentiation), we can fit fˆ by vector-field regression N
Lvf (w, ϑ) =
1 X 2 yi − fˆ(xi ) . N i=1
(4)
If only states are observed, we fit by one-step prediction. Let T be the set of observed one-step transitions (xt , xt+1 , ∆tt ) (potentially with irregular ∆tt ) and let Ψ∆t (·; fˆ) denote a one-step ODE solver (e.g., explicit Euler Ψ∆t (x; fˆ) = x + ∆t fˆ(x) or a Runge–Kutta method). The one-step loss is Lstep (w, ϑ) =
1 |T |
2 xt+1 − Ψ∆tt (xt ; fˆ) .
X
(5)
(xt ,xt+1 ,∆tt )∈T
In the following, Lfit denotes either Lvf or Lstep . We train with Lvf (state-only training is also compatible with OrthoReg; see Section I), as derivative-space training was empirically more stable, and evaluate in both derivative and state space. OrthoReg uses the empirical inner product ⟨·, ·⟩D over observed states D = {xi }N i=1 , defined in Section 3. 4.1
OrthoReg Objective
OrthoReg encourages a non-redundant decomposition by penalizing correlation between the neural augmentation and the symbolic library. Concretely, we penalize the (squared) empirical inner products
L⊥ reg (ϑ) = λ
M D E2 X fˆaug (·; ϑ), ϕj , D
j=1
(6)
ˆ so that L⊥ reg = 0 implies faug ⊥ span{ϕ1 , . . . , ϕM } w.r.t. ⟨·, ·⟩D . Let µ ≥ 0 control symbolic sparsity and λ ≥ 0 the orthogonality strength, then the full OrthoReg training objective is min Lfit (w, ϑ) + µ∥w∥1 + L⊥ reg (ϑ), w,ϑ
5
(7)
Algorithm 1: OrthoReg training (vector-field regression). M Input: samples {(xi , yi )}N i=1 , library {ϕj }j=1 , weights λ, µ ≥ 0. Initialize w ∈ RM , neural parameters ϑ. For each epoch and minibatch B ⊂ {1, . . . , N } of size B: PM 1. fˆ(xi ) = j=1 wj ϕj (xi ) + fˆaug (xi ; ϑ). P 2. Lfit = B1 i∈B ∥yi − fˆ(xi )∥2 . P PM 2 3. gj = 1 g . fˆaug (xi ; ϑ)⊤ ϕj (xi ), Lorth = λ B
j=1 j
i∈B
4. Lsparse = µ∥w∥1 . 5. Update (ϑ, w) by autodiff on Lfit + Lorth + Lsparse . The guarantee in Theorem 4.1 is tied to the additive objective in Equation (7): the orthogonality penalty must act on the additive residual fˆaug at the same empirical inputs as the data-fit. Staged or alternating updates of the same objective fall within the same formal scope when they reach a solution satisfying the theorem assumptions, whereas unconstrained residual fitting after shrinkage-based symbolic learning and compositional architectures of the form Tw + Nϑ ◦ Tw do not (Theorem C.1 and Section C). 4.2
Implementation and Computational Considerations
Implementation. Algorithm 1 summarises training. We optimise the continuous L1 -regularised objective with Adam and report sparsity after thresholding small coefficients. The orthogonality penalty costs O(B n M ) per batch; in mini-batch training, finite-batch squaring of the inner-product estimate introduces a small bias, so our experiments use full- or large-batch estimates unless stated otherwise. Hyperparameter selection, feature scaling, and the autonomous-augmentation caveat for non-autonomous data are described in Sections F and G.8. 4.3
Theoretical Guarantees
We first show that L2 regularization does not prevent – and can in fact induce – overlap between symbolic and neural components under sparse selection (Theorem A.2, Section A). The OrthoReg objective restores control via two structural results: a penalty bound showing that empirical overlap is controlled at the optimum and vanishes as λ → ∞, and a triangular error decomposition (the in-library / orthogonal-complement split) that, under the empirical inner product, serves as a formal lens for the symbolic–neural decomposition. Theorem 4.1 (Penalty bound for empirical orthogonality). Assume λ > 0 and that Equation (7) admits a global minimizer (w∗ , ϑ∗ ). Let wsym ∈ arg minw Lfit (w, ϑ0 ) + µ∥w∥1 denote the best pure-symbolic fit at fˆaug ≡ 0 (parameters ϑ0 ). Then M X ⟨fˆaug (·; ϑ∗ ), ϕj ⟩2D ≤ λ1 Lfit (wsym , ϑ0 ) + µ∥wsym ∥1 , j=1
so empirical orthogonality is enforced as λ → ∞. Proof idea. The orthogonality term is a quadratic penalty: comparing the optimum to the baseline fˆaug ≡ 0 yields a 1/λ bound on constraint violation. A full proof is given in Theorem B.1. Error decomposition. Beyond the penalty bound, orthogonality has a concrete geometric consequence for approximation. In the idealized vector-field regression setting, and when the learned augmentation is empirically orthogonal to Fphy , the squared error admits an orthogonal decomposition: ∥f − (fˆphy + fˆaug )∥2D = ∥PFDphy (f ) − fˆphy ∥2D + ∥f − PFDphy (f ) − fˆaug ∥2D , (8) Thus, the symbolic error is the gap to the in-library projection PFDphy (f ), while the neural error is the gap to the empirical orthogonal-complement component f − PFDphy (f ); see Theorem B.2. Under the empirical inner product, this decomposition serves as a formal lens for the symbolic–neural split: each component’s error is the gap to its half of the orthogonal direct sum. 6
Comparison to L2 regularization. L2 regularization controls only the magnitude of the augmentation through the decomposition ∥fˆaug ∥2D = ∥PFDphy (fˆaug )∥2D + ∥fˆaug − PFDphy (fˆaug )∥2D , where the equality follows from orthogonal decomposition in inner product spaces [Rudin, 1987]. Even when the total is small, the component PFDphy (fˆaug ) can be non-zero, allowing neural–symbolic overlap. When Fphy is learned through sparsity constraints, this overlap can occur even in realizable settings, motivating explicit orthogonality regularization (Section A). Empirical-to-population transfer. OrthoReg controls overlap under the empirical inner product used for training. A formal population guarantee for the learned augmentation would require uniform convergence over the augmentation hypothesis class; we leave this to future work (Section I).
5
Experiments
Table 1: Damped Pendulum results. Scale factors ek inside parentheses are powers of ten; x/ẋ are state/derivative spaces. Metric
Pure
L2
OrthoReg
In-Distribution Performance MSE(ẋID , e−2) (↓) MSE(xID , e−2) (↓) MSE(ẋID,ext , e2) (↓)
1.27 ±0.07 0.74 ±0.43 19.45 ±12.89
1.41 ±0.05 0.31 ±0.02 5.00 ±2.42
1.89 ±0.01 3.65 ±0.36 0.13 ±0.07
Out-of-Distribution Performance MSE(xOOD,T2 , e0) (↓) MSE(ẋOOD,T3 , e2) (↓)
1.02 ±0.07 34.00 ±22.22
1.07 ±0.02 7.90 ±4.02
1.01 ±0.01 0.02 ±0.00
System Identification Quality F1 (↑) #Terms (↓)
0.43 ±0.09 11.4 ±2.9
0.61 ±0.03 6.8 ±0.4
0.93 ±0.15 3.6 ±1.3
We evaluate OrthoReg across four dynamical systems of different complexities: a modified damped pendulum, a Lotka–Volterra predator-prey system, a time-modulated SIR epidemiological model, and a Duffing oscillator. Our evaluation focuses on two complementary metrics: (i) trajectory accuracy measured by normalized mean-squared error (MSE) on derivatives and integrated states2 on in-distribution, out-of-distribution, and extrapolation states; and (ii) symbolic recovery quality measured by F1 score and the number of symbolic terms. We compare three approaches: pure symbolic regression (SINDy), an L2 -regularized hybrid, and an OrthoReg-regularized hybrid. We report along four axes consistent with multi-axis evaluation principles in scientific ML [Wyder et al., 2026]: initial-condition extrapolation, parametric extrapolation, regime/basin extrapolation, and robustness under irregular sampling and observation noise. 5.1
Damped Pendulum: Missing Dynamics
The modified damped pendulum extends the classical damped pendulum [Kharkongor and Mahato, 2018] with three forcing terms that lie outside the symbolic dictionary: 2
θ̈ + αθ̇ + ω02 sin θ = β1 cos(3θ) + β2 e−θ + β3 tanh(θ̇),
(9)
with (β1 , β2 , β3 ) = (0.3, 0.25, 0.15). The dictionary used by the symbolic component is PolynomialLibrary(degree = 2) ⊕ FourierLibrary(nfreq = 1), which contains sin θ, cos θ, sin(2θ), cos(2θ) and monomials of (θ, θ̇) up to total degree two; the three forcing terms above lie outside this dictionary by construction. Table 1 reports performance under medium missing dynamics (β1 = 0.3, β2 = 0.25, β3 = 0.15 in Equation (9)). Derivative and state MSE quantify trajectory fit, and F1 measures symbolic recovery against ground-truth terms. OrthoReg yields the strongest symbolic recovery (highest F1, fewest 2 MSE values are normalized by the squared norm of the target signal for scale invariance.
7
terms) and the largest reductions in OOD-derivative error, at the cost of higher in-distribution state MSE. Quantitatively, OOD-derivative error drops by more than two orders of magnitude relative to the L2 baseline while in-distribution state MSE is higher than the L2 baseline. This trade-off is consistent with the symbolic / orthogonal-residual split established by Theorem B.2: enforcing the split sacrifices some in-distribution fit in exchange for components that transfer to OOD evaluations. 5.2
Comparison with Pure Neural Baselines
Table 2: Baseline comparison on the modified damped pendulum. OrthoReg uniquely provides symbolic recovery while remaining competitive on prediction. Metric MSE(ẋID , e−2) (↓) MSE(xOOD,T2 , e0) (↓) MSE(ẋOOD,T3 , e2) (↓) F1 (↑)
PINN
Universal ODE
OrthoReg
1.58 ±0.00 1.02 ±0.01 0.94 ±0.00
1.61 ±0.06 1.07 ±0.07 1.02 ±0.18
1.89 ±0.01 1.01 ±0.01 0.02 ±0.00
–
–
0.93 ±0.15
We also compare to pure neural approaches in Table 2, including PINNs [Raissi et al., 2019] and Universal Differential Equations [Rackauckas et al., 2020]. While these baselines achieve competitive trajectory fitting, they cannot recover symbolic components. OrthoReg recovers a sparse symbolic component and is competitive on OODT 2 and OODT 3 derivative MSE. Implementation details are in Section G.6. Canonical PySINDy STLSQ (same library and splits as the main experiments) is summarised in Section G.6: it improves symbolic support recovery relative to the Adam L1 puresymbolic rows when thresholding is the sparsity mechanism, but unstable rollouts in partial mismatch show why residual capacity remains necessary.
Figure 2: Ablations. (a) OrthoReg is most effective under partial library mismatch. (b) Irregular sampling degrades all methods, while OrthoReg retains competitive OOD behavior. (c) Intermediate λ gives the best trade-off between symbolic recovery and OOD error; right panel shows a residual– symbolic cosine diagnostic for OrthoReg to visualize separation. 5.3
Ablations
We ablate three factors that probe when and why orthogonal regularization is beneficial: (i) the degree to which the true dynamics exceed the symbolic library, (ii) robustness to non-ideal (irregular) observation times, and (iii) the strength of the orthogonality penalty. Figure 2 summarizes the main trends. Complete quantitative results and additional ablation analyses are reported in Section E. 8
Figure 3: (a) Damped pendulum in (θ, ω), with ω := θ̇. (b) Duffing in (x, ẋ) with shaded basins. Solid: inferred; dashed: ground truth. OrthoReg captures global trends; pure symbolic and L2 distort the dynamics. Rightmost panels show the OrthoReg decomposition, highlighting the neural residual. Dataset difficulty (Figure 2a). OrthoReg’s gains on stress-test generalization (xOOD,T3 ) are largest in the intermediate regime where the system partially exceeds the library; in low-missing regimes (β1 = 0.1, β2 = 0.08, β3 = 0.05) differences shrink on simpler evaluations, and in high-missing regimes (β1 = β2 = β3 = 2.0) all methods degrade with OrthoReg retaining a modest advantage. Full results in Section E.1. Sampling scheme (Figure 2b). Under irregular sampling (∆t ∼ Uniform[0.01, 0.05] vs. regular ∆t = 0.02), absolute performance drops across all methods, but OrthoReg retains a marginal F1 advantage and a competitive xOOD,T3 . Full breakdown in Section E.2. Regularization strength (Figure 2c). OrthoReg has an intermediate operating range: weak λ under-enforces the orthogonality penalty, while overly strong λ constrains the augmentation. We select λ and µ by validation prediction loss on held-out in-distribution trajectories; F1 is tracked as a diagnostic only (Section G.8). Right: cosine diagnostic visualizing residual–symbolic alignment within each method; regularization strength is assessed by F1 and OOD error. 5.4
Duffing oscillator: regime-dependent generalization
Table 3: Duffing oscillator. Training on the positive basin; OOD T2 tests cross-basin (Regime-Ext), OOD T3 tests parameter shifts (Param-Ext). Metric
Pure
L2
OrthoReg
Predictive Performance MSE(xID , e−1) (↓) MSE(xOOD,T2 , e0) (↓) MSE(xOOD,T3 , e−1) (↓)
4.75 ±0.13 11.03 ±6.09 3.70 ±0.21
5.36 ±0.14 7.46 ±6.20 3.44 ±0.12
5.59 ±0.09 4.93 ±0.35 3.46 ±0.05
System Identification Quality F1 (↑) #Terms (↓)
0.40 ±0.05 11.0 ±1.4
0.59 ±0.04 3.8 ±0.4
0.63 ±0.05 3.4 ±0.5
We study the unforced Duffing oscillator, a multistable system whose cross-basin behavior probes whether models capture global structure beyond the observed basin [Göring et al., 2024], and train on trajectories from the positive basin only. Here orthogonality is helpful because the symbolic library carries the global polynomial structure of the dynamics; discouraging the neural residual from absorbing those terms lets the symbolic component transfer across basins. Table 3 and Figure 3 show that in-distribution all methods are comparable, while on cross-basin OOD T2 OrthoReg achieves the lowest error and strongest symbolic recovery; under parameter shifts (OOD T3) differences are within noise (details in Section F.2). The gain is largest where evaluation stresses global state-space structure rather than raw fit, consistent with the orthogonal split established by Theorem B.2. 5.5
Cross-System Validation
As auxiliary stress tests beyond the pendulum and Duffing, we evaluate OrthoReg on a Lotka– Volterra predator-prey system (coupled polynomial dynamics) and a time-modulated SIR model 9
(severe library mismatch and state-dependent timescales); full results are in Section F. We read LV as a low-headroom sanity check (Table 6: methods are comparable, with a small F1 gain for OrthoReg) and SIR as a misspecification stress test (Table 7: OrthoReg trades ID and OOD-T2 accuracy for the strongest OOD-T3 extrapolation and the sparsest symbolic recovery – highest F1 and fewest terms). The pattern is consistent with the symbolic–neural split. 5.6
Summary
Across systems, OrthoReg is most effective in the regime it is designed for: partial library misspecification, where the symbolic library captures part of the dynamics but leaves structured residual effects. In this setting, OrthoReg improves symbolic recovery and reduces out-of-distribution derivative error, while accepting some loss in in-distribution fit. The gains persist under irregular sampling, but diminish under heavy observation noise and when the missing dynamics are strongly entangled with library features (Section E.4). These results support the view of OrthoReg as a decomposition regularizer: it is not a universal accuracy booster, but a mechanism for preserving a meaningful symbolic component while assigning non-symbolic effects to the residual.
6
Conclusion
Hybrid symbolic–neural models promise interpretable structure together with flexible residual modeling, but this promise depends on a clean division of labor between the two components. We showed that this division becomes fragile when the symbolic component is discovered through sparse library learning: standard L2 regularization controls the magnitude of the neural residual, but not whether it relearns symbolic structure. We introduced OrthoReg, which makes complementarity explicit by penalizing empirical overlap between the neural augmentation and the symbolic library. Our analysis shows that projection-based regularization can induce symbolic–neural overlap under sparse discovery (Theorem A.2), whereas OrthoReg controls empirical orthogonality violation at global minimizers (Theorem 4.1) and yields an in-library / orthogonal-complement error split (Theorem B.2). Empirically, OrthoReg improves symbolic recovery and out-of-distribution behavior in regimes of partial library mismatch, where the library captures part of the dynamics but leaves structured residual effects. Several limitations remain. Our guarantees are empirical-inner-product statements for additive models with fixed differentiable libraries and continuous sparsity relaxations; derivative supervision, feature scaling, and noise sensitivity also matter in practice (Section I). A key next step is to extend orthogonality-regularized residual learning to broader symbolic-discovery pipelines, including canonical SINDy selection–refit, genetic-programming, and transformer-based symbolic regression. More broadly, OrthoReg points toward hybrid scientific models in which approximate mechanistic vocabularies and flexible residuals can coexist without erasing the symbolic explanation.
Acknowledgements We thank Manuel Lubetzki and Alessandro Palma for helpful discussions and constructive feedback on this manuscript.
10
References Hananeh Aliee, Till Richter, Mikhail Solonin, Ignacio Ibarra, Fabian Theis, and Niki Kilbertus. Sparsity in continuous-depth neural networks. Advances in Neural Information Processing Systems, 35:901–914, 2022. Sören Becker, Michal Klein, Alexander Neitz, Giambattista Parascandolo, and Niki Kilbertus. Predicting ordinary differential equations with transformers. In International conference on machine learning, pages 1978–2002. PMLR, 2023. Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neural symbolic regression that scales. In International conference on machine learning, pages 936–945. Pmlr, 2021. Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3932–3937, 2016. Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. Zhao Chen, Yang Liu, and Hao Sun. Physics-informed learning of governing equations from scarce data. Nature communications, 12(1):6136, 2021. Edward Choi, Mohammad Taha Bahadori, Jimeng Sun, Joshua Kulas, Andy Schuetz, and Walter Stewart. Retain: An interpretable predictive model for healthcare using reverse time attention mechanism. Advances in neural information processing systems, 29, 2016. Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maziar Raissi, and Francesco Piccialli. Scientific machine learning through physics–informed neural networks: Where we are and what’s next. Journal of Scientific Computing, 92(3):88, 2022. Ramzi Dakhmouche, Ivan Lunati, and Hossein Gorji. Robust symbolic regression for dynamical system identification. Transactions on Machine Learning Research, 2025. Nazanin Ahmadi Daryakenari, Khemraj Shukla, and George Em Karniadakis. Representation meets optimization: Training pinns and pikans for gray-box discovery in systems pharmacology. Computers in Biology and Medicine, 201:111393, 2026. Stéphane d’Ascoli, Sören Becker, Philippe Schwaller, Alexander Mathis, and Niki Kilbertus. ODEFormer: Symbolic regression of dynamical systems with transformers. In International Conference on Learning Representations, volume 2024, pages 21943–21976, 2024. Veronika Eyring, William D Collins, Pierre Gentine, Elizabeth A Barnes, Marcelo Barreiro, Tom Beucler, Marc Bocquet, Christopher S Bretherton, Hannah M Christensen, Katherine Dagon, et al. Pushing the frontiers in climate modelling and analysis with machine learning. Nature Climate Change, 14(9):916–928, 2024. Niclas Alexander Göring, Florian Hess, Manuel Brenner, Zahra Monfared, and Daniel Durstewitz. Out-of-domain generalization in dynamical systems reconstruction. In Forty-first International Conference on Machine Learning, 2024. Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. Advances in neural information processing systems, 32, 2019. Gevik Grigorian, Sandip V George, Sam Lishak, Rebecca J Shipley, and Simon Arridge. A hybrid neural ordinary differential equation model of the cardiovascular system. Journal of the Royal Society Interface, 21(212), 2024. Zhongkai Hao, Songming Liu, Yichi Zhang, Chengyang Ying, Yao Feng, Hang Su, and Jun Zhu. Physics-informed machine learning: A survey on problems, methods and applications. arXiv preprint arXiv:2211.08064, 2022. 11
Konstantin Hess, Valentyn Melnychuk, Dennis Frauen, and Stefan Feuerriegel. Bayesian neural controlled differential equations for treatment effect estimation. In International Conference on Learning Representations, volume 2024, pages 52586–52622, 2024. Herbert W Hethcote. The mathematics of infectious diseases. SIAM review, 42(4):599–653, 2000. Jiao Hu, Jiaxu Cui, and Bo Yang. Learning interpretable network dynamics via universal neural symbolic regression. Nature Communications, 16(1):6226, 2025. Pierre-Alexandre Kamienny, Stéphane d’Ascoli, Guillaume Lample, and François Charton. End-toend symbolic regression with transformers. Advances in Neural Information Processing Systems, 35:10269–10281, 2022. D Kharkongor and Mangal C Mahato. Resonance oscillation of a damped driven simple pendulum. European Journal of Physics, 39(6):065002, 2018. Elham Kiyani, Khemraj Shukla, George Em Karniadakis, and Mikko Karttunen. A framework based on symbolic regression coupled with extended physics-informed neural networks for graybox learning of equations of motion from data. Computer Methods in Applied Mechanics and Engineering, 415:116258, 2023. John R Koza. Genetic programming as a means for programming computers by natural selection. Statistics and computing, 4(2):87–112, 1994. Adam J Kucharski, Timothy W Russell, Charlie Diamond, Yang Liu, John Edmunds, Sebastian Funk, Rosalind M Eggo, Fiona Sun, Mark Jit, James D Munday, et al. Early dynamics of transmission and control of covid-19: a mathematical modelling study. The lancet infectious diseases, 20(5): 553–558, 2020. Guillaume Lample and François Charton. Deep learning for symbolic mathematics. In International Conference on Learning Representations, 2020. Peter D Lax. Functional analysis. John Wiley & Sons, 2014. Senwei Liang and Haizhao Yang. Finite expression method for solving high-dimensional partial differential equations. Journal of Machine Learning Research, 26(138):1–31, 2025. Takashi Matsubara and Takaharu Yaguchi. FINDE: Neural differential equations for finding and preserving invariant quantities, 2022. S Chandra Mouli, Muhammad Alam, and Bruno Ribeiro. Metaphysica: Improving OOD robustness in physics-informed machine learning. In The Twelfth International Conference on Learning Representations, 2024. Brenden K Petersen, Mikel Landajuela Larma, Terrell N. Mundhenk, Claudio Prata Santiago, Soo Kyung Kim, and Joanne Taery Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. In International Conference on Learning Representations, 2021. Zhaozhi Qian, Krzysztof Kacprzyk, and Mihaela van der Schaar. D-code: Discovering closed-form odes from observed trajectories. In International Conference on Learning Representations, 2022. Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scientific machine learning. arXiv preprint arXiv:2001.04385, 2020. Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686–707, 2019. Till Richter, Weixu Wang, Alessandro Palma, and Fabian J Theis. Generative models of cell dynamics: from neural odes to flow matching. Communications Biology, 2026. 12
David Rolnick, Priya L Donti, Lynn H Kaack, Kelly Kochanski, Alexandre Lacoste, Kris Sankaran, Andrew Slavin Ross, Nikola Milojevic-Dupont, Natasha Jaques, Anna Waldman-Brown, et al. Tackling climate change with machine learning. ACM Computing Surveys (CSUR), 55(2):1–96, 2022. Walter Rudin. Real and complex analysis. McGraw-Hill, Inc., 1987. Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Data-driven discovery of partial differential equations. Science advances, 3(4):e1602614, 2017. Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 324(5923):81–85, 2009. Nabeel Seedat, Fergus Imrie, Alexis Bellot, Zhaozhi Qian, and Mihaela van der Schaar. Continuoustime modeling of counterfactual outcomes using neural controlled differential equations. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 19497–19521. PMLR, 17–23 Jul 2022. Fangzheng Sun, Yang Liu, Jian-Xun Wang, and Hao Sun. Symbolic physics learner: Discovering governing equations via Monte Carlo tree search. In The Eleventh International Conference on Learning Representations, 2023. Wassim Tenachi, Rodrigo Ibata, and Foivos I Diakogiannis. Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws. The Astrophysical Journal, 959(2):99, 2023. Yuan Tian, Wenqi Zhou, Michele Viscione, Hao Dong, David S Kammer, and Olga Fink. Interactive symbolic regression with co-design mechanism through offline reinforcement learning. Nature Communications, 16(1):3930, 2025. Jean-François Toubeau, Jérémie Bottieau, François Vallée, and Zacharie De Grève. Deep learningbased multivariate probabilistic forecasting for short-term scheduling in power markets. IEEE Transactions on Power Systems, 34(2):1203–1215, 2018. Silviu-Marian Udrescu and Max Tegmark. Ai feynman: A physics-inspired method for symbolic regression. Science advances, 6(16):eaay2631, 2020. Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. SymbolicGPT: A generative transformer model for symbolic regression. arXiv preprint arXiv:2106.14131, 2021. Martin Vastl, Jonáš Kulhánek, Jiří Kubalík, Erik Derner, and Robert Babuška. Symformer: End-toend symbolic regression using transformer-based architecture. IEEE Access, 12:37840–37849, 2024. Alistair White, Niki Kilbertus, Maximilian Gelbrecht, and Niklas Boers. Stabilized neural differential equations for learning dynamics with explicit constraints. Advances in Neural Information Processing Systems, 36:12929–12950, 2023. Philippe Wyder, Judah Goldfeder, Alexey Yermakov, Yue Zhao, Stefano Riva, Jan Williams, David Zoro, Amy Rude, Matteo Tomasetto, Joe Germany, et al. Common task framework for a critical evaluation of scientific machine learning algorithms. Advances in Neural Information Processing Systems, 38, 2026. Yuan Yin, Vincent Le Guen, Jérémie Dona, Emmanuel De Bézenac, Ibrahim Ayed, Nicolas Thome, and Patrick Gallinari. Augmenting physical models with deep networks for complex dynamics forecasting. Journal of Statistical Mechanics: Theory and Experiment, 2021(12):124012, 2021. Bob Junyi Zou, Matthew E Levine, Dessi P. Zaharieva, Ramesh Johari, and Emily Fox. Hybrid$^2$ neural ODE causal modeling and an application to glycemic response. In Forty-first International Conference on Machine Learning, 2024.
13
A
Why L2 Regularization Fails
Here we analyze when and why L2 regularization guarantees orthogonality. We show that APHYNITY’s approach works for fixed symbolic components (Theorem A.1), but fails when symbolic structure is learned jointly with sparsity constraints (Theorem A.2). Problem Setup. Throughout this appendix, F is a linear function space (closed under sums and scalar multiplication) sufficiently expressive to contain f − fphy for the fphy ∈ Fphy at hand; this is what allows the projection arguments below. Consider learning a hybrid model fˆ = fˆphy + fˆaug where: PM • fˆphy (x; w) = j=1 wj ϕj (x) with w ∈ RM (symbolic component). • fˆaug (x; ϑ) is a neural network (neural component); the family {fˆaug (·; ϑ)}ϑ admits a parameter ϑ0 for which fˆaug (·; ϑ0 ) ≡ 0 (e.g. all final-layer weights set to zero for an MLP without an additive skip connection at the output). The standard L2 -regularized loss is: LL2 (w, ϑ) = ∥f − fˆphy − fˆaug ∥2D + µ∥w∥1 + λ2 ∥fˆaug ∥2D | {z } | {z } | {z } reconstruction
sparsity
(10)
L2 regularization
When L2 Regularization Guarantees Orthogonality. We first establish the positive result: L2 regularization works as intended when the symbolic component is fixed. This is a key result of [Yin et al., 2021]. Proposition A.1 (Fixed symbolic hypothesis class). Assume Fphy = span{ϕ1 , . . . , ϕM } is fixed, faug is optimized over F, and λ2 > 0. Consider min
fphy ∈Fphy , faug ∈F
∥f − fphy − faug ∥2D + λ2 ∥faug ∥2D .
1 Then any minimizer satisfies fphy = PFDphy (f ) and faug = 1+λ (f − fphy ), hence faug ⊥ Fphy . 2
Proof. For fixed fphy , the objective is quadratic in faug and is minimized, up to D-equivalence, by 1 (f − fphy ). Plugging back reduces the problem to faug = 1+λ 2 min
fphy ∈Fphy
∥f − fphy ∥2D ,
whose minimizer is the orthogonal projection PFDphy (f ). The residual f − fphy is orthogonal to Fphy , and scaling preserves orthogonality. When L2 Regularization Fails. The situation changes when we learn the symbolic structure itself. The following theorem shows that joint optimization breaks the orthogonality guarantee. Theorem A.2 (L2 failure with sparse symbolic learning). There exist f and a fixed symbolic library ˆ∗ {ϕj }M j=1 such that the minimizer of Equation (10) is not orthogonal, that is, faug ̸⊥ Fphy , even though f ∈ Fphy . Proof. Let ϕ1 , ϕ2 be orthonormal w.r.t. ⟨·, ·⟩D and set f = ϕ1 + ϕ2 , so f ∈ Fphy := span{ϕ1 , ϕ2 }. Consider the (function-space) version of Equation (10), with g standing in for fˆaug to lighten notation: min
w∈R2 , g∈F
∥f − (w1 ϕ1 + w2 ϕ2 + g)∥2D + µ∥w∥1 + λ2 ∥g∥2D .
1 For fixed w, the unique minimizer in g is g(w) = 1+λ (f − w1 ϕ1 − w2 ϕ2 ), and the objective reduces 2 to the Lasso problem
min λ2 ∥f − (w1 ϕ1 + w2 ϕ2 )∥2D + µ∥w∥1 . w∈R2 1+λ2 14
2λ2 Since ϕ1 , ϕ2 are orthonormal and f = ϕ1 + ϕ2 , the problem decouples and for 0 < µ < 1+λ the 2 unique minimizer is 2) w1∗ = w2∗ = 1 − µ(1+λ 2λ2 .
Thus g ∗ = g(w∗ ) = 2λµ2 (ϕ1 + ϕ2 ), so ⟨g ∗ , ϕj ⟩D = 2λµ2 ̸= 0 for j ∈ {1, 2}. Hence g ∗ ̸⊥ Fphy although f ∈ Fphy .
B
Theoretical Guarantees of OrthoReg
Having established that L2 regularization can fail under joint optimization with sparsity constraints, we now develop the theory of orthogonal regularization. We show that explicit orthogonality penalties yield quantitative control of neural–symbolic overlap and provide useful geometric structure. The Orthogonality Constraint.
OrthoReg solves:
min Lfit (w, ϑ) + µ∥w∥1 + λ w,ϑ
M X ⟨fˆaug (·; ϑ), ϕj ⟩2D ,
(11)
j=1
where Lfit is any non-negative data-fit term (e.g., Equations (4) and (5)). Unlike L2 regularization, OrthoReg penalizes the squared inner products ⟨fˆaug , ϕj ⟩2D rather than the total norm ∥fˆaug ∥2D . Penalty method: orthogonality at optimum. Theorem B.1 (Penalty bound at optimum). Assume λ > 0. Let (w∗ , ϑ∗ ) be a global minimizer of Equation (11) and let wsym ∈ arg min Lfit (w, ϑ0 ) + µ∥w∥1 w
be the best pure symbolic fit, where fˆaug (·; ϑ0 ) ≡ 0. Then M X
⟨fˆaug (·; ϑ∗ ), ϕj ⟩2D ≤
j=1
1 Lfit (wsym , ϑ0 ) + µ∥wsym ∥1 . λ
In particular, for any sequence of global minimizers with λ → ∞, the empirical orthogonality violations converge to zero. Proof. By optimality of (w∗ , ϑ∗ ), we have Lfit (w∗ , ϑ∗ ) + µ∥w∗ ∥1 + λ
M X
2 fˆaug (·; ϑ∗ ), ϕj D
j=1
⩽ Lfit (wsym , ϑ0 ) + µ∥wsym ∥1 , where fˆaug (·; ϑ0 ) ≡ 0 by definition. Dropping the non-negative first two terms on the left yields the stated bound. The final claim follows by sending λ → ∞. Remark (uniqueness up to D-equivalence). The empirical ⟨·, ·⟩D is a seminorm on F: functions agreeing on D are identified. The minimisers in Theorem A.1 and Theorem 4.1 are therefore unique only up to D-equivalence; this matters most when D does not separate the augmentation hypothesis class on X , and is mild in our reported regimes where the empirical samples are dense. Approximation Quality. Beyond enforcing orthogonality, we want to understand what orthogonality buys in terms of approximation structure. We decompose the error into symbolic and orthogonal-complement components with the orthogonal decomposition: 15
Theorem B.2 (Orthogonal error decomposition). Let Fphy be a closed linear subspace of (F, ⟨·, ·⟩D ) and let PFDphy denote the orthogonal projection onto Fphy with respect to ⟨·, ·⟩D . For any fphy ∈ Fphy and any faug ∈ F satisfying faug ⊥ Fphy , we have ∥f − (fphy + faug )∥2D = ∥PFDphy (f ) − fphy ∥2D + ∥f − PFDphy (f ) − faug ∥2D . Proof. Write f = PFDphy (f ) + (f − PFDphy (f )) with PFDphy (f ) ∈ Fphy and f − PFDphy (f ) ⊥ Fphy . Then f − (fphy + faug ) = PFDphy (f ) − fphy + f − PFDphy (f ) − faug , ⊥ where the first term lies in Fphy and the second term lies in Fphy , hence they are orthogonal. Expanding the squared norm yields the identity.
B.1
Finite-Sample Analysis
Our theoretical analysis so far has operated under the empirical inner product. In finite-sample settings we enforce orthogonality empirically, and here we quantify how well empirical inner products approximate their population counterparts. Empirical vs Population Orthogonality. Definition B.3 (Orthogonality Gap). Define the empirical inner product N 1 X f (xi )⊤ g(xi ), ⟨f, g⟩D = N i=1 and the population inner product ⟨f, g⟩ν = Ex∼ν [f (x)⊤ g(x)]. The orthogonality gap is |⟨f, g⟩D − ⟨f, g⟩ν |. The following concentration result applies to fixed functions f and g, that is, functions that do not depend on the specific sample D. Theorem B.4 (Finite-Sample Guarantee). Let f, g : X → Rn be fixed functions and assume |f (x)⊤ g(x)| ≤ R almost surely. Let D = {x1 , . . . , xN } be i.i.d. samples from distribution ν. Then with probability at least 1 − δ: r 2 log(2/δ) |⟨f, g⟩D − ⟨f, g⟩ν | ≤ R . N Proof. Define random variables Zi = f (xi )⊤ g(xi ) for i = 1, . . . , N , where xi ∼ ν are i.i.d. and satisfy Zi ∈ [−R, R] almost surely. Then N
⟨f, g⟩D =
1 X Zi , N i=1
⟨f, g⟩ν = E[Zi ].
Hoeffding’s inequality states that for i.i.d. random variables Zi ∈ [a, b], ! N 1 X 2N ϵ2 P Zi − E[Zi ] ≥ ϵ ≤ 2 exp − . N i=1 (b − a)2 Applying this with a = −R, b = R and solving 2 exp(−2N ϵ2 /4R2 ) = δ yields p ϵ = R 2 log(2/δ)/N .
16
Corollary B.5 (Population Orthogonality). If we enforce ⟨fˆaug , ϕj ⟩D = 0 empirically, then under the assumption that fˆaug is fixed independently of D, we have r 2 log(2/δ) ˆ |⟨faug , ϕj ⟩ν | ≤ R N with probability 1 − δ. Intuition.
Each inner product N
⟨fˆaug , ϕj ⟩D =
1 Xˆ faug (xi )⊤ ϕj (xi ) N i=1
is an average of N bounded random variables. Hoeffding’s inequality√for bounded variables Zi ∈ [−R, R] yields concentration around the population mean at rate O(1/ N ). Remark. Since fˆaug is trained on the same sample used to evaluate ⟨fˆaug , ϕj ⟩D , the independence assumption is violated. A population-level guarantee for the learned augmentation therefore requires uniform convergence over the augmentation hypothesis class (e.g. Rademacher or VC-style bounds), and not single-function concentration; we discuss this in Section I. B.2
Monte Carlo Approximation
In training we estimate the empirical inner products in the orthogonality penalty using minibatches. We quantify the resulting stochastic approximation error. Minibatch estimator. Let D = {xi }N i=1 and fix a constraint ϕj . For a minibatch B ⊆ {1, . . . , N } of size B, define 1 Xˆ ⟨fˆaug , ϕj ⟩B := faug (xi )⊤ ϕj (xi ), B i∈B N
⟨fˆaug , ϕj ⟩D :=
1 Xˆ faug (xi )⊤ ϕj (xi ). N i=1
Assume |fˆaug (x)⊤ ϕj (x)| ≤ R for all x ∈ D. Lemma B.6 (Minibatch deviation). If B is sampled uniformly (with replacement) from D, then for any δ ∈ (0, 1), Bϵ2 P ⟨fˆaug , ϕj ⟩B − ⟨fˆaug , ϕj ⟩D ≥ ϵ ≤ 2 exp − 2 , 2R and equivalently, with probability at least 1 − δ, r 2 log(2/δ) ˆ ˆ ⟨faug , ϕj ⟩B − ⟨faug , ϕj ⟩D ≤ R . B The same scaling holds for sampling without replacement (up to a finite-population correction factor). Implication for optimization. Theorem B.6 shows that the minibatch inner products concentrate around their full-data counterparts at rate O(B −1/2 ). In the default configuration the orthogonality penalty is evaluated on the full training set each step (all trajectories fit into one batch). For completeness, we record a standard concentration bound for minibatch estimation.
C
When Residual Pipelines Can Reintroduce Symbolic Overlap
We revisit the setting of Theorem 1.1. Let ϕ1 , ϕ2 ∈ F be orthonormal with respect to ⟨·, ·⟩D , let f = ϕ1 + ϕ2 , and let the Stage-1 symbolic model be fitted with an L1 penalty of weight µ ∈ (0, 2). The orthonormal Lasso solution is ŵj = 1 − µ/2, 17
j = 1, 2,
and therefore leaves the empirical residual r := f −
2 X
ŵj ϕj =
j=1
µ (ϕ1 + ϕ2 ) ∈ span{ϕ1 , ϕ2 }. 2
Thus, the residue created by sparsity regularization lies exactly in the symbolic span. The question is whether a neural residual fitted after, or through, the symbolic component is forced to remain orthogonal to this span. Proposition C.1 (Unconstrained residual learning and compositional pipelines can retain symbolic overlap). In the setting above, the following hold. 1. Unconstrained residual learning. If w is fitted by Lasso and an unconstrained augmentation g ∈ F is subsequently fitted to the residual, then ĝ = r,
⟨ĝ, ϕj ⟩D = µ/2,
j = 1, 2.
2. Compositional architecture with an input-space penalty. Consider a compositional architecture fˆ(x) = (I + Nϑ )(Tw (x)) = Tw (x) + Nϑ (Tw (x)),
Tw (x) =
2 X
wj ϕj (x).
j=1
For this simplified example, Tw denotes the symbolic component itself; in compositional models with a residual structure fˆ = Tw + Nϑ ◦ Tw , Tw is the physics map onto which Nϑ composes. Define the effective additive augmentation at the Lasso solution ŵ eϑ (x) := Nϑ (Tŵ (x)). N Assume that, at the Lasso solution ŵj = 1 − µ/2, the 2N points N {xi }N i=1 ∪ {Tŵ (xi )}i=1
are all distinct. Then there exists a sufficiently expressive Nϑ such that fˆ(xi ) = f (xi ), for all i = 1, . . . , N . Consequently, X Nϑ (xi )⊤ ϕj (xi ) = 0,
Nϑ (xi ) = 0,
j = 1, 2,
but
eϑ (xi ) = r(xi ) N
eϑ , ϕj ⟩D = ⟨r, ϕj ⟩D = µ/2. ⟨N
i
Thus, an OrthoReg-style penalty applied to Nϑ at the original inputs {xi } does not control the additive overlap of the effective augmentation Nϑ ◦ Tŵ . Proof. For unconstrained residual learning, orthonormality gives the usual soft-thresholding solution ŵj = (1 − µ/2)+ . Since µ ∈ (0, 2), this gives ŵj = 1 − µ/2, and an unconstrained augmentation fits the residual exactly: 2 X µ ĝ = f − ŵj ϕj = (ϕ1 + ϕ2 ). 2 j=1 Taking empirical inner products with ϕ1 and ϕ2 gives ⟨ĝ, ϕj ⟩D = µ/2. For the compositional case, exact fit at the Lasso coefficients requires Nϑ (Tŵ (xi )) = f (xi ) − Tŵ (xi ) = r(xi ). This constrains Nϑ at the transported points Tŵ (xi ), not at the original data points xi . Since the 2N points are distinct, a sufficiently expressive Nϑ can interpolate both Nϑ (Tŵ (xi )) = r(xi )
and
Nϑ (xi ) = 0.
The empirical penalty on Nϑ (xi ) is then zero, while the effective additive augmentation satisfies eϑ (xi ) = Nϑ (Tŵ (xi )) = r(xi ), N eϑ , ϕj ⟩D = µ/2. and hence ⟨N 18
Orthogonality-regularised residual learning. The failure above is not caused by staging itself, but by fitting the residual without an orthogonality constraint. If Stage 2 instead minimises X ∥r − g∥2D + λ ⟨g, ϕj ⟩2D , j
then, on the orthonormal toy, ⟨ĝ, ϕj ⟩D = (µ/2)/(1+λ) → 0 as λ → ∞. Orthogonalisation therefore suppresses symbolic–neural overlap whether imposed jointly or in an explicitly regularised residual stage; it does not remove the L1 shrinkage bias on w, so the shrinkage residue is left as fit error rather than absorbed by the neural component.
D
Computational Complexity
We briefly discuss the additional computational overhead introduced by the orthogonality penalty. Per-iteration cost. Let B be the batch size, n the state dimension, and M the number of library terms. For the vector-field regression loss used in the main experiments, computing the penalty requires evaluating ϕj (xi ) and forming inner products fˆaug (xi )⊤ ϕj (xi ) for all i ∈ [B] and j ∈ [M ], which costs O(B n M ) arithmetic operations in addition to the usual network forward/backward passes; backpropagation has the same order. Rollout or sequence-prediction losses multiply this cost by the number of evaluated transitions or time steps. Timing protocol. Wall-clock per epoch is measured as the median over a fixed window of training epochs after a short warm-up, on the modified damped pendulum (B = 1000, n = 2, M = 20), with identical optimizer state and data loaders across configurations. The reported 5–15% overhead is the relative increase of the OrthoReg median over the Hybrid L2 median across five seeds, isolating the cost of the orthogonality term gj ; overhead grows linearly in M at fixed (B, n). Memory requirements. The additional memory is dominated by storing the evaluated library features on the batch, which scales as O(B n M ) floats. For example, with B = 1000, n = 2, and M = 20, this corresponds to approximately 4×104 scalars (about 160 KB in float32); on the SIR system (M = 45) the same calculation gives roughly 360 KB, well below the network activations cached for backpropagation. P Mini-batch gradient bias. The penalty j ⟨fˆaug , ϕj ⟩2D is a squared mean of bounded random variables, so the standard one-batch unbiased Monte-Carlo estimator of ∇ϑ ⟨fˆaug , ϕj ⟩D (which would multiply two independent mini-batch estimates) is replaced in our implementation by the squared single-batch average. Under stochastic gradient descent this introduces an additive bias of order O(1/B) relative to the full-data quantity (exactly the variance of the per-batch inner-product estimate); the full-batch evaluation used in the main experiments (B equal to the full training set, see Section G) makes this bias zero by construction, and the mini-batch ablation in Section E.3 confirms that the ranking and quantitative trends carry over for B ≥ 200. A two-batch unbiased estimator is a drop-in alternative for memory-constrained settings. Hardware and total compute. All reported experiments were run on an internal HPC cluster, allocating one general-purpose GPU per job, ≈ 64 GB of RAM, and a small CPU allocation (typically ≤ 8 cores). The workloads are intentionally small—autonomous ODEs in n ≤ 4 states, MLPs of width 64–128, 2000 Adam epochs, and batch sizes ≤ 1000—so that OrthoReg remains CPUtractable; we used GPUs primarily to accelerate the neural ODE component, which dominates wall-clock time. Each individual training run completes within a small number of GPU-minutes; the full reported grid (main systems, ablations, and cross-system validation, each repeated over five seeds) accounts for the bulk of the project’s compute, and preliminary tuning and discarded runs add a similar order-of-magnitude overhead on top.
19
E
Additional Ablations
E.1
Full Dataset Difficulty Ablation
Table 4: Dataset difficulty ablation across missing dynamics regimes (mean effect strength absent from the symbolic library). OrthoReg’s gains are largest at medium difficulty; in the low-difficulty regime it trades ID and OOD-T2 accuracy for higher F1 and OOD-T3 generalization. Metric
Pure
L2
OrthoReg
Low Difficulty (β1 = 0.1, β2 = 0.08, β3 = 0.05) MSE(ẋID , e−4) (↓) MSE(ẋOOD,T2 , e−3) (↓) MSE(ẋOOD,T3 , e0) (↓) F1 (↑)
7.77±1.64 1.53±0.26 161.60±166.18 0.50±0.09
10.99±0.83 2.01±0.00 5.86±0.01 0.72±0.05
26.82±5.00 2.79±0.09 0.58±0.06 0.86±0.00
Medium Difficulty (β1 = 0.3, β2 = 0.25, β3 = 0.15) MSE(ẋID , e−2) (↓) MSE(ẋOOD,T2 , e−2) (↓) MSE(ẋOOD,T3 , e0) (↓) F1 (↑)
1.27±0.07 1.15±0.05 940.60±614.71 0.43±0.09
1.41±0.05 1.27±0.04 218.49±111.32 0.61±0.03
1.89±0.01 1.54±0.02 2.45±0.43 0.93±0.15
High Difficulty (β1 = 2.0, β2 = 2.0, β3 = 2.0) MSE(ẋID , e−2) (↓) MSE(ẋOOD,T2 , e−2) (↓) MSE(ẋOOD,T3 , e−1) (↓) F1 (↑)
3.91±0.13 4.49±0.17 4.49±2.61 0.46±0.05
3.90±0.11 4.46±0.10 4.06±0.75 0.43±0.04
3.86±0.09 4.43±0.09 2.64±0.48 0.52±0.07
† For Medium, OOD T3 reports derivative MSE instead of state MSE due to numerical instabilities during trajectory
integration for baseline methods. Low and High report state MSE for OOD T3.
Table 4 reports the full dataset difficulty ablation across regimes of increasing missing dynamics (controlled by the mean effect strength β of terms absent from the symbolic library; cf. βi in Equation (9)). In the main paper, Figure 2a summarizes the central trends using stress-test generalization (xOOD,T3 ) and symbolic recovery (F1). The complete table additionally includes in-distribution derivative accuracy (ẋID ) and a milder OOD evaluation without parameter shifts (xOOD,T2 ). Low difficulty (β1 = 0.1, β2 = 0.08, β3 = 0.05). When the true dynamics are close to the symbolic library, all methods achieve low errors on in-distribution and mild OOD evaluations (Table 4, first block), and the hybrid component is not strongly required. Interestingly, the most flexible baselines can remain competitive on these simpler evaluations, consistent with the observation that regularization is least critical when the library is largely sufficient. Under the stress-test setting (xOOD,T3 ), however, OrthoReg remains substantially more robust (Figure 2a), indicating that enforcing a clean split can mitigate extrapolation failures even when the library mismatch is small. Medium difficulty (β1 = 0.3, β2 = 0.25, β3 = 0.15). This regime shows the most pronounced separation between methods. OrthoReg simultaneously improves OOD performance and symbolic recovery, yielding markedly lower xOOD,T3 error and substantially higher F1 than both the unregularized baseline and L2 regularization (Figure 2a; Table 4, second block). Conceptually, the orthogonality penalty constrains the augmentation to capture residual structure without redundantly re-learning library terms, stabilising identification and improving generalisation. High difficulty (β1 = 2.0, β2 = 2.0, β3 = 2.0). When the missing dynamics strongly exceed the representational capacity of the library, all methods deteriorate in symbolic recovery and OOD performance (Table 4, third block), reflecting that the symbolic–neural split becomes intrinsically harder to recover. Nevertheless, OrthoReg retains a modest advantage on xOOD,T3 and F1 (Figure 2a). In this regime, once the residual dominates, orthogonality alone is insufficient to recover interpretable structure and the symbolic–neural split becomes harder to identify. Overall, Table 4 complements Figure 2a by showing that OrthoReg’s strongest gains arise when the system partially exceeds the library: small mismatches do not demand strong regularization, whereas extreme mismatches limit recoverability for all methods. 20
Table 5: Sampling scheme ablation (regular vs. irregular). All methods degrade under irregular sampling; OrthoReg retains a small relative advantage on F1. Metric
L2
OrthoReg
1.41±0.05 1.07±0.02 7.90±4.02 0.61±0.03
1.89±0.01 1.01±0.01 0.02±0.00 0.93±0.15
Pure Regular Sampling
MSE(ẋID , e−2) (↓) MSE(xOOD,T2 , e0) (↓) MSE(ẋOOD,T3 , e2) (↓) F1 (↑)
1.27±0.07 1.02±0.07 34.00±22.22 0.43±0.09
Irregular Sampling MSE(ẋID , e0) (↓) MSE(xOOD,T2 , e0) (↓) MSE(ẋOOD,T3 , e1) (↓) F1 (↑)
E.2
3.52±0.00 3.85±0.00 4.04±0.09 0.30±0.01
3.52±0.00 3.85±0.00 3.90±0.10 0.30±0.01
3.52±0.00 3.85±0.00 3.74±0.05 0.31±0.01
Full Sampling Scheme Ablation
Table 5 reports the full sampling scheme ablation across regular and irregular observation settings. Irregular sampling substantially increases derivative and OOD prediction errors and reduces symbolic recovery across all methods, consistent with increased solver error and gradient variance under nonuniform time steps (see Section G). Despite this degradation, OrthoReg retains a marginal advantage in symbolic recovery (F1) compared to baseline methods. These results complement the trends shown in Figure 2b and indicate that orthogonal regularization remains effective for symbolic recovery beyond idealized sampling assumptions. E.3
Monte Carlo Sampling Ablation
We investigate the impact of Monte Carlo sampling on model performance by varying the number of training samples from 100 to 5000. Figure 4 shows the performance across different sample sizes for the medium missing dynamics regime (β1 = 0.3, β2 = 0.25, β3 = 0.15).
Figure 4: Monte Carlo sampling ablation study (medium missing dynamics). Performance is shown across different sample sizes (100–5000) for F1 score, OOD, T2 MSE, and a cosine diagnostic ⟨fˆaug , fˆphy ⟩D /(∥fˆaug ∥D ∥fˆphy ∥D ) between the learned augmentation and the assembled symbolic P prediction fˆphy = j ŵj ϕj on the empirical training distribution. The cosine takes values in [0, 1]: smaller values mean fˆaug and fˆphy are more directionally separated under the empirical inner product, larger values mean they share direction. We report it to visualise the concentration result of Theorem B.6: at rate O(N −1/2 ) the empirical inner products that drive the OrthoReg penalty concentrate around their full-data values, and the cosine – built from those same inner products – stabilises along with them as N grows. The cosine is a within-method stability diagnostic for a single OrthoReg run; it is built from the empirical inner products that appear in the OrthoReg training P penalty j ⟨fˆaug , ϕj ⟩2D but is not the penalty itself. F1 and the cosine diagnostic both stabilise as the number of training samples grows for OrthoReg, consistent with the O(N −1/2 ) concentration of the empirical inner products that the OrthoReg penalty operates on (Theorem B.6). 21
E.4
Noise Robustness Ablation
We investigate the robustness of OrthoReg, L2 hybrid, and Pure Symbolic methods to measurement noise in the training data. Gaussian noise with standard deviation ε ∈ {0.0, 0.01, 0.05, 0.1, 0.2, 0.5} is added to the state observations after trajectory generation, simulating realistic measurement uncertainty. OOD T2 Deriv MSE (Fixed reg_scale=0.85)
101 100
F1 Score (Fixed reg_scale=0.85)
OrthoReg L2 Pure Symbolic
OrthoReg L2 Pure Symbolic
0.9 0.8 0.7
102
F1 Score
OOD T3 Deriv MSE
OOD T2 Deriv MSE
102
OOD T3 Deriv MSE (Fixed reg_scale=0.85)
103
OrthoReg L2 Pure Symbolic
0.6 0.5 0.4
101
0.3 10 1
0.2 0.0
0.1
0.2
0.3
0.4
0.5
0.0
0.1
0.2
0.3
Noise Level ( )
0.4
0.5
0.0
0.1
0.2
0.3
0.4
0.5
Figure 5: Noise robustness study on the damped pendulum system. Performance degrades similarly across all methods as noise increases. OrthoReg maintains its advantage at zero noise (superior OOD, T2 and OOD, T3 performance, higher F1 score) but does not provide special robustness to noise compared to baselines. All methods converge to similar performance levels at high noise (ε ≥ 0.1), indicating that measurement noise fundamentally limits symbolic recovery regardless of regularization approach. Curves correspond, in order, to OrthoReg, L2 -regularised hybrid, and Pure Symbolic. Figure 5 shows that all three methods suffer similarly from increasing noise levels. OrthoReg maintains its performance advantages at zero noise (lower OOD, T2 and OOD, T3 derivative MSE, higher F1 score), but these advantages diminish rapidly as noise increases. At noise levels ε ≥ 0.1, all methods converge to similar performance, with OOD, T2 derivative MSE increasing from ∼ 0.05 (zero noise) to ∼ 0.6 (high noise) and F1 scores dropping from ∼ 0.9 to ∼ 0.2 − 0.4 across all methods. Overall, orthogonal regularization does not materially improve robustness to measurement noise in this setting; performance converges across methods at higher noise levels. Improving robustness likely requires changes to the estimation procedure (e.g., denoising or noise-aware differentiation) rather than additional separation constraints. At high noise (ε ≳ 0.1), the smoothed-finite-difference targets used to compute yi ≈ ẋ(ti ) are themselves dominated by amplified noise (Savitzky–Golay window of length 5), so the convergence of all methods at this regime partly reflects target corruption rather than failure of the regularizer itself; an analytic-derivative ablation would disentangle the two and is left to future work (Section I).
F
Cross-System Validation
The main experiments (Section 5) evaluate OrthoReg on the modified pendulum and Duffing systems. Here we assess transfer to additional systems that introduce different sources of library mismatch, focusing on Lotka–Volterra (temporal modulation) and SIR (state dependence and memory). Together, these systems test three distinct mismatch regimes: structured feature mismatch (Duffing), nonautonomous perturbations (Lotka–Volterra), and strongly misspecified state-dependent residual structure (SIR).
Autonomous augmentation with non-autonomous data. Both LV and SIR include explicit timedependent factors in the ground-truth dynamics, while the augmentation fˆaug (x; ϑ) in Equation (7) takes only the state x as input. An autonomous augmentation can therefore fit such factors only as a state-dependent envelope on the training horizon [0, T ]; the OOD splits below test state extrapolation rather than time extrapolation. Concatenating t with x in the network input is a one-line architectural change that would close this gap. 22
F.1
Lotka–Volterra System: Temporal Coupling
We evaluate OrthoReg on a modified predator-prey system with temporally modulated and statedependent interactions. The dynamics are: dx = αx − βxy + ε1 x sin(ωfast t) cos(ωfast xy) × sin ωslow (x + y) , dt x dy . = δxy − γy + ε2 y sin(ωfast t) cos(ωfast xy) × sin ωslow (x + y) sin dt y + creg Here, ε1 and ε2 control the strength of missing dynamics not captured by the symbolic feature library, while creg is a small constant that prevents division by zero. The runs reported below use (ε1 , ε2 ) = (3.0, 2.0) (see the per-system reproducibility values in Section G). The augmented terms introduce high-frequency temporal modulation, state-dependent coupling, and asymmetric predatorprey interactions; the explicit sin(ωfast t) factors make the system non-autonomous (Section F). Table 6: Lotka–Volterra results including derivative metrics. Scale factors shown inside the metric parentheses indicate the power of ten applied to the reported values (e.g., e−2). Uncertainty terms reported as ±0.0 fall below the displayed precision after rescaling. Abbreviations: MSE = mean squared error; ID = in-distribution; OOD = out-of-distribution. Metric
L2
Pure
OrthoReg
Predictive Performance MSE(ẋID , e−2) (↓) MSE(ẋOOD,T2 , e−2) (↓) MSE(ẋOOD,T3 , e−1) (↓)
1.64 ±0.00 1.21 ±0.00 1.74 ±0.00
1.64 ±0.00 1.22 ±0.00 1.73 ±0.00
1.64 ±0.00 1.22 ±0.00 1.71 ±0.00
System Identification Quality F1 (↑) #Terms (↓)
0.22 ±0.01 16.6 ±0.9
0.22 ±0.00 16.0 ±0.0
0.24 ±0.01 14.8 ±0.4
Results and Analysis. Across all predictive metrics in Table 6, the three methods perform similarly in this configuration, indicating that the chosen library and training distribution already capture most of the dynamics relevant to our evaluation protocol. OrthoReg yields a small but consistent gain in identification quality (F1: 0.24 vs. 0.22) while selecting slightly fewer terms. We therefore view Lotka–Volterra as a “low headroom” setting: it serves as a sanity check that OrthoReg does not degrade performance when the decomposition problem is comparatively easy, rather than as a regime with large expected improvements. F.2
Duffing oscillator: multistability and cross-basin generalization
This appendix complements the main-text results in Table 3 by detailing the experimental setup and the origin of the cross-basin generalization challenge. Our Duffing configuration and the motivation for the basin split follow Göring et al. [2024]. System and regimes.
We consider the unforced Duffing oscillator in first-order form, ẋ = y,
(12) 2
ẏ = ay − x(b + cx ),
(13)
1 with parameters [a, b, c] = [− 21 , −1, 10 ], placing the dynamics in a bistable regime with two p symmetric basins of attraction at x = ± −b/c (and y = 0). Göring et al. [2024] use this multistable setting to study out-of-domain generalization for learned dynamical-system models.
Why cross-basin generalization is challenging. When training data are restricted to a single basin, the model never observes trajectories from the opposite basin. For multistable flows that are not topologically transitive on the test domain, Göring et al. [2024] formalize that the resulting OOD generalization problem is not strictly learnable for universal-approximator hypothesis classes (Theorem 4.2), and empirically illustrate this behavior on the Duffing system. 23
Library mismatch. The shared symbolic library (polynomial degree 2 ⊕ Fourier(n = 1); Section G) does not contain the cubic term cx3 , so any method using this library must either (i) approximate the missing nonlinearity through other terms (pure symbolic) or (ii) represent it through the neural augmentation (hybrid models). Effect of orthogonal regularization. OrthoReg encourages the neural augmentation to be (approximately) orthogonal to the span of library terms on the training distribution, encouraging the symbolic component to retain the structure expressible by the library while pushing the neural component toward dynamics not captured by it. Empirically, this is associated with improved cross-basin performance and more stable identification compared to L2 -regularised hybrids (Table 3). Numerical note on basin crossings. Although the unforced system is basin-confining in the idealized continuous dynamics, numerical integration and initial conditions close to the separatrix can yield occasional apparent basin crossings. In the qualitative visualization we therefore treat trajectories as illustrative, and we generate/inspect trajectories to ensure the intended train-test basin split is respected (cf. Section G for integration details). F.3
SIR: State Dependence and Memory
System Design. We extend the classical normalised SIR model with a multiplicative tweaked term per compartment that combines a state-dependent bilinear coupling, an exponentially decaying time envelope, and a state-modulated frequency. Concretely, with ωS = ωfast (1 + S), ωI = ωfast (1 + 2I), ωR = ωslow (1 + R), SI dS = −β + ε1 S I sin(ωS t) e−t/τmemory tanh ωslow (I + R) , dt S+I +R dI SI = +β − γ I + ε2 I R sin(ωI t) e−t/τmemory tanh ωslow (S + R) , dt S+I +R dR = +γ I + ε3 S R sin(ωR t) e−t/τmemory tanh ωslow (S + I) , dt
(14) (15) (16)
with parameters β = 2.0, γ = 1.0, (ε1 , ε2 , ε3 ) = (0.3, 0.25, 0.2), ωfast = 3.0, ωslow = 1.5, τmemory = 2.0 (Section G). Each tweaked term is added (sign +) to the corresponding standard SIR derivative, so εi → 0 recovers the classical normalised SIR model exactly. The base rates β, γ are constant; all additional state dependence enters through the tweaked residual terms [Hethcote, 2000, Kucharski et al., 2020]. As in LV, the explicit time factors make this a non-autonomous stress test for the autonomous augmentation (Section F). Table 7: SIR results under severe library mismatch. OrthoReg improves OOD-T3 extrapolation and symbolic separation, but trades off ID and OOD-T2 accuracy. Metric
Pure
L2
OrthoReg
Predictive Performance MSE(ẋID , e−5) (↓) MSE(ẋOOD,T2 , e−4) (↓) MSE(ẋOOD,T3 , e−4) (↓)
2.12±0.73 2.40±0.11 72.68±19.04
31.42±4.93 7.75±0.53 15.33±3.33
90.64±7.56 16.54±0.80 8.76±0.52
System Identification Quality F1 (↑) #Terms (↓) Orth. (↑)
0.14±0.03 40.2±3.1 –
0.15±0.03 20.4±3.6 0.35±0.16
0.20±0.06 6.8±3.3 0.55±0.23
The SIR system demonstrates a challenging regime where the symbolic library has a severe mismatch with the true dynamics, a large number of candidate terms are present, and the tweaked term combines state-modulated frequencies sin(ωk (·)t) with a decaying time envelope e−t/τmemory that the autonomous augmentation can absorb only as a state-dependent envelope (Section F). While all methods are challenged by this setup, OrthoReg retains an advantage in the out-of-distribution setting (OOD-T3) and in system identification (highest F1, fewest terms). 24
Summary. Taken together, the cross-system experiments show that OrthoReg remains most useful when the library captures part of the dynamics but leaves structured residuals. The gains are modest in low-headroom settings such as LV, clearer for cross-basin Duffing generalisation, and mixed under the strongly misspecified SIR variant, where OrthoReg improves OOD-T3 and sparsity but sacrifices ID and OOD-T2 accuracy.
G
Experimental Design and Implementation
Overview. All experiments follow a common protocol: we generate synthetic trajectories from known ground-truth dynamical systems, train models to predict the vector field from state observations, and evaluate predictive accuracy, symbolic recovery, and neural–symbolic separation. To isolate the effect of orthogonal regularization, OrthoReg and the hybrid baselines share the same model structure and training pipeline; they differ only in the form of the regularization applied to the neural augmentation. G.1
Evaluation metrics
We report two complementary result metrics: • Predictive performance. Mean-squared error (MSE) on derivatives, computed against the ground-truth vector field (available for all synthetic systems), and, when numerically stable, MSE of long-horizon trajectory rollouts in state space. • System identification. F1 score comparing the recovered symbolic support against the ground-truth library terms after coefficient thresholding at 10−3 . A useful diagnostic is the exploratory cosine |⟨fˆaug , fˆphy ⟩D |/(∥fˆaug ∥D ∥fˆphy ∥D ) between the learned P augmentation and the assembled symbolic prediction fˆphy = j wj ϕj . This quantity is an indirect P proxy for the OrthoReg training penalty j ⟨fˆaug , ϕj ⟩2D (the two coincide up to a library reweighting only when w concentrates on a single term) and is numerically unstable when either component has small empirical norm, so we do not report it as a result and do not use it to rank methods. G.2
In-distribution and out-of-distribution evaluation
We distinguish four evaluation regimes: • In-distribution (ID): test trajectories drawn from the same distribution of initial conditions and time points as training; • In-distribution, time-extrapolated (ID,ext): same initial-condition distribution as training, evaluated past the training horizon T on the autonomous continuation of the trajectory; this stresses temporal extrapolation while keeping the IC distribution fixed; • OOD-T2 (initial-condition extrapolation, “IC-Ext”). trajectories initialized outside the training range of initial conditions; • OOD-T3 (parametric extrapolation, “Param-Ext”). trajectories generated from perturbed system parameters while keeping the symbolic library fixed, probing extrapolation under model mismatch. Tables retain the OOD-T2/OOD-T3 column labels for compactness; the regime/basin axis corresponds to the Duffing cross-basin split (Section F.2). Hyperparameters, including λ, are selected on a validation split disjoint from these test splits and fixed across systems. G.3
Data generation and derivative targets
Trajectories are generated by numerically integrating the ground-truth vector field using scipy.integrate.odeint (LSODA adaptive solver) and then sampled to obtain state observations {x(ti )}. Models are trained using the vector-field regression objective (Equation (4)) with derivative targets yi ≈ ẋ(ti ) computed via smoothed finite differences (PySINDy SmoothedFiniteDifference, second-order accuracy with Savitzky–Golay smoothing, window length 5). 25
Temporal discretization. For the pendulum experiments we simulate trajectories over a horizon T = 6.0 with nT = 100 time points, yielding a step size ∆t = T /(nT − 1) ≈ 0.0606. Other systems use their respective simulation horizons and discretizations as specified by the experimental configuration. G.4
Model classes and symbolic library
Neural augmentation. Hybrid models use a multilayer perceptron with three hidden layers and tanh activations. The hidden width is 128 for hybrid models in the main experiments. The pureneural baselines (PINN and Universal ODE) reported in Table 2 are trained at the same hidden width 128, so the neural-capacity comparison is at parity. Weights are initialized using standard variance-preserving schemes as implemented in the respective models. Symbolic library. We use a fixed feature library combining a polynomial basis (degree 2, with cross-terms for multivariate systems) and a Fourier basis (nfrequencies = 1, i.e. sin and cos at the unit frequency). Throughout, M counts scalar feature–component pairs (basis features times state dimension n), so for the pendulum (n = 2) with poly-2 ⊕ Fourier-1 we have M = 20, and for SIR (n = 3) the same construction gives M = 45. Symbolic coefficients are optimized with an L1 penalty and thresholded at 10−3 when computing support recovery metrics. G.5
Optimization and effective batch size
All models are optimized using Adam. The derivative regression phase is trained for 2000 epochs with learning rate 0.0089. In the default Monte Carlo setup we use nsamples = 1000 trajectories, all of which fit into a single optimization batch; thus, each epoch effectively processes the full training set. When minibatching is used, the batch size is stated explicitly. G.6
Baselines
Pure symbolic (SINDy-like). Sparse regression over the symbolic library with an L1 penalty, trained using the same optimizer infrastructure as the hybrid methods to ensure comparability. The canonical SINDy baseline uses sequential thresholded least squares (STLSQ). L2 -regularized hybrid. Joint optimization over symbolic coefficients w and neural parameters ϑ with loss L = Lfit (w, ϑ) + λ1 ∥w∥1 + λ2 ∥fˆaug ∥2D . OrthoReg.
Identical to the L2 hybrid baseline except that the neural regularizer is replaced by λ
M X
⟨fˆaug , ϕj ⟩2D .
j=1
PINN baseline. A physics-informed neural network that encodes known components of the pendulum dynamics as a soft constraint, while leaving unmodeled residual dynamics unconstrained. Universal ODE. A fully neural model without a symbolic component; uses the same depth and activations as the hybrid networks. Canonical PySINDy STLSQ (selection–refit baseline). The Pure Symbolic rows in the main tables (Tables 1, 3, 6 and 7) use the same continuous L1 +Adam sparse-regression pipeline as the hybrid models (Section G.6), which gives a controlled comparison within one optimization stack. As a canonical SINDy-style selection–refit reference, we additionally run PySINDy STLSQ on the same symbolic libraries and train/ID/OOD splits: SequentialThresholdedLeastSquares(threshold=0.045, normalize_columns=False) with second-order finite differences (FiniteDifference(order=2)). For pendulum ẋOOD,T3 under STLSQ, we fix the support selected on train and refit the active coefficients by ordinary least 26
Table 8: PySINDy STLSQ versus Pure (Adam L1 ) and OrthoReg (Tables 1, 3, 6 and 7). †: divergent rollout or clamp. Metric
STLSQ
Pure
OrthoReg
Pendulum (Equation (9)) 1.17 ±0.00 1.27 ±0.07 1.89 ±0.01 † 1.02 ±0.07 1.01 ±0.01 0.70 ±0.00 0.43 ±0.09 0.93 ±0.15 9.0 ±0.0 11.4 ±2.9 3.6 ±1.3
MSE(ẋID , e−2) MSE(xOOD,T2 , e0) F1 #Terms
Duffing (Section F.2) † †
4.75 ±0.13 5.59 ±0.09 11.03 ±6.09 4.93 ±0.35 0.62 ±0.00 0.40 ±0.05 0.63 ±0.05 7.0 ±0.0 11.0 ±1.6 3.4 ±0.5
MSE(xID , e−1) MSE(xOOD,T2 , e0) F1 #Terms
Lotka–Volterra (Section F.1) MSE(ẋOOD,T2 , e−2) 41.8 ±0.0 1.21 ±0.00 1.22 ±0.00 F1 0.30 ±0.00 0.22 ±0.01 0.24 ±0.01 SIR (Section F.3) MSE(ẋOOD,T2 , e−4) 233 ±0.00 2.40 ±0.11 16.54 ±0.80 F1 0.18 ±0.00 0.14 ±0.03 0.20 ±0.06
squares on the OOD-T3 trajectories before scoring. Where trajectory rollouts diverged or hit the numerical clamp, state MSE is marked †; STLSQ standard deviations are omitted when identical across all five seeds (deterministic splits and bitwise-identical fits). The comparison underscores why OrthoReg targets partial library misspecification rather than purely symbolic reconstruction. STLSQ reduces shrinkage bias and improves support recovery relative to the continuous L1 +Adam Pure baseline on the pendulum and Duffing systems. However, because STLSQ has no residual component, out-of-library dynamics must be approximated within the fixed dictionary, which can lead to unstable rollouts even when derivative-fit or F1 diagnostics look reasonable. OrthoReg instead couples sparse symbolic recovery with a flexible residual discouraged from using library directions, yielding stronger pendulum recovery and finite hybrid rollouts in the regimes above. The LV and SIR results show the complementary limitation: when the library mismatch is severe or the evaluation is dominated by time-dependent perturbations, F1 alone is not a reliable summary of dynamical performance. G.7
Per-system reproducibility values
For each system we list (i) the ground-truth parameters used for training trajectories, (ii) the perturbation used to generate the OOD-T3 split, and (iii) the simulation horizon and step size. Final regularization hyperparameters are shared across systems and listed at the end of the block; the same values are used for both the L2 baseline and OrthoReg, applied to the corresponding regularization term. • Modified damped pendulum (Equation (9), dataset theoretical_pendulum.yaml): – parameters ω0 = 1.0, α = 0.2, (β1 , β2 , β3 ) = (0.3, 0.25, 0.15); – OOD-T3 perturbation: all parameters scaled by 1.2; – T = 6.0, nT = 100, ∆t ≈ 0.0606. • Duffing oscillator (Section F.2, dataset duffing.yaml): – parameters a = −0.5, b = −1.0, c = 0.1 (cubic term cx3 absent from the dictionary); – OOD-T3 perturbation: (a, b, c) 7→ (1.2 a, 1.2 b, 2.0 c); – T = 40, ∆t = 0.01. • Lotka–Volterra (Section F.1, dataset complex_orthogonal_lv.yaml): – base parameters α, β, γ, δ as in lv.py (Blasius-style scaled rates 0.1·12, 0.005·12, 0.04· 12, 0.00004·12); missing-dynamics weights (ε1 , ε2 ) = (3.0, 2.0) at the reported difficulty 2.0; 27
– OOD-T3 perturbation: parameter ranges (0.2, 0.3) · 12, (0.01, 0.015) · 12, (0.08, 0.12) · 12, (0.00008, 0.00012) · 12 for (α, β, γ, δ) respectively, with initial conditions (N0 , P0 ) ∈ (100, 200) × (10, 20); – T = 10, ∆t = 0.1. • Tweaked SIR (Section F.3, dataset tweaked_sir.yaml): – base parameters β = 2.0, γ = 1.0; tweaked-term parameters (ε1 , ε2 , ε3 ) = (0.3, 0.25, 0.2), ωfast = 3.0, ωslow = 1.5, τmemory = 2.0; – OOD-T3 perturbation: β ∈ (8, 12), γ ∈ (0.8, 1.2) (parameter rescaling applied to the standard SIR backbone); – T = 10, ∆t = 0.1. • Final regularization hyperparameters (selected per Section G.8 and shared across systems): – L2 baseline: λ2 = 0.005 (on ∥fˆaug ∥2 ), symbolic sparsity µ = 0.001; D
– OrthoReg: λ = 0.005 (on the inner-product penalty), symbolic sparsity µ = 0.003. G.8
Regularization selection
We select regularisation strengths (λ, µ) for each method by a small validation sweep, using prediction loss on a held-out in-distribution trajectory split as the selection criterion. Symbolic recovery is tracked over the same sweep but reported only as a diagnostic. The selected (λ, µ) are fixed across in-distribution and out-of-distribution evaluations; the same protocol is applied to the L2 baseline with its own grid. P Penalty scale dependence on library norms. The orthogonality penalty λ j ⟨fˆaug , ϕj ⟩2D is scaledependent in the library: rescaling ϕj → cj ϕj rescales the per-feature contribution by c2j . The reported runs operate on bounded state ranges where the empirical norms ∥ϕj ∥D are comparable across j; for libraries with very different per-feature scales, standardising features to unit empirical norm (and rescaling λ accordingly) is a robust drop-in.
H
LLM Usage Disclosure
Large language models were used as writing and programming assistants throughout the project, including for text polishing, code development, and improving the presentation of mathematical arguments. The research idea, theoretical results, experimental design, analyses, figures, and manuscript were produced by the authors.
I
Limitations
Scope of the guarantee. Our guarantees apply to the additive decomposition f = fphy +faug under the empirical inner product used for training. They therefore do not directly cover compositional architectures such as fˆ = Tw + Nϑ ◦ Tw , unconstrained residual-learning pipelines, or population correlations of the learned augmentation without an additional uniform-convergence argument (Theorem C.1 and Section C). A full population-level theory is left for future work. Empirical-to-population transfer. The penalty bound (Theorem 4.1) and the error decomposition (Theorem B.2) are stated for the empirical inner product ⟨·, ·⟩D . A bound on the population correlation ⟨fˆaug , ϕj ⟩ν for the learned augmentation cannot be obtained by single-function concentration alone (the function depends on D); it would require uniform convergence over the augmentation hypothesis class (e.g. Rademacher complexity or VC-style arguments). A uniform-convergence statement is left to future work. Symbolic-discovery regime. OrthoReg is formulated for a fixed, finite, differentiable library and a continuous sparsity relaxation optimised jointly with the neural augmentation. Selection–refit methods such as STLSQ, L0 search, genetic programming, and transformer-based symbolic search operate in a different regime: after support selection, coefficients are typically refit by unregularised least squares, 28
avoiding shrinkage bias on the retained terms. Combining such solvers with orthogonality-regularised residual learning is a natural extension, but it is not covered by the joint-objective certificate analysed here. Dictionary-dependent interpretability and OOD scope. The recovered symbolic support is mechanistically meaningful only when the chosen library spans a scientifically relevant component of the true vector field. Theorem B.2 separates in-library and orthogonal-complement error under the empirical inner product; consequently, OrthoReg is most useful when missing dynamics are well separated from the library span, and less beneficial when the residual is strongly correlated with library functions on the observed distribution. OrthoReg therefore provides a complementary empirical decomposition, not recovery of governing laws absent from the library. Derivative supervision and noise. The experiments use vector-field regression with derivative targets estimated from smoothed finite differences. This isolates the decomposition problem but makes high-noise settings sensitive to derivative-estimation error, as seen in Section E.4; in particular, at ε ≳ 0.1 the smoothed-FD targets (Savitzky–Golay window of length 5) themselves limit recoverability. State-prediction losses (Equation (5)) or noise-aware differentiation are natural extensions for sensorlevel data. Implementation. The reported implementation optimises the continuous L1 -regularised objective with Adam and obtains exact sparsity by post-hoc thresholding. The orthogonality penalty is scaledependent in the library features, so feature standardisation is advisable for libraries with very different per-feature scales (Section G.8). The Pure-Symbolic rows use the same continuous sparsity pipeline as the hybrid models for controlled comparison; STLSQ-style pure-symbolic baselines are discussed in Section G.6.
29