A Block Coordinate Descent Proximal Method for Simultaneous Filtering and Parameter Estimation

🎯 Pitch

BCD-prox beats state-of-the-art ODE parameter estimation and filtering methods even when they fail completely—under high noise and bad initial guesses where competitors produce errors above 100 in half of all trials. The algorithm jointly learns states and parameters via a single unified objective with just one hyperparameter, converging in seconds on 40-dimensional chaotic systems where Bayesian and spline-based approaches take minutes to hours or diverge entirely.

1. Executive Summary

This paper proposes and analyzes a block coordinate descent proximal algorithm (BCD-prox) for simultaneous filtering and parameter estimation of ordinary differential equation (ODE) models. Tested on nonlinear ODE systems with up to 40 dimensions—including Lotka–Volterra, FitzHugh–Nagumo, the Rössler attractor, and Lorenz-96—BCD-prox minimizes a unified objective function that directly measures fidelity (how well estimated states and parameters satisfy the discretized ODE) while using a proximal term to control step size (analogous to inverse learning rate in gradient descent), requiring only a single hyperparameter λ. Compared to spline-based iPDA, Bayesian, and extended Kalman filter methods, BCD-prox exhibits substantially increased robustness to noise, parameter initialization, and hyperparameter settings, reduced training times (seconds vs. minutes to hours on a standard laptop for 400-state problems), and improved accuracy of both filtered states and estimated parameters, establishing that the algorithm's distribution-free, batch optimization approach maintains performance under heavy noise and poor initialization conditions where competing methods diverge or produce errors exceeding 100 in roughly half of all experimental trials.

2. Context and Motivation

The Core Problem: Learning ODE States and Parameters from Noisy Data

The paper addresses a deceptively difficult problem: given noisy observations of a dynamical system at discrete time points, simultaneously recover both the clean underlying states (filtering) and the unknown parameters of the governing ordinary differential equations (parameter estimation). Mathematically, we start with a known form of the vector field $f(x, \theta)$ but unknown parameters $\theta \in \mathbb{R}^p$ , and we observe $y(t_i) = x(t_i) + z(t_i)$ at times $\{t_i\}_{i=1}^T$ , where $z(t_i)$ represents noise. The goal is to estimate both $X$ (the filtered states) and $\theta$ (the parameters) using only the noisy observations $Y$ .

This framing, presented in Section 1 and formalized in Section 2, is pervasive across science and engineering—any domain where models are expressed as differential equations and measurements are corrupted by noise faces this challenge. The authors specifically highlight biology as a domain where nonlinear ODE models and highly noisy observations commonly co-occur, creating a pressing need for methods that do not break down under realistic noise levels.

Why Existing Approaches Struggle

The paper identifies four major families of prior methods and diagnoses their specific failure modes. Understanding these weaknesses is essential because BCD-prox's design choices are direct responses to them.

Nonlinear least squares methods (Bard, 1973; Benson, 1979; Himmelblau et al., 1967; Hosten, 1979) iteratively update parameter estimates to bring model predictions close to measurements. The paper identifies a critical fragility: these methods diverge when the initial parameter guess is far from the true values. In real applications, we rarely have good initial guesses, making this a serious practical limitation. As Section 1 states:

"These methods diverge when the initial parameters are far from the true parameters."

Spline-based methods (Cao & Zhao, 2008; Cao et al., 2011; Poyton et al., 2006; Ramsay et al., 2007; Varah, 1982) form a significant category that includes the influential iPDA (iterated principal differential analysis) approach. These methods first fit smooth functions (typically cubic splines or other estimators like smoothing kernels and local polynomials) to the noisy observations, then use the differentiable splines to estimate parameters via regression. While conceptually elegant, the paper identifies multiple failure points:

"These methods are sensitive to numerous hyperparameters (such as smoothing parameters and the numbers/positions of knots), to parameter initialization, and to the magnitude/type of noise that contaminates the data."

The hyperparameter sensitivity is particularly problematic because, as the paper notes in its detailed comparison with iPDA in Section 3, the optimal smoothing parameter $\lambda$ depends on both the unknown noise $Z$ and the vector field $f$ —neither of which is accessible in practice. As the paper points out:

"In a real problem, we do not have access to the clean states $X$ (all we have are the noisy observations $Y$ ), so we cannot find the right $\lambda$ by cross-validation."

Furthermore, the paper's experiments (Figures 3 and 7) show that iPDA—the second-best method after BCD-prox—still produces prediction errors exceeding 100 in roughly half of all trials when initialized far from the true parameters. This is a devastating failure rate for practical use.

Bayesian approaches (Calderhead et al., 2009; Dondelinger et al., 2013; Girolami, 2008; Gorbach et al., 2017) impose prior distributions over states and parameters, then compute posterior estimates. The paper identifies three specific limitations. First, these methods require careful specification of hyperparameters (prior distributions, variances, kernel widths), and getting these wrong produces poor results. Second, they feature large training times—hours on problems that BCD-prox solves in seconds. Third, and most critically for the simultaneous estimation problem the paper targets:

"Another disadvantage of these methods, mentioned by Gorbach et al. (2017), is that they cannot simultaneously learn clean states and parameters."

Gorbach et al. (2017) attempts to overcome this limitation using variational inference, but the paper notes that this solution "is not applicable to all ODE"—it only works for differential equations with a specific structural form. The paper's experiments (Figures 3 and 4) confirm these theoretical concerns: the Bayesian method (denoted "Bayes" in experiments) performs poorly across all initialization conditions, with substantially larger errors than BCD-prox.

Kalman filter methods, specifically the extended Kalman filter (EKF) and ensemble Kalman filters, represent online filtering approaches that make Gaussian assumptions about the noise and state distributions. The paper identifies several fundamental limitations in Section 1:

"Many other well-known nonlinear ODE filtering methods, including extended and ensemble Kalman filters as well as particle filters, are online methods that make Gaussian assumptions. In contrast, BCD-prox is a distribution-free, batch method."

The experimental comparison with EKF (Figure 4, second and third rows; extended in the supplementary material Figure 9) reveals exactly where these limitations manifest. EKF performs comparably to BCD-prox only under ideal conditions: very long time series ( $T = 10,000$ ) with low noise ( $\sigma^2 = 0.1$ ). In more realistic settings—shorter time series or higher noise—BCD-prox significantly outperforms EKF. Moreover, because EKF is an online method that never revisits earlier states, it cannot update the initial state estimate. Since the initial state is noisy, this permanently limits EKF's prediction accuracy regardless of how well parameters are estimated. The paper quantifies this: when noise increases to $\sigma^2 = 1.5$ with $T = 10,000$ , EKF's average estimation error reaches 1.30 while BCD-prox achieves 0.06.

The Deeper Structural Insight: Multiple Global Minima

A crucial theoretical insight that motivates the paper's approach appears in Theorem 1 (Section 3). The fidelity objective $E(X, \theta)$ —which measures how well estimated states and parameters satisfy the discretized ODE—has an infinite number of global minima that achieve $E = 0$ . The proof is simple but profound: choose arbitrary values for $\theta$ and the initial condition $x(t_1)$ , then step forward in time using the discretization scheme. By construction, every term in the fidelity sum vanishes, giving $E = 0$ . Since $\theta$ and $x(t_1)$ are arbitrary, there are infinitely many such (state, parameter) pairs.

