PaperlayerView exact revision on arXiv

Source-linked AI summary

Distributed Linear Programming on GPU Clusters at Extreme Scale

Arnaud Deza, Santanu Dey, Pascal Van Hentenryck

arXiv:2609.09108v1math.OCcs.DC

TL;DR

Large LPs may exceed a single node’s memory, and solver phases beyond the matrix products can restore that limit. SHARDLP keeps the matrix and primal–dual state partitioned throughout the solve; it meets the published criterion on nine of eleven Google PDLP benchmarks and validates larger multi-node solves.

  • Problem

    Although PDHG’s GPU-suited operations distribute naturally, scaling a complete sparse LP solver across nodes remains less developed, especially when the matrix and primal–dual state cannot fit on one node.

    Source: §1 · p. 1

  • Method

    SHARDLP keeps the matrix and primal–dual state partitioned through solver stages, from prepartitioned input shards through distributed recovery and output.

    Source: §1 · p. 2, §2.2 · p. 3

  • Results

    Nine of eleven Google PDLP benchmark instances meet the published criterion, versus eight in the published CPU PDLP study.

    Source: §1 · p. 2, §3.1 · p. 4

  • Takeaways & Limitations

    The results identify memory ownership and communication as distinct scale bottlenecks, with support-aware communication reducing modeled data movement on Design Match.

    Source: §5 · p. 6

    The LP must already be partitioned into shards, and support-aware communication is implemented only for column partitions with fixed matrix support.

    Source: §4 · p. 6

Abstract

from arXiv · show

Large linear programs can exceed the memory of a single compute node. Although first-order methods replace sparse factorizations with GPU-suited matrix-vector products, other solver phases can reintroduce a single-node memory limit. We present SHARDLP, a distributed GPU LP solver that keeps the matrix and primal-dual state partitioned from sharded input through solution output. On the Google PDLP benchmark, SHARDLP reaches the published criterion on nine of eleven instances, compared with eight in the published CPU PDLP study. On the largest benchmark, eight H200 GPUs solve a 1.185-billion-variable, 6.338-billion-nonzero LP in 9.9 minutes; the published CPU experiment reports 21.06 hours on different hardware. Beyond this benchmark, separately checked multi-node solves reach up to 13.604 billion variables and 40.807 billion nonzeros, while validated executions span up to 76 GPUs across 29 compute nodes. For column-partitioned solves, support-aware communication skips GPUs that store no coefficients for a row; on an LP with 2.76 billion nonzeros, it cuts modelled communication by 92.97% and improves solver time by 1.27x-1.52x

1. Introduction

Large LPs challenge factorization-based methods with memory-intensive linear algebra, while PDHG offers GPU-suited sparse operations. SHARDLP addresses the gap between distributing PDHG iterations and building a complete solver for problems whose matrix and primal–dual state cannot fit on one node.

Motivation

  • PDHG-based first-order methods use sparse matrix–vector products, projections, and vector operations that map naturally to GPUs and distributed memory.By contrast, simplex and interior-point methods can become memory intensive and difficult to parallelize at very large scale because of factorization-based linear algebra.Source: §1 · p. 1

Background

  • PDHG, also known as the Chambolle–Pock method, underlies PDLP, which adds diagonal scaling, presolve, adaptive step sizes, restart, and feasibility polishing.The principal benchmark uses eleven released instances from the large-scale Google PDLP study; its largest has 1.185 billion variables and 6.338 billion nonzeros.Source: §1 · p. 1

Research gap

  • Recent GPU work demonstrates that PDHG maps well to accelerators, but scaling a general sparse LP solver across multiple compute nodes remains less developed.D-PDLP distributes matrix products across eight H100 GPUs within one node, while MPAX leaves efficient distributed sparse-data sharding as future work.Source: §1 · pp. 1–2

Contributions

  • SHARDLP reaches the published criterion on nine of eleven Google PDLP benchmarks and separately validates multi-node solves up to 13.604 billion variables and 40.807 billion nonzeros.Its support-aware communication cuts modeled communication by 92.97% and improves solver time by 1.27×–1.52× on Design Match.Source: §1 · p. 2

2. From a distributed iteration to a distributed solver

SHARDLP builds on PDHG and a two-dimensional matrix decomposition, then preserves distributed ownership through the rest of the solve. This avoids requiring any solver phase to gather the full problem on one node.

