---
title: 'cuPDLPx: GPU-Accelerated First-Order LP Solver'
url: https://www.emergentmind.com/topics/cupdlpx
type: topic
---

# cuPDLPx: GPU-Accelerated First-Order LP Solver

Searching arXiv for recent 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 [2507.14051][2506.02174].

## 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\) [2506.02174][2604.23951].

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
\[
\tau=\frac{\eta}{\omega},\qquad \sigma=\eta\,\omega,
\]
where \(\eta>0\) is the overall step size and \(\omega>0\) is the primal weight. A safe choice stated in the survey literature is \(\eta<1/\|K\|_2\) with \(\omega=1\) [2506.02174].

A complementary formulation used in the distributed implementation writes the LP in saddle form
\[
\min_{x\in\mathcal X}\max_{y\in\mathcal Y}\mathcal L(x,y)=c^T x+y^T A x-p(-y;\ell_c,u_c),
\]
with
\[
x^{k+1}= \operatorname{proj}_{\mathcal X}\bigl(x^k-\tau(c-A^T y^k)\bigr),
\]
\[
z^{k+1}=A(2x^{k+1}-x^k),
\]
\[
y^{k+1}=y^k-\sigma z^{k+1}-\sigma \operatorname{proj}_{-\mathcal S}\bigl(\sigma^{-1}y^k-z^{k+1}\bigr),
\]
where \(\mathcal X=\{x:\ell_v\le x\le u_v\}\) and \(\mathcal S=\{s:\ell_c\le s\le u_c\}\). Convergence is monitored via three relative KKT errors \(r_{\text{primal}}, r_{\text{dual}}, r_{\text{gap}}\), and the solver stops when \(\max\{r_{\text{primal}},r_{\text{dual}},r_{\text{gap}}\}\le \epsilon\) [2601.07628].

## 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
\[
z^{k+1}=\frac{k+1}{k+2}\,\mathrm{PDHG}(z^k)+\frac{1}{k+2}z^0,
\]
and with reflection parameter \(\gamma\in[0,1]\),
\[
z_{\rm ref}=(1+\gamma)\,\mathrm{PDHG}(z^k)-\gamma z^k,\qquad
z^{k+1}=\frac{k+1}{k+2}z_{\rm ref}+\frac{1}{k+2}z^0.
\]
This is the restarted Halpern PDHG core described for cuPDLPx [2507.14051].

A defining enhancement is a restart criterion based on the fixed-point residual
\[
r(z)=\|z-\mathrm{PDHG}(z)\|_P.
\]
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 [2507.14051].

cuPDLPx also updates the primal weight \(\omega\) via a discrete PID control loop on \(\log \omega\). The tracking error is
\[
e^n=\log\!\Bigl(\frac{\sqrt{\omega^n}\,\|x^{n,t}-x^{n,0}\|_2}{(1/\sqrt{\omega^n})\,\|y^{n,t}-y^{n,0}\|_2}\Bigr),
\]
followed by
\[
\log\omega^{n+1}=\log\omega^n-\Bigl[K_P e^n+K_I\sum_{i=1}^n e^i+K_D(e^n-e^{n-1})\Bigr].
\]
The intended effect is to rebalance primal and dual progress at restarts [2507.14051].

The survey literature additionally describes diagonal scaling
\[
\tilde K=D_1 K D_2
\]
with \(D_1,D_2\succ 0\) chosen by a combination of Ruiz equilibration and the diagonal preconditioner of Pock–Chambolle, and an adaptive line-search rule for \(\eta\). 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 [2506.02174].

## 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 [2601.07628][2312.14832].

Its principal data structures are sparse matrix storage for \(A\) in CSR format and dense vectors for primal and dual variables, objective coefficients, bounds, and auxiliary state. In the distributed description, \(A\) 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 [2601.07628].

The kernel inventory is consistent across descriptions: sparse matrix–vector multiply for \(A x\) and \(A^T y\), 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 [2507.14051][2601.07628].

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 [2312.14832].

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 \(\gamma=1\), one chooses
\[
\mathcal T_1=\lambda_A I_m-AA^*,\qquad \lambda_A=1/\eta^2\ge \|A\|^2,\qquad \sigma=\eta/\omega,
\]
and identifies the relevant iterates so that the HPR-LP and cuPDLPx sequences coincide exactly [2509.23903].

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 \(O(1/k)\) convergence guarantees of HPR-LP apply directly to cuPDLPx’s base algorithm. In particular, the paper presents a KKT-residual estimate of the form
\[
\|\mathcal R(\bar w^{k+1})\|\le C\,\frac{R_0}{k+1}.
\]
This situates cuPDLPx within a broader operator-splitting lineage rather than treating it as an isolated solver design [2509.23903].

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 \(\eta\approx 1/\|A\|\), no adaptive penalty updates or restarts, possible stalling on poorly scaled problems, and open questions around degenerate steps \(\eta=1/\|A\|\). These caveats help explain why the enhanced cuPDLPx line places significant emphasis on restart and control mechanisms [2509.23903].

## 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 [2507.14051].

| Solver | MIPLIB \(\epsilon=10^{-4}\) SGM\(_{10}\) | Mittelmann \(\epsilon=10^{-4}\) SGM\(_{10}\) |
|---|---:|---:|
| 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 \(\epsilon=10^{-4}\), while remaining approximately 3x-5x at \(\epsilon=10^{-8}\) [2507.14051].

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 [2601.07628].

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 [2604.23951].

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 [2604.23951].

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 \(10^2\), \(10^3\), \(10^4\), and \(10^5\) iterations and, after each block, streams the current primal \(x^k\), dual \(y^k\), and reduced costs \(c-A^T y^k\) 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 [2605.05086].

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 [2605.05086].

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.

Source: https://www.emergentmind.com/topics/cupdlpx