Due to its desirable convexity and convergence features, Wasserstein
distance is considered as an important measure of dissimilarity
in many applications, ranging from computer vision and machine learning to Seismic and
Bayesian inversion;
see e.g. [1]}, [2]}, [3]}, [4]}.
Unfortunately, the numerical computation of Wasserstein metric is in
general infeasible in high dimensions and is mainly limited to low dimensions; see e.g. [1]}.
For instance, in the case \(c({\bf x}, {\bf y}) = || {\bf x}- {\bf y} ||^2\) , the quadratic Wasserstein distance between two
continuous probability density functions \(f: {\mathcal {X}}\rightarrow {\mathbb {R}}+\)  and \(g: {\mathcal {Y}}\rightarrow {\mathbb {R}}+\)  is given by [6]}
\(W_2^2(f, g) = \int _{\mathcal {X}} | {\bf x} - M({\bf x}) |^2 \, f(\bf x)\, d{\bf x},\) 
