Nicolas Tremblay
email: nicolas.tremblay@grenoble-inp.fr

Andreas Loukas
email: andreas.loukas@epfl.ch

Approximating Spectral Clustering via Sampling: a Review

Spectral clustering refers to a family of well-known unsupervised learning algorithms. Rather than attempting to cluster points in their native domain, one constructs a (usually sparse) similarity graph and computes the principal eigenvectors of its Laplacian. The eigenvectors are then interpreted as transformed points and fed into a k-means clustering algorithm. As a result of this non-linear transformation, it becomes possible to use a simple centroid-based algorithm in order to identify non-convex clusters, something that was otherwise impossible. Unfortunately, what makes spectral clustering so successful is also its Achilles heel: forming a graph and computing its dominant eigenvectors can be computationally prohibitive when dealing with more that a few tens of thousands of points. In this chapter, we review the principal research efforts aiming to reduce this computational cost. We focus on methods that come with a theoretical control on the clustering performance and incorporate some form of sampling in their operation. Such methods abound in the machine learning, numerical linear algebra, and graph signal processing literature and, amongst others, include Nyström-approximation, landmarks, coarsening, coresets, and compressive spectral clustering. We present the approximation guarantees available for each and discuss practical merits and limitations. Surprisingly, despite the breadth of the literature explored, we conclude that there is still a gap between theory and practice: the most scalable methods are only intuitively motivated or loosely controlled, whereas those that come with end-to-end guarantees rely on strong assumptions or enable a limited gain of computation time.

Introduction

Clustering is a cornerstone of our learning process and, thus, of our understanding of the world. Indeed, we can all distinguish between a rose and a tulip precisely because we have learned what these flowers are. Plato would say that we learned the Idea -or Form [120]-of both the rose and the tulip, which then enables us to recognize all instances of such flowers. A machine learner would say that we learned two classes: their most discriminating features (shape, size, number of petals, smell, etc.) as well as their possible intra-class variability.

Mathematically speaking, the first step on the road to classifying objects (such as flowers) is to create an abstract representation of these objects: with each object i we associate a feature vector p i ∈ R d , where the dimension d of the vector corresponds to the number of features one chooses to select for the classification task. The space R d in this context is sometimes called the feature space. The choice of representation will obviously have a strong impact on the subsequent classification performance. Say that in the flower example we choose to represent each flower by only d = 3 features: the average color of each RGB channel (level of red, green and blue) of its petals. This choice is not fail-proof: even though the archetype of the rose is red and the archetype of the tulip is yellow, we know that some varieties of both flowers can have very similar colors and thus a classification solely based on the color will necessarily lead to confusion. In fact, there are many different ways of choosing features: from features based on the expert knowledge of a botanist, to features learned by a deep learning architecture from many instances of labeled images of roses and tulips, via features obtained by hybrid methods more-or-less based on human intelligence (such as the first few components of a Principal Component Analysis of expert-based features).

The second step on the road to classifying n objects is to choose a machine learning algorithm that groups the set of n points P = (p 1 , . . . , p n ) in k classes (k may be known in advance or determined by the algorithm itself). Choosing an appropriate algorithm depends on the context:

• Availability of pre-labeled data. Classifying the points P in k classes may be seen as assigning a label (such as "rose" or "tulip" in our k = 2 example) to each of the points. If one has access to some pre-labeled data, we are in the case of supervised learning: a more-or-less parametrized model is first learned from the pre-labeled data and then applied to the unlabeled points that need classification. If one does not have access to any pre-labeled data, we are in the case of unsupervised learning where classes are typically inferred only via geometrical consideration of the distribution of points in the feature space. If one has only access to a few labeled data, we are in the in-between case of semi-supervised learning where the known labels are typically propagated in one form or another in the feature space. • Inductive vs transductive learning. Another important characteristic of a classification algorithm is whether it can be used to classify only the set of points P at hand (transductive), or if it can also be directly used to classify any never-seen data point p n+1 (inductive).

This chapter focuses on the family of algorithms jointly referred to as spectral clustering. These algorithms are unsupervised and transductive: no label is known in advance and one may not naturally 1 extend the results obtained on P to never-seen data points. Another particularity of spectral clustering algorithms is that the number of classes k is known in advance. Spectral clustering algorithms have received a large attention in the last two decades due to their good performance on a wide range of different datasets, as well as their ease of implementation. In a nutshell, they combine three steps:

1. Graph construction. A sparse similarity graph is built between the n points. 2. Spectral embedding. The first k eigenvectors of a graph representative matrix (such as the Laplacian) are computed. 3. Clustering. k-means is performed on these spectral features, to obtain k clusters.

For background information about spectral clustering, such as several justifications of its performance, out-of-sample extensions, as well as comparisons with local methods, the interested reader is referred to the recent book chapter [144].

One of the drawbacks of spectral clustering is its computational cost as n, d, and/or k become large (see Sec. 2.3 for a discussion on the cost). Since the turn of the century, a large number of authors have striven to reduce the computational cost while keeping the high level of classification performance. The majority of such accelerating methods are based on sampling: they reduce the dimension of a subproblem of spectral clustering, compute a low-cost solution in small dimension, and lift the result back to the original space.

The goal of this chapter is to review existing sampling methods for spectral clustering, focusing especially on their approximation guarantees. Some of the fundamental questions we are interested in are: where is the sampling performed and what is sampled precisely? how should the reduced approximate solutions be lifted back to the original space? what is the computational gain? what is the control on performances-if it exists? Given the breadth of the literature on the subject, we do not try to be exhaustive, but rather to illustrate the key ways that sampling can be used to provide acceleration, paying special attention on recent developments on the subject.

Chapter organization. We begin by recalling in Sec. 2 the prototypical spectral clustering algorithm. We also provide some intuitive and formal justification of why it works. The next three sections classify the different methods of the literature depending on where the sampling is performed with respect to the three steps of spectral clustering:

• Sec. 3 details methods that sample directly in the original feature space.

• Sec. 4 assumes that the similarity graph is given and details methods that sample nodes and/or edges to approximate the spectral embedding. • Sec. 5 assumes that the spectral embedding is given and details methods to accelerate the k-means step.

Finally, Sec. 6 gives perspective on the limitations of existing works and discusses key open problems.

Notation. Scalars, such as λ or d, are written with low-case letters. Vectors, such as u, z or the all-one vector 1, are denoted by low-case bold letters. Matrices, such as W or L are denoted with bold capital letters. Ensembles are denoted by serif font capital letters, such as C or X. The "tilde" will denote approximations, such as in z or Ũk . We use so-called Matlab notations to slice matrices: given a set of indices S of size m and an n × n matrix W, W(S, :) ∈ R m×n is W reduced to the lines indexed by S, W(:, S) ∈ R n×m is W reduced to the columns indexed by S, and W(S, S) ∈ R m×m is W reduced to the lines and columns indexed by S. The equation U k = U(:, : k) defines U k as the reduction of U to its first k columns. Also, C is the transpose of matrix C and C + its Moore-Penrose pseudo-inverse.

The operator X = diag(x) takes as an input a vector x ∈ R n and returns an n × n diagonal matrix X featuring x in its main diagonal, i.e., X(i, j) = x(i) if i = j and X(i, j) = 0, otherwise. Finally, we will consider graphs in a large part of this paper.

We will denote by G = (V, E, W) the undirected weighted graph of |V| = n nodes interconnected by |E| = e edges: e i j ∈ E is the edge connecting nodes v i and v j , with weight W(i, j) ≥ 0. Matrix W is the adjacency matrix of G. As G is considered undirected, W is also symmetric. In general, W can be any symmetric matrix with positive entries, but we usually prefer to work with sparse graphs without self-loops, in which case the matrix is also sparse and has a zero diagonal.

Spectral clustering

The input of spectral clustering algorithms consists of (i) a set of points P = (p 1 , p 2 , . . . , p n ) (also called featured vectors) representing n objects in a feature space of dimension d, and (ii) the number of classes k in which to classify these objects. The output is a partition of the n objects in k disjoint clusters. The prototypical spectral clustering algorithm [121,102], dates back in fact to fundamental ideas by Fiedler [START_REF] Fiedler | Algebraic connectivity of graphs[END_REF] and entails the following steps:

Algorithm 1. The prototypical Spectral Clustering algorithm Input. A set of n points P = (p 1 , p 2 , . . . , p n ) in dimension d and a number of desired clusters k.

1. Graph construction (optional) a. Compute the kernel matrix K ∈ R n×n : ∀(i, j), K(i, j) = κ( p i -p j 2 ).

b. Compute W = s(K), a sparsified version of K. c. Interpret W as the adjacency matrix of a weighted undirected graph G. 2. Spectral embedding a. Compute the eigenvectors u 1 , u 2 , • • • , u k associated with the k smallest eigenvalues of a graph representative matrix R (usually a Laplacian)

computed from W. b. Set U k = [ u 1 | u 2 | • • • | u k ] ∈ R n×k .
c. Embed the i-th node to x i = U k (i,:) q( U k (i,:) 2 ) , with q(•) a normalizing function. 3. Clustering a. Use k-means on x 1 , . . . , x n in order to identify k centroids c 1 , . . . , c k . b. Voronoi tesselation: construct one cluster per centroid c and assign each object i to the cluster of the centroid closest to x i .

Output: A partition of the n points in k clusters.

A few comments are in order:

• A common choice of kernel in step 1a is the radial basis function (RBF) kernel κ( p i -p j 2 ) = exp -p i -p j 2 2 /σ2 for some user-defined sparsification function σ (step 1b). The sparsification s of K usually entails setting the diagonal to 0 and keeping only the k largest entries of each column (i.e., set all others to 0). The obtained matrix K sp is not symmetric in general and a final "symmetrization" step W = K sp + K sp is necessary to obtain a matrix W interpretable as the adjacency matrix of a weighted undirected graph 2 G = (V, E, W). This graph is called the k nearest neighbour (k-NN) similarity graph (note that the k used in this paragraph has nothing to do with the number of clusters). Other kernel functions κ and sparsification methods are possible (see Section 2 of [138] for instance).

• There are several possibilities for choosing the graph representative matrix R in step 2a. We consider three main choices [138]: Let us denote by D the diagonal degree matrix such that D(i, i) = ∑ j W(i, j) is the (weighted) degree of node v i . We define the combinatorial graph Laplacian matrix L = D-W, the normalized graph Laplacian matrix L n = I -D -1/2 WD -1/2 , and the random walk Laplacian L rw = I -D -1 W. Other popular choices include3 the non-backtracking matrix [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], degree-corrected versions of the modularity matrix [START_REF] Ali | Improved spectral community detection in large heterogeneous networks[END_REF], the Bethe-Hessian matrix [114] or similar deformed Laplacians [START_REF] Dall'amico | Optimized Deformed Laplacian for Spectrumbased Community Detection in Sparse Heterogeneous Graphs[END_REF]. • The normalizing function q(•) used in step 2c depends on which representative matrix is chosen. In the case of the Laplacians, experimental evidence as well as some theoretical arguments [138] support using a unit norm normalization for the eigenvectors of L n (i.e. q is the identity function), and no normalization for the eigenvectors of L and L rw (i.e. q(•) = 1). • Step 1 of the algorithm is "optional" in the sense that in some cases the input is not a set of points but directly a graph. For instance, it could be a graph representing a social network between n individuals, where each node is an individual and there is an edge between two nodes if they know each other. The weight on each edge can represent the strength of their relation (for instance close to 0 if they barely know each other, and close to 1 if they are best friends). The goal is then to classify individuals based on the structure of these social connections and is usually referred to as community detection in this context [START_REF] Fortunato | Community detection in graphs[END_REF]. Given the input graph, and the number k of communities to identify, one can run spectral algorithms starting directly at step 2. Readers only interested in such applications can skip Sec. 3, which is devoted to sampling techniques designed to accelerate step 1.

After the spectral embedding X = ( x 1 , . . . , x n ) has been identified, spectral clustering uses k-means in order to find the set of k centroids C = ( c 1 , . . . , c k ) that best represents the data. Formally, the k-means cost function to minimize reads:

f (C; X) = ∑ x∈X min c∈C x -c 2 2 . (1)
We would ideally hope to identify the set of k centroids C * minimizing f (C; X).

Solving exactly this problem is NP-hard [START_REF] Drineas | Clustering in Large Graphs and Matrices[END_REF], so one commonly resorts to approximation and heuristic solutions (see for instance [128] for details on different such heuristics). The most famous is the so-called Lloyd-Max heuristic algorithm:

Algorithm 2. The Lloyd-Max algorithm [START_REF] Lloyd | Least squares quantization in PCM[END_REF] Input. Set of n points X = (x 1 , x 2 , . . . , x n ) and number of desired clusters k.

1. Start from an initial guess C ini of k centroids 2. Iterate until convergence: a. Assign each point x i to its closest centroid to obtain a partition of X in k clusters. b. Move each centroid c to the average position of all points in cluster .

Output: A set of k centroids C = (c 1 , . . . , c k ).
When the clusters are sufficiently separated and C ini is not too far from the optimal centroids, then the Lloyd-Max algorithm converges to the correct solution [START_REF] Kumar | Clustering with spectral norm and the k-means algorithm[END_REF]. Otherwise, it typically ends up in a local minimum.

A remark on notation. Two quantities of fundamental importance in spectral clustering are the eigenvalues λ i and especially the eigenvectors u i of the graph Laplacian matrix. We adopt the graph theoretic convention of sorting eigenvalues in nondecreasing order: 0 = λ 1 ≤ λ 2 ≤ . . . ≤ λ n . Also, for reasons of brevity, we overload notation and use the same symbol for the spectrum of the three Laplacians L, L n and L rw . Thus, we advise the reader to rely on the context in order to discern which Laplacian gives rise to the eigenvalues and eigenvectors. Finally, the reader should keep in mind that the largest eigenvalue if always bounded by 2 for L n and L rw .

An illustration of spectral clustering

The first two steps of the algorithm can be understood as a non-linear transformation from the initial feature space to another feature space (that we call spectral feature space or spectral embedding): a transformation of features p i in R d to spectral features x i in R k . The first natural question that arises is why do we run k-means on the spectral features X = (x 1 , . . . , x n ) that are subject to parameter tuning and costly to compute, rather than directly run k-means on the original P? Figures 1 and2 illustrate the answer.

In Figure 1, we show the result of k-means directly on a set of artificial features P known as the two-half moons dataset. In this example, the intuitive ground truth is that each half-moon corresponds to a class, that we want to recover. Running kmeans directly in this 2D feature space will necessarily output a linear separation between the two obtained Voronoi cells and will thus necessarily fail, as no straight line can separate the two half-moons. Spectral clustering, via the computation of the spectral features of a similarity graph, transforms these original features P in spectral features X that are typically linearly separable by k-means: the two half-moons are successfully recovered! We illustrate this in Figure 2. In the next section, we will examine a theoretical argument aiming to justify this phenomenon.

