Gradient-Free Kernel Stein Discrepancy
Authors: Matthew A Fisher
Chris J. Oates
Venue: NeurIPS 2023 poster (2023)
Keywords: Bayesian
discrepancy
kernel
sampling
Stein's method

Abstract: Stein discrepancies have emerged as a powerful statistical tool, being applied to fundamental statistical problems including parameter inference, goodness-of-fit testing, and sampling.  The canonical Stein discrepancies require the derivatives of a statistical model to be computed, and in return provide theoretical guarantees of convergence detection and control.  However, for complex statistical models, the stable numerical computation of derivatives can require bespoke algorithmic development and render Stein discrepancies impractical.  This paper focuses on posterior approximation using Stein discrepancies, and introduces a collection of non-canonical Stein discrepancies that are gradient-free, meaning that derivatives of the statistical model are not required.  Sufficient conditions for convergence detection and control are established, and applications to sampling and variational inference are presented.

Gradient-Free Kernel Stein Discrepancy
Matthew A. Fisher', Chris. Oates' 1.2
Newcastle University, UK Alan Turing Institute, UK
Abstract
Stein discrepancies have emerged as a powerful statistical tool, being applied to fundamental statistical problems including parameter inference, goodness-of-ff testing, and sampling. The canonical Stein discrepancies require the derivatives of a statistical model be computed, and return provide theoretical guarantees of convergence detection and control. However, for complex statistical models, the stable numerical computation of derivatives can require bespoke algorithmic development and render Stein discrepancies impractical. This paper focuses on posterior approximation using Stein discrepancies, and introduces a collection of non-canonical Stein discrepancies that are gradient-free, meaning that derivatives of the statistical model are not required. Sufficient conditions for convergence detection and control are established, and applications to sampling and variational inference are presented.
Introduction
Stein discrepancies were introduced in Gorham and Mackey [2015], way measure the quality of an empirical approximation to a continuous statistical model involving an intractable normalisa- tion constant. Rooted in Stein's 's method [Stein, 1972], the idea is consider empirical averages of a large collection test functions, each which is known to integrate to zero under the statisti- cal model. To date, test functions have been constructed by combining derivatives of the statistical model with reproducing kernels [Chwialkowski et al., 2016, Liu et al., 2016, Gorham and Mackey, 2017, Gong et al., 2021a,b], random features [Huggins and Mackey, 2018], diffusion coefficients and functions with bounded derivatives [Gorham et al., 2019], neural networks [Grathwohl et al., 2020], and polynomials [Chopin and Ducrocq, 2021]. The resulting discrepancies have been shown to be powerful statistical tools, with diverse applications including parameter inference [Barp et al., 2019, Matsubara et al., 2022], goodness-of-fe fit testing [Jitkrittum al., 2017, Fernandez et al., 2020], and sampling [Liu and Lee, 2017, Chen et al., 2018, 2019, Riabiz et al., 2022, Hodgkinson et al., 2020, Fisher al., 2021]. However, one of the main drawbacks of these existing works is the requirement that derivatives both exist and can be computed.
The use of non-differentiable statistical models is somewhat limited but includes, for example, Bayesian analyses where Laplace priors are used [Park and Casella, 2008, Ročková and George, 2018]. Much more common is the situation where derivatives exist but cannot easily be computed. In particular, for statistical models with parametric differential equations involved, one often re- quires different, more computationally intensive numerical methods to be used if the sensitivities (i.c. derivatives of the solution with respect to the parameters) are to be stably computed [Cockayne and Duncan, 2021]. For large-scale partial differential equation models, as used in finite element simulation, computation of sensitivities can increase simulation times by several orders of magni- tude, if is practical at all.
The motivation and focus of this paper is on computational methods for posterior approximation, and this end propose a collection of non-canonical Stein discrepancies that are gradient free, meaning that computation of the derivatives of the statistical model is not required. Gradient-free
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Stein operators were introduced in Han and Liu [2018] in the context of Stein variational gradient descent [Liu and Wang. 2016], but the theoretical properties of the corresponding discrepancy have yet to be investigated. General classes of Stein discrepancies were analysed in Huggins and Mackey [2018], Gorham et al. [2019]. but their main results do not cover the gradient-free Stein discrepancies developed in this work, for reasons that will be explained. The combination of gradient-free Stein operators and reproducing kernels studied in detail, to obtain discrepancies that can be explicitly computed. The usefulness of these discrepancies depends crucially on their ability to detect the convergence and non-convergence of sequences probability measures to the posterior target, and in both directions positive results are established.
Outline Gradient-free KSD (GF-KSD)) proposed and theoretically analysed in Section The proposed discrepancy involves certain degrees of freedom, including probability density denoted q in the sequel, and strategies for specifying these degrees of freedom are empirically assessed in Section 3. Two applications are then explored in detail; Stein importance sampling (Section 4.1) and Stein variational inference (Section 4.2). Conclusions are drawn in Section 5.
2 Methods
This section contains our core methodological (Section 2.1) and theoretical (Section 2.2) develop- ment. The following notation will be used:
Real Analytic Notation For a twice differentiable function f: B2 + R. let 0 denote the partial derivative of with respect to its ith argument, let f denote the gradient vector with entries a,f, and let V denote the Hessian matrix with entries 0, For a sufficiently regular bivariate function f 12d x _d + I. let (0, @ 0,)f indicate the application of 0 to the first argument of f, followed by the application of 0 to the second argument. (For derivatives of other orders, the same tensor notation @ will be used.)
Probabilistic Notation Let P(RR) denote the of probability distributions RRd. Let ((x) P(RR)) denote atomic distribution located at E Rdd. For ,, E P(RR), let < no indicate that Ir is absolutely continuous with respect to For π E P(R and (Tn)eeN C P(R"), write Th + æ to indicate weak convergence of m, to n. The symbols p and q are reserved for probability density functions on Rd, while is reserved for generic element of P(RR). For convenience, the symbols p and q will also be used to refer to the probability distributions that these densities represent.
2.1 Gradient-Free Kernel Stein Discrepancy
The aim of this section is to explain how GF-KSD can be constructed. Let p P(IR)) be target distribution of interest. Our starting point is a gradient-free Stein operator, introduced in Han and Liu [2018] in the context of Stein variational gradient descent [Liu and Wang, 2016]:
Definition (Gradient-Free Stein Operator). For p,q E P(Rd) with < p and V log well- defined, the gradient-fre free Stein operator defined as
Sp.r == 9 (V-h. h V logg),
acting on differentiable functions h : Rd - Rd
The Langevin Stein operator of Gorham and Mackey [2015] is recovered when p q. but when f p the dependence on the derivatives of is removed. The operator Sp.g can still be recognised as a diffusion Stein operator, being related the infinitesimal generator of a diffusion process that leaves invariant: however, it falls outside the scope of the theoretical analysis of Huggins and Mackey [2018]. Gorham al. [2019], for reasons explained in Remark 2. can also be viewed as a non-standard instance of the density method of Diaconis et al. [2004]; see Section 2 of Anastasiou et al. [2023]. The inclusion of q introduces an additional degree of freedom, specific choices for which are discussed in Section 3.
Remark 1. The ratio g/ p in Definition! could be viewed as an importance weight, but the construc- tion does not fully correspond to importance sampling due the V log q term, which is q q-dependent.
2

---

The Stein operator nomenclature derives from the vanishing integral property in Proposition be- low, which is central to Stein's 's method [Stein, 1972]:
Proposition 1. In the setting of Definition 1, assume that ||=]|<=) q(x) + as ||=||| + oo and [ all dq < oc. Then, for any function h : R + R whose first derivatives exist and are bounded, it holds that SSSp.gh dp = 0.
All proofs are contained in Appendix A. From Proposition 1, the expectation of Sp.gh with respect to π E ((Rd) will be zero when π and p are equal: conversely, the value of such an expectation can be used to quantify the extent to which π and p different. Consideration of multiple test functions increases the number and nature of the differences between и and that may be detected. A discrepancy is obtained by specifying which test functions h are considered, and then taking supremum over the expectations associated to this set. For computational convenience we take h to be contained in the unit ball of reproducing kernel Hilbert space, as described next.
For symmetric positive definite function Rd x Rdd -> RR called a kernel, denote the associated reproducing kernel Hilbert space H(k). Let H(k)) denote the Cartesian product of d copies of H(kk), equipped with the inner product (h. (h. == (h, H)))))))
Proposition 2. Let π E P(RT). In the setting of Definition 1, assume there is an a > 1 such that [(///)> dr < co and [ log a/[a//a-1) dx < cc. Let k: IRd x IR + R be a continuously dif- ferentiable kernel such that both is and its first derivatives x + (0; 0.)k(H, x). i 1,... d, are bounded. Then Sp.g is bounded linear operator from (k)) to L (TT)
For discrete distributions n, supported on a finite subset of Iad the moment conditions in Proposi- tion 2 are automatically satisfied. For general distributions π on RRd, the exponent a can be taken arbitrarily close to 1 to enable the more stringent moment condition J (q/p)² da < oo hold. An immediate consequence of Proposition 2 is that Definition 2 below is well-defined.
Definition 2 (Gradient-Pie Frredient-Free Kernel Stein Discrepancy). For p. q, k and æ satisfying the precondi- tions of Proposition 2, the gradient-free KSD is defined as
(p.g(T) suu { Spigh h du [[bl|x(() S $1).
(1)
The GF-KSD coincides with the canonical kernel Stein discrepancy (KSD) when p === q, and is thus strictly more general. Note that (T. is precisely the operator norm of the linear functional h + J S dx, which exists due to Proposition 2. Most common kemnlls satisfy the assumptions of Proposition 2, and a particularly important example is the inverse multi-quadric kernel
k(z, v) (a² + || - -
00), , E (0, 1),
which has bounded derivatives of all orders; see Lemma 4 of Fisher et al. [2021].
The use of reproducing kernels ensures that GF-KSD can be explicitly computed; see Proposi- tion 5 and Corollary 1 in Appendix A. A general spectral characterisation of GF-KSD is provided in Proposition of the supplement, inspired by the recent work of Wynne et al. [2022]. Note that the scale of Equation (6) does not matter when one is interested in the relative performance of dif- ferent π E P(R as approximations of a fixed target p E P(R.).. In this sense, GF-KSD may be employed with ğ in place of where p x p/Z Z and 2 is a normalisation constant. This feature makes GF-KSD applicable to problems of posterior approximation, and will be exploited for both sampling and variational inference in Section 4. On the other hand, GF-KSD is not applicable to problems in which the target distribution py involves a parameter 0, such as estimation and composite hypothesis testing, since then the normalisation term Zp cannot be treated as constant.
Going beyond the discrepancies discussed in Section several non-canonical discrepancies have recently been proposed based on Stein's 's method, for example sliced [Gong et al., 2,,1a, stochas- tie [Gorham et al., 2020], and conditional [Singhal et al., 2019] Stein discrepancies (in all cases derivatives of p are required). However, only a subset of these discrepancies have been shown to enjoy important guarantees of convergence detection and control. Convergence control is critical for the posterior approximation task considered in this work, since this guarantees that minimisation of Stein discrepancy will produce consistent approximation of the posterior target. The aim of the next section is to establish such guarantees for GF-KSD.
3

---

2.2 Convergence Detection and Control
The canonical KSD benefits from theoretical guarantees of convergence detection and control [Gorham and Mackey, 2017]. The aim of this section is to establish analogous guarantees for GF-KSD. To set the scene, for a Lipschitz function f R2d + IR& denote its Lipschitz constant L(f) := sup (z) - f())) | yll. Then, for measurable g - IR, we denote the tilted Wasserstein distance as ( p) == sup supLifist f fg dv - [ fg dpl whenever this expression is well-defined [Huggins and Mackey, 2018]. Note that the standard I-Wasserstein distance W WI(T,, p) is recovered when g 1. There is no dominance relation between Wi(+::9) for different g: the topologies they induce are different. Tilted Wasserstein distance induces a much weaker topology than, for example, divergences such as Kullback-Leibee or Hellinger, since it does not require ab- solute continuity of measures.
Theorem 1 (Convergence Detection). Let p.q P(Rd) with & p. V Vlogg q Lipschitz and V log q/|² dg < oc. Assume there an a > 1 such that the sequence (Tn) C((T..) sat- E (0,) || V glla/(a-1)
< 00, [[ log q|| q|]///a-1] (g/p) dr < oc, and I fq/p d 0o with f(x)= = ||=]|, for each n E N. Let i be a kernel such that each of k. (0, @ 0.)k and (0,0, @ 0, ) exist, are continuous, and are bounded, for i,j E (1.....d). Then W (Ta,
isfies ()
implies D (Tn
Thus the convergence of n. to p. in the sense of the tilted Wasserstein distance with g - q p/p, is detected by the GF-KSD D Dp.g. Note that the conditions on (Ta)eeN in Theorem are automatically satisfied when each Ta, has finite support.
Despite being a natural generalisation of KSD, this GF-KSD does not in general provide weak convergence control in the equivalent theoretical context. Indeed, Gorham and Mackey [2017] es- tablished positive results on convergence control for distributions in Q(E4), the set of probability distributions on Rd with positive density function g R + (0, oo) for which V log q is Lipschitz and y is distantly dissipative, meaning that
(V log g(x) - V log (()) v) ||:0 - yll = r > 0.
in inf P-
1x y/
Under the equivalent assumptions, convergence control fails for GF-KSD in general:
Proposition 3 (Convergence Control Fails in General). Let k be a (non-identically zero) radial kernel of the form k(æ, y) o(x y) for some twice continuously differentiable ó ; Rd + R. for which the preconditions of Proposition 2 are satisfied. Then there exist p E P(R"),, Q(R)) and a sequence (TR)nEN _ (((R"), also satisfying the preconditions of Proposition 2, such that Dp.((TT) - 0 and yet Th ++
The purpose of Proposition
to highlight that GF-KSD is not a trivial extension of canonical KSD; it requires bespoke treatment. This is provided in Theorem 2, next. Indeed, to ensure that GF-KSD provides convergence control, additional condition on the tails of q are required:
Theorem 2 (Convergence Control). Let p E P(RR), q E Q(Rd) be such that is continuous and inf q(x)/p(x)> > O. Assume there is an > 1 such that the sequence (T (nnnEN C ((RR) sat- dm < so and f | log ql|" (=/(a=1)(//)) dā,, < cc, for each n3 E N. Let k be the inverse mullti-quadric kernel in Equation (2). Then Dp. (T.) - 0 implies nn -
isfies (g/p)) dm,
The proof of Theorem 2 is based on carefully re-casting GF-KSD between π and p canonical KSD between q and transformed distribution मे (see Proposition 7 in the supplement), then appeal- ing the analysis of Gorham and Mackey [2017]. Compared to the analysis of KSD in Gorham and Mackey [2017], the distant dissipativity condition now appears on & (a degree of freedom), rather than on p distribution determined by the task at hand), offering realistic opportunity for this condition to be verified. For example, one could take a Gaussian measure that dominates the target distribution p for in GF-KSD. Nevertheless, the conditions of Theorem rule out distributions p that are heavy-tailed. Suitable choices for q are considered in Section 3.
Remark 2 (Related Work). Convergence control was estublished for discrepancies based on general classes of Stein operator earlier work, but the required assumptions are too stringent when applied in our context. In particular, to use Huggins and Mackey [2018] it is required that the gradient

---

Convergence
Converpence
Prior
Laplace
lop GI-KSD
16 P.LL HTPT
Inward Convergence
Inward Non-Convergence
GMM
9= KDB
Oblique Convergence
Oblique Non- Convergence
lng GP&SD
100
(a)
(b)
Figure Empirical assessment of gradient-free kernel Stein discrepancy. (a) Test sequences (Ta) defined in Appendix C.I. The first column displays sequences (solid) that converge to the distributional target p (black), while the second column displays sequences (dashed) which con- verge instead to fixed Gaussian target. (b) Performance of gradient-free kernel Stein discrepancy, when approaches to selecting q (described the main text) are employed. The colour and style of each curve in (b) indicates which of the sequences in (a) is being considered. [Here we fixed the kernel parameters a and 1/2.]
V llg(g/p) is bounded', while use Gorham al. [2019] is required that q/p is Lipschit??. In our context, where must be specified in ignorance of p. such conditions, which require that is almost as light as p the tail, cannot be guaranteed to hold. The present paper therefore instead contributes novel analysis for the regime where q may be appreciably heavier than in the tail.
This completes our theoretical assessment, but the practical performance GF-KSD remains to be assessed. Suitable choices for both q and the kernel parameters and 8 are proposed and investi- gated Section 3, and practical demonstrations of GF-KSD are presented in Section 4.
Implementation Detail
The purpose of this section is to empirically explore the effect of varying q. a and B, aiming to arrive at reasonable default settings. In the absence of an application-speciñe optimality criterion, we aim to select values that perform well (in a sense to specified) over range of scenarios that may be encountered. Here, assess performance several sequences (m, hee are considered, some of which converge to specified limit p and the rest of which converge to an alternative Gaussian target; see Figure (a)) An effective discrepancy should clearly indicate which of these sequences are convergent and which are not. On this basis, recommendations for are considered in Section 3.1, and recommendations for and b in Section 3.2. Of course, cannot expect default settings to perform universally well, so in Section 3.3 we highlight scenarios where our defaults may fail. Python code reproduce the experiments reported below can be downloaded at [blinded].
3.1
Choice of
In what follows we cast q as an approximation of p. aiming to inherit the desirable performance of canonical KSD for which and p are equal. The task to which the discrepancy is being applied will,
'oo see this, take A = q/p in Theorem 3.2 of Huggins and Mackey [2018].
To see this, take (9/p)/ in Theorem and Proposition of Gorham al. [2019].

---

in practice, constrain the nature and form of the distributions q that can be implemented. For expos-
itory purposes (only), the following qualitatively distinct approaches choosing q are considered:
Prior In Bayesian settings where p is posterior distribution, selecting q to be the prior distri- bution ensures that the condition gef ((z)/p(x) > is satisfied.
Laplace If the target p can be differentiated, albeit a possibly high computational cost, it may be practical to construct a Laplace approximation to the target [Gelman et al., 2013].
GMM One could take q to be Gaussian mixture model fitted to approximate samples from p. representing more flexible alternative to Laplace.
KDE Additional flexibility can be obtained by employing a kernel density estimator as a non- parametric alternative to GMM.
Of course, there is circularity to GMM and KDE which renders these methods impractical in general. The specific details of how each of the q were constructed are contained in the code that accompanies this paper, but the resulting are displayed as insets in Figure 1(b). The performance of GF-KSD with these different choices of q is also displayed in Figure 1(b). It was observed that all four choices of produced discrepancy that could detect convergence of Tn to p. though the detection of convergence was less clear for Prior due to slower convergence of the discrepancy 0 as n was increased. On the other hand, all approaches were able to clearly detect non-convergence to the target. That Laplace performed comparably with GMM and KDE was surprising, given that the target is not -ell-approximated by single Gaussian component. These results are for a specific choice of target p. but in Appendix C.5.1 a range of p are considered and similar conclusions are obtained. Section 3.3 explores how challenging p must be before the convergence detection and control properties associated to Laplace fail.
3.2 Choice of and |3
For the investigation in Section 3.1 the parameters of the inverse multi-quadric kernel (2) were fixed to a = 1 and 3 1/2, the latter being the midpoint of the permitted range 8 E (0, 1). In general, care in the selection of these parameters may required. The parameter a captures the scale of the data, and thus standardisation of the data may be employed to arrive at o = 1 as a natural default. In this paper (with the exception of Section 4.2) the standardisation I + C- '2 was performed, where C' is the covariance matrix of the approximating distribution g being used. In Appendix C.5.2 we reproduce the investigation of Section 3.1 using range of values for and |3: these results indicate the performance of GF-KSD is remarkably insensitive to perturbations around (a, () = (1, 1/2).
3.3 Avoidance of Failure Modes
GF-KSD is not a silver bullet, and there are a number of specific failure modes that care may be required to avoid. The four main failure modes are illustrated in Figure 2. These are as follows: (a) is substantially heavier than in a tail; (b) q is substantially lighter than p in a tail; (c) the dimension is too high; (d) p has well-separated high-probability regions. Under both (a) and (c), convergence detection can fail, either because theoretical conditions are violated or because the terms म, must be extremely close to p before convergence begins be detected, Under (b), the values of GF-KSD at small n can mislead. Point (d) is a well-known pathology of all score- methods; see Wenliang and Kanagawa [2021]. Liu et al. [2023]. These four failure modes inform our recommended usage of GF-KSD, summarised next.
Summary of Recommendations Based on the investigation just reported, the default settings we recommend are Laplace with = (post-standardisation)) and / = 1/2. Although not uni- versally applicable, Laplace does not require samples from p. and has no settings that must be user-specified. Thus we recommend the use of Laplace in situations where a Laplace approxima- tion can be justified and computed. If Laplace not applicable, then one may attempt to construct an approximation q using techniques available for the task at hand (e.g. in a Bayesian setting, one may obtain q via variational inference, or via inference hased on an approximate likelihood). These recommended settings will be used for the application presented in Section 4.1, next.

---

q p
d
P=p
log GFKSD
log. GFKSD
log GIRSD
log OFKSD
de
120
190
18
(a)
(h)
(d)
Figure 2: Failure modes: (a) is substantially heavier than p in tail; (b) q is substantially lighter than in a tail; (c) the dimension d is too high; (d) p has separated high-probability regions. [For each of (a), (b) and (d), the colour and style of the curves refers the mon--onvergence (outvardinwardlobiqqee presented in Figure 1, and we plot the logarithm of the gradient-free kernel Stein discrepancy function of the index n of the sequence (Tn)nE.. For (c)
sense of convergence
consider convergent sequences (Mn)EE of distributions IRd
Applications
To demonstrate potential uses of GF-KSD, two applications posterior approximation presented: Stein importance sampling (Section 4.1) and Stein variational inference using measure transport (Section 4.2). In each case, extend the applicability of existing algorithms to statistical models for which certain derivatives of p are either expensive or non-exisie non-existent.
4.1 Gradient-Free Stein Importance Sampling
Stein importance sampling [Liu and Lee, 2017, Hodgkinson et al., 2020] operates by first sampling independently from a tractable approximation of the target and then correcting the bias the sam- ples so-obtained. To date, applications of Stein importance sampling have been limited to instances where the statistical model p can be differentiated; our contribution to remove this requirement. In what follows analyse Stein importance sampling in which independent samples (z,, (Zn))eN generated from the same approximating distribution q that is employed within GF-KSD:
Theorem (Gradient-Free Stein Importance Sampling). Let p P(R q Q(R:) be such that is continuous and inf q(x)/p(a) O. Suppose that cxp(y/l V log all?) dq < oo for some y Let be the inverse multi quilii-quaric kernel Equation (2). Let (In)neN be independent samplex from To the sumple, assign optimal weights
w² E arg min { (D Dp.a Eurd(a:)) :0 WT Wny wi + + Wn = =1
Then Tn ;= EE wf 8(x,) satisfies nn 4 almost surely as n + cc.
The proof builds on earlier work in Riabiz et al. [2022]. Note that the optimal weights w can be computed without the normalisation constant of p. by solving a constrained quadratic programme at cost (
As an illustration, we implemented gradient-free Stein importance sampling to approximate a pos- terior arising from discretely observed
model
i(t) = au(() Bu(t)o(t),
i(t) -yv(t) + Su(t)u(t),
(n(0), u(0)) (ug, vo),
with independent log- normal observations with covariance matrix diag(of, oz). The parameters to be inferred are (a. ,, 5, up Up (..0.1.) σ2) and therefore = The data analysed are due Hewitt [1921], and full details are contained in Appendix C.3. The direct application of Stein importance sampling this task requires the numerical calculation of sensitivities of the differential equation at each of the samples that are to be re-weighted. Aside from simple cases where automatic differentiation or adjoint methods can be used, the stable computation sensitivities can form major computational bottleneck; see [Riabiz et al., 2022]. In contrast, approach required

---

GF-KSD Self Nomalised
loga
(h)
Figure 3: Gradient-Free Stein Importance Sampling: (a) The lower triangular panels display n = 20 independent (biased) samples from the Laplace approximation q. while the upper triangular panels display the same number of re-weighted samples obtained using gradient-free free Stein importance sampling. [Samples are shown in blue, with their size proportional to the square of their weight, to aid visualisation. The shaded background indicates the high probability regions of p. the target.] ] (b) The approximation quality, the energy distance between the approximation and the target. [Thh solid line corresponds to the output of gradient-free Stein importance sampling, while the dashed line corresponds to the output of nelf-normalised importance sampling. Standard error regions are shaded.]
function of the number of samples from q. is measured as
fixed number of gradient computations to construct a Laplace approximation', independent of the number of samples required; Appendix C.3 for detail. In the lower triangular portion of Figure 3()), biased samples from the Laplace approximation 7 (# p) are displayed, while in the upper triangular portion the same samples are re-weighted using gradient-free Stein importance sampling to form a consistent approximation of p. A visual reduction in bias and improvement in approximation quality can be observed. As baseline against which to assess the quality of our approximation, we consider norm-nord importance sampling; i.e. the approximation with weights such that W; p(xx)//(x:; Figure 3(b) reports the accuracy of the approximations p as quantified using energy distance [Cramér, 1928]. These results indicate that the approximations produced using gradient-free Stein importance sampling improve on those constructed using self- normalised importance sampling. This may be explained by the fact that the optimal weights w/ attempt to mitigate both bias due to g * p and Monte Carlo error, while the weights w only address the bias due to # and do not attempt Additional experiments in Appendix C.5.3 confirm that gradient-free Stein importance sampling achieves comparable performance with gradient-based Stein importance sampling in regimes where the Laplace approximation can justified.
mitigate error due to the randomness in Monte Carlo.
Although our recommended default settings for KSD were successful in this example, an inter- esting theoretical question would be to characterise an optimal choice of q this context. This ap- pears to be a challenging problem but we hope to address it in future work. In addition, although we focused on Stein importance sampling, our methodology offers the possibility to construct gradient- free versions of other related algorithms, including the Stein points algorithms of Chen et al. [2018, 2019], and the Stein thinning algorithm of Riabiz et al. [2022].
4.2 Stein Variational Inference Without Second-Order Gradient
Stein discrepancy proposed variational objective in Ranganath et al. [2016] and has demon- strated comparable performance with the traditional Kullback-Leibler objective in certain applica- tion areas, whilst abolishing the requirement that the variational family is absolutely continuous with
In this case, 49 first order derivatives were computed, of which 48 were on the optimisation path and was used to construct finite difference approximation to the Hessian.

---

GP-KSD O denntions
GF-KSD 2000 decations
GF KSD 10,000 Berritionss
GF-KSD 20,00 inceations
KSD 20,00 dentiions
Figure 4: Stein Variational Inference Without Second Order Gradient: The top row concerns ap- proximation of a distributional target that is "banana" shaped, while the bottom row concerns a "sinusoidal" target. The first four columns depict the variational approximation Te,., constructed using gradient descent applied to gradient-free kernel Stein discrepancy (i.e. first order derivatives of y required) along the stochastic optimisation sample path (m E {0,2 2 10S, 10 2 x 10*)), while the final column reports the corresponding approximation (m = 2 x 10*) constructed using stan- dard kernel Stein discrepancy (i.e. second order derivatives of p required).
respect the statistical model [Fisher al., 2021]. This offers an exciting, as yet largely unexplored opportunity to construct flexible variational families outside the conventional setting of normalising flows (which are constrained be diffeomorphisms of IRd). However, gradient-basee stochastic op- timisation of the canonical Stein discrepancy objective function means that second-order derivatives of the statistical model are required. GF-KSD reduces the order of derivatives that are required from second-order to first-order, and in many cases (such as when differential equations appear in the statistical model) this will correspond to a considerable reduction in computational cost.
In what follows we fix a reference distribution R E P(IR)) and a parametric class of maps To Rd Rdd 0 E IRP, and consider the variational family (To) P(Rd) whose elements To = T R are the pushforwards of R through the maps 7 0 E IRP. The aim is to use stochas- tic optimisation to minimise 0 + (Te) (or, equivalently, any strictly increasing transformation thereof). For this purpose a low-cost unbiased estimate of the gradient of the objective function is required: the details and sufficient theoretical conditions are contained in Appendix B.
As an interesting methodological extension, that departs from the usual setting of stochastic op- timisation, here consider interlacing stochastic optimisation over 0 and the selection of q, leveraging the current value Am on the optimisation path to provide a natural candidate Te... for q in this context'. Thus, for example, 0m+1- = On VoDp. (TT)² (To)²|e=__) le=0., for some learning rate > O. (In this iterative setting, to ensure the variational objective remains fixed, do not perform the standardisation of the data described in Section 3.2.)
vanilla stochastic gradient descent routine becomes
To assess the performance of GF-KSD in this context, re-instantiated an experiment from Fisher et al. [2021]. The results, in Figure 4, concern the approximation of "banana" and "sinusoidal" dis- tributions in dimension d 2, and were obtained using the reference distribution R = (0, and taking T to be the inverse autoregressive flow of Kingma et al. [2016]. These are both toy prob- lems, which do not themselves motivate our methodological development, but do enable us to have an explicit ground truth to benchmark performance against. Full experimental detail is contained in Appendix C.4; we highlight that gradient clipping was used, both to avoid extreme values of q/p encountered on the optimisation path, and to accelerate the optimisation itself [Zhang et al., 2019]. The rightmost panel depicts the result of performing variational inference with the standard KSD objective functional. It is interesting to observe that GF-KSD leads to similar performance in both examples, with the caveat that stochastic optimisation was more prone to occasional failure when GF-KSD was used. The development of a robust optimisation technique in this context requires care and a detailed empirical assessment, and is left as a promising avenue for further research.
*If the density of te,, is not available, for example because T is a complicated mapping, it can be consis- tently estimated using independent samples from 7

---

Conclusion
In this paper GF-KSD was proposed and studied. Theoretical and empirical results support the use of GF-KSD in settings where initial approximation to the distributional target can readily be constructed, and where the distributional target itself does not contain distant high probability regions, but poor performance can occur outside this context. Nevertheless, for many statistical analyses the principal challenge the cost evaluating the statistical model and its derivatives, rather than the complexity of the target itself, and in these settings the proposed diserepancy has the potential to be usefully employed. The focus of this work was on posterior approximation, with illustrative applications sampling and variational inference being presented. A natural extension to this work would involve a systematic empirical assessment of the performance of GF-KSD across broad range of applied contexts. However, we note that GF-KSD is applicable to problems such as estimation and composite hypothesis testing, where the target ranges over a parametric model class, and alternative strategies will be required to circumvent gradient computation in that context.
Acknowledgements MAF was supported by EP/W52238711 CJO was supported by EP/W019590/1. The authors grateful to François-Xave Xavier Briol, Jon Cockayne, Jeremias Knoblauch, Lester Mackey, Marina Riabiz, Rob Salomone, Leah South, and George Wynne for insightful comments early draft of the manuscript, as well as the Reviewers NeurlPS.
References
A. Anastasiou, A. Barp, F.-. Briol, Ebner, R. E. Gaunt, F. Ghaderinezhad, Gorham, A. Gret- ton, C. Ley, Q. Liu, L. Mackey, C. J. Oates, G. Reinert, and Y. Swan. Stein's method meets statistics: A review of some recent developments. Statistical Science, 38(1): 2023.
A. Barp, F.-X Briol, A. Duncan, M. Girolami, and L. Mackey. Minimum Stein discrepancy timators. In Proceedings of the 33rd Conference on Neural Information Processing Systems, volume 32, pages 12964-12978 12976, 2019.
W. Y. Chen, L. Mackey, J. Gorham, F.-X, Briol, and C. J. Oates. Stein points. In Proceedings of the 35th 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. J. Oates. Stein point Markov chain Monte Carlo. In Proceedings of the 36th International Conference on Machine Learning, pages 1011-1021. PMLR, 2019.
N. Chopin and G. Ducrocq. Fast compression of MCMC output. Entropy, 23(8):1017, 2021.
K. Chwialkowski, H. Strathmann, and A. Gretton. A kernel test of goodness of fit. In Proceedings of the 33rd International Conference on Machine Learning, pages 2606-2615. PMLR, 2016.
Cockayne and A. Duncan. Probabilistic gradients for fast calibration of differential equation models. SIAM/ASA Journal on Uncertainty Quantification, 9(4):1643-1673 2021.
H. Cramér. On the composition of elementary errors. Scandinavian Actuarial Journal, (1):141-180, 1928.
Diaconis, C. Stein, Holmes, and G. Reinert. Use of exchangeable pairs the analysis of simulations. In Stein's 's Method, pages 1-25. Institute of Mathematical Statisties, 2004.
T. Fernandez, N. Rivera, W. Xu, and A. Gretton. Kernelized Stein discrepancy tests of goodness- of-fit for time-to-event data. In Proceedings of the 37th International Conference on Machine Learning. pages 3112-3122, PMLR, 2020.
M. Fisher, T. Nolan, M. Graham, D. Prangle, and C. J. Oates. Measure transport with kernel Stein discrepancy. In Proceedings of the 24th International Conference Artificial Intelligence and Statistics, pages 1054-1062. PMLR, 2021. (Here we refer to the error-corrected version arXiv.2010.11779),
10

---

A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin. Bayesian Data Analysis. CRC Press, 2013.
W. Gong, Y. Li, and J. M. Hernández-Lobato. Sliced kemeliied Stein discrepancy. In Proceedings of the 9th International Conference Learning Representations, 2021a.
W. Gong, K. Zhang, Y. Li, and J. M. Bernández-Lobato Active slices for sliced Stein discrepancy. In Proceedings of the 38th International Conference on Machine Learning, pages 3766-3776, 2021b.
Gorham and L.. Mackey. Measuring sample quality with Stein's method. In Proceedings of the 29th Conference Neural Information Processing Systems, pages 226-234, 2015.
Gorham and Mackey, Measuring sample quality with kernels. In Proceedings of the 34th International Conference Machine Learning, pages 1292-1301. PMLR, 2017.
J. Gorham, A. B. Duncan, S. Vollmer, and L. Mackey, Measuring sample quality with diffusions. The Annals of Applied Probability, 29(5):2884-29) 2928, 2019.
J. Gorham, A. Raj, and L. Mackey. Stochastic Stein discrepancies. In Proceedings of the 34th Conference on Neural Information Processing Systems, 2020.
W. Grathwohl, K.-C Wang. J.-H. Jacobsen, D. Duvenaud, and R. Zemel. Learning the Stein dis- crepancy for training and evaluating mergy-baser models without sampling. In Proceedings of the 37th International Conference on Machine Learning, pages 3732-3747. PMLR, 2020.
J. Han and Q. Liu, Stein variational gradient descent without gradient. In Proceedings of the 35th International Conference on Machine Learning, pages 1900-1908. PMLR, 2018.
C. G. Hewitt. The Conservation of the Wild Life of Canada. Charles Scribner's Sons, 1921.
L. Hodgkinson, R. Salomone, and F. Roosta. The reproducing Stein kernel approach for post-hoc corrected sampling. arXiv:2001.091 092666 2020.
H. Huggins and L.. Mackey. Random feature Stein discrepancies. In Proceedings of the 32nd International Conference Neural Information Processing Systems, pages 1903-1913, 2018.
W. Jitkrittum, W. Xu, Z. Szabó, K. Fukumizu, and A. Gretton. A linear-time kernel goodness-of-fi test. In Proceedings of the 31st Conference on Neural Information Processing Systems. NIPS Foundation, 2017.
D. Kingma and Ba. Adam: A method for stochastic optimization. Proceedings of the 3rd International Conference Learning Representations, 2015.
D. Kingma, T. Salimans, R. Jozefowicz,
X. Chen, Sutskever, and M. Welling. Improved vari- ational inference with inverse autoregressive flow. In Proceedings of the 30th Conference on Neural Information Processing Systems, 2016.
D. C. Liu and Nocedal. On the limited memory BFGS method for large scale optimization. Mathematical Programming, 45(3):503-528, 1989.
Q Liu and J. Lee. Black-box importance sampling. In Proceedings of the 20th International Con- ference on Artificial Intelligence and Statistics, pages 952-961. PMLR, 2017.
Q. Liu and D. Wang. Stein variational gradient descent: A general purpose Bayesian inference algorithm. Proceedings of the 30th Conference on Neural Information Processing Systems, 2016.
Q. Liu, J. Lee, and M. Jordan. A kemelized Stein discrepancy for goodness-of-frt tests. In Proceed- ings of the 33rd International Conference on Machine Learning, pages 276-284, PMLR, 2016.
X. Liu, A. Duncan, and A. Gandy. Using perturbation to improve goodness-of-fit tests based on kernelized Stein discrepancy. In Proceedings of the 40th International Conference on Machine Learning, 2023.
11

---

T. Matsubara, J. Knoblauch, F.-X. Briol, and C. J. Oates. Robust generalised Bayesian inference for intractable likelihoods. Journal of the Royal Statistical Society, Series B, 84(3):997-0022, 2022.
B. O' O'Donoghue, E. Chu, N. Parikh, and S. Boyd. Conic optimization via operator splitting and homogeneous self-dual embedding. Journal of Oplimization Theory and Applications, 169(3); 1042-1068, June 2016.
Park and G. Casella. The Bayesian lasso. Journal of the American Statistical Association, 103 (482):681-688) 2008.
D. Prangle and C. Viscardi. Distilling importance sampling for likelihood free inference. Journal of Computational and Graphical Statistics, 2023. To appear.
R. Ranganath, D. Tran, J. Altosaar, and D. Blei. Operator variational inference. In Proceedings of the 30th Conference on Neural Information Processing Systems, pages 496-504, 2016.
M. Riabiz, W. Chen, Cockayne, P. Swietach, A. Niederer, L. Mackey, and Oates. Optimal thinning of MCMC output. Journal of the Royal Statistical Sociery, Series B, 1404),1059-1081 2022.
V. Ročková and E. I. George. The Spike-and- LASSO. Journal of the American Statistical Association, 1136521):411444,, 2018.
G. Schwarz. Estimating the Dimension of a Model. The Annals of Statistics, 6(2):461 -664, 1978.
464,
B. W. Silverman. Density Estimation for Statistics and Data Analysis. Chapman & Hall, 1986.
R. Singhal, X. Han, Lahlou, and R. Ranganath. Kernelized complete conditional Stein discrep- ancy. J00::9404.4478 2019.
Stan Development Team. RStan: the R interface to Stan, 2022. URL httppss///mcsstaa.orr/// R package version 2.21.5.
C. Stein. bound for the error in the normal approximation to the distribution of sum of dependent random variables. In Proceedings of the 6th Berkeley Symposium on Mathematical Statistics and Probability, volume 6.2, pages 583-603. University of California Press, 1972.
I. Steinwart and A. Christmann. Support Vector Machines. Springer Science & Business Media, 2008.
K. Wenliang and H. Kanagawa, Blindness score-based methods to isolated components and mixing proportions. In Proceedings of the NeurlPS workshop "Your Model is Wrong", 2021.
G. Wynne, M. Kasprzak, and A. B. Duncan. A spectral representation of kernel Stein discrep- ancy with application to goodness-offff tests for measures on infinite dimensional Hilbert spaces. arXiv:2206.04552 2022.
J. Zhang, T. He, S. Sra, and A. Jadbabaie. Why gradient clipping accelerates training: A theoretical justification for adaptivity. In International Conference on Learning Representations, 2019.
12

---

Appendices
These appendices supplement the paper Gradient Free Kernel Stein Discrepancy. The proofs for the- oretical results stated in the main text are contained in Appendix A. Theoretical analysis of stochastic gradient estimators for use in the context of Stein Variational Inference contained in Appendix B, as advertised in Section 4.2 of the main text. All additional details related to empirical assessment are contained in Appendix C.
A Proof of Results in the Main Text
This appendix contains proofs for all novel theoretical results reported in the main text.
Proof of Proposition 1. First we show that the integral in the statement is well-defined. Since and its first derivatives bounded, we can define the constants Co sup(])((x]] x R and
C := sup h(z)l= R.. Then
I hp dp /
? V. h h- Vloggl dp = dp=
V h h Vlog dq
S + Co / |V qll dq < co,
so indeed the integral is well-defined. Now, let = fz Rd (a r), and t(r) sup (g(x) x E S,), so that by assumption -d 0 as oc. Let 1 if x E B and 0 $ B. Then J Sp dp- = lim, J 1 B, Sp.gh dp and, from the divergence theorem,
f S1 dp =
9 V h + h- v loggl dp
1
p
/ 1 - Vgl dr
B,
/ v (gh) dx
n,
f qh-n dr st(r) x r(d/2)
4/2
r 1 200 0,
S.
as required.
Proof of Proposition 2. Since ls and its first derivatives are bounded, can set
d Σ (a, @ 0,)
c/ = sup / x) aERd
C == sup acRd
From Cauchy-Schwan and the reproducing property, for any E H(kk it holds that [f(x)) =)) S [lfll=(s)ll(((( [|ll / /k(x,2). x). Furthermore, using the fact that & is contimuously differentiable, it can be shown that (x)) < 11/lbaa ((a, 0z)k((zr æ); see Corollary 4.36 of Steinwart and Christmann [2008]. As a consequence, for all h e H(AA)d d we have that, for all E Rdd
(A(, (k(,,)),,,)( x))(()
[h(z)]| h((1) S
k(x,z) 1he c/ [|]])
V --((x)) =
Σ Dz, 1z(2))
Σ (a, 2 00) S cfllllaae
13

---

To use analogous notation as in the proof of Proposition 1, set Co := Co | and C1 ==
Cf [blluaaa Then, using Hölder's inequality and the fact that (a + b)8 v 20 (a + b with
3 == aaaaaa- 1), we have
|Sp.g/llesss>> = / (ppigh do
o V h h. Vlogal dı
! dv
(V. h h.
log () - da
+ +00 / | == dx
< 2-
dæ
Ci
(ʃ
a
2=1
(c+))=++(0)) +Ct))=²
/ /IV log qll log
$ 2 [
dm
- dr
as required.
To prove obtain explicit computable formulae for the GF-KSD, two intermediate results are required: Proposition 4. Let k and π satisfy the preconditions of Proposition 2. Then the function
Rdd 2 2 f(x) = g(x))
k + k( Voo
log q(x)]
(3)
p(x)
takes values in H((k),, is Bochner r-inttgrrable and, thus, & : f dr E H()k)
Proof. Since k has continuous first derivatives + (0, 0,) Lemma 4.34 of Steinwart and Christmann [2008] gives that (a, @ 1)k(2.)) H((), and, thus f E H(k)". Furthermore, f R H(k)d is Bochner r-integrable since
f (xlaa) dx(x)) =
((x))
V V.k(x,)) + k(x, . log (z)) (2)]
p(z)
V dx (ʃ [IVaA(e + (( Vog (ogg p
( بير dr(x)
g dm 1
2
(ct +C+)) / | log a
dn
00
where we have employed the same Cf and C+ notation as used in the proof of Proposition 2. Thus, from the definition of the Bochner integral, & = dr exists and is an element of H(k)d.
Proposition 5. Let k and n satisfy the preconditions of Proposition Then
D (T)( = //
((x) q(y),
p(x) p(y) k kg v) dT(x)dr(v)
(4)
where
(g, = V. Vyk(I, v) + (Z, w). V, logg()))) + (V V,, v), Ve log (x))
+ klz 3) (V (I)) V, (ogg
Proof. Let ff be as in Equation (3). From Proposition 4, f dm H(k) Moreover, since f is Bochner T-integrable and Tf (h, H' continuous linear functional on H(k) from basic properties of Bochner integrals we have T T f dm = TS dr. In partieular,
((x))
(h, H(K)
h,
[uk(r,)) + k( V log Vlg ((x) dr(x))
((x)
H(k)
(xx)
V k(z, ))00 (h, k(2, 1) M(aj V log q(z)] dr(x)
p(x)
q(x)
V Vh(x)) + h(x) v log [(x) dx(x)= = / Spalt dr(z)
p(x)
14

---

which shows that & is the Riesz representer of the bounded linear functional h + f SS..h dæ on H(k) It follows from Cauchy-Schwarz
Schwarz that the (squared) operator norm of this functional is
( g(xr)
(T) = ( = (5, (E,E))(R))
Vak(a, -) + k(x, V Voga(()) dx(x).
p(æ)
/
((U) Vyk(y. k(y, Vlog ((v)l dx(x)
p(y)
H(k)
// ((x) q(y)
ka y) dz(z))a())
p()) p(y)
as claimed.
For concreteness, we instantiate Proposition 5 in the specific case of the inverse multi-quadric kernel, since this is the kernel that we recommend in Section 2.2:
Corollary (Explicit Form). For p. q. and satisfying the preconditions of Proposition and the inverse multi-quadric kernel in Equation (2), we have that
(T. = =1
((x))())
43(3 1)1l0 - y11²
d + (V log q(x)) log q(y), 2 y)
p(z)p(v) (o + o vl 2| + 2.9
(o2 + 20 211
V log ((z), V log {(v))
dx(r)dr(v)
(6)
(a² + æ y 2)3
Proof of Corollary First we compute derivatives of the kernel & in Equation (2):
233
V vak(,,
(x v)
(o² + || w112)8+1
2
V (,,(,,
(x y)
(o2 + |r 2
8+1
48(3 + 1) œ y
20d
V vyk(z,y) (o + || -11)~ 842
(a² + | v112)*+1
Letting uxx)= == V logg(2) form a convenient shorthand, we have that
(ig u) V vy v) + (V 2=, v). u(y)) + (v, u(x)) + k(,, v) (u(xx,,
4. 1) 8//²
d (u(x) - - u(v), 9) (o2 || y 12
(u(z), u(y))
(o2 + || 8|2 2/3
(o2 + || y
which, combined with Proposition 5, gives the result.
In addition to the results in the main text, here we present a spectral characterisation of GF-KSD. The following result was inspired by impressive recent contribution to the literature on kernel Stein discrepancy due to Wynne et al. [2022], and our (informal) proof is based on an essentially identical argument:
Proposition 6 (Spectral Characterisation). Consider positive definite isotropic kernel k. and recall that Bochner's theorem guarantees k(x, y) = f -(s,, dy((s) for some E ((RR)). Then, under regularity conditions that we leave implicit.
(pp( = /// p(x) {e (e Vg(z) - ise ise ite- g(x) dr(x))
du(s).
(7)
The Fourier transform Vg of Vq is defined as fe- f Vx(x) dr, and basic property of the Fourier transform that the transform of derivative be computed using the expression is (ef(((((())))))) dr. This implies that the inner integral in (7) vanishes when and p are equal. Thus we can interpret GSD as quantification the uniformity of dz dd/d,, with a weighting function based on the Fourier derivative identity with regard to Vg.
15

---

Proof of Proposition 6. From direct calculation, and assuming derivatives and integrals can be in- terchanged, we have that
V Vak(a, v) = / ise isc--l) dj((s)
(8)
V Vok(z, v) / iss- du((s)
(9)
Vx Vyk(2, v) = / -115|2--(60000) du((s).
(10)
Now, let
næ 8) p(x) { (e, Vq(x) ise 1(, q(x) iss
((x) { - v log (xr)(()) - ise 1(8,z)
p(x)
and note through direct calculation and Equations (8) to (10) that
/ 112,0 8) 71. 8) du(ss)
g(x) ((v)
/ ||B|12 + log
} -o du(s)
p(x) p(v)
is V log (1v) + V log V Vlog g(v)
g(z) ((v)
4 V,k(r, V Fg y) Vlogg(2) Vak(z.v) log (v) - k)r,y) Vlogg(>)
p(x) p(y)
V Vlog((v)
((x) ((v)
kn (xr v))
p(x) p(y)
Thus, integrating with respect to n, and assuming that we may interchange the order of integrals, we have that
/1/ m(o,s) dr<(0)) du(s)= / /ʃ n(x, 8) nly, s) dr(x)dr(y)) du(ss
/] / n(xx, s) mly, 8) dp(s) dπ(z)d=(y)
((xr) g(y), kg v) dm dr(x)dr(y) (T) p(z) ply)
where the final equality Proposition 5. This establishes the result.
To prove Theorem 1, two intermediate results required:
Proposition 7. For an element π E (((Rd), assume Z == f(//p) dr E (0, oc). Assume that k and satisfy the preconditions of Proposition 2, and that f | log [|] da < oc. Let := == (p)) Then it E ((Rd) and
Dp (5).
Proof. The assumption Z E (0, co) implies that P(Rd). Furthermore, the assumption I log glla/(a-1((q/)) dı < co implies that f | log (|] da < cc. Thus the assump- tions of Proposition 2 are satisfied for both and म. and thus both Dp (T) and ( are well- defined. Now, with in Proposition 4, notice that
//
(xx)
D (T) = ||El|x(())=
Vk(z, + k(x, Vog
p(z)
H(()"
/ Vk(2, .) + k(x, logg(()) Zdz
ZD (T),
H(A)d
as claimed.
Proposition Let f RRd - |0, oc) and स E P(Rd). Then f f dr > 0 ʃs dr > for all (0, oo).
16

---

Proof. From the definition of the Lebesgue integral, we have that f dæ = sup s dr s simple function with 0 f s f > 0. Thus there exists simple function s = شi=1 s,ls, with 0 S & S f and f s dm > O. Here the 8; E IR and the measurable sets S; C R are disjoint. In particular, it must be the case that at least one of the coefficients s; positive; without loss of
generality suppose $1 > 0. Then & == $1 1/"1s, is a simple function with O I 8 I f fndd f dr > 0.
follows that J' f dæ = s dæ a simple function with S S f} > 0.
Proof of Theorem 1. Since (((/p) drn. (0, oo) and q/p 2 from Proposition 8 we have, for each n, that Zn := J q/p dm,, > O. Thus the assumptions of Proposition 7 are satisfied by & and each Wu, which guarantees that D P..g(TT)) Zn ZnD g. where म. ( E P(R.).
Now, since W (Ta, pt q/p) - O, taking f 1 we obtain Z, fq/p dı, + fq/p dp dp In addition, note that
| fd d - / f dq
W1 (T = sup
1(f)5)
/ f(0) + +f(xx) - f(0)] din(z)) - / f()) + (f()))(()) f(0) dg(x)
sup L(f)s1
| |[r(s) s(0) d (1(()) f(0)) dg(xr)
sup L(fis1
/ d - Idu
sup L(NS1 /(0==0
Thus, from the triangle inequality, we obtain the bound
f dā - f dg
W2(nnn 9)
sup L(f)s1 =(0)=0
fq dT - pZn
lq dp
sup L(FiS1 =(0)=0
p
Jq dma
14 d
/
1q dmn -
-
I1 sa dp
sup L(FiSI
sup Lifisi /(0)=0
pZn
p
p
p
=(0)=0
Zm
fq dmn
q
sup u =(0)=0
Wi
n p;
Zn
P
p
+0
-0
oc. For (+), since q/p 2 the supremum is realised by f(x) = ||>|| and
f ||r
((x) dr, (x) oc. p(x)
Thus we have established that W ( g) + 0. Since V log is Lipschitz with / V log g||² dg < Do and & has continuous and bounded second derivatives, the standard kernel Stein discrepancy has 1- I-Wasserstein convergence detection [Proposition 9 of Gorham and Mackey, 2017], meaning that
W, W1( (m. q) > 0 implies that D9.9 (āx) + 0 and thus. since Zn + Dp.g Pop.((TT) (ma) - 0. This completes the proof.
To prove Proposition 3, an intermediate result is required:
Proposition 9. Let k(x,v) v) = d(x y) be kernel with & twice differentiable and let q E P(Rd) with V log well-defined. Then k. (x,,)) = - AØ((0) + ((0][V log q(x)),? where A V. and kg was defined in Proposition 5.
17

---

Proof. First, note that we must have Vo(0) = 0. else the symmetry property of k would be violated. Now, Vak(z. (Vo)(= y) V V,k(Z,
(Vo) - y) and V V V,,(x, y) -ĀØ( - y)) Thus V Vak(a, ووا = V v,k(,, y) z=y = 0 and V Vyk(o, v) z=v (0)).
Plugging these expressions into Equation (5) yields the result.
Proof of Proposition 3. Let kn be defined as in Proposition 5. From Cauchy-Schward we have that (2 v) S / /kgl22 2) (ug v), and plugging this into Proposition 5 we obtain the bound
(T) - /
((2)) p(z) V ka(ar, dr(())
(11)
For a radial kernel k<,, y) = 6xx - y) with & twice differentiable, we have (0) > 0 (else k must be the zero kernel, since by 1Rd), and kg (2x,2) .r) - -)( + ( log ((x)] (from Proposition 9). Plugging this expression
[k( y)1 S x), y) o(0) for all æ, E
into Equation (11) and applying Jensen's inequality gives that
g(x))
(T)((T) <
p(z)² -2((0) + ((0) V logg(z))1 dz(z).
Now we may pick a choice of p. and (nm)neN (mn - p) for which this bound can be made arbi- trarily small. One example is NOO, 1), N(0, o²) (any fixed >> 1), for which we have
D2.( - / 0² ><<p(>>|a=]] [- (6( 600) [²] dr(x))
where 7 = o 0= > 0. Then it is clear that, for example, the sequence म = δ(ne1) (where ej = 11,0..... satisfies the assumptions of Proposition 2 and, for this choice,
(Tn( S o² exp((=>n²) [- ( + (0))n +
and yet Th + as claimed.
Proof of Theorem 2. Since infrers (<x)/p(x) > for each have Zn = f dm, and, furthermore, the assumption If Vl all /lolg//////)///////// dr, 00 implies that f V log q||a/(a-1) dmn < oc. Thus the assumptions of Proposition 7 are satisfied by and each Tn, and thus we have (T.)(Tn) = Z,, (z. where म, == (gmm)/ppZn) E P(R.).
From assumption, Z, 2 inf ((x)/p(x) is bounded away from O. Thus if D (TT) + 0 then Da.q(Tr, (T,)) + O. Furthermore, since q E Q(E:)) and the inverse multi-quadric kernel & is used, the standard kernel Stein discrepancy has convergence control, meaning that (q.g(#Tn) - implies
nn d q [Theorem of Gorham and Mackey, 2017]. It therefore suffices show that än 4 implies n. - p.
From the Portmanteau theorem, Tn 4 p equivalent to f dmn - g dp for all functions g which are continuous and bounded. Thus, for arbitrary continuous and bounded function g. consider
f= -gp/g. which is also continuous and bounded. Then, since म. 4 q. we have (again from the Portmanteau theorem) that Z,- f dr., = f dæ, + f dg = f g dp. Furthermore, the specific choice g shows that Z, 1, and thus f g dr,, f dp in general. Since y was arbitrary,
have established that म. 4 p. completing the proof.
To prove Theorem 3, an intermediate result is required:
Proposition 10. Let Q E P(Rd)) and let kg IRd - R be reproducing kernel with (x, dg 0 for all Rd. Let (In)neN be a sequence of random variables independently sampled from and assume that (x, x)} dg(x)) < 2o for some y >> 0. Then
D 9-2
Σ ((xz)
- 0
n
1=1
almost surely as n 20.
18

---

Proof. This is Lemma 4 in Riabiz et al. [2022], specialised to the case where samples are indepen- dent and identically distributed. Although identical the statement in Riabiz al. [2022], one obtains this result by following an identical argument and noting that the expectation of (21, Tz)) identically O when * j (due independence of z; and σz), so that bounds on these terms are not required.
L
Proof of Theorem 3. Since has finite support, all conditions of Theorem 2 are satisfied. Thus it is sufficient to show that almost surely pp.g( (Tx) ) - To this end, we follow Theorem of Riabiz et [2022] and introduce the classical importance weights
p(x;)// (/())) which are well-defined since Q > 0. The normalised weights ú; = W// W, W, == ブÿ= "y satisfy S W1,..., th and W1 + +.... W = 1, and thus the optimality of w", together with the integral form of the GF-KSD in Equation (4), gives that
Dp.g Σ w b(xx))
s <Dpa Drn
Σ (.:((z))
kig (g)))))))) zzj)
Wn,
im
j=l
Wn,
D
Σ 6(1))
(12)
From the strong law of large numbers, almost surely n=W,, + f dq Thus it suffices to show that the final term Equation (12) converges almost surely to To achieve this, we can check the conditions of Proposition 10 are satisfied.
Since q and f-)) Spak(z,) satisfy the conditions of Proposition the condition I (2, dq O is satisfied. Let ((=) = + ||= so that k(x, = 6(x - y) is the inverse multi-quadric kernel. Note that (0)) = 1 and Ad(0) = -d. Then, from Proposition 9, we have that (g x) (0( + ((0) log log all = d + ||V log zl12. Then
/ exp(>kk x}} dg(x) = exp(nd] / /epp/1llV log d go
dg < co.
which establishes that the conditions of Proposition 10 satisfied and completes the proof.
B Stein Variational Inference Without Second-Order Gradient
This section contains sufficient conditions for unbiased stochastic gradient estimators to exist in the context of Stein Variational Inference; see Section 4.2 the main text. The main result that we prove is as follows:
Proposition 11 (Stochastic Gradients). Let p, q, R P(Rd) and T Rd - Rd for each 0 E IRP Let f + V T be bounded. Assume that for each 9 E IRP there is an open neighbourhood No C
IP such that
sup /
(1))
dR(x) < co,
BeNo (TT)
I sup
(TT
V log (x)) dR(x) < co.
oE P(T (x))
I Be P(T (x))
(T (x))
2 log (7 (x))] dR(x) < oo,
sup
for each of r E {..q). Let le be the inverse multi-quadric kernel in Equation (2) and let u(ar y) denote the integrand in Equation (6). Then
V VoD p (m)) E
Σ Vou(T (x:), T (2;))
n(n
where the expectation is taken with respect independent samples æ1, Z1...., In R.
19

---

The role of Proposition 11 is to demonstrate how an unbiased gradient estimator may be constructed, whose computation requires first-order derivatives p only, and whose cost is O(n²). Although the assumption that + Va2 (x) is bounded seems strong, it can typically be satisfied by re- parametrisation of E IRP.
To prove Proposition 11, we exploit the following general result due to Fisher et al. [2021]: Proposition 12. Let R E P(RP). Let è S IP be an open set and let u : IRd x R + R. T R Rdd 8 O, be functions such that, for all 0 E E
(A1) If [u(TT (x), (y)) dR(x)dR(y) < oo;
(A2) there exists an open neighbourhood No C & of 0 such that
sup || OEN.
(Vou (1),, T (v)) dR(x)dR(y) < oo.
Then F(0) = JJ ul (x), dR(z)dR(y) is well-defined for all Q and
VaF(O) = E
Σ (x)
T
n
if3
where the expectation taken with respect to independent samples #1,... In ~ R.
Proof. This is Proposition 1 in Fisher et al. [2021].
Proof of Proposition 11. what follows we aim to verify the conditions of Proposition 12 hold for the choice u(x, = kg (ar, y). where ki was defined in Equation (5).
(Al): From the first line in the proof of Proposition 4, the functions Sp.g((,)) +), Rd, are in M(k)². Since (x, y) + u(x, y) (z, -), (g:k)((()))) is positive semi-definite, from Cauchy-Schwarz, (u(x,)) I Vu(z, m) u(v.v). v). Thus
1 (c), (v)) dR(x)dR(y) - - [u([,,)] aT R(x)dT R(y)
$( ( a R(x)))
Since k(2,)) $(x - v). we have from Proposition 9 that
x(x = (x(x)
[- [-))(0] + Voo)
p(x))
and
J Vu(z,a) R() s
v /
2
/11/11/-0000 -0(0) <0]] log dr
dTg
y By/
p
which is finite by assumption.
(A2): Fix x. y E Rd and let R2 (0) = q(T (x))/ (v)). From repeated application rule of differentiation, we have that
appicttall of the product
Vou (x), T (v)) (T (x), T (v)) [R (R) + R. (O)
Let bp(x))= := Vogg(w)) (x)) V log ((r), b(x) = -,()) and (V.7 (x))k, (1) (z). In what follows, we employ a matrix norm dxd which is consistent with the Eu- clidean norm on Rd meaning that (x) (x))] <
for each 0 e e and æ E Rdd Considering the first term (*), further applications of the chain rule yield that
(z)
V. (0)R, = R. (0) (z) + V (T)b(T"()) (v))
20

---

and from the triangle inequality we obtain a bound
V (O) I Rs (0)
(
IV.7 (v)]l (v))1l (v))lI]
Let S denote inequality up an implicit multiplicative constant. Since we assumed that (oTT())) is bounded, and the inverse multi-quadrie kernel k is bounded, we obtain that
[(*) < R2 (,) [[]6(T (z)))I
(2))1l.
Similarly, from Equation (5), and using also the fact that the inverse multi-quadric kernel & has derivatives or all orders [Lemma 4 of Fisher et al., 2021], we obtain a bound
(T (y)il & [1 + | (T (z)) + IV (T (x))i] [1+ (T6 (x)II
[ + b (T (y))1 Vby (T (w))l] [1 + (T (z))
which we multiply by R, (0) to obtaín bound on (++). Thus we have an overall bound
[ (x). (v)) || & R. (,) (0) { + (T, (z)) + (T)) [] 1 (TT
+[1+ 16, ( Vb, [y) [1 (TT
Substituting this bound into Jf supseN, oT (z), T || dl(()d/B(y)) and factoring terms into products of single integrals, we obtain an explicit bound this double integral in terms of the following quantities (where E {p.
sup (0) dR(x)
BEN,
| geN,
sup R2 (1)0,,, dR(z)
[ sup BEN,
(T)) (R(x)
which we have assumed exist.
Thus the conditions of Proposition hold, and the result immediately follows.
C
Experimental Details
These appendices contain the additional empirical results referred to in Section 3, together with full details required to reproduce the experiments described in Section 4 of the main text.
C.1 Detection of Convergence and Non-Convergence
This appendix contains full details for the convergence plots of Figure 1, In Figure 1, we considered the target distribution
p(z) w.N(z;p
where N(x; /4, d²) is the univariate Gaussian density with mean and variance a². The parame- ter choices used (wi, W2, Wz) (0.375, 0.5625, 0.0625), (21, us) = (-0.4, 0.3, 0.06) and (of, oz. oz) =
The approximating sequences considered were location-scale sequences of the form L where L"(x) a, b, for some (an)neN and (bı)nεN and P(R). For the converging sequences, we set u = p and for the non-converging sequences, we set u = (0, We considered three different choices of (a, heN and (b,)neN, for each colour. The sequences (a,)neR and (b,)) used are shown in Figure S1. The specification of our choices of q the following:
Prior: We took g NO 0.75².
Laplace: The Laplace approximation computed A/033, 0.2041
21

---

3.0
dg
dn
b
bn
20
2.0
1.5
0.5
0.5
F1
71
(a)
(b)
Figure S1: The sequences On and b, used the location-scale sequences. (a) The sequences an and bn used in the location-scale sequences in Figure 1. (b) The sequences and b, used in the location-scale sequences in Figure 2. In each case, the colour used of each curve indicates which of the sequences (I neN they correspond to.
GMM: The Gaussian mixture model was computed using 100 samples from the target p. The number of components used was 2, since this value minimised the Bayes information criterion [Schwarz, 1978].
The kernel density estimate was computed using 100 samples from the target p. We utilised a Gaussian kernel k(r,,) v) exp(-(e (x: /(²) with the lengthscale or bandwidth parameter & determined by Silverman's 's rule of thumb [Silverman, 1986].
KDE:
The values of GF-KSD reported in Figure 1 were computed using a quasi Monte Carlo approxi- mation to the integral (6), utilising a length 300 low-discrepancy sequence. The low discrepancy sequences were obtained by first specifying a uniform grid over [0, 1] and then performing the in- verse CDF transform for each member of the sequence Tp.
C.2 Avoidance of Failure Modes
This appendix contains full details of the experiment reported in Section 3.3 and an explanation of the failure mode reported in Figure 2a. The sequences considered displayed in Figure S2. Each sequence was location-scale sequences of the form L u, where L"(x) On + b, for some hee) and (b,), and u E (IR). For the converging sequences, we set u = p. The specification of the settings of each failure mode are follows:
Failure mode (a) [Figure 2a]: We took p as the target used in Figure and detailed in Appendix C.I and took q (0, 1.5"). The a, and b, sequences used are displayed in Figure Sla. The values GF-KSD reported were computed using quasi Monte Carlo approximation, using length 300 low discrepancy sequence. The low discrepancy sequences were obtained by first specifying a uniform grid over (0, 1] and then performing the inverse CDF transform for each member of the sequence n..
Failure mode (b) [Figure 2b]: We took p ~~ NO, 1) and g ~ N(--0.,0.13). The a, and by sequences used are displayed in Figure Slb. The values of GF-KSD reported were computed using a quasi Monte Carlo approximation, using length 300 low discrepancy sequence. The low discrepancy sequences were obtained by first specifying a uniform grid over [0, 1] and then performing the inverse CDF transform for each member of the sequence π.
Failure Mode (c) [Figure 2c]: each dimension considered, we took p ~ N/O, I) and q ~ (0, 1.17). The an and b, sequences used are displayed in Figure Slb. The values of GF-KSD reported were computed using a quasi Monte Carlo approximation, using a length 1, 024 Sobol sequence in each dimension d.
NNot that for d 1, still considered location-scall sequences of the form (x) an b.
22

---

12 12
13
0.
0.15
044
0.00
0.85
0.00
2.00
1.75
Consrrgigg Sequene
1.90
Sequiiee
Figure S2: Test sequences (T,, ne used in Figure 2. The colour and style of each sequence indicates which of the curves in Figure is being considered. In the second row from the top, the sequence used when d Figure 2c shown in the final column.
Failure mode (d) [Figure 2d]: We took p= with p(z) 0,5.( 0.5. 0.5.V(z; 1, 0,13), where N(x; 14, o²) is the univariate Gaussian density with mean and vari- ance d². The a, and b, sequences used are displayed in Figure S1b. For the non-converging sequences we took u = (11, 0.1 and used the a, and b,, sequences specified in Figure S1b. The values GF-KSD reported were computed using quasi Monte Carlo approximation, using a length 300 low discrepancy sequence. The low discrepancy sequences were obtained by first specifying uniform grid over [0, 1] and then performing the inverse CDF transform for each member of the sequence T..
In Figure S3, provide account of the degradation of convergence detection between g Prior considered in Figure and g = (0, 1.52) of Failure mode (a). In Figure S3a, it can be seen that the value of the integrals J (q/p)² dm,, are finite for each element of the pink sequence πn- However, in Figure S3b. it be seen that the values of the integrals f (g/pp) dıy, are infinite for the last members of the sequence ,., thus violating a condition of Theorem 1.
C.3 Gradient-Free Stein Importance Sampling
This appendix contains full details for the experiment reported in Section 4.1. We considered the following Voltera dynamical system:
i(t) = a'u a'u(t) ('u(t)o(t),
a(t)
5'' ('u(t))(tt)
+
(u(0), v(0) = (u6, v%).
23

---

Prior
4= VV/00 1.52)
/3/p
M
4=M0.
4=
1.39
==M0.14)
0
x
(b)
(c)
Figure S3: Explanation of Failure Mode (a). (a) Values of (q/p) for g Prior and for the
converging pink sequence displayed in Figure (b) Values of (g " for g = / 1.5%) and for the converging pink sequence displayed in the first column and first of Figure S2. (c) Values of q², /p for different choices of q.
Using 21 observations U1,..., U21 and U1 U...... U21 over times t <... < 22 we considered the probability model
U; ~ Log- Looggnorml((lou((t)).
0 Log-normal(log v(f.), (o1)²).
In order to satisfy positivity constraints, we performed inference on the logarithm of the parameters (a, 3, 7, 6, up Uo. 11. σ2) (log a', log 3', log log 8', log u/. log 1/1 log o' log o2). We took the following independent priors on the constrained parameters:
a' Log-normallogg(0))...00 0.6²), g' Log-normal(lo((..000 0.3²),
?/ (oggnnoral((log(0.)))00 0.6?), 8' Log-normal((og(0002)) 0.3?).
" Log-normalllog (10). 1), v Log-normal(oog(0); 1),
of Log-normal(log(12222 0.02?), oe2 (og-normal(og(( 0.02?).
In order to obtain independent samples from the posterior for comparison, we utilised Stan [Stan Development Team, 2022] to obtain 8, 000 posterior samples using four Markov chain Monte Carlo chains. Each chain was initialised at the prior mode. The data analysed are due to Hewitt [1921] and can be seen, along with posterior predictive check, in Figure $4.
The Laplace approximation was obtained by the use of 48 iterations of the L-BFGS optimisation algorithm [Liu and Nocedal, 1989] initialised at the prior mode. The Hessian approximation was obtained using Stan's default numeric differentiation the gradient.
Finally, the quadratic programme defining the optimal weights of gradient-free Stein importance sampling (refer to Theorem [2016].
solved using the splitting conic solver of O' O'Donoghue et al.
C.4 Stein Variational Inference Without Second-Order Gradient
This appendix contains full details for the experiment reported in Section 4.2. We considered the following bivariate densities
(i := : (x; 0, ni) (1y (ax), m22)
(2 u) (x;; 0, of) ( lxz ag)). ,x²,,?))
where N(x; μl, o²) is the univariate Gaussian density with mean μ and variance σ². The parameter choices for the sinusoidal experiment P1 7/Î 1.32 n7 0.092 and choices for the banana experiment P2 were of = 1, 0.22 and = 0.5.
1.2. The parameter
The development of robust stochastic optimisation routine for measure transport with GF-KSD beyond the scope this work, and what follows we simply report one strategy that was suc- cessfully used in the setting of the application reported in the main text. This strategy was based on

---

Lymx
80-
60-
0-
Hade
80-
80-
40-
20-
0-
1500
1908
1910 Year
1915
1920
Data
Prodiction
Figure S4: Posterior predictive check for the Lotka- Volterra model. The shaded blue region indi- cates the 50% interquartile range of the posterior samples.
0.4
0.2
2500
3000
7500
10000
12500
15000
17500
20000
Figure S5: Tempering sequence (Em)meN used each the variational inference experiments.
tempering of p. the distributional target, to reduce a possibly rather challenging variational optimisa- tion problem into sequence of easier problems to be solved. Specifically, we considered tempered distributions Pm ((IRR)) with log density
log Pm (z) cm Pog (1 cm. logp(z))
where (cm)mes MEN E [0, 1 the tempering sequence and Po E P(Rd) is fixed. In this case Po was taken to be (0, 27) in both the banana and sinusoidal experiment. Then, at iteration of stochastic optimisation, we considered the variational objective function
π + log D (TT
where y = mo... as explained the main text. Tempering has been applied in the context of normal- ising flows in Prangle and Viscardi [2023]. The tempering sequence used (Em)mEN for each of the experiments is displayed in Figure S5.
For each experiment, the stochastic optimisation routine used was Adam [Kingma and Ba, 2015] with learning rate 0.001. Due to issues involving exploding gradients due the q/p term in GF- KSD, we utilised gradient clipping in each of the variational inference experiments, with the max- imum 2- 2-norm value taken be 30. both the banana and sinusoidal experiment, the parametric class of transport maps T was the inverse autoregressive flow of Kingma et al. [2016]. In the banana experiment, the dimensionality of the hidden units in the underlying autoregressive neural network was taken 20. In the sinusoidal experiment, the dimensionality of the hidden units in the underlying autoregressive neural network was taken as 30. For the comparison with standard kernel Stein discrepancy, the same parametric class 8 and the same initialisations of 9 were used.
25

---

1.0
0.6
2
0.2
-1. -1.0
(a)
(b)
Figure S6: The To and tempering sequences used in additional convergence detection experi- ments. The colour of each curve indicates which of the sequences in Figure S7 and Figure S8 they correspond with. (a) The mi choices each tempered sequence of distributions. (b) The tempering
sequences mEN considered.
C.5 Additional Experiments
Appendix C.5.1 explores the impact of the conclusions drawn in the main text. Appendix C.5.2 investigates the sensitivity of the proposed discrepancy to the choice of the parameters and 3 that appear in the kernel. Appendix C.5.3 compares the performance of GF-KSD importance sampling. KSD importance sampling and self-normalised importance sampling.
C.5.1 Exploring the Effect of
In this section investigate the robustness of the convergence detection described in Figure subject to different choices of the target We consider two further choices of
(1
C (x;
P2(x) Σ d, Student- u,
where N(r; d is the univariate Gaussian density with mean μ and variance d² and Student-T(x; M, m) is the univariate Student-T density with degrees of freedom location pa- rameter m and scale parameter s. The parameter choices for Pl were
(c1, C2, c3, 01) (0.3125, 0.3125, 0.3125, 0.0625),
(Mi, U2, u3, p4) = (-0.3, 0.3, 0),
(of, o2 o oi) =
The parameter choices for P2 were v = 10 and
(d,, d2, d3, d4) = 0.21, 0.3,0 0.4), (m,, m2, M3, m4) = 0..2 (81, 82, 83, 84) = (0.05, 0.1,0.1..111
Instead of using the location-scale sequences of Figure 1, we instead considered tempered sequences of the form
log n (z) En log Tu (2 + (,, log uxz).
For the converging sequences considered we to be the target (either Pı P2) and set u = A (x;; O, 0.42) for each of the non-converging sequences. The different sequences vary in choice of mo and tempering sequence (en)neN. These choices displayed in Figure S6 and are taken the same for both of the targets considered.
The specification of our choices of is the following:
26

---

Prior: For P1, we took q=0,,= For we took q ~ Mondent-T(11,00)) 0.5))
Laplace: For the Laplace approximation computed was (0, 0.051²). For p2, the Laplace approximation computed was q (O, 0.125").
GMM: For both targets, the Gaussian mixture model was computed using 100 samples from the target. both cases, the number of components used was 3, since this value minimised the Bayes information criterion [Schwarz, 1978].
KDE; For both targets, the kernel density estimate was computed using 100 samples from the target. both cases, we utilised a Gaussian kernel k(x,, y) = exp(-(x y) / with the lengthseale or bandwidth parameter é determined by Silverman's 's rule of thumb [Silverman, 1986).
Results for P1 are displayed Figure S7 and results for Ps are displayed in Figure S8. can be seen that for both target distributions and the different sequences considered, GF-KSD conrectly detects convergence in each case. For both targets and for q = Laplace, can be seen that GF-KSD exhibits the same behaviour of Failure Mode (b), displayed in Figure 2b.
The values of GF-KSD reported in Figure and Figure S8 were computed using quasi Monte Carlo approximation to the integral (6), utilising a length 300 low-discrepancy sequence. Due to the lack of an easily computable inverse CDF, we performed an importance sampling estimate of GF-KSD as follows
Dp-(TT) = // (Spa & Spa)k((,, v) dr(x)) dr(v)
= // (Sp. Spg) w(r)w(y)
(()))(v) dw(x) dwlyy,
where is the proposal distribution. For each element of a sequence Tne we used Gaussian proposal wn, of the form:
log Wh, (I:) = G, log (x) + (1 - c.) log [o O, 0.4").
Since To is Gaussian for each sequence, this construction ensures that each W,, both Gaussian and good proposal distribution for π. The low-discrepancy sequences were then obtained by first specifying uniform grid over [0, 1] and the performing an inverse CDF transformation using wy-
C.5.2 Exploring the Effeet of a and 3
In this section we investigate the effect convergence detection that results from changing the parameters and | in the inverse multi-quadric kernel (2). Utilising the same test sequences and choices of used in Figure I, we plot the values of GF-KSD in Figure S9. It can be seen that the convergence detection is robust to changing values of a and B.
C.5.3 GF-KSD KSD Importance Sampling
In this section we investigate the performance of gradient-free Stein importance sampling, standard Stein importance sampling, and nelf-normllised importance sampling, as the distribution q varies in quality approximation We consider two different regimes:
1. - -(,, and = -N AI) for 0.7 s A S 1.3.
2. -N. I) and ()(11 for -0.6 S S 0.6, where = (I)
In both cases, we consider the performance of each approach for varying dimension d and number of samples Results are reported in Figure S10 and Figure S11 for each regime respectively. The quadratic programme defining the optimal weights of gradient-free Stein importance sampling and Stein importance sampling (refer to Theorem 3) was solved using the splitting conic solver of O' Donoghuue et al. [2016].
27

---

Convergence
q= Prior
q= Lapaee
20
log GF-KSD
Octward Convergmne
Ouneanl CoreCeeeeet
-10
GMM
q= KDE
log GF-KSD
Converpence
Oblique Convergeace
25
75 25
f1 100
(a)
(b)
Figure S7: Additional empirical assessment of gradient-free kernel Stein discrepancy using the target P1 defined in Appendix C.5.1. (a) Test sequences (mn))EN, defined in Appendix C.5.1. The first column displays sequences (solid) that converge the distributional target p (black), while the second column displays sequences (dashed) which converge instead to fixed Gaussian target. (b) Performance of gradient-free kernel Stein discrepancy, when different approaches to selecting g are employed. The colour and style of each curve in (b) indicates which of the sequences in (a) is being considered. [Here we fixed the kernel parameters o and 3 = 1/2.]
Inward Convergence
Iowaed Nou. Convergence
q= Prior
9 Laplace
20
log GF-KSD
Octoorli Consergnne
Ouricenl Nan-Ceneergeter
-20
GMM
q= KDE
20
log GF-KSD
Oblique Converpence
Obligue Converpence
-10
-20
11
(a)
(b)
Figure S8: Additional empirical assessment of gradient-free kernel Stein discrepancy using the target p2 defined in Appendix C.5.I. (a) Test sequences (TT))E,, defined in Appendix C.5.1. The first column displays sequences (solid) that converge to the distributional target p (black), while the second column displays sequences (dashed) which converge instead to fixed Gaussian target. (b) Performance of gradient-free kernel Stein discrepancy, when different approaches to selecting are employed. The colour and style of each curve in (b) indicates which of the sequences in (a) is being considered. [Here we fixed the kernel parameters o = and 3 = 1/2.]
28

---

q= Prior
9= Laplace
9 GMM
q= KDE
(a,f)= (0.5.0.5)
(o.0)=
(o.p)= (o.p)= (o.0)=
Figure S9: Comparison of different values of a and in the inverse multiquadric kernel. Here the vertical axis displays the logarithm of the gradient free kernel Stein discrepancy. The colour and style of each of the curves indicates which of the sequences in Figure is being considered.
29

---

m 10
100
500
-1.75
GF-KSD
log ED (d= 1)
Self Normalised
55
-2.75
-60
-3.00
-6.
-3.25
-2.75
log ED (d=
-08
-10
-12
-5.00
-14
-1.
-6.00
log ED (d= 20)
0.90
0.00
-0.25
-3
1.0
1.2
1.2
0.8
1.5
log HD (d==0)
1.0
1.0
1.2
A
Figure S10: Comparison of the performance of importance sampling methodologies in varying di- mension d and number of sample points considered under the regime AT),, The approx- imation quality is quantified as the logarithm of the Energy Distance (ED).
30

---

100
500
-2.00
4.0
log ED (d=1) 11=1)
-2.25
KSD
-2.50
-50
Self Normalised
-3.00
-3.25
-3.50
-65
-10
0.0
log ED (d=5)
4.0
-2.5
4.5
-30
-5.0
-1.0
-5.5
-60
15
log ED (d= 20)
0.0
-0.5
0.0
0.5
-0.5
0.0
-05
0.0
0.5
2.0
log ED (d= 50)
0.5
-03
Figure S11: Comparison of the performance of importance sampling methodologies in varying di- mension d and number of sample points considered n under the regime NN I). The approxi- mation quality is quantified as the logarithm the Energy Distance (ED).
31


Reviews:
[{"review_id": "krK9MxBxu9", "reviewer": "Reviewer_SmKF", "summary": "The authors study the kernel Stein discrepancy based on a Stein operator, introduced previously in [Liu 2018], which does not require access to the gradient of the target density. As their main theoretical result, they prove that, under certain assumptions on the target distribution (and the auxiliary distribution q) and the approximating sequence, the discrepancy controls weak convergence. The authors provide recommendations for the choice of the auxiliary distribution q in the construction of the discrepancy and present an experimental study which identifies certain failure modes of the discrepancy. They also discuss two applications to posterior approximation.", "questions": "Related to what I wrote above, could the authors state clearly what class of targets $p$ their Theorem 2 applies to? Could they also try to characterise the class of targets $p$ for which one can easily construct a useful auxiliary distribution $q$ in practice (using one of the methods of section 3.1), such that $q$ and $p$ satisfy the assumptions of Theorem 2?", "limitations": "I believe the real applicabilty of the theoretical results presented by the authors is somewhat limited, to an extent that is not clearly acknowledged in the paper. But I am looking forward to reading the authors' response as perhaps I'm missing something.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The idea of the authors to use the gradient-free Stein operator in order to define a new kind of KSD is interesting and worth studying. There are certainly a lot of situations where the evaluation of the gradient of the model is costly in practice. Giving the practitioners an option to avoid evaluating it when computing the KSD is certainly useful. The discussion of the properties of the new GF-KSD is very thorough and based on both theoretical results and a substantial experimental study. I particularly liked the the experiments revealing failure modes of GF-KSD as those are not necessarily captured by the theoretical results.", "weaknesses": "The main theoretical contribution of the paper seems to be provided by Theorem 2. However, I am a bit concerned about its real applicability. The authors require that the auxiliary distribution $q$ have tails heavier than the target $p$. On the other hand, they require that $q$ is distantly dissipative, which precludes the use of heavy-tailed $q$. This all means that $p$ must not be heavy-tailed if the assumptions of Theorem 2 are to be satisfied. \n\nMoreover, the authors suggest four potential choices of $q$ (Prior, Laplace, GMM and KDE) but note that GMM and KDE are impractical in general. At the same time, the Prior method does not satisfy the assumptions of Theorem 2 for heavy-tailed priors. The Laplace method, on the other hand, seems to satisfy the assumptions of Theorem 2 only for targets which are sub-Gaussian.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful feedback, especially at short notice, which is greatly appreciated.\n\n> This all means that $p$ must not be heavy-tailed if the assumptions of Theorem 2 are to be satisfied.\n\nThis is completely correct - but we emphasise that the standard kernel Stein discrepancy (KSD) also has this requirement.  That is, the standard KSD is only guaranteed to have convergence control in settings where the target distribution $p$ is distantly dissipative (Theorem 8 of Gorham and Mackey, 2017), meaning that $p$ cannot be heavy tailed.  Despite this limitation, standard KSD has been widely used for applications such as sampling and variational inference, for which convergence control is required.  This is encouraging, in the sense that there are a wide variety of important problems for which $p$ is not heavy tailed and KSD has been successfully used.  There is active research on Stein's method for heavy tailed targets in the Probability community [for example, UB2022], but (as far as we are aware) the application of these techniques to the problem of posterior approximation has yet to be attempted.\n\nThank you for raising this point - we will explicitly highlight the non-applicability of both KSD and GF-KSD to heavy-tailed $P$ in the revised manuscript, when discussing the preconditions of Theorem 2.\n\n> Related to what I wrote above, could the authors state clearly what class of targets $p$ their Theorem 2 applies to?  Could they also try to characterise the class of targets $p$ for which one can easily construct a useful auxiliary distribution $q$ in practice (using one of the methods of section 3.1), such that $q$ and $p$ satisfy the assumptions of Theorem 2?\n\nThank you for the opportunity to clarify this point -- we first provide a mathematical answer, and then also a practical answer which addresses the fact that an appropriate and actionable choice of $q$ is required.\n\nMathematical answer:\n\nIf we choose $q=p$ we obtain the standard KSD, and in this case the class of $p$  for which KSD provided convergence control is known to include all distributions $p(x)$  on $\\mathbb{R}^d$  with Lipschitz continuous $\\nabla \\log p(x)$  and _distantly dissipative_ tails (see Section 2.2). The latter condition is implied when $p$  is strongly log-concave outside a compact set (i.e. sub-Gaussian tails).  \n\nFor any $p$ (not necessarily distantly dissipative) with sub-Gaussian tails there exists a dominating Gaussian distribution $q$, i.e. with $\\inf_{x \\in \\mathbb{R}^d} q(x) / p(x) > 0$, and since Gaussians are distantly dissipative Theorem 2 can be applied.  Thus GF-KSD is in principle _more widely applicable_ than standard KSD.\n\nWe will make sure to emphasise this point in the revised manuscript.\n\nPractical answer:\n\nOf course, if $p$ is implicitly defined then it may not be clear how $q$ should be selected.  If $p$ itself is distantly dissipative then $q$ could be taken to be a tempered version of $p$ (i.e. $q(x) \\propto p(x)^t$ for some $0 < t < 1$), which would automatically guarantee that $\\inf_{x \\in \\mathbb{R}^d} q(x) / p(x) > 0$.  But this \"solution\" raises the question of how to take the gradient of $q$, if taking the gradient of $p$ is presumed to be difficult.  As a more practical strategy, we propose to take $q$ to be a Laplace approximation of $p$, in the hope that the curvature of $p$ at the mode is an accurate proxy for the tail decay of $p$.  This strategy was seen to perform well in the experiments reported in Sections 3.1 and 4.1, but one can of course construct examples where this strategy will fail.  \n\nIn preparatory work, we also considered taking $q$ to be a Laplace approximation with an \"inflated\" covariance matrix; a conservative choice aimed at guaranteeing the $\\inf_{x \\in \\mathbb{R}^d} q(x) / p(x) > 0$ condition is satisfied.  This did not lead to improved performance in any of our experiments, so was omitted from the manuscript.\n\nPlease allow us to emphasise that this is an instance of the fundamental problem that if $p$ is implicitly defined via Bayes' theorem, such that the mathematical properties (tail decay, curvature, etc) of $p$ are unknown to the user, then we cannot hope to predict how well _any_ sampling algorithm will work.  In specific applications where the prior and the likelihood are mathematically tractable, we may be able to _derive_ appropriate choices of $q$ using similar methodology that is used for constructing dominating measures for rejection sampling; this would be an interesting direction for further work.  \n\nAll of the above discussion will be incorporated into the revised manuscript.\n\nThank you once again for your comments and questions, which we hope we have adequately addressed.  To summarise; the theoretical applicability of GF-KSD is strictly greater than standard KSD in terms of convergence control, but neither KSD nor GF-KSD are appropriate when $p$ is heavy-tailed.  The practical choice of $q$ requires some care, but a promising strategy based on Laplace approximation was found to work well in the experiments we performed.  We hope that our commitment to clarify these points in the manuscript will allow you to increase your score for our work.\n\nReferences:\n\n[UB2022] Upadhye NS and Barman K, 2022. A unified approach to Stein’s method for stable distributions. Probability Surveys, 19, pp.533-589."}, {"review_id": "nbccl0MuWj", "reviewer": "Reviewer_9g9f", "summary": "This paper proposed a posterior approximation using a new Stein discrepancy, which does not require derivatives of the statistical model. For that purpose, the authors derived the new discrepancy, called gradient-free KSD, and studied its statistical and convergence behaviors theoretically.\nThen the authors developed algorithms for differential equations, for which stable computation of derivatives is difficult, and a new sampling algorithm that bypass the Hessian calculations.", "questions": "I would appreciate it if the authors would answer the concerns described in Weakness.\n\nAs for minor questions;\n- What do the dotted and solid lines in Figure 1(b) correspond to?\n\n- Is it required to adjust parameters of the Laplace distribution, KDE, and GMM  for $q$ in some way ? If so, what is the recommended method ?\n\n- Looking at Figure 3 (b), it seems that the number of samples ($n$) must be very large ($\\log n=5$, i.e., $n=150$) even for low-dimensional problems such as d=8 in order for there to be any difference in energy distance. Is my understanding correct?", "limitations": "The limitation of the proposed method is discussed in detail.", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- A new KSD that does not require gradients, similar to the idea of importance sampling, is proposed. This leads to new algorithms for differential equations where the gradient is difficult to compute and for problems requiring Hessian calculation.\n\n- The authors studied the theoretical property of the proposed GF-KSD by extending the existing KSD theory.\n\n- Not only theoretical analysis but also detailed numerical investigations of the choice of parameters and $q$ are carried out with the actual use in mind.\n\n", "weaknesses": "- The writing style is such that the main paper alone is not complete, and it is assumed that the reader will read the Appendix. For example, Eq. 6 of Line 115 does not appear in the main text, and the tilted Wasserstein distance defined in Theorem 1 is introduced without any explanation of its properties in the main text.\n\n- The writing style could be improved since the discussion about existing research and the explanation of the proposed method are mixed, making the paper difficult to read.\n\n-Some parts are mathematically undefined or under-discussed\n - In Definition 2, sup is undefined\n - In Line 158, at last, $\\not \\to$  is undefined.\n - I don't know how widely the tilted Wasserstein distance (TWD) in Theorem 1 is known to the general public, but there is no discussion of the properties of TWD. Therefore I could not understand how important Theorem 1 is, that is, how important it is when it is said that convergence of TWD leads to convergence of GF-KSD; even after reading the proof of Theorem 1, I could only understand that TWD is a convenient form of the usual Wasserstein distance, which is obtained after applying the triangle inequality.\n\n- I could not understand the importance of the proposed method because I am not sure for what problems the proposed method is effective. \nI agree that it may be useful for differential equation problems, but the authors only applied the method to very small models of Lotka-Volterra in the experiments. In such a setting, MCMC is the standard approach, and even if the model is high-dimensional, we can solve it efficiently by variational inference. \nAlso, although the combination with Stein variational inference seems interesting, I wondered if GF-KSD is really flexible enough to generate samples for complex real data under the two restrictions suggested in Section 3. \nOne restriction is that the tail of $q$ should not be far from the target distribution; the other is that it must not be high-dimensional. \nI think the application to differential equations seems promising, so it would be better to find a problem setting where GF-KSD is more useful than MCMC and standard variational inference.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking a look at our manuscript.\n\nPlease see our global response to all Reviewers on the topic of how our manuscript has been presented.\n\n> In Definition 2, sup is undefined\n\nThe symbol \"$\\sup$\" is the supremum; our understanding is that we can assume familiarity with \"$\\sup$\" for NeurIPS, but we respectfully defer to the Area Chair for guidence on this point.\n\n> In Line 158, at last, $\\nrightarrow$ is undefined.\n\nThe symbol \"$\\nrightarrow$\" is the logical opposite of $\\rightarrow$, i.e. \"does not converge\"; our understanding is again that we can assume familiarity with \"$\\nrightarrow$\" for NeurIPS, but we again respectfully defer to the Area Chair for guidence on this point.\n\n> I don't know how widely the tilted Wasserstein distance (TWD) in Theorem 1 is known to the general public, but there is no discussion of the properties of TWD. Therefore I could not understand how important Theorem 1 is, that is, how important it is when it is said that convergence of TWD leads to convergence of GF-KSD; even after reading the proof of Theorem 1, I could only understand that TWD is a convenient form of the usual Wasserstein distance, which is obtained after applying the triangle inequality.\n\nThank you for the opportunity to clarify this point:  The TWD was introduced quite recently, in Proposition 3.3 of Huggins and Mackey 2018, so it is not widely known at present.  The special case where the tilting function $g(x) = 1$ is constant corresponds to the standard 1-Wasserstein distance.  Loosely speaking, for a tilting function $g(x)$ that is bounded away from $0$ on a compact subset $K \\subset \\mathbb{R}^d$, the topologies induced by TWD and 1-Wasserstein distances will be identical (i.e. convergence in one implies convergence in the other, and vice versa).  This is because the tilting function $g(x)$ can be viewed as importance weights, so that instead of comparing two measures $p(x)$ and $q(x)$ directly we compare instead their importance re-weighted measures, proportional to $g(x)p(x)$ and $g(x)q(x)$, and it is clear that convergence of either implies convergence of the other under compactness.  However, if $g(x)$ can approach either $0$ or $\\infty$, which is what happens in our work concerning measures on $\\mathbb{R}^d$, then the topologies of TWD and 1-Wasserstein distances do not coincide in general.  \n\nOn the other hand, we note that TWD induces a much weaker topology than, for example, divergences such as Kullback--Leibler or Hellinger, since it does not require absolute continuity of measures.  Further, we emphasise that Theorem 2 (convergence control) is really the main result of our manuscript, as opposed to Theorem 1 (convergence detection).  This is because convergence control justifies the design of algorithms that seek to minimise GF-KSD, guaranteeing the consistency of the approximations that are generated.\n\nBased on your useful feedback, we commit to expand our discussion of TWD in the manuscript to include the above points, which should broaden the accessibility of our manuscript.\n\n> I could not understand the importance of the proposed method because I am not sure for what problems the proposed method is effective [...] I wondered if GF-KSD is really flexible enough to generate samples for complex real data.\n\nThis is an excellent question; our manuscript takes only a first step toward answering it, establishing the theoretical foundations of Gradient Free KSD (GF-KSD), enumerating its possible failure modes, and providing positive proofs of concept.  Subsequent follow-up work will be required to understand in detail the scenarios where GF-KSD can be effective and when it will encounter failure modes.  We will be sure to emphasise the need for a subsequent detailed empirical investigation in the conclusion section of our revised manuscript.\n\n> What do the dotted and solid lines in Figure 1(b) correspond to?\n\nIn the caption of Figure 1 we state that \"the colour and style of each curve in (b) indicates which of the sequences in (a) is being considered\", but we will rephrase this to be explicit that the (e.g.) dashed blue curves in (b) correspond to the sequence of approximating distributions displayed as dashed blue curves in (a).\n\n> Is it required to adjust parameters of the Laplace distribution, KDE, and GMM for $q$ in some way? If so, what is the recommended method?\n\nPlease allow us to emphasise that KDE and GMM are not part of any method that we propose, we simply used these to generate some examples of distributions $q$, for the purpose of illustration in Figure 1.  We will add additional emphasis on this point in the manuscript.  The Laplace approximation has no degrees of freedom to be specified.\n\n> Looking at Figure 3 (b), it seems that the number of samples ($n$) must be very large ($\\log n = 5$, i.e., $n = 150$) even for low-dimensional problems such as $d=8$ in order for there to be any difference in energy distance. Is my understanding correct?\n\nThat is completely correct, but we would argue that $n = 150$ is not a \"very large\" number of samples if we aim to accurately represent an 8-dimensional target.  Of course, \"very large\" will be application-specific, but if we were to take a Cartesian product of a grid of size just 3 in each dimension then we would need $n = 3^8 = 6561$ samples, which is much greater than $n = 150$.\n\nThank you once again for your constructive suggestions on how our manuscript can be improved; we hope that our commitment to resolve them will be reflected in your updated reviewer scores."}, {"review_id": "LhsvWmDL1e", "reviewer": "Reviewer_DoGh", "summary": "This paper explores the use of Stein discrepancies in scenarios where the score function of the target distribution is unavailable or computationally impractical to evaluate. The authors introduce a novel approach called gradient-free kernelized Stein discrepancy (GF-KSD), which leverages a Stein operator developed by Han and Liu (2018) that does not rely on the score function. The authors establish sufficient conditions for the resulting divergence to control and detect convergence. The empirical evaluation of this divergence extends its application to Stein importance sampling and Stein variational inference, surpassing the scope of the previous work by Han and Liu (2018).", "questions": "It would be helpful to elaborate on the numerical instability issue noted in Section 4.2. Does this issue only occur when using GF-KSD for Stein variational inference, or does it also occur in the experiments in Section 4.1? Is that an artefact of the specific form of the target densities chosen for this experiment? Can this be avoided by a judicious choice of $q$?", "limitations": "Examples where the proposed approach can fail are extensively discussed in Section 3.3.", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "**Theories**: Whilst the gradient-free Stein operator has been previously discussed in the community (Han and Liu, 2018), its theoretical properties have not received sufficient attention. This paper established conditions under which the GF-KSD can detect and control convergence of a sequence of empirical distributions to a target probability measure (Theorem 1 and 2). To my knowledge, this contribution fills an important gap in the existing literature and significantly enhances our understanding of the topic.\n\n**Discussions**: Key results on both the theoretical and empirical sides are sufficiently discussed. Limitations of the GF-KSD are thoroughly examined and supported by empirical evidence (Section 3.3). Discussions on the choice of the degree of freedom, the density, are also included (Section 3.1).\n\n**Structure**: This paper exhibits excellent writing with clear motivations throughout.\n\nJ. Han and Q. Liu. Stein variational gradient descent without gradient. In Proceedings of the 35th International Conference on Machine Learning, pages 1900–1908. PMLR, 2018.", "weaknesses": "**Experiments**: My only major concern is over the empirical results. The problems examined in the experiments are primarily toy examples with relatively low dimensionalities. where the one with the highest dimension is a 8-dimensional inference problem for a Lotka-Volterra model. The highest-dimensional problem considered is an 8-dimensional inference problem for a Lotka-Volterra model. To further validate the applicability of the proposed method in more realistic scenarios, additional empirical evidence on higher-dimensional problems would be beneficial.\n\nSpecifically, it would be interesting to investigate the performance of the default choice of Laplace approximation when the dimensionality of the target distribution is high. This analysis could shed light on whether the default approach continues to yield satisfactory results in non-toy scenarios.\n\nMoreover, the reported numerical instability issue in Section 4.2 warrants further investigation, particularly in high-dimensional problems. Understanding the extent to which this instability becomes prominent in higher dimensions is crucial as it may impact the practical viability of the GF-KSD approach.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your kind comments and for your eloquent assessment of our manuscript.  As the reviewer most familiar with the literature on KSD, we sincerely hope that you will champion our work for acceptance, in light of the rather mixed scores we have received.\n\n> To further validate the applicability of the proposed method in more realistic scenarios, additional empirical evidence on higher-dimensional problems would be beneficial.\n\nWe completely agree; as we mentioned to Reviewer 9g9f, our manuscript is limited to establishing the theoretical foundations of Gradient Free KSD (GF-KSD), enumerating its possible failure modes, and providing positive proofs of concept.  Subsequent follow-up work will be required to understand in detail the scenarios where GF-KSD can be effective and when it will encounter failure modes.  We will be sure to emphasise the need for a subsequent detailed empirical investigation in the conclusion section of our revised manuscript, including to problems that are higher-dimensional.  \n\n> It would be helpful to elaborate on the numerical instability issue noted in Section 4.2. Does this issue only occur when using GF-KSD for Stein variational inference, or does it also occur in the experiments in Section 4.1? Is that an artefact of the specific form of the target densities chosen for this experiment? Can this be avoided by a judicious choice of $q$?\n\nThank you for the opportunity to clarify this point:  The numerical issues that arose in Section 4.2 were encountered due to our simultaneous learning of both $p$, the distributional target, and $q$, the axuiliary distribution used to set up GF-KSD.  At early stages along the optimisation path, where the distribution $q$ is not yet an accurate reasonable of $p$, extreme values of the ratios $q/p$ were encountered and gradient clipping was required to regularise the stochastic gradient descent.  This problem was not present in Section 4.1, where $q$ was fixed as a Laplace approximation to $p$, such that extreme values of $q/p$ were not encountered.  These results are encouraging; GF-KSD is numerically stable in the context of Stein Importance Sampling when the Laplace approximation is used.  The use of GF-KSD in Stein Variational Inference requires further investigation; for example, there is at present no theoretical justification for the joint learning of $q$ along the optimisation path - we simply report that this numerical strategy was found to work well and represents \"a promising avenue for further\nresearch\".\n"}, {"review_id": "Tp5bNF9BCJ", "reviewer": "Reviewer_WN2C", "summary": "The paper provides a method to estimate the kernel Stein discrepancy without gradients. ", "questions": "How limited is the class of distributions that allow gradient free estimation of the kernel Stein discrepancy?\n\nAnd the paper probably cherry picked good experiments, how often do the mentioned failure modes occur?", "limitations": "The paper is not clearly written for someone who is not familiar with the field. After reading, I'm still confused about whether someone else has published a similar method before. ", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper provides a method to estimate the kernel Stein discrepancy without gradients. ", "weaknesses": "The paper is not clearly written for someone who is not familiar with the field. After reading, I'm still confused about whether someone else has published a similar method before. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking a look at our manuscript.\n\nPlease see our global response to all Reviewers on the topic of how our manuscript has been presented.\n\n> After reading, I'm still confused about whether someone else has published a similar method before.\n\nThank you for the opportunity to clarify this point.  A gradient-free Stein operator had been described in the context of Stein variational gradient descent (Han and Liu, 2018), but a Stein _discrepancy_ based on this operator had not been proposed or analysed.  Our work contributes the first such gradient-free Stein discrepancy, which is based on a reproducing kernel Hilbert space Stein set.  Reviewer DoGh has summarised the situation well, explaining that \"whilst the gradient-free Stein operator has been previously discussed in the community (Han and Liu, 2018), its theoretical properties have not received sufficient attention. This paper established conditions under which the GF-KSD can detect and control convergence of a sequence of empirical distributions to a target probability measure (Theorem 1 and 2). To my knowledge, this contribution fills an important gap in the existing literature and significantly enhances our understanding of the topic\".\n\n> How limited is the class of distributions that allow gradient free estimation of the kernel Stein discrepancy?\n\nFor concreteness, we have interpreted this question as asking \"for which $p$ does gradient-free kernel Stein discrepancy (KSD) offer convergence control?\", since convergence control is the principal theoretical justification for using gradient-free KSD (GF-KSD) as an objective to be minimised.  Please let us know in the discussion if this was not what was being asked.\n\nIf we choose $q = p$ we obtain the standard KSD, and in this case the class of $p$ for which KSD provided convergence control is known to include all distributions $p(x)$ on $\\mathbb{R}^d$ with Lipschitz continuous $\\nabla \\log p(x)$ and _distantly dissipative_ tails (see Section 2.2).  The latter condition is implied when $p$ is strongly log-concave outside a compact set (i.e. sub-Gaussian tails).  For $q \\neq p$, the distant dissipativity condition applies to $q$ rather than $p$, but then our $\\inf_{x \\in \\mathbb{R}^d} q(x) / p(x) > 0$ condition in Theorem 2 implies that $p$ must also have sub-Gaussian tails in order for the preconditions of Theorem 2 to hold.  Thus, GF-KSD applies to distributions $p$ with sub-Gaussian tails -- this is more general than standard KSD since the requirement for $p$ to be distantly dissipative is removed.  We will make sure to emphasise this point in the revised manuscript.\n\nThank you for your feedback; we hope that our clarifications, our commitment to improve the manuscript, and the remarks of Reviewer DoGh will have given you a more positive impression of our work."}]

Author rebuttal:
__On the Presentation of our Manuscript__

Our manuscript is theoretical in nature and is written for readers who have some familiarity with the concept of kernel Stein discrepancy (KSD).  This is a small but important subsection of the NeurIPS community, who we believe will be the readers most interested in this work (e.g. KSD is currently being used in the community for applications such as goodness-of-fit testing, sampling, gradient estimation, and variational inference).  Excellent introductions to KSD already exist, such as "A Short Introduction to Kernelized Stein Discrepancy" on Qiang Liu's website, or "Measuring Sample Quality with Kernels" by Gorham and Mackey (ICML 2017); there is also now a Wikipedia page on "Stein Discrepancy".  For readers with this background, we believe our mansucript strikes an appropriate balance between precision and concision.  For example, Reviewer DoGh comments that "this paper exhibits excellent writing with clear motivations throughout".

Meta-review: This is paper makes nice contribution to the KSD literature by developing the properties of a KSD that can be applied when it is impractical to compute model log density gradients. However, in addition to the points brought up by the reviewers, the authors should address a two other points in their camera-ready version:

1. Clarify that they did not "introduce" this gradient-free KSD. Rather, it was first proposed in Han and Liu [2018].
2. Clarify that the gradient-free KSD is essentially an importance-sampled KSD using p as the importance distribution.
Decision: Accept (poster)