GPU-Parallel Interval Lower Bound Solver

• The paper presents a GPU-parallel lower bounding solver that partitions branch-and-bound nodes into subdomains and uses rigorous interval arithmetic to compute tight lower bounds.
• It leverages the Mean Value Form with adaptive partitioning strategies to improve bound accuracy while reducing iterations in nonconvex global optimization.
• Empirical results demonstrate orders-of-magnitude speedups over CPU implementations and competitive performance compared to traditional McCormick relaxations.

An interval-based GPU-parallel lower bounding solver is a lower-bounding module for branch-and-bound methods that uses interval arithmetic, temporary domain partitioning, and massive GPU parallelism to compute rigorous lower bounds over many subdomains simultaneously. In the formulation developed for deterministic global optimization, each branch-and-bound node with domain box $X \subseteq \mathbb{R}^n$ is partitioned into numerous subdomains, interval enclosures of the objective and constraints are evaluated on the GPU, and infeasible subdomains are discarded by interval tests; the surviving subdomains define a conservative node lower bound. The design described in "Accelerating Deterministic Global Optimization via GPU-parallel Interval Arithmetic" integrates this mechanism into the open-source solver MAiNGO through CUDA, using either a single monolithic kernel or a CUDA Graph, and reports speedups of three orders of magnitude against CPU interval arithmetic without domain partitioning (Zhang et al., 28 Jul 2025).

1. Position within spatial branch-and-bound

Spatial Branch-and-Bound for nonconvex NLPs hinges on bounding: tighter lower bounds prune nodes earlier and reduce iterations, but computing tight bounds is expensive if done serially on CPUs. The interval-based GPU-parallel lower bounding solver addresses this bottleneck by moving lower-bound construction from a single-box CPU computation to a subdomain-parallel GPU computation. The target problem is a nonconvex global optimization problem with objective $f:\mathbb{R}^n \to \mathbb{R}$ over a decision box $X \subseteq \mathbb{R}^n$, together with inequality constraints $g_j(x)\le 0$ for $j=1,\dots,n_g$ and equality constraints $h_j(x)=0$ for $j=1,\dots,n_h$ (Zhang et al., 28 Jul 2025).

At a branch-and-bound node, the node domain is a box $X=(X_1,\dots,X_n)$ with $X_i=[\underline{X}_i,\overline{X}_i]$. The solver temporarily partitions this box into $S$ subdomains $f:\mathbb{R}^n \to \mathbb{R}$0, evaluates interval bounds for the objective and constraints on each $f:\mathbb{R}^n \to \mathbb{R}$1, removes certainly infeasible subdomains through interval feasibility tests, and aggregates the remaining objective enclosures into a node lower bound. The central motivation is not only parallel throughput. The design also decouples GPU workload from the number of open nodes in the branch-and-bound tree, thereby avoiding the startup bottleneck of GPU-parallel node-level bounding; it reduces warp divergence relative to McCormick relaxations; and it exploits the Mean Value Form’s higher convergence order with respect to domain width (Zhang et al., 28 Jul 2025).

A common misconception is that GPU acceleration of branch-and-bound must operate across many tree nodes at once. The subdomain lower bounding design instead parallelizes inside a single node by evaluating many subdomains in parallel. In the reported framework, this internal parallelism is the mechanism that preserves rigor while producing tighter bounds than regular interval arithmetic without partitioning (Zhang et al., 28 Jul 2025).

2. Interval arithmetic, feasibility tests, and lower-bound construction

The interval component is rigorous because directed rounding ensures outward rounding and thus inclusion. For basic operations, interval addition and multiplication are performed with $f:\mathbb{R}^n \to \mathbb{R}$2 and $f:\mathbb{R}^n \to \mathbb{R}$3 so that lower endpoints are rounded downward and upper endpoints upward. Division is implemented via multiplication by an outward-rounded reciprocal, handling zero-crossings robustly. For transcendental functions such as $f:\mathbb{R}^n \to \mathbb{R}$4, $f:\mathbb{R}^n \to \mathbb{R}$5, and $f:\mathbb{R}^n \to \mathbb{R}$6, the implementation uses double-direction rounding and careful enclosures to guarantee inclusion despite 1–3 ulp inaccuracies of CUDA math functions (Zhang et al., 28 Jul 2025).

