%\documentstyle[11pt]{article %\documentclass{article} \documentstyle{smjrnl} % fonts removed by Henry %\textheight22cm %\textheight16cm %\textheight18.5cm %\textwidth13cm %\usepackage{latexsym} \input epsf \input psfig \newtheorem{prop}{Proposition}[section] \newtheorem{lem}{Lemma}[section] \newtheorem{thm}{Theorem}[section] \newtheorem{cor}{Corollary}[section] \newtheorem{rem}{Remark}[section] \newtheorem{alg}{Algorithm}[section] \newcommand{\1}{\;\;\; } \newcommand{\diag}{\rm {diag} } \newcommand{\Diag}{\rm {Diag} } \newcommand{\Sn}{{\cal S}_n } \newcommand{\E}{{\cal E} } \newcommand{\Se}{{\cal S}_e } \newcommand{\Sd}{{\cal S}_d } \newcommand{\Sc}{{\cal S}_C } \newcommand{\Sh}{{\cal S}_H } \newcommand{\snn}{{\cal S}_{n-1} } \newcommand{\KK}{{\cal K} } \newcommand{\LL}{{\cal L} } \newcommand{\DD}{{\cal D} } \newcommand{\BB}{{\cal B} } \newcommand{\PP}{{\cal P} } \newcommand{\TT}{{\cal T} } \newcommand{\A}{{\cal A} } \newcommand{\p}{{\cal P} } \newcommand{\trace}{{\rm trace\,}} \newcommand{\kvec}{{\rm vec\,}} \newcommand{\bpr}{{\bf Proof.} \hspace{1 em}} \newcommand{\epr}{ \\ \hspace*{4.5in} $\Box$ } \newcommand{\beq}{ \begin{equation} } \newcommand{\eeq}{ \end{equation} } \newcommand{\bt}{ \begin{tabular} } \newcommand{\et}{ \end{tabular} } \begin{document} \begin{article} \bibliographystyle{plain} \journame{Small Journal Name} \volnumber{?} \issuenumber{?} \issuemonth{???} \volyear{199?} %% Do not delete either of the following two commands. %% Please supply facing curly brackets for the part you %% are not using for this article. %\received{May 1, 1991}\revised{} \authorrunninghead{A. Alfakih, A. Khandani, H. Wolkowicz} \titlerunninghead{Euclidean Distance Matrix Completion Problems} \setcounter{page}{1} %% This command is optional. %% May set page number only for first page in %% issue, if desired. %% <<== End of commands to be entered at Kluwers %% Authors, start here ==>> \title{SOLVING EUCLIDEAN DISTANCE MATRIX COMPLETION PROBLEMS VIA\\ SEMIDEFINITE PROGRAMMING} \authors{Abdo Y. Alfakih} \email{aalfakih@orion.math.uwaterloo.ca} \affil{ University of Waterloo, Department of Combinatorics and Optimization, Waterloo, Ontario N2L 3G1, Canada} \authors{Amir Khandani} \email{khandani@claude.uwaterloo.ca} \affil{ University of Waterloo, Department of Electrical \& Computer Engineering, Waterloo, Ontario N2L 3G1, Canada} \authors{Henry Wolkowicz} \email{henry@orion.math.uwaterloo.ca.} \affil{ University of Waterloo, Department of Combinatorics and Optimization, Waterloo, Ontario N2L 3G1, Canada. Research supported by Natural Sciences Engineering Research Council Canada. } \abstract{ Given a partial symmetric matrix $A$ with only certain elements specified, the Euclidean distance matrix completion problem (EDMCP) is to find the unspecified elements of $A$ that make $A$ a Euclidean distance matrix (EDM). In this paper, we follow the successful approach in \cite{JoKrWo:95} and solve the EDMCP by generalizing the completion problem to allow for approximate completions. In particular, we introduce a primal-dual interior-point algorithm that solves an equivalent (quadratic objective function) semidefinite programming problem (SDP). Numerical results are included which illustrate the efficiency and robustness of our approach. Our randomly generated problems consistently resulted in low dimensional solutions when no completion existed.} \keywords{Euclidean Distance Matrices, Semidefinite Programming, Completion Problems, Primal-Dual Interior-Point Method.} \noindent {\bf AMS Subject Classifications:} 90C25 90C27 15A48 \footnotetext{ This report is available by anonymous ftp at orion.math.uwaterloo.ca in the directory\\ pub/henry/reports or with URLs: ftp://orion.uwaterloo.ca/pub/henry/reports/distmat.ps.gz\\ or http://orion.math.uwaterloo.ca:80/\~{ }hwolkowi/henry/reports/distmat.ps.gz.} %\tableofcontents \begin{quote} \begin{center} Dedicated to Olvi Mangasarian \end{center} The first time that I came across Olvi's work was as a graduate student in the 70's when I studied from his book on Nonlinear Programming (now a SIAM classic) and also used the {\em Mangasarian-Fromovitz constraint qualification}. This is {\em the} constraint qualification (CQ) in nonlinear programming (NLP), since it guarantees the existence of Lagrange multipliers and is equivalent to stability of the NLP. This CQ has since been extended to various generalizations of NLP and plays a crucial role in perturbation theory. In 1983 I was a visitor at The University of Maryland, College Park, and was teaching a course in the Business College. While walking through the halls one day I noticed the name Fromovitz on one of the doors. I could not pass this by and knocked and asked if this was {\em the} Fromovitz. The reply was ``yes''; and, this is the story of the now famous CQ. Stan Fromovitz had just received his PhD from Stanford and was working at Shell Development Company in the Applied Math Dept. Olvi needed a special Theorem of the Alternative for his work on a CQ. Stan went digging into the Math Library at Berkeley and came up with exactly what was needed: Motzkin's Theorem of the Alternative. The end result of this was the MF CQ. I have followed Olvi's work very closely throughout my career. His work is marked by many beautiful and important results in various areas. Some of the ones I am aware of are: condition numbers for nonlinear programs; generalized convexity; complementarity problems; matrix splittings; solution of large scale linear programs. It is a pleasure and an honour to be able to contribute to this special issue. \flushright Henry \end{quote} \section{INTRODUCTION} \label{intro} Euclidean distance matrices have received a lot of attention in recent years both for their elegance and for their many important applications. Two other research areas of high interest currently, are semidefinite programming and interior-point methods. In this paper we solve the Euclidean distance matrix completion problem by generalizing the completion problem to allow for approximate completions, i.e. we find a weighted, closest Euclidean distance matrix. In particular, we introduce a primal-dual interior-point algorithm that solves an equivalent (quadratic objective) semidefinite programming problem. An $n \times n$ symmetric matrix $D=(d_{ij})$ with nonnegative elements and zero diagonal is called a {\em pre-distance matrix} (or dissimilarity matrix). In addition, if there exist points $x^1,x^2,\ldots,x^n$ in $\Re^r$ such that \beq d_{ij} = {\| x^i- x^j\|}^2, \1 \1 \1 i,j=1,2,\ldots,n, \eeq then $D$ is called a {\em Euclidean distance matrix} (EDM). The smallest value of $r$ is called {\em the embedding dimension} of $D$. Note that $r$ is always $\leq n-1$. Given a partial symmetric matrix $A$ with certain elements specified, the Euclidean distance matrix completion problem (EDMCP) consists in finding the unspecified elements of $A$ that make $A$ a EDM. Alternatively, for the approximate EDMCP, let $A$ be a pre-distance matrix, $H$ be an $n \times n$ symmetric matrix with nonnegative elements, and let $\|A\|_F= \sqrt{ \trace A^tA}$ denote the {\em Frobenius norm} of $A.$ Consider the objective function \[ f(D) := {\| H \circ (A - D) \|}^2_F, \] where $\circ$ denotes {\em Hadamard product}, e.g. \cite{HoJo:85}. The weighted, {\em closest Euclidean distance matrix problem} is \[ (CDM_0) \bt{ccc} $\mu^*$ := & $\min$ & $f(D)$ \\ & $\mbox{ subject to }$ & $D \in \E, $ \et \] where $\E$ denotes the cone of EDMs. Applications of EDMs abound: e.g. molecular conformation problems in chemistry \cite{ch88,PardShallXue}; multidimensional scaling and multivariate analysis problems in statistics \cite{lh82,MR94c:86008}; genetics, geography, and others, \cite{Alhom:93}. Many of these applications require a low embedding dimension, e.g. $r=3.$ Theoretical properties of EDMs can be found in e.g. \cite{MR96a:15025,MR88k:15023,gow85,hwlt91,jt95,infieldsLaur:97,sch35}. This includes characterizations as well as graph theoretic conditions for existence of completions. More information can be found in the recent survey article by Laurent \cite{infieldsLaur:97}. (Generalizations of EDM arise in \cite{MR96m:51020}.) An interesting discussion on algorithms for EDMCP appears in \cite{Trosset:97}. The point is made that there is no definitive general algorithm for EDMCP, i.e. one cannot provide an efficient decision rule for the question of whether a completion exists or not. However, there are many algorithms that find approximate completions. In \cite{Trosset:97,Trosset:97a,Trosset:97b}, the author presents results on finding EDM completions based on spectral decompositions. In particular, the computationally hard problem of fixing the rank (the embedding dimension) is discussed. Some work on finding the closest EDM to a given symmetric matrix appears in \cite{MR91j:90074,MR89b:15040,MR1459991}. Another approach based on global optimization and allowing for intervals for the distances is given in \cite{MoreWu:96,MoreWu:97} and also in \cite{ZoBySch:96}. We build a convex (tractable) model by relaxing the constraint on the embedding dimension (rank). Semidefinite programming, SDP, is an extension of linear programming where nonnegativity constraints on vector variables is replaced by positive semidefiniteness constraints on matrix variables. This area has attracted a lot of interest recently because of the many applications as well as the elegant mathematics involved, see e.g. the survey papers: \cite{Goem:97,VanBoy:94,VanBoy:94b,Al:94}. A lot of the interest in SDP is from the interior-point community who have completed so much successful work on linear programming. At the moment, interior-point methods are the most successful algorithms for general SDP problems, see e.g. the above survey articles as well as the books \cite{int:Nesterov5,Yinye:97b} and the recent theses \cite{Al:91,Alhom:93,KLERK:97,Sturm:97,Helmberg:94,Pat:96,Ram:93}. The above references provide some evidence of the current high level of research activity in these areas. The main contribution of this paper is a new approach to solving EDMCP. This approach is different than those in the literature in two ways. First we change the EDMCP into an approximation problem. This latter problem is a convex problem, i.e. our model is a tractable (polynomial time) model. Moreover, we further relax the original problem since we do not fix the embedding dimension. Thus, we do not solve problems that require a given dimension, though, we hope that this approach could be used as a first approximation to such problems. (A discussion on reducing the dimension is given in the conclusion, Section \ref{sect:concl}.) Our randomly generated tests consistently resulted in optimal solutions with low dimension, e.g. $r \leq 3,$ when a completion did not exist. This is in contrast to the case where completions do exist, where the interior point method found the completion of maximum rank, since interior-point methods find optimal solutions in the relative interior of the optimal face. Second, we use a semidefinite programming formulation and a primal-dual interior-point algorithm to solve the approximation problem. And, we prove that the Slater constraint qualification holds for our model if and only if the graph of the matrix of weights is connected, see Corollary \ref{cor:slater}. Usually, the lack of the constraint qualification results in numerical difficulties due to unbounded optimal sets. However, in our case we can take advantage of a disconnected graph to replace the original problem by two smaller simpler problems. As a side issue, we point out that our algorithm uses a new search direction for semidefinite programming introduced in \cite{KrReWo:96}. This search direction is based on applying a Gauss-Newton approach to the optimality conditions. Our purpose is not to compare different search directions, and other public domain packages may be used for the problems that we have solved; though our tests show that our approach is comparable. We use this approach, rather than the more standard approaches already in the literature, since it is very well suited for our particular application. And, we think that having a program that is specific for this problem has many advantages for e.g. exploiting structure. Moreover, the details of the algorithm are self-contained in this paper; and, we provide a MATLAB program for those interested in our tests. Numerical results are included which illustrate the efficiency and robustness of the interior-point approach. The paper is organized as follows. In Section~\ref{sect:distgeom} we introduce the basic results of EDM's. In Section~\ref{sect:dualnopt} we present the optimality conditions for the problem. In Section~\ref{sect:pdalgor} we derive the algorithm and the Slater constraint qualification result. We conclude with several remarks and numerical tests in Section~\ref{sect:concl}. In addition, we include Section \ref{app:A} with some technical details on the SDP algorithm. \section{DISTANCE GEOMETRY} \label{sect:distgeom} It is well known, e.g. \cite{sch35,gow85,hwlt91,MR14:889k}, that a pre-distance matrix $D$ is a EDM if and only if $D$ is negative semidefinite on \[ M:=\left\{ x \in \Re^n : x^t e = 0 \right\}, \] the orthogonal complement of $e$, where $e$ is the vector of all ones. Thus the set of all EDMs is a convex cone, which we denote by $\E.$ We exploit this result to translate the cone $\E$ to the cone of semidefinite matrices in ${\cal S}_{n-1},$ the space of symmetric matrices of order $n-1.$ Define the $n \times n$ orthogonal matrix \beq \label{defV} Q:= \left[ \frac{1}{\sqrt{n}}e \; | \; V \right], \1 \1 Q^tQ = I. \eeq Thus $V^te=0$ and $~V^tV=I$. Moreover, the subspace $M$ can be represented as the range of the $n \times (n-1)$ matrix $V$ and \beq \label{eq:Vmp} J := V V^t= I- \frac{e e^t}{n} \eeq is the orthogonal projection onto $M$. Now define the {\em centered} and {\em hollow} subspaces of $\Sn$ \beq \begin{array}{rcl} \Sc &:=& \{ B \in \Sn : Be = 0 \}, \\ \Sh& := & \{ D \in \Sn : \diag(D) = 0 \}, \end{array} \eeq where diag($D$) denotes the column vector formed from the diagonal of $D$. Following \cite{cri88}, we define the two linear operators acting on $\Sn$ \beq \begin{array}{rcl} \label{KK} \KK(B)& := & \mbox{diag}(B)\,e^t + e \, \mbox{diag}(B)^t - 2B, \end{array} \eeq and \beq \begin{array}{rcl} \label{T} \TT(D)& := & -\frac 12 JDJ. \end{array} \eeq The operator $- 2 \TT$ is an orthogonal projection onto $\Sc;$ thus it is a self-adjoint idempotent. \begin{thm} \label{KT} The linear operators satisfy \begin{eqnarray*} \KK ( \Sc) = \Sh, \\ \TT ( \Sh) = \Sc, \end{eqnarray*} and $\KK_{|\Sc}$ and $\TT_{|\Sh}$ are inverses of each other. \end{thm} \bpr See \cite{gow85,jt95}. \epr It can easily be verified that \beq \begin{array}{rcl} \label{K*T*} \KK^*(D)& = & 2 \left( \Diag(De)-D \right) \\ \end{array} \eeq is the adjoint operator of $\KK,$ where $\Diag(De)$ denotes the diagonal matrix formed from the vector $De$. In addition, a hollow matrix $D$ is EDM if and only if $B=\TT(D)$ is positive semidefinite, denoted by $B \succeq 0$. (We denote positive definiteness by $B \succ 0.$) Equivalently, $D$ is EDM if and only if $D=\KK(B),$ for some $B$ with $Be=0$ and $B \succeq 0$. In this case the embedding dimension $r$ is given by the rank of $B$. Moreover, if $B=XX^t$, then the coordinates of the points $x^1,x^2,\ldots,x^n$ that generate $D$ are given by the rows of $X$ and, since $Be=0,$ it follows that the origin coincides with the centroid of these points. For these and other basic results on EDM see e.g. \cite{gow85,hwlt91,jt95,infieldsLaur:97,sch35}. We now introduce the composite operators \beq \begin{array}{rcl} \label{KV} \KK_V(X)& := & \KK( V X V^t), \end{array} \eeq and \beq \begin{array}{rcl} \label{TV} \TT_V(D)& := & V^t\TT( D)V= - \frac 12 V^t D V, \end{array} \eeq where $V$ is defined in (\ref{defV}). \begin{lem} \label{KVTV} \begin{eqnarray*} \KK_V ( \snn) =\Sh, \\ \TT_V ( \Sh) =\snn, \end{eqnarray*} and $\KK_V$ and $\TT_V$ are inverses of each other on these two spaces. \end{lem} \bpr This immediately follows from Theorem \ref{KT} and the definition of $V.$ \epr From (\ref{K*T*}), we get that \beq \begin{array}{rcl} \label{K*VT*V} \KK^*_V(D)& = & V^t \KK^*(D) V \\ \end{array} \eeq is the adjoint operator of $\KK_V$. The following corollary summarizes the relationships between $\E$, the cone of Euclidean distance matrices of order $n$, and $\p$, the cone of positive semidefinite matrices of order $n-1$, which we will need for our model problem. \begin{cor} \label{co} Suppose that $V$ is defined as in (\ref{defV}). Then: \begin{eqnarray*} \KK_V(\p) & = & \E , \\ \TT_V(\E) & = &\p. \end{eqnarray*} \end{cor} \bpr We saw earlier that $D$ is EDM if and only if $D=\KK(B)$ with $Be=0$ and $B \succeq 0$. Let $X=V^tBV$, then since $Be=0$ we have $B=VXV^t$. Therefore, $VXV^t \succeq 0$ if and only if $X \succeq 0$; and the result follows using (\ref{KV}) and Lemma \ref{KVTV}. \epr Note that the $n \times (n-1)$ matrix $V$ as defined in (\ref{defV}) is not unique. In our code we use \beq \label{eq:Vorth} V := \left[ \begin{array}{cccc} y & y& \ldots & y \\ 1+x & x & \ldots & x \\ x & 1+x & \ldots & x \\ ... &... & \ddots &... \\ x & x & \ldots & 1+x \\ \end{array} \right], \eeq where $x = \frac{-1}{n + \sqrt{n} }$ and $y= \frac{-1}{\sqrt{n}}$. With this choice, it can be easily verified that $V^te=0$, $V^tV=I$, and $VV^t=J$ as required by (\ref{defV}). \subsection{ Program Formulations} \label{sect:formul} Since diag($A$) = diag($D$) = 0, we can assume without loss of generality that diag($H$) = 0. Note that $H_{ij}=0$ means that $D_{ij}$ is free, while $H_{ij}>0$ forces $D_{ij}$ to be approximately equal to $A_{ij}$. i.e., $A_{ij}$ is approximately fixed. If we want $D_{ij}=A_{ij}$ exactly, then we can add a linear constraint to the program. (See below.) Recall that the graph of $H$ is {\em connected} if for all indices $i \neq j$ there is a {\em path} of indices $i_1, i_2, ..., i_k$ such that $H_{i,i_1}\neq 0, H_{i_1,i_2}\neq 0, \ldots, H_{i_{k-1},i_{k}}\neq 0, H_{i_k,j}\neq 0,$ see e.g. \cite{BrualRys:91}. Thus, we can assume that the graph of $H$ is connected or the problem can be solved more simply as two smaller problems, see Lemma~\ref{lem:connect} below. In particular, we can assume that $H$ does not have a row (hence a column) of all zeros; otherwise the corresponding row and column in $A$ and $D$ are free (independent) and the problem can be posed in a lower dimensional space. By abuse of notation, let the function \[ f(X) := {\| H \circ (A - \KK_V ( X)) \|}^2_F = {\| H \circ \KK_V(B - X) \|}^2_F, \] where $B = \TT_V(A)$. We now apply Corollary \ref{co} and get the following problem, equivalent to $(CDM_0).$ \[ (CDM) \bt{ccc} $\mu^*$ := & $\min$ & $f(X)$ \\ & subject to & $\A X=b $ \\ & & $X \succeq 0.$ \et \] We allow for an additional constraint using the linear operator $\A : \snn \longrightarrow \Re^m$, i.e., $b \in \Re^m.$ The addition of this linear operator could represent some of the fixed elements in the given matrix $A$, e.g. adding the constraint $\left(K_V(X)\right)_{ij}=A_{ij}$ fixes the $ij$ element of $D.$ Also, note that $X \in \snn$. It is in this lower dimensional space that we solve the problem. We can recover the optimal distance matrix using the optimal $X$ and the relation \[ D = \KK_V(X). \] Using finite precision, we can never solve the approximation problem exactly. In addition, we need to calculate the embedding dimension. The following lemma shows we lose little in the objective function if we choose a small embedding dimension using a numerical rank approach, i.e. if we only discard very small eigenvalues, then the objective function changes very little. \begin{lem} Suppose that $X^*$ solves (CDM). Let $\bar{X}$ be the closest symmetric matrix to $X^*$ with rank $k$, i.e. we set the smallest $n-k$ eigenvalues of $X^*$ to 0, $\lambda_{k+1}= \ldots \lambda_n=0.$ Then, \beq \label{eq:error} \sqrt{f(\bar{X})}\leq \sqrt{f(X^*)}+ 2 \gamma \left(\sqrt{n}+1\right) \sqrt{ \sum_{i=k+1}^n \lambda_i^2}, \end{equation} where $\gamma:= \max_{ij}H_{ij}.$ \end{lem} \bpr By the Cauchy-Schwartz inequality, \[ \sqrt{f(\bar{X})} \leq \sqrt{f(X^*)}+ \gamma ||{\cal K}(V(X^*-\bar{X})V^t)||_F. \] The result now follows from the definition of $\cal K.$ More precisely, let $B= V(X^*-\bar{X})V^t.$ Then \begin{eqnarray*} || {\cal K}(V(X^*-\bar{X})V^t)||_F &=& || {\cal K}(B)||_F \\ &\leq & 2||B||_F +2||\diag(B)e^t||_F\\ & \leq & 2||B||_F +2\sqrt{n} ||B||_F. \end{eqnarray*} \epr \section{DUALITY and OPTIMALITY} \label{sect:dualnopt} We now derive the optimality conditions and duality theory needed for a primal-dual interior-point approach. For $\Lambda \in \snn$ and $y \in R^m$, let \beq L(X,y,\Lambda) = f(X) + \langle y, b - \A(X) \rangle - \trace \Lambda X \eeq denote the {\em Lagrangian} of (CDM). It is easy to see that the primal program (CDM) is equivalent to \beq \mu^* = \min_X \max_{\stackrel{y}{\Lambda \succeq 0}} L(X,y,\Lambda) = \min_{X \succeq 0} \max_{\stackrel{y}{\Lambda \succeq 0}} L(X,y,\Lambda). \eeq We assume that the generalized Slater's constraint qualification, \[ \exists X \succ 0 \mbox{ with } \A(X)= b, \] holds for (CDM). Slater's condition implies that strong duality holds, i.e., \beq \mu^* = \max_{\stackrel{y}{\Lambda \succeq 0}} \min_{X} L(X,y,\Lambda) = \max_{\stackrel{y}{\Lambda \succeq 0}} \min_{X\succeq 0} L(X,y,\Lambda), \eeq and $\mu^*$ is attained for some $y$ and $\Lambda \succeq 0$, see e.g. \cite{w11}. Since the semidefinite constraint on $X$ can be treated as redundant, the inner minimization of the convex, in $X$, Lagrangian is unconstrained and we can differentiate to get the equivalent problem \beq \mu^*= \max_{\stackrel{\nabla f(X) - \A^*y=\Lambda}{\Lambda \succeq 0}} f(X) + \langle y, b - \A(X) \rangle - \trace \Lambda X. \eeq We can now state the dual problem \beq \label{eq:dualprob} \bt{cccc} (DCDM) & $\mu^*$ :=& max & $f(X)+\langle y, b - \A(X) \rangle - $ \trace $ \Lambda X$ \\ & & subject to & $\nabla f(X) - \A^*y- \Lambda = 0$ \\ & & & $ \Lambda \succeq 0, (X \succeq 0). $ \et \eeq We keep the semidefinite constraint on $X$ in brackets to emphasize that it is a hidden constraint in the dual, though it is a given constraint in the primal. The above pair of primal and dual programs, (CDM) and (DCDM), provide an optimality criteria in terms of feasibility and complementary slackness. This provides the basis for many algorithms including primal-dual interior-point algorithms. In particular, we see that the duality gap, in the case of primal and dual feasibility, is given by the complementary slackness condition: \beq \trace X (2\KK^*_V( H^{(2)} \circ \KK_V( \bar{X}-B))- \A^*\bar{y}) = 0, \eeq or equivalently \[ X (2\KK^*_V( H^{(2)} \circ \KK_V( \bar{X}-B))- \A^* \bar{y}) = 0, \] where $ H^{(2)}=H \circ H.$ \begin{thm} \label{opt} Suppose that Slater's condition holds. Then $\bar{X} \succeq 0$, and $\bar{y}$, $\bar{ \Lambda} \succeq 0$ solve (CDM) and (DCDM), respectively, if and only if the following three equations hold. \[ \bt{cc} $\A (\bar{X}) = b $ & primal feasibility \\ $ 2\KK^*_V( H^{(2)} \circ \KK_V( \bar{X}-B)) - \A^* \bar{y} - \bar{\Lambda} =0$ & dual feasibility \\ \trace $\bar{ \Lambda} \bar{X}=0$ & complementary slackness \et \] \end{thm} In the case that there is no linear operator $\cal A,$ we get a simplified variational principle. The statement that no feasible direction is a descent direction translates into the following characterization of an optimal solution $\bar{X}$ of (CDM): \beq \nabla f(\bar{X}) \in (\PP - \bar{X})^+, \label {Ps} \eeq where $\PP$ is the cone of positive semidefinite matrices and \[ S^+ = \{ P : \trace QP \geq 0, \1, \forall Q \in S\} \] is the {\em polar cone} of the set $S$. This yields the following characterization of optimality. \begin{thm} Suppose that (CDM) has no linear constraint $\cal A.$ The matrix $\bar{X} \succeq 0 $ solves (CDM) if and only if \beq \trace \KK^*_V( H^{(2)} \circ \KK_V ( \bar{X} - B))\; (X-\bar{X}) \geq 0, \1 \forall X \succeq 0. \eeq \end{thm} \bpr Note that the gradient acting on the symmetric matrix $h$, is \beq \bt{ll} $\langle \nabla f(X),h \rangle $ & = 2 \trace $(H^{(2)} \circ \KK_V(\bar{X} - B)) \KK_V(h)$ \\ & = 2 \trace $ \KK^*_V(H^{(2)} \circ \KK_V(\bar{X}-B))h.$ \et \eeq Therefore the gradient of the objective function is \beq \nabla f(X) = 2 \KK^*_V(H^{(2)} \circ \KK_V( \bar{X} - B)). \eeq The result follows upon replacing $h$ by the direction $X - \bar{X}$, for $X \succeq 0$, and applying the so-called Pshenichnyi condition (\ref{Ps}). \epr \section{PRIMAL-DUAL INTERIOR-POINT\\ALGORITHM} \label{sect:pdalgor} We now derive a primal-dual interior point method using the log-barrier approach \cite{HeReVaWo:93}. This is an alternative way of deriving the optimality conditions in Theorem \ref{opt}. For simplicity, we treat the unconstrained problem. i.e., we consider (CDM) without the linear equality constraint $\A X = b$. In this case, the Slater constraint qualification holds for both primal and dual problems, as can be seen from the following lemma. (In the case that some of the elements of $A$ are definitely fixed, then one needs to rederive the algorithm and include the constraint $\A X = b$.) \begin{lem} \label{lem:connect} Let $H$ be an $n \times n$ symmetric matrix with nonnegative elements and 0 diagonal such that the graph of $H$ is connected. Then \[ \KK_V^*(H^{(2)} \circ \KK_V(I)) \succ 0 , \] where $I \in \snn$ is the identity matrix. \end{lem} \bpr A simple calculation shows that $\KK_V(I) = 2 (ee^t-I)$. Thus $H^{(2)} \circ \KK_V(I)= 2 H^{(2)}$ and $ \KK_V^*(H^{(2)} \circ \KK_V(I))= 4V^t $(Diag$(H^{(2)}e) - H^{(2)} )V$. Note that $e$ is an eigenvector of Diag$(H^{(2)}e) - H^{(2)}$ with zero eigenvalue. We show that 0 is a simple eigenvalue. Assume to the contrary that there exists another nonzero eigenvector $u$, i.e. $u$ is not a multiple of $e$. Then \[u^t (\mbox{Diag}(H^{(2)}e)-H^{(2)})u= \sum_{i 0$ can be chosen such that \beq \label{eq:Xfeas} X=B+ \alpha I \succ 0; \end{equation} we can then apply the previous Lemma \ref{lem:connect} to get \beq \label{eq:Lamfeas} \KK_V^*(H^{(2)} \circ \KK_V(X-B)) \succ 0. \end{equation} \begin{cor} \label{cor:slater} Slater's constraint qualification holds for the dual problem (\ref{eq:dualprob}) if and only if the graph of $H$ is connected. \end{cor} \bpr Sufficiency follows directly from the above Lemma \ref{lem:connect}. To prove necessity, suppose that the graph is not connected. Then (CDM) can be solved as two smaller disconnected problems where the distances joining points for these two smaller problems are completely free, i.e. the set of optimal distances matrices for (CDM$_0$) is unbounded. Therefore, the set of optimal solutions $X$ of (CDM) is unbounded. This implies that Slater's condition cannot hold in the dual, see e.g. the proof of \cite[Theorem 4]{Golds:72} or \cite{LempMaur:80}. In our case, it is rather simple to prove this result and we include it for completeness, i.e. we prove that Slater's condition for the dual implies that the optimal set for the primal is bounded. Suppose that Slater's condition holds for the dual, i.e. there exists $\bar{\Lambda} \succ 0$ and $\bar{X} \succeq 0$ such that $\nabla f(\bar{X}) - \bar{\Lambda} = 0.$ Equivalently, this means that $\nabla f(\bar{X}) \succ 0.$ Let $X^*$ be optimal for (CDM). Then convexity of $f$ implies that \[ \left< \nabla f(\bar{X}),X^*-\bar{X} \right> \leq 0. \] Therefore, we get \[ \begin{array}{rcl} \left< \nabla f(\bar{X}),\bar{X} \right> & \geq & \left< \nabla f(\bar{X}),X^* \right>\\ & \geq & \lambda_{\min}(\nabla f(\bar{X})) \trace X^* \\ & \geq & \lambda_{\min}(\nabla f(\bar{X})) \lambda_{\max}(X^*), \end{array} \] i.e. the norm of $X^*$ is bounded by $ \frac n{\lambda_{\min}(\nabla f(\bar{X}))} \left< \nabla f(\bar{X}),\bar{X} \right>.$ \epr The log-barrier problem for (CDM) is \[ \min_{X \succ 0} B_{\mu}(X) := f(X) - \mu \log \det(X), \] where $\mu \downarrow 0$. For each $\mu > 0$ we take one Newton step for solving the stationarity condition \beq \nabla B_{\mu}(X) = 2 \KK^*_V(H^{(2)} \circ \KK_V(X - B)) - \mu X^{-1}=0. \eeq Let \beq C := 2 \KK^*_V(H^{(2)} \circ \KK_V( B)) = 2 \KK^*_V(H^{(2)} \circ A). \eeq Then the stationarity condition is equivalent to \beq \nabla B_{\mu}(X) = 2 \KK^*_V\left(H^{(2)} \circ \KK_V(X)\right) - C - \mu X^{-1}=0. \eeq By equating $\Lambda = \mu X^{-1},$ and multiplying through by $X,$ we get the optimality conditions, $F:=\left( \begin{array}{c} F_d \\ F_c \end{array} \right)=0,$ \beq \label{cs} \begin{array}{llccl} F_d&:= 2 \KK^*_V\left(H^{(2)} \circ \KK_V(X)\right) - C - \Lambda &=&0 & \mbox{dual feas.} \\ F_c&:= \Lambda X - \mu I &=&0 & \mbox{perturbed compl. slack.}, \end{array} \eeq and the estimate of the barrier parameter \beq \mu = \frac{1}{n-1} \trace \Lambda X. \eeq Following, we present the p-d i-p framework we used. This is a common framework for both linear and semidefinite programming, see e.g. \cite{SWright:96}. We include a centering parameter $\sigma_k$ (rather than the customary predictor-corrector approach) and let ${\cal F}^0$ denote the set of strictly feasible primal-dual points; $F^\prime$ denotes the derivative of the function of optimality conditions. \begin{alg} \label{alg:pdip} (p-d i-p framework:)\\ {\bf Given} $(X^0,\Lambda^0) \in {\cal F}^0$\\ {\bf for} $k=0,1,2 \ldots $\\ \hspace{.5in}{\bf solve} for the search direction (in a least squares sense) \\ \[ ~~~~F^{\prime}(X^k,\Lambda^k) \left( \begin{array}{ccc} \delta X^k \\ \delta \Lambda^k \end{array} \right) = \left( \begin{array}{ccc} -F_d \\ \-\Lambda^k X^k+ \sigma_k \mu_k I \end{array} \right) \] \hspace{.5in}where $\sigma_k$ centering, $\mu_k=\trace X^k\Lambda^k/(n-1)$ \[ (X^{k+1},\Lambda^{k+1}) = (X^k,\Lambda^k) + \alpha_k (\delta X^k, \delta \Lambda^k), ~~ \alpha_k > 0, \] \hspace{.5in}so that $(X^{k+1},\Lambda^{k+1}) \succ 0$\\ {\bf end (for)}. \end{alg} We use two approaches for solving the least squares problem. First we solve the large least squares problem; we call this the Gauss-Newton (GN) method. In the second approach, we restrict dual feasibility and substitute for $\delta \Lambda^k$ in the second equation; we call this the restricted Gauss-Newton (RGN) method. The numerical tests show that the second approach is significantly more efficient. \section{CONCLUSION and COMPUTATIONAL RESULTS} \label{sect:concl} In this paper we have presented an interior-point algorithm for finding the weighted, closest Euclidean distance matrix; this provides a solution for the approximate Euclidean distance matrix completion problem. The algorithm has been extensively tested and has proven to be efficient and very robust, i.e. it has not failed on any of our test problems. In addition, an important observation is that the ranks of the optimal solutions $X,$ for sparse problems where no completion existed, were typically between 1 and 3, i.e. a very small embedding dimension was obtained without any original rank restrictions. However, when a completion existed, then typically the embedding dimension was high. In this case, we can apply the technique presented in \cite{AlWo:98} to ``purify'', i.e. to iteratively move to an optimal solution of smaller rank on the optimal face. (See also Pataki \cite{Pat:96}.) We discuss more details of the Gauss-Newton approach in Section \ref{app:A}. The program was written in MATLAB. The tests were done on randomly generated problems; we used s SPARC 20 with 64 megs of RAM and SPECint 65.3, SPECfp 53.1. The documented MATLAB code, as well as ongoing test results, can be obtained with URL (or anonymous ftp) \\ ftp://orion.math.uwaterloo.ca/pub/henry/software/distance.d, or\\ http://orion.math.uwaterloo.ca/\~{ }hwolkowi/henry/software/distance.d. In Table \ref{table:closest} we present a sample of test results. These results include problems with matrices up to dimension $n=42.$ We would like to emphasize that these are just preliminary test results. For example, we do not use a predictor-corrector approach, which has become the standard approach in interior-point methods, but rather use a centering parameter. Using the predictor-corrector approach should reduce the number of iterations from 20-30\%. Tests with larger sized matrices are in progress. \begin{table}[htb] \begin{center} \begin{tabular}{|c|c|c|c|c|c|c|c|} \hline dim & toler & $H$ dens. & $rank(X)$ & {iterations} & \multicolumn{2}{c|} {lss cpu-time} \\ \cline{6-7} & & & & & GN& RGN\\ \hline 8 & $10^{-13}$& .8 & 2 & 25&.16 &.1 \\ 9 & $10^{-13}$& .8 & 2 & 23 & .24 &.13 \\ 10 & $10^{-13}$& .8 & 3 & 25 &.34 &.18 \\ 12 & $10^{-9}$& .5 & 3 & 17 &.73 &.32 \\ 15 & $10^{-9}$& .5 & 2 & 20 &2.13 &.79 \\ 18 & $10^{-9}$& .5 & 4 & 20 &6.15 & 1.9 \\ 20 & $10^{-9}$& .3 & 2 & 20 & 11.35 &3.3 \\ 24 & $10^{-9}$& .3 & 2 & 20 &34.45 & 8.4 \\ 30 & $10^{-9}$& .3 & 4 & 20 & 138. & 31.5 \\ 35 & $10^{-9}$& .2 & 3 & 19 & 373. & 77. \\ 38 & $10^{-9}$& .2 & 3 & 19 & 634. & 127 \\ 40 & $10^{-8}$& .1 & 2 & 20 & 845.9 & 181.7 \\ 42 & $10^{-8}$& .1 & 4 & 18 & 1118.& 232.02 \\ \hline \end{tabular} \caption{\label{table:closest} data for closest distance matrix: dimension; tolerance for duality gap; density of nonzeros in $H$; rank of optimal $X$; number of iterations; cpu-time for one least squares solution of the GN and restricted GN directions. } \end{center} \end{table} We conclude with a specific example and its graph: the dimension is 11 generated with sparsity $0.5;$ the optimal value is $0.9256;$ the rank of $X$ is 3; and the number of iterations is 25 to get 13 decimals accuracy. The matrices $H$ and $A$ and the optimal distance matrix $D$ corresponding to the graph in Figure 1 are, respectively: \begin{small} \begin{verbatim} matrix H = 0 0 0 0 0 5 2 5 0 7 0 0 0 0 3 4 4 0 0 1 0 5 0 0 0 0 0 0 1 0 0 0 2 0 3 0 0 0 1 0 0 3 0 3 0 4 0 0 0 2 7 0 3 0 0 5 4 0 1 2 0 2 0 0 0 7 2 0 1 0 7 2 0 6 0 1 2 5 0 0 0 0 0 6 0 0 0 0 0 1 0 3 3 0 0 0 0 3 0 7 0 0 0 0 0 1 0 3 0 0 0 5 2 3 0 7 2 0 0 0 0 -------------------------- matrix A = 0 0 0 0 0 4 0 0 6 6 2 0 0 0 8 4 0 2 0 0 0 0 0 0 0 5 6 0 6 0 0 0 7 0 8 5 0 0 4 1 4 5 4 0 0 4 6 0 0 0 0 0 5 0 0 4 0 0 4 0 0 1 0 0 0 0 0 2 6 1 0 1 0 0 0 0 3 0 0 0 4 0 0 0 0 0 0 0 6 0 0 5 5 0 0 0 0 0 0 6 0 0 4 0 0 0 0 0 0 0 2 0 7 0 0 0 3 0 0 0 0 \end{verbatim} \end{small} \begin{small} \begin{verbatim} -------------------------- matrix D = Columns 1 through 7 0 6.9200 7.2655 6.9713 0.7190 3.9912 0.5989 6.9200 0 8.1310 6.6123 3.8224 0.6523 3.5381 7.2655 8.1310 0 9.3511 5.9827 6.4399 6.0000 6.9713 6.6123 9.3511 0 3.7085 3.6117 4.7459 0.7190 3.8224 5.9827 3.7085 0 1.4994 0.0775 3.9912 0.6523 6.4399 3.6117 1.4994 0 1.4981 0.5989 3.5381 6.0000 4.7459 0.0775 1.4981 0 0.1684 4.9773 6.4452 5.6419 0.2296 2.5198 0.1321 8.9782 0.4804 9.0262 5.0302 4.8968 0.9976 4.9394 5.9737 0.3183 7.3517 4.0355 2.7767 0.1996 2.7897 5.2815 1.0844 7.0000 2.3417 2.1414 0.2465 2.3892 Columns 8 through 11 0.1684 8.9782 5.9737 5.2815 4.9773 0.4804 0.3183 1.0844 6.4452 9.0262 7.3517 7.0000 5.6419 5.0302 4.0355 2.3417 0.2296 4.8968 2.7767 2.1414 2.5198 0.9976 0.1996 0.2465 0.1321 4.9394 2.7897 2.3892 0 6.6871 4.1359 3.5984 6.6871 0 0.3050 0.7459 4.1359 0.3050 0 0.2315 3.5984 0.7459 0.2315 0 \end{verbatim} \end{small} \begin{figure}[htb] \centering % \centerline{\ \psfig{figure=fig11513.ps,width=5.3in} \label{fig2} \caption{Approximate Completion Problem} \end{figure} \noindent \subsection{Gauss-Newton Direction} \label{app:A} The linear system for the search direction in Algorithm \ref{alg:pdip} is overdetermined. Therefore, it is not clear what is meant by solving this system. There are many different search directions that have been used for SDP, see e.g. \cite{Todd:97}. Our problem is not a standard SDP since it has a quadratic objective and no constraints. Therefore, the standard public domain packages do not apply directly. \footnote{It has been pointed out to us that the package SDPpack, \cite{AHNO97}, will solve our problem if we replace our original objective function with the variable $t$ and add a constraint $||W||_F \leq t,$ where $W$ represents the appropriate quantity in our original objective function. The numerical tests with this approach appear to be comparable to our approach. We thank Mike Overton for this observation.} In our approach we use the Gauss-Newton direction introduced in \cite{KrReWo:96}, i.e. we linearize (\ref{cs}) and find the least squares solution of the resulting overdetermined linear system. This direction has many good properties, e.g.: it always exists; the linear systems for both GN and RGN are nonsingular at each iteration and in the limit; quadratic convergence steps are expected since the optimal value of the nonlinear least squares problems are 0. (For the details of these and other properties see \cite{KrReWo:96}.) \bibliography{.psd,.master,.publs} \end{article} \end{document}