Spectral convergence of probability densities for forward problems in uncertainty quantification

The remainder of the paper is organized as follows: Sec. 2 provides the preliminaries on the gPC method, pushforward measures and a more detailed account of previous results. Secs. 3 and 4 detail the paper’s main theoretical results and their proofs, respectively. We then turn in Sec. 5 to provide an alternative, optimal-transport based analysis of the measure pushforward problem. Finally, in Sec. 6 we present numerical experiments and outline open problems.

We present a brief review of the gPC method, see [50] for a thorough introduction. For clarity, we first review the gPC method for the one-dimensional case. The Legendre polynomials $\{p_{n}(\alpha)\}_{n=0}^{\infty}$ are the set of orthogonal polynomials with respect to the Lebesgue measure on $\Omega=~{}[-1,1]$, see [38]

The Legendre polynomials constitute an orthonormal basis of $L^{2}(\Omega,d\alpha)$, i.e.,

The Galerkin-gPC expansion of a function $f\in L^{2}(\Omega)$ of order $n$ is the projection of $f$ to its first $n+1$ modes, i.e., $g_{n}(\alpha)=\sum_{j=0}^{n}\hat{f}(n)p_{n}(\alpha)$. In the multidimensional case $\Omega=[-1,1]^{d}$, the expansion of order $n$ of $f$ is defined as the $L^{2}(\Omega)$ projection of $f$ into the space of polynomials of maximal degree $n$ [8, 21, 52]

where $p_{\bf j}$ is the tensor-product Legendre polynomial of multi-index ${\bf j}=(j_{1},\ldots,j_{d})$ and $\|{\bf j}\|_{\infty}=\max_{i}|j_{i}|$. If $f$ can only be evaluated at a discrete set of points, the integrals in expansion coefficients $\hat{f}(n)$ cannot be exactly computed. The collocation-gPC approach is to approximate these coefficients using the Gauss quadrature formula: Let $q:\Omega\to\mathbb{R}$ be an integrable function, then the Gauss-Legendre quadrature formula of order $n$ is $\int_{-1}^{1}q(\alpha)\,d\alpha\approx\sum_{k=1}^{n}g(\alpha_{k})w_{k}$, where $\{\alpha_{k}\}_{k=1}^{N}$ are the distinct and real roots of $p_{N}(\alpha)$, $w_{k}:=\int_{\Omega}l_{k}(\alpha)\,d\mu(\alpha)$ are the weights, and $l_{k}(\alpha)$ are the Lagrange interpolation polynomials with respect at the quadrature points [14]. Taking $q=p_{j}(\alpha)f(\alpha)$ yields the gPC collocation method in one dimension [51],

The generalization of (3) to multiple dimensions is a direct result of integration by tensor-grid Gauss quadratures, see [8, 14, 51] for details. We further note that the gPC-collocation polynomial (3) is also the polynomial interpolant of degree $N-1$ in the Guass-Legendre points, i.e., $g_{n}(\alpha_{\bf k})=f(\alpha_{\bf k})$ for all quadrature points [12].

Recall that for $A=[-1,1]^{d}$ or $\mathbb{R}^{d}$, and any pair of integers $k,p\geq 1$, the corresponding Sobolev space is defined as [2, 18]

where ${\bf j}$ is a multi-index and $D^{\bf j}u=(\Pi_{k=1}^{d}\partial_{\alpha_{k}}^{j_{k}})u(\alpha)$. For $p=2$ we will adopt the conventional notation $W^{k,2}(A)=H^{k}(A)$. Both the Galerkin (2) and the collocation expansion (3) converge spectrally in $L^{2}(\Omega)$, where the underlying measure is the Lebsegue measure [24, 40, 47, 50]. By this we mean that if $f\in H^{k}(\Omega)$ for some $k\geq 0$, then [8, 24, 50]

and if $f$ is analytic (in the multivariate sense), then $\|f-g_{n}\|_{L^{2}(\Omega)}$ decays exponentially in $n$ [13], See [53] and the references therein for details. From a probability and UQ point of view, the convergence in $L^{2}(\Omega)$ guarantees that the moments of $g_{n}$ converge to the moments of $f$. Indeed, if $\varrho$ is the Lebesgue measure on $\Omega$, a simple application of the Cauchy-Schwartz inequality shows that

and therefore for e.g., a smooth quantity of interest $f$, the gPC method provides an accurate approximation of $\mathbb{E}f$ for a relatively low order $n$. The convergence of other moments follow similarly, see e.g., [15].¹¹1For measures $\varrho$ which are not the Lebesgue measure, see [16]. Both gPC expansions (2) and (3) also converge in higher-regularity Sobolev spaces, as given by the now classical result:

