Daniele Calandriello

Ioannis Koutis

Alessandro Lazaric

Michal Valko

Improved large-scale graph learning through ridge spectral sparsification

The representation and learning benefits of methods based on graph Laplacians, such as Laplacian smoothing or harmonic function solution for semi-supervised learning (SSL), are empirically and theoretically well supported. Nonetheless, the exact versions of these methods scale poorly with the number of nodes n of the graph. In this paper, we combine a spectral sparsification routine with Laplacian learning. Given a graph G as input, our algorithm computes a sparsifier in a distributed way in O(n log 3 (n)) time, O(m log 3 (n)) work and O(n log(n)) memory, using only log(n) rounds of communication. Furthermore, motivated by the regularization often employed in learning algorithms, we show that constructing sparsifiers that preserve the spectrum of the Laplacian only up to the regularization level may drastically reduce the size of the final graph. By constructing a spectrally-similar graph, we are able to bound the error induced by the sparsification for a variety of downstream tasks (e.g., SSL). We empirically validate the theoretical guarantees on Amazon co-purchase graph and compare to the state-of-the-art heuristics.

Introduction

Graphs are a very effective data structure to represent relationships between entities (e.g., social and collaboration networks, influence graphs). Over the years, many machine learning problems have been defined and solved exploiting the graph representation, such as graph-regularized least squares (LAPRLS, [START_REF] Belkin | On manifold regularization[END_REF], Laplacian smoothing (LAPSMO, [START_REF] Sadhanala | Graph sparsification approaches for laplacian smoothing[END_REF]) graph semi-supervised learning (SSL, [START_REF] Chapelle | Semi-supervised learning[END_REF][START_REF] Zhu | Semi-Supervised Learning Using Gaussian Fields and Harmonic Functions[END_REF], laplacian embedding (LE, Belkin & Niyogi 2001, and spectral Proceedings of the 35 th International Conference on Machine Learning, Stockholm, Sweden, PMLR 80, 2018. Copyright 2018 by the author(s).

clustering (SC, Von Luxburg 2007). The intuition behind graph-based learning is that the information expressed by the graph helps to capture the underlying structure of the problem (e.g., a manifold), thus improving the learning. For instance, LAPSMO and SSL rely on the assumption that nodes that are close in the graph are more likely to have similar labels. Similarly, LE and SC try to find a lowdimensional representation of the nodes using the eigenvectors of the Laplacian of the graph. In general, given a graph G of n nodes and m edges, most of graph-based learning tasks require computing the minimum of a cost function based on the associated n × n Laplacian matrix L G , which contains m non-zero entries. Solving exactly such optimization problems amounts to O(n 3 ) time and O(n 2 ) space complexity in the worst case and they become infeasible even for mildly large/dense graphs.

A complete review of the literature on large-scale graph learning is beyond the scope of this paper and we only consider methods that reduce learning space and time complexity starting from a given graph received as input. 1 We identify mainly three possible approaches. We can (1) reduce runtime replacing the pseudo-inverse operator L + G with an iterative solver, (2) reduce time and space complexity replacing the large graph G with a sparser approximation H, or (3) reduce runtime and increase memory capacity by distributing the computation across multiple machines.

Iterative solvers. Iterative methods can solve a number of learning problems without explicitly constructing L + G (e.g., gradient descent, GD, for LAPSMO, iterative averaging for SSL, and the power method for SC). In this case we only need O(m) time per iteration. Unfortunately, all simple iterative methods (e.g., GD) converge in a number of iterations proportional to the condition number of the Laplacian, κ = λ max (L G )/λ min (L G ), which may grow linearly with the number of nodes n, thus removing the advantage of the iterative method, whose complexity tends to O(n 3 ) in the worst case. Advanced iterative methods, such as the preconditioned conjugate gradient, use preconditioning to find an accurate solution in a number of iterations independent of κ. [START_REF] Koutis | A nearly-m log n time solver for SDD linear systems[END_REF] gives a nearly-linear solver for Laplacians or strongly diagonally dominant (SDD) matri-ces, that using a chain of preconditioners, converges in only O(m log(n)) time. As space and time costs scale with the number of edges, a natural desire is to reduce m by sparsifying and distributing the graph.

Graph sparsification. The objective of sparsification methods is to remove redundant edges, so that the resulting sparse sub-graph can be easily stored in memory and efficiently manipulated to compute final solutions. A simple graph-sparsification technique is to sample nq (with q > 1) edges from G with probabilities proportional to the edge weights with replacement. While computationally very efficient, uniform sampling requires sampling a number of edges proportional to O(nµ(G)) (i.e., q ∝ µ(G)), where µ(G) is the coherence of the Laplacian matrix, and it can grow as large as n when the graph is highly structured (e.g., if there is a single edge e connecting two components of the graph we need to sample all of the edges of the graphpotentially O(n 2 )-to guarantee that we do not exclude e and generate an inappropriate H). A more refined approach is the k-neighbors (kN) sparsifier [START_REF] Sadhanala | Graph sparsification approaches for laplacian smoothing[END_REF], which performs local sparsifications node-by-node by keeping all edges at nodes with degree smaller than q, and samples them proportionally to their weights whenever the degree is bigger than q. While in certain structured graphs, this method may perform much better than uniform (Von [START_REF] Luxburg | Hitting and commute times in large random neighborhood graphs[END_REF], in the general case q, still needs to scale with the coherence µ(G). A more effective method is to sample edges proportionally to their effective resistance, which intuitively measures the importance of an edge in preserving the minimum distance between two nodes. As a result, only relevant edges are kept and the sparsified graph could be reduced to O(n polylog(n)) edges. Nonetheless, computing effective resistances also requires the pseudo-inverse L + G , thus being as expensive as solving any graph-Laplacian learning problem.

Distributed computing. When the number of edges m is too large to fit the whole graph in a single machine, we are forced to distribute the edges across multiple machines. At the same time, if the sparsifier construction or the downstream inference can be parallelized, we can also reduce their runtime. Unfortunately, distributing data and computation across multiple machines can cause large communication costs. For example, simple GD or label propagation methods require O(κ) iterations (and communication rounds) to converge and access to non-local (e.g., neighbors in a graph) data. While preconditioned solvers reduce the number of iterations, almost none of their memory access is local, thus making difficult to have efficient distributed implementations.

Contribution. In this paper, we propose a new approach that aims at integrating the benefits of the three different methods above. Using the large memory and computational capacity of distributed computing and leveraging the sequen-tial sparsification methods of [START_REF] Kelner | Spectral sparsification in the semi-streaming setting[END_REF] and [START_REF] Calandriello | Distributed sequential sampling for kernel matrix approximation[END_REF], we show how to compute an accurate sparsifier H of graph G in O(n log 3 (n)) time, O(n log 2 (n)) work and O(n log(n)) memory, using only log(n) rounds of communication. Afterwards, learning tasks can be solved directly on L H on a single machine using near-linear time solvers, resulting in an overall O(n log 3 (n)) runtime. Moreover, we show that the regularization used in some graphbased learning algorithms allows using even sparser graphs. In particular, we introduce the notion of ridge effective resistance to obtain sparsifiers that are better adapted to solve Laplacian-regularized learning tasks (e.g., LAPSMO, SSL) and are smaller than standard spectral sparsifiers without compromising the performance of downstream tasks.

Background

We use lowercase letters a for scalars, bold lowercase letters a for vectors and uppercase bold letters A for matrices. We use A B to denote that B -A is positive semidefinite (PSD), [A] i,j to indicate the (i, j)-th entry of A, and ordered the eigenvalues as λ 1 (A) ≤ . . . ≤ λ n (A).

Graphs and graph Laplacian

We denote with G = (V, E), an undirected weighted graph with n nodes V and m edges E. Each edge e i,j ∈ E has a weight a ei,j measuring the "similarity" between nodes i and j. Given graphs G and G over the same set of nodes V, G + G denotes the graph obtained by summing the weights of their edges. For graph G, we introduce the weighted adjacency matrix A G with entries [A G ] i,j = a ei,j , the total weights A = e a e , and the diagonal degree matrix D G with entries [D G ] i,i j a ei,j . The Laplacian of G is the PSD matrix L G D G -A G . Furthermore, we assume that G is connected and thus L G has only one eigenvalue equal to 0 and Ker(L G ) = 1. Let L + G be the pseudoinverse of L G and L -1/2 G = (L + G ) 1/2 . For any node i = 1, . . . , n, we denote with χ i ∈ R n , the indicator vector, so that b e √ a e (χ iχ j ) is the "edge" vector. If we denote with B G the m × n signed edge-vertex incidence matrix, then the Laplacian can be written as

L G = e b e b T e = B T G B G .

Learning on graphs

Given graph G and its Laplacian L G , we denote with f ∈ R n , a labeling of its nodes, where [f ] i is the value associated with the i-th node. Many graph learning algorithms assume that the optimal labeling f is smooth w.r.t. the graph, i.e., the quantity e a e ([f ]

ei -[f ] ej ) 2 = f T L G f is small.
In the following, we review examples from the supervised, semi-supervised and unsupervised learning with graphs.

Laplacian smoothing (LAPSMO) with Gaussian noise.

Given a graph G on n nodes, let y f + ξ be a noisy measurement of f with [ξ] i ∼ N (0, σ 2 ). The goal of LAPSMO is to find a vector f that accurately reconstructs f under the graph smoothness assumption by solving

f arg min f ∈R n (f -y) T (f -y) + λf T L G f = (λL G + I) -1 y, (1)
where λ is a regularization parameter.

Graph semi-supervised learning (SSL). In SSL, the input f is a partial observation of the labels f for a subset S ⊂

[n] of nodes. The goal is to predict the labels f u of the unrevealed nodes. The harmonic function solution (HFS) by [START_REF] Zhu | Semi-Supervised Learning Using Gaussian Fields and Harmonic Functions[END_REF] solves the optimization problem

f HFS arg min f ∈R n 1 (f -y) T S (f -y) + λf T L G f = (λ L G + I S ) + y S , (2)
where |S| is the number of labeled nodes received as input, I S ∈ R n×n is the identity matrix with zeros at nodes not in S, and y S I S y ∈ R n . Similarly, in local transductive regression (LTR) [START_REF] Cortes | Stability of transductive regression algorithms[END_REF], the optimization problem is

f LTR arg min f ∈R n (f -y) T C(f -y) + f T (L G + λI)f = (C -1 (L G + λI) + I) -1 y S , (3)
where C is a diagonal matrix with entries c for nodes in S, c u for entries not in S, and c ≥ c u > 0.

Spectral clustering (SC).

Applying the Laplacian smoothness assumption, the goal of SC is to find k disjoint subset assignments such that the clusters are smooth w.r.t. the Laplacian. Let {f c } k c=1 be the cluster indicator vectors such that [f c ] i 1 if node i is in the c-th cluster and [f c ] i 0 otherwise. Denote with F ∈ R n×k , the matrix containing the assignments, and let C be the space of feasible clustering, such that all f c are binary and each row of F contains only one non-zero entry. Since computing the minimum ratio-cut is NP-hard (Von Luxburg, 2007;[START_REF] Lee | Multiway spectral partitioning and higher-order cheeger inequalities[END_REF], even under constraints [START_REF] Cucuringu | Simple and scalable constrained clustering: a generalized spectral method[END_REF], SC defines instead the relaxed problem

F arg min F:F T F=I k ,fc⊥1
Tr(F T LF).

Once the relaxed solution is computed, we can use different heuristics to recover the clustering, such as thresholding or performing a k-means clustering on the F matrix.

Computational complexity. The problems above require either to compute an eigendecomposition of the Laplacian L G or to solve a linear system involving L G . Computing these exactly is not feasible when the number of nodes n and edges m grows. In particular, (a) storing L G in memory requires O(m) space, and it is not feasible when m is large, (b) even if L G is sparse and m is small, the pseudo-inverse L + G might be dense, and thus computing and storing L + G exactly requires up to O(n 3 ) time and O(n 2 ) space.

Distributed Spectral Sparsification

In this section, we describe a new, sequential, distributed, and efficient algorithm for graph sparsification that can be used as a preprocessing step to solve a large variety of downstream learning tasks, without significantly affecting their performance. We point out that while distributing data-agnostic sparsifiers (e.g. uniform sampling) is straightforward, distributing the computation of sparsifiers based on effective resistances requires a careful merging procedure to guarantee satisfactory memory vs. accuracy tradeoff, which is what we provide in this section.

(ε, γ)-spectral sparsifiers

We start with the introduction of the notion of (ε, γ)sparsifier that is adapted for the learning tasks that use sparsified graph Laplacian.

Definition 1. A (ε, γ)-spectral sparsifier of G is a re- weighted sub-graph H ⊆ G whose Laplacian L H satisfies (1 -ε)L G -εγI L H (1 + ε)L G + εγI. (4)
For γ = 0, this definition reduces to the standard notion of ε-spectral sparsifier [START_REF] Spielman | Spectral sparsification of graphs[END_REF]. The main difference is that an (ε, γ)-spectral sparsifier allows for an extra additive error of order εγ. This change is directly motivated by the fact that the sparsifier H may be used in learning tasks whose solution may not be sensitive to small (additive) errors. As a result, (ε, γ)-spectral sparsifiers are able to further reduce the size of H w.r.t. (ε, 0)sparsifiers, without significantly affecting the final learning performance. Formally, an ε-sparsifier preserves all the quadratic forms up to a small multiplicative (constant) error, and thus can be used to provide an accurate approximation to many important quantities such as graph cuts or eigenvalues. In fact, for all i ∈ [n], an ε-sparsifier guarantees that

(1 -ε)λ i (L G ) ≤ λ i (L H ) ≤ (1 + ε)λ i (L G ).
Nonetheless, in many learning tasks (e.g., LTR) the noise level in the signal f requires regularizing the solution so that the Laplacian L G itself is eventually replaced by L G + λI (e.g., Eq. 1). This corresponds to soft-thresholding the eigenvalues of the Laplacian, so that eigenvalues below λ are partially ignored. If λ is properly tuned w.r.t. the noise, the regularization increases stability and improves the learning performance. Therefore, constructing a sparsifier that accurately reconstructs all eigenvalues of L G may be wasteful, as it may require keeping most of the edges. As a result, in tasks where L G is regularized, it is better to use (ε, γ)-sparsifiers, as their additive error γI is homogeneous with the regu-larization λI and their smaller size allows scaling to up. 2We now extend the results of [START_REF] Spielman | Graph sparsification by effective resistances[END_REF] for the construction of ε-spectral sparsifiers to the general case of (ε, γ)-sparsifiers. We redefine the edge effective resistance to account for the regularization. Definition 2. The γ-effective resistance of an edge e in graph G is defined as

r e (γ) b T e L G + γI -1 b e . (5)
The "effective dimension" of the graph is the total sum of the γ-effective resistances, d eff (γ) e r e (γ).

We can now construct a sparsifier H by sampling q times each edge with a probability proportional to its γ-effective resistance. More formally, the resulting (random) graph contains q e ∼ B(r e (λ); q) copies of each edge, where B is the Binomial distribution, and its associated Laplacian is L H = e∈H q e /(qr e (γ))b e b T e , which is an unbiased estimator of L G . We can then apply existing results from sketching of PSD matrices [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] to prove that H is a valid (ε, γ)-sparsifier. Proposition 1 [START_REF] Cohen | Input sparsity time low-rank approximation via ridge leverage score sampling[END_REF]. Let ε > 0 and γ ≥ 0 be the accuracy parameters and 0 ≤ δ ≤ 1 the probability of error. Let H be the graph obtained by sampling edges in G with a probability proportional to their γ-effective

resistances. If q ≥ 4 log(4n/δ)/ε 2 , then w.p. 1 -δ, H is an (ε, γ)-sparsifier with O(d eff (γ)q) edges.
We first notice that this result reduces to the one of [START_REF] Spielman | Graph sparsification by effective resistances[END_REF] for γ = 0. In fact, d eff (0) = n -1 for all graphs, thus matching the space requirement q for ε-sparsifiers. Nonetheless, as γ increases, the size of H reduces significantly. Using L G = B T G B G , the effective dimension d eff (γ) can be conveniently rewritten as

d eff (γ) = Tr B T G B G (B T G B G +γI) -1 = n i=2 λ i (L G ) λ i (L G ) + γ ,
thus showing that d eff (γ) is the "soft" rank of the Laplacian, where γ significantly reduces the contribution of small eigenvalues to the total sum. While in the worst case d eff (γ) can be as large as n -1, for a variety of graphs with rapidly decaying spectrum [START_REF] Jamakovic | The Laplacian spectrum of complex networks[END_REF][START_REF] Samukhin | Laplacian spectra of, and random walks on, complex networks: Are scale-free architectures really important[END_REF][START_REF] Zhan | On the distributions of Laplacian eigenvalues versus node degrees in complex networks[END_REF][START_REF] Akoglu | Graph based anomaly detection and description: a survey[END_REF], d eff (γ) may be significantly smaller than n -1, thus reducing the number of edges q required to obtain an (ε, γ)-sparsifier.

The algorithm

As pointed out in the introduction, the main limitation of effective-resistance-based sparsification is that the computation of r e requires inverting the Laplacian matrix, thus resulting in a computational cost that already matches the cost H 1, ← G ← {(e i,j , q e = 1, p 1,e = 1)}

3: Initialize set S 1 = {H 1, } k =1 4: for h = 1, . . . , k -1 do 5: Pick two sparsifiers H h,i , H h,i from S h 6: H ← Merge-Resparsify(H h,i , H h,i ) 7:
Place H back into S h+1 8: end for 9: Return H G , the last sparsifier in S k Algorithm 2 Merge-Resparsify Require:

(ε, γ)-sparsifiers H h,i , H h,i of graphs G h,i , G h,i Ensure: H, an (ε, γ) sparsifier of G h,i + G h,i 1: Initialize H = H h,i + H h,i 2:
For all e ∈ H, use a fast SDD solver to compute

r h+1,e (γ) ← (1 -ε)b T e (L H + (1 + ε)γI) -1 b e
3: Set probabilities p h+1,e ← min{ r h+1,e (γ), p h,e } 4: Sample q h+1,e from B( p h+1,e / p h,e , q h,e ) 5: Return H ← {(e i,j , q h+1,e , p h+1,e )} for all q h+1,e > 0 of the learning tasks themselves. Moreover, large graphs cannot be stored in memory, and multiple passes over the graph would result in a disk access overhead larger than the computational cost. In order to avoid these problems, we adapt our previous work [START_REF] Calandriello | Distributed sequential sampling for kernel matrix approximation[END_REF] in online sparsification and randomized linear algebra (see a thorough discussion and comparison at the end of the section) to obtain the distributed sequential resparsification (DiSRe) algorithm (Alg. 1).3

The structure. We represent a sparsifier H as a collection of weighted edges H {(e i,j , q e , p e )}, and the Laplacian can be reconstructed as L H e∈H 1/ p e (q e /q)b e b T e . Intuitively, each edge e has an associated weight based on its probability p e , and a number of included copies q e . Keeping multiple copies of each edge helps the random L H to concentrate towards L G , where the maximum number of copies q for an edge trades-off success probability and the size of H. We assume we have k machines. DiSRe begins by partitioning the graph G into k sub-graphs G on n vertices and m ≥ n edges, such that G = {G } k i= In other words, it splits the matrix B G into submatrices B Gi by arbitrarily selecting a subset of rows. The sub-graphs are small enough that they can be stored in memory,4 and they are also obviously sparsifiers of themselves, therefore we can define an initial set of sparsifiers

H {4,1} H {3,1} + H {3,2} H {3,1} H {2,1} H {1,1} + H {1,2} H {1,1} G1 H {1,2} G2 H {3,2} H {2,2} + H {2,3} H {2,2} H {1,3} G3 H {2,3} H {1,4} G4 h = 4 h = 3 h = 2 h = 1
S 1 {H 1, } k =1 , with H 1,
{(e i,j , q 1,e = q, p 1,e = 1)} e∈G l . With this definition, H 1, contains edges e i,j with unit weight p 1,e = 1 and H 1, = G . Starting from these initial sparsifiers, DiSRe proceeds through a sequence of merge and sparsify operations where two sparsifiers are first combined and then sparsified again to keep having manageable-size graphs at each step. While DiSRe can run on any arbitrary sequence of merges, we consider the most (computationally) effective scheme, where sparsifiers are merged two-by-two in parallel, thus inducing a balanced full binary merge tree (see Fig. 2). For notational convenience, we consider that at each iteration h, the inner loop of Alg. 1 only merges two arbitrary sparsifiers from the pool of available sub-graphs S h and merges them into a new sparsifier. In practice, multiple merge-and-sparsify operations can be executed in a parallel and asynchronous way. The size of S h , number of sparsifiers present at layer h, is |S h | = kh + 1. Therefore, a node in the tree corresponding to a sparsifier is uniquely identified by two indices {h, } where h is the height of the layer and ≤ |S h | is the index of the node in the layer. We also define the graph G {h, } as the union of all sub-graphs G that are reachable from node {h, } as leaves (descendants of {h, }). For example, in Fig. 2, sparsifier H 3,1 in node {3, 1} approximates the graph G {3,2} = G 3 + G 4 , where we highlight in red the descendant tree.

The resparsification. In Alg. 2 we detail how two arbitrary sparsifiers are combined to obtain a temporary graph H. While the merge operation simply combines H h,i and H h,i by summing their weights, the resparsification aims at generating a valid sparsifier from the "original" sub-graph (G h,i + G h,i ), as if it was directly sparsified at the beginning. We first compute estimates r(γ) of the γ-effective resistance by using fast solvers to invert the strongly diagonal dominant L H + γI matrix. Instead of sampling edges in H directly proportionally to r(γ) (more precisely p h+1,e ), we perform a "resampling" scheme where an edge e is preserved with a "reweighted" probability p h+1,e / p h,e . Intuitively, the overall stored on a single machine), we can simply apply the same merging scheme of DiSRe by loading small enough chunks of the graph and sparsifying them sequentially. sequence of resampling guarantees that at each step h + 1, an edge e ∈ (G h,i + G h,i ) has the "correct" probability p h+1,e of being included in the sparsifier.

Performance. We now study the performance of DiSRe and its complexity. Time complexity refers to the amount of time necessary to compute the final solution and work complexity refers to the total amount of operations carried out by all machines to compute the final solution.

Theorem 1. Let ε > 0 be the accuracy, 0 ≤ δ ≤ 1 the probability of error, and ρ (1 + 3ε)/(1ε). Given an arbitrary graph G and an arbitrary merge tree structure, if DiSRe is run with parameter q 26ρ log(3n/δ)/ε 2 , then each sub-graphs H {h, } is an (ε, γ)-sparsifier of G {h, } with at most 3qd eff (γ) edges with probability 1δ. Whenever the merge tree is balanced and k is big enough such that m/k ≤ 3qd eff (γ),5 then merge operations can be run in parallel across the machines with an overall time complexity of O(d eff (γ) log 3 (n)), a total work O(m log 3 (n)), and O(log(n)) rounds of communication.

Discussion. [START_REF] Kelner | Spectral sparsification in the semi-streaming setting[END_REF] proposed a sequential algorithm for graph sparsification that closely emulates the batch sampling of [START_REF] Spielman | Graph sparsification by effective resistances[END_REF] in a semi-streaming setting and incrementally constructs an ε-sparsifier. However, their proof had a flaw since they treated dependent variables as independent [START_REF] Calandriello | Analysis of kelner and levin graph sparsification algorithm for a streaming setting[END_REF]. Kyng et al. (2016) resolved the issues in the proof of [START_REF] Kelner | Spectral sparsification in the semi-streaming setting[END_REF] and showed that a slightly modified algorithm can construct a sparsifier with O(n log(n)/ε 2 ) edges in O(m log 2 (n)/ε 2 ) time, matching the space complexity of batch sampling. The method proposed by Kyng et al. (2016) can be further improved by parallelizing its computation over multiple machines. Using the parallel sparsification algorithm of [START_REF] Koutis | Simple parallel and distributed algorithms for spectral graph sparsification[END_REF], the time complexity can be reduced up to O(log 6 (n)). Nonetheless, since these methods require random access to the edges, they cannot be easily distributed (it would have O(m polylog(n)) communication cost) and scaled to graphs that cannot be stored on a single machine. Furthermore, the algorithm of Kyng et al. (2016) accurately reconstructs the whole spectrum of the Laplacian, which leads to sparsifiers whose number of edges scales linearly with n. On the other hand, in regularized learning tasks, the presence of multiplicative and additive spectral error allows creating smaller sparsifiers whose size scales with d eff (γ). Notice that, for γ large enough, this possibly means sparsifiers with less than n -1 edges, necessarily leading to disconnected graphs. Finally, note that merging two traditional ε-sparsifiers gives an ε-sparsifier, merging two (γ, ε)sparsifiers produces a less accurate (2γ, ε)-sparsifier. Therefore simple merge-and-reduce strategies [START_REF] Feldman | Turning big data into tiny data: Constant-size coresets for k-means, pca and projective clustering[END_REF], which address every resparsification as independent, would either cumulate errors or require multiple passes over the data. Similarly to Kyng et al. (2016), DiSRe's sequential Merge-Resparsify solves this problem (Appendix B).

Mixed additive-multiplicative reconstruction is studied more extensively in randomized matrix algebra [START_REF] Drineas | Lectures on randomized numerical linear algebra[END_REF]. [START_REF] Cohen | Online row sampling[END_REF] developed an efficient method to spectrally sparsify generic matrices up to (1 ± ε) multiplicative and γ-additive errors using an incremental sampling method based on ridge leverage scores (i.e., the analog of γ-effective resistances for matrices). If applied to graph Laplacians, their method adds edges incrementally and returns an (ε, γ)-sparsified graph with [START_REF] Cohen | Online row sampling[END_REF] provided only ε-sparsifiers, suggesting to set γ as small as possible, and did not explore the advantages possible in machine learning. Moreover, no existing (ε, γ)-sparsifier construction method can leverage both distribution and fast solvers. [START_REF] Cohen | Online row sampling[END_REF] can only add edges (but not remove them as DiSRe), preventing repeated merge-and-resparsify. Other streaming RLS sampling methods, such as by [START_REF] Cohen | Input sparsity time low-rank approximation via ridge leverage score sampling[END_REF], use dense intermediate sketches, such as frequent directions [START_REF] Ghashami | Frequent directions: Simple and deterministic matrix sketching[END_REF], that are not Laplacians of a subgraph and cannot be easily paired with near-linear solvers for Laplacians.

O(d eff (γ) log 2 (n)) edges in O(m log(n)) time. Nonethe- less,

Downstream Guarantees

We now show how the spectral reconstruction guarantees provided by (ε, γ)-sparsifiers translate into guarantees on the quality of the approximate solutions computed using H instead of G. We first introduce a result for ε-sparsifiers in SSL and then show how for regularized problems, (ε, γ)sparsification can further improve computational performance without loss in accuracy in LAPSMO.

Generalization bounds for SSL

Given the closed form solutions of HFS (Eq. 2) and LTR (Eq. 3), we simply replace L G with L H and then run a nearlylinear time solver to obtain approximate solutions f HFS and f LTR . We compare approximate solutions to their exact counterparts in the context of algorithmic stability. Definition 3. Let L be a transductive learning algorithm. We denote by f and f the solutions obtained by running L on datasets V (S, T ) and V (S , T ) respectively. L is uniformly β-stable w.r.t. the squared loss if there exists β ≥ 0 such that for any two partitions (S, T ) and (S , T ) that differ by exactly one training (and test) point and for all

i ∈ [n], we have |([f ] i -[y] i ) 2 -([f ] i -[y] i ) 2 | ≤ β.
The stability of LTR was proven by [START_REF] Cortes | Stability of transductive regression algorithms[END_REF]. On the other hand, the singularity of the Laplacian may lead to unstable behavior in HFS due to the (γ L G + I S ) + pseudo-inverse, with drastically different results for small perturbations of the dataset. For this reason, we take the Stable-HFS algorithm by [START_REF] Belkin | Regularization and Semi-Supervised Learning on Large Graphs[END_REF], where an additional regularization term is introduced to restrict the space of admissible solutions to the space F {f : f , 1 = 0} of solutions orthogonal to the null space of L G (i.e., centered functions). As shown by [START_REF] Belkin | Regularization and Semi-Supervised Learning on Large Graphs[END_REF], to satisfy the constraint, it is sufficient to set an additional regularization parameter µ to µ ((γ L G + I S ) + y S ) T 1/((γ L G + I S ) + 1) T 1, and compute the solution f STA as f STA (γ L G + I S ) + (y S -µ1). While Stable-HFS is more stable and thus more suited for theoretical analysis, its space and time requirement remains O(m) and cannot be applied to graphs with a large number of edges. Therefore, we again replace f STA with an approximate solution f STA computed using

L H . Define R(f ) 1 i=1 (f (x i ) -y(x i )) 2 as the empirical error and R(f ) 1 u u i=1 (f (x i ) -y(x i )) 2 as the generalization. Theorem 2. Let G be a fixed (connected) graph with eigen- values 0 = λ 1 (G) < λ 2 (G) ≤ . . . ≤ λ n (G)
, and H an ε-sparsifier of G. Let y ∈ R n be the labels of the nodes in G with |y(x)| ≤ c and F be the set of centered functions such that |f (x)y(x)| ≤ 2c. Let S ⊂ V be a random subset of labeled nodes, if the labels y S are centered, then w.p. at least 1δ (w.r.t. the random generation of the sparsifier H and the random subset of labeled points S) the resulting Stable-HFS solution satisfies

R( f ) ≤ R( f ) + β + 2β + 4c 2 ( + u) u π( , u) ln 1 δ 2 + 1 1 -ε 2(1 + ε)ε γλ 2 (G)c ((1 -ε) γλ 2 (G)-1) 2 2 , (6)
where f and f are computed on H and G,

π( , u) u + u -0.5 2 max{ , u} 2 max{ , u} -1 and β ≤ 3c √ ((1 -ε) γλ 2 (G) -1) 2 + 4c (1 -ε) γλ 2 (G) -1 • Thm. 2 (full proof in Appendix A)
shows how approximating G with H impacts the generalization error as the number of labeled samples increases. If we set ε = 0, we recover the bound of [START_REF] Cortes | Stability of transductive regression algorithms[END_REF], which depends only on R( f ) and β. When ε > 0, we see from Eq. 6 that the two terms already present in the exact case are either unchanged ( R( f )) or increase only by a constant factor β. Because of the approximation, a new error term (the last one in Eq. 6) is added to the bound, but we can see that it is negligible compared to β. In fact, it converges to zero as O(ε 2 / 2 (1ε) 4 ) as grows and it is dominated by β for any constant value of ε. This means that increasing ε corresponds to a constant increase in the bound, regardless of the size of the problem. Consequently, ε can be freely chosen to trade off accuracy and space complexity (Thm. 1) depending on the problem constraints. Finally, because the eigenvalues present in the bound are the ones of the original graph, any additional knowledge on the spectral properties of the input graph can be easily included in the analysis. Therefore, it is straightforward to provide stronger guarantees for Sparse-HFS when combined with assumptions on the graph generating model. Finally, we remark the level of generality of this result that holds for the integration between HFS and any ε-accurate spectral sparsification method. We postpone computational considerations to the following subsection.

Generalization bounds for LAPSMO

Starting from the closed form solution of LAPSMO (Eq. 1) we can replace the L G matrix with a sparsified Laplacian L H and using a fast linear solver, compute an approximate solution f = (λL

H + I) -1 y in O(n log 2 (n)) time and O(n log(n)) space. Finally, we can decompose the error as f -f 2 2 ≤ f -f 2 2 + f -f 2 2 .
The first term can be bounded using classical results from empirical process theory [START_REF] Bühlmann | Statistics for highdimensional data: methods, theory and applications[END_REF]. We bound the second term in the following theorem.

Theorem 3. For an arbitrary graph G and its (ε, γ)sparsifier, let f be the LAPSMO solution computed using L G and f the solution computed using L H . Then,

f -f 2 2 ≤ ε 2 1 -ε (0.25 + λγ) λ f T L G f + λγ f 2 2 ,
where λ is the regularization of LAPSMO.

For ε-sparsifiers, [START_REF] Sadhanala | Graph sparsification approaches for laplacian smoothing[END_REF] derive a similar bound f -

f 2 2 ≤ O(λ f T L G f ).
Setting γ = 0, we recover their bound up to constants. When γ > 0 instead, additional error terms emerge due to the introduced bias. In particular, the term λγ f 2 2 depends on the norm of the exact solution f , which in turn depends on the value of λ. Nonetheless, when f 2 2 is small, as is the case in our experiments, setting γ = 1/λ makes this term a constant, which is reflected by the good empirical performance. Computationally, for both Stable-HFS and LAPSMO, passing from computing a solution on the full graph to computing a solution on the sparsifier reduces the number of edges, which makes the memory and runtime plummet. Moreover, carefully distributing the sparsification process across multiple machines allows computing a final solution in a time independent from the number of edges, since the preprocessing sparsification step takes only O(n log 3 (n)) time, and the solution step only O(n log 2 (n)). Up to logarithmic terms, this results in an overall O(n) near-linear runtime, without any assumptions on the input graph. For graphs with a particularly favorable spectrum and problems with enough regularization, this is only O(d eff (γ)), resulting in a potentially sub-linear runtime. This result, only possible due to a particular structure of learning problems, opens up unexplored possibilities that would not be possible for general graph problems.

Bounds for other problems

Many other problems can be well approximated using (ε, γ)sparsifiers. For example, the cost of a SC solution evaluated on L H is very close to the cost evaluated on L G .

Proposition 2. For any rank k orthogonal projection

F T F, if H is an (ε, γ)-sparsifier of G, we have Tr(F T L H F) ≤ (1 + ε) Tr(F T L G F) + εγk.
Therefore, a clustering that well separates the sparsifier will also separate well the true graph. Similarly, we can obtain strong approximation guarantees for a variety of other Laplacian-based algorithms. Regularized problems such as LTR [START_REF] Cortes | Stability of transductive regression algorithms[END_REF], Laplacian-regularized least squares, and Laplacian SVM [START_REF] Belkin | On manifold regularization[END_REF] are of particular interest since the additive γ error is absorbed by the regularization and it is possible to provide strong generalization guarantees.

Experiments

We empirically validate our theoretical findings by testing how (ε, γ)-sparsifiers improves computational complexity without sacrificing final accuracy.

Dataset. We run experiments on the Amazon co-purchase graph [START_REF] Sadhanala | Graph sparsification approaches for laplacian smoothing[END_REF]. This graph fits our setting: It cannot be generated from vectorial data and is only artificially sparse, since the crawler that created it had no access to the true private co-purchase network held by Amazon. To compensate, [START_REF] Gleich | Using local spectral methods to robustify graph-based learning algorithms[END_REF] use a densification procedure that given the graph adjacency matrix

A G , com- putes all k-step neighbors A G,k k s=1 A s G .
We make the graph unweighted for numerical stability. The final graph has n = 334, 863 nodes and m = 98, 465, 352 edges, with an average degree of 294. We followed an approach similar to [START_REF] Sadhanala | Graph sparsification approaches for laplacian smoothing[END_REF] and introduce a hand-designed smooth signal as a target. We then perform 2000 iterations of the power method to compute an approximation of the smallest eigenvector v min , which is used as a smooth function over the graph.

Baselines. For all setups, we compute an "exact" solution (up to convergence error) using a fast linear solver. Computing this EXACT baseline requires O(m log(n)) time and O(m) space and achieves the best performance. Afterwards, we compare three different sparsification procedures to eval- Experimental procedure. We repeat each experiment 10 times with different sparsifiers and report the average performance of f on the specific task and its standard deviation. More details on experiments are given in the Appendix C.

Laplacian smoothing with Gaussian noise

We set f = v min and test different levels of noise, log 10 (σ) ∈ {-3, -2, -1, 0}. After constructing the sparsifier H, we compute an approximate solution f using LAPSMO (Eq. 2) with λ ∈ {10 -3 , 10 -2 , 10 -1 , 1, 10}. We measure the performance by the squared error

D( f ) = f -f 2 2 . As f 2 2 = v min 2 2 = 1, good values of D( f ) should be below 1.
Accuracy. In the interest of space, in Tab. 1, we report results for σ = {0.001, 0.01} and the best regularization λ for each method. We first notice that all sparsifiers are considerably smaller than the original graph, keeping only a small fraction of its edges. The smallest sparsifiers are obtained by DiSRe when γ is large. The comparison with DiSRe with γ = 0 (i.e., ε-sparsifier) confirms that the additive error translates into an extra compression of the resulting sparsifier. This also impacts the accuracy which degrades as γ increases. Nonetheless, we notice that while ε-sparsifiers perfectly match the accuracy of the exact method, even for large γ (and thus much smaller graph), DiSRe still outperforms kN, which has a significantly worse accuracy. Finally, we note that for γ = 0, the impact of q is as expected: Increasing q increases the size of the sparsifier and slightly improves the performance.

Computational complexity. All algorithms require 90s to load the graph from disk. The preprocessing phase of kN takes slightly less than 1min, while DiSRe's takes 12min on 4 machines. For the solving step, EXACT is unsurprisingly the slowest, requiring 12min to compute an f solution. Both kN and (ε, γ)-sparsifiers require 1-2min, depending on the number of edges preserved. Overall, preprocessing the graph with DiSRe before computing a solution does not introduce any overhead compared to EXACT (both take roughly 12min). We notice that while kN is overall faster, the time for DiSRe could be easily reduced by increasing the number of parallel processes when computing effective resistances or with a better network topology allowing pointto-point communication. Moreover, once we have access to an accurate ε-sparsifier, it is easier to solve problem repeatedly, e.g., to cross-validate regularization. For example, computing a solution for 4 different values of λ (see the appendix) is crucial for good performance and requires 48min for EXACT and only 20min for DiSRe. Finally, memory usage is reduced by a factor of 3 as EXACT requires over 30GB of memory to execute while DiSRe never exceeds 10GB. We expect these advantages to only grow larger as we scale to larger graphs.

SSL with harmonic function solution

We also test DiSRe on a SSL problem. The labels are generated taking the sign of f = v min and ∈ {20, 346, 672, 1000} labels are revealed. The labeled nodes are chosen at random so that 0 and 1 labels are balanced in the dataset. We run Stable-HFS with λ ∈ {10 -6 , 10 -4 , 10 -2 , 1}. In Tab. 1, we report results for = {346, 672} and the best λ for each method. We run DiSRe with γ = 0 as Stable-HFS does not have any regularization and ε-sparsifiers are preferable. The average size of the sparsifiers is the same as before as they are agnostic to the learning task. Similar to the smoothing case, DiSRe achieves a performance that closely approximates the exact solution, despite the significant compression of the original graph. Furthermore, the effectiveness of the ε-sparsifier returned by DiSRe is confirmed by its comparison with kN, whose error is significantly worse. Finally, we notice that the computational analysis in the previous section holds for SSL as well. In fact, although the learning task is different, we use the same SSD solver to compute the HFS and thus the running time are comparable in the two tasks.

In the proofs, we make several uses the following reformulation of Def. 1.

Proposition 3. A sub-graph H is a (ε, γ)-sparsifier of G if and only if (1 -ε)L G -εγI L H (1 + ε)L G + εγI ⇐⇒ (L G + γI) -1/2 (L H -L G )(L G + γI) -1/2 2 2 ≤ ε.
Theorem 2. Let G be a fixed (connected) graph with eigenvalues 0 = λ 1 (G) < λ 2 (G) ≤ . . . ≤ λ n (G), and H an ε-sparsifier of G. Let y ∈ R n be the labels of the nodes in G with |y(x)| ≤ c and F be the set of centered functions such that |f (x)y(x)| ≤ 2c. Let S ⊂ V be a random subset of labeled nodes, if the labels y S are centered, then w.p. at least 1δ (w.r.t. the random generation of the sparsifier H and the random subset of labeled points S) the resulting Stable-HFS solution satisfies

R( f ) ≤ R( f ) + β + 2β + 4c 2 ( + u) u π( , u) ln 1 δ 2 + 1 1 -ε 2(1 + ε)ε γλ 2 (G)c ((1 -ε) γλ 2 (G)-1) 2 2 , ( 6
)
where f and f are computed on H and G, π( , u) u + u -0.5 2 max{ , u} 2 max{ , u} -1 and

β ≤ 3c √ ((1 -ε) γλ 2 (G) -1) 2 + 4c (1 -ε) γλ 2 (G) -1 • Proof of Thm. 2.
Step 1 (generalization of stable algorithms). Let β be the stability of Stable-HFS when using the sparsified Laplacian L H in place of L G . Then using the result of [START_REF] Cortes | Stability of transductive regression algorithms[END_REF], we have that with probability at least 1δ (w.r.t. the randomness of the labeled set S) the solution f satisfies

R( f ) ≤ R( f ) + β + 2β + c 2 ( + u) u π( , u) log(1/δ) 2 •
In order to obtain the final result, we first derive an upper bound on the stability β and relate the empirical error of f STA to the one of f STA . Furthermore, it can be shown that if we center the vector of labels y S y Sy S , with y 1 y T S 1, then the solution of Stable-HFS can be rewritten in closed form as f STA = (γ L G + I S ) + ( y S -µ1) = (P(γlL G + I S ))

+ y S .

Step 2 (stability). The bound on the stability follows similar steps as in the analysis of Stable-HFS by [START_REF] Belkin | Regularization and Semi-Supervised Learning on Large Graphs[END_REF] integrated with the properties of spectral sparsifiers reported in Def. 1. Let S and S be two labeled sets only that only differ by one element and f and f be the solutions obtained by running Stable-HFS using L H and S and S respectively.

Without loss of generality, we assume that I S ( , ) = 1 and I S ( + 1, + 1) = 0, and the opposite for I S . The original proof of [START_REF] Cortes | Stability of transductive regression algorithms[END_REF] showed that the stability β can be bounded as β ≤ ff . In the following, we show that the difference between the solutions f and f and thus the stability of the algorithm, is strictly related to eigenvalues of the sparse graph H. Let A P( γL H + I S ) and B P( γL H + I S ). We remind that if the labels are centered, the solutions of Stable-HFS can be conveniently written as f = A -1 y S and f = B -1 y S . As a result, the difference between the solutions can be written as

f -f = A -1 y S -B -1 y S ≤ A -1 ( y S -y S ) + A -1 y S -B -1 y S . (7)
Let us consider any vector f ∈ F, since the null space of a Laplacian L H is the one vector 1 and P = L H L + H , then Pf = f . Thus, we have

P( γL H + I S )f (1) ≥ P γL H f -PI S f (2) ≥ PlγL H f -f (3) ≥ ( γλ 2 (H)-1) f , (8)
where (1) follows from the triangle inequality and (2) follows from the fact that PI S f ≤ f , since the largest eigenvalue of the projection matrix P is one and the norm of f restricted on S is smaller than the norm of f . Finally, (3) follows from the fact that

PL H f = L H L + H L H f = L H f and since f is orthogonal to the null space of L H then L H f ≥ λ 2 (H) f
, where λ 2 (H) is the smallest non-zero eigenvalue of L H . At this point we can exploit the spectral guarantees of the sparsified Laplacian L H and we have that λ 2 (H) ≥ (1ε)λ 2 (G). As a result, we have an upper bound on the spectral radius of the inverse operator (P( γL H + I S )) -1 and thus

A -1 (y S -y S ) ≤ 4M γ(1 -ε)λ 2 (G) -1 ,
where the first step follows from Eq. 8 since both y S and y S are centered and thus (y Sy S ) ∈ F and the second step is obtained by bounding y Sy S ≤ y Sy S + y Sy S ≤ 4M . The second term in Eq. 7 can be bounded as

A -1 y S -B -1 y S = B -1 (B -A)A -1 y S = B -1 P(I S -I S )A -1 y S ≤ 1.5M √ ( γ(1 -ε)λ 2 (G) -1) 2 ,
where we used y S ≤ y S + y S ≤ 2M √ , P(I S -I S ) ≤ √ 2 < 1.5 and we applied Eq. 8 twice. Putting it all together we obtain the stated bound.

Step 3 (empirical error). The other element effected by the sparsification is the empirical error R( f ). We first recall that

P = L + G L G = L -1/2 G L G L -1/2 G
(and equivalently with G replaced by H) and we introduce

P = L -1/2 G L H L -1/2 G Let A P( γL H + I S ), A P( γL G + I S )
, then we can rewrite the empirical error as

R( f ) = 1 I S f -I S f + I S f -y S 2 ≤ 1 I S f -y S 2 + 1 I S f -I S f 2 ≤ R( f ) + 1 I S ( A -1 -A -1 ) y S 2 ≤ R( f ) + 1 A -1 ( A -A) A -1 y S 2 = R( f ) + 2 γ 2 A -1 (P(L G -L H )) A -1 y S 2 = R( f ) + γ 2 A -1 (P(L G -L H )P) A -1 y S 2 ,
where in the last passagen we use L G P = L G and L H P = L H . To bound the second term, we derive

A -1 P(L G -L H )P A -1 y S 2 (1) = A -1 L 1/2 G L -1/2 G (L G -L H )L -1/2 G L 1/2 G A -1 y S 2 (2) = A -1 L 1/2 G PL -1/2 G (L G -L H )L -1/2 G PL 1/2 G A -1 y S 2 (3) = A -1 L 1/2 G P(P -P)PL 1/2 G A -1 y S 2 ≤ A -1 L 1/2 G P 2 P -P 2 PL 1/2 G A -1 2 y S 2 (4) ≤ M 2 ε 2 A -1 L 1/2 G P 2 A -1 L 1/2 G P 2 ,
where in (1) and (2), we use the definition of P, in (3) we use the definition of P, while in (4) we use the fact that Def. 1 implies that (1ε)P P (1 + ε)P and thus the largest eigenvalue of P -P is ε and P -P 2 ≤ ε 2 . We need now to bound

A -1 L 1/2 G P 2 = A -1 L 1/2 G 2
. From the definition of spectral norm,

A -1 L 1/2 G 2 = max x =1 x T L 1/2 G A -1 A -1 L 1/2 G x = max x =1 x T A -1 L G A -1 x ≤ 1 1 -ε max x =1 x T A -1 L H A -1 x = 1 1 -ε A -1 L 1/2 H 2 = 1 1 -ε A -1 L 1/2 H P 2 .
Similarly to Eq. 8, finding a lower bound on AL

-1/2 H
Px for all x is equivalent to find a lower bound for all f ∈ F to

P( γL H + I S )L -1/2 H f ≥ P γL H L -1/2 H f -PI S L -1/2 H f ≥ P γL H L -1/2 H f -L -1/2 H f ≥ γ λ 2 (H)- 1 λ 2 (H) f ≥ 1 λ 2 (H) ( γλ 2 (H)-1) f ≥ 1 (1 + ε)λ 2 (G) ( γ(1 -ε)λ 2 (G)-1) f .
Similarly, we can show that

P(lγL G + I S )L -1/2 G f ≥ P γL G L -1/2 G f -PI S L -1/2 G f ≥ P γL G L -1/2 G f -L -1/2 G f ≥ γ λ 2 (G) λ 2 (G) - 1 λ 2 (G) f ≥ 1 λ 2 (G) ( γλ 2 (G)-1) f ≥ 1 (1 + ε)λ 2 (G) ( γ(1 -ε)λ 2 (G)-1) f .
Taking this and putting all together gives

R( f ) ≤ R( f ) + 1 1 -ε (1 + ε) 2 ε 2 2 γ 2 λ 2 (G) 2 M 2 ( γ(1 -ε)λ 2 (G)-1) 4 •
Combining the three steps above concludes the proof.

Theorem 3. For an arbitrary graph G and its (ε, γ)-sparsifier, let f be the LAPSMO solution computed using L G and f the solution computed using L H . Then,

f -f 2 2 ≤ ε 2 1 -ε (0.25 + λγ) λ f T L G f + λγ f 2 2 ,
where λ is the regularization of LAPSMO.

Proof. We need to bound the distance between f and f . Using the definition, the fact that

A -1 -B -1 = B -1 (B -A)A -1 and collecting (L G + γI) 1/2 we have f -f 2 2 = (λL H + I) -1 -λL G + I) -1 )y 2 2 = (λL H + I) -1 (λL H -λL G )(λL G + I) -1 )y 2 2 = λ 2 (λL H + I) -1 (λL H -L G )(λL G + I) -1 )y 2 2 = λ 2 (λL H + I) -1 (L G + γI) 1/2 (L G + γI) -1/2 (L H -L G )(L G + γI) -1/2 (L G + γI) 1/2 f 2 2 .
Then, using Prop. 3,

