Sparse Quadratic Programming with Ipoptr

This post is a follow up to my last post on quadratic programming facilities in R. A commenter pointed me to the ipoptr project which exposes an R interface to the COIN-OR optimization routine Ipopt. COIN-OR is a suite of optimization utilities implemented in C++ and supported by a back-end of configurable FORTRAN linear system solvers. ipoptr may be a good solution for users wishing to solve sparse nonlinear constrained optimization problems through an R frontend. Some highlights of this solver are:

  1. It is a general interior point solver that can handle nonlinear objectives with nonlinear constraints. In particular, no convexity assumptions are required to obtain local solutions.
  2. It offers a flexible and lightweight specification of the objective function and a sparse matrix representation for the constraints and other problem data.

In this post, I’ll explain how ipoptr can be applied to solve quadratic programs and I’ll compare the performance of this solver to other quadratic program solvers (quadprog, ipop) available in R. We’ll see that ipoptr is very fast and efficient on large sparse quadratic programs, seemingly an order of magnitude faster than quadprog on the demonstration problem considered in my previous post. Because the Ipopt backend is a bit tricky to install, the last section provides a detailed overview of how I successfully built this solver under Ubuntu Linux.

The ipoptr interface

Before applying ipoptr to quadratic programming, it may be helpful to present a high-level picture of the ipoptr interface. The following is a summary of the ipoptr interface drawn from the documentation. ipoptr stages and invokes an interior point algorithm from within R to find a local solution of the following constrained optimization problem:

    \[\begin{aligned}\underset{\alpha \in \mathbb{R}^n}{\text{Minimize}}: \qquad & f(x) \\ \text{Subject to:} \qquad & g_L \leq g(x) \leq g_U \\ & x_L \leq x \leq x_U \end{aligned}\]

where:

  • f:\mathbb{R}^n \to \mathbb{R} the objective function is twice continuously differentiable
  • g:\mathbb{R}^n \to \mathbb{R}^m defines the constraint set and is twice continuously differentiable
  • x_U, and x_L are fixed vectors in \mathbb{R}^n
  • g_U, and g_L are fixed vectors in \mathbb{R}^m

Note that the solver does not need to make any assumptions about the convexity of either f or g. The solver makes use of the following ingredients to find a solution:

  1. An initial guess x_0 for the solution
  2. The objective function gradient \nabla f
  3. the Jacobian of the constraint map J_g
  4. Optionally, the Hessian of the problem Lagrangian in the form:

        \[H(x, \sigma_f, \vec{\lambda}) = \sigma_f \nabla^2 f(x) + \sum_{i=1}^m \lambda_i \nabla^2 g_i(x).\]

When n is large, the dense matrix representations of g, \nabla f, and J_g will have a substantial footprint. In application problems, however, the Jacobian and/or the Hessian object will often be sparse.  For efficiency, ipoptr utilizes a sparse matrix representation for H and J_g.  This format is defined as follows. Let A_{\mathrm{values}} be the list of all non-zero values in A (read into the list from left to right along the rows of A).  Then let A_{\mathrm{mask}} be a list of lists, where the list at position i in A_{\mathrm{mask}} provides the column indexes of the nonzero elements in the i-th row of A. Then (A_{\mathrm{values}}, A_{\mathrm{mask}}) is a sparse representation of A.

Generally, it could be quite difficult to come up with the explicit representations of the Jacobian and Hessian that ipoptr consumes. Although the Hessian is an optional argument, including it can significantly improve convergence behaviour (try taking the Hessian arguments out in the examples below). For a quadratic program, however, it is simple to compute the Jacobian and Hessian directly. We’ll do this in the next section and then provide an R wrapper function that can transform any symmetric positive definite quadratic program in standard matrix form into a quadratic program that can be solved by ipoptr.

Solving Quadratic Programs with ipoptr

Consider the quadratic program:

    \[\begin{aligned} \underset{\alpha \in \mathbb{R}^n}{\text{Minimize}}: \qquad & f(x) = -d^Tx + \frac{1}{2}x^T D x \\ \text{Subject to:} \qquad & g(x) = Ax \leq b \end{aligned}\]

where

  • D is an n \times n matrix (which we’ll assume is symmetric)
  • A is an m \times n matrix.
  • d is a fixed vector in \mathbb{R}^n and b is a fixed vector in \mathbb{R}^m.

Then using some basic identities of vector calculus and the symmetry of D:

    \[\begin{aligned} \nabla f(x) &= -d + \frac{1}{2}(D^T + D) x = -d + Dx \\ \nabla^2 f(x) &= D \\ J_g(x) &= A^T \\ \nabla^2 g_i(x) &= 0 \mbox{ for all } i \\ H(x, \sigma_f, \vec{\lambda}) &= \sigma_f D. \end{aligned}\]

Because D is positive definite, the global minimizer \tilde{x} of f can be computed as:

    \[\nabla f(\tilde{x}) = 0 \rightarrow \tilde{x} = D^{-1} d.\]

Goldfarb and Idnani, whose QP algorithm is implemented in the quadprog package, use the global minimizer \tilde{x} as the initial point for their interior point primal/dual method. This choice is based on empirical evidence cited by the authors that this selection can substantially reduce the number of iterations required to find a solution. We’ll follow suit with ipoptr by setting x_0 =\tilde{x} for our initial guess.

The following R function uses these facts to solve a quadratic program in standard form with ipoptr:

Performance Comparison

In my previous post, I compared the performance of quadprog and kernlab’s ipop solvers on the circus tent demonstration problem. For this post, I repeated the same experiment with the ipoptr solver and the results were very impressive. Ipoptr was substantially faster at solving this sparse demonstration problem than quadprog. Here are the solve times I obtained:

Here is the experiment code:

It’s not so surprising that the ipoptr package is substantially faster in this example because the matrix Dmat that defines the objective function for the tent problem is very sparse. Quadprog requires dense matrix representations to describe the QP that we wish to solve, while in ipoptr we need only specify the functional form the quadratic objective function. Here is the plot of Dmat’s sparsity structure:

To investigate the performance of ipoptr on dense QPs, I generated some random positive definite QPs of differing sizes and compared solve times against the solver from quadprog and the ipop solver from the kernlab package. In this dense problem case, quadprog appears to have a significant advantage over ipoptr. Still, ipoptr shows solid performance and is considerably faster than the pure R implementation of the interior point solver found in kernlab’s ipop.

Here is the code for the randomly generated QP experiment:

Conclusions and Comments

In this post, we looked at how the ipoptr interface can be used to solve quadratic programming problems and compared the Ipopt solver to other quadratic program solvers available in R. We found:

  • ipoptr is very efficient and well-suited for solving large, sparse quadratic programs, even outperforming the quadprog package solver in this case.
  • However, for dense symmetric positive definite QPs, quadprog is much faster than ipoptr.
  • Both solvers are substantially faster than the pure R interior point solver ipop from the kernlab package.

Some additional remarks:

  • For the experiments in this post, Ipopt was configured to use the MUMPS linear solver. It is possible to configure Ipopt to use other linear solvers and this could have an impact on performance.
  • Initializing Ipopt with x_0 as the global minimum of the objective function (as done internally in the quadprog package) seems to be somewhat beneficial for reducing runtimes. This effect should be investigated more carefully, however.
  • The timings presented above do not include the time required to translate from a QP specified in matrix format (e.g. Dmat, Amat, etc.) to the input format required by ipoptr. In practice, however, this step is rather expensive and users will probably want to generate the ipoptr input data directly instead of forming dense matrix inputs and then converting to the ipoptr format. The specifics of how this should be done will depend on the particular problem which the user is trying to solve.

 

Installing ipoptr for Ubuntu Linux

Because of complex licensing issues in COIN-OR suite and its many dependencies, it’s not possible to bundle the backend ipopt solver and its R interface into a simple package installable directly from R. Thus, installing ipoptr won’t be as easy as installing most other R pacakges. Fortunately, the CoinOpt documentation and the R interface documentation are fairly complete.

After a bit of work, I was able to install everything on a Ubuntu-Linux box using the following steps.

Now, we can fire up R and install the package from source using:

4 thoughts on “Sparse Quadratic Programming with Ipoptr

  1. Great Post and helpful explanation! Quick comment: It seems that your r-code was deleted somehow after the line "The following R function uses these facts to solve a quadratic program in standard form with ipoptr: ". Would you mind re-posting?

  2. Thanks for noticing that Chandler! I recently changed my github user name which broke all the embedded gists in this post. I've fixed all the links and the blog should render correctly now.

Leave a Reply

Your email address will not be published. Required fields are marked *