# Bayes Hilbert Spaces for Posterior Approximation

George Wynne  University of Bristol

###### Abstract

Performing inference in Bayesian models requires sampling algorithms to draw
samples from the posterior. This becomes prohibitively expensive as the size
of data sets increase. Constructing approximations to the posterior which are
cheap to evaluate is a popular approach to circumvent this issue. This begs
the question of what is an appropriate space to perform approximation of
Bayesian posterior measures. This manuscript studies the application of Bayes
Hilbert spaces to the posterior approximation problem. Bayes Hilbert spaces
are studied in functional data analysis in the context where observed
functions are probability density functions and their application to
computational Bayesian problems is in its infancy. This manuscript shall
outline Bayes Hilbert spaces and their connection to Bayesian computation, in
particular novel connections between Bayes Hilbert spaces, Bayesian coreset
algorithms and kernel-based distances.

###### Contents

  1. 1 Introduction
  2. 2 Background
  3. 3 Measure Discrepancies and Bayes Hilbert Spaces
  4. 4 Bayesian Coresets and Bayes Hilbert Spaces
  5. 5 Maximum Mean Discrepancy and Bayes Hilbert Spaces
  6. 6 Hilbert Bayesian Coresets
  7. 7 Kullback-Leibler Bayesian Coresets
  8. 8 Conclusion
  9. 9 Appendix

##  1 Introduction

The aim of this manuscript is to advocate for, and demonstrate the utility of,
Bayes Hilbert spaces as tools to investigate computational Bayesian problems.
Bayesian statistics is a common statistical modelling paradigm which involves
the fusion of prior knowledge about the unknown quantity to be inferred with
observed data. The result is a posterior measure over the unknown quantity.
The posterior measure is the object of interest since it can be used to
perform prediction and uncertainty quantification. Practically though, using
the posterior measure for these tasks involves deploying sampling algorithms
on the posterior measure. The practical issue is that it is expensive to
sample from the posterior measure when the observed data set is large. The
typical cost when one has observed $N$ data points is $O(N)$ per iteration of
the sampling algorithm, which becomes prohibitive since many thousands of
iterations could be required to obtain the desired number of samples. For more
on Bayesian approaches to modelling see Gelman et al. (2013) and for sampling
in Bayesian methods see the handbook Kroese et al. (2011).

There are a range of solutions to this problem. The two main approaches are
designing algorithms which avoid the $O(N)$ cost each iteration and forming an
approximation to the target posterior and sampling from that approximation
instead of the target posterior. Examples of the former approach involve
sampling algorithms which employ sub-sampling (Bierkens et al., 2019; Dang et
al., 2019), variational inference (Blei et al., 2017; Hoffman et al., 2013)
and divide and conquer methods (Scott et al., 2016; Srivastava et al., 2015;
Vyner et al., 2023). The latter approximation approach will be the focus of
this paper. This approach typically involves keeping the prior the same and
approximating the likelihood function with something cheaper to evaluate than
$O(N)$, resulting in a cheaper per-iteration cost when sampling. The area of
this methodology which shall be the focus in the sequel is Bayesian coresets
(Huggins et al., 2016; Campbell and Broderick, 2019). This is where the an
approximation of the likelihood is formed based on $M$ points and weights,
with $M\ll N$. This approximation results in an $O(M)$ per-iteration cost once
inserted into a sampling algorithm. Many variants of this method exist
(Campbell and Broderick, 2018, 2019; Campbell and Beronov, 2019; Manousakas et
al., 2020; Naik et al., 2022).

When performing this approximation, indeed when performing any posterior
approximation, one must ask themselves what is an appropriate space to form
the approximation? This choice plays a critical role in the analysis and
performance of posterior approximations. It will impact both the notion of
distance between the approximate posterior and the target posterior and how
one performs optimisation with respect to this distance to form the
approximation. The aim of this paper is to show Bayes Hilbert spaces are
appropriate spaces to perform such analysis and facilitate a novel
perspectives on Bayesian coreset algorithms and posterior approximation in
general.

Bayes Hilbert spaces (van den Boogart et al., 2011; van den Boogaart et al.,
2014; Barfoot and D’Eleuterio, 2023) are spaces of measures which have notions
of addition and scalar multiplication which are coherent with Bayes’s theorem.
This makes them particularly appealing spaces to form posterior
approximations. The notion of a Bayes Hilbert space has origins in
compositional data analysis (Aitchison, 1982; Pawlowsky-Glahn and Buccianti,
2011) and is now a common tool used in distributional data analysis (Mateu and
Giraldo, 2021; Petersen et al., 2022). This is a sub-field of functional data
analysis where a user observes a data set which is a collection of probability
density functions and wishes to perform analysis. To deal with the specific
structure of probability density functions, different versions of addition and
scalar multiplication are used in the Bayes Hilbert space than typical
function spaces. The use of Bayes Hilbert spaces in distributional data
analysis is mature with many applications (Mateu and Giraldo, 2021). On the
other hand, the use of Bayes Hilbert spaces in computational Bayesian problems
is in its infancy. The primary reference of this application is the innovative
contribution of Barfoot and D’Eleuterio (2023) which made clear how
variational inference procedures can be viewed in terms of projections in
Bayes Hilbert spaces. This connection was made by showing variational
inference methods can be viewed as optimising coefficients of an approximating
function in a Bayes Hilbert space.

The main technical contributions of this manuscript are relating the Bayes
Hilbert space norm to commonly used distances on measures, providing a novel
connection between Bayes Hilbert spaces and kernel methods and framing
multiple Bayesian coreset algorithms in terms of Bayes Hilbert spaces.

Section 2 will provide background content for this manuscript. In particular,
Section 2.1 will introduce the mathematical details of Bayes Hilbert spaces,
Section 2.2 will describe the sense in which Bayes Hilbert spaces are coherent
with Bayes’ theorem and Section 2.3 will outline existing literature related
to Bayes Hilbert spaces. After the background content, Section 3 will provide
a novel result which relates the distance in a Bayes Hilbert space to common
discrepancies for measures, Section 4 will outline how Bayesian coresets can
be expressed in a Bayes Hilbert space, Section 5 will provide a novel link
between maximum mean discrepancy and Bayesian coresets, Section 6 will
leverage this link to provide novel connections between maximum mean
discrepancy algorithms and Hilbert coreset algorithms, Section 7 will describe
a novel connection between Bayesian coresets constructed using the Kullback-
Leibler divergence and Bayes Hilbert spaces. Concluding remarks and future
avenues of research are outlined in Section 8.

##  2 Background

This section shall provide background information to motivate and frame the
novel results that will occur in the later sections. Section 2.1 shall outline
Bayes Hilbert spaces which are the central focus of this manuscript, Section
2.2 shall outline Bayes’ theorem and Section 2.3 shall discuss lines of
research related to the Bayes Hilbert spaces methodology, with a focus on
methods which map measures into Hilbert spaces.

###  2.1 Bayes Hilbert Spaces

The aim of this section is to motivate and introduce the mathematical details
of Bayes Hilbert spaces. This will include defining the elements of a Bayes
Hilbert space, the notions of addition and scalar multiplication and the
Hilbert inner product.

Before beginning with the mathematical details of a Bayes Hilbert space it is
helpful to understand how the ideas developed. The ideas behind a Bayes
Hilbert space may be traced back to compositional data analysis (Aitchison,
1982). Mathematically, compositional data is non-negative multivariate data
with a sum constraint that represents relative, rather than absolute,
information. This data type is appropriate when one wishes to analyse the
composition of certain objects, hence the name compositional. For example, one
may be studying a collection of chemicals which are each made up from a fixed
number of other chemicals. Then compositional data analysis can be used to
perform analysis on the relative proportions of each of the chemicals. For a
long list of other examples see Aitchison (1986).

This type of data has many quirks. For example, a coherent notion of addition
and subtraction is not straightforward. This is because the typical notion of
vector addition would lead to a notion of subtraction which could result in
negative values, and negative values cannot represent proportions. Issues
therefore also occur if one were to try and use the typical notion of scalar
multiplication with this type of data.

The pioneering paper by Aitchison (1982) derived a mathematical framework to
deal with the quirks of compositional data outlined above. These tools are
known as the Aitchison geometry and form the bedrock for the field of
compositional data analysis. The geometry revolves around deriving particular
notions of equality, addition, scalar multiplication and subtraction for
compositional data along with a logarithmic transform which maps the data from
the simplex into a nicer space. The application areas of this mathematical
framework are legion, covering areas such as geology, ecology and microbiome
research (Pawlowsky-Glahn and Buccianti, 2011; Filzmoser et al., 2018;
Aitchison, 1986; van den Boogaart and Tolosana-Delgado, 2013) as well as an
insightful survey paper marking 40 years since Aitchison’s original
publication (Greenacre et al., 2022).

The idea of a Bayes Hilbert space is to generalise the Aitchison geometry used
in compositional data analysis to the situation where the objects of interest
are measures. The easiest way to conceptualise this move is to view
compositional data as a histogram, with each the $n$-th entry of a
compositional vector representing the probability of the $n$-th event. Then,
the trick is to view a probability density function as a continuous limit of a
histogram and to then adapt the Aitchison geometry to this continuous limit,
which in practice means moving from finite sums over the entries of a
compositional vector to integrals with respect to probability density
functions. Then the move from probability density functions to measures is
made by using Radon-Nikodym derivatives.

The rest of this section will be spent making these ideas and intuition
concrete. All technical content is taken from existing sources (Hron et al.,
2022; van den Boogaart et al., 2014; van den Boogart et al., 2011; Maier et
al., 2021).

Let $(\Theta,\mathcal{A})$ be a measurable space and $\mu$ a probability
measure on $\Theta$. This measure $\mu$ will act as the base measure. What
follows can be generalised to finite measures that are not probability
measures, the difference is simply additional normalisation terms. Elements of
$\Theta$ will be denoted by $\theta$. Define $M(\mu)$ as the set of measures
on $(\Theta,\mathcal{A})$ that are $\sigma$-finite and mutually absolutely
continuous with respect to $\mu$, which means that they have the same null
sets as $\mu$. Discussion regarding the choice of $\mu$ for Bayesian
computation problems is given in Section 4.

Define an equivalence relation $=_{B}$ between two measure $\eta,\nu\in
M(\mu)$ as $\eta=_{B}\nu$ if and only if there exists a constant $c>0$ such
that $\eta(A)=c\nu(A)\>\forall\>A\in\mathcal{A}$. It can be easier to
understand this notion of equivalence in terms of the Radon-Nikodym
derivatives with respect to $\mu$. For a measure $\eta\in M(\mu)$ the Radon-
Nikodym derivative of $\eta$ with respect to $\mu$ is the non-negative
function $\frac{\mathrm{d}\eta}{\mathrm{d}\mu}$ such that
$\eta(A)=\int_{A}\frac{\mathrm{d}\eta}{\mathrm{d}\mu}\mathrm{d}\mu\>\forall
A\in\mathcal{A}$. The existence of such a function is guaranteed by the Radon-
Nikodym theorem, see for example Bogachev (2007, Theorem 3.2.2). The function
is unique $\mu$-almost everywhere. The interpretation of this result is simply
that measures in $M(\mu)$ possess what can be viewed as densities with respect
to $\mu$ which makes interpretation of results to come more straightforward.
To this end, for $\eta\in M(\mu)$ set
$p_{\eta,\mu}\coloneqq\frac{\mathrm{d}\eta}{\mathrm{d}\mu}$ as the Radon-
Nikodym derivative of $\eta$ with respect to $\mu$ for notational convenience.
Note that $p_{\mu,\mu}=1$ in the sense that $p_{\mu,\mu}$ is the function
constantly equal to one. Then $\eta=_{B}\nu$ if and only if there exists a
constant $c>0$ such that $p_{\eta,\mu}=cp_{\nu,\mu}$ $\mu$-almost everywhere.
Indeed, anytime that these Radon-Nikodym derivatives are written as equal it
is to be interpreted in the $\mu$-almost everywhere sense.

Define $B(\mu)$ to be the set of equivalence classes with respect to $=_{B}$
within $M(\mu)$. With this notion of equality in place, which is helpful when
studying Bayesian problems as described in Section 4, the next step is to
define addition and scalar multiplication. For $\eta,\nu\in B(\mu)$ addition
is defined

$\displaystyle(\eta\oplus\nu)(A)\coloneqq\int_{A}p_{\eta,\mu}(\theta)\cdot
p_{\nu,\mu}(\theta)\mathrm{d}\mu(\theta)$

for every $A\in\mathcal{A}$, For $\eta,\nu\in B(\mu)$ and
$\alpha\in\mathbb{R}$ scalar multiplication is defined as

$\displaystyle(\alpha\odot\eta)(A)\coloneqq\int_{A}p_{\eta,\mu}(\theta)^{\alpha}d\mu(\theta)$

for every $A\in\mathcal{A}$. The combination of addition and scalar
multiplication facilitates the identification of an additive inverse, namely
for $\nu\in B(\mu)$ define $\ominus\nu\coloneqq(-1\odot\nu)$ and for
$\eta,\nu\in B(\mu)$ define $\eta\ominus\nu\coloneqq\eta\oplus(\ominus\nu)$.

While these definitions may look abnormal at first, all of the operations
$\oplus,\odot,\ominus$ can be interpreted straightforwardly through the Radon-
Nikodym derivatives. Namely, for any $\eta,\nu\in B(\mu)$ the operation
$\oplus$ can be written as $p_{\eta,\mu}\oplus p_{\nu,\mu}\coloneqq
p_{\eta,\mu}\cdot p_{\nu,\mu}$ where $\cdot$ means standard multiplication so
that $p_{\eta\oplus\nu,\mu}=p_{\eta,\mu}\oplus p_{\nu,\mu}$. Scalar
multiplication $\odot$ can be written as $\alpha\odot p_{\eta,\mu}\coloneqq
p_{\eta,\mu}^{\alpha}$ so that $p_{\alpha\odot\eta,\mu}=\alpha\odot
p_{\eta,\mu}$. Finally, $\ominus$ can be written as
$p_{\eta\ominus\nu,\mu}=p_{\eta,\mu}/p_{\nu,\mu}$ so that
$p_{\eta\ominus\nu,\mu}=p_{\eta,\mu}-p_{\eta,\mu}$. These operations defined
directly on the Radon-Nikodym derivatives make the natural identification
$\eta$ with $p_{\eta,\mu}$ for all $\eta\in B(\mu)$ coherent in the sense that
if one writes $\eta=_{B}p_{\eta,\mu}\>\forall\eta\in B(\mu)$ then
$\eta\oplus\nu=_{B}p_{\eta\oplus\nu,\mu}$ for all $\eta,\nu\in B(\mu$),
similarly for $\odot,\ominus$. Therefore, for ease of notation $\eta$ may be
replaced with $p_{\eta,\mu}$ at certain places in the sequel as at times the
focus shall be on the Radon-Nikodym derivatives rather than the measures
themselves.

