Daniele Calandriello

Alessandro Lazaric

Michal Valko

Distributed Adaptive Sampling for Kernel Matrix Approximation

 show that sampling points with replacement according to their ridge leverage scores (RLS) generates small dictionaries of relevant points with strong spectral approximation guarantees for K n . The drawback of RLS-based methods is that computing exact RLS requires constructing and storing the whole kernel matrix. In this paper, we introduce SQUEAK, a new algorithm for kernel approximation based on RLS sampling that sequentially processes the dataset, storing a dictionary which creates accurate kernel matrix approximations with a number of points that only depends on the effective dimension d eff (γ) of the dataset. Moreover since all the RLS estimations are efficiently performed using only the small dictionary, SQUEAK never constructs the whole matrix K n , runs in linear time O(nd eff (γ) 3 ) w.r.t. n, and requires only a single pass over the dataset. We also propose a parallel and distributed version of SQUEAK achieving similar accuracy in as little as O(log(n)d eff (γ) 3 ) time.

Introduction

One of the major limits of kernel ridge regression (KRR), kernel PCA [START_REF] Schölkopf | Kernel principal component analysis[END_REF], and other kernel methods is that for n samples storing and manipulating the kernel matrix K n requires O(n 2 ) space, which becomes rapidly infeasible for even a relatively small n. For larger sizes (or streams) we cannot even afford to store or process the data on as single machine. Many solutions focus on how to scale kernel methods by reducing its space (and time) complexity without compromising the prediction accuracy. A popular approach is to construct low-rank approximations of the kernel matrix by randomly selecting a subset (dictionary) of m columns from K n , thus reducing the space complexity to O(nm). These methods, often referred to as Nyström approximations, mostly differ in the distribution used to sample the columns of K n and the construction of low-rank approximations. Both of these choices significantly affect the accuracy of the resulting approximation [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF]. Bach [START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF] showed that uniform sampling preserves the prediction accuracy of KRR (up to ε) only when the number of columns m is proportional to the maximum degree of freedom of the kernel matrix. This may require sampling O(n) columns in datasets with high coherence [START_REF] Gittens | Revisiting the Nyström method for improved large-scale machine learning[END_REF], i.e., a kernel matrix with weakly correlated columns. On the other hand, Alaoui and Mahoney [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] showed that sampling columns according to their ridge leverage scores (RLS) (i.e., a measure of the influence of a point on the regression) produces an accurate Nyström approximation with only a number of columns m proportional to the average degrees of freedom of the matrix, called effective dimension. Unfortunately, the complexity of computing RLS requires storing the whole kernel matrix, thus making this approach infeasible. However, Alaoui and Mahoney [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] proposed a fast method to compute a constant-factor approximation of the RLS and showed that accuracy and space complexity are close to the case of sampling with exact RLS at the cost of an extra dependency on the inverse of the minimal eigenvalue of the kernel matrix. Unfortunately, the minimal eigenvalue can be arbitrarily small in many problems. Calandriello et al. [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF] addressed this issue by processing the dataset incrementally and updating estimates of the ridge leverage scores, effective dimension, and Nyström approximations on-thefly. Although the space complexity of the resulting algorithm (INK-Estimate) does not depend on the minimal eigenvalue anymore, it introduces a dependency on the largest eigenvalue of K n , which in the worst case can be as big as n, thus losing the ad-vantage of the method. In this paper we introduce an algorithm for SeQUEntial Approximation of Kernel matrices (SQUEAK), a new algorithm that builds on INK-Estimate, but uses unnormalized RLS. This improvement, together with a new analysis, opens the way to major improvements over current leverage sampling methods (see Sect. 6 for a comparison with existing methods) closely matching the dictionary size achieved by exact RLS sampling. First, unlike INK-Estimate, SQUEAK is simpler, does not need to compute an estimate of the effective dimension for normalization, and exploits a simpler, more accurate RLS estimator. This new estimator only requires access to the points stored in the dictionary. Since the size of the dictionary is much smaller than the n, SQUEAK needs to actually observe only a fraction of the kernel matrix K n , resulting in a runtime linear in n. Second, since our dictionary updates require only access to local data, our algorithm allows for distributed processing where machines operating on different dictionaries do not need to communicate with each other. In particular, intermediate dictionaries can be extracted in parallel from small portions of the dataset and they can be later merged in a hierarchical way. Third, the sequential nature of SQUEAK requires a more sophisticated analysis that take into consideration the complex interactions and dependencies between successive resampling steps. The analysis of SQUEAK builds on a new martingale argument that could be of independent interest for similar online resampling schemes. Moreover, our SQUEAK can naturally incorporate new data without the need of recomputing the whole resparsification from scratch and therefore it can be applied in streaming settings. We note there exist other ways to avoid the intricate dependencies with simpler analysis, for example by resampling [START_REF] Musco | Provably useful kernel matrix approximation in linear time[END_REF], but with negative algorithmic side effects: these methods need to pass through the dataset multiple times. SQUEAK passes through the dataset only once 1 and is therefore the first provably accurate kernel approximation algorithm that can handle both streaming and distributed settings.

Background

In this section, we introduce the notation and basics of kernel approximation used through the paper.

Notation. We use curly capital letters A for collections. We use upper-case bold letters A for matrices and operators, lower-case bold letters a for vectors, and lower-case letters a for scalars, with the exception of f, g, and h which denote functions. We denote by [A] ij and [a] i , the (i, j) element of a matrix and ith element of a vector respectively. We denote by I n ∈ R n×n , the identity matrix of dimension n and by Diag(a) ∈ R n×n the diagonal matrix with the vector a ∈ R n on the diagonal. We use e n,i ∈ R n to denote the indicator vector for element i of dimension n. When the dimension of I and e i is clear from the context, we omit the n. We use A B to indicate that A -B is a Positive Semi-Definite (PSD) matrix or operator. Finally, the set of integers between 1 and n is denoted as [n] := {1, . . . , n}, and between i and j as [i : j] := {i, . . . , j}.

Kernel. We consider a positive definite kernel function K : X × X → R and we denote with H its induced Reproducing Kernel Hilbert Space (RKHS), and with ϕ : X → H its corresponding feature map. Using ϕ, and without loss of generality, for the rest of the paper we will replace H with a high dimensional space R D where D is large and potentially infinite. With this notation, the kernel evaluated between to points can be expressed as K(x, x ) = K(x, •), K(x , •) H = ϕ(x), ϕ(x ) H = ϕ(x) T ϕ(x ). Given a dataset of points D = {x t } n t=1 , we define the (empirical) kernel matrix K t ∈ R t×t as the application of the kernel function on all pairs of input values (i.e., [K t ] ij = k i,j = K(x i , x j ) for any i, j ∈ [t]), with k t,i = K t e t,i as its i-th column. We also define the feature vectors φ i = ϕ(x i ) ∈ R D and after introducing the matrix Φ t = [φ 1 , φ 2 , . . . , φ t ] ∈ R D×t we can rewrite the kernel matrix as K t = Φ T t Φ t . Kernel approximation by column sampling. One of the most popular strategies to have low space complexity approximations of the kernel K t is to randomly select a subset of its columns (possibly reweighted) and use them to perform the specific kernel task at hand (e.g., kernel regression). More precisely, we define a column dictionary as a collection

I t = {(i, w i )} t i=1
, where the first term denotes the index of the column and w i its weight, which is set to zero for all columns that are not retained. For the theoretical analysis, we conveniently keep the dimension of any dictionary I t to t, while in practice, we only store the non-zero elements. In particular, we denote by |I t | be the size of the dictionary corresponding to the elements with non-zero weights w i . Associated with a column dictionary, there is a selection matrix S t = Diag( √ w 1 . . . √ w t ) ∈ R t×t such that for any matrix A t ∈ R t×t , A t S t returns a t × t matrix where the columns selected by I t are properly reweighted and all other columns are set to 0. Despite the wide range of kernel applications, it is possible to show that in most of them, the quality of a dictionary can be mea-sured in terms of how well it approximates the projection associated to the kernel. In kernel regression, for instance, we use K t to construct the projection (hat) matrix that projects the observed labels y t to y t . In particular, let P t = K t K + t be the projection matrix (where K + t indicates the pseudoinverse), then y t = P t y t . If K t is full-rank, then P t = I t is the identity matrix and, we can reconstruct any target vector y t exactly. On the other hand, the only sampling scheme which guarantees to properly approximate a full rank P t requires all columns to be represented in I t . In fact, all columns have the same "importance" and no low-space approximation is possible. Nonetheless, kernel matrices are often either rank deficient or have extremely small eigenvalues (exponentially decaying spectrum), as a direct (and desired) consequence of embedding low dimensional points x i into a high dimensional RKHS. In this case, after soft thresholding the smaller eigenvalues to a given value γ, K t can be effectively approximated using a small subset of columns. This is equivalent to approximating the γridge projection matrix

P t def =(K t + γI) -1/2 K t (K t + γI) -1/2 .
We say that a column dictionary is accurate if the following condition is satisfied. Definition 1. A dictionary I t = {(i, w i )} t i=1 and its associated selection matrix S t ∈ R t×t are ε-accurate w.r.t. a kernel matrix

K t = K 1/2 t K 1/2 t if 2 P t -P t ≤ ε, (1)
where for a given γ > 0, the approximated projection matrix is defined as

P t def =(K t + γI t ) -1 2 K 1/2 t S t S T t K 1/2 t (K t + γI t ) -1 2 .
Notice that this definition of accuracy is purely theoretical, since P t is never computed. Nonetheless, as illustrated in Sect. 5, ε-accurate dictionaries can be used to construct suitable kernel approximation in a wide range of problems.

Ridge leverage scores sampling. Alaoui and Mahoney [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] showed that an ε-accurate dictionary can be obtained by sampling columns proportionally to their γ-ridge leverage scores (RLS) defined as follows. Definition 2. Given a kernel matrix

K t ∈ R t×t , the γ-ridge leverage score (RLS) of column i ∈ [t] is τ t,i = e T t,i K t (K t + γI t ) -1 e t,i , (2)
Furthermore, the effective dimension d eff (γ) t of the kernel matrix K t is defined as

d eff (γ) t = t i=1 τ t,i (γ) = Tr K t (K t + γI t ) -1 . ( 3
)
2 the matrix norm we use is the operator (induced) norm

I = I t-1 ∪ {(t, p t-1,t = 1, q t-1,t = q)} Expand 5:
I t = Dict-Update(I) using Eq. 4

if q t-1,i = 0 then 4:
Compute τ t,i using I

5:

Set p t,i = min{ τ t,i , p t-1,i }

6:

Set q t,i ∼ B( p t,i / p t-1,i , q t-1,i ) p t,i = p t-1,i and q t,i = q t-1,i

9:

end if 10: end for The RLS can be interpreted and derived in many ways, and they are well studied [START_REF] Michael | Input sparsity time low-rank approximation via ridge leverage score sampling[END_REF][START_REF] Michael | Online row sampling[END_REF][START_REF] David | Sketching as a tool for numerical linear algebra[END_REF] in the linear setting (e.g. φ t = x t ). Patel et al. [START_REF] Patel | oA-SIS: Adaptive column sampling for kernel matrix approximation[END_REF] used them as a measure of incoherence to select important points, but their deterministic algorithm provides guarantees only when K t is exactly low-rank. Here we notice that τ t,i = e T t,i K t (K t + γI t ) -1 e t,i = e T t,i P t e t,i , which means that they correspond to the diagonal elements of the P t itself. Intuitively, this correspond to selecting each column i with probability p t,i = τ t,i will capture the most important columns to define P t , thus minimizing the approximation error P t -P t . More formally, Alaoui and Mahoney [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] state the following.

Proposition 1. Let ε ∈ [0, 1] and I n be the dictionary built with m columns randomly selected proportionally to RLSs {τ n,i } with weight

w i = 1/(mτ n,i ). If m = O( 1 ε 2 d eff (γ) n log( n δ )
), then w.p. at least 1δ, the corresponding dictionary is ε-accurate.

Unfortunately, computing exact RLS requires storing K n and this is seldom possible in practice. In the next section, we introduce SQUEAK, an RLS-based incremental algorithm able to preserve the same accuracy of Prop. 1 without requiring to know the RLS in advance. We prove that it generates a dictionary only a constant factor larger than exact RLS sampling.

In the previous section, we showed that sampling proportionally to the RLS {τ t,i } leads to a dictionary such that P t -P t ≤ ε. Furthermore, since the RLS correspond to the diagonal entries of P t , an accurate approximation P t may be used in turn to compute accurate estimates of τ t,i . The SQUEAK algorithm (Alg. 1) builds on this intuition to sequentially process the kernel matrix K n so that exact RLS computed on a small matrix (K t with t n) are used to create an ε-accurate dictionary, which is then used to estimate the RLS for bigger kernels, which are in turn used to update the dictionary and so on. While SQUEAK shares a similar structure with INK-Estimate [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF], the sampling probabilities are computed from different estimates of the RLS τ t,i and no renormalization by an estimate of d eff (γ) t is needed. Before giving the details of the algorithm, we redefine a dictionary as a collection I = {(i, p i , q i )} i , where i is the index of the point x i stored in the dictionary, p i tracks the probability used to sample it, and q i is the number of copies (multiplicity) of i. The weights are then computed as w i = q i /(q p i ), where q is an algorithmic parameter discussed later. We use p i to stress the fact that these probabilities will be computed as approximations of the actual probabilities that should be used to sample each point, i.e., their RLS τ i .

SQUEAK receives as input a dataset D = {x t } n t=1 and processes it sequentially. Starting with an empty dictionary I 0 , at each time step t, SQUEAK receives a new point x t . Adding a new point x t to the kernel matrix can either decrease the importance of points observed before (i.e., if they are correlated with the new point) or leave it unchanged (i.e., if their corresponding kernel columns are orthogonal) and thus for any i ≤ t, the RLS evolves as follows.

Lemma 1. For any kernel matrix K t-1 at time t -1 and its extension K t at time t, we have that the RLS are monotonically decreasing and the effective dimension is monotonically increasing,

1 τ t-1,i + 1 τ t-1,i ≤ τ t,i ≤ τ t-1,i , d eff (γ) t ≥ d eff (γ) t-1 .
The previous lemma also shows that the RLS cannot decrease too quickly and since τ t-1,i ≤ 1, they can at most halve when τ t-1,i = 1. After receiving the new point x t , we need to update our dictionary I t-1 to reflect the changes of the τ t,i . We proceed in two phases. During the Expand phase, we directly add the new element x t to I t-1 and obtain a temporary dictionary I, where the new element t is added with a sampling probability p t-1,t = 1 and a number of copies q t-1,t = q, i.e., I = I t-1 ∪ {(t, p t-1,t = 1, q t-1,t = q)}. This increases our memory usage, forcing us to update the dictionary using Dict-Update, in order to decrease its size. Given as input I, we use the following estimator to compute the approximate RLS τ t,i ,

τ t,i = (1 -ε)φ T i (Φ t SS T Φ T t + γI) -1 φ i = 1-ε γ (k i,i -k T t,i S(S T K t S + γI t ) -1 S T k t,i ), ( 4
)
where ε is the accuracy parameter, γ is the regularization and S is the selection matrix associated to I. This estimator follows naturally from a reformulation of the RLS. In particular, if we consider φ i , the RKHS representation of x i , the RLS τ t,i can be formulated as τ t,i = φ T i (Φ t I t Φ T t + γI) -1 φ i , where we see that the importance of point x i is quantified by how orthogonal (in the RKHS) it is w.r.t. the other points. Because we do not have access to all the columns (SS

T = I t ),
similarly to what [START_REF] Michael | Online row sampling[END_REF] did for the special case φ i = x i , we choose to use τ t,i ≈ φ T (Φ t SS T Φ T t + γI) -1 φ i , and then we use the kernel trick to derive a form that we can actually compute, resulting in Eq. 4. The approximate RLSs are then used to define the new sampling probabilities as p t,i = min{ τ t,i , p t-1,i }. For each element in I, the Shrink step draws a sample from the binomial B( p t,i / p t-1,i , q t-1,i ), where the minimum taken in the definition of p t,i ensures that the binomial probability is well defined (i.e., p t,i ≤ p t-1,i ). This resampling step basically tracks the changes in the RLS and constructs a new dictionary I t , which is as if it was created from scratch using all the RLS up to time t (with high probability). We see that the new element x t is only added to the dictionary with a large number of copies (from 0 to q) if its estimated relevance p t,t is high, and that over time elements originally in I t-1 are stochastically reduced to reflect the reductions of the RLSs. The lower p t,i w.r.t. p t-1,i , the lower the number of copies q t,i w.r.t. q t-1,i . If the probability p t,i continues to decrease over time, then q t,i may become zero, and the column i is completely dropped from the dictionary (by setting its weight to zero). The approximate RLSs enjoy the following guarantees.

Lemma 2. Given an ε-approximate dictionary I t-1 of matrix K t-1 , construct I by adding element (t, 1, q) to it, and compute the selection matrix S. Then for all i in I such that q t-1,i = 0, the estimator in Eq. 4 is α-accurate, i.e., it satisfies τ t,i /α ≤ τ t,i ≤ τ t,i , with α = (1 + ε)/(1ε). Moreover, given RLS τ t-1,i and τ t,i , and two α-accurate RLSs, τ t-1,i and τ t,i , the quantity min { τ t,i , τ t-1,i } is also an α-accurate RLS.

This result is based on the property that whenever I t-1 is ε-accurate for K t-1 , the projection matrix P t can be approximated by P t constructed using the temporary dictionary I and thus, the RLSs can be accurately estimated and used to update I t-1 and obtain a new ε-accurate dictionary for K t . Since τ t,i is used to sample the new dictionary I t , we need each point to be sampled almost as frequently as with the true RLS τ t,i , which is guaranteed by the lower bound of Lem. 2. Since RLSs are always smaller or equal than 1, this could be trivially achieved by setting τ t,i to 1. Nonetheless, this would keep all columns in the dictionary. Consequently, we need to force the RLS estimate to decrease as much as possible, so that low probabilities allow reducing the space as much as possible. This is obtained by the upper bound in Lem. 2, which guarantees that the estimated RLS are always smaller than the exact RLS. As a result, Shrink sequentially preserves the overall accuracy of the dictionary and at the same time keeps its size as small as possible, as shown in the following theorem.

Theorem 1. Let ε > 0 be the accuracy parameter, γ > 1 the regularization, and 0 < δ < 1 the probability of failure. Given an arbitrary dataset D in input together with parameters ε, γ, and δ, we run SQUEAK with

q = 39α log (2n/δ) ε 2 , where α = (1 + ε)/(1 -ε).
|I t | ≤ 3qd eff (γ) n .
As a consequence, on a successful run the overall complexity of SQUEAK is bounded as space complexity = max t=1,...,n

|I t | 2 ≤ (3qd eff (γ) n ) 2 , time complexity = O nd eff (γ) 3 n q 3 .
We show later that Thm. 1 is special case of Thm. 2 and give a sketch of the proof with the statement of Thm. 2. We postpone the discussion about this result and the comparison with previous results to Sect. 6 and focus now on the space and time complexity. Note that while the dictionaries I t always contain t elements for notational convenience, Shrink actually never updates the probabilities of the elements with q t-1,i = 0. This feature is particularly important, since at any step t, it only requires to compute approximate RLSs for the elements which are actually included in I t-1 and the new point x t (i.e., the elements in I) and thus it does not require recomputing the RLSs of points x s (s < t) that have been dropped before! This is why SQUEAK computes an ε-accurate dictionary with a single pass over the dataset. Furthermore, the estimator in Eq. 4 does not require computing the whole kernel column k t,i of dimension t. In fact, the components of k t,i , corresponding to points which are no longer in I, are directly set to zero when computing Algorithm 2 The distributed SQUEAK algorithm Input: Dataset D, parameters γ, ε, δ Output:

I D 1: Partition D into disjoint sub-datasets D i 2: Initialize I Di = {(j, p 0,i = 1, q 0,i = q) : j ∈ D i } 3: Build set S 1 = {I Di } k i=1 4: for h = 1, . . . , k -1 do 5: if |S h | > 1 then Dict-Merge 6:
Pick two dictionaries I D , I D from S h 7:

I = I D ∪ I D 8:
I D,D = Dict-Update(I) using Eq. 5 end if 13: end for 14: Return I D , the last dictionary in S k k T t,i S. As a result, for any new point x t we need to evaluate K(x s , x t ) only for the indices s in I. Therefore, SQUEAK never performs more than n(3qd eff (γ) n ) 2 kernel evaluations, which means that it does not even need to observe large portions of the kernel matrix. Finally, the runtime is dominated by the n matrix inversions used in Eq. 4. Therefore, the total runtime is O(n max t=1,...,n |I t |

3 ) = O(nd eff (γ) 3 n q 3 ). In the next section, we introduce DISQUEAK, which improves the runtime by independently constructing separate dictionaries in parallel and then merging them recursively to construct a final ε-accurate dictionary.

Distributed RLS Sampling

In this section, we show that a minor change in the structure of SQUEAK allows us to parallelize and distribute the computation of the dictionary I n over multiple machines, thus reducing even further its time complexity. Beside the computational advantage, a distributed architecture is needed as soon as the input dimension d and the number of points n is so large that having the dataset on a single machine is impossible. Furthermore, distributed processing can reduce contention on bottleneck data sources such as databases or network connections. DIstributed-SQUEAK (DISQUEAK, Alg. 2) partitions D over multiple machines and the (small) dictionaries that are generated from different portions of the dataset are integrated in a hierarchical way. The initial dataset is partitioned over k disjoint sub-datasets D i with i = 1, . . . , k and k dictionaries I Di = {(j, p 0,i = 1, q 0,i = q) : j ∈ D i } are initialized simply by placing all samples in D i into I with weight 1 and multiplicity q. Alternatively, if the datasets D i are too large to fit in memory, we can run SQUEAK to generate the initial dictionaries. The dictionaries I Di are added to a dictionary collection S, and following any predefined merge tree as in Fig. 1, these dictionaries are progressively merged together with Dict-Merge. Given two input dictionaries I D and I D , we combine them into a single dictionary I (the equivalent of the Expand phase in SQUEAK) and then Dict-Update is run on the merged dictionaries to create an updated dictionary I D∪D , which is placed back in the dictionary collection S. Since Dict-Merge only takes the two dictionaries as input and does not require any information on the dictionaries in the rest of the tree, separate branches can be run simultaneously on different machines, and only the resulting (small) dictionary needs to be propagated to the parent node for the future Dict-Merge. Unlike in SQUEAK, Dict-Update is run on the union of two distinct dictionaries rather than one dictionary and a new single point. As a result, we need to derive the "distributed" counterparts of Lemmas 1 and 2 to analyze the behavior of the RLSs and the quality of the estimator.

