Polynomial chaos expansions for dependent random variables

TL;DR

This paper introduces a novel Gram-Schmidt orthogonalization method for constructing polynomial chaos expansions with dependent variables, significantly improving accuracy over existing mapping and dominating support approaches.

Contribution

The paper proposes a new GSO-based approach for PCE with dependent variables, enhancing accuracy and robustness in uncertainty quantification.

Findings

GSO-based PCE outperforms existing methods in accuracy

Weighted Leja sequences improve polynomial interpolation stability

Proposed method is effective for non-intrusive stochastic collocation

Abstract

Polynomial chaos expansions (PCE) are well-suited to quantifying uncertainty in models parameterized by independent random variables. The assumption of independence leads to simple strategies for evaluating PCE coefficients. In contrast, the application of PCE to models of dependent variables is much more challenging. Three approaches can be used. The first approach uses mapping methods where measure transformations, such as the Nataf and Rosenblatt transformation, can be used to map dependent random variables to independent ones; however we show that this can significantly degrade performance since the Jacobian of the map must be approximated. A second strategy is the class of dominating support methods which build PCE using independent random variables whose distributional support dominates the support of the true dependent joint density; we provide evidence that this approach appears to produce approximations with suboptimal accuracy. A third approach, the novel method proposed here, uses Gram-Schmidt orthogonalization (GSO) to numerically compute orthonormal polynomials for the dependent random variables. This approach has been used successfully when solving differential equations using the intrusive stochastic Galerkin method, and in this paper we use GSO to build PCE using a non-intrusive stochastic collocation method. The stochastic collocation method treats the model as a black box and builds approximations of model output from a set of samples. Building PCE from samples can introduce ill-conditioning which does not plague stochastic Galerkin methods. To mitigate this ill-conditioning we generate weighted Leja sequences, which are nested sample sets, to build accurate polynomial interpolants. We show that our proposed approach produces PCE which are orders of magnitude more accurate than PCE constructed using mapping or dominating support methods.

Equations

$$f(z)=\sum _{n=1}^{\infty }{\alpha }_n{\varphi }_n(z),$$

$$({\varphi }_i(z),{\varphi }_j(z){)}_{L_{\omega }^2(\Omega )}:={\int }_{\Omega }{\varphi }_i(z){\varphi }_j(z)\omega (z)={\delta }_{i,j},$$

$$\Omega$$

$${\varphi }_{\lambda }(z)=\prod _{i=1}^d{\varphi }_{{\lambda }_i}^i(z_i),$$

$${\int }_{{\Omega }_i}{\varphi }_j^i(z_i){\varphi }_k^i(z_i){\omega }_i(z_i)\mathrm{d}{z}_i$$

$$f(z)$$

$$\Lambda$$

$${\pi }_{\Lambda }$$

$$\sum _{\lambda \in \Lambda }{\alpha }_{\lambda }{\varphi }_{\lambda }(z)=\sum _{n=1}^N{\alpha }_n{\varphi }_n(z).$$

$$\displaystyle\pi_{\Lambda}\coloneqq\mathrm{span}\left\{\prod_{i=1}^{d}z_{i}^{\lambda_{i}}\;\;\big{|}\;\;\lambda\in\Lambda\right\}.$$

$$f({\mathcal{T}}^{-1}(u))=:g(u)$$

$$u_1$$

$$U_1-F_1(z_1)$$

$${\eta }_{R^V}(u)=\frac{1}{\sqrt{(2\pi {)}^d\det (R^V)}}\exp \left(-\frac{1}{2}{u}^T{R^V}^{-1}u\right),$$

$$\omega (z)=\frac{{\eta }_{R^V}(\hat{u})}{{\prod }_{i=1}^d\eta ({\hat{u}}_i)}\prod _{i=1}^d{\omega }_i(z_i).$$

$$R_{ij}^Z={\int }_{\mathbb{R}}{\int }_{\mathbb{R}}\left(\frac{z_i-\mathbb{E}(Z_i)}{\sqrt{\text{Var}(Z_i)}}\right)\left(\frac{z_j-\mathbb{E}(Z_j)}{\sqrt{\text{Var}(Z_j)}}\right)\frac{{\eta }_{R^{V_i,V_j}}({{\Phi }_g}^{-1}(F_i(z_i)),{{\Phi }_g}^{-1}(F_j(z_j)))}{\eta ({{\Phi }_g}^{-1}(F_i(z_i)))\eta ({{\Phi }_g}^{-1}(F_j(z_j)))}{\omega }_i(z_i){\omega }_j(z_j)\mathrm{d}{z}_i\mathrm{d}{z}_j,$$

$$\Vert f{\Vert }_{L_g^p}$$

$${ϵ}_N$$

$$\Vert f-f_N{\Vert }_{L_{\omega }^p}$$

$${\varphi }_n(Z)$$

$${\int }_{\Omega }f(z)\omega (z)\mathrm{d}z\approx \sum _{j=1}^J{w}_{j}f(z^{(j)}),$$

$$(\Psi ({\mathcal{Z}}_J,N){)}_{j,n}={\Psi }_{j,n}$$

$$\forall p\in \mathrm{span}\{{\psi }_n{\}}_{n=1}^N\backslash \{0\},\sum _{j=1}^J{w}_j{p}^2(z^{(j)})>0.$$

$$\sqrt{W}\Psi =QR,$$

$${\varphi }_n(z)=\sum _{k=1}^n{\psi }_k(z)(R^{-1}{)}_{kn},n\in [N]$$

$$({\psi }_n,{\varphi }_k{)}_{L_{\omega }^2}$$

$$\Phi \alpha =y,⟺V\Phi \alpha =Vy,$$

$$z^{(N+1)}=\underset{z\in \Omega }{\mathrm{argmax}}v(z)\prod _{n=1}^N|z-z^{(n)}|$$

$$z^{(N+1)}=\underset{z\in \Omega }{\mathrm{argmax}}|\det v(z)\Phi (\mathcal{Z},z^{(N+1)})|.$$

$$\mathrm{span}\{{\varphi }_n{\}}_{n=1}^N$$

Full text

1 Introduction

All models of modern scientific applications are subject to various sources of uncertainty. The effect of these uncertainties on model predictions can be assessed by viewing the model as an input-output map, where the inputs are random variables with known (but possibly complicated) distribution, and the output is a set of model quantities of interest (QoI). Uncertainty quantification (UQ) then refers to the process of computing the output statistics that result from this input-output map.

UQ of high-fidelity models typically requires large numbers of simulations. Building an approximation or surrogate for the model input-output map is an effective and popular approach to reduce the computational burden of UQ, and numerous techniques have been developed to build such approximations. Within the computational science and engineering community, some of the most widely adopted methods for approximating models parameterized by random variables are those based on generalized polynomial chaos expansions [25, 64], sparse grid approximation [63, 46], Gaussian process models [50] and low-rank tensor decompositions [48]. These methods can be very efficient when building approximations of models subject to independent random variables, and as such, most computational tools are built assuming the input variables are independent. However, there is a dearth of algorithmic options when the variables are dependent.

In this paper we present sampling strategies for building polynomial chaos expansions (PCE) of models influenced by a high-dimensional random vector that are accurate and stable when the components of the random vector are not independent. PCE essentially seeks to build a polynomial approximation of a function (model) whose inputs are the random variables; our goal will be the generation of a PCE surrogate that is accurate in a norm weighted by the probability density function of the dependent variables. PCE are ideally suited to approximating functions of random variables because they employ basis functions which are orthonormal to the probability measure $\omega$ of the variables, and this relationship can be exploited to construct stable approximation schemes that are accurate in regions of high probability.

The stochastic Galerkin [25] and stochastic collocation [2, 63] methods are the two main approaches for generating a PCE surrogate, which amounts to computing a set of PCE coefficients. In this paper we focus on stochastic collocation because it allows the computational model to be treated as a black box. Stochastic collocation proceeds in two steps: (i) running the computational model with a set of realizations of the random parameters and (ii) constructing an approximation of the corresponding model output. In the relevant situation when the model is expensive, step (i) is the most time-consuming portion. Stochastic collocation methods include pseudospectral projection [12, 14], sparse grid interpolation [4, 8], least orthogonal interpolation [43], least squares [40, 59], compressive sensing [17, 65, 29] and low-rank tensor decompositions [10, 18, 26].

