Graph Learning via Spectral Densification
Authors: Zhuo Feng
Yongyu Wang
Zhiqiang Zhao
Venue:  (2021)
Keywords: Spectral Graph Theory
Undirected Graphical models
Gaussian Markov Random Fields

Abstract: Graph learning plays important role in many data mining and machine learning tasks, such as manifold learning, data representation and analysis, dimensionality reduction, data clustering, and visualization, etc. For the first time, we present a highly-scalable spectral graph densification approach (GRASPEL) for graph learning from data. By limiting the precision matrix to be a graph-Laplacian-like matrix in graphical Lasso, our approach aims to learn ultra-sparse undirected graphs from potentially high-dimensional input data. A very unique property of the graphs learned by GRASPEL is that the spectral embedding (or approximate effective-resistance) distances on the graph will encode the similarities between the original input data points. By interleaving the latest high-performance nearly-linear
time spectral methods, ultrasparse yet spectrally-robust graphs can be learned by identifying and including the most spectrally-critical edges into the graph. Compared with prior state-of-the-art graph learning approaches, GRASPEL is more scalable and allows substantially improving computing efficiency and solution quality of a variety of data mining and
machine learning applications, such as manifold learning, spectral clustering (SC), and dimensionality reduction.

Under review as a conference paper at ICLR 2021
GRAPH LEARNING VIA SPECTRAL DENSIFICATION
Anonymous authors
Paper under double-blind review
ABSTRACT
Graph learning plays an important role in many data mining and machine learn- ing tasks, such as manifold learning, data representation and analysis. dimen- sionality reduction, data clustering, and visualization, etc. For the first time, we present a highly-scalable spectral graph densification approach (GRASPEL) for graph learning from data. By limiting the precision matrix to be graph- Laplacian-like matrix in graphical Lasso, our approach aims to learn ultra-sparse undirected graphs from potentially bigh-dimensional input data. A very unique property of the graphs learned by GRASPEL is that the spectral embedding (or approximate effective-resistance)) distances on the graph will encode the simi- larities between the original input data points. By interleaving the latest high- performance nearly-linear time spectral methods, ultrasparse yet spectrally-robus: graphs can be learned by identifying and including the most spectrally critical edges into the graph. Compared with prior ofate-of-thh-aar art graph learning ap- proaches, GRASPEL is more scalable and allows substantially improving com- puting efficiency and solution quality of variety of data mining and machine learning applications, such as manifold learning, spectral clustering (SC), and di- mensionality reduction.
INTRODUCTION
Graph learning is playing increasingly important roles in many machine learning and data mining applications. For example, key step of many existing machine learning methods requires convert- ing potentially nigh-dimensional data sets into graph representations: represent each (high-dimensional) data point as a node, and assign each edge a weight to encode the similarity between the two nodes (data points). The constructed graphs can be efficiently leveraged to represent the underlying structure of a data set or the relationship between data points (Jebara et al. 2009 Maier et al.) 2009. Liu et al. 2018). However, how to learn meaningful graphs from large data set at scale still remains a challenging problem.
a common practice to
Several recent graph learning methods leverage emerging graph signal processing (GSP) techniques for estimating sparse graph Laplacians, which show very promising results (Dong al. 2016 Egilmez et al. 2017] et al. 2019 Kalofolias Perraudin, 2019). For example, (Egilmez et al, 2017) addresses the graph learning problem by restricting the precision matrix to be graph Lapla- cian and maximizing a posterior estimation of attractive Gaussian Markov Random Field (GMRF) [ while an (I-reuularization term is used to promote graph sparsity; (Rabbat) 2017) provides error analysis for inferring sparse graphs from smooth signals; (Kalofolias & Perraudin, 2019) leverages approximate nearest-neighhor (ANN) graphs to reduce the number of variables for optimization; (Kumar al. 2019) introduces graph Laplacian learning method by imposing Laplacian spectral constraints. However, even the state-of-the-ae not scale well for large data set due to their extremely high algorithm complexity. For example. solving the optimization problem for Laplacian estimation in (Dong et al.) 2016 Kalofolias] [2066] Egilmez et al., 2017: Dong et al., 2019) requires O(N2) time complexity per iteration for N data entities and nontrivial parameters tuning for controlling graph sparsity which limits their applica- tions to only very small data sets (e. g. with up to a few thousands of data points); the method introduced in (Carey 2017) leverages Isomap manifold embedding (Tenenbaum et al., 2000 for graph construction, which requires ( time for manifold construction and thus does not scale to large data set; the latest graph learning approach (Kalofolias & Perraudin 2019) takes advantages of ANN graphs, but still runs very slowly for large data sets; the Laplacian estimation method with
Laplacian estimation methods for graph learning do
'If' the precision matrix of GMRF is M-matrix with all non-negative off-diagona elements, call it attractive GMRF (Slawski & Hem 2015 Dong al. 2019.

---

Under review as a conference paper at ICLR 2021
spectral constraints requires a good graph structure to be provided in advance (Kumar et al. 2019)
which otherwise can be very costly when going through exhaustive graph structure searches.
This work for the first time introduces a spectral graph densification approach (GRASPEL) for learning ultra-sparse graphs from data by leveraging the latest results in spectral graph theory (Feng)
2016 2018: Zhao et al. 2018). GRASPEL has close connection with prior GSP-based Lapla-
cian estimation methods (Dong et al., 2016; Kalofolias) 2016, Egilmez et al., 2017 Kalofolias & Perraudin 2019 Dong et al., 2019). and the graphical Lasso method (Friedman et al. 2008). By treating M -dimensioonal data points as M graph signals, GRASPEL allows elficiently solving con- vex problem by iterative identifying and including the most spectrally-critical edges into the latest graph leveraging recent nearly-lined Compared with prior spectral graph sparsification algorithms (Spiclman & Srivastava 2011 Feng 2016) that aim remove edges from a given graph while preserving key graph spectral properties, GRASPEL aims to add edges into the graph such that the learned graphs will have spectral embed- ding (or effective-resistance distances encoding the distances between the original input data points. Comparing with state-of-the-att graph learning methods, GRASPEL is more scalable for estimation of attractive Gaussian Markov Random Fields (GMRFs) for even very large data set. We summarize the contribution of this work as follows:
time spectral methods (Fengl 2016; 2018: Zhao et al. 2018).
We propose a spectral graph densification approach (GRASPEL) that allows efficient estimation of attractive Gaussian Markov Random Fields (GMRFs) leveraging the latest spectral graph theory.
We show that the graphical Lasso problem with Laplacian-like precision matrix can be efficiently solved by including crritrall--riticall edges to dramatically reduce spectral embedding distortions.
The key to achieving high efficiency is spectral emhedding scheme for finding spectrally-critical edges, allowing each GRASPEL itcration to be completed in O(N log N) instead of O(N²) time.
For the first time, we introduce a novel convergence criterion for graph learning tasks based on graph spectral stability: when the maximum embedding distortion becomes small enough, or equivalently when graph spectra become sufficiently stable, GRASPEL iterations can be terminated.
Our experiment results show that the graphs learned from high-dimensional data using GRASPEL can lead to more efficient and accurate spectral clustering (SC) as well as dimensionality reduction.
BACKGROUND OF GRAPH LEARNING VIA LAPLACIAN ESTIMATION
Given M observations on N data entities in a data matrix X [21 ..., Eygl E R N each column vector of X can be considered as a signal on a graph. For example, the USPS data set including 9, 298 images of handwritten digits with each image having 256 pixels will result in feature matrix X E RN X with N = 9, 298 and M = 256. The recent GSP-based graph learning methods (Dong al. 2016 estimate graph Laplacians from X for achieving the following desired characteristics:
Smoothness of Graph Signals. The graph signals corresponding to the real-world data should be sufficiently smooth on the learned graph structure: the signal values will only change gradually across connected neighboring nodes. The smoothness of a signal æ over an undirected graph (V, E, w) can be measured with the following Laplacian quadratic form:
x Lx =
Σ (P.))E
Wp.g(1 r (p) x (4))²
where = D W denotes the Laplacian matrix of graph G with D and W denoting the degree and the weighted adjacency matrices of G, and Wp.q denotes the weight for edge (p, q). The smaller value of (I) indicates the smoother signals across the graph. To quantify the smoothness (Q) of set of signals X over graph G, the following matrix trace can be computed (Kalofolias 2016):
Q(x, L.) = Tr(X LX),
where Tr( denotes the matrix trace.
Sparsity of the Estimated Graph (Laplacian). Graph sparsity is another critical consideration in graph learning. One of the most important motivations of learning graph is to use it for downstream data mining or machine learning tasks. Therefore, desired graph learning algorithms should allow better capturing and understanding the global structure (manifold) of the data set, while producing
2

---

Under review as a conference paper at ICLR 2021
sufficiently sparse graphs that can be easily stored and efficiently manipulated in the downstream algorithms, such as graph clustering, partitioning, dimension reduction, data visualization, etc. To this end, the graphical Lasso algorithm (Friedman et al. 2008) has been proposed to learn the struc- ture in an undirected Gaussian graphical model using /1 regularization control the sparsity of the precision matrix. Given a sample covariance matrix S and regularization parameter 8, graphical Lasso targets the following convex optimization task:
max log det(@) - Tr(eS) - 3 ||
(3)
e
over all non-negative definite precision matrices O The first two terms together can be interpreted the log-likelihood under Gaussian Markov Random Field (GMRF). || denotes the entry- 11 norm, so B||@||, becomes the sparsity promoting regularization term. This model tries to learn the graph structure by maximizing the penalized log-likelihood. When the sample covariance matrix S obtained from M i.id independent and identically distributed) samples X [21,..., FM where X ~ N (0, S) has an N N-dimensional Gaussian distribution with zero mean, each element in the pre- cision matrix L, encodes the conditional dependence between variables X; and Xj. For example, O = implies that the corresponding variables X, and X, are conditionally independent, given the rest. However, the log-determinant problems are very computationally expensive. The emerg- ing GSP-based methods infer the graph by adopting the criterion of signal smoothness (Kalofolias) 2016 Dong et al. 2016 Egilmez et al. 2017] Kalofolias & Perraudin, 2019). However, their ex- tremely high complexities do not even allow for learning graphs with even a few thousands of nodes Furthermore, these methods require nontrivial parameters tuning for achieving good performance.
GRASPEL: GRAPH SPECTRAL LEARNING AT SCALE
Similar to recent GSP-based Laplacian estimation methods, GRASPEL aims to more efficiently
solve the following convex problem for estimation of attractive GMRFs (Dong et al. 2019 Lake & Tenenbaum 2010 that is also similar to the graphical Lasso problem (Friedman et al. 2008):
1 Tr(X TT(XTex) ex) - 0|
max: log det
M
where θ - = 1,,, L denotes the set of valid graph Laplacian matrices, / denotes the identity matrix, and 0² > denotes prior feature variance. It can be shown that the three terms in 4 are corresponding to log (e). Tr(eS) and @||@|| in 3 respectively. When each column vector in the data matrix X is treated as graph signal vector, there is a close connection between formulation and the graphical Lasso problem. Since e = L + / I correspond to symmetric and positive definite (PSD) matrices (or M matrices) with non-positive off-diagonal entries, this formu- lation will lead to the estimation of attractive GMRFs (Dong et al. 2019:; Slawski & Hein 2015). In case that X is non- Gaussian, (4) can be understood Laplacian estimation based on minimizing the Bregman divergence between positive definite matrices induced by the function θ + log det (9)
(llawski & Hein 2015).
3.1 THEORETICAL BACKGROUND
Express the Laplacian matrix
L
Σ (P.q)E
wp.cep.00p.
(5)
where Ep E R denotes the standard basis vector with all zero entries except for the p-th entry being 1, and ep.q = Ep- Eg. Consider the objective function in (
F
log (
Tr(X exx B1 11,,
(6)
M
which can be further simplified as follows by substituting 5 into C
F Σ log(A: + 1/02)
1
Tr(X 0²
"p cmalli
43
Σ
Wpar (7)
M
(F-G)ER
(F-2)EF
For each of the N row vectors :()::) E R x where 1,.., N. the following two-step data pre- processing will be performed: (1) ((i,)) ((i:)) ns, where μa, denotes the sample mcan of (i, then
8 xx which results in Tr(OS) # Tr(X eX). (2) Feature normalization by X X/||
3

---

Under review as a conference paper at ICLR 2021
where the Laplacian eigenvectors corresponding to the ascending eigenvalues X denoted by uu for I,..., N, satisfying:
Lus A,U;
(8)
Taking the partial derivative with respect to the weight upo of edge (p, () leads to:
OF
N
X T ep,91é 43, M
1
(9)
1=2 X 1/02 Prcp.
Since the last two terms in are all fixed (constant) values for a given data matrix X where B can be considered as additional offset added to all data pairs (candidate edges), can drop the third term by simply setting = O, which will not impact the ranking of candidate edges in graph learning. The above simplification implies the second term alone will effectively penalize graph density for estimating Laplacian-like like precision matrix: including more edges will result in greater trace Pr(XTeX) The spectral perturbation analysis in Theorem 1 (Appendix) will lead to:
(u. (p.a)
(10)
00p.
we construct a subspace matrices U, for spectral graph embedding by using the first r- weighted nontrivial Laplacian cigenvectors as follows:
U2
Ury
U+
(11)
1e
1/022
VAr 1/02
then by setting 3 -0, can be rewritten as follows:
oF
1 M
1
10,1 enall
/ ||x ppall2 =
.emb
data "P.
(12)
@wp.p
=p.q
M
||,," Cp.g 12 and Fp. ||X ep.gll denote the &g distances in the spectral embed- ding space and the original data space, respectively. We note that using larger r for constructing U, will allow more accurate estimation of aF but at a higher computational cost.
where zp.4
awy.."
3.2 GRAPH LEARNING VIA SPECTRAL DENSIFICATION
OF value, it is considered as spectrally-critical If a candidate edge (p.g) 9) has relatively large Swp. edge, since will significantly impact and have large spectral embedding distortions defined as
rmmb p.g
Up.
M
-data
(13)
pp.4
Subsequently, 12 can be further simplified as follows:
OF
(1 (1-
emb
(14)
Owpa
Tpp.
"p.q
which implies that including the candidate edges with greater 7p.g and smmb values into the latest graph will allow a faster convergence of (4) in gradient descent (GD), as long as 7p.4 >> holds.
The proposed approach. When a² + too and r + N, the spectral embedding distance
emb becomes the effectiiv-resistance distance Reff and the scaled spectral embedding distortion Wp.g Reff becomes the edge leverage score for spectral graph sparsification (Spielman & Pig
P-g
M
Srivastava 2011 if each edge weight is computed by Wp.g = - Prior research proves that for
every undirected graph a spectrally sparsified graph with O(N log N) edges can be computed by sampling each edge according its effective resistance (Spielman & Srivastava 2011); the trary, GRASPEL can be regarded as a spectral graph densification procedure that aims to include O log N) spectrally-critical edges with large spectral embedding distortions. Therefore, the graphical Lasso problem with Laplacian-lik, precision matrix can be efficiently solved via a spec- tral graph densification procedure that iteratively includes the most spectrally critticll edges into the latest graph that the distortion in the graph spectral emhedding space be drastically mitigated.
4

---

Under review as a conference paper at ICLR 2021
Step 2. Spectral graph embedding
Node embedding vectors
New edge
Data vectors
Gk
If |
P
Step 3A, Check embedding distortion
GO
P.
Step 1. Initial graph construction
Step 3B. Add
Step 4. Check embedding stability
spectrally-critical edge
Figure 1: The overview of the proposed GRASPEL framework.
Consequently, unique feature of GRASPEL is that the spectral embedding (effective-resistance) distances the learned graph will encode the l2 distances between the original data points. which key to many manifold learning and dimensionality reduction problems (Belkin & Niyogi 2003 Carey 2017). We note that the proposed GRASPEL algorithm shares some similar features as the general stagewise algorithm (Tibshirani 2015): as the step size goes to zero the sequence of for- ward stagewise estimates will exactly coincide with the lasso path. Consequently, GRASPEL will produce an approximate solution the original graphical Lasso problem when using rather small step size (e.g. adding only one edge with a small edge weight in each GRASPEL iteration).
Convergence analysis. The global optimal solution can be obtained when 14) becomes zero or there exists no edge with 7 of GRASPEL iterations can be determined based on graph spectral stability: when the maximum embedding distortion 1) becomes small enough (e.g. 1/maz S tol), equivalently when graph spectra become sufficiently stable, the GRASPEL iterations can be terminated.
for inclusion to the latest graph. Consequently. the convergence
Complexity analysis. Comparing with the state-of of-thh--rttt-ttto graph construction methods (Dngg et al. et Dong 2019 Kalofolias & Perraudin 2019) which
2 Kalofolias (0166] [Egilmez et
Dong
require at least O(N2 time in each iteration, (14) allows each GRASPEL iteration to Identify the
most spectrally-criteee
edges in O log N) time: given the subspace projection matrix U, in 12
that can be computed in nearly-linear time leveraging our latest Laplacian eigensolver (Zhao et al,
2018), the spectral embedding distortion (17p.) of each candidate edge (p, q) can be estimated in constant time; then (T] allows identifying the most spectrally-critin
edges almost nearly-linear time by limiting the search within only the candidate edge connections between small number of the top and bottom sorted nodes according to the 1-D spectral embedding using the Fiedler vector.
DETAILED STEPS GRASPEL
Overview, GRASPEL will iteratively identify and add the most spectrally critiial edges into the latest graph that the spectral embedding distortion be greatly mitigated, until such edges can be found (as illustrated Figure] Tp. The detailed GRASPEL algorithm flow for graph learning has been described in Algorithm in the Appendix, and summarized into the following key steps:
Step 1: Initial graph construction. GRASPEL first constructs an ultra-sparse kNN graph to approximate the local manifold of given data set. An extra step of spectral sparsification (Spiclman & Srivastava 2011 Feng 2019) may be adopted to further simplify the graph.
Step 2: Spectral graph embedding. GRASPEL extracts now-dimensional vector repre-
sentations (e.g. v, for node S in Figure [ using scalable spectral graph embedding.
Step 3: Spectrally-critical edge identification. GRASPEL identifies the most spectrally- critical edges with the greatest distortions, and include them into the latest graph.
5

---

Under review as a conference paper at ICLR 2021
Step 4: Spectral stability checking. The GRASPEL iterations are terminated when the embedding distortions become sufficiently small, or equivalently when the graph spectra (e.g. the first few Laplacian eigenvalues and cigenvectors) become adequately stable.
Step I: Initial graph construction. As aforementioned, (approximate) kNN graphs can be con- structed as the initial graph, since they can be obtained efficiently (Muja & Lowe 2009), while being able to approximate the local data proximity (Roweis & Saul 2000). However, the optimal & value (the number of nearest neighbors) is usually problem dependent and can be very difficult to find. In this work, GRASPEL will start with creating an (approximate) kNN graph using a relatively small k value (e.g. h = 2) ), which will suffice for approximating the local structure of the manifold, and strive to iteratively improve the approximation of the global manifold structure by adding small portion of spectrally-criticn
edges through solving the proposed convex problem in 4 In addi-
tion, spectral sparsification (Spielman & Srivastava simplify the initial KNN graph (Wang & 2e11)
2011 Feng 2019) can be leveraged to further
Step 2:
Spectral graph embedding. Spectral graph embedding directly leverages the first few nontrivial cigenvectors for projecting nodes onto -iw-dimeniional space (Belkin & Niyogi [003))
The eigenvalue decomposition of Laplacian matrix is usually the computational bottleneck spec-
tral graph embedding. especially for large graphs (Shi & Malik 2000, Von Luxburg 2007 Chen
ct al., 2011). To achieve good scalability, nearly-lineal time Laplacian solvers (Koutis et 2010)
or multilevel Laplacian solvers (Zhao et al., 2018) can be exploited for much faster eigenvector (eigenvalue) computations without loss of accuracy.
Step 3: Spectrally-critical edge identification. Once Laplacian eigenvectors are available for the latest graph, through the following phases GRASPEL will identify spectrally-critical edges by looking each candidate edge's embedding distortion defined in (3)), while the following theorem can be derived for quantifying each candidate edge's Sage's impact on the first few Laplacian eigenvalues.
Theorem
The perturbation of the i-th Laplacian eigenvalue de and the total relutive spectral perturbation of the first r Laplacian eigenvalues due to the inclusion a candidate edge (p. q) can be estimated by 8λ, = Supp ep.))² and A.
Σ Abs Ae = =dup.all I1,, ep.al12 x fp... respectively.
Proof: See the Appendix.
Phase A: candidate edge identification with Fiedler vectors. Our approach for identifying spectrally- critical edges starts with sorting nodes according to the Fiedler vector that can be computed in nearly-linear time leveraging fast Laplacian solvers (Koutis et al. 2010 Spielman & Teng 2014). This scheme is equivalent to including only the first nontrivial Laplacian eigenvector into U, (r 2) in (IT) for spectral graph embedding. Subsequently, we can search candidate edge connections between the top and bottom few nodes in the I ID sorted node array. According to (14), only a small portion of node pairs with large embedding distances needs to be examined candidate edges.
Phase B: embedding distortion estimation with multiple eigenvectors. With the first r Laplacian eigenvectors computed in the previous spectral embedding step, each node of the latest graph will be associated with an r-dimensional embedding vector, which allows the spectral embedding distor- tion of each candidate sneectrally-critica edge to be quickly estimated. As discussed in Section 3.2, the spectral embedding distances computed with the first eigenvectors can well approximate the effective-resistance distance thus the gradient the proposed optimization task (9) Only the can- didate edges with top embedding distortions will be added into the latest graph. Since the distances approximated using the first few eigenvectors will be the lower bounds of the effective-resistance distances, the lower bound of embedding distortions due to such an approximation can be estimated.
Step 4: Spectral stability checking. In this work, we propose to evaluate the edge embedding distortions with (3) for checking the spectral stability of the learned graph. there exists no ad- ditional edge that has an embedding distortion greater than a given tolerance level (lol), GRASPEL iterations can be terminated. It should be noted that choosing different tolerance levels will result in graphs with different densities. For example, choosing a smaller distortion tolerance will require more edges to be included so that the resultant spectral embedding distances on the learned graph can more precisely encode the distances between the original data points. In practice, even for very large data set GRASPEL converges very quickly when starting with an initial 2NN or uNN (ultra-sparse nearest-neighhor) graph (Wang& & Feng 2017).

---

Under review as a conference paper at ICLR 2021
5 EXPERIMENTS
In this section, extensive experiments have been conducted to evaluate the performance of GRASPEL for variety of public domain data sets (see the Appendix A.. for detailed setting and evaluation metrics). In Sections [A66T and A.662 (see Appendix), we report additional experimen- tal results for spectral clustering and dimensionality reduction (ISNE) applications leveraging the proposed GRASPEL approach. Since this work primarily focuses learning graphs from high- dimensional data points, the proposed method can be orthogonal to existing research related to deep learning based spectral clustering methods: an autoencoder can be first applied to transform the input data into more optimal features that can subsequently become the input of GRASPEL for learning graphs in spectral clustering. For fair comparisons with other state-of-the-art graph learn- ing methods, we directly use the raw data as input without any additional pre-processing steps. The following experiments are performed using MATLAB R2020b running on a Laptop with 10th Intel(R) Core(TM) i5 CPU and 8GB RAM.
5.1 EXPERIMENT SETUP
When applying Algorithm| D to our data sets for the graph learning tasks shown in this section, we randomly sample candidate edges that connect between the top and bottom 0.05 (c = 0.05) nodes in the ID sorted array according to the Fiedler vector, which allows GRASPEL to quickly identify the most spectrally-critical edges. Note that choosing smaller value will allow more efficient edge sampling for estimating global graph (manifold) structural properties, while choosing a greater e value will require more samples but may lead to better preservation of mid-to-short range graph (manifold) structural properties. When estimating the spectral distortion of each candidate edge we compute the first step. We set the edge sampling ratio to be Ç = 0.001 and add one edge (s 1) to the latest graph in each GRASPEL iteration. Clustering Accuracy (ACC) and Normalized Mutual Information (NMI) values (see definitions in the Appendix) are used for evaluating the spectral clustering results.
to 10 Laplacian eigenvectors for the spectral graph embedding
1E3 in (TT) has been used for computing O in all experiments.
5.2 EXPERIMENT RESULTS
Spectral stability checking. In Figure 2 we
show how the objective function defined in (6)
GRASPEL w/ uNN Graph
would change during the GRASPEL iterations
-300
when starting with (a) a 2NN graph, and (b)
Fge-343
333
an ultra-sparse nearest neighbor (uNN) graph
for the USPS data set. Only the first 50 Lapla-
-400
GRASPEL w/ an 2NN Graph
dellt))
cian eigenvalues are computed for evaluating (). The uNN graph is obtained by spectrally sparsifying 5NN graph using the GRASS al- gorithm (with relative condition number of 30) (Feng 2019). We set the distortion toler- ance to be tol = for both cases and demon- strate the results for the first 30 GRASPEL it- erations. As observed in Figure 2 6 grows rapidly within the first 10 iterations when start- ing with the 2NN graph, which indicates the GRASPEL iterations allows quickly converg- ing graph with relatively stable graph spec- tra (b) achieves a greater objective function value after 30 iterations when comparing with (a), which is mainly due to the much sparser structure of the initial uNN graph. Therefore, for very large data sets a spectral sparsifica- tion procedure can be performed before going through GRASPEL iterations for further im- proving graph sparsity and thus runtime effi-
4- o-e 0---0-0---------------0---------
-150
420
Fgm-482
15
20
25
GRASPEL Iterations (tol=1)
10
log max
8
6
2
6
8
10
12
Iterations
Figure 2: Objective function changes (top figure) and maximum spectral distortions (bottom figure) during GRASPEL iterations.
*GRASPEL starts to experience negligible gradient values after 10 iterations and can not identify additional spectrally-critical candidate edges with large emhedding distortions.
7

---

Under review as a conference paper at ICLR 2021
subset the MNIST dataset including 24,462 handwritten digits from to
[8]
(c)
ta|=100 (10 Iterations)
(e)
tol=20 (29 iterations)
0.02
0.8
0.11
02
Large gap
Large
(b)
(d)
500
-34
-30
124
-38
NMI-0.012
-130
ACC-16.15%
ACC-915%
NM1+-.873
ACC-90.7%
2NN Graph
2NN Graph 10 Edges
2NN Graph 29 Edges
The test set of the Fashion MNIST dataset including 10,000 article images from style 0 to
031
2.60
(c)
tol=100
[e)
tol=10
2.01
0.009
Large gap
Large gap
2.02
(f)
(d)
(b)
-282
-24
-2.4
-18
NM 599
58:
ACC-43.9%
NMI-0.611
NM1-0.62 626
0.22
ACC-55.2%
ACC-55.5%
405
2NN Graph
2NN Graph Edges
2NN Graph + 66 Edges
Figure 3: The first few Laplacian eigenvalues and 3D spectral drawings of the 2NN graph in figures (a) and (b), and the GRASPEL- graphs in figures (c) (f).
ciency in graph learning tasks. However, it may be challenging to choose a good Ie for creating that
the initial kNN graph will be further spectrally sparsified: choosing too large & will result in very dense initial graph even after spectral sparsification, whereas choosing a small & may lead to slower convergence.
Spectral embedding distortion. We also show how the proposed scheme for spectral stability checking can be applied based on the embedding distortion metric defined in (I3). GRASPEL iterations can be terminated when there exists no candidate edge that has spectral embedding distortion greater than a given tolerance level (e.g. 7) 2 tol). As shown in Figure 2 for the USPS data set GRASPEL adds one additional edge in each iteration and requires only I1 iterations to effectively mitigate the maximum embedding distortion by over 2, 600x.
Spectral embedding results. In Figure we show the first few Laplacian eigenvalues and 3D spectral drawings of the graphs learned with different distortion tolerance levels for a subset (in- cluding all the 24, 462 handwritten digits from 0 to 6) of the MNIST data set (see Appendix for details) and the test set of the Fashion MNIST data set including 10, 000 article images from style 0 to (Xiao et al. 2017). When creating the 3D spectral drawing layouts, each entry of the first
8

---

Under review as a conference paper at ICLR 2021
0.24
NMI=0.834
0.4
NMI--.882
NMI-0 894
0.28
ACC-75.3%
ACC-72.4%
ACC-92.5%
02
01R
3.2
-642
-0.00
Resistance Correlation: R 0.063
Resistance Correlation: 0.093
Resistance Correlation: R 0.352
12
10
D.1
0.28
0.2
3NN Graph
4NN Graph
2NN Graph 100 Edges
Figure 4: The Pearson correlation coefficients computed with the effective-resistance distances ob- tained on the 3NN graph (left), 4NN graph (middle) and the graph learned by GRASPEL (right).
three nontrivial Laplacian eigenvectors (u2 Ug U.) corresponds to the x, y and & coordinates of each node (data point), respectively. The ground-truth label of each data point is shown using its corresponding color. The edges have been omitted in the graph layouts to better reveal the structure of the data points (manifolds). Different embedding distortion tolerance levels are considered in our experiment. Starting from initial 2NN graph, by adding has dramatically mitigated the spectral embedding distortions from by adding only few additional edges into the initial 2NN graph. By examining the first few Laplacian eigenvalues, we notice that the initial 2NN graph of the MNIST data set has nine connected components (that equals the number of zero cigenvalues), while with 10 extra edges added via GRASPEL iterations, a well-connected graph can be formed for approximately preserving the structure of the original data set. For the MNIST data set we also observe relatively large gaps between the 6th and 7th eigenvalues, indicat- ing that the intrinsic dimensionality of the GRASPEL-learned graphs (manifolds) is approximately five, whereas for the Fashion MNIST data set the intrinsic dimensionality is approximately six.
edge in each iteration GRASPEL
Resistance-distance correlation. In the last, for the full USPS data set we evaluate graph learning quality by checking if the effective- resistancee distances on the learned graph will properly encode the lg distances between the original data points. To this end, we randomly pick up 1, 000 node pairs and compute their effective-resistance distances. To avoid picking up nearby nodes, we first sort nodes according the Fiedler vector and then choose the node pairs (Nopp. Nboor) from the node sets Nzop and Noot formed with the top and bottom 5%|V V| nodes, respectively. Then, these resistance distances
compared with the corresponding & distances between the original data points by checking the Pearson correlation coefficient. In Figure 4 we observe that the effective- resistance distances on the graph learned via 100 GRASPEL iterations (27 seconds) have the highest correlation (R 0.352) with the lg distance between the original data points, while the 3NN or 4NN graphs with much greater edge density only achieve R = 0.063 and It = 0.093, respectively. It is also observed that increasing & from 3 to 4 for constructing the kNN graph can improve the resistance correlation but not necessarily the spectral clustering quality (c.g. NMI and ACC metrics), whereas the graph learned by GRASPEL achieves the best results in all aspects.
6 CONCLUSION
In this work, we present a highly-scalable spectral graph densification approach (GRASPEL) for graph learning from data. By limiting the precision matrix to be a graph-Laplacian-like matrix in graphical Lasso, GRASPEL can always efficiently learns sparse undirected graphs from high- dimensional input data by identifying and including the most spectrally-criticn graph. A unique feature of the graphs learned by GRASPEL. is that the reisstive distances will encode the similarities of the original data points. Compared with offte-of-the--e art graph learning approaches, GRASPEL is more scalable and always leads to substantially improved computing efficiency and solution quality for spectral clustering (SC) and dimensionality reduction tasks.
edges into the latest

---

Under review as a conference paper at ICLR 2021
REFERENCES
Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation, 15(6):1373-139, 2003.
CJ Carey. Graph Construction for Manifold Discovery. PhD thesis, University of Massachusetts Amherst, 2017.
Yen Chen, Yangqiu Song. Hongjie Bai, Chih-Jen Lin, and Edward Y Chang. Parallel spectral clustering in distributed systems. IEEE transactions on pattern analysis and machine intelligence, 33(3):668- 586, 2011.
Xiaowen Dong, Dorina Thanou, Pascal Frossard, and Pierre Vandergheynst. Learning laplacian matrix in smooth graph signal representations. IEEK Transactions Signal Processing, 64(23): 6160-6173, 2016.
Xiaowen Dong. Dorina Thanou, Michael Rabbat, and Pascal Frossard. Learning graphs from data: A signal representation perspective. IEEE Signal Processing Magazine, 36(3):44 63, 2019.
Hilmi E Egilmez, Eduardo Pavez, and Antonio Ortega, Graph learning from data under laplacian and structural constraints. IEEE Journal of Selected Topics in Signal Processing, 11(6):825-841, 2017.
Zhuo Feng. Spectral graph sparsification in nearly-linear time leveraging efficient spectral perturba- tion analysis. In Proceedings of the 53rd Annual Design Automation Conference, pp. 57. ACM, 2016.
Zhuo Feng. Similarity-aware spectral sparsification by edge filtering. ACM/ESDA/IEEE Design Automation Conference (DAC), pp. 1-6. IEEE, 2018.
In 2018 55th
Zhuo Feng. Grass: Graph spectral sparsification leveraging scalable spectral perturbation analysis.
arXiv preprint 19X1::111..04533 2019.
Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432-441, 2008.
Shivam Garg, Vatsal Sharan, Brian Zhang, and Gregory Valiant. A spectral view of adversarially robust features. In Advances Neural Information Processing Systems, pp. 10138-10148, 2018.
Xiaojie Guo. Robust subspace segmentation by simultaneously learning data representations and their affinity matrix. In IJCAI, pp. 3547-3553, 2015.
Tony Jebara, Jun Wang, and Shih-Fu Chang. Graph construction and b-matching for semi-supervised learning. In Proceedings of the 26th annual international conference on machine learning, pp- 441-448. ACM, 2009.
Vassilis Kalofolias. How to learn graph from smooth signals. In Artificial Intelligence and Statis- tics, pp. 920-929, 2016.
Vassilis Kalofolias and Nathanaël Perraudin. Large scale graph learning from smooth signals. In- lernational Conference on Learning Representations (ICLR 2019), 2019.
I. Koutis, G. Miller, and R. Peng. Approaching Optimality for Solving SDD Linear Systems. In Proc. IEEE FOCS, pp. 235-244, 2010.
Sandeep Kumar, Jiaxi Ying, Jose Vinicius de Miranda Cardoso, and Daniel Palomar. Structured graph learning via laplacian spectral constraints. In Advances in Neural Information Processing Systems, pp. 1651111663,, 2019.
Brenden Lake and Joshua Tenenbaum. Discovering structure by learning sparse graphs. 2010.
George C Linderman and Stefan Steinerberger. Clustering with I-sne, provably. arXiv e-print, 170.00222.. 2017.
Yang Liu, Quanxue Gao, Zhaohua Yang, and Shujian Wang. Learning with adaptive neighbors for image clustering. In IJCAI, pp. 2483-2488 2018.
10

---

Under review as a conference paper at ICLR 2021
Laurens van der Maaten and Geoffrey Hinton. Visualizing Data using t-SNE. Journal of machine learning research, 9NNv))2359--6655, 2008.
Markus Maier, Ulrike Luxburg, and Matthias Hein. Influence of graph construction on graph- based clustering measures, In Advances in neural information processing proccessn systems, 1025-1032, 2009.
Marius Muja and David G Lowe. Fast approximate nearest neighbors with automatic algorithm configuration. VISAPP (1), 2(331-340):2 2009.
Christos H Papadimitrou and Kenneth Steiglitz. Combinatorial optimization: algorithms and com- plexity. 1982.
Vittal Premachandran and Ramakrishna Kakarala. Consensus of k-nns for robust neighborhood selection on graph-based manifolds. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 1594-1601, 2013.
Michael G Rabbat. Inferring sparse graphs from smooth signals with theoretical guarantees. In 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp- 6537. IEEE, 2017.
Sam T Roweis and Lawrence Saul. Nonlinear dimensionality reduction by locally linear embed- ding. science, 290(5500):2333-2255 2326, 2000.
Jianbo Shi and Jitendra Malik. Normalized cuts and image segmentation. IEEE Transactions on pattern analysis and machine intelligence, 22(8):888-906 905, 2000.
Martin Slawski and Matthias Hein. Estimation of positive definite m-matrices and structure learning for attractive gaussian markov random fields. Linear Algebra and its Applications, 473:145-111 179, 2015.
Daniel Spielman and Nikhil Srivastava. Graph Sparsification by Effective Resistances. SIAM Jour- nal on Computing, 40(6):1913- 2011.
Daniel Spielman and Shang-Hua Teng. Nearly linear time algorithms for preconditioning and solving symmetric. diagonally dominant linear systems. SIAM Journal on Matrix Analysis and Applications, 35(3):835-885, 2014.
Alexander Strehl and Joydeep Ghosh. Cluster ensembles knowledge reuse framework for com- bining multiple partitions. Journal of machine learning research, 3(Dec):583-61) 617, 2002.
Joshua Tenenbaum, Vin De Silva, and John C Langford. A global geometric framework for nonlinear dimensionality reduction. science, 29005500))2311-233338 2000.
Ryan J Tibshirani. A general framework for fast stagewise algorithms. The Journal of Machine Learning Research, 16(1):254332555) 2015.
Laurens Van Der Maaten. Accelerating t-SNE Using Tree-based Algorithms. The Journal of Ma- chine Learning Research, 15(1):3221-3245 2014.
Ulrike Von Luxburg. A tutorial on spectral clustering. Statistics and computing, 17(4)399-4166 2007.
Yongyu Wang and Zhuo Feng. Towards scalable spectral clustering via spectrum-preeseringg sparsi- fication. arXiv e-print, 1XXi::7110045000 2017.
Han Xiao, Kashif Rasul, and Roland Vollgraf. Pashion-mnist. a novel image dataset for benchmark- ing machine learning algorithms, 2017.
Zhiqiang Zhao, Yongyu Wang, and Zhuo Feng. Nearly-linear time spectral graph reduction for scalable graph partitioning and data visualization. arXiv e-print,
2018.
11

