5. The engine and what moved
Draft
Sections 2 to 4 gave the mathematics of a global solver: relaxations that bound, a tree that searches, and formulations that decide how tight a relaxation can be. This section is about the solver as a piece of engineering and about how its progress is measured. It is needed now for two reasons. The first is that every leading global MINLP solver in use today is built on a mixed-integer linear programming engine. That engine, not the nonlinear machinery, is where most of the running time goes. It is also where most of the measured progress of the last twenty-five years was made. A reader who intends to design a parallel or GPU solver needs to know which components of the engine carry the weight, because accelerating the wrong component buys nothing. The second reason is that the performance numbers in this series come from benchmarks, and benchmarks have conventions. A solved count, a shifted geometric mean, a virtual best and a vendor's "twice as fast" are four different kinds of statement. This section says what each of them means and what it does not.
The six subsections go from measurement to machinery. Section 5.1 reports how much faster the MILP engine became between 2001 and 2020, and which components the ablation studies credit. An ablation is an experiment in which one component of a fixed solver is switched off and the slowdown is measured. Section 5.2 defines the benchmark arithmetic the field uses, with a live figure. Section 5.3 names the global solvers that return certificates, gives their lineage, and prints the dated scoreboard. Section 5.4 covers the primal heuristics that need no linear program and the alternating-direction heuristics for separable nonconvex problems. It states the Shapley–Folkman theorem and proves the Udell–Boyd bound built on it, which explains why the relaxations of those problems are tight. Section 5.5 asks what machine learning has changed inside the tree and what has shipped. Section 5.6 explains what a floating-point solver's answer means, what an exact solver does differently, and why exactness stops at the linear case.
A thousandfold
This subsection establishes two numbers and one order: how much faster the MILP engine became between 2001 and 2020 with the hardware held fixed, how many instances became solvable at all, and which components of the engine the ablation studies credit. The GPU design of Section 7 competes against this engine as it stands, so the order among the components matters more than the factor. The cleanest measurement of progress in MILP solving is the study of Koch, Berthold, Pedersen and Vanaret. They ran the solvers of 2001 and the solvers of 2020 on the same instances and separated the contribution of hardware from the contribution of algorithms.T. Koch, T. Berthold, J. Pedersen and C. Vanaret, "Progress in mathematical programming solvers from 2001 to 2020", EURO Journal on Computational Optimization 10 (2022), 100031; arXiv 2206.09787 (2022). Every number attributed to this paper below is taken from its text. Their design rests on one notion that recurs in every benchmark table of this section.
Definition 5.1.1 (virtual best, virtual worst). Given solvers \(S_1, \dots, S_k\) run on the same instances under the same limits, the virtual best is the fictitious solver whose time on each instance is the smallest of the \(k\) times. It counts an instance as solved if any of the \(k\) solvers solved it. The virtual worst takes the largest time on each instance and counts an instance as solved only if all \(k\) solvers solved it.
(The factor, and how it was measured) Koch et al. compared the virtual best of three 2001 codes (CPLEX 7.0, Xpress-MP 14.10 and MOSEK 3.2.1.8) with the virtual best of five 2020 codes (CPLEX 12.10, Gurobi 9.0, Xpress 8.11, MOSEK 8.1 and COPT 1.4). In the terms of the definition, \(k = 3\) on the 2001 side and \(k = 5\) on the 2020 side: an instance's 2001 time is the smallest of the three codes' times on it, it counts as solved in 2001 if any of the three solved it, and likewise for 2020 with five. Comparing virtual bests rather than two single codes makes the measured progress a property of the state of the art in each year, not of one vendor's choices. The 2001 machine was an 870 MHz Pentium III. The 2020 machine was an Intel Core i7-9700K at 3.6 GHz with 64 GB of memory. The measured hardware factor varied between about 7 and more than 45 depending on the load on the modern machine, and the authors adopt "about 20". With the hardware held fixed, the algorithms alone improved by a factor of about 9 for LP and about 50 for MILP. The combined speed-ups are about 180 and about 1,000. The MILP factor is computed on the 105 of 339 instances that both generations solved within the 24-hour limit, because a speed-up can only be measured where both codes finish. Spread over the nineteen years between the two releases, a factor of 50 is a compound rate of 22.9% a year. Over twenty years it is 21.6%. That is the "22% every year" of the paper.
(Counts are not speed-ups) The second fact is about what the old codes could not do at all. Of the 240 instances of the MIPLIB 2017 benchmark set, 149 could not be solved by any of the 2001 codes within a day, even on the modern machine. The 2020 codes solve those 149 in a geometric mean of 104 seconds, the benchmark average that Section 5.2 defines.Koch et al. (2022), Conclusion. The paper's Section 3.4 gives the same share as 156 of 240; the paper calls both 62% and does not reconcile them. This post quotes the Conclusion's figure. The authors warn against dividing one number by the other. A ratio "1 day / 104 sec = 830-fold" would be a speed-up measured on instances where one side never finished. Their Section 3.4 shows that the choice of time limit and of penalty for a time-out moves a single aggregate ratio on the same data from 37 to 8,125. The reason is that an instance which timed out has no running time, only a limit, and a ratio that includes it has to invent one: the limit itself, some multiple of it, or a penalty of the analyst's choosing. With 149 of 240 instances unsolved on one side, the invented numbers dominate the aggregate, and the choice among them is in effect the choice of the answer. A speed-up on the instances both codes solve and a count of newly solvable instances are two separate facts. This series reports them separately throughout.
(What the paper credits) What moved is listed in the same paper. On the LP side the gains came from refined implementations, 64-bit addressing and parallel barrier codes, and the authors note that "little progress has been made on the theoretical side". On the MILP side they list new heuristics such as RINS and local branching, and new or improved cutting planes such as multi-commodity-flow cuts. They add "a large number of tricks for the bag": conflict analysis, symmetry detection, solution polishing and dynamic search, together with parallel tree search. RINS and local branching are the large-neighbourhood searches of Section 3.2. Conflict analysis, which turns an infeasible or pruned node into an inequality valid for the whole tree, is in Section 3.6. Symmetry detection finds variables the problem does not distinguish, so that of several branches that differ only by a renaming one is explored; solution polishing is a search for better feasible points near the incumbent, run late in the solve; and dynamic search is the name under which CPLEX ships its default tree search. They add that the features overlap, so that a two-point comparison cannot assign a factor to each. Per-component factors come from ablation studies.
(Three ablation studies, two decades) Three ablation studies span two decades. Bixby, Fenelon, Gu, Rothberg and Wunderling measured CPLEX 8.0 against CPLEX 5.0 on 758 instances: a geometric-mean speed-up of about 12, and about 528 on the instances that took the old code more than 27 hours.R. E. Bixby, M. Fenelon, Z. Gu, E. Rothberg and R. Wunderling, "Mixed-integer programming: a progress report", in M. Grötschel (ed.), The Sharpest Cut (MPS-SIAM, 2004), 309–325; the 758, 12 and 528 are as summarized in Koch et al. (2022), Section 1.1. Their ablation table for CPLEX 8.0 is the one most often quoted. The usual summary of it is that cutting planes were by far the most valuable component, then presolve, then the choice of branching variable, then the primal heuristics, with factors of about 54, 11, 3 and 1.4.This order and these four factors are recalled from the table as reproduced in R. Bixby and E. Rothberg, "Progress in computational mixed integer programming: a look back from the other side of the tipping point", Annals of Operations Research 149 (2007), 37–41, and in A. Lodi, "Mixed integer programming computation", in 50 Years of Integer Programming 1958–2008 (Springer, 2010), 619–645. Neither of those two papers nor the chapter itself was re-read for this post, and no open copy of the chapter exists (Unpaywall, 5 October 2026). The four numbers stay unverified until they are copied from one of the three. Achterberg and Wunderling repeated the exercise for CPLEX 12.5 on 2,928 models. They found a geometric-mean speed-up of 4.71 over version 8.0, up to 78.6 on the instances that were hard for 8.0, and 732 time-outs at 10,000 seconds for the old code against 87 for the new.T. Achterberg and R. Wunderling, "Mixed integer programming: analyzing 12 years of progress", in M. Jünger and G. Reinelt (eds), Facets of Combinatorial Optimization (Springer, 2013), 449–481; the figures are as summarized in Koch et al. (2022), Section 1.1. The chapter's own ablation table is paywalled and was not re-read for this post. Their component ablation is reported in the same qualitative terms, with conflict analysis and symmetry handling among the smaller factors. Lodi's survey covers the longer arc. CPLEX 1.2 of 1991 against CPLEX 11.0 of 2007 on 1,734 instances gives a geometric-mean speed-up of 67.9. On 1,852 instances with a 30,000-second limit the old code solved 15.0% and the new 67.1%.Lodi (2010), as summarized in Koch et al. (2022), Section 1.1. For the LP side see R. E. Bixby, "Solving real-world linear programs: a decade and more of progress", Operations Research 50 (2002), 3–15, and R. E. Bixby, "A brief history of linear and mixed-integer programming computation", Documenta Mathematica, Extra Volume "Optimization Stories" (2012), 107–121, which puts the machine-independent LP improvement from 1988 to 2004 at a factor of 3,300 and the machine improvement at 1,600. The one ablation with current numbers that this post could verify is Mexi's, for SCIP 10. Switching presolving off costs a factor of 2.60 in shifted geometric mean time, the benchmark average of Section 5.2, and replacing the branching rule by a random one 2.55. Switching cutting planes off costs 2.10, and switching off primal heuristics or conflict analysis 1.31 each, on 349 instances with five seeds.G. Mexi, The Two Faces of Mixed-Integer Programming: Primal and Dual Progress, doctoral thesis, Technische Universität Berlin (2026), Section 2.7.1 and Table 2.1; read 5 October 2026 from the university's repository depositonce.tu-berlin.de. The test set is 349 instances from the MIPLIB 2003, 2010 and 2017 benchmark sets and the COR@L collection that SCIP 10.0 or an earlier release can solve, run with five seeds (1,745 instance-seed pairs) on Intel Xeon Gold 6338 machines with a two-hour limit. The thesis describes its findings as "broadly consistent with the results of Achterberg and Wunderling (2013) for CPLEX 12.5". Two remarks keep the three studies in proportion. "Random branching" is not the same ablation as switching a component off, since a solver cannot run without some branching rule. And the studies agree on what is large and what is small, cuts, presolve and branching on one side and heuristics and conflict analysis on the other, but not on the order among the three large factors. The CPLEX 8.0 table as recalled puts cuts first by a wide margin. The SCIP 10 table puts presolve first and branching second, with cuts third.
(The order of components) The lesson for a GPU design is direct. The components that carry the engine are the cut loop and presolve, which are mostly sequential passes over rows and columns, and branching, which is a batch of child LPs. The heuristics, which are the easiest to parallelize, carry the least. Section 7 returns to this order.
(The nonconvex curve, self-reported) Nonconvex MINLP has its own, shorter curve, measured by the solver teams themselves. In November 2022 GAMS published both measurements. For BARON the statement is: "Over the last approximately 20 years, BARON's performance on a set of 87 MINLPLib instances has improved by about 10x in speed, and 3x in the number of problems solvable." For SCIP the code was recompiled at every version from 3.0.2 (October 2013) to 8.0.1 (June 2022) with the same CPLEX and Ipopt underneath, so that only SCIP's own code changed. On 183 MINLPLib instances "the number of solvable models increased from 54 with SCIP 3.0.2, to 108 with SCIP 8.0.1, and at the same time the mean solve time was reduced 3.1-fold".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 BARON figure is Sahinidis's own and the SCIP figure is from the SCIP team's anniversary talk; both are self-reported. The blog adds that "results of any benchmark need to be taken with a grain of salt". The SCIP 8 paper states the same measurement as twice as many instances solved and a speed-up of about three.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), 287–310; arXiv 2301.00587 (2023), Section 4. A factor of 10 over nineteen years is 12.9% a year and a factor of 3.1 over nine years is 13.4% a year, against the MILP engine's 22%. The block below the next figure computes these rates and prints the numbers quoted in this subsection.
# Compound annual rates behind the progress figures of this subsection.
#
# A speed-up by a factor F over Y years is a compound rate of
# F^(1/Y) - 1 a year. The factors are the ones quoted in the text:
# Koch et al. (2022) for MILP and LP, the GAMS blog for BARON and SCIP.
FIGURES = [
("MILP algorithms 2001 to 2020, factor 50, over 19 years", 50, 19),
("the same factor over 20 years", 50, 20),
("LP algorithms, factor 9, over 19 years", 9, 19),
("BARON 2003 to 2022, factor 10, over 19 years", 10, 19),
("SCIP 3.0.2 to 8.0.1, factor 3.1, over 9 years", 3.1, 9),
]
for name, factor, years in FIGURES:
rate = factor ** (1 / years)
root = f"{factor}^(1/{years})"
print(f"{name}:")
print(f" {root:<9} = {rate:.3f}, "
f"i.e. {100 * (rate - 1):.1f}% a year")
MILP algorithms 2001 to 2020, factor 50, over 19 years:
50^(1/19) = 1.229, i.e. 22.9% a year
the same factor over 20 years:
50^(1/20) = 1.216, i.e. 21.6% a year
LP algorithms, factor 9, over 19 years:
9^(1/19) = 1.123, i.e. 12.3% a year
BARON 2003 to 2022, factor 10, over 19 years:
10^(1/19) = 1.129, i.e. 12.9% a year
SCIP 3.0.2 to 8.0.1, factor 3.1, over 9 years:
3.1^(1/9) = 1.134, i.e. 13.4% a year
(The engine a GPU design competes against) The block computes five scalar roots, and nothing in it is parallel. Both factors, 9 for LP and 50 for MILP, were measured with the 2001 codes single-threaded, apart from two barrier LP codes, and the 2020 codes on one thread or eight, whichever was faster on each instance. The algorithms that gave the factor of 50, the dual simplex with warm starts (Section 3.1), the cut loop (Section 3.3), presolve (Section 2.5) and reliability branching (Section 3.1), are the ones Section 7 finds hardest to move to a GPU. A GPU design for nonconvex MINLP competes against this engine as it stands in 2026, not against the codes of 2001, and the ablation studies say which components decide that competition.
How the field measures itself
Every number in Section 5.3 comes from a benchmark, and a benchmark is a set of conventions: which instances, which time limit, which machine, how the times are averaged, how a failure is scored, and whose results are published. This subsection fixes those conventions, because the same raw times can be made to rank the same solvers in different orders. It names the three instance libraries. It defines the shifted geometric mean and shows in a figure what the shift does. It computes the virtual best on the current MINLP benchmark and defines the trusted gap. It then records three facts that make a benchmark number fragile: performance variability, the vendors' withdrawal from public benchmarks, and the tuning of solvers to known test sets. It ends with the rules this series follows whenever it quotes one.
(The three instance libraries) MINLPLib is the library behind every MINLP number in this series. Bussieck, Drud and Meeraus assembled it in 2003 as a collection of GAMS models, and it is maintained today at minlplib.org. Each instance has a page recording the best known point and its infeasibility (points are accepted at an infeasibility of at most \(10^{-8}\)), the dual bounds reported by individual solvers, and the problem's classification.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), 114–119; the library is at minlplib.org. The library's instance export of 5 October 2026 lists 1,633 instances, the most recent added on 18 March 2026. Of these, 376 are classified convex, 1,218 nonconvex and 39 unclassified, and 1,058 have discrete variables. The largest type classes are mixed-binary NLP (423), mixed-binary QCP (275), NLP (239), QCP (217), general MINLP (108), binary QP (83) and QP (74). By the bounds recorded in the export, 1,001 instances have a relative gap between the best known point and the recorded dual bound of at most \(10^{-6}\), and 1,095 a gap of at most \(10^{-4}\). The remaining 538 are open in the library's own sense, with a gap above \(10^{-4}\) against the dual bound the library records. That recorded bound is the library's best bound, not the trusted bound of Definition 5.2.3 below.MINLPLib instance export (CSV) from minlplib.org, read 5 October 2026. The counts depend on the export date and should be re-read with it. QPLIB is the quadratic counterpart. Furini and twelve co-authors selected 453 instances, 319 discrete and 134 continuous, from 8,164 submissions, keeping those that were hard for the solvers of the time and diverse in structure. Each instance carries a three-letter code. The first letter classifies the objective (among others, L linear, C convex quadratic, Q nonconvex quadratic). The second classifies the variables (among others, C continuous, B binary, M mixed, I integer), and the third the constraints (among others, N none, B box, L linear, C convex quadratic, Q nonconvex quadratic). A QBL instance therefore has a nonconvex quadratic objective over binary variables with linear constraints.F. Furini, E. Traversi, P. Belotti, A. Frangioni, A. Gleixner, N. Gould, L. Liberti, A. Lodi, R. Misener, H. Mittelmann, N. V. Sahinidis, S. Vigerske and A. Wiegele, "QPLIB: a library of quadratic programming instances", Mathematical Programming Computation 11 (2019), 237–265; the library page qplib.zib.de states "Starting from 8,164 submitted instances, the final version of QPLIB contains 319 discrete and 134 continuous instances" (read 5 October 2026). The full alphabet of the code is in the paper's Section 3 and was not re-read for this post. MIPLIB 2017 is the MILP library: 1,065 instances in the collection and a benchmark set of 240. The benchmark set was chosen by solving a mixed-integer program over the candidates that asked for diversity across feature clusters and for solvability by at least one solver within the limit.A. Gleixner, G. Hendel, G. Gamrath, T. Achterberg, M. Bastubbe, T. Berthold, P. Christophel, K. Jarck, T. Koch, J. Linderoth, M. Lübbecke, H. D. Mittelmann, D. Ozyurt, T. K. Ralphs, D. Salvagnin and Y. Shinano, "MIPLIB 2017: data-driven compilation of the 6th mixed-integer programming library", Mathematical Programming Computation 13 (2021), 443–490; miplib.zib.de. Mittelmann runs the 240 after they have been "preprocessed by using SCIP to randomly perturb rows/columns and performing a quick presolve", a defence against solvers tuned to the public files.H. D. Mittelmann, "The MIPLIB2017 benchmark instances (preprocessed data)", run of 7 July 2026, plato.asu.edu/ftp/milp.html, read 5 October 2026.
Run times are aggregated by one statistic, which the next definition fixes.
Definition 5.2.1 (shifted geometric mean). For run times \(t_1, \dots, t_n > 0\) and a shift \(s \ge 0\),
\[\operatorname{sgm}_s(t) \;=\; \exp\Big( \frac{1}{n} \sum_{i=1}^n \ln \max\big(1,\ t_i + s\big) \Big) \;-\; s .\]Mittelmann's tables use \(s = 10\) seconds and count a run that fails or reaches the time limit as having taken the time limit. In a "scaled" row they divide each solver's mean by the smallest mean in the table, so that the best solver scores 1.H. D. Mittelmann, "Shifted geometric mean", plato.asu.edu/ftp/shgeom.html; the guard \(\max(1, \cdot)\) is the convention of T. Achterberg, Constraint Integer Programming, PhD thesis, TU Berlin (2007), Appendix A. The SCIP 10 report uses \(s = 1\) s for times and \(s = 100\) for node counts in its exact-mode study (C. Hojny et al., "The SCIP Optimization Suite 10.0", arXiv 2511.18580 (2025), Section 3.1.9).
(The formula on three times) On the three times 1, 10 and 100 s, which are solver A's first three in the toy benchmark drawn below, the formula with \(s = 10\) reads \(\exp\big(\tfrac13(\ln 11 + \ln 20 + \ln 110)\big) - 10 = (11 \cdot 20 \cdot 110)^{1/3} - 10 = 18.9\) s. That sits between the geometric mean of the three, \(10.0\) s, and their arithmetic mean, \(37.0\) s, which is the proposition below in one instance; the block after the figure prints all three numbers. The guard \(\max(1, \cdot)\) only acts when \(t_i + s < 1\), so it is irrelevant for \(s \ge 1\), and the figure below omits it. The shift has a precise effect, which is worth stating as a proposition because the figure is built on it.
Proposition 5.2.2 (what the shift does). Let \(\tau_i = \max(1, t_i + s)\). (a) \(\min_i \tau_i - s \le \operatorname{sgm}_s(t) \le \max_i \tau_i - s\), and \(\operatorname{sgm}_s(t)\) is nondecreasing in every \(t_i\). (b) If all \(t_i \ge 1\), then \(\operatorname{sgm}_0(t)\) is the geometric mean of the \(t_i\). (c) For fixed \(t\) with all \(t_i \ge 1\), \(\operatorname{sgm}_s(t)\) is nondecreasing in \(s\) and converges to the arithmetic mean \(\frac{1}{n} \sum_i t_i\) as \(s \to \infty\). The shifted geometric mean therefore interpolates between the geometric mean, which weights every instance by its relative speed-up, and the arithmetic mean, which weights every instance by its absolute time.
Proof. (a) The geometric mean of positive numbers lies between their minimum and their maximum, and each \(\tau_i\) is nondecreasing in \(t_i\). (b) With \(s = 0\) and \(t_i \ge 1\) every \(\tau_i = t_i\). (c) When all \(t_i \ge 1\) the guard is inactive and \(\operatorname{sgm}_s(t) = G(s) - s\) with \(G(s) = \prod_i (t_i + s)^{1/n}\) the geometric mean of the shifted times. Differentiating, \(G'(s) = G(s) \cdot \frac{1}{n} \sum_i (t_i + s)^{-1} = G(s)/H(s)\), where \(H(s)\) is the harmonic mean of the shifted times. Since the geometric mean is at least the harmonic mean, \(G'(s) \ge 1\) and \(\frac{d}{ds}\operatorname{sgm}_s(t) = G'(s) - 1 \ge 0\). For the limit, \(\operatorname{sgm}_s(t) = s\big[\prod_i (1 + t_i/s)^{1/n} - 1\big]\) and \(\ln(1 + t_i/s) = t_i/s + O(s^{-2})\), so \(\operatorname{sgm}_s(t) = s\big[\exp\big(\frac{1}{n}\sum_i t_i/s + O(s^{-2})\big) - 1\big] = \frac{1}{n}\sum_i t_i + O(s^{-1})\). ∎
(Three consequences for reading a table) Three consequences govern how a table is read. A shift of 10 seconds makes an instance solved in 0.1 s and one solved in 5 s nearly indistinguishable. This is intended: both are easy, and a geometric mean without the shift would let a solver that is a hundred times faster on trivial instances dominate the ranking. The second consequence is that the ratio of two solvers' shifted means is not invariant under rescaling all times, since a faster machine changes the ratios unless the shift is rescaled with them. The shift is a number of seconds that does not scale with the times: doubling every time doubles a geometric mean exactly, but a shifted mean by less, and by a different amount for each solver, since the shift matters more to a solver whose times are short. This is why the "scaled" row divides by the best mean on the same machine, and why numbers from different machines are never divided by one another. The third consequence is not about the shift but about the failures. A time-out is censored data: the true time is known only to exceed the limit, and whatever number stands in for it is a convention. The mean of a column that charges failures at the limit depends on the limit, and raising the limit can only raise the same solver's mean. A newly solved instance is charged its true time, which exceeds the old limit, and an instance still unsolved is charged the new, larger limit. Every mean therefore comes with its solved count.
(The toy benchmark, and what to look for) The figure makes these statements concrete on a table small enough to read in full. Six instances are solved by three solvers, A, B and C. Each solve time is a dot on a logarithmic time axis, the fastest solver on each row carries a ring, and a run that reached the two-hour limit sits on the dashed line as a red cross. The three colours are only labels for the three solvers here. Red keeps its meaning as a failure. Below the table, each solver's geometric mean, shifted geometric mean and arithmetic mean are drawn on the same axis, with the virtual best underneath. Proposition 5.2.2 is visible as geometry. The filled dot of the shifted mean sits between the hollow dot of the geometric mean and the bar of the arithmetic mean. It slides from one toward the other as the slider moves \(s\) from 0 to 100 seconds. The last control decides how a failure enters the mean. The thing to look for is that the same six columns of times can rank the solvers in different orders, and that the shift decides which order the table reports.
(The three scenarios, in numbers) In the default table solver A takes 1, 10, 100, 4 and 50 seconds on the first five instances and fails on the sixth. B takes 5 seconds on each of the first five and 300 on the sixth. C takes 2, 10, 20, 30, 50 and 80 seconds. With the failure charged at 7,200 seconds the arithmetic means are 1,228, 54.2 and 32.0 s, which rank the solvers C, B, A. The geometric means are 33.6, 9.9 and 19.1 s, which rank them B, C, A. The shifted geometric means at \(s = 10\) are 62.6, 14.8 and 24.0 s, the geometric order again. Moving the slider shows where the two orders meet. B and C change places at \(s = 80\) s, and at \(s = 100\) the shifted means 151.5, 31.2 and 29.5 s reproduce the arithmetic order. B and C solve six of six, A five. The virtual best takes 1, 5, 5, 4, 5 and 80 seconds, has means 16.7, 5.8 and 9.0 s, and leads every column. The three solvers' shifted means are 7.0, 1.7 and 2.7 times its own. That ratio is the quantity Mittelmann's "scaled" row reports, taken there against the best real solver rather than the virtual best. The second scenario gives A 120 seconds on the sixth instance, so that everyone finishes. The three means then give three different orders. The arithmetic mean ranks the solvers C, A, B (32.0, 47.5, 54.2 s). The geometric mean ranks them B, A, C (9.9, 17.0, 19.1 s). The shifted mean at \(s = 10\) ranks them B, C, A (14.8, 24.0, 27.2 s), with A and C changing places at \(s = 2.2\) s and B and C at 80 s. The third scenario lets A and B fail one instance each. They are then nearly tied, 62.6 against 59.5 s, and change places at \(s = 21\) s, while C leads every mean and solves the most. Switching the last control to "not at all" removes a failure from its solver's mean. A's shifted mean then falls from 62.6 to 18.9 s although A has solved nothing more. In the third scenario the order by the shifted mean becomes B, A, C while C is the solver that finishes the most instances. That is the case against the convention, and it is why Mittelmann charges a failure at the time limit.
The default table as s grows (a failure charged at 7,200 s)
s = 0 s = 10 s = 100 s -> inf
geometric shifted shifted arithmetic
A 33.6 ---> 62.6 ---> 151.5 ---> 1,228
B 9.9 ---> 14.8 ---> 31.2 ---> 54.2
C 19.1 ---> 24.0 ---> 29.5 ---> 32.0
virtual best 5.8 ---> 9.0 ---> 13.9 ---> 16.7
order B, C, A B, C, A C, B, A C, B, A
^
B and C change places
at s = 80
every row climbs from its geometric mean toward its arithmetic
mean (Proposition 5.2.2); the shift decides which order is shown
The block below recomputes the default table at two shifts, the crossing point, and the two worked examples that the figure's readout also prints. The first example takes the first three instances alone. Solver A's first three times, 1, 10 and 100 s, have arithmetic mean 37.0, geometric mean 10.0 and shifted mean 18.9 s, while B's 5, 5, 5 give 5.0 three times. Had A failed the third and been charged at the limit, its shifted mean would be \((11 \cdot 20 \cdot 7210)^{1/3} - 10 = 106.6\) s. Every number the text quotes appears in the output.
# Scoring a benchmark table.
#
# The toy table of the figure (six instances, three solvers) under the
# arithmetic, geometric and shifted geometric means, with a failure
# charged at the two-hour limit; then the shift at which B and C change
# places, and the two worked examples on the first three instances.
import numpy as np
LIMIT = 7200.0
table = {
"A": [1, 10, 100, 4, 50, None],
"B": [5, 5, 5, 5, 5, 300],
"C": [2, 10, 20, 30, 50, 80],
}
def charged(times):
"""A failure (None) enters the mean at the time limit."""
return np.array([LIMIT if t is None else t for t in times], float)
def sgm(t, s):
"""(prod (t_i + s))^(1/n) - s, computed in log space."""
return np.exp(np.mean(np.log(t + s))) - s
def score(times, s):
"""Arithmetic, geometric and shifted means, and the solved count."""
t = charged(times)
return dict(am=t.mean(), gm=sgm(t, 0.0), sgm=sgm(t, s),
solved=sum(x is not None for x in times))
def crossover(ta, tb):
"""The shift at which two solvers' shifted means agree, by bisection."""
d = lambda s: score(ta, s)["sgm"] - score(tb, s)["sgm"]
grid = [0.0] + [10 ** (k / 20) for k in range(-20, 81)]
for lo, hi in zip(grid[:-1], grid[1:]):
if np.sign(d(lo)) != np.sign(d(hi)):
for _ in range(60):
mid = 0.5 * (lo + hi)
if np.sign(d(mid)) == np.sign(d(lo)):
lo = mid
else:
hi = mid
return 0.5 * (lo + hi)
# The virtual best: the smallest time on each instance (Definition 5.1.1).
vb = [min(t for t in col if t is not None) for col in zip(*table.values())]
print("virtual best, instance by instance:", vb)
for s in (10, 100):
rows = ({k: score(v, s) for k, v in table.items()}
| {"virtual best": score(vb, s)})
print(f"\nshift s = {s} s: solver arithmetic geometric"
f" shifted solved")
for k, q in rows.items():
print(f" {k:<14}{q['am']:9.1f}{q['gm']:11.1f}{q['sgm']:9.1f}"
f" {q['solved']} of 6")
def order(key):
return ", ".join(sorted(table, key=lambda k: rows[k][key]))
print(f" order by arithmetic {order('am')}; "
f"by geometric {order('gm')}; by shifted {order('sgm')}")
vb_sgm = rows["virtual best"]["sgm"]
print(" ratio of shifted means to the virtual best:",
", ".join(f"{k} {rows[k]['sgm'] / vb_sgm:.1f}x" for k in table))
print(f"\nB and C change places at s = "
f"{crossover(table['B'], table['C']):.1f} s")
print("\nThe first three instances alone, s = 10:")
A3, B3 = np.array([1.0, 10.0, 100.0]), np.array([5.0, 5.0, 5.0])
print(f" A: arithmetic {A3.mean():.1f}, geometric {sgm(A3, 0):.1f}, "
f"shifted {sgm(A3, 10):.1f};"
f" B: {B3.mean():.1f}, {sgm(B3, 0):.1f}, {sgm(B3, 10):.1f}")
print(" A failing the third, charged at the limit:")
print(" (11 * 20 * 7210)^(1/3) - 10 =",
f"{sgm(np.array([1.0, 10.0, LIMIT]), 10):.1f}")
virtual best, instance by instance: [1, 5, 5, 4, 5, 80]
shift s = 10 s: solver arithmetic geometric shifted solved
A 1227.5 33.6 62.6 5 of 6
B 54.2 9.9 14.8 6 of 6
C 32.0 19.1 24.0 6 of 6
virtual best 16.7 5.8 9.0 6 of 6
order by arithmetic C, B, A; by geometric B, C, A; by shifted B, C, A
ratio of shifted means to the virtual best: A 7.0x, B 1.7x, C 2.7x
shift s = 100 s: solver arithmetic geometric shifted solved
A 1227.5 33.6 151.5 5 of 6
B 54.2 9.9 31.2 6 of 6
C 32.0 19.1 29.5 6 of 6
virtual best 16.7 5.8 13.9 6 of 6
order by arithmetic C, B, A; by geometric B, C, A; by shifted C, B, A
ratio of shifted means to the virtual best: A 10.9x, B 2.2x, C 2.1x
B and C change places at s = 80.3 s
The first three instances alone, s = 10:
A: arithmetic 37.0, geometric 10.0, shifted 18.9; B: 5.0, 5.0, 5.0
A failing the third, charged at the limit:
(11 * 20 * 7210)^(1/3) - 10 = 106.6
The shifted mean step by step: solver A, first three instances
A solved all three, s = 10
t_i 1 10 100
t_i + s 11 20 100 + 10
sgm_10 = (11 * 20 * (100 + 10))^(1/3) - 10 = 18.9
A fails the third, charged at the limit of 7,200 s
t_i 1 10 7200
t_i + s 11 20 7210
sgm_10 = (11 * 20 * 7210)^(1/3) - 10 = 106.6
A on the three it solved:
geometric 10.0 <= shifted 18.9 <= arithmetic 37.0
B (5, 5, 5): 5.0 under all three means
The cost of the computation is one logarithm per entry, and the crossing is sixty bisection steps on a one-dimensional function. Nothing here needs parallelism. The block's purpose is that the reader can change a time and watch the ranking move.
(The virtual best of the 2026 MINLP benchmark) The virtual best of Definition 5.1.1 is the one benchmark number that Mittelmann's pages never print and that is quoted anyway. For the MINLP benchmark of 26 February 2026 it can be computed from his per-instance table. Count an instance as solved when its printed time is below 7,200 seconds and is not marked as a time-out. This rule reproduces the four published counts (BARON 158, SCIP 154, SHOT 96, and 119 for LINDO in the table's later run) and gives 183 of 200 instances, 91.5%, solved by at least one of the four published solvers. The table's own row "optimal auto settings" prints the same 183 with 17 time-outs. The number is therefore Mittelmann's own, printed in the table rather than on the page.H. D. Mittelmann, per-instance comparison table for the MINLP benchmark, plato.asu.edu/ftp/compare.txt, file dated 6 March 2026, read and recomputed on 5 October 2026. Under the stricter reading that excludes the two rows flagged "#" (an infeasible instance that every solver proved infeasible, and one instance SCIP proved infeasible while BARON timed out) the union is 181. The 17 instances no published solver solves are bayes2_20, bayes2_30, camshape100, casctanks, chp_shorttermplan2b, faclay80, feedtray, gasnet, gasprod_sarawak16, gastrans582_mild11, multiplants_mtg1a, parabol5_2_3, supplychainp1_022020, topopt-mbb_60x40_50, transswitch0009r, wastepaper6 and water4. Without LINDO the union is 180. Forty-nine instances are solved by all four, and 21 by exactly one (BARON 10, SCIP 6, LINDO 3, SHOT 2). Two caveats attach. Gurobi and Xpress were run but their results are suppressed, so this is a virtual best over four solvers, not six. And the table's LINDO column (119 solved, from a run of 6 March 2026 under GAMS 53.1) is a later run than the page's (116 solved), so even one published page mixes runs. The number this series uses is "183 of 200", with the date and the method, and never "92%" alone.
The virtual best over four solvers of Mittelmann's MINLP benchmark
200 instances; Mittelmann's page of 26 February 2026 and his
per-instance table of 6 March 2026; solved = printed time below
7,200 s, not marked as a time-out; LINDO from the table's later
run (the page has 116); Gurobi and Xpress were run, but their
results are suppressed: a union over four solvers, not six
BARON 158 SCIP 154 LINDO 119 SHOT 96
| | | |
'--------------+------+-------+-------------'
|
| union, instance by instance
v
solved by at least one of the four 183 of 200 (91.5%),
as the table's own row
"optimal auto settings"
prints it
by all four 49
by exactly one 21: BARON 10, SCIP 6,
LINDO 3, SHOT 2
solved by none of the four 17
the union without LINDO: 180
the union without the two rows flagged "#": 181
Where the solvers do not agree on the optimum, the bound itself needs a convention.
Definition 5.2.3 (trusted dual bound, trusted gap). For a minimization instance with best known value \(z_{\mathrm{inc}}\) and dual bounds \(d_1, \dots, d_k\) reported by \(k\) independent solvers, the trusted dual bound is \(d_{\mathrm{tr}} = \min_j d_j\), the one bound that every solver's report implies. The trusted gap is \((z_{\mathrm{inc}} - d_{\mathrm{tr}}) / |z_{\mathrm{inc}}|\). An instance is called open at tolerance \(\varepsilon\) when its trusted gap exceeds \(\varepsilon\).
Each solver's bound is a claim that the optimum is at least \(d_j\). A solver claiming \(d_3 = 99.2\) also claims every smaller bound, so the minimum is the claim all of them share.
Example 5.2.4 (the trusted gap). With best known value 100 and dual bounds 98.0, 98.5 and 99.2 from three solvers, the trusted dual bound is 98.0 and the trusted gap is 2%. The strongest single claim of 0.8% is not trusted until a second solver reaches it. Requiring agreement of two solvers out of three would give 98.5 and a gap of 1.5%.
This series uses the unanimous convention, and it is the definition under which it calls an instance open. MINLPLib's instance pages list the dual bounds per solver, which is what makes the test possible. A bound that a known feasible point contradicts is evidence of a bug in the solver that reported it. A bound that only one solver has reached is a conjecture until a second solver confirms it.
(A column is a software stack) A benchmark column is a software stack, not an algorithm. BARON solves its LP, MIP and QP relaxations with CPLEX, Xpress, CLP/CBC or HSL LA04, chosen automatically by default, and its NLP subproblems with MINOS, SNOPT, IPOPT, FilterSD or FilterSQP. SCIP under GAMS runs with CPLEX or SoPlex as the LP solver and Ipopt as the NLP solver. SHOT is an outer-approximation scheme, the convex-case method of Section 3.4, whose work is done by Gurobi, CPLEX or Cbc. In Mittelmann's MINLP run it handles only the quadratic instances of the set, so its 96 is a count on a subset.GAMS, BARON solver manual, options LPSol and NLPSol, gams.com/latest/docs/S_BARON.html, read 5 October 2026; A. Lundell, J. Kronqvist and T. Westerlund, "The supporting hyperplane optimization toolkit for convex MINLP", Journal of Global Optimization 84 (2022); H. D. Mittelmann, "Latest progress in optimization software", INFORMS Annual Meeting, Atlanta, 28 October 2025, slide 21 ("SHOT can only handle quadratic instances"), plato.asu.edu/talks/informs2025.pdf. The result files behind the February 2026 page come from four different GAMS distributions, with the sub-solver versions that each distribution shipped. BARON ran under GAMS 52.3, SCIP under 53.1, SHOT under 49.6, and LINDO under 52.1 for the page and 53.1 for the table.H. D. Mittelmann, result and trace files of the MINLP benchmark, plato.asu.edu/ftp/minlp_res/ and plato.asu.edu/ftp/minlp_trc/, file headers read 5 October 2026. A column therefore measures a solver together with its sub-solvers, its thread setting and its interface. A change in any of them moves the number.
Software stacks behind Mittelmann's MINLP benchmark, February 2026
GAMS 52.3 --> BARON --+--> LP, MIP, QP: CPLEX, Xpress, CLP/CBC
| or HSL LA04 (chosen automatically
| by default)
'--> NLP: MINOS, SNOPT, IPOPT, FilterSD
or FilterSQP
GAMS 53.1 --> SCIP ---+--> LP: CPLEX or SoPlex
'--> NLP: Ipopt
GAMS 49.6 --> SHOT ------> Gurobi, CPLEX or Cbc do the work of
its outer approximation; in this run
only the quadratic instances, so its
96 is a count on a subset
GAMS 52.1 (page), 53.1 (table) --> LINDO
a column measures the solver with its sub-solvers, its thread
setting and its interface: a change in any of them moves it
(Performance variability) Run times also move when nothing mathematical changes. Permuting the rows or columns of an instance, changing the random seed, or running on a different machine changes the path a solver takes, because ties are broken differently and floating-point sums come out differently. The same instance can then take seconds or hours. Koch and co-authors defined a variability score for MIPLIB 2010 as the coefficient of variation of the running times over permutations, the standard deviation of the times divided by their mean, and Lodi and Tramontani surveyed the phenomenon.T. Koch et al., "MIPLIB 2010", Mathematical Programming Computation 3 (2011), 103–163; A. Lodi and A. Tramontani, "Performance variability in mixed-integer programming", INFORMS TutORials in Operations Research (2013), 1–12. The practical consequence is that a change of a few percent in a mean is within the noise. A feature is therefore evaluated over several seeds or permutations and on subsets by difficulty. The SCIP release reports use three seeds and the brackets of instances needing at least 100 and at least 1,000 seconds, and the Xpress Global paper reports "affected" instances separately.C. Hojny et al., "The SCIP Optimization Suite 10.0", arXiv 2511.18580 (2025), Section 2; P. Belotti, T. Berthold, T. Gally, L. Gottwald and I. Pólik, "Solving MINLPs to global optimality with FICO Xpress Global", Optimization Online (July 2025), Section 4. A speed-up of 5% measured on one seed and one run is within the noise and should not be reported.
(Who publishes, since 2018) Whose results are published is itself a convention, and in MILP it changed in 2018. At the INFORMS Annual Meeting of that year Gurobi presented benchmark claims against CPLEX and Xpress. It retracted them in a published statement two weeks later: "we published analytics claiming Gurobi was faster, as compared to CPLEX and Xpress, than it actually is … We apologize to Prof. Mittelmann for this misleading characterization of his involvement".Gurobi Optimization, announcement of 7 November 2018, plato.asu.edu/ftp/apology.pdf; H. D. Mittelmann, benchmark index page, plato.asu.edu/bench.html, read 5 October 2026, for the withdrawal history quoted next. Mittelmann's index page records the sequel. IBM and FICO demanded that results for CPLEX and Xpress be removed. In August 2024 Gurobi withdrew from the benchmarks as well, and on 24 December 2024 MindOpt followed. The consequence is that no independent public comparison of Gurobi, CPLEX and Xpress exists today. On the MIPLIB 2017 page of 7 July 2026 the leading published column is COPT 8.0.3, with scaled mean 1.00 and 219 of 240 solved. HiGHS 1.15.1 in its parallel build scores 5.44 with 192 solved, and SCIP 10.0.0 scores 9.93 with 136 solved. The machine is an AMD Ryzen 9 5900X with 12 threads and a two-hour limit.Mittelmann, MIPLIB 2017 benchmark page, run of 7 July 2026, as above. The other published columns are Optverse 2.0.1 (scaled mean 1.72, 210 solved), a FiberSCIP-with-HiGHS portfolio (5.15, 174), a forthcoming concurrent SCIP 11 (6.59, 153) and HiGHS 1.15.0 (7.55, 158). Gurobi's own statement that version 13.0 is "~16% faster" on difficult MIP models and "more than 2X faster on non-trivial MINLPs" than 12.0 is a vendor claim on an unpublished test set. This series labels every such number as one.Gurobi Optimization, "What's New in Gurobi 13.0", gurobi.com/whats-new-gurobi-13-0/, vendor page, read 5 October 2026.
(Tuning to the instances) The newest threat to a benchmark is tuning to its instances. Mittelmann's INFORMS talk of 28 October 2025 ends its "Recent Developments" section with a slide headed "The Benchmarks at a Turning Point". Its first line is "We are in the age of AI". The slide then lists the possible use of over-tuning and machine learning, and notes that the benchmarks are exposed to it except those with undisclosed instances. It announces that Optverse's numbers will be removed after the conference because the code is unavailable, and asks of the coming MIPLIB 2024, "will it be AI-proof?". Two further bullets concern GPUs and the growing number of solvers that handle nonconvexity globally.Mittelmann, INFORMS Annual Meeting talk of 28 October 2025, slide 8; the comparison on known and unknown datasets is slide 14 ("This comparison attempts to sense any tuning effort"). His defences are the ones already mentioned: perturbed and presolved copies of the MIPLIB instances, and undisclosed instances (16 of the 65 in the LP feasibility benchmark). He also compares each code's score on the public and on the undisclosed instances side by side, which "attempts to sense any tuning effort". Section 5.5 returns to the same worry from the learning side.
The rules under which this series quotes a benchmark number follow from the above and are stated once.
- The aggregate is the shifted geometric mean of run times with shift 10 s, scaled by the smallest mean in the same table, with time-outs and failures charged at the time limit.
- A solved count accompanies every mean, because a time-out is censored data and a mean of censored times depends on the limit.
- Every number carries the page date, the machine, the thread count and the time limit, because the pages differ in all four. Solver versions and instance counts are part of the number.
- Comparisons of versions or features report several seeds or permutations and subsets by difficulty, because performance variability is of the order of a few percent.
- A comparison across solver generations states the virtual-best construction and never divides a time limit by a mean time.
- A virtual best is quoted only when it has been computed from per-instance results and the computation is on record.
- A number from a vendor's page is a vendor claim and is labelled as one.
An alternative to any mean is the performance profile of Dolan and Moré, the empirical distribution over instances of each solver's time ratio to the best. It shows the whole distribution instead of one number and is the usual complement to a table. The PAVER tool of Bussieck, Dirkse and Vigerske automates both for GAMS trace files.E. D. Dolan and J. J. Moré, "Benchmarking optimization software with performance profiles", Mathematical Programming 91 (2002); M. R. Bussieck, S. P. Dirkse and S. Vigerske, "PAVER 2.0: an open source environment for automated performance analysis of benchmarking data", Journal of Global Optimization 59 (2014).
For GPU work two further conventions matter, and both are set out in Section 7. Mittelmann's GPU tracks (LP feasibility, convex QP, sparse SDP) give GPU codes a shorter time limit than CPU codes, 1,000 against 15,000 seconds in the LP feasibility benchmark. They report the device and its memory, and mark instances that failed for memory. The measure of a heuristic, GPU or not, is the primal integral of Section 2.3 rather than a solved count, since a heuristic proves nothing. No benchmark of GPU codes on MINLPLib exists, and Section 7.8 states the protocol one would need.
Global solvers that finish
A solver finishes a nonconvex MINLP when it returns a feasible point together with a certificate that no feasible point is better by more than the tolerance. In 2001 that was a research result for any but the smallest instances. Today a dozen codes do it routinely, and a benchmark ranks them. This subsection names the codes, says where each came from, records what the benchmarks of 2026 show, and says what each code does at a node. The codes fall into three families, and the table at the end of the lineage sorts them. They are the convex-only descendants of the methods of Section 3.4, the global spatial branch-and-bound codes, and the extensions that the MILP vendors have added to their engines since 2013. Four guarantee classes recur, and the table's guarantee column uses them, qualified by the problem class where a solver's guarantee is narrower.
Definition 5.3.1 (guarantee classes). A solver is global for a class of problems if, given finite bounds on the variables in nonconvex terms and a tolerance \(\varepsilon > 0\), it terminates on every instance of the class with an \(\varepsilon\)-optimal point or a proof of infeasibility. \(\varepsilon\)-optimality is that of Definition 1.5.20, with the tolerances tabulated in Section 3.5. It is convex-only if it has that guarantee for convex MINLP and, on a nonconvex instance, returns a feasible point and a dual bound that may not close. It is local if it returns a point satisfying first-order conditions (Section 1.3) and no bound, and a heuristic if it returns a feasible point and no bound.
On the table at the end of the lineage the four classes read as follows. BARON is global: on any instance with finite bounds it stops with a point and a bound within \(\varepsilon\) of each other, or a proof that there is no point. SHOT is convex-only: on a convex MINLP it has the same guarantee, on a nonconvex one it returns a feasible point and a bound that may stay apart, since the cuts it adds can be invalid there and are repaired rather than proved. DICOPT on a nonconvex instance is local: its answer is a point at which the first-order conditions hold and nothing is said about any other point. And cuOpt's GPU side is a heuristic: a feasible point, with the bound left to the CPU tree. The word global is therefore a promise about termination on a class, not about speed on an instance.
(The convex-case lineage) The first MINLP codes were for the convex case and grew out of the methods of Section 3.4. DICOPT, from Grossmann's group at Carnegie Mellon, implements the outer-approximation algorithm of Viswanathan and Grossmann and ships in GAMS. In that algorithm nonlinear equalities are relaxed to inequalities according to the signs of their multipliers, and the linearizations are enforced through a penalty term. AlphaECP implements Westerlund's extended cutting plane method and its extension to pseudoconvex constraints. SBB in GAMS and Leyffer's MINLP_BB are NLP-based branch-and-bound codes.GAMS solver manuals for DICOPT, AlphaECP and SBB, gams.com/latest/docs/S_DICOPT.html and the sibling pages; J. Viswanathan and I. E. Grossmann, "A combined penalty function and outer-approximation method for MINLP optimization", Computers & Chemical Engineering 14 (1990), 769–782; M. A. Duran and I. E. Grossmann, "An outer-approximation algorithm for a class of mixed-integer nonlinear programs", Mathematical Programming 36 (1986); T. Westerlund and F. Pettersson, "An extended cutting plane method for solving convex MINLP problems", Computers & Chemical Engineering 19 (1995); T. Westerlund and R. Pörn, "Solving pseudo-convex mixed integer optimization problems by cutting plane techniques", Optimization and Engineering 3 (2002); S. Leyffer, "Integrating SQP and branch-and-bound for mixed integer nonlinear programming", Computational Optimization and Applications 18 (2001). Bonmin (2008) packaged the four strategies of Section 3.4, NLP-based branch and bound, outer approximation, the single-tree LP/NLP method and a hybrid, as open source on COIN-OR. Minotaur (2021) has both families, and Juniper (2018) is an NLP-based branch and bound in Julia. MindtPy and GDPopt bring outer approximation, the extended cutting plane method (both Section 3.4) and generalized Benders decomposition (Section 3.7) to Pyomo models.P. Bonami, L. T. Biegler, A. R. Conn, G. Cornuéjols, I. E. Grossmann, C. D. Laird, J. Lee, A. Lodi, F. Margot, N. Sawaya and A. Wächter, "An algorithmic framework for convex mixed integer nonlinear programs", Discrete Optimization 5 (2008); A. Mahajan, S. Leyffer, J. Linderoth, J. Luedtke and T. Munson, "Minotaur: a mixed-integer nonlinear optimization toolkit", Mathematical Programming Computation 13 (2021); O. Kröger, C. Coffrin, H. Hijazi and H. Nagarajan, "Juniper: an open-source nonlinear branch-and-bound solver in Julia", CPAIOR 2018, LNCS 10848; D. E. Bernal, Q. Chen, F. Gong and I. E. Grossmann, "Mixed-integer nonlinear decomposition toolbox for Pyomo (MindtPy)", Computer Aided Chemical Engineering 44 (2018). Pajarito and Pavito add conic outer approximation (Section 4.8).C. Coey, M. Lubin and J. P. Vielma, "Outer approximation with conic certificates for mixed-integer convex problems", Mathematical Programming Computation 12 (2020), 249–293. SHOT (2022) is the extended supporting hyperplane method inside a MILP solver's single tree. It was the fastest code in the sixteen-solver comparison of Kronqvist, Bernal, Lundell and Grossmann quoted in Section 3.4, and the review attributes the lead of SHOT and AOA to "a single-tree approach closely integrated with the MILP solver".A. Lundell, J. Kronqvist and T. Westerlund, "The supporting hyperplane optimization toolkit for convex MINLP", Journal of Global Optimization 84 (2022); J. Kronqvist, D. E. Bernal, A. Lundell and I. E. Grossmann, "A review and comparison of solvers for convex MINLP", Optimization and Engineering 20 (2019), whose solved counts are quoted in Section 3.4. Knitro's and MOSEK's mixed-integer modes are of the same kind: branch and bound over convex NLP or conic relaxations, exact for convex problems and local or heuristic otherwise.Artelys, Knitro User Guide, "Mixed-integer nonlinear programming" and "Multi-start", artelys.com/docs/knitro; MOSEK appears only in the convex and conic columns of Mittelmann's tables cited below.
(The global line: BARON) The global line begins with BARON. Sahinidis released it in the 1990s as the implementation of branch and reduce, the name Ryoo and Sahinidis gave to spatial branch and bound with domain reduction (Section 2.6) applied at every node.N. V. Sahinidis, "BARON: a general purpose global optimization software package", Journal of Global Optimization 8 (1996); H. S. Ryoo and N. V. Sahinidis, "A branch-and-reduce approach to global optimization", Journal of Global Optimization 8 (1996). BARON is sold by The Optimization Firm; Sahinidis moved from Illinois to Carnegie Mellon and then to Georgia Tech. Tawarmalani and Sahinidis made polyhedral relaxations, linear outer approximations of the envelopes solved by an LP code, the default in 2005. Khajavirad and Sahinidis added in 2018 a hybrid that switches between the LP relaxation and a convex NLP relaxation under learned rules. They report that it raised the number of problems solvable within 500 seconds by 30% over four instance libraries. Kılınç and Sahinidis added in 2018 the MILP devices that exploit integrality.M. Tawarmalani and N. V. Sahinidis, "A polyhedral branch-and-cut approach to global optimization", Mathematical Programming 103 (2005); A. Khajavirad and N. V. Sahinidis, "A hybrid LP/NLP paradigm for global optimization relaxations", Mathematical Programming Computation 10 (2018); M. R. Kılınç and N. V. Sahinidis, "Exploiting integrality in the global optimization of mixed-integer nonlinear programming problems with BARON", Optimization Methods and Software 33 (2018). BARON leads Mittelmann's MINLP benchmark today and has for most runs over the past decade. The exception is Octeract, which topped the then 87-instance table in July 2022 and April 2023 before leaving GAMS and the benchmark.H. D. Mittelmann, INFORMS Annual Meeting 2023 talk, slide 17, plato.asu.edu/talks/informs2023.pdf; M. Miltenberger's interactive history of the benchmarks, mattmilten.github.io/mittelmann-plots/.
(The open-source global codes) Couenne (2009) is the open-source spatial branch and bound that Belotti, Lee, Liberti, Margot and Wächter built on COIN-OR's Bonmin, Cbc, Clp and Ipopt. Their paper fixed the vocabulary of bound tightening and branching points used in Section 3.5.P. Belotti, J. Lee, L. Liberti, F. Margot and A. Wächter, "Branching and bounds tightening techniques for non-convex MINLP", Optimization Methods and Software 24 (2009). SCIP has solved nonconvex MINLPs by spatial branch and bound for more than a decade, since versions 1.2 and 2.0 of 2009–10 and Vigerske's thesis of 2013. Version 8.0 (2022) rebuilt the nonlinear machinery around one expression framework with "nonlinear handlers" that detect structure (convex and concave pieces, quadratics, perspective structure, products) and attach the matching relaxation to each piece.S. Vigerske, Decomposition of Multistage Stochastic Programs and a Constraint Integer Programming Approach to Mixed-Integer Nonlinear Programming, PhD thesis, Humboldt-Universität zu Berlin (2013); S. Vigerske and A. Gleixner, "SCIP: global optimization of mixed-integer nonlinear programs in a branch-and-cut framework", Optimization Methods and Software 33 (2018); Bestuzheva et al. (2025), as above. Version 10.0 (24 November 2025) is 6% faster than 9.0 on SCIP's MINLP test set and 20% faster on the instances that need at least 1,000 seconds (9% and 22% against 9.2.4). It brings the exact solving mode of Section 5.6 and symmetry handling extended to reflections. It also adds a form of conflict analysis (Section 3.6) that reasons on cuts rather than on clauses, and pseudocosts inherited from ancestor nodes when a variable has none of its own. A probabilistic lookahead rule decides when strong branching can stop; pseudocosts and strong branching are the two branching devices of Section 3.1.Hojny et al. (2025), Sections 2 and 3. The release also adds implied-integer detection, flower inequalities for multilinear terms and kernel-search heuristics, which this post does not describe. Version 10.1.0 (18 September 2026) adds a -t option that runs several SCIP configurations concurrently, and the development branch for version 11 adds Feasibility Jump and local-search heuristics.SCIP releases on GitHub, v10.0.0 (24 November 2025) to v10.1.0 (18 September 2026), github.com/scipopt/scip/releases; SCIP CHANGELOG, master branch, read 4 October 2026. GloMIQO (2013) and ANTIGONE (2014) are Misener and Floudas's global solvers: a reformulation to a standard factorable form, term-wise relaxations with RLT cuts (Section 4.7) and edge-concave cuts, the cuts for functions concave along every coordinate direction (Section 2.4).R. Misener and C. A. Floudas, "GloMIQO: global mixed-integer quadratic optimizer", Journal of Global Optimization 57 (2013); R. Misener and C. A. Floudas, "ANTIGONE: Algorithms for coNTinuous / Integer Global Optimization of Nonlinear Equations", Journal of Global Optimization 59 (2014). LINDO's global solver combines branch and bound on convex relaxations with multistart local search.Y. Lin and L. Schrage, "The global solver in the LINDO API", Optimization Methods and Software 24 (2009). MAiNGO, from Mitsos's group at RWTH Aachen, builds McCormick relaxations in the original variable space through the MC++ library, offers the reduced-space formulations of Section 2.4, and runs its tree in parallel with MPI.D. Bongartz, J. Najman, S. Sass and A. Mitsos, "MAiNGO – McCormick-based Algorithm for mixed-integer Nonlinear Global Optimization", technical report, Process Systems Engineering (AVT.SVT), RWTH Aachen University (2018), permalink.avt.rwth-aachen.de/?id=729717; the MAiNGO documentation lists C. Witte among the software's authors. Alpine solves nonconvex MINLPs by the adaptive piecewise McCormick partitioning of Section 4.6, and RAPOSa is a global solver for polynomial problems built on the reformulation-linearization technique of Section 4.7.H. Nagarajan, M. Lu, S. Wang, R. Bent and K. Sundar, "An adaptive, multivariate partitioning algorithm for global optimization of nonconvex programs", Journal of Global Optimization 74 (2019); B. González-Rodríguez, J. Ossorio-Castillo, J. González-Díaz, Á. M. González-Rueda, D. R. Penas and D. Rodríguez-Martínez, "Computational advances in polynomial optimization: RAPOSa, a freely available global solver", Journal of Global Optimization 85 (2023). Octeract solved all 87 instances of the MINLP benchmark of 2023 (BARON 77, SCIP 64, ANTIGONE 53, LINDO 42, Couenne 24, on an Intel i7-11700K with a 7,200-second limit). It was frozen at version 4.7.1 when it left GAMS, was removed from GAMS in version 46 (February 2024), and is sold through AMPL and AIMMS. Its vendor describes a distributed engine, a claim that no benchmark tests.Mittelmann, INFORMS 2023 talk, slide 17 ("Since Octeract will be removed from GAMS, it will be frozen at version 4.7.1"); GAMS release notes, distribution 46 (February 2024); Octeract Ltd, "Octeract Engine" product page, octeract.com/octeract-engine/, vendor claims.
(The MILP vendors arrive) The larger change of the last decade is that the MILP vendors arrived, in the following order. CPLEX 12.6 (2013–14) solved nonconvex mixed-integer quadratic programs to global optimality, and Gurobi 9.0 (November 2019) added nonconvex quadratic constraints and objectives, translated into bilinear form and solved by spatial branching (Section 4.7 gives the sources for both). FICO's Xpress Global (October 2022, with Xpress 9.0) was the first general MINLP solver from a MILP vendor. It is one of the two columns whose publication Mittelmann's MINLP page suppresses.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 (July 2025), optimization-online.org/2025/07/solving-minlps-to-global-optimality-with-fico-xpress-global/; the paper's introduction dates the solver's release to October 2022. Gurobi 11.0 (28 November 2023) added spatial branching with dynamic outer approximation for its function constraints as an option, through the new parameter FuncNonlinear, with the static piecewise-linear approximation still the default. Version 12.0 (12 November 2024) made the dynamic approach the default and added general nonlinear constraints \(y = f(x)\) given as expression trees, solved "to global optimality using a branch-and-bound algorithm". Version 13.0 (11 November 2025) deprecates the function constraints in favour of the general nonlinear constraints (Section 4.6), adds a nonlinear barrier method, as a preview, for local optima of continuous nonconvex models, and adds a PDHG method for LP that can run on a GPU (Section 7.2).Gurobi Optimizer Reference Manual, release notes for 11.0, 12.0 and 13.0, docs.gurobi.com; the release dates are the upload dates of gurobipy 11.0.0, 12.0.0 and 13.0.0 on PyPI, read 5 October 2026. COPT 8.0 (January 2026) added global solving of nonconvex quadratically constrained problems.D. Ge, Q. Huangfu, Z. Wang, J. Wu and Y. Ye, "Cardinal Optimizer (COPT) user guide", arXiv 2208.14314 (v4, January 2026). Gurobi's table of 13.0 against 12.0 (2.5 times faster on MINLP, 54.7% on nonconvex MIQCP) is a vendor claim on an unpublished test set.Gurobi Optimization, "What's New in Gurobi 13.0", vendor page, as above.
The MILP vendors arrive, in order of release
2013-14 o CPLEX 12.6: nonconvex MIQP to global optimality
|
Nov 2019 o Gurobi 9.0: nonconvex quadratic constraints and
| objectives, in bilinear form, spatial branching
|
Oct 2022 o Xpress Global, with Xpress 9.0: the first general
| MINLP solver from a MILP vendor
|
28 Nov 2023 o Gurobi 11.0: spatial branching with dynamic outer
| approximation for its function constraints, as
| an option (FuncNonlinear)
|
12 Nov 2024 o Gurobi 12.0: the dynamic approach by default;
| general nonlinear constraints y = f(x)
|
11 Nov 2025 o Gurobi 13.0: function constraints deprecated;
| nonlinear barrier (preview, local); PDHG for LP,
| can run on a GPU
|
Jan 2026 o COPT 8.0: nonconvex quadratically constrained
problems solved globally
(The Xpress Global ablations) The Xpress Global paper is the one published account of how a MILP vendor's engine becomes a global solver, and its ablations extend the measurements of Section 5.1 to the nonlinear case. The solver is the Xpress MILP code with a nonlinear presolve, convexification cuts and spatial branching added, and FICO's successive linear programming code as a local heuristic. Successive linear programming is a local NLP method that linearizes the objective and the constraints at the current point and solves a sequence of LPs inside a trust region, a bound on the step that keeps each LP's answer where the linearization is still a fair model. The convexification cuts are RLT cuts (Section 4.7) and lifted tangent inequalities, which are tangent cuts on convex pieces strengthened by the variable bounds. On the authors' internal test set, disabling presolve loses about 8% of the solved instances with a 40% slowdown on the rest, and formula simplification alone accounts for 82 fewer solved instances. Removing the RLT cuts loses 34 solved instances with only a 1% effect on time. Evaluating strong-branching candidates without convexification cuts loses 583 solved instances, and without propagation and cuts 689. Using them costs 15 to 30% in time per candidate, "which leads to much better performance, especially in the number of solved instances". The local heuristic matters on hard instances, with 102 fewer solved without it.Belotti, Berthold, Gally, Gottwald and Pólik (2025), Sections 4.5 to 4.7 and Table 5. A deterministic parallel run follows the same path whatever the thread timing, and an opportunistic one does not (Section 6.2 makes this precise). Sixteen deterministic threads give a speed-up of 2.59 over all solvable instances and 3.97 on the instances that need at least 100 nodes, and 32 threads add nothing over all instances (3.97 to 4.05 on the node-heavy subset, a figure Section 6.3 returns to). The opportunistic mode is only 2 to 9% faster than the deterministic one, which is why deterministic parallelism is the default.Belotti et al. (2025), Tables 7 and 8 and Section 4.9. Section 6 takes up determinism and the scaling of trees. Against the MILP measurements of Section 5.1 the nonlinear solver differs in one place. The nonlinear presolve contributes more than the linear one, because a factorable reformulation (Section 2.4) is only as good as the expression it is given.
Xpress Global: a MILP engine made global, and its ablations
Xpress MILP code
+ nonlinear presolve off: about 8% of the solved instances
lost, 40% slower on the rest; formula
simplification alone: 82 fewer solved
+ convexification cuts
RLT cuts off: 34 fewer solved, 1% in time
lifted tangent cuts
+ spatial branching
+ SLP local heuristic off: 102 fewer solved (hard instances)
= Xpress Global
strong-branching candidates evaluated
without convexification cuts ............. 583 fewer solved
without propagation and cuts ............. 689 fewer solved
(using them costs 15 to 30% in time per candidate)
ablations on the authors' internal test set
deterministic threads (the default)
16 threads speed-up 2.59 over all solvable instances,
3.97 on those that need at least 100 nodes
32 threads nothing more over all instances; 4.05 on the
node-heavy subset
opportunistic mode: only 2 to 9% faster than deterministic
The table collects the field. The guarantee column uses Definition 5.3.1. The sub-solver column names the codes that do the LP, MILP and NLP work inside each solver, because, as Section 5.2 said, a benchmark column measures the stack.
| solver | origin, licence | guarantee | relaxation and search at a node | sub-solvers |
|---|---|---|---|---|
| BARON | Sahinidis; The Optimization Firm; commercial | global | branch and reduce: FBBT, OBBT, probing; polyhedral LP relaxation of the envelopes, convex NLP relaxation when it pays; local NLP for incumbents; violation-transfer branching | LP/MIP/QP: CPLEX, Xpress, CLP/CBC, HSL LA04 (automatic); NLP: MINOS, SNOPT, IPOPT, FilterSD, FilterSQP |
| SCIP | Zuse Institute Berlin; Apache 2.0 since 8.0.3 | global (MINLP); exact mode for MILP | constraint integer programming: nonlinear handlers add tangents, secants, McCormick and perspective cuts; OBBT at the root; hybrid spatial branching; forty heuristics | LP: SoPlex or CPLEX; NLP: Ipopt; FiberSCIP/UG for parallel trees |
| Couenne | COIN-OR; EPL, open | global | spatial B&B on Bonmin/Cbc: FBBT, OBBT with a depth schedule, reduced-cost tightening; rounding heuristics | LP: Clp; NLP: Ipopt |
| ANTIGONE | Misener and Floudas; via GAMS | global | standard form; term-wise relaxations, RLT and edge-concave cuts; OBBT | LP/MILP: CPLEX; NLP: CONOPT or SNOPT |
| LINDO Global | LINDO Systems; commercial | global | convex relaxations in a tree; multistart local search | own |
| Octeract | Octeract; commercial; frozen at 4.7.1 when it left GAMS, removed from GAMS 46 (February 2024); sold through AMPL and AIMMS | global | spatial B&B; vendor claims a distributed engine | own (left Mittelmann's benchmark after April 2023) |
| MAiNGO | RWTH Aachen; EPL-2.0 (maingopy) | global | McCormick relaxations in the original or reduced space (MC++); interval bounds | LP: CPLEX or CLP; MPI tree; GPU interval bounder as a prototype (Section 6.4) |
| Alpine | Julia, open | global | adaptive piecewise McCormick partitioning | MILP and NLP solvers of the user's choice |
| RAPOSa | Univ. Santiago de Compostela; free | global (polynomial) | RLT relaxations with spatial branching | LP solver of the user's choice |
| Gurobi 13 | Gurobi Optimization; commercial | global (nonconvex MIQCP, expression trees) | bilinear form with McCormick and spatial branching; dynamic outer approximation of univariate pieces; nonlinear barrier (local) | own LP, QP, barrier, PDHG (GPU optional) |
| Xpress Global | FICO; commercial | global | Xpress MILP engine plus nonlinear presolve, RLT and lifted tangent cuts, spatial branching | own LP and MILP; SLP local solver |
| COPT 8 | Cardinal Operations; commercial | global (nonconvex (MI)QCQP) | spatial B&B over the COPT MILP engine | own LP (CPU simplex, barrier; GPU PDLP and barrier for LP) |
| CPLEX 22.1 | IBM; commercial | global (nonconvex MIQP only) | MIQP branch and bound; a learned classifier decides whether to linearize | own |
| SHOT | Åbo Akademi; COIN-OR, EPL-2.0 | convex-only; heuristic with repair on nonconvex problems | extended supporting hyperplanes in a MILP solver's tree; invalid cuts repaired | MILP: Gurobi, CPLEX, Cbc, HiGHS; NLP: Ipopt |
| Bonmin | COIN-OR; EPL, open | convex-only | NLP-B&B, OA, LP/NLP-B&B and a hybrid | LP: Clp or Cbc; NLP: Ipopt |
| DICOPT | Grossmann group; GAMS | convex-only (local on nonconvex) | outer approximation with equality relaxation and augmented penalty | MILP and NLP solvers of GAMS |
| Knitro | Artelys; commercial | convex-only (local on nonconvex) | NLP-based and LP/NLP-based branch and bound; parallel deterministic multistart | own NLP, MILP |
| MOSEK | MOSEK ApS; commercial | convex conic only | mixed-integer conic branch and bound | own conic interior point |
| Minotaur | Argonne, Wisconsin; open | convex-only | NLP-based and LP/NLP-based branch and bound | Ipopt, filterSQP; CLP or CPLEX |
| Juniper | Julia, open | convex-only (local on nonconvex) | NLP-based branch and bound | Ipopt or any JuMP NLP solver |
| MindtPy, GDPopt | Pyomo; open | convex-only (local on nonconvex) | OA, ECP, GBD; logic-based OA for GDP | any Pyomo MILP and NLP solvers |
| cuOpt 26.08 | NVIDIA; Apache 2.0 | LP and MILP only (no nonlinear terms) | GPU primal heuristics with a CPU branch and bound; QP, with conic constraints in beta | own: GPU PDLP and barrier, CPU dual simplex |
The scoreboard follows the rules of Section 5.2. Every number is a count of instances solved within the limit, dated by its page, with the machine, limit and thread count where the page gives them.H. D. Mittelmann, benchmark pages read on 4 and 5 October 2026: MINLP (page dated 26 February 2026), plato.asu.edu/ftp/minlp.html; binary nonconvex QPLIB (9 May 2026), plato.asu.edu/ftp/qplib.html; discrete non-binary nonconvex QPLIB (12 May 2026), plato.asu.edu/ftp/nonbinary.html; convex discrete QPLIB (7 September 2026), plato.asu.edu/ftp/convex.html; MISOCP (10 September 2026), plato.asu.edu/ftp/misocp.html; QUBO (21 September 2026), plato.asu.edu/ftp/qubo.html. The virtual best of the MINLP row is the computation of Section 5.2. One page is left out on purpose: the continuous nonconvex QPLIB page of 17 May 2026, plato.asu.edu/ftp/cnconv.html, lists its solvers in one order and heads its per-instance table in another, so that the solved counts 35, 31, 22, 15 and 28 attach to different solvers under the two readings; a recount from its per-solver logs favours the table-header order, but its SCIP column comes from logs of April 2025 with SCIP 9.2.1 while the page names SCIP 10.0.1, so nothing from that page is quoted here. The last row is the QUBO benchmark, quadratic unconstrained binary optimization, \(\min\{x^\top Q x : x \in \{0, 1\}^n\}\), which Section 7.6 takes up. The table after the figure prints the same counts with the same conditions, so that a count can be quoted with everything that qualifies it.
| benchmark, instances; page date; machine; limit; threads | solved |
|---|---|
| MINLP, 200 MINLPLib instances; 26 Feb 2026; AMD Ryzen 9 5900X; 2 h; through GAMS; feasibility tolerance 1e-6 | BARON 158, SCIP 154, LINDO 116, SHOT 96 (quadratic instances only); at least one of the four 183 (from compare.txt, 6 Mar 2026); Gurobi and Xpress run, publication suppressed |
| binary nonconvex QPLIB, 128; 9 May 2026; Ryzen 9 5900X; 1 h; 12 threads (SCIP and Couenne single-threaded); zero MIP gap | SHOT with Gurobi 91, COPT 8.0.4 83, BARON 25.11.17 74, RAPOSa 4.4.1 69, SCIP 10.0.0 36, ANTIGONE 1.1 16, Couenne 0.5 6 |
| discrete non-binary nonconvex QPLIB, 160; 12 May 2026; Intel Xeon Gold 6230; 3 h; 8 threads | SHOT with Gurobi 99, COPT 88, BARON 25.3.19 66, SCIP 10.0.0 40 |
| convex discrete QPLIB, 31; 7 Sep 2026; machine and limit not recorded for this post | SHOT 25, COPT 24, BARON 22, MOSEK 20, Knitro 15, SCIP 14, Minotaur 8, Bonmin 7 |
| MISOCP, 47; 10 Sep 2026; Intel i7-11700K; 1 h; MIP gap 0 (as quoted in Section 4.8) | COPT 47, SCIP 39, MOSEK 37, Knitro 33 |
| QUBO (QPLIB), 23; 21 Sep 2026; machine not recorded for this post; 1 h; 12 threads; solved globally | QuBowl 22, QUPLANE 22, BARON 13, COPT 13, McSparse 12, SHOT 12, BiqBin 9, SCIP 8 |
(Dating the MINLP counts) The MINLP numbers move between runs and must be dated. The INFORMS slide of 2023 (87 instances, Intel i7-11700K) read Octeract 87, BARON 77, SCIP 64. The INFORMS slide of 28 October 2025, on the current 200 instances and the current machine, read BARON 161, SCIP 153, LINDO 116 and SHOT 96. The page of 26 February 2026 reads 158, 154, 116 and 96.Mittelmann, INFORMS 2023 talk, slide 17; INFORMS 2025 talk, slide 21; MINLP benchmark page of 26 February 2026, as above. The 200 instances are the test set of the SCIP 8 paper. That paper selected them from the 1,505 MINLPLib instances every solver could read so that at least one solver solved each, the set was varied in integrality and nonlinearity, and no family of similar names dominated. It ran each instance in five variants, the original and four permutations of its rows and columns, for 1,000 serial runs under GAMS 41.2 with a two-hour limit. BARON solved 790, SCIP 776, Octeract 671 and Lindo API 538, the virtual best 967 and the virtual worst 368. SCIP's 41 failures were 16 wrong optimal values, 23 infeasible solutions and 2 aborts. BARON's 27 were 26 wrong values and one infeasible solution. With threads, on the 200 unpermuted instances, BARON solved 161, 160, 160 and 158 at 1, 4, 8 and 16 threads and FiberSCIP 161, 145, 147 and 152.Bestuzheva et al., Journal of Global Optimization 91 (2025), Sections 3.1 to 3.2, Table 2 (p. 305) and Table 3; Table 2 is identical to Table 2 of arXiv 2301.00587 v1 (January 2023). The conclusion that "SCIP's performance is currently on par with the state-of-the-art commercial solver BARON" is the authors'. The thread columns are the first sign of the theme of Section 6: more threads did not mean more solved instances for either code.
The SCIP 8 paper's test set and runs (Bestuzheva et al.)
1,505 MINLPLib instances that every solver could read
| selected so that at least one solver solves each, the
| set is varied in integrality and nonlinearity, and no
| family of similar names dominates
v
200 instances
| five variants each: the original and four permutations
| of its rows and columns
v
1,000 serial runs under GAMS 41.2, two-hour limit
|
v
solved virtual best 967 BARON 790 SCIP 776
Octeract 671 Lindo API 538 virtual worst 368
failures BARON 27 = 26 wrong values + 1 infeasible solution
SCIP 41 = 16 wrong optimal values + 23 infeasible
solutions + 2 aborts
with threads, on the 200 unpermuted instances (solved):
threads 1 4 8 16
BARON 161 160 160 158
FiberSCIP 161 145 147 152
Read as a statement about the state of the art, the February 2026 page says this. On a curated set of 200 instances with up to about 100,000 variables, the best single code finishes 79% within two hours, the best of four published codes 91.5%, and 17 instances, 8.5% of the set, are beyond all four. Two commercial codes were run on the same set and their results are not public.Mittelmann, INFORMS 2025 talk, slide 21 (instance sizes); MINLP benchmark page of 26 February 2026 ("*: publication suppressed").
(What each code does at a node) What each code does at a node is Section 3.6's subject, and only the differences that matter for the rest of this series are listed here. BARON runs domain reduction first (FBBT, optimality-based tightening and probing, the three devices of Section 2.6, all on by default). It then solves a polyhedral LP relaxation of the envelopes, or a convex NLP relaxation when its rules say the latter pays, runs a local NLP search for an incumbent, and branches by violation transfer, the spatial branching score of Section 3.5. Its tolerances are in the table of that section.BARON user manual (The Optimization Firm) and GAMS/BARON documentation, options LBTTDo, OBTTDo, PDo, DoLocal, BrVarStra, read 5 October 2026; Y. Puranik and N. V. Sahinidis, "Domain reduction techniques for global NLP and MINLP optimization", Constraints 22 (2017). SCIP's nonlinear handlers separate convex pieces by tangents and concave pieces by the tightest plane valid at the box's vertices, products by McCormick planes, quadratics by their eigen-decomposition, and terms over semicontinuous variables by perspective cuts (Section 4.3). OBBT runs at the root with an iteration budget of ten times the root LP's. Spatial branching scores violation and pseudocosts with equal weight, pulls the branching point three quarters of the way to the midpoint, a pull that Section 3.5 shows shrinking with the relative domain width deeper in the tree, and keeps it a fifth of the width from either end. Its tolerances are also in that table.SCIP source, src/scip/set.c (SCIP 10.0.0), parameters branching/midpull, branching/clamp; Bestuzheva et al. (2025), Section 2. Couenne runs three rounds of FBBT, OBBT on a depth schedule, aggressive bound tightening and reduced-cost tightening. It places the branching point a quarter of the way from the midpoint to the relaxed value, and uses the rounding heuristics of Nannicini and Belotti.Couenne source (COIN-OR, master), options max_fbbt_iter, log_num_obbt_per_level, branch_midpoint_alpha, branch_lp_clamp; G. Nannicini and P. Belotti, "Rounding-based heuristics for nonconvex MINLPs", Mathematical Programming Computation 4 (2012). ANTIGONE relaxes term by term, with piecewise McCormick and edge-concave relaxations and RLT cuts, and applies OBBT with a probability that decreases with depth.Misener and Floudas (2014); Puranik and Sahinidis (2017), Section 7. Gurobi solves nonconvex quadratics in bilinear form with McCormick planes and spatial branching. It solves univariate pieces of expression trees by dynamic outer approximation: tangents where the piece is convex on the node's domain, secants where it is concave, refined by branching. Its gap tolerance is in the same table.Gurobi Optimizer Reference Manual 13.0, parameters NonConvex, FuncNonlinear, and the "Nonlinear constraints" page, docs.gurobi.com. Gurobi's branching and domain-reduction rules are not documented in detail. Xpress Global does what the previous paragraph said. MAiNGO evaluates McCormick relaxations through the expression graph at each node and solves the resulting LP with CPLEX or CLP. SHOT adds supporting hyperplanes inside the MILP solver's tree. When a cut from a constraint it could not prove convex makes the master infeasible, it repairs the master by adding slack to that cut and forces progress with a primal objective cut. On nonconvex problems it keeps a valid dual bound but has no guarantee of closing the gap, which is why its 96 of 200 is not comparable with BARON's 158.A. Lundell and J. Kronqvist, "Polyhedral approximation strategies for nonconvex mixed-integer nonlinear programming in SHOT", Journal of Global Optimization 82 (2022): "there is no theoretical guarantee that the gap can be reduced to the globally optimal solution for nonconvex problems".
SHOT on a nonconvex problem: the repair step
+--------------------------------------------------+
| MILP master: supporting hyperplanes as cuts, |<-----.
| inside the MILP solver's tree | |
+--------------------------------------------------+ |
| |
v |
does a cut from a constraint SHOT could not |
prove convex make the master infeasible? |
| | |
| no | yes |
v v |
continue in the tree add slack to that cut; |
force progress with a ------'
primal objective cut
dual bound: valid; closing the gap: no guarantee, which is why
its 96 of 200 is not comparable with BARON's 158
(The message for GPU work) For GPU work the table has one message. Every solver in the global class runs its tree on the CPU, and all but LINDO and Alpine bound with a polyhedral relaxation solved by a simplex code, BARON switching to a convex NLP relaxation when its rules say so. No production solver runs a spatial branch and bound for general nonconvex MINLP on a GPU. cuOpt is an LP and MILP solver whose GPU side is heuristic. Gurobi 13 and COPT 8 use the GPU for the LP method, not for the tree. MAiNGO's GPU interval bounder and the generated McCormick kernels of ParBB are research prototypes that Section 6.4 and Section 7.8 describe.NVIDIA, cuOpt User Guide 26.08, "Introduction", docs.nvidia.com/cuopt; 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); 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). The codes that finish are the baseline a GPU design has to beat, on the benchmark rules of Section 5.2.
Heuristics that do not need the LP
The incumbent is half of pruning. A node is discarded when its bound is no better than the best feasible point known, so a good point found early shrinks the tree, and the primal integral of Section 2.3 measures how early. Section 3.2 catalogued the heuristics that solvers run, almost all of which solve a linear program or a sub-MIP. This subsection is about the ones that do not. There are two kinds. The first is local search over integer assignments with a violation measure as objective, of which Feasibility Jump is the clearest case. It needs no LP at all and parallelizes in a way nothing built on the simplex method can, which is why it was the first MIP heuristic to run on a GPU. The second is the alternating direction method of multipliers applied to problems whose nonconvexity is separable. It needs a factorization but no tree, and it comes with a bound on how far from optimal its answer can be. The bound, the Shapley–Folkman theorem of Starr and its sharpening by Udell and Boyd, is the reason the portfolio problems of Section 4.4 and Section 9 are nearly convex, and it is stated and proved here.
Example 5.4.1 (one move of Feasibility Jump). The smallest example shows what one move of a local search is. Take the single row \(3x_1 + 2x_2 \ge 5\) over integers \(x_1, x_2 \in [0, 5]\), at the point \(x = (0, 0)\), where it is violated by 5. Moving \(x_1\) to 2 removes the whole violation, and so does moving \(x_2\) to 3. Each move is scored by the violation it removes, which is 5 for both, and either variable could be moved. The paragraph and the definition below make "the best move of one variable" precise, and the example continues after them with a second row that breaks the tie.
(The objects, recalled from Section 3.2) Section 3.2 fixed the objects. The weighted violation \(F^w(\bar x) = \sum_{i \in M} w_i \max\{0,\ a_i^\top \bar x - b_i\}\) of a point \(\bar x\) under row weights \(w \ge 0\), for rows \(a_i^\top x \le b_i\) with an equality written as two inequalities, is the display before Proposition 3.2.10. The one-variable violation \(G_j(t)\), obtained by fixing every variable but \(x_j\) at \(\bar x\), is the function of that proposition, convex and piecewise linear in \(t\) with its breakpoints at the ratios \(d_i / a_{ij}\). The jump value \(\mathrm{Jump}_j(\bar x)\) is the smallest minimizer of \(G_j\) over the admissible values of \(x_j\) other than \(\bar x_j\), and Algorithm 3.2.11 computes it exactly. One quantity is named here because the rest of this subsection and Section 7.7 rank moves by it.
Definition 5.4.2 (score of a move). The score of \(x_j\) at \(\bar x\) is \(s_j = G_j(\bar x_j) - G_j(\mathrm{Jump}_j(\bar x))\), the violation that the best single move of \(x_j\) removes. A variable is a candidate for a move when \(s_j > 0\).
(The score on the example, and the cost per step) In Example 5.4.1 with its single row, \(G_1(0) = 5\) and \(G_1(\mathrm{Jump}_1) = G_1(2) = 0\), so \(s_1 = 5\), and \(s_2 = 5\) in the same way; both variables are candidates. Excluding the current value \(\bar x_j\) from the minimization matters only when \(\bar x_j\) is itself a minimizer of \(G_j\). Then every other value has \(G_j(t) \ge G_j(\bar x_j)\), the score is nonpositive, and the variable is not moved. Proposition 3.2.10 gives the cost of the sequential scheme: \(O(1)\) to choose a move from a sample of at most 25 candidates, \(O(\eta \log \eta)\) to recompute the mover's jump value when it appears in \(\eta\) rows, \(O(\eta \mu)\) to update the scores of the variables that share a row with it, where \(\mu\) is the largest number of variables in a row, and \(O(m \mu)\) for a weight update.
(The algorithm and its record) Feasibility Jump (Algorithm 3.2.12) starts from any assignment and repeatedly moves one variable with positive score to its jump value. At a local minimum of \(F^w\), where no score is positive, it raises by one the weight of every violated row and continues. Luteberget and Sartor describe the method as "a guided local search over a Lagrangian relaxation". The weights play the role of multipliers for the rows, and the weight update is a step of dual ascent on the rows that remain violated: the multiplier of a constraint that is still violated is raised, which is the direction in which the dual function of Section 2.2 increases.B. Luteberget and G. Sartor, "Feasibility Jump: an LP-free Lagrangian MIP heuristic", Mathematical Programming Computation 15 (2023), 365–388; reference implementation at github.com/sintef/feasibilityjump. It won the computational competition of the 2022 MIP Workshop against LP-based entries, with Salvagnin's entry as runner-up.MIP Workshop 2022, computational competition results, mixedinteger.org/2022/competition/, read 5 October 2026. Its record is in Section 3.2: feasible solutions for 123 of the 240 MIPLIB 2017 benchmark instances without an LP solve, a default heuristic in Xpress 9 whose effect shows in the primal integral rather than in the solved count, and adoption by HiGHS, by SCIP's development branch and, on the GPU, by NVIDIA's cuOpt.
Example 5.4.1 continues. Add the row \(x_1 + x_2 \le 2\), also with weight 1. The move \(x_2 \to 3\) now violates it by 1, as much as \(x_2 \to 2\) leaves in the first row, and the tie goes to the smaller value, so the jump value of \(x_2\) becomes 2, with a score of 4. The move \(x_1 \to 2\) still removes everything and keeps its score of 5, so \(x_1\) is chosen. Each score touched only the rows containing the variable. The block below computes the jump values and scores by evaluating \(G_j\) on the integer domain, which is exact here and is what the breakpoint scan of Algorithm 3.2.11 computes without enumeration.
# One Feasibility Jump move (Example 5.4.1).
#
# Rows a_i x <= b_i with weights w_i. The jump value of x_j is the value
# t in its domain, different from the current one, that minimizes the
# weighted violation of the rows containing x_j with the other variables
# held fixed; the score is the violation removed by the move.
import numpy as np
def violation(rows, w, x):
return sum(wi * max(0.0, a @ x - b) for (a, b), wi in zip(rows, w))
def jump(rows, w, x, j, lo, hi):
"""Minimizer of G_j(t) = sum_i w_i max(0, a_ij t - d_i).
Taken over the integers of [lo, hi] other than x_j, the smallest on
ties; here by direct evaluation on the (small) integer domain.
"""
best = None
for t in range(lo, hi + 1):
if t == x[j]:
continue
y = x.copy()
y[j] = t
g = violation(rows, w, y)
if best is None or g < best[0] - 1e-12:
best = (g, t)
return best[1], best[0]
x = np.array([0.0, 0.0])
dom = (0, 5)
# Each case is a title and its rows (a_i, b_i); the row 3x1 + 2x2 >= 5
# is written as -3x1 - 2x2 <= -5.
cases = [("one row: 3x1 + 2x2 >= 5",
[(np.array([-3.0, -2.0]), -5.0)]),
("two rows: 3x1 + 2x2 >= 5 and x1 + x2 <= 2",
[(np.array([-3.0, -2.0]), -5.0),
(np.array([1.0, 1.0]), 2.0)])]
for title, rows in cases:
w = np.ones(len(rows))
cur = violation(rows, w, x)
print(title)
print(f" at x = {x.astype(int).tolist()} "
f"the weighted violation is {cur:.0f}")
for j in range(2):
t, g = jump(rows, w, x, j, *dom)
print(f" jump value of x{j + 1}: {t} "
f"(violation after the move {g:.0f}, score {cur - g:.0f})")
one row: 3x1 + 2x2 >= 5
at x = [0, 0] the weighted violation is 5
jump value of x1: 2 (violation after the move 0, score 5)
jump value of x2: 3 (violation after the move 0, score 5)
two rows: 3x1 + 2x2 >= 5 and x1 + x2 <= 2
at x = [0, 0] the weighted violation is 5
jump value of x1: 2 (violation after the move 0, score 5)
jump value of x2: 2 (violation after the move 1, score 4)
(The same computation, shaped for a device) The performance-shaped version of the same computation is the one a GPU runs. The listing below stores the rows by column, computes every variable's jump value by the breakpoint rule of Proposition 3.2.10, and treats the variables as independent tasks dealt to threads. On the two-row example it prints the moves above. It compiles with clang++ -std=c++23.
// Feasibility Jump: the jump value of every variable, all at once.
//
// The jump values are taken at the current point. Rows a_i x <= b_i are
// stored by column (column j lists the rows it touches). With the other
// variables fixed,
//
// G_j(t) = sum_i w_i max(0, a_ij t - d_i)
//
// is convex and piecewise linear in t, so its minimum over the integers
// of [lo_j, hi_j] is at a breakpoint floor(d_i/a_ij), ceil(d_i/a_ij) or
// at a bound. Excluding the current value x_j changes the minimizer only
// when x_j is itself the minimizer; the score is then nonpositive and the
// variable is not moved, so the candidate set below suffices.
//
// One task per variable; here the tasks are spread over std::jthreads.
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <limits>
#include <thread>
#include <vector>
// One nonzero of a column: its row and its coefficient a_ij.
struct Entry {
int row;
double a;
};
// The rows by column, with right-hand sides b, row weights w and the
// variable domains [lo_j, hi_j].
struct Problem {
std::vector<std::vector<Entry>> col;
std::vector<double> b, w, lo, hi;
};
// Weighted violation of the rows containing x_j when x_j moves from xj
// to t (act = row activities).
static double G(const Problem& P, const std::vector<double>& act, int j,
double xj, double t) {
double g = 0.0;
for (const Entry& e : P.col[j]) {
double excess = act[e.row] - e.a * xj + e.a * t - P.b[e.row];
g += P.w[e.row] * std::max(0.0, excess);
}
return g;
}
// The jump value: the candidate minimizing G, different from x_j, the
// smallest on ties.
static double jump(const Problem& P, const std::vector<double>& act,
const std::vector<double>& x, int j) {
std::vector<double> cand{P.lo[j], P.hi[j]};
for (const Entry& e : P.col[j]) {
// slack of row i without x_j
double d = P.b[e.row] - (act[e.row] - e.a * x[j]);
for (double c : {std::floor(d / e.a), std::ceil(d / e.a)})
if (c >= P.lo[j] && c <= P.hi[j])
cand.push_back(c);
}
std::sort(cand.begin(), cand.end());
double best = std::numeric_limits<double>::infinity();
double arg = x[j];
for (double t : cand) {
if (t == x[j])
continue;
if (double g = G(P, act, j, x[j], t); g < best) {
best = g;
arg = t;
}
}
return arg;
}
int main() {
// Example 5.4.1: rows 3x1 + 2x2 >= 5 (as -3x1 - 2x2 <= -5) and
// x1 + x2 <= 2, weights 1, x = (0, 0), domains [0, 5]
Problem P;
P.col = {{{0, -3.0}, {1, 1.0}},
{{0, -2.0}, {1, 1.0}}};
P.b = {-5.0, 2.0};
P.w = {1.0, 1.0};
P.lo = {0.0, 0.0};
P.hi = {5.0, 5.0};
std::vector<double> x{0.0, 0.0}, act(P.b.size(), 0.0);
for (std::size_t j = 0; j < P.col.size(); ++j)
for (const Entry& e : P.col[j])
act[e.row] += e.a * x[j];
const std::size_t n = P.col.size();
const std::size_t cores = std::thread::hardware_concurrency();
const std::size_t T = std::min<std::size_t>(n, cores);
std::vector<double> v(n), score(n);
// One independent task per variable, dealt round-robin to T threads;
// a GPU gives each its own thread.
{
std::vector<std::jthread> pool;
for (std::size_t t = 0; t < T; ++t)
pool.emplace_back([&, t] {
for (std::size_t j = t; j < n; j += T) {
v[j] = jump(P, act, x, static_cast<int>(j));
score[j] = G(P, act, j, x[j], x[j]) -
G(P, act, j, x[j], v[j]);
}
});
} // jthreads join here
for (std::size_t j = 0; j < n; ++j)
std::printf("x%zu: jump value %g, score %g\n", j + 1, v[j],
score[j]);
return 0;
}
x1: jump value 2, score 5
x2: jump value 2, score 4
(What a GPU variant trades) The cost of one pass is \(O(\mathrm{nnz}(A))\) for the activities plus \(O(\eta_j \log \eta_j)\) per variable with the breakpoint scan of Algorithm 3.2.11 (the listing evaluates \(G_j\) at each of its candidates and costs \(O(\eta_j^2)\)), all variables independent. That independence is the whole point for a device. The sequential algorithm updates one variable per step and keeps the scores of its neighbours current lazily, which is cheap per step but serial. A GPU variant recomputes every violation (one thread per row) and every jump value and score (one thread per column) after each move, or after a batch of moves on variables that share no row. It trades \(O(\mathrm{nnz}(A))\) work per step for full parallelism. The weight update at a local minimum is a reduction over the violated rows. Independent runs from different starting points and seeds are a trivial portfolio on any hardware.
One step of the GPU variant of Feasibility Jump
a move, or a batch of moves on variables that share no row
|
v
rows [ 1 ][ 2 ][ 3 ] . . . [ m ] one thread per row:
| every violation
v
columns [ 1 ][ 2 ][ 3 ] . . . [ n ] one thread per column:
| every jump value, score
v
the next move; at a local minimum, the weight update is a
reduction over the violated rows
O(nnz(A)) work per step, for full parallelism; the sequential
algorithm moves one variable per step and rescores its
neighbours lazily: cheap per step, but serial
(Relatives, and the fused GPU heuristics) Feasibility Jump has relatives. Local-MIP, by Lin, Cai, Zou and Lin, is a stand-alone local-search solver for MIP with its own move operators and weighting scheme. ViolationLS, by Davies, Didier and Perron, is the constraint-based local search inside Google's CP-SAT, which works on the violations of general constraints rather than linear rows.P. Lin, S. Cai, M. Zou and J. Lin, "Local-MIP: efficient local search for mixed integer programming", Artificial Intelligence 348 (2025), 104405 (a conference version appeared at CP 2024); T. O. Davies, F. Didier and L. Perron, "ViolationLS: constraint-based local search in CP-SAT", CPAIOR 2024, LNCS, 243–258. On the GPU the heuristics have been fused. Çördük, Sielski, Boucher and Aatish, the cuOpt team, run PDLP, the GPU first-order LP method of Section 7.2, as an approximate LP solver. They keep a probing cache, a store of the bound implications of fixing single variables, for early infeasibility detection. On the device they combine a feasibility pump, Feasibility Jump and fix-and-propagate, the heuristic of Section 3.2 that fixes a variable, propagates bounds and repeats. They report "221 feasible solutions and 22% objective gap in the MIPLIB2017 benchmark on a presolved dataset".A. Çördük, P. Sielski, A. Boucher and K. Aatish, "GPU-accelerated primal heuristics for mixed integer programming", arXiv 2510.20499 (2025). The 2026 MIP Workshop's computational competition was on exactly this, GPU primal heuristics that call no MIP solver. It was won by CHAP, a hybrid of a GPU tabu search, a local search that forbids recently reversed moves, with cuPDLPx, a GPU implementation of restarted PDHG (Section 7.2), as its LP and fix-and-propagate on the CPU. CHAP found solutions on 47 of the 50 hidden instances within five minutes, against 43 for cuOpt in heuristics-only mode and 44 for Gurobi in default mode.G. K. Tjusila, A. Hoen, N.-C. Kempke, G. Mexi, T. Berthold, A. Gleixner, T. Koch and S. Pokutta, "CHAP: a hybrid GPU-CPU heuristic for MIP", arXiv 2605.05086 (2026); MIP Workshop 2026, computational competition page, mixedinteger.org/2026/competition/. What these heuristics produce is an incumbent, never a bound. Section 7.7 gives the GPU engineering and the measured primal integrals, and Section 7.8 asks what the nonlinear analogue would be, a local search whose violation function is evaluated through the expression graph.
(The second family: separable nonconvexity) The second family needs a factorization but no tree, and it is the family the tax problem of Section 9 belongs to. Its functions are allowed to take the value \(+\infty\), and one word of vocabulary is needed for them. A function \(f : \mathbb{R} \to \mathbb{R} \cup \{+\infty\}\) is closed if it is lower semicontinuous, that is, if its value at a point is at most the limit of its values along any sequence approaching the point; equivalently, if its epigraph, the set of points on or above its graph (Definition 1.2.1), is a closed set. Its domain is the set where it is finite. The value \(+\infty\) is how a constraint on one variable is written as part of its cost, so that a lot that must be sold whole or not at all, or an integer number of shares, is a function and not a constraint.
Definition 5.4.3 (separable–affine problem). A separable–affine problem is
\[\min_{x \in \mathbb{R}^n}\ \sum_{i=1}^n f_i(x_i) \quad \text{subject to} \quad A x = b,\]with \(A \in \mathbb{R}^{m \times n}\) and each \(f_i : \mathbb{R} \to \mathbb{R} \cup \{+\infty\}\) closed, with a finite union of intervals as its domain. Any constraint on a single variable (a bound, an integrality requirement, a minimum trade size) is written as an infinite value of \(f_i\), so that \(m\) counts only the constraints that couple variables.N. Moehle, J. Gindi, S. Boyd and M. J. Kochenderfer, "Portfolio construction as linearly constrained separable optimization", Optimization and Engineering 24 (2023), 1667–1687; arXiv 2103.05455. Their equation (5).
The portfolio problem with a factor risk model \(V = F F^\top + D\) (Section 4.4) is a separable–affine problem in the holdings, the cash and the \(k\) factor exposures \(y = F^\top h\). The only coupling rows are the \(k\) exposure equalities and the budget, so \(m = k + 1\), and every nonconvex term (the tax liability of each asset, a fixed charge, an integer number of shares) sits in its own \(f_i\).
(ADMM: the split and the three steps) The alternating direction method of multipliers, ADMM, is the standard way to split such a problem into a step on the \(f_i\) and a step on the coupling rows. Introduce a copy \(z\) of \(x\) that carries the coupling, and write the problem as
\[\min_{x,\, z}\ \sum_{i=1}^n f_i(x_i) + \mathbb{I}_{\{Az = b\}}(z) \quad \text{subject to} \quad x = z,\]where \(\mathbb{I}_C\) is the indicator of the set \(C\), zero on \(C\) and \(+\infty\) off it. With a penalty \(\rho > 0\) and a scaled multiplier \(u \in \mathbb{R}^n\) for the constraint \(x = z\) (the letter \(\lambda\) is kept for the multipliers of the rows \(Ax = b\) below), the augmented Lagrangian, the Lagrangian of the constraint \(x = z\) with a quadratic penalty on its violation added (Section 1.3 met it in the local methods of LANCELOT and ALGENCAN), is
\[L_\rho(x, z, u) \;=\; \sum_{i=1}^n f_i(x_i) + \mathbb{I}_{\{Az = b\}}(z) + \frac{\rho}{2}\,\lVert x - z + u \rVert^2 - \frac{\rho}{2}\,\lVert u \rVert^2 .\]ADMM minimizes \(L_\rho\) over \(x\) with \(z\) and \(u\) fixed, then over \(z\) with \(x\) and \(u\) fixed, and then takes the multiplier step \(u \leftarrow u + x - z\). The multiplier step is dual ascent on the constraint \(x = z\), whose gradient with respect to the multiplier is the residual \(x - z\). The \(x\)-step separates over the coordinates. For a scalar function \(f\) the proximal point of \(v\) is
\[\operatorname{prox}_{f/\rho}(v) \;=\; \arg\min_{x}\ \Big[ f(x) + \frac{\rho}{2}\,(x - v)^2 \Big],\]the point that trades a decrease of \(f\) against the squared distance from \(v\). Section 7.2 uses the same operator inside the first-order LP method. Step 1 of the algorithm below is \(x_i \leftarrow \operatorname{prox}_{f_i/\rho}(z_i - u_i)\), one scalar minimization per variable, and for a nonconvex \(f_i\) it is computed exactly, piece by piece. The \(z\)-step is the Euclidean projection of \(x + u\) onto the affine set \(\{Az = b\}\), one linear solve with a matrix that never changes.
</figure>
Algorithm 5.4.4 ADMM for a separable-affine problem
(Moehle, Gindi, Boyd and Kochenderfer 2023)
Input closed f_1, ..., f_n on R, each with a finite union of intervals
as domain; A in R^{m x n}; b in R^m;
a penalty rho > 0; tolerances eps_res and eps_obj; a patience N;
a starting pair (z^0, u^0), by default the solution of the
convexified problem in which every f_i is replaced by its
convex envelope.
Output a point x with A x = b and dist(x, dom f) < eps_res, and the
value z_admm = f(Proj_{dom f}(x)).
0. Factor the KKT matrix [[I, A^T], [A, 0]] once and cache the
factorization.
Repeat for k = 0, 1, 2, ...
1. x-update, one proximal step per variable:
x_i^{k+1} = argmin_x f_i(x) + (rho / 2) (x - z_i^k + u_i^k)^2 .
For a piecewise-quadratic f_i: add the quadratic to every piece,
clip each piece's vertex to its interval, take the least value over
the pieces.
2. z-update, the coupling step:
z^{k+1} = argmin { || z - x^{k+1} - u^k ||^2 : A z = b },
one solve with the cached factorization.
3. u^{k+1} = u^k + x^{k+1} - z^{k+1}.
4. Every 10 iterations take the candidate x = z^{k+1} (it satisfies
A x = b); measure r(x) = dist(x, dom f) and o(x) = f(Proj_{dom f}(x));
record (x, o(x)) when r(x) < eps_res and o(x) improves on the best so
far; stop when the best value has not improved by eps_obj in N
iterations.
Invariant
A z^k = b at every iteration. For convex f_i the iterates converge
to a primal-dual optimal pair (Boyd et al. 2011). For nonconvex f_i
there is no guarantee, so the stopping rule is on the best feasible
objective seen, not on residuals.
Cost
step 1 is O(sum_i pieces_i) scalar work with no communication;
step 2 is one back-substitution with the cached factor; step 3 is
O(n). With a factor model the factor is a dense (k + 1) x (k + 1)
Schur complement (Definition 4.7.1; the block-elimination form is
Proposition 7.5.2) after eliminating the diagonal block, so step 2
costs O(n k) per iteration.
Parallel
step 1 is a kernel over variables, and over accounts; step 2 is a
small dense solve plus two products with the n x k exposure matrix,
a batched matrix product across accounts that share one risk model;
the test in step 4 is a reduction.
The split x = z of Algorithm 5.4.4, and what crosses it
sum_i f_i(x_i): separable I_{Az = b}(z): the coupling
+--------------------------+ x + u +--------------------------+
| 1. x_i <- prox_{f_i/rho} | ------> | 2. z <- the projection |
| (z_i - u_i) | | of x + u onto |
| for i = 1, ..., n: | z - u | {A z = b}: one solve |
| one scalar problem | <------ | with the KKT factor |
| each, exact piece by | | cached in step 0 |
| piece, in parallel | | |
+--------------------------+ +--------------------------+
\ /
+------> 3. u <- u + x - z <-----+
dual ascent on x = z
4. every 10 iterations the candidate x = z^{k+1}, which
satisfies A x = b, is measured and the best one recorded
(What the heuristic achieves on portfolio instances) Moehle, Gindi, Boyd and Kochenderfer apply this to portfolio construction with separable nonconvex terms (capital-gains taxes by lot, integer shares, minimum trade sizes) with \(m = k + 1\) coupling rows. On monthly instances with about 1,000 securities and 72 risk factors they report a mean solve time of 251 milliseconds (standard deviation 159) for the nonconvex problem and 152 (standard deviation 67) for the convexified one. The gap \(z_{\mathrm{admm}} - d^\star\) between the heuristic's value and the certified bound of the convexified problem was 0 to 10 basis points of account value (a basis point is a hundredth of a percent) over 692 monthly instances, with mean 0.6 and standard deviation 1.1.Moehle, Gindi, Boyd and Kochenderfer (2023), Sections 5 to 7; the abstract rounds \(n = 998\) and \(k = 72\) to "around 1000 securities and 100 risk factors". The same group's two-stage method for the tax-aware problem of Section 9 solves the convexified problem, rounds the direction of every trade and re-solves the convex restricted problem. It matched the MIQP optimum found by CPLEX 12.9 to within 0.05 basis points on 549 of the 565 instances (of 744) that CPLEX finished within 300 seconds, and was "several hundred times faster".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), 364–383, Section 6; the paper gives the instance count as 744 in its text and as 720 in one figure caption, and this series uses 744. The method descends from the ADMM heuristic for embedded mixed-integer QP of Takapoui, Moehle, Boyd and Bemporad and from the NCVX system of Diamond, Takapoui and Boyd. Both keep every step of ADMM for a convex problem and replace one projection onto a convex set by a projection onto a nonconvex one.R. Takapoui, N. Moehle, S. Boyd and A. Bemporad, "A simple effective heuristic for embedded mixed-integer quadratic programming", International Journal of Control 93 (2020); S. Diamond, R. Takapoui and S. Boyd, "A general system for heuristic minimization of convex functions over non-convex sets", Optimization Methods and Software 33 (2018); S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein, "Distributed optimization and statistical learning via the alternating direction method of multipliers", Foundations and Trends in Machine Learning 3 (2011), for the convex convergence theory. Geißler, Morsi, Schewe and Schmidt showed that the feasibility pump of Section 3.2 is itself a penalty alternating-direction method, so the pump and these heuristics are one family. All of them are block coordinate descent on a nonconvex penalty, that is, they minimize over one block of variables at a time with the others fixed. Such a method converges to partial minima, points that no single block can improve, and therefore needs restarts.B. Geißler, A. Morsi, L. Schewe and M. Schmidt, "Penalty alternating direction methods for mixed-integer optimization: a new view on feasibility pumps", SIAM Journal on Optimization 27 (2017).
(The certificate: the convexified problem's value) What makes the ADMM heuristic more than a heuristic is the bound that comes with it. Run the same algorithm on the convexified problem, in which every \(f_i\) is replaced by its convex envelope \(\hat f_i = \operatorname{vex}_{\operatorname{conv}(S_i)} f_i\) on the interval \(\operatorname{conv}(S_i)\), where \(S_i\) is the domain of \(f_i\) (Section 2.4). That problem is convex, ADMM converges on it, and its value \(\hat z\) is a lower bound on \(z^\star\) because \(\hat f_i \le f_i\): it is a relaxation in the sense of Definition 1.1.3, with the same feasible set and a smaller objective. The same value is the Lagrangian dual value of the original problem, and that is what certifies it. Write \(a_i\) for the \(i\)-th column of \(A\) and \(f^*(s) = \sup_x [\,s x - f(x)\,]\) for the conjugate of a scalar function, as defined in Section 2.2 before Proposition 2.2.3. With a multiplier \(\lambda \in \mathbb{R}^m\) on the rows \(Ax = b\), the dual function of Section 2.2 is
\[\begin{aligned} q(\lambda) \;&=\; \inf_{x} \Big[ \sum_{i=1}^n f_i(x_i) + \lambda^\top (A x - b) \Big] \;=\; -\,b^\top \lambda + \sum_{i=1}^n \inf_{x_i} \big[ f_i(x_i) + (a_i^\top \lambda)\, x_i \big] \\ &=\; -\,b^\top \lambda \;-\; \sum_{i=1}^n f_i^{*}(-a_i^\top \lambda) . \end{aligned}\]The infimum separates over the coordinates because the objective is a sum and the coupling is linear. The next proposition says that the same expression is the dual function of the convexified problem, and that strong duality holds for that problem.
Proposition 5.4.5 (the convexified problem is the Lagrangian dual). Let each \(f_i\) be closed with compact domain \(S_i \subset \mathbb{R}\) and bounded on it, and let \(\hat f_i\) be its convex envelope on \(\operatorname{conv}(S_i)\), extended by \(+\infty\) outside. (i) \(\hat f_i = f_i^{**}\), the biconjugate of \(f_i\), and \(\hat f_i^{\,*} = f_i^{*}\). (ii) The dual function \(q\) above is also the dual function of the convexified problem \(\hat z = \min\{\sum_i \hat f_i(x_i) : Ax = b\}\). (iii) If the convexified problem is feasible, then \(d^\star = \sup_\lambda q(\lambda) = \hat z \le z^\star\). The same holds with inequality rows \(Ax \le b\) and multipliers \(\lambda \ge 0\).For (i): that the biconjugate is the closed convex envelope is stated in Section 2.2 before Proposition 2.2.3 and used in its clause (iii); that for a closed function with compact domain the convex hull is already closed is R. T. Rockafellar, Convex Analysis (Princeton University Press, 1970), Section 12. (ii) and (iii) are the one-variable case of Theorem 3.7.3, whose proof is in Section 3.7.
Proof sketch. (i) is the statement of Section 2.2 that \(f^{**}\) is the largest closed convex function below \(f\), its closed convex envelope, which on a compact domain is the convex envelope itself, together with \(f^{***} = f^{*}\). (ii) By (i), \(\inf_x [\hat f_i(x) + s x] = -\hat f_i^{\,*}(-s) = -f_i^{*}(-s)\), so replacing \(f_i\) by \(\hat f_i\) leaves \(q\) unchanged. (iii) The convexified problem has a closed convex objective, a compact convex domain and affine constraints, so its dual value equals its optimal value (Theorem 3.7.3, where the argument goes through the perturbation function). Weak duality, Theorem 2.2.2, gives \(d^\star \le z^\star\). ∎
The ADMM heuristic with its certificate (Moehle and co-authors)
f_i ------- convex envelope -------> hat f_i
| |
| v
| ADMM on the convexified problem converges;
| with each f_i closed, of compact domain and
| bounded on it, and the convexified problem
| feasible, its value hat z = d* is a lower
| bound on z* (Proposition 5.4.5)
| |
| the starting pair (z^0, u^0) |
| <-------------------------------------+
v
Algorithm 5.4.4 on the nonconvex problem: the heuristic's
value z_admm
reported by Moehle, Gindi, Boyd and Kochenderfer (2023): the
gap z_admm - d* was 0 to 10 basis points of account value over
692 monthly instances, mean 0.6, standard deviation 1.1; mean
solve time 152 ms convexified, 251 ms nonconvex
(Three statements on the size of the gap) The next three statements say how far above \(d^\star\) the optimum can be, and why the answer does not depend on the number of variables. Section 3.7 already used the first of them, the nonconvexity \(\rho(f_i)\) of a block, and pointed here for its definition. For a scalar closed \(f\) with compact domain \(S\) it reads as follows.
Definition 5.4.6 (nonconvexity of a function). Let \(f : \mathbb{R} \to \mathbb{R} \cup \{+\infty\}\) be closed with compact domain \(S\), so that \(f = +\infty\) off \(S\), and let \(\hat f\) be its convex envelope on \(\operatorname{conv}(S)\), which is finite there. The nonconvexity of \(f\) is
\[\rho(f) \;=\; \sup_{x \in \operatorname{conv}(S)} \big( f(x) - \hat f(x) \big) \;\in\; [0, +\infty] .\]\(\rho(f) = 0\) exactly when \(f\) is convex and closed. \(\rho(f) = +\infty\) whenever \(S\) is not an interval, because \(f = +\infty\) at a point of \(\operatorname{conv}(S) \setminus S\) where \(\hat f\) is finite.M. Udell and S. Boyd, "Bounding duality gap for separable problems with linear constraints", Computational Optimization and Applications 64 (2016), 355–378, Section 1, whose functions are stated on convex domains; the remark on nonconvex domains follows from the definition with \(f = +\infty\) off \(S\).
Two instances fix the definition. For the fixed charge of Example 5.4.9 below, \(f = 0\) at \(x = 0\) and \(1 + x\) on \((0, 1]\), the envelope on \([0, 1]\) is the chord \(2x\) and \(\rho(f) = \sup_{x \in (0, 1]} (1 + x - 2x) = 1\), approached as \(x \to 0^+\). For a binary variable, \(S = \{0, 1\}\) with any finite costs at the two points, the envelope is finite at \(\tfrac12\) and \(f\) is not, so \(\rho(f) = +\infty\): the definition charges an integrality requirement an infinite nonconvexity, and the theorem below inherits that.
Theorem 5.4.7 (Shapley–Folkman; Starr 1969). Let \(S_1, \dots, S_n \subset \mathbb{R}^m\) be compact. Every point of \(\operatorname{conv}(S_1 + \cdots + S_n)\) can be written as \(x_1 + \cdots + x_n\) with \(x_i \in \operatorname{conv}(S_i)\) for every \(i\) and \(x_i \in S_i\) for all but at most \(m\) indices.R. M. Starr, "Quasi-equilibria in markets with non-convex preferences", Econometrica 37 (1969), 25–38, where the lemma is attributed to Shapley and Folkman; proof also in I. Ekeland and R. Témam, Convex Analysis and Variational Problems, SIAM Classics 28 (1999), Appendix I.
(What the theorem says) The theorem says that a sum of many nonconvex sets in a space of fixed dimension is nearly convex. The convex hull of the sum differs from the sum itself in at most \(m\) of the summands, however many summands there are. The smallest picture is \(S_i = \{0, 1\}\) on the line, with \(m = 1\). The sum of \(n\) copies is \(\{0, 1, \dots, n\}\) and its convex hull is \([0, n]\). Any \(c \in [0, n]\) is a sum of \(n\) terms from \([0, 1]\) with all but one term in \(\{0, 1\}\): take \(\lfloor c \rfloor\) ones, one term equal to \(c - \lfloor c \rfloor\), and zeros for the rest. Example 5.4.9 below is this picture with a cost attached to each term, the epigraphs of the costs playing the role of the sets and the sum constraint being the one coupling row. Applied to the epigraphs of the \(f_i\) the theorem gives the following.
Theorem 5.4.8 (duality gap of separable problems; Udell and Boyd 2016). Consider the feasible problem
\[z^\star \;=\; \min\Big\{ \sum_{i=1}^n f_i(x_i) \ :\ A x \le b,\ \ G x = h \Big\},\]with each \(f_i\) closed with compact domain \(S_i \subset \mathbb{R}\) and bounded on it, and every constraint on a single variable absorbed into \(S_i\). Let \(\tilde m\) be the number of equality rows plus the largest number of inequality rows that can be active at one point. Let \(\hat z\) be the value of the convexified problem in which each \(f_i\) is replaced by \(\hat f_i\), and order the nonconvexities of Definition 5.4.6 as \(\rho_{(1)} \ge \rho_{(2)} \ge \cdots \ge \rho_{(n)}\), infinite values first. Then the convexified problem has a solution \(x^\star\) with \(x^\star_i \in S_i\) and \(f_i(x^\star_i) = \hat f_i(x^\star_i)\) for all but at most \(\tilde m\) indices, and
\[\hat z \;\le\; z^\star \;\le\; \sum_{i=1}^n f_i(x^\star_i) \;\le\; \hat z + \sum_{i=1}^{\min(\tilde m,\, n)} \rho_{(i)} ,\]an inequality in \((-\infty, +\infty]\). Since \(\hat z = d^\star\) by Proposition 5.4.5, the duality gap satisfies \(z^\star - d^\star \le \sum_{i=1}^{\min(\tilde m, n)} \rho_{(i)}\). The right-hand side is finite only when every \(S_i\) is an interval. The bound is attained on some instances. Udell and Boyd write \(p^\star\) and \(\hat p\) for \(z^\star\) and \(\hat z\).Udell and Boyd (2016), Theorems 1 and 2 and the display after Theorem 1, which are stated for convex domains \(S_i\). The statement with \(f_i = +\infty\) off a possibly nonconvex \(S_i\) is the same inequality read in \((-\infty, +\infty]\); the structural clause about the indices is what the proof gives in either case.
Proof sketch. Write the convexified problem in epigraph form: minimize \(\sum_i t_i\) over \((x_i, t_i) \in K_i := \operatorname{conv} \operatorname{epi}(f_i)\), with each epigraph capped above by \(\max_{S_i} f_i\) so that \(K_i\) is compact, subject to the coupling rows. A linear function over a compact convex set attains its minimum at an extreme point \((x^\star, t^\star)\) of the feasible set, which is \(\prod_i K_i\) intersected with the rows. Let \(L\) be the affine set defined by the equality rows and by the inequality rows active at \((x^\star, t^\star)\). Its codimension is at most \(\tilde m\). Suppose \(\tilde m + 1\) components \((x^\star_i, t^\star_i)\) were not extreme in their \(K_i\). Each such component is the midpoint of a segment inside \(K_i\), so it admits a direction \(d_i \ne 0\) in the \((x_i, t_i)\)-plane with \((x^\star_i, t^\star_i) \pm \epsilon d_i \in K_i\) for small \(\epsilon\). Embedded in the full space, with zeros in every other block, these \(\tilde m + 1\) directions span a subspace of dimension \(\tilde m + 1\), because their supports are disjoint. The linear part of \(L\) has codimension at most \(\tilde m\). Two subspaces whose dimensions add up to more than the dimension of the space meet in a nonzero vector \(v\). Moving from \((x^\star, t^\star)\) by \(\pm \epsilon v\) stays in every \(K_i\) by convexity, stays in \(L\), and keeps the inactive rows satisfied for small \(\epsilon\). So \((x^\star, t^\star)\) is the midpoint of two feasible points and is not extreme, a contradiction. Hence at most \(\tilde m\) components are non-extreme. An extreme point of \(\operatorname{conv} \operatorname{epi}(f_i)\) belongs to \(\operatorname{epi}(f_i)\) itself, so on an extreme component \(x^\star_i \in S_i\) and \(t^\star_i \ge f_i(x^\star_i)\), while optimality gives \(t^\star_i = \hat f_i(x^\star_i) \le f_i(x^\star_i)\). Hence \(f_i(x^\star_i) = \hat f_i(x^\star_i)\) there. On the at most \(\tilde m\) other components, \(f_i(x^\star_i) - \hat f_i(x^\star_i) \le \rho(f_i)\), which is \(+\infty\) when \(x^\star_i \notin S_i\), and summing the worst \(\tilde m\) gives the right-hand inequality. The left inequalities hold because \(\hat f_i \le f_i\) and because \(x^\star\) satisfies the original rows. Udell and Boyd exhibit an instance on which the bound holds with equality. They also give a constructive version: minimizing a random linear objective over the optimal set of the convexified problem yields, with probability one, a point satisfying the inequality. ∎
(The proof in a picture, and the earlier bounds) In a picture, the convexified feasible set is a product of convex sets cut by a few affine constraints. A vertex of such a set can be "fractional", off the original nonconvex pieces, in at most as many coordinates as there are cutting constraints. Each fractional coordinate costs at most its nonconvexity. Adding variables adds terms to the sum but never raises the number of fractional coordinates, so the absolute gap is bounded independently of \(n\) and the relative gap goes to zero as \(n\) grows. The earlier bound of Aubin and Ekeland is \((m + 1)\,\rho_{(1)}\), which uses the largest nonconvexity \(m + 1\) times. Bertsekas gave the Lagrangian form, and Bertsekas, Lauer, Sandell and Posbergh used it to show that the relative duality gap of unit-commitment problems vanishes as the number of generating units grows.J.-P. Aubin and I. Ekeland, "Estimates of the duality gap in nonconvex optimization", Mathematics of Operations Research 1 (1976), 225–245; D. P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods (Academic Press, 1982); D. P. Bertsekas, G. S. Lauer, N. R. Sandell and T. A. Posbergh, "Optimal short-term scheduling of large-scale power systems", IEEE Transactions on Automatic Control 28 (1983), 1–11.
Example 5.4.9 (a fixed charge under one coupling row). Let \(f(x) = 0\) at \(x = 0\) and \(f(x) = 1 + x\) on \((0, 1]\): a fixed charge of one plus a linear cost, the shape of a transaction cost or of a lot that must be sold whole or not at all. Its domain \([0, 1]\) is an interval and \(f\) is bounded on it, so Theorem 5.4.8 applies with a finite right-hand side. Its convex envelope on \([0, 1]\) is the chord \(2x\) through \((0, 0)\) and \((1, 2)\), and its nonconvexity is \(\rho = \sup_{x \in (0, 1]} (1 + x - 2x) = 1\), approached as \(x \to 0^+\). Minimize \(\sum_{i=1}^n f(x_i)\) subject to \(\sum_i x_i = c\) on \([0, 1]^n\). There is one coupling equality, so \(\tilde m = 1\) and Theorem 5.4.8 bounds the gap by \(\rho_{(1)} = 1\) for every \(n\). With \(n = 3\) and \(c = 2.5\) the true optimum must use all three variables and costs \(3 + 2.5 = 5.5\), the envelope relaxation costs \(2c = 5.0\), and the gap is \(0.5\). The gap approaches 1 as \(c\) approaches an integer from above (at \(c = 7.01\) with \(n = 10\) it is \(0.99\)) and vanishes as \(c\) approaches an integer from below. The \((m + 1)\,\rho_{(1)}\) count of Aubin and Ekeland gives 2 for this example, and Theorem 5.4.8 gives 1. Both hold.
# The Shapley-Folkman bound on a separable example (Example 5.4.9).
#
# f(x) = 0 at x = 0 and 1 + x on (0, 1], with convex envelope 2x;
# minimize sum_i f(x_i) subject to sum_i x_i = c on [0, 1]^n.
import itertools
import numpy as np
def f(x):
return 0.0 if x <= 0 else 1.0 + x
def env(x):
return 2.0 * x
grid = np.linspace(0, 1, 2001)[1:]
rho = max(f(t) - env(t) for t in grid)
print(f"nonconvexity rho = sup (f - envelope) = {rho:.4f}, "
"approached as x -> 0+")
def optimum(n, c):
"""The true optimum, using the fewest blocks that can carry c.
Every nonzero block costs 1 + x_i.
"""
k = int(np.ceil(c - 1e-12))
return k + c if k <= n else np.inf
for n, c in [(3, 2.5), (10, 2.5), (10, 7.01), (10, 7.99)]:
p_star, p_hat = optimum(n, c), 2.0 * c
print(f"n = {n:2d}, c = {c:4.2f}: optimum {p_star:5.2f}, "
f"envelope relaxation {p_hat:5.2f}, gap {p_star - p_hat:.2f}")
print(" <= rho_(1) = 1 (one coupling equality)"
" <= (m + 1) rho = 2")
# brute force for n = 3, c = 2.5 on a grid of step 0.05
best = np.inf
for x1, x2 in itertools.product(np.arange(0, 1.0001, 0.05), repeat=2):
x3 = 2.5 - x1 - x2
if -1e-9 <= x3 <= 1 + 1e-9:
best = min(best, f(x1) + f(x2) + f(x3))
print(f"brute force, n = 3, c = 2.5: best cost on the grid {best:.2f}")
nonconvexity rho = sup (f - envelope) = 0.9995, approached as x -> 0+
n = 3, c = 2.50: optimum 5.50, envelope relaxation 5.00, gap 0.50
<= rho_(1) = 1 (one coupling equality) <= (m + 1) rho = 2
n = 10, c = 2.50: optimum 5.50, envelope relaxation 5.00, gap 0.50
<= rho_(1) = 1 (one coupling equality) <= (m + 1) rho = 2
n = 10, c = 7.01: optimum 15.01, envelope relaxation 14.02, gap 0.99
<= rho_(1) = 1 (one coupling equality) <= (m + 1) rho = 2
n = 10, c = 7.99: optimum 15.99, envelope relaxation 15.98, gap 0.01
<= rho_(1) = 1 (one coupling equality) <= (m + 1) rho = 2
brute force, n = 3, c = 2.5: best cost on the grid 5.50
The block costs nothing to run. Its brute-force check is a sum over a grid and would vectorize trivially, but the point is the arithmetic, not the speed.
</figure>
(What the bound says for the tax problem) For the tax problem the count is \(\tilde m = k + 1\): 73 coupling rows for a 72-factor model. Whether the a priori bound says anything then depends on the blocks. When every block is a continuous holding on a direction-split box, each \(S_i\) is an interval and the bound sums the 73 largest kinks of the liability functions. It is loose against the measured gaps of 0 to 10 basis points quoted above. When a block must hold a whole lot or an integer number of shares, its \(S_i\) is not an interval and its nonconvexity is infinite. The bound is then vacuous as soon as the convexified optimum places such a block off its domain. What survives in every case is the structural statement of the proof: an extreme optimal solution of the convexified problem lies off the original domain, or off the original cost, in at most \(\tilde m\) blocks. The bound's value is therefore structural rather than numerical. It says that the gap is governed by the number of risk factors and the size of the tax kinks, not by the number of assets or lots. The tree Section 9 designs therefore branches on a few directions, not on \(2^n\) patterns. One more caveat belongs with the theorem. Not every solution of the convexified problem is near-optimal. Udell and Boyd give a symmetric instance in which ADMM returns the symmetric, non-extreme solution, off by an amount proportional to \(n\), while an extreme point satisfies the bound. They prescribe a second solve with a random linear objective over the optimal set to reach one. Moehle and co-authors avoid the issue by using the convexified solution only as the starting point of Algorithm 5.4.4.Udell and Boyd (2016), Sections 2 and 7; Moehle, Gindi, Boyd and Kochenderfer (2023), Section 6.
(What parallelizes, and what does not) What parallelizes in this subsection is almost everything, which is why its methods were the first to reach the GPU. In Feasibility Jump the jump values and scores are one task per variable and the violations one task per row. In ADMM the proximal step is one task per variable and the coupling step a small dense solve. Across accounts, scenarios or seeds every run is independent. What does not parallelize is the thing none of these methods produces: the dual bound. A heuristic's incumbent shrinks the tree, and the envelope bound \(d^\star\) certifies a separable problem to within the sum of a few nonconvexities. Closing the last few basis points still takes the branch and bound of Section 3.5 on the directions that the relaxation left fractional.
Learning inside the tree
Every decision in the tree of Section 3 that is not a bound computation is a rule written by hand and tuned on a test set. The rules decide which variable to branch on, which cuts to keep, which node to take next, which heuristic to run and how often. Since 2016 the question has been whether a model trained on solved instances can make these decisions better. This subsection says what has been tried, what has shipped, and what standard of evidence separates the two. The exactness of branch and bound does not depend on these decisions, only the size of the tree does (Section 3.1). A learned branching rule or cut selector can waste time but cannot return a wrong answer. Nobody learns the one decision that could, the pruning test, which stays a comparison of a certified bound with the incumbent.
The decisions in the tree, and the one that no model makes
open nodes --> which node to take next? ............... [rule]
^ |
| v
| the relaxation: which cuts to keep? ..... [rule]
| which heuristic to run, how often? ...... [rule]
| |
| v
| bound no better than the incumbent? ..... [test]
| | no | yes: pruned
| v
| which variable to branch on? ............ [rule]
| |
+-------------+ the children
[rule] written by hand and tuned on a test set, or learned:
a poor choice changes the size of the tree, never
the answer
[test] the pruning test, a comparison of a certified bound
with the incumbent: nobody learns it
(Learning to branch) The first line of work learned to branch by imitating strong branching, the rule of Section 3.1 that solves both child LPs of every candidate and is the best branching rule by node count and the worst by time. Khalil, Le Bodic, Song, Nemhauser and Dilkina trained, during the solve of one instance, a ranking model on strong-branching scores and then used it in place of strong branching deeper in the same tree.E. B. Khalil, P. Le Bodic, L. Song, G. Nemhauser and B. Dilkina, "Learning to branch in mixed integer programming", Proceedings of the AAAI Conference on Artificial Intelligence 30 (2016). Gasse, Chételat, Ferroni, Charlin and Lodi trained a graph convolutional network, a network whose layers pass information along the edges of a graph, offline on the bipartite graph of variables and constraints, with the LP solution as node features, to imitate strong branching's choice. On the four instance families they trained on, the network came close to strong branching's decisions at a fraction of its cost. Its trees are several times larger than strong branching's (5 to 50 times in their Table 2) and are solved in a fraction of the time, and it beat SCIP's default reliability branching on instances of the same families.M. Gasse, D. Chételat, N. Ferroni, L. Charlin and A. Lodi, "Exact combinatorial optimization with graph convolutional neural networks", Advances in Neural Information Processing Systems 32 (2019); arXiv 1906.01629. Ecole, the library that exposes SCIP's tree as a learning environment, is A. Prouvost, J. Dumouchelle, L. Scavuzzo, M. Gasse, D. Chételat and A. Lodi, "Ecole: a gym-like library for machine learning in combinatorial optimization solvers", NeurIPS 2020 workshop; arXiv 2011.06069. Nair and co-authors at DeepMind added two components. Neural diving is a network that predicts a partial integer assignment and hands the rest to a sub-MIP. Neural branching is an imitation of strong branching whose expert they made affordable by computing the child LP scores for all candidates at once, with a batched variant of the alternating direction method on a GPU. The batched expert "generates 1.4× and 12× more training data in the same running time" on their two largest data sets. The instances ranged from thousands to a million variables.V. Nair, S. Bartunov, F. Gimeno, I. von Glehn, P. Lichocki, I. Lobov, B. O'Donoghue, N. Sonnerat, C. Tjandraatmadja, P. Wang, R. Addanki, T. Hapuarachchi, T. Keck, J. Keeling, P. Kohli, I. Ktena, Y. Li, O. Vinyals and Y. Zwols, "Solving mixed integer programs using neural networks", arXiv 2012.13349 (2020). Scavuzzo and co-authors left imitation behind and learned to branch by reinforcement, that is, from the outcome of the rule's own decisions rather than from an expert's, with a tree-shaped Markov decision process in which a branching decision is credited with the size of the subtree it creates. The learned rule beats imitation exactly where strong branching is itself a poor expert, on problems whose LP bound barely moves when a variable is fixed.L. Scavuzzo, F. Chen, D. Chételat, M. Gasse, A. Lodi, N. Yorke-Smith and K. Aardal, "Learning to branch with tree MDPs", NeurIPS 2022; arXiv 2205.11107.
Khalil et al. (2016): learning to branch during one solve
root
/ \ strong branching scores
o o the candidates, and a
/ \ / \ ranking model is trained
o o o o on its scores
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
/ \ / \ / \ / \ deeper in the same tree:
o o o o o o o o the model, in place of
strong branching
Gasse et al. (2019): a network trained to imitate strong branching
at a node constraints o o o o the bipartite graph
|\ /|\ /|\ /| of variables and
| X | X | X | constraints, with
|/ \|/ \|/ \| the LP solution as
variables o o o o node features
|
v
graph convolutional network
|
v
the branching variable: close to strong branching's
choice, at a fraction of its cost
trained offline on four instance families; on instances of the
same families its trees are 5 to 50 times larger than strong
branching's, solved in a fraction of the time, and it beat
SCIP's default reliability branching
(Learning to select cuts) Cut selection is the second component. The cut loop of Section 3.3 scores candidates by efficacy, parallelism to the objective, integer support (the share of a cut's nonzero coefficients that fall on integer variables) and mutual parallelism, the scores of Definition 3.3.1, with fixed weights and keeps the best. Turner, Koch, Serrano and Winkler showed that the best fixed weights depend on the instance and trained a network to predict them per instance, which is possible because SCIP exposes cut selection as a plugin.M. Turner, T. Koch, F. Serrano and M. Winkler, "Adaptive cut selection in mixed-integer linear programming", Open Journal of Mathematical Optimization 4 (2023); arXiv 2202.10962. Tang, Agrawal and Faenza learned by reinforcement which Gomory cuts to add from the rows of the tableau. Wang and co-authors learned both which cuts to add and how many with a hierarchical sequence model. Deza and Khalil surveyed the line.Y. Tang, S. Agrawal and Y. Faenza, "Reinforcement learning for integer programming: learning to cut", ICML 2020; arXiv 1906.04859; Z. Wang et al., "Learning cut selection for mixed-integer linear programming via hierarchical sequence model", ICLR 2023; arXiv 2302.00244; A. Deza and E. B. Khalil, "Machine learning for cutting planes in integer programming: a survey", IJCAI 2023; arXiv 2302.09166. The survey of Scavuzzo, Aardal, Lodi and Yorke-Smith is the map of the whole field, organized by solver component (primal heuristics, branching, cutting planes, node selection, configuration). Lodi and Zarpellon's earlier survey covers branching alone.L. Scavuzzo, K. Aardal, A. Lodi and N. Yorke-Smith, "Machine learning augmented branch and bound for mixed integer linear programming", Mathematical Programming 217 (2026), 123–166 (online 2024); arXiv 2402.05501; A. Lodi and G. Zarpellon, "On learning and branching: a survey", TOP 25 (2017).
(What has shipped) What has shipped in a production solver is narrower than what has been published, and it has one shape. CPLEX contains the classifier of Section 4.7, trained offline, that decides at the root whether to linearize the products of binaries in a mixed-integer quadratic program or to keep the quadratic. Bonami, Lodi and Zarpellon describe it in the paper cited there, and it is the first deployment of learning in a commercial MIP solver documented in a refereed paper. FICO Xpress has two more: a learned decision on whether to scale a model, and a learned decision on whether to separate local cuts in the tree. The latter is "a practical implementation inside Xpress on a large, diverse set of real-world industry MIPs".T. Berthold and G. Hendel, "Learning to scale mixed-integer programs", Proceedings of the AAAI Conference on Artificial Intelligence 35 (2021); T. Berthold, M. Francobaldi and G. Hendel, "Learning to use local cuts", Mathematical Programming Computation 17 (2025), 437–450. For nonconvex MINLP, Berthold and Geis learned, inside Xpress Global, which of the existing branching rules to use on a given instance. They report a reduction of 8 to 9% in geometric-mean time and over 10% on hard instances. Ghaddar and co-authors did the same for the spatial branching rules of a polynomial solver by offline algorithm selection from instance features, at no cost at solve time.T. Berthold and F. Geis, "Learning to choose branching rules for nonconvex MINLPs", arXiv 2602.09996 (2026); 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). Whether Gurobi uses learned components inside its solver is not documented in any primary source. Its release notes for versions 11 to 13 contain no such statement, and this series makes none. The shape is the same in every verified case: a decision about configuration, trained offline, cheap at solve time, where a wrong prediction costs time rather than correctness. Learned branching in the inner loop, the subject of most of the research, has not shipped.
(The evidence standard) The reason is the evidence standard, and it has three parts. The first is distribution. Every imitation and reinforcement result above is reported on instances drawn from the family the model was trained on, and the gains do not transfer to a new family without retraining. Gasse and co-authors say so, and the survey of Scavuzzo and co-authors treats generalization as the open problem.Gasse et al. (2019) report results family by family and note the dependence on the training distribution; Scavuzzo et al. (2026) list generalization across distributions among the open questions. The second is the baseline. A learned rule must beat the solver's tuned default on time, not on node count, with several seeds, on instances the solver has not seen. The 2026 results that meet that standard point away from large models. Bayramoğlu, Nemhauser and Sahinidis obtain sparse branching models with "fewer than 4% of the parameters of a state-of-the-art graph neural network" that are "faster than the default solver and the GPU-accelerated GNN". Zhang and co-authors evolve CPU-only branching rules with a language model that outperform both SCIP's default and learned GPU policies.S. Bayramoğlu, G. L. Nemhauser and N. V. Sahinidis, "Speeding up mixed-integer programming solvers with sparse learning for branching", arXiv 2604.00094 (2026); C. Zhang, B. Zhang, Z. Xu, H. Chen, X. Lu, S. Fan, Y. Teng and G. Fan, "BiFE: search-efficient discovery of CPU-only branching policies via LLM-based bi-fidelity evolution", arXiv 2609.36735 (2026). The third is the benchmark itself. A model trained on MIPLIB cannot be evaluated on MIPLIB, which is Mittelmann's worry of Section 5.2 in another form. A benchmark on known instances measures tuning as much as algorithms, and only unknown instances measure the solver. For the reader whose problem is solved thousands of times a day on accounts that look alike, the first part of the standard is a point in favour. The distribution is stable and the labels are cheap, which is the setting in which learned configuration is most likely to pay, and nothing published yet says that it does.
(Which component belongs on the device) For a GPU the question is which learned component belongs on the device. Inference for a network at every node loses to a sparse model on the CPU, as the 2026 results show, because the node is small and the launch latency is not. The device-shaped form is the one Nair and co-authors used: a batched computation of scores for all candidates at once. That is the strong-branching expert itself. If strong branching becomes a batched kernel over the children of a frontier, as Section 6.4 and Section 7.4 describe, the student loses its reason to exist, and what remains to learn is the configuration layer, offline. Section 7.8 lists on-device learned scoring as an open problem for that reason.
Exact arithmetic
(What the tolerances are for, and what they permit) A solver's answer is a statement with tolerances attached. Section 3.5 defined the three tolerances, feasibility, integrality and gap, and tabulated the defaults of BARON, SCIP, Couenne and Gurobi. This subsection says what they permit, how an exact solver removes them, and why removing them stops at the linear case. Section 8 rewrites the three tolerances for a problem stated in dollars. The tolerances are not there to absorb floating-point rounding, which is far smaller: adding 0.1 to itself ten times in binary floating point gives 0.9999999999999999, a relative error of \(10^{-16}\). They are there because the simplex method's factorizations, the cut coefficients and the comparisons of bounds are all computed in floating point and need slack to be stable. Without it a solver declares feasible problems infeasible and cycles on degenerate vertices. The problem is what the tolerances permit, which the indicator example shows. An indicator \(z = 10^{-6}\) is within an integrality tolerance of \(10^{-6}\) of zero, so the solver treats the switch as off, yet \(x \le M z\) with \(M = 10^6\) then allows \(x = 1\), one share held by a position the solver reports as closed. With Gurobi's tolerance of \(10^{-5}\) the same \(M\) allows ten, and its documentation gives this very example.Gurobi Optimizer Reference Manual, "Guidelines for numerical issues", with the example \(z = 0.0000099999\) and \(x = 9.999\) for \(M = 10^6\); the remedy the guidelines give is a bound on \(x\) as small as the model allows, not a smaller tolerance. The floating-point rounding error of \(10^{-16}\) quoted above is not the solver tolerance and plays no part in this arithmetic. A feasibility tolerance of \(10^{-6}\), applied after the solver has scaled a row with coefficients of order \(10^7\) shares down to order one, allows a violation of ten shares in the original units. And a Gomory mixed-integer cut computed in floating point from a tableau row is not guaranteed valid, because the row is only approximately an identity. Cook, Dash, Fukasawa and Goycoolea had to show how to make the cut provably valid by computing it in directed rounding and treating the error terms conservatively.W. Cook, S. Dash, R. Fukasawa and M. Goycoolea, "Numerically safe Gomory mixed-integer cuts", INFORMS Journal on Computing 21 (2009), 641–649.
| tolerance | applied to | multiplied by | what it permits |
|---|---|---|---|
| integrality, \(10^{-6}\) | an indicator \(z\) in \(x \le M z\), treated as zero | \(M = 10^6\) | \(x = 1\): one share, held by a position the solver reports as closed |
| integrality, \(10^{-5}\) (Gurobi's; its documentation gives this example) | the same indicator | \(M = 10^6\) | \(x = 10\): ten shares |
| feasibility, \(10^{-6}\) | a row with coefficients of order \(10^7\) shares, after the solver scales it to order one | \(10^7\), scaling back | a violation of ten shares in the original units |
(How an exact solver removes them) An exact solver removes the tolerances without abandoning floating point. The idea rests on two facts that Section 7.3 proves. By weak duality (Theorem 2.2.2) any dual-feasible vector gives a lower bound on an LP's value. A floating-point dual vector is only approximately feasible, and its bound is repaired in one of two ways. The bound shift of Theorem 7.3.1 charges each reduced cost, the cost of a variable net of what the multipliers pay for its column, the worst it can do over the variable's finite interval and evaluates the result with directed rounding, the rounding toward the safe side that Section 2.6 introduced for interval arithmetic. The repaired number is a rigorous bound, at the cost of one matrix–vector product. Project and shift instead moves the vector onto the exact dual-feasible set and evaluates it in rational arithmetic. The design, due to Cook, Koch, Steffy and Wolter, keeps the floating-point LP for guidance and makes every inference that affects correctness safe or exact. On the one-variable LP \(\min\{x : 3x \ge 2\}\) of this section, the one multiplier \(y \ge 0\) is dual feasible when the reduced cost \(1 - 3y\) is nonnegative, that is, when \(y \le 1/3\), and the raw value \(2y\) is a bound only there.
Two repairs of a floating-point dual vector, one guarantee
the dual vector of the floating-point LP:
only approximately dual feasible
/ \
v v
bound shift (Theorem 7.3.1) project and shift
each reduced cost charged the the vector moved onto
worst it can do over the the exact dual-feasible
variable's finite interval, set, evaluated in
evaluated with directed rounding; rational arithmetic
one matrix-vector product
\ /
v v
a rigorous lower bound on the LP value
(weak duality, Theorem 2.2.2)
Theorem 5.6.1 (exact rational mixed-integer programming; Cook, Koch, Steffy and Wolter 2013). Let a MILP have rational data and finite bounds on all variables. Consider a branch and bound in which:
- Each node's LP relaxation is solved in floating point. The floating-point solution is used only to choose the branching variable and the order in which the bounds below are tried.
- A node is pruned by bound only on a lower bound that is valid in exact arithmetic: a bound shift (Theorem 7.3.1) evaluated with directed rounding, a project-and-shift bound evaluated in rational arithmetic, or an exact rational LP solve.
- A node is declared infeasible only on a certificate verified in exact arithmetic, a Farkas ray (a vector \(y \ge 0\) with \(y^\top A = 0\) and \(y^\top b < 0\) for rows \(Ax \le b\), whose existence certifies infeasibility and whose form with variable bounds is Corollary 7.3.4) checked in rationals or passed through the bound shift, or on an exact LP solve.
- The exact LP is solved when the safe bound is too weak to prune while the floating-point value says that pruning should be possible.
- A floating-point solution that looks integral is verified in rational arithmetic. If it fails, the exact LP is solved and the search continues from the exact solution.
- Every candidate incumbent is verified in rational arithmetic before it is accepted.
- Branching splits the range of an integer variable at an integer.
Then the method terminates after finitely many nodes with an exactly optimal rational solution or an exact proof of infeasibility.
Proof sketch. Theorem 3.1.5 gives the invariant of branch and bound: the smallest node bound over the open list is at most \(z^\star\), the incumbent value is at least \(z^\star\), and every feasible point lies in an open node or is no better than the incumbent. The invariant survives each step here because every number it uses is exact or safe. A node pruned by (2) has a certified lower bound that is at least the incumbent value, so it contains no better point. A node pruned by (3) is certified empty. An incumbent accepted by (6) is exactly feasible, so its value bounds \(z^\star\) from above. A node whose floating-point solution looks integral but fails verification is not lost. By (5) its LP is solved exactly. The exact solution is integral and becomes the incumbent after verification, or it is fractional and is branched on, or its exact value prunes the node. Branching by (7) partitions the node's feasible set, and since every integer range is finite the tree is finite (Theorem 3.1.6). The floating-point LP therefore decides only where to branch and which safe bound to try first. An inaccurate LP can enlarge the tree but cannot change the answer. ∎
Theorem 5.6.1: one node, guided in floating point, decided exactly
the node LP, solved in floating point: guidance only (1)
|
+-- looks infeasible --> a Farkas ray checked in rationals or
| through the bound shift, or an exact
| LP solve: the node is infeasible (3)
|
+-- looks prunable ----> safe bound >= incumbent value: pruned
| (2); safe bound too weak: an exact LP
| solve (4)
|
+-- looks integral ----> verified in rationals before it is
| accepted as incumbent (5), (6); if it
| fails, an exact LP solve, and the
| search continues from the exact
| solution (5)
|
+-- otherwise ---------> branch on the variable it chose,
splitting its range at an integer
(1), (7)
an inaccurate LP can enlarge the tree, not change the answer
(Guidance and inference, and the certificates) The proof is short because the design puts all of the floating-point work on the side of guidance and all of the certified work on the side of inference. Eifler and Gleixner revised the framework a decade later with presolve in rational arithmetic, safe Gomory mixed-integer cuts and a repair heuristic that rounds floating-point solutions to exactly feasible ones. They also made reliability branching exact, since strong branching's bound tightenings are inferences and must be safe too. They report the revised code to be 10.7 times faster than the original framework and to solve 2.9 times as many instances within the limit.W. Cook, T. Koch, D. E. Steffy and K. Wolter, "A hybrid branch-and-bound approach for exact rational mixed-integer programming", Mathematical Programming Computation 5 (2013); D. E. Steffy and K. Wolter, "Valid linear programming bounds for exact mixed-integer programming", INFORMS Journal on Computing 25 (2013), which compares the bound shift, project and shift and the exact LP; the exact LP solver is QSopt_ex, D. Applegate, W. Cook, S. Dash and D. Espinoza, "Exact solutions to linear programming problems", Operations Research Letters 35 (2007). L. Eifler and A. Gleixner, "A computational status update for exact rational mixed integer programming", Mathematical Programming 197 (2023), 793–812; the factors 10.7 and 2.9 are from the abstract of the journal version, whereas the abstract of the preprint arXiv 2101.09141 reports 6.6 and 2.8 for an earlier state of the code. L. Eifler and A. Gleixner, "Safe and verified Gomory mixed-integer cuts in a rational mixed-integer program framework", SIAM Journal on Optimization 34 (2024), 742–763. The solver can also write down why its answer is right. The VIPR certificate format of Cheung, Gleixner and Steffy records every bound as a nonnegative combination of constraints, a rounding, or a branching disjunction, and an independent program checks the file. A checker formally verified in the HOL4 proof assistant, a program that checks every step of a mathematical proof, exists through CakeML, a compiler whose own correctness is proved in the same system.K. K. H. Cheung, A. Gleixner and D. E. Steffy, "Verifying integer programming results", IPCO 2017, LNCS 10328; the verified checker is reported in Hojny et al. (2025), Section 3.1.
(The price, measured) SCIP 10.0 ships this as a solving mode. Switched on by exact/enabled, it is restricted to mixed-integer linear programs and uses the GMP, MPFR and Boost libraries for rational and multiprecision arithmetic. It presolves in rational arithmetic through PaPILO, in parallel, to recover part of the overhead. It computes safe node bounds by bound shift before falling back to an exact LP solve, separates safe Gomory cuts, repairs floating-point heuristic solutions, and emits VIPR certificates. The price was measured on the MIPLIB 2017 benchmark set with a two-hour limit and three seeds, 720 instance-seed runs in all. Exact SCIP solved 161 of the 720 runs, a floating-point SCIP restricted to the same features 235, and default SCIP 342. The slowdown is a factor of 2.6 to 3.6 against the reduced floating-point code and 6.8 to 10.8 against the default.Hojny et al. (2025), Section 3.1, in particular 3.1.3 (presolve), 3.1.4 (safe dual bounding), 3.1.5 (safe cuts) and 3.1.9 (the computational study, shifts 1 s for times and 100 for nodes); A. Gleixner, L. Gottwald and A. Hoen, "PaPILO: a parallel presolving library for integer and linear optimization with multiprecision support", INFORMS Journal on Computing 35 (2023). The price of an exact answer is therefore about an order of magnitude in time. The restriction to the linear case is not an engineering gap but a mathematical one, and the proposition after the figure is the reason: a nonlinear problem's optimum need not be a rational number, so there is nothing for rational arithmetic to return.
SCIP 10's exact mode (exact/enabled, for MILP only)
rational presolve tree search
(PaPILO, in ----> safe node bounds by bound shift,
parallel) an exact LP solve as the fallback;
safe Gomory cuts; floating-point
heuristic solutions repaired
|
v
VIPR certificate: every bound a
nonnegative combination of
constraints, a rounding, or a
branching disjunction
|
v
an independent checker (one is
formally verified in HOL4,
through CakeML)
Proposition 5.6.2 (irrational optima). The problem \(\min\{(x^2 - 2)^2 : 0 \le x \le 2\}\) has optimal value \(0\), attained only at \(x = \sqrt 2\). No rational point attains the optimum, so no finite computation in rational arithmetic can return an optimal solution, and the only certifiable statement about a returned point \(\bar x \in \mathbb{Q}\) is that it is \(\varepsilon\)-optimal for some \(\varepsilon > 0\). By contrast \(\min\{x : 3x \ge 2\}\) has the rational optimum \(2/3\), and every feasible and bounded MILP with rational data has a rational optimal solution (Theorem 2.1.4).
Proof. \((x^2 - 2)^2 \ge 0\) with equality if and only if \(x^2 = 2\), and \(\sqrt 2\) is irrational. The second statement is Meyer's theorem, Theorem 2.1.4 in Section 2.1, used already in Proposition 1.5.6. The convex hull of the feasible set of a rational MILP is a rational polyhedron, and a linear function bounded below on it attains its minimum on a rational face, which contains a rational feasible point. ∎
Example 5.6.3 (the two doubles around \(\sqrt 2\)). This example is Proposition 5.6.2 with numbers. The two doubles adjacent to \(\sqrt 2\), 1.4142135623730949 and 1.4142135623730951, both give the floating-point objective \(2 \cdot 10^{-31}\). In exact arithmetic the second, which is the nearest double, has the smaller residual, \(x^2 - 2 = 2.7 \cdot 10^{-16}\) against \(-3.5 \cdot 10^{-16}\) for the first. Neither is the minimizer, and no decimal with 2, 4, 8 or 16 digits makes \(q^2 - 2\) zero. The optimal solutions of nonlinear problems are algebraic numbers in general, and once \(\exp\), \(\log\) or \(\sin\) appear they need not even be that (Section 1.5). Rational arithmetic represents neither. The honest form of exactness for a nonlinear problem is therefore different in kind: a rigorous enclosure. That is an interval \([\underline z, \overline z]\) computed with directed rounding (Section 2.6) that is guaranteed to contain the optimum, together with a point whose value lies in it. That is what interval branch and bound certifies, what the safe bounds of Section 7.3 give for a relaxation, and what MINLPLib's "solved" means when a dual bound and a primal bound agree to \(10^{-6}\). The tolerance \(\varepsilon\) is part of the problem statement for a nonlinear problem, not a defect of the solver. The right way to choose it is in the units of the application, and for a tax problem that means dollars (Section 8). The block below prints both examples; its indicator line quotes the tolerance of \(10^{-5}\), which \(z = 10^{-6}\) is also within.
# Tolerances and exactness: the two examples of the text.
#
# The first part prints the arithmetic of the tolerances: ten 0.1s summed
# in floating point and in rationals, the indicator example, a scaled
# row. The second is Example 5.6.3, min (x^2 - 2)^2 on [0, 2], whose
# minimizer sqrt 2 no double and no decimal attains.
import math
from fractions import Fraction
s = 0.0
for _ in range(10):
s += 0.1
r = sum(Fraction(1, 10) for _ in range(10))
print("0.1 added ten times")
print(f" in binary floating point: {s!r} (equal to 1: {s == 1.0})")
print(f" in rationals: {r} (equal to 1: {r == 1})")
print(f"relative error of the floating-point sum: {abs(s - 1.0):.1e},")
print(" far below any solver tolerance")
print("an indicator z = 1e-6, within an integrality tolerance of 1e-5 "
"of zero,")
print(f" times M = 1e6 allows x = {1e-6 * 1e6:g}: one share")
print("a feasibility tolerance of 1e-6, applied after a row of order 1e7 "
"shares")
print(f" is scaled to order one, allows {1e-6 * 1e7:g} shares")
print("\nExample 5.6.3: min (x^2 - 2)^2 on [0, 2]; the minimizer is sqrt 2")
lo, hi = 1.0, 2.0
# bisection on the sign of x^2 - 2 ends on two adjacent doubles
for _ in range(60):
mid = 0.5 * (lo + hi)
lo, hi = (mid, hi) if mid * mid < 2 else (lo, mid)
print(" the two doubles adjacent to sqrt 2:")
print(f" {lo:.16f} and {hi:.16f}")
print(f" math.sqrt(2) = {math.sqrt(2.0):.16f}")
for x in (lo, hi):
# the residual of the double itself, in exact arithmetic
exact = Fraction(x) ** 2 - 2
print(f" x = {x:.16f}: "
f"floating-point objective {(x * x - 2) ** 2:.1e},")
print(f" exact x^2 - 2 = {float(exact):+.1e}")
for d in (2, 4, 8, 16):
q = Fraction(round(hi * 10 ** d), 10 ** d)
print(f" the decimal with {d:2d} digits, q = {float(q):.16f}:")
print(f" q^2 - 2 = {float(q * q - 2):+.1e}, not zero")
print(f" compare min x subject to 3x >= 2: optimum {Fraction(2, 3)}, "
"exact in rationals")
0.1 added ten times
in binary floating point: 0.9999999999999999 (equal to 1: False)
in rationals: 1 (equal to 1: True)
relative error of the floating-point sum: 1.1e-16,
far below any solver tolerance
an indicator z = 1e-6, within an integrality tolerance of 1e-5 of zero,
times M = 1e6 allows x = 1: one share
a feasibility tolerance of 1e-6, applied after a row of order 1e7 shares
is scaled to order one, allows 10 shares
Example 5.6.3: min (x^2 - 2)^2 on [0, 2]; the minimizer is sqrt 2
the two doubles adjacent to sqrt 2:
1.4142135623730949 and 1.4142135623730951
math.sqrt(2) = 1.4142135623730951
x = 1.4142135623730949: floating-point objective 2.0e-31,
exact x^2 - 2 = -3.5e-16
x = 1.4142135623730951: floating-point objective 2.0e-31,
exact x^2 - 2 = +2.7e-16
the decimal with 2 digits, q = 1.4099999999999999:
q^2 - 2 = -1.2e-02, not zero
the decimal with 4 digits, q = 1.4141999999999999:
q^2 - 2 = -3.8e-05, not zero
the decimal with 8 digits, q = 1.4142135600000001:
q^2 - 2 = -6.7e-09, not zero
the decimal with 16 digits, q = 1.4142135623730951:
q^2 - 2 = +4.3e-16, not zero
compare min x subject to 3x >= 2: optimum 2/3, exact in rationals
The block is sixty bisection steps and a few rational operations. Nothing in it is parallel, and the Fraction type is the whole point, since it is what a GPU does not have.
</figure>
(Where exact mode is used, and what it means for the tax problem) Exact mode is used where the answer is a theorem: a proof in combinatorics that rests on a MILP being infeasible, a verification task, a contract in which a rounding error is a wrong answer rather than an approximation. For the tax problem of Section 9 the lesson is narrower and still useful. The tax part of the model is linear in the lot variables, so a lot-selection subproblem can in principle be solved and certified to the cent in rational arithmetic. The risk term is quadratic and its optimum is a point whose coordinates are, in general, irrational, so the problem as a whole is certified to a tolerance in dollars. And the wash-sale indicator is a fact that a position of one share can get wrong: the position whose indicator the solver reports as zero but which \(x \le M z\) lets hold one share, as in the indicator example above. The three statements are different, and a solver's log does not distinguish them.
Where this is used
Every component of this section runs in a named code. The MILP engine of Section 5.1, with the cut loop, presolve and branching that the ablations credit, is the engine inside BARON (through CPLEX, Xpress or CLP/CBC), SCIP (through SoPlex or CPLEX), Xpress Global, Gurobi and COPT, and the landscape table of Section 5.3 names the stack under every global solver. The benchmark arithmetic of Section 5.2 is the arithmetic of Mittelmann's pages, of the SCIP release reports and of the Xpress Global paper, and the seven rules stated there govern every number quoted in Sections 6 to 9. Feasibility Jump runs by default in Xpress 9 and HiGHS 1.11, on the GPU in cuOpt, and in SCIP's development branch. The ADMM heuristic of Algorithm 5.4.4 runs in the portfolio systems of Moehle and co-authors and is the method Section 9 adapts to the tax problem. Learning has shipped as offline configuration, in CPLEX as the linearization classifier and in Xpress as the decisions on scaling, local cuts and the branching rule, and not as in-loop branching. Exact arithmetic ships in SCIP 10 as a mode for mixed-integer linear programs with VIPR certificates, at a cost of about an order of magnitude in time.
What parallelizes
The components that carry the engine, the cut loop, presolve and branching, are sequential passes over rows and columns or batches of child LPs solved by a simplex code, and every global solver of Section 5.3 runs its tree on the CPU. The heuristics that reached the GPU first are the LP-free ones of Section 5.4. In Feasibility Jump the jump values and scores are one task per variable and the violations one task per row. In ADMM the proximal step is one task per variable and the coupling step a small dense solve. Across accounts, scenarios or seeds every run is independent. None of these methods produces a dual bound. Exact arithmetic splits across host and device. Directed rounding is available per operation on NVIDIA hardware (__dadd_rd, __dmul_ru and their relatives), so the safe bound of Section 7.3 translates to the device as it stands, and a frontier of node bounds can be certified in one kernel. Rational arithmetic does not translate: variable-length integers, data-dependent branching and a few verifications per node are host work. Floating-point guidance in bulk on the device, with safe bounds computed there and the rare exact verification on the host, is the exact-MIP architecture of Theorem 5.6.1 drawn across two processors. Its cost on the CPU is a factor of about ten today, and whether the device changes that factor is one of the questions of Section 7.8. The MILP engine became about fifty times faster, with the hardware held fixed, in twenty years, and the ablation studies agree on which of its components are large and which are small, though not on the order among the large ones. Section 6 asks how much of that engine, and of the tree around it, runs on many processors at once, and what the measured answer has been.
Theorem 5.6.1 drawn across two processors
device host
+--------------------------+ +----------------------+
| floating-point guidance, | | rational arithmetic: |
| in bulk | the rare | variable-length |
| | exact | integers and data- |
| safe bounds by directed | verification | dependent branching |
| rounding per operation | -----------> | |
| (__dadd_rd, __dmul_ru | | a few verifications |
| and their relatives) | | per node |
| | | |
| a frontier of node | | |
| bounds certified in one | | |
| kernel | | |
+--------------------------+ +----------------------+
its cost on the CPU is a factor of about ten today; whether
the device changes that factor is a question of Section 7.8