0
Research Papers

[+] Author and Article Information
Panagiotis A. Tsilifis

CSQI,
Institute of Mathematics,
School of Basic Sciences,
École Polytechnique Fédérale de Lausanne,
Lausanne CH-1015, Switzerland
e-mail: panagiotis.tsilifis@epfl.ch

Manuscript received November 10, 2016; final manuscript received July 5, 2018; published online July 30, 2018. Assoc. Editor: Sez Atamturktur.

J. Verif. Valid. Uncert 3(1), 011005 (Jul 30, 2018) (11 pages) Paper No: VVUQ-16-1030; doi: 10.1115/1.4040802 History: Received November 10, 2016; Revised July 05, 2018

## Abstract

The recently introduced basis adaptation method for homogeneous (Wiener) chaos expansions is explored in a new context where the rotation/projection matrices are computed by discovering the active subspace (AS) where the random input exhibits most of its variability. In the case where a One-dimensional (1D) AS exists, the methodology can be applicable to generalized polynomial chaos expansions (PCE), thus enabling the projection of a high-dimensional input to a single input variable and the efficient estimation of a univariate chaos expansion. Attractive features of this approach, such as the significant computational savings and the high accuracy in computing statistics of interest are investigated.

<>

## Introduction

While new technological advancements and the rapid increase of computational power simply seem to motivate the need for solving even more complex physical problems, uncertainty quantification tasks seem to suffer from the never-ending challenge of the “relatively” limited computational resources available to the experimenter. To overcome this challenge of developing problem- and computer code-independent efficient methodologies that will enable experimentation and predictive capabilities, one necessarily focuses on further advancing the mathematical and algorithmic tools for quantifying the various forms of uncertainties and their effects on physical models.

Most standard uncertainty quantification approaches involve the use of Monte Carlo (MC) methods [1] where several samples are drawn according to the model input distribution and their corresponding outputs are used to compute certain statistics of interest. The use of such methods typically varies from uncertainty propagation problems [2] to calibration problems [3,4] to stochastic optimization [5,6] and experimental design problems [7,8]. The slow convergence of the method along with the usually expensive computational models can easily make the approach unaffordable and infeasible, causing one to resort to cheaper alternatives. These consist of substituting the computational model with a surrogate that can be repeatedly evaluated at almost no cost, therefore accelerating computations. Such surrogates can be, for instance, polynomial chaos expansions (PCE) [912], Gaussian processes [13,14], or adaptive sparse grid collocation [15]. The common characteristic in all the above constructions is that an initial set of forward evaluations at preselected points is used to construct a functional approximation of the original computational model. After such an approximation becomes available, one no longer relies on the expensive-to-evaluate computer code but instead uses the surrogate. Although this characteristic generally outperforms MC methods, it can still fail to address the issue of computational efficiency due to the curse of dimensionality effect, that is the increase of the dimensionality of the uncertain input results is an exponentially increasing number of evaluations required to compute the surrogate parameters.

Several ways to achieve dimensionality reduction have been proposed in the literature, among which are sensitivity analysis methods [16] that rank the inputs according to their influence, thus allowing to neglect the components with respect to which, the output is insensitive. These methods involve, for example, variance decomposition methods such as the Sobol' indices [17,18]. Similar approaches might consist of spectral methods, such as the Karhunen–Loéve expansion (KLE) [9,19,20], for random vectors and random fields, where one expands in a series of scalar variables, the importance of which is determined by the eigenvalues of the covariance matrix. The series can be truncated accordingly, leaving out the unimportant scales of fluctuation of the random quantity. For an empirical version of the KLE, one might also consider the principal component analysis [21] or the kernel principal component analysis [22]. At last, one can proceed with the active subspace (AS) method [2325] of discovering the low-dimensional projection of the random input space where the model exhibits maximal variability. This makes use of a spectral decomposition of the covariance of the gradient vector, the eigenvectors of which reveal the projection mapping to the active subspace.

In this work, our attention is focused on reduction methods that are applicable to polynomial chaos expansions. PCEs provide an analytic and compact representation of the model output in terms of a number of scalar random variables, typically referred to as the germs, through polynomials that are pairwise orthogonal with respect to inner product (expectation) of the underlying Hilbert space of square integrable random variables. Their use, pioneered by Ghanem and Spanos [9], has been explored throughout several contexts and applications such as flow-through heterogeneous porous media [26,27], fluid dynamics [2830], and aeroelasticity [31] for either forward propagation [32,33] or inference problems [3436].

Most of the aforementioned dimensionality reduction methods can be applied to PCEs once the expansion becomes known. In Refs. [37] and [38], explicit formulas for the Sobol' sensitivity indices were derived with respect to the coefficients of chaos expansions, thus enabling the fast computation of variance decomposition factors whenever such expansions are available. The idea was used in the context of stochastic differential equations [39] to separate the influence of uncertainty due to random parameters from that caused by the driving white noise, and in a similar fashion, it has also been applied to large-eddy simulations [40] for sensitivity analysis due to turbulent effects. However, the key challenge of reducing the dimensionality of the input, in order to make the computation of the coefficients a feasible task, still remains unsolved and most of the works still focus their attention on simply obtaining sparse representations by adaptive finite elements [41], 1-minimization [42,43], and compressive sensing [44].

Recently [45] it was observed that for Hermite expansions with Gaussian input variables, the original input can be rotated using isometries onto the underlying Gaussian Hilbert space to obtain expansions with respect to the new basis that concentrate their dependence on only a few components. Certain choices for the rotation matrix were analytically seen not only to guarantee certain levels of sparsity in the expansion but even to concentrate the probabilistic behavior of scalar quantities of interest (QoI) on low-dimensional Gaussian subspaces. Computational schemes for efficiently computing the coefficients of the adapted expansion were also developed and the numerical results were impressive. The basis adaptation concept was further extended from scalar QoIs to random fields [46], where explicit formulas for computing the coefficients with respect to any rotation were derived and the case of a parametric family of rotations was discussed that gives rise to expansions with a Gaussian process input. Even though the idea of rotating the basis has been further applied to design optimization problems [47] and has been used to develop efficient schemes coupled with compressing sensing methods [48,49], it is still quite restricted to the homogeneous chaos expansions, where the Gaussian distribution remains invariant under rotations, a property that is not valid in other cases.

The key difference when applying the above idea in random vectors other than Gaussians is that the distribution or the rotated variables is, in general, not analytically available and the construction of polynomials, orthogonal with respect to this distribution, is not an easy task. In this paper, we are exploring the special case where an isometry can be used to rotate the basis in the same fashion as in the basis adaptation methodology, such that only an one-dimensional (1D) component of the new variables will suffice to build an accurate polynomial chaos expansion. In this quite restrictive case, the input can easily be mapped to a uniform distribution via its own cumulative distribution function, and therefore, the Legendre polynomials can be used to expand the series. To find the proper rotation matrix, we make use of the active subspace methodology that successfully discovers the rotation such that maximal variability can be achieved in one only variable. To do so, unlike the traditional active subspace method where the covariance matrix of the gradient vector is computed via MC methods, we make use of an initial chaos expansion that is obtained using low-order quadratures, that is with very few model evaluations. In this case, the matrix can be analytically computed at no additional cost. Once the active subspace is revealed, a new one-dimensional Legendre chaos expansion of relatively high order can be easily constructed as again it does not demand many evaluations. Apart from the fact that the approach enables surrogate construction of expensive models that would be otherwise infeasible, we also discover some attractive features: The capability of obtaining high order one-dimensional expansions allows for very accurate predictions of the probabilistic behavior of QoIs where a lower order full-dimensional expansion fails.