If $f$ is analytic (for $d=1$) on an ellipse in $\mathbb{C}$ with foci at $\pm 1$ and with both axis summing to $r>1$, then in one dimension $\hat{f}(n)=O(n^{-1/2}r^{-n})$ [46, 47]. By Parseval identity (with respect to Legendre polynomials),

and so the $L^{2}$ approximation error is almost exponential as well. Furthermore, by uniform convergence of $g_{n}$ to $f$ on $\Omega$, one can differentiate $g_{n}$ term-wise, yielding faster than polynomial (exponential for all practical purposes) convergence of $g_{n}$ to $f$ in all Sobolev spaces $H^{k}(\Omega)$ with $r\geq 0$.

The pushforward of a Borel measure $\varrho$ on $\Omega$ by a measurable $f:\Omega\to\mathbb{R}$ is a measure $\mu:\,=f_{*}\varrho$ defined by $\mu(A)=\varrho(f^{-1}(A))$ for any Borel subset $A\subseteq\mathbb{R}$. If a measure $\mu$ on $\mathbb{R}$ is absolutely continuous with respect to the Lebesgue measure, its probability density function (PDF) is its Radon-Nykodim derivative, i.e., $p_{\mu}\in L^{1}(\mathbb{R})$ which satisfies $\mu(A)=\int_{A}p_{\mu}(y)\,dy$ for any Borel set $A\subseteq\mathbb{R}$. Alternatively, if the cumulative distribution function (CDF) $F_{\mu}(y)=\mu(f^{-1}(-\infty,y))$ is differentiable, then $p_{\mu}(y)=dF_{\mu}(y)/dy$. In the one-dimensional case, if $f$ is piecewise $C^{1}$ and piecewise monotonic, and $\varrho$ is an absolutely continuous probability measure, then [15]

This relation is the source of many of the difficulties and peculiarities in understanding PDFs of pushforward measures. For example, if $f(J)=c$ for some constant $c$ on an interval $J$ with $\varrho(J)>0$, then $\mu$ has a singular part at $c$ and therefore has no PDF. Even for a non-constant smooth monotonic function such as $f(\alpha)=\alpha^{2}$ and a simple $\varrho$ such as the uniform measure on $[0,1]$, then $p_{\mu}(y)\sim 1/\sqrt{y}$ is singular. Analogously for $d>1$, then $p_{\mu}\sim\int_{f^{-1}(y)}p_{\varrho}|\nabla f|^{-1}d\sigma$, where $d\sigma$ is the $(d-1)$-dimensional surface elements.

The practical goal of surrogate models in the context of PDF approximation is to approximate $\mu$ by $\nu_{n}=(g_{n})_{\sharp}\varrho$ with small error terms $\|p_{\mu}-p_{\nu_{n}}\|_{L^{q}(\mathbb{R})}$ for some $q\geq 1$, while maintaining $n$ small, as $n$ is a good proxy to the computational cost. To see that, first note that since sampling arbitrarily many times from the polynomial $g_{n}$ is computationally cheap, approximating $p_{\nu_{n}}$ to arbitrary precision given $g_{n}$ is relatively cheap as well, see the analysis in [15]. Therefore, it is constructing $g_{n}$ which is costly. In the collocation gPC (3), the larger the degree of approximation $g=g_{n}$ is, the more evaluations of $f$ are required. If for example $f$ models the response of a partial differential equation (PDE), each such evaluation amounts to a solution of the PDE, which is usually computationally expansive. In the Galerkin-type gPC (2), one usually projects the original PDE with random parameters to several PDEs with deterministic parameters, the number of which also grows with the degree $n$, see [50] for details. Therefore, in either case guaranteeing good accuracy for a small degree $n$ should be the goal of our analysis.

As noted in the introduction, $L^{2}(\Omega)$ convergence alone does not guarantee convergence of the pushforward PDFs. Indeed, one can construct a sequence such that $g_{n}\to f$ in $L^{2}(\Omega)$ but $\|p_{\nu_{n}}-p_{\mu}\|_{L^{q}(\mathbb{R})}>{\rm const}(q)$ for all $n,q\geq 1$ [15]. Previously however [32], we proved that $L^{2}$ convergence is sufficient to establish convergence in the weaker Wasserstein metric, or more precisely, that $\rm Wass_{2}(f_{\sharp}\varrho,g_{\sharp}\varrho)\leq\|f-g\|_{L^{2}(\Omega,\varrho)}$ (see relevant definitions in Sec. 5). Furthermore, we showed that uniform boundedness of $\|f-g_{n}\|_{\infty}$ combined with $L^{2}(\Omega)$ convergence guarantees Wasserstein-$p$ convergence for any $1\leq p<\infty$, see (26) below. Since $\rm Wass_{1}(\mu,\nu)$ is equal to the $L^{1}(\mathbb{R})$ convergence of the CDFs [34, 43], the convergence of the CDFs in $L^{1}$ is a corollary of our result. A result in the direction of PDF approximation by surrogate methods was obtain by Ditkowski, Fibich, and the author in [15]. If $g_{n}$ is taken to be the spline interpolant of $f$ of order $m$, on e.g., a uniform grid with a total of $N$ grid points, then $\|p_{\mu}-p_{\nu_{n}}\|_{L^{q}(\mathbb{R})}\lesssim N^{-m/d}$ for any $1\leq q<\infty$. The main characteristic of splines that is useful to establish this result is the pointwise $C^{1}(\Omega)$ convergence of $g_{n}$ to $f$, see Theorem 4 below. We will use this tool tool later to prove Theorem 1.

Let $\Omega=\Omega(d)=[-1,1]^{d}$ be equipped with an absolutely continuous probability measure $\varrho$. Consider a smooth function of interest $f:\Omega\to\mathbb{R}$ and let $g_{n}:\Omega\to\mathbb{R}$ be its generalized polynomial chaos (gPC) approximation - either its $L^{2}(\Omega)$ projection to the space of Legendre polynomials of order $\leq n$ (Galerkin type) or its polynomial interpolant on the Gauss-Legendre quadrature points of order $n$ (collocation type). Consider the pushforward measures $\mu:\,=f_{*}\varrho$ and $\nu_{n}:\,=(g_{n})_{*}\varrho$ and denote their respective probability density functions (PDF) by $p_{\mu}$ and $p_{\nu_{n}}$. In what follows, we see under what conditions $p_{\nu_{n}}$ converges to $p_{\mu}$ in $L^{q}(\mathbb{R})$ for $1\leq q<\infty$, and find the convergence rates.

Local $C^{1}$ convergence as established in Lemma 3 implies convergence of PDFs by the following result:

Indeed, (13a) is guaranteed explicitly and (9) guarantees (13b) with an explicit convergence rate $\tau=\sigma-\sigma_{\min}-2>0$.

For brevity, we omit the $n$ subscripts, denoting $g_{n}=g$ and $\nu=\nu_{n}$. First, we treat the case where $|f^{\prime}|>\kappa_{f}>0$. Suppose $f$ is monotonic increasing, i.e., $f^{\prime}>\kappa_{f}>0$. By the embedding theorem (10), that $f\in H^{6}(\Omega)$ implies that $f\in C^{2}(\Omega)$ and so

for every $y\in{\rm range}(f)$ [15, Lemma 4.1]. Since $g$ is a polynomial, $g\in C^{1}(\Omega)$. Furthermore, by Lemma 3 we have that $\|f-g\|_{C^{1}}\lesssim n^{-5/2}\|f\|_{H^{4}(\Omega)}$, and so for sufficiently large $n$, $g^{\prime}(\alpha)>\kappa_{f}/2>0$ for all $\alpha\in\Omega$. Therefore $p_{\nu}(y)=r(g^{-1}(y))/g^{\prime}(g^{-1}(y))$ for every $y\in{\rm range}(g)$. It might be that the ranges of $f$ and $g$, which are the supports of $p_{\mu}$ and $p_{\nu}$, respectively, do not overlap. Assume for simplicity that $f(-1)=g(-1)$ but $g(1)>f(1)$, the other cases can be treated similarly. Then

We begin with the second integral -

where the last inequality is due to the Sobolev-Morrey embedding (10). Therefore, we need only to consider the first integral in (15), and we therefore assume without loss of generality that ${\rm range}(f)={\rm range}(g)$, and so

Denote $\alpha=f^{-1}(y)$ and $\alpha_{*}:\,=\alpha_{*}(\alpha)=g^{-1}(f(\alpha))$, then by change of variables

Since $g^{\prime}(\alpha)$ and $r(\alpha)$ are differentiable, then for any $\alpha\in\Omega$

