Many datasets in applications adopt dissimilarities that are non-Euclidean and/or nonmetric. Examples of such popular dissimilarity measures [19,7] include Minkowski distance (L p ), cosine similarity, Hamming, Jaccard, Mahalanobis, Chebyshev, and KL-divergence. Studies in psychology have long recognized that human perception of similarity is not a metric [42]. Further, dissimilarity matrices with negative entries (e.g.: cosine similarity, correlation, signed distance) have also been widely used in various problems. Negative inner product norms also have deep connections to hyperbolic spaces as well as the study of spacetime in special relativity theory.
In many machine learning practices, embedding in low dimensional vector space is often explicitly or implicitly done within the data processing pipeline. Such embedding may already need to consider more general dissimilarities. For example, one branch of graph learning adopts embedding in non-Euclidean spaces (e.g., hyperbolic spaces [9]), and machine learning for physics model and data (AI4Sicence) needs to consider more general inner product norms [22]. In transformer models [45,13], the attention mechanism can also be viewed as learning a general bilinear form on tokens. Despite the wide adoption of general dissimilarity measures in practice, theoretical study of embedding and dimension reduction for non-Euclidean non-metric data appears to be still very limited.
In this paper we consider one of the most classical algorithms for data embedding and dimension reduction -multidimensional scaling -and develop a non-Euclidean, non-metric version with theoretical performance guarantee.
Background on MDS. Classical multidimensional scaling (cMDS) takes as input a Euclidean distance matrix (EDM), i.e., a symmetric matrix D 2 R n⇥n where each entry is the squared Euclidean distance between two points in Euclidean space, and recovers the Euclidean coordinates. Using a standard double centering trick, D can be turned into a ⇤ Rutgers University. {cd751,jg1555,kll160,fluo,hongbin.sun,cx122}@rutgers.edu Gram matrix B = X T X where X encodes the Euclidean coordinates. For the purpose of producing a low-dimensional vector, classical MDS takes k eigenvectors corresponding to top k largest eigenvalues of the Gram matrix B. This minimizes the strain, the di↵erence (Frobenius norm) in terms of the Gram matrix.
When the input distance matrix is not a Euclidean distance matrix, this problem is called metric MDS [1]. Metric MDS considers the minimization of STRESS [28], defined as the sum of squared di↵erence of pairwise embedding distances to the input dissimilarities. Minimizing STRESS makes the problem to be non-linear and there is no closed-form solution, though one can use either gradient descent or Newton's method [35,26]. Nevertheless, in practice, cMDS is often applied for non-Euclidean distance matrix. In this case, the centered matrix B is no longer positive semi-definite. The common practice is keep the top positive eigenvalues and throw away the negative eigenvalues.
Two recent papers [41,38] pointed out that classical MDS produces suboptimal solutions on non-Euclidean distance matrix when considering STRESS. This is not a surprise, since cMDS, minimizing strain, does not minimize STRESS. However, a more problematic issue is that when using more dimensions in cMDS (i.e., increasing k), the STRESS error first drops and then increases. We call this phenomena Dimensionality Paradox. It is theoretically unsatisfactory and counter-intuitive that embeddings by classical MDS using more dimensions could yield worse results.
The error analysis in [38] sheds light on this issue. When the input matrix is Non-Euclidean, the negative eigenvalues carry crucial information. cMDS, keeping only positive eigenvalues, is intrinsically biased -the more positive eigenvalues used the more it deviates from the input data. A real eradication of this issue must address the root cause, i.e., applying an algorithm meant for Euclidean geometry on non-Euclidean data.
Our Contributions. We extend multidimensional scaling to non-Euclidean geometry, by generalizing the dissimilarity function from the standard inner product (which defines Euclidean geometry) to the broader family of symmetric bilinear forms (u, v) = u T Av, where the symmetric matrix A does not have to be positive semi-definite. For dimension reduction, we look for both a low dimensional vector representation and an associated bilinear form, that together approximate the input dissimilarity matrix with minimum STRESS error. Specifically, the key contributions are as follows:
We conduct an in-depth analysis on STRESS error for any chosen subset of k < n eigenvalues of the input centered dissimilarity matrix. We propose Neuc-MDS, an e cient algorithm that finds the best subset of eigenvalues to minimize a lower bound of STRESS.
Beyond the constraints of eigenvalue subsets, we extend our findings to the general linear combinations of eigenvalues. Our advanced algorithm, Neuc-MDS + , finds the best linear combination of eigenvalues to minimize the lower bound objective.
We provide theoretical analysis for the asymptotic behavior of cMDS and Neuc-MDS on random symmetric matrices. First, both necessarily produce large STRESS if the target dimension k = o(n) -on completely unstructured data aggressive dimension reduction shall not be expected. Further, when k = ⇥(n), the STRESS of Neuc-MDS monotonically decreases to 0 while the STRESS of cMDS increases and eventually reaches a plateau.
Empirically we evaluate Neuc-MDS and Neuc-MDS + on ten diverse datasets encompassing di↵erent domains. The experiment results show that both methods substantially outperform previous baselines on STRESS and average distortion, and fully resolve the issue of dimensionality paradox in cMDS. Our codes are available on Githubfoot_0 .
Our work is in the general family of similarity learning [25,3] with dimension reduction, going beyond metric learning and embedding. Due to the huge amount of literature on this topic we only mention those that are most relevant.
Let P denote a dataset of size n. Let D 2 R n⇥n be the dissimilarity matrix of dataset P , D ij = D ji is a real-valued symmetric dissimilarity measure between pair p i , p j 2 P , i 6 = j 2 [n], and all diagonal entries D ii = 0 (i.e., D is a hollow matrix). D is the analog of squared Euclidean distance matrix in classical MDS. But here D is not necessarily Euclidean, may not be a metric (e.g. violating triangle inequality) and may have negative entries. Our goal is to obtain (1) a low-dimensional vector representation for each element in P , and (2) a function f that computes a dissimilarity measure using the calculated low dimensional vectors. Since the input dissimilarities are not necessarily Euclidean nor a metric, we look for the function f beyond Euclidean distances, but stay within a broader family of inner products or bilinear forms.
A bilinear form on a vector space V is a function : V ⇥ V ! R which is linear in each variable when the other variable is fixed. More precisely, (au + v, w) = a (u, w) + (v, w) and (w, au + v) = a (w, u) + (w, v) for all u, v, w 2 V and any scalar a. We only consider symmetric bilinear forms , i.e., (u, v) = (v, u). A bilinear form is positive definite (or positive semi-definite) if (u, u) > 0 for 8u 6 = 0 (or (u, u) 0). Symmetric matrices and symmetric bilinear forms are two sides of a coin. Namely, fix a basis W = {w 1 , ..., w n } of the vector space, there is a one-to-one correspondence between them. That is, A = [ (w i , w j )] n⇥n is a symmetric matrix. Conversely, give a symmetric matrix A, one defines a symmetric bilinear form (u, v) = u T W Av W where u W is the coordinate of vector u in the basis W .
Formally, we have the following problem definition.
Definition 1 (Non-Euclidean Dimension Reduction). Given a symmetric dissimilarity matrix D of a dataset P of size n and a natural number k n, find a collection of n k-
F for the dissimilarity matrix D of P given by Dij = f ( pi , pj ) is minimized.
A special case of a symmetric bilinear form is the standard inner product h, i on R n . Indeed the inner product hu, vi = u T v is a symmetric positive definite bilinear form. The inner product of u v and u v is precisely the squared Euclidean distance ||u v|| 2 . The Euclidean space is R n equipped with the standard inner product. Thus metric geometry of Euclidean space is governed by the inner product. On the other hand, a basic theorem of linear algebra states that a finite dimensional vector space equipped with a symmetric positive definite bilinear form is isometric to the Euclidean space of the same dimension.
Thus geometry of a positive definite symmetric bilinear form, or symmetric positive definite matrices, is nothing but Euclidean geometry.
When D represents the inner products of pairwise di↵erences (i.e., squared Euclidean distances) of n points P in R d , D is called a Euclidean distance matrix (EDM). It can be shown that by taking B = 1 2 CDC, where C = I 1 n 1 n 1 T n is the centering matrix and 1 n is a vector of ones, one obtains the Gram matrix B = X T X where X is d⇥n dimensional matrix of the n coordinates of dimension d. The matrix B is a symmetric positive semi-definite matrix. Thus using eigendecomposition of B one can recover the coordinates X. This procedure is classical multidimensional scaling (cMDS) [40]. Furthermore, if one would like to use k-dimension coordinates with k < d, classical MDS suggests to take the eigenvectors corresponding to the k largest positive eigenvalues of the Gram matrix B. This minimizes the strain, the di↵erence (Frobenius norm k • k F ) in terms of the Gram matrix.
X cmds , arg min
Now consider a general symmetric dissimilarity matrix D of size n⇥n, it naturally associates R n with a symmetric bilinear form . In many real-world situations, the dissimilarity matrix is not positive, nor negative definite, i.e., the bilinear form is indefinite. This means the intrinsic geometry (R n , ) is non-Euclidean. The relationship between the Gram matrix and square distance matrix still holds for indefinite bilinear forms. See Appendix A
In practice, classical MDS is often the default choice even when the input dissimilarity matrix D is not an EDM, i.e., the centered matrix B = 1 2 CDC is not positive semidefinite. Classical MDS, which simply drops negative eigenvalues of B to produce positive semi-definiteness does not respect the geometry well. Indeed, multiple researchers have observed suboptimal and counter-intuitive performance [41,38]. For example, increasing k may lead to increased STRESS error -keeping more dimensions makes the approximation worse! On a second thought, such results are not surprising. If the input data does not carry Euclidean geometry, forcing it through a procedure for Euclidean geometry is fundamentally problematic.
Our main observation is that a vector space with an indefinite symmetric bilinear form has its own intrinsic geometry. This geometry, even though may not be metrical, carries the most accurate information about the dissimilarity matrix and the datasets. Suppose the centered matrix B = 1 2 CDC has p > 0 positive eigenvalues and q > 0 negative eigenvalues. The associated indefinite bilinear form has signature (p, q). One such example is
The geometry of a finite dimensional vector space with a symmetric bilinear form of signature (p, q) where p, q > 0 is much less developed compared to Euclidean geometry. When q = 1, this is the Minkowski geometry and is closely related to relativity theory in physics and hyperbolic geometry [33] in mathematics. For general (p, q), one should probably abandon the notion of distance for indefinite spaces. Namely, the expression (u v, u v), even if it is positive, should not be considered as the square of the distance between two points u, v.
According to [33], one calls p (u v, u v) the Lorenzian distance between u, v. Despite the term distance in the name, the Lorenzian distance does not satisfy triangle inequality in general.
We propose Non-Euclidean MDS, a novel linear dimension reduction technique using bilinear forms. For a vector v, we use Diag(v) as the diagonal matrix with v along the main diagonal and zero everywhere else. For any given symmetric dissimilarity matrix D 2 R n⇥n , let the eigen decomposition of the centralization of D be given as follows:
where U 2 R n⇥n is the orthogonal matrix of eigenspace and ⇤ 2 R n⇥n is the diagonal matrix ⇤ =: Diag( ) with eigenvalues ,
By using an algorithm that will be discussed in Section 4 we will choose k of the eigenvalues, represented by a binary indicator vector w = (w 1 , • • • , w n ) 2 {0, 1} n with a value of 1 (or 0) indicating the corresponding eigenvalue is chosen (or not chosen). Let
Note that p ⇤ contains complex numbers for non-PSD matrix D. Since w has only k nonzero values, we can drop the n k zero rows in X (corresponding to the eigenvalues not selected) and have a k-dimensional vector representation of the data. Now we can derive dissimilarities by defining Di,j , (X i X j ) T (X i X j ), D = DIS(X) := Di,j .
See 1 for details.
In the description above, we allow the coordinates in X to take complex numbers, such that we can take Di,j to be the standard dot product of X i X j with X i X j . Alternatively, we can keep X to take real coordinates, i.e., X = p
Again we drop the n k zero rows in X and take a bilinear form f (u, v) = u T Av, where A is a k by k diagonal matrix, with the element at (i, i) to be 1 (or 1) if the corresponding eigenvalue chosen is positive/negative. Notice that the bilinear form takes precisely the format of (p, q)-distance as in Equation (2).
We have a few remarks in place: First we are not throwing away the negative eigenvalues. As will be explained later we actually keep eigenvalues of largest magnitude and some could be negative. As a consequence D may have negative real values, which is expected as we are moving away from Euclidean distances and input entries in D may even start to be negative. Second, when D is an EDM, i.e., all eigenvalues are non-negative, Neuc-MDS reduces to classical MDS. Last, similar to cMDS, our method also starts with computing the eigenvalues of the Gram matrix. For large datasets, fast (approximation) algorithm for partial SVD can also be applied to our methods. For example, on n ⇥ n symmetric matrices, using power methods one can iteratively compute partial SVD up to k largest/smallest eigenvalues in O(kn 2 ). With randomness introduced, it can be reduced to O(log k • n 2 ) [18,27]. For really large datasets, the dissimilarity matrix with size O(n 2 ) might already be too large to be acquired or stored, one may extend MDS through local embedding methods like Landmark MDS [12] or Local MDS [10]. This approach can also be applied to neuc-MDS to substantially speed up computation without su↵ering too much on performance (See Section 6 for empirical results).
To establish the foundation of theoretical analysis, we first analyze the STRESS error of Neuc-MDS and decompose it into three terms. Next we show an e cient algorithm that minimizes the first two terms that dominate. Last, we examine Neuc-MDS and classical MDS on random Gaussian matrices.
Inspired by the analysis of STRESS of classical MDS [38], we adopt a similar approach and decompose the STRESS error into three terms. Let 2 R n be the vector of all eigenvalues
w be the indicator vector of dropped eigenvalues, w 2 {0, 1} n , and be the Hadamard product. We have the following result with proof in Appendix B.
Note that C 1 /4 is the sum of squared eigenvalues that are dropped, and C 2 /4 is the square of the sum of eigenvalues that are dropped. Individually, we can minimize C 1 by keeping eigenvalues of large absolute value; and C 2 by balancing the dropped eigenvalues such that the summation has a small magnitude. For term C 3 , from Equation 21in Proof B, we know C 3 0. In [38] it is argued that if one takes the approximation
Although it is empirically observed in [38] that C 3 is roughly constant and hence negligible for optimization, there are some cases in our experiments in which C 3 is not negligible.
In light of Theorem 2, we would like to approximately optimizing the STRESS by minimizing the lower bound C 1 + C 2 , which can be formulated as a quadratic integer programming problem: Given a set of n values L = { i } and a positive integer k > 0, choose a k-subset S ✓ L such that min
When i 's are all positive, the best choice is to take the k largest eigenvalues of L, i.e., the cMDS' solution. However, with a mixture of positive and negative eigenvalues, taking the top k largest eigenvalues is no longer optimal -specifically, as k increases, the first error term is monotonically reduced but the second error term could start going up. In the following subsection, we discuss an optimal algorithm to solve Equation (5).
The optimization problem described in Equation ( 5) is a special case of the family of quadratic integer programming problems. Though in general, quadratic integer programming is NP-hard [31], we show that this particular one is actually solvable in polynomial time. Formally, we have: Theorem 3. For the optimization problem defined in (5), there exits an optimal solution with r largest positive values and s smallest negative values in L, r + s = k. And, there is an O(n)-algorithm that outputs an optimal solution.
Our algorithm ( 2) is greedy and iteratively selects the eigenvalue with the highest absolute value. Specifically, let S be the set of selected eigenvalues and H(S) be the sum of eigenvalues not selected. Initially S = ;. In each iteration, if H(S) < 0, select the negative eigenvalues remained of the greatest magnitude and add it to S; if H(S) > 0, pick the largest positive one. If H(S) = 0, pick the eigenvalue with the greatest magnitude. The complete proof of Theorem 3 is delayed to Appendix C.
Here, we analyze the important error term C 1 + C 2 in Theorem 2 for cMDS and Neuc-MDS when the input (centered) dissimilarity matrix is a symmetric random matrix. Our analysis is established upon Wigner's famous Semicircle Law [46,47], which states the following. Proposition 4 (Semicircle Law [47]). Suppose B 2 R n⇥n is a symmetric random matrix where every element is independently distributed with equal densities and second moments 2 . Let S a,b (B) be the number of eigenvalues of B that lie in the interval Our results are formally stated as follows with proof in Appendix D.
Theorem 5. Suppose a random symmetric matrix B 2 R n⇥n where B ij is independently distributed with equal densities and second moments 2 , is taken as the centered matrix 3 to classical MDS and Neuc-MDS, both selecting k eigenvalues with k n. Let e C denote the C 1 + C 2 error for cMDS and e N for Neuc-MDS, we have:
. Therefore the STRESS error of both classical MDS and Neuc-MDS cannot be very small if the target dimension k is o(n). Lastly, when k = cn, with c ! 1, e N is monotonically decreasing and eventually reaches 0. On the other hand, for cMDS, when c 1/2, e C reaches a plateau at about 0.1801 • n 3 2 . Notice that this error has an extra factor of n compared to the error for small k.
For real world data the Gram matrix is likely far from a random matrix. The analysis above points out that aggressive dimension reduction can be indeed only a luxury for structured data, even if we use inner products that are not limited to Euclidean distances. Second, any real world data carries some random measurement noise. When the scale of such random noise becomes non-negligible, STRESS error introduced by such noise cannot be small with aggressive dimension reduction. We would recommend practitioners to examine the spectrum of eigenvalues to gain insights on the power or limit in reducing dimensions.
Both cMDS and Neuc-MDS choose a subset of k eigenvalues from the input Gram matrix. This can be considered as applying Diag(w) to the eigenvalue vector ⇤ (Equation 4) to filter some eigenvalues out. Such operators can be viewed as a special case of more general low-rank linear maps T : R n ! R n with rank(T ) = k n.
Here we consider a new family of embeddings X = U ⇤1/2 with ⇤ being a rank-k diagonal matrix whose diagonal entries are given by some ˜ 2 R n with only k nonzero entries. Note that there always exists a rank-k linear operator T such that ˜ = T ( ). Then we ask if we can improve Neuc-MDS further with a rank-k linear map of ? This question can be answered by the following theorem: Theorem 6. Given a dissimilarity matrix D 2 R n⇥n with eigenvalues 2 R n of its Gram matrix CDC/2 = U ⇤U T , for any rank-k (k n) diagonal matrix ⇤ = Diag( ˜ ) with ˜ 2 R n with k nonzero entries only on coordinates given by some index k-set W ✓ [n], let w := 1 W 2 {0, 1} n be the indicator vector of W , and let D be the dissimilarity matrix reconstructed from X = ⇤1/2 U T . Then the STRESS error of D can be expressed as: 3 Here we take B as the centered matrix. Therefore the diagonal entries of B do not have to be 0. Further, the distribution of centered matrix (ignoring the scaling factor) CAC is the same as the distribution of (A) plus an additional zero eigenvalue. Therefore, one can consider (A) sampled from random matrices. When n ! 1, they are asymptotically the same. See more discussion in the appendix (Lemma 12).
where := ˜ . The first two terms, C1 + C2 , as a lower bound of the STRESS, is minimized as
Remark 7. Recall the decomposition of STRESS in Theorem 2, the lower bound 4 wT (2) + 4( wT
Therefore, it has a better lower bound of STRESS compared to Neuc-MDS.
Theorem 6 provides a constructive way to obtain the optimal ˜ that (approximately) minimizes the STRESS error for a prefixed w which is determined by an indicator set W served as the constraints of nonzero entries on ˜ . The next question is how to find an optimal k-set W on which the optimal ˜ ⇤ has the lowest lower bound. Essentially, we need to solve the following optimization problem:
Proposition 8. For the optimization problem (8), there exits an optimal solution with r largest positive values and s smallest negative values in L, r + s = k. And, there is an O(n)-algorithm that outputs the optimal solution.
The algorithm follows a similar greedy manner as EV-Selection, with only one adjustment: in each step, compare the two marginal gains provided by the largest positive eigenvalue and the lowest negative eigenvalue among the remained ones, and choose the positive/negative eigenvalue of largest magnitude if the corresponding gain is smaller. We name the algorithm Neuc-MDS + , and present the complete algorithm with its correctness proof in Appendix C.
This section presents experimental results of Neuc-MDS and Neuc-MDS + . First, we evaluate the performance on dissimilarity error of two proposed algorithms comparing with closelyrelated baselines on three metrics: STRESS, distortion and additive error. Then we show that the dimensionality paradox issue observed on cMDS is fully resolved by Neuc-MDS and Neuc-MDS + .
Synthetic Data. We introduce two synthetic datasets: Random-simplex and Euclideanball, both with non-Euclidean dissimilarities. See details in Appendix E.1. On a high level, suppose the dataset has size n, we construct a Random-simplex such that for each vertex, the first n 1 coordinates virtually form a simplex, while the last coordinate almost dominates the distances between the other points, which creates a large negative eigenvalue for the Gram matrix. The Euclidean-ball dataset (similar to Delft's balls [14]) considers n balls of di↵erent radii with the distance of two balls defined as the smallest distance of two points on the two respective balls. The dissimilarities by this definition no longer satisfy triangle inequality.
Real-world Data. We consider both genomics and image data. For genomics data, we include 5 datasets from the Curated Microarray Database (CuMiDa) [17], each indicating a certain type of cancer. Following the practice mentioned in [43], pairwise dissimilarities are generated with entropic a nity with the diagonal as zero. We also test three celebrated image datasets: MNIST, Fashion-MNIST and CIFAR-10. The dissimilarity matrix for each dataset captures 1000 images randomly sampled for each class. We use three measures in [38] to calculate the dissimilarities.
An essential guidance of our choice of datasets is the existence of substantial negative eigenvalues. Otherwise, Neuc-MDS and cMDS are equivalent. Figure 1 illustrates the eigenvalue distribution of the Renal dataset. Table 1 shows the basic statistics of each dataset and the number of negative eigenvalues of the Gram matrix. All datasets are non-Euclidean and non-metric. Baselines. We include cMDS, Lower-MDS [38] and SMACOF (Scaling by MAjorizing a COmplicated Function) [29] as baselines. Lower-MDS [38] looks for a symmetric, lowrank, trace-zero, and positive semi-definite matrix. SMACOF minimizes STRESS using majorization and is one of the best nonlinear optimization algorithms for MDS. All methods are deterministic therefore variance is not concerned. Performance on Dissimilarity Error. In addition to STRESS as the primary metric, we also test the average distortion (i.e. multiplicative error) and scaled additive error on all datasets. For all metrics, a smaller value indicates a more favorable performance. Limited by space, we leave scaled additive error together with other observations in Appendix E. With k as the target dimension, for synthetic datasets and images, we set k = 100, for genomics data, k = 20. The details are presented in Table 2.
The results 4 show that in terms of STRESS, Neuc-MDS and Neuc-MDS + outperform cMDS and Lower-MDS consistently by a large margin. SMACOF has comparable performance with Neuc-MDS in a couple data sets but can go out of bound in others. For average distortion, the genomics datasets di↵erentiate di↵erent methods drastically while the rest datasets produce comparable results. Neuc-MDS + occasionally gives slightly higher STRESS than Neuc-MDS. Recall that both methods focus on optimizing for a lower bound of STRESS (in Theorem 2 and Theorem 6). This shows that C 3 may play a role in practice.
Neuc-MDS Addresses 'Dimensionality Paradox'. Dimensionality paradox of classical MDS refers to the observation that STRESS increases as the dimension goes up. When raising this concern, [38] proposes a Lower-MDS algorithm as mitigation. We show that Neuc-MDS and Neuc-MDS + address this issue even better. For Lower-MDS the target dimension cannot be larger than the number of positive eigenvalues. Our methods do not have this limitation. Lower-MDS has a shorter curve because the target dimension k is limited to be smaller than the number of positive eigenvalues.
Landmark MDS Landmark MDS [12] is a heuristic to speed up classical MDS. One chooses a small number of landmarks and apply MDS on the landmarks first. The coordinates of the remaining points are obtained through a triangulation step with respect to the landmarks. We can use the same heuristic to speed up Neuc-MDS. On the random-simplex dataset, with only 25% points randomly chosen as landmarks the STRESS is only a factor of 1.0644 of the STRESS obtained by Neuc-MDS. If we use only 10% points as landmarks, the final STRESS is only 1.0898 of the STRESS of Neuc-MDS. This shows that Neuc-MDS can also be significantly accelerated using the landmark idea, achieving nearly the same STRESS.
This paper presents an extension of classical MDS to non-Euclidean non-metric settings. We would like to mention a few future directions. Since we step out of the domain of Euclidean embedding, both the input dissimilarity matrix and the one obtained after dimension reduction can have negative values. Therefore if one would like to feed the output dissimilarity matrix to another data processing module that by default requires non-negative values, special care must be taken to address the negative values. In experiments, we discovered that Neuc-MDS + achieves similar stress as Neuc-MDS and produces much fewer negatives values in the output dissimilarity matrix (Appendix E.3). How to e↵ectively use such dissimilarities in downstream learning and inference tasks would be a major future work. Note that the geometry of general bilinear forms is a largely unexplored territory.
We will prove the following proposition relating hollow symmetric matrices and squared "distance" matrices.
Proposition 9.
(1) Given any symmetric matrix G = [g ij ] n⇥n , there exists a symmetric bilinear form h, i on R n and n vectors v 1 , ..., v n such that g ij = hv i , v j i.
(2) Given any hollow symmetric matrix M = [m ij ] n⇥n , there exists a symmetric bilinear form h, i on R n and n vectors v 1 , ..., v n such that m ij = hv i v j , v i v j i for all i, j.
Note that the hollow condition m ii = 0 is necessary.
Proof. Part ( 1) is a well known fact from linear algebra. We take v 1 , ..., v n to be the standard basis of R n . From the symmetric matrix G, one defines the symmetric bilinear form by the formula hu, vi = u T Gv and vice versa.
To see part (2), consider the centered symmetric matrix B = CM C where C = I J n is the double centering matrix and J is the matrix whose entries are 1. One can see that its (i, j)-th entry b ij of matrix B is given by
One can easily verify from the formula that the centering condition holds:
By part (1) of the proposition, one finds n vectors v 1 , ..., v n in R n and a symmetric bilinear form h, i on R n such that 2b ij = hv i , v j i. Now we claim that
Indeed, by dropping the constant term 1
P n l=1 m kl in (9), we have
The last step uses the fact that M is a hollow matrix, i.e., m ii = m jj = 0.
In summary, we call G = [hv i , v j i] the Gram matrix and M = [hv i v j , v i v j i] and the squared distance matrix of a symmetric bilinear form h, i. Obviously, we can determine M from G by the formula
We can also recover the Gram matrix G from M from the double centering construction CM C if the center of the vectors v 1 , ..., v n is 0, i.e., P n i=1 v i = 0. To see this, given a squared distance matrix M associated to n vectors v 1 , .., v n , we can replace v i 's by v i w for a fixed vector w without changing M . Now by taking w = 1 n P n i=1 v i to be the center of v 1 , ..., v n , we can normalize P n i=1 v i = 0. In this case, formula (9) states that
where
In this section, we want to prove our main Theorem 6 for the error analysis of STRESS error. Note that Theorem 2 is a special case of Theorem 6 for w = 0, which is exactly the case in MDS and Neuc-MDS with ˜ given by w. Theorem 6. Given a dissimilarity matrix D 2 R n⇥n with eigenvalues 2 R n of its Gram matrix CDC/2 = U ⇤U T , for any rank-k (k n) diagonal matrix ⇤ = Diag( ˜ ) with ˜ 2 R n with k nonzero entries only on coordinates given by some index k-set W ✓ [n], let w := 1 W 2 {0, 1} n be the indicator vector of W , and let D be the dissimilarity matrix reconstructed from X = ⇤1/2 U T . Then the STRESS error of D can be expressed as:
where := ˜ . The first two terms, C1 + C2 , as a lower bound of the STRESS, is minimized as
The proof will essentially follow the pipeline in [38]. Our results and proofs can be viewed as an extension to a more general family of problems. We first show that the STRESS error can be decomposed into three terms C1 + C2 + C3 . Then we analyze the approximate lower bound C1 + C2 and figure out the minimum.
We first introduce some notations and lemmas we need to use later. Recall that C = I 1 n 1 n 1 T n is the centering matrix. First note the following fact: Lemma 10. The n-dimensional vector 1 n is an eigenvector with eigenvalue 0 of the matrices C, U and XT X.
Proof. By definition of C = I 1 n 1 n 1 T n , it is easy to check C1 n = 0. Also by definitions of U and XT X we can get similar results.
As a corollary, we have the following lemma: Let Q be a Householder reflector matrix defined as follows:
For any symmetric matrix A 2 R n⇥n , let (A) 2 R (n 1)⇥(n 1) , f(A) 2 R n 1 , ⇠(A) 2 R be given by QAQ through the following equation:
We need the following properties of Q: 32]). For a symmetric matrix A, we have
Lemma 13 ([21, 38]). For any matrix M , we use diag(M ) = (M i,i ) to denote the column vector of M 's diagonal entries. For a symmetric hollow matrix
we are ready to prove our main theorem.
Proof of Theorem 6.
For the first term:
For the second term:
For the third term:
Taking the di↵erence and then using a theorem on diagonalization and the Hadamard product from [24], we have:
Finally, we have that
The above arguments prove the first part of the Theorem. Now we want to analyze the lower bound of
Denote C := wT , Z := w T . By the Cauchy-Schwartz inequality, we know that
The equality is obtained when has the same value on all coordinates of W . Now check that w T ( ) (2)
The last inequality achieves the lower bound with Z = k 1+k C through the standard analysis of univariate quadratic equations. Based on the above arguments, we have that
Since the lower bounds from the Cauchy-Schwartz inequality and the univariate quadratic equation are all tight, we can conclude that the final lower bound is tight and it can be achieved with
This finishes the proof.
Proof of Theorem 3. For any subset S ✓ L, define F (S) =
Note that this is the value we want to minimize (i.e., C 1 + C 2 ). Without loss of generality, we may assume that none of the chosen eigenvalues in S is 0, since it does nothing to change F (S).
We first prove the first claim of Theorem 3.
Proof. First, we show that an optimal solution S of size k exists. Assume we have some optimal solution S of size less than k. Let b 2 L \ S. Then we have:
Now it is clear that if H(S) is 0, we can add any eigenvalue to S and F (S) would not change. If H(S) is not 0, then adding an eigenvalue to S that matches the sign of H(S), would reduce F (S) and the opposite would increase F (S). We will continue to make use of this fact. Adding an eigenvalue that matches the sign of H(S) is also always possible since H(S) is the sum of the eigenvalues we are choosing from. Thus, we can just continue to add eigenvalues to S until it is of size k.
Next, assuming S is of size k, we show that it must be the one as described in Theorem 3. We can prove this by contradiction. Without loss of generality, we may assume that S contains a positive eigenvalue b such that there exists some a 2 L \ S satisfying a > b. First, we observe that since S is optimal:
which reduces to H(S) b. Then, we have:
We know that b a < 0 and b > 0, H(S) b, so we get that the expression is negative. Thus, S [ {a} \ {b} has a lower value under F and contradicts our assumption that S is optimal.
Let S be the optimal solution from Theorem 3. Let T be the solution obtained by EV-Selection. We will show that F (S) = F (T ). We use contradiction again. Without loss of generality, we assume that S has more positive eigenvalues than T . Now consider H(S \ T ). From a similar argument used in the proof of Theorem 3 (1), due to S \ T containing only positive eigenvalues, H(S \ T ) 0. Now if H(S \ T ) = 0, we know that S \ T contains one eigenvalue, -otherwise, removing a positive eigenvalue from S would decrease F (S). Then, we know that F (S) = F (T ) and that is a contradiction. If H(S \ T ) > 0, then we have reached a contradiction in the definition of T because EV-Selection will only choose eigenvalues in T \ S when all negative eigenvalues in T \ S, Q, has been chosen. However, for all subsets P of positive eigenvalues of T , H(P [ Q) > 0, so no eigenvalue in T \ S will ever be chosen by EV-Selection.
Note that this is the value we want to minimize (i.e., C 1 + 1 |S|+1 C 2 ). Without loss of generality, we may assume that none of the chosen eigenvalues in S is 0, since it does nothing to change F (S).
We first prove the first claim in Proposition 8.
Proof. First, we show that an optimal solution S of size k exists. Assume we have some optimal solution S of size less than k. Let b 2 L \ S. Then we have:
Now it is clear that if H(S) is 0, adding any eigenvalue to S would not increase F (S). If H(S) is not 0, then adding an eigenvalue to S that matches the sign of H(S), would reduce F (S). We will continue to make use of this fact. Adding an eigenvalue that matches the sign of H(S) is also always possible since H(S) is the sum of the eigenvalues we are choosing from. Thus, we can just continue to add eigenvalues to S until it is of size k.
Next, we will show there must be an S as described in Proposition 8. We can prove this by contradiction. Without loss of generality, we may assume some optimal S contains a positive eigenvalue b such that there exists some a 2 L \ S satisfying a > b. We can directly observe since S is optimal
Since we know b a < 0, we get that (|S| + 2)b + |S|a + 2H(S) 0. First, for the case where it is equal to 0, we can just take out b and put in a. If S still does not have the greatest positive eigenvalues, we can repeat our analysis and examine another candidate positive eigenvalue to be replaced. In the case that it is strictly negative, then we have H(S) < 0 since |S|, b, a > 0. Then there must exist some negative eigenvalue c that has not been chosen and if we replace a with c in our analysis, we get that b c > 0 and (|S| + 2)b + |S|c + 2H(S) < 0 which gives F (S [ {c} \ {b}) < F (S) which contradicts the optimality of S.
Let S be the optimal solution from Proposition 8. Let T be the solution obtained by ?? 3. We will show that F (S) = F (T ). Note that now that we know |S| = k, we now define
We use contradiction again. Without loss of generality, we assume that S has more positive eigenvalues than T . Let a be the largest eigenvalue in magnitude from T \ S. Let T 0 be the set chosen by ?? 3 right before choosing a. Then let b be the largest eigenvalue in S \ T 0 . By definition of T , we know that:
Now let's compare that with:
Now clearly H(S) + b H(T 0 ). In the case that H(S) + b = H(T 0 ), we clearly have S [ {a} \ {b} = T , so we get F (S) = F (T ). In the case that H(S) + b < H(T 0 ), we also know that a b < 0, so we get
That violates the fact that S is an optimal set, so we have a contradiction.
We compare the classical multidimensional scaling and our non-Euclidean multidimensional scaling under the random matrix context, and give the proof of Theorem 5.
We start with a symmetric random matrix B 2 R n⇥n as in Proposition 4, and want to select k eigenvalues by using the classical and non-Euclidean multidimensional scaling.
We start with classical multidimensional scaling. To select the largest k eigenvalues, we need to select all eigenvalues greater than (2 r) p n such that
). (38) There is not an explicit way to write r as a function of k. For a fixed r, since we want to drop all eigenvalues smaller than (2 r) p n , we have e
For non-Euclidean multidimensional scaling, to select the k eigenvalues with largest magnitude, we need to select all eigenvalues with magnitude greater than (2 r)
For a fixed r, since we want to drop all eigenvalues with magnitude smaller than (2 r) p n , we have
We first consider the case that k = o(n), i.e. we want to reduce the dimension from n to a much smaller k. Under this assumption, r is a very small positive number. In the classical case, by applying the Taylor series to equation (38), we get r ⇡ ( 3⇡k 2n )
2 3 . Then we plug it into equation (39) to get
Similarly, in the non-Euclidean case, we apply the Taylor series to equation (40), we get r ⇡ ( 3⇡k 4n )
2 3 . Then we plug it into equation (41) to get
This finishes the proof of Theorem 5 (1). Now we turn to the case that k = cn for a constant c 2 [0, 1]. We first solve c = k n in equations (38) and (40) to get the corresponding r, then plug the r-values in equations (39) and (41) to get the corresponding e C and e N , respectively. Since equations (38) and (40) can not be solved explicitly, we can only get numerical values of e C and e N . For c 2 (0, 0.5], numerical values of e C and e N as shown in Table 4. For c 2 [0.5, 1], the classical multidimensional scaling stabilizes, with e C ⇡ (0.5 + 0.1801 • n)n 2 2 . Meanwhile, the error e N for non-Euclidean multidimensional scaling decreases to zero as c increases to 1. More precisely, if c = 1 ✏ with very small ✏ > 0, equations (40) and (41) give
Some numerical values of e N with c 2 [0. Table 5: This table shows the error term e N of di↵erent choices of c 2 [0.5, 1) with c = k/n, normalized by dividing n 2 2 .
In this appendix we fill the missing details in the experiment section. Our experiments are implemented with Intel Core i9 CPU of 32GB memory, no GPU is required and the execution time is no longer than 30 seconds. In the rest of this section, we first introduce the generation process of synthetic datasets in Appendix E.1; Next we define the evaluation metrics and perturbation approaches in Appendix E.2; More results on dissimilarities regarding the scaled additive error, negative distances and negative eigenvalues selected, are demonstrated in Appendix E.3. With respect to the dimensionality paradox, we provide illustrations similar to Figure 2 on other datasets in Appendix E.4.
The random-simplex dataset is generated with the idea of creating a simplex with one added dimension that generates all of the variance in the the distances such that the added dimension corresponds with a negative eigenvalue. This way, we should have a dimension with a large negative eigenvalue that is impactful on the stress of the embedding. Specifically, this dataset was generated with points x 1 , ...x 1000 , each point having 1000 dimensions. For each point x i , the first 100 coordinates were chosen uniformly randomly between 0 and .01. The next 899 coordinates were chosen uniformly randomly between 0 and p 0.5/899. The last coordinate is i ⇥ 0.3/1000. Then, the distance between x i and x j was found in the following way:
where x i (k) is the kth coordinate of x i . The first 999 coordinates essentially form a simplex, because selecting a large number of coordinates randomly from the same distribution causes all of the distances to converge to a normal distribution with low variance. Thus, all of the distances between points using just the first 999 coordinates are very similar. Then, the last coordinate has "large" discrepancies between all of the points, causing most of the di↵erences in the distances.
The Euclidean-ball metric was inspired by [49]. The set of distances between objects in space does not follow triangle inequality, so it becomes a non-euclidean dissimilarity measure. The particular dataset we use generates the distances by first randomly selecting 1000 points uniformly in a 10 dimensional hypercube of length 100. Each of these points is now the center of a ball. Then, for each ball, with probability 0.9, we select a radius uniformly randomly between 0 and 5. With probability 0.1, we let the radius be 0.8 times the distance between the center of this ball to the closest other ball. This closest other ball is based on the radius of the other ball if it has already been decided. The final distance matrix is just filled by the distances between the 1000 balls. In this way, we attempt to create some balls with large radius, so that the triangle inequality is heavily violated, therefore causing the dataset to deviate more from the Euclidean setting.
In addition to STRESS, our evaluation metric include scaled additive error and average geometric distortion. Scaled additive error by definition is allowing scaling of the distance matrix of our embedding before calculating stress. To do this, we first ran the embedding method (e.g., cMDS, Neuc-MDS, or other methods) on the given dissimilarity squared matrix D to get the output dissimilarity squared matrix D. Then, we would flatten both D and D to vectors, and project D onto the line through D to get E. Then we calculate kD Ek to get scaled additive error. Clearly, this is equivalent to allowing scaling of D to find the minimum stress. For average geometric distortion, we again input D to get D, but this time, we evaluate the distortion of each dissimilarity by dividing entries of p D by entries of p D. Since distortion is not well defined when dissimilarities become complex, we skip those. Then, with the remaining valid distortions, we scale the distortions so that there's an equal number of them greater than 1 and less than 1, then for those that remain less than 1, we take their reciprocal. Then, we finally take the geometric average of all of these distortions. In this way, we basically again allowed scaling of the dissimilarities to find the minimum possible average distortion.
Next, we had to decide on metrics to use for the image datasets. We took the following 3 metrics from [38]. For the first metric, the distance between two images is first finding the Euclidean distance between their coordinate representations, and then adding Gaussian noise to the distances. Because the noise is added directly onto the distances, the end result is highly likely to be non-Euclidean. For the second metric, we build a k nearest neighbor graph of the images based on their Euclidean distances. Since the end result is a graph structure, the shortest path distances are also highly likely to be non-Euclidean. For the third metric, we randomly removed entries from the coordinate representation of the images. Then, the distance between two images is decided by taking the Euclidean distance between coordinates that both images still included. Here, the Euclidean property breaks down again because triangle inequality would no longer have to hold. In Table 2, we choose the k-NN metric with k = 2. In the next Appendix section, we show results for higher k and other metrics applied to the image datasets. As for genomics datasets, we use the entropic a nities commonly used in t-SNE methods. The main idea is to apply an adaptive kernel, whose bandwidth depends on a parameter termed as perplexity, to pairwise data entries.
In addition to STRESS and average distortion already reported in Table 2, we further show the third metric: scaled additive error as defined in Appendix E.2. Note that our algorithms may produce negative distances, we also report number of negative distances and number of negative eigenvalues selected. All original results (without scaling) on 10 datasets are demonstrated in Table 6. Recall we mentioned the results in Table 2 contains some scaling on the synthetic data and images, we also provide the original values in Table 7.
On scaled additive error, our proposed methods still consistently outperform cMDS and Lower-MDS considerably. Thus, we have shown that on all metrics our methods yield more favorable results. It is interesting to observe that Neuc-MDS + produces much fewer negative distances than Neuc-MDS, though the number of negative eigenvalues selected are quite close. Further, two methods have similar performance on STRESS and other metrics. Therefore in practice, if there is a considerable concern on negative distances, Neuc-MDS + should become a better choice. Until now all results reported on images are using kNN metric with k = 2 to produce non-Euclidean dissimilarities. There could be other perturbations for the same purpose, as mentioned in Appendix E.2. We further test with k = 10 and the other two metrics: adding noise and random removal on the STRESS. We add noise sampled from Gaussian distribution with variance as the maximum instance-wise di↵erence scaled by 500; For random removal, we skip 50 images for each. The results are shown in Table 8.
The purpose of this section is to show the consistent performance of Neuc-MDS and Neuc-MDS + on all datasets, with respect to mitigating the dimensionality paradox issue. Therefore, in addition to Figure 2, we provide the same plots for other datasets, and the parameters all follow the main setup. Figure 3 gives a clear illustration that Neuc-MDS and Neuc-MDS + always have lower STRESS than cMDS and Lower-MDS as dimension grows.
https://github.com/KLu9812/MDSPlus
For the sake of clarity we scaled some results by a factor of 10 3 . We did so on all methods to ensure we did not engender bias on analysis. We provide the original results in Appendix E.3.