Sparse grid interpolations and pseudospectral approximations are by construction ideally suited to approximation when the input random variables are independent. Polynomial-based sparse grid methods utilize univariate high-order (e.g., Gaussian) quadrature rules that are optimal for tensor product approximation [46, 41]. Recently, sub-sampled tensor-product quadrature has been used to generate sample sets for independent random variables in the context of interpolation [35], least squares [54] and sparse regression [58]. There has also been extensive work on using random sampling for least squares [40, 44] and compressive sensing [51, 27, 31] as well as some work on generating deterministic sequences for interpolation [43, 42].

The aforementioned methods usually rely on tensor product structure and are most effective for independent random variables. If one is only interested in computing moments, recent advances in polynomial quadrature for dependent measures can be utilized [1, 15, 32, 16]. The polynomial accuracy of these methods is inadequate for efficient pseudo spectral projection. Consequently regression-based methods are the only viable alternative.

If $f$ is the model under consideration and $Z$ is a vector of random parameters that are input to the model, a PCE approach posits the representation

[TABLE]

where $z$ are realizations of the variables $Z$, $\alpha_{n}$ are the PCE coefficients that must be computed, and the basis functions $\phi_{n}$ are polynomial basis functions that are pairwise orthonormal under an inner product weighted by the probability density of $Z$. Above we assume that $f$ is scalar-valued, but the procedures we describe carry over to vector- or function-valued outputs. When the components of $Z$ are independent, one can generate the multivariate polynomials $\phi_{n}$ from univariate orthogonal polynomials, but such a construction is not easily accomplished when $Z$ has dependent components. One approach for dealing with dependent variables is to build an approximation for a set of independent variables whose tensor-product measure dominates the dependent measure [28, 9]. Such an approach introduces an error which is proportional to the “distance” between the tensor-product and dependent measures.

In this paper we use Gram-Schmidt orthogonalization to generate a set of polynomials orthonormal to arbitrary probability measures. This approach was first proposed in the multivariate setting in [62] for stochastic Galerkin methods and used for solving time-dependent PCE using stochastic Galerkin projection [24]. However, the use of Gram-Schmidt orthogonalization has received comparatively less attention for stochastic collocation. We will observe in this paper a well-known phenomenon, that the use of a Gram-Schmidt procedure to build basis functions is poorly conditioned. However we propose a method that can reduce the amount of ill-conditioning by using specialized preconditioning and sampling approaches. Specifically, in this paper we adapt and improve the strategies for building weighted Leja sequences for interpolation developed for independent random variables in [42].

This paper is devoted to computational studies of novel and recently-developed algorithms for computing PCE expansions for dependent variables. The novelty and main outcomes of this study are described below:

• •

We provide a systematic study, on practical problems, which compares the use of mapping methods (such as the Rosenblatt transformation), measure domination methods, and Gram-Schmidt Orthogonalization (GSO) methods. In short, we observe that GSO provides the most accurate procedure, but may be ill-conditioned in some cases. However, this ill-conditioning frequently does not degrade the accuracy of the resulting emulator.

• •

We propose the use of weighted Leja sequences (where the weight is related to the dependent density) when constructing PCE approximations for dependent measures. An ingredient in this approach is the use of GSO for constructing an orthonormal basis for the dependent measure. We observe that this approach frequently performs much better than the alternatives described above.

In the remainder of this paper we introduce polynomial chaos expansions and their construction for both tensor-product and dependent probability measures. We then propose stable and accurate sampling schemes for regression and interpolation, using PCE, which can be used for arbitrary measures. We then provide an investigation of the Gram-Schmidt procedure, compare its performance with probabilistic transformation maps for multivariate approximation, and conclude with a number of numerical examples to highlight the strengths of our proposed approach.

2 Polynomial chaos expansions for independent random variables

Let $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a function (model) of a $d$-variate random variable $Z=(Z_{1},\ldots,Z_{d})$. The random variable has associated probability density function $\omega(z)$ for $z\in\Omega\subset\mathbb{R}^{d}$. Polynomial chaos expansions represent the model output $f(z)$ as an expansion in orthonormal polynomials, as in (1). The basis functions $\phi_{n}$ are typically constructed to be orthonormal with respect to the density $\omega$, that is

[TABLE]

where $\Omega$ is the support of the density $\omega$ and $\delta_{i,j}$ is the Kronecker delta function. Under mild conditions on the distribution $\omega(z)$, any function $f(z)$ with finite variance, i.e. $f\in L^{2}_{\omega}(\Omega)$ can be represented by a PCE that converges in $L^{2}_{\omega}$ to the true function asymptotically [20]. In this section we describe the “canonical” construction when $Z$ has independent components, which is well-known and standard. The next section discusses how one can construct a PCE basis when $Z$ has dependent components.

Polynomial chaos expansions are most easily constructed when the components of $Z$ are independent. Under the assumption of independence, we have

[TABLE]

where $\omega_{i}$ are the marginal densities of the variables $Z_{i}$, which completely characterizes the distribution of $Z$. This allows us to express the basis functions $\phi$ as tensor products of univariate orthonormal polynomials. That is

[TABLE]

where $\lambda=(\lambda_{1}\ldots,\lambda_{d})\in\mathbb{N}_{0}^{d}$ is a multi-index, and the univariate basis functions $\phi^{i}_{j}$ are defined uniquely (up to a sign) for each $i=1,\ldots,d$, as

[TABLE]

In practice the PCE (1) must be truncated to some finite number of terms, say $N$, defined by a multi-index set $\Lambda\subset\mathbb{N}_{0}^{d}$:

[TABLE]

Frequently the PCE is truncated to retain only the multivariate polynomials whose associated multi-indices have norm at most $p$, i.e.,

[TABLE]

Taking $q=1$ results in a total-degree space having dimension $\text{card}\;\Lambda^{d}_{p,1}\equiv N={d+p\choose d}$. The choice of $\Lambda$ identifies a subspace in which $f_{N}$ has membership:

[TABLE]

Under an appropriate ordering of multi-indices, the expression (3), and the expression (1) truncated to the first $N$ terms, are identical. Defining $[N]:=\{1,\ldots,N\}$, for $N\in\mathbb{N}$, we will in the following frequently make use of a linear ordering of the PCE basis, $\phi_{k}$ for $k\in[N]$ from (1), instead of the multi-index ordering of the PCE basis $\phi_{\lambda}$ for $\lambda\in\Lambda$ from (3). Therefore,

[TABLE]

Any bijective map between $\Lambda$ and $[N]$ will serve to define this linear ordering, and the particular choice of this map is not relevant in our discussion.

3 PCE for dependent random variables

In this section we present three of the most popular strategies for constructing a PCE surrogate when the distribution of $Z$ is not of tensor-product form. The first strategy is a mapping method, where the dependent coordinates are mapped to a set of independent coordinates. The second algorithm identifies a dominating measure of tensor-product form, and a PCE is built with respect to the dominating measure. The final strategy attempts to explicitly compute an orthonormal polynomial basis in $L^{2}_{\omega}$ and to subsequently construct a PCE approximation.

Much of our notation from the previous section carries over, except that we require more caution when speaking of multi-index sets $\Lambda$. Without a tensor-product construction, the definition of the polynomial $\phi_{\lambda}$ is not unique, and thus the meaning of the subspace $\pi_{\Lambda}$ is unclear. However, the following definition mitigates part of this issue:

Definition 3.1**.**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be a downward-closed index set, i.e., $\lambda\in\Lambda$ implies that $\nu\in\Lambda$ for all $\nu\in\mathbb{N}_{0}^{d}$ satisfying $\nu\leq\lambda$. Here, $\leq$ is the partial lexicographic ordering on the $d$-dimensional integer lattice. Then we define

[TABLE]

This definition is consistent with the meaning of $\pi_{\Lambda}$ for independent variables introduced in Section 2 (assuming $\Lambda$ is downward-closed).

3.1 Mapping methods

Let $\mathcal{T}:\Omega\rightarrow{\tilde{\Omega}}$ denote an invertible transformation which maps possibly dependent random variables $Z$ to independent random variables $\mathcal{T}(Z)=U=(U_{1},\ldots,U_{d})$ with marginal densities $\rho_{i}(u_{i})$, $i\in[d]$. Typically, mapping methods are defined so that $U$ has a uniform distribution on ${\tilde{\Omega}}=[0,1]^{d}$ (and hence is a tensor-product distribution). In this case, we can form the approximation $f_{N}$ by first forming an approximation in $U$ space [19, 60]:

[TABLE]

