General bounds on the quality of Bayesian coresets

Large-scale data is now commonplace in scientific and commerical applications of Bayesian statistics. But despite its prevalence, and the corresponding wealth of research dedicated to scalable Bayesian inference, there are still suprisingly few general methods that provably provide inferential results, within some reasonable tolerated error, at a significant computational cost savings. Exact Markov chain Monte Carlo (MCMC) methods require many full passes over the data [1, Ch. 6–12, 2, Ch. 11–12], limiting the utility of these methods when even a single pass is expensive. A wide range of MCMC methods that access only a subset of data per iteration, e.g., via delayed acceptance [3, 4, 5, 6], pseudomarginal or auxiliary variable methods [7, 8, 9], and basic subsampling [10, 11, 12, 13], provide at most a minor improvement over full-data MCMC [14, 15, 16]. On the other hand, methods including carefully constructed log-likelihood function control variates can provide substantial gains [17, 18, 19]. However, black-box control variate constructions for large-scale data often rely on assumptions such as posterior density differentiability and unimodality that do not hold in many popular models, e.g., those with discrete variables or multimodality. See [15, 20] for a survey of scalable MCMC methods. Parametric approximations via variational inference [21] or the Laplace approximation [22, 23] can be obtained scalably using stochastic optimization methods, but existing general theoretical guarantees for these methods again typically rely on posterior normality assumptions [24, p. 141–144,25, 26, 27, 28, 29, 30] (see [21, 31] for a review).

Although many existing methods rely on asymptotic normality or unimodality in the large-scale data regime, the problem of handling large-scale data in Bayesian inference does not fundamentally require this structure. Instead, one can more generally exploit redundancy in the data (i.e., the existence of good approximate sufficient statistics), which can be used to draw principled conclusions about a large data set based only on a small fraction of examples. Indeed, while approximate posterior normality often does not hold in models with latent discrete or combinatorial objects, weakly identifiable or unidentifiable parameters, persisting heavy tails, multimodality, etc., such models can and regularly do exhibit significant redundancy in the data that can be exploited for faster large-scale inference. Bayesian coresets [32]—which involve replacing the full dataset during inference with a sparse weighted subset—are based on this notion of exploiting data redundancy. Empirical studies have shown the existence of high-quality coreset posterior approximations constructed from a small fraction of the data, even in models that violate posterior normality assumptions and for which standard control variate techniques work poorly [33, 34, 35, 36, 37]. However, existing theoretical support for Bayesian coresets in the literature is limited. There exist no lower bounds on Bayesian coreset approximation error, and while upper bounds do exist, they currently impose restrictive assumptions. In particular, the best available theoretical upper bounds to date apply to exponential family models [36, 38] and models with strongly log-concave and locally smooth log-densities [37].

This article presents new theoretical techniques and results regarding the quality of Bayesian coreset approximations. The main results are two general large-data asymptotic lower bounds on the KL divergence (Theorems 3.3 and 3.5), as well as a general upper bound on the KL divergence (Theorem 5.3) under the assumption that the log-likelihoods satisfy a multivariate generalization of subexponentiality (Definition 5.2). The main general results in this paper lead to various novel insights about specific Bayesian coreset construction methods. Under mild assumptions,

• •

  common importance-weighted coreset constructions (e.g. [32]) require a coreset size $M$ proportional to the dataset size $N$ (Corollary 4.1), even with post-hoc optimal weight scaling (Corollary 4.2), and thus yield a negligible improvement over full-data inference;

• •

  any construction algorithm requires a coreset size $M>d$ when the log-likelihood function is determined by $d$ parameters locally around a point of concentration (Corollary 4.3);

• •

  subsample-optimize coreset construction algorithms (e.g. [38, 36, 37, 39]) achieve an asymptotically bounded error with a coreset size $\mathrm{polylog}N$ in a wide variety of models (Corollary 6.1).

The paper includes empirical validation of the main theoretical claims on two models that violate common assumptions made in the literature: a multimodal, unidentifiable Cauchy location model with a heavy-tailed prior, and an unidentifiable logistic regression model with a heavy-tailed prior and persisting posterior heavy tails. Experiments were performed on a computer with an Intel Core i7-8700K and 32GB of RAM. Proofs of all theoretical results may be found in Appendix A.

Notation. We use standard asymptotic growth symbols $O,\Omega,\Theta,o,\omega$ (see, e.g., [40, Sec. 3.3]), and their probabilistic variants $O_{p},\Omega_{p},\Theta_{p},o_{p},\omega_{p}$ (see, e.g., [24, Sec. 2.2]). We use the same symbol to denote a measure $\pi$ and its density $\pi(\cdot)$ with respect to a specified dominating measure. We also regularly suppress integration variables and differential symbols in integrals throughout for notational brevity when these are clear from context; for example, $\int\pi\exp(\ell)$ is shorthand for $\int\pi(\mathrm{d}\theta)\exp(\ell(\theta))$. Finally, the pushforward of a measure $\pi$ by a map $\eta$ is denoted simply $\eta\pi$.

Define a target probability distribution $\pi$ on a space $\Theta$ comprised of a sum of $N$ potentials $\ell_{n}:\Theta\to\mathbb{R}$, $n=1,\dots,N$ and a base distribution $\pi_{0}(\mathrm{d}\theta)$,

where the normalization constant $Z$ is not known. In the Bayesian context, this distribution corresponds to a Bayesian posterior distribution for a statistical model with prior $\pi_{0}$ and conditionally i.i.d. data $X_{n}$, where $\ell_{n}(\theta)=\log p(X_{n}|\theta)$. The goal is to compute or approximate expectations under $\pi$; but the likelihood $\ell$ (and its gradient) becomes expensive to evaluate when $N$ is large. To avoid this cost, Bayesian coresets [32, 33, 34, 35, 36, 37] involve replacing the target with a surrogate density

where $w\in\mathbb{R}^{N}$, $w\geq 0$ are a set of weights, and $Z(w)$ is the new normalizing constant. If $w$ has at most $M\ll N$ nonzeros, the $O(M)$ cost of evaluating $\sum_{n}w_{n}\ell_{n}$ (and its gradient) is a significant improvement upon the original $O(N)$ cost. In this work, the problem of coreset construction is formulated in the data-asymptotic limit; a coreset construction method should

• •

  run in $o(N)$ time and memory (or at most $O(N)$ with a small leading constant),

• •

  produce a small coreset of size $M=o(N)$,

• •

  produce a coreset with $O(1)$ posterior forward/reverse KL divergence as $N\to\infty$.

These three desiderata ensure that the effort spent constructing and sampling from the coreset posterior is worthwhile: the coreset provides a meaningful reduction in computational cost compared with standard Markov chain Monte Carlo algorithms, and has a bounded approximation error.

This section presents lower bounds on the KL divergence of coreset approximations for general models and data generating processes. The first key steps in the analysis are to write all expectations in terms of distributions that do not depend on $w$, and to remove the difficult-to-control influence of the tails of $\pi$ and $\pi_{w}$ by restricting certain integrals to some small subset $B\subseteq\Theta$ of the parameter space. Lemma 3.1, the key theoretical tool used in this section, achieves both of these two goals; note that the result has no major assumptions and applies generally in any setting that a Bayesian coreset can be used. For convenience, define

and the decreasing, nonnegative function $f:\mathbb{R}_{+}\to\mathbb{R}_{+}$,

Note that while the integrals in the fraction denominator in $J_{B}(w)$ range over the whole $\Theta$ space, a further lower bound on $\underline{\operatorname{KL}}(w)$ can be obtained by restricting their domains arbitrarily. Also, crucially, the bound in Lemma 3.1 does not depend on $\pi_{w}(B^{c})$, which would be difficult to analyze without detailed knowledge of the tail behaviour of $\pi_{w}$ as a function of the coreset weights $w$. Although the bound in Lemma 3.1 applies generally, it is most useful when $B$ is small (so that simple local approximations of $\ell$ and $\ell_{w}$ can be used), $\pi$ concentrates on $B$ (so that $\pi(B^{c})\approx 0$), and $\pi$ and $\pi_{w}$ are very different when restricted to $B$; the behaviour of the bound in this case is roughly (see the proof in Appendix A) $f(J_{B}(w))\approx-\log(1-\operatorname{TV}(\pi,\pi_{w}))$. Finally, note that Lemma 3.1 remains valid if one replaces $\ell_{w}$ with $\ell_{w}-c$ and $\ell$ with $\ell-c^{\prime}$ for any constants $c,c^{\prime}$ that do not depend on $\theta$ but may depend on the data and coreset weights $w$.

For the remainder of this section, consider the setting where $\Theta$ is a measurable subset of $\mathbb{R}^{d}$ for some $d\in\mathbb{N}$, fix some $\theta_{0}\in\Theta$, and assume each $\ell_{n}$ is differentiable in a neighbourhood of $\theta_{0}$. Let

Theorems 3.3 and 3.5 characterize KL divergence lower bounds in terms of the sum of the coreset weights $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ and the log-likelihood gradients $g,g_{w}$. Intuitively for the full data set where all $w_{n}=1$ and $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}=N$, and an i.i.d. data generating process from the likelihood with parameter $\theta_{0}$, the central limit theorem asserts under mild conditions that $g_{w}/\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\overset{p}{\to}0$ at a rate of $N^{-1/2}$. Theorems 3.3 and 3.5 below provide KL lower bounds when the coreset construction algorithm does not match this behavior. In particular, Theorem 3.3 provides results that are useful when $g_{w}/\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\overset{p}{\to}0$ occurs reasonably quickly but slower than $N^{-1/2}$, while Theorem 3.5 strengthens the conclusion when $g_{w}/\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\overset{p}{\to}0$ very slowly or not at all. The major benefit of Theorems 3.3 and 3.5 for analyzing coreset construction methods is that they reduce the problem of analyzing posterior KL divergence to the much easier problem of analyzing the 2-norm $\|\cdot\|_{2}$ of a weighted sum of random vectors in $\mathbb{R}^{d}$.

Consider a sequence $r\to 0$ as $N\to\infty$, and for a fixed matrix $H\succ 0$ let

be a sequence of neighbourhoods around $\theta_{0}$; these will appear in Assumptions 3.2, 3.4, 3.3 and 3.5 below. Note that throughout, all asymptotics will be taken as $N\to\infty$, and various sequences (e.g., $r$ and $B$) are implicitly indexed by $N$. To simplify notation, this dependence is suppressed. Assumption 3.2 makes some weak assumptions about the model and data generating process: it intuitively asserts that the potential functions are sufficiently smooth around $\theta_{0}$, that $r\to 0$ slowly, and that $\pi$ concentrates at $\theta_{0}$ at a usual rate. Note that Assumption 3.2 does not assume data are generated i.i.d. and places no conditions on the coreset construction algorithm.

Two additional assumptions related to the coreset construction algorithm—namely, that it works well enough that $\frac{1}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}\sum_{n}w_{n}\nabla^{2}\ell_{n}(\theta)\overset{p}{\to}H$ and $g_{w}/\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\overset{p}{\to}0$ at a rate faster than $r\to 0$—lead to asymptotic lower bounds on the best possible quality of coresets produced by the algorithm, as well as lower bounds even after optimal post-hoc scaling of the weights.

Theorem 3.3 is restricted to the case where the coreset algorithm is performing reasonably well. Theorem 3.5 extends the bounds to the case where the algorithm is performing poorly, in the sense that it is unable to make $\frac{g_{w}}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}\overset{p}{\to}0$ at a rate faster than $r\to 0$ (or perhaps $\frac{g_{w}}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}$ does not converge to 0 at all). In order to draw conclusions in this setting, we need a weak global assumption on the potential functions. A function $f:\Theta\to\mathbb{R}$ is $L$-smooth below at $\theta_{0}$ if

Note that $L$-smoothness below is weaker than Lipschitz smoothness and does not imply concavity; Eq. 15 restricts the growth of the function only in the negative direction, and only when the expansion is taken at $\theta_{0}$. Assumption 3.4 asserts that the potential functions are smooth below.

Theorem 3.5 uses Assumptions 3.2 and 3.4 and additional assumptions related to the coreset construction algorithm to obtain lower bounds in a setting that relaxes the “performance” conditions in Theorem 3.3: $-\frac{1}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}\sum_{n}w_{n}\nabla^{2}\ell_{n}(\theta)$ no longer needs to converge to $H$ in probability, and $g_{w}/\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ can converge to 0 slowly or not at all.

An important final note in this section is that while Theorems 3.3 and 3.5, as stated, require choosing $\Theta$ to be some measurable subset of $\mathbb{R}^{d}$ and that the posterior $\pi$ concentrates around some point of interest $\theta_{0}\in\mathbb{R}^{d}$, these results can be generalized to a wider class of models and spaces. In particular, Corollary 3.6 demonstrates that if $\Theta$ is arbitrary, but the potential functions $\ell_{n}$ only depend on $\theta$ through some other function $\eta:\Theta\to\mathbb{R}^{d}$, that the conclusions of Theorems 3.3 and 3.5 still hold.

In this section, the general theoretical results from Section 3 are applied to specific algorithms, Bayesian models, and data generating processes to explain previously observed empirical behaviour of coreset construction, as well as to place fundamental limits on the necessary size of coresets. Consider a setting where the data $X_{n}$ arise as an i.i.d. sequence drawn from some probability distribution $\nu$, $\ell_{n}(\eta(\theta))=\log p(X_{n}|\eta(\theta))$ for $\eta:\Theta\to\mathbb{R}^{d}$, $\eta_{0}=\eta(\theta_{0})$, and the following technical criteria hold (where $\mathbb{E}$ denotes expectation under the data generating process):

• (A1)

  $\mathbb{E}\left[\nabla\ell_{n}(\eta_{0})\right]=0$ and $H=\mathbb{E}\left[-\nabla^{2}\ell_{n}(\eta_{0})\right]=\mathbb{E}\left[\nabla\ell_{n}(\eta_{0})\nabla\ell_{n}(\eta_{0})^{T}\right]\succ 0$.

• (A2)

  $\mathbb{E}\left[\|\nabla\ell_{n}(\eta_{0})\|_{2}^{2+\delta}\right]<\infty$ for some $\delta>0$ and $\mathbb{E}\left[\|\nabla^{2}\ell_{n}(\eta_{0})\|^{2}_{F}\right]<\infty$.

• (A3)

  On a neighbourhood of $\eta_{0}$, $\|\nabla^{2}\ell_{n}(\eta)-\nabla^{2}\ell_{n}(\eta_{0})\|_{2}\leq R(X_{n})\|\eta-\eta_{0}\|_{2}$, $\mathbb{E}\left[R(X_{n})\right]<\infty$.

• (A4)

  $\eta\pi_{0}$ is twice differentiable a neighbourhood of $\eta_{0}$, and $\pi(\eta_{0})>0$.

• (A5)

  For all $r\to 0$ such that $r^{2}=\omega\mathopen{}\mathclose{{}\left(\log N/N}\right)$, $-\log\eta\pi\mathopen{}\mathclose{{}\left(\|\eta-\eta_{0}\|>r}\right)=\Omega_{p}(Nr^{2})$.

These conditions apply to a wide range of models, e.g., an unidentifiable, multimodal location model posterior with heavy tails on $\Theta=\mathbb{R}$, where the Bayesian model is specified by

and the data are generated from the likelihood with parameter $\theta_{0}=5$, and an unidentifiable logistic regression posterior with heavy tails on $\mathbb{R}^{2}$, where the Bayesian model is specified by

the covariates are generated via $X_{n}\overset{\text{iid}}{\sim}{\mathrm{Unif}}(\{x\in\mathbb{R}^{2}:\|x\|_{2}\leq 1\})$, and the observations $Y_{n}$ are generated from the likelihood with parameter $\theta_{0}=\begin{bmatrix}1&6\end{bmatrix}^{T}$. See Proposition A.6 in Appendix A for the verification of (A1-5) for these two models. Example posterior log-densities for these models are displayed in Fig. 1.

This section presents upper bounds on the KL divergence of coreset approximations. As in Section 3, the first step is to write all expectations in terms of distributions that do not depend on $w$. Lemma 5.1 does so without imposing any major assumptions; the result again applies generally in any setting that a Bayesian coreset can be used. For convenience, define

The upper bound in Lemma 5.1 is nonvacuous (i.e., finite) as long as there exists a $\alpha>1$ such that the $\alpha$ Rényi divergence $D_{\alpha}(\pi_{w}||\pi)$ [49, p. 3799] is finite. Note that as in Lemma 3.1, the bound in Lemma 5.1 remains valid if one replaces $\ell_{w}$ with $\ell_{w}-c$ and $\ell$ with $\ell-c^{\prime}$ for any constants $c,c^{\prime}$ that do not depend on $\theta$ but may depend on the coreset weights $w$ and data.

More practical bounds necessitate an assumption about the behaviour of the potentials $(\ell_{n})_{n=1}^{N}$. Definition 5.2 below asserts that the multivariate moment generating function of $(\ell_{n})_{n=1}^{N}$ is bounded when the vector is close to 0. This definition is a generalization of the usual definition of subexponentiality for the univariate setting (e.g., [50, Sec. 2.7]). Theorem 5.3 subsequently shows that Definition 5.2 is sufficient to obtain simple bounds on $\overline{\operatorname{KL}}$.

Definition 5.2, the key assumption in Theorem 5.3, is satisfied by a wide range of models when choosing $h(x)=x$ and $A\propto\operatorname{Cov}_{\pi}\mathopen{}\mathclose{{}\left((\ell_{n})_{n=1}^{N}}\right)$, as demonstrated by Proposition 5.4. Because this case applies widely, let $A$-subexponential be shorthand for $(h,A)$-subexponentiality with $h(x)=x$.

In other words, intuitively, if a coreset construction algorithm produces weights such that $\operatorname{Var}_{\pi}(\bar{\ell}_{w}-\bar{\ell})$ is small, then $\overline{\operatorname{KL}}(w)$ is also small. That being said, the generality of Definition 5.2 to allow arbitrary $h,A$ is still helpful in obtaining upper bounds in specific cases; see, e.g., Propositions A.1 and A.2.

A strategy to construct Bayesian coresets that has recently emerged in the literature, shown in Algorithm 3, is to first subsample the data to select $M$ data points, and subsequently optimize the weights for those selected data points [36, 38, 37]. The subsampling step serves to pick a reasonably flexible basis of log-likelihood functions for coreset approximation, and avoids the slow greedy selection routines from earlier work [33, 35, 34]. The optimization step tunes the weights for the selected basis, avoiding the poor approximations of importance-weighting methods. Indeed, Algorithm 3 creates exact coresets $\pi_{w^{\star}}=\pi$ with high probability in Gaussian location models [36, Prop. 3.1] and finite-dimensional exponential family models [37, Thm. 4.1], and near-exact coresets with high probability in strongly log-concave models [37, Thm. 4.2] and Bayesian linear regression [38, Prop. 3].

Corollary 6.1 generalizes these results substantially, and demonstrates that coresets of size $M=O(\mathrm{polylog}(N))$ produced by the subsample-optimize method in Algorithm 3 maintain a bounded KL divergence as $N\to\infty$. Two key assumptions are subexponentiality of the potentials and a polynomial (in $N$) growth of $\operatorname{Var}_{\pi}(\ell(\theta))$; these conditions are not stringent and should hold for a wide range of Bayesian models and i.i.d. data generating processes. The last key assumption in Eq. 30 is that a randomly-chosen potential function $\ell_{I}$, $I\sim{\mathrm{Categorical}}(p_{1},\dots,p_{N})$ (with probabilities as in Algorithm 3) is well-aligned with the residual coreset error function. Similar alignment conditions have appeared in past results for more restrictive settings (see, e.g., $J(\delta)$ in [37, Thm. 4.1]).

Fig. 3 confirms that subsample-optimize coreset construction methods applied to the logistic regression and Cauchy location models in Eqs. 18 and 19 (which both violate the conditions of past upper bounds in the literature) are able to provide high-quality posterior approximations for very small coresets—in this case, $M\propto\log N$.

This article presented new general lower and upper bounds on the quality of Bayesian coreset approximations, as measured by the KL divergence. These results were used to draw novel conclusions regarding importance-weighted and subsample-optimize coreset methods, which align with simulation experiments on two synthetic models that violate the assumptions of past theoretical results. Avenues for future work include general bounds on the subexponentiality constant $\beta$ in Proposition 5.4, as well as the alignment probability in Eq. 30, in the setting of Bayesian models with i.i.d. data generating processes. A limitation of this work is that both quantities currently require case-by-case analysis.

The author gratefully acknowledges the support of an NSERC Discovery Grant (RGPIN-2019-03962).