arXiv CS open access

B$^3$-PWL: GPU-Batched Branch-and-Bound for Piecewise-Linear Optimization with SOS2 Constraints

2026 · arxiv_cs

arXiv CS · Papers · 2026

cloud distributed-computing parallel-computing

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

B3-PWL: GPU-Batched Branch-and-Bound for Piecewise-Linear Optimization with SOS2 Constraints

arXiv:2608.28988v1 [math.OC] 29 Aug 2026

Yilin Guan 1 3

Shuqing Luo 2

Pingzhi Li 2

Tianlong Chen 2,†

Kaidi Xu 3,4†

1 Johns Hopkins University 2 University of North Carolina, Chapel Hill City University of Hong Kong 4 The Hong Kong Institute of AI for Science

Abstract Piecewise-linear (PWL) optimization problems arise in many mixed-integer programming (MIP) optimization applications, including portfolio optimization, workforce scheduling, and resource allocation. But solving them to global optimality remains computationally expensive because branch-and-bound repeatedly solves LP relaxation subproblems. Existing solvers are largely CPU-centric, leaving the scalability of modern GPUs underutilized. Few prior GPU-accelerated branchand-bound either targets neural network which is not suitable for general PWL optimization, or accelerates only auxiliary subroutines such as strong branching heuristics within CPU-centric MIP solvers. To bridge this gap, we propose B3 -PWL, a GPU-centric batched branch-and-bound framework for piecewise-linear optimization with Special Ordered Set of type 2 (SOS2) constraints. Our method solves batches of LP relaxation subproblems concurrently on the GPU using a first-order primal-dual solver, enabled by a specialized batched block-tiled sparse matrix kernel. To complement bound computation, we further introduce a unified feasibility search module that combines an SOS2 repair primal heuristic with a batched feasibility pump to rapidly obtain feasible incumbents and improve pruning efficiency. On a benchmark of 43 PWL-MIP instances, B3 -PWL achieves a 9.25× geometric-mean speedup over NVIDIA cuOpt while reaching high-quality feasible incumbents on every tested instance. On a public valve-point unit-commitment benchmark, it further outperforms NVIDIA cuOpt and the open-source CPU solvers SCIP and HiGHS, demonstrating the potential of first-order LP methods as the central engine of GPU-accelerated branch-and-bound.

Introduction

Combinatorial optimization with piecewise-linear (PWL) structure underlies many real-world decision-making tasks, including portfolio optimization [Lobo et al., 2007], workforce scheduling [Ernst et al., 2004], and large-scale resource allocation [Karimi et al., 2003]. In mixed-integer programming (MIP), PWL formulations are commonly used to approximate nonlinear relationships while preserving compatibility with optimization frameworks. However, solving such problems to global optimality remains computationally expensive, because the integer constraints induce a combinatorial search space that grows exponentially with problem size. The standard branch-and-bound (B&B) method tackles this by recursively partitioning the search space and solving an LP relaxation at every node it explores, accumulating thousands to millions of LP solves on a single instance [Land and Doig, 2009, Nemhauser and Wolsey, 1988]. In conventional MIP solvers, these LP relaxations are handled by simplex-based methods such as dual simplex [Dantzig, 1990], whose basis re-optimization † Corresponding authors. Contact: yguan29@jh.edu, tianlong@cs.unc.edu, kaidixu@cityu.edu.hk.

Preprint.

and sequential pivoting are difficult to parallelize. As a result, although modern GPUs have utilized massive parallel computing cores to revolutionize large-scale machine learning, the same accelerated computing potential remains largely unexploited for PWL optimization problems. GPU-accelerated branch-and-bound has been studied in two adjacent settings. The first line of work targets neural network verification, where frameworks such as αβ-CROWN Xu et al. [2021], Wang et al. [2021] run B&B almost entirely on the GPU. Their efficiency relies on two properties specific to neural network verification. First, neural networks have a layered structure that supports GPU-friendly bound propagation methods such as LiRPA Xu et al. [2020], which compute linear lower bounds on the network output layer by propagating linear relaxations through the layers, avoiding LP solves at each node. Second, verification only asks whether the worst-case network output violates a safety property, not what its exact value is. Each LP subproblem therefore does not need to be solved exactly: a loose lower bound is enough to prune or further branch the subproblem. This leads to efficient pruning that drastically reduces the number of nodes explored in search tree. However, PWL optimization is more general in structure and it requires solving each LP subproblem precisely to get the globally minimal objective value. Obtaining tight enough bounds and preventing the search tree from exploding are therefore both substantially harder in PWL optimization. The second line of work uses GPU-accelerated LP solving as an auxiliary component of MIP solvers. NVIDIA cuOpt NVIDIA [2025], Blin et al. [2026] apply GPU-batched first-order LP methods to specific subroutines such as full strong branching (FSB) [Achterberg et al., 2005] and optimizationbased bound tightening (OBBT) Gleixner et al. [2017]. However, the LP relaxation at each node of the main search tree is still solved on the CPU using simplex-based methods such as dual simplex. The GPU here only accelerates expensive auxiliary procedures rather than the core B&B loop. In this paper, we propose B3 -PWL, a GPU-batched branch-and-bound framework for piecewise-linear optimization problems with Special Ordered Set of type 2 (SOS2) constraints. Unlike prior work that relegates GPU LP solving to auxiliary subroutines, B3 -PWL puts batched LP relaxations of subproblems on the GPU as the central compute primitive of the B&B loop. To support this paradigm, we extend the first-order method cuPDLPx Lu et al. [2025] to concurrently solve multiple nodespecific LP relaxations within a single batched dispatch, backed by a custom block-tiled sparse matrix kernel that handles the ragged constraint matrices. In addition, we introduce a unified feasibility search module that integrate an SOS2 repair primal heuristic [Keha et al., 2006] and batched feasibility pump strategy tightly and seamlessly into our solving pipeline. It effectively and rapidly repairs and finds feasible solutions to update the best incumbent and tighten the global upper-bound for efficient search tree pruning. Our contributions are as follows: • We present B3 -PWL, a GPU-batched branch-and-bound framework for piecewise-linear optimization with SOS2 constraints, in which routine node LP relaxations are batched and solved on the GPU via a modified first-order solver for batched LP problems. • We develop a unified feasibility search module combining an SOS2 repair primal heuristic and a batched feasibility pump to rapidly obtain feasible incumbents for efficient pruning. • On a 43-instance benchmark, B3 -PWL achieves a 9.25× geometric-mean speedup and obtain higher quality solutions over cuOpt. We further evaluate on a public valve-point unitcommitment benchmark, where B3 -PWL outperforms cuOpt, and open-source CPU solvers SCIP and HiGHS, highlighting the potential of first-order LP methods for GPU-accelerated branch-and-bound.

Related Work

First-order methods for linear programming. Simplex algorithm [Dantzig, 1990] and interiorpoint methods [Wright, 1997, Karmarkar, 1984] have long dominated linear programming solving. While these algorithms remain the formidable foundation of state-of-the-art CPU solvers, their reliance on matrix factorizations and sequential computations limits their scalability on GPUs. To unlock GPU acceleration for large-scale optimization, first-order methods (FOMs) have emerged as a scalable alternative. PDLP [Applegate et al., 2021, 2025] reformulates LP as a saddle-point problem and solve it via primal-dual hybrid gradient (PDHG) iterations [Chambolle and Pock, 2011] with a range of practical enhancements, which transforms FOMs into a practical, numerically stable 2

paradigm for large-scale LPs. To leverage the massive memory bandwidth of modern accelerators, cuPDLP [Lu and Yang, 2025, Lu et al., 2023] extends the CPU-based PDLP to GPUs by offloading core sparse matrix computations and introducing hardware-aware heuristics. cuPDLPx [Lu et al., 2025] further accelerates this paradigm by tailoring the algorithm, achieving multi-fold speedups over its predecessor on large-scale benchmarks. Nevertheless, these works focus primarily on accelerating individual LP solves, rather than integrating batched LP computation into B&B. GPU-accelerated branch-and-bound. GPU acceleration of B&B has been particularly successful in neural network verification, where αβ-CROWN [Xu et al., 2021, Wang et al., 2021] pioneered this problem by combining branch-and-bound with GPU-accelerated linear bound propagation (e.g., LiRPA) [Xu et al., 2020] and tightening loose bounds with batched gradient descent on GPUs. In traditional mixed-integer programming (MIP), however, the use of GPUs within exact LP-based branch-and-bound remains more limited. Early efforts incorporating GPUs, such as Neural Branching [Nair et al., 2020] uses a GPU-scalable ADMM-based variant of full strong branching to generate offline imitation-learning targets for branching policies. More recently, Blin et al. [2026] extended batched primal-dual hybrid gradient (PDHG) [Chambolle and Pock, 2011] methods to execute online, exact LP batches on GPUs. However, their application is limited to accelerating specific heuristic sub-routines, including full strong branching (FSB) [Achterberg et al., 2005] and optimization-based bound tightening (OBBT) [Gleixner et al., 2017]. In contrast, we use batched GPU LP solving as the core engine for routine node processing in exact B&B for PWL optimization with SOS2 constraints.

Preliminaries

3.1

PWL Problem Formulation with SOS2 Constraints

