    # Approximate Bayesian Computation with the Sliced-Wasserstein Distance

Approximate Bayesian Computation (ABC) is a popular method for approximate inference in generative models with intractable but easy-to-sample likelihood. It constructs an approximate posterior distribution by finding parameters for which the simulated data are close to the observations in terms of summary statistics. These statistics are defined beforehand and might induce a loss of information, which has been shown to deteriorate the quality of the approximation. To overcome this problem, Wasserstein-ABC has been recently proposed, and compares the datasets via the Wasserstein distance between their empirical distributions, but does not scale well to the dimension or the number of samples. We propose a new ABC technique, called Sliced-Wasserstein ABC and based on the Sliced-Wasserstein distance, which has better computational and statistical properties. We derive two theoretical results showing the asymptotical consistency of our approach, and we illustrate its advantages on synthetic data and an image denoising task.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Consider the problem of estimating the posterior distribution of some model parameters

given data points

. This distribution has a closed-form expression given by the Bayes’ theorem up to a multiplicative constant:

 π(θ|y1:n)∝π(y1:n|θ)π(θ).

For many statistical models of interest, the likelihood cannot be numerically evaluated in a reasonable amount of time, which prevents the application of classical likelihood-based approximate inference methods. Nevertheless, in various settings, even if the associated likelihood is numerically intractable, one can still generate synthetic data given any model parameter value. This generative setting gave rise to an alternative framework of likelihood-free inference techniques. Among them, Approximate Bayesian Computation (ABC, [1, 2]) methods have become a popular choice and have proven useful in various practical applications (e.g., [3, 4, 5]

). The core idea of ABC is to bypass calculation of the likelihood by using simulations: the exact posterior is approximated by retaining the parameter values for which the synthetic data are close enough to the observations. Closeness is usually measured with a discrepancy measure between the two datasets reduced to some ‘summary statistics’ (e.g., empirical mean or empirical covariance). While summaries allow a practical and efficient implementation of ABC, especially in high-dimensional data spaces, the quality of the approximate posterior distribution highly depends on them and constructing sufficient statistics is a non-trivial task. Summary statistics can be designed by hand using expert knowledge, which can be tedious especially in real-world applications, or in an automated way, for instance see

.

Recently, discrepancy measures that view data sets as empirical probability distributions to eschew the construction of summary statistics have been proposed for ABC. Examples include the Kullback-Leibler divergence (KL,

), maximum mean discrepancy , and Wasserstein distance 

. This latter distance emerges from the optimal transport (OT) theory and has attracted abundant attention in statistics and machine learning due to its strong theoretical properties and applications on many domains. In particular, it has the ability of making meaningful comparisons even between probability measures with non-overlapping supports, unlike KL. However, the computational complexity of the Wasserstein distance rapidly becomes a challenge when the dimension of the observations is large. Several numerical methods have been proposed during the past few years to speed-up this computation. Wasserstein-ABC (WABC,



) uses an approximation based on the Hilbert space-filling curve and termed the Hilbert distance, which is computationally efficient but accurate for small dimensions only. Besides, under a general setting, the Wasserstein distance suffers from a curse of dimensionality in the sense that the error made when approximating it from samples grows exponentially fast with the data space dimension

. These computational and statistical burdens can strongly affect the performance of WABC applied to high-dimensional data.

The Sliced-Wasserstein distance (SW, [11, 12]) is an alternative OT distance and leverages the attractive property that the Wasserstein distance between one-dimensional measures has an analytical form which can be efficiently approximated. SW is defined as an average of one-dimensional Wasserstein distances, and therefore has a significantly better computational complexity. Several recent studies have reported empirical success on generative modeling with SW [13, 14, 15] as well as nice asymptotic and statistical properties [16, 17, 18], making this alternative distance increasingly popular. In this paper, we develop a novel ABC framework that uses SW as the data discrepancy measure. This defines a likelihood-free method which does not require choosing summary statistics and is efficient even with high-dimensional observations. We derive asymptotical guarantees on the convergence of the resulting ABC posterior, and we illustrate the superior empirical performance of our methodology by applying it on a synthetical problem and an image denoising task.

