Convex Optimization Learning of Faithful Euclidean Distance Representations in Nonlinear Dimensionality Reduction

Classical multidimensional scaling only works well when the noisy distances observed in a high dimensional space can be faithfully represented by Euclidean distances in a low dimensional space. Advanced models such as Maximum Variance Unfolding (MVU) and Minimum Volume Embedding (MVE) use Semi-Definite Programming (SDP) to reconstruct such faithful representations. While those SDP models are capable of producing high quality configuration numerically, they suffer two major drawbacks. One is that there exist no theoretically guaranteed bounds on the quality of the configuration. The other is that they are slow in computation when the data points are beyond moderate size. In this paper, we propose a convex optimization model of Euclidean distance matrices. We establish a non-asymptotic error bound for the random graph model with sub-Gaussian noise, and prove that our model produces a matrix estimator of high accuracy when the order of the uniform sample size is roughly the degree of freedom of a low-rank matrix up to a logarithmic factor. Our results partially explain why MVU and MVE often work well. Moreover, we develop a fast inexact accelerated proximal gradient method. Numerical experiments show that the model can produce configurations of high quality on large data points that the SDP approach would struggle to cope with.


Introduction
The chief purpose of this paper is to find a complete set of faithful Euclidean distance representations in a low-dimensional space from a partial set of noisy distances, which are supposedly observed in a higher dimensional space.The proposed model and method thus belong to the vast field of nonlinear dimensionality reduction.Our model is strongly inspired by several high-profile Semi-Definite Programming (SDP) models, which aim to achieve a similar purpose, but suffer two major drawbacks: (i) theoretical guarantees yet to be developed for the quality of recovered distances from those SDP models and (ii) the slow computational convergence, which severely limits their practical applications even when the data points are of moderate size.Our distinctive approach is to use convex optimization of Euclidean Distance Matrices (EDM) to resolve those issues.In particular, we are able to establish theoretical error bounds of the obtained Euclidean distances from the true distances under the assumption of uniform sampling, which has been widely used in modelling social networks.Moreover, the resulting optimization problem can be efficiently solved by an accelerated proximal gradient method.In the following, we will first use social network to illustrate how initial distance information is gathered and why the uniform sampling is a good model in understanding them.We then briefly discuss several SDP models in nonlinear dimensionality reduction and survey relevant error-bound results from matrix completion literature.They are included in the first three subsections below and collectively serve as a solid motivation for the current study.We finally summarize our main contributions with notation used in this paper.

Distances in Social Network and Their Embedding
The study of structural patterns of social network from the ties (relationships) that connect social actors is one of the most important research topics in social network analysis (Wasserman, 1994).To this end, measurements on the actor-to-actor relationships (kinship, social roles, affective, transfer of information, etc) are collected or observed by different methods (questionnaires, direct observation, experiments, written records, etc) and the measurements on the relational information are referred as the network composition.In other words, without appropriate network measurements, we are not able to study any structural property.The measurement data usually can be presented as an n×n measurement matrix, where the n rows and the n columns both refer to the studied actors.Each entry of these matrices indicates the social relationship measurement (e.g., presence/absence or similarity/dissimilarity measure) between the row and column actors.In this paper, we are only concerned with symmetric relationships, i.e., the relationship from actor i to actor j is the same as that from actor j to actor i.Furthermore, there exist standard ways to convert the measured relationships into Euclidean distances, see (Cox and Cox, 2000, Section 1.3.5)and (Borg and Groenen, 2007, Chapter 6).
However, it is important to note that in practice, only partial relationship information are observed, which means the measurement matrix is usually incomplete and noisy.The observation processes are often assumed to follow certain random graph model.One simple but wildly used model is the Bernoulli random graph model (Solomonoff and Rapoport, 1951;Erdős and Rényi, 1959).Let n labeled vertices be given.The Bernoulli random graph is obtained by connecting each pair (or not) of vertices independently with the common probability p (or 1 − p) and it reproduces well some principal features of the real-world social network such as the "small-world" effect (Milgram, 1967;de Sola Pool and Kochen, 1979).Other properties such as the degree distribution, the connectivity, the diameters of the Bernoulli random graph, can be found in (e.g., Bollobás, 2001;Janson et al., 2011).For more details on the Bernoulli as well as other random models, one may refer to the review paper (Newman, 2003) and references therein.In this paper, we mainly focus on the Bernoulli random graph model.Consequently, the observed measurement matrix follows the uniform sampling rule, which will be described in Section 2.3.
In order to examine the structural patterns of a social network, the produced images (e.g., embedding in 2 or 3 dimensional space for visualization) should preserve the structural patterns as much as possible, as highlighted by Freeman (2005), "the points in a visual image should be located so the observed strengths of the inter-actor ties are preserved."In other words, the designed dimensional reduction algorithm has to assure that the embedding Euclidean distances between points (nodes) fit in the best possible way the observed distances in a social space.Therefore, the problem now reduces to whether one can effectively find the best approximation to the true social measurement matrix, which has a low embedding dimension, from the observed incomplete and noisy data.The classical Multidimensional Scaling (cMDS) (see Section 2.1) provides one of the most often used embedding methods in using distance information.However, cMDS alone is often not adequate to produce satisfactory embedding, as rightly observed in several high-profile embedding methods in manifold learning.

