An adaptive local reduced basis method for solving PDEs with uncertain inputs and evaluating risk

Accepted Manuscript An adaptive local reduced basis method for solving PDEs with uncertain inputs and evaluating risk Zilong Zou, Drew Kouri, Wilkins Aquino

PII: DOI: Reference:

S0045-7825(18)30532-2 https://doi.org/10.1016/j.cma.2018.10.028 CMA 12133

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date : 12 July 2018 Revised date : 17 October 2018 Accepted date : 18 October 2018 Please cite this article as: Z. Zou, et al., An adaptive local reduced basis method for solving PDEs with uncertain inputs and evaluating risk, Comput. Methods Appl. Mech. Engrg. (2018), https://doi.org/10.1016/j.cma.2018.10.028 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights (for review)

Research Highlights • We develop a methodology for adaptively constructing local reduced bases. • We provide rigorous derivations of a posteriori errors for risk measures as QoI. • We introduce an implicit Voronoi-partitioning that underpins the concept of local bases. • We demonstrate the performance of the method using a damped Helmholtz problem and an advection-diffusion problem. • We compare the efficiency and accuracy of our method to local and global sparse grid algorithms.

*Manuscript Click here to download Manuscript: localRB.pdf

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Click here to view linked References

An Adaptive Local Reduced Basis Method for Solving PDEs with Uncertain Inputs and Evaluating Risk Zilong Zoua , Drew Kourib , Wilkins Aquinoa,∗ a Department b Sandia

of Civil and Environmental Engineering, Duke University, Durham, NC 27708, USA National Laboratories, P.O. Box 5800, MS 0380, Albuquerque, NM 87185-9999, USA

Sandia National Laboratories is a multimission 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-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

Abstract Many physical systems are modeled using partial differential equations (PDEs) with uncertain or random inputs. For such systems, naively propagating a fixed number of samples of the input probability law (or an approximation thereof) through the PDE is often inadequate to accurately quantify the “risk” associated with critical system responses. In this paper, we develop a goal-oriented, adaptive sampling and local reduced basis approximation for PDEs with random inputs. Our method determines a set of samples and an associated (implicit) Voronoi partition of the parameter domain on which we build local reduced basis approximations of the PDE solution. The samples are selected in an adaptive manner using an a posteriori error indicator. A notable advantage of the proposed approach is that the computational cost of the approximation during the adaptive process remains constant. We provide theoretical error bounds for our approximation and numerically demonstrate the performance of our method when compared to widely used adaptive sparse grid techniques. In addition, we tailor our approach to accurately quantify the risk of quantities of interest that depend on the PDE solution. We demonstrate our method on an advection-diffusion example and a Helmholtz example. Keywords: stochastic partial differential equation, adaptive sampling, reduced basis, a posteriori error estimation, coherent risk measure, conditional value-at-risk

Preprint submitted to Elsevier

October 17, 2018

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1. Introduction Uncertainty is ubiquitous in virtually all engineering applications and, for a wide range of problems, it is inadequate to simulate the underlying physics without quantifying the uncertainty induced by random forces and coefficients, unknown boundary and initial conditions, and unverifiable modeling assumptions. Furthermore, the need to accurately quantify the risks associated with critical system responses is omnipresent. To numerically approximate these problems is computationally expensive since one often must solve a system of partial differential equations (PDEs) with uncertain inputs, which requires discretization in physical space, time and the uncertain input space. A measure of risk is a functional acting on a random quantity of interest (QoI) that assigns a numerical value for the overall “hazard” associated with large outcomes. A common approach to quantifying risk is to use the expected value of the QoI plus a measure of statistical deviation such as the standard deviation (see [35] and the references within). In many high consequence engineering applications, the mean plus standard deviation does not sufficiently characterize risk preference because it does not accurately characterize large magnitude, low probability events. In this paper, we restrict our attention to a specific class of risk measures called coherent risk measures [1]. This class is sufficiently broad and includes many risk measures that characterize tail events such as the conditional value-at-risk [36]. A difficulty with evaluating such risk measures is that they are typically not differentiable and thus approximation can be challenging. Indeed, accurately approximating the risk often requires an enormous number of samples of the QoI and hence an enormous number of PDE solves. Significant effort has been devoted to developing efficient methods for solving PDEs with random inputs. There are two prevalent approaches: projection methods and sampling methods. Projection methods seek to approximate the PDE solution by projecting onto a specified, cheaply computable basis (typically polynomial). Popular projection methods include polynomial chaos and stochastic Galerkin [4, 17, 24, 41]. On the other hand, sample-based approaches generate a finite set of realizations of the PDE solution ∗ Corresponding

author Email addresses: [email protected] (Zilong Zou), [email protected] (Drew Kouri), wa20@duke (Wilkins Aquino)

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

that are then used to compute statistics or are interpolated to build an approximation of the full solution. These approaches include Monte Carlo, stochastic collocation [3, 32, 40], stochastic reduced order models (SROMs) [26, 25]. Both polynomial chaos and stochastic collocation methods are based on polynomial approximation and obtain exponential convergence rates provided the PDE solution is sufficiently regular with respect to the random inputs. However, these methods typically suffer from the curse of dimensionality. Adaptive methods are promising candidates to mitigate the curse of dimensionality. Since the sensitivity of the PDE solution with respect to the random inputs is typically not uniform, one can achieve computational savings by biasing the computational effort towards parameter regions where the solution varies rapidly [6, 22, 23, 29, 30, 31, 39]. In fact, stochastic collocation on adaptive sparse grids has been successfully employed to solve optimization problems constrained by PDEs with random inputs [27, 28]. In addition to the above approaches for uncertainty quantification (UQ), reduced basis (RB) methods have proven to be powerful tools for efficiently solving time-dependent and parameter-dependent PDEs [7, 8, 9, 11, 12, 13, 14, 16, 34, 37]. Typical RB methods for parametrized PDEs construct a global basis using PDE solutions at select atoms. This basis is then used to approximate the PDE solution over the entire parameter domain. For these methods, a key concept is the offline-online computing strategy. During the offline stage, the basis and other large computational ingredients are computed and stored. In the online stage, only a reduced approximation problem needs to be solved at each new parameter. Since a large number of samples are typically needed to accurately evaluate QoI statistics, it is essential that one controls the cost of online evaluation. When the target PDE is affinely parametrized, the online cost can be made independent of the dimension of the spatial discretization. Detailed explanations of this approach can be found in [34, 37]. Fortunately, many non-affine problems can be approximated by affine counterparts using methods like empirical interpolation [5]. The main contribution of this work is the development of an adaptive local RB approximation for solving PDEs with random inputs and approximating risk measures evaluated at a QoI that depends on the PDE solution. Our method adaptively constructs a set of parameter samples which gives rise to an implicit Voronoi partition of the parameter domain. Using this decomposition, we construct local RB approximations of 3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

the PDE solution within each Voronoi cell. In contrast to existing RB methods, we form local bases for each Voronoi cell using solutions only at a fixed number of proximal atoms as well as the gradient information at the generating atom. Due to the local nature of our proposed method, the online computational cost does not increase as more atoms are added in the adaptive process. In addition, thanks to the gradient information, the solution of our reduced problem is guaranteed to converge almost everywhere over the parameter domain to the solution of the full problem provided certain smoothness assumptions hold. We build our sample set on the fly by adapting the greedy approach utilized in [34, 37] and guiding the adaptivity using residual-based a posteriori error indicators. Our error indicators require the evaluation of the stability (inf-sup) constant of the parametrized PDE operator for each realization of the random input. To circumvent the need to solve an eigenvalue problem for the exact stability constant at each sample, we develop a surrogate model of the constant using dimension-adaptive sparse grids [22, 23]. Evaluating the error associated with the adaptive sparse grid surrogate is nontrivial and problem dependent. Hence, our computable error indicator is, in general, a heuristic. Our adaptive sampling approach permits a natural extension to approximate coherent risk measures of QoIs that depend on the PDE solution. To this end, we develop a rigorous bound for the error in the risk approximation and subsequently derive an error indicator that is compatible with our adaptive sampling framework. With the risk error indicator, we can easily perform adaptive sampling to specifically target risk approximation. To conclude, we demonstrate the efficiency of our method through numerical examples and comparison with existing methods. This paper is organized as follows. First, we describe the problem formulation including a discussion of PDEs with uncertain inputs and local Taylor-like approximations. Following this discussion, we introduce our adaptive local RB approach including a priori and a posteriori error analysis. We then extend this analysis to handle the approximation of coherent risk measures. We conclude with numerical examples that confirm our results.

4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2. Background 2.1. Problem Formulation Let (Ω, F, P) be a probability space, where Ω is the sample space, F ⊆ 2Ω is a σ-