---

Under review as a conference paper at ICLR 2021
A APPENDIX
A.1 PROOF OF THEOREM 1
Let L RNNN denote the Laplacian matrix of an undirected graph, and us denote the i-th eigen- vector of I. corresponding to the i-th eigenvalue X that satisfies:
Lus X,
(15)
then we have the following eigenvalue perturbation analysis:
(I. 5L. (1+ 8mz) (A: dA:)) (u: + dus)
(16)
where perturbation 8T. = dw pup..gpp..ppp that implies new edge connection is applied I, result- ing in perturbed eigenvalues and eigenvectors X; + 8λ; and u, 5u, for i 1,,., N, respectively.
Keeping only the first-order terms leads to:
Lou, + 5Lu = Aqou, + 87,U4
(17)
Write du, in terms of the original eigenvectors u, for for I,.., N:
N
Su, a;u
(18)
Substituting (18) into L7 leads to:
N a + 8L44 = A Σ
L
a;4 + 67
(19)
2m1
Multiplying of both sides of 19 results in:
N a,u, + u 8L4 = Aru
u L
Σ
a:44 + 8X,2 Us.
(20)
1=1
Since un for for i l,... N are unit-length, mutuaaly--rthogoonal eigenvectors, we have:
u, L Σ ag4 a.u! Lug,
Ahay Σ a,u; aquu Aru.
(21)
Substituting 115 into 21
we have:
au Lus andt Aqui.
(22)
According to
we have:
0g44 Auu Σ
(23)
1=1
1=1
Substituting
23 into 20) leads to:
u 8LL; 8X.u7 u; 67
(24)
Then the eigenvalue perturbation 8A, due to 8 is given by:
8), Supp (u. ena)
(25)
12

---

Under review as a conference paper at ICLR 2021
The total relative perturbation of the first eigenvalues due the inclusion of edge (p, () becomes:
5\ 8up.allU, Cp.a 12 8wp.. 2p
.emb
(26)
where 2 denotes the spectral embedding distance between and when U, is formed with the first r 1 weighted nontrivial eigenvectors as defined in (TT). If each edge weight computed by
Swp.g with data
|| cp.g [. as long as we can find a spectrally-critical edge with large
Pp-
P.
1 M up.g. or equivalently spectral embedding distortion, including this edge into the current graph will significantly perturb the first Laplacian eigenvalues.
A. adata
A.2 ALGORITIM FLOW
Algorithm The GRASPEL Algorithm Flow
Input: A data matrix (X [x]....IM] E RV with N data points in M M-dimensional, embedding distortion tolerance (1 I tol), window size for edge sampling (0 < S 50%), edge sampling ratio (0 < s 1), and the number of edges to be selected in each iteration (0 < 3). Output: The spectrally-learned graph G.
1: Construct an initial 2NN graph G using approximate kNN algorithms.
2: while Mmaz 2 tol do
3: Embed the latest graph G using its Fiedler vector and sort the nodes into a ID array Inode'
4: Obtain node set Noop ( (Noot) by including only the top (bottom) [eN]] nodes in Inode:
5: Sample each of the [8/[] edges by randomly choosing node from Noop and another node
from Nooti
6: Form edge set Eset using edges with large distortions (1) tol) and set the largest edge embedding distortion as Mjaar.
7:
if Eset 2 then
8: Add the top edges with largest from Ezet into G; else
9:
Add all the edges in Eset into G; end if
11:
end while
Return the learned graph G.
Algorithm 2 Spectral Clustering Algorithm
Input: A graph (V, w.) and the number of clusters Output: Clusters C,..C,,.
1: Compute the adjacency matrix A, and diagonal matrix D;
2: Obtain the unnormalized Laplacian matrix L=D-A;
3: Compute the eigenvectors uy that correspond the bottom nonzero eigenvalues of L;
4: Construct U, E I2N xr with r eigenvectors of L stored as columns;
5: Perform k-means algorithm to partition the rows of U, into r clusters and return the result.
A.3 DATA SETS DESCRIPTION
COIL20: the data set contains 4.4 gray-scale images 20 objects, and each object on turntable has 72 normalized gray-scale images taken from different degrees. The image size is 32x 32 pixels.
PenDigits: the data set consists of 7.494 images of handwritten digits from 44 writers, using the sampled coordination information. Each digit is represented by 16 attributes.
USPS: the data set includes 9, 298 scanned hand-written digits from 0 on the envelops from U.S. Postal Service with 256 attributes.
MNIST: the data set consists of 70,000 images of handwritten digits. Each image has 28-by-28 pixels in size. This database can be found at website (http://yann.leeen wtt..//ccu..cc//e////iiiise
13

---

Under review as a conference paper at ICLR 2021
A.4 COMPARED ALGORITTM
Standard KNN: the most widely used affinity graph construction method. Each node is connected to its A nearest neighbors.
Consensus of KNN (cons-kNN) (Premachandran & Kakarala 2013): adopts the fute-of-the-an neighborhood selection methods to construct the affinity graphs. It selects strong neighborhoods to improve the robustness of the graph by using the consensus information from different neighbor- hoods in a given kNN graph.
LSGL Kalofolias & Perraudin 2019:: a method to automatically select the parameters of the model introduced in (Kalofolias) 2016) given desired graph sparsity level. The default settings have been used in our experiments.
A.5 EVALUATION METRIC
(1) The ACC metric measures the agreement between the clustering results generated by clustering algorithms and the ground-truth labels. A higher value of ACC indicates better clustering quality. The ACC can be computed by:
Σ Sys, map(c;))
ACC
(27)
where n is the number of samples in the data set, y; is the ground-truth label provided by the data sets, and is clustering result obtained from the algorithm. 5(x, y) is delta function defined as: ((x, y)=1 for x y. and 6(2, y/)0., otherwise. map(o) permutation function that maps each cluster index e to ground truth label, which can be realized using the Hungarian algorithm (Papadimitrou & Steiglitz 1982).
(2) The NMI metric the range of [0, 1J, while a higher NMI value indicates a better matching between the algorithm generated result and ground truth result. For two random variables and Q, normalized mutual information is defined as (Strehl & Ghosh] 2002
NMI I(P,Q)
(28)
[H(P)H(Q)
where T(P, Q) denotes the mutual information between and Q. while H(P) and H(Q) are en- tropies of P and Q. In practice, the NMI metric can be calculated follows (Strehl & Ghosh 2002
Σ 1=1==1
ni.j log( m.mis.) ny.ny
NMI
(29)
k
n: no Te )
E ny nslog %
where n is the number of data points in the data set, k is the number of clusters, n, is the number of data points in cluster C; according to the clustering result generated by algorithm, n; is the number of data points in class C, according to the ground truth labels provided by the data set, and n.. is the number of data points in cluster C; according to the clustering result as well as in class C according to the ground truth labels.
A.6 ADDITIONAL EXPERIMENTAL RESULTS
A. GRAPH LEARNING FOR SPECTRAL CLUSTERING (SC)
The classical spectral clustering (SC) algorithm (see Algorithm 2 in the Appendix) first constructs kNN graph where each edge weight encodes similarities between different data points (entities); then SC calculates the eigenvectors of the graph Laplacian matrix and embeds data points into low- dimensional space (Belkin & Niyogi 2003); in the last, k-means algorithms are used to partition the data points into multiple clusters. The performance of SC strongly depends on the quality of the underlying graph (Guol 2015). In this section. we apply GRASPEL for graph construction, and show the learned graphs can result in drastically improved efficiency and accuracy in SC tasks.
14

---

Under review as a conference paper at ICLR 2021
Since the SC algorithm has intrinsic randommess, the clustering result in each run is different. So the reported numbers in our results have been averaged over 20 runs.
We show comprehensive results of SC using four graph learning (construction) methods in Table [] includes the total time for eigendecomposition of the Laplacian matrix and k-means clustering. The graph density (EE)/V) results are also shown in Table In Table the runtime of the consensus method includes the time for consensus information calculation and edge pruning, while the runtime of GRASPEI. includes the total time for spectral graph densification. As observed, GRASPEL can consistently achieve the state-of-the-eee results in SC: the GRASPEL learned graphs have led to significantly better SC accuracy (ACC), much lower graph densities and much shorter graph learning time.
and Table 2 The runtime reported in Table 2
Table 1: ACC and NMI results
ACC(%) NMI
Data Set
Standard K-NN
Consensus
LSGL
GRASPEL
COIL20
75.72/0.86
81 81.60/0.90 90
85/49/0
86.46/0.94
PinDigits
74.36/0.79
71.08/0.79
74.53/0.77
82.4070.79
USPS
64.3170.79
68.5470.8T
81.50/0.84
91.50/0,89
MNIST
64.20/0.74
74.63/0.78
Table 2: SC Runtime Results
Spectral clustering time (seconds)
Data Set
Standard k-NN Consensus
LSGL GRASPEL
COIL20
0.03
0.03
0.08
0.02
Penl Digits
0.18
0.16
4.42
0.17
USPS
0.72
0.56
7.05
0.28
MNIST
252.59
3.06
Table 3: Graph density results
Graph density (ET/V)
Data Set
Standard K-NN Consensus
LSGL.
GRASPEL.
COIL20
6.12
5.06
11.99
1.39
PenDigits
6.76
6.70
186.52
2.96
USPS
7.30
6.58
29.97
MNIST
7.46
1.72
Table 4: Graph learning (construction) time results
Graph construction time (seconds)
Data Set
Consensus
LSGL
GRASPEL
COIL20
2.43
13.56
0.29
PenDigits
172.51
1085.43
2.04
USPS 574.28
2074.78
3.37
MNIST
208.89
In Figures and] additional SC results have been provided for the USPS and the Pendigit data sets by comparing two embedding distortion tolerance levels (tol = 100 and tol =10). Not surpris- ingly, when starting with initial 2NN graphs few GRASPEL iterations have already dramatically improved SC results: the normalized mutual information (NMI) has been improved from 0.612 to 0.888 for the USPS data set, and from 0.04 to 0.840 for the Pendigit data set; The spectral clustering accuracy (ACC) has been improved from 40.4% 91.5% for the USPS data set, and from 15.2% to 89.4% for the Pendigit data set.
Since the number of zero Laplacian eigenvalues equals to the connected components in the learned graph, we observe that GRASPEL can always identify O(g) spectrally-critical edges added the initial 2NN graph so that its q connected components immediately get stitched into connected graph for both test cases.
15

---

Under review as a conference paper at ICLR 2021
The full USPS dataset including 9,298 handwritten digits from to
0023
(8)
(c) tol=100 (9 iterations)
0.02
(e)
tol=10 (20 lterations)
0.01
0.005
10
15
2015
(b)
(d)
4.02
NMI=0.612
NMI-0.883
NMI--0888
ACC-40.4%
ACC-90.7%
031
-001
2NN Graph
2NN Graph 9 Edges
2NN Graph Edges
Figure 5: The first 20 Laplacian cigenvalues (top) and 3D spectral drawings (bottom) of the 2NN graph in figures (a) and (b), and the GRASPEL-leamedd graphs in figures (c) (f).
subset of the Pendigit dataset including 5,278 handwritten digits (0 to 6)
2000
2000
(a)
[c) tol=100 (17 iterations)
(e) tol=10 (22 Hteraiions)
1500
1200
1000
500
500
30
0
10
20
30
50
0.6
0.00
ACC-7A.2%
0.04
NMIn0 ACC-89.4%
0.4
(b)
ACC-15.2%
0.21
(d)
(f)
4000
062
-22
02
0.02
0.4
-012
2NN Graph
2NN Graph Edges
2NN Graph Edges
Figure 6: The first 50 Laplacian eigenvalues (top) and 3D spectral drawings (bottom) of the 2NN graph in figures (a) and (b). and the GRASPEL-leamed graphs in figures (c) (f).
16

---

Under review as a conference paper at ICLR 2021
USPS data set
standard t-SNE: 99.8s
5X reduction t-SNE: 14.1s
MNIST data set
standard t-SNE: 1921.7s
5X reduction t-SNE: 149.3s
Figure 7: Multilevel t-SNE visualization results.
Discussion. The superior performance of GRASPEL is due the following reasons: 1) In tradi- tional KNN graphs, all the nodes have the same degrees; as a result, the clustering may strongly favor balanced cut, which may lead to improper cuts in high-density regions of the graph. In con- trast, GRASPEL always learns ultra-sparse graphs that only include edges with the largest impact to graph spectral (structural) properties; as a result, the corresponding cuts will always occur in proper regions of the graph, which enables handle even unbalanced data. 2) Recent work (Garg al., 2018) shows the fundamental connections between spectral properties graphs associated with data and the inherent robustness to adversarial examples. Since GRASPEL identifies candidate edges by leveraging spectral graph properties, the learned graph structure will also be robust input noises (perturbations).
A.6.2 GRAPH LEARNING FOR DIMENSIONALITY REDUCTION (DR)
The t-Distributed Stochastic Neighbor Embedding (t-SNE) has become one of the most popular vi- sualization tools for high-dimensional data analytic tasks (Maaten & Hinton 2008 Linderman & Steinerberger] 2017). However, its high computational cost limits its applicability to large scale problems. An substantially improved t-SNE algorithm has been introduced based on tree approxi- mation (Van Der Maaten 2014). However, for large data set the computational cost can still very high.
A multilevel t-SNE algorithm has been proposed in (Zhao al. 2018 coarsening as a pre- pre-processing step applied to the original kNN graph. A much smaller set of representative data points can be then selected from the coarsened graph for t-SNE visualization. In this work, we use GRASPEL to learn sparse graphs that can be further reduced into much smaller ones using spectral graph coarsening. Then more efficient t-SNE visualization can be achieved based on the sampled data points corresponding to the nodes in the coarsened graphs.
leveraging spectral graph
17

---

Under review as a conference paper at ICLR 2021
In experiments, first construct initial graphs by applying a spectral sparsification procedure to the 5NN graphs of both the MNIST and USPS data sets. Then spectral graph coarsening procedure (Zhao et al. (0188) has been applied to create hierarchy of coarse-level graphs. The be obtained by directly using the data points corresponding to the nodes
I-SNE visualization the coarsest graph. Figure shows the visualization and runtime results the standard t-SNE (with tree-based acceleration) (Van Der Maaten, 2014) and the multilevel t-SNE algorithm (Zhao al., 2018) based on graphs learned by GRASPEL. When using a 5X graph reduction ratio, I-SNE can be dramatically accelerated (12.8X and 7X speedups for MNIST and USPS data sets, respectively) without loss of visualization quality.
18


Reviews:
[{"review_id": "dV5cuVxttBa", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary and significance: Learning a graph from data is an important, yet less studied, problem. The proposed algorithm (GRASPEL) is based on a graphical Lasso formulation with the precision matrix restricted to be a graph Laplacian. The algorithm starts with a sparse kNN graph, and recursively adds critical edges (identification of these critical edges based on Lasso and spectral perturbation analysis is the main contribution of the paper). \nThe outcome is a highly scalable that learns a graph in nearly linear time (ignoring log factors and number of recursions). The scalability of the algorithm makes the contributions significant.\n\nOriginality: The basic formulation and idea of selecting spectrally critical edges seem original and interesting, although the reviewer is not an expert in related methods. The authors should note that the there are graph learning methods based on solution of graphical lasso [Pavez, Ortega, ICASSP 2016; Kumar et al, Neurips 2019]. Beyond this step, the authors employ several existing techniques to make GRASPEL scalable (although this part is not highly novel, the overall method is original).\n\nQuality and clarity: The theory in the paper is technically sound. The only exception (this is also an issue about clarity) is the assumption that $U_N^T e_{pq} \\approx U_r^T e_{pq}$. In general, this is not valid and hence it should be clarified when this assumption is reasonable. The experimental section is well executed. \nThe paper presented is mostly good. It would help to include Algorithm 1 and Table 1 in the main paper. \n\nTypo: On page 6 (Phase A), perhaps the authors meant (r=2) instead of (r=1). Also reference of Carey is currently in all caps\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "YKyB6Xkffbf", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper studies ways of adding edges to graphs to improve the result of spectral embedding / clustering. It refines existing embeddings using by measuring edges' effect on Laplacian eigenvalues, and adjusting such edges to reduce the distortions. The performance of the algorithm is justified using developments of worst-case efficient algorithms for Laplacian matrices, and experimentally, the algorithm converges quickly when starting with nearest neighbor graphs, and leads to significant increases in accuracy.\n\nStrengths:\n+ the paper puts together a lot of different ideas, many of which have solid theoretical foundations.\n+ full experimental evaluation on a moderate sized data set that demonstrates both good results and good performance.\n\nWeaknesses\n- the presentation could use significant improvement\n- the ideas are a bit disconnected, at least when one try to follow the ideas mathematically.\n\nI find the approach taken by this paper quite interesting: node embeddings are now widely used to preprocess graphs into vector data more friendly to learning pipelines. Many of these methods add additional edges imperatively / via local methods, without taking the overall data set into account. Iterating this process based on the overall embedding is an interesting, but algorithmically more intensive approach. To carry this out, the authors combined a variety of interesting tools, as well as some high performance packages that they developed. I feel the novelty of this combination, as well as the gains obtained, should make this result of broad interest to those working on spectral clustering/embeddings. However, many of the other reviews point out gaps in both the ideas and presentation, and I'm inclined to agree with them that this paper can benefit from a thorough revision before appearing at a major conference.", "extra_scores": {}, "extra_text": {}}, {"review_id": "QTPuBCK0QVH", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes a graph learning method for spectral embedding and associated problems such as clustering and dimension reduction. What differentiates the method from much the existing literature is that it focuses approximating an optimal densification of a very sparse initial graph rather than on sparsification of an initial graph, as is more common. The method is based on iteratively identifying edges to add to the graph so as to best improve the corresponding spectral embedding, so called \"spectrally critical\" edges. The authors motivate spectral criticality in relation to the partial derivatives of an objective function inspired by the log-likelihood of a Gaussian graphical model. In particular, those with the highest partial derivatives will tend to be those which, through their addition to the graph, lead to the greatest increase in this objective. The authors go on to discuss a close connection between spectral criticality and distance distortion when comparing the spectral embedding and the original input space. Since the initial graph is very sparse it can be efficiently determined, and the relatively small number of additional edges which need to be added by the proposed method to obtain a high quality embedding means that the entire procedure can be implemented efficiently.\n\nThe paper is well written for the most part, and the method is intuitive and persuasive. The connection between the partial derivatives and the distance distortions is very pleasing and provides a commonsense interpretation of the operations of the method. In addition the empirical performance of the method on some important benchmarks seems to be good.\n\nMy main concern about the paper is how the method is presented. The connection to the graphical LASSO seems at first to be pivotal, but then the LASSO component is dropped (\\beta set to 0, or effectively ignored in the actual algorithm). Without this connection, then, it isn't clear what is the motivation for the objective being optimised? It seems far more natural to me to motivate the method from the point of view of the distance distortions being used to select edges, and to treat the connection to the log-likelihood objective as an interesting theoretical point. Another connection beyond the graphical LASSO is given in relation to Bregman distances, but there is far too little discussion given for this connection to motivate the use of the objective. Furthermore, although I genuinely do appreciate the motivation of distance preservation, when we think about manifold learning it is only really the local distance structure which is of great importance, and so some of the justification is lost. It is intuitively the case that since the distance distortions are determined as ratios that this will implicitly pick up on deviations in smaller distances over larger ones in any case, so this is practically not a problem, but the authors do not mention this fact and so the direct justification for why this approach works for manifold learning is somewhat lacking. Finally, while the empirical performance shows some promise of the method, it unfortunately leaves a few important questions unanswered. Notably, if we consider Figure 2 we see that the proposed approach rapidly improves on the 2NN graph, and \"converges\" after relatively few iterations. Using the authors' interpretation of \"convergence\" in this case, however, it looks like their method does not improve appreciably on the uNN graph. Furthermore the initial uNN graph is far superior in objective value to the modified 2NN graph. I understand that the uNN graph used for initialisation may be far denser than the final GRASPEL graph which started from the 2NN graph, but this isn't discussed, nor is it clear that applying GRASPEL to the 2NN graph is computationally superior to just starting with the uNN solution and not using GRASPEL at all.\n\nThe paper and the proposed method clearly have some strong points, however in its current form I am concerned it leaves too much not adequately clear for the reader.\n\nIn addition to the points above, below find a few minor comments/questions/typos.\n- How does the discussion of smooth signals on graphs connect to the proposed approach? Algebraic connectivity seems like a more natural connection to the way the method is posed.\n- what are \"attractive\" Gaussian Markov Random Fields?\n- typo pg 2 \"a undirected\" -> \"an undirected\"\n- I'm not used to seeing matrices divided by scalars. While it isn't ambiguous, I'd recommend \\frac{1}{\\sigma^2}I as opposed to \\frac{I}{\\sigma^2}\n- I am confused by the dimensionality. You mention X is M observations on N data entities, which I interpret as sample size  = N and dimension = M. But this doesn't match the description of the sample covariance matrix formulation in the footnote on pg 3. It is also discussed that one makes M i.i.d. observations for this connection to a covariance matrix.\n- Why are the embedding distortions proportional to M and not M/r? The latter makes better intuitive sense since we would expect distances to scale roughly with the square root of dimensionality.\n- typo (?) pg 4: \"there exists no edge with \\eta > 1 can be found...\"\n- In Phase A of step 3 why only search the extreme points in the first non-trivial eigenvector? It seems critical points can show up in subsequent eigenvectors as well\n- Apologies if I missed this point, but is it the case that the columns of the data matrix are scaled to have unit norm? If not, wouldn't it be that the distortion distances could be arbitrarily large, since the eigenvectors are normalised?\n- In the experimental setup, is it the case that you sample 1/1000 edges which connect two points which lie in the extrema of the Fiedler vector? this is an extremely small number. Not a problem, but if the Fiedler vector is clearly indicative of criticality, why not simply use fewer than 5% each end and look at all pairwise distances among these potentially critical points.\n- The discussion on page 7 relating to Figure 2: \"As observed in Figure 2 (a) achieves a much greater objective function value after 30 iterations when comparing with (b)\" Doesn't the figure show the opposite of this? The light blue line is (a) and the red line is (b)\n\n#########################################################\n\nFinal Recommendation:\nI have considered the authors' responses to my comments, as well as the assessments given by other reviewers. I still feel as though, while the method looks as though it may offer a potentially useful practical alternative to other graph learning methods, in its current form I do not think it is presented in a manner which warrants acceptance at a prestigious conference such as ICLR. If the decision overall is that the paper is not to be accepted, then I wish the authors well with their work and hope they take into consideration the comments of the reviewers as I do believe the work has potential.", "extra_scores": {}, "extra_text": {}}, {"review_id": "KLzb5uiNYQf", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Update: In the revised version, the authors have addressed some of my technical concerns, therefore I slightly increased my score.\n\n\nSummary:\n\nThe paper proposes a method to learn a sparse graph from data by optimizing an objective similar to the graphical Lasso over the set of valid Laplacian matrices. The method starts with an initial sparse kNN graph which has been further sparsified using spectral sparsification algorithms. It then applies an approach referred to as \"spectral graph densification\", where spectrally-critical edges are identified which have large embedding distortions and thus significantly impact the objective function. The approach is terminated based on spectral stability checking, where the iterations are terminated when the embedding distortions become sufficiently small.\n\n\nStrengths:\n\n- The paper proposes a new approach for learning sparse graphs from data. \n- Compared to prior spectral graph sparsification algorithms which remove edges from a given graph, the proposed approach starts with a sparse graph and adds edges into the graph iteratively. It is shown in the experiments that this approach leads to a fast improvement in the objective function values in a few iterations.\n- The proposed approach has running time in $O(N \\log N)$ instead of $O(N^2)$ for each iteration and is thus faster than previous graph learning methods.\n\n\nConcerns:\n\n- The discussion in the beginning of Section 3 seems to suggest that the formulation of the objective (4) is an individual contribution of the paper. However, the convex optimization problem in (4) has been previously proposed in [Lake & Tenenbaum, 2010] (see also the discussion in [Dong et al, 2019]). The paper mentions [Dong et al, 2019] but also a reference to [Lake & Tenenbaum, 2010] should be provided.\n- The proposed approach relies on a number of other techniques: it requires as input a kNN graph which is then further sparsified using other spectral sparsification algorithms ([Spielman & Srivastava, 2011, Feng, 2018]) and then a lower-dimensional vector representation is computed by using a nearly-linear-time spectral graph embedding procedure ([Zhao et al, 2018]). The only technical contribution of the paper is then to provide a greedy optimization scheme for the objective in (4) by selecting the spectrally critical edges.\n- It is stated in the introduction that the experimental results show that the graphs learned by the proposed technique can lead to more efficient and accurate spectral clustering (SC) as well as dimensionality reduction. However, the results for spectral clustering and dimensionality reduction are only in the appendix and not discussed in the main part of the paper. In order to make the paper self-contained, at least parts of these results should be moved to the main paper. To make space for these results, one suggestion is to merge the sections 3.3 and 4 into one as they are a bit redundant.\n\n\nMinor Comments:\n\n- After (4) it says that \"It can be shown that the three terms in (4) are corresponding to $\\log \\det (\\Theta)$,  $Tr(\\Theta S)$ and $\\beta||\\Theta||_1$ in (3), respectively\". What is meant by \"correspond\" in this context?\n- In equation (11), what is the index $r$? It has not been defined before.\n- As (14) is directly obtained from (12) by plugging (13) into the left hand side of (12), why is there $ \\approx$ in (14) instead of $=$ ? Also, what is the purpose of the right hand side of (12) (after the $\\geq$)? \n- The claim that the spectral embedding distances on the learned graph will encode the $l_2$ distances between the original data points should be further clarified.\n- Theorem 1 is a bit imprecise, as it makes a statement about the \"spectral criticality\" of a candidate edge, which has not been formally defined before (spectral critical edges are only introduced informally before as \"the ones which can most effectively perturb the graph spectral properties\"). The proof in the appendix then computes the relative spectral perturbation of the first $r$ eigenvalues due to the inclusion of the edge $(p,q)$ but does not precisely relate it to the statement in the Theorem.\n- Is the notion of \"graph spectral stability\" novel? If not, please provide a reference.\n\n\n\nConclusion:\n\nThe idea of using spectral densification to optimize the graph learning objective is nice and elegant. However, I vote for rejection as the performance of the method is reliant on a number of other techniques for the construction of the initial sparse graph and fast computation of the spectral embedding. Moreover, the presentation of the paper could be improved (see comments above).\n\n\n\nReferences:\n\n[Dong et al, 2019] Xiaowen Dong, Dorina Thanou, Michael Rabbat, and Pascal Frossard. Learning graphs from data: A signal representation perspective. IEEE Signal Processing Magazine, 36(3):44–63, 2019.\n\n[Feng, 2018] Zhuo Feng. Similarity-aware spectral sparsification by edge filtering. In 2018 55th ACM/ESDA/IEEE Design Automation Conference (DAC), pp. 1–6. IEEE, 2018.\n\n[Lake & Tenenbaum, 2010] B. Lake and J. Tenenbaum. Discovering structure by learning sparse graph. Proceedings of the Annual Cognitive Science Conference, 2010.\n\n[Spielman & Srivastava, 2011] Daniel Spielman and Nikhil Srivastava. Graph Sparsification by Effective Resistances. SIAM Journal on Computing, 40(6):1913–1926, 2011.\n\n[Zhao et al, 2018] Zhiqiang Zhao, Yongyu Wang, and Zhuo Feng. Nearly-linear time spectral graph reduction for scalable graph partitioning and data visualization. arXiv e-print, arXiv:1812.08942, 2018.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers generally like the paper, in particular the scalability of the proposed approach. The author response and revised version clarified some questions of the reviewers, however, it didn't fully mitigate their concerns.
Decision: Reject