2.1. PDHG and the D-PDLP matrix decomposition

  • PDHG alternates a primal update using A^T y and a dual update using Ax; SHARDLP follows cuPDLPx optimization logic and uses a reflected-Halpern update within each restart epoch.Scaling, step selection, primal-weight adaptation, and restart follow Lu et al. [22], with control logic using scalar reductions.Source: §2.1 · p. 2
  • Rank (r,c) stores A_rc in an R×C process grid; row and column reductions form Ax and A^T y while ownership remains sharded through input and output.The grid is chosen before sharding so estimated matrix blocks and replicated vector slices fit each GPU; per-rank storage is O(nnz(A_rc) + n/C + m/R).Source: §2.1 · pp. 2–3
  • The two-dimensional decomposition follows D-PDLP [19]; SHARDLP’s contribution is preserving that ownership beyond the matrix products and throughout the solve.The figure’s panel (b) depicts row/column reductions while ownership remains sharded; panel (a) concerns centralized solver phases.Source: §2.1 · p. 3, Figure 1 · p. 3

2.2. Persistent ownership beyond the iteration

  • Each rank reads a prepartitioned shard, and scaling, restart, termination, polishing, recovery, and output remain distributed; converting a monolithic model into shards is offline.General presolve is disabled; a distributed pass fixes variables forced to zero by globally singleton zero-equality rows, removing 1.123 billion nonzeros (40.69%) in both RQ3 Design Match arms, but not in RQ1.Source: §2.2 · p. 3

3. Computational Results

The experiments test SHARDLP against the Google PDLP benchmark, at larger scales, and under multi-GPU communication. They report benchmark results, separately checked large solves, and cases where communication limits solver-loop speed or support-aware communication helps.

3. Computational Results

  • For solves beyond RQ1, a separate checker requires nine original-space quantities to be finite and at most 10^-6; timing excludes offline shard construction and the final check.RQ1 uses the Google study’s published scaled-LP criterion.Source: §3 · pp. 3–4

3.1. RQ1: Google PDLP benchmark

  • SHARDLP reaches the published PDLP criterion on 9 of 11 benchmark instances, versus 8 in the published CPU PDLP study.The study also reports 592.6 seconds of solver time and 933.4 seconds end to end for TSP-Gaia-100M on eight H200 GPUs; the published 32-core CPU experiment reports 21.06 hours.Source: §3.1 · p. 4

3.2. RQ2: Scaling beyond the Google benchmark

  • 13.604 billion variables and 40.807 billion nonzeros were separately validated on 76 GPUs across 29 nodes; that solve took 1,905 seconds.The larger-scale experiments also include validated solves with 25.000 billion nonzeros on 56 GPUs, and MCF5B was validated in all three runs on 40 GPUs across ten nodes.Source: §3.2 · p. 4, Table 2 · p. 5

3.3. RQ3: Does communication limit multi-GPU scaling?

  • On Mediterranean Shipping, 4,000 iterations took 107.2 seconds on one Blackwell GPU and 112.7 seconds on four, while end-to-end time fell from 281.5 to 168.1 seconds.The fixed-iteration speedup relative to one GPU is the Figure 2a comparison.Source: §3.3 · p. 5, Figure 2 · p. 5
  • On Design Match, support-aware communication cut modeled communication by 92.97% and was 1.269× faster on V100s and 1.518× faster on Quadro RTX 6000s.Figure 2b reports the geometric-mean speedup of support-aware over dense communication.Source: §3.3 · pp. 5–6, Figure 2 · p. 5
  • On QAP-WIL, the modeled communication reduction was 38.35% and the speedup was 1.024×; support-aware communication skips GPUs with no nonzeros in a row.Because the skipped GPUs contribute zero, the exact-arithmetic PDHG update is unchanged.Source: §3.3 · pp. 5–6

4. Limitations and Future Work

SHARDLP assumes prepared input shards, and support-aware communication is currently limited to column partitions with fixed matrix support. Future work targets larger problems, reduced communication and memory use, broader distributed solver stages, and applications to other optimization problems.

4. Limitations and Future Work

  • Monolithic-model conversion to SHARDLP shards is currently an offline preprocessing step, so the solver assumes the LP is already partitioned.The paper identifies this input-format requirement as a limitation.Source: §4 · p. 6
  • Support-aware communication is implemented only for column partitions with fixed matrix support.The authors describe these choices as sufficient for the paper’s experiments, while leaving room for a more general distributed solver pipeline.Source: §4 · p. 6
  • Future work includes better partitioning, broader distributed presolve and postsolve, and extending support-aware communication to more general process grids.The authors also identify quadratic and semidefinite programs and large discrete-optimization algorithms as possible applications.Source: §4 · p. 6

5. Conclusion

SHARDLP keeps the problem and solution distributed throughout the solve, addressing both memory ownership and communication as scale bottlenecks.

5. Conclusion

  • Persistent ownership keeps the matrix and primal-dual state distributed throughout every solver stage, with no phase requiring a full copy on one node.The paper reports that this approach meets the published criterion on nine of eleven Google benchmarks and separately checked LPs with up to 13.604 billion variables and 40.807 billion nonzeros.Source: §5 · p. 6
  • Support-aware communication follows matrix sparsity by skipping GPUs with no coefficients for a row; on Design Match it cuts modeled data movement by 92.97%.The reported solver-time speedups are 1.27×–1.52×.Source: §5 · p. 6
  • The results identify memory ownership and communication as distinct scale bottlenecks; smaller per-GPU blocks alone do not ensure speed.Source: §5 · p. 6

Appendix A. Computational environment and hardware

The experiments ran on Georgia Tech’s PACE Phoenix cluster with one MPI process per GPU, using multiple node and GPU configurations. Reported device-memory peaks are sampled measurements, not exact allocator high-water marks.

Appendix A. Computational environment and hardware

  • Experiments ran on PACE Phoenix under Slurm with one MPI process per GPU; solver GPU counts could be lower than the number installed in an allocated node.Jobs reserved GPUs, CPU cores, and host memory without assuming whole-node exclusivity.Source: Appendix A · p. 9
  • GPU communication used NCCL 2.26.5, with MPI managing the process grid and control reductions; NVLink connected H200 GPUs, while other listed GPU nodes used PCIe paths.The distributed C/CUDA runs used CUDA 12.9.1 and CUDA-aware OpenMPI 4.1.8.Source: Appendix A · p. 9
  • Inter-node transport varied by run: several benchmarks used NCCL’s socket path, while the two-node QAP-THO-150 run enabled InfiniBand.These settings do not establish a common physical-link bandwidth or cluster-wide RDMA performance.Source: Appendix A · p. 9
  • Table 4 reports physical storage and device-memory use; peak GPU memory is the largest periodic nvidia-smi sample matched to participating GPU UUIDs, not an exact allocator high-water mark.Input and primal–dual output sizes are totals across shards.Source: Appendix A · p. 10, Table 4 · p. 10

Appendix B. Instance origins

The appendix summarizes the dimensions and construction of the Google benchmark instances and distinguishes generated application-inspired models from instances built from named public data sources.

Appendix B. Instance origins

  • Table 5 summarizes the dimensions and construction of the eleven Google PDLP benchmark instances; m, n, and nnz(A) denote constraints, variables, and explicitly stored matrix coefficients.The benchmark instances were released in the large-scale Google PDLP study [5].Source: Appendix B · p. 10, Table 5 · p. 11
  • Synthetic instances are generated, application-inspired models, while LPs built from named public data sources are identified separately.Generators and construction notes are available for a subset with released construction code.Source: Appendix B · p. 10
  • MCF8.35B uses 8,303 commodities and 1,000 factories per commodity; MCF13.60B uses 34 commodities and 20,000 factories per commodity, and has more nonzeros but fewer rows, reducing row-vector communication.Both instances use 1,000 or 20,000 warehouses, respectively, and five stores.Source: Appendix B · p. 10
  • Table 6 summarizes the origins and validation status of the beyond-benchmark experiments.Source: Table 6 · p. 12

C.1. Data layout across solver stages

SHARDLP begins from prepared matrix shards and retains distributed ownership through the solver lifecycle, recovery, and output; its prepared partition remains fixed during a solve.

C.1. Data layout across solver stages

  • Each process starts from a prepared shard; converting a monolithic input into shards happens offline and is excluded from reported times.In an R × C grid, each rank holds its matrix block and vectors replicated only within its grid column or row.Source: §C.1 · p. 13
  • Table 7 describes the distributed lifecycle, where a global reduction is scalar or partition-sized and never gathers the complete matrix or primal–dual state.The same partition is retained through the solve.Source: Table 7 · p. 13, §C.1 · p. 13
  • Singleton-pivot recovery updates the primal variable, row dual, and reduced cost in original units; infeasibility and unboundedness rays are not implemented.Prepared blocks use contiguous row and column intervals in a fixed coordinate order.Source: §C.1 · p. 13

C.2. Validation criteria and timing

The section specifies scaled-LP acceptance criteria, an independent original-coordinate check, and timing conventions that distinguish solver time from end-to-end time.