algebra of events, and P : F → [0, 1] is a probability measure. We denote the expectation R of a random variable X : Ω → R by E[X] := Ω X(ω) dP(ω) and by almost surely (a.s.)

we mean “up to a set of P-measure equal to zero.” We denote the space of random variables on (Ω, F, P) with p finite moments, p ∈ [1, ∞), by Lp (Ω, F, P) and the space

of essentially bounded random variables on (Ω, F, P) by L∞ (Ω, F, P). Now, let W be

an arbitrary Banach space. We denote the norm on W by k · kW and the topological

dual space of W by W ∗ (i.e., W ∗ is the Banach space consisting of all bounded linear functionals acting on elements of W ). We denote the action of a functional w∗ ∈ W ∗

on w ∈ W by w∗ (w) = hw∗ , wiW ∗ ,W . If W is a Hilbert space, then we denote the inner product on W by (·, ·)W .

Let D ⊂ Rd , d = 1, 2, 3, denote the physical domain and let U = U (D) and V = V (D)

denote Hilbert spaces of (deterministic) functions mapping D into R. We consider the following parametric equation: For fixed f ∈ V ∗ , find u : Ω → U such that b L(ω)u(ω) = f,

a.s.

(1)

b Here L(ω) is a bounded linear operator from U into V ∗ for P-almost all (a.a.) ω ∈

Ω. To facilitate numerical approximation, we assume that the finite-dimensional noise assumption holds [40]. That is, we assume there exists a random vector, ξ : Ω → Ξ with

Ξ := ξ(Ω) ⊆ RM , and a bounded linear operator L(ξ) mapping U into V ∗ a.s. satisfying b L(ω) = L(ξ(ω))

a.s.

We denote the joint distribution of ξ by Fξ and assume that the probability law of ξ, P ◦ ξ −1 , is absolutely continuous with respect to the Lebesgue measure with density ρ : Ξ → [0, ∞). Here and throughout, we abuse notation and denote the random vector and its realizations by ξ. We note here that uncertainty can also be incorporated in the right-hand side f , but for ease of presentation, we restrict our attention to uncertainty in the operator. Through a change of variables, (1) becomes: Find u : Ξ → U such that L(ξ)u(ξ) = f, 5

∀ξ ∈ Ξ.

(2)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In the subsequent analysis, it will be convenient to define the parametrized bilinear form a(·, ·; ξ) : U × V → R associated with L(ξ) for ξ ∈ Ξ by a(u, v; ξ) := hL(ξ)u, viV ∗ ,V

∀ u ∈ U, v ∈ V.

We can then rewrite (2) as the variational problem: Find u : Ξ → U such that a(u(ξ), v; ξ) = hf, viV ∗ ,V ∀v ∈ V, ξ ∈ Ξ.

(3)

We assume that the spaces U and V are chosen so that problem (3) is well-posed. In particular, we assume that there exists positive constants 0 < m1 ≤ m2 , independent of ξ ∈ Ξ, such that for all ξ ∈ Ξ, hL(ξ)u, viV ∗ ,V ≤ m2 kukU kvkV inf

sup

06=u∈U 06=v∈V

∀u ∈ U, v ∈ V

(4a)

|hL(ξ)u, viV ∗ ,V | =: γ(ξ) > m1 kukU kvkV

(4b)

hL(ξ)u, viV ∗ ,V = hL(ξ)∗ v, uiU ∗ ,U = 0 ∀ u ∈ U

=⇒

v = 0.

(4c)

Note that when U and V are finite-dimensional, (3) represents, e.g., the discretization of a PDE. Our goal is to build an approximation u : Ξ → U of u that can be evaluated at ξ ∈ Ξ efficiently without the need to solve the full problem (3). Once we construct u, we can approximate the statistics of u using the statistics of u. For example, given a function g : U → R, the expectation of g(u(ξ)) can be approximated using Monte Carlo as E[g(u)] ≈ E[g(u)] =

Z

Ξ

g(u(ξ))ρ(ξ)dξ ≈

N 1 X g(u(ξi )) N i=1

(5)

where ξi are independent and identically distributed random vectors. 2.2. Deterministic Samples and Local Taylor Approximation Our proposed surrogate, u, relies on a set of deterministic samples distributed in Ξ. These samples generate an implicit Voronoi partition of the parameter domain on which we construct local RB approximations. This method is inspired by local Taylor approximation on a finite set of (deterministic) samples of ξ. In particular, we draw motivation from SROMs [25, 26]. We will review some critical facts about Taylor approximation that motivate our approximation method. 6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Any finite set of (deterministic) samples of ξ, denoted {ξ k }m k=1 (e.g., computed as a SROM of ξ), generates an implicit Voronoi partition of the parameter domain Ξ = ∪m k=1 Ξk where {Ξk } are the Voronoi cells generated by the atoms {ξ k }, i.e., Ξk := {ξ ∈ Ξ : kξ − ξ k k2 ≤ kξ − ξ j k2 , j 6= k}. Here, k · k2 denotes the Euclidean norm on Rm . Using this partition, we can assign each

sample ξ k the probability pk := P({ω ∈ Ω : ξ(ω) ∈ Ξk }) for k = 1, . . . , m. Clearly, we Pm S have that k=1 pk = 1 since Ξ = Ξk and P ◦ ξ −1 (Ξi ∩ Ξj ) = 0 for i 6= j, which follows from the absolute continuity of P ◦ ξ −1 with respect to the Lebesgue measure. Given the

samples {ξ k } and assuming sufficient regularity of the random field PDE solution u, we can approximate u using a Taylor expansion over Ξk anchored at the sample ξ k in the usual way. For example, the local linear approximation of u(ξ) is u0 (ξ) :=

m X

k=1

Here,





1Ξk (ξ) u(ξ k ) + ∇ξ u(ξ k ) · (ξ − ξ k ) .

(6)

1Ξk (ξ) denotes the characteristic function of the set Ξk . That is, 1Ξk (ξ) is one

if ξ ∈ Ξk and is zero otherwise. The computational cost of building u0 is completely dictated by the number of atoms m. The following theorem provides an error bound for the above Taylor approximation (6). It is adapted from [25, Theorem 2]. Theorem 1. Let Ξ ⊇ Ξ be an open set in RM and suppose u : Ξ → U is continuously Fr´echet differentiable. Moreover, we assume that there exists L ≥ 0 such that max k∂i u(ξ) − ∂i u(ξ 0 )kU ≤ Lkξ − ξ 0 k

i=1,...,M

∀ ξ, ξ 0 ∈ Ξ

(7)

where ∂i denotes the partial derivative with respect to the ith component of ξ. Then, m

ku(ξ) − u0 (ξ)kU ≤

LX 1Ξk (ξ)kξ − ξ k k22 2 k=1

∀ ξ ∈ Ξ.

(8)

In particular, if Ξ is bounded then sup ku(ξ) − u0 (ξ)kU ≤ ξ∈Ξ

L max diam(Ξk )2 2 k=1,...,m

where the diameter of Ξk is defined by diam(Ξk ) := sup kξ − ξ 0 k2 ξ,ξ 0 ∈Ξk

7

(9)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Proof. Notice that subtracting u0 (ξ) from u(ξ) and applying the integral mean value theorem gives u(ξ) − u0 (ξ) =

m X

k=1

1Ξk (ξ)

Z

0

1

(∇ξ u(ξ k + t(ξ − ξ k )) − ∇ξ u(ξ k )) · (ξ − ξ k )dt.

Computing the norm of both sides and applying Jensen’s and the Cauchy-Schwarz inequalities as well as the Lipschitz continuity bound (7) yields (8). When Ξ is bounded, (9) follows trivially from (8) and H¨ older’s inequality. As can be seen, the Taylor approximation is a special discontinuous collocation method. The collocation points, {ξ k }m k=1 , are selected according to the probability law of the input random parameters ξ. Within some vicinity of each ξ k , the truncated Taylor expansion is accurate. In the context of SROMs, the authors in [21] demonstrate the efficacy of (6) to solve stochastic problems and its specific advantages over the traditional stochastic collocation and stochastic Galerkin methods. Local approximation based on Taylor expansion is only accurate when the size of the cells Ξk is sufficiently small, as can be seen in (8). Moreover, for problems in which the gradient of u varies rapidly with the random parameters, the truncated Taylor expansion is likely to yield a large error that can only be mitigated by refining the local partition with additional atoms. The number of samples required to generate a sufficiently refined partition can be quite large, which leads to excessive computational cost. On the other hand, SROM samples {ξ k } are constructed to approximate the input random parameters and do not necessarily yield a good approximation of the PDE solution. In the subsequent sections, we will develop a method that improves upon the local Taylor approach. In particular, we will replace the truncated Taylor expansion with a RB approximation and employ adaptivity to select the samples, and hence to partition the parameter domain. The sample locations are chosen using an a posteriori error indicator. For certain problems, these changes yield significant improvement in terms of the approximation quality of the PDE solution.

