{"id":"4ef8f1b9-a406-49ec-af29-167c3b893b93","arxiv_id":"2506.21776","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"rodeo is a JAX-based Python library that implements probabilistic ODE solvers and several Bayesian parameter inference methods with linear scaling in system size.","lead":"This paper introduces rodeo, a Python library for probabilistic ODE solving and parameter inference built on JAX. It combines several existing methods and demonstrates fast, accurate inference on biological and epidemiological models.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The linear-in-d claim rests on the untested block-diagonal Jacobian approximation; for strongly coupled ODEs this may lose accuracy at fixed N, undermining the headline speed-accuracy trade-off.","rationale":"The central claim has two parts: the linear O(d) solver and fast, accurate parameter inference. The O(d) claim is true only under the block-diagonal structure of Section 2.3, and the interrogation method that preserves that structure while remaining first-order accurate is Krämer et al.'s J*_f, with off-diagonal Jacobian blocks set to zero. The reader's weakest assumption identifies exactly this: no quantitative accuracy comparison against the full Jacobian is reported. This is load-bearing because the block structure is not a cosmetic optimization; it changes the filter's posterior geometry (block-diagonal covariance, no cross-variable linear feedback), so the approximation's validity depends on the strength of coupling and step size. The paper does provide indirect support: in Sections 5.8, 5.10, and 5.12, Krämer-interrogation posteriors match RKDP posteriors for FN, Hes1, and SEIRAH at sufficiently small dt, and Table 3 shows real speed-ups. Those examples are smooth and moderately coupled, however, so they do not demonstrate that 'minimal loss in accuracy' holds broadly. This justifies keeping the reader's CONDITIONAL verdict rather than ACCEPT or REJECT: the library is credible and the examples support the main functionality, but the headline scaling claim should be qualified by the J* approximation or backed by a direct head-to-head test. The proposed linear-system experiment with tunable off-diagonal coupling is a minimal, decisive check of whether this concern lands.","tokens_in":45420,"tokens_out":13310,"duration_ms":140793,"concrete_test":"Run a two-variable linear ODE x' = A x with A = [[-1, κ], [κ, -1]] and κ ranging from 0 to 100, using the same IBM prior, N = 100, and a Gaussian observation model as in Section 5.2. Compare rodeo with Krämer interrogation against rodeo with the full Tronarp Jacobian (or a dense Kalman implementation) and an RKDP reference: record the maximum error of the posterior mean µ_{0:N|N} and the empirical coverage of 95% credible intervals for the state. If error grows materially with κ or coverage collapses below nominal, the unconditional O(d) claim needs an explicit accuracy caveat for strongly coupled systems.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 2.3 obtains O(d) complexity by requiring block-diagonal W and, for the interrogation used in the benchmark examples, Krämer et al.'s modified Jacobian J*_f, which zeros every off-diagonal block of ∂f/∂X. The paper asserts 'extensive evidence ... minimal loss in accuracy' but reports no direct comparison against the full Jacobian. The consequence is structural: the Kalman update for variable k sees no linearized sensitivity to any other variable, so the posterior covariance stays block-diagonal and cross-variable uncertainty is discarded. For strongly coupled systems (large off-diagonal Jacobian entries), the filter mean can drift from the full-linearization solution at any fixed N, and the speed-accuracy trade-off in Table 3, obtained at the smallest N that visually matches RKDP, would not transfer. The paper's examples are smooth and moderately coupled, so they do not stress this assumption.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents rodeo, a Python library built on JAX for probabilistic parameter inference in ODEs. The library implements a Bayesian filtering framework with Gaussian Markov process priors, several model interrogation methods (Chkrebtii et al., Schober et al., Tronarp et al., Krämer et al.), and multiple inference algorithms: a basic plug-in likelihood, Fenrir, DALTON (Gaussian and non-Gaussian), marginal MCMC, and a Markov-prior variant of MAGI. The central claims are that the solver scales linearly in the number of discretization points N and system variables d (via block-diagonal interrogation and a modified Jacobian), and that inference is fast, accurate, and scalable, as demonstrated on four examples (a second-order univariate ODE, FitzHugh-Nagumo, Hes1, and SEIRAH) with comparisons to LSODA and RKDP. The paper is written primarily as a software and methods showcase, with extensive code listings and pseudocode for each algorithm.","tokens_in":45587,"tokens_out":10596,"duration_ms":101816,"significance":"If the claims hold, rodeo is a valuable contribution: it provides a JAX-native, AD-compatible, JIT-compiled implementation of several published probabilistic ODE inference methods in one package, making them accessible to Python users. The paper includes detailed pseudocode, runnable examples, and benchmarks against independent deterministic solvers (LSODA and RKDP via difrax), which are welcome strengths. The generalization of MAGI to arbitrary-order ODEs with Markov priors is a useful extension. However, the headline \"linear scaling\" and \"fast\" claims rest on a block-diagonal Jacobian approximation whose accuracy is asserted but not quantitatively demonstrated, and the speed summary in Section 5.13 is not consistent with Table 3. These issues are fixable but should be addressed before the paper can be recommended for publication.","major_comments":[{"comment":"The linear-in-d complexity claim depends on the block-diagonal modified Jacobian J*_f (Krämer et al. 2021), which zeroes all off-diagonal blocks of the Jacobian. The paper states that \"there is extensive evidence in our experiments and those of Krämer et al. (2021) suggesting there is minimal loss in accuracy compared to using the full Jacobian,\" but no quantitative comparison is reported in this manuscript. This is load-bearing because the speed-accuracy trade-off in the numerical examples is obtained only under this approximation. Please add a numerical study on a strongly coupled ODE (e.g., with large off-diagonal Jacobian entries) comparing the posterior produced by J* versus the full Jacobian at fixed N, or explicitly state the limitation that the linear-scaling claim holds only when the block-diagonal approximation is accurate.","section":"§2.3"},{"comment":"The text says the last three examples show rodeo is \"2 to 12 times faster than RKDP and 2-12 times faster than LSODA,\" but the RKDP column for those rows is 3.94, 3.53, and 2.86, i.e., a range of 2.86-3.94, not 2-12. The 12.64 figure comes from the univariate Chkrebtii example (N=30), which is not in the emphasized multivariate set. Additionally, the column labeled \"rodeo (no blocking)\" is ambiguous: the values 3.49, 2.77, 4.42 appear to indicate the speedup factor from blocking (consistent with the text that \"blocking is 3-4 times faster\"), but the header suggests it reports the speed of the no-blocking solver itself. Please correct the summary ranges and clarify what the last column measures.","section":"§5.13, Table 3"},{"comment":"The speed comparisons are based on the smallest N for which a Laplace posterior is \"indistinguishable\" from the RKDP posterior. This criterion is not formalized or quantified. Without a concrete measure (e.g., Wasserstein distance between marginal posteriors, KL divergence, or a check that credible interval coverages match), the reported speedups are not reproducible and the \"fast, accurate\" claim is not fully supported. Please specify the comparison metric used to select N.","section":"§5.13, Table 3"}],"minor_comments":[{"comment":"In the blocking definitions, the blocks a^(k)_n, B^(k)_n, and V^(k)_n are described as having dimensions p_k×1, p_k×p_k, and p_k×p_k, respectively. However, these quantities must have the same number of rows as the ODE equations for variable k (r_k); for example, in a first-order univariate ODE with three state derivatives, B_n is 1×3, not 3×3. Please correct the notation to make the dimensions consistent with the observation model W^(k) + B^(k) having shape r_k × p_k.","section":"§2.3, Eq. (14)"},{"comment":"The table caption says \"The proposed method in the accompanying reference(s)...\" but the table does not list rodeo itself; it lists only existing methods. Please clarify whether the table is intended to include rodeo or to survey prior work.","section":"§1, Table 1"},{"comment":"There is a typo: \"An significant contribution\" should be \"A significant contribution.\"","section":"§3.5"},{"comment":"The notation pβ(Θ, Ũ1:N | Z1:N, Y0:M) omits the conditioning value Z1:N = 0; the posterior should be written pβ(Θ, Ũ1:N | Z1:N = 0, Y0:M) for consistency with Equation (38).","section":"§3.5, Eq. (39)"},{"comment":"The statement that unobserved components are handled by setting the observation variance to zero and using \"the log density of Normal(0; 0, 0)\" is not mathematically well-defined, since a normal distribution with zero variance is degenerate. Please clarify how the code actually handles these entries (e.g., by masking or by using a small positive variance).","section":"§5.10"},{"comment":"The formula for the MAGI prior temperature β in the text (β = η^{-2} Δt^{2-2q} Δt') and the code (beta = dt_obs * dt_sim ** (2 - 2 * n_deriv) * sigma[0] ** (-2)) should be cross-checked for notational consistency, particularly the role of q versus n_deriv, so that readers can reproduce the choice.","section":"§5.7"},{"comment":"There are numerous typographical and formatting issues in the code listings and acknowledgments (e.g., \"Cananda\" for \"Canada\", inconsistent spacing in comments, stray characters). A careful proofreading pass would improve the presentation.","section":"Throughout"}],"recommendation":"major_revision","confidential_remarks":"The paper is a software contribution in scope for stat.CO. The core algorithms appear sound and the library has real utility. The main concerns are (1) an unsupported accuracy claim for the block-diagonal Jacobian approximation that underlies the linear-scaling headline, and (2) an inconsistent and ambiguous speed summary in Section 5.13. Both are fixable with additional experiments and text corrections. I would not recommend reject; the revisions are within the manuscript's scope."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is a solid software paper, not a methodological breakthrough. The genuinely new artifact is rodeo, a JAX/Python library that unifies several probabilistic ODE inference methods (Basic, Fenrir, DALTON, marginal MCMC, MAGI) under one filtering framework and adds a modest but useful generalization: MAGI on arbitrary-order ODEs with a Markov prior. It fills a real gap: the main alternative, ProbNumDiffEq.jl, is Julia-first, and existing Python/R tools are slower or narrower. The code looks real, the examples actually run, and the comparisons against LSODA, RKDP, and diffrax are the right kind of sanity check.\n\nWhat the paper does well: the API is clean, the algorithms are traceable to the cited literature, and the empirical section covers three realistic inference settings, including partially observed systems and a Poisson measurement model. The linear-in-d blocking is not invented here; it is Krämer et al. (2021), but the implementation in JAX with square-root Kalman variants is a useful engineering contribution. The examples support the central claim that these methods recover parameters at moderate N.\n\nSoft spots, in proportion. First, the speed claim in Section 5.13 is worded too strongly: “2 to 12 times faster than RKDP” appears right after presenting Table 3, but for the three multivariate examples the RKDP column is 3.94, 3.53, and 2.86. The 12.64 is the univariate Chkrebtii row. That range is misread from the table and should be fixed. Second, the linear-in-d claim rests on the block-diagonal modified Jacobian J*; Section 2.3 says there is “extensive evidence” of minimal accuracy loss but the paper reports no direct comparison against the full Jacobian. The stress-test note is right that this matters for strongly coupled systems. It is not a fatal flaw—the examples are moderately coupled and the speed-accuracy trade-off at least holds there—but a quantitative ablation would make the headline claim much safer. Third, the benchmark comparisons are visual: “indistinguishable from the true posterior” is asserted, not measured. A small table of posterior distances (e.g., Wasserstein or KL) would turn this from impression into evidence. Finally, the repository link is given, but the preprint lacks a commit hash and the data files (hes1data.npy, seirahdata.npy) are referenced but not included; for a library paper, ship those or at least pin the version.\n\nBottom line: this deserves a serious referee. The issues are fixable and mostly about precise reporting, not about a broken method. The likely audience is applied statisticians and ML researchers who want Python-native probabilistic ODE inference, and for them this paper is useful. I would ask for the speed-language correction, a Jacobian ablation, and quantitative posterior comparisons before accepting, but I would not desk-reject it.\n\nRecommendation: send to peer review, conditional on those revisions.","headline":"A useful JAX/Python library paper that deserves peer review, provided the authors tighten the speed claims, add a Jacobian ablation, and quantify 'indistinguishable'.","tokens_in":46123,"tokens_out":2379,"would_cite":false,"duration_ms":28620,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"rodeo, a JAX-based Python library, claims that probabilistic ODE parameter inference can scale linearly in both evaluation points and system variables, matching deterministic-solver accuracy while running 2–12x faster in its tested…","keywords":["probabilistic ODE solvers","parameter inference","Bayesian filtering","Kalman filter","JAX","integrated Brownian motion","block-diagonal computation","dynamical systems"],"falsifier":"Run the same parameter-inference benchmark twice at a fixed step size—once with the Krämer block-diagonal Jacobian interrogation and once with the full Tronarp Jacobian—and measure the distance between the two resulting parameter posteriors (e.g., Wasserstein distance). If the posteriors separate by more than the solver's own uncertainty at moderate coupling strength, the 'minimal loss in accuracy' premise fails and the linear-scaling speed advantage no longer comes free.","tokens_in":45198,"feed_emoji":"🧮","tokens_out":9757,"duration_ms":90085,"temperature":0.7,"pith_summary":"This paper presents rodeo, a Python library that makes probabilistic ODE solvers practical for parameter estimation. Its central claim is that solving an ODE with quantified numerical uncertainty—by treating the solver as a Bayesian filter—can be done in time linear in both the number of grid points and the number of system variables, and that the resulting parameter posteriors match those of a high-accuracy deterministic solver at a fraction of the compute. Correct, this would let practitioners report uncertainty from numerical integration as part of their parameter estimates without paying the cost usually associated with probabilistic methods. The library unifies five inference approaches (Basic, Fenrir, DALTON, marginal MCMC, and a Markov-process version of MAGI) under one interface, supporting partially observed systems and non-Gaussian measurement models.","feed_headline":"Fit ODE parameters 2–12x faster with quantified uncertainty","feed_subtitle":"rodeo's probabilistic solver quantifies numerical error while matching deterministic posteriors at lower cost.","key_machinery":"The central object is the Kalman-filter probabilistic ODE solver with block-diagonal 'model interrogation.' A Gaussian Markov process prior (by default integrated Brownian motion, giving transition matrix $Q$ and covariance $R$) is placed on the solution and its derivatives, and the ODE residual $W\\mathbf{X}_n - f(\\mathbf{X}_n, t_n)$ is treated as a zero-mean observation $Z_n = 0$ in the state-space model, so that posterior draws and means come from Kalman filtering and smoothing. The claim of linear-in-$d$ scaling rests on the block-diagonal structure of the weight matrix $W$ together with the modified Jacobian $J^*_f$, which zeroes out the off-diagonal blocks of the full Jacobian so that each variable's Kalman recursions run independently. That device, plus JAX's automatic differentiation and JIT compilation, is what carries the speed-accuracy results.","core_discovery":"The paper claims that the traditional cost gap between probabilistic and deterministic ODE solvers can be closed by combining a Kalman-filter formulation of the solver with variable-wise block-diagonal structure. The ODE residual is treated as a zero-mean pseudo-observation in a nonlinear state-space model; a first-order Taylor 'model interrogation' linearizes it; and by keeping only the block-diagonal part of the resulting Jacobian (the Krämer-modified $J^*_f$), the Kalman recursions decompose per system variable, reducing complexity from $O(d^3)$ to $O(d)$ in the number of variables $d$. On top of this solver, several likelihood approximations from the literature—the plug-in Basic method, Fenrir, DALTON, the Chkrebtii marginal MCMC, and a Markov-prior version of MAGI—are expressed uniformly and implemented with JAX automatic differentiation and just-in-time compilation. In the FitzHugh-Nagumo, Hes1, and SEIRAH benchmarks, the paper shows these posteriors becoming indistinguishable from a high-accuracy RKDP deterministic solver's posterior at moderate step sizes, while running 2–12 times faster than LSODA and RKDP.","pith_inferences":["A direct accuracy comparison between the block-diagonal Jacobian and the full Jacobian is the missing experiment; a natural follow-up would measure posterior divergence versus coupling strength to map where linear scaling holds.","If the insensitivity to the IBM prior scale noted in the examples holds generally, the solver's reported uncertainty bands are effectively determined by discretization error and the data rather than the user's prior choice—a testable calibration property.","The same blocking device could extend to time-varying parameters, which the paper names as future work, with no change to the $O(d)$ complexity.","If the speed advantage replicates on stiff systems—listed by the authors as an open direction—the practical case for replacing deterministic integrators in routine model calibration would be much stronger."],"forward_implications":["Parameter posteriors from the Basic, Fenrir, and DALTON approximations converge to the deterministic high-accuracy posterior as the solver grid refines, so numerical integration error can be treated as part of the statistical model instead of being hidden.","Because the block-diagonal smoother keeps memory linear in the number of grid points, rodeo can backpropagate directly through the solver steps, making gradient-based inference faster than adjoint-method alternatives.","The same machinery covers partially observed systems (Hes1 with unobserved $H(t)$, SEIRAH with two observed compartments) and non-Gaussian noise (Poisson count data), so the probabilistic-solver benefit extends beyond textbook Gaussian settings.","Blocking itself is worth 3–4× over the unblocked $O(d^3)$ solver in the reported timings, and the full solver is 2–12× faster than LSODA and RKDP at matched posterior accuracy.","As the number of evaluation points $N$ grows, the solver's posterior mean approaches the true ODE solution while its uncertainty bands shrink, meaning the uncertainty output is a usable diagnostic for under-resolved discretizations."],"supporting_citations":[{"why":"Supplies the block-diagonal modified Jacobian $J^*_f$ and the blocking scheme that reduces the multivariate solver from $O(d^3)$ to $O(d)$.","marker":"Krämer et al. (2021)"},{"why":"Establishes the Bayesian nonlinear filtering formulation of probabilistic ODE solvers, including the state-space model and first-order Taylor model interrogation.","marker":"Tronarp et al. (2018)"},{"why":"Provides the integrated Brownian motion prior with its $Q, R$ matrices, plus the simplest zeroth-order interrogation used as a baseline.","marker":"Schober et al. (2019)"},{"why":"Supplies the marginal MCMC algorithm for parameter inference and the motivating second-order example used to benchmark solver accuracy.","marker":"Chkrebtii et al. (2016)"},{"why":"Contributes the MAGI posterior derivation and the Hes1 benchmark model; the paper extends their method to a Markov prior.","marker":"Yang et al. (2021)"},{"why":"Defines the data-adaptive DALTON likelihood approximation, including the non-Gaussian variant used for the SEIRAH model.","marker":"Wu and Lysy (2024)"},{"why":"Defines the Fenrir likelihood approximation that yields a closed-form Gaussian likelihood within the filtering framework.","marker":"Tronarp et al. (2022)"},{"why":"Provides JAX, the implementation backbone whose automatic differentiation and JIT compilation deliver the reported speedups.","marker":"Bradbury et al. (2018)"}],"fun_headline_variants":["rodeo: ODE parameter inference up to 12x faster with uncertainty","Linear-scaling probabilistic solver speeds ODE fitting 2-12x","Fast probabilistic parameter inference for ODEs with rodeo","rodeo library brings probabilistic ODE solvers to production speed","Quantified uncertainty for ODE params at 2-12x lower cost"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that zeroing out the off-diagonal blocks of the ODE Jacobian (the Krämer modification) causes minimal accuracy loss; the paper asserts 'extensive evidence' for this but reports no quantitative comparison, so for strongly coupled systems the linear-time solver's posterior could drift from the full-Jacobian answer.","fun_headline_variants_meta":{"raw":{"variants":["rodeo: ODE parameter inference up to 12x faster with uncertainty","Linear-scaling probabilistic solver speeds ODE fitting 2-12x","Fast probabilistic parameter inference for ODEs with rodeo","rodeo library brings probabilistic ODE solvers to production speed","Quantified uncertainty for ODE params at 2-12x lower cost"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000504,"raw_usage":{"total_tokens":2455,"prompt_tokens":935,"completion_tokens":1520,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":551,"completion_tokens_details":{"reasoning_tokens":1427}},"tokens_in":551,"tokens_out":1520,"duration_ms":11835,"temperature":1.0,"reasoning_tokens":1427,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T22:19:20.151575+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the same parameter-inference benchmark twice at a fixed step size—once with the Krämer block-diagonal Jacobian interrogation and once with the full Tronarp Jacobian—and measure the distance between the two resulting parameter posteriors (e.g., Wasserstein distance). If the posteriors separate by more than the solver's own uncertainty at moderate coupling strength, the 'minimal loss in accuracy' premise fails and the linear-scaling speed advantage no longer comes free.","supporting_citations":[{"cited_title":"Probabilistic Solutions To Ordinary Differential Equations As Non-Linear Bayesian Filtering: A New Perspective","cited_arxiv_id":"1810.03440","evidence_quote":"Establishes the Bayesian nonlinear filtering formulation of probabilistic ODE solvers, including the state-space model and first-order Taylor model interrogation."},{"cited_title":"Bayesian solution uncertainty quantification for differential equations","cited_arxiv_id":null,"evidence_quote":"Supplies the marginal MCMC algorithm for parameter inference and the motivating second-order example used to benchmark solver accuracy."},{"cited_title":"Inference of dynamic systems from noisy and sparse data via manifold-constrained Gaussian processes","cited_arxiv_id":null,"evidence_quote":"Contributes the MAGI posterior derivation and the Hes1 benchmark model; the paper extends their method to a Markov prior."},{"cited_title":"Data-Adaptive Probabilistic Likelihood Approximation for Ordinary Differential Equations","cited_arxiv_id":null,"evidence_quote":"Defines the data-adaptive DALTON likelihood approximation, including the non-Gaussian variant used for the SEIRAH model."},{"cited_title":"JAX : composable transformations of P ython+ N um P y programs","cited_arxiv_id":null,"evidence_quote":"Provides JAX, the implementation backbone whose automatic differentiation and JIT compilation deliver the reported speedups."}],"review_version":1}