The local lower-bound computation is organized around subdomains. For each subdomain $f:\mathbb{R}^n \to \mathbb{R}$7, the interval objective enclosure defines

$f:\mathbb{R}^n \to \mathbb{R}$8

Feasibility is determined by interval tests on the constraint enclosures. For an inequality constraint $f:\mathbb{R}^n \to \mathbb{R}$9, if $X \subseteq \mathbb{R}^n$0, the subdomain is certainly infeasible and is discarded; if $X \subseteq \mathbb{R}^n$1, the subdomain is certainly feasible for that constraint; otherwise it remains potentially feasible. For an equality constraint $X \subseteq \mathbb{R}^n$2, the subdomain is retained only if $X \subseteq \mathbb{R}^n$3 (Zhang et al., 28 Jul 2025).

The node lower bound is the minimum over feasible subdomains:

$X \subseteq \mathbb{R}^n$4

The framework also allows a refinement via the global interval hull

$X \subseteq \mathbb{R}^n$5

where $X \subseteq \mathbb{R}^n$6 is an inclusion function such as the Mean Value Form or Natural Interval Extension. In that case, the node lower bound is $X \subseteq \mathbb{R}^n$7, while feasibility is integrated through the analogous hulls of the constraint enclosures (Zhang et al., 28 Jul 2025).

The key interval extension is the Mean Value Form centered at the midpoint $X \subseteq \mathbb{R}^n$8. For differentiable $X \subseteq \mathbb{R}^n$9 on $g_j(x)\le 0$0,

$g_j(x)\le 0$1

and the practical enclosure is

$g_j(x)\le 0$2

Equivalently, the enclosure can be expressed coordinatewise as

$g_j(x)\le 0$3

where the interval gradient is obtained by forward-mode AD through CuTangent on each subdomain (Zhang et al., 28 Jul 2025).

The theoretical benefit of splitting is explicit. If each dimension is uniformly partitioned into $g_j(x)\le 0$4 subintervals, then $g_j(x)\le 0$5 and $g_j(x)\le 0$6. The excess width of the refined enclosure satisfies

$g_j(x)\le 0$7

for Natural Interval Extension, while for the Mean Value Form it satisfies

$g_j(x)\le 0$8

This quadratic dependence on subdomain width is the central reason that Mean Value Form benefits more strongly from subdomain splitting than Natural Interval Extension (Zhang et al., 28 Jul 2025).

3. Partitioning strategies and integration into node processing

The node-level workflow begins with a choice of subdomain count $g_j(x)\le 0$9, partitioning strategy, and interval extension. Three partitioning strategies are described.

Strategy	Definition	Intended effect
Uniform	Same number $j=1,\dots,n_g$0 per dimension, so $j=1,\dots,n_g$1	Regular refinement of all coordinates
Largest-dimension	Split only $j=1,\dots,n_g$2	Control growth of $j=1,\dots,n_g$3
Adaptive	Start with the largest uniform $j=1,\dots,n_g$4 not exceeding GPU core count, then increment $j=1,\dots,n_g$5 in dimensions with largest widths until $j=1,\dots,n_g$6 available core count	Match GPU resources while targeting wide dimensions

After partitioning, each subdomain uses $j=1,\dots,n_g$7, computes interval bounds for the objective and constraints through the Mean Value Form or Natural Interval Extension, and applies the feasibility tests. The node lower bound is then obtained by aggregating the feasible subdomains, while optional refinement computes the hulls of the objective and constraints for pruning (Zhang et al., 28 Jul 2025).

Several tightness heuristics are explicit. Increasing $j=1,\dots,n_g$8 improves bound tightness; using the midpoint as $j=1,\dots,n_g$9 is empirically robust and matches the Mean Value Form definition; adaptive partitioning reduces dependency and wrapping effects without excessive serialization; and Mean Value Form is favored over Natural Interval Extension. The Jacobian or gradient intervals are computed per subdomain, rather than reused from a parent box, because parent-box enclosures would be looser (Zhang et al., 28 Jul 2025).

