Stochastic Galerkin reduced basis methods for parametrized linear elliptic partial differential equations

TL;DR

This paper introduces a stochastic Galerkin reduced basis (SGRB) method to efficiently estimate parameter-dependent statistics of elliptic PDE solutions, significantly reducing computational costs for multiple parameter queries.

Contribution

The paper develops a novel SGRB approach that combines stochastic Galerkin finite element methods with reduced basis techniques, enabling fast statistical output estimation for parametrized elliptic PDEs.

Findings

SGRB achieves similar accuracy with fewer basis functions compared to standard methods.

A single SGRB run estimates statistical outputs for new parameters, unlike Monte Carlo approaches.

Numerical tests demonstrate the efficiency of SGRB in a convection-diffusion-reaction problem.

Abstract

We consider the estimation of parameter-dependent statistics of functional outputs of elliptic boundary value problems (BVPs) with parametrized random and deterministic inputs. For a given value of the deterministic paremeter, a stochastic Galerkin finite element (SGFE) method can estimate the corresponding expectation and variance of a linear output at the cost of a single solution of a large block-structured linear system of equations. We propose a stochastic Galerkin reduced basis (SGRB) method as a means to lower the computational burden when statistical outputs are required for a large number of deterministic parameter queries. Our working assumption is that we have access to the computational resources necessary to set up such a reduced order model for a spatial-stochastic weak formulation of the parameter-dependent BVP. To derive an SGRB model, we project the spatial-stochastic weak solution of a parameter-dependent SGFE model onto a POD reduced basis generated from snapshots of SGFE solutions at representative values of the parameter. We propose residual-corrected estimates of the parameter-dependent expectation and variance of linear functional outputs and provide respective computable error bounds.~We test the SGRB method numerically for a convection-diffusion-reaction problem, choosing the convective velocity as a deterministic parameter and the parametrized reactivity field as a random input. Compared to a standard reduced basis model embedded in a Monte Carlo sampling procedure, the SGRB model requires a similar number of reduced basis functions to meet a given tolerance requirement. However, only a single run of the SGRB model suffices to estimate a statistical output for a new deterministic parameter value, while the standard reduced basis model must be solved for each Monte Carlo sample.

Figures

Tables

Table 1. Table 1: Model parameters of the test problem.

symbol	value	description
${\kappa }_{\_}0$	$-1000$	expected value of reactivity
$\sigma$	$200$	standard deviation of reactivity
$L$	$1$	correlation length
$K$	$5$	Karhunen-Loève truncation index
$\Omega$	${(-0.5,0.5)}^2$	spatial domain with boundary $\partial \Omega$
$𝒫$	${[-200,200]}^2$	deterministic parameter domain
${\Xi }_{\_}$	${[-\sqrt{3},\sqrt{3}]}^K$	random parameter domain

Table 2. Table 2: Discretization parameters of the test problem. The default values are used in the snapshot simulations. The reference values are used to assess the accuracy of the snapshot discretizations.

symbol	default	reference	description
$M_{\_}\text{FE}$	$225$	$961$	number of FE degrees of freedom
$N_{\_}{\xi }_{\_}$	$1024$	$16384$	number of MC samples of $y_{\_}\in {\Xi }_{\_}$
$d$	$2$	$3$	degree of SG polynomials
$M_{\_}\text{SG}$	$243$	$1024$	number of SG degrees of freedom: ${\left(d+1\right)}^K$

Equations

$$E[g]:=\frac{1}{N_{\_}{\xi }_{\_}}\sum _{\_}{n=1}^{N_{\_}{\xi }_{\_}}g({\xi }_{\_}{}^n),\underline{E}[g]:=\frac{1}{N_{\_}{\xi }_{\_}-1}\sum _{\_}{n=1}^{N_{\_}{\xi }_{\_}}g({\xi }_{\_}{}^n),V[g]:=\underline{E}[g^2]-\underline{E}[g]E[g],$$

$$u(y_{\_},{\mu }_{\_})\in X:a(u(y_{\_},{\mu }_{\_}),v;y_{\_},{\mu }_{\_})=f(v;y_{\_},{\mu }_{\_})\forall v\in X.$$

$$u^R(y_{\_},{\mu }_{\_})\in X^R:a(u^R(y_{\_},{\mu }_{\_}),v;y_{\_},{\mu }_{\_})=f(v;y_{\_},{\mu }_{\_})\forall v\in X^R.$$

$$X_{\_}^R(1)$$

$$X_{\_}^R(2)$$

$$X_{\_}^R(3)$$

$$X_{\_}^R(4)$$

$$u_{\_}(i)\in X$$

$$u_{\_}{(i)}^R\in X_{\_}^R(i)$$

$$r(\cdot )$$

$$r_{\_}(i)(\cdot )$$

$$\alpha (y_{\_},{\mu }_{\_}):=\underset{\_}{\inf }v\in X\setminus \{0\}\frac{a(v,v;y_{\_},{\mu }_{\_})}{\Vert v{\Vert }_{\_}{X}^2}\forall (y_{\_},{\mu }_{\_})\in {\Xi }_{\_}{}^{N_{\_}{\xi }_{\_}}\times \mathcal{P},$$

$$\Vert F{\Vert }_{\_}{X}^{\prime }:=\underset{\_}{\sup }v\in X\setminus \{0\}\frac{|F(v)|}{\Vert v{\Vert }_{\_}X}\forall F\in X^{\prime }.$$

$$\Vert u-u^R{\Vert }_{\_}X\le \frac{\Vert r{\Vert }_{\_}{X}^{\prime }}{\alpha }.$$

$$\alpha \Vert e{\Vert }_{\_}{X}^2\stackrel{\text{eqref}eq:coercivityStrong}{\le }a(e,e)\stackrel{\text{eqref}eq:strongPrimal}{=}f(e)-a(u^R,e)\stackrel{\text{eqref}eq:residualStrong}{=}r(e)\stackrel{\text{eqref}eq:dualNormStrong}{\le }\Vert r{\Vert }_{\_}{X}^{\prime }\Vert e{\Vert }_{\_}X.$$

$$|l(u)-l(u^R)+r(u_{\_}{(1)}^R)|$$

$$|l(u)-l(u^R)+r(u_{\_}{(1)}^R)|\stackrel{\text{eqref}eq:residualStrong}{=}|l(e)+f(u_{\_}{(1)}^R)-a(u^R,u_{\_}{(1)}^R)|\stackrel{\text{eqref}eq:strongPrimal}{=}|l(e)+a(e,u_{\_}{(1)}^R)|$$

$$\stackrel{\text{eqref}eq:residualStrongDual}{=}|r_{\_}(1)(e)|\stackrel{\text{eqref}eq:dualNormStrong}{\le }\Vert r_{\_}(1){\Vert }_{\_}{X}^{\prime }\Vert e{\Vert }_{\_}X\stackrel{\text{eqref}eq:errorReducedStrong}{\le }\frac{\Vert r_{\_}(1){\Vert }_{\_}{X}^{\prime }\Vert r{\Vert }_{\_}{X}^{\prime }}{\alpha }.$$

