arXiv CSopen access

Large-Scale Regularized Matching on GPU Clusters

Large-Scale Regularized Matching on GPU Clusters Aida Rahmattalabi∗1 , Gregory Dexter∗†2 , Sanjana Garg∗1 , Qinquan Song†3 , Shenyinying Tu1 , Yuan Gao1 , Zhipeng Wang1 , and Rahul Mazumder1

arXiv:2606.07777v1 [cs.DC] 5 Jun 2026

1

LinkedIn Corporation, USA 2 Nubank, USA 3 OpenAI, USA

Abstract Production decision systems such as ad allocation and content matching involve millions of users and thousands of items, reducing to large-scale linear programs with a common structure: sparse constraint matrices with diagonal block structure across sources (users), solved repeatedly on recurring cadences over slowly evolving inputs. Three systems gaps stand out. Scale: production instances routinely exceed the memory capacity of GPU solvers such as cuPDLP and D-PDLP under fixed hardware budgets. Temporal instability: solution variability across consecutive runs induces downstream churn and complicates SLAs, a failure mode for which existing solvers provide no explicit control. Extensibility: CPU-based solvers such as DuaLip-Scala converge slowly and couple problem formulation to fixed schemas, making new constraint families expensive to incorporate. We present a distributed multi-GPU LP solver built natively in PyTorch with systems–algorithm co-design for this structure. It adopts a column-sharded execution model with fused Triton kernels and batched operations to reduce per-iteration overhead. As the number of users grows, only local computation increases, while communication is limited to a single reduction of item-level dual variables to rank 0, yielding near-linear scaling with GPU count at fixed item size. Second, we adopt ridge-regularized LPs from prior work to improve solver stability, introducing a tunable γ that provably bounds run-to-run solution drift. This control is not exposed in existing GPU LP solvers. A continuation schedule over γ further balances convergence speed and solution fidelity. Finally, we introduce an operator-centric programming model that replaces DuaLip-Scala’s schema-bound interface with three composable primitives (ObjectiveFunction, ProjectionMap, Maximizer), enabling new formulations through local changes without modifying the solve loop or distributed infrastructure. On synthetic matching workloads representative of production usage, our system achieves order-of-magnitude wall-clock speedup over DuaLip-Scala, near-linear multi-GPU scaling (3.86× on 4 GPUs), while scaling to problem sizes beyond the reach of GPU-based solvers.

Keywords: linear programming, GPU computing, distributed optimization, dual ascent, matching, ridge regularization

∗ Three authors contributed equally to this research. † Work completed at Linkedin

1

1

Introduction

Large-scale linear programs (LPs) underpin production decision systems, where instances with billions of decision variables and hundreds of millions of constraints are solved repeatedly under tight latency and stability requirements. Classical LP solvers such as simplex or interior-point methods do not scale well to this regime: both are bottlenecked by sparse matrix factorizations, which are memory-intensive and difficult to parallelize efficiently at scale. This leads to poor utilization of modern parallel hardware and limits their applicability in large-scale production settings. Another challenge is stability. In these settings, small perturbations in the data can lead to large changes in the solution. To address this, regularization terms are often introduced to encourage proximity to a reference solution (e.g., the prior-day solution) and reduce solution drift. However, in the linear setting, this effectively introduces additional constraints into an already large-scale problem, further increasing computational burden. The LPs arising in these settings share a structure that existing solvers do not fully exploit. Decision variables partition into per-source blocks (one per user), each subject to simple local constraints—such as simplex or box constraints—that decompose independently across blocks. Coupling constraints (e.g., budgets, pacing limits, and frequency caps) interact with each source block in an element-wise manner, yielding a constraint matrix with diagonal block structure across sources. The resulting matrices are highly sparse: nonzeros occur only at eligible source–destination pairs, and within each block the interaction is diagonal. Recently, first-order methods have emerged as a scalable approach for large linear programs. These are iterative primal–dual algorithms, notably Primal-Dual Hybrid Gradient (PDHG) [2], whose per-iteration cost reduces to sparse matrix–vector multiplication, making them well-suited to GPUs. This has enabled solvers such as cuPDLP.jl [11] and cuPDLPx [9], the latter improving convergence via restarted Halpern PDHG. However, these methods are strictly single-device, with memory-bound problem sizes and no path to scale-out. D-PDLP [7] extends PDHG to multiple GPUs via 2D matrix partitioning and NCCL AllReduce, achieving up to 6× speedup on NVLink-connected H100s, but remains confined to a single NVLink domain and requires synchronous communication at every iteration, with no mechanism for stability control across recurring solves. Crucially, these methods ignore structural properties of production LPs, where variables decompose into independent per-source blocks with simple local constraints, and coupling constraints act element-wise, inducing a diagonal block structure. Existing solvers treat the system as unstructured, missing this decomposition. A complementary line of work, including ECLIPSE [3] and DuaLip [14], exploits this structure via ridge-regularized dual ascent and scales to CPU clusters, but remains limited to fixed schemas and CPU-centric designs, making extensions costly and limiting accelerator efficiency. We build on DuaLip [14] and present a distributed multi-GPU LP solver built natively in PyTorch that co-designs algorithm and system components around the diagonal block structure of production matching LPs. Our contributions are: 1. Column-sharded SpMV and projection via NCCL. We partition the constraint matrix across GPUs by source columns. Because constraints decompose per source, each GPU computes local primal updates and gradient contributions independently, with no cross-device dependencies during computation. Each iteration requires only a reduction and broadcast of destination-level dual variables to and from rank 0, whose size is independent of the number of sources and the column partitioning. As a result, scaling the number of sources increases only local computation while communication remains fixed, enabling near-linear scaling on instances that exceed single-device memory capacity. 2. Tunable regularization with stability control. We adopt a ridge-regularized formulation from prior work and expose the parameter γ as a first-class control, which bounds run-to-run primal drift and directly addresses stability requirements in recurring production solves. A continuation schedule initializes with larger γ for fast, stable convergence and anneals it toward the target fidelity. 3. Operator-centric programming model. We replace DuaLip-Scala’s rigid schemas with three composable primitives: ObjectiveFunction, Maximizer and ProjectionMap. New LP formulations require only local implementation of the objective and gradient, while the solve loop and distributed infrastructure remain unchanged, improving extensibility without affecting performance. On synthetic matching workloads calibrated to production-scale sizes, our system achieves a 9.1× wall-clock