This paper is structured as follows: Sec. 2 reviews the polynomial chaos (PC) and active subspace approaches and then presents how the first can be used to analytically compute the latter along with an algorithm for efficient construction of 1D Legendre chaos expansions. In Sec. 3, we apply the methodology to two numerical examples. First, a simple quadratic polynomial function with analytically known active subspace is used for validation and second, a multiphase flow problem that simulates transport of ammonium and its oxidation to nitrite along a one-dimensional rod is used to draw my conclusions.

## Methodology

To provide a proper setting upon which our methodology is developed, we first consider a general response function $f:D⊂ℝd→ℝ$ where the points $ξ∈D$ will be thought of as inputs and their mappings $f(ξ)$ will the model outputs or the QoI. Typically, D is also equipped with some probability measure $P$; therefore, its elements are D-valued random variables from a space of events Ω to D. A σ-algebra $F$ that consists of the $P$-measurable subsets of Ω or in other words, the inverse mappings of the Borel subsets of D is naturally defined and thus the probability triplet $(Ω,F,P)$ is the space on which we will be working.

###### Generalized Polynomial Chaos.

Throughout this paper, we will assume that the quantity $f(ξ)$ has finite variance and therefore is a square integrable function, that is $f∈L2(Ω,F,P)$. It is known then [10] that f admits a series expansion in terms of orthogonal polynomials of $ξ$, given as Display Formula

(1)$f(ξ)=∑α,|α|=0∞fαψα(ξ)$

where $α=(α1,…,αd)∈J:=ℕ∪{0}$ are finite-dimensional multi-indices with norm $|α|=α1+⋯+αd$ and $ψα$ are d-dimensional polynomial functions of $ξ$ that are orthogonal with respect to the measure defined by the density function $p(ξ)$, that is Display Formula

(2)$E{ψα(ξ)ψβ(ξ)}=∫Dψα(ξ)ψβ(ξ)p(ξ)dξ=||ψα||2δα,β$

where $||ψα||=(E{ψα(ξ)2})1/2$ and $δα,β$ is the Dirac delta function which is 1 if $α=β$ and 0 otherwise. Without loss of generality, we assume here that $||ψα||=1$, that is the polynomials are normalized. We will refer to Eq. (1) as the PC expansion of f. In the case where $ξ=(ξ1,…,ξd)$ consists of independent and identically distributed random variables, the polynomials are expanded as Display Formula

(3)$ψα=∏i=1dψαi(ξi)$

where $ψαi(ξi)$ are univariate polynomials of ξi of order αi, i = 1,…, d. Common choices of the density function $p(ξ)=∏i=1dp(ξi)$ give rise to known forms of the polynomials, for instance the Gaussian, uniform, gamma, and beta distributions correspond to Hermite, Legendre, Laguerre, and Jacobi polynomials. respectively [10].

In practice, we work with truncated versions of Eq. (1), that is for $Q∈ℕ, JQ:={α∈J: |α|≤Q}$, we will assume that f can be accurately approximated by Display Formula

(4)$f(ξ)≈∑α∈JQfαψα(ξ)$

The above-mentioned expansion consists of Display Formula

(5)$NQ=(d+QQ)=(d+Q)!d!Q!$
basis terms, and the estimation of the corresponding coefficients is typically a challenging task. Several approaches have been developed in the literature for the estimation of the coefficients ${fα}α∈JQ$, divided into two main categories: Intrusive and nonintrusive methods. The first, pioneered by Ghanem and Spanos [9], treats the solution of a differential equation as a random field that can be written as in Eq. (4), where the coefficients vary as functions of the spatial or time parameters associated with the computational domain. The expression is then substituted in the equation satisfied by f, in order to derive the governing equations satisfied by the coefficients, which then need to be solved once. Nonintrusive approaches involve estimation of the coefficients using the relation Display Formula
(6)$fα=E{f(ξ)ψα(ξ)}$
which is computed using numerical integration techniques. As mentioned in Sec. 1, both approaches suffer by the curse of dimensionality which in the first case implies that the system of equations to be solved increases exponentially as a function of the dimensionality d and in the second case, the exponential increase is on the number of collocation points where the model f needs to be evaluated.

###### Discovering the Active Subspace.

Let us assume that the function f can be well approximated by Display Formula

(7)$f(ξ)≈g(WTξ)$

where W is a $d×d′$ matrix with orthonormal columns, that is it satisfies Display Formula

(8)$WTW=Id′$

and $g:ℝd′→ℝ$ is a $d′$-dimensional function that will be called the link function. Intuitively, such an assumption implies that the domain D can be rotated using a d × d orthonormal matrix Display Formula

(9)$A=[WV]$
with $ATA=Id$ and W, V being $d×d′$ and $d×(d−d′)$ matrices, respectively, such that f exhibits most of its variation on the space spanned by ${WTξ: ξ∈D}$ while it remains insensitive to variations on its orthogonal complement, thus motivating the term active subspace [23].

As it is clear from the above, the main challenge in the above construction is the determination of the rotation matrix A such that the active subspace defined by W will have the minimum possible dimensionality $d′$. To explore the directions along which f exhibits greatest sensitivity, we define the d × d matrix Display Formula

(10)$C=E{∇f(ξ)∇f(ξ)T}$

where $∇f(ξ)=((∂f(ξ)/∂ξ1),…,(∂f(ξ)/∂ξd))T$ is the gradient vector of f. This is essentially the uncentered covariance matrix of $∇f(ξ)$ and is a symmetric and positive semidefinite matrix; therefore, it admits a decomposition as Display Formula

(11)$C=UΛUT$

where $Λ$ is the diagonal matrix with entries the eigenvalues of C in decreasing order, namely λi, i = 1,…, d with λ1 ≥… ≥ λd and U is the unitary matrix whose ith column is the eigenvector of C corresponding to λi. By setting A = U and $η=ATξ$ it can be shown [23] that Display Formula

(12)$E{∇ηf)T(∇ηf)}=λ1+⋯+λd$

and by decomposing A as in Eq. (9) and setting $ηw=WTξ$ gives us that the mean-squared gradients of f with respect to $ηw$ will be the sum of the $d′$ largest eigenvalues. This construction provides us a way of rotating the input space D and separating it into two subspaces namely the active subspace and its orthogonal complement by simply selecting the $d′$ most dominant eigenvalues of C and their corresponding eigenvectors. Then the function f can be approximated as in Eq. (7).

###### Active Subspace Computation for Polynomial Chaos Expansions and Reduction to Univariate Legendre Chaos.

The key challenge in the active subspace methodology is the computation of the gradient matrixC, the standard way of which is by estimating the expectation using Monte Carlo sampling. Recently a gradient-free approach for determining the projection matrix when f is approximated by a Gaussian Process was developed [50]. For this purpose, assume that a PC expansion is available for f as in Eq. (4). Then, C can be written as Display Formula

(13)$C≈E{∇(∑fαψα(ξ))·∇(∑fβψβ(ξ))T}$

and the entries Cij will be given by Display Formula

(14)$Cij≈E{∂∂ξi(∑fαψα(ξ))·∂∂ξj(∑fβψβ(ξ))T}=E{(∑fα∂ψα(ξ)∂ξi)·(∑fβ∂ψβ(ξ)∂ξj)T}=∑α∑βfαfβE{∂ψα(ξ)∂ξi·∂ψβ(ξ)∂ξj}=fTKijf$

where Kij is the stiffness matrix with entries Display Formula

(15)$(Kij)αβ=E{∂ψα(ξ)∂ξi·∂ψβ(ξ)∂ξj}$

and $f={fα}α∈JQ$ is the vectorized representation of the chaos coefficients. The values of the entries of the stiffness matrices Kij, i, j = 1,…, d depend solely on the polynomials used in the PC expansion and their corresponding probability measures with respect to which the expectation is taken and can be computed independently of the nature of the function f. Detailed computation of Kij of the case of Legendre polynomials is provided in the Appendix. Computation of C generally requires full knowledge of the coefficients f and can be an expensive task. However, as we describe in Sec. 2.3, a relatively cheap estimation of the coefficients, based on low-level quadrature rules can suffice of this purpose.