C.2. Validation criteria and timing

  • The published rule requires scaled-LP absolute primal and stationarity residuals at most 10^-8, a relative primal–dual gap at most 10^-2, and bound and sign membership.Finite-precision membership tests use a 10^-8 tolerance [5].Source: §C.2 · p. 13
  • Sign-admissible projected dual quantities are used for stationarity and dual-objective evaluation, while distances from the reported signs are checked separately.The interval-based multiplier signs depend on finite lower or upper bounds, and each interval contributes to the dual objective.Source: §C.2 · p. 14
  • Nine criteria cover variable-bound feasibility, row feasibility, stationarity, dual and reduced-cost signs, and the relative primal–dual objective gap.For each paired feasibility or stationarity test, one quantity is normwise and the other is the worst coordinatewise normalized residual; Table 8 prints three tests and the maximum over all nine.Source: §C.2 · pp. 14–15
  • The separate checker evaluates original coordinates with long-double accumulation; solved runs use a 10^-6 check for primal feasibility, stationarity, admissible signs, and primal–dual objective agreement.Source: §C.2 · pp. 13, 15
  • Solver time includes solver iterations, residual evaluations, restart, and enabled polishing; end-to-end time also includes input, initialization, recovery, and output, but excludes offline conversion and the separate final check.Source: §C.2 · p. 15
  • The Google study used a six-day limit [5], so Table 1 compares coverage and reports raw runtimes rather than a controlled hardware comparison between the published systems.The published CPU experiments used either a 16-core AMD EPYC 7302 node with 256 GB of memory or a 32-core Intel Xeon Platinum 8352Y node with 1 TB.Source: §C.2 · p. 15

Appendix D. Support-aware update equivalence and communication model

The support-aware path assigns each row to only the column shards containing its nonzeros, while reconstructing required dual state at synchronization points. Its reported communication reduction is modeled, and the equivalence argument assumes fixed support and matching solver rules.

Appendix D. Support-aware update equivalence and communication model

  • Each row’s participant set contains exactly the GPUs whose column shards store a nonzero coefficient for that row.The dense path stores updated dual coordinates on every rank; the participant path returns updates only to ranks whose transpose product uses that coordinate.Source: Appendix D · p. 16, Figure 3 · p. 16
  • The communication plan is built once after support-changing preprocessing, and structural changes or repartitioning require rebuilding it.Rows sharing a participant set are grouped; full-participation buckets use NCCL Reduce/Broadcast, while partial buckets use point-to-point transfers through the lowest participating rank.Source: Appendix D · p. 16
  • Whole-solver equality at synchronization points is proved under exact arithmetic, identical solver rules, fixed exact row support, and authoritative dual-coordinate ownership.Required dual coordinates are reconstructed before operations needing replicated state, while global statistics count each authoritative coordinate once.Source: Appendix D · p. 16
  • At shared V100 Design Match checkpoints, the largest reported dense-versus-participant differences range from 1.30 × 10−15 for relative gap to 1.98 × 10−11 for relative stationarity.The reported maximum for relative primal residual is 1.87 × 10−12.Source: Appendix D · p. 17
  • The communication model counts logical scalar-hops and excludes packing, imbalance, topology, and latency, so its reductions need not predict elapsed time.It counts reduction and dissemination legs plus one-way full-state synchronization boundaries.Source: Appendix D · p. 17
  • 92.9748% lower modeled communication is reported for Design Match through solver-specific integration of support-aware updates.The contribution combines a fixed producer/consumer plan with the nonlinear update, empty-row ownership, solver state, polishing, checking, and output boundaries.Source: Appendix D · p. 17

Appendix E. Related work and comparison scope

The comparison focuses on general sparse GPU LP solvers and counts runs that export a complete primal–dual solution passing an independent checker. SHARDLP’s scale claim is not a matched-speedup comparison, and one larger capacity run was not independently checked.

Appendix E. Related work and comparison scope

  • The comparison includes general sharded sparse-matrix GPU LP solvers only when runs export complete primal–dual solutions that pass an independent checker.The 70-GPU capacity experiment did not export complete vectors for a separate primal–dual check, so its runtime is not a checked solve or time-to-solution result.Source: Appendix E · pp. 17–18
  • D-PDLP distributes two PDHG matrix products across eight H100 GPUs within one node, while MPAX reports dense multi-GPU PDHG and leaves efficient distributed sparse-data sharding as future work.D-PDLP’s largest-variable instance has 127.5 million variables and 257.5 million nonzeros; MPAX reports approximately 900 million represented nonzeros on up to four H100 GPUs.Source: Appendix E · p. 17
  • Other systems include a matching-specialized solver spanning sixteen H100 GPUs across two nodes and ECLIPSE’s distributed matrix–vector products for a structured web LP.The matching solver reports models up to one billion nonzeros; ECLIPSE optimizes the smooth dual of a ridge-perturbed model.Source: Appendix E · p. 17
  • MCF13.60B is reported as the largest validated solve among the discussed general sparse GPU LP solvers, with 40.807 billion nonzeros.The authors explicitly distinguish this scale claim from a claim of being first or demonstrating matched speedup over the compared systems.Source: Appendix E · pp. 17–18
Loading 2609.09108v1…