The next result links all these operations together and assures us that these
operations result in a valid real vector space.

######  Theorem 1.

(van den Boogart et al., 2011, Theorem 5) $B(\mu)$ equipped with
$\oplus,\odot$ is a real vector space with $\mu$ the additive zero element.

This vector space is known as a Bayes linear space (van den Boogart et al.,
2011) and provides the basic structure needed to proceed to Bayes Hilbert
spaces. The next step is defining the following set

$\displaystyle B^{2}(\mu)=\left\\{\eta\in
B(\mu)\colon\mathbb{E}_{\mu}\left[\left\lvert\log\frac{\mathrm{d}\eta}{\mathrm{d}\mu}\right\rvert^{2}\right]<\infty\right\\}.$
(1)

Note that $\eta\in B^{2}(\mu)$ if and only if $\log p_{\eta,\mu}\in
L^{2}(\mu)$, where $L^{2}(\mu)$ is the typical space of (equivalence classes
of) functions on $\Theta$ that are square integrable with respect to $\mu$.

As previously mentioned, the central tool for compositional data analysis was
a logarithm transform to map the compositional data into a nice space. The
analogy to this for the present case is the centred log-ratio (CLR) transform
which is defined (Egozcue et al., 2006; van den Boogaart et al., 2014)

$\displaystyle\Psi_{\mu}(\eta)=\log p_{\eta,\mu}-\mathbb{E}_{\mu}[\log
p_{\eta,\mu}],$

where the expectation is being taken with respect to the input argument of
$\log p_{\eta,\mu}$. A few remarks are in order. First, the expectation is
finite given the assumption that $\mu$ is a finite measure and that $\eta\in
B^{2}(\mu)$. Second, since $p_{\eta,\mu}$ accepts as argument an element of
$\Theta$ so to does the CLR $\Psi_{\mu}(\eta)$. Third, $\Psi_{\mu}(\mu)=0$, in
the sense that it is the function constantly zero, since $p_{\mu,\mu}=1$
meaning that the base measure, which is the additive identity in $B^{2}(\mu)$
is mapped to the zero function. The CLR maps into $L^{2}_{0}(\mu)=\\{f\in
L^{2}(\mu)\colon\mathbb{E}_{\mu}[f]=0\\}$, the elements of $L^{2}(\mu)$ that
have zero mean, which is the analogy to the nice space that is used in
compositional data analysis. This space is also used in information geometry
as a tangent space to manifolds of measures, see Section 2.3 for more
discussion. Finally, the integrability condition in $B^{2}(\mu)$ is very weak
and means that $B^{2}(\mu)$ can contain infinite measures. Discussion on this
is given in Section 8.

The operations $\oplus,\odot$ relate addition with multiplication and scalar
multiplication with exponentiation, these relations are maintained by
logarithms which is why the CLR transform is helpful. Also, the subtraction of
the expectation of the logarithm ensures the scale invariance property in the
definition of $=_{B}$ holds. Overall, this means for any $\eta,\nu\in
B^{2}(\mu)$ and $\alpha\in\mathbb{R}$ the following linearity holds
$\Psi_{\mu}(\alpha\odot(\eta\oplus\nu))=\alpha\cdot\left(\Psi_{\mu}(\eta)+\Psi_{\mu}(\nu)\right)$
and that for any $\eta,\nu\in B^{2}(\mu)$ such that $\eta=_{B}\nu$ the CLR
maps $\eta,\nu$ to the same output, meaning
$\Psi_{\mu}(\eta)=\Psi_{\mu}(\nu)$. The inverse function of the CLR map is
simply $\exp$, to see this take any $\eta\in B^{2}(\mu)$ then

$\displaystyle\exp(\Psi_{\mu}(\eta))=_{B}p_{\eta,\mu}\cdot\exp(-\mathbb{E}_{\mu}[\log
p_{\eta,\mu}])=_{B}p_{\eta,\mu}$

due to the scale invariance of $=_{B}$ (van den Boogart et al., 2011), where
the natural identification between $\eta$ and $p_{\eta,\mu}$, as described
immediately before Theorem 1, has been used.

With the CLR established an inner product structure can be defined. For
$\eta,\nu\in B^{2}(\mu)$ define

$\displaystyle\langle\eta,\nu\rangle_{B^{2}(\mu)}\coloneqq\langle\Psi_{\mu}(\eta),\Psi_{\mu}(\nu)\rangle_{L^{2}(\mu)}=\int_{\Theta}\Psi_{\mu}(\eta)(\theta)\Psi_{\mu}(\nu)(\theta)\mathrm{d}\mu(\theta).$
(2)

with corresponding norm

$\displaystyle\lVert\eta-\nu\rVert_{B^{2}(\mu)}=\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{2}(\mu)}.$
(3)

The next result shows the inner product on $B^{2}(\mu)$ provides the desired
Hilbertian structure.

######  Theorem 2.

(van den Boogaart et al., 2014, Theorem 1) The map $\Psi_{\mu}$ is an isometry
between $B^{2}(\mu)$ and $L^{2}_{0}(\mu)=\\{f\in
L^{2}(\mu)\colon\mathbb{E}_{\mu}[f]=0\\}$. The inverse of $\Psi_{\mu}$ is
$\exp$ and $B^{2}(\mu)$ is a Hilbert space.

The space $B^{2}(\mu)$ equipped with the inner product (2) is called a Bayes
Hilbert space. Theorem 2 reveals critical structure of $B^{2}(\mu)$. This
Hilbertian structure will facilitate approximations based on dictionaries of
functions and optimisation methods based on orthogonal projections. Such
structure is typically not present in common representations of measures and
therefore Bayes Hilbert spaces have great potential for applications in
measure approximation. For example, depending on the choice of $\mu$ it can be
straightforward to derive an orthonormal basis for $B^{2}(\mu)$ by using
orthogonal polynomials. The case of $\mu$ being a Gaussian and using Hermite
polynomials is studied in Barfoot and D’Eleuterio (2023) in the context of
Gaussian variational inference. More discussion regarding the potential use of
these basis approximations for computational Bayesian problems is given in
Section 8.

###  2.2 Bayes’ Theorem and Bayes Hilbert Spaces

The aim of this section is to explain how Bayes Hilbert spaces have structure
which is coherent with Bayes’ theorem, justifying their name. In particular,
it will be shown that the notions of addition and scalar multiplication used
in $B^{2}(\mu)$ are coherent with using Baye’s theorem to update a prior
belief using likelihoods to obtain a posterior. This section contains no novel
technical content with exposition taken from van den Boogart et al. (2011);
van den Boogaart et al. (2014).

One may go back to compositional data analysis, the genesis of Bayes Hilbert
spaces as outlined in Section 2.1, to see how the connection to Bayes’s
theorem occurs. Aitchison noted that the simplex, the canonical sample space
for compositional data analysis, is “familiar in other areas of statistics
$\ldots$ as the operation of Bayes’s formula to change a prior probability
assessment into a posterior probability assessment through the perturbing
influence of the likelihood function” (Aitchison, 1986; van den Boogart et
al., 2011).

These allusions will now be made concrete. The following recites content from
van den Boogart et al. (2011) and for more on Bayes’ theorem consult Stuart
(2010). The sample space will be $\Theta$ and the prior measure $\pi_{0}$.
Suppose there are $N$ observations $\\{x_{n}\\}_{n=1}^{N}\subset\mathcal{X}$,
the data space, that are independently and identically distributed. Assume the
likelihood for each observation is $l$ which takes as input an observation
from $\mathcal{X}$ and a parameter from $\Theta$. Then, assuming the
observations are independent and identically distributed, the likelihood for
all the data is $L(\theta)\coloneqq\prod_{n=1}^{N}l(\theta,x_{n})$, or using
shorter notation $L=\prod_{n=1}^{N}l_{x_{n}}$ where
$l_{x_{n}}(\theta)\coloneqq l(\theta,x_{n})$. Bayes’ theorem states the
posterior measure $\pi$ over $\Theta$, given the prior, likelihood and
observations, satisfies

$\displaystyle\frac{\mathrm{d}\pi}{\mathrm{d}\pi_{0}}=Z^{-1}L,$ (4)

where $Z=\mathbb{E}_{\pi_{0}}[L]$ is a constant and is known as the evidence.
It is implicitly assumed that the prior, likelihood and data are such that
$Z<\infty$ which is a very mild and common assumption.

Mathematically, (4) shows that the Radon-Nikodym derivative of the posterior
with respect to the prior is proportional to the likelihood. Therefore, using
the definitions of $=_{B}$ and $p_{\pi,\pi_{0}}$ given in Section 2.1

$\displaystyle p_{\pi,\pi_{0}}=_{B}L=_{B}\bigoplus_{n=1}^{N}l_{x_{n}}.$ (5)

This means shows that Bayes’ theorem is simply a sum over the likelihood terms
when written in Bayes Hilbert space notation. This is the sense in which the
Bayes Hilbert space operations are coherent with Bayes’ theorem. Another
helpful property of a Bayes Hilbert space is the definition of $=_{B}$ as
being equality up to scalar constant. This is useful since in practice
posterior sampling algorithms are agnostic of scalar constants and therefore a
mathematical framework for posterior approximation should also be invariant to
scalar constants.

######  Example 1.

An example of Bayesian inference that will recur is logistic regression.
Suppose one observes input output pairs $x_{n}=\\{u_{n},y_{n}\\}$ with
$u_{n}\in\mathbb{R}^{d}$ and $y_{n}\in\\{0,1\\}$. For example, $u_{n}$ might
represent medical information about the $n$-th person in a population and
$y_{n}$ indicates whether they do or do not have a certain illness. A model
for the regression problem of predicting $y_{n}$ given $u_{n}$ is logistic
regression (Gelman et al., 2013, Chapter 3.7). The parameter space is
$\Theta=\mathbb{R}^{d}$ and for parameter choice $\theta$ the response $y_{n}$
is modelled as
${y_{n}\sim\emph{Bernoulli}(1/1+e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}})}$
meaning the likelihood is

$\displaystyle
l_{x_{n}}(\theta)=\frac{y_{n}}{1+e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}+\frac{(1-y_{n})e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}{1+e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}.$

A commonly used prior for $\theta$ is $\pi_{0}=\emph{N}(0,I_{d})$, the
standard multivariate Gaussian.

Under the assumption that $\pi\in B^{2}(\mu)$ the CLR for $\pi$ is

$\displaystyle\Psi_{\mu}(\pi)=\log p_{\pi,\mu}-\mathbb{E}_{\mu}[\log
p_{\pi,\mu}].$

Under the additional assumption $\pi_{0}\in B^{2}(\mu)$, using the chain rule
for Radon-Nikodym derivatives $p_{\pi,\mu}=p_{\pi,\pi_{0}}p_{\pi_{0},\mu}$ the
CLR transform can be written

$\displaystyle\Psi_{\mu}(\pi)$ $\displaystyle=\log
p_{\pi,\pi_{0}}-\mathbb{E}_{\mu}[\log p_{\pi,\pi_{0}}]+(\log
p_{\pi_{0},\mu}-\mathbb{E}_{\mu}[\log p_{\pi_{0},\mu}])$
$\displaystyle=\mathcal{L}-\mathbb{E}_{\mu}\left[\mathcal{L}\right]+(\log
p_{\pi_{0},\mu}-\mathbb{E}_{\mu}[\log p_{\pi_{0},\mu}]),$ (6)

where $\mathcal{L}\coloneqq\log L=\sum_{n=1}^{N}\log l_{x_{n}}$ is the log-
likelihood. Note that if $\mu=\pi_{0}$ then the final term in brackets is zero
since $p_{\pi_{0},\pi_{0}}=1$ and the CLR transform becomes equal to the
centred log-likelihood function.

######  Example 2.

Taking $\mu=\pi_{0}$ then continuing Example 1 the CLR for the logistic
regression model is

$\displaystyle\Psi_{\pi_{0}}(\pi)$
$\displaystyle=\sum_{n=1}^{N}\log\left(\frac{y_{n}}{1+e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}+\frac{(1-y_{n})e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}{1+e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}\right)$
$\displaystyle-\sum_{n=1}^{N}\mathbb{E}_{\pi_{0}}\left[\log\left(\frac{y_{n}}{1+e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}+\frac{(1-y_{n})e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}{1+e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}\right)\right].$

The importance of these representations from a computational Bayesian point of
view should be considered. The main focus in computational Bayesian statistics
is the posterior measure. A computational Bayesian typically wishes to draw
samples from the posterior and will almost always have to use a sampling
algorithm. The typical cost is $O(N)$ per iteration of the sampling algorithm,
with many thousands of iterations typically required, therefore the overall
cost becomes prohibitive for large data sets. In this large data scenario one
usually employs either a sampling algorithm which makes cheap approximations,
such as subsampling data during the execution of the sampling algorithm, or
will try and form a cheap approximation of the posterior and sample from that
instead. The difficulty in the latter approach is that measures are not
typically viewed as lying in a space that is applicable to standard
approximation and optimisation theory.

The Bayes Hilbert space provides appealing structure in which to view the
posterior. First, it possesses an appropriate notion of equality as it is
agnostic to scaling constants, much like common sampling algorithms. Second,
the Hilbertian structure facilitates the use of approximation theory,
concentration inequalities and projections which will all play a role in the
rest of this manuscript to construct approximations to the posterior. Later,
Theorem 3 in Section 3 shows that the norm in a Bayes Hilbert space can bound
commonly used distances between measures and therefore approximation in a
Bayes Hilbert space results in minimising commonly used distances.

###  2.3 Related Work

The central idea of Bayes Hilbert spaces is to construct a Hilbert space into
which measures shall be mapped. The Hilbert space provides nice structure to
form distances, construct approximations and perform optimisation. The
application of Bayes Hilbert spaces to computational Bayesian problems beyond
expressing Bayes rule in terms of Bayes Hilbert space operations (van den
Boogart et al., 2011) is currently modest. This is because the application of
Bayes Hilbert spaces has so far focused on distributional data analysis
(Petersen et al., 2022; Mateu and Giraldo, 2021). The sole contribution to the
problem of posterior approximation using a Bayes Hilbert space is Barfoot and
D’Eleuterio (2023) which frames variational inference in terms of a Bayes
Hilbert space.

Other frameworks exists which embed measures into other function spaces. In
information geometry, the most similar topic are exponential manifolds
(Pistone and Sempi, 1995; Cena and Pistone, 2006; Pistone, 2013). Informally,
this is where a manifold is constructed over the space of probability measures
such that at a given measure, call it the reference measure, a Banach space is
constructed and a subset found which is homeomorphic to the space of measures
that are absolutely continuous with respect to the reference measure. The
connection between this method and Bayes Hilbert spaces is that $B^{2}(\mu)$
plays the roles of the reference measure space and $L^{2}_{0}(\mu)$ the Banach
space, in this context known as the tangent space. In fact, the CLR transform
also features in this method, see Pistone and Sempi (1995, Equation 23). In
Pistone and Sempi (1995) the Banach space was an Orlicz space which enforces
stronger integrability conditions than $L^{2}_{0}(\mu)$ to ensure that the
corresponding measures are all finite, unlike in $B^{2}(\mu)$. Connections
between information geometry and compositional data analysis have been made by
Erb and Ay (2021) but extending these to Bayes Hilbert spaces was left for
future work.