2

speedup over DuaLip-Scala on a single GPU and exhibits near-linear scaling that improves with the number of sources.

2

Related Work

Large-scale matching and allocation. A substantial literature addresses weighted bipartite matching and network-flow problems at scale, with specialized solvers that exploit degree and capacity constraints [17, 6, 16]. However, these methods do not naturally support heterogeneous linear side constraints such as pacing, fairness, and frequency caps that arise in production allocation systems [8, 4]. The formulation inherited from ECLIPSE preserves general LP expressivity while exploiting matching structure in the dual, enabling new constraint families to be added without modifying the core solver. Production applications include notification and email allocation [3], ad allocation matching [12], opportunity assignment [14], and inactive-member promotion on social networks [1]. First-order methods. PDLP [2] established convergence guarantees for restarted PDHG on general LPs and underpins both cuPDLP variants [10, 9]. These methods reduce per-iteration computation to generic sparse matrix–vector products but do not exploit problem structure. In addition, their analysis does not address solution stability across repeated solves, a key requirement in production decision systems. Ridge-regularized dual ascent. Basu et al. [3] showed that ℓ2 regularization induces smoothness in the dual and enables scalable optimization for extreme-scale allocation, identifying the diagonal-block structure of webscale matching LPs as the key property for distributed decomposition. Ramanath et al. [14] operationalized these ideas in a Scala/Spark system (DuaLip), with efficient projections for simplex and box constraints. However, the system is CPU-centric and does not support accelerator-friendly sparse layouts or fine-grained collective communication. It also exposes only two fixed schemas, making new constraint families require non-local code changes. GPU-based LP solvers. Lu and Yang [10] showed that restarted PDHG can match commercial solvers when instances fit in GPU memory, but is fundamentally limited by a single-device memory ceiling. Lu et al. [9] improved convergence via restarted Halpern PDHG, achieving 2.5×–5× speedups over cuPDLP, while retaining the same memory constraint. Li et al. [7] extended PDHG to multiple GPUs via 2D matrix partitioning and NCCL-based synchronization, achieving up to 6× speedup on eight NVLink-connected H100s within a single server. However, these methods still fail to exploit the structure of production matching LPs. First, all reported experiments are confined to a single NVLink domain, leaving multi-node scalability unaddressed. Second, treating the constraint matrix as unstructured forces a generic 2D partitioning that requires two synchronous AllReduces per iteration along the row and column grid. In our setting, the column/source decomposition makes the per-iteration reduce volume scale with the constraint dimension (destinations) and independent of the decision-variable dimension (sources) which is the larger of the two in production matching workloads. High-performance tensor frameworks. Our implementation builds on PyTorch [13], leveraging sparse-dense and fused kernels, batched dense operations, and torch.distributed with an NCCL backend. This provides a flexible execution substrate for structured sparse linear algebra, bridging the gap between LP-specific algorithm design and general-purpose GPU programming. Further, it enables a unified implementation in which the solve loop and projection operators are shared across formulations while remaining efficient on distributed GPU hardware.

3

Problem Setting

This section formalizes the ridge-regularized LP used throughout the paper, derives its dual, and specializes both to the matching workloads that motivate our system. The resulting computational object consists of two sparse matrix–vector products and a blockwise projection over a constraint matrix with diagonal block structure, which forms the basis for the design in Section 4.

3

3.1

Ridge-Regularized LP

We consider the regularized LP problem min x∈C

c⊤ x +

γ ∥x∥22 2

s.t.

Ax ≤ b,

(1)

with primal variable x ∈ Rn , sparse coupling matrix A ∈ Rm×n , right-hand side b ∈ Rm , and regularization parameter γ > 0. The set C ⊆ Rn encodes simple constraints that decompose blockwise across sources and admits an efficient projection oracle ΠC ; the inequality Ax ≤ b encodes the coupling constraints that bind blocks together. This distinction is operationally central: simple constraints are absorbed into the primal subproblem via projection and do not appear in the dual, while coupling constraints are the only ones associated with explicit dual variables. The ridge term induces strong convexity in the primal, yielding a smooth dual with a fixed step size that can be derived analytically from A and γ and it preserves the blockwise separability of the inner minimization, which the column-sharded execution layer of Section 4 relies on.

3.2

Dual Formulation

Let λ ∈ Rm ≥0 denote the dual variables associated with the constraint Ax ≤ b. The Lagrangian dual is given by n o γ g(λ) = min c⊤ x + ∥x∥22 + λ⊤ (Ax − b) . (2) x∈C 2 Strong convexity of the primal objective (due to γ > 0) implies strong duality; in other words, solving the dual recovers the primal optimum with no duality gap. The inner minimization decomposes over the blocks of C and admits a closed-form expression as a projection of an unconstrained quadratic minimizer: 
 xγ∗ (λ) = ΠC − γ1 A⊤ λ + c . (3) By Danskin’s theorem, the dual is differentiable with gradient ∇g(λ) = Axγ∗ (λ) − b.

(4)

Equations (3)–(4) isolate the three dominant per-iteration operations: the sparse matrix–vector products A⊤ λ and Axγ∗ , and the projection ΠC . Their cost is fully determined by the structure of A, which we exploit in the next section.

3.3

Matching Workloads

The LPs arising in our target workloads share a structure that the rest of the paper exploits. Let xij denote the assignment of source i ∈ [I] (a user) to destination j ∈ [J] (a campaign, item, or slot), and let E ⊆ [I] × [J] denote the set of eligible pairs. We define xij = 0 for (i, j) ∈ / E and stack the per-source blocks xi = (xi1 , . . . , xiJ )⊤ into vector x = (x1 ; . . . ; xI ). P Simple constraints. Each source has a local constraint, typically a unit simplex constraint j xij ≤ 1 with xij ≥ 0 for (i, j) ∈ E, or a box variant, which decomposes independently across sources. These correspond to the simple constraints x ∈ C in (1), where each block xi lies in a separate polytope Ci and the projection operator ΠC decomposes as a Cartesian product of per-block projections [14]. Crucially, these constraints do not introduce dual variables; only the coupling constraints do. As a result, the number of explicit dual variables scales as O(J) rather than O(I + J), with J ≪ I in production settings, which directly determines the communication volume in Section 4. Coupling constraints. Cross-source requirements, including destination-level budgets, pacing, frequency caps, and fairness limits, couple all sources eligible for a given destination. A canonical example is destination capacity, X aij xij ≤ dj ∀j ∈ [J], (5) i∈[I]

which interacts with each source block in an element-wise manner. With x stacked in source-major order, this yields a coupling matrix formed by horizontally concatenated diagonal blocks. We allow multiple such constraint families to coexist, producing the structure exploited in Section 4. 4

Definition 1 (Matching coupling matrix). For m coupling-constraint families indexed by k ∈ [m] and x = (x1 ; . . . ; xI ) ∈ RIJ , the coupling matrix A ∈ RmJ×IJ is   D11 D12 · · · D1I  D21 D22 · · · D2I    (6) A= . .. ..  , ..  .. . . .  Dm1 Dm2 · · · DmI where Dki ∈ RJ×J diagonal for every (k, i). Each block Dki acts element-wise on xi , and is zero on coordinates j with (i, j) ∈ / E. This object generalizes the ECLIPSE formulation [3], which supported a single matching block, to an arbitrary number of constraint families, each sharing a diagonal block-across-sources structure. It is sparse in two ways: structurally, since off-diagonal entries within each Dki are zero by construction, and structurally in the data, since aij = 0 for ineligible pairs. The rest of the paper treats this dual sparsity as the central design constraint. Section 4 encodes it in a compact CSC layout that materializes only the J diagonal entries of each block and shards A across GPUs along the source axis, so per-iteration communication scales with the dual dimension mJ rather than the number of nonzeros or sources.

4

GPU System Design

The dual ascent of Section 3 reduces each iteration to two sparse matrix–vector products and a blockwise projection. This section describes how we implement these primitives on a multi-GPU cluster. Our goal is to approximate the performance of a custom C++/CUDA solver while retaining a lightweight implementation built on PyTorch sparse tensors and torch.distributed, so that the same code path applies across all formulations supported by the programming model of Section 5. We make four design choices, each driven by structural properties of the matching LPs in Definition 1: a sparse tensor layout preserving block-diagonal and within-block sparsity (§4.1); a bucketed batching scheme that maps per-block projections to a small number of high-occupancy kernel launches (§4.2); a fused Triton kernel for simplex projection that collapses a multi-op PyTorch pipeline into a single kernel (§4.3); and a column-sharded execution layer whose per-iteration communication depends only on the number of coupling constraints (§4.4).

4.1

Sparse Tensor Layout

The constraint matrix of Definition 1 exhibits two sources of sparsity: structured sparsity from the diagonal block form across sources, and unstructured sparsity within each block, since only a subset of source–destination pairs are eligible. We exploit both by collapsing the diagonal block structure into a compact representation. With a single matching constraint family, A = [D1 , . . . , DI ], where each Di ∈ RJ×J is diagonal, so the nonzero pattern is fully captured by the J diagonal entries per block. We store these in a sparse tensor T in Compressed Sparse Column (CSC) format with columns indexed by source i, so that T [:, i] = diag(Di ) and all variables for a given source are stored contiguously in memory. CSC stores exactly the required information at fixed precision: one column pointer, one row index per nonzero, and one value per nonzero, with no representation of structural zeros within blocks. This layout directly supports the rest of the pipeline. The sparse products Ax and A⊤ λ use PyTorch optimized sparse-dense kernels without format conversion, and columns of T align with the destination partitioning of Basu et al. [3], making block-local operations column-local.

4.2

Batched Projection

The closed-form primal step (3) requires projecting each per-source slice of the candidate primal onto its feasible set Ci , a unit simplex or box depending on the formulation. A direct GPU implementation encounters two bottlenecks. Per-source kernel launches incur prohibitive dispatch overhead at scale, while a single dense slab representation of size [smax × num_sources] eliminates launch overhead but introduces substantial 5

zero-padding waste due to heterogeneous slice lengths. We address both issues with a logarithmic bucketing scheme. Slices are grouped by length into ranges [2t−1 , 2t ). Within each bucket, we pack slices into a dense slab padded only to the bucket upper bound, apply a single batched projection kernel, and scatter results back. This bounds padding overhead to at most a factor of two within each bucket, while reducing the number of GPU launches to 1 + ⌊log2 smax ⌋. The scheme integrates directly with PyTorch batched dense kernels. Figure 2 reports per-iteration runtime and peak memory as a function of source count, comparing bucketing against the single-slab baseline.

4.3

Fused Simplex Projection

P Within each bucket, the dominant projection in our target workloads is onto the simplex {w ≥ 0, i wi ≤ z}. The standard algorithm of Duchi et al. [5]—sorting, prefix sums, threshold selection, recovery of the cutoff index ρ and threshold θ, followed by the final subtract-and-clamp projection—maps naturally to high-level tensor frameworks but compiles into multiple GPU kernel launches per invocation. For block sizes L ranging from hundreds to a few thousand, launch overhead and intermediate memory traffic dominate arithmetic, as temporary tensors (sorted values, prefix sums, and masks) are repeatedly materialized in and read back from global memory. We fuse this pipeline into a single Triton kernel [15]. Each program instance processes one column: it loads the slice with masked bounds handling, performs the sort and inclusive scan in registers, recovers the cutoff index ρ via a boolean reduction exploiting the monotone prefix structure of the Duchi condition, computes the threshold θ through a masked reduction, and applies the final subtract-and-clamp projection. This eliminates intermediate materialization and global memory traffic between stages. For the inequality variant, we add an in-kernel early exit that returns the input unchanged when already feasible, using the same numerical tolerance as the reference implementation to preserve numerical equivalence. The kernel runs in fp32, which best matches Triton’s register-resident sort and scan primitives, and supports column lengths up to 8192. Beyond this limit, execution falls back to the multi-launch implementation to avoid register spilling. The fused kernel serves as the default projection primitive in the inner loop, reducing each primal update to a single kernel launch and eliminating intermediate materialization in global memory. Figure 1 compares per-iteration runtime and peak memory against the PyTorch-eager baseline across workloads ranging from 1M to 50M sources.

