Distributed Adaptive Sampling for Kernel Matrix Approximation
Daniele Calandriello Alessandro Lazaric Michal Valko

















Subroutinethplop
 



















 











 


 #1#1 #1 




















































































































assumptionAssumption
theoremTheorem
corollaryCorollary
definitionDefinition
lemmaLemma
propositionProposition

remark[1][Remark]

  #1



 





[

Distributed Adaptive Sampling for Kernel Matrix Approximation

 Daniele Calandriello Alessandro Lazaric Michal Valko 
SequeL team, INRIA Lille - Nord Europe ]


Most kernel-based methods, such as kernel or Gaussian process regression, kernel PCA, ICA, or k-means clustering,
do not scale to large datasets, because constructing and storing
the kernel matrix _n requires at least (n^2) time and space for
n samples.
Recent works <cit.> 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 
_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 , 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 γ of the dataset.
Moreover since all the RLS estimations are efficiently performed using only the small dictionary,
is the first RLS sampling algorithm that never constructs the whole matrix _n, runs in linear time (nγ^3) w.r.t. n,
and requires only a single pass over the dataset. We also propose a parallel and distributed version of that linearly scales across multiple machines, achieving similar accuracy in as little as (log(n)γ^3) time.







§ INTRODUCTION



One of the major limits of kernel ridge regression (KRR), kernel PCA <cit.>, and other kernel methods is that for n samples storing and
manipulating the kernel matrix _n requires (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 _n, thus reducing the space complexity to
(nm). These methods, often referred to as Nyström
approximations, mostly differ in the distribution used to sample the columns of
_n and the construction of low-rank approximations. Both of these
choices significantly affect the accuracy of the resulting
approximation <cit.>. <cit.> 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 (n) columns in
datasets with high coherence <cit.>,
i.e., a kernel matrix with weakly correlated
columns. On the other hand, <cit.> 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 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, <cit.>
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. <cit.> addressed this
issue by processing the dataset incrementally and updating estimates of
the ridge leverage scores, effective dimension, and approximations
on-the-fly. Although the space complexity of the resulting algorithm
() does not depend on the minimal eigenvalue anymore, it introduces
a dependency on the largest eigenvalue of _n, which in the worst
case can be as big as n, thus losing the advantage of the
method. In this paper we introduce an algorithm for SeQUEntial Approximation
of Kernel matrices (), a new algorithm that builds on ,
but uses unnormalized RLS.
This improvement, together with a new analysis, opens the way to 
major improvements over current leverage sampling methods (see Sect. <ref> for a comparison with existing methods)
closely matching the dictionary size achieved by exact RLS sampling.
First, unlike , 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, needs to actually observe only a fraction of the kernel matrix _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 requires a more sophisticated
analysis that take into consideration the complex interactions and dependencies between successive resampling steps.
The analysis of builds on a new martingale argument that could be of independent interest for similar online resampling schemes.
Moreover, our 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 
<cit.>, but with negative algorithmic side effects: these methods 
need to pass through the dataset multiple times. passes 
through the dataset only once[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.]
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 𝒜 for collections. We use upper-case
bold letters A for matrices and operators, lower-case bold letters a
for vectors,
lower-case letters a for scalars, with the exception of f,g, and h which
denote functions, and [n] := {1,…,n} for the set of
integers between 1 and n. 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∈^n× n, the identity matrix of dimension n and by
( a)∈^n× n the diagonal matrix with the vector
a∈^n on the diagonal. We use e_n,i∈^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) operator.

Kernel.
We consider a positive definite kernel function
: ×→ and we denote with
 its induced Reproducing Kernel Hilbert Space (RKHS),
and with φ: → its corresponding feature map.
Using φ, and without loss of generality,
for the rest of the paper we will replace  with a high dimensional
space ^D where D is large and potentially infinite.
With this notation, the kernel evaluated between to points can be expressed as
( x,  x') = ( x,·)( x', ·) = φ( x)φ( x') = φ( x)^φ( x').
Given a dataset of points ={ x_t}_t=1^n, we define the (empirical) kernel matrix
_t∈^t× t as the
application of the kernel function on all pairs of input values (i.e.,
[_t]_ij = k_i,j = ( x_i,  x_j) for any i,j∈ [t]),
with k_t,i = _t  e_t,i as its i-th column.
We also define the feature vectors _i = φ( x_i) ∈^D and after
introducing the matrix 
_t = [ _1,_2, …, _t ] ∈^D × t                                     we can rewrite the kernel matrix as _t = _t^_t. 

 
Kernel approximation by column sampling.
One of the most popular strategies to have low space complexity approximations
of the kernel _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 _t = {(i,w_i)}_i=1^t, 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 _t to t, while in practice, we only
store the non-zero elements. In particular, we denote by |_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
_t = (√(w_1)…√(w_t))∈^t× t such
that for any matrix A_t∈^t × t, A_t _t returns a
t× t matrix where the columns selected by _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 measured 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 _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 _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  (_t + γ I)^-1/2_t(_t + γ I)^-1/2.

We say that a column dictionary is accurate if the following condition is satisfied.

A dictionary _t = {(i,w_i)}_i=1^t and its associated selection matrix _t ∈^t × t are ε-accurate w.r.t. a kernel matrix _t = _t^1/2_t^1/2
if [the matrix norm we use is the operator (induced) norm]

    P_t -  P_t ≤ε,

