Stein Π-Importance Sampling
Congye Wang^1, Wilson Chen^2, Heishiro Kanagawa^1, Chris. J. Oates^1

^1 Newcastle University, UK

^2 University of Sydney, Australia

Stein discrepancies have emerged as a powerful tool for retrospective improvement of Markov chain Monte Carlo output.
However, the question of how to design Markov chains that are well-suited to such post-processing has yet to be addressed.
This paper studies Stein importance sampling, in which weights are assigned to the states visited by a $\Pi$-invariant Markov chain to obtain a consistent approximation of $P$, the intended target.
Surprisingly, the optimal choice of $\Pi$ is not identical to the target $P$; we therefore propose an explicit construction for $\Pi$ based on a novel variational argument.
Explicit conditions for convergence of Stein $\Pi$-Importance Sampling are established.
For $\approx 70\%$ of tasks in the  benchmark, a significant improvement over the analogous post-processing of $P$-invariant Markov chains is reported.

§ INTRODUCTION

Stein discrepancies are a class of statistical divergences that can be computed without access to a normalisation constant.
Originally conceived as a tool to measure the performance of sampling methods [27], these discrepancies have since found wide-ranging statistical applications <cit.>.
Our focus here is the use of Stein discrepancies for retrospective improvement of mcmc, and here two main techniques have been proposed:
(i) Stein importance sampling [40, 33], and (ii) Stein thinning [46].
In Stein importance sampling (also called black box importance sampling), the samples
are assigned weights such that a Stein discrepancy between the weighted empirical measure and the target $P$ is minimised.
Stein thinning constructs a sparse approximation to this optimally weighted measure at a lower computational and storage cost.
Together, these techniques provide a powerful set of post-processing tools for mcmc, with subsequent authors proposing a range of generalisations and extensions [52, 16, 31, 24, 9].

The consistency of these algorithms has been established in the setting of approximate, $\Pi$-invariant mcmc, motivated by challenging inference problems where only approximate sampling can be performed.
In these settings, $\Pi$ is implicitly an approximation to $P$ that is as accurate as possible subject to computational budget.
However, the critical question of how to design Markov chains that are well-suited to such post-processing has yet to be addressed.
This paper provides a solution, in the form of a specific construction for $\Pi$ derived from a novel variational argument.
Surprisingly, we are able to demonstrate a substantial improvement using the proposed $\Pi$, compared to the case where $\Pi$ and $P$ are equal.
The paper proceeds as follows:
<Ref> presents an abstract formulation of the task and existing results for optimally-weighted empirical measures are reviewed.
<Ref> derives our proposed choice of $\Pi$ and establishes that Stein post-processing of samples from a $\Pi$-invariant mala provides a consistent approximation of $P$.
The approach is stress-tested using the recently released  suite of benchmark tasks in <Ref>, before concluding with a discussion in <Ref>.

§ BACKGROUND

To properly contextualise our discussion we start with an abstract mathematical description of the task.
Let $P$ be a probability measure on a measurable space $\mathcal{X}$.
Let $\mathcal{P}(\mathcal{X})$ be the set of all probability measures on $\mathcal{X}$.
Let $D_P : \mathcal{P}(\mathcal{X}) \rightarrow [0,\infty]$ be a statistical divergence for measuring the quality of an approximation $Q$ to $P$, meaning that $D_P(Q) = 0$ if and only if $Q=P$.
In this work we consider approximations whose support is contained in a finite set $\{ x_1,\dots,x_n \} \subset \mathcal{X}$, and in particular we consider optimal approximations of the form

\begin{align*}
P_n^\star = \sum_{i=1}^n w_i^\star \delta(x_i), \qquad w^\star \in \argmin_{w \geq 0, \; 1^\top w = 1} D_P \left( \sum_{i=1}^n w_i \delta(x_i)  \right) .
\end{align*}

In what follows we restrict attention to statistical divergences for which such approximations can be shown to exist and be well-defined.
The question that we then ask is which states $\{x_1,\dots,x_n\}$ minimise the approximation error $D_P(P_n^\star)$?
Before specialising to Stein discrepancies, it is helpful to review existing results for some standard statistical divergences $D_P$.

§.§ Wasserstein Divergence

Optimal quantisation focuses on the $r$-Wasserstein ($r \geq 1$) family of statistical divergences $D_P(Q) = \inf_{\gamma \in \Gamma(P,Q)} \int \|x - y\|^r \mathrm{d}\gamma(x,y)$, where $\Gamma(P,Q)$ denotes the set of all couplings[A coupling $\gamma \in \Gamma(P,Q)$ is a distribution $\gamma \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d)$ whose marginal distributions are $P$ and $Q$. ] of $P,Q \in \mathcal{P}(\mathbb{R}^d)$, and the divergence is finite whenever $P$ and $Q$ have finite $r$-th moment.
Assuming the states $\{x_1,\dots,x_n\}$ are distinct, the corresponding optimal weights are $w_i^\star = P(A_i)$ where $A_i$ is the Voronoi neighbourhood[The Voronoi neighbourhood of $x_i$ is the set $A_i = \{ x \in \mathbb{R}^d : \|x-x_i\| = \min_{j=1,\dots,n} \|x-x_j\| \}$.] of $x_i$ in $\mathbb{R}^d$.
Optimal states achieve the minimal quantisation error for $P$;

\begin{align*}
e_{n,r}(P) & = \inf_{x_1,\dots,x_n \in \mathbb{R}^d } D_P \left( \sum_{i=1}^n w_i^\star \delta(x_i) \right) ,
\end{align*}

the smallest value of the divergence among optimally-weighted distributions supported on at most $n$ states.
Though the dependence of optimal states on $n$ and $P$ can be complicated, we can broaden our perspective to consider asymptotically optimal states, whose asymptotic properties can be precisely characterised.
To this end, for $A \subset \mathbb{R}^d$, let $\mathcal{U}(A)$ denote the uniform distribution on $A$, and define the universal constant $\smash{C_r([0,1]^d) = \inf_{n \geq 1} n^{r/d} \;  e_{n,r}(\mathcal{U}([0,1]^d))}$.
Suppose that $P$ admits a density $p$ on $\mathbb{R}^d$.
Then the $r$th quantisation coefficient of $P$ on $\mathbb{R}^d$, defined as

C_r(P) = C_r([0,1]^d)  \left( \int p(x)^{d/(d+r)} \; \mathrm{d}x \right)^{(d+r)/d} ,

plays a central role in the classical theory of quantisation, being the rate constant in the asymptotic convergence of the minimal quantisation error; $\smash{\lim_{n \rightarrow \infty} n^{r/d} e_{n,r}(P) = C_r(P)}$; see Theorem 6.2 of 30.
This suggests a natural definition; a collection $\{x_1,\dots,x_n\}$ is called asymptotically optimal if

\lim_{n \rightarrow \infty} n^{r/d}  D_P \left( \sum_{i=1}^n w_i^\star \delta(x_i) \right) = C_r(P) ,

which amounts to $P_n^\star$ asymptotically attaining the minimal quantisation error $e_{n,r}(P)$.
The main result here is that, if $\{x_1,\dots,x_n\}$ are asymptotically optimal, then $\smash{\frac{1}{n} \sum_{i=1}^n \delta(x_i) \rightarrow \Pi_r}$, where convergence is in distribution and $\Pi_r$ is the distribution whose density is $\smash{\pi_r(x) \propto p(x)^{d/(d+r)}}$; see Theorem 7.5 of 30.
This provides us with a key insight; optimal states are over-dispersed with respect to the intended distributional target.
The extent of the over-dispersion here depends both on $r$, a parameter of the statistical divergence, and the dimension $d$ of the space on which distributions are defined.

The $r$-Wasserstein divergence is, unfortunately, not well-suited for use in the motivating Bayesian context.
In particular, computing the optimal weights $w_i = P(A_i)$ requires knowledge of $P$, which is typically not available when $P$ is implicitly defined via an intractable normalisation constant.
On the other hand, the optimal sampling distribution $\Pi$ is explicit and can be sampled (for example using mcmc); for discussion of random quantisers in this context see <cit.>, <cit.> and <cit.>.
The simple form of $\Pi$ is a feature of the classical approach to quantisation that we will attempt to mimic in the sequel.

§.§ Kernel Discrepancies

The theory of quantisation using kernels is less well-developed.
A kernel is a measurable, symmetric, positive-definite function $k : \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}$.
From the Moore–Aronszajn theorem, there is a unique Hilbert space $\mathcal{H}(k)$ for which $k$ is a reproducing kernel, meaning that $k(\cdot,x) \in \mathcal{H}(k)$ for all $x \in \mathcal{X}$ and $\langle f , k(\cdot,x) \rangle_{\mathcal{H}(k)} = f(x)$ for all $f \in \mathcal{H}(k)$ and all $x \in \mathcal{X}$.
Assuming that $\mathcal{H}(k) \subset L^1(P)$, we can define the weak (or Pettis) integral

\begin{align}
\mu_P(\cdot) = \int k(\cdot,x) \; \mathrm{d}P(x) ,
\label{eq: pettis}
\end{align}

called the kernel mean embedding of $P$ in $\mathcal{H}(k)$.
The kernel discrepancy is then defined as the norm of the difference between kernel mean embeddings

\begin{align}
D_P(Q) = \|\mu_Q - \mu_P\|_{\mathcal{H}(k)} = \sqrt{\iint k(x,y) \; \mathrm{d}(Q-P)(x) \mathrm{d}(Q-P)(y)} \label{eq: MMD}
\end{align}

where, to be consistent with our earlier notation, we adopt the convention that $D_P(Q)$ is infinite whenever $\mathcal{H}(k) \not\subset L^1(Q)$.
The second equality in (<ref>) follows immediately from the stated properties of a reproducing kernel.
To satisfy the requirement of a statistical divergence, we assume that the kernel $k$ is characteristic, meaning that $\mu_P = \mu_Q$ if and only if $P=Q$.
In this setting, the properties of optimal states are necessarily dependent on the choice of kernel $k$, and are in general not well-understood.
Indeed, given distinct states $\{x_1,\dots,x_n\}$, the corresponding optimal weights $w^\star = (w_1^\star,\dots,w_n^\star)^\top$ are the solution to the linearly-constrained quadratic program

\begin{align}
\argmin_{w \in \mathbb{R}^d} \; w^\top K w - 2 z^\top w \qquad \text{s.t.} \qquad w \geq 0 , \; 1^\top w = 1  \label{eq: constrained program}
\end{align}

where $K_{i,j} = k(x_i,x_j)$ and $z_i = \mu_P(x_i)$.
This program does not admit a closed-form solution, but can be numerically solved.
To the best of our knowledge, the only theoretical analysis of approximations based on (<ref>) is due to
32, who established rates for the convergence of $P_n^\star$ to $P$ in the case where states are independently sampled from $P$.
The question of an optimal sampling distribution was not considered in that work.

Although few results are available concerning (<ref>), relaxations of this program have been well-studied.
The simplest relaxation of (<ref>) is to remove both the positivity ($w \geq 0$) and normalisation ($1^\top w = 1$) constraints, in which case the optimal weights have the explicit representation $w^* = K^{-1} z$.
The analysis of optimal states in this context has developed under the dual strands of kernel cubature and Bayesian cubature, where it has been theoretically or empirically demonstrated that (i) if states are randomly sampled, the optimal sampling distribution will be $n$-dependent [3] and over-dispersed with respect to the distributional target [10], and (ii)  space-filling designs are asymptotically optimal for typical stationary kernels on bounded domains $\mathcal{X} \subset \mathbb{R}^d$ [11].
Analysis of optimal states on unbounded domains appears to be more difficult; see e.g. 35.
Relaxation of either the positivity or normalisation constraints results in approximations that behave similarly to kernel cubature <cit.>.
However, relaxation of either constraint can result in the failure of $P_n^\star$ to be an element of $\mathcal{P}(\mathcal{X})$, limiting the relevance of these results to the posterior approximation task.

Despite relatively little being known about the character of optimal states in this context, kernel discrepancy is widely used.
The application of kernel discrepancies to an implicitly defined distributional target, such as a posterior distribution in a Bayesian analysis, is made possible by the use of a Stein kernel; a $P$-dependent kernel $k = k_P$ for which $\mu_P(x) = 0$ for all $x \in \mathcal{X}$  [45].
The associated kernel discrepancy

\begin{align}
D_P(Q) = \|\mu_Q\|_{\mathcal{H}(k_P)} = \sqrt{\iint k_P(x,y) \; \mathrm{d}Q(x) \mathrm{d}Q(y)} \label{eq: KSD}
\end{align}

is called a ksd [18, 41, 28], and this will be a key tool in our methodological development.
The corresponding optimally weighted approximation $P_n^\star$ is the Stein importance sampling method of 40.
To retain clarity of presentation in the main text, we defer all details on the construction of Stein kernels to <Ref>.

§.§ Sparse Approximation

If the number $n$ of states is large, computation of optimal weights can become impractical.
This has motivated a range of sparse approximation techniques, which aim to iteratively construct an approximation of the form $P_{n,m} = \frac{1}{m} \sum_{i=1}^m  \delta(y_i)$, where each $y_i$ is an element from $\{x_1,\dots,x_n\}$.
The canonical example is the greedy algorithm which, at iteration $j$, selects a state

\begin{align}
y_j \in \argmin_{y \in \{x_1,\dots,x_n\} } D_P\left( \frac{1}{j} \delta(y) + \frac{1}{j} \sum_{i=1}^{j-1} \delta(y_i) \right)  \label{eq: greedy alg}
\end{align}

for which the statistical divergence is minimised.
In the context of kernel discrepancy, the greedy algorithm (<ref>) has computational cost $O(m^2 n)$, which compares favourably[It is difficult to quantify the complexity of numerically solving (<ref>), since this will depend on details of the solver and the tolerance that are used.
On the other hand, if we ignore the non-negativity and normalisation constraints, then we can see that the computational cost of solving the $n$-dimensional linear system of equations is $O(n^3)$.] with the cost of solving (<ref>) when $m \ll n$.
Furthermore, under appropriate assumptions, the sparse approximation converges to the optimally weighted approximation; $\smash{D_P(P_{n,m}) \rightarrow D_P(P_n^\star)}$ as $m \rightarrow \infty$ with $n$ fixed.
See 52 for full details, where non-myopic and mini-batch extensions of the greedy algorithm are also considered.
The greedy algorithm can be viewed as a regularised version of the Frank–Wolfe algorithm (also called herding, or the conditional gradient method), for which a similar asymptotic result can be shown to hold [15, 4, 14].
Related work includes 21, 22, 50.
Since in what follows we aim to retrospectively improve mcmc output, where it is not unusual to encounter $n \approx 10^4$–$10^6$, sparse approximation will be important.

