cuPDLPx: GPU-Accelerated First-Order LP Solver

• cuPDLPx is a GPU-based linear programming solver that employs the primal–dual hybrid gradient method to map LP solutions primarily to sparse matrix–vector products.
• It integrates adaptive restart policies, discrete PID control for primal weight updates, and diagonal scaling to balance progress and enhance convergence.
• Its design leverages CUDA optimizations, kernel fusion, and scalable multi-GPU strategies to achieve significant speedups in benchmark LP problems.

Searching arXiv for papers on cuPDLPx and related work. cuPDLPx is a GPU-based first-order solver for linear programming built on the primal–dual hybrid gradient (PDHG) method. Across its published variants, it is described as a further enhanced GPU-based solver for LP, a C/CUDA implementation of restarted PDHG, and the computational substrate underlying later developments such as distributed multi-GPU LP solving, lightweight presolve pipelines, and hybrid MIP heuristics. Its central design premise is to map LP solution almost entirely to sparse matrix–vector products, elementwise projections, restart logic, and reductions that are well matched to modern GPU architectures (Lu et al., 18 Jul 2025, Lu et al., 2 Jun 2025).

1. Problem formulation and core iteration

cuPDLPx is formulated for general LPs in bound-constrained or extended conic form. One common presentation is

$$\min c^T x \quad \text{subject to} \quad A x = b,\; Gx \ge h,\; \ell \le x \le u,$$

with saddle-point formulation

$$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$

where $K^T=[G^T\,|\,A^T]$, $q^T=[h^T\,|\,b^T]$, $X=\{x:\ell\le x\le u\}$, and $Y=\{y:y_{1:m_1}\ge 0\}$. Another equivalent presentation uses row-wise bounds $\underline b \le Ax \le \overline b$ and variable-wise bounds $\underline x \le x \le \overline x$ (Lu et al., 2 Jun 2025, Cederberg et al., 27 Apr 2026).

The core PDHG iteration in cuPDLPx is

$$x^{k+1}=\operatorname{proj}_X\!\Bigl(x^k-\tau(c-K^T y^k)\Bigr),$$

$$y^{k+1}=\operatorname{proj}_Y\!\Bigl(y^k+\sigma\bigl(q-K(2x^{k+1}-x^k)\bigr)\Bigr),$$

with step sizes parameterized as

$$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$0

where $$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$1 is the overall step size and $$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$2 is the primal weight. A safe choice stated in the survey literature is $$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$3 with $$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$4 (Lu et al., 2 Jun 2025).

A complementary formulation used in the distributed implementation writes the LP in saddle form

$$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$5

with

$$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$6

$$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$7

$$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$8

where $$\min_{x\in X}\max_{y\in Y} L(x,y)=c^T x-y^T Kx+q^T y,$$9 and $K^T=[G^T\,|\,A^T]$0. Convergence is monitored via three relative KKT errors $K^T=[G^T\,|\,A^T]$1, and the solver stops when $K^T=[G^T\,|\,A^T]$2 (Li et al., 12 Jan 2026).

2. Restarted Halpern PDHG, weight control, and conditioning

Later cuPDLPx papers place restarted Halpern PDHG at the center of the algorithm. The Halpern interpolation step is written as

$K^T=[G^T\,|\,A^T]$3

and with reflection parameter $K^T=[G^T\,|\,A^T]$4,

$K^T=[G^T\,|\,A^T]$5

This is the restarted Halpern PDHG core described for cuPDLPx (Lu et al., 18 Jul 2025).

A defining enhancement is a restart criterion based on the fixed-point residual

$K^T=[G^T\,|\,A^T]$6

The implementation evaluates three restart triggers: sufficient decay, necessary decay combined with no local progress, and an artificial trigger based on maximum inner-loop length. The paper states that this residual-based rule is lighter to compute on the GPU than gap-based and KKT-based restarts, and is better aligned with the theoretical Halpern rates (Lu et al., 18 Jul 2025).

cuPDLPx also updates the primal weight $K^T=[G^T\,|\,A^T]$7 via a discrete PID control loop on $K^T=[G^T\,|\,A^T]$8. The tracking error is

$K^T=[G^T\,|\,A^T]$9

followed by

$q^T=[h^T\,|\,b^T]$0

The intended effect is to rebalance primal and dual progress at restarts (Lu et al., 18 Jul 2025).

The survey literature additionally describes diagonal scaling

$q^T=[h^T\,|\,b^T]$1

with $q^T=[h^T\,|\,b^T]$2 chosen by a combination of Ruiz equilibration and the diagonal preconditioner of Pock–Chambolle, and an adaptive line-search rule for $q^T=[h^T\,|\,b^T]$3. This suggests that cuPDLPx should be understood not as a bare PDHG implementation, but as a restart- and conditioning-intensive solver family whose practical behavior depends on scaling, restart policy, and primal-weight control (Lu et al., 2 Jun 2025).

