Simultaneous model discovery and state estimation under high data corruption

TL;DR

The paper addresses learning ODEs from highly corrupted data by jointly inferring state trajectories and sparse models using a hybrid likelihood-based regression framework. It introduces a sparse, library-based regression with a smooth $L_0$-like penalty and BIC-based model selection, solved with a second-order Levenberg–Marquardt method that exploits Hessian sparsity via graph coloring and automatic differentiation. Extensive experiments on the Van der Pol, Lorenz, Lorenz-96, and Colpitts systems demonstrate robustness to noise and data gaps and competitive performance against WSINDy. The approach is general to ODEs and potentially extensible to DAEs and PDEs, offering interpretable, data-driven models with limited tuning requirements.

Abstract

This paper proposes a sparse regression strategy for discovery of ordinary differential equations from incomplete and noisy data. Inference is performed over both equation parameters and state variables using a statistically motivated likelihood function. Sparsity is enforced by a selection algorithm which iteratively removes terms and compares models using statistical information criteria. Large scale optimization is performed using a second-order variant of the Levenberg-Marquardt method, where the gradient and Hessian are computed via automatic differentiation. The proposed method is illustrated and tested on several systems with varying levels of noisy and incomplete data. Comparisons are made to a state-of-the-art algorithm for system identification, demonstrating competitiveness of the proposed approach.

Figures (15)

• Figure 1: Visualization of subset selection methodology to hone in on the correct terms under two different scenarios with parameter drop count $k = 5$. Iteration numbers are shown above each point, while both cases show the first model is rejected in iteration 8. Left shows acceptance of small forward steps resulting in no need to go to more complex models. Right shows a small forward step is rejected, resulting in the need for backward steps or testing more complex models.
• Figure 1: Comparison of CPU time (seconds) between the proposed LM method and other standard optimization algorithms.
• Figure 1: Comparison of our algorithm vs. WSINDy for the Van der Pol system, for monomial libraries with maximum orders $p=3,6$ with $18$ and $54$ terms, respectively. Each boxplot summarizes 50 independent noise realizations per noise level. A line at TPR = 1 indicates that most realizations achieve perfect recovery for our method.
• Figure 1: State estimation using recovered models for the Van der Pol oscillator using an insufficient library of polynomial terms (up to degree 2).
• Figure 2: Validation error and relative $\ell_2$ error in the recovered coefficients for the Lorenz system under 10% noise (left) and 30% noise (right).