DuaLip: GPU-Accelerated LP Solver

• DuaLip is a large-scale linear programming solver that uses ridge-regularized dual ascent and first-order methods to efficiently solve ranking, allocation, and matching problems.
• It decouples problem specification from the optimization engine with composable primitives like ObjectiveFunction, ProjectionMap, and Maximizer, ensuring system extensibility.
• Its GPU-oriented architecture and algorithmic refinements, including Jacobi preconditioning and coordinate scaling, deliver robust convergence and over 10x wall-clock speedup compared to CPU implementations.

DuaLip is a large-scale linear programming solver family built around ridge-regularized dual ascent with first-order methods for decision systems such as ranking, allocation, and matching that must be solved repeatedly at massive scale. In its GPU-oriented redesign, DuaLip is presented as a solver architecture that decouples problem specification from the optimization engine, expresses LP formulations through composable primitives for dual objective evaluation and blockwise projection operators, and targets modern accelerators while maintaining convergence guarantees. On extreme-scale matching workloads, the GPU implementation achieves at least a 10x wall-clock speedup over the prior distributed CPU DuaLip solver under matched stopping criteria (Dexter et al., 4 Mar 2026).

1. Problem setting and solver formulation

DuaLip-GPU starts from the primal LP

$$\text{minimize } c^\top x \qquad \text{subject to } A x \le b,\; x \in C,$$

and adds an $\ell_2$-penalty $\gamma/2\|x\|_2^2$. Its Lagrangian dual is a smooth concave function

$$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$

with

$$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$

The report states the corresponding evaluation identities as

$$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$

and

$$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$

This formulation places the computational burden on sparse matrix-vector products and blockwise projections. The architectural consequence is that all sparse matrix-vector products $A x^*$, $A^\top \lambda$ and blockwise $\Pi_{(C)}$ are exposed as kernels; changing constraints or objective simply means registering a new ProjectionMap or ObjectiveFunction, without touching the solver loop. This suggests that extensibility is treated as a first-class systems property rather than an afterthought (Dexter et al., 4 Mar 2026).

2. Operator-centric architecture

The redesigned system adopts an operator-centric interface comprising three minimal primitives, explicitly intended to decouple problem formulation, projection logic, and solver internals. The style is compared to tensor frameworks such as PyTorch, where the solve loop orchestrates only sparse-dense operators and blockwise projections; new LP formulations register only the necessary kernels locally (Dexter et al., 4 Mar 2026).

Interface	Role	Method
Maximizer	Implements the dual ascent loop	maximize(obj, λ₀) → {λ\*, history}
ObjectiveFunction	Encapsulates $\ell_2$0 and the simple-constraint map	calculate(λ, γ) → {g(λ), ∇g(λ)}
ProjectionMap	Encodes per-block simple constraints such as simplex and box	project(block_id, v) → Π_Cblock(v)

In practice, the solve loop is described as $A^\top \lambda$4

The system-level significance of this decomposition is explicit in the technical report: new formulations can be added locally while reusing a shared optimization loop, diagnostics, and distributed infrastructure. That design directly addresses the limitation of the original implementation, which was tightly coupled to a small number of schemas and built on a CPU-centric Scala/Spark stack (Dexter et al., 4 Mar 2026).

3. Algorithmic refinements to ridge-regularized dual ascent

The report identifies three enhancements that stabilize and accelerate first-order maximization: Jacobi-style row normalization, primal coordinate scaling, and a continuation scheme for the regularization parameter (Dexter et al., 4 Mar 2026).

Jacobi-style row normalization defines

$\ell_2$1

over nonzero rows and replaces

$\ell_2$2

Since $\ell_2$3 has positive diagonal, $\ell_2$4. The dual Hessians satisfy

$\ell_2$5

so scaling is exactly Jacobi preconditioning of the dual Hessian. In the matching case, where $\ell_2$6 is block-diagonal by source, this clusters eigenvalues near one. The report further states that, under mild i.i.d. row-row correlation $\ell_2$7,