Let us now assume that the matrix C is available and that λ1 ≫ λ2. This implies that the rotation matrix A can be decomposed in a one-dimensional vector w with wTw = 1 and a d × (d – 1) matrix V such that f exhibits most of its variation on the space spanned by $η=wTξ$ and can be accurately approximated as a function of η as in Eq. (7). The square integrability condition of f directly implies that f can be expanded in a polynomial chaos expansion in terms of polynomials that are orthogonal with respect to the probability measure of η. To avoid the structure of such polynomials, since the probability measure p(η) can be arbitrarily complex, we introduce the uniform $U(−1,1)$ germ Display Formula

(16)$ζ=2Fη(η)−1$

where Fη(⋅) is the cumulative probability distribution of η and we write Display Formula

(17)$f(ξ)≈g(η)=g(Fη−1(ζ+12)):=g(ζ)$

Finally g can be expanded as Display Formula

(18)$g(ζ)=∑n=0Nζgnψn(ζ)$

where ψn(ζ) are the normalized univariate Legendre polynomials and thus we can achieve a one-dimensional chaos decomposition of f.

###### Efficient Basis Reduction Using Pseudo-Spectral Projections.

Using Eq. (6), the coefficients of the chaos expansion of f can be estimated after approximating the integral with Display Formula

(19)$fα=∫Df(ξ)ψα(ξ)p(ξ)dξ≈∑i=1sf(ξ(i))ψα(ξ(i))wi, α∈JQ$
using a quadrature rule ${ξ(i),wi}i=1sd$, where ${ξ(i)}i=1sd$ are quadrature points and ${wi}i=1sd$ are the corresponding weights. As mentioned previously, such a procedure can be prohibitive for relatively large d or for computationally expensive models f that exhibit high nonlinearity. It is desirable, in such cases, to develop a computational strategy that reduces the computational resources by limiting the number of model evaluations. Provided that a one-dimensional active subspace exists, writing f in the form Eq. (18) and computing the coefficients ${gn}n=0Nζ$ would be an efficient alternative. Since this requires the knowledge of the projection vector w, one could perform the following steps: First, a low-level quadrature rule that consists of an affordable number of points can be performed in order to compute the low-order coefficients ${fα}α∈JQ0, Q0 and an estimate of C and therefore of w can be obtained. Then, the PC coefficients of $g$ can be computed in a similar manner by Display Formula
(20)$gn≈∑i=1s1g(ζ(i))ψn(ζ(i))wi$

where ${ζ(i)}i=1s1$ are one-dimensional quadrature points on the [–1,1] space. In practice, the evaluation of $g(ζ(i))$ will require transforming ζ(i) to $η(i)=Fη−1(ζ(i)+1/2)$ and $η(i)=Aξ(i)$ where $η(i)$ is a d-dimensional arbitrary completion of η(i), where the remaining d – 1 entries will essentially play no role in the model output as they span the subspace with respect to which f is approximately invariant. This procedure is summarized in the pseudo-algorithm described in Algorithm 1.

For the implementation of the above-mentioned algorithm in the numerical examples that are presented in this paper, the python package chaos_basispy [51] has been used, which is equipped with polynomial chaos basis function evaluations capabilities, quadrature points generation, and computation of the stiffness matrices Kij derived above.

###### Error Analysis.

The procedure described so far includes computation of the coefficients of a low-order polynomial chaos expansion using an efficient numerical integration scheme, typically a low-level sparse quadrature rule. Although the estimation of the coefficients can be satisfactorily accurate, the computation of the gradient vector covariance matrix and subsequently the computation of the active subspace are still subject to the truncation error introduced by using a low-order chaos expansion. Below we attempt to quantify this error be deriving some rather standard upper bounds for error of the covariance estimate and its eigenvalues.

Let us denote with $f̂(ξ)$ the truncated version of $f(ξ)$ to be used for computation of CDisplay Formula

(21)$f̂(ξ)=∑α∈JQ0fαψα(ξ)$

where Q0 < Q and write Display Formula

(22)$f(ξ)=f̂(ξ)+ε(ξ)$

The error of approximating the gradient can then be written as Display Formula

(23)$∇ξε(ξ)=∇ξf(ξ)−∇ξf̂(ξ)$

By defining Display Formula

(24)$γQ0:=||∇ξε(ξ)||$

where $||·||$ is the Euclidean norm, and we can write Display Formula

(25)$||∇ξf̂(ξ)−∇ξf(ξ)||≤γQ0$

where clearly $γQ0→0$ as Q0Q. Define also Display Formula

(26)$Ĉ=E{∇f̂(ξ)∇f̂(ξ)T}$

By referring to the spectral norm (induced by the Euclidean norm) when $||·||$ is applied to matrices, the following Theorem can be stated:

Theorem 1. The norm of the difference between C and$Ĉ$is bounded by

$||Ĉ−C||≤E[γQ02]+2E[LγQ0]$
where$L=||∇ξf(ξ)||$.

Proof. First note that Display Formula

(27)$||∇ξf̂+∇ξf||=||∇ξf̂−∇ξf+2∇ξf||≤||∇ξf̂−∇ξf||+2||∇ξf||≤γQ0+2L$

Then we have Display Formula

(28)$||∇ξf̂∇ξf̂T−∇ξf∇ξfT||=12||(∇f̂−∇f)(∇f̂+∇f)T+(∇f̂+∇f)(∇f̂−∇f)T||≤||(∇f̂+∇f)(∇f̂−∇f)T||≤(γQ0+2L)γQ0$

Finally we get Display Formula

(29)$||Ĉ−C||=||E[∇ξf̂∇ξf̂T]−E[∇ξf∇ξfT]||=||E[∇ξf̂∇ξf̂T−∇ξf∇ξfT]||≤E[||∇ξf̂∇ξf̂T−∇ξf∇ξfT||]≤E[γQ02]+2E[LγQ0] ◻$

The following corollary also holds:

Corollary. The difference between the kth true eigenvalue λk and the corresponding estimate θk is also bounded as

$|λk−θk|≤E[γQ02]+2E[LγQ0]$

Proof. Simply observe that

$|λk−θk|≤||Ĉ−C||≤E[γQ02]+2E[LγQ0]$

where the first inequality follows from Corollary 8.1.6 in Ref. [52] and the second from Theorem 1.

In the above expressions for the error bounds, we can further write Display Formula

(30)$γQ0={∑i=1d(∑|α|=Q0+1Qfα∂ψα(ξ)∂ξi)2}1/2$
which gives Display Formula
(31)$E[γQ02]=∑i=1d∑|α|=Q0+1Q∑|β|=Q0+1QfαfβE[∂ψα(ξ)∂ξi∂ψβ(ξ)∂ξi]$
Display Formula
(32)$=∑i=1d∑|α|,|β|=Q0+1Qfαfβ(Kij)αβ$
Display Formula
(33)$=∑i=1d∑|α|,|β|=Q0+1α−i=β−iQfαfβ(∏j≠ij=1dE[ψ′αj(ξj)ψβj(ξj)])$

In the above, $α−i, β−i$ indicate the multi-indices $α, β$ where the entries αi, βi, respectively, are excluded and the last line follows by using the expressions for $(Kij)αβ$, derived in Appendix. For the second term of the error bound, one can simply use the Cauchy–Schwarz inequality to write $E[LγQ0]≤E[L2]1/2E[γQ02]1/2$ and derive a similar expression for $E[L2]$, as above. Further numerical investigation of the behavior of the above error bound falls beyond the scope of this work.

## Numerical Examples

###### Polynomial Function.

Consider the function