We consider an optimization problem in variables x ∈ Rn , where the first J components {xj }Jj=1 are linked to PWL functions. Each PWL function fj (xj ) is defined over the interval [ℓj , uj ] with Tj + 1 ordered breakpoints {dj,0 , dj,1 , . . . , dj,Tj }, partitioning the domain into Tj linear segments with corresponding values fj (dj,i ). Using the λ-formulation with Special Ordered Sets of type 2 (SOS2), we formulate the problem as: min

c⊤ x + e ⊤ z

s.t.

Ax + Gz ≤ b,

x,z,λ

(1) ℓ ≤ x ≤ u,

λj,i = 1,

∀ j ∈ [J],

(2)

i=0

xj =

Tj X

dj,i λj,i ,

zj =

i=0

Tj X

fj (dj,i ) λj,i ,

∀ j ∈ [J],

(3)

i=0

λj,i ≥ 0, ∀ i, j, SOS2(λj,· ), ∀ j ∈ [J]. T +1

Here, zj ∈ R denotes the value of fj (xj ), and λj,· = (λj,0 , . . . , λj,Tj ) ∈ R+j are convex combination weights. The SOS2 constraint requires that λj,· to have at most two nonzero components at adjacent indices. That is, if λj,i > 0 and λj,i′ > 0 with i ̸= i′ , then |i − i′ | = 1. This ensures that (xj , zj ) lies on a single linear segment or exactly on a breakpoint of the PWL graph of fj . 3.2

Branch and Bound (B&B) for SOS2 Constraints

The combinatorial difficulty of the PWL optimization arises from the SOS2 constraints. Each PWL function fj admits several candidate active segments, and identifying the correct one for every fj requires searching over an exponentially large combination of segment choices. Branch-andbound (B&B) solves this problem by recursively restricting the admissible breakpoint range for each PWL function. At the root node, all breakpoints are admissible. Each node D is characterized by per-function index intervals [Lj , Uj ] ⊆ [0, Tj ] and adds the restriction λj,i = 0,

∀ j ∈ [J], i ∈ / {Lj , . . . , Uj }. 3

(4)

CPU

GPU ③ Batched cuPDLPx

↻ outer loop

① Branch-and-Bound Search Tree

Stage 1 Batched PSLP Presolve + Ruiz Scaling

explored pruned active leaf

input batch

homogeneous constraint matrices

✖

② Pack into Batch LP2

LP3

⋯

LPB

batch

LP LB ≥ global UB, or LP infeasible

Update ↑

best feasible solution → update global UB

� {�(t) � }�=1

B²-SpMV �(t)⊤ ⋅ �(t) �

� {�(t) � }�=1

B²-SpMV �(t) ⋅ �̄(t) �+1

�̄(t) �+1

Π[ℓ, u]

2 �(t) �+1

↻ k ← k+1 iters

Π�

−

�(t) �

� {�(t) �+1 }�=1 � {�(t) �+1 }�=1

Stage 3 Batched PSLP Postsolve + Rescaling

(�(t)∗ , �(t)∗ ) ← postsolve(�̃(t) , ỹ(t) )

④ Verify Solutions & Update Tree

Prune ✖

heterogeneous constraint matrices

Stage 2 Batched PDHG

N active leaves selected

LP1

presolved batch

presolve + scaling

Branch ⋎

for � = 1, … , �

Batch of LP Solutions

{ ( �(t)∗ , �(t)∗ ) }

split active breakpoint range into two children

for � = 1, … , �

Figure 1: Overall workflow for B3 -PWL. Each round alternates between CPU-side B&B management and GPU-side batched LP solving. The CPU (left) selects N active subproblems from the search 1 packs their LP relaxations into a batch (⃝), 2 and after the GPU returns solutions, verifies tree (⃝), 4 The GPU feasibility and updates the tree by pruning, updating the incumbent, or branching (⃝). 3 per-task PSLP presolve and Ruiz scaling produce (right) executes a batched cuPDLPx pipeline (⃝): heterogeneous constraint matrices, batched PDHG iterations (built on our B2 -SpMV kernel) run until convergence, and postsolve and rescaling yield primal-dual solutions for the entire batch. The outer loop continues until the active queue is empty or the global gap closes. The subproblem at D is therefore the original SOS2 formulation (1)–(3) together with the interval restriction (4). Its LP relaxation, obtained by dropping the SOS2 constraints, provides a lower bound for D. Three cases may arise after the LP problem is solved: (i) if the LP is infeasible or its optimum is no better than the incumbent, D is fathomed; (ii) if the LP solution satisfies all SOS2 constraints, it updates the incumbent if improved and D is fathomed; (iii) otherwise, D is branched by splitting one active range into two children. The algorithm terminates when the active queue is empty or the global gap falls within tolerance.

Method

We introduce B3 -PWL, a GPU-batched branch-and-bound framework for piecewise-linear optimization with SOS2 constraints. In our framework, batched LP solving is the central compute primitive: the CPU organizes the global search state and assembles task batches, while the GPU executes batches of LP tasks and returns primal and dual solutions. We describe the end-to-end workflow in Section 4.1, then present the batched node-LP solver in Section 4.2, the block-tiled sparse backend for heterogeneous matrices in Section 4.3, and the unified feasibility-search module in Section 4.4. 4.1

Overall Workflow

Different node subproblems within a single PWL B&B tree descend from a common root formulation and differ only in their variable bounds (the active breakpoint ranges), which is naturally suited to batching. Exploiting this, we extend the GPU LP solver cuPDLPx [Lu et al., 2025] to concurrently solve multiple node LP relaxations within a single batched dispatch. Building on this batched solver, our framework proceeds in rounds of three stages. (1) In the assembly stage, the CPU selects active nodes from the active nodes queue according to the current B&B policy and packs their LP relaxations into a single GPU batch. (2) In the solve stage, the GPU solves all LP relaxations in the batch in parallel through batched matrix operations, and returns their primal-dual solutions. (3) In the postprocess stage, the CPU uses the returned LP solutions to update the global lower and upper bounds, prune fathomed nodes, and branch unresolved nodes 4

into new subproblems that are enqueued into the active node queue. Beyond these standard B&B operations, we further exploit the LP solutions through a unified feasibility-search module which aims at recovering high-quality feasible incumbents that tighten the upper bound (Section 4.4). The next round then assembles its batch from the updated queue, and the loop continues until the active queue is empty or the global gap falls within tolerance. 4.2

Heterogeneous Batched LP Engine on GPU

Each batch consists of B tasks T = {1, . . . , B}, each corresponding to one node LP relaxation. (t) (t) SOS2 constraints are dropped, and the active breakpoint range [Lj , Uj ] is encoded through the (t)

(t)

fixings λj,i = 0 for i ∈ / [Lj , Uj ]. The objective c(t) is either the original node objective for regular B&B search tasks or a feasibility-pump objective for feasibility tasks (see Section 4.4). Each task is independently transformed by task-local presolve and Ruiz rescaling into the reduced-space boxed LP form (t) (t) (t) min c(t)⊤ x(t) s.t. ℓA ≤ A(t) x(t) ≤ uA , ℓ(t) ≤ u(t) (5) x ≤x x , x(t)

yielding reduced matrices A(t) that differ across tasks in dimensions and sparsity patterns. Unlike the prior batched LP engine for heuristic subroutines in cuOpt, where approximate estimation allows presolve to be skipped so that all tasks share a common constraint matrix A, our setting requires presolve and therefore motivates a dedicated batched SpMV implementation for heterogeneous matrices A(t) (Section 4.3). Building on cuPDLPx, our engine batches its inner PDHG operator across all B tasks in a single GPU dispatch, while task-local stepsizes, primal weights, projections, restart states, and stopping criteria are maintained independently. The surrounding reflected Halpern updates and adaptive-restart logic follow cuPDLPx and are described in Appendix C. The inner PDHG operator is dominated by two sparse matrix–vector products. (t)

(t)

B Let X (t) = [ℓx , ux ] and S (t) = [ℓA , uA ], and let Xk = {xk }B t=1 and Yk = {yk }t=1 denote the batch-wise primal and dual iterates at PDHG iteration k. For each task t, the inner PDHG operator takes the form (t) (t) (t) (t) xk+1 = ΠX (t) xk − τk c(t) − A(t)⊤ yk , (6) (t)

(t)

x̄k+1 = 2xk+1 − xk ,

(7)

(t) (t) (t) (t) (t) yk+1 = yk − σk A(t) x̄k+1 − σk Π−S (t) (t)

(t)

(t) (t) (t) (σk )−1 yk − A(t) x̄k+1

(8)

(t)

Here, τk = ηk /ωk and σk = ηk ωk are the task-local primal and dual stepsizes induced by (t) (t) the cuPDLPx stepsize ηk and primal weight ωk . The row-bound proximal operator projects onto (t) (t) −S (t) = [−uA , −ℓA ]. Full details of cuPDLPx, including the reflected Halpern scheme, adaptive restarting, and primal-weight updates, are provided in Appendix C. 4.3

Accelerating Ragged SpMV with Batched Blocking

flag matrices (bool) index matrices (uint8) A key challenge in our approach is that different B&B nodes may produce LP relaxations with different numbers of active variables and constraints after pre-process presolve, delivering a ragged collection of sparse matrices {A(t) }B t=1 with heterogeneous dimensions. stacked tiles Standard GPU-accelerated sparse linear algebra li(float64) braries such as cuSPARSE NVIDIA Corporation [2024] provide no batched SpMV primitives. Sequentially launching an SpMV kernel per subproblem fails e siz k oc to fully utilize the GPU parallel resources, while flatbl block size tening the batch into a single block-diagonal sparse matrix leads to poor load balance and inefficient in- Figure 2: Building stacked tiles for dexing overhead. To support efficient batching on B2 -SpMV. such heterogeneous workloads, we introduce a batched blocking scheme for efficient and parallel 5 4 3 2 1 0

