Extending bluff-and-fix estimates for polynomial chaos expansions

Highlights

• An iterative polynomial chaos algorithm is applied to approximate the solution to a nonlinear PDE with stochastic parameters.

• The solution to a smaller polynomial chaos $M_0$ system is iteratively used to determine the solution to corresponding $M$ systems for a consecutive range of $M>M_0$.

• The algorithm accommodates choice around which component variables in a PDE system are higher priority to estimate accurately, allowing for efficient performance when the objective is to evaluate a function of the PDE solution that does not necessarily weigh the importance of each variable equally.

Abstract

Polynomial chaos methods, which are part of a broader class known as stochastic Galerkin schemes, can be used to approximate the solution to a PDE with uncertainties represented by stochastic inputs or parameters. The stochastic solution is expressed as an infinite polynomial expansion truncated to $M+1$ terms. The approach is then to derive a resulting system of coupled, deterministic PDEs and solve this system with standard numerical techniques. Some challenges with conventional numeric techniques applied in this context are as follows: (1) the solution to a polynomial chaos $M$ system cannot easily reuse an already existing computer solution to an $M_0$ system for some $M_0<M,$ and (2) there is no flexibility around choosing which variables in an $M$ system are more important or advantageous to estimate accurately. This latter point is especially relevant when, rather than focusing on the PDE solution itself, the objective is to approximate some function of the PDE solution that weights the solution variables with relative levels of importance. In Lyman and Iaccarino (2020) [1], we first present a promising iterative strategy (bluff-and-fix) to address challenge (1); we find that numerical estimates of the accuracy and efficiency demonstrate that bluff-and-fix can be more effective than using monolithic methods to solve a whole $M$ system directly. This paper is an extended version of Lyman and Iaccarino (2020) [1] that showcases how bluff-and-fix successfully addresses challenge (2) as well by allowing for choice in which variables are approximated more accurately, in particular when estimating statistical properties such as the mean and variance of an $M$ system solution.

Introduction

Uncertainty quantification (UQ) in physical models, including those governed by systems of partial differential equations, is important for building confidence in the resulting predictions. Uncertainties can originate from imperfect knowledge of boundary or initial conditions in the system of interest or variability in the operating conditions. A common approach is to represent the sources of uncertainty as stochastic variables; in this context the solution to the original differential equation becomes a random quantity. Stochastic Galerkin schemes (SGS) are used to approximate the solution to parametrized differential equations [2], [3], [4]. In particular, they utilize a functional basis in the parameter space to express the solution and then derive and solve a deterministic system of PDEs with standard numerical techniques [5], [6]. A Galerkin method projects the randomness in a solution onto a finite-dimensional basis, making deterministic computations possible. SGS are part of a broader class known as spectral methods [7], [8], [9].

The most common UQ strategies involve Monte Carlo (MC) algorithms, which suffer from a slow convergence rate proportional to the inverse square root of the number of samples [10]. If each sample evaluation is expensive — as is often the case for the solution of a PDE — this slow convergence can make obtaining tens of thousands of evaluations computationally infeasible [8]. Initial spectral method applications to UQ problems showed orders-of-magnitude reductions in the cost needed to estimate statistics with comparable accuracy [11].

In the conventional approach for SGS, the unknown quantities are expressed as infinite series of orthogonal polynomials in the space of the random input variables. This representation has its roots in the work of Wiener [12], who expressed a Gaussian process as a countably infinite series of Hermite polynomials. Ghanem and Spanos [7] truncated Wiener’s representation and used the resulting finite series as a key ingredient in a stochastic finite element method. SGS based on polynomial expansions are often referred to as polynomial chaos approaches. In contrast to sampling methods (e.g. Monte Carlo simulations), polynomial chaos approaches are intrusive methods i.e. they require the solution of a mathematical problem that is different from the one originally considered. Uncertain quantities are represented as expansions that separate deterministic coefficients from a chosen random orthogonal basis [13].

Let $\mathcal{D}=[0,1]\times [0,T]$ be a subset of the spatial and time domain ${\mathbb{R}}_x\times {\mathbb{R}}_{t\ge 0}$.¹ Then let $u:\mathcal{D}\to \mathbb{R}$ be continuous and differential in its space and time components; further, let $u\in {\mathcal{L}}^2\left(\mathcal{D}\right)$. This $u$ represents the solution to a differential equation,