$$\displaystyle|E[l(u)]-E[l(u^{R})]+E[r(u_{\_}{(1)}^{R})]|\leq E\Bigg{[}\frac{\|r\|_{\_}{X^{\prime}}\|r_{\_}{(1)}\|_{\_}{X^{\prime}}}{\alpha}\Bigg{]}.$$

$$|E[l(u)]-E[l(u^{R})]+E[r(u_{\_}{(1)}^{R})]|\leq E[|l(u)-l(u^{R})+r(u_{\_}{(1)}^{R})|]\stackrel{{\scriptstyle\eqref{eq:functionalBoundsDualBased}}}{{\leq}}E\Bigg{[}\frac{\|r\|_{\_}{X^{\prime}}\|r_{\_}{(1)}\|_{\_}{X^{\prime}}}{\alpha}\Bigg{]}.$$

$$\displaystyle\big{|}V[l(u)]-V[l(u^{R})]+V[r(u_{\_}{(1)}^{R})]+\underline{E}[r(u_{\_}{(2)}^{R})]-\underline{E}[r(u_{\_}{(3)}^{R})]-E[r(u_{\_}{(4)}^{R})]\big{|}$$

$$\displaystyle\qquad\leq\underline{E}\Bigg{[}\frac{\|r\|_{\_}{X^{\prime}}^{2}\|r_{\_}{(1)}\|_{\_}{X^{\prime}}^{2}}{\alpha^{2}}\Bigg{]}+E\Bigg{[}\frac{\|r\|_{\_}{X^{\prime}}\|r_{\_}{(1)}\|_{\_}{X^{\prime}}}{\alpha}\Bigg{]}\underline{E}\Bigg{[}\frac{\|r\|_{\_}{X^{\prime}}\|r_{\_}{(1)}\|_{\_}{X^{\prime}}}{\alpha}\Bigg{]}$$

$$\displaystyle\qquad\qquad+\underline{E}\Bigg{[}\frac{\big{\|}r_{\_}{(2)}-r_{\_}{(3)}-\frac{N_{\_}{\xi_{\_}{}}-1}{N_{\_}{\xi_{\_}{}}}r_{\_}{(4)}\big{\|}_{\_}{X^{\prime}}\big{\|}r\big{\|}_{\_}{X^{\prime}}}{\alpha}\Bigg{]}.$$

$$\displaystyle\big{|}\underline{E}[l(u)^{2}]-\underline{E}[l(u)]E[l(u)]-\underline{E}[l(u^{R})^{2}]+\underline{E}[l(u^{R})]E[l(u^{R})]$$

$$\displaystyle\qquad+\underline{E}[r(u_{\_}{(1)}^{R})^{2}]-\underline{E}[r(u_{\_}{(1)}^{R})]E[r(u_{\_}{(1)}^{R})]+\underline{E}[r(u_{\_}{(2)}^{R})]-\underline{E}[r(u_{\_}{(3)}^{R})]-E[r(u_{\_}{(4)}^{R})]\big{|}$$

$$\displaystyle\qquad\qquad\leq\underbrace{\Big{|}\underline{E}[(l(u)-l(u^{R})+r(u_{\_}{(1)}^{R}))^{2}]\Big{|}}_{\_}{\stackrel{{\scriptstyle\eqref{eq:functionalBoundsDualBased}}}{{\leq}}\underline{E}\Big{[}\frac{\|r\|_{\_}{X^{\prime}}^{2}\|r_{\_}{(1)}\|_{\_}{X^{\prime}}^{2}}{\alpha^{2}}\Big{]}}$$

$$\displaystyle\qquad\qquad\qquad+\underbrace{\Big{|}\underline{E}[l(u)-l(u^{R})+r(u_{\_}{(1)}^{R})]\Big{|}}_{\_}{\stackrel{{\scriptstyle\eqref{eq:functionalBoundsDualBased}}}{{\leq}}\underline{E}\Big{[}\frac{\|r\|_{\_}{X^{\prime}}\|r_{\_}{(1)}\|_{\_}{X^{\prime}}}{\alpha}\Big{]}}\underbrace{\Big{|}E[l(u)-l(u^{R})+r(u_{\_}{(1)}^{R})]\Big{|}}_{\_}{\stackrel{{\scriptstyle\eqref{eq:functionalBoundsDualBased}}}{{\leq}}E\Big{[}\frac{\|r\|_{\_}{X^{\prime}}\|r_{\_}{(1)}\|_{\_}{X^{\prime}}}{\alpha}\Big{]}}$$

$$\displaystyle\qquad\qquad\qquad+\Big{|}\;\underline{E}[\;\underbrace{2(l(u^{R})-r(u_{\_}{(1)}^{R}))l(e)+r(u_{\_}{(2)}^{R})}_{\_}{\stackrel{{\scriptstyle\eqref{eq:residualStrong},\eqref{eq:strongPrimal},\eqref{eq:residualStrongDual}}}{{=}}-r_{\_}{(2)}(e)}\;]-\underline{E}[\;\underbrace{E[l(u^{R})-r(u_{\_}{(1)}^{R})]l(e)+r(u_{\_}{(3)}^{R})}_{\_}{\stackrel{{\scriptstyle\eqref{eq:residualStrong},\eqref{eq:strongPrimal},\eqref{eq:residualStrongDual}}}{{=}}-r_{\_}{(3)}(e)}\;]$$

$$\displaystyle\qquad\qquad\qquad\qquad-E[\;\underbrace{\underline{E}[l(u^{R})-r(u_{\_}{(1)}^{R})]l(e)+r(u_{\_}{(4)}^{R})}_{\_}{\stackrel{{\scriptstyle\eqref{eq:residualStrong},\eqref{eq:strongPrimal},\eqref{eq:residualStrongDual}}}{{=}}-r_{\_}{(4)}(e)}\;]\;\Big{|},$$

$$\big{|}\underline{E}[r_{\_}{(2)}(e)]-\underline{E}[r_{\_}{(3)}(e)]-E[r_{\_}{(4)}(e)]\big{|}\stackrel{{\scriptstyle\eqref{eq:dualNormStrong},\eqref{eq:errorReducedStrong}}}{{\leq}}\underline{E}\Bigg{[}\frac{\big{\|}r_{\_}{(2)}-r_{\_}{(3)}-\frac{N_{\_}{\xi_{\_}{}}-1}{N_{\_}{\xi_{\_}{}}}r_{\_}{(4)}\big{\|}_{\_}{X^{\prime}}\big{\|}r\big{\|}_{\_}{X^{\prime}}}{\alpha}\Bigg{]}.$$

$$\bar{u}({\mu }_{\_})\in \bar{X}:\bar{a}(\bar{u}({\mu }_{\_}),v;{\mu }_{\_})=\bar{f}(v;{\mu }_{\_})\forall v\in \bar{X}.$$

$${\bar{u}}^R({\mu }_{\_})\in {\bar{X}}^R:\bar{a}({\bar{u}}^R({\mu }_{\_}),v;{\mu }_{\_})=\bar{f}(v;{\mu }_{\_})\forall v\in {\bar{X}}^R.$$

$${\bar{u}}_{\_}(1)\in \bar{X}$$

