Kernel Stein Discrepancy thinning: a theoretical perspective of pathologies and a practical fix with regularization
Authors: Clement Benard
Brian Staber
Sébastien Da Veiga
Venue: NeurIPS 2023 poster (2023)
Keywords: Bayesian inference
Markov chain Monte Carlo
kernelized Stein discrepancy
Stein thinning
kernel methods
TL;DR: Stein thinning, an algorithm to post-treat MCMC outputs, often exhibits strong pathologies, that are carefully analyzed to propose an improved regularized algorithm.

Abstract: Stein thinning is a promising algorithm proposed by (Riabiz et al., 2022) for post-processing outputs of Markov chain Monte Carlo (MCMC). The main principle is to greedily minimize the kernelized Stein discrepancy (KSD), which only requires the gradient of the log-target distribution, and is thus well-suited for Bayesian inference. The main advantages of Stein thinning are the automatic remove of the burn-in period, the correction of the bias introduced by recent MCMC algorithms, and the asymptotic properties of convergence towards the target distribution. Nevertheless, Stein thinning suffers from several empirical pathologies, which may result in poor approximations, as observed in the literature. In this article, we conduct a theoretical analysis of these pathologies, to clearly identify the mechanisms at stake, and suggest improved strategies. Then, we introduce the regularized Stein thinning algorithm to alleviate the identified pathologies. Finally, theoretical guarantees and extensive experiments show the high efficiency of the proposed algorithm. An implementation of regularized Stein thinning as the kernax library in python and JAX is available at https://gitlab.com/drti/kernax.