Justification of spectral clustering

A popular approach -and by no means the only one, see Sec. 2.2.3-to justify spectral clustering algorithms stems from its connection to graph partitioning. Suppose Fig. 2 Illustration of the spectral clustering algorithm on the two half-moons dataset (n = 500, d = 2, k = 2). The graph is created with a RBF kernel and via a sparsification done with k-nearest neighbours (with k = 5). The spectral embedding is done with the two eigenvectors associated to the two smallest eigenvalues of the combinatorial Laplacian matrix L. The embedding X is here in practice in 1D as the first eigenvector of L is always constant and thus not discriminative (to confirm this, first show that L is a PSD matrix and then prove that Lc = 0 for any constant vector c). Observe how the two clusters are now linearly separable in the spectral feature space. k-means on these features successfully recovers the two half-moons.

that the similarity graph G = (V, E, W) has been obtained and we want to compute a partition 4

P = {V 1 , V 2 , . . . , V k } of the nodes V in k groups.
Intuitively, a good clustering objective function should favor strongly connected nodes to end up in the same subset, and nodes that are far apart in the graph to end up in different subsets. This intuition can be formalized with graph cuts.

Considering two groups V 1 and

V 2 , define w(V 1 , V 2 ) = ∑ i∈V 1 ∑ j∈V 2 W(i, j
) to be the total weight of all links connecting V 1 to V 2 . Also, denote by V the complement of V in V, such that w(V , V ) is the total weight one needs to cut in order to disconnect V from the rest of the graph. Given these definitions, the simplest graph cut objective function, denoted by cut, is:

cut (P = {V 1 , . . . , V k }) = 1 2 k ∑ =1 w(V , V ). (2)
The best partition according to the cut criterion is P * = argmin P cut(P). For k = 2, solving this problem can be done exactly in O(ne + n 2 log(n)) amortized time using the Stoer-Wagner algorithm [126] and approximated in nearly linear time [START_REF] Karger | Minimum cuts in near-linear time[END_REF].

Nevertheless, this criterion is not satisfactory as it often separates an individual node from the rest of the graph, with no attention to the balance of the sizes or volumes of the groups. In clustering, one usually wants to partition into groups that are "large enough". There are two famous ways to balance the previous cost in the machine learning literature5 : the ratio cut [143] and normalized cut [121] cost functions, respectively defined as:

rcut(P) = 1 2 k ∑ =1 w(V , V ) |V | and ncut(P) = 1 2 k ∑ =1 w(V , V ) vol(V ) , (3)
where |V | is the number of nodes in V and vol(V ) = ∑ i∈V ∑ j∈V W(i, j) is the so-called volume of V . The difference between them is that ncut favors clusters of large volume, whereas rcut only considers cluster size-though for a d-regular graph with unit weights the two measures match (up to multiplication by 1/d). Unfortunately, it is hard to minimize these cost functions directly: minimizing these two balanced costs is NP-hard [139, 121] and one needs to search over the space of all possible partitions which is of exponential size.

A continuous relaxation. Spectral clustering may be interpreted as a continuous relaxation of the above minimization problems. Without loss of generality, in the following we concentrate on relaxing the rcut minimization problem (ncut is relaxed almost identically). Given a partition P = (V 1 , . . . , V k ), let us define

C = z 1 |V 1 | | . . . | z k |V k | ∈ R n×k , (4)
where z ∈ R n is the indicator vector of V :