$${\bar{u}}_{\_}(2)\in \bar{X}$$

$${\bar{u}}_{\_}(3)\in \bar{X}$$

$${\bar{u}}_{\_}^R(1)\in {\bar{X}}_{\_}^R(1)$$

$${\bar{u}}_{\_}^R(2)\in {\bar{X}}_{\_}^R(2)$$

$${\bar{u}}_{\_}^R(3)\in {\bar{X}}_{\_}^R(3)$$

Full text

To derive an SGRB model, we project the spatial-stochastic weak solution of a parameter-dependent SGFE model onto a POD reduced basis generated from snapshots of SGFE solutions at representative values of the parameter. We propose residual-corrected estimates of the parameter-dependent expectation and variance of linear functional outputs and provide respective computable error bounds. We test the SGRB method numerically for a convection-diffusion-reaction problem, choosing the convective velocity as a deterministic parameter and the parametrized reactivity field as a random input. Compared to a standard reduced basis model embedded in a Monte Carlo sampling procedure, the SGRB model requires a similar number of reduced basis functions to meet a given tolerance requirement. However, only a single run of the SGRB model suffices to estimate a statistical output for a new deterministic parameter value, while the standard reduced basis model must be solved for each Monte Carlo sample.

††journal: Computer Methods in Applied Mechanics and Engineering

{highlights}

ROMs to estimate parametrized statistics of functional outputs of elliptic PDEs.

Residual-corrected estimation of expectation and variance of statistical outputs.

Evaluation of stochastic Galerkin reduced basis (SGRB) model is sampling-free.

SGRB method achieves similar reduction of degrees of freedom as Monte Carlo RB.

1 Introduction

We consider linear elliptic boundary-value problems subject to a finite number of random and deterministic input parameters. Our goal is to compute the parameter-dependent expected value and variance of a functional output of interest. In this context, a reduced basis model provides a computationally inexpensive map between the deterministic input parameter and the corresponding output statistics. Moreover, it provides a computable a posteriori bound for the error between the reduced basis statistical estimate and a corresponding high-fidelity estimate. Reduced basis methods for linear elliptic boundary value problems with affinely parametrized deterministic data are well-understood [1, 2, 3]. We consider two approaches to include stochastic parameters:

Monte Carlo reduced basis (MCRB) method: The underlying equations are formulated weakly regarding the physical space, which means that the problem depends on both the deterministic and stochastic parameters. Monte Carlo sampling is used to estimate the parameter-dependent expected value and variance of a functional output of interest. An MCRB method for linear elliptic problems with error bounds for the expectation and variance of a linear functional output is derived in [4]. Improved error bounds are provided by [5]. Further advances are the introduction of a weighted error estimator [6] and the embedding in a multi-level procedure [7]. MCRB methods have also been applied to parabolic problems [8], saddle point problems [9], Bayesian inverse problems [10, 11, 12] and the assessment of rare events [13].

• 2.

Stochastic Galerkin reduced basis (SGRB) method: The underlying equations are formulated weakly regarding the spatial and stochastic dimensions, so that the problem depends on the deterministic parameters only. Parameter-dependent estimates of the expected value and variance of a functional output are obtained by direct integration of the reduced solution. The principle of SGRB methods is introduced in [14] for stochastic time-dependent incompressible Navier-Stokes problems, formulated weakly regarding the spatial and stochastic dimensions, with time acting as a parameter. Applications to linear dynamical systems are studied in [15, 16]. SGRB methods can be related to space-time reduced basis methods [17, 18], which rely on a weak formulation with respect to space and time. The idea of using SGRB methods to estimate parameter-dependent expected values is discussed in [19, section 8.2.1].

Our main contribution is the derivation of a stochastic Galerkin reduced basis method to compute residual-corrected parameter-dependent estimates of the expectation and variance of statistical outputs together with corresponding error bounds. A Monte Carlo reduced basis method will be used as a benchmark to assess the accuracy of the SGRB method and to compare the error bounds.

The creation of an SGRB model requires an underlying stochastic Galerkin finite element (SGFE) model or an equivalent high-fidelity Galerkin approximation. At least a few snapshots of the SGFE solution for different values of the deterministic parameter are needed to provide a suitable reduced basis. Moreover, evaluations of the SGFE linear and bilinear forms are necessary to derive the respective reduced-order Galerkin model and error bounds. We assume that the resources necessary for these computations are available in the setup phase. We do not discuss the associated costs explicitly as it is not the focus of this work but refer to [20] for a detailed analysis of the costs of stochatic Galerkin methods. We point out that the SGRB method presented in this paper is not a tool to reduce the computational burden associated with a single solution of an SGFE model as it is the objective in, e.g., [21, 22] using proper generalized decomposition, in [23, 24] using a rational Krylov method and a low-rank tensor approximation, respectively, and in [25, 26] using problem-tailored preconditioned iterative solvers. Instead, the SGRB approach targets the situation where a certain number of SGFE simulations are feasible in an expensive pre-processing step to create a reduced-order model (ROM) which can be evaluated cheaply for any given deterministic parameter. As an extreme scenario, one could imagine having supercomputer resources available in the setup phase whereas the ROM shall be evaluated on a microcontroller in real time. Therefore, SGRB models can be particularly useful in settings where statistical estimates are required for many values of the deterministic parameter, like in robust optimal control or the real-time exploration of parameter-dependent statistics.

Compared to Monte Carlo reduced basis methods, stochastic Galerkin reduced basis methods can substantially decrease the computational cost of estimating the expectation and variance for a given deterministic parameter value. The reason is that MCRB methods require sampling the reduced-order solution, which may lead to a large number of reduced-order simulations for a single query of the deterministic parameter. The same issue arises when a stochastic collocation method is applied instead of Monte Carlo [27]. SGRB methods overcome this drawback by evaluating the stochastic integrals in the offline stage, i.e., during the setup of the ROM. As a result, the cost of solving an SGRB model is similar to the cost of a single solution of a comparable RB model within an MC loop. At the same time, the SGRB model directly delivers a statistical estimate without sampling. Therefore, one can expect a speed-up factor in the order of magnitude of the number of MC samples.

2 Monte Carlo reduced basis method

