Year: 2021

Abstract:
In distribution compression, one aims to accurately summarize a probability distribution P using a small number of representative points. Near-optimal thinning procedures achieve this goal by sampling n points from a Markov chain and identifying √ n points with Õ(1/ √ n) distributional discrepancy to P. Unfortunately, these same algorithms suffer from quadratic or super-quadratic runtime in the sample size n. To address this deficiency, we introduce a simple metaprocedure—Compress++—for speeding up any input thinning algorithm while suffering at most a factor of four in error. When combined with the quadratictime kernel halving and kernel thinning algorithms of Dwivedi and Mackey (2021), Compress++ delivers √ n points with O( √ log n/n) integration error and better-than-Monte-Carlo maximum mean discrepancy in O(n log n) time and O( √ n log(n)) space. Moreover, Compress++ enjoys the same near-linear runtime given any quadratic-time input and reduces the runtime of super-quadratic algorithms by a square-root factor. In our benchmarks with high-dimensional Monte Carlo samples and long-running Markov chains targeting challenging differential equation posteriors, Compress++ matches or nearly matches the accuracy of its input algorithm in orders of magnitude less time.

text: √ n points with Õ(1/ √ n) distributional discrepancy to P. Unfortu-
nately, these same algorithms suffer from quadratic or super-quadratic runtime in the sample size n. To address this deficiency, we introduce a simple metaprocedure—Compress++—for speeding up any input thinning algorithm while suffering at most a factor of four in error. When combined with the quadratictime kernel halving and kernel thinning algorithms of Dwivedi and Mackey (2021), Compress++ delivers √ n points with O( √ log n/n) integration error and better-than-Monte-Carlo maximum mean discrepancy in O(n log2 n) time and O( √ n log2(n)) space. Moreover, Compress++ enjoys the same near-linear runtime given any quadratic-time input and reduces the runtime of super-quadratic algorithms by a square-root factor. In our benchmarks with high-dimensional Monte Carlo samples and long-running Markov chains targeting challenging differential equation posteriors, Compress++ matches or nearly matches the accuracy of its input algorithm in orders of magnitude less time.

heading: 1 INTRODUCTION
text: Distribution compression—constructing a concise summary of a probability distribution—is at the heart of many learning and inference tasks. For example, in Monte Carlo integration and Bayesian inference, n representative points are sampled i.i.d. or from a Markov chain to approximate expectations and quantify uncertainty under an intractable (posterior) distribution (Robert & Casella, 1999). However, these standard sampling strategies represent a bottleneck in computationally-demanding settings due to their slow root-n Monte Carlo error rate. For instance, the Monte Carlo estimate Pnf , 1n ∑n i=1 f(xi) of an unknown expectation Pf , EX∼P[f(X)] based on n i.i.d. points has Θ(n− 1 2 ) integration error |Pnf − Pf |, requiring n = 10000 points for 1% relative error and n = 106 points for 0.1% error. Such bloated sample representations preclude downstream applications with critically expensive function evaluations like computational cardiology, where a 1000-CPU-hour tissue or organ simulation is required for each sample point (Niederer et al., 2011; Augustin et al., 2016; Strocchi et al., 2020), or expert data annotation which can be monetarily expensive Beaugnon et al. (2017), and is sometimes coined as the human bottleneck of machine learning.
To restore the feasibility of such critically expensive tasks, it is common to thin down the initial sequence of points to a produce a much smaller coreset. The standard thinning approach, select every t-th sample point (Owen, 2017), is simple to implement but often leads to an substantial increase in error: e.g., standard thinning n points from a fast-mixing Markov chain yields Ω(n− 1 4 ) error when n 1 2 points are returned. Recently, Dwivedi & Mackey (2021b;a) introduced a more effective alternative, kernel thinning, that provides near optimal Õ(n− 12 ) error when compressing n points down to size n 1 2 . While practical for moderate sample sizes, the runtime of this algorithm scales quadratically with the input size n, making its execution prohibitive for very large input sizes. Our goal is to significantly improve the runtime of such compression algorithms while providing comparable error guarantees.
Problem setup Given a sequence Sin of n input points summarizing a target distribution P, our aim is to identify a high quality coreset Sout of size √ n in time nearly linear in n. We measure coreset
quality via its integration error |Pf −Poutf | , |Pf − 1|Sout| ∑ x∈Sout f(x)| for functions f in the reproducing kernel Hilbert space (RKHS)H induced by a given kernel k (Berlinet & Thomas-Agnan, 2011). We consider both single function error and kernel maximum mean discrepancy (MMD, Gretton et al., 2012), the worst-case integration error over the unit RKHS norm ball:
MMDk(P,Pout) , sup‖f‖k≤1|Pf − Poutf | = ‖Pk− Poutk‖k.
Our contributions We introduce a new simple meta procedure—COMPRESS++—that significantly speeds up a generic thinning algorithm while simultaneously inheriting the error guarantees of its input up to a factor of four. A direct application of COMPRESS++ to kernel thinning improves its quadratic O(n2) runtime to near linear O(n log2 n) time while maintaining the error guarantees up to a factor four. Since the Õ(n− 12 ) KT MMD guarantees of Dwivedi & Mackey (2021b) match the Ω(n− 1 2 ) minimax lower bounds of Tolstikhin et al. (2017); Phillips & Tai (2020) up to factors
of √
log(n) and constants depending on d, KT-COMPRESS++ also provides near-optimal MMD compression for a wide range of kernels and distributions P. Moreover, the practical gains from applying COMPRESS++ are substantial: KT thins 65, 000 points in 10 dimensions in 20m, while KTCOMPRESS++ needs only 1.5m; KT takes more than a day to thin 250, 000 points in 100 dimensions, while KT-COMPRESS++ takes less than an hour (a 32× speed-up). COMPRESS++ can also be directly combined with any thinning algorithm, even those that have suboptimal guarantees but often perform well in practice, like kernel herding (Chen et al., 2010), support points (Mak & Joseph, 2018), Stein points (MCMC) (Chen et al., 2018; 2019), and Stein thinning (Riabiz et al., 2020a), all of which run in Ω(n2) time. As a demonstration, we combine COMPRESS++ with the popular kernel herding algorithm and observe 10-60× speed-ups. In all of our experiments, COMPRESS++ leads to minimal loss in accuracy and, surprisingly, even improves upon herding accuracy for lower-dimensional problems.
Overview of COMPRESS++ To define COMPRESS++, we first introduce COMPRESS: A simple and elegant recursive strategy takes in a halving algorithm and provides an intermediate thinned coreset in significantly faster runtime. COMPRESS divides the coreset into four parts, recursively applies COMPRESS to each, combines the output of each and halves the resulting coreset to get the output. Then COMPRESS++ with a thinning algorithm ALG is defined as follows: In stage one, use COMPRESS with ALG instantiated for 2-thinning, to obtain an intermediate coreset of size 2g √ n for a suitable parameter g, and then in stage two apply ALG directly to further thin it down to √ n points. In this manner, ALG is ever applied to coresets of size 2g √ n or smaller thereby improving the runtime. The parameter g is chosen based on the run-time and known guarantees of the underlying thinning algorithm, but a default choice of 4 can be shown to be optimal in several theoretical settings, and it provided competitive performance across all our experiments. Overall, with n input points and a thinning algorithm with runtime nα, COMPRESS++ uses a outputs √ n
points in time Õ(nα/2) with errors similar to those provided by the input thinning algorithm. Notation We write SALG for the coreset outputted by an algorithm ALG and extend our MMD definition for point sequences (S1,S2) with empirical distributions (Q1,Q2) via MMDk(S1,S2) , MMDk(Q1,Q2) and MMDk(P,S1) , MMDk(P,Q1). We use a - b and a % b to mean a = O(b) and a = Ω(b) and use -k to denote that the underlying constants depend on k.
Basic definitions We use H to denote an inner-product space endowed with the inner-product 〈·, ·〉H. We start with the definition of sub-gamma random variables (Boucheron et al., 2013).
Definition 1 (f -Sub-Gamma) A random variable X is said to be sub-gamma on the right with parameters ( σ2, c ) denoted by Γ+ ( σ2, c ) if for all 0 < λ < 1c , we have E [ exp(λX) ] ≤
exp ( λ2σ2
2(1−cλ)
) . For any f ∈ H, we say that a random variable G on H is f -sub-gamma on the
right, denoted by Γf+ ( σ2, c ) , if the random variable 〈f,G〉H ∈ Γ+ ( σ2, c ) .
As a reminder, we note that c = 0 yields a sub-Gaussian tails on the right. For X ∈ Γ+ ( σ2, c ) , Boucheron et al. (2013, Section 2.4) shows that for t ≥ 0, and δ ∈ (0, 1], we have
P [ X > √ 2σ2t+ ct ] ≤ e−t, or equivalently P [ X ≤ √ 2σ2 log( 1δ ) + c log( 1 δ ) ] ≥ 1− δ. (1)
A property that we frequently use is that the sub-gamma property is suitably closed under multiplication and addition of sub-gamma random variables. A discussion is deferred to App. A.
Definition 2 (Thinning and halving algorithms) Consider an algorithm ALG that takes as input a point sequence Sin of size n and returns a (possibly random) subsequence SALG ⊂ Sin of size nout. We say that ALG is αn-thinning if nout = n/αn whenever n/αn ∈ N. When αn = 2, we say that ALG is a Halving algorithm, and when αn = √ n, we say ALG is root-thinning.
For a thinning algorithm ALG, we associate a kernel discrepancy embedding that measures the approximation quality for the n input points Sin provided by the nout points SALG output by ALG:
ψALG(Sin) , ∑ x∈Sin k(x, ·)− n nout ∑ x∈SALG k(x, ·). (2)
Letting PS denote the empirical distribution of S, the reproducing property of k yieldsthat
PSinf − PSALGf = 1n 〈f, ψALG〉k, for any f ∈ H, and MMD(Sin,SALG) = 1 n‖ψALG‖k. (3)
Next, we define the notion of a sub-gamma thinning algorithm via the object ψALG (2).
Definition 3 (Sub-Gamma Thinning Algorithm) Given the functions f ∈ H, σ : N→ R+ and c : N→ R+, an αn-thinning algorithm ALG is called f -sub-gamma on the right with parameters (σ2, c), if ψALG(Sin) ∈ Γf+(σ2(n), c(n)) conditioned on Sin of size n.
When the algorithm is both Γf+ and Γ −f + -sub-gamma, (1) and (3) immmediately imply a high probability tail bound on the integration error |PSinf − PSALGf |.