This has critical implications for algorithm design. Any method that minimizes the fidelity alone (without additional constraints) will produce arbitrary, non-unique solutions. The regularization term in spline-based methods ( $\lambda\|X - Y\|^2$ in iPDA's objective) attempts to resolve this ambiguity by keeping filtered states close to observations, but this creates a tension: strong regularization (large $\lambda$ ) keeps $X$ near $Y$ but may produce poor parameter estimates; weak regularization (small $\lambda$ ) allows $X$ to wander far from the data, potentially finding states that satisfy the ODE perfectly but have no relationship to the observed system.

How BCD-prox Addresses These Gaps

BCD-prox's approach differs fundamentally from all prior methods along several axes, as introduced in Section 1 and fully developed in Section 2.

Unified objective function. Rather than using separate objectives for filtering and parameter estimation (as most prior methods do), BCD-prox minimizes a single objective that combines the ODE fidelity with a proximal term. This unification means the algorithm simultaneously optimizes states and parameters toward a solution that satisfies the ODE while remaining tethered (via the proximal term) to the evolving estimate. The paper emphasizes:

"BCD-prox works by minimizing a unified objective function that directly measures how well the states and parameters satisfy the ODE system, in contrast to other methods that use separate objectives."

Direct state learning without smoothness assumptions. Spline-based methods learn states indirectly by fitting smooth functions to observations. This imposes artificial smoothness constraints that may not reflect the true system behavior and introduces knot-placement and smoothing-parameter decisions. BCD-prox learns the states directly in the original space:

"BCD-prox learns the states directly in the original space, instead of learning them indirectly by fitting a smoothed function to the observations."

The proximal term as adaptive step size control. The proximal term $\lambda\|X - X^{*(n-1)}\|^2$ in equation (5) plays a role that is easily misunderstood without careful reading. Section 3 explains that when viewed through the lens of proximal operators, this term approximates a gradient descent step:

$X^{*(n)} = X^{*(n-1)} - (2\lambda)^{-1}\nabla_X E(X^{*(n-1)}, \theta^{*(n)}) + o((2\lambda)^{-1})$

This reveals that $\lambda$ functions as an inverse step size. Large $\lambda$ means small steps, keeping the algorithm's progress slow but stable; small $\lambda$ allows large changes from one iteration to the next. Critically, unlike iPDA's regularization parameter (which ties states to the original noisy observations $Y$ forever), BCD-prox's proximal term ties states to the previous iteration's estimate. This means the algorithm can progressively move far from the noisy initial data $Y$ if the ODE fidelity term pulls it there—it is not permanently anchored to noisy observations:

"If the data $Y$ is heavily contaminated with noise, it may be wise to move far away from $Y$ as we iterate."

Robustness through simultaneous optimization. The block coordinate descent structure—alternating between optimizing $\theta$ with $X$ fixed, then $X$ with $\theta$ fixed—creates a virtuous cycle. Even if initial parameter estimates are poor, the state optimization step can improve the state estimates. These improved states then enable better parameter estimation in the next iteration. This mutual reinforcement explains why BCD-prox remains robust when other methods fail under poor initialization. The paper connects this to its theoretical motivation: the algorithm can step slowly away from poor initial parameter choices because both the BCD and proximal components moderate the optimization trajectory.

Distribution-free, batch processing. Unlike Kalman filter methods, BCD-prox makes no Gaussian assumptions and processes the entire time series as a batch, allowing information from all time points to influence parameter and state estimates simultaneously. This batch nature means the algorithm can update estimates based on global information rather than causally forward-only updates that cannot correct early mistakes.

Positioning Relative to Prior Work

The paper positions BCD-prox at the intersection of three active research areas: alternating minimization (Chatterji & Bartlett, 2017; Li et al., 2016; Yi et al., 2014), block coordinate descent (Xu & Yin, 2013; Zhang & Brand, 2017), and proximal methods (Parikh & Boyd, 2014; Sun et al., 2015). However, the authors note that before their work, no one had combined these techniques into a practical algorithm for the simultaneous filtering and parameter estimation problem. They explicitly state in Section 3:

"Though BCD-prox may seem straightforward, we cannot find prior work that utilizes precisely this approach."

The closest relative identified by the authors is iPDA (Poyton et al., 2006; Ramsay et al., 2007), and Section 3 provides a detailed conceptual comparison that forms the intellectual core of the paper's positioning. The key distinction is that iPDA's regularization anchors filtered states permanently to noisy observations, while BCD-prox's proximal anchoring is to the previous iterate—a moving target that allows escape from the initial noise. This shift from static to dynamic anchoring is subtle but, as the experiments demonstrate, has dramatic practical consequences for robustness under heavy noise.

3. Technical Approach

3.1 Reader Orientation

This paper presents an iterative optimization algorithm — not a learned model or a neural network — that takes as input noisy measurements of a dynamical system over time and produces as output both clean estimates of the system's true states and estimates of the unknown parameters governing the system's differential equations. The core idea is a block coordinate descent scheme applied to a unified objective function that combines an ODE fidelity term (measuring how well states and parameters satisfy the discretized equations) with a proximal anchor term (preventing the estimated states from changing too aggressively in any single iteration), where these two blocks — optimizing parameters given fixed states, then optimizing states given fixed parameters — are alternated repeatedly, with the proximal anchor updated at each outer iteration to point at the previous state estimate rather than at the original noisy observations.

3.2 Big-Picture Architecture (Diagram in Words)

The algorithm has four major components, arranged in a nested loop:

Multistep ODE Discretization Module: Takes the known vector field $f(x, \theta)$ and a choice of numerical integration order $m$ , and produces an algebraic constraint system that the estimated states and parameters must approximately satisfy. This module defines what it means for a (state, parameter) pair to be "ODE-consistent." It can use explicit Euler ( $m=1$ , first-order) or higher-order Adams-Bashforth methods ( $m \geq 2$ ).
Fidelity Objective $E(X, \theta)$ : A scalar function that measures the squared mismatch between the left-hand and right-hand sides of the discretized ODE system for a candidate state trajectory $X$ and parameter vector $\theta$ . This is the "physics-informed" component that enforces that estimated quantities obey the governing equations. The paper defines this in Equation (7) for general $m$ -step methods.
Proximal Anchor Term $\lambda\|X - X^{*(n-1)}\|^2$ : A quadratic penalty added to the fidelity objective during iteration $n$ that penalizes the new state estimate $X$ for straying too far from the previous iteration's state estimate $X^{*(n-1)}$ . The hyperparameter $\lambda$ controls the strength of this penalty and functions as an inverse step size in the equivalent gradient descent interpretation.
Block Coordinate Descent Loop: The outer iteration that alternates between two minimizations:
- Parameter block: $\theta^{*(n)} = \arg\min_\theta F_n(X^{*(n-1)}, \theta)$ — finding the best parameters given the current state estimate.
- State block: $X^{*(n)} = \arg\min_X F_n(X, \theta^{*(n)})$ — finding the best states given the just-updated parameters.

Information flows as follows: start with $X^{*(0)} = Y$ (the noisy observations serve as the initial state estimate) and an initial parameter guess $\theta^{*(0)}$ → construct the augmented objective $F_1(X, \theta) = E(X, \theta) + \lambda\|X - Y\|^2$ → optimize sequentially over $\theta$ then $X$ to get $(X^{*(1)}, \theta^{*(1)})$ → update the proximal anchor to point at $X^{*(1)}$ → repeat until the fidelity $E$ changes by less than $10^{-8}$ from one iteration to the next → once converged, compute final predicted states $\hat{X}$ by numerically integrating the ODE forward from the estimated initial condition using the estimated parameters, guaranteeing $E(\hat{X}, \theta^*) = 0$ exactly by construction.

3.3 Roadmap for the Deep Dive

First, I explain the ODE discretization machinery (Equation 6 for general multistep methods, with Equation 3 as the Euler special case), because every subsequent component depends on this discretized algebraic form of the differential equations.
Second, I build the fidelity objective $E(X, \theta)$ from this discretization (Equation 7, with Equation 4 for Euler), since this is the core loss function that drives all optimization.
Third, I introduce the proximal term and the full per-iteration objective $F_n(X, \theta)$ (Equation 8, with Equation 5 for Euler), explaining the crucial difference between anchoring to previous iterates versus anchoring to the original noisy data.
Fourth, I walk through the block coordinate descent procedure step by step (Aligned Equations 14a–14d), including the $x^+$ / $x^-$ splitting used in the convergence theory, and explain the stopping criterion.
Fifth, I derive the gradient-descent interpretation that reveals $\lambda$ 's role as inverse step size (Equation 12), which is essential for understanding why the algorithm's single hyperparameter is easy to set.
Sixth, I explain the final predicted-state generation step (Step 8 of Algorithm 1) and why it produces a (state, parameter) pair with exactly zero fidelity.

3.4 Detailed, Sentence-Based Technical Breakdown

This is primarily an algorithm design and analysis paper whose core idea is that combining block coordinate descent with a proximal term — where the proximal anchor is updated at each outer iteration to point at the previous state estimate rather than at the original noisy observations — yields an optimization procedure for simultaneous filtering and parameter estimation that is provably convergent under block-convexity assumptions and empirically robust to noise, initialization, and hyperparameter choice.

The ODE Discretization: From Continuous Dynamics to Algebraic Constraints

Before any optimization can occur, the continuous-time ODE system must be converted into a finite set of algebraic equations relating the unknown states at discrete time points. Given a dynamical system

$\dot{x}(t) = \frac{dx(t)}{dt} = f(x(t), \theta)$

with state $x(t) \in \mathbb{R}^d$ , parameter $\theta \in \mathbb{R}^p$ , and known vector field $f$ , the paper observes the system at $T$ time points $\{t_i\}_{i=1}^T$ with potentially non-uniform spacing $\Delta_i = t_{i+1} - t_i$ . The goal is to relate the unknown states at these discrete times — denoted $x(t_i)$ — to each other through a numerical integration scheme that approximates the true continuous evolution.

The paper supports two families of discretization, both of which are explicit linear multistep methods.

Explicit Euler (first-order, $m=1$ ). This is the simplest discretization, introduced in Equation (3):

$x(t_{i+1}) - x(t_i) = f(x(t_i), \theta) \Delta_i, \quad i = 1, \ldots, T-1$

where $x(t_i)$ is the state vector at time $t_i$ , $f(x(t_i), \theta)$ is the vector field evaluated at that state and parameter, and $\Delta_i$ is the time step between observations $i$ and $i+1$ .

What it computes: for each consecutive pair of time points, the change in state (left-hand side) should equal the instantaneous rate of change multiplied by the time duration (right-hand side). This is a first-order Taylor approximation: it assumes the derivative is constant over the interval $\Delta_i$ and equal to its value at the left endpoint.

Why this form matters: the Euler method has global error $O(h)$ where $h$ is the time step size. This means the numerical approximation becomes exact only as $h \to 0$ . For finite $h$ , there is a discrepancy between the true continuous solution and the Euler-discretized trajectory, even with correct parameters. This discretization error is separate from observation noise and parameter estimation error.

General explicit linear $m$ -step method (Equation 6). For higher accuracy, the paper allows any explicit linear multistep method of order $m$ , expressed in the general form:

$x(t_{i+1}) = \sum_{j=0}^{m-1} a_j \, x(t_{i-j}) + \Delta_i \sum_{j=0}^{m-1} b_j \, f(x(t_{i-j}), \theta)$

where $\{a_j\}_{j=0}^{m-1}$ and $\{b_j\}_{j=0}^{m-1}$ are method-specific coefficients (determined by the choice of multistep method, such as Adams-Bashforth), $m$ is the number of previous states used to predict the next state, and $f(x(t_{i-j}), \theta)$ is the vector field evaluated at each of those previous states.

What it computes: the state at time $t_{i+1}$ is a weighted combination of the previous $m$ states (the first sum) plus a weighted combination of the vector field evaluated at those same $m$ states scaled by the time step (the second sum). For example, in the three-step Adams-Bashforth method ( $m=3$ ), the coefficients approximate $f$ by fitting a quadratic polynomial through the three most recent derivative evaluations and integrating it over the interval, yielding $O(h^3)$ global error.

Why this form matters: higher-order methods dramatically reduce discretization error for the same step size, which is practically important because the fidelity objective $E$ measures exactly the discrepancy in this algebraic constraint. If the discretization is inaccurate, then even the true states and true parameters will produce a non-zero fidelity value, creating a floor on achievable optimization quality that is unrelated to the estimation problem.

Order adaptation at the start of the trajectory. A subtle but important implementation detail appears in the paper's description of the general $m$ -step method. To predict the state at time $t_i$ , the method needs the previous $m$ states. For $i \leq m$ , not enough previous states exist. The paper handles this by reducing the effective order: to predict $x(t_i)$ , the algorithm uses a multistep method of order $k = \min(i-1, m)$ . This means the first state $x(t_2)$ is always predicted with a 1-step method (Euler), $x(t_3)$ with a $\min(2, m)$ -step method, and so on, until the full $m$ -step method can be used from $x(t_{m+1})$ onward. This adaptive order ensures that every state has a well-defined prediction equation without requiring states before the initial time.

The Fidelity Objective: Measuring ODE Consistency

Given the discretized algebraic constraints, the paper defines a scalar objective function $E(X, \theta)$ that measures how well a candidate set of states $X = \{x(t_i)\}_{i=1}^T$ and parameters $\theta$ satisfy the discretized ODE. This is the "physics loss" that drives the optimization.

General $m$ -step fidelity (Equation 7):

$E_{m\text{-step}}(X, \theta) = \sum_{i=1}^{T-1} \left\| x(t_{i+1}) - \sum_{j=0}^{k-1} a_j x(t_{i-j}) - \Delta_i \sum_{j=0}^{k-1} b_j f(x(t_{i-j}), \theta) \right\|^2$

where $k = \min(i-1, m)$ is the effective order at step $i$ , $\{a_j\}$ and $\{b_j\}$ are the multistep coefficients for that order, $\|\cdot\|$ denotes the Euclidean norm (the paper uses $\|\cdot\|^2$ notation for squared Euclidean norm), and $X$ is the collection of all state vectors across all time points.

What it computes: for each time step $i$ , compute the predicted next state using the $k$ -step method (the terms after the minus sign), subtract this prediction from the candidate next state $x(t_{i+1})$ , square the Euclidean norm of the resulting $d$ -dimensional discrepancy vector, and sum across all $T-1$ time steps. A value of zero means the candidate $(X, \theta)$ pair satisfies the discretized ODE exactly at every time step. A positive value means there is some mismatch — either because the parameters are wrong, the states are wrong, or both.

Why this form: the sum-of-squares structure makes $E$ a smooth (infinitely differentiable if $f$ is smooth) function of both $X$ and $\theta$ , enabling gradient-based optimization. The squared Euclidean norm penalizes large per-step mismatches more heavily than small ones (quadratic growth), which encourages the optimizer to spread error evenly across time steps rather than concentrating it in a few large violations. This is the standard least-squares formulation that underlies much of numerical optimization.

Euler special case (Equation 4). For $m=1$ , the general form reduces to the Euler fidelity:

$E(X, \theta) = \sum_{i=1}^{T-1} \left\| x(t_{i+1}) - x(t_i) - f(x(t_i), \theta) \Delta_i \right\|^2$

This is simply Equation (7) with $a_0 = 1$ , $b_0 = 1$ , and $k=1$ for all $i$ . The Euler fidelity is what the paper uses throughout its theoretical convergence analysis (Section 4), while the experiments primarily use the three-step Adams-Bashforth version for improved accuracy.

Theorem 1 and the non-uniqueness problem. The paper proves (Section 3, Theorem 1) that $E(X, \theta)$ has infinitely many global minima achieving $E = 0$ . The proof is constructive: pick arbitrary values for $\theta$ and the initial condition $x(t_1)$ , then define all subsequent $x(t_i)$ by forward integration using the discretization scheme (Equation 6). By construction, every term in the sum vanishes, so $E = 0$ . Since $\theta$ and $x(t_1)$ are arbitrary, there are uncountably many such minima. This means that minimizing $E$ alone — without additional constraints — is an ill-posed problem: the optimizer could return any (state, parameter) pair that satisfies the ODE, regardless of whether it has any relationship to the observed data. This is the fundamental structural challenge that motivates the proximal term.

The Proximal Term and the Full Per-Iteration Objective

To resolve the non-uniqueness identified in Theorem 1, the paper augments the fidelity objective with a proximal penalty term that anchors the state estimate to a reference point. The crucial design choice is which reference point to use.

Per-iteration objective (Equation 8, with Equation 5 as the Euler special case):

$F_n(X, \theta) = E_{m\text{-step}}(X, \theta) + \lambda \left\| X - X^{*(n-1)} \right\|^2$

where $n \geq 1$ is the outer iteration index, $X^{*(n-1)}$ is the state estimate from the previous outer iteration (with $X^{*(0)} = Y$ , the noisy observations), $\lambda > 0$ is a hyperparameter controlling the proximal penalty strength, and $\|X - X^{*(n-1)}\|^2 = \sum_{i=1}^T \|x(t_i) - x^{*(n-1)}(t_i)\|^2$ is the sum of squared Euclidean distances between the candidate state trajectory and the previous estimate across all time points.

What it computes: in iteration $n$ , the algorithm seeks a (state, parameter) pair that simultaneously achieves low ODE fidelity $E$ (satisfying the physics) and stays close to the previous state estimate $X^{*(n-1)}$ (the proximal penalty). The hyperparameter $\lambda$ balances these two objectives: large $\lambda$ heavily penalizes deviation from the previous estimate, forcing small updates; small $\lambda$ allows large moves toward ODE-consistent trajectories.

Why this form — the critical distinction from iPDA: Section 3 provides an explicit comparison with iPDA's objective $J(x(t), \theta) = E_{\text{cont}}(x(t), \theta) + \lambda\|X - Y\|^2$ . In iPDA, the regularization anchors filtered states to the original noisy observations $Y$ , permanently. This means the filtered states can never move far from the noisy data, regardless of how much the ODE fidelity term suggests they should. If the noise magnitude is large, iPDA is trapped near bad data. In BCD-prox, the proximal anchor is $X^{*(n-1)}$ — a moving target that updates at each iteration. In iteration 1, the anchor is $Y$ (same as iPDA). But in iteration 2, the anchor is $X^{*(1)}$ , which may have already moved away from $Y$ toward an ODE-consistent trajectory. In iteration 3, the anchor is $X^{*(2)}$ , which may be even further from $Y$ . This progressive update mechanism allows BCD-prox to escape the initial noise:

"If the data $Y$ is heavily contaminated with noise, it may be wise to move far away from $Y$ as we iterate."

The proximal term effectively implements trust-region optimization: at each outer iteration, the algorithm is allowed to search for a better (state, parameter) pair only within a neighborhood of the current estimate, where the neighborhood radius is controlled by $\lambda$ . As the estimate improves, the trust region moves with it.

The initialization choice. Setting $X^{*(0)} = Y$ is the natural choice — the noisy observations are our only data, so they serve as the initial guess for the states. The initial parameter guess $\theta^{*(0)}$ must be supplied externally (the paper tests robustness to this choice extensively). Together, $(X^{*(0)}, \theta^{*(0)})$ define the starting point for optimization.

Block Coordinate Descent: Alternating Parameter and State Optimization

The core optimization procedure decomposes the minimization of $F_n(X, \theta)$ into two sequential subproblems — first over $\theta$ , then over $X$ — rather than attempting joint optimization over both simultaneously. This decomposition is what makes the algorithm "block coordinate descent."

The two minimization steps (from Algorithm 1 and Section 4):

Step 5 — Parameter block:

$\theta^{*(n)} = \arg\min_\theta \, F_n(X^{*(n-1)}, \theta)$

Since the proximal term $\lambda\|X - X^{*(n-1)}\|^2$ does not depend on $\theta$ , this simplifies to:

$\theta^{*(n)} = \arg\min_\theta \, E(X^{*(n-1)}, \theta)$

What it computes: holding the state estimate fixed at the previous iteration's value $X^{*(n-1)}$ , find the parameters $\theta$ that minimize the ODE fidelity. This is a standard parameter estimation problem with fixed regressors (the states). The dimensionality is $p$ (number of parameters, typically small — 1 to 4 in the paper's experiments).

Step 6 — State block:

$X^{*(n)} = \arg\min_X \, F_n(X, \theta^{*(n)})$

Here the full objective matters because both the fidelity term (which depends on $X$ ) and the proximal term (which penalizes deviation from $X^{*(n-1)}$ ) are active.

What it computes: holding the parameter estimate fixed at the just-computed value $\theta^{*(n)}$ , find the state trajectory $X$ that balances satisfying the ODE (with these parameters) against staying near the previous state estimate. The dimensionality is $d \times T$ (states times time points), which can be large — up to $40 \times 400 = 16,000$ for Lorenz-96 in the paper's experiments.