solving, dubbed B2 -SpMV. It partitions each sparse matrix A(t) into fixed-size tiles, removes all-zero tiles, and packs the remaining nonzero tiles from all subproblems into a single contiguous GPU tensor. For each subproblem, we build a tile index map that identifies, for every tile position, either the location of the corresponding nonzero tile in the packed tensor or an empty marker, enabling both sparse computing for each LP sub-problem and coalesced memory access for efficient GPU execution. B2 -SpMV executes both A(t) x(t) and A(t)⊤ y (t) with a single batched kernel over the entire ragged batch, where each thread block processes one tiled row or column for one subproblem, with empty tiles skipped based on the index map. 4.4

Unified Feasibility Search for Efficient Upper-Bound Tightening

Efficient pruning in B&B depends on tightening either the lower bounds for each node or the global upper bound from the best feasible incumbent. Therefore, high-quality incumbents are critical for pruning by bounding and help drastically reducing the search space. However, LP solutions typically violate SOS2 constraints. We therefore introduce a feasibility-search module that attempts to convert LP primal solutions into feasible incumbents through two-stage procedure: a fast and direct SOS2 repair primal heuristic, followed by a batched feasibility pump when direct repair fails. SOS2 repair primal heuristic. Consider an LP solution (x̂, ẑ, λ̂) at node D. We attempt to repair it into an SOS2-feasible one by re-interpolating λ̂j,· for each PWL function. The idea is to use x̂j to identify the active segment it lies on and reset λ̂j,· to the convex-combination weights of that segment’s two endpoints. L

Concretely, for each PWL function j with active breakpoint range [Lj , Uj ], we clip x̂j to [dj j , dj j ], k

k +1

locate the current segment kj ∈ {Lj , . . . , Uj − 1} with dj j ≤ x̂j ≤ dj j

, and reset λ̃j,· as the

k k +1 convex-combination weights of dj j and dj j : k +1

λ̃j j =

dj j

− x̂j

k +1

− dj

dj j

k +1

λ̃j j

, kj

x̂j − dj j k +1

dj j

− dj j

P k k with λ̃kj = 0 for all other k ̸= kj , kj + 1. The repaired outputs follow as x̃j = k dj λ̃j and P z̃j = k fj (dkj ) λ̃kj , yielding a candidate (x̃, z̃, λ̃) that satisfies SOS2 by construction. Since the repaired candidate satisfies SOS2 by construction, it remains to verify the original linear constraints Ax̃ + Gz̃ ≤ b and variable bounds. If satisfied, (x̃, z̃, λ̃) is feasible for the original PWL problem and updates the incumbent upper bound. Batched Feasibility-Pump (FP). When first-stage repair fails, the repaired vector λ̄(r) still encodes which SOS2 segment each x̂j most naturally suggests. The idea of objective feasibility pump (OFP) Achterberg and Berthold [2007] is to treat λ̄(r) as a soft target and re-solve the LP under a combined objective that balances the original PWL minimization objective with proximity to λ̄(r) , so that successive LP solves are pulled toward an SOS2-feasible point while remaining close to optimality, until repair eventually succeeds. We adapt OFP to our pipeline in two ways: the objective is tailored to the SOS2 formulation, and FP tasks are integrated into the same GPU batch with our B&B search LP tasks. Specifically, each pump round solves the same node LP with the original constraints intact but a modified objective: (r)

cFP = αr c(r) prox + (1 − αr ) γr corig . (r)

(r)

where cprox has entries (cprox )j,k = −λ̄j,k on the λ block and zero elsewhere, so that minimizing it rewards λ-mass on the two breakpoints selected by repair, thus pulls the LP solution toward the SOS2 target λ̄(r) . corig is the original objective matrix and keeps the solution close to optimal. αr ∈ [0, 1] (r) trades off feasibility against quality, and the rescaling factor γr = ∥cprox ∥2 /(∥corig ∥2 + ε) aligns their ℓ2 magnitudes so αr remains effective. We set α0 = 1 and decay αr+1 = 0.9 αr on each failure. Since FP tasks differ from search tasks only in their objective vectors and share the same constraint structure, our batched cuPDLPx solver, which supports per-task objectives, can dispatch them in the 6

same GPU batch as regular B&B search tasks at no additional scheduling cost. After each FP solve, the returned primal is repaired and checked: accepted candidates update the incumbent, otherwise the repaired point becomes λ̄(r+1) for the next round.

Experiments

5.1

Experimental Setup

Benchmarks. Public PWL-MIP instances that keep SOS2 constraints in native form are limited (detailed in Appendix G), e.g. MIPLIB 2017, the most widely used MIP benchmark library, contains no instance with SOS or piecewise-linear constraints, even though the MPS format supports them [Gleixner et al., 2021]. Recent PWL work also relies on its own instance generators [Hübner et al., 2026]. We therefore evaluate B3 -PWL on two benchmarks. The first is a suite of 43 synthetic PWL-MIP instances from three application families: portfolio optimization, workforce scheduling, and resource allocation. The suite contains 18 portfolio, 9 workforce, and 16 resource-allocation instances, with 100–1,200 PWL functions and 15–25 breakpoints per function. We can control LP size, breakpoint count, sparsity, and SOS2 density through generator to characterize the behavior of our GPU B&B framework across different scales and structures. Appendix G reports the formulations, generator parameters, seeds, and per-instance size statistics. The second is the public valve-point unit commitment benchmark of Pedroso et al. [2014] built on real-world cost data. This benchmark considers thermal unit commitment, where the objective is to determine unit commitment and generation schedules that satisfy multi-period demand while minimizing valve-point fuel costs. Hardware Settings. All end-to-end solver comparisons are run on the same Blackwell node, equipped with NVIDIA RTX PRO 6000 Blackwell Server Edition GPUs (96 GB each) and a dualsocket x86-64 host with 192 physical cores (384 logical threads) and 718 GB RAM. Solver Configurations. Both solvers are run with a 1% relative MIP gap and a 300-second wallclock limit. B3 -PWL uses a PDLP feasibility/optimality tolerances of 10−4 . The feasibility-pump module is enabled in the main configuration. cuOpt is run with v26.02.00 which combines CPUside dual simplex for node LPs, PaPILO presolve, MIR / Gomory / knapsack cutting planes, and GPU-accelerated batched LP for FSB and OBBT. Metrics.

Our evaluation focuses on two aspects: solution quality and solve speed.

