% 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}