z (i) = 1 if node i ∈ V , 0 otherwise. (5)
It will prove useful in the following to remark that, independently of how the partitions are chosen, we always have that C C = I, the identity matrix in dimension k. With this in place, the problem of minimizing rcut can be rewritten as (see discussion in

To understand why this equivalence holds, one should simply note that

tr C LC = k ∑ =1 1 |V | z Lz = k ∑ =1 1 |V | ∑ i> j W(i, j)(z (i) -z ( j)) 2 = k ∑ =1 w(V , V ) |V | = 2 rcut(P).
Solving ( 6) is obviously still NP-hard as the only thing we have achieved is to rewrite the rcut minimization problem in matrix form. Yet, in this form, it is easier to realize that one may find an approximate solution by relaxing the discreteness constraint "C as in (4)". In the absence of the hard-to-deal-with constraint, the relaxed problem is not only polynomially solvable but also possesses a closed-form solution! By the Courant-Fischer-Weyl (min-max) theorem, the solution is given by the first k eigenvectors

U k = [u 1 , u 2 , . . . , u k ] of L: U k = arg min C∈R n×k tr C LC subject to C C = I.
This relaxation is not unique to the combinatorial Laplacian. In the same spirit, the minimum ncut optimization problem can be formulated in terms of the normalized Laplacian matrix L n , and the relaxed problem's solution is given by the first k eigenvectors of L n .

A difficulty still lies before us: how do we go from a real-valued U k to a partition of the nodes? The two next subsections aim to motivate the use of k-means as a rounding heuristic. The exposition starts from the simple case when there are only two clusters (k = 2) before considering the general case (arbitrary k).

The case of two clusters: thresholding suffices

For simplicity, we first consider the case of two clusters. If one constructs a partitioning P t with V 1 = {v i : u 2 (i) > t} and V 2 = {v i : u 2 (i) ≤ t} for every level set t ∈ (-1, 1), then it is a folklore result that

rcut(P * ) ≤ min t rcut(P t ) ≤ 2 rcut(P * ) d max - λ 2 2 , (7)
with P * = arg min P rcut(P) being the optimal partitioning, d max is the maximum degree of any node in V, and λ 2 the second smallest eigenvalue of L. The upper bound is achieved by the tree-cross-path graph constructed by Guattery and Miller [START_REF] Guattery | On the quality of spectral separators[END_REF]. In an analogous manner, if P * = arg min P ncut(P) is the optimal partitioning w.r.t. the ncut cost and every P t has been constructed by thresholding the second eigenvector of L n , then

ncut(P * ) ≤ min t ncut(P t ) ≤ 2 ncut(P * ). (8)
Inequality ( 8) can be derived as a consequence of the Cheeger inequality, a key result of spectral graph theory [START_REF] Chung | Spectral graph theory[END_REF], which for the normalized Laplacian reads:

λ 2 2 ≤ ncut(P * ) ≤ min V w(V, V) min{w(V), w( V)} ≤ min t ncut(P t ) ≤ 2λ 2 .
As a consequence, we have

ncut(P * ) ≤ min t ncut(P t ) ≤ 2λ 2 ≤ 4ncut(P * ) = 2 ncut(P * ),
as desired. The derivation of the rcut bound given in [START_REF] Arya | An Optimal Algorithm for Approximate Nearest Neighbor Searching Fixed Dimensions[END_REF] follows similarly.

More than two clusters: use k-means

As the number of clusters k increases, the brute-force approach of testing every level set becomes quickly prohibitive. But why is k-means the right way to obtain the clusters in the spectral embedding? Though a plethora of experimental evidence advocate the use of k-means, a rigorous justification is still lacking. The interested reader may refer to [START_REF] Lee | Multiway spectral partitioning and higher-order cheeger inequalities[END_REF] for an example of an analysis of spectral partitioning without k-means.

More recently, Peng et al.

[107] came up with a mathematical argument showing that, if G is well clusterable and we use a k-means algorithm (e.g., [START_REF] Kumar | A simple linear time (1+ε)-approximation algorithm for k-means clustering in any dimensions[END_REF]) which guarantees that the identified solution C abides to

f ( C; X) ≤ (1 + ε) f (C * ; X),
where C * is the optimal solution of the k-means problem, then the partitioning P produced by spectral clustering when using L n has ncut cost provably close to that of the optimal partitioning P * . In particular, it was shown that, as long as λ k+1 ≥ ck 2 ncut(P * ), then

ncut(P * ) ≤ ncut( P) ≤ ζ ncut(P * ) 1 + ε k 3 λ k+1 ,
for some constants c, ζ > 0 that are independent of n and k (see also [START_REF] Kolev | A note on spectral clustering[END_REF]). Note that, using the higher-order Cheeger inequality [START_REF] Lee | Multiway spectral partitioning and higher-order cheeger inequalities[END_REF] λ k /2 ≤ ncut(P * ), the condition λ k+1 ≥ ck 2 ncut(P * ) implies

λ k+1 λ k ≥ ck 2 2 = Ω (k 2 ).
Though hopefully milder than this one 6 , such gap assumptions are very common in the analysis of spectral clustering. Simply put, the larger the gap λ k+1λ k is, the stronger the cluster structure and the easier it is to identify a good clustering. Besides quantifying the difficulty of the clustering problem, the gap also encodes the robustness of the spectral embedding to errors induced by approximation algorithms [START_REF] Davis | The rotation of eigenvectors by a perturbation[END_REF]. The eigenvectors of a perturbed Hermitian matrix exhibit an interesting property: instead of being arbitrary, induced changes are localized w.r.t. the eigenvalue axis, following an inverse square eigenvalue-distance law [START_REF] Loukas | How close are the eigenvectors of the sample and actual covariance matrices[END_REF]. More precisely, if ũi is the i-th eigenvector after perturbation, then the inner products ( ũ i u j ) 2 decrease proportionally with |λ iλ j | 2 . As such, demanding that λ k+1λ k is large is often helpful in the analysis of spectral clustering algorithms in order to ensure that the majority of useful information (contained within U k ) is preserved (in Ũk ) despite approximation errors 7 .

Choice of relaxation

The presented relaxation approach is not unique and other relaxations could be equally valid (see for instance [START_REF] Bie | Fast sdp relaxations of graph cut clustering, transduction, and other combinatorial problems[END_REF][START_REF] Bresson | Multiclass Total Variation Clustering[END_REF]112]). This relaxation has nevertheless the double advantage of being theoretically simple and computationally easy to implement. Also, justification of spectral clustering algorithms does not only come from this graph cut perspective and in fact encompasses several approaches that we will not detail here: perturbation approaches or hitting time considerations [138], a polarization theorem [START_REF] Brand | A unifying theorem for spectral embedding and clustering[END_REF], consistency derivations [135,[START_REF] Lei | Consistency of spectral clustering in stochastic block models[END_REF], etc. Interestingly, recent studies (for instance [START_REF] Bordenave | Non-backtracking Spectrum of Random Graphs: Community Detection and Non-regular Ramanujan Graphs[END_REF]) on the Stochastic Block Models have shown that spectral clustering (on other matrices than the Laplacian, such as the non-backtracking matrix [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], or the Bethe-Hessian matrix [114] or other similar deformed Laplacians [START_REF] Dall'amico | Optimized Deformed Laplacian for Spectrumbased Community Detection in Sparse Heterogeneous Graphs[END_REF]) perform well up to the detectability threshold of the block structure.

Computational complexity considerations

What is the computational complexity of spectral clustering as a function of the number of points n, their dimension d and the number of desired clusters k? Let us examine the three steps involved one by one.

The first step entails the construction of a sparse similarity graph from the input points, which is dominated by the kernel computation and costs O(dn 2 ). In the second step, given the graph G consisting of n nodes and e edges 8 , one needs to compute the spectral embedding (step 2 of Algorithm 1). Without exploiting the special structure of a graph Laplacian -other than its sparsity that is-there are two main options:

• Using power iterations, one may identify sequentially each non-trivial eigenvector u in time O(e/δ ), where δ = λλ -1 is the -th eigenvalue gap and e is the number of edges of the graph [136]. Computing the spectral embedding there-fore takes O(ke/δ ) with δ = min δ . Unfortunately, there exist graphs 9 such that δ = O(1/n), bringing the overall worst-case complexity to O(kne). • The Lanczos method can be used to approximate the first k eigenvectors in roughly O(ek + nk 2 ) time. This procedure is often numerically unstable resulting to a loss of orthogonality in the computed Krylov subspace basis. The most common way to circumvent this problem is by implicit restart [START_REF] Calvetti | An implicitly restarted lanczos method for large symmetric eigenvalue problems[END_REF], whose computational complexity is not easily derived. The number of restarts, empirically, depend heavily on the eigenvalue distribution in the vicinity of λ k : if λ k is in an eigenvalue bulk, the algorithms takes longer than when λ k is isolated. We decide to write the complexity of restarted Arnoldi as O(t(ek + nk 2 )) with t modeling the number of restarts. Note that throughout this paper, t will generically refer to a number of iterations in algorithm complexities. We refer the interested reader to [START_REF] Bai | Templates for the solution of algebraic eigenvalue problems: a practical guide[END_REF] for an in-depth discussion of Lanczos methods.

The third step entails solving the k-means problem, typically by using the Lloyd-Max algorithm to converge to a local minimum of f (C; X). Since there is no guarantee that this procedure will find a good local minimum, it is usually rerun multiple times, starting in each case from randomly selected centroids C ini . The computational complexity of this third step is O(tnk 2 ), where t is a bound on the number of iterations required until convergence multiplied by the number of retries (typically 10).

A taxonomy of sampling methods for spectral clustering

For the remainder of the chapter, we propose to classify sampling methods aiming at accelerating one or more of these three steps according to when they sample. If they sample before step 1, they are detailed in Sec. 3. Methods that assume that the similarity graph is given or well-approximated and sample between steps 1 and 2 will be found in Sec. 4. Finally, methods that assume that the spectral embedding has been exactly computed or well-approximated and sample before the k-means step are explained in Sec. 5. This classification of methods, like all classification systems, bears a few flaws. For instance, Nyström methods can be applied to both the context of Sections 3 and 4 and are thus mentioned in both. Also, we decided to include the pseudo-code of only a few chosen algorithms that we think are illustrative of the literature. This choice is of course subjective and debatable. Notwithstanding these flaws, we hope that this taxonomy clarifies the landscape of existing methods.

Sampling in the original feature space

This section is devoted to methods that ambitiously aim to reduce the dimension of the spectral clustering problem even before the graph has been formed. Indeed, the naive way of building the similarity graph (step 1 of spectral clustering algorithms) costs O(dn 2 ) and, as such, is one of the the main computational bottlenecks of spectral clustering. It should be remarked that the present discussion fits into the wider realm of kernel approximation, a proper review of which cannot fit in this chapter: we will thus concentrate on methods that were in practice used for spectral clustering.

Nyström-based methods

The methods of this section aim to obtain an approximation Ũk of the exact spectral embedding U k via a sampling procedure in the original feature space.

The Nyström method is a well known algorithm for obtaining a rough low rank approximation of a positive semi-definite (PSD) matrix A.

Here is a high level description of the steps entailed:

Algorithm 3. Nyström's method Input. PSD matrix A ∈ R n×n ,
B k = Q k Σ k Q k , where Q k = Q(:, : k) ∈ R n×k and Σ k = Σ(: k, : k).
Possible outputs:

• A low-rank approximation à = CB + C ∈ R n×n of A • A rank-k approximation Ãk = CB + k C ∈ R n×n of A • The top k eigenvectors of Ãk , stacked as columns in matrix Ṽk ∈ R n×k , obtained by orthonormalizing the columns of Qk = CQ k Σ -1 k ∈ R n×k
Various guarantees are known for the quality of à depending on the type of sampling utilized (i.e., how the indices in S are selected in step 1) and the preferred notion of error (spectral . 2 vs frobenius . F vs trace . * norm) [START_REF] Gittens | Revisiting the nyström method for improved large-scale machine learning[END_REF][START_REF] Kumar | Sampling techniques for the nystrom method[END_REF][START_REF] Frieze | Fast monte-carlo algorithms for finding low-rank approximations[END_REF]148]. For instance:

Theorem 1 (Lemma 8 for q = 1 in [START_REF] Gittens | Revisiting the nyström method for improved large-scale machine learning[END_REF]). Let ε ∈ (0, 1) and δ ∈ (0, 1) and suppose that S contains the indices of m columns drawn i.i.d. uniformly at random (with or without replacement). Then:

A -Ã 2 ≤ 1 + n (1 -ε)m A -A k 2
holds with probability at least 1 -3δ , provided that m ≥ 2ε -2 µk log (k/δ ); where

µ = n k max i=1,...,n V k (i, :) 2 2
is the coherence associated with the first k eigenvectors V k of A, and A k is the best rank-k approximation of A.

Guarantees independent of the coherence can be obtained for more advanced sampling methods. Perhaps the most well known method is that of leverage scores, where one draws m samples independently by selecting (with replacement) the i-th column with probability p i = V k (i, :) 2 2 /k. Theorem 2 (Lemma 5 for q = 1 in [START_REF] Gittens | Revisiting the nyström method for improved large-scale machine learning[END_REF]). Let ε ∈ (0, 1) and δ ∈ (0, 1) and suppose that S contains the indices of m columns drawn i.i.d. with replacement from such a probability distribution. Then:

A -Ã 2 ≤ A -A k 2 + ε 2 A -A k * holds with probability at least 0.8 -2δ provided that m ≥ O(ε -2 k log(k/δ )).
Computing leverage scores exactly is computationally prohibitive since it necessitates a partial SVD decomposition of A, which we are trying to avoid in the first place. Nevertheless, it is possible to approximate all leverage scores with a multiplicative error guarantee in time roughly O(ek log(e)) if A has O(e) non-zero entries. (see Algorithms 1 to 3 in [START_REF] Gittens | Revisiting the nyström method for improved large-scale machine learning[END_REF]). Many variants of the above exist [START_REF] Kumar | Sampling techniques for the nystrom method[END_REF][START_REF] Kumar | Sampling methods for the nyström method[END_REF], but to the best of our knowledge, the fastest current Nyström algorithm utilizes ridge leverage scores with a complex recursive sampling scheme and runs in time nearly linear in n [100].

Nyström for spectral clustering. Though initially conceived for low-rank approximation, Nyström's method can also be used to accelerate spectral clustering. The key observation is that U k , the tailing k eigenvectors of the graph representative matrix R, can be interpreted as the top k eigenvectors of the PSD matrix A = R 2 I -R. As such, the span of the k top eigenvectors of Ãk obtained by running Algorithm 3 on A is an approximation of the span of the exact spectral embedding. Different variants of this idea have been considered for the acceleration of spectral clustering [START_REF] Fowlkes | Spectral grouping using the Nystrom method[END_REF]141,[START_REF] Li | Time and space efficient spectral clustering via column sampling[END_REF][START_REF] Bouneffouf | Sampling with minimum sum of squared similarities for nystrombased large scale spectral clustering[END_REF][START_REF] Mohan | Exploiting sparsity to improve the accuracy of nyström-based large-scale spectral clustering[END_REF][START_REF] Li | Towards Large Scale Spectral Problems via Diffusion Process[END_REF].

Following our taxonomy, we hereby focus on the case where we have at our disposal n points p i in dimension d, and the similarity graph has yet to be formed. The case where the graph is known is deferred to Sec. 4.

In this case, we cannot run Algorithm 3 on A = R 2 I -R as the graph, and a fortiori its representative matrix R, has not yet been formed. What we can have access to efficiently is B = s(K(S, S)) and C = s(K(:, S)), as these require only a partial computation of the kernel and cost only O(dnm). Note that s is a sparsification function that is applied on a subset of the kernel matrix.

The following pseudo-code exemplifies how Nyström-based techniques can be used to approximate the first k eigenvectors U k associated with the normalized Laplacian matrix (i.e., here R = L n ): Algorithm 3b. Nyström for spectral clustering [START_REF] Li | Time and space efficient spectral clustering via column sampling[END_REF] Input. The set of points P, the number of desired clusters k, a sampling set S of size m ≥ k

1. Compute the sub-matrices B = s(K(S, S)) ∈ R m×m and C = s(K(:, S)) ∈ R n×m , where s is a sparsification function. 2. Let D r = diag(B1) be the m × m degree matrix. 3. Compute the top k eigenvalues Σ k and eigenvectors Q k of D -1 /2 r BD -1 /2 r . 4. Set Qk = CD -1 /2 r Q k Σ -1 k . 5. Let D l = diag( Qk Σ k Q k 1) be the n × n degree matrix. 6. Compute Ũk obtained by orthogonalizing D -1 /2 l

Qk .

Output: Ũk , an approximation of the spectral embedding U k . This algorithm runs in O(nm max(d, k)) time, which is small when m depends mildly on the other parameters of interest. Nevertheless, the algorithm (and others like it) suffers from several issues:

• Alg. 3b attempts to use Nyström's method on

A = 2I -L n = I + D -1 2 s(K)D -1 2
via the exact computation of two sub-matrices of K. In doing so, it makes two strong (and uncontrolled) approximations: First of all, the sparsification step (step 1 in Alg. 3b) is applied to the sub-matrices K(S, S) and K(:, S), deviating from the correct sparsification procedure that takes into account the entire kernel matrix K. Second, the degree matrix D is never exactly computed as knowing it exactly would entail computing exactly s(K), which is precisely what we are trying to avoid. Existing methods thus rely on heuristic approximations of the degree in order to bypass this difficulty (see steps 2 and 5 of Alg. 3b). • Since we don't have direct access to the kernel matrix, we cannot utilize advanced sampling methods such as leverage scores to draw the sampling set S. This is particularly problematic if (due to sparsification), matrices B and C are sparse, as for sparse matrices uniform sampling is known to perform poorly [START_REF] Mohan | Exploiting sparsity to improve the accuracy of nyström-based large-scale spectral clustering[END_REF].

Techniques that rely on distances between columns do not fair much better. Landmark-based approaches commonly perform better in simple problems but suffer when the clusters are non-convex [START_REF] Bouneffouf | Sampling with minimum sum of squared similarities for nystrombased large scale spectral clustering[END_REF]. We refer the reader to the work by Mohan et al. [START_REF] Mohan | Exploiting sparsity to improve the accuracy of nyström-based large-scale spectral clustering[END_REF] for more information on landmark-based methods. The latter work also describes an involved sampling scheme that is aimed at general (i.e., non-convex) clusters.

For the reasons highlighted above, the low-rank approximation guarantees accompanying the classical Nyström method cannot be directly used here. A fortiori, it is an open question how much the quality of the spectral clustering solution is affected by using the centroids obtained by running k-means on Ũk .

Column sampling. Akin in spirit to Nyström methods, an alternative approach to accelerating spectral clustering was inspired by column sampling low-rank approximation techniques [START_REF] Drineas | Fast monte carlo algorithms for matrices ii: Computing a low-rank approximation to a matrix[END_REF][START_REF] Deshpande | Matrix approximation and projective clustering via volume sampling[END_REF].

An instance of such algorithms was put forth under the name of cSPEC (column sampling spectral clustering) by Wang et al. [141].

Let C = U C Σ C V C be the singular value decomposition of the n × m matrix C = s(K(:, S)). Then, matrices Σ = n m Σ C and Ũ = CV C Σ + C
are interpreted as an approximation of the actual eigenvalues and eigenvectors of K and thus U k can be substituted by the first k columns of Ũ. This algorithm runs in O(ndm + nm 2 ).

Authors in [START_REF] Chen | Large Scale Spectral Clustering with Landmark-Based Representation[END_REF] propose a hybrid method, between column sampling and the representative-based methods discussed in Sec. 3.3, where they propose the following approximate factorization of the data matrix:

(p 1 | . . . |p n ) FZ ∈ R d×n , (9)
where F ∈ R d×m concatenates the feature vectors of m sampled points and Z ∈ R m×n represents all unsampled points as approximate linear combinations of the representatives, computed via sparse coding techniques [START_REF] Lee | Efficient sparse coding algorithms[END_REF] 10 . The SVD of D-1/2 Z, with D the row-sum of Z, is then computed to obtain an approximation Ũk of U k . The complexity of their algorithm is also O(ndm + nm 2 ). In these methods, the choice of the sample set S is, of course, central and has been much debated. Popular options are uniformly at random or via better-tailored probability distributions, via a first k-means (with k = m) pass on P, or via other selective sampling methods. Also, as with most extensions of Nyström's method to spectral clustering, column sampling methods for spectral clustering do not come with end-to-end approximation guarantees on U k .

In the world of low-rank matrix approximation the situation is somewhat more advanced. Recent work in column sampling utilizes adaptive sampling with leverage scores in time O(e + npoly(k)), or uniformly i.i.d. after preconditioning by a fast randomized Hadamard transform [145,[START_REF] Drineas | Lectures on randomized numerical linear algebra[END_REF]. Others have also used a correlated version called volume sampling to obtain column indices [START_REF] Deshpande | Matrix approximation and projective clustering via volume sampling[END_REF]. Nevertheless, this literature extends beyond the scope of this chapter and thus we invite the interested reader to consider the aforementioned references for a more in-depth perspective.

Random Fourier features

Out of several sketching techniques one could a priori use to accelerate spectral clustering, we focus on random Fourier features (RFF) [110]: a method that samples in the Fourier space associated to the original feature space. Even though RFFs have originally been developed to approximate a kernel matrix K in time linear in n instead of the quadratic time necessary for its exact computation, they can in fact be used to obtain an approximation Ũk of the exact spectral embedding U k .

Let us denote by κ the RBF kernel, i.e., κ(t) = exp(-t 2 /σ 2 ), whose Fourier transform is:

κ(ω) = R d κ(t) exp -iω t dt. ( 10
)
The above takes real values as κ is symmetric. One may write:

κ(p, q) = κ(p -q) = 1 Z R d κ(ω) exp iω (p-q) dω, (11)
where, in order to ensure that κ(p, p) = 1, the normalization constant is set to Z = R d κ(ω)dω. According to Bochner's theorem, and due to the fact that κ is positivedefinite, κ/Z is a valid probability density function. κ(p, q) may thus be interpreted as the expected value of exp iω (p-q) provided that ω is drawn from κ/Z:

κ(p, q) = E ω exp iω (p-q) ( 12
)
Drawing ω from the distribution κ/Z is equivalent to drawing independently each of its d entries according to the normal law of mean 0 and variance 2/σ 2 . Indeed:

κ(ω) = π d/2 σ d exp(-σ 2 ω 2 /4) and Z = R d κ(ω)dω = (2π) d , leading to κ(ω) Z = σ 2 √ π d exp -σ 2 ω 2 /4 .
In practice, we draw independently m such vectors ω to obtain the set of sampled frequencies Ω = (ω 1 , . . . , ω m ). For each data point p i , and given this set of samples Ω , we define the associated random Fourier feature vector:

ψ i = 1 √ m [cos(ω 1 p i )| • • • | cos(ω m p i )| sin(ω 1 p i )| • • • | sin(ω m p i )] ∈ R 2m , ( 13
)
and call 11 , degree-normalize Ψ to obtain a low-rank approximation of the normalized Laplacian L n and perform an SVD to directly obtain an approximation Ũk of the spectral embedding U k . The total cost to obtain this approximation is O(ndm + nm 2 ). These ideas were developed in Refs. [START_REF] Chitta | Efficient Kernel Clustering Using Random Fourier Features[END_REF]146] for instance.

Ψ = (ψ 1 | • • • |ψ n ) ∈
As in Nyström methods however, the concentration guarantees of RFFs for K do not extend to the degree-normalized case; moreover, the sparsification step 1b of spectral clustering is ignored. Note that improving over RFFs in terms of efficiency and concentration properties is the subject of recent research (see for instance [START_REF] Le | Fastfood-approximating kernel expansions in loglinear time[END_REF]).

The paradigm of representative points

The methods detailed here sample in the original feature space and directly obtain a control on the misclustering rate due to the sampling process. They are based on the following framework:

1. Sample m so-called representatives. b. a correspondence table to associate each p i to its nearest representative 2. Run spectral clustering on Y to get the cluster membership of each y i . 3. Lift the cluster membership to each p i by looking up the cluster membership of its representative in the correspondence table.

Output: k clusters 11 an approximation of the degree d i of node v i is ψ i ψ where ψ = ∑ j ψ j . All degrees can thus be estimated in time O(nm 2 ).

The complexity of KASP is bounded by 12 O(mdnt + m 3 ). For a summary of the analysis given in [147], let us consider the cluster memberships given by exact spectral clustering on P as well as the memberships given by exact spectral clustering on P = (p 1 + ε 1 , . . . , p n + ε n ) where the ε i are any small perturbations on the initial points. Let us denote by L n (resp. Ln ) the normalized Laplacian matrix of the similarity graph on P (resp. P). The analysis concentrates on the study of the miss-clustering rate ρ:

ρ = # of points with different memberships n . ( 14
)
The main result, building upon preliminary work in [START_REF] Huang | Spectral clustering with perturbed data[END_REF], stems from a perturbation approach and reads:

Theorem 3. Under the assumptions of Theorem 3 in [147]: ρ ≤ O k g 2 0 L n -Ln F ,
where g 0 is a value depending on the spectral gap. Also, under the assumptions of Theorem 6 in [147], one has, with high probability:

L n -Ln F ≤ O σ (2) ε + σ (4) ε , (15) with σ (2)
ε and σ

ε the 2nd and 4th moments of the perturbation's norms ε i 2 . Combining both bounds, one obtains an upper bound on the misclustering rate that depends on the second and fourth moments of the perturbation's norms ε i 2 . The "collapse" of points onto the m representative points, interpreted as a perturbation on the original points, should thus tend to minimize these two moments, leading the authors to propose distortion-minimizing algorithms, such as KASP. A very similar algorithm, eSPEC, is described in [141].

Other methods

Approximate nearest neighbour search algorithms

The objective here is to approximate the nearest neighbour graph efficiently. Even though these methods are not necessarily based on sampling, we include them in the discussion as they are frequently used in practice.

Given the feature vectors p 1 , . . . , p n ∈ R d and a query point q ∈ R d , the exact nearest neighbour search (exact NNS) associated to P and q is p * = argmin p∈P dist(q, p) where dist stands for any distance. Different distances are possible depending on the choice of kernel κ. We will here consider the Euclidean norm as it enters the definition of the classical RBF kernel. Computing the exact NNS costs O(dn). The goal of the approximate NNS field of research is to provide faster algorithms that have the following control on the error.

Definition 1. Point p * is an ε-approximate nearest neighbor of query q ∈ R d , if ∀p ∈ P dist(q, p * ) ≤ (1 + ε) dist(q, p).
For ε = 0, this reduces to exact NNS.

Extensions of this objective to the k-nearest neighbour goal are considered in the NNS literature. A k-nearest neighbour graph can then be constructed simply by running an approximate k-NNS query for each object p i . Thus, approximate NSS algorithms are interesting candidates to approximate the adjacency matrix of the nearest-neighbour affinity graph, that we need in step 1 of spectral clustering. Many algorithms exist, their respective performances depending essentially on the dimension d of the feature vectors. According to [START_REF] Avrithis | High-dimensional approximate nearest neighbor: kd Generalized Randomized Forests[END_REF], randomized k-d forests as implemented in the library FLANN [START_REF] Muja | Scalable nearest neighbor algorithms for high dimensional data[END_REF] are considered state-of-the-art for dimension of around 100, whereas methods based on Balanced Box Decomposition (BBD) [START_REF] Arya | An Optimal Algorithm for Approximate Nearest Neighbor Searching Fixed Dimensions[END_REF][START_REF] Anagnostopoulos | Low-Quality Dimension Reduction and High-Dimensional Approximate Nearest Neighbor[END_REF] are known to perform well for d roughly smaller than 100. In high dimensions, to avoid the curse of dimensionality, successful approaches are for instance based on hashing methods (such as Locality Sensitive Hashing (LSH) [START_REF] Andoni | Near-Optimal Hashing Algorithms for Approximate Nearest Neighbor in High Dimensions[END_REF], Product Quantization (PQ) [START_REF] Kalantidis | Locally Optimized Product Quantization for Approximate Nearest Neighbor Search[END_REF]) or k-d generalized random forests [START_REF] Avrithis | High-dimensional approximate nearest neighbor: kd Generalized Randomized Forests[END_REF]. Finally, proximity graph methods, that sequentially improve over a first coarse approximation of the k-NN graph (or other graph structures such as navigable graphs) have received a large attention recently and are becoming state-of-the-art in regimes where quality of approximation primes (see for instance [START_REF] Malkov | Efficient and robust approximate nearest neighbor search using Hierarchical Navigable Small World graphs[END_REF][START_REF] Dong | Efficient k-nearest neighbor graph construction for generic similarity measures[END_REF][START_REF] Fu | Fast Approximate Nearest Neighbor Search With The Navigating Spreading-out Graph[END_REF][START_REF] Aumüller | ANN-Benchmarks: A Benchmarking Tool for Approximate Nearest Neighbor Algorithms[END_REF]). Such tools come with various levels of guarantees and computation costs, the details of which are not in the scope of this chapter. Experimentally, to obtain an approximate k-NN graph with a typical recall rate13 of 0.9, these algorithms are observed to achieve a complexity of O(dn α ) with α close to 1 (α 1.1 in [START_REF] Dong | Efficient k-nearest neighbor graph construction for generic similarity measures[END_REF] for instance).

Feature selection and feature projection

Some methods work on reducing the factor d of the complexity O(dn 2 ) of the kernel computation via feature selection, i.e., the sampling of features deemed more useful for the underlying clustering task, or feature projection, i.e., the projection on usually random subspaces of dimension d < d. Feature selection methods are usually designed to improve the classification by removing features that are too noisy or useless for the classification. We thus do not detail further these methods as they are not approximation algorithms per se. The interested reader will find some entries in the literature via references [START_REF] Dash | Feature Selection for Clustering[END_REF][START_REF] He | Laplacian score for feature selection[END_REF]149,[START_REF] Cai | Unsupervised feature selection for multi-cluster data[END_REF]. Projection methods use random projections of the original points P on spaces of dimension d ∼ log n in order to take advantage of the Johnson-Lindenstrauss lemma of norm conservation: the kernel computed from the projected features in dimension d is thus an approximation of the true kernel with high probability. We refer to the works [116,[START_REF] Hunter | Compressive Spectral Clustering[END_REF] for more details.

Sampling given the similarity graph

We now suppose that the similarity graph is either given (e.g., in cases where the original data is a graph) or has been well approximated (by approximate k-NN search for instance) and concentrate on sampling-based methods that aim to reduce the cost of computing the first k eigenvectors of R.

These methods predominantly aim to approximate R by a smaller matrix R of size m. The eigen-decomposition is done in R m which can be significantly cheaper when m n. In addition, each method comes with a fast way of lifting vectors from R m back to R n (this is usually a linear transformation). After lifting, the eigenvectors of R are used as a proxy for those of R.

Unlike the previous section where a strong approximation guarantee of the exact embedding U k by an efficiently computed Ũk was a distant and difficult goal to achieve in itself; we will see in this section that the knowledge of the similarity graph not only enables to obtain such strong approximation guarantees, but also enables to control how the error on U k transfers as an error on the k-means cost.

To be more precise, recall (1) defining the k-means cost f (C; X) associated to the n points X = (x 1 , . . . , x n ) and a centroid set C. Now, suppose that we have identified a set of n points X = ( x1 | . . . | xn ) that are meant to approximate the exact spectral embedding X. Moreover, let C * (resp. C * ) be the optimal set of k centroids minimizing the k-means cost on X (resp. X). We will see that several (not all) approximation methods of this Section achieve an end-to-end approximation guarantee of the form

f (C * ; X) 1 /2 -f ( C * ; X) 1 /2 ≤ ε,
for some small ε with -at least-constant probability. Such an end-to-end guarantee is indeed more desirable than a simple guarantee on the distance between U k and Ũk : it informs us on the approximation quality of the attained clustering.

Nyström-based methods

The Nyström-based methods discussed in Sec. The first step towards an end-to-end guarantee relies on the following result:

Lemma 1 (see the proof of Theorem 6 in [START_REF] Boutsidis | Spectral clustering via the power methodprovably[END_REF]). Denote by C * the optimal centroid set obtained by solving k-means on the rows of Ũk . It holds that f (C * ; X)

1 /2 -f ( C * ; X) 1 /2 ≤ 2 E F , (16)
where

E = U k U k -Ũk Ũ k .
This means that the error made by considering the optimal k-means solution based on Ũk (instead of U k ) is controlled by the Frobenius norm of the projector difference

E = U k U k -Ũk Ũ k . Furthermore, since 14 E F ≤ √ 2k E 2 and E 2 = sin(Θ (U k , Ũk )) 2 ,
we can apply the Davis-Kahan sinΘ perturbation theorem (see for instance Section VII of [START_REF] Bhatia | Matrix analysis[END_REF]) and, provided that σ k -σk+1 > 0, obtain:

E F ≤ √ 2k E 2 ≤ √ 2k A -Ã 2 σ k -σk+1 ,
where {σ i } (resp. { σi }) are the singular values of A (resp. Ã) ordered decreasingly 15 . The final bound is obtained by combining the above with the leverage score sampling bound given by Theorem 2:

Theorem 4. Let Ũk be the eigenvectors obtained by running Alg. 3 on A = 2I -L n (with the leverage score sampling scheme for the m samples S of step 1). Denote by C * the optimal centroid set obtained by solving k-means on the rows of Ũk . Then, for some constant C > 1, we have

f (C * ; X) 1 /2 -f ( C * ; X) 1 /2 ≤ 2 √ 2k σ k -σk+1 σ k+1 (A) + Ck log (k/δ ) m n ∑ j=k+1 σ j
with probability at least 0.8 -2δ .

Examining the above bound one notices that 2 √ 2k

σ k+1 (A) σ k -σk+1
is independent of the number of samples. The incompressibility of this error term emanates from A being (in general) different from its best low-rank approximation. On the other hand, all remaining error terms can be made independent of k and n by setting

m = O k √ k log k n ∑ j=k+1 σ j σ k -σk+1 .
This end-to-end guarantee is not satisfactory for several reasons. First of all, it relies on the assumption σ k > σk+1 , which is not necessarily true. Moreover, the Davis-Kahan theorem could in theory guarantee E 2 ≤ A k -Ãk 2 /σ k and E 2 ≤ A -Ãk 2 /σ k , which are stronger than the bound depending on A -Ã 2 that we used. Unfortunately, Nyström approximation theorems do not give controls on A k -Ãk 2 nor on A -Ãk 2 , impeding tighter end-to-end bounds.

Graph coarsening

Inspired by the algebraic multi-grid, researchers realized early on that a natural way to accelerate spectral clustering is by graph coarsening [START_REF] Hendrickson | A multi-level algorithm for partitioning graphs[END_REF][START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF][START_REF] Dhillon | Weighted graph cuts without eigenvectors a multilevel approach[END_REF]. Here, instead of solving the clustering problem directly on G, one may first reduce it to a coarser graph G c consisting of m n nodes using a multi-level graph coarsening procedure. The expensive eigen-decomposition computation is done at a lower cost on the representative matrix of the small graph and the final spectral embedding is obtained by inexpensively lifting and refining the result.

In the notation of [START_REF] Loukas | Spectrally approximating large graphs with smaller graphs[END_REF], coarsening involves a sequence of c + 1 graphs

G = G 0 = (V 0 , E 0 , W 0 ) G 1 = (V 1 , E 1 , W 1 ) • • • G c = (V c , E c , W c ) (17)
of decreasing size

n = n 0 > n 1 > • • • > n c = m,
where each vertex of G represents one of more vertices of G -1 . To express coarsening in algebraic form, we suppose that L(G 0 ) = L is the combinatorial Laplacian associated with G. We then obtain L(G c ) by applying the following repeatedly

L(G ) = P ∓ L(G -1 )P + , (18)
where P ∈ R n ×n -1 is a matrix with more columns than rows, = 1, 2, . . . , c is the level of the reduction and symbol ∓ denotes the transposed pseudoinverse. An eigenvector ũ ∈ R m of L(G c ) is lifted back to R n by backwards recursion

ũ -1 = P ũ ,
where ũc = ũ. Matrices P 1 , P 2 , . . . , P c are determined by the transformation performed at each level. Specifically, one should define for each level a surjective map ϕ : V -1 → V between the original vertex set V -1 and the smaller vertex set V . We refer to the set of vertices V (r) -1 ⊆ V -1 mapped onto the same vertex v r of V as a contraction set:

V (r) -1 = {v ∈ V -1 : ϕ (v) = v r }
It is easy to deduce from the above that contraction sets induce a partitioning of V -1 into n subgraphs, each corresponding to a single vertex of V .

Then, for any v r ∈ V and v i ∈ V -1 , matrices P ∈ R n ×n -1 and P + ∈ R n -1 ×n are given by:

P (r, i) =    1 |V (r) -1 | if v i ∈ V (r) -1 0 otherwise and P + (i, r) = 1 if v i ∈ V (r) -1 0 otherwise.
The preceding construction is the only one that guarantees that every L(G ) will be the combinatorial Laplacian associated with G [START_REF] Loukas | Graph reduction with spectral and cut guarantees[END_REF].

Note that from a computational perspective the reduction is very efficient and can be carried out in linear time: each coarsening level entails multiplication by a sparse matrix, meaning that O(e) and O(n) operations suffice, respectively, to coarsen L and lift any vector (such as the eigenvectors of L(G c )) from R m back to R n .

Coarsening for spectral clustering

Using coarsening effectively boils down to determining for each how to partition G -1 into n contraction sets V (1) , . . . ,V (n ) , such that, after lifting, the first k eigenvectors Ũk of L(G c ) approximate the spectral embedding U k derived from L. Alternatively, one may also solve the k-means problem in the small dimension and only lift the resulting cluster assignments [START_REF] Dhillon | Weighted graph cuts without eigenvectors a multilevel approach[END_REF]. This scheme is computationally superior but we will not discuss it here as it does not come with any guarantees. Perhaps the most simple (and common) method of forming contraction sets is by the heavy edge matching heuristic-originally developed in the multi-grid literature and first considered for graph partitioning in [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF]. This method is derived based on the intuition that, the larger the weight of an edge, the less likely it will be that the vertices it connects will reside in different clusters. We should therefore aim to contract pairs of vertices connected by a heavy edge (i.e., of large weight) first. Let us consider this case further. By focusing on edges, we basically constrain ourselves by enforcing that every contraction set V (r) contains either two nodes connected by an edge, or a single node, signifying that said node is chosen to remain as is in the coarser graph. As such, we can reformulate the problem of selecting contraction sets at each level as that of selecting the largest number of edges (to attain the largest reduction), while also striving to make the cumulative sum of selected edge weights as large as possible (giving preference to heavy edges). This is exactly the maximum weight matching problem, which can be approximated in linear time [START_REF] Duan | Linear-time approximation for maximum weight matching[END_REF].

A plethora of numerical evidence motivates the use of matching-based coarsening methods, such as the heavy-edge heuristic, for accelerating spectral clustering [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF][START_REF] Dhillon | Weighted graph cuts without eigenvectors a multilevel approach[END_REF]115]. From a theoretical perspective, the approximation quality of matching-based methods was characterized in [START_REF] Loukas | Spectrally approximating large graphs with smaller graphs[END_REF]. Therein, the matching was constructed in the following randomized manner: Algorithm 5. Randomized edge contraction (one level) [START_REF] Loukas | Spectrally approximating large graphs with smaller graphs[END_REF] Input. A graph G = (V, E) 1. Associate with each e i j ∈ E a probability p i j > 0. 2. While |E| > 0:

a. Draw a sample e i j from E with probability ∝ p i j . b. Remove from E both e i j as well as all edges sharing a common endpoint with it. c. Construct contraction set (v i , v j ).

Output: Contraction sets

The following approximation result is known: Theorem 5 (Corollary 5.1 in [START_REF] Loukas | Spectrally approximating large graphs with smaller graphs[END_REF]). Consider a graph with bounded degrees d i n and λ k ≤ min e i j ∈E

d i +d j 2
. Suppose that the graph is coarsened by Algorithm 5, using a heavy-edge potential such that p i j ∝ w i j . For sufficiently large n, a single level, and δ > 0,