8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3. The Adaptive Local Reduced Basis Method 3.1. Local Reduced Basis Approximation Our approximation approach is motivated by the observation that the local first-order Taylor approximation (6) can be rewritten as u0 (ξ) =

m X

1Ξk (ξ) Φk,0 wk (ξ)

(10)

k=1

where   Φk,0 = u(ξ k ), ∂1 u(ξ k ), ∂2 u(ξ k ), . . . , ∂M u(ξ k ) >  wk (ξ) = 1, (ξ1 − ξ k,1 ), (ξ2 − ξ k,2 ), . . . , (ξM − ξ k,M ) .

(11) (12)

Here, ξi and ξ k,i denote the ith components of ξ and ξ k , respectively. Thus, for fixed ξ ∈ Ξk , the surrogate in (6) lies in the linear span of Φk,0 with coefficients wk (ξ). As discussed above, this approximation is only accurate when Ξk is “sufficiently small” and is not optimal in terms of fully exploiting the range of Φk,0 . We address this deficiency using Galerkin projection as a means of computing the coefficient vector wk (ξ). To provide more flexibility and to make the proposed method more robust, we provide the option to extend the range of Φk,0 with additional sampled solutions at proximal points. In particular, we append to Φk,0 the solutions u(ξ j ) at the closest N (N ≤ m − 1) atoms. That is,

  Φk = Φk,0 , u(ξ k1 ), u(ξ k2 ), . . . , u(ξ kN ) .

(13)

An example of the proposed basis is shown in Figure 1, where a partition of the parameter domain Ξ (with two random parameters) is generated by a group of atoms (red circles and dots). The surrogate solution at the blue dot is computed using a basis consisting of the solution at the red dots as well as the solution and gradient at the black dot. Using Φk , we define our surrogate, u, as u(ξ) :=

m X

1Ξk (ξ)uk (ξ)

(14)

k=1

where uk (ξ) solves the following variational problem: Find uk : Ξk → Uk such that a(uk , vk ; ξ) = hf, vk iV ∗ ,V ∀vk ∈ Vk (ξ), ξ ∈ Ξk . 9

(15)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

Figure 1: Example of the local reduced basis method with two random parameters and N = 3. The surrogate solution at the blue dot is computed using a basis consisting of the solution at the red dots as well as the solution and gradient at the black dot. We plot 2,000 Monte Carlo samples of ξ in the background.

Here, Uk := span(Φk ) is the reduced approximation space and Vk (ξ) ⊂ V is the reduced test space that is possibly dependent on ξ. In general, (15) is not automatically well-posed for an arbitrary choice of Vk (ξ). To guarantee the well-posedness of (15), we extend the concept of optimal testing. Similar ideas have also been used in [19] in the context of discontinuous Petrov-Galerkin method and in [34, 37] for RB applications. To this end, we first introduce the operator T (ξ) := RV−1 ◦ L(ξ) where RV−1 : V ∗ → V is the inverse Riesz isomorphism. It is straightforward to verify that (T (ξ)u, v)V = a(u, v; ξ) ∀u ∈ U, v ∈ V

(16)

we then define the test space as Vk (ξ) := T (ξ)Uk . Due to the properties of a(·, ·; ξ) (4b), T (ξ) is injective. Obviously, Vk (ξ) ⊂ V since Uk ⊂ U . With this test space, the following theorem demonstrates the well-posedness of (15). Theorem 2. If Vk (ξ) = T (ξ)Uk and (4b) holds. Then, we have that |a(u, v; ξ)| ≥ γ(ξ) > m1 06=u∈Uk 06=v∈Vk (ξ) kukU kvkV inf

sup

10

∀ ξ ∈ Ξ.

(17)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Proof. We have that |a(u, v; ξ)| |a(u, T (ξ)u; ξ)| |(T (ξ)u, T (ξ)u)V | ≥ inf = inf 06=u∈Uk kukU kT (ξ)ukV 06=u∈Uk kukU kT (ξ)ukV 06=u∈Uk 06=v∈Vk (ξ) kukU kvkV inf

sup

=

inf

|(T (ξ)u, v)V | |a(u, v; ξ)| = inf sup 06=u∈Uk 06=v∈V kukU kvkV kukU kvkV

sup

06=u∈Uk 06=v∈V

≥ inf

|a(u, v; ξ)| = γ(ξ) > m1 kukU kvkV

sup

06=u∈U 06=v∈V

where the first and third equalities follow from the definition of T (ξ) and the second equality follows from kT (ξ)uk2V |(T (ξ)u, T (ξ)u)V | |(T (ξ)u, v)V | = = kT (ξ)ukV = kL(ξ)ukV ∗ = sup . kT (ξ)ukV kT (ξ)ukV kvkV 06=v∈V

We have the following a priori error bound for our surrogate u(ξ). Theorem 3. Let the assumptions of Theorem 1 hold. Then,   m m2 L X 1Ξk (ξ)kξ − ξ k k22 ku(ξ) − u(ξ)kU ≤ 1 + m1 2 k=1

∀ ξ ∈ Ξ.

(18)

Proof. We have that ku(ξ) − u(ξ)kU ≤ ≤ ≤







1+ 1+

m2 m1 m2 m1

m2 1+ m1

X m

k=1 X m



k=1

inf kw − u(ξ)kU

w∈Uk

ku0 (ξ) − u(ξ)kU m

LX 1Ξk (ξ)kξ − ξ k k22 2 k=1

for all ξ ∈ Ξ where the first inequality follows from Theorem 2.2 in [2], the second inequality is due to the fact that u0 (ξ) ∈ Uk and the final is due to Theorem 1.

Note that Vk (ξ) is dependent on ξ, which is essential to establish the stability of the reduced problem (15) for arbitrary ξ ∈ Ξ. However, for uniformly coercive elliptic problems, we can set Vk = Uk . Due to the coercivity of the bilinear form a(·, ·; ξ), the restriction argument implies that the problem in (15) is automatically well posed. A similar error bound can be derived as in Theorem 3. 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Finally, note that since N is a fixed number that does not grow with the introduction of more atoms, the extended basis Φk is still local and the Galerkin projection (15) has at most M + N + 1 dimensions which are independent of the number of atoms m. When N = 0, we recover the Taylor basis Φk,0 that is used by the SROM method. On the other hand, when N = m − 1, we recover a global reduced basis method similar to that used in [37], but with additional local gradient enrichment. We define Nk = dim(Φk ) as the reduced problem dimension. Obviously, Nk ≤ M + N + 1. We achieve a reduction in the computational cost since the dimension of the full problem (3) is much larger than Nk . In many cases, the cost of solving the latter is negligible in comparison to the cost of solving the full problem. We want to point out that due to the addition of gradients, the size of the basis scales with M . In cases where M is very large, the computational cost of the reduced problem may render our method impractical. In this case, basis compression methods, such as Proper Orthogonal Decomposition (POD), can be used to construct lower dimensional bases and maintain computational efficiency. However, the main drawback of these approaches is the associated difficulty in guaranteeing the error bound in Theorem 3.

3.2. A Posteriori Error Estimation and Adaptive Sampling We now describe our adaptive approach for generating the atoms set Θ := {ξ k }m k=1 . Our sequential selection of Θ is guided by the current surrogate solution u. Given a current set of atoms, we select the next atom ξ k+1 by identifying the region of Ξ with the largest error, which we estimate using a reliable a posteriori error estimate u (ξ). In this work, we employ residual-based error estimation as the foundation for adaptivity. Residual-based error estimators have been used with great success in adaptive finite element methods [10] and in the adaptive construction of reduced bases [34]. We first define the surrogate error as e(ξ) = u(ξ) − u(ξ) and note that e satisfies a(e(ξ), v; ξ) = r(u(ξ), v; ξ) ∀v ∈ V

(19)

where r(u(ξ), v; ξ) = hf, viV ∗ ,V − a(u(ξ), v; ξ) is the residual associated with u(ξ). Note that r(u(ξ), ·; ξ) ∈ V ∗ is a bounded linear functional on V . Hence, from the inf-sup 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

condition (4b), we have γ(ξ)ke(ξ)kU ≤ sup v∈V

a(e(ξ), v; ξ) = kr(u(ξ), ·; ξ)kV ∗ . kvkV

(20)

That is, ke(ξ)kU ≤

kr(u(ξ), ·; ξ)kV ∗ =: 0u (ξ). γ(ξ)

(21)