$\ell_2$8

Primal coordinate scaling is introduced for the regime in which primal entries $\ell_2$9 vary widely and a uniform $\gamma/2\|x\|_2^2$0 distorts geometry. With a diagonal scaling $\gamma/2\|x\|_2^2$1 and substitution $\gamma/2\|x\|_2^2$2, the equivalent problem becomes

$\gamma/2\|x\|_2^2$3

so that

$\gamma/2\|x\|_2^2$4

Choosing $\gamma/2\|x\|_2^2$5 balances the curvature of $\gamma/2\|x\|_2^2$6 across dual coordinates.

The continuation scheme addresses the fact that very small $\gamma/2\|x\|_2^2$7 gives accurate recovery of the original LP but poor smoothness, with Lipschitz constant

$\gamma/2\|x\|_2^2$8

The schedule initializes $\gamma/2\|x\|_2^2$9 at a moderate value and decays it geometrically:

$$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$0

for example $$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$1 every $$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$2 iterations. Simultaneously, the AGD step-size is scaled as $$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$3 to maintain stable accelerated updates. Early iterations progress rapidly on a smoother problem, while later stages refine the solution to within $$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$4 of the true LP optimum.

4. GPU realization for sparse matching constraints

The GPU implementation is specialized to sparse matching constraints. Matching LPs are described as having $$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$5 assembled from diagonal $$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$6 blocks $$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$7 for $$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$8 constraint families and sources $$g(\lambda)=\min_{x\in C}\left[c^\top x+\frac{\gamma}{2}\|x\|^2+\lambda^\top(Ax-b)\right],$$9, with further sparsity from ineligible $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$0 pairs (Dexter et al., 4 Mar 2026).

The first execution technique is a constraint-aligned CSC layout. The system stores the $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$1 matrix $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$2 in GPU memory as a block-CSC tensor whose columns correspond to destination coordinates $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$3 across all sources $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$4. Each nonzero entry in column $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$5 thus relates to the $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$6 variables $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$7 contiguously. At $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$8 storage, the layout incurs one pointer per column, one index per nonzero, and one value per nonzero, enabling coalesced SpMV via PyTorch’s sparse-dense kernels.

The second technique is a batched projection kernel for per-slice projections $$\nabla g(\lambda)=A x^*_{(\gamma)}(\lambda)-b,\qquad x^*_{(\gamma)}(\lambda)=\Pi_C\!\left(-(A^\top\lambda+c)/\gamma\right).$$9. Rather than launching one CUDA kernel per slice, the implementation buckets sources by slice length into ranges $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$0, gathers the corresponding source vectors into a dense $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$1 tensor with padding, launches a single high-occupancy batched-projection kernel on $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$2 vectors, and scatters results back. The report states that this requires $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$3 kernels, amortizes launch overhead, and keeps padding waste under $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$4.

The third technique is distributed dual communication. In multi-GPU or multi-node setups, columns of $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$5 are partitioned evenly across devices, while $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$6 and $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$7 are replicated. Each iteration consists of local computation of $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$8, local $$x^* = \Pi_C (-(A^\top \lambda + c)/\gamma), \qquad \nabla g(\lambda) = A x^* - b,$$9, and local $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$0; a single NCCL All-Reduce aggregating $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$1 and two scalars; an AGD update on GPU 0; and two NCCL Broadcasts sending updated $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$2 and momentum states to all GPUs. Per iteration communication is therefore one $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$3-vector reduce and two $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$4-vector broadcasts, independent of $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$5 or the per-GPU partition, and the report characterizes the resulting scaling as near-ideal (Dexter et al., 4 Mar 2026).

5. Empirical behavior and performance

The experimental evaluation benchmarks synthetic matching LPs with $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$6 up to 100 M sources, $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$7 K destinations, and sparsity $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$8, comparing Scala/Spark DuaLip on CPU with PyTorch/GPU DuaLip on 1–4 GPUs. All runs use Accelerated Gradient Descent with fixed hyperparameters: initial step $$g(\lambda) = c^\top x^* + (\gamma/2)\|x^*\|^2 + \lambda^\top(A x^* - b).$$9, max step $A x^*$0, and $A x^*$1. The stopping criterion is a fixed dual-objective gap (Dexter et al., 4 Mar 2026).