Kernel Stein Discrepancy thinning: a theoretical perspective of pathologies and a practical fix with regularization
Clément Bénard"
Brian Staber' Sébastien Da Veiga*
Safran Tech, Digital Sciences Technologies, 78114 Magny-Les-Hameatl France
ENSAL, CREST, F-35000 Rennes, France
{clement.c benard, brian. staber) (safrangrou com]
sebastien, davveigadennsai.fff.
Abstract
Stein thinning is promising algorithm proposed by Riabiz et al. [202 for post- processing outputs of Markov chain Monte Carlo (MCMC). The main principle is to greedily minimize the kernelized Stein discrepancy (KSD), which only requires the gradient of the log-target distribution, and is thus well-suited for Bayesian inference. The main advantages of Stein thinning are the automatic remove of the burn-in period, the correction of the bias introduced by recent MCMC algorithms, and the asymptotic properties of convergence towards the target distribution. Nevertheless, Stein thinning suffers from several empirical pathologies, which may result observed in the literature. In this article, we conduct a theoretical analysis of these pathologies, to clearly identify the mechanisms at stake, and suggest improved strategies. Then, we introduce the regularized Stein thinning algorithm to alleviate the identified pathologies. Finally, theoretical guarantees and extensive experiments show the high efficiency of the proposed algorithm. An implementation of regularized Stein thinning as the kernax library in python and JAX available at hhtps://gilaabcom//tttkkkaaaa
poor approximations,
1 Introduction
Bayesian inference a powerful approach to solve statistical tasks, and is especially efficient to incorporate prior expert knowledge of the studied system, to provide uncertainties of the estimated quantities. Bayesian methods have thus demonstrated high empirical performance for a wide range of applications, in particular in the fields of physics and computational biology, to just name few. However, the Bayesian framework often leads to the evaluation of expectations with respect to posterior distribution, which is not tractable [Green et al., 2015), except in the specific case of conjugate prior distribution and likelihood, which hardly occurs in practice. To overcome this issue, Markov chain Monte Carlo (MCMC) is one of the most commonly used computational methods to estimate these integrals. Indeed, MCMC algorithms iteratively generate a sample, which follows the targeted posterior distribution, the Markov chain converges to its stationary state (Robert and Casella 1999 Brooks et al. 2011]. Consequently, the quality of the resulting estimates strongly depends on the convergence of the MCMC and how its output post-processed Standard post- processing procedures of MCMC outputs consist in removing the first iterations, called the burn-in period, and thinning the Markov chain with a constant frequency. Burn-in removal aims at reducing the bias introduced by the random initialization of the Markov chain. The R convergence diagnosis of Gelman et al. (1995) is, for instance, a well known method for determining the burn-in period. On the other hand, thinning the Markov chain allows for compressing the MCMC output and may also reduce the correlation between the iteratively selected points. More recently, promising kernel-based
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

procedures were proposed to automatically remove the burn-in period, compress the output, and reduce the asymptotic bias (South et al., 2022]. These approaches consist in minimizing a kemel- based discrepancy measure D(P, Qm) between the empirical distribution Qm of a subsample of the MCMC output of size m, and the target distribution P.. In this respect, minimization of the maximum mean discrepancy (MMD) was investigated by several authors, but these strategies require the full knowledge of the target distribution P. whose density is not tractable in non-conjugate Bayesian inference.
Based on the previous works of Chen al. 2018 and Chen et al. [209 Riabiz et al. 2022
propose to minimize the kernelized Stein discrepancy (KSD), to design an effciintt kernel based algorithm to thin MCMC outputs score-based discrepancy measure, i.e., it only requires the knowledge of the score function of the target P. which is readily available in our Bayesian framework. Importantly, Gorham and Mackey] [2017) showed that under suitable mild conditions, the KSD enjoys good convergence properties. More precisely, the KSD) provided that the sample size is large enough. Therefore, KSD thinning is highly promising tool for postppsoes and measuring the quality of MCMC outputs. This article thus focuses on the Stein thinning algorithm proposed by Riabiz et al. 20022]. which consists in selecting m points amongst the n iterations of the MCMC output, by greedily minimizing the KSD distance. Thanks to the convergence properties of the KSD, the empirical measure of the selected points weakly converges towards the posterior law P. However, on the practical side, several articles [Wenliang and Kanagawal 2020 Korba et al. 2021 have noticed empirical limitations of KSD- sampling algorithms, especially for multimodal target distributions. In fact, these limitations happen to be quite problematic, even in simple experiments, and have been slightly overlooked in the literature far, in our opinion. Therefore, this article first focuses on the analysis of KSD pathologies in Section 2. taking both an empirical and theoretical point of view. Then, we propose strategies mitigate the identified problems, and introduce the regularized Stein thinning in Section 3 We show the efficiency of our algorithm through both theoretical analysis and extensive experiments in Section 4 Notice that proofs and additional experiments are gathered in Appendices 1-7 in the Supplementary Material. In the remaining of this initial section, we mathematically formalize the KSD distance and the associated Stein thinning algorithm.
non non-tractable Bayesian setting, The KSD (Liu et al. 2016 is a
valid distance to detect samples drawn form the target distribution,
Kernelized Stein discrepancy. Kernelized Stein discrepancy was independently introduced by Chwialkowski et al. (2016), Liu et al. (2016), Gorham and Mackey (2017) as a promising tool for measuring dissimilaritie between two distributions IP and C on Rd with d 2 1, whenever IP admits a continuously differentiable density p. and the normalization constant of p is not tractable. Let k: Rdd Rd + R be a positive semi-definite kernel and let H(k) be the associated reproducing kernel Hilbert space (RKHS) with inner product (., ()())) and norm || [| Kernelized Stein discrepancy belongs to the family of maximum mean discrepancies (MMD) [Gretton et al., 2006 defined as
MMD Q) = sup
[[ - [[)]))
where X P. Q. If the kernel l is characteristic, then the MMD is distance between probability distributions. In practice, the MMD may not be computable as it involves mathematical expectations with respect to P. whose density is not tractable. To circumvent this issue, Gorham and Mackey (2015) proposed the Stein discrepancy which relies on Stein's method [Stein 1972]. It consists in defining an operator Tp that maps functions g RRd + Rd to real-valuee functions such that B[T,o(()) [X)] O. with X ~ P. for all in (k) = [g: Rd - Rd Ex |9: (())) < 1). The probability measure on Rd is assumed to admit a continuously differentiable Lebesgue density p E C (Rd), such that [[| [log [12] < co. The Stein discrepancy is then defined as SD(P, Q) = suppec(() [E( (T)g))))) [Z)]) where Z ~ Q. If the Stein operator T, is chosen as (x)x). V p(x) + (x))))), then Stein's 's discrepancy has
the Langevin operator (x) closed-form expression known as kernelized Stein discrepancy [Chwialkowski et al., 2016 ( [Lu et al. 2016), KSD² Q) = [l,,,Z, Z'), where Z Q, Z' Q, and hip denotes the Langevin Stein kernel defined from the score function sp(x) = log for x, E Rdd as
(x,xx')) x') =Vx, k x')) + (x), Vxrk(x,')) + (x,(x'), ) Vzk(x,x)) + (xp( sp(xx (x'))))))))))) (x)))(xxxx))
(2)
The main advantage of the KSD is that it only requires the knowledge of the score function, and does not involve any integration with respect to P. Gorham and Mackey [0017 also established

---

convergence guaranties when the kernel k is chosen as the inverse multi-quadraate (IMQ) kernel function k(x,, x = (c ||x -x11i)--- with c > 0, 3 (0. the positive definite matrix the identity matrix, and the density p is distantly dissipative as defined below. Log-concave distributions outside of a compact set are a typical example of such probability densities.
Definition 1.1 (Distant dissipativity Gorham and Mackey
(017)) The density p C² (RR) is
distantly dissipative if lim inf x(r) > 0, where x((r) inf 2 (1)())_s,
||x yll r
xx x
y-200
Stein thinning algorithm. Let P be a target probability measure that admits density p. and let {x² }=z c. R be MCMC output. The Stein thinning algorithm [Riabiz et al. 2022] selects m S n particles X Xxm by greedily minimizing the kernelized Stein discrepancy. Given t - < m particles Xmp...., X the t-th particle is defined
Tt E argmin kp (x;; X;) + (E()...))
kip (Xx, j=1
where the KSD of an empirical distribution has been used to simplify the objective function. The kernel function k is usually chosen as the IMQ kernel function, defined above. for both its good theoretical properties and empirical efficiency. Indeed, several articles [Chen et al., 2018 Riabiz et al., 2022] have led extensive experiments to show the better practical performance of the IMQ kernel over other choices. Also notice that the bandwidth parameter E is quite influential on the algorithm performance, but happens to be very difficult to tune, as highlighted by CChopin and Ducrocq 2021] Indeed, since the normalization constant of the target distribution is unknown, no additional metric available to assess the precise performance of the thinning procedure when varies. Furthermore, the sample quality output by Stein thinning varies in an erratic fashion with respect to f, making the design of heuristic procedures for the choice of notoriously difficult. Following the literature recommendations [Riabiz et al. 2022], we use the median heuristic to set in our experiments, and refer to Garreau et al. 2017 for an extensive analysis of this approach for kernel methods.
2 Analysis of KSD Pathologies
Although kernelized Stein discrepancy is a highly promising approach to thin MCMC outputs, several
empirical studies have highlighted that KSD-based algorithms may suffer from strong pathologies in simple experiments [Wenliang and Kanagawa 2020 Korba et al.] [02T Riabiz et al. 2022 Liu et al.
2023). The most established KSD pathology is that Stein thinning ignores the weights of distant modes of the target distribution, leading to the selection of samples of poor quality by Stein thinning.
This problem, called Pathology I throughout the article, is analyzed in Subsection 2.T Additionally, Korba et al. 2021) also notice that KSD thinning may result in samples concentrated in regions of low probability of p. As opposed Pathology I, the mechanism leading to this problematic hehavior is not well understood in the literature, to our best knowledge. Subsection 2.2 is thus dedicated to the theoretical characterization and illustration of Pathology II. Throughout the article, we illustrate KSD thinning using the running example of a Gaussian mixture, defined in Example below, where initial particles are directly sampled from p to better highlight pathologies. We will come back the thinning of MCMC outputs in detail in Section 4
Example 1. Let the density p be a Gaussian mixture model of two components, respectively centered in (-p, u, 0a-1) 1) and (u, Oa- 1. of weights and 1 - w, and of variance oaIa. The initial particles {x; In are drawn from p. The KSD thinning algorithm selects m < n points to approximate p.
2.1 Pathology I: mode proportion blindness
We first focus on Pathology I, which states that Stein thinning is blind to the relative weights of multi- ple distant modes of target distribution. Indeed. Weeliang and Kanagawa [220)] show that the score sp is insensitive to distant mode weights. Consequently. the KSD distance is unable to properly iden- tify samples with different weights than those of the target, in finite sample settings, as long as samples are accurately distributed within each mode. To be more specific, we illustrate this pathology with our Example of Gaussian mixture in dimension 2. We set μ = 3 and a = to enforce the two modes to be well separated, and take an unbalanced proportion = 0.2 for the left mode, and 1 - = 0.8 for the right mode. We generate n = 3000 observations and select en = 300 particles with Stein thinning. Clearly, the red selected sample displayed in Figure has wrong proportions, with about half of the
3

---

particles in each mode, instead of the expected 20 - 80%, reflected by the initial black particles sam- pled from p. More precisely, over 100 repetitions of the Stein thinning algorithm, we obtain an average proportion of 0.53 particles in the left mode, with standard deviation of 0.08 across the 100 runs.
Although Wenliang and Kanagawa [2020] clearly show that the KSD distance insensitive to the mode weights in the specific case of Gaussian mixtures, the mechanism leading to the selection of about half of the particles in each mode by Stein thinning, as in Example remains unexplained in the literature, to our best knowledge. Therefore, we conduct a theoretical analysis in the general case of any mixture distribution with two distant modes, stated Assumption 2.1 below. For the sake of clarity, we only study the case of a number of modes of two, without loss of generality. Importantly, notice that a finite sample drawn from distribution with distant modes, takes the form of clusters of particles around each mode, as illustrated in Figure [ Then, Stein thinning selects particles among these clusters to approximate p. and these particles define empirical law of a density 9 with compact support around each mode. [Wenliang and Kanagawa [2020) explain that the score Sp is especially insensitive to the mode weights in these compact areas around modes, which is the root cause of the generation of samples with wrong proportions, in Figure I Therefore, Assumption 2. below defines this observed setting, required have Pathology I occur, where density q has compact supports around each distant mode. Additionally, we also need to formalize Assumption 2.2 which tells that the distributions of the two modes of the mixture g have close KSD distance with respect to the target p. In particular, this assumption can be easily verified when both p and q have symmetric mode distributions, since the KSD distance is insensitive to the weights of
Stein Tiinning
5.0
2.5
0.0
2.5
5500
-5
0
5
20
Figure 1: Illustration of Pathology
with the Gaussian mixture of Ex-
ample [ (d = 2, μ 3,
1, w = 0.2. n = 3000, m = 300). Initial particles are in black, and the Stein thinning output is red.
Assumption 2.1 (Distant bimodal mixture distributions). Let p and q be two mixture distributions in R4, made of two modes centered in (-,,, 0a-1) and (u, 02-1), with μ > 0. The distribution of each mode of p E C² (TR)) has Rd as support, whereas each mode distribution of q have a compact support, included in ball of radius r > 0. with r < μ.. The left mode of p has weight "p * 1/2, and the right mode has weight 1 - Wp. Similarly, w and 1 - w are the mode weights of q. Let Q and Qu be the probability measures that respectively admit the density of the left and right modes of q. and P and Qu. be also the probability laws for p and q.
Assumption 2.2. For distant bimodal mixture distributions q and p satisfying Assumption 2.1 and for 17 E (0, 1), we have (PP Qt)/KSD² (P, Qx) - < 17.
Theorem 2.3. Let kp be the Stein kernel associated with the radial kernel k(x, x') = o [| x 1l2 where x, x' E Rd, E > O, and & E C²(R), such that 6(2) - O. 6'(2) + O, and dz(2) 2) - for 00. Let p and be two bimodal mixture distributions satisfying Assumptions 2.J and 2.2 for
any 17 E (0, 1). We define w" as the optimal mixture weight of q with respect to the KSD distance, i.e., -(-"
w= -argmin KSD(P, Q..). Then, for p large enough, we have |w=
=E(0.1]
Theorem 2.3 proved in Appendix B states that the weight wt of the optimal mixture g, which minimizes the KSD distance to the target p. is close to 1/2 regardless of the true target weight wpp. whenever the distributions of the two modes of the mixture have a close KSD distance to p. and provided that the two modes are distant. In particular, this is the case in the experiment of Example and Figure T where the two modes are symmetric and well separated. Additionally, a more specific empirical illustration of Theorem 2.3 can be found in Appendix[A.]] In Section 3 we will propose strategies improving Stein thinning to recover samples with accurate mode proportions.
2.2 Pathology II: spurious minimum
The core of this section is dedicated to the theoretical characterization of Pathology II. We first need to introduce additional notations to formalize our analysis. We thus define Ms the region of the input space where the score norm is lower than the threshold so 2 0, formally Ma, {x

---

Stein Thinning
5.0
20
2.5
0.0
1.0
I-p Io
-2.5
-5.0
0.
-5
5
()
2
Figure 2: Illustration of Pathology II for the Gaussian mixture of Example I (d 2. 1 p= = a= 1,
W = 0.5, = 3000, m 300): many particles are selected around the line x²1) = 0 (left panel), because of the squared first component of the score sp(x) along z²!) (for x(2) = 0 in the right panel).
R² [x)]( < $o). We also introduce an independent and identically distributed (iid) sample X Xm of P., with P,, Pm the associated empirical measure for positive integer m. and X the j-th component of X. Then, Theorem below shows that samples concentrated in regions of the input space where the norm of the score is low, have smaller KSD than samples drawn from the true target distribution p. for small sample sizes. Additionally, the score norm is low around stationary points of p. including local minimum and saddle points, shown in Corollary 2.5 below. However, samples concentrated at local minimum of p are bad approximations of the target distribution by definition. Therefore, pathological samples may be generated by Stein thinning. which minimizes the empirical KSD, and thus explains Pathology II observed by Korba et al. [222]] and shown in Figure 2 For the sake of clarity, we formalize our result for the IMQ kernel used in practice, and set c 1 without loss of generality, since it is equivalent to tune e or in the Stein thinning algorithm.
Theorem 2.4 (KSD spurious minimum). Let kp be the Stein kernel associated with the IMQ kernel with & > 0, 3 E (0, 1), and e = 1. Let {x: 1i
v C Mas {x E IR ||sp (x)(x)]l2 S 8 be fixed
- set of points of empirical measure Qm m Li=1 8(x:), with 8o 2 and m 2 2. We have KSD2 (P, Qm) < E[KSD2 (P. m.)) if the score threshold so and the sample size m are small - 32)/(23d/2 + 23ss// + sQ).
enough to satisfy m < + (E[llsp(X)
Corollary 2.5 (Low KSD samples at density minimumm. Let kp be the Stein kernel associated with the IMQ kernel with f > 0, 3 E (0, 1), and == 1. Let be density with at least one local minimum or saddle point. For m 2 2, if {x: m C RRd is a set of points, all located at locgl minimum or saddle points of p. then we have KSD P. Qm) < LEKSD* (m,Pm)], if m < 1 252 (1)pp
The proofs of Theorem 2.4 and Corollary 2.5 reported in Appendix C. are built on the idea that the KSD of the empirical law of (x; ) has bias of the form E (xx))l//m² Consequently, when m is small, the bias has a strong influence on KSD estimates, which favor samples concentrated in regions of low score norm, as stationary points of p. This mechanism is illustrated in Figure 2 and Corollary [266 for Gaussian mixtures. In this case, Stein thinning aligns a large number of particles around the line of saddle points defined by ())) 0. an area of low probability of the targeted mixture distribution, because of the variations of the score function, if the sample size m is small enough, From another perspective, for any sample size m. exists a Gaussian mixture with /1/a a large enough, such that Pathology II occurs. Therefore, Pathology II can appear for arbitrarily large samples m, depending on the target distribution properties.
Corollary 2.6 (KSD spurious minimum for Gaussian mixtures). Let kip be the Stein kernel associated with the IMQ kernel with & > O, 3 E (0, 1). and c = L. Let the density be a Gaussian mixture model of two components with equal weights, respectively centered in (-p, 0d-1) and (μ. 02-1). of variance o² I'Ia, and let v = u/o. If 1 and 0 S < [uv222 - In(v + vv² 1)) p,, then m C Ma. of empirical measure Qm, we have for any {x,}m
(i) KSD2 (IP, Qm) < E[KSD? Pm.)) if m and so satisfy < 1+ 28d;P4 42300/
E EXIS
(ii) there exists three disjoint intervals T-us To, Tp C TR, respectively centered around -u, 0, and p.
such that x (1) (1)
...., Zi
I-p U Ul U
5

---

3 Regularized Stein Thinning
Stein thinning suffers from two main pathologies, analyzed in Section 2 In a word, Pathology I comes from the insensitivity of the score to the relative weights of distant modes, whereas Pathology Il originates from the variations of the score norm, which do not differentiate local minimum from local maximum of the target distribution. We propose to regularize the KSD distance to fix these two problems, using terms that are highly sensitive the type of stationary point and the relative weights of modes. The proposed algorithm is first introduced in Subsection 3.1 then theoretical properties are discussed in Subsection 3.2 and finally the good empirical performance will be shown in Section [
3.1 Algorithm
Entropic regularization. In order to compensate the blindness of the KSD to mode proportions in multimodal distributions, we introduce the following entropic regularized KSD, denoted by KSDx, and defined as KSD (P, Q) = E k (2, Z') - and IP admits the density p. In our Bayesian setting, [[loo(p(Z)) is known up an additive constant since the normalization factor of is not tractable. However, it is possible to use KSD (P, Q) as the objective function of the Stein thinning algorithm, as the greedy selection of particles to optimize this quantity does not rely on the unknown additive constant. The main idea of this entropic regularization that - log log(((x)) takes higher values in modes of smaller probability, and therefore provides the relative mode weight information, which is missing in the KSD) distance. More precisely, modes with smaller weights take smaller density values, and are therefore more penalized than modes of higher weights. Therefore, with such entropic penalization, regularized Stein thinning tends to select particles in modes of higher weights more frequently than in modes of smaller weights, and we recover appropriate proportions.
- - AlEloog(p (p(g)((((()))))) [))] where Z and Z' have probability law Q
Laplacian correction. Chen et al. 2018 and Riabiz al.) [002 have noticed that the term kp(x,, Xx), which naturally appears in the empirical kernelized Stein discrepancy with the Langevin operator, can be interpreted as a regularization term. For example, Stein thinning does not select particles in the burn-in period of an MCMC output thanks to this regularization. However, this term k,(x, x2) also responsible for Pathology II, of samples concentrated in stationary points of p, as shown in Theorem 2.4 Therefore, we add second regularization term to compensate the weaknesses of (p(x,.,,) by penalizing particles located at local minimum and saddle points of the density p. Such points are located in areas of convexity of the target distribution, which can thus be detected with the positive values of the Laplacian of the density. Therefore, using the truncated Laplacian operator A [xx Ej=ı (x² (()))(r())) for function f E (R(R),, we propose the L.-KID) estimate with a Laplacian correction for densities E (IRd), defined by
L-KSD2 Qm)
(x,,
m² (x,, A log(p(x;))]
m
Regularized Stein thinning. Overall, we obtain the following estimate for the entropic regularized KSD with Laplacian correction L-KSD Qm. = L-KSD² Q... Eim log(p(x)).. Then, at each iteration & E {I,...,, .m), the regularized Stein thinning selects the particle index me e {1, [I....,n] to greedily minimize
m
(X Xa,) + A log(p(x..)) - At log(p((x,, +2 (xp Xx,)
Finally, Figure 3 illustrates the performance of regularized Stein thinning to fix the two pathologies analyzed in Section 2 in the case of Example [ with Gaussian mixtures. Indeed, the top panel of Figure [aaure/]shoe shows that the majority of particles are selected in the right mode, expected from the target distribution with w = 0.2. More precisely, an average proportion of 0.11 of the particles are located in the left mode over 100 repetitions of the procedure (0.89 in the right mode), with standard deviation of 0.03. For the value choice of X, we refer to the next subsection and the experimental Section 4 On the bottom panel of Figure 3 we observe that no particle is now selected on the line
a(1) = 0, as expected from the target Gaussian mixture distribution.
Remark 3.1. The truncated Laplacian operator is simply given by the trace of the Hessian matrix, where negative components are set to It follows that the computational cost the regularized algorithm is similar to the original Stein thinning.

---

The Laplacian correction of ky introduces recond-orde derivatives of p in the Stein discrepancy, and therefore enables differentiate local minimum and saddle points of the density p from its local maximum. A natural approach to introduce second-order derivatives of p in KSD) estimates, is to define the Stein discrepancy using second-order operators. Laplacian Stein operator [Oates et al. 2017] is derived in Appendix G but experiments show that this strategy is not efficient to fix Pathologies & II.
Remark 3.2.
Reguluized Steiu Thiming
Reguluiied Stein Thizning
5.0
5.0
2.5
2.5
0.0
0.00
2.5
-2.5
-500
Figure 3: Pathologies fixed by the regularized Stein thinning.
3.2 Theoretical properties
This subsection is dedicated to the theoretical analysis of regularized Stein thinning. First, we show that the proposed algorithm now enjoys good properties regarding Pathologies I and II, and thus mitigates the identified problems of the original Stein thinning. Secondly, we extend the convergence analysis of Riabiz al. [2022] for the post-treatment of MCMC output, show the convergence of the empirical law output by regularized Stein thinning towards the target probability measure.
Entropic regularization. In the previous section, Theorem 2.3 highlights how Pathology I of mode proportion blindness originates from the score insensitivity to mode weights. On the other hand, the entropic regularization is directly built on the target density, and therefore strongly depends on the mode weights. In the same setting of Assumption 2.1 required for Pathology to occur with the original algorithm, the following Theorem 3.3 shows that the entropic regularized KSD is minimized for the appropriate target weight, with the suitable regularization strength A. Notice that Theorem 3.3 proved in Appendix D is valid if [log [)] * [B[lo(p(ZZx)) with Z Qr and ZR- QR. otherwise the impact of the entropic regularization on wj vanishes. However, as wp _ 1/2 is required occur, p asymmetric, and [[log(p(22 (Z E[log(p(ZR)) (ZR)) is
in Assumption 2.1 for for Pathology only possible in very specific cases. Theorem 3.3 clearly shows that the regularized entropic KSD is sensitive to the weights of distant modes. Efficient strategies to choose the regularization strength will be first discussed in the asymptotic analysis below, and then in the experiments of the next section.
Theorem 3.3. Let kp be the Stein kernel associated with the radial kernel k(x, x') 6[]xx-x|12/8) where x, x' E Rdd E> O, and & C² (IR). Let p and q be two bimodal mixture distributions satisfying Assumption 2.7 We define Wx, the optimal mixture weight of q with respect to the entropic regularized KSD distance, i.e., wx = argmin KSDx Qw). If ([llo(p((Zz) * ([llo(p(ZR)) w0(
where Z ~ Qr and ZR= Qn it exists X E IR. such that w. = Wpp-
Laplacian correction. First, stress that the L-KSD is strongly consistent estimate of the KSD distance, where the proof follows from the law of large numbers. Therefore, the Laplacian correction introduced in the L-KSD estimate does not undermine the good asymptotic properties of the KSD distance. Secondly, the following theorem shows that samples concentrated in local minimum or saddle points of the target distribution and of low density values, are well identified by the L-KSDD as samples of worse quality than those truly sampled from the target. Consequently, the Laplacian correction fixes Pathology II, previously formalized in Theorem 2.4
Theorem 3.4. Let kip be the Stein kernel associated with the IMQ kernel with E. > 0, 3 E (0, 1), and 1. For m 2 2, let {x+)% C IR be a set of points located at Xo. local minimum saddle point of p. and of empirical measure Qm- Then, we have L--SSD2 (TP, Qm)) >> ELL.-DD* (TP P..)), if
the density at xo satisfies p(xo) < A Atp(x))/ (Elll X) [2 + E A log Ixx)).
7

---

Convergence of regularized Stein thinning, While regularized Stein thinning fixes finite sample size pathologies, the asymptotic properties of Stein thinning are also preserved. Indeed, if the initial set of particles is drawn from a different distribution than the target using a Markov chain Monte Carlo, Theorem 3t6 that the empirical measure of the sample obtained with regularized Stein Thinning, converges towards the target measure P. and thus extends the results of Riabiz et al. [2022). Notice that the weak convergence of a sequence of probability measure is denoted by =, and that distantly dissipative distributions are defined in Definition essentially states mild integrability conditions, and enforces that the MCMC output is not too far from sample drawn from p Appendix additional details.
I.I The required assumption below,
Assumption 3.5. Let Q be a probability distribution on R., such that absolutely continuous with respect to Q. Let {Z, HEEN C [Rd be Q-invariant, time-homogeeneous Markov chain, generated using a V Vuniformly ergodic transition kernel, such that V Vxx) 2 dE v 1/3d//2 + [xp(x)]] Suppose
that, for some y > 0, sup log(((2)))) iEN
< 00, sup
EE log log)(Z,) Z
< 00,
iEN
de sup m1x(11 2 aq (z.) ( +1ss, SEN
< CO, sup
E
dP
23d
+ (11s,(22)1BV(z.) < 00.
F2
EEN
Theorem 3.6. Let IP be distantly dissipative probability measure, that admits the density p E (IRIR), kp be the Stein kernel ussociated with the IMQ kernel where E, > 0,3 E (0, 1). Let (Zi)ieN C Ted be a Markov chain satisfying Assumption 3.5 π be the index sequence of length mn generated by regularized Stein thinning, and Qm, be the empirical measure of (Zx, J If log(n)) < mn < n, with any a > 1, and Amn. = =o(log(m) (m)) \/m)), then we have almost surely Q
P.
Theorem 3.6 proved in Appendix F provides us with interesting insights about the entropic reg- ularization strength X. We already know that À should be chosen with rate at least as fast as 0(1/m), to avoid the introduction of a higher a bias in the L--SDD than the original KSD. Indeed, for sample drawn from the true target distribution p. this bias E[L-KSD² (IP, Pm)]] takes the form (p()) Therefore, for slower rates of À than O(1/m). trivial samples concentrated a local maximum of the target p. can have smaller L-KSD than samples drawn from p. Then, Theorem 3.6 states that, for our ultimate application of MCMC post-processing, the Stein thinning sample distribution converges towards the target for such À rates of (1/mm) or faster. In practice, in our Bayesian setting, it is not possible to fine tune this parameter A because metric is available to assess the Stein thinning quality for various values of X, as already mentioned in the case of the bandwidth parameter f. In addition, we cannot theoretically determine which exact range of values of À leads to good thinned samples in a finite sample regime. However, we will see in the experiments of the following section that both slower and faster À rates than (11/m) lead to samples of degraded quality. Therefore, we set 1/m in the regularized Stein thinning, to ensure good empirical performance and the algorithm convergence.
(X.)) + E Log(p(X))]/m m
Empirical Assessment
This section shows how regularized Stein thinning outperforms the original algorithm through three batches of experiments: mixtures of standard distributions using exact or MCMC sampling, and Bayesian logistic regression real datasets. For the experiments considered in Sections 4.2 and 4.3 two Hatropolis-Hasttng samplers are considered with the Metropolis-Adiusee (MALA) and the No-U- sampler (NUTS). We use the IMQ kernel with l set with the median heuristic, 3 = 1/2, and c = 1, as recommended in Chen et al. (2018). J. Riabiz et al. (0022) We also set the regularization parameter with the default value of A = I m.. Notice that additional experiments are provided in Appendix A and that the code is available at
Langevin A Algorithm
When the target distribution is known, the efficiency of the Stein thinning algorithms are assessed by computing the MMD distance (see Equation (I) between a large sample drawn from the target distribution and the thinned samples. More specifically,
use the following closed-form expression of the MMD [rretton et al. 20066 with X, X' ~ P and Z, Z' ~ Q.
MMD (P, Q = [,,, X')] + E|k Blk, Z')] 21E[k(X, Z)]
(3)
where the kernel function & is chosen as the distance-induced kernel studied by Sejdinovic et al. (2013) and given by l(x, x') = ||x||2 + ||x'||2 - ||x - x2|12, for x, x' E Rd. In this setting, the MMD) reduces to the well known energy distance, as shown by SSedinovic et al. (2013.
8

---

4.1 Gaussian mixtures with exact sampling
As a first batch of experiments, we build on Example | and consider more complicated two- dimensional Gaussian mixtures to further illustrate the correction of Pathologies & II. The first Gaussian mixture made of four modes located at f4 = (-3,3) 22 (-3, 3), M3 = (3,3). and
M4 (3, -3), and with weights w1 W12 = 001,, and W3 = wi = 0.4, respectively. The second mixture is taken from [Quu and Wang 2023 It is made of 6 equally weighted Gaussian distributions centered at μ4 = (3 (2 - 1)/6), 3 sin (2r(i - 1) 1)/6)). for i ......6 6. For both experiments, we rely on exact Monte Carlo sampling to generate n 3000 observations, and apply Stein thinning and its regularized variant to select m 300 particles. The observed samples and the selected particles are shown in Figure The first example shows that vanilla Stein thinning does not capture the right proportions, while the regularized variant appropriately penalizes modes with lower weights. The second example illustrates Pathology II, which is corrected by the regularized Stein thinning.
Thimnin:
Rosiieiad
Thimnins
Requsaiianl Saiis Thiniig
25
2.5
0.0
B.
-24
22
-25
-25
-50
-.0
Figure 4: Gaussian mixtures with exact Monte Carlo sampling. Solutions (red dots) obtained by Stein thinning and its regularized variant.
4.2 Banana-shaped and Gaussian mixtures with MCMC sampling
We consider mixture of two distant modes of d-dimensional banana-shaped distributions with t-tails and unbalanced weights [Haario et al. 1999 Pompe et al. 2020). illustrated in Figure 5 and precisely defined in Appendix [ We sample this target banana mixture with both MALA and NUTS using three different step sizes E and iterations. The generated samples are post-processed with the Stein thinning and regularized Stein thinning algorithms, and their performances compared with the MMD between the post-processed samples and large samples drawn from the known target banana mixture. This experiment is for various thinning sizes m and dimensions d, with 20 repetitions to quantify uncertainties. The results obtained with the MALA sampler are shown in Figure the regularized Stein thinning clearly generates samples of higher quality than the vanilla Stein thimning. Additionally, an example of post-processed MALA output is depicted in Figure [] together with a heatmap of the Laplacian correction. On the left panel of Figure [ we see that pathologies are especially strong in this experiment, with large number of particles lying between the two modes in a region of low probability. On the right panels, we observe that regularized Stein thinning fix pathologies. Similar results were obtained with NUTS and are reported in Appendix A.2 for brevity. Next, we conduct the same experiments for d-dimensional Gaussian mixture of four modes with different variances, as detailed in Appendix A.2 Again, Figure shows the better
Stein Thinning
Regularized Stein Thininn
A log(p)
21
10-
0
-10
25
-20
Figure 5: t-banana-shaped mixture (d = 10). From left to right: solutions obtained with standard and regularized Stein thinning with contour lines of p, and heatmap of the Laplacian correction.

---

100
1.20
Thinin
Stso Thrccin
Regglarted Ston Tinrin
Haucarion Thinning
20
11:23
100
300
300
700
500
1100)
21
31
03
80
00
110
11.00
100
300
900
200
5001
1100
1100
21
31
10
51
50
73
80 50
00
m
Figure 6: (MALA) Graphs of the MMD with respect to the thinning size m (for d = 2) and with respect to (for m = 300). Left two panels: banana mixture. Right two panels: Gaussian mixture.
performance of regularized Stein thinning. Besides, we take advantage of this last experiment to explore other regularization rates than our default À = 1/m. Figure 12 in Appendix A.2 shows that slower rate of À 1/ log(m), which violates the convergence assumptions of Theorem 3.6, has significantly worse performance than the original Stein thinning. On the other hand, with a faster rate than 1/m such as 1/m² the effect of the entropic regularization disappears, and we recover similar results than the original Stein thinning. This supports that the default value of = 1/m is an efficient heuristic, since slower and faster rates of À strongly degrade the algorithm performance.
4.3 Bayesian logistic regression
We now compare the two Stein thinning algorithms in the Bayesian logistic regression setting for binary classification, since such problem usually involves multimodal posterior-see. e.g., Gersh-
man et al. [2012]. Liu and Wang 2016), Fong et al. 2019 Korba et al. [021] Given a dataset DN (X;, Ho made of N pairs of features X, I and labels Y, E fo. 1}, the probabil- ity that Y; is of class 1 is given by - 3, Bb) = 1/(11 + exp(- - 33" X,)), for some
parameters 0 (So, 3) E Rd+1 The prior distributions of the weight vector θ assumed to be and a Gamma prior with parameters (a, b) is chosen for the precision 00) Following [Fng al. 2019). the hyperparameters are chosen as a =
Gaussian, p( (2)) = =N(3G))
The posterior distribution of the weights is sampled with both MALA and NUTS using 48 independent chains, of respectively 104 and 10. iterations, and four step sizes & are considered along with three thinning sizes m. Each MCMC sample is post-processed with the two Stein thinning algorithms. For new input x*, the resulting thinned samples are used to approxi- mate the posterior predictive distribu- tion p(Y
Table 1: AUCs obtained with NUTS sampler for Stein Thinning (ST) and Regularized Stein Thinning (RST).
m
50
m
300
Dataset
ST
RST
ST
RST
Breast W. 0.88 (0.02) 0.96 (0.00) 0.93 (0.01) 0.96 (0.00)
Diabetes 0,52 (0.01) 0.50 (0.02) 0.53 (0.02) 0.57 (0.02)
Haberman 0.51 (0.04) 0.53 (0.02) 0.53 (0.03) 0.58 (0.02)
Liver 0.53 (0.04) 0.69 (0.01) 0.61 (0.04) 0.70 (0.01)
1|x", DN), defined by ())(O)DN) )dO. The per-
Sonar 0.80 (0.02) 0.81 0.001)0.81 (0.01) (0.01)
I J = 1|x*, 1|x*, formance of Stein thinning algorithms are assessed using the standard AUC metric for classification problems, estimated with 10-fold cross-validation and 10 repetitions for uncertainties. Table gathers the results for five public datasets from the UCI repository [Dua and Graff] 2017], and described in Appendix A..3 where the best AUC obtained for each algorithm over the four MCMC step sizes are reported. Clearly, regularized Stein thinning significantly improves the performance of Bayesian logistic regression.
5 Conclusion
Stein thinning has raised a high interest in recent years, as a powerful tool to post-process MCMC outputs, by the greedy minimization of the kernelized Stein discrepancy. Unfortunately, empirical studies have shown that KSD- algorithms suffer from strong pathologies. We have conducted an in-depth theoretical analysis to identify the mechanisms at stake. From this understanding, we propose an improved Stein thinning algorithm relying on entropic regularization and Laplacian correction. This approach exhibits relevant theoretical properties regarding pathologies, as well as highly improved empirical performance. Finally, the analysis of these regularization terms for other types of KSD-based algorithms, such a KSD descent, seems promising route for future work.
10

---

References
C. Bishop and N.M. Nasrabadi. Pattern Recognition and Machine Learning, volume Springer, 2006.
S. Brooks, A. Gelman, G. Jones, and X-I.. Meng. Handbook of Markov Chain Monte Carlo. CRC press, 2011.
W.Y. Chen, L. Mackey, J. Gorham, F. Briol, and C. Oates. Stein points. In International Conference on Machine Learning, pages 844-853, PMLR, 2018.
W.Y. Chen, A. Barp, F-X. Briol, J. Gorham, M. Girolami, L. Mackey, and C. Oates. Stein point markov chain monte carlo. In International Conference Machine Learning, pages 1011- 1021. PMLR, 2019.
N. Chopin and G. Ducrocq. Fast compression of MCMC output. Entropy, 23:1017, 2021.
K. Chwialkowski, H. Strathmann, and A. Gretton. A kernel test of goodness of fit. In International Conference on Machine Learning, pages 2606-2615. PMLR, 2016.
D. Dua and C. Graff. UCI machine learning repository, 2017. URL http::/ uci edu/ml
E. Fong, Lyddon, and C. Holmes. Scalable nonparametric sampling from multimodal posteriors with the posterior bootstrap. In International Conference on Machine Learning, pages 1952-1962. PMLR, 2019.
D. Garreau, W. Jitkrittum, and M. Kanagawa. Large sample analysis of the median heuristic. arXiv preprint arXiv: 1eXXvv:707..0m 1707.07267 2017.
A. Gelman, J.B. Carlin, H.S. Stem,, and D.B. Rubin. Bayesian Data Analysis. Hall/CRC, 1995.
Chapman and
S. Gershman, M. Hoffman, and D. Blei. Nonparametric variational inference. arXiv preprint arXiv:1206. 4665, 2012.
Gorham and Mackey, Measuring sample quality with stein's method. Advances in Neural Information Processing Systems, 28, 2015.
J. Gorham and L. Mackey. Measuring sample quality with kernels. In International Conference Machine Learning, pages 1292-1300 PMLR, 2017.
P.J. Green, K. Latuszyriski M. Pereyra, and C.P. Robert. Bayesian computation: summary of the current state, and samples backwards and forwards. Statistics and Computing, 25:835-866 862, 2015.
A. Gretton, K. Borgwardt, M. Rasch, B. Schölkopf, and A. Smola. kernel method for the two- sample-problem. Advances in Neural Information Processing Systems, 19, 2006.
H. Haario, F. Saksman, and J. Tamminen. Adaptive proposal distribution for random walk metropolis algorithm. Computational Statistics, 14:375-395, 1999.
A. Korba, P-C. International Conference on Machine Learning, pages 5719-5730. PMLR, 2021.
S. Majewski, and P. Ablin. Kernel stein discrepancy descent. In
Q Liu and D. Wang. Stein variational gradient descent: A general purpose bayesian inference algorithm. In Advances in neural information processing systems, pages 2378-237811 22-62 2016.
Q Liu, J. Lee, and M. Jordan. A kemeliizd stein discrepancy for goodness-of-fit tests. In International Conference on Machine Learning, pages 276-284. PMLR, 2016.
X. Liu, A.B. Duncan, and A. Gandy. Using perturbation to improve goodness-of-fii tests based on kernelized stein discrepancy, arXiv preprint arXiv:2304. 147622 2023.
C. Oates, Barp, and M. Girolami. Posterior integration on a riemannian manifold. arXiv preprint 1111v:1712.01193, 2017.

---

E. Pompe, C. Holmes, and K. Latuszyúísk. A framework for adaptive meme targeting multimodal distributions. The Annals of Statistics, 48:2930-29522 2020.
Yixuan Qiu and Xiao Wang. Efficient multimodal sampling via tempered distribution flow, Journal of the American Statistical Association, 000):--15, 2023. doi: 10. 10.1000006211599.0000 011144592222299805555 2198059.
M. Riabiz, W.Y. Chen, J. Cockayne, P. Swietach, S.A. Niederer, L. Mackey, and C.J. Oates. Optimal thinning MCMC output. Journal of the Royal Statistical Sociery: Series B, in press, 2022.
CP.. Robert and G. Casella. Monte Carlo Statistical Methods, volume 2. Springer, 1999.
D. Sejdinovic, B. Sriperumbudur, A. Gretton, and K. Fukumizu. Equivalence of distance-basee and rkhs-based statistics hypothesis testing. The Annals of Statistics, pages 2263-2291, 2013.
L.F. South, M. Riabiz, O. Teymur, and C.J. Oates. Postprocessing of meme. Annual Review of Statistics and Its Application, 9:529-555, 2022.
C. Stein. A bound for the error in the normal approximation to the distribution of sum of dependent random variables. In Proceedings of the sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 2: Probability Theory, volume 6, pages 583-603. University of California Press, 1972.
L.K. Wenliang and H. Kanagawa. Blindness of score- based methods to isolated components and mixing proportions. arXiv preprint arXiv:2008. 10087, 2020.
12

---

Appendix
A
Additional Experiments
A.1 Illustration of Theorem 2.3
To better illustrate Theorem we run an additional experiment, where is still defined as in Figure 7 from Example recalled below, with unbalanced mode weights of 0.2 and 0.8. The density is distributed as but each mode is truncated outside circle of two standard deviation radius, and y has weight w. Next, for various values of wr E 00 0.9), we draw two samples of size = 3000 from p and q. and compute (SD(P, Qw)) (with 30 repetitions for each value). The result is displayed Figure 8 and shows that the optimal weight is close to 1/2, as predicted by Theorem 2.3 since JKSD* Qt)/KSD² (P, QR) 1 estimated as 0.01 this case, implying that I 1 < 0.005.
Example 2. Let the density be Gaussian mixture model of two components, respectively centered (-1, 0g-1) and (μ, 0a-1 1), of weights w and 1 w. and of variance o The initial particles [x: in are drawn from The KSD thinning algorithm selects < n points to approximate p.
Stein Thinning
5.0
2.5
0.0
2.5
5.0
5
Figure 7: Illustration of Pathology with the Gaussian mixture of Example (d = 2, 1 3, 1, w 0.2,
3000, m = 300). Initial particles are in black, and the Stein thinning output is red.
A.2 Gaussian and Banana-shaped Mixtures
This appendix gathers additional results and details for the Gaussian and banana-shaped mixtures, as well as the MMD distance used to evaluate thinning performance, and the regularization parameter
Gaussian mixture. The second batch of experiments in Section considers d-dimensional Gaussian mixture of four modes of equal weight, with d 2 2, illustrated in Figure 9 The center of modes are chosen as (-2, 0). (2, 0). (-3,4) and (3, 4), and null values for higher dimension coordinates. The first two modes have an identity covariance matrix, while the remaining two modes have a diagonal covariance matrix with variance equal to 2. The results for regularized Stein thinning and the original Stein thinning are provided in Figure 10 for MALA sampler, and in Figure [I [ for NUTS sampler. both figures, the three tested step size are displayed, with a small impact on the resulting performance. Figures 9 10 and II show the high performance improvement of regularized Stein thinning over the original algorithm.
Regularization parameter X. Figure 12 displays the MMD obtained with regularization parame- ters X, set as À 1/m² and X = 1/ log m) These results should be compared with the ones shown
13

---

0.12
0.11
0.10
KSD Qu)
0.09
0.08
0.07
0.06
0.1
0.2
33
0.4
0.5
0.6
0.7
0.8
0.9
Figure 8: (S Qw) for as defined in Example with μ = 3. o = 1, wp 0.2, and q
truncation of and with weight w. The KSD estimated with 3000 and 30 repetitions for each w value.
Stein Thinning
Regularized Stein Thinning
A log(p)
7.5
3.0
2.5
0.0
-25
Figure 9: (Gaussian mixture with d = 2. MALA) First two panels: solutions obtained with Stein thinning and regularized Stein thinning with contour lines the target distribution. Last panel: heatmap of the Laplacian correction A log(p).
in Figures 10 and 11 which were obtained with regularization parameter À 1/m. These addi- tional experiments show the importance of choosing the regularization parameter as X = (11mm). suggested by Theorem 3.6 Indeed, slower rates of X give poor quality samples, and faster rates than À = O(1/m) tend to remove the effect of the entropic regularization, and we then recover similar performance than the original Stein thinning. On the other hand, À 1/m provides high improvement over the standard thinning, as shown in Figures 10 and 11
Banana-shaped mixture with t-tails. The first batch of experiments in Section considers a banana- -hanaa mixture with t-tails, defined as follows. Let şp : Fd Rd be the transformation defined by (z) = x; if * 2, and 42 (x) = 22 + bx² 1006. Let Z be random variable that follows the multivariate t-Student distribution with degrees of freedom Then, the random variable X (Z) + follows a t-banana-shaped distribution centered at μ.. We consider mixture of two s.banana-shapee distributions centered in O₂ and (0. 8, 0g-2), with weights w 0.25 and
0.75, respectively, which is illustrated in Figure 13 for d = 10. The results for regularized Stein thinning and the original Stein thinning are provided in Figure 14 for MALA sampler, and
in Figure 15 for NUTS sampler. In both figures, the three tested step size are displayed, which
confirms higher performance regularized Stein thinning. Figures [3 and T5 show the high performance improvement of regularized Stein thinning over the original algorithm.
14

---

0.1
Steun Thanann
0.2
Regularroed Statn Thiming
0.1
100
300
300
700
500
1100
1300
700
900
170
130
300
500
700
300
1100
306
G.5
MMDDP:Q.. P.Q..
MMN P. QQ.
MMDXP
Stein Thaning
Regularzed Stetn Thaming
33
40
10
81
33
100
10
100
20
10
20
100
Figure 10: (Gaussian mixture, MALA) Graphs of the MMD distance with respect to the thinning size m (with d 2) and with respect to d (with 300) for various step sizes
BLOS
0.2
00-
MMIEF..CCM (2.)
MMDDX
Stoin Thiming
0.2
Regglaariedd Thinming
2.0
MMDDE. Qc)
MSSIXP Qea
Stoin Thiming
MMDXE
1.0
Regularieed Thiming
8.6
0.0
Figure 11: (Gaussian mixture, NUTS) Graphs of the MMD distance with respect to the thinning size m (with 2) and with respect to (with m 300) for various step sizes e.
A.3 Bayesian Logistic Regression
This appendix gathers additional results for Bayesian logistic regression. In particular, Table provides a description of the tested datasets. Table gives the resulting AUC, for m = 50, 100, 300, using NUTS or MALA sampler. We recall that only the best AUC over the four tested MCMC step size e reported.
Recall that the Bayesian logistic regression defines the probability that Y, is of class 1 as p(Y, 1 1X,, (6) 1/(1 cxp(-90 - x,)), for some vector of parameters 0 (Bu. B3 E Rd+i The prior distributions of the weight vector O is assumed to be Gaussian, p(30)) ()) N 0, 1/), and a Gamma prior with parameters (a, b) chosen for precision 7.) Upon marginalizing, it is found that 366 is distributed as the non-standardized t-distribution Student-t(2a, 0, b/a) [Bishop and Nasrabadi 2006). Following (Fong et al. 2019 the hyperparameters chosen =b=l.
15

---

0,06
01
0.5
MMDXP-00
MMDXP.... Qu.
MMDDP Go.
0.70
Stein Thenning
A77
Megularze Sten Timiiii
0.72
110
300
3.0
7.00
1000
102
300
588
700
1100
100
300
0.0
710
1000
c- 0.05
0.2
-65
1.2
Qn))
MMDEP Gr.
MMDXP
1.0
Stein Thrminn
(egnlariied Stein Thimingg
088
Figure 12: (Gaussian mixture, MALA) Graphs of the MMD distance with respect the thinning size m (with d = 10) for various step sizes e. For the first row, we set À = 1/m", and we observe that the effect of entropic regularization almost vanishes, since the performance close the original Stein thinning. For the second row, we set À = 1/ log(m), violating the convergence assumption, and resulting in bad thinned samples.
Stein Thinning
Regularized Stein Thinning
A log(p)
30
20
0-
-20
20
20
20
20
0
20
(1)
x: (1)
11)
Figure 13: (t-banana-shapel mixture with d = 10, MALA) First two panels: solutions obtained with Stein thinning and regularized Stein thinning with contour lines of the target distribution. Last panel: heatmap the Laplacian correction + log(p) for (3) = (10) =
Table 2: Description of UCI datasets
Dataset
Sample size
Dimension
Breast Wisconsin
569
30
Diabetes
768
8
Haberman
306
3
Liver Disorders
345
6
Sonar
208
60
16

---

Stdn Thiming
1.6
Regultiool Siii Thinning
NMDDPP (g.
MMIXP
1:0
0.5
0.
80
000
300
230
200
500
1100
(U)
100
200
500
100
9000
1000
100
500
703
903
1800
LD
Nh-m Thrming
Regelaried Steis Thirping
MMDDP Qu)
22 Q
Q..)
11
30
30
10
50
100
13
201
10
150
101
70
50
50
70
60
00
100
Figure 14: (Mixture of t-banana-shaped distributions, MALA) Graphs of the MMD distance with respect to the thinning size m (with d 2) and with respect (with 300) for various step sizes E.
01
0.5
1.0
Sten Thrming
MMDXP Q-)
Thiznine
MMDDP..0 Q-.
MMD
100
100
230
200
500
1100
100
200
300
700
1100
100
300
1000
703
900
1900
B
= LD
St-in Thuming
Reppbaried Seett Thirning
MNDDP. Q..
MMDXP
MADDP
Figure 15: (Mixture of t-banana- shaped distributions, NUTS) Graphs of the MMD distance with
respect to the thinning size m (with d = 2) and with respect to (with m 300) for various step sizes E.
B Proof of Theorem 23
Assumption 2.1 (Distant bimodal mixture distributions). Let and q be two mixture distributions in Rdd made of two modes centered in (-14, 0a-1) and (j,, 02-1), with > 0. The distribution of each mode of p C² (IR)) has Rd as support, whereas each mode distribution of have compact support, included in a ball of radius r > 0, with r < p. The left mode of p has weight wp * 1/2, wy. Similarly, w and
w are the mode weights of Let Qr and Qn be the probability measures that respectively admit the density of the left and right modes of q. and IP and Qu be also the probability laws for p and q.
and the right mode has weight
Assumption 2.2. For distant bimodal mixture distributions q and satisfying Assumption 2.1 and for 17 E (0, 1), we have KSD? QL)/KSD² (P, Qx) - -1|
Theorem 2.3. Let kp be the Stein kernel associated with the radial kernel x(x,x') = o([[x x' e)) where x, x' E RRd, E. O, and & E (IR), such that Φ(z) > ('(z) + O, and c²" co. Let p and be two bimodal mixture distributions satisfying Assumptions [ 2.2 for
for
any E (0, 1). define w² as the optimal mixture weight of q with respect to the KSD distance, iee. 21177
w argmin KSD(P, Qw). Then, for p large enough, we have wn
w=(0,1]
17

---

Table 3: AUCs obtained by Stein Thinning (ST) and Regularized Stein Thinning (RST). A 10-fold cross-validation is performed and the experiments are repeated 10 times to provide uncertainties.
NUTS Sampler
m
50
m
100
300 m
Dataset
ST
RST
ST
RST
ST
RST
Breast W. 0.88 (0.020) 0.96 (0.004) 0.91 (0.023) 0.96 (0.003) 0.93 (0.008)
0.96 (0.004)
Diabetes 0.52 (0.009) 0.50 (0.019) 0.52 (0.021)
0.55 (0.018) 0.53 (0.015)
0.57 (0.019)
Haberman 0.51 (0.038) 0.53 (0.023) 0.54 (0.033)
0.58 (0.035) 0.53 (0.034)
0.58 (0.017)
Liver
0.53 (0.044) 0.69 (0.014) 0.56 (0.038)
0.70 (0.013) 0.61 (0.039) 0.70 (0.011)
Sonar
0.80 (0.021) 0.81 (0.007) 0.81 (0.009) 0.82 (0.011) 0.81 (0.011)
0.81 (0.009)
MALA Sampler
50
100
300
Dataset
ST
RST
ST
RST
ST
RST
Breast W. 0.68 (0.044) 0.93 (0.010) 0.72 (0.048) 0.93 (0.007) 0.72 (0.037)
0.88 (0.026)
Diabetes 0.51 (0.012) 0.48 (0.010) 0.53 (0.028)
0.51 (0.014) 0.53 (0.016)
0.56 (0.016)
Haberman 0.52 (0.034) 0.60 (0.027) 0.53 (0.024) 0.58 (0.017) 0.55 (0.024)
0.61 (0.013)
Liver
0.54 (0.033) 0.70 (0.008) 0.55 (0.034) 0.69 (0.005) 0.57 (0.024)
0.62 (0.032)
Sonar
0.80 (0.019) 0.80 (0.019) 0.80 (0.010) 0.80 (0.010) 0.81 (0.013) 0.80 (0.010)
Lemma 1. If kip the Stein kernel associated with the radial kernel k(x, y) <[[[x yl12/e). where& l O, Φ E C² (RR), and x,y R2 such that sp(x), sy < So, then we have
d-1
2so B' [1| yl12/E)
[x,((,)) y)1 5
f|xx - y112 o (l|x y112/0) 72 0 ((1x-xxle/E))
+ S x - yl12//).
Proof of Theorem B We consider the mixture distributions p and satisfying Assumption 2. for
J 0 and r > and assume that Assumption 2.2 is satisfied for v) E (0, 1). More precisely, denote by qr the distribution of the left mode of the mixture and similarly, OR is the distribution of the right mode of q. The probability measures Qr and QR respectively admits the densities qr and
By definition of the KSD, write
KSD (w) = (ip(x, x"
Additionally, given the above notations, g takes the form y wgr. (1 -w)gn- Then, we can develop the KSD expression to get
/ kp (x,xxx)) (wwq + (1 ())n(x))( (x') + (1 w))n (xx) dxxdx =w² / Xx (xdxxxx' + (1 (1-w)² / kip(x, x" (x))g((())))))))) ddxdx
KSD (P. Qw)
2wo(1 -w) w) / kp(x, x') )nn ( (x))n(()))))))())))(((())()) )dxdx',
where the last term follows from the symmetry of kp. Finally, we have
KSD2 Qw) = Qr. + (1 K)D QR) + 2wwl w) /
x (xxx ddxdxx
18

---

and we denote by AL.R the last term of this equation, which now writes
(4) We first focus on the last term LR R of this expression, which can be shown be arbitrarily small when μt gets large. According to Assumption [.1] the distance between the centers of the two modes 2μ., and both qr. and On have compact support included
KSD (IP, Qw.. KSD (P, @2) + (1 k) (I, Qn) 2wo(1 w)
ball of radius r. Consequently, for
x, x' E IR such that (x) > 0 and gn((x)) > 0, then ||x x'I22 2(μ r). Additionally, since the score Sp continuous, Sp bounded on a compact set, and exists Sg > 0 such that sp(X) < 8o and sp(x')
Then, from Lemma] 1 we have
d 1-1 20( r)
1 03 (1|x x x1112/E) + 12 p00 1ox-xx11 x11222)) 8
280
(x,(x, x'I S
o (|[1x - x x(12/0)
+ 1/0011x x 2212/)),
and since (m.xx)n((x* = 0 for || - x112 < 2( - r), we get
AL. S
sup
d-1
(=)) + p (2) e xx(z) + }o(=))
280
(2 r)
By assumption, ((2) + 0, Ø (z) + 0, and (z) + 0 for = xc, and then, we have
lim Ar. = 0. 10000
Next, we reorder the terms of equation
to get second-order polynomial in follows
KSD2 (FP, Qw) = 2 QL) + KSD2 (F. QR) 2
2w KSD 2 (P. QR) R + KSD2 (IP, QR).
Notice that the coefficient of w² is KSD* (IP, (1/2)/ 4, and therefore positive, Then, KSD2(I, Qw)
admits unique minimum with respect to given by
KSD (P, QR) AL.R
KS + (P, QR) - 22
We rewrite w* as follows,
1/2KSD* (P, Qn) + 11/2KSD" Q1) - AL, + 1/2KKD (P, Qn) 1/2KSD (P, Qz)
w*
KSD2 (P, Q2) + KSD2 (P. Qn) 24L.R
1 KSD2 Qn) (SP
+ 2 KSD Q2) - KSDD((P QR) - 2AL.R
2
1 1
1 KSD2 (T, (t)/kK() Qn)
2 21 + KFD Qt)/KSD QR) - 2AL.R QR)
1 1
(FS"(F/QQ)/KKDD' Qn)
+
2
2 2(1 AL,R/ ALLR/KSD' (IP, Qn)) (KSD) QL)/KSDD (I, QR) 1)
We can deduce the following bound
1
(PD Qr), QR) 1
w
(5)
2 2 LL. KSD (P. QR)) + (KSD (I,, Qt) QR) - 1)]
According to Assumption 2.2
with
17 1,
KSD2 (P, Qr.)
< <m
KSD 2 (IP Qn)
which gives an upper bound for the numerator of the right hand side of inequality ( Additionally, for large enough, RR is arbitrarily small, and in particular, we can have AR. < KSD2(P, QR)/2, and then 2(1 - AL. @R)) > 1. Next, we use the triangle inequality to get [2(1 Az, At.R/KSD* QR)) (KSD? Q (P, QR) - -1)1 > 2(1 At.R/KSD (P, QR)) - (I (r)/ QR) -
> 1 7
19

---

where the last inequality is obtained using 2(1 - AL.R/ AAR/ (P, QR)) > 1 for je large enough, and Assumption 2.2 again. Finally, this lower bound on the denominator and the upper bound on the numerator combined with inequality (5) give
< 77
2
2(1 11
Proof of Lemmal 7 We consider x, E Rd, such that y, (x)) < 8o and sp(y) < Su- From Equation (2) of the main article, we derive the Stein kernel obtained for a radial kernel @[| ||x hee-yles 12/(),, where & 0 and (R), and get
1 - d
xx = [|x y112
6'( - x-112/1) x2 02 000111-y11//)
(sp(x)) sply)) (x - y) 0'(1|x - -1112/0)
y) Using Cauchy-Schwartz Schwartz inequality, we have sp(x) -sply) S sg. and
+ (sp(xx) -pp()))((1xx - yyll//).
f||x y112
(8p(x) sp())) (x -y) I p
290
f|(x y112
Overall, we obtain the following bound
1xp(x,, y)1 S E||x y112 Φ (l|x yl12/0) (e(1x-ylle/11 + 280
1 d
d'([x yl122/0
+ yl1l/())
C Proofs of Theorem 2.4 Corollary 2.5 and Corollary 2.6
Theorem 2.4 (KSD spurious minimum). Let ki be the Stein kernel associated with the IMQ kernel
with 8 >> 0, 3 E (0,1), and e = 1. Let {x; 1 Mag {x E IR [|sp(x)]lz S 80 be fixed and m 2 We have small
set of points of empirical measure Qm m
- 2i=1 8(x:), with 8p
KSD² (P, Qm) < E[KSD2 (...P)), the score threshold so and the sample size /2)/(23d/f + 23so/ +.).
enough to satisfy m < + (E[ls,(())
Corollary 2.5 (Low KSD samples at density minimum). Let kp be the Stein kernel associated with the IMQ kernel with > 0,
(0, 1), and 1. Let p be a density with at least one local minimum or saddle point. For m 2 2, if x; 1m C R is set of points, all located at local minimum or saddle points of p. then we have KSD* (P. E[KSD* (F, (P,Pm)), if m < 11+ E[1s, (X)]
Corollary 2.6 (KSD spurious minimum for Gaussian mixtures). Let k be the Stein kernel associated with the IMQ kernel with & O, 8 E (0, 1), and e = Let the density p be a Gaussian mixture model of two components with equal weights, respectively centered in (-H, 0a-1) and (j,, Oa-1) 1), of variance o² Ia,,, and v = 1//o. Ifv
and 0 S 80 < [UVV22 1-In VV2 1)] pp, then for any )m C. .M., of empirical measure Qmi we have
(i) KSD? (IP Qm)) < E[KSD² (P,, P..)). )1 if m and sp satisfy m <1+ 28dJe" 22800/8+st
(ii) there exists three disjoint intervals I- Iu, 1 C R. respectively centered around -4. 0, and p.,
such that x (1)
(11 E I_ Uo U Up-
Lemma 2. Let kp be the Stein kernel associated to the IMQ kernel, with parameters & E> 0, 3 E (0, 1).
and c 1. For so 2 minxeRs [|s(x)llz and x, y E Man we have
2,3d 2,300 + s
(x,
Proof of Theorem 2.4 By definition of the kernelized Stein discrepancy between the target distribu- tion and the empirical measure P.,. -/ Σx 6(X;), we have
m
KSD (P. P..))
Σ
(X, X,)
m
JJ=1
20

---

In what follows, hpp denotes the Stein kernel obtained for the inverse multi-quadratrie kernel function, i.e.,
4,(33 + | + ||x - 8-2
kx,,,,)
23 + (d (x,(x) sy(r)) (x - y)) ||x y112/89) 8-1
+ (x) s1 + ||x y13/
Given that X; and Xj are independent random variables that follow the distribution IP, one has (X,, X,)] for any i * j. Using this property and the closed-form expression of the Stein
kernel Kp, it is found that
E IKSD
E (X, X
m 2=1
Eipp (X, X)]
m
2
E E11ls((X)1
ml2
where X ~ P. On the other hand, the kernelized Stein discrepancy between the target IP and the empirical distribution Qm 1/2 >x= 8(x;) is given by
KSD2 Qm)
1
Σ
(x, X3)
m
i,
1
1 m
kip [x, x;)
kp (xi, X3)
m
m
ifj
23dd
m
8p
(x;
kip( (Xi, xj)-
mf2
m
(
Next, it can be shown that the difference m [KSD (PP Qm)
E|KSD2 (I. P )l) takes the form
m(KSD Qm) (P,
(x,
ky (x,, x2) (X)
m
Since Xy M can Lemma to bound the terms kp (x,, X,). and then obtain
mKSDD (P, Qm) EP S o + (m - 1) ()(22)) 233d
2,38011
82 +
$8) - (X)) (x)
El|
where the right hand side is always negative if
Eilss((X) [1 $
< 1+
23 + 2350 $o
Proof of Lemma
2
The Stein kernel kp obtained for the inverse multi quati-quuadratic kernel function, > O. 3 (0, 1), and 1, is given by
with parameters
(x, = 43(B 1)
|xx - (1 + ||x - y113/19 -8-2
8
20 +
(d (,, - sy(r)) - x)) ||x y113/8)-P0 -8-1 8-1
+ (x) (1 + ||x ~11//8)--
21

---

Since the first term is always negative and (1 + |xx y11//e2 - - 1 for a 8, + 1, we obtain
23dd
23
kp(X, y) -
+
(xp(x)
sp(y)) (x y)](1 + 1 y12// (xpp(x) sp(y))
p2
82
23d
+
23
((p(())-8p()) -p((y)
x x-y
xx - y y122/E
12/ + 1()p()) sp(x))
||x y|12 (1 +
We define the function 9 for 2 2 0 as
g(2) (1+ 22)
and simple function analysis shows that
def
sup =20
()) V23
1
23 + 1 8+1
2,9 + 2
Since 3 E 1), we have S 1. We combine the last two inequalities for ( yx and e1 to get
yx. -
23d
+ (<p((x)
Sa0x))-8pp0)) -sp(())-sp()))
x y
+ [<xp(x)) sp(y))]-
f2
||x yl12
We can apply Cauchy-Schwand inequality, and since x, .May we get
[ssp(x) - sy(y))
y [1(sp(x) - sp(y))]l2
||x yll2
< 290+
xx y112
||x - 1|2
and also
(x)p(x) sp(y)) S S
Overall, we obtain
2,dd (2
2,300
kp (x.
+ sg. L
2
Proof of Corollary 2.5
As minimum and saddle points are stationary points of p, we have
{x C Mo,
and we can apply Theorem 2.4 for so = 0 to get the final result.
Proof of Corollary 2.6 Let the density p be Gaussian mixture model of two components with
equal weights, respectively centered in (-p, 0a-1) and (n, 08-1). and of variance I'Ia. and let
u/o. We assume that > and 0 < So < [VVE2- In(v + v2² 1)]
1)]/p. /p. Then, according
to Theorem 2.4 for any {x; C Ma., of empirical measure Qm.
8(x,), we have
(i) KSD* Qm. < E (P, P PP)) if m and 80 satisfy m < 1 + BIIl, (XII- 23d/ 2380/E44
To prove statement (ii), we need to characterize the shape of the set Map C Rdd, given by the level lines of the squared score norm (x) The density p is Gaussian mixture, i.e.,
p(x) 2(2x)d/20d 1|2/200 (e (w)) /20²
where (-1) is the vector x without the first component. Then, the score is also given by an explicit formula.
(2)
tanh (
sp(x)
where tanh is the standard hyperbolic tangent function. An important property of this score function is that the j-th component of Sp only depends on 0,( which makes (x) invariant by any translation orthogonal the j-th axis. Then, can compute the squared score norm
(sp(x)]12 s ())) ()) +
04
22

---

75
(x
5
Lo
80
6
2
min
2max
(1)
"max
min
Figure 16: Squared first component of the score for a Gaussian mixture with 3 and 1.
where 8/ (2)²
2/r tanh( (zz )=² simple function analysis of this univariate function, illustrated in Figure 16 shows that sp (z)? has two local maximum in 2max and 2 and three local minimum in 2 min O, and miin provided that 11/a > We also get that 69
=
(2) (z)² grows to
to0 when (1) + +//- 20. The extreme values are ordered as follows
1u mmin max 0 < 2 < = <
J..
(1) The values of mmn Fmaxx Z and z given by the zeros of the first derivative of sp ())( defined by
dsp (2) 2( 1
-
H
tanh
cosh
This derivative vanishes when one of the two factors null. Since u/o 1, the first term is null when
4²²
1
1,
cosh
2)2
which leads to
o2
o²
arcosh
P
and
+
arcosh P
max
mmax
The second factor is null when
tanh
O.
0
1
Obviously, 0 is solution. Since u/o > 1, equation they do not have closed form, have
(6) has two additional solutions. Although
Z E (-14, ²
(2 1)).
Also notice that, as u/a gets larger, amin is closer to 11 and z to μ. For example in Figure 16 set u1/o 3, and we hardly gap between - and Z μ and a
23

---

By definition, for x E Ma we have
(100)) < |sp(x)] - s
Therefore, given the variations of sy (1) (x(1) detailed above and illustrated in Figure 16 if s < Sp (=max)² 8p (=max)² 2 it exists three disjoint intervals Tpp To la R, respectively centered around -M, 0, and μ4, such that
x Mep
20 E
To conclude, compute the value of (2) (mmaxx) that is
(Emaa)) = aaacosh ( + 1 o²
tanh (arcosh( ())coss())))
10
Using the formulas (anh(arcosh(x)) = 2x2-1 for |a|| 1, arcosh((r) = --(nn vx 1), and with V u/o oo, we get
(1) (= = =0 (1/0// - 1/0 inn//o+ // (1//0)=-1) -1/1
[[VV22 - 1-nnn( 1)]11)
which is always strictly positive since > 1. By assumption, 0 S 8D [WV - -nnv + v22 - T) /u, and therefore, we have 80 < (2) (=ttax = op (zmax). which concludes the proof of statement (ii).
D Proof of Theorem 3.3
Theorem 3.3. Let kp be the Stein kernel associated with the radial kernel k(x, [x[|xxx\\l2// where x, x' E Rdd E> O, and & (RR). Let p and q be two bimodal mixture distributions satisfying Assumption 2.7 We define W³, the optimal mixture weight of q with respect to the entropic regularized KSD distance, i.e., Wx = argmin KSDx Qw). If [[og(((((Z)) + ([llo(((ZRR)) wE(0.1]
where ZL ~ Qr and ZR On it exists X E RR such that w² Wpp
consider the mixture distributions p and y satisfying Assumption 2.1 for >> and > O, and denote by qr. the distribution of the left mode of the mixture q, and similarly, aR is the distribution of the right mode of q. The probability measures QL and QR respectively admits the densities qL and QR- As in the proof of Theorem 2.3 get that
Proof of Theorem 3.3
KSD2(IP Qw) = KSDD (IP, Qu. + KSD (P, QR) - 2AI.
2w KSD [IP QR) A + KKSD2 (R.
Next, we define Qui Zr Qr. and Z QR. and develop the entropic regularization term, E[log(p(Z))
/ log(p(x)) + (1 (x)))((())ss
/ loglo(x))) (1- w) / [og(((()))Rxx))dd
wB[log(p(Z [)] (l)og((((ZR)))
Combining these two results, have
(PD Qw. = [KSD2 Qz) + KSD2 Qn) 24
2w [KSD (I, Qn) AL. n 3/2(E((og(((Zz))) - ([log(p(Zn))))
KS, QR) (p]log(p(ZzR))
We recall that p does not depend on w, but only on the fixed weight wp and the two mode distributions. Consequently, KSD (IP, Qw)) is second-order polynomial with respect to w. As for Theorem 2.3 the
24

---

coefficient of w² is KSDD Q 21/22/44 and is therefore positive. Then, the polynomial is minimized with respect to w, at wn defined by
KSD (P, QR) - AL. + 1/1] [[])g(((ZR))
"
KSD Qt) + KSD (IP, QR) L
Since [[to)) ([log(p(Znn)) + by assumption, we get that w up for lambda defined by
X = 2 u, (P,
w. k, (P, + 2wR)
E Bloo(p E (pg(p(ZR)) ZR))]
E Proof of Theorem 3.4
Theorem 3.4. Let ki be the Stein kernel associated with the IMQ kernel with 8 0, 3 E (0, 1), and c= 1. For m 2 2, lel (x:)m C IRd be set of points concentrated at Xo, local minimum or suddle point of p. and of empirical measure Qm- Then, we have L-KSD* (IP Qm) > ELL-KSD* (I,, Pm)], ) if
the density Xy satisfies p(xo) >(x(0))) [Ell) (X) E A log ((X)))
Proof of Theorem 3.4 Following the same approach as in the proof of Theorem 2.4 we have
(FLKSDD((F,IP.))) = 23d
E Elll
pp()s((X) (X)
E + lxg)(X))
m82
where X I., and, with the empirical distribution Qm
1 E. 0(x:).
23d
122
m
L-KSD2(P,.OA
(x,))
(x,,
+
log(p(x:)).
mf2
m
m
m
As (x: (x:) are concentrated at local minimum saddle point Xo, the score is null for all particles, as well the distances between them, and we get
23d
(m
1)23d
A log(p(xo))
L-KSDD Qm)
mf2
ml2
m
Next, the difference m (L-KSD (P, Qm)
E LSLS (P,, (T)Pmm)] writes
m(L-KSD² Qm) E|L-KSD (P))) A log(p(xo)) + (m 1) 28d
(E|l + B log (())).
By definition,
A log)(x) def Σ
82 log p(x) Oxli)2
with
8 logp(x)
= 8p((x)
(11 Op(x) 2
Ozli22
p(x) Orli2
p(x) azols)
As Xy is stationary point of p, 0p(xo) (zz = 0, and obtain
82 log px))
p2x(x)
Oz0i2
p(xo)
Oxl12
leading to
A logg((x))
Σ
1
@'p(xx)
A Atp(X))
pl p(xo)
p(xo)
j=1
25

---

Finally, we get
p+p(xo)
23d
m (LL-KSD((P,Q.))
(P,
(m {
p(xo)
(B[I)sp(X) + E|A (x))(
which ensures that L.KKSD (I, Qm)
TE [L-L-S (TP, Pm]]
> 0 for 2 provided that
A 1 (xo)
p(xa) <
E8
s
2 + E| [xg](X))
F Proof of Theorem 3.6
Assumption 3.5. Let Q be probability distribution IR such that IP is absolutely continuous with respect to Q Let {Z, hie C I be a Q-invariant, time-homogeneous Markov chain, generatd
using V- -uniformly ergodic transition kernel, such that V(x) 2 dl 2/3d//2 + [sp(x)]l Suppose that, for some 7 > sup Ble log(p(Z,))]
oo, sup El + (Z))
< 00.
VEN
ieN
2 sup E(e E max(1, (1, aq (z,))) 2 sEN
< 00, sup EEN
E
dP
(Z.)
23d
+ [1sp(Zz))BBV
< 00.
L2
Assumption 3.5 is close to the assumption made in Riabiz et al. 2022 For the last two finite expectations of the assumption, notice that we have just plugged the formula (x, x) = 2,d/E2 [xs](x)] into the integrability conditions, since this formula is quite straightforward. But in addition to Riabiz al. [2022). require two integrability assumptions, for each regularization term. We give additional insights about these quite complex conditions below. Overall, our assumptions are hardly stronger than those of Riabiz et al. [022]].
The condition E|e21 log()(Z.))] < oo is satisfied if f(x)p(x)-ddz < oc, where f is the distribution of Z; (since p(x) 1 only on compact set and p is continuous). Therefore, there exists y satisfying such condition, provided that the tails of density f are not too heavy compared to the tails of p. This has be verified for all iterations of the Markov Chain, and happens be quite mild assumption.
Regarding the condition involving the Laplacian term, (2022) (Appendix S2.4). relies on a Lipschitz condition for the score function V log p. to transform the first finite expectation of Assumption
the analysis from Riabiz al.
amenable and practical integrability condition which writes [[ < oo for some K E (0, oo). In the same spirit, if we add Lipschitz condition on A log p. second additional integrability assumption reduces a similar amenable condition. Finally, even though Riabiz et al (2022) do not elaborate further on this condition, it can be noticed that it is satisfied if all Z; have subgaussian tails, for example. Besides, one of the main example of distantly dissipative distributions are log-concave functions outside of a compact set. In this case, the Laplacian regularization is null outside of a compact set, since the Laplacian is negative for concave functions. Then, the integrability condition is automatically satisfied for this type of distributions.
3.5 into
Theorem 3.6. Let TP be a distantly dissipative probability measure, that admits the density p E C² kp be the Stein kernel associated with the IMQ kernel wh f,c > 0, E (0, Ler (Zihie C. Rdd be Markov chain satisfying Assumption 3.5 be the index sequence of length mn generated by regularized Stein thinning, and Qm., be the empirical measure of (Z mm If log(n)a < Mh < n, with any >> 1, and Am., = ooloooommmmm (onsmmmm//mma then we have almost surely
Qm.. x
Theorem 3.6 extends Theorem from Riabiz et al. 2022] to regularized Stein Thinning, us-
ing Lemmas and assuming the following convergence rate of the regularization parameter Am, olog m., /.,). LemmaBJalso extends Theorem from Riabiz et al. [ 1],
whereas Lemma is slight modification of Lemma 5 from Riabiz et al. 2022
26

---

Lemma 3. Let IP be a probability measure on Re that admits density (RR)), kp be reproducing Stein kernel, and {x;}i=1 C. I a fixed set of points. If π is an index sequence of length m produced by regularized Stein thinning, then we have for O.
wi ((x;)) 1 + log(m)
KSD Σ o(x,,)) = KSD2
max kp xx) i=
2=
+ log(m)
A + log(p(xi)) 2A max
[log(p(x;))1
max El....
m
==,...
where the weights w are defined
w" E arg
min
KSD*
Σ wx((x:))
Σ; w. =1 W; 20
Lemma 4. Let be non-nes non-negative function on R Consider a sequence of random variables
(Xi)ee C Bd such that, for some 7 > 0,
b def
sup E(e oo.
iEN
Iflog(n) < mm < n, with any a > 1, then we have almost surely,
log(ma)
lim
max
f(X.) 0.
12 12-000
mh.
The proofs of Lemmas 3 and are reported the end of this section. We first proceed with the proof
of Theorem 3.6
Proof of Theorem 3.6
From Lemma
we have
mm
log(mn)
KSD?
8(Zz,)
< KSD?
i ((z)))
+
max EI.....
kp (Z,, Z,)
My
j=1
+ log(m,)
+
max (=.....
log(p(Z,)) + 2Am.. max
log(r(Z.))
mp
Riabiz et al. [2022 [Proof of Theorem 3, p.12] showed that the term (+) converges towards 0 almost surely as > oc. For the remaining terms, from Assumption 3.5 we have
sup 1/c1l///*+11,,222)0 oc, EeN
sup E[e Ef (oggp(2g)(21) < <00; teN
and
sup 7A log(p(Z.))] 00. ieN
We can use Lemma 4 with fxx = (x, x) and f(x) = A log(x) to deduce that (**) + 0 and (k (****) 0, respectively, The remaining term (o) be rewritten as
27m max log(p(Z.)) 2mmdmm, log(mn) log(mnn)
max log(p(Z.))) i= I......
i=.....
ma
Using the assumption that Am.. = o(log(mn)) ) /m,) and Lemma with f(x) = log(p(x))), we
conclude that (o) - 0. It follows that KSD I Ey 8(x,,)) + 0 almost surely as n - co.
m., Given that p assumed be distantly dissipative, we apply Theorem from Sam al. 2019 to obtain that Qm.. * P almost surely, +
27

---

Proof of Lemma 3 We consider an iteration t E f2....,mi m. of regularized Stein thinning, where 2 is the final length of the thinned sample. We define a PKSD وj= 8(x,:) and ft 1 - kp (x(xx,,)). We also denote Sf = max,=1...... ky,,, x;) + max:=1.... A log (),))( and S2 2 max log(p(x)))|- Using the definition of the squared KSD of an empirical measure, we have
a: = 0c-1 + (x(xxx Xn.) + Σ Kip xxxx
j=1
Let x = arg m ized Stein thinning algorithm at iteration and we have
(x. 1._, fy- By definition, xx, minimizes the cost function of the regular-
(x )xx) + 2 (xx, xx.) + At log(p(x=:)) - At log(p(xxx))
1=1
$ xx,,,xx)) + 2
_ (xx, xf) + A log(p(x))) At log(p(x)))
We combine this last inequality with the first equation to obtain
a: a + kp x;)
ki (xa x,) + A+ log(p(x;)) A log(p(x<,))
Xt(log(((x))) - log(p(x...)
and then,
6-1
at S a xi) +
Σ j=1
kp x() ) + - log(p(x;) A log(p(x..))
+ At( log(p( x )) + log(ppx>;,)))
By definition, hence, we have
S log(p(x:)) S S2 and (x x
+ A log(p(x;)) - A log(p(x=()) S S.P.
a: S at-1 + St tAS2 +
min fy-1
y (
As in the proof of Riabiz al. 20222 we have miny Pin J S /ot- A' ((()))))) where h* the element in the RKHS (k,,) of the form Li=1 kp (xi: .). As result, a, is bounded as follows:
a S a + S S + +1S²- -2 h (12(0)))))
We then show by induction that
ag S (11111112 + Cc + +S2))
where
C2 def (s- |1 ((())))
t
With such as a result, we will have KSD2 (P. - 1= 8(x,,) )) S KSD2(IP, (IP, Σn wfd(x )) + C + XS2 and obtain the advertised result in Theorem 3 for the last iteration = m.
For t = 1, we have a1 = kp( xx,) S S7 and thus - h- l( +C1+AS For fixed 2 2, 1-1 1
assume that a,-1 <
(1(e,) +O where C,-1 = [-1 (Sp - 1 )
クj=1 ]
We then have
a, Sa,- + S9 + t2S22 + 22 12" (10-111(111((()))))
<(t 1)2( 11 + CC-1 A82) + ST LAS2
+ 2 1)
|1h- /11~112200,) + Cc-1 + XS21
-((( + C. + AS2) + R4
28

---

where
Re= (1=111))0-11-P
21) (P( X52) + sी
+ 1AS1 + 2(1 - v) ||+ 11 + C.-1 + AS2
-1- - -. + (1 (0)]11"11 + S
21 - 1) + 2(t - 1), || /11+11 C--1 252111"
Using Riabiz et al.
2022 Lemma 1], we have
h (((k))) | Thack + Cc-1 A82 S h 22 + Cc-1 AS2
It follows from Equation
9 that we need R, S 0, i...
P. (t (t- 5t 1-11h* 111h1| (((p)) A82 + + (1
2 2| (e(,,)) C--1 + AS S
t-1
The above inequality is always satisfied long
2 (R(a,) + C--It PCC, 1)-0-1
-7-11h"1
+ +211h"1
--
t
which is equivalent
1C, - -t- 11- ( |1 [1(
and always true by definition of Ct. Hence we have shown that at S 82 (1 ((( + + AS2).
Given that / (eeck,))
KSDD(IP, Σi (x;)(, we have
KSD (IP, Σ 8(xx5)) & KSD (P, Σ wii wjo(xx)) + AS2,
where C. < 1+log(t)
maxx (x; x;) + max A+ log(((())))) (see Riabiz et al.
2022 [emma 2)).
We follow the proof Lemma from Riabiz et al. 2022 the last step the proof, we essentially need to show that
Proof of Lemma
ca (m 00
and c2(mm) < <000
where
(m def 2 2log(mm) log(nb)
and 02 (2(m)) def 4
log(ma) log(n((m, + 1)2)6)
m²
m
With the assumption that log(n)3 S mn < n with 3 1, it is deduced that ei (m?) - and cg(m,) + as n + co.
G Laplacian Stein Operator
Instead of the standard Langevin operator, we can use Tx(g) = Appg)/p mentioned in Oates et al. [2017, Appendix A.2]. However, such Stein operator introduces similar problems Pathology I, since the Stein kernel associated with (g) also has spurious minimum in regions where second derivatives of p vanish, illustrated below. Therefore, it is more appropriate to taylor specific Laplacian correction as proposed in this paper, which cannot directly be derived from a Stein kernel, study the operator Tμ defined Tpg A(pg)/p and such that (Oates et al. 2017 Appendix A.2]
since it is not differentiable. In this appendix,
E (T,g)) (Z)] = 0,
29

---

for all g belonging to G, and with Z ~ IP. For each dimension j E let T be the operator defined as
1 p(x)
(Tx9)(x)
V2 p(x)((x) + 2V 2; (x) p(x) V 9(x))
for g : Rd + R. and where žj is the j-th coordinate x to simplify notations. The operator 7p can then be rewritten as (x)g)(x) ンラー (x)g)(x) Gorham and Mackey [ 2] establishes the closed-form expression of the KSD in the case of the multidmmnniooall Langevin operator. We generalize the proof of Gorham and Mackey
2] for the Laplacian operator (x) = A (((x)g(())) (x))) and establish closed-form expression of the Stein kernel kp R& x R + R. For a given kernel function k Rd Rd R, k E C the Stein kernel kp is given by [Gorham and Mackey 2017]
d (x, = Σ 5=1
1 (),,
(10)
where
x) = (TB (Tz (1:(x, ))s T3 (x(:,)))))
(11)
After a few developments, it is found that
=
+ 2V ,,,,,,
((x.))
x(x,y)
+ 2V 2V
+ 2p(y V - V.. kxx,) (12) k(x,y)
P(x - V ,x, + p V M(x,y)
+ 2p(x)
A closed-form expression can then be obtained in the case of, e.g., an inverse multiquadratic kernel of the form k(x,y) = (1 + |x y / ())12 where denotes the bandwidth. this case, one has
112
V v, ,(, k(x,)) 41), Ve,k(x,,) - k
V M, k(x,)) V2 k(x,y) = - k(x,y)3 3 (xx vz)) 1
2 Va, k(x,y) = Ø y)3 e (xx, (I) - 03)²
(13)
(I, -13) 15k(x,y) WA, k(x y) = -9((x,)) y)
(zz y/
2
2 By k(x,y) =-V, By Vy, kxx,)
M k(x,y) = y)( 90k(x,y) y) (x, v;)
105 y)
(xz Ui
84
g6
ps
A closed-form expression of kp can then be obtained by combining Equations
10 13). In contrast
the Langevin Stein kernel (see Equation (2) of the main article), the above Stein kernel involves second-order derivatives of the density p.
Example of Gaussian mixtures for this new Stein operator Tog = (pg)/p. We sequentially set = 2 and μ 5, with a = 1 and w 0.5. Results are displayed in Figure 17 In the left panel with μ 2, see that Pathology II does not occur, as opposed to Figure 2 of the main article with the Langevin Stein operator. However, when μ is set in the right panel, all particles are concentrated in spurious minimum again. Therefore, introducing higher-order derivatives of the target density through the Stein operator does not seem to be a promising route.
We run experiments based
30