I {4,1} I {3,1} + I {3,2} I {3,1} I {2,1} + I {2,2} I {2,1} I {1,1} D 1 I {2,2} I {1,2} + I {1,3} I {1,2} D 2 I {1,3} D 3 I {3,2} I {2,3} I {2,3} I {1,4} + I {1,5} I {1,4} D 4 I {1,5} D 5 h = 4 h = 3 h = 2 h = 1
Lemma 3. Given two disjoint datasets D, D , for ev- ery i ∈ D ∪ D , τ i,D ≥ τ i,D∪D and 2d eff (γ) D∪D ≥ d eff (γ) D + d eff (γ) D ≥ d eff (γ) D∪D .
While in SQUEAK we were merging an ε-accurate dictionary I t and a new point, which is equivalent to a perfect, 0-accurate dictionary, in DISQUEAK both dictionaries used in a merge are only ε-accurate. To compensate for this loss in accuracy, we introduce a new estimator,

τ D∪D ,i = 1-ε γ (k i,i -k T i S(S T KS + (1 + ε)γI) -1 S T k i ), (5)
where S is the selection matrix associated with the temporary dictionary I = I D ∪ I D . Eq. 5 has similar guarantees as Lem. 2, with only a slightly larger α.

Lemma 4. Given two disjoint datasets D, D , and two ε-approximate dictionaries I D , I D , let I = I D ∪ I D and S be the associated selection matrix. Let K be the kernel matrix computed on D ∪ D , k i its i-th column, and τ D∪D ,i the RLS of k i . Then for all i in I such that q i = 0, the estimator in Eq. 5 is α-accurate, i.e. it satisfies

τ D∪D ,i /α ≤ τ D∪D ,i ≤ τ D∪D ,i , with α = (1 -ε)/(1 + 3ε).
Given these guarantees, the analysis of DISQUEAK follows similar steps as SQUEAK. Given ε-accurate dictionaries, we obtain α-accurate RLS estimates τ D∪D ,i that can be used to resample all points in I and generate a new dictionary I D,D that is ε-accurate. To formalize a result equivalent to Thm. 1, we introduce additional notation: We index each node in the merge tree by its height h and position l. We denote the dictionary associated to node {h, l} by I {h,l} and the collection of all dictionaries available at height h of the merge tree by S h = {I {h,l} }. We also use K {h,l} to refer to the kernel matrix constructed from the datasets D {h,l} , which contains all points present in the leaves reachable from node {h, l}. For instance in Fig. 1, node {3, 1} is associated with I {3,1} , which is an εapproximate dictionary of the kernel matrix Theorem 2. Let ε > 0 be the accuracy parameter, γ > 1 the regularization factor, and 0 < δ < 1 the prob. of failure. Given an arbitrary dataset D and a merge tree structure of height k as input together with parameters ε, γ, and δ, we run DISQUEAK with

K {3,1} constructed
q = 39α log (2n/δ) ε 2 • where α = (1 + 3ε)/(1 -ε). Then, w.p. at least 1 -δ,
DISQUEAK generates a sequence of collections of dictionaries {S h } k h=1 such that each dictionary I {h,l} in S h is ε-accurate (Eq. 1) w.r.t. to K {h,l} , and that at any node l of height h the size of the dictionary is bounded as |I {h,l} | ≤ 3qd eff (γ) {h,l} . The cumulative (across nodes) space and time requirementsof the algorithm depend on the exact shape of the merge tree.

Theorem 2 gives approximation and space guarantees for every node of the tree. In other words, it guarantees that each intermediate dictionary processed by DISQUEAK is both an ε-accurate approximation of the datasets used to generate it, and requires a small space proportional to the effective dimension of the same dataset. From an accuracy perspective, DISQUEAK provides exactly the same guarantees of SQUEAK. Analysing the complexity of DISQUEAK is however more complex, since the order and arguments of the Dict-Update operations is determined by the merge tree. We distinguish between the time and work complexity of a tree by defining the time complexity as the amount of time necessary to compute the final solution, and the work complexity as the total amount of operations carried out by all machines in the tree in order to compute the final solution. We consider two special cases, a fully balanced tree (all inner nodes have either two inner nodes as children or two leaves), and a fully unbalanced tree (all inner nodes have exactly one inner node and one leaf as children). For both cases, we consider trees where each leaf dataset contains a single point D i = {x i }.

In the fully unbalanced tree, we always merge the current dictionary with a new dataset (a single new point) and no Dict-Merge operation can be carried out in parallel. Unsurprisingly, the sequential algorithm induced by this merge tree is strictly equivalent to SQUEAK. Computing a solution in the fully unbalanced tree takes O(nd eff (γ) 3 D q 3 ) time with a total work that is also O(nd eff (γ) 3 D q 3 ), as reported in Thm. 1. On the opposite end, the fully balanced tree needs to invert a d eff (γ) {h,l} dimensional matrix at each layer of the tree for a total of log(n) layers. Bounding all d eff (γ) {h,l} with d eff (γ) D , gives a complexity for computing the final solution of O(log(n)q 3 d eff (γ) 3 D ) time, with a huge improvement on SQUEAK. Surprisingly, the total work is only twice O(nq 3 d eff (γ) 3 D ), since at each layer h we perform n/2 h inversions (on n/2 h machines), and the sum across all layers is log(n) h=1 n/2 h ≤ 2n. Therefore, we can compute a solution in a much shorter time than SQUEAK, with a comparable amount of work, but at the expense of requiring much more memory across multiple machines, since at layer h, the sum

|S h | l=1 d eff (γ) {h,l} can be much larger than d eff (γ) D .
Nonetheless, this is partly alleviated by the fact that each node {h, l} locally requires only d eff (γ) 2 {h,l} ≤ d eff (γ) 2 D memory. Proof sketch: Although DISQUEAK is conceptually simple, providing guarantees on its space/time complexity and accuracy is far from trivial. The first step in the proof is to carefully decompose the failure event across the whole merge tree into separate failure events for each merge node {h, l}, and for each node construct a random process Y that models how Alg. 2 generates the dictionary I {h,l} . Notice that these processes are sequential in nature and the various steps (layers in the tree) are not i.i.d. Furthermore, the variance of Y is potentially large, and cannot be bounded uniformly. Instead, we take a more refined approach, inspired by Pachocki [START_REF] Pachocki | Analysis of resparsification[END_REF], that 1) uses Freedman's in-equality to treat W, the variance of process Y, as a random object itself, 2) applies a stochastic dominance argument to W to reduce it to a sum of i.i.d. r.v. and only then we can 3) apply i.i.d. concentrations to obtain the desired result.

Applications

In this section, we show how our approximation guarantees translate into guarantees for typical kernel methods. As an example, we use kernel ridge regression. We begin by showing how to get an accurate approximation K n from an ε-accurate dictionary. Lemma 5. Given an ε-accurate dictionary I t of matrix K t , and the selection matrix S t , the regularized Nyström approximation of K t is defined as

K n = K n S n (S T n K n S n + γI m ) -1 S T n K n , (6)
and satisfies

0 K t -K t γ 1 -ε K t (K t + γI) -1 γ 1 -ε I. (7)
This is not the only choice of an approximation from dictionary I n . For instance, Musco and Musco [START_REF] Musco | Provably useful kernel matrix approximation in linear time[END_REF] show similar result for an unregularized Nyström approximation and Rudi et al. [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF] for a smaller K n , construct the estimator only for the points in

I n . Let m = |I n |, W = (S T n K n S n +γI m ) ∈ R m×m is nonsingu- lar and C = K n S n ∈ R n×m ,
applying the Woodbury formula we compute the regression weights as

w n =( K n + µI n ) -1 y n = (CW -1 C T + µI n ) -1 y n = 1 µ y n -C C T C + µW -1 C T y n . (8)
Computing (C T C + µW) -1 takes O(nm 2 ) time to construct the matrix and O(m 3 ) to invert it, while the other matrix-matrix multiplication take at most O(nm 2 ) time. Overall, these operations require to store at most an n × m matrix. Therefore the final complexity of computing a KRR using the dictionary is reduced from O(n 3 ) to O(nm 2 + m 3 ) time, and from O(n 2 ) to O(nm) space. The following corollary provides guarantees for the empirical risk of the solution w n in a fixed design setting.

Corollary 1 ([1, Thm. 3]). For an arbitrary dataset D, let K be the kernel matrix constructed on D. Run SQUEAK or DISQUEAK with regularization parameter γ. Then, the solution w computed using the regularized Nyström approximation K satisfies

R D ( w) ≤ 1 + γ µ 1 1 -ε 2 R D ( w),
where µ is the regularization of kernel ridge regression and R D ( w) is the empirical risk on D. For the random design setting, [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF] provides a similar bound for batch RLS sampling, under some mild assumption on the kernel function K(•, •) and the dataset D. Deriving identical guarantees for SQUEAK and DISQUEAK is straightforward.

Time |I n | (Total space = O(n|I n |)) Increm. Exact n 3 n - Bach [2] nd max 2 n ε + d max 3 n ε d max,n ε No Alaoui and Mahoney [1] n(|I n |) 2 λ min + nµε λ min -nµε d eff (γ) n + Tr(K n ) µε No Calandriello et al. [3] λ 2 max γ 2 n 2 d eff (γ) 2 n ε 2 λ max γ d eff (γ) n ε 2 Yes SQUEAK n 2 d eff (γ) 2 n ε 2 d eff (γ) n ε 2 Yes RLS-sampling nd eff (γ) 2 n ε 2 d eff (γ) n ε 2 -
Other applications. The projection matrix P n naturally appears in some form in nearly all kernel-based methods. Therefore, in addition to KRR in the fixed design [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] and in the random design [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF], any kernel matrix approximation that provides ε-accuracy guarantees on P n can be used to provide guarantees for a variety of kernel-based methods. As an example, Musco and Musco [START_REF] Musco | Provably useful kernel matrix approximation in linear time[END_REF] show this is the case for kernel PCA [START_REF] Schölkopf | Kernel principal component analysis[END_REF], kernel CCA with regularization, and kernel K-means clustering.

Discussion

Tab. 1 compares several kernel approximation methods w.r.t. their space and time complexity. For all methods, we omit O(log(n)) factors. We first report RLS-sampling, a fictitious algorithm that receives the exact RLSs as input, as an ideal baseline for all RLS sampling algorithms. The space complexity of uniform sampling [START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF] scales with the maximal degree of freedom d max . Since

