Is there a constrained nonlinear optimization library like IPOPT that runs on GPUs?

Question

Somebody on my team wants to paralelize IPOPT. (at least some of functions of it). I have not been able to find a GPU implementation of it or a similar package. Nor I've found anything on their docs.

So the question is, is there an alternative already implemented on GPU? or at least somebody working on porting it to GPU so we can work toguether?

Answer 1

tl;dr: My general impression from the literature is that speedups are modest (if they exist). The main kernel you'll see in these methods is a sparse-direct method (e.g., sparse LU, sparse LDLT), and memory accesses are irregular; these characteristics don't favor use of GPUs. Also, parallel IPMs are in their infancy. I still suspect people will work on GPU implementations, but I am skeptical much will come of them. (However, distributed-memory IPMs seem slightly more promising and useful.)

---

A few people have worked on interior-point methods (IPMs) for GPUs:

• Smith, Gondzio, and Hall developed an IPM for linear programs (LPs) that achieved 4-10x speedups
• Jung and O'Leary look at some of the linear algebra kernels in IPMs for LPs and see modest speedups for GPUs on larger problems in their test set

In general, the papers don't compare their work to highly tuned LP solvers. Jung and O'Leary compare to linprog, which wouldn't be the choice of most practitioners, and the impression I get from scanning the Smith, et al paper is that a serial interior-point method was written from scratch. Given that Hall is very involved in open-source LP solvers, I take that work slightly more seriously. One of the better open-source LP solvers out there is Clp, which Hall maintains, and if they had used that code, it would be mentioned by name. So I'd pay more attention if these codes were accelerating already highly tuned solvers, instead of one-off serial comparisons that aren't state-of-the-art.

That said, the existing work is for LPs, and you're probably asking for a nonlinear programming (NLP) solver, because that's what IPOPT is. Here is what I know about the NLP solver case:

• parallelizing IPMs is generally challenging (for the LP & SOCP case, see Elemental; you'd use something like the SOCP version (or a QP version) for a general IPM
• if your NLP doesn't have obvious structure (e.g., isn't a stochastic program, doesn't have obvious bordered-block-angular structure, isn't a PDE-constrained optimization problem where preconditioners for the PDE are well-known), the best linear algebra methods involve sparse-direct methods (e.g., LU, LDLT) with irregular memory accesses that aren't typically suited to GPU architectures (as mentioned before)
• if your NLP does have obvious structure, you'd want to use an inertia-free line search algorithm that doesn't require an inertia-revealing factorization (building upon recent work by Chiang and Zavala), so you could use known preconditioners to help solve the Karush-Kuhn-Tucker system
• the best parallel matrix interior-point methods available that I know of are those by Gondzio(Implementation paper, theory paper), and by Waechter, et al. (Large-scale interior point overview paper, paper going into more detail about implementation of large-scale interior point methods); there are also algorithms by Petra, Zavala, et al that exploit structure (structured nonconvex optimization paper, Schur complement for interior-point methods paper).

Most research seems to be working on parallelism for the distributed-memory case (which is hard enough). There is some work on the shared-memory case for LP & QP solvers, because that work makes it into commercial codes (e.g., Gurobi, CPLEX, Xpress). Existing work is interesting, but there doesn't seem to be anything compelling for GPUs yet, excepting special applications (e.g., machine learning, for which you can use different algorithms better suited to GPUs).

Answer 2

I'm a little late to the party, but the short answer is that yes it's possible to parallelize an interior point method for GPUs, but whether or not that is successful depends on the structure of the problem. In terms of existing software, Optizelle can do it. Grab the development branch until a new release occurs in the near future.

The situations differ slightly depending on whether or not the original problem contains equalities or inequalities. There's a variety of ways to do this, but, in my opinion, the best way to do this for problems with only inequalities constraints is using an inexact trust-region method Newton method combined with an primal dual interior point method.

For inequalities only, the basic inexact trust-region Newton method can be found in Nocedal and Wright's Numerical Optimization on page 171 or on Conn, Gould, and Toint's Trust-Region Methods on page 205. This algorithm can be combined successfully with a primal-dual interior point method by essentially using the modified truncated-CG method from page 890 of the paper An Interior Point Method for Large-Scale Nonlinear Programming by Byrd, Hribar, and Nocedal. Personally, I don't like how they setup their interior point system, so I wouldn't use their interior point formulation, but that's preference. NITRO is a good algorithm. As far as the interior point details, Optizelle's manual explains how to do this in its manual. I probably ought to post an updated manual, but the development branch is current.

For the case with both inequality and equality constraints, I believe the best algorithm is combining the inexact trust-region composite-step SQP method from Heinkenschoss and Ridzal in a paper titled A Matrix-Free Trust-Region SQP Method for Equality Constrained Optimization. Basically, the process of tacking on an interior point method works pretty much the same as the unconstrained case except that the quasinormal step needs to be safeguarded as well.

