Abstract:
Kernelized Stein discrepancy (KSD), though being extensively used in goodness-offit tests and model learning, suffers from the curse-of-dimensionality. We address this issue by proposing the sliced Stein discrepancy and its scalable and kernelized variants, which employ kernel-based test functions defined on the optimal one-dimensional projections. When applied to goodness-of-fit tests, extensive experiments show the proposed discrepancy significantly outperforms KSD and various baselines in high dimensions. For model learning, we show its advantages over existing Stein discrepancy baselines by training independent component analysis models with different discrepancies. We further propose a novel particle inference method called sliced Stein variational gradient descent (S-SVGD) which alleviates the mode-collapse issue of SVGD in training variational autoencoders.

heading: 1 INTRODUCTION
text: Discrepancy measures for quantifying differences between two probability distributions play key roles in statistics and machine learning. Among many existing discrepancy measures, Stein discrepancy (SD) is unique in that it only requires samples from one distribution and the score function (i.e. the gradient up to a multiplicative constant) from the other (Gorham & Mackey, 2015). SD, a special case of integral probability metric (IPM) (Sriperumbudur et al., 2009), requires finding an optimal test function within a given function family. This optimum is analytic when a reproducing kernel Hilbert space (RKHS) is used as the test function family, and the corresponding SD is named kernelized Stein discrepancy (KSD) (Liu et al., 2016; Chwialkowski et al., 2016). Variants of SDs have been widely used in both Goodness-of-fit (GOF) tests (Liu et al., 2016; Chwialkowski et al., 2016) and model learning (Liu & Feng, 2016; Grathwohl et al., 2020; Hu et al., 2018; Liu & Wang, 2016).
Although theoretically elegant, KSD, especially with RBF kernel, suffers from the ”curseof-dimensionality” issue, which leads to significant deterioration of test power in GOF tests (Chwialkowski et al., 2016; Huggins & Mackey, 2018) and mode collapse in particle inference (Zhuo et al., 2017; Wang et al., 2018). A few attempts have been made to address this problem, however, they either are limited to specific applications with strong assumptions (Zhuo et al., 2017; Chen & Ghattas, 2020; Wang et al., 2018) or require significant approximations (Singhal et al., 2019). As an alternative, in this work we present our solution to this issue by adopting the idea of “slicing”. Here the key idea is to project the score function and test inputs onto multiple one dimensional slicing directions, resulting in a variant of SD that only requires to work with one-dimensional inputs for the test functions. Specifically, our contributions are as follows.
• We propose a novel theoretically validated family of discrepancies called sliced Stein discrepancy (SSD), along with its scalable variant called max sliced kernelized Stein discrepancy (maxSKSD) using kernel tricks and the optimal test directions.
• A GOF test is derived based on an unbiased estimator of maxSKSD with optimal test directions. MaxSKSD achieves superior performance on benchmark problems and restricted Boltzmann machine models (Liu et al., 2016; Huggins & Mackey, 2018).
∗Work done at Microsoft Research Cambridge
• We evaluate the maxSKSD in model learning by two schemes. First, we train an independent component analysis (ICA) model in high dimensions by directly minimising maxSKSD, which results in faster convergence compared to baselines (Grathwohl et al., 2020). Further, we propose a particle inference algorithm based on maxSKSD called the sliced Stein variational gradient descent (S-SVGD) as a novel variant of the original SVGD (Liu & Wang, 2016). It alleviates the posterior collapse of SVGD when applied to training variational autoencoders (Kingma & Welling, 2013; Rezende et al., 2014).

heading: 2 BACKGROUND