λ 2 (λL H + I) -1 (L G +γI) 1/2 (L G + γI) -1/2 (L H -L G )(L G + γI) -1/2 (L G + γI) 1/2 f 2 2 ≤ ε 2 λ 2 (λL H + I) -1 (L G + γI) 1/2 2 2 f T (L G + γI) f = ε 2 1 -ε λ 2 (λL H + I) -1 ((1 -ε)L G -εγI + γI)(λL H + I) -1 2 f T (L G + γI) f ≤ ε 2 1 -ε λ 2 (λL H + I) -1 (L H + γI)(λL H + I) -1 2 f T (L G + γI) f = ε 2 1 -ε λ 2 max i λ i (L H ) + γ (λλ i (L H ) + 1) 2 f T (L G + γI) f = ε 2 1 -ε λ max i λλ i (L H ) (λλ i (L H ) + 1) 2 + λγ (λλ i (L H ) + 1) 2 f T (L G + γI) f ≤ ε 2 1 -ε λ max i λλ i (L H ) (λλ i (L H ) + 1) 2 + max i λγ (λλ i (L H ) + 1) 2 f T (L G + γI) f ≤ ε 2 1 -ε λ (0.25 + λγ) f T (L G + γI) f ≤ ε 2 1 -ε ((0.25 + λγ) λ f T L G f + (0.25 + λγ) λγI),
which concludes the proof.