---

3
2.
2
3
4
5
10
Figure Stein thinning with Stein operator Tpg = A Appg) /p. for Gaussian mixtures of Example with μ = 2 (left panel), and = 5 (right panel), a 1, and w 0.5.
31


Reviews:
[{"review_id": "nAAScXZuLo", "reviewer": "Reviewer_4eLX", "summary": "The goal of Stein thinning is to post-process the outputs of Markov chain Monte Carlo (MCMC) methods by minimizing the kernelized Stein discrepancy (KSD) between the produced chain and the target distribution. It is useful and has become fast a quite popular method in Bayesian inference because it automatically removes the burn-in period, corrects the bias, and has asymptotic properties of convergence towards the target.\n\nHowever it suffers from two pathologies, namely Pathology I of mode proportion blindness and Pathology II of mode collapse. Pathology I originates from the score insensitivity to mode weights, while Pathology II arises from the over-representation of a few modes. To mitigate these pathologies, the paper proposes a regularization technique called regularized Stein thinning based on adding an entropy and Laplacian term to the KSD. ", "questions": "Do the authors have any intuition on a general recipee or dependence of lambda on the parameters of the problem?", "limitations": "The authors reasonably discuss the limitations of their work. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper does a good pedagogic job on the issues related to the optimization of the KSD, and offers a complete picture by providing theoretical negative results about the regular KSD (eg Theorem 2.3, 2.4).\nFor instance Th 2.4 (through Cor 2.5) says one can attain a better KSD than empirical samples of the true target mixture by putting all the points between the modes because this is a region where the score is close to zero.\n\nMoreover the proposed regularised KSD is shown to enjoy good properties regarding Pathologies I and II and has theoretical guarantees and extensive experiments to demonstrate its efficiency.", "weaknesses": "Theorem 2.4 is still limited to the IMQ kernel.\n\nThe Laplacian correction involves second-order derivatives of the target.\n\nTheir regularised KSD introduces an additional hyperparameter lambda that a user should tune, however it is not clear how to fix this parameter in advance (depending on the dimension, the target etc). It is not clear in the theoretical result Theorem 3.2 neither what lambda recovers the true mixture weights and fixes Pathology I. Still, the authors conducted a reasonable number of experiments with different values of lambda. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you again for the review and the positive comments. We hope that we adress the identified weaknesses in the global answer above."}, {"review_id": "a2WXipT7oY", "reviewer": "Reviewer_zNe7", "summary": "This paper proposes a regularized version of Stein thinning for post-processing the output of Markov Chain Monte Carlo algorithms. The regularization addresses two common challenges of Stein thinning: insensitivity to mode proportions and samples concentrating at stationary points. The paper analyzes the above two pathologies under a theoretical framework which is then utilized to justify the proposed regularization technique. Numerical experiments demonstrate gains over standard Stein thinning in both toy examples and in a logistic regression setting.", "questions": "- Why do we see an initial sharp increase of MMD with respect to $d$ followed by gradual decrease in the Gaussian mixture example (Fig. 5) for RST, but strict decrease in MMD for ST?\n  \n- Is there a systematic way to set $\\lambda$ beyond the heuristics provided in the paper?", "limitations": "Limitations are not discussed.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The paper is well-written and easy to follow. The structure of the paper is logical.\n  \n- I found the study of pathologies illuminating and interesting. The theoretical framework appears to be sound.\n  \n- The proposed technique shows promise both in the toy case and in the simple logistic regression setup.", "weaknesses": "- More intuition on the assumptions of the theoretical formulation would be helpful, with some discussion on the limitations of the results with more focus on practical applications of the method.\n  \n- Extra hyperparameter is introduced, for which a heuristic tuning method is provided, however a systematic approach is missing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you again for the review, and the questions raised. Notice that we tackle the main points in the global answer above.\n\n**MMD variations.** The variations of the MMD with respect to dimension d in Figure 5 indeed depend on the considered examples, with different observed behaviors for the Gaussian mixture and the banana-shaped one. This is quite difficult to interpret, since these variations are due to the intrinsic complexity of summarizing a high-dimensional distribution with just m=300 points (as in Fig. 5), and is thus case dependent. But we agree that it would be interesting to investigate more precisely what distribution properties could imply such differences."}, {"review_id": "lbGKzyDf6O", "reviewer": "Reviewer_XEkq", "summary": "In this article the author(s) studied the Stein thining method, an algorithm for post-processing outputs of MCMC based on the kernelized Stein discrepancy (KSD). This article first theoretically analyzed two pathologies of KSD, and then proposed methods to mitigate the two issues by regularizing the KSD objective function. Both theoretical analysis and numerical experiments show that the regularized Stein thinning method has improved the sampling quality.", "questions": "1. Both Theorem 2.3 and Theorem 3.2 are based on radial basis function kernels, whereas Theorem 2.4 and Theorem 3.3 are using  inverse multi-quadratic kernels. What considerations is this based on?\n2. Assumption 2.2 seems too hard to verify. Can you give an example for it to hold? Considering Example 1 of this article with only one dimension, is it possible to given an expression for $\\eta$? I would expect $\\eta$ to be a function of $\\mu$ and $w$, and $\\eta\\rightarrow 0$ as $\\mu\\rightarrow\\infty$.\n3. As mentioned in the weakness section, the computing time may be taken into account and reported in the experiments.", "limitations": "After rebuttal: my previous concern on computational cost is no longer a major problem, and other limitations may include the tuning of $\\lambda$, and some technicalities such as the choice of kernel functions.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The blindness of score-based methods to multimodal distributions is a long-standing issue, and this paper presents a careful theoretical anaysis in the context of Stein thining. The author(s) also point out another pathology of KSD caused by spurious minimum of the objective function. Both analyses deepen our understanding of KSD-based methods. The proposed regularzed Stein thinning method also shows promising performance.", "weaknesses": "1. For the theoretical analyses, some questions need to be clarified. See the questions section.\n2. For the numerical experiments, the target densities explored so far are somewhat simplified. I would suggest considering the following additional scenarios:\n    - Sampling from energy-based models, e.g. [1].\n    - Gaussian mixtures with more than two distant modes, e.g. the experiments in [2].\n    - Bayesian neural networks, as an extension of the linear logistic regression.\n\nThese are all \"difficult\" target distributions, and I would not expect regularized Stein thinning solving them perfectly. But exploring at least  some of the cases above would be helpful.\n\n3. If my understanding is correct, the truncated Laplacian operator requires computing the eigen decomposition of the Hessian matrix for every point in the MCMC sample, and this seems to be a huge computational cost. (After author rebuttal: **this point was my misunderstanding**, and the actual computing cost was clarified by the author(s). The truncated Laplacian operator only requires computing the diagonal elements of the Hessian matrix.)\n\n[1] Che, T., Zhang, R., Sohl-Dickstein, J., Larochelle, H., Paull, L., Cao, Y., & Bengio, Y. (2020). Your GAN is secretly an energy-based model and you should use discriminator driven latent sampling. Advances in Neural Information Processing Systems.\n\n[2] Qiu, Y., & Wang, X. (2023). Efficient Multimodal Sampling via Tempered Distribution Flow. Journal of the American Statistical Association.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you again for the detailed review, and the many suggestions provided. We hope that we tackle the main points in the global answer above. We answer to the remaining specific points below.\n\n**Assumption 2.2.** Notice that Assumption 2.2 is satisfied when the two modes of $q$ have a similar KSD distance with respect to $p$. In particular, this is true for symmetric distributions. Unfortunately, it is not possible to give an expression for $\\eta$, since it strongly depends on the distributions $p$ and $q$, and the result is stated for any continuous distributions $p$ and $q$.\n\n**Additional experiments.** Actually, we conducted unreported experiments with Gaussian mixtures with several distant modes as in [2] and in a setting similar to the first energy model considered in [1]. They show thinned samples of good quality even for such difficult target distributions, as you can see in the attached pdf in the global rebuttal. We will add such cases in the article to show the good behavior of the regularized Stein thinning. Thank you for this suggestion.\n\nAnother suggestion is to conduct experiments with Bayesian neural networks, which is definitely an exciting idea. However, the scope of such application is very large, involves ultra-high dimensional distributions, and thus deserves a full separate paper in our opinion, since this setting has not been considered yet in the previous KSD literature."}, {"review_id": "0mwMqIHx0l", "reviewer": "Reviewer_GM6g", "summary": "This paper explores two pathologies associated with KSD thinning: the inability to distinguish mixing weights and concentration on low-probability regions of the target. The authors provide formalisations of when these pathologies can occur (Theorem 2.3 and Theorem 2.4), supported by empirical evidence. To address these issues, they propose the addition of two penalty terms: entropy regularisation and a Laplacian correction. The effectiveness of the proposed regularised KSD approach in mitigating these pathologies is demonstrated through both theoretical studies and numerical experiments.", "questions": "**Q1. Entropic regularisation**: Could you elaborate on why the entropic regularisation helps alleviate Pathology I? In particular, L242 says “The main idea of this entropic regularization is that $− \\log(p(x))$ takes higher values in modes of smaller probability, and therefore provides the relative mode weight information”. Whilst I agree with the first half of this claim, I am unsure how this sensitivity to low-density regions relates to the ability to detect mis-specification of the mixing weights, and why this regularisation would lead to small values when the mixing weights are correctly specified.\n\n**Q2. Figure 3**: The proportion of particles selected by L-KSD that lie in the left mode, while better complying with the true proportion of $w=0.2$ compared to the result of KSD, still differs from the true value by a statistically significant gap (since the reported estimated proportion is $0.11 \\pm 0.03$). Could you explain whether this gap is expected? Does it mean that the regularised penalty suffers from a significant bias in the finite sample scheme, despite being consistent in the infinite limit?", "limitations": "More discussions on any limitations of the approach should be discussed; see the question in \"Weaknesses\".", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "**Novelty**: Although the pathologies of blindness to mixing proportions and spurious concentration to low-probability regions have been acknowledged in the community, this paper stands out by providing formalizations and solutions to these issues. Solutions to Pathology I has been previously proposed in other applications of KSD (e.g., Zhang et al., 2022; Liu et al., 2023), but no solution nor formalisation have been provided regarding Pathology II. The proposed regularisation for Stein thinning represents a novel and principled solution, supported by strong theoretical guarantees.\n\n**Applicability**: The paper thoroughly discusses and provides experimental evidence on the impact of the regularisation parameter, the degree of freedom, on the performance of the proposed approach. Practical guidelines are given to aid in choosing appropriate parameter values in practice.\n\n**Writing**: This paper is well motivated and clearly written. \n\nZhang, M., Key, O., Hayes, P., Barber, D., Paige, B., & Briol, F. X. (2022). Towards healing the blindness of score matching. *arXiv preprint arXiv:2209.07396*.\n\nLiu, X., Duncan, A. B., & Gandy, A. (2023). Using Perturbation to Improve Goodness-of-Fit Tests based on Kernelized Stein Discrepancy. ICML 2023.", "weaknesses": "This paper is well written overall. I only found a couple of minor areas that could potentially be improved:\n\n**Limitations**: It would be beneficial to include more in-depth discussions on any limitations of the L-KSD approach. E.g., whilst a convergence guarantee was given when $n \\to \\infty$ (Theorem 3.5), for finite $n$ and $m_n$, how prominent is the bias introduced due to the addition of regularisation?\n\n**Clarity**: Whilst this paper is mostly clearly written, some parts deserve elaboration and clarification. E.g.,\n\n1. In Section 3.2, the intuition behind why the entropic regularisation solves Pathology I is not very clearly explained (also see Q1).\n2. In Theorem 3.3, it would be helpful to clarify the precise meaning of “concentrated at $x_0$”.\n3. Also in Theorem 3.3, it might be beneficial to provide an interpretation of the inequality for $p(x_0)$ on L313. Specifically, is this condition satisfied by the mixture of Gaussian example used in the paper?\n4. In Section 4, it would be helpful to clarify the choice of the kernel (and any hyper-parameters) in the MMD used to evaluate the results.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you again for the detailed review, and the questions raised. We anwser to the main points below. See also the global rebuttal for additional insights.\n\n**Theorem 3.3.** Regarding Theorem 3.3, “concentrated at $x_0$” means that all particles are located at $x_0$. We will improve clarity of this statement in the final version of the article. Thank you for the comment. The inequality of $p(x_0)$ tells us that the Laplacian regularization penalizes particles located around stationary points of the distribution in low probability regions. Since several intricate terms in this inequality depend on the distribution $p$, it is unfortunately not possible to give a precise interpretation of the threshold. However, many experiments show that such regularization is indeed efficient to avoid Stein thinning solutions with particles located in low probability regions. In particular, this is true for Example 1 with Gaussian mixtures. While it is possible to derive formula in this specific case (in the same spirit than Corollary 2.6), their interpretation is difficult.\n\n**Entropic regularization.** Modes with smaller weights take smaller density values, and are therefore more penalized by $–\\log(p(x))$ than modes of higher weights. Therefore, with such entropic penalization, regularized Stein thinning tends to select particles in modes of higher weights more frequently than in modes of smaller weights, and we recover appropriate proportions.\n\n**Mode proportions in Figure 3 and finite sample bias.** The mode proportion achieved by L-KSD, although much better than the standard KSD one, is indeed not the true proportion in this example. This is not due to a significant bias from the regularized penalty, but this is related to our heuristic for $\\lambda$. As mentioned in the global answer, we only have the guarantee that an optimal $\\lambda$ exists to recover the true proportion, but unfortunately, a precise tuning is impossible, and our rule of thumb does not ensure such exact recovery. It only gives a default value, which consistently outperforms the standard KSD, with much closer proportions to the target distribution, more representative of the true balance between the modes.\n\n**MMD settings.** All details about the MMD settings are given in Appendix 1.2. We will add these details in the main paper thanks to the additional page if the paper is accepted."}]

