8.2 Quadratic Programming

qp( P, q, [, G, h [, A, b[, solver]]])
Solves a convex quadratic program

\begin{displaymath}
\begin{array}{ll}
\mbox{minimize} & (1/2) x^TPx + q^T x \\
\mbox{subject to} & Gx \preceq h \\ & Ax = b.
\end{array}\end{displaymath}

P is a square dense or sparse real matrix, representing a symmetric matrix in 'L' storage, i.e., only the lower triangular part of P is referenced. G and A are dense or sparse real matrices. Their default values are sparse matrices with zero columns. q, h and b are single-column real dense matrices. The default values of h and b are matrices of size (0,1).

The default CVXOPT solver is used when the solver argument is absent or None. The MOSEK solver (if installed) can be selected by setting solver='mosek'.

qp() returns a dictionary with keys 'status', 'x', 's', 'y', 'z'. The possible values of the 'status' key are as follows.

'optimal'
In this case the 'x' entry is the primal optimal solution, the 's' entry is the corresponding slack in the inequality constraints, the 'z' and 'y' entries are the optimal values of the dual variables associated with the linear inequality and linear equality constraints. These values (approximately) satisfy the optimality conditions

\begin{displaymath}
Px + q + G^T z + A^T y = 0, \qquad Gx + s = h, \qquad
Ax = b, \qquad s \succeq 0, \qquad z \succeq 0, \qquad s^T z = 0.
\end{displaymath}

'primal infeasible'
This only applies when solver is 'mosek', and means that a certificate of primal infeasibility has been found. The 'x' and 's' entries are None, and the 'z' and 'y' entries are vectors that approximately satisfy

\begin{displaymath}
G^T z + A^T y = 0, \qquad h^T z + b^T y = -1, \qquad z \succeq 0.
\end{displaymath}

'dual infeasible'
This only applies when solver is 'mosek', and means that a certificate of dual infeasibility has been found. The 'z' and 'y' entries are None, and the 'x' and 's' entries are vectors that approximately satisfy

\begin{displaymath}
Px = 0, \qquad q^Tx = -1, \qquad Gx + s = 0, \qquad Ax=0, \qquad
s \succeq 0.
\end{displaymath}

'unknown'
This means that the algorithm reached the maximum number of iterations before a solution was found. The 'x', 's', 'y', 'z' entries are None.

As an example we compute the trade-off curve on page 187 of the book Convex Optimization, by solving the quadratic program

\begin{displaymath}
\begin{array}{ll}
\mbox{minimize} & -\bar p^T x + \mu x^T S ...
...ox{subject to} & {\bf 1}^T x = 1, \quad x \succeq 0
\end{array}\end{displaymath}

for a sequence of positive values of mu. The code below computes the trade-off curve and produces two figures using the Matplotlib package.
\includegraphics[width=10cm]{figures/portfolio1.eps} \includegraphics[width=10cm]{figures/portfolio2.eps}

from math import sqrt
from cvxopt.base import matrix
from cvxopt.blas import dot 
from cvxopt.solvers import qp
import pylab

# Problem data.
n = 4
S = matrix([[ 4e-2,  6e-3, -4e-3,    0.0 ], 
            [ 6e-3,  1e-2,  0.0,     0.0 ],
            [-4e-3,  0.0,   2.5e-3,  0.0 ],
            [ 0.0,   0.0,   0.0,     0.0 ]])
pbar = matrix([.12, .10, .07, .03])
G = matrix(0.0, (n,n))
G[::n+1] = -1.0
h = matrix(0.0, (n,1))
A = matrix(1.0, (1,n))
b = matrix(1.0)

# Compute trade-off.
N = 100
mus = [ 10**(5.0*t/N-1.0) for t in xrange(N) ]
portfolios = [ qp(mu*S, -pbar, G, h, A, b)['x'] for mu in mus ]
returns = [ dot(pbar,x) for x in portfolios ]
risks = [ sqrt(dot(x, S*x)) for x in portfolios ]

# Plot trade-off curve and optimal allocations.
pylab.figure(1, facecolor='w')
pylab.plot(risks, returns)
pylab.xlabel('standard deviation')
pylab.ylabel('expected return')
pylab.axis([0, 0.2, 0, 0.15])
pylab.title('Risk-return trade-off curve (fig 4.12)')
pylab.yticks([0.00, 0.05, 0.10, 0.15])

pylab.figure(2, facecolor='w')
c1 = [ x[0] for x in portfolios ] 
c2 = [ x[0] + x[1] for x in portfolios ]
c3 = [ x[0] + x[1] + x[2] for x in portfolios ] 
c4 = [ x[0] + x[1] + x[2] + x[3] for x in portfolios ]
pylab.fill(risks + [.20], c1 + [0.0], '#F0F0F0') 
pylab.fill(risks[-1::-1] + risks, c2[-1::-1] + c1, '#D0D0D0') 
pylab.fill(risks[-1::-1] + risks, c3[-1::-1] + c2, '#F0F0F0') 
pylab.fill(risks[-1::-1] + risks, c4[-1::-1] + c3, '#D0D0D0') 
pylab.axis([0.0, 0.2, 0.0, 1.0])
pylab.xlabel('standard deviation')
pylab.ylabel('allocation')
pylab.text(.15,.5,'x1')
pylab.text(.10,.7,'x2')
pylab.text(.05,.7,'x3')
pylab.text(.01,.7,'x4')
pylab.title('Optimal allocations (fig 4.12)')
pylab.show()