%======================= main.tex ===========================================
\documentclass{slides}
\usepackage{latexsym}
\input epsf
\input psfig
\newtheorem{exam}{Example}
\newtheorem{prop}{Proposition}
\newtheorem{lem}{Lemma}
\newtheorem{thm}{Theorem}
\newtheorem{cor}{Corollary}
\newcommand{\1}{\;\;\; }
\newcommand{\Sn}{{\cal S}_n }
\newcommand{\Se}{{\cal S}_e }
\newcommand{\Sd}{{\cal S}_d }
\newcommand{\Sc}{{\cal S}_C }
\newcommand{\Sh}{{\cal S}_H }
\newcommand{\snn}{{\cal S}_{n-1} }
\newcommand{\KK}{{\cal K} }
\newcommand{\DD}{{\cal D} }
\newcommand{\BB}{{\cal B} }
\newcommand{\PP}{{\cal P} }
\newcommand{\TT}{{\cal T} }
\newcommand\T{{\cal T}}
\newcommand{\bpr}{{\bf Proof.} \hspace{1 em}}
\newcommand{\epr}{ \\ \hspace*{4.5in} $\Box$ }
\newcommand{\beq}{ \begin{equation} }
\newcommand{\eeq}{ \end{equation} }
\newcommand{\bt}{ \begin{tabular} }
\newcommand{\et}{ \end{tabular} }
\newcommand\A{{\cal A}}
\newcommand{\Diag}{{\rm Diag\,}}
\newcommand{\diag}{{\rm diag\,}}
\newcommand{\tr}{{\rm trace\,}}
\newcommand{\trace}{{\rm trace\,}}
\newcommand{\rank}{{\rm rank\,}}
\newcommand{\p}{{\cal P}}
\newcommand{\kvec}{{\rm vec\,}}
\newcommand{\trian}{{\rm trian\,}}
\newcommand{\Trian}{{\rm Trian\,}}
\newcommand{\BoDiag}{{\rm B^0Diag\,}}
\newcommand{\OoDiag}{{\rm O^0Diag\,}}
\newcommand{\arrow}{{\rm arrow\,}}
\newcommand{\Arrow}{{\rm Arrow\,}}
\newcommand{\bodiag}{{\rm b^0diag\,}}
\newcommand{\oodiag}{{\rm o^0diag\,}}
\newcommand{\ck}{{\cal C}_k}
\newcommand{\dk}{{\cal D}_k}
\newcommand{\uu}{{\cal U}}
\newcommand{\uk}{{\cal U}_k}
\newcommand{\vk}{{\cal V}_k}
\newcommand{\wk}{{\cal W}_k}
\newcommand{\wm}{{\cal W}_m}
\newcommand{\w}{{\cal W}}
\newcommand{\wks}{{\cal W}^s_k}
\newcommand{\wms}{{\cal W}^s_m}
\newcommand{\zk}{{\cal Z}_k}
\newcommand{\zks}{{\cal Z}^s_k}
\newcommand{\z}{{\cal Z}}
\newcommand{\n}{{\cal N}}
\newcommand{\ra}{{\cal R}}
\newcommand{\q}{{\cal Q}}
\newcommand{\s}{{\cal S}_n}
\newcommand{\m}{{\cal M}_{n}}
\newcommand{\req}[1]{(\ref{#1})}
\newcommand{\adj}{{\rm adj\,}}
\newcommand{\relint}{{\rm relint\,}}
\pagestyle{plain}
\bibliographystyle{plain}
\begin{document}
%======================= 
%======================= 
\begin{slide}{}
\begin{center}
{\bf  SEMIDEFINITE PROGRAMMING\\ and\\ APPLICATIONS}
\end{center}
 
~~\\
~~\\
~~\\
~~\\
~~\\
~~\\
Henry Wolkowicz \\
University of Waterloo

\end{slide}

\begin{slide}{}

\begin{center}
Central Theme:
\end{center}
\begin{description}
\item
SDP relaxations arise through the dual of the Lagrangian dual;
provides recipe for relaxations,
geometry explanations, and approximations of optimum;
\end{description}
\end{slide}
\begin{slide}{}
\begin{center}
OUTLINE
\end{center}
\begin{description}
\item[$\bullet$]
Introduction:\\
What is semidefinite programming, SDP, and why use it.
(quick)
\item[$\bullet$]
SDP relaxation:\\
 from dual of Lagrangian dual
\item[$\bullet$]
several examples, e.g. QAP, GP,\\
\item[$\bullet$]
specific example: Euclidean distance matrix problems
\end{description}
\end{slide}
\begin{slide}{}
\begin{large}
\begin{center}
INTRODUCTION
\end{center}
\end{large}
~~\
~~\
~~\
~~\
~~\
\begin{center}
Semidefinite Programming\\
 looks just like\\
Linear Programming
\end{center}
\[ {\bf (PSDP)}
\begin{array}{cccc}
    p^*=  & \max &\tr CX & (\left< C,X \right>) \\
 &  \mbox{s.t.} & {\cal A}X = b & \mbox{(linear)}\\
  && X \succeq 0,~~(X \in \p)& \mbox{(nonneg)}
    \end{array}
\]
\end{slide}
\begin{slide}{}
$\preceq$ denotes the L{\"{o}}wner partial order\\
$A\preceq B$ if $B-A \succeq 0$\\

${\cal S}_n$ denotes  $n \times n$ symmetric matrices

\[ {\cal A} :{\cal S}_n \rightarrow \Re^m
\]
\[({\cal A}X)_i =\tr (A_iX),~
  \mbox{for given}~ A_i \in {\cal S}_n 
\]

$\p$ - 
cone of positive semidefinite matrices
\begin{center}
replaces
\end{center}
$\Re^n_+$ - nonnegative orthant

\end{slide}
\begin{slide}{}

{\bf Why use SEMIDEFINITE PROGRAMMING?}

Quadratic approximations are better than linear approximations. 

{\bf Quadratic approximations are too hard to solve!}

But, we
can solve relaxations of
quadratic approximations efficiently using semidefinite programming.
\end{slide}
\begin{slide}{}
\begin{center}
APPLICATIONS
\end{center}
Finding bounds and good feasible solutions
for \underline{\bf hard combinatorial problems} such as:
max-cut; graph partitioning; quadratic assignment problem;
max-clique.

Unconstrained and constrained \underline{\bf NLP}
 techniques, e.g.:
quasi-Newton updates that preserve positive definiteness;
Trust region algorithms for large scale minimization;
Extended SQP techniques for constrained minimization.

Partial Hermitian \underline{\bf matrix completion problems} and
Euclidean distance matrix completion problems.

\underline{Engineering problems} such as:
Ricatti equations; 
min-max eigenvalue problems; matrix norm minimization; eigenvalue
localization.

\end{slide}
\begin{slide}{}
\begin{center}
{\bf SIMILARITIES WITH LP}\\
~\\
~\\
~\\
~\\
~\\
DUALITY
\end{center}
payoff function, player $Y$ to player $X$ (Lagrangian)
\[ L(X,y) :=  \tr (CX) +y^t(b-{\cal A}X)
\]

Optimal (worst case) strategy for player $X$:
\[p^* = 
      \max_{X \succeq 0 }  \min_{y }  L(X,y) 
\]
Using the {\em hidden constraint} $b-{\cal A}X=0$, recovers primal problem.

\end{slide}
\begin{slide}{}
\[ 
\begin{array}{rcl}
L(X,y) &=&  \tr (CX) +y^t(b-{\cal A}X)\\
       &=&  b^ty + \tr \left(C -{\cal A}^*y \right) X
\end{array}
\]

~\\
~\\
adjoint operator,~~ ${\cal A}^*y= \sum_i y_i A_i $
\[
\left<{\cal A}^*y,X \right>= \left<y,{\cal A}X \right>, ~~~ \forall X,y
\]


\end{slide}
\begin{slide}{}

\[p^* = 
       \max_{X \succeq 0 } \min_{y }  L(X,y) 
\leq d^*:=\min_y \max_{X \succeq 0} L(X,y) 
\]
The dual is obtained from the optimal strategy of the competing player,
Y. \\
The {\em hidden constraint} $C-{\cal A}^*y \preceq 0$\\
 yields the dual
\[ {\bf (DSDP)}
\begin{array}{ccc}
    d^*=& \min &b^ty \\
 &  \mbox{s.t.} & {\cal A}^*y \succeq  C \\
    \end{array}
\]

for the primal
\[ {\bf (PSDP)}
\begin{array}{ccc}
    p^*=  & \max &\tr CX \\
 &  \mbox{s.t.} & {\cal A}X = b\\
  && X \succeq 0
    \end{array}
\]
\end{slide}
\begin{slide}{}
Characterization of optimality for the\\
   dual pair $X,y$~~(slack $Z\succeq 0$)
 \[ 
\begin{array}{cc} 
    {\cal A}^*y -Z = C  & \mbox{dual feasibility}\\
~\\
    AX = b & \mbox{primal feasibility}\\
~\\
    ZX  = 0 & \mbox{complementary slackness}
\end{array}
\]
\[
    ZX = \mu I ~~~~~ \mbox{perturbed}
\]

Forms the basis for:\\ ~~\\
primal simplex method\\
dual simplex method\\
interior point methods
\end{slide}
\begin{slide}{}
\begin{center}
Direct SDP Relaxation
\end{center}
\begin{center}
Max-Cut Problem
\end{center}
\[
 \begin{array}{c}
    \max ~ \frac 12 \sum_{i<j} w_{ij}(1-x_ix_j),~~~x \in  \{ \pm 1 \}^n.
\end{array}
\]
Equate $x_i=1$ with $i \in \cal I$ and -1 otherwise.  \\
Let
\[ q(x) := x^tQx, \]
where $Q$ is an $n \times n$ symmetric matrix.
An equivalent problem is  the homogeneous
{\em $(\pm 1)$-quadratic programming problem}
\[
\mu^*:=\max ~ q(x),~~~x \in \{ \pm 1 \}^n.
\]
Replace $x \in \{ \pm 1 \}^n$ constraints $x_i^2=1.$

Note that for
\[ X=xx^t,  \]
\[ X \succeq 0,~ \diag (X) = e, ~q(x)=\tr XQ.  \] 
Relax the rank-1 condition on $X$ to get SDP.
\end{slide}
\begin{slide}{}
How does SDP arise from general quadratic approximations?

Let 
\[q_i(y)=\frac 12 y^tQ_iy+y^tb_i + c_i,~y\in \Re^n \]
\[ {\bf (QQP)}
\begin{array}{ccc}
    q^*=  & \min &q_0(y) \\
 &  \mbox{s.t.} & q_i(y)\leq 0\\
     &&  i=1,\ldots m
    \end{array}
\]
Lagrangian:\\
\[ 
\begin{array}{ccc}
  L(y,x) &=& q_0(y) + \sum_{i=1}^m x_iq_i(y)
\end{array}
\]
or equivalently
\[ 
\begin{array}{ccc}
  L(y,x) &=& \frac 12 y^t (Q_0 +\sum_{i=1}^m x_iQ_i)y  \\
 &&   +y^t(b_0 +\sum_{i=1}^m x_ib_i)\\
   &&+ (c_0 +\sum_{i=1}^m x_ic_i)
\end{array}
\]
\end{slide}
\begin{slide}{}

Weak duality follows from ``hidden constraint''
\[d^*=\max_{x\geq 0} \min_y L(y,x) \leq q^* = \min_y \max_{x\geq 0} L(y,x).
\]

Now homogenize; multiply linear term by new variable $y_0$
\[ y_0y^t(b_0 +\sum_{i=1}^m x_ib_i), ~~ y_0^2=1.  \]

\end{slide}
\begin{slide}{}
\[
\begin{array}{cccc}
  d^*=\\
  =\max\limits_{x \geq 0} \min_y &L(y,x)\\
  = \max\limits_{x \geq 0} 
         \min\limits_{y_0^2=1}& \frac 12 y^t (Q_0 +\sum_{i=1}^m x_iQ_i)y 
                            ~~~~~+ty_0^2 \\
 &   +y_0y^t(b_0 +\sum_{i=1}^m x_ib_i)\\
 &+ (c_0 +\sum_{i=1}^m x_ic_i)
                            ~~~~~-t\\
  = \max\limits_{x \geq 0,t} 
         \min\limits_{y}& \frac 12 y^t (Q_0 +\sum_{i=1}^m x_iQ_i)y 
                            ~~~~~+ty_0^2 \\
 &   +y_0y^t(b_0 +\sum_{i=1}^m x_ib_i)\\
 &+ (c_0 +\sum_{i=1}^m x_ic_i)
                            ~~~~~-t\\
\end{array}
\]

The hidden semidefinite constraint yields the semidefinite program. 
\end{slide}
\begin{slide}{}
We get
\[ 
B:=\left( \begin{array}{cc}
      0 & b_0^t \\ b_0 &Q_0
   \end{array}  \right),
\] 
 $A: \Re^{m+1} \rightarrow {\cal S}_{n+1}$
\[ 
A \left( \begin{array}{c}
      t \\ x
   \end{array}  \right)
   := -\left[ \begin{array}{cc}
        t &  \sum_{i=1}^m x_ib_i^t \\
      \sum_{i=1}^m x_ib_i  & \sum_{i=1}^m x_i Q_i
        \end{array}   \right]
\]
and the SDP constraint
\[
B-A \left( \begin{array}{c}
      t \\ x
   \end{array}  \right)
 \succeq 0.
\]

{\bf NOTE There is NO hidden constraint on the $Q_i$ 
if all $q_i$ are convex.} There is no need for SDP in the convex case.
Better algorithms exist, e.g. proximal methods, using quadratic cones,
...

\end{slide}
\begin{slide}{}
The dual program is equivalent to the SDP (with $c_0=0$)
\[ {\bf (D)}
\begin{array}{ccc}
    d^*=  & \mbox{sup} & -t +\sum_{i=1}^m x_ic_i \\
 &  \mbox{s.t.} & A\left( \begin{array}{c}
      t \\ x
   \end{array}  \right)
 \preceq B\\
  && x \in \Re^m, t \in \Re
    \end{array}
\]
As in linear programming, the dual of the dual is obtained from the optimal
strategy of the competing player:
\[ {\bf (DD)}
\begin{array}{ccc}
    d^*=& \inf &\tr BU \\
 &  \mbox{s.t.} & A^*U = \left( \begin{array}{c}
      -1 \\ c
   \end{array}  \right) \\
  && U \succeq 0.
    \end{array}
\]
\end{slide}
\begin{slide}{}
In some sense, the Lagrangian relaxation is the {\bf best tractable
relaxation}. 


(Is there a better one?)

\end{slide}
\begin{slide}{}

geometrical description of SDP relaxations: approximate feasible set $F$ 
by valid (not necessarily convex) inequalities:
\[ 
\begin{array}{ccc}
   F \subset \left\{ y: \sum_{i=1}^m x_iq_i(y) \leq 0 \right \}
\end{array}
\]
where $x \geq 0$ and
\[
   Q_0+\sum_{i=1}^m x_iQ_i \succeq 0.
\]


---------------------------------

Recover optimal approximation from optimal dual $x*$
\[
  \bar{y} \in \arg \min L(y,x^*)
\]
Obtained from column 1 of optimal SDP matrix $Y^*.$
Then shift to guarantee complementary slackness with valid optimal
inequality.
\end{slide}
\begin{slide}{}
{\bf The Trust Region Subproblem:}

Let
\[ q(x):= x^tAx - 2a^tx,
\]
\begin{eqnarray*}
(TRS)~~~~~ \mu^* := &\min& q(x)\\
&\mbox{s.t.}&
x^tx = s^2~~(\leq s^2).
\end{eqnarray*}

where

$A=A^t$, not necessarily semidefinite\\
$a \in \Re^n$, $s>0.$
~~\\
(application: quadratic model for unconstrained minimization)
\end{slide}
\begin{slide}{}
{\bf Homogenization of TRS}
\begin{eqnarray*}
\mu^* &=& \min\limits_{||x||=s,~y_0^2=1}  x^tAx - 2y_0a^tx \\
&=& \max\limits_t \min\limits_{||x||=s}  x^tAx - 2y_0a^tx +ty_0^2-t \\
&=& \max\limits_t \min\limits_{||x||=s,~y_0^2=1}  x^tAx - 2y_0a^tx
+ty_0^2-t \\
&=& \max\limits_t \min\limits_{||x||^2+y_0^2=s^2+1}  x^tAx - 2y_0a^tx
+ty_0^2-t
\end{eqnarray*}

\[ =  \max\limits_t (s^2+1)\lambda_1(D(t)) -t\]

$$D(t) = \left(
\begin{array}{cc}
t & -a^t \\
-a & A
\end{array}
\right).
$$
\end{slide}
\begin{slide}{}
%%This is slide 14
{\bf unconstrained dual problem to TRS}

$$D(t) = \left(
\begin{array}{cc}
t & -a^t \\
-a & A
\end{array}
\right)
$$

$y=\left( \begin{array}{c} y_0\\ x \end{array} \right)$ 
normalized eigenvector for $\lambda_{\min} D(t)$
\[
k(t) =  (s^2+1)\lambda_{\min}(D(t)) -t,
\]

\[  \mbox{***   }~~~ \max_t k(t)
\]

Note
\[  k^\prime(t) = (s^2+1)y_0^2 -1=0  \]
is feasibility for $x$
\end{slide}
\begin{slide}{}
%%This is slide 15
{\bf SDP Primal-Dual Pair}
\[
\max_t k(t) =  (s^2+1)\lambda_{\min}(D(t)) -t,
\]


add the variable $\lambda$
\[
\begin{array}{cc}
\max & (s^2+1)\lambda - t \\
\mbox{s.t.} & D(t) \succeq \lambda I
\end{array}
(DSDP)
\]

Lagrangian dual of this dual is:
\[
\begin{array}{cc}
\min & \tr D(0)X \\
\mbox{s.t.} & \tr X = s^2+1 \\
      &  X_{11} = 1\\
      & X \succeq 0
\end{array}(PSDP)
\]
\end{slide}
\begin{slide}{}
%%This is slide 16
primal-dual interior point method:

approx. solve perturbed optimality conditions
using Newton's method:
\[
\begin{array}{c}
  \tr X = s^2+1 \\
        X_{11} = 1\\
  D(t) - \lambda I - Z = 0\\
  \mu Z^{-1} - X = 0  \\
       X \succ 0, Z \succ 0
\end{array}
\]


\end{slide}
\begin{slide}{}
\begin{center}
{\bf QUADRATIC ASSIGNMENT PROBLEM\\ QAP} 
\end{center}

\[
\begin{array}{ccc}
\mu^*:= &\min\limits_{X \in \Pi} & \tr AXBX^t - 2CX^t \\
\end{array}
\]
~\\
~\\
$A, B$ and $C$ are real $n\times n$ matrices\\
$\Pi$ is the set of permutaion matrices.

~~\\
~~\\
~~\\
Rewrite as
\[
(QAP_E)~~
\begin{array}{ccl}
\mu^*:=
&\min & \tr AXBX^t - 2CX^t \\
&\mbox{~s.t.~} & XX^t = I, \left( X^tX = I\right) \\
         && \left(Xe = X^te = e \right)\\
        && X_{ij}^2 - X_{ij} = 0,~~ \forall i,j.
\end{array}
\]

ignore $Xe = X^te = e$ for now
\end{slide}
\begin{slide}{}
Find the semidefinite relaxation by taking the dual of the Lagrangian
dual.

-----------------------------------------------

We first add the (0,1)-constraints to the objective function using
Lagrange multipliers $W_{ij}$
\[
\begin{array}{ccc}
\mu_{\cal O} &=& \min\limits_{XX^t=X^tX = I} \max\limits_W
\tr AXBX^t - 2CX^t \\
       &&+ \sum_{ij} W_{ij}(X_{ij}^2 - X_{ij}).
\end{array}
\]
~~\\
~~\\
We now homogenize the objective function by multiplying by a constrained
scalar $x_0$
\[
\begin{array}{cc}
\mu_{\cal O} \geq \mu_R = \\
\max\limits_W  \min\limits_{\stackrel{XX^t=X^tX=I}{x_0^2 =1}}&
\tr \left[ AXBX^t + \right.\\ 
  &  \left.W(X \circ X)^t
    -x_0(2C+ W)X^t \right].
\end{array}
\]

\end{slide}
\begin{slide}{}
Introducing a Lagrange multiplier $w_0$ for the constraint on $x_0$ and
Lagrange multipliers $S_b$ for $XX^t=I$ and $S_o$ for $X^tX=I$ we
get
\[
\begin{array}{ll}
\mu_{\cal O} \geq \mu_R := \\
   \max\limits_W  \min\limits_{X,~x_0} &
\tr \left[ AXBX^t +  W(X \circ X)^t + w_0 x_0^2 \right. \\
         & \left. + S_b XX^t + S_o X^tX
\right] \\ 
& - \tr x_0(2C+ W)X^t\\
 &  - w_0 - \tr S_b - \tr S_o.
\end{array}
\]

We have grouped the quadratic, linear, and constant terms together.
We now define $x:= \kvec X$, $y^t:=(x_0,x^t)$ and $w^t:= (w_0,\kvec W^t)$ and 
get
\[
\begin{array}{ll}
\mu_R = \\
\max\limits_W \min\limits_{y} & y^t \left[ L_Q+Arrow(w)+\BoDiag(S_b)+
                 \right. \\
& \left. \OoDiag(S_o)  \right] y \\
& - w_0 - \tr S_b - \tr S_o
\end{array}
\]
\end{slide}
\begin{slide}{}
We used the $(n^2+1) \times (n^2+1)$ matrix
\[
L_Q := \left[ \begin{array}{cc}
0 & - \kvec (C)^t \\
-\kvec (C) & B \otimes A
\end{array} \right],
\]
and the (interesting) linear operators
\[
\Arrow (w) := \left[ \begin{array}{cc}
w_0 & - \frac 12 w_{1:n^2}^t \\
-\frac 12 w_{1:n^2} & \Diag \left(w_{1:n^2}\right)
\end{array} \right],
\]
\[
\BoDiag (S) := \left[
\begin{array}{cc}
0 & 0 \\
0 & I \otimes S_b
\end{array}
\right]
\]
and
\[
\OoDiag (S) := \left[
\begin{array}{cc}
0 & 0 \\
0 & S_o \otimes I
\end{array}
\right].
\]

\end{slide}
\begin{slide}{}
The hidden semidefinite constraint yields the equivalend SDP:
\[
(D_{\cal O})~~
\begin{array}{llc}
\max & - w_0 - \tr S_b - \tr S_o \\
\mbox{~s.t.~}& L_Q +\Arrow(w) + \\
        & \BoDiag(S_b) + \OoDiag(S_o) \succeq 0,
    \end{array}
\]
The dual of this dual yields the semidefinite relaxation.

$Y \succeq 0$ is $(n^2+1) \times ( n^2+1)$\\
     the dual matrix variable

\[
(SDP_{\cal O})~~
\begin{array}{cllcll}
\min &\tr L_QY \\
\mbox{~s.t.~}
  & \bodiag(Y) = I && \oodiag(Y) = I \\
  & \arrow(Y) = e_{0} && Y \succeq 0
\end{array}
\]
\end{slide}
\begin{slide}{}
adjoint operators are:
\[\arrow (Y) := \diag (Y) - (0, (Y_{0,1:n^2})^t. \]

$$ \bodiag(Y) := \sum\limits_{k=1}^n Y_{(k-1)n+1:kn,(k-1)n+1:kn} $$

$$ [\oodiag(Y)]_{ij} := \tr Y_{(i-1)n+1:in,(j-1)n+1:jn}  $$

\end{slide}
\begin{slide}{}
\begin{center}
{\bf Direct Approach to SDP Relaxation}
\end{center}

Let\\
 $X \in \Pi_n$ be a permutation matrix\\
$x=\kvec(X),~ c=\kvec(C).$
\begin{eqnarray*}
  q(X) &=& \tr AXBX^t - 2CX^t \\
       &=& x^t (B \otimes A) x -2c^tx \\
       &=& \trace xx^t (B \otimes A)  -2c^tx \\
       &=& \trace L_Q Y_X,
\end{eqnarray*}
where $L_Q$ is as above and
\[
Y_X := \left[
\begin{array}{cc}
1& x^t \\
x & xx^t
\end{array}
\right].
\]
\end{slide}
\begin{slide}{}

{\def\baselinestretch{0.8}\large\small
\begin{verbatim}
function [phi, X, y] = psd_ip( L);
% solves: max trace(LX) s.t. X psd, diag(X) = b;  b = ones(n,1)/4
%         min b'y       s.t. Diag(y) - L psd, y unconstrained,
% input:  L ... symmetric matrix
% output: phi ... optimal value of primal, phi =trace(LX)
%         X   ... optimal primal matrix
%         y   ... optimal dual vector
% call:   [phi, X, y] = psd_ip( L);

digits = 6;                     % 6 significant digits of phi
[n, n1] = size( L);             % problem size
b = ones( n,1 ) / 4;            % any b>0 works just as well
X = diag( b);                   % initial primal matrix is pos. def.
y = sum( abs( L))' * 1.1;       % initial y is chosen so that
Z = diag( y) - L;               % initial dual slack Z is pos. def.
phi = b'*y;                     % initial dual
psi = L(:)' * X( :);            % and primal costs
mu = Z( :)' * X( :)/( 2*n);     % initial complementarity
iter=0;                         % iteration count

disp(['      iter    alphap    alphad      gap     lower      upper']);

while phi-psi > max([1,abs(phi)]) * 10^(-digits)

      iter = iter + 1;          % start a new iteration
      Zi = inv( Z);             % inv(Z) is needed explicitly
      Zi = (Zi + Zi')/2;
      dy =  (Zi.*X) \ (mu * diag(Zi) - b);      % solve for dy
      dX = - Zi * diag( dy) * X + mu * Zi - X;  % back substitute for dX
      dX = ( dX + dX')/2;       % symmetrize

\end{verbatim}
}

\end{slide}
\begin{slide}{}


{\def\baselinestretch{0.8}\large\small
\begin{verbatim}
% line search on primal
      alphap = 1;               % initial steplength
      [dummy,posdef] = chol( X + alphap * dX );  % test if pos.def
      while posdef > 0,
              alphap = alphap * .8;
              [dummy,posdef] = chol( X + alphap * dX );
              end;
      if alphap < 1, alphap = alphap * .95; end; % stay away from boundary
% line search on dual; dZ is handled implicitly: dZ = diag( dy);
      alphad = 1;
      [dummy,posdef] = chol( Z + alphad * diag(dy) );
      while posdef > 0;
              alphad = alphad * .8;
              [dummy,posdef] = chol( Z + alphad * diag(dy) );
              end;
      if alphad < 1, alphad = alphad * .95; end;
% update
      X = X + alphap * dX;
      y = y + alphad * dy;
      Z = Z + alphad * diag(dy);
      mu = X( :)' * Z( :) / (2*n);
      if alphap + alphad > 1.8, mu = mu/2; end; % speed up for long steps
      phi = b' * y; psi = L( :)' * X( :);
% display current iteration
        disp([ iter alphap alphad (phi-psi) psi  phi ]);

        end;            % end of main loop
\end{verbatim}
}
\end{slide}
\begin{slide}{}
 486 66 Mhz, PC using a C encoded version;\\
\# of iterations does not seem to depend significantly on $n$\\
$hh$ gives the number of hours, $mm$ the number of minutes, 
$ss$ the seconds.\\
random unweighted graphs with edge probability $p=0.5$.\\
stopping condition to 6 digits of accuracy\\
\end{slide}
\begin{slide}{}
\begin{center}
\begin{tabular}{|r|r|r|}
\hline
\multicolumn{1}{|c}{size} $n$ &
\multicolumn{1}{|c}{iterations}&
\multicolumn{1}{|c|}{hh:mm:ss}\\
\hline
100 & 14 &      48\\
150 & 12 &    2:30\\
200 & 12 &    6:30\\
250 & 13 &   14:24\\
300 & 14 &   30:24\\
400 & 14 & 1:24:36\\
500 & 14 & 2:13.18\\ \hline
\end{tabular}
\end{center}
\end{slide}
\begin{slide}{}
included triangle constraints;\\
favorable to use Mehrotra's predic\-tor -- corrector approach\\
system is solved twice with the same 
factorization but with different right hand sides.\\
The triangle
inequalities were included successively according to their 
amount of violation\\
programs have $n$ equality and 4${n\choose 3}$
inequality constraints in the dual program.\\
 The final solution satisfies
all these constraints.\\
\# of iterations astonishingly small\\
e.g.  problems with $n= 100$
satisfies $4{100\choose 3}= 646,800$ inequality constraints.\\
\end{slide}
\begin{slide}{}
\begin{center}
\begin{tabular}{|r|r|r@{ (}c@{) }|r|r|}\hline
\multicolumn{6}{|l|}{random graphs, edge-weights $\in\{-10,\ldots,10\}$}\\ \hline
\multicolumn{1}{|c}{size} $n$ &
\multicolumn{1}{|c}{cut}&
\multicolumn{1}{|c@{ (}}{upper bnd}&
\%&
\multicolumn{1}{|c}{hh:mm:ss}&
\multicolumn{1}{|c|}{iters}\\ \hline\hline
 30 &  330 &  330.0001 &0.0&       54 &  59 \\
 40 &  619 &  619.0000 &0.0&     2:02 &  60 \\
 50 &  708 &  708.0001 &0.0&  1:04:22 & 166 \\
 70 & 1215 & 1248.7142 &2.8&  4:42:21 & 220 \\
100 & 2440 & 2540.2528 &4.1& 37:56:05 & 266 \\ \hline\hline
\hline
\hline
\multicolumn{6}{|l|}{ random unweighted graphs, edge probability $0.5$}\\ \hline
\multicolumn{1}{|c}{size} $n$ &
\multicolumn{1}{|c}{cut}&
\multicolumn{1}{|c@{ (}}{upper bnd}&
\%&
\multicolumn{1}{|c}{hh:mm:ss}&
\multicolumn{1}{|c|}{iters}\\ \hline\hline
 30 &  143 &  143.0000 &0.0&       56 &  62 \\
 40 &  245 &  245.6864 &0.3&    59:19 & 247 \\
 50 &  374 &  375.3689 &0.4&  1:54:18 & 250 \\
 70 &  723 &  726.0399 &0.4&  5:57:22 & 257 \\ \hline
\end{tabular}
\end{center}
\end{slide}
\begin{slide}{}
\begin{center}
{\bf Euclidean Distance Matrix Completion Problem}
\end{center}
$n \times n$ symmetric nonnegative (elementwise)  
matrix $D=(d_{ij})$ with zero 
diagonal is called a {\em pre-distance matrix}. 

if there 
exist points $x^1,x^2,\ldots,x^n$ in $\Re^r$ such that
\beq  
d_{ij} = {\| x^i- x^j\|}^2, \1 \1 \1 i,j=1,2,\ldots,n,
\eeq
then $D$ is called a {\em Euclidean distance matrix} (EDM). 

The smallest value of $r$ is called {\em the embedding dimension} of $D$.  
Note that $r$ is always $\leq n-1$. 

\end{slide}
\begin{slide}{}
a pre-distance matrix 
$D$ is a EDM if and only if $D$ is negative semidefinite on 
\[ M:=\left\{ x \in \Re^n : x^t e = 0 \right\},
\]
the orthogonal complement of $e$,
where $e$ is the vector of all ones.

Define the $n \times n$ orthogonal matrix 
\beq \label{defV} 
   Q:= \left[ \frac{1}{\sqrt{n}}e \; | \; V \right], \1 \1 Q^tQ = I. 
\eeq
Thus $V^te=0$ and $~V^tV=I$. Moreover, the subspace $M$ can be represented 
as the range of the $n \times (n-1)$ matrix $V$ and 
\beq \label{eq:Vmp}
J := V V^t= I- \frac{e e^t}{n}
\eeq
is the orthogonal projection onto $M$.

\end{slide}
\begin{slide}{}
Now define the {\em centered} and {\em hollow} subspaces
\beq \begin{array}{rcl}
\Sc &:=&  \{ B \in \Sn :  Be = 0 \}, \\ 
\Sh& := & \{ D \in \Sn :  \diag(D) = 0 \}, 
\end{array}
\eeq 
where diag($D$) denotes the column vector formed from the diagonal of $D$.

define the two linear operators acting on
$\Sn$
\beq \begin{array}{rcl} \label{KK} 
\KK(B)& := &  \mbox{diag}(B)\,e^t + e \, \mbox{diag}(B)^t - 2B,
\end{array} \eeq 
and
\beq \begin{array}{rcl} \label{T} 
\TT(D)& := &  -\frac 12 JDJ.
\end{array} \eeq
The operator $- 2 \TT$ is an orthogonal projection onto $\Sc;$ 
thus it is a self-adjoint idempotent.

\end{slide}
\begin{slide}{}
\begin{thm} \label{KT} 
The linear operators satisfy
\begin{eqnarray*}
\KK (  \Sc) = \Sh, \\
\TT (  \Sh) = \Sc, 
\end{eqnarray*} 
and $\KK_{|\Sc}$ and $\TT_{|\Sh}$ are inverses of each other.
\end{thm} 
~~ \epr

\begin{center}
a hollow matrix $D$ is EDM\\
if and only if\\
 $B=\TT(D)$ is positive
semidefinite, ( $B \succeq 0$) 
\end{center}

Equivalently, $D$ is EDM
 if and only if
$D=\KK(B),$
for some $B$ with $Be=0$ and $B \succeq 0$.   
In this case the embedding dimension $r$ is given
by the rank of $B$. 


\end{slide}
\begin{slide}{}

composite operators
\beq \begin{array}{rcl} \label{KV} 
\KK_V(X)& := &  \KK( V X V^t),
\end{array} \eeq 
and
\beq \begin{array}{rcl} \label{TV} 
\TT_V(D)& := &  V^t\TT( D)V= - \frac 12 V^t D V,
\end{array} \eeq 
where $V$ is defined in (\ref{defV}).  

\begin{lem} \label{KVTV} 
\begin{eqnarray*} 
\KK_V ( \snn) =\Sh, \\
\TT_V ( \Sh) =\snn, 
\end{eqnarray*} 
and $\KK_V$ and $\TT_V$ are inverses of each other on these two spaces.
\end{lem}
~~\epr

\end{slide}
\begin{slide}{}

\begin{cor} \label{co} 
Suppose that $V$ is defined as in (\ref{defV}).  Then:
\begin{eqnarray*} 
\KK_V(\p)  & = &  EDM , \\
\TT_V(EDM) & = &\p.
\end{eqnarray*} 
\end{cor}
~~\epr

Note that $V$ as defined in (\ref{defV}) is not unique. In our code
we use  
\beq
\label{eq:Vorth}
  V := \left[ \begin{array}{cccc}  y & y& \ldots  & y \\
                        1+x & x & \ldots, &  x \\  
                         x & 1+x & \ldots, &  x \\  
                           ... &...  & \ddots &...   \\
                         x & x & \ldots, & 1+x \\  
                      \end{array} \right],
\eeq  
where $x = \frac{-1}{n + \sqrt{n} }$ and $y= \frac{-1}{\sqrt{n}}$. With this
choice, it can be 
easily verified that $V^te=0$, $V^tV=I$, and $VV^t=J$ as required by
(\ref{defV}). 

\end{slide}
\begin{slide}{}
{\bf  Program Formulation} 

Let $A$ be a predistance matrix, i.e. symmetric, nonnegative
(elementwise) and with 0 diagonal. Let
$H$ be an $n \times n$ symmetric nonnegative (elementwise) matrices and
let $\|A\|_F= \sqrt{A^tA}$ denote the {\em Frobenius norm} of $A.$ Consider
the objective function
\[ f(D) := {\| H \circ (A - D) \|}^2_F,   \] 
where $\circ$ denotes {\em Hadamard product}. 
The weighted, {\em closest Euclidean distance
matrix problem} is 
\[ 
(CDM_0)
 \bt{ccc}
          $\mu^*$ := & $\min$   &   $f(D)$  \\
                  & $\mbox{ subject to }$ & $D \in EDM. $ 
  \et
\]

WLOG:\\
diag($H$) = 0\\
graph of $H$ is {\em connected}\\

\end{slide}
\begin{slide}{}
define
\[ f(X) := {\| H \circ (A - \KK_V ( X)) \|}^2_F
            \]
\[
(CDM)
 \bt{ccc}
          $\mu^*$ := & $\min$   &   $f(X)$  \\
                  &  subject to  &  $X \succeq 0.$
  \et
\]

We can recover the optimal distance matrix using the optimal $X$ and the 
relation 
\[ D = \KK_V(X).  \]  

\end{slide}
\begin{slide}{}
We can now state the dual problem 
\beq 
\label{eq:dualprob}
\bt{cccc} 
& & max & $f(X)+\langle y, b - \A(X) \rangle -
     $ trace $ \Lambda X$ \\  
       &           &  subject to & $\nabla f(X) - \A^*y- \Lambda = 0$  \\
       &           &             & $  \Lambda \succeq 0, (X \succeq 0). $  
\et  \eeq 

Strong duality holds (Slater's holds for the primal)

Slater's holds for the dual if and only if the graph of $H$ is
connected.
\end{slide}
\begin{slide}{}
REMARKS:
\begin{description}
\item
An initial stepsize of 1.5 was taken; tests were also done with 2 and
improved convergence was noted.
\item
A centering parameter $\sigma$ changes adaptively
depending on the length of the stepsize, i.e. it decreases if the
stepsize is large and increases to 1 if the stepsize is small.
\item
The barrier parameter $\mu$ also changes adaptively; there is an
aggressive decrease $\mu/4n$ if the stepsize is large; this decrease is
smaller when small stepsizes are chosen.
\end{description}
\end{slide}
\begin{slide}{}
\begin{description}
\item
Backtracking was initially
done with the value .95 to be conservative, as is
usually the case with SDP. However, tests were done with the aggressive
.997 as is done in LP. This did not result in faster convergence.
However, robustness was maintained.
\item
The upper part of $F^{\prime}$, i.e. $[ F^{\prime}_{u1}~ F^{\prime}_{u2}],$ 
need only be evaluated once
at the beginning; while the lower (simpler) part of the operator needs to be
evaluated at each iteration.
Both evaluations take very little time compared to the solution of the
least squares problem.
\end{description}
\end{slide}
\begin{slide}{}
\begin{description}
\item
No predictor-corrector approach was used. Rather, centering was used to
speed up convergence.
However, the step lengths were
consistently larger than 1 which suggests that no corrector step is
needed, i.e. centering keeps 
the steps very close to the central path.
\item
Though the initial starting points were feasible, dual feasibility was
lost after the first iteration, i.e. unlike Newton's method, the
Gauss-Newton direction does not maintain feasibility. However,
complementarity and dual feasibility converged together to 0, i.e. at a
similar rate. This might account for the strong convergence obtained
with this direction.
\end{description}
\end{slide}
\begin{slide}{}
\begin{description}
\item
The search direction was very robust, i.e. machine accuracy was
attainable whenever asked for. This was the case even when the
aggressive initial stepsize and reduction of $\mu$ were used.
\item
Problem dimension was restricted to $n \leq 42$. 
This was because of the high cost of the least squares solution
in MATLAB. (Done
using qls, Version 1.6, from Pontus Matstoms, Linkoping University.)
However, if one solves the
least squares problem using, for example, a Lanczos approach, then one
need never form the matrices which represent $F^{\prime}$ and thus
can exploit sparsity and solve much larger problems. These tests are
ongoing.
\end{description}
\end{slide}
\begin{slide}{}
\begin{description}
\item
We have assumed, without loss of generality, that the graph of the
matrix $H$ was connected, i.e. if this was not the case then we can
solve two smaller problems independently and more simply. In fact, one
should always solve the smaller problems independently. Note that
Slater's constraint qualification is always satisfied for the primal and
so there is no duality gap. However, if Slater's constraint
qualification fails for the dual (which is equivalent to the graph being
disconnected), then the set of optimal
solutions for the primal is unbounded and the algorithm has great
difficulty converging. Here is an example of a class of problems where
strong duality holds but interior-point methods have difficulty.
\end{description}
\end{slide}
\begin{slide}{}
\end{slide}
\end{document}