heading: 2.1 KERNELIZED STEIN DISCREPANCY
text: For two probability distributions p and q supported on X ⊆ RD with continuous differentiable densities p(x) and q(x), we define the score sp(x) = ∇x log p(x) and sq(x) accordingly. For a test function f : X → RD, the Stein operator is defined as
Apf(x) = sp(x)T f(x) +∇Txf(x). (1) For a function f0 : RD → R, the Stein class Fq of q is defined as the set of functions satisfying Stein’s identity (Stein et al., 1972): Eq[sq(x)f0(x) +∇xf0(x)] = 0. This can be generalized to a vector function f : RD → RD where f = [f1(x), . . . , fD(x)]T by letting fi belongs to the Stein class of q for each i ∈ D. Then the Stein discrepancy (Liu et al., 2016; Gorham & Mackey, 2015) is defined as
D(q, p) = sup f∈Fq Eq[Apf(x)] = sup f∈Fq Eq[(sp(x)− sq(x))T f(x)]. (2)
When Fq is sufficiently rich, and q vanishes at the boundary of X , the supremum is obtained at f∗(x) ∝ sp(x)− sq(x) with some mild regularity conditions on f (Hu et al., 2018). Thus, the Stein discrepancy focuses on the score difference of p and q. Kernelized Stein discrepancy (KSD) (Liu et al., 2016; Chwialkowski et al., 2016) restricts the test functions to be in a D-dimensional RKHS HD with kernel k to obtain an analytic form. By defining up(x,x′) = sp(x)Tsp(x′)k(x,x′) + sp(x) T∇x′k(x,x′) + sp(x′)T∇xk(x,x′) + Tr(∇x,x′k(x,x′)) the analytic form of KSD is:
D2(q, p) = ( sup
f∈HD,||f ||HD≤1 Eq[Apf(x)]
)2 = Eq(x)q(x′)[up(x,x′)]. (3)

heading: 2.2 STEIN VARIATIONAL GRADIENT DESCENT
text: Although SD and KSD can be directly minimized for variational inference (VI) (Ranganath et al., 2016; Liu & Feng, 2016; Feng et al., 2017), Liu & Wang (2016) alternatively proposed a novel particle inference algorithm called Stein variational gradient descent (SVGD). It applies a sequence of deterministic transformations to a set of points such that each of mappings maximally decreases the Kullback-Leibler (KL) divergence from the particles’ underlying distribution q to the target p.
To be specific, we define the mapping T (x) : RD → RD as T (x) = x+ φ(x) whereφ characterises the perturbations. The result from Liu & Wang (2016) shows that the optimal perturbation inside the RKHS is exactly the optimal test function in KSD. Lemma 1. (Liu & Wang, 2016) Let T (x) = x + φ(x) and q[T ](z) be the density of z = T (x) when x ∼ q(x). If the perturbation φ is in the RKHSHD and ||φ||HD ≤ D(q, p), then the steepest descent directions φ∗q,p is φ∗q,p(·) = Eq[∇x log p(x)k(x, ·) +∇xk(x, ·)] (4) and ∇ KL[q[T ]||p]| =0 = −D2(q, p).
The first term in Eq.(4) is called drift, which drives the particles towards a mode of p. The second term controls the repulsive force, which spreads the particles around the mode. When particles stop moving, the KL decrease magnitude D2(q, p) is 0, which means the KSD is zero and p = q a.e.

heading: 3 SLICED KERNELIZED STEIN DISCREPANCY
text: We propose the sliced Stein discrepancy (SSD) and kernelized version named maxSKSD. Theoretically, we prove their correctness as discrepancy measures. Methodology-wise, we apply maxSKSD to GOF tests, and develop two ways for model learning.