The error estimator (21) requires explicit knowledge of the inf-sup constant for each ξ ∈ Ξ, which is not readily available and can be computationally expensive to obtain. In particular, the exact inf-sup constant, γ(ξ), for a single ξ requires the solution of a generalized eigenvalue problem with size equal to that of the full problem. To circumvent this issue, we build an efficient surrogate model for γ(ξ), denoted by γ b(ξ). Assuming that γ b(ξ) satisfies the approximation property ∃0<δ<1

such that

sup ξ∈Ξ

|γ(ξ) − γ b(ξ)| ≤ δ, γ(ξ)

(22)

we have that (1 − δ)γ(ξ) ≤ γ b(ξ) ≤ (1 + δ)γ(ξ) and we can bound the error as ke(ξ)kU ≤ u (ξ) :=

kr(u(ξ), ·; ξ)kV ∗ θb γ (ξ)

(23)

where θ = (1 + δ)−1 . Of course, accurately determining δ is generally a nontrivial task, because evaluating the error associated with γ b(ξ) requires solving additional eigenvalue problems.

Upon building the error indicator u (ξ), we develop an adaptive algorithm that re-

duces the expected error E [ke(ξ)kU ] in a greedy manner. Thanks to the Voronoi partition of the parameter domain, we can decompose the error into a local error indicator in each cell as E [ke(ξ)kU ] =

m X

k=1

E [ke(ξ)kU 1Ξk (ξ)] ≤

m X

k=1

E [u (ξ)1Ξk (ξ)] =

m X

ηk

(24)

k=1

where ηk := E [u (ξ)1Ξk (ξ)] is the local error indicator. In a greedy method, the next atom is selected from the cell with the largest ηk . To compute ηk , an intuitive approach would be to employ Monte Carlo integration via a set of training samples that captures the probability law of ξ. However, such an approach has two drawbacks. First, as the dimension of ξ becomes large, a significant number of training samples is required. Second, a predefined, fixed training set may 13

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

not provide adequate representation of the solution behavior in some regions of the parameter space. To overcome these limitations, we propose an approach to generate training samples on the fly. To this end, notice that ηk = pk E [u (ξ)|ξ ∈ Ξk ]

(25)

where pk = P({ω ∈ Ω : ξ(ω) ∈ Ξk }) is the probability of Ξk and E [u (ξ)|ξ ∈ Ξk ] is the conditional expectation of u (ξ) given that ξ ∈ Ξk . Hence we can estimate ηk by estimating pk and E [u (ξ)|ξ ∈ Ξk ] individually.

To estimate pk , we maintain a population (Nbkg ) of background Monte Carlo samples Θbkg . Provided that Nbkg is sufficiently large, we have pk ≈

|Θbkg ∩ Ξk | Nbkg

(26)

where |A| denotes the number of samples in A. This computation is realized by the implicit Voronoi tessellation. This step is computationally inexpensive since no PDE or reduced problem needs to be solved. The conditional expectation E [u (ξ)|ξ ∈ Ξk ] can be estimated in a similar fashion. Let Θtr ⊂ Ξ be a dynamic training set that adapts to our need to explore the parameter

domain with the largest local error. We can estimate E [u (ξ)|ξ ∈ Ξk ] as P ξ∈Θtr ∩Ξk u (ξ) E [u (ξ)|ξ ∈ Ξk ] ≈ |Θtr ∩ Ξk | and we hence have |Θbkg ∩ Ξk | ηk ≈ × Nbkg

P

ξ∈Θtr ∩Ξk u (ξ)

|Θtr ∩ Ξk |

(27)

(28)

which is the computable local error indicator for each cell Ξk . Upon computing ηk , the next atom ξ m+1 is selected from Θtr ∩ Ξk as the one that maximizes u (ξ) over this set, i.e., ξ m+1 = arg maxξ∈Θtr ∩Ξk u (ξ) and Θ is updated to Θ ∪ {ξ m+1 }. A new partition is then implicitly formed using the updated Θ. The next step is to update the training set Θtr so that there are enough training min samples for further exploration. In particular, we require at least Ntr training samples

inside each Voronoi-cell of the new partition. To this end, we verify the cell identity of the training samples by implicit Voronoi tessellation again. We identify cells with insufficient number of training samples and draw additional training samples in those cells using 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

min rejection sampling to ensure that mink=1,...,m+1 |Ξtr ∩ Ξk | ≥ Ntr . To initialize the

init algorithm, we start with a small (Ntr ) of training samples in Θtr and only one atom

Θ = {E[ξ]}. Algorithm 1: Adaptive sequential construction of Θ and u(ξ) Initialization: • Specify the maximum number of atoms m and the error tolerance Errtol . N

bkg • Form a background set Θbkg := {ξi }i=1 consisting of a sufficiently large

number of Monte Carlo samples of ξ. N init

tr • Form an initial training set Θtr := {ξi }i=1

consisting of a few random

samples of ξ. • Select an initial atom Θ = {E[ξ]} and build the initial surrogate model u(ξ) based on the solution u(E[ξ]) and gradient ∇ξ u(E[ξ]). • Set Natom = 1 and Errmax = +∞. while Natom < m and Errmax > Errtol do Evaluate the error indicator u (ξ) of u(ξ) at each ξ ∈ Θtr ; Compute ηk , k = 1, 2, . . . , Natom via implicit Voronoi tessellation of Ξ; Set k = arg max ηk , Errmax = k=1,...,Natom

max

k=1,2,...,Natom

ηk ;

Set ξmax = arg max u (ξ), Θ = Θ ∪ {ξmax } and Natom = Natom + 1; ξ∈Θtr ∩Ξk

Compute u(ξmax ) and gradient ∇ξ u(ξmax ); Incorporate the new information at {ξmax } into the surrogate u(ξ); Draw additional training samples to ensure

min

k=1,...,Natom

min |Ξtr ∩ Ξk | ≥ Ntr .

end Store the atoms Θ and surrogate model u(ξ) for future computations.

3.3. Adaptive Local Estimation for Coherent Risk Measures In this section, we extend the Local RB method and the adaptive sampling procedure to approximate coherent measures of risk evaluated at a quantity of interest (QoI) that depends on the PDE solutions. 15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

As mentioned, many applications require the accurate approximation of the “risk” associated with critical system responses. We will introduce an axiomatic definition of a measure of risk and describe how our adapted reduce basis approximation can be extended to approximate risk. In general, we are interested in quantifying the risk of some QoI g : U → R evaluated at the PDE solution u. Suppose the random variable

G = g(u(ξ)) ∈ X where X := Lp (Ω, F, P) with 1 ≤ p < ∞ and let R : X → (−∞, +∞] denote our target measure of risk. A measure of risk, R, is coherent [1] if it satisfies:

(R1) Subadditivity: R(X + X 0 ) ≤ R(X) + R(X 0 ) for all X, X 0 ∈ X ;

(R2) Translation Equivariance: R(X + t) = R(X) + t for all X ∈ X and for all t ∈ R; (R3) Monotonicity: X, X 0 ∈ X with X ≤ X 0 a.s. implies R(X) ≤ R(X 0 ); (R4) Positive Homogeneity: R(tX) = tR(X) for all X ∈ X and t ≥ 0. Properties (R1) and (R4) guarantee that R is convex. Moreover, Theorem 6.7 in [38] ensures that R(X) = sup E[ϑX]

(29)

ϑ∈A

where the risk envelope, A, is a convex bounded and weak∗ closed subset of X ∗ consisting of probability density functions. That is, A ⊆ { θ ∈ X ∗ : E[θ] = 1, θ ≥ 0 a.s. } . In particular, property (R3) implies θ ≥ 0 a.s. for all θ ∈ A and property (R2) ensures if θ ∈ A then E[θ] = 1.

One popular coherent measure of risk is the conditional value-at-risk (CVaR)   Z 1 1 1 qα (X) dα = inf t + E[(X − t)+ ] R(X) = CVaRβ (X) := t∈R 1−β β 1−β

(30)

for 0 ≤ β < 1, where qα (X) := inf{η ∈ R : P(X ≤ η) ≥ α} is the upper α-quantile of X and (·)+ := max{·, 0}. When the distribution function of the random variable X is continuous, CVaRβ (X) is the expected value of X conditioned on the event that X is larger than its β-quantile [36]. Additionally, the risk envelope A for CVaRβ is given by   1 ∗ A = θ ∈ X : E[θ] = 1, 0 ≤ θ ≤ a.s. . 1−β Another popular coherent measure of risk is the mean-plus-upper-semideviation of order p ∈ [1, ∞),

1

R(X) = E[X] + cE[(X − E[X])p+ ] p , 16