This completes our overview of background material.
In what follows we seek to mimic classical quantisation by deriving a choice for $\Pi$ that is straight-forward to sample using mcmc and is appropriately over-dispersed relative to $P$.
This should be achieved while remaining in the framework of kernel discrepancies, so that optimal weights can be explicitly computed, and coupled with a sparse approximation that has low computational and storage cost.

§ METHODOLOGY

The methods that we consider first sample states $\{x_1,\dots,x_n\}$ using $\Pi$-invariant mcmc, then post-process these states using kernel discrepancies (<Ref>) and sparse approximation (<Ref>), to obtain an approximation to the target $P$.
A variational argument, which we present in <Ref>, provides a suitable $n$-independent choice for $\Pi$ (which agrees with our intuition from <Ref> that $\Pi$ should be in some appropriate sense over-dispersed with respect to $P$).
Sufficient conditions for strong consistency of the approximation are established in <Ref>.

§.§ Selecting $\Pi$

Here we present a heuristic argument for a particular choice of $\Pi$; rigorous theoretical support for Stein $\Pi$-Importance Sampling is then provided in <Ref>.
Our setting is that of <Ref>, and the following will additionally be assumed:

It is assumed that

* $C_1^2 := \inf_{x \in \mathcal{X}} k(x,x) > 0$

* $C_2 := \int \sqrt{k(x,x)} \; \mathrm{d}P(x) < \infty$.

Note that (A2) implies that $\mathcal{H}(k) \subset L^1(P)$, and thus (<ref>) is in fact a strong (or Bochner) integral.

A direct analysis of the optimal states associated to the optimal weights $w^\star$ appears to be challenging due to the fact that the components of $w^\star$ are strongly inter-dependent.
Our solution here is to instead consider optimal states associated with weights that are near-optimal and whose components are only weakly dependent.
Specifically, we will be assuming that $P$ is absolutely continuous with respect to $\Pi$ (denoted $P \ll \Pi$), and study convergence of snis, i.e. the approximation $\smash{P_n = \sum_{i=1}^n w_i \delta(x_i)}$, $w_i \propto (\mathrm{d}P/\mathrm{d}\Pi)(x_i)$, where $x_1,\dots,x_n \sim \Pi$ are independent.
Since $w \geq 0$ and $\smash{1^\top w = 1}$, from the optimality of $w^\star$ under these constraints we have that $D_P(P_n^\star) \leq D_P(P_n)$.
It is emphasised that the snis weights are a theoretical device only, and will not be used for computation; indeed, we can demonstrate that the snis weights $w$ perform substantially worse than $w^\star$ in general.

The analysis of snis weights $w$ is tractable when viewed as approximation of the kernel mean embedding $\mu_P$ in the Hilbert space $\mathcal{H}(k)$.
Indeed, recall that $D_P(P_n) = \|\xi_n / \sqrt{n}\|_{\mathcal{H}(k)}$ where $\smash{\xi_n = \sqrt{n}(\mu_{P_n} - \mu_P)}$.
Then, following Section 2.3.1 of 1, we observe that

\begin{align}
\xi_n
= \sqrt{n} \left( \sum_{i=1}^n w_i k(\cdot,x_i) - \mu_P \right) = \frac{\frac{1}{\sqrt{n}} \sum_{i=1}^n \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x_i) \left[ k(\cdot,x_i) - \mu_P \right] }{ \frac{1}{n} \sum_{i=1}^n \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x_i) } . \label{eq: Sergios frac}
\end{align}

The idea is to seek $\Pi$ for which the asymptotic variance of $\xi_n$ is small.
Supposing that

\begin{align*}
\int \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x)^2 \; \mathrm{d}\Pi(x) < \infty , \tag{S1}
\end{align*}

from the weak law of large numbers the denominator in (<ref>) converges in probability to 1.
Further supposing that

\begin{align}
\int \left\| \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x) [k(\cdot,x) - \mu_P ] \right\|_{\mathcal{H}(k)}^2 \; \mathrm{d}\Pi(x) < \infty ,
\tag{S2}
\end{align}

from the Hilbert space central limit theorem the numerator in (<ref>) converges in distribution to a Gaussian
$\smash{\frac{1}{\sqrt{n}} \sum_{i=1}^n (\mathrm{d}P / \mathrm{d}\Pi)(x_i) \left[ k(\cdot,x_i) - \mu_P \right] \stackrel{\text{d}}{\rightarrow} \mathcal{N}(0,\mathcal{C})}$ where $\mathcal{C} : \mathcal{H}(k) \rightarrow \mathcal{H}(k)$ is the covariance operator defined via

\langle f , \mathcal{C} g \rangle_{\mathcal{H}(k)} = \int \left\langle f , \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x) [ k(\cdot,x) - \mu_P ] \right\rangle_{\mathcal{H}(k)} \left\langle g , \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x) [ k(\cdot,x) - \mu_P ] \right\rangle_{\mathcal{H}(k)} \; \mathrm{d}\Pi(x) ,

see Section 10.1 of 39.
Thus, from Slutsky's lemma applied to (<ref>), we conclude that $\smash{\xi_n \stackrel{\text{d}}{\rightarrow} \mathcal{N}(0,\mathcal{C})}$.
Recalling that $\smash{n D_P(P_n)^2 = \|\xi_n\|_{\mathcal{H}(k)}^2}$, and noting that the mean square of the limiting Gaussian random variable is $\text{tr}(\mathcal{C})$, a natural idea is to select the sampling distribution $\Pi$ such that $\text{tr}(\mathcal{C})$ is minimised.

Fortunately, the trace of $\mathcal{C}$ can be explicitly computed.
It simplifies presentation to restrict attention to a Stein kernel $k = k_P$, for which $\mu_P = 0$, giving
$\text{tr}(\mathcal{C}) = \int (\mathrm{d}P / \mathrm{d}\Pi)(x)^2 k_P(x) \; \mathrm{d}\Pi(x)$, where for convenience we have let $k_P(x) := k_P(x,x)$.
Assuming that $P$ and $\Pi$ admit densities $p$ and $\pi$ on $\mathcal{X} = \mathbb{R}^d$, the variational problem we wish to solve is

\begin{align}
\argmin_{\pi \in \mathcal{Q}} \int \frac{p(x)^2}{\pi(x)} k_P(x) \; \mathrm{d}x \quad \text{s.t.} \quad \int \pi(x) \; \mathrm{d}x = 1,
\label{eq: variational}
\end{align}

where $\mathcal{Q}$ be the set of positive measures on $\mathbb{R}^d$ for which (S1-2) are satisfied.
To solve this problem, we first relax the constraints (S1-2) and solve the relaxed problem using the Euler–Lagrange equations, which yield

\begin{align}
\pi(x) \propto p(x) \sqrt{k_P(x)} . \label{eq: optimal Q}
\end{align}

Note that the normalisation constant of $\pi$ is $C_2$ from (<ref>), whose existence we assumed.
Then we verify that (S1-2) in fact hold for this choice of $\Pi$.

\begin{align*}
\text{(S1)} & = \int \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x)^2 \; \mathrm{d}\Pi(x) = C_2 \int \frac{1}{k_P(x)} \; \mathrm{d}\Pi(x) \leq \frac{C_2}{C_1^2} < \infty \\
\text{(S2)} & = \int \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x)^2 k_P(x) \; \mathrm{d}\Pi(x) = C_2 \int \sqrt{k_P(x)} \; \mathrm{d}P(x) = C_2^2 < \infty ,
\end{align*}

which shows that we have indeed solved (<ref>).
The sampling distribution $\Pi$ we have obtained is characterised up to a normalisation constant in (<ref>), so just like $P$ we can sample from $\Pi$ using techniques such as mcmc.
The Stein kernel $k_P$ determines the extent to which $\Pi$ differs from $P$, as we illustrate next.

§.§ Illustration

Illustrating our choice of $\Pi$ in 1D.  (a)  The univariate target $P$ (black), and our choice of $\Pi$ based on the Langevin–Stein kernel (purple), the KGM1–Stein kernel (green), and the KGM3–Stein kernel (blue).  (b)  The mean kernel Stein discrepancy (KSD) for Stein $\Pi$-Importance Sampling using the Stein kernels from (a); in each case, KSD was computed using the same Stein kernel used to construct $\Pi$.
Solid lines indicate the baseline case of sampling from $P$, while dashed lines indicate sampling from $\Pi$.
(The experiment was repeated 100 times and standard error bars are plotted.)

For illustration, consider the univariate target $P$ (black curve) in <Ref>, a 3-component Gaussian mixture model.
Our recommended choice of $\Pi$ in (<ref>) is shown for both the Langevin–Stein kernel (purple curve) and the KGM$s$–Stein kernels with $s \in \{1,3\}$ (green and blue curves).
The Stein discrepancy corresponding to the Langevin–Stein kernel provides control over weak convergence (i.e. convergence of integrals of functions that are continuous and bounded), while the KGM$s$–Stein kernel provides additional control over the convergence of polynomial moments up to order $s$; full details about the construction of Stein kernels are contained in <Ref>.
The Langevin and KGM1–Stein kernels have $k_P(x) \asymp x^2$, while the KGM3–Stein kernel has $k_P(x) \asymp x^6$, in each case as $|x| \rightarrow \infty$, and thus greater over-dispersion results from use of the KGM3–Stein kernel.
This over-dispersion is less pronounced[The same holds for classical quantisation; c.f.  <Ref>.] in higher dimensions; see <Ref>.

To illustrate the performance of Stein $\Pi$-Importance Sampling, we generated a sequence $(x_n)_{n \in \mathbb{N}}$ of independent samples from $\Pi$.
For each $n \in \{1,\dots,100\}$, the samples $\{x_1,\dots,x_n\}$ were assigned optimal weights $w^\star$ by solving (<ref>), and the associated ksd was computed.
As a baseline, we performed the same calculation using independent samples from $P$.
<Ref> indicates that, for both Stein kernels, substantial improvement results from the use of samples from $\Pi$ compared to the use of samples from $P$.
Interestingly, the KGM3–Stein kernel demonstrated a larger improvement compared to the Langevin–Stein kernel, suggesting that the choice of $\Pi$ may be more critical in settings where ksd enjoys a stronger form of convergence control.

To illustrate a posterior approximation task,
consider a simple regression model $y_i = f_i(x) + \epsilon_i$ with $f_i(x) = x_1(1 + t_i x_2)$, $t_i = i - 5$, $i = 1,\dots,10$, with $\epsilon_i$ independent $\mathcal{N}(0,1)$.
The parameter $x = (x_1,x_2)$ was assigned a prior $\mathcal{N}(0,I)$.
Data were simulated using $x = (0,0)$.
The posterior distribution $P$ is depicted in the leftmost panel of <Ref>, while our choice of $\Pi$ corresponding to the Langevin (centre left), KGM3 (centre right) and Riemann–Stein kernels (right) are also displayed.
For the Langevin and KGM3 kernels, the associated $\Pi$ target their mass toward regions where $P$ varies the most.
The reason for this behaviour is clearly seen for the Langevin–Stein kernel since $k_P(x) = c_1 + c_2 \|\nabla \log p(x)\|^2$ for some $c_1, c_2 > 0$; see <Ref> for detail.
The Riemann–Stein kernel can be viewed as a preconditioned form of the Langevin–Stein kernel which takes into account the geometric structure of $P$; see <Ref> for full detail[The use of geometric information may be beneficial, in the sense that the associated diffusion process may mix more rapidly, and rapid mixing leads to sharper bounds from the perspective of convergence control [26].
However, the Riemann–Stein kernel is associated with a prohibitive computational cost; it is included here only for academic interest.].
Results in <Ref> demonstrate that Stein $\Pi$-Importance Sampling improves upon the default Stein importance sampling method (i.e. with $\Pi$ and $P$ equal) for all choices of kernel.

An additional illustration involving a GARCH model with $d=4$ parameters is presented in <Ref>, where the effect of varying the order $s$ of the KGM–Stein kernel is explored.

Illustrating our choice of $\Pi$ in 2D.  The bivariate target $P$ (left), together with our choice of $\Pi$ based on the Langevin–Stein kernel (centre left), the KGM3–Stein kernel (centre right), and the Riemann–Stein kernel (right).

§.§ Theoretical Guarantees

The aim of this section is to establish when post-processing of $\Pi$-invariant mcmc produces a strongly consistent approximation of $P$, for our recommended choice of $\Pi$ in (<ref>).
Our analysis focuses on mala [48], leveraging the recent work of [20] to present explicit and verifiable conditions on $P$ for our results to hold.
In fact, we consider the more general preconditioned form of mala, where the symmetric positive definite preconditioner matrix $M$ is to be specified.
Our results also allow for (optional) sparse approximation, to circumvent direct solution of (<ref>) (c.f. <Ref>).
The resulting algorithms, which we call Stein $\Pi$-Importance Sampling (S$\Pi$IS-MALA) and Stein $\Pi$-Thinning (S$\Pi$T-MALA), are quite straight-forward and contained, respectively, in <Ref>.