4.4

Distributed Execution

We launch one OS process per GPU via torchrun and form a single torch.distributed process group over NCCL. The same runtime supports both single-node multi-GPU and multi-node deployments through standard rendezvous environment variables. The CSC tensor T and corresponding entries of c are partitioned across devices in a balanced column split, while the dual variable λ and constraint vector b are replicated on every rank. Each rank reads the shared instance directly from the network filesystem and materializes only its local column block, so the full matrix is never assembled on a single device and no startup scatter is required. The column partition aligns with the decomposition induced by the simple constraints. For a fixed λ, each rank computes its local slice of xγ∗ (λ) via (3) and its contribution to ∇g(λ) independently, with no cross-device dependencies during computation. Each iteration proceeds as follows. First, each rank computes its local gradient contribution, objective value, and regularization term on-device. Next, a reduce-to-rank-0 (SUM) aggregates the gradient vector of size |λ|, while scalar reductions aggregate the objective and regularization values. Rank 0 then performs the accelerated gradient update, producing the next iterate and momentum vector, followed by projection onto the nonnegative orthant for inequality coordinates. Finally, the updated dual variables are broadcast from rank 0 to all ranks for the next iteration. The per-iteration communication volume therefore consists of one |λ|-sized reduction, one |λ|-sized broadcast, and O(1) scalar reductions, independent of the constraint-matrix nonzero count and the number of source partitions. Increasing the number of sources or GPUs therefore increases only local computation while communication remains fixed. Although the dual update is serialized on rank 0, its cost is negligible relative to the parallel sparse matrix–vector products and projections.

6

Component

Interface / Semantics

ObjectiveFunction

Encodes (A, b, c). calculate(λ, γ) returns (g(λ), ∇g(λ), xγ∗ (λ)), where xγ∗ (λ) = ΠC (− γ1 (A⊤ λ + c)) and ∇g = Axγ∗ − b. Implements only tensor-level ops; reuses CSC layout (Sec. 4.1), fused projection (Sec. 4.3), and block structure (Def. 1). Blockwise projection operator ΠC . Provides batched primitives for simplex / box / box-cut (Sec. 4.2). User-defined only for new constraint families; execution and batching infrastructure are reused. Runs dual ascent on λ ≥ 0: iterates ObjectiveFunction.calculate, applies accelerated update, conditioning, and continuation (Sec. 6), and hides distributed execution (Sec. 4.4).

ProjectionMap

Maximizer

Table 1: Programming model with shared ObjectiveFunction, projection, and optimization components.

5

Programming Model

The system is reusable across formulations only if new LPs can be expressed without modifying the solve loop or distributed runtime. We achieve this by exposing a boundary at the three primitives from Section 3.2: sparse products Ax and A⊤ λ, and the blockwise projection ΠC . All algorithmic logic (dual ascent, gradient evaluation, column-sharded execution, and NCCL communication) is shared; only formulation-specific operators are implemented locally. DuaLip-Scala instead places the boundary at the LP-family level: users select one of two declarative schemas, and gradient computation is dispatched via schema-specific implementations. This makes even P small changes non-local. For example, adding a global-count constraint (i,j)∈E xij ≤ M only augments A and λ, leaving the dual algorithm unchanged, yet requires modifications to the schema, parser, gradient code, Spark execution plan, and dispatch logic. The abstraction therefore tracks problem families rather than computation. Our system lowers the boundary to the gradient computation itself, implemented as tensor operations over sparse and dense structures. New coupling constraints correspond to a local change: one additional dual coordinate, one term in A⊤ λ, and one gradient contribution from x. The solver, projections, and distributed execution remain unchanged, as they operate purely on x and λ.

6

Algorithmic Improvements

Our solver builds on the ridge-regularized dual ascent framework of Section 3.1, optimizing the smoothed dual g(λ) via first-order updates with gradient access ∇g(λ) = Axγ∗ (λ) − b. Prior Scala/Spark implementations required instance-specific tuning of AGD variants and regularization strength to obtain stable convergence. We instead adopt a single optimization strategy based on continuation over the ridge parameter γ. The method solves a sequence of regularized problems, starting from a large γ to improve conditioning and stabilize early iterations, and gradually decreasing it to the target objective. Each stage is warm-started from the previous dual iterate, ensuring continuity across the schedule. This removes the need for instance-specific optimizer tuning and improves robustness across heterogeneous matching workloads of Definition 1. In practice, it combines fast initial progress under strong regularization with accurate solutions at small γ. Within each stage, we apply Jacobi preconditioning to improve conditioning by rescaling the complex constraints using the diagonal of the Hessian proxy induced by A. This improves the effectiveness of first-order updates, particularly in poorly conditioned regimes at intermediate values of γ. Together, continuation and Jacobi preconditioning yield a stable optimization trajectory from heavily smoothed problems to the original objective, improving convergence without manual hyperparameter tuning. Further details are provided in Appendix B.2.

7

Experiments

We evaluate our system along three axes: (i) system-level performance and multi-GPU scaling; (ii) comparisons with state-of-the-art GPU-solvers; and (iii) the effect of algorithmic enhancements, including preconditioning 7

DuaLip (PyTorch) Sources

Scala

1 GPU

2 GPUs

3 GPUs