Implementation parity is evaluated by plotting dual objective versus iteration for Scala versus PyTorch on 1–4 GPUs. The trajectories overlap almost exactly, with relative error below $A x^*$2 after 100 iterations, which the report interprets as confirming numerical fidelity.

Wall-clock results are given as average time per AGD iteration in seconds. For 25 M sources, Scala CPU requires 2.46 s, while 1 GPU requires 0.27 s and 4 GPUs require 0.07 s. For 50 M sources, Scala CPU requires 3.44 s and 4 GPUs require 0.13 s. For 75 M sources, Scala CPU requires 2.63 s and 4 GPUs require 0.21 s. For 100 M sources, Scala CPU requires 3.33 s and 4 GPUs require 0.27 s. The report summarizes these measurements by stating that even on a single GPU, iteration time is $A x^*$3 lower than Scala, and that sharding across 4 GPUs yields nearly linear speedup, specifically $A x^*$4 versus ideal $A x^*$5, while allowing instances beyond CPU memory limits.

The impact of algorithmic enhancements is also quantified. With Jacobi preconditioning, early-stage convergence is $A x^*$6–$A x^*$7 faster on a 25 M-source problem; without preconditioning, heterogeneity in row norms slows descent. A fixed $A x^*$8 is compared with geometric decay in which $A x^*$9 is halved every 25 iterations; continuation accelerates early progress on the smoother problem and converges to the same final gap, delivering $A^\top \lambda$0 speedup in total iterations (Dexter et al., 4 Mar 2026).

6. Convergence properties and feasibility control

The report links dual progress to primal feasibility through Appendix Lemma A.2. If dual suboptimality is defined as

$A^\top \lambda$1

then primal infeasibility satisfies

$A^\top \lambda$2

This bound is important because DuaLip is motivated by operational decision systems in which feasibility tolerances have direct practical meaning (Dexter et al., 4 Mar 2026).

Empirically, both preconditioning and continuation reduce $A^\top \lambda$3 rapidly, ensuring feasibility within economically meaningful tolerances in under 200 iterations. In that sense, the solver’s convergence story is not limited to objective ascent speed; it is tied explicitly to constraint satisfaction in the original LP.

The report’s overall conclusion is correspondingly systems-oriented rather than purely algorithmic: DuaLip-GPU attains at least 10x wall-clock acceleration over the CPU-based DuaLip under matched stopping criteria, with predictable single- and multi-GPU scaling and robust convergence across extreme-scale matching workloads (Dexter et al., 4 Mar 2026).

7. Nomenclature and potential ambiguity

A common source of confusion is the near-homography between DuaLip and DualLip. The former denotes the LP solver lineage discussed above, including the GPU redesign for ridge-regularized dual ascent on large-scale matching-style LPs (Dexter et al., 4 Mar 2026). The latter denotes a different 2020 system, “DualLip: A System for Joint Lip Reading and Generation,” concerned with lip reading, lip generation, and talking face generation rather than linear programming (Chen et al., 2020).

That distinction is substantive, not merely orthographic. “DualLip” in the 2020 paper treats lip reading and lip generation as dual tasks in a semi-supervised learning framework, using pseudo pairs derived from unlabeled text and lip video data, and reports results on GRID and TCD-TIMIT (Chen et al., 2020). By contrast, “DuaLip-GPU Technical Report” addresses sparse LP optimization, operator-centric solver design, GPU kernels, NCCL-based distributed communication, and convergence acceleration for ridge-regularized dual ascent (Dexter et al., 4 Mar 2026).

For technical discourse, especially in bibliographic and systems contexts, maintaining this distinction avoids conflating an optimization infrastructure for large-scale LPs with a multimodal audiovisual learning system.