💻 computer science

Simultaneous model discovery and state estimation under high data corruption

This paper proposes a sparse regression strategy that simultaneously discovers ordinary differential equations and estimates state variables from incomplete, noisy data by leveraging a statistically motivated likelihood function, iterative model selection, and second-order optimization via automatic differentiation.

Original authors: Teddy Meissner, Karl Glasner

Published 2026-02-18

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Teddy Meissner, Karl Glasner

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are a detective trying to solve a mystery, but the crime scene is a mess. The clues (data) are scattered, some are missing entirely, and the ones you do have are covered in mud and scribbles (noise). Your goal? To figure out the exact set of rules (the laws of physics) that govern how the suspects moved, even though you can't see them clearly.

This paper introduces a new, super-smart detective tool designed to solve exactly this kind of problem. Here is how it works, broken down into simple concepts:

1. The Problem: The "Messy Crime Scene"

In science, we often try to write down equations that describe how things change over time (like how a virus spreads or how a pendulum swings).

The Old Way: Traditional methods usually assume the data is perfect. If the data is noisy or missing chunks, these methods get confused and give up. It's like trying to solve a puzzle when half the pieces are missing and the rest are glued together.
The Challenge: Real-world data is rarely perfect. It's often incomplete (we couldn't measure everything) and noisy (our sensors make mistakes).

2. The Solution: The "Dual-Task Detective"

The authors propose a method that does two things at the same time, like a detective who is both reconstructing the crime and writing the rulebook simultaneously.

Task A: State Estimation (Reconstructing the Crime): Instead of just looking at the muddy clues, the method tries to guess what the true path of the suspect actually was. It fills in the missing gaps and smooths out the scribbles, creating a clean, smooth story of what happened.
Task B: Model Discovery (Writing the Rulebook): While it's cleaning up the story, it's also trying to figure out the mathematical laws that caused that story. It asks, "What simple equation could produce this smooth path?"

