{"id":"56b035c8-fa26-4e00-bea4-5efbd3e7b2f5","arxiv_id":"2507.03159","paper_version":2,"verdict":"ACCEPT","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"low","formal_verification":"none","parameter_count":0,"one_line_summary":"A Julia package that embeds trained ML predictors into JuMP optimization models, uniquely offering GPU-accelerated gray-box nonlinear formulations.","lead":"MathOptAI.jl is a new open-source Julia library that lets you embed trained machine learning models, such as neural networks and decision trees, directly into mathematical optimization problems built with JuMP. Its gray-box formulation is the first in a general-purpose library to offload function and derivative evaluations of large PyTorch models to a GPU during optimization.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Gray-box Hessian derivation in §2.3.1 covers only the full-space constraint path; the advertised reduced-space gray-box path needs a different weight vector and is not specified or tested.","rationale":"The reader's ACCEPT is reasonable if the code archive is trusted, but the correctness of the gray-box derivative plumbing is the exact place where the central claim could break. The paper's written derivation in §2.3.1 only handles the full-space constraint multiplier case, while §2.3 advertises a reduced-space gray-box path. In reduced-space, the second-order weight is the adjoint of the outer expression, not Ipopt's multiplier for y=F(x), and the paper does not specify how that weight is passed through JuMP's operator interface. This is not an exotic corner case: substituting F(x) into a nonlinear objective or constraint is precisely where reduced-space formulations are most useful. The proposed targeted test—comparing gray-box solutions to an algebraic formulation and checking the reported Hessian against finite differences—would settle whether the gap is real. If it passes, the concern is cleared and the paper can be accepted as written; if it fails, the gray-box claim must be narrowed or corrected. Hence CONDITIONAL rather than a flat rejection.","tokens_in":10041,"tokens_out":14714,"duration_ms":182905,"concrete_test":"With F(x)=tanh(Wx+b) (one input, one output, smooth), solve three Ipopt models from the same start: (i) algebraic full-space min (y-1)^2+x^2 s.t. y=F(x); (ii) MathOptAI gray-box full-space; (iii) MathOptAI gray-box reduced-space min (F(x)-1)^2+x^2. If (iii) is not supported, the paper should say so; if it is supported, it must match (i). Additionally, at a non-stationary (x,λ), compare the gray-box operator's reported ∇²(λᵀF) with second-order finite differences of λᵀF. This test isolates whether the Jacobian/Hessian ordering and weight semantics are correct.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The gray-box contribution rests on JuMP/Ipopt receiving correct F(x), ∇F(x), and ∇²(λᵀF)(x). Section 2.3.1 derives ∇²G for the full-space constraint y=F(x), where λ is the vector of multipliers of the P equality constraints. Section 2.3 also states that the gray-box predictor can be added via the reduced-space form y:=F(x). In that form there is no constraint F(x)-y=0; the correct weight multiplying ∇²F_k in the Hessian of the Lagrangian is the total derivative of the surrounding objective and constraint functions with respect to y_k, not Ipopt's multiplier for y=F(x). The paper does not explain how JuMP's user-defined operator interface supplies this weight to the G(x,λ) construction. If MathOptAI reuses the Ipopt multiplier vector in the reduced-space path, the supplied Hessian is wrong and Ipopt can converge to a point that is not stationary for the intended model. Even in the full-space path, correctness depends on unstated sign and ordering conventions for the Jacobian of F(x)-y=0. No numerical validation of the gray-box path appears in the paper; the code archive is the only evidence, so this derivative plumbing is the least secured link in the central claim.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"MathOptAI.jl is an open-source Julia package for embedding trained machine-learning predictors into JuMP optimization models. The paper describes three embedding formulations, full-space, reduced-space, and gray-box, and presents the gray-box approach as the novel contribution: the predictor is exposed to Ipopt through user-defined nonlinear operator callbacks for function, Jacobian, and Hessian-of-the-Lagrangian evaluations, with PyTorch-based oracle evaluations offloaded to a GPU through Julia's Python interface. The paper also gives a usage example, a feature comparison with OMLT, gurobi-machinelearning, PySCIPOpt-ML, and GAMSPy, and a set of design principles. The authors explicitly state that the paper contains no computational experiments and refer readers to a companion paper for performance evidence.","tokens_in":10252,"tokens_out":9949,"duration_ms":128431,"significance":"If the gray-box implementation is correct, this is a genuinely useful software contribution: it is the only package in the Table 1 comparison with a general-purpose gray-box formulation and GPU acceleration, and it fills a gap between Python machine-learning tooling and the JuMP optimization ecosystem. The algebra in Section 2.3.1 is standard and the library is transparently disclosed, with source code archived and prior work by Casas and coauthors credited. I see no circularity issue, since the formulations rest on standard nonlinear programming constructions rather than fitted parameter claims. However, the significance is conditional on the correctness of the derivative plumbing described in Sections 2.3 and 2.3.1, which is the least secured part of the paper: the reduced-space gray-box Hessian derivation is incomplete, and no in-paper validation is provided. The stress-test concern therefore lands and should be addressed before the central gray-box claim is accepted.","major_comments":[{"comment":"The gray-box Hessian derivation is written for the full-space constraint F(x)-y=0, where λ is the vector of Lagrange multipliers of the P equality constraints. The same section, however, states that the predictor can be added via the reduced-space form y := F(x). In the reduced-space problem there are no constraints y=F(x), so the correct weight multiplying ∇²F_k in the Hessian of the Lagrangian is the gradient of the original Lagrangian with respect to y_k, evaluated at y=F(x), not an Ipopt multiplier of a nonexistent constraint. The manuscript never explains how MathOptAI obtains this weight for the reduced-space path, and if it reuses the full-space multiplier vector, the Hessian is wrong and Ipopt may converge to a non-stationary point. Please either give the explicit reduced-space derivation and the exact operator-callback semantics, or restrict the gray-box claim to the full-space formulation.","section":"Sections 2.3 and 2.3.1"},{"comment":"Even for the full-space path, the identity G(x,λ) = (Affine{λ,0}∘F)(x) is correct only under a matching sign and ordering convention. The paper writes the constraint as F(x)-y=0 but does not state whether λ is the multiplier of F(x)-y=0 or y-F(x)=0, and it does not specify the row/column ordering of ∇F and of the Hessian block. Because JuMP and Ipopt also have their own conventions for how user-defined operator Hessians are weighted, an unexplained mismatch can silently produce a wrong search direction. Please state the conventions explicitly and verify them on a small example with a finite-difference or KKT-residual check.","section":"Section 2.3.1"},{"comment":"The paper explicitly disclaims computational experiments and directs readers to a companion paper for performance evidence. This disclosure is honest, but for the gray-box claim it leaves the central novel mechanism unverified within the manuscript itself: the correctness of the derivative callbacks, not just their speed, is exactly what the algebraic description cannot establish with certainty. If the archived code contains derivative-consistency tests or correctness tests, the manuscript should say so and summarize them; otherwise, please add a minimal numerical example showing that Ipopt converges to the correct KKT point in both full-space and reduced-space gray-box modes.","section":"Sections 1.2, 2.4, and 6"}],"minor_comments":[{"comment":"The display for BinaryDecisionTree mixes an output expression with a set of side constraints; please define the path set P and leaf values y_p precisely, and avoid overloading y for both the predictor output and the leaf values.","section":"Section 1.1"},{"comment":"The default non-smooth nonlinear ReLU formulation is used with Ipopt in the main code example, but Ipopt is a smooth NLP solver; the paper should either note which nonsmooth-capable solvers are expected to work with this default or explain how the formulation is regularized when Ipopt is selected.","section":"Section 3"},{"comment":"For consistency and readability, use explicit 'No' entries instead of blanks in the Formulations and GPU-acceleration rows, so that readers can distinguish 'not supported' from 'not listed'.","section":"Table 1"}],"recommendation":"major_revision","confidential_remarks":"This is a worthwhile software-description paper, and I would not reject it. My recommendation of major revision rests on one load-bearing point: the reduced-space gray-box Hessian derivation is incomplete and potentially wrong as written. The authors should show the exact JuMP/MOI user-defined operator semantics that supply the multiplier weights, or restrict the gray-box claim to full-space. A small KKT/finite-difference test for both gray-box modes would also substantially increase confidence. If the implementation is already correct, this should be straightforward to document."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"This is a genuinely useful software paper, and the authors are honest about what it is: a general-purpose library for embedding trained ML predictors into JuMP, with the gray-box idea attributed to Casas et al. rather than claimed as new. The design principles are sensible, the comparison table is informative, and the code is archived. If you work on optimization with learned surrogates, this is a practical tool worth knowing about.\n\nWhat the paper does well: it packages the gray-box formulation in a way that no other compared library does, supports GPU evaluation of the required oracles through PyTorch, and spells out the full-space, reduced-space, and gray-box formulations clearly. The composition abstraction (everything is a predictor F(x) mapping vectors to vectors) is a good design choice and is explained with concrete code. The authors also explicitly point to Parker et al. for numerical scaling results, so they are not overclaiming inside this paper.\n\nThe soft spots are real but not fatal. First, the stress-test note about reduced-space gray-box Hessians has force. Section 2.3.1 derives the Hessian contribution as the derivative of G(x, lambda) = lambda^T F(x), where lambda is the multiplier on the constraint F(x) - y = 0. That is correct for the full-space path. But the reduced-space path y := F(x) has no such equality constraint, so the correct weight on the Hessian of F is not an Ipopt multiplier for F(x)-y=0; it is the total derivative of the Lagrangian with respect to the operator's output. The paper does not say how JuMP's user-defined operator interface supplies that weight in the reduced-space case. The implementation might well be correct—the archive is the only evidence, and no numerical check appears in the paper—but the missing explanation leaves a hole in the central claim. A small verification example, even in an appendix, would close it.\n\nSecond, the paper has no standalone numerical validation. The authors disclose this and point to a companion paper, which is fair, but for a software contribution a minimal sanity check (e.g., a small problem where the gray-box solution matches the full-space solution) would strengthen the submission at low cost.\n\nOverall: this is a solid engineering artifact, clearly written, with honest attribution and a useful feature comparison. The math that is shown is standard; the gap is in documenting the derivative plumbing for one advertised path. A serious referee should be asked to examine that plumbing, but this should not be desk-rejected—it deserves review.\n\nRecommendation: send to peer review, and ask the authors to clarify or test the reduced-space gray-box Hessian path.","headline":"A useful, honestly-scoped software paper that deserves peer review, but the gray-box Hessian story is under-specified for the reduced-space path—worth asking the authors to pin down.","tokens_in":10806,"tokens_out":2039,"would_cite":true,"duration_ms":27298,"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":"This paper presents MathOptAI.jl, a Julia library that embeds trained machine-learning predictors into JuMP optimization models, including a gray-box mode that offloads derivative computation to a GPU.","keywords":["MathOptAI.jl","JuMP","Julia","machine learning predictors","gray-box formulation","GPU acceleration","PyTorch","nonlinear optimization"],"falsifier":"Take a small pretrained network with known weights, embed it in JuMP once in full-space and once in gray-box, and solve identical optimization problems with Ipopt; a disagreement between the two solutions beyond solver tolerance, or a mismatch between the gray-box Jacobian and Hessian-of-the-Lagrangian and finite-difference values at a random feasible point, would show that the oracle bridge does not represent y = F(x).","tokens_in":9813,"feed_emoji":"🧠","tokens_out":10419,"duration_ms":103769,"temperature":0.7,"pith_summary":"MathOptAI.jl is a Julia library whose purpose is to take a trained machine-learning predictor—a neural network, decision tree, random forest, or Gaussian process—and make it part of a mathematical optimization model built in JuMP. The paper's central claim is that this embedding can be done in three ways: a full-space formulation that adds variables and constraints for every piece of the predictor, a reduced-space formulation that inlines the predictor as one nonlinear expression, and a gray-box formulation that leaves the predictor's internals outside the solver and supplies function, Jacobian, and Hessian-of-the-Lagrangian evaluations through callbacks. The gray-box path is the paper's new contribution: among the five packages compared in Table 1, MathOptAI.jl is the only one with this formulation, and it is what allows PyTorch models to be evaluated on a GPU while the rest of the model remains a Julia-side JuMP problem. A sympathetic reader should care because the library removes the need for a closed algebraic form: if a predictor has trained weights and a differentiation engine, it can be optimized over directly.","feed_headline":"Library embeds trained ML models into JuMP optimization","feed_subtitle":"Gray-box mode keeps the solver model small and moves network derivative math to the GPU.","key_machinery":"The load-bearing object is the AbstractPredictor interface, which forces every supported model to act as a vector-in, vector-out function y = F(x), and the Pipeline wrapper that composes predictors so a neural network is just a sequence of layers. On top of this, the gray-box formulation rests on JuMP's user-defined nonlinear operator support: F(x) is registered as an operator, automatic differentiation supplies ∇F(x), and the Hessian-of-the-Lagrangian is computed by forming the composed predictor G(x, λ) = (Affine{λ,0} ∘ F)(x) and differentiating it once, avoiding the need to build P separate Hessians. These pieces together let the solver treat the predictor as a black-box oracle while still receiving exact derivative information.","core_discovery":"The discovery the paper puts forward is an engineering capability: a general-purpose, solver-independent modeling library can embed a trained predictor without requiring the predictor to be written out algebraically, as long as the modeling layer can provide oracle evaluations. The paper shows the full mechanism for JuMP and Ipopt, and it argues that the gray-box formulation is what makes very large neural networks tractable, because the optimization model's size scales with the predictor's inputs and outputs rather than with the number of layers or parameters. On the paper's own terms, MathOptAI.jl establishes that the gray-box technique, previously demonstrated only in custom implementations, can be packaged in a general library and combined with GPU acceleration through Python.","pith_inferences":["Inference: the gray-box design is portable—any modeling layer with user-defined nonlinear operators and any automatic-differentiation engine with a GPU backend could reproduce it, so the Julia-plus-Python combination in the paper is an implementation choice rather than a requirement.","Inference: the authors leave a numerical comparison of the three formulations to future work, so a natural next test is to benchmark full-, reduced-, and gray-box on the same set of problems to map where each wins in iterations, wall time, and solution quality.","Inference: the vector-only input and output convention, which the paper calls its biggest limitation, points to an obvious extension—native support for n-dimensional tensors such as images, which would remove the reshaping burden the library currently places on users.","Inference: in an iterative training-and-optimization loop, retraining a surrogate would not require rewriting the optimization model, only reloading the predictor, so closed-loop pipelines could become simpler to maintain."],"forward_implications":["Users can embed PyTorch, Flux, Lux, GLM, and DecisionTree models into a JuMP model and solve with any of the many solvers JuMP supports, not just one vendor's optimizer.","The gray-box formulation keeps the number of optimization variables and constraints tied to the predictor's input and output dimensions, so networks with very large parameter counts become embeddable; the authors point to related work demonstrating scales of hundreds of millions of parameters.","GPU-accelerated oracle evaluation is available for PyTorch models, moving function, Jacobian, and Hessian-of-the-Lagrangian computation to the GPU while the rest of the nonlinear oracles remain on the CPU.","Choosing a reformulation is a runtime option: the same embedded network can be solved with a nonlinear solver like Ipopt or, by changing the ReLU reformulation, with a mixed-integer solver like HiGHS.","Because predictors compose and third-party Julia packages plug in through package extensions, adding a new model type is a small amount of code, as demonstrated by the GLM adapter."],"supporting_citations":[{"why":"Provides JuMP, the modeling language whose user-defined nonlinear operator interface carries the full-space, reduced-space, and gray-box formulations.","marker":"Lubin et al. 2023"},{"why":"Provides Ipopt, the nonlinear solver that the gray-box formulation feeds with function, Jacobian, and Hessian-of-the-Lagrangian callbacks.","marker":"Wächter and Biegler 2006"},{"why":"Provides PyTorch, the source of the trained models and the automatic differentiation used to evaluate gray-box oracles, including on a GPU.","marker":"Paszke et al. 2019"},{"why":"Provides OMLT, the main prior general-purpose embedding toolkit and the baseline against which MathOptAI.jl's reduced-space and gray-box support are compared.","marker":"Ceccon et al. 2022"},{"why":"Demonstrated the gray-box embedding technique in a custom implementation, establishing the precedent that MathOptAI.jl packages into a general library.","marker":"Casas et al. 2025"},{"why":"Supplies the related numerical results on GPU-accelerated neural-network constraints that the paper cites for the scaling and performance benefits of gray-box.","marker":"Parker et al. 2025"},{"why":"Provides gurobi-machinelearning, whose decision-tree formulation and feature set are used as reference points in the comparison table.","marker":"Gurobi Optimization, LLC 2026a"}],"fun_headline_variants":["Embed ML predictors in JuMP models with MathOptAI.jl","Gray-box trick keeps JuMP model size independent of network layers","MathOptAI.jl: GPU offload for neural nets in JuMP","Optimization model size independent of ML parameter count"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The gray-box formulation is only correct if the Julia–Python bridge returns the function value and the function's first- and second-order derivatives in exactly the order and convention Ipopt expects; if the derivatives are misordered, mis-scaled, or computed incorrectly, the solver is optimizing a different problem than y = F(x).","fun_headline_variants_meta":{"raw":{"variants":["Embed ML predictors in JuMP models with MathOptAI.jl","Gray-box trick keeps JuMP model size independent of network layers","MathOptAI.jl: GPU offload for neural nets in JuMP","Optimization model size independent of ML parameter count"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001401,"raw_usage":{"total_tokens":5620,"prompt_tokens":855,"completion_tokens":4765,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":471,"completion_tokens_details":{"reasoning_tokens":4693}},"tokens_in":471,"tokens_out":4765,"duration_ms":34832,"temperature":1.0,"reasoning_tokens":4693,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T20:16:35.378230+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a small pretrained network with known weights, embed it in JuMP once in full-space and once in gray-box, and solve identical optimization problems with Ipopt; a disagreement between the two solutions beyond solver tolerance, or a mismatch between the gray-box Jacobian and Hessian-of-the-Lagrangian and finite-difference values at a random feasible point, would show that the oracle bridge does not represent y = F(x).","supporting_citations":[],"review_version":1}