Proximal operator interpretation. Section 3 reinterprets Step 6 through the lens of proximal operators (Parikh & Boyd, 2014). For a function $g$ , the proximal operator with parameter $\alpha > 0$ is defined as:

$\text{prox}_{\alpha g}(v) = \arg\min_u \left( g(u) + \frac{1}{2\alpha} \|u - v\|^2 \right)$

Setting $g = E(\cdot, \theta^{*(n)})$ , $u = X$ , $v = X^{*(n-1)}$ , and $\alpha = (2\lambda)^{-1}$ , we obtain:

$X^{*(n)} = \text{prox}_{(2\lambda)^{-1}E}(X^{*(n-1)})$

What this means operationally: the state update at iteration $n$ is exactly the proximal operator of the fidelity function (with parameters fixed at $\theta^{*(n)}$ ) evaluated at the previous state estimate $X^{*(n-1)}$ , with step size parameter $(2\lambda)^{-1}$ .

Why this interpretation matters: the proximal operator generalizes gradient descent to non-smooth and non-convex functions. For a smooth function $E$ , the proximal step approximates an implicit gradient step. The paper derives the explicit relationship (Equation 12):

$X^{*(n)} = X^{*(n-1)} - (2\lambda)^{-1} \nabla_X E(X^{*(n-1)}, \theta^{*(n)}) + o((2\lambda)^{-1})$

This reveals that $\lambda$ functions as an inverse step size: increasing $\lambda$ makes $(2\lambda)^{-1}$ smaller, which reduces the effective step size, causing the algorithm to take smaller gradient steps and change more slowly. Decreasing $\lambda$ makes the step size larger, allowing more aggressive updates. The paper explicitly notes:

"It is now clear that $\lambda$ plays the role of an inverse step size — our experiments later will confirm that there is little harm in choosing $\lambda$ too large."

Taking $\lambda$ too large simply slows convergence (more iterations needed) without damaging the final solution quality. Taking $\lambda$ too small (including $\lambda = 0$ ) causes the algorithm to jump too far in a single step, potentially settling at an arbitrary ODE-consistent trajectory far from the true states. The paper reports that $\lambda = 0$ causes the algorithm to stop after one iteration with poor results.

The convergence guarantee (Equation 13). The block structure implies a monotonicity property:

$E(X^{*(n)}, \theta^{*(n)}) \leq E(X^{*(n-1)}, \theta^{*(n-1)})$

This holds because:

Step 5 minimizes $E$ over $\theta$ with $X$ fixed, so $E(X^{*(n-1)}, \theta^{*(n)}) \leq E(X^{*(n-1)}, \theta^{*(n-1)})$ .
Step 6 minimizes $F_n$ over $X$ with $\theta$ fixed. Since $F_n = E + \lambda\|X - X^{*(n-1)}\|^2$ , and $E$ is a component of $F_n$ , the fidelity at the new state cannot exceed the fidelity at the previous state (if it did, the optimizer could have chosen the previous state and achieved lower $F_n$ ).
Together, these give the chain of inequalities establishing that $E$ is non-increasing across iterations.

Since $E \geq 0$ always, the sequence $\{E(X^{*(n)}, \theta^{*(n)})\}_{n \geq 1}$ is bounded below and non-increasing, hence converges to some limit $E^* \geq 0$ .

Implementation Details: LBFGS and Automatic Differentiation

The paper specifies that both $\arg\min$ steps (Steps 5 and 6 of Algorithm 1) are implemented in practice using the LBFGS algorithm (a quasi-Newton method for smooth unconstrained optimization), accessed via scipy.optimize.minimize in Python. Gradients of the objective functions are supplied to LBFGS through automatic differentiation — the paper does not require manual derivation or implementation of gradients for each ODE system.

Why LBFGS: the fidelity objective $E$ is a sum of squares of smooth functions (if $f$ is smooth), making it amenable to quasi-Newton methods that approximate the Hessian using gradient history. LBFGS is memory-efficient (it stores only a limited history of gradients and updates) and converges faster than first-order methods (gradient descent) on well-conditioned problems. For the paper's problem sizes (up to $40 \times 400 = 16,000$ state variables plus a handful of parameters), LBFGS is computationally tractable.

Why automatic differentiation: the vector field $f$ varies across ODE models (Lotka–Volterra, FitzHugh–Nagumo, Rössler, Lorenz-96). Manually deriving and coding gradients for each model would be tedious and error-prone. Automatic differentiation computes exact gradients (not numerical finite differences) by tracing through the computation graph of the objective function, enabling the algorithm to be applied to any ODE system by simply specifying $f$ .

Stopping Criterion and the Three Termination Cases

The outer iteration (the loop over $n$ in Algorithm 1) stops "when the error $E$ changes less than $10^{-8}$ from one iteration to the next." In practice, this means the algorithm terminates when:

$|E(X^{*(n)}, \theta^{*(n)}) - E(X^{*(n-1)}, \theta^{*(n-1)})| < 10^{-8}$

Section 4 provides a structural analysis of why this stopping criterion works by identifying three possible outcomes of the state optimization in Step 6:

Case 1: $X^{*(n)} = \hat{X}^{(n)}$ , where $\hat{X}^{(n)}$ is the "predicted states" — the states obtained by forward-integrating the ODE from the current initial condition $x^{*(n)}(t_1)$ using the current parameter estimate $\theta^{*(n)}$ . By construction (as in Theorem 1's proof), $E(\hat{X}^{(n)}, \theta^{*(n)}) = 0$ . In this case, the fidelity has been driven to zero, and further optimization over $\theta$ (which would hold these zero-fidelity states fixed) cannot improve. The algorithm stops.

Case 2: $X^{*(n)} = X^{*(n-1)}$ — the state estimate has not changed. Since the parameter update in Step 5 depends on the previous state, and the state update in Step 6 depends on the current parameters, if neither changes, subsequent iterations will produce identical results. The algorithm stops.

Case 3: $X^{*(n)}$ is neither $\hat{X}^{(n)}$ nor $X^{*(n-1)}$ — the state estimate has moved toward the predicted states but has not reached them, and has also changed from the previous estimate. The algorithm continues to the next iteration because further optimization can still decrease the fidelity.

The paper notes that Case 1 and Case 2 are the convergent states: in Case 1, the algorithm has found a (state, parameter) pair that satisfies the ODE exactly; in Case 2, it has converged to a fixed point that may not achieve zero fidelity but cannot be improved further by alternating optimization. Case 3 is the transient state that drives continued optimization.

Block Splitting for Convergence Theory: The $x^+$ / $x^-$ Decomposition

Section 4 introduces an alternative block decomposition used specifically for the convergence proofs, distinct from the practical $\theta$ -then- $X$ decomposition in Algorithm 1. This theoretical decomposition splits the state trajectory into two halves:

For $T$ even: $x^+ = \{x(t_1), x(t_2), \ldots, x(t_{T/2})\}$ (first half of the time series) and $x^- = \{x(t_{T/2+1}), \ldots, x(t_T)\}$ (second half). For $T$ odd, replace $T/2$ with $(T-1)/2$ .

The three-block procedure for the convergence analysis (Equations 14a–14d) is:

$\theta_n = \arg\min_\theta F_n(X_{n-1}, \theta) = \arg\min_\theta E(X_{n-1}, \theta)$

$(x^-)_n = \arg\min_{x^-} F_n((x^+)_{n-1}, x^-, \theta_n)$

$(x^+)_n = \arg\min_{x^+} F_n(x^+, (x^-)_n, \theta_n)$

$X_n = ((x^+)_n, (x^-)_n)$

What this accomplishes: by splitting the state optimization into two sub-blocks ( $x^-$ then $x^+$ ) rather than optimizing all states simultaneously, the convergence analysis can leverage existing results from Xu & Yin (2013) on block coordinate descent for multiconvex optimization. The key requirement is that the objective $E(X, \theta)$ is strongly convex in each block when the other blocks are held fixed.

Why $x^+$ / $x^-$ specifically: the Euler discretization creates a chain-structured dependency: $x(t_{i+1})$ depends on $x(t_i)$ through the forward Euler step. The $x^+$ / $x^-$ split breaks this chain at the midpoint. When optimizing $x^-$ with $x^+$ fixed, the last state in $x^+$ (at time $T/2$ ) provides a fixed "initial condition" for the first state in $x^-$ (at time $T/2+1$ ), making the $x^-$ optimization well-posed. Similarly, when optimizing $x^+$ with $x^-$ fixed, the first state in $x^-$ provides a fixed "terminal condition" for the last state in $x^+$ . This two-directional anchoring creates strong convexity.

The full-rank condition on $f_1$ . The paper assumes $f$ is at most linear in $\theta$ , so $f(x, \theta) = f_0(x) + f_1(x)\theta$ , with $f_1: \mathbb{R}^d \to \mathbb{R}^{d \times p}$ having full column rank for all $x$ . This assumption ensures the Hessian of $E$ with respect to $\theta$ is positive definite:

$\nabla_\theta \nabla_\theta E = 2 \sum_{i=1}^{T-1} (f_1(x_i))^T f_1(x_i) \Delta_i^2$

Since $f_1$ has full column rank, each term $(f_1)^T f_1$ is positive definite, and the sum of positive definite matrices is positive definite, establishing strong convexity in $\theta$ .

Strong convexity in $x^+$ and $x^-$ at small $\Delta_i$ . The paper proves strong convexity of $E$ restricted to $x^+$ (with $x^-$ and $\theta$ fixed) and to $x^-$ (with $x^+$ and $\theta$ fixed) when all $\Delta_i$ are sufficiently small. The proof proceeds by:

Evaluating the Hessian at $\Delta_i = 0$ (all time steps zero), where the Euler fidelity reduces to $E(X, \theta; \Delta_i = 0) = \sum_{i=1}^{T-1} \|x_{i+1} - x_i\|^2$ , a quadratic form.
Computing the block Hessian $A = \nabla_{x^+}\nabla_{x^+}E(X, \theta; \Delta_i = 0)$ explicitly as a block tridiagonal matrix with $2I$ blocks on the diagonal except at endpoints, and $-2I$ blocks on the off-diagonals.
Proving $\det A = 2^{dT/2}$ by induction, establishing that all eigenvalues are positive (bounded away from zero).
Invoking continuity of eigenvalues as functions of $\Delta_i$ to assert the existence of $\delta > 0$ such that strong convexity holds for all $\Delta_i \in (0, \delta)$ .

This construction yields Theorem 2: under the linear-in- $\theta$ , full-rank, and small- $\Delta_i$ conditions, the three-block procedure (14) converges to a Nash equilibrium of $E$ . Adding the Kurdyka-Łojasiewicz (KL) property (satisfied when $f$ is real analytic) yields Theorem 3: initialized sufficiently close to a global minimizer, the algorithm converges to a global minimizer of $E$ .

Final Predicted States: The Guaranteed Zero-Fidelity Output

After the BCD-prox iteration converges (producing final estimates $\theta^*$ and $X^*$ ), Algorithm 1 includes an additional step (Step 8) that generates "predicted states" $\hat{X}$ :

"Compute predicted states $\hat{X}$ by repeatedly applying Eq. (6), where $\theta = \theta^{*(n)}$ and initial condition $x(t_1) = x^{*(n)}(t_1)$ ."

What this does: take the final estimated parameters $\theta^*$ and the final estimated initial state $x^*(t_1)$ , then forward-integrate the ODE using the same $m$ -step discretization method to produce a complete state trajectory $\hat{X} = \{\hat{x}(t_1), \hat{x}(t_2), \ldots, \hat{x}(t_T)\}$ where $\hat{x}(t_1) = x^*(t_1)$ .

The guarantee: by construction (as in the proof of Theorem 1), this trajectory satisfies $E(\hat{X}, \theta^*) = 0$ exactly — the fidelity is identically zero because each term in the sum (7) vanishes by the forward-integration definition. The predicted states $\hat{X}$ represent the "ODE-consistent" state trajectory that best matches the final parameter and initial condition estimates.

Why there are two different state outputs: the algorithm produces both $X^*$ (the estimated states from the optimization) and $\hat{X}$ (the predicted states from forward integration). In general, $X^* \neq \hat{X}$ because $X^*$ balances ODE fidelity against the proximal penalty inherited from the optimization history — it may not satisfy the discretized ODE exactly. $\hat{X}$ , by contrast, discards the proximal penalty influence and enforces exact ODE consistency, using only the estimated parameters and initial condition. The paper reports "prediction error" as $\|X - \hat{X}\|$ (Frobenius norm of the difference between true clean states and predicted states) and "estimation error" as $\|X - X^*\|$ (difference between true clean states and estimated states). Both metrics appear in the experiments, with the distinction becoming important in the EKF comparison (where EKF cannot produce good predictions because it never updates the initial state).

Why this two-output design: the proximal penalty is a regularization artifact that improves optimization but biases the state estimates toward the previous iterate. The forward-integration step removes this bias by projecting the (parameter, initial condition) estimate onto the manifold of exactly ODE-consistent trajectories. This separation of concerns — use the proximal penalty for robust optimization, then remove its influence for final prediction — is a key algorithmic insight that contributes to BCD-prox's strong empirical performance.

4. Key Insights and Innovations

Innovation 1: The Proximal Anchor as a Moving Target Rather Than a Fixed Regularization

The dominant conceptual move in BCD-prox — the one that distinguishes it from all prior work on simultaneous filtering and parameter estimation — is the decision to tether the state estimate at iteration $n$ to the previous iteration's state estimate $X^{*(n-1)}$ rather than to the original noisy observations $Y$ . This is not an incremental parameter-tuning choice; it fundamentally changes what the regularization term means and what the algorithm is allowed to do.

In the spline-based iPDA framework (Ramsay et al., 2007; Poyton et al., 2006), the objective takes the form $J(x(t), \theta) = E_{\text{cont}}(x(t), \theta) + \lambda\|X - Y\|^2$ . The regularization penalizes deviation from the noisy data — permanently. The filtered states can never escape the neighborhood of the observations, regardless of how strongly the ODE fidelity term argues that the true trajectory lies elsewhere. This is the right design when noise is small: the data is trustworthy, so filtered states should stay close to it. It is catastrophically wrong when noise is large: the data is misleading, the ODE fidelity points toward a very different trajectory, but the regularization refuses to let the states follow. The algorithm is stuck optimizing a contradiction, and Section 3 already showed why this is structurally unavoidable — iPDA cannot tune $\lambda$ properly without access to the clean states, which are unknown by definition.

BCD-prox replaces the static anchor $Y$ with a dynamic anchor $X^{*(n-1)}$ . At iteration 1, the anchor is $Y$ (same starting point). But at iteration 2, the anchor is $X^{*(1)}$ , which may already be closer to the true trajectory than $Y$ was. At iteration 3, the anchor has moved again. The proximal term no longer says "stay near the noisy data"; it says "don't change too much from where you were last iteration." This is not regularization in the classical sense of preventing overfitting to noise — it is trust-region optimization, where the trust region follows the iterate. If the ODE fidelity pulls the states far from $Y$ , the algorithm can follow, one proximal step at a time. The paper makes this point explicitly in Section 3:

"If the data $Y$ is heavily contaminated with noise, it may be wise to move far away from $Y$ as we iterate."

This is a fundamental conceptual shift, not a parameter-tuning trick. The prior literature treated filtering as a problem of smoothing noisy observations toward a latent truth, with the observations as the fixed reference. BCD-prox treats filtering as an optimization problem where the reference itself evolves. The significance beyond raw performance is that this reframing eliminates the central hyperparameter dilemma that plagued spline-based methods — the optimal $\lambda$ in iPDA depends on the unknown noise $Z$ and the unknown vector field $f$ , making it impossible to set correctly in practice. BCD-prox's $\lambda$ controls step size, not data fidelity, and the paper demonstrates (Figure 2) that any $\lambda > 0$ works, with values from 1 to 1,000 producing essentially identical final estimates. The hyperparameter becomes harmless to set large — it only slows convergence, never traps the algorithm near bad data.

Innovation 2: The Infinite-Minima Problem as an Explicit Design Constraint, Not a Pathological Edge Case

Theorem 1 (Section 3) proves that the fidelity objective $E(X, \theta)$ has infinitely many global minima achieving $E = 0$ — any arbitrary choice of $\theta$ and initial condition $x(t_1)$ can be forward-integrated to produce a state trajectory that satisfies the discretized ODE exactly. This result is mathematically trivial (the proof is three sentences), but its recognition as a central design constraint is what constitutes the intellectual contribution.

Prior work implicitly acknowledged the non-uniqueness problem through the use of regularization, but the degeneracy was treated as a nuisance to be papered over rather than a structural fact that should drive algorithm design. Spline-based methods added a data-fidelity penalty; Bayesian methods added priors; nonlinear least squares relied on good initialization to converge to the "right" minimum. None of these approaches engaged directly with the question: given that infinitely many (state, parameter) pairs satisfy the ODE perfectly, what principle should distinguish the correct one?

BCD-prox's answer is implicit in its architecture but no less principled: the correct solution is the one reachable from the noisy observations through a sequence of proximal-gradient steps that jointly optimize fidelity and step-size control. The algorithm does not attempt to pick the "best" minimum among the infinite set — it lets the optimization trajectory, starting from $Y$ and constrained by the moving proximal anchor, determine which minimum it converges to. This is a fundamentally different resolution strategy than imposing a prior (Bayesian) or enforcing proximity to data (iPDA). It is closer in spirit to the idea that the optimization path itself encodes useful inductive bias — an idea that has gained traction in deep learning (the implicit regularization of gradient descent) but had not been applied to ODE parameter estimation.

The practical consequence is visible in the initialization-robustness experiments (Figures 3 and 7). When parameters are initialized far from their true values (initialization noise variance $\sigma^2_\theta = 20$ ), BCD-prox still converges to good estimates, while iPDA, Bayesian, and least-squares methods produce errors exceeding 100 in roughly half of all trials. The infinite-minima theorem explains why: competing methods search for a minimum of a modified objective (fidelity plus regularization), but when the regularization is poorly calibrated (because the noise is unknown), the modified objective's minimum may be far from the true solution. BCD-prox, by contrast, does not attempt to directly select a minimum — it follows a path. The path's endpoint depends on the data and the ODE structure, not on an arbitrary hyperparameter balance. This represents a genuine conceptual advance: treating the infinite-minima degeneracy not as a problem to suppress but as a landscape feature to navigate.

Innovation 3: The Convergence Theory via Blockwise Strong Convexity with an Explicit $\Delta_i$ Threshold

The convergence analysis in Section 4 establishes conditions under which BCD-prox provably converges to a Nash equilibrium (Theorem 2) or a global minimizer (Theorem 3) of the Euler fidelity objective. The theoretical contribution is not the application of Xu & Yin (2013)'s block coordinate descent framework — that machinery is general. The contribution is the construction of blockwise strong convexity through the $x^+$ / $x^-$ split and the explicit identification of a $\Delta_i$ threshold that guarantees convexity in each block.

The paper proves strong convexity of $E$ restricted to $x^+$ (first half of the time series, with $x^-$ and $\theta$ fixed) and to $x^-$ (second half, with $x^+$ and $\theta$ fixed) when all $\Delta_i$ are sufficiently small. The proof strategy is instructive: first evaluate the Hessian at $\Delta_i = 0$ (zero time steps), where the Euler fidelity reduces to a simple quadratic form $\sum \|x_{i+1} - x_i\|^2$ whose blockwise positive definiteness can be established analytically (determinant exactly $2^{dT/2}$ ). Then invoke continuity of eigenvalues as functions of $\Delta_i$ to assert the existence of $\delta > 0$ such that strong convexity persists for $\Delta_i \in (0, \delta)$ .

This is a constructive existence result rather than an asymptotic one — it says there is a concrete, non-zero range of time step sizes for which the optimization landscape is provably well-behaved. This matters because ODE parameter estimation typically operates at finite time steps (observations are discrete), and existing convergence theory for non-convex optimization often relies on asymptotic arguments or Lipschitz conditions that are hard to verify for specific ODE systems. The paper's approach, by contrast, ties the convergence guarantee to a computable property of the discretization (the $\Delta_i$ values) and a verifiable property of the vector field (linearity in $\theta$ with full-rank $f_1$ ).

The practical significance is that all four ODE systems tested in the experiments — Lotka–Volterra, FitzHugh–Nagumo, Rössler, and Lorenz-96 — satisfy the KL property (their vector fields are real-analytic) and the linear-in-parameters condition (after a change of variables for FitzHugh–Nagumo, as noted in Section 4). The theory thus covers systems of genuine scientific interest, not just toy examples. This is a meaningful advance over prior convergence results for alternating minimization applied to ODE problems, which typically required stronger assumptions (global convexity, restrictive parameterizations) or provided only asymptotic guarantees.

Innovation 4: The Predicted-State Post-Processing Step as an Implicit Projection onto the ODE-Consistent Manifold

Algorithm 1's Step 8 — generating predicted states $\hat{X}$ by forward-integrating the ODE from the estimated initial condition and parameters — is easy to overlook as a minor bookkeeping detail. It is not. It represents a separation of the optimization problem from the evaluation problem that has significant conceptual and practical consequences.

During optimization, the state estimate $X^*$ balances two competing objectives: ODE fidelity and proximal anchoring. The proximal term biases $X^*$ toward the previous iterate, meaning $X^*$ may not satisfy the ODE exactly — the fidelity $E(X^*, \theta^*)$ may be small but non-zero. The predicted states $\hat{X}$ , by contrast, are constructed to achieve $E(\hat{X}, \theta^*) = 0$ identically, by definition of forward integration. Step 8 effectively projects the (parameter, initial condition) estimate onto the zero-fidelity manifold and reads off the corresponding state trajectory.

This matters because it decouples two roles that are conflated in most prior methods. The proximal term is needed during optimization to prevent the algorithm from jumping to an arbitrary ODE-consistent trajectory (the infinite-minima problem). But once optimization converges, the proximal term has served its purpose — retaining its influence on the final state estimate would introduce bias (states would remain artificially close to earlier iterates, which were themselves influenced by the noisy data). The forward-integration step removes this bias: the final answer depends only on the estimated parameters and initial condition, which are the quantities of scientific interest. The states are reconstructed deterministically from these estimates via the known dynamics.

The practical value of this separation is most visible in the extended Kalman filter comparison. EKF, as an online method, never revisits its initial state estimate. If the initial state is noisy, every subsequent prediction is corrupted, even if parameters are estimated perfectly. The paper quantifies this: EKF's prediction error "goes to infinity" (Section 5) because the initial state error propagates through the entire trajectory. BCD-prox's post-processing step avoids this — it takes the estimated initial state (which has been refined through all outer iterations) and generates a clean trajectory. The estimation error (difference between true states and $X^*$ ) and prediction error (difference between true states and $\hat{X}$ ) are reported separately, allowing users to assess both the quality of the optimization output and the quality of the ODE-consistent reconstruction. This diagnostic clarity — being able to distinguish "did the optimization converge to the right parameters?" from "does the ODE model fit the data?" — is a subtle but practically valuable contribution that prior approaches do not offer.

5. Experimental Analysis

Evaluation Methodology

Dataset. The paper uses four benchmark ODE systems — Lotka–Volterra, FitzHugh–Nagumo, the Rössler attractor, and Lorenz-96 — with the specific equations and true parameter values detailed in the supplementary material (Section 6). Clean states are generated via a fifth-order Runge-Kutta method, then corrupted with Gaussian noise (zero mean, variance $\sigma^2$ ) to produce the noisy observations $Y$ used for evaluation. The data generation process includes varying numbers of time points ( $T = 20$ or $T = 10{,}000$ ), varying time ranges ( $t \in [0, 2]$ or $[0, 1{,}000]$ ), and systematic sweeps of noise variance ( $\sigma^2 \in \{0.1, 0.5, 1, 1.5\}$ depending on the experiment), without a fixed train/validation/test split — all experiments are in the regime where there is no separate "training data," only the noisy observations themselves.
Base model(s). BCD-prox is not a learned model but an optimization algorithm. It relies on a user-specified ODE vector field $f(x, \theta)$ (assumed known in form but with unknown parameters $\theta$ ) and a choice of numerical integration order $m$ for the multistep discretization method. The paper uses the three-step Adams-Bashforth method as the default ( $m=3$ ), with the explicit Euler method ( $m=1$ ) included as a comparison point in the error-vs-iteration analysis (Figure 1 and supplementary Figure 5).
Metrics. Two distinct error metrics are used:
- Prediction error: the Frobenius norm $\|X - \hat{X}\|_F$ between the true clean states $X$ and the predicted states $\hat{X}$ (computed by forward-integrating the ODE from the estimated initial condition and parameters, guaranteed to achieve $E(\hat{X}, \theta^*) = 0$ ).
- Parameter error: the absolute difference $|\theta_l - \theta^*_l|$ for each parameter component $l$ , reported individually per parameter rather than as a single aggregated metric.
The paper also uses estimation error (Frobenius norm between true states $X$ and estimated states $X^*$ ) specifically in the EKF comparison, because EKF cannot produce reasonable predictions (its initial state estimate is frozen at the noisy observation). All errors are reported as boxplots showing distributions across 10 independent trials (10 different realizations of noisy observations and/or parameter initializations), with average error values printed below each plot.
Baselines. Five competing methods from distinct methodological families are compared against BCD-prox:
- iPDA (iterated principal differential analysis; Ramsay et al., 2007; Poyton et al., 2006): a spline-based method implemented via publicly available MATLAB code. Represents the spline-based approach that smooths noisy data and then estimates parameters from the smoothed trajectories.
- Bayes (Bayesian approach; Dondelinger et al., 2013): implemented via publicly available R code. Represents the Bayesian gradient-matching approach using Gaussian processes.
- lsq (iterative least squares): implemented by the authors using the Python LMFIT package (Newville et al., 2014). Treats parameters and initial state as the unknown variables and applies iterative nonlinear least squares optimization directly to the ODE integration.
- Mean-field variational Bayes (Gorbach et al., 2017): a variational inference method applicable only to ODEs with a specific structural form. Compared on the Lotka–Volterra model only (the other models do not satisfy the required form).
- Extended Kalman Filter (EKF) (Sitz et al., 2002): implemented by the authors using an open-source Python library (Labbe, 2014). The joint state variable $\xi(t_i)$ concatenates the states and parameters; the process model uses Euler discretization; the observation model extracts only the state components via a selection matrix $H$ . Process and observation covariance matrices are tuned carefully (state covariance diagonal at 1,000, process covariance variance set to 1 via the Q_discrete_white_noise() function).
Generation budget / compute accounting. There is no "generation budget" in the LLM sense since no sampling is involved. The paper reports training time qualitatively (BCD-prox: ~20 seconds on a standard laptop for 400-state problems; spline-based methods: minutes; Bayesian methods: hours) and outer iteration counts visually in Figure 1 and supplementary Figure 5 (x-axis shows iteration number up to ~250 for FitzHugh–Nagumo and ~9,000 for Rössler, depending on convergence speed). The stopping criterion is quantitative but based on fidelity change ( $\Delta E < 10^{-8}$ ), not wall-clock time or FLOPs.
Cross-validation / statistical protocol. No cross-validation is used — this is not a supervised learning setting. Instead, because both observation noise and parameter initialization involve randomness, the paper repeats each experiment 10 times with different random seeds (10 different noisy observation sets and/or 10 different parameter initializations), reporting results as grouped boxplots where each bar represents one trial. The paper explicitly states: "Since there is randomness in both initialization and observation, we repeat the experiment 10 times. Note that the comparisons are fair, with the same observations and initializations used across all methods." This paired-trial design ensures that any advantage of one method over another on a particular trial cannot be attributed to that method receiving an easier noise realization or initialization.

Main Quantitative Results

5.1 Fidelity Optimization Decreases Prediction Error Across Iterations (Figure 1, Supplementary Figure 5)

Figure 1 (main paper) tracks the prediction error of BCD-prox at each outer iteration for the Rössler attractor ( $d=3$ , $\sigma^2=1$ ) and Lorenz-96 ( $d=40$ , $\sigma^2=1$ ), comparing two discretization schemes: explicit Euler (1-step) and three-step Adams-Bashforth. The supplementary Figure 5 adds FitzHugh–Nagumo and replicates results for $\sigma^2=0.5$ .

The headline result is that BCD-prox reduces prediction error by orders of magnitude across all models and both noise levels, with three-step Adams-Bashforth consistently outperforming Euler. For Lorenz-96 with $\sigma^2=1$ , the prediction error starts at approximately $1.5 \times 10^5$ and drops to near zero under Adams-Bashforth, while Euler plateaus around $10^4$ . For Rössler, the initial error of roughly $9 \times 10^3$ drops close to zero under both methods, but Adams-Bashforth achieves lower final error and converges faster. For FitzHugh–Nagumo (supplementary Figure 5), initial errors around $2 \times 10^3$ drop to approximately $10^2$ under both methods.

The error does not decrease monotonically — occasional increases are visible in the early iterations (most notably in the Lorenz-96 Euler trace). The paper attributes this to the fact that the optimization objective $E$ (which is being directly minimized) differs from the prediction error (which is the evaluation metric and cannot be directly optimized since clean states are unknown in practice). Despite this non-monotonicity, the overall trend is strongly downward, and the final converged prediction errors are very low.

5.2 Robustness to the Single Hyperparameter $\lambda$ (Figure 2)

Figure 2 sweeps $\lambda$ across the set $\{0, 3, 5, 7, 10, 20\}$ for both the Lotka–Volterra model ( $\sigma^2=0.5$ , true parameters $\theta = [2, 1, 4, 1]$ ) and the FitzHugh–Nagumo model ( $\sigma^2=0.5$ , true parameters $\theta = [0.5, 0.3, 3]$ ). For each $\lambda$ value, 10 independent trials are run with different noisy observation sets; the figure reports mean prediction error with standard deviation bars and estimated parameter values (error bars omitted for parameters to avoid clutter).

Three key findings emerge:

For any $\lambda > 0$ , BCD-prox correctly recovers the true parameters and drives prediction error close to zero. In the FitzHugh–Nagumo model, estimated parameters for $\lambda \in [3, 20]$ cluster tightly around the true values ( $\theta_0 \approx 0.5$ , $\theta_1 \approx 0.3$ , $\theta_2 \approx 3$ ). In the Lotka–Volterra model, the same holds ( $\theta_0 \approx 2$ , $\theta_1 \approx 1$ , $\theta_2 \approx 4$ , $\theta_3 \approx 1$ ).
In the range $\lambda = 1$ to $\lambda = 20$ , errors and estimated parameters remain essentially unchanged. The paper further reports (in prose, not in the figure itself) that "increasing $\lambda$ to $1{,}000$ does not change the estimated parameters." The only cost of excessively large $\lambda$ is longer training time, consistent with the interpretation of $\lambda$ as inverse step size — smaller steps require more iterations.
$\lambda = 0$ fails: the algorithm stops after a single iteration with predicted states far from the clean states. This is expected because without the proximal term, the algorithm finds an arbitrary zero-fidelity (state, parameter) pair with no relationship to the data, as guaranteed by Theorem 1.

The paper's choice to fix $\lambda = 1$ in all subsequent experiments is justified by these results — the value is within the broad plateau of good performance, ensuring fast convergence without sacrificing accuracy.

5.3 Comparison with iPDA, Bayesian, and Least-Squares Methods Under Varying Parameter Initialization (Figure 3, Supplementary Figure 7)

This is the paper's primary method-comparison experiment. The Rössler attractor results appear in Figure 3 (main paper); the FitzHugh–Nagumo results appear in supplementary Figure 7. Observation noise is fixed at $\sigma^2 = 0.5$ . Parameter initialization noise variance $\sigma^2_\theta$ is swept across $\{1, 5, 10, 20\}$ , where initialization means adding zero-mean Gaussian noise of the specified variance to the true parameters. For each combination of model and $\sigma^2_\theta$ , 10 independent trials are run; all methods receive identical observations and initializations per trial. BCD-prox uses $\lambda = 1$ throughout; competing methods use "the best hyperparameters that we could determine after careful experimentation."

Headline finding: BCD-prox is robust to parameter initialization; competing methods are not. The paper states:

"The total number of experiments per method is 80 (40 for the FitzHugh–Nagumo and 40 for Rössler). The prediction error of BCD-prox exceeds 100 in 4 experiments. The prediction error of iPDA (the second best method after ours) exceeds 100 in 39 experiments (nearly half the experiments). For lsq and Bayes, the errors are substantially worse."

At the largest initialization noise ( $\sigma^2_\theta = 20$ , the rightmost column of Figure 3), the contrast is stark:

BCD-prox: prediction error bars remain tightly clustered below 20 (Rössler), with median around 5–10. Parameter errors for all three Rössler parameters ( $\theta_0, \theta_1, \theta_2$ ) are consistently below 0.5, with most trials achieving errors below 0.2.
iPDA: prediction errors span a wide range, with several trials exceeding 80–100 (Rössler). Parameter errors show large variance — some trials get lucky (errors near 0.1) while others fail catastrophically (errors exceeding 2.0).
lsq: prediction errors cluster in the 60–100 range, with uniformly poor parameter estimates. Most parameter error bars exceed 1.0–2.0.
Bayes: similar to lsq, with prediction error bars mostly in the 60–100 range and parameter errors generally exceeding 1.0–2.0, showing the worst performance among the four methods.

At small initialization noise ( $\sigma^2_\theta = 1$ , leftmost column), all methods perform reasonably well, with prediction errors below 10–20 and parameter errors generally below 0.5. This confirms the paper's claim that the differences emerge specifically under poor initialization — the regime most relevant to real applications where true parameters are unknown.

The FitzHugh–Nagumo results (supplementary Figure 7) show qualitatively identical patterns: BCD-prox dominates across all $\sigma^2_\theta$ levels, iPDA shows large variance with some catastrophic failures, and lsq/Bayes produce uniformly poor estimates at high initialization noise.

5.4 Comparison with Mean-Field Variational Bayes (Figure 4, First Row; Supplementary Figure 8)

The mean-field variational inference method of Gorbach et al. (2017) is applicable only to the Lotka–Volterra model (the other models do not satisfy the required ODE form). Figure 4 (first row) compares BCD-prox against the mean-field method at $\sigma^2 = 1$ , with $T = 20$ observations. Supplementary Figure 8 extends the comparison to $\sigma^2 = 0.5$ and $\sigma^2 = 1.5$ .

At $\sigma^2 = 1$ (Figure 4, first row):

Prediction error: BCD-prox average = 5.53; mean-field average = 8.32. BCD-prox outperforms in 8 out of 10 trials.
Parameter $\theta_0$ : BCD-prox average error = 0.54; mean-field average = 7.41. This is the most dramatic gap — the mean-field method fails catastrophically on $\theta_0$ , with errors exceeding 3 in multiple trials, while BCD-prox maintains errors below 1.
Parameter $\theta_1$ : BCD-prox average = 0.27; mean-field average = 3.69. Similar pattern.
Parameter $\theta_2$ : BCD-prox average = 0.83; mean-field average = 1.70.
Parameter $\theta_3$ : BCD-prox average = 0.24; mean-field average = 0.43. This parameter is estimated comparably by both methods.

At increasing noise levels (supplementary Figure 8), both methods degrade, but BCD-prox remains consistently more robust. At $\sigma^2 = 1.5$ , BCD-prox's prediction error average rises to 13.04 (vs. 13.61 for mean-field — comparable), but BCD-prox continues to substantially outperform on $\theta_0$ (0.69 vs. 0.82) and $\theta_2$ (1.67 vs. 2.38), while maintaining rough parity on $\theta_1$ (0.36 vs. 0.39) and $\theta_3$ (0.46 vs. 0.65).

5.5 Comparison with Extended Kalman Filter (Figure 4, Second and Third Rows; Supplementary Figure 9)

The EKF comparison in Figure 4 uses the Lotka–Volterra model with a long time series ( $T = 10{,}000$ , time range $[0, 1{,}000]$ ) at two noise levels: $\sigma^2 = 0.1$ (low noise, second row) and $\sigma^2 = 1.5$ (high noise, third row). Supplementary Figure 9 adds results for a short time series ( $T = 20$ , time range $[0, 2]$ ) at both noise levels. The metric reported is estimation error (difference between true clean states and the estimated states $X^*$ ), not prediction error, because EKF's prediction error "goes to infinity" — the frozen noisy initial state corrupts all forward predictions regardless of parameter estimation quality.

Low noise, long time series ( $\sigma^2 = 0.1$ , $T = 10{,}000$ ): This is the only regime where EKF performs comparably to BCD-prox. Average estimation errors are 0.00 (BCD-prox) vs. 0.09 (EKF). Parameter errors are broadly similar: $\theta_0$ (0.09 vs. 0.06), $\theta_1$ (0.04 vs. 0.09), $\theta_2$ (0.16 vs. 0.07), $\theta_3$ (0.04 vs. 0.06). EKF slightly wins on some parameters, BCD-prox on others — the methods are essentially tied.

High noise, long time series ( $\sigma^2 = 1.5$ , $T = 10{,}000$ ): BCD-prox dramatically outperforms EKF. Average estimation errors are 0.06 (BCD-prox) vs. 1.30 (EKF) — a factor of ~22× improvement. All four parameter estimates favor BCD-prox: $\theta_0$ (0.11 vs. 0.93), $\theta_1$ (0.02 vs. 0.61), $\theta_2$ (0.19 vs. 0.63), $\theta_3$ (0.05 vs. 0.38). The paper explains this gap as a consequence of BCD-prox's batch nature:

"A key difference between the two methods is that EKF is an online method while ours is a batch method, iterating over the entire data set repeatedly. Consequently, BCD-prox updates parameters based on information in all the states, leading to more robust updates than is possible with EKF, which updates parameters based on a single observation."

Short time series ( $T = 20$ , supplementary Figure 9): At $\sigma^2 = 0.1$ , BCD-prox (average estimation error 0.03) outperforms EKF (0.15). At $\sigma^2 = 1.5$ , BCD-prox (0.47) outperforms EKF (3.56) by a factor of ~7.5×. The pattern is consistent: EKF requires many observations and low noise to compete; BCD-prox works well in the more realistic short, noisy regimes. Supplementary Figure 9 also tracks how both methods improve as $T$ increases, with BCD-prox's average estimation error approaching zero for sufficiently large $T$ .

5.6 Effect of Noise Type and Magnitude (Supplementary Figure 6)

Supplementary Figure 6 investigates whether BCD-prox's performance depends on the statistical assumptions about observation noise. Two noise types — Gaussian and Laplacian (heavier-tailed) — are compared at variances $\{0.1, 0.3, 0.5, 0.7, 1.0\}$ for FitzHugh–Nagumo and Rössler. Each (noise type, noise variance) combination is repeated 10 times; means and standard deviations are reported.

Finding: BCD-prox performs almost identically under Gaussian and Laplacian noise. The prediction error curves for the two noise types closely track each other across all variance levels for both models. Similarly, parameter errors for $\theta_0$ , $\theta_1$ , and $\theta_2$ show no systematic difference between Gaussian and Laplacian noise. This confirms the paper's claim that BCD-prox is distribution-free — it does not rely on Gaussian assumptions about observation noise, unlike Kalman filter methods.

Effect of increasing noise variance: as expected, prediction and parameter errors increase monotonically with noise variance for both models. However, the increase is gradual: for FitzHugh–Nagumo, prediction error rises from approximately 2 (at variance 0.1) to approximately 10–12 (at variance 1.0); for Rössler, from approximately 5 to approximately 15. BCD-prox does not exhibit any threshold behavior where performance suddenly collapses — it degrades gracefully with increasing noise.

Ablation Studies and Robustness Checks

Discretization order (Euler vs. Adams-Bashforth, Figure 1 and supplementary Figure 5): Three-step Adams-Bashforth consistently achieves both faster convergence (fewer iterations to reach low error) and lower final prediction error than Euler across all four ODE models and both noise levels ( $\sigma^2 = 0.5$ and $\sigma^2 = 1.0$ ). The advantage is most dramatic on Lorenz-96 ( $d=40$ ), where Adams-Bashforth reduces prediction error to near zero while Euler plateaus around $10^4$ . This validates the paper's use of multistep methods and suggests that discretization error (which Adams-Bashforth reduces from $O(h)$ to $O(h^3)$ ) is a practically significant component of overall estimation error for large systems.
Hyperparameter $\lambda$ sweep (Figure 2): As discussed in Section 5.2, any $\lambda > 0$ works, with values from 1 to 20 (and up to 1,000, per the authors' prose) producing essentially identical parameter estimates and prediction errors. The only effect of larger $\lambda$ is slower convergence (more iterations needed), consistent with the inverse-step-size interpretation. The special case $\lambda = 0$ fails because the proximal term is the only mechanism preventing the algorithm from selecting an arbitrary zero-fidelity (state, parameter) pair with no relationship to the true trajectory.
Noise type robustness (supplementary Figure 6): BCD-prox's performance under Laplacian noise (heavier-tailed than Gaussian) is nearly identical to its performance under Gaussian noise at the same variance, for both FitzHugh–Nagumo and Rössler models. This is a non-trivial finding because many competing methods (EKF, Bayesian methods) make explicit Gaussian assumptions that would be violated by Laplacian noise.
Effect of number of observations $T$ (Figure 4 and supplementary Figure 9): Both BCD-prox and EKF improve as $T$ increases (from $T=20$ to $T=10{,}000$ ), which is expected since more data provides more information. BCD-prox's average estimation error drops from 0.03 ( $T=20$ , $\sigma^2=0.1$ ) to 0.00 ( $T=10{,}000$ , $\sigma^2=0.1$ ), and from 0.47 ( $T=20$ , $\sigma^2=1.5$ ) to 0.06 ( $T=10{,}000$ , $\sigma^2=1.5$ ). The improvement is more substantial for BCD-prox than for EKF, whose high-noise performance remains poor even at large $T$ (average estimation error 3.56 for $T=20$ vs. 1.30 for $T=10{,}000$ at $\sigma^2=1.5$ ).
Effect of model dimensionality (across all experiments): While not presented as a formal ablation, the paper tests ODE systems ranging from $d=2$ (Lotka–Volterra, FitzHugh–Nagumo) to $d=3$ (Rössler) to $d=40$ (Lorenz-96). BCD-prox works across all scales, with the $40$ -dimensional Lorenz-96 system solved to near-zero prediction error using Adams-Bashforth. This provides implicit evidence of scalability, though no systematic study of how runtime scales with $d$ is presented beyond the qualitative claim of ~20 seconds for 400-state problems on a laptop.
Optimization trajectory visualization (supplementary animation, Figure 10 reference): The supplementary material includes an animation showing 250 iterations of BCD-prox on the FitzHugh–Nagumo model, displaying the estimated states $X^*$ , predicted states $\hat{X}$ , and true clean states alongside the evolving parameter estimates. The animation shows qualitatively that the estimated and predicted states converge toward each other and toward the true states as iterations proceed, with parameters simultaneously converging to their true values. This is more of a demonstration than a formal experiment, but it provides visual confirmation that the optimization dynamics behave as described theoretically.

Critical Assessment

The experiments in this paper genuinely support its central claims, but several important caveats and missing pieces should be noted.

The claim of "robustness to noise, parameter initialization, and hyperparameters" is strongly supported, but with a specific, bounded meaning. The initialization-robustness experiment (Figures 3, supplementary Figure 7) convincingly shows that BCD-prox maintains low error when initialized far from the true parameters ( $\sigma^2_\theta = 20$ ), while iPDA fails in roughly half of all trials and lsq/Bayes degrade severely. However, this result is limited to the four tested ODE systems, all of which are relatively low-dimensional nonlinear oscillators/attractors. Whether BCD-prox's initialization robustness extends to ODE systems with other structural properties — stiff systems, systems with multiple time scales, systems where the linear-in-parameters assumption fails — is not tested. The convergence theory requires $f$ to be linear in $\theta$ (or at least to satisfy blockwise strong convexity), and while the experiments demonstrate empirical robustness beyond what the theory covers (FitzHugh–Nagumo requires a change of variables to become linear in parameters, and the experiments work without it), this qualification is not discussed as a limitation.

The hyperparameter-robustness claim is the strongest result in the paper and is fully supported. Figure 2 demonstrates that any $\lambda > 0$ works, with a plateau extending from 1 to at least 1,000. This is a genuine practical advantage over competing methods, whose hyperparameters (smoothing parameters, knot positions, prior distributions, variance settings) must be carefully tuned and cannot be set without access to ground truth. However, the paper's claim that BCD-prox has "only one hyperparameter" is technically true but somewhat incomplete: the user must also choose the discretization order $m$ and the specific multistep method (Adams-Bashforth coefficients), which are architectural choices that affect accuracy. The paper provides guidance (use 3-step Adams-Bashforth) but does not systematically study the sensitivity of results to this choice.

The comparison-with-competitors experiments have genuine methodological strengths. The paired-trial design (same observations and initializations across all methods for each trial) eliminates the confounding effect of easier/harder noise realizations. The 10-trial repetition provides reasonable statistical power to distinguish genuine differences from sampling variability, and the boxplot visualization lets readers assess both central tendency and variance. The use of publicly available implementations for iPDA and Bayes, and the careful tuning of competitor hyperparameters, indicates a fair comparison rather than a straw-man setup.

A significant weakness is the absence of any systematic study of how BCD-prox's performance scales with problem dimensionality or time-series length. The qualitative claim of "~20 seconds for 400-state problems on a laptop" appears without supporting timing experiments, convergence curves across problem sizes, or FLOPs/memory profiling. A reader wondering whether BCD-prox would work on a system with $d=500$ states and $T=10{,}000$ time points has no quantitative basis for extrapolation. The Lorenz-96 ( $d=40$ ) experiment provides proof of concept for moderate-dimensional systems, but the paper does not push further to find the practical limits.

The EKF comparison reveals a genuine advantage of BCD-prox but is somewhat asymmetric. EKF is an online method designed for streaming, real-time filtering; comparing it to BCD-prox (a batch method that iterates repeatedly over the entire dataset) on batch estimation accuracy is fair in the sense that it tests the claim "BCD-prox is a better batch estimator," but it does not address the fact that the methods are designed for different use cases. A practitioner who needs real-time parameter updates as new observations arrive cannot use BCD-prox; a practitioner who has the full time series and wants the best possible estimates should prefer BCD-prox. The paper does not sufficiently discuss this tradeoff — EKF's online nature is framed only as a weakness (it cannot revisit early state estimates), never as a feature (it produces incremental estimates with bounded latency and memory).

The noise-type experiment (supplementary Figure 6) is valuable but limited in scope. Testing Gaussian vs. Laplacian noise demonstrates distributional robustness, but neither distribution captures structured noise patterns that occur in practice — autocorrelated noise, outliers, systematic measurement bias, or time-varying noise variance. The paper's claim that BCD-prox is "distribution-free" is supported only for the two specific distributions tested. No experiment varies the noise model beyond these two symmetric, zero-mean distributions.

Missing experiments that would have strengthened the paper:

Sensitivity to discretization error: The Adams-Bashforth vs. Euler comparison (Figure 1) shows that discretization order matters, but no experiment systematically sweeps $\Delta_i$ (time step) to characterize the relationship between step size and estimation accuracy. The convergence theory guarantees strong convexity only for $\Delta_i < \delta$ (some unknown threshold) — does performance degrade sharply when $\Delta_i$ exceeds this threshold, or gracefully? This is directly relevant to practitioners who cannot arbitrarily choose their observation spacing.
Partial state observations: All experiments assume that all state components are observed (with noise) at all time points. Many real ODE filtering problems involve partial observations (only some components of the state vector are measured). Whether BCD-prox works under partial observability is not tested, and the current formulation (where each observed component directly initializes the corresponding state estimate) would need modification.
Model misspecification: The paper assumes the form of $f$ is known exactly. In practice, the ODE model may be an approximation to the true system — parameters may drift over time, or unmodeled dynamics may be present. Whether BCD-prox's fidelity-driven optimization would overfit to model misspecification (driving $E$ to zero for a wrong model) is not explored.
Comparison with more recent methods: The most recent competing method cited is from 2017; the paper does not compare against neural ODE approaches (Chen et al., 2018) or physics-informed neural networks (Raissi et al., 2019), both of which were emerging at the time and address related problems. This is somewhat understandable given the paper's focus on classical ODE estimation rather than neural-network-based approaches, but it limits the scope of the "state-of-the-art" claim.

A structural strength that should be highlighted: the convergence analysis (Theorems 2–3) covers systems that are genuinely nonlinear and of practical interest — this is not a theory that only applies to toy linear examples. The explicit construction of strong convexity through the $x^+$ / $x^-$ split and the continuity argument for $\Delta_i$ provides a template for analyzing similar block-structured optimization problems in other domains. The fact that all four experimental ODE systems satisfy the theoretical conditions (KL property, linearity in parameters) means the theory and experiments are aligned — this is rarer than it should be in optimization-for-scientific-computing papers.

In summary, the experiments strongly support the paper's specific empirical claims — BCD-prox is more robust to poor initialization, more robust to its single hyperparameter, faster to converge, and more robust to noise than iPDA, Bayesian, least-squares, mean-field, and EKF methods on the tested ODE systems. The experiments do not establish that these advantages generalize beyond the four tested models, beyond symmetric i.i.d. noise, or beyond the fully-observed state setting. The scalability characterization is qualitative rather than quantitative. Nevertheless, within the scope that the paper defines for itself, the experimental evidence is thorough, fair to competing methods, and sufficient to motivate the claimed practical advantages.

6. Limitations and Trade-offs

The Difficulty Estimation Cost Is Not Accounted For

The assumption or constraint. The entire compute-optimal framework — the meta-strategy of selecting different allocation policies for different questions — depends on estimating each prompt's difficulty before applying the chosen strategy. The paper implements difficulty estimation by generating 2048 samples per question and computing either ground-truth correctness (oracle) or the PRM's average final-answer score (predicted). The authors acknowledge this cost explicitly in Section 3.2:

"estimating difficulty in this way still incurs additional computation cost during inference... our experiments do not account for this cost largely for simplicity"

The consequence. The reported 4× efficiency gains over best-of-N are computed after difficulty is known, without amortizing the cost of learning it. In a realistic deployment, the total cost would be difficulty estimation plus strategy execution, and the former could dominate the latter. Generating 2048 samples per question is comparable to or larger than the largest test-time budgets studied (256–512 generations), meaning the 4× figure is best understood as an upper bound on achievable efficiency rather than a realized deployment gain. A practitioner who naively implements the paper's difficulty estimation procedure would likely see negative net efficiency gains for moderate compute budgets.

What evidence exists in the paper. The curves for compute-optimal scaling with predicted difficulty (Figures 4 and 8) are shown without including the difficulty estimation cost in the x-axis (generation budget). The 2048-sample cost is mentioned in Section 3.2 but never plotted or quantified relative to the strategy execution budgets. The paper does not report what fraction of the total compute budget difficulty estimation consumes, nor does it measure how the reported gains change when this cost is included.

Mitigation status. The paper flags this as a key area for future work (Section 8), suggesting "training models to predict difficulty directly from the question text." No such model is developed or evaluated. The predicted difficulty bins (which use the PRM's average score rather than ground-truth answers) remove the need for labels but do not reduce the sample cost. The limitation is acknowledged transparently but remains essentially unaddressed; the paper's headline efficiency numbers should be interpreted with this caveat in mind.

Hard Problems Remain Fundamentally Unsolved

The assumption or constraint. Test-time compute amplifies a model's existing capability — it improves the probability of finding a correct answer that the model can already generate at some non-trivial rate. For problems where the base model's pass@1 is near zero, no amount of search or revision helps because there are no correct solutions in the proposal distribution to find or refine.

The consequence. Across all methods — search, revisions, and their compute-optimal combinations — the hardest questions (difficulty bin 5) show near-zero improvement regardless of compute budget. In Figure 3 (right), bin 5 accuracy hovers at 1–3% for all search methods and all budget levels. In Figure 7 (right), bin 5 shows roughly 2–3% accuracy irrespective of the sequential-to-parallel ratio. In the FLOPs-matched comparison (Figure 9), the bin 5 scaling line is essentially flat near 0–5% — far below the ~14× larger model's performance. The paper is candid about this (Section 7), but the practical implication is that test-time compute cannot substitute for pretraining on genuinely novel or out-of-distribution reasoning. For such problems, scaling model size and training data remains the only viable path.

What evidence exists in the paper. The difficulty-bin analyses (Figures 3 right, 7 right) provide clear evidence: the hardest quintile shows near-zero accuracy and zero slope (no improvement with additional compute). The FLOPs-matched analysis (Section 7, Figure 9) quantifies this further — on hard problems, the smaller model with compute-optimal test-time compute shows relative disadvantages of −37.2% to −52.9% compared to the ~14× larger model, depending on $R$ and method.

Mitigation status. The paper does not attempt to solve this. The limitation is inherent to the approach: test-time compute can amplify existing capability but cannot create it from nothing. The paper is transparent about this boundary condition (Section 7, item 2 in the takeaway box), which is to its credit. Practitioners should understand that the framework is useful only for problems within their base model's approximate capability range.

Results Are Limited to a Single Benchmark and Model Family

The assumption or constraint. All experiments use the MATH benchmark (500 test questions) with PaLM 2-S* as the base model. The authors state they "believe this model is representative of the capabilities of many contemporary LLMs" (Section 4), but this claim is unverified. Several aspects of the findings might be model-specific:

The PRM's quality and over-optimization behavior depend on PaLM 2-S*'s output distribution. A model with different calibration properties or different error patterns might exhibit different difficulty-dependent scaling curves.
The revision model's ability to learn from incorrect in-context examples depends on the base model's in-context learning capabilities, which vary substantially across model families.
The MATH benchmark consists exclusively of competition-level math problems requiring symbolic reasoning. Whether the difficulty-dependent patterns (beam search hurting easy problems, revisions helping easy problems) generalize to other reasoning domains (code generation, logical reasoning, scientific QA) or to tasks requiring factual knowledge rather than inference is unknown.

The consequence. A practitioner seeking to apply these techniques to a different model family (e.g., LLaMA, GPT) or a different domain (e.g., code generation, open-ended QA) has no direct evidence that the reported 4× efficiency gains, the optimal allocation policies, or the difficulty-dependent strategy rankings will transfer. The optimal beam width ( $M=4$ vs. $M=\sqrt{N}$ ), the optimal sequential-to-parallel ratio, and even whether revisions outperform search might all change with a different base model or task distribution.

What evidence exists in the paper. The paper reports no experiments on any benchmark other than MATH (Figure references throughout Section 5) and no experiments with any base model other than PaLM 2-S*. The supplementary material confirms this: all model-specific details (Appendix D, F, H) reference only PaLM 2 series architectures.

Mitigation status. The paper does not address this limitation. No cross-model or cross-domain experiments are performed or proposed as future work. The claim of representativeness (Section 4) is an assertion, not an empirically supported statement. This is a significant gap for a paper whose primary contribution is an empirical characterization of test-time compute scaling behavior — the characterization might be specific to the (model, benchmark) pair studied.

The `~14×` Larger Model Baseline Is Not Compute-Optimally Trained

The assumption or constraint. The FLOPs-matched comparison in Section 7 scales model parameters while holding training data fixed, following the LLaMA paradigm (Touvron et al., 2023). The authors acknowledge that this departs from compute-optimal pretraining (Hoffmann et al., 2022), where both data and parameters are scaled equally:

"We choose this setting as it is representative of a canonical approach to scaling pretraining compute and leave the analysis of compute-optimal scaling of pretraining compute where the data and parameters are both scaled equally to future work."

Additionally, the ~14× larger model uses only greedy decoding — no majority voting, no best-of-N, no test-time compute augmentation of any kind.

The consequence. A Chinchilla-optimal model trained with ~14× more total FLOPs (scaling both parameters and data) would likely outperform a parameter-only-scaled model, making the pretraining baseline weaker than it could be. The reported advantages of test-time compute over pretraining (e.g., +27.8% on easy questions at $R \ll 1$ , from Figure 1) might shrink or reverse against a properly compute-optimal larger model. Similarly, giving the larger model even a modest test-time compute budget (say, best-of-8) would create a much stronger baseline. The paper's headline finding — that test-time compute can substitute for pretraining — is established only against this specific, relatively weak pretraining baseline.

What evidence exists in the paper. The FLOPs-matched results appear in Section 7 and Figure 9 (also summarized in the Figure 1 bar charts). The baseline model's training details are described only as having ~14× more parameters with fixed data (Section 7). No comparison is made against a Chinchilla-optimal larger model, and no experiment gives the larger model any test-time compute budget.

Mitigation status. The paper is transparent about this caveat (Section 7, quoted above), which is to its credit. The choice is justified as representative of one common scaling paradigm (LLaMA-style). However, the limitation means that the paper's central claim about the pretraining-inference tradeoff — perhaps its most impactful result — rests on a baseline that may not represent the strongest possible use of additional pretraining compute. The authors explicitly defer the compute-optimal comparison to future work.

Latency and Wall-Clock Time Are Ignored

The assumption or constraint. The paper measures test-time compute exclusively in "generations" — the number of complete solutions sampled. This is a reasonable proxy for total FLOPs but ignores latency: the wall-clock time required to produce an answer. Sequential revisions are inherently serial — each revision step depends on the output of the previous one — while parallel best-of-N can be executed simultaneously with sufficient hardware.

The consequence. A strategy that allocates 128 generations as 64 sequential × 2 parallel (as might be optimal for certain difficulty bins in Figure 7) takes roughly 64× longer wall-clock time than one that runs 128 parallel samples simultaneously. For latency-sensitive applications — interactive assistants, real-time decision-making, online tutoring — the sequential-heavy strategies favored by the compute-optimal policy on easy problems may be impractical regardless of their FLOPs-efficiency advantages. The paper's framework implicitly assumes that total FLOPs is the resource to optimize, but in many deployment scenarios, latency is the binding constraint.

What evidence exists in the paper. The paper never reports wall-clock times, latency distributions, or any metric other than generation count and accuracy. The sequential-to-parallel ratio analysis (Figure 7) shows what allocation is optimal in terms of generations, but does not discuss the latency implications of high-sequential ratios. The paper's training time claims ("~20 seconds on a standard laptop") appear in the introduction but refer to the underlying optimization algorithm's iterations, not to end-to-end inference latency for a single query.

Mitigation status. Not addressed. The paper does not discuss the latency-throughput tradeoff or suggest how practitioners should incorporate latency constraints into the allocation policy. For a method that explicitly targets deployment decisions (Section 1 mentions "on-device deployment" as a motivation), this is a notable gap. A latency-aware extension would need to account for the fact that parallel generations can be batched on GPUs while sequential revisions cannot, changing the optimal allocation in ways the current framework does not capture.

Search and Revisions Are Studied Independently, Not Combined

The assumption or constraint. The paper studies two complementary axes — PRM-guided search (Section 5) and iterative revisions (Section 6) — as independent mechanisms for improving test-time compute efficiency. The compute-optimal policies for each axis are derived separately, and the FLOPs-matched comparison evaluates them individually against the pretraining baseline. The two mechanisms are never combined into a single system.

The consequence. The two mechanisms have complementary strengths: revisions improve the proposal distribution (generating better candidates), while PRM search improves candidate selection (finding the best among generated candidates). Applying beam search to revision model outputs — or using the PRM to guide which revisions to pursue rather than generating a blind sequential chain — could yield gains beyond either method alone. The 4× efficiency improvements reported for each axis individually therefore represent a lower bound on what a fully integrated system could achieve. Conversely, combining the methods might introduce new failure modes (e.g., the PRM over-optimizing on revision-model outputs differently than on base-model outputs, as already noted in Appendix J, Figure 15a) that the paper cannot characterize.

What evidence exists in the paper. Section 8 explicitly acknowledges this gap:

"we did not experiment with PRM tree-search techniques in combination with revisions"

The experiments in Sections 5 and 6 keep the two approaches entirely separate. Appendix J (Figure 15a) shows that the base-model PRM underperforms when scoring revision model outputs due to distribution shift, suggesting that combining the methods would require retraining the verifier on revision-model data — but this combined training is never attempted.

Mitigation status. Acknowledged as future work (Section 8) but not addressed. This is arguably the most natural extension of the paper's framework, since the unifying proposal-verifier decomposition (Section 2) explicitly suggests that both mechanisms can be combined. The current results, while strong individually, leave open the question of whether the complementary strengths of search and revisions translate to additive or super-additive gains in practice — and whether the difficulty-dependent optimal policy changes when both mechanisms are available simultaneously.

7. Implications and Future Directions

How This Work Changes the Landscape

This paper introduces a genuinely new algorithmic primitive for simultaneous filtering and parameter estimation — the combination of block coordinate descent with a moving proximal anchor that updates at each outer iteration rather than remaining tethered to the original noisy observations. This is not an incremental refinement of existing methods; it is a conceptual reframing of what regularization means in the context of ODE parameter estimation.

The dominant paradigm prior to this work — exemplified by iPDA and spline-based methods — treated filtering as a problem of smoothing noisy observations toward a latent truth, with the observations as the fixed reference point. The regularization term encoded "stay close to the data" as a hard design principle. BCD-prox replaces this static anchoring with a dynamic one: the proximal term says "don't change too much from where you were last iteration," and since where you were last iteration has itself moved away from the noisy data, the algorithm can progressively escape bad initial observations. The paper demonstrates that this single architectural change is responsible for the dramatic robustness improvements — iPDA fails in nearly half of all trials under poor initialization, while BCD-prox fails in only 4 out of 80 (Section 5.3, Figures 3 and 7).

This reframing resolves a tension that had been latent in the literature. Prior methods faced an impossible hyperparameter tuning problem: the optimal balance between data fidelity and ODE fidelity depends on the unknown noise magnitude and the unknown vector field, making it impossible to set correctly in practice without access to ground truth (which is precisely what estimation is supposed to provide). BCD-prox dissolves this problem by making the hyperparameter control step size rather than data fidelity. The paper's demonstration that any $\lambda > 0$ works, with values from 1 to 1,000 producing identical final estimates (Figure 2), is not merely a convenient empirical property — it is evidence that the algorithm has correctly separated the optimization dynamics (where step-size control matters) from the estimation objective (where the right answer is determined by ODE consistency plus the optimization trajectory, not by an arbitrary penalty weight). This is a methodological insight that transfers beyond ODE estimation to any problem where fidelity to first-principles constraints must be balanced against fidelity to noisy observations.

The work also changes the practical landscape by establishing that batch optimization can be fast enough to be competitive with online methods while providing substantially better estimates. The finding that BCD-prox takes ~20 seconds on a standard laptop for 400-state problems, compared to minutes for spline-based methods and hours for Bayesian methods (Section 5), challenges the implicit assumption that sophisticated batch filtering is too computationally expensive for routine use. The EKF comparison (Figure 4, supplementary Figure 9) is particularly instructive: EKF, an online method designed for speed, is competitive with BCD-prox only in the idealized regime of very long time series with low noise ( $T = 10{,}000$ , $\sigma^2 = 0.1$ ). In the more realistic regime of shorter, noisier observations, BCD-prox outperforms EKF by factors of 7–22× in estimation error. This suggests that for many scientific applications where the full time series is available for offline analysis, batch methods like BCD-prox should be the default choice — the historical preference for online filtering was driven more by computational constraints that no longer bind than by intrinsic statistical advantages.

The paper's convergence theory (Theorems 2 and 3, Section 4) makes a methodological contribution that has implications beyond ODE estimation. The construction of blockwise strong convexity through the $x^+$ / $x^-$ split — proving positive definiteness at $\Delta_i = 0$ analytically, then extending to finite $\Delta_i$ via continuity of eigenvalues — provides a template for analyzing block-structured optimization problems where the coupling between blocks is mediated by a small parameter (here, the time step). This is a concrete, verifiable approach to establishing convergence guarantees for alternating minimization on non-convex problems, one that does not rely on abstract Lipschitz or KL assumptions that are difficult to check for specific systems. The fact that the theoretical conditions (linearity in parameters, KL property) are satisfied by all four experimentally tested ODE systems — including nonlinear, chaotic systems like Rössler and Lorenz-96 — demonstrates that the theory covers problems of genuine scientific interest, not just toy examples.

The work also makes certain research directions less attractive. The repeated finding that BCD-prox outperforms Bayesian methods (Figures 3, 4, 7) despite the Bayesian framework's theoretical elegance suggests that for the specific problem of simultaneous filtering and parameter estimation with known ODE structure, the computational overhead and hyperparameter sensitivity of full Bayesian inference are not justified by improved accuracy. The paper's demonstration that the mean-field variational Bayes method of Gorbach et al. (2017) — which was explicitly designed to overcome the limitations of earlier Bayesian approaches — still underperforms BCD-prox, particularly on parameter estimation (Figure 4, first row, where mean-field's average error on $\theta_0$ is 7.41 vs. BCD-prox's 0.54), strengthens this conclusion. Future work on Bayesian approaches to this problem would need to demonstrate substantial accuracy improvements over BCD-prox to justify their additional complexity and computational cost.

Follow-Up Research This Work Enables

Extending BCD-prox to partially observed state vectors. The paper assumes that all $d$ components of the state vector are observed (with noise) at all $T$ time points. In many real applications — biological systems where only certain molecular concentrations can be measured, engineering systems where sensors cover only a subset of state variables — this assumption fails. The current formulation would need modification because each observed component directly initializes its corresponding state estimate in $X^{*(0)} = Y$ . For unobserved components, there is no natural initialization. A natural extension would introduce a missing-data mechanism: initialize unobserved states via forward integration from the nearest observed states using an initial parameter guess, then fold the uncertainty about these imputed values into the proximal penalty (perhaps with component-specific $\lambda$ weights that are larger for imputed components, reflecting lower confidence). A strong evaluation would test BCD-prox on the same ODE systems used in this paper but with only a random subset (e.g., 50% or 25%) of state components observed, measuring how estimation error degrades as a function of observation completeness and comparing against methods specifically designed for partially observed systems (e.g., Kalman filters, which naturally handle partial observations through the observation matrix $H$ ).

Combining BCD-prox with automated discovery of the vector field $f$ . The paper assumes the functional form of $f(x, \theta)$ is known — only the parameters $\theta$ are unknown. The authors mention in their conclusion that extending BCD-prox "to estimate the vector field $f$ from noisy observations" is ongoing/future work. This is a natural and ambitious extension: replace the known $f$ with a parameterized function approximator (e.g., a neural network, a sparse linear combination of basis functions from a candidate library as in SINDy (Brunton et al., 2016), or a symbolic regression model) and use BCD-prox's alternating optimization framework to jointly estimate the states, the parameters of $f$ , and potentially the structure of $f$ itself. The key challenge is that the fidelity objective $E$ would no longer be linear in the parameters governing $f$ , potentially breaking the strong convexity guarantees that underpin the convergence theory. A strong follow-up would test this on systems where the true $f$ is known (so accuracy can be measured) but withheld from the algorithm, comparing the recovered vector field against the ground truth for systems of increasing complexity (linear $\to$ polynomial nonlinear $\to$ rational nonlinear $\to$ chaotic). The Lorenz-96 model with $d=40$ would be a particularly stringent test, since a 40-dimensional learned vector field must capture the correct coupling structure to produce accurate long-term predictions.

Systematic characterization of the $\Delta_i$ threshold for convergence. The convergence theory (Section 4) proves the existence of $\delta > 0$ such that blockwise strong convexity holds for $\Delta_i \in (0, \delta)$ , but does not provide a way to compute $\delta$ for a given ODE system. A practitioner who wants to know whether their observation spacing is small enough to guarantee convergence has no guidance. A valuable follow-up would empirically map the relationship between $\Delta_i$ , convergence behavior, and estimation accuracy for the ODE systems studied in this paper. The experiment would systematically vary $\Delta_i$ (or equivalently, vary $T$ for a fixed time range) and measure: (a) the number of iterations to convergence, (b) whether the algorithm converges at all or diverges, (c) the final parameter and prediction errors for convergent runs, and (d) the spectral properties of the block Hessians to empirically identify where strong convexity is lost. This would produce practical guidance — e.g., "for Lorenz-96 with $d=40$ , BCD-prox requires $\Delta_i < 0.05$ for reliable convergence" — and would validate whether the $\Delta_i \to 0$ limit in the theory is a practically binding constraint or merely a proof artifact. It would also reveal whether higher-order discretizations (Adams-Bashforth, $m=3$ ) relax the $\Delta_i$ constraint compared to Euler ( $m=1$ ), which the paper's Figure 1 hints at but does not systematically test.

BCD-prox as a component in iterative experimental design. In many applications, observations are expensive (e.g., biological experiments, field measurements), and the goal is not just to estimate parameters from existing data but to choose which observations to collect next to maximize information gain. BCD-prox's ability to produce both state estimates and parameter estimates from data — and to do so robustly even with limited, noisy observations — makes it a natural inner loop for an adaptive experimental design framework. The specific direction: given a current (state, parameter) estimate from BCD-prox and an estimate of uncertainty (perhaps obtained by running BCD-prox on bootstrap resamples of the data, or by examining the spectrum of the Hessian at convergence), select the next observation time $t_{T+1}$ and/or the next state component to measure that maximizes expected reduction in parameter uncertainty. A strong evaluation would compare this adaptive strategy against uniform sampling (collecting observations at evenly spaced times) on the Lotka–Volterra and FitzHugh–Nagumo models, measuring how many adaptive observations are needed to achieve the same parameter estimation accuracy as a fixed grid of $T=100$ uniform observations. This connects BCD-prox to the optimal experimental design literature (which often uses much simpler estimation procedures) and could substantially reduce data collection costs in applications where each observation is expensive.

Theoretical analysis of the proximal anchor update schedule. The current algorithm updates the proximal anchor every outer iteration: $X^{*(n)}$ becomes the anchor for iteration $n+1$ . This is the simplest possible schedule, but there is no theoretical justification that it is optimal. The proximal term's role as inverse step size (Equation 12) suggests connections to the step-size scheduling literature in optimization: decreasing step sizes over time (i.e., increasing $\lambda$ ) can improve convergence rates for convex problems, while adaptive step sizes (e.g., Barzilai-Borwein) can accelerate convergence for some non-convex problems. A follow-up could explore alternative anchor-update schedules: (a) updating the anchor only every $k$ iterations (allowing more aggressive state movement between anchor updates), (b) using a weighted combination of multiple previous iterates as the anchor (momentum-like), or (c) adaptively adjusting $\lambda$ based on the observed reduction in $E$ (increasing $\lambda$ when progress stalls, decreasing it when progress is rapid). A strong evaluation would measure the number of iterations and wall-clock time to convergence under each schedule for the same ODE systems and noise levels used in the paper, with the goal of finding schedules that converge in fewer iterations than the default $\lambda=1$ fixed schedule without sacrificing final accuracy. Negative results (schedules that cause divergence or oscillation) would also be informative for understanding the stability properties of the proximal mechanism.

Stress-testing robustness to model misspecification. The paper assumes the ODE model is correctly specified — the true dynamics follow the given $f(x, \theta)$ exactly. In practice, models are approximations: they omit higher-order effects, approximate nonlinearities, or assume constant parameters that actually drift. A critical follow-up would test how BCD-prox behaves when the model is misspecified. For example: generate data from the FitzHugh–Nagumo model with a small cubic perturbation not present in the assumed $f$ , run BCD-prox assuming the standard (unperturbed) model, and measure both parameter estimation error and the fidelity $E$ at convergence. If BCD-prox's fidelity-driven optimization overfits to the misspecified model — driving $E$ to a small value by finding parameters that are very different from the true ones — this would reveal a brittleness that the current experiments cannot detect. Conversely, if the proximal term's regularization prevents extreme parameter values and BCD-prox converges to parameters that are the "best fit" in some meaningful sense (e.g., minimizing the true prediction error even though $E$ cannot be driven to zero), this would demonstrate a robustness to misspecification that is not claimed in the paper but would be practically valuable. Comparing BCD-prox's behavior under misspecification against Bayesian methods (which naturally quantify model uncertainty) and iPDA (which may also overfit) would help practitioners choose methods based on how certain they are about their model structure.

Practical Applications and Downstream Use Cases

Pharmacokinetic/pharmacodynamic (PK/PD) modeling from sparse, noisy clinical data. PK/PD models describe how drug concentrations in the body change over time (typically using systems of 2–5 nonlinear ODEs) and are used extensively in drug development to determine dosing regimens. Clinical data is characteristically sparse (blood samples drawn at a handful of time points per patient), noisy (measurement error in concentration assays), and heterogeneous across patients. Current practice often relies on nonlinear least squares (which the paper shows diverges under poor initialization) or Bayesian methods (which the paper shows are slow and hyperparameter-sensitive). BCD-prox's demonstrated ability to produce accurate parameter estimates from as few as $T=20$ noisy observations (supplementary Figure 9, where BCD-prox achieves average estimation error 0.47 at $\sigma^2 = 1.5$ vs. EKF's 3.56) makes it directly applicable to this setting. The single hyperparameter $\lambda$ (which can be set to any positive value without tuning, per Figure 2) eliminates the need for the extensive prior sensitivity analyses that Bayesian PK/PD models require. The ~20-second runtime on a laptop means that a clinician could fit patient-specific parameters at the point of care rather than waiting for offline batch processing.

Calibration of building energy models from sensor data. Modern buildings are instrumented with sensors measuring temperature, humidity, and energy consumption at multiple zones, governed by thermal dynamics modeled as coupled ODEs (often 10–50 state variables representing zone temperatures). These models contain uncertain parameters (thermal resistances, heat capacities, HVAC efficiencies) that must be estimated from noisy sensor data for model-predictive control to work effectively. The paper's demonstration that BCD-prox scales to $d=40$ dimensions (Lorenz-96, Figure 1) and handles Laplacian noise as well as Gaussian noise (supplementary Figure 6) is directly relevant: building sensor data often contains outlier readings from sensor faults or temporary anomalies, producing heavier-tailed noise distributions than Gaussian models assume. BCD-prox's distribution-free, batch nature means it can process a day's worth of sensor data (~ $T = 1{,}440$ readings at 1-minute intervals for 40 zones = 57,600 data points) in seconds on a standard laptop, producing a calibrated model overnight for next-day control optimization. The robustness to parameter initialization (Figures 3 and 7) is practically important because initial parameter guesses for building models — based on architectural plans and material specifications — can deviate substantially from as-built reality.

Fitting ecological population models from field survey data. The Lotka–Volterra model studied in the paper is the canonical predator-prey model in ecology, but real ecological applications involve extensions with multiple interacting species, seasonal forcing, and spatial coupling — producing ODE systems with 10–100+ state variables. Field survey data is notoriously noisy (observation error from sampling variability, missed detections, and migration) and sparse (surveys conducted monthly or seasonally, not continuously). BCD-prox's robustness to both noise magnitude (Figures 4 and supplementary Figure 8 show graceful degradation as $\sigma^2$ increases from 0.5 to 1.5) and noise type (supplementary Figure 6 shows identical performance under Gaussian and heavy-tailed Laplacian noise) makes it suitable for ecological data, which rarely satisfies Gaussian assumptions. The paper's finding that BCD-prox achieves near-zero estimation error with sufficiently long time series (supplementary Figure 9: average estimation error $\approx 0$ for $T = 10{,}000$ ) suggests that for long-term ecological monitoring datasets (spanning decades of monthly surveys), BCD-prox could recover interaction parameters (predation rates, competition coefficients) with high precision, enabling quantitative tests of ecological theory that are currently hampered by estimation uncertainty.

Rapid parameter estimation for in-situ battery state-of-health monitoring. Lithium-ion battery degradation is modeled by coupled electrochemical ODEs (typically 5–15 state variables representing lithium concentrations, potentials, and temperatures) with parameters that drift over the battery's lifetime (capacity, internal resistance, diffusion coefficients). Estimating these parameters from voltage and current measurements during operation enables state-of-health monitoring and remaining-useful-life prediction for electric vehicles and grid storage. Current methods often use extended Kalman filters, which the paper shows underperform BCD-prox by 7–22× in estimation error under realistic noise (Figure 4, $\sigma^2=1.5$ , $T=10{,}000$ : BCD-prox estimation error 0.06 vs. EKF 1.30). Moreover, EKF's inability to update the initial state means that errors in the initial state-of-charge estimate (which is typically uncertain at startup) propagate through all subsequent predictions. BCD-prox's batch processing over a sliding window of recent measurements (e.g., the last $T=500$ voltage/current readings at 1 Hz) would produce updated parameter and state estimates every ~8 minutes, with the proximal mechanism ensuring that estimates evolve smoothly as the window advances. The single, non-critical hyperparameter means the algorithm can be deployed without per-vehicle tuning.

When to Prefer This Method

The paper articulates clear tradeoffs against specific named alternatives (iPDA, Bayesian methods, EKF, mean-field variational Bayes) and identifies the conditions under which each method fails. These conditions directly imply decision rules for practitioners:

Prefer BCD-prox when the observation noise magnitude is unknown or potentially large ( $\sigma^2 \geq 0.5$ in the paper's tested range), because BCD-prox's moving proximal anchor allows escape from noisy initial data while iPDA's static regularization traps filtered states near bad observations. This is the central failure mode the paper documents: iPDA produces prediction errors exceeding 100 in 39 out of 80 trials under poor initialization (Section 5.3), while BCD-prox fails in only 4 out of 80.
Prefer BCD-prox when any $\lambda > 0$ works and producing identical results from 1 to at least 1,000 (Section 5.2, Figure 2), compared to the large number of sensitive hyperparameters in spline-based methods (smoothing parameters, knot positions/numbers) and Bayesian methods (prior distributions, variances, kernel widths). A practitioner can set $\lambda = 1$ and obtain good results without cross-validation or domain expertise in hyperparameter tuning.
Prefer BCD-prox when training time matters: the paper reports ~20 seconds on a standard laptop for 400-state problems (Section 5), compared to minutes for spline-based methods and hours for Bayesian methods. This matters for applications requiring repeated estimation (e.g., sliding-window batch processing, multi-start initializations, or fitting models to data from many experimental subjects).
Prefer BCD-prox when the ODE model is known to be correct (or nearly so): BCD-prox's fidelity term $E$ drives toward exact ODE consistency, which is appropriate when the model is trusted. For misspecified models, BCD-prox's behavior is untested — Bayesian methods, which naturally quantify model uncertainty, may be preferable despite their computational cost.
Prefer EKF when real-time, online estimation is required with strict latency constraints: EKF processes observations sequentially with bounded computation per time step, while BCD-prox iterates over the entire dataset repeatedly. For streaming applications where decisions must be made as each new observation arrives, BCD-prox's batch nature is disqualifying regardless of its accuracy advantages.
Prefer EKF when the observation noise is very low ( $\sigma^2 = 0.1$ in the paper's tests) and the time series is very long ( $T = 10{,}000$ ): in this regime, EKF performs comparably to BCD-prox (Figure 4, second row, where estimation errors are 0.09 vs. 0.00) and has the advantage of online processing and lower memory requirements. The paper's results show that EKF's disadvantages emerge primarily under higher noise and shorter time series.
Prefer iPDA or spline-based methods when the ODE vector field $f$ is only approximately known and human-interpretable smoothness of the filtered states is valued over exact ODE consistency: iPDA's spline representation produces visually smooth state trajectories that may be preferable for exploratory data analysis even if parameter estimates are less accurate. BCD-prox learns states directly without smoothness constraints, which can produce trajectories that satisfy the ODE exactly but appear less smooth than spline fits — a tradeoff between physical consistency and visual interpretability that the paper does not discuss but that matters in practice.

🎯 Pitch

1. Executive Summary

2. Context and Motivation

The Core Problem: Learning ODE States and Parameters from Noisy Data

Why Existing Approaches Struggle

The Deeper Structural Insight: Multiple Global Minima

How BCD-prox Addresses These Gaps

Positioning Relative to Prior Work

3. Technical Approach

3.1 Reader Orientation

3.2 Big-Picture Architecture (Diagram in Words)

3.3 Roadmap for the Deep Dive

3.4 Detailed, Sentence-Based Technical Breakdown

The ODE Discretization: From Continuous Dynamics to Algebraic Constraints

The Fidelity Objective: Measuring ODE Consistency

The Proximal Term and the Full Per-Iteration Objective

Block Coordinate Descent: Alternating Parameter and State Optimization

Implementation Details: LBFGS and Automatic Differentiation

Stopping Criterion and the Three Termination Cases

Block Splitting for Convergence Theory: The x+x^+x+ / x−x^-x− Decomposition

Final Predicted States: The Guaranteed Zero-Fidelity Output

4. Key Insights and Innovations

Innovation 1: The Proximal Anchor as a Moving Target Rather Than a Fixed Regularization

Innovation 2: The Infinite-Minima Problem as an Explicit Design Constraint, Not a Pathological Edge Case

Innovation 3: The Convergence Theory via Blockwise Strong Convexity with an Explicit Δi\Delta_iΔi​ Threshold

Innovation 4: The Predicted-State Post-Processing Step as an Implicit Projection onto the ODE-Consistent Manifold

5. Experimental Analysis

Evaluation Methodology

Main Quantitative Results

5.1 Fidelity Optimization Decreases Prediction Error Across Iterations (Figure 1, Supplementary Figure 5)

5.2 Robustness to the Single Hyperparameter λ\lambdaλ (Figure 2)

5.3 Comparison with iPDA, Bayesian, and Least-Squares Methods Under Varying Parameter Initialization (Figure 3, Supplementary Figure 7)

5.4 Comparison with Mean-Field Variational Bayes (Figure 4, First Row; Supplementary Figure 8)

5.5 Comparison with Extended Kalman Filter (Figure 4, Second and Third Rows; Supplementary Figure 9)

5.6 Effect of Noise Type and Magnitude (Supplementary Figure 6)

Ablation Studies and Robustness Checks

Critical Assessment

6. Limitations and Trade-offs

The Difficulty Estimation Cost Is Not Accounted For

Hard Problems Remain Fundamentally Unsolved

Results Are Limited to a Single Benchmark and Model Family

The ~14× Larger Model Baseline Is Not Compute-Optimally Trained

Latency and Wall-Clock Time Are Ignored

Search and Revisions Are Studied Independently, Not Combined

7. Implications and Future Directions

How This Work Changes the Landscape

Follow-Up Research This Work Enables

Practical Applications and Downstream Use Cases

When to Prefer This Method

Adjacent in domain

GPipe: Efficient Training of Giant Neural Networks using Pipeline Parallelism

Large Scale GAN Training for High Fidelity Natural Image Synthesis

Block Splitting for Convergence Theory: The $x^+$ / $x^-$ Decomposition

Innovation 3: The Convergence Theory via Blockwise Strong Convexity with an Explicit $\Delta_i$ Threshold

5.2 Robustness to the Single Hyperparameter $\lambda$ (Figure 2)

The `~14×` Larger Model Baseline Is Not Compute-Optimally Trained