How can classical multidimensional scaling go wrong?
Year: 2022

Abstract:
Given a matrix D describing the pairwise dissimilarities of a data set, a common task is to embed the data points into Euclidean space. The classical multidimensional scaling (cMDS) algorithm is a widespread method to do this. However, theoretical analysis of the robustness of the algorithm and an in-depth analysis of its performance on non-Euclidean metrics is lacking. In this paper, we derive a formula, based on the eigenvalues of a matrix obtained from D, for the Frobenius norm of the difference between D and the metric Dcmds returned by cMDS. This error analysis leads us to the conclusion that when the derived matrix has a significant number of negative eigenvalues, then kD DcmdskF , after initially decreasing, will eventually increase as we increase the dimension. Hence, counterintuitively, the quality of the embedding degrades as we increase the dimension. We empirically verify that the Frobenius norm increases as we increase the dimension for a variety of non-Euclidean metrics. We also show on several benchmark datasets that this degradation in the embedding results in the classification accuracy of both simple (e.g., 1-nearest neighbor) and complex (e.g., multi-layer neural nets) classifiers decreasing as we increase the embedding dimension. Finally, our analysis leads us to a new efficiently computable algorithm that returns a matrix Dl that is at least as close to the original distances as Dt (the Euclidean metric closest in `2 distance). While Dl is not metric, when given as input to cMDS instead of D, it empirically results in solutions whose distance to D does not increase when we increase the dimension and the classification accuracy degrades less than the cMDS solution.

heading: 1 Introduction
text: Multidimensional scaling (MDS) refers to a class of techniques for embedding data into Euclidean space given pairwise dissimilarities [Carroll and Arabie, 1998, Borg and Groenen, 2005]. Apart from the general usefulness of dimensionality reduction, MDS has been used in a wide variety of applications including data visualization, data preprocessing, network analysis, bioinformatics, and data exploration. Due to its long history and being well studied, MDS has many variations such as non-metric MDS [Shepard, 1962a,b], multi-way MDS [Kroonenberg, 2008], multi-view MDS [Bai et al., 2017], confirmatory or constrained MDS [Heiser and Meulman, 1983], etc. (See France and Carroll [2010], Cox and Cox [2008] for surveys).
The basic MDS formulation involves minimizing an objective function over a space of embeddings. There are two main objective functions associated with MDS: STRESS and STRAIN.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
The STRAIN objective (Equation 1 below) was introduced by Torgerson [1952], whose algorithm to solve for this objective is now commonly referred to as the classical MDS algorithm (cMDS).
Xcmds := argmin X2Rr⇥n
X T X
✓ V DV
2
◆ . (1)
Here V is the centering matrix given by V := I 1nJ, and I is the identity matrix and J is the matrix of all ones. cMDS first centers the squares of the given distance matrix and then uses its spectral decomposition to extract the low dimensional embedding. cMDS is one of the oldest and most popular methods for MDS, and its popularity is in part due to the fact that this decomposition is fast and can scale to large matrices. The point set produced by cMDS, however, is not necessarily the point set whose Euclidean distance matrix is closest under say Frobenius norm to the input dissimilarity matrix. This type of objective is instead captured by STRESS, which comes in a variety of related forms. In particular, in this paper we consider the SSTRESS objective (see Equation 2).
Specifically, given an embedding X , let EDM(X) be the corresponding Euclidean distance matrix, that is EDM(X)ij = kXi Xjk2F , where Xi, Xj are the ith and jth columns of X . If D is a dissimilarity matrix whose entries are squared, then we are interested in the matrix,
Dt := argmin D0=EDM(X),X2Rr⇥n
kD0 Dk2F . (2)
Note that the reason we assume our dissimilarity matrix has squared entries is because the standard EDM characterizations uses squared entries (see further discussion of EDMs below). Equation 2 is a well studied objective [Takane et al., 1977, Hayden and Wells, 1988, Qi and Yuan, 2014].
There are a number of similarly defined objectives. If one considers this objective when the matrix entries are not squared (i.e. q D0ij p Dij), then it is referred to as STRESS. If one further normalizes each entry of the matrix difference by the input distance value (i.e. ( q D0ij p Dij)/ p Dij) then it
is called Sammon Stress. In this paper, we are less concerned with the differences between different types of Stress, and instead focus on how the cMDS solution behaves generally under a Stress type objective. Thus for simplicity we focus on SSTRESS. It is important to note that there are algorithms to solve the SSTRESS objective, but the main drawback is that they are slow in comparison to cMDS [Takane et al., 1977, Hayden and Wells, 1988, Qi and Yuan, 2014]. Thus, many practitioners default to using cMDS and do not optimize for SSTRESS.
In this paper, we shed light on the theoretical and practical differences between optimizing for these two objectives. Let Dcmds := EDM(Xcmds), where Xcmds is the solution to Equation 1, and let Dt be the solution to Equation 2. We are interested in understanding the quantity
err := kD Dcmdsk2F . (3)
Doing so will provide practitioners with multiple advantages and will guide the development of better algorithms. In particular,
1. Understanding err is the first step in rigorously quantifying the robustness of the cMDS algorithm. 2. If err is guaranteed to be small, then we can use the cMDS algorithm without having to worry about loss in quality of the solution. 3. If err is big, we can make an informed decision about the benefits of the speed of the cMDS algorithm versus the quality of the solution. 4. Understanding when err is big helps us design algorithms to approximate Dt that perform better when cMDS fails.
Contributions. Our main theorem, Theorem 1, decomposes err into three components. This decomposition gives insight into when and why cMDS can fail with respect to the SSTRESS objective. In particular, for Euclidean inputs, err naturally decreases as the embedding dimension increases. For non-Euclidean inputs, however, our decomposition shows that after an initial decrease, counterintuitively err can actually increase as the embedding dimension increases. In practice one may not know a priori what dimension to embed into, though one might assume it suffices to embed into some sufficiently large dimension. Importantly, these results demonstrate that when using cMDS to embed, choosing a dimension too large can actually increase error.
This degradation of the cMDS solution is of particular concern in relation to the robustness in the presence of noisy or missing data, as may often be the case for real world data. Several authors [Cayton and Dasgupta, 2006, Mandanas and Kotropoulos, 2016, Forero and Giannakis, 2012] have proposed variations to specifically address robustness with cMDS. However, our decomposition of err, suggests a novel approach. Specifically, by attempting to directly correct for the problematic term in our decomposition (which resulted in err increasing with dimension) we produce a new lower bound solution. We show empirically that this lower bound corrects for err increasing, both by itself and when used as a seed for cMDS. Crucially the running time of our new approach is comparable to cMDS, rather than the prohibitively expensive optimal SSTRES solution. Finally, and perhaps more importantly, we show that if we add noise or missing entries to real world data sets, then our new solution outperforms cMDS in terms of the downstream task of classification accuracy, under various classifiers. Moreover, our decomposition can be used to quickly predict the dimension where the increase in err might occur.
The main contributions of our paper are as follows.
1. We decompose the error in Equation 3 into three terms that depend on the eigenvalues of a matrix obtained from D. Using this analysis, we show that there is a term that tells us that as we increase the dimension that we embed into, eventually, the error starts increasing. 2. We verify, using classification as a downstream task, that this increase in the error for cMDS results in the degradation of the quality of the embedding, as demonstrated by the classification accuracy decreasing. 3. Using this analysis, we provide an efficiently computable algorithm that returns a matrix Dl such that if Dt is the solution to Equation 2, then kDl DkF  kDt DkF , and empirically we see that kDl DtkF  kDcmds DtkF . 4. While Dl is not metric, when given as input to cMDS instead of D, it results in solutions that are empirically better than the cMDS solution. In particular, this modified procedure results in a more natural decreasing of the error as the dimension increases and has better classification accuracy.

heading: 2 Preliminaries and Background
text: In this section, we lay out the preliminary definitions and necessary key structural characterizations of Euclidean distance matrices.
cMDS algorithm. For completeness, we include the classical multidimensional scaling algorithm in Algorithm 1.
Algorithm 1 Classical Multidimensional Scaling. 1: function CMDS(D, r) 2: X = V ⇤D ⇤ V/2 3: Compute µ1 . . . µr > 0, U as the eigenvalues and eigenvectors of X 4: return U ⇤ diag(pµ1, . . . , p µr, 0, . . . , 0)
EDM Matrices
Definition 1. D 2 Rn⇥n is a Euclidean Distance Matrix (EDM) if and only if there exists a d 2 N such that there are points x1, . . . , xn 2 Rd with
Dij = kxi xjk22.
Note that unlike other distance matrices, an EDM consists of the squares of the (Euclidean) distances
between the points.
This form permits several important structural characterizations of the cone of EDM matrices.
1. Gower [1985], Schoenberg [1935] show that a symmetric matrix D is an EDM if only if
F := (I 1sT )D(I s1T )
is a positive semi-definite matrix for all s such that 1T s = 1 and Ds 6= 0.
2. Schoenberg [1938] showed that D is an EDM if and only if exp( D) is a PSD matrix with 1s along the diagonal for all > 0. Note here exp is element wise exponentiation of the matrix.
3. Another characterization is given by Hayden and Wells [1988] in which D is an EDM if and only if D is symmetric, has 0s on the diagonal (i.e., the matrix is hollow), and D̂ is negative semi-definite, where D̂ is defined as follows
QDQ =
 D̂ f
f T ⇠
. (4)
Here f is a vector and
Q = I 2 vT v vv T for v = [1, . . . , 1, 1 + p n]T . (5)
Note Q is Householder reflector matrix (in particular, it’s unitary) and it reflects a vector about the hyperplane span(v)?.
The characterization from Hayden and Wells [1988] is the main characterization that we shall use. Hence we establish some important notation.
Definition 2. Given any symmetric matrix A 2 Rn⇥n, let us define Â 2 Rn 1⇥n 1, f(A) 2 Rn 1, and ⇠(A) 2 R as follows
QAQ =
 Â f(A)
f(A)T ⇠(A)
.
In addition to characterizations of the EDM cone, we are also interested in the dimension of the EDM. Definition 3. Given an EDM D, the dimensionality of D is the smallest dimension d, such that there exist points x1, . . . , xn 2 Rd with Dij = kxi xjk22. Let E(r) be the set of EDM matrices whose dimensionality is at most r.
Using Hayden and Wells [1988]’s characterization of EDMs, Qi and Yuan [2014] show D 2 E(r) if and only if D is symmetric, hollow (i.e., 0s along the main diagonal), and D̂ in Equation 4 is negative semi-definite with rank at most r.
Conjugation matrices: Q and V . Conjugating distance matrices either by Q (in Equation 5) or by the centering matrix V is an important component of understanding both EDMs and the cMDS algorithm. We observe that V is also (essentially) a Householder matrix like Q. Let w = [1, . . . , 1]T and observe that J = wwT so that V = I 1wTwww T .
Qi and Yuan [2014] establishes an important connection between Q and V that we make use of in our analysis in Section 3. Specifically, for any symmetric matrix A, we have that
V AV = Q  Â 0 0 0 Q. (6)
Here Â is the matrix given by Definition 2. This connection gives a new framework in which we can understand the cMDS algorithm. We know that given D, cMDS first computes V DV/2. Using Equation 6, this is equal to Q  D̂ 0 0 0 Q/2. Thus, when cMDS computes the spectral
decomposition of V DV/2, this is equivalent of computing the spectral decomposition of D̂. Then using the characterization of E(r) from Qi and Yuan [2014], setting all but the largest r eigenvalues to 0 might seem like the optimal solution. However, this procedure ignores condition that the matrix must be hollow. As we shall show, this results in sub-optimal solutions.

heading: 3 Theoretical Results
text: Throughout this section we fix the following notation. Let D be a distance matrix with squared entries. Let r be the dimension into which we are embedding. Let 1  . . .  n 1 be the eigenvalues of D̂ and U be the eigenvectors. Let be an n-dimensional vector where i = i i>0 or i>r for
i = 1, . . . , n 1 and n = 0. Let S = Q ⇤  U 0 0 1 . Let Dt be the solution to Problem 2 and Dcmds
the resulting EDM from the solution to Problem 1. Let
C1 = n 1X
i=1
2i , C2 = n 1X
i=1
i, C3 = nk(S S) k2F C22
2 .
Then the main result of the paper is the following spectral decomposition of the SSTRESS error. Theorem 1. If D is a symmetric, hollow matrix then, kDcmds Dk2F = C1 + C22 + C3.
The idea behind the proof of Theorem 1 is to decompose kDcmds Dk2F = kQDcmdsQ QDQk2F
= 4kD̂/2 D̂cmds/2k2F + (⇠(D) ⇠(Dcmds))2 + 2kf(D) f(Dcmds)k2F , which follows from Definition 2. We relate each of these three terms to C1, C2, and C3. The following lemmas will work towards expressing each of these terms separately. In the following discussion, let Yr := Xcmds T Xcmds and recall that Xcmds is the solution to the classical MDS problem given in Equation 1. The proofs for the following lemmas are in the appendix. Lemma 1. If G is a positive semi-definite Gram matrix, then 12V EDM(G)V = Q  Ĝ 0 0 0 Q
Lemma 2. The value of the objective function obtained by Xcmds in Equation 1 is C1 4 . Specifically, we have that 4kYr ( V DV )/2k2F = n 1X
i=1
2i =: C1.
Lemma 3. 12D̂cmds = Ŷr.
Lemma 4. If Tr(D) = 0, then (⇠(D) ⇠(Dcmds))2 =
n 1X
i=1
i !2 =: C22 .
Lemma 5. If D is hollow, then 2kf(D) f(Dcmds)k2F = nk(S S) k2F C22
2 =: C3.
The first term in the error is C1. From the definition, we can see C1 is the sum of the squares of the eigenvalues of D̂ corresponding to eigenvectors that are not used (or are discarded) in the cMDS embedding. As we increase the embedding dimension, we use more eigenvectors. Hence as the embedding dimension increases, we see that C1 monotonically decreases. In the case when D is an EDM, we can use all of the eigenvectors so this term will go to zero. On the other hand, if D is not an EDM and D̂ has positive eigenvalues (i.e. negative eigenvalues of V DV ), then these are eigenvalues that correspond to eigenvectors that cannot be used. Thus, cMDS will always exhibit this phenomenon for C1 regardless of the input D (for both STRAIN and SSTRESS).
The second term is C22 . This term looks similar to C1, but instead of summing the eigenvalues squared, we first sum the eigenvalues and then take the square. This subtle difference has a big impact. First, we note that as r increases, C2 becomes more negative. Suppose that D is not an EDM, (i.e., D̂ has positive eigenvalues) and let K be the number of negative eigenvalues of D̂. Then, since D has trace 0, when r = K, C2 is negative. Hence, C2 decreases and is eventually negative as r increases which implies that as r increases, there is a a certain value of r after which C22 increases. This results in the quality of the embedding decreasing. As we will see, this term will be the
dominant term in the total error.
While C1 and C2 are simple to understand, C3 is more obtuse. To simplify it, we consider the following. If is an entry of a random n by n unitary matrix, then as n goes to infinity the distribution for p n converges to N (0, 1) and the total variation between the two distributions is bounded above by 8/(n 2) [Easton, 1989, Diaconis and Freedman, 1987]. Therefore, we can assume that the variance of an entry of a random n by n orthogonal matrix is about 1/n. So, heuristically, S S ⇡ 1n11 T and
nk(S S)~ k2F ⇡ n n2 k11T k2F = n n2 k1C2k2F = C22 .
Then since C3 = nk(S S) k2F C 2 2 2 ⇡ C22 C 2 2
2 , we see that, at least heuristically, the overall behavior of C3 is dominated by C2.
Algorithm 2 Lower Bound Algorithm. 1: function LOWER(D, r) 2: Compute D̂, f(D) and ⇠(D). 3: Compute 1  . . .  n 1, U as the eigenvalues and eigenvectors of D̂ 4: Initialize ci = i and cn = ⇠(D) 5: Let negative_C2 = 0. 6: for i = 1 . . . n 1 do 7: if i > 0 or i > r then 8: negative_C2 += i 9: Set ci to 0 10: E := number of c’s not equal to 0 11: sub := negative_C2/E 12: for i = 1 . . . n 1 do 13: if cn i! = 0 then 14: if cn i + sub  0 then 15: cn i += sub, E = 1, negative_C2 = sub 16: else 17: E = 1, negative_C2 += cn 1, sub = negative_C2/E 18: cn 1 = 0 19: cn += sub, negative_C2 = sub 20: Let D̂ = U ⇤ diag(c1, . . . , cn 1) ⇤ UT
21: return Q ⇤ 
D̂ f(D) f(D)T cn ⇤Q
From the previous discussion, we see that the C22 term in the decomposition is the most vexing due to the term (⇠(D) ⇠(Dcmds))2. This term is from the excess in the trace; that is, the result of discarding the eigenvalues changes the value of the trace from 0 to non-zero. From Lemma 3, we see that cMDS projects D̂/2 onto the cone of PSD matrices (which do not necessarily have trace 0). To retain the trace 0 condition, we perform the following adaptation. Let (r) be the space of symmetric, trace 0 matrices A, such that Â is negative semi-definite of rank at most r and project D onto (r) instead. That is, we seek a solution the following problem.
Dl := argmin A2(r)
kA Dk2F . (7)
Clearly this problem is a strict relaxation of the SSTRESS MDS problem and hence we get that
kD DtkF kDl DkF Because we no longer require the solution to be hollow, conjugating by S permits us to rewrite Problem 7 as:
c := argmin c2Rn,cT 1=0, 8n>i, ci0,
8n>j>r, cj=0
n 1X
i=1
(ci i)2 + (cn ⇠(D))2 (8)
Theorem 2. If c is the solution to Problem 8 and Dl is the solution to problem 7 and if we let M be the n 1 by n 1 diagonal matrix with the first n 1 terms of c on the diagonal, then
S T DlS =
 M f(D)
f(D)T cn
.
Unlike Problem 2, Problem 8 can be solved directly. We do so via Algorithm 2. We see that in the first loop (lines 6-9), we set ci to be 0 if i > 0 or i > r to ensure the proper constraints on cis (i.e., the eigenvalues). Doing so, we incur a cost of C1. That is, we sum the squares of those eigenvalues. Now the solution after line 9 does not necessarily satisfy cT 1 = 0. In fact, at this stage of the algorithm, we have that cT 1 = C2, and we have E many c’s that are non zero that we can
adjust. That is, we modify these entries on average by C2/E. Because we use a Frobenius norm to measure error, this incurs a cost of C22E/E2 = C22/E. Finally, we note that the major computational step of the algorithm is the spectral decomposition. Hence it has a comparable run time to cMDS. Theorem 3. If D is a symmetric, hollow matrix, then Algorithm 2 computes Dl the solution to 7 in O(n3) time.
Proposition 1. If D is a symmetric, hollow matrix, then
kDl Dk2F C1 + C
2 2
r + 1
It is important to note that Dl is not a metric. However, if we want an EDM and hence an embedding, we can use Dl as the input to cMDS algorithm. Noting that D̂l is negative semi-definite matrix of rank at most r, we see that if Dlcmds is the metric returned by cMDS on input Dl, then we have that kDl Dlcmdsk2F = 2kf(Dl) f(Dlcmds)k2F = 2kf(D) f(Dlcmds)k2F . Similar to Lemma 5, we get the following result.
Lemma 6. 2kf(D) f(Dlcmds)k2F = nk(S S)ck2F
C22 r+1
2 =: C4.
Using Lemma 6, we see that kDlcmds Dk2F  kDlcmds Dlk2F + kDl Dk2F
 2kf(Dlcmds) f(Dl)k2F + kDt Dk2F .
) kDlcmds Dk2F kDt Dk2F  nk(S S)ck2F
C22 r+1
2 = C4
On the other hand, we upper bound kDcmds Dk2F kDt Dk2F using Proposition 1. Taking the difference, we get
kDcmds Dk2F kDt Dk2F  r
r + 1 C
2 2 + C3
Using our heuristics for C3 and C4, we expect Dcmds to be be much further from Dt, than Dlcmds. These claims are experimentally verified in the next section.
Solving Equation 2. Another approach would be to solve Equation 2 directly. However, solving this problem is much more computationally expensive and works such as Qi and Yuan [2014] present new algorithms to do so. Algorithm 2 can be used to compute the solution as well. Here we Dykstra’s method of alternating projections as Hayden and Wells [1988] do, but instead of only projecting D̂ on the cone of negative semi definite matrices, without any rank constraint, we project onto (r) using Algorithm 2. The second, projection is then onto the the space of hollow matrices. Alternating these projections gives us an iterative method to compute Dt.

heading: 4 Experiments
text: In this section, we do two things. First, we empirically verify all of the theoretical claims. Second, we show that on the downstream task of classification, if we use cMDS to embed the data, then as the
embedding dimension increases, the classification accuracy gets worse. Thus, suggesting that the embedding quality of cMDs degrades.
Verifying theoretical claims. To do this, we look at three different types of metrics. First, are metrics that come from graphs and for these we use Celegans Rossi and Ahmed [2015] and Portugal Rozemberczki et al. [2019] datasets. Second, are metrics obtained in an intermediate step of Isomap. Third, are perturbations of Euclidean metrics. For both of these metrics we use the heart dataset Detrano et al. [1989]. In all cases, we show that the classical MDS algorithm does not perform as previously expected; instead, it matches our new error analysis. In each case, we demonstrate that our new algorithm Lower + cMDS outperforms cMDS in relation to SSTRESS.
Results. First, let us see what happens when we perturb an EDM. Here we consider two different measures. Let D be the original input, and Dp be the perturbed. Then we measure
kDp Dcmdsk2F kDpk2F and kD Dcmdsk2F kDk2F .
As we can see from Figure 2c and 2d, as the embedding dimension increases both quantities eventually increase. The increase in the relative SSTRESS is as we predicted theoretically. The increase in the second quantity suggests that cMDS does not do a good job of denoising either. This suggests that the cMDS objective is not a good one.
Next, we see how cMDS does on graph datasets and with Isomap on Euclidean data. Here we plot the relative squared error (kD Dcmdsk2F /kDk2F ). Figure 2 shows that, as predicted, as the embedding dimension increases, the cMDS error eventually increases as well. For both the Portugal dataset alone and with Isomap on the heart dataset, this error eventually becomes worse than the error when embedding into two dimensions! As we can see from Figure 1, this increase is exactly due to the C22 term. Also, as heuristically predicted, the C3 term is roughly constant.
Let us see how the true SSTRESS solution and new approximation algorithm perform. We look at the relative squared error again. First, Figure 2 shows us that our lower bound tracks the true SSTRESS solution extremely closely. We can also see from Figure 2d that the true SSTRESS solution and the Lower + cMDS solution are closer to the original EDM compared to the perturbed matrix. Thus, they are better at denoising than cMDS. Finally, we see that our new algorithm Lower + cMDS performs better than just cMDS and, in most cases, fixes the issue of the SSTRESS error increasing with dimension. We also see that kDlcmds Dk2F kDt Dk2F is also roughly constant.
We point out cMDS and Lower+cMDS have comparable running times as they can respectively be computed with one or two spectral decompositions. Computing Dt requires computing a spectral decomposition in every round, with larger datasets requiring hundreds of rounds till they appear to converge. This is a significant advantage of Lower+cMDS over True, and also the reason we were not able to compute True for the Portugal dataset above, and the classifications experiments below.
Classification. For the classification tasks, we switch to more standard benchmark datatsets: MNIST, Fashion MNIST, and CIFAR10. If we treat each pixel as a coordinate then these datasets are Euclidean and we cannot demonstrate the issue with cMDS. We obtain non-Euclidean metrics in three ways. First, we compute the Euclidean metric and then perturb it with a symmetric, hollow matrix whose off diagonal entries are drawn from a Gaussian. Second, we construct a k nearest neighbor graph and then compute the all pair shortest path metric on this graph. This is the metric that Isomap tries to embed using cMDS. Third, we imagine each image has missing pixels and then compute the distance between any pair of images using the pixels that they have in common (as done in Balzano et al. [2010], Gilbert and Sonthalia [2018]). Thus, we have 9 different metrics to embed. For each dataset, we constructed all 3 metrics for the first 2000 images. We embedded each of the nine metrics into Euclidean space of dimensions from 1 to 1000 using cMDS and our Lower + cMDS algorithm. We do not embed by solving Equation 2 as this is computationally expensive. For each embedding, we then took the first 1000 images are training points and trained a feed-forward 3 layer neural network to do classification. Results with a nearest neighbor classifier are in the appendix. We then tested the network on the remaining 1000 points. We can see from Figure 3, that as the embedding dimension increases, the classification accuracy drops significantly when trained on the cMDS embeddings. On the other hand, when trained on the Lower + cMDS embeddings, the classification accuracy does not drop, or if it degrades, it degrades significantly less. Thus, showing that the increase in SSTRESS for cMDS on non Euclidean metrics, does result in a degradation of the quality of the embedding and that our new method fixes this issues to a large extent.

heading: Acknowledgements
text: Work on this paper by Gregory Van Buskirk and Benjamin Raichel was partially supported by NSF CAREER Award 1750780.

Review summaries:
The work shows a scenario that the classical MDS algorithm fails. "When the derived matrix has a significant number of negative eigenvalues," the classical MDS algorithm tends to fail when the number of dimensions increases sufficiently large. Then the work designs an algorithm with improved performance bound and shows the potential on practical applications.

The authors study the problem of multidimensional scaling (MDS) for embedding data into Euclidean space and perform a theoretical error analysis of the classical muldimensional scaling (cMDS) algorithm. Their main result decomposes the error of cMDS into 3 components, which they use to explain the behaviour and degradation of the cMDS error with respect to the embedding dimension. Then, they propose an algorithm, based on their theoretical findings, and show that it outperforms cMDS on standard image classification datasets.

A new decomposition of the classical multidimensional scaling (cMDS) error term is generated and used to get additional insights of the method.  As I first came across the draft, I naively thought; “Well, cMDS is extremely well known, how can one innovate on such topic?”. Moreover, the departing justification stating cMDS does not excel in non-linear mappings wasn’t that great; I mean cMDS ASSUMES Euclidean space and linearity so it was no wonder that used out of its natural scope it may be suboptimal. But I am very glad to say that the authors proved my naivety wrong very quickly! There is a very clear contribution, the paper is strongly supported, and reading enjoyable.  STRENGTHS + The unusual “late comer”. Innovating in a method which is so well known is difficult.  WEAKNESSES + There are some unclear decisions and justifications but it is somewhat contradictory that I’m criticizing this when, I in the place of the authors, I would very likely have proceed equally… Notwithstanding, I’m indicating this below, in case a second thought can be given + There is an unconventional mixture of theoretical claims with experimental testing whereby some “theoretical proofs” are complemented with heuristics verified experimentally, but actually not theoretically proven. 

The paper considers the classic embedding method known as multidimensional scaling and focus on the so called STRAIN  and SSTRESS objectives, both designed for embedding a given squared distance matrix into a lower dimensional space, often for visualization, but potentially also for reducing problem size for downstream tasks.   The idea is that the SSTRESS is the better objective to optimize as it captures how well the distances of the embedded data points approximate the input distances, but computationally hard, so maybe it is enough to just use the solution for the simpler STRAIN objective, that boils down to Singular Value Decomposition and that may work fine.  The paper analyzes the solution of STRAIN in terms of the more relevant SSTRESS objective, and shows a decomposition of the error of the STRAIN solution as measured with the STRESS objective and analyzes this objective to understand what happens. The punchline is that if the input distance matrix is not Euclidian (no data set with these pairwise distances exist in any dimension), then the SSTRESS error of the embedding achieved by optimizing STRAIN can  actually increase with increased embedding dimension. A property that seems counterintuitive as a higher dimension should allow improved embeddings. The paper shows  experimentally that this phenomena actually occurs in practice as well and that the performance is nicely explained by the error decomposition proved and discussed in the paper. In particular the experiments show that the classic Isomap embedding has this issue.  Furthermore, if these STRAIN embeddings are used for a simple downstream task like standard classification, the results become increasingly worse as the embedding size is increased.