0 ≤ c ≤ 1.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The risk envelope A for this risk measure is given by A = { θ ∈ X ∗ : θ = 1 + θ0 − E[θ0 ], kθ0 kX ∗ ≤ c, θ0 ≥ 0 a.s.} . See [38, § 6.3] for more examples of coherent risk meaures. For general coherent risk measures, using the subadditivity (R1), the monotonicity (R3) and the biconjugate representation (29), we have that for any X, X 0 ∈ X , R(X) = R(X 0 + (X − X 0 )) ≤ R(X 0 ) + R(X − X 0 ) ≤ R(X 0 ) + R(|X − X 0 |) = R(X 0 ) + sup E[ϑ|X − X 0 |]. ϑ∈A

The same bound holds if we exchange X with X 0 and thus we have the error bound |R(X) − R(X 0 )| ≤ sup E[ϑ|X − X 0 |].

(31)

ϑ∈A

Returning to our original goal of approximating R(G), we have that (31) holds with

X = G and X 0 = G = g(u(ξ)). To bound the right-hand side of (31), we employ our a posteriori error indicator. Given the current atom set Θ, the surrogate, u, is locally defined over the implicit Voronoi cells Ξk , k = 1, . . . , m. Using the explicit representation (14), we define Gk = g(uk (ξ)) and note that |R(G) − R(G)| ≤ sup

ϑ∈A

m X

k=1

E[ϑ1Ξk (ξ)|G − Gk |] ≤

m X

k=1

sup E[ϑ1Ξk (ξ)|G − Gk |].

ϑ∈A

Hence, the error in the risk of our surrogate model is bounded by the expected worst case error over the Voronoi tessellation. This bound provides a natural means to balance the error. Intuitively, we should refine cells whose local error is large. To arrive at a practical local error estimate, we assume the functional g is H¨ older continuous. That is, there exists K > 0 and α > 0 such that |g(w) − g(w0 )| ≤ Kkw − w0 kα U

∀ w, w0 ∈ U.

Under this assumption, we arrive at E[ϑ1Ξk (ξ)|G − Gk |] ≤ KE[ϑ1Ξk (ξ)ku(ξ) − uk (ξ)kα U] 17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

for all θ ∈ A and all k = 1, . . . , m. Therefore, we can bound the right-hand side of (31) using our a posteriori error indicator (21) as |R(G) − R(G)| ≤ K

m X

k=1

sup E[ϑ1Ξk (ξ)ke(ξ)kα U] ≤ K

ϑ∈A

m X

k=1

sup E[ϑ1Ξk (ξ)0u (ξ)α ].

ϑ∈A

Our local cell error estimate is thus δk0 := supϑ∈A E[ϑ1Ξk (ξ)0u (ξ)α ] = R(1Ξk (ξ)0u (ξ)α )

whereas our computable (heuristic) local error indicator is δk := supϑ∈A E[ϑ1Ξk (ξ)u (ξ)α ] = R(1Ξk (ξ)u (ξ)α ). We note that when approximating a QoI and an associated risk measure, we can employ a dual-weighted residual method to derive an error indicator for G that is independent of the inf-sup constant [15, 20]. However, the dual-weighted error indicator would require the additional expense of a reduced adjoint solve at each sample. Recall that in Section 3.2 we have a cell error indicator ηk which is used in the greedy adaptive sampling approach in Algorithm 1. It is straight forward to modify Algorithm 1 to account for the error associated with the risk computation by replacing ηk with δk . This substitution provides an adaptive reduced basis approach for efficiently evaluating the risk associated with QoIs that depends on the PDE solution. 3.4. Computational Cost In this subsection, we derive a coarse estimate of the computational complexity associated with each stage of our method. To this end, we assume that U and V are finite dimensional with dimension Nspace and that M + N  Nspace . In our estimate, we omit constants that are neither an algorithmic parameter nor a dimension of the discretization. When adding estimates, we omit terms that are clearly comparable to or dominated by another term. In practice, the indicator u (ξ) is a conservative upper bound of the true error ke(ξ)kU and only serves as a guide for the adaptive algorithm to select the next atom. Therefore, a somewhat loose approximation u (ξ) may be sufficient to construct an accurate surrogate

u. To build a satisfactory surrogate model γ b(ξ), we employ dimension-adaptive sparse

grid interpolation [22, 23, 30]. The dimension-adaptive sparse grid method attempts to find the important stochastic dimensions along which the stability constant varies

most rapidly and places more interpolation points along such important dimensions [23]. 18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Of course, error bounds of the form (22) are not readily available for adaptive sparse grid interpolation and thus our error indicator (23) is heuristic. We limit the number of collocation points to be Nγ . Given that we only need moderate accuracy of the surrogate model, Nγ can be chosen to be a relatively small number. For instance, in all our examples, we chose Nγ = max(100, 10M ). We assume a conservative cost of the solution 3 of the smallest generalized eigenvalue as O(Nspace ) [33], although more refined estimates 3 are available for sparse systems. Hence γ b(ξ) is built with O(Nγ Nspace ) operations. We

would like to point out that a similar or higher computational cost is required in existing RB methods to construct a lower bound for the stability constant [34]. Once γ b(ξ) is constructed, evaluating it at a new parameter takes negligible cost.

Our adaptive procedure is similar to the greedy sampling of conventional RB methods

[34]. For m atoms, the number of high fidelity solutions is m(M + 1). Hence, the 3 operation count is O(m(M + 1)Nspace ). The key to RB methods is the assumption of an

(approximate) affine parametrization of the weak formulation. Given this assumption, the online cost to solve a reduced problem and evaluate the a posteriori error estimator (23) can be made independent of Nspace by precomputing and storing offline the affine components of the reduced system and those of the dual norm of the residual. Using the test basis as in Theorem 2 induces little extra cost since the affine components are readily available as byproducts of the residual norm. The details of this offline preprocessing are similar to those in the existing RB methods [37]. The most expensive part in the preprocessing steps is the repeated application of an inverse Riesz map RV−1 : V ∗ → V . However, note that RV is ξ-independent, which means that we only need

2 to factorize RV once. After the factorization, the cost of applying RV−1 is O(Nspace ). 3 2 Overall, an additional cost of O(Nspace ) + O(mM (M + N + 1)Nspace ) is required for the

offline preprocessing. The last term can be absorbed into the former cost estimate given that Nspace is the dominant dimension, yielding a combined offline construction cost of 3 O((1 + mM + m)Nspace ).

Under the assumption of an affine parametrization, the reduced problem and error indicator u (ξ) can be computed with O(Nk3 + M 2 Nk2 ) operations at each sample, where Nk ≤ M + N + 1 is the size of the reduced problem and is by assumption much smaller than Nspace . If a training set Θtr , as in the existing RB methods, is used in the adaptive 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

sampling process, an overhead of O(m |Θtr | (Nk3 + M 2 Nk2 )) is added to the construction cost of the surrogate model. The latter cost is not negligible if the cardinality of Θtr , which is assumed to be an exhaustive set to capture the input uncertainty, is large. In contrast, our approach uses a dynamic training set of much smaller cardinality. The init min total number of our training samples is bounded by Ntr +mNtr , which is substantially

smaller than the cardinality of an exhaustive training set |Θtr |. Hence the construction

3 init min cost for our method is O (1 + mM + m)Nspace +O m(Ntr + mNtr )(Nk3 + M 2 Nk2 ) . In addition, a potential saving for our method can be realized by exploiting the

gradient computations. At each atom, we perform one high fidelity solve for the state u and M high fidelity solves for the partial derivatives. However, these 1 + M solutions share the same operator. If a direct solver were to be used, only one factorization per atom would be required. The latter would reduce our construction cost to

3 init min O (1 + m)Nspace + O m(Ntr + m2 Ntr )(Nk3 + M 2 Nk2 ) . However, we do not exploit this computational advantage in our numerical examples when comparing to other collocation based methods. In summary, the overall computational cost to build a surrogate model using our adap

3 init min tive approach is O (Nγ + 1 + mM + m)Nspace +O m(Ntr + m2 Ntr )(Nk3 + M 2 Nk2 ) .

Given that our choices of the algorithmic parameters typically results in a very small

training set, the last term is negligible in comparison to the first one in large scale prob
3 . lems. Therefore, the overall cost can be approximated as O (Nγ + 1 + mM + m)Nspace

Once a surrogate model is constructed, the online computational cost to obtain a

surrogate solution at a new sample is simply O(Nk3 + M 2 Nk2 ). Given a fixed N , this number is independent of Nspace as well as of the number of atoms used to construct the surrogate model m.

4. Numerical Examples In this section, we present two numerical examples to showcase the capability of the proposed local reduced basis method (local RB) in terms of accurately approximating the solution to PDEs with random inputs as well as quantifying the risk (CVaR in particular) of QoIs of the solution. We compare the performance of our local RB method with two popular adaptive 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

