Algorithm¶

cuPDLPx solves linear programs using restarted reflected Halpern primal–dual hybrid gradient (r²HPDHG)^[1]. Each iteration updates the primal and dual iterates using sparse matrix–vector products and projections onto the bound constraints.

LP formulation¶

Primal problem¶

cuPDLPx solves linear programs in bound form:

\[ \begin{aligned} \operatorname*{minimize}_{x\in\mathbb R^n}\quad & c^\top x + c_0 \\ \text{subject to}\quad & \ell_c \le Ax \le u_c, \\ & \ell_v \le x \le u_v. \end{aligned} \]

Here $A\in\mathbb R^{m\times n}$ is the constraint matrix, $c\in\mathbb R^n$ is the objective vector, and $c_0$ is the objective constant. The vectors $\ell_c,u_c$ and $\ell_v,u_v$ are the constraint and variable bounds, respectively. Bounds may be infinite.

Define the variable and constraint sets

\[ \mathcal X=\{x:\ell_v\le x\le u_v\},\qquad \mathcal S=\{s:\ell_c\le s\le u_c\}. \]

Dual problem¶

\[ \begin{aligned} \operatorname*{maximize}_{y\in\mathcal Y,\,r\in\mathcal R}\quad & -p(-y;\ell_c,u_c)-p(-r;\ell_v,u_v)+c_0 \\ \text{subject to}\quad & c-A^\top y=r. \end{aligned} \]

The support function of the interval $[\ell,u]$ is

\[ p(z;\ell,u)=u^\top z^+-\ell^\top z^-, \qquad z^+=\max(z,0),\quad z^-=\max(-z,0). \]

Let $\mathbb R_+=[0,+\infty)$ and $\mathbb R_-=(-\infty,0]$. The dual and dual-slack domains are the Cartesian products $\mathcal Y=\prod_{i=1}^m\mathcal Y_i$ and $\mathcal R=\prod_{j=1}^n\mathcal R_j$, where

\[ \mathcal Y_i= \begin{cases} \{0\}, & \ell_{c,i}=-\infty,\ u_{c,i}=+\infty,\\ \mathbb R_-, & \ell_{c,i}=-\infty,\ u_{c,i}\in\mathbb R,\\ \mathbb R_+, & \ell_{c,i}\in\mathbb R,\ u_{c,i}=+\infty,\\ \mathbb R, & \text{otherwise}, \end{cases} \qquad \mathcal R_j= \begin{cases} \{0\}, & \ell_{v,j}=-\infty,\ u_{v,j}=+\infty,\\ \mathbb R_-, & \ell_{v,j}=-\infty,\ u_{v,j}\in\mathbb R,\\ \mathbb R_+, & \ell_{v,j}\in\mathbb R,\ u_{v,j}=+\infty,\\ \mathbb R, & \text{otherwise}. \end{cases} \]

Here $r$ is the dual-slack variable.

Saddle-point formulation¶

The LP is equivalent to

\[ \min_{x\in\mathcal X}\max_{y\in\mathcal Y} L(x,y) =c^\top x+c_0-y^\top Ax-p(y;-u_c,-\ell_c). \]

The objective constant $c_0$ shifts both objective values without changing the feasible set or the PDHG iterates. The paper^[1] uses $c_0=0$.

Restarted reflected Halpern PDHG¶

The algorithm applies reflection and Halpern anchoring to the PDHG update, with adaptive restarts^[1].

Algorithm 1. Restarted reflected Halpern PDHG

Input: An LP, an initial point $(x^{0,0},y^{0,0})$, and a termination tolerance $\varepsilon$.

Precondition. Apply geometric mean scaling, Ruiz equilibration, Pock–Chambolle scaling, and objective and bound scaling.
Initialize. Set $$ \eta=\frac{0.998}{\lVert A\rVert_2},\qquad \omega^0=1,\qquad \tau=\frac{\eta}{\omega^0},\qquad \sigma=\eta\omega^0. $$
for $n=0,1,\ldots$ do
Set $k\leftarrow 0$.
repeat
PDHG update. $$ \begin{aligned} \widehat x^{n,k+1} &=\operatorname{proj}_{\mathcal X} \left(x^{n,k}-\tau(c-A^\top y^{n,k})\right),\\[0.4em] \widehat y^{n,k+1} &=y^{n,k} -\sigma A(2\widehat x^{n,k+1}-x^{n,k}) -\sigma\operatorname{proj}_{[-u_c,-\ell_c]} \left(\sigma^{-1}y^{n,k} -A(2\widehat x^{n,k+1}-x^{n,k})\right). \end{aligned} $$
Termination check.
if $\operatorname{KKT}(\widehat x^{n,k+1},\widehat y^{n,k+1}) \le\varepsilon$ then

return $(\widehat x^{n,k+1},\widehat y^{n,k+1})$.

end if
Reflection and Halpern update. $$ \begin{aligned} x^{n,k+1} &=\frac{k+1}{k+2} \left(\gamma\bigl(2\widehat x^{n,k+1}-x^{n,k}\bigr) +(1-\gamma)x^{n,k}\right) +\frac{1}{k+2}x^{n,0},\\ y^{n,k+1} &=\frac{k+1}{k+2} \left(\gamma\bigl(2\widehat y^{n,k+1}-y^{n,k}\bigr) +(1-\gamma)y^{n,k}\right) +\frac{1}{k+2}y^{n,0}. \end{aligned} $$
Set $k\leftarrow k+1$.
until a restart condition holds.
Restart.
Set the next anchors:
$$ (x^{n+1,0},y^{n+1,0}) =(\widehat x^{n,k},\widehat y^{n,k}). $$
Update $\omega^{n+1}$ with the PID controller, then update $\tau$ and $\sigma$.
end for

The algorithm uses the following notation:

$\operatorname{proj}_{\mathcal X}$: projection onto the variable bounds, computed by clipping each coordinate to its interval;
$\eta$: the step size, based on the scaled matrix's spectral norm; the active-set step size boost may increase it at a restart;
$\omega^n$: the primal weight for epoch $n$; the initial value $\omega^0=1$ assumes the default objective and bound scaling described under Preconditioning;
$\tau$ and $\sigma$: the primal and dual step sizes, respectively;
$\gamma$: the reflection coefficient (reflection_coefficient, default 1), which weights the reflected point $2\widehat x^{n,k+1}-x^{n,k}$ against the current iterate $x^{n,k}$;
$n$ and $k$: the restart epoch and the iteration within that epoch (here $n$ is an epoch index, rather than the variable count in the LP);
$(x^{n,0},y^{n,0})$: the primal and dual anchors for epoch $n$;
$(\widehat x^{n,k+1},\widehat y^{n,k+1})$: the PDHG iterate before the reflection and Halpern update, at which the KKT residual is evaluated.

The KKT (Karush–Kuhn–Tucker) check tests primal feasibility, dual feasibility, and the primal–dual gap. See Termination criteria for the residuals and tolerances.

The following pages describe the components: Base algorithm, Presolve and postsolve, Preconditioning, Step size, Adaptive restart, Primal weight, and Termination criteria. Optional feasibility polishing can reduce primal and dual feasibility residuals after the main solve.

References¶

[1] Haihao Lu, Zedong Peng, and Jinwen Yang. cuPDLPx: A Further Enhanced GPU-Based First-Order Solver for Linear Programming, 2025.