=============================
Semidefinite Programming
=============================
A *positive semidefinite* matrix is a symmetric real matrix whose eigenvalues are all nonnegative.
A semidefinite programming problem is simply a linear program where the solutions are positive
semidefinite matrices instead of points in Euclidean space.
Primal and Dual formulations
=============================
A typical semidefinite program (SDP for short) in the primal form is the following optimization problem:
.. math::
\left\lbrace
\begin{array}{lll}
\min & \sum_{i=1}^m b_i x_i & \\
\textrm{subject to} & & \\
& \sum_{i=1}^m A_{ij}x_i - C_j \succeq 0 & j=1,\dots,k.
\end{array}\right.
The dual program associated to the above SDP will be the following:
.. math::
\left\lbrace
\begin{array}{lll}
\max & \sum_{j=1}^k tr(C_j\times Z_j) & \\
\textrm{subject to} & & \\
& \sum_{j=1}^k tr(A_{ij}\times Z_j) = b_i & i=1,\dots,m,\\
& Z_j \succeq 0 & j=1,\dots,k.
\end{array}\right.
For convenience, we use a block representation for the matrices as follows:
.. math::
C = \left(
\begin{array}{cccc}
C_1 & 0 & 0 & \dots \\
0 & C_2 & 0 & \dots \\
\vdots & \dots & \ddots & \vdots \\
0 & \dots & 0 & C_k
\end{array}
\right),
and
.. math::
A_i = \left(
\begin{array}{cccc}
A_{i1} & 0 & 0 & \dots \\
0 & A_{i2} & 0 & \dots \\
\vdots & \dots & \ddots & \vdots \\
0 & \dots & 0 & A_{ik}
\end{array}
\right).
This simplifies the :math:`k` constraints of the primal form in to one constraint
:math:`\sum_{i=1}^m A_i x_i - C \succeq 0` and the objective and constraints of the
dual form as :math:`tr(C\times Y)` and :math:`tr(A_i\times Z_i) = b_i` for :math:`i=1,\dots,m`.
The ``sdp`` class
=============================
The ``sdp`` class provides an interface to solve semidefinite programs using various range of
well-known SDP solvers. Currently, the following solvers are supported:
``CVXOPT``
----------------------------
This is a python native convex optimization solver which can be obtained from `CVXOPT `_.
Beside semidefinite programs, it has various other solvers to handle convex optimization problems.
In order to use this solver, the python package ``CVXOPT`` must be installed.
``DSDP``
----------------------------
If `DSDP `_ and ``CVXOPT`` are installed and ``DSDP`` is callable from command line,
then it can be used as a SDP solver. Note that the current implementation uses ``CVXOPT`` to call ``DSDP``, so ``CVXOPT`` is a
requirement too.
``SDPA``
----------------------------
In case one manages to install `SDPA `_ and it can be called from command line, one can use
``SDPA`` as a SDP solver.
``CSDP``
----------------------------
Also, if `csdp `_ is installed and can be reached from command, then it can be used to solve
SDP problems through ``sdp`` class.
To initialize and set the solver to one of the above simply use::
SDP = sdp('cvxopt') # initializes and uses `cvxopt` as solver.
.. note::
In windows, one can provide the path to each of the above solvers as the second parameter of the constructor::
SDP = sdp('csdp', {'csdp':"Path to executable csdp"}) # initializes and uses `csdp` as solver existing at the given path.
Set the :math:`b` vector:
----------------------------
To set the vector :math:`b=(b_1,\dots,b_m)` one should use the method ``sdp.SetObjective`` which takes a list or a numpy array of
numbers as :math:`b`.
Set a block constraint:
----------------------------
To introduce the block of matrices :math:`A_{i1},\dots, A_{ik}` associated with :math:`x_i`, one should use the method
``sdp.AddConstraintBlock`` that takes a list of matrices as blocks.
Set the constant block `C`:
----------------------------
The method ``sdp.AddConstantBlock`` takes a list of square matrices and use them to construct :math:`C`.
Solve the input SDP:
----------------------------
To solve the input SDP simply call the method ``sdp.solve()``. This will call the selected solver on the entered SDP and
the output of the solver will be set as dictionary in ``sdp.Info`` with the following keys:
+ ``PObj``: The value of the primal objective.
+ ``DObj``: The value of the dual objective.
+ ``X``: The final :math:`X` matrix.
+ ``Z``: The final :math:`Z` matrix.
+ ``Status``: The final status of the solver.
+ ``CPU``: Total run time of the solver.
Example:
----------------------------
Consider the following SDP:
.. math::
\left\lbrace
\begin{array}{lll}
\min & x_1 - x_2 + x_3 \\
\textrm{subject to} & \\
& \left(\begin{array}{cc}7 & 11\\ 11 & -3 \end{array}\right)x_1 +
\left(\begin{array}{cc}-7 & 18\\ 18 & -8 \end{array}\right)x_2 +
\left(\begin{array}{cc} 2 & 8\\ 8 & -1 \end{array}\right)x_3
\succeq\left(\begin{array}{cc} -33 & 9\\ 9 & -26 \end{array}\right) \\
& \left(\begin{array}{ccc}21 & 11 & 0\\ 11 & -10 & -8\\ 0 & -8 & -5\end{array}\right)x_1 +
\left(\begin{array}{ccc}0 & -10 & -16\\ -10 & 10 & 10\\ -16 & 10 & -3\end{array}\right)x_2 +
\left(\begin{array}{ccc} 5 & -2 & 17\\ -2 & 6 & -8\\ 17 & -8 & -6\end{array}\right)x_3
\succeq\left(\begin{array}{ccc} -14 & -9 & -40\\ -9 & -91 & -10\\ -40 & -10 & -15\end{array}\right) \\
\end{array}
\right.
The following code solves the above program::
from numpy import matrix
from Irene import sdp
b = [1, -1, 1]
C = [matrix([[-33, 9], [9, -26]]),
matrix([[-14, -9, -40], [-9, -91, -10], [-40, -10, -15]])]
A1 = [matrix([[7, 11], [11, -3]]),
matrix([[21, 11, 0], [11, -10, -8], [0, -8, -5]])]
A2 = [matrix([[-7, 18], [18, -8]]),
matrix([[0, -10, -16], [-10, 10, 10], [-16, 10, -3]])]
A3 = [matrix([[2, 8], [8, -1]]),
matrix([[5, -2, 17], [-2, 6, -8], [17, -8, -6]])]
SDP = sdp('cvxopt')
SDP.SetObjective(b)
SDP.AddConstantBlock(C)
SDP.AddConstraintBlock(A1)
SDP.AddConstraintBlock(A2)
SDP.AddConstraintBlock(A3)
SDP.solve()
print SDP.Info