Scaling Decision Optimization to 100 Million Variables and Beyond with mPDLP in NVIDIA cuOpt

extract

Supply chain problems are expanding across more SKUs, lanes, and constraints than ever before, while energy grids are balancing more distributed sources in real time. To meet these challenges, teams need to evaluate larger models and more uncertainty within a practical planning window. NVIDIA cuOpt GPU-accelerated decision optimization can already deliver speedups of more than 10x over CPU solvers on a single GPU for large-scale linear programming (LP) problems. However, today’s largest planning problems can take hours to converge or exceed the memory capacity of a single GPU. These constraints cause a bottleneck for the workflows that drive business decisions. Historically, large and notoriously challenging LP problems have been a consistent subject of interest for researchers and practitioners. One such problem, zib03, introduced by Thorsten Koch in 2008 and featuring over 104 million nonzeros, has become a widely studied benchmark for new algorithms and solvers. Figure 1 shows how much progress has been achieved in the field with new algorithms and new hardware for solving this problem. cuOpt is now pushing the limits even further by solving the problem on multiple GPUs with mPDLP, achieving a solve rate nearly 10x faster than a year before. This post introduces the new cuOpt Multi-GPU Primal-Dual hybrid gradient for Linear Programming (mPDLP) solver. It distributes LP problems across NVLink-connected GPUs, enabling two major outcomes: - Reducing solve time for large problems that are impractical to solve within the time window available for a planner - Up to 6x lower peak memory usage per GPU compared to single-GPU PDLP (with LP problems capped at 2.1B nonzeros) The post also highlights the benefits of this approach through two partner applications: Kinaxis in globally constrained supply and production planning, and PSR in large-scale energy-system capacity-expansion models. How does the PDLP algorithm solve LPs? Solving an LP consists in finding \(x\) for the following equation: \begin{aligned} \min_{x} \quad & c^\top x \\ \text{s.t.} \quad & Ax = b \\ & l \leq x \leq u \end{aligned} Primal-Dual hybrid gradient for LP (PDLP) is a first order method for solving LPs. It is based on a gradient-descent-like method, making it highly parallelizable and GPU-friendly, enabling significant speedups. The hot loop of the algorithm is a Sparse Matrix-Vector Multiplication (SpMV), followed by several map operations. This method iterates on two vectors: - Primal solution (\(x\)) - Dual solution (\(y\)) The constraints of the LP are abstracted as a matrix (\(A\)). This is the matrix that is used in the SpMV loop. There is also the transpose of this matrix (\(A\)T). The iterative algorithm consists of these operations: \begin{aligned} x^{k+1} &= \text{proj}_{[l,\, u]}\!\left(x^k – \tau(c – A^\top y^k)\right) \\ y^{k+1} &= y^k + \sigma\left(b – A(2x^{k+1} – x^k)\right) \end{aligned} Here, τ and σ are the primal and dual step lengths, proj[l,u] denotes projection onto the variable bounds \(l ≤ x ≤ u\), \(b\) is the right-hand side of the linear equality constraints and \(c\) is the cost vector of the objective function we are trying to minimize. This algorithm is primarily composed of element-wise operations: projections and basic arithmetic. These operations are trivial to distribute on multiple GPUs, so they can be ignored while developing the distributed algorithm. Stripped down to its nontrivial communication steps, the algorithm is: \begin{aligned} x^{k+1} &= A^\top y^k \\ y^{k+1} &= Ax^{k+1} \end{aligned} The last value of \(y\) is used to compute the next \(x\) with a SpMV, and the last value of \(x\) is used to compute the next \(y\) with another SpMV. Distributing PDLP The only nontrivial task to distribute in PDLP is the SpMV, which is a memory-bound operation. This means that increasing the available bandwidth increases the speed of the operation. To increase the bandwidth, you can choose to either wait for the next generation of NVIDIA GPUs to be released or use more GPUs, effectively multiplying the total available bandwidth by the number of GPUs. This, however, comes with the costs associated with distributed algorithms. When distributing an algorithm, two main overheads must be minimized: - Communication overhead: Time spent waiting on information from other GPUs - Load imbalance: Uneven distribution of work across GPUs, causing some GPUs to wait for others NVIDIA provides the communication tools needed for efficient distributed computation for both the hardware and software stack. For the hardware stack, NVLink and NVSwitch provide the computation backbone. NVLink enables fast, optimized data exchanges between two GPUs, while the NVSwitch efficiently orchestrates the communication among several GPUs. For the software stack, NCCL is a C/C++ library that gives access to GPU-direct point-to-point communication and collective operations (AllReduce, AllGather, Broadcast) among GPUs. These operations natively use the NVLink and NVSwitch to their full capacity. With NVLink, NVSwitch, and NCCL together, multiple distinct GPUs become one fast shared machine. Data moves between devices at high bandwidth and low overhead, so distributed algorithms can scale without leaving the GPU. 2D partitioned D-PDLP Previous efforts to distribute PDLP across multiple GPUs include D-PDLP. First, the matrix can be reordered to improve load balancing. It is then partitioned along its rows and columns into several sub-blocks, which are distributed across the GPUs. Finally, the SpMV of the original matrix is computed by performing an SpMV on each sub-block and summing the resulting partial outputs corresponding to the same row to reconstruct the final result. This means each result of the output vector needs to be constructed from communication among multiple GPUs. This technique works well and D-PDLP used it to demonstrate the potential of multi-GPU acceleration. However, 2D partitioning in D-PDLP computes each SpMV independently, while cuOpt mPDLP exploits the fact that consecutive SpMVs share dependencies through the \(x\) and \(y\) vectors. Min-cut-partitioned PDLP Min-cut-partitioning leverages the dependencies shared between the two SpMVs performed in each PDLP iteration. Unlike the 2D partitioning approach, which treats the SpMVs independently, min-cut partitioning considers both operations jointly when determining how to distribute the constraint matrix across GPUs. The key idea is to partition the matrix so that tightly connected dependencies between the input and output vectors are kept on the same GPU whenever possible. Minimizing the communication across consecutive SpMVs by exploiting the sparsity pattern of \(A\). When performing \(y^{k+1} = Ax^{k+1}\) the underlying operations are: \begin{aligned} y_1 &= 1x_1 + 2x_2 + 4x_3 + \textcolor{gray}{0x_4} \\ y_2 &= 3x_1 + 6x_2 + \textcolor{gray}{0x_3} + \textcolor{gray}{0x_4} \\ y_3 &= \textcolor{gray}{0x_1} + 7x_2 + 1x_3 + 2x_4 \\ y_4 &= \textcolor{gray}{0x_1} + \textcolor{gray}{0x_2} + 3x_3 + 4x_4 \end{aligned} When computing \(y\), it’s unnecessary to use the full vector \(x\) for the different components \(y\)1, \(y\)2, \(y\)3, \(y\)4 in \(y\). This is due to the sparsity of \(A\). When performing the operations to compute \(y\)1: \(x\)4 is not necessary. For \(y\)2: \(x\)3 and \(x\)4 are not necessary. This is illustrated in the bipartite graph in Figure 3. Figure 3 shows a general overview of what data is needed to compute each row (this graph also represents the dependencies on SpMV for AT if read from right to left). This also means that if the graph gets partitioned in k-ways, k being the number of GPUs, the edges that are inside a partition can be computed locally, and the edges that are split among two partitions will require communication between the two GPUs associated with the edge. Back to the example, in Figure 4, the graph is partitioned with GPU 1 owning (that is, being responsible for the computation of) \(y\)1, \(y\)2, \(x\)1, \(x\)2 in the upper group and GPU 2 owning the other rows and columns in the lower group. This split allows for increased locality when performing the SpMVs. On GPU 1, when computing \(y\)1, data from \(x\)1, \(x\)2, and \(x\)3 is needed. However, \(x\)1 and \(x\)2 are local and were already computed in the previous step. It’s only necessary to pull \(x\)3 from GPU 2, which is remote. Communication is tightly linked with edge cuts. The \(y\)1-\(x\)3 edge is cut between the partitions, which results in one additional communication needed to compute \(y\)1. This means that partitioning the graph while reducing the number of edge cuts—depending on the sparsity pattern of \(A\)—might greatly reduce the need for communication. This is exactly what min-cut-partitioning does: a graph partitioning that reduces edge cuts. To summarize: - \(A\) is the sparse matrix representing the constraints of the linear program - First, materialize the bipartite graph defined by the matrix \(A\) - The graph is then partitioned k-ways (for k GPUs) to reduce edge cuts - Each GPU is assigned a subset of rows and columns according to the nodes in the partition - PDLP iterations are performed; work is mostly GPU-local, and communication occurs on edges that were cut in the partition This method is highly dependent on the number of edge cuts in the graph, resulting in performance variations depending on the sparsity pattern of matrix \(A\). Benchmarking results The benchmarking results presented here compare mPDLP with two other PDLP implementations: single-GPU cuOpt PDLP and multi-GPU-based D-PDLP. We ran the benchmarks on more than 100 LP instances spanning several datasets, including: - Hans Mittleman LPFeas benchmark set - Hans Mittleman LPFeas Addendum of large problems - PDLP benchmark dataset - D-PDLP benchmark dataset - Open Energy Benchmark - Problems from NVIDIA industry partners For each instance, we used the target tolerance 10−6, with a time limit of one hour. The distributed solvers, mPDLP and D-PDLP, are evaluated on a node equipped with NVIDIA DGX B200 GPUs connected through an all-to-all NVLink topology. The single-GPU cuOpt PDLP implementation is evaluated using a single B200 GPU. Figure 5 compares mPDLP cuOpt and single-GPU cuOpt. It shows that speedup strongly correlates with problem size. Speedups are noticeable when the number of nonzero entries (nnz) exceeds 107. Beyond this threshold, speedups increase with problem size, meaning larger problems benefit more from the multi-GPU implementation than from the single-GPU version. Note that these run times include several iteration-independent overheads: presolving, host-side transpose of A, graph partitioning and postsolving. Only the PDLP steps are distributed across workers, so only they benefit from the additional GPUs. Measuring the speedup of the PDLP steps alone, mPDLP reaches up to 11.4x on tsp-gaia-10m, well above the 4.2x measured end-to-end. Figure 6 compares mPDLP with D-PDLP. Again, lower performance is observed on small problems, while mPDLP becomes more competitive once the number of nonzero entries exceeds 107. On most of these large-scale instances (nnz>107), mPDLP achieves speedups between 1.2x and 2.5x over D-PDLP. However, on the three ultra-large-scale instances, mPDLP is slower than D-PDLP. It is worth noting that these instances come from established optimization benchmarks and may not be representative of the sparsity patterns found in the real-world problems that cuOpt mPDLP is tuned on. This could contribute to the observed slowdown, although the small number of ultra-large-scale instances also makes it difficult to draw a general conclusion. Multi-GPU PDLP delivers speedups only when problems are large enough, where communication overhead becomes negligible compared to compute time. On smaller problems, the constant cost of GPU synchronization and data transfer at every iteration outweighs the parallelism benefit. Sparsity structure also matters: problems with high edge-cut ratios after min-cut partitioning incur heavy cross-GPU communication that can neutralize or reverse the speedup. Integrating multi-GPU PDLP across the optimization ecosystem NVIDIA is working with partners across the optimization ecosystem to bring multi-GPU PDLP solver capabilities into the tools and platforms developers already use. - Kinaxis achieved a 3.3x speedup on a CPG supply chain LP model with over 135M variables using NVIDIA cuOpt mPDLP on their Maestro platform using eight NVLink-connected NVIDIA H100 GPUs. - PSR demonstrated more than 5x speedup on a stochastic energy expansion LP model with 185M variables using NVIDIA cuOpt mPDLP on eight NVLink-connected NVIDIA B200 GPUs. Continuing to improve cuOpt mPDLP NVIDIA is also working on the following cuOpt mPDLP improvements and features: - Min-cut partitioning: The current partitioning doesn’t account for computational load per GPU. Weight-aware min-cut and hypergraph partitioning will make the distribution load-aware, reducing imbalance and communication overhead across GPUs. - Hiding communications: The current SpMV hot loop communicates and computes sequentially; splitting the matrix into local and remote components will allow local computation and remote data transfer to run in parallel, reducing idle GPU time. - Feasibility polishing: Large-scale LPs often require millions of PDLP steps to converge. Feasibility polishing trades optimality for speed, delivering feasible solutions in fewer iterations and expanding the set of problems mPDLP can solve within practical time limits. Get started with mPDLP in NVIDIA cuOpt Run your LP models across multiple GPUs with the cuOpt mPDLP tutorial. Follow the notebook to install cuOpt, configure the solver, and compare solve times across GPU counts. Start with the provided examples, then load your own MPS file to measure performance on your workload. Explore the cuOpt source code on the NVIDIA/cuopt GitHub repo and refer to the product documentation for APIs, solver settings, and deployment guidance. Engage with the developer community on the cuOpt GitHub Discussion Forum. Stay up to date on NVIDIA cuOpt by subscribing to NVIDIA news and following NVIDIA AI on LinkedIn, X, Discord, and YouTube.