%======================= main.tex ===========================================
\documentclass{slides}
\newtheorem{exam}{Example}
\newtheorem{prop}{Proposition}
\newtheorem{lem}{Lemma}
\newtheorem{thm}{Theorem}
\newtheorem{cor}{Corollary}
\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}
%\blackandwhite{example}
%======================= example.tex
%  the file example.tex must exist - it is empty though
%======================= 
%======================= 
\begin{slide}{}
\begin{center}
{\bf Newton Directions and Symmetrization\\
in\\
 Primal-Dual  Interior-Point Methods\\
 for\\
Semidefinite Programming}
\end{center}
 

~~\\
~~\\
~~\\
~~\\
~~\\
~~\\
Henry Wolkowicz \\
University of Waterloo\\
(with Serge Kruk and Franz Rendl)

\end{slide}

\begin{slide}{}
\begin{center}
Outline
\end{center}
\begin{description}
\item
What is semidefinite programming, SDP
(quick)
\item
Complementarity Conditions\\
 and Symmetrization
\item
Ill-conditioning
\item
Search directions
\item
Lots of questions\\
 (from a numerical analysis point of view)
\end{description}
\end{slide}
\begin{slide}{}
\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 constr.}\\
  && X \succeq 0,~~(X \in \p)& \mbox{nonneg constr.}
    \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}.

\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 -requires 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
\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{complementary slackness}
\end{array}
\]

\[
    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.\\
(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}{}
{\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^{\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$

\[ (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)
\]
\hspace{.5in}so that $(X^{k+1},Z^{k+1}) \succ 0$\\
{\bf end (for)}.

\end{slide}
\begin{slide}{}
Complications:

Overdetermined system, i.e.
how to 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}{}

{\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}

