Quadratic Programming

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.

A quadratic program (QP) is written in standard form as:

\begin{equation*} \begin{array}{rl} \mathrm{minimize} & (1/2) x^T P x + q^T x \\ \mathrm{subject\ to} & G x \leq h \\ & A x = b \end{array} \end{equation*}

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:

QP is about finding the minimum of a quadratic function (circles) over a linear set (polygon)

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.

Solving QPs in Python

The most readily-available QP solver in Python is quadprog. It can be installed by simply running sudo pip install quadprog if you have the pip package manager on your system. To keep the same matrix names as above, let us wrap it with the following function:

import quadprog

def solve_qp(P, q, G=None, h=None, A=None, b=None):
    qp_G = .5 * (P + P.T)   # make sure P is symmetric
    qp_a = -q
    if A is not None:
        qp_C = -numpy.vstack([A, G]).T
        qp_b = -numpy.hstack([b, h])
        meq = A.shape[0]
    else:  # no equality constraint
        qp_C = -G.T
        qp_b = -h
        meq = 0
    return quadprog.solve_qp(qp_G, qp_a, qp_C, qp_b, meq)[0]

For a small example, let us see how to solve the following QP:

\begin{align*} \begin{array}{rl} \mathrm{minimize} & \left\| \left[\begin{array}{ccc} 1 & 2 & 0 \\ -8 & 3 & 2 \\ 0 & 1 & 1 \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ x_3\end{array}\right] - \left[\begin{array}{c} 3 \\ 2 \\ 3\end{array}\right] \right\|^2 \\ \mathrm{subject\ to} & \left[\begin{array}{ccc} 1 & 2 & 1 \\ 2 & 0 & 1 \\ -1 & 2 & -1 \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ x_3\end{array}\right] \leq \left[\begin{array}{c} 3 \\ 2 \\ -2 \end{array} \right] \end{array} \end{align*}

First, we write our QP matrices in proper format:

M = array([[1., 2., 0.], [-8., 3., 2.], [0., 1., 1.]])
P = dot(M.T, M)
q = dot(array([3., 2., 3.]), M).reshape((3,))
G = array([[1., 2., 1.], [2., 0., 1.], [-1., 2., -1.]])
h = array([3., 2., -2.]).reshape((3,))

Finally, we compute the solution using the function above:

In [1]: solve_qp(P, q, G, h)
Out[1]: array([-0.49025721, -1.57755261, -0.66484801])

To go further, you can check out this blog post on Quadratic Programming in Python, as well as the qpsolvers which implements the solve_qp function for a variety of QP solvers.

Content on this website is under the CC-BY 4.0 license.