Solution quality. we use Gurobi as an external oracle reference. We run Gurobi with tight tolerance ⋆ gap (= 10−6 ) and denote its reference objective by zref . We evaluate solvers on (1) Exact Match counts instances on which the solver returns the objective value that exactly matches the solution of Gurobi. (2) Eff-Opt counts instances on which the solver’s objective value zsolver satisfying ⋆ |zsolver − zref | −2 . ⋆ |} ≤ 10 max{1, |zref Thus, Eff-Opt corresponds to finding a feasible incumbent within 1% of the Gurobi reference objective. Solve time. All solvers terminate upon reaching a 1% MIP gap, exhausting the active search queue, or hitting the 300 s wall-clock budget, whichever comes first; total time is the wall-clock time from solver start to termination. For aggregate runtimes, we report shifted geometric means with shift s = 1 s following standard MIP benchmarking practice [Achterberg, 2007, 2009]: ! n 1X GeoTime(t) = exp log(ti + s) − s. n i=1 The shift reduces the influence of trivially easy instances and avoids overweighting sub-second timing differences. Arithmetic means report the average of wall-clock times. Geo-Spd is the geometric mean of per-instance speedup ratios: ! n cuOpt X 1 t log Bi3 -PWL . GeoSpd(cuOpt/B3 -PWL) = exp n i=1 ti 7

Table 1: Solution quality and solve time of B3 -PWL vs. cuOpt on 43 PWL-MIP instances. Eff-Opt Category

Exact

Arith. Time (s)

Geo. Time (s)

Geo-Spd

n B3 -PWL cuOpt B3 -PWL cuOpt B3 -PWL cuOpt B3 -PWL cuOpt (cu/B3 -PWL)

B3 -PWL Wins

Portfolio 18 Workforce 9 Resource Alloc. 16

18 9 16

10 8 12

15 9 13

0 1 8

2.5 2.6 82.3

142.2 68.9 107.9

2.1 2.2 39.0

57.0 35.7 67.3

29.35× 17.69× 1.75×

18/18 9/9 11/16

Total

32.2

114.1

7.0

55.0

9.25×

38/43

Table 2: Comparison with established CPU MIP solvers on the same 43 instances. Solver 3

B -PWL Gurobi 13.0.2 SCIP 10.0.1 HiGHS 1.15.1 cuOpt 26.02

Eff-Opt

Shifted Geo. Time (s)

43/43 43/43 37/43 40/43 30/43

7.0 1.0 27.1 25.2 55.0

with values > 1 indicating B3 -PWL is faster on average. We additionally report Time to Best Incumbent (TTB): the wall-clock time at which a solver first reaches the incumbent it ultimately returns, censored at 300 s for runs that fail to reach Eff-Opt. 5.2

Main Results

Table 1 reports the performance comparisons between B3 -PWL and cuOpt across our PWL benchmark. Detailed per-instance objectives and wall-clock times are provided in Appendix E. Solution quality. B3 -PWL reaches Eff-Opt on all 43 instances, compared with 30 / 43 for cuOpt, and matches the Gurobi reference objective (Exact Match) on 37 / 43 instances, compared with 9 / 43 for cuOpt. These results indicate that B3 -PWL consistently returns higher-quality feasible incumbents than cuOpt on our benchmarks. Solve speed. In total wall-clock time, B3 -PWL achieves a shifted geometric mean of 7.0 seconds, compared with 55.0 seconds for cuOpt, and an arithmetic mean of 32.2 seconds, compared with 114.1 seconds for cuOpt. Under the reported Geo-Spd metric, this corresponds to a 9.25× aggregate speedup. We say B3 -PWL wins on an instance if it reaches Eff-Opt while cuOpt does not, or if both reach Eff-Opt and B3 -PWL terminates faster. Under this criterion B3 -PWL wins on 38 / 43 instances. Comparison with established CPU solvers. Table 2 places these results in the context of mature CPU MIP solvers. Gurobi remains the fastest solver on this benchmark, reaching Eff-Opt on all 43 instances with a shifted geometric-mean time of 1.0,s, compared with 7.0,s for B3 -PWL. B3 -PWL nevertheless reaches Eff-Opt on all 43 instances and achieves substantially lower aggregate runtime than the open-source CPU solvers SCIP and HiGHS, which reach Eff-Opt on 37/43 and 40/43 instances, respectively. These results position B3 -PWL as a proof of concept that batched firstorder LP solving can serve as the core engine of GPU-centric branch-and-bound, rather than as a replacement for mature commercial MIP solvers. 5.3

Public valve-point UC benchmark

We further evaluate B3 -PWL on the public valve-point unit commitment benchmark of Pedroso et al. [2014]. Since valve-point costs are nonconvex and nonsmooth, the original paper solves the problem through an adaptive piecewise-linear approximation, yielding a sequence of PWL-MIP instances. We use the public 40-unit eld40/ucp40 data with the original horizons T ∈ {1, 3, 6, 12, 24}. Unlike the adaptive breakpoint-refinement procedure of the original paper, we fix each cost curve to a static breakpoint count (20 unless noted) to obtain deterministic solver-comparison instances. As shown in Table 3, Gurobi remains fastest and we do not claim otherwise. Against the open-source simplex-based solvers B3 -PWL is faster on 6 of 7 instances (geomean 10.3× over SCIP, 5.1× over 8

Table 3: Time to a 1% relative MIP gap (s) on the valve-point UC benchmark. Rows 1-5 sweep the horizon T at 20 breakpoints; rows 5-7 sweep breakpoints at T = 24. Geo-Spd is the geometric mean of per-instance ratios against B3 -PWL. B3 -PWL

Gurobi

SCIP

HiGHS

cuOpt

T =1 T =3 T =6 T = 12 T = 24, bp=20

1.13 3.47 2.12 2.92 5.93

4.02 0.45 0.61 1.19 3.04

0.12 37.8 21.1 36.4 >300

0.06 211.4 11.0 14.2 39.9

3.33 245.9 10.9 209.0 >300

T = 24, bp=15 T = 24, bp=30

5.73 9.02

2.36 3.95

>300 >300

34.6 254.0

>300 >300

—

0.47×

10.34×

5.13×

25.40×

Instance

Geo-Spd

Table 4: Final relative gap (%) and best objective on the valve-point UC benchmark. “n/a”: no feasible solution within the budget. B3 -PWL

Gurobi Obj.

Gap

SCIP Obj.

Gap

HiGHS Obj.

Gap

cuOpt

Instance

Gap

Obj.

Gap

Obj.

T =1 T =3 T =6 T = 12 T = 24, bp=20

0.052 122,356 0.485 122,889 0.258 122,609 0.010 122,305 0.086 122,401 0.999 293,101 0.990 293,375 0.929 293,394 0.460 292,890 0.876 293,266 0.544 587,487 0.371 586,675 0.634 588,350 0.288 586,440 0.436 586,981 0.327 1,192,044 0.673 1,196,353 0.827 1,198,228 0.289 1,191,898 0.905 1,199,171 0.369 2,401,839 0.182 2,397,989 6.781 2,556,152 0.262 2,400,451 5.721 2,538,496

T = 24, bp=15 0.345 2,404,897 0.201 2,402,118 6.242 2,547,243 0.233 2,403,133 1.682 2,437,913 T = 24, bp=30 0.414 2,398,461 0.189 2,393,431 n/a 0.115 2,391,921 1.894 2,434,827

HiGHS). The exception is T = 1, too small to amortize GPU setup. The margin widens with horizon length. At T = 24, SCIP and cuOpt exhaust the budget at every breakpoint setting while B3 -PWL certifies the gap in under 10 s. 5.4

Time to Best Incumbent

Total wall-clock time conflates two phases of branch-and-bound: finding a high-quality feasible incumbent and closing the final dual-bound gap. We separate them by reporting Time to Best Incumbent (TTB) (defined in Section 5.1), censored at 300 s for runs that do not reach Eff-Opt. Since B3 -PWL reaches Eff-Opt on every instance, its TTB is observed in all 43 cases. cuOpt’s TTB is censored on the 13 instances where it fails to reach Eff-Opt within the budget. Table 5 reports the comparison. Across all 43 instances, B3 -PWL achieves a 12.28× geometric-mean TTB speedup over cuOpt overall (4.22 s vs. 45.77 s shifted geomean), and a strictly lower TTB on every instance. Per-category speedups follow the same pattern as total runtime in Table 1: 28.66× on portfolio and 17.56× on workforce, narrowing to 3.87× on resource allocation, where cuOpt finds incumbents quickly via dual simplex. Comparing Table 5 with the total solve times in Table 1 reveals B3 -PWL’s end-to-end advantage stems primarily from rapid primal progress: high-quality incumbents are typically found within the first few GPU batches, after which most of the remaining wall-clock budget is spent on dual-bound tightening to certify optimality. This profile suggests that further end-to-end gains will come from accelerating dual-side convergence such as tighter LP relaxations. 5.5

Unified Feasibility Search Module Ablation

We ablate the two stages of the feasibility-search module in Section 4.4, the SOS2 repair primal heuristic and the batched Feasibility Pump (FP), under three configurations: (i) no feasibility heuristic, (ii) repair only, and (iii) repair + FP. Table 6 reports solve time across the 43-instance benchmark. 9

Table 5: Time to best incumbent (TTB) on the 43-instance benchmark. Geo. TTB is the shifted geometric mean with shift s = 1 s; Geo-Spd is the geometric mean of per-instance cuOpt/B3 -PWL TTB ratios. Geo. TTB (s)

Arith. TTB (s)

n B -PWL cuOpt Geo-Spd B3 -PWL cuOpt

Category

Feas

Geo. (s)

Arith. (s)

Geo. (s)

Arith. (s)

(ii)/(iii)

Portfolio Workforce Resource Alloc.

18 9 16

0/18 0/9 0/16

300.0 300.0 39.8

4.97 6.10 38.98

36.2 37.2 82.3

2.09 2.17 38.98

2.5 2.6 82.3

2.18× 2.61× 1.00×

Total

0/43

142.1

11.56

53.6

7.05

32.2

1.69×

SOS2 repair primal heuristic is necessary for PDLP outputs. Unlike simplex-based solvers, which return vertex solutions that naturally satisfy integrality at leaf nodes, PDHG iterates in PDLP converge from the interior of the LP polytope. Its postsolved primal iterate may satisfy the LP relaxation within global tolerance while still leaving small nonzero mass on SOS2 coefficients that should be exactly zero. Consequently, even leaf nodes whose branching restrictions identify a single admissible SOS2 segment may fail the strict incumbent check. Column (i) of Table 6 confirms this empirically: no feasible incumbent is recorded on any of the 43 instances. The SOS2 repair step resolves the mismatch by projecting each PWL component onto an admissible segment. With repair enabled (column (ii)), feasibility is recovered on all 43 instances at low overhead (Table 7). FP provides category-dependent acceleration. The Feasibility Pump injects proximity-objective LP tasks into the GPU batch to accelerate discovery of high-quality incumbents. Its benefit depends on whether early incumbent quality is the bottleneck for pruning. On portfolio and workforce instances, where the first high-quality incumbent dominates total solve time, FP yields geometricmean speedups of 2.18× and 2.61× respectively over repair alone. On resource allocation, however, the branching heuristic already discovers incumbents cheaply at nearby leaves, so FP yields no measurable acceleration (1.00×). Accelerating Ragged SpMV with Batched Blocking 2

We benchmark the efficiency of B -SpMV with varied batch size, where the size of the constraint matrix per sample is randomly generated within 8000 ∼ 25000. We use the default SpMV interface from cuSPARSE as baseline, with looped execution for batching. It scales with batch size and saturates beyond 50× at larger batches, confirming that the kernel effectively amortizes per-task launch and memory overhead across the batch.

Speedup (B 2 −SpMV vs. Looped cuSPARSE)

5.6

57.51×60.94×

55.52×

Time Breakdown

29.25×

30 20

Table 7 reports the time-weighted breakdown across the 43 instances. Batched LP solving (PDHG loop, task-local PSLP presolve and rescaling) consumes over 85% of B&B 10

48.70×

39.76×

16.23× 8.66×

5.7

53.66× 56.06×

1.21×1.79×

2.51×

4.42×

64 128 256 512 1024 2048 4096 8192

Batch Size (Benchmarking on NVIDIA H200 GPU)

Figure 3: Speedup Visualization for B2 -SpMV vs. Looped cuSPARSE.

3 Table 7: Time-weighted P component P BaB breakdown of B -PWL’s branch-and-bound phase. Each entry reports / i ti within the corresponding category. i ti

Component PDHG loop & GPU setup PSLP presolve/rescale SOS2 check & primal heuristic Branch selection & node expansion Batch construction Other postprocessing

Portfolio Workforce Resalloc Overall 82.1% 9.1% 2.0% 4.4% 3.8% 0.1%

70.4% 15.9% 1.0% 4.9% 4.7% 3.2%

73.9% 11.4% 6.0% 3.7% 3.3% 2.0%

74.0% 11.3% 5.9% 3.8% 3.3% 2.0%

time, while all CPU-side operations (SOS2 checking, primal heuristics, branching, batch construction, and postprocessing) account for less than 15% combined. This pattern holds consistently across all three categories, indicating that B3 -PWL is bound by GPU-side LP throughput rather than CPU-side tree management.

Conclusion

We introduced B3 -PWL, a GPU-batched branch-and-bound framework for PWL optimization with SOS2 constraints that solves batched node LP relaxations on the GPU as the main engine of the B&B loop. On a 43-instance benchmark, B3 -PWL achieves a 9.25× geometric-mean speedup over NVIDIA cuOpt while reaching high-quality incumbents on every tested instance. We see this work as a first step toward using batched first-order method for LP as the core engine of B&B, beyond its current role as auxiliary acceleration. Extending the paradigm to broader PWL and MIP settings, including the incorporation of cutting planes for tighter dual bounds and more general branching heuristics for efficient tree exploration, remains an open direction for future work.

References T. Achterberg. Constraint integer programming. PhD thesis, 2007. T. Achterberg. Scip: solving constraint integer programs. Mathematical Programming Computation, 1(1):1–41, 2009. T. Achterberg and T. Berthold. Improving the feasibility pump. Discrete Optimization, 4(1):77–86, 2007. T. Achterberg, T. Koch, and A. Martin. Branching rules revisited. Operations Research Letters, 33 (1):42–54, 2005. D. Applegate, M. Díaz, O. Hinder, H. Lu, M. Lubin, B. O’Donoghue, and W. Schudy. Practical large-scale linear programming using primal-dual hybrid gradient. Advances in Neural Information Processing Systems, 34:20243–20257, 2021. D. Applegate, M. Díaz, O. Hinder, H. Lu, M. Lubin, B. O’Donoghue, and W. Schudy. Pdlp: A practical first-order method for large-scale linear programming. arXiv preprint arXiv:2501.07018, 2025. N. Blin, S. Gualandi, C. Maes, A. Lodi, and B. Stellato. Batched first-order methods for parallel lp solving in mip. arXiv preprint arXiv:2601.21990, 2026. A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of mathematical imaging and vision, 40(1):120–145, 2011. G. B. Dantzig. Origins of the simplex method, page 141–151. Association for Computing Machinery, New York, NY, USA, 1990. ISBN 0201508141. URL https://doi.org/10.1145/87252. 88081. A. T. Ernst, H. Jiang, M. Krishnamoorthy, and D. Sier. Staff scheduling and rostering: A review of applications, methods and models. European journal of operational research, 153(1):3–27, 2004. A. Gleixner, G. Hendel, G. Gamrath, T. Achterberg, M. Bastubbe, T. Berthold, P. Christophel, K. Jarck, T. Koch, J. Linderoth, et al. Miplib 2017: data-driven compilation of the 6th mixedinteger programming library. Mathematical Programming Computation, 13(3):443–490, 2021. A. M. Gleixner, T. Berthold, B. Müller, and S. Weltge. Three enhancements for optimization-based bound tightening. Journal of Global Optimization, 67(4):731–757, 2017. T. Hübner, A. Gupte, and S. Rebennack. Spatial branch-and-bound for nonconvex separable piecewise linear optimization. INFORMS Journal on Computing, 38(2):645–675, 2026. B. Karimi, S. F. Ghomi, and J. Wilson. The capacitated lot sizing problem: a review of models and algorithms. Omega, 31(5):365–378, 2003. N. Karmarkar. A new polynomial-time algorithm for linear programming. In Proceedings of the sixteenth annual ACM symposium on Theory of computing, pages 302–311, 1984. A. B. Keha, I. R. de Farias Jr, and G. L. Nemhauser. A branch-and-cut algorithm without binary variables for nonconvex piecewise linear optimization. Operations research, 54(5):847–858, 2006. B. Knueven, J. Ostrowski, and J.-P. Watson. On mixed-integer programming formulations for the unit commitment problem. INFORMS Journal on Computing, 32(4):857–876, 2020. A. H. Land and A. G. Doig. An automatic method for solving discrete programming problems. In 50 Years of Integer Programming 1958-2008: From the Early Years to the State-of-the-Art, pages 105–132. Springer, 2009. M. S. Lobo, M. Fazel, and S. Boyd. Portfolio optimization with linear and fixed transaction costs. Annals of Operations Research, 152(1):341–365, 2007. H. Lu and J. Yang. cupdlp. jl: A gpu implementation of restarted primal-dual hybrid gradient for linear programming in julia. Operations Research, 73(6):3440–3452, 2025. 12

H. Lu, J. Yang, H. Hu, Q. Huangfu, J. Liu, T. Liu, Y. Ye, C. Zhang, and D. Ge. cupdlp-c: A strengthened implementation of cupdlp for linear programming by c language. arXiv preprint arXiv:2312.14832, 2023. H. Lu, Z. Peng, and J. Yang. cupdlpx: A further enhanced gpu-based first-order solver for linear programming. arXiv preprint arXiv:2507.14051, 2025. V. Nair, S. Bartunov, F. Gimeno, I. Von Glehn, P. Lichocki, I. Lobov, B. O’Donoghue, N. Sonnerat, C. Tjandraatmadja, P. Wang, et al. Solving mixed integer programs using neural networks. arXiv preprint arXiv:2012.13349, 2020. G. L. Nemhauser and L. A. Wolsey. Integer and combinatorial optimization, volume 18. Wiley New York, 1988. NVIDIA. Nvidia cuopt: GPU accelerated decision optimization, 2025. URL https://github. com/NVIDIA/cuopt. NVIDIA Corporation. cuSPARSE Library, 2024. URL https://developer.nvidia.com/ cusparse. Included in CUDA Toolkit. J. P. Pedroso, M. Kubo, and A. Viana. Unit commitment with valve-point loading effect. arXiv preprint arXiv:1404.4944, 2014. S. Wang, H. Zhang, K. Xu, X. Lin, S. Jana, C.-J. Hsieh, and J. Z. Kolter. Beta-crown: Efficient bound propagation with per-neuron split constraints for neural network robustness verification. Advances in neural information processing systems, 34:29909–29921, 2021. S. J. Wright. Primal-dual interior-point methods. SIAM, 1997. K. Xu, Z. Shi, H. Zhang, Y. Wang, K.-W. Chang, M. Huang, B. Kailkhura, X. Lin, and C.-J. Hsieh. Automatic perturbation analysis for scalable certified robustness and beyond. Advances in Neural Information Processing Systems, 33:1129–1141, 2020. K. Xu, H. Zhang, S. Wang, Y. Wang, S. Jana, X. Lin, and C.-J. Hsieh. Fast and complete: Enabling complete neural network verification with rapid and massively parallel incomplete verifiers. In International Conference on Learning Representations, 2021. URL https://openreview.net/ forum?id=nVZtXBI6LNn.

Broader Impacts

B3 -PWL is a foundational systems contribution: a GPU-accelerated framework for solving piecewiselinear mixed-integer optimization problems to global optimality. It is not tied to any specific deployed application, but PWL-MIP formulations underlie a wide range of decision-making tasks, including portfolio optimization, workforce scheduling, energy and resource allocation, and supplychain planning. Faster and more accessible exact solvers for this class of problems reduce the computational and energy cost of solving real-world planning problems, and lower the barrier for smaller organizations and academic users who lack access to commercial-grade CPU solvers. We do not foresee direct negative societal impacts from this work. The framework solves classical optimization problems whose objectives and constraints are specified by the user; it does not learn from data, generate content, or make autonomous decisions. As with any general-purpose optimization tool, downstream impact depends on the application, and we view such considerations as properties of mathematical optimization as a whole rather than risks specific to our contribution.

Limitations and Future Work

While B3 -PWL demonstrates that batched first-order LP solvers can serve as the central engine of branch-and-bound for PWL optimization, the current framework has several limitations that point to natural directions for future work. Slow dual-bound certification. First-order methods such as PDHG converge sublinearly in highaccuracy regimes, making dual-bound tightening substantially slower than with simplex-based solvers. Empirically (Section 5.4), B3 -PWL typically locates a high-quality incumbent within the first few GPU batches, but the remaining wall-clock budget is dominated by closing the final dual gap to certify optimality, and on a small number of instances (e.g., the largest resource-allocation cases in Appendix E) the 1% MIP gap is not closed within the 300 s budget despite a near-optimal incumbent being found early. Hybridizing first-order LP solves with simplex- or interior-point-based crossover for high-accuracy phases is a promising direction for narrowing this gap. No warm starting between subproblems. Simplex-based MIP solvers exploit the structural similarity between parent and child node LPs through dual-simplex warm starting from the parent basis, which often resolves a child node in only a handful of pivots. PDHG, in contrast, lacks a basis representation, and our batched first-order method solvers re-solve each subproblem essentially from scratch despite the high similarity between sibling LPs. Designing warm-start strategies for batched first-order B&B (e.g., transferring iterate scaling, primal weights, or partial restart history across related subproblems) could substantially reduce per-node work and is, to our knowledge, an open problem in the first-order MIP literature. Simple branching and no cutting planes. PWL optimization is a broad problem class with highly heterogeneous structure across applications. Our current pipeline focuses primarily on GPU-batched node-LP solving and SOS2-aware feasibility search, but it does not yet incorporate the full range of solver machinery used by highly engineered commercial MIP solvers, such as sophisticated branching rules, cut generation, bound tightening, and advanced presolve strategies. Future work could explore branching heuristics that are robust across diverse PWL structures, as well as GPU-compatible cutting-plane and bound-tightening routines that can be integrated into the batched first-order pipeline. By tightening LP relaxations and guiding the search toward more informative regions of the B&B tree, these extensions may substantially improve pruning efficiency, reduce dual-bound certification time, and narrow the gap between first-order GPU-centric B&B and mature commercial solvers.

cuPDLPx: GPU-Accelerated First-Order LP Solver

Our framework utilizes cuPDLPx [Lu et al., 2025] as the linear programming solver of subproblem relaxations in the B&B search tree, as it is the state-of-the-art, GPU-accelerated first-order solver designed to solve LP problems as an equivalent primal-dual saddle point problem. 14

Primal-Dual Hybrid Gradient (PDHG) The foundational algorithm of cuPDLPx is the PrimalDual Hybrid Gradient (PDHG). Given primal variables x, dual variables y, objective vector c, and constraint matrix A, PDHG iteratively updates the primal-dual iterate z k = [xk ; y k ] by alternating between a primal projection step and a dual update step: xk+1 = projX xk − τ (c − A⊤ y k ) , (9) y k+1 = y k − σA(2xk+1 − xk ) − σprojS σ −1 y k − A(2xk+1 − xk ) , (10) where X and S denote the primal and constraint feasible domains respectively, and τ, σ are the primal and dual step sizes. cuPDLPx further reparameterizes τ, σ with a base step size η and a primal weight ω (i.e., τ = η/ω, σ = ηω). Because operations in (9) and (10) are dominated by matrix-vector multiplications, PDHG are well suited for GPU acceleration. Algorithmic Enhancements in cuPDLPx To overcome the slow tail convergence of vanilla PDHG algorithm, cuPDLPx augments this base routine with four advanced mechanisms, achieving robust and accelerated convergence for large-scale LPs. Reflected Halpern Scheme Instead of standard average schemes, cuPDLPx employs a Reflected Halpern update. It interpolates between an over-relaxed PDHG step and an anchor point z 0 : k + 1 1 z0, (11) (1 + γ)PDHG(z k ) − γz k + z k+1 = k+2 k+2 where γ ∈ [0, 1] is the reflection parameter. This mechanism enables larger effective step sizes and strictly bounds the fixed-point residual. Adaptive Restart The anchor z 0 is periodically dynamically reset. cuPDLPx monitors the fixed-point error metric r(z) = ∥z − PDHG(z)∥P . When specific decay conditions are met (e.g., sufficient decay or lack of local progress), it restarts by resetting z 0 to the current iterate, keeping the search tightly focused on the optimal neighborhood. Constant Stepsize To maximize parallel efficiency on GPUs, cuPDLPx eschews the sequential adaptive stepsize heuristics typical of CPU solvers. Instead, it adopts a constant stepsize η ≈ 1/∥A∥2 , efficiently approximated via power iteration. PID-Controlled Primal Weight To balance the progress between primal and dual spaces, the primal weight ω is dynamically regulated. At each restart occurrence, cuPDLPx updates ω using a Proportional-Integral-Derivative (PID) controller based on the logarithmic gap between the primal and dual distances to optimality.

Numerical Robustness and Stopping-Criterion Sensitivity

We further examine the numerical robustness of B3 -PWL with respect to the PDLP feasibility/optimality tolerance and the global B&B stopping criterion. Specifically, we evaluate the feasibility and stability of the returned incumbents, the validity and stability of the final lower bounds, and sensitivity to tighter target MIP gaps. Incumbent feasibility and stability. We rerun all 43 instances with the PDLP feasibility and optimality tolerances tightened from 10−4 to 10−5 and 10−6 . Under all three tolerance settings, every instance returns an incumbent within 1% of the Gurobi reference objective. Relative to the 10−4 baseline, the maximum change in the solver-reported incumbent objective is only 0.0213% under the tighter 10−5 and 10−6 tolerances. We additionally validate each returned incumbent by fixing its integer/SOS2 pattern and re-solving the remaining continuous LP with HiGHS at a 10−9 feasibility tolerance. Across the three PDLP tolerance settings, this yields 43 × 3 = 129 incumbent patterns. All 129 patterns yield feasible solutions, with a maximum constraint violation of 2.91 × 10−11 . The maximum objective adjustment introduced by this validation is 0.0180% at 10−4 , 0.000457% at 10−5 , and 0.000933% at 10−6 . Lower-bound stability and validity. For each tighter PDLP tolerance, we compare the final lower bounds against the 10−4 baseline. We also verify ⋆ LB ≤ zGurobi ≤ U Bvalidated 15

Table 8: Validity and stability of final lower bounds under tighter PDLP feasibility/optimality tolerances. LB shifts are instance-wise relative changes from the 10−4 baseline. PDLP tolerance

Valid finite bounds

Median LB shift

95th-percentile LB shift

Maximum LB shift

43/43 43/43 43/43

Baseline 0.00193% 0.00329%

Baseline 0.167% 0.371%

Baseline 0.206% 0.989%

−4

10 10−5 10−6

Table 9: Sensitivity to the global B&B target MIP gap. Incumbent and final-LB shifts are measured relative to the 1% configuration. Target MIP gap

Eff-Opt

Maximum incumbent shift

Maximum final-LB shift

Valid finite bounds

Time-limit runs

1% 0.5% 0.1%

43/43 43/43 43/43

Baseline 0.000217% 0.000217%

Baseline 0.0989% 0.102%

43/43 43/43 43/43

3 7 12

for every available finite lower bound, where U Bvalidated denotes the objective of the independently validated incumbent. Table 8 reports instance-wise relative changes in the final lower bounds with respect to the 10−4 baseline. The final lower bounds remain valid for all 43 instances under each tolerance setting and exhibit only small changes relative to the 10−4 baseline. Sensitivity to the B&B stopping criterion. We further tighten the target relative MIP gap from 1% to 0.5% and 0.1%, while keeping the PDLP tolerance and 300-second wall-clock limit fixed. Table 9 reports incumbent and final-LB shifts relative to the 1% baseline. Tightening the MIP-gap criterion does not materially change the returned incumbents or final lower bounds. Its primary effect is to increase the cost of dual-bound certification: the number of runs reaching the 300-second limit increases from 3 at a 1% target gap to 12 at a 0.1% target gap. Thus, the main solution-quality conclusions remain stable under stricter stopping criteria, while closing a tighter global optimality gap requires substantially more computation.

Per-Instance Detailed Results

Tables 10–12 give per-instance wall-clock time and returned objective for PWL-BaB and cuOpt on every canonical instance, alongside the Gurobi certified optimum (MIPG AP=10−6 , separate 32-thread Xeon CPU run; Gurobi runtime is omitted because it is not directly comparable to the GPU-accelerated PWL/cuOpt times). TO marks the 300 s budget exhaustion; — marks no feasible incumbent returned. Each instance is tagged pwl{N}_s{seed}, where N is the number of PWL functions in the instance and seed is the random-generator seed. Within resource allocation, two pairs of size groups share the same N (pwl300 and pwl500); we disambiguate them with the suffixes a/b, whose generator parameters are listed in Appendix G.

Implementation Details

Multi-level branching for batch filling. GPU throughput is maximized when all B slots are occupied, but the B&B queue may contain fewer than B nodes, particularly in early iterations. Inspired by the split-depth mechanism in αβ-CROWN [Xu et al., 2021], we preemptively branch nodes multiple levels before solving. Specifically, when the queue size |Q| < 0.9B, we compute a split depth 0.9B dsplit = log2 , (12) max(1, |Q|) and expand each seed node into up to 2dsplit leaf nodes via BFS-style recursive branching, using the branching heuristic to select variables at each level. This ensures the GPU batch is consistently saturated. 16

Table 10: Per-instance results: Portfolio (18 instances). PWL-BaB Instance

cuOpt

Gurobi

Time (s)

Objective

Time (s)

Objective

0.6 0.6 0.5 1.2 1.7 1.2 1.4 1.3 1.3 2.6 3.0 2.1 5.8 2.9 3.2 3.6 3.4 8.6

-1636.0085 -1591.3689 -1702.2300 -3329.2898 -3493.2999 -2820.2763 -5132.8143 -5312.8877 -4326.8367 -6729.8250 -6930.5113 -6186.0892 -8793.2572 -8844.0370 -8396.4138 -9917.6555 -10257.0782 -10361.3325

TO TO TO TO TO TO 7.8 8.8 10.1 12.9 12.0 13.8 17.8 18.7 TO TO 27.1 30.3

— — — -3103.7171 — -2752.9085 -5092.7688 -5270.0022 -4314.2420 -6717.2791 -6908.1034 -6150.4538 -8752.8616 -8822.4380 — — -10242.9353 -10332.3481

-1636.0084 -1591.3686 -1702.2300 -3329.2898 -3493.2999 -2820.2763 -5132.8144 -5313.4625 -4326.8434 -6729.8211 -6930.5117 -6189.1684 -8793.2572 -8844.0369 -8398.0055 -9918.4741 -10257.0767 -10361.3306

pwl100_s42 pwl100_s123 pwl100_s999 pwl200_s42 pwl200_s123 pwl200_s999 pwl300_s42 pwl300_s123 pwl300_s999 pwl400_s42 pwl400_s123 pwl400_s999 pwl500_s42 pwl500_s123 pwl500_s999 pwl600_s42 pwl600_s123 pwl600_s999

Table 11: Per-instance results: Workforce (9 instances). PWL-BaB Instance pwl300_s42 pwl300_s123 pwl300_s999 pwl600_s42 pwl600_s123 pwl600_s999 pwl1200_s42 pwl1200_s123 pwl1200_s999

cuOpt

Gurobi

Time (s)

Objective

Time (s)

Objective

0.8 0.8 0.8 2.1 1.8 2.1 4.8 5.5 4.2

281.7799 279.7645 277.9249 563.5598 559.5290 555.8498 1121.4303 1115.9629 1121.4227

12.5 11.6 12.2 27.9 14.8 35.8 107 TO 98.3

283.5340 279.7645 279.3438 567.8821 563.2236 559.9900 1127.8207 — 1125.9802

281.7799 279.7645 277.9249 563.5598 559.5290 555.8498 1121.4303 1115.9629 1121.4227

Adaptive batch sizing. We dynamically adjust B based on GPU memory utilization. After each solving round, we query the available GPU memory and compute a capacity factor γ = Mtarget /Mused , where Mtarget = 0.85 · Mtotal . The batch size is updated as B ← min(2B, ⌊γB⌋), subject to a minimum increase threshold, ensuring stable growth without out-of-memory failures.

Benchmark Details

To our knowledge, there is no broadly adopted, MIPLIB-like benchmark suite specifically for largescale nonconvex PWL optimization that preserves native SOS2 structure. In particular, the MIPLIB 2017 paper reports that, although its input format permits SOS and piecewise-linear constraints, none of the submitted instances used these structures. Existing PWL studies therefore often rely on generated application-oriented families. Recent external studies have released reproducible generators and instances for nonconvex PWL network-flow and knapsack problems, but these remain generated benchmarks rather than a broadly adopted standard [Hübner et al., 2026]. Public application-specific libraries also exist, notably PGLib-UC, which combines convex piecewise-linear production costs with ordinary binary commitment, startup, and shutdown decisions, and therefore targets a different problem scope [Knueven et al., 2020]. Given this landscape, we construct a controlled systems benchmark to characterize the behavior of our GPU B&B framework, rather than to claim broad industrial representativeness. We selected three classical PWL application classes, portfolio optimization, workforce scheduling, and resource allocation, whose generators allow us to systematically vary LP size, breakpoint count, sparsity, and 17

Table 12: Per-instance results: Resource Allocation (16 instances). PWL-BaB Instance pwl300a_s42 pwl300a_s123 pwl300a_s999 pwl300b_s123 pwl400_s42 pwl400_s123 pwl400_s999 pwl500a_s42 pwl500a_s123 pwl500a_s999 pwl500b_s42 pwl500b_s123 pwl500b_s999 pwl1200_s42 pwl1200_s123 pwl1200_s999

cuOpt

Gurobi

Time (s)

Objective

Time (s)

Objective

64.1 64.5 64.0 43.4 9.0 8.9 8.9 40.0 40.0 39.8 11.3 11.4 11.4 TO TO TO

316.2431 316.2431 316.2431 316.2431 421.6067 421.6067 421.6067 527.0718 527.0718 527.0718 527.0083 527.0083 527.0083 1278.4221 1280.2254 1287.2700

28.5 38.4 25.4 23.2 31.0 TO 31.4 69.0 37.4 81.0 48.2 60.7 52.1 TO TO TO

316.2431 316.2431 316.2431 316.2431 421.6067 — 425.3778 527.0718 527.0718 527.7716 527.0083 531.4298 531.7667 — — —

316.2431 316.2431 316.2431 316.2431 421.6067 421.6067 421.6067 527.0718 527.0718 527.0718 527.0083 527.0083 527.0083 1278.1789 1280.0844 1287.0155

SOS2 density and thereby characterize the behavior of our GPU B&B framework across different scales and structures. We remove root-trivial instances, since they primarily benchmark a single LP relaxation rather than the batched B&B mechanism studied here. Synthetic generation also allows us to control instance difficulty within an informative range, avoiding cases that are either too easy or too difficult to meaningfully distinguish our framework from CPU simplex-based solvers. Our benchmark suite spans three problem families that exercise complementary structural regimes for branch-and-bound with batched LP relaxations. The portfolio benchmark (Section G.1) features dense SOS2 coverage—one SOS2 substructure per asset—coupled through a small number of global linear constraints. The workforce scheduling benchmark (Section G.2) features a multi-period grid of SOS2 substructures coupled through skill-coverage and budget constraints. The resource allocation benchmark (Section G.3) concentrates combinatorial structure on a sparse subset of variables under a block-angular linear coupling. All three families share the same univariate PWL cost-curve generator (described in Section G.1) so that differences in solver behavior across families reflect coupling structure rather than curve shape. G.1

Portfolio Benchmark with Piecewise-Linear Transaction Costs

Problem Formulation. We consider a portfolio optimization benchmark with separable PWL transaction-cost penalties. This problem models capital allocation across N assets under a global budget, sparse factor-risk limits, and sector exposure caps. The objective balances expected return against a per-asset transaction penalty. For each asset j ∈ {1, . . . , N }, let wj ≥ 0 denote the dollar holding and let zj denote the transactioncost penalty associated with wj . The objective is min − w,z

N X

µj wj + γ

j=1

N X

zj ,

(13)

j=1

where µj is the expected return of asset j, and γ > 0 controls the trade-off between return and transaction cost. The portfolio is subject to three classes of linear coupling constraints. First, the total capital invested cannot exceed the budget C: N X

wj ≤ C.

j=1

(14)

Second, for each risk factor i ∈ {1, . . . , F }, the factor exposure is bounded by a risk budget Birisk : N X

σi,j wj ≤ Birisk ,

(15)

j=1

where σi,j is the loading of asset j on factor i. Third, sector allocations are capped. Let s(j) denote the sector of asset j; then for each sector s ∈ {1, . . . , Ssec }, X wj ≤ Bssec . (16) j: s(j)=s

Transaction cost is modeled independently for each asset by a univariate PWL function. For asset j, let {(dj,k , ϕj,k )}Jk=1 be the breakpoints of its cost curve, where dj,k is a breakpoint on the holding axis and ϕj,k is the corresponding penalty value. We use the standard SOS2 formulation with convex-combination variables λj,k ∈ [0, 1]: J X

λj,k = 1,

(17)

k=1

wj =

J X

dj,k λj,k ,

(18)

ϕj,k λj,k ,

(19)

k=1

zj =

J X k=1

together with an SOS2 constraint on λj,· , so that at most two adjacent breakpoints are active. This enforces that (wj , zj ) lies on the graph of the PWL penalty function. The complete formulation is min

w,z,λ

s.t.

−

N X

µj wj + γ

j=1

j=1 N X

N X

wj ≤ C,

j=1 N X

σi,j wj ≤ Birisk ,

i = 1, . . . , F,

j=1

wj ≤ Bssec ,

s = 1, . . . , Ssec ,

j: s(j)=s J X

(20)

λj,k = 1,

j = 1, . . . , N,

k=1

wj =

J X

dj,k λj,k ,

j = 1, . . . , N,

ϕj,k λj,k ,

j = 1, . . . , N,

k=1

zj =

J X k=1

λj,· is SOS2, wj ∈ [wj , wj ], Data Generation. point mj :

j = 1, . . . , N, j = 1, . . . , N.

The holding interval for each asset is restricted to a tight box around a reference [wj , wj ] = [mj − 0.35 spanj , mj + 0.35 spanj ]. 19

(21)

This narrow box is chosen deliberately so that no single SOS2 segment is admissible at the LP relaxation, ensuring that branch-and-bound must explore multiple segment combinations rather than terminating at the root LP. The breakpoints (dj,k , ϕj,k ) are sampled from a univariate U-shaped curve generator that combines a low-volume penalty region, a flat economical mid-range, and a steeply rising tail; the same generator is reused across all three benchmark families to isolate the effect of coupling structure. Expected returns are sampled independently as µj ∼ U(0.02, 0.15), and the global transaction-cost weight is fixed to γ = 0.005. Asset-wise scales are drawn as wjmax ∼ U(50, 500). Assets are assigned to sectors in round-robin order, s(j) = j mod Ssec . For each factor i, the loading matrix σi,j has approximately 30% of entries nonzero (each entry independently masked with probability 0.7), with nonzero values drawn from U(0.01, 0.10). Let wjmax,band = wj denote the upper end of the tight box. The total budget is set to C = αC

wjmax,band ,

αC ∼ U(0.6, 0.8).

Risk budgets are set by Birisk = max αi σi⊤ wmax,band , 1 ,

αi ∼ U(0.8, 1.2),

and sector budgets by Bssec = αs

wjmax,band ,

αs ∼ U(0.75, 0.95).

j: s(j)=s

Problem Size.

For N assets and J breakpoints per asset, the model contains

N holding variables + N cost variables + N J SOS2 weights = N (J + 2) variables in total. The number of global coupling constraints is 1 + F + Ssec , while the PWL representation contributes 3N local constraints. This benchmark is therefore characterized by many homogeneous SOS2 substructures coupled through a relatively small set of sparse linear constraints.

G.2

Workforce Scheduling with Overtime Costs

Problem Formulation. We consider a multi-period workforce scheduling problem with skill coverage and budget constraints. Let w = 1, . . . , W index workers, p = 1, . . . , P index periods, and t = 1, . . . , T index skills. The decision variable hw,p denotes the number of hours assigned to worker w in period p, and zw,p denotes the corresponding labor cost. Each worker has a skill coefficient κw,t ∈ [0, 1] for skill t, and each period p requires at least mt,p effective hours of each skill t. In addition, each period is subject to a payroll budget Bp . To model overtime and other nonlinear labor costs, each worker–period pair (w, p) is associated with a piecewise-linear cost function zw,p = fw,p (hw,p ). We represent this function using J breakpoints {(dw,p,k , ϕw,p,k )}Jk=1 and SOS2 variables λw,p,k , yielding J X k=1

λw,p,k = 1,

hw,p =

J X

dw,p,k λw,p,k ,

k=1

zw,p =

J X k=1

together with an SOS2 constraint on λw,p,· . 20

ϕw,p,k λw,p,k ,

The full optimization problem is min

h,z,λ

s.t.

W X P X

zw,p

w=1 p=1 W X

κw,t hw,p ≥ mt,p ,

∀t, p,

rw hw,p ≤ Bp ,

∀p,

λw,p,k = 1,

∀w, p,

w=1 W X w=1 J X k=1

hw,p =

J X

dw,p,k λw,p,k ,

∀w, p,

ϕw,p,k λw,p,k ,

∀w, p,

k=1

zw,p =

J X k=1

hw,p ∈ [hw , hw ], λw,p,k ∈ [0, 1], SOS2(λw,p,1 , . . . , λw,p,J ),

∀w, p, k, ∀w, p.

Data Generation. Worker-specific maximum hours are sampled as hmax ∼ U(8, 16), wage rates w as rw ∼ U(15, 50), and staffing requirements as mt,p ∼ U(2, 8). The skill matrix κw,t has approximately 40% of entries nonzero, with nonzero entries sampled from U(0.3, 1.0); if a worker draws no nonzero skill, one randomly chosen entry is resampled from U(0.3, 1.0) to ensure every worker is qualified for at least one skill. Period budgets are set to 0.7–0.9 times a reference full-capacity payroll level. Each worker–period pair receives an independently sampled PWL cost curve from the same generator used in Section G.1, with curves associated with the same worker sharing a common hour range determined by that worker’s feasible shift length. Problem Size.

For W workers, P periods, and J breakpoints per cost curve, the model contains

W P hour variables + W P cost variables + W P J SOS2 weights = W P (J + 2) variables in total, plus T P skill-coverage constraints, P payroll-budget constraints, and 3W P local linear equations from the PWL encoding. In our benchmark suite, the number of local PWL functions (W P ) typically ranges from 200 to 2000. G.3

Block-Angular Resource Allocation with Piecewise-Linear Costs

Problem Formulation. We consider a multi-line production-planning benchmark with blockangular resource structure and separable PWL operating-cost penalties on a subset of products. This problem models a manufacturer operating G product groups (e.g., plants or product lines), in which each group draws on its own dedicated equipment (local resources) but also competes for a shared pool of global resources (e.g., raw-material budget, warehouse capacity, logistics throughput). A subset of products carries a non-convex production-cost curve—capturing setup amortization at low volume, economies of scale at medium volume, and overtime or congestion at high volume—which we model as a univariate PWL penalty. For each group g ∈ {1, . . . , G}, let Jg denote the set of products in that group, with |Jg | = K and N = GK. For each product j ∈ {1, . . . , N } the decision variable xj ≥ 0 is the production level. A subset P ⊆ {1, . . . , N } of size Npwl carries a PWL cost variable zj , and the objective minimizes the total PWL operating cost, X min zj . (22) x,z

j∈P

The constraint matrix exhibits a block-angular structure: G diagonal local blocks coupled through one shared global block. For each group g, L local capacity constraints couple only products in Jg : X loc aloc r = 1, . . . , L, (23) g,r,j xj ≤ cg,r , j∈Jg

while M global capacity constraints couple all N products: N X glob aglob m,j xj ≤ cm ,

m = 1, . . . , M.

(24)

j=1

All linear coupling constraints are inequalities. For each j ∈ P, the PWL operating cost is parameterized by J breakpoints {(dj,k , ϕj,k )}Jk=1 , encoded with SOS2 convex-combination variables λj,k ∈ [0, 1]: J J J X X X λj,k = 1, xj = dj,k λj,k , zj = ϕj,k λj,k , (25) k=1

k=1

together with an SOS2 constraint on λj,· . The complete formulation is X min zj x,z,λ

s.t.

j∈P

loc aloc g,r,j xj ≤ cg,r ,

g = 1, . . . , G, r = 1, . . . , L,

j∈Jg N X glob aglob m,j xj ≤ cm ,

m = 1, . . . , M,

j=1 J X

λj,k = 1,

j ∈ P,

k=1

xj =

J X

dj,k λj,k ,

j ∈ P,

ϕj,k λj,k ,

j ∈ P,

k=1

zj =

J X k=1

SOS2(λj,1 , . . . , λj,J ), j ∈ P, xj ∈ [xj , xj ], j = 1, . . . , N, where [xj , xj ] is a tight box for j ∈ P (specified below) and the wider default range [0, 100] for j∈ / P. Data Generation. For each group g, the local block Aloc ∈ RL×K is generated with target g sparsity ρloc = 0.02; nonzero positions are drawn uniformly at random and nonzero values are sampled independently from U(−2, 2). The global block Aglob ∈ RM ×N is generated with sparsity ρglob = 0.008 and the same value distribution; allowing both signs in Aglob models global rows that aggregate signed contributions across products (e.g., net demand vs. supply of a shared resource). glob Capacities are sampled as cloc ∼ U(1, 50). g,r , cm The PWL-active subset P is drawn uniformly at random without replacement from {1, . . . , N }; breakpoint values (dj,k , ϕj,k ) are obtained from the same univariate U-shaped curve generator described in Section G.1, applied over the nominal production range [xmin , xmax ] = [0.01, 20]. To force a non-trivial branch-and-bound search, every j ∈ P is restricted to a tight box of width 0.7(xmax − xmin ) around the midpoint of the PWL range, [xj , xj ] = [m̄ − 0.35 span, m̄ + 0.35 span], 22

with m̄ = (xmin +xmax )/2 and span = xmax −xmin ; products outside P take the wider default range xj ∈ [0, 100]. After the boxes are imposed, the right-hand sides of both blocks are post-processed by adding a uniform slack s ∼ U(5, 30) on top of the worst-case left-hand side minx∈[x,x] Ai,: x for each row i, so that the linear feasibility region {x : Ax ≤ c, x ≤ x ≤ x} is non-empty by construction. Problem Size. For N = GK products, an active subset of size Npwl , and J breakpoints per PWL function, the model has N production variables + Npwl cost variables + Npwl J SOS2 weights = N + Npwl (J + 1) variables in total. The number of linear coupling constraints is GL + M (block-diagonal local plus dense global), while the PWL representation contributes 3Npwl local linear equations. In contrast to the portfolio and workforce benchmarks—where every primary variable carries an SOS2 substructure—this benchmark concentrates combinatorial structure on only a fraction of the variable set while keeping a sparse, hierarchical linear coupling.

Problem Characteristics and Solver Behavior

We analyze how structural differences across the three problem families affect the behavior of our GPU B&B framework. PWL coverage and sparsity. Resource-allocation instances embed a relatively small fraction of PWL variables in a much larger sparse block-angular LP. Because branching on one PWL variable modifies only a local part of this structure, its improvement to the global lower bound can be small relative to the cost of solving the full node LP. Consequently, more branching may be required before effective pruning occurs. Meanwhile, this sparse structure favors mature simplex-based solvers through basis reuse, sparse factorization, cuts, and other B&B enhancements. LP size. The node LPs in the resource-allocation family are also substantially larger, reducing the maximum batch size that fits within the same GPU memory and increasing the cost of each batched iteration. Search-tree structure. These properties tend to produce narrower, deeper search trees because each local PWL branch provides limited global tightening. In contrast, portfolio and workforce expose many independent PWL decisions and tend to generate shallower, wider frontiers that better match GPU batching. Branching quality is consequently more important for resource allocation: an unproductive early decision can send the search deep into a poor direction and require many weak or infeasible node LPs to be solved before pruning. Unlike simplex methods, first-order methods do not return basic vertex solutions and converge slowly at high dual accuracy, weakening lower-bound-based pruning, especially in deeper search trees.

B$^3$-PWL: GPU-Batched Branch-and-Bound for Piecewise-Linear Optimization with SOS2 Constraints

Analyze & explore related research in Conceptio

Related documents