The computational trade-off is equally explicit. Per node, interval Mean Value Form evaluation over all subdomains costs $h_j(x)=0$0, where $h_j(x)=0$1 depends on DAG size. Aggregation on the CPU costs $h_j(x)=0$2 for minima and $h_j(x)=0$3 for hulls if computed by a single pass, although the implementation uses sorting for refinements, which can be $h_j(x)=0$4. As $h_j(x)=0$5 grows, lower bounds tighten and branch-and-bound iterations decrease, but GPU work, data transfer, and CPU aggregation costs all increase (Zhang et al., 28 Jul 2025).

This also clarifies another frequent confusion: the temporary subdomain partition is not the branch-and-bound branching operation itself. It is a lower-bounding procedure local to one tree node, discarded after the node bound is computed. That distinction is central to why GPU work can be saturated independently of the evolving number of open nodes (Zhang et al., 28 Jul 2025).

4. CUDA realization and MAiNGO integration

The implementation is integrated into MAiNGO through CUDA, with the MAiNGO DAG supplied via MC++. Two GPU realizations are provided.

Implementation	Structure	Reported effect
Single monolithic kernel	One kernel performs partitioning and per-subdomain interval evaluation of the DAG	Simpler, but DAG branches do not execute concurrently inside the kernel
CUDA Graph	A root partition kernel is followed by kernels for DAG nodes, connected by dependencies, with dependent memory-copy nodes at the end	Higher efficiency; reduced launch overhead; potential parallel execution of independent DAG branches

In the monolithic-kernel realization, each thread ideally handles one subdomain $h_j(x)=0$6. The thread evaluates the simplified DAG using CuInterval and CuTangent with directed outward rounding. The data flow is: copy node box and DAG constants to the GPU, launch the kernel, copy per-subdomain bounds back to the CPU, then compute the hull and minimum on the CPU. The DAG evaluation order is sequential within the kernel, so independent DAG branches do not execute concurrently there (Zhang et al., 28 Jul 2025).

In the CUDA Graph realization, the graph is constructed once and then instantiated and repeatedly launched. A partition kernel creates the array of subdomains; each DAG node is replaced by a kernel and connected according to DAG dependencies; and dependent memory-copy nodes transfer bounds back only after objective and constraint kernels finish. The reported benefits are reduced kernel-launch overhead, potential parallel execution of independent DAG branches, and improved memory scheduling (Zhang et al., 28 Jul 2025).

The mapping from subdomains to GPU threads is one-to-one when possible, so the number of threads equals $h_j(x)=0$7. Memory layout is arranged for coalesced access to subdomain endpoints and per-node intermediate values. Two C++ classes manage configuration and memory across varying $h_j(x)=0$8 and interval-extension choices. Shared-memory usage is limited by the per-kernel design; main storage is in global memory; and graph orchestration is used to reduce transfers and launch overhead (Zhang et al., 28 Jul 2025).

Rigor and reproducibility depend on the interval library choices. The implementation uses FP64 arithmetic and CUDA intrinsic rounding modes for $h_j(x)=0$9, $j=1,\dots,n_h$0, $j=1,\dots,n_h$1, $j=1,\dots,n_h$2, $j=1,\dots,n_h$3, and $j=1,\dots,n_h$4. For transcendental functions such as $j=1,\dots,n_h$5, $j=1,\dots,n_h$6, $j=1,\dots,n_h$7, and $j=1,\dots,n_h$8, CuInterval uses double-direction rounding and specialized enclosures to guarantee inclusion. Given a fixed graph and device, the combination of outward rounding and dependency-respecting graph execution makes the computation deterministic (Zhang et al., 28 Jul 2025).

The main hardware bottlenecks are also explicit. Performance saturates when $j=1,\dots,n_h$9 exceeds available CUDA cores and FP64 throughput becomes the bottleneck. On the reported RTX A1000 Laptop GPU, FP64/FP32 is approximately $X=(X_1,\dots,X_n)$0. As $X=(X_1,\dots,X_n)$1 grows, serialization and CPU-side reduction overheads increase, and the benefit of further splitting may become small once nodes are already narrow (Zhang et al., 28 Jul 2025).