where for a given γ > 0, the approximated projection matrix is defined as

    P_t  (_t + γ I_t)^-1/2_t^1/2_t _t^_t^1/2(_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. <ref>, ε-accurate dictionaries can be used to construct suitable kernel approximation in a wide range of problems.

Ridge leverage scores sampling. <cit.> showed that an ε-accurate dictionary can be obtained by sampling columns proportionally to their γ-ridge leverage scores (RLS) defined as follows.


Given a kernel matrix K_t∈^t× t, the γ-ridge leverage score (RLS) of column i∈ [t] is

    τ_t,i =  e_t,i^_t(_t + γ I_t)^-1 e_t,i,

Furthermore, the effective dimension γ_t of the kernel matrix _t is defined as

    γ_t = ∑_i=1^t τ_t,i(γ) = (_t(_t + γ I_t)^-1).



The RLS can be interpreted and derived in many ways, and they are well studied
<cit.> in the linear setting (e.g. _t =  x_t).
<cit.> used them
as a measure of incoherence to select important points,
but their deterministic algorithm provides guarantees only
when _t is exactly low-rank.
Here we notice that

    τ_t,i = e_t,i^ K_t(_t + γ I_t)^-1 e_t,i
    =  e_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, <cit.> state the following.

Let ε∈[0,1] and _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 = (1/ε^2_nlog(n/δ)), then w.p. at least 1-δ, the corresponding dictionary is ε-accurate.


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


 



§ SEQUENTIAL RLS SAMPLING





[t!]
       
    The algorithm
    


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  algorithm
(Alg. <ref>) 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 shares a similar structure with
 <cit.>, the sampling probabilities are
computed from different estimates of the RLS τ_t,i and no renormalization by
an estimate of _t is needed. Before giving the details of the
algorithm, we redefine a dictionary as a collection =
{(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/(qp_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.

receives as input a dataset ={ x_t}_t=1^n and processes it sequentially. Starting with an empty dictionary _0, at each time step t, 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.


    For any kernel matrix _t-1 at time t-1 and its extension
    _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,      γ_t≥γ_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 _t-1
to reflect the changes of the τ_t,i. We proceed in two phases.
During the phase, we directly add the new element
x_t to _t-1 and obtain a
temporary dictionary , 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., = _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 ,
in order to decrease its size.
Given as input , we use the following estimator to compute
the approximate RLS _t,i,
    
    _t,i   = (1-ε)_i^(_t S S^_t^ + γ I)^-1_i
        =1-εγ(k_i,i -  k_t,i^(^_t + γ I_t)^-1^ k_t,i),

where ε is the accuracy parameter, γ is the regularization
and  is the selection matrix associated to
. 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 =
_i^(_t I_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
(^≠ I_t),
similarly to what <cit.> did for the special case _i =  x_i,
we choose to use
_t,i≈^(_t^_t^ +
γ I)^-1_i, and then we use the kernel trick
to derive a form that we can actually compute, resulting in Eq. <ref>.
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 , the step draws a sample from
the binomial ℬ(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
_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 _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.


    Given an ε-approximate dictionary _t-1
    of matrix _t-1, construct  by adding element
    (t, 1, q) to it, and compute the selection matrix .
    Then for all i in  such that q_t-1,i≠ 0, the estimator in Eq. <ref> 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 _t-1 is
ε-accurate for _t-1, the projection matrix P_t
can be approximated by P_t constructed using the temporary dictionary
 and thus, the RLSs can be accurately estimated and used to
update _t-1 and obtain a new ε-accurate dictionary for
_t. Since τ_t,i is used to sample the new dictionary
_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. <ref>. 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. <ref>, which
guarantees that the estimated RLS are always smaller than the exact RLS. As a
result, 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.

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

    q =  39αlog(2n/δ)/ε^2

where  α = (1+ε)/(1-ε). Then, w.p. at least 1-δ,
 generates a sequence of random dictionaries {_t}_t=1^n that are ε-accurate (Eq. <ref>) w.r.t. any of the intermediate kernels _t, and the size of the dictionaries is bounded as max_t=1,…,n |_t| ≤ 3qγ_n.

As a consequence, on a successful run the overall complexity of is bounded as

    space complexity   =  (max_t=1,…,n |_t|)^2 ≤(3qγ_n)^2,
    time complexity   =(n γ_n^3q^3).



We show later that Thm. <ref> is special case of Thm. <ref>
and give a sketch of the proof with the statement of Thm. <ref>.
We postpone the discussion about this result and the comparison with previous
results to Sect. <ref> and focus now on the space
and time complexity. Note that while the dictionaries _t always
contain t elements for notational convenience, 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
_t-1 and the new point x_t (i.e., the elements in
) and thus it does not require recomputing the RLSs of points
x_s (s < t) that have been dropped before!
This is why computes an ε-accurate dictionary
with a single pass over the dataset. Furthermore, the estimator in
Eq. <ref> 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 , are directly set
to zero when computing k_t,i^. As a result, for any new
point x_t we need to evaluate ( x_s, x_t) only for the indices
s in .
Therefore, never performs more than n(3qγ_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. <ref>.
Therefore, the total runtime is (n(max_t=1,…,n |_t|)^3)
= (n γ_n^3q^3). In the next section, we introduce
, 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 allows us to parallelize and distribute the computation of the dictionary
_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 i guess this is some distributed lingo? on bottleneck data sources
such as databases or network connections.
DIstributed-(, Alg. <ref>) partitions  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 _i with i=1,…,k and k
dictionaries __i = {(j,p_0,i = 1, q_0,i = q) : j ∈_i}
are initialized simply by placing all samples in _i into 
with weight 1 and multiplicity q. Alternatively, if the datasets
_i are too large to fit in memory, we can run to generate the initial dictionaries.
The dictionaries __i are added to a dictionary collection
_1, and at each step h ∈ [k] Alg. <ref>
arbitrarily chooses two dictionaries _ and _'
from _h, and merges them.
first combines them into a single dictionary  (the equivalent of
the phase in ) and then is run on the merged dictionaries
 to create an updated dictionary _∪', which
 is placed back in the dictionary collection .
This sequence of merges can be represented using a binary merge tree,
as in Fig. <ref>.
Since 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 .
Unlike in , 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 <ref> and <ref> to analyze the
behavior of the RLSs and the quality of the estimator.

    Given two disjoint datasets , ',
    for every i ∈∪',
    τ_i,≥τ_i,∪'
    and
    
    2γ_∪'≥γ_ + γ_'≥γ_∪'.




need to reflect the fact that we don't need bound on RLS decrease anymore
While in we were merging an ε-accurate dictionary _t and a new
point, which is equivalent to a perfect, 0-accurate dictionary,
in both dictionaries used in a merge are only ε-accurate.
To balance this change, we introduce a new
estimator,
    
    _∪',i = 1-εγ(k_i,i-  k_i^(^ + (1+ε)γ I)^-1^ k_i),




where  is the selection matrix associated with the temporary dictionary = _∪_'.
Eq. <ref> has similar guarantees as Lem. <ref>,
with only a slightly larger α.



    Given two disjoint datasets , ',
    and two ε-approximate dictionaries
    _, _',
    let = _∪_'
    and  be the associated selection matrix.
    Let  be the kernel matrix computed on ∪',
    k_i its i-th column, and τ_∪',i
    the RLS of k_i.
    Then for all i in  such that q_i≠ 0, the estimator in Eq. <ref> is
    α-accurate, i.e. it satisfies τ_∪',i/α≤_∪',i≤τ_∪',i, with α = ( 1-ε )/( 1+3ε ).
                


Given these guarantees, the analysis of follows similar steps as . Given
ε-accurate dictionaries, we obtain α-accurate RLS estimates
_∪',i that can be used to resample all points in
 and generate a new dictionary _,'
that is ε-accurate.
To formalize a result equivalent to Thm. <ref>, 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 _{h,l} and the collection of
all dictionaries available at height h of the merge tree by _h = {_{h,l}}. We also use
_{h,l} to refer to the kernel matrix constructed from the
datasets _{h,l}, which contains all points present in the leaves reachable from node {h,l}.
For instance in
Fig. <ref>, node {3,1} is associated with _{3,1},
which is an ε-approximate dictionary of the kernel matrix _{3,1}
constructed from the dataset _{3,1}. _{3,1} contains _1,
_2, _3 (descendent nodes are highlighted in red) and it has dimension
(|_1|+|_2|+|_3|).
Theorem <ref> summarizes the guarantees for .

 
Let ε>0 be the
accuracy parameter, γ > 1 the regularization factor, and 0 < δ
< 1 the prob. of failure. Given an arbitrary dataset  and a merge tree structure of height k as
input together with parameters ε, γ, and δ, we run
with

    q =  39αlog(2n/δ)/ε^2·

where α = (1+3ε)/(1-ε). Then, w.p. at least 1-δ, generates a sequence of collections of dictionaries {_h}_h=1^k such that each dictionary _{h,l} in _h is ε-accurate (Eq. <ref>) w.r.t. to _{h,l}, and that at any node l of height h the size of the dictionary is
bounded as |_{h,l}| ≤ 3qγ_{h,l}.
The cumulative (across nodes) space and time requirementsof the algorithm depend on the exact shape of the merge tree.




Theorem <ref> gives approximation and space guarantees
for every node of the tree. In other words, it guarantees that
each intermediate dictionary processed by 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, provides exactly the same
guarantees of .
Analysing the complexity of is however more complex,
since the order and arguments of the 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 _i = { x_i}.
In the fully unbalanced tree, we always merge the current dictionary with a new dataset
(a single new point) and no operation can be carried
out in parallel. Unsurprisingly, the sequential algorithm induced by this
merge tree is strictly equivalent to .
Computing a solution in the fully unbalanced tree takes (nγ_^3q^3) time
with a total work that is also (nγ_^3q^3),
as reported in Thm. <ref>.
On the opposite end, the fully balanced tree needs to invert
a γ_{h,l} dimensional matrix at each layer of the tree for a total of log(n) layers. Bounding all γ_{h,l}
with γ_, gives a complexity for computing the final solution of (log(n)q^3γ_^3) time, with a huge
improvement on .  Surprisingly, the total work is only
twice (nq^3γ_^3), since
at each layer h we perform n/2^h inversions (on n/2^h machines),
and the sum across all layers is ∑_h=1^log(n) n/2^h≤ 2n.
Therefore, we can compute a solution in a much shorter time than ,
with a comparable amount of work, but at the expense of requiring much more memory
across multiple machines, since at layer h, the sum ∑_l=1^|_h|γ_{h,l}
can be much larger than γ_. Nonetheless,
this is partly alleviated by the fact that each node {h,l}
locally requires only γ_{h,l}^2 ≤γ_^2
memory.

Proof sketch:
Although 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. <ref>
generates the dictionary _{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 <cit.>, that 1) uses
Freedman's inequality 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 _n
from an ε-accurate dictionary.

    Given an ε-accurate dictionary _t
    of matrix _t, and the selection
    matrix S_t, the regularized
Nyström approximation of _t is defined as
    
    _n = _n _n(_n^_n _n + γ I_m)^-1_n^_n,

    and satisfies
    
    0 ≼_t - _t ≼γ/1-ε_t(_t + γ I)^-1≼γ/1-ε I.





This is not the only choice of an approximation from dictionary _n. For instance,
<cit.> show similar result for an unregularized approximation
and <cit.> for a smaller _n, construct the estimator only for the points in
_n.
Let
C = _n _n ∈^n × m, W = (_n^_n _n + γ I_m) ∈^m × m,
with m = |_n|,
and 
using the Woodbury formula define the regression weights as 
    w_n =    (_n + μ I_n)^-1 y_n = ( C W^-1 C^ + μ I_n)^-1 y_n
    
            =   1/μ( y_n -  C(  C^ C + μ W)^-1 C^ y_n).




Computing (  C^ C + μ W)^-1, inverting it, and the other matrix-matrix multiplication take (nm^2 + m^3) time,
and require to store at most an n × m matrix. Therefore the final
complexity of computing w_n is reduced from (n^3) to (nm^2 + m^3) time,
and from (n^2) to (nm) space.
We now provide guarantees
for the empirical risk of w_n in a fixed design setting.

    For an arbitrary dataset , let  be the kernel matrix
    constructed on .
    Run or with regularization
    parameter γ. Then, the solution
    w computed using the regularized approximation
     satisfies
    
    ℛ_( w) ≤   (1 + γ/μ1/1-ε)^2ℛ_( w),

    where μ is the regularization of kernel ridge regression
    and ℛ_( w) is the empirical risk on .

Using tools from <cit.>,
and under some mild assumption on the kernel function (·,·)
and the dataset ,
similar results can be derived for the random design setting.


Other applications. 
The projection P_n naturally appears in some form across nearly
all kernel-based methods. Therefore, in addition to KRR in the fixed <cit.>
and random <cit.> design setting, any kernel matrix approximation that provides ε-accuracy guarantees on P_n
can be used to provide guarantees for a variety of other kernel problems.
As an example, <cit.> show this is the case for
kernel PCA <cit.>, kernel CCA with regularization,
and kernel K-means clustering. Similarly, inducing points methods
for Gaussian Processes <cit.> can benefit from fast and provably accurate
dictionary construction.

 





§ DISCUSSION




Tab. <ref> compares several kernel approximation methods w.r.t. the time complexity
required to compute an ε-accurate dictionary, as well as the size of the final dictionary
|_n|. Note that for all methods the final complexity to construct
an approximate _n from _n (e.g. using Eq. <ref>) scales as (n|_n|^2 + |_n|^3) time
and (n|_n| + |_n|^2) space.
For all methods, we omit (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 <cit.> scales with the maximal degree
of freedom d_max. Since
d_max = n max_i τ_n,i≥∑_iτ_n,i = _n,
uniform sampling is often outperformed by RLS sampling.
Moreover, uniform sampling needs to know d_max in advance to guarantee
ε-accuracy, and this quantity is expensive to compute <cit.>.
While <cit.> 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 (nγ/(λ_min - nγ))
columns, which quickly grows above n^2 when λ_min becomes small.
Finally, <cit.> require that the maximum dictionary
size is fixed in advance, which implies some information about  the effective
dimensions _n, and requires estimating both _t,i and
γ_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 n^3 in the worst case. 
Therefore, we can see from the table that achieves the
same space complexity (up to constant factors) as knowing the RLS in advance and hence outperforms previous methods.
Moreover, when parallelized across k machines, can reduce this
linear runtime by a factor of k (linear scaling) with little communication cost.


A recent method by <cit.> achieves comparable space and time guarantees as .[The technical report of <cit.> was developed independently from our work.] While they rely on a similar estimator, the two approaches are very different. Their method is batch in nature, estimating leverage scores using repeated independent sampling from the whole dataset, and it requires multiple passes on the data. On the other hand, is intrinsically sequential and it only requires one single pass on , 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 <cit.> 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 . Furthermore, <cit.> requires centralized coordination after each of the sampling passes, and cannot leverage distributed architectures to match 's sub-linear runtime.


Future developments
Both and 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 _n and
indefinitely process new data beyond the initial dataset .
It is also interesting to see whether could be used in conjuntion with existing meta-algorithms (e.g., <cit.> 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.



 



  
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) 

plainnat
18


urlstyle

[Alaoui and Mahoney(2015)]alaoui2014fast
Ahmed El Alaoui and Michael W. Mahoney.
Fast randomized kernel methods with statistical guarantees.
In Neural Information Processing Systems, 2015.

[Bach(2013)]bach2013sharp
Francis Bach.
Sharp analysis of low-rank kernel matrix approximations.
In Conference on Learning Theory, 2013.

[Calandriello et al.(2016)Calandriello, Lazaric, and
  Valko]calandriello2016analysis
Daniele Calandriello, Alessandro Lazaric, and Michal Valko.
Analysis of Nyström method with sequential ridge leverage
  scores.
In Uncertainty in Artificial Intelligence, 2016.

[Cohen et al.(2016)Cohen, Musco, and Pachocki]cohen2016online
Michael B. Cohen, Cameron Musco, and Jakub Pachocki.
Online row sampling.
International Workshop on Approximation, Randomization, and
  Combinatorial Optimization, 2016.

[Cohen et al.(2017)Cohen, Musco, and Musco]cohen2017input
Michael B. Cohen, Cameron Musco, and Christopher Musco.
Input sparsity time low-rank approximation via ridge leverage score
  sampling.
In Symposium on Discrete Algorithms, 2017.

[Gittens and Mahoney(2013)]gittens2013revisiting
Alex Gittens and Michael W Mahoney.
Revisiting the Nyström method for improved large-scale machine
  learning.
In International Conference on Machine Learning, 2013.

[Kumar et al.(2012)Kumar, Mohri, and Talwalkar]kumar_sampling_2012
Sanjiv Kumar, Mehryar Mohri, and Ameet Talwalkar.
Sampling methods for the Nyström method.
Journal of Machine Learning Research, 130
  (Apr):0 981–1006, 2012.

[Levy(2015)]levy2015stochastic
Haim Levy.
Stochastic dominance: Investment decision making under
  uncertainty.
Springer, 2015.

[Musco and Musco(2016)]musco2016provably
Cameron Musco and Christopher Musco.
Provably useful kernel matrix approximation in linear time.
Technical report, 2016.
URL <http://arxiv.org/abs/1605.07583>.

[Pachocki(2016)]pachocki2016analysis
Jakub Pachocki.
Analysis of resparsification.
Technical report, 2016.
URL <http://arxiv.org/abs/1605.08194>.

[Patel et al.(2015)Patel, Goldstein, Dyer, Mirhoseini, and
  Baraniuk]patel_oasis:_2015
Raajen Patel, Thomas A. Goldstein, Eva L. Dyer, Azalia Mirhoseini, and
  Richard G. Baraniuk.
oASIS: Adaptive column sampling for kernel matrix approximation.
Technical report, 2015.
URL <http://arxiv.org/abs/1505.05208>.

[Quiñonero-Candela and
  Rasmussen(2005)]quinonero-candela_unifying_2005
Joaquin Quiñonero-Candela and Carl Edward Rasmussen.
A unifying view of sparse approximate Gaussian process regression.
Journal of Machine Learning Research, 60
  (Dec):0 1939–1959, 2005.
URL <http://www.jmlr.org/papers/v6/quinonero-candela05a.html>.

[Rasmussen and Williams(2006)]rasmussen_gaussian_2006
Carl Edward Rasmussen and Christopher K. I. Williams.
Gaussian processes for machine learning.
Adaptive computation and machine learning. MIT Press, Cambridge,
  Mass, 2006.
ISBN 978-0-262-18253-9.
OCLC: ocm61285753.

[Rudi et al.(2015)Rudi, Camoriano, and Rosasco]rudi2015less
Alessandro Rudi, Raffaello Camoriano, and Lorenzo Rosasco.
Less is more: Nyström computational regularization.
In Neural Information Processing Systems, 2015.

[Schölkopf et al.(1999)Schölkopf, Smola, and
  Müller]scholkopf1999kernel
Bernhard Schölkopf, Alexander J. Smola, and Klaus-Robert Müller.
Kernel principal component analysis.
In Advances in kernel methods, pages 327–352. MIT Press
  Cambridge, MA, USA, 1999.

[Tropp(2011)]tropp2011freedman
Joel Aaron Tropp.
Freedman's inequality for matrix martingales.
Electronic Communications in Probability, 16:0
  262–270, 2011.

[Tropp(2015)]tropp2015an-introduction
Joel Aaron Tropp.
An introduction to matrix concentration inequalities.
Foundations and Trends in Machine Learning, 80
  (1-2):0 1–230, 2015.

[Woodruff et al.(2014)]woodruff2014sketching
David P Woodruff et al.
Sketching as a tool for numerical linear algebra.
Foundations and Trends in Theoretical Computer Science,
  100 (1–2):0 1–157, 2014.


 








§ 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.

For any symmetric matrix A ∈^n × m and any γ > 0

    A( A^ A + γ I_m)^-1 A^ =  A A^( A A^ + γ I_n)^-1.

For any symmetric matrix A ∈^n × n and diagonal matrix B ∈^n × n
such that B has n-s zero entries, and s non-zero entries,
define C ∈^n × s as the matrix obtained by removing all
zero columns in B.
Then

    A B( B A B + γ I_m)^-1 B A
    = A C( C^ A C + γ I_s)^-1 C^ A.

For any appropriately shaped matrix A, B, C, with A and B invertible,
the Woodbury matrix identity states

    ( A +  C B C^)^-1 =  A^-1 -  A^-1 C( C^ A^-1 C +  B^-1)^-1 C^ A^-1.




Kernel matrix.
Since _t is real symmetric matrix, we can eigendecompose it as

    _t =  U_t Λ_t U_t^,

where U_t ∈^t × t is the eigenvector matrix and Λ_t ∈^t × t is the diagonal eigenvalue matrix, with all non-negative elements since _t is PSD. Considering the feature mapping from  to the RKHS, we can write the kernel matrix as

    _t = _t^_t,

where _t ∈^D × t is the feature matrix, whose SVD decomposition is

    _t =  V_t Σ_t U_t^,

where V_t ∈^D × D contains the left singular vectors, Σ_t ∈^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 _t. Furthermore, we have 

    Λ_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.

The following identity holds.

    P_t -  P_t2   =(_t + γ I_t)^-1/2_t^1/2( I_t -  S_t S_t^)_t^1/2(_t + γ I_t)^-1/22
       = (_t_t^ + γ I_D)^-1/2_t( I_t -  S_t S_t^)_t^(_t_t^ + γ I_D)^-1/22.



Using the SVD decomposition,

    (_t + γ I_t)^-1/2_t^1/2( I_t -  S_t S_t^)_t^1/2(_t + γ I_t)^-1/22
       =( U_t Λ_t U^ + γ U_t I_t U_t^)^-1/2 U_t Λ_t^1/2 U^( I_t -  S_t S_t^) U_t Λ_t^1/2 U^( U_t Λ_t U^ + γ U_t I_t U_t^)^-1/22
       = U_t( Λ_t + γ I_t)^-1/2 Λ_t^1/2 U^( I_t -  S_t S_t^) U_t Λ_t^1/2( Λ_t + γ I_t)^-1/2 U^2
       = U_t( Λ_t + γ I_t)^-1/2 Λ_t^1/2 U^( I_t -  S_t S_t^) U_t Λ_t^1/2( Λ_t + γ I_t)^-1/2 U^2
       =( Λ_t + γ I_t)^-1/2 Λ_t^1/2 U^( I_t -  S_t S_t^) U_t Λ_t^1/2( Λ_t + γ I_t)^-1/22
       =( Σ_t Σ_t^ + γ I_D)^-1/2 Σ_t U^( I_t -  S_t S_t^) U_t Σ_t^( Σ_t Σ_t^ + γ I_D)^-1/22
       = V_t( Σ_t Σ_t^ + γ I_D)^-1/2 Σ_t U^( I_t -  S_t S_t^) U_t Σ_t^( Σ_t Σ_t^ + γ I_D)^-1/2 V_t^2
       = V_t( Σ_t U_t U_t^ Σ_t^ + γ I_D)^-1/2 V_t^ V_t Σ_t U^( I_t -  S_t S_t^) U_t Σ_t^ V^ V( Σ_t U_t U_t^ Σ_t^ + γ I_D)^-1/2 V_t^2
       = (_t_t^ + γ I_D)^-1/2_t( I_t -  S_t S_t^)_t^(_t_t^ + γ I_D)^-1/22.







 


§ RIDGE LEVERAGE SCORES AND EFFECTIVE DIMENSION (PROOF OF LEMMA <REF> AND <REF>)


From <cit.> we know that τ_t,i≤τ_t-1,i and
γ_t ≥γ_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   = _i^(_t_t^ + γ I)^-1_i
    = _i^(_t-1_t-1^ + _t_t^ + γ I)^-1_i
       = _i^(_t-1_t-1^  + γ I)^-1_i
    - _i^(_t-1_t-1^  + γ I)^-1_t_t^(_t-1_t-1^  + γ I)^-1_i/1+_t^(_t-1_t-1^  + γ I)^-1_t
       = τ_t-1,i
    - _i^(_t-1_t-1^  + γ I)^-1_t_t^(_t-1_t-1^  + γ I)^-1_i/1+_t^(_t-1_t-1^  + γ I)^-1_t·

Let

    x = (_t-1_t-1^  + γ I)^-1/2_i
      and   y = (_t-1_t-1^  + γ I)^-1/2_t.

Then τ_t,i/τ_t-1,i is equal to

    τ_t,i/τ_t-1,i = 1
    - (_t^(_t-1_t-1^  + γ I)^-1_i)^2/(1+_t^(_t-1_t-1^  + γ I)^-1_t)_i^(_t-1_t-1^  + γ I)^-1_i   = 1 -  ( y^ x)^2/(1 +  y^ y) x^ x·

Defining the cosine between y and x
as cos( y, x) =  y^ x/( x y),
we have that

    1 -  ( y^ x)^2/(1 +  y^ y) x^ x
    = 1 -   y^ y x^ xcos( y, x)^2/(1 +  y^ y) x^ x
    = 1 -  y^2/ 1 +  y^2cos( 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^ x)^2/(1 +  y^ y) x^ x
    = 1 -  y^2/ 1 +  y^2cos( y, x)^2
    ≥ 1 -  x^2/ 1 +  x^2
    = 1/ 1 +  x^2
    = 1/1+τ_t,i

which concludes the proof of Lem. <ref>.

For Lem. <ref>, the first point τ_i,≥τ_i,∪' can be easily proven by choosing  to construct a kernel matrix _, and then invoke
Lem. <ref> as we add one sample at a time
from '.
Also as easily for γ_ + γ_'
we have

    γ_ + γ_'≤ 2max{γ_; γ_'}≤ 2max{γ_∪'; γ_∪'}
        = 2γ_∪'.

Finally, we prove the other side of the inequality for γ_ + γ_.
Let _, _' be the matrices
constructed using the feature vectors of the samples in  and
' respectively. Then,

    γ_ + γ_'   = ∑_i ∈τ_, i + ∑_i ∈'τ_',i
        = ∑_i ∈ ϕ_i^(__^ + γ I_D)^-1 ϕ_i + ∑_i ∈' ϕ_i^(_'_'^ + γ I_D)^-1 ϕ_i
       ≥∑_i ∈ ϕ_i^(_∪'_∪'^ + γ I_D)^-1 ϕ_i + ∑_i ∈' ϕ_i^(_∪'_∪'^ + γ I_D)^-1 ϕ_i
       = ∑_i ∈∪' ϕ_i^(_∪'_∪'^ + γ I_D)^-1 ϕ_i
        =∑_i ∈∪'τ_∪',i
        = γ_∪'.







§ RIDGE LEVERAGE SCORES ESTIMATION (PROOF OF LEMMA <REF> AND <REF>)



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

    τ_t,i   =  e_t,i^_t(_t + γ I_t)^-1 e_t,i
    =  e_t,i^_t^_t(_t^_t + γ I_t)^-1 e_t,i
       =  e_t,i^_t^(_t_t^ + γ I_D)^-1_t e_t,i
    =  ϕ_i^(_t_t^ + γ I_D)^-1 ϕ_i.

This formulation, combined with Def. <ref>, suggests
ϕ_i^(_t S_t S_t^_t^ + γ I_D)^-1 ϕ_i
as an estimator for τ_t,i. However, at step t, we only have access to an ε-accurate dictionary
_t-1 and not _t. Therefore, we augment it with (t,1,q) to construct
_t and the corresponding _t, which will have
[_t]_t,t = 1.
We will now show how to implement this estimator efficiently.
From Prop. <ref>, we apply Woodbury matrix identity
with A = γ I, B =  I and C = _t S_t,

    _t,i   = (1-ε) ϕ_i^(_t S_t S_t^_t^ + γ I_D)^-1 ϕ_i
       = (1-ε) ϕ_i^(_t S_t I_t S_t^_t^ + γ I_D)^-1 ϕ_i
    
    (Prop. <ref>)    = (1-ε) ϕ_i^(1/γ I_D - 1/γ^2_t S_t(1/γ S_t^_t^_t S_t +  I_t)^-1 S_t^_t^) ϕ_i
       = (1-ε)/γ ϕ_i^( I_D - _t S_t( S_t^_t^_t S_t + γ I_t)^-1 S_t^_t^) ϕ_i
       = (1-ε)/γ( ϕ_i^ ϕ_i -  ϕ_i^_t S_t( S_t^_t^_t S_t + γ I_t)^-1 S_t^_t^ ϕ_i)
       = (1-ε)/γ(k_i,i -  k_t,i S_t( S_t^_t S_t + γ I_t)^-1 S_t^ k_t,i),

which is the estimator defined in Eq. <ref>.

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 _k.
Given disjoint datasets {_i}_i=1^k with associated
feature matrices _i. From each dataset, construct
an ε-accurate dictionary _i, with its associated
selection matrix _i.

To estimate the RLS τ_i of point i w.r.t. the whole
dataset = ∪_j=1^k_j, and corresponding feature matrix , we set the estimator to be

    _i   = (1-ε) ϕ_i^(∑_j=1^k_j S_j S_j^_j^ + (1+(k-1)ε)γ I_D)^-1 ϕ_i.




Part 1: accuracy of the RLS estimator τ_i. 
Since each of the dictionaries _i used to generate S_i is ε-accurate,
we can use the equivalence
from Lem. <ref>, 

    P_i -  P_i2   = (_i_i^ + γ I_D)^-1/2(_i_i^ - _i S_i S_i^_i^)(_i_i^ + γ I_D)^-1/22≤ε,

which implies that 

    (1-)_i_i^ - εγ I_D ≼_i S_i S_i^_i^≼ (1+)_i_i^ + εγ I_D.

Therefore, we have

    _i   = (1-ε) ϕ_i^(∑_j=1^k_j S_j S_j^_j^ + (1+(k-1)ε)γ I_D)^-1 ϕ_i
       ≤ (1-ε) ϕ_i^(∑_j=1^k((1-ε)_j_j^-εγ I_D) + (1+(k-1)ε)γ I_D)^-1 ϕ_i
       ≤ (1-ε) ϕ_i^((1-ε)^ - kεγ I_D + (1+(k-1)ε)γ I_D)^-1 ϕ_i
       = (1-ε) ϕ_i^((1-ε)(^ + γ I_D))^-1 ϕ_i
    = (1-ε)/(1-ε) ϕ_i^(^ + γ I_D)^-1 ϕ_i = τ_i,

and

    _i   = (1-ε) ϕ_i^(∑_j=1^k_j S_j S_j^_j^ + (1+(k-1)ε)γ I_D)^-1 ϕ_i
       ≥ (1-ε) ϕ_i^(∑_j=1^k((1+ε)_j_j^ + εγ I_D) + (1+(k-1)ε)γ I_D)^-1 ϕ_i
       = (1-ε) ϕ_i^((1+ε)^ + kεγ I_D + (1+(k-1)ε)γ I_D)^-1 ϕ_i
       = (1-ε) ϕ_i^((1+ε)^ +(1 + (2k -1)ε) γ I_D))^-1 ϕ_i
       ≥ (1-ε) ϕ_i^((1+(2k -1)ε)^ +(1 + (2k -1)ε) γ I_D))^-1 ϕ_i
       = (1-ε)/(1+(2k-1)ε) ϕ_i^(^ - γ 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. <ref>,
and with k=2 to prove the accuracy claim in Lem. <ref>.


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. <ref>, we know that τ_t-1≥τ_t.
Given α-accurate _t and _t-1 we have the upper bound

    min{_t, _t-1}≤min{τ_t, τ_t-1}
        = τ_t,

and the lower bound,

    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.
 





§ PROOF OF THM. <REF> AND THM. <REF>



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

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. <ref>.

We partition  into k disjoint sub-datasets _i
of size n_i, such that = ∪_i=1^k_i.
For each dataset _i, we construct an initial dictionary
_{1,i} = {(j,p_0,i = 1, q_0,i = q) : j ∈_i}
by inserting all points from _i into __i with
weight p_0,i = 1 and number of copies q_0,i = q.
It is easy to see that _{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
_i, we can run Alg. <ref>
on _i to generate _{1,i}, and the following proof
will remain valid.
Regardless of construction, the initial dictionaries _{1,i} are included into the dictionary pool _1.

At iteration h, the inner loop of Alg. <ref>
arbitrarily chooses  two dictionaries from
_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 <ref> 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
_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 _h (number of dictionaries present at layer h)
is |_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≤ |_h| is the index of the node in the
layer. For example, in Figure <ref>, the node containing
_1,2,3 is indexed as {3,1}, and the highest
node containing _4 is indexed as {3,2}.

We also define the dataset
_{h,l} as the union of all sub-datasets _l' that are
 reachable from node {h,l} as leaves.
For example, in Fig. <ref>, dictionary _1,2,3
in node {3,1} is constructed starting from all points in _{3,1} = _1 ∪_2 ∪_3,
where we highlight in red the descendant tree.
We now define ^h as the block diagonal kernel matrix
where each diagonal block _{h,l} is constructed on
_{h,l}. Again, from Fig. <ref>, ^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 _3,1 constructed
on _{3,1} = _1 ∪_2 ∪_3, and a second n_4 × n_4 block
_3,2 constructed on _{3,2} = _4. Similarly, we can 
adapt Def. <ref> to define P^h as a block diagonal projection matrix,
where block P_{h,l} is defined using _{h,l},
and block diagonal P^h,
where block P_{h,l} is defined using _{h,l}
and _{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 _{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

    (∃  h∈{1,…,k}:  P^h -  P^h_2 ≥ε ∪ max_l=1,…,|_h||_{h,l}| ≥ 3qγ_{h,l})
       = (∃  h∈{1,…,k} : (max_l=1,…,|_h| P_{h,l} -  P_{h,l}_2) ≥ε_A_h ∪ (max_l=1,…,|_h||_{h,l}| ≥ 3qγ_{h,l})_B_h) ≤δ,

where event A_h refers to the case when some dictionary _{h,l} at an intermediate layer h fails to accurately approximate _{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 _{h,l} at a certain layer h). We can conveniently decompose the previous joint (negative) event into two separate conditions as

    (⋃_h = 1^k A_h ∪ B_h) =({⋃_h = 1^k A_h}∪{⋃_h = 1^k B_h})
    =({⋃_h = 1^k A_h}) +({⋃_h = 1^k B_h}∩{⋃_h = 1^k A_h}^∁) 
       =({⋃_h = 1^k A_h}) +({⋃_h = 1^k B_h}∩{⋂_h = 1^k A_h^∁})
    =({⋃_h = 1^k A_h}) +(⋃_h = 1^k{B_h ∩{⋂_h' = 1^k A_h'^∁}}).

Applying this reformulation and a union bound we obtain

    (∃  h∈{1,…,k}:  P^h -  P^h_2 ≥ε ∪ max_l=1,…,|_h||_{h,l}| ≥ 3qγ_{h,l})
       = (∃  h∈{1,…,k} : (max_l=1,…,|_h| P_{h,l} -  P_{h,l}_2) ≥ε)
           + (∃  h∈{1,…,k} : max_l=1,…,|_h||_{h,l}| ≥ 3qγ_{h,l}∩{∀  h' ∈{1, …, h} :   P^h' -  P^h'_2 ≤ε})
       ≤∑_h = 1^k∑_l=1^|_h|( P_{h,l} -  P_{h,l}_2 ≥ε)
           +∑_h = 1^k∑_l=1^|_h|(|_{h,l}| ≥ 3qγ_{h,l}∩{∀  h' ∈{1, …, h} :   P^h' -  P^h'_2 ≤ε}) ≤δ.

As discussed in Sect. <ref>, 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. <ref> 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 |_h| = k - h + 1,
we have

    ∑_h = 1^k∑_l=1^|_h|δ/2k^2 = ∑_h = 1^k(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. <ref>
is that we can now analyze the processes that generated each P_{h,l} -  P_{h,l}
(and each dictionary _{h,l}) separately.
Focusing on a single node {h,l} restricts our
problem on a well defined dataset _{h,l},
where we can analyze the evolution of _{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. fixed vs random tree This prevents us from using concentration
inequalities for i.i.d. processes that are at the basis of the analysis of
uniform sampling <cit.> and the method proposed
by <cit.>. 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 _{h',l'}. This requires a careful definition
of the martingale process to still use matrix concentration inequalities for
fixed-size matrices (see Sect. <ref>). 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. <ref> 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. <ref> (step 2)), which
finally allows us to use an i.i.d. matrix concentration inequality to bound
the total variance (Sect. <ref>-(step 5)). This finally leads
to the bound on the accuracy. The bound on the space complexity (Sect. <ref>)follows similar
(but simpler) steps.






 §.§ Bounding the projection error 𝐏_{h,l} - 𝐏_{h,l}



The sequential process. Thanks to the union bound in Eq. <ref>,
instead of having to consider the whole merge tree followed by
Alg. <ref>, we can focus on each individual node
{h,l} and study the sequential process that generated its dictionary _{h,l}.
We will now map more clearly the actions
taken by Alg. <ref> 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 _{h,l} and the random dictionary _{h,l},
where the randomness influences both which points are included in _{h,l},
and the weight with which they are added.
fixed vs random tree

Note that since the merge tree is decided in advance, the dataset _{h,l} is
not a random object, and is fixed for the whole process.
Consider now a point
i ∈_{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 _{s,l'}, it means that either (1) Alg. <ref>
used _{s,l'} to compute p_s,i, and p_s,i to compute q_s,i,
or (2) at layer h Alg. <ref> 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} = |_{h,l}| so that we can use index i ∈ [_{h,l}] to index all points in
_{h,l}. Given the symmetric matrix = (_{h,l} + γ I)^-1/2_{h,l}^1/2 with its i-th column _i = (_{h,l} + γ I)^-1/2_{h,l}^1/2 e__{h,l},i, we can rewrite the projection matrix as that
P_{h,l} = ^ = ∑_i=1^_{h,l}_i_i^.
Note that

    _i_i^
    = _i^_i
    =  e__{h,l},i^^ e__{h,l},i
    =  e__{h,l},i^^ e__{h,l},i
    =  e__{h,l},i^ P_{h,l} e__{h,l},i = τ__{h,l},i,

or, in other words, the norm _i_i^ is equal to the RLS of
the i-th sample w.r.t. to dataset _{h,l}. Note that since i is present only in node l on layer h, its
RLS is uniquely defined w.r.t. _{h,l} and can be shortened as
τ_h,i = τ__{h,l},i.
Using _i, we can also introduce the random
matrix P^{h,l}_s as

    P^{h,l}_s = ∑_i=1^_{h,l}q_s,i/qp_s,i_i_i^
    = ∑_i=1^_{h,l}∑_j=1^qz_s,i,j/qp_s,i_i_i^.

where z_s,i,j are {0,1} r.v. such that q_s,i = ∑_j=1^qz_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. <ref>.
But, for
a general s ≠ h P^{h,l}_s does not have a direct interpretation in the context
of Alg. <ref>. It combines the vectors _i,
which are defined using _{h,l} at layer h, with the weights
p_s,i computed by Alg. <ref> 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. <ref>.

fixed vs random merge tree
Taking into account that we are now considering a specific node {h,l},
we can drop the index from the
dataset _{h,l} =, RLS τ__{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}_s=1^h defined as

    Y_s =  P_{h,l} -  P^{h,l}_s = 1/q∑_i=1^∑_j=1^q(1 - z_s,i,j/p_s,i)_i_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.

 


 §.§ Bounding 𝐘_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. <ref> 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}_s=1^h:
rewrite item (1)


  (1) On the one hand, the overall algorithm may fail in generating an accurate
dictionary at some intermediate iteration, and yet return an accurate dictionary at
the end. On the other hand, whenever one of the intermediate Y_s is larger than ε (inaccurate intermediate dictionary) we lose our guarantees for p_s,i for the whole process. This is because an inaccurate
p_s,i that underestimates too much p_s.i will influence all successive p_s',i through the minimum p_s',i = min{_s',ip_s'-1,i}.
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 _{h,l} is an ε-accurate dictionary of
_{h,l}, knowing Y_s =  P_{h,l} -  P^{h,l}_s≤ε
for s < h does not guarantee that all descendant _{s,l'} are ϵ-accurate their _{s,l'}.
Nonetheless, we will show that Y_s≤ε
is enough to guarantee that the intermediate estimate p_s,i
computed in Alg. <ref>, 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{ Y_s-1≤ε}
        + Y_s-1{ 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,

    (  Y_h≥ε)
    ≤(  Y_h≥ε),

and if we can bound (  Y_h≥ε)
we will have a bound for the failure probability of Alg. <ref>,
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 = { u_s,i,j≤p_s,i,j/p_s-1,i,j}z_s-1,i,j,

where u_s,i,j∼𝒰(0,1) is a [0,1] uniform random variable and p_s,i,j is defined as

    p_s,i,j =p_s,i{ Y_s-1≤ε∩ z_s-1,i,j = 1}
        +p_s-1,i,j{ 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 (  Y_h≥ε)
≤(  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,…,^2q], 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_i^,

which allows us to write the cumulative matrix as
Y_{s,i,j} = ∑_r=1^{s,i,j}  X_{s,i,j}
where the checkpoints {s,,q} correspond to Y_s,

    Y_{s,,q} =  Y_s = 
        1/q∑_i=1^∑_j=1^q(1-z_s,i,j/p_s,i,j)_i_i^.

Let _s be the filtration containing all the
realizations of the uniform random variables u_s,i,j up to the step s, that is _s = { u_s',i',j', ∀{s',i',j'}≤ s}. Again, we notice that _s defines the state of
the algorithm after completing iteration s because, unless a “freezing” happened,
Alg. <ref> 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 _{s,l'} (for some l' unique at layer s), they are fully determined by _s-1.
Furthermore, since _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 _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| _{s,i',j'} = z_s,i,j| _{s-1,,q}∼ℬ( p_s,i,j/p_s-1,i,j)

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| _{s-1,,q}⊥z_s,i',j'| _{s-1,,q}.

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.

    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/α.


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 _{h,l}, and simply
refer to it as . Similarly, we will refer with
_i to _{s,l'} (as in, the dataset used to
compute _s,i), and with _i to the samples
in  not contained in _i (complement of _i).
here we use a square S, have to recheck if it is introduced in main text or appendix
Define A as the || × |_i| matrix that contains
the columns of _s related to points in _i,
and similarly define B as the || × |_i| matrix that contains
the columns of _s related to points in _i,
where _s can be reconstructed by interleaving columns of A
and B.
From its definition in Eq.<ref>,
we know that _s,i is computed by Alg. <ref> as

    _s,i = (1-ε) ϕ_i^(_ A A^_^ + (1+ε)γ I_D)^-1 ϕ_i,

using only the points in A that are available at node {s,l'}.
From Lem.<ref>
we know that 

    Y_s
    = P_{h,l} -  P^{h,l}_s2   = (__^ + γ I_D)^-1/2(__^ - _ S_s S_s^_^)(__^ + γ I_D)^-1/22
       = (__^ + γ I_D)^-1/2(__^ - _ A A^_^ - _ B B^_^)(__^ + γ I_D)^-1/22≤ε

and we know that this implies

    _ A A^_^≼__^ + ε(__^ + γ I_D) - _ B B^_^≼__^ + ε(__^ + γ I_D).

Plugging it in the initial definition,

    _s,i   = (1-ε) ϕ_i^(_ A A^_^ + (1+ε)γ I_D)^-1 ϕ_i
       ≥ (1 - ) ϕ_i^(__^ + ε(__^ + γ I_D) + (1+ε)γ I_D)^-1 ϕ_i
       = (1 - )1/1+2ε ϕ_i^(__^ + γ I_D)^-1 ϕ_i
    ≥1-ε/1+2ετ_h,i≥τ_h,i/α.


This result is weaker then Lemmas <ref> and <ref>,
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}_s=1^h 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

    [ Y_s _s-1]
        = Y_s-1⇔[ X_{s,i,j}_s-1] =  0.

We begin by inspecting the conditional random variable X_{s,i,j} | _s-1. Given the definition of X_{s,i,j}, the conditioning on _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 _s-1 which are contained in _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
[ X_{s,i,j}_s-1] =  0
is trivially satisfied.
On the other hand, if Y_s-1≤ε we have

    _u_s,i,j   [1/q(z_s-1,i,j/p_s-1,i,j - z_s,i,j/p_s,i,j)_i_i^_s-1]
       = 1/q( z_s-1,i,j/p_s-1,i,j -z_s-1,i,j/p_s,i,j[{ u_s,i,j≤p_s,i,j/p_s-1,i,j}_s-1])_i_i^
       = 1/q(z_s-1,i,j/p_s-1,i,j - z_s-1,i,j/p_s,i,jp_s,i,j/p_s-1,i,j)_i_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_i^≤1/q1/p_s,i,j_i_i^
    = 1/q1/p_s,i,jτ_h,i≤1/qα/τ_h,iτ_h,i
    =α/q R,

where we used Lemma <ref> 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 <cit.> to bound the norm of Y.


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 ∑_j=1^k [  X_j^2 { X_s}_s=0^j-1],
      for k = 1, 2, 3, ….

Then, for all ε≥ 0 and σ^2 > 0,

    ( ∃ 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

    (  Y_h≥ε)
    ≤(  Y_h≥ε)
       = (  Y_h≥ε∩ W_h≤σ^2)
    + (  Y_h≥ε∩ W_h≥σ^2)
       ≤(  Y_h≥ε∩ W_h≤σ^2)_
        + (  W_h≥σ^2)_.

Using the bound on  X_{s,i,j}_2, we can directly apply Proposition <ref> to bound  for any fixed σ^2.
To bound the part , we use the following lemma, proved later in Sec. <ref>.


    
    (  W_h≥6α/q)
            ≤ n ·exp{ - 2q/α}



Combining Prop. <ref> with σ^2 = 6α/q, Lem <ref>, the fact that 2ε/3 ≤ 1 and the value used by Alg. <ref> q = 39αlog(2n/δ)/ε^2 we obtain

    (  P_{h,l} -  P_{h,l}_2 ≥ε)
        = (  Y_h≥ε)
    ≤(  Y_h≥ε∩ W_h≤σ^2)
    + (  W_h≥σ^2)
       ≤ 2ν·exp{-ε^2q/α(1/12+2ε/3)}
    + n ·exp{- 2q/α}
       ≤ 3n ·exp{ - ε^2/13αq}
    = 3n ·exp{ - 3log(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.

 




 §.§ Proof of Lemma <ref> (bound on predictable quadratic variation)




cnt-lem-quad-variation



Step cnt-lem-quad-variationcnt-lem-quad-variation (a preliminary bound).
We start by writing out W_r for the process Y_s,

    W_r =1/q^2∑_{s,i,j}≤ r[(z_s-1,i,j/p_s-1,i,j - z_s,i,j/p_s,i,j)^2_{s,i,j}-1]_i_i^_i_i^.

We rewrite the expectation terms in the equation above as

    [ 
    (z_s-1,i,j/p_s-1,i,j - z_s,i,j/p_s,i,j)^2_{s,i,j}-1]
       = [z_s-1,i,j^2/p_s-1,i,j^2 -2 z_s-1,i,j/p_s-1,i,jz_s,i,j/p_s,i,j +z_s,i,j^2/p_s,i,j^2_{s,i,j}-1]
       (a)=[z_s-1,i,j^2/p_s-1,i,j^2 -2 z_s-1,i,j/p_s-1,i,jz_s,i,j/p_s,i,j +z_s,i,j^2/p_s,i,j^2_s-1]
       = z_s-1,i,j^2/p_s-1,i,j^2 -2 z_s-1,i,j/p_s-1,i,j1/p_s,i,j[z_s,i,j_s-1] +1/p_s,i,j^2[z_s,i,j^2 _s-1]
       (b)=z_s-1,i,j/p_s-1,i,j^2
    -2 z_s-1,i,j/p_s-1,i,jz_s-1,i,j/p_s-1,i,j
    +1/p_s,i,j^2[z_s,i,j_s-1]
       =1/p_s,i,j^2[z_s,i,j_s-1] - z_s-1,i,j/p_s-1,i,j^2
       (c)=1/p_s,i,jz_s-1,i,j/p_s-1,i,j - z_s-1,i,j/p_s-1,i,j^2
    =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∑_j=1^q∑_i=1^∑_s=1^h z_s-1,i,j/p_s-1,i,j(1/p_s,i,j - 1/p_s-1,i,j)_i_i^_i_i^.


We can now upper-bound W_h as

    W_h   ≼1/q^2∑_j=1^q∑_i=1^∑_s=1^h z_s-1,i,j/p_s-1,i,j(1/p_s,i,j - 1/p_s-1,i,j)_i_i^_i_i^
       = 1/q^2∑_j=1^q∑_i=1^(z_h,i,j/p_h,i,j^2 - z_h,i,j/p_h,i,j^2 + ∑_s=1^h z_s-1,i,j/p_s-1,i,j(1/p_s,i,j - 1/p_s-1,i,j))_i_i^_i_i^
       = 1/q^2∑_j=1^q∑_i=1^(z_h,i,j/p_h,i,j^2 + (∑_s=1^h -z_s,i,j/p_s,i,j^2 + z_s-1,i,j/p_s,i,jp_s-1,i,j) - z_0,i,j/p_0,i,j^2)_i_i^_i_i^
       ≼1/q^2∑_j=1^q∑_i=1^(z_h,i,j/p_h,i,j^2 + (∑_s=1^hz_s-1,i,j/p_s,i,jp_s-1,i,j - z_s,i,j/p_s,i,jp_s-1,i,j))_i_i^_i_i^
       = 1/q^2∑_j=1^q∑_i=1^(z_h,i,j/p_h,i,j^2 + ∑_s=1^hz_s-1,i,j(1 - z_s,i,j)/p_s,i,jp_s-1,i,j)_i_i^_i_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 z_s-1,i,j(1 - z_s,i,j)/p_s,i,jp_s-1,i,j = z_s-1,i,j(1 - z_s,i,j)/p_s-1,i,j^2, 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(1-z_s,i,j) terms can be active for s ∈ [h] and thus

    W_h   ≼1/q^2∑_j=1^q∑_i=1^(z_h,i,j/p_h,i,j^2 + ∑_s=1^hz_s-1,i,j(1 - z_s,i,j)/p_s-1,i,j^2)_i_i^_i_i^
       =1/q^2∑_j=1^q∑_i=1^(max{max_s=1,…,h{z_s-1,i,j(1 - z_s,i,j)/p_s-1,i,j^2} ; z_h,i,j/p_h,i,j^2})_i_i^_i_i^
       = 1/q^2∑_j=1^q∑_i=1^_i_i^_i_i^(max_s=0,…,h{z_s,i,j/p_s,i,j^2})·




Step cnt-lem-quad-variationcnt-lem-quad-variation (introduction of a stochastically dominant process).
We want to study max_s=0,…,h{z_s,i,j/p_s,i,j^2}.
To simplify notation, we will consider max_s=0,…,h{z_s,i,j/p_s,i,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,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 z_s,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}_s=0^h 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,[
Notice that unlike z_s,i,j, w_s,i,j is no longer _s-measurable but it is '_s-measurable, where 

    '_{s,i,j} = { u_s',i',j', ∀{s',i',j'}≤{s,i,j}}∪{  w_s,i,j}
    = _{s,i,j}∪{  w_s,i,j}.

] 

    (1/w_s,i,j≤ a _s)
     = 
    0     for    a < 1/p_s,i,j
    
    1-1/p_s,i,ja    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. <ref>
to guarantee that 1/p_s,i,j≤ a < α/p_h,i.
Note that the distribution of 1/w_s,i,j conditioned on _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
_s-1 (computed using the previous iteration),
we have that  

    (1/w_s,i,j≤ a _s)
    = (1/w_s,i,j≤ a _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
_s-1. It follows that given  _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. <ref>, which illustrates the dependence between
the various variables.

Finally for  the special case w_0,i,j the definition above reduces to

    (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
    ,

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 cnt-lem-quad-variationcnt-lem-quad-variation 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,

    (A ≥ a) ≥(B ≥ a) ⇔(A ≤ a) ≤(B ≤ a).

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

    (A ≥ a) ≥(B ≥ a) [A] ≥[B],

while the reverse (A dominates B, if [A] ≥[B]) is not
true in general.
Following this definition of  stochastic dominance, our goal is to prove

    (max_s=0^hz_s,i,j/p_s,i,j≤ a)
    ≥(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., _h. From the definition of w_h,i,j, we know that,
w.p. 1, 1/p_h,i≤ 1/w_h,i,j.
Therefore,

    (max_s=0,…,hz_s,i,j/p_s,i,j≤ a)
    ≥(max{max_s=0,…,h-1z_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 _k-1. Since u_k,i,j and
w_k,i,j are independent given _k-1, we have

    ( z_k,i,j/w_k,i,j≤ a _k-1)
       =( {u_k,i,j≤p_k,i,j/p_k-1,i,j}1/w_k,i,j≤ a _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,ja) = 1-1/p_k-1,i,ja    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,ja    for     1/p_k-1,i,j≤ a < 1/p_k,i,j
    
    1-1/p_k-1,i,ja    for    1/p_k,i,j≤ a < α/p_h,i
    
    1     for   α/p_h,i≤ a 
    
       =(1/w_k-1,i,j≤ a _k-2)
    =(1/w_k-1,i,j≤ a _k-1) ,

where the inequality is also represented in Fig. <ref>.

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. <ref> 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.

    (max{max_s=0 … k-1z_s,i,j/p_s,i,j ; z_k,i,j/w_k,i,j}≤ a )
    =__k-1[(max{max_s=0 … k-1z_s,i,j/p_s,i,j ; z_k,i,j/w_k,i,j}≤ a _k-1)]
       (a)≥__k-1[(max{max_s=0 … k-1z_s,i,j/p_s,i,j ; z_k-1,i,j/w_k-1,i,j}≤ a _k-1)]
       =__k-1[(max{max_s=0 … k-2z_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 _k-1)]
       =__k-1[(max{max_s=0 … k-2z_s,i,j/p_s,i,j ;z_k-1,i,jmax{1/p_k-1,i,j; 1/w_k-1,i,j}}≤ a _k-1)]
       (b)=__k-1[(max{max_s=0 … k-2z_s,i,j/p_s,i,j ;z_k-1,i,j/w_k-1,i,j}≤ a _k-1)]
    = (max{max_s=0 … k-2z_s,i,j/p_s,i,j ; z_k-1,i,j/w_k-1,i,j}≤ a ),
where in (a), given _k-1, everything is fixed except u_k,i,j and
w_k,i,j and we can use the stochastic dominance in (<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,

    (max_s=0,…,hz_s,i,j/p_s,i,j≤ a)
    ≥(max{z_0,i,j/p_0,i,j ; z_0,i,j/w_0,i,j}≤ a )
    ≥(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 (<ref>).

Step cnt-lem-quad-variationcnt-lem-quad-variation (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_h',i,j^max
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_h',i,j^max·

We first notice that while Y_h is not
necessarily PSD, W_h is a sum of PSD matrices.
Introducing the function Λ({1/p_h,i,j^max}_i,j) we can restate Eq. <ref> as

    W_h
    =λ_max( W_h)
    ≤Λ({1/p_h,i,j^max}_i,j) λ_max(1/q^2∑_j=1^q∑_i=1^(1/p_h,i,j^max)^2_i_i^_i_i^).

In Step 4, we showed that 1/p_h,i,j^max is stochastically
dominated by 1/w_0,i,j for every i and j. In order to bound
Λ({1/p_h,i,j^max}_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_h,i,j^max}_i,j) as

    λ_max(1/q^2∑_j=1^q∑_i=1^(1/p_h,i,j^max)^2_i_i^_i_i^)
    = max_ x :  x = 1 x^(1/q^2∑_j=1^q∑_i=1^(1/p_h,i,j^max)^2_i_i^_i_i^) x
       = max_ x :  x = 11/q^2∑_j=1^q∑_i=1^(1/p_h,i,j^max)^2_i_2^2 x^_i_i^ x
    = max_ x :  x = 11/q^2∑_j=1^q∑_i=1^(1/p_h,i,j^max)^2(_i_2_i^ x)^2.

From this reformulation, it is easy to see that, because 1/p_h,i,j^max
is strictly positive,
the function Λ({1/p_h,i,j^max}_i,j) is
monotonically increasing w.r.t. the individual 1/p_h,i,j^max, or
in other words that increasing an 1/p_h,i,j^max 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) max_ x :  x = 11/q^2∑_j=1^q∑_i=1^(1/w_0,i,j)^2(_i_2_i^ x)^2,

we now need to prove the stochastic dominance of Λ({1/w_0,i,j}_i,j) over
Λ({1/p_h,i,j^max}_i,j).
Using the definition of 1/p_h,i,j^max, w_h,i,j, and the
monotonicity of Λ we have

    (Λ({1/p_h,i,j^max}_i,j) ≤ a)
       =(Λ({max{max_s=0,…,h-1z_s,i,j/p_s,i,j ; z_h,i,j/p_h,i,j}}_i,j) ≤ a )
       ≥(Λ({max{max_s=0,…,h-1z_s,i,j/p_s,i,j ; z_h,i,j/w_h,i,j}}_i,j) ≤ a ).

Now pick 1≤ k ≤ h, for a fixed _k-1, 1/p_k-1,i,j^max is a constant
and max{1/p_k,i,j^max; x} is a monotonically
increasing function in x, making Λ(max{1/p_k,i,j^max; x}) also an increasing function. Therefore, we have

    (Λ({max{1/p_k-1,i,j^max ; z_k,e,j/w_k,i,j}}_i,j) ≤ a )
    =__k-1[(Λ({max{1/p_k-1,i,j^max ; z_k,e,j/w_k,i,j}}_i,j) ≤ a _k-1)]
       (a)≥__k-1[(Λ({max{1/p_k-1,i,j^max ; z_k-1,i,j/w_k-1,i,j}}_i,j) ≤ a _k-1)]
       (b)=__k-1[(Λ({max{1/p_k-2,i,j^max ; z_k-1,i,j/w_k-1,i,j}}_i,j) ≤ a _k-1)],

where inequality (a) follows from the fact that stochastic dominance is preserved
by monotonically increasing functions <cit.>, such as
Λ, combined with the fact that for
a fixed _k-1 the variables z_k,i,j and w_k,i,j are all
independent and (b) from the definition of 1/p_k-1,i,j^max 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

    ( W_h≥σ^2)
       ≤(Λ({1/p_h,i,j^max}_i,j) ≥σ^2)
    ≤(λ_max(1/q^2∑_j=1^q∑_i=1^(1/w_0,i,j)^2_i_i^_i_i^) ≥σ^2).


Step cnt-lem-quad-variationcnt-lem-quad-variation (concentration inequality).
Since all w_0,i,j are (unconditionally) independent from each other, we can apply the following theorem.

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([ V]) = λ_max(∑_k[ X_k]).

Then, for all h ≥ 0,

    (  λ_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^21/w_0,i,j_i_i^_i_i^≼1/q^2α^2/p_h,i^2_i_i^_i_i^≼1/q^2α^2/p_h,i^2_i_i^^2 I
    ≼α^2/q^2 I,

where the first inequality follows from the definition of w_0,i,j in Eq. <ref>, the second
from the PSD ordering, and the third from the definition of _i_i^.

Therefore, we can use L α^2/q^2 for the purpose of Prop. <ref>. We need now to compute [ 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_0,i,j^2 as

    F_1/w_0,i,j^2(a) = (1/w_0,i,j^2≤ a) = (1/w_0,i,j≤√(a))
     = 
    0     for    a < 1
    
    1-1/√(a)    for    1 ≤ a < α^2/p_h,i^2
    
    1     for   α^2/p_h,i^2 ≤ a
    .

Since (1/w_0,i,j^2 ≥ 0) = 1, we have that 


    [1/w_0,i,j]
    = ∫_a=0^∞[1-F_1/w_0,i,j(a)] a
       = ∫_a=0^1 (1 - F_1/w_0,i,j(a)) a + ∫_a=1^α^2/p_h,i^2(1 - F_1/w_0,i,j(a)) a + ∫_a=α^2/p_h,i^2^∞(1 - F_1/w_0,i,j(a)) a
       = ∫_a=0^1 (1 - 0) a + ∫_a=1^α^2/p_h,i^2(1 - (1 - 1/√(a)) ) a +∫_a=α^2/p_h,i^2^∞(1-1) a 
       = ∫_a=0^1  a + ∫_a=1^α^2/p_h,i^21/√(a) a
    = 1+[2√(a)]_1^α^2/p_h,i^2
    = 2α/p_h,i - 1.

Therefore,

    μ_max   = λ_max([ V]) = λ_max(∑_{i,j}[ X_{i,j}])
    = λ_max(1/q^2∑_j=1^q∑_i=1^[1/w_0,i,j^2]_i_i^_i_i^)
       = λ_max(1/q∑_i=1^(2α/p_h,i - 1)p_h,i_i_i^)
    ≤λ_max(2α/q∑_i=1^_i_i^)
    = 2α/qλ_max( P)
    ≤2α/q L.

Therefore, selecting h = 2, σ^2=6α/q and applying Prop. <ref> we have

    ( W_h≥σ^2)
       ≤(λ_max(1/q^2∑_j=1^q∑_i=1^1/w_0,i,j^2_i_i^_i_i^) ≥ (1+2)2α/q)
       ≤·exp{ - 2α/qq^2/α^2(3log(3) - 2) }≤ n ·exp{ - 2q/α}·


 




 §.§ Space complexity bound



Denote with A the event A = {∀  h' ∈{1, …, h} :   P^h' -  P^h'_2 ≤ε}, and again  = |_{h,l}|. Letting q = |_{h,l}| = ∑_i=1^ q_h,i = ∑_j=1^q∑_i=1^  z_h,i,j be the random number
of points in _{h,l}, we reformulate

    (|_{h,l}| ≥ 3qγ_{h,l}∩{∀  h' ∈{1, …, h} :  ( P^h' -  P^h'_2 ≤ε) ≤ε})
       = (|_{h,l}| ≥ 3qγ_{h,l}∩ A)
    = (∑_j=1^q∑_i=1^  z_h,i,j≥ 3qγ_{h,l}∩ A)
       =(  ∑_j=1^q∑_i=1^ z_h,i,j≥ 3qγ_{h,l} A )
    (A ).

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. <ref>,
we are going to find a stochastic variable to dominate z_h,i,j.
Denoting with u'_s,i,j∼𝒰(0,1) a uniform random variable,
we will define w'_s,i,j as

    w'_s,i,j|_{s,i',j'} = w'_s,i,j|_s-2{u'_s,i,j≤p_h,i/p_s-1,i}∼ℬ(p_h,i/p_s-1,i)

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

    ( ∑_j=1^q∑_i=1^ z_h,i,j≥ g  A)
    = __t-1|A[( ∑_j=1^q∑_i=1^{u_h,i,j≤p_h,i/p_t-1,i}z_h-1,i,j≥ g _t-1∩ A)]
       ≤__t-1|A[( ∑_j=1^q∑_i=1^{u'_h,i,j≤p_h,i/p_t-1,i}z_h-1,i,j≥ g _t-1∩ A)]
    = ( ∑_j=1^q∑_i=1^ w'_h,i,jz_h-1,i,j≥ g  A),

where we used the fact that conditioned on A, _{h,l} is accurate _{h,l}, which guarantees that p_h,i≤ p_h,i.
Plugging this in the previous bound,

    ( ∑_j=1^q∑_i=1^ z_h,i,j≥ g  A )
    (A )
    ≤( ∑_j=1^q∑_i=1^ w'_h,i,jz_h-1,i,j≥ g ∩ A)
    ≤( ∑_j=1^q∑_i=1^ w'_h,i,jz_h-1,i,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.

    ( ∑_j=1^q∑_i=1^ w'_s,i,jz_s-1,i,j≥ g )
    = __s-2[( ∑_j=1^q∑_i=1^{u'_s,i,j≤p_h,i/p_s-1,i}z_s-1,i,j≥ g _s-2)]
       = __s-2[( ∑_j=1^q∑_i=1^{u'_s,i,j≤p_h,i/p_s-1,i}{u_s-1,i,j≤p_s-1,i/p_s-2,i}z_s-2,i,j≥ g _s-2)]
       = __s-2[( ∑_j=1^q∑_i=1^{u'_s-1,i,j≤p_h,i/p_s-2,i}z_s-2,i,j≥ g _s-2)]
    = ( ∑_j=1^q∑_i=1^ w'_s-1,i,jz_s-2,i,j≥ g )

Applying this repeatedly from s=h to s=2 we have,

    ( ∑_j=1^q∑_i=1^ w'_h,i,jz_h-1,i,j≥ g )
    =( ∑_j=1^q∑_i=1^ w'_1,i,jz_0,i,j≥ g )
    =( ∑_j=1^q∑_i=1^ w'_1,i,j≥ g ).

Now, all the w'_1,i,j are independent Bernoulli random variables,
and we can bound their sum with a Hoeffding-like bound using Markov inequality,

    ( ∑_j=1^q∑_i=1^ w'_1,i,j≥ g ) = inf_θ > 0( e^∑_j=1^q∑_i=1^θ w'_1,i,j≥ e^θ g)
       ≤inf_θ > 0[ e^∑_j=1^q∑_i=1^θ w'_1,i,j]/e^θ g
    = inf_θ > 0[ ∏_j=1^q∏_i=1^ e^θ w'_1,i,j]/e^θ g
    = inf_θ > 0∏_j=1^q∏_i=1^[ e^θ w'_1,i,j]/e^θ g
       =inf_θ > 0∏_j=1^q∏_i=1^ (p_h,i e^θ + (1-p_h,i))/e^θ g
    = inf_θ > 0∏_j=1^q∏_i=1^ (1+p_h,i( e^θ -1))/e^θ g
       ≤inf_θ > 0∏_j=1^q∏_i=1^ e^p_h,i( e^θ -1)/e^θ g≤inf_θ > 0e^q( e^θ -1)∑_i=1^ p_h,i/e^θ g
    = inf_θ > 0 e^(γ_{h,l}q(e^θ - 1) - θ g)≤inf_θ > 0 e^(γ_{h,l}q(e^θ - 1) - θ g),

where we use the fact that 1 + x ≤ e^x, w'_1,i,j∼ℬ(p_h,i) and by Def. <ref>, ∑_i=1^ p_h,i = ∑_i=1^τ_h,i = γ_{h,l}.
The choice of θ minimizing the previous expression is obtained as

    d/dθe^(qγ_{h,l}(e^θ - 1) - θ g )
    = e^(qγ_{h,l}(e^θ - 1) - θ g )(qγ_{h,l}e^θ - g) = 0,

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

    inf_θexp{qγ_{h,l}(e^θ - 1) - θ g)}   = exp{ g  - qγ_{h,l} -  g log(g/qγ_{h,l})}
       = exp{ -g(log(g/qγ_{h,l}) - 1)} e^- qγ_{h,l},

and choosing g = 3qγ_{h,l}, we conclude our proof.
 

§ APPLICATIONS


    
    _n = _n _n(_n^_n _n + γ I_n)^-1_n^_n
       = _n^_n _n(_n^_n^_n _n + γ I_n)^-1_n^_n^_n
       = _n^_n _n_n^_n^(_n _n_n^_n^ + γ I_D)^-1_n
       = _n^(_n _n_n^_n^ + γ I_D - γ I_D)(_n _n_n^_n^ + γ I_D)^-1_n
       = _n^( I_D - γ(_n _n_n^_n^ + γ I_D)^-1)_n
       = _n^ I_D_n - γ_n^(_n _n_n^_n^ + γ I_D)^-1_n
       = _n - γ_n^(_n _n_n^_n^ + γ I_D)^-1_n

    From Lem. <ref>, and the fact that _n is ε-approximate we have that
    
    (_t_t^ + γ I_D)^-1/2_t( I_t -  S_s S_s^)_t^(_t_t^ + γ I_D)^-1/22≤ε,

    which implies
    
    _t_t^ - _t S_s S_s^_t^≼ε(_t_t^ + γ I_D)

    and
    
    _t_t^ - ε(_t_t^ + γ I_D) ≼_t S_s S_s^_t^.

    Therefore,
    
    _n - _n    = γ_n^(_n _n_n^_n^ + γ I_D)^-1_n
         ≼γ_n^(_t_t^ - ε(_t_t^ + γ I_D) + γ I_D)^-1_n
       = γ/1-ε_n^(_t_t^ + γ I_D)^-1_n
         = γ/1-ε_t(_t + γ I)^-1≼γ/1-ε I.


 