f (C * ; X) 1 /2 -f ( C * ; X) 1 /2 = O   1 -m n δ ∑ k =2 λ λ k+1 -λ k  
with probability at least 1-δ . Above, C * is the optimal k-means solution when using the lifted eigenvectors of L c as a spectral embedding.

We deduce that coarsening works better when the spectral clustering problem is easy (as quantified by the weighted gap ∑ k =2 λ /(λ k+1λ k )) and the achieved error is linear on the reduction ratio 1m/n.

There also exist more advanced techniques for selecting contraction sets that come with stronger guarantees w.r.t. the attained reduction and quality of approximation, but feature running time that is not smaller than that of spectral clustering [START_REF] Loukas | Graph reduction with spectral and cut guarantees[END_REF]. In particular, these work also with the normalized Laplacian and can be used to achieve multi-level reduction. Roughly, their strategy is to identify and contract sets S ⊂ V for which x(i) ≈ x( j) for all vectors x ∈ U k and v i , v j ∈ S. This strategy ensures that the best partitionings of G are preserved by coarsening. We will not expand on these methods here as they do not aim to improve the running time of spectral clustering.

Other approaches

In the following, we present two additional approaches for approximately computing spectral embeddings. The former can be interpreted as a sampling-based method (but in a different manner than the techniques discussed so far), whereas the latter is only vaguely linked to sampling. Nevertheless, we find that both techniques are very interesting and merit a brief discussion.

Spectral sparsification

This approach is best suited for cases when the input of spectral clustering is directly a graph 16 . Differently from the methods discussed earlier, here the aim is to identify a Laplacian matrix L of the same size as L but with fewer entries. Additionally, it should be ensured that

1 1 + ε x Lx ≤ x Lx ≤ (1 + ε)x Lx for all x ∈ R n (19)
for some small constant ε > 0 There are different ways to use sparsification in order to accelerate spectral clustering. The most direct one is to exploit the fact that the eigenvalues λk and eigenvectors Ũk of L approximate, respectively, the eigenvalues and eigenvectors of L up to multiplicative error. This yields the same flavor of guarantees as in graph coarsening and ensures that the computational complexity of the partial eigen-decomposition will decrease when e = ω(n log n). A variation of this idea was considered in [142], though the latter did not provide a complete error and complexity analysis. Alternative approaches are also possible. We refer the interested reader to [136] for a rigorous argument that invokes a Laplacian solver.

Despite these exciting developments, we should mention that the overwhelming majority of graph sparsification algorithms remain in the realm of theory. That is, we are currently not aware of any practical and competitive implementation and thus retain a measure of skepticism with regards to their utility in the setting of spectral clustering.

Random eigenspace projection

There also exists approaches that do not explicitly rely on sampling. The key starting point here is that, with regards to spectral clustering, one does not need the eigenvectors exactly-any rotation of U k suffices (indeed, k-means is an algorithm based on distances and rotations conserve distances). Even more generally, consider Ũk ∈ R n×m with m ≥ k and denote:

ε = min Q∈Q U k I k×m Q -Ũk F ,
where Q is the space of m × m unitary matrices and I k×m consists of the first k rows of an m × m identity matrix.

The following lemma (which is a generalization of Lemma 1) shows how ε can be used to provide control on the k-means error: Lemma 2 (Lemma 3.1 in [START_REF] Martin | Fast approximate spectral clustering for dynamic networks[END_REF]). Let C * be the optimal solution of the k-means problem on Ũk . It holds that 17

f (C * ; X) 1 /2 -f ( C * ; X) 1 /2 ≤ 2ε. ( 20
)
There exists (at least) two approaches to efficiently compute Ũk while controlling ε [START_REF] Boutsidis | Spectral clustering via the power methodprovably[END_REF]133] (see also related work in [START_REF] Halko | Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions[END_REF]). We will consider here a simple variant of the one proposed in [133] and further analyzed in [START_REF] Martin | Fast approximate spectral clustering for dynamic networks[END_REF]: Let G ∈ R n×m be a random Gaussian matrix with centered i.i.d. entries, each having variance 1 m . Further, suppose that we project G onto span(U k ) by multiplying each one of its columns by an ideal projector P k defined as

P k = U I k 0 0 0 U . (21)
Theorem 6 ([133, 95]). Let C * be the optimal solution of the k-means problem on the rows of Ũk = P k G. For every δ ≥ 0, one has f (C * ; X)

1 /2 -f ( C * ; X) 1 /2 ≤ 2 k m ( √ k + δ ), (22)
with probability at least 1exp(-δ 2 /2).

This result means that for an ideal projector P k , dimension m = O(k 2 ) suffices to guarantee good approximation (since the error becomes independent of k and n)! A similar argument also holds when the entries of G, instead of being Gaussian, 17 A remark on the definition of the k-means cost. Note that, here, the lines X of Ũk are points in dimension m ≥ k, such that the optimal centroid set C * minimizing the k-means cost on X is a set of k points in dimension m ≥ k. In this context, the notation f ( C * ; X) is ill-defined: it is a sum of distances between points that do not necessarily have the same dimension. We abuse notations and give the following meaning to f ( C; X). First, consider the matrix form of the k-means cost, as used in the proofs of Lemmas 1 and 2:

f (C; X) = X -CC X 2 F , where X = (x 1 | . . . |x n ) ∈ R n×k and C ∈ R n×k
is the (weighted) cluster indicator matrix associated to the Voronoi tesselation of X given C: C i = 1/ √ s if data point i belongs to cluster , and 0 otherwise; where s is the size of cluster . Now, let C ∈ R n×k be the cluster indicator matrix associated to the Voronoi tesselation of X given C. One writes:

f ( C; X) = X -C C X 2 F .
are selected i.i.d. from {-√ 3, 0, + √ 3} with probabilities {1/6, 2/3, 1/6}, respectively [START_REF] Achlioptas | Database-friendly random projections: Johnson-lindenstrauss with binary coins[END_REF]. This construction has the benefit of being sparser and, moreover, is reminiscent of sampling. It should be noted that in [133], m = O(log n) was deemed enough because one only wanted that the distance between two rows of U k was approximated by the distance between the same two rows of Ũk . There was in fact no end-to-end control on the k-means error.

The discussion so far assumed that P k is an ideal projector onto span(U k ). However, in practice one does not have access to this projector as we are in fact in the process of computing U k . One may choose to approximate the action of P k by an application of a matrix function h on the representative matrix R [132, 111]. Assuming a point λ * in the interval [λ k , λ k+1 ) is known, one may select a polynomial [123] or rational function [START_REF] Isufi | Autoregressive moving average graph filtering[END_REF][START_REF] Loukas | Distributed autoregressive moving average graph filters[END_REF] that approximates the ideal low-pass response, i.e., h(λ ) = 1 if λ ≤ λ * and h(λ ) = 0, otherwise. The approximated projector Pk = h(R) can be designed to be very close to P k . For instance, in the case of Chebychev polynomials of order c using the arguments of [80, Lemma 1] it is easy to prove that w.h.p. using h(R) instead of P k does not add more than O(c -c √ n) error in [START_REF] Boutsidis | Unsupervised Feature Selection for the kmeans Clustering Problem[END_REF]. Furthermore, the operation Pk G can conveniently be computed in O(mce) time via this polynomial approximation.

The last ingredient needed for this approximation is λ * , i.e., a point in the interval [λ k , λ k+1 ). Finding efficiently a valid λ * is difficult. An option is to rely on eigencount techniques [START_REF] Di Napoli | Efficient estimation of eigenvalue counts in an interval[END_REF]105,109] to find one in 18 O(ck 2 (log n)(e + n log(λ n /(λ k+1λ k )))) time, which features similar complexity as the Lanczos method (see discussion in Sec. 2.3). Another option is to content oneself with values of λ * known only to be close to the interval [λ k , λ k+1 ), but thereby loosing the end-to-end guarantee [133].

Sampling in the spectral feature space

Having computed (or approximated) the spectral embedding X = (x 1 , x 2 , . . . , x n ), what remains is to solve the k-means problem on X, in order to obtain k centroids together with the associated k classes obtained after Voronoi tessellation.

The usual heuristic used to solve the k-means problem, namely the Lloyd-Max algorithm, is already very efficient as it runs in O(nk 2 t) time as seen in Sec. 2.3. Nonetheless, this section considers ways to accelerate k-means even further. In the 18 Proof sketch: Given λ ∈ (0, λ n ], denote by j the largest integer such that λ j ≤ λ and by P j the orthogonal projector on U j . Let G ∈ R n×m be a random Gaussian matrix with centered i.i.d. entries, each having variance 1 m and denote by ĵ = P j G F . Relying on Theorem 4.1 (and the following discussion in Section 4.2) of [109] with E λ = 0, one has with prob. at least 1ε that

(1 -δ ) j ≤ ĵ ≤ (1 + δ ) j for all j = 1, . . . , n provided m ≥ 1 δ 2 log n ε . Setting δ = 1/(2k + 3), gives w.h.p. that 2k+2 2k+3 j ≤ ĵ ≤ 2k+4 2k+3 j for all j = 1, . . . , n provided m ≥ O(k 2 log n).
This implies that w.h.p. for every j ≤ k + 1 it must be that round( ĵ) = j, whereas when j > k + 1 we have round( ĵ) > k + 1. Note that round( ĵ) is the closest integer to ĵ. By dichotomy on λ ∈ (0, λ n ], one thus finds a λ * in time O(ck 2 (log n)(e + n log(λ n /(λ k+1λ k )))).

following, we classify the relevant literature in five categories and point towards representative references for each case. In our effort to provide depth (as well as breadth) of presentation, the rest of the section details only methods that belong to the first and last categories.

• Exact acceleration of Lloyd-Max. There exists exact accelerated Lloyd-Max algorithms, some of them based on avoiding unnecessary distance calculations using the triangular inequality [START_REF] Hamerly | Accelerating Lloyd's Algorithm for k-Means Clustering[END_REF]101], or on optimized data organization [START_REF] Kanungo | An Efficient k-Means Clustering Algorithm: Analysis and Implementation[END_REF], and others concentrating on clever initializations [START_REF] Arthur | k-means++: The advantages of careful seeding[END_REF]104]. The latter concern sampling and are discussed in Sec. : it is based on mini-batches sampled uniformly at random from X. We will not discuss further this method as it does not come with guarantees on the cost of the obtained solution. • Methods involving sampling in the Fourier domain. There are a few samplingbased heuristics to solve the k-means problem, that are different from the Lloyd-Max algorithm. For instance, the work in [START_REF] Keriven | Compressive K-means[END_REF] proposes to sample in the frequency domain to obtain a sketch from which one may recover the centroids with an orthogonal matching pursuit algorithm specifically tailored to this kind of compressive learning task [START_REF] Gribonval | Compressive Statistical Learning with Random Feature Moments[END_REF]. These methods are reminiscent of the random Fourier features sketching approach introduced in Sec. 3.2. We will not discuss them further. • Methods involving sampling features. Similarly to ideas presented in Sec. 3.4.2 but here specific to the k-means setting, some works reduce the ambient dimension of the vectors, either by selecting a limited number of features [START_REF] Boutsidis | Unsupervised Feature Selection for the kmeans Clustering Problem[END_REF][START_REF] Altschuler | Greedy Column Subset Selection: New Bounds and Distributed Algorithms[END_REF], or by embedding all points in a lower dimension using random projections [START_REF] Boutsidis | Randomized Dimensionality Reduction for $k$ -Means Clustering[END_REF][START_REF] Cohen | Dimensionality reduction for k-means clustering and low rank approximation[END_REF][START_REF] Makarychev | Performance of Johnson-Lindenstrauss Transform for k-Means and k-Medians Clustering[END_REF].

The tightest results to day are a (1 + ε) multiplicative error on the k-means cost f either by randomly selecting O(ε -2 k log k) features or by projecting them on a random space of dimension O(ε -2 log (k/ε)) (sublinear in k!). The sampling result is useless in the spectral clustering setting as the ambient dimension of the spectral features is already k. The projection result could in principle be applied in our setting, to reduce the cost of the k-means step to O(tnk log k). We will nevertheless not discuss it further in this chapter. • Methods involving sampling points. Finally, the last group of existing methods are the ones that solve k-means on a subset S of X, before lifting back the result on the whole dataset. We classify such methods in two categories. In Sec. 5.2, we detail methods that are graph-agnostic, meaning that they apply to any k-means problem; and in Sec. 5.3 we discuss methods that explicitly rely on the fact that the features x were in fact obtained from a known graph. We argue that the latter are better suited to the spectral clustering problem.

Clever initialization of the Lloyd-Max algorithm

Recall that the k-means objective on X is to find the k centroids C = (c 1 , . . . , c k ) that minimize the following cost function:

f (C; X) = ∑ x∈X min c∈C x -c 2 2 . ( 23
)
and that C * = arg min C f (C; X) is the optimal solution attaining cost f * = f (C * ; X).

Recall also that the Lloyd-Max algorithm (see Algorithm 2) converges to a local minimum of f , that we will denote by C lm , for which the cost function equals f lm = f (C lm ; X). It is crucial to note that the initialization of centroids C ini in the first step of the Lloyd-Max algorithm, which usually is done by randomly selecting k points in X, is what determines the distance | f *f lm | to the optimal. As such, significant efforts have been devoted to smartly selecting C ini by various sampling schemes.

As usual, we also face here the usual trade-off between sampling effectively and efficiently. The fastest sampling method is of course uniformly at random, but it does not come with any guarantee on the quality of the local minimum C lm it leads to. An alternative sampling scheme, called k-means++ initialization, is based on the following more general D 2 -sampling algorithm. k-means++ initialization boils down to running Alg. 6 with m = k to obtain a set of k initial centroids. Importantly, when the Lloyd-Max heuristic is run with this initialization, the following guarantee holds: Theorem 7 ( [START_REF] Arthur | k-means++: The advantages of careful seeding[END_REF]). For any set of data points, the cost f lm obtained after Lloyd-Max initialized with k-means++ is controlled in expectation:

E( f lm ) ≤ 8(log k + 2) f * .
In terms of computation cost, D 2 -sampling with m = k runs in O(nkd), that is, O(nk 2 ) in our setting of a spectral embedding X in dimension k. This work inspired other initialization techniques that come with similar guarantees and are in some cases faster [START_REF] Bahmani | Scalable K-means++[END_REF][START_REF] Bachem | Fast and Provably Good Seedings for k-Means[END_REF]. The interested reader is referred to the review [START_REF] Celebi | A comparative study of efficient initialization methods for the k-means clustering algorithm[END_REF] for further analyses on the initialization of k-means.

Graph agnostic sampling methods: coresets

The rest of Sec. 5, considers sampling methods that fall in the following framework: (i) sample a subset S of X, (ii) solve k-means on S, (iii) lift the result back on the whole dataset X. Sec. 5.2, focuses on coresets: general sampling methods designed for any arbitrary k-means problem; whereas in Sec. 5.3, we will take into account the specific nature of the spectral features encountered in spectral clustering algorithms.

