Fast Bayesian Coresets via Subsampling and Quasi-Newton Refinement
Authors: Cian Vasanttilak Naik
Judith Rousseau
Trevor Campbell
Venue: NeurIPS 2022 Accept (2022)
Keywords: Bayesian coresets
Markov chain Monte Carlo
data subsampling

Abstract: Bayesian coresets approximate a posterior distribution by building a small weighted subset of the data points. Any inference procedure that is too computationally expensive to be run on the full posterior can instead be run inexpensively on the coreset, with results that approximate those on the full data. However, current approaches are limited by either a significant run-time or the need for the user to specify a low-cost approximation to the full posterior. We propose a Bayesian coreset construction algorithm that first selects a uniformly random subset of data, and then optimizes the weights using a novel quasi-Newton method. Our algorithm is a simple to implement, black-box method, that does not require the user to specify a low-cost posterior approximation. It is the first to come with a general high-probability bound on the KL divergence of the output coreset posterior. Experiments demonstrate that our method provides significant improvements in coreset quality against alternatives with comparable construction times, with far less storage cost and user input required. 

Fast Bayesian Coresets via Subsampling and Quasi-Newton Refinement
Cian Naik
Judith Rousseau
Department of Statistics
Department of Statistics
University of Oxford
University of Oxford
cian, nwi.kstats ox. ac. uk
judith. rousseadaausss ox ac. uk
Trevor Camphell
Department of Statisties University of British Columbia trevor@stat. .ubistttt.bb..... ca
Abstract
Bayesian coresets approximate posterior distribution by building small weighted subset of the data points. Any inference procedure that is too computationally expensive to be run on the full posterior can instead be run inexpensively on the coreset, with results that approximate those on the full data. However, current approaches limited by either significant run-time the need for the user specify a low-cost approximation to the full posterior. We propose a Bayesian coreset construction algorithm that first selects a uniformly random subset of data, and then optimizes the weights using novel quasi-Newton method. Our algorithm simple to implement, black-box method, that does not require the user to specify a low-cost posterior approximation. It is the first to come with a general high-probability bound on the KI. divergence of the output coreset posterior. Experiments demonstrate that our method provides significant improvements in coresct quality against alternatives with comparable construction times, with far less storage cost and user input required.
Introduction
Bayesian methods are key tools for parameter estimation and uncertainty quantification, but exact inference is rarely possible for complex models, Currently, the gold standard method for approximate inference is Markov chain Monte Carlo (MCMC) [1; 2: 3, Ch. 11,12], which involves simulating Markov chain whose stationary distribution the Bayesian posterior. However, modern applications often concerned with very large datasets; in this setting, MCMC methods typically have ((NT) complexity- for samples and dataset size N wwhic quickly becomes intractable as N increases. Motivated by stochastic methods in variational inference [4, 5], this cost can be reduced by involving only random subsample M < data points in each Markov chain iteration [6-13]. This reduces the per-iteration computation time of MCMC, but distribution of the resulting Markov chain and cause slow mixing [11, 14-17j.
create substantial error in the stationary
Bayesian coresets (18-23) provide an alternative for reducing the cost of MCMC and other inference methods. The key idea is that in a large-scale data setting, much of the data is often redundant. In particular, there often exists fixed small, weighted subset of the data-a coreset- -tat suffices to capture the dataset as whole in some sense [21]. Thus, one can find such a coreset, it can be used in place of the full dataset in the MCMC algorithm, providing the simplicity, generality, and per-iteration speed of subsammee MCMC without the statistical drawbacks. Current state-of-
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

the-art Bayesian coreset approaches based on sparse variational inference [21] empirically provide high-quality coresets, but are limited by their significant run- run-time and lack of theoretical convergence guarantees. Methods based on sparse regression [19, 20, 23] are significantly faster and come with some limited theoretical guarantees, but require low-cost approximation to the full posterior. This approximation is often impractical to find and difficult to tune. Moreover, it can fundamentally limit the quality of the coreset obtained. Aside from [21], these methods also require ((V) storage, which is problematic in the large-data setting. There is third class of constructions based on importance sampling [18], but these do not provide reliable posterior approximations in practice.
In this work, we develop novel coreset construction algorithm that is both faster and easier to tune than the current state of the art, and has theoretical guarantees on the quality of its output. The key insight in our work is that we can split the construction of coreset into two stages: first we select the data in the coreset via uniform subsampling an idea developed concurrently by [24, 25j-and then optimize the weights using few steps of a novel guasi-Newton method. Selecting points uniformly randomly avoids the slow inner-outer loop optimization of [21], and guarantees that the optimally- weighted coreset has low KL divergence to the posterior (Theorems 4.1 and 4.2). Weighting these uniformly selected points correctly is crucial, and our quasi-Newton method is guaranteed to converge to point close to the optimally-
coreset (Theorem 4.3) in significantly fewer iterations than the method of [21]. Finally, our method is easy to tune and does not require the low-cost posterior approximation of [19, 20]. Our experiments show that the algorithm exhibits significant improvements in coreset quality against alternatives with comparable construction times.
2 Bayesian coresets
The problem we study is as follows. Suppose we are given target probability density x(0) for variable 0 E Q that is comprised of N potentials (fn(0)) and base density (o))..
N
π (0)
Z exp Σ (0)
To (0).
(I)
where Z is the normalizing constant. This setup corresponds to a Bayesian statistical model with prior wo and i.i.d. data X, conditioned on 0, where (0) = log pl X,, (0). Computing expectations under æ exactly is often intractable. MCMC methods can be employed here, but typically have a computational complexity of ((NT) to obtain T samples, since Σn, fn (0) needs to be evaluated at each step. Bayesian coresets [18] offer an alternative, by constructing a small, weighted subset of the dataset, on which any MCMC algorithm can then be run.
To do this, we find set of weights w E R. corresponding each data point, with the constraint
that only M < N of these nonzero, i.e. ||w|| = Em= Tu M. The weighted subset of points corresponding to the strictly positive weights is called a coreset, and we can then create the coreset posterior approximation,
Σ wnfn(()) To (0)
To
(w) exp
where Z Z(w) is the new normalizing constant, and π = ग corresponds the full posterior. Running MCMC on this approximation has complexity O(MT), a considerable speedup if M << N.
Constructing a sparse set of weights w that T.. is as close to possible is the key challenge here. Methods based on importance sampling [18] and sparse regression [19, 20, 23] have been developed, but these provide poor posterior approximations or require first finding low-cost posterior approximation. The current state-of-the-art approach resolves both of the aforementioned issues by formulating the coreset construction problem as variational inference in the family of coresets [21]:
w" arg min WERN
DKL. (mwllmm) s.t. 2 0, |/ llo S M,
where (-I--ll-) denotes the KL divergence. Noting that coresets are a sparse subset of an exponential
family-the the weights form the natural parameter w E R. the component potentials (fn (0)) n=
form the sufficient statistic, log Zww) is the log partition function, and xo is the base density.
Tw (0) == exp (wT f(0) log Z(w)) To()), f(0) = (fi ... (() E R,,
(4)
2

---

one can obtain a formula for the gradient
W.. (Tw [ = - Vu log Z (1)(1 - w) = Cov (1-w)) (1,J"(1-))]
(5)
where Cov denotes covariance under Tw. Here, and throughout, we take derivatives with respect w the ambient space w R The method of [21] then involves selecting a data point add to the coreset, optimizing the coreset weights using stochastic gradient estimates based on Eq. (5), and then iterating. In practice, this method infeasibly slow outside of small-data problems. For every data point selected it requires sampling from Tы a number of times equal to the number of optimization steps (set to 100 in (21 This can easily lead to tens of thousands of sampling steps in total, each of which typically necessitates the use of an MCMC algorithm. Conversely, our method requires tens of sampling steps in total.
3 Construction via subsampling and quasi-Newton refinement
In this section, we provide new Bayesian coreset construction algorithm (Algorithm 1) to solve the sparse variational inference problem Eq. (3). This method involves first uniformly subsampling the data, and then optimizing the weights the subsample. Here, we give a detailed explanation of how the algorithm works. We defer to Section the theoretical analysis of the method, which demonstrates that, with high probability, it finds near-optimal coreset in small number of optimization iterations.
3.1 Uniform subsampling
The first key insight to our approach is that select the M data points that comprise the coreset by simple uniform subsampling. From an algorithmic standpoint, the benefit of this approach is clear: it makes the selection of coreset points fast and easy to implement. It also decouples the subset selection from the optimizing of the coreset weights, leaving us with a simpler optimization problem to solve.' 1 Intuitively, this is reasonable from a coreset quality standpoint because uniform subsample of the data will, with high probability, create functions (fm))= with which to approximate the full data log-likel Of course, we still need to optimize the coreset weights; but we have not limited the flexibility of the coreset family significantly by selecting points via uniform subsampling. Theorem 4.2 provides precise statement of these ideas, and shows that the optimal weighted coreset posterior built using a uniformly random subset of data is good posterior approximation.
"basis" of log-likeliiee ligeelkedoood
Since the coreset data points are chosen uniformly randomly, we will w.l.o.g. re-label the selected subset to have indices I,.... M, and the remaining data points to have indices M + +1,... 1,.... N from now on. We can specify a new weight constraint set W C N of vectors w with 0 entries for every index beyond M, and a further subset Wy where the weights sum to N,
W =
{ w E N w 2 0, n >> M = Wn = =0},
Wy = fwe 1s == N}.
... is also useful to denote g(0) (M())) E R-M as the first M components of (0).
3.2 Quasi-Newton optimization
Given the choice of M data points to include in the coreset, it is crucial to optimize their weights correctly. This is the second step of our algorithm. The KL divergence gradient for the M weights is
V DKL (T)) (- v², log ((w)(1 - w))
()) 1:M = - Cov [g. f (1(1w]]] w)]
(6)
To reduce the number of weight optimization iterations each of which generally involves MCMC sampling from Tw. where w is the current set of weights-we develop a second-order optimization method. Again using the properties of exponential families, we can derive the Hessian of the KL divergence with respect to the M weights
2 DEL (Tw m = (V2 log 7 (w:) V3 log (1w) (1-w) 1:M.1MM
'The benefits of this initial uniformly random selection step also noted by [24, 25] in concurrent work.
3

---

This matrix is not guaranteed to be positive (semi-)definite. However, writing out the second term explicitly as
V log (1) - w) - Er [ Er. Er. (o" -ew ET. (1 w)))
we see that this term will be small near the optimum, where the coreset approximation ideally
good, i.e., w² f(0) ~ 1² The first term does not contain this f w) term, and thus should dominate the expression.
This motivates the use of
V (Tw||m)) R (V log (w)) 1:M,1:M = Cov [9.9]
(7)
scale gradient steps in the optimization method rather than the true Hessian, thus creating quasi-Newton method [26]. This heuristic motivates our approach, but it is rigorously justified by Theorem 4.3, which proves that this scaling results in optimization iterations that are guaranteed to converge exponentially to coreset which is nearly globally optimal. Note that when M is larger than the inherent dimension of the space of log-likelihood functions, Cov w 9 9] will have zero eigenvalues and be noninvertible. Therefore we add regularization r >> 0 prior to inversion. Theorem 4.3 shows how the regularization influences the optimization, but in general we want 7 to be as small as possible while still ensuring the numerical stability of inverting Cov.. [9. gl.
Our optimization method is as follows. First, we initialize at the uniformly weighted coreset, i.e.,
N
Wi W, W0m
m 1,....,, M.
M
Next, given step-size tuning sequence Nh, k 2 O, we update the weights each iteration & E NUU0) using the (T > O)-reuularized quasi-Newton step
i Wk + 7x (G((w) + wk)
(8)
G(w)= = Cov [g. H(w) = Cov [g. fl.
(9)
Theorem 4.3 suggests that the step size 7rk should be constant: 7 = 7 E (0, 1]; but the analysis in that section assumes that we can compute G(w) and !! (w) exactly, whereas in practice we will have to estimate them. Hence we allow Yh to depend on the iteration number in general. Finally, we project ux++1 back onto the constraint set obtain next iterate Wk+1, i.e. for m = l,..., M set
Wk+ 1.
max (wx++ 1,, 0).
(10)
3.3 Algorithm
The pseudocode for the quasi- Newsi- coreset construction method is shown in Algorithm 1. There are a number of practical considerations needed to use this method; we discuss these here.
The first consideration is that we cannot calculate G(w) and Hww) in closed form. Instead, at step k of the algorithm, we use a Monte Carlo estimate; first taking S samples (0.). ~ Twg and calculating 9. = g(0.) 1 Σ s(Or), so that we may estimate G(w))
sid.
G - S Σ N.0 E MRM
(11)
Using a Monte Carlo estimate in this way introduces one source of error. Furthermore, we often cannot sample exactly from Tue. and instead use MCMC to do so. This could lead to additional errors, for example if our samples are from unconverged chains. However, we find that this method works well practice, for a reasonable number of samples
As H H(ws) (w)) is an M x N matrix, estimating directly will incur a ((MN) storage cost. Instead, we note that we only require H (wx)(1 wx)) = Cov wx [g. f (1 wx)] = Cov we [g. h]. where h = f-1- wk- = J-g" wk. which requires O(N) time, but only O(S) space. We can then estimate H(ww) - wxk as
-u We can estimate h using h. = Σ (0.) - E fn [0)) 0f wk.
fr wx)
1 S
Σ
h.
(12)
mm

---

Algorithm QNC (QUASI-NEWTON CORESET)
Require: (Jn(0)) No, S. K, K Kunner 7s T, M.
1: Uniformly sample M data points, w.l.o.g. these correspond to indices 1,..... M. N
2: Set W0m
m M.
M
3: for k = 0, 0...., K do
4: Sample (0s),=1 Lid.
Twx x <exp (Ww [(0) To(
5: Set 9o - g(0)) - Er= ((r)) and h2 + En=1 (0()))- - fn(0.)] - -g/we wk
6: Set H4 w&) 1 Ich, and Gk E 9.0%
7:
if I Kume then
8:
Choose %/k via line search with starting value y
9:
else
Set 7k 6
11: Take a quasi-Newtone step: wx+1 +- wk + "k (Gs + +11 Hpr(1 - ws)
Project: for all m [M], Wk+1mm < max (18x+1m,
13:
return w = WK
Thus, we can find wk+1 from wi by taking the stochastic Newton step
0x41 WK. + 7he (Gs+r1) H1 wx),
(13)
and projecting onto the constraint set: Wk+11m 0k dx 1.. 20. We set the regularization parameter r by examining the condition number of G2 - and keeping it below reasonable value. We tune 7 using line search method. As we do not have access the objective function (i.e. the KL-divergence between the coreset and full posteriors), we use the curvature part of the Wolfe conditions [27] to tune this. In practice, this line search expensive. Thus, we only tune 7 for k Kume. and leave it as a constant thereafter. The intuition here is that we may start with quite poor coreset weights, and so need choose the initial steps carefully. Appendix C, we perform sensitivity analysis for the parameters S, Krune and т. We see that our results are generally not sensitive the choice of these parameters, within reasonable ranges.
Computing G + T/ involves taking S samples, forming the product 9.02. and then inverting the resulting matrix. Its time complexity is thus OSSM++M²) Computing H (wxx) wyx similarly has complexity O (SM + SN), and computing the product with (G& + +-)=- has complexity (M²). For K Newton steps Algorithm 1, the overall time complexity thus O(K(M3 + SM² + +N)), which is linear in N. The space complexity is ((M2) (to store GR). which is sublinear in N, unlike [20, 23].
In theory, In particular, can instead calculate he only for indices in I {1, (1...., M) UT, where T is uniformly selected sample S (1...., N}. Caleulating Á (wx) (1 wx) using this subsampled vector h then has complexity O(SM + ST) in the worst case, where T = [T). The overall time complexity is then OKK + SM² + ST)), which is sublinear in N if T = o(N). This could give significant improvements in coreset construction time in the large-data regime. However, this approach requires further study, as we find that it leads to a degradation in the performance of our algorithm. We do not use it in our experiments.
reduce the time complexity by using further stochastic estimate of fdx - wx).
The final practical consideration the design of Algorithm is that we cannot actually calculate the objective function that we are trying to minimize, namely the KL divergence DKL (T l between the coreset and full posteriors. Thus, it may be hard tell if the optimization is actually making progress. In practice, we monitor the norm of the gradient of the KI. divergence (which we do have access to, as given by Eq. (6)). We can terminate our algorithm early if we do not see significant enough decrease in this measure.
5

---

4 Theoretical analysis
In this section we provide a theoretical analysis of the proposed method. We demonstrate (in Theorems 4.1 and 4.2) that the optimal coreset posterior built using a uniformly random subset of data is, with high probability, an exact or near-exact approximation of the true posterior. Second, we show (in Theorem 4.3) that the proposed optimization algorithm converges exponentially quickly to point near the optimal coreset, thus realizing the guarantee Theorem 4.2 practice, and theoretically justifying our proposed algorithm. The assumptions required for our theory are somewhat technical. We give some intuition here on settings where they hold, but discuss this further Appendix A. Proofs of all results may be found in Appendix B.
4.1 Subset selection via uniform subsampling
Theorem 4.1 and Theorem 4.2 provide theoretical foundation for selection via uniform subsampling. Both recommend setting M Z d log N, where d is essentially the "dimension" of the span of the log-likelihood functions. Theorem 4.1 states that if the span of the log-likelihood functions is indeed mitte-dimensionll with dimension d' e.g., in usual exponential family model- thonel-then with high probability, the optimal coreset using uniformly random selection of M points is exact. Theorem 4.2 extends this to the general setting where O C RR N is large relative to D, and where we assume the posterior concentrates as N -> oc. In particular, Theorem 4.2 states that with high probability, the optimal coreset using uniformly random selection of M data points provides low error.
In what follows, assume that (-) are i..dd from some true data generating distribution p. and we define f(-) = E, (.)). We also introduce the new notation .) where we consider the potentials f as functions of both and 0. This enables us to make the data-dependence explicit when needed, i.e., (-) = f((nn :); for instance, in the standard context of Bayesian inference
lid. fa(+) = f( -) = log p(znl.), with In ~ Po for n E [N].
Theorem 4.1. Let So be the vector space of functions on O spanned by (f(x, -);;));; E X Assume that So is finite dimensional, with dimension d. Let be any measure on θ such that So C L²)
and denote r = Ep (Ep((fi-
- -)))). For >> define
2r
J(Ø) = inf Pr (aaJa--- ag86; [all
(14)
VNO
Then for 5 > O such that ((6) + 1, a universal constant C₁, and M E N such that
M
2d
Log.N N + log
((8) + +01)
48
J(8)
we have that with probability 2 1 8 e- MJ(6)/4
min KT((TT(( = 0. WEWN
Note that the condition that So has finite dimension < Do is satisfied for all d-dimensional exponential families when f0) log p(.,, Further, note that Theorem 4.1 still applies in situations where increases with N. To understand the behaviour of J(0), we consider in Appendix A full rank exponential family. We can show that, under certain conditions, J(8) converges strictly positive value as 8N + co. In this case, we can set ó = (1/N), such that with probability (1/N), indeed have that min (T( 0 for M 2 d dlog N.
In the interest of obtaining a more general approximation result, we now consider the case where So is infinite dimensional, but where we can find dimenssoia space SJ which approximates (0) - f(0) reasonably well. A typical setup where this is valid is when è C "R the posterior concentrates (as N + oo) at a point On θ that maximizes f(0), and each f-(-) is smooth near 8g and globally well-behaved in some sense. Theorem 4.2 below states, roughly. that for model satisfying these assumptions, min (L(T=II(T) will be small with high probability. We give the exact conditions in Appendix A. Although they are perhaps strong, the result is still indicative of when a good posterior approximation obtained by our method.
Intuitively, we require that (x, 0) be differentiable, strongly concave in 0, and "smooth enough". Letting O be the unique value such that f(0)= S f(06) vo E O, we also require that we split

---

(z; 0) into two parts,
f(,,0) (2, - (z,
E is a finiteddimensional vector space Sr. This split needs to be made such that f( is negligible. in the sense that f(z, 0) = f() 8) locally around 0g. A usual case where these conditions hold is when f(x; (2) is C in A. In this case, can obtain (x) 0) by second-order Taylor expansion,
where the span of (z.)
(x: 0) f(=;00) + (0 0o) Vo)
(2; 0o) (0 00-00 (o)'' 0o)
which has dimension bounded by 1 + ++ D 1)/2. ((2) 0) is then the remainder term, which is indeed negligible around 0g in the desired sense. In Appendix A, we provide further discussion on the conditions of Theorem 4.2. Recalling that a quantity is o(1) as N oo if and only if it converges to O, we can then state the result as follows.
Theorem 4.2. Suppose the assumptions A AI-A4 in Appendix A hold, and we set M 2 D(log N + 1).
Then as N oo, with probability least o(1), we have that
min wEWN
KL (L(TT|(=) o(1)
In the above example, we can examine the proof to obtain a bound on (L(T=ll)) of order 1/v/N with probability converging to I, showing that this result does indeed lead useful bound.
4.2 Quasi-Newton weight optimization
Theorems 4.1 and 4.2 guarantee the existence of high-quality coreset, but do not provide any insight into whether it is possible to find it tractably. Theorem 4.3 addresses this remaining gap. In particular, intuitively, asserts that as long as (1) the regularization parameter т > nonzero but small enough that it does not interfere with the optimization (in particular, the minimum eigenvalue of G(w) over the space), and (2) the optimal coreset is good approximation to the full dataset (which already guaranteed by Theorem 4.2), then the weights wa at iteration k of the approximate Newton method converge exponentially to point close to a global optimum. In the result below, let G(w) and H (w) be defined as in Eq. (9).
Theorem 4.3. Let W C W be closed convex set, W C argminwew KL(T be maximal closed convex subset, and fix the regularization parameter T > >0. Define
A (
g = inf min weW
E eigvals G(w), A 2>0 0
Suppose ] E [0, 1) and 8 2 0 such that for all w E W
ll
w0\ll S cllw-w" w1
(15)
where w² = prrpi-- (w). Then the wtth iterate we E W and its projection wt projw- (we) of the projected Newton method defined in Section 3.2 by Eqs. (8), (9) and (10), initialized at wo E W with fixed step size 7 E [0, 1] and regularization satisfies
1 (1-1)
|| - wffl x Ilwo will + 70
17 1 1 1=1-7(8-0 -))
For this result to make sense, we require that 17 E [0. 1]. We have that 1 2 O, since y E [0, 1]. c 2 and & < 1 by definition. It may be the case that 77 1, in which case result still holds, but is vacuous.
In order to use this result, the minimum positive eigenvalue of G(w) for E W C W needs to be bounded, as does the error of the optimal coresct weights w" the sense of Eq. (15). notable corollary of Theorem 4.3 occurs when there exists w" such that
M
AB E θ
Σ m=l
winn fm (6) (m).
7

---

i.e., the optimal coreset is identical to the true posterior distribution, as is guaranteed by Theorem 4.1. For any model where the conditions of Theorem 4.1 hold, such as any exponential family model, Eq. (15) holds with 8 = 0. We can thus take
to find that
1R will S (1 - - |100 - will-
This provides some intuition on Theorem 4.3; as long the optimal coreset posterior is reasonable approximation to т,, and we use a small regularization T > 0 (such that Ş = 1), then we generally need only small number of optimization steps to find the optimal coreset. The conditions of Theorem 4.1 are sufficient, but not necessary, for Eq. (15) to hold. We present further discussion in Appendix A.
We note that if O, then there is no convergence, even in the setting where Theorem 4.1 holds. However, it is important to emphasises that, for any given dataset, we are free to choose т. We can (and should) choose T to be sufficiently small that is roughly 1. In particular, this means setting significantly smaller than the minimum positive eigenvalue of G(w).
5 Experiments
In this section, we compare our proposed quasi-Newton coreset (QNC) construction against existing constructions-unifoe subsampling (UNIF), greedy iterative geodesic ascent (GIGA) [20] and itera- tive hard thresholding (IHT) [23]-as as well as the Laplace approximation (LAP), which represents what one obtains by assuming posterior normality in the large-data setting. Experiments were per- formed on a machine with a 2.6GHz 6-Core Intel Core i7 processor, and 16GB memory; code is avail- cwm/ wwww/t/tooteeeeeeeeeeeeeeeeeeeee
able at https:://itthul
-nweiineewooo--- wwweeeooorreeseeeeeeeeeeeeee wxwrssets-epeeimeeeeee experiments.
In each case, we use S = 500 Monte Carlo samples during coreset construction. We compute error metrics using 1000 samples from each method's approximate posterior. We also compare to the baseline of sampling from the full posterior (FULL) to establish a noise floor for the given comparison sample size of 1000. In the synthetic Gaussian experiment, which is the simplest one we consider, we see that sparse variational inference (SVI) [21] is prohibitively slow, taking 4.2 x 10Ss to construct a coreset of size (and scaling at best linearly with coreset size). Thus, we do not compare against here. In Appendix C we run a smaller, lower dimensional, synthetic Gaussian experiment. Here, we compare against SVI, and confirm that it is prohibitively slow for the size of datasets that we consider.
For GIGA and IHT, we need to supply a low-cost approximation to the posterior. To ensure these methods apply as generally QNC and SVI, use a uniformly sampled coreset approximation of size M with weights N/M (where M is the same as the desired coreset size we are constructing). We also use these weights for UNIF. In Appendix C we provide additional results for GIGA and IHT with a Laplace approximation used for +. This does not lead to a significant improvement, and we provide further discussion there.
Experiments are performed in settings standard in the Bayesian coresets literature, but with larger dataset sizes than those that could be considered previously. For each experiment we plot the approximate reverse KL divergence obtained by assuming posterior normality, and the build time. In Appendix C we provide comparisons on the additional metries of relative mean and covariance error, forward KL. divergence and per sample time to sample from the respective posteriors. For the experiments in Sections 5.2 and 5.3 with heavy tailed priors, we also compare using the maximum mean discrepancy [28] and kernel Stein discrepancy [29, 30J, with inverse multi-quadratic (IMQ) kernel.
Throughout the experiments, we see that QNC outperforms the other subsampling methods for most coreset sizes consider. GIGA and IHT in particular are limited by fixed choice of म and finite projection dimension defined by our choice of S = 500-we we discuss this further in Appendix C. Furthermore, outside of the synthetic Gaussian experiment, we see that QNC outperforms LAP for coreset sizes above certain thresholds, which represent a small fraction of the full dataset.
In Appendix C, we also perform a sensitivity analysis for the parameters S, Keune and that we use in Algorithm 1. We see that our results are generally not sensitive to the choice of these parameters, within reasonable ranges.
8

---

10
20
BAN Low 14
10
10
FULL
GISA
GISA
FULL GIGA
OIGA
IH
IHT
10
LAP
10"
LAP
10
10
10
Dovaser
Coreser
Connact
(a) Gaussian
(b) Sparse linear regression
100
10"
10"
10%
10
10
10
Auh Drew 92
10
10
10
FULL GIGA IHH
10
HH
10
GIGA
IH
10
CIGA
IHT
10
QNC
107
10
101
10'
Conooet
Doveser
Goreer
Conoeet
(c) Logistic regression
Basis function regression
Figure 1: Reverse KI. divergence (left) and build time in seconds (right) for each experiment. We plot the median and shaded area between the 25th percentiles over 10 random trials. Our algorithm (QNC) provides improvement in coreset quality, with comparable run-time and less user input.
5.1 Synthetic Gaussian location model
Our first comparison on Gaussian location model, with prior 0 ~ NN oo.), and likelihood
(X,) iid. (0, ol) in D dimensions. Here, we take p4 = 0 and on = Closed form expressions are available for the subsampled posterior distributions [21, Appendix BJ, and we can sample from them without MCMC. We compare the methods on synthetic dataset with N 100, 000,
N(,, ol) and set (m)) iid. NN 100) and a = 100. From Fig. la we see that our method outperforms the other subsampling methods for all coresct sizes. This example is largely illustrative; we expect LAP to provide the exact posterior here
D 100, where we generate the (X,),
by design, and this is indeed what we find (the reverse KL plots for LAP and FULL overlap).
5.2 Bayesian sparse linear regression
Next, we study a Bayesian sparse linear regression problem, where the data (I,, yr) consists of feature on RDD and an outcome Y/. E TR.N The posterior distribution is that of TR. in the model
= 0 + En.
n n= L,..,NN
(16)
where en iid N( 0), n I....,N N. We place independent Cauchy priors on the coefficients 0: 0, lid (O, rλ,), = 1,...,,, D,
where the hyperpriors are A; iid. Half-Cauchy 1), Half-Cauchy(0.. ou). with ou 2, and or I(ao, be) with ag = 1, by = 1. We perform posterior inference on the 2D + 2 dimensionnl set of parameters (0 ........ 0p. ....... .D. Sampling is performed using STAN [31]. The flight delays dataset, with N = 100, 000 and D = 13 (so that the overall inference problem is 28 dimensional). The response variable yn is the delay in the departure time of a flight, and the features are meteorological and flight-specific information. We see in Fig. 1b that our method outperforms the other subsampling methods for all coreset sizes, and LAP for sizes above 500-represnnting 0. 0.5% of the data. In Appendix C we see that the target posterior in this case has heavy tails, which makes the Laplace approximation particularly unsuited this problem. We can see the effect this has even more clearly in the additional results presented there.
dataset we study
TTis dataset was constructed by merging airport on time data from the US Bureau of Transportation Statis- tics https: //www. transtats. bts. ets.goo/DDD gov/DL_See Sew.SeeeeeFe... a records from https: wttp///nwnndddd..oo
asp? go- VQ-FGJ with historical weather

---

5.3 Heavy-tailed Bayesian logistic regression
For this comparison, we perform Bayesian logistic regression with parameter 0 E R having heavy-tailed Cauchy prior 0g iid. ~ Cauchy(O, o), = I.... D. The data (In, M/ N consists of feature In E RP and label Yn E 11,, The relevant posterior distribution is that of 0 E R which governs the generation of yn given In via
indep
Yn Fny 0
Bern
(17)
1+e -+
Again, sampling is performed using STAN, As in Section 5.2, we use the flight delays dataset, so we have that N = 100, 000 and D = 13. The response variables yn are binarized, so that yn = 1 if flight was cancelled or delayed by more than an hour, and yn -1 it was not. In Fig. Ic we see that our method outperforms other subsampling methods for all coreset sizes, and LAP for sizes above 500 representing 0.5% of the data.
5.4 Bayesian radial basis function regression
Our final comparison is on a Bayesian basis function regression example. This is a larger version of the same experiment performed by [21], with D 301 before, but N -100, 000. Full details are given in Appendix C. From Fig. 1d we see that our method provides a significant improvement over LAP for all coreset sizes, which performs particularly poorly in this experiment. despite the Gaussianity of the target posterior [21, Appendix B]. It also outperforms the other subsampling methods for coreset sizes of 1000 and above (representing 1% of the data) but this difference is not as marked as in the other examples. The fact that QNC requires large coreset to start providing benefit beyond UNIF stems from the fact that, in this higher dimensional example, there are sometimes data points that are very influential for certain basis, and both UNIF and QNC can miss these. This indicates the need (in some settings) for secondary method that can check for and include these, without the significant cost that e.g. SVI entails. We defer this study to future work. However, we do note that for this example QNC has the lowest build time, even outperforming LAP. This because we only perform small number of quasi-Newton steps before no further improvement is possible.
6 Conclusion
This paper introduces a novel method for data summarization prior to Bayesian inference. In particular, the method first selects small subset of data uniformly randomly, and then optimizes the weights on those data points using a novel quasi-Newton method. Theoretical results demonstrate that the method is guaranteed find near-optimal coreset, and that the optimal coreset has low KL. divergence to the posterior with high probability. Future work includes investigating the performance of the method in more complex models, and studying conditions under which the method can provide compression in time sublinear in dataset size.
Acknowledgments and Disclosure of Funding
Campbell was supported by a National Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant and Discovery Launch Supplement. C. Naik was supported by the Engineering and Physical Sciences Research Council and Medical Research Council [award reference 1930478]. J. Rousseau was supported by the European Research Council (ERC) under the European Union's 's Horizon 2020 research and innovation programme (grant agreement No 834175).
10

---

References
[1] Christian Robert and George Casella. Monte Carlo Statistical Methods. Springer, 2nd edition, 2004.
[2] Christian Robert and George Casella. A short history of Markov Chain Monte Carlo: subjective recollections from incomplete data. Statistical Science, 26(1):102-115, 2011.
[3] Andrew Gelman, John Carlin, Hal Stern, David Dunson, Aki Vehtari, and Donald Rubin. Bayesian data analysis. CRC Press, 3rd edition, 2013.
[4] Matthew D. Hoffiman, David M. Blei, Chong Wang, and John Paisley, Stochastic variational inference. Journal of Machine Learning Research, 14: 14:1333--1333 1347, 2013.
[5] Rajesh Ranganath, Sean Gerrish, and David Blei. Black box variational inference. In Interna- tional Conference on Artificial Intelligence and Statistics, 2014.
[6] Rémi Bardenet, Arnaud Doucet, and Chris Holmes. On Markov chain Monte Carlo methods for tall data. The Journal of Machine Leurning Reseurch, 185))11155 2017.
[7] Anoop Korattikara, Yutian Chen, and Max Welling. Austerity in MCMC land: cutting the Metropolis-Hastings budget. In International Conference on Machine Learning, 2014.
[8] Dougal Maclaurin and Ryan Adams. Firefly Monte Carlo: exact MCMC with subsets of data. In Conference on Uncertainty in Artificial Intelligence, 2014.
[9] Max Welling and Yee Whye Teh. Bayesian learning via stochastic gradient Langevin dynamics. In International Conference on Machine Learning, 2011.
[10] Sungin Ahn, Anoop Korattikara, and Max Welling. Bayesian posterior sampling via stochastic gradient Fisher scoring. In International Conference on Machine Learning, 2012.
[11] Matias Quiroz, Robert Kohn, and Khue-Dung Dang. Subsampling MCMC- an introduction for the survey statistician. Sankhya: The Indian Journal of Statistics, 80-A:S33-S69, 2018.
[12] Joris Bierkens, Paul Fearnhead, and Gareth Roberts. The zig-zag process and super-efficient sampling for bayesian analysis of big data. The Annals of Statisties, 47(3): 2019.
[13] Murray Pollock, Paul Fearnhead, Adam M Johansen, and Gareth O Roberts. Quasi-stationary monte carlo and the scale algorithm. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82((): 10n--116712222 2020.
[14] James Johndrow, Natesh Pillai, and Aaron Smith. No free lunch for approximate MCMC. arXiv:2010. 12514, 2020.
[15] Tigran
Nagapetyan, Andrew Duncan, Leonard Hasenclever, Sebastian Vollmer, Lukasz Szpruch, and Konstantinos Zygalakis. The true cost of stochastic gradient Langevin dynamics. (Tf:::170602692 2017.
[16] Michael Betancourt.
The fundamental incompatibility of Hamiltonian Monte Carlo and data subsampling. In International Conference on Machine Learning, 2015.
[17] Matias Quiroz, Robert Kohn, Mattias Villani, and Minh-Ngoc Tran. Speeding up MCMC by efficient data subsampling. Journal of the American Statistical Association, 114(526):831-844 2019.
[18] Jonathan Huggins, Trevor Campbell, and Tamara Broderick. Coresets for scalable Bayesian
logistic regression. In Advances in Neural Information Processing Systems, 2016.
[19] Trevor Campbell and Tamara Broderick. Automated scalable Bayesian inference via Hilbert coresets. Journal of Machine Learning Research, 20(15), 2019.
[20] Trevor Campbell and Tamara Broderick. Bayesian coreset construction via greedy iterative
geodesic ascent. In International Conference on Machine Learning, 2018
[21] Trevor Campbell and Boyan Beronov. Sparse variational inference: Bayesian coresets from scratch. In Advances in Neural Information Processing Systems, 2019.
[22] Dionysis Manousakas, Zuheng Xu, Cecilia Mascolo, and Trevor Campbell. Bayesian pseudo- coresets. In Advances in Neural Information Processing Systems, 2020.
[23] Jacky Zhang, Rajiv Khanna, Anastasios Kyrillidis, and Oluwasanmi Koyejo. Bayesian coresets: revisiting the nonconvex optimization perspective. In Artificial Intelligence in Statistics, 2021.
11

---

[24] Martin Jankowiak and Du Phan. Surrogate likelihoods for variational annealed importance sampling, arXiv preprint 1.Ki::211..1294, 2021.
[25]
Naitong Chen, Zuheng Xu, and Trevor Campbell. Bayesian inference via sparse hamiltonian flows. arXiv preprint arXiv:2203.05723. 2022.
[26] Charles Broyden. Quasi-Newton methods. In W Murray, editor, Numerical Methods for Unconstrained Optimization, pages 87-106. Academic Press, 1972.
[27] Philip Wolfe. Convergence conditions for ascent methods. SIAM review, 11(2):226-23 235, 1969.
[28] Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. The Journal of Machine Learning Research, 13(1):723-773, 2012.
[29] Qiang Liu, Jason Lee, and Michael Jordan. A kemelized Stein discrepancy for goodness-of-fi tests. In International conference machine learning, pages 276 784. PMLR, 2016.
[30] Anna Korba, Pierre-Cyril Frank-Frankwwski, Szymon Majewski, and Pierre Ablin. Kernel Stein discrepancy descent. arXiv preprint arXiv:2105. 09994, 2021.
[31] Bob Carpenter, Andrew Gelman, Matthew Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of statistical software, 76(1):1-32, 2017.
12

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] The claims made in the abstract are supported by the theoretical and experimental results in Section 4 and Section 5
(b) Did you describe the limitations of your work? [Yes] Methodological limitations discussed in Section 3.3 mainly the fact that we need to used Monte Carlo estimates of some quantities of interest. The theoretical results in Section 4 show the setting in which our method can be expected to work. We give some intuition for what settings these are in the main text, with further discussion in Appendix A. our fourth experiment in Section 5.4 we also discuss some of the limitations of our work in high dimensions.
(c) Did you discuss any potential negative societal impacts of your work? [No] Our work provides generic preprocessing algorithm for sampling from probability distribu- tions whose log density involves many summed terms. Our method is generic and foundational nature: has many downstream applications, including Bayesian infer- ence, which may have a societal impact depending on the particular data and model under consideration. But we do not speculate on the impacts of potential downstream applications in this work.
(d) Have you read the ethics review guidelines and ensured that your paper conforms them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] The assump- tions for Theorem 4.1 and Theorem 4.3 are given in the theorem statements. Before the statement of Theorem 4.2 we give intuitive explanations of its assumptions, but the full assumptions are explicitly laid out in Appendix A.
(b) Did you include complete proofs of all theoretical results? [Yes] Proofs of all theoretical results are included in Appendix
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main ex- perimental results (either in the supplemental material URL)? [Yes] Code, data, and instructions needed to reproduce all experimental results is provided in the supplemental material.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they chosen)? [Yes] Important details provided in Section 5, and all further details are included in the code provided as part of the supplemental material.
(c) Did you report error bars (e.g.. with respect the random seed after running experi- ments multiple times)? [Yes] Error bars representing the 25th percentiles over 10 random trials presented for each plot of the results.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] The details of the machine used to perform all experiments are given in Section
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] Creators of all existing assets are cited.
(b) Did you mention the license the assets? [Yes] All the datasets we use are publicly available, and some do include a license. Where relevant, the license is mentioned.
(c) Did you include any new assets either the supplemental material or as URL? [Yes] We include new code in the supplemental material.
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [No] The datasets we use are publicly available, and do not pertain to identifiable individuals. We obey the license terms of the data where applicable, do not discuss whether consent was obtained from the data creators.
13

---

(c) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [No] None of the datasets use contain personally identifiable information or offensive content.
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
14


Reviews:
[{"review_id": "m8pQBwOn6C", "reviewer": "Reviewer_LhtQ", "summary": "This paper presents a novel algorithm for approximating posterior distributions via weighted subsets of the data. The proposed algorithm samples points uniformly at random and then computes the weights using a quasi-Newton method; crucially it scales well with respect to the number of data points, both in terms of time and storage. The authors provide a rigorous theoretical analysis of their proposed algorithm and empirically evaluate it against commonly-used methods for posterior approximation. ", "questions": "My understanding of the empirical evaluation is that all four settings satisfy the assumptions of either Theorem 4.1 or 4.2: is this correct and if so, have the authors considered experimenting on settings where those assumptions do not hold? It would be illuminating to see the empirical performance when the theoretical guarantees do not apply. ", "limitations": "Yes, I believe the authors' discussion of their work's limitations is sufficient; related to my question above, I would have found it helpful if the authors had provided some intuition about settings in which their theoretical results would not hold much like how they did for settings in which the assumptions do hold. ", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, I found this paper to be an excellent read and believe it to be a strong candidate for acceptance. \n\n1. The proposed method is not particularly creative but frankly, I think the simplicity of the ideas being presented here are a strength of this work. The sum of the contributions in my mind are sufficient to get this paper above the originality bar for acceptance. \n\n2. The theoretical and empirical analysis are compelling although I do have a few suggestions/concerns regarding the experiments that if addressed, could raise my score of this submission even higher: \n    - I would have liked to see some comparison against SVI, even if that meant scaling down the experiments to settings where SVI is computationally tractable. It would be very telling if QNC was competitive with or even outperformed SVI in small-data regimes; if not, then finding some sort of cost-accuracy trade-off between the two methods could be useful for practitioners. \n   - I would have liked to see some sort of sensitivity analysis to the model hyperparameters, specifically S, K_tune (about which very little is said) and tau. \n\n3. The manuscript is very well-written; I genuinely found this paper to be an enjoyable read and applaud the authors' effort to include intuition wherever possible. \n\n4. I believe this work will have significant impact on a non-trivial percentage of the NeurIPS population: the proposed algorithm is simple and easy-to-implement and thus, accessible to people who might want to use/expand upon it. Furthermore, the scalability makes this algorithm broadly applicable, in settings where previously-proposed methods might not be. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "NHJfVuTkGoT", "reviewer": "Reviewer_nwJJ", "summary": "This paper addresses the problem of approximating a Bayesian posterior with a large number of datapoints. They do so by selecting a sparse, weighted subset of datapoints (a coreset). Previous methods have addressed this problem, but the authors note that previous methods suffer from requiring some subset of: (1) potentially expensive loops over the whole dataset, (2) user input or parameter tuning, (3) a large number of potentially expensive MCMC simulations. The authors address (1) by noting that a coreset can be formed on a random subset of the data, (2) by proposing a quasi-Newton method that requires one MCMC simulation per iteration, and (3) just by virtue of providing a different algorithm without the same tuning / user input requirements. \n\nThe authors prove conditions under which (1) still allows for exact or approximate recovery of the posterior. And they prove conditions under which their optimization algorithm (2) does not require many iterations. In experiments, the authors show that their method is both faster and more accurate than other coreset methods.", "questions": "- Can the authors clarify what's going on with $\\xi$ in Theorem 4.3 (see weaknesses above)?\n\n- (minor question): On lines 118-122, the authors note that when M is larger than the dimension of the space of log-likelihood functions, $G$ is not invertible, which motivates adding some small $\\tau$ to its diagonal. If $G$ has a zero eigenvalue in some directions, this implies the weight vector $w$ has an orthogonal decomposition into directions along which its variation does and does not matter. Could a reparameterization allow the optimization algorithm to instead just work along the directions for which the variation of $w$ does matter? This might avoid the use of $\\tau$. ", "limitations": "I think the authors have adequately discussed the social impact of their work.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n- The writing in the paper is good; the material is pretty technical, but I found it relatively easy to follow. \n- The method proposed by the authors intuitively fixes clearly explained problems with previous work and has desirable theoretical properties (mostly -- see weaknesses below) and this intuition and theory actually seem to hold up in practice in non-exhaustive, but at least detailed and careful experiments. That's a pretty solid combination. Based on this, I vote for the paper to be accepted. \n\n**Weaknesses**\n1. The method is claimed to be \"black box\" with \"far less ... user input required\" than other algorithms. I think this is overselling things a little bit. I see five different hyperparameters (at least five, since one of them is a sequence): the number of coreset points $M$, the step size sequence $\\gamma_k$, the Hessian regularization $\\tau$ (or, equivalently, the condition number threshold mentioned on line 142), the number of iterations to tune the stepsize for $K_{tune}$, and the size of the subsampling set $T$. Since there are no sensitivity studies for these hyperparameters, I don't think it's fair to suggest that little user input is required. Relatedly, some competing methods are criticized for needing an initial, user-specified posterior approximation $\\hat\\pi$ as an input, such as the Laplace approximation. I don't think this requires much user input; given access to an optimization library and automatic differentiation tools, a Laplace approximation can be constructed without user input.\n2. Theorem 4.3, which studies the convergence rate of the authors' proposed algorithm (an important selling point of the work), essentially assumes the smallest eigenvalue of the covariance $G(w)$ is not too small. For example if $\\xi$, is zero, and the optimal coreset perfectly approximates the posterior, then Theorem 4.3 predicts no convergence of their algorithm. But it seems -- and maybe the authors can correct me on this -- that as $N \\to \\infty$ we should expect $\\xi \\to 0$. For example, say we expect the posterior to concentrate to a point mass as $N \\to \\infty$ and some $w^* \\in W$ perfectly approximates the posterior. If the posterior is contracting to a point mass, won't $G(w^*)$ approach the matrix of 0's? (I'm writing this as if the inf over $W$ is actually the inf over $W_N$; if it's the inf over $W$, then the inf contains arbitrarily large $w$'s even for $N = 1$).\n\n**Smaller things**\n\n- There are a lot of derivatives with respect to $w$ taken throughout the paper. But $w$ sits in a constrained set $W$. How are these derivatives defined when $w$ is on the boundary of $W$?\n- Should the left hand side of equations 6 and 7 be indexed $1:M$ and $1:M, 1:M$, respectively?\n- Line 113: it's not intuitively clear to me why the first term should dominate the expression here.\n- Algorithm 1: is input $K$ supposed to be $K_{tune}$?\n- \"This is sublinear in $N$\". I usually take \"sublinear\" to mean $o(N)$, whereas this runtime is $O(N)$.\n- In the equation right after line 188, I don't think $A$, $K, X$, and $X_1$ were ever defined.\n- In Theorem 4.3: could the inf over $w \\in W$ be replaced by the inf over ${ w \\in W :  \\|\\| w - w^* \\|\\| \\leq \\|\\|w_0 - w^* \\|\\| \\}$? Since there's monotone improvement towards $w^*$, it seems like you'll never leave this ball, and so you can ignore any bad behavior outside of it.\n- In the experiments, it would be good to include the runtime of FULL to see how marginal / non-marginal the gains here are.\n- In the proof of Theorem 4.3, there are a couple places where $G(w) + \\lambda I$ is written. Should this $\\lambda$ be a $\\tau$?\n- It would be helpful if the equations in the appendix were contiguously numbered with those in the main text. Otherwise it's not immediately obvious which equation \"Eq. (1)\" refers to.\n- The letter $t$ shows up in a few places in proofs. It seems like this is a number between 0 and 1, but I don't think this is ever stated.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "HvsCI1xLnZC", "reviewer": "Reviewer_hdj8", "summary": "In the area of Bayesian Coresets, this paper presents an algorithm that proposes:\n1. To uniformly select $M$ data points that will be used to  approximate the posterior distribution, and,\n1. to use a quasi-Newton method to refine the weights that will be assigned to each of the selected data points.\n\nFor this quasi-Newton method, the authors propose using an estimated covariance matrix to approximate the objective’s Hessian. \nFinally, the authors provide theoretical justification for their algorithm via three different theorems.", "questions": "The most important thing I would like the authors to clarify is the derivation of equation (6) in the Appendix because it is critical for the proof of both, Thm. 4.1 and Thm. 4.2.\n\nIn the proof of Thm. 4.3, I think it should be that $\\mathbf u \\in \\mathbb R^M$ is in the null space of $G(\\mathbf w)$ iff $\\mathbf{u}^{\\mathrm{T}}f(\\theta) = c$ holds   $\\pi_w$-almost everywhere. \n\nFor the following limitations, I'm not asking the authors to solve them. Instead, I suggest the authors that, if my assessment is correct, they should make the limitations more explicit as open problems.\n\nFirst, I see the following limitations:\n\n- The use of Monte Carlo estimations. \nThis has two implications:\n    1. The obvious one is the Monte-Carlo error, which might imply a bias.\n    2. A not-so-obvious one is that the method assumes that we have an oracle that allows us to sample from the corresponding distribution. However, often times we can only get approximated samples from the distribution (or samples from an unconverged chain).\n\nPart of point (1) was briefly mentioned, but there is no discussion about (2). To be concrete, it should be pointed out that a version of Thm. 4.3 with the two sources of errors, is still an open problem, extending the discussion at the end of section 3.2 \n\nAnother limitation is the fact that we don't know if the optimization is actually making progress. A quick look at Algorithm 1 will show that there is no function of the objective being used as stopping criteria. I understand why this happens. Nevertheless, I think the authors should be, again, more explicit about this. This is related to the fact that the Armijo condition is not being used.", "limitations": "The authors mention some of the paper’s limitations. Above I'm suggesting mentioning a few more.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I think this is a very solid submission.\nThe authors propose an algorithm and prove that their approximation might behave as the full posterior and provide a convergence rate.\nA notable feature is that assumptions are deeply discussed in the body of the paper (Section 4) or the appendix.\nIn general, I would say that it is a very precise paper.\nOn the weak side, I found the proof in the appendix quite hard to follow. Many inequalities were not evident at all and limited the possibility of doing a complete verification of the paper’s contributions [cf. Questions below].\nIt is important to highlight that when I was able to derive the bounds, the proofs are correct, as in Thm. 4.3.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "Ajqq-WHI-Ad", "reviewer": "Reviewer_No1L", "summary": "This paper proposes a method for approximating a posterior distribution in the framework of Bayesian statistics. Specifically, this posterior distribution involves a sum over $N$ potentials $\\sum_{n=1}^N f_n(\\theta)$ which has to be evaluated each time an expectation wrt the posterior has to be calculated. It is proposed to replaced this sum of another weighted sum over $M \\ll N$ terms.  This corresponds to selecting a subset of data, called a coreset. \nThese points are sampled uniformly at random and their weights are computed by minimizing the KL divergence between the full posterior and the approximated posterior. A Quasi-Newton optimization method is proposed to calculate these weights.\nStatistical guarantees are given for the size of the coreset in order to reach a small KL divergence. The convergence of the Quasi-Newton method is also discussed. Numerical simulations illustrate the accuracy of the posterior approximation.", "questions": "I have several questions about Thm 4.1 and its proof. \n- First, the statement of this result is difficult to interpret since $J(\\delta)$ depends on $\\delta$ in a non-trivial way. After Thm 4.1, the authors try to explain the behaviour of $J(\\delta)$ at the lines 185-191.\nI do not find this explanation very clear. Although we might guess what are the random variable $X$ and $X_1$ from the context, it is worth defining them explicitly.\nThe limiting behaviour for $J(\\delta)$ at line 189 is not obvious. Can you give a few lines of explanation?\n- Second, it is difficult to follow the proof of Thm 4.1 in Supplementary Material. Could you give a roadmap for this proof? What are the main steps? It would help the reader to read this proof.\n\nHere are a few technical questions about the proof of Thm 4.1. \n- At line 61 (SM), the authors introduce a variable $t$, which is presumably in the $[0,1]$ interval. It is worth defining $t$.\nNext, at line 62 (SM), there is a formula for the KL divergence. Can you give a few words of explanation? Why is this integral equal to the KL divergence?\n\n\n- Also, at line 75 (SM), the authors claim they use a union bound. This step necessitates a clarification. Could you give a few details?\nWhat are the events that you consider in the union? How do you bound $\\phi$?\n\n- At line 79 (SM), there is a typo. At the right-hand side of first inequality, there should be $\\Pr( \\dots < t_N)$ but $<t_N$ is missing. The same thing happen at the second inequality. \n\n- At line 82 (SM), I presume there is another typo. In my opinion, $\\log(A_d^{-d/2})$ should be replaced by $\\log(A_d^{d/2})$. Also, the last equality should be an inequality. \n\nAlso, I have a few questions about the proof of Thm 4.2\n\n- At line 98 (SM), the authors upper bound $I_{\\epsilon, 1}$. The last inequality is unclear. I presume that the authors claim that \n$[g - \\sum_m v_m (f_m - \\bar{f})-E_{w(v),t}]^2\\leq [g - \\sum_m v_m (f_m - \\overline{f})]^2$. \nWhy is this correct? \nIs the quantity $g - \\sum_m v_m (f_m - \\bar{f}) - E_{w(v),t}$ positive?\n\n- At line 101 and at several other places (e.g. line 113), it is written that the probability of an event is bounded from below by $1 - f(N)$, or some decreasing function $f(N)$. To be complete, it should be mentioned that $N$ has to be large enough so that $f(N)<1$.\n\n- At line 100 (SM), the authors write *We now assume $\\|S_N(\\theta)\\| \\leq (2- \\sqrt{2}) \\epsilon L_0 \\sqrt{N}$*. This is not an assumption since the authors calculate right after the probability that this event does not happen. See also line 120 for a similar use of the verb *assume*.\n\n- At line 102 (SM), another inequality is unclear. At the second lower bound, it seems that the authors take the square of an inequality. Are both sides of this inequality positive numbers? \n\n- At line 127 (SM), the first equality is unclear. Could you give some details?\n\nHere are some questions about Thm 4.3.\n\n- At line 229, I think that the assumptions are incomplete, namely, about the parameter $\\eta = 1 - \\gamma(\\xi - \\epsilon)$ which influences the convergence rate. What are the conditions on these parameters such that $\\eta >0$? Shouldn't there be conditions on $\\epsilon$ so that $\\xi - \\epsilon >0$?\n\n- About the convergence guarantees, there is an upper bound on $\\| w_k - w_k^\\star\\|$ where $w_k^\\star$ is the projection of $w_k$ on $W^\\star$. Here, $w_k$ is the weight at iteration $k$. What does this bound tell about the convergence to a global optimum?\n\n\nTypos:\n\n104 unfinished sentence\n\nAlgorithm 1: indicies &larr; indices\n\nSM 34 extra full stop\n\nSM 87 $S_N$ &larr; $S_N(\\theta_0)$\n\nSM 122 missing ) on the left-hand side\n\n=========================================================================\n\nNote: I have read the authors' response which I found convincing. I update my score.\n", "limitations": "The limitations of this work are summarized in the checklist and are briefly discussed in section 5.4. I do not see any potential negative societal impact.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Originality: This paper discusses interesting ideas. Being familiar with coresets in general but less with this specific application, I cannot completely assess the originality of this paper.\n    \nQuality: From a global perspective, the main paper is well-written and the results are well presented. However, from the technical viewpoint, the interpretation of some mathematical claims is not straightforward, for instance, in the case of Thm 4.1. The numerical simulations are convincing.\n    \nClarity: The main paper is easy to follow. Concerning the technical details, I have several questions about the proofs of the theoretical results in Supplementary Material. These proofs are long (e.g., the proof of Thm 4.2 takes more than 8 pages). In principle, this is not an issue. However, these proves deserve extra clarifications at several places and are not easy to read; see the questions below. To guide the reader, it would also be beneficial if a brief description of the proof strategy is given before each proof.\n    \nSignificance: This method is certainly interesting and probably useful in the context of Bayesian inference.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept