{"id":"76747166-25e1-4674-92cd-53219b200556","arxiv_id":"2501.15767","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"Verification of Markov chains with ML-computed probabilities and rewards becomes a global optimization problem that a decomposition and bound propagation scheme solves, often hundreds of times faster than direct solving.","lead":"The paper introduces a way to prove worst-case properties of Markov chains whose transition probabilities and rewards are computed by machine learning models. It encodes these checks as optimization problems and accelerates global solving with bound propagation, backed by an open-source tool and a healthcare cost example.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Formal global-optimality claim rests on MILP-representability, but the experiments and supported classifiers use logistic/softmax models handled only by numerical spatial branch-and-bound, so the formal guarantee is not established for those models.","rationale":"The paper's contribution is genuinely valuable for exactly MILP-representable learned parameters: the bilinear formulation is valid, the bound-propagation lemmas are standard and correct, and adding valid bounds on π, P, r, and v cannot cut off an optimal feasible point. The decomposition and interval Gauss-Seidel tightening are reasonable and the released code is a plus. The reader's weakest assumption correctly identifies the soft spot: the formal guarantee is only proven under MILP representability, yet the experiments that support the abstract's 'global optimality' and '100x faster' claims use logistic regression and softmax classifiers, which are not exactly MILP-representable. The paper's own text in Section 3.1 concedes this by appealing to Gurobi's spatial branch-and-bound, and Appendix F.2's claim of an exact MILP formulation for softmax is not mathematically possible for a finite MILP. Since the bounds on θ are the root of all downstream propagation, an invalid or merely approximate bound there undermines the formal status of the final global optimum. This is not a rejection of the exact-MILP core; it is a demand to either certify the nonlinear handling to a rigorous tolerance or restrict the formal claims. The reader's CONDITIONAL verdict already captures this, so no verdict change is needed. Secondary concerns, such as the runtime reclassification rule in Appendix G, are worth reporting but are less load-bearing than the exactness gap.","tokens_in":22839,"tokens_out":7046,"duration_ms":70264,"concrete_test":"Run the smallest state-space experiment (n=5) with logistic-regression probabilities and compare the θ-bounds produced by markovml/Gurobi's spatial branch-and-bound against rigorous enclosures computed by an independent interval-arithmetic tool (e.g., INTLAB or Arb) for the same feature set X=[-1,1]^5. If any reported θmin or θmax lies outside the rigorous enclosure by more than the solver tolerance, the bound-propagation chain is formally invalid for logistic models. Repeat once for a 5-neuron softmax MLP classifier. If the rigorous enclosures confirm the bounds, the concern is resolved for the tested cases; if not, the formal claims must be restricted to exactly MILP-representable models.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim of exact global verification depends on the Section 3.1 assumption that every f_i is MILP-representable. That assumption is not satisfied by the models used in the headline experiments: the state-space and model-count experiments use logistic regression, and supported neural classifiers include softmax layers (Appendices F.1, F.2, G). The paper itself states in Section 3.1 that logistic/softmax are handled 'using Gurobi’s built-in spatial branch-and-bound techniques for nonlinear functions – essentially dynamic piecewise linear approximations.' A finite MILP cannot exactly encode the transcendental softmax or logistic function, and a dynamic piecewise-linear approximation with default numerical tolerances is not an exact MILP encoding. Consequently, the bounds on θ computed in Section 4 are not formally guaranteed lower/upper bounds for these models, and all propagated bounds on π, P, r, and v inherit that uncertainty. If the numerical global solver reports an optimum only to tolerance, the 'formal verification' guarantee is not established for exactly the models used in the experimental evaluation and case study. Appendix F.2's statement that the authors 'implemented the MILP formulation ourselves' for softmax masks rather than resolves this gap. The exact-MILP core (trees, ReLU networks without softmax, linear models) is sound; the load-bearing weakness is the external validity of the formal claim for the broader model class.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper introduces a framework for formally verifying properties of finite-state Markov chains whose transition probabilities, initial distribution, and rewards are determined by machine learning models. The verification task (reachability, hitting time, total discounted reward) is cast as a mixed-integer bilinear program under an assumption that each ML model is MILP-representable. The paper's main algorithmic contribution is a decomposition and bound-propagation scheme: compute elementwise bounds on the learned outputs by solving small MILPs, propagate these to π, P, and r via affine equalities, derive coarse bounds on the value function v using Markov-chain theory, tighten those bounds with interval Gauss-Seidel, and then solve the original bilinear program with an off-the-shelf global solver, using the propagated bounds as constraints. Experiments on synthetic instances compare this approach with direct solving, reporting large speedups and better scaling with state count, model count, tree depth, and neural-network width, and a case study applies the method to a cardiac-arrest cost-effectiveness model with a tree-based discharge classifier. The paper releases an open-source Python package, markovml, and scripts to reproduce the experiments.","tokens_in":22936,"tokens_out":8645,"duration_ms":78116,"significance":"If scoped appropriately, this is a useful and timely contribution: the idea of embedding learned parameters directly into an optimization formulation for Markov-chain verification is natural but nontrivial, and the decomposition into small MILPs plus interval linear-system tightening is a practically effective way to accelerate global optimization. The theoretical core for exactly MILP-representable models (trees, ReLU networks without softmax, and linear models) is sound: Lemma 4.1 propagates interval bounds correctly, Lemma 4.3 gives valid value-function bounds, and Theorem 4.4 provides a clean condition for the interval-Gauss-Seidel hull property. The reported speedups for decision-tree instances (tens to hundreds of times, with a 1000x speedup on some instances) are substantial, and the ablation in Table 1 isolates the contribution of the v-bounds. The release of documented code (markovml) and reproducible experimental scripts is a genuine strength. However, the paper's headline claims of 'formal verification' and 'global optimality' over 'a broad class of machine learning models, including ...","major_comments":[{"comment":"The core theoretical framework is conditional on the assumption that every f_i is MILP-representable, but the experimental evaluation uses logistic regression (Experiments 1 and 2) and softmax-based classifiers (Experiment 4 and the torch support in Section 5), which do not satisfy this assumption. Section 3.1 itself acknowledges that logistic and softmax functions are handled 'using Gurobi’s built-in spatial branch-and-bound techniques for nonlinear functions – essentially dynamic piecewise linear approximations.' A finite MILP cannot exactly encode a transcendental function, and the bounds on θ computed in Section 4 are therefore not certified lower/upper bounds for these models. Consequently, the propagated bounds on π, P, r, and v are not formally valid for these instances, and the abstract's claim that the method 'solves the problem to global optimality' is not established for exactly the model classes used in the headline experiments. The fix is to restrict the formal claims (abstract, introduction, Section 4) to exactly MILP-representable models, and to either re-run the affected experiments with such models, remove them from the formal-claim experiments, or provide a rigorous analysis of the spatial-branch-and-bound tolerances and their effect on the final optimum.","section":"Section 3.1, Section 4, Appendix G"},{"comment":"The statement that the authors 'implemented the MILP formulation ourselves' for softmax is misleading: a piecewise-linear approximation of softmax encoded with binary variables is not an exact MILP representation. Unless the paper provides a certified error bound for this approximation (e.g., a formal statement of how the piecewise approximation relates to the true softmax, with a guarantee that the final objective value is bracketed within a stated tolerance), the results for softmax-based networks should be presented as numerical or approximate, not as exact formal verification. This is not a mere presentational point; it determines whether the advertised 'formal guarantees' hold for a large portion of the supported model zoo.","section":"Appendix F.2, Section 5"}],"minor_comments":[{"comment":"The sentence listing 'linear and logistic regression' among MILP-representable functions is inaccurate, and the paper's own caveat immediately below contradicts it. Rephrase to avoid implying that logistic regression is exactly encoded by a finite MILP.","section":"Section 3.1"},{"comment":"There is a duplicated phrase: 'from a set T ⊆ [n] of transient states from a set T ⊆ [n] of transient states equal to the complement of S.' Remove the repetition.","section":"Definition 3.2"},{"comment":"The runtime comparison reclassifies suboptimal solver statuses as 'optimal' when the objectives agree to 1e-12. This is disclosed, but it conflates cases where the direct solver proved optimality with cases where it only found a feasible solution with the same objective. Please report the runtime comparison separately for instances where Gurobi proved optimality, and present the reclassified analysis only as a sensitivity check.","section":"Section 6.1 and Appendix G"},{"comment":"The error bars are described as 'standard deviation' of geometric means, but the multiplier-based display is ambiguous. State explicitly whether the standard deviation is computed on log-transformed runtimes and whether the bars are multiplicative (e.g., geometric mean times division by the geometric standard deviation).","section":"Figure 2"},{"comment":"The NeurIPS checklist item on Limitations is left as '[TODO]' in the submitted manuscript. Please complete it and include a discussion of the MILP-representability scope in the Limitations section of the main text.","section":"Checklist, Item 2"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is competent and the core algorithmic idea is sound for MILP-representable models, with genuinely impressive speedups on decision-tree instances. The difficulty is scope control: the abstract and introduction claim formal guarantees over a broad class including neural networks, while the experiments rely heavily on logistic/softmax models that do not fit the MILP-representability assumption. This is a fixable issue — the authors should either narrow the formal claims or provide a rigorous error analysis for the nonlinear case — rather than a fatal technical error. I would not reject, and I see no circularity or fabrication in the results. I would ask the authors to also tighten the runtime-comparison methodology and to complete the checklist. Given the paper's current overclaim, major revision is the appropriate recommendation."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Quick take: the paper defines a genuinely new verification problem — Markov chains whose transition probabilities, rewards, and initial distribution are outputs of ML models — and reduces it to a mixed-integer bilinear program with a decomposition and bound propagation scheme. For models that are exactly MILP-representable (trees, ReLU nets without softmax, linear regressions), the core math is sound: the bilinear encoding is valid, the bound propagation lemmas check out, and the interval Gauss-Seidel tightening is standard and correctly applied. The interval M-matrix condition (Theorem 4.4) is a nice touch, and the released code and software are a real plus. The 100x speedup claim is credible for the exact-MILP cases tested.\n\nThe soft spot is the gap between the formal claim and the experiments. Under the paper's own Section 3.1 assumption, logistic regression and softmax networks are not MILP-representable; they are handled with Gurobi's spatial branch-and-bound — essentially dynamic piecewise-linear approximation with numerical tolerances. That means the 'formal verification' guarantee does not actually extend to the logistic models used in the state-space and model-count experiments, or to the softmax classifiers. This is not a fatal flaw in the core method, but the abstract and conclusion overstate the scope. The fix is straightforward: either restrict the claims to exactly MILP-representable models, or present the nonlinear-case results as sound numerical bounds with explicit tolerances, not formal guarantees.\n\nOne more caution: the runtime analysis in Appendix G only keeps instances solved to optimality and then reclassifies the direct solver's 'suboptimal' results as optimal when the objective values match to 1e-12. That rule can only help their method, so the speedup ratios should be read with that caveat. A dedicated baseline against a specialized bilinear solver would strengthen the comparison.\n\nReproducibility is adequate: code is released, but no commit hash, incomplete seeds, and no patient data for the case study — minor.\n\nWho should read this: anyone working on probabilistic verification with learned parameters, healthcare cost-effectiveness analysts who want exact subgroup analyses instead of Monte Carlo, and people building solvers for bilinear programs with large variable bounds. It deserves a serious referee; with a scope-clarifying revision and tighter runtime reporting, it's a solid contribution.","headline":"Useful new problem class and a sound bound-propagation scheme for the exact-MILP core, but the 'formal' guarantee is overclaimed for the logistic/softmax models actually used in the experiments.","tokens_in":23609,"tokens_out":3566,"would_cite":true,"duration_ms":29127,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["90C26","90C11","60J10"],"pacs":[],"model":"deepseek-v4-flash","headline":"Formal verification of Markov processes with ML-learned parameters reduces to a mixed-integer bilinear program, solvable to global optimality by a decomposition and bound-propagation scheme that is orders of magnitude faster than direct…","keywords":["formal verification","Markov chains","machine learning verification","bilinear programming","bound propagation","interval Gauss-Seidel","global optimization","markovml"],"falsifier":"Run markovml on a two-state Markov reward process whose transition probability is a logistic regression in one feature $x \\in [0,1]$ with a fixed reward vector, and compare its claimed global minimum of total discounted reward against an independent dense scan of $x$ using the closed-form value $\\pi^\\top(I-\\lambda P(x))^{-1}r$; if the scan finds a strictly better objective than the certified optimum by more than the solver tolerance, the global-optimality claim for nonlinear models is false.","tokens_in":22459,"feed_emoji":"🧮","tokens_out":7010,"duration_ms":65034,"temperature":0.7,"pith_summary":"Markov chains used in high-stakes settings increasingly take their transition probabilities and rewards from machine-learning models, but proving worst-case properties of the resulting system has been out of reach. The paper claims that reachability, expected hitting time, and total discounted reward for such chains can be formulated as one mixed-integer bilinear program, and that a decomposition scheme—bound each ML output with small MILPs, propagate bounds through the affine parameter maps, then tighten value-function bounds with interval Gauss-Seidel—makes the program tractable to global optimality. In numerical experiments the method solves instances orders of magnitude faster than direct bilinear solving, with the gap growing with model complexity, and it is packaged in an open-source tool called markovml. If the claim holds, it turns formal verification of ML-parameterized Markov processes from a Monte Carlo approximation into an exact optimization problem.","feed_headline":"Exact verification for ML-driven Markov chains","feed_subtitle":"A decomposition scheme solves the bilinear program globally, beating direct solvers by up to 1000x.","key_machinery":"The load-bearing mechanism is a four-stage bound-propagation decomposition of the bilinear program. Stage one solves $2\\ell$ small MILPs to bound each component of the ML output $\\theta$ over the feature set $X$; stage two pushes these intervals through the affine maps $\\pi = A_\\pi\\theta + b_\\pi$, $\\mathrm{vec}(P) = A_P\\theta + b_P$, $r = A_r\\theta + b_r$ using Lemma 4.1; stage three exploits the fact that every row of $(I-\\lambda P)^{-1}$ sums to $1/(1-\\lambda)$ to get initial bounds on the value vector $v$; stage four runs interval Gauss-Seidel on the linear system $(I-\\lambda P)v = r$, which is guaranteed by Theorem 4.4 to return the tightest rectangular enclosure (the hull) exactly when the upper-bound transition matrix satisfies $\\rho(P_{\\max}) \\le 1/\\lambda$. These bounds are then added as constraints to the original bilinear program, and any out-of-the-box global solver finishes the job.","core_discovery":"The central claim is that formal verification of Markov processes whose parameters are outputs of ML models is a solvable optimization problem: for linear models, tree ensembles, and ReLU networks, computing the worst- or best-case reachability probability, hitting time, or total reward can be written as a mixed-integer bilinear program. The paper's method solves that program by first bounding each ML output through small MILPs, propagating the resulting intervals through the affine equalities that define $\\pi$, $P$, and $r$, deriving initial value-function bounds from the Markov structure via $(I-\\lambda P)^{-1}$, and tightening those bounds with interval Gauss-Seidel. The final bilinear program with the propagated bounds is passed to a global solver, so the result is a certificate of the extremal property, not a simulation estimate. The paper reports that this pipeline solves to global optimality with speedups ranging from tens to roughly a thousand times over direct solving, and demonstrates the pipeline on a cardiac-arrest cost-effectiveness model.","pith_inferences":["The $\\theta$-bounding stage is independent of the Markov chain, so plugging in stronger neural-network verifiers—for example, better branch-and-bound for ReLU networks—should tighten bounds and extend scalability to deeper architectures, a direction the paper itself flags.","The spectral-radius condition in Theorem 4.4 doubles as a predictability diagnostic: instances with $\\rho(P_{\\max}) > 1/\\lambda$ are exactly the ones where interval Gauss-Seidel may not reach the hull, so users could estimate in advance which verification problems will be hard.","For logistic regression and softmax classifiers, the exact MILP-representability assumption is not met; the formal guarantee is then inherited from the numerical global solver, so the strongest formal statement should be scoped to piecewise-linear models, with nonlinear models carrying a numerical-global certificate.","The framework suggests a natural extension to subgroup certification: since bounds on $\\theta$ and $v$ are cheap, one could search over feature regions to certify performance across a whole subgroup without enumerating patients."],"forward_implications":["Worst-case and best-case answers to reachability, hitting time, and total discounted reward become provable for Markov chains whose parameters come from linear models, tree ensembles, and small ReLU networks, rather than approximated by Monte Carlo.","The same pipeline handles optimization and feasibility versions, so threshold questions—such as whether the worst-case failure probability stays below a bound—become exact yes-or-no checks.","The speedup over direct bilinear solving grows with model complexity: about 100x for deep trees, sometimes over 1000x, and neural-network instances that time out for direct solving are solved by the decomposition.","Because the final optimization still runs on an out-of-the-box global solver, future improvements in bilinear or MILP solvers should translate directly into faster verification."],"supporting_citations":[{"why":"Defines bilinear programs as NP-hard non-convex quadratic problems solved by branch-and-bound, the problem class the paper targets.","marker":"[46]"},{"why":"Provides McCormick's envelope, the convex relaxation that global solvers use for bilinear terms.","marker":"[47]"},{"why":"The out-of-the-box global solver used for all optimization and for the spatial branch-and-bound handling of nonlinear model outputs.","marker":"[48]"},{"why":"Supplies the interval matrix analysis, hull definition, and interval Gauss-Seidel theory used to tighten value bounds.","marker":"[54]"},{"why":"Provides the M-matrix results used in Lemma 4.2 and in the spectral-radius condition of Theorem 4.4.","marker":"[55]"},{"why":"Provides the MILP formulations for several sklearn models used in the markovml implementation.","marker":"[58]"},{"why":"Gives the interval linear systems exposition that the paper follows for its Gauss-Seidel implementation.","marker":"[61]"},{"why":"Demonstrates SMT/MILP encoding of ReLU networks, motivating how ML models can be embedded into the optimization framework.","marker":"[32]"}],"fun_headline_variants":["Exact ML-Markov verification solves globally, up to 1000x faster","Formal proof for Markov chains with ML parameters, 1000x speedup","Bilinear decomposition verifies ML-driven Markov chains exactly","Global-optimal reachability for ML-parameterized Markov processes"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"Everything stands on the assumption that each ML model's input-output relation can be exactly encoded with linear constraints and binary variables; for logistic and softmax models, which cannot, the paper relies on the numerical global solver's spatial branch-and-bound, so the formal guarantee for those models is only as strong as that numerical solve.","fun_headline_variants_meta":{"raw":{"variants":["Exact ML-Markov verification solves globally, up to 1000x faster","Formal proof for Markov chains with ML parameters, 1000x speedup","Bilinear decomposition verifies ML-driven Markov chains exactly","Global-optimal reachability for ML-parameterized Markov processes"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000295,"raw_usage":{"total_tokens":1685,"prompt_tokens":887,"completion_tokens":798,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":503,"completion_tokens_details":{"reasoning_tokens":720}},"tokens_in":503,"tokens_out":798,"duration_ms":7949,"temperature":1.0,"reasoning_tokens":720,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T13:57:57.780368+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run markovml on a two-state Markov reward process whose transition probability is a logistic regression in one feature $x \\in [0,1]$ with a fixed reward vector, and compare its claimed global minimum of total discounted reward against an independent dense scan of $x$ using the closed-form value $\\pi^\\top(I-\\lambda P(x))^{-1}r$; if the scan finds a strictly better objective than the certified optimum by more than the solver tolerance, the global-optimality claim for nonlinear models is false.","supporting_citations":[{"cited_title":"Handbook of global optimization, volume 2","cited_arxiv_id":null,"evidence_quote":"Defines bilinear programs as NP-hard non-convex quadratic problems solved by branch-and-bound, the problem class the paper targets."},{"cited_title":"Computability of global solutions to factorable nonconvex programs: Part i—convex underestimating problems","cited_arxiv_id":null,"evidence_quote":"Provides McCormick's envelope, the convex relaxation that global solvers use for bilinear terms."},{"cited_title":"New techniques for the analysis of linear interval equations","cited_arxiv_id":null,"evidence_quote":"Supplies the interval matrix analysis, hull definition, and interval Gauss-Seidel theory used to tighten value bounds."},{"cited_title":"Gurobi machine learning, 2024","cited_arxiv_id":null,"evidence_quote":"Provides the MILP formulations for several sklearn models used in the markovml implementation."},{"cited_title":"Interval linear and nonlinear systems","cited_arxiv_id":null,"evidence_quote":"Gives the interval linear systems exposition that the paper follows for its Gauss-Seidel implementation."},{"cited_title":"Reluplex: An efficient smt solver for verifying deep neural networks","cited_arxiv_id":null,"evidence_quote":"Demonstrates SMT/MILP encoding of ReLU networks, motivating how ML models can be embedded into the optimization framework."}],"review_version":1}