where $D:\,=\left[\|r\|_{\infty}\|g^{\prime\prime}\|_{\infty}+\|r^{\prime}\|_{\infty}\|f\|_{\infty}\right]$. In general, $D=D_{n}$ as $g=P_{n}f$ depends on $n$, and $\|g^{\prime\prime}\|_{\infty}$ might not be bounded. However, as in Lemma 3, $\|g^{\prime\prime}\|_{\infty}$ and therefore $D$ are uniformly bounded for all $n\geq 1$. Since $g^{\prime}\geq\kappa_{f}/2$ for sufficiently large $n$, then substituting (17) in the integral yields

Since $r(\alpha),(f^{\prime}-g^{\prime})\in L^{2}(\Omega)$, then by the Cauchy-Schwartz inequality ${\rm I}$ is bounded from above by

To bound ${\rm II}$ in (20) from above, we first note that by Lagrange’s mean-value theorem, there exists $\beta$ between $\alpha$ and $\alpha_{*}$ such that $g^{\prime}(\beta)(\alpha-\alpha_{*})=g(\alpha)-g(\alpha_{*})=g(\alpha)-f(\alpha)$, and therefore $|\alpha-\alpha_{*}|\leq|g(\alpha)-f(\alpha)|/\kappa_{\rm f}$. From here, the process of bounding ${\rm II}$ from above is the same as bounding ${\rm I}$, which yields ${\rm II}\leq D/\kappa_{f}\|f-g\|_{L^{2}(\Omega)}$. Therefore $\|p_{\mu}-p_{\nu}\|_{L^{1}(\mathbb{R})}\lesssim\|f-g\|_{H^{1}(\Omega)}$. Applying the relevant Sobolev approximation theorem (4), settles the case where $|f^{\prime}|>\kappa_{f}>0$.

We now turn to the case where $g^{\prime}$ and $f^{\prime}$ vanish at finitely many points. As we show, it does not matter whether these zero points coincide or not, and we will therefore treat the case where $f^{\prime}(-1)=0$, $g^{\prime}(-1)>0$ and $f^{\prime}(\alpha),g^{\prime}(\alpha)>0$ for all $\alpha\in(-1,1]$. Assume without loss of generality that $f(-1)=g(-1)$.³³3We showed how to treat the case where the ranges of $f$ and $g$ do not overlap above. Fix $\varepsilon>0$ and divide the integral for $\|p_{\mu}-p_{\nu}\|_{1}$ on the left hand side of (16) into two domains - (1) isolating the singular point in the PDF, i.e., $y\in[g(-1),g(-1+\varepsilon)]$ and (2) the rest of the domain, $y\in[g(-1+\varepsilon),g(1)]$.

1. (1)

   On the first domain $[g(-1),g(-1+\varepsilon)]$, we take a crude estimation

   (21a)		$$\int\limits_{g(-1)}^{g(-1+\varepsilon)}|p_{\mu}(y)-p_{\nu}(y)|\,dy\leq\int\limits_{g(-1)}^{g(-1+\varepsilon)}p_{\mu}(y)+p_{\nu}(y)\,dy\,.$$
   Taking, for example $p_{\mu}(y)$, by a change of variables $f(\alpha)=y$ we get
   (21b)		$$\int\limits_{g(-1)}^{g(-1+\varepsilon)}p_{\mu}(y)=\int\limits_{g(-1)}^{g(-1+\varepsilon)}\frac{r(f^{-1}(y))}{f^{\prime}(f^{-1}(y))}\,dy=\int\limits_{-1}^{-1+\varepsilon}r(\alpha)\,d\alpha\leq\|r\|_{\infty}\varepsilon\,.$$

2. (2)

   The integral on $[g(-1+\varepsilon),g(1)]$ reads the same as (20), where $\kappa_{f}$ is replaced by $\kappa_{f}^{\varepsilon}:\,=\min\limits_{x\in[-1+\varepsilon,1]}|f^{\prime}|$, yielding

   (22)

   $$\int\limits_{g(-1+\varepsilon)}^{g(1)}|p_{\mu}(y)-p_{\nu}(y)|\,dy\leq\frac{\tilde{D}}{(\kappa_{f}^{\varepsilon})^{2}}\|f-g\|_{H^{1}}\,,$$

   where $\hat{D}$ is $\varepsilon$-independent. How does $\kappa_{f}^{\varepsilon}$ depend on $\varepsilon$? Since in the worst case $f(-1)=f^{\prime}(-1)=\cdots f^{(k)}(-1)=0$ but $f^{(k+1)}(-1)=0$, then by Taylor expansion, $\kappa_{f}^{\varepsilon}\approx\varepsilon^{k}$ for sufficiently small $\varepsilon$.