7.1 Linear Cone Programs

conelp(c, G, h[, dims[, A, b[, primalstart[, dualstart[, kktsolver]]]]])

Solves a pair of primal and dual cone programs

minimize   cTx                        maximize  - hT z - bTy
subject to Gx + s = h                 subject to GTz + ATy + c = 0
          Ax = b                               z ≽ 0.
          s ≽ 0

The primal variables are x and s. The dual variables are y and z. The inequalities are interpreted as s ∈ C, z ∈ C, where C is a cone defined as a Cartesian product of a nonnegative orthant, a number of second-order cones, and a number of positive semidefinite cones:

C = C0 × C1 × ⋅⋅⋅× CM × CM+1 × ⋅⋅⋅× CM+N

with

                                                                                                     {             }
C0 = {u ∈ Rl | uk ≥ 0, k = 1,...,l}, Ck+1 = {(u0,u1) ∈ R×Rrk -1 | u0 ≥ ∥u1∥2}, k = 0,...,M - 1, Ck+M+1 = vec(u) | u ∈ Stk+ , k = 0,...,N - 1.

In this definition, vec(u) denotes a symmetric matrix u stored as a vector in column major order. The structure of C is specified by dims. This argument is a dictionary with three fields.

dims[’l’]:
l, the dimension of the nonnegative orthant (a nonnegative integer).
dims[’q’]:
[r0,,rM-1], a list with the dimensions of the second-order cones (positive integers).
dims[’s’]:
[t0,,tN-1], a list with the dimensions of the positive semidefinite cones (nonnegative integers).

The default value of dims is {’l’: G.size[0], ’q’: [], ’s’: []}, i.e., by default the inequality is interpreted as a componentwise vector inequality.

The arguments c, h and b are real single-column dense matrices. G and A are real dense or sparse matrices. The number of rows of G and h is equal to

       M∑-1    N∑- 12
K = l+     rk +    tk.
        k=0     k=0

The columns of G and h are vectors in

                         2         2
Rl × Rr0 × ⋅⋅⋅× RrM -1 × Rt0 × ⋅⋅⋅× RtN-1,

where the last N components represent symmetric matrices stored in column major order. The strictly upper triangular entries of these matrices are not accessed (i.e., the symmetric matrices are stored in the ’L’-type column major order used in the blas and lapack modules). The default values for A and b are matrices with zero rows, meaning that there are no equality constraints.

primalstart is a dictionary with keys ’x’ and ’s’, used as an optional primal starting point. primalstart[’x’] and primalstart[’s’] are real dense matrices of size (n,1) and (K,1), respectively, where n is the length of c. The vector primalstart[’s’] must be strictly positive with respect to the cone C.

dualstart is a dictionary with keys ’y’ and ’z’, used as an optional dual starting point. dualstart[’y’] and dualstart[’z’] are real dense matrices of size (p,1) and (K,1), respectively, where p is the number of rows in A. The vector dualstart[’s’] must be strictly positive with respect to the cone C.

The role of the optional argument kktsolver is explained in section 7.7.

conelp() returns a dictionary that contains the result and information about the accuracy of the solution. The most important fields have keys ’status’, ’x’, ’s’, ’y’, ’z’. The ’status’ field is a string with possible values ’optimal’, ’primal infeasible’, ’dual infeasible’ and ’unknown’. The meaning of the ’x’, ’s’, ’y’, ’z’ fields depends on the value of ’status’.

’optimal’
In this case the ’x’, ’s’, ’y’ and ’z’ entries contain the primal and dual solutions, which approximately satisfy
                         T    T                               T
Gx+s = h,    Ax = b,    G z+A  y+c = 0,    s ≽ 0,   z ≽ 0,   s z = 0.

The other entries in the output dictionary summarize the accuracy with which these optimality conditions are satisfied. The fields ’primal objective’, ’dual objective’, and ’gap’ give the primal objective cTx, dual objective -hTz - bTy, and the gap sTz. The field ’relative gap’ is the relative gap, defined as

         sTz                     T     T     T
max-{- cTx,--hTz --bTy} if max{- c x,- h z - b y} > 0

and None otherwise. The fields ’primal infeasibility’ and ’dual infeasibility’ are the residuals in the primal and dual equality constraints, defined as

     ∥Gx + s- h∥   ∥Ax - b∥      ∥GT z +AT y+ c∥
max {max-{1,∥h∥} ,max-{1,∥b∥} },  ---max{1,∥c∥}--,

respectively.

’primal infeasible’
The ’x’ and ’s’ entries are None, and the ’y’, ’z’ entries provide an approximate certificate of infeasibility, i.e., vectors that approximately satisfy
  T     T          T    T
G  z + A y = 0,   h z +b y = - 1,   z ≽ 0.

The field ’residual as primal infeasibility certificate’ gives the residual

  T     T
∥G-z-+-A-y∥.
max {1,∥c∥}

’dual infeasible’
The ’y’ and ’z’ entries are None, and the ’x’ and ’s’ entries contain an approximate certificate of dual infeasibility
Gx + s = 0,   Ax = 0,    cTx = - 1,   s ≽ 0.

The field ’residual as dual infeasibility certificate’ gives the residual

    --∥Gx-+s∥-- ---∥Ax∥---
max{max {1,∥h∥},max{1,∥b∥}}.