sparse grid approaches [22, 29]. The first method, denoted by ASG1 [29] from here on, uses locally supported piecewise multi-linear functions for interpolation. The second approach, denoted by ASG2 [22], uses global Lagrange polynomials for interpolation and a dimension adaptive strategy. We refer the readers to [42] for a comparison of the local RB method with the SROM approach presented in [25]. In all the comparisons presented herein, we fix the number of high-fidelity solves for all methods and compare the accuracy of their approximation. The comparison is particularly relevant in solving large-scale stochastic PDEs with limited computational budget. That is, we want to achieve a relatively low approximation error with only a small number of full PDE solves. It is important to point out that our reported number of full PDE solves does not include the Nλ = max{100, 10M } eigenvalue solves to build the surrogate model for the stability constants. Also, we did not exploit the fact that using a direct solver we could reduce the number of high-dimension matrix factorizations by a factor of M + 1. However, the advantage of local RB method is apparent from our results. 4.1. 1D Helmholtz equation with damped resonance Let D = (0, 1). We first present results for a damped Helmholtz problem, which is defined as

∂ − ∂x



∂u(x, ω) ν(x, ω) ∂x



− icτ u(x, ω) − τ 2 u(x, ω) = f (x), u(0, ω) = u(1, ω) = 0,

x ∈ D, a.s.

(32)

a.s.

where ν(x, ω) is a piecewise constant random field, c is a deterministic damping factor, τ is a deterministic angular frequency, and f is a deterministic source. Stability results for this type of problem can be found in [18]. We model the random field ν(x, ω) as ν(x, ω) = [0.3 + 1.2ξ1 (ω)] 1[0,0.5] (x) + [0.5 + 3.5ξ2 (ω)] 1[0.5,1] (x)

(33)

where ξ1 ∼ Beta(1, 3) and ξ2 ∼ Beta(3, 2) with correlation equal to 0.5. We set the angular frequency to τ = 4π. These parameters are chosen so that damped resonances 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

occur for certain realizations of the random field ν(x, ω). To see this, we plot the L2 (D) norm of the solution ku(·, ξ)kL2 (D) over the entire parameter domain Ξ with two different values of damping constants c = 0.01 and c = 0.2 in Figure 2. From Figure 2, it is clear that resonance occurs in three different regions of Ξ. It is well known that as c decreases towards zero, the solution norm at resonances increases unboundedly. This behavior makes Helmholtz-type problems particularly challenging for constructing reduced-order models in the parameter space. For this example, we fix c = 0.2, which is shown in Figure 2 (b).

(a)

(b)

Figure 2: The L2 (D) norm of u(·, ξ) over the parameter domain Ξ with (a) c = 0.01 (b) c = 0.2. Notice that resonances occur for certain realizations of ξ.

In Figure 3, we show the joint density of ξ as well as Nbkg = 2, 000 background samples of these parameters used in our local RB construction. The samples clearly concentrate on the left part of the domain. The other algorithmic parameters for this min init = 50 and Ntr = 5. The resulting reduced example are: Nspace = 200, N = 3, Ntr

basis dimension in this example is equal to 6, which is significantly smaller than Nspace . In Figure 4, we show the samples and corresponding Voronoi cells constructed using the local RB method in three different stages of Algorithm 1. We show the atoms with the corresponding training set Ξtr in the background. Notice how the adaptively selected atoms are placed close to resonance regions (Figure 2), where we expect the error u (ξ) to be large. It is interesting to notice that the atoms and training samples also cluster around the resonance area on the top right of the domain, even though this area is in the region of low probability density (Figure 3). The latter indicates how important 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

Figure 3: (a) Joint density of ξ. (b) Background samples of ξ.

it is to accurately capture the solution behavior around resonance regions in stochastic Helmholtz problems. 1

1

1

0.5

0.5

0.5

0

0

0.5

1

0

0

0.5

1

0

0

0.5

1

Figure 4: Atoms, partition and training set generated by the local RB method when Natom = 10, 25, 50, respectively.

In Figure 5, we compare the relative L2 (D) error eu (ξ) =

ku(ξ) − u(ξ)kL2 (D) ku(ξ)kL2 (D)

(34)

computed for the background samples of Figure 3. As can be observed, the mean relative L2 (D) error of the local RB method is almost three orders of magnitude smaller than that of the two adaptive sparse grid methods. The bad performance for the dimension adaptive method ASG2 is anticipated due to the localized nature of resonance regions in parameter space. We show the mean and variance of the solutions obtained using the different algorithms in Figure 6 and Figure 7, respectively. From all comparisons, we see the local RB method has better accuracy than the sparse grid algorithms given the same number of 23

0.4

1 0

0.3 -1

Distribution

log10 (eu ), mean and std

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

-2 -3

-5

0

0.2 0.1

ASG1 ASG2 Local RB

-4

ASG1 ASG2 Local RB

0

200 400 Number of full PDE solves

600

-6

-4

-2

0

log10 (eu )

(a)

(b)

Figure 5: Comparison of the relative L2 (D) error of the adaptive sparse grid methods and the local RB method. (a) Mean and standard deviation for different number of full PDE solves, (b) Distribution of the relative error when 600 full PDE solves are used to construct the surrogate models.

full PDE solves. The statistics of the solution can be accurately captured using the local RB surrogate constructed with only 60 full PDE solves (20 atoms). Finally, we investigate the performance of our adaptive algorithm for estimating the CVaR of the QoI G(ξ) = g(u(ξ)) = ku(·, ξ)kL2 (D) .

(35)

The functional g is H¨older continuous with α = 1 and K = 1. We aim to approximate CVaRβ (G) with β = 0.9. We refer to the QoI and CVaRβ computed from full PDE solutions the reference values and compare these values with the local RB method as well as the two sparse grid methods. Figure 8(a) shows the true QoI over the parameter domain Ξ, Figure 8(b) shows the 100, 000 Monte Carlo samples that we use to compute the risk, the red dots indicate the samples that contribute to the risk, i.e., the 10% largest samples of G. The cumulative distribution of G(ξ) as well as CVaRβ (G) are shown in Figure 8(c). In Figure 9(a), we plot the Voronoi tessellations from the first six steps of Algorithm 1 modified to approximate the risk CVaR0.9 (G). The color of each cell corresponds to the logarithm of the local error indicator δk := supϑ∈A E[ϑ1Ξk (ξ)u (ξ)α ] as defined in Section 3.3. At each iteration of Algorithm 1, the cell with the largest error, δk , is selected and a 24

#10 -3

5

6

#10 -3

5

4

0

#10 -3

6

2

2 real

imag

real

imag

-5 0

-10

0

-10 -2 Full ASG1 ASG2 Local RB

-15

-20

#10 -3

4

0

-5

0

0.5 x

-2 -15

-4

1

-6

Full ASG1 ASG2 Local RB

0

0.5 x

-20

1

0

0.5 x

(a)

-4 -6

1

0

0.5 x

1

(b)

Figure 6: Comparison of the mean of the PDE solution using ASG1 (blue), ASG2 (cyan) and local RB (red). The black line corresponds to the full PDE solution using 100, 000 Monte Carlo samples. (a) depicts the real and imaginary parts of the mean generated by surrogate models constructed using 60 full PDE solves, (b) depicts the real and imaginary parts of the mean generated by surrogate models constructed using 600 full PDE solves.

#10 -3

1

Full ASG1 ASG2 Local RB

0.8 0.6 0.4 0.2 0

#10 -3 Full ASG1 ASG2 Local RB

0.8

Variance

1

Variance

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.6 0.4 0.2

0

0.2

0.4

0.6

0.8

0

1

x

0

0.2

0.4

0.6

0.8

1

x

(a)

(b)

Figure 7: Comparison of the variance of the PDE solution using ASG1 (blue), ASG2 (cyan) and local RB (red). The black line corresponds to the full PDE solution using 100, 000 Monte Carlo samples. (a) depicts the real and imaginary parts of the variance generated by surrogate models constructed using 60 full PDE solves, (b) depicts the real and imaginary parts of the variance generated by surrogate models constructed using 600 full PDE solves.

25

CDF

1

CVaR

0.5

0

(a)

VaR 90%

CVaR = 0.05277

0

0.05 QoI

(b)

0.1

(c)

Figure 8: (a) The true QoI, (b) the parameter samples that contribute to CVaRβ (G) are shown in red, (c) the cumulative distribution of G(ξ) and CVaRβ (G).

new atom is found within the cell, which are depicted as the red dots. As more atoms are added, the local errors, δk , decrease resulting in an accurate approximation of CVaRβ (G). In Figure 9(b), we plot the relative error of the risk, ˆ − CVaRβ (G) CVaRβ (G) eR = , CVaRβ (G)

(36)