## 2 Background

Consider a probability space with associated expectation operator

, on which all the random variables are defined. Let

be a sequence of independent and identically distributed (i.i.d.) random variables associated with some observations valued in . Denote by the common distribution of and by the set of probability measures on . For any , denotes the empirical distribution corresponding to observations. Consider a statistical model parametrized by . We focus on parameter inference for purely generative models: for any , we can draw i.i.d. samples from , but the numerical evaluation of the likelihood is not possible or too expensive. For any , is the empirical distribution of i.i.d. samples generated by , . We assume that , endowed with the Euclidean distance , is a Polish space (i.e., complete and separable), , endowed with the distance , is a Polish space, parameters are identifiable, i.e.  implies . denotes the Borel -field of . Approximate Bayesian Computation. ABC methods are used to approximate the posterior distribution in generative models when the likelihood is numerically intractable but easy to sample from. The basic and simplest ABC algorithm is an acceptance-rejection method , which iteratively draws a candidate parameter from a prior distribution , and ‘synthetic data’ from , and keeps if is close enough to the observations . Specifically, the acceptance rule is , where is a data discrepancy measure taking non-negative values, is a tolerance threshold, and with small is a summary statistics. The algorithm is summarized in Pseudo-code 1 and returns samples of that are distributed from:
(1)
The choice of directly impacts the quality of the resulting approximate posterior: if the statistics are sufficient statistics, converges to the true posterior as , otherwise, the limiting distribution is at best [19, 20]. Wasserstein-ABC has been proposed to avoid this loss of information. Wasserstein distance and ABC. For , consider for some and the Wasserstein distance of order defined for any by,
(2)
where is the set of probability measures on that verifies: , , .
Evaluating the Wasserstein distance between multi-dimensional probability measures turns out to be numerically intractable in general, and solving (2) between empirical distributions corresponding to samples leads to computational costs in . Nevertheless, between one-dimensional measures has a closed-form expression [22, Theorem 3.1.2.(a)], given by:
(3)
where and

denote the quantile functions of

and respectively. For empirical one-dimensional distributions, (3) can be efficiently approximated by simply sorting the samples drawn from each distribution and computing the average cost between the sorted samples. This amounts to operations at worst.
Wasserstein-ABC  is a variant of ABC (1) that uses , between the empirical distributions of the observed and synthetic data, in place of the discrepancy measure between summaries. To make this method scalable to any dataset size,  introduces a new approximation of (2), the Hilbert distance, which extends the idea behind the computation of in 1D to higher dimensions, by sorting samples according to their projection obtained via the Hilbert space-filling curve. This alternative can be computed in , but yields accurate approximations only for low dimensions . They also use a second approximation, the swapping distance, based on an iterative greedy swapping algorithm. However, each iteration requires operations, and there is no guarantee of convergence to . Input: observations , number of iterations , data discrepancy measure , summary statistics , tolerance threshold . for  do        repeat              and i.i.d.       until ;        return Pseudo-code 1 Vanilla ABC.

## 3 Sliced-Wasserstein ABC

Sliced-Wasserstein distance. The analytical expression of in (3) motivates the formulation of an alternative OT distance, called the Sliced-Wasserstein distance [11, 12]. SW is obtained by reducing multi-dimensional distributions to one-dimensional representations through linear projections, and then by averaging 1D Wasserstein distances between these projected distributions. More formally, we denote by the -dimensional unit sphere, and by the Euclidean inner-product. For any , is the linear form associated with such that for any . For , the Sliced-Wasserstein distance of order between is defined as,

 SWpp(μ,ν)=∫Sd−1Wpp(u⋆♯μ,u⋆♯ν)dσ(u), (4)

where

is the uniform distribution on

and for any measurable function and , is the push-forward measure of by , defined as: , , with . is a distance on  with significantly lower computational requirements than the Wasserstein distance: in practice, the integration in (4) is approximated with a finite-sample average, using a simple Monte Carlo scheme.