Definition

Let S ⊂ X be a subset of X of size m. To each element s ∈ S associate a weight ω(s) ∈ R + . Define the estimated k-means cost associated to the weighted set S as:

f (C; S) = ∑ s∈S ω(s) min c∈C s -c 2 2 . ( 24
)
Definition 2 (Coreset). Let ε ∈ (0, 1 2 ). The weighted subset S is a ε-coreset for f on X if, for every set C, the estimated cost is equal to the exact cost up to a relative error:

∀C f (C; S) f (C; X) -1 ≤ ε. (25)
This is the so-called "strong" coreset definition 19 , as the ε-approximation is required for all C. The great interest of finding a coreset S comes from the following fact. Writing C * the set minimizing f , the following inequalities hold

(1 -ε) f (C * ; X) ≤ (1 -ε) f ( C * ; X) ≤ f ( C * ; S) ≤ f (C * ; S) ≤ (1 + ε) f (C * ; X).
The first inequality comes from the fact that C * is optimal for f , the second and last inequality are justified by the coreset property of S, and the third inequality comes from the optimality of C * for f . This has two consequences:

1. First of all, since ε < 1 2 :

f (C * ; X) ≤ f ( C * ; X) ≤ (1 + 4ε) f (C * ; X),
meaning that C * is a well controlled approximation of C * with a multiplicative error on the cost. 2. Estimating C * can be done using the Lloyd-Max algorithm on the weighted subset 20 S, thus reducing the computation time from O(nk 2 ) to O(mk 2 ).

Coreset methods for k-means thus follow the general procedure: Coreset methods compete with one another on essentially two levels: the coreset size m should be as small as possible in order to decrease the time of Lloyd-Max on S, and the coreset itself should be sampled efficiently (at least faster than running k-means on the whole dataset!), which turns out in fact to be a strong requirement. The reader interested in an overview of coreset construction techniques is referred to the recent review [START_REF] Munteanu | Coresets-Methods and History: A Theoreticians Design Pattern for Approximation and Streaming Algorithms[END_REF], as well as chapter 2 of this book.

An instance of coreset-sampling algorithm

We focus on a particular coreset algorithm proposed in [START_REF] Bachem | Practical Coreset Constructions for Machine Learning[END_REF] that builds upon results developed in [START_REF] Langberg | Universal ε-approximators for integrals[END_REF][START_REF] Feldman | A unified framework for approximating and clustering data[END_REF]: it is not state-of-the-art in terms of coreset size, but has the advantage of being easy to implement and fast enough to compute. It reads: Algorithm 8: a coreset sampling algorithm [START_REF] Bachem | Practical Coreset Constructions for Machine Learning[END_REF]. Input. X, m the number of required samples, t an iteration number 1. Repeat t times: draw a set of size k using D 2 -sampling. Out of the t sets obtained, keep the set B that minimizes f (B; X).

2. α ← 16(log k + 2) 3. For each b ∈ B, define B the set of points in X in the Voronoi cell of b 4. Set φ = 1 n f (B; X). 5. For each b ∈ B and each x ∈ B , define s(x) = α φ x -b 2 2 + 2α φ |B | ∑ x ∈B x -b 2 2 + 4n |B | 6.
Define the probability of sampling x i as p i = s(x i )/ ∑ x s(x) 7. S ← sample m nodes i.i.d. with replacement from p and associate to each sample s the weight ω s = 1 mp s .

Output: A weighted set S of size m.

Theorem 2.5 of [START_REF] Bachem | Practical Coreset Constructions for Machine Learning[END_REF] states:

Theorem 8. Let ε ∈ (0, 1/4) and δ ∈ (0, 1). Let S be the output of Alg. 8 with t = O(log 1/δ ). Then, with probability at least 1δ , S is a ε-coreset provided that:

m = Ω k 4 log k + k 2 log 1/δ ε 2 . ( 26
)
The computation cost of running this coreset sampling algorithm, running Lloyd-Max on the weighted coreset, and lifting the result back to X is dominated, when 21 n k, by step 1 of Alg. 6 and thus sums up to O(nk 2 log 1/δ ).

Remark 1. The coreset sampling strategy underlying this algorithm relies on the concept of sensitivity [START_REF] Langberg | Universal ε-approximators for integrals[END_REF]. Many other constructions of coresets for k-means are possible [START_REF] Munteanu | Coresets-Methods and History: A Theoreticians Design Pattern for Approximation and Streaming Algorithms[END_REF] with better theoretical bounds then [START_REF] Calvetti | An implicitly restarted lanczos method for large symmetric eigenvalue problems[END_REF]. Nevertheless, as the coreset line of research has been essentially theoretical, practical implementations of coreset-sampling algorithms are scarce. A notable exception is for instance the work in [START_REF] Frahling | A fast k-means implementations using coresets[END_REF] that proposes a scalable hybrid coreset-inspired algorithm for k-means. Other exceptions are the sampling algorithms based on the farthest-first procedure, a variant of D 2 -sampling that chooses each new sample to be arg max i d i instead of drawing it according to a probability proportional to d i . Once S of size m is drawn, then ∀s ∈ S, each weight ω s is set to be the cardinal of the Voronoi cell associated to s. Authors in [113] show that such weighted sets computed by different variants of the farthest-first algorithm are ε-coresets, but for values of ε that can be very large. For a fixed ε, the number of samples necessary to have a ε-coreset with this type of algorithm is unknown (see also chapter 3 of this book).

Graph-based sampling methods

The methods discussed so far in this section are graph agnostic both for the sampling procedure and the lifting: they do not take into account that, in spectral clustering, X are in fact spectral features of a known graph.

A recent line of work [133, [START_REF] Martin | Fast approximate spectral clustering for dynamic networks[END_REF][START_REF] Gadde | Active Semi-supervised Learning Using Sampling Theory for Graph Signals[END_REF][START_REF] Gadde | Active learning for community detection in stochastic block models[END_REF] based on Graph Signal Processing (GSP) [122,118] leverages this additional knowledge for accelerating both the sampling and the lifting steps. For the purpose of the following discussion, define by z ∈ R n the ground truth indicator vector of cluster , i.e., z (i) = 1 if node v i is in cluster , and 0 otherwise. The goal of spectral clustering is, of course, to recover {z } =1,...,k . 21 To be precise, the statement holds if n ≥ O k In both cases, P S should be set to 1 N I m if uniform sampling was chosen, to diag(p * s 1 , . . . p * s m ) if leverage score sampling was chosen, and to diag(π s 1 , . . . π s m ) if DPP sampling was chosen. 4. Assign each node j to the cluster for which ẑ ( j)/ ẑ 2 is maximal.

Output: A partition of X in k clusters

To aid understanding, let us start by a high-level description of Algorithm 9. The indicator vectors z are interpreted as graph signals that are (approximately) bandlimited on the similarity graph G (see Sec. 5.3.1 for a precise definition). As such, there is no need to measure these indicator vectors everywhere: one can take advantage of generalized Shannon-type sampling theorems to select the set S of m nodes to measure (step 1). Then k-means is performed on S to obtain the indicator vectors z r ∈ R m on the sample set S (step 2). These reduced indicator vectors are interpreted as noisy measurements of the global cluster indicator vectors z on S. The solutions z r are lifted back to X as ẑ via solving an inverse problem taking into account the bandlimitedness assumption or via label-propagation on the graph structure reminiscent of semi-supervised learning techniques (step 3). As the lifted solutions ẑ do not have a binary structure as true indicator vectors should have, an additional assignment step is necessary: assign each node j to the class for which ẑ ( j) ẑ 2 is maximal (step 4). The rest of this section is devoted to the discussion of the three sampling schemes as well as the two lifting procedures considered in this framework. To this end, we will first introduce a few graph signal processing (GSP) concepts in Sec. 5.3.1 before discussing in Sec. 5.3.2 several examples of graph sampling theorems appropriate to the spectral clustering context.

A brief introduction to graph signal processing (GSP)

Denote by U = (u 1 | . . . |u n ) ∈ R n×n the matrix of orthonormal eigenvectors of the Laplacian matrix L, with the columns ordered according to their associated sorted eigenvalues: 0 = λ 1 ≤ λ 2 ≤ . . . ≤ λ n . In the GSP literature [122, 118], these eigenvectors are interpreted as graph Fourier modes for two main reasons:

• By analogy to the ring graph, whose Laplacian matrix is exactly the (symmetric) double derivative discrete operator, and is thus diagonal in the basis formed by the classical 1D discrete Fourier modes. • A variational argument stemming from the Dirichlet form can be exploited to express eigenvectors u i of L as the basis of minimal variation x Lx =

1 2 ∑ i j W i j [x(i) -x( j)]
2 on G and eigenvalues λ i as a sum of local variations of u i , i.e., a generalized graph frequency.

A graph signal z ∈ R n is a signal that is defined on the nodes of a graph: its i-th element is associated to node v i . Given the previous discussion, the graph Fourier transform of z, denoted by z, is its projection on the graph Fourier modes: z = U z ∈ R n . The notion of graph filtering naturally follows as a multiplication in the Fourier domain. More precisely, define a real-valued filter function h(λ ) defined on [0, λ n ]. The signal x filtered by h reads Uh(Λ )U x, where we use the convention

h(Λ ) = diag(h(λ 1 ), h(λ 2 ), . . . , h(λ n )).
In the following, we will use the following notation for graph filter operators:

h(L) = Uh(Λ )U . (27)
For more details on the graph Fourier transform and filtering, their various definitions and interpretations, we refer the reader to [131].

Of interest for the discussion in this chapter, one may define bandlimited graph signals as linear combinations of the first few low-frequency Fourier modes. Writing U k = (u 1 | . . . |u k ) ∈ R n×k , we have the formal definition:

Definition 3 (k-bandlimited graph signal). A graph signal z ∈ R n is k-bandlimited if z ∈ span(U k ), i.e., ∃ α ∈ R k such that z = U k α.
To grasp why the notion of k-bandlimitedness lends itself natually to the approximation of spectral clustering, consider momentarily a graph with k disconnected components and z ∈ R n the indicator vector of cluster . It is a well known property of the (combinatorial) Laplacian that {z } =1,...,k form a set of orthogonal eigenvectors of L associated to eigenvalue 0: that is, the set of indicator vectors {z } =1,...,k form a basis of span(U k ). Understanding arbitrary graphs with block structure as a perturbation of the ideal disconnected component case, the indicator vectors {z } =1,...,k of the blocks should live close to span(U k ) (in the sense that the difference between any z and its orthogonal projection onto span(U k ) is small). This in turn implies that every z should be approximately k-bandlimited.

As we will see next, the bandlimitedness assumption is very useful because it enables us to make use of generalized versions of Nyquist-Shannon sampling theorems, taking into account the graph.

Graph sampling theorems

The periodic sampling paradigm of the Shannon theorem for classical bandlimited signals does not apply to graphs without specific regular structure. In fact, a number of sampling schemes have been recently developed with the purpose of generalizing sampling theorems to graph signals [117,[START_REF] Chen | Discrete Signal Processing on Graphs: Sampling Theory[END_REF]134,109] (see [START_REF] Lorenzo | Chapter 9 -sampling and recovery of graph signals[END_REF] for a review of existing schemes).

Let us introduce some notations. Sampling entails selecting a set S = (s 1 , . . . , s m ) of m nodes of the graph. To each possible sampling set, we associate a measurement matrix M = (δ s 1 |δ s 2 | . . . |δ s m ) ∈ R m×n where δ s i ( j) = 1 if j = s i , and 0 otherwise. Now, consider a k-bandlimited signal z ∈ span(U k ). The measurement of z on S reads:

y = Mz + n ∈ R m , (28)
where n models measurement noise. The sampling question boils down to: how should we sample S such that one can recover any bandlimited z given its measurement y? There are three important components to this question: (i) how many samples m do we allow ourselves (m = k being the strict theoretical minimum)? (ii) how much does it cost to sample? (iii) how do we in practice recover z from y and how much does that inversion cost? There are a series of works that propose greedy algorithms to find the "best" set S of minimal size m = k that embed all k-bandlimited signals (see for instance [129] and references therein). These algorithms cost O(nk 4 ) and are thus not competitive in our setting 22 . Moreover, in our case, we don't really need to be that strict on the number of samples and can allow more than k samples. A better choice is to use random graph sampling techniques. In the following we consider two types of independent sampling (uniform and leverage-score sampling) as well as a more involved method based on determinantal point processes. Independent sampling. In the i.i.d. setting, one defines a discrete probability distribution p ∈ R n over the node set V. The sampling set S is then generated by drawing m nodes independently with replacement from p. At each draw, the probability to sample node v i is denoted by p i . We have ∑ i p i = 1 and write P = diag(p). Under this sampling scheme, the following Restricted Isometry Property holds for the associated measurement matrix M [109].

Theorem 9. For any δ , ε ∈ (0, 1), with probability at least 1δ :

(1 -ε) z 2 2 ≤ 1 m MP -1/2 z 2 2 ≤ (1 + ε) z 2 2 ( 29
)
for all z ∈ span(U k ) provided that m ≥ 3 ε 2 (ν k p ) 2 log 2k δ ( 30
)
where ν k p is the so-called graph weighted coherence:

ν k p = max i p -1/2 i U k δ i 2 . ( 31
)
This property is important as it says, in a nutshell, that any two different bandlimited signals will be identifiable post-sampling provided the number of samples is large enough. The concept of large enough depends on (ν k p ) 2 : a measure of the interplay between the probability distribution and the norms of the rows of U k . In the uniform i.i.d. case since

p i = 1/n, one has (ν k p ) 2 = n max i U k δ i 2 2
, which stays under control only for very regular graphs, but can be close to n in irregular graphs such as the star graph. The good news is that there exists an optimal sampling distribution (in the sense that it minimizes the right-hand side of inequality [START_REF] Chen | Large Scale Spectral Clustering with Landmark-Based Representation[END_REF]) that adapts to the graph at hand:

p * i = U k δ i 2 2 k (32)
In fact, in this case, (ν k p * ) 2 matches its lower bound k and the necessary number of samples m to embed all bandlimited signals drops to O(k log k). The distribution p * is also referred to by the name "leverage scores" in parts of the literature (see discussion in Sec. 3.1) [START_REF] Drineas | Lectures on randomized numerical linear algebra[END_REF]. As such, i.i.d. sampling under p * will be referred to as leverage score sampling. Now, for lifting, there are several options. Theorem 10. Let S be the i.i.d. nodes sampled with distribution p and M be the associated sampling matrix. Let ε, δ ∈ (0, 1) and suppose that m satisfies [START_REF] Chen | Large Scale Spectral Clustering with Landmark-Based Representation[END_REF]. With probability at least 1δ , for all z ∈ span(U k ) and n ∈ R m , the solution ẑ of (33) verifies:

ẑ -z 2 ≤ 2 m(1 -ε) P -1/2 S n 2 .
This means that a noiseless measurement of a k-bandlimited signal yields a perfect reconstruction. Also, this quantifies how increasing m reduces the error of reconstruction due to a noisy measurement. Note that this error may be large if there is a significant measurement noise on a node that has a low probability of being sampled. However, by definition, this is not likely to happen. • One can also use a label-propagation decoder reminiscent to semi-supervised learning techniques [START_REF] Bengio | Label Propagation and Quadratic Criterion[END_REF][START_REF] Chapelle | Semi-Supervised Learning[END_REF]: ẑ = arg min

w∈R n P -1/2 S (Mw -y) 2 2 + γ w g(L)w, (34)
where γ is a regularization parameter, g(L) a graph filter operator as in [START_REF] Celebi | A comparative study of efficient initialization methods for the k-means clustering algorithm[END_REF] with g(λ ) a non-decreasing function. As g is non-decreasing, the regularization term of (34) penalizes high frequency solutions, that is, solutions that are not smooth along paths of the graph. Theorems controlling the error of reconstruction are more involved and we refer the reader to Section 3.3 of [109] for details. • Other decoders [START_REF] Bellec | A sharp oracle inequality for graph-slope[END_REF]106] are in principle possible, replacing for instance the 2 Laplacian-based regularization w g(L)w by 1 -regularizers ∇w 1 , but they come with an increased computation cost, lesser guarantees, and have not been used for spectral clustering: we will thus not detail them further.