4 GPUs

25M 50M 75M 100M

2.46 3.44 2.63 3.33

0.27 -

0.14 0.27 0.43 -

0.09 0.18 0.29 0.37

0.07 0.13 0.21 0.27

Table 2: Average time per AGD iteration (seconds). Multi-GPU sharding enables larger instances and yields substantial speedups over the Spark-based Scala implementation. and regularization continuation. We further demonstrate numerical parity with the Scala-based solver in the Appendix Section B. We use synthetic matching data to enable controlled scaling of problem size and sparsity. The data generation procedure and complete experimental setup are described in Section A of the Appendix.

7.1

System Performance and Scaling

Cross-Platform (CPU–GPU) Runtime and Scaling Efficiency We compare per-iteration runtime between Scala (Spark-based) and the PyTorch-based GPU implementation. Using a fixed random seed, we generate identical problem instances compatible with the Scala solver input schema and evaluate both systems under equal configurations. Table 2 reports the average time per algorithm iteration across 1000 iterations. For moderate problem sizes (25M sources), a single GPU already provides close to an order-of-magnitude improvement over Scala. As the problem size increases, model sharding across multiple GPUs becomes necessary to satisfy memory constraints. Beyond enabling larger instances, our multi-GPU parallelization significantly reduces iteration time, achieving more than 35× speedup relative to Scala and near-linear scaling within the GPU setting. Kernel-Level Optimizations. We compare the fused projection kernel of Section 4.3 against the multilaunch PyTorch reference implementation that it replaces. Figure 1 reports per-iteration time and peak GPU memory on a single H100, sweeping problem sizes from 1M to 50M sources. The fused kernel achieves a 2.5–5× per-iteration speedup for 10M–50M sources, and over 20× on the 1M instance, where kernel launch and dispatch overhead dominates the relatively small arithmetic workload on per-column slices. Peak memory is reduced by a consistent ∼ 20% across all settings above 1M, reflecting the removal of materialized intermediate tensors (sorted values, prefix sums, and masks) from global memory. Both implementations scale linearly with problem size, so the memory reduction translates directly into an increase in the maximum single-GPU feasible instance before saturating the 80 GiB device budget. Projection Batching. Figure 2 validates the bucketing scheme described in Section 4.2. We compare the bucketed projection against the single-slab baseline (batching=False in our implementation, corresponding to the second failure mode in Section 4.2) on top of the Triton-fused kernel, sweeping problem size from 10M to 50M sources. Bucketing delivers a consistent ∼1.2× per-iteration speedup by avoiding wasted compute on the zero-padded entries of the single slab, and a consistent ∼24% reduction in peak memory by replacing the global [smax × num_sources] slab with much smaller per-bucket slabs. Both gains are governed by the skew of the per-source slice-length distribution: heavier tails enlarge the fraction of arithmetic and memory the unbucketed path spends on padding. Multi-GPU and Multi-Node Scaling. We evaluate the distributed execution scheme of Section 4.4 over 1, 2, 4, 8, and 16 GPUs. The 1–8 GPU configurations run on a single H100 node, while 16 GPUs span two nodes coordinated via a torchrun-based launcher. All runs use the Triton-fused projection with batching enabled. Figure 3 reports end-to-end solve time and speedup across problem sizes. With a fixed destination count of 10,000, per-iteration communication is independent of source size (Section 4.4), so scaling is dominated by local gradient computation. At 16 GPUs, we observe 13.1× speedup on 75M sources (82% efficiency), 12.0× on 50M (75%), and 9.7× on 25M (61%). Smaller problems saturate 8

Figure 1: Triton fused projection vs. PyTorch-eager Duchi (single H100). (Top) p95 per-iteration time; (Bottom) peak GPU memory. earlier because the fixed reduce-and-broadcast overhead becomes a larger fraction of each iteration (6.0× at 10M, 37% efficiency). The scaling remains smooth across the 8 → 16 GPU transition, indicating that the per-iteration communication—one |λ|-reduce and two |λ|-broadcasts— is well within NCCL bandwidth on the H100 interconnect. Beyond speedup, distributed execution enables instances that exceed single-GPU memory. The 100M-source problem OOMs on a single device (working set exceeds 80 GiB), but runs on ≥ 2 GPUs due to linear column partitioning of the constraint matrix. At 16 GPUs, the 100M-source instance completes in under 3 minutes.

7.2

Comparison with D-PDLP

We compare our system against the state-of-the-art GPU primal–dual solver D-PDLP [7] on matching LPs generated with varying number of sources from 50M to 100M. Both solvers use float64 data type. The base formulation yields LPs with up to one billion nonzeros. We also consider an ℓ1 -regularized variant of the problem, where the ℓ2 -regularization in Problem 1 is replaced by a term γ∥x∥1 in the objective. This formulation can be equivalently expressed as a linear program by introducing auxiliary variables, and is useful for enforcing stability across consecutive solver runs. All experiments run on a single node with NVIDIA H100 80 GB GPUs (either 4× or 8× on the same host). The solver is invoked in single-node multi-GPU mode with AGD updates, a six-stage geometric γ-schedule γ ∈ {103 , 102 , 10, 1, 10−1 , 10−2 } (10,000 iterations per stage; 60,000 total), Jacobi preconditioning, and AGD step-size range [10−5 , 10−1 ]. D-PDLP is run on the

9

