Linear Optimization Library for C++ with GPU Support

Question

Does anyone know any linear optimization libraries for C++ supporting GPUs for parallelization? If multiple, which do you recommend?

The GPU support is important to me since I am dealing with large matrices in my problem constraints where the solving section is pretty much time-consuming without parallelization.

I myself found 1-2 PhD theses but they hadn't published their implementation as a library publicly.

Answer 1

For context: most (if not all) major LP solvers are built on 2 algorithms: the simplex method, and the interior-point method.

The simplex method is intrinsically sequential: you're doing a lot of (cheap) operations called pivots, and the matrices involved are usually sparse. At each pivot, you essentially perform a rank-one update of a sparse LU factorization, plus a few matrix-vector products. Unless your problem exhibits some particular structure (as in this paper), the potential for parallelism is relatively low (even on CPU), especially when dealing with problems that exhibit hyper sparsity

In the interior-point method, the main cost at each iteration is computing the Cholesky factorization of a matrix of the form $A D A^{T}$, where $A$ is the constraint matrix of your problem (in standard form) and $D$ is diagonal. Again, the matrices involved are sparse, but Cholesky factorization is more amenable to parallelism. Here, modern LP solvers will use multiple (CPU) threads when available.

From there, here are the 3 main reasons that LP solvers typically don't run on a GPU:

1. Because the matrices are sparse, the potential GPU speed-up is less dramatic than, e.g., in machine learning (where everything mostly boils down to dense matrix-matrix multiplication). Even just handling memory requirements on a GPU can quickly become a nightmare.
2. The available GPU libraries for sparse linear algebra are not yet as numerically robust and broad as existing CPU ones. For instance, as far as I know, there is no sparse Cholesky factorization on CuSparse.
3. Solver development teams are usually small, and building/maintaining GPU support is a lot of work, which may simply not be worth it at the moment.

One further reason: LP solvers perform computations in double-precision floating-point arithmetic. However, recent trends in GPU development are moving towards lower-precision arithmetic, such as single-precision or even lower (see NVIDA's recent "Tensor32" arithmetic). When looking at GPUs' flop-count for double-precision, the potential speedups are less dramatic.