8.5 Nonlinear Convex Programming

cp( F[, G, h[, A, b]])

Solves an optimization problem

$\begin{displaymath} \begin{array}{ll} \mbox{minimize} & f_0(x) \\ \mbox{subj... ... k=1,\ldots,m \\ & G x \preceq h \\ & A x = b, \end{array}\end{displaymath}$

(8.3)

with $f=(f_0,\ldots,f_m)$ convex and twice differentiable.

F is a function that evaluates the objective and nonlinear constraint functions. It must handle the following calling sequences.

F() returns a tuple (m, x0), where m is the number of nonlinear constraints and x0 is a point in the domain of f. x0 is a dense real matrix of size (n,1).
F(x), with x a dense real matrix of size (n,1), returns a tuple (f, Df). f is a dense real matrix of size (m+1,1), with f[k] equal to f_k(x). (If m is zero, f can also be returned as a number.) Df is a dense or sparse real matrix of size (m+1,n) with Df[k,:] equal to the transpose of the gradient of f_k at x. If x is not in the domain of f, F(x) returns None or a tuple (None,None).
F(x,z), with x a dense real matrix of size (n,1) and z a positive dense real matrix of size (m+1,1) returns a tuple (f, Df, H). f and Df are defined as above. H is a square dense or sparse real matrix of size (n, n), whose lower triangular part contains the lower triangular part of

$\begin{displaymath} z_0 \nabla^2f_0(x) + z_1 \nabla^2f_1(x) + \cdots + z_m \nabla^2f_m(x). \end{displaymath}$

If F is called with two arguments, it can be assumed that x is in the domain of f.

G and A are dense or sparse real matrices with n columns. Their default values are matrices of size (0,n). h and b are dense real matrices with one column, and the same number of rows as G and A, respectively. Their default values are matrices of size (0,1).

cp() returns a dictionary with keys 'status', 'x', 'snl', 'sl', 'y', 'znl', 'zl'. The possible values of the 'status' key are:

'optimal'

In this case the 'x' entry of the dictionary is the primal optimal solution, the 'snl' and 'sl' entries are the corresponding slacks in the nonlinear and linear inequality constraints, and the 'znl', 'zl' and 'y' entries are the optimal values of the dual variables associated with the nonlinear inequalities, the linear inequalities, and the linear equality constraints. These vectors approximately satisfy the Karush- Kuhn-Tucker (KKT) conditions

$\begin{displaymath} \nabla f_0(x) + D\tilde f(x)^T z_\mathrm{nl} + G^T z_\mat... ... k=1,\ldots,m, \qquad Gx + s_\mathrm{l} = h, \qquad Ax = b, \end{displaymath}$

where $\tilde f = (f_1,\ldots, f_m)$ ,

$\begin{displaymath} s_\mathrm{nl}\succeq 0, \qquad s_\mathrm{l}\succeq 0, \qquad... ...\mathrm{nl}^T z_\mathrm{nl} + s_\mathrm{l}^T z_\mathrm{l} = 0. \end{displaymath}$

'unknown'

This indicates that the algorithm reached the maximum number of iterations before a solution was found. The 'x', 'snl', 'sl', 'y', 'znl' and 'zl' entries are None.

Example: equality constrained analytic centering

The equality constrained analytic centering problem is defined as

$\begin{displaymath} \begin{array}{ll} \mbox{minimize} & -\sum_{i=1}^m \log x_i \\ \mbox{subject to} & Ax = b. \end{array}\end{displaymath}$

The function acent() defined below solves the problem, assumping it is solvable.

from cvxopt import solvers 
from cvxopt.base import matrix, spmatrix, log

def acent(A, b):
    m, n = A.size
    def F(x=None, z=None):
        if x is None: return 0, matrix(1.0, (n,1))
        if min(x) <= 0.0: return None
        f = -sum(log(x))
        Df = -(x**-1).T 
        if z is None: return f, Df
        H = z[0] * spmatrix(x**-2, range(n), range(n))
        return f, Df, H
    return solvers.cp(F, A=A, b=b)['x']

Example: robust least-squares

The function robls() defined below solves the unconstrained problem

$\begin{displaymath} \begin{array}{ll} \mbox{minimize} & \sum_{k=1}^m \phi((Ax-b)... ...{\mbox{\bf R}}^{m\times n}, \quad \phi(u) = \sqrt{\rho + u^2}. \end{displaymath}$

from cvxopt import solvers 
from cvxopt.base import matrix, spmatrix, sqrt, div

def robls(A, b, rho): 
    m, n = A.size
    def F(x=None, z=None):
        if x is None: return 0, matrix(0.0, (n,1))
        y = A*x-b
        w = sqrt(rho + y**2)
        f = sum(w)
        Df = div(y, w).T * A 
        if z is None: return f, Df 
        H = A.T * spmatrix(z[0]*rho*(w**-3), range(m), range(m)) * A
        return f, Df, H
    return solvers.cp(F)['x']

Example: floor planning

This example is the floor planning problem of section 8.8.2 in the book Convex Optimization:

$\begin{displaymath} \begin{array}{ll} \mbox{minimize} & W + H \\ \mbox{subjec... ...gamma \leq w_k \leq \gamma h_k, \quad k=1,\ldots,5. \end{array}\end{displaymath}$

This problem has 22 variables

$\begin{displaymath} W, \qquad H, \qquad x\in{\mbox{\bf R}}^5, \qquad y\in{\mbox{... ...}^5, \qquad w\in{\mbox{\bf R}}^5, \qquad h\in{\mbox{\bf R}}^5, \end{displaymath}$

5 nonlinear inequality constraints, and 26 linear inequality constraints. The code belows defines a function floorplan() that solves the problem by calling cp(), then applies it to 4 instances, and creates a figure.

import pylab
from cvxopt import solvers
from cvxopt.base import matrix, spmatrix, mul, div

def floorplan(Amin):

    #     minimize    W+H
    #     subject to  Amink / hk <= wk, k = 1,..., 5 
    #                 x1 >= 0,  x2 >= 0, x4 >= 0
    #                 x1 + w1 + rho <= x3  
    #                 x2 + w2 + rho <= x3 
    #                 x3 + w3 + rho <= x5  
    #                 x4 + w4 + rho <= x5
    #                 x5 + w5 <= W
    #                 y2 >= 0,  y3 >= 0,  y5 >= 0 
    #                 y2 + h2 + rho <= y1 
    #                 y1 + h1 + rho <= y4 
    #                 y3 + h3 + rho <= y4
    #                 y4 + h4 <= H  
    #                 y5 + h5 <= H
    #                 hk/gamma <= wk <= gamma*hk,  k = 1, ..., 5
    #
    # 22 Variables W, H, x (5), y (5), w (5), h (5).
    #
    # W, H:  scalars; bounding box width and height
    # x, y:  5-vectors; coordinates of bottom left corners of blocks
    # w, h:  5-vectors; widths and heigths of the 5 blocks

    rho, gamma = 1.0, 5.0   # min spacing, min aspect ratio

    # The objective is to minimize W + H.  There are five nonlinear 
    # constraints 
    #
    #     -wk + Amink / hk <= 0,  k = 1, ..., 5

    def F(x=None, z=None):
        if x is None:  return 5, matrix(17*[0.0] + 5*[1.0])
        if min(x[17:]) <= 0.0:  return None 
        f = matrix(0.0, (6,1))
        f[0] = x[0] + x[1]  
        f[1:] = -x[12:17] + div(Amin, x[17:]) 
        Df = matrix(0.0, (6,22))
        Df[0, [0,1]] = 1.0
        Df[1:,12:17] = spmatrix(-1.0, range(5), range(5))
        Df[1:,17:] = spmatrix(-div(Amin, x[17:]**2), range(5), range(5))
        if z is None: return f, Df
        H = spmatrix( 2.0* mul(z[1:], div(Amin, x[17::]**3)), range(17,22), range(17,22) )
        return f, Df, H

    G = matrix(0.0, (26,22)) 
    h = matrix(0.0, (26,1))
    G[0,2] = -1.0                                       # -x1 <= 0
    G[1,3] = -1.0                                       # -x2 <= 0 
    G[2,5] = -1.0                                       # -x4 <= 0
    G[3, [2, 4, 12]], h[3] = [1.0, -1.0, 1.0], -rho     # x1 - x3 + w1 <= -rho 
    G[4, [3, 4, 13]], h[4] = [1.0, -1.0, 1.0], -rho     # x2 - x3 + w2 <= -rho 
    G[5, [4, 6, 14]], h[5] = [1.0, -1.0, 1.0], -rho     # x3 - x5 + w3 <= -rho 
    G[6, [5, 6, 15]], h[6] = [1.0, -1.0, 1.0], -rho     # x4 - x5 + w4 <= -rho 
    G[7, [0, 6, 16]] = -1.0, 1.0, 1.0                   # -W + x5 + w5 <= 0
    G[8,8] = -1.0                                       # -y2 <= 0 
    G[9,9] = -1.0                                       # -y3 <= 0 
    G[10,11] = -1.0                                     # -y5 <= 0 
    G[11, [7, 8, 18]], h[11] = [-1.0, 1.0, 1.0], -rho   # -y1 + y2 + h2 <= -rho 
    G[12, [7, 10, 17]], h[12] = [1.0, -1.0, 1.0], -rho  #  y1 - y4 + h1 <= -rho 
    G[13, [9, 10, 19]], h[13] = [1.0, -1.0, 1.0], -rho  #  y3 - y4 + h3 <= -rho 
    G[14, [1, 10, 20]] = -1.0, 1.0, 1.0                 # -H + y4 + h4 <= 0  
    G[15, [1, 11, 21]] = -1.0, 1.0, 1.0                 # -H + y5 + h5 <= 0
    G[16, [12, 17]] = -1.0, 1.0/gamma                   # -w1 + h1/gamma <= 0 
    G[17, [12, 17]] = 1.0, -gamma                       #  w1 - gamma * h1 <= 0
    G[18, [13, 18]] = -1.0, 1.0/gamma                   # -w2 + h2/gamma <= 0 
    G[19, [13, 18]] = 1.0, -gamma                       #  w2 - gamma * h2 <= 0
    G[20, [14, 18]] = -1.0, 1.0/gamma                   # -w3 + h3/gamma <= 0  
    G[21, [14, 19]] = 1.0, -gamma                       #  w3 - gamma * h3 <= 0
    G[22, [15, 19]] = -1.0, 1.0/gamma                   # -w4  + h4/gamma <= 0 
    G[23, [15, 20]] = 1.0, -gamma                       #  w4 - gamma * h4 <= 0
    G[24, [16, 21]] = -1.0, 1.0/gamma                   # -w5 + h5/gamma <= 0 
    G[25, [16, 21]] = 1.0, -gamma                       #  w5 - gamma * h5 <= 0.0

    # solve and return W, H, x, y, w, h 
    sol = solvers.cp(F, G, h)  
    return  sol['x'][0], sol['x'][1], sol['x'][2:7], sol['x'][7:12], sol['x'][12:17], sol['x'][17:] 

pylab.figure(facecolor='w')
pylab.subplot(221)
Amin = matrix([100., 100., 100., 100., 100.])
W, H, x, y, w, h =  floorplan(Amin)
for k in xrange(5):
    pylab.fill([x[k], x[k], x[k]+w[k], x[k]+w[k]], 
               [y[k], y[k]+h[k], y[k]+h[k], y[k]], '#D0D0D0')
    pylab.text(x[k]+.5*w[k], y[k]+.5*h[k], "%d" %(k+1))
pylab.axis([-1.0, 26, -1.0, 26])
pylab.xticks([])
pylab.yticks([])

pylab.subplot(222)
Amin = matrix([20., 50., 80., 150., 200.])
W, H, x, y, w, h =  floorplan(Amin)
for k in xrange(5):
    pylab.fill([x[k], x[k], x[k]+w[k], x[k]+w[k]], 
               [y[k], y[k]+h[k], y[k]+h[k], y[k]], '#D0D0D0')
    pylab.text(x[k]+.5*w[k], y[k]+.5*h[k], "%d" %(k+1))
pylab.axis([-1.0, 26, -1.0, 26])
pylab.xticks([])
pylab.yticks([])

pylab.subplot(223)
Amin = matrix([180., 80., 80., 80., 80.])
W, H, x, y, w, h =  floorplan(Amin)
for k in xrange(5):
    pylab.fill([x[k], x[k], x[k]+w[k], x[k]+w[k]], 
               [y[k], y[k]+h[k], y[k]+h[k], y[k]],'#D0D0D0')
    pylab.text(x[k]+.5*w[k], y[k]+.5*h[k], "%d" %(k+1))
pylab.axis([-1.0, 26, -1.0, 26])
pylab.xticks([])
pylab.yticks([])

pylab.subplot(224)
Amin = matrix([20., 150., 20., 200., 110.])
W, H, x, y, w, h =  floorplan(Amin)
for k in xrange(5):
    pylab.fill([x[k], x[k], x[k]+w[k], x[k]+w[k]], 
               [y[k], y[k]+h[k], y[k]+h[k], y[k]],'#D0D0D0')
    pylab.text(x[k]+.5*w[k], y[k]+.5*h[k], "%d" %(k+1))
pylab.axis([-1.0, 26, -1.0, 26])
pylab.xticks([])
pylab.yticks([])

pylab.show()

$\includegraphics[width=15cm]{figures/floorplan.eps}$