Here, $\psi$ is a tensor product basis of univariate polynomials orthonormal with respect to the univariate densities $\rho_{i}$, in this case the univariate orthonormal Legendre polynomials on $[0,1]$.

This approach is straightforward to implement if $\mathcal{T}$ is available, but the construction of $\mathcal{T}$ is often the bottleneck in implementations. In the following we describe two popular approaches for computing the transformation $\mathcal{T}$.

3.1.1 Rosenblatt Transformation

When the joint cumulative distribution function $F_{Z}$ of $Z$ is continuous, an explicit construction of $\mathcal{T}$ is given by the Rosenblatt transformation [53]. The Rosenblatt transformation defines the following components for the transformation $\mathcal{T}$:

[TABLE]

where $F_{i|i-1,\dots,1}(\cdot\mid Z_{1}=z_{1},\ldots,Z_{i-1}=z_{i-1})$ is the conditional distribution of $Z_{i}$ given that $(Z_{1},\ldots,Z_{i-1})=(z_{1},\ldots,z_{i-1})$. The marginalization needed to compute the conditional distributions generally requires multivariate integration. For example, computing $u_{1}$ requires integration over $d-1$ variables, and hence a quadrature approach suffers the curse of dimensionality.

The inverse Rosenblatt transformation can be used to obtain the original sample $Z$ from the sample $U$ by solving the following optimization problem

[TABLE]

For a fixed order of the one-dimensional transformations, there exists a unique solution to this problem, since the $F_{i}$ are strictly monotonic increasing functions. We can solve this problem by applying standard root-finding algorithms, such as the bisection method.

The Rosenblatt transformation essentially requires a closed form for the joint density $\omega(Z)$. This requirement can be relaxed by utilizing density estimation techniques, such as Gaussian kernel or sparse grid density estimation to construct an approximation to the joint density [22], which can help mitigate the curse of dimensionality.

3.1.2 Nataf Transformation

When the dependence between the random variables $Z$ is linear then the Rosenblatt transformation can be simplified. The resulting Nataf transformation [38] requires the availability of the correlation matrix $R^{Z}\in\mathbb{R}^{d\times d}$, of $Z$, the marginal distributions $F_{i}(z_{i})$, and the marginal density of $Z_{i}$, denoted $\omega_{i}(z_{i})$.

The Nataf transformation assumes that the joint density $\omega(Z)$ can be expressed in terms of the multivariate Gaussian density

[TABLE]

with correlation matrix $R^{V}$ that must be computed. The non-tensorial density $\omega$ is assumed to have the form

[TABLE]

Here $\omega_{i}$ are the marginal distributions of the joint density $\omega$ and $\hat{u}_{i}={\Phi_{g}}^{-1}(F_{i}(z_{i}))$, with $\eta$ and ${\Phi_{g}}$ respectively denoting the univariate standard normal density and cumulative distribution functions. Each entry $R^{Z}_{ij}$ of the correlation matrix $R^{Z}$ is related to $R^{V}_{ij}$ by

[TABLE]

where $R^{V_{i},V_{j}}=\begin{pmatrix}1&R^{V}_{ij}\\ R^{V}_{ij}&1\;\end{pmatrix}\in\mathbb{R}^{2\times 2}.$ The integral in (10) is two-dimensional and so can be solved with Gaussian quadrature and standard root finding algorithms, such as bisection [36].

Applying the Nataf transformation consists of three steps. First the “corrected” correlation matrix ${R^{V}}$ must be constructed by solving (10). This matrix is independent of the marginals of the variable $Z$ and so can be computed once and stored for any subsequent use. The second step creates a set of intermediate correlated Gaussian random variables $\hat{U}=(\hat{U}_{1},\ldots,\hat{U}_{d})$ created from $Z$ via $\hat{U}_{i}={\Phi_{g}}^{-1}(F_{i}(Z_{i}))$. Finally these intermediate variables are decorrelated using $U=L^{-1}\hat{U}$, where $L$ is the lower triangular Cholesky factor of ${R^{V}}$. In this paper we use this Nataf trafansformation $\mathcal{T}_{\text{nataf}}^{\text{gauss}}$ to build PCE using univariate Hermite polynomials. We will also use the transformation $\mathcal{T}_{\text{nataf}}^{\text{unif}}=\mathcal{T}_{n2u}\circ\mathcal{T}_{\text{nataf}}^{\text{gauss}},$ where $\mathcal{T}_{n2u}$ maps independent normally-distributed random variables to independent uniform variables on $[-1,1]$. When using this transformation, we will build PCE based upon tensor-products of univariate Legendre polynomials.

There is a close connection between the Nataf transformation and copula dependence modeling. Sklar’s theorem states that any joint distribution can be expressed as the product of univariate marginal distributions and a copula that describes the dependency structure. The Nataf density in (9) does exactly this, using the Gaussian copula $\eta_{R^{V}}(V)/\prod_{i=1}^{d}\eta_{i}(V_{i})$. We will utilize this connection between the Nataf transformation and Gaussian copula modeling to efficiently generate dependent multivariate samples to test the performance of our proposed method in Section 5. However we will show in the same section that the usage of Nataf transformations for constructing approximations via (5) produces poor approximations when compared to the other approaches considered in this paper.

3.2 Domination methods

A domination method can be used to compute a PCE when $Z$ has dependent components. This approach uses PCE’s consisting of tensor product basis functions to approximate functions of dependent random variables, but this approximation comes at a cost [64, 28, 9]. Given $\Omega$ and $\omega$, the basic idea is to identify a tensor-product density $g$ with support $G\subset\mathbb{R}^{d}$ so that $G\supseteq\Omega$. For example, if $\Omega$ is compact then one may identify $G$ as the smallest bounding box for $\Omega$, and define $g$ as the uniform probability density over $G$. One then constructs a PCE using the density $g$ with the strategy from Section 2. The error committed by this strategy is essentially well-understood.

Denote the $g$-weighted $L^{p}$ norm as

[TABLE]

We have the following lemma, which characterizes accuracy in a $\omega$-weighted norm given an approximation that is accurate in the $g$-weighted norm.

Lemma 3.2** (Strong convergence [9]).**

Let $p$ and $q$ be conjugate exponents, i.e., $1/p+1/q=1$ with $p,q\geq 1$. Assume the error of the approximation $f_{N}$ of $f$ satisfies

[TABLE]

Then,

[TABLE]

provided $C_{r}<\infty$.

This lemma states that the error induced by approximating using a polynomial construction that is accurate in a dominating measure induces an error that grows as the “difference” between the optimal and non-optimal basis increases. This difference is quantified as the maximum $C_{r}$ of the ratio of the two densities over $\Omega$. Thus if we use a PCE based upon the tensor product of univariate orthogonal polynomials for approximating a function of highly dependent variables then this can induce a substantial increase in error for a fixed sample size, when compared to the error obtained using a polynomial approximation that is designed to be accurate in the original $\omega$-weighted norm. We demonstrate in Section 5 that this happens in practice.

3.3 Orthogonalization methods

The most straightforward way to construct a PCE basis for dependent variables uses Gram-Schmidt orthogonalization [62]. Assume any linearly independent set of polynomials $\{\psi_{n}\}_{n=1}^{N}$ is given (a common choice might be polynomials that are orthonormal under a related tensorial measure). Then the polynomials $\phi_{n}$ orthonormal under $\omega$ are constructed numerically by setting $\phi_{1}(z)\equiv 1$ and computing

[TABLE]

The $L^{2}_{\omega}$ norms and inner products above can be computed approximately using quadrature. Let $\mathcal{Z}_{J}=(z^{(1)},\ldots,z^{(J)})$, with $z^{(m)}=(z_{1}^{(m)},\ldots,z_{d}^{(m)})$, and $w=(w_{1},\ldots,w_{J})$ be a set of quadrature samples and weights satisfying

[TABLE]

and let $\Psi$ denote the $J\times N$ Vandermonde-type matrix

[TABLE]

Note that designing a quadrature rule with small $J$ so that (15) holds is a very difficult task in general. Depending on the approach used, we expect the requisite $J$ to scale with dimension $d$ and with the desired accuracy of (15).

We assume that the quadrature rule $(\mathcal{Z}_{J},w)$ forms a proper discrete $\ell^{2}$ norm on the span of $\psi_{n}$, or equivalently

[TABLE]

Under this assumption, the columns of $\Psi$ weighted by $w$ are linearly independent, and so there is a unique QR factorization

[TABLE]

which is effectively performing the operations in (14). In the remainder of the paper, we will refer to $\sqrt{W}\Psi=QR$ as the moment matrix. Under these assumptions