heading: 3.1 SLICED STEIN DISCREPANCY
text: Before moving to the details, we give a brief overview of the intuition on how to tackle the curse-offimensionality issue of SD (The right figure of Figure 1). For detailed explanation, refer to appendix B.1. This issue of Stein discrepancy (Eq.2) comes from two sources: the score function sp(x) and the test function f(x) defined on X ⊂ RD. First, we notice that comparing sp and sq is equivalent to comparing projected score srp = s T p r and s r q for all r ∈ SD−1 on an hyper-sphere (Green square in Figure 1 (Right)). This operation reduces the test function’s output from RD to R (Green circle in Figure 1 (Right)). However, its input dimension is not affected. Reducing the input dimension of test functions is non-trivial, as directly removing input dimensions results in the test power decrease. This is because less information is accessed by the test function (see examples in appendix B.1). Our solution to this problem uses Radon transform which is inspired by CT-scans. It projects the original test function f(x) in Stein discrepancy (Eq. 2) (as an RD → R mapping) to a group of R → R functions along a set of directions (g ∈ SD−1). Then, this group of functions are used as the new test functions to define the proposed discrepancy. The invertibility of Radon transform ensures that testing with input in the original space RD is equivalent to the test using a group of low dimensional functions with input in R. Thus, the above two steps not only reduce the dimensions of the test function’s output and input, but also maintain the validity of the resulting discrepancy as each step is either equivalent or invertible.
In detail, assume two distributions p and q supported on RD with differentiable densities p(x) and q(x), and define the test functions f(·; r, g) : RD → R such that f(x; r, g) = frg ◦ hg(x) = frg(x
Tg), where hg(·) is the inner product with g and frg : R→ R. One should note that the r and g in f(·; r, g) should not just be treated as parameters in a test function f . In fact, they are more like the index to indicate that for each pair of r, g, we need a new f(·; r, g), i.e. new frg, which is completely independent to other test functions. The proposed sliced Stein discrepancy (SSD), defined using two uniform distributions pr(r) and pg(g) over the hypersphere SD−1, is given by the following, with frg ∈ Fq meaning f(·; r, g) ∈ Fq:
S(q, p) = Epr,pg
[ sup
frg∈Fq Eq[srp(x)frg(xTg) + rTg∇xT gfrg(xTg)]
] . (5)
We verify the proposed SSD is a valid discrepancy measure, namely, S(q, p) = 0 iff. q = p a.e. Theorem 1. (SSD Validity) If assumptions 1-4 in appendix A are satisfied, then for two probability distributions p and q, S(q, p) ≥ 0, and S(q, p) = 0 if and only if p = q a.e.
Despite this attractive theoretical result, SSD is difficult to compute in practice. Specifically, the expectations over r and g can be approximated by Monte Carlo but this typically requires a very
large number of samples in high dimensions (Deshpande et al., 2019). We propose to relax such limitations by using only a finite number of slicing directions r from an orthogonal basis Or of RD, e.g. the standard basis of one-hot vectors, and the corresponding optimal test direction gr for each r. We call this variant maxSSD, which is defined as follows and validated in Corollary 1.1:
Smax(q, p) = ∑ r∈Or sup frgr∈Fq,gr∈SD−1 Eq[srp(x)frgr (xTgr) + rTgr∇xT grfrgr (xTgr)]. (6)
Corollary 1.1. (maxSSD) Assume the conditions in Theorem 1, then Smax(q, p) = 0 iff. p = q a.e.

heading: 3.2 CLOSED FORM SOLUTION WITH THE KERNEL TRICK
text: The optimal test function given r and g is intractable without further assumptions on the test function families. This introduces another scalability issue as optimizing these test functions explicitly can be time consuming. Fortunately, we can apply the kernel trick to obtain its analytic form. Assume for each test function frg ∈ Hrg, where Hrg is a scalar-valued RKHS equipped with kernel k(x,x′; r, g) = krg(xTg,x′Tg) that satisfies assumption 5 in appendix A and frg(xTg) = 〈frg, krg(xTg, ·)〉Hrg . We define the following quantities:
ξp,r,g(x, ·) = srp(x)krg(xTg, ·) + rTg∇xT gkrg(xTg, ·), (7) hp,r,g(x,y) = s r p(x)krg(x Tg,yTg)srp(y) + r Tgsrp(y)∇xT gkrg(xTg,yTg)+
rTgsrp(x)∇yT gkrg(xTg,yTg) + (rTg)2∇2xT g,yT gkrg(xTg,yTg). (8) The following theorem describes the optimal test function inside SSD (Eq.(5)) and maxSSD (Eq.(6)). Theorem 2. (Closed form solution) If Eq[hp,r,g(x,x)] <∞, then
D2rg(q, p) = || sup frg∈Hrg,||frg||≤1 Eq[srp(x)frg(xTg) + rTg∇xT gfrg(xTg)]||2
= ||Eq[ξp,r,g(x)]||2Hrg = Eq(x)q(x′)[hp,r,g(x,x′)]. (9)
Next, we propose the kernelized version of SSD with orthogonal basis Or, called SKSD. Theorem 3. (SKSD as a discrepancy) For two probability distributions p and q, given assumptions 1,2 and 5 in appendix A and Eq[hp,r,g(x,x)] <∞ for all r and g, we define SKSD as
SKo(q, p) = ∑ r∈Or ∫ SD−1 pg(g)D 2 rg(q, p)dg, (10)
which is equal to 0 if and only if p = q a.e.
Following the same idea of maxSSD (Eq.6), it suffices to use optimal slice direction gr for each r ∈ Or, resulting in a slicing matrix G ∈ SD×(D−1). We name this discrepancy as maxSKSD, or maxSKSD-g when we need to distinguish it from another variant described later. Corollary 3.1. (maxSKSD) Assume the conditions in Theorem 3 are satisfied. Then
SKmax(q, p) = ∑ r∈Or sup gr D2rgr (q, p) (11)
is equal to 0 if and only if p = q a.e.
Figure 1 (Left) clarifies the connections between the mentioned discrepancies. We emphasise that using a single projection g in maxSKSD may be insufficient when no single projected feature xTg is informative enough to describe the difference between p and q. Instead, in maxSKSD, for each score projection r ∈ Or, we have a corresponding gr. One can also use the optimal r to replace the summation over Or, which provides additional benefits in certain GOF tests. We call this discrepancy maxSKSD-rg, and its validity can be proved accordingly. Interestingly, in appendix G, we show under certain scenarios maxSKSD-g can have inferior performance due to the noisy information provided by the redundant dimensions. Further, we show that such limitation can be efficiently addressed by using maxSKSD-rg.
Kernel choice and optimalG RBF kernel with median heuristics is a common choice. However, better kernels, e.g. deep kernels which evaluate a given kernel on the transformed input φ(x), might be preferred. It is non-trivial to directly use such kernel on SKSD or maxSKSD. We propose an adapted form of Eq.(10) to incorporate such kernel and maintain its validity. We include the details in appendix D and leave the experiments for future work.
The quality of sliced directionG is crucial for the performance of both maxSKSD-g or maxSKSD-rg. Indeed, it represents the projection directions that two distributions differ the most. The closed-form solutions of G is not analytic in general, in practice, finding the optimal G involves solving other difficult optimizations as well (projection r and test function frg). For the scope of this work, we obtained G by optimizing maxSKSD-g or maxSKSD-rg using standard gradient optimization, e.g. Adam, with random initialization. Still in some special cases (e.g. p, q are full-factorized), analytic solutions of optimalG exists, which is further discussed in appendix E.

heading: 3.3 APPLICATION OF MAXSKSD
text: Goodness-of-fit Test Assume the optimal test directions gr ∈ G are available, maxSKSD (Eq.(11)) can then be estimated using U-statistics (Hoeffding, 1992; Serfling, 2009). Given i.i.d. samples {xi}Ni=1 ∼ q, we have an unbiased minimum variance estimator:
SK ∧ max(q, p) = 1 N(N − 1) ∑ r∈Or ∑ 1≤i 6=j≤N hp,r,gr (xi,xj). (12)
The asymptotic behavior of the estimator is analyzed in appendix F.1. We use bootstrap (Liu et al., 2016; Huskova & Janssen, 1993; Arcones & Gine, 1992) to determine the threshold for rejecting the null hypothesis as indicated in algorithm 1. The bootstrap samples can be calculated by
SK ∧∗ m = ∑
1≤i 6=j≤N
(wmi − 1
N )(wmj −
1 N ) ∑ r∈Or hp,r,gr (xi,xj) (13)
where (wm1 , . . . , w m N ) M m=1 are random weights drawn from multinomial distributions Multi(N, 1N , . . . , 1 N ).
Algorithm 1: GOF Test with maxSKSD U-statistics Input :Samples {xi}Ni=1 ∼ q(x), score function sp(x), Orthogonal basis Or, optimal test
direction gr for each r ∈ Or, kernel function krg , significant level α, and bootstrap sample size M .
Hypothesis :H0: p = q v.s. H1: q 6= p Compute SK ∧ max(q, p) using U-statistic Eq.(12); Generate M bootstrap samples {SK ∧∗
m}Mm=1 using Eq.(13); Reject null hypothesis H0 if the proportion SK ∧∗ m > SK ∧ max(q, p) is less than α.
Model Learning The proposed maxSKSD can be applied to model learning in two ways. First, it can be directly used as a training objective, in such case q is the data distribution and p is the model to be learned, and the learning algorithm performs minp SKmax(q, p). The second model learning scheme is to leverage the particle inference for latent variables and train the model parameters using an EM-like (Dempster et al., 1977) algorithm. Similar to the relation between SVGD and KSD, we can derive a corresponding particle inference algorithm based on maxSKSD, called sliced-SVGD (S-SVGD). In short, we define a specific form of the perturbation as φ(x) = [φgi(x Tgi), . . . , φgD (x TgD)]
T and modify the proofs of Lemma 1 accordingly. The resulting S-SVGD algorithm uses kernels defined on one dimensional projected samples, which sidesteps the vanishing repulsive force problem of SVGD in high dimensions (Zhuo et al., 2017; Wang et al., 2018). We illustrate this in Figure 2 by estimating the variance
of a standard Gaussian with the particles obtained by SVGD or S-SVGD (see appendix J.1). We see that as the dimension increases, SVGD severely under-estimates the variance of p, while the S-SVGD remains robust. Furthermore, its validity is justified since in such case the KL gradient equals to maxSKSD which is a valid discrepancy. Readers are referred to appendix F.2 for the derivations. We also give an analysis of their memory and computational cost for both GOF and model learning in appendix H.

heading: 4 EXPERIMENTS

heading: 4.1 GOODNESS OF FIT TEST
text: We evaluate maxSKSD (Eq.(11)) for GOF tests in high dimensional problems. First, we demonstrate its robustness to the increasing dimensionality using the Gaussian GOF benchmarks (Jitkrittum et al., 2017; Huggins & Mackey, 2018; Chwialkowski et al., 2016). Next, we show the advantage of our method for GOF tests on 50-dim Restricted Boltzmann Machine (RBM) (Liu et al., 2016; Huggins & Mackey, 2018; Jitkrittum et al., 2017). We included in comparison extensive baseline test statitics for GOF test: Gaussian or Cauchy random Fourier features (RFF) (Rahimi & Recht, 2008), KSD with RBF kernel (Liu et al., 2016; Chwialkowski et al., 2016), finite set Stein discrepancy (FSSD) with random or optimized test locations (Jitkrittum et al., 2017), random feature Stein discrepancy (RFSD) with L2 SechExp and L1 IMQ kernels (Huggins & Mackey, 2018), and maximum mean discrepancy (MMD) (Gretton et al., 2012) with RBF kernel. Notice that we use gradient descent to obtain the test directions gr (and potentially the slicing directions r) for Eq.(11).