Figure 2: Geometric bucketing vs. single-slab baseline (Triton kernel, single H100). (Top) p95 per-iteration time; (Bottom) peak GPU memory. same instance with presolve disabled and a 10−4 relative tolerance on the optimality and feasibility residuals. Presolve is disabled to allow more memory headroom for D-PDLP solver. Dualip-DPU materializes the LP directly on the GPU; D-PDLP first loads a pre-generated MPS file the size of which grows from 63 GB (s50M base) to 306 GB (s100M ℓ1 -reformulated). Table 3 reports end-to-end solve time on the eight problem instances and two GPU counts. Two observations stand out. First, Dualip completes every base instance up to 109 nonzeros at both 4 and 8 GPUs, while D-PDLP runs out of memory on the largest base instance (s100M, 1.0 ×109 nonzeros) and on every ℓ1 -reformulated instance, including the smallest (∼ 109 nonzeros). Second, the runtime growth on the instances Dualip solves is near-linear in problem size and in hardware: doubling the nonzero count from 0.50B to 1.00B increases 8-GPU solve time by 1.9×, and doubling the GPU count from 4 to 8 yields a 1.88–1.92× speedup across the four base sizes. Where D-PDLP succeeds it is markedly faster per instance, but this efficiency does not translate into a usable solver at the scales required by the ℓ1 reformulation, which is the operating point of interest for our application. Note that the GPU-based DuaLip solver can solve such instances within a 2-hour time budget, making it practical for production use cases where this level of latency is acceptable. Table 4 compares the solutions returned by the two solvers at 8 GPUs on the instances where D-PDLP completes, together with the largest instance (s100M, 1.0×109 nonzeros) which only Dualip solves. D-PDLP terminates adaptively once both the relative primal and dual residuals fall below 10−4 . Dualip runs the full 60,000-iteration schedule described above; we report the primal and dual objectives at the end of the schedule, the absolute primal–dual gap, and the maximum positive constraint violation (“slack”) of the 10

Figure 3: Scaling from 1–16 GPUs (16-GPU spans 2 H100 nodes). (Top) Solve time; (Bottom) speedup vs. 1 GPU. recovered primal. On instances where both solvers succeed, the dual objectives agree to four significant figures (∆dual/|dual| < 10−6 across all three problem sizes), indicating that both methods converge to the same optimum when the regularization parameter is sufficiently small (0.01 in the final stage). This result is important because it shows that the regularized formulation provides a controllable trade-off: the regularization strength can be tuned to adjust solution fidelity while also improving numerical stability when needed. Additionally, our system achieves a substantially smaller primal–dual gap (down to 10−12 ) compared to D-PDLP, reflecting the improved conditioning of the smoothed problem formulation.

7.3

Effect of Algorithmic Enhancements

We now isolate the impact of two algorithmic improvements introduced in Section 6: diagonal preconditioning and regularization continuation. Preconditioning. Figure 4 plots log(|L − L̂|), where L is the dual objective and L̂ is the converged reference value. With diagonal preconditioning, convergence accelerates substantially, particularly in early iterations. This confirms that scaling the dual variables mitigates ill-conditioning arising from heterogeneous constraint magnitudes.

11

Table 3: Solver runtime on synthetic statistical-matching LPs at scale. Each cell is solve time in seconds on N NVIDIA H100 80GB GPUs (single node). “OOM” = the solver exhausted GPU/CPU memory and produced no output. “–” = configuration not run. Dualip generates the LP on the fly; D-PDLP loads pre-built MPS files (135–327 GB on disk for the largest cases). n×m

Instance

Dualip (s)

NNZ

D-PDLP (s)

4 GPU

8 GPU

4 GPU

8 GPU

2384 3522 3738 4619

1269 1841 1944 2441

6.8 10.2 10.9 OOM

4.5 6.6 7.0 OOM

ℓ1 -reformulated instances (≈ 2× NNZ of base) s50M d10K + L1 50M × 10K ∼1.0B NA s75M d10K + L1 75M × 10K ∼1.5B NA s80M d10K + L1 80M × 10K ∼1.6B NA s100M d10K + L1 100M × 10K ∼2.0B NA

NA NA NA NA

OOM OOM OOM OOM

OOM OOM OOM OOM

Base instances s50M d10K s75M d10K s80M d10K s100M d10K

50M × 10K 75M × 10K 80M × 10K 100M × 10K

0.50B 0.75B 0.80B 1.00B

Table 4: Solution quality at 8 GPUs. D-PDLP terminates when relative primal and dual residuals fall below 10−4 . Dualip runs 60,000 iterations with a six-stage γ schedule (γ ∈ {103 , 102 , 10, 1, 10−1 , 10−2 }, 10,000 iterations per stage). All experiments completed within 2 hours time limit set for both solvers. “Gap” denotes relative objective gap. Slack is the maximum constraint violation by each solver. Instance

Solver

Primal

Dual

Gap

Constraint Slack

s50M–d10K

D-PDLP Dualip

−2.510 × 106 −2.444 × 106

−2.511 × 106 −2.444 × 106

7.2 × 10−5 −6.9 × 10−7

4.6 × 10−7 † 5.3 × 10−2

s75M–d10K

D-PDLP Dualip

−3.829 × 106 −3.724 × 106

−3.829 × 106 −3.724 × 106

7.9 × 10−5 3.5 × 10−10

4.4 × 10−7 † 1.2 × 10−5

s80M–d10K

D-PDLP Dualip

−4.095 × 106 −3.983 × 106

−4.096 × 106 −3.983 × 106

8.0 × 10−5 5.6 × 10−12

4.1 × 10−7 † 2.7 × 10−6

s100M–d10K

D-PDLP Dualip

−4.879 × 106

OOM −5.012 × 106 1.6 × 10−6

1.4 × 10−1

† D-PDLP reports a relative primal residual instead of max constraint slack.

Regularization Continuation. Figure 5 evaluates the continuation strategy for the ridge parameter γ. Starting from a larger γ stabilizes and accelerates early optimization, while gradual decay ensures the final solution closely approximates the unregularized LP optimum. Decaying γ from 0.16 to 0.01 (halved every 25 iterations) yields faster convergence compared to using a fixed regularization level.

8

Discussion and Limitations

In this paper, we described the main architectural, algorithmic, and systems design choices for extreme-scale matching LPs in production systems. A natural next step is to evaluate the same dual-ascent framework on broader classes of linear programming benchmarks beyond matching. A key limitation is the lack of publicly available datasets at the scale and structure of production matching workloads. As a result, evaluation is based primarily on large synthetic instances calibrated to observed industrial distributions. Developing standardized benchmark suites for extreme-scale matching LPs would significantly improve reproducibility