Embedding Methods in Manifold Learning
The cMDS and its variants have found many applications in data dimension reduction and have been well documented in the monographs (Cox and Cox, 2000;Borg and Groenen, 2007;Pȩkalska and Duin, 2005).When the distance matrix (or dissimilarity measurement matrix) is close to a true Euclidean Distance Matrix (EDM) with the targeted embedding dimension, cMDS often works very well.Otherwise, a large proportion of unexplained variance has to be cut off or it may even yield negative variances, resulting in what is called embedding in a pseudo-Euclidean space and hence creating the problem of unconventional interpretation of the actual embedding (see, e.g., Pȩkalska et al., 2002).
cMDS has recently motivated a number of high-profile numerical methods, which all try to alleviate the issue mentioned above.For example, the ISOMAP of Tenenbaum et al. (2000) proposes to use the shortest path distances to approximate the EDM on a low-dimensional manifold.The Maximum Variance Unfolding (MVU) of Weinberger and Saul (2006) through SDP aims for maximizing the total variance and the Minimum Volume Embedding (MVE) of Shaw and Jebara (2007) also aims for a similar purpose by maximizing the eigen gap of the Gram matrix of the embedding points in a low-dimensional space.The need for such methods comes from the fact that the initial distances either are in stochastic nature (e.g., containing noises) or cannot be measured (e.g., missing values).The idea of MVU has also been used in the refinement step of the celebrated SDP method for sensor network localization problems (Biswas et al., 2006).
It was shown in (Tenenbaum et al., 2000;Bernstein et al., 2000) that ISOMAP enjoys the elegant theory that the shortest path distances (or graph distances) can accurately estimate the true geodesic distances with a high probability if the finite points are chosen randomly from a compact and convex submanifold following a Poisson distribution with a high density, and the pairwise distances are obtained by the k-nearest neighbor rule or the unit ball rule (see Section 2.3 for the definitions).However, for MVU and MVE (both have enjoyed a numerical success), there exist no theoretical guarantee as to how good the obtained Euclidean distances are.At this point, it is important to highlight two observations.(i) The shortest path distance or the distance by the k-nearest neighbor or the unit-ball rule is often not suitable in deriving distances in social network.This point has been emphasized in the recent study on E-mail social network by Budka et al. (2013).(ii) MVU and MVE models only depend on the initial distances and do not depend on any particular ways in obtaining them.They then rely on SDP to calculate the best fit distances.From this point of view, they can be applied to social network embedding.This is also pointed out in Budka et al. (2013).Due to the space limitation, we are not able to review other leading methods in manifold learning, but refer to (Burges, 2009, Chapter 4) for a guide.
Inspired by their numerical success, our model will inherit the good features of both MVU and MVE.Moreover, we are able to derive theoretical results in guaranteeing the quality of the obtained Euclidean distances.Our results are the type of error bounds, which have attracted growing attention recently.We review the relevant results below.

Error Bounds in Low-Rank Matrix Completion and Approximation
As mentioned in the preceding section, our research has been strongly influenced by the group of researches that related to the MVU and MVE models, which have natural geometric interpretations and use SDP as their major tool.Their excellent performance in data reduction calls for theoretical justification, which in any type seems in nonexistence.Our model also enjoys a similar geometric interpretation, but departs from the two models in that we deal with EDM directly rather than reformulating it as SDP.This key departure puts our model in the category of matrix approximation problems, which have attracted much attention recently from machine learning community and motivated our research.
The most popular approach to recovering a low-rank matrix solution of a linear system is via the nuclear norm minimization (Mesbahi, 1998;Fazel, 2002), which is of SDP.What makes this approach more exciting and important is that it has a theoretically guaranteed recoverability (recoverable with a high probability).The first such a theoretical result was obtained by Recht et al. (2010) by employing the Restricted Isometric Property (RIP).However, for the matrix completion problem the sample operator does not satisfy the RIP (see, e.g., Candès and Plan, 2010).For the noiseless case, Candès and Recht (2009) proved that a low-rank matrix can be fully recovered with high probability provided that a small number of its noiseless observations are uniformly sampled.See Candès and Tao (2010) for an improved bound and the near-optimal bound on the sample number.By adapting the techniques from quantum information theory developed in Gross et al. (2010) to the matrix completion problem, a short and intelligible analysis of the recoverability was recently proposed by Recht (2011).
The matrix completion with noisy observations was studied by Candès and Plan (2010).Recently, the noisy case was further studied by several groups of researchers including Koltchinskii et al. (2011), Negahban andWainwright (2012) and Klopp (2014), under different settings.In particular, the matrix completion problem with fixed basis coefficients was studied by Miao et al. (2012), who proposed a rank-corrected procedure to generate an estimator using the nuclear semi-norm and established the corresponding non-asymmetric recovery bounds.
Very recently, Javanmard and Montanari (2013) proposed a SDP model for the problem of (sensor network) localization from an incomplete set of noisy Euclidean distances.Using the fact that the squared Euclidean distances can be represented by elements from a positive semidefinite matrix: where x i ∈ d are embedding points and X defined by X ij = x T i x j is the Gram matrix of those embedding points, the SDP model aims to minimize Tr(X) (the nuclear norm of X).Equivalently, the objective is to minimize the total variance x i 2 of the embedding points because it is commonly assumed that the embedding points are already centered.This objective obviously contradicts the main idea of MVU and MVE, which aim to make the total variance as large as possible.It is important to point out that making the variance as big as possible seems to be indispensable for SDP to produce high quality of localization.This has been numerically demonstrated in Biswas et al. (2006).
The impressive result in Javanmard and Montanari (2013) roughly states that the error bound reads as O((nr d ) 5 ∆ r 4 ) containing an undesirable term (nr d ) 5 , where r is the radius used in the unit ball rule, d is the embedding dimension, ∆ is the bound on the measurement noise and n is the number of embedding points.As pointed out by Javanmard and Montanari (2013) that the numerical performance suggested the error seems to be bounded by O( ∆ r 4 ), which does not match the derived theoretical bound.This result also shows tremendous technical difficulties one may have to face in deriving similar bounds for EDM recovery as those for general matrix completion problems, which have no additional constraints other than being low rank.
To summarize, most existing error bounds are derived from the nuclear norm minimization.When translating to the Euclidean distance learning problem, minimizing the nuclear norm is equivalent to minimizing the variance of the embedding points, which contradicts the main idea of MVU and MVE in making the variance as large as possible.Hence, the excellent progress in matrix completion/approximation does not straightforwardly imply useful bounds about the Euclidean distance learning in a low-dimensional space.Actually one may face huge difficulty barriers in such extension.In this paper, we propose a convex optimization model to learn faithful Euclidean distances in a low-dimensional space.We derive theoretically guaranteed bounds in the spirit of matrix approximation and therefore provide a solid theoretical foundation in using the model.We briefly describe the main contributions below.

Main Contributions
This paper makes two major contributions to the field of nonlinear dimensionality reduction.One is on building a convex optimization model with guaranteed error bounds and the other is on a computational method.
(i) Building a convex optimization model and its error bounds.Our departing point from the existing SDP models is to treat EDM (vs positive semidefinite matrix in SDP) as a primary object.The total variance of the desired embedding points in a low-dimensional space can be quantitatively measured through the so-called EDM score.The higher the EDM score is, the more the variance is explained in the embedding.Therefore, both MVU and MVE can be regarded as EDM score driven models.Moreover, MVE, being a nonconvex optimization model, is more aggressive in driving the EDM score up.However, MVU, being a convex optimization model, is more computationally appealing.Our model strikes a balance between the two models in the sense that it inherits the appealing features from them.It results in a convex optimization model whose objective consists of three parts (see Subsection 3.1 for more details).
Each part in the objective has its own purpose.The first part is a least-square term that governs the deviation from the observed distances.The second is the nuclear norm minimization term, which, as we have argued before, contradicts the idea of maximizing the total variance.Hence, the second term is corrected by the third term, which involves a set of orthogonal axes of approximating the true embedding space.
The last term is crucial to our analysis and is also to accommodate situations where the initial (valuable) information is available about the embedding space.To illustrate such a situation, we may simply think that the leading eigenvectors of the distance matrix estimator from ISOMAP form a good approximation to the true embedding space.When the three parts are combined, the model drives the EDM score up.
What makes our model more important is that it yields guaranteed error bounds under the uniform sampling rule.More precisely, we show that for the unknown n × n Euclidean distance matrix with the embedding dimension r and under mild conditions, the average estimation error is controlled by Crn log(n)/m with high probability, where m is the sample size and C is a constant independent of n, r and m.It follows from this error bound result that our model will produce an estimator with high accuracy as long as the sample size is of the order of rn log(n), which is roughly the degree of freedom of a symmetric hollow matrix with rank r up to a logarithmic factor in the matrix size.It is interesting to point out that with special choice of model parameters, our model reduces to one of the subproblems solved by MVE.Hence, our theoretical result partially explains why the MVE often leads to configurations of high quality.To our knowledge, it is the first such theoretical result that shed lights on the MVE model.
(ii) An efficient computational method.Treating EDM as a primary object not only benefits us in deriving the error-bound results, but also leads to an efficient numerical method.Previously, both MVU and MVE models have numerical difficulties when the data points are beyond 1000.They may even have difficulties with a few hundreds of points when their corresponding slack models are to be solved.This probably explains why most of publications in using the two models do not report cpu time.On the contrary, our model allows us to develop a very fast inexact accelerated proximal gradient method (IAPG) even for problems with a few thousands of data points.Our method fully takes advantages of recent advances in IAPG, saving us from developing the corresponding convergence results.We are also able to develop theoretical optimal estimates of the model parameters.This gives a good indication how we should set the parameter values in our implementation.Numerical results both on social networks and the benchmark test problems in manifold learning show that our method can fast produce embeddings of high quality.

Organization and Notation
The paper is organized as follows.Section 2 provides necessary background with a purpose to cast the MVU and MVE models as EDM-score driven models.This viewpoint will greatly benefit us in understanding our model, which is described in Section 3 with more detailed interpretation.We report our error bound results in Section 4, where only the main results are listed with all the technical proofs being moved to Appendix.Section 5 contains an inexact accelerated proximal gradient method as well as the theoretical optimal estimates of the model parameters.We report our extensive numerical experiments in Section 6 and conclude the paper in Section 7.
Notation.Let S n be the real vector space of n × n real symmetric matrices with the trace inner product X, Y := trace(XY ) for X, Y ∈ S n and its induced Frobenius norm • .Denote S n + the symmetric positive semidefinite matrix cone.We use I ∈ S n to represent the n by n identity matrix and 1 ∈ n to represent the vector of all ones.Let e i ∈ n , i = 1, . . ., n be the vector with the i-th entry being one and the others being zero.For a given X ∈ S n , we let diag(X) ∈ n denote the vector formed from the diagonal of X.
Below are some other notations to be used in this paper: • For any Z ∈ m×n , we denote by Z ij the (i, j)-th entry of Z.We use O n to denote the set of all n by n orthogonal matrices.
• For any Z ∈ m×n , we use z j to represent the j-th column of Z, j = 1, . . ., n.Let J ⊆ {1, . . ., n} be an index set.We use Z J to denote the sub-matrix of Z obtained by removing all the columns of Z not in J .
• Let I ⊆ {1, . . ., m} and J ⊆ {1, . . ., n} be two index sets.For any Z ∈ m×n , we use Z IJ to denote the |I| × |J | sub-matrix of Z obtained by removing all the rows of Z not in I and all the columns of Z not in J .
• We use " • " to denote the Hardamard product between matrices, i.e., for any two matrices X and Y in m×n the (i, j)-th entry of • For any Z ∈ m×n , let Z 2 be the spectral norm of Z, i.e., the largest single value of Z, and Z * be the nuclear norm of Z, i.e., the sum of single values of Z.The infinity norm of Z is denoted by Z ∞ .

Background
This section contains three short parts.We first give a brief review of cMDS, only summarizing some of the key results that we are going to use.We then describe the MVU and MVE models, which are closely related to ours.Finally, we explain three most commonly used distance-sampling rules.

cMDS
cMDS has been well documented in Cox and Cox (2000); Borg and Groenen (2007).In particular, Section 3 of Pȩkalska et al. (2002) explains when it works.Below we only summarize its key results for our future use.A n × n matrix D is called Euclidean distance matrix (EDM) if there exist points p 1 , . . ., p n in r such that D ij = p i − p j 2 for i, j = 1, . . ., n, where r is called the embedding space and r is the embedding dimension when it is the smallest such r.
An alternative definition of EDM that does not involve any embedding points {p i } can be described as follows.Let S n h be the hollow subspace of S n , i.e., Define the almost positive semidefinite cone K n + by where , 1935;Young and Householder, 1938) that D ∈ S n is EDM if and only if Moreover, the embedding dimension is determined by the rank of the doubly centered matrix JDJ, i.e., r = rank(JDJ) and where J is known as the centering matrix.Since −JDJ is positive semidefinite, its spectral decomposition can be written as where P T P = I and λ 1 ≥ λ 2 ≥ • • • ≥ λ n ≥ 0 are the eigenvalues in the nonincreasing order.Since rank(JDJ) = r, we must have λ i = 0 for all i ≥ (r + 1).Let P 1 be the submatrix consisting of the first r columns (eigenvectors) in P .One set of the embedding points are cMDS is built upon the above result.Suppose a pre-distance matrix D (i.e., D ∈ S n h and D ≥ 0) is known.It computes the embedding points by (2).Empirical evidences have shown that if the first r eigenvalues are positive and the absolute values of the remaining eigenvalues (they may be negative as D may not be a true EDM) are small, then cMDS often works well.Otherwise, it may produce mis-leading embedding points.For example, there are examples that show that ISOMAP might cut off too many eigenvalues, hence failing to produce satisfactory embedding (see e.g., Teapots data example in Weinberger and Saul ( 2006)).Both MVU and MVE models aim to avoid such situation.
The EDM score has been widely used to interpret the percentage of the total variance being explained by the embedding from leading eigenvalues.The EDM score of the leading k eigenvalues is defined by It is only well defined when D is a true EDM.The justification of using EDM scores is deeply rooted in the classic work of Gower (1966), who showed that cMDS is a method of principal component analysis, but working with EDMs instead of correlation matrices.
The centering matrix J plays an important role in our analysis.It is the orthogonal projection onto the subspace 1 ⊥ and hence J 2 = J.Moreover, we have the following.Let S n c be the geometric center subspace in S n : Let P S n c (X) denote the orthogonal projection onto S n c .Then we have P S n c (X) = JXJ.That is, the doubly centered matrix JXJ, when viewed as a linear transformation of X, is the orthogonal projection of X onto S n c .Therefore, we have It is also easy to verify the following result.
Lemma 1 For any X ∈ S n h , we have