[1]
$x_0$ (initial state), $\epsilon$ (step size), $M$ (preconditioner matrix), $n$ (chain length), $k_P$ (Stein kernel)
$\textstyle x' \gets \underbrace{ \textstyle x_{i-1} + \epsilon M^{-1} \nabla \log p(x_{i-1}) + \frac{\epsilon}{2} M^{-1} \nabla \log k_P(x_{i-1}) }_{=: \nu(x_{i-1})} + \sqrt{2 \epsilon} M^{-1/2} Z_i$   $Z_i \stackrel{\text{IID}}{\sim} \mathcal{N}(0,I)$
$\textstyle L \gets \log\left(\frac{p(x')}{p(x_{i-1})}\right) + \frac{1}{2} \log\left(\frac{k_P(x')}{k_P(x_{i-1})}\right) - \frac{1}{4 \epsilon} \|x_{i-1} - \nu(x')\|_{M^{-1}}^2 + \frac{1}{4 \epsilon} \| x' - \nu(x_{i-1}) \|_{M^{-1}}^2 $

if $\log(U_i) < L$ then $x_i \gets x'$; else $x_i \gets x_{i-1}$; end if
$U_i \stackrel{\text{IID}}{\sim} \mathcal{U}([0,1])$

$\Pi$-Invariant Metropolis-Adjusted Langevin Algorithm (MALA)

[1]
$\{x_1,\dots,x_n\}$ from <Ref>, $k_P$ (Stein kernel)
$w^\star \in \argmin_{w \in \mathbb{R}^d} \{ \langle w , K_P w \rangle \; : \; w \geq 0, \; 1^\top w = 1 \} $  $[K_P]_{i,j} = k_P(x_i,x_j)$

Stein $\Pi$-Importance Sampling (S$\Pi$IS-MALA)

[1]
$\{x_1,\dots,x_n\}$ from <Ref>, $m$ (number of samples to retain), $k_P$ (Stein kernel)
$\textstyle y_i \gets \argmin_{y \in \{x_1,\dots,x_n\}} \frac{1}{2} k_P(y) + \sum_{j=1}^{i-1} k_P(y,y_j)$

Stein $\Pi$-Thinning (S$\Pi$T-MALA)

Let $A \preceq B$ indicate that $A - B$ is a positive semi-definite matrix for $A,B \in \mathbb{R}^{n \times n}$.
For a symmetric positive definite matrix $A$ let $\smash{\|z\|_A := \sqrt{z^\top A^{-1} z} }$ for $\smash{z \in \mathbb{R}^d}$.
Let $\smash{C^s(\mathbb{R}^d)}$ denote the set of $s$-times continuously differentiable real-valued functions on $\mathbb{R}^d$.

Let <Ref> hold where $k = k_P$ is a Stein kernel, and let $D_P : \mathcal{X} \times \mathcal{X} \rightarrow [0,\infty]$ denote the associated ksd.
Assume also that

* $\nabla \log p \in C^2(\mathbb{R}^d)$ with $\sup_{x \in \mathbb{R}^d} \|\nabla^2 \log p(x)\| < \infty$

* $\exists$ $b_1 > 0, B_1 \geq 0$ such that $- \nabla^2 \log p(x) \succeq b_1 I$ for all $\|x\| \geq B_1$

* $k_P \in C^2(\mathbb{R}^d)$.

* $\exists$ $0 < b_2 < 2 b_1 C_1^2$, $B_2 \geq 0$ such that $\nabla^2 k_P(x) \preceq b_2 I$ for all $\|x\| \geq B_2$

Let $P_n^\star = \sum_{i=1}^n w_i^\star \delta(x_i)$ be the result of running <Ref> and let $\smash{P_{n,m} = \frac{1}{m} \sum_{i=1}^m  \delta(y_i)}$ be the result of running <Ref>.
Let $m \leq n$ and $\smash{m = \Omega((\log n)^{\delta})}$ for some $\delta > 2$.
Then there exists $\epsilon_0 > 0$ such that, for all step sizes $\epsilon \in (0,\epsilon_0)$ and all initial states $x_0 \in \mathbb{R}^d$, $\smash{D_P(P_n^\star)\rightarrow 0 , D_P (P_{n,m}) \rightarrow 0}$ almost surely as $m,n \rightarrow \infty$.

The proof is in <Ref>.
Compared to earlier authors[These earlier results required high-level assumptions on the convergence of mcmc, for which explicit sufficient conditions had yet to be derived.], such as 13, 46, a major novelty here is that our assumptions are explicit and can often be verified  <cit.>.
ass: p concave is strong log-concavity of $P$ when $B_1=0$, while for $B_1 > 0$ this condition is slightly stronger than the related distant dissipativity condition assumed in earlier work [28, 46].
ass: kP convex holds for the Langevin–Stein kernel (i.e. weak convergence control) and for the KGM1–Stein kernel (i.e. weak convergence control + control over first moments), but not for the higher-order KGM–Stein kernels.
Extending our proof strategy to the higher-order KGM–Stein kernels would require further research into the convergence properties of mala, and this is expected to be difficult.

§ BENCHMARKING ON

The area of Bayesian computation has historically lacked a common set of benchmark problems, with classical examples being insufficiently difficult and case-studies being hand-picked [17].
To introduce objectivity into our assessment, we exploited the recently released  benchmark  [43].
This project is an attempt toward standardised benchmarking, consisting of a collection of posteriors to be numerically approximated.
Here, we systematically compared the performance of S$\Pi$IS-MALA against the default Stein importance sampling algorithm (i.e. $\Pi = P$; denoted SIS-MALA), and also against unprocessed $P$-invariant mala (i.e. uniform weights), reporting results across the breadth of .
The test problems in  are defined in the  probabilistic programming language, and so  [49] was used to directly access posterior densities and their gradients as required.
For all instances of mala, an adaptive algorithm was used to learn a suitable preconditioner matrix $M$ during the warm-up period; see <Ref>.
All experiments that we report can be reproduced using code available at <https://github.com/congyewang/Stein-Pi-Importance-Sampling>.

3|c||Langevin–Stein Kernel     3|c|KGM3–Stein Kernel

Task     $d$     MALA     SIS -MALA     S$\Pi$IS -MALA     MALA     SIS -MALA     S$\Pi$IS -MALA

earnings-earn_height     3     1.41     0.0674     0.0332     5.33     0.656     0.181

gp_pois_regr-gp_regr     3     0.298     0.0436     0.0373     1.22     0.385     0.223

kidiq-kidscore_momhs     3     1.04     0.109     0.0941     4.66     0.848     0.476

kidiq-kidscore_momiq     3     5.03     0.516     0.358     25.3     4.86     1.55

mesquite-logmesquite_logvolume     3     1.10     0.179     0.156     4.97     1.70     0.844

arma-arma11     4     4.47     1.09     1.01     26.0     8.91     6.03

earnings-logearn_logheight_male     4     9.46     1.96     1.59     53.9     15.4     8.65

garch-garch11     4     0.543     0.159     0.130     4.70     1.16     1.01

kidiq-kidscore_momhsiq     4     5.21     0.982     0.897     29.3     7.25     5.05

earnings-logearn_interaction_z     5     3.09     1.36     1.33     19.3     10.4     8.94

kidiq-kidscore_interaction     5     7.74     1.65     1.79     47.8     13.2     10.1

kidiq_with_mom_work-kidscore_interaction_c     5     1.35     0.659     0.711     7.92     4.05     4.17

kidiq_with_mom_work-kidscore_interaction_c2     5     1.38     0.689     0.699     8.09     4.24     4.25

kidiq_with_mom_work-kidscore_interaction_z     5     1.11     0.500     0.499     6.62     2.63     3.25

kidiq_with_mom_work-kidscore_mom_work     5     1.07     0.507     0.545     6.70     2.63     3.04

low_dim_gauss_mix-low_dim_gauss_mix     5     5.51     1.87     1.76     37.5     14.7     11.3

mesquite-logmesquite_logva     5     1.83     0.821     0.818     12.6     5.73     5.59

hmm_example-hmm_example     6     1.99     0.578     0.523     11.6     4.13     3.40

sblrc-blr     6     479     154     134     3300     1100     854

sblri-blr     6     201     66.7     60.3     1340     514     595

arK-arK     7     6.87     3.39     3.16     60.4     26.4     23.0

mesquite-logmesquite_logvash     7     1.89     1.18     1.23     15.5     8.88     10.1

bball_drive_event_0-hmm_drive_0     8     1.15     0.679     0.698     8.55     4.72     3.99

bball_drive_event_1-hmm_drive_1     8     42.9     11.9     12.4     285     85.6     67.8

hudson_lynx_hare-lotka_volterra     8     4.62     2.29     2.15     47.4     18.8     18.9

mesquite-logmesquite     8     1.46     1.00     1.06     13.3     8.28     9.14

mesquite-logmesquite_logvas     8     2.02     1.31     1.35     19.2     10.8     12.2

mesquite-mesquite     8     0.429     0.268     0.235     3.71     2.17     2.42

eight_schools-eight_schools_centered     10     0.526     0.100     0.182     7.53     2.15     215

eight_schools-eight_schools_noncentered     10     0.210     0.137     0.137     43.6     28.7     27.5

nes1972-nes     10     6.16     3.89     3.45     72.9     36.2     34.4

nes1976-nes     10     6.67     3.86     3.53     77.5     35.5     34.4

nes1980-nes     10     4.34     2.68     2.57     49.8     25.4     25.7

nes1984-nes     10     6.18     3.75     3.43     71.3     34.9     33.6

nes1988-nes     10     7.40     3.70     3.27     81.4     34.6     32.4

nes1992-nes     10     7.52     4.32     3.84     89.1     39.7     37.3

nes1996-nes     10     6.44     3.87     3.53     74.1     36.4     34.3

nes2000-nes     10     3.35     2.22     2.20     38.6     21.3     22.8

diamonds-diamonds     26     196     157     143     5120     2990     2620

mcycle_gp-accel_gp     66     11.3     8.25     9.79     960     623     815

Benchmarking on .
Here we compared raw output from MALA with the post-processed output provided by the default Stein importance sampling method of 40 (SIS-MALA) and the proposed Stein $\Pi$-Importance Sampling method (S$\Pi$IS-MALA).
Here $d = \text{dim}(P)$ and the number of MALA samples was $n=3 \times 10^3$.
The Langevin and KGM3–Stein kernels were used for SIS-MALA and S$\Pi$IS-MALA and the associated ksd are reported.
Ten replicates were computed and statistically significant improvement is highlighted in bold.

Results are reported in <Ref> for $n = 3 \times 10^3$ samples from mala.
These focus on the Langevin–Stein kernel, for which our theory holds, and the KGM3–Stein kernel, for which it does not.
There was a significant improvement of S$\Pi$IS-MALA over SIS-MALA in 73% of test problems for the Langevin–Stein kernel and in 65% of test problems for the KGM3–Stein kernel.
Compared to unprocessed mala, a significant improvement occurred in 100% and 97% of cases, respectively for each kernel.
However, the extent of improvement decreased when the dimension $d$ of the target increased, supporting the intuition that we set out earlier and in <Ref>.
An in-depth breakdown of results, including varying the number $n$ of samples that were used, and the performance S$\Pi$T-MALA, can be found in <Ref>.

Our focus is on the development of algorithms for minimisation of ksd; the properties of ksd themselves are out of scope for this work[The interested reader is referred to 28, 8, 34.].
Nonetheless, there is much interest in better understanding the properties of ksd, and we therefore also report performance of S$\Pi$IS-MALA in terms of 1-Wasserstein divergence in <Ref>.
The main contrast between these results and the results in <Ref> is that, being score-based, ksd suffer from the blindness to mixing proportions phenomena which has previously been documented in 53, 38, 42.
Caution should therefore be taken when using algorithms based on Stein discrepancies in the context of posterior distributions with multiple high probability regions that are spatially separated.
This is also a failure mode for mcmc algorithms such as mala, and yet there are still many problems for which mala has been successfully used.

The alternative choice $\Pi_1$, with $\smash{\pi_1(x) \propto p(x)^{d/(d+1)}}$, which provides a generic form of over-dispersion and is optimal for approximation in 1-Wasserstein divergence (c.f. <Ref>), was also considered.
Results in <Ref> indicate that, while $\Pi_1$ yields an improvement compared to the baseline of using $P$ itself, $\Pi_1$ may be less effective than our proposed $\Pi$ when $P$ is skewed.

§ DISCUSSION

This paper presented Stein $\Pi$-Importance Sampling; an algorithm that is simple to implement, admits an end-to-end theoretical treatment, and achieves a significant improvement over existing post-processing methods based on ksd.
On the negative side, second order derivatives of the statistical model are required.
For models for which access to second order derivatives is impractical, our methodology and theoretical analysis are directly applicable to gradient-free ksd [24], and this would be an interesting direction for future work.

CW was supported by the China Scholarship Council.
HK and CJO were supported by EP/W019590/1.

[1]
Agapiou, S., Papaspiliopoulos, O., Sanz-Alonso, D., and Stuart, A. M. (2017).
Importance sampling: Intrinsic dimension and computational cost.
Statistical Science, 32(3):405–431.

[2]
Anastasiou, A., Barp, A., Briol, F.-X., Ebner, B., Gaunt, R. E., Ghaderinezhad,
F., Gorham, J., Gretton, A., Ley, C., Liu, Q., Mackey, L., Oates, C. J.,
Reinert, G., and Swan, Y. (2023).
Stein's method meets statistics: A review of some recent
Statistical Science, (38):120–139.

[3]
Bach, F. (2017).
On the equivalence between kernel quadrature rules and random feature
The Journal of Machine Learning Research, 18(1):714–751.

[4]
Bach, F., Lacoste-Julien, S., and Obozinski, G. (2012).
On the equivalence between herding and conditional gradient
In Proceedings of the 29th International Conference on Machine

[5]
Barbour, A. D. (1988).
Stein's method and poisson process convergence.
Journal of Applied Probability, 25(A):175–184.

[6]
Barbour, A. D. (1990).
Stein's method for diffusion approximations.
Probability Theory and Related Fields, 84(3):297–322.

[7]
Barp, A., Oates, C. J., Porcu, E., and Girolami, M. (2022a).
A Riemann–Stein kernel method.
Bernoulli, 28(4):2181–2208.

[8]
Barp, A., Simon-Gabriel, C.-J., Girolami, M., and Mackey, L. (2022b).
Targeted separation and convergence with kernel discrepancies.
arXiv preprint arXiv:2209.12835.

[9]
Bénard, C., Staber, B., and Da Veiga, S. (2023).
Kernel stein discrepancy thinning: A theoretical perspective of
pathologies and a practical fix with regularization.
arXiv preprint arXiv:2301.13528.

[10]
Briol, F.-X., Oates, C. J., Cockayne, J., Chen, W. Y., and Girolami, M. (2017).
On the sampling problem for kernel quadrature.
In Proceedings of the 34th International Conference on Machine
Learning, pages 586–595.

[11]
Briol, F.-X., Oates, C. J., Girolami, M., Osborne, M. A., and Sejdinovic, D.
Probabilistic integration: A role in statistical computation (with
discussion and rejoinder).
Statistical Science, 34(1):1–22.

[12]
Carmeli, C., De Vito, E., and Toigo, A. (2006).
Vector valued reproducing kernel Hilbert spaces of integrable
functions and mercer theorem.
Analysis and Applications, 4(04):377–408.

[13]
Chen, W. Y., Barp, A., Briol, F.-X., Gorham, J., Girolami, M., Mackey, L., and
Oates, C. J. (2019).
Stein point Markov chain Monte Carlo.
In Proceedings of the 36th International Conference on Machine
Learning, pages 1011–1021.

[14]
Chen, W. Y., Mackey, L., Gorham, J., Briol, F.-X., and Oates, C. J. (2018).
Stein points.
In Proceedings of the 35th International Conference on Machine
Learning, pages 844–853.

[15]
Chen, Y., Welling, M., and Smola, A. (2010).
Super-samples from kernel herding.
In Proceedings of the 26th Conference on Uncertainty in
Artificial Intelligence, pages 109–116.

[16]
Chopin, N. and Ducrocq, G. (2021).
Fast compression of MCMC output.
Entropy, 23(8):1017.

[17]
Chopin, N. and Ridgway, J. (2017).
Leave pima indians alone: Binary regression as a benchmark for
Bayesian computation.
Statistical Science, 32(1):64–87.

[18]
Chwialkowski, K., Strathmann, H., and Gretton, A. (2016).
A kernel test of goodness of fit.
In Proceedings of the 33rd International Conference on Machine
Learning, pages 2606–2615.

[19]
Cohort, P. (2004).
Limit theorems for random normalized distortion.
The Annals of Applied Probability, 14(1):118–143.

[20]
Durmus, A. and Moulines, É. (2022).
On the geometric convergence for MALA under verifiable conditions.
arXiv preprint arXiv:2201.01951.

[21]
Dwivedi, R. and Mackey, L. (2021).
Kernel thinning.
In Proceedings of 34th Conference on Learning Theory, pages

[22]
Dwivedi, R. and Mackey, L. (2022).
Generalized kernel thinning.
In Proceedings of the 10th International Conference on Learning

[23]
Ehler, M., Gräf, M., and Oates, C. J. (2019).
Optimal Monte Carlo integration on closed manifolds.
Statistics and Computing, 29(6):1203–1214.

[24]
Fisher, M. A. and Oates, C. J. (2022).
Gradient-free kernel Stein discrepancy.
arXiv preprint arXiv:2207.02636.

[25]
Girolami, M. and Calderhead, B. (2011).
Riemann manifold Langevin and Hamiltonian Monte Carlo
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 73(2):123–214.

[26]
Gorham, J., Duncan, A. B., Vollmer, S. J., and Mackey, L. (2019).
Measuring sample quality with diffusions.
The Annals of Applied Probability, 29(5):2884–2928.

[27]
Gorham, J. and Mackey, L. (2015).
Measuring sample quality with Stein's method.
In Proceedings of the 29th Conference on Neural Information
Processing Systems, pages 226–234.

[28]
Gorham, J. and Mackey, L. (2017).
Measuring sample quality with kernels.
In Proceedings of the 34th International Conference on Machine
Learning, pages 1292–1301.

[29]
Gotze, F. (1991).
On the rate of convergence in the multivariate CLT.
The Annals of Probability, pages 724–739.

[30]
Graf, S. and Luschgy, H. (2007).
Foundations of Quantization for Probability Distributions.

[31]
Hawkins, C., Koppel, A., and Zhang, Z. (2022).
Online, informative MCMC thinning with kernelized Stein
arXiv preprint arXiv:2201.07130.

[32]
Hayakawa, S., Oberhauser, H., and Lyons, T. (2022).
Positively weighted kernel quadrature via subsampling.
In Proceedings of the 35th Conference on Neural Information
Processing Systems.

[33]
Hodgkinson, L., Salomone, R., and Roosta, F. (2020).
The reproducing Stein kernel approach for post-hoc corrected
arXiv preprint arXiv:2001.09266.

[34]
Kanagawa, H., Gretton, A., and Mackey, L. (2022).
Controlling moments with kernel Stein discrepancies.
arXiv preprint arXiv:2211.05408.

[35]
Karvonen, T., Oates, C. J., and Girolami, M. (2021).
Integration in reproducing kernel Hilbert spaces of Gaussian
Mathematics of Computation, 90(331):2209–2233.

[36]
Karvonen, T., Oates, C. J., and Sarkka, S. (2018).
A Bayes–Sard cubature method.
In Proceedings of the 32nd Conference on Neural Information
Processing Systems.

[37]
Kent, J. (1978).
Time-reversible diffusions.
Advances in Applied Probability, 10(4):819–835.

[38]
Koehler, F., Heckett, A., and Risteski, A. (2022).
Statistical efficiency of score matching: The view from isoperimetry.
In Proceedings of the 36th Conference on Neural Information
Processing Systems.

[39]
Ledoux, M. and Talagrand, M. (1991).
Probability in Banach Spaces: Isoperimetry and Processes,
volume 23.
Springer Science & Business Media.

[40]
Liu, Q. and Lee, J. (2017).
Black-box importance sampling.
In Proceedings of the 20th International Conference on
Artificial Intelligence and Statistics, pages 952–961.

[41]
Liu, Q., Lee, J., and Jordan, M. (2016).
A kernelized Stein discrepancy for goodness-of-fit tests.
In Proceedings of the 33rd International Conference on Machine
Learning, pages 276–284.

[42]
Liu, X., Duncan, A., and Gandy, A. (2023).
Using perturbation to improve goodness-of-fit tests based on
kernelized Stein discrepancy.
In Proceedings of the 40th International Conference on Machine

[43]
Magnusson, M., Bürkner, P., and Vehtari, A. (2022).
: A set of posteriors for Bayesian inference
and probabilistic programming.

[44]
Meyn, S. P. and Tweedie, R. L. (2012).
Markov Chains and Stochastic Stability.
Springer Science & Business Media.

[45]
Oates, C. J., Girolami, M., and Chopin, N. (2017).
Control functionals for Monte Carlo integration.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 79(3):695–718.

[46]
Riabiz, M., Chen, W., Cockayne, J., Swietach, P., Niederer, S. A., Mackey, L.,
and Oates, C. J. (2022).
Optimal thinning of MCMC output.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 84(4):1059–1081.

[47]
Roberts, G. O. and Rosenthal, J. S. (1998).
Optimal scaling of discrete approximations to Langevin diffusions.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 60(1):255–268.

[48]
Roberts, G. O. and Stramer, O. (2002).
Langevin diffusions and Metropolis–Hastings algorithms.
Methodology and Computing in Applied Probability, 4:337–357.

[49]
Roualdes, E., Ward, B., Axen, S., and Carpenter, B. (2023).
BridgeStan: Efficient in-memory access to Stan programs through
Python, Julia, and R.

[50]
Shetty, A., Dwivedi, R., and Mackey, L. (2022).
Distribution compression in near-linear time.
In Proceedings of the 10th International Conference on Learning

[51]
Sonnleitner, M. (2022).
The Power of Random Information for Numerical Approximation and
PhD thesis, University of Passau.

[52]
Teymur, O., Gorham, J., Riabiz, M., and Oates, C. J. (2021).
Optimal quantisation of probability measures using maximum mean
In Proceedings of the 24th International Conference on
Artificial Intelligence and Statistics, pages 1027–1035.

[53]
Wenliang, L. K. and Kanagawa, H. (2021).
Blindness of score-based methods to isolated components and mixing
In Proceedings of “Your Model is Wrong” @ the 35th Conference
on Neural Information Processing Systems.

This supplement contains supporting material for the paper Stein $\Pi$-Importance Sampling.
The mathematical background on Stein kernels is contained in <Ref>.
The proof of <Ref> is contained in <Ref>.
For implementation of Stein $\Pi$-Importance Sampling without the aid of automatic differentiation, various explicit derivatives are required; the relevant calculations can be found in <Ref>.
The empirical protocols and additional empirical results are presented in <Ref>.

§ MATHEMATICAL BACKGROUND

This appendix contains mathematical background on reproducing kernels and Stein kernels, as used in the main text.
<Ref> introduces matrix-valued reproducing kernels, while <Ref> specialises to Stein kernels by application of a Stein operator to a matrix-valued kernel.
A selection of useful Stein kernels are presented in <Ref>.

§.§ Matrix-Valued Reproducing Kernels

A matrix-valued kernel is a function $K : \mathbb{R}^d \times \mathbb{R}^d \rightarrow \mathbb{R}^{d \times d}$, that is both

* symmetric; $K(x,y) = K(y,x)$ for all $x,y \in \mathbb{R}^d$, and

* positive semi-definite; $\sum_{i=1}^n \sum_{j=1}^n \langle c_i , K(x_i,x_j) c_j \rangle \geq 0$ for all $x_1,\dots,x_n \in \mathbb{R}^d$ and $c_1,\dots,c_n \in \mathbb{R}^d$.

Let $K_x = K(\cdot,x)$.
For vector-valued functions $g,g' : \mathbb{R}^d \rightarrow \mathbb{R}^d$, defined by $g = \sum_{i=1}^n K_{x_i} c_i$ and $g' = \sum_{j=1}^m K_{x_j'} c_i'$, define an inner product

\begin{align}
\langle g , g' \rangle_{\mathcal{H}(K)} = \sum_{i=1}^n \sum_{j=1}^m \langle c_i , K(x_i , x_j') c_j' \rangle . \label{eq: HK inn prod}
\end{align}

There is a unique Hilbert space of such vector-valued functions associated to $K$, denoted $\mathcal{H}(K)$; see Proposition 2.1 of 12.
This space is characterised as

\begin{align*}
\mathcal{H}(K) = \overline{\mathrm{span}}\{ K_x c : x, c \in \mathbb{R}^d \}
\end{align*}

where here the closure is taken with respect to the inner product in (<ref>).
It can be shown that $\mathcal{H}(K)$ is in fact a rkhs which satisfies the reproducing property

\begin{align*}
\langle g , K_x c \rangle_{\mathcal{H}(K)} = \langle g(x) , c \rangle
\end{align*}

for all $g \in \mathcal{H}(K)$ and $x, c \in \mathbb{R}^d$.
Matrix-valued kernels are the natural starting point for construction of ksd, as described next.

§.§ Stein Kernels

A general construction for Stein kernels is to first identify a matrix-valued rkhs $\mathcal{H}(K)$ and an operator $S_P : \mathcal{H}(K) \rightarrow L^1(P)$ for which $\int S_p h \; \mathrm{d}P = 0$ for all $h \in \mathcal{H}(K)$.
Such an operator will be called a Stein operator.
The collection $\{S_p h : h \in \mathcal{H}(K)\}$ inherits the structure of an rkhs, whose reproducing kernel

\begin{align}
k_P(x,y) = \langle S_P K_x , S_P K_y \rangle_{\mathcal{H}(K)}  \label{eq: stein kernel}
\end{align}

is a Stein kernel, meaning that $\mu_P = 0$ where $\mu_P$ is the kernel mean embedding from (<ref>); see 8.
Explicit calculations for the Stein kernels considered in this work can be found in <Ref>.

For univariate distributions, 5 proposed to obtain Stein operators from infinitessimal generators of $P$-invariant continuous-time Markov processes; see also 6, 29.
The approach was extended to multivariate distributions in 27.
The starting point is the $P$-invariant Itô diffusion

\begin{align}
\mathrm{d}X_t = \frac{1}{2} \frac{1}{p(X_t)} \nabla \cdot [ p(X_t) M(X_t) ] \mathrm{d}t + M(X_t)^{1/2} \mathrm{d}W_t ,
\label{eq: ito}
\end{align}

where $p$ is the density of $P$, assumed to be positive, $M : \mathbb{R}^d \rightarrow \mathbb{R}^{d \times d}$ is a symmetric matrix called the diffusion matrix, and $W_t$ is a standard Wiener process [37, 48].
Here the notation $[\nabla \cdot A]_i = \nabla \cdot (A_{i,\cdot}^\top)$ indicates the divergence operator applied to each row of the matrix $\smash{A(x) \in \mathbb{R}^{d \times d}}$.
The infinitessimal generator is

\begin{align*}
(A_P u)(x) = \frac{1}{2} \frac{1}{p(x)} \nabla \cdot [p(x) M(x) \nabla u(x)] .
\end{align*}

Substituting $h(x)$ for $\frac{1}{2} \nabla u(x)$, we obtain a Stein operator

\begin{align}
(S_P h)(x) = \frac{1}{p(x)} \nabla \cdot [p(x) M(x) h(x)] \label{eq: diff stein op}
\end{align}

called the diffusion Stein operator [26].
This is indeed a Stein operator, since under mild integrability conditions on $K$, the divergence theorem gives that $\int S_p h \; \mathrm{d}P = 0$ for all $h \in \mathcal{H}(K)$; for full details and a proof see 8.

§.§ Selecting a Stein Kernel

There are several choices for a Stein kernel, and which we should use depends on what form of convergence we hope to control [28, 26, 33, 8, 34].
<Ref> describes the Langevin–Stein kernel for weak convergence control, <Ref> describes the KGM–Stein kernels for additional control over moments, and <Ref> presents the Riemann–Stein kernel, whose convergence properties have to-date been less well-studied.

All of the kernels that we consider have length scale parameters that need to be specified, and some also have location parameters to be specified.
As a reasonably automatic default we define

x_\star \in \argmax p(x) , \qquad
\Sigma^{-1} = - \nabla^2 \log p(x_\star)

as a location and a matrix of characteristic length scales for $P$ that will be used throughout.
These values can typically be obtained using gradient-based optimisation, which is usually cheaper to perform compared to full approximation of $P$.
It is assumed that $\nabla^2 \log p(x_\star)$ is positive definite in the sequel.

§.§.§ Weak Convergence Control with Langevin–Stein Kernels

The first kernel we consider, which we called the Langevin–Stein kernel in the main text, was introduced by 28.
This Stein kernel was developed for the purpose of controlling the weak convergence of a sequence $(Q_n)_{n \in \mathbb{N}} \subset \mathcal{P}(\mathbb{R}^d)$ to $P$.
Recall that a sequence $(Q_n)_{n \in \mathbb{N}}$ is said to converge weakly (or in distribution) to $P$ if $\smash{\int f \dd Q_n \to \int f \dd P}$ for all continuous bounded functions $\smash{f : \mathbb{R}^d \rightarrow \mathbb{R}}$.
This convergence is denoted $\smash{Q_n \stackrel{\text{d}}{\rightarrow} P}$ in shorthand.

The problem considered in 28 was how to select a combination of matrix-valued kernel $K$ (and, implicitly, a diffusion matrix $M$) such that the Stein kernel $k_P$ in (<ref>) generates a ksd $D_P(Q)$ in (<ref>) for which $D_P(Q_n) \rightarrow 0$ implies $\smash{Q_n \stackrel{\text{d}}{\rightarrow} P}$.
Their solution was to combine the inverse multi-quadric kernel with an identity diffusion matrix;

K(x,y) = (1 + \|x-y\|_\Sigma^2 )^{-\beta} I , \qquad M(x) = I

for $\beta \in (0,1)$.
Provided that $P$ has a density $p$ for which $\nabla \log p(x)$ is Lipschitz, and that $P$ is distantly dissipative <cit.>, the associated ksd enjoys weak convergence control.
Technically, the results in 28 apply only when $\Sigma = I$, but Theorem 4 in 13 demonstrated that they hold also for any positive definite $\Sigma$.
Following the recommendation of several previous authors, including 14, 13, 46, we take $\beta = \frac{1}{2}$ throughout.

§.§.§ Moment Convergence Control with KGM–Stein Kernels

Despite its many elegant properties, weak convergence can be insufficient for applications where we are interested in integrals $\int f \; \mathrm{d}P$ for which the integrand $f : \mathbb{R}^d \rightarrow \mathbb{R}$ is unbounded.
In particular, this is the case for moments of the form $\smash{f(x) = x_1^{\alpha_1} \dots x_d^{\alpha_d}}$, $\smash{0 \neq \alpha \in \mathbb{N}_0^d}$.
In such situations, we may seek also the stronger property of moment convergence control.
The development of ksd for moment convergence control was recently considered by 34, and we refered to their construction as the KGM–Stein kernels in the main text.
(For convenience, we have adopted the initials of the authors in naming the KGM–Stein kernel.)

A sequence $(Q_n)_{n \in \mathbb{N}} \subset \mathcal{P}(\mathbb{R}^d)$ is said to converge to $P$ in the $s$th order moment if
$\int \lVert x \rVert^s \dd Q_n (x) \to \int \lVert x \rVert^s \dd P(x)$.
To establish convergence of moments, we need an additional condition on top of weak convergence control: uniform integrability control.
A sequence of measures $(Q_n)_{n \in \mathbb{N}}$ is said to have uniformly integrable $s$th moments if for any $\varepsilon > 0$, we can take $r>0$ such that

\sup_{n \in \mathbb{N}} \; \int_{\lVert x \rVert > r} \lVert x  \rVert^s \; \dd Q_n(x) < \varepsilon.

This condition essentially states that the tail decay of the measures is well-controlled (so that it has a convergent moment).
The ksd convergence $D_P(Q_n) \to 0$ implies uniform integrability if for any $\varepsilon > 0$, we can take $r_{\varepsilon} > 0$ and $f_{\varepsilon} \in \mathcal{H}(K)$ such that

\begin{align}
S_P f_{\varepsilon} (x) \geq \lVert x \rVert^s 1\{\lVert x \rVert > r_\varepsilon\}- \varepsilon, \label{eq: Heish req}
\end{align}

i.e., the Stein-modified rkhs can approximate the (norm-weighted) indicator function arbitrarily well.
Such a function $f_{\varepsilon}$ can be explicitly constructed (while not guaranteed to be a member of the rkhs).
Specifically, the choice $f_{\varepsilon} = (1-\iota_{\varepsilon})g$ satisfies (<ref>) under an appropriate dissipativity condition, where $\iota_\varepsilon$ is a differentiable indicator function vanishing outside a ball, and $\smash{g(x)= - x / \sqrt{1 + \lVert x \rVert^2}}$.
This motivated 34 to introduce the $s$th order KGM–Stein kernel, which is based on the matrix-valued kernel and diffusion matrix

K(x,y) = \left[ \phi(\|x- y\|_\Sigma) + \kappa_{\rm lin}(x, y) \right] I , \qquad M(x) = (1 + \|x - x_\star\|_\Sigma^2)^{\frac{s-1}{2}} I ,

where $\smash{(x,y) \mapsto \phi(\|x-y\|_\Sigma)}$ is a $C_0^1$ universal kernel <cit.>.
For comparability of our results, we take $\phi$ to be the inverse multi-quadric $\smash{\phi(r) = (1 + r^2)^{-1/2}}$, and

\[ \kappa_{\rm lin}(x,y) = \frac{1 + (x - x_\star)^\top \Sigma^{-1} (y-x_\star) }{\sqrt{1 + \lVert x - x_\star \rVert_\Sigma^2} \sqrt{1 + \lVert y - x_\star \rVert_\Sigma^2} }.\]

Here the normalised linear kernel $\kappa_{\rm lin}$ ensures $g \in \mathcal{H}(K)$, while the $C_0^1$ universal kernel $\phi$ allows approximation of $S_P \iota_\varepsilon g$; see 34.

§.§.§ Exploiting Geometry with Riemann–Langevin–Stein Kernels

For academic interest only, here we describe the Riemann–Stein kernel that featured in <Ref> of the main text.
This Stein kernel is motivated by the analysis of 26, who argued that the use of rapidly mixing Itô diffusions in Stein operators can lead to sharper convergence control.
The Riemann–Stein kernel is based on the class of so-called Riemannian diffusions considered in 25, who proposed to take the diffusion matrix $M$ in (<ref>) to be $M = (\mathcal{I}_{\text{prior}} + \mathcal{I}_{\text{Fisher}})^{-1}$, the inverse of the Fisher information matrix, $\mathcal{I}_{\text{Fisher}}$, regularised using the Hessian of the negative log-prior, $\mathcal{I}_{\text{prior}}$.
For the two-dimensional illustration in <Ref>, this leads to the diffusion matrix

M(x) = \left( I + \sum_{i=1}^n [\nabla f_i(x)] [\nabla f_i(x)]^\top \right)^{-1} ,

where we recall that $y_i = f_i(x) + \epsilon_i$, where the $\epsilon_i$ are independent with $\epsilon_i \sim \mathcal{N}(0,1)$, and the prior is $x \sim \mathcal{N}(0,1)$.
For the presented experiment we paired the above diffusion matrix with the inverse multi-quadric kernel $K(x,y) = (1 + \|x-y\|_\Sigma^2)^{-\beta}$ for $\beta = \frac{1}{2}$.
The Riemann–Stein kernel extends naturally to distributions $P$ defined on Riemannian manifolds $\mathcal{X}$; see 7 and Example 1 of 33.

Unfortunately, the Riemann–Stein kernel is prohibitively expensive in most real applications, since each evaluation of $M$ requires a full scan through the size-$n$ dataset.
The computational complexity of Stein $\Pi$-Thinning with the Riemann–Stein kernel is therefore $O(m^2 n^2)$, which is unfavourable compared to the $O(m^2 n)$ complexity in the case where the Stein kernel is not data-dependent.
Furthermore, the convergence control properties of the Riemann–Stein kernel have yet to be established.
For these reasons we included the Riemann–Stein kernel for illustration only; further groundwork will be required before the Riemann-Stein kernel can be practically used.

§ PROOF OF <REF>

This appendix is devoted to the proof of <Ref>.
The proof is based on the recent work of 20, on the geometric convergence of mala, and on the analysis of sparse (greedy) approximation of kernel discrepancies performed in 46; these existing results are recalled in <Ref>.
An additional technical result on preconditioned mala is contained in <Ref>.
The proof of <Ref> itself is contained in <Ref>.

§.§ Auxiliary Results

To precisely describe the results on which our analysis is based, we first need to introduce some notation and terminology.
Let $V: \mathcal{X} \to [1,\infty)$ and, for a function $f : \mathcal{X} \rightarrow \mathbb{R}$ and a measure $\mu$ on $\mathcal{X}$, let

\|f\|_V := \sup_{x \in \mathcal{X}} \frac{|f(x)|}{V(x)}, \qquad \|\mu\|_V := \sup_{\|f\|_V \leq 1} \left| \int_{\mathcal{X}} f \mathrm{d}\mu \right| .

Recall that a $Q$-invariant Markov chain $(x_i)_{i \in \mathbb{N}} \subset \mathcal{X}$ with $n^{\text{th}}$ step transition kernel $\mathrm{Q}^n$ is $V$-uniformly ergodic <cit.> if and only if $\exists R \in [0,\infty), \rho \in (0,1)$ such that

\begin{align}
\| \mathrm{Q}^n(x,\cdot)-Q \|_{V} \leq R \rho^n V(x)  \label{eq: V unif ergod}
\end{align}

for all initial states $x \in \mathcal{X}$ and all $n \in \mathbb{N}$.

Although mala (<Ref>) is classical [48], until recently explicit sufficient conditions for ergodicity of mala had not been obtained.
The first result we will need is due 20, who presented the first explicit conditions for $V$-uniform convergence of mala.
It applies only to standard mala, meaning that the preconditioning matrix $M$ appearing in <Ref> is the identity matrix.
The extension of this result to preconditioned mala will be handled in <Ref>.

Let $Q \in \mathcal{P}(\mathbb{R}^d)$ admit a density, $q$, such that

* there exists $x_0$ with $\nabla \log q(x_0) = 0$

* $q$ is twice continuously differentiable with $\sup_{x \in \mathbb{R}^d} \| \nabla^2 \log q(x-x_0) \| < \infty$

* there exists $b > 0$ and $B \geq 0$ such that $- \nabla^2 \log q(x-x_0) \succeq b I$ for all $\|x-x_0\| \geq B$.

Then there exists $\epsilon_0 > 0$ such that for all step sizes $\epsilon \in (0,\epsilon_0)$, standard $Q$-invariant mala (i.e. with $M = I$) is $V$-uniformly ergodic for $V(x) = \exp(\frac{b}{16} \|x-x_0\|^2)$.

This is Theorem 1 of 20.

The next result that we will need establishes consistency of the greedy algorithm applied to samples from a Markov chain that is $Q$-invariant.

Let $P,Q \in \mathcal{P}(\mathcal{X})$ with $P \ll Q$.
Let $k_P : \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}$ be a Stein kernel and let $D_P : \mathcal{X} \times \mathcal{X} \rightarrow [0,\infty]$ denote the associated ksd.
Consider a $Q$-invariant, time-homogeneous Markov chain $(x_i)_{i \in \mathbb{N}} \subset \mathcal{X}$ such that