$f:ℝd→ℝ$

with Display Formula

(34)$f(ξ)=a+bwTξ+cξTwwTξ$

with constants a, b, c, and $w∈ℝd$ such that $wTw=1$. For the variables, we assume that $ξ=(ξ1,…,ξd)$ with ξi independent $U(−1,1)$. This gives Display Formula

(35)$∇f(x)=w(b+2cwTξ)$
which gives Display Formula
(36)$E[∇f(ξ)∇f(ξ)T]=wE[(b+2cwTξ)(b+2cwTξ)]wT=w(b2+4c2E[ξTwwTξ])wT=w(b2+43c2)wT$

and it is clear that $b2+(4/3)c2$ is the only nonzero eigenvalue corresponding to the eigenvector w, giving a one-dimensional active subspace. Setting $η=wTξ$ we get Display Formula

(37)$f(ξ)=g(wTξ)$

where the link function g(η) is Display Formula

(38)$g(η)=a+bη+cη2$

Letting Fη be the cumulative distribution function of η, we are interested in constructing a PC expansion with respect to ζ = (2Fη(η) – 1). In our numerical implementation, we use arbitrary values for a, b, c, and w that were randomly generated (and are easily reproducible by fixing the random seed) and we take d = 10. We have Display Formula

(39)$w=(0.1404,−0.3574,0.4267,−0.0931,−0.21460.2642,0.2560,−0.1895,0.0046,−0.6680)T$

(40)$a=1.1500,b=0.9919,c=0.9533$

Figure 1 (left) shows the empirical pdf of η for the above choice of w and Fig. 1 (right) shows the quadrature points of a level-5 Clenshaw–Curtis rule and their transformed values on the η-space through the mapping $Fη−1(ζ+1/2)$. Note here that the expression (34) can be rearranged as a series of Legendre polynomials; therefore, its PC expansion is essentially known and has exact order 2. The same expansion can be computed numerically using sparse quadrature rule and specifically, a level-2 rule suffices for an accurate estimation, corresponding to 221 evaluations. Estimation of the coefficients of a 1d PC expansion with respect to ζ using the known orthogonal transformation w can be achieved using a level-5 rule that corresponds to 33 evaluations of f. Note that the PC expansion with respect to ζ is no longer a second order series. The order of polynomials necessary to achieve convergence here depends on the nature of the inverse cdf $Fη−1$. In our case, we find that a 20th order polynomial chaos expansion suffices to achieve a low error. Figure 2 shows the resulting expansions evaluated at 1000 sample points $ξ$ or their corresponding $ζ=2Fη(wTξ)−1$ accordingly and the density function of the output. Note that in order to obtain a common input domain, the plots show the dependence of the expansions with respect to the η sample points.

###### Multiphase Flow: One-Dimensional Transport and Decay of Ammonium.

Multiphase-multicomponent flow and transport models are known for their particular challenges associated with the uncertainty propagation problem which are mainly due to the highly nonlinear nature of the coupled systems of equations that describe the complex physical process and the varying role of the many random parameters involved in the system. Here we consider the problem of a one-dimensional transport of the microbially mediated first-order oxidation of ammonium $(NH4+)$ to nitrite $(NO2−)$ and the subsequent oxidation of nitrite to its ion $(NO3−)$ that has been previously investigated in Refs. [5355]. The flow domain under consideration where the flow of ammonium takes places is taken here to be a one-dimensional column of 2 m length.

For the numerical solution of the transport and decay model, we employ the multiphase-multicomponent simulator TOUGH2 [56] that solves the integral form of the system of governing equations using the integral finite difference method. More specifically, the EOS7r module [57] that specializes in modeling radionuclide transport is used here, allowing the description of the oxidation effect in a similar manner as specifying half-life parameters to describe radioactive or biodegradation effects. Detailed description of the balance equations and capillary pressure model can be found in Ref. [56]. The mass components considered here are water, air, ammonium, and nitrite, and the fluid phases are aqueous and gas. We consider a first-order decay law for the mass components [57].

In the scenario considered here, the domain is discretized into 400 blocks, each of length 0.005 m. and a source of ammonium is placed on the one boundary (y = 0), while no concentration of nitrite is initially present. The domain is shown in Fig. 3. Uncertainty is introduced in the model through six model parameters, namely the decay parameters of ammonium and nitrite, $T1/2am$ and $T1/2n$, respectively, the constant flux of ammonium at the boundary or more precisely the constant injected concentration $χ0am$, the distribution coefficient $Kdwam$ of ammonium that affects its adsorption onto the immobile solid grains, and the porosity ϕ of the medium. All other model parameters were assigned fixed values which are shown in Table 1. All random input parameters, denoted in a vector form as $θ=(T1/2am,T1/2am,χ0am,Kdwam,ϕ)T$ are assumed to follow uniform distributions $θi∼U(θli,θui)$ where the intervals $[θli,θui], i=1,…,5$ are given in Table 2. Next, for the construction of a polynomial chaos approximation of model output, the random input parameters are rescaled to $U(−1,1)$ random variables via the transformation Display Formula

(41)$ξi=2θi−θliθui−θli−1$

###### Results.

The transport model runs for a time period corresponding to t = 200 h and the mass fraction of ammonium in the aqueous phase is observed. The quantity of interest considered here is the mass fraction $χtam(y)$ at the middle point of the column (y = 1 m distance from the source), where neglect the index w for simplicity. First, to illustrate the motivation for applying a basis reduction procedure on the QoI, we compute the polynomial chaos expansion of the mass fraction of ammonium over the whole domain y ∈ [0, 2]. Figure 4 shows a comparison of 10 Monte Carlo outputs of the model with the polynomials chaos expansions evaluated at the same MC inputs. The four figures correspond to the cases where the coefficients have been estimated using a Clenshaw–Curtis quadrature rule of level 2, 3, 4, and 5, respectively. Qualitatively we observe that for level 3 and higher, the surrogate is a very good approximation of the model. However, at y = 1 we can see that the chaos expansions of $χtam(y)$ fall sometimes below zero which is not an acceptable value since $χt(y)am≥0$. This immediately makes the chaos expansion inaccurate not only as a surrogate for fast computation of the model output but will also result in false empirical distributions of the QoI. Next, we use the estimated coefficients for estimation of the gradient matrix G and we denote with G the gradient matrix estimated using the coefficients from level quadrature rule. The eigenvalues and the dominating eigenvector w are shown in Fig. 5. For the cases  = 3, 4, and 5, the eigenvalues and the vector w almost coincide, indicating convergence of the eigenvalues and existence of a one-dimensional active subspace. In the  = 2 case, the eigenvalues seem to be inaccurate and perhaps the projection vector does not define an optimal rotation for dimension reduction. In addition, even though a “gap” can be seen between the first and the remaining eigenvalues, their values indicate that they are not negligible. Figure 6 (left) shows the 56 estimated coefficients of the expansions for the different quadrature rule levels, and similar conclusions can be drawn regarding the convergence of the numerical integration and the accuracy of the expansion, as the level-2 results are completely “off” comparing to the remaining cases and even several coefficients of the level-3 case display some discrepancy before being stabilized in level-4 and level-5 rules.

After defining the random active variables $ηℓ=wℓTξ$ with cdf $Fηℓ(·)$ corresponding to the active subspaces revealed by the chaos expansions estimated above, we construct the PC approximations with respect to the germs Display Formula