SW also seems to have better statistical properties than the Wasserstein distance and its approximations. We illustrate it with the task of estimating the scaling factor of the covariance matrix in a multivariate Normal model, as in the supplementary material of . Denote for any , the

-dimensional Gaussian distribution with zero-mean and covariance matrix

. Observations are assumed i.i.d. from with , and we draw the same number of i.i.d. data from for 100 values of equispaced between 0.1 and 9. We then compute and between the empirical distributions of the samples, and the swapping and Hilbert approximations presented in , for and 1000 observations. between two Gaussian measures has an analytical formula, which boils down in our setting to: , and we approximate the exact SW using a Monte Carlo approximation of:

 SW22(μσ⋆,μσ)=W22(μσ⋆,μσ)∫Sd−1uTu dσ(u),

This formula is derived from (4) and the exact between one-dimensional Gaussian distributions. We also compute KL with the estimator proposed for KL-based ABC (KL-ABC, ). Figure 1 shows the distances plotted against for each . When the dimension increases, we observe that (i) as pointed out in , the quality of the approximation of empirical Wasserstein returned by Hilbert and swapping rapidly deteriorates, and (ii) SW, approximated using densities or samples, is the only approximate metric that attains its minimum at . This curse of dimensionality can be a limiting factor for the performance of WABC and KL-ABC in high dimensions.

Sliced-Wasserstein ABC. Motivated by the practical success of SW regardless of the dimension value in the previous experiment, we propose a variant of ABC based on SW, referred to as the Sliced-Wasserstein ABC (SW-ABC). Our method is similar to WABC in the sense that it compares empirical distributions, but instead of , we choose the discrepancy measure to be , . The usage of SW allows the method to scale better to the data size and dimension. The resulting posterior distribution, called the SW-ABC posterior, is thus defined in (1) with replaced by . Figure 1: Comparison of OT distances and KL between data generated from d-dimensional Gaussian distributions μσ vs. μσ⋆, σ2⋆=4, with 1000 i.i.d draws. SW is approximated with 100 random projections.

## 4 Theoretical Study

In this section, we analyze the asymptotic behavior of the SW-ABC posterior under two different regimes. Our first result concerns the situation where the observations are fixed, and goes to zero. We prove that the SW-ABC posterior is asymptotically consistent in the sense that it converges to the true posterior, under specific assumptions on the density used to generate synthetic data.

###### Proposition (Asymptotic consistency when ε→0, y1:n fixed).

Let . Suppose that has a density w.r.t. the Lebesgue measure such that is continuous and there exists satisfying and . In addition, assume that there exists such that , where . Then, with fixed, the SW-ABC posterior converges to the true posterior as goes to 0, in the sense that, for any measurable , , where is defined by (1).

###### Proof of Section 4.

The proofs consists in applying [9, Proposition 3.1], which establishes the conditions for the data discrepancy measure to yield an ABC posterior that converges to the true posterior in the asymptotic regime we consider. This amounts to verify that:

1. [wide, labelwidth=!, labelindent=0pt,label=()]

2. For any and , with respective empirical distributions and , if and only if .

3. is continuous in the sense that, if converges to in the metric , then, for any empirical distribution ,  , where is the empirical measure of .

Condition 1 follows from the fact that is a distance [16, Proposition 5.1.2]. Now, let and be a continuous function such that for any with . Since converges to in the metric and is continuous, we get that . This implies that weakly converges to in [23, Definition 6.7], which, by [23, Theorem 6.8], is equivalent to . By applying the triangle inequality and [16, Proposition 5.1.3], there exists such that, for any empirical measure ,

 |SWp(^μn,^μkθ,m)−SWp(^μn,^μθ,m)| ≤SWp(^μkθ,m,^μθ,m) ≤C1/p Wp(^μkθ,m,^μθ,m).

We conclude that , making condition 2 applicable. ∎