MVU and MVE Models
The input of MVU and MVE models is a set of partially observed distances Let {p i } n i=1 denote the desired embedding points in r .They should have the following properties.The pairwise distances should be faithful to the observed ones.That is, and those points should be geometrically centered in order to remove the translational degree of freedom from the embedding: Let K := n i=1 p i p T i be the Gram matrix of the embedding points.Then the conditions in ( 5) and ( 6) are translated to To encourage the dimension reduction, MVU argues that the variance, which is Tr(K), should be maximized.In summary, the slack model (or the least square penalty model) of MVU takes the following form: where ν > 0 is the penalty parameter that balances the trade-off between maximizing variance and preserving the observed distances.See also Weinberger et al. (2007); Sun et al. (2006) for more variants of this problem.The resulting EDM D ∈ S n from the optimal solution of ( 7) is defined to be and it satisfies K = −0.5JDJ.Empirical evidence shows that the EDM scores of the first few leading eigenvalues of K are often large enough to explain high percentage of the total variance.MVE seeks to improve the EDM scores in a more aggressive way.Suppose the targeted embedding dimension is r.MVE tries to maximize the eigen gap between the leading r eigenvalues of K and the remaining eigenvalues.This gives rise to max There are a few standard ways in dealing with the constraints corresponding to (i, j) ∈ Ω 0 .We are interested in the MVE slack model: where ν > 0. The MVE model ( 8) often yields higher EDM scores than the MVU model ( 7).However, ( 7) is a SDP problem while ( 8) is nonconvex, which can be solved by a sequential SDP method (see Shaw and Jebara, 2007).

Distance Sampling Rules
In this part, we describe how the observed distances indexed by Ω 0 are selected in practice.We assume that those distances are sampled from unknown true Euclidean distances d ij in the following fashion.
where ξ ij are i.i.d.noise variables with E(ξ) = 0, E(ξ 2 ) = 1 and η > 0 is a noise magnitude control factor.We note that in (9) it is the true Euclidean distance d ij (rather than its squared quantity) that is being sampled.There are three commonly used rules to select Ω 0 .
(i) Uniform sampling rule.The elements are independently and identically sampled from Ω with the common probability 1/|Ω|.
(iii) Unit ball rule.For a given radius > 0, (i, j) ∈ Ω 0 if and only if d ij ≤ .
The k-NN and the unit ball rules are often used in low-dimensional manifold learning in order to preserve the local structure of the embedding points, while the uniform sampling rule is often employed in some other dimensionality reductions including embedding social network in a low-dimensional space.

A Convex Optimization Model for Distance Learning
Both MVU and MVE are trusted distance learning models in the following sense.They both produce a Euclidean distance matrix, which is faithful to the observed distances and they both encourage high EDM scores from the first few leading eigenvalues.However, it still remains a difficult (theoretical) task to quantify how good the resulting embedding is.In this part, we will propose a new learning model, which inherit the good properties of MVU and MVE.Moreover, we are able to quantify by deriving error bounds of the resulting solutions under the uniform sampling rule.Below, we first describe our model, followed by detailed interpretation.

Model Description
In order to facilitate the description of our model and to set the platform for our subsequent analysis, we write the sampling model ( 9) as an observation model.Define two matrices D and D (1/2) respectively by A sampled basis matrix X has the following form: for some (i, j) ∈ Ω.
For each (i, j) ∈ Ω 0 , there exists a corresponding sampling basis matrix.We number them as X 1 , . . ., X m .Define the corresponding observation operator O : That is, O(A) samples all the elements A ij specified by (i, j) ∈ Ω 0 .Let O * : m → S n be its adjoint, i.e., Thus, the sampling model ( 9) can be re-written as the following compact form where y = (y 1 , . . ., y m ) T and ξ = (ξ 1 , . . ., ξ m ) T are the observation vector and the noise vector, respectively.Since −JDJ ∈ S n + , we may assume that −JDJ ∈ S n + has the following single values decomposition (SVD): where P ∈ O n is an orthogonal matrix, λ = (λ 1 , λ 2 , . . ., λ n ) T ∈ n is the vector of the eigenvalues of −JDJ arranged in nondecreasing order, i.e., Suppose that D is a given initial estimator of the unknown matrix D, and it has the following single value decomposition −J DJ = P Diag( λ) P T , where P ∈ O n .In this paper, we always assume the embedding dimension r := rank(JDJ) ≥ 1.Thus, for any given orthogonal matrix P ∈ O n , we write P = [P 1 P 2 ] with P 1 ∈ n×r and P 2 ∈ n×(n−r) .For the given parameters ρ 1 > 0 and ρ 2 ≥ 0, we consider the following convex optimization problem This problem has EDM as its variable and this is in contrast to MVU, MVE and other learning models (e.g., Javanmard and Montanari, 2013) where they all use SDPs.The use of EDMs greatly benefit us in deriving the error bounds in the next section.Our model (13) tries to accomplish three tasks as we explain below.

Model Interpretation
The three tasks that model ( 13) tries to accomplish correspond to the three terms in the objective function.The first (quadratic) term is nothing but corresponding to the quadratic terms in the slack models ( 7) and ( 8).Minimizing this term is essentially to find an EDM D that minimizes the error rising from the sampling model (11).
The second term I, −JDJ is actually the nuclear norm of (−JDJ).Recall that in cMDS, the embedding points in (2) come from the spectral decomposition of (−JDJ).Minimizing this term means to find the smallest embedding dimension.However, as argued in both MVU and MVE models, minimizing the nuclear norm is against the principal idea of maximizing variance.Therefore, to alleviate this confliction, we need the third term − P 1 P T 1 , −JDJ .
In order to motivate the third term, let us consider an extreme case.Suppose the initial EDM D is close enough to D in the sense that the leading eigenspaces respectively spanned by { P 1 } and by {P 1 } coincide.That is P 1 = P 1 .Then, Hence, minimizing the third term is essentially maximizing the leading eigenvalues of (−JDJ).Over the optimization process, the third term is likely to push the quantity up, and the second term (nuclear norm) forces the remaining eigenvalues down.The consequence is that the EDM score EDMscore(r) = f (t, s) := t t + s gets higher.This is because ∀ t 2 > t 1 and s 2 < s 1 .
Therefore, the EDM scores can be controlled by controlling the penalty parameters ρ 1 and ρ 2 .The above heuristic observation is in agreement with our extensive numerical experiments.Model (13) also covers the MVE model as a special case.Let ρ 2 = 2 and D to be one of the iterates in the MVE SDP subproblems.The combined term is just the objective function in the MVE SDP subproblem.In other words, MVE keeps updating D by solving the SDP subproblems.Therefore, our error-bound results also justify why MVE often leads to higher EDM scores.Before we go on to derive our promised error-bound results, we summarize the key points for our model (13).It is EDM based rather than SDP based as in the most existing research.The use of EDM enables us to establish the error-bound results in the next section.It inherits the nice properties in MVU and MVE models.We will also show that this model can be efficiently solved.

Error Bounds Under Uniform Sampling Rule
Suppose that X 1 , . . ., X m are m independent and identically distributed (i.i.d.) random matrices over Ω with the common1 probability 1/|Ω|, i.e., for any 1 ≤ i < j ≤ n, Thus, for any A ∈ S n h , we have Moreover, we assume that the i.i.d.noise variables in ( 9) have the bounded fourth moment, i.e., there exists a constant γ > 0 such that E(ξ 4 ) ≤ γ.
Let D be the unknown true EDM.Suppose that the positive semidefinite matrix −JDJ has the singular value decomposition (12) and P = [P 1 , P 2 ] with P 1 ∈ n×r .We define the generalized geometric center subspace in S n by (compare to (3)) Let T ⊥ be its orthogonal subspace.The orthogonal projections to the two subspaces can hence be calculated respectively by It is clear that we have the following orthogonal decomposition A = P T (A) + P T ⊥ (A) and P T (A), P T ⊥ (B) = 0 ∀ A, B ∈ S n . (15) Moreover, we know from the definition of P T that for any A ∈ S n , which implies that rank(P T ⊥ (A)) ≤ 2r.This yields for any A ∈ S n For given ρ 2 ≥ 0, define Let ζ := (ζ 1 , . . ., ζ m ) T be the random vector defined by The non-commutative Bernstein inequality provides the probability bounds of the difference between the sum of independent random matrices and its mean under the spectral norm (see, e.g., Recht, 2011;Tropp, 2012;Gross, 2011).The following Bernstein inequality is taken from (Negahban and Wainwright, 2012, Lemma 7), where the independent random matrices are bounded under the spectral norm or bounded under the ψ 1 Orlicz norm of random variables, i.e., Lemma 2 Let Z 1 , . . ., Z m ∈ S n be independent random symmetric matrices with mean zero.Suppose that there exists M > 0, for all l, Z l 2 ≤ M or Z l 2 ψ 1 ≤ M .Denote σ 2 := E(Z 2 l ) 2 .Then, we have for any t > 0, Now we are ready to study the error bounds of the model ( 13).Denote the optimal solution of ( 13 The second major technical result below shows that the sampling operator O satisfies the following restricted strong convexity (Negahban and Wainwright, 2012) in the set C(τ ) for any τ > 0, where Lemma 4 Let τ > 0 be given.Suppose that m > C 1 n log(2n), where C 1 > 1 is a constant.Then, there exists a constant C 2 > 0 such that for any A ∈ C(τ ), the following inequality holds with probability at least 1 − 1/n.
Next, combining Proposition 3 and Lemma 4 leads to the following result.
Proposition 5 Assume that there exists a constant b > 0 such that D ∞ ≤ b.Let κ > 1 be given.Suppose that ρ 1 ≥ κη 1 m O * (ζ) 2 and ρ 2 ≥ 0. Furthermore, assume that m > C 1 n log(2n) for some constant C 1 > 1.Then, there exists a constant C 3 > 0 such that with probability at least 1 − 1/n, This bound depends on the model parameters ρ 1 and ρ 2 .In order to establish an explicit error bound, we need to estimate ρ 1 (ρ 2 will be estimated later), which depends on the quantity  18).To this end, from now on, we always assume that the i.i.d.random noises ξ l , l = 1, . . ., m in the sampling model ( 11) satisfy the following sub-Gaussian tail condition.
Assumption 1 There exist positive constants K 1 and K 2 such that for all t > 0, By applying the Bernstein inequality (Lemma 2), we have Proposition 6 Let ζ = (ζ 1 , . . ., ζ m ) T be the random vector defined in (18).Assume that the noise magnitude control factor satisfies η < ω := O(D (1/2) ) ∞ .Suppose that there exists C 1 > 1 such that m > C 1 n log(n).Then, there exists a constant C 3 > 0 such that with probability at least 1 − 1/n, This result suggests that ρ 1 can take the particular value: where κ > 1.Our final step is to combine Proposition 5 and Proposition 6 to get the following error bound.
Theorem 7 Suppose that there exists a constant b > 0 such that D ∞ ≤ b, and the noise magnitude control factor satisfies η < ω = O(D (1/2) ) ∞ .Assume the sample size m satisfies m > C 1 n log(2n) for some constant C 1 > 1.For any given κ > 1, let ρ 1 be given by ( 22) and ρ 2 ≥ 0.Then, there exists a constant C 4 > 0 such that with probability at least 1 − 2/n, The only remaining unknown parameter in ( 23) is ρ 2 though α(ρ 2 ).It follows from ( 17) that (α(ρ 2 ))2 = 1 2r Since P 1 P T 1 2 = P 1 P T 1 2 = r and P 1 P T 1 , P 1 P T 1 ≥ 0, we can bound α(ρ 2 ) by This bound seems to suggest that ρ 2 = 0 (corresponding to the nuclear norm minimization) would lead to a lower bound than other choices.In fact, there are better choices.The best choice ρ * 2 for ρ 2 is when it minimizes the right-hand side bound and is given by ( 25) in Subsection 5.1, where we will show that ρ 2 = 1 is a better choice than both ρ 2 = 0 and ρ 2 = 2 (see Proposition 8).
The major message from Theorem 7 is as follows.We know that if the true Euclidean distance matrix D is bounded, and the noises are small (less than the true distances), in order to control the estimation error, we only need samples with the size m of the order r(n − 1) log(2n)/2, since |Ω| = n(n − 1)/2.Note that, r = rank(JDJ) is usually small (2 or 3).Therefore, the sample size m is much smaller than n(n − 1)/2, the total number of the off-diagonal entries.Moreover, since the degree 2 of freedom of n by n symmetric hollow matrix with rank r is n(r − 1) − r(r − 1)/2, the sample size m is close to the degree of freedom if the matrix size n is large enough.