3. GPU implementation and systems design

The single-GPU architecture of cuPDLPx is described as a pure-CUDA implementation of PDHG enriched with Halpern extrapolation, adaptive restarts, and PID-controlled primal weight, while the earlier C-language line is described as a ground-up rewrite in C/C++ with generic CUDA APIs (Li et al., 12 Jan 2026, Lu et al., 2023).

Its principal data structures are sparse matrix storage for $q^T=[h^T\,|\,b^T]$4 in CSR format and dense vectors for primal and dual variables, objective coefficients, bounds, and auxiliary state. In the distributed description, $q^T=[h^T\,|\,b^T]$5 is stored as rowPtr[m+1], colIdx[nnz], and val[nnz] in global GPU memory, while dense vectors are allocated via cudaMalloc and host-to-device transfers use pinned DMA (Li et al., 12 Jan 2026).

The kernel inventory is consistent across descriptions: sparse matrix–vector multiply for $q^T=[h^T\,|\,b^T]$6 and $q^T=[h^T\,|\,b^T]$7, fused vector updates such as axpy or scale-add, coordinatewise projection by clamping to bounds, Halpern extrapolation, and reductions for norms or residuals. The implementation literature emphasizes kernel fusion, persistent kernels for small-to-medium problems, warp-level and block-level reductions, CUDA streams for overlap, in-place updates, and reuse of auxiliary buffers (Lu et al., 18 Jul 2025, Li et al., 12 Jan 2026).

The C/CUDA reimplementation reports several systems-level changes relative to the Julia prototype: removal of dynamic Julia dispatch and garbage-collection overhead, reorganization of data into tightly packed C arrays for maximum memory throughput, fusion of vector updates and projections into custom kernels, tailor-tuned sparse-matrix-vector routines rather than high-level wrappers, and more efficient incorporation of standard presolve/scaling modules from CLP/COPT. On the same GPU hardware, this version achieves roughly 50% faster iteration rates (Lu et al., 2023).

A recurring engineering theme is that the algorithmic bottlenecks are intentionally reduced to GPU-favorable primitives: sparse MV, branch-free clamp, fused saxpy-style updates, and small reductions. This architecture is what makes cuPDLPx a representative example of GPU-native first-order LP solving rather than a GPU port of a factorization-based simplex or interior-point codebase.

4. Relationship to HPR-LP and theoretical interpretation

A later analytical paper argues that the base algorithm of cuPDLPx is a special case of the base algorithm of HPR-LP, and therefore that cuPDLPx is another concrete implementation instance of HPR-LP. In that formulation, with $q^T=[h^T\,|\,b^T]$8, one chooses

$q^T=[h^T\,|\,b^T]$9

and identifies the relevant iterates so that the HPR-LP and cuPDLPx sequences coincide exactly (Chen et al., 28 Sep 2025).

The same analysis states that, once the active sets have been identified, HPR-LP and EPR-LP become equivalent under the same initialization, and that all non-ergodic $X=\{x:\ell\le x\le u\}$0 convergence guarantees of HPR-LP apply directly to cuPDLPx’s base algorithm. In particular, the paper presents a KKT-residual estimate of the form

$X=\{x:\ell\le x\le u\}$1

This situates cuPDLPx within a broader operator-splitting lineage rather than treating it as an isolated solver design (Chen et al., 28 Sep 2025).

At the same time, the cuPDLPx literature emphasizes enhancements absent from the minimal base-algorithm comparison, notably restart policies, PID-controlled weight updates, and GPU-tailored execution. A plausible implication is that there are two distinct levels at which cuPDLPx is discussed in the literature: an algebraic level, where it is related to HPR-LP, and an implementation level, where restart heuristics, weight control, preconditioning, and fused kernels materially shape runtime behavior.

The same analytical paper also notes limitations of the base form: fixed step size $X=\{x:\ell\le x\le u\}$2, no adaptive penalty updates or restarts, possible stalling on poorly scaled problems, and open questions around degenerate steps $X=\{x:\ell\le x\le u\}$3. These caveats help explain why the enhanced cuPDLPx line places significant emphasis on restart and control mechanisms (Chen et al., 28 Sep 2025).

5. Benchmark performance and multi-GPU scaling

The enhanced single-GPU solver reports substantial gains over earlier PDLP implementations. Across benchmark datasets, cuPDLPx achieves 2.5x-5x speedups on MIPLIB LP relaxations and 3x-6.8x on Mittelmann’s benchmark set, with particularly strong improvements in high-accuracy and presolve-enabled settings (Lu et al., 18 Jul 2025).

Solver	MIPLIB $X=\{x:\ell\le x\le u\}$4 SGM$X=\{x:\ell\le x\le u\}$5	Mittelmann $X=\{x:\ell\le x\le u\}$6 SGM$X=\{x:\ell\le x\le u\}$7
cuPDLP (Julia)	11.07 s	67.38 s
HPR-LP (Julia)	6.50 s	24.28 s
cuPDLPx (Julia)	5.24 s	22.78 s
cuPDLPx (C)	4.39 s	19.11 s
Gurobi-PDLP	7.14 s	32.47 s