$$\mathcal{F}(u,x,t)=0.$$

Here $\mathcal{F}$ is a general differential operator that contains both spatial and temporal derivatives. Often $\mathcal{F}$ is assumed to be nonlinear.

Let $\xi$ be a zero mean, square-integrable, real random variable. We assume uncertainty is present in the initial condition $u(\cdot ,0)$ and represent it by setting

$$u(x,0;\xi ):{\mathbb{R}}_x\to \mathbb{R}u(x,0;\xi )=f(x,\xi )$$

where $f$ is a known function of $x$ and $\xi$. Accordingly, the solution $u(x,t;\xi )$ to $\mathcal{F}(u,x,t)$ is now a random variable indexed by $(x,t)\in \mathcal{D}$, meaning $u(x,t;\xi )$ is a stochastic process.

As statistics of $\xi$, we require that both $u(x,t;\xi )$ and $f(x,\xi )$ have existing second moments. Note that these are the only restrictions; namely, even though $f$ is chosen as sinusoidal in the example of Section 2, we do not require $f$ to be periodic, bounded over the real line, zero on the whole boundary $\partial \mathcal{D}$, etc.

We consider a polynomial chaos expansion (PCE) and separate the deterministic and random components of $u$ by writing

$$u(x,t;\xi )=\sum _{k=0}^{\infty }{u}_k(x,t){\mathrm{\Psi }}_k\left(\xi \right).$$

The $u_k:\mathcal{D}\to \mathbb{R}$ are deterministic coefficients, while the ${\mathrm{\Psi }}_k$ are orthogonal polynomials with respect to the measure $d\xi$ induced by $\xi$. Let $〈\cdot ,\cdot 〉$ denote the inner product mapping $(\varphi ,\psi )\mapsto \int \varphi \left(\xi \right)\psi \left(\xi \right)d\xi$. The triple-product notation $〈\varphi \psi \phi 〉$ is understood as $〈\varphi ,\psi \phi 〉=〈\varphi \psi ,\phi 〉=\int \varphi \left(\xi \right)\psi \left(\xi \right)\phi \left(\xi \right)d\xi$ and the singleton $〈\varphi 〉$ as $\int \varphi \left(\xi \right)d\xi$. Then we require the ${\mathrm{\Psi }}_k$ to satisfy the properties [5]

$$〈{\mathrm{\Psi }}_k〉={\mathbb{E}}_{\xi }\left({\mathrm{\Psi }}_k\right)={\delta }_{k0}〈{\mathrm{\Psi }}_i{\mathrm{\Psi }}_j〉=c_i{\delta }_{ij}$$

where the $c_i\in \mathbb{R}$ are nonzero and ${\delta }_{ij}$ is the Kronecker delta. Note that the Hermite ${\mathrm{\Psi }}_k$ are necessarily independent as random variables in $\xi$ from their orthogonality.

Let $S_M(x,t;\xi )={\sum }_{k=0}^M{u}_k(x,t){\mathrm{\Psi }}_k\left(\xi \right)$ be a stochastic process indexed by $(x,t)\in \mathcal{D}$. By the Cameron–Martin theorem² the truncated PCE denoted by $S_M$ converges in mean square,³

$$S_M(x,t;\xi )\underset{{\mathcal{L}}^2}{\stackrel{M\to \infty }{\stackrel{}{⟶}}}u(x,t;\xi ),$$

at every fixed $(x,t)$ in the domain $\mathcal{D}$. This justifies the PCE and its truncation to a finite number of terms for the sake of computation.

Stochastic properties of the solution $u$ can be readily computed [8], [15], [17], [18] given the attributes in (2). Namely,

$$\mathbb{E}(x,t)≔{\mathbb{E}}_{\xi }\left(u(x,t;\xi )\right)=\sum _{k=0}^{\infty }{u}_k(x,t){\mathbb{E}}_{\xi }\left({\mathrm{\Psi }}_k\left(\xi \right)\right)=u_0(x,t).$$