(42)$ζℓ=2Fηℓ(ηℓ)−1, ℓ=2,…,5$
of order 15 using the level-5 one-dimensional CC quadrature rule. The estimated coefficients are displayed in Fig. 6 (right), where except the one corresponding to the level 2 rule, they seem to be in good agreement, as expected from the similarity of their respective projection vectors w. Figure 7 shows values of 1000 MC outputs of the true model and the 1D PC expansion as a function the common input η. Evaluation of the full model is carried out at the points $ξ(i)=wℓηℓ(i)$ while evaluation of the 1D PC expansions is obtained at $ζ(i)=Fηℓ(ηℓ(i))$. It is observed again that for ≥ 3, the one-dimensional representation of $χtam$ provides a quite accurate approximation of the true QoI whose scatter around the 1D curve due to the influence of the orthogonal subspace $span{wTξ}T$ is very small.

Finally, we compute the empirical pdfs based on 105 MC samples of all 5D and 1D PC expansions that we have and compare their histograms with that of the 1000 MC samples that are available directly from the TOUGH2 model. The results are shown in Fig. 8. It is observed again that those based on the 1D expansions provide far better approximations of the true histogram which is by definition bounded from below at 0 while those based on the full expansions fail dramatically to capture the lower tail behavior. The computational resources required for the estimation of the expansions described are also in favor of the basis reduction methodology: The 1D expansion based on an orthogonal projection w that was computed using  = 3 quadrature rule required a total of 241 model evaluations for the estimation of the third-order 5D PC plus 33 model evaluations for a level-5 1D quadrature rule resulting in a total of 274 model evaluations and eventually provides a more accurate density function than a  = 4, or 5 rule for the full PC that require 801 and 2433 model evaluations, respectively!

## Conclusions

I have presented a methodology that combines the active subspace approach for dimensionality reduction with the construction of polynomial chaos surrogates for the purpose of efficient exploration of response surfaces and uncertainty propagation problems. This approach, serving at the same time as an extension of the basis adaptation framework [45] for homogeneous chaos expansions, due to its applicability to generalized polynomial chaos, provides an explicit expression of the gradient covariance matrix in terms of the chaos coefficients that allows fast computation of the active subspace that is not only useful for dimension reduction, but even for sensitivity analysis purposes. Besides the attractive features in terms of computational efficiency and improved accuracy of the surrogate that were mentioned throughout the paper and demonstrated in the numerical examples, the method shows the direction to future challenges: Precisely, the extension of the current approach to higher-dimensional active subspaces, either by finding a way to map the projected variables to a vector of independent ones or by constructing orthogonal polynomials with respect to arbitrary joint distributions, faces several obstacles as detailed in Ref. [58] but it would provide a new ground on polynomial chaos-based dimensionality reduction methods.

## Appendices

###### Appendix: Computation of the Stiffness Matrix Kij

First write the partial derivatives Display Formula

(A1)$∂ψα(ξ)∂ξi=∂∂ξi∏m=1dψαm(ξm)=ψ′αi(ξi)∏m=1m≠idψαm(ξm)$

Then for any $α,β∈JP$, we have Display Formula

(A2)$(Kij)αβ=E{∂ψα(ξ)∂ξi·∂ψβ(ξ)∂ξj}=E{(ψ′αi(ξi)∏m=1m≠idψαm(ξm))·(ψ′βj(ξj)∏m=1m≠jdψβm(ξm))}$
which for i ≠ j gives Display Formula
(A3)$(Kij)αβ=(∏m=1m≠i,m≠jdE{ψαm(ξm)ψβm(ξm)})·E{ψ′αi(ξi)ψβi(ξi)}·E{ψαj(ξj)ψ′βj(ξj)}=(∏m=1m≠i,m≠jdδαm,βm)·E{ψ′αi(ξi)ψβi(ξi)}·E{ψαj(ξj)ψ′βj(ξj)}$

and for i = j gives Display Formula

(A4)$(Kii)αβ=(∏m=1m≠idE{ψαm(ξm)ψβm(ξm)})·E{ψ′αi(ξi)ψ′βi(ξi)}=(∏m=1m≠idδαm,βm)·E{ψ′αi(ξi)ψ′βi(ξi)}.$

The computation of the terms $E{ψ′α(ξ)ψβ(ξ)}$ and $E{ψ′α(ξ)ψ′β(ξ)}$ depends on the choice of polynomials. The case of Hermite polynomials was derived in Ref. [48]. Here I explore the case of Legendre polynomials.

In the case where $ψα(ξ)=(ℓα(ξ)/||ℓα(ξ)||)$ are the normalized Legendre polynomials, the expectations in the last two terms of the above relation can be computed recursively using the relation Display Formula

(A5)$(2α+1)ℓα(ξ)=ddξ[ℓα+1(ξ)−ℓα−1(ξ)]$
that gives Display Formula
(A6)$E{ℓ′α(ξ)ℓβ(ξ)}=δα−1,β+E{ℓ′α−2(ξ)ℓβ(ξ)}=δα−1,β+δα−3,β+δα−5,β+⋯$

(A7)$E{ℓ′α(ξ)ℓ′β(ξ)}=(2α−1)δα−1,β−1+(2α−1)E{ℓα−1(ξ)ℓ′β−2(ξ)}+(2β−1)E{ℓ′α−2(ξ)ℓβ−1(ξ)}+E{ℓ′α−2(ξ)ℓ′β−2(ξ)}=(2α−1)(δα−1,β−1+δα−1,β−3+⋯)+(2α−5)(δα−3,β−1+δα−3,β−3+⋯)+(2α−9)(δα−5,β−1+δα−5,β−3+⋯)+⋯.$

Using Eqs. (A6) and (A7), finally I compute Display Formula

(A8)$E{ψ′α(ξ)ψβ(ξ)}=(2α+1)(2β+1)E{ℓ′α(ξ)ℓβ(ξ)}$

(A9)$E{ψ′α(ξ)ψ′β(ξ)}=(2α+1)(2β+1)E{ℓ′α(ξ)ℓ′β(ξ)}$

## References

