% semidef. relaxations in matlab

\documentstyle[11pt]{article}
\topmargin 0pt
\textheight 8in
\textwidth 5.4in

\begin{document}

\newcommand{\eps}{\varepsilon}
\newcommand{\diag}{\mbox{diag}}
\newcommand{\Diag}{\mbox{Diag}}
\newcommand{\tr}{\em{tr}}
\newcommand{\ms}{\medskip}
\newcommand{\bs}{\bigskip}
\newcommand{\opa}{{\cal A}}
\renewcommand{\d}{\delta}

\title{
A Matlab toolbox for semidefinite programming}

\author{
Franz Rendl\\
Technische Universit\"at Graz\\
Institut f\"ur Mathematik\\
Kopernikusgasse 24, A-8010 Graz, Austria \\
e-mail: rendl@fmatbds01.tu-graz.ac.at
}

\maketitle


\begin{abstract}
Semidefinite programs have found increasing
interest in the recent past, see \cite{AHO,HRVW}.
This note describes an
easy to use toolbox in Matlab to solve semidefinite programs
by interior point methods.
The user has to provide the data of the problem and write a
routine that provides an initial feasible interior point to
the problem. Then she can use the preset routines to solve
the problem.
\end{abstract}

\section{The problem}

Let $C$ be a given symmetric $n \times n$ matrix, $a \in \Re^m$ and
$\opa~$ be a linear operator, that maps symmetric matrices of size $n$
into vectors in $\Re^m$.
We provide a Matlab toolbox to solve the following pair of primal and
dual semidefinite problems.

\medskip
{\bf primal problem.}
maximize $~\tr CX~$ such that $~a - \opa(X) = 0,~~X~ psd.$
\ms

{\bf dual problem.}
minimize $~a^ty~$ such that $~\opa^t(y) - C~ psd.$
\ms

{\bf interior point assumption.}
It is assumed throughout, that there exists $(X,y)$ such that $X$ is positive
definite, $\opa(X)=a, \opa^t(y) - C$ positive definite.
We call these points {\sl feasible interior points}.

\bs
Under this assumption strong duality holds. For the theoretical backbone of
the algorithm presented, we refer to \cite{HRVW}.
The operator $\opa~$ can be represented by $m$ symmetric matrices $A_i$
as follows.
$$\opa(X) = \left(
\begin{array}{c} \tr A_1 X \\ \vdots \\ \tr A_m X \end{array}
\right)$$
This implies that $$\opa ^t(y) = \sum_i y_i A_i.$$

\section{The algorithm}

The method used in the Matlab routines is based on \cite{HRVW}. The
reader is referred to this paper for any theoretical details.

It is convenient to introduce a slack matrix
$$Z = \opa ^t(y) -C$$
for the dual. The problems are solved using a predictor-corrector
primal-dual interior point approach. Here is the key step.
Given an interior feasible point $(X,y,Z)$, and a value
$\mu$, a new interior feasible point is obtained by solving the
following linear system, corresponding to a Newton step on the
KKT conditions at $(X,y,Z)$.
\begin{eqnarray*}
\d Z &=& \opa ^t(\d y) \\
\opa ( \d X) &=& 0 \\
Z \d X + \d Z X & =& \mu I - ZX
\end{eqnarray*}
The last condition comes from linearizing
$$(Z+ \d Z) (X + \d X) = \mu I.$$
Note that the first and third of these equations can be used to express $\d X$
and $\d Z$ in terms of $\d y$.
This system is solved using a predictor-corrector approach, that provides
an improved solution to the nonlinear third equation. Here are the details.
\begin{eqnarray*}
\d Z^p + \d Z^c & =& \opa^t( \d y^p + \d y^c)\\
\opa ( \d X^p + \d X^c)& =& 0 \\
ZX + \d Z^pX + \d Z^c X + Z \d X^p + Z \d X^c + \d Z^p \d X^p& =& \mu I.
\end{eqnarray*}
The last equation is obtained from approximating
$(Z + \d Z^p + \d Z^c)(X + \d X^p + \d X^c)$ by the term on the left hand
side. These equations are decoupled as follows.
\begin{eqnarray*}
\d Z^p & = & \opa (\d y^p) \\
\d Z^c & = &\opa (\d y^c) \\
\opa( \d X^p) & = & 0 \\
\opa( \d X^c) & = & 0 \\
\d Z^p X + Z \d X^p & = & -ZX \\
\d Z^c X + Z \d X^c & = & \mu I - \d Z^p \d X^p \\
\end{eqnarray*}
These systems can be solved first for the predictor $(\d X^p, \d y^p, \d
Z^p)$,
which in turn is used to solve for the corrector $(\d X^c, \d y^c, \d Z^c)$.

\ms
{\bf predictor system.}
After substitution one obtains:
\begin{eqnarray*}
\opa ( Z^{-1} \opa ^t(\d y ^p) X) &=& -a \\
\d Z^p & = & \opa ^t( \d y^p) \\
\d X^p &= &- X - Z^{-1} \d Z^p X
\end{eqnarray*}
This is used to find the correctors.

\bs
{\bf corrector system.}
\begin{eqnarray*}
\opa ( Z^{-1} \opa ^t(\d y^c) X) &=& \mu \opa (Z^{-1}) -\opa (Z^{-1} \d Z^p
\d X^p) \\
\d Z^c & = & \opa ^t( \d y^c) \\
\d X^c& = &\mu Z^{-1}  - Z^{-1} \d Z^p \d X^p - Z^{-1} \d Z^c X
\end{eqnarray*}

{\bf actual search direction.}
The actual search direction is now the sum of the two steps.
$$dy = \d y^p + \d y^c, ~~dX = \d X^p + \d X^c,~~ dZ = \d Z^p + \d Z^c.$$

To conclude one iteration one finds primal and dual step sizes
$\alpha_p$ and $\alpha_d$, so that the new point
$$X + \alpha_p dX,~~ y + \alpha_d dy,~~Z + \alpha_d dZ$$
is again a feasible interior point. The parameter $\mu$ is updated, and one
iterates, until the duality gap is sufficiently small.

\section{User supplied data}

The user has to provide the data defining the problem, and write
a routine, that provides an interior feasible starting point.
The operator $\opa~$ is defined through the $m$ matrices $A_i$. These are
stored in two arrays.

\begin{verbatim}
nonz(i), i = 1, ... , m
\end{verbatim}

contains the number of nonzeros of $A_i$ in the upper triangle, including
the main diagonal. The second array stores the indices and actual entries
of the nonzeros of each $A_i$. This array consists of 3 rows.

\begin{verbatim}
store( 3, :)
\end{verbatim}
A typical column of this array contains
$$ \left(
\begin{array}{c} i\\ j\\ A_k(i,j) \end{array} \right).$$
The order of the nonzeros for a specific matrix $A_k$ is arbitrary, but
it is assumed that the first $nonz(1)$ entries  correspond to $A_1$,
the next $nonz(2)$ entries to $A_2$, and so on. Thus
the number of columns of
the array $store$ is equal to the sum of the entries of the array $nonz$.

Thus an instancs of a semidefinite program is defined by the quadruple

\begin{verbatim}
C, a, store, nonz
\end{verbatim}

where $C$ is a symmetric matrix of size $n$, $a$ and $nonz$ are
column vectors of size $m$ and $store$ contains the indices of the $A_i$.

Finally the user has to provide a Matlab function, that generates an interior
feasible starting point. Here is the call to this function.
\begin{verbatim}
[X, y, Z]  = psd_strt( C, a, store, nonz)
\end{verbatim}
Once this is done the problem can be solved by  calling the matlab function
$psd.$
\begin{verbatim}
[X, y, psi] = psd( C, a, store, nonz);
\end{verbatim}

\section{The Matlab functions}

The toolbox consists presently of the following Matlab functions, which
need not be changed by the user.
\begin{itemize}

\item psd.m: This is the main routine. It has the following calling sequence.
On successful termination, $X$ is primal optimal, $y$ is dual optimal, and
$psi$ is the optimal value of the primal.
\begin{verbatim}
[X, y, psi] = psd( C, a, store, nonz);
\end{verbatim}

\item op\_a.m: This function evaluates $\opa(X)$. Calling sequence:
\begin{verbatim}
vector = op_a( X, store, nonz);
\end{verbatim}

\item op\_at.m: This function evaluates $\opa^t(y)$. Calling sequence:
Note that here the primal matrix size $n$ is a fourth input argument.
\begin{verbatim}
matrix = op_at( y, store, nonz, n);
\end{verbatim}

\item op\_biga.m: This function evaluates the system matrix corresponding
to the operator $\opa( Z^{-1} \opa^t(.) X)$.
Calling sequence:
\begin{verbatim}
matrix = op_biga( X, Zinv, store, nonz);
\end{verbatim}

\item psd\_feas.m: This function checks strict feasibility
of a point $(X, y, Z)$.
Calling sequence:
\begin{verbatim}
okay = psd_feas(X, y, Z, C, a, store, nonz);
\end{verbatim}
If $okay=0$ then the solution is strictly feasible, within the tolerances
set in the program. Otherwise a message is printed, that indicates which of
the conditions is violated.

\item pd\_check.m: This function checks whether a given matrix is positive
definite or not.

\item pd\_solve.m: This function solves a positive definite system, once the
cholesky factor is given.
\end{itemize}



\section{An example}

Consider the following artificial example.

$$\max \tr CX ~~~\mbox{such that}$$
$$x_{ii} + x_{i+1, i+1} - x_{i,i+1} - x_{i+1,i} = 2,~~~i=1, \ldots, n-1$$
$$\tr(X) = n,$$
$$X \mbox{~~positive semidefinite.}$$
Here there are $n$ constraints. The first $n-1$ constraints
are zero matrices with a 2 by 2 block
of the form
$$ \left( \begin{array}{r r} 1& -1 \\ -1& 1 \end{array} \right)$$
on the main diagonal starting in row and column $i$.
The last matrix $A_n = I_n$, the identity matrix of size $n$.

Suppose, $n=5$ and $C$ is the matrix of all ones.
Note that $X=I_n$ is a feasible interior primal point and
$y_i=0, ~~i=1, \ldots, n-1, ~~y_n = n+1$ leads to an interior dual solution.

Here is the matlab function that defines the starting point.

{\def\baselinestretch{0.8}\large\small
\begin{verbatim}
function [ X, y, Z] = psd_strt( C, a, store, nonz);
  [n, n1] = size( C);
  X = eye(n);
  y = zeros(n,1);
  y(n)=max(sum(abs(C)))*1.1;
  Z = op_at( y, store, nonz, n) - C;
\end{verbatim} }

Note that the definition of $Z$ is done using the supplied function
$op\_at$, which evaluates $\opa ^t(y)$.

The data can be set as follows.
{\def\baselinestretch{0.8}\large\small
\begin{verbatim}
n=5;
a = 2 * ones(n,1);
a(n) = n;
C = ones( n);
nonz = 3*ones(n,1);
nonz(n) =n;
for i = 1:n-1
store(:,3*i-2:3*i)= [i i i+1
                     i i+1 i+1
                     1 -1  1];
                     end;
store(:, 3*(n-1)+1: 3*(n-1)+n)=[1:n; 1:n; ones(1,n)];
\end{verbatim}  }
This is contained in the m-file {psd\_expl.m}. Here is a sample
output.

{\def\baselinestretch{0.8}\large\small
\begin{verbatim}
[X,y,psi]=psd(C,a,store,nonz)
   iter    alphap    alphad   log(gap)  primal     dual     unsym
      0    1.0000    1.0000    1.3522    5.0000   27.5000
 1.0000    0.7600    1.0000    1.1093   10.1044   22.9663  -99.0000
 2.0000    1.0000    1.0000    0.5072   11.6531   14.8686  -15.3525
 3.0000    1.0000    1.0000   -0.3958   12.7525   13.1545  -14.6122
 4.0000    1.0000    1.0000   -1.2989   12.9722   13.0224  -13.2453
 5.0000    1.0000    1.0000   -2.2020   12.9963   13.0025  -12.7682
 6.0000    1.0000    1.0000   -3.1051   12.9995   13.0003  -12.0412
 7.0000    1.0000    1.0000   -4.0082   12.9999   13.0000  -10.9620
 8.0000    1.0000    1.0000   -4.9113   13.0000   13.0000   -9.9340
 9.0000    1.0000    1.0000   -5.8144   13.0000   13.0000   -9.0309
10.0000    1.0000    1.0000   -6.7175   13.0000   13.0000   -8.1278
elapsed time:
 4.9400
X =
 0.8400    0.0400    0.8400    0.0400    0.8400
 0.0400    1.2400    0.0400    1.2400    0.0400
 0.8400    0.0400    0.8400    0.0400    0.8400
 0.0400    1.2400    0.0400    1.2400    0.0400
 0.8400    0.0400    0.8400    0.0400    0.8400
y =
-2.0000
-1.0000
-1.0000
-2.0000
 5.0000
psi =
13.0000
\end{verbatim}}

\section{Comments}

There are several observations to make..
\begin{enumerate}
\item I do not take any {\bf responsibility} or {\bf garantuee}
whatsoever on the performance of the code. It works well on the problems that
I tested it, and I would encourage any feedback on the code.
\item I realize that it is limiting to have the user provide an initial
strictly interior solution. This may be handled automatically in the update
to this version.
\item It is possible to include linear inequalities into the model. This
feature definitely will be included in the next update.
\item For specific problems, it may be advantagous to rewrite the function
files describing the action of the operator $\opa$. For instance in the
max-cut case, one has
$$\opa(X) = \diag(X),$$
and the system matrix simply is
$$X .* Z^{-1},$$
i.e. the elementwise product of the two matrices.
\item This version does not exploit any features of version 4 of Matlab.
\end{enumerate}


\bibliographystyle{plain}
\begin{thebibliography}{10}

\bibitem{AHO}
F.~ALIZADEH, J.A.~HAEBERLY and M.L.~OVERTON.
\newblock A new primal-dual interior point method for semidefinite
programming.
\newblock Technical Report, August 1994.

\bibitem{HRVW}
C.~HELMBERG, F.~RENDL, R.~J. VANDERBEI, and H.~WOLKOWICZ.
\newblock An interior-point method for semidefinite programming.
\newblock Technical Report CORR 93-20, SOR 93-15, Department of Combinatorics
  and Optimization, Waterloo, Ont, 1993.
\newblock to appear in: SIAM optimization.

\end{thebibliography}
\end{document}