heading: 2 COMPRESS
text: We now introduce our first meta procedure COMPRESS (which is also a building block for COMPRESS++). COMPRESS takes as input a halving algorithm HALVE and a oversampling factor g, and then given any input of size n, it outputs at thinned coreset of size 2g √ n.
The algorithm is extremely simple to implement: first, divide the input sequence into four subsequences of size n4 (in any manner the user chooses); second, recursively call COMPRESS on each subsequence to produce four coresets of size 2g−1 √ n; finally, call HALVE on the concatenation of those coresets to produce the final output of size 2g √ n. We denote this algorithm as COMPRESS(HALVE, g), and describe it formally in Alg. 1. COMPRESS can also be implemented in a streaming fashion but we defer the discussion to App. G. We remind that COMPRESS with g = 0 is a root thinning algorithm.
Algorithm 1: COMPRESS Input: Halving algorithm HALVE, oversampling factor g, point sequence Sin of size n if n = 4g then return Sin else
Partition Sin into four subsets {Si}4i=1 each of size n/4 for i = 1, 2, 3, 4 do S̃i ← COMPRESS(Si,HALVE, g) // return coresets of size 2g · √ n 4 end S̃ ← CONCATENATE(S̃1, S̃2, S̃3, S̃4) // coreset of size 2 · 2g · √ n S ← HALVE (k, S̃) // coreset of size 2g √ n
return S end

heading: 2.1 SUB-GAMMA PROPERTY AND RUNTIME GUARANTEES FOR COMPRESS
text: We now establish the sub-gamma property and runtime of COMPRESS in terms of the input halving algorithm HALVE. We measure computational complexity in terms of the dominant operations done by the halving algorithms (e.g., kernel computations with kernel thinning), and ignore the time taken for reading and writing to memory. See App. B for the proof. In the result, our notation highlights explicitly the dependence of σ2, c and rCP on g.
Theorem 1 (Sub-gamma property and runtime of COMPRESS) Let HALVE denote a Γf+halving algorithm with parameters (σ2HALVE, cHALVE), with running time rHALVE(n) on inputs of size n. Then with inputs of size n, COMPRESS(HALVE, g) outputs a coreset of size 2g √ n, and is Γf+thinning with parameters
σ2CP(n, g) , n 4g+1 ∑log4 n−g−1 i=0 σ 2 HALVE ( 2g+1−i √ n ) ≤ n4g+1 (log4 n− g)σ 2 HALVE ( 2g+1 √ n ) , (4)
cCP(n, g) , max0≤i≤log4 n−g−1 [ √ n 2g+i+1 cHALVE (√ n2g+1−i )] ,
and has the computational complexity bounded by rCP(n, g) = ∑log4 n−g−1 i=0 4 irHALVE(2 g+1 √ n4−i).
Remark 1 (Worst-case guarantee for COMPRESS) For any monotone σ2HALVE and cHALVE, we have σ
2 CP(n,0) n2 ≤ 1 4nσ 2 HALVE ( 2 √ n ) (log4 n) and cCP n ≤ 1 2 √ n cHALVE ( 2 √ n ) . Recall from (1), this
gives us a high probability bound on 1n 〈f, ψCP〉H. From (3), this bounds the integration error PSinf − PSCPf . Comparing this with the corresponding expressions for HALVE for compressing an input of size 2 √ n to √ n, we find that in the worst case the error from COMPRESS is only a logarithmic factor larger. In particular, this shows that the COMPRESS algorithm gives guarantees that are only a logarithmic factor increase in error compared other methods of combining HALVE to get a root-thinning algorithm (such as RECHALVE discussed in App. H).
Remark 2 (Explicit speed-up of runtimes) Let us demonstrate how COMPRESS speeds up a thinning algorithm. For ease of discussion set g = 0. First, we note that rCP(n) . log n · rHALVE ( 2 √ n ) . Thus, the running time is upper bounded by a logarithmic factor over the running HALVE on a much smaller input of size 2 √ n. Consequently, COMPRESS provides a significant improvement for any halving algorithm with super-linear runtime. Unrolling the recursion for explicit polynomial runtime, we find that
• rHALVE(n) = Cn2 =⇒ rCP(n) = C(log n− g)4g+1n
• rHALVE(n) = Cnk for k > 2 =⇒ rCP = Cnk/22kg+k
Note that in both these setting, we get an approximately quadratic speed up. In particular, for kernel thinning the speed up is from O(n2) to O(n log n), which also provides significant speed up in practice as our experiments later demonstrate.
Remark 3 (Explicit bounds on sub-gamma parameters) Unrolling the recursion (4), we now derive a explicit expressions for σ2CP(n, g) under the following growth conditions on σ 2 HALVE:
• σ2HALVE(n) = C log4 n =⇒ σ2CP(n, g) = C n4g+1 · 1 4 (log4 n+ 3g + 3)(log4 n− g)
• σ2HALVE(n) = C =⇒ σ2CP(n, g) = Cn4−g−1(log4 n− g) .
Example 1 (KT-SPLIT meets COMPRESS.) We now illustrate the gains of COMPRESS with the KT-SPLIT algorithm introduced in Dwivedi & Mackey (2021a, Alg1a.). Given a target kernel k, KT-SPLIT uses non-uniform randomness to recursively thin an input sequence Sin of size n to a sequence Sout of size n/α. Its runtime is bounded as O(n2), and thus Rem. 2 immediately yields that COMPRESS leads to an O( nlogn ) speed-up in runtime—thereby reducing the runtime of KT from days to 10 minutes for a million samples in 100 dimensions as seen in Fig. 1. Next, we show that the single function error of COMPRESS is at most √
log n times worse than that of KT-SPLIT. With α = 2, Dwivedi & Mackey (2021a, Thm 1, Rem. 1a) then implies that for any f in the RKHS of k, we have
1 n ∣∣∣∑x∈Sin f(x)− 2∑x∈Sout f(x)∣∣∣ ≤ 1n√ 32 log 43 ‖k‖∞‖f‖k(√log(2n) log( 1δ )+ log( 1δ ))
with probability 1 − δ. Inverting this bound, we obtain that
Pr [∣∣∣∑x∈Sin f(x)− 2∑x∈Sout f(x)∣∣∣ ≥ cnt+√2σ2nt] ≤ e−t for cn = √ 32 log 43 ‖k‖∞‖f‖k and σ2n = 4 log 4 3 log(2n) · ‖f‖ 2 k‖k‖∞. Applying (Boucheron et al., 2013, Theorem 2.3), we thus
conclude that KT-SPLIT (in the halving mode) is Γf+ ( 32σ2n + 64c 2 n, 8cn ) -sub-gamma.
Having established the sub-gamma property of KT-SPLIT halving, we can directly apply Theorem 1 to show that COMPRESS(KT-SPLIT, g) is also Γf+ sub-gamma with parameters
σ2CP ≤ 237 · n4g ‖f‖ 2 k‖k‖∞(log n+ 3g + 3)(log n− g), and cCP ≤ 22 √ ‖k‖∞‖f‖k,
which in turn implies that for g = 0, we have∣∣∣ 1n∑x∈Sin f(x)− 1√n∑x∈SCP f(x)∣∣∣ ≤ 22‖f‖k√‖k‖∞√n [√(log2 n+ 3 log n) log( 1δ )+ log( 1δ )] with probability 1− δ, which is at most O( √ log n) worse than the KT-SPLIT guarantee itself.