Next, we study the limiting SW-ABC posterior when the value of is fixed and the number of observations increases, i.e. . We suppose the size of the synthetic dataset grows to with , such that can be written as a function of , , satisfying . We show that, under this setting and appropriate additional conditions, the resulting approximate posterior converges to the prior distribution on restricted to the region .

###### Proposition (Asymptotic consistency when ε fixed, n→∞, m/n→α>0).

Let , and be an increasing sequence satisfying , for . Assume that the statistical model is well specified, i.e. there exists such that , and that almost surely the following holds:

 limn→∞SWp(^μn,^μθ,m(n))=SWp(μθ⋆,μθ), (5)

where , denote the empirical distributions of the observations and synthetic data respectively. Then, the SW-ABC posterior converges to the restriction of the prior on the region as , i.e. for any ,

 limn→∞πεy1:n(θ) =π(θ|SWp(μθ⋆,μθ)≤ε) ∝π(θ)1{SWp(μθ⋆,μθ)≤ε}.
###### Proof of Section 4.

The result follows from the application of [7, Theorem 1] to and the required conditions. ∎

###### Remark 1.

Condition (5) is a mild assumption, e.g. is fulfilled when is compact and separable: in this case, for any and its empirical instantiation over samples, -almost surely , and by [16, Proposition 5.1.3], -almost surely. Equation 5 then follows by applying the triangle inequality.

## 5 Experiments Figure 2: Comparison of SMC-ABC strategies in the multivariate Gaussians problem. Each strategy uses 1000 particles and are run for 3 hours max. First row shows ABC posteriors against true posterior for σ2, second row reports W1-distance to true posterior versus time.

Multivariate Gaussians. As a first set of experiments, we investigate the performance of SW-ABC on a synthetical setting where the posterior distribution is analytically available. We consider observations i.i.d. from a -dimensional Gaussian , with and . The parameter is

for which the prior distribution is assigned to be an inverse gamma distribution

. Therefore, the posterior distribution of given and is .

We compare SW-ABC against ABC using the Euclidean distance between sample variances (Euclidean-ABC), WABC with the Hilbert distance, WABC with the swapping distance and KL-ABC. Each ABC approximation was obtained using the sequential Monte Carlo sampler-based ABC method

, which is more computationally efficient than vanilla ABC (1) and implemented in the package pyABC . We provide our code in . Figure 2 reports for , the resulting ABC posteriors and to the true posterior vs. time. Due to the poor performance of the estimator of KL between two empirical distributions proposed in  (see Fig. 1), KL-ABC fails at approximating well the posterior in these experiments. Hence, we excluded those results from Fig. 2 for clarity. Euclidean-ABC provides the most accurate approximation as expected since it relies on statistics that are sufficient in our setting. WABC performs poorly with high-dimensional observations, contrary to SW-ABC, which approximates well the posterior for each dimension value and is as fast.

Image denoising. We now evaluate our approach on a real application, namely image denoising. We consider a widely used algorithm for this task, the Non-Local Means algorithm (NL-means, ), and we present a novel variant of it derived from SW-ABC.

We formally define the denoising problem: let , denote a clean gray-level image. We observe a corrupted version of this image, , where is some noise in . The goal is to restore from . We focus on denoising methods that consider ‘patch-based representations’ of images, e.g. NL-means. Specifically, let be a patch size and the set of pixel positions. For , denotes the pixel value at position in image , and is a window in centered at : for ,  , where is extended to by periodicity. Let be a dictionary of positions, and such that, for , , i.e.  is the position in of the most similar patch to . For , an estimator of is given by , being the uniform distribution on . In practice, it is approximated with a Monte Carlo scheme:

 ^Pj≈(Tn)−1∑Tt=1∑ns=1Pi(t)+l(t,s), (6)

where , , and , are mutually independent. Finally, we construct an estimate of as follows: for any ,   . The classical NL-means estimator corresponds to the case where (thus ) and for any and , , where is a search window.

In our work, we assume that the likelihood is not available, but we observe for , i.i.d. from . By replacing in (6) by the SW-ABC posterior, we obtain the proposed denoising method, called the SW-ABC NL-means algorithm. We provide our Python implementation in .