Robert, C. , and Casella, G. , 2013, Monte Carlo Statistical Methods, Springer Science & Business Media, New York.
Morokoff, W. , and Caflisch, R. , 1995, “ Quasi-Monte Carlo Integration,” J. Comput. Phys., 122(2), pp. 218–230.
Tarantola, A. , 2005, Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics, Philadelphia, PA.
Mosegaard, K. , and Tarantola, A. , 1995, “ Monte Carlo Sampling of Solutions to Inverse Problems,” J. Geophys. Res.: Solid Earth, 100(B7), pp. 12431–12447.
Spall, J. C. , 2005, Introduction to Stochastic Search and Optimization: Estimation, Simulation, and Control, Vol. 65, Wiley, Hoboken, NJ.
Spall, J. , 1992, “ Multivariate Stochastic Approximation Using a Simultaneous Perturbation Gradient Approximation,” IEEE Trans. Autom. Control, 37(3), pp. 332–341.
Huan, X. , and Marzouk, Y. , 2013, “ Simulation-Based Optimal Bayesian Experimental Design for Nonlinear Systems,” J. Comput. Phys., 232(1), pp. 288–317.
Tsilifis, P. , Ghanem, R. , and Hajali, P. , 2017, “ Efficient Bayesian Experimentation Using an Expected Information Gain Lower Bound,” SIAM/ASA J. Uncertainty Quantif., 5(1), pp. 30–62.
Ghanem, R. , and Spanos, P. , 2012, Stochastic Finite Elements: A Spectral Approach, Revised Edition, Dover Publications, Mineola, NY.
Xiu, D. , and Karniadakis, G. , 2002, “ The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations,” SIAM J. Sci. Comput., 24(2), pp. 619–644.
Babuska, I. , Nobile, F. , and Tempone, R. , 2007, “ A Stochastic Collocation Method for Elliptic Partial Differential Equations With Random Input Data,” SIAM J. Numer. Anal., 45(3), pp. 1005–1034.
Reagan, M. , Najm, H. , Ghanem, R. , and Knio, O. , 2003, “ Uncertainty Quantification in Reacting-Flow Simulations Through Non-Intrusive Spectral Projection,” Combust. Flame, 132(3), pp. 545–555.
Bilionis, I. , and Zabaras, N. , 2012, “ Multi-Output Local Gaussian Process Regression: Applications to Uncertainty Quantification,” J. Comput. Phys., 231(17), pp. 5718–5746.
Bilionis, I. , Zabaras, N. , Konomi, B. , and Lin, G. , 2013, “ Multi-Output Separable Gaussian Process: Towards an Efficient, Fully Bayesian Paradigm for Uncertainty Quantification,” J. Comput. Phys., 241, pp. 212–239.
Ma, X. , and Zabaras, N. , 2009, “ An Adaptive Hierarchical Sparse Grid Collocation Algorithm for the Solution of Stochastic Differential Equations,” J. Comput. Phys., 228(8), pp. 3084–3113.
Saltelli, A. , Ratto, M. , Andres, T. , Campolongo, F. , Cariboni, J. , Gatelli, D. , Saisana, M. , and Tarantola, S. , 2008, Global Sensitivity Analysis: The Primer, Wiley, Chichester, UK.
Sobol', I. , 1990, “ On Sensitivity Estimation for Nonlinear Mathematical Models,” Mat. Model., 2, pp. 112–118.
Owen, A. , 2013, “ Variance Components and Generalized Sobol' Indices,” SIAM/ASA J. Uncertainty Quantif., 1(1), pp. 19–41.
Karhunen, K. , 1946, “ Über Lineare Methoden in Der Wahrscheinlichkeits-Rechnung,” Ann. Acad. Sci. Fennicade Ser. A1, Math. Phys., 37, pp. 3–79.
Loéve, M. , 1955, Probability Theory, D. Van Nostrand, Princeton, NJ.
Pearson, K. , 1901, “ Liii. on Lines and Planes of Closest Fit to Systems of Points in Space,” London, Edinburgh, Dublin Philos. Mag. J. Sci., 2(11), pp. 559–572.
Ma, X. , and Zabaras, N. , 2011, “ Kernel Principal Component Analysis for Stochastic Input Model Generation,” J. Comput. Phys., 230(19), pp. 7311–7331.
Constantine, P. , Dow, E. , and Wang, Q. , 2014, “ Active Subspace Methods in Theory and Practice: Applications to Kriging Surfaces,” SIAM J. Sci. Comput., 36(4), pp. A1500–A1524.
Constantine, P. , Emory, M. , Larsson, J. , and Iaccarino, G. , 2015, “ Exploiting Active Subspaces to Quantify Uncertainty in the Numerical Simulation of the Hyshot Ii Scramjet,” J. Comput. Phys., 302, pp. 1–20.
Lukaczyk, T. , Palacios, F. , Alonso, J. , and Constantine, P. , 2014, “ Active Subspaces for Shape Optimization,” AIAA Paper No. AIAA 2014-1171.
Ghanem, R. , 1998, “ Scales of Fluctuation and the Propagation of Uncertainty in Random Porous Media,” Water Resour. Res., 34(9), pp. 2123–2136.
Ghanem, R. , and Dham, S. , 1998, “ Stochastic Finite Element Analysis for Multiphase Flow in Heterogeneous Porous Media,” Transp. Porous Media, 32(3), pp. 239–262.
Najm, H. , 2009, “ Uncertainty Quantification and Polynomial Chaos Techniques in Computational Fluid Dynamics,” Annu. Rev. Fluid Mech., 41(1), pp. 35–52.
Le Maître, O. , Reagan, M. , Najm, H. , Ghanem, R. , and Knio, O. , 2002, “ A Stochastic Projection Method for Fluid Flow—II: Random Process,” J. Comput. Phys., 181(1), pp. 9–44.
Xiu, D. , and Karniadakis, G. , 2003, “ Modeling Uncertainty in Flow Simulations Via Generalized Polynomial Chaos,” J. Comput. Phys., 187(1), pp. 137–167.
Arnst, M. , Ghanem, R. , and Soize, C. , 2010, “ Identification of Bayesian Posteriors for Coefficients of Chaos Expansions,” J. Comput. Phys., 229(9), pp. 3134–3154.
Ghanem, R. , Doostan, A. , and Red-Horse, J. , 2008, “ A Probabilistic Construction of Model Validation,” Comput. Methods Appl. Mech. Eng., 197(29–32), pp. 2585–2595.
Ghanem, R. , and Red-Horse, J. , 1999, “ Propagation of Probabilistic Uncertainty in Complex Physical Systems Using a Stochastic Finite Element Approach,” Phys. D: Nonlinear Phenom., 133(1–4), pp. 137–144.
Marzouk, Y. M. , Najm, H. N. , and Rahn, L. , 2007, “ Stochastic Spectral Methods for Efficient Bayesian Solution of Inverse Problems,” J. Comput. Phys., 224(2), pp. 560–586.
Marzouk, Y. , and Najm, H. , 2009, “ Dimensionality Reduction and Polynomial Chaos Acceleration of Bayesian Inference in Inverse Problems,” J. Comput. Phys., 228(6), pp. 1862–1902.
Ghanem, R. , and Doostan, R. , 2006, “ Characterization of Stochastic System Parameters From Experimental Data: A Bayesian Inference Approach,” J. Comput. Phys., 217(1), pp. 63–81.
Sudret, B. , 2007, “ Global Sensitivity Analysis Using Polynomial Chaos Expansions,” Reliab. Eng. Syst. Saf., 93(7), pp. 964–979.
Crestaux, T. , Le Maître, O. , and Martinez, J. , 2009, “ Polynomial Chaos Expansion for Sensitivity Analysis,” Reliab. Eng. Syst. Saf., 94(7), pp. 1161–1172.
Le Maître, O. , and Knio, O. , 2015, “ Pc Analysis of Stochastic Differential Equations Driven by Wiener Noise,” Reliab. Eng. Syst. Saf., 135, pp. 107–124.
Lucor, D. , Meyers, J. , and Sagaut, P. , 2007, “ Sensitivity Analysis of Large-Eddy Simulations to Subgrid-Scale-Model Parametric Uncertainty Using Polynomial Chaos,” J. Fluid Mech., 585, pp. 255–279.
Blatman, G. , and Sudret, B. , 2008, “ Sparse Polynomial Chaos Expansions and Adaptive Stochastic Finite Elements Using a Regression Approach,” C. R. Méc., 336(6), pp. 518–523.
Peng, J. , Hampton, J. , and Doostan, A. , 2014, “ A Weighted ℓ1-Minimization Approach for Sparse Polynomial Chaos Expansions,” J. Comput. Phys., 267, pp. 92–111.
Yang, X. , and Karniadakis, G. , 2013, “ Reweighted ℓ1 Minimization Method for Stochastic Elliptic Differential Equations,” J. Comput. Phys., 248, pp. 87–108.
Hampton, J. , and Doostan, A. , 2015, “ Compressive Sampling of Polynomial Chaos Expansions: Convergence Analysis and Sampling Strategies,” J. Comput. Phys., 280, pp. 363–386.
Tipireddy, R. , and Ghanem, R. , 2014, “ Basis Adaptation in Homogeneous Chaos Spaces,” J. Comput. Phys., 259, pp. 304–317.
Tsilifis, P. , and Ghanem, R. , 2017, “ Reduced Wiener Chaos Representation of Random Fields Via Basis Adaptation and Projection,” J. Comput. Phys., 341, pp. 102–120.
Thimmisetty, C. , Tsilifis, P. , and Ghanem, R. , 2017, “ Homogeneous Chaos Basis Adaptation for Design Optimization Under Uncertainty: Application to the Oil Well Placement Problem,” Ai Edam, 31(3), pp. 265–276.
Yang, X. , Lei, H. , Baker, N. , and Lin, G. , 2016, “ Enhancing Sparsity of Hermite Polynomial Expansions by Iterative Rotations,” J. Comput. Phys., 307, pp. 94–109.
Tsilifis, P. , Huan, X. , Safta, C. , Sargsyan, K. , Lacaze, G. , Oefelein, J. , Najm, H. , and Ghanem, R. , 2018, “ Compressive Sensing Adaptation for Polynomial Chaos Expansions,” Preprint .
Tripathy, R. , Bilionis, I. , and Gonzalez, M. , 2016, “ Gaussian Processes With Built-in Dimensionality Reduction: Applications to High-Dimensional Uncertainty Propagation,” J. Comput. Phys., 321, pp. 191–223.
Tsilifis, P. , 2017, “ A Polynomial Chaos Basis Reduction Framework in Python,” PyPI: Python Package Index, accessed July 17, 2018,
Golub, G. , and Van Loan, C. , 2012, Matrix Computations, Vol. 3, JHU Press, Baltimore, MD.
Cho, C. , 1971, “ Convective Transport of Ammonium With Nitrification in Soil,” Can. J. Soil Sci., 51(3), pp. 339–350.
McNab, W. , and Narasimhan, T. , 1993, “ A Multiple Species Transport Model With Sequential Decay Chain Interactions in Heterogeneous Subsurface Environments,” Water Resour. Res., 29(8), pp. 2737–2746.
Van Genuchten, M. , 1985, “ Convective-Dispersive Transport of Solutes Involved in Sequential First-Order Decay Reactions,” Comput. Geosci., 11(2), pp. 129–147.
Pruess, K. , Oldenburg, C. , and Moridis, G. , 1999, “ TOUGH2 User's Guide, Version 2,” Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, Report No. LBNL-43134.
Oldenburg, C. , and Pruess, K. , 1995, “ EOS7R: Radionuclide Transport for TOUGH2,” Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, Report No. LBL-34868.
Soize, C. , and Desceliers, C. , 2010, “ Computational Aspects for Constructing Realizations of Polynomial Chaos in High Dimension,” SIAM J. Sci. Comput., 32(5), pp. 2820–2831.
View article in PDF format.

