A geometrical connection between sparse and low-rank matrices and its application to manifold learning
Authors: Lawrence K. Saul
Venue: Accepted by TMLR ()

Abstract: We consider when a sparse nonnegative matrix $\mathbf{S}$ can be recovered, via an elementwise nonlinearity, from a real-valued matrix~$\mathbf{L}$ of significantly lower rank. Of particular interest is the setting where the positive elements of $\mathbf{S}$ encode the similarities of nearby points on a low dimensional manifold. The recovery can then be posed as a problem in manifold learning---in this case, how to learn a norm-preserving and neighborhood-preserving mapping of high dimensional inputs into a lower dimensional space. We describe an algorithm for this problem based on a generalized low-rank decomposition of sparse matrices. This decomposition has the interesting property that it can be encoded by a neural network with one layer of rectified linear units; since the algorithm discovers this encoding, it can also be viewed as a layerwise primitive for deep learning. The algorithm regards the inputs $\mathbf{x}_i$ and $\mathbf{x}_j$ as similar whenever the cosine of the angle between them exceeds some threshold $\tau\in(0,1)$. Given this threshold, the algorithm attempts to discover a mapping $\mathbf{x}_i\mapsto\mathbf{y}_i$ by matching the elements of two sparse matrices; in particular, it seeks a mapping for which $\mathbf{S}=\max(0,\mathbf{L})$, where $S_{ij} = \max(0,\mathbf{x}_i\!\cdot\!\mathbf{x}_j\! -\! \tau\|\mathbf{x}_i\|\|\mathbf{x}_j\|)$ and $L_{ij} = \mathbf{y}_i\!\cdot\!\mathbf{y}_j\! -\! \tau\|\mathbf{y}_i\|\|\mathbf{y}_j\|$. We apply the algorithm to data sets where vector magnitudes and small cosine distances have interpretable meanings (e.g., the brightness of an image, the similarity to other words). On these data sets, the algorithm is able to discover much lower dimensional representations that preserve these meanings.