5. Empirical behavior and comparison with alternative bounding methods

The reported experiments were conducted on an Intel i7-13700H paired with an NVIDIA RTX A1000 Laptop GPU with 2560 CUDA cores and FP64/FP32 $X=(X_1,\dots,X_n)$2. On an ANN surrogate of the Peaks function, increasing the number of subdomains tightened root-node lower bounds significantly. Mean Value Form became tighter than Natural Interval Extension as $X=(X_1,\dots,X_n)$3 increased, consistent with the quadratic excess-width improvement. Branch-and-bound iterations decreased monotonically with $X=(X_1,\dots,X_n)$4, and Mean Value Form yielded fewer iterations than Natural Interval Extension across all tested $X=(X_1,\dots,X_n)$5 values (Zhang et al., 28 Jul 2025).

Wall-clock behavior reflected both interval accuracy and hardware effects. A pure-CPU serial Subdomain Lower Bounding scheme improved wall-clock time by approximately two orders of magnitude up to 64 subdomains relative to 1 subdomain, after which CPU performance plateaued. The heterogeneous GPU framework achieved over three orders of magnitude speedup against CPU interval arithmetic without partitioning. With $X=(X_1,\dots,X_n)$6, gains saturated; increasing to 4096 further tightened bounds, but wall-clock time stayed similar because of serialization and CPU sorting overhead (Zhang et al., 28 Jul 2025).

The two CUDA realizations were not equivalent in practice. The CUDA Graph implementation reduced lower-bounding time by approximately 20–90% compared with the single monolithic kernel, depending on $X=(X_1,\dots,X_n)$7 and problem complexity. The reported reason is the combination of reduced launch overhead and concurrency across independent DAG branches (Zhang et al., 28 Jul 2025).

Comparison with MAiNGO’s default McCormick-based solver is more nuanced than a simple dominance claim. For ANN surrogates of Alpine02 in 2–6 dimensions, GPU-parallel Subdomain Lower Bounding produced tighter lower bounds and fewer branch-and-bound iterations for 2–5 dimensions. The default solver sometimes had lower wall-clock time for $X=(X_1,\dots,X_n)$8 variables; however, at 5D it failed within a 2-hour limit whereas the proposed method converged within approximately 1 hour; at 6D both timed out, but the proposed method obtained a much tighter final lower bound. On constrained 2D Alpine02 validation problems with various linear and nonlinear constraints, interval bounds of constraints were used to prune infeasible subdomains and nodes, the method handled all constraint types, and it consistently required fewer branch-and-bound iterations. Without bound tightening, the proposed method was competitive or superior on several problems, including Styblinski–Tang minimization and kinetic ODE parameter estimation. At the same time, MAiNGO’s default solver gained substantial speedups from Duality-Based Bound Tightening and Constraint Propagation; Constraint Propagation remained compatible with Subdomain Lower Bounding, but its effect was often modest in the tested set (Zhang et al., 28 Jul 2025).

These results address a standard criticism of interval methods. Natural interval arithmetic on a single unsplit box can indeed be weak, but the reported solver does not rely on that regime. It couples splitting, midpoint-centered Mean Value Form, forward-mode AD, and GPU parallelism, and in some case studies it delivers performance competitive with or better than a McCormick-relaxation default solver (Zhang et al., 28 Jul 2025).

6. Relation to broader GPU-parallel lower-bounding research

The interval-based GPU-parallel lower bounding solver belongs to a broader family of GPU-accelerated bounding modules in exact optimization. In linear domain propagation for mixed-integer and pseudo-boolean optimization, a GPU-parallel feasibility-based bounds tightening method processes propagation rounds entirely on the GPU over CSR-stored sparse matrices, and reports geometric mean speed-ups around $X=(X_1,\dots,X_n)$9–$X_i=[\underline{X}_i,\overline{X}_i]$0 on reasonably large instances and up to $X_i=[\underline{X}_i,\overline{X}_i]$1 on favorably large instances (Sofranac et al., 2020). In complete search for bound-constrained nonlinear minimization, a GPU-based interval method uses single program, single data execution and variable cycling, and reports complete, rigorous enclosures of multimodal benchmark minima up to 10,000 dimensions on one GPU (Zhang et al., 2 Jul 2025).

Other lower-bounding architectures target different relaxations while retaining the same offloading principle. Batched first-order LP solving for MIP uses GPU-parallel PDHG to evaluate many related LP relaxations at once, with reported speedups of $X_i=[\underline{X}_i,\overline{X}_i]$2, $X_i=[\underline{X}_i,\overline{X}_i]$3, $X_i=[\underline{X}_i,\overline{X}_i]$4, and up to $X_i=[\underline{X}_i,\overline{X}_i]$5 in strong-branching case studies, and an average OBBT speedup of $X_i=[\underline{X}_i,\overline{X}_i]$6 against sequential dual simplex on neural-network verification instances (Blin et al., 29 Jan 2026). For exact sparse linear regression with an $X_i=[\underline{X}_i,\overline{X}_i]$7–$X_i=[\underline{X}_i,\overline{X}_i]$8 penalty, a GPU-accelerated nonlinear branch-and-bound framework solves an interval relaxation at each node by ADMM with closed-form coordinate-wise updates and reports GPU speedups of $X_i=[\underline{X}_i,\overline{X}_i]$9–$S$0 versus a CPU implementation of the same algorithm (Meng et al., 4 Feb 2026).

Constraint-programming lower bounds exhibit a similar pattern. For the Bin Packing constraint, GPU-accelerated dual feasible function lower bounds are used inside interval-based propagation, with the GPU implementation solving more instances or matching stronger baselines on several large-capacity benchmark families, and solving the Scholl 3 instances that the compared L2 and Arc-Flow propagators did not solve within the reported limit (Tardivo et al., 2024). In permutation flow-shop branch-and-bound, earlier GPU work had already shown the value of pooling many subproblems and parallelizing only the lower-bound evaluation, reaching peaks near $S$1 on a Tesla C2050 and $S$2 in a two-GPU adaptive variant [(Nouredine et al., 2012); (Chakroun et al., 2012)].

Taken together, these works suggest that GPU-parallel lower bounding is best understood as a design pattern rather than a single algorithmic template. The recurring principles are homogeneous kernel work, explicit management of memory hierarchy and transfer cost, and the retention of irregular search control on the CPU unless the bound itself admits a GPU-resident fixed-point or batched formulation.

7. Limitations and future directions

The main limitations of the subdomain interval solver are structural rather than incidental. Very high-dimensional problems make uniform partitioning infeasible because $S$3 grows too quickly; even adaptive strategies may not permit enough splitting to obtain strong tightness. Highly non-smooth or discontinuous functions degrade Mean Value Form tightness; the interval enclosures remain conservative, but the bounds may become weak. On consumer GPUs, FP64 throughput limitations induce serialization, and CPU-side aggregation of per-subdomain results can become the bottleneck as $S$4 increases. Large DAGs combined with large $S$5 can also create memory pressure when storing per-subdomain intermediates (Zhang et al., 28 Jul 2025).

The reported future directions are correspondingly targeted. Adaptive partitioning may be guided by bound sensitivity or error indicators rather than width alone; refinement may continue until a target bound tightness is reached; mixed CPU/GPU parallelization may process multiple nodes on CPU in parallel while each node invokes GPU Subdomain Lower Bounding; GPU-native reductions for refinement hulls and minima may replace CPU aggregation; interval libraries on GPUs may be extended with correctly rounded transcendental implementations and efficient interval AD; and hybrid strategies may combine Mean Value Form intervals with McCormick or semidefinite relaxations, selecting the bounding technique per node from DAG structure and estimated tightness (Zhang et al., 28 Jul 2025).

A plausible implication is that the long-term competitiveness of interval-based GPU bounding will depend as much on systems issues as on enclosure theory. The paper already shows that CUDA Graphs, memory scheduling, and reduction placement materially affect end-to-end performance. If those systems choices are combined with stronger adaptive partitioning and tighter GPU interval primitives, interval lower bounding could occupy a broader role in deterministic global optimization than its traditional single-box CPU implementations have allowed (Zhang et al., 28 Jul 2025).