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

PUBLISHER ARTICLE PREVIEW

From the original article

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:

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.

begin{aligned}min_{x} quad & c^top x \text{s.t.} quad & Ax = b \& l leq x leq uend{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.

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).

begin{aligned}x^{k+1} &= text{proj}_{[l,, u]}!left(x^k – tau(c – A^top y^k)right) \y^{k+1} &= y^k + sigmaleft(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.

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.