[TABLE]

defines a polynomial basis orthonormal with respect to the discrete inner product defined by the quadrature rule in (15). For large $J$, this approximates the arbitrary measure $\omega$. The polynomials $\phi_{n}$ have the following properties, which follow from (16), (17), and (18):

• •

$\mathrm{span}\{\phi_{n}\}_{n=1}^{P}=\mathrm{span}\{\psi_{n}\}_{n=1}^{P}$ for any $P\leq N$.

• •

If there is a multi-index set $\Lambda$ such that $\{\psi_{n}\}_{n=1}^{N}$ is a basis for $\pi_{\Lambda}$, then $\{\phi_{n}\}_{n=1}^{N}$ is also a basis for $\pi_{\Lambda}$.

• •

$\sum_{j=1}^{J}w_{j}\phi_{n}(z^{(j)})\phi_{m}(z^{(j)})=\delta_{m,n}$.

From (18) it is clear that the basis $\phi_{n}$ depends on the ordering of the polynomials $\{\psi_{n}\}_{n=1}^{N}$. In the univariate setting, the technique described above is referred to as the Stieltjes procedure [23].

The numerical conditioning of the aforementioned Gram-Schmidt procedure is dependent on the accuracy of the quadrature rule used in (15). In low dimensions when $\Omega$ is tensorial and when the density $\omega$ is known explicitly, we can use tensor-product Gauss quadrature rules to form $\mathcal{Z}_{J}$ and $w$. Specifically, to compute a quadrature rule for a dependent measure, we construct a tensor-product quadrature rule for an independent measure $\nu$ that dominates the multivariate dependent measure $\omega$111By dominate we mean that the dependent measure is absolutely continuous with respect to the tensor-product measure.. Given the tensor-product rule with points and weights $\{z^{(q)},v_{q}\}_{q=1}^{Q}$, we can compute the weighted $L^{2}_{\omega}$ inner product necessary to orthogonalize the tensor product basis by a change of variable such that

[TABLE]

The use of a tensor-product quadrature rule is of course inefficient in general, however it only requires evaluations of polynomials in this context.

In higher-dimensions, one can use custom polynomial quadrature rules (e.g.  [30]) in the Gram-Schmidt procedure. However, often in these high dimensional settings or when Bayesian inference is used to condition prior estimates of uncertainty on observational data, the joint density of $Z$ is not known explicitly, but is instead characterized by samples drawn from the unknown distribution. In these situations we must use Monte Carlo or Psuedo Monte Carlo quadrature, where in (15) we use $w_{i}=1/J$.

4 Sampling schemes for dependent probability measures

Numerous methods exist for computing the coefficients of a polynomial chaos expansion. In this paper we focus on interpolation-based collocation methods. In the collocation setting, we have $N=M$, where $N$ is the dimension of the approximation space and $M$ is the number of samples. Given a set of $N$ realizations $\mathcal{Z}_{N}=\{z^{(1)},\ldots,z^{(N)}\}$, with corresponding model outputs $y=(f(z^{(1)}),\ldots,f(z^{(N)}))^{T}$, we would like to find a solution that satisfies $\Phi\alpha=y$, where $\alpha=(\alpha_{1},\ldots,\alpha_{N})^{T}$ denotes the vector of PCE coefficients and $\Phi\in\mathbb{R}^{N\times N}$ denotes the Vandermonde matrix with entries $\Phi_{mn}=\phi_{n}(z^{(m)}),\quad m\in[N],\;n\in[N]$.

Given data $y$ on $\mathcal{Z}_{N}$, we wish to construct a unique polynomial interpolant $f_{N}(Z)$ from $\pi_{\Lambda}$ that interpolates $y$. We assume unisolvence of the interpolation problem for $\mathcal{Z}_{N}$ on $\pi_{\Lambda}$; we will later prescribe a method that numerically guarantees this condition. Assuming unisolvence, there is a unique solution to the linear system:

[TABLE]

where the latter expression is true when $v(\cdot)$ is a non-vanishing weight function on $\mathcal{Z}_{N}$. The linear algebraic formulation above can be used to directly motivate a sampling scheme: build $\mathcal{Z}_{M}$ sequentially so that the determinant of $V\Phi$ is maximized. Such a procedure is a particular kind of D-optimal design [21], but differs from standard D-optimal design constructions in that (i) the procedure is sequential, i.e., greedy maximization is performed so that a sequence, instead of a non-nested grid, is constructed, and (ii) the choice of weight $v$ determines special properties of the sequence, as discussed in the following section. The reason we formulate the interpolation problem (20) with the additional weight $v$ is both to improve numerical stability, and to form connections with theoretical asymptotic results. For example, the results in [41] show that choosing $v=\sqrt{\omega}$ for the construction of weighted Leja sequences produces nodal sets that are asymptotically optimal.

4.1 Weighted Leja sequences via LU factorization

A Leja sequence (LS) is essentially a doubly-greedy computation of a determinant maximization procedure. Given an existing set of nodes $\mathcal{Z}_{M}$, a Leja sequence update chooses a new node $z^{(M+1)}$ by maximizing the determinant of a new Vandermonde-like matrix with an additional row and column: the additional column is formed by adding a single predetermined new basis element, $\phi_{M+1}$, and the additional row is defined by the newly added point. Hence a LS is both greedy in the chosen interpolation points, and also assumes some a priori ordering of the basis elements. The introduction of a row-based weighting makes this a weighted Leja sequence; as the previous section suggests, our weight will be the function $v$; we leave this function undefined for now, but make an explicit choice in Section 4.3.

In one dimension, a weighted LS can be understood without linear algebra: Let $\mathcal{Z}_{N}$ be a set of nodes on $\Omega$ with cardinality $N\geq 1$. We will add a new point $z^{(N+1)}$ to $\mathcal{Z}$ determined by the following:

[TABLE]

We omit notation indicating the dependence of $z^{N+1}$ on $\mathcal{Z}_{N}$. By iterating (21), one progressively builds up the Leja sequence $\mathcal{Z}$ by recomputing and maximizing the objective function for increasing $N$. The literature contains many instances of using the above objective function to add nodes [34, 52, 3].

In multiple dimensions, formulating a generalization of the univariate procedure is challenging. The following linear algebra formulation [55, 6] greedily maximizes the weighted Vandermonde-like determinant

[TABLE]

The above procedure is an optimization with no known explicit solution, so constructing a Leja sequence is challenging. In [42], gradient based optimization was used to construct weighted Leja sequences. However a simpler procedure based upon LU factorization can also be used [6]. The simpler approach comes at a cost of slight degradation in the achieved determinant of the LS. We adopt the LU-based approach here due to its ease of implementation.

The algorithm for generating weighted Leja sequences using LU factorization is outlined in Algorithm 1. The algorithm consists of 5 steps. First a polynomial basis must be specified. The number of polynomial basis elements must be greater than or equal to the number of desired samples in the Leja sequence, i.e. $N\geq M$. The input basis must also be ordered, and the Leja sequence is dependent on this ordering. Unless otherwise specified, in this paper we only consider total-degree polynomial spaces, that is we have

[TABLE]

for some polynomial degree $k$. We use lexigraphical ordering on $\Lambda$ to define the basis. The second step consists of generating a set of $S$ candidate samples $\mathcal{Z}_{S}$; ideally, $S\gg M$. Our candidate samples will be generated as independent and identically-distributed realizations of a random variable. The precise choice of the random draw will be discussed in the next section. For now we only require that the measure of the draw have support identical with the measure of $Z$. Once candidates have been generated we then form the $S\times N$ Vandermonde-like matrix $\Phi$, precondition this matrix with $V$, and compute a truncated LU factorization. (Computing the full LU factorization is expensive and unnecessary.) We terminate the LU factorization algorithm after computing the first $M$ pivots. These ordered pivots correspond to indices in the candidate samples that will make up the Leja sequence. If we assume that there is any size-$M$ subset of $\mathcal{Z}_{S}$ that is unisolvent for interpolation, then by the pivoting procedure, a Leja sequence is always chosen so that the interpolation problem is unisolvent.

Once a Leja sequence $\mathcal{Z}_{M}$ has been generated one can easily generate a polynomial interpolant with two simple steps. The first step evaluates the function at the samples in the sequence, i.e. $y=f(\mathcal{Z})$. The coefficients of the PCE interpolant can then be computed via

[TABLE]

where the matrices $P$, $L$, and $U$ are identified in Algorithm 1.

