8. What remains
Draft
The machinery of Sections 2 to 7 is what a global solver runs. This section lists what still stands between that machinery and a solver that finishes on nonconvex problems at the sizes practitioners bring. For each item it says what the evidence is and when it was read. The date matters. Several of the items are measurements that move between benchmark runs, and each statement carries its date, which is 5 October 2026 unless a sentence says otherwise. Eight items follow: the bound, the linear program inside the node, the shape of the tree, the arithmetic, the learned components, the formulation, the benchmarks themselves, and, as the counterweight, what is already solved and at what size. The last item is the bridge to Section 9. In the ordinary monthly trade the tax problem is a convex mixed-integer quadratic program that current solvers finish. When the wash-sale rule enters the model it becomes nonconvex.
The tax problem of Section 9, without and with the wash-sale rule
the tax problem in the the wash-sale rule
ordinary monthly trade: ----- enters the model -----> nonconvex
a convex mixed-integer
quadratic program that (the rule as Section 9 sketches
current solvers finish it, from memory and unverified)
Bounds are the bottleneck
(The benchmark evidence) A MILP solver's run usually ends with a proof. A nonconvex MINLP solver's run usually ends with a time limit, and the gap it then reports is the measure of what its relaxations could not see. The benchmark evidence is Mittelmann's MINLP table of 26 February 2026. It runs 200 MINLPLib instances through GAMS with a two-hour limit on an AMD Ryzen 9 5900X with 12 cores and 128 GB. The feasibility tolerance is \(10^{-6}\), failures are charged at the time limit, and every instance has to be solved globally. BARON solved 158 of the 200, SCIP 154, LINDO 116 and SHOT 96. Gurobi and Xpress were run as well but do not permit their columns to be published.H. D. Mittelmann, "Mixed Integer Nonlinear Programming Benchmark (MINLPLIB)", page dated 26 February 2026, plato.asu.edu/ftp/minlp.html, read 5 October 2026. The header lists BARON, GUROBI, LINDO, SCIP, SHOT and XPRESS with the footnote ": publication suppressed". The unscaled shifted geometric means of the four published columns are BARON 82.5, LINDO 398.1, SCIP 117.8 and SHOT 414.3 seconds. The 200 instances are the test set of the SCIP 8 paper (K. Bestuzheva, A. Chmiela, B. Müller, F. Serrano, S. Vigerske and F. Wegscheider, "Global optimization of mixed-integer nonlinear programs with SCIP 8", *Journal of Global Optimization 91 (2025), online December 2023; arXiv 2301.00587), to which the page links. SHOT handles only the quadratic instances of the set (Mittelmann, INFORMS Annual Meeting, 28 October 2025, slide 21, plato.asu.edu/talks/informs2025.pdf). Section 5.2 computes from the per-instance table behind the page that at least one of the four published solvers solved 183 of the 200, which is 91.5 per cent, and that 17 instances were solved by none of them within the two hours.H. D. Mittelmann, per-instance comparison table plato.asu.edu/ftp/compare.txt, last modified 6 March 2026, read 5 October 2026; the counting rule and the caveat on the LINDO column are in Section 5.2. The earlier version of this series rounded that union to "92%". The number is a computation from the per-instance files and not a figure on the page, so it is quoted here as 183 of 200 with the method in Section 5.2.
(The reason is the envelope) The reason is the envelope, the tightest convex function that lies below a given function on a box (Section 2.4), because a global solver's bound is the value of a relaxation built from envelopes, and the bound can be no better than the envelopes are tight. On the bilinear running example the McCormick relaxation claims \(0.40\) against a true maximum of \(0.18\), a gap of 122 per cent of that maximum. The Al-Khayyal–Falk theorem of Section 2.4 says why. Each of the four planes misses the product by \((x^U - x^L)(y^U - y^L)/4\) at the centre of the box, so the two envelopes are half the box's area apart there, and they agree with the product only on the box's edges.F. A. Al-Khayyal and J. E. Falk, "Jointly constrained biconvex programming", Mathematics of Operations Research 8 (1983); the planes are G. P. McCormick's, "Computability of global solutions to factorable nonconvex programs: Part I. Convex underestimating problems", Mathematical Programming 10 (1976). Only two things close that gap. Shrinking the box narrows the band by the square of the count: on a box of side \(1/k\) each plane misses the product by \(1/(4k^2)\) at the centre. The relaxed maxima the McCormick figure of Section 2.4 reports are \(0.400\), \(0.233\), \(0.193\), \(0.192\), \(0.187\) and \(0.185\) for \(k = 1\) to \(6\) boxes per axis, at the price of \(k^2\) boxes in two dimensions and \(k^d\) in \(d\). The bound itself does not fall by the square of the count, because the position of the constraint line relative to the grid decides where the relaxed maximum sits. Lifting closes the gap without branching when the structure allows. The full level-1 RLT of Section 4.7, the reformulation–linearization technique that multiplies pairs of the constraints and replaces each product of variables by a new variable, brings the bound to \(0.30\), and one convexity fact about the square, \(X \ge x^2\), brings it to \(0.18\), the exact answer. On Haverly's pooling problem the same picture appears at a scale a refinery would recognize. The McCormick relaxation of the p-formulation claims a profit of 500 on the base case, whose true optimum is 400. It does so because it lets the single pool carry 3 per cent sulphur toward one output and 1 per cent toward the other at once. The pq-formulation gives the same 500 on this instance, since with one pool and no pool capacity its RLT rows reproduce the McCormick rows. A single spatial branch on the pool quality at 2 per cent gives 400 on the lower half and 100 on the upper half, which proves the optimum. The second and third Haverly cases have root bounds of 1,000 against 600 and 800 against 750.C. A. Haverly, "Studies of the behavior of recursion for the pooling problem", ACM SIGMAP Bulletin 25 (1978). The bounds are those the pooling figure of Section 3.5 computes in the page; the p- and q-formulation vocabulary is that of M. Tawarmalani and N. V. Sahinidis, Convexification and Global Optimization in Continuous and Mixed-Integer Nonlinear Programming (Kluwer, 2002). The spatial figure of Section 3.5 shows the price of closing a gap by boxes alone. With midpoint splits it needs 13 nodes to a tolerance of 0.01, 21 to 0.001 and 27 to 0.0001, on a problem with two variables, where the uniform grid needs 25 boxes for the first of these.
Haverly's pooling problem, base case: one spatial branch
root: profit bound 500, true optimum 400
(the McCormick relaxation of the p-formulation;
the pq-formulation gives the same 500)
/ \
pool quality <= 2 per cent pool quality >= 2 per cent
/ \
the lower half: bound 400 the upper half: bound 100
at the root the relaxation lets the single pool carry 3 per cent
sulphur toward one output and 1 per cent toward the other at
once; the bounds of the two halves, 400 and 100, prove the
optimum 400
the second and third cases, root bound against the optimum:
1,000 against 600, and 800 against 750
(What the bound is made of: two ablations) What a global solver spends its time on is therefore the bound, and the ablations of one modern global solver, an ablation being a run with one component switched off and measured against the full solver, say which parts of the bound matter. FICO Xpress Global was tested by its authors on an internal set of 1,308 instances with three permutations each. Switching off presolve loses about 8 per cent of the solved instances and costs 40 per cent in time. Removing the RLT cuts loses 34 solved instances. Strong branching in this solver, the branching rule that evaluates each candidate variable by a trial branch before choosing one (Section 3.1), applies bound propagation and convexification cuts to each candidate before evaluating it. Removing the convexification cuts from that evaluation loses 583 solved instances, removing bound propagation loses 96, and removing both loses 689. Removing the SLP local heuristic, where SLP is sequential linear programming, a local NLP method that linearizes at the current point, loses 102, mostly on hard instances.P. Belotti, T. Berthold, T. Gally, L. Gottwald and I. Pólik, "Solving MINLPs to global optimality with FICO Xpress Global", Optimization Online, entry 2025/07 (published 31 July 2025, updated 4 August 2025), optimization-online.org/2025/07/solving-minlps-to-global-optimality-with-fico-xpress-global/, Sections 4.5 to 4.9 and Tables 5 and 6. Table 5 reads "Propagation −96", "Cuts −583", "Both propagation and cuts −689" in solved instances. The test set is described in its Section 3. Every one of those components except the last acts on the bound. The lesson usually drawn from twenty years of MILP ablations is that cutting planes and presolve are the two largest factors and heuristics a small one. It carries over to the nonconvex case with the envelope rows and their tightening in the role of the cuts.R. E. Bixby, M. Fenelon, Z. Gu, E. Rothberg and R. Wunderling, "Mixed-integer programming: a progress report", in The Sharpest Cut (MPS-SIAM, 2004); T. Achterberg and R. Wunderling, "Mixed integer programming: analyzing 12 years of progress", in Facets of Combinatorial Optimization (Springer, 2013). The two chapters are usually summarized as ranking cutting planes first and presolve second. Their per-component factors could not be read for this series, since neither chapter has an open copy, and Section 5.1 quotes only what could be read. The one open ablation with numbers, G. Mexi, The Two Faces of Mixed-Integer Programming: Primal and Dual Progress, doctoral thesis, Technische Universität Berlin (2026), doi 10.14279/depositonce-26588, read 5 October 2026, ranks presolve above cutting planes for SCIP 10 on 349 instances with five seeds: disabling presolve slows the solver by a factor of 2.60, random branching by 2.55, disabling cutting planes by 2.10, and disabling primal heuristics or conflict analysis by about 1.3.
| part removed | solved instances lost |
|---|---|
| bound propagation, in strong branching's evaluation of each candidate | 96 |
| convexification cuts, in strong branching's evaluation of each candidate | 583 |
| both of the two, in strong branching | 689 |
| RLT cuts | 34 |
| SLP local heuristic | 102, mostly on hard instances |
| presolve | about 8 per cent of the solved instances; costs 40 per cent in time |
Reformulation moves the bound more than search does, and Section 8.6 collects the evidence from Section 4. There is no general procedure for finding the good formulation. The 17 unsolved instances of the benchmark are the instances for which nobody has found one, or for which none is known to exist.
(The bound on a GPU, as of 5 October 2026) For the GPU the state of the bound is the first item on the agenda of Section 7.8, and as of 5 October 2026 it stands as follows. Two bounding engines for deterministic global optimization have run on a GPU. The interval lower bounder of MAiNGO partitions a node's box into \(m^d\) sub-boxes, evaluates the natural interval extension and the mean value form on each in one kernel, and takes the minimum. Its authors report bounds that tighten with the partition, fewer branch-and-bound iterations, and wall-clock times three orders of magnitude below CPU interval arithmetic without partitioning.H. Zhang, T. Kerkenhoff, N. Kichler, M. Dahmen, A. Mitsos, U. Naumann and D. Bongartz, "Accelerating deterministic global optimization via GPU-parallel interval arithmetic", arXiv 2507.20769 (2025). The mean value form's gap falls as the square of the sub-box width (Section 6.4; Section 7.8, problem 2), which is why a cheap partition turns a useless bound into a sharp one; the exponent \(d\) is why it stops being cheap beyond a handful of nonconvex variables per node. The pointwise McCormick kernels generated by SourceCodeMcCormick.jl evaluate relaxations, interval extensions and subgradients for thousands of boxes at once. The branch and bound ParBB built on them is reported at about 9 nanoseconds per relaxation evaluation on the GPU against 237 on the CPU, and 11 to 22 times faster than EAGO on the authors' test problems. These are the authors' numbers, and the supported operation set was small when the package was read.R. X. Gottlieb, P. Xu and M. D. Stuber, "Automatic source code generation for deterministic global optimization with parallel architectures", Optimization Methods and Software 41 (2026), online September 2024; software at github.com/PSORLab/SourceCodeMcCormick.jl. The timings are recorded from the package's README; the paper's own tables were not read for this series. Neither engine has been evaluated on MINLPLib under a benchmark protocol, and no GPU spatial branch and bound for general factorable mixed-integer problems exists. Everything else a GPU can do for a tree, which Sections 7.2 to 7.7 described, is downstream of a bound that is tight enough to prune.
The price grows exponentially in \(d\) while the gain does not depend on \(d\), which is the whole trade: a cheap partition makes a useless bound sharp on a node with two or three nonconvex variables, and stops being cheap beyond a handful.
Two GPU bounding engines: many boxes at once
MAiNGO's interval lower bounder: the sub-boxes of one node
the node's box one kernel, on every sub-box:
partitioned into --------> the natural interval extension
m^d sub-boxes and the mean value form
[][][][] ... [] |
v
the minimum over the sub-boxes:
the node's lower bound
SourceCodeMcCormick.jl: thousands of boxes at once
[][][][][][] ... [][] --> the generated McCormick kernels:
relaxations, interval extensions
and subgradients
ParBB, built on them: about 9 ns per relaxation evaluation on
the GPU against 237 on the CPU (the authors' numbers)
the mean value form's gap falls as the square of the sub-box
width; the exponent d is why the partition stops being cheap
beyond a handful of nonconvex variables per node
The LP is still sequential
(The warm-started simplex, and the batched alternative) The relaxation a MILP solver solves at a node is a linear program that differs from its parent's by one bound. The dual simplex method starts from the parent's optimal basis, the set of columns that determines the parent's optimal vertex, which is dual feasible for the child, meaning that the reduced costs keep their signs and only the primal side is upset by the changed bound, and reaches the child's optimum in a few pivots. Each pivot is a sparse LU update, an incremental change to the triangular factors of the basis matrix, that depends on the last, and the whole solve runs on one core. This is the right tool for the thousands of small, similar LPs a tree needs, and nothing built on first-order methods has displaced it inside the tree. Blin, Gualandi, Maes, Lodi and Stellato put the question in the form a GPU needs. A batch of \(K\) node LPs that share the constraint matrix and differ in their bounds is a matrix of primal iterates \(X \in \mathbb{R}^{n \times K}\). One step of PDHG on all of them is two sparse-matrix–dense-matrix products. They identify the problem sizes at which that batch beats a sequence of dual simplex solves, and the sizes at which it does not.N. Blin, S. Gualandi, C. Maes, A. Lodi and B. Stellato, "Batched first-order methods for parallel LP solving in MIP", arXiv 2601.21990 (2026). Their experiments are strong branching and bound tightening at the root, which is where cuOpt applies the same idea (mip_batch_pdlp_strong_branching, and from release 26.04 an option for batched PDLP in reliability branching; NVIDIA cuOpt release notes, docs.nvidia.com/cuopt/user-guide/latest/release-notes.html, read 5 October 2026).
One node LP on one core, or K node LPs as one matrix
dual simplex: the child from the parent's optimal basis, which
is dual feasible for the child; a few pivots, on one core
parent's basis --> pivot --> pivot --> ... --> child's optimum
each pivot a sparse LU update that depends on
the last
batched PDHG: K node LPs that share the constraint matrix and
differ in their bounds
node LP 1 node LP 2 node LP K
+-----------+-----------+-- ... --+-----------+
| | | | |
X | primal | primal | | primal | n
| iterate | iterate | | iterate | rows
| | | | |
+-----------+-----------+-- ... --+-----------+
K columns
one PDHG step on all K: two sparse-matrix-dense-matrix products
(The accuracy trade, measured) The trade is in the accuracy. A first-order iterate that is \(10^{-4}\) infeasible is not a solution, and a dual iterate that is \(10^{-4}\) from dual feasibility is not a bound. The repair is the safe bound of Section 7.3: for any \(y \ge 0\) the Lagrangian dual function with the box kept as a hard constraint is a valid bound, loosened by the bound-weighted dual infeasibility.A. Neumaier and O. Shcherbina, "Safe bounds in linear and mixed-integer linear programming", Mathematical Programming 99 (2004). What the repair costs can be measured. Take a frontier of 256 node relaxations of a 40-item knapsack LP, on which an exact LP would prune 230 of the nodes against the incumbent. After 10 PDHG iterations the safe bound already prunes 228 of them, and its worst looseness is 23.5 objective units. After 50 iterations it prunes all 230 with a worst looseness of 2.56, while the primal iterates are still infeasible by 0.34 on average. After 200 iterations the looseness is 0.26 and the infeasibility 0.03.Mode A of the frontier listing of Section 7.4 (batched_pdhg_frontier.py), whose table prints 228 pruned with a worst looseness of 23.47 at iteration 10, 230 pruned with 2.56 and a mean primal infeasibility of 0.337 at iteration 50, and 230 with 0.26 and 0.030 at iteration 200. The PDHG figure of Section 7.2 shows the same correction on the running example: at iteration 120 the raw dual value 4.03 is not a bound, the safe bound 6.21 is one, and the true value is 5.56. Most pruning decisions need very low accuracy and a few need near-exact solves, and the few are the nodes whose bound is within the gap. The tolerance figure of Section 7.4 shows what happens to the tree when every bound is loosened by a relative \(\delta\). On its 12-item knapsack the exact search takes 37 nodes. The search with \(\delta = 10^{-4}\) takes 39, with \(\delta = 10^{-2}\) 81, and with \(\delta = 10^{-1}\) 1,865. Loosening the bound enlarges the tree, as these numbers show. On this instance, and on most of the figure's instances, the growth is faster than \(\delta\) itself once \(\delta\) approaches one integer unit of the objective.
(What the solvers' tolerances are) A solve to \(10^{-4}\) serves a heuristic and not a proof. cuOpt's PDLP solves to \(10^{-4}\) relative accuracy by default and its barrier method to \(10^{-8}\). Gurobi 13's PDHG stops at its own tolerances and hands the point to crossover, the step that moves an interior point to a vertex and a basis from which the simplex can finish exactly. The condensed interior-point methods of Section 7.5 run at \(10^{-4}\) to \(10^{-6}\).NVIDIA cuOpt User Guide, "LP and MILP settings" (release 25.12 and later): "PDLP solves to 1e-4 relative accuracy by default. Barrier solves to 1e-8 relative accuracy by default. Dual Simplex solves to 1e-6 absolute accuracy by default." Gurobi Optimizer Reference Manual 13.0, parameters Method (value 6, PDHG), PDHGGPU, PDHGRelTol, PDHGAbsTol, PDHGConvTol, Crossover, docs.gurobi.com/projects/optimizer/en/current/reference/parameters.html; the 13.0.2 release notes say the GPU PDHG "has matured from the beta state and now is a fully supported feature". S. Shin, M. Anitescu and F. Pacaud, "Accelerating optimal power flow with GPUs: SIMD abstraction of nonlinear programs and condensed-space interior-point methods", Electric Power Systems Research 236 (2024), 110651; F. Pacaud, S. Shin, A. Montoison, M. Schanen and M. Anitescu, "Condensed interior-point methods for scalable nonlinear programming on GPUs", Mathematical Programming Computation (2026, online August 2026); arXiv 2405.14236. Mittelmann's GPU QP benchmark shows the cost of asking for more. On 21 QPLIB quadratic programs on a B200, cuOpt solves 19 at tolerance \(10^{-4}\), 16 at \(10^{-6}\) and 13 at \(10^{-8}\), and its scaled time grows from 1 to 3.48 to 8.29.H. D. Mittelmann, "Benchmarking Optimization Software: a (Hi)Story", slides of a 2019 EURO talk updated and extended for a talk at PolyU Hong Kong, March 2026, slide 35 (GPU QP benchmark of 18 February 2026; HPR-QP on the same table: 19, 15 and 15 solved at scaled 3.05, 8.89 and 10.1), plato.asu.edu/talks/hongkong26.pdf. The file is hongkong26.pdf; the address printed on the slides' last page, hongkong2026.pdf, returns no document. For a heuristic that is fine. For a proof it is not, because a node bound that is off by \(10^{-4}\) of a large objective can exceed the gap being proved.
(cuOpt, and the MIPFEAS benchmark) The one general solver built around a GPU shows where the frontier is. cuOpt runs its primal heuristics on the device and its branch and bound on the CPU. Its documentation says of the MIP solver that it "excels at finding high-quality feasible solutions quickly" while "proving feasible solutions optimal remains under active development".NVIDIA cuOpt User Guide 26.08, "Introduction", docs.nvidia.com/cuopt/user-guide/latest/introduction.html, read 5 October 2026: the MIP solver combines "GPU-accelerated primal heuristics for improving the primal bound with traditional CPU algorithms, including branch and bound, to improve the dual bound", and "primal heuristics such as local search, feasibility pump, and feasibility jump run on the GPU". The repository github.com/NVIDIA/cuopt is Apache-2.0, created 8 April 2025; Papilo-based presolve is on by default from release 25.10 (October 2025), root cuts from 26.02 (February 2026), work-stealing branch-and-bound workers from 26.06. The MIPFEAS benchmark of March 2026 measures exactly that split. On the 233 feasible instances of the MIPLIB 2017 benchmark set, with 600 seconds per instance, cuOpt 26.02 on an H100 found a feasible solution on 225 instances and proved 57 optimal, with a time-averaged primal integral of 0.0651, the primal integral being the area under the incumbent's gap against time, which Section 2.3 defines and Section 8.7 returns to, so that a lower value means good solutions found sooner. HiGHS 1.12 found 212 and proved 97 with 0.0567, and SCIP 10.0.1 found 204 and proved 77 with 0.0755.The figures are those of slide 36 of Mittelmann's Hong Kong talk (March 2026): cuOpt 26.02 on an H100, HiGHS 1.12.0, SCIP 10.0.1 with SoPlex, CBC 2.10.11, all through GAMS 53.1 on the 233 feasible MIPLIB 2017 instances; the "virtual mean commercial solver" on the same slide found 227 and proved 180 with 0.0132. The benchmark is M. Bussieck and S. Dirkse, "Expanding the focus: introducing the MIPFEAS benchmark", GAMS blog, posted 17 March 2026 and updated through 16 September 2026, gams.com/blog/2026/03/expanding-the-focus-introducing-the-mipfeas-benchmark/, read 5 October 2026; its tables have been revised several times and differ from the slide in solver versions and in the virtual solvers' counts. Its table of 16 September 2026 (600 s, 24 threads) reads cuOpt 26.08 on a B200 226 found, 78 proved, 0.0286; HiGHS 1.15.1 214, 107, 0.0453; SCIP 10.0.1 204, 77, 0.0755; virtual mean commercial solver 227, 183, 0.0105. The integral is the time average \(\bar P(T)\) of Section 2.3 with the benchmark's own integrand. Mittelmann's MIPFEAS page of 16 September 2026, plato.asu.edu/ftp/mipfeas.html, prints it as 2 if \(z(t) = \infty\), 1 if \(z(t)\,z^\star < 0\), and otherwise \(|z(t) - z^\star| / \max(|z(t)|, |z^\star|, 1)\). The post's own convention is Berthold's, stated in Sections 2.3 and 7.7, which uses 1 in both of the first two cases. A MIPFEAS integral is therefore never smaller than the post's for the same run, and it is larger whenever the run spends time without an incumbent. The three primal integrals are quoted in the benchmark's own convention, which scores the time before the first incumbent at 2 where the post's convention, Berthold's, scores it at 1. The GPU side finds more solutions than any open-source CPU solver and proves fewer: 57 against 97 for HiGHS and 77 for SCIP. On the LP side the picture is better and still not decisive. On Mittelmann's LP feasibility benchmark of 16 September 2026 the best GPU code, HPR-LP-C, has scaled shifted geometric mean 1 with 63 of 65 instances solved. cuOpt follows at 1.17 with 62 and COPT's GPU barrier at 1.29 with 64. The CPU codes are COPT at 1.67 with 65 and MOSEK at 5.90 with 56. The decisive wins are on instances with \(10^8\) nonzeros, where the CPU codes do not finish and the GPU codes run out of memory before they run out of time.H. D. Mittelmann, "LPfeas Benchmark (find a PD feasible point; also for GPUs)", page dated 16 September 2026, plato.asu.edu/ftp/lpfeas.html: B200 for the GPU codes with a 1,000 s limit against 15,000 s for the CPU codes. On the six additional huge instances of the October 2025 talk (up to 126 million variables and 253 million nonzeros) cuPDLP-C was the only GPU code to solve all six, in 1,166 to 7,288 s; cuOpt failed on all six for memory (Mittelmann, INFORMS 2025 slides, slides 12 and 23; H. Lu, J. Yang, H. Hu, Q. Huangfu, J. Liu, T. Liu, Y. Ye, C. Zhang and D. Ge, "cuPDLP-C: a strengthened implementation of cuPDLP for linear programming by C language", arXiv 2312.14832 (2023)).
| solver | feasible solution found | proved optimal | time-averaged primal integral |
|---|---|---|---|
| cuOpt 26.02 (H100) | 225 | 57 | 0.0651 |
| HiGHS 1.12 | 212 | 97 | 0.0567 |
| SCIP 10.0.1 | 204 | 77 | 0.0755 |
(A division of labour) For a tree these numbers suggest a division of labour. Node relaxations that share a matrix batch well, and siblings share one. A column of the batch can be retired the moment its safe bound crosses the incumbent, which for most columns is early. The columns that cannot be retired are the ones that need accuracy, and those are the ones a dual simplex with a warm start solves best. No published rule says when to stop a column, and no measurement on MINLPLib says how much larger the tree becomes when every bound is a safe bound from a loose dual. That is the fourth item of the agenda of Section 7.8, and it is open.
A division of labour in the tree, as the numbers suggest
node relaxations that share a matrix (siblings share one)
|
v
+----> a PDHG step on every column of the batch
| |
| v
| has the column's safe bound --yes--> the column is
| crossed the incumbent? retired: for most
| | no columns, early
| v
+--yes-- iterate on? (no published rule says when to stop)
| no
v
the columns that cannot be retired, the ones that need accuracy:
a dual simplex warm-started from the parent's basis, which
reaches the child's optimum in a few pivots
open, item 4 of the agenda of Section 7.8: no measurement on
MINLPLib of how much larger the tree becomes when every bound
is a safe bound from a loose dual
Trees are irregular
(What the shape costs) The shape of a branch-and-bound tree is unknown until the search is over, and the efficiency of running it in parallel is bounded by that shape. The parallel figure of Section 6.3 shows the three phases, ramp-up, a primary phase and ramp-down, on a knapsack of about a thousand nodes. Eight workers with incumbents shared reach 7.86 times the speed of one at 98 per cent efficiency, and sixty-four workers reach 42.63 times with incumbents shared and 26.51 without, because each worker then explores what the others' solutions would have pruned. Lai and Sahni showed in 1984 that adding processors can make a best-first branch and bound slower, or more than proportionally faster, depending on which node is evaluated when. Li and Wah gave the conditions under which the detrimental anomaly cannot occur, which amount to distinct bounds, so that the best-first order is a total order (Theorem 6.2.5). These are the anomaly theorems of Section 6.2.T.-H. Lai and S. Sahni, "Anomalies in parallel branch-and-bound algorithms", Communications of the ACM 27 (1984); G.-J. Li and B. W. Wah, "Coping with anomalies in parallel branch-and-bound algorithms", IEEE Transactions on Computers C-35 (1986).
The three phases of a parallel run (the figure of Section 6.3)
workers
busy
^ primary phase
| +-------------------------------+
| / \
| / \
| / \
+------+---------------------------------------+------> time
ramp-up ramp-down
the figure's knapsack of about a thousand nodes, speed against
one worker:
8 workers, incumbents shared 7.86 times, 98 per cent
efficiency
64 workers, incumbents shared 42.63 times
64 workers, incumbents not shared 26.51 times: each worker
explores what the others'
solutions would have pruned
(The measured scaling is modest) The measured scaling of real solvers is modest. Koch, Ralphs and Shinano observed in 2012 that the average speedup of MIP solvers from 1 to 12 threads on MIPLIB 2010 was roughly a factor of 3.T. Koch, T. Ralphs and Y. Shinano, "Could we use a million cores to solve an integer program?", Mathematical Methods of Operations Research 76 (2012). The one thread-scaling table published for a global MINLP solver is Xpress Global's, which Section 6.3 reproduces: about 2.6 times on sixteen threads over all solvable instances, about 4 times on the instances that needed at least 100 nodes, and no further gain from 32 threads. The deterministic MILP frameworks of 2025 and 2026 report the same shape, Para-B&B for HiGHS reaching 2.17 times at eight threads with its threads idle a third of the time (Section 6.3). ParaSCIP's record of twenty-one previously unsolved MIPLIB instances, with up to 80,000 cores for one run (Section 6.3), is a record of what sustained effort over seven years and several machines can do for individual instances. It is not a scaling curve.
(The GPU results are on other trees) The GPU results are on trees unlike a MINLP's. Gmys's flow-shop solver and the Chapel codes of Section 6.4 solved open permutation flow-shop instances on hundreds of GPUs and reach 44 per cent efficiency relative to one node at 1,024 GPUs, on N-Queens and flow-shop trees. Those trees have bounds that cost a few hundred arithmetic operations, and they are as regular as trees get. A spatial branch-and-bound node is an LP of thousands of rows whose bound tightening alone may cost two LPs per variable. The frontier-batch pattern of Section 6.4 inherits a price that is independent of the problem: a batch committed before the incumbents it will itself produce are known evaluates nodes a serial search would have pruned. In the batch figure of Section 6.4 the best batch size for its device gives a speedup of 9.30 while wasting 116 of its 1,145 bounded nodes, and a larger batch is slower again because most of its rounds are not full. Lane efficiency is the other cost. A round that waits for its slowest node keeps its lanes busy for a fraction of the time, and that fraction falls with the batch size and with the spread of node work.
A frontier batch, one round: where the time and the nodes go
lane 1 [######]....................
lane 2 [#############].............
lane 3 [###].......................
lane 4 [##########################] the slowest node
|<-------- one round ----->| then the barrier
# bounding a node . idle, waiting for the slowest
idle lanes: the busy fraction of a round falls with the batch
size and with the spread of node work
wasted nodes: the batch is committed before the incumbents it
will itself produce are known, so it bounds nodes a serial
search would have pruned
the batch figure of Section 6.4, at the best batch size for its
device: a speedup of 9.30, with 116 of its 1,145 bounded nodes
wasted; a larger batch is slower again, most of its rounds not
full
(Two absences, dated) Two absences have to be stated as absences, with their date. As of 5 October 2026 no published study measures the strong scaling of a global MINLP solver across more than one node, in processes or in devices, on MINLPLib or QPLIB instances. The nearest evidence is Xpress Global's one-node table, FiberSCIP's shared-memory experiments, ParaSCIP's MILP records and the Chapel results on combinatorial trees. And no benchmark measures any GPU code on MINLPLib instances. Mittelmann's three GPU tracks are LP feasibility, convex continuous QPLIB and sparse SDP, the MIPFEAS benchmark covers MILP heuristics, and his MINLP page has no GPU column and no thread-scaling arm.H. D. Mittelmann, "Benchmarks for Optimization Software", index page with entries dated through 1 October 2026, plato.asu.edu/bench.html, read 5 October 2026, labels three tracks "also for GPUs". Y. Shinano, S. Heinz, S. Vigerske and M. Winkler, "FiberSCIP — a shared memory parallelization of SCIP", INFORMS Journal on Computing 30 (2018), reports MIP and MINLP results and a deterministic mode. The protocol under which the two gaps could be filled is item 12 of the agenda of Section 7.8. One paper would falsify either statement, which is why both carry a date.
(Determinism on a CPU) Determinism needs two statements, one for each kind of processor. On a CPU, determinism is a documented product property with a known mechanism and a measured price, and Sections 6.2 and 6.3 give the record. Gurobi is deterministic by default, in parallel as well, with the time-dependent parameters, the default LP method and the concurrent MIP mode as the documented exceptions, and its work limit stops a run at the same point every time. Xpress runs its tree search deterministically by default, its authors having measured the opportunistic alternative at 2 to 9 per cent faster over 2 to 32 threads (Section 6.2). CPLEX's default automatic mode applies as much parallelism as possible while still achieving deterministic results (Section 6.3). FiberSCIP's deterministic mode counts communication points and circulates a token in a fixed order, and ParaSCIP is opportunistic by design. The mechanism in every case is a logical clock, as in Section 6.2. Workers compute without communicating until a counter of work reaches a common value, then exchange in a fixed order, and a wall-clock limit is the one thing that breaks it.
Determinism on a CPU: a logical clock, not the wall clock
wall-clock time ------------------------------------------->
worker 1 [==== W units of work ====].....|1|[==== W units ...
worker 2 [=== W units of work ===].......|2|[=== W units ...
worker 3 [====== W units of work ======].|3|[====== W units ...
^
the exchange, in the fixed order
1, 2, 3
W: the common value that each worker's counter of work reaches
before it communicates. The exchanges fall at the same counts of
work in every run, whatever the wall-clock times; a wall-clock
limit is the one thing that breaks it
(Determinism on a GPU: the arithmetic) On a GPU the sentence "the same answer twice" is not true today, for two reasons specific to the device. The first is arithmetic. Floating-point addition is not associative. A reduction accumulated into one word by many threads with atomicAdd is a left comb (recursive summation, Section 6.2), and its leaf order is the order in which the memory system serializes the threads. That order is not specified. NVIDIA's own floating-point guide draws the conclusion: "seemingly identical computations can produce different results even if all basic operations are computed in compliance with IEEE 754."NVIDIA, Floating Point and IEEE 754, CUDA 13.4 documentation, page dated 13 September 2026, docs.nvidia.com/cuda/floating-point/index.html; NVIDIA, CUDA Programming Guide 13.4.2, appendix "C++ language extensions", section "Atomic functions". The worst-case bound for two summation orders is Higham's, \(2(n-1)u\sum_i |x_i|\) to first order: N. J. Higham, "The accuracy of floating point summation", SIAM Journal on Scientific Computing 14 (1993). The listing of Section 6.2 measures the spread on one vector of a thousand doubles summed in a thousand arrival orders. The recursive sum takes 273 distinct values, a fixed pairwise tree whose leaves follow the arrival order still takes 45, the same tree over the operands in a canonical order takes one, and accumulation into a scaled 64-bit integer, the single-bin case of Demmel and Nguyen's reproducible summation (Proposition 6.2.19), takes one in every order.J. Demmel and H. D. Nguyen, "Parallel reproducible summation", IEEE Transactions on Computers 64 (2015); W. Ahrens, J. Demmel and H. D. Nguyen, "Algorithms for efficient reproducible floating point summation", ACM Transactions on Mathematical Software 46 (2020), whose abstract gives the cost of the full method: one read-only pass, one parallel reduction, and about \(9n\) floating-point and \(3n\) bitwise operations for a six-word accumulator. The fixed tree costs no extra arithmetic. It costs the freedom to accumulate in any order, so a kernel that scattered partial sums into one word with atomics needs a segmented layout or a second pass. A pruning test that compares a reduced bound with an incumbent to the last bit therefore depends on the order unless the schedule is fixed, and the dual objective \(b^\top y\) at a near-optimal dual, a cancelling sum with terms of both signs, is where the spread is largest (Theorem 6.2.17). The pairwise tree in canonical order is the schedule a reduction kernel must implement to be reproducible, and it parallelizes as a balanced tree of depth ten.
Summing on a GPU: the order of the additions decides the bits
atomicAdd into one word: a fixed pairwise tree,
a left comb, its leaves in the its leaves in arrival order
order in which the memory or in a canonical order
system serializes the threads
(an order not specified)
+ +
/ \ / \
+ x + +
/ \ / \ / \
+ x + + + +
/ \ / \ / \ / \ / \
+ x x x x x x x x x
/ \
x x
a thousand doubles summed in a thousand arrival orders (the
listing of Section 6.2): distinct results
recursive sum, the comb ............................. 273
pairwise tree, leaves in arrival order .............. 45
pairwise tree, operands in a canonical order ........ 1
accumulation into a scaled 64-bit integer, the
single-bin case of Demmel and Nguyen ................ 1
over a thousand operands the canonical tree has depth ten
(Determinism on a GPU: the batch) The second source is the composition of a batch. Which nodes form a round, and whether an incumbent found during the round is applied before the round ends, depend on arrival times unless the schedule is fixed. Two results of Section 6 separate what is at stake. Proposition 6.4.5 says that a fixed batch schedule, with canonical batches and identifiers, fixed-schedule bounds, incumbents merged at round boundaries and no pruning on the device, makes the state after every round a function of the input and of \(b\) alone, whatever the number of lanes, blocks and devices. Proposition 6.2.21 says that under directed rounding every reduction order gives a valid bound, so pruning is correct under any order of accumulation, while reproducibility needs the fixed schedule.
A fixed batch schedule (Proposition 6.4.5): rounds and barriers
round r round r + 1 round r + 2
----|==================|==================|==================|-->
^ ^ ^ ^
the barriers: an incumbent found in a round is merged at the
round boundary, and prunes from the next round on
each round: canonical batches and identifiers, fixed-schedule
bounds, no pruning on the device
the state after every round is a function of the input and of b
alone, whatever the number of lanes, blocks and devices
(What is reproducible today) The two propositions separate a correctness requirement from a product requirement. Directed rounding, which CUDA offers as per-operation intrinsics, makes pruning correct whatever the hardware does with the order, at the cost of one intrinsic per operation. A fixed schedule makes the run reproducible, at the cost described above. The price of the fixed schedule in a tree has three parts: a barrier per round, idle lanes while the slowest node of the round finishes, and the one-round lag before an incumbent prunes. What is reproducible on a GPU today is this. A single kernel with fixed-schedule reductions is bitwise reproducible on one device, and cuBLAS documents that its routines are, on GPUs with the same architecture and the same number of streaming multiprocessors, and that the guarantee is lost across toolkit versions and when several streams share one workspace (Section 6.2). cuOpt's deterministic parallel branch-and-bound mode, added in release 26.02, is experimental and excludes the GPU heuristics. No GPU MINLP solver with a deterministic mode exists.NVIDIA, cuBLAS Library 13.4, Section 2.1.4 "Results reproducibility", page dated 13 September 2026: "all cuBLAS API routines from a given toolkit version, generate the same bit-wise results at every run when executed on GPUs with the same architecture and the same number of SMs", lost "when multiple CUDA streams are active". NVIDIA cuOpt release notes, entry 26.02: "experimental support for determinism in the parallel branch-and-bound solver. GPU heuristics are not supported yet in this mode." MAiNGO's documentation states the validity half for its GPU interval bounder: "Due to the outward rounding, interval arithmetic will keep the true lower bound enclosed in the interval no matter which precision is used" (avt-svt.pages.rwth-aachen.de/public/maingo/special_uses.html, read 5 October 2026). The CPU standard, a deterministic path on the same machine with the same thread count, is not yet available for any GPU branch and bound. Its cost in wall-clock time and tree size is item 8 of the agenda of Section 7.8, and it has not been measured by anyone.
Numerics
Every solver works to tolerances. There are three, and they are different in kind. A feasibility tolerance says how far a constraint may be violated. An integrality tolerance says how far from an integer an integer variable may sit. A gap tolerance says how close the bound must come to the incumbent before the search stops. The defaults, read from the documentation and source code of the solvers on 4 and 5 October 2026, are these.Gurobi Optimizer Reference Manual 13.0, parameters FeasibilityTol, IntFeasTol, MIPGap, MIPGapAbs, and the Help Center article "What is the MIPGap?"; IBM ILOG CPLEX Optimization Studio 22.1.2 parameter pages EpRHS, EpInt, EpGap, EpAGap; FICO Xpress Optimizer controls FEASTOL, MIPTOL, MIPRELSTOP, MIPABSSTOP; SCIP source src/scip/def.h and src/scip/set.c on the master branch (SCIPsetIsFeasIntegral tests integrality against feastol; limits/gap 0; numerics/epsilon \(10^{-9}\)); BARON User Manual, version 2026.9.10, and the GAMS/BARON page, which sets EpsR and EpsA from GAMS optCR (default \(10^{-4}\)) and optCA (default 0); BARON options at dev.ampl.com/solvers/baron/options.html; NVIDIA cuOpt User Guide, "MIP settings". Couenne's default gap is 0.
| solver | feasibility | integrality | gap (relative; absolute) | gap formula |
|---|---|---|---|---|
| Gurobi 13 | FeasibilityTol 1e-6 | IntFeasTol 1e-5 | MIPGap 1e-4; MIPGapAbs 1e-10 | \(\vert z_{\mathrm{inc}} - \mathrm{bound} \vert / \vert z_{\mathrm{inc}} \vert\) |
| CPLEX 22.1 | EpRHS 1e-6 | EpInt 1e-5 | EpGap 1e-4; EpAGap 1e-6 | \(\vert \mathrm{bound} - z_{\mathrm{inc}} \vert / (10^{-10} + \vert z_{\mathrm{inc}} \vert)\) |
| Xpress | FEASTOL 1e-6 | MIPTOL 5e-6 | MIPRELSTOP 1e-4; MIPABSSTOP 0 | \(\vert z_{\mathrm{inc}} - \mathrm{bound} \vert \le \mathrm{tol} \cdot \max(\vert \mathrm{bound} \vert, \vert z_{\mathrm{inc}} \vert)\) |
| SCIP 10 | feastol 1e-6 | feastol 1e-6 (shared) | limits/gap 0; limits/absgap 0 | \(\vert p - d \vert / \min(\vert p \vert, \vert d \vert)\); infinite if signs differ |
| BARON 2026 | AbsConFeasTol 1e-6 | AbsIntFeasTol 1e-5 | EpsR 1e-6; EpsA 1e-6 standalone; under GAMS: optCR 1e-4, optCA 0 | \(\vert U - L \vert \le \mathrm{EpsR} \cdot \max(\vert L \vert, \vert U \vert)\) |
| cuOpt 26.08 | MIP absolute 1e-6 | 1e-5 | 1e-4; 1e-10 | \(\vert z_{\mathrm{inc}} - \mathrm{bound} \vert / \vert z_{\mathrm{inc}} \vert\); PDLP itself 1e-4 relative |
(Three readings of the table) Three readings of the table matter. First, "optimal" from a MILP or MINLP solver means a relative gap of \(10^{-4}\) unless the user changes it, with two exceptions. SCIP's gap limit is zero, so it finishes when its bound meets the incumbent to within its \(10^{-9}\) epsilon and the constraints hold to \(10^{-6}\). Standalone BARON stops at a relative or absolute gap of \(10^{-6}\), but under GAMS it inherits optCR and stops at \(10^{-4}\) like the others. Second, the gap formulas differ in the denominator, so "a 1 per cent gap" is not one number: Gurobi, CPLEX and cuOpt divide by the incumbent, SCIP by the smaller of the two bounds, Xpress and BARON by the larger. Third, the integrality tolerance is \(10^{-5}\) in Gurobi, CPLEX, BARON and cuOpt, \(5 \times 10^{-6}\) in Xpress and \(10^{-6}\) in SCIP, where it is the feasibility tolerance under another name.
(The tolerance is part of the problem statement) The tolerance is part of the problem statement, not a detail of the implementation, and the tax problem shows why in two ways. Nonlinear optima are irrational in general. The minimum of \((x^2 - 2)^2\) on \([0, 2]\) is at \(\sqrt 2 = 1.41421356\ldots\), and no finite decimal is exact (Proposition 5.6.2), so some \(\varepsilon > 0\) has to be in the specification before "solved" means anything. Compare the linear problem \(\min\{x : 3x \ge 2\}\), whose optimum \(2/3\) is exactly representable in rationals. And a relative tolerance is a sum of money. A relative gap of \(10^{-4}\) on an account whose after-tax objective is of the order of 720,000 dollars is 72 dollars. That is more than the 20.40 dollars by which the rules of thumb miss the optimum at 55,000 dollars of cash in Section 9.4. A tax is computed to the cent. The right stopping rule for the tax problem is an absolute gap of a few dollars, which is a few times \(10^{-6}\) relative on that account and would be a few times \(10^{-9}\) on an account a thousand times larger. Dividing the objective by \(10^5\) leaves a relative gap tolerance unchanged and makes an absolute one \(10^5\) times looser in dollars, so the stopping rule has to be stated in dollars and kept in dollars when the model is scaled.
(Tolerances compose with the model) The tolerances compose with the model, and the composition is where real errors come from. The big-M example of Section 5.6 is the case, a big-M constraint being a row \(x \le Mz\) in which a binary \(z\) switches a continuous \(x\) off through a large constant \(M\): the indicator constraint \(x \le M z\) with \(M = 10^6\) accepts \(z = 10^{-6}\) as the integer 0 and permits \(x = 1\), a position of one share that the solver believes is zero, and an absolute feasibility tolerance on a row whose coefficients have been divided by \(10^6\) is a million times looser in the original units. The remedy is a bound on \(x\) as small as the model allows, or an indicator constraint handled by the solver, rather than a smaller tolerance. A wash sale either occurred or did not, and a lot is sold whole or not at all. A model that returns a fractional lot, or a trade that violates a rule by a rounding error, has produced a wrong answer, not an approximate one. The script below prints the table's two nonzero relative gap tolerances, \(10^{-4}\) and \(10^{-6}\), and SCIP's \(10^{-9}\) epsilon in dollars on the account above.
The big-M example of Section 5.6: a tolerance composed with M
the indicator constraint x <= M z, with M = 10^6
the relaxed value z = 10^-6
|
| the integrality test
v
accepted as the integer 0
|
| the row: x <= M z = 10^6 * 10^-6 = 1
v
x = 1 permitted: a position of one share that the solver
believes is zero
and on a row whose coefficients have been divided by 10^6, an
absolute feasibility tolerance is a million times looser in the
original units
# Relative gap tolerances in dollars.
#
# The table's relative gap tolerances, 1e-4 and 1e-6, and SCIP's
# 1e-9 epsilon, in dollars on an account of $720,000. Plain Python,
# no imports.
for tol in (1e-4, 1e-6, 1e-9):
print(f"relative gap {tol:.0e} on a bound of $720,000 "
f"is ${720_000 * tol:,.4f}")
relative gap 1e-04 on a bound of $720,000 is $72.0000
relative gap 1e-06 on a bound of $720,000 is $0.7200
relative gap 1e-09 on a bound of $720,000 is $0.0007
(What the script shows, and the dual bound) The script's work is three multiplications. There is nothing to parallelize, and that is the point, since a tolerance is a decision made once per problem rather than once per node. The same arithmetic applies to a dual bound. A bound computed in floating point can lie on the wrong side of the true bound by a few units in the last place (Section 6.2). The safe-bound theorem of Section 7.3 is the remedy: evaluate the bound with every operation rounded toward \(-\infty\) and it is valid, at the cost of a hair of tightness. The C++23 listing of Section 7.3 returns, at the root of the two-variable example, \(-5.555838995037160\) from the exact dual, about \(10^{-15}\) below the exact value, and for the four children bounds between \(6 \times 10^{-7}\) and \(3 \times 10^{-6}\) below the exact child values from duals perturbed in the seventh digit, all valid.
A bound computed in floating point, and the same bound rounded down
a few units in the last place,
on the wrong side
|<--->|
-----------------------+----+-----+----------------------> value
| | |
| | computed in floating point:
| | may land here, not valid
| the true bound
every operation rounded toward -inf:
valid, at the cost of a hair of tightness
the C++23 listing of Section 7.3, on the two-variable example:
the root, from the exact dual: -5.555838995037160, about
1e-15 below the exact value
the four children, from duals perturbed in the seventh digit:
between 6e-7 and 3e-6 below the exact child values, all valid
(Exact arithmetic, for the linear case only) Exact arithmetic exists for the linear case only, and Section 5.6 gives the design and its price. Theorem 5.6.1 is the hybrid design of Cook, Koch, Steffy and Wolter: a floating-point LP solver for guidance, a safe or exact bound for every inference, and rational verification of every incumbent. SCIP 10.0 ships it as a solving mode (exact/enabled = TRUE), restricted to mixed-integer linear programs and able to write a certificate in the VIPR format that an independent checker verifies, at a slowdown that Section 5.6 measures on the MIPLIB 2017 benchmark set. Almost all of the work of an exact solver is still floating point. Exactness comes from making every inference safe and verifying every claim, not from computing everything in rationals. Nothing comparable exists for nonlinear problems. Their optima are irrational in general (Proposition 5.6.2), algebraic for polynomial data, as \(\sqrt 2\) above, and transcendental once \(\exp\), \(\log\) or \(\sin\) appear, so their certificates would need interval or algebraic arithmetic throughout the tree. For the tax problem the linear part, which is the tax itself, is in integer cents and can be made exact. The quadratic part, the risk, cannot, and the right stopping rule is a gap of a few dollars on bounds evaluated with directed rounding.
Exact MILP in SCIP 10.0 (exact/enabled = TRUE)
the hybrid design of Theorem 5.6.1 (Cook, Koch, Steffy, Wolter)
+---------------------------------------------------------------+
| floating-point LP solver ---> guidance |
| every inference ---> made on a safe or exact bound |
| every incumbent ---> verified in rational arithmetic |
+---------------------------------------------------------------+
|
| able to write
v
a certificate in the VIPR format ---> an independent checker
almost all of the work is still floating point
mixed-integer linear programs only: nonlinear optima are
irrational in general (Proposition 5.6.2), so their certificates
would need interval or algebraic arithmetic throughout the tree
Learning does not transfer
(The record has one shape) The question since 2016 has been whether a learned component can replace a hand-written one inside the tree, and Section 5.5 gives the record. The record has one shape. A graph neural network trained to imitate strong branching comes close to its decisions at a fraction of the cost, neural diving and branching reach instances with \(10^3\) to \(10^6\) variables, and learning has touched branching, node selection, cut selection, heuristics and configuration. The gains hold on instance distributions that resemble the training distribution. A learned branching rule is good on the distribution it was trained on and unpredictable off it, and the expert it imitates is itself a heuristic, which fails on problems where the relaxation's value does not move with branching. For the spatial case the one published result learns offline which branching rule to apply to a polynomial problem from instance features.B. Ghaddar, I. Gómez-Casares, J. González-Díaz, B. González-Rodríguez, B. Pateiro-López and S. Rodríguez-Ballesteros, "Learning for spatial branching: an algorithm selection approach", INFORMS Journal on Computing 35 (2023). The evidence of 2026 points away from a per-node network on the GPU and toward smaller models and learned configuration: sparse models with a few per cent of a graph network's parameters, evolved CPU-only branching rules, and learned rule selection inside Xpress Global, whose numbers Section 7.8, problem 7, records. What has shipped in a production solver with a published account is learned configuration, where a wrong prediction costs time and never correctness: CPLEX 12.10's classifier that decides whether to linearize a mixed-integer quadratic program before solving it.P. Bonami, A. Lodi and G. Zarpellon, "A classifier to decide on the linearization of mixed-integer quadratic problems in CPLEX", Operations Research 70 (2022). Gurobi's internal use of learned components is not described in a publication the author could read, and is not claimed here. The evidence standard the field has settled on is the right one for the reader to apply. A learned component is compared against a tuned default on held-out instances, in time and not in node counts, over several seeds, with the training distribution named. A decision made by a graph network cannot be explained to a compliance officer, which a tax model must allow for.
Learned configuration in CPLEX 12.10: linearize the MIQP or not
a mixed-integer quadratic program
|
v
the classifier: linearize it?
yes / \ no
v v
linearize, then solve solve it as an MIQP
\ /
v v
a correct answer either way: a wrong
prediction costs time, never correctness
(A target application, not a result) For a problem solved ten thousand times a day on accounts that look alike, the stronger forms may be exactly right, and nobody has published that they are. The tax problem of Section 9 supplies cheap labels, since the direction of every asset's trade is known after a convex solve that takes milliseconds, and a fixed structure across months. The research question is whether a predictor of the directions from lot features confines the tree to the few names whose direction is uncertain. That is a target application and not a result.
The tax problem as a learning target: the research question
accounts that look alike, a fixed structure across months:
a convex solve, in milliseconds ---> the direction of every
asset's trade: a cheap
label
|
v
lot features ----------------------> a predictor of the
directions
|
v
does it confine the tree to the few
names whose direction is uncertain?
a target application, not a result
Formulation is the user's job
(What formulation moved, in numbers) The relaxation figure of Section 2.1 showed two polygons around one set of integers, and the solver solves the polygon it is given. At the drawn angle the LP relaxation over the tighter polygon gives 5.56 against an integer optimum of 4.95, a gap of 10.9 per cent of the bound. The weaker polygon gives 5.63 and 12.0 per cent, for the same 22 integer points. Section 4 was the same lesson with curves. The hull of a two-box disjunction, the convex hull of the union of two boxes and so the tightest convex set that contains both (Section 4.2), is exact at the root, where the tightest big-M formulation gives 1.444 against a true optimum of 0.947, a gap of 52.5 per cent of that optimum, and loosening \(M\) only widens it. The perspective of a fixed charge plus a quadratic, the function \(z\,f(x/z)\) of Section 4.3 built from the cost \(f\) and the binary \(z\) that switches it on, is the convex envelope of the true cost: 1.47 against a true 1.86 at the drawn point, where the big-M relaxation charges 0.58. The cardinality figure's big-M relaxation returns the unconstrained portfolio unchanged, 11.12 per cent volatility against a true two-asset optimum of 12.45 per cent. The perspective relaxation reaches 11.71 per cent, and the factor-model split of the same covariance closes 72 to 100 per cent of the gap. On the bilinear example the lifted relaxations climb the ladder 0.40, 0.30, 0.18 with no branching at all. Reformulation as algorithm, the Bertsimas–Cory-Wright line of Section 4.9, is the same lesson taken to its conclusion. The ridge term is a modelling choice that buys convexity in the binaries, and the problem size certified went from about 400 securities, where the perspective methods had stopped, to about 3,200.D. Bertsimas, R. Cory-Wright and J. Pauphilet, "A unified approach to mixed-integer optimization problems with logical constraints", SIAM Journal on Optimization 31 (2021); D. Bertsimas and R. Cory-Wright, "A scalable algorithm for sparse portfolio selection", INFORMS Journal on Computing 34 (2022), whose Table 1 lists the largest instance each earlier method certified (Frangioni and Gentile 2009 and Zheng, Sun and Li 2014 at 400 securities) and whose Tables 7 to 9 run the S&P 500, Russell 1000 and Wilshire 5000 universes with cardinalities 10, 50, 100 and 200 under a 600-second limit on one thread (table numbers of the arXiv version, v5). The perspective numbers are those of the cardinality figure of Section 4.4 and of the factor-split computation of Section 4.5.
(Two negative results) Two negative results say why the solver cannot do this for you. Deciding whether a polynomial of degree four is convex is NP-hard, and so are strict, strong, quasi- and pseudo-convexity for even degrees of at least four. For quadratics convexity is the positive semidefiniteness of the Hessian and is polynomial.A. A. Ahmadi, A. Olshevsky, P. A. Parrilo and J. N. Tsitsiklis, "NP-hardness of deciding convexity of quartic polynomials and related problems", Mathematical Programming 137 (2013). So "convex MINLP" is a syntactic class, recognized by rules walked over the expression tree. A solver detects the convexity it has been taught to see. The rules are the tree-walk rules of Fourer and coauthors, the PSD test for quadratics, and, since SCIP 10.0.0, the detection of a rotated second-order cone in "simple bilinear constraints, e.g., x*y >= 1".R. Fourer, C. Maheshwari, A. Neumaier, D. Orban and H. Schichl, "Convexity and concavity detection in computational graphs: tree walks for convexity assessment", INFORMS Journal on Computing 22 (2010); SCIP CHANGELOG, release 10.0.0 (24 November 2025), "Features and Performance Improvements, Nonlinearity": "extended SOC detection to simple bilinear constraints, e.g., x*y >= 1", github.com/scipopt/scip/blob/master/CHANGELOG. A nonconvex term the solver cannot recognize as convex is relaxed as nonconvex, with the envelope's gap and the tree that follows. That is why the modeller declares convexity by writing the model in a recognizable form. The second negative result is representability, Section 4.10. Some sets cannot be written as mixed-integer convex sets at all, the parabola on the integers among them, and some rules of a tax model may fall outside what a MILP or a mixed-integer convex program can express.M. Lubin, J. P. Vielma and I. Zadik, "Mixed-integer convex representability", Mathematics of Operations Research 47 (2022); the midpoint test is the worked example of Section 4.10.
One constraint, x*y >= 1, as the solver's rules see it
>=
/ \
* 1
/ \
x y
the rules, walked over the expression tree
/ \
v v
not recognized as recognized: since SCIP 10.0.0,
convex a rotated second-order cone in
| a simple bilinear constraint
v |
relaxed as nonconvex: v
the envelope's gap and handled as convex
the tree that follows
the rules: the tree walks of Fourer and coauthors, the PSD test
for quadratics, and this cone detection; deciding convexity in
general is NP-hard already for polynomials of degree four
(What the solver does do) What the solver does do for the formulation is presolve and strengthening, and the ablations measure it. In Xpress Global, formula simplification in the nonlinear presolve is worth 82 solved instances on the internal test set. SCIP 8 strengthens every estimator of a convex term in a semicontinuous variable by the perspective, which reproduces Frangioni and Gentile's cuts without the modeller writing them.Belotti, Berthold, Gally, Gottwald and Pólik (2025), Section 4.5; K. Bestuzheva, A. Chmiela, B. Müller, F. Serrano, S. Vigerske and F. Wegscheider (2025), Section 2.6, and the computational study K. Bestuzheva, A. Gleixner and S. Vigerske, "A computational study of perspective cuts", Mathematical Programming Computation 15 (2023). The large choices remain the modeller's: which disjunction to write as a hull and which as a big-M, which diagonal to extract for the perspective, whether to add the regularizer that makes the binary problem convex, and whether to keep the wash sale out of the model by the trading calendar or to put it in. Formulation contributed as much as search to the progress of the last decade (Section 4). The choice of formulation is the part of the work that cannot be delegated to a solver.
Benchmarks are fragile
(A number with its date) A benchmark number is only interpretable with its date, machine, thread count, time limit and tolerances, and the MINLP benchmark of this series is the example. Section 5.3 records its three dated tables: 87 instances on an Intel i7-11700K in 2023, with Octeract solving all of them before it left GAMS and was frozen, then 200 instances on the Ryzen 9 5900X in October 2025 and again in February 2026, with BARON's count moving from 161 to 158 and SCIP's from 153 to 154 between the two runs, and a per-instance table whose LINDO column comes from a third run. None of these differences is an error. A benchmark is a set of runs on a day, and its numbers move when the runs are repeated.
Mittelmann's MINLP benchmark: three dated tables (Section 5.3)
2023 October 2025 February 2026
--+---------------------+----------------------+------------>
| | |
87 instances on an 200 instances on the the same 200 on the
Intel i7-11700K; Ryzen 9 5900X same machine
Octeract solves
all of them, then BARON 161 ------------> BARON 158
leaves GAMS and SCIP 153 ------------> SCIP 154
is frozen
and a per-instance table whose LINDO
column comes from a third run
as published: his INFORMS slides of 2023 and of 28 October 2025,
his benchmark page of 26 February 2026 and its per-instance table
(The aggregation is a choice) The aggregation is a choice, and it can reverse an order. Section 5.2 defined the shifted geometric mean and showed on a six-instance table that the shift alone decides the order of two solvers. Leaving failures out of the mean rewards a solver for giving up, and a "virtual best" row dominates every column by construction. A virtual best is quoted in this series only when it has been computed from the per-instance files and the computation is on record, as the 183 of 200 of Section 5.2 is. Section 5.1 carries the same warning across generations, in Koch, Berthold, Pedersen and Vanaret's refusal to divide a day by 104 seconds and call the ratio a speed-up.
(The vendors are absent) The commercial vendors are absent from the public tables, and the reason is on record in Section 5.2: Gurobi's retraction of 7 November 2018 of a comparison it had attributed to Mittelmann, and the withdrawals of CPLEX, Xpress, Gurobi and MindOpt that followed through December 2024. The consequence for the reader is that no independent current comparison of Gurobi, CPLEX and Xpress exists in public. Every figure of the form "2.5 times faster on MINLP than the previous version" is a vendor claim on the vendor's own test set, and Section 5.3 labels Gurobi 13's table of such figures as one.
(Four smaller fragilities) Four smaller fragilities are worth naming because each one caught a draft of this series. A table can be misread. Mittelmann's continuous nonconvex QPLIB page of 17 May 2026 lists its five solvers in one order above the table and prints the summary rows in the order of the table's header, which differs. Read in list order the page says BARON 35 and SCIP 22, read in table order it says ANTIGONE 35 and SCIP 15. The per-solver logs settle it for the table order, and they also show that the SCIP column comes from runs of April 2025 with SCIP 9.2.1 where the page names 10.0.1, so this series quotes nothing from that table.H. D. Mittelmann, "Continuous Non-Convex QPLIB Benchmark", page dated 17 May 2026, plato.asu.edu/ftp/cnconv.html, and its logs at plato.asu.edu/ftp/cnconv_logs/, read 5 October 2026: all 102 BARON, MINOTAUR and COPT logs were parsed; BARON has 33 normal completions of which the page credits 31, MINOTAUR 21 or 22 against the page's 22, COPT 28 as on the page. Hardware moves the number. The GAMS blog's update of September 2026 says that cuOpt's MIPFEAS primal integral "dropped from 0.0651 (on H100) to just 0.0383". It attributes the drop to "both software improvements and the Blackwell architecture", so the two runs did not use the same code. Its table of 16 September 2026 lists cuOpt 26.08 on the B200 at 0.0286. Mittelmann's LP feasibility table gives the GPU codes 1,000 seconds against 15,000 for the CPU codes.Bussieck and Dirkse, MIPFEAS blog, update of 16 September 2026, read 5 October 2026; Mittelmann, LPfeas page, 16 September 2026. Progress curves are self-reported. The GAMS figures of 2022 for BARON and SCIP that Section 5.1 quotes were supplied by the solver teams themselves, and the blog adds that "results of any benchmark need to be taken with a grain of salt".S. Mann, "Progress in MINLP solver technology", GAMS blog, 23 November 2022, gams.com/blog/2022/11/progress-in-minlp-solver-technology/, read 5 October 2026; the figures are quoted in Section 5.1. And the primal integral has two normalizations in use, Berthold's \(P(T) = \int_0^T \gamma(t)\,dt\) in seconds and the time average \(\bar P(T) = P(T)/T\) that MIPFEAS and the 2026 competition report, so every quoted value must say which. In either, a heuristic that never proves anything can score three times better than an exact solver that proves optimality in the same budget, which is why a primal-integral win is not a solver win.T. Berthold, "Measuring the impact of primal heuristics", Operations Research Letters 41 (2013); the two-run example of Section 7.7 (a GPU heuristic with \(\bar P(T) = 0.052\) that never reaches the optimum against an exact solver with \(0.175\) that reaches it at 180 seconds).
(What a number needs) What a benchmark number needs is therefore a list, and the rule of this series is to carry the list with the number. The list is the date of the page, the machine and its memory, the thread count, the time limit, the feasibility and gap tolerances with the gap formula, whether failures are charged at the limit, and the per-instance files. For the two measurements that do not yet exist, a GPU code on MINLPLib and a scaling curve of a global MINLP solver across devices, the protocol is written down in item 12 of the agenda of Section 7.8:
- the 200 instances by name, with the MINLPLib reference bounds on the day of the run;
- one time limit for every code, on the CPU and on the GPU alike;
- both the primal and a dual integral, with the symbol stated;
- three permutations per instance;
- deterministic and opportunistic modes as separate rows;
- every configuration run twice, with the agreement of node counts and incumbent sequences reported as a column;
- no virtual best without the computation.
The dual integral will be 1 for every heuristic that produces no bound, which is a result and not a missing value.
What is already solved
The list has a positive side, and it is worth stating at the same precision. The table collects what is solved to a certificate today, at what size, by what machinery and on whose evidence. Every row carries its date because the sizes move.Sources by row. LP: Mittelmann, LPfeas page (16 September 2026) and INFORMS 2025 slides 12 and 23; D. Applegate, M. Díaz, O. Hinder, H. Lu, M. Lubin, B. O'Donoghue and W. Schudy, "PDLP: a practical first-order method for large-scale linear programming", Mathematical Programming Computation (2026), arXiv 2501.07018 (the eleven-instance study; Gurobi's barrier solved three of them and exceeded 1 TB of memory on others). Convex QP: N. Moehle, J. Gindi, S. Boyd and M. J. Kochenderfer, "Portfolio construction as linearly constrained separable optimization", Optimization and Engineering 24 (2023), Section 7. Cardinality: Bertsimas and Cory-Wright (2022), Table 1 and Tables 7 to 9 of the arXiv version (v5). Tax: N. Moehle, M. J. Kochenderfer, S. Boyd and A. Ang, "Tax-aware portfolio construction via convex optimization", Journal of Optimization Theory and Applications 189 (2021), Section 6 and Table 1. QPLIB: H. D. Mittelmann, "Binary Non-Convex QPLIB Benchmark" (9 May 2026, AMD Ryzen 9 5900X, zero MIP gap), plato.asu.edu/ftp/qplib.html, and "Discrete Non-Convex QPLIB Benchmark (non-binary)" (12 May 2026, Intel Xeon Gold 6230), plato.asu.edu/ftp/nonbinary.html. MINLP: the page of 26 February 2026 and compare.txt, as above; the instance sizes from Mittelmann, INFORMS 2025, slide 21. MINLPLib: the library's CSV export, read 5 October 2026 (1,633 rows, the latest addition dated 18 March 2026; the library reports primal bounds at infeasibility \(10^{-8}\) and dual bounds per solver), minlplib.org; M. R. Bussieck, A. S. Drud and A. Meeraus, "MINLPLib—a collection of test models for mixed-integer nonlinear programming", INFORMS Journal on Computing 15 (2003). QUBO: H. D. Mittelmann, "Nonconvex QUBO-QPLIB Benchmark", 21 September 2026, plato.asu.edu/ftp/qubo.html.
| problem class | size and evidence | machinery |
|---|---|---|
| LP | Mittelmann LPfeas, 16 Sep 2026: 65 instances; the six huge instances of the Oct 2025 talk (up to 126 M variables, 253 M nonzeros) solved by cuPDLP-C in 1,166 to 7,288 s; PDLP (journal version) solves 8 of 11 instances with 125 M to 6.3 billion nonzeros to feasibility 1e-8 and gaps near 1% within six days | simplex, barrier; PDHG on the GPU at 1e-4 to 1e-8 with crossover |
| convex QP with a factor model | \(n = 998\) assets, \(k = 72\) factors: a convexified ADMM solve in a mean of 152 ms; the nonconvex ADMM heuristic in a mean of 251 ms, its value 0 to 10 bp from the certified bound of the convexified problem over 692 instances (Moehle, Gindi, Boyd and Kochenderfer 2023) | KKT system of size \(O(n + k)\) via the Woodbury identity; interior point |
| convex MIQP: cardinality-constrained mean-variance | about 400 securities with perspective and SDP-chosen splits (Frangioni and Gentile 2009; Zheng, Sun and Li 2014); about 3,200 (Wilshire 5000) with certificates, \(k \in \{10, 50, 100, 200\}\), 600 s, one thread (Bertsimas and Cory-Wright 2022) | perspective reformulation, diagonal split, ridge dual, outer approximation in one tree with lazy cuts |
| tax-aware MIQP at the lot level | \(n = 998\), \(k = 72\), share-level lots: CPLEX 12.9 solved 565 of 744 monthly instances within 300 s; the two-stage heuristic matched the relaxation bound, and so was certified optimal to the solver tolerance of about 0.05 bp, on 678 of 744 (Moehle, Kochenderfer, Boyd and Ang 2021, Table 1; the authors' numbers) | MIQP branch and bound; two convex solves (relaxation, then the direction-fixed QP) |
| nonconvex binary QP (QPLIB) | Mittelmann, page of 9 May 2026: 128 instances, 1 h, 12 threads (SCIP single-threaded): SHOT with Gurobi 91 solved, COPT 83, BARON 74, RAPOSa 69, SCIP 36 | RLT, SDP-derived cuts, spatial and integer branching |
| nonconvex discrete QP, non-binary (QPLIB) | Mittelmann, page of 12 May 2026: 160 instances, 3 h, 8 threads: SHOT with Gurobi 99, COPT 88, BARON 66, SCIP 40 | as above |
| general MINLP (MINLPLib selection) | Mittelmann, page of 26 Feb 2026: 200 instances up to about 100,000 variables, 2 h: BARON 158, SCIP 154, LINDO 116, SHOT 96 (quadratic instances only); at least one of the four 183, none 17 (his per-instance table compare.txt, 6 Mar 2026; Section 5.2) | spatial branch and bound over polyhedral or LP relaxations |
| MINLPLib as a whole | 1,633 instances (CSV export read 5 Oct 2026; latest addition 18 Mar 2026): 1,001 with a primal-dual gap at most 1e-6, 1,095 at most 1e-4, 538 open | the library's reported bounds, with the dual bound per solver |
| QUBO (QPLIB, unconstrained binary) | Mittelmann, page of 21 Sep 2026: 23 instances, 1 h, 12 threads: QuBowl 22 of 23 solved exactly, QUPLANE 22, BARON 13, COPT 13 | SDP bounds in branch and bound |
(Reading down the table) Read down the table, the solved part is the convex part and the small nonconvex part. Linear programs are solved at any size memory allows, on either kind of processor, to a tolerance that is a choice. Convex quadratic programs with a factor model are solved in milliseconds, and the direction-fixed subproblems of the tax problem are such programs. Convex mixed-integer quadratic programs with a few hundred to a few thousand binaries are solved with certificates when the formulation is right, and the right formulation took the field from 1996 to 2022 to find. The nonconvex rows are where two hours buy a certificate on about three quarters of a curated set, and where the hardest instances have resisted every code for years. The tax problem, as Moehle and coauthors formulated it, is in the second and fourth rows. Its only nonconvexity is the direction of each trade, its relaxation certifies the heuristic on 678 of 744 instances, and the exact MIQP solver finishes 565 of 744 in five minutes. It crosses into the nonconvex rows when the wash-sale rule enters the model as a product of a buy and a sale, and Section 9 is about both sides of that line.
The tax-aware MIQP at the lot level: two convex solves
Moehle, Kochenderfer, Boyd and Ang (2021): n = 998, k = 72,
share-level lots, 744 monthly instances
convex solve 1: the relaxation ---------------> its bound
| |
| the direction of each trade |
v |
convex solve 2: the direction-fixed QP |
| |
v v
the heuristic's value, compared with the relaxation bound:
matched it on 678 of 744, and so was certified optimal to
the solver tolerance of about 0.05 bp
the exact alternative, MIQP branch and bound in CPLEX 12.9:
565 of 744 solved within 300 s
bp: basis point, one hundredth of one per cent of account value
Where this is used
The rows of the table are finished by the solvers of Section 5.3. The LP row belongs to the simplex and barrier codes on the CPU and to the GPU first-order codes, the convex QP and convex MIQP rows to the branch and bound of CPLEX, Gurobi, Xpress and MOSEK over quadratic or conic relaxations, with SCIP adding the perspective cuts on its own, and the nonconvex rows to BARON, SCIP, COPT and SHOT with Gurobi, with SDP-based codes on the QUBO row. Moehle and coauthors' tax instances were finished by CPLEX 12.9.
What parallelizes
For the GPU the table is a map of what to accelerate. The rows that are already solved on the CPU in milliseconds gain nothing from a device unless there are thousands of them, which for accounts there are. The rows that end at a time limit gain nothing from a device unless the bound improves, which is the first item of the agenda. The sizes of the tax problem, hundreds of lots in a thousand names across possibly hundreds of thousands of accounts, are in the first category for every account whose relaxation is exact and in the second for the rest. The post that follows this series has to be built for both.