heading: 2.2 MMD GUARANTEES FOR COMPRESS
text: Next, we turn attention to establishing the worst-case function error, namely the MMD guarantees for COMPRESS. For brevity, we give simplified upper bounds on the error parameters. For the exact expressions and the proof, see App. D. Here SCP denotes the output of COMPRESS.
Theorem 2 (MMD Guarantees for COMPRESS) Let HALVE be a halving algorithm such that for any size n input it satisfies
P [ MMD(Sin,SHALVE) ≥ an + √ 2vnt+ cnt ] ≤ e−t for all t ≥ 0, and EψHALVE = 0.
Then, COMPRESS(HALVE, g) with size n input satisfies Pr [ MMD(Sin,SCP) ≥ σM,CP · √ 2t+ 2tRM,CP ] ≤ 2(n+ 1)e−t for all t ≥ 0, (5)
where RM,CP, √ 8vq+2a2q+4cq , and σM,CP,RM,CP √ log4 n−g with q, √ n2g+1.
Remark 4 (Symmetrization) We can convert any halving algorithm into one that satisfies the mean zero condition EψHALVE = 0 without impacting integration error, by symmetrization, i.e., by returning either the outputted half or its complement with equal probability.
Example 2 (Kernel thinning meets COMPRESS.) For MMD guarantees, we analyze COMPRESS with kernel thinning of Dwivedi & Mackey (2021b, Alg 1.) as the halving algorithm with g = 0, so that the overall algorithm is root-thinning, We call this instantiation KT-Compress, and compare and contrast it to KT for root thinning. 1
First, we note that KT has a quadratic runtime of order n2, so that COMPRESS once again provides a near linear O( nlogn ) speed-up by Rem. 2. Next, we claim that the MMD error for KT COMPRESS is at most O(log n) worse compared to KT, a mild degradation compared to the linear speed up. To prove this, first we note that for an output SKT of size n2m from KT, Dwivedi & Mackey (2021a, Thm. 2) guarantees that
MMD(Sin,SKT) -k 2 m n + 2m n √ log n · Mk(B2(R), 2 m n ) √ log(1/δ) + 2 m n √ log n log(1/δ), (6)
with probability at least 1 − δ, whereMk(A, ε) denotes the ε-L∞ cover of the unit ball of RKHS with domain restricted to the set A (Dwivedi & Mackey (2021a, Def. 1)), and R , maxx∈Sin‖x‖2. Inverting this bound, we conclude that KT as a halving algorithm (i.e., m = 1) satisfies the assumption of Theorem 2 with
an - 1n , vn -k logn n2 · Mk(B2(R), 2 n ), and cn -k √ logn n .
1Throughout this work, we use KT to refer to the KT variant which uses target kernel k in both KT-SPLIT and KT-SWAP steps, a variant referred to as target KT by Dwivedi & Mackey (2021a).
Now applying Theorem 2, we find that KT compress admits the tail bound (5) with RM,KT−CP -k √ log n · Mk(B2(R), 1√n ) and σM,KT−CP -k log n · √ Mk(B2(R), 1√n ).
Choosing t?=log( 2(n+1)δ ), we conclude that the √ n sized KT-compress output SKT−CP satisfies
MMD(Sin,SKT−CP) -k √ log n (√ log n · log(nδ ) + log( n δ ) ) 1√ n √ Mk(B2(R), 1√n ), (7)
with probability at least 1 − δ. On the other hand, with high probability the output from KT (when root thinning) itself satisfies 1√
n √ log n · Mk(B2(R), 1√n ) √ log(1/δ) where we note that
this bound follows from (6) by substituting m = 12 log2 n. Comparing (7) to the above bound, we conclude that KT-COMPRESS provides at most O(log n) worse MMD error than KT itself for general k, and Sin.

heading: 3 COMPRESS ++
text: In this section, we introduce another simpleyet effective meta procedure—COMPRESS++—that uses a generic thinning algorithm along with COMPRESS in two stages. Overall, it provides a significantly improved error at the cost of a mild logarithmic overhead in the running time when compared to COMPRESS. Algorithmically, COMPRESS++ takes as input an α-thinning algorithm (Def. 2) THIN(α) for α ∈ N.2 In the first stage, it runs COMPRESS with THIN(2) as the halving procedure and oversampling factor g set to a non-zero quantity using THIN(2) as the halving algorithm to reduce the size from n to √ n2g. Next, we run the THIN(2g) on the output of the first stage to thin the sequence down to size √ n. We formally present this algorithm in Alg. 2.
Algorithm 2: COMPRESS++ Input: Thinning algorithm THIN, Oversampling parameter g, Input point sequence Sin SCOMPRESS ← COMPRESS(THIN(2), g,Sin) // Coreset of size 2g √ n
SCOMPRESS++ ← THIN(2g,SCP) // Coreset of size √ n return SCP++