Orlicz spaces are rather technical and other work has focused on constructing
similar manifolds which use Hilbert spaces instead. This includes Newton
(2012) which also uses $L^{2}_{0}(\mu)$ as its Hilbert space but uses a map
slightly different from the CLR transform, see Newton (2012, Equation 6). A
reproducing kernel Hilbert space approach, see Section 5, was performed by
Fukumizu (2009) in which the Banach space is a reproducing kernel Hilbert
space. The CLR transform appears again in Fukumizu (2009, Theorem 18.1). This
approach is appealing since reproducing kernel Hilbert spaces are very nice
spaces to study and manipulate. However, one main difference is that the set
of measures absolutely continuous with respect to a given base measure is only
homeomorphic to a subset of a Hilbert space. This means it does not inherit
the vector space structure from the Hilbert space. This is in contrast to how
Bayes Hilbert spaces have vector space structure.

The purpose of all these methods in information geometry, as the name of the
field suggests, is to better understand the geometry of spaces of measures
rather than to directly approximate them. The closest to tackling an
approximation problem is Sriperumbudur et al. (2017) which uses the
reproducing kernel Hilbert space approach of Fukumizu (2009) to construct
density estimators from samples of a distribution, which is still a distinct
task to posterior approximation. Another distinguishing factor between the
Bayes Hilbert space approach and the exponential manifold approach is that in
the latter the space of measures is bijective with a subset of a function
space, not a subspace like the former, since they are tangent spaces to the
manifold. More discussion is given in Section 8 regarding the potential
application of exponential manifolds to the problems discussed in this
manuscript.

Another area which uses maps from measures into Hilbert spaces are kernel mean
embeddings (Muandet et al., 2017; Gretton et al., 2012; Berlinet and Thomas-
Agnan, 2004; Guilbart, 1978). This method is discussed at length in Section 5.
The idea is to represent a measure using a kernel and then reason with the
representation of the measure rather than the measure itself, for example
forming distances and performing hypothesis testing. This method has also been
applied to Bayesian computation (Fukumizu et al., 2013). The crucial
difference between this approach and the Bayes Hilbert space and exponential
manifold approach is that the map from measures to the Hilbert space is not
easy to invert, both theoretically and practically. This makes the kernel mean
embedding approach difficult to apply to the posterior approximation. The
kernel mean embedding distance will appear in Theorem 4 as a way to relate
distances in the Bayes Hilbert space to distances between the empirical
measures of observed data in posterior approximation.

##  3 Measure Discrepancies and Bayes Hilbert Spaces

The aim of this section is to prove a novel result which shows the Bayes
Hilbert space norm between two measures upper bounds typical notions of
discrepancy between measures. Specifically, Theorem 3 provides bounds for the
Hellinger, Kullback-Leibler and Wasserstein-$1$ distances in terms of a Bayes
Hilbert space norm. Such bounds are important if the Bayes Hilbert space is to
be used for posterior approximation as they show that the approximation is
meaningful. The bounds are obtained by a simple application of a recent result
from the Bayesian inverse problems literature regarding the stability of
posterior measures (Sprungk, 2020). This result is possible since the notion
of error of posterior approximation employed within Bayesian inverse problem
literature coincides with the distance in Bayes Hilbert spaces.

Before the result is stated some notation needs to be outlined. First,
$\Theta$ is assumed to be a separable, complete metric space with metric
$d_{\Theta}$ and let $\mathcal{B}(\Theta)$ denote the set of finite Borel
measures on $\Theta$ and $\mathcal{P}(\Theta)\subset\mathcal{B}(\Theta)$ the
subset of probability measures. In all three following definitions,
$\eta,\nu\in\mathcal{B}(\Theta)$.

Assuming $\eta,\nu$ are both absolutely continuous with respect to some
measure $\lambda$, such a $\lambda$ always exists for example
$\lambda=\eta+\nu$, then the square of the Hellinger distance is defined

$\displaystyle\text{H}(\eta,\nu)^{2}=\frac{1}{2}\int_{\Theta}\left(\sqrt{\frac{\mathrm{d}\eta}{\mathrm{d}\lambda}}-\sqrt{\frac{\mathrm{d}\nu}{\mathrm{d}\lambda}}\right)^{2}\mathrm{d}\lambda,$
(7)

where this definition is invariant to the choice of $\lambda$.

Assuming $\eta$ is absolutely continuous with respect to $\nu$ the Kullback-
Leibler (KL) divergence is defined

$\displaystyle\text{KL}(\eta\;\|\;\nu)=\int_{\Theta}\log\frac{\mathrm{d}\eta}{\mathrm{d}\nu}\mathrm{d}\eta,$

for more discussion on the KL divergence see Section 7.

Let $\text{Lip}(1)$ denote the set of functions
$f\colon\Theta\rightarrow\mathbb{R}$ with Lipschitz constant less than or
equal to one, then the Wasserstein-$1$ distance is defined

$\displaystyle\text{W}_{1}(\eta,\nu)=\sup_{f\in\text{Lip}(1)}\left\lvert\int_{\Theta}f\mathrm{d}\eta-\int_{\Theta}f\mathrm{d}\nu\right\rvert.$

For $\eta\in\mathcal{P}(\Theta)$ define

$\displaystyle\lVert\eta\rVert_{\mathcal{P}^{2}}$
$\displaystyle=\inf_{\theta_{0}\in\Theta}\left(\int_{\Theta}d_{\Theta}(\theta,\theta_{0})^{2}\mathrm{d}\eta(\theta)\right)^{1/2}$
$\displaystyle\mathcal{P}^{2}(\Theta)$
$\displaystyle=\\{\eta\in\mathcal{P}(\Theta)\colon\lVert\eta\rVert_{\mathcal{B}^{2}}<\infty\\},$

which is a quantity which will occur in the Wasserstein-$1$ bound that
provides a notion of size of the sample space with respect to the measure and
the metric on the space. Finally, for $f\in L^{2}(\mu)$,
$\operatorname*{ess\,sup}_{\mu}f$ is the essential supremum of $f$ and if
there is an $\tilde{f}$ in the $L^{2}(\mu)$-equivalence class of $f$ such that
$\tilde{f}(\theta)\leq B\>\forall\theta\in\Theta$ then
$\operatorname*{ess\,sup}_{\mu}f\leq B$. For a more detailed description of
essential supremum consult Bogachev (2007, Section 2.11).

######  Theorem 3.

Let $\Theta$ be a complete, separable metric space and
$\mu\in\mathcal{P}^{2}(\Theta),\eta,\nu\in\mathcal{B}(\Theta)$. Let
$\eta,\nu\in B^{2}(\mu)$ with

$\displaystyle\frac{\mathrm{d}\eta}{\mathrm{d}\mu}$
$\displaystyle=Z_{\eta}^{-1}\exp(\Psi_{\mu}(\eta))\qquad
Z_{\eta}=\mathbb{E}_{\mu}[\exp(\Psi_{\mu}(\eta))]$ (8)
$\displaystyle\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$
$\displaystyle=Z_{\nu}^{-1}\exp(\Psi_{\mu}(\nu))\qquad
Z_{\nu}=\mathbb{E}_{\mu}[\exp(\Psi_{\mu}(\nu))].$ (9)

Assume that
$\operatorname*{ess\,sup}_{\mu}\Psi_{\mu}(\eta),\operatorname*{ess\,sup}_{\mu}\Psi_{\mu}(\nu)\leq
B$ then

$\displaystyle\emph{H}(\eta,\nu)$
$\displaystyle\leq\frac{1}{2}\left(e^{B}+e^{2B}\right)^{1/2}\lVert\eta-\nu\rVert_{B^{2}(\mu)}$
$\displaystyle\emph{KL}(\eta\;\|\;\nu)$ $\displaystyle\leq
2e^{B}\lVert\eta-\nu\rVert_{B^{2}(\mu)}$ $\displaystyle\emph{W}_{1}(\eta,\nu)$
$\displaystyle\leq(e^{B}+e^{2B})\lVert\mu\rVert_{\mathcal{P}^{2}}\lVert\eta-\nu\rVert_{B^{2}(\mu)}.$

The proof is contained in the Appendix and is a straightforward application of
the results in Sprungk (2020). A bound using an $L^{1}(\mu)$ norm also holds
for the KL and Wasserstein-$1$ case and $\mu$ being a probability measure can
be relaxed to $\mu$ being a finite measure as the cost of an extra constant,
as discussed in the proof. So long as the upper bound $B$ is controlled,
Theorem 3 shows that the Bayes Hilbert space norm converging to zero implies
that the other three notions of distance converge to zero too. The most
similar result is Huggins et al. (2018, Proposition 6.3) which provides a
bound on the Wasserstein-$1$ and Wasserstein-$2$ distances, the difference is
that an $L^{2}$ norm is used which involves derivatives of the log-likelihoods
rather than the Bayes Hilbert space norm. Corollary 1 in Section 4 deals with
the specific case of approximating a posterior measure.

Implicit in all of this is the choice of the base measure $\mu$. The choice
does not impact the left hand side of the bounds and will impact the right
hand side by changing $B$ and the value of the norm. Theorem 3 provides a
guiding principal of choosing $\mu$ to make $B$ as small as possible.

######  Example 3.

Continuing the logistic regression example, since $l_{x_{n}}\leq 1$ the log-
likelihood is bounded as $\log l_{x_{n}}\leq 0$ therefore the CLR is bounded
as

$\displaystyle\Psi_{\pi_{0}}(\pi)\leq-\sum_{n=1}^{N}\mathbb{E}_{\pi_{0}}\left[\log\left(\frac{y_{n}}{1+e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}+\frac{(1-y_{n})e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}{1+e^{-\langle\theta,u_{n}\rangle_{\mathbb{R}^{d}}}}\right)\right],$

which can be used as a bound in Theorem 3.

The notation for the $Z_{\eta},Z_{\nu}$ terms in the Radon-Nikodym derivatives
in (8),(9) is atypical since they do not factor out all the scalar factors.
Specifically, since the CLR map $\Psi_{\mu}$ involves centering by a constant,
see (6), there is still a multiplicative factor which could be absorbed into
the normalising terms. This was done when introducing the posterior in (4)
with the $Z$ term. The reason this is not done in Theorem 3 is to retain the
full expressions for Bayes Hilbert space norm in the right hand side.

######  Example 4.

Continuing the logistic regression example, the $Z$ term in (4) would be

$\displaystyle
Z=\mathbb{E}_{\pi_{0}}\left[L\right]=\mathbb{E}_{\pi_{0}}\left[\exp\left(\mathcal{L}\right)\right]$

where $L=\prod_{n=1}^{N}l_{x_{n}},\mathcal{L}=\log L$. In contrast to this the
$Z_{\pi}$ term in Theorem 3 is

$\displaystyle
Z_{\pi}=\mathbb{E}_{\pi_{0}}\left[\exp(\Psi_{\pi_{0}}(\pi))\right]=\mathbb{E}_{\pi_{0}}\left[\exp(\mathcal{L}-\mathbb{E}_{\pi_{0}}[\mathcal{L}])\right].$

Typically, bounds of the type used in Theorem 3, see Sprungk (2020), assume
that $\mathcal{L}\leq 0$ and involve $Z^{-1}$ as a constant, where
$Z=\mathbb{E}_{\pi_{0}}[L]$. Note $\mathcal{L}\leq 0$ is a very mild
assumption and can be obtained when assuming first $\mathcal{L}\leq C$ for
constant $C$ and then using $\mathcal{L}-C$, see Sprungk (2020, Section 2).
The task of upper bounding $Z^{-1}$ is closely related to the upper bound
assumption on $\Psi_{\pi_{0}}$ in Theorem 3. To see this note that
$Z^{-1}=\mathbb{E}_{\pi_{0}}[\exp(\mathcal{L})]^{-1}\leq\exp(-\mathbb{E}_{\pi_{0}}[\mathcal{L}])$
by Jensen’s inequality. When $\mathcal{L}\leq 0$ then the bound $B$ only needs
to satisfy $-\mathbb{E}_{\pi_{0}}[\mathcal{L}]\leq B$. Hence $Z^{-1}\leq
e^{B}$. This shows how despite using different normalising constants in the
Radon-Nikodym derivative representation, to ensure the bounds still involve
the CLR transform, Theorem 3 recovers constants in the upper bounds typical of
results in the Bayesian inverse problems literature, for example Sprungk
(2020, Theorem 5), with the difference being the slack in Jensen’s inequality.

##  4 Bayesian Coresets and Bayes Hilbert Spaces

The aim of this section will be to introduce Bayesian coreset posterior
approximations (Campbell and Broderick, 2019; Huggins et al., 2016) and
provide a novel relationship between them and Bayes Hilbert spaces. This
connection is crucial for the rest of this manuscript and will be analysed
further in later sections.

As mentioned in the previous section, for a likelihood that is a product of
$N$ terms, the typical cost per iteration of a sampling algorithm is $O(N)$.
This is because typical sampling algorithms have to “touch” each data point in
the likelihood during each of their iterations. The idea of a Bayesian coreset
is to approximate the likelihood with a weighted product of $M$ terms. The
resulting approximation will have $O(M)$ cost which can be substantially
cheaper if $M\ll N$. The first paper on Bayesian coresets was focused on
logistic regression (Huggins et al., 2016). Since then the area has developed
widely and has produced multiple innovations, both in terms of how the
distance to the full posterior is evaluated and in the optimisation methods
used to form the coreset (Campbell and Broderick, 2018, 2019; Manousakas et
al., 2020; Zhang et al., 2021; Naik et al., 2022).

Retaining the same setting as Section 2.2, the prior is $\pi_{0}$, the
observed data $\\{x_{n}\\}_{n=1}^{N}\subset\mathcal{X}$, the individual
likelihood functions $l_{x_{n}}(\theta)\coloneqq l(\theta,x_{n})$ and the full
likelihood $L=\prod_{n=1}^{N}l_{x_{n}}$. For the rest of this manuscript, a
Bayesian coreset is a collection of points
$\\{z_{m}\\}_{m=1}^{M}\subset\mathcal{X}$ and non-negative weights
$\\{w_{m}\\}_{m=1}^{M}\subset\mathbb{R}_{\geq 0}$ for some $M\in\mathbb{N}$
and the Bayesian coreset posterior approximation corresponding to the Bayesian
coreset is the measure $\pi_{w,z}$ satisfying

$\displaystyle\frac{\mathrm{d}\pi_{w,z}}{\mathrm{d}\pi_{0}}=Z_{w,z}^{-1}L_{w,z},$
(10)