d max = n max i τ n,i ≥ i τ n,i = d eff (γ) n ,
uniform sampling is often outperformed by RLS sampling. While Alaoui and Mahoney [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] also sample according to RLS, their two-pass estimator does not preserve the same level of accuracy. In particular, the first pass requires to sample O (nµε/(λ minnµε)) columns, which quickly grows above n 2 when λ min becomes small. Finally, Calandriello et al. [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF] require that the maximum dictionary size is fixed in advance, which implies some information about the effective dimensions d eff (γ) n , and requires estimating both τ t,i and d eff (γ) t . This extra estimation effort causes an additional λ max /γ factor to appear in the space complexity. This factor cannot be easily estimated, and it leads to a space complexity of n3 in the worst case. Therefore, we can see from the table that SQUEAK achieves the same space com-plexity (up to constant factors) as knowing the RLS in advance and hence outperforms previous methods.

A recent method by Musco and Musco [START_REF] Musco | Provably useful kernel matrix approximation in linear time[END_REF] achieves comparable space and time guarantees as SQUEAK. 3 While they rely on a similar estimator, the two approaches are very different. Their method is batch in nature, as it processes the whole dataset at the same time and it requires multiple passes on the data for the estimation of the leverage scores and the sampling. On the other hand, SQUEAK is intrinsically sequential and it only requires one single pass on D, as points are "forgotten" once they are dropped from the dictionary. Furthermore, the different structure requires remarkably different tools for the analysis. While the method of Musco and Musco [START_REF] Musco | Provably useful kernel matrix approximation in linear time[END_REF] can directly use i.i.d. concentration inequalities (for the price of needing several passes), we need to rely on more sophisticated martingale arguments to consider the sequential stochastic process of SQUEAK. Furthermore, DISQUEAK smoothly extends the ideas of SQUEAK over distributed architectures, while it is not clear whether [START_REF] Musco | Provably useful kernel matrix approximation in linear time[END_REF] could be parallelized or distributed.

Future developments Both SQUEAK and DIS-QUEAK need to know in advance the size of the dataset n to tune q. An interesting question is whether it is possible to adaptively adjust the q parameter at runtime. This would allow us to continue updating I n and indefinitely process new data beyond the initial dataset D. It is also interesting to see whether SQUEAK could be used in conjuntion with existing meta-algorithms (e.g., [START_REF] Kumar | Sampling methods for the Nyström method[END_REF] with model averaging) for kernel matrix approximation that can leverage an accurate sampling scheme as a black-box, and what kind of improvements we could obtain.

Notation summary

Kernel matrix at time

t K t ∈ R t×t , i-th column k t,i ∈ R t , (i, j)-th entry k i,j ∈ R Eigendecomposition K t = U t Λ t U T t ∈ R t×t , eigenvector matrix U t ∈ R t×t and diagonal SDP eigenvalues matrix Λ t ∈ R t×t .
Kernel matrix at time t in the RKHS,

K t = Φ T t Φ t , with Φ t ∈ R D×t SVD decomposition Φ t = V t Σ t U T
t ∈ R D×t with left singular vectors, V t ∈ R D×D , right singular vector U t (same as K t ), and singular values Σ t ∈ R D×t (σ i on the main diagonal and zeros under it) SVD decomposition Φ t = V t Σ t U T t ∈ R D×t with left singular vectors, V t ∈ R D×D , right singular vector U t (same as K t ), and singular values Σ t ∈ R D×t (σ i on the main diagonal and zeros under it)

P t = Ψ t Ψ T t with Ψ t = (K t + γI t ) -1/2 K 1/2 t ∈ R t×t Column dictionary at time t, I t = {(i, p t,i , q t,i )} t i=1 Selection matrix S t ∈ R t×t with { qt,i q pt,i } t i=1 on the diagonal P t = Ψ t S t S T t Ψ T t with S t ∈

A Preliminaries

In this section, we introduce standard matrix results and equivalent definitions for the kernel matrix and the the projection error which use convenient representations exploiting the feature space.

Matrix identity. We often use the following identity.

Proposition 2. For any symmetric matrix A ∈ R n×m and any γ > 0

A(A T A + γI m ) -1 A T = AA T (AA T + γI n ) -1 .
For any symmetric matrix A ∈ R n×n and diagonal matrix B ∈ R n×n such that B has ns zero entries, and s non-zero entries, define C ∈ R n×s as the matrix obtained by removing all zero columns in B. Then

AB(BAB + γI m ) -1 BA = AC(C T AC + γI s ) -1 C T A.
For any appropriately shaped matrix A, B, C, with A and B invertible, the Woodbury matrix identity states

(A + CBC T ) -1 = A -1 -A -1 C C T A -1 C + B -1 -1 C T A -1 .
Kernel matrix. Since K t is real symmetric matrix, we can eigendecompose it as

K t = U t Λ t U T t ,
where U t ∈ R t×t is the eigenvector matrix and Λ t ∈ R t×t is the diagonal eigenvalue matrix, with all non-negative elements since K t is PSD. Considering the feature mapping from X to the RKHS, we can write the kernel matrix as K t = Φ T t Φ t , where Φ t ∈ R D×t is the feature matrix, whose SVD decomposition is

Φ t = V t Σ t U T t ,
where V t ∈ R D×D contains the left singular vectors, Σ t ∈ R D×t has the singular values σ i on the main diagonal (followed by zeros), and U t contains the right singular vector, which coincide with the eigenvectors of K t . Furthermore, we have Λ t = Σ T t Σ t , and each eigenvalue λ i = σ i , with i = 1, . . . , t.

Projection error. We derive a convenient lemma on the formulation of the projection error both in terms of kernels and feature space. Lemma 6. The following identity holds.

P t -P t 2 = (K t + γI t ) -1/2 K 1/2 t (I t -S t S T t )K 1/2 t (K t + γI t ) -1/2 2 = (Φ t Φ T t + γI D ) -1/2 Φ t (I t -S t S T t )Φ T t (Φ t Φ T t + γI D ) -1/2 2 .
Proof of Lemma 6. Using the SVD decomposition,

(K t + γI t ) -1/2 K 1/2 t (I t -S t S T t )K 1/2 t (K t + γI t ) -1/2 2 = (U t Λ t U T + γU t I t U T t ) -1/2 U t Λ 1/2 t U T (I t -S t S T t )U t Λ 1/2 t U T (U t Λ t U T + γU t I t U T t ) -1/2 2 = U t (Λ t + γI t ) -1/2 Λ 1/2 t U T (I t -S t S T t )U t Λ 1/2 t (Λ t + γI t ) -1/2 U T 2 = U t (Λ t + γI t ) -1/2 Λ 1/2 t U T (I t -S t S T t )U t Λ 1/2 t (Λ t + γI t ) -1/2 U T 2 = (Λ t + γI t ) -1/2 Λ 1/2 t U T (I t -S t S T t )U t Λ 1/2 t (Λ t + γI t ) -1/2 2 = (Σ t Σ T t + γI D ) -1/2 Σ t U T (I t -S t S T t )U t Σ T t (Σ t Σ T t + γI D ) -1/2 2 = V t (Σ t Σ T t + γI D ) -1/2 Σ t U T (I t -S t S T t )U t Σ T t (Σ t Σ T t + γI D ) -1/2 V T t 2 = V t (Σ t U t U T t Σ T t + γI D ) -1/2 V T t V t Σ t U T (I t -S t S T t )U t Σ T t V T V(Σ t U t U T t Σ T t + γI D ) -1/2 V T t 2 = (Φ t Φ T t + γI D ) -1/2 Φ t (I t -S t S T t )Φ T t (Φ t Φ T t + γI D ) -1/2 2 .
B Ridge Leverage Scores and Effective Dimension (Proof of Lemma 1 and 3)

From Calandriello et al. [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF]Lem. 1] we know that τ t,i ≤ τ t-1,i and d eff (γ) t ≥ d eff (γ) t-1 . We will now prove the lower bound τ t,i ≥ τ t-1,i /(τ t-1,1 + 1).

Considering the definition of τ t,i in terms of φ i and Φ t , and applying the Sherman-Morrison formula we obtain

τ t,i = φ T i (Φ t Φ T t + γI) -1 φ i = φ T i (Φ t-1 Φ T t-1 + φ t φ T t + γI) -1 φ i = φ T i (Φ t-1 Φ T t-1 + γI) -1 φ i - φ T i (Φ t-1 Φ T t-1 + γI) -1 φ t φ T t (Φ t-1 Φ T t-1 + γI) -1 φ i 1 + φ T t (Φ t-1 Φ T t-1 + γI) -1 φ t = τ t-1,i - φ T i (Φ t-1 Φ T t-1 + γI) -1 φ t φ T t (Φ t-1 Φ T t-1 + γI) -1 φ i 1 + φ T t (Φ t-1 Φ T t-1 + γI) -1 φ t • Let x = (Φ t-1 Φ T t-1 + γI) -1/2 φ i and y = (Φ t-1 Φ T t-1 + γI) -1/2 φ t .
Then τ t,i /τ t-1,i is equal to

τ t,i τ t-1,i = 1 - (φ T t (Φ t-1 Φ T t-1 + γI) -1 φ i ) 2 (1 + φ T t (Φ t-1 Φ T t-1 + γI) -1 φ t )φ T i (Φ t-1 Φ T t-1 + γI) -1 φ i = 1 - (y T x) 2 (1 + y T y) x T x •
Defining the cosine between y and x as cos(y, x) = y T x/( x y ), we have that

1 - (y T x) 2 (1 + y T y) x T x = 1 - y T yx T xcos(y, x) 2 (1 + y T y) x T x = 1 - y 2 1 + y 2 cos(y, x) 2 ,
where y 2 1+ y 2 depends only on the norm of y and not its direction, and cos(y, x) depends only on the direction of y and is maximized when y = x. Therefore,

τ t+1,i τ t,i = 1 - (y T x) 2 (1 + y T y) x T x = 1 - y 2 1 + y 2 cos(y, x) 2 ≥ 1 - x 2 1 + x 2 = 1 1 + x 2 = 1 1 + τ t,i ,
which concludes the proof of Lem. 1.

For Lem. 3, the first point τ i,D ≥ τ i,D∪D can be easily proven by choosing D to construct a kernel matrix K D , and then invoke Lem. 1 as we add one sample at a time from D . Also as easily for d eff (γ) D + d eff (γ) D we have

d eff (γ) D + d eff (γ) D ≤ 2 max{d eff (γ) D ; d eff (γ) D } ≤ 2 max{d eff (γ) D∪D ; d eff (γ) D∪D } = 2d eff (γ) D∪D .
Finally, we prove the other side of the inequality for d eff (γ) D +d eff (γ) D . Let Φ D , Φ D be the matrices constructed using the feature vectors of the samples in D and D respectively. Then,

d eff (γ) D + d eff (γ) D = i∈D τ D,i + i∈D τ D ,i = i∈D φ T i (Φ D Φ T D + γI D ) -1 φ i + i∈D φ T i (Φ D Φ T D + γI D ) -1 φ i ≥ i∈D φ T i (Φ D∪D Φ T D∪D + γI D ) -1 φ i + i∈D φ T i (Φ D∪D Φ T D∪D + γI D ) -1 φ i = i∈D∪D φ T i (Φ D∪D Φ T D∪D + γI D ) -1 φ i = i∈D∪D τ D∪D ,i = d eff (γ) D∪D .
C Ridge Leverage Scores Estimation (Proof of Lemma 2 and 4)

We begin with a convenient reformulation of the ridge leverage scores,

τ t,i = e T t,i K t (K t + γI t ) -1 e t,i = e T t,i Φ T t Φ t (Φ T t Φ t + γI t ) -1 e t,i = e T t,i Φ T t (Φ t Φ T t + γI D ) -1 Φ t e t,i = φ T i (Φ t Φ T t + γI D ) -1 φ i .
This formulation, combined with Def. 1, suggests φ T i (Φ t S t S T t Φ T t + γI D ) -1 φ i as an estimator for τ t,i . However, at step t, we only have access to an ε-accurate dictionary w.r.t. Φ t-1 and not w.r.t. Φ t . Therefore, we augment it with (t, 1, q) to construct I t and the corresponding S t , which will have [S t ] t,t = 1. We will now show how to implement this estimator efficiently. From Prop. 2, we apply Woodbury matrix identity with A = γI, B = I and

C = Φ t S t , τ t,i = (1 -ε)φ T i (Φ t S t S T t Φ T t + γI D ) -1 φ i = (1 -ε)φ T i (Φ t S t I t S T t Φ T t + γI D ) -1 φ i (P rop. 2) = (1 -ε)φ T i 1 γ I D - 1 γ 2 Φ t S t 1 γ S T t Φ T t Φ t S t + I t -1 S T t Φ T t φ i = (1 -ε) γ φ T i I D -Φ t S t S T t Φ T t Φ t S t + γI t -1 S T t Φ T t φ i = (1 -ε) γ φ T i φ i -φ T i Φ t S t S T t Φ T t Φ t S t + γI t -1 S T t Φ T t φ i = (1 -ε) γ k i,i -k t,i S t (S T t K t S t + γI t ) -1 S T t k t,i ,
which is the estimator defined in Eq. 4.

We can generalize this estimator to the case where instead of using a single dictionary and fresh data to estimate τ t,i , we are using k ε-accurate dictionaries I k . Given disjoint datasets {D i } k i=1 with associated feature matrices Φ i . From each dataset, construct an ε-accurate dictionary I i , with its associated selection matrix S i .

To estimate the RLS τ i of point i w.r.t. the whole dataset D = ∪ k j=1 D j , and corresponding feature matrix Φ, we set the estimator to be

τ i = (1 -ε)φ T i   k j=1 Φ j S j S T j Φ T j + (1 + (k -1)ε)γI D   -1 φ i .
Part 1: accuracy of the RLS estimator τ i . Since each of the dictionaries I i used to generate S i is ε-accurate, we can use the equivalence from Lem. 6,

P i -P i 2 = (Φ i Φ T i + γI D ) -1/2 (Φ i Φ T i -Φ i S i S T i Φ T i )(Φ i Φ T i + γI D ) -1/2 2 ≤ ε, which implies that (1 -ε)Φ i Φ T i -εγI D Φ i S i S T i Φ T i (1 + ε)Φ i Φ T i + εγI D .
Therefore, we have

τ i = (1 -ε)φ T i   k j=1 Φ j S j S T j Φ T j + (1 + (k -1)ε)γI D   -1 φ i ≤ (1 -ε)φ T i   k j=1 (1 -ε)Φ j Φ T j -εγI D + (1 + (k -1)ε)γI D   -1 φ i ≤ (1 -ε)φ T i (1 -ε)ΦΦ T -kεγI D + (1 + (k -1)ε)γI D -1 φ i = (1 -ε)φ T i (1 -ε)(ΦΦ T + γI D ) -1 φ i = (1 -ε) (1 -ε) φ T i ΦΦ T + γI D -1 φ i = τ i , and
τ i = (1 -ε)φ T i   k j=1 Φ j S j S T j Φ T j + (1 + (k -1)ε)γI D   -1 φ i ≥ (1 -ε)φ T i   k j=1 (1 + ε)Φ j Φ T j + εγI D + (1 + (k -1)ε)γI D   -1 φ i = (1 -ε)φ T i (1 + ε)ΦΦ T + kεγI D + (1 + (k -1)ε)γI D -1 φ i = (1 -ε)φ T i (1 + ε)ΦΦ T + (1 + (2k -1)ε)γI D ) -1 φ i ≥ (1 -ε)φ T i (1 + (2k -1)ε)ΦΦ T + (1 + (2k -1)ε)γI D ) -1 φ i = (1 -ε) (1 + (2k -1)ε) φ T i ΦΦ T -γI D -1 φ i = (1 -ε) (1 + (2k -1)ε) τ i .
Then, we can instantiate this result with k = 1 to prove the accuracy claim in Lem. 2, and with k = 2 to prove the accuracy claim in Lem. 4.

Part 2: accuracy of min { τ t , τ t-1 }. To simplify the notation, for this part of the proof we indicate with τ t ≤ τ t-1 that for each i ∈ {1, . . . , t -1} we have τ t,i ≤ τ t-1,i . From Lem. 1, we know that τ t-1 ≥ τ t . Given α-accurate τ t and τ t-1 we have the upper bound and the lower bound,

min { τ t , τ t-1 } ≤ min {τ t , τ t-1 } = τ t ,
min { τ t , τ t-1 } ≥ 1 α min {τ t , τ t-1 } = 1 α τ t ,
which combined gives us 1 α τ t ≤ min { τ t , τ t-1 } ≤ τ t as required by the definition of α-accuracy.

D Proof of Thm. 1 and Thm. 2

From the discussion of Thm. 2, we know that running SQUEAK is equivalent to running DISQUEAK on a specific (fully unbalanced) merge tree. Therefore, we just prove Thm. 2, and invoke it on this tree to prove Thm. 1.

We begin by describing more in detail some notation introduced in the main paper and necessary for this proof.

Merge trees

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

We partition D into k disjoint sub-datasets D i of size n i , such that D = ∪ k i=1 D i . For each dataset D i , we construct an initial dictionary I {1,i} = {(j, p 0,i = 1, q 0,i = q) : j ∈ D i } by inserting all points from D i into I Di with weight p 0,i = 1 and number of copies q 0,i = q. It is easy to see that I {1,i} is an ε-accurate dictionary, and we can split the data in small enough chunks to make sure that it can be easily stored and manipulated in memory. Alternatively, if we want our initial dictionaries to be small and cannot choose the size of D i , we can run Alg. 1 on D i to generate I {1,i} , and the following proof will remain valid. Regardless of construction, the initial dictionaries I {1,i} are included into the dictionary pool S 1 .

At iteration h, the inner loop of Alg. 2 arbitrarily chooses two dictionaries from S h and merges them into a new dictionary. 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. Figure 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 dictionaries I 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 depth of the tree minus 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 dictionary merge, and the size of S h (number of dictionaries present at layer h) is

|S h | = k -h + 1.
Therefore, a node corresponding to a dictionary is uniquely identified with two indices {h, l}, where h is the height of the layer and l ≤ |S h | is the index of the node in the layer. For example, in Figure 2(a), the node containing I 1,2,3 is indexed as {3, 1}, and the highest node containing I 4 is indexed as {3, 2}.

We also define the dataset D {h,l} as the union of all sub-datasets D l that are reachable from node {h, l} as leaves. For example, in Fig. 2(a), dictionary I 1,2,3 in node {3, 1} is constructed starting from all points in D {3,1} = D 1 ∪ D 2 ∪ D 3 , where we highlight in red the descendant tree. We now define K h as the block diagonal kernel matrix where each diagonal block K {h,l} is constructed on D {h,l} . Again, from Fig. 2, K 3 is a n × n matrix with two blocks on the diagonal, a first (n

1 + n 2 + n 3 ) × (n 1 + n 2 + n 3 ) block K 3,1 constructed on D {3,1} = D 1 ∪ D 2 ∪ D 3 ,
and a second n 4 × n 4 block K 3,2 constructed on D {3,2} = D 4 . Similarly, we can adapt Def. 1 to define P h as a block diagonal projection matrix, where block P {h,l} is defined using K {h,l} , and block diagonal P h , where block P {h,l} is defined using K {h,l} and I {h,l} .

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 l P {h,l} -P {h,l} ≤ ε ⇒ P {h,l} -P {h,l} ≤ ε for all blocks l on the diagonal, and since each block corresponds to a dictionary I {h,l} , this means that if P h -P h ≤ ε, all dictionaries at layer l are ε-accurate approximation of their respective represented datasets. Our goal is to show

P ∃h ∈ {1, . . . , k} : P h -P h 2 ≥ ε ∪ max l=1,...,|S h | |I {h,l} | ≥ 3qd eff (γ) {h,l} = P ∃h ∈ {1, . . . , k} : max l=1,...,|S h | P {h,l} -P {h,l} 2 ≥ ε A h ∪ max l=1,...,|S h | |I {h,l} | ≥ 3qd eff (γ) {h,l} B h ≤ δ, (9)
where event A h refers to the case when some dictionary I {h,l} at an intermediate layer h fails to accurately approximate K {h,l} and event B h considers the case when the memory requirement is not met (i.e., too many points are kept in one of the dictionaries I {h,l} at a certain layer h). We can conveniently decompose the previous joint (negative) event into two separate conditions as

P k h=1 A h ∪ B h = P k h=1 A h ∪ k h=1 B h = P k h=1 A h + P   k h=1 B h ∩ k h=1 A h C   = P k h=1 A h + P k h=1 B h ∩ k h=1 A C h = P k h=1 A h + P k h=1 B h ∩ k h =1 A C h .
Applying this reformulation and a union bound we obtain

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

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 | l=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 that we can now analyze the processes that generated each P {h,l} -P {h,l} (and each dictionary I {h,l} ) separately. Focusing on a single node {h, l} restricts our problem on a well defined dataset D {h,l} , where we can analyze the evolution of I {h,l} sequentially.

Challenges. Due to its sequential nature, the "theoretical" structure of the process generating P {h,l} -P {h,l} is considerably complicated, where each step in the sampling process (layer in the merge tree) is highly correlated with the previous steps. This prevents us from using concentration inequalities for i.i.d. processes that are at the basis of the analysis of uniform sampling [START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF] and the method proposed by Alaoui and Mahoney [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF]. As a result, we first show that the project error is a martingale process. The main difficulty technical difficulty in analyzing how the projection error evolves over iterations is that the projection matrices change dimension every time a new point is processed. In fact, all of the matrices P {h ,l } -P {h ,l } descending from P {h,l} -P {h,l} have potentially a different size since they are based on different kernel matrices K {h ,l } . This requires a careful definition of the martingale process to still use matrix concentration inequalities for fixed-size matrices (see Sect. D.1). Another major technical challenge is to control the variance of the martingale increments. In fact, at each the projection error may increase by a quantity whose cumulative variance can be arbitrarily large. As a result, a direct use of the Freedman matrix inequality in Sect. D.2 would not return an accurate result. In order to provide a tighter bound on the total variance of the martingale process of the projection error, we need to introduce an i.i.d. stochastically dominant process (Sect. D.3 (step 2)), which finally allows us to use an i.i.d. matrix concentration inequality to bound the total variance (Sect. D.3-(step 5)). This finally leads to the bound on the accuracy. The bound on the space complexity (Sect. D.4)follows similar (but simpler) steps.

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

The sequential process. Thanks to the union bound in Eq. 10, instead of having to consider the whole merge tree followed by Alg. 2, we can focus on each individual node {h, l} and study the sequential process that generated its dictionary I {h,l} . We will now map more clearly the actions taken by Alg. 2 to the process that generated P {h,l} -P {h,l} . We begin by focusing on P {h,l} , which is a random matrix defined starting from the fixed kernel matrix K {h,l} and the random dictionary I {h,l} , where the randomness influences both which points are included in I {h,l} , and the weight with which they are added.

Consider now a point i ∈ D {h,l} . Since the starting datasets in the leaves are disjoint, there is a single path in the tree, with length h, from the leaves to {h, l}. This means that for all s < h, we can properly define a unique p s,i and q s,i associated with that point. More in detail, if at layer s point i is present in D {s,l } , it means that either (1) Alg. 2 used I {s,l } to compute p s,i , and p s,i to compute q s,i , or (2) at layer h Alg. 2 did not have any merge scheduled for point i, and we simply propagate p s,i = p s-1,i and q s,i = q s-1,i . Consistently with the algorithm, we initialize p 0,i = 1 and q 0,i = q.

Denote ν {h,l} = |D {h,l} | so that we can use index i ∈ [ν {h,l} ] to index all points in D {h,l} . Given the symmetric

matrix Ψ = (K {h,l} + γI) -1/2 K 1/2 {h,l} with its i-th column ψ i = (K {h,l} + γI) -1/2 K 1/2
{h,l} e ν {h,l} ,i , we can rewrite the projection matrix as that P {h,l} = ΨΨ T = ν {h,l} i=1 ψ i ψ T i . Note that

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

Taking into account that we are now considering a specific node {h, l}, we can drop the index from the dataset D {h,l} = D, RLS τ D {h,l} ,i = τ h,i , and size ν {h,l} = ν. Using this shorter notation, we can reformulate our objective as bounding P {h,l} -P {h,l} 2 = P {h,l} -P {h,l} h 2 , and reformulate the process as a sequence of matrices {Y s } h s=1 defined as

Y s = P {h,l} -P {h,l} s = 1 q ν i=1 q j=1 1 - z s,i,j p s,i ψ i ψ T i ,
where Y h = P {h,l} -P {h,l} h = P {h,l} -P {h,l} , and Y 1 = P {h,l} -P {h,l} 0 = 0 since p 0,i = 1 and q 0,i = q.

D.2 Bounding Y h

We transformed the problem of bounding P {h,l} -P {h,l} into the problem of bounding Y h , which we modeled as a random matrix process, connected to Alg. 2 by the fact that both algorithm and random process Y h make use of the same weight p s,i and multiplicities q s,i . There are two main issues in analyzing the process {Y s } h s=1 :

(1) While the overall algorithm may fail in generating an accurate dictionary at some intermediate iteration and yet return an accurate dictionary at the end, whenever one of the intermediate Y s is larger than ε we lose our guarantees for p s,i for the whole process, since an inaccurate p s,i that underestimates too much p s.i will influence all successive p s ,i through the minimum. To solve this, we consider an alternative (more pessimistic) process which is "frozen" as soon as it constructs an inaccurate dictionary. Freezing the probabilities at the first error gives us a process that fails more often, but that provides strong guarantees up to the moment of failure.

(2) While Y h = P {h,l} -P {h,l} h ≤ ε guarantees that I {h,l} is an ε-accurate dictionary of K {h,l} , knowing Y s = P {h,l} -P {h,l} s ≤ ε for s < h does not guarantee that all descendant I {s,l } are ε-accurate w.r.t. their K {s,l } . Nonetheless, we will show that Y s ≤ ε is enough to guarantee that the intermediate estimate p s,i computed in Alg. 2, and used as weights in P {h,l} s , are never too small. The frozen process. 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. 2, 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,i,j and p s,i,j . We recursively define the indicator ({0, 1}) random variable z s,i,j as z s,i,j = I u s,i,j ≤ p s,i,j p s-1,i,j z s-1,i,j , where u s,i,j ∼ U(0, 1) is a [0, 1] uniform random variable and p s,i,j is defined as

p s,i,j = p s,i I Y s-1 ≤ ε ∩ z s-1,i,j = 1 + p s-1,i,j I Y s-1 ≥ ε ∪ z s-1,i,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,i,j and enforcing that p s+1,i,j = p s,i,j when z s,i,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,i,j / p s,i 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, i, and j into a linear index r = s in the range [1, . . . , ν 2 q], obtained as r = {s, i, j} = (s -1)νq + (i -1)q + j. We also define the difference matrix as

X {s,i,j} = 1 q z s-1,i,j p s-1,i,j - z s,i,j p s,i,j ψ i ψ T i ,
which allows us to write the cumulative matrix as Y {s,i,j} = {s,i,j} r=1

X {s,i,j} where the checkpoints {s, ν, q}

correspond to Y s , Y {s,ν,q} = Y s = 1 q ν i=1 q j=1 1 - z s,i,j p s,i,j ψ i ψ T i .
Let F s be the filtration containing all the realizations of the uniform random variables u s,i,j up to the step s, that is F s = {u s ,i ,j , ∀{s , i , j } ≤ s}. Again, we notice that F s defines the state of the algorithm after completing iteration s because, unless a "freezing" happened, Alg. 2 and Y s flip coins with the same probability, and generate the same dictionaries. Since τ s,i and p s,i,j are computed at the beginning of iteration s using the dictionary I {s,l } (for some l 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 ,i,j up to z s-1,i,j for any i and j, we have that all the Bernoulli variables z s,i,j at iteration s are conditionally independent given F s-1 . In other words, we have that for any i , and j such that {s, 1, 1} ≤ {s, i , j } < s the following random variables are equal in distribution,

z s,i,j F {s,i ,j } = z s,i,j F {s-1,ν,q} ∼ B p s,i,j p s-1,i,j , (11)
and for any i , and j such that {s, 1, 1} ≤ {s, i , j } ≤ {s, ν, q} and s = {s, i , j } we have the independence

z s,i,j F {s-1,ν,q} ⊥ z s,i ,j F {s-1,ν,q} . (12)
While knowing that Y s ≤ ε is not sufficient to provide guarantees for the approximate probabilities p s,i , we can show that it is enough to prove that the frozen probabilities p s,i,j are never too small. Lemma 7. Let α = (1 + 3ε)/(1ε) and p s,i,j be the sequence of probabilities generated by the freezing process. Then for any s, i, and j, we have p s,i,j ≥ p h,i /α = τ h,i /α.

Proof of Lemma 7. 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,i,j , we have that

p s,i,j ≥ p s,i = p s,i = max {min { τ s,i , p s-1,i } , p s-1,i /2} ≥ min { τ s,i , p s-1,i } = min { τ s,i , p s-2,i } = min { τ s,i , p s-3,i } . . . = min { τ s,i , p 0,i } = τ s,i
, and therefore p s,i,j ≥ τ s,i . Now let {s, l } be the node where τ s,i was computed. We will again drop the {h, l} index from D {h,l} , and simply refer to it as D. Similarly, we will refer with D i to D {s,l } (as in, the dataset used to compute τ s,i ), and with D i to the samples in D not contained in D i (complement of D i ). Define A as the |D| × |D i | matrix that contains the columns of S s related to points in D i , and similarly define B as the |D| × |D i | matrix that contains the columns of S s related to points in D i , where S s can be reconstructe by interleaving columns of A and B. From its definition in Eq.5, we know that τ s,i is computed by Alg. 2 as

τ s,i = (1 -ε)φ T i Φ D AA T Φ T D + (1 + ε)γI D -1 φ i ,
using only the points in A that are available at node {s, l }. From Lem.6 we know that

Y s = P {h,l} -P {h,l} s 2 = (Φ D Φ T D + γI D ) -1/2 (Φ D Φ T D -Φ D S s S T s Φ T D )(Φ D Φ T D + γI D ) -1/2 2 = (Φ D Φ T D + γI D ) -1/2 (Φ D Φ T D -Φ D AA T Φ T D -Φ D BB T Φ T D )(Φ D Φ T D + γI D ) -1/2 2 ≤ ε
and we know that this implies

Φ D AA T Φ T D Φ D Φ T D + ε(Φ D Φ T D + γI D ) -Φ D BB T Φ T D Φ D T D + ε(Φ D Φ T D + γI D ).
Plugging it in the initial definition,

τ s,i = (1 -ε)φ T i Φ D AA T Φ T D + (1 + ε)γI D -1 φ i ≥ (1 -ε)φ T i (Φ D Φ T D + ε(Φ D Φ T D + γI D ) + (1 + ε)γI D ) -1 φ i = (1 -ε) 1 1 + 2ε φ T i (Φ D Φ T D + γI D ) -1 φ i ≥ 1 -ε 1 + 2ε τ h,i ≥ τ h,i /α.
This result is weaker then Lemmas 2 and 4, since we do not provide an upper bound, and only show that p s,i ≥ p h,i /α and not p s,i ≥ p s,i /α, but it guarantees that the probabilities used at any intermediate layer s are bigger than a fraction 1/α of the exact probabilities that we would use at layer h, which will suffice for our purpose.

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,i,j} F s-1 = 0.
We begin by inspecting the conditional random variable X {s,i,j} |F s-1 . Given the definition of X {s,i,j} , the conditioning on F s-1 determines the values of z s-1,i,j and the approximate probabilities p s-1,i,j and p s,i,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,i,j} is the variable z s,i,j . Specifically, if Y s-1 ≥ ε, then we have p s,i,j = p s-1,i,j and z s,i,j = z s-1,i,j (the process is stopped), and the martingale requirement

E X {s,i,j} F s-1 = 0 is trivially satisfied. On the other hand, if Y s-1 ≤ ε we have E us,i,j 1 q z s-1,i,j p s-1,i,j - z s,i,j p s,i,j ψ i ψ T i F s-1 = 1 q z s-1,i,j p s-1,i,j - z s-1,i,j p s,i,j E I u s,i,j ≤ p s,i,j p s-1,i,j F s-1 ψ i ψ T i = 1 q z s-1,i,j p s-1,i,j - z s-1,i,j p s,i,j p s,i,j p s-1,i,j ψ i ψ T i = 0,
where we use the recursive definition of z s,i,j and the fact that u s,i,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,i,j} = 1 q z s-1,i,j p s-1,i,j - z s,i,j p s,i,j ψ i ψ T i ≤ 1 q 1 p s,i,j ψ i ψ T i = 1 q 1 p s,i,j τ h,i ≤ 1 q α τ h,i τ h,i = α q def = R,
where we used Lemma 7 to bound p s,i,j ≤ τ h,i /α. 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] Aaron | Freedman's inequality for matrix martingales[END_REF] to bound the norm of Y.

Proposition 3 (Tropp [START_REF] Aaron | Freedman's inequality for matrix martingales[END_REF], 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, . . . .
Define the predictable quadratic variation process of the martingale as

W k def = 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,i,j} 2 , we can directly apply Proposition 3 to bound (a) for any fixed σ 2 . To bound the part (b), we use the following lemma, proved later in Sec. D.3.

Lemma 8 (Low probability of the large norm of the predictable quadratic variation process).

P W h ≥ 6α q ≤ n • exp -2 q α
Combining Prop. 3 with σ 2 = 6α/q, Lem 8, the fact that 2ε/3 ≤ 1 and the value used by Alg. 2 q = 39α log(2n/δ)/ε 2 we obtain

P P {h,l} -P {h,l} 2 ≥ ε = P ( Y h ≥ ε) ≤ P Y h ≥ ε ∩ W h ≤ σ 2 + P W h ≥ σ 2 ≤ 2ν • exp - ε 2 q α 1 12 + 2ε/3 + n • exp -2 q α ≤ 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 since at most we can split our dataset in n parts, concludes this part of the proof.

D.3 Proof of Lemma 8 (bound on predictable quadratic variation)

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

W r = 1 q 2 {s,i,j}≤r E z s-1,i,j p s-1,i,j - z s,i,j p s,i,j 2 F {s,i,j}-1 ψ i ψ T i ψ i ψ T i .
We rewrite the expectation terms in the equation above as

E z s-1,i,j p s-1,i,j - z s,i,j p s,i,j 2 F {s,i,j}-1 = E z 2 s-1,i,j p 2 s-1,i,j -2 z s-1,i,j p s-1,i,j z s,i,j p s,i,j + z 2 s,i,j p 2 s,i,j
F {s,i,j}-1

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

W h = W {h,m,q} = 1 q 2 q j=1 ν i=1 h s=1 z s-1,i,j p s-1,i,j 1 p s,i,j - 1 p s-1,i,j ψ i ψ T i ψ i ψ T i .
We can now upper-bound W h as

W h 1 q 2 q j=1 ν i=1 h s=1 z s-1,i,j p s-1,i,j 1 p s,i,j - 1 p s-1,i,j ψ i ψ T i ψ i ψ T i = 1 q 2 q j=1 ν i=1 z h,i,j p 2 h,i,j - z h,i,j p 2 h,i,j + h s=1 z s-1,i,j p s-1,i,j 1 p s,i,j - 1 p s-1,i,j ψ i ψ T i ψ i ψ T i = 1 q 2 q j=1 ν i=1 z h,i,j p 2 h,i,j + h s=1 - z s,i,j p 2 s,i,j + z s-1,i,j p s,i,j p s-1,i,j - z 0,i,j p 2 0,i,j ψ i ψ T i ψ i ψ T i 1 q 2 q j=1 ν i=1 z h,i,j p 2 h,i,j + h s=1 z s-1,i,j p s,i,j p s-1,i,j - z s,i,j p s,i,j p s-1,i,j ψ i ψ T i ψ i ψ T i = 1 q 2 q j=1 ν i=1 z h,i,j p 2 h,i,j + h s=1 z s-1,i,j (1 -z s,i,j ) p s,i,j p s-1,i,j ψ i ψ T i ψ i ψ T i ,
where in the inequality we use the fact p s,i,j ≤ p s-1,i,j . From the definition of p s,i,j , we know that when z s,i,j = 0, p s,i,j = p s-1,i,j . Therefore zs-1,i,j (1-zs,i,j ) p s,i,j p s-1,i,j = zs-1,i,j (1-zs,i,j )

p 2 s-1,i,j
, since the term is non-zero only when z s,i,j = 0.

Finally, we see that only one of the z s-1,i,j (1z s,i,j ) terms can be active for s ∈ [h] and thus

W h 1 q 2 q j=1 ν i=1 z h,i,j p 2 h,i,j + h s=1 z s-1,i,j (1 -z s,i,j ) p 2 s-1,i,j ψ i ψ T i ψ i ψ T i = 1 q 2 q j=1 ν i=1 max max s=1,...,h z s-1,i,j (1 -z s,i,j ) p 2 s-1,i,j ; z h,i,j p 2 h,i,j ψ i ψ T i ψ i ψ T i = 1 q 2 q j=1 ν i=1 ψ i ψ T i ψ i ψ T i max s=0,...,h z s,i,j p 2 s,i,j • (13)
Step 2 (introduction of a stochastically dominant process). We want to study max s=0,...,h zs,i,j p 2 s,i,j

.

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

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

Unfortunately, due to the dependence between different copies of the point 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,i,j /p s,i,j } h s=0 with a set of (independent) random variables w 0,i,j that will stochastically dominate them.

We define the random variable w s,i,j using the following conditional distribution,4

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

1
for α/p h,i ≤ a .

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

P 1 w s,i,j ≤ a F s = P 1 w s,i,j ≤ a F s-1 .
Notice that in the definition of w s,i,j , none of the other w s ,i ,j (for any different s , i , or j ) appears and p s,i,j is a function of F s-1 . It follows that given F s-1 , w s,i,j is independent from all other w s ,i ,j (for any different s , i , 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,i,j the definition above reduces to

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

{z 0,i,j } {p 0,i,j } {w 0,i,j } {u 1,i,j } {z 1,i,j } {p 1,i,j } {w 1,i,j } {u 2,i,j } {z 2,i,j } {p 2,i,j } {w 2,i,j } {u 3,i,j } {z 3,i,j } {p 3,i,j } {w 3,i,j } {u 4,i,j } {z 4,i,j } {p 4,i,j } {w 4,i,j } {u 5,i,j } {z 5,i,j } {p 5,i,j } {w 5,i,j }
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,i,j p s,i,j ≤ a ≥ P 1 w 0,i,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,i,j , we know that, w.p.

1 1/ pk-1 1/ pk α 2 /p h 0 1-pk pk-1 1-p h α 2 pk 1 P max zk-1 pk-1 ; zk w k ≤a | F k-1 P zk-1 w k-1 ≤a|F k-1 zk-1 =1
Figure 4: C.d.f. of max z k-1,i,j /p t,i,j ; z k,i,j /p k,i,j and z k-1,i,j /w k-1,i,j conditioned on F {k-1} . For conciseness, we omit the i, j indices.

given F k-1 , we have

P z k,i,j w k,i,j ≤ a F k-1 = P I u k,i,j ≤ p k,i,j p k-1,i,j 1 w k,i,j ≤ a F k-1 =            0 for a ≤ 0 1 - p k,i,j p k-1,i,j for 0 ≤ a < 1/p k,i,j 1 - p k,i,j p k-1,i,j + p k,i,j p k-1,i,j 1 -1 p k,i,j a = 1 - 1 p k-1,i,j a for 1/p k,i,j ≤ a < α/p h,i 1 for α/p h,i ≤ a ≥            0 for a < 1/p k-1,i,j 1 - 1 p k-1,i,j a for 1/p k-1,i,j ≤ a < 1/p k,i,j 1 - 1 p k-1,i,j a for 1/p k,i,j ≤ a < α/p h,i 1 for α/p h,i ≤ a (15) = P 1 w k-1,i,j ≤ a F k-2 = P 1 w k-1,i,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,i,j = 1 (note that the peeling is trivially true when z k-1,i,j = 0 since the whole chain terminated at step z k-1,i,j ). We show how to move from an iteration k ≤ h to k -1.

P max max

s=0...k-1 z s,i,j p s,i,j ; z k,i,j w k,i,j ≤ a = E F k-1 P max max s=0...k-1 z s,i,j p s,i,j ; z k,i,j w k,i,j ≤ a F k-1 (a) ≥ E F k-1 P max max s=0...k-1 z s,i,j p s,i,j ; z k-1,i,j w k-1,i,j ≤ a F k-1 = E F k-1 P max max s=0...k-2 z s,i,j p s,i,j ; z k-1,i,j p k-1,i,j ; z k-1,i,j w k-1,i,j ≤ a F k-1 = E F k-1 P max max s=0...k-2 z s,i,j p s,i,j ; z k-1,i,j max 1 p k-1,i,j ; 1 w k-1,i,j ≤ a F k-1 (b) = E F k-1 P max max s=0...k-2 z s,i,j p s,i,j ; z k-1,i,j w k-1,i,j ≤ a F k-1 = P max max s=0...k-2 z s,i,j p s,i,j ; z k-1,i,j w k-1,i,j ≤ a ,
where in (a), given F k-1 , everything is fixed except u k,i,j and w k,i,j and we can use the stochastic dominance in [START_REF] Aaron | An introduction to matrix concentration inequalities[END_REF], and in (b) we use the fact that the inner maximum is always attained by 1/w k,i,j since by definition 1/w k-1,i,j is lower-bounded by 1/p k-1,i,j . Applying the inequality recursively from k = h to k = 1 removes all z s,i,j from the maximum and we are finally left with only w 0,i,j as we wanted, P max s=0,...,h z s,i,j p s,i,j ≤ a ≥ P max z 0,i,j p 0,i,j ; z 0,i,j w 0,i,j ≤ a ≥ P 1 w 0,i,j ≤ a ,

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

Step 4 (stochastic dominance on W h ). Now that we proved the stochastic dominance of 1/w 0,i,j , we plug this result in the definition of W h . For the sake of notation, we introduce the term p max h ,i,j to indicate the maximum over the first h step of copy i, j such that max s=0,...,h z s,i,j p s,i,j = 1 p max h ,i,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,i,j } i,j ) we can restate Eq. 13 as

W h = λ max (W h ) ≤ Λ({1/p max h,i,j } i,j ) def = λ max   1 q 2 q j=1 ν i=1 1 p max h,i,j 2 ψ i ψ T i ψ i ψ T i   .
In Step 4, we showed that 1/p max h,i,j is stochastically dominated by 1/w 0,i,j for every i and j. In order to bound Λ({1/p max h,i,j } i,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,i,j } i,j ) as

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

Λ({1/w 0,i,j } i,j ) def = max x: x =1 1 q 2 q j=1 ν i=1 1 w 0,i,j 2 ψ i 2 ψ T i x 2 ,
we now need to prove the stochastic dominance of Λ({1/w 0,i,j } i,j ) over Λ({1/p max h,i,j } i,j ). Using the definition of 1/p max h,i,j , w h,i,j , and the monotonicity of Λ we have 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,i,j and w k,i,j are all independent and (b) from the definition of 1/p max k-1,i,j and the fact that by definition 1/w k-1,i,j is lower-bounded by 1/p k-1,i,j . We can iterate this inequality to obtain the desired result Step 5 (concentration inequality). Since all w 0,i,j are (unconditionally) independent from each other, we can apply the following theorem.

P   Λ   1 p max h,i,j i,j   ≤ a   = P Λ max
µ max def = λ 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 {i,j} = 1 q 2 1 w 0,i,j ψ i ψ T i ψ i ψ T i 1 q 2 α 2 p 2 h,i ψ i ψ T i ψ i ψ T i 1 q 2 α 2 p 2 h,i ψ i ψ T i 2 I α 2 q 2 I,
where the first inequality follows from the definition of w 0,i,j in Eq. 14, the second from the PSD ordering, and the third from the definition of ψ i ψ T i .

Therefore, we can use L def = α 2 /q 2 for the purpose of Prop. 4. 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,i,j . Let us denote the c.d.f. of 1/w 2 0,i,j as F 1/w 2 0,i,j (a) = P Therefore, selecting h = 2, σ 2 = 6α/q and applying Prop. 4 we have While we do know that the z h,i,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,i,j , and in addition there could be dependencies between z h,i,j and z h,i ,j . Similarly to Lem. 8, we are going to find a stochastic variable to dominate z h,i,j . Denoting with u s,i,j ∼ U(0, 1) a uniform random variable, we will define w s,i,j as w s,i,j |F {s,i ,j } = w s,i,j |F s-2 def = I u s,i,j ≤ p h,i p s-1,i ∼ B p h,i p s-1,i

P W h ≥ σ 2 ≤ P   λ max   1 q 2 q j=1 ν i=1 1 w 2 0,i,j ψ i ψ T i ψ i ψ T i   ≥ (1 + 2) 2α q   ≤ ν • exp - 2α q q 2 α 2 (3 log(3) -2) ≤ n • exp -
for any i and j such that {s, 1, 1} ≤ {s, i , j } < {s, i, j}. Note that w s,i,j , unlike z s,i,j , does not have a recursive definition, and its only dependence on any other variable comes from p s-1,i . First, we peel off the last step

P   q j=1 ν i=1 z h,i,j ≥ g A   = E Ft-1|A   P   q j=1 ν i=1 I u h,i,j ≤ p h,i p t-1,i z h-1,i,j ≥ g F t-1 ∩ A     ≤ E Ft-1|A   P   q j=1 ν i=1 I u h,i,j ≤ p h,i p t-1,i z h-1,i,j ≥ g F t-1 ∩ A     = P   q j=1 ν i=1 w h,i,j z h-1,i,j ≥ g A   ,
where we used the fact that conditioned on A, I {h,l} is accurate w.r.t. K {h,l} , which guarantees that p h,i ≤ p h,i . Plugging this in the previous bound, 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. ν i=1 e p h,i (e θ -1) e θg ≤ inf θ>0 e q(e θ -1) ν i=1 p h,i e θg = inf θ>0 e (d eff (γ) {h,l} q(e θ -1)-θg) ≤ inf θ>0 e (d eff (γ) {h,l} q(e θ -1)-θg) ,

Figure 1 :

 1 Figure 1: Merge tree for Alg. 2 with an arbitrary partitioning and merging scheme.

  from the dataset D {3,1} . D {3,1} contains D 1 , D 2 , D 3 (descendent nodes are highlighted in red) and it has dimension (|D 1 |+|D 2 |+|D 3 |). Theorem 2 summarizes the guarantees for DISQUEAK.

I 1 , 2 , 3 , 4 , 5 I 1 , 2 , 3 + I 4 I 1 , 2 , 3 I 1 + I 2 , 3 I 1 I 1 D 1 I 2 , 3 I 2 + I 3 I 2 D 2 I 3 D 3 I 4 I 4 I 4 D 4 s = 4 s = 3 s = 2 s = 1 (a) arbitrary tree I 1 , 2 , 3 , 4 I 1 , 2 , 3 + I 4 I 1 , 2 , 3 I 1 , 2 + I 3 I 1 , 2 I 1 + I 2 I 1 D 1 I 2 D 2 I 3 I 3 D 3 I 4 I 4 I 4 D 4 s = 4 s = 3 s = 2 s = 1 (b) sequential tree I 1 , 2 , 3 , 4 I 1 , 2 + I 3 , 4 I 1 , 2 I 1 , 2 I 1 + I 2 I 1 D 1 I 2 D 2 I 3 , 4 I 3 + I 4 I 3 I 3 D 3 I 4 I 4 D 5 s = 4 s = 3 s = 2 s = 1 (

 123451234123123111232322334444432112341234123123121211223334444432112341234121212112234343334454321 c) minimum depth tree

Figure 2 :

 2 Figure 2: Merge trees for Algorithm 2.

Figure 3 :

 3 Figure 3: The dependence graph of the considered 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.

P( W h ≥ σ 2

 2

Proposition 4 (

 4 Tropp[START_REF] Aaron | An introduction to matrix concentration inequalities[END_REF], 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 def = k X k ,and the maximum eigenvalue of its expectation

1 - 1 a=0 1 -α 2 /p 2 h,i a=1 1 - 1 a=0( 1 - 1 =

 111221111 F 1/w0,i,j (a) da= F 1/w0,i,j (a) da + F 1/w0,i,j (a) da + ∞ a=α 2 /p 2 h,i 1 -F 1/w0,i,j (a) da = 2α/p h,i -1.Therefore,µ max = λ max (E [V]) = λ max {i,j} E X {i,j} = λ max

2q α • D. 4 z

 4 Space complexity boundDenote with A the event A = ∀h ∈ {1, . . . , h} : P h -P h 2 ≤ ε , and again ν = |D {h,l} |. Lettingq = |I {h,l} | = ν i=1 q h,i = q j=1 ν i=1 z h,i,j be the random number of points in I {h,l} , we reformulateP |I {h,l} | ≥ 3qd eff (γ) {h,l} ∩ ∀h ∈ {1, . . . , h} : P h -P h 2 ≤ ε ≤ ε = P |I {h,l} | ≥ 3qd eff (γ) {h,l} ∩ A = P h,i,j ≥ 3qd eff (γ) {h,l} ∩ A i,j ≥ 3qd eff (γ) {h,l} A   P (A) .

  i,j z h-1,i,j ≥ g ∩ A i,j z h-1,i,j ≥ g   .

III

  u s,i,j ≤ p h,i p s-1,i z s-1,i,j ≥ g F s-2 u s,i,j ≤ p h,i p s-1,i I u s-1,i,j ≤ p s-1,i p s-2,i z s-2,i,j ≥ g F s-2 u s-1,i,j ≤ p h,i p s-2,i z s-2,i,j ≥ g F s-2 1,i,j z s-2,i,j ≥ g Applying this repeatedly from s = h to s = 2 we have, the w 1,i,j are independent Bernoulli random variables, and we can bound their sum with a Hoeffding-like bound using Markov inequality, h,i e θ + (1p h,i ))e θg = infθ>0 q j=1 ν i=1 (1 + p h,i (e θ -1)) e θg ≤ inf θ>0 q j=1

Table 1 :

 1 Comparison of Nyström methods. λmax and λmin refer to largest and smallest eigenvalues of Kn.

  R t×tFor node {h, l} in the merge tree, given ν = |D {h,l} | Similarly defined K {h,l} ∈ R ν×ν , P {h,l} ∈ R ν×ν Similarly defined I {h,l} , S {h,l} ∈ R ν×ν , P {h,l} ∈ R ν×ν For layer h in the merge tree, block diagonal matrices K h , P h , P h with K {h,l} , P {h,l} , P {h,l} on the diagonalΨ = (K {h,l} + γI ν ) -1/2 K i ψ i ψ T i ∈ R ν×ν constructedusing weights p s,i from step s of Alg. 2 and ψ i from {h, l} Y s = P {h,l} -P {h,l} = P {h,l} -P {h,l} h ∈ R ν×ν sampling process based on K {h,l} , at final step h Y s = P {h,l} -P {h,l} s ∈ R ν×ν sampling process based on K {h,l} , at intermediate step s Y s ∈ R ν×ν sampling process based on K {h,l} , with freezing W h ∈ R ν×ν total variance of sampling process based on K {h,l} , with freezing

	{h,l} s	=	ν {h,l} i=1	qs,i q ps,

1/2 {h,l} ∈ R ν×ν , with ψ i its i-th column P

  D {h,l} and can be shortened as τ h,i = τ D {h,l} ,i . Using ψ i , we can also introduce the random matrix P

	{h,l} s	as
	P {h,l} s

,l} ,i ΨΨ T e ν {h,l} ,i = e T ν {h,l} ,i P {h,l} e ν {h,l} ,i = τ D {h,l} ,i , or, in other words, the norm ψ i ψ T i is equal to the RLS of the i-th sample w.r.t. to dataset D {h,l} . Note that since i is present only in node l on layer h, its RLS is uniquely defined w.r.t.

  1, 1/p h,i ≤ 1/w h,i,j . Therefore,

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

Now focus on an arbitrary intermediate step 1 ≤ k ≤ h, where we fix F k-1 . Since u k,i,j and w k,i,j are independent

  Now pick 1 ≤ k ≤ h, for a fixed F k-1 ,

				max s=0,...,h-1	z s,i,j p s,i,j	;	z h,i,j p h,i,j	i,j	≤ a
	≥ P Λ	max	max s=0,...,h-1	z s,i,j p s,i,j	;	z h,i,j w h,i,j	i,j	≤ a .
	1 p max k-1,i,j	is a constant and max

Proceedings of the 20 th International Conference on Artificial Intelligence and Statistics (AISTATS) 2017, Fort Lauderdale, Florida, USA. JMLR: W&CP volume 54. Copyright 2017 by the author(s).

Note that there is an important difference in whether the method passes through kernel matrix only once or through the dataset only once, in the former, the algorithm may still need access one data point up to n times, thus making it unsuitable for the streaming setting and less practical for distributed computation.

The technical report of Musco and Musco[START_REF] Musco | Provably useful kernel matrix approximation in linear time[END_REF] was developed independently from our work.

Notice that unlike zs,i,j, ws,i,j is no longer Fs-measurable but it is F s -measurable, where F {s,i,j} = u s ,i ,j , ∀{s , i , j } ≤ {s, i, j} ∪ {ws,i,j} = F {s,i,j} ∪ {ws,i,j} .

Acknowledgements The research presented was supported by French Ministry of Higher Education and Research, Nord-Pas-de-Calais Regional Council and French National Research Agency projects ExTra-Learn (n.ANR-14-CE24-0010-01) and BoB (n.ANR-16-CE23-0003)

where we use the fact that 1 + x ≤ e x , w 1,i,j ∼ B(p h,i ) and by Def. 2, ν i=1 p h,i = ν i=1 τ h,i = d eff (γ) {h,l} . The choice of θ minimizing the previous expression is obtained as d dθ e (qdeff(γ){h,l}(e θ -1)-θg) = e (qdeff(γ){h,l}(e θ -1)-θg) qd eff (γ) {h,l} e θg = 0,

and thus θ = log(g/(qd eff (γ) {h,l} )). Plugging this in the previous bound,

-1 e -qd eff (γ) {h,l} , and choosing g = 3qd eff (γ) {h,l} , we conclude our proof.

E Applications

Proof of Lemma 5.

From Lem. 6, and the fact that I n is ε-approximate we have that

and