B. Proof of Thm. 1

Theorem 1. Let ε > 0 be the accuracy, 0 ≤ δ ≤ 1 the probability of error, and ρ (1 + 3ε)/(1ε). Given an arbitrary graph G and an arbitrary merge tree structure, if DiSRe is run with parameter q 26ρ log(3n/δ)/ε 2 , then each sub-graphs H {h, } is an (ε, γ)-sparsifier of G {h, } with at most 3qd eff (γ) edges with probability 1δ. Whenever the merge tree is balanced and k is big enough such that m/k ≤ 3qd eff (γ),6 then merge operations can be run in parallel across the machines with an overall time complexity of O(d eff (γ) log 3 (n)), a total work O(m log 3 (n)), and O(log(n)) rounds of communication.

The proof of Thm. 1 assembles the techniques from two sources. First, it is based on the analysis of the algorithm of [START_REF] Kelner | Spectral sparsification in the semi-streaming setting[END_REF], that shows it produces a spectral sparsifier in high probability. We have previously published this analysis as a technical report [START_REF] Calandriello | Analysis of kelner and levin graph sparsification algorithm for a streaming setting[END_REF] for the case of γ = 0. In this alternative proof, we rigorously take into account the dependencies across subsequent resparsifications using martingale inequalities, fixing a flaw in the original analysis. Second, it also copies the steps for the analysis of SQUEAK [START_REF] Calandriello | Distributed sequential sampling for kernel matrix approximation[END_REF] which is an algorithm for general kernel sparsificatioN. In particular, Alg. 1 is an instantiation of SQUEAK to the special case of graph sparsification. While SQUEAK's analysis by [START_REF] Calandriello | Distributed sequential sampling for kernel matrix approximation[END_REF] holds in general, in DiSRe we can exploit the specific structure of graph Laplacians to perform a few optimizations. For completeness and ease of verification, we restate here SQUEAK's original proof with the necessary modifications to the notation, constants, and relevant quantities to target the graph learning setting.