We end this section by noting that (approximate) Fekete points are an alternative determinant-maximizing choice for interpolation points [55, 6, 5]. We opt to use Leja sequences here because they are indeed a sequence, whereas a Fekete point construction is not hierarchical.

4.2 The induced measure

Generating Leja sequences with large determinants using LU factorization requires generating a large number $S$ of candidate samples $\mathcal{Z}_{S}$. The only theoretical requirement on the distribution of these samples is that they are sampled over the domain of the random variables $\Omega$. When sampling directly from $\omega$ is not feasible, one could instead sample uniformly over $\Omega$. However, when a significant portion of the probability is concentrated within a small region of $\Omega$, then $S$ needs to be prohibitively large. Computationally, this manifests as an ill-conditioned matrix $\Phi$. Generation of Leja sequences does not require the evaluation of the expensive simulation model $f$, so that in principle candidate sets can be quite large. However, we wish to avoid LU factorizations of matrices with millions of rows. To reduce the number of candidate samples needed, we instead propose to sample from specific measures.

Let the samples $z^{(s)}$ be generated as iid realizations from a density $\nu$. In the context of maximizing determinants, an attractive choice of the measure $\nu$ for generating the samples is the so-called induced measure

[TABLE]

where $k(z)$ is known as the Christoffel function. This measure is a property only of $\omega$ and $\pi_{\Lambda}$, and not of the individual basis elements $\phi_{\lambda}$. Therefore, this biased density is well-defined even when $Z$ has dependent components. This measure has been shown to generate well-conditioned matrices $\Phi$ [11].

Unfortunately, generating samples from the so called induced measure $v$ requires explicit knowledge of the probability measure $\omega$ of the random variables. In some settings, only samples from the measure are available. To mitigate this issue in this paper, we recommend sampling from a large $\Lambda$-asymptotic distribution called the equilibrium measure. When $\Lambda$ is a total-degree space, the induced measure converges to the equilibrium measure as $N\rightarrow\infty$ [44]. In contrast to the induced measure, a simple closed form of the equilibrium measure is known for some distributions. For random variables with bounded densities $\omega$ on a hypercube, the equilibrium measure is the Chebyshev distribution. Unfortunately when the probability density $\omega$ has areas of high concentration in a hypercube, even generating samples from the equilibrium measure can fail to generate a good candidate set. This is because the equilibrium measure is asymptotically optimal, but may not be actually optimal for finite $\Lambda$. However, we found that enriching the equilibrium set with samples from the probability distribution $\omega$ produced excellent candidate sets, and we adopted this approach for all numerical examples in our paper.

Let the joint density of $Z$ be given by (9), where the marginal distributions $\omega_{i}$ are each univariate Beta random variables with parameter $(\alpha,\beta)=(2,5)$ and the correlation matrix $R^{V}$ with $R^{V}_{11}=R^{V}_{22}=1$ and $R^{V}_{ij}=-0.9,i\neq j$. The resulting PDF $\omega$ and the induced measure density $\omega_{\Lambda}$ for $\Lambda$ corresponding to degree-3 and degree-20 total degree spaces are depicted in Figure 1.

The induced measure assigns significant non-zero probability in regions where the probability of the random variables $\omega$ is small, and the difference between the joint and induced densities appears to increases with polynomial degree. This is consistent with theory that states that for variables on bounded convex domains, the induced distribution will converge (with degree) to a Chebyshev-like distribution [7], which concentrates samples on the boundary of the variable domain $\Omega$.

4.3 Weight function

The properties of a Leja sequence are greatly influenced by the choice of the weight function $v(Z)$. In [42] Leja sequences were generated by setting $v=\sqrt{\omega}$. However in this paper we use the root inverse of the Christoffel function:

[TABLE]

It is shown in [33] that the Christoffel function is the optimal choice of weight function for generating Fekete nodes in one dimension.

5 Numerical Results

In this section we demonstrate the efficacy of our proposed approach on a number of numerical examples.

To measure the performance of an approximation, we will use the $\omega$-weighted $\ell^{2}$ error on a set of test nodes. We generate a set of $S=10,000$ random samples $\{z^{(j)}\}_{j=1}^{S}\subset\Omega$ drawn from the density $\omega$. The error is computed as

[TABLE]

where $f$ is the exact function and $f_{N}$ is the interpolative approximation.

When considering dependent random variables with probability measures concentrated in a small region of the variable domain, we are careful not to generate misleading results. Functions that vary significantly in regions of high-probability are much harder to approximate in a weighted $L^{2}_{\omega}$ norm than functions that vary in regions of low probability. Here we consider the oscillatory Genz function

[TABLE]

This function has strong variation throughout the domain. However, in a further attempt to avoid constructing a function that only varies strongly in regions of low-probability, we set the coefficients $c$ and $d$ randomly. Specifically we draw $e$ and $c_{i},i\in[d]$ randomly from the uniform distribution on $[0,1]$ and set

[TABLE]

To facilitate notation, let $B(\cdot;\alpha,\beta)$ denote a tensor-product density function of a Beta random variable with parameters $\alpha$ and $\beta$:

[TABLE]

where $K$ is a normalizing constant to ensure that $B$ is a probability density. Note that we use the same parameters $\alpha,\beta$ for each dimension.

We will consider three approaches in our numerical results:

• •

Nataf – This is the approach outlined in (5) which uses the transformation $\mathcal{T}_{\text{nataf}}^{\text{gauss}}$ defined in Section 3.1.2. The approximation $g_{N}$ from (5) in $U$-space is constructed using weighted Leja sequences as discussed in Section 4.

• •

DOM $(\alpha,\beta)$ – This is a domination strategy, where the dominating density is $B(\cdot;\alpha,\beta)$. Again, we build approximations with respect to the dominating methods using weighted Leja sequences.

• •

GS $(\alpha,\beta)$ – This is the proposed strategy in this paper. We first construct a basis that is (approximately) $L^{2}_{\omega}$ orthonormal via the technique in Section 3.3. We then use weighted Leja sequences as described in Section 4 to construct a PCE approximation in this basis.

Unless otherwise stated we use $10,000$ candidate samples to build the Leja sequence, where half of these samples are from a tensor-product Chebyshev density, and the remaining half are drawn from the density $\omega$.

We note that all three approches use weighted Leja sequences to construct the approximation. Therefore, our examples are a direct comparison for our three strategies to handle dependent variables: mapping methods, domination methods, and orthogonalization methods.

5.1 Leja sequence for domination measures

In this section we investigate the impact of constructing a polynomial approximation using the dominating measure strategy of Section 3.2. Thus given the measure associated with the random variables, we first identify a dominating measure, construct a PCE approximation using an interpolation sequence and basis elements from that dominating measure, and then compute the error with respect to the random variable measure. Although this paper is an exposition on approximation strategies for dependent random variables, here we choose the measure of the random variables to have independent components to facilitate comparison with known optimal strategies for tensor-product approximation, which are not applicable when variable dependencies exist.

In Figure 2 (left) we consider a one-dimensional case. We compare interpolants of a univariate oscillatory Genz function parameterized by a single Beta random variable $Z$ with the density $\omega=B(\cdot;10,10)$. Our dominating measure $g$ will be the uniform measure over the same domain, i.e., $g=B(\cdot;1,1)$. We generate two interpolants. The first one, $f_{N}^{\omega}$, uses Gauss quadrature nodes $\mathcal{Z}_{N}^{\omega}$ from $\omega$. The second, $f_{N}^{g}$, uses Gauss quadrature nodes $\mathcal{Z}_{N}^{g}$ from $g$. However, the error for both is measured with the same formula and samples in (5) that are generated from the Jacobi measure. The Gauss-Jacobi quadrature samples are optimal in the univariate setting and will produce a unitary condition number. The sampling scheme of the Jacobi basis places samples in a way that balances stability with sampling in high-probability regions. The procedure with $g$ does this as well, but for the uniform density. Consequently the samples for the uniform measure appear more frequently in regions of low $\omega$-probability. This manifests itself in a larger error in regions of high $\omega$ probability than obtained using the optimal $\omega$ construction. It is also apparent from Figure 2 (left) that the approximation $f_{N}^{\omega}$ peforms poorly in regions of very low $\omega$ probability. However, as measured in the $L^{2}_{\omega}$ norm, this is allowable.