## References

Robert, C. , and Casella, G. , 2013, Monte Carlo Statistical Methods, Springer Science & Business Media, New York.
Morokoff, W. , and Caflisch, R. , 1995, “ Quasi-Monte Carlo Integration,” J. Comput. Phys., 122(2), pp. 218–230.
Tarantola, A. , 2005, Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics, Philadelphia, PA.
Mosegaard, K. , and Tarantola, A. , 1995, “ Monte Carlo Sampling of Solutions to Inverse Problems,” J. Geophys. Res.: Solid Earth, 100(B7), pp. 12431–12447.
Spall, J. C. , 2005, Introduction to Stochastic Search and Optimization: Estimation, Simulation, and Control, Vol. 65, Wiley, Hoboken, NJ.
Spall, J. , 1992, “ Multivariate Stochastic Approximation Using a Simultaneous Perturbation Gradient Approximation,” IEEE Trans. Autom. Control, 37(3), pp. 332–341.
Huan, X. , and Marzouk, Y. , 2013, “ Simulation-Based Optimal Bayesian Experimental Design for Nonlinear Systems,” J. Comput. Phys., 232(1), pp. 288–317.
Tsilifis, P. , Ghanem, R. , and Hajali, P. , 2017, “ Efficient Bayesian Experimentation Using an Expected Information Gain Lower Bound,” SIAM/ASA J. Uncertainty Quantif., 5(1), pp. 30–62.
Ghanem, R. , and Spanos, P. , 2012, Stochastic Finite Elements: A Spectral Approach, Revised Edition, Dover Publications, Mineola, NY.
Xiu, D. , and Karniadakis, G. , 2002, “ The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations,” SIAM J. Sci. Comput., 24(2), pp. 619–644.
Babuska, I. , Nobile, F. , and Tempone, R. , 2007, “ A Stochastic Collocation Method for Elliptic Partial Differential Equations With Random Input Data,” SIAM J. Numer. Anal., 45(3), pp. 1005–1034.
Reagan, M. , Najm, H. , Ghanem, R. , and Knio, O. , 2003, “ Uncertainty Quantification in Reacting-Flow Simulations Through Non-Intrusive Spectral Projection,” Combust. Flame, 132(3), pp. 545–555.
Bilionis, I. , and Zabaras, N. , 2012, “ Multi-Output Local Gaussian Process Regression: Applications to Uncertainty Quantification,” J. Comput. Phys., 231(17), pp. 5718–5746.
Bilionis, I. , Zabaras, N. , Konomi, B. , and Lin, G. , 2013, “ Multi-Output Separable Gaussian Process: Towards an Efficient, Fully Bayesian Paradigm for Uncertainty Quantification,” J. Comput. Phys., 241, pp. 212–239.
Ma, X. , and Zabaras, N. , 2009, “ An Adaptive Hierarchical Sparse Grid Collocation Algorithm for the Solution of Stochastic Differential Equations,” J. Comput. Phys., 228(8), pp. 3084–3113.
Saltelli, A. , Ratto, M. , Andres, T. , Campolongo, F. , Cariboni, J. , Gatelli, D. , Saisana, M. , and Tarantola, S. , 2008, Global Sensitivity Analysis: The Primer, Wiley, Chichester, UK.
Sobol', I. , 1990, “ On Sensitivity Estimation for Nonlinear Mathematical Models,” Mat. Model., 2, pp. 112–118.
Owen, A. , 2013, “ Variance Components and Generalized Sobol' Indices,” SIAM/ASA J. Uncertainty Quantif., 1(1), pp. 19–41.
Karhunen, K. , 1946, “ Über Lineare Methoden in Der Wahrscheinlichkeits-Rechnung,” Ann. Acad. Sci. Fennicade Ser. A1, Math. Phys., 37, pp. 3–79.
Loéve, M. , 1955, Probability Theory, D. Van Nostrand, Princeton, NJ.
Pearson, K. , 1901, “ Liii. on Lines and Planes of Closest Fit to Systems of Points in Space,” London, Edinburgh, Dublin Philos. Mag. J. Sci., 2(11), pp. 559–572.
Ma, X. , and Zabaras, N. , 2011, “ Kernel Principal Component Analysis for Stochastic Input Model Generation,” J. Comput. Phys., 230(19), pp. 7311–7331.
Constantine, P. , Dow, E. , and Wang, Q. , 2014, “ Active Subspace Methods in Theory and Practice: Applications to Kriging Surfaces,” SIAM J. Sci. Comput., 36(4), pp. A1500–A1524.
Constantine, P. , Emory, M. , Larsson, J. , and Iaccarino, G. , 2015, “ Exploiting Active Subspaces to Quantify Uncertainty in the Numerical Simulation of the Hyshot Ii Scramjet,” J. Comput. Phys., 302, pp. 1–20.
Lukaczyk, T. , Palacios, F. , Alonso, J. , and Constantine, P. , 2014, “ Active Subspaces for Shape Optimization,” AIAA Paper No. AIAA 2014-1171.
Ghanem, R. , 1998, “ Scales of Fluctuation and the Propagation of Uncertainty in Random Porous Media,” Water Resour. Res., 34(9), pp. 2123–2136.
Ghanem, R. , and Dham, S. , 1998, “ Stochastic Finite Element Analysis for Multiphase Flow in Heterogeneous Porous Media,” Transp. Porous Media, 32(3), pp. 239–262.
Najm, H. , 2009, “ Uncertainty Quantification and Polynomial Chaos Techniques in Computational Fluid Dynamics,” Annu. Rev. Fluid Mech., 41(1), pp. 35–52.
Le Maître, O. , Reagan, M. , Najm, H. , Ghanem, R. , and Knio, O. , 2002, “ A Stochastic Projection Method for Fluid Flow—II: Random Process,” J. Comput. Phys., 181(1), pp. 9–44.
Xiu, D. , and Karniadakis, G. , 2003, “ Modeling Uncertainty in Flow Simulations Via Generalized Polynomial Chaos,” J. Comput. Phys., 187(1), pp. 137–167.
Arnst, M. , Ghanem, R. , and Soize, C. , 2010, “ Identification of Bayesian Posteriors for Coefficients of Chaos Expansions,” J. Comput. Phys., 229(9), pp. 3134–3154.
Ghanem, R. , Doostan, A. , and Red-Horse, J. , 2008, “ A Probabilistic Construction of Model Validation,” Comput. Methods Appl. Mech. Eng., 197(29–32), pp. 2585–2595.
Ghanem, R. , and Red-Horse, J. , 1999, “ Propagation of Probabilistic Uncertainty in Complex Physical Systems Using a Stochastic Finite Element Approach,” Phys. D: Nonlinear Phenom., 133(1–4), pp. 137–144.
Marzouk, Y. M. , Najm, H. N. , and Rahn, L. , 2007, “ Stochastic Spectral Methods for Efficient Bayesian Solution of Inverse Problems,” J. Comput. Phys., 224(2), pp. 560–586.
Marzouk, Y. , and Najm, H. , 2009, “ Dimensionality Reduction and Polynomial Chaos Acceleration of Bayesian Inference in Inverse Problems,” J. Comput. Phys., 228(6), pp. 1862–1902.
Ghanem, R. , and Doostan, R. , 2006, “ Characterization of Stochastic System Parameters From Experimental Data: A Bayesian Inference Approach,” J. Comput. Phys., 217(1), pp. 63–81.
Sudret, B. , 2007, “ Global Sensitivity Analysis Using Polynomial Chaos Expansions,” Reliab. Eng. Syst. Saf., 93(7), pp. 964–979.
Crestaux, T. , Le Maître, O. , and Martinez, J. , 2009, “ Polynomial Chaos Expansion for Sensitivity Analysis,” Reliab. Eng. Syst. Saf., 94(7), pp. 1161–1172.
Le Maître, O. , and Knio, O. , 2015, “ Pc Analysis of Stochastic Differential Equations Driven by Wiener Noise,” Reliab. Eng. Syst. Saf., 135, pp. 107–124.
Lucor, D. , Meyers, J. , and Sagaut, P. , 2007, “ Sensitivity Analysis of Large-Eddy Simulations to Subgrid-Scale-Model Parametric Uncertainty Using Polynomial Chaos,” J. Fluid Mech., 585, pp. 255–279.
Blatman, G. , and Sudret, B. , 2008, “ Sparse Polynomial Chaos Expansions and Adaptive Stochastic Finite Elements Using a Regression Approach,” C. R. Méc., 336(6), pp. 518–523.
Peng, J. , Hampton, J. , and Doostan, A. , 2014, “ A Weighted ℓ1-Minimization Approach for Sparse Polynomial Chaos Expansions,” J. Comput. Phys., 267, pp. 92–111.
Yang, X. , and Karniadakis, G. , 2013, “ Reweighted ℓ1 Minimization Method for Stochastic Elliptic Differential Equations,” J. Comput. Phys., 248, pp. 87–108.
Hampton, J. , and Doostan, A. , 2015, “ Compressive Sampling of Polynomial Chaos Expansions: Convergence Analysis and Sampling Strategies,” J. Comput. Phys., 280, pp. 363–386.
Tipireddy, R. , and Ghanem, R. , 2014, “ Basis Adaptation in Homogeneous Chaos Spaces,” J. Comput. Phys., 259, pp. 304–317.
Tsilifis, P. , and Ghanem, R. , 2017, “ Reduced Wiener Chaos Representation of Random Fields Via Basis Adaptation and Projection,” J. Comput. Phys., 341, pp. 102–120.
Thimmisetty, C. , Tsilifis, P. , and Ghanem, R. , 2017, “ Homogeneous Chaos Basis Adaptation for Design Optimization Under Uncertainty: Application to the Oil Well Placement Problem,” Ai Edam, 31(3), pp. 265–276.
Yang, X. , Lei, H. , Baker, N. , and Lin, G. , 2016, “ Enhancing Sparsity of Hermite Polynomial Expansions by Iterative Rotations,” J. Comput. Phys., 307, pp. 94–109.
Tsilifis, P. , Huan, X. , Safta, C. , Sargsyan, K. , Lacaze, G. , Oefelein, J. , Najm, H. , and Ghanem, R. , 2018, “ Compressive Sensing Adaptation for Polynomial Chaos Expansions,” Preprint .
Tripathy, R. , Bilionis, I. , and Gonzalez, M. , 2016, “ Gaussian Processes With Built-in Dimensionality Reduction: Applications to High-Dimensional Uncertainty Propagation,” J. Comput. Phys., 321, pp. 191–223.
Tsilifis, P. , 2017, “ A Polynomial Chaos Basis Reduction Framework in Python,” PyPI: Python Package Index, accessed July 17, 2018,
Golub, G. , and Van Loan, C. , 2012, Matrix Computations, Vol. 3, JHU Press, Baltimore, MD.
Cho, C. , 1971, “ Convective Transport of Ammonium With Nitrification in Soil,” Can. J. Soil Sci., 51(3), pp. 339–350.
McNab, W. , and Narasimhan, T. , 1993, “ A Multiple Species Transport Model With Sequential Decay Chain Interactions in Heterogeneous Subsurface Environments,” Water Resour. Res., 29(8), pp. 2737–2746.
Van Genuchten, M. , 1985, “ Convective-Dispersive Transport of Solutes Involved in Sequential First-Order Decay Reactions,” Comput. Geosci., 11(2), pp. 129–147.
Pruess, K. , Oldenburg, C. , and Moridis, G. , 1999, “ TOUGH2 User's Guide, Version 2,” Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, Report No. LBNL-43134.
Oldenburg, C. , and Pruess, K. , 1995, “ EOS7R: Radionuclide Transport for TOUGH2,” Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, Report No. LBL-34868.
Soize, C. , and Desceliers, C. , 2010, “ Computational Aspects for Constructing Realizations of Polynomial Chaos in High Dimension,” SIAM J. Sci. Comput., 32(5), pp. 2820–2831.