Going forward, the symbol $\mathbb{E}$ denotes a deterministic ${\mathcal{C}}^2\left(\mathcal{D}\right)$ function in $x$ and $t$ in addition to the usual expectation operator.

Similarly, the variance of the solution is

$$\mathbb{V}(x,t)≔{\mathbb{V}}_{\xi }\left(u(x,t;\xi )\right)=\sum _{k=1}^{\infty }{u}_k{(x,t)}^2〈{\mathrm{\Psi }}_k{\mathrm{\Psi }}_k〉=\sum _{k=1}^{\infty }{c}_k{u}_k{(x,t)}^2.$$

Note that variance is defined, because we assume $u(x,t;\xi )$ has existing second moments over all $(x,t)$ indexes. Using the truncated $u$ expansion $S_M$, we can approximate the solution’s true variance ${\mathbb{V}}_{\xi }\left(u(x,t;\xi )\right)$ as⁴

$${\mathbb{V}}_M(x,t)≔{\mathbb{V}}_{\xi }\left(S_M(x,t;\xi )\right)={\mathbb{V}}_{\xi }\left(\sum _{k=0}^M{u}_k(x,t){\mathrm{\Psi }}_k\left(\xi \right)\right)\approx \mathbb{V}\left(u(x,t;\xi )\right).$$

As mentioned in the abstract, this paper is an extended version of [1] that considers a PDE system obtained by explicitly representing the PDEs in terms of ${\mathbb{V}}_M(x,t)$ and $\mathbb{E}(x,t)$. In Section 3.4, our objective will not be to compute $S_M$ but instead to estimate $\mathbb{E}$ and ${\mathbb{V}}_M$. The beauty is that, unlike the case of approximating $u$ by $S_M$, not all of the coefficients functions $u_k(x,t)$ need to be determined perfectly in order to approximate $\mathbb{E}$ and ${\mathbb{V}}_M$ accurately. In this sense, the algorithm we introduce (bluff-and-fix) is preferable to conventional methods that do not have the flexibility to choose which $u_k(x,t)$ are more important to estimate well.

Substituting the truncation $S_M$ into Eq. (1), we have

$$\mathcal{F}\left(\sum _{k=0}^M{u}_k(x,t){\mathrm{\Psi }}_k\left(\xi \right),x,t\right)=0.$$

Furthermore, we can determine the initial conditions for the deterministic component functions. Multiplying $u(x,0;\xi )$ by any ${\mathrm{\Psi }}_k$ and integrating with respect to the $\xi$-measure $d\xi$ yields

$$\sum _{i=0}^{\infty }{u}_i(x,0)\int {\mathrm{\Psi }}_i\left(\xi \right){\mathrm{\Psi }}_k\left(\xi \right)d\xi =\int f(x,\xi ){\mathrm{\Psi }}_k\left(\xi \right)d\xi ={\mathbb{E}}_{\xi }\left[f(x,\xi ){\mathrm{\Psi }}_k\left(\xi \right)\right]$$$$u_k(x,0)=\frac{1}{c_k}{\mathbb{E}}_{\xi }\left[f(x,\xi ){\mathrm{\Psi }}_k\left(\xi \right)\right].$$

The scalars $c_k$ of course are dependent on the choice of the orthogonal ${\mathrm{\Psi }}_k$ polynomials. Consequently, the initial conditions for $\mathbb{E}(x,t)$ and ${\mathbb{V}}_M(x,t)$ are

$$\mathbb{E}(x,0)=u_0(x,0)=\frac{1}{c_0}{\mathbb{E}}_{\xi }\left[f(x,\xi ){\mathrm{\Psi }}_0\left(\xi \right)\right],$$$${\mathbb{V}}_M(x,0)=\sum _{k=1}^M{c}_k{u}_k{(x,0)}^2=\sum _{k=1}^M\frac{1}{c_k}{\mathbb{E}}_{\xi }{\left[f(x,\xi ){\mathrm{\Psi }}_k\left(\xi \right)\right]}^2.$$

In a similar manner, we can “integrate away” the randomness in Eq. (7) by projecting onto each basis polynomial. This is discussed in detail in the next section.