Sparse Orthogonal Regression Technique: A Spectral Framework for Equation Discovery, Approximation, and Integration

AI summary

Overview

Research area: Machine learning for scientific computing, specifically data-driven equation discovery, numerical integration, and function approximation using sparse regression in orthonormal bases.

Technical level: Intermediate. The paper assumes familiarity with orthonormal bases, linear regression, L1 regularization, and ordinary differential equations, but it is written as a workflow-level contribution rather than a new optimization theory.

Scope: The paper introduces the Sparse Orthogonal Regression Technique (SORT), a framework that learns sparse expansions in prescribed orthonormal bases from noisy, irregularly sampled data, and reuses the resulting coefficients for vector-field reconstruction, integral estimation, and order-consistent approximation.

What This Paper Is About

Scientific computing and machine learning both need to reconstruct an unknown function from finite, noisy, often irregularly sampled observations and then reuse that representation downstream. SORT addresses this by fitting a sparse coefficient vector in a fixed orthonormal basis using L1-regularized regression, so coefficients are estimated directly from data without quadrature or analytic inner products. The central application is discovering governing ordinary differential equations, where each vector-field component becomes a sparse spectral expansion rather than a symbolic formula or a selection from a finite candidate library.

Key Contributions

1. A reusable sparse spectral workflow for dynamical systems. SORT represents each vector-field component as a sparse expansion in a chosen orthonormal basis (Legendre or Fourier in the experiments), estimated by L1-regularized regression, positioning it as a complementary route to symbolic regression, grammar-based discovery, and SINDy-style sparse identification.

2. Integration as coefficient readout. The paper shows that if a basis direction is chosen to represent the integration functional, the integral is recovered from a single estimated coefficient up to a known normalization factor. The construction uses the normalized Riesz representer of the integration functional and a factorization trick (the term D(t) = product over j of (1 + t_j²)) that gives a closed-form norm.

3. Order-consistent model growth as a practical diagnostic. Because increasing truncation order in an orthonormal basis refines rather than reparameterizes the model, low-order coefficients retain their meaning across model sizes, allowing models to be compared by active coefficient sets, coefficient persistence, sparsity level, and predictive accuracy.

4. An explicit reframing of basis design as a modeling decision. The paper demonstrates through controlled comparisons that when a finite library misses the problem-specific nonlinearity, library-based models fail, and argues that basis choice belongs to the scientific modeling problem rather than to implementation detail.

Main Findings

• Competitive accuracy with more stable degradation under coarse sampling. In the Dynobench-derived experiment, both SORT and SINDy behave similarly at fine sampling with reliable derivatives. As the sampling interval increases (panel (b) illustrates Δt_sample = 0.75), Lotka–Volterra and Van der Pol show abrupt SINDy rollout failure while the orthogonal sparse model remains bounded and degrades more gradually. The paper states the main effect is not better fine-sampling accuracy but slower, less catastrophic degradation.

• Neutral and well-matched cases exist. Stuart–Landau is reported as a neutral case, and Bar Magnets tests the Fourier setting. On the Thomas attractor, both methods perform comparably and trigonometric features improve accuracy for both.

• Robustness under representation mismatch. On a three-dimensional cyclic system where coupling uses the Bessel function J₁(·), which is not sparse in low-degree polynomial or trigonometric libraries, the SINDy-style baseline degrades as the training fraction decreases or noise increases, while SORT remains robust. Panels in Figure 2 vary training fraction at fixed noise σ = 0.1 and vary noise at fixed training fraction of 20%.

• All comparisons are representation-controlled. SORT and the SINDy baseline share the same sampled states, derivative estimates, noise model, subsampling strategy, rollout horizon, and evaluation metric. In the noise experiments, derivatives are estimated by applying a Savitzky–Golay filter to each observed coordinate.

• Integrals track closed-form values in one dimension. SORT accurately tracks closed-form values for Gaussian integrals ∫₀^L e^(−αx²) dx and Fresnel-type integrals ∫₀^L cos(αx²) dx across a wide range of α, including strongly oscillatory Fresnel cases where accurate quadrature would require resolving rapid oscillations.