3. The Secret Sauce: "Sparsity" (The Occam's Razor Principle)

When trying to find the rulebook, there are infinite possibilities. You could write a 100-page equation that fits the data perfectly, but it would be nonsense.

The Analogy: Imagine you are trying to guess the ingredients of a soup. You could say it's made of 50 different spices, but the chef only used salt, pepper, and garlic.
The Method: This tool uses a "sparsity" strategy. It aggressively hunts for the simplest explanation. It tries to remove unnecessary ingredients (terms in the equation) until only the essential ones remain. It uses a statistical "scorecard" (called BIC) to decide: "Is this extra ingredient actually helping, or is it just clutter?"

4. How It Handles the Mess

The paper highlights three specific ways this detective is better than the old ones:

Noise: Instead of trying to filter out the noise (which can sometimes remove good clues), the method infers the truth. It realizes, "Okay, this data point looks weird, but if I adjust my guess of the true path slightly, everything makes sense."
Missing Data: If you didn't measure the suspect's speed for 10 minutes, the method doesn't panic. It uses the laws of physics it's discovering to "fill in the blanks" logically.
Complexity: It can handle equations that aren't just simple straight lines. It can find complex, curved relationships (like exponential growth) that other methods miss.

5. The Engine: The "Super-Optimizer"

To do all this math, the paper uses a powerful engine called Automatic Differentiation combined with a Levenberg-Marquardt algorithm.

The Analogy: Imagine trying to find the lowest point in a foggy mountain valley.
- Old methods take small, blind steps (like feeling the ground with a stick).
- This method has a helicopter that can see the whole shape of the valley and the slope of the ground instantly. It knows exactly which way to jump to get to the bottom (the best solution) much faster.

6. The Results: A Winning Detective

The authors tested this method on famous chaotic systems (like the Lorenz system, which models weather, and the Van der Pol oscillator, which models electrical circuits).

The Test: They took perfect data, added huge amounts of noise (up to 50% error!), and removed chunks of the timeline.
The Outcome: Their method successfully reconstructed the original equations and the true paths of the systems, even when the data was terrible. It beat the current "gold standard" (a method called WSINDy) in almost every test.

Summary

Think of this paper as a new, high-tech toolkit for scientists. It allows them to look at a messy, broken, and incomplete set of observations and say, "I can tell you exactly what the underlying laws of nature are, and I can tell you exactly what happened, even though you gave me terrible data." It's about finding the simple, elegant truth hidden inside a chaotic mess.

1. Problem Statement

The paper addresses the challenge of model discovery (identifying the governing Ordinary Differential Equations, or ODEs) and state estimation (reconstructing the true system trajectory) simultaneously from data that is both noisy and incomplete.

Context: Traditional methods often struggle when data is corrupted.
- Residual Regression (e.g., SINDy): Highly sensitive to noise because it requires differentiating noisy data to compute residuals.
- Parameter Estimation (Nonlinear Least Squares): Often ignores model uncertainty, leading to non-convex optimization problems that are difficult to solve.
- Incomplete Data: Many real-world scenarios have missing time steps or unmeasured state variables, which breaks methods relying on dense differentiation.
Goal: To develop a robust framework that infers sparse ODE parameters ( $\theta$ ) and the underlying state variables ( $u$ ) simultaneously, even when the data is sparse and corrupted by significant noise.

2. Methodology

The authors propose a hybrid sparse regression strategy grounded in statistical likelihood and second-order optimization.

A. Statistical Formulation (Likelihood Function)

Instead of minimizing the residual of the data directly, the method maximizes a likelihood function that accounts for uncertainty in both the data and the model.

The objective function $L_\lambda(\theta, u)$ $L_{λ} (θ, u)$ is defined as:
$L_\lambda(\theta, u) = \|N(u; \theta)\|_2^2 + \lambda \| \hat{u} - u|_D \|_2^2$
Where:
- $N(u; \theta) = 0$ represents the discretized ODE system.
- $\hat{u}$ is the noisy, incomplete observed data.
- $u|_D$ is the model state restricted to the observed data points.
- $\lambda$ is a ratio of variances balancing model fidelity vs. data fidelity.
Interpretation: This formulation interpolates between nonlinear least squares (fitting data) and residual minimization (fitting the equation). By optimizing over both $u$ and $\theta$ , the method effectively "denoises" the data by finding a smooth trajectory that satisfies the ODE constraints.

B. Sparsity Enforcement via Bayesian Information Criteria (BIC)

To discover the correct equation structure (i.e., which terms exist), the method enforces sparsity without using hard thresholds or $L_1$ penalties (which can introduce bias).

Smooth Approximation: A smooth approximation of the $L_0$ norm is used to allow gradient-based optimization.
Selection Algorithm: The authors employ an Adaptive Pruning algorithm:
1. Start with a full library of candidate terms.
2. Iteratively remove terms based on a Bayesian Information Criterion (BIC) score, which balances goodness-of-fit against model complexity.
3. Forward/Backward Refinement: The algorithm adaptively adjusts the number of terms removed. If a step is rejected, it switches to a "backward" mode to add terms back, ensuring the search does not get stuck in local minima of the model space.

C. Optimization Strategy

The optimization problem is high-dimensional and non-convex.

Second-Order Method: The authors utilize a Levenberg-Marquardt (LM) variant using the exact Hessian.
Efficiency: To handle the computational cost of the Hessian:
- Automatic Differentiation (AD): Used to compute gradients and Hessians.
- Graph Coloring: Exploits the sparsity structure of the Hessian (derived from the time-discretization stencil) to minimize the number of Hessian-Vector Products (HVPs) required.
Result: This approach converges significantly faster (approx. 10x) than first-order methods (like L-BFGS or Adam) and other second-order methods tested.

3. Key Contributions

Hybrid Regression Framework: A unified likelihood-based approach that simultaneously infers parameters and states, making it robust to high noise and missing data where residual-based methods fail.
BIC-Driven Sparsity: A novel selection algorithm that uses BIC rather than fixed thresholds or $L_1$ regularization, allowing for a more statistically rigorous model selection process.
Second-Order Optimization with AD: Demonstrates the practicality of using exact second-order derivatives (Hessian) for system identification by leveraging automatic differentiation and graph coloring to handle sparsity efficiently.
General Applicability: The formulation is not restricted to equations linear in parameters, allowing for the discovery of models with nonlinear parameter dependencies (e.g., exponential terms with unknown exponents).

4. Results

The method was tested on several benchmark systems with varying levels of noise (up to 50%) and missing data (up to 50% of points removed).

Van der Pol Oscillator:
- Outperformed the state-of-the-art WSINDy algorithm across all noise levels.
- Successfully recovered the system structure and state even with 50% noise and continuous gaps in data.
Lorenz System (Chaotic):
- Demonstrated the ability to recover the "ground truth" chaotic trajectory despite significant noise, a task where standard integration of discovered models usually fails due to sensitivity to initial conditions.
- Achieved high True Positive Rates (TPR > 0.9) for model recovery even at 40% noise.
Lorenz-96 (High Dimensional):
- Successfully scaled to a 5-variable chaotic system, maintaining high accuracy in parameter estimation and state reconstruction.
Colpitts Oscillator (Nonlinear Parameters):
- Successfully identified a system containing an exponential term with an unknown parameter in the exponent ( $e^{ax}$ ), a scenario where linear-in-parameter methods (like standard SINDy) cannot be directly applied.
Comparison: In all cases, the proposed method showed lower parameter error ( $RE(\theta)$ ) and higher model recovery rates compared to WSINDy.

5. Significance

This paper represents a significant advancement in data-driven scientific discovery:

Robustness: It solves the "noise sensitivity" problem that has plagued many sparse regression techniques (like SINDy) by treating the state as a latent variable to be inferred rather than a fixed input.
Handling Incompleteness: It provides a viable solution for systems where sensors fail or data is sparse, a common issue in biological and physical experiments.
Computational Efficiency: By proving that second-order optimization with automatic differentiation is feasible and superior for this class of problems, it opens the door for more accurate and faster model discovery.
Generality: The framework moves beyond the restrictive assumption that models must be linear in parameters, broadening the scope of discoverable physical laws.

In summary, Meissner and Glasner present a statistically motivated, computationally efficient, and highly robust framework for discovering differential equations from imperfect real-world data, outperforming current benchmarks in noisy and sparse regimes.