Quadratic Programming in Python

Quadratic programs are a particular kind of mathematical optimization problems
that can be applied to solve a variety of problems, for instance: in statistics
for curve fitting, in machine learning to compute support vector machines
(SVMs), in robotics to
solve inverse kinematics, etc. They are
the first step beyond linear programming (LP) in convex optimization. We will
now see how to solve quadratic programs in Python using a number of available
solvers: CVXOPT, CVXPY, Gurobi, MOSEK, qpOASES and quadprog.

Here, \(x\) is the vector of optimization variables \(x_1, \ldots,
x_n\). The matrix \(P\) and vector \(q\) are used to define any
quadratic objective function on these variables, while the matrix-vector
couples \((G, h)\) and \((A, b)\) are used to define inequality and
equality constraints, respectively. Vector inequalities apply coordinate by
coordinate.

This mathematical construction means the following: a QP finds the minimum of a
quadratic function over a linear set:

In the 2D illustration above, the level sets of the quadratic function are
represented by ellipses, and the linear set is the blue polygon. Since the
global optimal of the objective function is outside of the polygon, the
solution \(x^*\) of the QP lies at the boundary of the linear set. The set
of linear constraints that are saturated at \(x^*\) is called the active
set, but that's a story for
another post...

Back to the standard form, notice that there is no constant term in the
objective function. Indeed, it would have no effect on the result of the
optimization, which is the location of the solution \(x^*\). Expressions
like \(\| A x - b \|^2\) are written in standard form as \(x^T (A^T A)
x - 2 (A^T b)^T x\).

The standard form also assumes without loss of generality that the matrix
\(P\) is symmetric. Any matrix \(M\) can be decomposed as sum of its
symmetric part \(M^+\) and antisymmetric part \(M^-\), and the latter
yields zero in \(x^T M^- x\). Note that some QP solvers assume that you
feed them a symmetric cost matrix: they won't check this, and will return wrong
results if you don't.

Setting up QP solvers

The two readily-available QP solvers in Python are CVXOPT and quadprog. They
can be installed by:

CVXOPT uses its own matrix type, and it requires the matrix \(P\) of the
objective function to be symmetric. To be on the safe side, you can wrap it as
follows:

defcvxopt_solve_qp(P,q,G=None,h=None,A=None,b=None):P=.5*(P+P.T)# make sure P is symmetricargs=[matrix(P),matrix(q)]ifGisnotNone:args.extend([matrix(G),matrix(h)])ifAisnotNone:args.extend([matrix(A),matrix(b)])sol=cvxopt.solvers.qp(*args)if'optimal'notinsol['status']:returnNonereturnnumpy.array(sol['x']).reshape((P.shape[1],))

The quadprog module works directly on NumPy arrays so there is no need for
type conversion. Its matrix representation is equivalent but with different
names:

defquadprog_solve_qp(P,q,G=None,h=None,A=None,b=None):qp_G=.5*(P+P.T)# make sure P is symmetricqp_a=-qifAisnotNone:qp_C=-numpy.vstack([A,G]).Tqp_b=-numpy.hstack([b,h])meq=A.shape[0]else:# no equality constraintqp_C=-G.Tqp_b=-hmeq=0returnquadprog.solve_qp(qp_G,qp_a,qp_C,qp_b,meq)[0]

In case you want to try out other solvers, I implemented similar wrappers for
those I could get my hands on (like Gurobi or MOSEK) in the qpsolvers module.

Comparing solver performances

The three others are symbolic, meaning that if you dig into their API they
allow you to construct your problem formally (with variable names) rather than
using the matrix-vector representation. This is convenient for big sparse
problems, but slower and small problems such as the one we are looking at here.
The three symbolic solvers I tested are:

The Toeplitz matrix used to generate inequalities is just an upper-tridiagonal
matrix with coefficients 1, 2, 3, all other coefficients being zero. This
matrix is sparse but represented by (dense) NumPy arrays here. Using the
function above, I generated a benchmark for problem sizes ranging from 10 to
2,000, averaging computation times over 10 runs for each point. Here are the
results:

The bottom line of this small comparison is that quadprog, which implements
the Goldfarb-Idnani dual algorithm, simply rocks. More generally,
active-set solvers (quadprog and qpOASES) perform best on these dense problems.
To see the benefit of symbolic solvers (CVXPY or Gurobi), one would have to use
sparse matrix representation, which I didn't do here.

You can try for yourself on your own machine, all scripts and solver wrapping
functions are in the qpsolvers
repository.