• Integration accuracy declines gradually with dimension. Higher-dimensional separable integrals remain accurate in moderate dimensions, with error increasing with dimension, reflecting growth of the basis and the sample complexity of sparse coefficient recovery.

• SORT suppresses unstable high-frequency modes. For f(x) = x on [−π, π], which has Fourier sine coefficients |b_k| = 2/k, OLS fits unstable high-frequency coefficients under noise and subsampling and produces oscillatory test error, while SORT preserves the dominant low-order spectral structure. This used 20% of samples for training and 80% for testing.

• Lowest test error and stable coefficient hierarchy in 2D approximation. For f(x₁, x₂) = exp(x₁ cos(2x₂)) on [−1, 5]², using 20% noisy training data (σ = 0.1), 80% clean test data, tensor-product Legendre bases with total-degree cutoff D_max = 50, and parameter counts k ∈ {60, 80, 100}, SORT achieves the lowest test error among the compared methods. OLS in a monomial basis changes substantially with order, OLS in a Legendre basis benefits from orthogonality but remains noise-sensitive, and RBF/RFF ridge models produce dense weights without an interpretable hierarchy.

• Systems excluded and why. The coupled phase oscillator is excluded as non-autonomous; bacterial respiration, predator–prey, glider, and shear flow contain rational, singular, or ratio-trigonometric terms such as 1/x or cot(y); Lorenz is excluded because its standard polynomial form gives a polynomial SINDy library a representational advantage.

• Theoretical background cited, not newly derived. The paper notes the bound ||ĉ − c||_∞ ≤ 2λ in the ideal empirically orthonormal case, the Fourier coefficient decay |f̂_k| ≤ V/(2π|k|), and the K-mode Fourier approximation error bound ||f − S_K f||_L² ≤ V/√(πK), presenting these as background intuition rather than new results.

• Computational scaling stated but not benchmarked. Constructing the design matrix costs O(NK), and coordinate-descent LASSO solvers scale with repeated passes over this matrix. The paper states SORT does not remove the curse of dimensionality and that severe basis mismatch can degrade performance. No explicit runtime measurements, dataset sizes, or numeric RMSE values are reported in the provided content.

Methodology in Plain English

The approach rests on a simple observation: if you pick a set of basis functions that are orthonormal with respect to some domain and measure, then any function in that space is a weighted sum of those basis functions, and the weights (coefficients) are inner products. Classical spectral methods compute those weights by quadrature or exact projection. SORT instead treats the weights as unknowns and fits them by regression on the sample values, adding an L1 penalty that drives weak or noise-supported coefficients to zero.

The regularization strength λ is selected empirically, for example by validation error or cross-validation, with a preference for the sparsest stable model. To turn this into an equation-discovery method, the researchers take observed trajectories, estimate derivatives (by finite differences in one experiment, by Savitzky–Golay smoothing and differentiation in the noisy experiment), and regress those derivative values onto the basis functions evaluated at the corresponding states. Each vector-field component gets its own sparse coefficient vector.

For integration, the trick is to transform the integrand so that the integration functional itself becomes a basis direction: the normalized constant basis element is set equal to D(t) divided by its L² norm, so the integral reduces to reading off the zero-th coefficient and multiplying by (b−a)/2 raised to the dimension d and by the norm of D. For higher dimensions, sparse index sets such as hyperbolic-cross truncations limit basis size.

Across all experiments, SORT and the baselines are given identical data, derivative estimates, noise, subsampling, and evaluation metrics, so any difference is attributable to the representation rather than the preprocessing.

Why This Matters

Impact on research. The paper argues that sparse regression on orthogonal features is not algorithmically new, and does not claim otherwise. Its contribution is representational: it positions a learned sparse orthonormal expansion as an intermediate object that can be reused across equation discovery, integration, and approximation, and that can precede a symbolic search stage. This reframes basis selection as an explicit part of scientific modeling, comparing it to the use of Bessel functions for problems with radial or cylindrical symmetry. It also supplies an order-consistent diagnostic that most machine-learning workflows lack, since increasing model complexity there usually changes the whole parameterization rather than refining an existing structure.