Published in Transactions on Machine Learning Research (12/2022)
A geometrical connection between sparse and low-rank matrices and its application to manifold learning
Lawrence K. Saul
bonl@ffationninsttutt.oy
Center for Computational Mathematics Flatiron Instiutt
162 Fifth Avenue New York, NY 10010
Reviewed on OpenReview: Attps: // openr/veewwwww.. net/ formm? id- p8gncJ.Mit
//tps://owe
Abstract
We consider when a sparse nonnegative matrix S can be recovered, via an elementwise non- linearity, from real-valued matrix of significantly lower rank. Of particular interest the setting where the positive elements of S encode the similarities of nearby points on a low dimensional manifold. The recovery can then be posed problem in manifold learning in this case, how to learn norm-preserving and neighborhood- preserving mapping of high dimensional inputs into a lower dimensional space. We describe an algorithm for this prob- lem based on a generalized low-rank decomposition of sparse matrices. This decomposition has the interesting property that can be encoded by a neural network with one layer of rectified linear units; since the algorithm discovers this encoding, it can be viewed as a layerwise primitive for deep lcarning. The algorithm regards the inputs æ; and x, as similar whenever the cosine of the angle between them exceeds some threshold E (0, 1). Given this threshold, the algorithm attempts to discover mapping æ, + y, by matching the elements of two sparse matrices; in particular, it seeks mapping for which S = max(0, L), where | []) and Lij = yı yı -T|lBilll 3,ll- We apply the algorithm to
Sig = max(0, x, z - || data sets where vector magnitudes and small cosine distances have interpretable meanings (e.g., the brightness an image, the similarity to other words). On these data sets, the algorithm discovers much lower dimensional representations that preserve these meanings.
1 Introduction
Many algorithms for high dimensional data analysis succeed by capitalizing on particular types of structure in large matricos. One important type of structure sparsity; for example, when a matrix is sparse, with large number of zero elements, it can be stored a compressed format [Dulf et al., 2017). Another type of structure is linear dependence; when matrix is rank deficient, with one or more singular values equal zero, it can be expressed as the product of two smaller matrices (Horn & Johnson 2012 It is well known that neither one of these structures implies the other. For example, very sparse matrix can be full rank (e.g., the identity matrix), and conversely, very low rank matrix can be dense (e.g., the matrix of all ones).
It would be incorrect, however, to say that these types of structure are unrelated. One can find more subtle connections between rank and sparsity by looking beyond the canonical decompositions algebra (Udell et al., 2016: Wright & Ma (021). One such connection was explored recently for sparse monnegative matrices (Saul 2022). In particular, let. S denote such a matrix. Then one can find a real-valued matrix L of equal or lower rank such that
linear
S x max L),
(1) applied independently to all elements of L. Note that a
where the piccewise linear operation in eq. trivial result of this form is obtained by setting L S (where the approximation is exact). A more intriguing possibility emerges from the following observation: the sparser the matriz S, the larger the space that can be

---

Published in Transactions on Machine Learning Research (12/2022)
max(0, I.)
I WH
10
0
0
15
10
-11
4
-11
-44
15
max(O. (0, Wh;)
0
13
0
16
24
0
-11
13
-1
16
24
-15
S = 81
34
0
-1
34
-19
-15
0
29
16
-19
29
-10
24
24
16
64
-24
24
24
16
64
-40
H = [h h2 hm
13
0
0
25
15
-15
-13
-10
-40
25
rank(S) =
rank(L)=2
Figure 1: Left: A sparse nonnegative matrix can be expressed via eq. in terms of matrix L of equal lower rank. The matrices share the same positive elements (shown in blue), but the zeros of S may correspond to either zeros or negative values in L (shown in red). Large geps in rank may become possible when S is very sparse. Right: The columns of can be encoded by neural network with top layer of rectified linear units and a bottom layer of d hidden units, where d= rank(L).
explored to discover matriz L of significantly lower rank. The intuition behind this observation is illustrated in the left panels of Fig. [ An even more striking example is the following. Let denote the n x n identity matrix; then for all n, there exists a matrix L such that I = max(0, L) and rank(L) S 3 (Saul 2022).
The computation in eq. (1) is interesting for another reason namely, that it can be implemented by a neural network with one layer rectified linear units (ReLUs) (Mazumdar & Rawat 2019a). To see this, suppose that the low-rank matrix in eq. (1) factored as WH. Then we can view the columns H as hidden units that encode the columns of S, and we can view the rows of W weight vectors that reconstruct the rows of S from these encodings. This idea is illustrated in the right panel of Fig. [
Eq. suggests many possibilities for high dimensional data analysis, especially for data that is naturally represented as a sparse nonnegative matrix (e-g,, n-gram counts (Michel et al. 2011), affinities in a social network Newman
2003 But eq. by itself has several limitations, of which we briefly mention three:
1. It is not obvious how to interpret the nonlinear encodings of data that obtained in this way. In particular, they suggest neither a subspace model of the data, as might be obtained by singular value decomposition (Eckart & Young 1936 Deerwester et al. 1990 Turk & Pentland 1991). nor a
whole- parlee-ppssssee decomposition, as might be obtained by nonnegative matrix factorization
Lee & Seung
1990 Gillas 2021) or sparse dictionary method (Olshausen Field
1996: Chen al. 1998 Mairal
2010
Rubinstein al. 2010). For many data sets, a geometrical picture is lacking.
2. It is not obvious how to apply the transformation in eq. 1 to data that not sparse, even though many studies have shown such data can be modeled by neural networks with rectified linear units (Hinton Ghahramani 1997 Radford et al. 2016
3. Based the single-layer network in Fig. I one might be tempted to imagine eq. [ greedy, layerwise strategy for unsupervised lcarning in deep networks (Hinton & Salakhutdinov 2006 Hinton et al. 2006 Mazundar & Rawal 20195). But Fig. T also reveals an immediate complication: unlike the original matrix S, the matrix H of hidden encodings is neither sparse nor nonnegative. Thus it not obvious how to compose this method in a way that generalizes to deep architectures.
In this paper we describe model, also based on eq. [ that overcomes the first two of these limitations and suggests a way around the third. The model is designed to learn faithful low-dimensional representations of nigh-dimensional inputs in an unsupervised and highly interpretable manner. The inputs are not required be sparsc. This the so-called problem of nonlinear dimensionality roduction, or manifold learning, and
it has been widely studied (Tenenbaum et al. 20000] Roweis & Saul 2000 Belkin & Niyogi 2003 Coifman
et al. 2005 Weinberger & Saul| 2006 van der Maaten & Hinton
2008 Mclnnes al. 2018 Agrawal et al.
2

---

Published in Transactions on Machine Learning Research (12/2022)
2021 Nevertheless, the model in this paper offers novel solution, powered by eigenvector methods, based on the neural network and nonlinear matrix decomposition in Fig. [ When applied to images, it is also designed to learn -euuiiaraaiv representations (Nair & Hinton 2010 Hinton et al. 2011), and thus it differs from previous approaches that do not proceed with this goal in mind.
Most of this paper devoted to describing this model and evaluating the simplest algorithm that suggests to implement. But the main contribution of this paper is a conceptual one: it is to attach a geometrical picture to eq. ① To begin, suppose that data set of high-dimensional sensory inputs has been sampled from low-dimensional perceptual manifold (Seung & Lee 2000). The problem of manifold learning can be viewed, at a high level, as one of translating between the two different types of matrices that arise in this setting. On one hand, there are sparse matrices, whose nonzero elements reveal at a glance which high dimensional inputs are to be regarded as similar or nearby on the manifold (Sengupta et al., 2018 Chen ct al. 2018); on the other hand, there are low-rank matrices whose nonzero eigenvalues reveal at a glance the dimensionality of some Euclidean space which the manifold can be embodded (Tenenbaum et al., 2000; Weinberger & Saul 2006). It. is in this sense that certain sparse and low-rank matrices can be viewed "two sides" of the same manifold, each revealing different aspects of its structure. The main contribution of this paper is the following: when properly formulated, these different matrices- sparse, one low- rank-are are mathematically connected by eg. 4 This connection is both surprisingly simple, involving only a piecewise linearity, and also deeply familiar, involving the same piecewise linearity that is most widely used in applications of deep learning. One might even say that the goal of this paper is to reveal a geometrical connection between sparse and low-rank matrices that was hiding in plain sight.
The organization of this paper is as follows. In section we develop model for nonlinear dimensionality reduction in which input similarities are encoded by sparse matrix S of thresholded inner products. Given this sparse matrix S, the model's embeddings are obtained by solving a sequence of eigenvalue problems for the corresponding low-rank matrix L in eq. [ In section 3 we describe one way to extend this model to larger data sets where it is too expensive to construct and/or or solve these eigenvaluc problems directly. This exlension requires only the solution of a sparse least-squares problem. In both of these sections, we also apply the model to images of handwritten digits (LeCun et al.) 1998) and pre-trained word vectors (Mikolov ct al. 2018 Finally, in section 4 conclude by discussing multiple avenues for future work.
2 Thresholded similarity matching
We begin by stating the problem that will ultimately lead us back to eq. [ Let X' = (z), ......,,,
denote a data set of digeniiiensional inputs in RD, Our goal to discover a corresponding set of outputs
y {31, ./...... 8n in Rdd that provide a faithful but much lower-dimensional reprosentation (or embedding) of the data (with d D). Section 2.1 defines what we mean by faithful in this context, section 2.2 describes algorithm to compute the embedding, section 2.3 presents the results of this algorithm two data sets of interest, and section 2.4 relates our approach to previous work on this problem.
2.1 Model
What does it mean for an embedding to be a faithful representation? We start from a familiar place, measuring similarities between vectors via inner products. As shorthand, let cos(x, | denote
ar.e"
the cosine of the angle between the vectors æ and w'. Then we define an embedding to be a 7-faithful representation of the data (with 0<+<1) if it meets the following three conditions:
(i) Magnitudes are preserved:
|| = ||>||| for all i. cos yj) = cos(2,, mz) whenever cos(a;; > T. cos( y,) S r whenever cos(2,, mz) ST.
Small angles are preserved:
Large angles remain large:
Likewise we shall say (informally) that an embodding an approximately r-fatthful representation if these conditions are approximately met. Note that equalities appear in the first two conditions, but not the third, In particular, large angles do not need to be preserved; they only need to remain larger than cos-
3

---

Published in Transactions on Machine Learning Research (12/2022)
We emphasize the following point: it is by design that these conditions refer to magnitudes and angles opposed to pairwise distances (or proximities), as is more commonly done in studies of nonlinear dimension- ality reduction (TTenenbamm et al. 2000 Belkin & Niyogi, 2003) Weinberger & Saul 2006 van der Maaten & Hinton 2008 McInnes et al. 2018 Agrawal et al.) 2021). The first condition is particularly apt for data sets which the magnitude of inputs convey meaninglul information their own right. In this paper, for example, we apply the model to data sets of images and word vectors. Note that the magnitude of pixel vector measures the overall brightness of an image, while the magnitude of word vector reflects a word's specificity (Schakel & Wilson 2015) (i.e., whether it can be used in one or many different contexts). Thus the origins in these spaces have definite meaning (e.g., a blank image, a nondescript. word), and the first condition requires that the origin in the input space is mapped uniquely to the origin in the output space, thus preserving this meaning. In fact, the second condition requires further that parallel inputs are mapped to parallel outputs. This is necessary to learn intensity-equiivrian mappings (Nair & Hinton 2010, Hinton ct al. 2011), as is desirable in sensory domains where the perceived strength of a signal is coded by its mag- nitude (e.g.: the brightness of images, the volume of sounds). It should also be noted that these conditions are not appropriate for all data sets. For example, they would thwart the "unfoldings" of many popular but artificial data sets (e.g., Swiss rolls) where the magnitudes of inputs do not have this interpretation.
The name of our model is derived from the following observation. Note how all three of the above conditions are compactly captured by the single equation:
max x. - x.. (0yl)) = max y: y, - 1/y: bs)l))
(2)
Now suppose in this equation that the threshold E 10) exceeds the cosine of the angle between most pairs of imputs in X. Then the left side will define highly sparse matrix, and r-faithful embedding will only be obtained if the same thresholded matrix is produced on the right side by the outputs in V Next consider the dense matrix yi-Ui - y ||ااا on the right side of eq. (2). If the outputs {y: live in Red then the rank of this matrix can be at most d+1. Hence if there does exist a f-ittthful embedding into R then there also exists a matrix of rank at most d+1 from which the sparse matrix on the left side can be recovered. These are the essential idcas behind our approach, which we call Thresholded Similarity Matching (TSM).
The sparsity of the matrices in eq- 2 is governed by the threshold t, which must be chosen manually. Given data set X, reasonable value for can be determined by examining the histograms of angles between nearest neighbors (using cosine distance) versus other pairs of inputs; the goal is to idcntify a threshold above which inputs can be meaningfully regarded as similar. Some of these histograms are shown in Fig. 2 for two data sets-one of images of MNIST handwritten digits (LeCun et al. 1998), and one of pre-trained word vectors from a FastText embedding (Mikolov et al.) 2018). (More details on these data sets are given in section 2.3) ) The figure highlights particular choices of T for these data sets; for these choices, the matrix the left side of eq- (2) is 99% sparse but has nonzero entries for nearly all pairs of one-nearest neighbors (INN). Put another way, r-faithful embeddings with these values of r will preserve the angles between the top one-percent of closely aligned inpuls. Note in these embeddings that angles are preserved not only between 1NNs, but also between all pairs of imputs whose cosines exceed the threshold T.
2.2 Algorithm
Next we describe how to compute an approximately r-faithful embedding of some specified dimensionality d. Our algorithm proceeds in three steps.
Step 1. Compute the sparse matrix S. To begin, we identify the left side of eq. (2) with the sparse matrix S in eq. 1
Siy = max(0, x; 1 - ||==>||||=,,|||| |||=,||)-
(3)
The algorithm starts simply by computing this matrix. There only one possible complication in this step; it arises when row or column of contains only a single nonzero element along the diagonal. This situation occurs when input a; (typically outlier) is not within an angle of cos cos) from any other inputs
4

---

Published in Transactions on Machine Learning Research (12/2022)
DIGITS
WORDS
INN other
IMN
other
15
30
45 angle (*)
60
75 90 30
75
90
angle (*)
Histograms of angles between inputs that are one-ncarest neighbors (INNs) blue versus those that are not in orange. The histograms are separately normalized as there are many fewer pairs of inputs in the first calegory. The histograms were used to set the threshold, т, that identifies pairs of similar inputs. We set 0.75 for images of handwritten digits (left) and 0.47 for word vectors (right); the angles corresponding to these values are indicated by the red vertical lines. For these choices, the matrix on the left side of eq- (2) was 99% sparse. Note that the overlap between histograms is small but nonzero, and for both data sets there remain a large number of kNN (with E> 1 but small) that are regarded as similar.
Figure/ 2:
besides itself. In this case, the matrix S may not provide sufficient information to co-locate the output ye in the low-dimensional embedding (indicaling only that it should be far from every other ouput). In practice, when this situation occurs, we augment the data set X' to include virtual inputs that interpolate between æ; and its nearest neighbor and that differ by angles less than cos- Note that if this situation occurs for more than a small fraction of inputs, then it suggests to lower the threshold T in eqs. (23
Step& 2. Compute the low-rank matrix L.
Next the algorithm attempts to find a low-rank matrix L from which the sparse matrix S in eq. (3) can be recovered via eq- 1 In our experience, the algorithm works best when the rank of L is equated to the target dimensionality (d) of the desired embedding, but higher choices are also possible.
For this step we adapt an alternating minimization procedure that was recently proposed for this prob- lem [Saul 2022 This procedure alternates between two updates-on one for the matrix L, and another for auxiliary matrix Z of the same size. The updates attempt to solve the constrained optimization:
min L - z | such that { d rank (L), S = max (0, Z).
Note that the objective function in eq. is bounded below by zero, and it only attains this minimum value when S = max(0, L) and L is a matrix of rank d. Though not explicitly indicated in eq. (4) both L and Z are understood to be symmetric matrices (like S).
The alternating updates for this optimization take an especially simple form. When the matrix L is fixed, the objective function minimizod by sctting
{
if Sig > 0, min Laj) otherwise.
Siy
Ziz
(5)
Likewise, when the matrix Z is fixed, the objective function minimized by setting L to the best low-rank approximation of Z measured by the Frobenius norm in eq. 4 Thus the update for L takes the form:
L arg min |M
2 such that rank(M) d.
(6)
M
Since is symmetric (though generally not positive definite), the solution eq. (6) obtained by setting
QDQ T where D is the d'd diagonal matrix whose diagonal elements specify the largest eigenvalues (in magnitude) of Z and Q is the corresponding x d matrix of normalized eigenvectors.
1 All baspd approaches noulinear dimensionality reduction are plagued some extent by outlicrs. Other alternatives (not considered here) are analyze only the largest connected component inputs (Tenenhaum et 2000) construct. neighhorhood graph via minimum spanning trees Zemel & Carreirss-Prriinánn 2004), thus ensuring that connected.
5

---

Published in Transactions on Machine Learning Research (12/2022)
The optimization in eq. 4 is not convex (due to the constraints on L and Z), and hence this procedure is not guaranteed to converge to global minimum. By construction, however, each of the updates in eqs. 5 reduces the value of the objective function (except at stationary points), and thus this property provides a
6
guarantec of monotonic improvement. In practice, the updates can be alternated as long as the objective
lunction is decreasing at each iteration by some small amount. The algorithm works on the same general principle as other alternating minimizations, such as for Appectation-Maximizzaiii
(Dempster et al. 1977 Meng & van Dyk 1997 McLachlan & Krishnan 2008) and nonnegative matrix factorization (Lee Seung, 1999 Gillas, 2021] and thus similar empirical strategies can be used to test for its convergence. We refer the reader to appendix A for further details of implementation; these details include the choice of initialization,
an additional linear constraint on the elements of Z, and the use of a momentum parameter, all of which lead to improved solutions and/or faster convergence.
Step 3. Compute the embedding V
The last step of the algorithm is to recover the d-dimensional outputs (2),.... Un.. To do so, we identify the matrix argument on the right side of eq. 2 with the low-rank matrix L computed in the previous step:
Laj U; =llslllles
The problem now is to recover the outputs y; E Rd from the known elements of L. A similar problem arises in multidimensional scaling (Cox & Cox 2000). where one recovers n vectors (centered on the origin) from n x n matrix of pairwise squared distances (o.g., |/y =,,l)))) The problem here is analogous, except that we are recovering the outputs from the norms and inner products on the right side of eq. 7 moreover, the outputs in our problem are not required to be centered.
We first derive an exact solution to this problem under the assumption that one docs, in fact, cxist. Assuming that eq. 77) holds for some outputs y E Rd, then it follows from the diagonal elements that Lu = (1-T)lly:1
or equivalently, || = VT (1-r) 11/2. Next we substitute this result into the right side of eq. 7 This yields an expression for the Gram matrix G whose elements Gij = yı yj store the inner products of the outputs. In particular, have:
Gig = Lig + 1-7 Linyi
(8)
The outputs are then determined (up to an arbitrary rotation) by the elements of this Gram matrix, thus completing the solution.
The above solution was based on assumption nannny, namely, that there exist outputs whose inner products are perfectly consistent with the matrix L computed in the previous step of the algorithm. When this is not the case, we compute an approximate solution that is optimal in a least-squares sense (Cox & Cox 2000 Tenenbaum et al. 2000). Let {Aald= denote the d largest eigenvalues the matrix G, and let {ea]a-) denote the corresponding eigenvectors. Then an approximate solution to eq.
given by
Yoi V mmax (0, An) Cin- ein.
When the matrix G from eq. is positive semidefinite, this approximation computes the outputs in g8d (up to global rotation) that minimize the error Eng (G:--y::)) y, When this is not the case, the approximation proceeds by ignoring the negative eigenvalues of G. (Large or multiple negative eigenvalues suggest that one should choose a higher value for the target dimensionality d.)
Practical considerations. The algorithm proceeds from choices of the threshold in step 1 and the target dimensionality d in step 2. It should be clear that these choices are interrelated. The smaller the value of r, the larger the number of inner products that the algorithm attempts to preserve in steps and 2. Likewise, the smaller the value of d, the less flexibility there is to discover a semilaritty-reserviin embedding in steps 2 and 3. The main question the following: is there a meaningful regime for these choices in which (a) the optimization in eq. 4 not eviscerated by local minima, and (b) the low-rank matrix L satisfying S = max L) yields sensible Gram matrix via eq. empirical evaluations, to which turn next.
( ? This question can be answered to some degree by
6

---

Published in Transactions on Machine Learning Research (12/2022)
DIGITS (14K)
WORDS (25K)
50
SVD
-0-TSM
-0-TSM
30
MEAN ANGULAR 15
DEVIATION (9)
10
16
64
128
256
4
16
64
128
256
dimensionality
dimensionality
Figure 3: Mean angular deviations in eq. (left) and 25K word vectors (right), as function of their dimensionality.
11 for the embeddings from TSM and SVD of 14K digit images
2.3
Evaluation and results
Data sets. We evaluated the algorithm of the previous section on two data sets where it is common to measure similarity using cosine distances:
Handwrillen digit images
The first data set was the MNIST data set of grayscale handwritten digits (LeCun et al., 1998). In this data set, each image encoded by a 784-dimensional vector of pixel values. The data contains total of 70K images, which we divided into five non-overlapping folds of 14K images. We ran the algorithm separately each of these folds.
2. Pre-trained word vectors
The second data set consisted of subset of pre-trained 300-dimensional word vectors from FastText embedding of 1M English tokens (Mikolov al. 2018). For the experiments in this section, we the algorithm on a subset of approximately 25K word vectors that were selected to obtain sub-vocabulary of mostly contentful ful words. This sub-vocabulary included the leading 18K words of the Text embedding that are valid tournament plays the game Serabble"] as well as all 905 unique words in the Semantic- Syntactic Word Relationship test for analogy completion (Mikolov et al., 2013a b). The latter were included so that we could use analogy completion to assess different dow-dimensional embeddings of word-vectors; to do this, it was necessary to supplement many proper nouns, such as state capitals, that appear in these analogies but are not valid plays in Scrabble. Finally, we also included the 15-ncarest neighbors (in cosine distance) for each of these 905 word-vectors so that they did not appear as isolated instances in our data set but were linked to other words with similar mcanings. This brought the total sub-vocabulary to 25124 words.
sub-vocabulary with (for instance) the
The data sets were restricted size for these initial evaluations because the updates in eqs. (-6) work with dense n matrices (namely, Z and L), and these matrices must fit comfortably into memory. We evaluate more memory-efficient way to handle larger data sets in section
Baselines. To obtain baseline embeddings, we performed a truncated singular value decomposition (SVD) of the above data sets; the inputs were then projected into the lowor-dimensional subspace spanned by the
top left singular vectors. More concretely this was done follows. Let X [ æ2 ... æ.,] denote the matrix of high dimensional inputs, and let U denote the matrix whose columns store the top d eigenvectors of XX We computed the lincar embedding + via the orthogonal projections
v;
U T.:.
(10)
2 An official tournament. word list is available htttps://wwu. http:////wwaaaaiaaaatdaerrrrreed

---

Published in Transactions on Machine Learning Research (12/2022)
DIGITS (14K)
WORDS (25K)
TSM -VD
-TTSMM
-VVVD
0.8
0.8
JACCARD INDEX
0.6
0.6
0.4
0.4
0.2
0.2
0
16
64
128
256
4
16
64
128
256
dimensionality
dimensionality
Figure 4: Jaccard indicos in eq. 25K word vectors (right), as function of their dimensionality.
12 for the embeddings from TSM and SVD of 14K digit images (left) and
Since these projections are orthogonal, the outputs in these embeddings always have smaller norms than their corresponding inputs: ||v/|| & || As discussed next, however, we chose metrics for evaluation that did not penalize the baseline embeddings for this property. These linear embeddings are optimal in a least-squares sense: if "U", then also the projection matrix of rank that minimizes the error X
Metrics. We computed two fairly intuitive metrics to evaluate the quality of embeddings 2; + As useful shorthand, let Ω2 {(i, Sy > 0} denote the set of indices where the matrix S in eq. (3) has positive (i.e., nonzero) elements, and let [9tz] denote the total number of these positive elements. The first metric we computed was the mean angular deviation
COS COS
-cos
cos
y/ Vj | y
(11)
P
|1 [|z,
SjED
Note that the sum in eq. [1 is restricted to those pairs of inputs whose angles the algorithm is trying to preserve. As another useful shorthand, let Py = ((i,j) ye yj - =|ly,|lll,,|| > 0} denote the set of indices
where the matrix on the right side of eq. (2) has positive elements. The second metric we computed was the
Jaccard index (Jaccard 1901 Levandowsky & Winter 1971):
f2 0 S2y
J(Zz, 2ny =
(12)
12 U.y
The Jaccard index lies between 0 and 1; it is equal to when exactly the same pairs of inputs and outputs have angles less than cos cr), and it is equal to 0 when none of the same pairs share this property.
Comparisons. Figures 3 and compare these metrics the embedding obtained from TSM and SVD. The embeddings were computed over wide range of target dimensionalities for the data sets of 14K digit image:f and 25K word vectors. The results are similar for both data sets: over a wide range of dimen- sionalities, the mean angular deviation is significantly lower for the TSM embeddings, while the Jaccard index is significantly higher. The results confirm the intuition suggested by Fig. the nonlinear matrix decomposition in eq. (3) is capable of deriving faithful low-dimensional embeddings by capitalizing on the sparsity of the matrix S.
Analogy completion. Analogies have been used test the meaningfulnes of adding and subtracting word-vectors in different embeddings (Mikolov al. 2013a b). Consider, for example, the following analogy:
3 AAppendix B presents an additional but less straightforward comparison to results from the Isomap algorithm. For the digit images, the results in the figures were obtained by averaging over the five different folds of the full MNIST data sct; the variation across folds, shown in the figure, was generally smaller than the marker sizes. safoguard against stereotypes and implicit biases when word-vectors used for real-world applica-
Sit is important tions (Bolukbasi et al. 2016 Gonen & Goldherg 2019 This especially when such blases might inadvertently exacerbated by dimensionality reduction.
8

---

Published in Transactions on Machine Learning Research (12/2022)
Table 1: Examples of analogies in different semantic and syntactic categories. These are also examples where the lower dimensional word-vectors from TSM (with target dimensionality d=32) predict the correct answer (shown in blue), while those from SVD predict incorrect answer (shown in red).
Type
Category
Prompt
Query
TSM
SVD
common capitals
(Athens, Greece) (Abuja, Nigeria)
Baghdad
Iraq
Afghanistan
other capitals
Amman
Jordan
Morocco
Semantic
currencies
(Algeria, dinar)
Cambodia
riel
rupee
city-state
(Chicago, Illinois)
Houston
Texas
Mississippi
family
(boy, girl)
father
mother
daughter
adjective-adverb
(amazing, amazingly)
calm
calmly
gently
opposites
(aware, unaware) (bad, worse)
honest
dishonest
cynical
comparisons
cold
colder
damp
superlatives
(bad, worst. (code, coding)
cool
coolest
awesome
Syntactic
participles
jump Austria
jumping
leap
nationalities past tense
(Albania, Albanian) (dancing, danced)
Austrian
Hungarian
describing
described
represented
plural nouns plural verbs
(banana, bananas) (decrease, decreases)
goat
goats estimates
cows
estimate
yields
"boy" to "girl" as "father" is to "mother." Given an embedding of word-vectors, one can attempt to complete the analogy automatically by computing the translated vector
tanslogy = vec vec("father") + vec(""grrl") vea("boy"))
(13)
then finding the word-vector (other than "father" or "girl") nearest to tanalogy in cosine distance. We compared these sorts of completions (limited to our sub-voaabulary of 25K words) from the original Fest embedding and the lover-dimensional embeddings obtained from TSM and SVD, The completions were evaluated on the Semantic-Syntactic Word Relationship data set (Mikolov et al., 2013a bj, which consists of 19558 analogies from 14 different calegories; of these 14 categories, are based on semantic relationships, and 9 are based on syntactic relationships. Table [ shows example from each category, and Fig. 5 plots the overall accuracy rates for analogy completion, broken down by type. For semantic relationships, we see that the lower-dimensiona embeddings from TSM approach the accuracy of the Fast'Text embedding much more quickly than those from SVD. For syntactic relationships, the embeddings from TSM are also better at small dimensionalities, though the difference is less pronounced. qualitative sense of these differences can be seen from the examples that appear in Table [ These examples were selected from analogies that only the embedding from TSM completed correctly in d2 dimensions.
2.4 Related work
The algorithm for TSM shares a familiar three-step recipe with many other eigenvector- based approaches to nonlinear dimensionality reduction. In a nutshell, the first of these steps computes a sparse matrix S from pairs of nearby inputs, the second constructs a related matrix of the same size (e.g., of graph dis- tances (Tenenbaum al. 2000), of linear reconstruction weights (Roweis & Saul 2000), of graph Laplacian elements (Belkin & Niyogi 2003), of "unfolded" inner products /Weinberger & Saul 2006 ), and the final step obtains an embedding from the top or bottom eigenvectors of this second matrix. Within this framework, however, there remain several important differences. First, previous approaches have not attempted to pre- serve the mcaningfulness of norms and angles, and most of them (for instance) do not map the origin in the input space to the origin in the output space. Second, the algorithm for TSM enforces a low-rank constraint in its alternating minimization, whereas other algorithms in this framework typically obtain a nested set of emboddings that are not optimized for a particular target dimensionality. Third, and most importantly, TSM is based on the connection in eq. ( between the matrices and L. and thus it provides a bridge between these eigenvector- baaeedd approaches and other neural-network models of unsupervised learning.

---

Published in Transactions on Machine Learning Research (12/2022)
SEMANTIC
SYNTACTIC
100
100
ANALOGY
COMPLETION
FastText
accuracy (%)
75
TSM
accuracy (%)
75
SVD
50
50
FastText -0- TSM
25
25
WORDS (25K)
SVD
0
16 32
64
128
16
32
64
128
dimensionality
dimensionality
Analogy completion rates for the embeddings from TSM and SVD of 25K word vectors, as function of their dimensionality, for semantic (left) and syntactic (right) relationships.
Figure 5:
Another line of work in dimensionality reduction has focused on embeddings for visualization. Well-known algorithms for this purpose include t-SNE (van der Maaten & Hinton 2008), UMAP (McInnes et al.) 2018). and minimum distortion embedding (Agrawal et al. (021)). These algorithms search for embeddings that balance attractive and repulsive forces between their low-dimensional outputs. In general, there are many possible choices for these forces, and they tend to be validated by the visualizations they produce. As such, they have been less explored for higher-dimensional emboddings. Unlike TSM, they rely on gradicnt- basod optimizations. With clever approximations (van der Maaten 2014 Linderman et al., 2019) and sampling- based methods (MeInnes et al., 2018 Agrawal et al. 2021 they have been scaled to very large data sets.
TSM is based on nonlinear decomposition originally proposed for sparse nonnegative matrices (Saul 2022 It builds on previous work by developing a geometrical picture for this decomposition. With this picture, eq. (I) can also be applied to data that is not intrinsically sparse; this is done by using the sparse matrix S to store thresholded similarities and the low-rank matrix L to recover a low-dimensional embedding. TSM also fits within a long line of work on generalized low-rank models that scarch for low dimensional structure
in nonnegative binary data (Lee & Sompolinsky 1999 Collins et al. 2002 Hoff 2005
Udell et al. 2016).
motivated in large part by a recent model of nonnegative similarity matching in biologically-inspieee neural networks (Sengupta et al. 2018 The authors of this work observed that localized receptive fields emerge when neural networks are trained to optimize a similarity-matchhn objective function. They also showed both numerically and analytically that the receptive fields these networks learn to tile the low-dimensional manifolds that are populated by their sensory inputs. The present work began with the realization that the responses of these receptive fields could be collected form a sparse matrix S, such as the that appears in eq- [1] The threshold in TSM also plays role analogous to the similarity threshold in their work.
Finally, note that approach
3 Locally linear extension
In this section we describe one way to extend the embeddings from TSM to larger data sets. The approach is based on familiar ideas from graph- learning (Weinberger et al.) 2005 Vladymyrov & Carreira- Parreira-Perppinn 2013 With this approach, we were able obtain results on larger data sets with only slightly more computation than the experiments in the previous section.
sumiriuper learning (Zhu et al. 2003) and large-scale manifold
3.1 Algorithm
We begin by reframing the problem. We suppose that the method of the previous section has been used to obtain an embodding of inputs, and now we wish to extend the embodding to an additional n-m inputs. We refer to the original m outputs as landmarks, and we propose two-step algorithm to compute locally linear exlension (LLX) of the smaller embedding.
10

---

Published in Transactions on Machine Learning Research (12/2022)
DIGITS (70K)
WORDS (108K)
50
TSMILLX
TSMILLX
30
MEAN ANGULAR 15
DEVIATION (9)
10
0
16
64
128
256
4
16
64 128
256
dimensionality
dimensionality
Figure 6: Mean angular deviations in eq. 11 for the embeddings from TSM/LLX and SVD of 70K digit
images (left) and 108K word vectors (right), as function of their dimensionality.
Step 1. Compute a sparse matrix W of reconstruction weights. This step similar to the first step of the algorithm known as locally linear embedding (LLE) (Roweis & Saul 2000). For each input æy we compute of & weights that reconstruct the input (approximately) from its k-nearest neighbors. Let K;; denote the set of k-nearest neighbors for input ws. Then these weights can be computed by minimizing the regularized sum of locally linear reconstruction errors,
2 + E
(x(W) = x;
Σ JEK,
Wiyays
Σ JEK,
||,,|1 W? E1
(14)
where e>> is a small value to avoid ill- ill-posed problems. There are two minor but important differences with the original algorithm for LLE: the first is that nearest neighbors are selected according cosine distance, and the second is that the rows of W are not constrained to sum one. Both differences arise from the desire to preserve the meaningfulness of norms and small cosine distances.
Step 2. Extend the embedding. The embedding is extended by assuming that nearby outputs should be linearly related in the same way as nearby inputs (Weinberger et al. 2005 Vladymyrov & Carreira- Carreira-Perpiián 2013 We therefore procoed by minimizing the analogous roconstruction error
Ew(O) = =||31 EWiyy |
(15)
over the localions of the n- unknown oulputs. To emphasize, the m landmarks in eq- 15 are fixed to their known values, while the weights Wij are determined previously by the values that minimize eq. (14] The minimization of eq. (15) reduces to d independent least-squares problems for the coordinates of the n-m unknown outputs. Most importantly, because the matrix W is sparse, these problems can be solved very efficiently using preconditioned conjugate gradients methods (Barrett et al. 1994). Further details of this approach can be found in the appendix.
3.2 Results
We experimented with this algorithm, which we denote by TSM/LLX, on the full data set of MNIST handwritten digits and a larger data set of word vectors. For both data sets, we used the algorithm to extend an embedding of m outputs (computed by TSM) to embedding of n 5m outputs (computed by LLX). For the digits, we used m 14000 and n = 70000 images, and for the text embedding, we used m m 21534 and n= 107670 word vectors. The vocabulary in this larger data set of word-vectors was expanded from the previous one to include all valid Scrabble plays that were among the leading one million tokens of the Fest embedding. For these experiments, we set 4=4d in eq. (14h and 10 in eq. (15)
Figures G and 7 show the results of these experiments, comparing the mean angular deviations and Jaccard indices of TSM/LLX embeddings versus those of SVD. These results were averaged over the five non-
11

---

Published in Transactions on Machine Learning Research (12/2022)
DIGITS (70K)
WORDS (108K)
TSMILLX
1
TTSSMILLX
-VDD
-SVD
0.8
0.8
JACCARD INDEX
0.6
0.6
0.4
0.4
0.2
0.2
0
0
16
64
128
256
4
16
64
128
256
dimensionality
dimensionality
Figure 7: Jaccard indices in eq. 12 for the embeddings from TSM/LLX and SVD of 70K digit images (left) and 108K word vectors (right), as a function of their dimensionality.
overlapping folds of each data set (where particular fold specified which m of examples were chosen as landmarks). Here again, we see that the TSM/LLX embeddings are much better as evaluated by these metrics. As to be expected, these TSM/LLX embeddings have somewhat smaller gains than the TSM embeddings in Figs. 3 and ! which do not involve an additional layer approximation. Nevertheless, the gains are still quite signíficant. Though not shown the figures, we also observed in these experiments that the results did not vary much across different folds.
4 Conclusion
In this paper have explored geometrical connection between sparse and low-rank matrices. The connec- tion emerged from study of low-dimensional embeddings that preserve the meaningfulness of norms and small angles. These embeddings were derived from the nonlinear matrix decomposition in eq. [ and they
were learned by solving sequence of cigenvalue problems to minimize the cost in eq. (4)
There are two clear directions for future work. The first is to investigale more scalable implementations of the algorithm in section 2 One such approach has been described section 3 but there are many
other (complementary) tools for this purpose, including the Nystrom approximation (Williams Seoger
2001
Drineas & Mahoney 2005 Gittens & Mahoney 2016), randomized algorithms (Mahoney 2010 Halko
et al. 2011 Tropp et al.
and alternating least-squares methods (Hastie et al. 2015 The second direction is to develop the ideas in this paper for deep learning, has been done for other models of
2017
manifold learning (Jia et al. 2015 Chen et al. 2018 Mishne et al. 2019 Du et al. 2021 Pai et al. 2022
In particular, we note that a simple modification of eq. (2) yields a parametric ReL.U model for learning
similarity- preeervip mappings x + y. Such a mapping is described by the implicit equation
max x- x x- +|| max(0. y.u, -- -
(16)
where the weight vectors μj and Dj prototypes that tile, respectively, the spaces of high-dimensional inputs and low-dimensional outputs (Sengupta et al.) 2018). Tt seems fruitful to investigate eq. (16)
a
layerwisc prescription for unsupervised lcarning in deep networks (Hinton & Salakhutdinoy 2006 Hinton
et al. 2006 Mazundar & Rawal, 2019b).
Acknowledgements
The author thanks the anonymous reviewers for their detailed and constructive feodback. He also acknow]- edges many illuminating discussions with S. Seung, K. Luther, R. Yang, A. Sengupta, D. Chklovski,, and
Y. Ma.
12

---

Published in Transactions on Machine Learning Research (12/2022)
References
A. Agrawal, A. Ali, and S. Boyd. Minimum distortion embedding. Foundations and Trends in Machine Learning Research, 14(3):211 378, 2021.
R. Barrett, M. Berry, T. F. Chan, J. Demmel, Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. van der Vorst. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. SIAM, 1994.
M. Belkin and P. Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural Computation, (5(6):1373-1396 2003.
T. Bolukbasi, K.-W. Chang, J. Zou, V. Saligrama, and A. Kalai. Man is to computer programmer as woman is to homemaker? debiasing word embeddings. In D. Lee, M. Sugiyama, U. Luxburg, Guyon, and R. Garnett (eds.), Advances in Neural Information Processing Systems 29, pp. 4349 357, 2016.
S. S. Chen, D. Donoho, and M. A. Saunders. Atomic decomposition by basis pursuit. SIAM Journal on Scientific Compuling, 20(1):33 61, 1998.
Y. Chen, D. M. Paiton, and B. A. Olshausen. The sparse manifold transform. In S. Bengio, H. Wallach, Н. Larochelle, K. Grauman, N. Cesa Biancii and R. Garnett (eds.), Advances in Neural Information Processing Systems 31, 10534 10545. Curran Associates, Inc., 2018.
R. R. Coifman, Lafon, A. B. Lee, M. Maggioni, В. Nadler, F. Warner, and Zucker. Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps. Proceedings of the National Academy of Sciences, 102((1):7426 7431, 2005.
M. Collins, S. Dasgupta, and R. F. Schapire. A generalization of principal components analysis to the expo- nential family. In T. G. Dictterich, S. Becker, and Z. Ghahramani (eds.), Advances in Newral Information Processing Systems 14, pp. 617 624, Cambridge, MA, 2002. MIT Press.
T. Cox and M. A. Cox. Multidimensional scaling. Chapman and Hall /RCC, 2000.
S. Deerwester, S. T. Dumais, G. W. Furnas, T. K. Landauer, and R. Harshman. Indexing by latent semantic analysis. Journal of the Association for Information Science and Technology, 41(6):391-400 1990.
A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion). Journal of the Royal Statistical Society B, 39:1 38, 1977.
L. Donoho and C. Grimes. Hessian eigemaps: locally linear embedding techniques for high-dimensional data. Proceedings of the National Academy of Sciences (USA), 100(10):5591-559)) 5596, 2003.
P. Drineas and M. W. Mahoney. On the Nyström method for approximating Gram matrix for improved kernel-based learning. Journal of Machine Learning Research, 6(72):2153 2175, 2005.
Y. Du, K. Collins, J. Tenenbaum, and V. Sitzmann. Learning signal-agnostic manifolds of neural fields. In M. Ranzato, A. Beygelzimer, Y. Dauphin, S. Liang, and J. W. Vaughan (ods.), Advances in Neural Information Processing Systems 34, pp. 8320 8331. Curran Associates, Inc., 2021.
S. Duff, A. M. Erisman, and K. Reid. Direct Methods for Sporse Motrices. Oxford University Press, 2017.
C. Eckart and G. Young. The approximation of one matrix by another of lower rank. Psychometrika, 1: 211 218, 1936.
Gillas. Nonnegative motria factorization. SIAM, 2021.
Alex Gittens and Michael W. Mahoney. Revisiting the nyström method for improved largo-scale machine learning. Journal of Machine Learning Research, 17(117):1-655 2016.
13

---

Published in Transactions on Machine Learning Research (12/2022)
H. Gonen and Y. Goldberg. Lipstick pig: Debiasing methods cover up systematic gender biases in word embeddings but do not remove them. In Proceedings of the North American Associalion for Compulational Linguistics (NAACL-19), pp. 609 614, 2019.
N. Halko, P. G. Martinsson, and J. A. Tropp. Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2)2217-288, 2011.
T. Hastie, R. Mazumder, J. D. Lee, and R. Zadeh. Matrix completion and low-rank SVD via fast alternating least squares. Journal of Machine Learning Research, 16(104):3367 3402, 2015.
G. E. Hinton and Z. Ghahramani. Generative models for discovering sparse distributed representations. Philosophical Transactions of the Royal Society of London B, 312(13588:1117-11990, 1997.
G. E. Hinton and R. R. Salakhutdinov. Reducing the dimensionality of data with neural networks. Science, 313(5786):504 507, 2006.
G. E. Himton, S. Osindero, and Y.-W. Teh. A fast learning algorithm for deep belief nets. Neural Computation, 18(7):15771154, 2006.
G. E. Hinton, A. Krizhovsky, and S. D. Wang. Transforming auto-encoders. In Proceedings of the Interna- tional Conference on Artificial Newral Networks (ICANN-11), pp. 44 51, 2011.
P. D. Hoff. Bilinear mixed-effocts models for dyadic data. Journal of the American Statistical Association, (00:286-295, 2005.
R. A. Horn and C. R. Johnson. Matric Analysis. Cambridge University Press, 2012.
P. Jaccard. Étude comparative de la distribution florale dans une portion des alpes et des jura. Bulletin de la Société Vaudoise des Sciences Naturelles, 37(142):547-579, 1901.
K. Jia, L. Sun, S. Gao, 7. Song, and B. E. Shi. Laplacian auto-encoders: explicit learning of nonlinear data manifold. Neurocomputing, 160:250 260, 2015.
Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied document recognition. Proceedings of the IEEE, 86(11):2278-2324, 1998.
D. D. Lee and H. S. Scung. Learning the parts of objects by non-negative matrix factorization. Nature, 401: 788 791, 1999.
D. D. Lee and Н. Sompolinsky. Learning a continuous hidden variable model for binary data. In M. J. Kearns, S. A. Solla, and D. A. Cohn (eds.), Advances in Neural Information Processing Systems 11, pp. 515 521, Cambridge, MA, 1999. MIT Press.
M. Levandowsky and D. Winter. Distance between sets. Nature, 234(5):34 35, 1971.
G. C. Linderman, M. Rachh, J. G. Hoskins, S. Steinerberger, and Y. Kluger. Fast interpolation-basee t-SNE for improved visualization of single-cell RNA-seq data. Nature Methods, 16:243-244 245, 2019.
M. W. Mahoney. Randomized algorithms for matrices and data. Foundations and Trends in Machine Learning, 3(2):123 224, 2010.
J. Mairal, F. Bach, J. Ponce, and G. Sapiro. Online learning for matrix factorization and sparse coding, Journal of Machine Learning Research, 11:19-60, 2010.
A. Mazumdar and A. S. Rawat. Learning and recovery in the ReLU model. In Proceedings of the 57th Anaual Allerton Conference Communicalion, Control, and Compuling, pp. 108 115, 2019a.
A. Mazundar and A. S. Rawat. Learning network parameters in the ReLU model. In Proceedings of the NewrIPS 2019 Workshop on Solving Inverse Problems with Deep Networks, 2019b.
L. Mclnnes, Healy, and Melville. UMAP: uniform manifold approximation and projection for dimension reduction. ArXiv e-print 1802.03426, 2018.
14

---

Published in Transactions on Machine Learning Research (12/2022)
G. J. McLachlan and T. Krishnan. The EM Algorithm and Eatensions. Wiley, 2008.
X.-L. Meng and D. Dyk. The FM algorithm the Royal Statistical Society B, 59(3):511 567, 1997.
old algorithm sung to fast tune. Journal of
J.-B. Michel, Y. K. Shen, A. P. Aiden, A. Veres, M. K. Gray, J. P. Pickett, D. Hoiberg, D. Clancy, P. Norvig, J. Orwant, S. Pinker, M. A. Nowak, and E. L. Aiden. Quantilative analysis of culture using millions of doi: /..116/scienne..
digitized books. Science, (81((11)) 182, 2011.
URL https://www
science org. htt//do/aabs/10/ wwt.:/111sscccccco science. 1199644
T. Mikolov, K. Chen, G. Corrado, and J. Dean. Efficient estimation of word representations in vector space. In Proceedings of the First International Conference on Learning Representations (ICLR-13). 2013a. URL http://arxix htt.aaawww....iaaaaaaaaa/aaeeeea aba/ 1301.3781
Mikolov, Sutskever, K. Chen, G. S. Corrado, and J. Dean. Distributed representations of words and phrases and their compasitionality. In C. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger (eds.), Advances in Neural Information Processing Systems 26, pp. 3111 3119. Curran Associates, Inc., 2013b.
T. Mikolov, E. Grave, P. Bojanowski, C. Puhrsch, and A. Joulin. Advances in pre-training distributed word representations. In N. Calzolari, K. Choukri, C. Cieri, T. Declerck, S. Goggi, K. Hasida, H. Isahara, B. Maegaard, J. Mariani, H. Mazo, A. Moreno, J. Odijk, S. Piperidis, and T. Tokunaga (eds.), Proceedings of the Eleventh International Conference on Longuage Resources and Evaluation (LREC-18). European Language Resources Association, 2018.
G. Mishne, U. Shaham, A. Cloninger, and 1. Cohen. Diffusion nets. Applied and Computational Harmonic Analysis, 47(2):259-22) 285, 2019.
V. Nair and G. E. Hinton. Rectified linear units improve restricted Boltzmann machines. In Proceedings of the 27th International Conference on Machine Learning (ICML-10), pp. 807-814, 2010.
M. E. J. Newman. The structure and function of complex networks. SIAM Review, 45:167-256 256, 2003.
B. A. Olshausen and D. J. Field. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature, (65833))) 709,, 1996.
G. Pai, A. Bronstein, R. Talmon, and R. Kimmel. Deep isometric maps. Image and Vision Computing, pp- 104461, 2022.
B. T. Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4(5):1-17, 1964.
Radford, L. Metz, and S. Chintala. Unsupervised representation learning with deep convolutional gen- erative adversarial networks. Y. Bengio and Y. LeCun (eds.), Proceedings of the Fourth International Conference on Learning Representations (ICLR-16). 2016. URL http:://arxiv..rg/bbbba
06434
S. T. Roweis and T.. K. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science, 290: 2323 2326, 2000.
R. Rubinstein, A. M. Bruckstein, and M. Elad. Dictionarics for sparse represcntation modeling. Proceedings of the IEEE, 98(6):1045-1088 2010.
L. K. Saul. A nonlinear matrix decomposition for mining the zeros of sparse data. STAM Journal on Mathematics of Data Science, 4(2):431 463, 2022.
A. M. J. Schakel and B. J. Wilson. Mcasuring word significance using distributed representations of words. Arxiv e-print 1508.02297, 2015.
15

---

Published in Transactions on Machine Learning Research (12/2022)
A. M. Sengupta, C. Pehlevan, M. Tepper, A. Genkin, and D. B. Chklovskii. Manifold-tiling localized receptive fields are optimal in similarity-pprserviin neural networks. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi and R. Garnett (eds.), Advances in Neural Information Processing Systems 31, pp. 7080-7090. Curran Associates, Inc., 2018.
H. S. Seung and D. D. Lec. The manifold ways of perception. Science, 290(3500):2268 221222226229,, 2000.
J. B. Tenenbaum, V. de Silva, and J. C. Langford. A global geometric framework for nonlinear dimensionality reduction. Science, 290:2319-2533 2323, 2000.
J. A. Tropp, A. Yurtsever, M. Udell, and V. Cevher. Practical sketching algorithms for low-rank matrix approximation. SIAM Journal on Malriz Analysis and Applicalions, 38(4):1454 11001145--185,, 2017.
M. Turk and A. Pentland. Eigenfaces for recognition. Journal of Cognitive Neuroscience, 3:71-86, 1991.
M. Udell, C. Horn, R. Zadeh, and Boyd. Generalized low rank models. Foundation and Trends in Machine Learning, 9(1):1-118 2016.
L. van der Maaten. Accelerating t-SNE using tree-based algorithms. Journal of Machine Learning Research, 15:3221 15:3221-3245, 3245, 2014.
L. van der Maaten and G. Hinton. Visualizing data using t-SNE. Journal of Machine Learning Research, 9: 2579 2605, 2008.
M. Vladymyrov and M. Á. Parreira-Perpinán. Locally linear landmarks for large-scale manifold learning. In Proceedings of the European Conference on Machine Learning and Knowledge Discovery in Databases (ECML/PKDD-B)) pp. 256-271, 2013.
K. Q. Weinberger and L. K. Saul. An introduction to nonlincar dimensionality reduction by maximum variance unfolding. In Proceedings of the Twenly Pirst Nutional Conference on Artificial Intelligence (AAAL-06), pp. 1683 1686, 2006.
K. Q. Weinberger, B. Packer, and L. K. Saul. Nonlinear dimensionality reduction by semidefinite program- ming and kernel matrix factorization. In Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics, pp. 381 388, 2005.
C. K. I. Williams and M. Seeger. Using the Nyström method to speed up kernel machines. In T. Leen, T. Dietterich, and V. Tresp (eds.), Advances in Neural Information Processing Systems 13, pp. pages 682-66 688. MIT Press, 2001.
J. Wright and Y. Ma. High- Digh-Dimensi Data Analysis with Now-Dimensional Models: Principles, Compu- tation, and Applications. Cambridge University Press, 2021.
R. Zemel and M. Carreira-Perpiñáán. Proximity graphs for clustering and manifold learning. In T.. Saul, Y. Weiss, and L. Bottou (eds.), Advonces in Neural Information Processing Systems 17, pp. 225 232, 2004.
X. Zhu, Z. Ghahramani, and J. Lafferty. Semi-supervised learning using Gaussian fields and harmonic func- tions. Proceedings of the Twenlüeth International Conference on International Conference on Machine Learning (ICML-03), pp. 912 919, 2003.
16

---

Published in Transactions on Machine Learning Research (12/2022)
A Supplementary details for alternating minimization
In this appendix we provide further details of the alternating minimization for TSM in section 2.2 specifically its initialization, regularization, and convergence.
A.1 Initialization
As discussed in section 2.2 the results from the alternating minimization of eq. 4 can depend on how the procedure is initialized. For the results this paper, we adopted a simple approach, initializing this procedure from the baseline embeddings computed by SVD. We describe this initialization next.
To perform the updates in eqs. 5 it is necessary to initialize either the matrix L or the matrix Z. We initialized the matrix L in the following way. Let d denote the desired rank of the matrix L. We began by setting
Lij U. By 1 e/ |lyl
(17)
where o, are the d-dimensional projections the imputs from eq. (10). Note that this initialization generally produces matrix L of rank d+1;; 1; nevertheless, it seemed to work well in practice. The minimization of
eq. { proceeds by updating Z from eq. (), then updating L from eq. 6 and so on. The rank of is
therefore equal to d after all subsequent updates.
As in all non-convex optimizations, it is problematic when results are sensitive to different initializations. In this respect, TSM similar to other popular algorithms for nonlinear dimensionality reduction, such as t-SNE fvan der Maaten & Hinton 2008) and UMAP (Melnnes et al.] 2018), whose solutions also depend on how their non-convex optimizations are initialized. These approaches generally achieve their best results by employing specialized heuristios, such early exaggeration, by initializing their optimizations from the bottom eigenvectors of graph Laplacian.
A.2 Regularization
There can be other ways besides initialization to address the problem of local minima in non-convex opti- mizations. For TSM, we imposed an additional linear constraint on the matrix to guide the optimization better minima. This constraint was based on the observation that distances between dissimilar inputs tend to increase when they are faithfully embedded into a lower-dimensional space (Weinberger & Saul 2006). Based on this intuition, we sought to encourage embeddings æ; y: satisfying
Σ 1j
y: - 1
Σ æ.- - |0 ||| ||=
(18)
Eq. (18) states that pairs of outputs are on average further apart in angle than their corresponding pairs of inputs. We encouraged solutions with this property by identifying the left side of eq. (18) with the elements of Z and constraining their overall sum. In particular, we replaced eq. 4 with the constrained optimization
(L) = d.
min ||
||| such that max (0, Z) S. Σ. Ziz 2
(L-ZZ)
(19)
- T T||0;||| |||,,1111111
The original update for requires only a slight modification to enforce this extra constraint. To proceed, rewrite the original update in eq. (5) with placcholder variable Z. In particular, let
2, =
il Sij > 0, min Lis) otherwise.
Syy
(20)
Additional coustraints narrow the scarch space, and thus in general they may restrict solutions minima. But well-ckosen constraint narrow the search space in way that guides non-convex optimizations better minima. This what appeared occur in practice.
17

---

Published in Transactions on Machine Learning Research (12/2022)
100
DIGITS WORDS
normalized loss
10
n=25k, d=32
n=14K, d=16
10
50
100
150
200
250
iteration
Figure 8: Behavior of the objective function in eq. 4 normalized by its initial value, over the course of 250 alternating updates to estimate and L. The updates were performed with the initialization, regularization, and acceleration described in eqs. 17 23
Only one further calculation required for the modified update. Let no denote the number of zero elements in S, and let 8 denote the amount, per zero clement of S, by which the original update for violates the new constraint in eq.
max
0, Σ
x.
=||
7 |,ا| 3.0])
(21)
no
Then in terms of these quantities, the update for the constrained optimization in eq. 19 takes the simple form:
if Sij > O, m0 Lij) + otherwise.
Sij
Zij =
(22)
Comparing eq. (20 and eq. [22 we see that the constrained update differs only by a uniform shift of those elements in Z that are not already determined by the nonzero elements of S.
A.3 Acceleration
The alternating updates in eq. (6) and eq. (22) are guarantoed to monotonically decrease the cost in eq. (99 at each step. This is especially useful guarantee for debugging first implementation of the algorithm. Tn practice, though, found it useful to modify the update for Z with additional momentum term (Polyak)
1964h. In particular, let Z, denote the value of Z after the 9th update, and let Z be the updated value that would have been computed originally from eq. (22 Instead, we set
Z = Z + <Z: -Z:-Zz-1)),
(23)
where 7 E (0,1) 1) is so-called momentum parameter. In practice, we found that the update in eq- 23 not
only converges more quickly, but also finds better solutions. All of the results in section 2 were obtained by setting == 0.9 and completing 250 of these alternating updates. Fig, g plots the objective function in eq. 44, normalized by initial value, over the course of these updates for the data sets section 2.3
Supplementary results for Isomap
Many algorithms have been proposod for learning nonlincar grighborhhood-reservvnn embeddings. These inchude several methods, such as Isomap [Tenenbaum 2000), Hessian LLE (Donoho & Grimes 2003
and maximum variance unfolding (Weinberger & Saul| 20061, that are similar in important respects to TSM.
18

---

Published in Transactions on Machine Learning Research (12/2022)
DIGITS (14K)
WORDS (25K)
TSM ISOMAP
TSM
ISOMAP
0.4
SVD
0.8
0.6
JACCARD INDEX
0.2
0.4
0.2
0
4
12 dimensionality
dimensionality
Figure 9: Jaccard indices in eq. 12 for the embeddings from TSM, Isomap, and SVD 14K digit images (left) and 25K word vectors (right), as function of their dimensionality.
In particular, the embeddings from these algorithms are designed (and in some settings, guaranteed) to preserve distances bet ween nearby imputs; they are also derived from the eigenvectors of specially constructed n x n matrices. It therefore instructive to compare embeddings with these properties to those of TSM.
In this appendix we compare embeddings from TSM and Isomap on the data sets of 14K images and 25K word vectors in section 2.3 We chose Isomap as benchmark because it does not share the numerical difficult ies of Hessian LLE for high dimensional data, and also because it scales better than maximum variance unfolding to data sets with large numbers of examples. It is important to emphasize that the embeddings of Isomap are not designed to preserve norms or angles between nearby inputs; this has important consequences, as we observe below. Nevertheless, one can still evaluate how well Tsomap succeeds in computing a neighborhood- preserving embedding when it operates on the same sparse neighborhood graph as TSM.
The Isomap algorithm consists of three steps. First, constructs sparse graph whose nodes represent inputs and whose edges (connccting ncarby inputs) are weighted by pairwise distances. Second, it uses dynamic programming to compute all shortest pairwise distances through this graph. Finally, from the n n matrix of pairwise distances, it uses multidimensional scaling to recover a oow-dimesiional embedding. More specifically, this is done by assuming that the oow-dimenioonal outputs are centered at the origin, then dcriving and diagonalizing their Gram matrix.
For the experiments this section, we seeded Isomap with a sparse graph whose edges were one-to-one correspondence with the nonzero elements of the sparse matrix S computed by TSM in eq. (3). In addition, weighted the edges of this graph by the Euclidean distances || between nearby normalized inputs, where è = æ //=//= The Isomap embeddings were therefore designed to co-locate pairs of outputs whose corresponding inputs have small cosine distances 11 < 1 Equivalently these are the pairs of inputs for which cos(æ;, >T
As before, we used a Jaccard index between sets Q2 and Gy measure the quality of the embeddings x; > Yr obtained in this way; see eq. 11). For Isomap, however, we constructed the set Sy in slightly different way; namely,
1-7). This was done to account for the fact that Isomap attempts to preserve distances rather than angles. Note that if the embedding by Isomap were perfectly neighborhodd-prsorvingg then this set Gy would contain exactly those pairs (i.j) that index nonzero elements of the matrix S in eq. 3 With this allowance, we can therefore compare the Jaccard indices from Isomap, TSM, and SVD in a meaningful way.
took Py (t,1)) 1 - ,, <
Figure g compares the Jaccard indices from these algorithms for embeddings of different dimensionality. There are two striking results. First, the embeddings from Isomap have higher (better) Jaccard indices than TSM low dimensionalities; second, after certain point, the higher-dimensional embeddings from Tsomap
19

---

Published in Transactions on Machine Learning Research (12/2022)
0.5
0.5
0
-0.5
-0.5
-1
-1
-0.5
0
0.5
-1
Figure 10: Two data sets that illustrate the additional requirements norm ceem-presevving embeddings. these data sets, the inputs (shown in red) all have unit norm, and the nearest neighbors are connected by black edges. The inputs can be faithfully mapped into one dimension if the only goal is preserve distances between nearest neighbors. But dimensionality reduction is possible if the mapping must also preserve the norms of the inputs.
yield progressively lower (worse) Jaccard indices seemingly counterintuitive result. Next we show how to
reconcile these observations.
First we consider why Isomap's embeddings have higher Jaccard indices for lower dimensionalities. One plausible explanation is that Isomap only aims to compute neighborhooddpreservii embedding; unlike TSM, it is not also required to compute norm-preserving embedding. The latter may be difficult impossible to accomplish when the target dimensionality is small. The additional requirements of norm- preserving embeddings are illustrated by the two panels of Fig. [0] Each panel shows data set whose inputs, all unit norm, are distributed like beads along a string with two open ends. In both cases, Isomap can faithfully embed the inputs into R²2 as uniformly spaced points along the real line. But while such mappings can preserve distances between nearest neighbors, they cannot also preserve distances the origin. Such effects may account for the lower Jaccard indices of TSM at lower dimensionalities, as its embeddings in this regime are subject to more stringent constraints.
Next we examine why Isomap's embeddings of these data sets appear to deteriorate- at least, measured by their Jaccard indices for higher dimensionalities. Here the most plausible explanation is that the data sets violate a basic assumption of the Isomap algorithm; the algorithm's guarantees hold only when the underlying manifold isometric to convex subset of Fuclidean space. In practice this breakdown signaled by the presence of negative cigenvalucs for the data sets of 14K digits and 25K word vectors. For each data set,, the first negative eigenvalue appears at dimensionality that slightly precedes the onset of decreasing Jaccard indices in Fig. 4 The negative eigenvalues also increase in frequency as the positive eigenvalues decrease in magnitude. This suggests a plausible explanation for the decreasing Jaccard indices of Isomap at higher dimensionalities.
the matrix that Isomap diagonalizes. Figure I plots these cigenvalues
It should be noted that we have applied Isomap in an unconventional way normalizing the inputs to unit length before computing distances- in order to permit a direct comparison to TSM. Recognizing this, we believe that the results in this section best serve to highlight the different requirements of norm- corm-preserving embeddings, well the different approaches needed to compute them.
20

---

Published in Transactions on Machine Learning Research (12/2022)
x10
DIGITS (14K)
x10
WORDS (25K)
positive negative
positive negative
6
X d
ISOMAP
EIGENVALUES
2
0
-2
15
20
d
Figure 11: Eigenvalues of the Gram matrix diagonalized by Isomap. Negative eigenvalues (shown in red) indicate that Isomap's graph distances do not correspond distances in any Euclidean space.
C
Supplementary details for LLX
Here we give further details to solve the sparse least-squares problem in section 3.1 As shorthand, let Y= [y y/2 yo] denote the outputs of the embedding, and let denote the n x n identity matrix. In this notation, we can rewrite eq. (15p
YY Trace (1--W]
(24)
To procoed, we distinguish the m known outputs Y (i.c., the landmarks) previously computed by TSM versus the n-m unknown outputs that remained be computed by LLX. Without loss of generality, and to simplify notation, assume that the data is ordered such that the m landmarks correspond to the last. examples of the data set. Then can write
{
Yu
y1 y2 Øn-m]
(unknown), (known).
Y Yu Yn where
(25)
Yz Un-m+1 Un-m+2
yn
We continue in this way by decomposing the x matrix Specifically, we write this matrix as
W inside the quadratic form of eq.
24
A E R (m)m))
A
B BI
B p)
1 1- W
where
(26)
C D
C (m-x(
D Rmxm
Note that on the right side of this equation, all of the submatrices A, B, C, D are sparse a property they inherit from the matrices and W on the left. With this notation, we can explicitly write out the lincar equations that locate the minimum of eq. 24h. In particular, we have:
AT -[ATB+CTD]
(27)
In this way obtain d systems of sparse lincar equations for the rows of Y... In practice, do not solve
these equations by inverting the (n-m) x (n-m) matrix on the left side of eq. (27 Instead, we solve each system independently using preconditioned conjugate gradients method (Barrett et al. 1994 This latter approach is much faster.
21


Reviews:
[{"review_id": "X5oldp9WF7", "reviewer": "Reviewer_WDzC", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\n## Strengths\n\n- The paper is very well written -- the method is described clearly and\n  completely, the work is very well situated with respect to prior work in the\n  various areas the work touches, and the experimental results are clearly\n  presented.\n- The approach described here has some interesting potential connections to\n  ReLU networks in deep learning, although the connection is left to future\n  work to be explored in detail.\n- The authors provide detailed guidance of various ways to tune the\n  hyperparameters (Figure 2, and around Step 3/Practical Considerations\n  sections on page 6), and initialize and speed up the algorithm (in Appendix\n  A).\n\n## Weaknesses\n\n- The method does not seem to be as robust as some other manifold learning\n  approaches. Although all these approaches have some amount of hyperparameter\n  tuning involved (e.g., the number of nearest neighbors or neighborhood size\n  to use in constructing adjacency graphs), the authors' method also has\n  some brittle-feeling aspects like the use of \"virtual inputs\" that\n  interpolate between inputs in the dataset to avoid having zero columns in\n  $\\boldsymbol{S}$ (bottom of page 4), the procedure used to factor the matrix\n  $\\boldsymbol{L}$ when the Gram matrix representation is invalid (Step 3 on\n  page 6), and the seeming importance of the initialization used for good\n  performance with the alternating minimization procedure (Appendix A). \n  The experimental evaluation of performance under these aspects of the method\n  (e.g., ablations) is mostly omitted.\n- In general, the experimental evaluation is limited. Since the method is\n  eventually a method for manifold learning, it would be very useful to compare\n  to other manifold learning methods (not just SVD, as is done here). This\n  makes it infeasible to evaluate the paper's contribution beyond simply\n  introducing a new method and showing that the associated nonconvex procedure\n  can be tuned to produce valid embeddings.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors apply an algorithm of Saul (2022) for computing a representation of\na sparse matrix as an elementwise thresholding (or a ReLU-ing, in deep learning\nparlance) of a matrix of (perhaps significantly) lower rank to the problem of\nnonlinear dimensionality reduction of a dataset in a high-dimensional space\n$\\mathbb{R}^D$. In this context, the primary contribution is to identify that\nSaul's algorithm can be applied in this setting (because the data matrix itself\nneed not be sparse), and to evaluate the suitability of the algorithm in this\nsetting. The authors consider data for which the norm is meaningful (e.g., word\nembeddings or images/image patches with varying illumination), and seek an\nembedding into a lower-dimensional space $\\mathbb{R}^d$ that preserves these norms,\nas well as local information (cosine similarities) and some global information\n(things that are far apart stay far apart). The embedding process has three\nstages. First, the data are used to construct a matrix $S_{ij} = \\mathrm{max}\n\\\\{ \\langle \\boldsymbol{x}_i, \\boldsymbol{x}_j \\rangle - \\kappa\n\\\\| \\boldsymbol{x}_i \\\\|\\\\| \\boldsymbol{x}_j \\\\|, 0 \\\\}$ of thresholded similarities;\nthis is similar to many procedures in manifold learning, but happens to be\nhomogeneous in the scale of the data. Next, the algorithm of Saul is used to\nobtain a rank-$d$ matrix $\\boldsymbol{L}$ which satisfies $\\boldsymbol{S}\n\\approx \\mathrm{max} \\\\{ \\boldsymbol{L}, \\boldsymbol{0} \\\\}$, via an alternating\nminimization procedure (nonconvex objective) with special initialization. \nFinally, the embedded data matrix $\\boldsymbol{Y}$ of $d$-dimensional vectors\nis computed from the matrix $\\boldsymbol{L}$ by a standard decomposition of an\nassociated Gram matrix. (From this description, one may infer that the\ncontribution of the authors' approach is to combine existing techniques into a\nnew approach to this problem.)\nThe hyperparameters $\\kappa$ and $d$ above should be tuned on the data.\nThe authors do not provide convergence guarantees for their nonconvex approach,\nbut rather demonstrate experimentally that it produces sensible embeddings,\nmeasured through mean absolute deviation of the angles of nearby points\n(quantifying the degree of preservation of the local input geometry) and\nthe Jaccard index of the supports of $\\mathrm{max}\\\\{\\boldsymbol{S},\n\\boldsymbol{0}\\\\}$ and $\\mathrm{max}\\\\{\\boldsymbol{L}, \\boldsymbol{0}\\\\}$\n(quantifying the quality of the solution to the nonconvex optimization\nprocedure of Saul described above), on an image dataset (MNIST) and a language\ndataset (pretrained word embeddings) versus an SVD-based baseline embedding.\nThe proposed method obtains embeddings whose degree of \"faithfulness\" to the\ninput geometry, measured by these metrics, grows more quickly as the embedding\ndimension increases than those of the SVD baseline. The authors discuss an\nincremental version of their algorithm in Section 3 that has the potential to\nscale to larger datasets (e.g., $n = 10^5$) with a mild loss of representation\nquality.\n\n", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "- Comparisons to other manifold learning methods (either small extensions to\n  other methods to allow them to operate in the present \"scale-homogeneous\"\n  setting, or the use of the proposed method in the setting where all data have\n  the same scale) seem essential to judge whether the proposed method is\n  suitable for use in other manifold learning applications, how the\n  representations it learns are different from those of other manifold learning\n  procedures, etc. The SVD baseline does not seem sufficient here, as it has no\n  ability to capture nonlinear low-dimensional structure in the input data.\n  This seems to be an essential addition for future work to build on this\n  submission.\n\n## Minor changes/clarifications\n\n- At the bottom of page 7, this error term is between matrices that do not have\n  the same size. Should it be $\\boldsymbol{U} \\boldsymbol{P}$ instead?\n- I was slightly confused reading how the Jaccard index is used for evaluation\n  on page 8. The authors refer to eqn (7) as the definition of the matrix\n  $\\boldsymbol{L}$, but this equation defines this matrix in terms of the\n  embeddings $\\boldsymbol{Y}$ that are constructed in the sequel rather than as\n  the output of \"Step 2\", and it is not clear to me that the \"if and only if\"\n  statements made at the end of the Jaccard index explanation hold when the\n  Gram factorization procedure needs to be done approximately, as described in\n  \"Step 3\" (in particular, when this occurs, eqn (7) only holds approximately,\n  if I follow correctly). Could the authors clarify what is being measured\n  here?\n\n", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "gSCT2BHO03", "reviewer": "Reviewer_ejvv", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n---\n\n1) I think the idea of building connections between sparse matrices and low rank matrices is very interesting. Both concepts are very important in learning low complexity representations. So developing methods that can convert one representation to the other is very interesting. \n\n2) The experiments are interesting in the following sense. To me whenever, we learn lower dimensional representations of data, the holy grail of statistics to preserve is some metric (need not be the Euclidean metric). In this case, however, the paper suggests, ignoring the Euclidean metric, but instead focusing on the cosine similarity measure and well as preserving the norms of the data points. \n\n**Weakness**\n---\n\nFor the me the major weakness revolves around theoretical considerations of the problem \n\n1) The problem is non-convex. While this is clear due to the constraints. This is not discussed until the appendix. I think a discussion about this would be needed in the paper. \n\n2) Convergence of the method. I think some more work needs to be done to talk about the convergence of the method. Specifically, at the bottom of page 5, it says \"By construction, each of the updates in eqs. (5–6) is guaranteed to reduce the value of the objective function in eq. (4). Hence they can be alternated until they reach a desired level of convergence.\" Either a proposition/theorem with this result or a reference would be good. \n\nIn particular, this method has the same flavor as many different projection based optimization techniques such a Dykstra method of alternate projections. While the space of rank at most d matrices is not convex, the projections (modulo having repeated singular values) is well defined. So some theory may apply to provide convergence (not necessarily to the optimal). \n\nAlso why look at $\\|L - Z\\|_F^2$ as the objective and not $\\|S - \\max(0,L)\\|_F^2$ as the objective. I feel like these might be equivalent. \n\n3) I think more discussion on Equation (1) would be appreciated. In particular, what does the approximation mean here and are there results formally connecting sparsity of $S$ to the rank of $L$. \n\n\n**Questions**\n---\n\n1) Do you know anything about the true optimal solution for the problem? Maybe some approach that looks at Lagrange multiplier might yield some results?\n\n2) While the metric is not objective for the reduction. How well does it perform with respect to metric preservation? Note SVD, might not actually do that well on metric preservation (even MDS might not, if we start with non-euclidean data, but here is all about euclidean data.)\n\n\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents a new optimization problem for learning low rank representations of data. In particular, the paper is interested in learning low rank representations of sparse data. The paper then presents a method for learning representations by attempting to solve the optimization problem. Finally, they demonstrate the usefulness of the learned representation on two datasets - MNIST and Word Embeddings. The paper demonstrates superiority to the SVD on two different statistics. ", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "I have two major change requests. \n\n1) I think convergence issues of the algorithm (weaknesses 1,2) should be talked about more. Hence expanding on this would be good. \n\n2) Providing more background on the equation $S \\approx \\max(0,L)$ would also be appreciated. ", "broader_impact_concerns": "No broader impact concerns. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "4Hy44B64w2", "reviewer": "Reviewer_KAJF", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I believe the paper is very clear, to the point, and (to a good extent) complete in its discussions, for the setup under study. However, I would still like to ask the authors to address the points I bring up, before publication, especially I would like to see clarifications on the *modeling assumptions* here; situations in which they are satisfied and places they show up in the paper (in the arguments, interpretations, etc). ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "For a set of $N$ high-dimensional points given as the columns of matrix $X\\in \\mathbb{R}^{d\\times N}$, define $\\bar{G}(X) \\coloneqq [ \\langle x_i , x_j \\rangle - \\kappa \\|x_i\\| \\|x_j\\| ]\\_{ij}$ and consider a *thresholded similarity matrix* given by $S_\\kappa(X) \\coloneqq \\max(0, \\bar{G}(X) ) $. \n\nThe authors propose to find low-rank matrices $L$ and $Y$ such that first, $S_\\kappa(X)  \\simeq_{\\rm ls} \\max(0, L)$, and second, $L \\simeq_{\\rm ls} \\bar{G}(Y)$, as a means for *embedding* $X$ into a lower-dimensional space (as $Y$).  Note that: \n\nConsidering $S \\simeq \\max(0, L)$, observe that ``the sparser the matrix S, the larger the space that can be explored to discover a matrix L of significantly lower rank.`` Therefore, for a given sparse matrix $S$, one could be hopeful that a low-dimensional matrix $L$ in the above exists, satisfying $S \\simeq \\max(0, L)$. The reasoning behind the ``second`` part in the above being possible is less complete. Note $L$ and $Y$ are not being estimated simultaneously, rather first $L$ then $Y$ from $L$. \n\nThe authors discuss the interrelated roles of $\\kappa$ and the rank chosen for the final solution when estimating $L$, but do not provide algorithmic proposals for their simultaneous selection; especially in relation to possible properties of data ($X$). Various useful discussions on the role of $\\kappa$ has been given. \n\nThe authors propose an alternating minimization procedure (non-convex, based on least-squares) for finding $L$ and a least-squares procedure for estimating $Y$ from $L$. The former does not have any formal guarantees on convergence or on the rate of convergence. \n\nThe authors provide useful numerical experiments, as well as an ad-hoc procedure (LLX) for applying the proposed method to large data sets (as the proposed method suffers from memory issues due to SVD for $N\\times N$ matrices). ", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "\n---\n---\n\nRemark 0: The descriptions and discussions in Section 2.1 (Model) are excellent; clear and complete. However, I have various comments about the summaries and interpretations of this model given before Section 2.1, as mentioned below; especially the content of page 2 starting from ``But eq. (1) by itself has several limitations, ...`` till the end of the page. \n\n\n---\n\n\nThe problem description in the abstract (``We consider when a sparse nonnegative matrix $S$ can be recovered from a real-valued matrix $L$ of significantly lower rank``) is rather vague and significantly more general than $S=\\max(0,L)$ in Eq 1, and more specifically than the setup studied in the paper in Eq 2; especially note ``can be recovered from``. In fact, if Eq 2 is the limit of the scope of the current study, I think it should be mentioned from the very beginning as the problem description. \n\n---\n\nInstead of ``Many problems in high dimensional data analysis can be formulated as a search for structure in large matrices``, and considering the two examples given right after it, I think it is more appropriate to say ``Many problems in high dimensional data analysis can be formulated as a search for large structured matrices.``\n\n---\n\nBottom of page 1: Considering ``The intuition behind this observation is illustrated in the left panels of Figure 1.`` I think the left and middle panels of Figure 1 provide an illustration for Eq 1 and not for the implication claimed in ``the sparser the matrix $S$, the larger the space ...``. For the latter, it would be informative to use a sparsified version of $S$ (in Fig 1) and provide a matrix $L$ of *least rank* (and the same for $S$ itself, to compare); where the *least rank* property has been insured by exhaustive search (going to a smaller dimension for $S$ might help with computations; both to find such a least rank matrix, and to possibly have a unique least rank matrix so that the comparison is possibly more informative.)\n\n---\n\nIn the middle of page 2, please elaborate on Item 2; if `sparsification' (by thresholding) of the true similarity matrix is meant (as discussed in the third paragraph of page 10), I think it should be stated as a step enabled by the properties of the data (e.g., see Remark 2 and Remark 3) rather than a general way of dealing with non-sparse data (of course identifying and exploiting such property of the data is a contribution, and I am pointing out that it would be helpful if it is clarified). \n\nIn addition, I am not sure if I understand what the authors mean by ``a geometrical picture`` for ``a dataset``, in Item 1. After the three items, could you clarify what the ``inputs`` mean in ``The inputs are not required to be sparse``?\n\n---\n\n\nAs discussed on top of page 4, in Section 2.1, the paper considers manifold learning when the origin and angles are meaningful. I suggest that this distinction with manifold learning in general be mentioned in lines 3-5 of the abstract (I understand that it is somewhat 'apparent' in the rest of the abstract), and especially in the paragraph on top of page 3 (the paragraph before last, in Section 1.)\n\n\n---\n\nIn the last paragraph of Section 1, in ``input similarities are encoded by a sparse matrix $S$ of thresholded inner products.`` I am not sure if ``input similarities are encoded by`` is exactly a correct description; maybe ``input data is in the form of thresholded input similarities`` instead? \n\n\n\n---\n\n\nIn the caption of Figure 2, adding *1NN' to ``nearest neighbors`` (other than the legend) might be helpful to the readers.\n\n\n\n---\n---\n\nLet me provide a definition, aligned with the goal of paper in studying manifold learning in setups where origin has a meaning (top of page 4 in Section 2.1). For any nonzero point $x\\in\\mathbb{R}^d$, define $\\bar{x}\\coloneqq x/\\|x\\| \\in S^{d-1}$ where $S^{d-1}$ is the unit sphere in $\\mathbb{R}^d$. \n\n**Definition 1:** Given a **similarity-encoding** function $P: \\mathbb{R}^{d}\\times \\mathbb{R}^{d} \\to \\mathbb{R}_+$, for each matrix $X\\in \\mathbb{R}^{d\\times N}$ with no zero columns, define a **$P$-refined similarity matrix** as\n$$\nS_P (X)  \\coloneqq  [ \\| x_i \\| \\cdot \\| x_j \\| \\cdot P( x_i , x_j  )  ]_\\{ij} \\in \\mathbb{R}^{N \\times N}\n$$\nHere $P$ indicates **a refined similarity** provided by **a pre-processing step**; e.g., `Step 1' of Section 2.2, or a 1-nearest-neighbor procedure discussed at the end of Section 1. \n\n**Example 1:** For a given $\\kappa\\in (0,1)$, consider $P^\\kappa$ defined as \n$$\nP^\\kappa ( x_i , x_j ) \\coloneqq \\max ( 0 , \\cos( x_i , x_j )  - \\kappa )  = ( \\cos( x_i , x_j )-\\kappa ) \\cdot \\delta\\_{ij} ,\n$$\nwhere $\\delta_{ij}$ is a 0/1 indicator function for whether $\\cos(x_i,x_j)$ exceeds $\\kappa$. \nNote that $P^\\kappa$ is the output of the pre-processing step in `Step 1' of Section 2.2. \n\n**Example 2:** For a given $\\kappa\\in (0,1)$, consider $P^{\\kappa, \\rm nn}$ defined as \n$$\nP^{\\kappa, \\rm nn}( x_i , x_j ) \\coloneqq (\\cos(x_i, x_j)-\\kappa) \\cdot p\\_{ij} \\,,\n$$ \nwhere $p_{ij}$ is a 0/1 indicator function for whether $x_j$ is the nearest neighbor to $x_i$ or not. Note that we are slightly abusing the notation, by allowing $P$ to have access to all of $X$. \n\n**Example 3:** For a given $\\kappa\\in (0,1)$, consider $P^{\\kappa, \\rm deb}$ defined as \n$$\nP^{\\kappa,\\rm deb} ( x_i , x_j ) \\coloneqq \\cos(x\\_i, x\\_j) \\cdot \\delta\\_{ij} \\,,\n$$\nwhere $\\delta_{ij}$ is a 0/1 indicator function for whether $\\cos(x_i,x_j)$ exceeds $\\kappa$. Note that $P^{\\kappa,\\rm deb}$ is a ``debiased`` version of $P^\\kappa$ in Example 1. \n\n\n---\n\nRemark 1: There are two separate notions that are being used interchangeably in parts of the text. \nOne is the **thresholded similarity matrix** given by $S_{P^\\kappa}$, \nand the other is the **local similarity matrix** given by $S_{P^{\\kappa, \\rm nn}}$. \nI think it would be helpful to distinguish between these two and refer to the appropriate one when needed. I will refer to this remark throughout my review. \n\n---\n\n``The algorithm regards the inputs $x_i$ and $x_j$ as similar whenever the cosine of the angle between them exceeds some threshold`` -- any relations (setup or method) to TSC (Thresholded Subspace Clustering; arXiv:1303.3716)?\n\n---\n\nTo bring up a possibility: Thresholded (an operation, namely $S_{P^\\kappa}$ as defined in Remark 1, and as performed in Eq 2) Similarity Matching in stead of Sparse (a property, which is not present for neither $x_i$'s or $y_i$'s, unless an operator acts on them) Similarity Matching. \n\n---\n\nNote: in this notation, the use of $S_{P^\\kappa}$ and the ``de-biased`` version $S_{P^{\\kappa, \\rm deb}}$ seem to be equivalent.\n\n---\n---\n\n\n\nRemark 2: \nConsider Remark 1. With $\\kappa$ chosen based on a `\\%99' threshold, for the two specific datasets mentioned here, $S_{P^\\kappa}$ is approximately the **weighted 1-NN graph** of the original data (i.e., $S_{P^\\kappa}(X) \\simeq S_{P^{\\kappa, \\rm nn}}(X)$), and the proposed method, namely \n\n> ``(Task 1:) seeking a low-dimensional embedding of the given high-dimensional data whose `thresholded similarity matrix' matches that of the original data``, \n\nis close to \n\n> ``(Task 2:) seeking a low-dimensional embedding of the given high-dimensional data whose `weighted 1-NN graph' matches that of the original data``. \n\nCould you please discuss \n\n- (Q1) whether other approaches for `Task 2` exist in the literature, \n\n- (Q2) what are the corresponding *geometric insights* behind those algorithms that might differ from the geometric picture (Eq 1, etc) behind your model (the three items in Section 2.1, or equivalently in Eq 2)? \nI am asking this because while the model of Section 2.1 for the two datasets and the choice of `\\%99` matches `Task 2` to a good extent, it would deviate from it in other scenarios (other data, other $\\kappa$), hence is expected to correspond to a different underlying *geometric picture* compared to other methods addressing `Task 2`. \n\n- (Q3) are there data and $\\kappa$ settings, of practical interest, where `Task 1` and `Task 2` diverge?\n\n---\n---\n\nRemark 3: \nThe following are some questions on the proposed method, at the end of Section 2.1, for choosing $\\kappa$. \n\n- (Q1) Given $N$ points, there are $N$ angles between each point and its nearest neighbor and $N(N-2)$ angles between each point and its non-nearest neighbors. As a side point, I suggest clarifying that the blue and the orange histograms in Figure 2 have been normalized `separately'. Do the authors have any arguments as to why these two separately-normalized histograms, whose underlying data have sizes of different orders of magnitude, are comparable in the sense that their dividing point can be meaningfully defined, depicted, or used?\n\n- (Q2) Do the authors believe that a threshold such as `\\%99` would generally be close to the dividing point of the two histograms? (say for data sets that have been already deemed `suitable' for the proposed method by confirming that their two histograms are mostly non-overlapping; such as the two datasets here)\n\n\n- (Q3) When a dataset has two nearly-non-overlapping histograms, can we say that the dataset is **easy-to-embed** in some sense and proceed with a different method than the proposed one which forms $S_{P^\\kappa}(X)$? For example, *does the information that the farthest nearest neighbor is in less distance than the closest second-nearest neighbor allow for a simpler embedding procedure?* (basically, I am curious about any contenders from the literature in such a `simple' setup.)\n\n- (Q4) Could you reconcile the other points in this remark with the proposal at the end of page 4; ``Note that if this situation occurs for more than a small fraction of inputs, then it suggests to use a smaller value of $\\kappa$ in eqs. (2–3).`` In other words, is there a conflict between a dataset satisfying the above and it requiring lots of `virtual points`, and if not, then for $\\kappa$ chosen using the `\\%99`-procedure for datasets such as those above, what would reducing $\\kappa$ mean in relation to the histograms.  \n\n- (Q5) Could you reconcile the other points in this remark with the interrelation of $\\kappa$ and $d$, discussed at the end of Section 2.2, especially in relation to the properties of the Gram matrix. (beyond the empirical approach in ``This question can be answered to some degree by empirical evaluations, to which we turn next`` and Section 2.3, and instead, possibly using the classification of regimes offered by the *non-overlapping* property or the `Task 1` $\\simeq$ `Task 2` property from the above). \n\n---\n--- \n\nSection 2.2, bottom of page 5: please provide a formal discussion on why the method converges; i.e., why the *reduction* is strict, etc. \n\n---\n\nSection 2.2, bottom of page 5: it would be helpful to clearly state that the authors do not have a guarantee on the rate of convergence and the proposed algorithm could take a long time to converge in some scenarios, in addition to pointing to convergence plots for the datasets and values of $\\kappa$ on which they have already experimented. \n\n\n--- \n\nSection 2.2, middle of page 6, on ``The above solution was based on an assumption ... When this is not the case``: I believe the non-existence happens if and only if $G$ in Eq 8 is not positive semidefinite, and, the result is not desirable when $G$ in Eq 8 has a high rank. Could you please discuss these situations further (linear algebraic and Eq 8; see my *Summary Of Contributions* and what is mentioned as ``second``), before moving on to a least-squares approach?  \n\n--- \n\nI think the nonlinearity present in the specification of the model makes SVD irrelevant as a baseline (as clear from Table 1). Possibly, existing algorithms for `Task 2` (see Remark 3, Q3) could be used as baselines. \n\n--- \n\nIn Section 3, how do we select the landmarks?\n", "broader_impact_concerns": "Probably a discussion considering the analogy completion part. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper considers the problem of learning a similarity-preserving mapping of high-dimensional data into a low-dimensional space. The problem is formulated as recovering a low-rank matrix whose positive elements match a sparse matrix. The authors propose an alternating minimization algorithm to solve the problem and present experiments on MNIST and pre-trained word vectors to demonstrate the performance for manifold embedding. 

The findings could be of interest to researchers in manifold learning, low-dimensional models, and beyond. All the reviewers acknowledged the contribution of the paper and the submitted recommendations are (at least leaning toward) acceptance. I agree with the reviewers to recommend accepting. 

There are some comments (e.g., from Reviewer WDzC) that could help improve the manuscript and the authors agree to incorporate them into the final version. The authors are advised to address them in the final version. 
Decision: Accept (poster)