Let us discuss the computation costs of the previous sampling and lifting techniques. In terms of sampling time, uniform sampling is obviously the most efficient and runs in O(k). Leverage score sampling is dominated by the computation of the optimal sampling distribution p * of (32), which takes O(nk) time 23 . In terms of lifting time, solving the decoder of (33) costs O(nk + mk 2 ). Solving the decoder of (34) costs O(et) via the conjugate gradient method, where t is the iteration number of the gradient solver (usually around 10 or 20 iterations suffice to obtain good accuracy when g(L) = L).

This discussion calls for a few remarks. First of all, these theorems are valid if we suppose that z is exactly k-bandlimited, which is in fact only an approximation if we consider z to be the ground truth indicator vectors of the k clusters to detect in the spectral clustering context. In this case, we can always decompose z as the sum of its orthogonal projection onto span(U k ) and its complement β : z = U k U k z + β . (28) becomes y = MU k U k z + n where n now represents the sum of a measurement noise and the distance-to-model term Mβ . The aforementioned theorems can then be applied to U k U k z. Moreover, note that the decoder of (34) is not only faster than the other ones in general, it also does not constrain the solution ẑ to be exactly in span(U k ), which is in fact desirable in the spectral clustering context: we thus advocate for the decoder of [START_REF] Dall'amico | Optimized Deformed Laplacian for Spectrumbased Community Detection in Sparse Heterogeneous Graphs[END_REF].

DPP sampling. Determinantal Point Processes are a class of correlated random sampling strategies that strive to increase "diversity" in the samples, based on a kernel K expliciting the similarity between variables. DPP sampling has been used successfully in a number of applications in machine learning (see for instance [START_REF] Kulesza | Determinantal Point Processes for Machine Learning[END_REF]).

Denote by [n] the set of all subsets of {1, 2, . . . , n}. An element of [n] could be the empty set, all elements of {1, 2, . . . , n} or anything in between. DPPs are defined as follows:

Definition 4 (Determinantal Point Process [START_REF] Kulesza | Determinantal Point Processes for Machine Learning[END_REF]). Consider a point process, i.e., a process that randomly draws an element S ∈

[n]. It is determinantal if, ∀ A ⊆ S, P(A ⊆ S) = det(K A ),
where K ∈ R n×n , a semi-definite positive matrix 0 K 1, is called the marginal kernel; and K A is the restriction of K to the rows and columns indexed by the elements of A.

The marginal probability π i of sampling an element i is thus K ii . Consider the following projective kernel:

K k = U k U k . ( 35
)
One can show that DPP samples from such projective kernels are necessarily of size k. After measuring the k-bandlimited signal z on a DPP sample S, one has the choice between the same decoders as before (see Eqs. ( 33) and ( 34)). For instance:

Theorem 11. For all z ∈ span(U k ), let y = Mz + n ∈ R k be a noisy measurement of z on a DPP sample obtained from kernel K k . The decoder of (33) with P = diag(π 1 , . . . , π n ) necessarily enables perfect reconstruction up to the noise level. Indeed, one obtains:

ẑ -z 2 ≤ 1 λ min U k M P -1 S MU k P -1/2 S n 2 . ( 36
)
Proof. The proof is only partly in [129] and we complete it here. Let us write z

= U k α. Solving (33) entails computing α ∈ R k s.t. P -1/2 S (MU k α -y) 2
2 is minimal. Setting the derivative w.r.t. α to 0, and replacing y by MU k α + n, yields:

U k M P -1 S MU k α = U k M P -1 S MU k α + U k M P -1 S n.
Recall that S is a sample from a DPP with kernel K k : det(MU k U k M ) is thus strictly superior to 0, which implies that MU k is invertible, which in turn implies that α = α + (MU k ) -1 n. One thus has ẑ -

z 2 = α -α 2 = (MU k ) -1 n 2 = (P -1/2 S MU k ) -1 P -1/2 S n 2
. Using the matrix 2-norm to bound this error yields

ẑ -z 2 ≤ λ max U k M P -1 S MU k -1 P -1/2 S n 2 ,
as claimed.

Several comments are in order:

• The particular choice of kernel

K k = U k U k implies that the marginal probabil- ity of sampling node v i , π i = U k δ i 2 2
, is proportional to the leverage scores p * i . The major difference between the i.i.d. leverage score approach and the DPP approach comes from the negative correlations induced by the DPP. In fact, the probability of jointly sampling nodes v i and v j in the DPP case is

π i π j -K 2 i j = π i π j -(δ i U k U k δ j ) 2 . The interaction term (δ i U k U k δ j )

Perspectives

Almost two decades have passed since spectral clustering was first introduced. Since then, a large body of work has attempted to accelerate its computation. So, has the problem been satisfactorily addressed? -or, despite all these works, is there still room for improvement and further research?

To answer, we must first define what "satisfactorily addressed" would entail. As we have seen, the prototypical spectral clustering algorithm can be divided in three sub-problems: the similarity graph computation runs in O(dn 2 ); the spectral embedding computation runs in O(t(ek + nk 2 )) using an Arnoldi algorithm with t implicit restarts and assuming that e is the number of edges; and the k-means step runs in O(tnk 2 ), with t now being a bound on the number of iterations of the Loyd-Max algorithm. Our criteria for evaluating an approximation algorithm aiming to accelerate one (or more) of these sub-problems are two-fold:

• We ask that the approximation algorithm's computation cost is effectively lighter than the cost of the sub-problem(s) it is supposed to accelerate! The ultimate achievement is an order-of-magnitude improvement w.r.t. n (or e), d and/or k, especially when the complexity has no hidden constants (i.e., the algorithm is practically implementable). When such a gain is not possible, a gain on the constants of the theoretical cost is also considered worthwhile. • The algorithm should come with convincing guarantees in terms of the quality of the found solution. Heuristics or partially motivated methods do not cut it. We require that, under mild assumptions, the proposed solution is provably close to the exact solution. Let us clarify two aspects of this statement further:

-It is difficult to concretely classify assumptions as mild, but a useful rule of thump is checking whether the theoretical results are meaningful for the significant majority of cases where spectral clustering would be used. -The control of the approximation error comes in different flavors, that we detail here from the tightest to the loosest. The best possible error control in our context is a control over the clustering solution itself, via error measures such as the misclustering rate. This is unfortunately unrealistic in many cases. An excellent alternative is the multiplicative error -considered as the gold standard in approximation theory-over the k-means cost 24 , ensuring that the cost of the approximation is not larger than 1 + ε times the cost of the exact solution. Next comes the additive error over the cost: ensuring that the cost difference between approximated and exact solutions is not larger than ε. All these error controls are referred to as end-to-end controls, and represent the limit of what we will consider a satisfactory error control.

Reviewing the literature, we were surprised to discover that there are rarely any algorithms meeting fully the proposed criteria: a faster algorithm with end-to-end control over the approximation error under mild assumptions. Let us revisit one by one the different approaches presented in Sections 3, 4 and 5 examining them in light of our criteria for success. In each category of approximation algorithms, we order the methods according to the power of their error control.

Sampling methods in the original feature space [Sec. 3].

• Representative points methods as described in [147,[START_REF] Huang | Spectral clustering with perturbed data[END_REF] allow for an end-to-end control on the miss-clustering rate ρ, which is unfortunately quite loose. The constants involved in Theorem 3 are in fact undefined -thus potentially large-, which is problematic knowing that ρ is by definition between 0 an 1. Also, the theorem's assumptions include independence of the ε i , which is hard to justify in practice. On the other hand, the computation gain of such methods is very appealing. • Feature projection methods, where the dimension d of the original feature space is reduced to a dimension d ≤ d based on Johnson-Lindenstrauss arguments, come with a multiplicative error control on the pairwise distances in the original feature space, thus providing a control on the obtained kernel matrix. The impact of this initial approximation on the final clustering result has not been studied. • Nyström-inspired methods [START_REF] Fowlkes | Spectral grouping using the Nystrom method[END_REF][START_REF] Li | Time and space efficient spectral clustering via column sampling[END_REF][START_REF] Bouneffouf | Sampling with minimum sum of squared similarities for nystrombased large scale spectral clustering[END_REF][START_REF] Mohan | Exploiting sparsity to improve the accuracy of nyström-based large-scale spectral clustering[END_REF] can be very efficient in practice especially because they do not need to build the graph. However, precisely because they do not build the graph, these methods cannot exactly perform two key parts of the prototypical spectral clustering algorithm: the k-NN sparsification and the exact degree computation. The partial knowledge and sparsity of the kernel matrix also makes sampling difficult, as using leverage scores sampling is not possible anymore, whereas most other sampling schemes do not work very well with sparse matrices and come with weak guarantees. To the extent of our knowledge, there is also no convincing mathematical argument proving that using these methods will yield a clustering that is of similar quality to that produced by the exact spectral clustering algorithm. • Random eigenspace projection is a very fast method and has been rigorously analyzed [START_REF] Martin | Fast approximate spectral clustering for dynamic networks[END_REF]133,105,111]. It is true that a successful application depends on obtaining a good estimate of the k-th eigenvalue, which is very hard when the k-th eigenvalue gap is relatively small. Nevertheless, our current understanding of spectral clustering suggests that it only works well when the gap is (at least) moderately large. As such, though there are definitely situations in which random eigenspace projection will fail to provide an acceleration, these correspond to cases where one should not be using spectral clustering in the first place. The same argumentation can also be used in defense of all methods that come with mild gap assumptions (see coarsening, and spectral sparsification). • Simple coarsening methods, such as the heavy-edge matching heuristic [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF], have nearly-linear complexity, seem to work well in practice, and are accompanied by end-to-end additive error control [START_REF] Loukas | Spectrally approximating large graphs with smaller graphs[END_REF]. Nevertheless, the current analysis of these heuristics only accounts for very moderate reductions (m ≥ n/2) and thus does not fully prove their success: in real implementations coarsening is used in a multi-level fashion resulting to a drastic decrease in the graph size (m = O(n/2 c ) for c levels), whereas the end-to-end control only works for a single level. • Advanced coarsening methods, such as local variation methods [START_REF] Loukas | Graph reduction with spectral and cut guarantees[END_REF], come with much stronger guarantees that allow for drastic size reduction and acceleration. Yet, thus far, all evidence suggests that finding a good enough coarsening is computationally as hard as solving the spectral clustering problem itself. As a consequence, it is at this point unclear whether these methods can be used to accelerate spectral clustering.

• Spectral sparsification techniques come with excellent guarantees in theory: one may prove that a spectral sparsifier can be computed in nearly-linear time and, moreover, the latter's spectrum will be provably close to the original one. Yet, we have reasons to doubt their practicality. Indeed, current algorithms are very complex, feature impractically large constants, and are only relevant for dense graphs.

In addition, spectral sparsifiers, by definition, approximate the entire spectrum of a graph Laplacian matrix. However, spectral clustering only needs an approximation of a tiny fraction of the spectrum. From that perspective, it is reasonable to conclude that without modification current approaches will not yield the best possible approximation. • Nyström-approximation applied directly to the Laplacian matrix is a good option, especially when combined with leverage score sampling. Nevertheless, an endto-end error control has only been partially derived and is not yet satisfactory.

Sampling to accelerate the k-means step [Sec. 5].

• Exact methods to accelerate the Lloyd-Max algorithm, may they be via avoiding unnecessary distance calculations or via a careful initialization are always useful and should be taken into account. • Coresets come with the strongest guarantees: the minimum number of samples to guarantee a (1 + ε) multiplicative error on the cost function has been well studied. Nevertheless, practical coreset sampling methods are scarce; and in the best cases, the sampling cost is of the same order of the Lloyd-Max running cost itself. • Graph-based sampling comes with strong guarantees, but not over the k-means cost: on the reconstruction error based on a k-bandlimited model that is only an approximation in practice. Moreover, we interpret the reduced indicator vectors z r obtained by running Lloyd-Max on the sampled set S as (possibly noisy) measurements of z on S. This interpretation currently lacks solid theoretical ground and impedes an end-to-end control of this approximation method. Nevertheless, the leverage-score based sampling allows for a reduction in order of magnitude of the Lloyd-Max running cost. • Other methods to accelerate k-means are not always appropriate to the spectral clustering context. Spectral feature dimension reduction is unnecessary in our context where d = k, sketching methods appropriate to distributed cases where n is very large are not appropriate neither as the spectral features need centralized data to be computed in any case.

In practice. The attentive reader will have remarked that, unsurprisingly, the tighter the error control, the more expensive the computation, and vice versa. Also, although we have put here an emphasis on the approximation error controls, it should not undermine the fact that methods from the whole spectrum are in practice useful, depending on the situation at hand, and specifically on the range of values of n, d and k. In very large d situations, a first step of random projection (or feature selection if some features are suspected to be too noisy) should be considered. Then, in situations where the exact computation of the proximity graph is too expensive, one may resort either to sketching methods or to Nyström-type methods to decrease the cost from quadratic to linear in n, and directly obtain an approximation of the spectral embedding without any explicit graph construction. These methods, however, do not take into account a sparsity constraint on the proximity graph and are usually rough on the degree correction they make. The role of the sparsity constraint is not well understood theoretically, but seems to be important in some practical cases [138]. In such instances, a better option is to use approximate nearest neighbours methods to create a sparse similarity graph, and work from there. In extremely large data, say n ≥ 10 8 , the only workable methods are the representative-based, with, if possible, a first k-means (or compressive kmeans [START_REF] Keriven | Compressive K-means[END_REF]) to reduce n to m ; or, in last resort, a uniform random sampling strategy.

In situations where one has to deal with such a large similarity graph that Arnoldi iterations are too expensive to compute the spectral embedding (either a graph created via approximate nearest neighbours, or if the original data is a graph), projection methods such as in [133,[START_REF] Boutsidis | Spectral clustering via the power methodprovably[END_REF], coarsening methods such as in [START_REF] Loukas | Graph reduction with spectral and cut guarantees[END_REF], or Nyströmbased methods are different possibilities.

Sampling methods to accelerate the last k-means step may seem to be a theoretical endeavour given that the Lloyd-Max algorithm is already very efficient. Due to the quadratic term in k, it is nevertheless in practice useful when k grows large. In this situation, hierarchical k-means [103] is a nice option. Coresets, because they are so stringent on the error control, have a hard time actually accelerating k-means, unless hybrid coreset-inspired methods are envisioned [START_REF] Frahling | A fast k-means implementations using coresets[END_REF]. Finally, graph-based methods, because they take into account that spectral features are in fact derived from the graph itself, enable significant acceleration and are well-suited to the spectral clustering context. Future research. Different directions of research could be envisioned to improve the state-of-the-art:

• For Nyström-inspired methods in the context of Sec. 3 (directly applied on the original data) as well as the other methods based on computing a low-rank approximation of the kernel matrix K, further work is needed to control both the sparsification and the degree correction, in order to bridge the gap between a provably good low-rank approximation of K to a provably good low-rank approximation of R. • For Nyström methods in the context of Sec. 4 (applied on a known or wellapproximated similarity graph), it would be interesting to extend Theorem 2 (for instance) to a control over A k -Ãk instead of A -Ã . This would enable a tighter use of Davis-Kahan's perturbation theorem in the discussion of Sec. 4.1 and, in fine, a better end-to-end guarantee. • Projection-based methods of Sec. 4.3.2 currently necessitate to compute a value λ * known to be in the interval [λ k , λ k+1 ). The algorithm used to do so is based on eigencount techniques that turn out to require as much computation time as the Lanczos iterations needed to compute U k exactly. One should relax this constraint to obtain end-to-end guarantees as a function of the distance between a coarsely estimated λ * and the target interval.

• The derivation and analysis of randomized multi-level coarsening schemes with end-to-end guarantees is very much an open problem. We suspect that, by utilizing spectrum-dependent sampling-schemes akin to leverage-scores one should be able to achieve results superior to heavy edge-matching in nearly linear time. • There is an interesting similarity between coreset techniques and the graph-based sampling strategies discussed in Sec. 5.3 and it would be interesting to investigate this link theoretically, maybe paving the way to coresets for spectral clustering?

Finally, accelerating the prototypical spectral algorithm depicted in Algorithm 1 should not be the sole objective of researchers in this field. Indeed, taking the graph cut point-of-view of Sec. 2.2, Algorithm 1 makes three insufficiently motivated choices: (i) To begin with, the sparsification step in Algorithm 1 is not well understood. Apart from the fact that it is always computationally more convenient to work with a sparse similarity graph then a dense one, the precise effect of sparsification on the clustering performance has not been analyzed. (ii) As mentioned in Sec. 2.2.3, the relaxation employed by spectral relaxation is not unique. Why should we focus our attention on this one versus another? See for instance [START_REF] Bresson | Multiclass Total Variation Clustering[END_REF]112] for recent alternative options. (iii) Finally, the use of k-means on the spectral features is not yet fully justified. Most of the end-to-end guarantees presented here compare the k-means cost of the exact solution to the k-means cost of the approximate solution. Given that the very use of k-means is not theoretically grounded, this choice of guarantee is debatable. Other options, such as a control over the rcut or ncut objectives are possible (as in [START_REF] Mohan | Beyond the nystrom approximation: Speeding up spectral clustering using uniform sampling and weighted kernel k-means[END_REF]) and should be further investigated.

Fig. 1

 1 Fig. 1 Left: the two half-moons synthetic dataset (n = 500, d = 2, k = 2). Right: k-means with k = 2 directly on P is unsuccessful to separate the two half-moons.

  [138]): min C∈R n×k tr C LC s.t. C C = I and C as in (4)

