Stein $\Pi$-Importance Sampling
Authors: Congye Wang
Wilson Ye Chen
Heishiro Kanagawa
Chris J. Oates
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Bayesian
discrepancy
kernel
sampling
Stein's method

Abstract: 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 PosteriorDB benchmark, a significant improvement over the analogous post-processing of $P$-invariant Markov chains is reported.

Stein I1-Importance Sampling
Congye Wang', Wilson Ye Chen², Heishiro Kanagawa', Chris. J. Oates' Newcastle University, UK
University of Sydney, Australia
Abstract
Stein discrepancies have emerged as 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 II-invariant Markov chain to obtain a consistent approximation P. the intended target. Surprisingly, the optimal choice of not identical to the target P; we therefore propose explicit construction for II based on novel variational argument. Explicit conditions for convergence of Stein [1-mmpotance Sampling are established. For 70% of tasks in the Posteri orDB benchmark, significant improvement over the analogous post- processocessing of P-invariant Markov chains is reported.
Introduction
Stein discrepancies are a class of statistical divergences that can be computed without access to normalisation constant. Originally conceived tool measure the performance of sampling methods (Gorham and Mackey, 2015), these discrepancies have since found wide-ranging statistical applications (see the review of Anastasiou al., 2023). Our focus here is the use of Stein discrepancies for retrospective improvement of Markov chain Monte Carlo (MCMC), and here two main techniques have been proposed: (i) Stein importance sampling (Liu and Lee, 2017: Hodgkinson al., 2020), and (ii) Stein thinning (Riabiz et al., 2022). In Stein importance sampling (also called black box importance sampling), the samples assigned weights such that Stein discrepancy between the weighted empirical measure and the target P is minimised. Stein thinning constructs sparse approximation this optimally weighted measure at lower computational and storage cost. Together, these techniques provide a powerful set of post-processing tools for MCMC, with subsequent authors proposing range of generalisations and extensions (Teymur et al., 2021: Chopin and Ducrocq, 2021; Hawkins et al., 2022; Fisher and Oates, 2023; Bénard et al., 2023).
The consistency of these algorithms has been established in the setting of approximate, --invariant MCMC, motivated by challenging inference problems where only approximate sampling can be performed. In these settings, IT is implicitly approximation P that is as accurate possible subject to computational budget. However, the critical question of how to design Markov chains that are well-suited to such the form of a specific construction for II derived from a novel variational argument. Surprisingly. we are able demonstrate a substantial improvement using the proposed 11, compared to the case where II and P are equal. The paper proceeds as follows: Section 2 presents an abstract formulation empirical measures are reviewed. Section
has yet be addressed. This paper provides a solution, in
of the task and existing results for aptimally-weighee derives our proposed choice of II and establishes that Stein post-processing of samples from II-invariant Metropolis--djustee Langevin algorithm (MALA) provides consistent approximation of P. The approach is stress-tested using the recently released PosteriorDy orDB suite of benchmark tasks in Section 4, before concluding with a discussion in Section
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

2 Background
To properly contextualise our discussion we start with abstract mathematical description of the task. Let P be a probability measure on a measurable space X. Let ((2') be the set of all probability measures on X. Let Dr P(x)) + O0, so] be statistical divergence for measuring the quality of an approximation Q to P, meaning that (Q()) = 0 if and only if Q = P. In this work we consider approximations whose support is contained
finite set ( Inf X, and in particular we
consider optimal approximations of the form
PR =
Σ o (zs))
w E
arg min Dp w20, mu
Σ (
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 (I1,...,fn) minimise the approximation error Dp(?*?? Before specialising to Stein discrepancies, it is helpful review existing results for some standard statistical divergences Dp.
2.1 Wasserstein Divergence
Optimal quantisation focuses on the r- Wasserstetin (r 2 1) family of statistical divergences Dp(Q)= imf.)r((..)) || yl''dr(z y). where I(P, Q) denotes the set of couplings' of P,Q E P(R and the divergence is finite whenever P and Q have finite r-th moment. Assuming the states are distinet, the corresponding optimal weights are w = P(A.) where A, is the Voronoi neighbourhood of x; in Rdd Optimal states achieve the minimal quantisation error for P:
In {
Cn (P)
inf ERR
Dp
w ((x;))
the smallest value of the divergence among optimally- weighted distributions supported on most states. Though the dependence of optimal states on n and P can be complicated, we can broaden our perspective consider asymptotically optimal states, whose asymptotic properties can be precisely characterised. To this end, for A C IRd, let ((A)) denote the uniform distribution on A, and define the universal constant (10, 11 infa21 n//d en. (2,002((10 1]()). Suppose that admits density Rd, Then the rth quantisation coefficient of P on IRd, defined as
(d+r)/d
(f (frlo)taato dæ
Cx(P) ((0)
plays a central role in the classical theory of quantisation, being the rate constant in the asymptotic convergence of the minimal quantisation error; lim,- n² nt"""n,r-P) C, C,(P); see Theorem 6.2 of Graf and Luschgy (2007). This suggests a natural definition; collection fa1,.... ;n) ) is called asymptotically optimal if
lim Dr n00
Σ ot ((x))) = C.()),
which amounts to P asymptotically attaining the minimal quantisation error (P). The main result here that, if {ai...., Fn are asymptotically optimal, then E 5(x,) II, where convergence is in distribution and II, is the distribution whose density is (x) x d(x)d/(d+r) see Theorem 7.5 of Graf and Luschgy (2007). This provides us with 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 Wassersteii divergence is. unfortunately, not well-suited for use in the motivating Bayesian context. In particular, computing the optimal weights P(A,) requires knowledge of P. which typically available when P is implicitly defined via an intractable normalisation constant. On
AA coupling 7 E T(P, Q) is distribution 7 E (RR x TRd) whose marginal distributions are P and Q.
The Voronoi neighbourhood of is the set A, (w Idd ||æ xill minj=1...... ||. 2; Il]

---

the other hand. the optimal sampling distribution II is explicit and can be sampled (for example using MCMC); for discussion of random quantisers in this context see Graf and Luschgy (2007, Chapter 9), Cohort (2004, p126) and Sonnleitner (2022, Section 4.5). The simple form of TT is feature of the classical approach to quantisation that we will attempt to mimic in the sequel.
2.2 Kernel Discrepancies
The theory of quantisation using kernels is less well-developed. A kernel is a measurable, symmetric, desitive-effinte function le : X' X R. From the Moore-Arome AronssaAonnnajn theorem, there is a unique Hilbert space H((k) for which k is a reproducing kernel, meaning that k(-, x) E H((k) for all æ E X and (f, k(,, 2)) R((a)) = f(x) for all H(k)) and all X. Assuming that H(kk) C LP)), we can define the weak (or Pettis) integral
uP(.) = ((,22) dP(x),
called the kernel mean embedding of P in (k).. The kernel discrepancy is then defined as the norm of the difference between kernel mean embeddings
VI k(z,g) d(Q () (y)
Dp(Q)= I pu (Pllx(()))
where, to be consistent with our carlier notation, we adopt the convention that Dp(Q) is infinite whenever HKk) & (Q).. The second equality (2) follows immediately from the stated properties of a reproducing kernel. To satisfy the requirement of a statistical divergence, we assume that the kernel characteristic, meaning that Pr μo 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 unde-understood. Indeed, given distinct states fz....., }. the corresponding optimal weights w" =w",,... wx) are the solution to the constraee quadratic program
arg min w T Kw - 22 To w
S.L.
w 2 0, 1ww=1
wERd
where Ki = k(z,,1,) xT)) and = (2,( This program does not admit closed-fome solution, but can be numerically solved. To the best of our knowledge, the only theoretical analysis of approxi- mations based on (3) is due to Hayakawa et al. (2022), who established rates for the convergence of Pt to 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 (3), relaxations of this program have been well-studied. The simplest relaxation of (3) is to remove both the positivity (w 2 0) and normalisation (1 w = 1) constraints, in which case the optimal weights have the explicit representation w" K-11 2. 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) states are randomly sampled, the optimal sampling distribution will be n-dependent (Bach, 2017) and over-dispersed with respect to the distributional target (Briol et al., 2017), and (ii) space-filling designs are asymptotically optimal for typical stationary kernels on bounded domains X' C TRd (Briol al., 2019). Analysis of optimal states on unbounded domains appears to be more difficult; see e.g. Karvonen et al. (2021). Relaxation of either the positivity or normalisation constraints results in approximations that behave similarly to kernel cubature (see, respectively, Ehler et al., 2019; Karvonen et al., 2018). However, relaxation of either constraint can result in the failure of P* to be an element of P(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 widely used. The application of kernel discrepancies to an implicitly defined distribu- tional target, such as posterior distribution in Bayesian analysis, is made possible by the use of Stein kernel; a P-dependent kernel ie kp for which yr (x) = 0 for all x E X' (Oates et al., 2017). The associated kernel discrepancy
II kp(z. v) dQ(z)dQ())
DP(Q)= [|
(4)
3

---

is called kernel Stein discrepancy (KSD) (Chwialkowski et al., 2016; Liu et al., 2016; Gorham and Mackey, 2017), and this will be a key tool in our methodological development. The corresponding optimally weighted approximation Pt is the Stein importance sampling method of Liu and Lee (2017). To retain clarity of presentation in the main text, we defer all details on the construction of Stein kermels to Appendix A.
2.3 Sparse Approximation
If the number 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 Pr.m 1 > d(y:), where each ye is an element from ( Fn.. The canonical example is the greedy algorithm which, at iteration j, selects a state
Dp (for0) + 1 8(w))
j-1
Mj arg min De (
(5)
for which the statistical divergence is minimised. In the context of kernel discrepancy, the greedy algorithm (5) has computational cost O(m²n), which compares favourably? with the cost of solving (3) when m < Furthermore, under appropriate assumptions, the sparse approximation converges to the optimally weighted approximation; Dp(Pn,m) + Dp(P,) as m - oo with n fixed. See Teymur et al. (2021) 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-Won Wolfe algorithm (also called herding, the conditional gradient method), for which similar asymptotic result can be shown to hold (Chen et al., 2010; Bach et al., 2012; Chen et al., 2018). Related work includes Dwivedi and Mackey (2021, 2022); Shetty al. (2022); Hayakawa al. (2022). Since in what follows we aim to retrospectively improve MCMC output, where it is not unusual encounter n. 100_ 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 II that is straight-forward to sample using MCMC and is appropriately over-dispersed relative P. This should be achieved while remaining in the framework of kernel discrepancies, that optimal weights be explicitly computed, and coupled with sparse approximation that has low computational and storage cost.
3 Methodology
The methods that consider first sample states -.1.....,,, ..,. using IT-invariant MCMC, then post- process these states using kernel discrepancies (Section 2.2) and sparse approximation (Section 2.3), obtain an approximation to the target P. A variational argument, which we present Section 3.1, provides a suitable n-indeppenent choice for II (which agrees with our intuition from Section 2.1 that should be in some appropriate sense over-dispersed with respect P). Sufficient conditions for strong consistency of the approximation are established in Section 3.3.
3.1 Selecting II
Here present heuristic argument for particular choice of II; rigorous theoretical support for Stein II-mmortance Sampling is then provided in Section 3.3. Our setting is that of Section 2.2, and the following will additionally be assumed:
Assumption 1. It is assumed that
(A1) C7 = infaex k(x,x) >
(A2) C2 = Vkz,,) dP(x)) < <0.
'It is difficult to quantify the complexity of numerically solving (3), since this will depend on details of the solver and the tolerance that are used. On the other hand, we ignore the non-negan non-negativity and normalisation constraints, then we can see that the computational cost of solving the n dimeeniional linear system of equations is O(r²).

---

Note that (A2) implies that H((k) C (P),, and thus (1) is in fact strong (or Bochner) integral.
direct analysis of the optimal states associated to the optimal weights w* appears to be challenging due to the fact that the components of w² are strongly inter-dependent. Our solution here is to instead consider optimal states associated with weights that, while not optimal, much better than alternatives, with the advantage that their components are only weakly dependent. Specifically, we will be assuming that / is absolutely continuous with respect to II (denoted / < II), and study convergence of self-normalised importance sampling (SNIS), i.e. the approximation Pn P wx;(x;;) w/ x (dP/dII) where 21...., In II are independent. Since w and 1 w = from the optimality of w" under these constraints we have that (p(P%) I (P)).) 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 wn in general.
be expected to perform
The analysis of SNIS weights w is tractable when viewed as approximation of the kernel mean embedding μP in the Hilbert space H(k). Indeed, recall that DP(P,) /IE./ vall//lllaaaa where E. = Vn(up, - uP). Then, following Section 2.3.1 of Agapiou et al. (2017), we observe that
1 Vn D== dl (2x) (b(,22)) prl
dP
En vn Σ x,) pr
(6)
1
E
dP am
(I,)
The idea is seek 1I for which the asymptotic variance of &n is small. Supposing that
(c) dll (x) oc,
(SI)
from the weak law of large numbers the denominator in (6) converges in probability to 1. Further supposing that
/ dP am
(+) - pr] al[(x) < oo,
(S2)
H(A)
from the Hilbert space central limit theorem the numerator in (6) converges in distribution to a Gaussian J Ei= (dP/dlI((z;) ()(, up] $ upl covariance operator defined via
C00,C where C H(()) + H(k) is the
(F,Cg)m(())) - (1: an (z) - up]
dp
9. an (2)(((( - up]
dllx))
H(())
H(k)
see Section 10.1 of Ledoux and Talagrand (1991). Thus, from Slutsky's lemma applied to (6), we conclude that &, (0,,C Recalling that (PP)())) I&a (Begg) and noting that the mean square of the limiting Gaussian random variable is tr(C). a natural idea is select the sampling distribution II such that tr(C) is minimised.
Fortunately, the trace of C can be explicitly computed. It simplifies presentation to restrict attention to a Stein kernel 1: = kp, for which (4p = O, giving tr(C) = I (z) for convenience have let kp(x) kp(2,, a). Assuming that P and IT admit densities p and on X Rd the variational problem wish to solve is
dll()), where
/
[x(x)²
kp(x) dr s.t. / T(2) dx 1, T(x)
arg min argmin #EQ
where be the set of positive measures on Rd 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
<(x) ((r()))))))
(8)
Note that the normalisation constant of n is C2 from (1), whose existence we assumed. Then we verify that (S1-2) in fact hold for this choice of II. Indeed,
(S1) / dP (x)² all(.) =c² /
1
< C2
all(.)
< 00
kp(x)
C
(S2) / बा kp))kp()) dll(x)) C2| / (kp(z) dP(x) C² 0c,
dP
5

---

100
(Langevia)
(KGM))
0.3
(KCM))
EjKD
10-1
Langevio)) (Langevin)
(KGM))
0.1
(KGMI)
(KGMS)
10-2
(KGM3)
3
10'
20
(1.1)
(1.2)
Figure 1: Illustrating our choice II in 1D. (a) The univariate target (black), and our choice of II 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 II-Importance Sampling using the Stein kernels from (a); in each case, KSD was computed using the same Stein kernel used to construct II. Solid lines indicate the baseline case of sampling from P, while dashed lines indicate sampling from II. (The experiment was repeated 100 times and standard error bars are plotted.)
which shows that we have indeed solved (7). The sampling distribution II we have obtained is characterised up to normalisation constant in (8), so just like P we can sample from II using techniques such MCMC. sampling (i.e. without self normaasaaaais see Lemma I of Adachi et al. (2022). The Stein kernel kp determines the extent to which TI differs from P, as we illustrate next.
interesting note that II is also optimal for standard importance
3.2 Illustration
For illustration, consider the univariate target P (black curve) in Figure 1.1, a 3-component Gaussian mixture model. Our recommended choice of IT in (8) is shown for both the Langevin-Stein kernel (purple curve) and the KGMs-Stein kernels with E {1,3) (green and blue curves). The Stein discrepancy corresponding the Langevin- Stein kernel provides control over weak convergence (i.e. convergence of integrals of functions that are continuous and bounded), while the KGMs-Stein kernel provides additional control over the convergence of polynomial moments up order s; full details about the construction of Stein kernels are contained in Appendix A. The Langevin and KGMI-Stein kernels have (x) 2 while the KGM3-Steer kernel has kp kp(z) x6 each case as |a| - oc, and thus greater over-dispersion results from use of the KGM3- Stein kernel. This over-dispersion is less p higher dimensions; see Appendix D.I.
To illustrate the performance of Stein II-Importanc Sampling, we generated a sequence (nn))et of independent samples from II. For each {1...., 100). the samples {x1,-.. (...-... .,. 2.,) were assigned optimal weights w" by solving (3), and the associated KSD was computed. As a baseline, we performed the same calculation using independent samples from P. Figure 1.2 indicates that, for both Stein kernels, substantial improvement results from the use of samples from II compared to the use of samples from P. Interestingly, the KGM3-Stein kernel demonstrated a larger improvement compared to the Langevin-Stent kernel, suggesting that the choice of II 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, = fi(x) + e; with (x) = (1 1,22), t i 5, i 1, ...... 10, with e; independent (000, 1). The parameter x = (I1, 222 was assigned a prior (0, 7). Data were simulated using x = (0, 0). The posterior distribution P is depicted in the leftmost panel of Figure 2, while choice of II 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 II target their mass toward regions where P varies the most. The reason for this behaviour is clearly seen for the Langevin-Stein kernel since
'The same holds for classical quantisation; c.f. Section 2.1.

---

P
(Langevin)
II (KGM3)
11 (Riemann)
22
2
2
Figure 2: Illustrating our choice of II in 2D. The bivariate target P (left), together with our choice of based the Langevin-Sue kernel (centre left), the SGM3-Stein kernel (centre right), and the Riemann-Stene kernel (right).
Algorithm II-Invariant Mctropollis-Ajjustee Langevin Algorithm (MALA) Require: Zo (initial state), (step size), M (preconditioner matrix), (chain length), kp (Stein kernel) 1: for i 1...., do
2: x' 21-1 + cM- ' log + SM log (o + 2 1/2
uD Zi ()
==((x_1)
3: L + log
p(a')
+ log
kn(x') (p((111)
- (11;;---(<e)1 + 1 / u( 11
p(
4: (og(()) < L then x; t- a': else F, f Ca-1:end
UD
u(10, 1))
5: end for
Algorithm 2 Stein I1-mmportance Sampling (SITI--MLAN
Require: {.1,..., J from Algorithm
1, kp (Stein kernel) I: w E arg win (ww, Kow) w 2 0, 1 --w=11 1}
[Kp]a, = ke (x(
kp(x) C1 + c2 log ((=) for some C1, C2 > O; Appendix C 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: Appendix A for full detail.. Results in Figure $2 demonstrate that Stein IT-Importance Sampling improves upon the default Stein importance sampling method (i.e. with IT and P equal) for all choices of kernel.
An additional illustration involving a GARCH model with d = parameters is presented Ap- pendix D.4, where the effect of varying the order of the KGM-Stein kernel is explored.
3.3 Theoretical Guarantees
The aim of this section is establish when post-processing of Il-invariant MCMC produces strongly consistent approximation of P, for our recommended choice II in (8). Our analysis focuses on MALA (Roberts and Stramer, 2002), leveraging the recent work of Durmus and Moulines (2022) to present explicit and verifiable conditions on P for our results to hold. 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 (3) (c.f. Section 2.3). The resulting algorithms, which call Stein IT-Importance Sampling -MTIS-MALA) and Stein I--Thmmning (SITT-MALA), are quite straight- forward and contained, respectively, in Algorithms and 3. The cinstrlycoossrrained quadratic programme in Algorithm 2 was sol using the Python v0. 4 packages qpsolvers vers v3 1..4.0 O and proxsuite te v0. 3...7 7. While it difficult analyse the computational complexity associated with these methods, we believe they are at worst 3
'The use geometric information may 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 (Gorham et al., 2019). However, the Riemann Stein kemel is associated with prohibitive computational cost: it included here only for academic interest.

---

Algorithm 3 Stein 11-Thinning (SUTT-MALA)
Require: {.1....,, Fn from Algorithm 1, m (number of samples to retain), kp (Stein kernel)
1: for ,...., do
2: t arg min se( k)((v) + j= (p(
3: end for
Let A < indicate that A - a positive semi-definite matrix for A, B E RTxn, For a symmetric positive definite matrix A let || ;= V/1A-12 for E R. Let C* (RR)) denote the set of s-times continuously differentiable real-valued functions on R
Theorem 1 (Strong consistency of SIISS- and SITT-MALA), Let Assumption I hold where kp is a Stein kernel, and let Dr (())) + [0, oc] denote the associated KSD. Assume also that
(A1) V logg E (RR) with SuPIEr V ((g](())] < 00
(A2) b > 0, Bı such that 7 log b for all ||>||| 2 Bı
(A3) kp (T..))
(A4) 30 < b2 < Cb, B2 such that 42kp(z) 3 b2l for all |||||| 2 B2
Let Pt - w. (z;)) be the result of running Algorithum 2 and let Pn,m= the result of running Algorithm Let m S n and m = (2(llg n) for some o > 2. Then there exists > 0 such that, for all step sizes e E (0,eo) and all initial states E Rd almost surely as m, n + 00.
أiー1 d(x) be
m
Dp(P=) - 0, Dp(Pn,m)
The proof is in Appendix B. Compared to earlier authors such as Chen et al. (2019); Riabiz et al. (2022), a major novelty here is that our assumptions are explicit and can often be verified (see also Hodgkinson al., 2020). (A2) is strong log-concavity of P when B O, while for B₁ 0 this condition is slightly stronger than the related distant dissipativity condition assumed in earlier work (Gorham and Mackey, 2017; Riabiz al., 2022). (A4) holds for the Langevin-Stenn kernel (i.e. weak convergence control) and for the KGM1-Stein kernel (i.e. weak convergence control + control over first moments), but for the higher-order KGM-Stein kernels. Extending our proof strategy the higher-order KGM-Stein kernels would require further research into the convergence properties of MALA, and this is expected be difficult.
4 Benchmarking on PosteriorDB
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 (Chopin and Ridgway, 2017). To introduce objectivity into our assessment, we exploited the recently re- leased PosteriorDB benchmark (Magnusson et al., 2022). This project an attempt toward standardised benchmarking, consisting of collection of posteriors to be numerically approxi- mated, Here, we systematically compared the performance SIII--ALA against the default Stein importance sampling algorithm (i.c. II = P: denoted SIS-MALA), and also against un- processed P-invariant MALA (i.e. uniform weights), reporting results across the breadth of defined in the Stan probabilistic pro-
PosteriorDB. The test problems in PosteriorDE gramming language, and so BridgeStan (Roualdes et al., 2023) was used to directly access posterior densities and their gradients as required. For all instances of MALA, adaptive algorithm was used to learn a suitable preconditioner matrix M during the warm-up period; Appendix D.3. All experiments that we report be reproduced using code available at https://githor com/ wwwgy//wag//teeeeeeeeeee Stwperannn/tteeiiiiieieie Ini-mmortancee--e Serpoortanne-Smmplii...
Results reported in Table for 3 x 10x samples from MALA. These focus the Langevin- Stein kernel, for which our theory holds, and the KGM3- kernel, for which it does not. There was significant improvement of SIIIS- 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
'These earlier results required high-level assumptions on the convergence of MCMC, for which explicit sufficient conditions had yet to be derived.
8

---

Langevin Stein Kenee
KGMS- Siein Keme!
SIS
SITS
SIS
SIIS
MALA
MALA
MALA
MALA
MALA
cariins jogr-g_reee
1,11
0.0874
0.0332
5.33
0.856
0.181
gp_pois, Jegr- kidiq kidscore, mombs
3
0.298
0.0436
0.0373
1.22
0.385
0.223
1.04
0.109
0.0941
4.66
0.848
0.476
kidig kidseore monig
0.516
1.55
cenequite-logmesquiie Jogvvlllueee Jogvolume
3
1.10
0.179
0.156
4.97
1.70
0.844
4.47
1.09
1.01
8,91
6.03
curnings- logeass--ogam..mmmm Joghhegg logheight_male
9.46
1.96
1.59
53.9
15.4
8.65
surch-garchll
4
0.543
0.159
0.130
4. 4.70
1.16
1.01
kidiq-kidkoore mombsiq
0.982
29.5
5,05
curnings logeam, interaction_
5
3,09
1.36
1.33
19.3
10.4
8.94
kidiq- bedig-kidscore_intee interrccionnnn
7.74
1.65
1.79
10.1
kidiq, with mom, work kidscore interaction c
5
1.35
0.639
0.711
7.92
4.05
4.17
kidiq, with mom workm--or-i kidscore_ interaction c2
1.38
0.689
0.699
8.09
4.24
4.25
kidbq with mom, work kork-iideorre, jateraction
1.11
0,500
0.499
6.62
2,63
3.25
kidiq_with, with mom, work kidscore, mom, work
1.07
0,5417
0.545
6.70
2,63
3.04
low dim, gauss, mix llo dau.ggussss._aooooooooooauuusssd
1.87
mesquise-logmesse meguuquite logya
5
1.83
0.821
0.818
12.6
5.73
5.59
hmmm, example-hrmm example
1.99
0.578
0.523
4.13
3.40
sbirc-blr
479
154
134
3300
1100
854
shiri-hir
6
201
66.7
60.3
13:40
514
595
mesquite-logmeeet logmesquite, Jogyashh
7
1.89
1,18
1.23
15.5
8.88
10.1
hhall, drive event 0-hmm drive
1,15
0.679
0.698
8.55
4.72
3.99
bhall, eventt_ I-hmm drive
8
42.9
11.9
12.4
285
85.6
67.8
hudson_ Jynx lode- volterra
4.62
2.29
2.15
47.4
18.8
18.9
mesquite- mmsquite-ognessuui
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-mssque
0.268
0.235
eight schools--ighh scchsols_heeeeeeeeee centered
10
0.526
0.100
0.182
7.53
2.15
215
eight, seohttschholseiiee schools,
schnols__nncenterred noncentered
0,210
0.137
0.137
43.6
28.7
27.5
nos 1972-nss
3.89
3.45
72.9
36.2
34.4
nex1976-nes
10
6.67
3,86
3,53
77.5
35.5
34.4
1980-n9
25.7
nes1984-nes
10
6.18
3.75
3,43
71.3
34,9
33.6
nes 1988
10
7.40
3.70
3.27
81.4
34,6
32.4
nes 1992-nes
10
7.52
4.32
3,84
89.1
39.7
37.3
nos 1996-925
10
6.44
3.87
3,53
74.1
36.4
34.3
nes2000-nes
3.35
2.22
2.20
21.3
22.8
caanondd-diamondd
26
196
157
143
5120
2990
2620
meycle sp-accel _pp
66
8.25
966
623
815
Table 1: Benchmarking PosteriorDB. Here we compared raw output from MALA with the post-processed output provided by the default Stein importance sampling method of Liu and Lee (2017) (SIS-MALA and the proposed Stein II-Importance Sampling method (STHS--MAL) Here dim( P) and the number MALA samples was n = 10 The Langevin and KGM3-Stein kernels were used for SIS-MALA and LIIS-MALA and the associated KSDs are reported. Ten replicates were computed and statistically significant improvement is highlighted in bold.
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 of the target increased, supporting the intuition that we set out earlier and in Appendix D.1. An in-depth breakdown results, including varying the number n of samples that were used, and the performance SITT-MALA, can be found in Appendices D.5 and D.6.
If P and its gradients are cheap evaluate, the computational SIS-MALA, and one could run more iterations of MALA for an equivalent computational cost. But for more complex P, the computational cost of all algorithms will be gated by the number of times P and its gradients need to be evaluated, making the direct comparison in Table 1 meaningful. Further, we aim for compressed representation of P, then some form of post-processing of MALA would be required, which would then entail an additional computational cost.
of MALA lower than that of
Our focus is the development of algorithms for minimisation of KSDs; the properties of KSDs themselves are out of scope for this work?. Nonetheless, there is much interest in better understanding the properties of KSDs, and therefore also report performance of -MLIA in terms of 1- Wassersteei divergence in Appendix D.7. The main contrast between these results and the results in Table 1 that, being score-based, KSDs suffer from the blindness mixing proportions phenomena which has previously been documented in Wenliang and Kanagawa (2021); Koehler et al. (2022); Liu al. (2023). Caution should therefore be taken when using algorithms based on Stein discrepancies
The interested reader is referred Gorham and Mackey (2017); Barp et al. (2022b); Kanagawa et al. (2022).

---

in the context of posterior distributions with multiple high probability regions that are spatially separated. This is also 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 11, with Ti (3:) x p(x)d/(d++1) which provides generic form of over- dispersion and is optimal for approximation in I- I-Wasserstein divergence (c.f. Section 2.1), was also considered. Results in Appendix D.8 indicate that, while II, yields an improvement compared the baseline of using P itself, 1I, may be less effective than our proposed II when P is skewed.
5
Discussion
This paper presented Stein II-Importance Sampling: an algorithm that is simple to implement, admits an end-to-end theoretical treatment, and achieves significant improvement over existing post- processing methods based on KSD. On the negative side, second order derivatives of the statistical model are required, and we are ultimately bound the performance of the KSD on which Stein II- I-Immortance Sampling is based. Our analysis focused on MALA, but there is in principle no barrier to deriving sufficient conditions for consistent approximation that are applicable to other sampling algorithms, such as the unadjusted Langevin algorithm. Of course, it remains to be seen whether SIS-MALA or any of its variants will stand the test of time compared to continued development in MCMC methodology, but we believe this line of research merits further investigation. For models for which access to second order derivatives impractical, our methodology and theoretical analysis are directly applicable to gradient-free KSD (Fisher and Oates, 2023), and this would be an interesting direction for future work. Similarly, alternatives to KSD that are better-suited to high-dimensional could be considered, such as the sliced KSD of Gong et al. (2021a,b).
Acknowledgements CW was supported by the China Scholarship Council. HK and CJO were supported by EP/W019590/1. The authors are grateful to Fravçois-Xavier Briol for feedback on an earlier draft of the manuscript, and to the anonymous Reviewers for their input.
References
Adachi, M., Hayakawa, S., Jørgensen, M., Oberhauser, H., and Osborne, M. A. (2022). Fast Bayesian inference with batch Bayesian quadrature via kernel recombination. In Proceedings of the 35th Conference Neural Information Processing Systems.
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.
Anastasiou, A., Barp, A., Briol, F.-X., Ebner; B., Gaunt, R. E., Ghaderinezhad, F., Gorham, J., Gretton, A., Ley, C., Liu, Q. Mackey, I.., Oates, C. J., Reinert, G., and Swan, Y. (2023). Stein's method meets statistics: A review of some recent developments. Statistical Science, (38):100-139.
Bach, F. (2017). On the equivalence between kernel quadrature rules and random feature expansions. The Journal of Machine Learning Research, 18(1):714-751.
Bach, F., Lacoste-Julien, S., and Obozinski, G. (2012). On the equivalence between herding and conditional gradient algorithms. In Proceedings of the 29th International Conference on Machine Learning.
Barbour, A. D. (1988). Stein's method and poisson process convergence. Journal of Applied Probability, 2554):175-184.
Barbour, A. D. (1990). Stein's method for diffusion approximations. Probability Theory and Related Fields, 84(3)2297-322.
Barp. A., Oates, C. J., Porcu, E.,, and Girolami, M. (2022a). Riemann-Stein kernel method. Bernoulli, 28(4):2111-2288.
Barp, A., Simon- Gabriel,, C.-J,, Girolami, M., and Mackey, L. (20226). Targeted separation and convergence with kernel discrepancies. arXiv preprint arXiv:2209.12833.
10

---

Bénard, C., Staber, B., and Da Veiga, S. (2023). Kernel stein discrepancy thinning: A theoretical perspective of pathologies and practical fix with regularization. arXiv preprint arX1v:2301.13528.
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 Machine Learning, pages 586-595.
Briol, F.-X., Oates, C. J., Girolami, M., Osborne, M. A., and Sejdinovic, D. (2019). Probabilistic integration: A role in statistical computation (with discussion and rejoinder). Statistical Science, 34(1):1-22.
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.
Chen, W. Y., Barp, A., Briol, F.-X., Gorham, J., Girolami, M., Mackey, L.,, and Oates, J. (2019). Stein point Markov chain Monte Carlo. In Proceedings of the 36th International Conference on Machine Learning, pages 1011-1021.
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.
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.
Chopin, N. and Ducroeq, G. (2021). Fast compression of MCMC output. Entropy, 23(8):0017.
Chopin, N. and Ridgway, J. (2017). Leave pima indians alone: Binary regression as a benchmark for Bayesian computation. Statistical Science, 32(1):64-87) 87.
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.
Cohort, P. (2004). Limit theorems for random normalized distortion. The Annals of Applied Probability, 14(1):118-143.
Durmus, A. and Moulines, É. (2022). On the geometric convergence for MALA under verifiable conditions. arXiv preprint STti:220100195//1
Dwivedi, R. and Mackey, (2021). Kernel thinning. In Proceedings of 34th Conference on Learning Theory, pages 1753-1753.
Dwivedi, R. and Mackey, L. (2022). Generalized kernel thinning. In Proceedings of the 10th International Conference on Learning Representations.
Ehler, M., Gräf, M., and Oates, C. J. (2019). Optimal Monte Carlo integration on closed manifolds. Statistics and Computing, 29())1223-112)
Fisher, M. A. and Oates, C. J. (2023). Gradient-fres kernel Stein discrepancy. In Proceedings of the 37th Conference on Neural Information Processing Systems.
Girolami, M. and Calderhead, B. (2011). Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 1332):1
Gong, W., Li, Y., and Hernández-Lobato, J. M. (2021a). Sliced kernelized Stein discrepancy. In Proceedings of the 9th International Conference on Learning Representations.
Gong, W., Zhang, K., Li. Y., and Hernández-Lobato. M. (2021b). Active slices for sliced Stein discrepancy. arXiv:2102.03111
Gorham, J., Duncan, A. B., Vollmer, J., and Mackey, L. (2019). Measuring sample quality with diffusions. The Annals of Applied Probability, 29(5):2884 2928.
Gorham, and Mackey, L.. (2015). Measuring sample quality with Stein's method. In Proceedings of the 29th Conference Neural Information Processing Systems, pages 226-234.
11

---

Gorham, J. and Mackey, L. (2017). Measuring sample quality with kernels. In Proceedings of the 34th International Conference on Machine Learning, pages 1292-1301.
Gotze, F. (1991). On the rate of convergence in the multivariate CLT. The Annals of Probability, pages 724-739.
Graf, S. and Luschgy, H. (2007). Foundations of Quantization for Probability Distributions. Springer.
Hawkins, C., Koppel, A., and Zhang, Z. (2022). Online, informative MCMC thinning with kernelized Stein discrepancy. arXiv preprint arXi::220..07.300
Hayakawa, S., Oberhauser, H., and Lyons, (2022). Positively weighted kernel quadrature via subsampling. In Proceedings of the 35th Conference on Neural Information Processing Systems.
Hodgkinson, I.., Salomone, R., and Roosta, F. (2020). The reproducing Stein kernel approach for post-hoc corrected sampling. arXiv preprint arXiv:2001.09266.
Kanagawa, H., Gretton, A., and Mackey, L. (2022). Controlling moments with kernel Stein discrep- ancies. arXiv preprint arXiv:2211.05400 (vl).
Karvonen, T., Oates, C. J., and Girolami, M. (2021). Integration in reproducing kernel Hilbert spaces of Gaussian kernels. Mathematics of Computation, 90(3311:2209-222) -2333.
Karvonen, T., Oates, C. J., and Sarkka, S. (2018). A Bayes-Sard cubature method. In Proceedings of the 32nd Conference Neural Information Processing Systems.
Kent, J. (1978). Tive-ieversible diffusions. Advances in Applied Probability, 10(4)) 835.
Koehler, F., Heckett, A., and Risteski, A. (2022). Statistical efficiency of score matching: The view from isoperimetry. In Proceedings of the 36th Conference Neural Information Processing Systems.
Ledoux, M. and Talagrand, M. (1991). Probability in Banach Spaces: Isoperimetry and Processes, volume 23. Springer Science & Business Media.
Liu, Q. and Lee, J. (2017). Black-box importance sampling. In Proceedings of the 20th International Conference on Artificial Intelligence and Statisties, pages 952-961.
Liu, Q., Lee, J., and Jordan, M. (2016). kernelized Stein discrepancy for goodness-of-fit tests. In Proceedings of the 33rd International Conference on Machine Learning, pages 276-284.
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 Learning.
Magnusson, M., Bürkner, P., and Vehtari, A. (2022). PosteriorDB: A set of posteriors for Bayesian inference and probabilistic programming. https: //t:://///ttbb http://ggtttbb.cco/tttaede/////e wtw.sttaddde/pttterrrrreeeee
Meyn, S. P. and Tweedie, R. (2012). Markov Chains and Stochastic Stability. Springer Science & Business Media.
Oates, C. J., Girolami, M., and Chopin, N. (2017). Control functionals for Monte Carlo integration. Journal of the Royal Statistical Society: Series (Statistical Methodology), 79(3):695-718.
Riabiz, M., Chen, W., Cockayne, J., Swietach, P., Niederer, S. A.. Mackey, L., and Oates, J. (2022). Optimal thinning of MCMC output. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 84(4): 1059-10111111
Roberts, G. O. and Rosenthal, J. S. (1998). Optimal scaling of discrete approximations Langevin diffusions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255 268.
Roberts, G. O. and Stramer, (2002). Langevin diffusions and M algorithms. Methodology and Computing in Applied Probability, 4:337-357.
12

---

Roualdes, E., Ward, B., Axen, S., and Carpenter, B. (2023). BridgeStan: Efficient in memory access to Stan programs through Python, Julia, and R. httpp:://githiiee wwt..//oww.eee/r.ieeeeeeeoeeeeee
Shetty, A., Dwivedi, R., and Mackey, L. (2022). Distribution compression near-linear time, In Proceedings of the 10th International Conference on Learning Representations.
Sonnleiterr M. (2022). The Power of Random Information for Numerical Approximation and Integration. PhD thesis, University of Passau.
Teymur, O., Gorham, J., Riabiz, M., and Oates, J. (2021). Optimal quantisation of probability
measures using maximum mean discrepancy. In Proceedings of the 24th International Conference on Artificial Intelligence and Statistics, pages 1027-1035.
Wenliang, L. K. and Kanagawa, H. (2021), Blindness of score- based methods to isolated components and mixing proportions. In Proceedings of "Your Model is Wrong" @ the 35th Conference on Neural Information Processing Systems.
13

---

Appendices
These appendices contains supporting material for the paper Stein 11-Importance Sampling. The mathematical background on Stein kernels is contained in Appendix A. The proof of Theorem is contained in Appendix B. For implementation of Stein I1-Importance Sampling without the aid of automatic differentiation, various explicit derivatives are required; the relevant calculations can be found in Appendix C. The empirical protocols and additional empirical results are presented in Appendix D.
A Mathematical Background
This appendix contains mathematical background on reproducing kernels and Stein kernels, as used in the main text. Appendix A.1 introduces matrix-valued reproducing kernels, while Appendix A.2 specialises to Stein kernels by application of a Stein operator to matrix-valued kernel. A selection of useful Stein kernels are presented in Appendix A.3.
A.1 Matrix- Valued Reproducing Kernels
A matrix-valued kernel is function K Rd x Rd Rdxd that is both
1. symmetric; (xx, v) = K(y, x) for all x, E IRd, and
2. positive semi-definite; Σ (emncc, (Z,, 2 for all 21,...,fn E Rd and CI...., E Rdd
Let K2 K(; r). For vector-valued functions g.g' Rdd + Rd. defined by = Ei K G and g' = 2 j= Kgz ci define an inner product
m ΣΣ
1
(g.g') (
Cc, 2)),))))))) =1
(9)
There unique Hilbert space of such vector-valued functions associated to K, denoted H(K); see Proposition 2.1 of Carmeli et al. (2006). This space is characterised as
H(K) span x, IR
where here the closure is taken with respect to the inner product (9). It can shown that H(K) in fact a reproducing kernel Hilbert space (RKHS) which satisfies the reproducing property
(g. Kc)) (g(z).c)
for all H(K) and x, c Rd Matrix-valued kernels are the natural starting point for construction of KSDs, as described next.
A.2 Stein Kernels
general construction for Stein kernels is to first identify matrix-valued RKHS H(K) and an operator Sp H(K) + L (P) for which f S,h dP = 0 for all h E H(K). Such an operator will be called a Stein operator. The collection {Sph : H(K)) inherits the structure of RKHS, whose reproducing kernel
kp(a, v) (Sp Sp3KK,SpK,K, Spky) (y)((())) H(K)
(10)
Stein kernel, meaning that μı where μaP is the kernel mean embedding from (1); see Barp et al. (2022b). Explicit calculations for the Stein kernels considered in this work can be found in Appendix C.
For univariate distributions, Barbour (1988) proposed to obtain Stein operators from infinitessimal generators of P--invariant continuous-time Markov processes; also Barbour (1990); Gotze (1991).
14

---

The approach was extended to multivariate distributions in Gorham and Mackey (2015). The starting point is the P-invariant Itō diffusion
1 V. (p(X.) M(.)) dx;
(11)
p(,))
where p is the density of P, assumed to be positive, M : + Rdxd symmetric matrix called the diffusion matrix, and W, is standard Wiener process (Kent, 1978; Roberts and Stramer, 2002). Here the notation VV Al; = V (A..) indicates the divergence operator applied each row of the matrix A(x) E Raxd The infinitessimal generator is
(Apu)(x) = 2 p(x)
V (=())(())Vv(r))
Substituting h(z) for Vu(ar), obtain Stein operator
(Sph)(z)
V [v()) p(x)
(12)
called the diffusion Stein operator (Gorham et al., 2019). This is indeed Stein operator, since under mild integrability conditions on K, the divergence theorem gives that Sph dP = 0 for all h HK));; for full details and proof see Barp et al. (2022b).
A.3 Selecting Stein Kernel
There are several choices for a Stein kernel, and which we should use depends on what form of convergence hope control (Gorham and Mackey, 2017; Gorham et al., 2019; Hodgkinson et al., 2020; Barp et al., 2022b; Kanagawa et al., 2022). Appendix A.3.1 describes the Langevin-Ster kernel for weak convergence control, Appendix A.3.2 describes the KGM-Stein kernels for additional control over moments, and Appendix A.3.3 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 be specified. As reasonably automatic default we define
F. arg max p(z), -- -22 log p(z.)
as location and a matrix of characteristic length scales for P that will be used throughout. These values can typically obtained using gradient-based optimisation, which is usually cheaper to perform compared to full approximation of P. It is assumed that 4 log p(x.) is positive definite in the sequel.
A.3.1 Weak Convergence Control with Langevin-Stein Kernels
The first kernel consider, which called the Langevin-Stein kernel in the main text, was introduced by Gorham and Mackey (2017). This Stein kernel was developed for the purpose of controlling the weak convergence of a sequence (Q.) P(Rd) to P. Recall that sequence (Qn)eeN is said to converge weakly (or in distribution) to P Pif fdQ,, + fdP for all continuous bounded functions f Rd + R.. This convergence is denoted Qn
in shorthand.
The problem considered in Gorham and Mackey (2017) was how to select combination of matrix- valued kemel K (and, implicitly, a diffusion matrix M) such that the Stein kemnel kp (10) generates a KSD (Q) in (4) for which (Q().) - 0 implies Qn 4 P. Their solution was to combine the inverse multi-quadric kernel with identity diffusion matrix;
K K(x, 3) = (1 + ||. - - v1
M(x)=1 =
for E (0, Provided that P has a density p for which Voog Vlogp(x:) is Lipschitz, and that P is distantly dissipative (see Definition 4 of Gorham and Mackey, 2017), the associated KSD enjoys weak convergence control. Technically, the results Gorham and Mackey (2017) apply only when == T. but Theorem 4 in Chen et al. (2019) demonstrated that they hold also for any positive definite Σ. Following the recommendation several previous authors, including Chen et (2018, 2019); Riabiz et al. (2022), we take B throughout.
15

---

A.3.2 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 f f dP for which the integrand f Rd - R. is unbounded. In particular, this is the case for moments of the form f((:) x ... xa 0 o may seek also the stronger property moment convergence control. The development of KSDs for moment convergence control was recently considered by Kanagawa et al. (2022), 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
Ng. In such situations, we
kernel.)
sequence (Qn)nEN C. ((RR) said to converge P in the sth order moment if [|l==|ldQQ((a) + [||==|| To establish convergence of moments, we need an additional condition on top of weak convergence control: uniform integrability control. A sequence of measures (Qn) is said to have uniformly integrable sth moments if for any > O, we can take >> 0 such that
sup n ENN ||||>>
| dQn(r)= <.
This condition essentially states that the tail decay of the measures is well-controlled (so that it has convergent moment). The KSD) convergence Dr(Qn) -> implies uniform integrability for any e > 0, we can take r. 0 and fa (H(K) such that
(p) 2 1 > rz) -E,
(13)
i.e., the Stein-modified RKHS can approximate the (norm-weighted) indicator function arbitrarily well. Such a function fa can be explicitly constructed (while not guaranteed to be a member of the RKHS). Specifically, the choice fa (1 'z)g satisfies (13) under an appropriate dissipativity condition, where & is a differentiable indicator function vanishing outside a ball, and g(x) = -x/ I | This motivated Kanagawa al. (2022) to introduce the sth order KGM-Stein kernel, which is based on the matrix-vasn valued kernel and diffusion matrix
K(x = [ || >/|x>) + Min((: v)] I, M(x) (1 + 1| x, (1:) 2 1
where (z + क1( )||x) is Cj universal kernel (see Barp et al., 2022b, Theorem 4.8). For comparability our results, we take & to be the inverse multi-quuadic ((r) = (1 + /²) 1/2, and
1+ ( x.) x.) x.)
Klin (2r, y)
v1 -z- 1 / 1.
Here the normalised linear kernel Kiin ensures g H(K), while the Co universal kernel Φ allows approximation of Spteg; see Kanagawa et al. (2022).
A.3.3 Exploiting Geometry with Riemann Laaneenn-nnnnnnn--nnnneen Kernels
For academic interest only, here we describe the Riemann-Ster Stein kernel that featured in Figure 2 of the main text. This Stein kernel is motivated by the analysis of Gorham et al. (2019), who argued that the use of rapidly mixing Itô diffusions in Stein operators can lead to sharper convergence control. The Rieman kernel is based on the class of so-called Riemannian diffusions considered in Girolami and Calderhead (2011), who proposed to take the diffusion matrix M in (11) to be M = (Iprior Iiisher the inverse of the Fisher information matrix. Fisherrs regularised using the Hessian the negative log-prior, Ipnorr For the two-dimensional illustration in Section 3.2, this leads to the diffusion matrix
M(x) = I+
Σ V((W)))))))) 1=1
where we recall that = f(()) + E,, where the are independent with ~ (O, 1), and the prior x (O 1). For the presented experiment we paired the above diffusion matrix with the inverse multi-quadric kernel K [(x, 3) (1 ||æ y v for B The Riemann-Ster kernel extends naturally to distributions P defined on Riemannian manifolds X; see Barp et al. (2022a) and Example I of Hodgkinson et al. (2020).
Unfortunately, the Riemann-Steen -Sein kernel prohibitively expensive most real applications, since each evaluation of M requires a full scan through the size-n dataset. The computational complexity
16

---

of Stein II- T--Thinning with the Riemann-Ster kernel is therefore O(m²n²), which is unfavourable compared the (m)2n) complexity in the case where the Stein kernel is not data-dependent. Fur- thermore, 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.
B Proof of Theorem 1
This appendix is devoted to the proof of Theorem 1. The proof is based the recent work of Durmus and Moulines (2022), on the geometric convergence of MALA, and on the analysis of sparse (greedy) approximation of kernel discrepancies performed in Riabiz et al. (2022); these existing results recalled in Appendix B.I. An additional technical result on preconditioned MALA contained in Appendix B.2. The proof of Theorem itself is contained in Appendix B.3.
B.1 Auxiliary Results
To precisely describe the results on which our analysis is based, we first need to introduce some notation and terminology. Let X [1, 2o) and, for a function f X - IR and a measure μ on X, let
[f(x]]
1 fdu 11livs
I == sup =x V(z)
Mully == sup
Recall that Q-invariant Markov chain (æi)ee C X' with n step transition kernel Q" is V-uniformly ergodic (see Theorem 16.0.1 of Meyn and Tweedie, 2012) if and only if 3R E (o), E (0, such that
(x, Qllv Rp"V(x))
(14)
for all initial states E X and all E N.
Although MALA (Algorithm 1) is classical (Roberts and Stramer, 2002), until recently explicit sufficient conditions for ergodicity of MALA had not been obtained. The first result will need is due Durmus and Moulines (2022), 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 Algorithm is the identity matrix. The extension of this result to preconditioned MALA will be handled in Appendix B.2.
Theorem 2. Let Q E ((Rd) admit density, q, such that
(DMI) there exists 2o with V log q(z0) =0
(DM2) q is twice continuously differentiable with supzerd 2 log æo)1l < 0
(DM3) there exists > 0 and B 2 0 such that -V2 22 zo) ł bl for all | - -zoll B.
Then there exists > 0 such that for all step sizes E (0, standard Q- D-invariant MALA (i.e. with M T) is V-uniformly ergodic for V(x) exp( || 1o)l²))
Proof. This Theorem of Durmus and Moulines (2022).
The next result that we will need establishes consistency of the greedy algorithm applied to samples from Markov chain that is Q-invariant.
Theorem 3. Let P. Q E P(X) with P < Q. kp X x X - IR be Stein kernel and let
Dp : X X + [0, oo] denote the associated KSD. Consider Q-invariant, time-homogeneout Markov chain (Ią)sEN C X such that
(R+1) (R.)een is V V-uniformly ergodic, such that V(x) 2 aQ (p)(k()))
(R+2) sup
<0
17

---

(R+3) there exists y such that supien exp max (1 dā()))) kp(z:)) cc.
Let Pn.m be the result of running the greedy algorithm in (5). If S n and log(n) = 0(m7/3) for
some y < 1, then Dp(Pn,m + 0 almost surely m, nt OC.
Proof. This Theorem of Riabiz et al. (2022).
B.2 Preconditioned MALA
In addition to the auxiliary results in Appendix B.1, 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 Metsopolis-Hastings algorithm whose proposal is the Euler-Maruyama discretisation of the Itô diffusion (11).
Proposition 1. Let (x) = M for a symmetric positive definite and position- independent matrix M E Rdxdd Let Q E P(R4) admit a probability density function (PDF) q for which the Q-invariant diffusion (X)) given by setting p q in (11), is well-defined. Then under the change of variables
Y, = M
dY, ( 2 (Y)))) dW,
1
(15)
where ((x) x q(M-1/2)) for all x E Rd.
Proof. From the chain rule,
(V log (v) = V log (M) = (ogg)//
and thus, substituting Y. 1/ x,, (15) is equal to
dx, = M-1
M 1/2 (1/22((oogg log (
dd,
2
M (Vlog)
log -M 1/2
M
which is identical to (11) in the case where M(x) M is constant.
Let Q and @ be the distributions referred to in Proposition 1, whose PDFs are respectively ((x) and ((x)) & q( qM-11/23) Proposition then implies that the M-preconditioned MALA algorithm applied Q (i.e. Algorithm I for II Q) is equivalent to the standard MALA algorithm (i.e. M = applied Q. This fact allows to generalise the result of Theorem 2 as follows:
Corollary 1. Consider a symmetric positive definite matrix M E Rddd Assume that conditions (DM1-3) Theorem 2 are satisfied. Then there exists e > and > such that for all step sizes E 6)), the M M-preconditioned Q-invariani MALA is V-uniformly ergodic for V V(x)
1' || - =0.l²).
exp To
Proof. From Theorem 2 and Proposition the result follows if we can establish (DM1-3) for Q since M preconditioned MALA is equivalent to standard MALA applied to Q. For matrix A Rdxd let Amin (A) and Amax (A) respectively denote the minimum and maximum eigenvalues of A. For
(DM1) we set yo = M1/220 and observe that
(V log (vo) = - 112 log (zp) = 0.
For (DM2) we have that
2 (y vo)l = sup
1 log
sup
yE3d
yERd
< Amin (M
sup
|| lo - -z0)l < co.
ZERR
18

---

For (DM3) we have that
-(V2 log (v y0) - (M (M (y W0))
-M-
1- M (bT) bM 1 > B'T
where 6' = bdmuux (M) which holds for all || - zoll 2 B, and in particular for all lly - yoll 2 B'
where B' = B Xmaxx (M)1/2 Thus (DM1-3) are established for Q.
Remark 1. The choice M- =-1, which sets the preconditioner matrix M equal to the inverse of the length scale matrix E used in the specification of the kernel K (c.f. Appendix A.3), leads to the elegant interpretation that Stein IT-Importance Sampling applied to M M-precoonitiioee MALA is equivalent to the Stein II-mmporaance Sampling applied to standard MALA (i.e. with M for the whitened target P with PDF p(z) & ppMpp1///p For our experiments, however, the preconditioner matrix M was learned during a warm-up phuse of MALA, since general the curvature of P (captured by E) and the curvature of TT (captured by M-1) may be different.
B.3 Proof of Theorem
The route to establishing Theorem 1 has three parts. First, establish (DM1-3) of Theorem 2 with Q = II, to deduce from Corollary 1 that II-ivarriant M -rpreccondttioeed MALA is V-uniformly ergodie. This in turn enables us to establish conditions (R* 1-3) of Theorem 3, again for from which the strong consistency Dr(Pa,m) - of SHIT-MALA is established. Finally, we note that 0 < (P(( < Dy(Pn.m since the support of Pr.m is contained in the support of PA and the latter is optimally weighted, whence also the strong consistency of -TIISMALA.
II.
Establish (DM1-3) First we establish (DM1-3) for Q = II. Fix xo E Rd. For (DM2), first recall that the range of kp is [C7, oo) where C1 >> from Assumption 1. Since log(-) has bounded second derivatives on [CP, 0o), there is constant C > 0 such that
Va E Rd.
2 log (x))l) S C()V))((a)))
Thus, using compactness of the set {a || zoll S B2).
sup FER"
V log kp)())] < C max
sup \p'tl( |0-2al <B
sup 1V=kp((t) [ B2
<c. (16)
<de by (A3)
<b2ll7] by (A4)
Now, π is twice differentiable as is the product of twice differentiable functions p and kpp from (A1) and (A3), and moreover
sup zERd
2 log (( - xo)1l S sup reRd
sup || || log kp(x)) < co, V logp(t))
+
aeRd
<00 by (A1)
< by (16)
so (DM2) is satisfied. For (DM3), first note from the chain and product rules that for all |||||| 2 Bg
V2kp(s mo)
Vkr(w [o) to)]
b2
2 kog x0)=
(17)
kp(a- - No)
kp(z - 20)²
C
3(b2/C9)1 by (A4)
20
Thus, for all ||0 xoll 2 B = ||00|| + max B1, B2).
1
by 20
-V -Vo (zz) - 2 log - xo)
2
log kp (z x0)
bi
I
(18)
2
Ib1 by (A2)
30b///9 (17)
=
as required. The same argument establishes (DM1); )) from (18) we have li T(x) = 0. and since π a continuously differentiable density there must exist an 2o at which π is locally minimised. Thus we have established (DM1-3) for Q 11 and we may conclude from Corollary that there is an é > 0 and b' > such that, for all e E (0, <6), the il-invariant M prreconditioned MALA chain (æa)EN is V --unffommly ergodic for V(x) = C2 exp f. || -=ol?) (since if Markov chain is V-uniformly ergodic, then it is also CV CV-unifomly ergodic).
19

---

Establish 1-3) The aim is now to establish conditions (R 1-3) of Theorem for Q II. By construction dP/dII C2/ ( < C2/C1 < oc, where C₁ and C2 were defined in Assump- tion 1, so that P < II. It has already been established that (αi)keN is V-uniformly ergodic, and further
dP V(x) C2 exp (16 zoll?) 2 2 dn (a) (r(x)
V/
for all z, which establishes (R+1) 1). Let R and denote constants for which the V-uniform ergodicity property (14) satisfied. From V-uniform ergodicity, the integral V dII exists and
dP dll
(a))k(()))Vas)) -Ca / dlI C2 [[](( / dnl
E
S CP + 0
which establishes (R+ Fix y 0. By construction dP/dIl C2/Ci, and thus
exp {
dP max (I. (x)² Ap(c)) exp (7kp(x))
where ÿ = maax(1 Since we have assumed that kp is continuous with, from (A4),
ba = sum kp(x) oo,
<
+ ||=|
may take y such that 70g // //16, that || | exp(ik((t))] < and in particular
E (eexp(ikp(za)) - exp(Skp(x)) dT(())) = =xp x Rx"V(ao) + 0
which establishes 3). Thus we have established (R+ 1-3) for 11. so from Theorem we have strong consistency of SITT-ALA (i.e. Dp(Pn.m) + 0) provided that m S with log(n) = 7/2 for some 7 < 1. The latter condition is equivalent to = ?((log n))) for some > which we used for the statement. Since (p(Pt) S (p(Pn, (Pn.)) the strong consistency of SIIIS MALA is also established.
C Explicit Calculation of Stein Kernels
This appendix contains explicit calculations for the Langevin-Stein and KGM-Stein kernels kp- which are sufficient to implement Stein II-mmporaance Sampling and Stein II-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, observe that the diffusion Stein operator Sp in (12) applied to matrix-valued kemel K is equivalent to the Langevin-Ster operator applied to the kernel C(x., y) = .M) [xx y) (y) In the case of the Langevin- Stein and KGM-Stein kernels have K (zr, y) s(x, y) for some
<(x,, v) and M (1 || x. where
||
(-- (--1)/2 for some
{0, 1,2..... Thus C(x,, v) c(x,
cxx,))== = (1 + | x. (- + 11l xx113)( 1) 1),2))((((
and
kp(x,y) = V.. Vyc(x,)) v) + [V (V, logp(y)) + [Vye(a,))) [V logp(x))
+ c(x, v) V log)(z)) IV, logp()),
20

---

following the calculations in Oates et al. (2017). To evaluate the terms in this formula we start by differentiating c(z, 3). obtain
Voc(E, - (1+ 1
Ix
13/2 11
88 [){((1 y) y) + V V2K(Z2 v)
x
1 | - x.13
vol(,,)) (1+ + || xz x111 (0- 1)/2 (1 + 1 - x. 113)
(s 1))( ))((,,))" 9)2 (g
2.)
[, v\
x
1 |1y - 1.1
V. V Voc(o,u) = (1 + || xz (13) (5- 11 + |18 x4 16)
(s 1)2x( 4) x.) _ 2 2.n)
(s 1)(y -- (xn(x,))) (1 + 11 x 111)
x I, [1 + 11y I.
(s 1) xx) _ VgN(1, y)
+ Vx Vy
(1 + - x.
These expressions involve gradients of n(x, y). and explicit formulae for these are presented for the choice of s(x,, y) corresponding the Langevin-Stein kernel in Appendix C.I, and the KGM-Stein kernel in Appendix C.2.
To implement Stein II-Thinning we require access to both kp(x) and Vkp(2)) 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 k(z, y) is continuously differentiable, that partial derivatives with respect to x and y can be interchanged. Then
cx(x): c(I, x))
(1 + || -
"x(xx x)
(x) = V.c(1, y) --م
(s 1)x(I, 1))--1(a Iaa
= (1 + || - xi 0-1
+ 7, 10yy
(1 + || - xz 1
c2(x) = Vx Vyc(x, 8)1 -
(8 - 1)²K(H, x) - -xx)) T xx))
(1 + || xz 113)
(1 |x xx 11:)2
2(s - 1)(x - xx,) x,)TTTT T5-1 Va((E,)) v) + Vx Vyn(x. U)l
+
(1+ || - x.
so that
kp(x) = kp(x,)) = c2(x) + (x) Va logp(z) +
(19)
Let (z)c()) = (r](c((z)] and V log O2 0z logp(x). Now we can differentiate
(19) to get
(20)
In what follows we also derive explicit formulae for co(x), (x) and (2))), and hence for (
V (x) and (.)) for the case of the Langevin kernel in Appendix C.I, and the KGM- Stein kernel in Appendix C.2.
21

---

C.1 Explicit Formulae for the Langevin-Stein Kernel
The Sogevvn--teein kemee from Appendix A.3.1 corresponds to the choice s = 1 and <(22, y) the inverse I multi-quadric kernel, so that
x(x,y) = (1 + | - y 8
V xx((x,y) = -23(1 + || - y)
(( 28(1+ v -- y)
V V g)((,)) -43(3 1) + || - (z y) - 2x v)
+ +2 + || - - $-1
Evaluating on the diagonal:
<(z, ((2,x)=1 =
V Vak(2,y)lynp Vyk(E,8)) y-o 0
Vx Vyt((,; y) y-vz = 2,Str(E-1)
so that co(x)) = 1, (I) = 0, c2(x) = 20t Differentiaaing these formulae, Vzco(z) = 0. (zci())= 0. [x0((x) = O.
C.2 Explicit Formulae for the KGM-Stein Kernel
The KGM kernel of order s from Appendix A.3.2 corresponds to the choice
1+ (a - x x.)
<(x, v) = (1 + || - v -
xx xz 2/22) + 11 x 1/)
for which we have
Ver(E, v) = -23(1 + || - ull3) 88 (013))8-11-1 (x y)
Σ- - v.)) s + (r x.) T Σ (y- 'x xI,(11 + 1-= x, 13:)
(1 |/ x, 12 / + 11 x, 11 1//2
V ())((,))== = -33(1 + || - --1 v|13) -111))--11--6
y)
(x - xx) - + (x - x. 110-+----011010--0))00000-----
11y - x. =1133--1
(1+ || x.
V.. VoH(z,8)= = -43(3 + 1)(1 || - --- + )1 + ||0 -113)
tr( - ---
X.)
- s(1 + p z.
(y 2.) 22.,)))222 (1-2,))>-'(-11)
( + (x - x.)
Σ (y - x2) + || +11++11 x. (1 11y x.
x(x - I.) E y x2)
(1+ | - 3/ ly z. ) 1/2
Evaluating on the diagonal:
n(x,x) 1 + (1 + x.)) -++1
Von(z,y)ly- y VYN(I, -(s () x, (1 + ||ا - 11,133 -1
V. Vyn(Ie, اy| = -atr ( |1=
+ - (0(s-2)(11 + ||I - 2, 12,3)) (x-x)) T
Σ - - xT.)
so that
co(x)) = 1 + (1 + || - x, 21112)
(x) = (s - 1)(1 + || - xx 112)) 25 (x - x..)
c2(x) = [s (ss- + || x. )-1 1][x x.)
( (x I.)
+ ir(E + 23 + || - x. 11
(1 [1+ x Iz 22
(1+ x x. 12
22

---

Differentiating these formulae:
Vaco(x) 2(s - 1) |a - (x [c
(x)((x)= = - cx
b 3 ( xx,]
+ (1 - دー
Vxc((x) - - 10) 2+) T 2 x., 2-)1-=1 (x x.)
-
+ - =2
x.,)
- -2 + | 1x=- 2- (x x+a)
+
+r(22 112 (x x.)]
+ +4(1 + || - -2) x. -x.,) -r.).
These complete the analytic calculations necessary to compute the Stein kernel kp and its gradient.
D Empirical Assessment
This appendix contains full details of the empirical protocols that were employed and the additional empirical results deseribed in the main text. Appendix D.1 discusses the effect of dimension our proposed II. Additional illuatrative results from Section 3.2 are contained in Appendix D.2.
The full details for how MALA was implemented contained in Appendix D.3. An additional illustration using a generalised rugo-regressvve moving average (GARCH) model is presented in Appendix D.4. The full results for the convergence of the sparse approximation provided by SITT-MALA to the optimal weighted approximation is investigated. Finally, the performance of KSDD is quantified using the I-Wasserstein divergence in Appendix D.7.
are contained in Appendix D.5, and in Appendix D.6
D.1 The Effect of Dimension on II
The improvement of Stein II-Importance Sampling over the default Stein importance sampling algorithm (i.e. II = P) can be expected to reduce the dimension d of the target P is increased. To see this, consider the Langevin-Stein kernel
kp(x) + log])oo)((>>)]
(21)
for some C1, C2 > 0; see Appendix C. Taking P = NO, laxa), for which the length scale matrix Σ appearing in Appendix A.3 Laxd, we obtain
kp(x) ci + c2||2/12.
However, the sampling distribution IT defined in (8) depends on kp only up unspecified normalisation constant; we may therefore equally consider the asymptotic behaviour of kp(x) := kp(x)/d. Let X ~ P. Then (X)) eg is d-independent constant, and
(k(()) - (ci + =odd]² kpp Elkp (P(()) = [kp(ee)
203 dP(x) = d
- 0
d co. This shows that kp converges to constant function in I,2(P), and thus for "typical" values of x in the effective support of P.
t Txx) x p(a)) kp(x) & p(x),
so that R P in the d - oc limit. This intuition is borne out in simulations involving both the Langevin-Stenn kernel (as just discussed) and also the KGM3-Stein kernel. Indeed, Figure SI shows that as the dimension is increased, the marginal distributions of become increasingly similar to those of P.
23

---

Langevin Stein Kernel
KGM3 Stein Kernel
0.00
0.30
14
0.31
18.55
1020
0.00
0000
0/0
(SI.1)
(S1.2)
Figure S1: The effect of dimension on II: Here P was taken to be the standard Gaussian distribution (0, laxd) in Rd and the proposed distribution II was computed. The marginal distribution of the first component II is plotted for d E 21,2, 10}, for both (a) the Langevin- Steiin kernel and (b) the KGM3 Stein kernel.
10
B[KSD]
P (Langovin) (Langevin)
P (KGM3)
11 (KGM5)
P (Riemanns)
10
I (Rieraamn)
Figure S2: Assessing the performance the sampling distributions IT shown in Figure The mean kernel Stein discrepancy (KSD) for computation performed using the Langevin-Ster 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 construct II. Solid lines indicate the baseline case of sampling from P. while dashed lines indicate the proposed approach of sampling from II. (The experiment was repeated 10 times and standard error bars are plotted.)
D.2 2D Illustration from the Main Text
Section 3.2 of the main text contained a 2-dimensional illustration of Stein II-Importance Sampling and presented the distributions II corresponding to different choices of Stein kernel. Here, in Figure S2, we present the mean KSDs for Stein Il-Importance Sampling performed using the Langevin-Steii Steii kernel (purple), the KGM3-Stein kernel (bluc), and the Riemann- kernel (red), corresponding to the sampling distributions displayed in Figure 2 of the main text.
For this experiment, exact sampling from both P and II was performed using fine grid on which all probabilities were calculated and appropriately normalised. Results are broad agreement with the I-dimensional illustration contained in the main text, in the sense that in all cases Stein I1-Importance Sampling provides a significant improvement over the default Stein importance sampling method with II equal to P.
24

---

Algorithm 4 Adaptive MALA
Require: 20.0 (initial state), en (initial step size), Mu (initial preconditioner matrix), (niheo (epoch lengths), {as fi=1 (learning schedule), h (number of epochs), kp (Stcin kemel)
[#0.1 150,1-....2,100 m0,mo) MALA( 6o, Mo, no, kp)
2: for i 1....,h do
3:
< -1-1.n;-1
4: Pi-1 + 1
n;-1 EØ=1 1 -
b Average acceptance rate for chain i
5:
6 f 6:-1 exp(pi-1 0.57
> Update step size
6: M, a,M, +1 (a.)ev((fx;))) (.-1--1))))
Update preconditioner matrix
7: (2i,1 (2i.1..., Rinn; t MALA( 6s, Mi, Mi, kr)
8: end for
D.3 Implementation of MALA
For implementation of MALA in Algorithm 4 are required specify step size e and preconditioner matrix M. In general, suitable values for both of these parameters will be problem- dependent. Standard practice to perform some form of manual or automated tuning arrive at parameter values for which the average acceptance rate is close to 0.57, motivated by the asymptotic analysis of Roberts and Rosenthal (1998). Adaptive MCMC algorithms, which seek optimise the parameters of MCMC algorithms such as MALA during the warm-up period, provide appealing solution, and was the approach taken in this work.
The adaptive MALA algorithm which we used is contained in Algorithm 4, where we have let MALA(x, e, M, n, kp) denote the output from the preconditioned MALA with initial state æ, step size preconditioner matrix M, and chain length n, described in Algorithm 1. In Algorithm 4, we use cov(-) 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 6g = 1, Mo la. h = 10, and onı = lengths were nin (Th-1,11..... 1,,; } from the final epoch are returned, and constituted output from MALA for experimental assessment.
0g = 0.3. The warm-up epoch ng = 1, 000 and the final epoch length was ng = 10 The samples
To sample from instead of II, used Algorithm formally set kp(x) for all x E Rd. which recovers 11 = P as the target.
D.4 Illustration on a GARCH Model
This appendix contains additional illustrative experiment, concerning GARCH model that is particular instance of a model from the PosteriorDE database discussed Section 4. The purpose of this illustration is to facilitate empirical investigation in slightly higher dimension (d 4) and to explore the effect of changing the order s of the KGM-Stein kernel defined in Appendix A.3.2.
First we describe the GARCH model that was used. These models are widely- in econometrics to describe time series data (ye ) in settings where the volatility process is assumed be time-varying (but stationary). In particular, we consider the GARCH(1, model
y/ Φ1 + at
a, oie,, Et ~ (0, 1),
of dz + Φ307 0407-15
where d2 d3 > Φ4 > 0, and Φ3 d4 1 are the model parameters, constrained to subset of Rd. For ease of sampling, a change of variables r (Ø1, $2, Ø3. 64) + is performed in such way that the parameter E IR4 is unconstrained. Assuming an improper flat prior on f, the log-posterior density for 0 is given up an additive constant by
lo Vi.....,,)
+
Σ
log(o?)
vē 20
log J, (,,1((1)) (8)),
where (0) is the Jacobian determinant of -1
25

---

1.0%
3.0
II (RGM2)
2.5
II (KGM3)
0.8
(KGM4)
Density
2.0
0.6
1.5
0.4
1.0
0:2
0.5
0.0
0.0
4.8
5.0
0.2
6.4
5.6
-1.0
-0.5
0.0
0.5
2.0
01
02
0.7
0.35
0.5
0.30
0.5-
Density
0.4
0.20
0.15
0.10
0.0
0.00
-5
6
10
8g
04
Figure S3: Illustrating the shape of IT based on the KGMs-Stein kernel for a GARCH(1,1) model, controlling convergence of moments up to order {2, 4). The marginal density functions of each distribution were approximated using one-million samples obtained using MCMC.
For this illustration, real data were provided within the model description of PosteriorDB, for which the estimated maximum a posteriori parameter ò (5.04, 1.36, 0.53, 0.31). The marginal distributions of II corresponding to the KGM-Stein kernels of orders {2, 3, 4] are compared the marginals of P in Figure S3. It can be seen that higher orders correspond to greater over- dispersion of II; this makes intuitive sense since larger corresponds to more stringent KSD how the
(controlling the convergence of moments up to order s) which places greater emphasis tails of P are approximated. Further, for the final skewed marginal of P, we note that the distribution exaggerates the skew, placing more of its mass in the tail of the direction which positively skewed. Further discussion of skewed targets is contained in Appendix D.8.
D.5 Stein II-Importance Sampling for PosteriorDB
To introduce objectivity into assessment, exploited the PosteriorDh iorDB benchmark (Magnusson et al., 2022). This ongoing project is an attempt toward standardised benchmarking, consisting of collection of posteriors to be numerically approximated. The test problems PosteriorDD are defined in the Stan probabilistic programming language, and so BridgeStan (Roualdes et al., 2023) was used directly access posterior densities and their gradients as required. The ambition of PosteriorDB is to provide extensive set of benchmark tasks; at the time we conducted our research, PosteriorDB was at Version 0.4.0 and contained 149 models, of which 47 came equipped with gold-standard sample of size n 103, generated from long run of Hamiltonian Monte Carlo (the No-U-Turn sampler in Stan). Of these 47 models, a subset of 40 were found to compatible with BridgeStan, which was at Version 1.0.2 at the time this research was performed. The version of Stan that we used was Stanc3 Version 2.31. 1.0 (Unix). Thus we used a total of 40 test problems for our empirical assessment.
26

---

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. Table 1 in the main text summarises the performance of SITIS- full empirical results are provided.
fixing the number of samples to be n = 3 x 10 In this appendix,
For sampling from MALA, we used the adaptive algorithm described in Appendix D.3 with a final epoch of length Rmax = 10 Then, whenever set of n < Rmmax 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", This ensures that the performance of unprocessed MALA that we report is not negatively affected by burn-in, in, in so far as is practical to control.
Full results are presented in Figure S4. These results broadly support the interpretation that SITIS- MALA usually outperforms SIS-MALA, or otherwise both methods provide similar level of performance, for the sufficiently large sample sizes n considered. The sample size threshold at which SIIIS- MUIS-MALA outperforms SIS-MALA appears to be dimension--eependent A notable exception is panel 29 of Figure S4, a d = 10 dimensional task for which SIII-MALA provided a substantially worse approximation in KSD for the range of values of n considered.
Figure S4: Benchmarking on PosteriorDE Here we compared raw output from MALA (dotted lines) with the post-processed output provided by the default Stein importance sampling method of Liu and Lee (2017) -MAL; solid lines) and the proposed Stein II-Ipportanc Sampling method ISTIS-MAAA dashed lines). The Langevin (purple) and KGM3 Stein kernels (blue) were used for SIS-MALA and SITIS-MALE and the associated KSDs 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 SITSS
given in parentheses. -Seegevin-Stein kernel: MALA.
-SISMMAA, MALA, SIS-MALA, ----- -TIIS-MALA.
MALA. KGM3-Stein kernel:
earnings carn Le)ght()))
8P pois regr 5P regs 3LD)
E|JsSD]
(S4.1)
(S4.2)
kidiy Inleeonnnoombss33
kidg - koneorrooooiigs
(S4.3)
(S4.4)
27

---

mesquite lognesggite hggnorgghee//ggwhhsss1111 J)
arma a
ElJSSD
(S4.5)
(S4.6)
surinng loguen loghightt madel (D)
garch L
(S4.7)
(S4.8)
kidic kidsoore www.orrrewooooiiii1000 (D)
earnings lygam inl erod 46D
3F
(S4.9)
(S4.10)
kidiy Inlaron (on (DD)
Rndig wiiggwthh www.. kirbecore Jn Jnpe wwt.ortctiooooii. (GD))
[XDD
(S4.11)
(S4.12)
28

---

lidig. wetghg.ttttmmmm momm.wrk lndcoore internetion, )
kicig with work lidsecre jnner ((5))
ElJKSD
(S4.13)
(S4.14)
lódiq with momn work- kidsoore mom w.D)
low din. garses mix -oww dimn gaiss mix
154.15)
(64.16)
uesquiee bagneequ (5g
Lamm.xxample wmw..ssaccplelee0mm
E/KSD
3/
3F
1S4117)
(S4.18)
Mbre- Ne
Junn Mr(GD)
(XSD)
XSD)
(S4.19)
(S4.20)
29

---

mesquite lagmesquite (ggmequitcoonnd(((())
EEJKSDD
(S4.21)
(S4.22)
bhlll drive event hram drive 0(822)
bhall drive evect Jum drive 1(82)
3P
(S4.23)
(S4.24)
budeon lyxx hare- 30 Aot vollerta (8D)
aite logmsquite (8
E|KSD
E0KSDY
(S4.25)
(S4.26)
mosquite lo Jogass
mesyille moquite/sD
(S4.27)
(S4.28)
30

---

right..chode right schoois pesseredd (10D)
right schools right 1o
19
199
E|JKSD
BIJKDD
(S4.29)
(S4.30)
nes nes 102)
msl977 sesl L0D)
(S4.3))
(S4.32)
100 1960 ned LGD)
I1s198 1 Dedd 10D)
EQKSD
(S4.33)
(S4.34)
1057588
(S4.35)
(S4.36)
31

---

1906 (00)
besl LBD)
ElJKSDD
(S4.37)
(S4.38)
diamonns dimmsodds(2D)
meycle gp op(
E(SSD)
E|3DD
(S4.39)
(S4.40)
D.6 Stein II-Thinning for PosteriorDB
The results presented in the main text concerned limit at which the optimal weights w* can be computed in few seconds values n, sparse approximation methods are likely to required. In the main text we presented Stein II-Thinning, which employs greedy optimisation perspective to obtain a sparse approximation to the optimal weights at cost O(m²n), where m are the number of greedy iterations performed. Explicit and verifiable conditions for the strong consisteney of the resulting SIIT-MALE algorithm were established Section 3.3. The purpose of this appendix is to empirically explore the convergence of SIIT- using the PosteriorDB test bed.
x 10 samples from MALA, which near the
laptop PC. For larger
In the experiments we report the number of MALA samples was fixed 10 The results, in Figure S5, indicate that
10 and the number
of greedy iterations was varied from m 1 to for most models PosteriorDB the minimum value of KSD approximately reached when m anywhere from TO to 7, representing modest but practically significant reduction in computational cost compared to STIÍS-MALA This agrees with the qualitative findings reported in the original Stein thinning paper of Riabiz et al. (2022).
32

---

Figure S5: Benchmarking on PosteriorDB. Here we investigate the convergence of the sparse approximation provided by the proposed Stein II-Thinning method (SIIT-MALA), The Langevin (purple) and KGM3- kernels (blue) were used for SITT-MALA and the associated KSDs 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. [LaagevvinStein SIIT-MALA, KGM3-Stein kernel:
kernel:
SIIT-MALA.]
earnings car (3))
pois reg regs (DD)
B|KSD]
30
12'
10
7
(S5.1)
(S5.2)
kidlg IndaoonnneenhsaS.
ledig - -omloommmmmml(33)
m
m
(S5.3)
(S5.4)
mesquile lognogg Be Jogvolune (3D)
40w.0100111100 (44D)
XSDD
BXSDD]
m
(S5.5)
(S5.6)
33

---

carnings lognarn lognen..nggiightts.ooooo hngisighttmmale4))
garch garchll(4D)
ElJKSDD
EIJSD
(S5.7)
(S5.8)
kidic kidscore mocnbsiq (D)
earnings logearn interaction a(5))
7h
(S5.9)
(S5.10)
kiliy kidsoure juter linn/5D)
kidag with moo work kidbe initer and a
E|KSD
107
Y12
m
(S5.11)
(S5.12)
Lodgg. w.t....whhhmmom. work lo Becero, (nlen
with www.hh.hh.ooooooo.oooooooo kirbecore Joper weterntioorroo... 2((D)
XSDD
10
m
m
(35.13)
(S5.14)
34

---

iádiq. with mom.,wrkk work lidsoore momm mook_wrkk(D))
low dim gomen puie kow dim www.ddiignsaaa gamss mix(SD)
39
BlJKSDD
(55.15)
(65.16)
mescuite logmesquite (0g
kamm example Ima eoamaple(GD)
TM
7h
(S5.17)
(S5.18)
able- - K.
sbli Mr16D)
EXKSD
101
(S5.19)
(S5.20)
meseuite log.oe (ppneli(7)))))
[SDD
(S5.21)
(S5.22)
35

---

bhall wolllddiveeeeeeeeee poent hmm drioe ((SD)
bball drive pveat hmm _
100
ElJiSD
ElJKSD
(S5.23)
(S5.24)
budson lyex hare- Jocdon volterra| (.D)
maseaate lognns pit (itee
7h
(S5.25)
(S5.26)
meqquie lognzequilee Jww.ee/..//1100
nesqule mesqiite
E|KSD
101
(S5.27)
(S5.28)
<ighi right www.rrree (10D)
right ye hoods agglt (1G)) ()
KSD)
(S5.29)
(S5.30)
36

---

nos 1es 1es(10
mest97 nes ()))
ElJKSD
BlJKsD
10
(S5.31)
(S5.32)
nes ne6 nes nes L02)
mes1911 sesl L0D)
E|J(DD
TM
7h
(S5.33)
(S5.34)
100 1961968- ned LGD)
1H619022 Ded LTD)
E/XSD
E|KSD
101
(S5.35)
(S5.36)
1001506 me (C))
1000 Dow (U))
XXSD
m
(S5.37)
(S5.38)
37

---

dianoonds (2
meyclego -ood__p(0)))
E|JSDD
(S5.39)
(S5.40)
D.7 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 MAI in terms of 1-Wasserstein divergence reported. This was made possible since PosteriorDE supplies set of posterior samples obtained from long run of Hamiltonian Monte Carlo (the No-U-Tume sampler in Stan) which we treat as gold standard.
Full results are presented in Figure S6 and Table 2. Broadly speaking, in most cases the minimisation of KSD seems to be associated with minimisation of 1- Wasserstein distance, and in particular a significant improvement of SIIIS MALA over -IS-MALA is reported for z= 63% of tasks in the PosteriorDB benchmark. However there are some scenarios for which minimisation of KSD is loosely, if at all, related to minimisation of 1 1-Wasserstein divergence. In these cases, we attribute this performance to combination of two factors: First, the blindness to mixing proportions phenomena, described in Wenliang and Kanagawa (2021); Koehler et al. (2022); Liu et al. (2023), which is pathology of KSDs in general. Second, the Langevin- kernel cannot be expected to control convergence in 1- 1-Wasserstein, since convergence in I-Wasserstein is equivalent to weak convergence plus convergence of the first moment. Focusing therefore on the KGM3-Stein kernel only, it is encouraging to note that SIIIS outperforms SIS-MALA on 83% of tasks in PosteriorDs in the 1- Wasserstein metric, as shown in Table 2. However, it is interesting to observe that MALA performed well in the 1- I-Wasserstein sense across the PosteriorDi orDB test bed.
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 KSDs, including potentially KSDs with better or more direct control over standard notions of convergence (such as 1- I-Wasserstein) that in the future may be developed.
38

---

Figure S6: Performance of Stein discrepancies PosteriorDB. Here we compared raw output from MALA (dotted lines) with the post-processed output provided by the default Stein importance solid lines) and the proposed Stein II- sampling method of Liu and Lee (2017) dashed lines). The Langevin (purple) and KGM3- Importance Sampling method (SIIIS-AAL
Stein kernels (blue) were used for SIS- SIS-MALA and SITIS- 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 SIS-MALA, MALL)
MALA. Langevin- kernel: SIS-MALA,
SIIIS- KGM3-Stein kernel:
carnings - hrigh
Ap-possregr (3
B[WASD
E[WASDI
3:33P
2x
3:33
33
(S6.1)
(S6.2)
kicig beieeree (30))
kidig kidsecur moniiq()))
loswwla
(S6.3)
(S6.4)
mesquite lopruesqunte www.asgiittooggsoonnnee.ooo Jognohnnee(3)))
armal ((47)
B/WASD
WAWAS
0%
3
10
(S6.5)
(S6.6)
39

---

garch g
earnings www.oooggghh.oll(((((..
1x35
1830
10
(S6.7)
(S6.8)
kidiq kidsoore mombsio((D)
ea carnings Ingrarm..nteraationn (nter ((.D)
1:30
BWWASD]
BWWAS
6%
1338
10
(S6.9)
(S6.10)
kidie_with
kiidg- htt.ssoor__nneeettoonoooo (5D))
las vavaa
(S6.11)
(S6.12)
40

---

kidig mom. kidsone.jutesattion Intecattion (
lidiq mom. wew.oorreeeteacticn..222222
BWASD
16
(3).13)
(S6.14)
wwthwwthhhhhht.tttommm kidexee_mom ceceeeeooonnwwtttteeet (5)
low din gamss juix kow dian gouss (.D)
6% 20-2
E[WASD
WAASD
2:35
1
1
2
10
(5).15)
(66.16)
mesquite lommeqqrite JemesquiteJogaa(5DD (opra(5DD)
hmmmpssamplo (nmm_ssample(6D)
RIWASD
2
10
(S6.17)
(S6.18)
41

---

ablre M)(
abári blr(GD)
B/WASD]
E[[AASD]
9%x
(S6.19)
(S6.20)
ark a1((DD)
meeqqutte lognessquitte Iogvashul 7D)
4
3
2:10
(S6.21)
(S6.22)
Jril harm (()
WAASS
3x
10
OF
(S6.23)
(S6.24)
42

---

huiison.kyre hare lotka.voberra(SD)
mesquiee logmesquiteSS))
6x
B[WASTI
E[WASD
3:
(S6.25)
(S6.26)
ineqquite lognesquite logeas (SD)
mesquiee mis(iit((S
E(WAADD)
BWWASD
1x
3x
(S6.27)
(S6.28)
sightt.cchols- wwwgt..chooosssoooooo
cigght.chooolss.cmeooo/ wwwteroe 192))
aiglit.adooos bggbrattoook.uommetteeeeed.
12833
(S6.29)
(S6.30)
43

---

1981972 neal (0D)
nes neal 102)
E[[ASSD]
12
E[WASD]
20 101
2
(S6.31)
(S6.32)
Bes1980 nest
1os19884 mpsi
E/WASD
4x10
4% 003
3:10
(S6.33)
(S6.34)
968
nes 10D)
1962 nesi 10/)
6x10" 107
4.000%
BAASD
WAASS
4%
2
(S6.35)
(S6.36)
44

---

D06 neal (0D)
nes20000 nesl 102)
E[WASSD)
E[WASD]
(S6.37)
(S6.38)
dianoods (2
mocyle..g acoe_gg((62))
18
EJWASD
B[WASD
60x11 18
184x
(S6.39)
(S6.40)
D.8 Investigation for Skewed Target
This final appendix contrasts the I-Wasserstein optimal sampling distribution 11; (c.f. Section 2.1), with the choice of 11 that recommended in (8). In particular, we focus the KGM3-Sterr kernel under a heavily skewed P. for which II, and II can be markedly different.
For this investigation a bivariate skew-normal target was constructed, where the density is given by xx)) = 40(x)) } (62) ) ( (--32)). with á and respectively denoting the density and distribution functions of standard Gaussian. The density of P. together with the marginal densities of TT, and II, are plotted in Figure S7. It can be seen that, while both I7 and are over-dispersed with respect to P. our recommended II assigns proportionally more mass to the tail that is positively skewed.
The performance of Stein II-Importance Sampling based on 1I, and II is compared in Figure S8. Though both choices lead to an improvement relative Stein importance sampling algorithm with II = P. the use of II leads to significant further reduction (on average) in KSD compared to I1 Based on our investigations, this finding seems general; the use of II does not realise the full potential of Stein II- II-Imporance sampling when the target is skewed.
45

---

Langevin
Stein
Task
MALA
MALA
SIS
MALA
SIISS
MALA
MALA
3IS-
MALA
SIIIS -
camiing-cam beight sp_ sp_pois, regrep regr kidiq-kidscone, mombs kidia-kidkooe momig
3
6420.0
7230.0
9440.0
6420.0
7280.0
7390,0
0.0779
0,0758
0.0724
0.0779
0.0665
0.0778
3
0.282
0.581
0,508
0.282
0,950
0.604
3
0.826
2.19
2.01
0.826
2.53
1.82
mesquite- lomquiite-lgmmsquitt logvalame arma armll
0.0311
4
0.0127
0.0132
0.0131
0.0127
0.0144
0.0133
logeeum Ingheight_ male
1.46
1.75
1.34
1.46
1,56
0.725
garch sarch-garchll
0.260
0.241
0.243
0.260
0.365
0.263
kidiq- siddsore_mommiig
1,86
2.39
1.72
1.86
2.51
1,54
carnings-logeee joteraction
0.0245
0.0240
0.0240
0.0245
0,0252
0.0251
kidig kidscore interaction
5
13.9
14.5
15.3
13.9
14.4
20.4
kidig, with mom, work bork-kidscore interrctiinntttnnn__.o
0.289
kidig with mom, work-kidseree, interaction c2
0.258
0.269
0.252
0.258
0.638
0.430
kidiy with
work kidscore, interaction kidiq with mom, work- kork-kisscore mom_ work
0.914
0.904
0.889
0.914
1,46
1.17
1,01
1.05
1.06
1.01
1.74
1.43
low dim, gouss low dim gowss_ mix
0.0215
0.0206
0.0207
0.0215
0.0214
0.0212
mesquite leesquire-logmsqqiiee logmesqiite logva
himm, cxample-hmme example
6
0.0708
0.102
0.105
0.0708
0.161
0.12
shirc bir
6
0.0613
0.0614
0.0581
0.0613
0.0613
0.0605
shiri-bir
6
0.0729
0.0733
0.0843
0.0729
0.0722
0.0926
arK-arK
7
0,0544
0.0533
0.0525
0.0544
0.0618
0,0557
mesquitt-loggeesquiie logvash
0.165
0.158
0.161
0.180
0.171
bball drive_vvent_____eeeeeee eventteeetaaieeeee (-hmm drive
8
0.203
0,183
0.186
0.203
0.242
0.216
bball, event I-hmm drive_
0.812
0.825
0.754
hudson, lynx, hare-lotka, voltern
8
0.113
0.112
0.111
0.113
0.135
020
mesquite logmesquite
0,212
0.204
0.209
0.212
0.214
0.212
mesquite-logmesqqee logvas
0.197
0.192
0.196
0,208
0.203
mesquite- mesquite
8
123.0
122.0
114.0
123.0
123.0
121.0
eight, setisssshoolsseihhe eighhtethhhhhhhhteeeeeeeeees schooll_centeeeeee
9.50
schh_schhols-cighe
nonooetteeedd
10
16.8
16.7
16.7
16.8
16.7
16.7
nes 1972-nes
10
0,193
0.189
0.172
193
0.198
0.172
nes1976-nes
10 0.190
0.189
0.179
0.190
0.202
0.177
nos 1980-92s
10 0,229
0.228
0.222
0.229
0.256
0.233
nex 1984-ns
0.184
mes 1988-188-.. nes
10 0.212
0.210
0,185
0.212
0.214
0.185
1992-nes
0.172
0.169
0.159
0.182
0.168
nes 1996 nes
10
0.191
0.187
0.173
0.191
0.199
0.179
nes2000-nez diamonds-dianodds
10
0.275
0.273
0.274
0.275
0.306
0.284
26
0.648
0.647
0.643
0.648
0.648
0.658
meycle. gp
66
15.5
15.5
15.6
15.5
15.5
15.8
Table 2: Benchmarking on PosteriorDE. Here we compared raw output from MALA with the post-processed output provided by the default Stein importance sampling method of Liu and Lee (2017) (SIS-Mc) (MML)) and the proposed Stein II-mportannc Sampling method (ELIS-MALA) Here d = dim( (P) and the number of MALA samples
3 10 The Langevin and KGM3- Stein kernels were used for SIS-MALA and STTIS- and the associated I-Wasserstein distances are reported. Ten replicates computed and statistically significant improvement is highlighted in bold.
0.8
0.7
0.7
A
(TWass.)
0.5
11 (KGM3)
0.5
Donsity
0.5
0.4-
0.4
0.3
0.2-
0.1
0.1
0.0
0.0
2
3
-5
-3
-2
.11
T2
Figure S7: Comparing the proposed distribution IT (KGM3; based on the KGM3-Stein kernel) to 11, Wasss.. the optimal choice for I-Wasserstein quantisation from Section 2.1) for bivariate skew-normal target (d = 2). The marginal density functions of each distribution were approximated using 10 samples from MCMC.
46

---

C (ass.))
DI (KGM3)
10
101
Figure S8: Comparing the performance of using the proposed distribution II (KGM3; based on the KGM3 Stein kemel) II, Wass., the optimal choice for W-Waseerstein quantisation from Section 2.1) for bivariate skew-normal target (d = 2). The mean kernel Stein discrepancy (KSD) for Stein II-Importance Sampling was estimated; 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 IT. (The experiment was repeated 10 times and standard error bars are plotted.)
47


Reviews:
[{"review_id": "ZYWR1PXT9E", "reviewer": "Reviewer_fuoL", "summary": "This work proposes a design of a probability density that is more over-dispersed than the target density, so that, somewhat surprisingly, the resulting MCMC samples, after optimally reweighted, can achieve lower KSD than MCMC samples from the true target density. Consistency of the two proposed algorithms, SPiIS-MALA and SPiT-MALA, is proved. These two algorithms are benchmarked on PosteriorDB dataset to demonstrate their superior performance over raw MALA and SIS-MALA (MALA plus optimal reweighting).\n", "questions": "- I would like to see a rigorous statement and proof of the heuristic argument from Sec 3.1. The sketched-out argument makes sense (other than one detail --- see below) so I'm wondering why it is not made into a complete proof. \n  * In the paragraph below (S2), it seems to assume that $E_{x \\sim \\Pi}[\\frac{dP}{d\\Pi}(x)(k(\\cdot, x) - \\mu_P)] = 0$. Why is this true? I think $E_{x \\sim \\Pi}[\\frac{dP}{d\\Pi}(x)k(\\cdot, x) - \\mu_P] = 0$ but not when $\\mu_P$ is multiplied by the importance weight. Of course, if $\\mu_P = 0$ then it does not matter.\n- A heuristic argument is given in D.1 for the Langevin-Stein kernel and $P$ is the standard $d$-dimensional Gaussian. Aside from this very simple $P$, is it true that $P \\approx \\Pi$ in general? Can we say anything theoretical about it? \n- How is the simplex-constrained minimization in Algorithm 2 implemented? What is the time complexity? This also seems like the computational bottleneck of Alg. 2 and the reason $n$ is ony a few thousands in all experiments.\n- For the experiment done in D.8, if the plot is for Wasserstein-1 distance, will $\\Pi(1Wass)$ result in better numbers than $\\Pi(KGM3)$?\n", "limitations": "The authors have adequately addressed the limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The paper tackles an interesting yet to my knowledge underexplored problem in sampling, which is how to design a density so that the resulting samples provide a good discrete support for the ensuing optimally reweighting step that finds a weight to minimize KSD.\n- The angle of the attack is, although not new (e.g. Graf and Luschgy 2007), quite surprising (i.e. the density needs to be over-dispersed). The motivational section on the optimal quantization of Wasserstein distance is well-put.\n- The consistency of the two proposed methods is proved which is nice.\n- Many kernels are studied (Langevin-Stein/KGM/Riemann-Langevin-Stein) and are used in producing the empirical results.\n- The writing of the paper is excellent and lots of intuition and toy examples are given to illustrate the points.\n", "weaknesses": "- The consistency proof of Theorem 1 seems like a rather straightforward application of Theorem 2 (Durmus and Moulines 2022) and Theorem 3(Riabiz et al. 2022). In particular, it seems to me the same proof should go through for quite generic $\\Pi$, not necessarily the one that takes the form in (8). Hence, it is not clear whether the design (8) is theoretically justified, other than the heuristic argument given in Sec. 3.1.\n- There is a considerable gap between the heuristic argument in Sec. 3.1 and the proposed algorithm, which is the weights used in Sec. 3.1 is not optimal but in the algorithms they are. The authors claimed the choice $dP/d\\Pi$ is \"near-optimal\" (without justification), but then noted that using $dP/d\\Pi$ will perform substantially worse than the optimal weight, which seems contradictory to the first claim. \n- The experimental results comparing SIS-MALA and SPiIS-MALA are somewhat mixed. For an end user, there is no provided criterion on whether they should use the proposed method or the baseline SIS-MALA. Moreover, I cannot find standard derivation of the reported numbers.\n  *  There seems to be missing experiments that benchmark the performance of SPiT-MALA compared to baseline (e.g. $\\Pi=P4). The only experiment I found for SPiT-MALA is in D.6 where the consistency is verified. \n- In many places it is hinted that $P$ will be close to $\\Pi$ as the dimension $d$ increases. This implies that the proposed method is only applicable to small dimensions and thus the application value is limited.  Further analysis on the relation to dimension could be helpful. Moreover, there is only one experiment for $d=66$ (last row in Table 1) that corroborate the point that the extent of improvement decreases when the dimension increases; more experiments could be used to strengthen this point.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed report.\n \n>  it is not clear whether the design (8) is theoretically justified, other than the heuristic argument given in Sec. 3.1\n \nOur proposed $\\Pi$ is not the only choice for which consistency can be established; consistent approximation is possible also for $\\Pi = P$, for example.  Rather, we heuristically motivate a specific choice of $\\Pi$ that is expected to out-perform alternatives, and then we verify that consistency occurs for this specific choice of $\\Pi$.\n \n> The authors claimed the choice $dP/d\\Pi$ is \"near-optimal\" (without justification), but then noted that using $dP/d\\Pi$ will perform substantially worse than the optimal weight\n \nYou are also correct that there is a performance gap between the weights that we analyse in Sec 3.1 and the Stein importance sampling weights; we will replace the phrase \"near optimal\" with a more nuanced explanation, which acknowledges the performance gap but explains that nevertheless the weights that we analyse in Sec 3.1 are expected to perform much better than simpler choices, such as uniform weights.\n \n> For an end user, there is no provided criterion on whether they should use the proposed method or the baseline SIS-MALA.\n \nWe see this as a fundamentally difficult question, akin to asking how to pick a MCMC method.  In the MCMC setting, it is typical to try one algorithm and, if it is not performing well, to try a different algorithm instead.  On the other hand, if we are talking in theoretical terms, then low-to-moderate dimensional posteriors that are not highly multi-modal are likely to be well-suited to S$\\Pi$IS-MALA, but otherwise plain MALA (rather than SIS-MALA) is likely to perform best.  We emphasise that this is due to the fundamental pathologies of KSD itself, rather than our algorithms to minimise it.  A discussion will be added to the manuscript.\n \n> I cannot find standard derivation of the reported numbers.\n \nThe full results, including standard error bars, are contained in supplemental Appendix D.5 - space considerations prevented us from including these in Table 1.  It can be verified that the error bars are all relatively small.\n \n> There seems to be missing experiments that benchmark the performance of SPiT-MALA compared to baseline (e.g. $\\Pi=P$).\n \nThank you for the opportunity to discuss this point - we are confident that S$\\Pi$T-MALA can outperform the baseline Stein Thinning algorithm with $\\Pi = P$, but this will occur only at large sample sizes $n$.  The reason for this is that Stein Thinning is a greedy algorithm which favours the inclusion of high probability samples in its initial phase - once the modes have been well described it will only then move on to sampling from the tail.  We wanted to demonstrate this phenomenon, but it appears that in many cases we need $n > 3,000$ to see the behaviour just described.  Due to the super-linear cost of S$\\Pi$T-MALA, we have so far not been able to deploy sufficient computational resources to thoroughly examine this effect.  This, together with the development of more efficient approximations to S$\\Pi$T-MALA, are actives area of ongoing work.  As such, we focused most the manuscript on S$\\Pi$IS-MALA, only noting in passing that we also obtain a consistency proof for S$\\Pi$T-MALA.\n \n> There is only one experiment for $d=66$ that corroborate the point that the extent of improvement decreases when the dimension increases; more experiments could be used to strengthen this point.\n \nThis is a good suggestion; at the time of writing this was the highest-dimensional model in PosteriorDB that complied, but we can explore the inclusion of other high-dimensional examples in the revised manuscript.\n \n> I would like to see a rigorous statement and proof of the heuristic argument from Sec 3.1\n \nWith respect, we believe everything here is rigorously stated and proven, albeit not within a theorem environment in latex.  \n \n> In the paragraph below (S2) it seems to assume that $\\mathbb{E}_{x \\sim \\Pi} [ \\frac{\\mathrm{d}P}{\\mathrm{d}\\Pi}(x) ( k(\\cdot,x) - \\mu_P(\\cdot) ) ] = 0$.  Why is this true?\n \nThe equality holds due to the following argument: \n\n$\\mathbb{E}_{x \\sim \\Pi} [ \\frac{\\mathrm{d}P}{\\mathrm{d}\\Pi}(x) ( k(\\cdot,x) - \\mu_P(\\cdot) ) ]$ \n\n$=\\mathbb{E}_{x \\sim P} [ k(\\cdot,x) - \\mu_P(\\cdot) ]$\n\n$=\\mathbb{E}_{x \\sim P} [ k(\\cdot,x) ] - \\mu_P(\\cdot) = \\mu_P(\\cdot) - \\mu_P(\\cdot) = 0$\n \n> Is it true that $P \\approx \\Pi$ in general?\n \nSince $\\frac{\\mathrm{d}\\Pi}{\\mathrm{d}P}(x) = k_P(x,x)^{1/2}$, the difference between $P$ and $\\Pi$ is driven by the Stein kernel $k_P$.  While $k_P(x,x)$ is usually an unbounded function as $\\|x\\| \\rightarrow \\infty$, the tail behaviour is controlled by the choice of Stein operator and base kernel.  There are of course moment-type constraints on $\\Pi$ for consistency of SPiIS-MALA than mean it cannot differ arbitrarily from $P$ for our theory to hold.  \n \n> How is the simplex-constrained minimization in Algorithm 2 implemented? What is the time complexity?\n \nThis is a linearly-constrained quadratic programme which we solved in Python 3.10.4 using the qpsolvers package version 3.4.0 as the frontend in conjunction with the ProxSuite package version 0.3.6 serving as the backend.  The full details are contained in the accompanying code, and we will add an explicit mention of the packages that can be used to run Algorithm 2 into the main text.  The time complexity is difficult to quantify, but we believe it is upper-bounded by $O(n^3)$.  These details will be included in the revised manuscript.\n \n> For the experiment done in D.8, if the plot is for Wasserstein-1 distance, will $\\Pi$ (1-Wass) result in better numbers than $\\Pi$ (KGM3)?\n \nWe can certainly investigate -- but Wasserstein-1 and KSD (KGM3) are quite different performance metrics, the former not capturing convergence of second and third moments, unlike the latter.  Our aim in this work was limited to designing an algorithm that is able to minimise a user-specified KSD.\n"}, {"review_id": "yflwkLGyZf", "reviewer": "Reviewer_X2c9", "summary": "The paper analyses which target distribution to use for MCMC in the situation where a stein-descrepency will be used to post-process it's output samples.\nThey propose to use a different target distribution for the MCMC to the distribution being approximated, and show this improves performance on a variety of posterior inference problems.", "questions": "How does using the stein-discrepency **during** sample generation (e.g. stein variational gradient descent) compare to using it for post-processing? \nI note that this question is not relevant for the contributions of the paper, but an answer would be useful for my understanding of the paper's general usefulness. ", "limitations": "The authors address the limitations of their method in the discussion, namely that it requires second order derivatives of the model", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- The paper identifies and clear question: how to choose the invariant distribution \\pi for MCMC if using a stein-descrepency to post-process the samples.\n- They provide a method for selecting $\\pi$ via a variational problem framing, where $\\pi$ is selected to minimize the variance in the post-processed approximation. This gives a closed form expression for $\\pi$ (up to a normalizing constant) such that it can be easily used within MCMC. Their analysis is agnostic to the choice of stein kernel making their results broadly applicable. \n- Figure 1 nicely illustrates the property of over-dispersion that their choice of $\\pi$ has, and shows that on a simple 1D problem that this results in lower error bars compared to using $P$ for the MCMC. \n- In their experiments their proposed method is shown to consistently improve (better results in 70% of  PosteriorDB tasks) upon the baseline of setting $\\pi$ to the target distribution being approximated. \n- Their presentation is generally very clear, with informative figures provided. ", "weaknesses": "I could not see any weaknesses in this paper, however the it's subfield is not within my area of expertise and I did not check the proofs. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your kind comments on our manuscript.\n \n> How does using the stein-discrepancy during sample generation (e.g. stein variational gradient descent) compare to using it for post-processing?\n \nThere have been some attempts to directly address this issue, in particular Stein Points (Chen et al, ICML 2018), Stein Point MCMC (Chen et al, ICML 2019), and Kernel Stein Discrepancy Descent (Korba et al, ICML 2021).  Whilst these algorithms do make use of Stein discrepancy for guiding sampling, it has to be acknowledged that these algorithms are not widely used.  SVGD is more widely used, but it is a gradient flow on the KL divergence rather than on the KSD.  At present it is unclear whether these algorithms will stand the test of time compared to MCMC, but we believe they are certainly worth investigating.\n"}, {"review_id": "yMemW29TKK", "reviewer": "Reviewer_yEHQ", "summary": "This paper proposes a proposal distribution $\\Pi$ to generate finite points \nsuch that a weighted version approximates the target distribution $P$ under \nStein discrepancy. ", "questions": "1. Can the authors comment more on the effect of dimensionality on the improvement of $\\Pi$ over $P$?\n\n2. Can the same proposal work for other kernel functions ( non Stein)? ", "limitations": "N/A", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The main strength lies in a new proposal for the sampling distribution $\\Pi$ that \nis more efficient for follow-up approximation to the original distribution $P$ in the \nsense of Stein discrepancy. \n\nAsymptotic consistency is established, and extensive simulations on Bayesian computation \nis also conducted to illustrate the benefit of the proposal.", "weaknesses": "The authors might consider expanding the section on the actual contribution (the proposal of $\\Pi$), \nand be more brief on the background. For instance, the section on Wasserstein distance doesn't seem necessary. \n\nThe second limitation is in the theoretical guarantees. Can one provide a non-asymptotic convergence guarantee \nfor the proposed method?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive comments on our manuscript.\n \n> the section on Wasserstein distance doesn't seem necessary\n \nThank you, we will think carefully about how to improve all aspects of our presentation, bearing in mind that Reviewers fuoL and pfRv appreciated this part of the manuscript specifically (\"the motivational section on the optimal quantization of Wasserstein distance is well-put\", \"a high level of clarity throughout, including an extensive review of related literature and methodologies\").\n \n> Can one provide a non-asymptotic convergence guarantee for the proposed method?\n \nWe believe this is possible in principle by exploiting Theorem 2 of Riabiz et al (2022).  However, this strategy would only provide a bound on $\\mathbb{E}[\\text{KSD}^2]$.  In contrast, our manuscript contains an almost sure (albeit asymptotic) convergence result.\n"}, {"review_id": "n3leBIuYy4", "reviewer": "Reviewer_m4Rk", "summary": "The paper studies the design of the MCMC algorithm which is well suited for post-processing to obtain consistent approximation $P_n^\\star$ of the target measure $P$ using Stein kernel discrepancies ($D_P( \\cdot)$). The authors suggest the following novel procedure: (1) choose a measure $\\Pi$ which differs from the target measure $P$ by the factor $\\sqrt{k_p}$, where $k_p$ is the Stein kernel  (solving variational problem) (2); Sample from $\\Pi$ using (pre-conditioned) MALA; (3) solve linearly-constrained quadratic problem or construct sparse approximation to obtain $P_n^\\star$.  \nTheorem 1 provides assumptions to ensure convergence of Stein kernel discrepancies $D(P_n^\\star)$ to zero.  Step (1) is based on SNIS procedure and the choice of $\\Pi$ which gives the smallest variance of the Stein kernel discrepancies between SNIS estimate P_n and target measure $P$. \nThe results are illustrated by numerical experiments.", "questions": "It is well known that MALA sampler is not good for mixtures of distributions. Is it possible to replace MALA by some other MCMC sampler as HMC or adaptive MCMC?\n", "limitations": "-", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "-Novel variational algorithm for explicit construction of measure \\Pi. \n-Theoretical analysis of the algorithm in the case of using MALA as a MCMC sampler\n", "weaknesses": "- Only asymptotic convergence in D_P metric is proved in Theorem 1. \n- Would be good to sketch ideas of construction of Stein's kernel in the main text. \n- It could be difficult to calculate Stein's kernel in more practical problems\n- Works in moderate dimensions\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for carefully considering our manuscript.\n \n> It is well known that MALA sampler is not good for mixtures of distributions. Is it possible to replace MALA by some other MCMC sampler as HMC or adaptive MCMC?\n \nThis can of course be done in practice, but ensuring consistency of the resulting algorithm could be considerably more difficult.  It is certainly a good suggestion, and one that we will take forward in further work.\n"}, {"review_id": "qgf85jkRY1", "reviewer": "Reviewer_pfRv", "summary": "This paper presents a novel approach for constructing an MCMC target that is specifically designed for post-processing using Stein importance sampling and Stein thinning, where the goal is to assign optimal weights to a subset of sample particles in order to construct the best possible approximations of a distribution $P$.\n\nThe proposed method introduces a new target distribution, denoted as $\\Pi$, which is obtained by tilting the original target density $p(x)$ with the square root of a Stein kernel $k_P(x)$. This construction is derived by solving a variational problem that minimises the trace of the variance of a limiting Gaussian distribution. The use of $\\Pi$ as a target for the Metropolis-Adjusted Langevin Algorithm (MALA) instead of the original target distribution $P$ is justified through an almost sure consistency guarantee (Theorem 1) and numerical experiments conducted on a set of benchmark problems.", "questions": "L174: Could you elaborate on what is the relaxation and what one loses due to this relaxation, if any?\n\nFigure 1.2: It seems that for small $n$, targeting $P$ is more beneficial than targeting $\\Pi$. Is it purely coincidental, or does it reflect a specific characteristic or trade-off associated with the methodology? Could you offer some insights into this observation?\n\nL229: Could you explain why the co-domain of $D_P$ is $\\mathcal{X} \\times \\mathcal{X}$ instead of $\\mathcal{P}(\\mathcal{X})$?\n\nTable 1: Are the results averages over multiple repetitions or obtained from a single repetition? t", "limitations": "Some limitations of Stein IS and the use of different Stein kernels are addressed in Section 4 and Section 3.2.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "**Originality**: Whilst previous works have extensively studied post-processing using Stein's discrepancy, this paper introduces a novel perspective on improving Stein importance sampling through the design of a target distribution that is distinct from the original target $P$ but more suitable for post-processing techniques. This novel approach of leveraging the target design to improve post-processing methods is both interesting and original.\n\n**Quality**: The construction of the proposed target distribution $\\Pi$ is clearly explained (Section 3.1). The paper presents both a theoretical guarantee (Theorem 1) and extensive numerical evidence to support the choice of $\\Pi$. The assumptions in Theorem 1 appear to be mild, and the authors provide comprehensive discussions comparing them with similar conditions in related literature. Overall, the significance of the proposed method is convincingly demonstrated.\n\n**Clarity**: This paper exhibits a high level of clarity throughout, including an extensive review of related literature and methodologies (Section 2).", "weaknesses": "**Motivation**: While the authors have effectively justified the use of the proposed target $\\Pi$ instead of the original target for constructing MALA samplers, the rationale for employing S$\\Pi$IS over running MALA **without** post-processing is slightly weak. Specifically:\n\n1. A main advantage of Stein importance sampling (SIS) lies in its ability to provide unbiased estimation when the MCMC sampler used to generate the sample is biased. Consequently, allocating computational resources to post-processing is crucial in such cases, as it is uncertain whether running the chain for a longer time would improve the quality due to the bias. However, for $\\Pi$ sampling, one must set up an unbiased sampler (e.g., MALA) that targets $\\Pi$, raising the question of whether conducting S$\\Pi$IS is more beneficial than simply allocating the same computational budget to running the chain for an extended period. The reported experimental results do not seem to address this question, as they solely compare S$\\Pi$IS with standard SIS. Including experiments comparing S$\\Pi$IS with MALA without post-processing would be valuable in addressing this concern.\n2. S$\\Pi$IS requires the setup of an MCMC sampler targeting a distribution different from $P$. Thus, the generated sample must be used in conjunction with the post-processing step to perform inference on $P$. This contrasts with standard SIS, where samples can be drawn from a sampler targeting $P$, which is typically what practitioners would do regardless of their intention to use post-processing methods. This raises the question of whether setting up a sampler that targets such a specialised distribution and relies on post-processing is more practically attractive than directly targeting $P$. Including discussions addressing these concerns would be beneficial.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful report.\n\n> A main advantage of Stein importance sampling (SIS) lies in its ability to provide unbiased estimation when the MCMC sampler used to generate the sample is biased.  However, for $\\Pi$ sampling, one must set up an unbiased sampler (e.g., MALA) that targets $\\Pi$\n\nThis is a good point; for example, the unadjusted Langevin algorithm (ULA) can generate approximate samples at lower computational cost compared to MALA, and Stein Importance Sampling can be used to retrospectively correct for the bias in ULA.  If it were clear how to design a ULA algorithm to target $\\Pi$ then this could be directly applied in Stein $\\Pi$ Importance Sampling, however in general the explicit characterisation of the invariant distribution for ULA is intractable, making it unclear how to proceed.  This would be an interesting direction for further research and we will add a discussion on this point to the conclusion section of the revised manuscript.\n \n> Including experiments comparing S$\\Pi$IS with MALA without post-processing would be valuable in addressing this concern.\n \nPlease allow us to point out that the columns labelled \"MALA\" in Table 1 correspond to classical $P$-invariant MALA, which we believe is what you have asked for?  It is true that, then $P$ and its gradients are cheap to evaluate, the computational cost of MALA is lower than that of S$\\Pi$IS-MALA, and one could run more iterations of MALA for an equivalent computational cost.  But for more complex $P$ the computational cost of all algorithms will be gated by the number of times $P$ and its gradients need to be evaluated, making the direct comparison in Table 1 meaningful.  Further, if we aim for a compressed representation of $P$, then some form of post-processing of MALA would be required, which would then entail an additional computational cost.  We aim to discuss these important points in more detail in the revised manuscript.\n \n> This raises the question of whether setting up a sampler that targets such a specialised distribution and relies on post-processing is more practically attractive than directly targeting $P$. Including discussions addressing these concerns would be beneficial.\n \nAs we also mentioned to Reviewer X2c9, at present it is unclear whether these algorithms will stand the test of time compared to MCMC, but we believe they are certainly worth investigating.  We will expand the conclusion section of the manuscript to address this broader point.\n \n> L174: Could you elaborate on what is the relaxation and what one loses due to this relaxation, if any?\n \nActually, nothing was lost due to relaxing the constraints (S1-2), since we verified in lines 177-178 that the solution to the relaxed problem also happens to satisfy the constraints (S1-2).\n \n> Figure 1.2: It seems that for small $n$, targeting $P$ is more beneficial than targeting $\\Pi$. Is it purely coincidental, or does it reflect a specific characteristic or trade-off associated with the methodology? Could you offer some insights into this observation?\n \nIn this particular example $P$ is uni-modal while $\\Pi$ is multi-modal, and we conjecture that there is a \"warm up\" period where samples are needed in each of the \"modes\" of $\\Pi$ before S$\\Pi$IS starts to perform well.  In contrast, SIS requires only samples from $P$, which is uni-modal, leading to a shorter \"warm up\" period.\n \n> L229: Could you explain why the co-domain of $D_P$ is $\\mathcal{X} \\times \\mathcal{X}$ instead of $\\mathcal{P}(\\mathcal{X})$?\n \nThank you for catching this typo!  Indeed, $D_P : \\mathcal{P}(\\mathcal{X}) \\rightarrow [0,\\infty]$.\n \n> Table 1: Are the results averages over multiple repetitions or obtained from a single repetition? \n \nThey are averages over ten replicates; in the caption we wrote that \"ten replicates were computed\", but we will make explicit that averages over the ten replicates are being reported.\n"}]

Author rebuttal:

Meta-review: This paper was well received by all the reviewers. The paper makes an elegant contribution the active field of using stein methods for MCMC processing. They consider the task of optimizing the target distribution of an MCMC where its known that stein-descrepency for target measure P will be used to post-process it's output samples. They show that surprisingly the optimal choice to target is not P, and provide a characterization of this choice. The paper is well written with numerous illustrations, solid theory, and a compelling experiment section, and thus a clear accept.
Decision: Accept (spotlight)