The Wasserstein Distances Between Pushed-Forward Measures with Applications to Uncertainty Quantification

TL;DR

This paper demonstrates that the Wasserstein metric effectively measures the difference between push-forward measures in uncertainty quantification, providing bounds and insights into how approximations of system responses impact the distribution of outputs.

Contribution

The authors establish bounds on Wasserstein distances between push-forward measures based on function approximations, highlighting the metric's suitability over traditional density-based measures.

Findings

Wasserstein distance bounds for $f ext{-} ho$ and $g ext{-} ho$ are derived.

Lower bounds for $p=1,2$ are provided.

Applications to numerical methods in uncertainty quantification are discussed.

Abstract

In the study of dynamical and physical systems, the input parameters are often uncertain or randomly distributed according to a measure $\varrho$. The system's response $f$ pushes forward $\varrho$ to a new measure $f\circ \varrho$ which we would like to study. However, we might not have access to $f$ but only to its approximation $g$. We thus arrive at a fundamental question -- if $f$ and $g$ are close in $L^q$, does $g\circ \varrho$ approximate $f\circ \varrho$ well, and in what sense? Previously, we demonstrated that the answer to this question might be negative in terms of the $L^p$ distance between probability density functions (PDF). Here we show that the Wasserstein metric is the proper framework for this question. For any $p\geq 1$, we bound the Wasserstein distance $W_p (f\circ \varrho , g\circ \varrho) $ from above by $\|f-g\|_{q}$. Furthermore, we provide lower bounds for the cases of $p=1,2$. Finally, we apply our theory to the analysis of common numerical methods in the field of computational uncertainty quantification.

Equations

$$\Vert f-g_h{\Vert }_{\infty },\Vert \nabla f-\nabla g_h{\Vert }_{\infty }\le K{h}^{\tau },$$

$$\Vert p_{\mu }-p_{\nu }{\Vert }_{L^p}\le \tilde{K}{h}^{\tau },$$

