\begin{filecontents}{seceqn.sty}
% seceqn.sty
% From TeXhax 89.37
%
% Substyle file for use with "article" to cause equations to be numbered
% within sections.
%
% 13 Apr 89	Jerry Leichter
\typeout{Document Option `seceqn':  13 Apr 89}
\@addtoreset{equation}{section}		%Make equation=0 when section steps
\def\theequation{\thesection.\arabic{equation}}
					%How an equation number looks

% From NUMINSEC.STY (with SIAM stuff).
\@addtoreset{figure}{section}
\def\thefigure{\thesection.\@arabic\c@figure}

\@addtoreset{table}{section}
\def\thetable{\thesection.\@arabic\c@table}
\end{filecontents}

\documentclass[a4,landscape]{seminar}
%\documentclass[a4,portrait]{seminar}
\usepackage{semcolor}
\input{seminar.bug}
%\usepackage{theorem,exscale}
%\usepackage{fullpage}        % was not available on PC?
%\usepackage{float}
%\usepackage{psfrag}
%\usepackage{algorithmic}     % was not available on PC?
%\usepackage{hyperref}
%\usepackage{seceqn}
\usepackage{epsfig}
\usepackage{amssymb}
%\usepackage{amssymb,epsfig}
%\usepackage{amsmath,amsthm,amsfonts,latexsym,amssymb,epsfig}
%\usepackage{fullpage}        % was not available on PC?
\usepackage{float}
%\usepackage{algorithmic}        % was not available on PC?
%\usepackage{epsf}
\usepackage{hyperref}
\usepackage{seceqn}


\newenvironment{pf}{\parindent=0pt{\textbf{Proof: }}}{\hfill
  $\blacksquare$ \\} 
\newcounter{ctr}
\newtheorem{thm}{Theorem}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{cor}{Corollary}[section]
\newtheorem{con}{Conjecture}[section]
\newtheorem{prop}{Proposition}[section]
\newtheorem{defi}{Definition}[section]
\newtheorem{example}{Example}[section]
\newtheorem{rem}{Remark}[section]
\newtheorem{alg}{Algorithm}[section]
\newtheorem{ex}{Exercise}[section]
\newcommand{\adj}{{\rm adj\,}}
\newcommand{\trace}{{\rm trace\,}}
\newcommand{\spanl}{{\rm span\,}}
\newcommand{\tr}{{\rm trace\,}}
\newcommand{\Rn}{\mathbb{R}^{n}}
\newcommand{\Sn}{{\mathcal S^n\,}}
\newcommand{\Mn}{{\mathcal M^n\,}}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\beqr}{\begin{eqnarray}}
\newcommand\C{\mathbb C}
\newcommand{\NN}{\mathcal N}
\newcommand{\A}{\mathcal A}
\newcommand{\RR}{\mathbb R}
\newcommand{\benum}{\begin{enumerate}}
\newcommand{\eenum}{\end{enumerate}}
\newcommand{\st}{\textnormal{s.t.}}
\newcommand{\TRSe}{TRS$_=\;$}
\newcommand{\disp}{\displaystyle}
\newcommand{\QED}{\hfill ~\rule[-1pt] {8pt}{8pt}\par\medskip ~~}
\newcommand{\bpr}{\textsc{proof: }}
\newcommand{\epr}{\QED}
\newcommand{\Diag}{{\rm Diag\,}}
\newcommand{\diag}{{\rm diag\,}}

\newcommand{\offDiag}{{\rm offDiag\,}}
\newcommand{\usMat}{{\rm us2Mat\,}}
\newcommand{\usvec}{{\rm us2vec\,}}
\newcommand{\svec}{{\rm s2vec\,}}
\newcommand{\kvec}{{\rm vec\,}}
\newcommand{\ZS}{{\mathcal Z_S} }
\newcommand{\XS}{{\mathcal X} }
%\newcommand{\XS}{{\mathcal X_s} }
\newcommand{\XSu}{{\mathcal X_{u}} }
\newcommand{\XSd}{{\mathcal X_{d}} }
\newcommand{\XX}{{\mathcal X} }
%\newcommand{\SO}{{\mathcal S^s} }
\newcommand{\SO}{{\mathcal S} }
\newcommand{\SOu}{{\mathcal S} }
%\newcommand{\SOu}{{\mathcal S^{su}} }
\newcommand{\Ss}{{\mathcal S} }
\newcommand{\Stn}{{{\mathcal S}^{\scriptsize{\pmatrix{n\cr 2}}}\,}}
\newcommand{\tn}{{\scriptsize{\pmatrix{n\cr 2}}}\,}
\newcommand{\Rtnp}{{\R^{\scriptsize{\pmatrix{n+1\cr 2}}}\,}}





% Nick's definitions.
\def\Rnbyn{\mathbb{R}^{n\times n}}
\def\R{\mathbb{R}}
\def\normF#1{\|#1\|_F}
\def\normt#1{\|#1\|_2}
\def\norm#1{\|#1\|}
\def\cm{correlation matrix}
\def\cms{correlation matrices}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All refs in roman.
\def\eqref#1{{\normalfont(\ref{#1})}}

\begin{document}
%\bibliographystyle{plain}





\begin{slide}


\begin{center}
{\bf 
Large Sparse Semidefinite Programming\\ with Applications to\\ the Nearest
Correlation Matrix Problem
}
\end{center}

\vspace{3mm}

     Henry Wolkowicz \\
Department of Combinatorics \& Optimization \\
University of Waterloo

\mbox{
\begin{figure}
\psfig{file=UWlogori.ps,height=20mm}
\end{figure}
}

(M. Anjos, N. Higham, M. Takouda; (H. Wei))





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

\begin{center}
{\bf OUTLINE}
\end{center}
\begin{description}
\item $\bullet$
Background on SDP; Motivation
\item $\bullet$
Robust, `non-interior', path-following algorithm for SDP (and LP)\\
\qquad (Gauss-Newton direction, preconditioned conjugate gradients)
\item $\bullet$
Nearest Correlation Matrix Problem: duality and algorithm
\item $\bullet$
Numerics
\item $\bullet$
Other SDP models
\end{description}




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



\section{What is SDP? And Motivation}
\label{sect:unc}
\[
\mbox{(SDP)}\quad
\begin{array}{lcl}
& \min  & f(X)\\
&  \mbox{subject to} &  \A X=b \\
 &&   X \succeq 0,
\end{array}
\]
~~\\
where:  $f:\Sn \rightarrow \RR$ convex function\\
$~\qquad  \Sn$ $n \times n$ real symmetric matrices\\
$~\qquad  \A:\Sn \rightarrow \RR^m \mbox{ linear operator, }$\\
$X (\succeq) \succ 0$ denotes positive (semi)definite

$\left(\quad  (\A X)_i = \left< A_i,X\right>
                  = \trace A_i X,\quad  A_i=A_i^T, i = 1\ldots n)
                \quad  \right)$

\subsection{SDP Arises Naturally in Optimization}
\label{sect:natur}
\subsubsection{Quadratic Model for $\min_x f(x)$}
\label{sect:uncquad}
A quadratic model at the current estimate $x_c$:
\[
\mbox{(Quad)}\qquad
\begin{array}{lll}
q^* = & \min & q(d):=  f(x_c)+ \nabla f(x_c)^T d + \frac 12 d^T\nabla^2 f(x_c) d\\
& \mbox{s.t.} & \Vert d \Vert \leq s,  \qquad (\|d\|^2\leq s^2 )
\end{array}
\] 
where the normalization/constraint avoids unboundedness
and steplengths that are {\em too} long.

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

Lagrangian dual (with Hessian $\frac 12 A$ and gradient $-2 a$):
\begin{eqnarray*}
q^* = \nu^*&:=&
      \max_{\lambda\leq 0} \;\min_x x^T(A-\lambda I)x -2a^Tx + \lambda s^2\\ 
      &=& \max_{\lambda\leq 0} \; h(\lambda)
\end{eqnarray*}
{\bf hidden constraint}: inner minimization bounded below\\
$ \mbox{dom}(h)$ restricted to $\nabla^2L(x,\lambda)=A-\lambda I \succeq 0$

{\bf (convex) Lagrangian} is
$L(x,\lambda) := x^T(A-\lambda I)x -2a^Tx + \lambda s^2$;\\ 
{\bf (concave) dual functional} is
$h(\lambda) :=  \min_x \; L(x,\lambda)$\\

(This can be exploited to solve large sparse trust region subproblems,
e.g. classical More-Sorensen (1983) using Cholesky factorization;
$n>10^6$, cf. C. Fortin and W. (2002) using Lanczos/eigenvalues)


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

SDP has surprisingly many applications arising from
\begin{center}
 {\em quadratic models}
\end{center}
\begin{enumerate}
\item
Engineering (control theory, design)
\item
Hard Combinatorial Problems
\item
Robust Optimization
\item
Mathematics of Finance (below)
\item
{\bf $\ldots $}
\end{enumerate}



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


\section{(Non) Interior Path-Following}
\subsection{Illustration/Motivation on LP Case}

\beq \label{prob:PLP}
(\mbox{LP}) \qquad
        \begin{array}{ccl}
        p^* :=&\min   & c^Tx \quad (\mbox{or }\left<c,x\right>)\\
        &\mbox{s.t.} & Ax = b \\
        &     &  x \geq  0   \quad (\mbox{or } x \succeq 0)
        \end{array}
\eeq
\beq \label{prob:DLP}
(\mbox{DLP}) \qquad
        \begin{array}{ccl}
        d^* :=&\max   & b^Ty\\
        &\mbox{s.t.} & A^T y +z = c \\
        &     &  z \geq  0  \quad (\mbox{or } z \succeq 0)
        \end{array}
\eeq
$A \in \Re^{m \times n}$ full rank (onto); LP and DLP strictly feasible


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

{\bf dual log-barrier problem} with parameter $\mu >0$ is
\[
(\mbox{Dlogbarrier}) \qquad
        \begin{array}{ccl}
        d_\mu^* :=&\max   & b^Ty+\mu \sum_{j=1}^n \log z_j \quad 
                                    (+\mu \log \det (z))\\
        &\mbox{s.t.} & A^T y +z = c \\
        &     &  z >  0 \qquad (z \succ 0).
        \end{array}
\]
stationary point of the Lagrangian / 
optimality conditions
\[
F_\mu(x,y,z)=
 \left(
\begin{array}{cl}
A^T y +z - c \\
A x -b    \\
X-\mu Z^{-1}
\end{array}
\right)  =0,\\ 
\begin{array}{cl}
\quad x,z>0, \quad (\succ 0)\\ X=\Diag(x),\\Z=\Diag(z)
\end{array}
 \]

{\em central path}:
set of these solutions
$(x_\mu,y_\mu,z_\mu), \mu > 0$



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


As $\mu \rightarrow 0$,  Jacobian $F_\mu^\prime(x,y,z)$
grows ill-conditioned near central path

{\bf Cure/Fix:} Make nonlinear equations {\em less nonlinear},
i.e. preconditioning for Newton type methods;\\
premultiply by block-diag matrix with blocks $(I,I,Z)$:
\[
\begin{array}{rcl}
F_\mu(x,y,z) \leftarrow
\pmatrix{I&0&0\cr 0&I&0 \cr 0&0&Z}
F_\mu(x,y,z)
&=&
 \left(
\begin{array}{cl}
A^T y +z - c \\
A x -b    \\
ZX-\mu I
\end{array}
\right) \\
& =:& 
 \left(
 \begin{array}{cl}
R_d \\
r_p    \\
R_{ZX}
\end{array}
\right) 
\end{array}
 \]



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

Special structure of linearized system can be exploited;
linearization for the Newton
direction $\Delta s=\pmatrix{\Delta x \cr \Delta y \cr \Delta z}$
is \beq \label{eq:newtoneqn} F^{\prime}_\mu(x,y,z) \Delta s=
\pmatrix{ 0 & A^T & I \cr A & 0 & 0    \cr Z & 0 & X } \Delta s =
-F_\mu(x,y,z). \eeq

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

overdetermined system in SDP case: \\
\[
\Sn \times \Re^m \times \Sn \rightarrow  \Sn \times \Re^m \times \Mn
\]
 apply symmetrization 'undoes preconditioning'


\beq \label{eq:newtoneqn2} 
\pmatrix{ I & 0 & 0 \cr 0 & I & 0    \cr 0 & 0 & {\cal S} } 
\pmatrix{ 0 & A^T & I \cr A & 0 & 0    \cr Z & 0 & X } 
\eeq

e.g. last equation is linearization of:\\
$ZX+XZ - 2 \mu I = 0$ (AHO search direction)



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


\subsection{Reduction/Block Elimination for the Normal Equations}
\label{sect:rednormal}
Step 1 (Eliminate $\Delta z$):
\beq \label{eq:elimdz}
 \pmatrix{
I & 0 & 0 \cr 0 & I & 0    \cr -X & 0 & I }
 \pmatrix{
0 & A^T & I \cr A & 0 & 0    \cr Z & 0 & X }
=
\left(
 \begin{array}{ccc}
0 & A^T & I \cr  A & 0 & 0    \cr Z & -XA^T & 0
\end{array}
\right).
\eeq
We let \beq \label{eq:PKmats} P_Z=
 \pmatrix{
I & 0 & 0 \cr 0 & I & 0    \cr -X & 0 & I }, \quad K= \left(
 \begin{array}{ccc}
 0 & A^T & I \cr
A & 0 & 0 \cr Z & -XA^T & 0
\end{array}
\right). \eeq


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

with right-hand side
\[
 -\pmatrix{
I & 0 & 0 \cr 0 & I & 0    \cr -X & 0 & I }
\left(
\begin{array}{cl}
R_d \\ r_p    \\ R_{ZX} -\mu e
\end{array}
\right)
=
\left(
\begin{array}{cl}
-R_d \\ -r_p    \\ XR_d  -R_{ZX}
\end{array}
\right)
 \]



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

Step 2 (Eliminate $\Delta x$):
\beq
\begin{array}{rcl}
\label{eq:elimdx}
\vspace{.1in}
F_n:= P_nK & :=  & \pmatrix{ I & 0 & 0 \cr 0 & I &
-AZ^{-1} \cr 0 & 0  & Z^{-1} } \left(
 \begin{array}{ccc}
0 & A^T & I \cr A & 0 & 0 \cr Z & -XA^T & 0
\end{array}
\right)  \\
&=&
\left(
 \begin{array}{ccc}
0 & A^T & I_n \cr 0 & AZ^{-1}XA^T & 0  \cr I_n & -Z^{-1}XA^T & 0
\end{array}
\right)
\end{array}
\eeq
$AZ^{-1}XA^T$ can have:\\
 $\bullet$ uniformly bounded condition number, e.g. G{\"u}ler et al 1993\\
 $\bullet$ structured singularity, M. Wright 1997\\
But ${\rm cond}(F_n) \rightarrow \infty$.



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

The right-hand side becomes
\[
\begin{array}{l}
\vspace{8mm}
-P_n P_Z
\left(
\begin{array}{cl}
R_d \\ r_p    \\ R_{ZX}
\end{array}
\right)
=\\
\qquad \qquad 
\left(
\begin{array}{cl}
-R_d \\
-r_p  +A(x-Z^{-1}XR_d  -\mu Z^{-1} e)  \\
Z^{-1}XR_d -x +\mu Z^{-1} e
\end{array}
\right)
\end{array}
 \]



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



\begin{prop}
The condition number of $F_n^TF_n$ diverges to infinity 
if  $x(\mu)_i/z(\mu)_i $ diverges to infinity,
for some $i$, as $\mu$ converges to 0. The condition number of 
$(F_\mu^{\prime})^TF_\mu^{\prime}$ 
is uniformly bounded if there exists a unique primal-dual solution.
\end{prop}
\bpr
Note that
\[
F_n^TF_n =
 \pmatrix{
I_n & -Z^{-1}XA^T & 0 \cr
-AXZ^{-1} & (AA^T+(AZ^{-1}XA^T)^2+AZ^{-2}X^2A^T) & A  \cr
 0 & A^T & I_n
}.
\]
By interlacing of eigenvalues, ...
\epr



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

{\bf EXAMPLE:} {\em getting too close to boundary} (worse for SDP)
$ A=\pmatrix{1&1}, c=\pmatrix{-1 \cr 1}, b=1.  $\\
$ x^*=\pmatrix{1 \cr 0}, y^*=-1, z^*=\pmatrix{0 \cr 2};$\\


initial points:
$
\begin{array}{lc}
x = \pmatrix{9.183012e-001\cr
             1.356397e-008},
z = \pmatrix{2.193642e-008\cr
             1.836603e+000},\\
y =         -1.163398e+000.
\end{array}
$

residuals and duality gap:
$
\|r_b\| = 0.081699, \quad
\|R_d\| = 0.36537, \quad \mu=x^Tz/n = 2.2528e-008
$

5 decimals rounding before/after arithmetic\\
centering with $\sigma = .1$

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


search directions 
 found using
full matrix $F^{\prime}_{\mu}$ and backsolve matrix $F_n$:
\[
\pmatrix{\Delta x\cr \Delta y\cr \Delta z}= \pmatrix{ 
   8.17000e-02  \cr
   -1.35440e-08   \cr
    1.63400e-01   \cr
   -2.14340e-08   \cr
    1.63400e-01   
           };  \quad = 
   \pmatrix{
       -6.06210e+01 \cr
       -1.35440e-08 \cr
         1.63400e-01 \cr
        0.00000e+00 \cr
         1.63400e-01 
             }
\]
error in $\Delta y$ is small; error after
backsubstitution for
$(\Delta x)_1$ is large.\\
\[
\pmatrix{AZ^{-1}XA^T \cr -Z^{-1}XA^T}=
\pmatrix{
    4.18630e+07 \cr
   -4.18630e+07 \cr
   -7.38540e-09
}
\]
%\end{example}

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

\subsection{Alternate Second Step; Stable Reduction}
{\bf Assuming!} $A=[I_m ~ E]$.\\
Partition diagonal matrix $Z,X$ using vectors $z=\pmatrix{z_m \cr z_v},
 x=\pmatrix{x_m \cr x_v}$, $XA^T=\pmatrix{X_m \cr X_v E^T}$
\[
\begin{array}{rcl}
\vspace{.1in}
 F_s:&=&P_sK=\pmatrix{
I_n & 0 & 0 & 0  \cr 0 & I_m & 0 & 0  \cr 0 & -Z_m & I_m &  0
\cr 0 & 0 & 0 & I_v }
 \pmatrix{
        0   &  0   &       A^T         &      I_n      \cr
I_m & E & 0    &     0  \cr
 Z_m & 0   &  -X_m &  0  \cr
 0 & Z_v   &  -X_vE^T &  0
}\\
&=&
 \pmatrix{
0  & 0  & A^T & I_n \cr I_m & E & 0 & 0 \cr 0 & -Z_mE &  -X_m  & 0
\cr 0  & Z_v  & -X_vE^T  & 0 }.
\end{array}
\]


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

The right-hand side becomes
\[
\begin{array}{rcl}
-P_s P_Z
\left(
\begin{array}{cl}
A^T y +z - c \\ A x -b    \\ ZXe -\mu e
\end{array}
\right)
=
\vspace{.5mm}
-P_s
\left(
\begin{array}{cl}
R_d \\ r_p    \\ -XR_d +ZXe -\mu e
\end{array}
\right)\\
\qquad=
\left(
\begin{array}{cl}
-R_d \\ 
-r_p\\
-Z_mr_p
-X_m(R_d)_m +Z_mX_me -\mu e\\
-X_v(R_d)_v +Z_vX_ve -\mu e
\end{array}
\right)
\end{array}
 \]

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

\subsection{Solving the Last Two Rows}
We can
change to a symmetric indefinite system by changing the rows and
scaling by the diagonal matrix
$\pmatrix{X_v & 0 \cr 0 & Z_m}$,
\beq
\label{eq:symindef}
\pmatrix{Z_vX_v & -X_vE^TZ_m \cr -Z_mEX_v & -X_mZ_m}.
\eeq
If $E$ is sparse, then this system stays sparse after the scaling.


But -- we use a conjugate gradient type method below.



\end{slide}
\begin{slide}{}
\section{Finding the Nearest Correlation Matrix using SDP}

{\bf Given} symmetric matrix $A\in\Sn$:

\[
  \mu^*= \min  \frac12  \normF{A-X}^2 \quad
    \mbox{s. t. \quad $\diag X =e$, $X \succeq 0$, $X \in \Sn$}.
\]


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


Applications  e.g. Finance:\\
 {\bf approximate correlation matrices} are constructed 
from vectors of stock returns; sample correlations are
taken only from days on which both stocks have data available
(from inconsistent data).\\
Ref: N. Higham (2002) uses a projection technique to take advantage of 
low rank.

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

\subsection{Mixed-Cone Formulation}

direct approach using
a mixed SDP and second-order (or Lorentz) cone problem:
\beq
\label{Prob-SOCP}
\begin{array}{rcl}
 \sqrt{2\mu^*} = &\min &\alpha \\
     &\mbox{s.t.}& \diag X = e\\
         &&  Y+X=A,~ \normF{Y} \le \alpha \\
       &&    X,Y \in \Sn, X \succeq 0.
\end{array}
\eeq
(Public domain software packages are available - numeric
comparison below.)


\end{slide}
\begin{slide}{}
\subsubsection{Operator Notation: $\usvec,  \usMat $}

\[x=\usvec(X) \in \R^{n \choose 2}
\qquad (a=\usvec A, s=\usvec S)
\]
$\sqrt{2}$ times vector (columnwise) from 
strict upper-triang of $X$.

${n\choose 2}=n(n-1)/2$;
$\sqrt{2}$ guarantees isometry.

$\usMat:=\usvec^{\dagger}$ mapping into $\Sn$\\
adjoint transformation $\usMat^*= \usvec$ since:
\begin{eqnarray*}
\left<\usMat(v),S\right> &=& \trace \usMat(v)S
= \trace \usMat(v) \offDiag(S)\\
&=&  v^T \usvec(S) = \left<\usvec(S),v\right>
\end{eqnarray*}
orthogonal projection: $\offDiag(S)=\usMat \usvec (S)$


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

\subsection{Duality and Optimality Conditions}

(using $X = \usMat(x)+I$)
an equivalent problem is:
\[
   \mu^* := \min  \frac 12 \normt{x-a}^2
     \quad\mbox{subject to \quad $\usMat(x) + I \succeq 0$,
     $x \in \R^{n\choose 2}$}.
\]
strong (Lagrangian) duality:
\[
\mu^*=\nu^*:= \max_{S \succeq 0}
\min_x  \frac 12  || x-a ||^2- \trace S(\usMat(x) + I ).
\]


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

change to {\bf Wolfe dual} using stationarity of inner minimization:
\[
0=(x-a) - \usMat^*(S) =(x-a) -  \usvec(S),
\]
and since $\trace S(\usMat(x) + I )= x^T\usvec(S)+\trace S$.


\begin{equation}
   \label{Prob-P2-D}
\begin{array}{rll}
  \mu^*=&\max & \frac 12  || x-a ||^2 - \trace S(\usMat(x) + I )\\
  & \mbox{subject to}& x - \usvec(S) = a\\
        &&   S \succeq 0.
\end{array}
\end{equation}




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

Slater's CQ holds for both primal and dual:
\begin{thm}
The optimal values $\mu^*=\nu^*$ and the
primal-dual pair $x,(y,s)$ are optimal if and only if
\[
\begin{array}{rcll}
X&:=&\usMat(x) +I  \succeq 0    \quad  \mbox{(primal feasibility)}\\
   x&=&a + s,~~S^y:=\usMat(s)+\Diag(y) \succeq 0  \quad  \mbox{(dual feasibility)}\\
XS^y &=& 0   \quad  \mbox{(complementary slackness)}\\
\end{array}
\]
\epr
\end{thm}

(cf $A=\left[I ~ E\right]$ for LP)



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

Substitute feasibility equations; exact primal-dual feasibility
maintained during iterations; full rank Jacobian at optimality.\\

{\em single bilinear (perturbed) equation} in $s,y$:
    $\quad F_{\mu}(s,y): \R^{n+1 \choose 2} \rightarrow \Mn$
\fbox{
\begin{minipage}{4.05in}
\[
\begin{array}{c}
\vspace{.5mm}
F_{\mu}(s,y):=\left[A+\usMat(s)+I \right] \left[\usMat(s)
                    +\Diag(y) \right]   - \mu I = 0
\end{array}
\]
\end{minipage}
}

typical SDP - 
overdetermined system of bilinear equations;\\
current approach is to symmetrize - which results in
ill-conditioning! from rank deficient Jacobian at optimality.\\

BUT, here, no symmetrization used;\\
solve using (an inexact) Gauss-Newton method. 



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


Linearization for search direction
$\Delta v=\pmatrix{\Delta s \cr \Delta y}$:
($S^y=\usMat(s)+\Diag(y)$)
\[
\begin{array}{lcl}
F^{\prime}_{\mu}(s,y) \Delta v =\\
=\left[A+\usMat(s)+I \right]\left(\usMat( \Delta s) +\Diag(\Delta y)
\right)
   + \usMat(\Delta s)S^y \\
=\left(\XSu+\SOu\right) (\Delta s) + \XSd (\Delta y).
\end{array}
\]
This is a linear, full rank, overdetermined system.\\
Our search direction $\Delta v$ is its least squares solution.



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

\subsubsection{Algorithm: p-d i-e-p framework}
{\bf $\bullet$ Initialization:}
\begin{quote}
{\bf $\bullet\bullet$  Input data:} a real symmetric $n\times n$ matrix $A$
(set $\diag(A)=0$)\\
{\bf $\bullet\bullet$ Positive tolerances:}
\[\epsilon_1 \mbox{ (stopping), } \epsilon_2 \mbox{ (lss accuracy), } \epsilon_3
\mbox{ (crossover), }
\]
\\
{\bf $\bullet\bullet$ Find initial strictly feasible points:}
 both
$S^0,X^0:=\left(\offDiag (S+A) +I \right)~\succ~0$; $\mu>0$\\
{\bf $\bullet\bullet$ Set initial parameters:}
\[
   {\rm gap}=\trace S^0X^0;~~ \mu={\rm gap}/n;~~{\rm objval}=.5||X^0-A||_F^2;~~k=0.
\]
\end{quote}
{\bf $\bullet$ while} 
$\min \{ \frac {\rm gap}{{\rm objval}+1} , {\rm objval} \} >
              \epsilon_1$
\begin{quote}
{\bf $\bullet$$\bullet$  solve lss for search direction}
(accuracy $\epsilon_2\min\{\mu,1\}$)
     \[
  ~~~~F^{\prime}_{\sigma \mu}(v^k)
\left( \Delta v^k \right)
= -F_{\sigma \mu}(v^k),
\]
where $\sigma_k$ centering,
        $\mu_k=\frac{1}{n}\trace S^k(\offDiag(S^{k+1}+A)+I)$
\[ S^{k+1} = S^k + \alpha_k \Delta S^k,  ~~ \alpha_k > 0,
\]
so that both $S^{k+1},\offDiag(S^{k+1}+A)+I \succeq 0$\\
($\alpha_k=1$ after the crossover.)\\
{\bf $\bullet$$\bullet$  update}
\[ k \leftarrow k+1 \quad \mbox{ and then}
\]
\[  \sigma_k \quad \left(\mbox{set } \sigma_k=0 \mbox{ if }
\min \{ \frac {\rm gap}{{\rm objval}+1} , {\rm objval} \} < \epsilon_3\right)
\]
\end{quote}
{\bf $\bullet$ end (while)}.\\
{\bf $\bullet$ Conclusion: } $X$ is approx. $\usMat (s)+A +I$

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

After the {\bf crossover}, centering $\sigma=0$ and steplength $\alpha=1$,
we get q-quadratic convergence; allows for {\em warm starts}.

{\bf Long steps} can be taken {\em beyond} the positivity boundary.
(tests show improved convergence rates)

Perturbations (remove smaller nonzeros)
to given $A$ are done at start to reduce {\em nnz($A$)}.
Nonzeros are added back to $A$ once the duality gap is {\em small}. 
Tests show that q-quadratic convergence is maintained.

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


\subsection{Preconditioning}

\[
\left(\XSu+\SOu\right) P_s^{-1}(\widehat{\Delta s}) + \XSd P_y^{-1}
            (\widehat{\Delta y})  = - F_{\mu}(s,y)
\]
where 
\[
\widehat{\Delta s}=P_s(\Delta s),\quad
\widehat{\Delta y}=P_y(\Delta y)
\]


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

\subsubsection{Diagonal Preconditioning}
Optimal scaling Dennis and W. (1993)
full rank matrix $A\in\R^{m \times n}$, $m \ge n$,
with condition number
$\omega (K) := n^{-1}\trace (K)/\det(K)^{1/n}$,
the optimal scaling
\begin{equation}
\min \omega ( (AD)^T(AD)) \quad \mbox{subject to: $D$ positive and diagonal}
\end{equation}
solution: $d_{ii}= 1/\normt{A_{:i}}, i=1,\ldots,n$


{\em explicit} expressions for preconditioner\\
inexpensive


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

two {\em diagonal} operators $P_s,P_y$\\
evaluate using columns of $F_\mu^\prime(s,y)$. \\
$k\cong (i,j),~1\leq i < j\leq n$,
strictly upper triangular part of $S$\\
$i=1,\dots,n$ for $y$\\
$X=A+\usMat(s)+I$ and $S=\usMat(s)+\Diag(y)$
\[
\begin{array}{rcl}
\XSu(e_k)
&=&
X \usMat(e_k)\\
&=&
\frac 1{\sqrt{2}}
X
\left( e_ie_j^T+e_je_i^T \right)
\\
&=&
\frac 1{\sqrt{2}} \left\{
\left( X_{:i}e_j^T+X_{:j}e_i^T \right)
\right\}.
\end{array}
\]
\[
\begin{array}{rcl}
\SOu(e_k)
&=&
\usMat(e_k)(S+\Diag(y))\\
&=&
\frac 1{\sqrt{2}}
\left( e_ie_j^T+e_je_i^T \right)(S+\Diag(y))
\\
&=&
\frac 1{\sqrt{2}} \left\{
\left( e_i(S+\Diag(y))_{j:}+e_j(S+\Diag(y))_{i:} \right)
\right\}.
\end{array}
\]


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


Therefore
\[
\begin{array}{lcl}
\| (\XSu+\SOu)(e_k) \|^2_F=\\
\qquad =
\frac 12 \left\{ \| (S+\Diag(y))_{:i} \|^2+ \| (S+\Diag(y))_{:j} \|^2
       + \right.\\
  \qquad\qquad
 \left. \| X_{:i} \|^2+ \| X_{:j} \|^2
+2(S+\Diag(y))_{jj}X_{ii}\right. \\
 \qquad \qquad
\left. +4(S+\Diag(y))_{ji}X_{ij}
+2(S+\Diag(y))_{ii}X_{jj}
\right\}.
\end{array}
\]
need three Hadamard products\\
$X \circ X,  (S+\Diag(y)) \circ (S+\Diag(y)),  (S+\Diag(y)) \circ X$\\
 and {\em vector} Kronecker
product $\diag((S+\Diag(y))) \otimes \diag(X)$.

\[
\XSd(e_i) = X\Diag(e_i)
\]
Therefore
\[
\|\XSd(e_i)\|^2_F = ||X_{i,:}||^2.
\]


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

\subsubsection{Block-Diagonal Incomplete Cholesky Preconditioner}
natural block structure:
\[
\left[ \left(\XSu+\SOu\right) ~|~   \XSd \right] 
   \left( \begin{array}{c}\Delta s\\\Delta y\end{array} \right)
                       = -F_{\mu},
\]
normal equations have the block structure
\[
\begin{array}{l}  
\left[\begin{array}{c|c}  
 \left(\XSu^*+\SOu^*\right) \left(\XSu+\SOu\right) &  
                                \left(\XSu^*+\SOu^*\right) \XSd \\
\hline
 \XSd^* \left(\XSu+\SOu\right) &  
                                 \XSd^*\XSd \\
                 \end{array}  \right] 
   \left( \begin{array}{c}\Delta s\\\Delta y\end{array} \right)
=\\
\qquad\qquad \qquad\qquad  =-\left( 
    \begin{array}{c} \XSu^*+\SOu^* \\   \XSd^*  \end{array}
     \right) 
                       F_\mu
    \end{array}
\]

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


preconditioner based on partial Cholesky factorizations
of the  block diagonal positive definite operator:
$$
\widetilde{P}^* \widetilde{P} = \left[\begin{array}{c|c}  
 \left(\XSu^*+\SOu^*\right) \left(\XSu+\SOu\right) &  0 \\
\hline
  0&                                  \XSd^*\XSd \\
                 \end{array}  \right].
$$
$\widehat{P}^* \widehat{P}$ has the approximate  factorization
$$
 \left[\begin{array}{c|c}  
 \widehat{R}^T \widehat{R} &   \\
\hline
  &         D^2  \\
                 \end{array}  \right] 
$$

can exploit the special structure again


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

\section{Numerical Tests}
Pentium 4; MATLAB 6.5; 1 GIG RAM.

crossover heuristic: relative duality gap $< .1$.

Initially , zero out all elements $|A_{ij}|< \epsilon_4=.5$.\\
In each iteration with  
(relative duality gap) $< .003$, set
$\epsilon_4 \leftarrow (\epsilon_4 -.1)$

Stopping criteria (relative duality gap) $< \epsilon_1 = 1e-10$.\\
(But - average accuracy attained $1e-13$, q-quadratic convergence.)



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


\begin{figure}[htb]
\epsfxsize=200pt
\centerline{\epsfbox{cpu_comparison_graph.eps}}
\caption{comparison: 2nd order cone (SeDuMi) (avrg 10 each density)}
\label{fig:cpu_comparison}
\end{figure}




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




\begin{figure}[htb]
\epsfxsize=220pt
\centerline{\epsfbox{dens0050103.ps}}
\caption{density .0005:.001:.003, CPU times and nnz(X), n=200}
\label{fig:densityn200}
\end{figure}


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


\begin{figure}[htb]
\epsfxsize=220pt
\centerline{\epsfbox{dens300.ps}}
\caption{density  .0005:.001:.003 vs CPU times and nnz(X), n=300 }
\label{fig:densityn300}
\end{figure}


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


\begin{figure}[htb]
\epsfxsize=220pt
\centerline{\epsfbox{dens350.ps}}
\caption{density  .0005:.001:.027 vs CPU times and nnz(X), n=350 }
\label{fig:densityn350}
\end{figure}



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

\tiny{
\begin{verbatim}
random A dim n =  100
density  0.008
number of nonzeros  80
noiter   -log10(relgapp)    step        sigmaa     optval        cputime.dir     cputimechol    lsqr.iter  nnz-X
   1      2.79e-02         0.95             1     5.2083e+01    4.2000e-01       4.0000e-02         9       322
   2      6.38e-01         0.95         0.715     5.0593e+01    5.5000e-01       1.3000e-01        10       322
   3      1.23e+00         0.95         0.715     5.0125e+01    7.8000e-01       1.2000e-01        18       322
*******CROSSover started at: 0.05856
   4      2.85e+00            1             0     5.0026e+01    1.0900e+00       1.1000e-01        33       322
   5      3.79e+00            1             0     5.0008e+01    1.5300e+00       1.1000e-01        49       322
   6      4.48e+00            1             0     5.0003e+01    1.6100e+00       1.2000e-01        52       322
   7      5.27e+00            1             0     5.0002e+01    1.5700e+00       1.1000e-01        51       322
   8      6.53e+00            1             0     5.0002e+01    1.4200e+00       1.0000e-01        44       322
   9      8.93e+00            1             0     5.0002e+01    6.1000e-01       1.1000e-01        12       322
  10      1.37e+01            1             0     5.0002e+01    8.7000e-01       9.0000e-02        23       322
 final norm(XSy) is: 1.8144e-14
 rel. min of eigs is: -4.2533e-15
\end{verbatim}
}


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

\tiny{
\begin{verbatim}
start with  elements abs(A) > 0.4
and number of nonzeros in A: 24
noiter   -log10(relgapp)    step        sigmaa        optval    cputime.dir      cputimechol    lsqr.iter  nnz-X
   1      2.52e-01         0.855            1     5.1840e+01    3.8000e-01       3.0000e-02         9       126
   2      5.39e-01         0.95        0.7435     5.0492e+01    2.7000e-01       3.0000e-02         8       126
   3      1.16e+00         0.95         0.715     5.0110e+01    4.3000e-01       2.0000e-02        17       126
*******CROSSover started at: 0.069726
   4      2.80e+00            1             0     5.0020e+01    5.7000e-01       3.0000e-02        25       126
   5      3.77e+00            1             0     5.0004e+01    6.0000e-01       2.0000e-02        27       126
perturbing A - using elements abs(A) > 0.3
increasing number of nonzeros in A to:30
   6      4.40e+00            1             0     5.0001e+01    6.7000e-01       2.0000e-02        29       132
perturbing A - using elements abs(A) > 0.2
increasing number of nonzeros in A to:46
   7      5.01e+00            1             0     5.0000e+01    1.0200e+00       3.0000e-02        45       166
perturbing A - using elements abs(A) > 0.1
increasing number of nonzeros in A to:66
   8      3.89e+00            1             0     5.0001e+01    1.1100e+00       5.0000e-02        43       234
perturbing A - using elements abs(A) > 2.7756e-17
increasing number of nonzeros in A to:80
   9      4.76e+00            1             0     5.0003e+01    8.6000e-01       8.0000e-02        26       322
  10      5.72e+00            1             0     5.0002e+01    9.7000e-01       9.0000e-02        26       322
  11      7.35e+00            1             0     5.0002e+01    1.3200e+00       1.1000e-01        38       322
  12      1.06e+01            1             0     5.0002e+01    6.0000e-01       9.0000e-02        12       322
 final norm(XSy) is: 2.6748e-11
 rel. min of eigs is: -4.9538e-13
\end{verbatim}
}

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

\section{Conclusion}
Gauss-Newton direction:
\begin{description}
\item
Advantages/Disadvantages:
\begin{description}
\item
Robust, warm starts are simple, longer steps
\item
exact primal and dual feasibility at each iteration
\item
Can apply CG-type approaches
\item
q-quadratic convergence
\item
scale-invariant on the right
\end{description}
\item
Future:
\begin{description}
\item
Need large sparse QR efficient as Cholesky
\item
predictor-corrector
\end{description}
\end{description}


\end{slide}{}


%\bibliography{.master,.psd,.publs,.bjorBOOK}
\end{document}