Author rebuttal:
We greatly thank the reviewers for their positive comments about our article and relevant suggestions. We explain below how we will improve the article clarity following the reviewer guidelines. We tackle the main points below, and also provide specific answers to each reviewer.

**Computational complexity.** The truncated Laplacian operator is simply the trace of the Hessian matrix, where negative components are set to $0$, and no eigendecomposition is needed. Therefore, the computational cost of the regularized Stein thinning is similar to the original Stein thinning. We will improve clarity of this important point in the final version of the article.

**Tuning of $\lambda$.** In our Bayesian setting, it is unfortunately not possible to tune the parameter $\lambda$, since there is not metric to assess the quality of the samples obtained over a grid of $\lambda$ values. This is the exact same limitation encountered in practice for the bandwidth of IMQ kernel $\ell$, for which, to the best of our knowledge, there is no satisfying and systematic tuning procedure. However, as opposed to the bandwidth parameter $\ell$, we are able to provide a simple heuristic for $\lambda$, directly guided by the careful analysis of Theorem 3.5. In fact, it turns out that setting $\lambda$ to $1/m$ is consistent in practice, and we show in various experiments that this heuristic is highly efficient. Additionally, we provide two theoretical guarantees for this setting of $\lambda$: the algorithm converges (Theorem 3.5), and the bias for samples truly sampled from the target does not increase with this regularization. Notice that Figure 6 in the Supplementary Material provides experiment results for lower and higher values of $\lambda$ than $1/m$: performance strongly decreases in both cases. Besides, the goal of Theorem 3.2 is to show that the regularized Stein thinning is now sensitive to the weights of multimodal distributions, as opposed to the original Stein thinning. However, we cannot theoretically determine which exact range of values of $\lambda$ lead to good thinned samples in finite sample regimes. As suggested by the reviewers, we will comment more on the limitation of this aspect of the L-KSD in the paper.

