\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)\\ &\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}