## Figures

Fig. 1

Left: Plot of the empirical pdf of η=wTξ for w given in Eq. (39). Right: level-5 Clenshaw–Curtis quadrature points on [–1,1] and their mappings on η-space.

Fig. 2

Left: predictions of the link function using 1D- and 10D-expansions. Right: comparison of the true and the estimated density function bases on 1D-20th-order chaos expansion.

Fig. 3

Domain of the ammonium transport problem. Ammonium is injected at y = 0.

Fig. 4

Ten Monte Carlo realizations of the model outputs and their corresponding PC approximations using quadrature level  = 2, 3, 4, and 5, respectively

Fig. 5

Eigenvalues of the gradient matrix (left) and the values of the dominant eigenvector w (right) for chaos coefficients corresponding to different quadrature levels

Fig. 6

Estimated coefficients of the chaos expansions using different quadrature rule levels

Fig. 7

Evaluation of the model output and the 1d PC expansions on 1000 MC samples of η for  = 2, 3, 4, and 5

Fig. 8

Histogram comparison based on MC samples from the 1D- (left column) and 5D- (right column) chaos expansions with that based on 1000 MC samples drawn from the TOUGH2 simulator (all histograms are normalized such that the bins integrate to 1)

## Tables

Algorithm 1: Nonintrusive implementation
Table 1 Material and initial parameters used in my simulations
Table 2 Input parameter value range

## Errata

Some tools below are only available to our subscribers or users with an online account.

### Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections