%======================= main.tex ===========================================
\documentclass{slides}
\newtheorem{exam}{Example}
\newtheorem{prop}{Proposition}
\newtheorem{lem}{Lemma}
\newtheorem{thm}{Theorem}
\newtheorem{cor}{Corollary}
\newcommand\T{{\cal T}}
\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}
\begin{document}
%======================= 
%======================= 
\begin{slide}{}
\begin{center}
{\bf  Newton Directions in Primal-Dual Interior-Point Methods for
Semidefinite Programming}
\end{center}
 

~~\\
~~\\
~~\\
~~\\
~~\\
~~\\
Henry Wolkowicz \\
University of Waterloo\\

\end{slide}
\begin{slide}{}
in collaboration with:\\
~~\\
~~\\
Serge Kruk,\\
Masakazu Muramatsu,\\
    Franz Rendl,\\
Robert J. Vanderbei

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

\begin{center}
Central Theme:
\end{center}
\begin{description}
\item
Primal-dual interior-point methods for SDP do not require
symmetrization when finding the Newton direction.

(In fact, the Newton direction can be considered as a Gauss-Newton
direction.)
\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$]
basic duality, complementarity conditions
\item[$\bullet$]
Differences with LP:\\
duality gap, strict complementarity, symmetrization
\item[$\bullet$]
Search directions:\\
different symmetrizations, Gauss-Newton, least squares
\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}{}
Other dual pairs:

\[ {\bf (P)}
\begin{array}{ccc}
    p^*=& \max &b^ty \\
 &  \mbox{s.t.} & {\cal A}^*y \preceq  C \\
    \end{array}
\]

\[ {\bf (D)}
\begin{array}{ccc}
    d^*=  & \min &\tr CX \\
 &  \mbox{s.t.} & {\cal A}X = b\\
  && X \succeq 0
    \end{array}
\]
\end{slide}
\begin{slide}{}

\begin{center}
However, there are several subtle \\
 {\bf DIFFERENCES with LP}.
\end{center}
1. 
Duality gaps can exist for SDP 
in the absence of strictly feasible solutions (Slater's constraint
qualification, CQ).

2. Strict complementarity can fail at the optimum.

$\Rightarrow$~~
3. The perturbed optimality conditions involve nondiagonal matrices and
map between different spaces\\
 (overdetermined -does it require symmetrization?).

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

1. Duality gaps:

\underline{Example}

If the primal is
\[ {\bf (P)}
\begin{array}{ccc}
    p^*=  & \max &x_2 \\
 &  \mbox{s. t.} & 
\left[ \begin{array}{ccc}
       x_2 & 0 & 0\\
       0 & x_1 & x_2\\
       0 & x_2 & 0
       \end{array}  \right]
\preceq 
   \left[ \begin{array}{ccc}    
       1 & 0 & 0\\
       0 & 0 & 0\\
       0 & 0 & 0
       \end{array}  \right] 
    \end{array}
\]
then the dual is
\[ {\bf (D)}
\begin{array}{ccc}
    d^*=& \min &\tr U_{11}\\
 &  \mbox{s. t.} & U_{22}=0  \\
 &   & U_{11}+2U_{23}=1  \\
  && U \succeq 0.
    \end{array}
\]
Then $p^*=0 < d^*=1.$
\end{slide}
\begin{slide}{}
2. Strict complementarity:

\underline{Example}

If the primal is
\[ {\bf (P)}
\begin{array}{ccc}
    p^*=  & \max &x_1 \\
 &  \mbox{s. t.} & 
\left[ \begin{array}{ccc}
       x_1 & x_3 & x_2\\
       x_3 & x_2 &  0 \\
       x_2 &  0  & x_3
       \end{array}  \right]
\preceq 
   \left[ \begin{array}{ccc}    
       0 & 0 & 0\\
       0 & 0 & 0\\
       0 & 0 & 1
       \end{array}  \right] 
    \end{array}
\]
then the dual is
\[ {\bf (D)}
\begin{array}{ccc}
    d^*=& \min &\tr U_{33}\\
 &  \mbox{s. t.} & U_{11}=1  \\
 &   & U_{22}+2U_{13}=0  \\
 &   & U_{33}+2U_{12}=0  \\
  && U \succeq 0.
    \end{array}
\]
Then the (unique) optimum pair is:\\
 $U=\Diag(1,0,0)$\\
 and $x=[0,0,0]$ with slack $Z=\Diag(0,0,1)$
\end{slide}
\begin{slide}{}
\begin{center}
INTERIOR-POINT METHODS\\ for\\ LP and SDP
\end{center}
started with Karmarkar 1984
~~\\
~~\\

primal-dual interior-point methods:
\begin{enumerate}
\item
potential reduction
\item
 path following
\item
 infeasible-interior-point
\item
predictor-corrector
\end{enumerate}

(see e.g. recent book by Steve Wright, or interior-point archives)

\end{slide}
\begin{slide}{}
Connection between Karmarkar's algorithm and log-barrier problems
(Stanford gang of 4: Gill, Murray, M. Wright, Saunders)

\[ {\bf (D)}
\begin{array}{ccc}
    d^*=& \min &b^ty \\
 &  \mbox{s.t.} & {\cal A}^*y -Z=  C \\
  & & Z \succeq 0
    \end{array}
\]

log-barrier problem for dual SDP, with fixed $\mu > 0$ 
\[  
\begin{array}{ccc}
&\min& b^ty - \mu \log \det Z\\
 &  \mbox{s.t.} & {\cal A}^*y -Z=  C \\
    \end{array}
\]
and $\mu \downarrow 0$

Lagrangian for log-barrier
\[  
\min b^ty  +\tr X\left(C -  {\cal A}^*y +Z \right) - \mu \log \det Z 
\]
\end{slide}
\begin{slide}{}
Optimality conditions for the\\
   dual log-barrier\\
solve $F_{\mu} (X,y,Z) =0$, with $X,Z \succ 0$
 \[ 
\begin{array}{ccc} 
\frac {\partial}{\partial X}:&  C-  {\cal A}^*y +Z  =0  & \mbox{dual feasibility}\\
~\\
   \frac {\partial}{\partial y}:&  b- {\cal A}X = 0 & \mbox{primal feasibility}\\
~\\
   \frac {\partial}{\partial Z}:& X-\mu Z^{-1}  = 0 & \mbox{pert. compl. 
          slack.}
\end{array}
\]

or, to avoid ill-conditioning,  use:
\[
    ZX-\mu I=0~~ \mbox{perturbed complementary slackness}
\]

Note:\\
the nonlinear system $F_{\mu} (X,y,Z) =0$ is 
overdetermined if (more linear) $ZX-\mu I=0$ is used.\\
(this is {\bf not} the case for LP)

\end{slide}
\begin{slide}{}
\begin{center}
CENTRAL PATH
\end{center}
For each $\mu > 0$, there exists a unique
$(X_{\mu},y_{\mu},Z_{\mu}),$ with $X_{\mu},Z_{\mu} \succ 0,$
that solves 
the nonlinear system $F_{\mu} (X,y,Z) =0.$ 

\[
{\cal C} =
\{ (X_{\mu},y_{\mu},Z_{\mu}): \mu >0  \}
\]

well defined,\\
 converges to optimum as $\mu \downarrow 0,$\\
long steps possible in neighbourhoods near ${\cal C}$

p-d i-p methods followed:\\
Megiddo 1987, Kojima, Mizuono, and Yoshise 1987, many others ....


\end{slide}
\begin{slide}{}
In numerical analysis, how do we
solve a nonlinear system of equations?

The standard approach is to solve the equivalent
 nonlinear least squares problem using
``Gauss-Newton'' method. 
\[  \min \frac 12 || F_{\mu} (X,y,Z) ||^2_2  \]
subject to $X,Z>0.$

\end{slide}
\begin{slide}{}
In the LP case,
when we differentiate we get a square system and Gauss-Newton
reduces to Newton's method.

i.e. solving
\[
\left( F_{\mu}^{\prime}\right)^t F_{\mu}^{\prime} (\Delta s)
   = - \left( F_{\mu}^{\prime}\right)^t F_{\mu}
\]
(where $F_{\mu}^{\prime}$ is the Jacobian)\\
 is equivalent to solving
\[
 F_{\mu}^{\prime} (\Delta s)
   = -  F_{\mu}
\]



\end{slide}
\begin{slide}{}
{\bf framework} p-d i-p:\\
{\bf Given} $(X^0,y^0,Z^0) \in {\cal F}^0$\\
{\bf for} $k=0,1,2 \ldots $\\
\hspace{.5in}{\bf solve} for the search direction\\
     \[
  ~~~~F_{\mu}^{\prime}(X^k,y^k,Z^k) 
\left( 
\begin{array}{ccc}
\Delta X^k \\ \Delta y^k \\ \Delta Z^k 
\end{array}
\right)
=  
\left( 
\begin{array}{ccc}
0 \\ 0 \\ -X^kZ^k + \sigma_k \mu_k I 
\end{array}
\right)
\]
\hspace{.5in}where $\sigma_k$ centering, $\mu_k=\tr X^kZ^k/n$

\[ \begin{array}{c}
(X^{k+1},y^{k+1},Z^{k+1}) =~~~~~~~~~~~~~~~~~~\\
~~~~~~~~~~~~~~~~ (X^k,y^k,Z^k) + \alpha_k
(\Delta X^k,  \Delta y^k,  \Delta Z^k)
\end{array}
\]
\hspace{.5in}so that $(X^{k+1},Z^{k+1}) \succ 0$\\
{\bf end (for)}.

\end{slide}
\begin{slide}{}
Complications in the SDP case:

Overdetermined system, i.e.
should we symmetrize $ZX$, e.g. $ZX+XZ$, many others ...

line search

strict complementarity may fail

duality gaps

\end{slide}
\begin{slide}{}
\begin{center}
SYMMETRIZATION
\end{center}
from log-barrier we get:\\
$X-\mu Z^{-1} =0$\\
with linearization:\\
$X + \Delta X - \mu \left(Z^{-1} - Z^{-1}\Delta ZZ^{-1} \right) = 0$

more linear version  - multiply by $Z$:\\
$ZX - \mu I=0$\\
with linearization:\\
$ZX  - \mu I +Z \Delta X + \Delta ZX = 0$

symmetrize:\\
$ZX+ZX - 2\mu I=0$\\
with linearization:\\
$ZX +ZX  - 2\mu I +Z \Delta X + \Delta XZ+ \Delta ZX  +X \Delta Z  = 0$

Others:\\
$Z^{1/2} XZ^{1/2} - \mu I=0$, etc...\\
$H_P(M):= \frac 12 \left( PMP^{-1} + (PMP^{-1})^t \right)$\\
$M=XZ$ and $P$ any nonsingular matrix

\end{slide}
\begin{slide}{}
$XZ+ZX$:\\
 advantage: fast convergence and stability;\\
 disadvantage: solve Lyapunov equation when finding Newton direction

Heuristic; symmetrize $\Delta X$ after it is calculated:\\
advantages: fast convergence  and inexpensive
Newton direction, especially for max-cut problem!\\
disadvantage: stability

Still many open questions, e.g. why symmetrize?\\
 Can apply
Gauss-Newton method to solving overdetermined system of nonlinear
equations.
\end{slide}
\begin{slide}{}
``best direction'':  Newton direction for least squares problem

best compromise:  Gauss-Newton direction

other directions: are approximations of the Gauss-Newton direction
\end{slide}
\begin{slide}{}
\begin{center}
Shrinking the Overdetermined System
\end{center}


We can find redundant equations and solve a
smaller system.  We can remove
the lower triangular part from the nonlinear equation and get
the equivalent equation
\[
\T \left(ZX  - \mu I \right)= 0,
\]
where the linear operator $\T: \m \rightarrow \s$
by ignoring the strictly lower triangular part of the matrix, i.e. the
$i,j$ components are
\[
\left(\T(W)\right)_{ij} = \left\{
\begin{array}{cc}
W_{ij} & \mbox{if}~ i \leq j \\
W_{ji} & \mbox{if}~ i > j \\
\end{array}
\right.
\]


The resulting system is now square.

$\T$ is an orthogonal projection
\end{slide}
\begin{slide}{}

The optimality conditions become
\[
\left(
\begin{array}{c}
Z + C - \A^{*}y \\
b-\A(X)\\
\T(ZX-\mu I)\\
\end{array}
\right)
 = 0.
\]
This does not introduce new nonlinearities but does guarantee that we
map between the same spaces. In fact, we can prove that we do not lose
information in the optimality conditions
when we only consider the upper triangular parts.

LEMMA\\
Suppose that $Z \succ 0.$ Then
the linear operator $\T_Z(\cdot) := \T (Z  \cdot )$ is a one-one
mapping on  $\s$. In addition, if
$\T_Z(X) = \mu I$, for some $X=X^t$, then $ZX=\mu I.$

\end{slide}
\begin{slide}{}
The Newton direction is obtained from
solving the system
\begin{eqnarray*}
        \Delta Z-\A^{*}(\Delta y)        &=&  -F_{d}  \label{shnewt1}\\
        -\A(\Delta X)                     &=&  -F_{p}  \label{shnewt2}\\
      \T(Z\Delta X+\Delta Z X)             &=&  -\T(Z F_{Z,X}) \label{shnewt4}.
\end{eqnarray*}
Equivalently,
\[
 \left[
\begin{array}{ccc}
0 & -\A^* & I \\
-\A & 0 & 0 \\
\T_Z (\cdot) & 0 & 
\T(\cdot X)
\end{array}
\right]
\left( \begin{array}{c}
    \Delta X \\
    \Delta y \\
\Delta Z \\
            \end{array} \right)
=
\left( \begin{array}{c}
     -F_d\\
     -F_p\\
     -\T_Z(X)+\mu I \\
            \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}{}
Simple convex program:
\[ {\bf (CP)}
\begin{array}{ccc}
    p^*=  & \min &f(x) \\
 &  \mbox{s.t.} & x \in \Omega\\
    \end{array}
\]
where:\\
 $f: X \rightarrow \Re $ -  differentiable, convex function\\
$\Omega \subset X$ is an abstract inner-product space\\
$\Omega$ -  convex set

Pshenichnyi Condition:\\
$x^* \in \Omega$ solves (CP) {\bf iff} $\nabla f(x^*) \in \left( \Omega
- x^* \right)^+$

polar cone
\[ K^+= \{ \phi \in X^*: \left< \phi, k \right> \geq 0, ~ \forall k \in K \}
\]

special case when $\Omega$ closed convex cone\\
$x^* \in \Omega$ solves (CP)\\
 {\bf iff}\\
 $\nabla f(x^*) \in \left( \Omega
 \right)^+$ and\\
  $\left<\nabla f(x^*),x^*\right> = 0$ (compl. slack.)
\end{slide}
\begin{slide}{}
special case when $\Omega = \Re^n_+$\\
if $z,x \in  \Re^n_+$\\
$z^tx=0$ iff $z_ix_i=0$ for all $i$

special case when $\Omega = \p$, psd matrices\\
if $Z,X \in  \p$\\
$\left<Z,X\right>=\tr ZX=0$ iff $ZX=0$

However, these special cases do not hold in general. The purpose of {\em
complicating} the complementarity condition is to obtain the {\em
powerful} central path concept. (Avoids convergence to nonfeasible -but
stationary - points.)



\end{slide}
\begin{slide}{}
for some $y$ and $X$ feasible for ${\cal A}X=b$
\[
\begin{array}{ccc}
 p^*=\max  \tr (CX) +y^t(b-{\cal A}X)
 &  \mbox{s.t.} & X \succeq 0\\
    \end{array}
\]

Optimality Conditions,\\
 $F(X,y,Z)=0$, $X,Z \succeq 0$:
 \[ 
   \nabla F: \left\{
\begin{array}{cc} 
C-  {\cal A}^*y  - Z=0, & \mbox{dual feasibility}\\
~\\
    b- Ax   =0& \mbox{primal feasibility}\\
~\\
    \left<Z,X\right>  = 0 & \mbox{complementary slackness}
\end{array}   \right.
\]

Recall:\\
$\left<Z,X\right>=\tr ZX=0$ iff $ZX=0$

The above system is {\bf underdetermined}\\
Why not use simpler underdetermined system?

\end{slide}
\begin{slide}{}
\begin{center}
STABILITY QUESTIONS
\end{center}
Let $f: \Re^n \rightarrow \Re$ and
$g: \Re^n \rightarrow \Re$ be  twice continuously differentiable
functions\\

  (e.g. $f$ is linear part of log-barrier
problem, i.e. it is the Lagrangian of
(PSDP); while $g(x)= \det (x),$ is the determinant.)\\

consider
\[
(NLP)~~~~~~~\begin{array}{lll}
        &  \min  & f(x) \\
&\mbox{subject to} & g(x) \geq  0,  \\
            \end{array}
\]

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

Optimality under CQ:
\[
\begin{array}{rcl}
\nabla f(x) - \lambda \nabla g(x) &=& 0\\
\left< \lambda , g(x) \right> &=& 0 \\
g(x) &\geq & 0 \\
\lambda &\geq & 0.
\end{array}
\]
$\mu > 0$ fixed log-barrier parameter\\
log-barrier problem is:
\[
\min_x B_{\mu}(x):= f(x) - \mu \log g(x).
\]
stationarity
\[
0=\nabla B_{\mu}(x_{\mu})= \nabla f(x_{\mu}) - 
                \frac {\mu}{g(x_{\mu})} \nabla g(x_{\mu}),
  ~ g(x_{\mu}) > 0.
\]
we can assume $x_{\mu} \rightarrow x^*$ and
$g(x_{\mu}) \rightarrow 0$\\
Therefore
\[
\frac {\mu}{g(x_{\mu})} \rightarrow \lambda^*
\]
\end{slide}
\begin{slide}{}

Hessian
\[
\nabla^2 P_{\mu}(x)= 
\]
\[\left\{ \nabla^2 f(x) - \frac {\mu}{g(x)} \nabla^2 g(x)
                      \right\}
           +\left[ \left( \frac {\mu}{g(x)}\right)
         \frac 1{g(x)}  \nabla g(x)\nabla g(x)^t \right].
\]
first \{ \} part converges\\
second [ ] part is rank one $\Rightarrow$ ill-conditioning\\

Fix: add new variable $\Delta v \in \Re$
\[\Delta v = -\frac {\mu}{g^2(x)} \nabla g(x)^t \Delta x.\]
Newton direction: $d=(\Delta x, \Delta v)^t,$
\[
\left[ \begin{array}{ll}
   \nabla^2 f(x_{\mu}) - \frac {\mu}{g(x_{\mu})} \nabla^2 g(x_{\mu}) &  
          \nabla g(x_{\mu}) \\
\nabla g(x_{\mu})^t & \frac {g(x_{\mu})^2}{\mu}   \\
            \end{array} \right]
\left( \begin{array}{c}
    \Delta x \\
\Delta v \\
            \end{array} \right)
\]
\[
=
\left( \begin{array}{c}
    - \nabla f(x_{\mu}) +
                \frac {\mu}{g(x_{\mu})} \nabla g(x_{\mu})\\
0 \\
            \end{array} \right).
\]
\end{slide}
\begin{slide}{}

As $\mu \rightarrow 0$ the Hessian becomes
\[
\left[ \begin{array}{cc}
   \nabla^2 f(x^*) - \lambda^* \nabla^2 g(x^*) &  \nabla g(x^*) \\
\nabla g(x^*) & 0  \\
            \end{array} \right].
\]

This matrix is nonsingular under the standard second order sufficient
optimality conditions. 
i.e. that the Hessian of the Lagrangian is
positive definite on the tangent space; equivalently,
$ \nabla^2 f(x^*) - \lambda^* \nabla^2 g(x^*)$ is positive definite on the
orthogonal complement of $\nabla g(x^*).$


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

alternative to overcome ill-conditioning

use perturbed KKT conditions\\
use substitution $\lambda =
\frac {\mu}{g(x)}$ or $\lambda g(x) = \mu$ 

Then becomes
\[
\begin{array}{rcl}
\nabla f(x) - \lambda \nabla g(x) &=& 0\\
\lambda g(x) &= & \mu \\
g(x) &> & 0 \\
\end{array}
\]

i.e. $\lambda g(x) = \mu$ avoids ill-conditioning that arises when
differentiating $\lambda = \frac {\mu}{g(x)}$.

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

Summary:\\
two ways to avoid ill-conditioning of
log-barrier problem. 

1. increase the dimension of the
problem by one -  still finds the Newton direction 

2.
change the optimality conditions (perturb)\\
direction based solely on the optimality conditions and not on a
log-barrier problem\\
(not a Newton direction)

\end{slide}
\begin{slide}{}
Alternate motivation for changing true
log-barrier direction:\\
$X-\mu Z^{-1} =0$\\
with linearization:\\
$X + \Delta X - \mu \left(Z^{-1} - Z^{-1}\Delta ZZ^{-1} \right) = 0$

is to get more linear version  - multiply by $Z$:\\
$ZX - \mu I=0$\\
with linearization:\\
$ZX  - \mu I +Z \Delta X + \Delta ZX = 0$

\begin{center}
MORE LINEARITY\\
 SPEEDS UP
 NEWTON'S METHOD\\
    but\\
results in overdetermined system
\end{center}

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