Merge trees

We first formalize the random process induced by Alg. 1.

We partition

G into k disjoint sub-graphs G i of size n i , such that G = ∪ k e=1 G i .
For each sub-graph G i , we construct an initial sparsifier H {1,i} {(j, p 0,i = 1, q 0,i = q) : j ∈ G i } by inserting all edges from G i into H 1,i with weight p 0,i 1 and number of copies q 0,i q. It is easy to see that H {1,i} is an (0, 0)-accurate sparsifier, and we can split the graph into small enough sub-graphs to make sure that it can be easily stored and manipulated in memory. Afterwards, the initial sparsifiers H {1,i} are included into the sparsifier pool S 1 .

At iteration h, the inner loop of Alg. 1 arbitrarily chooses two sparsifiers from S h and merges them into a new sparsifier. Any arbitrary sequence of merges can be described by a full binary tree, i.e., a binary tree where each node is either a leaf or has exactly two children. Fig. 2 shows several different merge trees corresponding to different choices for the order of the merges. Note that starting from k leaves, a full binary tree will always have exactly k -1 internal nodes. Therefore, regardless of the structure of the merge tree, we can always transform it into a tree of depth k, with all the initial sparsifiers H 1,i as leaves on its deepest layer. After this transformation, we index the tree nodes using their height (longest path from the node to a leaf, also defined as the depth of the tree minus the depth of the node), where leaves have height 1 and the root has height k. We can also see that at each layer, there is a single sparsifier merge, and the size of S h (number of sparsifiers

H 1,2,3,4,5 H 1,2,3 + H 4 H 1,2,3 H 1 + H 2,3 H 1 H 1 G 1 H 2,3 H 2 + H 3 H 2 G 2 H 3 G 3 H 4 H 4 H 4 G 4 h = 4 h = 3 h = 2 h = 1 (a) arbitrary tree H 1,2,3,4 H 1,2,3 + H 4 H 1,2,3 H 1,2 + H 3 H 1,2 H 1 + H 2 H 1 G 1 H 2 G 2 H 3 H 3 G 3 H 4 H 4 H 4 G 4 h 4 h = 3 h = 2 h = 1 (b) sequential tree H 1,2,3,4 H 1,2 + H 3,4 H 1,2 H 1,2 H 1 + H 2 H 1 G 1 H 2 G 2 H 3,4 H 3 + H 4 H 3 H 3 G 3 H 4 H 4 G 5 h 4 h = 3 h = 2 h = 1 (c) minimum depth tree Figure 2. trees for Alg. 1. present at layer h) is |S h | = k -h + 1.
Therefore, a node corresponding to a sparsifier is uniquely identified with two indices {h, }, where h is the height of the layer and ≤ |S h | is the index of the node in the layer. For example, in Fig. 2(a), the node containing H 1,2,3 is indexed as {3, 1}, and the highest node containing H 4 is indexed as {3, 2}.