12

Figure 4: Effect of diagonal preconditioning. We report log(|L− L̂|) for a 25M-source instance (10k destinations, 0.1% sparsity). Preconditioning significantly improves early-stage convergence.

Figure 5: Effect of regularization continuation. Decaying γ during optimization accelerates convergence while preserving solution fidelity. and enable more systematic comparison of methods in this regime.

13

A

Experiment setting

Synthetic LP construction. We construct synthetic instances by first generating a sparse bipartite interaction graph and then assigning values and constraint coefficients on its edges. Given a target number of “requests” I, “resources” J, and a target sparsity level, we draw a lognormal “breadth” parameter for each resource j, normalize these to obtain probabilities pj , and sample the number of incident requests Kj ∼ Poisson(pj Iν), truncated at I, where ν is the desired average number of nonzeros per row. For each resource j, we then select Kj distinct requests and create edges (i, j). On each edge, we draw a resourcespecific value scale vj , a request-specific responsiveness factor ui , and multiplicative noise εij , and define a nonnegative value coefficient
cij = min vj ui εij , cmax . Constraint coefficients are taken to be scaled versions of these values, aij = sj cij , where the per-resource scale sj is also drawn from a lognormal distribution. This construction yields a sparse matrix A whose rows differ both in support size and magnitude (often by several orders), and a matching value matrix C with the same sparsity pattern. Source capacities and right-hand side. Right-hand side source capacities bjPare chosen to make a nontrivial fraction of constraints active. Instead of taking bj proportional to the sum i aij , we approximate the maximum feasible load each resource could receive under the per-request simplex constraint (each request can allocate to at most one resource) by a greedy assignment: for each request i, we identify the incident edge with the largest aij and assign that amount to the corresponding resource. Summing these contributions over requests gives a “greedy load” ℓj for each resource, and we set
bj = ρj ℓj + ε , where ρj is drawn uniformly from [0.5, 1.0] and ε > 0 is a small constant. This ensures that some resource constraints are binding while others remain slack in the optimal solution. The final LP data passed to the solver is the sparse CSC representation of A, the corresponding value matrix (with signs adjusted to match our minimization convention), and the source capacity vector b.

B

Additional Experiments

B.1

Implementation parity.

We first verify numerical equivalence between the PyTorch implementation and the original Scala solver under identical problem instances and hyperparameters. Both systems are initialized from the same primal-dual state and execute the same accelerated gradient descent (AGD) update schedule, with differences limited to execution backend (CPU vs. GPU) and memory layout. Figure 6 compares dual objective trajectories across iterations in single-GPU and multi-GPU configurations. Across all tested problem sizes and constraint families, the trajectories closely overlap, indicating that distributed execution and kernel fusion do not alter the optimization dynamics. To quantify discrepancies, Figure 7 reports the relative error in the dual objective with respect to the Scala solver, computed as |gtorch (λt ) − gscala (λt )| . |gscala (λt )| In all configurations, the relative error drops below 1% within the first 100 iterations and continues to decay as the iterates converge, consistent with stabilization of the primal-dual dynamics. We further observe that error behavior is stable across single-GPU and multi-GPU runs, indicating that NCCL-based aggregation and column sharding preserve numerical consistency. These results confirm that the PyTorch implementation faithfully reproduces the optimization trajectory of the production Scala system while supporting distributed execution and GPU acceleration. 14

Figure 6: Scala–DuaLip (PyTorch) parity. Each panel shows the dual objective versus AGD iteration for the Scala and PyTorch implementations. The near-perfect overlap confirms numerical equivalence.

B.2

Algorithm enhancements

A central area of focus in improving the DuaLip solver was to improve the performance and robustness of the ridge-regularized dual ascent method that underpins this line of work. The basic framework of Section 3.1 admits many first-order instantiations; in practice, the dominant convergence issues we observed on production matching workloads stemmed from three sources: (i) poor conditioning of AA⊤ in the dual, (ii) tradeoff between convergence speed and solution quality with the ridge regularization term and (iii) heterogeneous scales in the primal variables that interact badly with the quadratic regularizer. We address these in turn via Jacobi-style row normalization and a regularization schedule. Jacobi preconditioning / row normalization. Our first enhancement is to improve the conditioning of the dual problem by rescaling the complex constraints. We assume A is full row rank. Intuitively, when some rows of A have much larger norms than others, gradient steps for the smoothed dual move too cautiously in some directions and too aggressively in others. When the projection operator is inactive (or the identity), the dual gradient reduces to ∇g(λ) = − so the Hessian is

1 AA⊤ λ + Ac − b, γ

1 ∇2 g(λ) = − AA⊤ . γ

The convergence of first-order methods on g is therefore governed by the conditioning of AA⊤ .

15

Figure 7: Relative error in dual objective compared to the Scala solver. The error drops below 1% within the first 100 iterations across all settings. We apply a standard row-scaling transform. Let −1 D = diag ∥A1∗ ∥−1 2 , . . . , ∥Am∗ ∥2

for all rows with nonzero norm (rows with ∥Ar∗ ∥2 = 0 are redundant and may be dropped or left unscaled with Drr = 1), and define the row-scaled system A′ = DA,

b′ = Db.

Because D has positive diagonal entries, row scaling preserves the feasible set exactly: { x : Ax ≤ b } = { x : A′ x ≤ b′ }. Moreover, A′ A′⊤ = D(AA⊤ )D,

D2 = diag(AA⊤ )−1 (on nonzero rows),

so row normalization is precisely Jacobi preconditioning of the dual Hessian −∇2 g(λ) = γ1 AA⊤ . For the matching constraint matrix of Definition 1, A is a horizontal concatenation of diagonal subblocks across sources, so AA⊤ is a sum of (nearly) diagonal matrices (one per source) and is therefore close to diagonal in practice. Enforcing diag(A′ A′⊤ ) = I tightly clusters the spectrum; in the ideal diagonal/orthogonal case it yields A′ A′⊤ = I and condition number 1. The following lemma formalizes this intuition under a simple statistical model for the matching blocks.

