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

來源原標題

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

作者

Tanya Lenz

原始標籤

Agentic AI / Generative AI · Developer Tools & Techniques · Simulation / Modeling / Design · Top Stories · Algorithms / Numerical Techniques · linear programming

原始連結

https://developer.nvidia.com/blog/scaling-decision-optimization-to-100-million-variables-and-beyond-with-mpdlp-in-nvidia-cuopt/

摘要原文

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. Figure 1. Progression of solve times by linear programming solvers on problem zib03 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\). Figure 2. Example Sparse-Matrix-Vector multiplication between the constraint matrix \(A\) and the primal solution \(x\) , yielding \(y\) 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. Bipartite dependency graph linking primal variables \(x\) to dual variables \(y\) based on matrix sparsity Figure 3 shows a general overview of what data is needed to compute each row (this graph also represents the dependencies on SpMV for A T if read from right to left).