heading: 3.1 SUB-GAMMA PROPERTY AND RUNTIME GUARANTEES FOR COMPRESS++
text: First, we note the quality of the solution given by COMPRESS++ by expressing the Γf+ parameters in terms of those of THIN. The proof can be found in App. I.
Theorem 3 (Sub-gamma property and runtime of COMPRESS++) Let THIN(α) denote a Γf+α-thinning algorithm with parameters σ2THIN(n, α) and cTHIN(n, α), and running time rTHIN(n, α) for inputs of size n. Then COMPRESS++(THIN, g) is root-thinning and Γf+ with parameters
σ2CP++(n, g) = n 4g+1 ∑log4 n−g i=0 σ 2 THIN ( 2g+1−i √ n, 2 ) + n4gσ 2 THIN ( 2g √ n, 2g ) ≤ n4g+1 (log4 n− g)σ 2 THIN ( 2g+1 √ n, 2 ) + n4gσ 2 THIN ( 2g √ n, 2g
) cCP++(n, g) = max (√ n 2g cTHIN ( 2g √ n, 2g ) ,max0≤i≤log4 n−g−1 ( √ n 2g+i+1 cTHIN (√ n2g+1−i, 2 ))) ,
with its runtime bounded as rCP++(n) = ∑log4 n−g−1 i=0 4 irTHIN ( 2g+1−i √ n, 2 ) + rTHIN ( 2g √ n, 2g ) .
With some algebra with expressions above, one can see the advantage of COMPRESS++ over COMPRESS. Setting g = log4 log4 n, we can see that
σ2CP++(n, g) ≤ nσ2THIN ( 2 √ n log n, 2 ) + nlognσ 2 THIN (√ n log n, √ log4 n ) 2One can construct an α-thinning algorithm given a halving algorithm by recursively calling the halving
algorithm on the output of the previous stage. See App. H for further discussion.
Normalizing by the size of the input, we get that the integration error of COMPRESS++ is comparable to error of just running THIN on sets of size √ n log n to get an output of size √ n and the error of halving a data set of size 2 √ n log n. Thus, the integration error is comparable to that of methods such as RecHalve (discussed in App. H). And, we note that COMPRESS++ provides a logarithmic factor improvement in error over COMPRESS.
Example 3 (KT-SPLIT meets COMPRESS++) Example 1 showed that KT-SPLIT algorithm is Γf+ ( 32σ2n,α + 64c 2 n,α, 8cn,α ) when the input size is n and the output
size is n/α with cn,α = α · √ 32 log 4 3 ‖k‖∞‖f‖k and σ 2 n,α = α 2 4 log 4 3 log(2n) · ‖f‖2k‖k‖∞. Thus, if we run COMPRESS++(KT-SPLIT, g), we get σ2CP++ ≤ 4096 log 4
3 ‖f‖ 2 k‖k‖∞ [ n4−g−1(log4 n− g) log(2g+2 √ n) + n log(2g √ n) ]
and cCP++ ≤√ 1024 log 4
3 ‖f‖ 2 k‖k‖∞ √ n. Setting, g = log4 log4 n, gives us σ 2 CP++ which can be seen to
be factor of two from the bound on the sub-gamma parameters for KT-SPLIT with α = √ n.
Optimal g: In this setting, we derive more precise bounds in App. J for optimal g (to minimize the error) and note there that large sample values of n = 107, a default choice of gopt = 4 suffices (which we also use in all our experiments). We highlight that with such a choice of g, the run-time is order On log2 n.
Next, we give MMD guarantees for COMPRESS++ in terms of the MMD guarantees of the thinning algorithm THIN. The proof can be found in App. K.
Theorem 4 (COMPRESS++ MMD Guarantees) Let THIN(α) is an α-thinning algorithm such that on input Sin of size n outputs STHIN(α) such that
Pr [ MMD ( Sin,STHIN(α) ) ≥ an,α + √ 2vn,αt+ cn,αt ] ≤ e−t and EψTHIN(2) = 0.
Then, the √ n sized output SCP++ of COMPRESS++(THIN, g) satisfies
Pr [ MMD(Sin,SCP++) ≥ as + √ 2(σ2M,CP + vs)t+ (cs + 2RM,CP)t ] ≤ 2(n+ 2)e−t,
where s , (2g √ n, 2g), and (RM,CP, σM,CP) were defined in Theorem 2.
Example 4 (Kernel thinning meets COMPRESS++.) We now illustrate the advantages of COMPRESS++, which we call KT-Compress++, defined as COMPRESS++ with kernel thinning Dwivedi & Mackey (2021b, Alg 1.) as the thinning algorithm. Let SKT−CP and SKT−CP++ denote the output of the two stages of COMPRESS++.
Recall from Example 2 that
MMD(Sin,SKT−CP) -k √ log n (√ log n · log(nδ ) + log( n δ ) ) 1 2g √ n √ Mk(B2(R), 12g√n ),
with probability at least 1−δ. Also, recall that for thinning a set of size 2g √ n to √ n, kernel thinning has the guarantee
MMD(Sin,SKT) -k 1√n + 1√ n
√ log(2g
√ n) · Mk(B2(R), 1√n )
√ log(1/δ) + 1√
n
√ log 2g √ n log(1/δ),
with probability at least 1 − δ, whereMk(A, ε) and R carry the same meaning as in Example 2. Thus, with probability at least 1− 2δ, we have
MMD(Sin,SKT−CP++) ≤ MMD(Sin,SKT−CP) + MMD(SKT−CP++,SKT−CP)
Setting 2g = log n, we get a bound of MMD(Sin,SKT−CP++) -k 2 (√ log(nδ ) + log( n δ ) ) 1√ n √ Mk(B2(R), 1√n logn ),
On the other hand, with high probability the direct KT output of size √ n satisfies
1√ n √ log n · Mk(B2(R), 1√n ) √ log(1/δ) where we note that this bound follows from (6) by substituting m such that n/2m = 2g √ n (the output size). Comparing the bounds, we conclude that KT-COMPRESS provides a constant worse MMD error than KT itself for k, and Sin, when Mk(B2(R), 1√n ) = Θ(Mk(B2(R), 1√ n logn
)). This condition is satisfied whenMk(B2(R), ε) -d Rd(log(1/ε))ω for suitable ω; and Dwivedi & Mackey (2021a, Prop. 3) shows that this bound on covering number is valid for analytic kernels like Gaussian and inverse multiquadrics.

heading: 4 EXPERIMENTS
text: We now provide experiments (with additional experiments and suppmentary details in App. L) that demonstrate the speed ups with minimal loss of error in several experiments. We compare KT and kernel herding with their COMPRESS and COMPRESS++ versions, and in all tasks we do root thinning, namely take n input points and compress it down to √ n points. For COMPRESS++ we set g = 4 in all experiments. We evaulate the procedures both in terms of MMD error (averaged across 10 runs), and runtime (averaged across 3 runs); the error bars are too small to be visible in almost all settings. All timing experiments were run on a single core of an Intel Xeon CPU with 32GB RAM.
Target distributions and experiment settings We consider three sets of target distributions: (a) Gaussian target P = N (0, Id) for d ∈ {2, 4, 10, 100}, (b) M -Mixture of Gaussian target P = 1 M ∑M j=1N (µj , I2) with M ∈ {4, 6, 8, 31} with component means µj ∈ R2 defined in App. L, and (c) MCMC target for 12 different Bayesian posterior inference tasks with data taken from Riabiz et al. (2020a) and additional details provided in App. L. All our experiments are run with Gaussian kernel k(x, y) = exp(− 12σ2 ‖x− y‖ 2 2) with σ
2 set to 2d for (a), and (b), and chosen using median heuristic over the entire MCMC chain Garreau et al. (2017) for (c). For error performance, we evaluate the different thinning algorithms with respect to the MMD metric MMDk for (a) and (b) since it can be computed in closed-form, and for (c) since the posterior is intractable, we measure the quantity we bound namely MMDk(Pin,Pout)—the error between the input and output empirical distributions. For (a) and (b), the input sequence is generated i.i.d. from the target P, and for (c) we first standard thin the chain to get n input points which is then fed to the different thinning algorithms (see App. L for details). For all error plots, we also provide the OLS fit on the log-log scale for the error and the input size, and its slope is shown as the emprical decay rate, e.g., a slope of −0.48 for the best fit is labeled as n−0.48 for the corresponding algorithm.
Kernel thinning experiments We evaluate standard thinning (ST), kernel thinning (KT), KTCOMPRESS (denoted as KT-Comp), KT-COMPRESS++ (denoted as KT-Comp++) for all the settings (a), (b) and (c), and we present the results in Figs. 1, 2, and 6 respectively. We note that standard thinning provides close to the Monte Carlo rate of n− 1 4 with n 1 2 points in almost all the experiments,
while KT and its COMPRESS++ variants provide the best performance often achieving the n− 1 2 rate as guaranteed by our theory in moderate dimensions. The runtimes for the procedures are compared in Fig. 1 where we note that COMPRESS++ provides a noticeable 10-25x speed up over the KT runtime (and compress even more), consistent with our theoretical guarantees. In fact, the faster runtime of COMPRESS and COMPRESS++ allowed us to run experiments with them to a significantly larger input size for which we KT took many days to finish a single run (and thereby was not run).
Kernel herding experiments As another demonstration of the practical advantages of COMPRESS and COMPRESS++, we apply these procedures with kernel herding Chen et al. (2010) for Gaussian and MoG targets and report the results in Figs. 3 and 7. The results are labeled as kernel herding (Herd), COMPRESS herding (Herd-Comp), and COMPRESS++ with herding (Herd-Comp++). In all cases, we observe that Herd-Comp++ provides a competitive performance to herding, sometimes improving upon it. Fig. 3 shows COMPRESS++ with herding is 10 − 60x faster than herding for large samples sizes thereby providing further evidence about the compatibility of our meta procedures with generic thinning procedures, and their practical impact on saving computation time with minimal loss in accuracy.