We compare our approach with the classical NL-means. We consider one of the standard image denoising datasets (cf. ), namely the CBSD68 dataset . It consists of colored images of size

. We first convert the images to gray scale, then manually corrupt each image with a Gaussian noise with standard deviation

, and finally try to recover the clean image. The quality of the output images is evaluated with the Peak Signal to Noise Ratio (

) measure:  . In all of our experiments, we use a dictionary of patches picked uniformly at random, and we set , , , .

We report the average PSNR values for different values of the noise level in Table 1. We observe that for small , NL-means provides more accurate results, whereas when becomes larger SW-ABC outperforms NL-means, thanks to the patch representation and the use of the SW distance. On the other hand, another important advantage of SW-ABC becomes prominent in the computation time: the proposed approach takes on a standard laptop computer per image whereas the classical NL-means algorithm takes . Indeed, the computational complexity of SW-ABC NL-means is upper-bounded by , whereas it is for the naïve implementation of NL-means, where denotes the cardinal number of the set . By taking , we can directly observe that SW-ABC has a lower computational complexity since in practice. We note that the computation time of NL-means can be improved by certain acceleration techniques, which can be directly used to improve the speed of SW-ABC NL-means as well. Finally, in Figure 3, we illustrate the performance of SW-ABC on two images for visual inspection. The results show that the injected noise is successfully removed by the proposed approach.

## 6 Conclusion

We proposed a novel ABC method, SW-ABC, based on the SW distance. We derived asymptotic guarantees for the convergence of the SW-ABC posterior to the true posterior under different regimes, and we evaluated our approach on a synthetical and an image denoising problem. Our results showed that SW-ABC provides an accurate approximation of the posterior, even with high-dimensional data spaces and a small number of samples. Future work includes extending SW-ABC to generalized SW distances .

## Acknowledgments

This work is partly supported by the French National Research Agency (ANR) as a part of the FBIMATRIX project (ANR-16-CE23-0014) and by the Data Science & AI chair from Télécom Paris. Alain Durmus acknowledges support from Polish National Science Center grant: NCN UMO-2018/31/B/ST1/00253.

## References

•  S. Tavaré, D. J. Balding, R. C. Griffiths, and P. Donnelly. Inferring coalescence times from dna sequence data. Genetics, 145(2):505–518, 1997.
•  M. A. Beaumont, W. Zhang, and D. J. Balding. Approximate bayesian computation in population genetics. Genetics, 162(4):2025–2035, 2002.
•  G. Peters and S. Sisson. Bayesian inference, monte carlo sampling and operational risk. Journal of Operational Risk, 1, 09 2006.
•  M. Tanaka, A. Francis, F. Luciani, and S. Sisson. Using approximate bayesian computation to estimate tuberculosis transmission parameters from genotype data. Genetics, 173:1511–20, 08 2006.
•  S. N. Wood. Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466(7310):1102–1104, 8 2010.
•  P. Fearnhead and D. Prangle. Constructing summary statistics for approximate bayesian computation: semi-automatic approximate bayesian computation. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(3):419–474, 2012.
•  B. Jiang, T.-Y. Wu, and W.-H. Wong. Approximate bayesian computation with kullback-leibler divergence as data discrepancy. In Amos Storkey and Fernando Perez-Cruz, editors,

Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics

, volume 84 of Proceedings of Machine Learning Research, pages 1711–1721, Playa Blanca, Lanzarote, Canary Islands, 09–11 Apr 2018. PMLR.
•  M. Park, W. Jitkrittum, and D. Sejdinovic. K2-abc: Approximate bayesian computation with kernel embeddings. In Arthur Gretton and Christian C. Robert, editors, Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, volume 51 of Proceedings of Machine Learning Research, pages 398–407, Cadiz, Spain, 09–11 May 2016. PMLR.
•  E. Bernton, P. E. Jacob, M. Gerber, and C. P. Robert. Approximate bayesian computation with the wasserstein distance. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 81(2):235–269, 2019.
•  J. Weed and F. Bach. Sharp asymptotic and finite-sample rates of convergence of empirical measures in wasserstein distance. Bernoulli, 25(4A):2620–2648, 11 2019.
•  J. Rabin, G. Peyré, J. Delon, and M. Bernot. Wasserstein barycenter and its application to texture mixing. In