We also define the graph G {h, } as the union of all sub-graph G that are reachable from node {h, } as leaves. For example, in Fig. 2(a), sparsifier H 1,2,3 in node {3, 1} is constructed starting from all edges in G {3,1} = G 1 ∪ G 2 ∪ G 3 , where we highlight in red the descendant tree. We now define L h as the block diagonal matrix where each diagonal block L G {h, } is the Laplacian constructed on G {h, } . Without loss of generality, we will assume that each of the sub-graphs G i is connected and spans all the n nodes in the graph. This simplifies the notation for the L G {h, } matrices, making them all n × n matrices. If this is not the case, the whole proof still follows through with different number of nodes n {h, } for each G {h, } . Again, from Fig. 2, L 3 is a 2n × 2n matrix with two blocks on the diagonal, a first n × n block L G {3,1} constructed on G {3,1} = G 1 ∪ G 2 ∪ G 3 , and a second n × n block L G {3,2} constructed on G {3,2} = G 4 . Similarly, we combine Def. 1 and Prop. 3 to define

P {h, } (L G {h, } + γI) -1/2 L G {h, } (L G {h, } + γI) -1/2 and P {h, } (L G {h, } + γI) -1/2 L H {h, } (L G {h, } + γI) -1/2 ,
and have P h as a block diagonal projection matrix, where each block P {h, } is defined using L G {h, } , and block diagonal P h , where each block P {h, } is defined using L H {h, } and H {h, } .

The statement. Since P h -P h is block diagonal, we have that a bound on its largest eigenvalue implies an equal bound on each matrix on the diagonal, i.e., P h -P h = max P {h, } -P {h, } ≤ ε =⇒ P {h, } -P {h, } ≤ ε for all blocks on the diagonal, and since each block corresponds to a sparsifier H {h, } , this means that if P h -P h ≤ ε, all sparsifiers at layer are (ε, γ)-sparsifiers of their respective graphs. Let d {h, } eff (γ) be the effective dimension of L G {h, } .

Our goal is to show

P ∃h ∈ {1, . . . , k} : P h -P h 2 ≥ ε ∪ max =1,...,|S h | |H {h, } | ≥ 3qd {h, } eff (γ) = P ∃h ∈ {1, . . . , k} : max =1,...,|S h | P {h, } -P {h, } 2 ≥ ε A h ∪ max =1,...,|S h | |H {h, } | ≥ 3qd {h, } eff (γ) B h ≤ δ, (9)
where event A h refers to the case when some sparsifier H {h, } at an intermediate layer h fails to accurately approximate L {h, } and event B h considers the case when the memory requirement is not met (i.e., too many edges are kept in one of the sparsifiers H {h, } at a certain layer h). After reformulating and a union bound we obtain

P ∃h ∈ {1, . . . , k} : P h -P h 2 ≥ ε ∪ max =1,...,|S h | |H {h, } | ≥ 3qd {h, } eff (γ) ≤ k h=1 |S h | =1 P P {h, } -P {h, } 2 ≥ ε + k h=1 |S h | =1 P |H {h, } | ≥ 3qd {h, } eff (γ) ∩ ∀h ∈ {1, . . . , h} : P h -P h 2 ≤ ε ≤ δ. ( 10
)
The accuracy of the sparsifier (first term in the previous bound) is guaranteed by the fact that given an (ε, γ)-accurate sparsifier we obtain γ-effective resistance estimates (i.e. RLS estimates) which are at least a fraction of the true ones, thus forcing the algorithm to sample each column enough. On the other hand, the space complexity bound is achieved by exploiting the fact that estimates are always upper-bounded by the true γ-effective resistance, thus ensuring that Alg. 1 does not oversample columns w.r.t. the sampling process following the exact γ-effective resistance.

In the reminder of the proof, we will show that both events happen with probability smaller than δ/(2k 2 ).

Since |S h | = k -h + 1, we have k h=1 |S h | =1 δ 2k 2 = k h=1 (k -h + 1) δ 2k 2 = k(k + 1) δ 4k 2 ≤ k 2 δ 2k 2 = δ 2 ,
and the union bound over all events is smaller than δ. The main advantage of splitting the failure probability as we did in Eq. 10 is we can now analyze the processes that generated each P {h, } -P {h, } (and each sparsifier H {h, } ) separately.

Focusing on a single node {h, } restricts our problem on a well defined graph G {h, } , where we can analyze the evolution of H {h, } sequentially.

B.1. Bounding the projection error P {h, } -P {h, }

The sequential process. Thanks to the union bound in Eq. 10, instead of having to consider the whole merge tree followed by Alg. 1, we can focus on each individual node {h, } and study the sequential process that generated its sparsifier H {h, } .

We will now map more clearly the actions taken by Alg. 1 to the process that generated P {h, } -P {h, } . We begin by focusing on P {h, } , which is a random matrix defined starting from the fixed graph Laplacian L G {h, } and the random sparsifier Laplacian L H {h, } , where the randomness influences both which edges are included in H {h, } , and the weight with which they are added.

Note that since the merge tree is decided in advance, the graph G {h, } is not a random object and is fixed for the whole process. Consider now an edge e ∈ G {h, } . Again for simplicity and without loss of generality 7 , we will assume that the starting graphs in the leaves are edge-disjoint. Therefore, there is a single path in the tree, with length h, from the leaves to {h, }. This means that for all s < h, we can properly define a unique p s,e and q s,e associated with that point. More in detail, if at layer s point i is present in G {s, } , it means that either (1) Alg. 1 used H {s, } to compute p s,e , and p s,e to compute q s,e , or (2) at layer h, Alg. 1 did not have any merge scheduled for point i, and we simply propagate p s,e = p s-1,i and q s,e = q s-1,i . Consistently with the algorithm, we initialize p 0,i = 1 and q 0,i = q.

Denote m {h, } = |G {h, } | so that we can use index

i ∈ [m {h, } ] to index all edges in G {h, } . Given the n × m {h, } matrix Q = (L G {h, } + γI) -1/2 B G {h, }
with its e-th column q i = (L G {h, } + γI) -1/2 B G {h, } e m {h, } ,e , we can rewrite the projection matrix as that P {h, } = QQ T = m {h, } e=1

q e q T e . Note that q e q T e = q T e q e = e T m {h, } ,e Q T Qe m {h, } ,e = e T m {h, } ,e QQ T e m {h, } ,e = e T m {h, } ,e P {h, } e m {h, } ,e = r e (γ),

or, in other words, the norm q e q T e is equal to the γ-effective resistance of the e-th edge w.r.t. to graph G {h, } . Note that since e is present only in node l on layer h, its γ-effective resistance is uniquely defined w.r.t. G {h, } and can be shortened as r h,e . Using q e , we can also introduce the random matrix P {h, } s as

P {h, } s = m {h, } e=1
q s,e q p s,e q e q T e = m {h, } e=1 q j=1 z s,e,j q p s,e q e q T e .

where z s,e,j are {0, 1} r.v. such that q s,e = q j=1 z s,e,j , or in other words z s,e,j are the Bernoulli random variables that compose the Binomial q s,e associated with edge e, with j indexing each individual copy of the edge. Note that when s = h, we have that P {h, } h = P {h, } and we recover the definition of the approximate projection matrix from Alg. 1. But, for a general s = h P {h, } s does not have a direct interpretation in the context of Alg. 1. It combines the vectors q e , which are defined using L G {h, } at layer h, with the weights p s,e computed by Alg. 1 across multiple nodes at layer s, which are potentially stored in different machines that cannot communicate. Nonetheless, P {h, } s is a useful tool to analyze Alg. 1.

Taking into account that we are now considering a specific node {h, }, we can drop the index from the graphs G {h, } = G, γ-effective resistances τ h,e , and size m {h, } = m. Using this shorter notation, we can reformulate our objective as bounding P {h, } -P {h, } 2 = P {h, } -P {h, } h 2 , and reformulate the process as a sequence of matrices {Y s } h s=1 defined as

Y s = P {h, } -P {h, } s = 1 q m e=1 q j=1
1 -z s,e,j p s,e q e q T e ,

where Y h = P {h, } -P {h, } h = P {h, } -P {h, } , and Y 1 = P {h, } -P {h, } 0 = 0 since p 0,i = 1 and q 0,i = q.

B.2. Bounding Y h

We transformed the problem of bounding P {h, } -P {h, } into the problem of bounding Y h , which we modeled as a random matrix process, connected to Alg. 1 by the fact that both algorithm and random process Y h make use of the same weight p s,e and multiplicities q s,e .

The frozen process. Inspired by [START_REF] Cohen | Online row sampling[END_REF], we will now replace the process Y s with an alternative process Y s defined as

Y s = Y s-1 I Y s-1 ≤ ε + Y s-1 I Y s-1 ≥ ε .
This process starts from Y 0 = Y 0 = 0, and is identical to Y s until a step s where for the first time Y s ≤ ε and Y s+1 ≥ ε. After this failure happen the process Y s is "frozen" at s and Y s = Y s+1 for all s + 1 ≤ s ≤ h. Consequently, if any of the intermediate elements of the sequence violates the condition Y s ≤ ε, the last element will violate it too. For the rest, Y s behaves exactly like Y s . Therefore,

P ( Y h ≥ ε) ≤ P Y h ≥ ε ,
and if we can bound P Y h ≥ ε we will have a bound for the failure probability of Alg. 1, even though after "freezing" the process Y h does not make the same choices as the algorithm.

We will see now how to construct the process Y s starting from z s,e,j and p s,e,j . We recursively define the indicator ({0, 1}) random variable z s,e,j as z s,e,j = I u s,e,j ≤ p s,e,j p s-1,e,j z s-1,e,j , where u s,e,j ∼ U(0, 1) is a [0, 1] uniform random variable and p s,e,j is defined as

p s,e,j = p s,e I Y s-1 ≤ ε ∩ z s-1,e,j = 1 + p s-1,e,j I Y s-1 ≥ ε ∪ z s-1,e,j = 0 .
This definition of the process satisfies the freezing condition, since if Y s+1 ≥ ε (we have a failure at step s), for all s ≥ s + 1 we have z s ,i,j = z s+1,i,j with probability 1 (p s+1,i,j /p s,i,j = p s,i,j /p s,i,j = 1), and the weights 1/(qp s+1,i,j ) = 1/(qp s,i,j ) never change.

Introducing a per-copy weight p s,e,j and enforcing that p s+1,i,j = p s,e,j when z s,e,j = 0 avoids subtle inconsistencies in the formulation. In particular, not doing so would semantically correspond to reweighting dropped copies. Although this does not directly affect Y s (since the ratio z s,e,j / p s,e is zero for dropped copies), and therefore the relationship

P ( Y h ≥ ε) ≤ P Y h ≥ ε
still holds, we will see later how maintaining consistency helps us bound the second moment of our process.

We can now arrange the indices s, e, and j into a linear index t = s in the range [1, . . . , m 2 q], obtained as t = {s, e, j} = (s -1)mq + (e -1)q + j. We also define the difference matrix as X {s,e,j} = 1 q z s-1,e,j p s-1,e,j -z s,e,j p s,e,j q e q T e , which allows writing the cumulative matrix as Y {s,e,j} =

{s,e,j} r=1

X {s,e,j} where the checkedges {s, m, q} correspond to Y s ,

Y {s,m,q} = Y s = 1 q m e=1 q j=1
1 -z s,e,j p s,e,j q e q T e .

Let F s be the filtration containing all the realizations of the uniform random variables u s,e,j up to the step s, that is F s = {u s ,e ,j , ∀{s , e , j } ≤ s}. Again, we notice that F s defines the state of the algorithm after completing iteration s because, unless a "freezing" happened, Alg. 1 and Y s flip coins with the same probability, and generate the same sparsifiers. Since r s,e and p s,e,j are computed at the beginning of iteration s using the sparsifier H {s, } (for some unique at layer s), they are fully determined by F s-1 . Furthermore, since F s-1 also defines the values of all indicator variables z s ,e,j up to z s-1,e,j for any i and j, we have that all the Bernoulli variables z s,e,j at iteration s are conditionally independent given F s-1 . In other words, we have that for any e , and j such that {s, 1, 1} ≤ {s, e , j } < s the following random variables are equal in distribution, z s,e,j F {s,e ,j } = z s,e,j F {s-1,m,q} ∼ B p s,e,j p s-1,e,j ,

and for any e , and j such that {s, 1, 1} ≤ {s, e , j } ≤ {s, m, q} and s = {s, e , j } we have the independence

z s,e,j F {s-1,m,q} ⊥ z s,e ,j F {s-1,m,q} . (12)
While knowing that Y s ≤ ε is not sufficient to provide guarantees for the approximate probabilities p s,e , we can show that it is enough to prove that the frozen probabilities p s,e,j are never too small.

Lemma 1. Let α (1 + 3ε)/(1ε) and p s,e,j be the sequence of probabilities generated by the freezing process. Then for any s, e, and j, we have p s,e,j ≥ p h,e /α = r h,e /α.