Real-world applications (as described in the paper):

• Surrogate modelling for expensive simulations, where a learned sparse expansion stands in for costly simulation runs.
• Numerical integration from sampled observations, useful when explicit quadrature is difficult, high-dimensional, or unavailable.
• Data-driven discovery of governing equations from time-series data, producing a sparse spectral vector field as a reusable model.
• Approximation with interpretable growth, where increasing truncation order refines accuracy without invalidating previously learned low-order structure.

Industry relevance. The framework targets settings where derivative information is degraded or where the natural basis is uncertain, which matches practical sensor-driven and simulation-derived data. Sparse Polynomial Chaos Expansion is identified as methodologically close, and that field is heavily used for uncertainty quantification and sensitivity analysis in engineering. The paper also reports that code for reproducing experiments and figures is available at a Zenodo DOI (10.5281/zenodo.21707070), and the work is supported by the European Union's Horizon Europe Marie Skłodowska-Curie programme under grant agreement No. 101081355, the Slovenian Research Agency Gravity project GC-0001, and research core funding No. P1-0188.

Future Directions

1. Adaptive and anisotropic basis selection, so that the orthonormal coordinate system can be tuned to the problem rather than prescribed in advance.

2. Empirical orthogonalization for irregular samples, addressing the fact that with irregular finite samples the empirical orthogonality condition (1/N)ΦᵀΦ = I_K is only approximate, so coefficient stability must currently be assessed numerically rather than guaranteed.

3. Larger-scale applications, extending the framework beyond the two- and three-dimensional systems and moderate-dimensional separable integrals tested here.

4. Coupling to a symbolic stage, where a recovered stable expansion is used to guide a subsequent search for simpler closed-form expressions consistent with it, and where a future extension might handle non-autonomous systems such as the coupled phase oscillator that the paper excluded.

Target Audience

Readers who will benefit most are researchers and practitioners in scientific machine learning and system identification who already work with SINDy-style sparse identification, sparse Polynomial Chaos Expansion, or symbolic regression and want a complementary representation-level alternative. It also suits numerical analysts interested in how classical spectral approximation behaves under data-driven coefficient estimation, and applied scientists who need integral estimates or compact surrogate models from noisy, irregularly sampled measurements. The paper is not aimed at readers seeking new sparse-recovery theory, since it explicitly does not propose any, nor at those who need runtime benchmarks or large-scale performance evidence, which are not reported.

Authors’ abstract

We develop the Sparse Orthogonal Regression Technique (SORT), a sparse spectral framework for learning orthonormal-basis expansions from noisy and irregularly sampled data. SORT estimates expansion coefficients directly from observations using L1-regularized regression, avoiding explicit quadrature or analytic inner-product evaluation. The central application is data-driven discovery of ordinary differential equations: vector fields are represented in chosen orthogonal bases and learned as sparse coefficient expansions. This provides a complementary route to symbolic regression, grammar-based discovery, and SINDy-style sparse identification by first recovering a compact spectral representation, which can later guide searches for simpler analytic forms. Across the dynamical-system experiments, SORT matches or improves upon library-based sparse-regression baselines when the basis is well adapted to the problem, and shows more stable degradation under sparse sampling, noisy derivative estimates, and representation mismatch. Specific examples illustrate why this representation is useful: if a finite library misses the problem-specific nonlinearity, the resulting model can fail. SORT is not immune to mismatch, but it shifts the problem away from brittle selection among generic terms to basis design adapted to the problem domain. The experiments also show that dominant low-order coefficients persist as model order increases, supporting order-consistent model growth. Beyond equation discovery, the same learned expansion supports nonlinear approximation and estimation of complex, high-dimensional integrals by coefficient readout. Overall, SORT provides a reusable intermediate representation for system identification, approximation, and integration, while making basis design an explicit part of the scientific modeling problem.