corresponding to the local RB method with different numbers of atoms. Clearly, as more atoms are added, the relative error drops.

0 -2

log10 (eR )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

-4 -6 -8

0

(a)

50 100 150 Number of atoms

200

(b)

Figure 9: (a) Voronoi tessellation from the first six iterations of Algorithm 1. The color of each cell corresponds the logarithm of the local error indicator, log10 (δk ). The red dots correspond to the selected atoms. (b) shows the relative error of risk obtained with the local RB method using varying numbers of atoms.

In Figure 10(a) and (c), we show G computed over the parameter space using local 26

PDE Solves

ASG1

ASG2

Local RB

60

4.04 × 10−1

4.35 × 10−1

3.78 × 10−4

2.84 × 10−2

600

1.70 × 10−2

6.94 × 10−8

Table 1: The relative error of CVaR0.9 (G) obtained using ASG1, ASG2 and local RB surrogate models with 60 and 600 PDE solves, respectively.

RB, ASG1 and ASG2 with 60 and 600 PDE solves, respectively. The comparisons of the corresponding CVaRβ (G) are shown in Figure 10(b) and (d) with the relative error eR shown in Table 1. The comparison clearly demonstrates that the CVaRβ (G) is captured much more accurately using the local RB method.

ASG2

ASG1

1

LocalRB

VaR

CVaR

90%

CDF

0.8

Full ASG1 ASG2 Local RB

0.6 0.4 0.2 0

0

0.02

0.04

(a) ASG1

ASG2

0.06 QoI

0.08

0.1

(b) 1

LocalRB

VaR

CVaR

90%

0.8 CDF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Full ASG1 ASG2 Local RB

0.6 0.4 0.2 0

(c) Figure 10:

0

0.02

0.04

0.06 QoI

0.08

0.1

(d)

(a) and (c): Comparison of the G(ξ) obtained using ASG1, ASG2 and

local RB surrogate models with 60 and 600 PDE solves, respectively. The (b) and (d) plot the corresponding risk generated by the surrogate models with 60 and 600 PDE solves, respectively.

27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

5. 2D advection-diffusion equation Let D = (0, 1)2 . We consider the 2D advection-diffusion problem −∇ · (κ(x, ω)∇u(x, ω)) + v(x, ω) · ∇u(x, ω) = f (x, ω)

x ∈ D a.s.

(37a)

u(x, ω) = 0

x ∈ Γd a.s.

(37b)

κ(x, ω)∇u(x, ω) · n = 0

x ∈ Γn a.s.

(37c)

where Γd := [0, 1] × {0} and Γn := ∂D \ Γd . The uncertainty in this example arises from diffusivity, advection and the forcing term. The diffusion coefficient κ is modeled as a random field using the approach in [26]. In particular, we characterize the diffusion coefficient as κ(x, ω) =

5 X 5 X

φi (x1 )φj (x2 )ξ5(i−1)+j (ω)

(38a)

i=1 j=1

ξk = 0.02 + 0.48βk

k = 1, 2, . . . , 25

(38b)

where φi (x) = Bi,5 (x) are Bernstein’s polynomials of order 5 and βk ∼ Beta(2, 5). The advection field is divergence free and is parametrized by two uniform random variables as     1 −x1 v(x, ω) = ξ26 + ξ27 0 x2 ξ26 ∼ U [10, 20]

ξ27 ∼ U [0, 10]

(39a) (39b)

Finally, the forcing term f is modeled by four Gaussians functions with uncertain magnitudes and located at [0.25, 0.25], [0.25, 0.75], [0.75, 0.25] and [0.75, 0.75]   −(x1 − 0.25)2 − (x2 − 0.25)2 ξ28 f (x, ω) = exp 0.5σ 2   −(x1 − 0.25)2 − (x2 − 0.75)2 + exp ξ29 0.5σ 2   −(x1 − 0.75)2 − (x2 − 0.25)2 + exp ξ30 0.5σ 2   −(x1 − 0.75)2 − (x2 − 0.75)2 + exp ξ31 . 0.5σ 2 The forcing magnitudes are chosen to be uniformly distributed, ξ28 , ξ31 ∼ U [0, 5]

ξ29 , ξ30 ∼ U [0, 10]. 28

(40)

For this example, we use the following algorithmic parameters: Nbkg = 5,000, Nspace = init min 1,089, N = 20, Ntr = 100 and Ntr = 5. Hence, the reduced basis dimension is equal to

52 in this case. We again compare the adaptive local RB method with ASG1 and ASG2. We compare the relative L2 (D) error (34) over 5,000 test samples of the parameters. As seen from Figure 11, the mean relative L2 (D) error of the local RB method is one order of magnitude smaller than that of the two adaptive sparse grid methods. Notice that the sharp computational advantage of the local RB method persists as more full PDE solves are included for all methods. 0.4

0.25

0.3

0.2 Distribution

eu , mean and std

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.2 0.1 0 -0.1

ASG1 ASG2 Local RB

0.15 0.1 0.05

320 640 960 1280 1600 Number of full PDE solves

(a)

0

-5

-4

-3

-2

-1

0

log10 (eu )

(b)

Figure 11: Comparison of the relative L2 (D) error of the adaptive sparse grid methods and the local RB method (a) Relative error decreases slowly as more full PDE solves are used to construct the surrogate models. The local RB method is one order of magnitude more accurate than the adaptive sparse grid methods. (b) Distribution of the relative error when 1, 600 full PDE solves are used to construct the surrogate models.

In Figure 12 we show estimates of the mean and variance of the solution for the different algorithms with 1,600 full PDE solves. All fields are computed using 50,000 Monte Carlo samples of the uncertain parameters. Since the surrogate approximation at each sample point is cheap to compute given that the surrogate model has already been constructed, statistics of the solution such as the mean and variance fields are much cheaper to obtain via a surrogate model. In Figure 12, it can be seen that the mean field seems to be well approximated by all surrogate models. However, the variance field is only accurately reproduced by the local RB method. The two adaptive sparse grid methods fail to deliver good approximations 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Mean Variance

ASG1

ASG2

Local RB

9.40 × 10−2

1.44 × 10−2

1.56 × 10−4

1.66 × 10−1

1.29 × 10−1

1.59 × 10−3

Table 2: Relative L2 (D) error of the predicated moments of u using surrogate models constructed with 1,600 full PDE solves.

to the variance field for the number of full PDE solves used. We also computed the

(a) Mean

(b) Variance Figure 12: The approximated mean and variance of u using Monte Carlo method with 50,000 independent parameter samples.

relative L2 (D) error of the mean and variance of the solution. These errors are shown in Table 2. We can see from the table that the mean and variance are approximated with much higher accuracy by the local RB surrogate. Finally, we use the local RB method and the two sparse grid methods to approximate the CVaRβ (ku(ξ)kL2 (D) ) with β = 0.9. Figure 13 shows the relative error of the risk obtained by the local RB method with varying numbers of atoms. Figure 14 shows the CDF of the QoI obtained using the different methods when 1,600 PDE solves (50 atoms 30

eR

ASG1

ASG2

Local RB

3.56 × 10−2

3.97 × 10−2

3.90 × 10−4

Table 3: The relative error of CVaR0.9 (G) obtained using ASG1, ASG2 and local RB surrogate models (constructed with 1,600 PDE solves).

for local RB) are used to construct the surrogate models. The risk is accurately captured by the local RB method with a relative error of less than 0.1%, whereas the relative error is 3.56% for ASG1 and 3.97% for ASG2, as shown in Table 3. -2

log10 (eR )

-2.5

-3

-3.5

Figure 13:

0

10

20 30 Number of atoms

40

50

The relative error of CVaR0.9 (ku(ξ)kL2 (D) ) obtained by the local RB

method with different number of atoms.

VaR CVaR

1

90% Full ASG1 ASG2 Local RB

0.8 CDF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.6 0.4 0.2 0

0

0.5

1

1.5

QoI

Figure 14:

Comparison of the CVaR0.9 (ku(ξ)kL2 (D) ) obtained using ASG1, ASG2

and local RB surrogate models.

31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

6. Conclusion In this work, we have presented an adaptive sampling approach for solving PDEs with random inputs based on a local reduced basis approximation. Our method relies on a set of parameter atoms to generate an implicit Voronoi partition of the parameter domain and employs a local reduced basis projection to approximate the solution within each cell. There are three main contributions: a novel local reduced basis method enriched with gradient information, a greedy adaptive sampling algorithm based on an effective a posteriori error indicator with a dynamic training set, and the application of the proposed framework to approximate measures of risk. In the numerical examples, we compared our proposed method with two adaptive sparse grid methods. We demonstrated that our method outperforms the others in terms of accurately approximating the PDE solution as well as the statistics of the solution when the PDE solution is highly sensitive to the random inputs. In addition, we showed that our method accurately approximates coherent risk measures evaluated at a quantity of interest depending on the PDE solution.