We introduce a complete probability space $(\Theta,\mathcal{F},\mathbb{P}_{\_}{})$ consisting of a set $\Theta$ of elementary events, a $\sigma$-algebra $\mathcal{F}$ on $\Theta$ and a probability measure $\mathbb{P}_{\_}{}$ on $\mathcal{F}$. For $k=1,\dots,K$ with $K\in\mathbb{N}$, we define independent random variables $\xi_{\_}{k}\colon\Theta\rightarrow\Xi_{\_}{k}$, where $\Xi_{\_}{k}\subset\mathbb{R}$ is the image of $\xi_{\_}{k}$. We introduce respective probability distributions $\mathbb{P}_{\_}{\xi_{\_}{k}}$ and probability densities $p_{\_}{\xi_{\_}{k}}:\Xi_{\_}{k}\rightarrow\mathbb{R}^{+}$, so that $\mathbb{P}_{\_}{\xi_{\_}{k}}(B)=\int_{\_}Bp_{\_}{\xi_{\_}{k}}(y_{\_}{})\,\mathrm{d}y_{\_}{}=\mathbb{P}_{\_}{}(\xi_{\_}{k}^{-1}(B))$ for all $B$ in the Borel $\sigma$-algebra of $\Xi_{\_}{k}$. We collect the random variables in a random vector $\xi_{\_}{}\colon\Theta\rightarrow\Xi_{\_}{}$, where $\xi_{\_}{}=(\xi_{\_}{1},\dots,\xi_{\_}{K})^{T}$ and $\Xi_{\_}{}=\Xi_{\_}{1}\times\dots\times\Xi_{\_}{K}$, with joint distribution $\mathbb{P}_{\_}{\xi_{\_}{}}$ and density $p_{\_}{\xi_{\_}{}}:\Xi_{\_}{}\rightarrow\mathbb{R}^{+}$. We denote the expectation of any $\mathbb{P}_{\_}{\xi_{\_}{}}$-measurable function $g\colon\Xi_{\_}{}\rightarrow\mathbb{R}$ with density $p_{\_}{\xi_{\_}{}}$ by $\mathbb{E}[g]=\int_{\_}{\Xi_{\_}{}}g(y_{\_}{})\,\mathrm{d}\mathbb{P}_{\_}{\xi_{\_}{}}(y_{\_}{})=\int_{\_}{\Xi_{\_}{}}g(y_{\_}{})p_{\_}{\xi_{\_}{}}(y_{\_}{})\,\mathrm{d}y_{\_}{}$. We define the variance $\mathbb{V}[g]=\mathbb{E}[(g-\mathbb{E}[g])^{2}]$ for any $g\in L^{2}_{\_}{\xi_{\_}{}}(\Xi_{\_}{})$, where $L^{2}_{\_}{\xi_{\_}{}}(\Xi_{\_}{}):=\{\,v:\Xi_{\_}{}\rightarrow\mathbb{R}\mid\int_{\_}{\Xi_{\_}{}}v(y_{\_}{})^{2}p_{\_}{\xi_{\_}{}}(y_{\_}{})\,\mathrm{d}y_{\_}{}<\infty\,\}$.

We introduce a deterministic parameter $\mu_{\_}{}\in\mathcal{P}$ for some parameter domain $\mathcal{P}\subset\mathbb{R}^{P}$ with $P\in\mathbb{N}$. The final statistical outputs are scalar-valued $\mu_{\_}{}$-dependent functions representing approximations to the expectation and variance of a linear functional of a PDE solution.

We let $\{\xi_{\_}{}^{1},\dots,\xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\}$ be a set of independent copies of the random vector $\xi_{\_}{}$. For some $g\in L^{2}_{\_}{\xi_{\_}{}}(\Xi_{\_}{})$, we define Monte Carlo estimators

[TABLE]

for which $\mathbb{E}[E[g]]=\mathbb{E}[g]$ and $\mathbb{E}[V[g]]=\mathbb{V}[g]$ hold. We let $\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}:=\{y_{\_}{}^{1},\dots,y_{\_}{}^{N_{\_}{\xi_{\_}{}}}\}$ be a realization of $\{\xi_{\_}{}^{1},\dots,\xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\}$. A realization of a Monte Carlo estimate is obtained after substituting $\xi_{\_}{}^{n}$ by $y_{\_}{}^{n}$ in (1), assuming that $g(y_{\_}{})$ is computable for any $y_{\_}{}\in\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}$. In our approach, we view $N_{\_}{\xi_{\_}{}}$ as a discretization parameter and fix $\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}$ before we build the reduced basis model.

We focus on the case where $g$ depends on its argument via a discretized PDE problem. In the following, we provide a full-order model (section 2.1) and a reduced-order model (section 2.2) to approximate the solution of the PDE for a given realization of the deterministic and random input parameters. The computation of linear outputs and the corresponding statistics are described in section 2.3, together with the respective error bounds. For the separation of the computation into an expensive offline phase and a inexpensive online phase, we refer to [4, 5].

2.1 Monte Carlo finite element model

We use a stochastic strong form of a parametrized PDE problem with random data to formulate an MCFE model. Samples of the solution of the PDE problem are characterized by a separable Hilbert space $X$ with inner product $(\cdot,\cdot)_{\_}X$ and norm $\|\cdot\|_{\_}X$. We introduce a parametrized bilinear form $a(\cdot,\cdot;y_{\_}{},\mu_{\_}{}):X\times X\rightarrow\mathbb{R}$ as well as a parametrized linear form $f(\cdot;y_{\_}{},\mu_{\_}{}):X\rightarrow\mathbb{R}$. This allows a stochastic strong formulation of a linear elliptic PDE problem:

Problem 2.1** (MCFE model).**

For given $(y_{\_}{},\mu_{\_}{})\in\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\times\mathcal{P}$, find

[TABLE]

We assume that $a(\cdot,\cdot;y_{\_}{},\mu_{\_}{})$ is $(y_{\_}{},\mu_{\_}{})$-uniformly bounded and coercive on $X$ and that $f(\cdot;y_{\_}{},\mu_{\_}{})$ is $(y_{\_}{},\mu_{\_}{})$-uniformly bounded on $X$. Then Problem 2.1 has a unique solution for any given $(y_{\_}{},\mu_{\_}{})\in\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\times\mathcal{P}$ according to the Lax-Milgram lemma.

It is a usual premise in reduced basis methods that the discretization space of the underlying full-order model is assumed to be large enough to capture the solution with sufficient precision. The assessment of the error of the full-order solution with respect to the infinite-dimensional exact solution is delegated to the choice of the full-order discretization. In this spirit, we assume $X$ to be a sufficiently well-resolving finite element space with $M_{\_}\text{FE}$ degrees of freedom. Similarly, we assume $\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}$ to be a large enough sample set, so that the error associated with the MC sampling is sufficiently small. Consequently, the errors associated with the MC sampling and the FE discretization are not represented in our error estimates.

2.2 Monte Carlo reduced basis model

Let $X^{R}\subset X$ be an $R$-dimensional subspace. An example is given in section 4.1. A reduced-order model of Problem 2.1 is

Problem 2.2** (MCRB model).**

For given $(y_{\_}{},\mu_{\_}{})\in\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\times\mathcal{P}$, find

[TABLE]

The unique solvability of Problem 2.2 is a direct consequence of Problem 2.1 being well-posed and $X^{R}$ being a subspace of $X$.

2.3 Output statistics and error estimates

We derive residual-corrected RB approximations of MCFE estimates of the expectation and variance of linear outputs of the parametrized PDE problem. We provide error bounds converging quadratically in terms of residual norms. In particular, we transfer the dual-based error bounds of [5], considering the true expectation and variance of RB outputs, to the setting of [4], considering MC approximations of the expectation and variance. This requires an additional dual problem as well as a careful handling of different MC discretizations of the expected value, namely $E[\cdot]$ and $\underline{E}[\cdot]$ according to (1). Throughout this section, we assume the same dependency on the deterministic and stochastic parameters as in sections 2.1 and 2.2, but often omit an explicit notation of the parameter dependence for clarity.

