8.4 Semidefinite Programming

We use the following notation for a pair of primal and dual semidefinite programs (SDPs):

$\begin{displaymath} \mbox{Primal:}\quad \begin{array}[t]{ll} \mbox{minimize} & c... ...\mathrm{l} \succeq 0, \quad Z_\mathrm{s} \succeq 0. \end{array}\end{displaymath}$

(8.2)

The dimensions of the primal and dual variables are

$\begin{displaymath} x\in {\mbox{\bf R}}^n, \qquad s_\mathrm{l} \in {\mbox{\bf R... ...{\bf S}}^{m_0} \times \cdots \times {\mbox{\bf S}}^{m_{N-1}}, \end{displaymath}$

where ${\mbox{\bf S}}^n$ is the set of real symmetric matrices of order n. The problem data are the matrices

$\begin{displaymath} c\in{\mbox{\bf R}}^n, \qquad G_\mathrm{l} \in{\mbox{\bf R}}^... ... {\mbox{\bf R}}^{p\times n}, \qquad b \in {\mbox{\bf R}}^{p}, \end{displaymath}$

and the linear mapping $G_\mathrm{s} : {\mbox{\bf R}}^n \rightarrow {\mbox{\bf S}}^{m_1} \times \cdots \times {\mbox{\bf S}}^{m_N}$ and its adjoint $G_\mathrm{s}^T$ . As for LPs we store vector variables as dense real matrices with one column. Block-diagonal symmetric matrices are stored as lists of square dense real matrices, with the lower triangular part of each matrix representing the lower triangular part of a diagonal block. Entries above the diagonal are not referenced.

sdp( c[, Gl, hl[, Gs, hs[, A, b[, solver[, primalstart[, dualstart]]]]]])

Solves the pair of primal and dual SDPs (8.2).

c is a dense real matrix with one column. Gl and A are dense or sparse real matrices. hl and b are dense real matrices with one column. The default values for Gl, hl, A and b are empty matrices, i.e., matrices with zero rows.

Gs and hs are lists of length N that specify the linear matrix inequality constraints. hs is a list of square dense real matrices hs[k] of order m_k. Gs is a list of dense or sparse real matrices Gs[k] with m_k*m_k rows and n columns, such that the product Gs[k]*x is the kth diagonal block of Gs(x), stored columnwise.

The solver argument is used to choose between two solvers: the default solver (used when solver is absent or equal to None) and the external solver DSDP5 (solver='dsdp'); see section 8.8. With the 'dsdp' option the code does not accept problems with equality constraints.

The optional argument primalstart is a dictionary with keys 'x', 'sl', and 'ss', used as an optional primal starting point. dualstart is a dictionary with keys 'y', 'zl', 'zs', used as an optional dual starting point. These two arguments are ignored when the DSDP solver is used.

sdp() returns a dictionary with keys 'status', 'x', 'sl', 'ss', 'y', 'zl', 'ss'. The possible values of the 'status' item are as follows.

'optimal'.

In this case the 'x', 'sl', 'ss', 'y', 'zl', 'zs' entries contain primal and dual optimal solutions, which approximately satisfy

$\begin{displaymath} G_\mathrm{l}x+s_\mathrm{l} = h_\mathrm{l}, \qquad G_\mathr... ...m{l}^Tz_\mathrm{l} + G_\mathrm{s}^T(Z_\mathrm{s}) + A^Ty+c = 0 \end{displaymath}$

and

$\begin{displaymath} s_\mathrm{l} \succeq 0, \qquad S_\mathrm{s} \succeq 0, \qqu... ... z_\mathrm{l} + \mathop{\bf tr}(S_\mathrm{s}Z_\mathrm{s}) = 0. \end{displaymath}$

'primal infeasible'.

The 'x', 'sl' and 'ss' entries are None, and the 'y', 'zl', 'zs' entries provide an approximate certificate of infeasibility:

$\begin{displaymath} G_\mathrm{l}^Tz_\mathrm{l} + G_\mathrm{s}^T(Z_\mathrm{s}) +A... ...\qquad z_\mathrm{l} \succeq 0, \qquad Z_\mathrm{s} \succeq 0. \end{displaymath}$

'dual infeasible'.

The SDP is dual infeasible. The 'y', 'zl' and 'zs' entries are None, and the 'x', 'sl', 'ss' entries contain an approximate certificate of dual infeasibility:

$\begin{displaymath} G_\mathrm{l}x+s_\mathrm{l} = 0, \qquad G_\mathrm{s}(x)+S_\ma... ...\qquad s_\mathrm{l} \succeq 0, \qquad S_\mathrm{s} \succeq 0. \end{displaymath}$

'unknown'.

The 'x', 'sl', 'ss', 'y', 'zl' and 'zs' entries are None.

We illustrate the calling sequence with a small example.

$\begin{displaymath} \begin{array}{ll} \mbox{minimize} & x_1 - x_2 + x_3 \\ \mbo... ... 9 & 91 & 10 \\ 40 & 10 & 15 \end{array} \right] \end{array}\end{displaymath}$

>>> from cvxopt.base import matrix
>>> from cvxopt import solvers
>>> c = matrix([1.,-1.,1.])
>>> G = [ matrix([[-7., -11., -11., 3.], 
                  [ 7., -18., -18., 8.], 
                  [-2.,  -8.,  -8., 1.]]) ]
>>> G += [ matrix([[-21., -11.,   0., -11.,  10.,   8.,   0.,   8., 5.], 
                   [ 0.,  10.,  16.,  10., -10., -10.,  16., -10., 3.], 
                   [ -5.,   2., -17.,   2.,  -6.,   8., -17.,  -7., 6.]]) ]
>>> h = [ matrix([[33., -9.], [-9., 26.]]) ]
>>> h += [ matrix([[14., 9., 40.], [9., 91., 10.], [40., 10., 15.]]) ]
>>> sol = solvers.sdp(c, Gs=G, hs=h)  
>>> print sol['x']
  -3.6775e-01
   1.8983e+00
  -8.8747e-01
>>> print sol['zs'][0]
   3.9613e-03   0.0000e+00
  -4.3390e-03   4.7526e-03
>>> print sol['zs'][1]
   5.5803e-02   0.0000e+00   0.0000e+00
  -2.4103e-03   1.0411e-04   0.0000e+00
   2.4214e-02  -1.0459e-03   1.0507e-02

Note that only the lower triangular parts of the dual variables are returned (in the example the returned values of the upper triangular elements happen to be zero, but this is not necessarily the case).

Only the entries in Gs and hs that correspond to lower triangular entries need to be provided, so in the example h and G can also be defined as follows.

>>> G = [ matrix([[-7., -11., 0., 3.], 
                  [ 7., -18., 0., 8.], 
                  [-2.,  -8., 0., 1.]]) ]
>>> G += [ matrix([[-21., -11.,   0., 0.,  10.,   8., 0., 0., 5.], 
                   [  0.,  10.,  16., 0., -10., -10., 0., 0., 3.], 
                   [ -5.,   2., -17., 0.,  -6.,   8., 0., 0., 6.]]) ]
>>> h = [ matrix([[33., -9.], [0., 26.]]) ]
>>> h += [ matrix([[14., 9., 40.], [0., 91., 10.], [0., 0., 15.]]) ]