On the 94 “hard” MIPLIB instances, the high-accuracy, presolve-on speedup reaches nearly 5x, and on Mittelmann the reported speedups are up to 6.8x with presolve at $X=\{x:\ell\le x\le u\}$8, while remaining approximately 3x-5x at $X=\{x:\ell\le x\le u\}$9 (Lu et al., 18 Jul 2025).

The multi-GPU extension is built upon cuPDLPx and distributes both memory and computation using a two-dimensional grid partitioning of the constraint matrix. The reported hardware scale is up to 8 NVIDIA H100 GPUs with 80 GB HBM3 each, interconnected by NVLink. On strong-scaling experiments, MIPLIB Large instances with more than 10 million nonzeros improve from 115.7 s on 1 GPU to 66.0 s on 8 GPUs, sdm_50k_500k improves from 377 s on 1 GPU to 92 s on 4 GPUs and 61 s on 8 GPUs, and mcf_2500_100_500 improves from 2943 s on 1 GPU to 935 s on 4 GPUs and 504 s on 8 GPUs. The implementation preserves full FP64 numerical accuracy, and the empirical communication model indicates that communication remains below 10–15% of total runtime at 8 GPUs for well-chosen partitions (Li et al., 12 Jan 2026).

These results place cuPDLPx in two performance regimes. On a single GPU it is a high-performance first-order LP solver competitive with other GPU-accelerated PDLP variants; in distributed form it becomes a memory-scaling and throughput-scaling mechanism for instances that exceed single-GPU capacity.

6. Presolve integration, use in CHAP, and open directions

A major systems development around cuPDLPx is the introduction of the lightweight LP presolver PSLP. PSLP is an open-source C implementation adopted by the community, with integrations in cuPDLPx, cuOpt, and HPR-LP. It implements fast and medium reductions such as singleton rows, doubleton rows, redundant constraints, column singletons, variable locks, parallel rows, parallel columns, primal propagation, and dual propagation, while omitting slow explorers such as rank-revealing row dependence and symmetry detection (Cederberg et al., 27 Apr 2026).

On presolve time alone, PSLP is reported as 6–12x faster than Gurobi’s built-in presolver while capturing about 90%–94% of the reduction in nonzeros. In end-to-end cuPDLPx experiments, the aggregate MIPLIB-383 result changes from 371/383 solved with SGM10 10.77 for cuPDLPx alone to 381/383 with SGM10 6.13 for cuPDLPx + PSLP; on all 48 Mittelmann instances, the result changes from 42/48 with SGM10 69.01 to 46/48 with SGM10 39.57. The paper also notes that presolve can consume up to 20–40% of reduced-problem solve time, yet overall total time is typically cut by factors of 2–9 on large problems (Cederberg et al., 27 Apr 2026).

cuPDLPx has also been used as an approximate LP component inside CHAP, a hybrid GPU-CPU primal heuristic framework for mixed-integer programming. In CHAP, cuPDLPx runs PDHG for $Y=\{y:y_{1:m_1}\ge 0\}$0, $Y=\{y:y_{1:m_1}\ge 0\}$1, $Y=\{y:y_{1:m_1}\ge 0\}$2, and $Y=\{y:y_{1:m_1}\ge 0\}$3 iterations and, after each block, streams the current primal $Y=\{y:y_{1:m_1}\ge 0\}$4, dual $Y=\{y:y_{1:m_1}\ge 0\}$5, and reduced costs $Y=\{y:y_{1:m_1}\ge 0\}$6 into a shared solution pool under an “LP solutions” slot. The CPU and GPU heuristics then consume these snapshots without waiting for full LP convergence (Tjusila et al., 6 May 2026).

The CHAP ablation results show a nuanced but important effect. Replacing CPU PDHG with cuPDLPx changes the “Found” count from 47 to 46, but increases “Wins” from 3 to 13, reduces the shifted geometric mean final gap from 19.57 to 12.17, and reduces the primal integral from 72.59 to 51.01. Adding GPU tabu search gives 47 found, 22 wins, 8.84 gap, and 41.84 primal integral. The CHAP paper explicitly characterizes cuPDLPx in this setting as only an approximate LP solver and notes that first-order methods can converge slowly; it identifies checkpoint tuning, step-size scheduling, and tighter embedding with higher-precision crossover or dual simplex refinements on the CPU as future work (Tjusila et al., 6 May 2026).

Taken together, these developments show cuPDLPx as more than a standalone LP code. It is a reusable GPU-native first-order LP engine that has been successively expanded along three axes: stronger single-GPU algorithms, distributed multi-GPU execution, and integration into broader optimization pipelines such as presolve-first solve stacks and CPU-GPU MIP heuristics.