16

Lemma B.1. Let A = [A1 · · · AI ] ∈ RmJ×IJ with user blocks Ai ∈ RmJ×J that are i.i.d. across i and
−1/2 e = Dexp A. Then diagonal by rows as in Definition 1. Let Dexp = diag E∥A1∗ ∥22 , . . . , E∥Am2 ∗ ∥22 and A
⊤ eA e ] = I. If, in addition, for r ̸= s, diag E[A q   E ⟨Ar∗ , As∗ ⟩ ≤ η E∥Ar∗ ∥22 E∥As∗ ∥22 then

with η ∈ [0, 1),

eA e ⊤ ] ≤ 1 + (m − 1)η . κ E[A 1 − (m − 1)η

Thus, under mild cross-row correlation assumptions, row normalization nearly equalizes the eigenvalues of AA⊤ in expectation for matching workloads, which stabilizes dual first-order updates. Regularization decay The preconditioning above addresses the dominant factor when the simple constraints are inactive: in that regime, the condition number of AA⊤ controls the convergence rate of first-order methods. However, the ridge parameter γ enters separately through the smoothed dual. Next, we consider how to set the regularization hyperparameter over the course of optimization. While γ does not affect the conditioning of AA⊤ , it has a substantial effect on the smoothness of the objective. On the one hand, one prefers small γ to avoid perturbing the original LP. Lemma 2 in Basu et al. [3] guarantees that there is always a range of sufficiently small γ values that allow exact recovery of a solution to the un-smoothed LP. On the other hand, the Lipschitz constant of the gradient is only bounded by ∥A∥22 /γ (see Lemma 3 in [3]), so very small γ leads to a poorly conditioned dual problem and slow progress in practice. To balance convergence speed and solution fidelity, we implement a simple continuation scheme in which γ is initialized at a moderately large value (for stable, fast early progress) and then decayed on a pre-specified schedule and rate as the algorithm approaches a good dual solution. Since γ directly affects the smoothness of the dual objective, we scale the maximum AGD step size proportionally with the decay of γ to maintain stability across transition points. Overall, this approach allows for faster convergence in early iterations, while ensuring the final solution is not greatly perturbed from the true LP solution.

17

References [1] Ayan Acharya, Siyuan Gao, Borja Ocejo, Kinjal Basu, Ankan Saha, Keerthi Selvaraj, Rahul Mazumdar, Parag Agrawal, and Aman Gupta. 2023. Promoting inactive members in edge-building marketplace. In Companion Proceedings of the ACM Web Conference 2023. 945–949. [2] David Applegate, Mateo Díaz, Oliver Hinder, Haihao Lu, Miles Lubin, Brendan O’Donoghue, and Warren Schudy. 2021. Practical large-scale linear programming using primal-dual hybrid gradient. Advances in Neural Information Processing Systems 34 (2021), 20243–20257. [3] Kinjal Basu, Amol Ghoting, Rahul Mazumder, and Yao Pan. 2020. ECLIPSE: An extreme-scale linear program solver for web-applications. In International Conference on Machine Learning. PMLR, 704–714. [4] Niv Buchbinder, Moran Feldman, Arpita Ghosh, and Joseph Naor. 2014. Frequency capping in online advertising. Journal of Scheduling 17, 4 (2014), 385–398. [5] John Duchi, Shai Shalev-Shwartz, Yoram Singer, and Tushar Chandra. 2008. Efficient projections onto the l 1-ball for learning in high dimensions. In Proceedings of the 25th international conference on Machine learning. 272–279. [6] Lisa Fleischer, Michel X Goemans, Vahab S Mirrokni, and Maxim Sviridenko. 2006. Tight approximation algorithms for maximum general assignment problems. In SODA, Vol. 6. 611–620. [7] Hongpei Li, Yicheng Huang, Huikang Liu, Dongdong Ge, and Yinyu Ye. 2026. D-PDLP: Scaling PDLP to Distributed Multi-GPU Systems. arXiv:2601.07628 [math.OC] https://arxiv.org/abs/2601.07628 [8] Elita Lobo, Justin Payan, Cyrus Cousins, and Yair Zick. 2024. Fair and welfare-efficient constrained multi-matchings under uncertainty. Advances in Neural Information Processing Systems 37 (2024), 74579–74616. [9] Haihao Lu, Zedong Peng, and Jinwen Yang. 2025. cuPDLPx: A Further Enhanced GPU-Based First-Order Solver for Linear Programming. arXiv:2507.14051 [math.OC] https://arxiv.org/abs/2507.14051 [10] Haihao Lu and Jinwen Yang. 2024. cuPDLP.jl: A GPU Implementation of Restarted Primal-Dual Hybrid Gradient for Linear Programming in Julia. arXiv:2311.12180 [math.OC] https://arxiv.org/ abs/2311.12180 [11] Haihao Lu and Jinwen Yang. 2025. cuPDLP.jl: A GPU Implementation of Restarted Primal–Dual Hybrid Gradient for Linear Programming in Julia. Operations Research 73, 6 (2025), 3440–3452. [12] Joseph Naor and David Wajc. 2018. Near-optimum online ad allocation for targeted advertising. ACM Transactions on Economics and Computation (TEAC) 6, 3-4 (2018), 1–20. [13] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. 2019. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems 32 (2019). [14] Rohan Ramanath, S Sathiya Keerthi, Yao Pan, Konstantin Salomatin, and Kinjal Basu. 2021. Efficient algorithms for global inference in internet marketplaces. arXiv preprint arXiv:2103.05277 (2021). [15] Philippe Tillet. 2025. Introducing Triton: Open-source GPU programming for neural networks. [16] Panagiotis Tsiakis and Lazaros G Papageorgiou. 2008. Optimal production allocation and distribution supply chain networks. International Journal of Production Economics 111, 2 (2008), 468–483. [17] Vijay V Vazirani. 2001. Approximation algorithms. Vol. 1. Springer.

18