In Figure 2 (right) we also plot the convergence in the median error (over $10$ samples of $b_{i}$ and $e$) of PCE approximations of the analytic function (24) in three dimensions, $d=3$. We set the joint density as $\omega=B(\cdot;10,10)$, and construct a PCE via the measure domination technique in Section 3.2, with the dominating density $g=B(\cdot;\beta,\beta)$ set to be another tensor-product Beta measure over the same support. We use various values of $\beta<10$. As $\beta$ approaches $10$ the dominating measure approximation will become more efficient. To limit the effect of the sampling scheme we use tensor-products of the univariate Gauss quadrature rule for $g$, and set the approximation space as associated with a tensor-product index set approximation,

[TABLE]

For each choice of dominating measure we compute the constant $C_{r}$ from Lemma 3.2 as a measure of the distance between the measure of orthogonality of each polynomial basis and the probability measure $\omega$. The plot clearly shows that constructing an approximation from a measure that is not orthogonal to $\omega$ results in a degradation of accuracy, and the penalty grows as $C_{r}$ grows. Using a dominating measure effects the constant of convergence but not the rate of convergence which is consistent with Lemma 3.2.

5.2 Leja sequences using probabilistic transformations

In this section we explore the performance of interpolants constructed using the transformations $\mathcal{T}_{\text{nataf}}^{\text{gauss}}$ and $\mathcal{T}_{\text{nataf}}^{\text{unif}}$, defined in Section 3.1.2, along with the mapping procedure (5) and weighted Leja sequences.

Probabilistic transformations, such as the Nataf transformation, are highly non-linear and often introduces steep gradients into a function via the map composition. To highlight the effect of the increased non-linearity introduced by the Nataf transform, consider the oscillatory Genz function as a function of only one variable shown in Figure 3. In this case the Rosenblatt and $\mathcal{T}_{\text{nataf}}^{\text{unif}}$ transformations are equivalent and consist of simply applying inverse transform sampling via the distribution function to define the $Z$ to $U$ transformation and back. Here we approximate the function using Gauss quadrature nodes of polynomials orthonormal with respect to the variables $Z$ and $U$. These are optimal in one-dimension and will outperform even Fekete and Leja sequences. Thus, we have two approximations: $f_{N}^{\omega}(z)$ which is a Gauss quadrature interpolant built from the Gauss nodes of $\omega$; and $g_{N}(\mathcal{T}(z))$, where $g_{N}$ is a Gauss quadrature interpolant built from Gauss nodes of a uniform random variable, and subsequently mapped to $z$-space via composition with $\mathcal{T}=\mathcal{T}_{\text{nataf}}^{\text{unif}}$.

We see the strong non-linearity introduced by the transformation degrades the accuracy of the interpolant. We can also easily see that the approximation error is largest in the regions in which the non-linearity introduced is strongest. These regions occur in regions of significant probability and so the resulting $L^{2}_{\omega}$ error is larger when using the Nataf transformation than without.

Figure 4 depicts the two-dimensional version of the algebraic function defined in (24). The same figure also plots the function under the transformations $\mathcal{T}_{\text{nataf}}^{\text{gauss}}$ and $\mathcal{T}_{\text{nataf}}^{\text{unif}}$. Clearly both transformations increase the non-linearity of the function, which means that we must also increase the degree of the PCE and thus the number of samples needed to approximate the function accurately.

In Figure 4 we also superimpose a degree 12 Leja sequence constructed in the i.i.d uniform space, along with the same samples mapped back to the original space $Z$ on the right and left plots respectively. This Leja sequence produces an interpolation matrix with a small condition number, however this benefit is outweighed by the degradation in accuracy of the interpolant when compared to the dominating measure approach (see Figure 5). Note that we can also build Leja sequence using univariate Hermite polynomials in the i.i.d. standard normal space coupled with the transformation $\mathcal{T}_{\text{nataf}}^{\text{gauss}}$. But again due to the non-linearity of the Nataf transformation (see Figure 4), the accuracy of the interpolant is poor relative to the approximation constructed using the dominating measure approach. We omit this comparison for brevity.

Another way to interpret the slow convergence of the interpolant built using non-linear probabilistic transformations is that we are no-longer approximating with polynomials in the original space of $Z$, but rather with mapped polynomials, and hence do not have the same guarantees on convergence that exist when approximating with polynomials. The poor performance we observe here is consistent with the univariate results documented in [64] and the multivariate results reported in [19].

5.3 Leja sequences using GSO bases

In this section we compare our three approaches for dependent variables when the approximation is constructed using interpolation via weighted Leja sequences. To consider the effect of dimension we consider two problems with $d=2$ and $d=10$. In each case we assume that the joint density of $Z$ is given by (9) and that the marginal distributions $\omega_{i}$ are each univariate Beta random variables with parameters $(\alpha,\beta)=(2,5)$.

In the top row of Figure 5 we compare the Nataf, DOM, and GS approaches described earlier for building PCE interpolants using the test function (24) in two dimensions, where we set the entries of the correlation matrix $R^{V}$ to $R^{V}_{11}=R^{V}_{22}=1$ and $R^{V}_{ij}=-0.9,i\neq j$. This density is plotted in the left of Figure 1 and right of Figure 4.

With the exception of the Nataf interpolant, the rates of convergence in the median error (over $10$ samples of $b_{i}$ and $e$) are similar. However the improved constant of convergence in (13) obtained by using the GS-orthogonalized basis results in significant improvements in accuracy. The Leja sequences built using the GSO basis produce interpolants that have errors which are orders of magnitude smaller for a fixed sample size. Unlike the other approaches, the convergence rate of the error in the Nataf interpolant deteriorates significantly. The poor performance we observed here is consistent with results observed in the univariate case [64] and the multivariate results observed in [19]. In short, this is caused by a Jacobian with steep gradients, as explained in Section 5.2.

To assess the conditioning of the interpolation problem and the procedure used to construct the GS basis we respectively use the following condition numbers

[TABLE]

where $\sigma_{\max}(A)$ and $\sigma_{\min}(A)$ are the maximum and minimum singular values of a matrix A. In all cases the weighted Vandermonde-like matrices $V\Phi$ built using Leja sequences are well conditioned. This is despite the fact that the condition numbers $\kappa^{\text{GS}}_{N}$ of the moment matrices (17) used to compute the GSO basis are large. The condition number of the moment matrices does eventually affect accuracy of the interpolant, but only at extremely high degrees, specifically $p\geq 23$ for the GS $(1,1)$ case.

The condition number of the moment matrix is dependent on how disparate the orthogonality measure of the tensor product basis is from $\omega$. The condition number can be reduced by using a tensor product basis which is orthogonal to a measure which is “closer” to the probability measure $\omega$. This is evident by the reduction in condition number in the GS $(2,5)$ example. The problem tested here is quite challenging, i.e. the correlation of the variables is high and the ratio of the densities is large, but our proposed algorithm still performs well.

The bottom row of Figure 5 plots the convergence of the various interpolation strategies for the algebraic equation with $d=10$. Each off diagonal entry of ${R^{V}}$ is set to 0.9, then for $i=1,3,\ldots,9$, $j=1,\ldots,d$ we set ${R^{V}}_{ij}=-{R^{V}}_{ij}$, and ${R^{V}}_{ji}=-{R^{V}}_{ji}$. We observe similar trends to those in the 2D example. Specifically the GS approaches perform notably better than the other approaches. Again, the Nataf transformation performs poorly: after 1000 samples the error is still $\mathcal{O}(1)$. As in the 2D case, the GS $(2,5)$ case produces better Gram-Schmidt condition numbers than the GS $(1,1)$ case.

5.4 Using Leja sequences for quadrature

The construction of Leja points is motivated largely by interpolation; however it is straightforward to compute the mean and variance from a polynomial chaos expansion constructed via interpolation. Due to orthonormality, the mean and variance can be computed analytically from the PCE

[TABLE]

In the case when an approximate grid is used to orthogonalize polynomials, such as in (15), then the expressions above are only true with respect to the discrete measure defined by this grid. If the discrete grid is accurate, then the above can be reasonable approximations to expectations with respect to the continuous dependent density. If the PCE interpolates using a set of nodes $\mathcal{Z}_{N}$, then we can also easily compute a quadrature rule. Because we can compute mean and variance analytically using (26) such rules are not often needed. However the weights of the quadrature rule can be used to quantify the ill-conditioning of the estimates of integrals using PCE interpolants.

Given function data $y$, the polynomial chaos coefficients $\alpha_{n}$ satisfy

[TABLE]

Since $\alpha=\Phi^{-1}y$, and $\phi_{1}\equiv 1$ because $\omega$ is a probability density function,

[TABLE]

then we immediately conclude that the first row of the matrix $\Phi^{-1}$ gives us quadrature weights $v_{n}$ defining the Leja polynomial quadrature rule

[TABLE]

The subscript $N$ on $Q_{N}$ indicates that the quadrature rule has $N$ nodes.

Given the weights $v_{n}$, the $\ell^{1}$ condition number of the quadrature operator $Q_{N}$ is given by

[TABLE]

where the last equality holds under the assumption that the $\phi_{n}$ are orthonormal with respect to a probability density function $\omega$. Large values of $\kappa_{Q_{N}}$ indicate the presence of negative weights, which makes the computation susceptible to catastrophic cancellation.

In Figure 6, we consider the same examples as in Figure 5, but now plot the relative error in the mean of the PCE obtained using different Leja sequences and the quadrature rule condition numbers. The DOM $(1,1)$ and DOM $(2,5)$ approaches cannot directly estimate the mean since the basis elements are not orthogonal to the probability measure, so we do not plot the errors for these approaches (integrals can be estimated numerically by sampling on the interpolant from the correct measure). As observed when computing weighted $L^{2}_{\omega}$ errors, the Nataf transformation performs poorly when used for quadrature and the Gram-Schmidt Leja sequences produce the most accurate estimates. The condition number of all quadrature rules are observed to be $O(1)$. However we must note that again the condition number of the moment matrices used to construct the Gram-Schmidt basis affect accuracy at very high-degrees. In two-dimensions, we compute the true mean of the function to machine precision using Gauss-Legendre tensor product quadrature, that is, we evaluate $\int_{\Omega}f(z)\mathrm{d}\omega(z)$ by integrating $f(z)\omega(z)$ with respect to the uniform measure on $[0,1]^{2}$. In ten dimensions, tensor product quadrature is not feasible so we use Monte Carlo (MC) quadrature with $10^{6}$ samples drawn from the probability measure. The non-monotonic decrease in the error of the Gram-Schmidt PCE rules results from the error in the MC estimate of the mean being larger than the error in the PCE estimates.

5.5 Approximation without a closed form for the joint density

When building an interpolant of a simulation model, a closed form of the joint density of the random variables is often not available. Instead one may only have samples from the joint density. For example, Bayesian inference [57] is often used to infer densities of random variables conditional on available observational data. The resulting so-called posterior densities are almost never tensor-products of univariate densities. Moreover it is difficult to compute analytical expressions for the posterior density and so Markov Chain Monte Carlo (MCMC) sampling is often used to draw a set of random samples from the posterior density.

Here we investigate the impact of computing the moment matrix (17), used with the Gram-Schmidt orthogonalization, using Monte Carlo (MC) quadrature when no closed form of the joint density is known. In the top row of Figure 7 we compare the accuracy of the interpolation procedure, when applied to the 2D algebraic function, as the number $J$ of MC samples is varied. From the figure we can see that decreasing the number of MC samples increases the ill-conditioning of both the interpolation matrix and the moment matrix. When the condition numbers become sufficiently large the ill-conditioning also affects the accuracy of the interpolant. Clearly the onset of the degradation of the accuracy of the interpolant is dependent on the number of MC samples.

The bottom row of Figure 7 highlights the effect of the number of MC samples on the estimation of the mean of the 2D algebraic function. The accuracy of the mean obtained from the PCE is limited by the number of samples used to construct the orthonormal basis. This is consistent with the observations made in the univariate setting in [47]. The effect of the sample size appears to have a much greater impact on the accuracy of the quadrature procedure than it does on interpolation accuracy measured in the $L^{2}_{\omega}$ norm. This fact can be leveraged to accurately estimate statistics using the GSO procedure by randomly sampling on the PCE approximation. This sampling can be done with no additional expensive evaluations of the true function $f$. Note the curve depicting approximation with 100 MC samples does not extend to large degrees because a requirement of the Gram-Schmidt orthogonalization is that the moment matrix is over-determined, which is not true for large degrees when only 100 MC samples are used.

5.6 Approximating a chemical reaction model

In this section we will demonstrate the utility of using our approach for Bayesian inference and dimension reduction using a model of competing species absorbing onto a surface out of a gas phase [61]. Consider the following ordinary differential equation, prescribing evolution of the mass fractions $(u_{1},u_{2},u_{3})$ of three chemical species:

[TABLE]

for some constants $a$, $b$, $c$, $d$, $e$, and $f$ and initial conditions $u_{1}(0),u_{2}(0),u_{3}(0)$.

5.6.1 Bayesian inference

Bayesian inference [57] is often used to infer densities of random variables conditional on available observational data. To make this precise, let $f_{\text{o}}({Z}):\mathbb{R}^{d}\rightarrow\mathbb{R}^{n_{o}}$ be an observable quantity, parameterized by the same $d$ random variables ${Z}$, which predicts a set of $n_{o}$ observable quantities. Bayes rule can be used to define the posterior density for the model parameters ${Z}$ given observational data ${{y}_{o}}$:

[TABLE]

where any prior knowledge on the model parameters is captured through the prior density $\pi(z)$. The function $\pi(y_{0}\mid z)$ is the likelihood function and dictates an assumed model-versus-data misfit.

The construction of the posterior is often not the end goal of an analysis. Instead, one is often interested in statistics on the unobservable QoI. Here we will focus on approximation of the data-informed predictive distribution of the mass fraction of the first species $u_{1}(t=50)$ at $t=50$ seconds, which we assume we cannot measure.

In the following we will assume that the rate parameters $a=2(Z_{1}+3)/3\in[0,4]$ and $b=30(Z_{2}+2)/7+5\in[5,35]$ in (28) are random variables where the posterior distribution (29) of $Z$ is

[TABLE]

where $C$ is a constant chosen to normalize $\omega$ so that it is a probability density. This density is called a “banana density” and is a truncated non-linear transformation of a bivariate standard normal Gaussian distribution. The response surface of the mass fraction over the posterior density is depicted in the left of Figure 8. The response has a strong non-linearity which makes it ideal for testing high-order polynomial interpolation.

To approximate the response surface we use the GS $(1,1)$ method. We compute the moments needed to orthogonalize the basis using Gauss-Legendre quadrature (including the weight in the integrand), and compare with the same approach using Monte Carlo quadrature. We use rejection sampling to draw independent samples from the banana density to avoid the effects of reduced effective sampling size that arise when using MCMC sampling. In the right of Figure 8 we plot the $L^{2}_{\omega}$-errors for various total-degree polynomial interpolants (up to degree 20) of the chemical reaction model using the banana density. Similar to previous results, the GS basis obtained using exact moments achieves the smallest error for a fixed number of samples. Using samples from the density to estimate moments achieves comparable accuracy for small sample sizes until ill-conditioning resulting from the low accuracy of the moments takes effect.

We remark that previous attempts at building an approximation to increase the efficiency of Bayesian inference can be found in [13, 37, 39]. Two of these methods [13, 39] use low-order localized surrogates to facilitate sampling in regions of high-probability. The use of localized surrogates results in small-rates of convergence. This is in contrast to the method proposed in this paper, which can exhibit spectral convergence rates. The method presented in [37] uses a sequence of global approximations that eventually concentrate on the posterior distribution. However unlike our proposed approach, the method in [37] does not use all model evaluations in the final approximation. Many samples are only used to construct the intermediate approximations.

5.6.2 Dimension reduction

Using the same chemical species problem (28), we show here that our algorithm can be used to integrate high-dimensional ridge functions. Ridge functions are multivariate functions that can be expressed as a function of a small number of linear combinations of the input variables variables [49]. We define a ridge function to be a function of the form $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that can be expressed as a function $g$ of $s<d$ rotated variables,

[TABLE]

In applications it is common for $d$ to be very large, but for $f$ to be (an approximate) ridge function with $s\ll d$. When approximating a ridge function one does not need to build a surrogate in $\mathbb{R}^{d}$ but rather can focus on the more tractable problem of approximating in $\mathbb{R}^{s}$. The difficulty then becomes generating an orthonormal basis with respect to the low-dimensional transformed probability density, which is typically unknown and challenging to compute. When the $d$-dimensional space is a hypercube, then the corresponding $s$-dimensional probability density is defined on a multivariate zonotope, i.e., a convex, centrally symmetric polytope that is the $s$-dimensional linear projection of a $d$-dimensional hypercube. The vertices of the zonotope are a subset of the vertices of the $d$-dimensional hypercube projected onto the $s$-dimensional space via the matrix $A$.

In the following we will consider the approximation of a ridge function of independent uniform random variables with uniform density $v(y)$. We set the ridge function to be a function of the mass fraction of the first species $u_{1}$ of the competing species model (28). Specifically we define $z=Ay\in\mathbb{R}^{2}$, where $y\in\mathbb{R}^{20}$ and $A\in\mathbb{R}^{2\times 20}$ is a randomly generated matrix with orthogonal rows, and we define the ridge function to be $f(y)=u_{1}(Ay)$. The projection of the independent high-dimensional density $v$ into the lower dimensional space induces a new dependent probability density $\omega$ on a zonotope $\Omega$. The resulting zonotope and probability density are shown in the left of Figure 9. The zonotope density is obtained from a KDE constructed using $10^{4}$ samples $Z=AY$, where $Y$ were drawn randomly from the $d=20$ uniform distribution.

To build an orthonormal basis on the zonotope we must select a tensor-product basis, a quadrature rule to compute the moments of that basis, and a procedure for generating candidate samples for the Leja sequence. Here we use a monomial basis and Sobol sequences to compute moments of the monomial basis with respect to $\omega$. Specifically we generate a Sobol sequence of $J=10^{4}$ $d=20$-samples $\{y^{(j)}\}_{j=1}^{J}$ and then project these samples onto the zonotope to generate a quadrature rule $\mathcal{z}^{J},w$, where $z^{(j)}=Ay^{(j)}$ and $w_{q}=1/J$. To generate candidate samples we use a randomized algorithm [56] to find the vertices of the zonotope and use rejection sampling to sample inside the convex hull of the zonotope vertices.

In the right of Figure 9 we plot the $L^{2}_{\omega}$-errors in various polynomial total-degree interpolants (up to degree 20) of the chemical reaction model ridge function. We compare the GS procedure against a DOM $(1,1)$ procedure. Again the GS basis achieves a significantly smaller error, for a fixed number of samples, than when using a tensor product basis. In this example the accuracy of the moments obtained by Sobol quadrature is sufficient to avoid noticeable effects of ill-conditioning for the ranges of degrees considered.

5.7 Partial differential equations with random input data

In this section we apply our proposed methodology to build an approximation of a QoI obtained from a model of diffusion in porous media. We utilize this example to explore the effect of dimensionality on our proposed methodology.

Consider the following model of diffusion

[TABLE]

where the random diffusivity $k$ is a random field represented by the expansion

[TABLE]

where

[TABLE]

Given a correlation length $L_{c}$, which controls the variability of the random field, we set $L_{p}=\max(1,2L_{c})$ and $L=\frac{L_{c}}{L_{p}}$.

The major challenge of polynomial approximation is the fast growth of the polynomial basis with the number of random variables. To mitigate this issue, we use an anisoptropic index set that consists of higher degree terms in important directions and lower-degree approximations in the less important directions. Specifically we use the following index sets presented in [45]:

[TABLE]

where $\alpha_{\text{min}}=\min_{k}\alpha_{k}$ and

[TABLE]

To study the convergence of the approximation methods discussed in this paper, we consider problems of varying number of variables $d$ and investigate the convergence of the $L^{2}_{\omega}$ error as the number of points used to build each approximation increases. With this goal, we set $L_{c}=0.5$ and the joint density of $Z$ to be the mixture of two independent tensor-product beta distributions222Note that, although the individual components of the mixture are independent, the mixture of the components is not. Indeed, the distribution is bi-modal., i.e. $\omega(z)=\frac{1}{2}B(z,10,4)+B(z,4,10)$. In Figure 10, we plot the error (computed using 1000 validation samples) as a function of the number of Leja samples $M$, which is dependent on $l$ in (34).

From the left plot, we see that the rate of convergence of our GS approach333We use sparse grid quadrature to produce an orthonormal basis. The index set (34) can be exactly integrated by sparse grids for any tensor-product probability measure using weighted Leja sequences [41]. When $\omega$ is a mixture of independent measures, we can simply use a separate sparse grid quadrature rule for each mixture component and estimate the true mixture integral using a weighted sum of the individual quadrature rules. does not degrade as the number random variables increases for $d>11$. This is consistent with the behavior observed for sparse grid approximations (with i.i.d uniform variables) in [45]. In the middle and left plots, we observe that the GS-based PCE can be significantly more accurate than the DOM-based PCE. The accuracy of the GS PCE is degraded, relative to DOM, only when the accuracy of the quadrature rule is insufficient to accurately othonormalize the tensor product basis used in (14) or when insufficient candidates are used to build the Leja sequence. There is a mild dependence on dimension for both the number of Monte Carlo samples needed to produce a well-conditioned GSO basis and the number of candidate samples used to construct a conditioned Leja sequence. These two effects are evident by the deterioration of the error in the GS (1,1) $J=2\cdot 10^{3}$ and $\texttt{GS}^{\star}$ approximations, respectively, when the number of variables is increased from 11 to 41.

6 Conclusions

In this paper, we have presented an algorithm for building polynomial interpolants of functions of dependent variables. Most existing literature focuses on strategies for building polynomial approximations under the assumption of independence. Our work provides two major contributions to the existing literature, namely an investigation of the disadvantages of existing polynomial interpolation methods for approximation with dependent variables and an algorithm for generating interpolants that minimizes the ill-conditioning of the interpolation matrix.

We demonstrate through extensive numerical examples that our results are almost always significantly more accurate than existing approaches. The improved performance is obtained by balancing the need to sample in high-dimensions with minimizing the ill-conditioning of the polynomial interpolation matrix. The efficacy of our approach is only limited by the ability to produce multivariate polynomials which are orthonormal to the joint distribution of the random variables. We use Gram-Schmidt orthogonalization to generate such polynomials. The stability of this procedure is dependent on the dependence structure between the random variables and the accuracy of the quadrature rules used to compute the inner-products in the Gram-Schmidt procedure.

The ill-conditioning of Gram-Schmidt introduced by using approximate numerical integration can be minimized by using highly accurate quadrature rules, provided an explicit form for the joint density is available. However, in a number of practical situations, one does not know the joint density explicitly, but rather only has samples drawn from its distribution. In these settings, for example when Bayesian inference is used to condition prior estimates of uncertainty on data, we can still compute highly accurate interpolants, but the accuracy of our approach is dependent on the number of available distribution samples. In future work, we will investigate alternative strategies for producing orthonormal multivariate polynomials that do not suffer from the ill-conditioning associated with the Gram-Schmidt procedure.

Our approach is extremely flexible and can exploit various types of structure in the function being approximated. We demonstrated this by leveraging dimension reduction and a priori estimates of anisotropy to build interpolants that can be used for high-dimensional approximation. In future work, we plan to develop adaptive strategies for determining the best polynomial index set that minimizes cost while maximizing accuracy, when anisotropy must instead by discovered.

7 Acknowledgements

J. Jakeman and M. Eldred were supported by DARPA EQUiPS and DOE SCIDAC. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525. The views expressed in the article do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

A. Narayan is partially supported by NSF DMS-1720416, AFOSR FA9550-15-1-0467, and DARPA EQUiPS N660011524053.

F. Franzelin and D. Pflüger were supported by the German Reseach Foundation (DFG) within the Cluster of Excellence in Simulation Technology (EXC 310).

8 Appendix

8.1 Generating samples using the Gaussian copula

The following algorithm can be used to generate dependent samples, from a Gaussian copula with correlation matrix ${R^{V}}$, with marginals, with marginal probability densities given by $f_{i}(z_{i})$. This algorithm can be used to generate random variables which satisfy the assumptions of the Nataf transformation 3.1.2.

Generate a sample $u=(u_{1},\ldots,u_{d})$ from the multivariate Gaussian distribution $N(0,I)$ with zero mean and unit variance. 2. 2.

Compute the Cholesky factorization $L$ of the correlation matrix such that ${R^{V}}=LL^{T}$. 3. 3.

Compute a correlated standard normal sample $v=Lu$. 4. 4.

Modify the sample $v$ to have uniform marginals in every dimension. I.e. compute $\hat{v}_{i}=\Phi(v_{i})$, where $\Phi$ is the CDF of the standard normal distribution. 5. 5.

Generate a sample $z=(z_{1},\ldots,z_{d})$ from the desired distribution. I.e. compute $z_{i}=F_{i}(\hat{v}_{i})$, where $F_{i}$ is the CDF of the desired marginal density.