* $(x_i)_{i \in \mathbb{N}}$ is $V$-uniformly ergodic, such that $V(x) \geq \frac{\mathrm{d}P}{\mathrm{d}Q}(x) \sqrt{k_P(x)}$

* $\sup_{i \in \mathbb{N}} \mathbb{E}\left[ \frac{\mathrm{d}P}{\mathrm{d}Q}(x_i) \sqrt{k_P(x_i)} V(x_i)\right] < \infty$

* there exists $\gamma > 0$ such that
$\sup_{i \in \mathbb{N}} \mathbb{E} \left[ \exp\left\{ \gamma \max\left( 1 ,  \frac{\mathrm{d}P}{\mathrm{d}Q}(x_i)^2 \right) k_P(x_i) \right\} \right] < \infty$.

Let $P_{n,m}$ be the result of running the greedy algorithm in (<ref>).
If $m \leq n$ and $\log(n) = O(m^{\gamma/2})$ for some $\gamma < 1$, then $D_P (P_{n,m}) \rightarrow 0$ almost surely as $m,n \rightarrow \infty$.

This is Theorem 3 of 46.

§.§ Preconditioned mala

In addition to the auxiliary results in <Ref>, which concern standard mala (i.e. with $M = I$), we require an elementary fact about mala, namely that preconditioned mala is equivalent to standard mala under a linear transformation of the state variable.
Recall that the $M$-preconditioned mala algorithm is a Metropolis–Hastings algorithm whose proposal is the Euler–Maruyama discretisation of the Itô diffusion (<ref>).

Let $M(x) \equiv M$ for a symmetric positive definite and position-independent matrix $M \in \mathbb{R}^{d \times d}$.
Let $Q \in \mathcal{P}(\mathbb{R}^d)$ admit a pdf $q$ for which the $Q$-invariant  diffusion $(X_t)_{t \geq 0}$, given by setting $p=q$ in (<ref>), is well-defined.
Then under the change of variables $Y_t := M^{1/2} X_t$,

\begin{align}
\mathrm{d}Y_t = \frac{1}{2} (\nabla \log \tilde{q})(Y_t) \mathrm{d}t + \mathrm{d}W_t ,
\label{eq: transformed ito}
\end{align}

where $\tilde{q}(x) \propto q(M^{-1/2} x)$ for all $x \in \mathbb{R}^d$.

From the chain rule,

\begin{align*}
(\nabla \log \tilde{q})(y) = \nabla_y \log q(M^{-1/2} y) = M^{-1/2} (\nabla \log q)(M^{-1/2} y) ,
\end{align*}

and thus, substituting $Y_t = M^{1/2} X_t$, (<ref>) is equal to

\begin{align*}
\mathrm{d}X_t & = M^{-1/2} \left[
\frac{1}{2} M^{-1/2} (\nabla \log q)(M^{-1/2} M^{1/2} X_t) + \mathrm{d}W_t  \right] \\
& = \frac{1}{2} M^{-1} (\nabla \log q)(X_t) + M^{-1/2} \mathrm{d}W_t ,
\end{align*}

which is identical to (<ref>) in the case where $M(x) = M$ is constant.

Let $Q$ and $\tilde{Q}$ be the distributions referred to in <Ref>, whose pdf are respectively $q(x)$ and $\smash{\tilde{q}(x) \propto q(M^{-1/2} x)}$.
<Ref> then implies that the $M$-preconditioned mala algorithm applied to $Q$ (i.e. <Ref> for $\Pi = Q$) is equivalent to the standard mala algorithm (i.e. $M = I$) applied to $\tilde{Q}$.
This fact allows us to generalise the result of <Ref> as follows:

Consider a symmetric positive definite matrix $M \in \mathbb{R}^{d \times d}$.
Assume that conditions (DM1-3) in <Ref> are satisfied.
Then there exists $\epsilon_0' > 0$ and $b' > 0$ such that for all step sizes $\epsilon \in (0,\epsilon_0')$, the $M$-preconditioned $Q$-invariant mala is $V$-uniformly ergodic for $V(x) = \exp(\frac{b'}{16} \|x-x_0\|^2)$.

From <Ref> and <Ref>, the result follows if we can establish (DM1-3) for $\tilde{Q}$, since $M$-preconditioned mala is equivalent to standard mala applied to $\tilde{Q}$.
For a matrix $A \in \mathbb{R}^{d \times d}$, let $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ respectively denote the minimum and maximum eigenvalues of $A$.
For (DM1) we set $y_0 = M^{1/2} x_0$ and observe that

\begin{align*}
(\nabla \log \tilde{q})(y_0) & = M^{-1/2} (\nabla \log q)(x_0) = 0 .
\end{align*}

For (DM2) we have that

\begin{align*}
\sup_{y \in \mathbb{R}^d} \|\nabla^2 (\log \tilde{q})(y-y_0)\| & = \sup_{y \in \mathbb{R}^d} \| M^{-1/2} (\nabla^2 \log q)(M^{-1/2}(y-y_0)) M^{-1/2} \| \\
& \leq \lambda_{\min}(M)^{-1} \sup_{x \in \mathbb{R}^d} \| (\nabla^2 \log q)(x-x_0) \| < \infty .
\end{align*}

For (DM3) we have that

\begin{align*}
- (\nabla^2 \log \tilde{q})(y-y_0) & = - M^{-1/2} (\nabla^2 \log q)(M^{-1/2}(y-y_0)) M^{-1/2} \\
& = - M^{-1/2} (\nabla^2 \log q)(x-x_0) M^{-1/2}
\succeq M^{1/2} (b I) M^{1/2} = b M^{-1} \succeq b' I
\end{align*}

where $b' = b \lambda_{\max}(M)^{-1}$, which holds for all $\|x-x_0\| \geq B$, and in particular for all $\|y-y_0\| \geq B'$ where $B' = B \lambda_{\max}(M)^{1/2}$.
Thus (DM1-3) are established for $\tilde{Q}$.

The choice $M = \Sigma^{-1}$, which sets the preconditioner matrix $M$ equal to the inverse of the length scale matrix $\Sigma$ used in the specification of the kernel $K$ (c.f. <Ref>), leads to the elegant interpretation that Stein $\Pi$-Importance Sampling applied to $M$-preconditioned mala is equivalent to the Stein $\Pi$-Importance Sampling applied to standard mala (i.e. with $M=I$) for the whitened target $\tilde{P}$ with pdf $\tilde{p}(x) \propto p(M^{-1/2} x)$.
For our experiments, however, the preconditioner matrix $M$ was learned during a warm-up phase of mala, since in general the curvature of $P$ (captured by $\Sigma$) and the curvature of $\Pi$ (captured by $M^{-1}$) may be different.