We introduce a parametrized linear form $l(\cdot;y_{\_}{},\mu_{\_}{}):X\rightarrow\mathbb{R}$, assumed to be $(y_{\_}{},\mu_{\_}{})$-uniformly bounded on $X$. We complement Problems 2.1 and 2.2 with auxiliary sets of dual problems to allow for residual-corrected output computations. For brevity, we provide the definitions and problems all at once, which results in some interconnections between the following statements:

Definition 2.3**.**

Subspaces ${X^{R}_{\_}{(1)}},\dots,{X^{R}_{\_}{(4)}}$ and linear forms $l_{\_}{(1)},\dots,l_{\_}{(4)}$ are given by

[TABLE]

Problem 2.4** (dual MCFE models).**

For given $(y_{\_}{},\mu_{\_}{})\in\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\times\mathcal{P}$, find

[TABLE]

Problem 2.5** (dual MCRB models).**

For given $(y_{\_}{},\mu_{\_}{})\in\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\times\mathcal{P}$, find

[TABLE]

Definition 2.6**.**

A primal residual $r$ and dual residuals $r_{\_}{(1)},\dots,r_{\_}{(4)}$ are given by

[TABLE]

The following error bounds require a coercivity factor

[TABLE]

and a dual space $X^{\prime}$ of $X$, with norm

[TABLE]

For efficiency, an offline/online decomposition of the dual norms of the encountered functional is possible and the coercivity factor can be replaced by a strictly positive lower bound [4, 5].

First, we provide a bound for the error of the RB solution to Problem 2.2 with respect to the FE solution to Problem 2.1, point-wise in $\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\times\mathcal{P}$, see [5, Proposition 3.1]:

Theorem 2.7** (solution bound).**

For given $(y_{\_}{},\mu_{\_}{})\in\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\times\mathcal{P}$,

[TABLE]

Proof.

We define $e:=u-u^{R}$ and derive

[TABLE]

Dividing by $\alpha\|e\|_{\_}{X}$ gives the result. ∎

We approximate a parameter-dependent linear output $l(u(y_{\_}{},\mu_{\_}{});y_{\_}{},\mu_{\_}{})$ point-wise with a residual-corrected reduced-order approximation, see [5, Theorem 3.6]:

Theorem 2.8** (output bound).**

For given $(y_{\_}{},\mu_{\_}{})\in\Xi_{\_}{}^{N_{\_}{\xi_{\_}{}}}\times\mathcal{P}$,

[TABLE]

Proof.

We define $e:=u-u^{R}$ and reformulate

[TABLE]

∎

We approximate the MCFE estimate $E[l(u(\cdot,\mu_{\_}{});\cdot,\mu_{\_}{})]$ of the parameter-dependent expected linear output as follows, see [5, Corollary 4.2.]:

Theorem 2.9** (expected output bound).**

For given $\mu_{\_}{}\in\mathcal{P}$,

[TABLE]

Proof.

By Jensen’s inequality

[TABLE]

∎

Finally, we approximate the MCFE estimate $V[l(u(\cdot,\mu_{\_}{});\cdot,\mu_{\_}{})]$ of the parameter-dependent variance of the linear output, see [5, Theorem 4.5]:

Theorem 2.10** (output variance bound).**

For given $\mu_{\_}{}\in\mathcal{P}$,

[TABLE]

Proof.

By (1), defining $e:=u-u^{R}$,

[TABLE]

where

[TABLE]

∎

3 Stochastic Galerkin reduced basis method

In the following, we replace the Monte Carlo sampling by a stochastic Galerkin procedure. We provide a full-order model (section 3.1) and a reduced-order model (section 3.2) to approximate the stochastic solution of the PDE problem for a given realization of the deterministic input parameters. The computation of statistics of linear outputs are described in section 3.3, together with the respective error bounds. The computation can be separated into an expensive offline phase and a inexpensive online phase by standard means [2].

3.1 Stochastic Galerkin finite element model

We introduce a stochastic Galerkin discretization space $S\subset L^{2}_{\_}{\xi_{\_}{}}(\Xi_{\_}{})$. An example is given in section 5.2. We define the product space $\bar{X}:=S\otimes X$, which is a Hilbert space with inner product $(\cdot,\cdot)_{\_}{\bar{X}}:=(\cdot,\cdot)_{\_}{L^{2}_{\_}{\xi_{\_}{}}(\Xi_{\_}{},X)}$ and norm $\|\cdot\|_{\_}{\bar{X}}:=\|\cdot\|_{\_}{L^{2}_{\_}{\xi_{\_}{}}(\Xi_{\_}{},X)}$ in terms of the Bochner-type space $L^{2}_{\_}{\xi_{\_}{}}(\Xi_{\_}{},X):=\{\,v:\Xi_{\_}{}\rightarrow X\mid\int_{\_}{\Xi_{\_}{}}\|v(y_{\_}{})\|_{\_}{X}^{2}p_{\_}{\xi_{\_}{}}(y_{\_}{})\,\mathrm{d}y_{\_}{}<\infty\,\}$.

We derive a spatial-stochastic weak formulation by taking the expectation of (2). Defining $\bar{a}(w,v;\mu_{\_}{}):=\mathbb{E}[a(w,v;\cdot,\mu_{\_}{})]$ and $\bar{f}(v;\mu_{\_}{}):=\mathbb{E}[f(v;\cdot,\mu_{\_}{})]$ provides

Problem 3.1** (SGFE model).**

For given $\mu_{\_}{}\in\mathcal{P}$, find

[TABLE]

As a consequence of the coercivity and boundedness properties associated with Problem 2.1, the bilinear form $\bar{a}(\cdot,\cdot;\mu_{\_}{})$ is $\mu_{\_}{}$-uniformly bounded and coercive on $\bar{X}$ and the linear form $\bar{f}(\cdot;\mu_{\_}{})$ is $\mu_{\_}{}$-uniformly bounded on $\bar{X}$. Therefore, Problem 3.1 has a unique solution for any given $\mu_{\_}{}\in\mathcal{P}$ according to the Lax-Milgram lemma.

3.2 Stochastic Galerkin reduced basis model

We introduce an $R$-dimensional reduced space $\bar{X}^{R}\subset\bar{X}$. Suitable reduced spaces are provided in section 4.2. A reduced form of Problem 3.1 is given as follows:

Problem 3.2** (SGRB model).**

For given $\mu_{\_}{}\in\mathcal{P}$, find

[TABLE]

The subspace property $\bar{X}^{R}\subset\bar{X}$ and the well-posedness of Problem 3.1 imply that Problem 3.2 has a unique solution.

3.3 Output statistics and error estimates

We derive SGRB approximations of the expectation and variance of linear outputs together with error bounds with respect to the corresponding SGFE approximations. The variance can be interpreted in terms of quadratic outputs. We follow the ideas of [28, 29] to derive the respective error bounds.

We introduce a linear form $\bar{l}(v;\mu_{\_}{}):=\mathbb{E}[l(v;\cdot,\mu_{\_}{})]$, which is $\mu_{\_}{}$-uniformly bounded on $\bar{X}$. We complement the primal problem of section 3.1 with corresponding dual problems:

Problem 3.3** (dual SGFE models).**

For given $\mu_{\_}{}\in\mathcal{P}$, find

[TABLE]

Letting $\bar{X}^{R}_{\_}{(1)}\subset\bar{X}$, $\bar{X}^{R}_{\_}{(2)}\subset\bar{X}$ and $\bar{X}^{R}_{\_}{(3)}\subset\bar{X}$ be $R$-dimensional subspaces, we introduce the following set of reduced dual equations:

Problem 3.4** (dual SGRB models).**

For given $\mu_{\_}{}\in\mathcal{P}$, find

[TABLE]

The error bounds will be provided in terms of dual norms of residuals:

Definition 3.5**.**

Based on Problems 3.1, 3.2, 3.3 and 3.4,

[TABLE]

We define, for any $\mu_{\_}{}\in\mathcal{P}$, the coercivity factor

[TABLE]

and the continuity factor

[TABLE]

It is possible to replace these factors by efficiently computable upper and lower bounds [28]. We introduce the dual space $\bar{X}^{\prime}$ of $\bar{X}$ with norm

[TABLE]

We can derive the following error bound for the error in the reduced-order approximation of the solution:

Theorem 3.6** (solution bound).**

For given $\mu_{\_}{}\in\mathcal{P}$,

[TABLE]

Proof.

Analog to the proof of Lemma 2.7. ∎

In view of the definition of the weak linear form $\bar{l}$, we obtain the following bounds for the expected value and variance of the output:

Theorem 3.7** (expected output bound).**

For given $\mu_{\_}{}\in\mathcal{P}$,

[TABLE]

Proof.

Analog to the proof of Lemma 2.8. ∎

Theorem 3.8** (output variance bound).**

For given $\mu_{\_}{}\in\mathcal{P}$,

[TABLE]

Proof.

Setting $e=\bar{u}-\bar{u}^{R}$, we rewrite

[TABLE]

After expressing the variance in terms of expectations, the left-hand side of (18) becomes

[TABLE]

The final result follows from the following bounds:

[TABLE]

∎

4 Reduced spaces

We introduce candidate reduced spaces $X^{R}$ and $\bar{X}^{R}$ to be used in Problems 2.2 and 3.2, respectively. For simplicity, we focus on spaces generated by snapshot-based proper orthogonal decomposition (POD), but the theory of sections 2 and 3 does not dependent on this choice. We point out that the availability of computable error bounds also allows the use of greedy snapshot sampling [4, 5, 30].

The procedures described in this section can also be applied to create the dual reduced spaces encountered in Problems 2.5 and 3.4, by applying the POD to snapshots of the corresponding discretized dual solutions. The creation of the dual reduced spaces must follow a certain sequence because some of the discretized dual problems contain reduced primal and dual solutions on their right-hand sides. For instance, creating $\bar{X}^{R}_{\_}{(3)}$ from samples of $\bar{u}_{\_}{(3)}$ requires the availability of $\bar{X}^{R}$ and $\bar{X}^{R}_{\_}{(1)}$ due to the right-hand side of the discretized dual problem that defines $\bar{u}_{\_}{(3)}$, see Problem 3.3.

We motivate the POD spaces by corresponding continuous minimization problems. We discretize these minimization problems using quadrature [3, sections 6.4 and 6.5]. The discrete minimization problems can be solved using a weighted singular value decomposition of a snapshot matrix, based on [31]. Algorithm 1 provides a definition of the POD algorithm in terms of linear algebra, assuming $N$ snapshot vectors of length $M$. The algorithm is formulated in a way that allows a general snapshot weighting and the maximum possible number of output vectors. Actual implementations can benefit from using a simpler (diagonal) snapshot weighting matrix and assuming a small number of output vectors. Sections 4.1 and 4.2 describe how to generate the input to the algorithm in order to compute POD basis vectors from available FE or SGFE snapshots.

4.1 Spatial POD

We provide a POD of snapshots of the solution $u$ of Problem 2.1, resulting in a spatial POD space $X^{R}=\operatorname{span}(\varphi_{\_}{1},\dots,\varphi_{\_}{R})$ for $R\leq M_{\_}\text{FE}$. One can define a POD basis as a set of functions which solve the continuous minimization problems

[TABLE]

for $R=1,\dots,M_{\_}\text{FE}$. We approximate the double integral using Monte Carlo quadrature.

Concerning the Monte Carlo quadrature of the $\Xi_{\_}{}$-integral on the one hand, we already know that the reduced-order model will only be evaluated at the random parameter points $y_{\_}{}^{1},\dots,y_{\_}{}^{N_{\_}{\xi_{\_}{}}}$, because the Monte Carlo discretization of the final stochastic estimates is fixed from the beginning. We use exactly these points for the discretization of the POD minimization problem, too, because with this choice our reduced basis will be optimal in a mean-square sense with respect to approximating the finite element solution at $y_{\_}{}^{1},\dots,y_{\_}{}^{N_{\_}{\xi_{\_}{}}}$.

Concerning the Monte Carlo quadrature of the $\mathcal{P}$-integral on the other hand, our model should be able to estimate the output statistics reasonably well at any point in $\mathcal{P}$. Having no further information about how the reduced-order model will ultimately be used, we discretize the deterministic parameter domain using a training set $\mathcal{P}^{N_{\_}{\mu_{\_}{}}^{\text{train}}}=\{\mu_{\_}{}^{1},\dots,\mu_{\_}{}^{N_{\_}{\mu_{\_}{}}^{\text{train}}}\}$, with $\mu_{\_}{}^{1},\dots,\mu_{\_}{}^{N_{\_}{\mu_{\_}{}}^{\text{train}}}$ distributed independently and uniformly over $\mathcal{P}$. When testing the performance of the resulting reduced-order model, we use a different set of points in the parameter domain in order to verify the robustness of the model with respect to the deterministic parameter.

The Monte Carlo quadrature of the double integral in (23) finally results in a set of discretized POD minimization problems

[TABLE]

for $R=1,\dots,M_{\_}\text{FE}$. For the POD computation in terms of Algorithm 1, we set $N=N_{\_}{\xi_{\_}{}}N_{\_}{\mu_{\_}{}}^{\text{train}}$ and $M=M_{\_}\text{FE}$ and let $U^{(i-1)N_{\_}{\xi_{\_}{}}+j}_{\_}{}\in\mathbb{R}^{M}$ be the coefficient vector corresponding to the expansion of $u(y_{\_}{}^{i},\mu_{\_}{}^{j})\in X$ in a basis of $X$ for $i=1,\dots,N_{\_}{\xi_{\_}{}}$ and $j=1,\dots,N_{\_}{\mu_{\_}{}}^{\text{train}}$. By substituting the finite element basis expansions into (24), we find

[TABLE]