heading: 5 CONCLUSIONS AND FUTURE DIRECTIONS
text: In this paper, we introduce an easy to implement and versatile meta-algorithms to speed up distribution compression tasks at almost no cost to the error guarantees. We have demonstrated, both theoretically and empirically that the performance on benchmark task are comparable to that of previous work. Furthermore, we have shown that these meta-algorithms even work well in practice with algorithms that don’t have theoretical guarantees but work well in practice e.g. herding. Analyzing the behavior of these meta-algorithms in other compression setting is an interesting direction for future work.
REPRODUCIBILITY STATEMENT
Our Python code reconstructing all experiments can be found in a supplementary zip file.

heading: 1 Introduction 1

heading: 2 Compress 3
text: 2.1 Sub-gamma property and runtime guarantees for COMPRESS . . . . . . . . . . . . . . . . . . 3
2.2 MMD guarantees for COMPRESS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

heading: 3 Compress ++ 6
text: 3.1 Sub-gamma property and runtime guarantees for COMPRESS++ . . . . . . . . . . . . . . . . 6

heading: 4 Experiments 8

heading: 5 Conclusions and Future Directions 9

heading: A Properties of Sub-Gamma Random Variables 12

heading: B Proof of Theorem 1: Sub-gamma property and runtime of COMPRESS 13

heading: C COMPRESS Calculations 14

heading: D Proof of Theorem 2: MMD Guarantees for COMPRESS 14

heading: E Proof of Lem. 2: Proof of Exponential Matrix Freedman Inequality 16

heading: F Proof of Lem. 3: Tail Bounds to Moments 16

heading: G Streaming Version of COMPRESS 17

heading: H Guarantees for RecHalve 18

heading: I Proof of Theorem 3: Sub-gamma property and runtime of COMPRESS++ 19

heading: J COMPRESS++ Calculations 19
text: K Proof of Theorem 4: COMPRESS++ MMD Guarantees 20

heading: L Additional experiments. and supplementary information 21

heading: A PROPERTIES OF SUB-GAMMA RANDOM VARIABLES
text: Lemma 1 (Summation and scaling preserve Γf+) Let f ∈ H and G1 ∈ Γ f + ( σ21 , c1 ) . Then, for all
β ≥ 0, we have β ·G1 ∈ Γf+ ( β2σ2, βc ) . Furthermore, if G2 ∈ Γf+ ( σ22 , c2 ) conditioned on G1, then
G1 +G2 ∈ Γf+ ( σ21 + σ 2 2 ,max(c1, c2) ) .
A.1 PROOF OF LEM. 1: SUMMATION AND SCALING PRESERVE Γf+
First, let us consider scaling. For any k > 0, since 〈f, ·〉H is linear, we have
〈f, k ·G1〉H = k · 〈f,G1〉H.
Thus, for any 0 < λ ≤ 1kc , E [ exp ( λ · 〈f, k ·G1〉H )] ≤ k 2σ21λ 2
2(1−kc1λ)
as required.
Furthermore, note that
〈f,G1 +G2〉H = 〈f,G1〉H + 〈f,G2〉H. Let 0 < λ ≤ min (
1 c1 , 1c2
) . Then, since G2 is conditionally sub-gamma given G1, we have
E [ exp ( λ · 〈f,G1 +G2〉H )] = E [ exp ( λ · 〈f,G1〉H + λ · 〈f,G2〉H )] = E [ E [ exp ( λ · 〈f,G1〉H + λ · 〈f,G2〉H ) |G1 ]]
= E [ E [ exp ( λ · 〈f,G1〉H ) exp ( λ · 〈f,G2〉H ) |G1 ]]
= E [ exp ( λ · 〈f,G1〉H ) · E [ exp ( λ · 〈f,G2〉H ) |G1 ]]
≤ exp (
λ2σ22 2(1−c2λ) ) · E [ exp ( λ · 〈f,G1〉H )] = exp ( λ2σ21
2(1−c1λ)
) · exp ( λ2σ22
2(1−c2λ) ) ≤ exp ( λ2(σ21+σ 2 2)
2(1−max{c1,c2}λ) ) as required.

heading: B PROOF OF THEOREM 1: SUB-GAMMA PROPERTY AND RUNTIME OF COMPRESS
text: Note that
ψCOMPRESS(Sin; 2g) = ∑4 j=1 ψCOMPRESS(Si; g) + √ n2−g−1ψHALVE(S̃)
where Si and S̃ are defined as in Alg. 1. Let ( σ2COMPRESS(n, g), cCOMPRESS(n, g) ) and(
σ2HALVE(n), cHALVE(n) )
respectively denote the ΓF+ parameters of ψCOMPRESS(·; g), and ψHALVE(·) on input of size n. Since ψCOMPRESS(Si; g) and ψHALVE(S̃) are independent conditioned on the input, we have, using Lem. 1,
σ2COMPRESS(n, 2 g) = 4σ2COMPRESS
( n 4 , 2g ) + n4−g−1σ2HALVE (√ n2g+1 ) and
cCOMPRESS(n, 2 g) = max ( √ n2−g−1cHALVE (√ n2g+1 ) , cCOMPRESS ( n 4 , 2g ))
Also, note that we can take σ2COMPRESS(4 g, 2g) = cCOMPRESS(4 g, 2g) = 0 as COMPRESS on input a point sequence of size 4g returns the same point sequence. Thus,
σ2COMPRESS(n, 2 g) = n4−g−1 log4 n−g−1∑ i=0 σ2HALVE ( 2g+1 √ n4−i ) and
cCOMPRESS(n, 2 g) = max
0≤i≤log4 n−g−1
(√ n2−g−i−1cHALVE (√ n2g+1−i )) as required.
Next, we bound the running time of COMPRESS. Let rCP(n) be the running time of the algorithm on input of size n. From the definition, COMPRESS makes four recursive calls to COMPRESS on inputs of size n/4. Then, HALVE is run on an input of size 2g+1 √ n. Thus, r satisfies the recursion
rCP(n) = 4rCP(n/4) + rHALVE( √ n2g+1).
Note that rCP(4g) = 1. Unrolling the recursion, we get
rCP(n) = n4 −g + log4 n−g−1∑ i=0 4irHALVE(2 g+1 √ n4−i)
as required.

heading: C COMPRESS CALCULATIONS
text: To get a sense for the guarantees given by the above algorithm, we will instantiate this for various settings for σ2HALVE. In particular, natural choices for σHALVE are constant and logarithmic growth.
• σ2HALVE(n) = C log4 n
σ2CP(n, g) = C n 4g+1 · ∑log4 n−g−1 i=0 log4 ( 2g+1 √ n4−i ) = C n4g+1 · 1 2 ∑log4 n−g−1 i=0 (log4 n+ g + 1− i)
= C n4g+1 · 1 2 ∑log4 n+g+1 j=2g+2 j
(i) = C n4g+1 · 1 4 (log4 n+ 3g + 3)(log4 n− g),
• σ2HALVE(n) = C σ2CP(n, g) = Cn4 −g−1(log4 n− g).