§.§ Proof of <Ref>

The route to establishing <Ref> has three parts.
First, we establish (DM1-3) of <Ref> with $Q = \Pi$, to deduce from <Ref> that $\Pi$-invariant $M$-preconditioned mala is $V$-uniformly ergodic.
This in turn enables us to establish conditions (R$^+$1-3) of <Ref>, again for $Q = \Pi$, from which the strong consistency $D_P(P_{n,m}) \stackrel{\text{a.s.}}{\rightarrow} 0$ of S$\Pi$T-MALA is established.
Finally, we note that $0 \leq D_P(P_n^\star) \leq D_P(P_{n,m})$, since the support of $P_{n,m}$ is contained in the support of $P_n^\star$, and the latter is optimally weighted, whence also the strong consistency of S$\Pi$IS-MALA.

Establish (DM1-3)

First we establish (DM1-3) for $Q = \Pi$.
Fix $x_0 \in \mathbb{R}^d$.
For (DM2), first recall that the range of $k_P$ is $[C_1^2,\infty)$ where $C_1 > 0$, from <Ref>.
Since $\log(\cdot)$ has bounded second derivatives on $[C_1^2,\infty)$, there is a constant $C > 0$ such that

\begin{align*}
\forall x \in \mathbb{R}^d, \qquad \| \nabla^2 \log k_P(x) \| \leq C \| \nabla^2 k_P(x) \| .
\end{align*}

Thus, using compactness of the set $\{x : \|x-x_0\| \leq B_2\}$,

\begin{align}
\sup_{x \in \mathbb{R}^d} \| \nabla^2 \log k_P(x) \| \leq C \max\Bigg( \underbrace{ \sup_{\|x-x_0\| \leq B_2 } \| \nabla^2 k_P(x) \| }_{ < \infty \text{ by \labelcref{ass: kp} } } ,  \underbrace{ \sup_{ \|x-x_0\| \geq B_2 } \| \nabla^2 k_P(x) \| }_{ < b_2 \|I\| \text{ by
\labelcref{ass: kP convex} } } \Bigg) < \infty .
\label{eq: bounded logs}
\end{align}

Now, $\pi$ is twice differentiable as it is the product of twice differentiable functions $p$ and $\smash{k_P^{1/2}}$ from ass: p and ass: kp, and moreover

\sup_{x \in \mathbb{R}^d} \| \nabla^2 \log \pi(x-x_0) \| \leq \underbrace{ \sup_{x \in \mathbb{R}^d} \| \nabla^2 \log p(x) \| }_{< \infty \text{ by \labelcref{ass: p} } } + \frac{1}{2} \underbrace{ \sup_{x \in \mathbb{R}^d} \| \nabla^2 \log k_P(x) \| }_{ < \infty \text{ by \eqref{eq: bounded logs} } } < \infty ,

so (DM2) is satisfied.
For (DM3), first note from the chain and product rules that for all $\|x\| \geq B_2$

\begin{align}
\nabla^2 \log k_P(x-x_0) = \underbrace{ \frac{\nabla^2 k_P(x-x_0)}{k_P(x-x_0)} }_{ \preceq (b_2/C_1^2) I \text{ by \labelcref{ass: kP convex} } } - \underbrace{ \frac{[\nabla k_P(x-x_0)][\nabla k_P(x-x_0)]^\top}{k_P(x-x_0)^2} }_{ \succeq 0 } \preceq \frac{ b_2}{C_1^2} I. \label{eq: d2logk}
\end{align}

Thus, for all $\|x-x_0\| \geq B := \|x_0\| + \max(B_1,B_2)$,

\begin{align}
- \nabla^2 \log \pi(x-x_0) = \underbrace{ - \nabla^2 \log p(x-x_0) }_{ \succeq b_1 I \text{ by \labelcref{ass: p concave} } }  - \frac{1}{2} \underbrace{ \nabla^2 \log k_P(x-x_0) }_{ \preceq (b_2/C_1^2) I \text{ by \eqref{eq: d2logk} } } \succeq \underbrace{ \left( b_1 - \frac{b_2}{2 C_1^2} \right) }_{ =: b >0} I
\label{eq: q log conc}
\end{align}

as required.
The same argument establishes (DM1); from (<ref>) we have $\lim_{\|x\| \rightarrow \infty} \pi(x) = 0$, and since $\pi$ is a continuously differentiable density there must exist an $x_0$ at which $\pi$ is locally minimised.
Thus we have established (DM1-3) for $Q = \Pi$ and we may conclude from <Ref> that there is an $\epsilon_0' > 0$ and $b' > 0$ such that, for all $\epsilon \in (0,\epsilon_0')$, the $\Pi$-invariant $M$-preconditioned mala chain $(x_i)_{i \in \mathbb{N}}$ is $V$-uniformly ergodic for $V(x) = C_2 \exp(\frac{b'}{16} \|x-x_0\|^2)$ (since if a Markov chain is $V$-uniformly ergodic, then it is also $C V$-uniformly ergodic).

Establish (R$^+$1-3)

The aim is now to establish conditions (R$^+$1-3) of <Ref> for $Q = \Pi$.
By construction $\smash{\mathrm{d}P/\mathrm{d}\Pi = C_2 / \sqrt{k_P(x)} < C_2 / C_1 < \infty}$, where $C_1$ and $C_2$ were defined in <Ref>, so that $P \ll \Pi$.
It has already been established that $(x_i)_{i \in \mathbb{N}}$ is $V$-uniformly ergodic, and further

