# 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:

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.

## 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:

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.