References [1] Philippe Artzner, Freddy Delbaen, Jean-Marc Eber, and David Heath. Coherent measures of risk. Mathematical finance, 9(3):203–228, 1999. [2] Ivo Babuska. Error-Bounds for Finite Element Method. 333:322–333, 1971. [3] Ivo Babuˇska, Fabio Nobile, and Raul Tempone. A stochastic collocation method for elliptic partial differential equations with random input data. SIAM Journal on Numerical Analysis, 45(3):1005– 1034, 2007. [4] Ivo Babuska, Ra´ ul Tempone, and Georgios E Zouraris. Galerkin finite element approximations of stochastic elliptic partial differential equations. SIAM Journal on Numerical Analysis, 42(2):800– 825, 2004. [5] Maxime Barrault, Yvon Maday, Ngoc Cuong Nguyen, and Anthony T Patera. An empirical interpolationmethod: application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus Mathematique, 339(9):667–672, 2004. [6] Fabian Bastin, Cinzia Cirillo, and Philippe L Toint. An adaptive monte carlo algorithm for computing mixed logit estimators. Computational Management Science, 3(1):55–79, 2006. [7] S´ ebastien Boyaval, Claude Le Bris, Tony Leli` evre, Yvon Maday, Ngoc Cuong Nguyen, and Anthony T Patera. Reduced basis techniques for stochastic problems. Archives of Computational methods in Engineering, 17(4):435–454, 2010.

32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[8] S´ ebastien Boyaval, Claude Le Bris, Yvon Maday, Ngoc Cuong Nguyen, and Anthony T Patera. A reduced basis approach for variational problems with stochastic parameters: Application to heat conduction with variable robin coefficient. Computer Methods in Applied Mechanics and Engineering, 198(41):3187–3206, 2009. [9] Tan Bui-Thanh, Karen Willcox, and Omar Ghattas. Model reduction for large-scale systems with high-dimensional parametric input space. SIAM Journal on Scientific Computing, 30(6):3270–3288, 2008. [10] Carsten Carstensen, Martin Eigel, Ronald HW Hoppe, and Caroline L¨ obhard. A review of unified a posteriori finite element error control. Numerical Mathematics: Theory, Methods and Applications, 5(04):509–558, 2012. [11] Peng Chen and Alfio Quarteroni. Accurate and efficient evaluation of failure probability for partial different equations with random input data. Computer Methods in Applied Mechanics and Engineering, 267:233–260, 2013. [12] Peng Chen, Alfio Quarteroni, and Gianluigi Rozza. A weighted reduced basis method for elliptic partial differential equations with random input data. SIAM Journal on Numerical Analysis, 51(6):3163–3185, 2013. [13] Peng Chen, Alfio Quarteroni, and Gianluigi Rozza. Comparison between reduced basis and stochastic collocation methods for elliptic problems. Journal of Scientific Computing, 59(1):187–216, 2014. [14] Peng Chen, Alfio Quarteroni, and Gianluigi Rozza. Reduced basis methods for uncertainty quantification. SIAM/ASA Journal on Uncertainty Quantification, 5(1):813–869, 2017. [15] Peng Chen and Christoph Schwab. Sparse-grid, reduced-basis bayesian inversion. Computer Methods in Applied Mechanics and Engineering, 297:84–115, 2015. [16] Tiangang Cui, Youssef M Marzouk, and Karen E Willcox. Data-driven model reduction for the bayesian solution of inverse problems. International Journal for Numerical Methods in Engineering, 102(5):966–990, 2015. [17] Manas K Deb, Ivo M Babuˇska, and J Tinsley Oden. Solution of stochastic partial differential equations using galerkin finite element techniques. Computer Methods in Applied Mechanics and Engineering, 190(48):6359–6372, 2001. [18] Leszek F. Demkowicz. Viscoelastic vibrations in fluid, analysis of a model problem. Technical Report TICAM REPORT 96-12, The University of Texas at Austin, March 1996. [19] Leszek F Demkowicz and Jay Gopalakrishnan. An overview of the discontinuous petrov galerkin method. In Recent Developments in Discontinuous Galerkin Finite Element Methods for Partial Differential Equations, pages 149–180. Springer, 2014. [20] Martin Drohmann and Kevin Carlberg. The romes method for statistical modeling of reducedorder-model error. SIAM/ASA Journal on Uncertainty Quantification, 3(1):116–145, 2015. [21] RV Field, M Grigoriu, and JM Emery. On the efficacy of stochastic collocation, stochastic galerkin, and stochastic reduced order models for solving stochastic problems. Probabilistic Engineering Mechanics, 41:60–72, 2015. [22] Baskar Ganapathysubramanian and Nicholas Zabaras. Sparse grid collocation schemes for stochastic

33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

natural convection problems. Journal of Computational Physics, 225(1):652–685, 2007. [23] Thomas Gerstner and Michael Griebel. Dimension–adaptive tensor–product quadrature. Computing, 71(1):65–87, 2003. [24] Roger G Ghanem and Pol D Spanos. Stochastic finite elements: a spectral approach. Courier Corporation, 2003. [25] Mircea Grigoriu. A method for solving stochastic equations by reduced order models and local approximations. Journal of Computational Physics, 231(19):6495–6513, 2012. [26] Mircea Grigoriu. Response statistics for random heterogeneous microstructures. SIAM/ASA Journal on Uncertainty Quantification, 2(1):252–275, 2014. [27] Drew P Kouri, Matthias Heinkenschloss, Denis Ridzal, and Bart G van Bloemen Waanders. A trust-region algorithm with adaptive stochastic collocation for pde optimization under uncertainty. SIAM Journal on Scientific Computing, 35(4):A1847–A1879, 2013. [28] Drew P Kouri, Matthias Heinkenschloss, Denis Ridzal, and Bart G van Bloemen Waanders. Inexact objective function evaluations in a trust-region algorithm for pde-constrained optimization under uncertainty. SIAM Journal on Scientific Computing, 36(6):A3011–A3029, 2014. [29] Xiang Ma and Nicholas Zabaras. An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations. Journal of Computational Physics, 228(8):3084–3113, 2009. [30] Xiang Ma and Nicholas Zabaras. An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations. Journal of Computational Physics, 229(10):3884–3915, 2010. [31] Fabio Nobile, Raul Tempone, and Clayton G Webster. An anisotropic sparse grid stochastic collocation method for partial differential equations with random input data. SIAM Journal on Numerical Analysis, 46(5):2411–2442, 2008. [32] Fabio Nobile, Ra´ ul Tempone, and Clayton G Webster. A sparse grid stochastic collocation method for partial differential equations with random input data. SIAM Journal on Numerical Analysis, 46(5):2309–2345, 2008. [33] Victor Y Pan and Zhao Q Chen. The complexity of the matrix eigenproblem. In Proceedings of the thirty-first annual ACM symposium on Theory of computing, pages 507–516. ACM, 1999. [34] Alfio Quarteroni, Gianluigi Rozza, and Andrea Manzoni. Certified reduced basis approximation for parametrized partial differential equations and applications. Journal of Mathematics in Industry, 1(1):1, 2011. [35] R Tyrrell Rockafellar and Stan Uryasev. The fundamental risk quadrangle in risk management, optimization and statistical estimation. Surveys in Operations Research and Management Science, 18(1):33–53, 2013. [36] R Tyrrell Rockafellar and Stanislav Uryasev. Optimization of conditional value-at-risk. Journal of risk, 2:21–42, 2000. [37] Gianluigi Rozza, Dinh Bao Phuong Huynh, and Anthony T Patera. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential

34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

equations. Archives of Computational Methods in Engineering, 15(3):229–275, 2008. [38] Alexander Shapiro, Darinka Dentcheva, and Andrzej Ruszczy´ nski. Lectures on stochastic programming: modeling and theory, Second edition. SIAM, 2014. [39] Michael D Shields. Adaptive monte carlo analysis for strongly nonlinear stochastic systems. Reliability Engineering & System Safety, 175:207–224, 2018. [40] Dongbin Xiu and Jan S Hesthaven. High-order collocation methods for differential equations with random inputs. SIAM Journal on Scientific Computing, 27(3):1118–1139, 2005. [41] Dongbin Xiu and George Em Karniadakis. The wiener–askey polynomial chaos for stochastic differential equations. SIAM journal on scientific computing, 24(2):619–644, 2002. [42] Zilong Zou, Drew Kouri, and Wilkins Aquino. An adaptive sampling approach for solving pdes with uncertain inputs and evaluating risk. In 19th AIAA Non-Deterministic Approaches Conference, page 1325, 2017.

35