Scaling Linear Programming to 100 Million Variables with NVIDIA cuOpt mPDLP

Supply-chain models now cover more SKUs, transportation lanes, and constraints, while energy grids must balance increasing numbers of distributed sources in real time. These workloads require larger optimization models and more uncertainty to be evaluated within practical planning windows.

NVIDIA cuOpt provides GPU-accelerated decision optimization and can deliver speedups of more than 10x over CPU solvers on a single GPU for large-scale linear programming (LP) problems. However, the largest planning models can still take hours to converge or exceed the memory capacity of one GPU.

A long-standing benchmark illustrates the challenge. The zib03 problem, introduced by Thorsten Koch in 2008, contains more than 104 million nonzeros and has been used to evaluate LP algorithms and solvers. New algorithms and hardware have steadily reduced its solve time. NVIDIA cuOpt now addresses the problem across multiple GPUs with its Multi-GPU Primal-Dual hybrid gradient for Linear Programming, or mPDLP.

What cuOpt mPDLP provides

The mPDLP solver distributes LP problems across NVLink-connected GPUs. This approach is designed to provide two main benefits:

• Reduced solve times for problems that are impractical to complete within a planner's available time window
• Up to 6x lower peak memory usage per GPU than single-GPU PDLP, for LP problems capped at 2.1 billion nonzeros

The approach has also been applied by Kinaxis to globally constrained supply and production planning, and by PSR to large-scale energy-system capacity-expansion models.

How the PDLP algorithm solves LPs

An LP seeks a vector x that solves the following problem:

minimize    cᵀx
subject to  Ax = b
            l ≤ x ≤ u

Primal-Dual hybrid gradient for LP, or PDLP, is a first-order method based on gradient-descent-like operations. Its structure makes it highly parallelizable and suitable for GPUs.

The algorithm's main computational loop is a sparse matrix-vector multiplication (SpMV), followed by several element-wise operations. It updates two vectors:

• The primal solution, x
• The dual solution, y

The LP constraints are represented by the sparse matrix A. PDLP uses both A and its transpose, Aᵀ, in each iteration:

xₖ₊₁ = proj[l,u](https://developer.nvidia.com/blog/scaling-decision-optimization-to-100-million-variables-and-beyond-with-mpdlp-in-nvidia-cuopt/x%E2%82%96 - τ(c - Aᵀyₖ))
yₖ₊₁ = yₖ + σ(b - A(2xₖ₊₁ - xₖ))

Here, τ and σ are the primal and dual step lengths. The projection function enforces the variable bounds l ≤ x ≤ u. The vector b contains the right-hand sides of the equality constraints, while c contains the objective-function costs.

Most of these operations are element-wise projections or basic arithmetic, which are straightforward to distribute across GPUs. The significant communication steps are the two SpMVs:

xₖ₊₁ = Aᵀyₖ
yₖ₊₁ = Axₖ₊₁

The updated y is used to compute the next x through one SpMV, and the updated x is used to compute the next y through another.

Distributing PDLP across GPUs

SpMV is the main operation that must be distributed. It is memory-bound, so additional memory bandwidth can improve performance. That bandwidth can come from newer GPUs or from using multiple GPUs, which increases the aggregate available bandwidth.

Multi-GPU execution introduces two primary overheads:

• Communication overhead: Time spent waiting for data from other GPUs
• Load imbalance: Uneven work distribution that causes some GPUs to wait for others

NVIDIA's hardware and software stack provides mechanisms for addressing these costs. NVLink enables high-speed data transfers between GPUs, while NVSwitch coordinates communication among multiple devices. NCCL, a C/C++ library, provides GPU-direct point-to-point communication and collective operations such as AllReduce, AllGather, and Broadcast.

Together, NVLink, NVSwitch, and NCCL allow several GPUs to operate as a shared machine. Data can move between devices at high bandwidth while remaining within the GPU memory and communication system.

2D-partitioned D-PDLP

Earlier work on distributing PDLP includes D-PDLP. Its method reorders the matrix to improve load balance, partitions the matrix by rows and columns into sub-blocks, and distributes those blocks across GPUs. Each GPU performs an SpMV on its assigned block, and partial outputs for the same row are combined to reconstruct the result.

This method requires output-vector components to be assembled through communication among multiple GPUs. It demonstrates the potential of multi-GPU acceleration, but its two-dimensional partitioning treats the two SpMVs independently. cuOpt mPDLP instead uses the dependencies shared by consecutive SpMVs through the x and y vectors.

Min-cut partitioning in mPDLP

Min-cut partitioning considers both SpMVs in a PDLP iteration when distributing the constraint matrix. The goal is to keep closely connected input and output data on the same GPU whenever possible, reducing communication across consecutive SpMVs.

The sparsity pattern of A defines a bipartite graph. One set of nodes represents the primal variables x, and the other represents the dual variables y. An edge connects a primal variable to a dual variable when that variable is required to compute the corresponding matrix operation.

For example, when calculating y = Ax, the operations might be represented as:

y₁ = 1x₁ + 2x₂ + 4x₃ + 0x₄
y₂ = 3x₁ + 6x₂ + 0x₃ + 0x₄
y₃ = 0x₁ + 7x₂ + 1x₃ + 2x₄
y₄ = 0x₁ + 0x₂ + 3x₃ + 4x₄

Computing y₁ does not require x₄, and computing y₂ does not require x₃ or x₄. These dependencies can be represented as edges in the bipartite graph. If the graph is divided into k partitions for k GPUs, edges within one partition can be processed locally. Edges crossing partitions require communication between GPUs.

In a two-GPU example, GPU 1 might own y₁, y₂, x₁, and x₂, while GPU 2 owns the remaining rows and columns. When GPU 1 computes y₁, it already has local access to x₁ and x₂, but must retrieve x₃ from GPU 2. The connection between y₁ and x₃ is an edge cut, and it creates an additional communication requirement.

Min-cut partitioning attempts to reduce the number of these cross-partition edges. The method can be summarized as follows:

• A represents the LP's sparse constraint matrix
• A bipartite graph is constructed from the nonzero structure of A
• The graph is partitioned into k sections for k GPUs
• Each GPU receives rows and columns associated with its partition
• PDLP iterations run primarily on local data, with communication required for cut edges

Performance therefore depends on the sparsity pattern of A and the number of edge cuts created by the partitioning.

Benchmark results

The benchmark compared mPDLP with two other PDLP implementations: single-GPU cuOpt PDLP and the multi-GPU D-PDLP implementation.

More than 100 LP instances were tested from several datasets:

• Hans Mittleman LPFeas benchmark set
• Hans Mittleman LPFeas Addendum of large problems
• PDLP benchmark dataset
• D-PDLP benchmark dataset
• Open Energy Benchmark
• Problems supplied by NVIDIA industry partners

Each instance used a target tolerance of 10⁻⁶ and a one-hour time limit. The distributed solvers were tested on a node equipped with NVIDIA DGX B200 GPUs connected through an all-to-all NVLink topology. Single-GPU cuOpt PDLP used one B200 GPU.

Compared with single-GPU cuOpt, mPDLP showed a strong relationship between speedup and problem size. Speedups became noticeable when the number of nonzero entries exceeded 10⁷, and larger problems generally benefited more from the multi-GPU implementation.

The measured runtimes include iteration-independent overheads such as presolving, host-side transposition of A, graph partitioning, and postsolving. Only the PDLP iterations are distributed across GPUs. When measuring the PDLP portion alone, mPDLP reached an 11.4x speedup on tsp-gaia-10m, compared with a 4.2x end-to-end speedup.

Against D-PDLP on eight B200 GPUs, mPDLP again performed less strongly on smaller problems and became more competitive above 10⁷ nonzeros. For most large instances with more than 10⁷ nonzeros, mPDLP achieved speedups between 1.2x and 2.5x over D-PDLP.

On three ultra-large-scale instances, however, mPDLP was slower than D-PDLP. These instances came from established optimization benchmarks and may have sparsity patterns that differ from the real-world problems for which cuOpt mPDLP is tuned. The small number of ultra-large-scale cases also limits broader conclusions.

The results show that multi-GPU PDLP is most useful when the problem is large enough for computation to outweigh synchronization and data-transfer costs. On smaller models, communication at every iteration can exceed the benefit of parallel execution. Sparsity structure is also important: a high edge-cut ratio can create enough cross-GPU communication to reduce or eliminate the speedup.

Applications in the optimization ecosystem

NVIDIA has worked with optimization partners to integrate multi-GPU PDLP capabilities into existing tools and platforms.

• Kinaxis: A CPG supply-chain LP model with more than 135 million variables achieved a 3.3x speedup using cuOpt mPDLP on the Maestro platform. The test used eight NVLink-connected NVIDIA H100 GPUs.
• PSR: A stochastic energy expansion LP model with 185 million variables achieved a speedup of more than 5x using cuOpt mPDLP on eight NVLink-connected NVIDIA B200 GPUs.

Planned cuOpt mPDLP improvements

Further work on cuOpt mPDLP includes several areas:

• Min-cut partitioning: The current partitioning method does not account for computational load on each GPU. Weight-aware min-cut and hypergraph partitioning are intended to make distribution more load-aware and reduce imbalance and communication overhead.
• Overlapping communication and computation: The current SpMV loop performs communication and computation sequentially. Separating local and remote matrix components could allow local computation and remote data transfers to run in parallel, reducing idle GPU time.
• Feasibility polishing: Large LPs can require millions of PDLP iterations. Feasibility polishing prioritizes obtaining feasible solutions more quickly by trading some optimality for speed, potentially expanding the set of problems that can be solved within practical time limits.

Using mPDLP in NVIDIA cuOpt

The cuOpt mPDLP tutorial demonstrates how to run LP models across multiple GPUs. It covers cuOpt installation, solver configuration, and solve-time comparisons across different GPU counts. The examples can also be adapted to load an MPS file and measure performance on a specific workload.

The cuOpt source code and documentation provide additional information about APIs, solver settings, and deployment.