where $L_{w,z}=\prod_{m=1}^{M}l_{z_{m}}^{w_{m}}$ and
$Z_{w,z}=\mathbb{E}_{\pi_{0}}[L_{w,z}]$. Non-negative weights are used so that
the resulting posterior is valid in the sense that it has finite measure. This
approximation keeps the prior the same and approximates the posterior by using
the likelihood $L_{w,z}$ instead of $L$. Note that if $M=N$ and
$w_{m}=1,z_{m}=w_{m}\>\forall m\>$ then $L=L_{w,z}$ and if $w_{m}=0$ then
$l_{z_{m}}$ plays no role in the approximation.

The idea behind this approximation is that some of the observed data points
may be redundant and do not need to be included, or they may be accounted for
by re-weighting a likelihood term based on a different data point. The fact
that $\\{z_{m}\\}_{m=1}^{M}$ does not have to be a subset of the observed data
$\\{x_{n}\\}_{n=1}^{N}$ facilitates a flexible and efficient approximation.
The canonical example in this case is a Gaussian mean inference task. In this
scenario the posterior can be written exactly with $M=1$, see Manousakas et
al. (2020, Section 3). The scenario where
$\\{z_{m}\\}_{m=1}^{M}\not\subset\\{x_{n}\\}_{n=1}^{N}$ was named a pseudo-
coreset in Manousakas et al. (2020). This terminology is not maintained in the
present discussion since the general case, where no assumptions are placed on
$\\{z_{m}\\}_{m=1}^{M}$ so it could or could not be a subset of
$\\{x_{n}\\}_{n=1}^{N}$, is investigated.

The CLR transform of $\pi_{w,z}$ is easily obtained by using (6). Let $\mu$ be
the base measure and assume that $\pi_{w,z},\pi_{0}\in B^{2}(\mu)$, then

$\displaystyle\Psi_{\mu}(\pi_{w,z})=\mathcal{L}_{w,z}-\mathbb{E}_{\mu}\left[\mathcal{L}_{w,z}\right]+(\log
p_{\pi_{0},\mu}-\mathbb{E}_{\mu}[\log p_{\pi_{0},\mu}]),$ (11)

where $\mathcal{L}_{w,z}\coloneqq\log L_{w,z}=\sum_{m=1}^{M}w_{m}\log
l_{z_{m}}$ is the weighted log-likelihood based on $\\{z_{m}\\}_{m=1}^{M}$.

The Bayesian coreset approximation is leveraging the structure inherent in the
Bayesian methodology that the Radon-Nikodym derivative of the posterior with
respect to the prior is equal, up to a scalar constant, to a product of
likelihoods. This means the full posterior can be recovered exactly with the
choice $M=N$ and $w_{m}=1,z_{m}=w_{m}\>\forall m\>$. Therefore, exact recovery
only involves a finite number of parameter choices. This is different to
typical function approximation methods where the target function can only be
recovered exactly with an infinite number of parameter choices, for example an
infinite basis expansion.

Some simple choices of $\mu$ in the Bayesian coreset scenario are are
$\mu=\pi_{0}$, setting $\mu$ to be the posterior based upon a random subset of
the data points and a Laplace approximation to the posterior (Campbell and
Broderick, 2019). It is currently an open question as to how to chose $\mu$
and what it means for $\mu$ to be a best choice. Intuitively one hopes that
$\mu$ has most mass where the posterior has most mass but this causes a
chicken and the egg problem as the entire task is to approximate the
posterior. as mentioned after Theorem 3 a guide can be to choose a $\mu$ which
minimises the bound $B$ using in Theorem 3.

The expressions for $\pi,\pi_{w,z}$ can be inserted into Theorem 3 to
understand how the Bayes Hilbert space norm implies a bound on the common
notion of distance for probability measures in this particular case. A typical
scenario for Bayesian coresets is having $M=N$ and $z_{n}=x_{n}$ and only
changing the weights, with some weights being zero to provide sparsity. The
next lemma applies Theorem 3 to this case.

######  Lemma 1.

Under the assumptions of Theorem 3, if $\mathcal{L}\leq C$ for some constant
$C$, $M=N$, $z_{n}=x_{n}\>\forall n$, and the weights $w_{n}$ are non-negative
with $\lVert w\rVert_{\infty}\leq W$ then $B$ can be set as
$B=-\max(1,W)\mathbb{E}_{\mu}[\mathcal{L}]+C$.

###### Proof.

Under the assumptions,
$\Psi_{\mu}(\pi)=\mathcal{L}-C-\mathbb{E}_{\mu}[\mathcal{L}]+C\leq-\mathbb{E}_{\mu}[\mathcal{L}]+C$
with the analogous argument for $\Psi_{\mu}(\pi_{w,z})$. Using the assumptions
on $w_{m},z_{m}$

$\displaystyle-\mathbb{E}_{\mu}[\mathcal{L}_{w,z}]+C=\sum_{n=1}^{N}-w_{m}\mathbb{E}_{\mu}[l_{z_{m}}]+C$
$\displaystyle\leq\max(1,W)\sum_{n=1}^{N}-\mathbb{E}_{\mu}[l_{x_{n}}]+C$
$\displaystyle=-\max(1,W)\mathbb{E}_{\mu}[\mathcal{L}]+C,$

which completes the proof. ∎

Lemma 1 shows that if one forms a Bayesian coreset using points that are a
subset of the observed data points, which is very common (Campbell and
Broderick, 2018, 2019; Naik et al., 2022), then as long as the weights are
bounded the term $B$ will be agnostic of the choice of weights. The
requirement that the weights are bounded is not typically explicitly enforced
in existing Bayesian coreset algorithms although it is reasonable. This is
because the norm of the full likelihood is
$\lVert\Psi_{\mu}(\pi)\rVert_{L^{2}(\mu)}$ so one would expect each weight to
not be much larger than this overall norm value. An exception where the
weights are bounded is in the importance sampling coreset methods (Huggins et
al., 2016).

Armed with the approximation scheme (10) and the ability to bound common
notions of discrepancy in terms of the Bayes Hilbert space norm using
Corollary 1 and Theorem 3, the question is now how to choose the weights and
the points in the Bayesian coreset approximation. Many methods have emerged in
the literature over the past five years or so and the subsequent sections will
be dedicated to viewing these through the lens of a Bayes Hilbert space. The
first step in this process begins in the next section which makes a novel
connection between the Bayes Hilbert space distance and maximum mean
discrepancy distance, a kernel-based distance for measures.

##  5 Maximum Mean Discrepancy and Bayes Hilbert Spaces

This aim of this section is to outline a novel relationship between maximum
mean discrepancy (MMD) and Bayesian coresets by relating the MMD distance to a
Bayes Hilbert space distance. The relationship is made via the kernel pre-
image problem (Schölkopf and Smola, 2002, Chapter 18). This relationship will
later be used in Section 6 to investigate a genre of Bayesian coreset methods
called Hilbert Bayesian coresets. The rest of this section will introduce MMD
then the kernel pre-image problem and finally make the relation to Bayesian
coresets.

MMD is a kernel-based discrepancy between measures and has been studied for
nearly two decades in statistical machine learning (Gretton et al., 2006,
2012; Muandet et al., 2017) and computational statistics (Teymur et al., 2021;
Dwivedi and Mackey, 2021; Mak and Joseph, 2018), with its origins being in the
70s in more abstract mathematical statistics under a different name (Guilbart,
1978). MMD has both theoretical and practical strengths. The theoretical
strengths include the elegant representation in terms of a reproducing kernel
Hilbert space (RKHS), explained below, which facilitates MMD being rewritten
in multiple helpful ways. The main practical strength of MMD is that it can be
easily estimated by simply evaluating the kernel on samples from the two
measures in question.

To begin the description of MMD one must first start with a kernel. A kernel
$k\colon\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ is a symmetric,
positive definite function. The point of a kernel is that it measures
similarity between its inputs. A simple way to construct a kernel is to take
any function $\phi\colon\mathcal{X}\rightarrow H$, which will be called a
feature map, where $H$ is some real Hilbert space, called the feature space,
then set $k(x,y)=\langle\phi(x),\phi(y)\rangle_{H}$. In this case the kernel
is measuring the similarity between $x$ and $y$ by measuring the similarity of
the features $\phi(x)$ and $\phi(y)$ using the inner product in the feature
space.

For any kernel $k$ there is a space of functions called the reproducing kernel
Hilbert space that is uniquely determined by the kernel which has specific
properties relating to the kernel. In particular, the RKHS, denoted
$\mathcal{H}_{k}$, is the unique Hilbert space of functions from $\mathcal{X}$
to $\mathbb{R}$ such that
$k(x,\cdot)\in\mathcal{H}_{k}\>\forall\>x\in\mathcal{X}$ and $\langle
f,k(x,\cdot)\rangle_{k}=f(x)\>\forall f\in\mathcal{H}_{k},x\in\mathcal{X}$
where $\langle\cdot,\cdot\rangle_{k}$ is the inner product of
$\mathcal{H}_{k}$, see Christmann and Steinwart (2008, Section 4.2). The
latter property is called the reproducing property and it is extremely useful
since, for functions in $\mathcal{H}_{k}$, it facilitates many quantities of
interest to be written purely in terms of the kernel, which is easy to
evaluate. For more details regarding kernels and their RKHS consult
(Christmann and Steinwart, 2008; Berlinet and Thomas-Agnan, 2004; Paulsen and
Raghupathi, 2016).

Before outlining MMD a point on notation. Latin letters will be used to denote
measures and random variables on $\mathcal{X}$, the space where observed data
for the Bayesian inference problem will lie, as is common in the MMD
literature (Sriperumbudur et al., 2010). This is to distinguish them from
measures on $\Theta$, the parameter space, which use Greek letters as was done
in Section 2.1.

MMD is now defined, for a more in depth discussion than what is presented here
consult Sriperumbudur et al. (2010). Let
$\mathcal{B}_{k}(\mathcal{X})=\\{P\in\mathcal{B}(\mathcal{X})\colon\mathbb{E}_{P}[\sqrt{k(X,X)}]<\infty\\}$
where in the expectation $X$ denotes the random variable in $\mathcal{X}$ with
law $P$ and $\mathcal{B}(\mathcal{X})$ is the set of finite Borel measures on
$\mathcal{X}$. Note that all finite sums of empirical measures are in
$\mathcal{B}_{k}$. For $P,Q\in\mathcal{B}_{k}(\mathcal{X})$ the MMD using
kernel $k$ is

$\displaystyle\text{MMD}_{k}(P,Q)\coloneqq\sup_{\lVert f\rVert_{k}\leq
1}\left\lvert\mathbb{E}_{P}[f]-\mathbb{E}_{Q}[f]\right\rvert,$ (12)

where the supremum is being taken over the unit ball of the RKHS. The key
property of MMD is that it may be re-written in different forms to make it
more interpretable and easily estimated.

The first of these involves kernel mean embeddings. This is a mapping of a
measure into the RKHS via the kernel. For a kernel $k$ the kernel mean
embedding
$\Phi_{k}\colon\mathcal{B}_{k}(\mathcal{X})\rightarrow\mathcal{H}_{k}$ is
defined

$\displaystyle\Phi_{k}(P)\coloneqq\int_{\mathcal{X}}k(x,\cdot)\mathrm{d}P(x).$

The idea of a kernel mean embedding is to provide a feature expansion in the
RKHS of a measure $P$, much like how standard kernel methods on
$\mathbb{R}^{d}$ involve feature expansions of the data.

Analogous to how the reproducing property represents pointwise evaluation of
functions in the RKHS, kernel mean embeddings represent integration of
functions in the RKHS. For $f\in\mathcal{H}_{k}$ and
$P\in\mathcal{B}_{k}(\mathcal{X})$

$\displaystyle\mathbb{E}_{P}[f]=\int_{\mathcal{X}}f(x)\mathrm{d}P(x)=\int_{\mathcal{X}}\langle
f,k(x,\cdot)\rangle_{k}\mathrm{d}P(x)=\langle
f,\int_{\mathcal{X}}k(x,\cdot)\mathrm{d}P(x)\rangle_{k}=\langle
f,\Phi_{k}(P)\rangle_{k},$ (13)

where the swapping of the integral and inner product is a result of the
integrability assumption on $k$ and $P$ made when assuming
$P\in\mathcal{B}_{k}(\mathcal{X})$ (Sriperumbudur et al., 2010, Theorem 1).

Kernel mean embeddings relate to MMD in the following elegant way

$\displaystyle\text{MMD}_{k}(P,Q)$ $\displaystyle=\sup_{\lVert f\rVert_{k}\leq
1}\left\lvert\mathbb{E}_{P}[f]-\mathbb{E}_{Q}[f]\right\rvert$
$\displaystyle=\sup_{\lVert f\rVert_{k}\leq 1}\left\lvert\langle
f,\Phi_{k}(P)-\Phi_{k}(Q)\rangle_{k}\right\rvert$ (14)
$\displaystyle=\lVert\Phi_{k}(P)-\Phi_{k}(Q)\rVert_{k},$ (15)

where (14) is by (13) and (15) is by the Cauchy-Schwarz theorem. This shows
that MMD is equal the the difference in the RKHS of the feature expansions of
the two measures.

A third representation of MMD can be obtained by starting with (15) and again
using (13)

$\displaystyle\text{MMD}_{k}(P,Q)^{2}$
$\displaystyle=\lVert\Phi_{k}(P)-\Phi_{k}(Q)\rVert_{k}^{2}$
$\displaystyle=\langle\Phi_{k}(P),\Phi_{k}(P)\rangle_{k}-2\langle\Phi_{k}(P),\Phi_{k}(Q)\rangle_{k}+\langle\Phi_{k}(Q),\Phi_{k}(Q)\rangle_{k}$
(16) $\displaystyle=\mathbb{E}_{P\times
P}[k(X,X^{\prime})]-2\mathbb{E}_{P\times Q}[k(X,Y)]+\mathbb{E}_{Q\times
Q}[k(Y,Y^{\prime})],$ (17)

where (16) is simply expanding the RKHS norm and (17) is by using (13). In the
derivation above $X$ has law $P$, $Y$ has law $Q$ and $X^{\prime},Y^{\prime}$
are i.i.d. copies of $X,Y$, respectively. This double expectation formula is
what makes MMD practical since it can be easily estimated empirically given
samples from $P,Q$ (Gretton et al., 2012).

MMD has found many uses in machine learning and computational statistics, for
example two-sample testing (Gretton et al., 2012), parameter inference
(Cherief-Abdellatif and Alquier, 2020), training generative models (Li et al.,
2017) and distribution compression (Dwivedi and Mackey, 2021). A long list of
other applications and references may be found in Muandet et al. (2017,
Section 3.5).

The use for MMD that is the focus of this section is the following pre-image
problem (Schölkopf and Smola, 2002, Chapter 18).

######  Problem 1 (The pre-image problem).

Given a kernel $k$, $\\{x_{n}\\}_{n=1}^{N}\subset\mathcal{X}$ and
$M\in\mathbb{N}$ find non-negative weights
$\\{w_{m}\\}_{m=1}^{M}\subset\mathbb{R}_{\geq 0}$ and
$\\{z_{m}\\}_{m=1}^{M}\subset\mathcal{X}$ which minimise

$\displaystyle\emph{MMD}_{k}(P_{w,z},P_{x})^{2}$
$\displaystyle=\left\lVert\sum_{m=1}^{M}w_{m}k(z_{m},\cdot)-\sum_{n=1}^{N}k(x_{n},\cdot)\right\rVert_{k}^{2}$
$\displaystyle=\sum_{m,m^{\prime}=1}^{M}w_{m}w_{m^{\prime}}k(z_{m},z_{m^{\prime}})-2\sum_{m=1}^{M}\sum_{n=1}^{N}w_{m}k(z_{m},x_{n})+\sum_{n,n^{\prime}=1}^{N}k(x_{n},x_{n^{\prime}}),$

where
$P_{w,z}=\sum_{m=1}^{M}w_{m}\delta_{z_{m}},P_{x}=\sum_{n=1}^{N}\delta_{x_{n}}$
are the empirical measures corresponding to the two sets of data.

The pre-image problem was studied heavily in the 90s and early 00s in the
context of support vector machines (Burges, 1996; Scholkopf et al., 1999) and
kernel principal component analysis (Kwok and Tsang, 2004; Mika et al., 1998).
It is called a pre-image problem since it is seeking to find weights and
points which map close to $\sum_{n=1}^{N}k(x_{n},\cdot)$ in the RKHS.

This problem relates to many topics within kernel methods and computational
statistics. Reduced set methods solve a similar pre-image problem to find
sparse representations of kernel-based algorithms (Burges, 1996). Kernel
principal component analysis also solves a similar pre-image problem
(Schölkopf et al., 1997). Kernel herding (Chen et al., 2010; Bach et al.,
2012) is a method of solving the pre-image problem by greedily picking points
and then choosing weights. The main two methods either use uniform weights or
line-search. More will be discussed about kernel herding in Section 6. The
scenario where $\\{z_{m}\\}_{m=1}^{M}\subset\\{x_{n}\\}_{n=1}^{N}$ is known as
distribution compression (Dwivedi and Mackey, 2021), or quantisation (Teymur
et al., 2021; Graf and Luschgy, 2000), where the weights are typically left as
uniform over the data points. The scenario where a user is not trying to
approximate an empirical distribution but rather a continuous distribution
using kernel-based approaches is covered within the kernel Stein discrepancy
literature (Riabiz et al., 2022; Anastasiou et al., 2023) with quasi-Monte
Carlo being a related field which aims to find discretisations of given
measures for integral approximation (Caflisch, 1998; Dick et al., 2013). In
statistical depth, the $h$-depth (Wynne and Nagy, 2021) is an instance of the
pre-image problem when $M=1$ so only one point is used to represent the data.

The pre-image problem immediately appears to be related to Bayesian coresets
since both involve starting with a data set and finding weights and points.
The subtlety is that the effectiveness of a Bayesian coreset is measured by
the quality of its corresponding posterior approximation. This means the
notion of quality of weights and points lying in $\mathcal{X}$ is expressed in
terms of a distance between measures on $\Theta$. This is in contrast to most
investigations of the pre-image problem and the related problems outlined
above. The difference is due to the common use of kernels which involve
expressions that measure distance between their inputs purely in terms of the
geometry of $\mathcal{X}$. The trick to link the pre-image problem (1) to
Bayesian coresets is to use a kernel which maps the data into a Bayes Hilbert
space. This will facilitate comparison of points in $\mathcal{X}$ in terms of
corresponding posteriors.

This is done by using a kernel of the form
$k(x,y)=\langle\phi_{\mu}(x),\phi_{\mu}(y)\rangle_{H}$, such kernels were
discussed at the start of this section. The same setting and notation for
Bayes’ theorem used in Section 2.2 is maintained. The likelihood function
based at a point $x\in\mathcal{X}$ is $l_{x}$. For some base measure $\mu$ the
feature map is $\phi_{\mu}(x)=\log l_{x}-\mathbb{E}_{\mu}[\log l_{x}]$ with
feature space $H=L^{2}(\mu)$ and the kernel is defined

$\displaystyle k(x,y)=\langle\phi_{\mu}(x),\phi_{\mu}(y)\rangle_{L^{2}(\mu)}.$
(18)

This kernel measures the similarity between two points in $\mathcal{X}$ by
comparing the similarity of the centred versions of the log-likelihoods based
at those points. This provides the crucial link between the data space and the
sample space.

######  Theorem 4.

Let $\mathcal{X}$ be a metric space equipped with its Borel $\sigma$-algebra.
Let $\mu$ be a finite measure on the measurable space $\Theta$ and $l$ a
likelihood function such that $\phi_{\mu}\colon\mathcal{X}\rightarrow
L^{2}(\mu)$, $\phi_{\mu}(x)=\log l_{x}-\mathbb{E}_{\mu}[\log l_{x}]$ is well-
defined and measurable. Let $\pi_{0}\in B^{2}(\mu)$ and $\pi\in B^{2}(\mu)$ be
the posterior with prior $\pi_{0}$, likelihood $l$ and observations
$\\{x_{n}\\}_{n=1}^{N}\subset\mathcal{X}$. Let $M\in\mathbb{N}$,
$\\{w_{m}\\}_{m=1}^{M}\subset\mathbb{R}$ and
$\\{z_{m}\\}_{m=1}^{M}\subset\mathcal{X}$. If $k$ is the kernel (18) then

$\displaystyle\emph{MMD}_{k}$
$\displaystyle(P_{w,z},P_{x})=\left\lVert\pi_{w,z}-\pi\right\rVert_{B^{2}(\mu)},$
(19)

where $\pi_{w,z}$ is the posterior based on the weighted likelihood (10).

###### Proof.

The assumptions ensure that $k$ is well-defined and measurable and that
empirical measures on $\mathcal{X}$ are measurable.

$\displaystyle\text{MMD}_{k}$
$\displaystyle(P_{w,z},P_{x})^{2}=\left\lVert\sum_{m=1}^{M}w_{m}k(z_{m},\cdot)-\sum_{n=1}^{N}k(x_{n},\cdot)\right\rVert_{k}^{2}$
(20)
$\displaystyle=\sum_{m,m^{\prime}=1}^{M}w_{m}w_{m^{\prime}}k(z_{m},z_{m^{\prime}})-2\sum_{m=1}^{M}\sum_{n=1}^{N}w_{m}k(z_{m},x_{n})+\sum_{n,n^{\prime}=1}^{N}k(x_{n},x_{n^{\prime}})$
(21)
$\displaystyle=\left\lVert\sum_{m=1}^{M}\phi_{\mu}(z_{m})-\sum_{n=1}^{N}\phi_{\mu}(x_{n})\right\rVert_{L^{2}(\mu)}^{2}$
(22)
$\displaystyle=\left\lVert\Psi_{\mu}(\pi_{w,z})-\Psi_{\mu}(\pi)\right\rVert_{L^{2}}^{2}(\mu)$
(23) $\displaystyle=\left\lVert\pi_{w,z}-\pi\right\rVert_{B^{2}(\mu)}^{2},$
(24)

where (20) is by (15), (21) is by the reproducing property, (22) is by the
definition of $k$, (23) is by the expression for the CLR transform of the
approximation posterior (11) and the full posterior (6) where the term in the
brackets in the expressions cancel out and finally (24) is by the definition
of the Bayes Hilbert space norm (3). ∎

Theorem 4 shows that the pre-image problem, Problem 1, is equivalent to
minimising a Bayes Hilbert space distance between the Bayesian coreset
posterior and the target posterior when using the kernel (18). An immediate
consequence of this result is that bounds on MMD can be translated into bounds
on the Bayes Hilbert space distance, and therefore by Theorem 3 bounds on
commonly used distances between measures. An simple example is the following
corollary which deals with the case of approximating a posterior by using a
likelihood based on a subset of the observed data uniformly sampled without
replacement. This random uniform subsampling without replacement is used as a
common basic benchmark for evaluating approximate posterior methods (Campbell
and Broderick, 2019).

######  Corollary 1.

Under the same assumptions as Theorem 4, let $M\leq N$ and set $\pi_{M}$ to be
the approximate posterior satisfying
$p_{\pi_{M},\pi_{0}}=_{B}(N/M)\odot\bigoplus_{m=1}^{M}l_{z_{m}}$ where
$\\{z_{m}\\}_{m=1}^{M}$ is randomly uniformly subsampled without replacement
from $\\{x_{n}\\}_{n=1}^{N}$. Assume that
$\lVert\phi_{\mu}(x_{n})\rVert_{L^{2}(\mu)}\leq\gamma\>\forall\>1\leq n\leq N$
then with probability at least $1-\delta$

$\displaystyle\left\lVert\pi-\pi_{M}\right\rVert_{B^{2}(\mu)}\leq\sqrt{\frac{8\gamma^{2}N(N-M+1)}{M}}\sqrt{2\log\left(\frac{2}{\delta}\right)}.$

The proof is a direct application of the kernel mean embeddings concentration
inequality of Schneider (2016). The bound is similar to other concentration
inequalities for Bayesian coresets, such as Campbell and Broderick (2019,
Theorem 4.1). The difference is that Corollary 1 focuses on the simple random
sub-sample case. More sophisticated sub-sampling methods, such as importance
sampling as was done in Campbell and Broderick (2019, Theorem 4.1) could also
be easily applied in the context of a Bayes Hilbert space due to the
Hilbertian structure.

The connection between MMD and Bayes Hilbert spaces will be used again in the
next section to outline a novel connection between Bayesian coreset algorithms
and kernel-based methods for solving the pre-image problem.

##  6 Hilbert Bayesian Coresets

The aim of this section is to provide a novel connection between Hilbert
Bayesian coreset algorithms and Bayes Hilbert spaces. The former are Bayesian
coreset algorithms which revolve around using a Hilbert norm and inner product
to define a notion of distance between the target posterior and the
approximation. The primary focus will be on showing that the Frank-Wolfe
Bayesian coreset algorithm (Campbell and Broderick, 2019) is equivalent to
Frank-Wolfe kernel-herding (Chen et al., 2010; Bach et al., 2012) when using
the kernel (18). This novel relationship will be shown to be a consequence of
Theorem 4. A secondary focus will be on the iterative hard thresholding
coreset algorithm (Zhang et al., 2021) which will be related to MMD descent
methods.

###  6.1 Frank-Wolfe Bayesian coresets

The goal of both the Frank-Wolfe Bayesian coreset algorithm and Frank-Wolfe
kernel-herding is essentially is to approximate an element of a Hilbert space
using a candidate set. In both cases, as the name suggests, the Frank-Wolfe
optimisation method is used to find the best approximation over the set. To
begin with, the approximation problem of concern and the Frank-Wolfe algorithm
solution is outlined only in the generality needed to make the desired
connections. For further discussion of the method consult Clarkson (2010).

For some Hilbert space $H$, fix some approximation target $f\in H$ and let
$\\{g_{n}\\}_{n=1}^{N}\subset H$ be a set of elements which will be used to
approximate the target. Let $\\{\sigma_{n}\\}_{n=1}^{N}\subset\mathbb{R}_{>0}$
be a set of positive real numbers, set $\sigma=\sum_{n=1}^{N}\sigma_{n}$ and
define the polytope $\mathcal{W}=\\{w\in\mathbb{R}^{N}\colon w\geq
0,\sum_{n=1}^{N}w_{n}\sigma_{n}\leq\sigma\\}$, noting it has vertices
$v_{n}=(\sigma/\sigma_{n})\mathbf{1}_{n}$ where $\mathbf{1}_{n}$ is the one-
hot vector of all zeros except a $1$ in the $n$-th entry. The goal is to find

$\displaystyle\operatorname*{arg\,min}_{w\in\mathcal{W}}\frac{1}{2}\lVert
f-g_{w}\rVert_{H}^{2},$ (25)

where $g_{w}=\sum_{n=1}^{N}w_{n}g_{n}$. The Frank-Wolfe method to solve this
problem involves iterative, conditional gradient updates. It looks at the
residual of the current approximation, finds the direction most aligned with
it, then performs line search in that direction. Specifically, at the $t$-th
iteration one performs the update of the current guess $u_{t}\in\mathcal{W}$
as follows

$\displaystyle\overline{u}_{t}$
$\displaystyle=\operatorname*{arg\,max}_{v_{n}}\langle
f-g_{u_{t}},g_{v_{n}}\rangle_{H}$ (26) $\displaystyle u_{t+1}$
$\displaystyle=(1-\rho_{t})u_{t}+\rho_{t}\overline{u}_{t},$

where $\rho_{t}$ is a scalar value calculated in closed form using
$\overline{u}_{t},u_{t},f$ (Bach et al., 2012) and
$v_{n}=(\sigma/\sigma_{n})\mathbf{1}_{n}$ are the vertices of $\mathcal{W}$,
meaning $g_{v_{n}}=(\sigma/\sigma_{n})g_{n}$. Only the vertices are optimised
over in (26) as this is equivalent to optimising over the whole polytope
$\mathcal{W}$ since the function to optimise in (26) is linear. This means
that only one member of the approximating dictionary is added at each
iteration, meaning after $T$ iterations the current solution is $T$-sparse.

The Frank-Wolfe Hilbert coreset algorithm derived in Campbell and Broderick
(2019) fits this template. A version which includes log-likelihood centering
is now outlined, which was not used in Campbell and Broderick (2018, 2019)
though was later advocated by Campbell and Beronov (2019); Zhang et al. (2021)
since it enforces the desired scale invariance property for the likelihood
approximation. Set $H=L^{2}(\mu)$, $f=\Psi_{\mu}(\pi),g_{n}=\phi_{\mu}(x_{n})$
and $\sigma_{n}=\lVert\phi_{\mu}(x_{n})\rVert_{L^{2}(\mu)}$. With these
substitutions the goal becomes finding

$\displaystyle\operatorname*{arg\,min}_{w\in\mathcal{W}}\frac{1}{2}\left\lVert\Psi_{\mu}(\pi)-\sum_{n=1}^{N}w_{n}\phi_{\mu}(x_{n})\right\rVert_{L^{2}(\mu)}^{2}=\operatorname*{arg\,min}_{w\in\mathcal{W}}\frac{1}{2}\left\lVert\Psi_{\mu}(\pi)-\Psi_{\mu}(\pi_{w})\right\rVert_{L^{2}(\mu)}^{2},$
(27)

where $\pi_{w}$ denotes the Bayesian coreset approximation (10) with weights
$w$ and $z_{n}=x_{n}\>\forall\>n$. The Frank-Wolfe iterations are