heading: 4.1.1 GOF TESTS WITH HIGH DIMENSIONAL GAUSSIAN BENCHMARKS
text: We conduct 4 different benchmark tests with p = N (0, I): (1) Null test: q = p; (2) Laplace: q(x) = ∏D d=1 Lap(xd|0, 1/ √ 2) with mean/variance matched to p; (3) Multivariate-t: q is fully factorized multivariate-t with 5 degrees of freedom, 0 mean and scale 1. In order to match the variance of p and q, we change the variance of p to 55−2 ; (4) Diffusion: q(x) = N (0,Σ1) where the variance of 1st-dim is 0.3 and the rest is the same as in I . For the testing setup, we set the significance level α = 0.05. For FFSD and RFSD, we use the open-sourced code from the original publications. We only consider maxSKSD-g here as it already performs nearly optimally. We refer to appendix I.1 for details.
Figure 3 shows the GOF test performances and the corresponding discrepancy values. In summary, the proposed maxSKSD outperforms the baselines in all tests, where the result is robust to the increasing dimensions and the discrepancy values match the expected behaviours.
Null The left-most column in Figure 3 shows that all methods behave as expected as the rejection rate is closed to the significance level, except for RFSD with L2 SechExp kernel. All the discrepancy values oscillate around 0, with the KSD being less stable.
Laplace and Multivariate-t The two middle columns of Figure 3 show that maxSKSD-g achieves a nearly perfect rejection rate consistently as the dimension increases, while the test power for all
baselines decreases significantly. For the discrepancy values, similar to the KL divergence between q and p, maxSKSD-g linearly increases with dimensions due to the independence assumptions..
Diffusion This is a more challenging setting since p and q only differ in one of their marginal distributions, which can be easily buried in high dimensions. As shown in the rightmost column of Figure 3, all methods failed in high dimensions except maxSKSD-g, which still consistently achieves optimal performance. For the discrepancy values, we expect a positive constant due to the one marginal difference between p and q. Only maxSKSD-g behaves as expected as the problem dimension increases. The decreasing value at the beginning is probably due to the difficulty in finding the optimal direction g in high dimensions when the training set is small.
4.1.2 RBM GOF TEST
We demonstrate the power maxSKSD for GOF tests on RBMs, but we now also include results for maxSKSD-rg. We follow the test setups in Liu et al. (2016); Jitkrittum et al. (2017); Huggins & Mackey (2018) where different amounts of noise are injected into the weights to form the alternative hypothesis q. The samples are drawn using block Gibbs samplers. Refer to appendix I.2 for details. Figure 4 shows that maxSKSD based methods dominate the baselines, especially with maxSKSD-rg significantly outperforming the others. At perturbation level 0.01, maxSKSD-rg achieves 0.96 rejection rate, while others are all below 0.5. This result shows the advantages of optimizing the slicing directions r.

heading: 4.2 MODEL LEARNING
text: We evaluate the efficiency of maxSKSD-based algorithms in training machine learning models. First, we use independent component analysis (ICA) which is often used as a benchmark for evaluating training methods for energy-based model (Gutmann & Hyvärinen, 2010; Hyvärinen, 2005; Ceylan & Gutmann, 2018). Our approach trains the ICA model by directly minimizing maxSKSD. Next, we evaluate the proposed S-SVGD particle inference algorithm, when combined with amortization (Feng et al., 2017; Pu et al., 2017), in the training of a variational autoencoder (VAE) (Kingma & Welling, 2013; Rezende et al., 2014) on binarized MNIST. Appendix J.5 also shows superior results for S-SVGD when training a Bayesian neural network (BNN) on UCI datasets (Dua & Graff, 2017).

heading: 4.2.1 ICA
text: ICA consists of a simple generative process z ∼ Lap(0, 1) andx =Wz, where the model parameters are a non-singular matrix W ∈ RD×D. The log density for x is log p(x) = log pz(W−1x) + C, where the normalization constant C can be ignored when training with Stein discrepancies. We train the models on data sampled from a randomly initialized ICA model and evaluate the corresponding test log likelihoods. We compare maxSKSD with KSD and the state-of-the-art LSD (Grathwohl et al., 2020). For more details on the setup, we refer the reader to appendix J.2.
Table 1 shows that both maxSKSD and LSD are robust to increasing dimensions, with maxSKSD being better when D is very large. Also at D = 200, maxSKSD converges significantly faster than LSD (see Figure 10 in appendix J.3). This faster convergence is due to the closed-form solution for the optimal test functions, whereas LSD requires adversarial training. While KSD is also kernel-based, it suffers from the curse-of-dimensionality and fails to train the model properly for D > 20. Instead the proposed maxSKSD can successfully avoid the problems of KSD with high dimensional data.
Table 3: Label entropy and accuracy for imputed images.
Method Entropy Accuracy
Vanilla VAE 0.297 0.718 SVGD VAE 0.538 0.691 S-SVGD VAE 0.542 0.728

heading: 4.2.2 AMORTIZED SVGD
text: Finally, we consider training VAEs with implicit encoders on dynamically binarized MNIST. The decoder is trained as in vanilla VAEs, but the encoder is trained by amortization (Feng et al., 2017; Pu et al., 2017), which minimizes the mean square error between the initial samples from the encoder, and the modified samples driven by the SVGD/S-SVGD dynamics (Algorithm 3 in appendix J.4).
We report performance in terms of test log-likelihood (LL). Furthermore we consider an imputation task, by removing the pixels in the lower half of the image and imputing the missing values using (approximate) posterior sampling from the VAE models. The performance is measured in terms of imputation diversity and correctness, using label entropy and accuracy. For fair comparisons, we do not tune the coefficient of the repulsive force. We refer to appendix J.4 for details.
Table 2 reports the average test LL. We observe that S-SVGD is much more robust to the increasing latent dimensions compared to SVGD. To be specific, with D = 16, SVGD performs the best where S-SVGD performs slightly worse than SVGD. However, when the dimension starts to increase, LL of SVGD drops significantly. For D = 64, a common choice for latent space, it performs even significantly worse than vanilla VAE. On the other hand, S-SVGD is much more robust. Notice that the purpose of this experiment is to show compare their robustness instead of achieving the state-ofthe-art performance. Still the performance can be easily boosted, e.g. running longer S-SVGD steps before encoder update, we leave it for the future work.
For the imputation task, we compute the label entropy and accuracy for the imputed images (Table 3). We observe S-SVGD has higher label entropy compared to vanilla VAE and better accuracy compared to SVGD. This means both S-SVGD and SVGD capture the muli-modality nature of the posterior compared to uni-modal Gaussian distribution. However, high label entropy itself may not be a good indicator for the quality of the learned posterior. One can think of a counter-example that the imputed images are diverse but does not look like any digits. This may also gives a high label entropy but the quality of the posterior is poor. Thus, we use the accuracy to indicate the “correctness” of the imputed images, with higher label accuracy meaning the imputed images are closed to the original image. Together, a good model should give a higher label entropy along with the high label accuracy. We observe S-SVGD has more diverse imputed images with high imputation accuracy.

heading: 4.3 SUMMARY OF THE EXPERIMENTS IN APPENDIX
text: We present further empirical results on GOF tests and model learning in the appendix to demonstrate the advantages of the proposed maxSKSD. As a summary glance of the results:
• In appendix G, we analyse the potential limitations of maxSKSD-g and show that they can be mitigated by maxSKSD-rg, i.e. optimising the slicing direction r;
• In appendix I.3, we successfully apply maxSKSD to selecting the step size for stochastic gradient Hamiltonian Monte Carlo (SGHMC) (Chen et al., 2014);
• In appendix J.5, we show that the proposed S-SVGD approach out-performs the original SVGD on Bayesian neural network regression tasks.