Model Parameter Estimation and the Algorithm
In general, the choice of model parameters can be tailored to a particular application.A very useful property about our model ( 13) is that we can derive a theoretical estimate, which serves as a guideline for the choice of the model parameters in our implementation.In particular, We set ρ 1 by ( 22) and prove that ρ 2 = 1 is a better choice than both the case ρ 2 = 0 (correponding to the nuclear norm minimization) and ρ 2 = 2 (MVE model).The first part of this section is to study the optimal choice of ρ 2 and the second part proposes an inexact accelerated proximal gradient method (IAPG) that is particularly suitable to our model.

Optimal Estimate of ρ 2
It is easy to see from the inequality (23) that in order to reduce the estimation error, the best choice ρ * 2 of ρ 2 is the minimum of α(ρ 2 ).We obtain from (24) that ρ * 2 ≥ 0 and The key technique that we are going to use to estimate ρ * 2 is the Löwner operator.We express both the terms P 1 P T 1 and P 1 P T 1 as the values from the operator.We then show that the Löwner operator admits a first-order apprximation, which will indicate the magnitude of ρ * 2 .The technique is extensively used by Miao et al. (2012).We briefly describe it below.
Immediately we have Φ(−JDJ) = P 1 P T 1 .We show it is also true for D. It follows the perturbation result of Weyl for eigenvalues of symmetric matrices (Bhatia, 1997, p. 63) that We must have λ i ≥ λ r − δ for i = 1, . . ., r and λ i ≤ δ for i = r + 1, . . ., n.
We therefore have Φ(−J DJ) = P 1 P T 1 .As a matter of fact, the scalar function defined by ( 26) is twice continuously differentiable (actually, φ is analytic) on (−∞, δ) ∪ (λ r − δ, ∞).Therefore, we know from (Bhatia, 1997, Exercise V.3.9) that Φ is twice continuously differentiable near −JDJ (actually, Φ is analytic near −JDJ).Therefore, under the condition that δ < λ r /2, we have by the derivative formula of the Löwner operator (see, e.g., Bhatia, 1997, Theorem V.3.3) that where H := D − D and W ∈ S n is given by We note that the leading r × r block of W is 0, which implies Therefore, we know from (25) that if D is sufficiently close to D, ).This shows that ρ 2 = 1 is nearly optimal if the initial estimator D is close to D. We will show that in terms of the estimation errors the choice ρ 2 = 1 is always better than the nuclear norm penalized least squares model (ρ 2 = 0) and the minimum volume embedding model (ρ 2 = 2).

Inexact Accelerated Proximal Gradient Method
Without loss of generality, we consider the following convex quadratic problem min where K n + is the almost positive semidefinite cone defined by ( 1), X, C ∈ S n , a ∈ m , b ∈ k , and A : S n → m , B : S n → k are two given linear operators.By setting 1 )J, one can easily verify that ( 28) is equivalent with the trusted distance learning model ( 13).
Being a convex quadratic problem, ( 28) can be solved in a number of ways.Based on our extensive numerical experiments, we found that the inexact accelerated proximal gradient (IAPG) method, which is recently studied by Jiang et al. (2012) for large scale linearly constrained convex quadratic SDP problems, is particularly suitable to our problem.To describe this method, let us denote the objective function by f (X) := 1 2 A(X)−a 2 + C, X and the corresponding gradient by ∇f (X) = A * (A(X) − a) + C, where A * : m → S n is the adjoint of A. The algorithm is described as follows.
Algorithm 1 Choose the starting point Y 1 = X 0 .Let t 1 = 1.Set k = 1.Perform the k-th iteration as follows: Step 1 Find an approximate minimizer where for all X ∈ S n , and Q k is a self-adjoint positive definite linear operator.
Step 2 Compute Step 3 It is noted that the major computational part is on the approximate solution of (29).There are three facts that make this computation very efficient.The first is that we can cheaply choose Q k ≡ I, the identical mapping on S n .This is because For different applications, one may want to choose different forms of Q k (see Jiang et al., 2012, for more details).The second fact is that the resulting problem of (29) with Q k ≡ I takes the following form: The geometric meaning of problem ( 30) is that it is to compute the nearest EDM from the matrix (Y k − ∇f (Y k )).This type of problems can be efficiently solved by the semismooth Newton method developed in Qi ( 2013).The third fact is that problem ( 30) is only solved approximately to meet a sufficient accuracy criterion (see (Jiang et al., 2012, eq. ( 32)) for the detailed formulations).This leads to a significant speed up of the already very fast semismooth Newton method.We note that the inexact APG inherits the computational 1/ √ ε complexity of the exact APG such as the FISTA developed by Beck and Teboulle (2009).Other complexity results on our inexact APG algorithm are similar as those of (Jiang et al., 2012, Theorem 3.1).We omit the details here to save space.

Numerical Experiments
In this section, we demonstrate the effectiveness of the proposed EDM Embedding (EDME) model ( 13) by testing Algorithm 1 on some real world examples.The examples are in two categories: one is of the social network visualization problem, whose initial link observation can be modelled by uniform random graphs.The other is from manifold learning, whose initial distances are obtained by the k-NN rule.The known physical features of those problems enable us to evaluate how good EDME is when compared to other models such as ISOMAP and MVU.It appears that EDME is capable of generating configurations of very high quality both in terms of extracting those physical features and of higher EDM scores.The test also raises an open question whether our theoretical results can be extended to this case.
For comparison purpose, we also report the performance of MVU and ISOMAP for most cases.The SDP solver used is the state-of-art SDPT3 package, which allows us to test problems of large data sets.We did not compare with MVE as it solves a sequence of SDPs and consequently it is too slow for our tested problems.Details on this and other implementation issues can be found in Subsection 6.3.

