The Polynomial Stein Discrepancy for Assessing Moment Convergence
Authors: Narayan Srinivasan
Matthew Sutton
Christopher Drovandi
Leah F South
Venue: ICML 2025 poster (2025)
Keywords: Stein Discrepancy
Approximate MCMC
Goodness-of-fit
Primary area: probabilistic_methods->monte_carlo_and_sampling_methods
TL;DR: Scalable linear-time alternative to the kernel Stein discrepancy for measuring sample quality

Abstract: We propose a novel method for measuring the discrepancy between a set of samples and a desired posterior distribution for Bayesian inference. Classical methods for assessing sample quality like the effective sample size are not appropriate for scalable Bayesian sampling algorithms, such as stochastic gradient Langevin dynamics, that are asymptotically biased. Instead, the gold standard is to use the kernel Stein Discrepancy (KSD), which is itself not scalable given its quadratic cost in the number of samples. The KSD and its faster extensions also typically suffer from the curse-of-dimensionality and can require extensive tuning. To address these limitations, we develop the polynomial Stein discrepancy (PSD) and an associated goodness-of-fit test. While the new test is not fully convergence-determining, we prove that it detects differences in the first $r$ moments for Gaussian targets. We empirically show that the test has higher power than its competitors in several examples, and at a lower computational cost. Finally, we demonstrate that the PSD can assist practitioners to select hyper-parameters of Bayesian sampling algorithms more efficiently than competitors.

The Polynomial Stein Discrepancy for Assessing Moment Convergence
Narayan Srinivasan 12
Matthew Sutton Christopher Drovandi 12
Leah South
Abstract
target distribution P. These discrepancies can be applied when P iis only known up to normalising constant, making them particularly appealing in the context of Bayesian in- ference. Unlike classical tools for assessing sample quality in Bayesian inference, Stein discrepancies are applicable to hyperparameter tuning (Gorham & Mackey, 2017) and goodness-of -faesss--fffi-- testing (Chwialkowski et al., 2016; Liu et al., 2016) for biased Monte Carlo algorithms like stochastic gra- dient Markov chain Monte Carlo (SG-MCMC, Welling & Teh, 2011; 1; Nemeth & Fearnhead, 2021).
We propose novel method for measuring the discrepancy between set of samples and de- sired posterior distribution for Bayesian inference. Classical methods for assessing sample quality like the effective sample size are not appropriate for scalable Bayesian sampling algorithms,, such as stochastic gradient Langevin dynamics, that are asymptotically biased. Instead, the gold stan- dard is the kernel Stein Discrepancy (KSD), which is itself not scalable given its quadratic cost in the number of samples. The KSD and its faster extensions also typically suffer from the curse-of- dimensionality and can require extensive tuning. To address these limitations, we develop the poly- nomial Stein discrepancy (PSD) and an associated goodness-of-fi test. While the new test is not fully convergence-determiig
Currently, the most widely adopted Stein discrepancy is the kernel Stein discrepancy (KSD, Gorham & Mackey, 2017; Chwialkowski et al., 2016; Liu et al., 2016), which applies a so-called Stein operator (Stein, 1972) ducing kernel Hilbert space. The KSD has an analytically tractable form and can perform well in tuning and goodness- of-fit tasks using the recommended inverse multiquadric (IMQ) base kernel. However, the computational cost of KSD is quadratic in the number of samples, n, from the distribution Q, so it is infeasible to directly apply KSD in applications where large numbers of MCMC iterations have been used, To address this bottleneck, Liu et al. (2016) introduced early linear-time alternative. However, the statistic suffers from poor statistical power compared to its quadratic-time counterpart, meaning that the probabil- ity of rejecting the null hypothesis when it is false is low. Several approaches have been disseminated in the litera- ture to tackle these problems and reduce the cost of using Stein discrepancies. Notably, the finite set Stein discrepancy (FSSD, Jitkrittum et al., 2017) and the random feature Stein discrepancy (RFSD, Huggins & Mackey, 2018) have been proposed as linear-time alternatives to KSD.
repro-
prove that it detects differences in the first moments for Gaussian tar- gets. We empirically show that the test has higher power than its competitors in several examples, and at a lower computational cost. Finally, we demonstrate that the PSD can assist practitioners to select hyper-paranne papere-praetters of Bayesian sampling algorithms more efficiently than competitors.
1. Introduction
Reliable assessment of posterior approximations was re- cently named one of three "grand challenges in Bayesian computation" (Bhattacharya et al., 2024). Stein discrep- ancies (Gorham & Mackey, 2015) can assist with this by providing ways to assess whether distribution Q, known via samples (2r00 p Q, is a good approximation to
The main idea behind FSSD is to use a finite set of test locations to evaluate the Stein witness function, which can be computed in linear-time. The FSSD test can effectively capture the differences between P and Q, by optimizing the test locations and kernel bandwidth. However, the test is sensitive to the test locations and other tuning parameters, all of which need to be optimized (Jitkrittum et al., 2017). The test requires that the samples be split for this optimisation and for evaluation. Additionally, the FSSD experiences a degradation of power relative to the standard KSD high dimensions (Huggins Mackey, 2018).
'School of Mathematical Sciences, Queensland Uni- versity of Technology (QUT), Brisbane, Australia 'QUTT Centre for Data Science, Brisbane. Australia "School of Mathematics and Physics, University of Queensland, Bris- bane, Australia. Correspondence to: Narayan Srini- vasan <narryansesaadrriiiiiiiaaaa &ou.hh qu..
www.drpu..cuu..uuu. Leah South
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
Random Fourier features (RFF) (Rahimi & Recht, 2007) a well-established technique to accelerate kernel-based methods. However, it is known (Chwialkowski al., 2015, Proposition 1) that the resulting statistic fails to distinguish a large class of probability measures. Huggins & Mackey (2018) alleviate this by generalising the RFF approach with their RFSD method, which is hear-linear-time. The authors identify a family of gnnvereenceddetermniig discrepancy measures that can be accurately approximated with impor- tance sampling. Like FSSD, tuning can be challenging for RESD because of the large number of example-seecific choices, including the optimal feature map, the importance sampling distribution, and the number of importance sam- ples. We have also found empirically that the power of the goodness-of-fit test based on RFSD is reduced when direct sampling from P infeasible.
and provides background on Stein discrepancies, KSD and goodness-of-ft tests. The PSD and goodness-of-ft tests based on the same are presented in Section 3, along with theoretical results about detecting moment discrepancies and the asymptotic power of the test. Section 4 contains simulation studies performed on benchmark examples. The paper is concluded in Section 5.
2. Background
Let the probability measure Q, known through the n sam- ples, {x(() H be supported in X. Suppose the target dis- tribution P has a corresponding probability density function p.
Utilising the notion of integral probability metries (IPMs, Müller, 1997), discrepancy between P Pnd Q can be de- fined
In addition to these limitations, KSD and current linear- time approximations can fail to detect moment convergence (Kanagawa et al., 2022). The standard Langevin KSD us- ing the IMQ kernel controls the bounded-Lipschitz metric, which determines weak convergence, but fails to control the convergence in moments. This significant draw- back because moments are often the main expectations of interest and, for many biased MCMC algorithms, they are where bias is likely to appear, as explained in Section 3.3. Kanagawa et al. (2022) propose an extended, quadratic-tima KSD that is able to control convergence in moments but this method is slower than KSD with an IMQ kernel.
d P. Q) == sup |Ex-r - [x()]] (1)
hEW
Various IPMs correspond to different choices of function class H (Anastasiou et al., 2023). However, since we only have access to the unnormalized density P, the first expec- tation in (I) is typically intractable.
2.1. Stein Discrepancies
Stein discrepancies make use of a so-called Stein operator (A) and an associated class of functions, ((A), such that
Motivated by the shortcomings of KSD and linear-time KSD methods, we propose linear-time variant of the KSD named the polynomial Stein Discrepancy (PSD). The method we propose in this article detects discrepancies in moments while still being computable in linear time. This approach, based on the use of r th order polynomials, is motivated by the zero variance control variates (Assaraf & Caffarel, 1999; Mira et al., 2013) used in variance re- duction of Monte Carlo estimates. While PSD is not fully convergenee-detrriinngg we show that when P is Gaussian, which is a reasonable approximation for big data applica- tions (Bardenet et al., 2017), the discrepancy is zero if and only if the first moments of P and Q match. We empiri- cally show that PSD has good statistical power for detecting discrepancies in moments, including in applications with non-Gaussian P. and we also demonstrate its effectiveness for tuning SG-MCMC algorithms. Importantly, the decrease in power with increasing dimension is considerably less than competitors, and the method requires no calibration beyond the choice of the polynomial order, which has a clear inter- pretation. Liu & Wang (2018) demonstrate the theoretical advantages of using similar kernels for moment matching in Stein variational gradient descent.
Ex-r (AA) = for all E G(A).
(2)
This allows us to define the so-called Stein discrepancy
S(Q.A.9) sup || [[ |l. (3)
SEG(A)
There is flexibility in the choice of Stein operator and we fo- cus operators that are appropriate when x E X C Rdd By considering the generator method of Barbour (1990) applied to overdamped Langevin diffusion (Roberts & Tweedie, 1996), one arrives (Gorham & Mackey, 2015) at the second order Langevin- -Stein operator defined for real-valued g as
A g (zg(x) + Veg(z) V logp(z).
(4)
Under mild conditions on and g. Ex P.x-p/la ((A9))
[X)] =
O as required. For example for unconstrained spaces where X' = Rdd if g is twice continuously differen- tiable, log p is continuously differentiable and [[Vg(z)] s C||==|| -p())l~p(x>)>>)))))<)))) for some constant C' > 0 and 0, then Ex-P 9)()] = 0 as required (South et al., 2022).
The paper is organized as follows. Section 2 sets notation
Alternatively, one can consider a first order Langevin- -Seii operator defined for vector-valued g (Gorham & Mackey,
2

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
2015)
and they provide lower-bound for the KSD in terms of the bounded Lipschitz metric.
A = V. g((x) + g(x) Va logp(r).
(5)
The discrepancy S(Q, T,)) can be estimated with its corre- sponding U-statistic (Serfling, 2009) as follows
Under regularity conditions on g and p similar to those for (4) we get that, Bp[(AE (p](AA9))(X)) (x)] = 0.
1 n(n 1) ko (z:(()) ,I()).
KSD
2.2. Kernel Stein Discrepancy
J
The key idea behind KSD is to write the maximum dis- crepancy between the target distribution and the observed sample distribution by considering & in (3) to be functions in reproducing kernel Hilbert space (RKHS, Berlinet & Thomas-Agnan 2004) corresponding to an appropriately chosen kernel.
This is the minimum variance unbiased estimator of KSD (Liu et al., 2016). Alternatively, one could consider the V-Statistic proposed in Gorham & Mackey (2017) and Chwialkowski et al. (2016), given by
2 = KSD
1 π
8o(x(0),
(8)
A symmetric, positive-defintte kernel A: Rdd x Rd + IR induces an RKHS K₂ of functions from T2 - RR. For any r R2d, k(x,-) E Kg. From the reproducing property of the RKHS, ff E KR then fx) = (f. -)x
The V-statistic, while no longer unbiased, is strictly non- negative and hence can be better suited as a discrepancy metric when compared the U-statistic.
The Stein set of functions 9 with the associated kernel is taken to be the set of vector-valued functions g, such that each component 9j belongs to K2 and the vector of their norms llg belongs to the unit ball, i.e. 9 == (g=
2.3. Goodness-of Fit Testing
In goodness-of testing, the objective is to test the null hypothesis Hp: Q P against an altermative hypothesis, typically that H₁: Q + P. Existing goodness-of-ff tests are based either asymptotic distributions of U-statisties (Liu et 2016; Jitkrittum et al., 2017; Huggins & Mackey, 2018) or bootstrapping (Liu et al., 2016; Chwialkowski et al., 2016).
(g.....,ga) 9a) : ||0|| < 1 for by = |1,
Under certain regularity conditions enforced on the choice of the kernel k. Gorham & Mackey (Proposition 2, 2017) and also Liu et al. (2016); Chwialkowski et al. (2016) arrive at the closed form representation for the Stein discrepancy given in (3) for this particular choice of Stein set and label
it the KSD
Asymptotic goodness-of-fi tests cannot be implemented di- rectly for KSD or its early linear-time alternatives (Liu et al., 2016; Chwialkowski al., 2016). Instead, bootstrapping is used to estimate the distribution of the test statistic under Ho. Chwialkowski et al. (2016) develop test for poten- tially correlated samples (a(() J using the wild bootstrap procedure (Leucht & Neumann, 2013). Liu et al. (2016) use the bootstrap procedure of Hušková & Janssen (1993); Arcones & Gine ( (1992) for degenerate U-statistics, in the setting where the samples are uncorrelated. As the number of samples and the number of bootstrap samples go to infinity, Liu et al. (2016) and Chwialkowski et al. (2016) show that their bootstrap methods have correct type I error rate (i.e. correct probability of rejecting the null hypothesis when is true) and power one.
KSD = S(Q,A,9) = <XDD==S(.AA1)-//aee-ohdtt.))))
(6)
where in the particular case of the first order operator (5),
ko(a, x²) Va Vyk(o, x') + v=)(x=))) u(x')
+ V V2<k(z,z') x') u(z) + k(2, a )u(x) u(z²).
Here, u(x)= V log p(x) and k(xr,x') is the chosen ker- nel.
Gorham & Mackey (Theorem 6, 2017) show that KSDs based on common kernel choices like the Gaussian ker- nel <', = cxp(- 2 ||r - x' (l3) and the Matérn ker- nel fail to detect non-convergence for d 2 3. Gaussian kernels are also known experience rapid decay in statis- tical power in increasing dimensions for common errors (Gorham & Mackey, 2017). As alternative, Gorham & Mackey (2017) recommend the IMQ kernel k(xx,²') = (c² + || - z' 8 with e = 1 and / = -0.5. They show that IMQ KSD detects convergence and non-convergence for c >>0 0 and 3 E (-1,0) 0) for the class of distantly dissipative P. with Lipschitz log p and Eg-r[|/V (x))()))] < co,
More recent linear-time alternatives, specifically the FSSD and RFSD, employ tests based the asymptotic distribu- tion of U-statistics. These methods have similar thcoretical guarantees as + oc and the number m of simulations from the tractable null distribution go to infinity but they
'Lu et al. (2016) do not appear to directly state the requirement that n oc in their theorem, but this follows from the proof on which is based (Hušková & Janssen, 1993).
3

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
can suffer from poor performance in practice. These tests require estimate for a covariance under P. When di- rect sampling from P is infeasible, which is typically the case when measuring sample quality for Bayesian inference, these methods estimate this covariance using samples from Q. While this is asymptotically correct as n - oo (Jitkrit- tum et al., 2017; Huggins& & Mackey, 2018), we have found empirically that the use of samples from Q can substan- tially reduce the statistical power of the tests, particularly for RFSD.
where E (x]A) This derivation is provided in more detail in Appendix A. For the sample {ax(0) j= from Q, we have ū = +o2 1 (a) Note that PSD implicitly depends on a polynomial order r through J. A simple frrm for (P)(2)) is available in Appendix A of South et al. (2023).
The squared finar-time solution for PSD can also be ex- pressed as V-statistic. Specifically
PSD
2
3. Polynomial Stein Discrepancy
Motivated by the practical effectiveness of polynomial fune- tions in MCMC variance reduction and post-processing (Mira et al., 2013; Assaraf & Caffarel, 1999), we develop PSD linear-time alternative to KSD.
1
A (z 2000)
n²
1=1
j=1
where A(x,)) (())T((V) and [T(]]] = (Px(
Ek= A.P(0),,)
3.1. Formulation
The U-statistic version of the squared PSD, which will be helpful in goodness-of fit testing, is
This section presents an intuitive and straightforward ap- proach to derive the PSD. The corresponding goodness-of-fit tests are presented in Section 3.2.
PSD
1
a
T3
(1) 2000))
Consider the class of rth order polynomials. That is, let E No. Σ 0g S r), where [i] denotes the äth dimension of x. Intuitively, 9 is the monomial terms that we will simply denote by (x) for = I,....J.. J. This is a valid Stein set in [(x)] 0 for all g E G under mild conditions depending on the distribution P. For example, when the density of P. p(x), is supported on an unbounded set then a sufficient condition is that the tails of P decay faster than an rth order polynomial. One can also consider boundary conditions for bounded spaces, as per Proposition 2 of Mira et al. (2013).
n(n
1)
G ( =1[i]=
j
1
J
span of J (*+")
nnn 1)
n1
71
%k
,
(10)
k=1
k=1
that [A (2)
where 2/ 1/ Ei (APa (2
The computational complexity of this discrepancy is (nJ)) = On (((())")) In very high dimensions, practi- tioners concerned about computational cost could run an approximate version of the test for ((ndr) by excluding in- teraction terms from the polynomial, for example for r = 2, monomial terms z would only be included for i = j. In the case of a Gaussian P with diagonal covariance matrix, such a discrepancy would detect differences in the marginal moments.
Henceforth we denote the second order operator A! as simply A. Using the linearity of the Stein operator (4), the aim is to optimize over different choices / with real coefficients 3y. for h: = 1,2...., J in (3). Analogous to the optimization in KSD, the optimization is constrained over the unit ball, that is || < 1. The result is
We note that while it would be possible to implement KSD with the conventional polynomial kernel k(æ, y) (1 + xx" Ty)", such discrepancies have not yet been exten- sively explored in the literature. We show the promise of polynomial kernels in terms of statistical power and linear- time complexity. Our PSD also differs from what one would obtain with a conventional polynomial kernel, offering a simpler formulation that may be more effective for identify- ing specific moments where discrepancies occur; in the case of a Gaussian P with independent components, monomial terms correspond directly to multi-index moments. We also present novel theory (Proposition 3.2) that enhances the (9) understanding and interpretability of PSD in the Bayesian big data limit.
PSD- sup (x()l( SEG
sup
EQ
Bx (x)]]
Be (x)]
sup
=
k=1
4

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
3.2. Goodness-of-fi Test
bootstrap procedure suggested by Liu et al. (2016), Arcones & Gine (1992), Hušková & Janssen (1993). This bootstrap test will itself be linear-time; such method was not avail- able in the linear- time KSD alternatives of Jitkrittum et al. (2017) and Huggins & Mackey (2018).
For the goodness-of-fit test based on PSD, we will be test- ing the null hypothesis Ho: Q P against a composite or directional alternative hypothesis. The form of the alter- native hypothesis depends on P. but we show in Section 3.3 that for Gaussian P(a suitable assumption under the Bayesian big data/Bernstein Von-Mises limit), H₃ is that the first r moments of Q do not match the first moments of P. Achieving high statistical power for these moments is important, as described in Section 3.3.
In each bootstrap replicate of this test, one draws weights (W1, w2..... w., Multinomial(n: 1 .....: A)) and com-
putes the bootstrap statistic
n
(wg 1 ( 2010) n
PSD u. boot
(w)
Observe that (10) is degenerate U-statistic, so asymp- totic test can be developed using [r()]] = (z).
1=1 j=1
n2
241
Corollary 3.1 (Jitkrittum et al. (2017)). Let Zi..... ZJ be i.idd random variables with Z4_ ~ NN 1). Let == [xzg]T(z)] and Σ, == [T))] E R for r E [P,Q}. Let {wili be the eigenvalues of the co- variance matrix Ep = Ex<plT(z)+ (z)]. Assume that
(
(II)
(w, n )
E B A2 (2r v) < co. Then, the following statements hold
If the observed test statistic, (10), is larger than the 100(1 a) percentile of the bootstrap statistics, (11), then Ho is rejected. The consistency of this test follows from Theorem 4.3 of Liu et al. (2016) and from Hušková & Janssen (1993), We recommend the bootstrap test because we empirically find that it has higher statistical power than the asymptotic test using samples from Q.
1. Under Ho Q P.
J
nPSD d
2
Σ t=1
[23
1)ws.
2. Under H1 (Q # P in way specified by PSD), if
The aforementioned tests are suitable for independent sam- ples; in the case of correlated samples, an alternate linear- time bootstrap procedure for PSD can be developed follow- ing the method of Chwialkowski et al. (2016).
oin == 4u T Eight > 0, then
PPDD 2 - PPDD*) 4, o
Proof. This follows immediately from Proposition 2 of Jitkrittum et al. (2017), which itself follows from Liu et al.
3.3. Convergence of Moments
(Theorem 4.1, 2016) and Serfling (Chapter 5.5, 2009). L
Since are using a finite-dimensional (i.e. non- characteristic) kernel, PSD is not fully convergence- determining. We do not view this as major disadvantage. Rather than aiming to detect all possible discrepancies be- tween P and Q and doing so with low statistical power, the method designed to achieve high statistical power when the discrepancy is in one of the first moments.
simple approach to testing Hg: Q = would then be to simulate from the stated distribution under the null hypothesis m times and reject Ho at a significance level of a if the observed test statistic (10) is above the Ta = 100(1 a) percentile of the m samples. How- ever, simulating from this distribution can be impractical since Σn requires samples from the target P. which is typ- ically not feasible. Following the suggestion of Jitkrittum et al. (2017), Ep can be replaced with the plug-in estimate computed from the sample, S.g == 1 D (I:)
Often, moments of the posterior distribution, such as the mean and variance, are the main expectations of inter- est. These moments can also be where differences are most likely to appear; for example, the posterior variance is asymptotically over-estimated in SG-MCMC methods (Nemeth & Fearnhead, 2021). Tuning SG-MCMC amounts to selecting the right balance between a large step-size and small step-size. Large step-sizes lead to over-estimation (bias) in the posterior variance, while small step-sizes may not sufficiently explore the space for a fixed n, and thus may lead to unttrnestimated posterior variance. Hence, it is critical to assess the performance of these methods in estimated second-order moments.
1 E.= T<x))] 3. tum et al. (2017) ensures that replacing the covariance ma- trix with E,, still renders consistent test.
T22,) Theorem 3 of Jitkrit-
In practice, however, there may be unnecessary degradation of power, especially in high dimensions due to the approx- imation of the covariance in Corollary 3.1. To circumvent this, we propose an alternative test which is follow the
5

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
Proposition 3.2. Pf is Gaussian with symmetric positive- definite covariance matrix Σ,, PSD> = 0 if and only if the multi-index moments of and Q match up to order
Huggins & Mackey (2018), we use the L1 IMQ base kernel and fix the number of features to 10.
The simulations are run using the settings and implementa- tions provided by the respective authors, with the exception that we sample from Q for all asymptotic methods since sampling from P is rarely feasible in practice. Goodness- of-fit testing results for PSD in the main paper are with the bootstrap test, which we recommend in general. Results for the PSD asymptotic test with samples from Q are shown in Appendix C. Following Jitkrittum al. (2017) and Hug- gins & Mackey (2018), our bootstrap implementations for KSD and PSD use V-statistics with Rademacher resampling. The performance is similar to the bootstrap described in Liu et al. (2016) and in Section 3.2.
The proof of Proposition 3.2 is provided in Appendix B. As consequence of Proposition 3.2, we have the following result:
Corollary 3.3. Suppose the conditions in Corollary 3.1 hold. Then, in the Bernstein-von Mises limit (i.e. the Bayesian big data limit), the asymptotic and bootstrap tests have power + for detecting discrepancies in the first moments of P and Q as n + DC.
Proof. This result follows from the consistency of these tests and the conditions under which PSD 0 as per Propo- sition 3.2.
Code to reproduce these results available at https: / //ithub. www.coo////ooooo haam/Nars98 /tt/s///PSDDD This code builds on ex- isting code (Huggins, 2018; Jitkrittum, 2019) by adding PSD as new method. All experiments were on high performance computing cluster, using a single core for each individual hypothesis test.
This is an important result given that biased MCMC al- gorithms are often used for big-data applications, with subsampling- based methods arguably being the most com- mon application of KSD. In this context, P is often close to Gaussian because typically the "Bernstein-von Mises ap- proximation of the target posterior distribution is excellent" (Bardenet et al., 2017). We also empirically show good per- formance for detecting discrepancies in moments in Section 4 and in supplementary results in the Appendices.
Further empirical investigations are available in Appendix C. In Appendix C.2, we show that PSD) with r is tracking second order moments well for the second order moment example of Kanagawa al. (2022). We also provide two logistic regression examples in Appendices C.3 and C.4, where we demonstrate that PSD can tracks discrepancies in moments for realistic Bayesian inference tasks. The for- mer is an example where P is far from normality, while the Bernstein-von Mises approximation reasonable for the latter. We investigate performance of PSD without interac- tions in Appendix C.5 and we find that the performance is remarkably similar to PSD with interactions. Finally, we provide an extreme example where moments are not well approximated and show what expectations are being tracked instead in Appendix C.6.
4. Experiments
This section demonstrates the performance of PSD on the current benchmark examples from Liu et al. (2016), Chwialkowski et al. (2016). Jitkrittum et al. (2017) and Hug- gins & Mackey (2018). The proposed PSD is compared to existing methods on the basis of runtime, power in goodness- of-fi testing and performance as a sample quality measure. The following methods are the competitors:
IMQ KSD: standard, quadratic time KSD (Gorham & Mackey, 2017; Liu et al., 2016; Chwialkowski et al., 2016) using the recommended IMQ kernel with e = and 8 = -0.5.
4.1. Goodness-of-Fit Tests
Following standard practice for assessing Stein goodness- of-fit tests, we begin by considering P = N NNOa. Ta) and assessing the performance for a variety of Q and using statistical tests with significance level a 0.05. We use m = 500 bootstrap samples estimate the rejection thresh- old for the PSD and KSD tests.
Gauss KSD: standard, quadratic time KSD using the common Gaussian kernel with bandwidth selected us- ing the median heuristic.
We investigate four cases: (a) type I error rate: Q N Ta) (b) statistical power for misspecified variance: Q = Aosa, E), where Enj = 0 for + j. Σn = 1.7 and
FSSD: The linear-time FSSD method with optimized test locations (Jitkrittum et al., 2017). We consider the optimized test locations (FSSID-opt), optimized according to power proxy detailed in Jitkrittum et al. (2017). We set the number of test locations to 10.
>n, = 1 for i = 2, ..,.. d, (c) statistical power for misspeci- fied kurtosis: Q == T(0, 5), a standard multivariate student-t distribution with degrees of freedom, and (d) statistical power for misspecified kurtosis: q(x) = Ie_, Lap(z, I 0, (1 the product of d independent Laplace distributions
RFSD: The near-linear-time RESD method (Huggins & Mackey, 2018). Following recommendations by
6

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
0.20
rejection rat
1.0
Full rejection rat
1.0
null rejection rat
1.0
RFSD
Gauss KSD
null rejection
0.15
IMQ KSD
0.10
0.5
0.5
Gauss FSSD-opt
0.05
PSD r1
PSD r2
0.00
10
21
0.0
10
20
0.0
10
20
0.0
10
20
PSD r3
dimension D
dimension D
dimension D
dimension
PSD 14
(a) Unit Gaussian (i.e. Q = (b)
Variance-perturtedd (c) Student-l, 250 repeats
(d) Laplace, 500 repeats
P). 500 repeats
Gaussian, 200 repeats
Figure 1. Type error rate (a) and statistical power (b,c,d) for detecting discrepancies between the unit Gaussian P and the sampling distribution Q (see the main text for details).
with variance 1. Following Huggins & Mackey (2018), all experiments use n 1000 except the multivariate t, which uses n = 2000.
Table I. Null rejection rates for testing methods at different pertur- bation levels the RBM example.
PERTURBATION:
0
0.02
0.04
0.06
As seen in Figure la, the type I error rate is generally close to 0.05 even with this finite n. Exceptions include RFSD which has decaying type eπror rate with increasing d, and PSD with r = 4, which has slightly over-inflated type error rate for d = 1.
RFSD
0.00
0.48
0.93
0.98
FSSD-opt
0,05
0.70
0.96
0.99
IMQ KSD
0.08
0.99
1.00
1.00
Gauss KSD
0.08
0.95
1.00
1.00
PSD rl
0.08
0.51
0.96
0.99
PSD r2
0.06
1.00
1.00
1.00
As expected, PSD with r = 1 is incapable of detecting discrepancies with the second (b) or fourth order moments (c,d). Similarly, PSD with r 2 and detecting discrepancies with the fourth moment (c,d). When the polynomial order is at least as high the order of the moment in which there are discrepancies, PSD outperforms linear-time methods and is competitive with quadratic-time KSD methods. Specifically, PSD with r = 4 is the only method to consistently achieve a power of R 1 in Figures Ic and 1d. In Figure 1b, PSD with r 2 and r 4 are the only methods consistently achieve a power of 1, while PSD with = 3 has higher statistical power than the competitors at d = 20. Overall, the new methods have a statsistical power up to double the statistical power of IMQ KSD and four times that of linear-time competitors for = 20.
PSD r3
0.09
0.97
1.00
1.00
3 are incapable of
4.2. Measure of Sample Quality
To demonstrate the advantages of PSD as a measure of dis- crepancy we follow the stochastic gradient Langevin dynam- ies (SGLD) hyper-yaraee pareee--rrmmterrseee selection setup from Gorham Mackey (Section 5.3, 2015). Since no Metropolis- Metsonolis-Hastings Hastings correction is used, SGLD with constant step size is bi- ased MCMC algorithm that aims to approximate the true posterior. Importantly, the stationary distribution of SGLD deviates more from the target as grows, leading to an in- flated variance. However, smaller e decreases the mixing speed of SGLD. Hence, an appropriate choice of i critical for accurate posterior estimation.
Next, following Liu et al. (2016) and Jitkrittum et al. (2017), consider the case where the target P is the non-normalized density of a restricted Boltzmann machine (RBM); the samples Qn, are obtained from the same RBM perturbed by independent Gaussian noise with variance o². For x² = 0, Hn : Q = P holds, and for o² > 0 the goal to detect that the n = 1000 samples come from the per- turbed RBM. Similar to the previous goodness-of fit test, the null rejection rate for range of perturbations using 100 repeats is given in Table 1. Notably, PSD with 2 and r = 3 both outperform tnear-time methods and are com- petitive with quadratic-time time methods, but for substantially reduced computational cost. This illustrates that PSD is a potentially valuable tool for goodness-of-fit testing in the non-Ganne on--aussian setting.
Similar to the experiments considered in Gorham & Mackey (2015), Chwialkowski et al. (2016) and Huggins & Mackey (2018), the target P is the bimodal Gaussian mixture model posterior of Welling & Teh (2011), We compare the step size selection made by PSD to that of RESD and IMQ KSD when n 10000 samples are obtained using SGLD. Figure shows the performance of SGLD for a variety of step- sizes in comparison with high quality samples obtained using MALA. Figure 3 shows that PSD with r r= 2, 3 and r = 4 agree with IMQ KSD, selecting e = 0.005 which is visually optimal as per Figure 2. Moreover, when utilized as a measure of discrepancy, PSD is around 70 times faster than KSD and around times faster than RFSD
7

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
0.0005
6 0.005
e=0.05
high high-quality samples
a
-2
-2
-2
-2
0
1
Figure 2. Approximate posterior for mixture example with SGLD for varying step sizes and when sampling from the true posterior using MALA.
discrepancy measure
RFSD
needing efficient, dependable measures of sample quality, especially in the context of complex Bayesian inference applications. For practitioners using PSD for biased MCMC samplers like SG-MCMC, we recommend using PSD with 2.
IMQ KSD
PSD r1
PSD r2
10-1
PSD
PSD M
It well known that probability distributions are completely determined by their moment generating functions, provided they exist (a necessary condition is that all moments are finite). This highlights the main drawback of our method, since we can potentially expect the KSD to outperform PSD for target distributions lacking well-defined moments, For example those with heavy tails such as the Cauchy distri- bution. However, we also note that the standard KSD (and the linear time variants) also fail for Cauchy distributions, due to Theorem 10 of Gorham & Mackey (2015). In such situations, Stein discrepancies based on diffusion Stein op- erators can be considered (Kanagawa et al., 2022). Similar extensions for PSD can be considered in future research
10-3
10"4
10-3
10"2
Figure Step size selection results for SGLD using various meth-
ods.
4.3. Runtime
We now compare the computational cost of computing PSD with that of RFSD and IMQ KSD. Datasets of dimension 10 with the sample size n ranging from 500 to 10000 were generated from P (O,,, Ia). As seen in Figure even for moderate dataset sizes, the PSD and RFSDs computed orders of magnitude faster than KSD. While RFSD is faster than PSD with = 3 or = 4, we have found in Sections 4.1 and 4.2 that PSD can have higher statistical power for detecting discrepancies in moments.
Another situation in which KSD outperforms PSD of order is when the discrepancies lie in moments higher than the rth order moment. For example, KSD outperforms PSD with r = 2, 3 for the student-t and Laplace examples in Figure 1 since the discrepancy is the kurtosis (4th order moment).
time (seconds)
102
RFSD
IMQ KSD
As practical heuristic to determine whether PSD may perform well in detecting moments specifically, one could consider plotting the marginal gradients versus the marginal samples. The closer this is to perfect linear relationship, the more one might expect the moments to be assessed. Further, if the gradicnts are available in analytic form, then one can also determine which expectations are being tracked by investigating the system of linear equations in equation (14) of Appendix B. We do this for the Rosenbrock target in Appendix C.6, and we able to determine exactly which expectations are being tracked. Regardless of whether P is Gaussian, the RBM, SGLD and logistic regression examples show that PSD can a useful measure of sample quality.
10a
PSD r1
PSD r2
10-2
PSD r3
PSD 1
5000
10000
sample size
Figure 4. Runtime for various testing methods where P (Oa,, Ia) with d = 10.
5. Conclusion
Our proposed PSD is powerful measure of sample quality, particularly for detecting discrepancies in moments. The method eliminates the need for extensive tuning, is linear- time and empirically provides high statistical power in high dimensions. This makes it a valuable tool for practitioners
The primary advantage of using the Langevin Stein opera- tor is that, for Gaussian target, the operator preserves the form of the monomials, and consequently we are tracking
8

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
convergence in moments. This may not hold true for other Stein operators. Nevertheless, applying the aforementioned diffusion Stein operators and other Stein operators differ- ent types of polynomials or function classes can be studied further in future work.
References
Anastasiou, A., Barp, A., Briol, F.-X., Ebner, B., Gaunt, R. E., Ghaderinezhad, F., Gorham, J., Gretton, A., Ley, C., Liu. Q.. Mackey, L.. Oates, C. J., Reinert, G., and Swan, Y. Stein's 's method meets computational statistics: A review some recent developments. Statistical Science, 38 (I):120 139, 2023. doi: 10.1214/22-STS865 URL https: /ttps://do... http:///doooog////// 10. 12w112144/22 422 STS863.
We show in Appendix D that the case of Gaussian P (by extension, in the Bernstein-von Mises limit), one can ap- ply any invertible linear transform the parameters and the method will still detect discrepancies between the original Q and P in the first moments. We suggest applying whiten- ing or simply standardization when the variances differ sub- stantially across dimensions. Future work could also con- sider using alternative norms, for example by maximising subject to the constraint that ||@l|| s or ||,3 W Bl12 I 1, the latter of which could be used weight different mono- mials and therefore different moments.
Arcones, M. A. and Gine, E. On the Bootstrap of U and V Statistics. The Annals of Statistics, (02)):55-674, 1992. ISSN 0090-5364.
Assaraf, R. and Caffarel, M. Zero-variance principle for Monte Carlo algorithms. Physical Review Letters, 83(23): 4682-4685, 11 1999. doi: 10.1 10.1103/hhysReeeett.3.....
Since the PSD) is not translation invariant, determining the optimal parameterisation for Stein discrepancies is an open problem and an interesting point for further research. We investigate this further in Appendix C.7.
Barbour, A. D. Stein's method for diffusion approxima- tions. Probability theory and related fields, 84(3):297- 322, 1990.
further extension could be the use of PSD as a tool de- termine the moments in which discrepancies between and could examine and its
Bardenet, R., Doucet, A., and Holmes, C. On Markov chain Monte Carlo methods for tall data, Journal of Machine Learning Research, 18(47):1-33, 2017. URL http: html.
P are occurring. In particular, distribution under Ho to obtain an ordering of which mono- mial terms contribute the most to the discrepancy. Finally, thcoretical investigations into the topological properties re- lating to the convergence of PSD can potential avenue for future research.
/ wtt..o//ppppps/////eeees
Berlinet, A. and Thomas-Agnan, C. Reproducing Kernel Hilbert Spaces in Probability and Statisties. Springer US, Boston, MA, 2004. ISBN 978- 978-1-4613-4792-7 978-1- 4419-9096-9, doi: 10.1007/978 1 4419-9096-999
Acknowledgments
Bhattacharya, A., Linero, A., and Oates, C. J. The ISBA bulletin. International Society for Bayesian Analysis, 31 (3):11-14, 2024.
Narayan Srinivasan is supported by the Australian Govern- ment through a research training program (RTP) stipend. LFS is supported by a Discovery Early Career Re- searcher Award from the Australian Research Council (DE240101190) CD is supported by Future Fellowship from the Australian Research Council (FT210100260). LFS would like to thank Christopher Nemeth and Rista Botha for helpful discussions. Computing resources and services used in this work were provided by the High Performance Com- puting and Research Support Group, Queensland University of Technology, Brishane, Australia.
Chwialkowski, K., Strathmann, H., and Gretton, A. A kernel test of goodness of fit. In Balcan, M. F. and Weinberger, The K. Q. (eds.), Proceedings The 33rd International Con- ference on Machine Learning, volume 48 of Proceed- ings of Machine Learning Research, pp. New York, New York, USA, 20-22 Jun 2016. PMLR. mww..press//lle
URL 1.t..k..
chwialkowskil
Chwialkowski, K. P., Ramdas, A., Sejdinovic, D., and Gretton, A. Fast two-sample testing with analytic representations of probability measures. In Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., and Garneut R. (eds.), Advances in Neural Information Process- ing Systems, volume 28. Curran Associates, Inc., 2015. URL http:::/pproceedingg.nnt neurips. 3S5Tcccal1a9924222eeeaaaa1188113 pdf.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none of which we feel must be specifically highlighted here.

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
George, E. and McCulloch, R. Variable selection via Gibbs sampling. Journal of The American Statistical Associa- tion J AMER STATIST ASSN, 88:881-889, 09 1993. doi: 111188001622159919999
H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), Advances in Neural Information Process- ing Systems, volume 30. Curran Associates, Inc., 2017. URL https:: meweeediiiinnnnnnii....e.. cc cc//ppper_iileee/ppper///////// 9 009d772288804b991647ccll833388 Papeee pdf.
Goodman, J. and Weare, J. Ensemble samplers with affine invariance. Communications in Applied Mathe- matics and Computational Science, 5(1))65-80, January doi: 10.21 40.cmmooccmmo 20010
2010. ISSN 2157-5452.
Kanagawa, H., Barp, A., Gretton, A., and Mackey, L. Con- trolling moments with kernel Stein discrepancies. arXiv preprint arXiv:221 TD0:22//.05400) 1.05408, 2022.
5.65. URL https: http///spporo. camcos/ /01
5-1/p04.. xxtml. Publisher: Mathematical Sciences Publishers.
Leucht, A. and Neumann, M. Dependent wild bootstrap for degenerate U- and V-statistics. Journal of Multivariate Analysis, 117:257-280 05 2013. 2013.03.003.
Gorham, J. and Mackey, L. Measuring sample qual- ity with Stein's method. In Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., and Garnett, R. (eds.), Advances in Neural Information Processing Sys- tems, volume 28. Curran Associates, Inc., 2015. URL
doi: 10.1016/j.jmm
Liu, Q. and Wang. D. Stein variational gradient descent moment matching. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa- N., and Garnett, R. (eds.), Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. URL cc cc/sapee_fiileeppppe////////// .ne.r.ereeddnnng............... neurips. 1 125533cc66577033edddd4eeeee pdf.
neurips.
cc htc/ppper_fiilleppppee////////// 99661aa198aal21cce11119911100000 pdf.
Paper.
Gorham, and Mackey, L. Measuring sample quality with kernels. In Precup, D. and Teh, Y. W. (eds.), Proceed- ings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learn- ing Research, pp. 1292-1301. PMLR, 06-11 Aug 2017. / URL https: http//poceedengge/ nggs ml htt.s//m.rr.rreee/0erm.1 gorham1 7a html.
Liu, Q. Lee, J., and Jordan, M. A kernelized Stein dis- crepancy for goodness-of-fe tests. In Balcan, M. and Weinberger, K. Q. (eds.), Proceedings of The 33rd Inter- national Conference of Proceedings of Machine Learning Research, pp. 276- 284, New York, New York, USA, 20-22 Jun 2016. PMLR. ngs ml r. www.psrs/see........... /
press/v70/
Machine Learning. volume 48
Huggins, J. discrepanciees/sc//mstttrr.
risd package, 2018.
feature--ent
URL /ntttps:/pproeee liubl6 6 html.
Huggins, J. and Mackey, L. Random feature Stein discrepancies. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Disa-Bianchi, N., and Garnett, R. (eds.), Advances in Neural Information Process- ing Systems, volume 31. Curran Associates, Inc., 2018. URL https::/proceedingss neurips. oo cccpppper_flllespppr////011 (rB400beet3cbbldcffod5ooooo pdf.
Mira, A., Solgi, R., and Imparato, D. Zero variance Markov chain Monte Carlo for Bayesian estimators. Statistics and Computing, 23(5):653-66 September 2013. ISSN 0960- 3174, 1573-1375. doi: 10.1007/s1 1222-222222244422122
Müller, A. Integral probability metrics and their generating classes of functions. Advances in Applied Probability, 29(2):429-443, 1997. ISSN 0001-8678. doi: 10.2307/ 1428011.
Hušková, M. and Janssen, Consistency of the gen- eralized bootstrap for degenerate u-statistics. The Annals of Statistics, 1114):811--,83, 1993. ISSN 00905364, 21688966. URL http:/ http://www. jstor. org/stable//24255
Nemeth, C. and Fearnhead, P. Stochastic gradient Markov chain Monte Carlo. Journal of the American Statistical Association, 116(533):433-550, January 2021. ISSN 0162-1459.
doi: 10.1080001621459.22222 210.559.22222111111111 1847120.
Pagani, F. Wiegand, M., and Nadarajah, S. An
Jitkrittum, W.
kerne-goof package, 2019. hutps://gthhb..co///eee howpa//bbb.mmmtttatttt//eeeeeoo.
dimensional Rosenbrock distribution for Markov chain Monte Carlo testing. Scandinavian Jour- nal of Statistics, 49(2):657-680, 2022.
Jitkrittum, W., Xu, W., Szabo, Z., Fukumizu, K., and Gretton, A. A linear-time kernel goodness-of-fit test. In Guyon, I., Luxburg, U. V., Bengio, S., Wallach,
ISSN
1467-9469. doi: 10.1111/$jos.1255 wi wiley com/
URL
10

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
doi/abs/10 http://001111///.....o
_eprint:
http://oniiiliiaaa.wwseeese//ooh 12532.
12532.
Rahimi, A. and Recht, B. Random features for large-scale kernel machines. In Platt, J., Koller, D., Singer, Y., and Roweis, S. (eds.), Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2007. URL https::/pprooee.nngg.o neurips.. cc cc/paper_faeee [oapperrfileeeppppe////////e
01 013a0066000bcc33@fffeefffffffff2 Paper. pdf.
Roberts, G. O. and Tweedie, R. L. Exponential convergence of Langevin distributions and their discrete approxima- tions. Bernoulli, 2(4):341-363, 1996. ISSN 1350-7265. doi: 10.2307/3318418.
Rosenbrock, H. H. An Automatic Method for Finding the Greatest or Least Value of Function. The Com- puter Journal, 3(3):175-184, January 1960. ISSN 0010- 4620. doi: 10.1099//comjn//33..5555 URL https: //doi.
3.175.
Serfling, R. J. Approximation theorems of mathematical statistics. John Wiley & Sons, 2009.
South, L. F., Karvonen, Nemeth, C.. Girolami, M., and Oates, C. Semi-exact control functionals from Sard's s method. Biometrika, 109(2):351-367, June 2022. ISSN 1464-3510. doi: 10. 10.109//biomeetsaab033..
South, L. F., Oates, C. J., Mira, A., and Drovandi, C. Regu- larized zero-variance control variates. Bayesian Analysis, 18(3):865 888, 2023, doi: 12.2211/22-BA2222 URL https: //t:i////i//////o/ org/o/og///oooo /1c 4 /22 2
Stein, C. A bound for the error in the normal approximation of dependent random vari-
to the distribution of ables. In Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 2:
Probubility Theory, volume 6.2. pp. 583-603. University of California Press, January 1972.
Welling, M. and Teh, Y. Bayesian learning via stochastic gradient Langevin dynamics. In Proceedings of the 28th International Conference on Machine Learning, pp. 681-
688, 2011.

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
A. Derivation of Closed Form Solution
Consider the optimisation problem
PSD sup SEG
[Q[Ag(()]]
sup
EQ
Σ (x))
|E 111211
km1
Bx [APb(T)]
sup BER1:)
k=
sup |
3x
sup BER:: k=1
Be
where Σ (x)) This can be written as an optimisation problem with Lagrange multiplier enforce the
constraint ||3||2 I (or equivalently || < 1). Specifically, we have PSD suppER L(3, A), where
L(B,A) Σ Bkēx 111111 - 1).
k=1
Computing the gradient of L(3, A) with respect to 8 and setting this to zero for the optimisation, we have
0 = V L A)
0 = x 2A3
2X'
Acknowledging that the supremum occurs when || = we have ||0||2 =
1 |2 21 1. so X | and thus
| Finally, substituting this solution for 3 back into the original objective, the solution for the supremum is
PSD =
sup BER D
BR t
SI
Z
Ex
||
2
k=1
B. Proof of Proposition 3.2
Proof. Using the notation x polynomial order r and
d
x we have g = spanfx" a E No. Ed ag S <r) for PSD with
PSD
Eq|Arap.
neNd 2 0.5r
Thus, PSD 0 if and only if Eg[Ax*] = 0 for all E No.
Σ rr. We will now proceed by proving that
Eg [Ax*] = 0 for all a E No. Σ 0: Sr (i.e. PSD 0) if and only if the moments of P and Q match up to order r.
12

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
The condition that A 0 for all a E No. E as can be written system of (d27) linear equations. For each a, we wish to understand the conditions under which
EQ [Ax = EqlAx* + V Vogg(()) Vx
Eq|Ax + (- (z- 14) V
Eq|Ax*1 BQ -ーー Vx BQ - V2
(12)
0.
This uses the property that P a Gaussian distribution as per the Bernstein-von Mises limit, so V logp(s) = -x-1 p.).
Since P decays faster than polynomially in the tails, we can also use the property that EE [Ax]]] = O. This leads system of (d=") linear equations, where for each
Ep[Ar*] Ep|Ar" + Vx (-2- =1)l
Ep|Ar"] Ep(Vx Mー Ep V== p1]
(13)
0.
Subtracting (13) from (12),
E r ( RQ- Vx x-1x] Eq-r(Ve" -
(14)
where [()) is used as a shorthand for Eo [x)) [p][())) Equation (14) combines the condition that Bo(Ax=]
=
with the property that Ep [Ax = 0 into a single system of equations. Our task to prove that this system of equations
holds if and only the moments of P and Q match up order r. We will so using proof by induction on the polynomial order r.
B.1. Base case (r 1)
For 1, 2 simplifies to x; for E {l,..., d]. is simpler this case to work directly with (12) than to work with (14). For each i {L...., d}, have
Eo[Ax; = 0
Eq Ar:]
Eq|E~
E x
Vx, + Eq Σ #
Vr, 0
Eq
Eq
E [4] 0 Eq
Er. x)] = ΣT. μ4,
where Σ,¹ denotes the ith row of E-1 Vectorising the above system of equations we have
Eoleol = Σ- Eq|æ]
This proves that PSD with r = 1 zero if and only the first order moments (means) of P and Q match.
B.2. Base case (r = 2)
For r = 2, rd simplifies to m, for i,j E f1,...,d] dd. Following (14), for each i, we have
Eq p Array (k--p(E-'a 'x Vayag) + P(=- Vaynyy 0.
The first term is EQ-p p|Ar,a,] = 0 because ATąI, = 21, does not depend on so the expectations under P and Q match. Similarly, the final term r|u Va,0,] disappears because first-order moments match. Thus
EQ- 10-p(E-1>= Vr,æ,
EQ p Er xw,l + Eg- ( xx,]] =
( BB-- - r (x.r 1) is 15-1RQ r [x: T) -0 (I- Eq p| /2 1) + (Eq- 10ppplan jax T ) = 0
( A)1, + (AE- = 0,
(15)
13

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
where A EQ-Plaa ( ]. Note that wxr gives all possible second order moments so if A then all second order moments match. We need show that A 0 is the only possible solution this equation.
Let Σー QDQT, where D is the diagonal matrix of eigenvalues with DEs A for &o E (1....,d) d). This is an dennppccoooooiiinnn so Q and QT are orthogonal (Q" = =Q-1). Now, (15) holds for all ž,j so the vectorised system of linear equations is
E -E-1 = QDQ A+ AQDQQ 0
DQ AQ + +QAAQD = 0
D3 + 3D = O,
where B Q AQ. The third line comes from multiplying by Q on the left-hand side and Q on the right-hand side. Consider now the ijth element:
(DB)y + (3D); (A+ 11) O.
We know that A for all {I,...,d] d) because - symmetric positive definite, so it must be that 0 and therefore A 0 is the only solution.
This proves that PSD with 2 if and only if the moments of P and Q match up to order two, i.e. the means and (co)variances match.
B.3. Inductive Step
Suppose that the moments up and including order match. Noting that Be-rlAue is of order 2 and EV- pp] is of order r - 1 and using (14), we have
Eq-r Σ 2 =0 0.
(16)
We must show this signifies that the moments of order must also match.
For general rth order monomial, xd can alternatively be written II Zim. for im (l,...,d). Therefore have
EQ-P Σ
Fim
Σ PE-PPLE
II
Fimml 1 = 0.
m=L.mn
This can be written in matrix form for all possible multi-indices of order r using Kronecker products:
terms
0 Eq-pP 1z ox 2r
[ E p|x (x E or (1
Eq- p(x 2 (x
2 x + Eq p|x ox Tx
....0 .... x-
terms
Writing A = EQ-PP X x the products AB and CD, this becomes
@ x
E ! and using
D) = AB @ CD provided that one can form
A- * -00
Once again using Σ-1 =QQQ_ we have
QDQ A+ A(QQQ 818...01) A( - "QDQ") 0.
14

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
terms
Multiplying by Q on the left and (Q 8 Q 8
on the right
0 = DD A(Q @ Q.. Q) + Q" A(QDQ 810 Q.Q
+ QT A(I @
@ 10010000000000011000000000000000
0 DQ" (Q...QQ..00.0000. (0000)000AQ000000000000
@
+ A(Q 8 000
0 = DQ A(Q 8 .Q...000+0)) 00000000000000000000000000000000 8000001000000000111111101100
0 = D3 + 3 0 10 01) + + D& 8 .... 3( Q Q Q Q) 00 .
+ +Q A(Q (Q)
D),
where Q Q 3 Q 8 Q). Now consider the ijth entries in these equations, where i E (1...., d) and
j [l,,...,d-11 The Kronecker product of diagonal matrices is diagonal, so the terms (D@I@- 10010-----111 -.1..-11.. (I @I@D@.. I) up to (1 @
D) all consist of diagonal matrices with elements of À on the diagonal. The explicit value of X for any given i and is not important, for now consider indices 2g {L...., df for I...., P which may not be unique? Then we have
(Az, + +2a +... )0s = O.
Since all > 0 by positive-definiteness of _-1, we have that Big = 0 for all i and j and therefore
B = 0
QTA(QQQ0 1 2 Q) 0
A = 0,
where the third line comes from multiplying by Q the left and (Q Q ..... .T...QQ) on the right. We have now shown that [][Ax*] = O for a E N6 Σ a S r if and only if the moments of P and Q match up to order r. Thus PSD if and only if the moments of P and Q match up to order r.
C. Additional Empirical Results
This section provides additional empirical results, as explained in Section 4.
C.1. Asymptotic Goodness-offfit Test
null rejection rate
0.20
rejection rat
1.0
Full rejection rat
1.0
Full rejection rat
1.0
RFSD
Gauss KSD
0.15
IMQ KSD
0.10
0.5
Gauss FSSD-opt
0.05
PSD r1
PSD r2
0.00
10
21
0.0
10
20
0.0
10
20
10
20
PSD r3
dimension
dimension D
dimension D
dimension
PSD r4
(a) Unit Gaussian (i.e. Q (b) Variance- Variance-perturbed (c) Student-t, 250 repeats
(d) Laplace, 500 repeats
P). 500 repeats
Gaussian, 200 repeats
Figure 5. Type error rate (a) and statistical power (b,c,d) for detecting discrepancies between the unit Gaussian P and the sampling distribution Q. using the test based on the asymptotic distribution of the U-statistic.
We reproduce Figure 1 in Section 4.1 of the main paper by using Corollary 3.1 to implement a test based on the asymptotic distribution of the U-statistic, PSD Specifically, we set P = N NOa,, Ia) and assess the performance for a variety of Q
2
2To he explicit in an example of fourth order polynomial, consider instead and j indices starting zero so i E (0.....d-1) and j E (O....,d d 13. Then the explicit form of the equation is (A₃ + Aassaa ALsajj%a + A5% where a%d = dla/d] the remainder operator.
15

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
and d using statistical tests with significance level a 0.05. We draw 5000 samples from the null distribution given in Corollary 3.1 using the plugin estimator E.g. We investigate the four cases described in Section 4.1 with identical settings in the main paper.
The results are shown in Figure 5. Primarily, while PSD displays good performance for the variance variance-pertuubed (r = 2 and r = 4) and Laplace (r = 4) experiments, PSD with r = 4 does not maintain high power in the case of the student-t. Also, we note that the type error is not well controlled. Further, the asymptotic test can be slower than the bootstrap lest due to the computation of the eigenvalues of the covariance matrix (complexity (JJ3)). This can be improved by considering approximations the covariance matrix, however this would introduce bias asymptotically.
We recommend the bootstrap procedure for goodness-of-fit fit testing over the asymptotic test due to its higher statistical power, better control over type I error and lower computational complexity.
C.2. Detecting Second Moment Discrepancies
Section 4.1 of Kanagawa al. (2022) illustrates the failure of the standard KSD with IMQ kernel for detecting non- convergence of the second moment of simple sequence, (Qu)n21 to its target P, set to be the standard Gaussian of d dimensions, P = (N,O The sequence On is defined as follows:
On = (1-
Pr +
6Rn:
n1 +
n
where Pa = Σ ox,- and In -n++++ 1.. 1 with (1,..., 1) (a vector of ones), and {X1, .,,,,,,,,)),)).. iid Pare i.d. with X; ~ P.
Kanagawa et al. (2022) demonstrate that the sequence Qn converges to P almost surely. However, it has the following (almost sure) biased limit:
lim By (Y & YJ = Ex-p X & x]+ 1181. n-000
Figure shows trace plots of KSD, PSD and the Euclidean distance between the true and estimated moments up to order While KSD appears to decay exponentially to zero, the trace plot for our recommended PSD with r 2 is remarkably similar to the trace plot for the Euclidean discrepancy in moments up to order 2. This is in contrast to the method from Kanagawa al. (2022), which quickly asymptotes statistical power, we believe PSD has practical advantages since it is aimed at directly tracking discrepancies in moments. Moreover, we achieve moment tracking in linear time rather than the quadratic time of Kanagawa et al. (2022).
non-zero value. While their approach is likely to have higher
16

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
2.0
IMQ KSD
0.5
0.2
200
Sequence index
PSD (1=2)
20
50
200
500
1000
Sequence index
Sq Euclddan Distance in Moments (r=2)
10.0
0000 Pooo
2.0
0.5
2
20
50
200
500
1000
Sequence index
Figure 6. Moment tracking capahilities of Example Kanagawa (2022)
C.3. Logistic Regression: Sparsity Prior
We extend the empirical investigations to consider a more challenging logistic regression task with a sparsity-inducing prior. We use the stochastic search variable selection (SSVS) prior of George McCulloch (1993), corresponding to mixture of Gaussians with standard deviations 0. and 10. The resulting posterior is multimodal with modes capturing the zero and non-zero estimates from the sparse prior. We simulate data with 100 observations and 20 variables, of which three non-zero. We multiple ULA runs with different step sizes and different initialisations to thoroughly investigate the performance of PSD.
Our results (averaged over runs) are presented in Table 2 and the posterior distributions corresponding to the different step sizes are given in figure 7. The normalised median discrepancy values are reported. KSD and PSD) with r have similar performance and fail to identify the poor performance of ULA for large step sizes. PSD with 2 identifies the step size with the best moment estimation. In particular, despite posterior being far from normal, our methods still attain superior performance for assessing moment discrepancies compared to IMQ KSD.
17

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
Step Size
KSD
PSD1
Eucl1
PSD2
Eucl2
PSD3
Euc13
PSD4
Eucl4
1 x 10
1.0e+00
1.0c+00
6.2c-02
1.0c+00
2.2c-02
1.0e+00
7.4c-03
2.3c-01
2.4c-03
5 X 10-5
5.7e-01
5.6e-01
5.8c-02
5.6c-01
2.1e-02
5.7e-01
7.2c-03
1.3e-01
2.3e-03
1 x 10
4.3e-01
4.3e-01
5.4e-02
4.5e-01
2.0e-02
4.8e-01
6.7e-03
1.2e-01
2.2e-03
5 x 10-4
2.0e-01
2.0e-01
3.9e-02
2.3e-01
1.Se-02
2.6e-01
5.0e-03
6.6e-02
1.6e-03
1 10-3
1.4c-01
1.4e-01
3.1c-02
1.6c-01
1.2c-02
1.9c-01
3.8c-03
5.1e-02
1.2c-03
x 10
7.2e-02
6.7e-02
1.56e-02
8.38e-02
4.88e-03
1.07e-01
1.79e-03
3.25e-02
5.75e-04
1 1
6.6e-02
4.7c-02
7.0e-02
1.0c-01
2.9c-02
1.5e-01
1.1e-02
5.4c-02
3.8c-03
5 x 10
2.9e-02
2.4e-02
7.6e-01
1.5e-01
6.5e-01
5.3e-01
5.4e-01
4.5e-01
4.5e-01
1x 10
2.78e-02
1.70c-02
1.0c+00
2.0c-01
1.0c+00
9.7e-01
1.0e+00
1.0e+00
1.0c+00
Table Normalised discrepancies for the logistic regression example with sparsity-inducing prior. Shown are the KSD, the PSD and the Euclidean distance between estimated and gold-standard moments of order most The lowest value shown in bold.
Step Size 0.05
Step Size 0.005
Step Size 5e-05
ULA Semple Tree Pastoroo
ULA Sample Tue Pastoro
UBA Sangle Tree Postrir
0.1
01
0.1
besa2
beta2
beta2
-0.1
-0.1
-011
-0.3
-0.3
-0.3
-02
0.0
0.1
02
0.
04
05
-02
00
01
0.
04
0.
-02
0.0
02
03
0.4
0.5
beta1
beta1
beta1
Figure ULA sample plots in blue approximating the true posterior (red) for logistic regression example with sparsity-inducing prior. PSD selects step size of 0.005 which visually provides good approximation the true posterior.
C.4. Logistic Regression: Big Data Example
We consider another logistic regression example looking the performance of PSD when the posterior close to normal. This example is logistic regression with 10t observations and d 5 variables using Gaussian priors with standard deviation of 10, thereby making the Bernstein-von Mises approximation suitable. Similarly to the previous example, we run each example five times with multiple ULA initialisations. The results are shown in Table and the posterior distributions corresponding to the different step sizes given in figure 8. We find that PSD performs well in identifying the best choice of step size in ULA, though the performance for PSD with r = 1 is not optimal. Looking at plots of the marginal gradients versus the marginal samples as heuristic, as per the suggestion made in Section 5, we find that when the samples are largely in the tails (e.g. for larger step-sizes), the Gaussian approximation is poor. Higher-order PSD performs excellently and in particular, PSD with = 2, which is what recommend in general, is both computationally cheap and capable of dealing with samples taken mainly in the tails.
18

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
Step Size
PSD1
Eucl1
PSD2
Euc12
PSD3
Eucl3
PSD4
Eucl4
1 x 10
7.44c-01
5.92e-05
1.97e-04
7.98e-07
1.80e-06
7.70c-09
1.1le-08
7.12e-11 11
5 x 10-5
2.80e-01
2.89e-05
7.07e-05
4.22c-07
6.75e-07
4.30c-09
4.26e-09
3.70e-11
1 x 10
2.02e-01
2.31e-05
5.14e-05
3.41e-07
4.92e-07
3.68e-09
2.76e-09
3.35e-11
5 10-1
8.19e-02
1.40e-05
2.30e-05
2.51e-07
2.35c-07
2.97c-09
1.61e-09
2.77e-11
1 10-3
6.25e-02
2.10e-05
3.20e-05
3.89e-07
3.74c-07
4.65e-09
2.71e-09
4.35c-11
x 10-3
5.89e-01
2.92e-02
5.89e-02
1.19e-03
2.21e-03
3.53e-05
5.52e-05
9.13e-07
1 10
7.32e-01
7.14e-02
1.19e-01
5.46e-03
7.98c-03
3,40e-04
4.71c-04
1.97e-05
5 x 10
9.23e-01
4.60e-01
5.03e-01
2.08e-01
2.10e-01
9.19e-02
9.63e-02
4.05e-02
1 10-1 1.00e+00
1.00e+00
1.00c+00
1.00e+00
1.00e+00
1.00e+00
1.00e+00
1.00c+00
Table Normalised discrepancies for the logistic regression example with big data. Shown are the PSD and the Euclidean distance between estimated and gold-standard moments of order most The lowest value is shown in bold.
Step Size 0.001
Step Size 5e-04
Step Size 1e-05
0.10
0.10
0.10
0,00
0.00
0,00
bebb2
beta2
beta2
-0.10
40100
-0100
-0.20
ULA Contbur
ULA Contour
-0.00
USA Gantoul
-020
-0 10
0.00
0.05
0.10
-020
10
0.00 0.05
0.10
-020
-0.10
0.00
0.05
0.10
beta1
beta1
beta1
Figure ULA sample plots in approximating the true posterior (red) for the logistic regression example with big data. selects step size of 5 10 which visually provides good approximation the true posterior.
C.5. PSD without interaction terms
In Section we proposed an alternative PSD methods that excludes interactions achieve computational complexity that scales linearly with dimension, specifically ((drr). In the context of Gaussian P with independent components, this means the inability to detect discrepancies in correlation.
We now empirically investigate the performance of PSD without interactions on the SGLD example from Section 4.2 and on the two additional logistic regression examples. The results are shown in Figure for the SGLD example and in Figure 10 for the logistic regression examples. The performance remarkably similar to the performance with interactions. This also true for the logistic regression examples when the covariates are simulated using an auto-regressive process with high positive negative autocorrelation.
We do not investigate this for the Gaussian goodness-of-fit examples since we know there is no disadvantage to removing interactions in that context, and the statistical power would therefore be similar to or better than the power with interactions.
We believe PSD) without interactions should generally perform similarly to PSD with interactions when applied to biased sampling methods like ULA and SGLD. Removing interactions may have more effect on bespoke sampling methods that treat multivariate relationships differently.
19

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
10
RFSD
discrepancy measure
IMQ KSD
PSD r1
PSD r2
10"
PSD r3
PSD
PSD no
PSD no
PSD no
10
10
10
10
E
Figure 9. PSD with and without interactions for the SGLD example from Section 4.2.
1.0
1.0
Scaled median discrepancy
PIDT
Scaled median discrepancy
PID
PIZUS HOSESR
PSOS notrarad
0.8
P rotnerra PSDS rnsideeat
0.8
PGCC3 POreen
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0
-5
4
-3
-2
-1
-5
-3
-2
-1
log step size
log step size
Figure 10. PSD with and without interactions for the logistic regression examples with (a) pparsity-inducing prior and (b) big data.
C.6. Rosenbrock Target
We conduct an empirical study with the target P, set to the two-dimensional Rosenbrock function (Rosenbrock, 1960; Goodman & Weare, 2010), given by:
p(x) exp (-(11- -0))-11)2--0 b - ==²=²))) =²))),
where b a constant and [p[r]] = fr- Inference for this target is challenging because there is high probability mass in narrow, curved region with complex dependencies and different scales across dimensions (Pagani al., 2022).
We use this non-Gaussian example show analytical reasoning and empirical investigations for case where we failing to track differences in the first moments. We provide the form of the expectations we are tracking instead when = 1.
By explicitly writing the form of the PSD, we can see that PSD with r = zero and only if 2b [EQ[21122] Eg (ai))
Eq[z1] | - 141 (associated with monomial 21), and Eqq [= = Eq [z2] (associated with monomial 22)). This can be derived by setting [Eq[V [.,, log p(x)] = 0 (associated with monomial æ ), and Eo s., log p(x)] (associated with monomial xr)) and rearranging, Unlike in PSD for Gaussian targets, this is not simply assessing the accuracy of the first moment. The equation associated with monomial Ii may hold when E2 [x:]] is correct (i.e. [q]=]] f1), ). but it could also hold when [ is incorrect and this is offset by having Eo [r) * [=[=]]] Similarly, the equation associated with monomial 22 is not directly assessing whether the moments of P and Q match, but rather whether certain relationships between moments are correct. This is different from the situation when we have Gaussian P.
The conditions under which PSD with higher r is zero can also derived. For example, to have zero PSD with r = would also require that Eg[2 (2 + 201 Vz, log p(x)) (monomial æ²)) R2 + 2r V.y, logp(x)) = (monomial x??) and
20

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
Bo Vaz log p())()) + 22 Va, log logp(x)] (monomial x1 However, these become increasingly complex and difficult interpret for higher> r.
To empirically investigate the performance of PSD, we provide time series plots of PSD and KSD with increasing using samples either from P (Figure 11) or from an altered Rosenbrock target with incorrect (Figure 12). For the former, all expectations converge to their true values as n - oo and we expect all discrepancies to go towards zero. This what we find empirically as well. For the latter, the first moment is correct but the variance incorrect. One might then hope that PSD with empirically in Figure 12. Importantl, since P is far from Gaussian, PSD well. Nevertheless, it may be useful discrepancy, even for this unusual P, for certain tasks like choosing sampler.
would convergence towards zero but this the case, as we have discovered analytically and
tracking the discrepancies in the moments
1e+02
5.00
Discrepancy
Discrepancy
0.50
1e+00
0.05
KSD
KSD
PSD1
1e-02
0,11
Euclidean
Euclidean
100
1000
10000
10
100
1000
10000
No. of Samples
No. of Samples
Figure 11. Time series plots of KSD, PSD and Euclidean distances between estimated and true moments up to order r when sampling from the correct Rosenbrock target. Results are shown for (a) r and
1e+02
5.00
Discrepancy
Discrepancy
0.50
1e+00
0,05
KSD
KSD
PSD1
PSD2
0.01
Euclidean
16-02
Euclidean
10
100
1000
10000
10
100
1000
10000
No. of Samples
No. of Samples
Figure 12. Time series plots of KSD, PSD and Euclidean distances between estimated and true moments up order when sampling from Rosenbrock distribution with misspecified Results are shown for (a) r and (b) 2.
21

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
C.7. Translation of Mean
The PSD is not translation invariant. Specifically, one could inflate the value PSD by considering arbitrary translations of samples by their mean. To further investigate, we have performed empirical investigations into the effect of the transformation i == MQ. where μQ is the mean of Q. We believe this is the most sensible and practical translation consider. However, using a mean-shift, described, does not affect which discrepancies PSD is theoretically capable of detecting but it may affect its statistical power in doing so. It can be observed that PSD with solely on the score function, which does not change with such a transformation. However, PSD with higher r uses both the samples and the score function in the discrepancy so it can be sensitive to the mean-shift.
is mean-shift invariant because it is based
We consider cases based on Gaussian P with unit covariance, N 100 and d 5. Further, are interested in the case where the mean of P (ju)) is not the same across all dimensions, since believe this is where the impact statistical power will have the most effect. For this reason, we consider two cases, (1) μp (a, ......... 0) and (2) μp (0, a, ................................ a). The discrepancy will be in the first dimension in both cases, so we are interested in how the relative scales of the means that are correctly specified (versus misspecified) affect the results. We consider the two cases mean translation for situations with misspecified mean as well as misspecified variance.
Table We consider both case (1) and case (2) in situations where the discrepancy is in the (first) mean (uq pp 0.5ee) or in the (first) variance (Eoql1 [p]] + 0.5). A "t" in front of the discrepancy indicates we have performed mean-shift reparametcrisation.
a
Case Misspecified
Discrepancy
-1.5
-1.0
-0.5
0.0
0.5
Mean
PSD₂
0.70
0.11
0.06
0.16
0.87
tPSD2
0.98
0.22
0.06
0.13
0.84
PSD3
0.71
0.26
0.06
0.18
0.76
tPSD3
0.92
0.30
0.05
0.19
0.74
2
Mean
PSD2
1.00
1.00
0.76
0.20
0.08
tPSD2
1.00
1.00
0.84
0.14
0.06
PSD3
1.00
0.96
0.72
0.26
0.06
(PSD3
1,00
1.00
0.82
0.25
0.10
1
Variance
PSD₂
1.00
0.96
0.78
0.55
0.32
tPSD2
1.00
0.98
0.84
0.54
0.32
PSD3
0.96
0.92
0.70
0.36
0.17
(PSD3
0.97
0.90
0.70
0.32
0.17
2
Variance
PSD2
0.97
0.83
0.62
0.39
0.27
tPSD)
1.00
0.98
0.88
0.56
0.34
PSD3
0.83
0.65
0.44
0.29
0.16
tPSDa
0.98
0.90
0.70
0.36
0.16
Table 4 shows the estimated statistical power based 200 independent simulations. The results demonstrate that the original mean scaling does affect the results, but the performance with the mean-shift reparamerisation is generally similar to, if not slightly better than, the performance with no reparameterisation. Importantly, the choice of parameterisation is a problem that affects Stein discrepancies more generally. KSD with radial kernels, i.c. kernels that are functions of the form || ylll|. like the Gaussian kernel, are mean-shift invariant. However, they are still sensitive to other types of reparameterization, such as whitening (considered in Appendix D).
D. Using an Invertible Linear Transform
Corollary D.1. Consider PSD applied on the transformed space y = WX. where W E Bdxd is an invertible matrix that is independent of { Denote the distribution of Wr by Q. Using a change of variables, V log p(y)
W-T V, p()). The new discrepancy,
PSDw
Σ
EglAny + V, logí(y) V y
GENg Rd
zero if and only if the moments of and Q match up to order r.
22

---

The Polynomial Stein Discrepancy for Assessing Moment Convergence
Proof. Under the Bernstein-von-Miies limit (or generally for Gaussian targets), P = (u, E) so in the transformed space P= Σ) where jı W and Σ WWWT.
Next, we will show that Σ is symmetric positive definite. Since
WEWT ΣT, Σ is symmetric. Since Σ is positive-definite = Σ> = 2 WEW > 0 only requires that we do not have W 2 = 0 and therefore by invertibility of W. that we do not have 2 = 0. Thus, z T E2 = 2 WEW > for all non-zero 2 E R which by definition means that positive-definte.
The assumptions of Proposition 3.2 (symmetric, positive-definite covariance) are met for this transformed space so for any E N we have
(r terms))
(r terms)
V O @ .... 15 = y B Gy 2 1y.
E E W
(W
8 (Wx) ExeP (W)( (W..)
WE Ex 8... @ x
W = WE W @ 8... 0z [( ....8WT)
= Ez- 0 0x
Therefore we have that PSD applied to the transformed match up to order r.
WW is equal to zero if and only if the moments of / and Q
23


Reviews:
[{"review_id": "YLWcWlfwQl", "reviewer": "Reviewer_iLVd", "summary": "The authors introduce a variant of Stein discrepancy which uses bounded degree polynomials as a Stein set. Computing the proposed Polynomial Stein Discrepancy (PSD) is straightforward using evaluations of the target score $\\nabla \\log P(x)$ and samples from the proposal distribution $Q$. For degree $r$, evaluating PSD given $n$ samples takes time $O(n {d+r \\choose d}) \\sim O(n d^r)$, whereas Kernel Stein Discrepancy requires $O(n^2 d)$ steps, which can be prohibitive for large sample sizes. \n\nTo validate PSD, the authors compare it to KSD with IMQ and Gaussian kernels, as well as two competing linear time discrepancy metrics (RFSD and Gauss FSSD-opt). For goodness of fit testing, the PSD demonstrates exceptional performance relative to competing methods in distinguishing $P=\\mathcal{N}(0, I_d)$ from a few different proof-of-concept $Q$ distributions, even outperforming KSD. The PSD is also competitive with KSD and outperforms other linear time methods for goodness-of-fit testing for the samples generated by an RBM. Finally, for parameter tuning in a simple test-case LMC method, the optimal parameter choice under PSD agrees with that of KSD, while being significantly cheaper to compute.\n\n## update after rebuttal \n\nI appreciate the authors engagement with my questions and I keep my score.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your time reviewing our manuscript and for your valuable comments and suggestions.\n\nWe will correct the typographical issues and the appendix as suggested.\n\n**When does KSD outperform PSD? In Figure 1, there are no such examples, and in Table 1 setting $r>1$ is sufficient to match KSD. What are some practical examples of non-moment-based perturbations that KSD can detect but not PSD?**\n\nProbability distributions are completely determined by their moment generating functions, provided it exists (a necessary condition is that all moments are finite). Therefore, we can potentially expect the KSD to outperform PSD for target distributions lacking well-defined moments (e.g. those with heavy tails). We comment on this further in response to Reviewer 8qrD. We can also expect KSD to outperform PSD of order $r$ when the discrepancies lie in moments higher than the $r$th order moment. For example, KSD outperforms PSD with $r=1,2,3$ for the student-t and Laplace examples in Figure 1 since the discrepancy is in the kurtosis ($4$th order moment). We will add details of when we can expect KSD to outperform PSD to the discussion.", "claims_and_evidence": "The experimental claims made in this work are supported by established benchmarks. The derivations and proofs in this work are correct to the best of my knowledge.", "methods_and_evaluation_criteria": "The experimental claims in this work are well supported. The authors compare PSD to other algorithms using established benchmarks, such as the goodness-of-fit tests used in Gorham & Mackey 2018, the RBM sampling task used in Liu et al. 2016 and Jitkrittum et al. 2017, and the bimodal Gaussian task introduced in Welling & Teh (2011).", "theoretical_claims": "I checked the derivation of PSD, the proof of Proposition 3.2, and the proofs detailed in Appendices A and B.", "experimental_designs_or_analyses": "I am unfamiliar with the literature on discrepancy metrics for testing sample quality. However, to the best of my knowledge, this work has used well-established experimental baselines to demonstrate the efficacy of the proposed approach. To the best of my knowledge the comparisons made between PSD and competitors are fair.", "relation_to_broader_scientific_literature": "This paper belongs to the literature on Stein discrepancies and goodness-of-fit testing. The proposed method is a conceptually simple instantiation of this idea which is shown to outperform other linear-time methods in the literature.", "essential_references_not_discussed": "N/a", "other_strengths_and_weaknesses": "One strength of this work is that its writing is extremely clear and easy to understand. Also, the authors give many practical insights, e.g. \"we recommend the bootstrap test in general because we empirically find that it has higher power than the asymptotic test using samples from $Q$\" \"we have also found empirically that the power of the goodness-of-fit test based on RFSD is reduced when direct sampling from $P$ is infeasible,\" etc.", "other_comments_or_suggestions": "In Appendix A, only the first sequence of equalities is necessary. The fact that $\\\\|\\\\bar{z}\\\\|\\_2^2 = \\\\sup\\_{\\\\beta \\\\in \\\\mathbb{R}^J : \\\\|\\\\beta\\\\|\\_2 \\\\leq 1} \\\\sum\\_{k=1}^J \\\\beta\\_k \\\\bar{z}\\_k$ is a well known property called 'self duality of $l^2$ norm' and it does not need to be re-proven. \n\nSome very minor typos: \n- Line 104 left 'PSD discrepancy' should just be 'PSD'\n- Line 95 right should read $\\mathbb{E}_{X \\sim P} [(\\mathcal{A}_x^{(2)} g)(X)]$, inside the expectation is a random variable", "questions_for_authors": "When does KSD outperform PSD? In Figure 1, there are no such examples, and in Table 1 setting $r > 1$ is sufficient to match KSD. What are some practical examples of non-moment-based perturbations that KSD can detect but not PSD?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "pSLu2vrEyt", "reviewer": "Reviewer_7ab6", "summary": "The paper proposes a Stein discrepancy metric (PSD) for hypothesis testing. Rather than attempt to target a broad set of functions, the paper argues that by restricting to polynomial moments (up to some order), one can obtain a metric that is light-weight for computation and also has better statistical power.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your time reviewing our manuscript and for your valuable comments and suggestions. We will correct the typographical error.\n\n**What if I wanted to test against a specific set of functions other than those given in the paper? Is there any other example of this generating an interesting set of statistics?**\n\nIt is possible to use functions other than polynomials, but the challenge with doing this is that the application of the Langevin Stein operator, which is required for the zero expectation property, modifies the form of the function. This means we may not obtain interpretable forms for other functions like we do for polynomial functions and Gaussian distributions. As an example, if one were to consider $g(x) = \\sin(x)$, then after the application of a second-order Langevin Stein operator for a one-dimensional unit Gaussian target, the function would become $\\mathcal{A}_x^{(2)} g(x) = -\\sin(x) - x \\cos(x)$. A discrepancy based on $g(x) = \\sin(x)$ would therefore determine whether $\\mathbb{E}_Q[\\sin(x) + x \\cos(x)] = \\mathbb{E}_P[\\sin(x) + x \\cos(x)]$, which is not the same as the desired function $\\sin(x)$.\n\n**For instance, if I wanted to test against the moments after applying some transform to my random variable, could I possibly concoct interesting tests this way?**\n\nThis is an interesting question. Appendix D shows that applying an invertible linear transformation does not change the fact that, in the Gaussian case, the test is for the moments up to order $r$ of the original distribution. However, the statistical power can be affected by the choice of transformation. In particular, such a transformation may improve power when the posterior variances differ substantially in scale. It would be interesting to consider other transformations as well, and we will mention this in the discussion.\n\nFollowing up on this and a comment by Reviewer JxPo on the effect of a mean shift on the PSD statistic, we aim to empirically investigate the effect of a mean shift during the review period. \n\n**My main concern is that the paper could have benefited from further exploration of some of the ideas contained therein.**\n\nWe have further explored and commented on several topics as part of this review. We would be happy to consider further exploring other aspects of the paper if you would like to provide further details.", "claims_and_evidence": "The authors claim some empirical benefits and provide appropriate benchmarks. They also provide a method for approximating their discrepancy metric from samples, and compute some asymptotics to characterize the quality of the approximation.", "methods_and_evaluation_criteria": "Yes, the benchmarks are thorough and the comparison between methods seems fair.", "theoretical_claims": "As noted in claims, the primary theoretical claim relates to the asymptotics of their approximate statistics. There are also detailed derivations to justify correctness of discrepancy metric.", "experimental_designs_or_analyses": "The paper provides some assessments of PSD against other common KSDs, for common/fairly standard data sets.", "relation_to_broader_scientific_literature": "The paper relates to known literature on the construction/selection of KSDs. More specifically, it discusses common KSDs in the context of their computational tractability. The relative merits and downsides of various prior schemes are mentioned, in order to highlight the need for KSDs which are easily computed and yet which still have statistical advantages.", "essential_references_not_discussed": "As far as I am aware, the most relevant references have been covered.", "other_strengths_and_weaknesses": "I find the main idea behind this work to be rather compelling; namely, I like this concept of testing against a smaller subset of functions. I believe this could be interesting if developed further, for instance by considering other classes of test functions.\n\nMy main concern is that the paper could have benefited from further exploration of some of the ideas contained therein.\n\nNonetheless, the statistic is simple, computationally cheap, and I believe would certainly be of use to practitioners. Thus, I would recommend that this paper be accepted as it stands.", "other_comments_or_suggestions": "Line 322: sfix -> fix", "questions_for_authors": "What if I wanted to test against a specific set of functions other than those given in the paper? Is there any other example of this generating an interesting set of statistics? For instance, if I wanted to test against the moments after applying some transform to my random variable, could I possibly concoct interesting tests this way?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ENUJ5nWFzJ", "reviewer": "Reviewer_8qrD", "summary": "This paper proposes a class of monomials as test functions in kernelised Stein discrepancies, in order to speed up computations. It shows that when the target is Gaussian then the method works well.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your time reviewing our manuscript and for your valuable comments and suggestions.\n\n**On the BvM limit**\n\nProposition 3.2 explicitly states the Gaussianity assumption without mentioning the BvM limit, but we agree that this can be clarified in other parts of the text. We will reword the text in the abstract and the introduction of the paper. We will replace the text \"we prove that it detects differences in the first $r$ moments in the Bernstein-von Mises limit\" in the abstract with \"we prove that it detects differences in the first $r$ moments for Gaussian targets.\"\n\nWe opt to keep the text about the BvM as a remark after our theoretical result. This is relevant to the practical application of the method because a major application of KSD and related methods is in fast but biased Monte Carlo methods, such as SG-MCMC, where posterior distributions are often close to Gaussian. Bardenet et al (2017) highlight this in their review paper on MCMC for tall data, stating ``we emphasize that we have only been able so far to propose subsampling-based methods which display good performance in scenarios where the Bernstein-von Mises approximation of the target posterior distribution is excellent. It remains an open challenge to develop such methods in scenarios where the Bernstein-von Mises approximation is poor.\"\n\n**Method is based on comparing high moments**\n\nAs mentioned in Section 3.3, lower-order moments are often of interest and they can be where discrepancies are most likely to appear. We give the example of SG-MCMC, where large step-sizes lead to over-estimation of the posterior variance and small step-sizes can lead to under-estimation of the posterior variance. We therefore do not necessarily agree that the method is based on comparing high moments.\n\n**Rosenbrock target context**\n\nWe will explain that it is challenging because there is high probability mass in a narrow, curved region with complex dependencies and different scales across dimensions. We will add references to the original work (Rosenbrock, 1960), the probabilistic version (Goodman and Weare, 2010) and a paper talking about its complexity (Pagani, 2022).\n\n**Gaunt (2020) and relations to theoretical underpinnings**\n\nThank you for the reference. Since the polynomial functions are not Lipschitz, we cannot directly compute the solution of Stein's equation and establish bounds on the derivative. The suggested paper seems to describe a technique to establish the required bounds for unbounded test functions which seems suitable for the PSD case. However, it is not clear how to establish this for our formulation of PSD.  Further, the main goal of the paper seems to focus on establishing the closeness of the distribution of $g(W)$ to the distribution of $g(Z)$, where $W$ is a sequence of random variables approximating a standard Gaussian. We are unsure of the relevance of this continuous mapping theorem to our problem - what would be a possible interpretation of the function $g$?\n\nInitial results on the convergence of PSD for a simple setting are given in response to Reviewer JxPo.\n\n**Relationship to Xu et al (2022)**\n\nThe resampling method is indeed similar to the method we have considered. Following Liu et al (2016), we are able to establish the asymptotic exactness of our bootstrap procedure. Theoretical guarantees for resampling under a normal approximation could be a valuable reference to extend the method. Upon our preliminary review of that paper, we feel that the results in Appendix $B.1$ are of interest and we will cite this as a potential extension for the bootstrap procedure. Please let us know if you were referencing a different theorem or section of the paper.\n\n**What if moments do not exist?**\n\nWe believe this could affect things in two ways. First, we might not have the critical zero expectation property for the Stein operator applied to the polynomial, hence we might not expect that PSD $= 0$ when $P = Q$. Second, we will not be assessing moments exactly because the gradients of the log target will not be linear in $x$. We will mention limitations around moments not existing in the discussion.\n\n**Dependence on Langevin Stein operator**\n\nThe advantage of using the Langevin Stein operator is that when it is applied to monomials for a Gaussian target then we obtain monomials back, and therefore we are tracking convergence in moments. This may not hold true for other Stein operators. Nevertheless, applying diffusion Stein operators and other Stein operators, such as those considered in Kanagawa et al (2022), to different types of polynomials is an interesting question and can be studied further in future work. We will mention this in the discussion.\n\n**Anonymous Github site**\n\nThe full set of code has been provided as a zipped folder in the supplementary material for (optional) review.  We would prefer to use the first author's Github account if the paper is accepted to encourage appropriate attribution.", "claims_and_evidence": "The claims are \n\n* the method detects differences in the first r moments in the Bernstein-von-Mises limit;\n* the corresponding KSD test has higher power than some of its competitors, in some empirical examples.\n\nIt is the first claim that I find misleading. The paper shows that when the target is Gaussian, the method detects differences in the first r moments. It somewhat artificially puts the result in a Bayesian framework by instead of clearly stating that the method works well on Gaussian target, the Gaussian is called the `Bernstein-von-Mises limit'.", "methods_and_evaluation_criteria": "The empirical evaluations look ok, although the code is not available (for blinded review, it states, but it is possible to put the code on an anonymous github site).", "theoretical_claims": "The proofs are correct from what I see. \n\nHowever the method is based on comparing high moments. If the underlying distribution is, say, Cauchy, and not Gaussian, then it is not clear that the method would work at all.\n\nWhile the limiting behaviour of posterior mean (or mode or median) often follows a Gaussian distribution, this is what I believe the authors call the Bernstein-von-Mises limit. This convergence only holds under some conditions; see for example papers by Spokoiny et al. and Kasprzak et al.  which show that the rate of convergence depends on the moments. Hence when the moments in the test are high, the posterior mean may still be far from normal and hence the asympotic regime is not warranted.", "experimental_designs_or_analyses": "The experimental setup in the main paper is standard and useful to compare against other distributions. In the supplementary material also a Rosenbrock target is used. I found it difficult to get more information about this distribution; it seems to have an unknown normalising constant and a narrow ridge around the mode. This target distribution is clearly not Gaussian and the simulations show that the proposed method does not work very well in this case.", "relation_to_broader_scientific_literature": "In the area of Stein's method there are results available for multivariate normal approximation using as test functions those of polynomial growth, \n\nGaunt, Robert E. \"Stein’s method for functions of multivariate normal random variables.\" (2020): 1484-1513.\n\nIt would be good to relate the approach in the paper to theoretical underpinnings; in particular bounds on convergence rates could be of interest. \n\nIt would also be good to relate the results to quantitative Bernstein-von-Mises theorems such as the ones mentioned above.\n\nThe resampling method is very related to \n\nXu, Wenkai, and Gesine D. Reinert. \"A kernelised Stein statistic for assessing implicit generative models.\" Advances in Neural Information Processing Systems 35 (2022): 7277-7289.\n\nThat paper includes normal approximation as a special case and gives explicit theoretical guarantees for re-sampling.", "essential_references_not_discussed": "None beyond the above mentioned papers.", "other_strengths_and_weaknesses": "The idea of using polynomials as test functions is new in the context of KSD. However in my view the paper is hiding the key results under the notion of Bernstein-von-Mises limit, when really it means that the target is Gaussian. This could be seen as misleading.", "other_comments_or_suggestions": "Whether or not the rth moments exist is not discussed in he paper.", "questions_for_authors": "* What happens if the target does not have any moments? \n\n* How much do the results depend on the Langevin Stein operator? What happens if a different Stein operator is used?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "QRr1XoAw5h", "reviewer": "Reviewer_JxPo", "summary": "In Bayesian statistics, it is quite common to want to integrate some\nfunctions with respect to the posterior distribution (e.g. the mean). When\nthe posterior is complicated, this has no closed form solution, and so\npractitioners often resort to using MCMC samplers or diffusion samplers. In\nthe past decade, the Stein discrepancy was developed as a way to quantify\nthe quality of a sample to a target distribution, even if the\nnormalization constant of the target distribution was unknown. The typical\nchoice of the Stein set is usually an RKHS associated with the Gaussian or\ninverse multiquadric kernel.\n\nThe authors propose using a different Stein set here, i.e., a Stein set\nwhere the preimage is the vector space of all monomials with degree less\nthan or equal to $r$. The authors argue the benefits are twofold: this can\nbe more computationally efficient when $r$ and $d$ are relatively small\ncompared to $n$ and in the case of the target is Gaussian, it focuses the\ndiscrepancy more on the similarity of the moments of the two\ndistributions. They offer some theory showing that this holds in the of a\nGaussian and they also illustrate the efficacy of this discrepancy on a many\ndifferent examples commonly used in this literature.\n\n## update after rebuttal\nBased on the authors feedback, I'm inclined to keep my score of a weak accept. I think there are good ideas in the paper, but there are still some lingering questions (choice of polynomial basis, lack of shift invariance, etc) which degrade the practicality of the method.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your time reviewing our manuscript and for your valuable comments and suggestions. We will create a log/log plot, replace \"Theorem 3.1\" with \"Corollary 3.1.\" and mention that property of the graph Stein discrepancy.\n\n**(Lack of) confidence about power of PSD in the wild**\n\nWe have highlighted in the paper the strong theoretical performance on Gaussian targets and the strong empirical performance of PSD on non-Gaussian targets. During this review period, we have also highlighted issues with PSD for distributions where moments do not exist (e.g. Cauchy).  However, the development of theoretical results for the non-Gaussian setting represents a substantial undertaking that would be better considered in a future paper.\n\n**When to use PSD over KSD**\n\nBroadly we suggest using PSD when one wishes to obtain high statistical power for detecting discrepancies in moments and the target is \"close to Gaussian\" (with checks that we have expanded on elsewhere). Fast but biased sampling methods like SG-MCMC, where the discrepancy is often in the variance and the target is often close to Gaussian, are good applications that we highlight in the paper. The linear-time complexity of PSD is particularly helpful for large $n$.\n\n**How do practitioners know the posterior is Gaussian?**\n\nWe propose a diagnostic plot in the discussion. If the gradients are available in analytic form, rather than indirectly through automatic differentiation procedures, then one can also determine which expectations are being tracked by analysing the system of linear equations in equation (14) of Appendix B. We do this for the Rosenbrock target in Appendix C.6, and we are able to determine exactly which expectations are being tracked. We will suggest this technique in the discussion to diagnose situations where the PSD might not behave as expected.\n\n**Simpler ways to learn the moments if the target is Gaussian**\n\nEstimating the first $r$ moments need not be our only goal. However, checking for discrepancies in moments can be important in practice. For example, in the common application of SG-MCMC algorithms, poor step-sizes can lead to under- or over-estimated variance. Selecting the step-size based on matching the second-order moment can lead to a reasonable posterior approximation, as explained in the paper.\n\n**Moment-matching guarantees for a larger class of target distributions**\n\nWe can guarantee PSD is assessing moments for Gaussian targets because $\\nabla \\log p(x)$ is a linear function in that context, so the Stein operator applied to an $r$th order polynomial is itself an $r$th order polynomial. In this sense, the closer $\\nabla \\log p(x)$ is to a linear function, the closer the PSD is to assessing discrepancies in the first moments.\n\n**Sensitivity of PSD to mean-shift**\n\nWe are investigating this further by observing the sensitivity of PSD to transformations. The moment-tracking property of PSD will continue to hold due to similar reasoning as Appendix D.\n\n**Topology of PSD and properties when $PSD(P, Q_n)\\to 0$**\n\nWhile higher order PSD becomes increasingly complex, we can explicitly consider the case of $P=N(\\mu_P, \\Sigma_P)$ and $r = 1$. Here, $PSD(P,Q_n) = ||n^{-1}\\sum_{i=1}^n\\Sigma_P^{-1}(\\mu_P-X_i) ||_2 = ||\\Sigma_P^{-1}(\\mu_P-\\bar{X}) ||_2$ is measuring the Mahalanobis distance between $\\mu_P$ and the sample mean $\\bar{X}$. From the law of large numbers as $n\\rightarrow \\infty$, $\\bar{X}\\rightarrow \\mu_Q$ and consequently, $PSD(P,Q)\\rightarrow  \\Sigma_P^{-1}(\\mu_P-\\mu_Q)$ is zero if and only if the moments of $P$ and $Q$ match. Assuming iid draws from $Q$, the CLT gives $\\sqrt{n}\\Sigma_P^{-1}(\\bar{X}_Q-\\mu_q)\\sim N(0, \\Sigma_P^{-1}\\Sigma_Q\\Sigma_P^{-T})$ so $PSD = ||Y +\\Sigma_P^{-1}(\\mu_q-\\mu_P)||_2$ where $Y\\sim N(0, n^{-1}\\Sigma_P^{-1}\\Sigma_Q\\Sigma_p^{-T})$. The PSD behaves as the square root of a noncentral $\\chi^2$ distribution.\n\nThis result for $r=1$ serves as an illustrative example for iid samples and Gaussian targets. More practically useful theory, including for correlated samples, would be substantially more challenging and we will mention it as future work in the paper.\n\n**Simpler formulation for identifying specific moments where discrepancies occur**\n\nIf the target is Gaussian with diagonal covariance, then for the form of PSD we consider, the individual terms in the PSD correspond exactly to the moments of interest, as explained in the paper. This makes identification of problematic moments simpler than for alternative polynomial kernels. \n\n**Benefits of other spanning polynomials?**\n\nUnfortunately the benefits of orthogonal polynomials such as Chebyshev polynomials would likely be lost with PSD. It is not clear that one obtains a Chebyshev polynomial after applying a Langevin Stein operator to a Chebyshev polynomial.", "claims_and_evidence": "Yes, the claims in this paper are convincing. The main theoretical claim\n(Proposition 3.2) is proved in the appendix, and all other claims about its\nefficacy are demonstrated empirically. This is ample empirical work, which\nboth displays the benefits of this approach, and also some areas where the\nGaussian assumption of $P$ can go wrong.", "methods_and_evaluation_criteria": "Yes, the proposed methods make intuitive sense for the problem at\nhand. The authors even show that the PSD can perform well for distributions\nthat are bimodal (Figure 2), but that it can also fail for some non-Gaussian\nand heavy tailed targets (Appendix C.6).\n\nThey also demonstrate both theoretically and empirically that for low $r$,\nthis method is computationally more efficient than alternatives.", "theoretical_claims": "The authors show that PSD = 0 enforces moment matching up to order $r$ when\n$P$ is Gaussian. This illustrates that the PSD does at least what\npractitioners want in the most simple case. This was fully argued in\nAppendix B.", "experimental_designs_or_analyses": "The experimental analysis here seem valid. Most of them are borrowed from\nprevious papers in the field and thus have been used as benchmarks before.", "relation_to_broader_scientific_literature": "The paper is motivated by tailoring the Stein discrepancy to be\neven more useful than it currently is for Bayesian analysis. Since often\npractitioners often want their MCMC samplers to approximate a Bayesian\nposterior in the mean and covariance, this discrepancy is more tailored to\nhelping identify when that has gone awry.", "essential_references_not_discussed": "The paper did a good job citing the relevant literature.\n\nOne thing I might suggest is that in \"Measuring sample quality with Stein’s method\" paper, the graph Stein discrepancy actually does control convergence with respect to the Wasserstein metric, which implies convergence in mean.", "other_strengths_and_weaknesses": "The main strength of this paper is its signficance in articulating an\nefficient discrepancy that can help identify differences in moments when the\ntarget distribution is approximately Gaussian. The idea of using a Stein set\nof monomials is not so novel, but the thorough empirical work and ability to\nwork in a broad range of examples is quite remarkable.\n\nThe simplicity of the method and clarity of the paper are also positives.\n\nThere are some limitations to the paper. In the general theory of Stein discrepancies,\none often defines the suitable test set $\\mathcal{H}$ which controls convergence.\nE.g., in another version of the Stein discrepancy, this would be the set currently used\nas the span of monomials $\\mathcal{G}$. Then one usually solves the Stein equation\n$\\mathcal{A} g  = h$ so they can guarantee there is always a $g\\in\\mathcal{G}$ that\nsatisfies the Stein equation for each $h$ for some set $\\mathcal{G}$.\n\nThe authors here do implicitly show a version of this for multivariate Gaussian $P$, which is a notable feat. However,\nas discussed in the questions below, this does not naturally generalize to non-Gaussian\ndistributions. The authors do present some solid evidence that the PSD can still perform\nwell even for non-Gaussian $P$, but there are no robustness results to illustrate when\nthe PSD might behave quite surprisingly. And while the point about the Bernstein-von Moses\nlimit is well taken, it still doesn't provide much confidence how powerful this statistic is when\nusing the PSD in the wild.", "other_comments_or_suggestions": "+ A few places in the paper mention \"Theorem 3.1\" but in the paper this is\n  written as \"Corollary 3.1.\"\n\n+ It would be nice if Figure 4 was a log/log plot so it was easier to see the\ndifferent growth rates.", "questions_for_authors": "[Q1] Under what conditions should a practitioner elect to use the PSD over the\nKSD (with say IMQ kernel)? As surfaced in the paper, the IMQ-KSD at least controls\nweak convergence (but not necessarily moments), whereas the PSD controls moments in the case when $P$ is Gaussian.\nIf practitioners know that the posterior distribution is Gaussian, they can surely use\nthe PSD, but how often do practitioners know this in advance? If practitioners know\nthat the target is Gaussian, aren't there simpler ways to learn the moments of the\ntarget distribution? If the moment matching were guaranteed for a larger class of\ntarget distributions (e.g. log concave), then perhaps the validity of this would be\nmore concrete.\n\n[Q2] The PSD given here is not translation invariant, i.e., if one shifts\nboth P and Q by some vector $\\alpha$, the PSD will have a different\nvalue. It appears that the implicit choice of the origin imposes different\nweights on each basis monomials of $\\mathcal{G}$. How sensitive is the power of the\nPSD to this implicit parameterization? Should one try to center the\ndistributions $P$ and $Q$ before using the PSD?\n\n[Q3] The authors show that if $P$ is a multivariate Gaussian and $PSD(P, Q)=0$,\nthen the first $r$ moments of $P$ and $Q$ match. Can the authors say anything\nabout the topology of the PSD, even in the case that $P$ is Gaussian? I.e.,\nif $PSD(P, Q_n)\\to 0$, does this imply that the first r moments of $Q_n$ must converge\nto $P$?\n\n[Q4] The authors mention that the PSD presented here is \"offering a simpler formulation\nthat may be more effective for identifying specific moments where discrepancies occur.\"\nCan the authors explain this in more detail? In the non-Gaussian case this seems \nnon-trivial, and even in the Gaussian case, how would this look exactly? The differentiation operator and a non-identity covariance\nwould result in some weighted sum of monomials, which are not exactly known without knowing the covariance. Are there any\nbenefits to using a different set of spanning polynomials, e.g. the Chebyshev polynomials?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper proposes a novel discrepancy metric, the Polynomial Stein Discrepancy (PSD), as a scalable alternative to the Kernel Stein Discrepancy (KSD) for assessing the quality of samples from approximate Bayesian inference algorithms. By leveraging bounded-degree polynomial test functions, the method enables linear-time computation, making it appealing for large-scale applications. The paper combines theoretical grounding for Gaussian targets with empirical validation on both synthetic and real-world scenarios.
The strengths of the paper are its:
* Theoretical Contributions: The PSD’s ability to test moment convergence under Gaussian targets is formally established (Proposition 3.2), and the paper offers a nuanced discussion of its limitations and domain of applicability.
* Practical Impact: Reviewers appreciate the computational efficiency of PSD and its utility in tasks like hyperparameter tuning for sampling algorithms such as SG-MCMC.
There were, however, some shared limitations and areas for improvement noted:
* Generality: The method is strongly grounded in the Gaussian case. While it performs well empirically on some non-Gaussian targets, its theoretical robustness outside this setting remains unclear both empirically (ie on more examples) and theoretically. This is clearly the biggest weakness of the paper. 
* Theoretical Framing: Reviewer 8qrD argues that the framing via the Bernstein-von Mises (BvM) limit is somewhat overstated and potentially misleading, and recommends grounding the claims more directly in the Gaussian setting. The authors agreed and clarified this in the rebuttal.
* Comparisons: Reviewer ExjZ and others note that additional comparisons to other discrepancy-based and moment-based evaluation strategies (e.g., variational inference diagnostics) could strengthen the practical recommendations.

All reviewers had a positive opinion. with three giving a score of 4 (accept) and one giving a 3 (weak accept). While its main results are for Gaussian targets, the combination of theoretical guarantees and empirical performance makes it a valuable addition to the toolbox for assessing sample quality in approximate inference. The authors' responses addressed most of the reviewers' concerns, and the discussion around limitations and extensions (e.g., some non-Gaussian targets, choice of Stein set) was thoughtful and transparent. This paper is likely to be of interest to both theoretical and applied audiences in the Monte Carlo and Bayesian communities. Final Recommendation: Accept
Decision: Accept (poster)