’unknown’
This indicates that the algorithm terminated early due to numerical difficulties or because the maximum number of iterations was reached. The ’x’, ’s’, ’y’, ’z’ entries contain the iterates when the algorithm terminated. Whether these entries are useful, as approximate solutions or certificates of primal and dual infeasibility, can be determined from the other fields in the dictionary. The fields ’primal objective’, ’dual objective’, ’gap’, ’relative gap’, ’primal infeasibility’, ’dual infeasibility’ are defined as when ’status’ is ’optimal’. The field ’residual as primal infeasibility certificate’ is defined as
        T     T
------∥G-z-+-A-y∥------.
- (hT z + bTy)max {1,∥h∥}

if hTz + bTy < 0, and None otherwise. A small value of this residual indicates that y and z, divided by -hTz - bTy, are an approximate proof of primal infeasibility. The field ’residual as dual infeasibility certificate’ is defined as

max {---∥Gx-+-s∥----,-----∥Ax-∥-----}
     - cTx max{1,∥h∥} - cT xmax {1,∥b∥}

if cTx < 0, and as None otherwise. A small value indicates that x and s, divided by -cTx are an approximate proof of dual infeasibility.

It is required that

                   [   ]
rank(A) = p,   rank(  G  ) = n,
                     A

where p is the number or rows of A and n is the number of columns of G and A.

As an example we solve the problem

minimize   - 6x - 4x - 5x
             1     2    3
subject to 16x1 - 14x2 + 5x3 ≤ - 3
          7x1 + 2x2 ≤ 5
          (                    2                        2                     2)1∕2
           (8x1 + 13x2 - 12x3 - 2) + (- 8x1 +18x2 + 6x3 - 14) + (x1 - 3x2 - 17x3 - 13) ≤ - 24x1 - 7x2 + 15x3 +12
          (x21 + x22 + x23)1∕2 ≤ 10
          ⌊                                                      ⌋   ⌊                ⌋
          ⌈  7x1 + 3x2 + 9x3  - 5x1 + 13x2 + 6x3   x1 - 6x2 - 6x3 ⌉   ⌈  68   - 30 - 19 ⌉
            - 5x1 + 13x2 + 6x3 x1 + 12x2 - 7x3   - 7x1 - 10x2 - 7x3 ≼   - 30  99    23   .
              x1 - 6x2 - 6x3  - 7x1 - 10x2 - 7x3 - 4x1 - 28x2 - 11x3    - 19  23    10

>>> from cvxopt import matrix, solvers  
>>> c = matrix([-6., -4., -5.])  
>>> G = matrix([[ 16., 7.,  24.,  -8.,   8.,  -1.,  0., -1.,  0.,  0.,   7.,  -5.,   1.,  -5.,   1.,  -7.,   1.,   -7.,  -4.],  
                [-14., 2.,   7., -13., -18.,   3.,  0.,  0., -1.,  0.,   3.,  13.,  -6.,  13.,  12., -10.,  -6.,  -10., -28.],  
                [  5., 0., -15.,  12.,  -6.,  17.,  0.,  0.,  0., -1.,   9.,   6.,  -6.,   6.,  -7.,  -7.,  -6.,   -7., -11.]])  
>>> h = matrix( [ -3., 5.,  12.,  -2., -14., -13., 10.,  0.,  0.,  0.,  68., -30., -19., -30.,  99.,  23., -19.,   23.,  10.] )  
>>> dims = {’l’: 2, ’q’: [4, 4], ’s’: [3]}  
>>> sol = solvers.conelp(c, G, h, dims)  
>>> sol[’status’]  
’optimal’  
>>> print sol[’x’]  
[-1.22e+00]  
[ 9.66e-02]  
[ 3.58e+00]  
>>> print sol[’z’]  
[ 9.30e-02]  
[ 2.04e-08]  
[ 2.35e-01]  
[ 1.33e-01]  
[-4.74e-02]  
[ 1.88e-01]  
[ 2.79e-08]  
[ 1.85e-09]  
[-6.32e-10]  
[-7.59e-09]  
[ 1.26e-01]  
[ 8.78e-02]  
[-8.67e-02]  
[ 8.78e-02]  
[ 6.13e-02]  
[-6.06e-02]  
[-8.67e-02]  
[-6.06e-02]  
[ 5.98e-02]

Only the entries of G and h defining the lower triangular portions of the coefficients in the linear matrix inequalities are accessed. We obtain the same result if we define G and h as below.

>>> G = matrix([[ 16., 7.,  24.,  -8.,   8.,  -1.,  0., -1.,  0.,  0.,   7.,  -5.,   1.,  0.,   1.,  -7.,  0.,  0.,  -4.],  
                [-14., 2.,   7., -13., -18.,   3.,  0.,  0., -1.,  0.,   3.,  13.,  -6.,  0.,  12., -10.,  0.,  0., -28.],  
                [  5., 0., -15.,  12.,  -6.,  17.,  0.,  0.,  0., -1.,   9.,   6.,  -6.,  0.,  -7.,  -7.,  0.,  0., -11.]])  
>>> h = matrix( [ -3., 5.,  12.,  -2., -14., -13., 10.,  0.,  0.,  0.,  68., -30., -19.,  0.,  99.,  23.,  0.,  0.,  10.] )