**IMQ kernel for Theorem 2.4.** The most efficient kernel for KSD-based algorithm is the IMQ kernel, as often stated in the literature, since using other types of kernels strongly degrades performance (Riabiz et al., 2022; Chen et al., 2018). Therefore, it was of utmost importance
that our theoretical results hold for the IMQ kernel used in practice. In our opinion, it was of secondary importance to extend the analysis to radial kernels. In the case of Theorem 2.3, it happens that the result holds for radial kernels without any additional assumption than for the IMQ kernel case. We thus decided to state the result in full generality. On the contrary, extending Theorem 2.4 to radial kernels requires intricate additional assumptions about the kernel properties. More importantly, it will strongly increase the complexity of the inequality, which relates the sample size m and the threshold s, leading to an obscure and unintuitive condition. For the sake of clarity and the practical scope of the theory, we thus believe that Theorem 2.4 should only be stated for IMQ kernels.

**Additional details in the theoretical analysis.** We will take advantage of the additional page, if the paper is accepted, to add details about assumptions of the theoretical results, to improve clarity. We separately answer to the points raised by the reviewers in dedicated posts.


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

Chen, W. Y., Mackey, L., Gorham, J., Briol, F. X., & Oates, C. (2018, July). Stein points. In International Conference on Machine Learning (pp. 844-853). PMLR.

Meta-review: All reviewers agreed that this work was above the acceptance threshold for NeurIPS.  The methodology is sound and supported by appropriate theoretical analysis.  I slightly question whether the specific methodological innovations were too focused on fixing specific examples where the original Stein thinning method performed poorly, and that in making these modifications to the algorithm the authors may have introduced a different collection of examples where the method performs poorly.  But nevertheless this is a very welcome contribution to NeurIPS!
Decision: Accept (poster)