Proof. Let s be the step where the process freezes (s = h if it does not freeze), or, in other words, Y s < ε and Y s+1 ≥ ε. From the definition of p s,e,j , we have that p s,e,j ≥ p s,i = p s,e = max {min { r s,e , p s-1,e } , p s-1,e /2}

≥ min { r s,e , p s-1,e } = min { r s,e , p s-2,e } = min { r s,e , p s-3,e } . . . = min { r s,e , p 0,e } = r s,e , and therefore p s,e,j ≥ r s,e . Now let {s, } be the node where r s,e was computed. From Alg. 1 we know it is computed using the sparsifier H generated by the union of the two children of node {s, }, which is an (ε, 2γ)-sparsifier of G {s, } .

From its definition, we know that r s,e is computed by Alg. 1 as

r s,e = (1 -ε)φ T e (L H + (1 + ε)γI) -1 φ i ,
using only the edges in H that are available at node {s, }. Let once again G = G {h, } by dropping the subscript, and define H c as the complement set of all edges in other sparsifiers at level s not in H. From the definition of P {h, } s and Prop. 3 we know that

Y s = P {h, } -P {h, } s 2 = (L G + 2γI) -1/2 (L G -(L H + L H c ))(L G + 2γI) -1/2 2 ≤ ε
and we know that this implies

L H L G + ε(L G + 2γI) -L H c L G + ε(L G + 2γI).
Plugging it in the initial definition,

r s,e = (1 -ε)b T e (L H + (1 + ε)γI) -1 b i ≥ (1 -ε)b T e (L G + ε(L G + 2γI) + (1 + ε)γI) -1 b e = (1 -ε) 1 1 + 3ε b T e (B G B T G + γI) -1 b e ≥ 1 -ε 1 + 3ε r h,e ≥ r h,e α •
We now proceed by studying the process {Y s } h s=1 and showing that it is a bounded martingale. In order to show that Y s is a martingale, it is sufficient to verify the following (equivalent) conditions,

E Y s F s-1 = Y s-1 ⇔ E X {s,e,j} F s-1 = 0.
We begin by inspecting the conditional random variable X {s,e,j} |F s-1 . Given the definition of X {s,e,j} , the conditioning on F s-1 determines the values of z s-1,e,j and the approximate probabilities p s-1,e,j and p s,e,j . In fact, remember that these quantities are fully determined by the realizations in F s-1 which are contained in F s-1 . As a result, the only stochastic quantity in X {s,e,j} is the variable z s,e,j . Specifically, if Y s-1 ≥ ε, then we have p s,e,j = p s-1,e,j and z s,e,j = z s-1,e,j (the process is stopped), and the martingale requirement E X {s,e,j} F s-1 = 0 is trivially satisfied. On the other hand, if Y s-1 ≤ ε we have E us,e,j 1 q z s-1,e,j p s-1,e,j -z s,e,j p s,e,j q e q T e F s-1 = 1 q z s-1,e,j p s-1,e,j -z s-1,e,j p s,e,j E I u s,e,j ≤ p s,e,j p s-1,e,j F s-1 q e q T e = 1 q z s-1,e,j p s-1,e,j -z s-1,e,j p s,e,j p s,e,j p s-1,e,j q e q T e = 0,

where we use the recursive definition of z s,e,j and the fact that u s,e,j is a uniform random variable in [0, 1]. This proves that Y s is indeed a martingale. We now compute an upper bound R on the norm of the values of the difference process as X {s,e,j} = 1 q z s-1,e,j p s-1,e,j -z s,e,j p s,e,j q e q T e ≤ 1 q 1 p s,e,j q e q T e = 1 q

1 p s,e,j τ h,i ≤ 1 q α τ h,i τ h,i = α q R,
where we used Lem. 1 to bound p s,e,j ≤ r h,e /α. If instead, Y s-1 ≥ ε, the process is stopped and X s = 0 = 0 ≤ R.

We are now ready to use a Freedman matrix inequality from [START_REF] Tropp | Freedman's inequality for matrix martingales[END_REF] to bound the norm of Y.

Proposition 4 (Tropp, 2011, Theorem 1.2). Consider a matrix martingale {Y k : k = 0, 1, 2, . . . } whose values are self-adjoint matrices with dimension d, and let {X k : k = 1, 2, 3, . . . } be the difference sequence. Assume that the difference sequence is uniformly bounded in the sense that X k 2 ≤ R almost surely for k = 1, 2, 3, . . . .

Improved Large-Scale Graph Learning through Ridge Spectral Sparsification

Define the predictable quadratic variation process of the martingale as

W k k j=1 E X 2 j {X s } j-1 s=0 , for k = 1, 2, 3, . . . .
Then, for all ε ≥ 0 and σ 2 > 0,

P ∃k ≥ 0 : Y k 2 ≥ ε ∩ W k ≤ σ 2 ≤ 2d • exp - ε 2 /2 σ 2 + Rε/3 •
In order to use the previous inequality, we develop the probability of error for any fixed h as

P ( Y h ≥ ε) ≤ P Y h ≥ ε = P Y h ≥ ε ∩ W h ≤ σ 2 + P Y h ≥ ε ∩ W h ≥ σ 2 ≤ P Y h ≥ ε ∩ W h ≤ σ 2 (a) + P W h ≥ σ 2 (b)
.

Using the bound on X {s,e,j} 2 , we can directly apply Prop. 4 to bound (a) for any fixed σ 2 . To bound the part (b), we use the following lemma, proved later in Sec. B.3.

Lemma 2 (Low probability of the large norm of the predictable quadratic variation process). We have that

P W h ≥ 6α q ≤ n • exp - 2q α •
Since P {h, } is defined at most on n nodes, combining Prop. 4 with σ 2 = 6α/q, Lem 2, the fact that 2ε/3 ≤ 1 and the value used by Alg. 1, q = 39α log(2n/δ)/ε 2 we obtain

P P {h, } -P {h, } 2 ≥ ε = P ( Y h ≥ ε) ≤ P Y h ≥ ε ∩ W h ≤ σ 2 + P W h ≥ σ 2 ≤ 2n • exp - ε 2 q α 1 12 + 2ε/3 + n • exp - 2q α ≤ 3n • exp - ε 2 13α q = 3n • exp -3 log 2n δ = 3n • exp -log 2n δ 3 = 3nδ 3 8n 3 ≤ δ 2n 2 •
This, combined with the fact that k ≤ n 2 /n ≤ n since at most we can split our graph into n parts each containing n edges, concludes this part of the proof.

B.3. Proof of Lem. 2 (bound on predictable quadratic variation)

Step 1 (a preliminary bound). We start by writing out W t for the process Y s ,

W t = 1 q 2
{s,e,j}≤t E z s-1,e,j p s-1,e,j -z s,e,j p s,e,j 2 F {s,e,j}-1 q e q T e q e q T e .

We rewrite the expectation terms in the equation above as E z s-1,e,j p s-1,e,j -z s,e,j p s,e,j 2 F {s,e,j}-1

= E z 2 s-1,e,j p 2 s-1,e,j
-2 z s-1,e,j p s-1,e,j z s,e,j p s,e,j + z 2

s,e,j p 2

s,e,j F {s,e,j}-1

(a) = E z 2 s-1,e,j p 2 s-1,e,j
-2 z s-1,e,j p s-1,e,j z s,e,j p s,e,j + z 2

s,e,j p 2

s,e,j

F s-1 = z 2 s-1,e,j p 2 s-1,e,j -2 z s-1,e,j p s-1,e,j 1 p s,e,j E [z s,e,j | F s-1 ] + 1 p 2 s,e,j E z 2
s,e,j F s-1

(b) = z s-1,e,j p 2 s-1,e,j -2 z s-1,e,j p s-1,e,j z s-1,e,j p s-1,e,j + 1 p 2 s,e,j E [z s,e,j | F s-1 ] = 1 p 2 s,e,j E [z s,e,j | F s-1 ] - z s-1,e,j p 2 s-1,e,j (c)
= 1 p s,e,j z s-1,e,j p s-1,e,j - z s-1,e,j p 2 s-1,e,j = z s-1,e,j p s-1,e,j 1 p s,e,j - 1 p s-1,e,j ,
where in (a) we use the fact that the approximate probabilities p s-1,e,j and p s,e,j and z s-1,e,j are fixed at the end of the previous iteration, while in (b) and (c) we use the fact that z s,e,j is a Bernoulli of parameter p s,e,j /p s-1,e,j (whenever z s-1,e,j is equal to 1). Therefore, we can write W t at the end of the process as

W h = W {h,m,q} = 1 q 2 q j=1 m e=1 h s=1
z s-1,e,j p s-1,e,j 1 p s,e,j -1 p s-1,e,j q e q T e q e q T e .

We can now upper-bound W h as

W h 1 q 2 q j=1 m e=1 h s=1
z s-1,e,j p s-1,e,j 1 p s,e,j -1 p s-1,e,j q e q T e q e q T e = 1 q 2 q j=1 m e=1 z h,e,j p 2

h,e,j -z h,e,j p 2

h,e,j + h s=1 z s-1,e,j p s-1,e,j 1 p s,e,j -1 p s-1,e,j q e q T e q e q T e = 1 q 2 q j=1 m e=1 z h,e,j p 2

h,e,j

+ h s=1 - z s,e,j p 2
s,e,j + z s-1,e,j p s,e,j p s-1,e,j -z 0,e,j p 2 0,e,j q e q T e q e q T e 1 q 2 q j=1 m e=1 z h,e,j p 2

h,e,j + h s=1 z s-1,e,j p s,e,j p s-1,e,j -z s,e,j p s,e,j p s-1,e,j q e q T e q e q T e = 1 q 2 q j=1 m e=1 z h,e,j p 2

h,e,j + h s=1 z s-1,e,j (1z s,e,j ) p s,e,j p s-1,e,j q e q T e q e q T e ,

where in the inequality we use the fact p s,e,j ≤ p s-1,e,j . From the definition of p s,e,j , we know that when z s,e,j = 0, p s,e,j = p s-1,e,j . Therefore zs-1,e,j (1-zs,e,j ) p s,e,j p s-1,e,j = zs-1,e,j (1-zs,e,j )

p 2 s-1,e,j
, since the term is non-zero only when z s,e,j = 0. Finally, we see that only one of the z s-1,e,j (1z s,e,j ) terms can be active for s ∈ [h] and thus

W h 1 q 2 q j=1 m e=1 z h,e,j p 2
h,e,j + h s=1 z s-1,e,j (1z s,e,j ) p 2 s-1,e,j q e q T e q e q T e = 1 q 2 q j=1 m e=1 max max s=1,...,h z s-1,e,j (1z s,e,j ) p 2 s-1,e,j ; z h,e,j p 2 h,e,j q e q T e q e q T e = 1 q 2 q j=1 m e=1 q e q T e q e q T e max s=0,...,h z s,e,j p 2

s,e,j

• (13)

Step 2 (introduction of a stochastically dominant process). We want to study max s=0,...,h zs,e,j p 2

s,e,j

. To simplify notation, we will consider max s=0,...,h zs,e,j p s,e,j

, where we removed the square, which will be re-added in the end. We know trivially that this quantity is larger or equal than, 1 because z 0,e,j /p 0,e,j = 1, but upper-bounding this quantity is not trivial as the evolution of the various p s,e,j depends in complex way on the interaction between the random variables z s,e,j . Nonetheless, whenever p s,e,j is significantly smaller than p s-1,e,j , the probability of keeping a copy of edge e at iteration s (i.e., z s,e,j = 1) is also very small. As a result, we expect the ratio zs,e,j p s,e,j to be still small with high probability.

Unfortunately, due to the dependence between different copies of the edge at different iterations, it seems difficult to exploit this intuition directly to provide an overall high-probability bound on W h . For this reason, we simplify the analysis by replacing each of the (potentially dependent) chains {z s,e,j /p s,e,j } h s=0 with a set of (independent) random variables w 0,e,j that will stochastically dominate them.

We define the random variable w s,e,j using the following conditional distribution,8

P 1 w s,e,j ≤ a F s =      0 for a < 1/p s,e,j 1 -1 p s,e,j a
for 1/p s,e,j ≤ a < α/p h,e 1 for α/p h,e ≤ a .

To show that this distribution is well defined, we use Lem. 1 to guarantee that 1/p s,e,j ≤ a < α/p h,e . Note that the distribution of 1 ws,e,j conditioned on F s is determined by only p s,e,j , p h,e , and α, where p h,e and α are fixed. Remembering that p s,e,j is a function of F s-1 (computed using the previous iteration), we have that

P 1 w s,e,j ≤ a F s = P 1 w s,e,j ≤ a F s-1 .
Notice that in the definition of w s,e,j , none of the other w s ,e ,j (for any different s , e , or j ) appears and p s,e,j is a function of F s-1 . It follows that given F s-1 , w s,e,j is independent from all other w s ,e ,j (for any different s , e , or j ). This is easier to see in the probabilistic graphical model reported in Fig. 3, which illustrates the dependence between the various variables.

Finally for the special case w 0,e,j the definition above reduces to

P 1 w 0,e,j ≤ a =      0 for a < 1 1 -1 a for 1 ≤ a < α/p h,e 1 for α/p h,e ≤ a , (14)
since p 0,e,j = 1 by definition. From this definition, w 0,e,j and w 0,e ,j are all independent, and this will allow us to use stronger concentration inequalities for independent random variables.

F {s,e,j} = u s ,e ,j , ∀{s , e , j } ≤ {s, e, j} ∪ {ws,e,j} = F {s,e,j} ∪ {ws,e,j}.