number of samples m, desired rank k 1 .

 1 Let S be m column indices chosen by some sampling procedure. 2. Denote by B = A(S, S) ∈ R m×m and C = A(:, S) ∈ R n×m the sub-matrices indexed by S. 3. Let B = QΣQ be the eigen-decomposition of B with the diagonal of Σ sorted in decreasing magnitude. 4. Compute the rank-k approximation of B as

2 .

 2 Run spectral clustering on the representatives. 3. Lift the solution back to the entire dataset. Let us illustrate this with the example of KASP: Algorithm 4. KASP: k-means-based approximate spectral clustering [147] Input. A set of n points P = (p 1 , p 2 , . . . , p n ) in dimension d, a number of desired clusters k, and a number of representatives m. 1. Perform k-means with k = m on P and obtain: a. the cluster centroids Y = (y 1 , . . . , y m ) as the m representative points.

3 . 1

 31 are also applicable here. Let us concentrate on the choice R = L n to illustrate the main ideas. As explained in Sec. 3.1, the tailing k eigenvectors U k of L n , can be interpreted as the top k eigenvectors of the PSD matrix A = 2I -L n . As such, the span of the top-k eigenvectors of Ãk , span( Ũk ), obtained by running Algorithm 3 on A should approximate the span of U k . Now, how does ones goes from Nyström theorems such as Theorem 2 to error bounds on the k-means cost function?

Algorithm 6 :

 6 D 2 -sampling. Input. X, m the number of required samples 1. Initialize B with any x chosen uniformly at random from X. 2. Iterate the following steps until B contains m elements: a. Compute d i = min b∈B x ib 2 2 . b. Define the probability of sampling x i as d i / ∑ i d i . c. Sample x new from this probability distribution and add it to B. Output: B a sample set of size m.

Algorithm 7 . 1 . 2 .

 712 Coresets to avoid k-means on X. Input. X, sampling set size m, and number of clusters k ≤ m. Compute a weighted coreset S of size m using a coreset-sampling algorithm. Run the Lloyd-Max algorithm on the weighted set S to obtain the set of k centroids C. 3. "Closest-centroid lifting": classify the whole dataset X based on the Voronoi cells of C. Output: A set of k centroids C = (c 1 , . . . , c k ).

4

 4

U k δ i 2 2

 2 -based methods can be summarized in the following general methodology [133]: Algorithm 9. Graph-based sampling strategies to avoid k-means on X. Input. X, m the number of required samples, k the number of desired clusters 1. Choose the random sampling strategy. Either: a. uniform (i.i.d.) Draw m i.i.d. samples uniformly. b. leverage score (i.i.d.) Compute ∀x i , p * i = /k. Draw m i.i.d. samples from p * . (optional:) set the weight of each sample s to 1/p * s . c. DPP Sample a few times independently from a DPP with kernel K k = U k U k . (optional:) set the weight of each sample s to 1/π s . 2. Run the Lloyd-Max algorithm on the (possibly weighted) set S to obtain the k reduced cluster indicator vectors z r ∈ R m . 3. Lift each reduced indicator vector {z r } =1,...,k to the full graph either with a. Least-square Solve (33) with y ← z r . b. Tikhonov Solve (34) with y ← z r .

•=

  If one uses the unbiased decoder ẑ = arg min w∈span(U k ) MP -1/2 M , then the following reconstruction result holds [109]:

  R 2m×n the RFF matrix. Other embeddings are possible in the RFF framework, but this one was shown to be the most appropriate to the Gaussian kernel[127]. As m increases, ψ i ψ j concentrates around its expected value κ(p i , p j ): ψ i ψ j κ(p i , p j ). Proposition 1 of [127] states the tightness of this concentration: it shows that the approximation starts to be valid with high probability for m ≥ O(d log d). The Gaussian kernel matrix is thus well approximated as K Ψ Ψ. With such a low-rank approximation Ψ of K, one can: estimate the degrees

  [125]. Most fast algorithms for spectral sparsification entail sampling O(n log n) edges from the total edges present in the graph. Different sampling schemes are possible[124,[START_REF] Koutis | Approaching optimality for solving sdd linear systems[END_REF], but the most popular ones entail sampling edges with replacement based on their effective resistance. It should be noted that though computing all effective resistances exactly can be computationally prohibitive, the effective resistance of edges can be approximated in nearly linear-time on the number of edges based on a Johnson-Lindenstrauss argument[124].

  5.1. • Approximate acceleration of Lloyd-Max. Approximately accelerating the Lloyd-Max algorithm has also received attention for instance via approximate nearest neighbour methods [108], via cluster closure [140], or via applying Lloyd-Max hierarchically (in the large k context) [103]. An approach involving sampling is introduced in [119]

  2 will be typically large if v i and v j are in the same cluster, and small if not. In other words, different from the i.i.d. leverage score case where each new sample is drawn regardless of the past, the DPP procedure avoids to sample nodes containing redundant information.• Whereas the leverage score approach only guarantees a RIP with high probability after O(k log k) samples, the DPP approach has a stronger deterministic guarantee: it enables perfect invertibility (up to the noise level) after precisely m = k samples. The reconstruction guarantee of (36) is nevertheless not satisfactory: even corrected by the marginal probabilities P S , the matrix U k M P -1 S MU k can still have a very small λ min , such that reconstruction may be quite sensitive to noise. Improving this control is still an open problem. In practice, sampling independently 2 or 3 times from a DPP with kernel K k , creates a set S of size 2k or 3k that is naturally more robust to noise. • Whereas independent sampling is straightforward, sampling from a DPP with arbitrary kernel costs in general O(n 3 ) (see Alg. 1 of [74] due to [62]). Thankfully, in the case of a projective kernel such as K k , one can sample a set in O(nk 2 ) based on Alg. 3 of [130].

•

  Sketching methods such as the Random Fourier Features[110] is yet another way of obtaining a pointwise multiplicative (1 + ε) error on the Gaussian kernel computation. RFF enable to compute a provably good low-rank approximation of the kernel. They nevertheless suffer from the same problems as Nyström-based techniques: without building the graph, sparsification and degree-normalization are uncontrolled. In addition, the guarantees on the low-rank approximation of the kernel do not transfer easily to guarantees of approximation of the spectral embedding U k .• Approximate nearest neighbours methods are numerous and varied, and come with different levels of guarantees. Practical implementations of algorithms, however, often set aside theoretical guarantees to gain on efficiency and performances; and comparisons are usually done on benchmarks rather than on theoretical performances. In the best of cases, there is a control on how close the obtained nearest neighbour similarity graph is to the exact one; but with no endto-end control.

Spectral embedding approximation methods

[Sec. 4]

.

Out-of-sample extensions of spectral clustering do exist (see for instance Section 5.3.6 of [144]), but they require additional work.

Each node v i of V represents a point p i , an undirected edge exists between nodes v i and v j if and only if W(i, j) = 0, and the weight of that connection is W(i, j).

In some of these examples, the k largest eigenvalues (instead of the k lowest in the Laplacian cases) of the representative matrix, and especially their corresponding eigenvectors, are of interest. This is only a matter of sign of the matrix R and has no impact on the general discussion.

By definition, a partitionP = {V 1 , V 2 , . . . , V k } of the nodes V is such that ∪ =1,...,k V = V and ∀ = , V ∩ V = /

The reader should note that in the graph theory literature, the measure of conductance is preferred over ncut. Conductance is max w(V , V )/w(V ). The two measures are equivalent when k = 2.

To construct an example possibly verifying such a strong gap assumption, consider k cliques of size k connected together via only k -1 edges, so as to form a loosely connected chain. Even though this is a straightforward clustering problem known to be easy for spectral clustering algorithms, the above theorem's assumption implies λ k+1 = Ω (k 2 ncut(P * )) = Ω (k) which, independently of n, can only be satisfied when k is a small (recall that the eigenvalues of L n are necessarily between 0 and 2).

Usually, one needs to ensure that ∑ i≤k, j>k ( ũ i u j ) 2 /k remains bounded.

The combinatorial Laplacian of a complete balanced binary tree on k ≥ 3 levels and n = 2 k -1 nodes has[START_REF] Achlioptas | Database-friendly random projections: Johnson-lindenstrauss with binary coins[END_REF] n ≤ λ 2 ≤ 2 n[START_REF] Guattery | On the performance of spectral graph partitioning methods[END_REF].

Authors in[116] have a very similar proposition as[START_REF] Chen | Large Scale Spectral Clustering with Landmark-Based Representation[END_REF], adding a projection phase at the beginning to reduce the dimension d (see Sec. 3.4.2). Similar ideas may also be found in[137].

It is in fact O(mdnt) for step 1, and bounded by O(dm 2 + m 2 k + mk 2 ) for step 2. As n ≥ m and m ≥ k, the total complexity is bounded by O(mdnt + m 3 ).

The recall rate for a node is the number of correctly identified k-NN divided by k. The recall rate for a k-NN graph is the average recall rate over all nodes.

Based on three arguments: (i) for any two matrices M 1 and M 2 of rank r 1 and r 2 it holds that rank(M 1 + M 2 ) ≤ r 1 + r 2 , (ii) for any matrix M or rank r, M F ≤ √ r M 2 , and (iii) both U k and Ũk are of rank k.

Note that, in our setting, A = 2I -L n and σ k = 2λ k .

When one starts from a set of points, it is preferable to sparsify the graph by retaining a constant number of nearest neighbors for each point. The resulting nearest neighbor graph has already O(n) edges, which is the smallest possible.

A weaker version of this definition exists in the literature where the ε-approximation is only required for C * .

Generalizing Algorithm 2 to a weighted set is straightforward: in step 2b, instead of computing the center of each cluster, compute the weighted barycenter.

It takes longer to find a good sample than to run k-means on the whole dataset!

Note that the complexity is different from the leverage score computation of the Nyström techniques of Sections 3.1 and 4.1 because, here, we suppose U k known whereas U k was not known in the previous sections. With U k known, computing the leverage scores only entails computing the normalized energy of each line of U k .

A control in terms of the k-means cost is usually considered as k-means is the last step of spectral clustering. Nevertheless, recalling the minimum cut perspective of Sec. 2.2, the control should arguably be in terms of rcut or ncut costs.

Acknowledgments

This work was partially supported by the CNRS PEPS I3A (Project RW4SPEC) and the Swiss National Science Foundation (Project DLEG, grant PZ00P2 179981).