$$f_{k}(\alpha)=\left\{\begin{array}[]{ll}0&x\in[0,\frac{1}{2}]\,,\\ (x-\frac{1}{2})^{k}&x\in[\frac{1}{2},1]\,,\par\end{array}\right.\qquad k\geq 1\,.$$

$$W_p(\mu ,\nu ):={\left[\underset{\gamma \in \Gamma }{\inf }\int |x-y{|}^{p}d\gamma (x,y)\right]}^{\frac{1}{p}},$$

$$\mu (x)=\underset{\mathbb{R}}{\int }\gamma (x,y)dy,\nu (y)=\underset{\mathbb{R}}{\int }\gamma (x,y)dx.$$

$$F_{\sigma }(y):=\sigma ([y,\infty ))=\underset{y}{\overset{\infty }{\int }}{p}_{\sigma }(t)dt,$$

$$W_1(\mu ,\nu )=\underset{\mathbb{R}}{\int }|F_{\mu }(x)-F_{\nu }(x)|dx.$$

$$W_p(\mu ,\nu )\le \Vert f-g{\Vert }_{\infty }.$$

$$W_p(\mu ,\nu )\le \Vert f-g{\Vert }_p.$$

$$W_p(\mu ,\nu )\lesssim \Vert f-g{\Vert }_{\infty }^{\frac{p}{q+p}}\Vert f-g{\Vert }_q^{\frac{q}{q+p}},$$

$$\left|{\mathbb{E}}_{\varrho }f-{\mathbb{E}}_{\varrho }g\right|\le W_1(\mu ,\nu )\le \Vert f-g{\Vert }_{L^1(\Omega ,\varrho )}.$$

$$W_1(\mu ,\nu )=\Vert f-g{\Vert }_{L^1(\Omega ,\varrho )}.$$

$$\Vert \sigma {\Vert }_{{\dot{H}}^{-1}(\mathbb{R})}:=\underset{\Vert q{\Vert }_{{\dot{H}}^1(\mathbb{R})}\le 1}{\sup }|⟨q,\sigma ⟩|,$$

$$\Vert \sigma {\Vert }_{L^2}^2={\int }_{\mathbb{R}}{\left|\hat{\sigma }(\xi )\right|}^{2}d\xi ,\Vert \sigma {\Vert }_{{\dot{H}}^{-1}}^2={\int }_{\mathbb{R}}{\left|\frac{\hat{\sigma }(\xi )}{|\xi |}\right|}^{2}d\xi ,$$

$$\Vert \mu -\nu {\Vert }_{{\dot{H}}^{-1}}\le \max {\left\{\Vert p_{\mu }{\Vert }_{\infty },\Vert p_{\nu }{\Vert }_{\infty }\right\}}^{\frac{1}{2}}{W}_2(\mu ,\nu ),$$

$$A^{\star }=\{\alpha \in \Omega |\varrho (B(0,1))\cdot |\alpha {|}^d\le \varrho (A)\},$$

$$f^{\star }(\alpha )=\underset{0}{\overset{\infty }{\int }}\text{mathbbm}{1}_{\{{\alpha }^{\prime }\in \Omega |f({\alpha }^{\prime })>t{\}}^{\star }}(r)dt,r:=|\alpha |,$$

$$\Vert f{\Vert }_p=\Vert f^{\ast }{\Vert }_p,\Vert \nabla f^{\ast }{\Vert }_p\le \Vert \nabla f{\Vert }_p,$$

$$W_2(\mu ,\nu )\ge A_k\left|{\mathbb{E}}_{\varrho }{f}^k-{\mathbb{E}}_{\varrho }{g}^k\right|,$$

$$A_k=A_k(f,g,r)=\frac{\sqrt{2k-1}}{k}{\left(\max (f,g{)}^{2k-1}-\min (f,g{)}^{2k-1}\right)}^{-\frac{1}{2}}{\tau }^{\frac{1}{2}}\Vert r{\Vert }_{\infty }^{-\frac{1}{2}},$$

$$\mu_{j}:\,=f_{*}\varrho\big{|}_{I_{j}}\,,\qquad\nu_{j}:\,=g_{*}\varrho\big{|}_{I_{j}}\,,$$

$$W_p^p(\mu ,\nu )\le \sum _{j=1}^M{W}_p^p({\mu }_j,{\nu }_j)+W_p^p({\mu }_{M+1},{\nu }_{M+1}),$$

$$W_p^p({\mu }_j,{\nu }_j):$$

$$\le \left(\underset{(x,y)\in f(I_j)\times g(I_j)}{\sup }|x-y{|}^p\right)\varrho (I_j),$$

$$\underset{(x,y)\in f(I_j)\times g(I_j)}{\sup }|x-y{|}^p\le {\left(|f(t)-g(t)|+2ϵ\right)}^p.$$

$$W_p^p(\mu ,\nu )\le \Vert f-g{\Vert }_{\infty }^p\sum _{j=1}^M\varrho (I_j)+o(ϵ)\sum _{j=1}^M\varrho (I_j)+{\left(\Vert f{\Vert }_{\infty }+\Vert g{\Vert }_{\infty }\right)}^p\varrho (I_{M+1}).$$

$$\underset{(x,y)\in f(I_j)\times g(I_j)}{\sup }|x-y{|}^p$$

$$=|f(t_j)-g(t_j){|}^p+o(ϵ).$$

$$W_p^p(\mu ,\nu )\le \sum _{j=1}^M|f(t_j)-g(t_j){|}^p\varrho (I_j)+o(ϵ)\sum _{j=1}^M\varrho (I_j).$$

$$W_p^p(\mu ,\nu )\le W_p^p({\mu }_{{\Omega }_r},{\nu }_{{\Omega }_r})+W_p^p({\mu }_{\Omega \setminus {\Omega }_r},{\nu }_{\Omega \setminus {\Omega }_r}).$$

$$\Vert f-g{\Vert }_{L^q(\Omega )}^q\ge \Vert f-g{\Vert }_{L^q({\Omega }_r)}^q\ge \varrho ({\Omega }_r)\cdot r^q,$$

$$W_p^p(\mu ,\nu )\le \Vert f-g{\Vert }_{\infty }^p\Vert f-g{\Vert }_q^q{r}^{-q}+r^p.$$

Full text

1. Introduction

1.1. Problem formulation

Suppose a domain $\Omega\subseteq\mathbb{R}^{d}$ is equipped with a Borel probability measure $\varrho$ and that a function $f:\Omega\to\mathbb{R}$ pushes forward $\varrho$ to a new measure $\mu:\,=f_{*}\varrho$, i.e., $f_{*}\mu(B)=~{}\varrho(f^{-1}(B))$ for every Borel set $B\subseteq\mathbb{R}$. We wish to characterize $\mu$, but only have access to a function $g$ which approximates $f$. If $\|f-g\|_{L^{q}(\Omega,\varrho)}$ is small, does $\nu:\,=g_{*}\varrho$ approximate $\mu$ well, and if so in what sense?

1.2. Motivation

To motivate this rather abstract question, consider the following toy example: a harmonic oscillator is described by the ordinary differential equation (ODE) $y^{\prime\prime}(t)+y=0$ with $y(0)=0$ and $y^{\prime}(0)=v$. Suppose we are interested in $f(v)=~{}y^{2}(\pi/2;v)$. By solving this ODE, we know that $f(v)=[v\sin(\pi/2)]^{2}=v^{2}$. In many other cases, however, we do not have direct access to $f$, but only to its approximation $g$. This could happen for various reasons – it may be that we can only compute $f(\alpha)$ numerically, or that we approximate $f$ using an asymptotic method. Following on the harmonic oscillator example, suppose we know $f(v)$ only at four given points $v_{1}$, $v_{2}$, $v_{3}$, and $v_{4}$. For any other value of $v$, we approximate $f(v)$ by $g(v)$, which linearly interpolates the adjacent values of $f$, see Fig. 2(a).

The parameters and inputs of physical systems are often noisy or uncertain. We thus assume in the harmonic oscillator example that the initial speed $v$ is drawn uniformly at random from $[1,2]$. In these settings, $f(v)$ is random, and we are interested in the distribution of $f(v)$ over many experiments. Even though $f$ and $g$ look similar in Fig. 2(a), the probability density functions (PDF) of $\mu=f_{*}\varrho$ and $\nu=g_{*}\varrho$, denoted by $p_{\mu}$ and $p_{\nu}$ respectively, are quite different, see Fig. 2(b). We would therefore like to have guarantees that $\nu$ approximates the original measure of interest $\mu$ well.

It might seem obvious that the distance between $f$ and $g$ controls the distance between $\mu$ and $\nu$. This hypothesis fails, however, when one estimates this distance using the PDFs $p_{\mu}$ and $p_{\nu}$. For example, let $f(\alpha)=\alpha$ and $g(\alpha)=\alpha+\delta\sin((10\delta)^{-1}\alpha)$, where $1\gg\delta>0$. Since $\|f-g\|_{\infty}=\delta$, the two functions are seemingly indistinguishable from each other, see Fig. 3(a). Consider the case where $\varrho$ is the Lebesgue measure on $[0,1]$ . Then, since both functions are monotonic, $p_{\mu}(y)=~{}1/f^{\prime}(f^{-1}(y))=~{}1$ and $p_{\nu}(y)=1/g^{\prime}(g^{-1}(y))$, see [14] for details. Hence, $p_{\nu}$ is onto $[1.1^{-1},0.9^{-1}]\approx[0.91,1.11]$ and so $\|p_{\mu}-p_{\nu}\|_{\infty}>0.1$, irrespectively of $\delta$, see Fig. 3(b). The lack of apparent correspondence between $\|f-g\|_{q}$ and $\|p_{\mu}-p_{\nu}\|_{p}$ for any pair of integers $p$ and $q$ suggests that the PDFs are not a well-suited metric for the problem depicted in Fig. 1. Instead, in this paper we propose the Wasserstein distance as the proper framework to measure the distance between $\mu$ and $\nu$.

1.3. Relevant literature

The harmonic oscillator example in Sec. 1.2 serves as a toy example for a broad class of problems. While the ODE $y^{\prime\prime}(t)+y=0$ can be solved explicitly, many other differential equations do not admit such closed-form solutions. Instead, we only have an approximation for the quantities of interest at our disposal. Indeed, the general settings presented above have spurred numerous papers in a field of computational science known as Uncertainty-Quantification (UQ), see e.g., [14, 22, 43, 53, 54, 55].

Perhaps surprisingly, the full approximation of $\mu$ (rather than its moments alone) in these particular settings received little theoretical attention in the literature, even though it is of practical importance in diverse fields such as ocean waves [1], computational fluid dynamics [8], hydrology [9], aeronautics [17], biochemistry [25], and nonlinear optics [30, 36]. Even though $\|f-g\|_{q}$ does not control $\|p_{\mu}-p_{\nu}\|_{p}$ in general (see e.g., Fig. 3), a previous result by Ditkowski, Fibich, and the author gives sufficient conditions for PDF approximation:

Theorem 1** (Ditkowski, Fibich, and Sagiv [14]).**

Let $f\in C^{2}([0,1]^{d})$ and let $g_{h}\in~{}C^{2}([0,1]^{d})$ be an interpolant of $f$ on a tensor grid of maximal spacing $h>0$ such that

[TABLE]

where $K=K(f,d)$ and $\tau>0$ is fixed. Then

[TABLE]

for every $1\leq p<\infty$, with a constant $\tilde{K}=\tilde{K}(f,d,q)$.

The conditions on $g$ are motivated by spline interpolation method, see Sec. 4 for further details. Theorem 1 is, to the best of our knowledge, a first result in the direction of this paper’s main question. Even so, Theorem 1 is limited in several ways

1. (1)

The demand $|\nabla f|\geq\tau_{\rm f}>0$ is an arbitrary condition from an application standpoint. 2. (2)

The differentability and the pointwise derivative-approximation conditions $\|\nabla f-\nabla g\|_{\infty}\lesssim~{}h^{\tau}$ are strong demands which many other approximation methods do not fulfill. 3. (3)

It is essential that the domain $\Omega$ is compact for the proof to hold. 4. (4)

Even when $d=1$, it is required that $d\varrho(\alpha)=c(\alpha)\,d\alpha$ with $c\in C^{1}(\bar{\Omega})$. For comparison, absolute continuity is a weaker condition, as it requires that $c\in C(\Omega)\cap L^{1}(\Omega)$.

The Wasserstein distance (see Sec. 1.4) is thus proposed to measure the distance $\mu$ and $\nu$ since it does not suffer from the drawbacks of the norms $\|p_{\mu}-p_{\nu}\|_{p}$. Admittedly, the $L^{p}$ distances between the PDFs are both natural in practice and are associated with rich statistical theory; for $p=1$, then $\|p_{\mu}-p_{\nu}\|_{1}$ is twice the total variation [13], and $\|p_{\mu}-p_{\nu}\|_{2}^{2}$ is the Integrated Square Error, which is a building block in non-parametric statistics [47]. Nevertheless, the analysis of the norms $\|p_{\mu}-p_{\nu}\|_{p}$ in terms of the functions $f$ and $g$ can be technically cumbersome; if e.g., $\varrho$ is the Lebesgue measure, then $p_{\mu}(y)$ is proportional to $\int_{f^{-1}(y)}1/|\nabla f|\,d\sigma$, where $d\sigma$ is the $(d-1)$ dimensional surface measure [14]. Moreover, the distance $\|p_{\mu}-p_{\nu}\|_{p}$ is difficult to work with since it assumes that $\mu$ and $\nu$ have distributions. This is not always the case. For example, let $\varrho$ be the Lebesgue measure on $[0,1]$ and let

[TABLE]

Although $f_{k}$ is in $C^{k}([0,1])$, the measure $\mu_{k}=(f_{k})_{*}\varrho$ is not a absolutely continuous measure and does not have a PDF since $\mu(\left\{0\right\})=1/2$. It is therefore natural to look for other ways to measure the distance between $\mu$ and $\nu$. There are many ways to define distances between probabilities and measures, such as total variation, mutual information, and Kullback-Leibler divergence. The equivalencies and relationships between these norms, metrics, and semi-metrics are the topics of many studies, see e.g., [19].

1.4. The Wasserstein distance

In order for us to choose the proper metric between $\mu$ and $\nu$, we revisit Fig. 3. While the two PDFs seem very different on a local scale, they are quite similar on a coarser scale. For example, $\mu([0.3,0.4])\approx\nu([0.3,0.4])$ and so, if we were to ask what is the probability that the results of many experiments are between $0.3$ and $0.4$, then both $\mu$ and $\nu$ would have provided similar answers. More loosely speaking, since $p_{\nu}$ is oscillatory, the regions where $p_{\nu}>p_{\mu}$ and the regions where $p_{\nu}<p_{\mu}$ are adjacent, and therefore cancel-out each other. The PDF, on the other hand, is the derivative of the measure, and it is therefore heavily affected by local differences. Another disadvantage of the norm $\|p_{\mu}-p_{\nu}\|_{q}$ is that it does not take geometry into account. Consider for example a family of standard Gaussian measures with mean $t\in\mathbb{R}$, i.e., $p_{\mu,t}(y)=\exp(-(y-t)^{2})/\sqrt{2\pi}$ (see Fig. 4). Then for every $t>2$, $\|p_{\mu,t}-p_{\mu,0}\|_{1}\approx 2$, regardless of whether $t=3$ or $t=10$ or $t=1,000$.

A widely-popular metric that overcomes some of the above issues is the Wasserstein metric. Given two probability measures $\mu$ and $\nu$ on $\mathbb{R}$ with $p\geq 1$ finite moments, the Wasserstein distance of order $p$ is defined as

[TABLE]

If the $p$-th moments of $\mu$ and $\nu$ are finite, then a minimizer exists, $W_{p}(\mu,\nu)$ is finite, and it is a metric [38, 50]. Intuitively, the Wasserstein distance with $p=1$ computes the minimal work (distance times force) by which one can transfer a mound of earth that “looks” like $\mu$ to a one that “looks” like $\nu$, and it is therefore referred to as the earth-mover’s distance.

As noted, some of the difficulties in approximating the PDFs arise from the inverse proportion between $p_{\mu}$ and $p_{\nu}$ and the gradients of $f$ and $g$, respectively. It is therefore natural to avoid these issues by considering the integral of the PDF, the cumulative distribution function (CDF)

[TABLE]

for any Borel measure $\sigma$. Indeed, the Wasserstein distance of order $p=1$ is related to the CDF by the following theorem.

Theorem** (Salvemini [37], Vallender [49]).**

For any two probability Borel measures $\mu$ and $\nu$ on $\mathbb{R}$,

[TABLE]

This theorem reinforces the notion that $W_{1}$ is not as sensitive to local effects as $\|p_{\mu}-p_{\nu}\|_{p}$. Indeed, Fig. 3(c) shows that the two CDFs of $\mu$ and $\nu$ are almost indistinguishable. Furthermore, in the previous Gaussians example (see Fig. 4), $W_{1}(p_{\mu,t},p_{\mu,0})=t$ by direct computation of the CDFs, then, and the same can be proven for $p=2$ as well [20, 28]. Hence, the geometric distance between the Gaussians matters in under the Wasserstein metric. Generally, Wasserstein distances are a central object in optimal transport theory [38, 50], and have also become an increasingly popular in such diverse fields as image processing [29, 35], optimization and neural networks [3], well-posedness proofs for partial differential equations with an associated gradient-flow [7], and numerical methods for conservation laws [41, 45].

1.5. Structure of the paper

The rest of the paper is organized as follows: Sec. 2 presents the main theoretical results of this paper. The upper bounds on $W_{p}$ (Theorems 2 and 3) are presented in Sec. 2.1, and the lower bounds on $W_{1}$ (Corollary 4) and $W_{2}$ (Theorem 5) are presented in Sec. 2.2. The proofs and some technical details of these results are presented in Sec. 3. Finally, in Sec. 4 the theoretical results are applied to the numerical analysis of uncertainty quantification methods, and a numerical example is presented.

2. Main Results

2.1. Upper Bounds

In what follows, $\Omega\subseteq\mathbb{R}^{d}$ is a Borel set, $\varrho$ is a Borel probability measure on $\Omega$, $f,g:\Omega\to\mathbb{R}$ are measurable, $\mu=f_{*}\varrho$, $\nu=g_{*}\varrho$, and $L^{p}=L^{p}(\Omega,\varrho)$ for any $1\leq p\leq\infty$ unless stated otherwise.

Theorem 2**.**

Let $f$ and $g$ be continuous on $\bar{\Omega}$. (i) If $f,g\in L^{\infty}(\Omega,\varrho)$, then for every $p\geq 1$

[TABLE]

(ii) If $\Omega$ is bounded and $f,g\in L^{p}(\Omega,\varrho)$ then

[TABLE]

This result is sharp. Let $\varrho$ be any probability measure on $[0,1]$ and let $f(\alpha)\equiv x_{0}$ and $g(\alpha)\equiv y_{0}$, for some $x_{0},y_{0}\in\mathbb{R}$. Then $\mu$ and $\nu$, are the Dirac delta distribution centered at $x_{0}$ and $y_{0}$, respectively, and the only distribution $\gamma\in\Gamma$ is $\gamma=\delta_{(x_{0},y_{0})}$. Hence, $W_{p}^{p}(\mu,\nu)=|x_{0}-y_{0}|^{p}=\|f-g\|_{\infty}^{p}$. Furthermore, as opposed to Theorem 1, this theorem does not even demand that $f$ and $g$ be differentiable, and puts no restrictions on the Borel measure $\varrho$. Though this theorem is only valid for domains in $\mathbb{R}^{d}$, a generalization of case (i) to (infinite-dimensional) Polish spaces has been achieved by Boussaid [6].

Item (ii) of Theorem 2 uses $L^{p}$ information to bound $W_{p}(\mu,\nu)$. In many cases, however, upper bounds on $f-g$ are known only in a specific $L^{q}$ space. The next theorem shows how $L^{q}$ error estimates can provide nontrivial upper bounds on $W_{p}(\mu,\nu)$ for any $p$, even if $p\neq q$.

Theorem 3**.**

Under the assumptions (i)+(ii) of Theorem 2, then for every $p,q\geq 1$,

[TABLE]

where $\lesssim$ denotes inequality up to a constant which depend only on $p$ and $q$.

This limit agrees with Theorem 2 in the limit $q\to\infty$ and when $q=p$ (up to a constant). Furthermore, for any $q\neq p$, the bound in Theorem 3 may improve the $L^{\infty}$ bound in Theorem 2, since $\varrho$ is a probability measure, $(f-g)\in L^{\infty}\cap L^{q}$, and so $\|f-g\|_{q}\leq\|f-g\|_{\infty}$.

2.2. Lower bounds

The $W_{1}$ lower bound is the direct result of the Monge-Kantorovich duality, see Sec. 3.3 for details and proof.

Corollary 4**.**

If $f,g\in C(\bar{\Omega})$ and $\Omega$ is bounded, then

[TABLE]

Moreover, if $f\geq g$ almost everywhere with respect to $\varrho$, then

[TABLE]

We note that since the upper bound is sharp (see discussion on Theorem 2) and since equality might hold, the lower bound is sharp too. We further note that in the case where $\Omega$ is the unit circle, lower bounds on $W_{1}$ in terms of the Fourier coefficients of $f$ were proved by Steinerberger [42].

Next, to bound $W_{2}(\mu,\nu)$ from below, we introduce two concepts: the Sobolev space $\dot{H}^{-1}$ and the symmetric decreasing rearrangement. For any Borel measure $\sigma$ on $\mathbb{R}$, define the semi-norm

[TABLE]

where $\|q\|_{\dot{H}^{1}}^{2}=\int|q^{\prime}(x)|^{2}\,dx$ [2]. Note that $\|\sigma\|_{\dot{H}^{-1}}<\infty$ only if $\sigma(\mathbb{R})=0$. Another way to understand the Sobolev semi-norm $\dot{H}^{-1}$ and to compare it to the more frequently used $L^{2}$ norm is through Fourier analysis. By Plancharel Theorem

[TABLE]

where $\hat{\sigma}$ is the Fourier transform of $\sigma$ [2]. Thus, if $\mu$ and $\nu$ are different only in high frequencies, then their $L^{2}$ difference might be much higher than their $\dot{H}^{-1}$ difference (due to the $1/|\xi|$ term in the integral). Intuitively, it means that highly local effects in $\sigma=\mu-\nu$ are “subdued” in the negative Sobolev semi-norm. This is analogous to the way local effects in the PDFs are subdued in the $W_{1}$ distance, i.e., in the CDFs (see Fig. 3). As noted, this property also characterizes the Wasserstein distance, and indeed Loeper [27] and Peyre [31] related $W_{2}(\mu,\nu)$ to $\|\mu-\nu\|_{\dot{H}^{-1}}$ in the following theorem:

Theorem** (Loeper [27], Peyre [31]).**

Let $\mu$ and $\nu$ be probability measures on $\mathbb{R}$ with densities $p_{\mu},p_{\nu}\in L^{\infty}(\mathbb{R})$, respectively. Then,

[TABLE]

We now introduce the Symmetric decreasing rearrangement by an absolutely-continuous Borel probability measure on $\Omega\subseteq\mathbb{R}^{d}$ [26]. The symmetric decreasing rearrangement of a measurable set $A$ is

[TABLE]

where $B(0,1)\subset\mathbb{R}^{d}$ is the unit ball around the origin. Next, for a measurable non-negative function $f:\mathbb{R}^{d}\to\mathbb{R}_{+}$, define the symmetric decreasing rearrangement as

[TABLE]

where $\mathbbm{1}_{B}$ is the identifier of a set $B\subseteq\mathbb{R}^{d}$. For a numerical example of the symmetric decreasing rearrangement, see Fig. 5. In more intuitive terms, $f^{*}$ is the unique monotonic decreasing function such that $\varrho(A(f,t))=\varrho(A(f^{*},t))$ for all $t\in\mathbb{R}$, where $A(f,t):\,=~{}\{\alpha~{}{\rm s.t.}~{}f(\alpha)\geq~{}t\}$ are the super-level sets of $f$. Moreover, since $f^{*}$ is monotonic decreasing, one also have that $A(f^{*},t)$ is the interval $[0,\varrho(A(f,t))]$. The symmetric decreasing rearrangement is an important object in real analysis [26], with notable properties such as the Pólya-Szego inequality [32]

[TABLE]

for all $p\geq 1$. Hence, the symmetric decreasing rearrangement $f^{*}$ minimizes all Sobolev $W^{1,p}$ norms among the class of functions with the same super-level sets, it can be said to be the “canonical” representative this class.

Theorem 5**.**

Let $I$ be a closed and bounded interval equipped with an absolutely-continuous probability measure $\varrho$ with a bounded and continuous weight function $r(\alpha)$, i.e., $d\varrho(\alpha)=r(\alpha)d\alpha$, and let $f,g\in C^{1}$ with $|(f^{*})^{\prime}|,|(g^{*})^{\prime}|>\tau>0$. Then, for every $k\in\mathbb{N}$

[TABLE]

where $A_{k}$ is a positive coefficient given by

[TABLE]

and the $\max$ and $\min$ are taken over all $x\in I$.

We remark that even though $A_{k}=A_{k}(f,g,r)$ depends on $f$ and $g$, it does not depend directly on $f-g$. Hence, for a sequence $(g_{n}(\alpha))_{n=1}^{\infty}$ which converges uniformly to $f$, for each $k\in\mathbb{N}$, then $A_{k}(f,g_{n},r)$ would converge to a positive constant as $n\to\infty$. A specific example to the computation of the coefficients $A_{k}$ can be found in Sec. 3.4.

3. Proofs of main results and technical discussion

3.1. Proof of Theorem 2

Proof.

We begin with the case where $f$ and $g$ are uniformly continuous in $\Omega$. Let $\epsilon>0$, then by uniform continuity there exists $\eta=\eta(\epsilon)>0$ such that $|f(\alpha)-f(\beta)|<\epsilon$ and $|g(\alpha)-g(\beta)|<\epsilon$ for every $\alpha,\beta\in\Omega$ such that $|\alpha-\beta|<\eta$. Let $L\in\mathbb{N}$ and partition $[-L,L]^{d}$ to $M$ equal-size boxes $\{\tilde{I}_{j}\}_{j=1}^{M}$ such that ${\rm diam}(\tilde{I}_{j})<\eta$. Let $I_{j}=\tilde{I}_{j}\cap\Omega$ for every $1\leq j\leq M$ and let $I_{M+1}:\,=\Omega\setminus[-L,L]^{d}$. Next, let

[TABLE]

i.e., the measures induced by $f(I_{j})$ and $g(I_{j})$ for every $1\leq j\leq M+1$. Since $\int_{\mathbb{R}}\mu_{j}=\int_{\mathbb{R}}\nu_{j}=\varrho(I_{j})$, we can transport $\mu$ to $\nu$ by transporting each $\mu_{j}$ to $\nu_{j}$. Even though this might not be the optimal transport between $\mu$ and $\nu$, since $W_{p}$ is defined as an infimum over all transports then

[TABLE]

where $\Gamma_{j}$ is the set of all measures whose marginals are $\mu_{j}$ and $\nu_{j}$. For $1\leq j\leq M$, since ${\rm diam}(I_{j})<\eta$, then by uniform continuity for any $t\in I_{j}$

[TABLE]

Here the proofs of the $L^{p}$ and $L^{\infty}$ bounds slightly diverge and we begin with proving that $W_{p}(\mu,\nu)\leq\|f-g\|_{\infty}$. For any $1\leq j\leq M$ then $(|f(t)-g(t)|+2\epsilon)^{p}\leq\left(\|f-g\|_{\infty}+2\epsilon\right)^{p}$ . Similarly, for $j=M+1$, the supremum in (2) is bounded from above by $\left(\|f\|_{\infty}+\|g\|_{\infty}\right)^{p}$. Combining these bounds together, we have that

[TABLE]

Since $\varrho$ is a probability measure $\sum_{j=1}^{M}\varrho(I_{j})=\varrho(\Omega)=1$ and as $L\to\infty$ the third term on the right-hand-side vanishes. Hence $W_{p}^{p}(\mu,\nu)\leq\|f-g\|_{\infty}^{p}+2o(\epsilon)$ for every $\epsilon>0$, and so $W_{p}\leq\|f-g\|_{\infty}$.

Next, consider the case where $f,g$ are continuous on $\Omega$, but not uniformly continuous. Then for any two sequences $\epsilon_{n}\to 0$ and $L_{n}\to\infty$, choose $\eta_{n}=\eta(\epsilon_{n},L_{n})$ which satisfies the uniform continuity condition on the compact domain $\bar{\Omega}\cap[-L_{n},L_{n}]$. Then, by partitioning this domain to sufficiently many boxes $M_{n}=M(\eta_{n})=M(\epsilon_{n},L_{n})$ such that ${\rm diam}(I_{j,n})\leq\eta_{n}$, the proof holds as $n\to\infty$.

Finally, we prove that $W_{p}(\mu,\nu)\leq\|f-g\|_{p}$. Here we require that $\Omega$ is bounded, and so we can choose $L$ such that $\Omega\subseteq~{}[-L,L]^{d}$. For $1\leq j\leq M$ we have that for some $t_{j}\in I_{j}$ then

[TABLE]

Substituting this inequality in (2) yields

[TABLE]

As the partition is refined (i.e., $M\to\infty$ and $\eta\to 0$), the first element on the right-hand-side converges to $\|f-g\|_{L^{p}(\varrho)}$. Since $\varrho$ is a probability measure, $\sum_{j=1}^{M}\varrho(I_{j})=1$, and so the second element on the right-hand-side is $o(\epsilon)$. Since this inequality is true for any $\epsilon>0$, the proof follows.

∎

3.2. Proof of Theorem 3

Proof.

Define $\Omega_{r}:\,=\{\alpha\in\Omega~{}|~{}|f(\alpha)-g(\alpha)|\geq r\}$ for any $r>0$, and let $\mu_{\Omega_{r}}$, $\mu_{\Omega\setminus\Omega_{r}}$, $\nu_{\Omega_{r}}$, and $\nu_{\Omega\setminus\Omega_{r}}$ be the measure induced by $f(\Omega_{r})$, $f(\Omega\setminus\Omega_{r})$, $g(\Omega_{r})$, and $g(\Omega\setminus\Omega_{r})$, respectively. For any $p\geq 1$,

[TABLE]

The fist term on the right-hand-side of (3) is bounded from above by $\|f-g\|_{\infty}^{p}\varrho(\Omega_{r})$, due to Theorem 2. To bound $\varrho(\Omega_{r})$, note that

[TABLE]

where the first inequality is due to monotonicity of $\varrho$, and the last inequality is due the continouity of $|f-g|^{q}$. Hence, $\varrho(\Omega_{r})\leq\|f-g\|_{q}^{q}r^{-q}$, and so the first term in the right-hand-side of (3) is bounded from above by $\|f-g\|_{\infty}^{p}\|f-g\|_{q}^{q}r^{-q}$. Since the $L^{\infty}$ upper bound of Theorem 2 is applicable to $f$ and $g$, and since $\varrho\left(\Omega\setminus\Omega_{r}\right)\leq 1$, then the second term on the right-hand-side of (3) is bounded from above by $\|f-g\|_{L^{\infty}(\Omega\setminus\Omega_{r})}^{p}\leq r^{p}$. Having bounded from above both terms on the right-hand-side of (3), then

[TABLE]

To minimize the right-hand-side of this inequality, we derive with respect to $r$ and get that the minimum is achieved at $r_{\min}=(qp^{-1}\|f-g\|_{q}^{q}\cdot\|f-g\|_{\infty}^{p})^{1/(p+q)}$, and so

[TABLE]

∎

3.3. Proof of Corollary 4

Proof.

The Monge-Kantorovich duality states that [50]

[TABLE]

where $L(w)$ is the Lipschitz constant of $w$. So, to prove a non-trivial lower bound for $\mu=f_{*}\varrho$ and $\nu=g_{*}\varrho$, it is sufficient to provide any function $w$ for which the integral is not zero. Let $w(y)=y$. Since $L(w)=1$, then $W_{1}(\mu,\nu)\geq|\int_{\mathbb{R}}y\,d\mu(y)-\int_{\mathbb{R}}y\,d\nu(y)|$, which, by change of variables, means that $W_{1}(\mu,\nu)\geq|\int_{\Omega}(f(\alpha)-g(\alpha))\,d\varrho(\alpha)|$. Combined with Theorem 2 we arrive at the corollary. ∎

3.4. Proof of Theorem 5

Proof.

By definition of the symmetric decreasing rearrangement, $\mu=f^{*}_{*}\varrho$ and $\nu=g^{*}_{*}\varrho$. Moreover, since the theorem requires that $|(f^{*})^{\prime}|,|(g^{*})^{\prime}|>\tau:\,>0$, we can assume without loss of generality that $f$ and $g$ are strongly monotonically decreasing. Next, we have the following standard lemma (for proof, see e.g., [14]):

Lemma**.**

Let $h\in C^{1}(I)$ be piecewise monotonic, let $d\varrho(\alpha)=r(\alpha)d\alpha$ where $r$ is continuous in $\Omega$. Then the PDF of the measure $\sigma=h_{*}\varrho$ is given by

[TABLE]

Hence, by definition and the above lemma

[TABLE]

Consider the first integral under the supremum. By change of variables $y=f(x)$, we have

[TABLE]

Doing the respective change of variable for the second integral under the supremum, we have

[TABLE]

For ease of notations, denote $M=\max_{x\in I}\left\{f(x),g(x)\right\}$ and $m=\min_{x\in I}\left\{f(x),g(x)\right\}$. Since $f$ and $g$ are continuous on a closed bounded interval, both $m$ and $M$ are finite. Fix $k\in\mathbb{N}$, and let $q_{k}(x)=c_{k}x^{k}$, where the normalization constant $c_{k}:\,=(\sqrt{2k-1}/k)(M^{2k-1}-m^{2k-1})^{-1/2}$ is chosen so that $\|q_{k}\|_{\dot{H}^{1}[m,M]}=1$.111It might seem that the choice of the interval $[m,M]$ is made ad-hoc. However, this proof can be carried out in the space $\dot{H}^{-1}(\mathbb{R})$ regardless, by the following construction: extend $q_{k}(y)$ to $\mathbb{R}$ by setting $q_{k}(y)=q_{k}(m)$ for $y<m$ and $q_{k}(y)=q_{k}(M)$ for $y>M$. Since outside $[m,M]$, $q_{k}^{\prime}\equiv 0$, then $\|q_{k}\|_{\dot{H}^{1}(\mathbb{R})}=\|q_{k}\|_{\dot{H}^{1}([m,M])}$, and $\langle q_{k},\mu-\nu\rangle$ is unchanged too since $\mu-\nu$ is supported only on $[m,M]$. Our choice is also consistent with the result by Peyre [31], since these also ”take place” on the supports of $\mu$ and $\sigma$. Hence, substituting $q_{k}$ in (4) for every $k\in\mathbb{N}$

[TABLE]

Finally, to bound $W_{2}$ from below we need Loeper and Peyre’s theorem, and so we need to compute $\|p_{\mu}\|_{\infty}$ and $\|p_{\nu}\|_{\infty}$. As noted, since $f=f^{*}$ is strictly decreasing, it is also continuously differentiable almost everywhere. Hence. by the result noted above, $p_{\mu}=r(f^{-1}(y))/|f^{\prime}(f^{-1}(y))|$ almost everywhere, and so $\|p_{\mu}\|_{\infty}\leq\tau^{-1}\|r\|_{\infty}$. Since the same holds for $g$ and $\nu$ as well, we substitute in Loeper’s and Peyre’s bound and get that

[TABLE]

∎

We complement the proof by an example of a direct computation of the coefficients $A_{k}$. Let $f(\alpha)=3\alpha-3$, $g(\alpha)=2\alpha-2$ and $\varrho$ is the Lebesgue measure on $[0,1]$, then by direct computation we have that $M=0$, $m=-3$, $\|r\|_{\infty}=1$, $\tau=2$, and so

[TABLE]

4. Convergence of uncertainty-quantification methods and numerical examples

We apply the main theoretical results of this paper to the analysis of uncertainty quantification (UQ) methods. In many applications, one can only compute the quantity of interest $f(\alpha)$ for a finite subset of $\alpha$ values $\{\alpha_{j}\}_{j=1}^{N}$. To compute $\mu=f_{*}\varrho$, we first use these sampled values $\{f(\alpha_{j})\}_{j=1}^{N}$ to construct an approximate function $g(\alpha)$, and then we approximate $\mu\approx\nu=g_{*}\varrho$, see Fig. 1. This measure-approximation problem is characterized by the following trade-off: The computational cost comes from direct computation of the samples $\{f(\alpha_{j})\}_{j=1}^{N}$, 222Since $g$ is given in closed form, e.g., by a polynomial, it is computationally cheap to estimate the measure $\nu=g_{*}\varrho$. Computing $f(\alpha_{j})$, on the other hand, might involve a full numerical solution of a PDE. and so it increases linearly with $N$. On the other hand, we expect the approximation error to decrease with the sample size $N$, i.e., as we improve the sampling resolution. The question is, therefore, how to construct $g$ such that $\mu$ is accurately approximated with a small sample size $N$.

In terms of numerical analysis, the main result of this paper is that upper bounds on $\|f-g\|_{q}$ do guarantee an upper bound on the Wasserstein distances $W_{p}(\mu,\nu)$. This in turn immediately implies an upper bound on the $L^{1}$ distance between the CDFs, due to the previously-noted Salvemini-Vallender identity $W_{1}(\mu,\nu)=\|F_{\mu}-F_{\nu}\|_{1}$ [49].

The upper bounds on the Wasserstein-error stand in sharp contrast to the $L^{q}$ errors between the PDFs, since in general an upper bound on $\|f-g\|_{q}$ does not guarantee an upper bound on $\|p_{\mu}-p_{\nu}\|_{L^{p}}$, for any finite $p$ and $q$ [14]. We therefore see that the way we define the approximation-error in this problem is not a mere technicality, but rather determines the results of the convergence analysis. Furthermore, we see that CDFs are “easier” to approximate than PDFs, in the sense that the it is easier to guarantee their efficient approximation.

We demonstrate the applicability of our theory for two approximation methods (surrogate models), spline interpolation and generalized Polynomial Chaos (gPC).

4.1. Spline interpolation.

Given an interval $\Omega=[\alpha_{\min},\alpha_{\max}]$ and grid-points $\alpha_{\min}=\alpha_{1}<\alpha_{2}<\cdots<\alpha_{N}=\alpha_{\max}$, an interpolating $m$-th order spline $g(\alpha)\in C^{m-1}(\Omega)$ is a piecewise polynomial of order $m$ that interpolates $f(\alpha)$ at the grid-points, endowed with some additional boundary conditions so that it is unique. See [12, 33] for comprehensive expositions on splines, see [34, 39] for their extension to multidimensional domains via tensor-products, and see [4, 22] for their applicability to UQ problems. Since Theorem 6 is directly applicable to spline interpolation [14], if $g$ is the spline interpolant of $f$, then the PDFs of $\mu$ and $\nu$ are close, i.e., $\|p_{\mu}-p_{\nu}\|_{L^{p}}$ is bounded from above for any $1\leq p<\infty$. We show that in these settings, the Wasserstein distance between the measures is also bounded from above.

Theorem 6**.**

Let $f\in C^{m+1}([0,1]^{d})$, let $g(\alpha)$ be its (tensor-product) spline interpolant of order $m$ on a (tensor-product) grid of maximal grid size $h$, and let $\varrho$ be a probability Borel measure. Then, for every $p\geq 1$,

[TABLE]

where $N$ is the total number of interpolation points, and where $\lesssim$ and $\approx$ denote inequality and equality up to constants independent of $h$ and $N$, respectively.

See Sec. LABEL:sec:spline_pf for the proof. Theorem 6 is stronger than Theorem 1 in three aspects. First, Theorem 6 holds for a broader function class than the application of Theorem 1 to splines, since it does not require that $|\nabla f|>\tau_{f}>0$, or even that the underlying measure $\varrho$ would be absolutely continuous. Second, Theorem 6 is non-trivial even for those functions for which Theorem 1 does apply. To obtain a “trivial” upper bound, note that for any two probability measures of $\mu$ and $\nu$ with PDFs $p_{\mu}$ and $p_{\nu}$, then

[TABLE]

where ${\rm diam}(\Omega)$ is the diameter of $\Omega={\rm supp}(\mu)\cup{\rm supp}(\nu)$ [19]. Since $f$ and $g$ are continuous on a compact set, they are bounded, and so the supports of $\mu$ and $\nu$ are bounded as well. Hence, ${\rm diam}(\Omega)<\infty$, and so by Theorem 1, $W_{1}(\mu,\nu)\leq Kh^{m}$. Theorem 6, however, guarantees an additional order of accuracy and so non-trivially improves the previous results.333Unfortunately, Theorem 6 cannot improve the $L^{1}$ bound in Theorem 1 since, in general, $\|p_{\mu}-p_{\nu}\|_{1}\lesssim W_{1}(\mu,\nu)$ only for finite spaces [19]. Finally, Theorem 6 applies not only for $p=1$ but for all $p\geq 1$.

Numerical example.

Let

[TABLE]

We use a cubic spline interpolant on a grid of $N$ uniformly-spaced points, with the not-a-knot boundary condition [12]. Theorem 6 guarantees that in this case $W_{p}(\mu,\nu)\lesssim N^{-4}$. Indeed, Fig. 6(a) shows the $W_{1}$ difference between the two measures as a function of $N$, and that the convergence rate is $N^{-4.59}$.

4.2. Generalized Polynomial Chaos (gPC)

Next, we turn to study $W_{p}$ convergence of $L^{2}$-spectral methods, for which PDF convergence is an open problem. We focus on the widely popular generalized Polynomial Chaos (gPC).

Review of the Collocation gPC method. For a more detailed exposition, see e.g., [21, 54]. Let the Jacobi polyomials $\{p_{n}(x)\}_{n=0}^{\infty}$ be the orthogonal polynomials with respect to $\varrho$, i.e., $p_{n}$ is a polynomial of degree $n$, and $\int_{-1}^{1}p_{n}(\alpha)p_{m}(\alpha)\,d\varrho(\alpha)=\delta_{n,m}$, see [44] for details. This family of orthogonal polynomials constitutes an orthonormal basis of the space $L^{2}(\Omega,\varrho)$, i.e., for every $f\in L^{2}$ one can expand

[TABLE]

This expansion converges spectrally, i.e., if $f$ is in $C^{r}$, then $\{\hat{f}(n)\}\lesssim n^{-r}$, and if $f$ is analytic in an ellipse $E\subseteq\mathbb{C}$ that contains $[-1,1]$, then $|\hat{f}(n)|\lesssim e^{-\gamma n}$, for some $\gamma>0$. Thus, one has that for such analytic functions

[TABLE]

The expansion coefficients $\{\hat{f}(n)\}$ can be approximated using the Gauss quadrature

[TABLE]

where $\{\alpha_{j}\}_{j=1}^{N}$ are the quadrature points, the distinct and real roots of $p_{N}(\alpha)$, and $w_{j}$ are the quadrature weights [11]. We define the gPC Collocation approximation $g_{N}$ to be the truncated expansion of $f$ with the quadrature-based coefficients $\hat{f}_{N}(n)$. We remark that this approximation method has a much simpler form – The gPC collocation approximation is also the unique interpolating polynomial of $f$ of order $N-1$ at the quadrature points [14]. We remark that our theory can also be applied to Galerkin-gPC methods [53].

Density estimation in UQ: The main appeal of the gPC method is its spectral $L^{2}$ convergence. As noted above, it is an open question whether this can be used to prove convergence of the PDFs, i.e., an upper bound on $p_{\mu}-p_{\nu}$ in some $L^{p}$. However, Theorem 3 implies that spectral $L^{2}$ convergence of $g_{N}$ to $f$ can yield fast convergence of $W_{p}(\mu,\nu)$ for any $1\leq p<\infty$.

Theorem 7**.**

Let $f$ be analytic in an ellipse in the complex plane that contains $[-1,1]$, and let $d\varrho(\alpha)=k(1-~{}\alpha)^{\beta_{1}}(1+~{}\alpha)^{\beta_{2}}d\alpha$, for any $\beta_{1},\beta_{2}\in\mathbb{R}$ and a proper normalization constant $k=k(\beta_{1},\beta_{2})$. Let $g(\alpha)$ be the collocation gPC approximation of $f$, i.e., the $N$-th order polynomial interpolant of $f$ at the respective Gauss quadrature points. Then, for every $p\geq 1$,

[TABLE]

where $\gamma$ does not depend on $N$.

Proof.

If $f$ is analytic, the truncated expansion has the exponential accuracy

[TABLE]

for some constant $\gamma>0$ [46, 51, 53]. Next, since the collocation gPC is a spectrally accurate approximation of the polynomial projection in $L^{2}$ [21], then $\|f-g\|_{2}\lesssim e^{-\gamma N}$ as well for $N\gg 1$. Finally, since $\|f-\pi_{N}(f)\|_{\infty}$ does not grow exponentially [21], Theorem 3 applies. ∎

Two particularly important cases of this theorem are when $\varrho$ is the Lebesgue measure, associated with the Legendre polynomials ($\beta_{1,2}=0$) and the measure associated with the Chebyshev polynomials ($\beta_{1,2}=-1/2$). By Theorem 7, the convergence of the Wasserstein metric stands in sharp contrast to that of the PDFs, i.e., of $\|p_{\mu}-p_{\nu}\|_{L^{q}}$. As previously noted, the convergence of the PDFs for the gPC method has not been proved, and might not be obtained at all for moderate values of $N$ [14]. It remains an open question whether Theorem 7 can be extended to measures with an unbounded support, such as the normal and the exponential distributions. Such a generalization might require a generalization of Theorem 3 to unbounded domains. We further note that Theorem 7 can be extended to measures $\varrho^{\prime}$ that are bounded from above by $\varrho$, see [15] for details.

Numerical example.

We approximate the same function $f$, as defined in (5), and approximate it using polynomial interpolation at Gauss-Legendre quadrature points (see Sec. LABEL:sec:gpc_pf). Since $f$ is analytic, Theorem 7 guarantees that the gPC-based $\nu$ converges exponentially in $N$ to that of $\mu$, see Fig. 6(b). The convergence of the respective PDFs, on the other hand, is polynomial at best (see Fig. 6(c)). Quantitatively, the $W_{1}$ error decreases by $8$ orders of magnitudes between $N=4$ and $N=120$, whereas the $L^{1}$ between the PDFs decreases by only $4$.

4.3. Comparison to the histogram method

This paper, as noted, is motivated by the following class of algorithms: to approximately characterize $\mu=f_{*}\varrho$, first approximate $f$ by $g$, and then approximate $\mu$ by $\nu=g_{*}\varrho$. How does this approach compare with more standard statistical methods? We focus on one common nonparametric statistical density estimation method, the histogram method; Given i.i.d. samples from $\mu$, denoted by $f(\alpha_{1})=y_{1},\ldots,f(\alpha_{N})=y_{N}$, and a partition of the range of $f(\boldsymbol{\alpha})$ into $L$ disjoint intervals (bins) $\{B_{\ell}\}_{\ell=1}^{L}$, the histogram estimator of the PDF is

[TABLE]

where $\mathbbm{1}_{B_{\ell}}$ is the characteristic function of bin $B_{\ell}$ [52]. The histogram methods is intuitive and easy to implement. What is then the advantage of approximation-based UQ methods? In Sec. 4.4, using results by Bobkov and Ledoux [5], we prove that

Corollary 8**.**

Under the conditions of Theorem 6, the $d$-dimensional, $m$-th order spline-based estimator of $\mu$ outperforms the histogram method on average in the $W_{p}$ sense when $d<2(m+1)$.

The average in this corollary refers to all i.i.d. realizations of $y_{1},\ldots,y_{N}$ from $\mu$. This corollary is an example of the so-called “curse of dimensionality”. To maintain a constant resolution and accuracy, the amount of data points (and hence the computational complexity) needs to increase exponentially with the dimension. Hence, above a certain dimension, it is preferable to ignore the underlying structure (i.e., the approximation of $f$ by $g$) and to consider only the empirical distribution of the i.i.d. samples $\{f(\alpha_{j})\}_{j=1}^{N}$.

Proof.

The error of spline interpolation is controlled by the following theorem

Theorem** (de Boor [12] and Hall and Meyer [23]).**

Let $f\in C^{m+1}\left(\left[\alpha_{\min},\alpha_{\max}\right]\right)$, and let $g(\alpha)$ be its ”not-a-knot”, clamped or natural $m$-th spline interpolant. Then

[TABLE]

where $C_{\rm spl}^{(j)}>0$ is a universal constant that depends only on the type of boundary condition and $j$, $m$, and $h=\max_{1<j\leq N}\lvert\alpha_{j}-\alpha_{j-1}\rvert$.

This result is extended for higher dimensions using the the construction of tensor-product grid and tensor-product splines. The definitions here become more technical, and we refer to Schultz [39] for further detail. We note that even in the multidimensional case, the error is still bounded by the spacing $h^{m+1-j}$. However, the number of grid points $N$ is proportional to $h^{-d}$ (this is the so-called curse of dimensionality which we previously mentioned). By the above error bounds, and by Theorem 2, we have that $W_{p}(\mu,\nu)\leq\|f-g\|_{\infty}\lesssim h^{m+1}$. ∎

4.4. Proof of Corollary 8

Proof.

Given $N$ i.i.d. from $\mu$, denoted by $y_{1},\ldots,y_{N}$, define the empirical distribution as

[TABLE]

where $\delta_{y}$ is the Dirac delta distribution centered at the point $y\in\mathbb{R}$. Under certain broad assumptions (see [5] for details), $\mathbb{E}W_{p}(\mu,\mu_{\rm emp})\lesssim N^{-1/2}$, where the expectancy in these bounds is over all realizations of $y_{1},\ldots,y_{N}$ with respect to the measure $\mu$ [5].

By the triangle inequality and linearity of expectation,

[TABLE]

where $d\mu_{\rm hist}(y)=p_{\rm hist}(y)\,dy$ is the measure defined by the histogram estimator. It is therefore sufficient to show that $\mathbb{E}W_{p}(\mu_{\rm emp},\mu_{\rm hist})\lesssim N^{-(1+1/p)}$ for any $p\geq 1$. We will prove a slightly stronger claim – that $W_{p}(\mu_{\rm emp},\mu_{\rm hist})\lesssim N^{-(1+1/p)}$ for every set of numbers $y_{1},\ldots,y_{N}$.

Let $\{B_{\ell}\}_{\ell=1}^{L}$ be the bins of the histogram estimator and let $\mu_{{\rm emp},\ell}$ and $\mu_{{\rm hist},\ell}$ be the restriction of the measures $\mu_{\rm emp}$ and $\mu_{\rm hist}$ to $B_{\ell}$, respectively, for every $1\leq\ell\leq L$. By definition, there are exactly $N\cdot\mu_{{\rm hist}}(B_{\ell})$ samples that fall into $B_{\ell}$, and so $\mu_{{\rm hist},\ell}(B_{\ell})=\mu_{{\rm emp},\ell}(B_{\ell})$. Hence, the two measures $\mu_{{\rm emp},\ell}$ and $\mu_{{\rm hist},\ell}$ are comparable in the Wasserstein metric and we can write that

[TABLE]

Since $\mu_{{\rm hist},\ell}$ is uniform on $B_{\ell}$ for any $\ell$, the Wasserstein distance is the greatest if all of the samples in $B_{\ell}$ are located on the extreme edge of the bin, i.e., if $y_{j}\in B_{\ell}$ then $y_{j}=a_{\ell}$, where we denote $B_{\ell}=[a_{\ell},b_{\ell}]$. Hence, for every $1\leq\ell\leq L$,

[TABLE]

and so

[TABLE]

where the second inequality is due to the partition, in which $(b_{\ell}-a_{\ell})\sim N^{-1}$, and the last equality holds since $\mu_{\rm emp}$ is a probability measure and since $\{B_{\ell}\}_{\ell=1}^{L}$ is a partition of its support. ∎

5. Acknowledgments

The author would like to thank S. Steinerberger for many useful comments and advice. This research was partially carried out during a stay of the author as a guest of R.R. Coifman and the Department of Mathematics at Yale University, whose hospitality is gratefully acknowledged.