\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[landscape,slideColor,colorBG]{seminar} \documentclass[autumn_m,distiller,nototal,slideColor,colorBG]{prosper} %\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} %\usepackage[dvips]{graphicx} \usepackage{graphics} \usepackage{graphicx} %\usepackage{amsmath,amssymb,amsthm} %\usepackage{subeqnarray} \usepackage{hyperref} \usepackage{epsfig} \usepackage{float} \usepackage{seceqn} \usepackage{amssymb} \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{\Mmn}{{\mathcal M^{n-1}\,}} \newcommand{\beq}{\begin{equation}} \newcommand{\eeq}{\end{equation}} \newcommand{\beqr}{\begin{eqnarray}} \newcommand\C{\mathbb C} \newcommand{\BB}{\mathcal B} \newcommand{\NN}{\mathcal N} \newcommand{\LL}{\mathcal L} \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{\sMat}{{\rm sMat\,}} \newcommand{\usvec}{{\rm us2vec\,}} \newcommand{\svec}{{\rm svec\,}} \newcommand{\kvec}{{\rm vec\,}} \newcommand{\ZS}{{\mathcal Z_S} } \newcommand{\XS}{{\mathcal X} } \newcommand{\WW}{{\mathcal W} } %\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{\Stn}{{\mathcal S^{n-1}\,}} \newcommand{\tn}{{\scriptsize{\pmatrix{n\cr 2}}}\,} \newcommand{\Rtnp}{{\R^{\scriptsize{\pmatrix{n+1\cr 2}}}\,}} \newcommand{\EDM}{{\bf{\rm EDM\,}}} \newcommand{\SDP}{{\bf{\rm SDP\,}}} % 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}{\href{http://www.math.uwaterloo.ca/navigation/CompMath/}{Computational Mathematics} } \begin{center} \mbox{ \begin{figure} \psfig{file=math-mast.ps,width=50mm,height=23mm} \end{figure} } \end{center} \begin{center} {\bf A Stable Iterative Method for Linear Programming } \end{center} \vspace{.5in} Monday, Feb. 14 2005 %% %%\end{slide} %%\begin{slide}{ %%\href{http://www.mathstat.uoguelph.ca/}{Department of Math. \& Stats.} %%} %%\begin{center} %%\mbox{ %%\begin{figure} %%\psfig{file=logobar-left-bg_uog.ps,width=50mm,height=23mm} %%\end{figure} %%} %%\end{center} %%\begin{center} %%{\bf A Stable Iterative Method for Linear Programming %%} %%\end{center} %%\vspace{.5in} %%Thursday, Oct. 14 2004 %% %% %% \end{slide} \begin{slide}{ A Stable Iterative Method \\ for Linear Programming } \vspace{3mm} Henry Wolkowicz \\ \epsfxsize=200pt \centerline{\epsfbox{newbanner.eps}} %\vspace{.1in} %\epsfxsize=80pt %\centerline{\epsfbox{UWlogo_pms.eps}} University of Waterloo \mbox{ \begin{figure} \psfig{file=UWlogori.ps,height=20mm} \end{figure} } (with: Maria Gonzalez-Lima, Hua Wei) \end{slide} \begin{slide}{ OUTLINE} \begin{description} \item $\bullet$ {\cyan Background} on LP and SDP;\\ ~~~~~~~~~~~~ {\cyan Notation} and {\cyan Motivation} \item $\bullet$ {\cyan Robust}, (`non-interior') path-following algorithm for LP \\ \hspace{2in} ~~~~~~~~~~~~~~~~~~(and SDP) \item $\bullet$ Details of algorithm on LP with {\cyan numerics} \vspace{1.1in} \item $\bullet$ (time permitting) Application to Nearest Euclidean Distance Matrix Problem, Numerics, Comparisons with a dual algorithm (time permitting) \end{description} \end{slide} \begin{slide}{Notation and Motivation} \[ {\cyan \mbox{(SDP)}\quad \begin{array}{lcl} & \min & f(x)\\ & \mbox{subject to} & \A x=b \\ && x \succeq 0, \end{array} } \] where: \begin{description} \item $f$ convex function, $x \in X$ (e.g. $\RR^n, \Sn$, where $\Sn$ is space of $n \times n$ real symmetric matrices \item $\A$ is a linear operator, $b \in \RR^m$ \item $x \succeq 0$ denotes nonnegativity (in $\Sn$ or in $\RR^n$) %\item %$\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)$ \end{description} \end{slide} \begin{slide}{LP Maximization Problem} \[ (\mbox{DLP}) \qquad \begin{array}{cclc} d^* :=&\max & b^Ty & \mbox{profit}\\ &\mbox{s.t.} & A^T y \leq c & \mbox{resource constraints}\\ && y \in \RR_+^m & \mbox{products} \end{array} \] $A=(a_{ij}) \in \Re^{m \times n}$ full rank; $a_{ij}$ - units resource $j$/unit product $i$;\\ $\sum_i a_{ij}y_i$ total units of resource $j$ used;\\ $x_j \geq 0$ { shadow price of unit resource } $j$\\ \[ L(y,x):= b^Ty + x^T (c-A^Ty), \qquad \mbox{{\cyan Payoff Function} $X$ to $Y$} \] Optimal (worst case) strategy for {\cyan player $Y$} is:\\ \[ d^*=\max_{y \geq 0} \min_{x \geq 0} L(y,x):= b^Ty + x^T (c-A^Ty) \] \end{slide} \begin{slide}{Linear Program - Feasible Set} \begin{center} \mbox{ \begin{figure} \psfig{file=feasoptlp.ps,height=60mm} \end{figure} } \end{center} \end{slide} \begin{slide}{LP Dual - Player $X$} {\cyan Payoff Function} \[ L(y,x)= b^Ty + x^T (c-A^Ty) = c^Tx + y^T(b-Ax) \] Optimal (worst case) strategy for {\cyan player $X$} is:\\ (max-min interchanged)\\ \[ d^* \leq p^*:= \min_{x \geq 0}\max_{y \geq 0} L(y,x)=c^Tx + y^T(b-Ax) \] {\cyan Hidden Constraint} $b-Ax \leq 0$. \[ (\mbox{LP}) \qquad \begin{array}{cclc} p^* =&\min & c^Tx & \mbox{resource costs}\\ &\mbox{s.t.} & Ax \geq b & \mbox{be competitive}\\ & & x \geq 0 & \mbox{shadow prices} \end{array} \] \end{slide} \begin{slide}{Complementary Slackness} For primal-dual feasible pair $y,x$\\ ($x\geq 0, y\geq 0, A^Ty \leq c, Ax \geq b$) \[ \begin{array}{rcl} b^Ty &\leq& L(y,x)= b^Ty + x^T (c-A^Ty)\\ &=& c^Tx + y^T(b-Ax)\\ &\leq& c^Tx \end{array} \] with optimality $b^Ty^*=c^Tx^*$ if and only if \[ x^T (c-A^Ty)= y^T(b-Ax)=0 \qquad \mbox{\cyan Complementary Slackness} \] \end{slide} \begin{slide}{History/Applications of LP} \begin{description} \item developed during 1940's \\ \item motivation: complex planning problems in wartime operations \item George B. Dantzig - simplex method 1947 (based on John von Neumann theory of duality ) \item Nobel prize econonmics 1975: Leonid Kantorovich (USSR) and Tjalling Koopmans (USA) for contributions to {\em theory of optimal allocation of resources} \end{description} \end{slide} \begin{slide}{History/Applications of LP cont...} \begin{description} \item important applications e.g.: airline crew scheduling, shipping or telecommunication networks, oil refining and blending, and stock and bond portfolio selection.\\ (pre-wordprocessing SIAM survey - 70\% of all computer time) \item Leonid Khaciyan 1979 ellipsoid method shows solvability in a number of steps which is a polynomial function of the amount of data. \item Narendra Karmarkar 1984 - interior-point method Interior-point methods now generally competitive with simplex method. \end{description} \end{slide} \begin{slide}{(Non) Interior Path-Following on LP} \[ (\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} \] \[ (\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} \] $A \in \Re^{m \times n}$ full rank (onto); LP and DLP strictly feasible $p^* \geq d^*$ {\cyan Weak Duality} $p^* = d^*$ {\cyan Strong Duality} $z \cdot x = 0$ Complementary Slackness (elementwise product is 0) \end{slide} \begin{slide}{Linear Primal-Dual Pair of SDPs} (looks/behaves like Linear Program, LP) \[ \mbox{(PSDP)}\quad \begin{array}{lcl} & \min & \left=\trace CX\\ & \mbox{subject to} & \A X=b \\ && X \succeq 0, \end{array} \] \[ \mbox{(DSDP)}\quad \begin{array}{lcl} & \max & b^Ty\\ & \mbox{subject to} & \A^* y+Z =C \\ && Z \succeq 0, \end{array} \] $A_i, \quad A_i \in \Sn$ fixed $(\A X)_i = \trace A_iX$; {\em adjoint operator}: $ \A^* y = \sum^m_{i=1} y_i $ $ZX=0$ Complementary Slackness (matrix product) \end{slide} \begin{slide}{(Perturbed) Optimality Conditions} For {\em \cyan barrier parameter} $\mu > 0$: \[ \begin{array}{rcl} F_{\mu}(X,y,Z):=\pmatrix{ \A^* y + Z - C \cr \A X -b \cr ZX-\mu I } = 0 \qquad \pmatrix{ \mbox{dual feasibility} \cr \mbox{primal feasibility} \cr \mbox{pert. compl. slack.} } \end{array} \] For SDP:\\ $\qquad\qquad F_\mu: \Sn \times \Re^m \times \Sn \rightarrow \Sn \times \Re^m \times \Mn$ \\ i.e. overdetermined nonlinear system \end{slide} \begin{slide}{Central Path/Path Following} {\bf \cyan 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)= \pmatrix{ A^T y +z - c \cr A x -b \cr X-\mu Z^{-1}} =0, \qquad x,z>0, \quad (\succ 0) \] $X=\Diag(x), Z=\Diag(z)$ {\cyan \em central path}: set of these solutions $(x_\mu,y_\mu,z_\mu), \mu > 0$ \end{slide} \begin{slide}{Ill-Conditioning} Due to $Z^{-1}$ and $Z^*$ singular: as $\mu \rightarrow 0$, Jacobian $F_\mu^\prime(x,y,z)$ grows {\cyan ill-conditioned} near central path {\bf Cure/Fix:} Make nonlinear equations {\cyan \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}{Linearization} \[ \begin{array}{rcl} F_\mu(x,y,z) = \left( \begin{array}{cl} A^T y +z - c \\ A x -b \\ ZX-\mu I \end{array} \right) = 0 \end{array} \] Special structure of linearized system can be exploited;\\ {\cyan linearization} for the {\cyan Newton direction} $\Delta s=\pmatrix{\Delta x \cr \Delta y \cr \Delta z}$ is \[ 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). \] \end{slide} \begin{slide}{Damped Newton Method} {\cyan Damped} Newton steps \[ x \leftarrow x+ \alpha_p \Delta x, \quad y \leftarrow y+ \alpha_d \Delta y, \quad z \leftarrow z+ \alpha_d \Delta z, \] {\em backtrack} from nonnegativity boundary to maintain positivity/interiority, $x>0,z>0$.\\ {\cyan But:} on central path, $F_\mu(x,y,z) = 0$, \[ \mu = \frac 1n \mu e^Te =\frac 1n e^TZXe = \frac 1n z^Tx = \frac 1n (\mbox{ duality gap}), \] \[ \begin{array}{rcl} \mbox{barrier parameter } \mu \cong \mbox{\cyan duality gap} &=& c^Tx-b^Ty \\&=& x^T\left(c-A^Ty\right) \\&=& x^Tz. \end{array} \] {\cyan if exact feasibility holds!} \end{slide} \begin{slide}{SDP Case} {\cyan overdetermined} system in SDP case: \\ \[ \Sn \times \Re^m \times \Sn \rightarrow \Sn \times \Re^m \times \Mn \] apply symmetrization 'undoes preconditioning' \[ \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 } \] e.g. last equation is linearization of:\\ $ZX+XZ - 2 \mu I = 0$ (AHO search direction) \end{slide} \begin{slide}{Reduction/Block Elimination for the Normal Equations, NEQ} {\cyan Step 1} (Eliminate $\Delta z$): \[ \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). \] We let \[ P_Z= \pmatrix{ I & 0 & 0 \cr 0 & I & 0 \cr -X & 0 & I }, \quad {\cyan K= \left( \begin{array}{ccc} 0 & A^T & I \cr A & 0 & 0 \cr Z & -XA^T & 0 \end{array} \right)}. \] \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}{Normal Equations, NEQ} {\cyan Step 2} (Eliminate $\Delta x$): \[ \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 & {\cyan AZ^{-1}XA^T } & 0 \cr I_n & -Z^{-1}XA^T & 0 \end{array} \right) \end{array} \] {\cyan $AZ^{-1}XA^T$} can have:\\ $\bullet$ uniformly bounded condition number, e.g. G{\"u}ler et al 1993\\ $\bullet$ structured singular values, e.g. S. Wright; M. Wright\\ But ${\rm cond}(F_n) \rightarrow \infty$. \end{slide} \begin{slide}{The right-hand side} \[ \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 \\ {\cyan -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}{Ill-Posed System} {\cyan \bf Proposition } The {\cyan 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.\\ \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}{Condition Number Growth} {\cyan observe}\\ $\bullet$ condition number of $F_n^T F_n$ is greater than the largest eigenvalue of the block $AZ^{-2}X^2A^T$;\\ $\bullet$ equivalently, $\frac 1{{\rm cond} (F_n^T F_n)}$ is smaller than the reciprocal of this largest eigenvalue.\\ $\bullet$ If $x,z$ stay in neighbourhood of central path, then $\min_i(z_i/x_i)$ is $O(\mu)$. \begin{quote} {\yellow THEN:}\\ {\cyan reciprocal of the condition number of $F_n$ is $O(\mu)$.} \end{quote} {\yellow Summary:}\\ {\cyan precondition} initial ill-conditioned optimality conditions from log-barrier problem\\ {\cyan Block Elimination} brings back ill-conditioning \end{slide} \begin{slide}{Ex. Catastrophic Roundoff Error} $ 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: $ {\cyan \|r_b\| = 0.081699, \quad \|R_d\| = 0.36537}, \quad {\yellow \mu=x^Tz/n = 2.2528e-008} $ 5 decimals rounding before/after arithmetic\\ centering with $\sigma = .1$\\ {\bf BUT:} residuals are NOT order $\mu$. \end{slide} \begin{slide}{Search Direction} search direction is found using: \\ \qquad\qquad (i) {\cyan full matrix} $F^{\prime}_{\mu}$; \qquad (ii) {\cyan backsolve matrix} $F_n$ \[ \pmatrix{\Delta x\cr \Delta y\cr \Delta z}= {\cyan \pmatrix{ \fbox{8.17000e-02} \cr -1.35440e-08 \cr 1.63400e-01 \cr -2.14340e-08 \cr 1.63400e-01 } }; \quad = {\cyan \pmatrix{ \fbox{-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}{Alternate Second Step;\\ Stable Reduction} {\bf Assuming!} $A=[I_m ~ E]$. Partition \\ $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_s{\yellow K}=\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{ \begin{array}{c} 0 \cr I_m \end{array} & \begin{array}{lr} 0~~ & ~~~A^T \cr E~~ & ~~0 \end{array} & \begin{array}{c} I_n \cr 0 \end{array} \cr \begin{array}{c} 0 \cr 0 \end{array} & \fbox{ \begin{minipage}{1.2in} {\cyan $\begin{array}{cc} -Z_mE & -X_m \cr Z_v & -X_vE^T \end{array}$ } \end{minipage} } & \begin{array}{c} 0 \cr 0 \end{array} }. \end{array} \] \end{slide} \begin{slide}{The right-hand side} \[ \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\\ {\cyan -Z_mr_p -X_m(R_d)_m +Z_mX_me -\mu e}\\ {\cyan -X_v(R_d)_v +Z_vX_ve -\mu e} \end{array} \right) \end{array} \] \end{slide} \begin{slide}{Equivalent View of Stable Linearization} find (hopefully sparse) representation\\ {\yellow range of $N$ is nullspace of $A$} \[ Ax=b \quad \mbox{ if and only if } \quad x=\hat{x}+Nv, \mbox{ for some } v \in \Re^{n-m}. \] e.g. symmetric form \[ Ex_v \leq b, x_v \geq 0, \quad E \in \Re^{m \times (n-m)}\qquad {\cyan \Leftrightarrow }\qquad x_m+Ex_v=b, x\geq 0 \] Therefore assume $E$ sparse and \[ A=\pmatrix{ I_{m}& E}, \quad N = \pmatrix{ -E \cr I_{n-m} }. \] \end{slide} \begin{slide}{Substitute for $z,x$; Eliminate} {\cyan \bf Theorem} The primal-dual variables $x,y,z$, with $x=\hat{x}+Nv \geq 0,~z=c-A^Ty \geq 0$, are optimal for (LP),(DLP) if and only if they satisfy the single bilinear optimality equation \[ {\cyan F(v,y):=\Diag(c-A^Ty) \Diag(\hat{x}+Nv)e = 0}. \] \epr A single (perturbed) optimality conditions to use for the primal-dual method, \[ {\cyan F_\mu(v,y):= \Diag(c-A^Ty) \Diag(\hat{x}+Nv)e -\mu e = 0}. \] \end{slide} \begin{slide}{Linearization for Search Direction, $\Delta s$} \[ {\cyan - F_{\mu}(v,y) = F^{\prime}_{\mu}(v,y) \Delta s} \quad \Delta s:=\pmatrix{\Delta v \cr \Delta y} \] Jacobian matrix is \[ F^{\prime}_{\mu}(v,y) = \pmatrix{ \Diag(c-A^Ty)N \quad - \Diag(\hat{x} +Nv)A^T } \] system to solve for search direction is \[ -F_{\mu}(v,y) = \Diag(c-A^Ty) N {\cyan \Delta v} - \Diag(\hat{x}+Nv) A^T{\cyan \Delta y} . \] first part usually large, $n-m$ variables second part usually small, only $m$ variables. \end{slide} \begin{slide}{Well Conditioned} {\cyan \bf Theorem} \label{thm:cond} Consider the primal-dual pair (LP),(DLP). Suppose that $A$ is onto (full rank), the range of $N$ is the null space of $A$, $N$ is full column rank, and $(x,y,z)$ is the {\em unique} primal-dual optimal solution. Then the matrix of the linear system \beq \label{eq:SDPoptcondlinearize} \begin{array}{rcl} - F_{\mu} &=& F^{\prime}_{\mu} \Delta s \\ &=& Z N\Delta v - X A^T\Delta y \end{array} \eeq ($F_\mu^\prime$ is Jacobian of $F_\mu$) is {\yellow nonsingular}. \end{slide} \begin{slide}{Proof of Theorem} \bpr Suppose that ${F}'_{\mu}(v,y)\Delta s = 0$. We need to show that $\Delta s= (\Delta v, \Delta y) = 0$. Let ${\cal B}$ and ${\NN}$ denote the set of indices $j$ such that $x_j=\hat{x}_j + (Nv)_j > 0$ and set of indices $i$ such that $z_i=c_i - (A^Ty)_i > 0$, respectively. Under the nondegeneracy (uniqueness) and full rank assumptions, we get ${\cal B} \bigcup {\NN} = \{1, ... n\}$, ${\cal B} \bigcap {\NN} = \emptyset$, and the cardinalities $|{\cal B}| = m$, $|{\NN}| = n-m$. Moreover, the submatrix $A_{\BB}$, formed from the columns of $A$ with indices in $\BB$, is nonsingular. By our assumption and the linearization definition, we get that \[ \left({F}'_{\mu}(v,y) \Delta s \right)_k = (c-A^Ty)_k (N\Delta v)_k - (\hat{x} + Nv)_k (A^T \Delta y)_k = 0, \quad \forall k. \] \end{slide} \begin{slide}{Proof of Theorem cont...} From the definitions of $\BB,\NN$, this implies that \beq \label{eq:BNzero} (A^T \Delta y)_j = 0, \forall j \in {\BB}, \qquad (N\Delta v)_i = 0, \forall i \in {\NN}. \eeq The left part of \eqref{eq:BNzero} implies $A^T_{\cal B} \Delta y = 0$, i.e. we obtain $\Delta y = 0$. It remains to show that $\Delta v = 0$. From the definition of $N$ we have $AN = 0$. Therefore, using the right part of \eqref{eq:BNzero} implies \[ \begin{array}{rcl} 0 &=& \pmatrix{A_\BB & A_\NN}\pmatrix{ (N\Delta v)_\BB \cr (N\Delta v)_\NN} \\&=& A_\BB (N\Delta v)_\BB + A_\NN (N\Delta v)_\NN \\&=& A_\BB (N\Delta v)_\BB. \end{array} \] \end{slide} \begin{slide}{Proof of Theorem cont...} By the right part of \eqref{eq:BNzero} and the nonsingularity of $A_\BB$, we get \[ N\Delta v = 0. \] Now, full rank of $N$ implies $\Delta v=0$. \epr \end{slide} \begin{slide}{Primal-Dual Algorithm} {\cyan follow usual primal-dual interior-point framework: Newton's method applied to perturbed system of optimality conditions with damped step lengths for maintaining nonnegativity constraints. } \vspace{.5in} {\cyan Differences}: \\$\bullet$ eliminate, primal-dual linear feasibility (exact feas. maintained) \\$\bullet$ search direction found using PCG (LSQR) \\$\bullet$ no backtracking to preserve sufficient positivity of $z,x$ \\$\bullet$ crossover step (affine scaling, perturbation parameter $\mu = 0$ and full Newton step close to optimum; exterior method; NLF - Newton Liberation Front) \\$\bullet$ identify zero values for $x$, {\em purification step}. \end{slide} \begin{slide}{Preconditioning Techniques} \[ Z:= Z(y)=\Diag(c-A^Ty),\quad X:= X(v)=\Diag(\hat{x} + Nv) \] \[ J:=F_{\mu}^\prime(v,y) = \pmatrix{ ZN & -XA^T} \quad \mbox{Jacobian} \] find a preconditioner {\em simple nonsingular} $M$ such that $JM^{-1}$ is well conditioned and solve better conditioned systems {\yellow $JM^{-1} \Delta q = -F_{\mu}$ and $M \Delta s = \Delta q$}. \[ \mbox{We look for:} \quad {\cyan M^TM \cong J^TJ} \] LSQR (Paige-Saunders) implicitly solves normal equations \[ {\cyan J^TJ ds = -J^T F^{\prime}_{\mu}} \] \end{slide} \begin{slide}{Optimal Diagonal Column Preconditioning} simplest of preconditioners; given square matrix $K$ \[ {\cyan \omega (K) = \frac{ \trace(K)/n}{\det(K)^{1/n}}} \quad \mbox{condition number} \] If $M = \arg \min \omega((JD)^T(JD))$ over all positive diagonal matrices $D$ then (Dennis-W. 1990) \[ {\cyan M_{ii} = 1/\|J_{:i}\|} \quad \mbox{$i$-th column norm} \] \end{slide} \begin{slide}{Partial (Block) Cholesky Preconditioner} \[ J^TJ= \pmatrix{ N^TZ^2N & -N^TZXA^T \cr -AXZN & AX^2A^T }. \] For $z,x$ near central path, i.e. $ZX \cong \mu I$, off diagonal terms $\cong 0$ block (partial) Cholesky preconditioning is good preconditioner ~\\ $Q$-less QR factorization $Q_ZR_Z=ZN,~Q_XR_X=XA^T$ \[ R_Z^T R_Z=N^TZ^2N, \quad R_X^TR_X=AX^2A^T. \] \[ J^TJ \cong M^TM, \quad {\cyan M=\pmatrix{R_Z & 0 \cr 0 & R_X}} \] expensive! \end{slide} \begin{slide}{Crossover Criteria/Quadratic Convergence} assume nonsingular Jacobian at optimality (so unique primal and dual solutions, $s^*$) standard theory for Newton's method:\\ {\cyan $\exists$ quadratic convergence neighbourhood of $s^*$ } {\cyan \bf Theorem} (Kantorovich) Let $r>0$, $s_0 \in \mathbb R^n$, $F: \mathbb R^n \rightarrow R^n$, and assume that $F$ is continuously differentiable in $\NN(s_0, r)$. Assume for a vector norm and the induced operator norm that $J \in \mbox{Lip}_\gamma(\NN(s_0, r))$ with $J(s_0)$ nonsingular, and that there exist constants $\beta$, $\eta\geq 0$ such that $$ \|J(s_0)^{-1} \| \leq \beta, \;\; \|J(s_0)^{-1} F(s_0)\| \leq \eta . $$ Define $\alpha = \beta \gamma \eta $. If $\alpha \leq {1\over 2}$ and $r\geq r_0 :=(1- \sqrt{1-2\alpha}\, ) / (\beta\gamma)$. \end{slide} \begin{slide}{Quadratic Convergence cont...} Then the sequence $\{s_k\}$ produced by \[ s_{k+1} = s_k - J(s_k)^{-1} F(s_k), \;\; k = 0, 1,\dots, \quad \mbox{\cyan NLF} \] is well defined and converges to $s_*$, a unique zero of $F$ in the closure of $\NN(s_0, r_0)$. If $\alpha < {1\over 2} $, then $s_*$ is the unique zero of $F$ in $\NN(s_0, r_1)$, where $r_1:= \min [ r, (1+ \sqrt{1-2\alpha }\,)/ (\beta \gamma) ] $ and $$ \| s_k- s_* \| \leq (2\alpha)^{2^k} {\eta \over \alpha}, \;\; k = 0, 1, \dots, $$ \epr \end{slide} \begin{slide}{Lipschitz Constant for Region of Convergence} {\cyan \bf Lemma} The Jacobian $$F'(v,y){\cyan \pmatrix{ \cdot \cr \cdot }} := \pmatrix{\Diag(c-A^Ty)N{\cyan \cdot } & \Diag(\hat x + Nv)A^T{\cyan \cdot } }$$ is Lipschitz continuous with constant \[ \gamma=\sqrt 2 \|A\|\|N\| \] with respect to $(v,y)$ \epr \end{slide} \begin{slide}{Proof} \bpr We let $\Delta s=\pmatrix{\Delta v \cr \Delta y}$. Since \begin{eqnarray*} \|F'(s) - F'(\bar s)\| &=& \max{ \|(F^\prime (s)-F^\prime (\bar{s})) \Delta s\| \over \|\Delta s\| }\\ &=& \max{ \| \Diag(A^T(y-\bar{y})) N \Delta v - \Diag(A^T \Delta y) N (v -\bar{v}) \| \over \|\Delta s\| }\\ &\leq& \max{\|A^T(y-\bar{y})\|\| N \Delta v\| + \| A^T \Delta y\| \| N (v - \bar{v}) \| \over \| \Delta s\|}\\ &\leq& \|A\|\|N\| \|y-\bar y\| + \|A\|\|N\|\|v - \bar v\|\\ &\leq & \sqrt 2 \|A\|\|N\|\|s-\bar s\|. \end{eqnarray*} Therefore a Lipschitz constant is $ \gamma = \sqrt 2 \|A\| \| N \|$. \epr \end{slide} \begin{slide}{Region of Quadratic Convergence} $\|[ZN\;-XA^T]^{-1}\| \leq \beta$, e.g. using smallest singular value \[ \|J^{-1} F_0(v,y)\| = \| [ZN \; -XA^T]^{-1} (-XZe)\| \leq \eta. \] {\cyan \bf Theorem} Suppose that \[ \alpha = \gamma \beta \eta < \frac 12. \] Then the sequence $s_k$ generated by $$s_{k+1} = s_k - J(s_k)^{-1} F(s_k)$$ {\yellow converges (quadratically)} to $s^*$, the unique zero of $F$ in the neighbourhood $\NN(s_0, r_1)$. \epr \end{slide} \begin{slide}{Purify Step} $\bullet$ detect zero variables/active constraints at optimality $\bullet$ use the Tapia indicators 1995, \[ \frac {(x_{k+1})_i}{(x_{k})_i} \quad \mbox{ratio of $i$-th component of iterates} \] $\bullet$ perform a pivot step to eliminate these variables (small indicators) \end{slide} \begin{slide}{Summary; Path-following and NOT Interior-point} $\bullet$ staying interior is a heuristic for staying within a neighbourhood of the central path\\ $\bullet$ staying interior (well-centered) is required for numerical accuracy when solving the {\em current} ill-conditioned reduced systems\\ $\bullet$ {\cyan Main Advantages} \begin{description} \item no loss of sparsity \item high accuracy solutions available if desired \item exact primal and dual feasibility at each iteration (true duality gap) \item warm starts \item fast convergence (no backtracking from boundary) \end{description} \end{slide} \begin{slide}{Numerical Tests} $\bullet$ {\cyan randomly} generated data with {\cyan known optimum}\\ $\bullet$ {\cyan well-conditioned} basis matrix \\ $\bullet$ stopping condition relative gap $10^{-12}$\\ $\bullet$ MATLAB 6.5, Pentium 3 733MHz, 256MB RAM\\ $\bullet$ iterative approach; {\cyan LSQR} (Paige-Saunders); different preconditioners $\bullet$ NEQ {\cyan stalls} with relative gap approximately $10^{-11}$ on many problems \end{slide} \begin{slide}{NEQ vs Stable Method-Direct Solver} \tiny %\begin{table*} \begin{center} \begin{tabular}{|c|c|c|r|r|r|c|c|c|c|}\hline data& $m$& $n $ & ${\rm nnz}(E)$ & cond($A_\BB$) & cond(J) &\multicolumn{2}{c}{NEQ} \vline & \multicolumn{2}{c}{Stable direct} \vline \\ \hline & & & & & & D\_time & its & D\_Time & its \\ \hline 1&$100$&$200$ & $1233$ & $51295$ & $32584$ & $0.03$ & $*$ & $0.06$ & $6$ \\ \hline 2&$200$&$400$ & $2526$ & $354937$ & $268805$ & $0.09$ & $6$ & $0.49$ & $6$ \\ \hline 3&$200$&$400$ & $4358$ & $63955$ & $185503$ & $0.10$ & $*$ & $0.58$ & $6$ \\ \hline 4&$400$&$800$ & $5121$ & $14261771 $ & $2864905$ & $0.61$ & $ *$ & $ 3.66$ & $6 $ \\ \hline 5&$400$&$800$ & $8939$ & $459727 $ & $256269$ & $0.64 $ & $6$ & $4.43 $ & $6$ \\ \hline 6&$800$&$1600$ & $10332$ & $11311945 $ & $5730600$ & $5.02 $ & $6$ & $26.43 $ & $6$ \\ \hline 7&$800$&$1600$ & $18135$ & $4751747$ & $1608389$ & $5.11 $ & $*$ & $33.10$ & $6$ \\ \hline \end{tabular} \end{center} %\end{table*} \normalsize \tiny{ ${\rm nnz}(E)$ - number of nonzeros in $E$; \\ cond($\cdot$) - condition number; $J=(ZN~-XA^T)$ at optimum;\\ D\_time - avg. for search direction per iter.;\\ its - for interior point method\\ {\cyan *} denotes NEQ stalls at $10^{-11}$ } \end{slide} \begin{slide}{Stable Method with LSQR\\ and Two Precond.} \tiny{ \begin{center} \begin{tabular}{|c|c|c|c|c|c|c|c|c|} \hline data set& \multicolumn{4}{c}{LSQR with ILU} \vline & \multicolumn{4}{c}{LSQR with Diag} \vline \\ \hline & D\_Time & its & L\_its & Pre\_time & D\_Time & its & L\_its & Pre\_time \\ \hline 1& $0.15$ & $6$ & $37$ & $0.06$ & $0.41$ & $6$ & $556$ & $0.01$ \\ \hline 2& $3.42$ & $6$ & $343$ & $0.28$ & $2.24$ & $6$ & $1569$ & $0.00$ \\ \hline 3& $2.11$ & $6$ & $164$ & $0.32$ & $ 3.18$ & $ 6$ & $1595$ & $0.00$ \\ \hline 4& NA & Stalling & NA & NA & $13.37$ & $6$ & $4576$ & $0.01$ \\ \hline 5& NA & Stalling & NA & NA & $21.58 $ & $6$ & $4207 $ & $ 0.01$ \\ \hline 6& NA & Stalling & NA & NA & $90.24 $ & $6$ & $9239$ & $ 0.02$ \\ \hline 7& NA & Stalling & NA & NA & $128.67 $ & $6$ & $8254$ & $ 0.02$ \\ \hline \end{tabular} \end{center} Same data sets as above;\\ two different preconditioners (diagonal and incomplete Cholesky with drop tolerance $0.001$);\\ D\_time - average time for search direction;\\ its - iteration number of interior point methods;\\ L\_its - average number LSQR iterations per major iteration;\\ Pre\_time - average time for preconditioner;\\ Stalling - LSQR cannot converge due to poor preconditioning. } \end{slide} \begin{slide}{LSQR with Block Cholesky preconditioner} \tiny{ \begin{center} \begin{tabular}{|c|c|c|c|c|} \hline data set& \multicolumn{4}{c}{LSQR with block Chol. Precond.} \vline \\ \hline & D\_Time & its & L\_its & Pre\_time \\ \hline 1& $0.09$ & $6$ & $4$ & $0.07$ \\ \hline 2& $0.57 $ & $6$ & $5$ & $0.48$ \\ \hline 3& $0.68$ & $6 $ & $5$ & $0.58$ \\ \hline 4& $5.55$ & $ 6$ & $6$ & $5.16 $ \\ \hline 5& $6.87 $ & $6 $ & $6 $ & $6.45 $ \\ \hline 6& $43.28 $ & $6 $ & $5$ & $41.85 $ \\ \hline 7& $54.80 $ & $6 $ & $5$ & $53.35$ \\ \hline \end{tabular} \end{center} } \end{slide} \begin{slide}{Iterations/Degeneracy} %\begin{figure}[htb] \epsfxsize=250pt \centerline{\epsfbox{deg.eps}} %\caption{Iterations for Degenerate Problem} %\end{figure} \end{slide} \begin{slide}{Sparse; Well conditioned $A_\BB$} $\bullet$ about {\cyan 3-4 nonzeros per row} in $E$ \\ $\bullet$ the Jacobian nonsingular at optimum. \\ $\bullet$ {\cyan well-conditioned basis matrix}, $A_\BB$\\ $\bullet$ the same dimensions and two dense columns, while total number of nonzeros increases\\ \vspace{.1in} $\bullet$ The loss in sparsity has essentially no effect on NEQ, since the $ADA^T$ matrix is dense due to the two dense columns. But we can see the negative effect that the loss of sparsity has on the stable direct solver.\\ However, we see that for these problem instances, using LSQR with the stable system can be up to twenty times faster then NEQ solver. \end{slide} \begin{slide}{Sparse; Well conditioned $A_\BB$, cont...} \begin{tiny} \begin{tabular}{|cccr|rc|cc|ccr|} \hline \multicolumn{4}{|c|}{data sets}& \multicolumn{2}{|c|}{NEQ} & \multicolumn{2}{|c|}{Stable Direct} & \multicolumn{3}{|c|}{LSQR} \\ \hline % using \hline to add a horizon line below a line. Name & cond($A_\BB$) & cond(J) & nnz(E) & D\_Time & its & D\_Time & its & D\_Time & its &L\_its\\ \hline nnz2 & 19 & 13558 & 4490 & 9.87 & 7 & 18.78 & 7 & 0.55 & 7 & 81\\ nnz4 & 21 & 19540 & 6481 & 10.09 & 7 & 20.74 & 7 & 0.86 & 7 & 106\\ nnz8 & 28 & 10170 & 10456 & 10.05 & 7 & 29.48 & 7 & 1.51 & 7 & 132\\ nnz16 & 76 & 11064 & 18346 & 10.08 & 7 & 35.48 & 7 & 3.65 & 7 & 210\\ nnz32 & 201 & 11778 & 33883 & 10.04 & 9 & 41.49 & 9 & 8.96 & 8 & 339 \\ \hline \end{tabular} ${\rm cond(\cdot)}$ - (rounded) condition number; \\ nnz($E$) - number of nonzeros in $E$;\\ D\_time - average time for search direction;\\ its - number of iterations;\\ L\_its - average number LSQR iterations per major iteration;\\ All data sets have the same dimension, $1000 \times 2000$, and have 2 dense columns. \end{tiny} \end{slide} \begin{slide}{Sparse; Well conditioned $A_\BB$, ... Size} \begin{tiny} The time for the {\cyan NEQ} solver is proportional to {\cyan $m^3$}. The stable direct solver is about twice that of NEQ. LSQR is the best among these 3 solvers on these instances. The computational advantage of LSQR becomes more apparent as the dimension grows. \begin{table} \begin{tabular}{|cccr|rc|rc|rc|} \hline \multicolumn{4}{|c|}{data sets}& \multicolumn{2}{|c|}{NEQ} & \multicolumn{2}{|c|}{Stable Direct} & \multicolumn{2}{|c|}{LSQR} \\ \hline % using \hline to add a horizon line below a line. name&size&cond($A_\BB$)&cond(J)& D\_Time & its & D\_Time & its & D\_Time & its \\ \hline sz1 & $400\times800$ & 20 & 2962 & 0.63 & 7 & 1.37 & 7 & 0.15 & 7 \\ sz2 & $400\times1600$ & 15 & 2986 & 0.63 & 7 & 1.36 & 7 & 0.26 & 7 \\ sz3 & $400\times3200$ & 13 & 2358 & 0.63 & 7 & 1.39 & 7 & 0.53 & 7 \\ sz4 & $800\times1600$ & 19 & 12344 & 5.08 & 7 & 9.60 & 7 & 0.32 & 7 \\ sz5 & $800\times3200$ & 15 & 15476 & 5.06 & 7 & 9.64 & 7 & 0.76 & 7 \\ sz6 & $1600\times3200$ & 20 & 53244 & 39.01 & 7 & 72.12 & 7 & 1.35 & 7 \\ sz7 & $1600\times6400$ & 16 & 56812 & 38.83 & 7 & 72.16 & 7 & 2.32 & 8 \\ sz8 & $3200\times6400$ & 19 & 218664 & 346.24 & 7 & 549.44 & 7 & 2.99 & 7 \\ \hline \end{tabular} \tiny{ ${\rm cond(\cdot)}$ - (rounded) condition number;\\ D\_time - average time for search direction;\\ its - number of iterations\\ } \end{table} \end{tiny} \end{slide} \begin{slide}{Sparse; Well conditioned $A_\BB$, ...\\ \# Dense Cols} \begin{tiny} \begin{table} \begin{tabular}{|cccc|rc|rc|rc|} \hline \multicolumn{4}{|c|}{data sets}& \multicolumn{2}{|c|}{NEQ} & \multicolumn{2}{|c|}{Stable Direct} & \multicolumn{2}{|c|}{LSQR} \\ \hline % using \hline to add a horizon line below a line. name& dense cols.& cond($A_\BB$)& cond(J) & D\_Time & its & D\_Time & its & D\_Time & its \\ \hline den0 & 0 & 18 & 45 & 1.21 & 6 & 2.96 & 6 & 0.41 & 6 \\ den1 & 1 & 19 & 13341 & 10.11 & 7 & 18.29 & 7 & 0.47 & 7 \\ den2 & 2 & 19 & 18417 & 9.98 & 7 & 19.35 & 7 & 0.60 & 7 \\ den3 & 3 & 19 & 19178 & 9.92 & 7 & 18.64 & 7 & 0.72 & 7 \\ den4 & 4 & 18 & 18513 & 9.89 & 7 & 18.72 & 7 & 0.97 & 7 \\ \hline \end{tabular} \tiny{ ${\rm cond(\cdot)}$ - (rounded) condition number; \\ D\_time - average time for search direction;\\ its - number of iterations.\\ } \end{table} \end{tiny} \end{slide} \begin{slide}{ LSQR iterations at different stages} \begin{figure}[htb] \epsfxsize=250pt \centerline{\epsfbox{lsqrits.eps}} \label{fig:lsqrits} \end{figure} \end{slide} \begin{slide}{No Backtracking - Complete Step to Boundary} \begin{figure}[htb] \epsfxsize=250pt \centerline{\epsfbox{backtrack.eps}} \end{figure} \end{slide} \begin{slide}{Warm Starts} {\yellow goal}: start from optimal solution as an initial starting point to a perturbed problem. \vspace{.2in} {\cyan stable method} particularly successful at performing warm starts for small perturbations, since the Jacobian (under nondegeneracy assumptions) is nonsingular \vspace{.2in} optimality conditions become \[ F(x,y,z):= \left( \begin{array}{cl} (A+\Delta A)^T y +z - c \\ (A+\Delta A) x -b \\ ZXe \end{array} \right) = \pmatrix{\Delta c \cr \Delta b \cr 0}. \] \end{slide} \begin{slide}{} {\cyan \bf Theorem} Suppose that the linear programming problem is nondegenerate and $x^*, (y^*,z^*)$ is a strictly complementary primal-dual optimal solution. Let $\BB$ and $\NN$ be the partition of the index set of variables: $$ \BB = \{ i: x_i^*>0\} ~\mbox{and }~ \NN = \{i: 1\leq i \leq n, i \not\in \BB\}. $$ Suppose we perturb $A$ to $A + \Delta A$, $b$ to $b+\Delta b$, and $c$ to $c+\Delta c$. Assume $(A+\Delta A)_\BB$ is nonsingular. \begin{enumerate} \item \label{item:starta} Then, a full step in the affine direction from the starting point $x^*, y^*, z^*$ {\yellow converges in one step} to a solution satisfying the optimality conditions. \item Furthermore, if $\Delta c, \Delta b, \Delta A$ are sufficiently small, then we get both $x\geq 0,z\geq 0$. \end{enumerate} \epr \end{slide} \begin{slide}{} \bpr The affine search direction $\Delta x, \Delta y, \Delta z$ with the starting point $x^*, y^*, z^*$ is the solution to the following system \begin{eqnarray*} (A + \Delta A) \Delta x &=& \Delta b - \Delta A x^*\\ (A+ \Delta A)^T \Delta y + \Delta z &=& \Delta c- \Delta A^T y^*\\ x^* \circ \Delta z + z^* \circ \Delta x &=& 0 . \end{eqnarray*} By noting that $x^*_\NN =0$ and $ z^*_\BB =0$, solving the above system yields \begin{eqnarray*} &\Delta x_\BB=(A+ \Delta A)_\BB^{-1}(\Delta b - \Delta A x^*), ~~ \Delta x_\NN = 0,&\\ &\Delta y = (A + \Delta A)_\BB^{-T} (\Delta c_\BB - \Delta A_\BB^T y^*),&\\ & \Delta z_\BB =0, ~~ \Delta z_\NN = \Delta c_\NN - \Delta A^T_\NN y^* - (A + \Delta A)_\NN^T (A+ \Delta A)_\BB^{-T} (\Delta c_\BB- \Delta A^T_\BB y^*).& \end{eqnarray*} \end{slide} \begin{slide}{} One can now verify that $x^* +\Delta x$, $y^* + \Delta y$, and $z^* + \Delta z$ is a solution to system equation {optper}. When $\Delta A$, $\Delta b$, and $\Delta c$ are sufficiently small, then $\Delta x_\BB$ and $\Delta z_\NN$ are small too, and thus $x^* + \Delta x \geq 0$ and $z^* + \Delta z \geq 0$. \epr \end{slide} \begin{slide}{Warm Start Convergence Radii for $A+r\Delta A$, $b+r\Delta b$, $c+r\Delta c$} \begin{tiny} \begin{table} \begin{center} \begin{tabular}{|c|c|c|}\hline Problem 1 with & Problem 2 with & Problem 3 with \\ $\|A\| = 10.5$, $\|b\|=465.6$, & $\|A\| = 15.1$, $\|b\|=582.0$, & $\|A\| = 15.5$, $\|b\|=755.8$, \\ $\|c\|=155.6$ & $\|c\|=217.7$ & $\|c\|=215.9$ \\ \hline $0.05$ & $0.01$ & $0.01 $ \\ \hline $0.05$ & $0.03$ & $0.01 $ \\ \hline $0.09$ & $0.03$ & $0.01 $ \\ \hline $0.05$ & $0.03$ & $0.01 $ \\ \hline $0.13$ & $0.01$ & $0.11 $ \\ \hline $0.09$ & $0.01$ & $0.07 $ \\ \hline $0.09$ & $0.01$ & $0.01 $ \\ \hline $0.07$ & $0.03$ & $0.01 $ \\ \hline $0.03$ & $0.02$ & $0.01 $ \\ \hline $0.11$ & $0.07$ & $0.01 $ \\ \hline \end{tabular} \end{center} \end{table} \end{tiny} \end{slide} \begin{slide}{} \section*{(Nearest) Euclidean Distance Matrix Completion using SDP} {\bf Given}:\\ {\em pre-distance matrix} $A\in\Sn$ (nonnegative with zero diagonal)\\ {\em weight matrix} $H\in \Sn$: \[ \mbox{(NEDM)}\quad \mu^*= \min \frac12 \normF{H\circ (A-D)}^2 ~ \mbox{subject to: } D \in \EDM \] \begin{tiny} $\EDM = \{ D=(d_{ij})\in \Sn : d_{ij} = \|x_i-x_j\|^2, \mbox{ for some } x_i \in \Re^k \}$, $k$ is {\em embedding dimension} \end{tiny} \end{slide} \begin{slide}{} Applications e.g. molecular conformation problems in chemistry; multidimensional scaling and multivariate analysis problems in statistics; genetics, geography, .... \end{slide} \begin{slide}{} \subsection*{Mixed-Cone Formulation} direct approach using a mixed SDP and second-order (or Lorentz) cone problem: \[ \begin{array}{rcl} &\min &\alpha \\ &\mbox{s.t.}& Y=H\circ(\LL(X)-A),~ \normF{Y} \le \alpha \\ && X \in \Stn, Y \in \Sn, X \in \SDP \end{array} \] where $X \in SDP \Rightarrow \LL(X) \in \EDM$\\ (Public domain software packages are available) \end{slide} \begin{slide}{} \[ B=[x_1~x_2~\ldots ~ x_n], \quad k \times n \] \[ D_{ij}=\|x_i-x_j\|^2 = -2x_i^Tx_j+ \|x_i\|^2 + \|x_j\|^2 \] \[ D= -2B^TB+ e \left(\diag (B^TB)\right)^T+ \left(\diag (B^TB)\right) e^T \] ~~\\ With $X=B^TB \succeq 0$\\ connection between \SDP and \EDM. \end{slide} \begin{slide}{} \subsubsection*{Operator Notation:} $\usvec, \usMat, \svec, \sMat $ \[x=\svec(X) \in \R^{n+1 \choose 2}, \quad X= \sMat(x) \] $\sqrt{2}$ times vector (columnwise) from upper-triang of $X$. ${n+1\choose 2}=n(n+1)/2$; $\sqrt{2}$ guarantees isometry. $\sMat:=\svec^{-1}$ mapping into $\Sn$\\ adjoint transformation $\sMat^*= \svec$: \begin{eqnarray*} \left<\sMat(v),S\right> &=& \trace \sMat(v)S \\ &=& v^T \svec(S) = \left \end{eqnarray*} \end{slide} \begin{slide}{} \begin{center} $D$ is \EDM \quad ($\subset \Sn$)\\ ~\\ \fbox{ {\em iff} } ~\\ \end{center} \[ D=\LL( X):=\pmatrix{0 & \diag ( X)^T \cr \diag ( X) & \diag(X)e^T+e\diag(X)^T -2X }, \] \begin{center} for some $X \succeq 0, X \in {\mathcal S}^{n-1}$ \end{center} ($e$ is vector of ones) \[ \LL: \Stn \rightarrow \Sn, \quad \LL({\mathcal S}^{n-1}_+) = \EDM \] \end{slide} \begin{slide}{} with partition: \[ D=\left[\matrix{\alpha & d^T\cr d& {\bar D} \cr}\right], \] where $\alpha \in \R$ \[ \LL^*(D)=2\left(\Diag(d)+\Diag(\bar De) -\bar D\right) \] \[ \LL^\dagger (D)= \frac 12 \left( de^T+ed^T - \bar D \right) \] \[ \LL^*,\LL^\dagger: \Sn \rightarrow \Stn, \quad \LL^\dagger (EDM)={\mathcal S}^{n-1}_+ \] \end{slide} \begin{slide}{} \subsection*{Duality and Optimality Conditions} (using $X = \sMat(x)+I$) an equivalent problem is: \[ \mu^* := \min \frac 12 \|H \circ (A-\LL(X)) \|_F^2 \quad\mbox{subject to \quad $X \succeq 0$} \] strong (Lagrangian) duality holds (Slater's holds for primal and holds for dual if the graph is complete) \[ \mu^*=\nu^*:= \max_{\Lambda \succeq 0} \min_{X} \frac 12 \|H \circ (A-\LL(X)) \|_F^2 - \trace \Lambda X \] \end{slide} \begin{slide}{} change to {\bf Wolfe dual} and obtain optimality conditions: \[ \begin{array}{rcll} X&:=&\sMat(x) \succeq 0 & \mbox{(primal feasibility)}\\ \Lambda &:=& \LL^* \left\{ H^{(2)} \circ( \LL(X) ) \right\}-C , ~~\Lambda \succeq 0 & \mbox{(dual feasibility)}\\ \Lambda X &:=& 0 & \mbox{(complementary slack.)}\\ \end{array} \] \end{slide} \begin{slide}{} eliminate $\Lambda$\\ exact primal-dual feasibility during iterations\\ full rank Jacobian at optimality. ~~\\ {\em single bilinear (perturbed) equation} in $x$; $\quad F_{\mu}(x): \R^{n \choose 2} \rightarrow \Mmn$ \fbox{ \begin{minipage}{4.05in} \[ \begin{array}{c} \vspace{.5mm} F_{\mu}(x):= \left[ \LL^* \left\{ H^{(2)} \circ( \LL(X) ) \right\}-C \right]X - \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 - with PCG \end{slide} \begin{slide}{} Let $\WW(x):= \LL^* \left\{ H^{(2)} \circ( \LL(x) ) \right\}$ ~\\ Linearization for search direction $\Delta x$ at current $x=\svec (X)$: \fbox{ \begin{minipage}{4.05in} \[ \begin{array}{lcl} F^{\prime}_{\mu}(x) \Delta x = \left[ \WW(x) -C \right] \Delta x + \left[ \WW( \Delta x) \right]X \end{array} \] \end{minipage} } ~\\ This is a linear, full rank, overdetermined system.\\ Our search direction $\Delta x$ is its (approx.) 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 pre-distance $n\times n$ matrix $A$ \\ {\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 $X^0,\Lambda^0:=\left(\WW (X) -C \right)~\succ~0$; $\mu>0$\\ {\bf $\bullet\bullet$ Set initial parameters:} \[ {\rm gap}=\trace \Lambda^0X^0;~~ \mu={\rm gap}/n;~~{\rm objval}=f(X^0);~~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}(x^k) \left( \Delta x^k \right) = -F_{\sigma \mu}(x^k), \] where $\sigma_k$ centering, $\mu_k=\frac{1}{n}\trace (\WW(X^{k})-C)X^k$ \[ X^{k+1} = X^k + \alpha_k \Delta X^k, ~~ \alpha_k > 0, \] so that both $X^{k+1},(\WW(X^{k+1})-C) \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: } $D=\LL(X) \in \EDM$ is approx. to $A$ \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) \end{slide} \begin{slide}{} \subsection*{Preconditioning} \[ \left(\Lambda+\XS \WW \right) P^{-1}(\widehat{\Delta x}) = - F_{\mu}(x), \] where \[ \widehat{\Delta x}=P(\Delta x) \] \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 \[ \min \omega ( (AD)^T(AD)) \quad \mbox{subject to: $D$ positive and diagonal} \] solution: $d_{ii}= 1/\normt{A_{:i}}, i=1,\ldots,n$ {\em explicit} expressions for preconditioner\\ inexpensive \end{slide} \begin{slide}{} {\em diagonal} operator $P$; evaluate using columns of $F_\mu^\prime(v)$. \\ $k\cong (i,j),~1\leq i < j\leq n$, strictly upper triangular part\\ \[ \begin{array}{rcl} \| (\Lambda+\XS\WW)(e_k) \|^2_F &=& \| \Lambda(e_k) \|^2_F + \| (\WW(e_k))X \|^2_F\\ && \qquad +\left<\Lambda(E_{ij}),(\WW(E_{ij}))X\right>, \end{array} \] where \[ \begin{array}{rcl} \Lambda (e_k) &=& \left\{ \begin{array}{cc} \frac 1{\sqrt{2}} \left( \Lambda_{:i}e_j^T+\Lambda_{:j}e_i^T \right), & \mbox{if } i < j\\ \left( \Lambda_{:i}e_i^T \right), & \mbox{if } i = j. \end{array} \right. \end{array} \] and $\XS \WW $ ..... inexpensive - 50\% reduction in LSQR iterations \end{slide} \begin{slide}{} \section*{Numerical Tests} Pentium 4; MATLAB 6.5; 1 GIG RAM. crossover heuristic: relative duality gap $< .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=220pt \centerline{\epsfbox{dens0050103.ps}} \caption{density .0005:.001:.003, CPU times and nnz($\Lambda$), n=200} \label{fig:densityn200} \end{figure} \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}