Social Network
Four real-world networks arising from the different applications are used to demonstrate the quality of our new estimator from EDME.
(SN1) Communication networks.Two sets of data are used in this experiment: the Enron email dataset (n = 183 users, Cohen ( 2009)) and the facebook-like social network (n = 1893 users, Opsahl and Panzarasa (2009)).For each dataset, we count the observed communications between the i-th and j-th users during a fixed time period from the mail/message logs, and denote this quantity by C ij .The social distances (or dissimilarities) between users are computed from the communication counts.It is natural to assume that larger communication count implies smaller social distance.Without loss of generality, we employ the widely used Jaccard dissimilarity (Klavans and Boyack, 2006) to measure the social distance of users: Thus, the observed measurement matrix D is incomplete, i.e., only a few entrances of the social distance matrix D are observed (the Enron email network: < 13% and the Facebook-like social network: < 0.8%).The corresponding two dimensional embedding results obtained by the MVU and our EDME model are shown in Figure 1 & 2 respectively.It can be seen from the presented eigenvalue spectrums3 of the Gramm matrices K = −1/2JDJ in Figure 1 & 2 that our EDME method is able to more accurately represent the high dimensional social network data in two dimensions than the MVU for both examples, since we capture all variance of the data in the top two eigenvectors, i.e., the EDMscore(2) = 100%.In fact, the MVU actually disperses the data into a much higher dimensional space (the rank of the Gramm matrix is much higher than the desired dimension) and only a small percentage of the variance can be explained from the top two eigenvectors.The details on the numerical performance on the MVU and EDME on these two examples will be reported in Table 1.(SN3) US airport network.In this example, we try to visualize the social network of the US airport network from the data of 2010 Opsahl (2011).There are n = 1572 airports under consideration.The number of the passengers transported from the i-th airport to the j-th airport in 2010 is recorded and denoted by C ij .Therefore, the social distance between two cities can be measured by the passenger numbers.For simplicity, we use the Jaccard dissimilarity defined in (31) to compute the corresponding distances between two cities.The observed distance matrix is also incomplete, and only very few entrances are observed (< 1.4%).The two dimensional embeddings obtained by the MVU and EDME methods are shown in Figure 4.The ten busiest US airports by total passenger traffic in 2010 are indicated by the red circles.Note that there are a large number of passengers transporting between them, which means the corresponding social distances among them should be relatively small.Thus, it is reasonable to expect that the embedding points of these top ten airports cluster around the zero point.Both MVU and EDME methods are able to show this important feature.However, it can be seen from the eigenvalue spectrums in Figure 4 that the MVU only captured 74.3% variance in the top two leading eigenvectors, while the EDME method captured all the variance in the two dimensional space.In this paper, we will use the data on the links between the blogs, which can be found from Freeman (2010) to visualize the corresponding social network.Similar to the communication network, we use the widely used Jaccard dissimilarity defined in (31) to measure the social distance of blogs.Without loss of generality, the 718 isolated blogs are removed from the original data, which means that we consider the remaining n = 1222 blogs with 586 left-leanings and 636 right-leanings.The social networks obtained by the MVU and the EDME are presented in Figure 5. From the results, we can see clearly that the embedding points generated by the MVU are concentrated near the zero point, and the rank of the corresponding Gram matrix is much higher than 2, which is 1135.However, our EDME method is able to capture all variance of the data in the two dimensions, providing a more accurate lower dimensional embedding.In fact, the embedding points in the visualizing network obtained by the EDME are naturally separated into two groups: the left-leaning blogs (the blue circles) and the right-leaning ones (the red circles).

Manifold learning
In this subsection, we test 4 widely used data sets in manifold learning.The initial distances used are generated by the k-NN rule.We describe them below with our findings.
(ML1) Teapots data In this example, we consider the two dimensional embedding of the images of a teapot (Weinberger and Saul, 2006).The teapots are rotated 360 degrees, and n = 400 images are taken from different angles.Each image has 76 × 101 pixels with 3 byte  In this example, we try to represent the high dimensional face image data (Tenenbaum et al., 2000) in a low dimension space.There are n = 698 images (64 pixel by 64 pixel) of faces with the different poses (up-down and left-right) and different light directions.Therefore, it is natural to expect that these high dimensional input data lie in the three dimensional space parameterized by the face poses and the light directions and that the equal importance of the three features can be sufficiently captured.Similar to the previous example, we use k = 5 to generate a connected graph.Both MVU and EDME methods successfully represent the data in the desired three dimensional space and their embedding results of the MVU and EDME are similar.For simplicity only the result of the EDME is shown in Figure 7.However, the Gram matrix learned by the ISOMAP has more than three nonzero eigenvalues.This is shown in the corresponding eigenvalue spectrums in Figure 7. Furthermore, for the ISOMAP, if we only compute the two-dimension embedding, then we only capture a smaller percentage of the total variance.It is interesting to observe that EDME is the only model that treats the three features equally important (the three leading eigenvalues are roughly equal).Moreovre, the EDME model performs much better than MVU in terms of the numerical efficiency.See Table 1 for more details.(ML3) The digits Data The data is from the MNIST database (LeCun et al., 1998).
We first consider the data set of digit "1", which includes n = 1135 8-bit grayscale images of "1".Each image has 28 × 28 pixels, which is represented as 784 dimensional vector.We note that the two most important features of "1"s are the slant and the line thickness.Therefore, the embedding results are naturally expected to lie in the two dimensional space parameterized by these two major features.In this example, we set k = 6. Figure 8 shows the two dimensional embeddings computed by ISOMAP, MVU and EDME.It can be clearly seen that EDME significantly outperforms the other two methods.In particular, EDME is able to accurately represent the data in the two dimensional space and captures the correct features.However, MVU returns an almost one dimensional embedding and only captures one of the major features, i.e., the slant of "1"s.For the ISOMAP, it only captures a small percentage of the total variance.Next, we consider the learning task of the mixed digits "1" and "9".We randomly choose 500 images of "1"s and 500 images of "9"s from the MNIST.For the mixed images, the major features clearly become the size of the top arches and the slant of the digits.However, both ISOMAP and MVU fail to obtain a two dimensional embeddings, which are clearly indicated in the eigenvalue spectrum in Figure 9.Moreover, two dimensional projections (along the first two eigenvectors) of ISOMAP and MVU embeddings can not capture the desired two major features of images.In contrast, EDME successfully represents the data in the two dimensional space with the correct features (the top arch size and the slant of digits).(ML4) Data of Frey face Finally, we consider the comparison of MVU and EDME on the Frey face images data (Roweis and Saul, 2000), which has n = 1965 images of faces.Each image has 28 × 20 gray scale pixels and is represented by a vector of 560 dimensions.We set k = 4.The representative faces are shown at some randomly chosen embedding points.From the eigenvalue spectrums shown in Figure 10, MVU returns a three dimensional embedding, while EDME is able to obtain a two dimensional embedding.Moreover, the two dimensional projection of the MVU embedding along the first two major eigenvectors only contain 86.2% variance of the data, which is smaller than that from EDME.

Numerical performance
We tested the ISOMAP, the MVU and our proposed EDME methods in MATLAB R2014a (version 8.3.0.532), and the numerical experiments are run in MATLAB under a MAC OS X 10.9.3 64-bit system on an Intel 4 Cores i5 2.7GHz CPU with 8GB memory.
In our numerical experiments, we use the SDPT3 (Toh et al., 1999), a Matlab software package for semidefinite-quadratic-linear programming, to solve the corresponding SDP problem of the MVU model, since the SDPT3 is more efficient than other SDP solves such as CSDP (Borchers, 1999) andSeDuMi (Sturm, 1999).This is particularly the case for the large dimensional SDP problems (e.g., the Facebook-like communication network and the US airport).The termination tolerance of the SDPT3 is tol = 10 −3 (tol = 10 −6 for the Enron email network and the Madrid train bombing).For our EDME model, the proposed inexact accelerated proximal gradient (IAPG) method (Algorithm 1) is employed.We terminate the algorithm if the primal and dual infeasibilities conditions are met, i.e., max{R p , R d } ≤ tol, where R p = B(X)−b and R d = (∇f (X)−B * y−Z)/(1+ C ) are the measurements of the infeasibilities of the primal and dual problems, respectively, where the tolerance is chosen by tol = 10 −3 (tol = 10 −6 for the Enron email network and the Madrid train bombing).The details on the numerical performance of the MVU and EDME methods can be found from Table 1, where we report the EDM scores from the leading two eigenvalues and cpu time in seconds.We observe that the performance of EDME is outstanding in terms of numerical efficiency.The developed IAPG is much faster than the SDP solver.Taking USairport2010 as example, MVU used about 12 hours while EDME only used about 80 seconds.For the  2013) that the shortest path distance is not suitable to measure the distances in social networks.We also like to point out that for all tested problems, EDME captured nearly 100% variance and it treats the local features equally important in terms of the leading eigenvalues being of the same magnitude.

Conclusions
The paper aimed to explain a mysterious situation regarding the SDP methodology to reconstruct faithful Euclidean distances in a low-dimensional space from incomplete set of noisy distances.The SDP models can construct numerical configurations of high quality, but they lack theoretical backups in terms of bounding errors.We took a completely different approach that heavily makes use of Euclidean Distance Matrix instead of positive semidefinite matrix in SDP models.This led to a convex optimization that inherits the nice features of MVU and MVE models.More importantly, we were able to derive errorbound results under the uniform sampling rule.The optimization problem can also be efficiently solved by the proposed algorithm.Numerical results in both social networks and manifold leading showed that our model can capture low-dimensional features and treats them equally important.
Given that our model worked very well for the manifold learning examples, an interesting question regarding this approach is whether the theoretical error-bound results can be extended to the case where the distances are obtained by the k-NN rule.It seems very difficult if we follow the technical proofs in this paper.It also seems that the approach of Javanmard and Montanari (2013) would lead to some interesting (but very technical) results.We plan to investigate those issues in future.
For the first term of the right hand side of (33), we have By noting that D * , D ∈ S n H , we know from Lemma 1 that the rank of D * − D − J(D * − D)J is no more than 2, which implies Moreover, it follows from (4) that J(D * − D)J, D * − D − J(D * − D)J = 0, which implies Thus, we have By noting that P Meanwhile, since for any A ∈ S n , P T (A) * = P T 2 AP 2 * , we know from the directional derivative formulate of the nuclear norm (Watson, 1992, Theorem 1) that By using the decomposition (15) and the notations defined in (17), we conclude from (35) that Thus, together with (37), we know from (33) that ρ 1 κ and κ > 1, we know from ( 16) and ( 35) that Since r ≥ 1, the desired inequality (19) follows from (40), directly.
Next we shall show that (20) also holds.By (39), we have Therefore, by combining with (36) and ( 16), we know from the decomposition (15) that Finally, since r ≥ 1, we conclude that This completes the proof.
Thus, it follows from Massart's concentration inequality (see, e.g., Bühlmann and Van De Geer, 2011, Theorem 14.2) that By applying the standard Rademacher symmetrization (see, e.g., Koltchinskii, 2011, Theorem 2.1), we obtain that where {ε 1 , . . ., ε m } is an i.i.d.Rademacher sequence.Again, since A ∞ = 1/ √ 2, we know that | A, X i | ≤ A ∞ < 1.Thus, it follows from the contraction inequality (see, e.g., Ledoux and Talagrand, 1991, Theorem 4.12) that For any A ∈ C(τ ; Υ), we have Thus, we obtain that It follows from (41) that By choosing Υ = 2 k ν, we obtain from the above inequality that By noting that log(x) < x for any x > 1, we conclude that P(B) ≤ Finally, the lemma then follows if we prove that for m > C 1 n log n with C 1 > 1, there exists a constant C 1 > 0 such that The following proof is similar with that of Lemma 7 Klopp (2011) (see, e.g., Klopp, 2014, Lemma 6).We include it again for seeking of completion.Denote Z l := ε l X l , l = 1, . . ., m.
By applying the Bernstein inequality (Lemma 2), we obtain the following tail bound for any t > 0, Let C 1 = e 1/2 2 1/ log 2 .Then, we know from (45) that the inequality (42) holds.Thus, we know from ( 14) and ( 19) in Proposition 3 that By substituting τ , we obtain that there exists a constant C 3 > 0 such that nm .
The result then follows by combining these two cases.
) by D * .The following result represents the first major step to derive our ultimate bound result.It contains two bounds.The first bound (19) is on the norm-squared distance betwen D * and D under the observation operator O.The second bound (20) is about the nuclear norm of D * − D. Both bounds are in terms of the Frobenius norm of D * − D. Proposition 3 Let ζ = (ζ 1 , . . ., ζ m ) T be the random vector defined in (18) and κ > 1 be given.Suppose that ρ 1 ≥ κη 1 m O * (ζ) 2 and ρ 2 ≥ 0, where O * is the adjoint operator of O.Then, we have 1 m O * (ζ) 2 , where ζ = (ζ 1 , . . ., ζ m ) T ∈ m with ζ l , l = 1, . . ., m are i.i.d.random variables given by (

Figure 1 :
Figure 1: The embedding networks of Enron email network.

Figure 3 :
Figure 3: The embedding networks of Madrid train bombing.

Figure 4 :
Figure 4: The embedding networks of USairport2010

Table 1 :
Numerical performance comparison of the MVU and the EDME