$\displaystyle\overline{u}_{t}=\frac{\sigma}{\sigma_{x_{n}}}\mathbf{1}_{n}\text{
where }n$
$\displaystyle=\operatorname*{arg\,min}_{n\in[N]}\left\langle\Psi_{\mu}(\pi)-g_{u_{t}},\frac{\sigma}{\sigma_{x_{n}}}\phi_{\mu}(x_{n})\right\rangle_{L^{2}(\mu)}$
(28) $\displaystyle u_{t+1}$
$\displaystyle=(1-\rho_{t})u_{t}+\rho_{t}\overline{u}_{t},$

where the full expression for the update of $\overline{u}_{t}$, see (26), has
been written for clarity and $[N]=\\{1,\ldots,N\\}$. Note Campbell and
Broderick (2019, Equation 4.2) matches (25) and Campbell and Broderick (2019,
Equation 4.4) matches (26) once the centred log-likelihoods are substituted
for the log-likelihoods. As (27) and (28) are written in terms of CLR
transforms, the $L^{2}(\mu)$ inner products can be written in terms of the
$B^{2}(\mu)$ inner product, see (2). This shows that the Frank-Wolfe Hilbert
coreset algorithm can be written in the language of Bayes Hilbert spaces.

The Frank-Wolfe coreset algorithm will now be related to Frank-Wolfe kernel-
herding. Kernel-herding is where the Hilbert space $H$ is set to be an RKHS
and the algorithm is used in the context of solving the kernel pre-image
problem, see Problem 1. Let $k$ be the kernel (18), $H=\mathcal{H}_{k}$ with
approximation target $f=\sum_{n=1}^{N}k(x_{n},\cdot)$ and retain the same
choice of $g_{n},\sigma_{n}$ as above. Kernel-herding aims to find

$\displaystyle\operatorname*{arg\,min}_{w\in\mathcal{W}}\frac{1}{2}\left\lVert\sum_{n=1}^{N}k(x_{n},\cdot)-\sum_{n=1}^{N}w_{n}k(x_{n},\cdot)\right\rVert_{k}^{2},$
(29)

which, using Theorem 4, is equal to (27). The kernel-herding Frank-Wolfe
updates are

$\displaystyle\overline{u}_{t}=\frac{\sigma}{\sigma_{x_{n}}}\mathbf{1}_{n}\text{
where }n$
$\displaystyle=\operatorname*{arg\,min}_{n\in[N]}\left\langle\sum_{n=1}^{N}k(x_{n},\cdot)-g_{u_{t}},\frac{\sigma}{\sigma_{x_{n}}}k(x_{n},\cdot)\right\rangle_{k},$
(30)
$\displaystyle=\operatorname*{arg\,min}_{n\in[N]}\left\langle\Psi_{\mu}(\pi)-g_{u_{t}},\frac{\sigma}{\sigma_{x_{n}}}\phi_{\mu}(x_{n})\right\rangle_{L^{2}(\mu)},$
(31) $\displaystyle u_{t+1}$
$\displaystyle=(1-\rho_{t})u_{t}+\rho_{t}\overline{u}_{t},$

where (31) is again by Theorem 4.

Due to the equivalences between the optimisation goals (29) and (27) and
updates (30) and (28) it can be concluded that the Frank-Wolfe Hilbert coreset
algorithm is equivalent to Frank-Wolfe kernel-herding when one uses the kernel
(18) whose features map uses the CLR to map data to likelihoods in a Bayes
Hilbert space. This identification immediately proposes some questions. For
example, when kernel-herding is typically employed it does not used the subset
of data requirement enforced in the Frank-Wolfe Bayesian coreset method,
instead it chooses points that are anywhere on the domain which are called
“super-samples” in Chen et al. (2010). This leads to the question of how
kernel-herding performs when used to construct Bayesian coresets without the
limitation of the coreset points being a subset of the observed data points.
Such a methodology was called pseudo-coresets in Manousakas et al. (2020) when
investigating coresets constructed using KL divergence rather than a Hilbert
norm.

###  6.2 Iterative hard thresholding Bayesian coresets

The iterative hard thresholding (IHT) coreset algorithm (Zhang et al., 2021)
proposes to find the weights of a Bayesian coreset using gradient descent and
then thresholds the weights to ensure they non-negative and sparse. More
specifically, the problem under consideration is

$\displaystyle\operatorname*{arg\,min}_{w\in\mathcal{W}_{M}}\frac{1}{2}\left\lVert\Psi_{\mu}(\pi)-\sum_{n=1}^{N}w_{n}\phi_{\mu}(x_{n})\right\rVert_{L^{2}(\mu)}^{2},$
(32)

where $\mathcal{W}_{M}=\\{w\in\mathbb{R}^{N}\colon w\geq 0,\lVert
w\rVert_{0}\leq M\\}$ is the set of vectors in $\mathbb{R}^{N}$ with non-
negative weights and at most $M$ non-zero values. The parameter $M$ dictates
the level of sparseness of the approximation. The IHT technique (Zhang et al.,
2021) solves this problem by performing gradient descent with respect to $w$
at each iteration and then projecting the weights to $\mathcal{W}_{M}$, to
maintain a valid, sparse choice of weights.

Using Theorem 4 is it made apparent that (32) is equivalent to finding a
sparse, non-negative weights solution to the kernel pre-image problem, Problem
1. Gradient descent methods were some of the first methods considered for the
kernel pre-image problem (Burges, 1996) although there were not studied
extensively due to the computational considerations at the time. The IHT
method utilises advances in accelerated optimisation methods to efficiently
find choices of weights that result in a good approximation to the posterior.

The link to the pre-image problem and the IHT method can be taken further by
realising that the IHT method is optimising an MMD criterion. Therefore, IHT
can be seen as performing MMD descent on the empirical measure that
corresponds to the choice of points and weights. MMD descent is a method of
empirical measure approximation that uses gradients to minimise an MMD
objective (Arbel et al., 2019). There are related methods which use other
kernel-based discrepancies (Korba et al., 2021; Xu et al., 2022). This begs
the question of the performance of these other kernel-based descent methods on
the task of constructing Bayesian coresets when using the kernel (18).

##  7 Kullback-Leibler Bayesian Coresets

The aim of this section is to provide a novel connection between Kullback-
Leibler (KL) Bayesian coreset algorithms and Bayes Hilbert spaces. The former
are Bayesian coreset algorithms which use the Kullback-Leibler divergence to
measure the discrepancy between the target posterior and the Bayesian coreset
posterior. The primary focus will be on showing that the quasi-newton KL
Bayesian coreset algorithm (Naik et al., 2022) is equivalent, up to hyper-
parameter choices, to variational inference in a Bayes Hilbert space. This
novel relationship will be shown to be a consequence of the description of
variational inference techniques in terms of Bayes Hilbert spaces provided by
Barfoot and D’Eleuterio (2023).

The KL divergence, see Section 2.1 is a divergence that is used often in
variational inference (Blei et al., 2017). This is a measure approximation
method where the KL divergence is minimised over a user chosen set of
measures, called the variational family. The variational family is often
parameterised using the parameters of a given distribution. For example, the
variational family could be a mixture of Gaussian and the variational
parameters would be the mixture weights and the means and covariances of each
Gaussian. The optimisation to find the best approximation within the
variational family typically involves gradient descent updates. This does not
require samples from the target measure, an essential property in the case
where the target measure is expensive to sample from, for example when it is
the posterior after observing a large number of data points.

Variational inference was framed in terms of Bayes Hilbert spaces by Barfoot
and D’Eleuterio (2023) and this perspective is now outlined. Let $B^{2}(\mu)$
be a Bayes Hilbert space, $\nu\in B^{2}(\mu)$ be the target measure one wishes
to approximate and $\\{b_{n}\\}_{n=1}^{N}\subset B^{2}(\mu)$ some set of $N$
elements which form the dictionary of the approximation. The variational
family is $\\{\oplus_{n=1}^{N}w_{n}\odot b_{n}\colon w\in\mathbb{R}^{N}\\}$
with corresponding optimisation problem

$\displaystyle\operatorname*{arg\,min}_{w\in\mathbb{R}^{N}}\text{KL}\left(\bigoplus_{n=1}^{N}w_{n}\odot
b_{n}\;\|\;\nu\right),$ (33)

where the notation $\oplus,\odot$ is defined in Section 2.1.

The key insight of Barfoot and D’Eleuterio (2023) is that solving this by a
quasi-Newton gradient descent method with respect to $w$ is equivalent to
doing iterative projections in a Bayes Hilbert space. Specifically, Barfoot
and D’Eleuterio (2023, Equation 41) provide the quasi-Newton update

$\displaystyle
w_{t+1}=w_{t}+G_{\nu_{t}}(\textbf{b})^{-1}\langle\textbf{b},\nu\ominus\nu_{t}\rangle_{B^{2}(\nu_{t})},$
(34)

where $\nu_{t}$ is the approximation at iteration $t$ using weights $w_{t}$,
$G_{\nu_{t}}(\textbf{b})$ is the $N\times N$ Gram matrix of the approximating
dictionary $\\{b_{n}\\}_{n=1}^{N}$ with $(n,m)$-th entry $\langle
b_{n},b_{m}\rangle_{B^{2}(\nu_{t})}$ and
$\langle\textbf{b},\nu\ominus\nu_{t}\rangle_{B^{2}(\nu_{t})}$ is the
$\mathbb{R}^{N}$ vector with $n$-th entry $\langle
b_{n},\nu\ominus\nu_{t}\rangle_{B^{2}(\nu_{t})}$. See Section 2.1 for the
definition of $\ominus$ and the Bayes Hilbert space inner product. It was
noted that this quasi-Newton gradient update can be written as

$\displaystyle
w_{t+1}=G_{\nu_{t}}(\textbf{b})^{-1}\langle\textbf{b},\nu\rangle_{B^{2}(\nu_{t})},$

see Barfoot and D’Eleuterio (2023, Equation 43), meaning that at each
iteration the target $\nu$ is being projected to the span of
$\\{b_{n}\\}_{n=1}^{N}$ with respect to the $B^{2}(\nu_{t})$ inner product.
This shows how the geometry of Bayes Hilbert spaces links with variational
inference methods.

The connection between the quasi-Newton Bayesian coreset algorithm derived by
Naik et al. (2022) and this Bayes Hilbert space perspective on variational
inference is now made. The method of Naik et al. (2022) starts by setting the
target measure to be $\pi$, the posterior, and choosing some $M\in\mathbb{N}$.
Then a size $M$ subset from $\\{x_{n}\\}_{n=1}^{N}$ is uniformly sampled,
which has its indices re-ordered and is denoted $\\{x_{m}\\}_{m=1}^{M}$. The
weights for the log-likelihood terms $\\{l_{x_{m}}\\}_{m=1}^{M}$ are then
optimised using a quasi-Newton iterations and thresholding the weights at each
iteration to be non-negative.

The connection to Barfoot and D’Eleuterio (2023) is made by setting $\nu=\pi$,
$N=M$ and $b_{m}=l_{x_{m}}$ in (33). This means the optimisation is

$\displaystyle\operatorname*{arg\,min}_{w\in\mathbb{R}^{M}_{\geq
0}}\text{KL}\left(\bigoplus_{m=1}^{M}w_{m}\odot
b_{m}\;\|\;\pi\right)=\operatorname*{arg\,min}_{w\in\mathbb{R}^{M}_{\geq
0}}\text{KL}(\pi_{w,x}\;\|\;\pi),$

which is equal to (33) except with an additional non-negative weight
constraint.

Reading off the quasi-Newton iteration (Naik et al., 2022, Equation 13) one
sees it is equal, up to hyper-parameter choices, to the update (34) since the
covariance terms in Naik et al. (2022, Equation 13) coincide with the
$B^{2}(\nu_{t})$ inner product. The method of Naik et al. (2022) then
thresholds the weights to ensure non-negativity.

Overall, this shows that the quasi-Newton KL coreset method can be viewed in
terms of iterative projections in a Bayes Hilbert space. This cross
pollination of perspectives is productive since it provides a Bayes Hilbert
space geometric interpretation to the quasi-Newton Bayesian coreset algorithm.
In turn, the numerical and theoretical investigations of Naik et al. (2022)
bolster the geometric perspective. Given that Bayes Hilbert spaces are tangent
spaces in the information geometry sense (Fukumizu, 2009; Pistone, 2013) the
quasi-Newton coreset method (Naik et al., 2022) can be seen as iterative
updates of projecting in a direction in a tangent space and then updating the
base location of the tangent space. This is also related to natural gradient
descent as noted by Barfoot and D’Eleuterio (2023). This perspective was
introduced for the greedy KL coreset method in Campbell and Beronov (2019) and
the results of this section show it applies in a broader sense via the Bayes
Hilbert space perspective.

##  8 Conclusion

This manuscript has described Bayes Hilbert spaces and shown how they are
appropriate spaces to perform posterior measure approximation. Along the way a
novel bound relating the Bayes Hilbert space norm to common discrepancies was
provided as well as a novel connection between kernel-based discrepancies and
Bayes Hilbert spaces in the context of posterior approximation. Multiple
Bayesian coreset methods were expressed in terms of Bayes Hilbert spaces which
provided insight and new interpretations of how they work.

Future avenues of research relating to Bayes Hilbert spaces, Bayesian coresets
and posterior approximation are legion. Three main avenues are now outlined.
First, an issue noted when defining Bayes Hilbert spaces was that they are in
a sense too big, because they can contain infinite measures. The related work
in information geometry (Fukumizu, 2009) defines a related set of measures
that is homeomorphic to a subset of a Hilbert space that has stronger
integrability conditions, but the set of measures is not a vector space
whereas the Bayes Hilbert space is. Therefore, a question is what is an
appropriate space with stricter integrability conditions to maintain finite
measures while also having a vector space structure. Second, as noted in van
den Boogaart et al. (2014) it is straight forward to construct a basis for the
Bayes Hilbert space by using orthogonal polynomials. The specific case of a
Gaussian base measure and Hermite polynomials was studied in Barfoot and
D’Eleuterio (2023). There have been many recent advances in the study of
sparse, high-dimensional approximation using sparse polynomials (Adcock et
al., 2022). This begs the question of how these methods can be applied to the
Bayes Hilbert space. Using polynomials for likelihood approximation in the
context of Bayesian computation was studied in Huggins et al. (2017) and
combining this approach with the aforementioned advances could provide further
advantages. Finally, there have been many recent innovations in distribution
compression and quantisation of empirical measures (Riabiz et al., 2022;
Dwivedi and Mackey, 2021). Given that Theorem 4 links this problem to the
Bayesian coreset problem the question immediately arises as to how to leverage
thee innovations in the context of posterior approximation, in particular in
the case where the points to form the likelihood $\\{z_{m}\\}_{m=1}^{M}$ are
not a subset of the observed data. This pseudo-coreset method has only been
explored for Bayesian coresets in the KL coreset setting (Manousakas et al.,
2020) but it is a natural thing to do in distribution compression setting
(Chen et al., 2010).

##  9 Appendix

###  9.1 Proof of Theorem 3

The proof of this result is simply an adaptation of Theorem 5, Theorem 11 and
Theorem 14 from Sprungk (2020). All the proofs leverage a local Lipschitz
continuity which follows from the assumption involving the bound $B$ on the
CLR transforms.

Recall that the measures in question $\eta,\nu$ are written as

$\displaystyle\frac{\mathrm{d}\eta}{\mathrm{d}\mu}$
$\displaystyle=Z_{\eta}^{-1}\exp(\Psi_{\mu}(\eta))\qquad
Z_{\eta}=\mathbb{E}_{\mu}[\exp(\Psi_{\mu}(\eta))]$ (35)
$\displaystyle\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$
$\displaystyle=Z_{\nu}^{-1}\exp(\Psi_{\mu}(\nu))\qquad
Z_{\nu}=\mathbb{E}_{\mu}[\exp(\Psi_{\mu}(\nu))].$ (36)

By Jensen’s inequality, this means that
$Z_{\eta}\geq\exp(\mathbb{E}_{\mu}[\Psi_{\mu}(\eta)])=1$ since
$\Psi_{\mu}(\eta)\in L^{2}_{0}(\mu)$, the same lower bound applies to
$Z_{\nu}$ too. Therefore, $Z_{\eta}^{-1},Z_{\nu}^{-1}\leq 1$.

#### Hellinger Distance

First the Hellinger distance bound is derived. Starting from the definition
(7) and using $\mu$ as the base measure

$\displaystyle 2\text{H}(\eta,\nu)^{2}$
$\displaystyle=\int_{\Theta}\left(e^{\frac{1}{2}\Psi_{\mu}(\eta)}Z_{\eta}^{-\frac{1}{2}}-e^{\frac{1}{2}\Psi_{\mu}(\nu)}Z_{\nu}^{-\frac{1}{2}}\right)^{2}\mathrm{d}\mu$
$\displaystyle\leq
2\int_{\Theta}\left(e^{\frac{1}{2}\Psi_{\mu}(\eta)}Z_{\eta}^{-\frac{1}{2}}-e^{\frac{1}{2}\Psi_{\mu}(\nu)}Z_{\eta}^{-\frac{1}{2}}\right)^{2}\mathrm{d}\mu$
$\displaystyle+2\int_{\Theta}\left(e^{\frac{1}{2}\Psi_{\mu}(\nu)}Z_{\eta}^{-\frac{1}{2}}-e^{\frac{1}{2}\Psi_{\mu}(\nu)}Z_{\nu}^{-\frac{1}{2}}\right)^{2}\mathrm{d}\mu$
$\displaystyle\eqqcolon I_{1}+I_{2},$

where the inequality $(a-b)^{2}\leq 2(a-c)^{2}+2(c-b)^{2}$ has been used which
is a result of the triangle inequality. The terms $I_{1},I_{2}$ will now be
upper bounded.

$\displaystyle I_{1}$
$\displaystyle=2Z_{\eta}^{-1}\int_{\Theta}\left(e^{\frac{1}{2}\Psi_{\mu}(\eta)}-e^{\frac{1}{2}\Psi_{\mu}(\nu)}\right)^{2}\mathrm{d}\mu$
$\displaystyle\leq
2e^{B}\int_{\Theta}\left(\frac{\Psi_{\mu}(\eta)}{2}-\frac{\Psi_{\mu}(\nu)}{2}\right)^{2}\mathrm{d}\mu$
(37)
$\displaystyle=\frac{1}{2}e^{B}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{2}(\mu)}^{2}=\frac{1}{2}e^{B}\lVert\eta-\nu\rVert_{B^{2}(\mu)}^{2},$

where (37) uses $\lvert e^{x}-e^{y}\rvert\leq e^{\max(x,y)}\lvert
x-y\rvert\>\forall\>x,y,\in\mathbb{R}$ and the assumed upper bound $B$ on
$\Psi_{\mu}(\eta),\Psi_{\mu}(\nu)$. To bound $I_{2}$ first note that by the
definition of the normalisation constants,

$\displaystyle I_{2}$
$\displaystyle=2\int_{\Theta}\left(e^{\frac{1}{2}\Psi_{\mu}(\nu)}Z_{\eta}^{-\frac{1}{2}}-e^{\frac{1}{2}\Psi_{\mu}(\nu)}Z_{\nu}^{-\frac{1}{2}}\right)^{2}\mathrm{d}\mu$
$\displaystyle=2Z_{\eta}^{-1}\left(Z_{\eta}^{\frac{1}{2}}-Z_{\nu}^{\frac{1}{2}}\right)^{2}\leq
2\left(Z_{\eta}^{\frac{1}{2}}-Z_{\nu}^{\frac{1}{2}}\right)^{2}.$

Now, $\lvert x^{1/2}-y^{1/2}\rvert\leq\frac{1}{2}\min(x,y)^{-1/2}\lvert
x-y\rvert\>\forall\>x,y>0$ gives

$\displaystyle I_{2}\leq\frac{1}{\min(Z_{\eta},Z_{\nu})}\lvert
Z_{\eta}-Z_{\nu}\rvert^{2}\leq\frac{1}{2}\lvert Z_{\eta}-Z_{\nu}\rvert^{2}.$

Finally,

$\displaystyle\lvert Z_{\eta}-Z_{\nu}\rvert^{2}\leq\int_{\Theta}\left\lvert
e^{\Psi_{\mu}(\eta)}-e^{\Psi_{\mu}(\nu)}\right\rvert^{2}\mathrm{d}\mu\leq
e^{2B}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{2}(\mu)}^{2}=e^{2B}\lVert\eta-\nu\rVert_{B^{2}(\mu)}^{2},$

where the same bounding technique as (37) was used. Putting this together
gives

$\displaystyle\text{H}(\eta,\nu)^{2}\leq\frac{1}{4}\left(e^{B}+e^{2B}\right)\lVert\eta-\nu\rVert_{B^{2}(\mu)}^{2},$

as desired.

#### KL Divergence

The case for KL divergence is now covered. First note

$\displaystyle\frac{\mathrm{d}\eta}{\mathrm{d}\nu}=\frac{\mathrm{d}\eta}{\mathrm{d}\mu}\frac{\mathrm{d}\mu}{\mathrm{d}\nu}=Z_{\eta}^{-1}Z_{\nu}e^{\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)},$

therefore

$\displaystyle\text{KL}(\eta\;\|\;\nu)=\int_{\Theta}\log\frac{\mathrm{d}\eta}{\mathrm{d}\nu}\mathrm{d}\eta$
$\displaystyle\leq\left\lvert\log Z_{\eta}-\log
Z_{\nu}\right\rvert+\int_{\Theta}\left\lvert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\right\rvert
Z_{\eta}^{-1}e^{\Psi_{\mu}(\eta)}\mathrm{d}\mu$ $\displaystyle\eqqcolon
I_{1}+I_{2}.$

Using $\lvert\log x-\log y\rvert\leq\min(x,y)^{-1}\lvert x-y\rvert\>\forall
x,y>0$ gives

$\displaystyle I_{1}$ $\displaystyle\leq\min(Z_{\eta},Z_{\nu})^{-1}\lvert
Z_{\eta}-Z_{\nu}\rvert$ $\displaystyle\leq
e^{B}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{1}(\mu)}$
$\displaystyle\leq
e^{B}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{2}(\mu)}=e^{B}\lVert\eta-\nu\rVert_{B^{2}(\mu)},$

where the local Lipschitz bound $\lvert e^{x}-e^{y}\rvert\leq
e^{\max(x,y)}\lvert x-y\rvert\>\forall\>x,y,\in\mathbb{R}$ was used to get the
$L^{1}(\mu)$ norm upper bound, similar to its use in (37). The $L^{2}(\mu)$
norm upper bounds the $L^{1}(\mu)$ norm since $\mu$ is a probability measure.
If $\mu$ was assumed only to be a finite measure then there would be an extra
factor of $\mu(\Theta)$ in the bounds. Bounding $I_{2}$ is done in a
straightforward way

$\displaystyle
I_{2}=\int_{\Theta}\left\lvert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\right\rvert
Z_{\eta}^{-1}e^{\Psi_{\mu}(\eta)}\mathrm{d}\mu$ $\displaystyle\leq
e^{B}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{1}(\mu)}$
$\displaystyle\leq
e^{B}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{2}(\mu)}=e^{B}\lVert\eta-\nu\rVert_{B^{2}(\mu)}.$

Putting $I_{1},I_{2}$ together gives the desired bound for the KL divergence.

#### Wasserstein-$1$ Distance

Finally, the Wasserstein-$1$ case is dealt with. Take any
$\theta_{0}\in\Theta$ then by simply shifting the values it suffices to
consider the supremum over functions in $\text{Lip}(1)$ such that
$f(\theta_{0})=0$. For such functions $\lvert f(\theta)\rvert\leq
d_{\Theta}(\theta,\theta_{0})$. Following Sprungk (2020, Theorem 14),

$\displaystyle\left\lvert\int_{\Theta}f\mathrm{d}\eta-\int_{\Theta}f\mathrm{d}\nu\right\rvert$
$\displaystyle=\left\lvert\int_{\Theta}f\cdot\left(Z_{\eta}^{-1}e^{\Psi_{\mu}(\eta)}-Z_{\nu}^{-1}e^{\Psi_{\mu}(\nu)}\right)\mathrm{d}\mu\right\rvert$
$\displaystyle\leq\left\lvert
Z_{\eta}^{-1}-Z_{\nu}^{-1}\right\rvert\left\lvert\int_{\Theta}fe^{\Psi_{\mu}(\eta)}\mathrm{d}\mu\right\rvert+\left|Z_{\nu}^{-1}\int_{\Theta}f\cdot\left(e^{\Psi_{\mu}(\eta)}-e^{\Psi_{\mu}(\nu)}\right)\mathrm{d}\mu\right\rvert$
$\displaystyle\eqqcolon I_{1}+I_{2},$

where the inequality is by the triangle inequality. To bound $I_{1}$ note

$\displaystyle\left\lvert
Z_{\eta}^{-1}-Z_{\nu}^{-1}\right\rvert=\frac{\left\lvert
Z_{\eta}-Z_{\nu}\right\rvert}{Z_{\eta}Z_{\nu}}\leq\left\lvert
Z_{\eta}-Z_{\nu}\right\rvert\leq
e^{B}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{1}(\mu)},$

where the first inequality is by the lower bound of $1$ on the normalising
constants and the second by using again the local Lipschitz result, as was
done in (37). Using this gives

$\displaystyle I_{1}\leq\int_{\Theta}\left\lvert
fe^{\Psi_{\mu}(\eta)}\right\rvert\mathrm{d}\mu\cdot
e^{B}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{1}(\mu)}\leq
e^{2B}\lVert\mu\rVert_{\mathcal{P}^{1}}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{1}(\mu)},$

where the second inequality is by the bounds on $\lvert f(\theta)\rvert$,
$\Psi_{\mu}(\eta)$ and the way the choice of $\theta_{0}$ was arbitrary, which
allows the infimum bound $norm{\mu}_{\mathcal{P}^{1}}$.

Bounding $I_{2}$ is largely similar,

$\displaystyle I_{2}\leq\int_{\Theta}\lvert
d_{\Theta}(\theta,\theta_{0})\rvert\left\lvert
e^{\Psi_{\mu}(\eta)}-e^{\Psi_{\mu}(\nu)}\right\rvert\mathrm{d}\mu\leq
e^{B}\lVert\mu\rVert_{\mathcal{P}^{2}}\int_{\Theta}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{2}(\mu)},$

where the second inequality is by the Cauchy-Schwarz inequality and the local
Lipschitz bound again. Overall, using the Cauchy-Schwartz inequality along
with how $\mu$ is a probability measure to bound the $L^{1}$ norm with the
$L^{2}$ norm gives

$\displaystyle\text{W}_{1}(\eta,\nu)\leq\left(e^{B}+e^{2B}\right)\lVert\mu\rVert_{\mathcal{P}^{2}}\lVert\Psi_{\mu}(\eta)-\Psi_{\mu}(\nu)\rVert_{L^{2}(\mu)}=\left(e^{B}+e^{2B}\right)\lVert\mu\rVert_{\mathcal{P}^{2}}\lVert\eta-\nu\rVert_{B^{2}(\mu)}.$

###  9.2 Proof of Corollary 1

The proof is a simple application of Schneider (2016, Theorem 2). First, note
that by Theorem 4

$\displaystyle\left\lVert\pi_{M}-\pi\right\rVert_{B^{2}(\mu)}=\text{MMD}_{k}(P_{M},P_{N})=\left\lVert\Phi_{k}P_{M}-\Phi_{k}P_{N}\right\rVert_{k},$

where $k$ is the kernel (18), $P_{M}=\frac{N}{M}\sum_{m=1}^{M}z_{m}$ and
$P_{N}=\sum_{n=1}^{N}\delta_{x_{n}}$. Therefore, the quantity to be bounded is
the same as in Schneider (2016, Theorem 2) up to an additional scaling of $N$
since the referenced result involves the normalised empirical measures
$\widetilde{P}_{M}=\frac{1}{N}P_{M},\widetilde{P}_{N}=\frac{1}{N}P_{N}$. This
means the bound in Schneider (2016, Theorem 2) implies for any
$\varepsilon_{0}$

$\displaystyle\mathbb{P}\left(\left\lVert\Phi_{k}\widetilde{P}_{M}-\Phi_{k}\widetilde{P}_{N}\right\rVert_{k}\geq\varepsilon_{0}\right)=\mathbb{P}(\left\lVert\Phi_{k}P_{M}-\Phi_{k}P_{N}\right\rVert_{k}\geq
N\varepsilon_{0})\leq
2\exp\left(-\frac{M\varepsilon_{0}^{2}}{8\gamma^{2}(1-(M-1)/N)}\right),$

and then the final step is simply substituting $\varepsilon=N\varepsilon_{0}$
and rearranging to obtain the bound purely in terms of $\delta$.

## References

  * Adcock et al. [2022] B. Adcock, S. Brugiapaglia, and C. G. Webster.  _Sparse Polynomial Approximation of High-Dimensional Functions_.  Society for Industrial and Applied Mathematics, 2022.
  * Aitchison [1982] J. Aitchison.  The statistical analysis of compositional data.  _Journal of the Royal Statistical Society. Series B. Methodological_ , 44(2):139–177, 1982.  With discussion.
  * Aitchison [1986] J. Aitchison.  _The Statistical Analysis of Compositional Data_.  Chapman and Hall, 1986.
  * Anastasiou et al. [2023] A. Anastasiou, A. Barp, F.-X. Briol, B. Ebner, R. E. Gaunt, F. Ghaderinezhad, J. Gorham, A. Gretton, C. Ley, Q. Liu, L. Mackey, C. J. Oates, G. Reinert, and Y. Swan.  Stein’s method meets computational statistics: A review of some recent developments.  _Statistical Science_ , 38(1):120 – 139, 2023\.
  * Arbel et al. [2019] M. Arbel, A. Korba, A. SALIM, and A. Gretton.  Maximum mean discrepancy gradient flow.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Bach et al. [2012] F. Bach, S. Lacoste-Julien, and G. Obozinski.  On the equivalence between herding and conditional gradient algorithms.  In _Proceedings of the 29th International Coference on International Conference on Machine Learning_ , page 1355–1362, 2012.
  * Barfoot and D’Eleuterio [2023] T. D. Barfoot and G. M. T. D’Eleuterio.  Variational inference as iterative projection in a Bayesian Hilbert space with application to robotic state estimation.  _Robotica_ , 41(2):632–667, 2023.
  * Berlinet and Thomas-Agnan [2004] A. Berlinet and C. Thomas-Agnan.  _Reproducing Kernel Hilbert Spaces in Probability and Statistics_.  Springer US, 2004.
  * Bierkens et al. [2019] J. Bierkens, P. Fearnhead, and G. Roberts.  The Zig-Zag process and super-efficient sampling for Bayesian analysis of big data.  _The Annals of Statistics_ , 47(3):1288 – 1320, 2019.
  * Blei et al. [2017] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American Statistical Association_ , 112(518):859–877, 2017.
  * Bogachev [2007] V. I. Bogachev.  _Measure Theory_.  Springer Berlin Heidelberg, 2007.
  * Burges [1996] C. J. C. Burges.  Simplified support vector decision rules.  In _International Conference on Machine Learning_ , 1996.
  * Caflisch [1998] R. E. Caflisch.  Monte Carlo and quasi-Monte Carlo methods.  _Acta Numerica_ , 7:1–49, 1998.
  * Campbell and Beronov [2019] T. Campbell and B. Beronov.  Sparse variational inference: Bayesian coresets from scratch.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Campbell and Broderick [2018] T. Campbell and T. Broderick.  Bayesian coreset construction via greedy iterative geodesic ascent.  In _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pages 698–706, 2018.
  * Campbell and Broderick [2019] T. Campbell and T. Broderick.  Automated scalable bayesian inference via hilbert coresets.  _Journal of Machine Learning Research_ , 20(15):1–38, 2019.
  * Cena and Pistone [2006] A. Cena and G. Pistone.  Exponential statistical manifold.  _Annals of the Institute of Statistical Mathematics_ , 59(1):27–56, 2006.
  * Chen et al. [2010] Y. Chen, M. Welling, and A. Smola.  Super-samples from kernel herding.  In _Proceedings of the Twenty-Sixth Conference on Uncertainty in Artificial Intelligence_ , page 109–116, 2010.
  * Cherief-Abdellatif and Alquier [2020] B.-E. Cherief-Abdellatif and P. Alquier.  MMD-Bayes: Robust bayesian estimation via maximum mean discrepancy.  In _Proceedings of The 2nd Symposium on Advances in Approximate Bayesian Inference_ , volume 118, pages 1–21, 2020.
  * Christmann and Steinwart [2008] A. Christmann and I. Steinwart.  _Support Vector Machines_.  Springer New York, 2008.
  * Clarkson [2010] K. L. Clarkson.  Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm.  _ACM Transactions on Algorithms_ , 6(4), 2010.
  * Dang et al. [2019] K.-D. Dang, M. Quiroz, R. Kohn, M.-N. Tran, and M. Villani.  Hamiltonian Monte Carlo with energy conserving subsampling.  _Journal of Machine Learning Research_ , 20(100):1–31, 2019.
  * Dick et al. [2013] J. Dick, F. Y. Kuo, and I. H. Sloan.  High-dimensional integration: The quasi-Monte Carlo way.  _Acta Numerica_ , 22:133–288, 2013.
  * Dwivedi and Mackey [2021] R. Dwivedi and L. Mackey.  Kernel thinning.  In _Proceedings of Thirty Fourth Conference on Learning Theory_ , volume 134, pages 1753–1753, 2021.
  * Egozcue et al. [2006] J. J. Egozcue, J. L. Díaz-Barrero, and V. Pawlowsky-Glahn.  Hilbert space of probability density functions based on Aitchison geometry.  _Acta Mathematica Sinica (English Series)_ , 22(4):1175–1182, 2006.
  * Erb and Ay [2021] I. Erb and N. Ay.  The information-geometric perspective of compositional data analysis.  In _Advances in Compositional Data Analysis_ , pages 21–43. Springer International Publishing, 2021.
  * Filzmoser et al. [2018] P. Filzmoser, K. Hron, and M. Templ.  _Applied Compositional Data Analysis_.  Springer International Publishing, 2018.
  * Fukumizu [2009] K. Fukumizu.  _Exponential manifold by reproducing kernel Hilbert spaces_ , page 291–306.  Cambridge University Press, 2009.
  * Fukumizu et al. [2013] K. Fukumizu, L. Song, and A. Gretton.  Kernel Bayes’ rule: Bayesian inference with positive definite kernels.  _Journal of Machine Learning Research_ , 14(118):3753–3783, 2013.
  * Gelman et al. [2013] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin.  _Bayesian Data Analysis_.  Chapman and Hall/CRC, 2013.
  * Graf and Luschgy [2000] S. Graf and H. Luschgy.  _Foundations of Quantization for Probability Distributions_.  Springer Berlin Heidelberg, 2000.
  * Greenacre et al. [2022] M. Greenacre, E. Grunsky, J. Bacon-Shone, I. Erb, and T. Quinn.  Aitchison’s compositional data analysis 40 years on: A reappraisal.  _arXiv:2201.05197_ , 2022.
  * Gretton et al. [2006] A. Gretton, K. Borgwardt, M. Rasch, B. Schölkopf, and A. Smola.  A kernel method for the two-sample-problem.  In _Advances in Neural Information Processing Systems_ , volume 19. MIT Press, 2006.
  * Gretton et al. [2012] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola.  A kernel two-sample test.  _Journal of Machine Learning Research_ , 13(1):723–773, 2012.
  * Guilbart [1978] C. Guilbart.  _Etude des Produits Scalaires sur l’Espace des Mesures: Estimation par Projections_.  PhD thesis, Université des Sciences et Techniques de Lille, 1978.
  * Hoffman et al. [2013] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley.  Stochastic variational inference.  _Journal of Machine Learning Research_ , 14(40):1303–1347, 2013.
  * Hron et al. [2022] K. Hron, J. Machalová, and A. Menafoglio.  Bivariate densities in Bayes spaces: Orthogonal decomposition and spline representation.  _Statistical Papers_ , 2022.
  * Huggins et al. [2016] J. Huggins, T. Campbell, and T. Broderick.  Coresets for scalable Bayesian logistic regression.  In _Advances in Neural Information Processing Systems_ , volume 29, 2016.
  * Huggins et al. [2017] J. Huggins, R. P. Adams, and T. Broderick.  PASS-GLM: Polynomial approximate sufficient statistics for scalable Bayesian GLM inference.  In _Advances in Neural Information Processing Systems_ , volume 30, 2017.
  * Huggins et al. [2018] J. H. Huggins, T. Campbell, M. Kasprzak, and T. Broderick.  Practical bounds on the error of Bayesian posterior approximations: A nonasymptotic approach.  _arXiv:1809.09505_ , 2018.
  * Korba et al. [2021] A. Korba, P.-C. Aubin-Frankowski, S. Majewski, and P. Ablin.  Kernel stein discrepancy descent.  In _Proceedings of the 38th International Conference on Machine Learning_ , volume 139, pages 5719–5730, 2021.
  * Kroese et al. [2011] D. P. Kroese, T. Taimre, and Z. I. Botev.  _Handbook of Monte Carlo Methods_.  Wiley, 2011.
  * Kwok and Tsang [2004] J.-Y. Kwok and I.-H. Tsang.  The pre-image problem in kernel methods.  _IEEE Transactions on Neural Networks_ , 15(6):1517–1525, 2004.
  * Li et al. [2017] C.-L. Li, W.-C. Chang, Y. Cheng, Y. Yang, and B. Poczos.  MMD GAN: Towards deeper understanding of moment matching network.  In _Advances in Neural Information Processing Systems_ , volume 30, 2017.
  * Maier et al. [2021] E.-M. Maier, A. Stöcker, B. Fitzenberger, and S. Greven.  Additive density-on-scalar regression in Bayes Hilbert spaces with an application to gender economics.  _arXiv:2110.11771_ , 2021.
  * Mak and Joseph [2018] S. Mak and V. R. Joseph.  Support points.  _The Annals of Statistics_ , 46(6A):2562–2592, 2018.
  * Manousakas et al. [2020] D. Manousakas, Z. Xu, C. Mascolo, and T. Campbell.  Bayesian pseudocoresets.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 14950–14960, 2020.
  * Mateu and Giraldo [2021] J. Mateu and R. Giraldo, editors.  _Geostatistical Functional Data Analysis_.  Wiley, 2021.
  * Mika et al. [1998] S. Mika, B. Schölkopf, A. Smola, K.-R. Müller, M. Scholz, and G. Rätsch.  Kernel PCA and de-noising in feature spaces.  In _Advances in Neural Information Processing Systems_ , volume 11, 1998.
  * Muandet et al. [2017] K. Muandet, K. Fukumizu, B. Sriperumbudur, and B. Schölkopf.  Kernel mean embedding of distributions: A review and beyond.  _Foundations and Trends® in Machine Learning_ , 10(1-2):1–141, 2017.
  * Naik et al. [2022] C. Naik, J. Rousseau, and T. Campbell.  Fast Bayesian coresets via subsampling and quasi-newton refinement.  In _Advances in Neural Information Processing Systems_ , 2022.
  * Newton [2012] N. J. Newton.  An infinite-dimensional statistical manifold modelled on Hilbert space.  _Journal of Functional Analysis_ , 263(6):1661–1681, 2012.
  * Paulsen and Raghupathi [2016] V. I. Paulsen and M. Raghupathi.  _An Introduction to the Theory of Reproducing Kernel Hilbert Spaces_ , volume 152 of _Cambridge Studies in Advanced Mathematics_.  Cambridge University Press, Cambridge, 2016.
  * Pawlowsky-Glahn and Buccianti [2011] V. Pawlowsky-Glahn and A. Buccianti, editors.  _Compositional Data Analysis_.  Wiley, 2011.
  * Petersen et al. [2022] A. Petersen, C. Zhang, and P. Kokoszka.  Modeling probability density functions as data objects.  _Econometrics and Statistics_ , 21:159–178, 2022.
  * Pistone [2013] G. Pistone.  Nonparametric information geometry.  In _Lecture Notes in Computer Science_ , pages 5–36. Springer Berlin Heidelberg, 2013.
  * Pistone and Sempi [1995] G. Pistone and C. Sempi.  An infinite-dimensional geometric structure on the space of all the probability measures equivalent to a given one.  _The Annals of Statistics_ , 23(5):1543–1561, 1995.
  * Riabiz et al. [2022] M. Riabiz, W. Y. Chen, J. Cockayne, P. Swietach, S. A. Niederer, L. Mackey, and C. J. Oates.  Optimal thinning of MCMC output.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 84(4):1059–1081, 2022.
  * Schneider [2016] M. Schneider.  Probability inequalities for kernel embeddings in sampling without replacement.  In _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics_ , volume 51, pages 66–74, 2016.
  * Schölkopf and Smola [2002] B. Schölkopf and A. J. Smola.  _Learning with Kernels_.  The MIT Press, 2002.
  * Schölkopf et al. [1997] B. Schölkopf, A. Smola, and K.-R. Müller.  Kernel principal component analysis.  In _Lecture Notes in Computer Science_ , pages 583–588. Springer Berlin Heidelberg, 1997.
  * Scholkopf et al. [1999] B. Scholkopf, S. Mika, C. Burges, P. Knirsch, K.-R. Muller, G. Ratsch, and A. Smola.  Input space versus feature space in kernel-based methods.  _IEEE Transactions on Neural Networks_ , 10(5):1000–1017, 1999.
  * Scott et al. [2016] S. L. Scott, A. W. Blocker, F. V. Bonassi, H. A. Chipman, E. I. George, and R. E. McCulloch.  Bayes and big data: The consensus Monte Carlo algorithm.  _International Journal of Management Science and Engineering Management_ , 11(2):78–88, 2016.
  * Sprungk [2020] B. Sprungk.  On the local lipschitz stability of Bayesian inverse problems.  _Inverse Problems_ , 36(5), 2020.
  * Sriperumbudur et al. [2017] B. Sriperumbudur, K. Fukumizu, A. Gretton, A. Hyvärinen, and R. Kumar.  Density estimation in infinite dimensional exponential families.  _Journal of Machine Learning Research_ , 18(57):1–59, 2017.
  * Sriperumbudur et al. [2010] B. K. Sriperumbudur, A. Gretton, K. Fukumizu, B. Schölkopf, and G. R. Lanckriet.  Hilbert space embeddings and metrics on probability measures.  _Journal of Machine Learning Research_ , 11:1517–1561, 2010\.
  * Srivastava et al. [2015] S. Srivastava, V. Cevher, Q. Dinh, and D. Dunson.  WASP: Scalable Bayes via barycenters of subset posteriors.  In _Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics_ , volume 38, pages 912–920, 2015.
  * Stuart [2010] A. M. Stuart.  Inverse problems: A Bayesian perspective.  _Acta Numerica_ , 19:451–559, 2010.
  * Teymur et al. [2021] O. Teymur, J. Gorham, M. Riabiz, and C. Oates.  Optimal quantisation of probability measures using maximum mean discrepancy.  In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130, pages 1027–1035, 2021.
  * van den Boogaart and Tolosana-Delgado [2013] K. G. van den Boogaart and R. Tolosana-Delgado.  _Analyzing Compositional Data with R_.  Springer Berlin Heidelberg, 2013.
  * van den Boogaart et al. [2014] K. G. van den Boogaart, J. J. Egozcue, and V. Pawlowsky-Glahn.  Bayes Hilbert spaces.  _Australian & New Zealand Journal of Statistics_, 56(2):171–194, 2014.
  * van den Boogart et al. [2011] K.-G. van den Boogart, E. Juan José, and P.-G. Vera.  Bayes linear spaces.  _SORT-Statistics and Operational Research Transasctions_ , 34(2):201–222, 2011.
  * Vyner et al. [2023] C. Vyner, C. Nemeth, and C. Sherlock.  SwISS: A scalable markov chain monte carlo divide-and-conquer strategy.  _Stat_ , 12(1), 2023.
  * Wynne and Nagy [2021] G. Wynne and S. Nagy.  Statistical depth meets machine learning: Kernel mean embeddings and depth in functional data analysis.  _arXiv:2105.12778_ , 2021.
  * Xu et al. [2022] L. Xu, A. Korba, and D. Slepcev.  Accurate quantization of measures via interacting particle-based optimization.  In _Proceedings of the 39th International Conference on Machine Learning_ , volume 162, pages 24576–24595, 2022.
  * Zhang et al. [2021] J. Zhang, R. Khanna, A. Kyrillidis, and S. Koyejo.  Bayesian coresets: Revisiting the nonconvex optimization perspective.  In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130, pages 2782–2790, 2021.