heading: 5 RELATED WORK
text: Stein Discrepancy SD (Gorham & Mackey, 2015) and KSD (Liu et al., 2016; Chwialkowski et al., 2016) are originally proposed for GOF tests. Since then research progress has been made to improve these two discrepancies. For SD, LSD (Grathwohl et al., 2020; Hu et al., 2018) is proposed to increase the capacity of test functions using neural networks with L2 regularization. On the other hand, FSSD (Jitkrittum et al., 2017) and RFSD (Huggins & Mackey, 2018) aim to reduce
the computation cost of KSD from O(n2) to O(n) where n is the number of samples. Still the curse-of-dimensionality issue remains to be addressed in KSD, and the only attempt so far (to the best of our knowledge) is the kernelized complete conditional Stein discrepancy (KCC-SD (Singhal et al., 2019)), which share our idea of avoiding kernel evaluations on high dimensional inputs but through comparing conditional distributions. KCC-SD requires the sampling from q(xd|x−d), which often needs significant approximations in practice due to its intractability. This makes KCC-SD less suited for GOF test due to estimation quality in high dimensions. On the other hand, our approach does not require this approximation, and the corresponding estimator is well-behaved asymptotically.
Wasserstein Distance and Score matching Sliced Wasserstein distance (SWD) (Kolouri et al., 2016) and sliced score matching (SSM) (Song et al., 2019) also uses the “slicing” idea. However, their motivation is to address the computational issues rather than statistical difficulties in high dimensions. SWD leveraged the closed-form solution of 1D Wasserstein distance by projecting distributions onto 1D slices. SSM uses Hutchson’s trick (Hutchinson, 1990) to approximate the trace of Hessian.
Particle Inference Zhuo et al. (2017); Wang et al. (2018) proposed message passing SVGD to tackle the well-known mode collapse problem of SVGD using local kernels in the graphical model. However, our work differs significantly in both theory and applications. Theoretically, the discrepancy behind their work is only valid if p and q have the same Markov blanket structure (refer to Section 3 in Wang et al. (2018) for detailed discussion). Thus, unlike our method, no GOF test and practical inference algorithm can be derived for generic cases. Empirically, the Markov blanket structure information is often unavailable, whereas our method only requires projections that can be easily obtained using optimizations. Projected SVGD (pSVGD) is a very recent attempt (Chen & Ghattas, 2020) which updates the particles in an adaptively constructed low dimensional space, resulting in a biased inference algorithm. The major difference compared to S-SVGD is that our work still updates the particles in the original space with kernel being evaluated in 1D projections. Furthermore, S-SVGD can theoretically recover the correct target distribution. There is no real-world experiments provided in (Chen & Ghattas, 2020), and a stable implementation of pSVGD is non-trivial, so we did not consider pSVGD when selecting the baselines.

heading: 6 CONCLUSION
text: We proposed sliced Stein discrepancy (SSD), as well as its scalable and kernelized version maxSKSD, to address the curse-of-dimensionality issues in Stein discrepancy. The key idea is to project the score function on one-dimensional slices and define (kernel-based) test functions on one-dimensional projections. We also theoretically prove their validity as a discrepancy measure. We conduct extensive experiments including GOF tests and model learning to show maxSKSD’s improved performance and robustness in high dimensions. There are three exciting avenues of future research. First, although validated by our theoretical study in appendix D, practical approaches to incorporate deep kernels into SSD remains an open question. Second, the performance of maxSKSD crucially depends on the optimal projection direction, so better optimization methods to efficiently construct this direction is needed. Lastly, we believe “slicing” is a promising direction for kernel design to increase the robustness to high dimensional problems in general. For example, MMD can be easily extended to high dimensional two-sample tests using this kernel design trick.

Review summaries:
This paper tries to solve the curse-of-dimensionality problem of KSD and corresponding mode-collapse problem of SVGD by projecting both the input and output of test function onto 1D slices. By doing so, the paper proposes the new discrepancies called SSD and maxSKSD, and a new variant of SVGD called S-SVGD. Experiments on goodness-of-fit test (synthetic high-dim Gaussian & RBM) and model learning (ICA on synthetic data & amortized SVGD on MNIST) are reported in the main body of the paper.