Improved Large-Scale Graph Learning through Ridge Spectral Sparsification A grey filling indicates that a random variable is observed or a function of observed variables.

Step 3 Proving the dominance. We remind the reader that a random variable A stochastically dominates random variable B, if for all values a the two equivalent conditions are verified, P(A ≥ a) ≥ P(B ≥ a) ⇔ P(A ≤ a) ≤ P(B ≤ a).

As a consequence, if A dominates B, the following implication holds,

P(A ≥ a) ≥ P(B ≥ a) =⇒ E[A] ≥ E[B],
while the reverse (A dominates B, if

E[A] ≥ E[B]
) is not true in general. Following this definition of stochastic dominance, our goal is to prove

P h max s=0
z s,e,j p s,e,j ≤ a ≥ P 1 w 0,e,j ≤ a .

We prove this inequality by proceeding backwards with a sequence of conditional probabilities. We first study the distribution of the maximum conditional to the state of the algorithm at the end of iteration h, i.e., F h . From the definition of w h,e,j , we know that, w.p. F k-1 , we have P z k,e,j w k,e,j ≤ a F k-1 = P I u k,e,j ≤ p k,e,j p k-1,e,j

1 w k,e,j ≤ a F k-1 =            0 for a ≤ 0 1 - p k,e,j p k-1,e,j
for 0 ≤ a < 1/p k,e,j 1 -p k,e,j p k-1,e,j + p k,e,j p k-1,e,j 1 -1 p k,e,j a = 1 -

1 p k-1,e,j a
for 1/p k,e,j ≤ a < α/p h,e

1 for α/p h,e ≤ a ≥            0 for a < 1/p k-1,e,j 1 - 1 p k-1,e,j a
for 1/p k-1,e,j ≤ a < 1/p k,e,j 1 -

1 p k-1,e,j a
for 1/p k,e,j ≤ a < α/p h,e 1 for α/p h,e ≤ a (15)

= P 1 w k-1,e,j ≤ a F k-2 = P 1 w k-1,e,j ≤ a F k-1 ,
where the inequality is also represented in Fig. 4. We now proceed by peeling off layers from the end of the chain one by one, taking advantage of the dominance we just proved. Fig. 4 visualizes one step of the peeling when z k-1,e,j = 1 (note that the peeling is trivially true when z k-1,e,j = 0 since the whole chain terminated at step z k-1,e,j ). We show how to move from an iteration k ≤ h to k -1. where in (a), given F k-1 , everything is fixed except u k,e,j and w k,e,j and we can use the stochastic dominance in (15), and in (b) we use the fact that the inner maximum is always attained by 1/w k,e,j since by definition 1/w k-1,e,j is lower-bounded by 1/p k-1,e,j . Applying the inequality recursively from k = h to k = 1 removes all z s,e,j from the maximum and we are finally left with only w 0,e,j as we wanted, P max s=0,...,h z s,e,j p s,e,j ≤ a ≥ P max z 0,e,j p 0,e,j ; z 0,e,j w 0,e,j ≤ a ≥ P 1 w 0,e,j ≤ a ,

P

where in the last inequality we used that z 0,e,j = 1 from the definition of the algorithm and p 0,e,j = 1 while w 0,e,j ≤ 1 by ( 14).

Step 4 (stochastic dominance on W h ). Now that we proved the stochastic dominance of 1/w 0,e,j , we plug this result in the definition of W h . For the sake of notation, we introduce the term p max h ,e,j to indicate the maximum over the first h step of copy e, j such that max s=0,...,h z s,e,j p s,e,j = 1 p max h ,e,j

•

We first notice that while Y h is not necessarily PSD, W h is a sum of PSD matrices. Introducing the function Λ({1/p max h,e,j } e,j ) we can restate Eq. 13 as

W h = λ max (W h ) ≤ Λ({1/p max h,e,j } e,j ) λ max   1 q 2 q j=1 m e=1 1 p max h,e,j
2 q e q T e q e q T e   .

In Step 4, we showed that 1/p max h,e,j is stochastically dominated by 1/w 0,e,j for every e and j. In order to bound Λ({1/p max h,e,j } e,j ), we need to show that this dominance also applies to the summation over all columns inside the matrix norm. We can reformulate Λ({1/p max h,e,j } e,j ) as

λ max   1 q 2 q j=1 m e=1 1 p max h,e,j 2 q e q T e q e q T e   = max x: x =1 x T   1 q 2 q j=1 m e=1 1 p max h,e,j 2 q e q T e q e q T e   x = max x: x =1 1 q 2 q j=1 m e=1 1 p max h,e,j 2 q e 2 2 x T q e q T e x = max x: x =1 1 q 2 q j=1 m e=1 1 p max h,e,j 2 q e 2 q T e x 2 .
From this reformulation, it is easy to see that, because 1/p max h,e,j is strictly positive, the function Λ({1/p max h,e,j } e,j ) is monotonically increasing w.r.t. the individual 1/p max h,e,j , or in other words that increasing an 1/p max h,e,j without decreasing the others can only increase the maximum. Introducing Λ({1/w 0,e,j } e,j ) as Λ({1/w 0,e,j } e,j ) max

x: x =1 1 q 2 q j=1 m e=1 1 w 0,e,j 2 q e 2 q T e x 2 , we now need to prove the stochastic dominance of Λ({1/w 0,e,j } e,j ) over Λ({1/p max h,e,j } e,j ). Using the definition of 1/p max h,e,j , w h,e,j , and the monotonicity of Λ we have ; x also an increasing function. Therefore, we have

P   Λ   1 p max h,
P   Λ   max 1 p max k-1,e,j ; z k,e,j w k,e,j e,j   ≤ a   = E F k-1   P   Λ   max 1 p max k-1,e,j ; z k,e,j w k,e,j e,j   ≤ a F k-1     (a) ≥ E F k-1   P   Λ   max 1 p max k-1,e,j ; z k-1,e,j w k-1,e,j e,j   ≤ a F k-1     (b) = E F k-1   P   Λ   max 1 p max k-2,e,j ; z k-1,e,j w k-1,e,j e,j   ≤ a F k-1     ,
where inequality (a) follows from the fact that stochastic dominance is preserved by monotonically increasing functions [START_REF] Levy | Stochastic dominance: Investment decision making under uncertainty[END_REF], such as Λ, combined with the fact that for a fixed F k-1 the variables z k,e,j and w k,e,j are all independent and (b) from the definition of 1/p max k-1,e,j and the fact that by definition 1/w k-1,e,j is lower-bounded by 1/p k-1,e,j . We can iterate this inequality to obtain the desired result

P( W h ≥ σ 2 ) ≤ P   Λ   1 p max h,e,j e,j   ≥ σ 2   ≤ P   λ max   1 q 2 q j=1 m e=1
1 w 0,e,j 2 q e q T e q e q T e   ≥ σ 2   .

Step 5 (concentration inequality). Since all w 0,e,j are (unconditionally) independent from each other, we can apply the following theorem. Proposition 5 (Tropp, 2015, Theorem 5.1.1). Consider a finite sequence {X k : k = 1, 2, 3, . . . } whose values are independent, random, PSD Hermitian matrices with dimension d. Assume that each term in the sequence is uniformly bounded in the sense that λ max (X k ) ≤ L almost surely for k = 1, 2, 3, . . . .

Introduce the random matrix V k X k , and the maximum eigenvalue of its expectation

µ max λ max (E [V]) = λ max k E [X k ] .
Then, for all h ≥ 0,

P (λ max (V) ≥ (1 + h)µ max ) ≤ d • e h (1 + h) 1+h µmax L ≤ d • exp - µ max L ((h + 1) log(h + 1) -h) •
In our case, we have X {e,j} = 1 q 2 1 w 0,e,j q e q T e q e q T e 1 q 2 α 2 p 2 h,e q e q T e q e q T e 1 q 2 α 2 p 2 h,e q e q T e 2 I α 2 q 2 I, where the first inequality follows from the definition of w 0,e,j in Eq. 14, the second from the PSD ordering, and the third from the definition of q e q T e . Therefore, we can use L α 2 /q 2 for the purpose of Prop. 5. We need now to compute E [X k ], that we can use in turn to compute µ max . We begin by computing the expected value of 1/w 0,e,j . Let us denote the c.d.f. of 1/w 2 0,e,j as

F 1/w 2 0,e,j (a) = P 1 w 2 0,e,j ≤ a = P 1 w 0,e,j ≤ √ a =      0 for a < 1 1 -1 √ a for 1 ≤ a < α 2 /p 2 h,e 1
for α 2 /p 2 h,e ≤ a .

Since P 1/w 2 0,e,j ≥ 0 = 1, we have that (1 -0) da +

E 1 w 0,e,j = ∞ a=0 1 -F 1/w0,e,j (
α 2 /p 2 h,e a=1 1 -1 - 1 √ a da + ∞ a=α 2 /p 2 h,e (1 -1) da = 1 a=0 da + α 2 /p 2 h,e a=1 1 √ a da = 1 + [2 √ a] α 2 /p 2 h,e 1 = 2α/p h,e -1.
Therefore,

µ max = λ max (E [V]) = λ max {e,j} E X {e,j} = λ max   1 q 2 q j=1 m e=1 E 1 w 2 0,e,j
q e q T e q e q T e   = λ max 1 q m e=1 2α p h,e -1 p h,e q e q T e ≤ λ max 2α q m e=1 q e q T e = 2α q λ max (P) ≤ 2α q L.

Therefore, selecting h = 2, σ 2 = 6α/q and applying Prop. 5 we have

P W h ≥ σ 2 ≤ P   λ max   1 q 2 q j=1 m e=1 1 w 2 0,e,j q e q T e q e q T e   ≥ (1 + 2) 2α q   ≤ m • exp - 2α q q 2 α 2 (3 log(3) -2) ≤ n • exp - 2q α • B.4.

Space complexity bound

Denote with A the event A = ∀h ∈ {1, . . . , h} : P h -P h 2 ≤ ε , and again m = |G {h, } |. To begin with, we can show that under event A, the approximate γ-effective resistances (i.e. RLS) are accurate (Lem. 2, [START_REF] Calandriello | Distributed sequential sampling for kernel matrix approximation[END_REF].

Let H = H {h, } be a (ε, 2γ)-sparsifier of G = G {h, } . Using Prop. 3 it is straightforward to see that under A we have

p h+1,e ≤ r h+1,e (γ) = (1 -ε)b T e (L H + (1 + ε)γI) -1 b e ≤ (1 -ε)b T e ((1 -ε)L G -2εγI + (1 + ε)γI) -1 b e = (1 -ε) (1 -ε) b T e (L G + γI) -1 b e = r h+1,e = p h+1,e .
Letting q = |H {h, } | = m e=1 q h,e = q j=1 m e=1 z h,e,j be the random number of edges in H {h, } , we reformulate

P |H {h, } | ≥ 3qd {h, } eff (γ) ∩ ∀h ∈ {1, . . . , h} : P h -P h 2 ≤ ε ≤ ε = P |H {h, } | ≥ 3qd {h, } eff (γ) ∩ A = P   q j=1 m e=1 z h,e,j ≥ 3qd {h, } eff (γ) ∩ A   = P   q j=1 m e=1 z h,e,j ≥ 3qd {h, } eff (γ) A   P (A).
While we do know that the z h,e,j are Bernoulli random variables (since they are either 0 or 1), it is not easy to compute the success probability of each z h,e,j , and in addition there could be dependencies between z h,e,j and z h,e ,j . Similarly to Lem. 2, we are going to find a stochastic variable to dominate z h,e,j . Denoting with u s,e,j ∼ U(0, 1) a uniform random variable, we will define w s,e,j as for any e and j such that {s, 1, 1} ≤ {s, e , j } < {s, e, j}. Note that w s,e,j , unlike z s,e,j , does not have a recursive definition, and its only dependence on any other variable comes from p s-1,e . First, we peel off the last step

w s,
P   q j=1 m e=1 z h,e,j ≥ g A   = E Ft-1|A   P   q j=1 m e=1 I u h,e,j ≤ p h,e p t-1,e z h-1,e,j ≥ g F t-1 ∩ A     ≤ E Ft-1|A   P   q j=1 m e=1 I u h,e,j ≤ p h,e p t-1,e z h-1,e,j ≥ g F t-1 ∩ A     = P   q j=1 m e=1 w h,e,j z h-1,e,j ≥ g A   ,
where we used the fact that conditioned on A, H {h, } is an (ε, γ)-sparsifier of G {h, } , which guarantees that p h,e ≤ p h,e . Plugging this in the previous bound,

P   q j=1 m e=1 z h,e,j ≥ g A   P (A) ≤ P   q j=1 m e=1 w h,e,j z h-1,e,j ≥ g ∩ A   ≤ P   q j=1 m e=1 w h,e,j z h-1,e,j ≥ g   .
We now proceed by peeling off layers from the end of the chain one by one. We show how to move from an iteration s ≤ h to s -1,

C. Further details on experiments

Software and hardware All our code is implemented in Julia and runs on a cluster of 4 machines with 128 GB of RAM and a 10-core Xeon E5-2630. The distributed computation is achieved using a suboptimal but simple producer-consumer queue.

Linear solver We use approximate Gaussian elimination scheme by Kyng & Sachdeva (2016) from Laplacians.jl9 package.

Details on DiSRe. To compute all (ε, γ) sparsifiers, we use DiSRe after splitting the input graph into 8 sub-graphs,10 resulting in 3 rounds of resparsifications using 4 machines. For each resparsification, we compute the effective resistance in parallel on each machine using 10 processes.11

Additional results. For completeness, in the following with provide the tables containing all the combinations of hyperparameters (ε, γ, q, k, σ, l) used in our experiments. As the supplementary files, we also provide two OpenDocument spreadsheets containing these results, one for the Laplacian smoothing experiment (extra_results_smoothing.ods) and one for the SSL experiment (extra_results_ssl.ods).

Figure 1 .

 1 Figure 1. Merge tree for Alg. 1.

Figure 3 .

 3 Figure 3. The dependence graph of the variables. Red variables are random. Black variables are deterministically computed using their input (a function of their input), with bold lines indicating the deterministic (functional) relation. Blue variables are constants.A grey filling indicates that a random variable is observed or a function of observed variables.

  Applying this repeatedly from s = h to s = 2, we have P e,j z h-1,e,j ≥ g e,j z 0,e,j ≥ g the w 1,e,j are independent Bernoulli random variables, and we can bound their sum with a Hoeffding-like bound using Markov inequality, q(e θ -1)-θg) ,where we use the fact that 1 + x ≤ e x , w 1,e,j ∼ B(p h,e ) and by Def. 2, The choice of θ minimizing the previous expression is obtained asd dθ e qd {h, } eff (γ)(e θ -1)-θg = e qd {h, } eff (γ)(e θ -1)-θg qd {h, } eff (γ)e θg = 0,and thus θ = log(g/(qd {h, } eff (γ))). Plugging this in the previous bound, (e θ -1)θg) = exp gqd we conclude our proof.

  Algorithm 1 The DiSRe algorithm.

	Input: G Output: H G 1: Partition G into k sub-graphs: 2:

Table 1 .

 1 Results for the SSL and the smoothing problems. uate if they can accelerate computation while preserving accuracy. We run DiSRe with different values of γ depending on the setting. For empirically strong heuristics, we attempted to uniformly subsample the edges, but at the sparsity level achieved by the other methods, the uniformly sampled sparsifier is disconnected and highly inaccurate.

	Alg. |E| (x10 EXACT Parameters 98.5	0.312 ± 0.022	0.286 ± 0.010	0.067 ± 0.0004	0.756 ± 0.006
	kN	k = 60	15.7	0.329 ± 0.0143	0.311 ± 0.027	0.172 ± 0.0004	0.822 ± 0.002
	kN	k = 90	21.2	0.334 ± 0.024	0.311 ± 0.024	0.125 ± 0.0002	0.811 ± 0.003
	DiSRe	γ = 0, q = 100	15	0.314 ± 0.0165	0.296 ± 0.015	0.068 ± 0.0003	0.758 ±0.005
	DiSRe	γ = 0, q = 150	22.8	0.314 ± 0.0158	0.310 ± 0.024	0.068 ± 0.0004	0.756 ± 0.005
	DiSRe	γ = 10 3 , q = 100	7.3	-	-	0.072 ± 0.0003	0.789 ± 0.005
	DiSRe	γ = 10 2 , q = 100	11.8	-	-	0.068 ± 0.0002	0.772 ± 0.004
	DiSRe	γ = 10, q = 100	14.4	-	-	0.068 ± 0.0004	0.760 ± 0.004

6

) Err. SSL ( = 346) Err. SSL ( = 672) Err. D( f )(σ = 10 -3 ) Err. D( f ) (σ = 10 -2 ) Instead, we compare to the state-of-the-art k-neighbors (kN) heuristic by

[START_REF] Sadhanala | Graph sparsification approaches for laplacian smoothing[END_REF]

, which is just as fast as uniform sampling and more accurate in practice.

  1, 1/p h,e ≤ 1/w h,e,j . Therefore, Now focus on an arbitrary intermediate step 1 ≤ k ≤ h, where we fix F k-1 . Since u k,e,j and w k,e,j are independent given Figure 4. C.d.f. of max z k-1,e,j /p t,e,j ; z k,e,j /p k,e,j and z k-1,e,j /w k-1,e,j conditioned on F {k-1} . For conciseness, we omit the e, j indices.

			z t,e,j =1
	1-p t,e / pt,e 1		P max 1 pt,e ; z t +1,e,j pt+1,e ≤a | F t,m,N
				P 1 w t,e,j ≤a | F t,m,N
	1 1/ pt,e 0			α 2 /p t,e
	P max s=0,...,h	z s,e,j p s,e,j	≤ a ≥ P max	max s=0,...,h-1	z s,e,j p s,e,j	;	z h,e,j w h,e,j	≤ a .

  e,j |F {s,e ,j } w s,e,j |F s-2 I u s,e,j ≤

	p h,e p s-1,e	∼ B	p h,e p s-1,e

  Improved Large-Scale Graph Learning through Ridge Spectral Sparsification Left: DiSRe, q = 100, γ = 100, Middle: DiSRe, q = 100, γ = 10, Right: kN, k = 60

	σ	λ	mean		std	σ	λ		mean		std	σ	λ		mean	std
	0.001 0.001	0.1905		0.0004	0.001 0.001		0.1907		0.0007	0.001 0.001		0.3024	0.0010
	0.001	0.01	0.0689		0.0003	0.001	0.01		0.0681		0.0004	0.001	0.01		0.1724	0.0004
	0.001	0.1	0.2814		0.0010	0.001	0.1		0.2841		0.0005	0.001	0.1		0.2838	0.0008
	0.001	1	0.7925		0.0010	0.001	1		0.7983		0.0002	0.001	1		0.7856	0.0006
	0.001	10	0.9716		0.0002	0.001	10		0.9727		0.0001	0.001	10		0.9706	0.0001
	0.01 0.001	19.0622		0.0340	0.01 0.001		19.0317		0.0489	0.01 0.001		30.1645	0.0703
	0.01	0.01	5.3400		0.0122	0.01	0.01		5.3175		0.0253	0.01	0.01		15.7825	0.0399
	0.01	0.1	0.7727		0.0047	0.01	0.1		0.7606		0.0046	0.01	0.1		1.9238	0.0107
	0.01	1	0.8086		0.0029	0.01	1		0.8127		0.0027	0.01	1		0.8227	0.0027
	0.01	10	0.9719		0.0004	0.01	10		0.9730		0.0006	0.01	10		0.9712	0.0005
	0.1 0.001	1906.6166		6.2333	0.1 0.001		1905.0821		4.0734	0.1 0.001	3022.1232	7.4520
	0.1	0.01	533.4745		2.2259	0.1	0.01		531.3563		1.5535	0.1	0.01	1575.6447	3.3189
	0.1	0.1	49.6360		0.3576	0.1	0.1		48.2788		0.3639	0.1	0.1		166.1367	0.6022
	0.1	1	2.3184		0.0440	0.1	1		2.1985		0.0424	0.1	1		4.5508	0.0500
	0.1	10	1.0092		0.0084	0.1	10		1.0191		0.0240	0.1	10		1.0276	0.0111
	1 0.001 190451.4911 600.4823	1 0.001 190205.0504 370.3744	1 0.001 302265.6761 375.8875
	1	0.01	53479.7617 251.9642	1	0.01	53118.5249 176.1894	1	0.01 157503.6242 369.7657
	1	0.1	4961.0662	44.4829	1	0.1		4809.2633	34.7541	1	0.1	16585.2588	43.9628
	1	1	152.0926		4.5554	1	1		140.9587		3.7546	1	1		380.2291	5.5260
	1	10	4.4540		2.1627	1	10		4.8436		1.5676	1	10		6.8515	1.1969
			LAPSMO, σ	λ		mean		std	σ	λ		mean		std
			0.001 0.001		0.2823	0.0007	0.001 0.001		0.1896		0.0002
			0.001	0.01		0.1259	0.0003	0.001	0.01		0.0676		0.0005
			0.001	0.1		0.2834	0.0008	0.001	0.1		0.2856		0.0002
			0.001	1		0.7929	0.0005	0.001	1		0.8000		0.0002
			0.001	10		0.9719	0.0001	0.001	10		0.9730		0.0000
			0.01 0.001		28.2517	0.0560	0.01 0.001		18.9408		0.0554
			0.01	0.01		11.0908	0.0309	0.01	0.01		5.2789		0.0147
			0.01	0.1		1.1154	0.0087	0.01	0.1		0.7566		0.0067
			0.01	1		0.8119	0.0029	0.01	1		0.8127		0.0028
			0.01	10		0.9722	0.0006	0.01	10		0.9737		0.0003
			0.1 0.001	2823.5932	6.3830	0.1 0.001	1894.0375		5.7404
			0.1	0.01	1108.9582	5.3674	0.1	0.01		527.5392		1.0052
			0.1	0.1		84.0189	0.4347	0.1	0.1		47.7877		0.4702
			0.1	1		2.8117	0.0394	0.1	1		2.1503		0.0282
	σ	λ	mean 0.1 1 0.001 282388.0211 754.7190 std σ λ 10 1.0124 0.0178	mean 0.1 1 0.001 189515.7627 293.8406 std σ λ 10 1.0083 0.0184	mean	std
	0.001 0.001	0.1908 1	0.0006 0.01 110729.9273 250.8183 0.001 0.001	0.1902 1	0.01	0.0004 52639.4113 346.7699 0.001 0.001	0.1910	0.0007
	0.001	0.01	0.0681 1	0.1	0.0003 8373.5778 0.001	0.01 55.8408	0.0681 1	0.1	0.0004 4724.3469 0.001	0.01 44.4066	0.0723	0.0004
	0.001	0.1	0.2845 1	1	0.0006	0.001 200.9318	0.1 3.7598	0.2846 1	1	0.0005	0.001 138.3458	0.1	2.9040	0.2622	0.0020
	0.001	1	0.7991 1	10	0.0002	0.001 4.8285	1 1.0723	0.7993 1	10	0.0003	0.001 4.1240	1	0.9370	0.7198	0.0034
	0.001	10	0.9728		0.0001	0.001	10		0.9729		0.0001	0.001	10		0.9096	0.0047
	0.01 0.001	19.0389		0.0464	0.01 0.001		18.9854		0.0510	0.01 0.001		19.0615	0.0551
	0.01	0.01	5.3282		0.0273	0.01	0.01		5.2975		0.0240	0.01	0.01		5.5206	0.0143
	0.01	0.1	0.7587		0.0059	0.01	0.1		0.7565		0.0053	0.01	0.1		0.9020	0.0088
	0.01	1	0.8129		0.0025	0.01	1		0.8124		0.0014	0.01	1		0.7890	0.0050
	0.01	10	0.9731		0.0003	0.01	10		0.9731		0.0004	0.01	10		0.9523	0.0056
	0.1 0.001	1904.4359		5.9497	0.1 0.001		1899.5487		8.1249	0.1 0.001	1906.4671	2.9946
	0.1	0.01	530.8727		2.7780	0.1	0.01		528.6970		1.6710	0.1	0.01		550.5067	1.7484
	0.1	0.1	48.0577		0.4268	0.1	0.1		48.0165		0.3385	0.1	0.1		64.5503	0.6398
	0.1	1	2.1887		0.0547	0.1	1		2.1773		0.0563	0.1	1		7.4914	0.3399
	0.1	10	1.0044		0.0136	0.1	10		1.0041		0.0086	0.1	10		5.3122	0.3776
	1 0.001 190532.1581 437.1137	1 0.001 189582.9265 413.1919	1 0.001 190769.1154 348.4353
	1	0.01	53197.4003 187.2996	1	0.01	52905.9901 169.3004	1	0.01	54985.9953 241.7914
	1	0.1	4770.8654	39.7099	1	0.1		4775.0883	18.4023	1	0.1	6484.4296	73.4393
	1	1	140.4328		3.7919	1	1		137.7268		2.7805	1	1		688.6995	39.1172
	1	10	5.2586		2.0544	1	10		4.7664		1.9795	1	10		416.7218	40.3694

LAPSMO, Left: DiSRe, q = 100, Middle: DiSRe, q = 150, Right: DiSRe, q = 100, γ = 1000 LAPSMO, Left: kN, k = 90, Right: EXACT

Many algorithms reduce the complexity of graph learning at construction time but they cannot be applied to natural graphs (e.g., social graphs) and therefore we do not review them.

Whenever no regularization is required in the learning task (i.e., HFS, SC), we set γ = 0 and consider "standard" ε-sparsifiers.

Whenever the original graph contains m ≤ O(deff(γ))edges, there is no need to run DiSRe as the (ε, γ)-sparsifiers would not reduce the size of the graph.

Whenever this is not possible (i.e., m/k is too large to be

This implies that there are enough machines so that the leaves in the merge tree already have relatively sparse sub-graphs.

This implies that there are enough machines so that the leaves in the merge tree already have relatively sparse sub-graphs.

Notice that unlike zs,e,j, ws,e,j is no longer Fs-measurable but it is F s -measurable, where

http://github.com/danspielman/Laplacians.jl

To guarantee that each of the sub-graphs is defined on the same set of nodes, we pre-construct a spanning tree of G and include it in each of the sub-graphs. Note that the analysis holds even if each sub-graph is disconnected from the others. We choose this approach to avoid complicating the code with additional searches of connected components in the sparsifiers.

Laplacian.jl is strictly single-threaded. Therefore, we use multiple processes on a single machine. Faster runtime and lower memory usage could be achieved by sharing the memory with threads.

Acknowledgements

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RE-NATER, and several universities as well as other organizations (see https://www.grid5000.fr). The research presented was also supported by European CHIST-ERA project DELTA, French Ministry of Higher Education and Research, Nord-Pas-de-Calais Regional Council, Inria and Otto-von-Guericke-Universität Magdeburg associated-team north-european project Allocate, and French National Research Agency projects ExTra-Learn (n.ANR-14-CE24-0010-01) and BoB (n.ANR-16-CE23-0003). I. Koutis is supported by NSF CAREER award CCF-1814563.