\begin{align*}
V(x) = C_2 \exp\left( \frac{b'}{16} \|x-x_0\|^2 \right) \geq C_2 = \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x) \sqrt{k_P(x)}
\end{align*}

for all $x$, which establishes (R$^+$1).
Let $R$ and $\rho$ denote constants for which the $V$-uniform ergodicity property (<ref>) is satisfied.
From $V$-uniform ergodicity, the integral $\int V \; \mathrm{d}\Pi$ exists and

\begin{align*}
\left| \mathbb{E}\left[ \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x_i) \sqrt{k_P(x_i)} V(x_i)\right] - C_2 \int V \; \mathrm{d}\Pi \right| & = C_2 \left| \mathbb{E}[V(x_i)] - \int V \; \mathrm{d}\Pi \right| \\& \leq C_2 R \rho^n V(x_0) \rightarrow 0
\end{align*}

which establishes (R$^+$2).
Fix $\gamma > 0$.
By construction $\mathrm{d}P / \mathrm{d}\Pi \leq C_2 / C_1$, and thus

\begin{align*}
\exp \left\{ \gamma \max\left( 1 , \frac{\mathrm{d}P}{\mathrm{d}\Pi}(x)^2 \right) k_P(x) \right\} < \exp \left\{ \tilde{\gamma} k_P(x) \right\}
\end{align*}

where $\tilde{\gamma} = \max(1, C_2 / C_1 ) \gamma$.
Since we have assumed that  $k_P$ is continuous with, from ass: kP convex,

\begin{align*}
b_3 := \limsup_{\|x\| \rightarrow \infty} \frac{ k_P(x) }{  \|x\|^2 } < \infty ,
\end{align*}

we may take $\gamma$ such that $\tilde{\gamma} b_3 < b' / 16$, so that $\|x \mapsto \exp\{ \tilde{\gamma} k_P(x) \} \|_V < \infty$ and in particular

\begin{align*}
\left| \mathbb{E}\left[ \exp\{ \tilde{\gamma} k_P(x_i) \} \right] - \int \exp\{ \tilde{\gamma} k_P(x) \} \; \mathrm{d}\Pi(x) \right|
& \leq \|x \mapsto \exp\{ \tilde{\gamma} k_P(x) \} \|_V \times R \rho^n V(x_0) \rightarrow 0
\end{align*}

which establishes (R$^+$3).
Thus we have established (R$^+$1-3) for $Q = \Pi$, so from <Ref> we have strong consistency of S$\Pi$T-MALA (i.e. $D_P(P_{n,m}) \stackrel{\text{a.s.}}{\rightarrow} 0$) provided that $m \leq n$ with $\smash{\log(n) = O(m^{\gamma/2})}$ for some $\gamma < 1$.
The latter condition is equivalent to $\smash{m = \Omega((\log n)^{\delta})}$ for some $\delta > 2$, which we used for the statement.
Since $0 \leq D_P(P_n^\star) \leq D_P(P_{n,m})$, the strong consistency of S$\Pi$IS-MALA is also established.

§ EXPLICIT CALCULATION OF STEIN KERNELS

This appendix contains explicit calculations for the Langevin–Stein and KGM–Stein kernels $k_P$, which are sufficient to implement Stein $\Pi$-Importance Sampling and Stein $\Pi$-Thinning.
These calculations can also be performed using automatic differentiation, but comparison to the analytic expressions is an important step in validation of computer code.

To proceed, we observe that the diffusion Stein operator $S_P$ in (<ref>) applied to a matrix-valued kernel $K$ is equivalent to the Langevin–Stein operator applied to the kernel $C(x,y) = M(x) K(x,y) M(y)^\top$.
In the case of the Langevin–Stein and KGM–Stein kernels we have $K(x,y) = \kappa(x,y) I$ for some $\kappa(x,y)$ and $M(x) = (1 + \|x - x_\star\|_\Sigma^2)^{(s-1)/2} I$ for some $s \in \{0,1,2,\dots\}$.
Thus $C(x,y) = c(x,y) I$ where

\begin{align*}
c(x,y) & := (1 + \|x - x_\star\|_\Sigma^2)^{(s-1)/2} (1 + \|y - x_\star\|_\Sigma^2)^{(s-1)/2} \kappa(x,y)
\end{align*}

\begin{align*}
k_P(x,y) & = \nabla_x \cdot \nabla_y c(x,y) + [\nabla_x c(x,y)] \cdot [\nabla_y \log p(y)] + [\nabla_y c(x,y)] \cdot [\nabla_x \log p(x)] \\
& \hspace{200pt} + c(x,y) [\nabla_x \log p(x)] \cdot [\nabla_y \log p(y)] ,
\end{align*}

following the calculations in 45.
To evaluate the terms in this formula we start by differentiating $c(x,y)$, to obtain

\begin{align*}
\nabla_x c(x,y) & = (1 + \|x - x_\star\|_\Sigma^2)^{(s-1)/2} (1 + \|y-x_\star\|_\Sigma^2)^{(s-1)/2} \\
& \qquad \qquad \times \left[ \frac{(s-1) \kappa(x,y) \Sigma^{-1} (x-x_\star)}{1 + \|x-x_\star\|_\Sigma^2} + \nabla_x \kappa(x,y) \right] \\
\nabla_y c(x,y) & = (1 + \|x-x_\star\|_\Sigma^2)^{(s-1)/2} (1 + \|y-x_\star\|_\Sigma^2)^{(s-1)/2} \\
& \qquad \qquad \times \left[ \frac{(s-1) \kappa(x,y) \Sigma^{-1} (y-x_\star)}{1 + \|y-x_\star\|_\Sigma^2} + \nabla_y \kappa(x,y) \right] \\
\nabla_x \cdot \nabla_y
c(x,y) & = (1 + \|x-x_\star\|_\Sigma^2)^{(s-1)/2} (1 + \|y-x_\star\|_\Sigma^2)^{(s-1)/2} \\
& \quad \times \left[ \frac{(s-1)^2 \kappa(x,y) (x-x_\star)^\top \Sigma^{-2} (y-x_\star)}{(1 + \|x-x_\star\|_\Sigma^2) (1 + \|y-x_\star\|_\Sigma^2)} + \frac{(s-1) (y - x_\star)^\top \Sigma^{-1} \nabla_x \kappa(x,y)}{(1 + \|y-x_\star\|_\Sigma^2)} \right. \\
& \qquad \qquad \qquad \left.  + \frac{(s-1) (x-x_\star)^\top \Sigma^{-1} \nabla_y \kappa(x,y)}{(1 + \|x-x_\star\|_\Sigma^2)} + \nabla_x \cdot \nabla_y \kappa(x,y) \right] .
\end{align*}

These expressions involve gradients of $\kappa(x,y)$, and explicit formulae for these are presented for the choice of $\kappa(x,y)$ corresponding to the Langevin–Stein kernel in <Ref>, and to the KGM–Stein kernel in <Ref>.

To implement Stein $\Pi$-Thinning we require access to both $k_P(x)$ and $\nabla k_P(x)$, the latter for use in the proposal distribution and acceptance probability in mala.
These quantities will now be calculated.
In what follows we assume that $\kappa(x,y)$ is continuously differentiable, so that partial derivatives with respect to $x$ and $y$ can be interchanged.

\begin{align*}
c_0(x) & := c(x,x) \\
& = (1 + \|x-x_\star\|_\Sigma^2)^{s-1} \kappa(x,x) \\
c_1(x) & := \left. \nabla_x c(x,y) \right|_{y \rightarrow x} \\
& = (1 + \|x-x_\star\|_\Sigma^2)^{s-1} \left[ \frac{(s-1) \kappa(x,x) \Sigma^{-1} (x-x_\star)}{(1 + \|x-x_\star\|_\Sigma^2)} + \left. \nabla_x \kappa(x,y) \right|_{y \rightarrow x} \right] \\
c_2(x) & := \left. \nabla_x \cdot \nabla_y c(x,y) \right|_{y \rightarrow x} \\& = (1 + \|x-x_\star\|_\Sigma^2)^{s-1} \left[ \frac{(s-1)^2 \kappa(x,x) (x-x_\star)^\top \Sigma^{-2} (x-x_\star)}{(1 + \|x-x_\star\|_\Sigma^2)^2}  \right. \\
& \hspace{100pt} \left. + \frac{2 (s-1) (x-x_\star)^\top \Sigma^{-1} \left. \nabla_x \kappa(x,y) \right|_{y \rightarrow x}}{(1 + \|x-x_\star\|_\Sigma^2)} + \left. \nabla_x \cdot \nabla_y \kappa(x,y) \right|_{y \rightarrow x}  \right]
\end{align*}

so that

\begin{align}
k_P(x) & := k_P(x,x) = c_2(x) + 2 c_1(x) \cdot \nabla_x \log p(x) + c_0(x) \|\nabla_x \log p(x) \|^2 . \label{eq: kp lin in c}
\end{align}

Let $[\nabla_x c_1(x)]_{i,j} = \partial_{x_i} [c_1(x)]_j$ and $[\nabla_x^2 \log p(x)]_{i,j} = \partial_{x_i} \partial_{x_j} \log p(x)$.
Now we can differentiate (<ref>) to get

\begin{align}
\nabla_x k_P(x) & = \nabla_x c_2(x) + 2 [\nabla_x c_1(x)] [\nabla_x \log p(x)] + 2 [\nabla_x^2 \log p(x)] c_1(x) \nonumber \\
& \hspace{30pt} + [\nabla_x c_0(x)] \|\nabla_x \log p(x)\|^2 + 2 c_0(x) [\nabla_x^2 \log p(x)] [\nabla_x \log p(x)] .  \label{eq: dkp lin in c}
\end{align}

In what follows we also derive explicit formulae for $c_0(x)$, $c_1(x)$ and $c_2(x)$, and hence for $\nabla_x c_0(x)$, $\nabla_x c_1(x)$ and $\nabla_x c_2(x)$, for the case of the Langevin–Stein kernel in <Ref>, and the KGM–Stein kernel in <Ref>.

§.§ Explicit Formulae for the Langevin–Stein Kernel

The Langevin–Stein kernel from <Ref> corresponds to the choice $s=1$ and $\kappa(x,y)$ the inverse multi-quadric kernel, so that

\begin{align*}
\kappa(x,y) & = (1 + \|x-y\|_\Sigma^2)^{-\beta} \\
\nabla_x \kappa(x,y) & = - 2 \beta (1 + \|x-y\|_\Sigma^2)^{-\beta-1} \Sigma^{-1}(x-y) \\
\nabla_y \kappa(x,y) & = 2 \beta (1 + \|x-y\|_\Sigma^2)^{-\beta-1} \Sigma^{-1}(x-y) \\
\nabla_x \cdot \nabla_y \kappa(x,y) & = - 4 \beta (\beta + 1) (1 + \|x-y\|_\Sigma^2)^{-\beta-2} (x-y)^\top \Sigma^{-2} (x-y) \\
& \hspace{100pt} + 2 \beta \text{tr}(\Sigma^{-1}) (1 + \|x-y\|_\Sigma^2)^{-\beta-1}.
\end{align*}

Evaluating on the diagonal:

\begin{align*}
\kappa(x,x) & = 1 \\
\left. \nabla_x \kappa(x,y) \right|_{y \rightarrow x} = \left. \nabla_y \kappa(x,y) \right|_{y \rightarrow x} & = 0 \\
\left. \nabla_x \cdot \nabla_y \kappa(x,y) \right|_{y \rightarrow x} & = 2 \beta \text{tr}(\Sigma^{-1}) ,
\end{align*}

so that $c_0(x) = 1$, $c_1(x) = 0$, $c_2(x) = 2 \beta \text{tr}(\Sigma^{-1})$.
Differentiating these formulae, $\nabla_x c_0(x) = 0$, $\nabla_x c_1(x) = 0$, $\nabla_x c_2(x) = 0$.

§.§ Explicit Formulae for the KGM–Stein Kernel

The KGM kernel of order $s$ from <Ref> corresponds to the choice

\begin{align*}
\kappa(x,y) = (1 + \|x-y\|_\Sigma^2)^{-\beta} + \frac{1 + (x-x_\star)^\top \Sigma^{-1} (y-x_\star)}{(1 + \|x-x_\star\|_\Sigma^2)^{s/2} (1 + \|y-x_\star\|_\Sigma^2)^{s/2} } ,
\end{align*}

for which we have

\begin{align*}
\nabla_x \kappa(x,y) & = - 2 \beta (1 + \|x-y\|_\Sigma^2)^{-\beta-1} \Sigma^{-1} (x-y)  \\
& \qquad + \frac{ \Sigma^{-1}(y-x_\star) - s[1 + (x-x_\star)^\top \Sigma^{-1} (y-x_\star)] \Sigma^{-1}(x-x_\star) (1 + \|x-x_\star\|_\Sigma^2)^{-1} }{ (1 + \|x-x_\star\|_\Sigma^2)^{s/2} (1 + \|y-x_\star\|_\Sigma^2)^{s/2}  } \\
\nabla_y \kappa(x,y) & = 2 \beta (1 + \|x-y\|_\Sigma^2)^{-\beta-1} \Sigma^{-1} (x-y)  \\
& \qquad + \frac{ \Sigma^{-1}(x-x_\star) - s[1 + (x-x_\star)^\top \Sigma^{-1} (y-x_\star)] \Sigma^{-1}(y-x_\star) (1 + \|y-x_\star\|_\Sigma^2)^{-1} }{ (1 + \|x-x_\star\|_\Sigma^2)^{s/2} (1 + \|y-x_\star\|_\Sigma^2)^{s/2}  } \\
\nabla_x \cdot \nabla_y \kappa(x,y) & = - 4 \beta (\beta + 1) (1 + \|x-y\|_\Sigma^2)^{-\beta-2} (x-y)^\top \Sigma^{-2} (x-y) + 2 \beta \text{tr}(\Sigma^{-1}) (1 + \|x-y\|_\Sigma^2)^{-\beta-1} \\
& \qquad + \frac{ \left[ \begin{array}{l} \text{tr}(\Sigma^{-1}) - s(1 + \|x-x_\star\|_\Sigma^2)^{-1} (x-x_\star)^\top \Sigma^{-2} (x-x_\star) \\
\qquad - s(1 + \|y-x_\star\|_\Sigma^2)^{-1} (y-x_\star)^\top \Sigma^{-2} (y-x_\star) \\
\qquad + s^2[1 + (x-x_\star)^\top \Sigma^{-1} (y-x_\star)] (1 + \|x-x_\star\|_\Sigma^2)^{-1} (1 + \|y-x_\star\|_\Sigma^2)^{-1} \\
\hspace{150pt} \times (x-x_\star)^\top \Sigma^{-2} (y-x_\star) \end{array} \right] }{ (1 + \|x-x_\star\|_\Sigma^2)^{s/2} (1 + \|y-x_\star\|_\Sigma^2)^{s/2} } .
\end{align*}

Evaluating on the diagonal:

\begin{align*}
\kappa(x,x) & = 1 + (1 + \|x-x_\star\|_\Sigma^2)^{-s+1} \\
\left. \nabla_x \kappa(x,y) \right|_{y \rightarrow x} = \left. \nabla_y \kappa(x,y) \right|_{y \rightarrow x} & = - (s-1) \Sigma^{-1} (x-x_\star) (1 + \|x-x_\star\|_\Sigma^2)^{-s} \\
\left. \nabla_x \cdot \nabla_y \kappa(x,y) \right|_{y \rightarrow x} & = 2 \beta \text{tr}(\Sigma^{-1}) + \text{tr}(\Sigma^{-1}) (1 + \|x-x_\star\|_\Sigma^2)^{-s} \\
& \qquad + s(s-2) (1 + \|x-x_\star\|_\Sigma^2)^{-s-1} (x-x_\star)^\top \Sigma^{-2} (x-x_\star)
\end{align*}

so that

\begin{align*}
c_0(x) & = 1 + (1 + \|x-x_\star\|_\Sigma^2)^{s-1} \\
c_1(x) & = (s-1) (1 + \|x-x_\star\|_\Sigma^2)^{s-2} \Sigma^{-1} (x-x_\star) \\
c_2(x) & = \frac{ [(s-1)^2 (1 + \|x-x_\star\|_\Sigma^2)^{s-1} - 1] (x-x_\star)^\top \Sigma^{-2} (x-x_\star) }{ (1 + \|x-x_\star\|_\Sigma^2)^2 } + \frac{\text{tr}(\Sigma^{-1})[ 1 + 2\beta (1 + \|x-x_\star\|_\Sigma^2)^s ]}{ (1 + \|x-x_\star\|_\Sigma^2) } .
\end{align*}

Differentiating these formulae:

\begin{align*}
\nabla_x c_0(x) & = 2(s-1) (1 + \|x-x_\star\|_\Sigma^2)^{s-2} \Sigma^{-1}(x-x_\star) \\
\nabla_x c_1(x) & = 2(s-1)(s-2)(1 + \|x-x_\star\|_\Sigma^2)^{s-3} [\Sigma^{-1} (x-x_\star)] [\Sigma^{-1}(x-x_\star)]^\top \\
& \qquad + (s-1) (1 + \|x-x_\star\|_\Sigma^2)^{s-2} \Sigma^{-1} \\
\nabla_x c_2(x) & = 2 (s-1)^2 (s-3) (1 + \|x-x_\star\|_\Sigma^2)^{s-4} [(x-x_\star)^\top \Sigma^{-2} (x-x_\star)] \Sigma^{-1} (x-x_\star) \\
& \qquad + 2(s-1)^2 (1 + \|x-x_\star\|_\Sigma^2)^{s-3} \Sigma^{-2}(x-x_\star) \\
& \qquad + 4 \beta \text{tr}(\Sigma^{-1}) (s-1) (1 + \|x-x_\star\|_\Sigma^2)^{s-2} \Sigma^{-1} (x-x_\star) \\
& \qquad - 2 (1 + \|x-x_\star\|_\Sigma^2)^{-2}  [ \Sigma^{-2} (x-x_\star) + \text{tr}(\Sigma^{-1}) \Sigma^{-1} (x-x_\star) ] \\
& \qquad + 4 (1 + \|x-x_\star\|_\Sigma^2)^{-3} [(x-x_\star)^\top \Sigma^{-2} (x-x_\star)] \Sigma^{-1} (x-x_\star) .
\end{align*}

These complete the analytic calculations necessary to compute the Stein kernel $k_P$ and its gradient.

§ EMPIRICAL ASSESSMENT

This appendix contains full details of the empirical protocols that were employed and the additional empirical results described in the main text.
<Ref> discusses the effect of dimension on our proposed $\Pi$.
Additional illuatrative results from <Ref> are contained in <Ref>.
The full details for how mala was implemented are contained in <Ref>.
An additional illustration using a garch model is presented in <Ref>.
The full results for S$\Pi$IS-MALA are contained in <Ref>, and in <Ref> the convergence of the sparse approximation provided by S$\Pi$T-MALA to the optimal weighted approximation is investigated.
Finally, the performance of ksd is quantified using the 1-Wasserstein divergence in <Ref>.

§.§ The Effect of Dimension on $\Pi$

The improvement of Stein $\Pi$-Importance Sampling over the default Stein importance sampling algorithm (i.e. $\Pi = P$) can be expected to reduce as the dimension $d$ of the target $P$ is increased.
To see this, consider the Langevin–Stein kernel

\begin{align}
k_P(x) = c_1 + c_2 \| \nabla \log p(x) \|_\Sigma^2 \label{eq: kp lang}
\end{align}

for some $c_1, c_2 > 0$; see <Ref>.
Taking $P = \mathcal{N}(0,I_{d \times d})$, for which the length scale matrix $\Sigma$ appearing in <Ref> is $\Sigma = I_{d \times d}$, we obtain

\begin{align*}
k_P(x) = c_1 + c_2 \|x\|^2  .
\end{align*}

However, the sampling distribution $\Pi$ defined in (<ref>) depends on $k_P$ only up to an unspecified normalisation constant; we may therefore equally consider the asymptotic behaviour of $\tilde{k}_P(x) := k_P(x) / d$.
Let $X \sim P$.
Then $\mathbb{E}[\tilde{k}_P(X)] = c_2$ is a $d$-independent constant, and

\begin{align*}
\left\| \tilde{k}_P - \mathbb{E}[\tilde{k}_P(X)] \right\|_{L^2(P)}^2 = \int \left[ \frac{k_P(x) - (c_1 + c_2 d)}{d}  \right]^2 \; \mathrm{d}P(x) = \frac{2c_2^2}{d} \rightarrow 0
\end{align*}

as $d \rightarrow \infty$.
This shows that $\tilde{k}_P$ converges to a constant function in $L^2(P)$, and thus for “typical” values of $x$ in the effective support of $P$,

\begin{align*}
\pi(x) \propto p(x) \sqrt{\tilde{k}_P(x)} \stackrel{\approx}{\propto} p(x) ,
\end{align*}

so that $\Pi \approx P$ in the $d \rightarrow \infty$ limit.
This intuition is borne out in simulations involving both the Langevin–Stein kernel (as just discussed) and also the KGM3–Stein kernel.
Indeed, <Ref> shows that as the dimension $d$ is increased, the marginal distributions of $\Pi$ become increasingly similar to those of $P$.

The effect of dimension on $\Pi$:
Here $P$ was taken to be the standard Gaussian distribution $\mathcal{N}(0,I_{d \times d})$ in $\mathbb{R}^d$ and the proposed distribution $\Pi$ was computed.
The marginal distribution of the first component of $\Pi$ is plotted for $d \in \{1,2,10\}$, for both (a) the Langevin–Stein kernel and (b) the KGM3–Stein kernel.

§.§ 2D Illustration from the Main Text

<Ref> of the main text contained a 2-dimensional illustration of Stein $\Pi$-Importance Sampling and presented the distributions $\Pi$ corresponding to different choices of Stein kernel.
Here, in <Ref>, we present the mean ksd for Stein $\Pi$-Importance Sampling performed using the Langevin–Stein kernel (purple), the KGM3–Stein kernel (blue), and the Riemann–Stein kernel (red), corresponding to the sampling distributions $\Pi$ displayed in <Ref> of the main text.

For this experiment, exact sampling from both $P$ and $\Pi$ was performed using a fine grid on which all probabilities were calculated and appropriately normalised.
Results are in broad agreement with the 1-dimensional illustration contained in the main text, in the sense that in all cases Stein $\Pi$-Importance Sampling provides a significant improvement over the default Stein importance sampling method with $\Pi$ equal to $P$.

Assessing the performance of the sampling distributions $\Pi$ shown in <Ref>.
The mean kernel Stein discrepancy (KSD) for computation performed using the Langevin–Stein kernel (purple), the KGM3–Stein kernel (blue), and the Riemann–Stein kernel (red); in each case, KSD was computed using the same Stein kernel used to construct $\Pi$.  Solid lines indicate the baseline case of sampling from $P$, while dashed lines indicate the proposed approach of sampling from $\Pi$.
(The experiment was repeated 10 times and standard error bars are plotted.)

§.§ Implementation of mala

For implementation of mala in <Ref> we are required to specify a step size $\epsilon$ and a preconditioner matrix $M$.
In general, suitable values for both of these parameters will be problem-dependent.
Standard practice is to perform some form of manual or automated tuning to arrive at parameter values for which the average acceptance rate is close to 0.57, motivated by the asymptotic analysis of 47.
Adaptive mcmc algorithms, which seek to optimise the parameters of mcmc algorithms such as mala during the warm-up period, provide an appealing solution, and was the approach taken in this work.

The adaptive mala algorithm which we used is contained in <Ref>, where we have let $\texttt{MALA}(x, \epsilon, M, n, k_P)$ denote the output from the preconditioned mala with initial state $x$, step size $\epsilon$, preconditioner matrix $M$, and chain length $n$, described in <Ref>.
In <Ref>, we use $\mathrm{cov}(\cdot)$ to denote the sample covariance matrix.
The algorithm monitors the average acceptance rate and increases or decreases it according to whether it is below or above, respectively, the 0.57 target.
For the preconditioner matrix, the sample covariance matrix of samples obtained from the penultimate tuning run of mala is used.
For all experiments that we report using mala, we set $\epsilon_0=1$, $M_0=I_d$, $h=10$, and $\alpha_1 = \cdots = \alpha_9 = 0.3$.
The warm-up epoch lengths were $n_0 = \dots = n_8 = 1,000$ and the final epoch length was $n_9 = 10^5$.
The samples $\smash{ \{x_{h-1,1} , \dots , x_{h-1,n_{i-1}}\} }$ from the final epoch are returned, and constituted output from mala for our experimental assessment.

[1]
$x_{0,0}$ (initial state), $\epsilon_0$ (initial step size), $M_0$ (initial preconditioner matrix), $\{n_i\}_{i=0}^{h-1}$ (epoch lengths), $\{\alpha_i\}_{i=1}^{h-1}$ (learning schedule), $h$ (number of epochs), $k_P$ (Stein kernel)

$\{x_{0,1} \ldots, x_{0,n_0}\} \leftarrow \texttt{MALA}(x_{0,0}, \epsilon_0, M_0, n_0, k_P)$
$x_{i,0} \leftarrow x_{i-1, n_{i-1}}$
$\rho_{i-1} \leftarrow \frac{1}{n_{i-1}} \sum_{j=1}^{n_{i-1}} 1_{x_{i-1,j} \neq x_{i-1,j-1}}$ Average acceptance rate for chain $i$
$\epsilon_{i} \leftarrow \epsilon_{i-1} \exp(\rho_{i-1} - 0.57)$ Update step size
$M_i \leftarrow \alpha_i M_i + (1 - \alpha_i) \mathrm{cov}(\{x_{i-1,1} \ldots, x_{i-1,n_{i-1}}\})$ Update preconditioner matrix
$\{x_{i,1} \ldots, x_{i,n_i}\} \leftarrow \texttt{MALA}(x_{i,0}, \epsilon_i, M_i, n_i, k_P)$

Adaptive MALA

To sample from $P$ instead of $\Pi$, we used <Ref> we formally set $k_P(x) = 1$ for all $x \in \mathbb{R}^d$, which recovers $\Pi = P$ as the target.

§.§ Illustration on a GARCH Model

This appendix contains an additional illustrative experiment, concerning a garch model that is a particular instance of a model from the  database discussed in <Ref>.
The purpose of this illustration is to facilitate an empirical investigation in a slightly higher dimension ($d=4$) and to explore the effect of changing the order $s$ of the KGM–Stein kernel defined in <Ref>.

First we describe the garch model that was used.
These models are widely-used in econometrics to describe time series data $\{y_t\}_{t=1}^{n}$ in settings where the volatility process is assumed to be time-varying (but stationary).
In particular, we consider the GARCH(1,1) model

\[
\begin{aligned}
y_t &= \phi_1 + a_t, \\
a_t &= \sigma_t\epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0,1), \\
\sigma^2_{t} &= \phi_2 + \phi_3 a^{2}_{t-1} + \phi_{4} \sigma^2_{t-1},
\end{aligned}
\]