As far as the parallelization opportunities, the algorithms I reference above work well because these algorithms can be implemented matrix-free. Specifically, Optizelle's implementation for the problem

$\min\limits_{x\in X}\{ f(x) : g(x)=0, h(x)\geq 0\}$

Requires that the user provide an implementation for

• $f(x), \nabla f(x), \nabla^2 f(x)\partial x$

• $g(x), g^\prime(x)\partial x, g^\prime(x)^*\partial y, (g^{\prime\prime}(x)\partial x)^* \partial y$

• $h(x), h^\prime(x)\partial x, h^\prime(x)^*\partial y, (h^{\prime\prime}(x)\partial x)^* \partial y$

It doesn't care where these implementations come from or how their parallelized. They can be done in shared memory, distributed memory, or GPUs. It does not matter. What works best for a particular problem, depends on the structure. In addition, it requires the user to provide linear algebra for

• init: Memory allocation and size setting
• copy: y <- x (Shallow. No memory allocation.)
• scal: x <- alpha * x
• axpy: y <- alpha * x + y
• innr: innr <- <x,y>
• zero: x <- 0
• rand: x <- random
• prod: Jordan product, z <- x o y
• id: Identity element, x <- e such that x o e = x
• linv: Jordan product inverse, z <- inv(L(x)) y where L(x) y = x o y
• barr: Barrier function, barr <- barr(x) where x o grad barr(x) = e
• srch: Line search, srch <- argmax {alpha \in Real >= 0 : alpha x + y >= 0} where y > 0
• symm: Symmetrization, x <- symm(x) such that L(symm(x)) is a symmetric operator

These operations can be done in serial, parallel, distributed memory, shared memory, or on GPUs. It does not matter. What's best depends on the problem structure.

Finally, there's the linear systems and there are three that can be provided:

• Preconditioner for $\nabla^2 f(x)$
• Left preconditioner for $g^\prime(x)g^\prime(x)^*$
• Right preconditioner for $g^\prime(x)g^\prime(x)^*$

Each of these preconditioners can be implemented in either serial or parallel, distributed memory or shared, or on GPUs. Basically, the first preconditioner is the preconditioner for the truncated-CG system associated with the optimalitity systems. The last two preconditioners are used for the augmented system solves associated with the composite step SQP algorithm. In general, this is where you're going to get your biggest performance boost. Imagine if the constraint $g$ represented some kind of PDE. Then, the preconditioner for $g^\prime(x)g^\prime(x)^*$ corresponds to a forward PDE solve followed by an adjoint PDE solve. Note, if they were square, $(g^\prime(x)g^\prime(x)^*)^{-1}=g^\prime(x)^{-*}g^\prime(x)^{-1}$. For a huge number of PDE formulations, such as finite difference methods with explicit time integrators, these solves can be very well parallelized on a GPU.

Finally, the algorithms in Optizelle do work on symmetric cone problems, which include bound, second-order cone, and semidefinite constraints. Nevertheless, in general, the linear cone solves will tend to out perform it. Basically, linear cone solves can reduce the feasibility and optimality solves done to a really compact system that can be Choleski factored. Since Optizelle works with nonlinear systems, it can't really do that. At least, I don't know how. On top of that, there are restrictions on the size of the SDP blocks that Optizelle can handle. The operator linv above requires the inverse of the SDP matrices and that inverse is really expensive for large blocks. In addition, there's an extra safe guard that requires a Choleski factorization. These factorizations don't really parallelize well on a GPU. At least, I don't know of an implementation that does parallelize well. Anyway, the bottom line is that if it's a linear cone program, use a linear cone solver like CSDP or SDPT3.

TLDR; Use Optizelle. It's free, open source, and BSD licensed. I've scaled it to something like half a billion variables and it worked fine. I've run it with GPUs and it worked fine. Whether or not it works well with a GPU depends on whether or not the operations above parallelize well on a GPU.

Answer 3

Generally nonlinear optimization is hard to parallelize because its stepping algorithm is very sequential. GPUs are slower than CPUs so you only get a speedup if you have your problem as something that is highly parallel (i.e. thousands of threads). Thus I wouldn't expect much of a speedup (or any, it might usually be slower) by putting the nonlinear solver on the GPU. That's most likely why no one has done it well, and most people wouldn't try it.

On the other hand, if your objective and constraint functions can be calculated in a highly parallel (or vectorized) fashion, you can get a massive speedup by solving your objective/constraint functions on the GPU. This is a common way to use GPUs with nonlinear optimization since it uses GPU acceleration on the hardest (and most parallel) part of the code, and thus is probably the most efficient.