where ${\mathcal{M}}_{\_}{X}$ denotes the mass matrix corresponding to $X$ and $I_{\_}{N}$ denotes the $N\times N$ identity matrix. The output of Algorithm 1 is a POD basis matrix $\Phi=(\Phi_{\_}1,\dots,\Phi_{\_}{M_{\_}\text{FE}})\in\mathbb{R}^{M_{\_}\text{FE}\times M_{\_}\text{FE}}$. The $i$-th POD basis function $\varphi_{\_}{i}$ can be obtained from the $i$-th POD basis vector $\Phi_{\_}i$ via an expansion in the available basis of $X$, using the elements of $\Phi_{\_}i$ as expansion coefficients. Finally, an $R$-dimensional POD reduced space is given by $X^{R}=\operatorname{span}(\varphi_{\_}{1},\dots,\varphi_{\_}{R})$ for any $R=1,\dots,M_{\_}\text{FE}$ and the trivial space $X^{0}\subset X$ contains only the zero function.

4.2 Spatial-stochastic POD

We introduce a reduced basis space that can be used to derive a stochastic Galerkin reduced basis method. It employs a simultaneous reduction of the spatial and stochastic degrees of freedom of a stochastic Galerkin finite element discretization.

Spatial-stochastic POD reduced basis functions can be defined as solutions to a set of $\mathcal{P}$-continuous POD minimization problems

[TABLE]

for $R=1,\dots,M_{\_}\text{FE}M_{\_}\text{SG}$. A Monte Carlo quadrature of the $\mathcal{P}$-integral raises the issue of choosing the sample points. Using the same training set $\mathcal{P}^{N_{\_}{\mu_{\_}{}}^{\text{train}}}=\{\mu_{\_}{}^{1},\dots,\mu_{\_}{}^{N_{\_}{\mu_{\_}{}}^{\text{train}}}\}$ as in section 4.1 leads to discrete POD minimization problems

[TABLE]

for $R=1,\dots,M_{\_}\text{FE}M_{\_}\text{SG}$. Regarding the POD computation in terms of Algorithm 1, we set $N=N_{\_}{\mu_{\_}{}}^{\text{train}}$ and $M=M_{\_}\text{FE}M_{\_}\text{SG}$ and denote the stochastic Galerkin finite element coefficient vector of $\bar{u}(\mu_{\_}{}^{n})$ by $U^{n}_{\_}{}$. By substituting the stochastic Galerkin finite element basis expansions into (25), we obtain

[TABLE]

where $I_{\_}{N_{\_}{\mu_{\_}{}}^{\text{train}}}$ is the $N_{\_}{\mu_{\_}{}}^{\text{train}}\times N_{\_}{\mu_{\_}{}}^{\text{train}}$ identity matrix and ${\mathcal{M}}_{\_}{S}$ is the mass matrix containing the mutual $S$-inner products of the basis functions used to represent $S$. In view of Algorithm 1, the $i$-th POD basis function $\bar{\varphi}_{\_}{i}$ can be obtained from the $i$-th POD basis vector $\Phi_{\_}i$ via an expansion in the available basis of $\bar{X}$, using the elements of $\Phi_{\_}i$ as expansion coefficients. Finally, an $R$-dimensional POD reduced space is given by $\bar{X}^{R}=\operatorname{span}(\bar{\varphi}_{\_}{1},\dots,\bar{\varphi}_{\_}{R})$ for any $R=1,\dots,M_{\_}\text{FE}M_{\_}\text{SG}$ and the trivial space $\bar{X}^{0}\subset\bar{X}$ contains only the zero function.

5 Numerical experiments

We assess the provided error bounds and compare the accuracy of the MCRB and SGRB models in terms of computing the expectation and variance of a linear output for a convection-diffusion-reaction problem.

5.1 Example problem

Let $y_{\_}{}=(y_{\_}{1},\dots,y_{\_}{K})^{T}\in\Xi_{\_}{}\subset\mathbb{R}^{K}$ denote the value of a sample of a random parameter vector, $\mu_{\_}{}=(\mu_{\_}{1},\mu_{\_}{2})^{T}\in\mathcal{P}\subset\mathbb{R}^{2}$ the value of a deterministic parameter vector and $x_{\_}{}=(x_{\_}{1},x_{\_}{2})^{T}\in\Omega\subset\mathbb{R}^{2}$ the coordinate in the computational domain $\Omega$. We model the random input by a second-order random field with expected value $\kappa_{\_}{0}$ and separable exponential covariance $c(x_{\_}{})=\sigma^{2}\exp(-|x_{\_}{1}|/L-|x_{\_}{2}|/L)$, where $\sigma$ is the standard deviation and $L$ is the correlation length. We approximate the random field using a truncated Karhunen-Loève expansion $\kappa_{\_}{}(x_{\_}{};y_{\_}{})=\kappa_{\_}{0}+\sigma\sum_{\_}{k=1}^{K}\sqrt{\lambda_{\_}{k}}\kappa_{\_}{k}(x_{\_}{})y_{\_}{k},$ where $\lambda_{\_}{k}$ denote the eigenvalues of the corresponding eigenvalue problem, ordered decreasingly by magnitude, and $\kappa_{\_}{k}(x_{\_}{})$ denote respective eigenfunctions. The covariance function allows for an analytical solution of the eigenvalue problem [32]. We assume that the parameters of the Karhunen-Loève expansion originate from independent uniformly distributed random variables. The truncation index $K$ can be interpreted as a modeling parameter, because it enters the definition of the bilinear form. The governing equations of our example problem are provided as follows:

Problem 5.1** (spatial strong form).**

For any $(y_{\_}{},\mu_{\_}{})\in\Xi_{\_}{}\times\mathcal{P}$, find $u(\cdot;y_{\_}{},\mu_{\_}{})\colon\Omega\rightarrow\mathbb{R}$ such that

[TABLE]

The deterministic parameter vector $\mu_{\_}{}\in\mathcal{P}\subset\mathbb{R}^{2}$ can be interpreted as a spatially uniform convective velocity. The random parameter vector $y_{\_}{}\in\Xi_{\_}{}\subset\mathbb{R}^{K}$ enters via a parametrized random reactivity. A concrete instance of the example problem is determined by the model parameters given in table 1. The output of the example problem is given by $l(u(y_{\_}{},\mu_{\_}{});y_{\_}{},\mu_{\_}{})$, where

[TABLE]

In order to express the example problem in terms of the spatial weak form of Problem 2.1, we set

[TABLE]

with

[TABLE]

and

[TABLE]

A spatial weak form of Problem 5.1 is provided in terms of the standard infinite-dimensional Sobolev space $H^{1}_{\_}0(\Omega)$ as follows:

Problem 5.2** (spatial weak form).**

For given $(y_{\_}{},\mu_{\_}{})\in\Xi_{\_}{}\times\mathcal{P}$, find

[TABLE]

By taking the expectation and using the notation of section 3.1, a spatial-stochastic weak form is given by

Problem 5.3** (spatial-stochastic weak form).**

For given $\mu_{\_}{}\in\mathcal{P}$, find

[TABLE]

5.2 Discretization

The MCFE and SGFE discretizations (Problems 2.1 and 3.1) provide necessary links between the infinite-dimensional test problems (Problems 5.2 and 5.3) and the respective reduced-order models (Problems 2.2 and 3.2). In the following, we describe the computational details of the MCFE and SGFE discretizations of the test problem. Table 2 lists our choice of the relevant discretization parameters.

Finite element method

We derive an instance of the stochastic strong finite element problem by replacing $H^{1}_{\_}0(\Omega)$ in Problem 5.2 with a finite-dimensional subspace. In particular, we employ the space $X\subset H^{1}_{\_}0(\Omega)$ formed by continuous piecewise linear finite elements corresponding to a regular graded simplicial triangulation of $\Omega$, characterized by the number $M_{\_}\text{FE}$ of finite element degrees of freedom. We estimate the spatial discretization error using simulations on a finer reference triangulation as a substitute for the exact solution.

Monte Carlo method

We provide estimates of the expectation and variance by discretizing the respective stochastic integrals using Monte Carlo quadrature in the sense of section 2.1. To this end, we generate random samples $y_{\_}{}^{1},\dots,y_{\_}{}^{N_{\_}{\xi_{\_}{}}}\in\Xi_{\_}{}$ according to the distribution $\mathbb{P}_{\_}{\xi_{\_}{}}$ with a standard pseudorandom number generator. A reference simulation with a higher number of samples delivers an estimate of the sampling error.

Stochastic Galerkin method

Stochastic Galerkin methods estimate the expectation and variance by directly evaluating the respective stochastic integrals, given a stochastic Galerkin solution based on a finite-dimensional subspace of $L^{2}_{\_}{\xi_{\_}{}}(\Xi_{\_}{})$. In general, a stochastic Galerkin finite element method applied to a linear elliptic problem with a random elliptic coefficient leads to a large, block-structured system of linear algebraic equations. In our case, however, the underlying random variables $y_{\_}{1},\dots,y_{\_}{K}$ are independent and enter the bilinear form linearly, see (27). Under these conditions, it is possible to find a double-orthogonal polynomial basis which decouples the blocks in the system matrix [33, 34]. The resulting block-diagonal system of equations can be solved efficiently due to the relatively small bandwidth and the ability to treat the blocks in parallel. To define a suitable double-orthogonal basis, we start with $K$ spaces of possibly different dimensions, spanned by univariate Legendre polynomials over the interval $[-\sqrt{3},\sqrt{3}]$. We normalize the polynomials regarding the underlying uniform distribution and rotate the bases such that they consist of double-orthogonal univariate polynomials, as described in [33]. Finally, a tensor product of these univariate double-orthogonal polynomial bases forms a basis of an $M_{\_}\text{SG}$-dimensional subspace $S\subset L^{2}_{\_}{\xi_{\_}{}}(\Xi_{\_}{})$. In our experiments, we use the same polynomial degree $d$ in all directions. We assess the error associated with the choice of $d$ by comparing with a reference solution using a higher degree.

Reduced basis

The considered reduced-order models rely on POD spaces generated from snapshots of the underlying discretized solution. We choose $N_{\_}{\mu_{\_}{}}^{\text{train}}=64$ as the number of training samples in the deterministic parameter domain. Consequently, section 4 specifies the creation of the reduced spaces.

5.3 Results

Figure 1 presents the parameter-dependent output statistics obtained with the default parameter-dependent SGFE model. The crosses in Figure 1 correspond to the snapshot training parameter values provided by the pseudo random number generator. Additionally, Figure 1 shows the test parameter values that are used to assess the model.

The reduced basis estimates of the output statistics together with the respective error bounds are provided by Theorems 2.9 and 2.10 for the MCRB method and Theorems 3.7 and 3.8 for the SGRB method. First, we validate the error bounds for a single random realization of the deterministic parameter, marked by a square in Figure 1. The convergence regarding the number of reduced basis functions $R$ is presented in Figure 2. The error components of the underlying discretized solution are provided as a reference. Looking at the discretization errors only, we see that number of MC samples is sufficient to approximate the expectation but actually too small to balance the FE error in case of the variance. The SG error, on the other hand, is smaller than the FE error in all cases, which provides evidence that the stochastic Galerkin discretization of the stochastic domain is sufficiently fine. Concerning the reduced basis models, we observe that $R\approx 16$ reduced basis functions are sufficient to obtain reduced-order estimates which are on a par with the full-order estimates in all considered cases. The plots suggest that all error bounds converge at the same rates as the respective errors. This is useful, because it implies that efficiency of the error bounds does not become significantly worse when the number of reduced basis functions is increased.

We assess the convergence globally over $\mathcal{P}$ in order to confirm that the point-wise observation in the deterministic parameter space provided by Figure 2 is not a lucky coincidence. To this end, we employ an $L^{2}_{\_}{}(\mathcal{P})$-norm, approximated using Monte Carlo quadrature with $N_{\_}{\mu_{\_}{}}^{\text{test}}=64$ samples shown as circles in Figure 1. The convergence results are presented in Figure 3. Since we have averaged over the parameter space, the plots appear less random than the plots in Figure 2. The convergence of the estimates and the corresponding bounds correspond quite well. Moreover, the MCRB and SGRB methods perform similar in terms of accuracy per number of basis functions.

In Figures 2 and 3, it appears that the SGRB error bound over-estimates the actual error more severely (by 4 orders of magnitude) than the MCRB error bound (2 orders of magnitude). A closer inspection of the individual components of the error estimate reveals that for larger $R$ the lower-order term involving the continuity factor becomes responsible for the major portion of the error estimate. In particular, for $R=64$ at the parameter point corresponding to Figure 2, the terms on the right-hand side of Theorem 3.8 amount to approximately $4.2\cdot 10^{-12}$, $1.2\cdot 10^{-19}$ and $5.2\cdot 10^{-14}$, respectively.

6 Conclusion

We have observed that the SGRB method can deliver estimates of the expectation and variance of linear outputs with an accuracy similar to the MCRB method. Also, the SGRB error bounds regarding the expected value were very close to the respective MCRB bounds in our experiments. Concerning the variance, the presented SGRB bounds overestimate the error more severely than the available MCRB bounds, which opens opportunities for future improvement of the SGRB variance bound. Nevertheless, the MCRB and SGRB variance bounds both converge at the same order depending on the number of reduced basis functions. This behavior is reflected by the theory, which predicts the same order of convergence in terms of dual norms of residuals.

The MCRB statistical output estimates and error bounds require a Monte Carlo sampling of the reduced quantities point-wise in the random parameter domain. In our tests, 1024 samples were sufficient to balance the finite element error for the expectation, but an accurate prediction of the variance would require even more samples. The SGRB estimates and bounds, on the other hand, are obtained by an exact integration of the corresponding reduced basis expansions in the setup phase of the reduced-order model, and, thus, do not rely on Monte Carlo sampling. As a consequence, the primal and dual SGRB problems need to be solved only once for each new deterministic parameter. This benefit comes at the cost of a more expensive offline phase. In our tests, the SGRB and the MCRB methods achieved a similar reduction of degrees of freedom for a given error tolerance. As a consequence, the possible online speedup of SGRB methods compared to MCRB methods is in the order of magnitude of the number of Monte Carlo samples. This is particularly attractive in scenarios where evaluating the reduced order model shall be as inexpensive as possible for a large number of parameter queries while the offline costs are not a primary concern.