where $\phi_{2} > 0$, $\phi_{3} > 0$, $\phi_{4} > 0$, and $\phi_{3} + \phi_{4} < 1$ are the model parameters, constrained to a subset of $\mathbb{R}^4$.
For ease of sampling, a change of variables $\tau: (\phi_1, \phi_2, \phi_3, \phi_4) \mapsto \theta$ is performed in such a way that the parameter $\theta \in \mathbb{R}^4$ is unconstrained.
Assuming an improper flat prior on $\theta$, the log-posterior density for $\theta$ is given up to an additive constant by

\[
\log p(\theta \:|\: y_1,\ldots,y_n) \stackrel{+C}{=} \sum_{t=1}^{n} \left[-\frac{1}{2}\log(\sigma^{2}_{t}) - \frac{y^{2}_{t}}{2\sigma^{2}_{t}}\right]  + \log |J_{\tau^{-1}}(\theta)|,
\]

where $|J_{\tau^{-1}}(\theta)|$ is the Jacobian determinant of $\tau^{-1}$.

For this illustration, real data were provided within the model description of , for which the estimated maximum a posteriori parameter is $\smash{\hat{\phi} = (5.04, 1.36, 0.53, 0.31)}$.
The marginal distributions of $\Pi$ corresponding to the KGM–Stein kernels of orders $s \in \{2,3,4\}$ are compared to the marginals of $P$ in <Ref>.
It can be seen that higher orders $s$ correspond to greater over-dispersion of $\Pi$; this makes intuitive sense since larger $s$ corresponds to a more stringent ksd (controlling the convergence of moments up to order $s$) which places greater emphasis on how the tails of $P$ are approximated.
Further, for the final skewed marginal of $P$, we note that the distribution $\Pi$ exaggerates the skew, placing more of its mass in the tail of the direction which is positively skewed.
Further discussion of skewed targets is contained in <Ref>.

Illustrating the shape of $\Pi$ based on the KGM$s$–Stein kernel for a GARCH(1,1) model, controlling convergence of moments up to order $s \in \{2, 3, 4\}$.
The marginal density functions of each distribution were approximated using one-million samples obtained using mcmc.

§.§ Stein $\Pi$-Importance Sampling for

To introduce objectivity into our assessment, we exploited the  benchmark  [43].
This ongoing project is an attempt toward standardised benchmarking, consisting of a collection of posteriors to be numerically approximated.
The test problems in  are defined in the  probabilistic programming language, and so  [49] was used to directly access posterior densities and their gradients as required.
The ambition of  is to provide an extensive set of benchmark tasks; at the time we conducted our research,  was at Version 0.4.0 and contained 149 models, of which 47 came equipped with a gold-standard sample of size $n=10^3$, generated from a long run of Hamiltonian Monte Carlo (the No-U-Turn sampler in ).
Of these 47 models, a subset of 40 were found to be compatible with , which was at Version 1.0.2 at the time this research was performed.
The version of  that we used was  Version 2.31.0 (Unix).
Thus we used a total of 40 test problems for our empirical assessment.

For each test problem, a total of 10 replicate experiments were performed and standard errors were computed.
A sampling method was defined as being significantly better for approximation of a given target, compared to all other methods considered, if had lower mean ksd and the standard error bar did not overlap with the standard error bar of any other method.
<Ref> in the main text summarises the performance of S$\Pi$IS-MALA, fixing the number of samples to be $n = 3 \times 10^3$.
In this appendix, full empirical results are provided.

For sampling from mala, we used the adaptive algorithm described in <Ref> with a final epoch of length $n_{\max} = 10^5$.
Then, whenever a set of $n \ll n_{\max}$ consecutive samples from mala are required for our experimental assessment, these were obtained by selecting at random a consecutive sequence of length $n$ from the total chain of length $10^5$.
This ensures that the performance of unprocessed mala that we report is not negatively affected by burn-in, in so far as is practical to control.

Full results are presented in Figure <ref>.
These results broadly support the interpretation that S$\Pi$IS-MALA usually outperforms SIS-MALA, or otherwise both methods provide a similar level of performance, for the sufficiently large sample sizes $n$ considered.
The sample size threshold at which S$\Pi$IS-MALA outperforms SIS-MALA appears to be dimension-dependent.
A notable exception is panel 29 of Figure <ref>, a $d=10$ dimensional task for which S$\Pi$IS-MALA provided a substantially worse approximation in ksd for the range of values of $n$ considered.

Benchmarking on .
Here we compared raw output from MALA (dotted lines) with the post-processed output provided by the default Stein importance sampling method of 40 (SIS-MALA; solid lines) and the proposed Stein $\Pi$-Importance Sampling method (S$\Pi$IS-MALA; dashed lines).
The Langevin (purple) and KGM3–Stein kernels (blue) were used for SIS-MALA and S$\Pi$IS-MALA and the associated ksd are reported as the number $n$ of iterations of MALA is varied.
Ten replicates were computed and standard errors were plotted.
The name of each model is shown in the title of the corresponding panel, and the dimension $d$ of the parameter vector is given in parentheses.
[Langevin–Stein kernel:  MALA, SIS-MALA, S$\Pi$IS-MALA.
KGM3–Stein kernel:  MALA, SIS-MALA, S$\Pi$IS-MALA.]

§.§ Stein $\Pi$-Thinning for

The results presented in the main text concerned $n = 3 \times 10^3$ samples from mala, which is near the limit at which the optimal weights $w^\star$ can be computed in a few seconds on a laptop PC.
For larger values of $n$, sparse approximation methods are likely to required.
In the main text we presented Stein $\Pi$-Thinning, which employs a greedy optimisation perspective to obtain a sparse approximation to the optimal weights at cost $O(m^2 n)$, where $m$ are the number of greedy iterations performed.
Explicit and verifiable conditions for the strong consistency of the resulting S$\Pi$T-MALA algorithm were established in <Ref>.
The purpose of this appendix is to empirically explore the convergence of S$\Pi$T-MALA using the  test bed.

In the experiments we report the number of mala samples was fixed to $n = 10^3$ and the number of greedy iterations was varied from $m = 1$ to $m = 10^3$.
The results, in <Ref>, indicate that for most models in  the minimum value of ksd is approximately reached when $m$ is anywhere from $\frac{n}{10}$ to $\frac{n}{2}$, representing a modest but practically significant reduction in computational cost compared to S$\Pi$IS-MALA.
This agrees with the qualitative findings reported in the original Stein thinning paper of 46.

Benchmarking on .
Here we investigate the convergence of the sparse approximation provided by the proposed Stein $\Pi$-Thinning method (S$\Pi$T-MALA).
The Langevin (purple) and KGM3–Stein kernels (blue) were used for S$\Pi$T-MALA and the associated ksd are reported as the number $m$ of iterations of Stein thinning is varied.
Ten replicates were computed and standard errors were plotted.
The name of each model is shown in the title of the corresponding panel, and the dimension $d$ of the parameter vector is given in parentheses.
[Langevin–Stein kernel:S$\Pi$T-MALA.
KGM3–Stein kernel:S$\Pi$T-MALA.]

§.§ Performance of Stein Discrepancies

The properties of Stein discrepancies was out of scope for this work.
Nonetheless, there is much interest in better understanding the properties of KSDs, and in this appendix the performance of S$\Pi$IS-MALA in terms of 1-Wasserstein divergence is reported.
This was made possible since  supplies a set of posterior samples obtained from a long run of Hamiltonian Monte Carlo (the No-U-Turn sampler in ) which we treat as a gold standard.

Full results are presented in <ref>.
Broadly speaking, for most models the minimisation of ksd seems to be associated with minimisation of 1-Wasserstein distance, however there are some models for which minimisation of ksd is loosely, if at all, related to minimisation of 1-Wasserstein divergence.
In these cases, we attribute this performance to the blindness to mixing proportions phenomena, described in 53, 38, 42.
Convergence in 1-Wasserstein is equivalent to weak convergence plus convergence of the first moment, so the KGM–Stein kernels of order $s \geq 1$ control convergence in 1-Wasserstein.
In <Ref> we proved that S$\Pi$IS-MALA is strongly consistent in ksd for the KGM–Stein kernel in the case $s=1$, so we can expect strong consistency in 1-Wasserstein divergence for S$\Pi$IS-MALA in this case as well.
It is interesting to observe that better 1-Wasserstein quantisations tend to be provided by S$\Pi$IS-MALA compared to SIS-MALA when either the Langevin–Stein or KGM–Stein kernel are used.

The development of improved Stein discrepancies is an active area of research, and we emphasise that the methodology developed in this work can be applied to any ksd, including potentially ksd with better or more direct control over standard notions of convergence (such as 1-Wasserstein) that in the future may be developed.

Performance of Stein discrepancies on .
Here we compared raw output from MALA (dotted lines) with the post-processed output provided by the default Stein importance sampling method of 40 (SIS-MALA; solid lines) and the proposed Stein $\Pi$-Importance Sampling method (S$\Pi$IS-MALA; dashed lines).
The Langevin (purple) and KGM3–Stein kernels (blue) were used for SIS-MALA and S$\Pi$IS-MALA, and the 1-Wasserstein divergence is reported as the number $n$ of iterations of MALA is varied.
Ten replicates were computed and standard errors were plotted.
The name of each model is shown in the title of the corresponding panel, and the dimension $d$ of the parameter vector is given in parentheses.
[Legend:  Raw MALA.
Langevin–Stein kernel:  SIS-MALA, S$\Pi$IS-MALA.
KGM3–Stein kernel:  SIS-MALA, S$\Pi$IS-MALA.]

§.§ Investigation for a Skewed Target

This final appendix contrasts the 1-Wasserstein optimal sampling distribution $\Pi_1$ (c.f. <Ref>), with the choice of $\Pi$ that we recommended in (<ref>).
In particular, we focus on the KGM3–Stein kernel under a heavily skewed $P$, for which $\Pi_1$ and $\Pi$ can be markedly different.

For this investigation a bivariate skew-normal target was constructed, where the density is given by $p(x_1, x_2) = 4 \phi(x_1)\Phi(6 x_1)\phi(x_2)\Phi(-3 x_2)$, with $\phi$ and $\Phi$ respectively denoting the density and distribution functions of a standard Gaussian.
The density $p$ of $P$, together with the marginal densities of $\Pi_1$ and $\Pi$, are plotted in <Ref>.
It can be seen that, while both $\Pi_1$ and $\Pi$ are over-dispersed with respect to $P$, our recommended $\Pi$ assigns proportionally more mass to the tail that is positively skewed.

The performance of Stein $\Pi$-Importance Sampling based on $\Pi_1$ and $\Pi$ is compared in <Ref>.
Though both choices lead to an improvement relative to Stein importance sampling algorithm with $\Pi = P$, the use of $\Pi$ leads to a significant further reduction (on average) in ksd compared to $\Pi_1$.
Based on our investigations, this finding seems general; the use of $\Pi_1$ does not realise the full potential of Stein $\Pi$-Imporance sampling when the target is skewed.

Comparing the proposed distribution $\Pi$ (KGM3; based on the KGM3–Stein kernel) to $\Pi_1$ (1Wass.; the optimal choice for 1-Wasserstein quantisation from <Ref>) for a bivariate skew-normal target ($d=2$).
The marginal density functions of each distribution were approximated using $10^6$ samples from mcmc.

Comparing the performance of using the proposed distribution $\Pi$ (KGM3; based on the KGM3–Stein kernel) to $\Pi_1$ (1Wass.; the optimal choice for 1-Wasserstein quantisation from <Ref>) for a bivariate skew-normal target ($d=2$).
The mean kernel Stein discrepancy (KSD) for Stein $\Pi$-Importance Sampling was estimated; in each case, the KSD based on the KGM3–Stein kernel was computed.
Solid lines indicate the baseline case of sampling from $P$, while dashed lines indicate sampling from $\Pi$.
(The experiment was repeated 10 times and standard error bars are plotted.)