6. Parallelism
Draft
The node pipeline of Section 3 was written for one processor. On a hard nonconvex instance it ends at a time limit with a gap. The obvious remedy is more processors. This section is about what more processors buy for branch and bound, what limits them, and which patterns the GPU work of Section 7 is built on. The subject has forty years of theory and systems behind it. The theory says when adding a processor can make a search slower, how much idle time no load balancer can remove, and what a batch of nodes evaluated together costs in extra work. We begin with a taxonomy of the levels at which an exact solver can run in parallel. We then state the theorems, describe the systems that hold the records, and end with branch and bound on a GPU, which is the design problem the rest of this series leads to.
Two facts organize everything that follows. The first is that the search tree is the one irregular part of a global solver. Its shape is unknown in advance, it depends on the incumbent found along the way, and two sibling subtrees can differ in size by many orders of magnitude. Everything else in a node is regular work: bounding (an LP, a convex NLP, or an interval evaluation), constraint propagation (one pass per constraint), optimization-based bound tightening (one small LP per variable bound), and primal heuristics (independent trials). The second fact is that parallelizing the tree changes the algorithm. With \(P\) workers the set of nodes that gets evaluated depends on the schedule. The running time is therefore not a function of the instance alone. Speedups above \(P\) and below \(1\) both occur, and node throughput is not the same thing as time to solution. Both facts recur in every subsection.
The first fact: an irregular tree whose nodes do regular work
N the tree: shape unknown in advance,
/ \ dependent on the incumbent found
/ \ along the way
/\ \
/ \ /\ two sibling subtrees: their sizes
/ \ /__\ small can differ by many orders of
/ \ magnitude
/ \
/____o_____\ large
|
v
inside each node, everything else is regular work:
bounding an LP, a convex NLP, or an interval
evaluation
constraint propagation one pass per constraint [ ][ ][ ] ...
optimization-based one small LP per
bound tightening variable bound [ ][ ][ ] ...
primal heuristics independent trials [ ][ ][ ] ...
What parallelizes
This subsection defines the measures of a parallel run and sorts the ways of running an exact solver in parallel into a table. The table is the reference for Section 7, where each row becomes a question about a GPU. The two measures, speedup and efficiency, are used loosely in the literature, and the definition below fixes them, together with a third, the isoefficiency function, so that the theorems of Section 6.2 and the measurements of Section 6.3 can be read on one scale. The caveat at the end of the subsection, that a parallel branch and bound does not explore the tree the sequential run explores, is the reason every speedup for a tree search has to be reported with a node count beside it.
Definition 6.1.1 (speedup, efficiency, isoefficiency). For a fixed instance and fixed settings, let \(T(P)\) be the wall-clock time to the stopping criterion with \(P\) processing elements (threads, processes or devices). The speedup is \(S(P) = T(1)/T(P)\) and the efficiency is \(E(P) = S(P)/P\). With \(W = T(1)\) the work, the overhead of the run is \(T_o(W, P) = P\,T(P) - W\), the processor time not spent on the sequential algorithm's work, and
\[E(P) \;=\; \frac{W}{W + T_o(W, P)} .\]The isoefficiency function \(W_E(P)\) is the work needed to hold the efficiency at \(E\), and it solves \(W = \frac{E}{1 - E}\,T_o(W, P)\). An algorithm is called scalable when \(W_E(P)\) grows slowly in \(P\); the lower bound, and the ideal, is \(\Theta(P)\), and \(P \log P\) is close to it.V. Kumar and V. N. Rao, "Parallel depth first search. Part II. Analysis", International Journal of Parallel Programming 16 (1987); A. Grama, A. Gupta, G. Karypis and V. Kumar, Introduction to Parallel Computing, 2nd edition (Addison-Wesley, 2003), Chapter 5.
The overhead \(T_o\) has four sources: communication, idle time during the ramp-up and ramp-down phases when the tree is too thin to feed every worker, idle time from latency and contention, and redundant work, meaning nodes a parallel run evaluates that the sequential run would have pruned. Every later theorem bounds one of them.T. Koch, T. Ralphs and Y. Shinano, "Could we use a million cores to solve an integer program?", Mathematical Methods of Operations Research 76 (2012). The four categories are theirs. For branch and bound the phrase "the work \(W\)" needs one caveat, which Section 6.2 makes precise. The parallel run and the sequential run explore different trees. So \(W\) is the sequential tree's size, and the parallel tree's size is a separate quantity that must be reported next to the speedup.
Definition 6.1.1: the processor time P T(P), worker by worker
0 T(P)
+---------------------------------------------+
worker 1 |=========c==========---=========c======......|
worker 2 |..=====c========*****========c=====---=======|
worker 3 |....======c===========c===****=======........|
: | |
worker P |.......==c=======---======c==============....|
+---------------------------------------------+
ramp-up ramp-down
= the sequential algorithm's work
c communication
. idle in the ramp-up and the ramp-down: the tree is too thin
to feed every worker
- idle from latency and contention
* redundant work: nodes the sequential run would have pruned
the box is P T(P) = W + T_o(W, P), with W = T(1)
S(P) = T(1)/T(P) E(P) = S(P)/P = W / (W + T_o(W, P))
Since this is the first time the definition is applied, it is worth writing the symbols out on one set of numbers, those of the example that Section 6.2 attaches to Theorem 6.2.9. There a tree of \(W = T(1) = 10^9\) node-units is run on \(P = 10^4\) workers, and the theorem's bound allows a wall-clock time of \(T(P) = 1.2 \times 10^5\) node-units. The speedup is \(S(P) = T(1)/T(P) = 8{,}333\), the efficiency is \(E(P) = S(P)/P = 0.833\), the processor time is \(P\,T(P) = 1.2 \times 10^9\), so the overhead is \(T_o = P\,T(P) - W = 2 \times 10^8\), and the overhead's share of the work is \(T_o/W = 0.2\), which is \((1 - E)/E\) as the definition's identity requires. Holding \(E = 0.833\) on more workers then needs \(W = 5\,T_o\): the work has to grow five times as fast as whatever the overhead grows to.
Definition 6.1.2 (the three types of parallelism; Gendron and Crainic, 1994). Type 1: the operations on one node are performed in parallel, for example the bound computation itself. Type 2: several nodes of one tree are evaluated at once, which changes the search and is where anomalies arise. Type 3: several branch-and-bound trees run concurrently with different parameters and exchange information (today called concurrent or portfolio solving). Within Type 2 the pool of open nodes is centralized (one list) or distributed (one list per worker with load balancing), and node exchange is synchronous or asynchronous.B. Gendron and T. G. Crainic, "Parallel branch-and-bound algorithms: survey and synthesis", Operations Research 42 (1994). Later surveys add the dimensions of knowledge sharing (incumbent, bounds, cuts) and of granularity (node or subtree): T. G. Crainic, B. Le Cun and C. Roucairol, "Parallel branch-and-bound algorithms", in E.-G. Talbi (ed.), Parallel Combinatorial Optimization (Wiley, 2006); D. A. Bader, W. E. Hart and C. A. Phillips, "Parallel algorithm design for branch and bound", in Tutorials on Emerging Methodologies and Applications in Operations Research (Springer, 2005).
Definition 6.1.2: the three types, and the pool of Type 2
Type 1 the operations on one node in N
parallel, e.g. the bound / | \
computation itself w_1 w_2 .. w_P
Type 2 several nodes (*) of one tree o
evaluated at once / \
o o
/ \ / \
* * * *
Type 3 several trees with different o <--> o ...
parameters, exchanging / \ / \
information o o o o
Type 2, the pool of open nodes:
centralized: one list distributed: one list per worker,
with load balancing
+---- load balancing ----+
| | |
+-----------------------+ +-----+ +-----+ +-----+
| L, all open nodes | | L_1 | | L_2 | ... | L_P |
+-----------------------+ +-----+ +-----+ +-----+
| | | | | |
w_1 w_2 ... w_P w_1 w_2 ... w_P
node exchange: synchronous or asynchronous
The anomalies that the definition's Type 2 mentions are the runs in which adding processors makes a search slower, or faster by more than the number of processors added; they are the subject of the first half of Section 6.2, which defines them and says when they can and cannot occur. The taxonomy below extends this classification by one level above the tree and one below the node, because both matter for the problems this series is heading toward. Above the tree sit independent problems: accounts, scenarios, seeds. Below the node sit the steps of the node pipeline of Section 3.6, each with a different shape of work. The column "speedup" says what the theory and the measured results of this section allow one to expect. The column "who does it" names the systems of Section 6.3 and Section 7 that do it today. Two terms in the table are defined later: a work-stealing deque in Definition 6.2.7, and the batch inflation \(N(b)/N(1)\) in Definition 6.4.1. Four were defined earlier and are used here as reminders only: the warm start, a simplex run that begins from the parent node's optimal basis rather than from scratch (Section 3.1); FBBT and OBBT, the feasibility-based and optimization-based bound tightening of Section 2.6, the first a pass over the constraints that shrinks variable bounds by interval arithmetic, the second one small LP per variable bound; the expression DAG, the graph of a factorable function's operations (Section 2.4); and PDHG, the primal-dual hybrid gradient method, a first-order method for linear programs that does matrix-vector products instead of pivots (Section 7.2). The hyper-heuristic named in the second row is cuOpt's controller over its GPU and CPU local searches, a population of candidate solutions with rules for recombining them that its parameters document.
| level | shape of the work | speedup | who does it |
|---|---|---|---|
| across problems | independent instances: accounts, scenarios, seeds | linear | processes, machines, anyone |
| across searches (portfolio) | one problem, several configurations, first to finish wins; racing ramp-up | a portfolio effect at the price of redundant work | Gurobi ConcurrentMIP, SCIP concurrent mode, UG racing ramp-up, cuOpt's hyper-heuristic |
| across tree nodes | irregular, unknown in advance, coupled through the incumbent | sublinear as a rule: ramp-up, ramp-down and anomalies in either direction; randomized work stealing gives an expected time of \(T_1/P + O(T_\infty)\) on the tree actually executed (Theorem 6.2.9) | UG/ParaSCIP and FiberSCIP, Xpress, Gurobi, CPLEX, cuOpt; Chase–Lev deques on one machine |
| across nodes, in a batch | same-shape relaxations evaluated together on a device | high throughput, but nodes bounded for nothing: the batch inflation \(N(b)/N(1)\) | GPU frontier batching (Gmys; the Chapel codes; B3-PWL; Liu and Lodi); batched PDHG (Blin et al.) |
| inside a node: the LP | simplex: sequential pivots, warm-started | about 1 on a device; a few-fold on a multicore CPU | one core; a few cores for a parallel dual simplex (Section 7.1) |
| PDHG: two matrix-vector products per iteration | large, at low accuracy; the bound needs a correction | GPU: cuOpt, Gurobi 13, COPT; Sections 7.2 and 7.3 | |
| inside a node: the NLP | KKT factorization with pivoting | large only once condensed, at \(10^{-4}\) to \(10^{-6}\) accuracy | GPU via condensed KKT: MadNLP with cuDSS; Section 7.5 |
| inside a node: tightening | FBBT: one pass per constraint; OBBT: \(2n\) LPs sharing one matrix | a map over constraints (Jacobi) that may need more sweeps; \(2n\) LPs batch | CPU today; the Jacobi-style GPU sweep and batched OBBT are open problems of Section 7.8 |
| inside a node: the bound by evaluation | interval or McCormick evaluation of the DAG on many sub-boxes or points | embarrassingly parallel over sub-boxes, at a cost \(m^d\) | MAiNGO's GPU interval bounder; SourceCodeMcCormick.jl and ParBB |
| heuristics | independent trials: multistart, local search moves, pump roundings | embarrassingly parallel | cuOpt and CHAP on the GPU; CPU threads in every solver |
One row needs comment before the theory. The row "inside a node: the LP" is the reason the GPU story is not simple. The simplex method's pivots are sequential and warm-started from the parent's basis, which is why the LP at a node has stayed on one core for thirty years. The first-order alternative of Section 7.2 gives up exactness for parallelism. Its bounds therefore have to be repaired before a tree can prune on them, which is the subject of Section 7.3. The other rows are taken up in order: the tree in Sections 6.2 and 6.3, and the batch in Section 6.4.
The row "inside a node: the LP": two ways to bound a node
simplex on one core: sequential pivots, warm-started
parent's +-----+ +-----+ +-----+
basis ---->|pivot|-->|pivot|-->...-->|pivot|--> exact bound --+
+-----+ +-----+ +-----+ |
v
a tree can prune on
the bound
^
PDHG on a GPU (Section 7.2): two matrix-vector |
products per iteration, low accuracy |
+-------+ +-------+ |
|mat-vec|-->|mat-vec|-->...--> inexact bound --> repaired ----+
|mat-vec| |mat-vec| (Section 7.3)
+-------+ +-------+
Measurement needs one caveat. The speedup of Definition 6.1.1 compares the same code at two processor counts, and for a branch and bound the two runs explore different trees. A parallel run that finds a good incumbent early can prune a subtree the sequential run spent an hour in. A parallel run can also expand nodes the sequential run never touched. Both effects are real, and papers that report node counts next to times show both. Gmys goes one step further and initializes his scaling experiments with the optimal value. No incumbent is then found during the run, and the parallel tree is the sequential tree. Helbecque and coauthors follow the same protocol. This measures the machine and the load balancer, not the solver, and the papers say so.J. Gmys, "Exactly solving hard permutation flowshop scheduling problems on peta-scale GPU-accelerated supercomputers", INFORMS Journal on Computing 34 (2022), who measures scaling "in the absence of speed-up anomalies" by initializing with the optimum; G. Helbecque, E. Krishnasamy, T. Carneiro, N. Melab and P. Bouvry, "A Chapel-based multi-GPU branch-and-bound algorithm", Euro-Par 2024 Workshops, LNCS 15385 (2025), follow the same protocol. Any reported speedup for a tree search should state which of the two was measured.
One instance, two runs: the sequential tree and a parallel tree
the sequential run a parallel run
root root
/ \ / \
o \ o x S pruned
/ \ \ / \ at its root
o x /\ o \
/ \ /\
/ S \ / \ R
/______\ /____\
x pruned
S a subtree the sequential run spent an hour in; a parallel run
that finds a good incumbent early can prune it at its root
R nodes the sequential run never touched, below a node it
pruned; a parallel run can expand them
started from the optimal value (Gmys; Helbecque and coauthors),
no incumbent is found during the run: the two trees are the same
Parallel branch and bound: theory
This subsection states what is proved about branch and bound on many processors. The questions are five: when adding processors can slow a search down, what a randomized work-stealing scheduler guarantees, how much idle time the shape of the tree forces, how a parallel search knows that it is finished, and what a deterministic parallel search costs. The theorems are old, mostly from the 1980s and 1990s, and they are stated for abstract trees. That is why they transfer without change from MILP to spatial branch and bound and to a GPU. We minimize throughout. The knapsack scripts maximize and say so.
The synchronous model and its anomalies
Definition 6.2.1 (branch-and-bound tree; essential, critical and non-essential nodes). Let a search produce a rooted tree whose node \(N\) carries a lower bound \(\bar z(N)\) such that \(\bar z(N) \le z^\star\) whenever \(N\) contains an optimal point, and \(\bar z(N') \ge \bar z(N)\) for every child \(N'\) of \(N\) (monotone bounds). A leaf is a node whose relaxation point is feasible for the original problem. Its bound is its objective value, and the incumbent \(z_{\mathrm{inc}}\) is the smallest leaf value found so far. A node is essential if \(\bar z(N) < z^\star\), critical if \(\bar z(N) = z^\star\), and non-essential if \(\bar z(N) > z^\star\). Every correct algorithm expands every essential node it generates, and with monotone bounds every essential node is generated.
The essential nodes are the work no algorithm can avoid, because no incumbent can prune them: an incumbent has value at least \(z^\star\), and an essential node's bound is strictly below that. The critical nodes are the ones a tie-breaking rule decides about, and the non-essential nodes are the ones a good incumbent prunes.
Definition 6.2.2 (the synchronous parallel model; Lai and Sahni). With \(n\) processors, each iteration selects the \(n\) open nodes of smallest bound (fewer if fewer are open), evaluates their children, updates the incumbent, prunes, and repeats. \(I(n)\) is the number of iterations until termination, and with unit time per node it is the running time. The sequential algorithm is the case \(n = 1\).
Lai and Sahni write \(n\) for the processor count, and this block and its script keep their symbol so that \(I(n)\) reads as in the paper. From Theorem 6.2.6 on the processor count is \(P\), as in Definition 6.1.1, and \(n\) returns to its usual meaning of a number of variables or summands.
Definition 6.2.3 (anomalies). For \(n_1 < n_2\), a detrimental anomaly is \(I(n_2) > I(n_1)\): more processors, more time. An acceleration anomaly is \(I(n_1)/I(n_2) > n_2/n_1\): a speedup larger than the processor ratio. Both are properties of a pair of schedules, not of any hardware.
Theorem 6.2.4 (Lai and Sahni, 1984: both anomalies exist). In the synchronous best-first model there are trees for which \(I(n_2) > I(n_1)\) with \(n_1 < n_2\), and trees for which \(I(n_1)/I(n_2) > n_2/n_1\).T.-H. Lai and S. Sahni, "Anomalies in parallel branch-and-bound algorithms", Communications of the ACM 27 (1984). The abstract states both existence results and reports experiments on 0/1 knapsack and travelling salesman instances. The two trees below are this post's constructions, not the paper's.
Proof (by construction). The script below builds both trees. Ties are broken first-in-first-out in both constructions. Detrimental: let the optimum \(z^\star = 5\) lie at depth four below the root, at the end of a chain of nodes with bounds \(1, 2, 3\). Let a second child of the root carry the bound exactly \(5\) and head a bushy subtree whose internal nodes all have bound \(5\) and whose leaves are worth \(6\). One processor expands the root and the three chain nodes in iterations 1 to 4. When the leaf worth \(5\) enters the open list, the root's bound-\(5\) child is older, so iteration 5 expands that child. Iteration 6 selects the leaf, and the run stops at iteration 7 because the smallest open bound is \(5 \ge z_{\mathrm{inc}}\). Two processors expand both children of the root in the second iteration. From then on each iteration takes one chain node and one bound-\(5\) node. When the optimal leaf appears in the open list, the bound-\(5\) nodes already there are older than it, so the FIFO rule selects them first: eight iterations. Acceleration: let the root's older child \(A\) have bound \(5\), let its only child \(A'\) have bound \(5\), and let the only child of \(A'\) be the optimal leaf, worth \(5\). Let the root's younger child head a chain with bounds \(1, 2, 3, 4\) that ends in a node \(Q\) of bound \(5\). Let \(Q\) head a ternary tree of depth two whose internal nodes have bound \(5\) and whose nine leaves are worth \(6\). One processor expands the root and the chain in iterations 1 to 5, then \(A\), \(Q\) and \(A'\) in FIFO order in iterations 6 to 8. It then expands the three children of \(Q\) in iterations 9 to 11, selects the optimal leaf at iteration 12 and stops at iteration 13. Two processors expand \(A\) and the chain's head together in iteration 2, \(A'\) and the next chain node in iteration 3, and the leaf worth \(5\) with the third chain node in iteration 4. Iteration 5 expands the last chain node, whose child \(Q\) has bound \(5 \ge z_{\mathrm{inc}}\) and is never pushed, and the list is empty: five iterations, a speedup of \(2.6\) with two processors. ∎
The two tables list the nodes expanded in each iteration of the two constructions, and the script below is the check.
| iteration | \(n = 1\) | \(n = 2\) |
|---|---|---|
| 1 | root | root |
| 2 | chain node, bound 1 | chain node, bound 1; \(B\) (bound 5) |
| 3 | chain node, bound 2 | chain node, bound 2; a bound-5 node below \(B\) |
| 4 | chain node, bound 3 | chain node, bound 3; a bound-5 node below \(B\) |
| 5 | \(B\) (bound 5), older than the leaf | two bound-5 nodes, both older than the leaf |
| 6 | the leaf worth 5; \(z_{\mathrm{inc}} = 5\) | two bound-5 nodes, both older than the leaf |
| 7 | stop: smallest open bound \(5 \ge z_{\mathrm{inc}}\) | the leaf worth 5; \(z_{\mathrm{inc}} = 5\) (the second node selected is pruned) |
| 8 | stop: smallest open bound \(5 \ge z_{\mathrm{inc}}\) |
The detrimental tree of Theorem 6.2.4: when each node is taken
(i / j): taken in iteration i when n = 1 and j when n = 2
root (1 / 1)
/ \
/ \
chain node, bound 1 (2 / 2) B, bound 5 (5 / 2)
| / | \
chain node, bound 2 (3 / 3) below B: internal nodes of
| bound 5, leaves worth 6
chain node, bound 3 (4 / 4) n = 1: none is taken
| n = 2: one in iteration 3,
leaf worth 5 = z* (6 / 7) one in 4, two in 5, two in
6, one pruned in 7
essential (bound < 5): the root and the three chain nodes
critical (bound = 5): B, internal nodes below B, the leaf worth 5
non-essential (bound > 5): the leaves worth 6
the run stops in iteration 7 when n = 1 and in 8 when n = 2
| iteration | \(n = 1\) | \(n = 2\) |
|---|---|---|
| 1 | root | root |
| 2 | chain node, bound 1 | chain node, bound 1; \(A\) (bound 5) |
| 3 | chain node, bound 2 | chain node, bound 2; \(A'\) (bound 5) |
| 4 | chain node, bound 3 | chain node, bound 3; the leaf worth 5; \(z_{\mathrm{inc}} = 5\) |
| 5 | chain node, bound 4 | chain node, bound 4; its child \(Q\) has bound \(5 \ge z_{\mathrm{inc}}\) and is not pushed; the list is empty: done |
| 6 | \(A\) (bound 5) | |
| 7 | \(Q\) (bound 5) | |
| 8 | \(A'\) (bound 5) | |
| 9 to 11 | the three children of \(Q\), one per iteration | |
| 12 | the leaf worth 5; \(z_{\mathrm{inc}} = 5\) | |
| 13 | stop: smallest open bound \(5 \ge z_{\mathrm{inc}}\) |
The acceleration tree of Theorem 6.2.4: when each node is taken
(i / j): taken in iteration i when n = 1 and j when n = 2
root (1 / 1)
/ \
/ \
A, bound 5 (6 / 2) chain node, bound 1 (2 / 2)
| |
A', bound 5 (8 / 3) chain node, bound 2 (3 / 3)
| |
leaf worth 5 = z* chain node, bound 3 (4 / 4)
(12 / 4) |
chain node, bound 4 (5 / 5)
|
Q, bound 5 (7 / never pushed)
/ | \
Q heads a ternary tree of depth
two: internal nodes of bound 5,
nine leaves worth 6; n = 1 takes
the three children of Q in 9 to 11
n = 2 finds z* = 5 in iteration 4, so Q is never pushed and the
run stops in iteration 5; n = 1 stops in 13. The tree has 30 nodes
(The two mechanisms) In the terms of Definition 6.2.1 the first tree has four essential nodes, the root and the three chain nodes, because their bounds lie below \(z^\star = 5\); the nodes with bound exactly \(5\) are critical, which is why the tie-breaking rule decides how long the run takes; and the leaves worth \(6\) are non-essential, pruned as soon as the incumbent is \(5\). The two mechanisms are the ones to remember. A parallel run evaluates nodes in a different order, so it may expand critical nodes the sequential run never touched (detrimental), or find the incumbent that prunes a large region far earlier (acceleration). Both are consequences of the incumbent coupling the subtrees. A search without pruning, such as the parity instance of the branch-and-bound figure in Section 3.1, has neither. The script builds the two trees and runs the synchronous model on them. It then builds a third tree that marks the limit of the next theorem, tests that theorem on six thousand random trees, and ends with the depth-first example of Theorem 6.2.6.
# The anomalies of Lai and Sahni on explicit trees.
#
# Minimization; a leaf's bound is its value. Synchronous best-first
# model: each iteration expands the n open nodes of smallest bound,
# oldest first among ties.
import heapq
import itertools
import random
class Tree:
def __init__(self):
self.bound, self.kids, self.leaf = [], [], []
def add(self, bound, leaf=False, parent=None):
i = len(self.bound)
self.bound.append(bound)
self.kids.append([])
self.leaf.append(leaf)
if parent is not None:
self.kids[parent].append(i)
return i
def chain(self, parent, bounds):
"""A path of internal nodes below parent; returns its last node."""
for b in bounds:
parent = self.add(b, parent=parent)
return parent
def bush(self, parent, bound, arity, depth, leaf):
"""A complete subtree of equal bounds, every leaf worth `leaf`."""
level = [parent]
for _ in range(depth):
level = [self.add(bound, parent=p)
for p in level for _ in range(arity)]
for p in level:
self.add(leaf, leaf=True, parent=p)
def iterations(T, n):
"""I(n): iterations of the synchronous model with n processors."""
stamp = itertools.count()
live = [(T.bound[0], next(stamp), 0)]
inc = float("inf")
it = 0
while live:
it += 1
batch = [heapq.heappop(live) for _ in range(min(n, len(live)))]
if batch[0][0] >= inc:
# the smallest open bound cannot beat the incumbent: done
break
for b, _, v in batch:
if b >= inc:
continue
if T.leaf[v]:
inc = min(inc, b)
continue
for c in T.kids[v]:
if T.bound[c] < inc:
heapq.heappush(live, (T.bound[c], next(stamp), c))
return it
# 1. Detrimental anomaly: a chain 1, 2, 3 to the optimum 5, beside a bushy
# subtree of bound-5 nodes with leaves worth 6.
T = Tree()
r = T.add(0)
T.add(5, leaf=True, parent=T.chain(r, (1, 2, 3)))
T.bush(T.add(5, parent=r), 5, 3, 4, 6)
print("detrimental anomaly, I(n) for n = 1, 2, 3, 4, 8:",
[iterations(T, n) for n in (1, 2, 3, 4, 8)])
# 2. Acceleration anomaly, same model: older child A(5) -> A'(5) -> leaf 5;
# younger child heads 1, 2, 3, 4 -> Q(5) -> ternary bush of 5s,
# leaves 6.
T = Tree()
r = T.add(0)
T.add(5, leaf=True, parent=T.chain(r, (5, 5)))
T.bush(T.chain(r, (1, 2, 3, 4, 5)), 5, 3, 2, 6)
print("acceleration anomaly, I(n) for n = 1, 2, 4:",
[iterations(T, n) for n in (1, 2, 4)],
"on a tree of", len(T.bound), "nodes")
# 3. Distinct monotone bounds, and still I(3) > I(2): an anomaly between
# two parallel runs, none against I(1).
T = Tree()
r = T.add(0)
a = T.add(0.2454, parent=r)
T.add(12.2162, leaf=True, parent=a)
b = T.add(1.9521, parent=a)
for bd, leaves in ((3.5135, (4.0236,)),
(4.0900, (18.31, 12.83)),
(4.5784, (7.46, 6.95, 11.12))):
c = T.add(bd, parent=b)
for l in leaves:
T.add(l, leaf=True, parent=c)
print("distinct bounds, I(n) for n = 1, 2, 3, 4:",
[iterations(T, n) for n in (1, 2, 3, 4)],
"on a tree of", len(T.bound), "nodes")
# 4. Random trees with distinct monotone bounds (parent < child, no ties):
# the hypothesis of Theorem 6.2.5.
random.seed(1)
bad_seq = bad_low = bad_I1 = pairs = total = trees = 0
for trial in range(6000):
T = Tree()
fr = [T.add(0.0)]
depth = random.randint(3, 9)
maxb = random.randint(2, 4)
for d in range(depth):
nxt = []
for p in fr:
for _ in range(random.randint(1, maxb)):
leaf = d == depth - 1 or random.random() < 0.2
step = (15 + 5 * random.random() if leaf
else random.random())
nxt.append(T.add(T.bound[p] + step + 1e-12 * len(T.bound),
leaf, p))
fr = [c for c in nxt if not T.leaf[c]]
if not fr or len(T.bound) > 4000:
break
for c in fr:
T.add(T.bound[c] + 15 + 5 * random.random(), True, c)
zstar = min(T.bound[i] for i in range(len(T.bound)) if T.leaf[i])
# essential nodes
E = sum(1 for i in range(len(T.bound))
if not T.leaf[i] and T.bound[i] < zstar)
I = {m: iterations(T, m) for m in (1, 2, 3, 4, 5, 8)}
trees += 1
bad_I1 += I[1] not in (E + 1, E + 2)
for m in (2, 3, 4, 5, 8):
total += 1
bad_seq += I[m] > I[1]
bad_low += I[m] < -(-E // m) + 1
for m1, m2 in itertools.combinations((2, 3, 4, 5, 8), 2):
pairs += I[m2] > I[m1]
print(f"{trees} random trees:")
print(f" I(1) outside {{E+1, E+2}}: {bad_I1}")
print(f" runs with I(n) > I(1): {bad_seq} of {total}")
print(f" runs with I(n) < ceil(E/n)+1: {bad_low}")
print(f" pairs 1 < n1 < n2 with I(n2) > I(n1): {pairs} of {10 * trees}")
# 5. Depth-first with private stacks (Theorem 6.2.6); an idle processor
# takes the oldest node of the fullest stack.
# Left: a head of bound 1 over a binary subtree, bounds 10 to 14,
# with 32 leaves worth 50.
# Right: a chain 2, 3, 4, 5, 6, then the optimum 9.
def dfs_private(T, n):
stacks = [[] for _ in range(n)]
stacks[0].append(0)
inc = float("inf")
it = 0
while any(stacks):
it += 1
for p in range(n):
if not stacks[p]:
donor = max(range(n), key=lambda q: len(stacks[q]))
if len(stacks[donor]) > 1:
stacks[p].append(stacks[donor].pop(0))
else:
continue
v = stacks[p].pop()
if T.bound[v] >= inc:
continue
if T.leaf[v]:
inc = min(inc, T.bound[v])
continue
for c in reversed(T.kids[v]):
if T.bound[c] < inc:
stacks[p].append(c)
return it
T = Tree()
r = T.add(0)
level = [T.add(1, parent=r)]
for d in range(5):
level = [T.add(10 + d, parent=p) for p in level for _ in range(2)]
for p in level:
T.add(50, leaf=True, parent=p)
T.add(9, leaf=True, parent=T.chain(r, (2, 3, 4, 5, 6)))
print("depth-first with private stacks, I(n) for n = 1, 2, 4:",
[dfs_private(T, n) for n in (1, 2, 4)])
print(" on a tree of", len(T.bound), "nodes")
detrimental anomaly, I(n) for n = 1, 2, 3, 4, 8: [7, 8, 7, 7, 6]
acceleration anomaly, I(n) for n = 1, 2, 4: [13, 5, 5] on a tree of 30 nodes
distinct bounds, I(n) for n = 1, 2, 3, 4: [6, 5, 6, 5] on a tree of 13 nodes
6000 random trees:
I(1) outside {E+1, E+2}: 0
runs with I(n) > I(1): 0 of 30000
runs with I(n) < ceil(E/n)+1: 0
pairs 1 < n1 < n2 with I(n2) > I(n1): 0 of 60000
depth-first with private stacks, I(n) for n = 1, 2, 4: [71, 10, 9]
on a tree of 102 nodes
(What the script reports) Two processors take eight iterations where one takes seven, and eight processors take six, a speedup of \(1.17\) with eight processors. On the second tree two processors are \(2.6\) times faster than one, in the best-first model itself. The script costs one heap operation per node per processor count and runs in two seconds. Nothing in it is parallel, since it simulates the schedule rather than running it. The fourth part is the empirical side of the next theorem, for the comparison the theorem makes. On 6,000 random trees with distinct monotone bounds, no run with \(n\) processors took more iterations than the sequential run. Every run respected the lower bound \(\lceil E/n \rceil + 1\), and every sequential run took \(E + 1\) or \(E + 2\) iterations, with \(E\) the number of essential nodes. The third tree shows what the theorem does not say. Its thirteen bounds are distinct and monotone, and still \(I(3) = 6\) while \(I(2) = 5\). The last line belongs to Theorem 6.2.6, where it is explained.
Theorem 6.2.5 (Lai and Sahni, 1984; Li and Wah, 1986: no slowdown against the sequential run). Suppose the bounds are monotone and no two nodes have the same bound, so that the best-first order is a total order. Let \(E\) be the number of essential nodes. Then for every \(n \ge 1\),
\[\Big\lceil \frac{E}{n} \Big\rceil + 1 \;\le\; I(n) \;\le\; I(1) \;=\; E + 1 \ \text{ or } \ E + 2,\]so no detrimental anomaly against the sequential run can occur, and the speedup \(I(1)/I(n)\) is at most \(n\) up to the rounding of one or two iterations.G.-J. Li and B. W. Wah, "Coping with anomalies in parallel branch-and-bound algorithms", IEEE Transactions on Computers C-35 (1986), whose abstract states "sufficient conditions to guarantee no degradation in performance due to parallelism and necessary conditions for allowing parallelism to have a speedup greater than the number of processors"; T.-H. Lai and A. Sprague, "Performance of parallel branch-and-bound algorithms", IEEE Transactions on Computers C-34 (1985), and "A note on anomalies in parallel branch-and-bound algorithms with one-to-one bounding functions", Information Processing Letters 23 (1986). The displayed inequalities are the ones the proof sketch supports and the ones the script checks. The primary texts of Li and Wah and of Lai and Sprague were not consulted for this post; what is attributed to them here is taken from their abstracts and from the survey of Gendron and Crainic (1994), Section 3.
Proof sketch. With distinct monotone bounds the essential nodes form a subtree containing the root. Every run generates and expands every one of them, and no non-essential node is ever the minimum of the open list while an essential node is open. The sequential run therefore expands exactly the \(E\) essential nodes, then selects the optimal leaf (one iteration), then stops when the next minimum has bound \(\ge z_{\mathrm{inc}}\) (one further iteration unless the list is empty). With \(n\) processors each iteration expands at most \(n\) essential nodes, which gives the lower bound. Each iteration expands at least one essential node while any is open, which gives \(I(n) \le E + 1\) plus the same final check. ∎
(What the theorem compares, and what it leaves open) With a strict order on bounds the parallel best-first run does the sequential run's essential work in fewer rounds. The rest of its work is speculative work on non-essential nodes, which costs nothing because it fills processors that would otherwise idle. The theorem compares every parallel run with the sequential run and with nothing else. Between two parallel runs a detrimental anomaly can occur even under its hypotheses. The third tree of the script is the example. With three processors the third one expands, in iteration 4, a non-essential node (bound \(4.5784\)) that two processors pop only in iteration 5, after the optimal leaf has set the incumbent, and prune. Both runs select the optimal leaf in iteration 5; the sixth iteration of the three-processor run pops the one child of that node still open, a leaf worth \(11.12\), and stops. Lai and Sprague's 1986 note treats exactly this setting of one-to-one bounding functions. It gives conditions under which the anomaly cannot occur when the processor count is doubled. The note itself was not read for this post. The two constructions of Theorem 6.2.4 show what the random trees of the script exclude: ties among critical nodes, and, in the depth-first example below, a selection rule that is not best-first. The practical consequence is cheap insurance. Break ties by a fixed total order (node identifier, then depth), and the detrimental anomaly against the sequential run is ruled out for best-first search. Section 6.4 returns to this point. A GPU reduction can change a bound by one unit in the last place (one ulp, the spacing of floating-point numbers at the value), and that is exactly a way of breaking ties differently on each run.
Theorem 6.2.6 (Rao and Kumar, 1988, 1993: superlinear speedup on average for depth-first search). Consider parallel depth-first search for one solution with \(P\) processors, each exploring a disjoint part of the tree. If the solutions are uniformly distributed over the leaves the expected speedup is linear in \(P\), and if their distribution is non-uniform it is superlinear. For heuristic depth-first search, with the children ordered by a heuristic, the expected speedup is at least linear and superlinear on a subset of instances.V. N. Rao and V. Kumar, "Superlinear speedup in parallel state-space search", Foundations of Software Technology and Theoretical Computer Science, LNCS 338 (1988); V. N. Rao and V. Kumar, "On the efficiency of parallel backtracking", IEEE Transactions on Parallel and Distributed Systems 4 (1993), whose abstract states the three cases. The asynchronous branch-and-bound version is A. de Bruin, G. A. P. Kindervater and H. W. J. M. Trienekens, "Asynchronous parallel branch and bound and anomalies", IRREGULAR 1995, LNCS 980.
Proof sketch. Sequential depth-first search visits the leaves in a fixed order and stops at the first solution. Its expected work is the expected position of the first solution in that order. \(P\) processors each scan a \(1/P\) share, and the search stops when any share reaches a solution, so the time is the minimum over \(P\) such positions. For a uniform distribution the minimum of \(P\) scaled copies has mean \(1/P\) of the sequential mean. When solutions cluster, the sequential order may have the cluster late while one of the \(P\) shares has it early, and the expected minimum falls faster than \(1/P\). ∎
(Why a depth-first speedup can exceed the processor count) A depth-first order is a guess about where the good leaves are. \(P\) processors make \(P\) guesses at once, and the best of \(P\) unequal guesses beats the first guess by more than a factor of \(P\). The last tree of the script is the mechanism in miniature, with one private depth-first stack per processor and an idle processor taking the oldest node of the fullest stack. One processor expands all 63 internal nodes of the left subtree, which nothing prunes, before it reaches the chain on the right: 71 iterations. With two processors the second takes the right subtree in the first iteration, pops the five chain nodes in iterations 1 to 5, and finds the optimum \(9\) in iteration 6. From then on every left node is pruned when popped: ten iterations, a speedup of \(7.1\) with two processors. For branch and bound this is the mechanism of the acceleration anomaly: the early incumbent prunes. It is also the reason a superlinear speedup is usually a statement about the sequential algorithm rather than about the parallel one. A sequential code given the parallel run's incumbent would have been just as fast. A sequential code that diversifies its own search, by restarts or by a portfolio of dives, captures most of the same effect on one processor.
Theorem 6.2.6 in miniature: private depth-first stacks, 102 nodes
root
______/ \______
/ \
o chain 2
/ \ |
10 10 3
/ \ / \ |
11 11 11 11 4
.. .. .. .. |
five levels, bounds 10 to 5
14, with 32 leaves worth 50 |
below them: 63 internal 6
nodes with their head o |
optimum 9
n = 1: the stack dives left first and expands all 63 internal
nodes before it reaches the chain: 71 iterations
n = 2: worker 2 steals chain 2 in iteration 1, pops the chain
in iterations 1 to 5 and finds 9 in iteration 6; from then
on every left node is pruned when popped: 10 iterations,
a speedup of 7.1
n = 4: 9 iterations
Load balancing and work stealing
The synchronous model has one open list and an implicit barrier at every iteration. Real systems let every worker own a list and move work between lists only when a worker runs dry. The analysis of that design has two parts: a theorem about how well random stealing balances the load, and a theorem about how much work must be present for the balancing to be cheap.
Definition 6.2.7 (work stealing; work and span). Each worker owns a double-ended queue of nodes. The owner pushes and pops at the bottom, which is depth-first, cache-warm and gives the newest nodes. An idle worker, a thief, steals from the top of a randomly chosen victim's queue, which holds the oldest node and therefore the root of the largest remaining subtree. A multithreaded computation is a directed acyclic graph of unit tasks. Its work \(T_1\) is the number of tasks and its span \(T_\infty\) is the length of its longest path. For a branch-and-bound tree with a fixed set of expanded nodes, \(T_1\) is the node count times the node time and \(T_\infty\) is the depth times the node time.The lock-free deque with one compare-and-swap in the contended case is D. Chase and Y. Lev, "Dynamic circular work-stealing deque", SPAA 2005; the memory fences it needs under the C11/C++11 model are given by N. M. Lê, A. Pop, A. Cohen and F. Zappa Nardelli, "Correct and efficient work-stealing for weak memory models", PPoPP 2013.
The deque the systems use is Chase and Lev's circular array. It is lock-free: the only operation that two threads contend for is a compare-and-swap. That is the atomic instruction that replaces a word by a new value only if the word still holds the value the caller expected, and reports whether it did.
Algorithm 6.2.8 (parallel branch and bound with per-worker deques and random work stealing). The pseudocode below gives the loop every worker runs, the two shared words it touches, and the invariant that no node is processed twice.
Algorithm 6.2.8 Parallel branch and bound with per-worker deques
and random work stealing
Input root node N0 with box B0; a bound oracle bar z(.); a branching
rule; P workers; tolerance eps
Output incumbent z_inc and its point, with
z_inc - min over open nodes of bar z(N) <= eps
State deque D_i per worker (the owner pops and pushes at the bottom,
thieves steal at the top); shared atomic z_inc (minimize);
shared counter G of open nodes (Proposition 6.2.13)
1. push N0 onto D_0; G <- 1;
z_inc <- +inf, or the value of a heuristic solution
2. every worker i, in parallel, repeats:
3. if D_i is non-empty: N <- pop the bottom of D_i
(depth-first: the newest node)
4. else: flush delta_i into G; if G = 0 stop;
pick a random victim j != i; N <- steal the top of D_j;
if the steal fails go to 3
5. if bar z(N) >= z_inc - eps: finish N (prune) and go to 3
6. if N is a leaf (integer feasible, box small enough):
z_inc <- min(z_inc, f(x_N)) by compare-and-swap; go to 3
7. branch N into N', N''; compute bar z(N'), bar z(N'');
push the children whose bound is below z_inc - eps onto the
bottom of D_i, the better child last
8. finish N; go to 3
Invariant
at every moment some open node contains an optimal point, unless
z_inc is already optimal; every node is popped by exactly one
worker (deque semantics), so no node is processed twice.
Cost per node
one bound (the dominant term for a MINLP node: an LP or an NLP)
and two deque operations (a store plus a release fence to push; a
sequentially consistent fence plus, in the one-element case, one
compare-and-swap to pop).
Parallel
all P workers run steps 3 to 8 independently; the only shared
state is z_inc (read at every node, written rarely) and G (flushed
in batches). Steals cost O(1) each and O(P T_inf) in total.
The algorithm's node order is a dive on each worker and a breadth-first spread across workers, since thieves take the oldest, highest nodes. That is the order the theory wants. The deep nodes are cheap to produce and consume locally, and the stolen nodes head large subtrees, so a steal moves a lot of work for one message.
Work stealing in Algorithm 6.2.8: worker 0's first dive and D_0
N0 D_0, before any steal
/ \ +------+
N1 S1 ----------------> | S1 | top: a thief steals the
/ \ |------| oldest node, the root of
N2 S2 -------------------> | S2 | the largest remaining
/ \ |------| subtree
N3 S3 ----------------------> | S3 | bottom: the owner pushes
^ +------+ and pops, depth-first
N3 is being processed
step 7 pushes the children, the better child last, and step 3
pops it next, so one deferred sibling per level stays in D_0
shared by all workers: z_inc (compare-and-swap) and G, into
which each worker flushes its delta_i
Theorem 6.2.9 (Blumofe and Leiserson, 1999: randomized work stealing). For a fully strict multithreaded computation with work \(T_1\) and span \(T_\infty\), executed on \(P\) processors by the randomized work-stealing scheduler, the expected execution time is \(T_1/P + O(T_\infty)\). With probability at least \(1 - \delta\) the time is \(T_1/P + O(T_\infty + \lg P + \lg(1/\delta))\). The space is at most \(S_1 P\), where \(S_1\) is the sequential space, and the expected total communication is \(O(P\,T_\infty\,(1 + n_d)\,S_{\max})\), with \(S_{\max}\) the largest activation record and \(n_d\) the maximum number of times a thread synchronizes with its parent. All three bounds are existentially optimal to within a constant factor.R. D. Blumofe and C. E. Leiserson, "Scheduling multithreaded computations by work stealing", Journal of the ACM 46 (1999); "existentially optimal" is the abstract's phrase, meaning that for each bound some computation exists on which no scheduler does better. The reference implementation is Cilk-5: M. Frigo, C. E. Leiserson and K. H. Randall, "The implementation of the Cilk-5 multithreaded language", PLDI 1998.
A multithreaded computation is fully strict when every thread synchronizes only with its parent. A branch-and-bound tree, in which a node spawns its children and never waits for them, is fully strict with \(n_d = 0\). An activation record is the memory a thread holds while it runs, which here is one node.
Proof sketch. Charge each processor step either to work, of which there are at most \(T_1\) in total, or to a steal attempt. A potential-function argument shows that \(O(P)\) steal attempts reduce the potential of the critical path by a constant factor in expectation, so \(O(P\,T_\infty)\) steal attempts suffice with high probability. Dividing the \(T_1 + O(P\,T_\infty)\) processor steps among \(P\) processors gives the time bound. The space bound follows because each processor's deque holds a stack whose depth is at most the sequential stack depth. ∎
(Parallel slackness, and one example at scale) Work stealing is efficient whenever there is parallel slackness, \(T_1/(P\,T_\infty) \gg 1\), because then steals are rare and each moves more work than it costs. When \(P\) approaches \(T_1/T_\infty\) the \(O(T_\infty)\) term is the ramp-up and ramp-down of the next sub-subsection in disguise. A tree of depth \(d\) with node time \(t\) cannot keep more than about \(T_1/(d\,t)\) workers busy, whatever the scheduler. One example fixes the scale. Take a tree of \(T_1 = 10^9\) node-units with span \(T_\infty = 2 \times 10^4\) node-units, which is a depth of 200 at 100 units per node, or of 20,000 at one. At \(P = 100\) the slackness is \(500\), and the bound \(T_1/P + T_\infty\) gives an efficiency of at least \(0.998\). At \(P = 10^4\) the slackness is \(5\) and the efficiency is at least \(0.833\). At \(P = 10^5\) the slackness is \(0.5\) and the guarantee falls to \(0.333\), and at \(P = 10^6\) it is \(0.048\). The constant hidden in the \(O(\cdot)\) makes the real numbers worse. The caveat of Section 6.1 also applies: \(T_1\) and \(T_\infty\) are properties of the tree actually executed, which the schedule itself changes.
The second theorem asks how much work is needed to keep \(P\) processors efficient when the balancing is done by requests rather than by a shared deque. That is the distributed-memory setting of the supercomputer codes of Section 6.3.
Proposition 6.2.10 (isoefficiency of random polling; Kumar, Grama and Vempaty, 1994; Sanders, 2002). Let \(W\) be the work of a depth-first search whose subtrees are split on request, with the guarantee that a donor keeps at least a fraction \(\alpha\) of its work, and let a request and a transfer cost \(O(1)\) messages. If idle processors poll random victims, every processor has been asked at least once after \(O(P \log P)\) requests in expectation, the total communication is \(O(P \log P \log W)\), and the isoefficiency function is \(W = \Theta(P \log^2 P)\) up to contention.V. Kumar, A. Y. Grama and N. R. Vempaty, "Scalable load balancing techniques for parallel computers", Journal of Parallel and Distributed Computing 22 (1994); P. Sanders, "Randomized receiver initiated load-balancing algorithms for tree-shaped computations", The Computer Journal 45 (2002), whose abstract states upper bounds for random polling that match lower bounds for a large class of algorithms and machines, including a fully asynchronous model. The comparison with asynchronous and global round robin, which both come out at \(O(P^2 \log P)\), and the lower bound \(\Omega(P \log P)\) for every request-based scheme, are reproduced from memory of Grama, Gupta, Karypis and Kumar (2003), Section 11.4; they were not re-checked against that text or against the 1994 paper, and they are not part of the statement above for that reason.
Proof sketch. After every processor has been asked once, the largest piece of work anywhere has shrunk by the factor \(1 - \alpha\), so \(\log_{1/(1 - \alpha)} W\) such rounds reduce the largest piece to a constant. For random polling the number of requests until every processor has been asked is the coupon-collector time, \(P \ln P\) in expectation. Setting the communication \(O(P \log P \log W)\) equal to a constant fraction of \(W\) gives \(W = \Theta(P \log^2 P)\). ∎
(What the proposition says in numbers) The content of the proposition is the exponent of \(\log P\). Take a run on \(1{,}024\) processors and multiply the processor count by \(8\), \(64\) or \(1{,}024\). To hold the efficiency, random polling needs \(13.5\), \(164\) or \(4{,}096\) times the work, the ratios of \(P \log_2^2 P\). That is close to linear, which is why random stealing is the default on one machine and random polling the default across machines. It is also a statement about problem size: a tree that is barely large enough for a thousand processors is far too small for a million.
Ramp-up, ramp-down and termination
Definition 6.2.11 (phases of a parallel tree search). The ramp-up phase runs from the start until enough open nodes exist to keep every worker busy for the first time. The primary phase runs from the first to the last moment at which every worker is busy. The ramp-down phase runs from the last such moment to the end.The vocabulary is from the ParaSCIP reports: Y. Shinano, T. Achterberg, T. Berthold, S. Heinz, T. Koch and M. Winkler, "Solving hard MIPLIB2003 problems with ParaSCIP on supercomputers: an update", ZIB-Report 13-66 (2013), also IPDPS Workshops 2014.
Idle time in the two outer phases is not removable by any load balancer, because no scheduler can create nodes faster than the tree produces them. The simplest bound makes this quantitative.
Proposition 6.2.12 (ramp-up bound for a binary tree). A binary tree of \(N\) unit-time nodes explored breadth-first by \(P\) synchronous workers needs at least
\[T(P) \;\ge\; \lceil \log_2 P \rceil \;+\; \Big\lceil \frac{N - 2^{\lceil \log_2 P \rceil} + 1}{P} \Big\rceil\]rounds, since in round \(k\) at most \(2^{k-1}\) nodes exist. Hence \(E(P) \le N/(P\,T(P))\).
Proof. Let \(K = \lceil \log_2 P \rceil\). In round \(k\) at most \(2^{k-1}\) nodes exist, so rounds \(1\) to \(K\) process at most \(1 + 2 + \dots + 2^{K-1} = 2^K - 1\) nodes. The remaining \(N - 2^K + 1\) nodes are processed at most \(P\) per round, which takes at least their count divided by \(P\) further rounds. ∎
For \(N = 10^6\) the bound gives an efficiency of at most \(0.998\) at \(P = 256\), \(0.954\) at \(P = 4{,}096\) and \(0.492\) at \(P = 65{,}536\). At the last of these the search needs at least \(31\) rounds, and \(16\) of them are spent waiting for the tree to widen. Real trees are not binary and not balanced, so the real numbers are worse. This is the arithmetic behind the ZIB group's observation that at tens of thousands of cores the idle time of ramp-up and ramp-down dominates the running time unless the tree is enormous. It is also the reason racing ramp-up exists (Section 6.3): it gives the processors something to do during ramp-up other than wait.Koch, Ralphs and Shinano (2012), cited above, identify idle time in ramp-up and ramp-down as the dominant overhead at tens of thousands of cores.
Proposition 6.2.12: a binary tree widens by doubling
round 1 o at most 1 node
round 2 o o at most 2 nodes
round 3 o o o o at most 4 nodes
round 4 o o o o o o o o at most 8 nodes
: : :
round k at most 2^(k-1)
rounds 1 to K = ceil(log2 P): fewer than P nodes exist, so
workers idle; after them, at most P nodes per round
N = 10^6, P = 65,536: at least 31 rounds, 16 of them waiting;
efficiency at most 0.492 (0.954 at P = 4,096, 0.998 at P = 256)
The script below runs the shared-pool model of Definition 6.2.2 on a thirty-item knapsack whose tree has about \(1.3 \times 10^5\) nodes under Dantzig's bound of Section 3.1. It works in synchronous rounds of unit node time. It reports the phases and the efficiency for best-first and depth-first order, with the incumbent shared the moment it is found or kept private to each worker. The instance maximizes.G. B. Dantzig, "Discrete-variable extremum problems", Operations Research 5 (1957). The bound fills the knapsack greedily by value per unit of weight and adds the fitting fraction of the first item that does not fit, which is the value of the LP relaxation.
# P workers draw from one shared pool for a 0/1 knapsack branch and bound.
#
# Dantzig bound, one time unit per node, synchronous rounds. Best-first
# (the pool is a heap) or depth-first (a stack); the incumbent is shared
# the moment it is found, or private to each worker. Prints rounds T(P),
# nodes, speedup, efficiency, ramp-up, ramp-down.
import heapq
import random
random.seed(5)
n = 30
# strongly correlated: hard for the bound
w = [random.randint(20, 80) for _ in range(n)]
v = [wi + 10 for wi in w]
cap = sum(w) // 2
order = sorted(range(n), key=lambda i: -v[i] / w[i])
w = [w[i] for i in order]
v = [v[i] for i in order]
def bound(k, wt, val):
for i in range(k, n):
if wt + w[i] <= cap:
wt += w[i]
val += v[i]
else:
return val + (cap - wt) * v[i] / w[i]
return val
def run(P, best_first=True, share=True):
pool = [(-bound(0, 0, 0), 0, 0, 0)]
pop = heapq.heappop if best_first else (lambda l: l.pop())
push = heapq.heappush if best_first else (lambda l, x: l.append(x))
inc = [0] * P
shared = 0
t = 0
nodes = 0
busy_hist = []
while pool:
t += 1
busy = 0
batch = [(p, pop(pool)) for p in range(min(P, len(pool)))]
for p, (nb, k, wt, val) in batch:
cur = shared if share else inc[p]
if -nb <= cur:
# pruned on pop: costs no bound evaluation
continue
busy += 1
nodes += 1
if k == n:
inc[p] = max(inc[p], val)
if share:
shared = max(shared, val)
continue
# depth-first pops the 1-branch first (pushed last)
for take in (0, 1):
if take and wt + w[k] > cap:
continue
nw, nv = wt + w[k] * take, val + v[k] * take
b = bound(k + 1, nw, nv)
if b > (shared if share else inc[p]):
push(pool, (-b, k + 1, nw, nv))
busy_hist.append(busy)
full = [i for i, b in enumerate(busy_hist) if b == P]
return (t, nodes, (full[0] if full else t),
(t - full[-1] - 1 if full else 0), max(inc))
for bf in (True, False):
T1, N1, _, _, best = run(1, bf)
if not bf:
print()
print(f"{'best-first' if bf else 'depth-first'}: "
f"knapsack n={n}, optimum {best};")
print(f" sequential rounds {T1}, nodes {N1}")
print(f"{'':>4} |{'':>46} | private")
print(f"{'P':>4} |{'T(P)':>7}{'nodes':>7}{'speedup':>8}{'eff':>6}"
f"{'ramp-up':>8}{'ramp-down':>10} |{'T(P)':>7}{'nodes':>7}"
f"{'speedup':>8}")
for P in (1, 4, 16, 64, 256):
T, N, ru, rd, _ = run(P, bf, True)
Tn, Nn, _, _, _ = run(P, bf, False)
print(f"{P:>4} |{T:>7}{N:>7}{T1 / T:>8.2f}{T1 / T / P:>6.3f}"
f"{ru:>8}{rd:>10} |{Tn:>7}{Nn:>7}{T1 / Tn:>8.2f}")
best-first: knapsack n=30, optimum 991;
sequential rounds 165694, nodes 132376
| | private
P | T(P) nodes speedup eff ramp-up ramp-down | T(P) nodes speedup
1 | 165694 132376 1.00 1.000 0 33318 | 165694 132376 1.00
4 | 41426 132377 4.00 1.000 2 8331 | 41426 132380 4.00
16 | 10365 132422 15.99 0.999 4 2086 | 10365 132437 15.99
64 | 2605 132711 63.61 0.994 6 527 | 2605 132774 63.61
256 | 669 134240 247.67 0.967 8 138 | 669 134495 247.67
depth-first: knapsack n=30, optimum 991;
sequential rounds 121615, nodes 121595
| | private
P | T(P) nodes speedup eff ramp-up ramp-down | T(P) nodes speedup
1 | 121615 121595 1.00 1.000 0 0 | 121615 121595 1.00
4 | 30405 121570 4.00 1.000 2 2 | 30464 121636 3.99
16 | 7631 121675 15.94 0.996 4 10 | 8071 123404 15.07
64 | 1935 122095 62.85 0.982 6 13 | 2069 125216 58.78
256 | 516 124068 235.69 0.921 8 5 | 569 131045 213.73
(What the pool script shows) Ramp-up lasts \(\lceil \log_2 P \rceil\) rounds, as Proposition 6.2.12 predicts for a binary tree: \(2, 4, 6, 8\) rounds for \(P = 4, 16, 64, 256\). This tree is wide, with thousands of open nodes for most of the run, so the efficiency is \(0.967\) at \(P = 256\) with best-first order and \(0.921\) with depth-first order. Under best-first order the incumbent hardly matters until the end, because best-first never selects a node that the final incumbent would prune until the essential nodes are exhausted. The \(33{,}318\) "ramp-down" rounds of the sequential best-first run are rounds whose popped node was pruned, which is how a sequential best-first run ends. Under depth-first order the incumbent matters at every moment. With private incumbents the 256 workers do \(7{,}000\) more nodes and lose nine points of efficiency (\(0.835\) against \(0.921\)), because each worker must rediscover the solution that prunes its own dives. The knapsack tree is as regular as trees get. On a tree with a thin frontier the losses are far larger, and the only way to know the efficiency on a given MINLP instance is to measure it. The script costs one bound evaluation per node and is sequential, since it simulates the rounds.
The next question is how a set of workers with private deques knows that the search is over. A worker with an empty deque cannot conclude anything from its own state, because a node it finished may have had children stolen by someone else.
Proposition 6.2.13 (termination detection by counting). Let every worker keep a private counter \(\delta_i\) of the nodes it has pushed minus the nodes it has finished, and let a global counter \(G\) receive \(\delta_i\) whenever worker \(i\) flushes. If a worker flushes before every read of \(G\) and never exits while its own deque is non-empty, then \(G = 0\) observed by an idle worker implies that every node it could still reach has been processed. The protocol never loses a node, and all workers eventually exit.
Proof sketch. The true number of open nodes is \(G + \sum_i \delta_i\), counting the unflushed increments. A worker with \(\delta_j > 0\) unflushed has pushed nodes it has not finished. They sit in its own deque or have been stolen and are then counted by the thief. Worker \(j\) cannot exit with a non-empty deque, so every node is eventually popped and processed. A worker that reads \(G = 0\) after flushing may exit while other workers still hold work. Correctness is unaffected and only efficiency is lost, which is why the implementation in Section 6.3 also flushes after every steal. ∎
(Termination on distributed memory) On distributed memory the classical algorithms circulate a token around a ring and colour it when work has moved, and the counting scheme above is the shared-memory special case. The UG framework of Section 6.3 detects termination centrally, when the coordinator's pool is empty and every solver reports idle.E. W. Dijkstra, W. H. J. Feijen and A. J. M. van Gasteren, "Derivation of a termination detection algorithm for distributed computations", Information Processing Letters 16 (1983); F. Mattern, "Algorithms for distributed termination detection", Distributed Computing 2 (1987); for UG, Y. Shinano, T. Achterberg, T. Berthold, S. Heinz and T. Koch, "ParaSCIP: a parallel extension of SCIP", in Competence in High Performance Computing 2010 (Springer, 2012), also ZIB-Report 10-27.
(Sharing the incumbent across devices) The incumbent is the one piece of information every worker needs from every other. The pool script shows the two regimes: under best-first order it can arrive late without harm, under depth-first order it must arrive at once. Two design rules follow for a machine with several devices. Broadcast the incumbent at round boundaries rather than per node, since a device cannot act on it in the middle of a kernel anyway. Use an all-reduce of the minimum for the broadcast. An all-reduce is the collective operation that combines one value from every participant and leaves the result on all of them, and the collective libraries provide it as a primitive that overlaps with the bounding kernels.NVIDIA NCCL User Guide, "Collective Operations" (ncclAllReduce with ncclMin, ncclBroadcast, ncclSend and ncclRecv), docs.nvidia.com/deeplearning/nccl, read 5 October 2026. And never let a worker prune on a value it received without the solution that produced it, so that a run restarted from a checkpoint can reproduce its proof. In UG the incumbent and its solution travel together through the coordinator.
Determinism and its cost
A production solver is expected to return the same answer by the same path when it is run twice on the same input. On one processor that is automatic. On many it is a design decision with a price. On a GPU it is, as of this writing, not available for any MINLP solver. This sub-subsection defines the terms, gives the CPU design and its measured cost, and then identifies the two sources of nondeterminism that are specific to a device: the order of a floating-point reduction and the composition of a batch. It ends with the fact that makes the problem tractable: correctness and reproducibility are separate requirements with separate remedies.
Definition 6.2.14 (deterministic, opportunistic, thread-independent, reproducible; logical clock). The path of a run is the sequence of solver states it passes through: the nodes selected, the bounds computed, the incumbents found, the cuts added, in order. A parallel solver is deterministic if two runs with identical input, parameters, thread count and machine produce identical paths and output. It is thread-independent if the path does not depend on the thread count either, and opportunistic if the path may depend on timing. A result is reproducible if the output, though not necessarily the path, is identical. A logical clock is a counter of work that is a function of the computation alone: LP iterations, nonzeros touched, propagation steps, or kernel launches, a kernel being a function launched on the device (Section 7.1). A deterministic synchronization point is a value of the logical clock at which the workers exchange information in a fixed order.IBM's documentation states the first notion as "multiple runs with the same model at the same parameter settings on the same platform will reproduce the same solution path and results", and the opportunistic alternative as "even slight differences in timing among threads or in the order in which tasks are executed in different threads may produce a different solution path": IBM ILOG CPLEX Optimization Studio 22.1.2, "Parallel mode switch", ibm.com/docs/en/icos/22.1.2, read 5 October 2026. Gurobi counts work units and its WorkLimit "will stop the optimization of a given model at the exact same point every time" on the same hardware; CPLEX counts ticks and offers dettimelimit; UG's deterministic clock counts communication-point calls. Any decision keyed to wall time, such as a time limit or a time-based heuristic budget, breaks determinism, which is why Gurobi lists TimeLimit and NoRelHeurTime among the exceptions to its guarantee: Gurobi Optimizer Reference Manual 13.0, "Parameters", read 5 October 2026.
(Deterministic, reproducible, correct) The three words name three different things. A deterministic run reproduces its bugs every time, and a nondeterministic run with valid bounds is correct every time. Two runs that return the same optimal value by different paths are reproducible but not deterministic. Xpress, FiberSCIP's deterministic mode and the recent Para-B&B framework share one CPU design. Workers compute without communicating until their logical clocks reach a common value, then exchange in a fixed order.
Algorithm 6.2.15 (deterministic parallel branch and bound with a logical clock). The pseudocode below is that design, with the invariant that makes the state after every round a function of the input alone.
Algorithm 6.2.15 Deterministic parallel branch and bound
with a logical clock
Input problem; P workers; a quantum Q of work units; a canonical node
order (bound, then identifier)
Output optimal value and solution, identical for every run with the
same input, P and Q
State open list L (shared); incumbent z_inc; per-worker clock c_w;
per-worker outbox O_w
1. root: presolve and bound the root once; L <- {root};
z_inc <- +inf; c_w <- 0 for all w
2. for round r = 1, 2, ...:
3. assign (sequential, a fixed rule of (L, r)): remove nodes from L
in canonical order and give each worker its share
4. work (parallel, no communication): worker w processes its nodes;
children, improved incumbents and globally valid cuts go to O_w;
c_w advances by the work units spent;
w stops when c_w >= r Q or its nodes are exhausted
5. barrier: wait until every worker has stopped
6. merge (sequential, w = 1, ..., P): apply O_w in that order:
z_inc <- min(z_inc, found values);
insert the children into L with identifiers assigned in order;
add cuts to the pool;
prune L against z_inc
7. until L is empty; return z_inc and its solution
Invariant
the state (L, z_inc, pool) after step 6 of round r is a function of
the input, P and Q alone, because every decision in steps 3 and 6
is keyed to the logical clock and to canonical orders, and step 4
reads only state fixed in step 3.
Cost per round
the idle time at the barrier, which is the spread of the workers'
finishing times, plus the sequential merge; a stopping rule inside
the loop must be a work limit, since a wall-clock limit would make
the stopping round depend on time.
Parallel
step 4 in full; steps 3 and 6 are sequential and bound the
efficiency as any loop with a sequential merge is bounded
(Amdahl's law, Proposition 6.4.4).
One round r of Algorithm 6.2.15: assign, work, barrier, merge
3. assign (sequential): nodes leave L in canonical order
|
v
4. work (parallel, no communication) 5. barrier
w = 1 ===========================.............|
w = 2 ==================......................|
: |
w = P ========================================|
| ----------------------------------------+---> time
| a worker stops when c_w >= r Q or when
| its nodes are exhausted
v
6. merge (sequential): O_1, O_2, ..., O_P, in that order, into
L, z_inc and the pool; then round r + 1, until L is empty
= working: c_w advances by the work units spent
. idle time at the barrier: the spread of the finishing times
(The price, measured twice) The price has been measured twice. The Xpress Global paper runs its deterministic default against an opportunistic variant on 1,308 MINLP instances. The opportunistic variant's time ratio is \(0.98, 0.98, 0.97, 0.95\) and \(0.93\) at \(2, 4, 8, 16\) and \(32\) threads over all instances, and \(0.99\) to \(0.91\) on the node-heavy subset. That is a saving of "2–9 %", which the authors decline to take because the opportunistic numbers "are not really reproducible". The same paper notes that the solution path is identical across Windows, Linux and macOS.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), Table 8 and footnote 18; the test set and machine are given in its Section 3 (FICO Xpress 9.5.3, two sixteen-core Intel Xeon Gold 5218 sockets, one hour, sixteen threads in the experiments; Section 4.9 states that the deterministic mode is the default). The design of the deterministic tree search is T. Berthold, J. Farmer, S. Heinz and M. Perregaard, "Parallelization of the FICO Xpress-Optimizer", Optimization Methods and Software 33 (2018): "MIP solvers like Xpress are expected to be deterministic. This inevitably results in synchronization latencies." Para-B&B is a deterministic shared-memory parallelization of HiGHS. It replicates the solver state across workers and removes every nondeterministic synchronization primitive. It reports a geometric-mean speedup of \(2.17\) at eight threads over 80 MIPLIB 2017 instances, \(5.12\) on the node-heavy ones, with "thread idle rates averaging 34.7 %".J. Zhang, D. Huang, Y. Liu, S. Wang, Z. Pu and Z. Liu, "Para-B&B: load-balanced deterministic parallelization of solving MIP", arXiv 2604.09556 (2026), abstract. The distributed-memory counterpart with a sliding-window deterministic mode is L. Wang, J. Liu, F. Zhang, J. Wei, Y. Tang, J. Sun and X. Luo, "N2N: a parallel framework for large-scale MILP under distributed memory", arXiv 2511.18723 (2025), which reports speedups of \(22.52\) and \(12.71\) with 1,000 MPI processes on two clusters in its nondeterministic mode and gives no number for the deterministic one in the abstract. A third of the threads idle at a barrier is the price of a barrier per round on an irregular tree. It is the number to keep in mind when a GPU design proposes a barrier per batch.
On a device two further sources of nondeterminism appear, and both are absent from a CPU solver that uses a logical clock. The first is arithmetic.
Proposition 6.2.16 (floating-point addition is not associative). In IEEE binary64 under round-to-nearest, with \(a = 1\) and \(b = c = 10^{-16}\), \(\mathrm{fl}(\mathrm{fl}(a + b) + c) = 1\) while \(\mathrm{fl}(a + \mathrm{fl}(b + c)) = 1 + 2^{-52}\).
Proof. The spacing of binary64 numbers at \(1\) is \(2^{-52} \approx 2.22 \times 10^{-16}\), so \(\mathrm{fl}(1 + 10^{-16}) = 1\) twice over. But \(\mathrm{fl}(b + c) \approx 2 \times 10^{-16}\) exceeds half the spacing, so \(1 + \mathrm{fl}(b + c)\) rounds up to \(1 + 2^{-52}\). ∎
Proposition 6.2.16: the same three numbers, two addition trees
a = 1, b = c = 10^-16 in binary64; the spacing at 1 is 2^-52,
about 2.22 x 10^-16
fl(fl(a + b) + c) = 1 fl(a + fl(b + c)) = 1 + 2^-52
+ [1] + [1 + 2^-52]
/ \ / \
[1] + c a + [about 2 x 10^-16]
/ \ / \
a b b c
[v]: the rounded value of the sum at that node
1 + 10^-16 rounds to 1, twice; about 2 x 10^-16 exceeds half the
spacing, so 1 plus it rounds up to 1 + 2^-52
(Schedules, fixed and dynamic) A sum of \(n\) terms is a binary tree of additions, and its schedule is the pair (tree shape, assignment of terms to leaves). Recursive summation is the left comb, pairwise summation the balanced tree. The schedule is fixed if it is a function of \(n\) and of the index order alone. It is dynamic if it depends on the order in which operands become available. That happens when many threads accumulate into one word with atomic additions, or when work is stolen. It also happens when a library picks its kernel by the number of streaming multiprocessors (the device's processing units, Section 7.1) or by the active streams (its independent command queues). An atomicAdd on a double is a left comb whose leaf order is the order in which the memory system serializes the threads. The hardware does not specify that order, so a dual objective \(b^\top y\) accumulated that way can differ between two runs of the same kernel.NVIDIA, CUDA Programming Guide 13.4.2, appendix "C++ language extensions", section "Atomic functions": the atomic functions perform "atomic read-modify-write operations on a 32-, 64-, or 128-bit word", have "a memory ordering of cuda::std::memory_order_relaxed", and atomicAdd() supports float and double; read 5 October 2026. NVIDIA's own floating-point guide gives the single-precision form of the example above and concludes that "the order in which operations are executed affects the accuracy of the result. The results are independent of the host system": N. Whitehead and A. Fit-Florea, Floating Point and IEEE 754, CUDA 13.4 documentation, read 5 October 2026.
Theorem 6.2.17 (Higham, 1993: the error of a summation schedule). Let \(\hat s\) be the sum of \(x_1, \dots, x_n\) computed by a schedule in which leaf \(i\) is \(h_i\) additions below the root, under round-to-nearest with unit roundoff \(u\) (\(u = 2^{-53}\) in binary64, \(2^{-24}\) in binary32). Then
\[|\hat s - s| \;\le\; u \sum_{i=1}^n h_i\,|x_i| + O(u^2), \qquad s = \sum_{i=1}^n x_i .\]For recursive summation (\(h_i \le n - 1\)) the bound is \((n - 1)\,u \sum_i |x_i|\), and for pairwise summation (\(h_i \le \lceil \log_2 n \rceil\)) it is \(\lceil \log_2 n \rceil\, u \sum_i |x_i|\), both to first order. Two schedules with heights \(h\) and \(h'\) therefore differ by at most \((h + h')\,u \sum_i |x_i|\), and relative to the exact sum by at most \((h + h')\,u\,\kappa\) with \(\kappa = \sum_i |x_i| / |\sum_i x_i|\) the condition number of the summation.N. J. Higham, "The accuracy of floating point summation", SIAM Journal on Scientific Computing 14 (1993), which analyses five methods and shows that "no one method is uniformly more accurate than the others"; the two first-order bounds are the standard ones.
Proof sketch. Each internal node computes \(\mathrm{fl}(p + q) = (p + q)(1 + \theta)\) with \(|\theta| \le u\). Unrolling from the root, \(\hat s = \sum_i x_i \prod_{k=1}^{h_i}(1 + \theta_{i,k})\), and \(|\prod_k (1 + \theta_{i,k}) - 1| \le h_i u + O(u^2)\). ∎
Theorem 6.2.17: two schedules for the sum of x_1, ..., x_n
recursive summation: pairwise summation:
the left comb the balanced tree
+ +
/ \ / \
+ x_n + +
/ \ / \ / \
... x_(n-1) ... ... ... ...
/ /\ /\
+ x_1 x_2 x_(n-1) x_n
/ \
x_1 x_2
leaf i sits h_i additions below the root
h_i <= n - 1 h_i <= ceil(log2 n)
(Why cancellation makes the order matter) The condition number \(\kappa\) is what makes the matter serious for a solver. A dual objective \(b^\top y\) with terms of both signs, and the reduced costs \(c - A^\top y\) near a dual optimum, are cancelling sums with large \(\kappa\). A pruning test that compares such a bound with the incumbent to the last few digits is therefore decided by the order of accumulation. Two remedies exist, and the proposition and the listing below give both.
Proposition 6.2.18 (a fixed schedule is bitwise reproducible). If the schedule of a reduction is a function of \(n\) and the index order alone, and every addition is an IEEE operation in one fixed precision and rounding mode, then the computed sum is a function of the operand values alone: identical in every run and on every conforming processor, whichever threads compute whichever internal nodes.
Proof. Induction on the tree. A leaf's value is its operand. An internal node's value is \(\mathrm{fl}(\mathrm{left} + \mathrm{right})\), and IEEE 754 specifies \(\mathrm{fl}\) uniquely for given operands, precision and rounding mode. The order of evaluation does not enter. ∎
(Three caveats on a real kernel) Three caveats decide whether the proposition applies to a real kernel. A compiler that contracts a product and a sum into a fused multiply-add changes the operation sequence. A change of precision changes \(\mathrm{fl}\). Library functions such as exp are not required to be correctly rounded, so two math libraries differ. Within one binary on one device, however, a block reduction satisfies the hypothesis: a block reduction is a sum over a fixed tree held in shared memory, the fast memory private to one group of threads. The vendor libraries state the library-level form of the proposition. cuBLAS routines "generate the same bit-wise results at every run when executed on GPUs with the same architecture and the same number of SMs". The guarantee is lost across toolkit versions and when several streams share one workspace, which is exactly the fixed-versus-dynamic distinction.NVIDIA, cuBLAS Library 13.4, Section 2.1.4 "Results reproducibility", read 5 October 2026. The design discipline for bit-reproducibility across compilers and machines is A. Arteaga, O. Fuhrer and T. Hoefler, "Designing bit-reproducible portable high-performance applications", IPDPS 2014.
Proposition 6.2.19 (one-bin integer accumulation). Let \(m = \max_i |x_i| > 0\), choose the integer \(s = 62 - \lceil \log_2(n m) \rceil\) and let \(q_i = \mathrm{rint}(x_i 2^s) \in \mathbb Z\). Then every partial sum of the \(q_i\) has modulus below \(2^{63}\), so their sum in 64-bit two's-complement arithmetic is exact and identical in every order. Moreover \(|2^{-s} \sum_i q_i - \sum_i x_i| \le n\,2^{-(s+1)}\), and the conversion of the integer total to binary64 rounds once, identically for every order.
Proof. \(|q_i| \le m 2^s + 1/2\), so a partial sum has modulus at most \(n m 2^s + n/2 < 2^{62} + n/2 < 2^{63}\), and integer addition without overflow is exact, associative and commutative. Each \(|q_i 2^{-s} - x_i| \le 2^{-(s+1)}\). The final conversion is one rounding of one integer. ∎
(The one-bin method's error, and its many-bin refinement) The grid \(2^{-s}\) is set by the largest summand, so the absolute error is about \(n^2 m\,2^{-63}\), a thousand times below the worst case of recursive summation. The relative error on a cancelling sum is poor, because it is measured against \(m\) and not against \(|s|\). The reproducible summation of Demmel and Nguyen removes that weakness by keeping several bins at consecutive exponent ranges and pre-rounding each summand into them. It is a "reproducible accumulator" that sums floating-point numbers "independent of summation order" in one read-only pass and one parallel reduction, at a cost of about \(9n\) floating-point and \(3n\) bitwise operations for a six-word accumulator.J. Demmel and H. D. Nguyen, "Fast reproducible floating-point summation", ARITH 2013; J. Demmel and H. D. Nguyen, "Parallel reproducible summation", IEEE Transactions on Computers 64 (2015); W. Ahrens, J. Demmel and H. D. Nguyen, "Algorithms for efficient reproducible floating point summation", ACM Transactions on Mathematical Software 46 (2020), article 22, whose abstract gives the operation counts. The GPU construction with a long fixed-point accumulator and error-free transformations, of which Proposition 6.2.19 is the single-word case, is C. Collange, D. Defour, S. Graillat and R. Iakymchuk, "Numerical reproducibility for the parallel reduction on multi- and many-core architectures", Parallel Computing 49 (2015). At \(n = 1{,}024\) the recursive bound is about a hundred times the pairwise one and about a thousand times the one-bin one, which is the sense of "a thousand times below" above.
The listing computes the sum of a batch of a thousand values with magnitudes spread over several decades and both signs, three ways. It sums recursively in arrival order, by a pairwise tree whose pairing is fixed by the index, and by the one-bin integer accumulation. It then shuffles the arrival order a thousand times and counts the distinct bit patterns each method produces.
// Two order-independent reductions over a batch of numbers (node bounds,
// dual products), C++23.
//
// fixed_tree_sum: a pairwise tree whose pairing is fixed by the index, so
// the result is a function of the operands in their canonical order;
// any thread may compute any pair, and it is what a deterministic GPU
// block reduction does. It is not invariant under permutation of the
// operands, so the batch must be put in a canonical order (by node id)
// before the reduction.
// int_bin_sum: one-bin integer accumulation, the single-bin case of
// Demmel and Nguyen's reproducible summation: round every summand to a
// common grid 2^-s chosen from n and max|x|, add exactly in int64 (no
// overflow by construction, so the sum is the same in every order),
// unscale once.
// Check: 1,000 random permutations of 1,000 values; count distinct bit
// patterns of each reduction.
// Compile check: clang++ -std=c++23 -fsyntax-only -Wall -Wextra reduce.cpp
#include <algorithm>
#include <bit>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <numeric>
#include <random>
#include <set>
#include <span>
#include <vector>
double recursive_sum(std::span<const double> x) {
double s = 0.0;
for (double v : x) {
s += v;
}
return s;
}
// n - 1 additions, depth ceil(log2 n)
double fixed_tree_sum(std::span<const double> x) {
std::vector<double> a(x.begin(), x.end());
while (a.size() > 1) {
const std::size_t h = a.size() / 2, odd = a.size() % 2;
// reads stay ahead of writes
for (std::size_t i = 0; i < h; ++i) {
a[i] = a[2 * i] + a[2 * i + 1];
}
if (odd) {
a[h] = a[a.size() - 1];
}
a.resize(h + odd);
}
return a.empty() ? 0.0 : a[0];
}
// n scale-and-round, n - 1 exact integer adds
double int_bin_sum(std::span<const double> x) {
double m = 0.0;
for (double v : x) {
m = std::max(m, std::fabs(v));
}
const int s = 62 - static_cast<int>(
std::ceil(std::log2(static_cast<double>(x.size()) * m)));
// n * max|x| * 2^s < 2^62: no overflow in any order
std::int64_t acc = 0;
for (double v : x) {
acc += static_cast<std::int64_t>(std::llrint(std::ldexp(v, s)));
}
return std::ldexp(static_cast<double>(acc), -s);
}
int main() {
const std::size_t n = 1000;
std::mt19937_64 gen(20261005);
std::lognormal_distribution<double> mag(0.0, 2.5);
std::bernoulli_distribution sign(0.5);
std::vector<double> x(n);
for (double& v : x) {
v = (sign(gen) ? 1.0 : -1.0) * mag(gen);
}
std::vector<std::size_t> idx(n);
std::iota(idx.begin(), idx.end(), std::size_t{0});
std::set<std::uint64_t> rec, tree_perm, tree_canon, bin;
std::vector<double> y(n);
for (int trial = 0; trial < 1000; ++trial) {
// the arrival order of the batch
std::shuffle(idx.begin(), idx.end(), gen);
for (std::size_t i = 0; i < n; ++i) {
y[i] = x[idx[i]];
}
rec.insert(std::bit_cast<std::uint64_t>(recursive_sum(y)));
tree_perm.insert(std::bit_cast<std::uint64_t>(fixed_tree_sum(y)));
bin.insert(std::bit_cast<std::uint64_t>(int_bin_sum(y)));
// restore the canonical order (by id) before the tree
std::vector<double> z = y;
std::vector<std::size_t> pos(n);
for (std::size_t i = 0; i < n; ++i) {
pos[idx[i]] = i;
}
for (std::size_t i = 0; i < n; ++i) {
z[i] = y[pos[i]];
}
tree_canon.insert(std::bit_cast<std::uint64_t>(fixed_tree_sum(z)));
}
std::printf("n = %zu values, 1000 arrival orders: distinct results\n",
n);
std::printf(" recursive sum in arrival order %zu\n",
rec.size());
std::printf(" fixed tree on the arrival order %zu\n",
tree_perm.size());
std::printf(" fixed tree on the canonical order %zu"
" value %.17g\n",
tree_canon.size(), fixed_tree_sum(x));
std::printf(" one-bin integer accumulation %zu"
" value %.17g\n",
bin.size(), int_bin_sum(x));
return 0;
}
n = 1000 values, 1000 arrival orders: distinct results
recursive sum in arrival order 273
fixed tree on the arrival order 45
fixed tree on the canonical order 1 value 106.09006385447151
one-bin integer accumulation 1 value 106.09006385447196
(What the listing shows) The recursive sum takes 273 different values over the thousand orders. The fixed tree still takes 45 when its leaves are assigned in arrival order. A fixed tree is therefore not enough: the operand order must be fixed too, which in a batch of nodes means sorting the batch by node identifier before the reduction. Restoring the canonical order makes the tree reproducible, and the integer accumulation is reproducible in any order. The two reproducible values differ in the sixteenth significant digit, as Proposition 6.2.19 allows. They are two different correctly computed quantities, each identical across all thousand orders. The fixed tree does the same \(n - 1\) additions as any balanced reduction at depth \(\lceil \log_2 n \rceil\). The integer method adds one scale-and-round per element and one exact maximum first, itself an order-free reduction. Every level of the tree and every element of the integer pass parallelizes, and the final integer reduction is itself a fixed tree.
The listing's three reductions over 1,000 arrival orders
1,000 values, shuffled 1,000 times
|
+--> recursive sum, in arrival order .......... 273 results
|
+--> fixed tree, leaves in arrival order ...... 45 results
|
+--> back to the canonical order (by id)
| |
| +--> fixed tree ....................... 1 result
| 106.09006385447151
|
+--> scale each value by 2^s, round it to int64
|
+--> exact integer sum in any order, one
conversion ....................... 1 result
106.09006385447196
A difference in the last bit of a bound matters to a tree search because of the ties.
Proposition 6.2.20 (marginal decisions). Consider a branch and bound that prunes node \(N\) when its computed bound \(\hat z(N)\), the floating-point value obtained for the exact relaxation value \(\bar z(N)\), satisfies \(\hat z(N) \ge z_{\mathrm{inc}} - \varepsilon\). Let two runs compute bounds that each differ from \(\bar z(N)\) by at most \(\eta\,|\bar z(N)|\) at every node they visit. Then the first node at which the two runs decide differently satisfies \(|\,\bar z(N) - (z_{\mathrm{inc}} - \varepsilon)\,| \le \eta\,|\bar z(N)|\) at the time of its test. Every difference between the two trees originates in such a marginal decision or in an earlier difference.
Proof. The statement restates the pruning test. ∎
(How a last-bit difference reaches the tree) Arithmetic nondeterminism of relative size \(\eta\) acts on the tree only through the decisions that were within \(\eta\) of the threshold. Once two runs expand different nodes the anomaly theorems take over. Theorem 6.2.4 says a different expansion order can raise or lower the node count, and Theorem 6.2.5 says detrimental anomalies against the sequential run need ties, non-monotone bounds, or a selection rule that is not best-first. A marginal decision is a near-tie, so the condition under which the no-anomaly theorem protects a run is exactly the absence of marginal decisions. Scheduling nondeterminism enters through the same door without any arithmetic difference, because it decides ties. The script measures both on the knapsack of the pool script. It runs a best-first search several times under three perturbations and under two batch schedules. The batch modes anticipate Section 6.4.
# Node counts of a best-first branch and bound under three kinds of
# nondeterminism, on the knapsack of the shared-pool script (30 items).
#
# D: deterministic, ties broken by node identifier.
# N: one ulp of relative noise on every bound before the pruning test,
# as an order-dependent reduction would produce.
# T: ties broken at random (scheduling noise, exact arithmetic).
# B1/B2: frontier batches of b nodes with a fixed, or a randomized,
# batch schedule.
# Every run returns the same optimum; only the node counts move.
import heapq
import random
random.seed(5)
n = 30
w = [random.randint(20, 80) for _ in range(n)]
v = [wi + 10 for wi in w]
cap = sum(w) // 2
order = sorted(range(n), key=lambda i: -v[i] / w[i])
w = [w[i] for i in order]
v = [v[i] for i in order]
def bound(k, wt, val):
for i in range(k, n):
if wt + w[i] <= cap:
wt += w[i]
val += v[i]
else:
return val + (cap - wt) * v[i] / w[i]
return float(val)
def serial(noise=0.0, random_ties=False, rng=None):
inc, evals, marginal, ident = 0, 0, 0, 0
key = lambda: rng.random() if random_ties else ident
pool = [(-bound(0, 0, 0), key(), 0, 0, 0)]
while pool:
negb, _, k, wt, val = heapq.heappop(pool)
bnd = -negb
if abs(bnd - inc) <= 1e-9 * max(1, abs(inc)):
# a decision noise could flip
marginal += 1
if noise:
bnd *= 1.0 + noise * (2 * rng.random() - 1)
if bnd <= inc:
continue
evals += 1
if k == n:
inc = max(inc, val)
continue
for take in (1, 0):
if take and wt + w[k] > cap:
continue
nw, nv = wt + w[k] * take, val + v[k] * take
bb = bound(k + 1, nw, nv)
if bb > inc:
ident += 1
heapq.heappush(pool, (-bb, key(), k + 1, nw, nv))
return inc, evals, marginal
def batched(b, randomized=False, rng=None):
inc, evals, ident, pending = 0, 0, 0, []
pool = [(-bound(0, 0, 0), 0, 0, 0, 0)]
while pool or pending:
if randomized:
# a random b-subset of the 2b best nodes forms the batch
cand = [heapq.heappop(pool)
for _ in range(min(2 * b, len(pool)))]
rng.shuffle(cand)
batch, rest = cand[:b], cand[b:]
for node in rest:
heapq.heappush(pool, node)
else:
batch = [heapq.heappop(pool) for _ in range(min(b, len(pool)))]
# the device knows only inc0 during the round
inc0, children, found = inc, [], []
for negb, _, k, wt, val in batch:
if -negb <= inc0:
continue
evals += 1
if k == n:
found.append(val)
continue
for take in (1, 0):
if take and wt + w[k] > cap:
continue
nw, nv = wt + w[k] * take, val + v[k] * take
bb = bound(k + 1, nw, nv)
if bb > inc0:
children.append((bb, k + 1, nw, nv))
# a randomized schedule delivers some incumbents a round late
now, later = [], []
for x in found:
if randomized and rng.random() < 0.5:
later.append(x)
else:
now.append(x)
for x in pending + now:
inc = max(inc, x)
pending = later
for bb, k, nw, nv in children:
if bb > inc:
ident += 1
heapq.heappush(pool, (-bb, ident, k, nw, nv))
return inc, evals
R = 12
def report(name, runs):
counts = sorted(r[1] for r in runs)
optimum = sorted(set(r[0] for r in runs))
print(f"{name} {str(optimum):<7} {counts[0]:>6} to {counts[-1]:>6}"
f" {len(set(counts)):>4} over {R} runs")
d = serial(rng=random.Random(0))
print(f"D deterministic: {d[2]} pruning decisions were exact ties "
f"(bound = incumbent)")
print()
print(f"{'':25} {'optimum':<7} {'evaluations':<16} distinct")
report("D deterministic ",
[serial(rng=random.Random(i)) for i in range(R)])
report("N one ulp of bound noise",
[serial(noise=2.0 ** -52, rng=random.Random(i)) for i in range(R)])
report("T random tie-breaking ",
[serial(random_ties=True, rng=random.Random(i)) for i in range(R)])
report("B1 fixed batches, b=1024 ",
[batched(1024, rng=random.Random(i)) for i in range(R)])
report("B2 random batches, b=1024",
[batched(1024, randomized=True, rng=random.Random(i))
for i in range(R)])
D deterministic: 2435 pruning decisions were exact ties (bound = incumbent)
optimum evaluations distinct
D deterministic [991] 127505 to 127505 1 over 12 runs
N one ulp of bound noise [991] 128085 to 128156 12 over 12 runs
T random tie-breaking [991] 121597 to 121906 12 over 12 runs
B1 fixed batches, b=1024 [991] 137215 to 137215 1 over 12 runs
B2 random batches, b=1024 [991] 139986 to 143010 12 over 12 runs
(What the nondeterminism script shows) The answer is \(991\) in all sixty runs. The node count is a constant only in the two deterministic modes. One ulp of relative noise in the bound, which is what an order-dependent reduction produces, moves the count by about half a percent and makes every one of twelve runs different. The reason is that \(2{,}435\) pruning decisions in this integer-valued instance sit exactly on the threshold, and the noise lets some of them survive, with their subtrees. Random tie-breaking with exact arithmetic moves the count more, and downwards here, which is the acceleration anomaly of Theorem 6.2.4 at work. The fixed batch schedule is deterministic, and its count exceeds the serial count by the inflation that Section 6.4 analyses. Randomizing the batch composition and delaying some incumbents by one round adds a spread of three thousand nodes on top. The script runs in ten to fifteen seconds and is sequential. The matching measurement on a real machine is in Section 6.3, where the work-stealing pool processes between \(537{,}013\) and \(537{,}123\) nodes across runs of one instance.
(Correctness and reproducibility are separate requirements) The last proposition is the one that makes GPU design tractable. It separates the requirement that pruning be correct from the requirement that a run be reproducible. It concerns the safe bound of Theorem 7.3.1, which for a node LP \(\min\{c^\top x : Ax \ge b,\ l \le x \le u\}\) is \(b^\top y + \sum_j \min(r_j l_j, r_j u_j)\) at any multiplier vector \(y \ge 0\) on the rows (the \(\lambda\) of Section 2.2, written \(y\) for a linear program as in Proposition 2.2.4(b)), with \(r = c - A^\top y\) the reduced costs. That number is a valid lower bound on the node for every \(y\), feasible or not. It is display (2.2.1), and Section 7.3 gives its attribution and its floating-point form.
Proposition 6.2.21 (validity is order-independent under directed rounding; reproducibility is not). Let a node bound be computed as the safe bound of Theorem 7.3.1 with every operation rounded toward \(-\infty\). Then every reduction schedule yields a valid lower bound on the node, different schedules may yield different valid bounds, and pruning against any of them is correct. Under round-to-nearest no such guarantee holds: a schedule may overestimate the exact bound by up to the amount in Theorem 6.2.17 and prune wrongly when the margin to the incumbent is smaller than that.
Proof. Each \(\mathrm{fl}_{\downarrow}(p + q) \le p + q\), so by induction on any tree the computed sum is at most the exact sum. The products with the box ends are handled by taking the minimum over the enclosure, as Section 7.3 does. Different trees round at different nodes and give different values, all below the exact one. ∎
(What a GPU can offer today) Directed rounding is available on NVIDIA devices as per-operation intrinsics rather than as a global mode. It makes pruning correct whatever the hardware does with the order, at the cost of one intrinsic per operation. A fixed schedule makes the run reproducible, at the cost of the freedom described above. The first is a correctness requirement and the second a product requirement, and a design can take the first without the second.The intrinsics are __dadd_rd, __dadd_ru, __dmul_rd, __dmul_ru, __fma_rd, __fma_ru and their division and square-root relatives: NVIDIA, CUDA Math API, "Double Precision Intrinsics", read 5 October 2026. MAiNGO's documentation states the first half for its GPU interval bounder: "Due to the outward rounding, interval arithmetic will keep the true lower bound enclosed in the interval no matter which precision is used" (MAiNGO documentation, "Special uses of MAiNGO", read 5 October 2026). What a GPU can offer today, then, is three things: a single kernel that is bitwise reproducible on one device, a batched bound that is valid in every order, and, in one solver, an experimental deterministic mode for the CPU tree that excludes the GPU heuristics (Section 6.3). The CPU standard, a deterministic path on the same machine with the same thread count, is not yet available for any GPU branch and bound. The sentence "the same answer twice" has to be qualified accordingly.
Systems
This subsection describes the systems that run branch and bound in parallel today: the supercomputer frameworks that hold the records, the commercial solvers and their deterministic modes, and what is known about the global MINLP solvers. It then places the figure that simulates a shared pool, and closes with a compile-checked implementation of Algorithm 6.2.8. The purpose is to give the reader the design decisions each system made, with the numbers that justify them. The GPU designs of Section 6.4 can then be read as variations on the same decisions.
Frameworks and supercomputer records
(UG and its five mechanisms) The Ubiquity Generator (UG) of the Zuse Institute Berlin is an external parallelization framework: a load coordinator process and solver processes that each run an unmodified base solver on a subproblem. Its instantiation with SCIP over MPI, the message-passing interface used between the processes of a cluster, is ParaSCIP, and over threads it is FiberSCIP.Y. Shinano, S. Heinz, S. Vigerske and M. Winkler, "FiberSCIP — a shared memory parallelization of SCIP", INFORMS Journal on Computing 30 (2018), also ZIB-Report 13-55; Y. Shinano, "The Ubiquity Generator framework: 7 years of progress in parallelizing branch-and-bound", Operations Research Proceedings 2017 (Springer, 2018). The current survey of MILP parallelism is T. Ralphs, Y. Shinano, T. Berthold and T. Koch, "Parallel solvers for mixed integer linear optimization", in Handbook of Parallel Constraint Reasoning (Springer, 2018). Five mechanisms carry the design. The first is layered presolving, presolve being the simplification pass of Section 2.5 that tightens bounds and removes redundant rows and columns before a solve: the root is presolved once, and every subproblem is presolved again by its solver. The second is ramp-up, which UG offers in three forms. In normal ramp-up, workers that already hold a subproblem return every second child to the coordinator until its pool holds enough good nodes. In racing ramp-up, all workers solve the whole problem from the root with different settings until a winner is chosen and its open nodes are redistributed. In the self-split ramp-up of UG 1.0, every worker generates the same deterministic tree prefix independently and keeps the leaves assigned to it, with no communication at all. The third mechanism is dynamic load balancing by a collecting mode, described below. The fourth is checkpointing of the primitive nodes, those with no ancestor in the coordinator's pool, so that a run can be resumed across jobs. The fifth, in FiberSCIP only, is a deterministic mode. It is built on a deterministic clock that counts communication-point calls and on a token circulated among the solvers in a fixed order, and it is intended for debugging and reproducible experiments. Both ParaSCIP and FiberSCIP run opportunistically by default.The ramp-up strategies and the collecting mode are described in Y. Shinano, T. Achterberg, T. Berthold, S. Heinz, T. Koch and M. Winkler, "Solving open MIP instances with ParaSCIP on supercomputers using up to 80,000 cores", IPDPS 2016, also ZIB-Report 15-53, Section 2.2, and in ZIB-Report 13-66; self-split ramp-up and the statement that UG 1.0 is "completely different from the previous versions internally" are in K. Bestuzheva et al., "The SCIP Optimization Suite 8.0", arXiv 2112.08872 (2021), Section 9. A UG 1.0 beta with gap-limit support ships with SCIP 9 (S. Bolusani et al., arXiv 2402.17702 (2024), Section 7) and a UG pseudo-Boolean application with SCIP 10 (arXiv 2511.18580 (2025), Section 8), where in the 2024 Pseudo-Boolean Competition SCIP solved 759 of 1,207 instances and FiberSCIP 776.
Algorithm 6.3.1 (racing ramp-up followed by node redistribution). The pseudocode below is the racing ramp-up of UG and of Gurobi's distributed MIP.
Algorithm 6.3.1 Racing ramp-up followed by node redistribution
(UG; Gurobi distributed MIP)
Input problem; P workers; P parameter settings s_1, ..., s_P
(branching, heuristics, seeds); a racing limit tau in time or
in open nodes; a threshold p of good nodes
1. the coordinator sends the presolved root to all workers;
worker i starts a complete solve with s_i
2. every worker reports periodically (its tree's lower bound, its
open-node count, its incumbent); incumbents are broadcast to all
workers at once, since they prune in every tree
3. if some worker solves the problem during the race:
broadcast termination and return its solution
4. at the limit tau: choose the winner w by a score combining its
tree's lower bound and its open-node count
5. collect the open nodes of w; every other worker discards its tree
6. if w has fewer open nodes than there are workers:
switch to normal ramp-up (workers return every second child
until the coordinator holds p good nodes);
else distribute w's nodes to the idle workers
7. continue with dynamic load balancing
Invariant
the global lower bound is the minimum over all workers' tree
bounds, and the incumbent is global at all times, so correctness
does not depend on which tree wins.
Cost
P full solves for time tau, of which P - 1 are discarded
(redundant work (P - 1) tau); one transfer of w's frontier.
Parallel
the race is embarrassingly parallel and gives a portfolio effect
(Type 3 of Definition 6.1.2).
Racing ramp-up as in Algorithm 6.3.1, on P workers
coordinator: the presolved root
| limit tau
+--> worker 1, s_1 ====== tree 1 ======| discarded
+--> worker 2, s_2 ====== tree 2 ======| discarded
: |
+--> worker w, s_w ====== tree w ======| the winner
: |
+--> worker P, s_P ====== tree P ======| discarded
^ |
incumbents: broadcast to all | score: the tree's
workers at once | lower bound and
| open-node count
v
the open nodes of w, to the coordinator
|
v
fewer open nodes than workers?
| yes | no
v v
normal ramp-up: workers return w's nodes go to the
every second child until the idle workers
coordinator holds p good nodes |
| |
+-----------+------------+
v
dynamic load balancing
if some worker solves the problem during the race (step 3):
termination is broadcast and its solution is returned
(Racing ramp-up elsewhere, and the doubt about it) Gurobi's distributed MIP documents the same scheme under the same name: "each distributed worker solves the problem concurrently. At a certain stage, Gurobi selects the most promising branch-and-bound tree created by these workers (also known as 'racing ramp-up')", after which "the different machines explore different parts of the tree and synchronize their progress".Gurobi Optimizer Reference Manual, "Distributed algorithms", docs.gurobi.com, read 5 October 2026; the page adds that distributed MIP is "especially effective for models that generate large but shallow search trees". Koch, Ralphs and Shinano are sceptical, however. In their experience racing approaches "have not proven to be effective enough to make up for the reduction in subproblems solved per thread per second in the initial parts of the algorithm". The race converts idle time into a portfolio, and a portfolio is not the same as a wider tree.T. Koch, T. Ralphs and Y. Shinano, "Could we use a million cores to solve an integer program?", Mathematical Methods of Operations Research 76 (2012), the discussion of ramp-up.
Algorithm 6.3.2 (the load coordinator of UG/ParaSCIP, abridged). The pseudocode below gives the coordinator's loop, including the collecting mode and the checkpoint rule.
Algorithm 6.3.2 Load coordinator of UG/ParaSCIP, abridged
State the coordinator's node pool; per-solver status (open nodes,
subtree bound); global incumbent and dual bound
1. solvers run the base solver on a subproblem and send incumbents
and subtree bounds to the coordinator periodically
2. a solver that empties its subtree asks for a node; the
coordinator sends the best unassigned node in its pool
3. the coordinator counts the "good" nodes in its pool: nodes whose
bound is within a threshold of the global bound, relative to
max(|global bound|, 1)
4. if fewer than p good nodes are unassigned:
enter COLLECTING MODE: ask the solvers whose subtrees contain
good nodes to send back every second child they generate,
until p have been collected
5. leave collecting mode when p good nodes are held; re-enter
whenever the pool runs thin
6. checkpoint: write only the primitive nodes (nodes with no
ancestor in the pool) to disk; on restart, re-presolve the root
and redistribute the saved nodes
7. terminate when the pool is empty and every solver is idle
Cost
O(1) messages per transferred node; the coordinator becomes a
serial bottleneck as the transfer rate grows.
Parallel
the solvers of step 1 run independently; the coordinator is serial.
The UG coordinator of Algorithm 6.3.2 and one solver, over time
solver i: base solver load coordinator: node pool,
on a subproblem global incumbent, dual bound
| |
|-- incumbent, subtree bound ->| periodically
| |
|-- subtree empty: a node? --->|
|<---- best unassigned node ---|
| |
| | fewer than p good nodes
|<-- send every second child --| unassigned:
|-- a child ------------------>| COLLECTING MODE
|-- a child ------------------>|
| ... | p good nodes held:
| | collecting mode left
| |
| |--> disk: the primitive
| | nodes (checkpoint)
| |
pool empty and every solver idle: terminate
only solvers whose subtrees contain good nodes are asked to
send children; the mode is entered again whenever the pool
runs thin
(The records) The records set with this design span a decade. In 2010 ParaSCIP solved two open MIPLIB 2003 instances, ds and stp3d, on up to 2,048 cores of HLRN II, in more than ten consecutive checkpoint-restarted runs. By 2013 it ran on up to 35,200 cores of HLRN II, HLRN III and Titan. The 2016 paper reports twelve previously unsolved MIPLIB instances with up to 80,000 cores on Titan. The 2020 report brings the count to twenty-one and states that ParaSCIP "can stably handle over 40,000 cores". It also records that the instance rmine10 took 48 jobs over about 75 days and about 5,660 CPU-core-years on HLRN III and Titan. One of those jobs was the single 80,000-core run of 13 hours, which by itself spent \(80{,}000 \times 13 = 1{,}040{,}000\) core-hours, about \(119\) core-years.Shinano et al., ZIB-Report 10-27 (2010; printed in Competence in High Performance Computing 2010, Springer, 2012); ZIB-Report 13-66 (2013); IPDPS 2016; and Y. Shinano, T. Achterberg, T. Berthold, S. Heinz, T. Koch and M. Winkler, "Solving previously unsolved MIP instances with ParaSCIP on supercomputers by using up to 80,000 cores", ZIB-Report 20-16 (May 2020). The 2020 report presents its results as "solving open instances over a period of seven years", on several machines (among them HLRN II, HLRN III, ISM and Titan), and the 80,000 cores were used in one run. Other UG instantiations solved Steiner tree instances with SCIP-Jack on up to 43,000 cores and ran ParaXpress on up to 43,344 cores. A shortest-lattice-vector solver on UG ran on up to 100,032 cores, and a UG solver with a Lagrangian doubly nonnegative bound solved the quadratic assignment instances tai30a and sko42 for the first time.Y. Shinano, D. Rehfeldt and T. Koch, "Building optimal Steiner trees on supercomputers by using up to 43,000 cores", CPAIOR 2019, LNCS 11494; Y. Shinano, T. Berthold and S. Heinz, "ParaXpress: an experimental extension of the FICO Xpress-Optimizer to solve hard MIPs on supercomputers", Optimization Methods and Software 33 (2018); N. Tateiwa, Y. Shinano, S. Nakamura, A. Yoshida, S. Kaji, M. Yasuda and K. Fujisawa, "Massive parallelization for finding shortest lattice vectors based on Ubiquity Generator framework", SC20 (2020); K. Fujii, N. Ito, S. Kim, M. Kojima, Y. Shinano and K.-C. Toh, "Solving challenging large scale QAPs", arXiv 2101.09629 (2021).
ParaSCIP's records, and rmine10 as a chain of checkpointed jobs
2010 2013 2016 2020
-+----------------+----------------+----------------+---->
| | | |
up to 2,048 up to 35,200 up to 80,000 "stably"
cores of cores of HLRN cores on over 40,000
HLRN II II, HLRN III Titan cores
and Titan
ds and stp3d 12 previously 21 in all
(MIPLIB 2003) unsolved MIPLIB
instances
rmine10, on HLRN III and Titan: 48 jobs over about 75 days
[job 1]--ckpt-->[job 2]--ckpt--> ... --ckpt-->[job 48]
|
checkpoint: the primitive nodes to disk; the next
job re-presolves the root and redistributes them
the 48 jobs: about 5,660 CPU-core-years
one of them: 80,000 cores x 13 hours = 1,040,000 core-hours,
about 119 core-years
(What limits the scaling, and the older frameworks) Koch, Ralphs and Shinano (2012) ask what limits this scaling, and their findings are the ones a GPU designer should also know. The average speedup of commercial MIP solvers from one to twelve threads on MIPLIB 2010 was roughly a factor of three. A distributed LP for a half-terabyte root relaxation would cap a whole machine at about fifteen branch-and-bound nodes per second. Node files at scale would mean petabytes of storage. And the phases, not the arithmetic, dominate at tens of thousands of cores.Koch, Ralphs and Shinano (2012), cited in Section 6.1. Older frameworks built the same ideas on earlier machines. Eckstein's centrally controlled search on the CM-5 reached "near-linear speedups using 64–128 processors", and a randomized successor followed it. Bixby, Cook, Cox and Lee wrote a distributed code. PICO and its successor PEBBL were built at Sandia. The ALPS/CHiPPS library hierarchy came with a computational study, which found that "properties of the problem class itself can have a substantial effect on the efficiency". Bob++ and DryadOpt followed. The latter ran branch and bound on a data-parallel execution engine and reported linear scaling in the number of machines.J. Eckstein, "Parallel branch-and-bound algorithms for general mixed integer programming on the CM-5", SIAM Journal on Optimization 4 (1994); J. Eckstein, "Distributed versus centralized storage and control for parallel branch and bound: mixed integer programming on the CM-5", Computational Optimization and Applications 7 (1997); R. E. Bixby, W. Cook, A. Cox and E. K. Lee, "Computational experience with parallel mixed integer programming in a distributed environment", Annals of Operations Research 90 (1999); J. Eckstein, C. A. Phillips and W. E. Hart, "PICO: an object-oriented framework for parallel branch and bound", in Inherently Parallel Algorithms in Feasibility and Optimization and their Applications (Elsevier, 2001); J. Eckstein, W. E. Hart and C. A. Phillips, "PEBBL: an object-oriented framework for scalable parallel branch and bound", Mathematical Programming Computation 7 (2015); T. K. Ralphs, L. Ladányi and M. J. Saltzman, "A library hierarchy for implementing scalable parallel search algorithms", The Journal of Supercomputing 28 (2004); Y. Xu, T. K. Ralphs, L. Ladányi and M. J. Saltzman, "Computational experience with a software framework for parallel integer programming", INFORMS Journal on Computing 21 (2009); A. Djerrah, B. Le Cun, V.-D. Cung and C. Roucairol, "Bob++: framework for solving optimization problems with branch-and-bound methods", HPDC 2006; M. Budiu, D. Delling and R. F. Werneck, "DryadOpt: branch-and-bound on distributed data-parallel execution engines", IPDPS 2011.
Commercial solvers
(Gurobi, Xpress, CPLEX and SCIP) The commercial solvers differ from the frameworks in one respect above all. They are deterministic by default, or offer determinism as a documented mode, and they pay for it. Gurobi's documentation states that the solver "is designed to be deterministic" and "is also deterministic for parallel optimization: you always get the same results with multiple threads". The exceptions are the time-dependent parameters, the default LP method, which "can run non-deterministic concurrent in some cases", and the concurrent MIP mode, which runs several independent MIP solves with different settings and "is not deterministic". Distributed MIP "maintains deterministic behavior when using identical input parameters, model data, and machine configurations and if no time-dependent behavior is used", whereas distributed concurrent does not.Gurobi, "Is Gurobi deterministic?", Help Center article 360031636051, and the Reference Manual pages "Parameters" (ConcurrentMIP, Method, WorkLimit, Seed) and "Distributed algorithms", all read 5 October 2026. The support article on threads states that more threads help models that need many nodes and "generally will not help much" for models solved at or near the root: "Does using more threads make Gurobi faster?", article 360013419951. FICO Xpress follows "a partial information approach and separating the concepts of simultaneous tasks and independent threads". Its authors report "almost linear" scaling on the CPUs of the time, with a solution path that is deterministic in a fixed environment and "to a certain extent, thread-independent".Berthold, Farmer, Heinz and Perregaard (2018), cited in Section 6.2. CPLEX exposes the choice as a parameter. Its parallel mode switch takes the values opportunistic (\(-1\)), automatic (\(0\), the default: "let CPLEX decide whether to invoke deterministic or opportunistic search") and deterministic (\(1\)). The same page says that "by default, CPLEX applies as much parallelism as possible while still achieving deterministic results", so the default automatic setting yields deterministic results. It adds that the parallel barrier method is "only deterministic", and that in deterministic mode the order in which user callbacks run across threads can still vary.IBM ILOG CPLEX Optimization Studio 22.1.2, "Parallel mode switch" (CPXPARAM_Parallel, ParallelMode), read 5 October 2026. The default value is automatic, not deterministic; the page's own description of the default behaviour is the sentence quoted. SCIP itself is sequential at the tree level. Its concurrent mode runs several SCIP instances with different settings and exchanges domain reductions between them, and tree parallelism is delegated to UG.R. L. Gottwald, S. J. Maher and Y. Shinano, "Distributed domain propagation", SEA 2017, LIPIcs 75.
(cuOpt) NVIDIA's cuOpt, the one production MIP solver with GPU heuristics, keeps its branch and bound on the CPU and has rebuilt that part release by release. Its release notes record "parallel branch and bound on the CPU: multiple best-first search and diving threads" in 25.10. In 26.02 they record "new parallel reliability branching" and "experimental support for determinism in the parallel branch-and-bound solver", with the note that "GPU heuristics are not supported yet in this mode". In 26.06 the workers "maintain their own heaps and steal nodes from other workers", and by NVIDIA's own account "on average, this increases the number of nodes explored by 3x". Its C API exposes the switch between an opportunistic and a deterministic mode, and a work limit as the deterministic substitute for a time limit. The dated release table of Section 7.4 gives the dates of these releases.NVIDIA cuOpt User Guide, "Release Notes" 25.05 to 26.08 and "MIP C API" (CUOPT_MIP_DETERMINISM_MODE, CUOPT_MODE_OPPORTUNISTIC, CUOPT_MODE_DETERMINISTIC), and the source header constants.h (mip_determinism_mode, work_limit, cudss_deterministic), github.com/NVIDIA/cuopt, all read 5 October 2026. The 3x figure is the vendor's. No published number measures the cost of the deterministic mode.
Global MINLP solvers
(What is documented for the global solvers) For the global MINLP solvers the record is thinner. The statement the documentation supports is that for most of them no tree-level parallelism has been verified from primary sources. BARON, Couenne, ANTIGONE, LINDO and SHOT should be treated as unknown in this respect rather than as sequential. MAiNGO documents a parallel version "which allows the use of multiple processors via MPI", in a manager–worker arrangement. Its "main use is for large problems where the B&B algorithm needs a long time to converge". It is absent from the PyPI and Julia packages, and no scaling number for it was found.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 four-author citation the MAiNGO documentation asks for; MAiNGO documentation, "Special uses of MAiNGO" and "Installation", and the repository README, read 5 October 2026. Xpress Global runs a deterministic parallel branch and cut, deterministic by default, and the experiments of its paper use sixteen threads. That paper contains the only table of thread speedups for a global MINLP solver that this post could find. The speedups relative to one thread are \(1.34, 1.72, 2.19, 2.59\) and \(2.59\) at \(2, 4, 8, 16\) and \(32\) threads over all solvable instances of its 1,308-instance set, and \(1.54, 2.21, 3.13, 3.97\) and \(4.05\) on the instances that needed at least a hundred nodes. The authors remark that "going from 16 to 32 threads, the number of solved instances no longer increases but actually goes down slightly".Belotti, Berthold, Gally, Gottwald and Pólik (2025), Table 7 and Section 4.9; Section 3 gives the setup, "using 16 threads in deterministic parallel mode", and Section 4.9 states that "the deterministic mode is the default". The set is all MINLPLib instances solvable by at least three solvers, an internal customer set, and QPLIB and MINLPLib instances solvable by earlier Xpress versions, each in three permutations. The other published thread-count evidence is a table of solved counts rather than of speedups: the SCIP 8 paper's thread table, quoted in Section 5.3, in which BARON solved 161, 160, 160 and 158 of the 200 unpermuted instances at 1, 4, 8 and 16 threads and FiberSCIP 161, 145, 147 and 152. More threads did not mean more solved instances for either code. FiberSCIP's own paper demonstrates "the current performance of FiberSCIP for solving mixed-integer linear programs (MIPs) and mixed-integer nonlinear programs (MINLPs) in parallel", and is the other published parallel MINLP result on shared memory. Octeract, whose record and freeze Section 5.3 reports, describes itself as a distributed engine. No benchmark tests the claim.The distributed-execution statement is on the vendor's product page, octeract.com/octeract-engine, read 5 October 2026, and is a vendor claim. Interval branch and bound has its own parallel literature. Casado and coauthors parallelized an interval algorithm on shared memory and found the per-thread memory allocator to be the bottleneck. Berenguel and coauthors estimate the work left in an interval search from the pruning observed so far, which is what a scheduler needs in order to decide how many processors to add.L. G. Casado, J. A. Martínez, I. García and E. M. T. Hendrix, "Branch-and-bound interval global optimization on shared memory multiprocessors", Optimization Methods and Software 23 (2008); J. L. Berenguel, L. G. Casado, I. García and E. M. T. Hendrix, "On estimating workload in interval branch-and-bound global optimization algorithms", Journal of Global Optimization 56 (2013).
| system | pool and balancing | scale on record | deterministic mode |
|---|---|---|---|
| UG / ParaSCIP | load coordinator; normal, racing or self-split ramp-up; collecting mode; checkpoints of primitive nodes | up to 80,000 cores in one run; "stably" over 40,000; 21 open MIPLIB instances solved, 2010–20 | no (opportunistic by design) |
| FiberSCIP | the same over threads | one node; 161, 145, 147 and 152 of the 200 unpermuted instances solved at 1, 4, 8 and 16 threads, in the SCIP 8 paper's thread table (quoted in Section 5.3) | yes, via a deterministic clock and a token; off by default |
| Gurobi | threads on one machine; distributed MIP with racing ramp-up | one machine, or distributed jobs | yes by default, per its documentation; exceptions: time-dependent parameters, the default LP method (Method = -1) in some cases, concurrent MIP, distributed concurrent; distributed MIP only with identical input parameters, model data and machine configurations and no time-dependent behavior |
| FICO Xpress (Xpress Global) | task-based tree search | 16 threads in its experiments | yes by default; 2–9 % saved by the opportunistic variant (Xpress Global paper, Table 8) |
| CPLEX | threads; distributed parallel MIP | one machine, or distributed | switch: opportunistic, automatic (default; the page documents deterministic results by default), deterministic; ticks |
| SCIP 10 | concurrent solvers sharing incumbents and domain reductions; tree parallelism via UG | one machine | tree search sequential; tree parallelism delegated to UG |
| cuOpt 26.08 | CPU workers with their own heaps and work stealing; GPU heuristics share a pool | one machine (several GPUs for PDLP, not for the tree) | experimental, from 26.02; excludes the GPU heuristics; work limit |
| MAiNGO | manager–worker MPI | not reported | not documented |
| Para-B&B (HiGHS) | replicated solver state, no nondeterministic primitives | 8 threads: 2.17x geometric mean over 80 MIPLIB 2017 instances, 5.12x on node-heavy ones, as its authors report | yes, the only mode; 34.7 % idle on average, as reported |
| N2N (SCIP) | distributed tasks | 1,000 MPI processes: 22.5x and 12.7x on two clusters, in its nondeterministic mode, as its authors report | sliding-window deterministic mode; no number given for it in the abstract |
The shared pool, simulated
The figure runs the shared-pool model on a knapsack tree of fifteen items and about a thousand nodes, the instance that the batch figure of Section 6.4 reuses. \(P\) workers draw from one pool, by best bound or depth-first. The incumbent is shared the moment it is found or kept private to each worker. Each node costs between one and two time units, with the cost attached to the node rather than to the worker, so that the same node costs the same in every run. A node pruned at the moment it is popped costs a tenth of a unit. The lower panel is a Gantt chart: one row per worker, one bar per node along a time axis, white where the worker idles. The phases of Definition 6.2.11 are the white margin at the left and the ragged edge at the right. In the default view eight workers share every incumbent and take the open node with the best bound. They finish in \(173.1\) time units against \(1{,}360.8\) for one worker, a speedup of \(7.86\) and an efficiency of \(98\,\%\), processing \(1{,}063\) nodes where the sequential run needed \(1{,}057\). The first \(4.6\) time units are ramp-up. The readout's table is the speedup curve. Sharing incumbents, the speedups are \(2.00, 3.98, 7.86, 15.22, 28.03\) and \(42.63\) for \(P = 2, 4, 8, 16, 32, 64\), with the node count rising from \(1{,}057\) to \(1{,}197\). With private incumbents they are \(1.95, 3.73, 6.75, 11.70, 18.04\) and \(26.51\), with the node count rising to \(2{,}707\) at \(P = 64\). Watch the speedup curve bend away from the diagonal as \(P\) grows, which is ramp-up and ramp-down on a tree of a thousand nodes. It bends further when the incumbents are private, because each worker then explores what the others' solutions would have pruned.
A work-stealing pool in C++23
The listing implements Algorithm 6.2.8 on a fifty-item knapsack with Dantzig's bound. Each node is packed into 64 bits (the next item, the weight so far, the value so far), which is also the representation one would use for a device-resident frontier. Each thread owns a Chase–Lev deque, and the incumbent is one shared atomic integer. Termination follows Proposition 6.2.13, with the per-thread count flushed to the global counter whenever a thread's deque runs dry and after every steal; the listing's batch flush, at a count of \(-1{,}024\), never fires, since a thread's count can fall only by the one or two nodes its deque held at its last flush. The instance maximizes.
The listing's data: a node in one 64-bit word, a deque per worker
63 56 55 28 27 0
+--------+------------------------+------------------------+
| k | weight so far | value so far |
+--------+------------------------+------------------------+
8 bits 28 bits 28 bits
top bottom
| |
v v
+--------+--------+--------+--------+--------+--------+
| oldest | | ... | | newest | |
+--------+--------+--------+--------+--------+--------+
| | ^
v v |
steal (a thief): pop (owner): push
the oldest node, the newest (owner)
the largest subtree node, a dive
k: the next item to decide, 0 at the root, 50 at a leaf
the deque: 2^20 slots, used circularly
// Shared-pool parallel branch and bound (C++23).
//
// Algorithm 6.2.8 on a 0/1 knapsack of fifty items with Dantzig's
// bound: one Chase–Lev deque per thread, random work stealing, a
// shared atomic incumbent, and a count of open nodes flushed in
// batches for termination. The optional argument is the number of
// threads P, 8 by default.
#include <algorithm>
#include <atomic>
#include <chrono>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <random>
#include <thread>
#include <vector>
// items, sorted by value per unit of weight
constexpr int N = 50;
static int W[N], V[N], CAP;
// A node packed into 64 bits:
// next item (8 bits) | weight so far (28) | value so far (28)
using Node = uint64_t;
Node pack(int k, int wt, int val) {
return Node(k) << 56 | Node(wt) << 28 | Node(val);
}
int k_of(Node x) { return int(x >> 56); }
int wt_of(Node x) { return int(x >> 28 & 0xFFFFFFF); }
int val_of(Node x) { return int(x & 0xFFFFFFF); }
// One worker's deque: Chase and Lev (2005); fences as in Le, Pop, Cohen
// and Zappa Nardelli (2013).
struct alignas(128) Deque {
// 2^20 slots, used circularly (index mod CAPB)
static constexpr int64_t CAPB = 1 << 20;
// thieves advance top: keep it off the owner's cache line
alignas(64) std::atomic<int64_t> top{0};
alignas(64) std::atomic<int64_t> bottom{0};
std::vector<std::atomic<Node>> buf =
std::vector<std::atomic<Node>>(CAPB);
// owner only, at the bottom
bool push(Node x) {
int64_t b = bottom.load(std::memory_order_relaxed);
int64_t t = top.load(std::memory_order_acquire);
// full: the caller expands the node itself
if (b - t >= CAPB) return false;
buf[b % CAPB].store(x, std::memory_order_relaxed);
std::atomic_thread_fence(std::memory_order_release);
bottom.store(b + 1, std::memory_order_relaxed);
return true;
}
// owner only, from the bottom: the newest node, a dive
bool pop(Node& x) {
int64_t b = bottom.load(std::memory_order_relaxed) - 1;
bottom.store(b, std::memory_order_relaxed);
std::atomic_thread_fence(std::memory_order_seq_cst);
int64_t t = top.load(std::memory_order_relaxed);
if (t > b) {
bottom.store(b + 1, std::memory_order_relaxed);
return false;
}
x = buf[b % CAPB].load(std::memory_order_relaxed);
// more than one element left: no race with a thief
if (t != b) return true;
bool won = top.compare_exchange_strong(t, t + 1,
std::memory_order_seq_cst,
std::memory_order_relaxed);
bottom.store(b + 1, std::memory_order_relaxed);
return won;
}
// any thief, from the top: the oldest node, the largest subtree
bool steal(Node& x) {
int64_t t = top.load(std::memory_order_acquire);
std::atomic_thread_fence(std::memory_order_seq_cst);
if (t >= bottom.load(std::memory_order_acquire)) return false;
x = buf[t % CAPB].load(std::memory_order_relaxed);
return top.compare_exchange_strong(t, t + 1,
std::memory_order_seq_cst,
std::memory_order_relaxed);
}
};
// Dantzig: fill greedily from item k, then the fitting fraction
double bound(int k, int wt, int val) {
for (int i = k; i < N; ++i) {
if (wt + W[i] <= CAP) {
wt += W[i];
val += V[i];
} else {
return val + (CAP - wt) * double(V[i]) / W[i];
}
}
return val;
}
// read at every node, written rarely
std::atomic<int> incumbent{0};
// pushed minus finished, flushed in batches; 0 means the search is over
alignas(64) std::atomic<int64_t> open_nodes{0};
// nodes processed, summed at exit
alignas(64) std::atomic<int64_t> processed{0};
// The loop every worker runs: Algorithm 6.2.8, steps 2 to 8.
void worker(int id, std::vector<Deque>& dq, int P) {
std::mt19937 rng(id * 7919 + 1);
std::uniform_int_distribution<int> pick(0, P - 1);
Node x;
int64_t local = 0;
// delta: my pushes minus my finished nodes since the last flush
int64_t delta = 0;
auto flush = [&] {
if (delta) {
open_nodes.fetch_add(delta, std::memory_order_acq_rel);
delta = 0;
}
};
// Process one node. Invariant: an optimal leaf lies below an open
// node, or the incumbent is optimal.
auto process = [&](auto&& self, Node nd) -> void {
int k = k_of(nd), wt = wt_of(nd), val = val_of(nd);
++local;
if (bound(k, wt, val) <=
incumbent.load(std::memory_order_relaxed)) {
return; // pruned
}
// a leaf (every item decided): incumbent <- max(incumbent, val)
if (k == N) {
int cur = incumbent.load();
while (val > cur &&
!incumbent.compare_exchange_weak(cur, val)) {
// a failed exchange reloads cur; retry while val is better
}
return;
}
// push the 0-branch first so the owner pops the 1-branch next
for (int take : {0, 1}) {
if (take && wt + W[k] > CAP) continue; // item k does not fit
Node c = pack(k + 1, wt + take * W[k], val + take * V[k]);
if (bound(k + 1, wt_of(c), val_of(c)) <=
incumbent.load(std::memory_order_relaxed)) {
continue; // the child cannot beat the incumbent
}
++delta;
// deque full: expand the child at once
if (!dq[id].push(c)) {
self(self, c);
--delta;
}
}
};
while (true) {
if (dq[id].pop(x)) {
process(process, x);
if (--delta <= -1024) flush();
continue;
}
// my deque is empty: publish my count, then read the global one
flush();
if (open_nodes.load(std::memory_order_acquire) == 0) break;
int victim = pick(rng);
if (victim != id && dq[victim].steal(x)) {
process(process, x);
--delta;
flush();
}
}
processed.fetch_add(local);
}
int main(int argc, char** argv) {
int P = argc > 1 ? std::atoi(argv[1]) : 8;
std::mt19937 g(11);
std::vector<int> idx(N);
// strongly correlated items; the capacity is half the total weight
for (int i = 0; i < N; ++i) {
W[i] = 10 + g() % 51;
V[i] = W[i] + 10;
CAP += W[i];
idx[i] = i;
}
CAP /= 2;
// sort the items by value per unit of weight, best first
std::sort(idx.begin(), idx.end(), [](int a, int b) {
return double(V[a]) / W[a] > double(V[b]) / W[b];
});
int W2[N], V2[N];
for (int i = 0; i < N; ++i) {
W2[i] = W[idx[i]];
V2[i] = V[idx[i]];
}
std::copy(W2, W2 + N, W);
std::copy(V2, V2 + N, V);
// the root: one open node, in the deque of worker 0
std::vector<Deque> dq(P);
open_nodes = 1;
dq[0].push(pack(0, 0, 0));
auto t0 = std::chrono::steady_clock::now();
{
std::vector<std::jthread> th;
for (int i = 0; i < P; ++i) {
th.emplace_back(worker, i, std::ref(dq), P);
}
} // the jthreads join here
auto t1 = std::chrono::steady_clock::now();
double ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
std::printf("P=%d optimum=%d nodes=%lld time=%.1f ms\n", P,
incumbent.load(), (long long)processed.load(), ms);
}
One worker's loop in the listing, with its count of open nodes
+--> pop the bottom of my deque --ok--> process it; --delta;
| | failed flush if delta <= -1024
| v |
| flush: open_nodes += delta; delta = 0 |
| | |
| open_nodes == 0 ? --yes--> leave the loop; |
| | no processed += |
| v my node count |
| steal the top of a random victim's deque |
| | won: process it; --delta; flush |
| | lost, or the victim is me: nothing |
| v |
+-------+-----------------------------------------+
process: pruned if its bound <= incumbent; at a leaf the
incumbent is raised by compare-and-swap; otherwise each child
that fits and whose bound beats the incumbent is pushed, ++delta
(The measurements, and two engineering facts) The timings that follow are this post's own measurements, on a ten-core laptop with eight performance cores and a build with clang++ -std=c++23 -O2 -pthread, and they vary with the machine's load. On quiet runs the code processes \(537{,}013\) nodes sequentially in \(5.8\) to \(6.0\) milliseconds, about \(11\) nanoseconds per node. This search is as fine-grained as the flow-shop bounds of the GPU literature and far finer than any MINLP node. Eight threads take \(1.0\) to \(1.6\) milliseconds on the best of those runs, a speedup of \(3.6\) to \(6\), and \(3\) to \(6\) milliseconds on others. A re-run on the same laptop while it was loaded gave \(9\) to \(22\) milliseconds at one thread and \(1.1\) to \(1.3\) milliseconds at eight. The node count varies between runs, from \(537{,}013\) to \(537{,}123\) over all the runs recorded for this post, because the incumbent arrives at different moments. That is Theorem 6.2.4 in practice. Two engineering facts from writing the code are worth more than the timings. A first version did a read-modify-write on the global counter at every node and had the deques' top and bottom on shared cache lines, a cache line being the 64-byte unit in which a core fetches memory, so that two words on one line are fetched and invalidated together and a write by one core to either word forces every other core to refetch the line. It ran seven to nine times slower with eight threads than with one. Aligning the atomics and flushing the count in batches gave the numbers above. At eleven nanoseconds of work per node, the cost of synchronization exceeds the cost of the work. The design rule is to compare the cost of a synchronization with the cost of a node: a deque operation costs tens of nanoseconds, a knapsack node tens of nanoseconds, an LP node milliseconds. The second fact is that a run occasionally shows no speedup at all. The root is pushed to one deque, and the thieves' early attempts fail until that worker has built a frontier, so ramp-up is itself a random variable. A production design seeds every deque before the threads start, which is the breadth-first prefix of the Chapel codes and the self-split ramp-up of UG 1.0. Everything in the listing parallelizes except the first node. The only shared writes are the rare incumbent updates and the batched counter flushes.
Shared words in the first version and in the listing, 8 threads
first version the listing
top and bottom on one line alignas(64): a line each
+-----------------------+ +-----------------------+
| top bottom | | top | thieves
+-----------------------+ +-----------------------+
thieves advance top, the | bottom | owner
owner moves bottom +-----------------------+
the global counter: a read- open_nodes, flushed in batches:
modify-write at every node when the deque runs dry, or
once delta <= -1024
8 threads: seven to nine on quiet runs:
times slower than 1 thread 1 thread 5.8 to 6.0 ms
8 threads 1.0 to 1.6 ms on
the best of them,
3 to 6 ms on others
Ramp-up in the listing, and every deque seeded first (P = 8)
the listing: the root is pushed to one deque
deque 0 [ root ] worker 0 builds a frontier from it
deque 1 [ ] \
deque 2 [ ] | workers 1 to 7 try to steal: their
: | early attempts fail until that
deque 7 [ ] / frontier exists
a production design: every deque seeded before the threads start
deque 0 [ nodes of the prefix ]
deque 1 [ nodes of the prefix ]
:
deque 7 [ nodes of the prefix ]
ramp-up in the listing is a random variable; seeding every
deque first is the breadth-first prefix of the Chapel codes and
the self-split ramp-up of UG 1.0
GPU branch and bound
This subsection is about running the tree itself on a device. A GPU executes threads in lockstep groups, of thirty-two on NVIDIA devices, called warps, and one thread of such a group is a lane. A thread that takes a different branch of a conditional stalls its group. A kernel launch costs a fixed latency whatever the kernel does, of the order of microseconds on current devices, and a pointer followed into memory costs more than hundreds of arithmetic operations (Definition 7.1.1 gives the model and quotes no figure). A branch-and-bound node as a CPU solver represents it, with its own LP basis, its own cut pool and its own child pointers, is the opposite of what the device wants. The two designs that have worked, and the records they hold, are both ways of giving device threads constant-shape work with no pointers to follow. We describe them, analyse the one quantity that every batched design trades away, place the figure that measures it, and end with what is regular and what is irregular in a spatial branch-and-bound node. That last question is the design question for the GPU programme of Section 7.8.
Two designs
(The device-resident design) The device-resident design puts the whole tree on the GPU as tens of thousands of independent depth-first explorers. Gmys and coauthors introduced the data structure that makes this possible for permutation problems, the integer-vector-matrix (IVM). It is a constant-memory encoding of a depth-first position in a permutation tree, so that a device thread can own a search without a linked list. Work is stolen inside the GPU and, in the multi-GPU version, between GPUs through a coordinator.J. Gmys, M. Mezmaz, N. Melab and D. Tuyttens, "A GPU-based branch-and-bound algorithm using integer–vector–matrix data structure", Parallel Computing 59 (2016); J. Gmys, R. Leroy, M. Mezmaz, N. Melab and D. Tuyttens, "Work stealing with private integer–vector–matrix data structure for multi-core branch-and-bound algorithms", Concurrency and Computation: Practice and Experience 28 (2016); J. Gmys, M. Mezmaz, N. Melab and D. Tuyttens, "IVM-based parallel branch-and-bound using hierarchical work stealing on multi-GPU systems", Concurrency and Computation: Practice and Experience 29 (2017), reporting "near-linear speed-up on up to four GPUs" for flow-shop. The record result is Gmys's 2022 solution of permutation flow-shop instances open since Taillard posed them in 1993. Six design facts are verified from the paper. There are \(K = 16{,}384\) explorers per GPU. Work units are intervals \([a_i, b_i)\) of the factoradic numbering of the \(n!\) leaves, so that a thief steals the right half of a victim's interval. The factoradic numbering is the mixed-radix numbering of permutations whose \(k\)-th digit runs from \(0\) to \(k\). Victims inside a GPU are chosen on a 4-ary 7-cube of the explorers, the graph on \(4^7 = 16{,}384\) vertices in which two explorers are adjacent when their base-4 labels differ by one, cyclically, in exactly one digit. Checkpoint messages carry lists of up to 16,384 intervals. Incumbents are supplied by CPU heuristic threads, because depth-first search alone finds poor solutions (Theorem 6.2.6 in the other direction). And the coordinator is "active nearly 100 % of the time", that is, saturated, at 32 to 128 GPUs on the smaller instances.
The device-resident design: one GPU of Gmys's solver (2022)
K = 16,384 depth-first explorers per GPU; explorer i owns an
interval [a_i, b_i) of the factoradic numbering of the n! leaves
explorer 1 explorer 2 explorer 3 . . . explorer K
[a_1, b_1) [a_2, b_2) [a_3, b_3) [a_K, b_K)
a steal inside the GPU: the thief takes the right half
victim, before [a_i --------------+-------------- b_i)
victim, after [a_i --------------)
thief, after [-------------- b_i)
victims are chosen on a 4-ary 7-cube of the explorers:
4^7 = 16,384 labels of 7 base-4 digits; two explorers are
adjacent when their labels differ by one, cyclically, in
exactly one digit
a label: d d d d d d d a digit: 0 -- 1 -- 2 -- 3
| | |
one digit moves by one +--------------+
CPU heuristic threads ------- incumbents -------> the explorers
several GPUs: GPU <-- work --> coordinator <-- work --> GPU
(Gmys's results) The paper reports three results. Eleven of the twenty-three open Taillard instances were solved and eight upper bounds improved. The scaling experiments ran on up to 384 V100 GPUs and 3,840 CPU cores. Parallel efficiency was near \(90\,\%\) on 16, 64 and 128 GPUs for instances needing one, four and twenty-seven single-GPU hours, and the speedup reached \(215\) on 384 GPUs for the twenty-seven-hour instance. The largest proof, Ta058, took 256 GPUs for 13 hours 17 minutes: a tree of \(339 \times 10^{12}\) nodes, 3,399 GPU-hours, and about 64 CPU-years of sequential work by the paper's own estimate.J. Gmys, "Exactly solving hard permutation flowshop scheduling problems on peta-scale GPU-accelerated supercomputers", INFORMS Journal on Computing 34 (2022); the instances are from E. Taillard, "Benchmarks for basic scheduling problems", European Journal of Operational Research 64 (1993). The companion sequential engine is J. Gmys, M. Mezmaz, N. Melab and D. Tuyttens, "A computationally efficient branch-and-bound algorithm for the permutation flow-shop scheduling problem", European Journal of Operational Research 284 (2020).
(The arithmetic of the record run) The arithmetic on those numbers is instructive. The Ta058 run averaged \(339 \times 10^{12} / (13 \times 3600 \times 256) \approx 28\) million node decompositions per second per GPU. One V100 in isolation reached \(37.6\) million on a smaller instance, so the 256-GPU run kept about three quarters of single-device throughput over thirteen hours. And 64 CPU-years is \(561{,}024\) core-hours, which over 13 hours is \(43{,}156\) sequential-core equivalents, about \(169\) per GPU.
(The pool-based design) The pool-based design keeps the tree on the host and sends the regular work to the device. Chakroun, Melab and coauthors introduced it for flow-shop: CPU threads manage pools of nodes and the GPU evaluates the bounds of a pool.N. Melab, I. Chakroun, M. Mezmaz and D. Tuyttens, "A GPU-accelerated branch-and-bound algorithm for the flow-shop scheduling problem", IEEE CLUSTER 2012; I. Chakroun, N. Melab, M. Mezmaz and D. Tuyttens, "Combining multi-core and GPU computing for solving combinatorial optimization problems", Journal of Parallel and Distributed Computing 73 (2013); I. Chakroun and N. Melab, "Towards a heterogeneous and adaptive parallel branch-and-bound algorithm", Journal of Computer and System Sciences 81 (2015); T.-T. Vu and B. Derbel, "Parallel branch-and-bound in multi-core multi-CPU multi-GPU heterogeneous environments", Future Generation Computer Systems 56 (2016). The earliest GPU branch-and-bound papers are T. Carneiro, A. E. Muritiba, M. Negreiros and G. A. L. de Campos, "A new parallel schema for branch-and-bound algorithms using GPGPU", SBAC-PAD 2011, and, for the knapsack, A. Boukedjar, M. E. Lalami and D. El Baz, "Parallel branch and bound on a CPU-GPU system", PDP 2012, and M. E. Lalami and D. El Baz, "GPU implementation of the branch and bound method for knapsack problems", IPDPS Workshops 2012. The current form is the Chapel code of Helbecque, Krishnasamy, Carneiro, Melab and Bouvry. It keeps one host pool per CPU thread, with work stealing between pools and one or more GPUs per pool. A breadth-first prefix search fills the pools before the device work starts. A chunk of \(q = \min(|\text{pool}|, M)\) nodes is offloaded only once a pool holds at least \(m\) nodes. These are the two thresholds \(m \le M\) of Algorithm 6.4.3 below, and the paragraph after that algorithm says what they trade. The paper reports strong-scaling efficiency, the efficiency \(E(P)\) of Definition 6.1.1 on a fixed instance. On a node of eight NVIDIA A100 GPUs it is up to \(63\,\%\), and on a node of eight AMD MI50 GPUs about \(75\,\%\) on average. The instances are 20×20 Taillard instances whose trees have \(9.5 \times 10^6\) to \(3.9 \times 10^9\) nodes. On the smaller instances \(77\,\%\) on two A100s falls to \(50\,\%\) on eight. Their 2025 paper takes the same code to 1,024 GPUs of the LUMI supercomputer. For N-Queens the speedup relative to one eight-GPU node is \(4.02\) at 8 nodes, \(8.00\) at 16, \(15.81\) at 32, \(31.12\) at 64 and \(56.30\) at 128 nodes, an efficiency of \(44\,\%\) at 1,024 GPUs. Within a node eight GPUs give \(6.26\times\) on N-Queens and \(5.15\times\) on flow-shop. The 20×20 flow-shop instances stop scaling beyond 16 nodes because their trees are too small for the machine, which is Proposition 6.2.12 at scale.G. Helbecque, E. Krishnasamy, T. Carneiro, N. Melab and P. Bouvry, "A Chapel-based multi-GPU branch-and-bound algorithm", Euro-Par 2024 Workshops, LNCS 15385 (2025), whose exact wording is "a strong scaling efficiency of up to 63 % and 75 % on average using a GPU-powered processing node including 8 NVIDIA A100 devices and AMD MI50 GPUs, respectively": two platforms, not a range. G. Helbecque, E. Krishnasamy, T. Carneiro, N. Melab and P. Bouvry, "Portable PGAS-based GPU-accelerated branch-and-bound algorithms at scale", Concurrency and Computation: Practice and Experience 37 (2025), e70321, is the reference for GPU branch and bound at scale; its instances are N-Queens and permutation flow-shop. It is distinct from the same group's CPU-only Chapel paper, G. Helbecque, J. Gmys, N. Melab, T. Carneiro and P. Bouvry, "Parallel distributed productivity-aware tree-search using Chapel", Concurrency and Computation: Practice and Experience 35 (2023), which makes no GPU claim. The portability study across up to 512 NVIDIA and 1,024 AMD GPUs is T. Carneiro, E. Kayraklioglu, G. Helbecque and N. Melab, "Investigating portability in Chapel for tree-based optimization on GPU-powered clusters", Euro-Par 2024; the earlier CPU-cluster paper is T. Carneiro, J. Gmys, N. Melab and D. Tuyttens, "Towards ultra-scale branch-and-bound using a high-productivity language", Future Generation Computer Systems 105 (2020), with "up to 84 % of the linear speedup" on 1,024 cores.
The pool-based design: the Chapel code of Helbecque et al.
a breadth-first prefix search fills the pools; then
<=============== work stealing between the pools ===============>
^ ^ ^
v v v
+----------+ +----------+ +----------+
| pool of | | pool of | | pool of |
| CPU | | CPU | | CPU | . . .
| thread 1 | | thread 2 | | thread 3 |
+----------+ +----------+ +----------+
| ^ | ^ | ^
| q | bounds | q | bounds | q | bounds
v | v | v |
+----------+ +----------+ +----------+
| GPU(s) | | GPU(s) | | GPU(s) | . . .
+----------+ +----------+ +----------+
q = min(|pool|, M) nodes go down only once the pool holds at
least m nodes: the thresholds m <= M of Algorithm 6.4.3
| year | who | problem | design | result |
|---|---|---|---|---|
| 2012 | Boukedjar, Lalami and El Baz | 0/1 knapsack | list processed on CPU or GPU by its size | one Tesla C2050; "substantial speedup" |
| 2013 | Chakroun, Melab et al. | permutation flow-shop | host pools; GPU bounds a pool | multi-core plus GPU |
| 2016 | Gmys et al. | flow-shop | IVM: the tree on the device; stealing inside the GPU | one GPU; the device-resident design |
| 2017 | Gmys et al. | flow-shop | IVM with hierarchical stealing across GPUs | near-linear on 4 GPUs |
| 2022 | Gmys | flow-shop (Taillard) | 16,384 explorers per GPU; factoradic intervals; CPU heuristics for incumbents | 11 of 23 open instances; Ta058 on 256 GPUs in 13 h, \(339 \times 10^{12}\) nodes |
| 2025 | Helbecque et al. (LNCS 15385) | flow-shop 20×20 | multi-pool Chapel; \(m\) and \(M\) thresholds; breadth-first prefix | up to 63 % on 8 A100s, about 75 % on 8 MI50s |
| 2025 | Helbecque et al. (CCPE 37) | N-Queens, flow-shop | the same, on LUMI | 44 % efficiency at 1,024 GPUs; flow-shop stops scaling at 16 nodes |
| 2015 | Xi, Li and Xi | MIQP for hybrid MPC | B&B on the GPU; several branches at once ("multi-point radiation") | dual neural network as the QP solver |
| 2025 | Verheijen, Elkady, Lazar, Goswami | mixed-binary MPC | batches of branches sized to memory; intermediate nodes skipped | four-tank system; "promising" against Gurobi |
| 2020–2022 | Bunel et al.; Xu et al.; Wang et al.; Zhang et al. | neural-network verification | thousands of subdomains per batch; bound propagation as the bound | "up to three orders of magnitude faster than LP-based BaB" |
| 2026 | Blin, Gualandi, Maes, Lodi, Stellato | MILP strong branching and bound tightening | batched PDHG: \(K\) LPs sharing \(A\) as one matrix-matrix iteration | sizes where first-order beats simplex identified |
| 2026 | Guan, Luo, Li, Chen and Xu (B³-PWL) | piecewise-linear MIPs with SOS2 | batched first-order LP engine inside a GPU branch and bound | 9.25× geometric mean over cuOpt on 43 instances |
| 2026 | Liu and Lodi | k-sparse GLMs | padded node tensors kept on the device; batches of nodes | one to two orders of magnitude at zero gap |
| 2026 | Lucas, Meng and Mazumder | sparse regression | perspective nodes by ADMM, coordinates fine-grained, nodes coarse-grained | beats a CPU B&B and commercial solvers |
| 2025 | Zhang et al. (MAiNGO) | factorable NLP/MINLP | Type 1: the node's interval bound on \(m^d\) sub-boxes, mean value form | three orders of magnitude over CPU intervals without partition |
| 2026 | Gottlieb, Xu and Stuber (ParBB) | factorable NLP | pointwise McCormick kernels generated from the expressions | 9 ns against 237 ns per evaluation; 11–22× over EAGO (authors' numbers) |
(What the table shows) Every entry in the result column is the authors' own report, and none has been reproduced independently for this post. The rows from Xi, Li and Xi (2015) downward are the frontier that bears on MINLP, and Section 7 takes each of them up in detail. Here they stand as evidence for one claim. Every GPU branch and bound to date either keeps the tree on the device in a constant-shape encoding or sends batches of same-shape relaxations to it. Every one exploits a fixed matrix shared by all nodes or a closed-form bound.W. Xi, D. Li and Y. Xi, "A mixed-integer quadratic programming solver based on GPU", 34th Chinese Control Conference (2015); P. C. N. Verheijen, M. Elkady, M. Lazar and D. Goswami, "A GPU-aware batched branch and bound method for solving mixed-binary MPC problems", IEEE Conference on Control Technology and Applications (2025); R. Bunel, J. Lu, I. Turkaslan, P. H. S. Torr, P. Kohli and M. P. Kumar, "Branch and bound for piecewise linear neural network verification", Journal of Machine Learning Research 21 (2020); K. Xu et al., "Fast and complete: enabling complete neural network verification with rapid and massively parallel incomplete verifiers", ICLR 2021; S. Wang et al., "Beta-CROWN: efficient bound propagation with per-neuron split constraints for complete and incomplete neural network robustness verification", NeurIPS 2021; H. Zhang et al., "General cutting planes for bound-propagation-based neural network verification", NeurIPS 2022; N. Blin, S. Gualandi, C. Maes, A. Lodi and B. Stellato, "Batched first-order methods for parallel LP solving in MIP", arXiv 2601.21990 (2026); Y. Guan, S. Luo, P. Li, T. Chen and K. Xu, "B³-PWL: GPU-batched branch-and-bound for piecewise-linear optimization with SOS2 constraints", arXiv 2608.28988 (2026); J. Liu and A. Lodi, "From sequential nodes to GPU batches: parallel branch and bound for optimal k-sparse GLMs", arXiv 2605.22188 (2026); R. Lucas, X. Meng and R. Mazumder, "A GPU-accelerated nonlinear branch-and-bound framework for sparse linear models", arXiv 2602.04551 (2026); 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), whose per-evaluation and over-EAGO figures are the authors' own, from the package documentation. The embedded-control literature is the one place where branch and bound for small mixed-integer quadratic programs has been studied as a real-time problem. Its CPU work is the context for the two GPU papers in the table. That work covers ADMM-based QP solvers inside branch and bound for embedded hybrid model predictive control, structure-exploiting trees warm-started across control steps, and certification of the worst-case node count of a branch and bound for MIQP.A. Bemporad and V. V. Naik, "A numerically robust mixed-integer quadratic programming solver for embedded hybrid model predictive control", IFAC-PapersOnLine 51 (2018); B. Stellato, V. V. Naik, A. Bemporad, P. Goulart and S. Boyd, "Embedded mixed-integer quadratic optimization using the OSQP solver", ECC 2018; P. Hespanhol, R. Quirynen and S. Di Cairano, "A structure exploiting branch-and-bound algorithm for mixed-integer model predictive control", ECC 2019; S. Shoja, D. Arnström and D. Axehill, "Overall complexity certification of a standard branch and bound method for mixed-integer quadratic programming", ACC 2022. No general GPU branch and bound for MIQP was found.
The frontier batch and its inflation
Both designs trade pruning for parallel slack, and the quantity they trade is the same. We define it, prove what it consists of, state the host loop, and measure it.
Definition 6.4.1 (frontier batch; batch inflation). A frontier batch is a set of \(b\) open nodes whose bounds are evaluated together, with pruning performed only against the incumbent known when the batch was formed. Let \(N(b)\) be the number of bound evaluations the search performs with batch size \(b\). The batch inflation is \(N(b)/N(1)\).
Proposition 6.4.2 (what inflation consists of). Consider best-first frontier batching on a tree with distinct monotone bounds. Then \(N(b) \ge N(1)\) for every \(b\). The extra evaluations \(N(b) - N(1)\) are of two kinds. The first kind is a node of a batch that an incumbent found earlier in the same round would have pruned. The second kind is a non-essential node bounded in a round that ended with an incumbent still above its bound. The batch reached it only because the essential frontier held fewer than \(b\) nodes. If the incumbent \(z_{\mathrm{inc}} = z^\star\) is known from the start, every batch size evaluates exactly the essential nodes of Definition 6.2.1, the same set as the sequential search.
Proof. With distinct bounds the sequential best-first run evaluates exactly the essential nodes and the optimal leaf, as in the proof of Theorem 6.2.5, so \(N(1) = E + 1\). Every correct run generates and evaluates every essential node, and it must evaluate the optimal leaf to find the optimum, so \(N(b) \ge E + 1\). Every extra node is therefore non-essential. Such a node enters a batch only if fewer than \(b\) open nodes have a smaller bound. Every open essential node has a smaller bound, so the essential frontier held fewer than \(b\) nodes at that moment. If the round ends with an incumbent at most the node's bound, an incumbent found earlier in the same round would have pruned it, which is the first kind. Otherwise it is the second. With \(z_{\mathrm{inc}} = z^\star\) fixed, the pruning test depends on the node alone, so the surviving set is the set of essential nodes whatever the order or the batching. ∎
(When the two kinds of waste arise) A batch assumes that the incumbent will not change during the round and that the frontier is wide enough to fill it with essential work. Early in a search the first assumption almost always holds, since the incumbent moves rarely, and the second holds once the frontier is wide. Late in the search both fail, since the frontier is thin and consists of nodes near \(z^\star\). Both kinds of extra evaluation have the same cause: the sequential run reaches the optimum first, and the batched run has already bounded nodes it would have pruned. That is why the GPU designs lean so hard on heuristics that supply a strong incumbent early. Gmys's design has CPU heuristic threads for this, and cuOpt and the 2026 competition entries of Section 7.7 have GPU heuristics. It is also why the batch size must follow the frontier width during the run. The batch figure below measures the first kind as its red share and the total of both kinds as the node count relative to the sequential run. Monotonicity of \(N(b)\) in \(b\) is not claimed, and it does not hold. On the knapsack of the script below \(N(6) = 132{,}385\) exceeds \(N(7) = 132{,}384\), because a change of batch size realigns which nodes share a round. The Lai–Sahni anomalies of Theorem 6.2.4 are the same phenomenon stated for processors instead of batches.
Algorithm 6.4.3 (frontier-batched branch and bound with a device bound). The pseudocode below is the pool-based design as a host loop around one batched kernel.
Algorithm 6.4.3 Frontier-batched branch and bound with a device bound
(the pool-based design)
Input root node; a bound kernel that evaluates up to M nodes at once;
a branching rule; thresholds m <= M; an incumbent z_inc from a
cheap heuristic, or +inf; tolerance eps
Output optimal value and proof
State host open list L (a heap for best-first order); device buffers
for a structure of arrays of boxes
1. L <- {root}
2. while L is non-empty:
3. if |L| < m: (a thin frontier cannot fill a device)
pop one node, bound it on the host, branch, push the children
4. else:
remove a batch of q = min(|L|, M) nodes
(the q best bounds; for depth-first the q newest);
copy their boxes l, u, and any warm-start data, to the device
as a structure of arrays
5. device:
branch the batch into 2q children;
bound all of them in one kernel (a batched first-order LP
with safe bounds, an interval or McCormick evaluation of
the DAG, or a Lagrangian bound);
optionally round or locally search the children for
candidate incumbents;
return the bounds and the best candidate
6. host:
z_inc <- min(z_inc, candidate);
push the children with bound < z_inc - eps;
when z_inc has improved by more than a threshold, sweep L
and drop every node whose bound is >= z_inc - eps
7. return z_inc and its solution
Invariant
the device never decides pruning alone, so a stale incumbent on the
device costs work but never correctness; every node in L has bound
below z_inc - eps after each sweep.
Cost per round
one launch latency; q bounds spread over the device's lanes;
O(q log |L|) host heap operations; O(q) bytes each way. The number
of rounds is at least the depth of the essential subtree.
Parallel
step 5 across the batch; the host bookkeeping of steps 4 and 6 stays
serial unless the node store lives on the device, as in the
device-resident design. At q of the order of 10^4 the host heap is
the bottleneck.
One round of Algorithm 6.4.3: host and device
HOST DEVICE
2. L non-empty? --no--> 7. return
| yes z_inc and
| its solution
v
3. |L| < m ? --yes--> pop one node,
| no bound it on the
| host, branch,
| push the children;
| back to 2
v
4. remove a batch of q = min(|L|, M)
nodes; copy their boxes l, u and
any warm-start data, as a structure
of arrays ------------------------> 5. branch the batch into
2q children; bound all
of them in one kernel;
optionally round or
locally search them
|
the bounds and |
the best |
candidate |
6. z_inc <- min(z_inc, candidate) <------------------+
push the children with
bound < z_inc - eps;
sweep L when z_inc has improved
by more than a threshold
|
+--> back to 2
The device buffers of Algorithm 6.4.3: a structure of arrays
q boxes [l, u] in n variables, one contiguous array per field
node 1 node 2 node 3 . . . node q
l_1 [ . | . | . | . . . | . ]
l_2 [ . | . | . | . . . | . ]
:
l_n [ . | . | . | . . . | . ] all lower
u_1 [ . | . | . | . . . | . ] bounds,
: then all
u_n [ . | . | . | . . . | . ] upper bounds
^ ^ ^
lane 1 lane 2 lane 3 . . .
the lanes of a warp, one node each, read consecutive addresses
(The two thresholds, and the script) The device buffers of the algorithm hold the boxes as a structure of arrays: one contiguous array per field, all lower bounds and then all upper bounds, so that a warp reads consecutive addresses. The algorithm has two parameters, \(m\) and \(M\), and the trade-off they control has a sharp knee. The script runs the loop on the knapsack of Section 6.2 with eight batch sizes from one to \(65{,}536\), counts rounds and bound evaluations, and prices a round under an illustrative device model. The model charges a launch latency \(L = 10\,\mu\mathrm{s}\) per round and \(2\,\mu\mathrm{s}\) per node spread over \(\min(b, 4096)\) lanes, against \(0.5\,\mu\mathrm{s}\) per node on a CPU core. The constants are illustrative and the instance maximizes. In the terms of Definition 6.4.1 the script's sequential run gives \(N(1) = 132{,}376\) bound evaluations on this instance, batches of \(1{,}024\) give \(N(1{,}024) = 141{,}311\), an inflation of \(1.07\), and batches of \(65{,}536\) give \(N(65{,}536) = 1{,}048{,}575\), an inflation of \(7.92\); the table it prints is these ratios for all eight sizes, with the rounds and the modelled times beside them.
# Frontier-batch branch and bound on the knapsack of the shared-pool
# script (Section 6.2).
#
# Each round takes up to b nodes from the best-first open list and bounds
# them together, as one device kernel would, pruning only against the
# incumbent known when the round began; the host then prunes the children
# against the updated incumbent. Counts rounds and bound evaluations, and
# applies an illustrative time model.
import heapq
import random
random.seed(5)
n = 30
w = [random.randint(20, 80) for _ in range(n)]
v = [wi + 10 for wi in w]
cap = sum(w) // 2
# the items in order of value per unit of weight, best first
order = sorted(range(n), key=lambda i: -v[i] / w[i])
w = [w[i] for i in order]
v = [v[i] for i in order]
def bound(k, wt, val):
"""Dantzig's bound of a node that has decided items 0 .. k-1.
wt and val are the weight and the value of the items taken so far.
"""
for i in range(k, n):
if wt + w[i] <= cap:
wt += w[i]
val += v[i]
else:
return val + (cap - wt) * v[i] / w[i]
return val
def run(b):
"""Best-first search in rounds of up to b nodes.
Returns (rounds, bound evaluations, optimum).
"""
inc = 0
rounds = 0
evals = 0
pool = [(-bound(0, 0, 0), 0, 0, 0)]
while pool:
rounds += 1
batch = [heapq.heappop(pool) for _ in range(min(b, len(pool)))]
# the device knows only this value during the round
inc0 = inc
children = []
for key, k, wt, val in batch:
if -key <= inc0:
continue
# one bound evaluation per surviving node
evals += 1
if k == n:
inc = max(inc, val)
continue
for take in (1, 0):
if take and wt + w[k] > cap:
continue
nw, nv = wt + w[k] * take, val + v[k] * take
bb = bound(k + 1, nw, nv)
if bb > inc0:
children.append((bb, k + 1, nw, nv))
# host: prune with the incumbent after the round
for bb, k, nw, nv in children:
if bb > inc:
heapq.heappush(pool, (-bb, k, nw, nv))
return rounds, evals, inc
R1, E1, opt = run(1)
L = 10e-6 # launch latency
LANES = 4096 # lanes
t_gpu = 2e-6 # per-node time on a lane
t_cpu = 0.5e-6 # per-node time on a CPU core
print(f"knapsack n={n}, optimum {opt}")
print(f"sequential best-first: {E1} bound evaluations, "
f"CPU time model {E1 * t_cpu * 1e3:.1f} ms")
print()
print(f"{'batch b':>8} {'rounds':>7} {'evals':>8} {'N(b)/N(1)':>9} "
f"{'GPU time model':>15}")
for b in (1, 16, 64, 256, 1024, 4096, 16384, 65536):
R, E, _ = run(b)
# per round: one launch, then the nodes spread over min(b, LANES) lanes
T = R * L + E * t_gpu / min(b, LANES)
print(f"{b:>8} {R:>7} {E:>8} {E / E1:>9.2f} {T * 1e3:>13.2f} ms")
knapsack n=30, optimum 991
sequential best-first: 132376 bound evaluations, CPU time model 66.2 ms
batch b rounds evals N(b)/N(1) GPU time model
1 165694 132376 1.00 1921.69 ms
16 10365 132431 1.00 120.20 ms
64 2604 132735 1.00 30.19 ms
256 669 134399 1.02 7.74 ms
1024 187 141311 1.07 2.15 ms
4096 69 168329 1.27 0.77 ms
16384 45 344063 2.60 0.62 ms
65536 40 1048575 7.92 0.91 ms
(What the inflation script shows) The inflation is negligible up to \(b = 256\), \(7\,\%\) at \(b = 1{,}024\), \(27\,\%\) at \(b = 4{,}096\), and then it explodes. At \(b = 65{,}536\) the search degenerates toward enumeration, \(2^{20} - 1\) evaluations, because whole levels of the tree are bounded before any leaf can prune them. The time model has its minimum at \(b \approx 16{,}384\). There an inflation of \(2.6\) is paid for by \(24\) fewer launches than at \(b = 4{,}096\): \(0.24\) ms of latency saved against \(0.09\) ms of extra bounding spread over \(4{,}096\) lanes. The minimum moves left as the per-node cost grows (an LP instead of a greedy bound) and right as the launch overhead grows (a multi-GPU reduction per round). It also moves during the run as the frontier narrows, so a fixed batch size is wrong at both ends of the search. The script is sequential and costs one bound evaluation per surviving node. In a real device loop step 5 of the algorithm is the kernel, and everything in the script's inner loop would run across lanes.
The batch figure
The figure runs Algorithm 6.4.3 on the fifteen-item knapsack of the parallel figure. Its device charges a latency of \(L\) time units per batch plus \(b/T\) for its \(b\) slots, whether or not they are filled. Its sequential reference is the same procedure one node at a time, which bounds \(893\) nodes in best-bound order. That count is not the parallel figure's one-worker run of \(1{,}057\) nodes and \(1{,}360.8\) time units, which prices node costs and pops differently, and the two figures are not to be compared on that number. In the default view, with \(b = 32\), \(L = 4\) and \(T = 64\), the search takes \(33\) batches and \(148.5\) time units against \(893\) sequential, a speedup of \(6.01\). It bounds \(897\) nodes, \(4\) of them wasted, at an occupancy of \(85\,\%\), occupancy being the share of the batch slots the device was charged for that actually held a node, and \(132.0\) of the \(148.5\) time units are latency. The first five batches hold \(1, 2, 4, 8\) and \(16\) nodes, which is ramp-up. The best batch size for that device is \(128\), at \(96.0\) time units and a speedup of \(9.30\), with \(116\) of its \(1{,}145\) nodes wasted. At \(256\) the time rises again to \(112.0\), with \(1{,}359\) nodes bounded and \(128\) wasted, and an occupancy of \(38\,\%\) because only \(3\) of \(14\) batches are full. The whole curve at the default device is U-shaped. It starts at \(3{,}586.0\) time units at \(b = 1\), where the device is four times slower than the sequential run because every node pays the latency alone. It then falls through \(1{,}802.0\), \(914.1\), \(470.3\), \(250.8\) and \(148.5\) with nothing wasted up to \(b = 16\), reaches \(105.0\) at \(b = 64\) with \(58\) wasted nodes and \(96.0\) at \(128\), and rises to \(112.0\) at \(256\). The cases say the rest. With no latency the best batch size is \(1\), since a single node already runs the device at full throughput. On a slow device with \(T = 8\) the best size is \(64\). In depth-first order with \(b = 64\) the run bounds \(2{,}617\) nodes against a depth-first sequential count of \(1{,}833\), a speedup of \(7.64\). The lower chart draws the two quantities Proposition 6.4.2 separates. In ink it draws the nodes bounded relative to the sequential run, which counts both kinds of extra evaluation. In red it draws the share bounded for nothing because an incumbent found earlier in the same batch would have discarded them, which is the first kind alone.
Two more propositions bound what a batched design can achieve. The first is Amdahl's law applied to the node.
Proposition 6.4.4 (Amdahl bound for GPU-accelerated branch and bound). If a fraction \(s \in (0, 1)\) of the sequential per-node work (node selection, branching, bookkeeping, host–device transfer) stays sequential and the remaining fraction \(1 - s\) (bounding) is accelerated by a factor \(g\), then the speedup of the whole search is at most \(1/(s + (1 - s)/g) < 1/s\), before counting the inflation of Proposition 6.4.2.G. M. Amdahl, "Validity of the single processor approach to achieving large scale computing capabilities", AFIPS Conference Proceedings 30 (1967), 483–485.
Proof. Divide the sequential time \(T\) into \(sT\) and \((1 - s)T\), and the second part becomes \((1 - s)T/g\). ∎
(Why the node data must live on the device) With \(s = 0.1\) and \(g = 100\) the bound is \(9.2\), and with \(s = 0.01\) it is \(50\). This is why every successful GPU branch and bound keeps the node data structure on the device rather than shipping nodes back and forth. Gmys's interval encoding, the padded node tensors of Liu and Lodi, and the block-tiled sparse kernels of B³-PWL all do this. Perumalla and Alam draw the same conclusion. The best prospects for GPU mixed-integer programming lie with problems whose relaxation matrix "fits entirely within one accelerator's memory" and whose trees "cannot be fully contained within a small number of computational nodes".K. S. Perumalla and M. Alam, "Design considerations for GPU-based mixed integer programming on parallel computing platforms", ICPP Workshops 2021.
The second returns to determinism, now for a batch.
Proposition 6.4.5 (a fixed batch schedule makes frontier batching deterministic). In Algorithm 6.4.3 suppose five things. Each round's batch consists of the first \(b\) nodes of \(L\) in a canonical order (bound, then identifier). Identifiers are assigned canonically (parent, child side). Bounds are computed by fixed schedules in one precision. The incumbents found in a round are merged at the end of that round in a fixed order. And the device never prunes on its own. Then the state after round \(r\) is a function of the input and \(b\) alone, and the run is deterministic for every number of lanes, blocks and devices.
Proof. Induction on rounds, exactly as for Algorithm 6.2.15: the five conditions fix the batch, its identifiers, every bound, the incumbent and the surviving children. ∎
(The price of a deterministic batch) The price has three parts. The first is a barrier per round, since the next batch cannot be formed until the merge is complete. The second is the one-round incumbent lag, which is the inflation already measured. The third is idle lanes, because a round in lockstep lasts as long as its slowest node. The lag exists in an opportunistic round as well. What determinism removes is the option of applying an incumbent in the middle of a round, and what it adds is the barrier. The determinism script of Section 6.2 showed the gain: the fixed schedule at \(b = 1{,}024\) gives one node count in every run, and the randomized one gives twelve counts in twelve runs. The idle lanes can be modelled in one line.
Two rounds under a fixed batch schedule (Proposition 6.4.5)
round r round r + 1
lane 1 #########....... | ######...... |
lane 2 #####*.......... | ############ |
lane 3 ################ | ########.... |
. . . | |
lane b ######.......... | #####....... |
^ ^
barrier barrier
# a node's work, in lockstep . an idle lane
* an incumbent found in round r: merged at the barrier in a
fixed order, it prunes from round r + 1 on
the price: one barrier per round, the one-round incumbent lag
and the idle lanes; the lane efficiency sum t_i / (b max t_i)
of Proposition 6.4.6 is the share of # in a round
Proposition 6.4.6 (lane efficiency of a lockstep round). If a round evaluates \(b\) nodes with work \(t_1, \dots, t_b\) in lockstep and ends when the slowest finishes, its lane efficiency is \(\sum_i t_i / (b \max_i t_i)\). For independent, identically distributed work with finite mean and an unbounded distribution, the expected efficiency tends to zero as \(b \to \infty\).
Proof. The efficiency lies in \([0, 1]\). By the strong law of large numbers \(\frac{1}{b}\sum_i t_i \to \mathbb E[t_1]\) almost surely, and \(\max_i t_i \to \infty\) almost surely because the distribution is unbounded. Their ratio therefore tends to zero almost surely, and dominated convergence carries the limit to the expectation. ∎
The script draws lognormal node work with median one and a chosen spread and forms 2,000 rounds of \(b\) nodes each. It prints the mean lane efficiency of Proposition 6.4.6 for three batch sizes and four spreads, together with the factor between the median node and the 84th percentile.
# Lane efficiency of a lockstep round (Proposition 6.4.6).
#
# b nodes are bounded together, and the round ends when the slowest
# finishes. Per-node work is lognormal with median 1 and
# log-standard-deviation sigma; efficiency = sum t / (b max t), averaged
# over 2,000 rounds for each batch size and each sigma.
import numpy as np
rng = np.random.default_rng(7)
rounds = 2000
print(f"{'sigma':>6} {'b=32':>8} {'b=256':>8} {'b=4096':>8} "
f"84th percentile of node work")
print(" " * 36 + "over its median")
for s in (0.1, 0.25, 0.5, 1.0):
row = []
for b in (32, 256, 4096):
t = rng.lognormal(0.0, s, size=(rounds, b))
row.append(np.mean(t.sum(axis=1) / (b * t.max(axis=1))))
effs = " ".join(f"{e:>8.3f}" for e in row)
print(f"{s:>6.2f} {effs} a factor {np.exp(s):.2f}")
sigma b=32 b=256 b=4096 84th percentile of node work
over its median
0.10 0.819 0.758 0.700 a factor 1.11
0.25 0.616 0.511 0.418 a factor 1.28
0.50 0.410 0.280 0.187 a factor 1.65
1.00 0.212 0.103 0.046 a factor 2.72
(What the lane-efficiency script shows) With a spread of node work of a factor \(1.65\) between the median and the 84th percentile, a round of \(4{,}096\) lanes keeps them busy less than a fifth of the time. Even at \(b = 256\) more than half the lanes idle. This is a model, not a measurement of any solver, and it says what the remedies must be. One remedy is batches of nodes with similar predicted work, which is a scheduling problem that a learned subtree-size estimator could serve. The other is a kernel with a fixed iteration budget per round, as in the batched first-order LP of Section 7.4. There every column does the same work by construction, and the cost moves to columns that converged early and iterate until the round ends or are compacted away at a round boundary. The script costs one random draw per lane and one maximum per round. The rounds are independent, so everything in it parallelizes across rounds, and the maximum within a round is itself a reduction.
Regular and irregular work in a spatial node
The spatial branch-and-bound node of Section 3.5 is a sequence of steps, and the question for a device is which of them have constant shape across a batch of boxes. Two abbreviations in the table are standard: SpMV is a sparse-matrix–vector product, and SpMM a sparse-matrix–dense-matrix product, which is \(K\) such products sharing one matrix. A Jacobi sweep, named in the FBBT row, propagates every constraint from the bounds of the previous sweep at once, instead of each constraint from the bounds its predecessors have just tightened, so that the constraints of one sweep are independent and can run as a map; it may need more sweeps to reach the same fixed point.
| step | shape across a batch of boxes | regular? | GPU status (October 2026) |
|---|---|---|---|
| node selection | a heap or a canonical order over the frontier | no | host; a device-side priority structure is open |
| FBBT | one pass per constraint, the same DAG per box, data-dependent number of sweeps | yes | a map over constraints and boxes (Jacobi); not yet built for MINLP; Section 7.8 |
| OBBT | \(2n\) LPs per box sharing one matrix | yes | batched LPs are Blin et al.'s setting; not built |
| McCormick relaxation (the LP's rows) | the same rows per box with box-dependent coefficients | yes, if padded | coefficients in one kernel; a shared-matrix SpMM becomes a batched SpMV: Section 7.8 |
| the LP bound | simplex: pivots, warm start from the parent | no | one core |
| PDHG: matrix-vector products, a fixed iteration budget per round | yes | batched PDHG (Blin et al.; cuOpt's batch PDLP strong branching); safe bound: Section 7.3 | |
| the interval or McCormick bound by evaluation | the same straight-line program for every box and every point | yes | MAiNGO's GPU interval bounder; ParBB's kernels |
| local NLP for an incumbent | an interior-point or SQP solve per box with its own factorization and iteration count | no | condensed-KKT GPU NLP exists for one problem at a time (Section 7.5); batching is open |
| branching | one split per box | yes | trivial |
| pruning and the incumbent | a comparison and a min-reduction | yes | exact with atomic minimum instructions, which are exact on integers; a dual objective needs a fixed schedule (Proposition 6.2.18) |
(The two bounding engines that run on a device) The regular steps are the ones with a GPU prototype, and the two bounding engines that have actually run on a device for factorable problems are both evaluations rather than LPs. The MAiNGO interval bounder partitions a node's box into \(m^d\) sub-boxes and evaluates the mean value form on all of them in one kernel with directed rounding. The form is an enclosure whose excess width is of second order in the width of the box, and its statement is short enough to give here.
Proposition 6.4.7 (the mean value form; Moore, 1966; Neumaier, 1990). Let \(f\) be continuously differentiable on a box \(X \subseteq \mathbb{R}^n\) with centre \(c\), and let \(\nabla F(X)\) be an inclusion-isotone interval extension of the gradient \(\nabla f\), in the sense of Definition 2.6.3, evaluated on \(X\). The mean value form of \(f\) on \(X\) is the interval
\[f_{\mathrm{MV}}(X) \;=\; f(c) + \nabla F(X)^\top (X - c) \;=\; f(c) + \sum_{i=1}^n \nabla F_i(X)\,(X_i - c_i).\](a) It encloses the range: \(f(X) \subseteq f_{\mathrm{MV}}(X)\). (b) If \(\nabla f\) is Lipschitz on \(X\) and each component \(\nabla F_i\) satisfies the first-order bound \(w(\nabla F_i(X)) \le 2L\,w(X)\) of Proposition 2.6.5(b), there is a constant \(C\) with \(w(f_{\mathrm{MV}}(X)) - w(f(X)) \le C\,w(X)^2\). The excess width is of second order in the width of the box, against first order for the natural interval extension.R. E. Moore, Interval Analysis (Prentice-Hall, 1966); A. Neumaier, Interval Methods for Systems of Equations (Cambridge University Press, 1990), Chapter 2; R. E. Moore, R. B. Kearfott and M. J. Cloud, Introduction to Interval Analysis (SIAM, 2009), Chapter 6. The theorem numbers in Neumaier and in Moore, Kearfott and Cloud were not verified for this post, so both are cited by chapter.
Proof sketch. (a) For \(x \in X\) the mean value theorem gives \(f(x) = f(c) + \nabla f(\xi)^\top (x - c)\) for some \(\xi\) on the segment from \(c\) to \(x\), and the segment lies in \(X\) because a box is convex. By Theorem 2.6.4, \(\nabla f(\xi) \in \nabla F(X)\), and \(x - c \in X - c\), so \(f(x) \in f_{\mathrm{MV}}(X)\). (b) Write \(g = \nabla f(c)\). Sub-distributivity of interval arithmetic gives \(f_{\mathrm{MV}}(X) \subseteq f(c) + g^\top (X - c) + (\nabla F(X) - g)^\top (X - c)\). The interval \(f(c) + g^\top (X - c)\) has width \(\sum_i |g_i|\,w(X_i)\). The third term has width at most \(\sum_i w(\nabla F_i(X))\,w(X_i) \le 2nL\,w(X)^2\), because \(g_i \in \nabla F_i(X)\) and \(X_i - c_i\) is symmetric about zero. Taylor's theorem at the two vertices with coordinates \(c_i \pm \tfrac12 w(X_i) \operatorname{sign}(g_i)\), with the Lipschitz constant of \(\nabla f\) bounding the remainder, shows that the range \(f(X)\) has width at least \(\sum_i |g_i|\,w(X_i) - C'\,w(X)^2\). Subtracting the two widths gives the bound with \(C = 2nL + C'\). ∎
(What the form buys, and its cost) The picture is that the form linearizes \(f\) at the centre and pays for the curvature only through the width of the gradient enclosure, which itself shrinks with the box. The minimum over the sub-boxes is the bound, and its wall-clock speedup over CPU interval arithmetic without partitioning is reported as three orders of magnitude. On some case studies it is competitive with MAiNGO's default McCormick relaxations. Its cost is the exponent \(d\), which confines it to nodes with few nonconvex variables. The generated McCormick kernels of SourceCodeMcCormick.jl evaluate pointwise relaxations, intervals and subgradients for thousands of boxes or points at once. The authors' figures are about 9 nanoseconds per evaluation on the GPU against 237 on the CPU, and 11 to 22 times over EAGO for the branch and bound built on them. Its operation library was small when this post read it.Zhang et al. (2025) and Gottlieb, Xu and Stuber (2026), cited above; Section 7.8 states the first experiments for both. Both are Type 1 parallelism in the sense of Definition 6.1.2, inside the bound, and both leave the tree on the host.
MAiNGO's GPU interval bounder on the box X of one node
X cut into m pieces along each of d coordinates: m^d sub-boxes
X_k with centres c_k (two of the d coordinates drawn)
1 2 . . . m
+--------+--------+--------+--------+
1 | | | | |
+--------+--------+--------+--------+
2 | | | | |
+--------+--------+--------+--------+
: | | | | |
+--------+--------+--------+--------+
m | | | | |
+--------+--------+--------+--------+
|
| one kernel, directed rounding
v
on every X_k: f_MV(X_k) = f(c_k) + grad F(X_k)^T (X_k - c_k)
|
v
the bound: the minimum over the sub-boxes, that is, over k
of the lower end of f_MV(X_k)
the cost grows with the exponent d; the tree stays on the host
What has been measured, and what has not
(The record, and its gaps) Every GPU branch-and-bound result verified for this post is for flow-shop, N-Queens, knapsack, sparse regression, piecewise-linear MIPs, neural-network verification, small mixed-integer QPs from control or factorable NLPs. Only the Chapel codes (flow-shop and N-Queens) and Gmys's flow-shop solver report multi-GPU scaling, and both initialize with the optimum. The knapsack results of Lalami and El Baz and of Boukedjar, Lalami and El Baz are on one GPU, as are the sparse-regression, piecewise-linear and verification results. For MINLPLib instances the published parallel results are CPU-only: FiberSCIP's MINLP experiments, ParaSCIP's MILP records, Xpress Global's thread-scaling table on its internal set, and MAiNGO's MPI mode without a reported curve. As of 5 October 2026 no published study measures the strong scaling of a global MINLP solver across more than one node on MINLPLib or QPLIB instances. No benchmark measures any GPU code on MINLPLib. Mittelmann's three GPU tracks are LP feasibility, convex continuous QPLIB and sparse SDP, and his MINLP page of 26 February 2026 has no GPU column and no thread-scaling arm.H. D. Mittelmann, "Benchmarks for Optimization Software", plato.asu.edu/bench.html, entries dated through 1 October 2026, and "Mixed Integer Nonlinear Programming Benchmark (MINLPLIB)", plato.asu.edu/ftp/minlp.html, page dated 26 February 2026, both read 5 October 2026. The statement is one of absence after the literature sweeps recorded for this post, and one paper falsifies it, which is why it carries a date. Section 5.2 gives the benchmark methodology. The protocol that would fill both gaps is the scaling and determinism measurement of Algorithm 7.8.2, and the dual integral it reports is defined there as Definition 7.8.1. Neither has numbers yet, because the experiment has not been run by anyone, and Section 7.8 ranks it among the open problems of the GPU programme.
Where this is used
Every commercial MILP solver runs a parallel tree on one machine by default. Gurobi and Xpress run it deterministically by default, and CPLEX's default automatic mode is documented to yield deterministic results. Gurobi's distributed MIP uses racing ramp-up, as UG does, and CPLEX's distributed parallel MIP is governed by the same parallel-mode switch. The records on open instances belong to the UG framework over SCIP and Xpress, and the GPU records to Gmys's flow-shop solver and the Chapel codes, on combinatorial trees whose bounds cost nanoseconds. Of the global MINLP solvers, Xpress Global and FiberSCIP are the ones with published parallel results, and MAiNGO has an MPI mode. cuOpt, the one production MIP solver with GPU heuristics, keeps its tree on the CPU and offers determinism only as an experimental mode that excludes those heuristics. The theory of this section is what those systems implement. Ties are broken by a total order (Theorem 6.2.5). Idle workers steal the oldest nodes at random (Theorem 6.2.9). Ramp-up strategies exist because of Proposition 6.2.12, and termination is detected by counting or by a token (Proposition 6.2.13). Logical clocks give determinism (Algorithm 6.2.15), and batch sizes follow the frontier (Proposition 6.4.2).
What parallelizes
The tree is irregular and the node is regular. The two designs that have worked assign them accordingly: the tree to the host, or to the device in a constant-shape encoding, and the regular work to the device in batches. For a spatial branch and bound the regular work is the McCormick LP with its box-dependent coefficients, FBBT as a map over constraints, OBBT as \(2n\) LPs with one matrix, interval and McCormick evaluation over sub-boxes, and multistart for incumbents. The irregular work is node selection, the simplex warm start, and the local NLP. The quantities to design against are the four of this section: the slackness \(T_1/(P\,T_\infty)\), the ramp-up bound, the batch inflation \(N(b)/N(1)\), and the lane efficiency of a lockstep round. Correctness under a nondeterministic order comes from directed rounding and safe bounds (Proposition 6.2.21 and Section 7.3). Reproducibility, if a product needs it, comes from fixed schedules and a fixed batch schedule (Propositions 6.2.18 and 6.4.5). Its price for the CPU is two to nine percent of time or a third of the threads, and for a GPU it has not been measured. Section 7 takes the device itself as its subject, and Section 7.8 turns the open items of this section into experiments.