Scale Space and Variational Methods in Computer Vision

, pages 435–446, Berlin, Heidelberg, 2012. Springer Berlin Heidelberg.
•  N. Bonneel, J. Rabin, G. Peyré, and H. Pfister. Sliced and Radon Wasserstein barycenters of measures. Journal of Mathematical Imaging and Vision, 51(1):22–45, 2015.
•  I. Deshpande, Z. Zhang, and A. G. Schwing. Generative modeling using the sliced Wasserstein distance. In

IEEE Conference on Computer Vision and Pattern Recognition

, pages 3483–3491, 2018.
•  A. Liutkus, U. Şimşekli, S. Majewski, A. Durmus, and F.-R. Stoter. Sliced-Wasserstein flows: Nonparametric generative modeling via optimal transport and diffusions. In International Conference on Machine Learning, 2019.
•  S. Kolouri, P. E. Pope, C. E. Martin, and G. K. Rohde. Sliced Wasserstein auto-encoders. In International Conference on Learning Representations, 2019.
•  N. Bonnotte. Unidimensional and Evolution Methods for Optimal Transportation. PhD thesis, Paris 11, 2013.
•  I. Deshpande, Y.-T. Hu, R. Sun, A. Pyrros, N. Siddiqui, S. Koyejo, Z. Zhao, D. Forsyth, and A. Schwing. Max-Sliced Wasserstein distance and its use for GANs. In IEEE Conference on Computer Vision and Pattern Recognition, 2019.
•  K. Nadjahi, A. Durmus, U. Şimsekli, and R. Badeau. Asymptotic Guarantees for Learning Generative Models with the Sliced-Wasserstein Distance. In Advances in Neural Information Processing Systems, 2019.
•  S. A. Sisson, Y. Fan, and M. A. Beaumont. Overview of Approximate Bayesian Computation. arXiv e-prints, page arXiv:1802.09720, Feb 2018.
•  D. T. Frazier, G. M. Martin, C.P. Robert, and J. Rousseau. Asymptotic properties of approximate Bayesian computation. Biometrika, 105(3):593–607, 06 2018.
•  G. Peyré, M. Cuturi, et al. Computational optimal transport. Foundations and Trends® in Machine Learning, 11(5-6):355–607, 2019.
•  S. T. Rachev and L. Rüschendorf. Mass transportation problems. Vol. I. Probability and its Applications (New York). Springer-Verlag, New York, 1998. Theory.
•  C. Villani. Optimal Transport: Old and New. Grundlehren der mathematischen Wissenschaften. Springer, 2009 edition, September 2008.
•  T. Toni, D. Welch, N. Strelkowa, A. Ipsen, and M. Stumpf. Approximate bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society, Interface / the Royal Society, 6:187–202, 03 2009.
•  E. Klinger, D. Rickert, and J. Hasenauer. pyABC: distributed, likelihood-free inference. Bioinformatics, 34(20):3591–3593, 05 2018.
•  K. Nadjahi, V. De Bortoli, A. Durmus, R. Badeau, and U. Şimşekli. Sw-abc software implementation.
•  A. Buades, B. Coll, and J.-M. Morel. A non-local algorithm for image denoising. In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), volume 2, pages 60–65. IEEE, 2005.
•  L. Fan, F. Zhang, H. Fan, and C. Zhang. Brief review of image denoising techniques. Visual Computing for Industry, Biomedicine, and Art, 2(1):7, 2019.
•  D. Martin, C. Fowlkes, D. Tal, J. Malik, et al. A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. Iccv Vancouver:, 2001.
•  S. Kolouri, K. Nadjahi, U. Simsekli, R. Badeau, and G. K. Rohde. Generalized Sliced Wasserstein Distances. In Advances in Neural Information Processing Systems, 2019.