heading: D PROOF OF THEOREM 2: MMD GUARANTEES FOR COMPRESS
text: By the definition of the COMPRESS algorithm, the COMPRESS discrepancy vector (2) satisfies
ψCP(S, g) = √ n2−g−1ψHALVE(S̃, g) + ∑4 i=1 ψCP(Si, g).
Unrolling the recursion, we find that HALVE is applied to 4k coresets of size √ n2g+12−k. If we
denote these coresets by S̃k,1, · · · , S̃k,4k , we have
ψCP = √ n2−g−1 ∑log4 n−g−1 k=0 ∑4k j=1 2
−kψHALVE(S̃k,j , g). (8) For each (k, j), let
ψ̃k,j = ψHALVE(S̃k,j , g) represent the HALVE discrepancy vector. Then, for each l, the variables (ψ̃l,j)j≥1 are conditionally independent and zero mean given (ψ̃k,j)k>l,j≥1.
Moreover, if we number the elements of Sin as (x1, . . . , xn), ‖ψ̃k,j‖2k = u>ψ̃k,jKuψ̃k,j = ∥∥∥K1/2uψ̃k,j∥∥∥22 for
K = (k(xi, xj))1≤i,j≥n and uψ̃k,j = ∑n i=1 ei(1− n nout I[xi ∈ Sout]) ∈ Rn.
For each u ∈ Rn, consider the symmetric dilation Mu ∈ Rn+1×n+1 given by( 0 uT
u 0
) .
We will abbreviate Muψ̃k,j as Mk,j . We will bound the spectral norm of this matrix in order to control ‖ψ̃k,j‖k. In order to do this, we will use the following subexponential version of the Matrix Freedman inequality of (Tropp, 2011, Thm. 1.2), a proof of which can be found in App. E.
Lemma 2 (Exponential Matrix Freedman Inequality) Consider a finite sequence of random matrices Xk ∈ Rd×d on that are self-adjoint. Let
E [ Xk|{Xi}i<k ] = 0
E [ Xpk |{Xi}i<k ] p!2 R p−2A2k.
Then, if σ2 = ∥∥∑A2k∥∥ Then, for all t > 0
Pr [ λmax (∑ kXk ) ≥ t ] ≤ d exp ( −t2 2(σ2+Rt) ) .
First note that ‖ψCP‖k = n ·MMD(Sin,SCP). So it suffice to bound ‖ψCP‖k. We will do this using Lem. 2. Note that for all the summands in (8) ,EMk,j = 0 and furthermore,
M2qk,j = ‖ψ̃k,j‖ 2q−2 k M 2 k,j ‖ψ̃k,j‖ 2q k I
M2q+1k,j = ‖ψ̃k,j‖ 2q−1 k Mk,j ‖ψ̃k,j‖ 2q+1 k I
Therefore, EMqk,j E‖ψ̃k,j‖ q kI . By assumption, we have that HALVE satisfies Pr [ MMD(Sin,SHALVE) ≥ an + √ 2vnt+ cnt ] ≤ e−t.
on input sequences of size n. This implies that Pr [ ‖ψ̃ell,j‖k ≥ a′` + √ 2v′`t+ c ′ `t ] ≤ e−t for a′` = `a`, v ′ ` = ` 2v`, c ′ ` = `c` for ψ̃k,j of size ` = 2
−k+g+1√n. We will use the following result that translates tail bounds into moment inequalities. A proof can be found in App. E.
Lemma 3 (Tail Bounds to Moments) Let Z be a non-negative random variable such that Pr [ Z > a+ √ 2vt+ ct ] ≤ e−t.
Then, EZq ≤ Γ ( q 2 + 1 ) (8v + 2a2)q/2 + (q)!(4c)q.
For ψ̃k,j , this implies E‖ψ̃k,j‖qk ≤ Γ(q/2 + 1)(8v′` + 2a ′2 ` ) q/2 + (q)!(4c′`) q for ` = 2g+1−k
√ n.
Thus, E‖wq ψ̃k,j Mqk,j‖k q! 2 R q−2 ` A 2 ` for R` = wψ̃k,j ( √ 8v′` + 2a ′2 ` + 4c ′ `) and A` =
2wψ̃k,j ( √ 8v′` + 2a ′2 ` + 4c ′ `)I where wψ̃k,j is the coefficient of ψ̃k,j in (8).
Thus, from Lem. 2, we have ψCP has exponential tails i.e. Pr [ ‖ψCP‖k ≥ t ] ≤ 2(n+ 1) exp ( −t2
2(σ2+Rt) ) with
σ2 = ∑log4 n−g−1 k=0 ∑4k j=1 n2
−2k−2g−2 (√
8v′√ n2g+1−k + 2a ′2√ n2g+1−k + 4c′√ n2g+1−k )2 = ∑log4 n−g−1 k=0 n2 −2g−2 (√ 8v′√ n2g+1−k + 2a ′2√ n2g+1−k + 4c′√ n2g+1−k
)2 ≤ (√
8v′√ n2g+1 + 2a ′2√ n2g+1 + 4c′√ n2g+1 )2∑log4 n−g−1 k=0 n2 −2g−2
≤ (√
8v′√ n2g+1 + 2a ′2√ n2g+1 + 4c′√ n2g+1 )2 n2−2g−2(log4 n− g)
and
R = max k≤log4 n−g−1
√ n2−g−1−k (√ 8v′√
n2g+1−k + 2a ′2√ n2g+1−k + 4c′√n2g+1−k ) ≤ √ n2−g−1 (√ 8v′√
n2g+1 + 2a ′2√ n2g+1 + 4c′√n2g+1
) .
We get the required bound by inverting expression in the exponential.

heading: E PROOF OF LEM. 2: PROOF OF EXPONENTIAL MATRIX FREEDMAN INEQUALITY
text: Lemma 4 Consider a finite sequence of random matrices Xk ∈ Rd×d that are self-adjoint. Let E [ Xk|{Xi}i<k ] = 0
E [ Xpk |{Xi}i<k ] p!2 R p−2A2k.
Then, if
σ2 = ∥∥∑A2k∥∥
Then, for all t > 0
Pr [ λmax (∑ kXk ) ≥ t ] ≤ d exp ( −t2 2(σ2+Rt) ) .
We first note the following lemmas about the moment generating functions for matrix valued random variables.
Lemma 5 (Master Tail Bounds (Tropp, 2011) ) Let Xk be an adaptive sequence of matrix valued random variables in Rd (with respect to the filtration Fk and let Vk be a sequence that is measureable with respect to Fk−1. Assume that
logE [ exp(θXk)|Fk−1 ] g(θ)Vk
Let Yk = ∑ i≤kXi and Wk = ∑ i≤k Vi. Then,
Pr [ ∃k > 0 : λmax(Yk) ≥ t and λmax(Wk) ≤ w ] ≤ d exp ( −θt+ g(θ)w ) .
Lemma 6 (Matrix MGF Exponential (Tropp, 2012)) LetX be a random self-adjoint matrix such that
EX = 0 EXq p!2 A 2 for q = 2, 3, . . .
Then,
E exp(θX) exp ( θ2 2(1−θ) ·A 2 ) for 0 ≤ θ ≤ 1.
Assume that R = 1. We get the general result by scaling. We apply Lem. 6 conditional on {Xi}i<k along with the operator monotonicity of log to get
logE [ exp(θXk)|{Xi}i<k ] θ 2 2(1−θ) ·A 2
By assumption, we have that λmax (∑ iA 2 i ) ≤ σ2. Thus, applying Lem. 5, we get
Pr [ λmax (∑ iXi )] ≤ d exp ( −θt+ θ 2σ2 2(1−θ) ) We get the result by choosing θ = t/(t+ σ2).

heading: F PROOF OF LEM. 3: TAIL BOUNDS TO MOMENTS
text: We bound the moments of the random variable X = Z − a. First note that X = X+ −X− where X± = max(±X, 0). Also, |X|q = Xq+ +X q −. Furthermore, X q − ≤ aq by assumption.
E [ |X|q ] ≤ aq + E [ Xq+ ] .
In order to bound the moments of X+, we use the tail bounds on X . Note that Pr[X+ > t] = Pr[X > t] for any t > 0.
E [ Xq+ ] = q ∫∞ 0 tq−1 Pr[X > t]dt
≤ q ∫∞ 0 ( √ 2vt+ ct)q−1 (√ 2vt+ct 2t ) e−t t dt
≤ q2 ∫∞ 0 ( √ 2vt+ ct)q e −t t dt
≤ q2q−2 ∫∞ 0 ( (2tv)q/2 + (2ct)q ) e−tdt
≤ 2q−2 ( Γ(q/2 + 1)2q/2vq/2 + (q)!2 (2c) q )
≤ Γ(q/2 + 1)(8v)q/2 + (q)!(4c)q.
Thus, we have EZq ≤ 2q−1(aq + E|X|q) which gives us the desired result.

heading: G STREAMING VERSION OF COMPRESS
text: COMPRESS can be efficiently implemented in a streaming fashion, by unrolling the recursion in the original COMPRESS and looking at the algorithm from a bottom up point of view. In this setup, the algorithm maintains memory in layers, and takes in inputs in batches of size 4g+1. Whenever there is more than 2i4g+1 points in layer i, the algorithm uses HALVE on the points in the layer and moves the points up to level i + 1. The points in the highest non-empty layer form the current approximation to the points seen so far. See Alg. 3 for formal description.
Algorithm 3: COMPRESS (Streaming) Input: Halving algorithm HALVE, oversampling factor g, input stream x1, x2, . . . S0 ← {} for t = 1, 2, . . . , do S0 ← S0 ∪ (xj)t·4 g+1 j=1+(t−1)·4g+1 // Process input in batches of size 4 g+1
if t == 4j then Sj+1 ← {} end for i = 0, . . . , log4 t+ 1 do
if |Si| == 2i4g+1 then S = HALVE (k,Si) // coreset of size 2i−14g+1 Si+1 ← Si+1 ∪ S Si ← {}
end end return Sblog4 tc+1 // Coreset stored in the highest level is the desired approximation at current stage
end
Next, we analyze the space complexity of the streaming version of the algorithm. The intuition for gains in storage is very similar to that for running time: we only maintain and run HALVE only on subsets of much smaller size than the original input. See App. G for the proof.
Proposition 1 (COMPRESS Streaming Memory Bound) Let Halve be an algorithm that that takes as input a point sequence Sin of size n and return a data set S with size n/2 with space complexity sHALVE(n). Then, at any time t, streaming version of COMPRESS with oversampling factor g has space complexity
sCP(t) = 4 g+3 √ td+ s(2g+1 √ t)d
real words of memory.
G.1 PROOF OF PROP. 1: COMPRESS STREAMING MEMORY BOUND
For convenience assume that n is divisible by 4g+1. At time t, we would like to estimate the space usage of the algorithm. At the ith level of memory, we can have at most 2i+24g − 1 data points.
Since we are maintaining a data set of size at most √ t4g at time t, there are at most log t2 levels. Thus, the maximum memory used up to time t is bounded by
0.5 log t∑ i=0 2i+24g ≤ 4g+3 √ t
data points stored. Furthermore, at any time up to time t, we have run HALVE on a point sequence of size at most √ t2g+1 which requires storing at most s (√ t2g+1 ) additional points. Each data point
is stored with d real words of memory giving the required result.
H GUARANTEES FOR RECHALVE
Algorithm 4: Recursive Halving (RecHalve) Input: Input point sequence Sin of size n, Halving algorithm HALVE, Thinning factor m if m = 0 then return Sin ; else S1 ← HALVE(Sin) // halved coreset size n/2 return RecHalve (S1,HALVE,m− 1) end
Proposition 2 (RecHalve Subgamma Guarantees) Let HALVE be a halving algorithm that is ΓF+ ( σ2HALVE, cHALVE ) . Then, RecHalve with depth parameter m is a 2−m √ n-root thinning algorithm and is ΓF+ with parameters
σ2RH(n,m) = ∑m−1 i=0 4 iσ2HALVE ( 2−in ) , and
cRH(n,m) = max0≤i≤m−1 ( 2i · cHALVE ( n2−i )) .
Denote by S1 the output of the halving algorithm after the first round. Then,
ψRH(Sin;m) = ψHALVE(Sin) + 2ψRH(S1;m− 1).
From the ΓF+ guarantees of HALVE and noting the independence conditioned on the inputs, we have
σ2RH(n,m) = σ 2 HALVE(n) + 4σ 2 RH(n/2,m− 1)
cRH(n,m) = max ( cHALVE ( n), 2 · cRH ( n/2,m− 1 ))) Also note that σ2RH(n, 0) = cRH(n, 0) = 0. Unrolling the recursion, we get
σ2RH(n,m) = ∑m−1 i=0 4 iσ2HALVE ( 2−in ) cRH(n,m) = max0≤i≤m−1 ( 2i · cHALVE ( n2−i
)) as required.
Proposition 3 (RecHalve Running Time) Let HALVE be a halving algorithm with running time h(n) on inputs of size n. Then, RecHalve with depth parameter m is a 2−m √ n-root algorithm with running time ∑m−1 i=0 h(n2 −i)
on inputs of size n.

heading: I PROOF OF THEOREM 3: SUB-GAMMA PROPERTY AND RUNTIME OF COMPRESS++
text: As in the definition of the algorithm, let SCOMPRESS denote the output of COMPRESS. Recall that from Proposition 1, we have that
ψCOMPRESS(Sin; g) , ∑ x∈Sin k(x, ·)− √ n 2g ∑ x∈SCOMPRESS k(x, ·)
is Γf+-sub-Gaussian with parameter
σCOMPRESS(n, g) = √ n4−g−1 ∑log4 n−g i=0 σ 2 THIN ( 2g+1−i √ n, 1/2 ) cCOMPRESS = max0≤i≤log4 n−g−1 (√ n2−g−i−1cTHIN (√ n2g+1−i, 1/2 )) .
Furthermore, by assumption, we have that ψTHIN(SCOMPRESS) , ∑ x∈SCOMPRESS k(x, ·)− 2 g ∑ x∈SCOMPRESS++ k(x, ·)
is Γf+ with parameters
σ2THIN ( 2g √ n, 2−g ) and
cTHIN ( 2g √ n, 2−g ) conditioned on the input SCOMPRESS. Also, note that
ψCOMPRESS++ = ψCOMPRESS + √ n 2g ψTHIN(SCOMPRESS)
Thus, we have that ψCOMPRESS++ is Γ f + with parameters
σ2THIN(n, g) = n4 −g−1∑log4 n−g i=0 σ 2 THIN ( 2g+1−i √ n, 1/2 ) + n4−gσ2THIN ( 2g √ n, 2−g ) and
cCP++(n, g) = max ( √ n2−gcTHIN ( 2g √ n, 2−g ) ,max0≤i≤log4 n−g−1 (√ n2−g−i−1cCP++ (√ n2g+1−i, 1/2 ))) as required.
The runnning time bound follows from the running time of COMPRESS ( THIN(1/2), g ) and
THIN ( 2g √ n, 2−g ) .

heading: J COMPRESS++ CALCULATIONS
text: The squared sub-Gaussian constant of KT with input size √ n · 2g and output size √ n is
σ2KT( √ n2g; √ n) = C · 4 g−1 3 log4( √ n2g)
where C is the constant floating with kernel-halving. The run-time is C ′n4g. Thus if we run COMPRESS (g) with kernel halving, and then KT on the output of size √ n · 2g, the overall subGaussian constant for
ψCg-KT(Sin) , ψCP(Sin; g) + √ n 2g ψKT(SCP)
is given by
C n4g+1 · 1 4 (log4 n+ 3g + 3)(log4 n− g) + n 4g · C · 4g−1 3 · log4(
√ n2g)
= C n4g+1 · 1 4 (log4 n+ 3g + 3)(log4 n− g) + n 4g · C · 4g−1 6 · (log4 n+ g).
Define α , log4 n and a , 4 g, and note that log4 x = c log x for c , 1/ log 4. Then the expression simplifies to
n 2 [ 1 8 · 1 4g ( α2 − αg + 3αg− 3g2 + 3α− 3g ) + 13 · 1 4g ( (4g − 1)α+ (4g − 1)g )] = n2 [ 1 8 · 1 4g ( α2 + 3α+ g(2α− 3)− 3g2 ) + 13 · 1 4g ( (4g − 1)α+ (4g − 1)g
)] = n2a [ 1 8 ( α2 + 3α+ c log a(2α− 3)− 3c2(log a)2 ) + 13 (a− 1)(α+ c log a) ] .
Note that dda [g(a)/a] = 0⇔ g(a) = ag ′(a). Taking the derivative and setting it to 0, we have
3α2 + 9α+ 3c(2α− 3) log a− 9c2(log a)2 + 8a(α+ c log a)− 8(α+ c log a) = 3c(2α− 3)− 18c2 log a+ 8a(α+ c log a) + 8ac− 8c
which implies
8ac = 3α2 + 9α+ 3c(2α− 3) log a− 9c2(log a)2 − 8(α+ c log a)− 3c(2α− 3) + 18c2 log a+ 8c ⇒ a = 18c [ 3α2 + α+ 17c− 6cα+ (6cα− 17c+ 18c2) log a− 9c2(log a)2 ] = 17−6α8 + 3α2+α 8c + (6α−17+18c) 8 log a− 9c 8 (log a) 2. (9) Let gopt denote the optimal g obtained by solving (9) (recall 4g = a). In Fig. 4, we plot the values computed numerically using fixed point iterations for (9). We note that even for large sample values of n = 107, a default choice of gopt = 3 suffices.
Also, note that from the above expression We highlight that with such a choice of g, the run-time (using Theorem 4) is bounded as n 4g + C(log4 n− g)4 g+1n+ Cn4g ≤ nlog4 n + Cn(4(log4 n) 2 + log4 n− 4 log4 n log4 log4 n)
≤ 5Cn(log4 n)2 with g = gopt.
K PROOF OF THEOREM 4: COMPRESS++ MMD GUARANTEES
Recall that
ψCP++ = ψCP + √ n 2g ψTHIN(SCP)
By assumption, we have that on input Sin of size n, THIN outputs Sout such that Pr [ MMD(Sin,Sout) > a`,α + √ 2v`,αt+ c`,αt ] ≤ e−t.
Again, recall that nMMD = ‖ψ‖k. Thus, we need to bound 1 n‖ψCP++‖k = 1 n‖ψCP + √ n
2g ψTHIN(SCP)‖k ≤ 1n‖ψCP‖k + 1√ n2g ‖ψTHIN‖k
= MMD(Sin,SCP) + MMD(SCP,Sout).
From Theorem 2, we have that
Pr [ MMD(Sin,SCP) > 2γRMMD,CP + √ 2γσ2MMD,CP ] ≤ 2(n+ 1) exp(−γ)
for σ2MMD,CP ≤ 2−2g−2 (√ 8v√n2g+1 + 2a 2√ n2g+1 + 4c√n2g+1 )2 (log4 n− g),
and RMMD,CP ≤ 2−g−1 (√
8v√n2g+1 + 2a 2√ n2g+1 + 4c√n2g+1
) .
Since, the guarantees of THINare conditioned on the inputs, we get
Pr [ MMD(Sin,SCP++) > 2tRMMD,CP + √ 2tσ2MMD,CP + a2g √ n,2g + √ 2v2g √ n,2gt+ c2g √ n,2−gt ] ≤
Pr [ MMD(Sin,SCP) > 2tRMMD,CP + √ 2tσ2MMD,CP ] + Pr [ MMD(Sin,Sout) > a2g√n,2g + √ 2v2g √ n,2gt+ c2g √ n,2−gt ] ≤ ( e−t ) + ( 2(n+ 1)e−t
) ≤ 2(n+ 2)e−t.
Thus, we get
Pr [ MMD(Sin,SCP++) ≥ a2g√n,2g + √ 2(σ2MMD,CP + v2g √ n,2g)t+ (c2g √ n,2g + 2RMMD,CP)t ] ≤ 2(n+ 2)e−t
as required.

heading: L ADDITIONAL EXPERIMENTS. AND SUPPLEMENTARY INFORMATION
text: MoG Target Our mixture of Gaussians target is given by P = 1M ∑M j=1N (µj , Id) for M ∈ {4, 6, 8} where
µ1 = [−3, 3]>, µ2 = [−3, 3]>, µ3 = [−3,−3]>, µ4 = [3,−3]>, µ5 = [0, 6] >, µ6 = [−6, 0]>, µ7 = [6, 0]>, µ8 = [0,−6]>.
Visualizing coresets In Fig. 5, we visualize the coresets for a 32-MoG target in with means located on two concntric circles, from where we observe the visible advantages of all the thinning strategies over standard thinning, and the comparable performance of the COMPRESS++ procedures. We see that even with 32 points all procedures succeed in locating the means for the MoG, thereby providing a better approximation to the underlying target with i.i.d. points.
In Figs. 6 and 7, we present the MMD errors of KT and herding experiments for the MOG target, and notice again that COMPRESS++ provides a competitive performance to the original algorithm.
MCMC Target Our set-up for MCMC experiments follows closely the set-up of Dwivedi & Mackey (2021b). For complete details on the targets and sampling algorithms we refer the reader to Riabiz et al. (2020a, Sec. 4).
Goodwin and Lotka-Volterra experiments The dataset of Riabiz et al. (2020b) consists of output of four distinct MCMC procedures targeting each of two d = 4-dimensional posterior distributions P: (1) a posterior over the parameters of the Goodwin model of oscillatory enzymatic control (Goodwin, 1965), and (2) a posterior for the Lotka-Volterra model of oscillatory predator-prey evolution (Lotka, 1925; Volterra, 1926)—both of these are differential equation based models and often used as benchmark datasets for Bayesian infere.ce Riabiz et al. (2020b) provide 2×106 sample points for each task from four MCMC chains: Gaussian random walk (RW), adaptive Gaussian random walk (adaRW, Haario et al., 1999), Metropolis-adjusted Langevin algorithm (MALA, Roberts & Tweedie, 1996), and pre-conditioned MALA (pMALA, Girolami & Calderhead, 2011).
Hinch experiments Riabiz et al. (2020b) also made available the output of two independent Gaussian random walk MCMC chains targeting each of two d = 38-dimensional posterior distributions P for a cardiac experiment that approximates a posterior over the parameters of the Hinch model of calcium signalling in cardiac cells (Hinch et al., 2004) and (2) a tempered version of the same posterior, as defined by Riabiz et al. (2020a, App. S5.4).
Burn-in times For all chains, the initial burn-in points of each chain using the maximum burn-in period reported in Riabiz et al. (2020a, Tabs. S4 & S6, App. S5.4). For each Hinch chain, we further normalized by subtracting the post-burn-in sample mean and dividing each coordinate by its post-burn-in sample standard deviation. To obtain an input sequence Sin of length n to be fed into a thinning algorithm, we downsample the remaining even indices of points using standard thinning (odd indices are held out). When applying standard thinning to any Markov chain output, we adopt the convention of keeping the final sample point.
The selected burn-in periods for the Goodwin task were 820,000 for RW; 824,000 for adaRW; 1,615,000 for MALA; and 1,475,000 for pMALA. The respective numbers for the Lotka-Volterra task were 1,512,000 for RW; 1,797,000 for adaRW; 1,573,000 for MALA; and 1,251,000 for pMALA.
MCMC Results From Fig. 8, we note that in some cases standard thinning itself provides a better than Monte Carlo approximation to the input, but KT and its COMPRESS++ variant are able to provide further improvement. For the Hinch experiments, we found that the symmetrization step as noted in Rem. 4 (and enforced in all the KT/Herding experiments) lead to worse performance, and not implementing it enabled COMPRESS++ (labeled as KT-Comp++(NS)) to once again provide competitive performance to KT.

Review summaries:
This paper gives a meta algorithm for speeding up coreset constructing algorithms for the distribution compression problem.  The benefit of this meta algorithm is that its running time is faster by a square-root factor (e.g. quadratic to linear) while keeping the error rate roughly the same: only a factor of 4 worse.  The method is very simple: split the input into four pieces of size n/4 each, recursively solve each of them, then combine all four answers into a set and take a coreset of this set.  The speed-up in running time is immediate (just because of how recursive formulas work) and the error bounds are also rather clear.

The topic of the study is distribution-compression/thinning algorithms for Monte-Carlo estimation of functions in RKHS, i.e., given a set of $n$ points, such that uniform distribution on these n points approximates an underlying distribution (with respect to integration over functions in a RKHS) with a certain error, output a set of $m$ points $(m \ll n)$ that has a similar error. The paper proposes three reductions/meta-algorithms to speed up existing distribution-compression algorithms while maintaining the error guarantee of the original algorithms. The main result is that the proposed reduction improves the runtime of existing algorithms (with quadratic runtime or more) by a quadratic factor without increasing the corresponding error by a polylog factor.

The paper introduces meta algorithms "Compress" and "Compress++" that take existing thinning algorithm as a subroutine and improves on their runtime while incurring marginally more error.   The algorithm "Compress" runs in recursive fashion: divides input set into four parts, runs Compress on each part independently, combines the resulting sets and halves the combined set using "HALVE" algorithm. As an example, the authors demonstrate that one can use KT-SPLIT as HALVE algorithm and get faster runtimes to compress albeit suffering extra error (upto log factors). Compress ensures that KT-SPLIT runs on sets of much smaller size.   They build on Compress to obtain better error rates while sacrificing on runtimes (upto log factors). For this, they use the thinning algorithm on much smaller set obtained after running Compress. The runtimes are quite easy to derive based on recursion of algorithms. They use sub-gamma property to prove the error guarantees of the algorithms. They demonstrate the faster runtimes of their algorithms using high dimensional Monte Carlo samples.