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.

  • 2
    $\begingroup$ Are you aware that BLAS implementations can use multiple cores when doing matrix operations? Actually for NVIDIA GPUs there is CuBLAS which shares the same interface as BLAS. Also see: or.stackexchange.com/questions/1024/… $\endgroup$ Commented May 22, 2021 at 17:54
  • 7
    $\begingroup$ Most LPs are very sparse, and the sparse linear algebra used in LP algorithms does not really work very well on GPUs. $\endgroup$ Commented May 23, 2021 at 1:45
  • 1
    $\begingroup$ @ErwinKalvelagen I know that newer A100 based card have hardware accelerated support for sparse matrices. $\endgroup$ Commented May 23, 2021 at 22:00
  • 2
    $\begingroup$ If your LP problems have a (near-)full constraint matrix, then simplex solvers based on dense numerical linear algebra have been developed, including GPU implementations. However, the requirement to solve dense LP problems has always been seen as almost non-existent, so the value of adding the necessary API for these codes - which are typically tested using random LPs - is not considered to be worthwhile. $\endgroup$ Commented May 24, 2021 at 8:49
  • 3
    $\begingroup$ SpMV alone wouldn't help much unless your LPs have very many more variables than constraints, since the SpMV operation in simplex generally takes only O(10%) of the solution time. So, once you've got your data to and from the GPU, any advantage may be lost. $\endgroup$ Commented May 24, 2021 at 8:52

2 Answers 2


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.

  • 2
    $\begingroup$ Hyper-sparsity isn't a property of the matrix. It's a property of solving linear systems and forming the matrix-vector product in the revised simplex method, and is observed when the solution is a sparse vector. $\endgroup$ Commented May 24, 2021 at 8:43
  • $\begingroup$ thanks for pointing this out :) $\endgroup$
    – mtanneau
    Commented May 24, 2021 at 12:23

As mentioned, both simplex method and barrier method don't suit GPUs. However, recently Google successfully gives a first-order method PDLP[1]. By using tricks like restarting, the performance of PDLP is improved but still not as good as other non-first-order-methods. But the difference is that this approach could be efficiently realized on GPU. With a little bit modification, a gpu version of PDLP in both Julia and C is given[2][3] and performance is compared. Sources are available at github[4][5].

[1] https://arxiv.org/abs/2106.04756

[2] https://arxiv.org/abs/2311.12180

[3] https://arxiv.org/abs/2312.14832

[4] https://github.com/jinwen-yang/cuPDLP.jl

[5] https://github.com/COPT-Public/cuPDLP-C


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.