---
abstract: 'In this work we introduce the Multi-Index Stochastic Collocation method (MISC) for computing statistics of the solution of a PDE with random data. MISC is a combination technique based on mixed differences of spatial approximations and quadratures over the space of random data. We propose an optimization procedure to select the most effective mixed differences to include in the MISC estimator: [such optimization is a crucial step and allows us to build a method that, provided with sufficient solution regularity, is potentially more effective than other multi-level collocation methods already available in literature]{}. We then provide a complexity analysis that assumes decay rates of product type for such mixed differences, [showing that in the optimal case the convergence rate of MISC is only dictated by the convergence of the deterministic solver applied to a one dimensional problem.]{} We show the effectiveness of MISC with some computational tests, comparing it with other related methods available in the literature, such as the Multi-Index and Multilevel Monte Carlo, Multilevel Stochastic Collocation, Quasi Optimal Stochastic Collocation and Sparse Composite Collocation methods.'
address:
- 'CEMSE, King Abdullah University of Science and Technology, Thuwal 23955-6900, Saudi Arabia. '
- 'CSQI - MATHICSE, École Polytechnique Fédérale de Lausanne, Station 8, CH 1015, Lausanne, Switzerland'
- 'Dipartimento di Matematica “F. Casorati”, Università di Pavia, Via Ferrata 5, 27100 Pavia, Italy'
author:
- 'Abdul-Lateef Haji-Ali'
- Fabio Nobile
- Lorenzo Tamellini
- Raúl Tempone
title: 'Multi-Index Stochastic Collocation for random PDEs'
---

Uncertainty Quantification ,Random PDEs ,Multivariate approximation ,Sparse grids ,Stochastic Collocation methods ,Multilevel methods ,Combination technique.

41A10 , 65C20 , 65N30 , 65N05

Introduction
============

Uncertainty Quantification (UQ) is an interdisciplinary, fast-growing research area that focuses on devising mathematical techniques to tackle problems in engineering and natural sciences in which only a probabilistic description of the parameters of the governing equations is available, due to measurement errors, intrinsic non-measurability/non-predictability, or incomplete knowledge of the system of interest. In this context, “parameters” is a term used in broad sense to refer to constitutive laws, forcing terms, domain shapes, boundary and initial conditions, etc.

UQ methods can be divided into deterministic and randomized methods. While randomized techniques, which include the Monte Carlo sampling method, are essentially based on random sampling and ensemble averaging, deterministic methods proceed by building a surrogate of the system’s response function over the parameter space, which is then processed to obtain the desired information. Typical goals include computing statistical moments (expected value, variance, higher moments, correlations) of some quantity of interest of the system at hand, typically functionals of the state variables (forward problem), or updating the statistical description of the random parameters given some observations of the system at hand (inverse problem). In any case, multiple resolutions of the governing equations are needed to explore the dependence of the state variables on the random parameters. The computational method used should therefore be carefully designed to minimize the computational effort.

In this work, we focus on the case of PDEs with random data, for which both deterministic and randomized approaches have been extensively explored in recent years. As for the deterministic methods, we mention here the methods based on polynomial expansions computed either by global Galerkin-type projections [@ghanem.spanos:book; @lemaitre:book; @matthies.keese:galerkin; @todor.schwab:convergence; @xiu.karniadakis:wiener] or collocation strategies based on sparse grids (see e.g. [@babuska.nobile.eal:stochastic2; @b.griebel:acta; @nobile.tempone.eal:aniso; @xiu.hesthaven:high]), low-rank techniques [@khoromskij.schwab:tensor; @khoromskij.oseledets:tensor; @nouy:2008a; @ballani.gras:tensor] and reduced basis methods (see e.g. [@boyaval.patera:RBforPDE; @chen.eal:RBvsSG]). All these approaches have been found to be particularly effective when applied to problems with a moderate number of random parameters (low-dimensional probability space) and smooth response functions. Although significant effort has been expended on increasing the efficiency of such deterministic methods with respect to the number of random parameters (see, e.g., [@cohen.devore.schwab:nterm2], the seminal work on infinite dimensional polynomial approximation of elliptic PDEs with random coefficients), Monte Carlo-type approximations remain the primary choice for problems with non-smooth response functions and/or those that depend on a high number of random parameters, despite their slow convergence with respect to sample size.

A very promising methodology that builds on the classical Monte Carlo method and enhances its performance is offered by the so-called *Multilevel Monte Carlo* (MLMC). It was first proposed in [@Heinrich:MLMC] for applications in parametric integration and extended to weak approximation of stochastic differential equations in [@giles:MLMC], which also provided a full complexity analysis. Let $\{h_\ell\}_{\ell=0}^L$ be a (scalar) sequence of spatial/temporal resolution levels that can be used for the numerical discretization of the PDE at hand and $\{{F}_\ell\}_{\ell=0}^L$ be the corresponding approximations of the quantity of interest, and suppose that the final goal of the UQ analysis is to compute the expected value of ${F}$, ${{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}}$. While a classic Monte Carlo approach simply approximates the expected value by using an ensemble average over a sample of independent replicas of the random parameters, the MLMC method relies on the simple observation that, by linearity of expectation, $$\label{eq:MLMC-telescopic}
  {{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}} \approx {{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}_L\right]}} =
{{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}_0\right]}} + \sum_{\ell=1}^L{{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}_\ell - {F}_{\ell-1}\right]}},$$ and computes by independent Monte Carlo samplers each expectation in the sum. Indeed, if the discretization of the underlying differential model is converging with respect to the discretization level, $\ell$, the variance of $({F}_\ell-{F}_{\ell-1})$ will be smaller and smaller as $\ell$ increases, i.e., when the spatial/temporal resolution increases. Dramatic computational saving can thus be obtained by approximating the quantities ${{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}_\ell-{F}_{\ell-1}\right]}}$ with a smaller and smaller sample size, since most of the variability of ${F}$ will be captured with coarse simulations and only a few resolutions over the finest discretization levels will be performed. The MLMC estimator is therefore given by $$\label{eq:MLMC-estim}
  {{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}} \approx \sum_{\ell=0}^L \frac{1}{M_\ell}\sum_{m=1}^{M_\ell} \left({F}_{\ell}(\omega_{m,\ell}) - {F}_{\ell-1}(\omega_{m,\ell})\right),
  \quad \text{with } {F}_{-1}(\cdot)=0,$$ where $\omega_{m, \ell}$ are the i.i.d. replicas of the random parameters. The application of MLMC methods to UQ problems involving PDEs with random data has been investigated from the mathematical point of view in a number of recent publications, see e.g. [@bsz11; @bls13; @scheichl.charrier:MLMC; @scheichl.giles:MLMC; @mss12]. Recent works [@teckentrup.etal:MLSC; @van-wyk:MLSC; @hps13; @kss12] have explored the possibility of replacing the Monte Carlo sampler on each level by other quadrature formulas such as sparse grids or quasi-Monte Carlo quadrature, obtaining the so-called Multilevel Stochastic Collocation (MLSC) or Multilevel Quasi-Monte Carlo (MLQCM) methods. See also [@tesei:MCCV] for a related approach where the Multilevel Monte Carlo method is combined with a control variate technique.

The starting point of this work is instead the so-called Multi-Index Monte Carlo method (MIMC), recently introduced in [@abdullatif.etal:MultiIndexMC], that differs from the Multilevel Monte Carlo method in that the telescoping idea presented in equations - is applied to discretizations indexed by a multi-index rather than a scalar index, thus allowing each discretization parameter to vary independently of the others. Analogously to what done in [@teckentrup.etal:MLSC; @van-wyk:MLSC; @hps13] in the context of stochastic collocation, here we propose to replace the Monte Carlo quadrature with a sparse grid quadrature at each telescopic level, obtaining in our case the Multi-Index Stochastic Collocation method (MISC). In other words, MISC can be seen as a multi-index version of MLSC, or a stochastic collocation version of MIMC. From a slightly different perspective, MISC is also closely related to the combination technique developed for the solution of (deterministic) PDEs in [@Bungartz.Griebel.Roschke.ea:pointwise.conv; @b.griebel:acta; @Griebel.schneider.zenger:combination; @Hegland:combination; @griebel.harbrecht:tensor]; in this work, the combination technique is used with respect to both the deterministic and stochastic variables.

One key difference between the present work and [@teckentrup.etal:MLSC; @van-wyk:MLSC; @hps13] is that the number of problem solves to be performed at each discretization level is not determined by balancing the spatial and stochastic components of the error (based, e.g., on convergence error estimates), but rather suitably extending the knapsack-problem approach that we employed in [@nobile.eal:optimal-sparse-grids; @back.nobile.eal:optimal; @back.nobile.eal:lognormal] to derive the so-called Quasi-Optimal Sparse Grids method (see also [@griebel.knapek:optimized]). A somewhat analogous approach was proposed in [@bieri:sparse.tensor.coll], where the number of solves per discretization level is prescribed *a-priori* based on a standard sparsification procedure (we will give more details on the comparison between these different methods later on). [In this work, we provide a complexity analysis of MISC and illustrate its performance improvements, comparing it to other methods by means of numerical examples.]{}

The remainder of this paper is organized as follows. In Section \[s:problem-setting\], we introduce the problem to be solved and the approximation schemes that will be used. The Multi-Index Stochastic Collocation method is introduced in Section \[s:method\], and our main theorem detailing the complexity of MISC for a particular choice of an index set is presented in Section \[s:complexity\]. Finally, Section \[s:numerics\] presents some numerical tests, while Section \[s:conclusions\] offers some conclusions and final remarks. The Appendix contains the technical proof of the main theorem. Throughout the rest of this work we use the following notation:

-   ${\mathbb{N}}$ denotes the set of integer numbers including zero;

-   ${\mathbb{N}}_+$ denotes the set of positive integer numbers, i.e. excluding zero;

-   ${\mathbb{R}}_+$ denotes the set of positive real numbers, ${\mathbb{R}}_+ = \{ r \in {\mathbb{R}}: r > 0\}$;

-   ${\bm{1}}$ denotes a vector whose components are always equal to one;

-   ${{\ensuremath{{\boldsymbol e}}}}^\kappa_\ell$ denotes the $\ell$-th canonical vector in ${\mathbb{R}}^\kappa$, i.e., $({{\ensuremath{{\boldsymbol e}}}}^\kappa_\ell)_i = 1$ if $\ell=i$ and zero otherwise; however, for the sake of clarity, we often omit the superscript $\kappa$ when obvious from the context. For instance, if ${{{\ensuremath{{\boldsymbol v}}}}}\in {\mathbb{R}}^N$, we will write ${{{\ensuremath{{\boldsymbol v}}}}}- {{{\ensuremath{{\boldsymbol e}}}}}_1$ instead of ${{{\ensuremath{{\boldsymbol v}}}}}- {{{\ensuremath{{\boldsymbol e}}}}}_1^N$;

-   given ${{{\ensuremath{{\boldsymbol v}}}}}\in {\mathbb{R}}^N$, $|{{{\ensuremath{{\boldsymbol v}}}}}| = \sum_{n=1}^N v_n$, $\max({{{\ensuremath{{\boldsymbol v}}}}}) = \max_{n=1,\ldots N}v_n$ and $\min({{{\ensuremath{{\boldsymbol v}}}}}) = \min_{n=1,\ldots N}v_n$;

-   given ${{{\ensuremath{{\boldsymbol v}}}}}\in {\mathbb{R}}^N$ and $f:{\mathbb{R}}\to {\mathbb{R}}$, $f({{{\ensuremath{{\boldsymbol v}}}}})$ denotes the vector obtained by applying $f$ to each component of ${{{\ensuremath{{\boldsymbol v}}}}}$, $f({{{\ensuremath{{\boldsymbol v}}}}}) = [f(v_1),f(v_2),\cdots,f(v_N)] \in {\mathbb{R}}^N$;

-   given ${{{\ensuremath{{\boldsymbol v}}}}}, {{{\ensuremath{{\boldsymbol w}}}}}\in {\mathbb{R}}^N$, the inequality ${{{\ensuremath{{\boldsymbol v}}}}}> {{{\ensuremath{{\boldsymbol w}}}}}$ holds true if and only if $v_n > w_n$ $\forall n=1,\ldots,N$.

-   given ${{{\ensuremath{{\boldsymbol v}}}}}\in {\mathbb{R}}^D$ and ${{{\ensuremath{{\boldsymbol w}}}}}\in {\mathbb{R}}^N$, $[{{{\ensuremath{{\boldsymbol v}}}}}, {{{\ensuremath{{\boldsymbol w}}}}}] = (v_1,\ldots,v_D,w_1,\ldots, w_N) \in {\mathbb{R}}^{D+N}$.

Problem setting {#s:problem-setting}
===============

Let $\mathscr{B} \subset {\mathbb{R}}^{d}$, $d=1,2,3$, be an open hyper-rectangular domain (referred to hereafter as the “physical domain”) and let ${{{\ensuremath{{\boldsymbol y}}}}}=(y_1,y_2,\ldots,y_N)$ be a $N$-dimensional random vector whose components are mutually independent and uniformly distributed random variables with support $\Gamma_n \subset {\mathbb{R}}$ and probability density function $\rho_n(y_n) = \frac{1}{|\Gamma_n|}$. Denoting $\Gamma = \Gamma_1 \times \Gamma_2 \ldots \times \Gamma_N$ (referred to hereafter as the “stochastic domain” or “parameter space”) and by $\sigma_B(\Gamma)$ the Borel $\sigma$-algebra over $\Gamma$, $\rho({{{\ensuremath{{\boldsymbol y}}}}})d{{{\ensuremath{{\boldsymbol y}}}}}= \prod_{n=1}^N \rho_n(y_n) dy_n$ is therefore a probability measure on $\Gamma$, due to the independence of $y_n$, and $(\Gamma,\sigma_B(\Gamma),\rho({{{\ensuremath{{\boldsymbol y}}}}})d{{{\ensuremath{{\boldsymbol y}}}}})$ is a complete probability space. Consider the following generic PDE, together with the assumption stated next:

\[pb:strong\_form\_yy\] Find $u:\mathscr{B} \times \Gamma \rightarrow {\mathbb{R}}$ such that for $\rho$-almost every ${{{\ensuremath{{\boldsymbol y}}}}}\in \Gamma$ $$  \begin{cases}
    \mathcal{L}(u;{{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}}) = {\mathscr{F}}({{{\ensuremath{{\boldsymbol x}}}}})  & {{{\ensuremath{{\boldsymbol x}}}}}\in \mathscr{B},\\
    u({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})=0  & {{{\ensuremath{{\boldsymbol x}}}}}\in \partial \mathscr{B}.
  \end{cases}$$

\[assump:wellposed\] Problem \[pb:strong\_form\_yy\] is well posed in some Hilbert space $V$ for $\rho$-almost every ${{{\ensuremath{{\boldsymbol y}}}}}\in \Gamma$.

The solution of Problem \[pb:strong\_form\_yy\] can be seen as an $N$-variate Hilbert-space valued function $u({{{\ensuremath{{\boldsymbol y}}}}}) : \Gamma \rightarrow V$. The random variables, $y_n$, can represent scalar values whose exact value is unknown, or they can stem from a spectral decomposition of a random field, like a Karhunen-Loève or Fourier expansion, possibly truncated after a finite number of terms, see, e.g., [@babuska.nobile.eal:stochastic2; @back.nobile.eal:lognormal]. It is also useful to introduce the Bochner space $
 L^2_\rho(\Gamma; V) =
 \left\{ u: \Gamma \rightarrow V \mbox{ strongly measurable s.t. } \int_{\Gamma} \| u({{{\ensuremath{{\boldsymbol y}}}}}) \|_{V}^2 \rho({{{\ensuremath{{\boldsymbol y}}}}})d{{{\ensuremath{{\boldsymbol y}}}}}< \infty\right\}.
$ Finally, given some functional of the solution $u$, $\Theta: V \rightarrow {\mathbb{R}}$, we denote by ${F}: \Gamma \rightarrow {\mathbb{R}}$ the $N$-variate real-valued function assigning to each realization ${{{\ensuremath{{\boldsymbol y}}}}}\in \Gamma$ the corresponding value of $\Theta[u]$ (quantity of interest), i.e., ${F}({{{\ensuremath{{\boldsymbol y}}}}}) = \Theta[u(\cdot,{{{\ensuremath{{\boldsymbol y}}}}})]$, and we aim at estimating its expected value, $${{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}} = \int_\Gamma {F}({{{\ensuremath{{\boldsymbol y}}}}}) \rho({{{\ensuremath{{\boldsymbol y}}}}}) d{{{\ensuremath{{\boldsymbol y}}}}}.$$

\[example:ell\] As a motivating example, consider the following elliptic problem: find $u:\mathscr{B} \times \Gamma \rightarrow {\mathbb{R}}$ such that for $\rho$-almost every ${{{\ensuremath{{\boldsymbol y}}}}}\in \Gamma$ $$\label{eq:ell_PDE_strong_yy_example}
  \begin{cases}
    -{\textnormal{div}}( a({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})\nabla u({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}}) ) ={\mathscr{F}}({{{\ensuremath{{\boldsymbol x}}}}})  & {{{\ensuremath{{\boldsymbol x}}}}}\in \mathscr{B},\\
    u({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})=h({{{\ensuremath{{\boldsymbol x}}}}}) & {{{\ensuremath{{\boldsymbol x}}}}}\in\partial \mathscr{B},
  \end{cases}$$ holds, where ${\textnormal{div}}$ and $\nabla$ denote differentiation with respect to the physical variables, ${{{\ensuremath{{\boldsymbol x}}}}}$, only, and the function $a: \mathscr{B} \times \Gamma \rightarrow {\mathbb{R}}$ is bounded away from $0$ and $\infty$, i.e., there exist two constants, $a_{min}, a_{max}$, such that $$\label{eq:bounding_a}
0< a_{min} \leq a({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}}) \leq a_{max} <\infty,
\quad \forall {{{\ensuremath{{\boldsymbol x}}}}}\in \mathscr{B} \mbox{ and for } \rho\mbox{-almost every } {{{\ensuremath{{\boldsymbol y}}}}}\in \Gamma.$$ This boundedness condition guarantees that Assumption \[assump:wellposed\] is satisfied, i.e. the equation is well posed for $\rho$-almost every ${{{\ensuremath{{\boldsymbol y}}}}}\in \Gamma$, thanks to a straightforward application of the Lax-Milgram lemma; moreover, the equation is well posed in $L^2_\rho(\Gamma; V)$, where $V$ is the classical Sobolev space $H^1_0(\mathscr{B})$, see, e.g., [@babuska.nobile.eal:stochastic2]. This is the example we will focus on in Section \[s:numerics\], where we will test numerically the performance of the Multi-Index Stochastic Collocation method that we will detail in Section \[s:method\].

The method that we present in the following sections can be also applied to more general problems than Problem \[pb:strong\_form\_yy\] in which the forcing terms, boundary conditions and possibly domain shape are also modeled as uncertain; the extension to time-dependent problems with uncertain initial conditions is also straightforward. Other probability measures can also be considered; the very relevant case in which the random variables, $y_n$, are normally distributed is an example.

As will be clearer in a moment, the methodology we propose uses tensorized solvers for deterministic PDEs. Although for ease of exposition we have assumed that the spatial domain, $\mathscr{B}$, is a hyper-rectangle, it is important to remark that the methodology proposed in this work can also be applied to non hyper-rectangular domains: this can be achieved by introducing a mapping from a reference hyper-rectangle to the generic domain of interest (with techniques such as those proposed in the context of Isogeometric Analysis [@hughes:IGA] or Transfinite Interpolation [@gordon.hall:transfinite.interp]) or by a Domain Decomposition approach [@quarteroni1999domain] if the domain can be obtained as a union of hyper-rectangles.

Approximation along the deterministic and stochastic dimensions {#s:interpolation}
---------------------------------------------------------------

In practice, we can only access the value of ${F}$ via a numerical solver yielding a numerical approximation of the solution $u$ of Problem \[pb:strong\_form\_yy\], which depends on a set of $D$ discretization parameters, such as the mesh-size, the time-step, the tolerances of the numerical solvers, and others, which we denote by $h_i, i=1,\ldots,D$; we remark that in general $D$, the number of parameters, might be different from $d$, the number of spatial dimensions. For each of those parameters, we introduce a sequence of discretization levels, $h_{i,\alpha}, \alpha=1,2,\ldots$, and for each multi-index ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}^D_+$, we denote by $u^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})$ the approximation of $u$ obtained from setting $h_i = h_{i,\alpha_i}$, with the implicit assumption that $u^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}}) \to u({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})$ as $\min_{1\leq i \leq D} \alpha_i \to \infty$ for $\rho$-almost every ${{{\ensuremath{{\boldsymbol y}}}}}\in \Gamma$; similarly, we also write ${F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}({{{\ensuremath{{\boldsymbol y}}}}}) = \Theta[u^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}(\cdot,{{{\ensuremath{{\boldsymbol y}}}}})]$. For instance, we could solve the problem stated in Example \[example:ell\] by a finite differences scheme with grid-sizes $h_{i,\alpha_i} = h_02^{-\alpha_i}$ in direction $i=1,\ldots,D$, for some $h_0 > 0$.

The discretization of ${F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$ over the random parameter space $\Gamma$ will consist of a suitable linear combination of tensor interpolants over $\Gamma$ based on Lagrangian polynomials. Observe that this approach is sound only if ${F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$ is at least a continuous function over $\Gamma$ (the smoother ${F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$ is, the more effective the Lagrangian approximation will be); for instance, for the problem stated in Example \[example:ell\], it can be shown under moderate assumptions on $a({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})$ that ${F}$ and ${F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$ are ${{{\ensuremath{{\boldsymbol y}}}}}$-analytic, see, e.g., [@back.nobile.eal:optimal; @cohen.devore.schwab:nterm2]; we will return to this point in Section \[s:numerics\].

To derive a generic tensor Lagrangian interpolation of ${F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$, we first introduce the set ${\mathcal{C}}^0(\Gamma_n)$ of real-valued continuous functions over $\Gamma_n$, and the subspace of polynomials of degree at most $q$ over $\Gamma_n$, ${\mathbb{P}}^q(\Gamma_n) \subset {\mathcal{C}}^0(\Gamma_n)$. Next, we consider a sequence of univariate Lagrangian interpolant operators in each dimension $Y_n$, i.e., $\{\mathscr U_n^{m(\beta_n)}\}_{\beta_n \in {\mathbb{N}}_+}$, where we refer to the value $\beta_n$ as the “interpolation level”. Each interpolant is built over a set of $m(\beta_n)$ collocation points, $\mathscr{H}_n^{m(\beta_n)} = \{y_n^1, y_n^2 \ldots y_n^{m(\beta_n)}\} \subset \Gamma_n$, where $m$ is a strictly increasing function, with $m(0)=0$ and $m(1)=1$, that we call the “level-to-nodes function”; thus, the interpolant yields a polynomial approximation, $$\mathscr U^{m(\beta_n)}_n : {\mathcal{C}}^0(\Gamma_n) \rightarrow {\mathbb{P}}^{m(\beta_n)-1}(\Gamma_n), \qquad
\mathscr U^{m(\beta_n)}_n[f](y_n) = \sum_{j=1}^{m(\beta_n)}
                    \left( f(y_n^j) \prod_{k=1,k \neq j}^{m(\beta_n)} \frac{y_n - y_n^k}{y^j_n - y_n^k} \right),$$ with the convention that $\mathscr U^{0}_n[f] = 0 \,\, \forall f \in {\mathcal{C}}^0(\Gamma_n)$. The $N$-variate Lagrangian interpolant can then be built by a tensorization of univariate interpolants: denote by ${\mathcal{C}}^0(\Gamma)$ the space of real-valued $N$-variate continuous functions over $\Gamma$ and by ${\mathbb{P}}^{{{\ensuremath{{\boldsymbol q}}}}}(\Gamma) = \bigotimes_{n=1}^N {\mathbb{P}}^{q_n}(\Gamma_n)$ the subspace of polynomials of degree at most $q_n$ over $\Gamma_n$, with ${{{\ensuremath{{\boldsymbol q}}}}}= (q_1,\ldots,q_N) \in {\mathbb{N}}^N$, and consider a multi-index ${{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}\in {\mathbb{N}}_+^N$ assigning the interpolation level in each direction, $y_n$; the multivariate interpolant can then be written as $$\begin{aligned}
\mathscr{U}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} : {\mathcal{C}}^0(\Gamma) \rightarrow {\mathbb{P}}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}) -{\bm{1}}}(\Gamma), \qquad
\mathscr{U}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}[{F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}] ({{\ensuremath{{\boldsymbol y}}}})
&= \left( \mathscr{U}^{m(\beta_1)}_1 \otimes \cdots \otimes \mathscr{U}^{m(\beta_N)}_N \right) [{F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}]({{\ensuremath{{\boldsymbol y}}}}).$$

The set of collocation points needed to build the tensor interpolant $\mathscr{U}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}[u] ({{\ensuremath{{\boldsymbol y}}}})$ is the tensor grid $\mathscr{T}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} = \times_{n=1}^N \mathscr{H}_n^{m(\beta_n)}$ with cardinality $ \#\mathscr{T}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} = \prod_{n=1}^N m(\beta_n)$. Observe that the Lagrangian interpolant immediately induces an $N$-variate quadrature formula, $\mathscr{Q}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} : {\mathcal{C}}^0(\Gamma) \rightarrow {\mathbb{R}}$, $$\begin{aligned}
& \mathscr{Q}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}[{F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}]
= {{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{\mathscr{U}}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}[{F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}]({{\ensuremath{{\boldsymbol y}}}})\right]}}
= \sum_{j=1}^{\# \mathscr{T}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} } {F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}(\widehat{{{{\ensuremath{{\boldsymbol y}}}}}}_j) \varpi_j,\end{aligned}$$ where $\widehat{{{{\ensuremath{{\boldsymbol y}}}}}}_j \in \mathscr{T}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}$ and the quadrature weights $\varpi_j$ are the expected values of the Lagrangian polynomials centered in $\widehat{{{{\ensuremath{{\boldsymbol y}}}}}}_j$, which can be computed exactly for most of the common interpolation knots and probability measures of the random variables.

It is recommended that the collocation points $\mathscr{H}_n^{m(\beta_n)}$ to be used in each direction are chosen according to the underlying probability measure, $\rho(y_n)dy_n$, to ensure good approximation properties of the interpolant and quadrature operators, $\mathscr{U}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}$ and $\mathscr{Q}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}$. Common choices are Gaussian quadrature points like Gauss-Legendre for uniform measures or Gauss-Hermite for Gaussian measures, cf. e.g., [@trefethen:book2013], which are however *not nested*, i.e., $\mathscr{H}_n^{m(\beta_n)} \not\subset \mathscr{H}_n^{m(\beta_n+1)}$. This means that they are not optimal for successive refinements of the interpolation/quadrature, and we will not consider them in this work. Instead, we will work with *nested* collocation points, and specifically with Clenshaw-Curtis points [@nobile.eal:optimal-sparse-grids; @trefethen:comparison], that are a classical choice for the uniform measure that we are considering here; other choices of nested points are available for uniform random variables, e.g., the Leja points [@nobile.eal:optimal-sparse-grids; @Chkifa:leja], whose performance is somehow equivalent to that of Clenshaw-Curtis for quadrature purposes, see [@nobile.etal:leja; @narayan:Leja]. Clenshaw-Curtis points are defined as $$\label{eq:CC_def}
  y_n^j=\cos\left(\frac{(j-1) \pi}{m(i_n)-1}\right), \quad 1\leq j \leq m(i_n),$$ together with the following level-to-nodes relation, $m(i_n)$, that ensures their nestedness: $$\label{eq:m_CC}
m(0)=0,\,\, m(1)=1,\,\, m(i_n)=2^{i_n-1}+1.$$ We conclude this section by introducing the following operator norm, which acts as a “Lebesgue constant” from ${\mathcal{C}}^0(\Gamma)$ to $L^2_\rho(\Gamma)$: $$\label{eq:leb_def}
\mathbb{M}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} = \prod_{n=1}^N \mathbb{M}_n^{m(\beta_n)},
\quad \mbox{with } \quad
\mathbb{M}_n^{m(\beta_n)} = \sup_{\| f\|_{L^\infty(\Gamma_n)}=1} \| \mathscr{U}^{m(\beta_n)}_n f \|_{L^2_\rho(\Gamma_n)}.$$ In particular, for the Clenshaw-Curtis points, it is possible to bound $\mathbb{M}_n^{m(\beta_n)}$ as: $$\label{eq:CC_leb_est}
  \mathbb{M}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} \leq \mathbb{M}_{est}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} = \prod_{n=1}^N \mathbb{M}_{n,est}^{m(\beta_n)},  \quad \quad
  \mathbb{M}_{n,est}^{q} =
  \begin{cases}
    1 & \mbox{ for } q=1 \\[2mm]
    \displaystyle \frac{2}{\pi}\log(q-1)+1 & \mbox{ for } q \geq 2.
  \end{cases}$$ See [@nobile.eal:optimal-sparse-grids] and references therein.

Nested collocation points have been studied also for other probability measures than uniform probability measures. In the very relevant case of a normal distribution, one possible choice is the Genz-Keister points [@genz.keister:kpnquad; @back.nobile.eal:lognormal]; we mention also the recent work [@narayan:Leja] on generalized Leja points that can be used for arbitrary measures on unbounded domains.

Multi-Index Stochastic Collocation {#s:method}
==================================

It is easy to see that an accurate approximation of ${{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}}$ by a direct tensor technique as the one just introduced, ${{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}} \approx \mathscr{Q}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} [{F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}]$, might require a prohibitively large computational effort even for moderate values of $D$ and $N$ (what is referred to as the “curse of dimensionality”). In this work, following the setting that was presented in [@nobile.eal:optimal-sparse-grids; @abdullatif.etal:MultiIndexMC], we propose the Multi-Index Stochastic Collocation as an alternative. It can be seen as a generalization of the telescoping sum presented in the introduction, see equations and . Denoting $\mathscr{Q}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})} [{F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}] = {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}$, the building blocks of such a telescoping sum are the first-order difference operators for the deterministic and stochastic discretization parameters, denoted respectively by $\Delta_{i}^\textnormal{det}$ with $1\le i\le D$ and $\Delta_{j}^\textnormal{stoc}$ with $1\le j\le N$: $$\begin{aligned}
\label{eq:misc_delta_operations}
& \Delta_{i}^\textnormal{det}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}] =
\begin{cases}
  {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} - {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}-{{{\ensuremath{{\boldsymbol e}}}}}_i,{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}, & \text{if }\alpha_i > 1,\\
  {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} & \text{if } \alpha_i=1,
  \end{cases} \\& \Delta_{j}^\textnormal{stoc}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]    =
\begin{cases}
  {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} - {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}-{{{\ensuremath{{\boldsymbol e}}}}}_j}}, & \text{if }\beta_j > 1,\\
  {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} & \text{if } \beta_j=1.
  \end{cases} $$ We then define the first-order tensor difference operators, $$\begin{aligned}
& {{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{det}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]
= \bigotimes_{i=1}^D \Delta_i^\textnormal{det}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]
= \Delta_1^\textnormal{det} \left[ \, \Delta_2^\textnormal{det} \left[ \, \cdots
  \Delta_D^\textnormal{det} \left[ {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} \right] \, \right] \, \right] = \sum_{{{{\ensuremath{{\boldsymbol j}}}}}\in \{0,1\}^D} (-1)^{|{{{\ensuremath{{\boldsymbol j}}}}}|} {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}-{{{\ensuremath{{\boldsymbol j}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}},
\label{eq:delta_det_def} \\
& {{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{stoc}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]
= \bigotimes_{j=1}^N \Delta_j^\textnormal{stoc}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]
= \sum_{{{{\ensuremath{{\boldsymbol j}}}}}\in \{0,1\}^N} (-1)^{|{{{\ensuremath{{\boldsymbol j}}}}}|} {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}-{{{\ensuremath{{\boldsymbol j}}}}}}} \, .
\label{eq:delta_stoc_def}\end{aligned}$$ with the convention that $F_{{{\ensuremath{{\boldsymbol v}}}}, {{\ensuremath{{\boldsymbol w}}}}}=0$ whenever a component of ${{\ensuremath{{\boldsymbol v}}}}$ or ${{\ensuremath{{\boldsymbol w}}}}$ is zero. Observe that computing ${{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{det}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$ actually requires up to $2^D$ solver calls, and analogously applying ${{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{stoc}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$ requires interpolating ${F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$ on up to $2^N$ tensor grids; for instance, if $D=N=2$ and ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}> {\bm{1}}$, we have $$\begin{aligned}
  &{{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{det} [{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]
  = \Delta^\textnormal{det}_2 \left[ \, \Delta^\textnormal{det}_1 \left[ \, {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} \, \right] \, \right]
   = \Delta^\textnormal{det}_2 [ {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} - {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}-{{{\ensuremath{{\boldsymbol e}}}}}_1,{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]
   = {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} - {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}-{{{\ensuremath{{\boldsymbol e}}}}}_1,{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}
   - {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}-{{{\ensuremath{{\boldsymbol e}}}}}_2,{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} + {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}-{\bm{1}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}, \\[2pt]  &{{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{stoc} [{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]
  = {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} - {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}-{{{\ensuremath{{\boldsymbol e}}}}}_1}}
   - {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}-{{{\ensuremath{{\boldsymbol e}}}}}_2}} + {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}-{\bm{1}}}}. \nonumber\end{aligned}$$ Finally, letting ${{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}] =
{{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{stoc}[ {{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{det}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]]$, we define the Multi-Index Stochastic Collocation (MISC) estimator of ${{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}}$ as $$\label{eq:misc_estimator}
    \mathscr{M}_{\mathcal{I}}[{F}]
  = \sum_{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in \mathcal I} {{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]
  = \sum_{ [{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in \mathcal I}  c_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} {F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}},$$ where $\mathcal{I} \subset {\mathbb{N}}_+^{D+N}$ and $c_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \in {\mathbb{Z}}$. Observe that many of the coefficients in , $c_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$, may be zero: in particular, $c_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ is zero whenever $[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}]+{{\ensuremath{{\boldsymbol j}}}} \in \mathcal{I}$ $\forall {{\ensuremath{{\boldsymbol j}}}}\in\{0,1\}^{N+D}$. Similarly to the analogous sparse grid construction [@nobile.eal:optimal-sparse-grids; @wasi.wozniak:cost.bounds; @b.griebel:acta], we shall require that the multi-index set $\mathcal{I}$ be downward closed, i.e., $$\begin{aligned}
  \forall \,[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in \mathcal{I}, \quad
  \begin{cases}
    {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}- {{\ensuremath{{\boldsymbol e}}}}_i \in \mathcal{I} \mbox{ for } 1 \leq i \leq D \mbox{ and } \alpha_i > 1,\\
    {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}- {{\ensuremath{{\boldsymbol e}}}}_j \in \mathcal{I} \mbox{ for } 1 \leq j \leq N \mbox{ and } \beta_j > 1.
  \end{cases}\end{aligned}$$

In theory, a MISC approach could also be developed to approximate the entire solution $u({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})$ and not just the expectation of functionals, considering differences between consecutive interpolant operators, $\mathscr{U}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}$, on the stochastic domain rather than differences of the quadrature operators, $\mathscr{Q}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}$, as a building block for the ${{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{stoc}$ operators, as well as considering the discretized solution $u^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$ rather than just the quantity of interest, ${F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$, in the construction of the ${{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{det}$ operators.

A knapsack-like construction of the set $\mathcal{I}$
-----------------------------------------------------

The efficiency of the MISC method in equation will heavily depend on the specific choice of the index set, $\mathcal{I}$; in the following, we will first propose a general strategy to derive quasi-optimal sets and then prove in Section \[s:complexity\] a convergence result for such sets under some reasonable assumptions.

To derive an efficient set, $\mathcal{I}$, we recast the problem of its construction as an optimization problem, in the same spirit of [@abdullatif.etal:MultiIndexMC; @nobile.eal:optimal-sparse-grids; @back.nobile.eal:optimal; @b.griebel:acta; @griebel.knapek:optimized]. We begin by introducing the concepts of “work contribution”, $\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$, and “error contribution”, $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$, for each hierarchical surplus operator, ${{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$. The work contribution measures the computational cost (measured, e.g., as a function of the total number of degrees of freedom, or in terms of computational time) required to add ${{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$ to $\mathscr{M}_\mathcal{I}[{F}]$, i.e., to solve the associated deterministic problems and to compute the corresponding interpolants over the parameter space, cf. equations and ; the error contribution measures instead how much the error $|{{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}} - \mathscr{M}_\mathcal{I}[{F}]| $ would decrease once the operator ${{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$ has been added to $\mathscr{M}_\mathcal{I}[{F}]$. In formulas, we define $$\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}
= {{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathscr{M}_{\mathcal{I} \cup \{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}]\}}[{F}]\right]}} - {{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathscr{M}_{\mathcal{I}}[{F}]\right]}}
= {{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[{{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]\right]}},$$ so that $$\label{eq:work_decomp}
  {{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathscr{M}_\mathcal{I}[{F}]\right]}} = \sum_{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in \mathcal{I}} \Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}},$$ Observe that this work definition is sharp only if we think of building the MISC estimator with an incremental approach, i.e., we assume that adding the multi-index $({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})$ to the index set $\mathcal I$ would not reduce the work that has to be done to evaluate the MISC estimator on the index set. This implies that one cannot take advantage of the fact that some of the coefficients in , $c_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$, that are non-zero when considering the set $\mathcal{I}$ could become zero if the MISC estimator is instead built considering the set $\mathcal{I} \cup \{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}]\}$, hence it would be possible not to compute the corresponding approximations ${F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}$. This approach is discussed in Section \[ss:impl\].

Similarly, we define $$\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}
 = \Big| {\mathscr{M}_\mathcal{I \cup \{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}]\}}[{F}]} - {\mathscr{M}_\mathcal{I}[{F}]} \Big|
 = \left| {{{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]} \right|.$$ Thus, by construction, the error of the MISC estimator can be bounded as the sum of the error contributions not included in the estimator $\mathscr{M}_\mathcal{I}[{F}]$, $$\begin{aligned}
\label{eq:error-decomp}
{{\ensuremath{\mathrm{Error}}\mspace{-2mu}\left[\mathscr{M}_\mathcal{I}[{F}]\right]}} =
|{{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}} - \mathscr{M}_\mathcal{I}[{F}]|
& = \left| {\sum_{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \notin \mathcal I} {{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]} \right| \nonumber \\[2mm]
 & \leq \sum_{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \notin \mathcal I} \left| { {{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}] } \right|
  = \sum_{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \notin \mathcal I} \Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}.\end{aligned}$$ Consequently, a quasi-optimal set $\mathcal{I}$ can be computed by solving the following “binary knapsack problem” [@martello:knapsack]: $$\begin{aligned}
  \mbox{maximize} &\sum_{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in {\mathbb{N}}_+^{D+N}}\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} x_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \nonumber \\
  \mbox{such that }& \sum_{[{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in {\mathbb{N}}_+^{D+N}}\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} x_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \leq W_{max}, \label{eq:knapsack} \\
  & \qquad x_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \in \{0,1\}, \nonumber\end{aligned}$$ and setting $\mathcal{I}=\{ [{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in {\mathbb{N}}_+^{D+N}: x_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},
  {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} = 1\}$. Observe that such a set is only “quasi” optimal since the error decomposition is not an exact representation but rather an upper bound. The optimization problem above is well known to be computationally intractable. Still, an approximate greedy solution (which coincides with the exact solution under certain hypotheses that will be clearer in a moment) can be found if one instead allows the variables $x_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ to assume fractional values, i.e., it is possible to include fractions of multi-indices in $\mathcal{I}$. For this simplified problem, the resulting problem can be solved analytically by the so-called Dantzig algorithm [@martello:knapsack]:

1.  compute the “profit” of each hierarchical surplus, i.e., the quantity $$    P_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} = \frac{\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}{\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}};$$

2.  sort the hierarchical surpluses by decreasing profit;

3.  add the hierarchical surpluses to $\mathscr{M}_\mathcal{I}[{F}]$ according to such order until the constraint on the maximum work is fulfilled.

Note that by construction $x_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}=1$ for all the multi-indices included in the selection except for the last one, for which $x_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}<1$ might hold; in other words, the last multi-index is the only one that might not be taken entirely. However, if this is the case, we assume that we could slightly adjust the value $W_{max}$, so that all $x_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ have integer values (see also [@b.griebel:acta]); observe that this integer solution is also the solution of the original binary knapsack problem with the new value of $W_{max}$ in the work constraint. Thus, if the quantities $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ and $\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ were available, the quasi-optimal index set for the MISC estimator could be computed as $$\label{eq:opt_set}
 \mathcal{I} = \mathcal{I}(\epsilon) =\left\{ [{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in {\mathbb{N}}_+^{D+N}\::\:
   \frac{\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}{\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} \geq \epsilon \right\},$$ for a suitable $\epsilon>0$.

\[remark:webster\_and\_bieri\] The MISC setting could in principle include the Multilevel Stochastic Collocation method proposed in [@teckentrup.etal:MLSC] as a special case, by simply considering a discretization of the spatial domain on regular meshes, and letting the diameter of each element (the mesh-size) be the only discretization parameter, i.e., $D=1$.

However, the sparse grids to be used at each level are determined in [@teckentrup.etal:MLSC] by computing the minimal number of collocation points needed to balance the stochastic and spatial error. This is done by relying on sparse grid error estimates; yet, since in general it is not possible to generate a sparse grid with a predefined number of points, some rounding strategy to the sparse grid with the nearest cardinality must be devised, which may affect the optimality of the multilevel strategy. In the present work, we overcome this issue by relying instead on profit estimates to build a set of multi-indices that simultaneously prescribe the spatial discretization and the associated tensor grid in the stochastic variables. Furthermore, only standard isotropic Smolyak sparse grids are considered in the actual numerical experiments in [@teckentrup.etal:MLSC] (although in principle anisotropic sparse grids could be used as well, provided that good convergence estimates for such sparse grids are available), while our implementation naturally uses anisotropic stochastic collocation methods at each spatial level.

The MISC approach also includes as a special case the “Sparse Composite Collocation Method” developed in [@bieri:sparse.tensor.coll], by considering again only one deterministic discretization parameter, i.e., $D=1$, and then setting $$\label{eq:iso-MISC}
    \mathcal{I} = \left\{ [\alpha,{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in {\mathbb{N}}_+^{1+N} : \alpha + \sum_{n=1}^N \beta_n \leq w \right\},$$ with $w \in {\mathbb{N}}_+$. In other words, the approach in [@bieri:sparse.tensor.coll] is based neither on profit nor on error balancing.

Complexity analysis of the MISC method {#s:complexity}
======================================

In this section, we assume suitable models for the error and work contributions, $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ and $\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ (which are numerically verified in Section \[s:numerics\] for the problem in Example \[example:ell\]) and then state our main convergence theorem for the MISC method built using a particular index set, $\mathcal{I}^*$, which can be regarded as an approximation of the quasi-optimal set introduced in the previous section.

\[assump:growth\_of\_dof\] The discretization parameters, $h_i$, for the deterministic solver depend exponentially on the discretization level $\alpha_i$, and the number of collocation points over the parameter space grows exponentially with the level $\beta_i$: $$h_{i,\alpha_i}=h_0 2^{-\alpha_i} \qquad \text{and}\qquad  C_{m,low} 2^{\beta_i} \leq m(\beta_i) \leq C_{m,up} 2^{\beta_i}.$$

\[assump:dW\_dE\_factor\] The error and work contributions, $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ and $\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$, can be bounded as products of two terms, $$    \Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \leq \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}  \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}  \qquad \text{and}\qquad
    \Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \leq \Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}  \Delta
    W_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc},  $$ where $\Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}$ and $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}$ denote the cost and the error contributions due to the deterministic difference operator, $\Delta^\textnormal{det}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$, and similarly $\Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}$ and $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}$ denote the cost and the error contribution due to the stochastic difference operator, $\Delta^\textnormal{stoc}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$, cf. equations -.

\[assump:dW\_dE\_model\] The following bounds hold true for the factors appearing in Assumption \[assump:dW\_dE\_factor\]: $$\begin{aligned}
    & \Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det} \leq C_{\mathrm{work}}^{\textnormal{det}} \prod_{i=1}^D (h_{i,\alpha_i})^{-\widetilde{\gamma}_i} , \label{eq:deltaW_alpha_est}\\
    & \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det} \leq C_{\mathrm{error}}^{\textnormal{det}} \prod_{i=1}^D (h_{i,\alpha_i})^{\widetilde{{{r_{}}}_i} } , \label{eq:dE_alpha_est}\\
    & \Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc} \leq \widetilde{C}_{\mathrm{work}}^{\textnormal{stoc}} \prod_{n=1}^N m(\beta_n)
    \leq C_{\mathrm{work}}^{\textnormal{stoc}}  \prod_{n=1}^N 2^{\beta_n},\label{eq:deltaW_beta_est}  \\
    & \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc} \leq C_{\mathrm{error}}^{\textnormal{stoc}}\, e^{- \sum_{i=1}^N \widetilde{g}_i m(\beta_i)}, \label{eq:deltaE_beta_est}  \end{aligned}$$ for some rates $\widetilde{\gamma}_i, \widetilde{{{r_{}}}}_i, \widetilde{g}_i > 0$.

With these assumptions, we are now ready to state our main theorem. The proof is technical and we therefore place it in the appendix. The proof is based on summing the error contributions outside a particular index set, $\mathcal I^*$, and the work contributions inside the same index set. This can be seen as a weighted cardinality argument in finite dimensions. See also [@griebel:infdim; @griebel:tensor] for similar arguments for different choices of finite and infinite dimensional index sets.

[theorem]{}[misccomplexity]{} \[thm:misc\_complexity\] Under Assumptions \[assump:growth\_of\_dof\] to \[assump:dW\_dE\_model\], the bounds for the factors appearing in Assumption \[assump:dW\_dE\_factor\] can be equivalently rewritten as

$$\begin{aligned}
    \label{eq:deltaW_ass}
    & \Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \leq {\ensuremath{C_{\text{\textnormal{work}}}}}e^{\sum_{i=1}^D \gamma_i \alpha_i} e^{ \delta |{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}|}, \\
    \label{eq:deltaE_ass}
    & \Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \leq {\ensuremath{C_{\text{\textnormal{error}}}}}e^{ - \sum_{i=1}^D {{r_{i}}}\alpha_i} e^{- \sum_{j=1}^N g_j \exp(\delta {\beta_j})},
  \end{aligned}$$

with $\gamma_i = \widetilde{\gamma}_i\log 2$, ${{r_{i}}} =\widetilde{{{r_{}}}}_i\log 2$, $\delta = \log 2$ and $g_i = \widetilde{g}_i C_{m,low}$. Define the following set $$\label{eq:opt_set_detailed_logform}
    \mathcal{I}^*(L) =
    \left\{ [{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in {\mathbb{N}}_+^{D+N}\::\:
      \sum_{i=1}^D ({{r_{i}}}+\gamma_i) \alpha_i + \sum_{i=1}^N (\delta \beta_i +  g_i e^{\delta \beta_i} ) \leq L \right\} \mbox{ with } L \in {\mathbb{R}}_+.$$ Then there exists a constant $\mathscr C_{\text{W}}$ such that, for any $W_{\max}$ satisfying $$\label{eq:maxW_ass}
    W_{\max} \geq \mathscr C_{\text{W}} \exp\left(\chi\right),
      $$ and choosing $L$ as $$\label{eq:opt-L}
  L = L(W_{\max}) = \frac{1}{\chi} \left(
    \log\left(\frac{W_{\max}}{\mathscr{C}_{\text{W}}}\right) - (\mathfrak{z}-1)
    \log\left(\frac{1}{\chi}\log\left(\frac{W_{\max}}{\mathscr{C}_{\text{W}}}\right)\right)
        \right),$$ with ${{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}= \left( \frac{\gamma_1}{\gamma_1 + r_1},\ldots,\frac{\gamma_D}{\gamma_D + r_D} \right)$, $\chi = \max({{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}})$, $\zeta=\min_{i=1,\ldots,D} \frac{{{r_{i}}}}{\gamma_i}$ and $\mathfrak{z} = \#\{i=1,\ldots D \::\: \frac{{{r_{i}}}}{\gamma_i} =
  \zeta\}$, the MISC estimator $\mathscr M_{\mathcal I^*(L(W_{\max}))}$ satisfies

\[eq:misc\_complexity\] $$\begin{aligned}
      &{{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathscr M_{\mathcal I^*(L(W_{\max}))}\right]}} \leq W_{\max}, \label{eq:misc_complexity_work}
      \\[2mm]
      &\limsup_{W_{\max} \uparrow \infty} \frac
{{{\ensuremath{\mathrm{Error}}\mspace{-2mu}\left[\mathscr M_{\mathcal I^*(L(W_{\max}))}\right]}}}{W_{\max}^{-\zeta}        \left(\log\left(W_{\max}\right)\right)^{\left(\zeta+1\right)\left(\mathfrak{z}-1\right)}}
        = \mathscr C_{\text{E}} < \infty. \label{eq:misc_complexity_error}
    \end{aligned}$$

\[rem:this-is-the-set\] The set $\mathcal{I}^*$ proposed in Theorem \[thm:misc\_complexity\] can be obtained by assuming that the bounds in equations and are actually equalities and by using the definition of the quasi-optimal set : $$\begin{aligned}
    \mathcal{I}^*
    & = \left\{ [{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in {\mathbb{N}}_+^{D+N} : \frac{\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}{\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}} \geq \epsilon \right \} \\
    & = \left\{ [{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in {\mathbb{N}}_+^{D+N} :
      \frac{e^{ - \sum_{i=1}^D {{r_{i}}}\alpha_i} e^{- \sum_{j=1}^N g_j \exp(\delta {\beta_j})}}
      {e^{\sum_{i=1}^D \gamma_i \alpha_i} e^{ \delta |{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}|}} \geq \epsilon \right \} \\
    & = \left\{ [{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}] \in {\mathbb{N}}_+^{D+N}\::\:
      \sum_{i=1}^D ({{r_{i}}}+\gamma_i) \alpha_i + \sum_{i=1}^N (\delta \beta_i +  g_i e^{\delta \beta_i} ) \leq L \right\},
  \end{aligned}$$ where the last equality holds with $L = - \log \epsilon$.

\[rem:only-space-matters\] Refining along the spatial or the stochastic variables has different effects on the error of the MISC estimator. Indeed, due to the double exponential $e^{- \sum_{j=1}^N g_j \exp(\delta {\beta_j})}$ in , the stochastic contribution to the error will quickly fade to zero, which in turn implies that most of the work will be used to reduce the deterministic error. This is confirmed by the fact that the error convergence rate in Theorem \[thm:misc\_complexity\] only depends on $\gamma_i$ and ${{r_{i}}}$, i.e., the cost and error rates of the deterministic solver, respectively. This observation coincides with that in [@bieri:sparse.tensor.coll page 2299]: “since the stochastic error decreases exponentially, the convergence rate should tend towards the algebraic rate of the spatial discretization \[...\]; see Proposition 3.8”. Compared with the method proposed in [@bieri:sparse.tensor.coll], MISC takes greater advantage of this fact, since it is based on an optimization procedure, cf. equation ; this performance improvement is well documented by the comparison between the two methods shown in the next section. Figure \[fig:MISC\_set\] shows the multi-indices included in $\mathcal{I}$ according to for increasing values of $L$, for a problem with $N=D=1$, $\gamma_i=1,r=2$, and $g=1.5$: as expected, the shape of $\mathcal{I}$ becomes more and more curved as $L$ grows, due to this lack of balance between the stochastic and deterministic directions.

![Index sets $\mathcal{I}(L)$ for $D=N=1$, according to equation .[]{data-label="fig:MISC_set"}](opt_MISC_set.pdf){width="50.00000%"}

Theorem \[thm:misc\_complexity\] is only valid in case Assumptions \[assump:wellposed\]– \[assump:dW\_dE\_model\] are true. In the next section we motivate these assumption for a specific elliptic problem that we use for numerically testing the MISC method. However, we stress that deriving bounds on the error and work contributions is problem-dependent and the corresponding analysis must be carried out in each case. Moreover, under a different set of assumption the complexity theorem would have to be rewritten accordingly.

Example and numerical evidence {#s:numerics}
==============================

In this section, we test the effectiveness of the MISC approximation on some instances of the general elliptic equation in Example \[example:ell\]; more precisely, we consider a problem with one physical dimension ($d=1$) as well as a more challenging problem with three dimensions ($d=3$); in both cases, we set $\mathscr{B} = [0,1]^d$, ${\mathscr{F}}({{{\ensuremath{{\boldsymbol x}}}}})=1$. As for the random diffusion coefficient, we set $$\label{eq:log-uniform-a}
      a({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})= e^{\gamma_N({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})}, \quad
      \gamma_N({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}}) = \sum_{n=1}^N \lambda_n \psi_n({{{\ensuremath{{\boldsymbol x}}}}}) y_n,$$ where $y_n$ are uniform random variables over $\Gamma_n = [-1,1]$, $\lambda_n = \sqrt{3} \exp(-n)$ and take $\psi_n$ to be a tensorization of trigonometric functions. More precisely, we define the function $$\phi_n(x) =
\begin{cases}
\displaystyle \sin\left(\frac{n}{2} \pi x  \right) & \text{if $n$ is even} \\[8pt]
\displaystyle \cos\left(\frac{n-1}{2} \pi x\right) & \text{if $n$ is odd}
\end{cases}$$ and set $\psi_n(x) = \phi_n(x)$ if $d=1$. If $d=3$, we take $\psi_n({{\ensuremath{{\boldsymbol x}}}}) = \phi_{i(n)}(x_1)\phi_{j(n)}(x_2)\phi_{k(n)}(x_3)$ for some indices ${i(n), j(n), k(n)}$ detailed in Table \[tab:modes\]. Observe that the boundedness of the supports of the random variables $y_n$ guarantees the existence of the two bounding constants in equation , $a_{min}$ and $a_{max}$, that in turn assures the well posedness of the problem. Finally, the quantity of interest is defined as $$\label{eq:qoi_def}
  {F}({{{\ensuremath{{\boldsymbol y}}}}}) = \int_\mathscr{B} u({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}}) Q({{{\ensuremath{{\boldsymbol x}}}}}) d{{{\ensuremath{{\boldsymbol x}}}}},
  \quad Q({{{\ensuremath{{\boldsymbol x}}}}}) = \frac{1}{(\sigma\sqrt{2\pi})^d}\exp \left( -\frac{\|{{{\ensuremath{{\boldsymbol x}}}}}-{{{\ensuremath{{\boldsymbol x}}}}}_0\|_2^2}{2\sigma^2} \right)$$ with $\sigma=0.16$ and locations ${{{\ensuremath{{\boldsymbol x}}}}}_0=0.3$ for $d=1$ and ${{{\ensuremath{{\boldsymbol x}}}}}_0=[0.3, 0.2, 0.6]$ for $d=3$. We also make the choice $h_0=1/3$ in Assumption \[assump:growth\_of\_dof\]. With these values and using the coarsest discretization, $h_0=1/3$, in all dimensions, the coefficient of variation of the quantity of interest can be approximated to be between $90\%$ and $100\%$ depending on the number of dimensions, $d$, and the number of random variables, $N$, that we consider below.

  $n$        1   2   3   4   5   6   7   8   9   10
  -------- --- --- --- --- --- --- --- --- --- ----
  $i(n)$     1   2   1   1   3   2   2   1   1    1
  $j(n)$     1   1   2   1   1   2   1   3   2    1
  $k(n)$     1   1   1   2   1   1   2   1   2    3

  : Included functions for $d=3$ in . Here $\psi_n({{\ensuremath{{\boldsymbol x}}}}) = \phi_{i(n)}(x_1)\phi_{j(n)}(x_2)\phi_{k(n)}(x_3)$.[]{data-label="tab:modes"}

Verifying bounds on work and error contributions
------------------------------------------------

In this subsection we discuss the validity of Assumptions \[assump:growth\_of\_dof\] to \[assump:dW\_dE\_model\], upon which the MISC convergence theorem is based. To this end, we analyze separately the properties of the deterministic solver and of the collocation method applied to the problem just introduced.

#### Deterministic solver

The deterministic solver considered in this work consists of a tensorized finite difference solver, with the grid size along each direction, $x_1,\ldots,x_d$, defined by $h_{i,\alpha_i} = h_0 2^{-\alpha_i}$ and no other discretization parameters are considered: therefore, $D=d$, Assumption \[assump:growth\_of\_dof\] is satisfied, and, due to the Dirichlet boundary conditions prescribed for $u$, the overall number of degrees of freedom of the corresponding finite difference solution is $\prod_{i=1}^D \left( \frac{1}{h_{i,\alpha_i}}-1 \right)
\leq \prod_{i=1}^D \left( \frac{1}{h_{i,\alpha_i}} \right)$. The associated linear system is solved with the GMRES method. We have numerically fitted the parameters, $\vartheta$ and $C_{\mathrm{GMRES}}$, in the model: $${{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[{F}_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\right]}} \leq C_{\mathrm{GMRES}}  \prod_{i=1}^D
  (h_{i,\alpha_i})^{-\vartheta},$$ for each individual tensor grid solve and found that $\vartheta=1$ gives a good fit in our numerical experiments. From this we can recover the rates $\{\tilde{\gamma}_i\}_{i=1}^D$ and the constant $C_{\mathrm{work}}^{\textnormal{det}}$ in with the following argument: since computing ${{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{det}[{F}^{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}]$ requires up to $2^D$ solver calls, each on a different grid (cf. equation ), we have $$\begin{aligned}
\Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}
& = {{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[{{\ensuremath{{\boldsymbol \Delta}}}}^\textnormal{det}[ {F}_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}]\right]}}
= \sum_{{{{\ensuremath{{\boldsymbol j}}}}}\in \{0,1\}^D} {{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[{F}_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}-{{{\ensuremath{{\boldsymbol j}}}}}}\right]}} \\
& \leq C_{\mathrm{GMRES}}\sum_{{{{\ensuremath{{\boldsymbol j}}}}}\in \{0,1\}^D} \prod_{i=1}^D \left( h_0 2^{-(\alpha_i-j_i)}\right)^{-\vartheta} \\
& = C_{\mathrm{GMRES}} \left( \prod_{i=1}^D \left( h_0 2^{-\alpha_i}\right)^{-\vartheta} \right) \sum_{{{{\ensuremath{{\boldsymbol j}}}}}\in \{0,1\}^D} \prod_{i=1}^D 2^{-j_i \vartheta} \\
& = C_{\mathrm{GMRES}} (1+ 2^{-\vartheta})^D \prod_{i=1}^D  (h_{i,\alpha_i})^{-\vartheta},\end{aligned}$$ i.e., bound is verified with $\widetilde{\gamma}_i = \vartheta, \forall i=1,\ldots,D$ and $C_{\mathrm{work}}^{\textnormal{det}} = C_{\mathrm{GMRES}} (1+
2^{-\vartheta})^D$. Hence, the sum of costs of the solver calls is proportional to the cost of the call on the finest grid.

Concerning the error contribution $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}$, we observe numerically that bound holds true in practice with $\widetilde{{{r_{}}}}_i=2$, $i=1,\ldots,D$, due to the fact that $a \in C^{\infty}(\mathscr{B})$ for $\rho$-almost every ${{{\ensuremath{{\boldsymbol y}}}}}\in \Gamma$, ${\mathscr{F}}\in C^{\infty}(\mathscr{B})$ and the function $Q$ appearing in the quantity of interest is also infinitely differentiable, confined in a small region inside the domain and zero up to machine precision on the boundary. In more detail, assuming for a moment that Assumption \[assump:dW\_dE\_factor\] is valid (we will numerically verify it later in this section), in Figure \[fig:check\_rates\_det\] we show the value of $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} = \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det} \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}$ for fixed ${{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}={\bm{1}}$ and variable ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}= j \bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}} + {\bm{1}}, j=1,2,\ldots$, as well as the corresponding value of the bound for $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}$. The line obtained by choosing $\bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}} = [1,\, 0,\, 0]$ confirms that the size of $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}$ indeed decreases exponentially fast with respect to $\alpha_1$, and by fitting the computed values of $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}$ we obtain that the convergence rate is $\widetilde{{{r_{}}}}_j=2$, as previously mentioned; analogous conclusions can be obtained for $\alpha_2$ and $\alpha_3$ by setting $\bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}} = [0,\, 1,\, 0]$ (shown in Figure \[fig:check\_rates\_det\]) and $\bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}} = [0,\, 0,\, 1]$ (not shown). Most importantly, confirmation of the product structure of $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}$ can be obtained by observing, e.g., the decay of $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}$ for $\bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}} = [1,\, 1,\, 0]$ and $\bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}} = [1,\, 1,\,1]$.

#### Stochastic discretization

The interpolation over the parameter space is based on the tensorized Lagrangian interpolation technique with Clenshaw-Curtis points explained in Section \[s:interpolation\], cf. eqs. and . In particular, due to the nestedness of the Clenshaw-Curtis points, adding the operator $\Delta^\textnormal{stoc}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$ to the MISC estimator will require $\Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc} = \prod_{j=1}^N \big( m(\beta_j) - m(\beta_j-1) \big)$ new collocation points, which in view of equation can be bounded as $$m(\beta_j) - m(\beta_j-1) =
  \begin{cases}
    1 \mbox{ if }  \beta_j=1 \\
    2 \mbox{ if } \beta_j=2 \\
    2^{\beta_j-2}, \mbox{ if } \beta_j>2,
  \end{cases}
  \leq 2^{\beta_j{-1}}, \quad \forall j=1,2,\ldots, $$ provided that the set $\mathcal{I}$ is downward closed: Assumption \[assump:growth\_of\_dof\] and bound in Assumption \[assump:dW\_dE\_model\] are thus verified. Observe that the nestedness of the Clenshaw-Curtis knots is a key property here: indeed, if the nodes are not nested $\Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}$ is not uniquely defined, i.e., it depends on the set $\mathcal{I}$ to which $\Delta^\textnormal{stoc}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$ is being added, see, e.g., [@nobile.eal:optimal-sparse-grids Example 1 in Section 3].

Finally, to discuss the validity of bound for $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}$, we rely on the theory developed in our previous works [@back.nobile.eal:optimal; @nobile.eal:optimal-sparse-grids]. We begin by introducing the Chebyshev polynomials of the first kind ${\Psi}_{q}(y)$ for $y \in [-1,1]$, which are defined by the relation $${\Psi}_{q}(\cos(\theta))
= \cos(q \theta), \quad 0 \leq \theta \leq \pi, \quad q \in {\mathbb{N}}.$$ Then, for any multi-index ${{{\ensuremath{{\boldsymbol q}}}}}\in {\mathbb{N}}^N$, we consider the $N$-variate Chebyshev polynomials ${\bm{\Psi}}_{{{\ensuremath{{\boldsymbol q}}}}}({{{\ensuremath{{\boldsymbol y}}}}}) = \prod_{n=1}^N {\Psi}_{q_n}(y_n)$ and introduce the spectral expansion of $f:[-1,1]^N\rightarrow {\mathbb{R}}$ over $\{{\bm{\Psi}}_{{{\ensuremath{{\boldsymbol q}}}}}\}_{{{{\ensuremath{{\boldsymbol q}}}}}\in {\mathbb{N}}^N}$, $$f({{{\ensuremath{{\boldsymbol y}}}}}) = \sum_{{{{\ensuremath{{\boldsymbol q}}}}}\in {\mathbb{N}}^N} f_{{{\ensuremath{{\boldsymbol q}}}}}{\bm{\Psi}}_{{{\ensuremath{{\boldsymbol q}}}}}({{{\ensuremath{{\boldsymbol y}}}}}),
\quad
f_{{{\ensuremath{{\boldsymbol q}}}}}= \int_\Gamma f({{{\ensuremath{{\boldsymbol y}}}}}) {\bm{\Psi}}_{{{\ensuremath{{\boldsymbol q}}}}}({{{\ensuremath{{\boldsymbol y}}}}}) \prod_{n=1}^N
\frac{1}{\sqrt{1-y_n^2}} d{{{\ensuremath{{\boldsymbol y}}}}},
$$ Next, given any ${\xi}_{1},{\xi}_{2},\ldots,{\xi}_{N}>1$ we introduce the Bernstein polyellipse ${\mathcal{E}}_{{\xi}_1,\ldots,{\xi}_N} =\prod_{n=1}^N {\mathcal{E}}_{n,{\xi}_n}$, where ${\mathcal{E}}_{n,{\xi}_n}$ denotes the ellipses in the complex plane defined as $$\begin{aligned}
  {\mathcal{E}}_{n,{\xi}_n} = \bigg\{ z_n \in \mathbb{C} : \,\,
  { \mathfrak{Re}\left( z \right)} \leq \frac{{\xi}_{n} + {\xi}_{n}^{-1}}{2}\cos{\phi}, \,\,\,
  {\mathfrak{Im}\left( z \right)} \leq \frac{{\xi}_{n} - {\xi}_{n}^{-1}}{2}\sin{\phi},
  \,\, \phi \in [0, 2\pi) \bigg\},\end{aligned}$$ and recall the following lemma (see [@nobile.eal:optimal-sparse-grids Lemma 2] for a proof).

\[lemma:inc\_cheb\_conv\] Let $f: [-1,1]^N \rightarrow {\mathbb{R}}$, and assume that there exist ${\xi}_{1},{\xi}_{2},\ldots,{\xi}_{N}>1$ such that $f$ admits a complex continuation $f^*:\mathbb{C}^N\rightarrow {\mathbb{R}}$ holomorphic in the Bernstein polyellipse ${\mathcal{E}}_{{\xi}_1,\ldots,{\xi}_N}$ with $\sup_{{{{\ensuremath{{\boldsymbol z}}}}}\in {\mathcal{E}}_{{\xi}_1,\ldots,{\xi}_N}} |f^*({{{\ensuremath{{\boldsymbol z}}}}})| \leq B$ and $B = B({\xi}_{1},{\xi}_{2},\ldots, {\xi}_{n}) < \infty$. Then $f$ admits a Chebyshev expansion that converges in $C^0([-1,1]^N)$, and whose coefficients $f_{{{\ensuremath{{\boldsymbol q}}}}}$ are such that $$\label{eq:polyell_cheb_conv}
    | f_{{{\ensuremath{{\boldsymbol q}}}}}| \leq C_\textnormal{Cheb}({{{\ensuremath{{\boldsymbol q}}}}}) \prod_{n=1}^N e^{-g^*_n q_n},  \quad g^*_n = \log {\xi}_{n}$$ with $C_\textnormal{Cheb}({{{\ensuremath{{\boldsymbol q}}}}}) = 2^{\|{{{\ensuremath{{\boldsymbol q}}}}}\|_0} B$, where $\|{{{\ensuremath{{\boldsymbol q}}}}}\|_0$ denotes the number of non-zero elements of ${{{\ensuremath{{\boldsymbol q}}}}}$.

The following lemma then shows that the region of analyticity of ${F}({{{\ensuremath{{\boldsymbol y}}}}})$ indeed contains a Bernstein ellipse, so that a decay of exponential type can be expected for its Chebyshev coefficients.

\[lemma:qoi\_polyellipse\] The quantity of interest ${F}({{{\ensuremath{{\boldsymbol y}}}}})=\Theta[u(\cdot,{{{\ensuremath{{\boldsymbol y}}}}})]$ is analytic in a Bernstein polyellipse with parameters ${\xi}_n = \tau_n + \sqrt{\tau_n^2+1}$, for any $\tau_n < \frac{\pi}{2 N \lambda_n}$.

Equation can be extended in the complex domain by replacing ${{{\ensuremath{{\boldsymbol y}}}}}$ with ${{{\ensuremath{{\boldsymbol z}}}}}\in \mathbb{C}^N$, and is analytic in the set $\Sigma = \left\{ {{{\ensuremath{{\boldsymbol z}}}}}\in \mathbb{C}^N : \Re e{[a({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol z}}}}})]} >0 \right\}$, see, e.g., [@babuska.nobile.eal:stochastic2]. By writing $z_n = b_n + i c_n$, we have $$\begin{aligned}
a({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol z}}}}})
& = \exp \left( \sum_n z_n \lambda_n \psi_n({{{\ensuremath{{\boldsymbol x}}}}}) \right)
  = \exp \left( \sum_n b_n \lambda_n \psi_n({{{\ensuremath{{\boldsymbol x}}}}}) \right) \exp \left( \sum_n i c_n \lambda_n \psi_n({{{\ensuremath{{\boldsymbol x}}}}}) \right) \\
& = \exp \left( \sum_n b_n \lambda_n \psi_n({{{\ensuremath{{\boldsymbol x}}}}}) \right)
    \left[ \cos \left(\sum_n c_n \lambda_n \psi_n({{{\ensuremath{{\boldsymbol x}}}}}) \right) + i \sin \left( \sum_n c_n \lambda_n \psi_n({{{\ensuremath{{\boldsymbol x}}}}}) \right) \right]\end{aligned}$$ so that the region $\Sigma$ can be rewritten as $$\Sigma =
\left\{ {{{\ensuremath{{\boldsymbol z}}}}}= {{{\ensuremath{{\boldsymbol b}}}}}+ i {{{\ensuremath{{\boldsymbol c}}}}}\in \mathbb{C}^N : \cos \left( \sum_n c_n \lambda_n \psi_n({{{\ensuremath{{\boldsymbol x}}}}}) \right) >0 , \forall {{{\ensuremath{{\boldsymbol x}}}}}\in \mathscr{B} \right\}.$$ Such a region includes the smaller region $$\Sigma_2 = \left\{ {{{\ensuremath{{\boldsymbol z}}}}}= {{{\ensuremath{{\boldsymbol b}}}}}+ i {{{\ensuremath{{\boldsymbol c}}}}}\in \mathbb{C}^N :
\left\| \sum_n c_n \lambda_n \psi_n \right\|_{L^\infty(\mathscr{B})} < \frac{\pi}{2} \right\},$$ which in turn includes $$\Sigma_3
= \left\{ {{{\ensuremath{{\boldsymbol z}}}}}= {{{\ensuremath{{\boldsymbol b}}}}}+ i {{{\ensuremath{{\boldsymbol c}}}}}\in \mathbb{C}^N :
\sum_n  \lambda_n |c_n|  < \frac{\pi}{2} \right\},$$ where the last equality is due to the fact that, by construction, $\left\|\psi_n\right\|_{L^\infty(\mathscr{B})}=1$, cf. equation . Next we let $\tau_n = \frac{\pi}{2 N \lambda_n}$ and define the following subregion of $\Sigma_3$: $$\Sigma_4 = \bigg\{ {{{\ensuremath{{\boldsymbol z}}}}}= {{{\ensuremath{{\boldsymbol b}}}}}+ i {{{\ensuremath{{\boldsymbol c}}}}}\in \mathbb{C}^N : |c_n|  < \tau_n \bigg\} \subset \Sigma_3.$$ $\Sigma_4$ is actually a polystrip in the complex plain that it in turn contains the Bernstein ellipse with parameters ${\xi}_n$ such that $$\frac{{\xi}_n - {\xi}_n^{-1}}{2} = \tau_n
\Rightarrow {\xi}_n^2 - 1 - 2\tau_n{\xi}_n = 0
\Rightarrow {\xi}_n = \tau_n + \sqrt{\tau_n^2+1}$$ in which $u({{{\ensuremath{{\boldsymbol x}}}}},{{{\ensuremath{{\boldsymbol y}}}}})$ is analytic. Finally, the quantity of interest, ${F}=\Theta[u]$, is also analytic in the same Bernstein polyellipse due to the linearity of the operator $\Theta$.

Incidentally, we remark that the choice of $\tau_n$ considered in Lemma \[lemma:qoi\_polyellipse\] degenerates for $N \to \infty$. In this case, if we know that $\sum_{n=0}^\infty (\lambda_n \|\psi_n\|_{L^\infty(\mathscr{B})})^p  < \infty$ for some $p<1$, then we could set $\tau_n = \frac{\pi}{2} (\lambda_n \|\psi_n\|_{L^\infty(\mathscr{B})})^{p-1}$, which does not depend on $N$.

$$\label{eq:deltaE_beta_gen}
    \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc} \leq C_E e^{-\sum_{n=1}^N g^*_n m(\beta_n-1)} \mathbb{M}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}$$

holds, where $C_E= 4^N B \prod_{n=1}^N \frac{1}{1-e^{-g^*_n}}$, $B$ as in Lemma \[lemma:inc\_cheb\_conv\], $g^*_n = \log {\xi}_n $ with ${\xi}_n$ as in Lemma \[lemma:qoi\_polyellipse\], and $\mathbb{M}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}$ has been defined in equation .

Combining Lemmas \[lemma:inc\_cheb\_conv\] and \[lemma:qoi\_polyellipse\], we obtain that the Chebyshev coefficients of ${F}$ can be bounded as $$| {F}_{{{\ensuremath{{\boldsymbol q}}}}}| \leq C_\textnormal{Cheb}({{{\ensuremath{{\boldsymbol q}}}}}) \prod_{n=1}^N e^{-g^*_n q_n},$$ with $g^*_n = \log {\xi}_n = \log (\tau_n + \sqrt{\tau_n^2+1})$ and $\tau_n$ as in Lemma \[lemma:qoi\_polyellipse\]. Then, the result can be obtained following the same argument of [@nobile.eal:optimal-sparse-grids Lemma 5].

To conclude, we first observe that $\mathbb{M}^{m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})}$ grows logarithmically with respect to $m({{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}})$, see eq. , so it is asymptotically negligible in the estimate above, i.e. we can write $$\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc} \leq C_{E_2}(\epsilon) \prod_{n=1}^N e^{-g^*_n (1-\epsilon_E) m(\beta_n-1)}$$ for an arbitrary $\epsilon_E > 0$ and with $C_{E_2}(\epsilon_E)>C_E$, and furthermore that the definition of $m(i)$ in implies that $m(i-1) \geq \frac{m(i)-1}{2}$. We can finally write $$\begin{aligned}
  \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}
  & \leq C_{E_2}(\epsilon) \prod_{n=1}^N e^{-g^*_n (1-\epsilon_E) \frac{m(\beta_n)-1}{2}}    = C_{\mathrm{error}}^{\textnormal{stoc}} \prod_{n=1}^N e^{-\widetilde{g}_n m(\beta_n)}, $$ with $C_{\mathrm{error}}^{\textnormal{stoc}} = C(\epsilon) \prod_{n=1}^N e^{\frac{g^*_n}{2} (1-\epsilon)}$ and $\widetilde{g}_n = \frac{g^*_n}{2} (1-\epsilon_E)$. The latter bound actually shows that bound in Assumption \[assump:dW\_dE\_model\] is valid for the test we are considering. Finally, we point out that in practice we work with the expression , whose rates $g_n$ are actually better estimated numerically, using the same procedure used to obtain the deterministic rates $\tilde {{{r_{}}}}_j=2$: we choose a sufficiently fine spatial resolution level ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$, consider a variable ${{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}= j\bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} + {\bm{1}}$ and fit the (simplified) model $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc} \leq C \prod_{n=1}^N e^{-g_n 2^{\beta_n}}$. The values obtained are reported in Table \[tab:g\_values\], and they are found to be equal for the case $d=1$ and $d=3$ (see also [@back.nobile.eal:comparison; @back.nobile.eal:optimal; @nobile.tempone.eal:aniso]). To make sure that the estimated value of $g_n$ does not depend on the spatial discretization, one could repeat the procedure for a few different values of ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$ and verify that the estimate is robust with respect to the spatial discretization: we note, however, that a rough estimate of $g_n$ will also be sufficient, since the convergence of the method is in practice dictated by the deterministic solver, as we have already discussed in Remark \[rem:only-space-matters\]. Figure \[fig:check\_rates\_stoc\] then shows the validity of the bound $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc} \leq C \prod_{n=1}^N e^{-g_n 2^{\beta_n}}$ comparing for fixed ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}={\bm{1}}$ and ${{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}= j\bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} + {\bm{1}}$ the value of $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}^\textnormal{det} \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}$ and the corresponding estimate.

      $g_1$    $g_2$    $g_3$    $g_4$    $g_5$    $g_6$    $g_7$    $g_8$    $g_9$    $g_{10}$
  -- -------- -------- -------- -------- -------- -------- -------- -------- -------- ----------
      2.4855   2.8174   4.5044   4.1938   4.7459   6.8444   7.1513   7.8622   8.6584    9.4545

  : Values of rates $g$ for the test cases considered.[]{data-label="tab:g_values"}

#### Stochastic-deterministic product structure

We conclude this section by verifying Assumption \[assump:dW\_dE\_factor\], i.e., the fact that the error contribution can be factorized as $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} = \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det} \Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}$ and that an analogous decomposition holds for $\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$. While the latter is trivial, to verify the former we employ the same strategy used to verify the models for $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det}$ and $\Delta E_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}$, this time letting both ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}$ and ${{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}$ change for every point, i.e., ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}= j \bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}} + {\bm{1}}$ and ${{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}= j {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}_0 + {\bm{1}}$. Figure \[fig:check\_mixed\_rates\] shows the comparison between the computed value of $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ and their estimated counterpart and confirms the validity of the product structure assumption.

![Comparison of $|\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}|$ for ${{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}=j\bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}+{\bm{1}}$ and ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}= j\bar{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}+{\bm{1}}$ for the test case with $D=3$ and $N=5$. The [*dashed*]{} lines are based on the model in with ${\tilde{{{r_{}}}}_i=2}$ for all $i=1,2,3$ and $g_j$ as in Table \[tab:g\_values\] for $j=1,\ldots, 5$. The [*solid*]{} lines are based on computed values.[]{data-label="fig:check_mixed_rates"}](COOK_d3_N5_rates.pdf){width="50.00000%"}

Test setup
----------

In our numerical tests, we compare MISC with the methods listed below. For each of them we show (for both test cases considered) plots of the convergence of the error in the computation of ${{\ensuremath{\mathbb{E}}\mspace{-2mu}\left[{F}\right]}}$ with respect to the computational work, taking as a reference value the result obtained using a well-resolved MISC solution. To avoid discrepancies in running time due to implementation details, the computational work is estimated in terms of the total number of degrees of freedom, i.e., using and . The names used here for the methods are also used in the legends of the figures showing the convergence plots.

“a-priori” MISC

:   refers to the MISC method with index set $\mathcal I$ defined by , where $\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ and $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ are taken to equal their upper bounds in and , respectively. The resulting set is explicitly written in . The convergence rate of this set is predicted by Theorem \[thm:misc\_complexity\], cf. Remark \[rem:this-is-the-set\]. Note that we do not need to determine the value of the constants ${\ensuremath{C_{\text{\textnormal{work}}}}}$ and ${\ensuremath{C_{\text{\textnormal{error}}}}}$ since they can be absorbed in the parameter $\epsilon$ in .

“a-posteriori” MISC

:   refers to the MISC method with index set $\mathcal I$ defined by , where $\Delta W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ is taken to equal its upper bound in , and $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ is instead computed explicitly as $\left| \Delta[{F_{\alpha,\beta}}] \right|$. Notice that this method is not practical since the cost of constructing set $\mathcal I$ would dominate the cost of the MISC estimator by far. However, this method would produce the best possible convergence and serve as a benchmark for both “a-priori” MISC and the bound .

MLSC

:   (only in the case $d>1$) refers to the Multilevel Stochastic Collocation obtained by setting $\alpha_1=\ldots=\alpha_D$ (i.e. considering the mesh-size as the only discretization parameter), as already mentioned in Remark \[remark:webster\_and\_bieri\]; we recall this is not exactly the MLSC method that was implemented in [@teckentrup.etal:MLSC], see again Remark \[remark:webster\_and\_bieri\]. Just as with MISC, we consider both the “a-priori” and “a-posteriori” version of MLSC, where $\Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}$ is taken to be equal to its upper in the former case and assessed by direct numerical evaluation in the latter case.

SCC

:   refers to the “Sparse Composite Collocation method” in Remark \[remark:webster\_and\_bieri\], see equation .

MIMC

:   refers to the Multi-Index Monte Carlo method as detailed in [@abdullatif.etal:MultiIndexMC], for which the complexity ${ {\ensuremath{ \mathcal O\left( W_{\max}^{-0.5} \right)}} }$ can be estimated for the test case at hand and as long as $d<4$.

SGSC

:   refers to the quasi-optimal Sparse Grids Stochastic Collocation (SGSC) with fixed spatial discretization as proposed in [@back.nobile.eal:optimal; @nobile.eal:optimal-sparse-grids]. To determine the needed spatial discretization for a given work and for a fair comparison against MISC, we actually compute the convergence curves of SGSC for all relevant levels of spatial discretizations and then show in the plots only the lower envelope of the corresponding convergence curves, ignoring the possible spurious reductions of error that might happen due to non-asymptotic, unpredictable cancellations, cf. Figure \[fig:envelope\]. In this way, we ensure that the error shown for such “single-level methods” has been obtained with the smallest computational error possible. Again, this is not a computationally practical method but is taken as a reference for what a sparse grids Stochastic Collocation method with optimal balancing of the space and stochastic discretization errors could achieve.

![Envelope of SGSC convergence curves for the test case with $d=3$ and $N=10$.[]{data-label="fig:envelope"}](COOK_d3_N10.pdf){width="70.00000%"}

Implementation details {#ss:impl}
----------------------

To implement MISC, we need two components:

1.  Given a profit level parameter, $\epsilon$, we build the quasi-optimal set $\mathcal I$ based on , and . One method to achieve this is to exploit the fact that this set is downward closed and use the following recursive algorithm.

        FUNCTION BuildSet(epsilon, multiIndex)
            FOR i = 1 to (D+N)
                IF Profit(multiIndex + e_i) > epsilon
                THEN
                    ADD multiIndex+e_i to FinalSet
                    CALL BuildSet(epsilon, multiIndex+e_i)
                END IF
            END FOR
        END FUNCTION
          

2.  Given the set, $\mathcal I(L)$, we evaluate . Here we have two choices:

    -   Evaluate the individual terms ${{\ensuremath{{\boldsymbol \Delta}}}}[{F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}]$ for every ${{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}\in \mathcal I$. To do so, we use the operator defined in  along each stochastic and spatial direction. By storing the values of these terms, we can evaluate the MISC with different index sets (contained in $\mathcal I(L)$), which might be required to test the convergence of the MISC method. Moreover, this implementation is suitable for adaptive methods that expand the index set based on some criteria and reevaluate the MISC estimator. On the other hand, this implementation has a computational overhead since most computed values of ${F_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}$ will actually not contribute to the final value of the estimator. However, this computational overhead is only a fraction of the minimum time required to evaluate the estimator.

    -   Use the combination form of  and only compute the terms that have $ c_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \neq 0$. This would remove the overhead of computing terms that make zero contribution to the estimator. This implementation is more efficient but less flexible as we cannot evaluate the estimator on sets contained in $\mathcal I(L)$ or build the set adaptively.

Test with $D=1$
---------------

Here we consider three different numbers of stochastic variables, namely $N=1,5,10$. Results are shown in Figure \[fig:res\_D=1\]. As expected, a-posteriori MISC shows the best convergence, with a-priori MISC being slightly worse and the single level methods following. Finally, we verify the accuracy of the estimated asymptotic convergence rate provided by Theorem \[thm:misc\_complexity\]: in this case, $\zeta = \frac{{{r_{}}}_1}{\gamma_1} = \frac{\widetilde{{{r_{}}}}_1 \log 2}{\widetilde{\gamma}_1 \log 2} = 2$ and $\mathfrak{z}=1$ holds. Hence, the predicted convergence rate is $W^{-\zeta}(\log W)^{(\zeta +1)(\mathfrak{z}-1)} = W^{-2}$, which appears to be in good agreement with the experimental convergence rate.

![Results for test $D=1$, case $N=1$ (top), $N=5$ (center) and $N=10$ (bottom).[]{data-label="fig:res_D=1"}](COOK_d1_N1.pdf "fig:"){width="80.00000%"}\
![Results for test $D=1$, case $N=1$ (top), $N=5$ (center) and $N=10$ (bottom).[]{data-label="fig:res_D=1"}](COOK_d1_N5.pdf "fig:"){width="80.00000%"}\
![Results for test $D=1$, case $N=1$ (top), $N=5$ (center) and $N=10$ (bottom).[]{data-label="fig:res_D=1"}](COOK_d1_N10.pdf "fig:"){width="80.00000%"}\

Test with $D=3$
---------------

In this case, we obtain the convergence curves shown in Figure \[fig:res\_D=3\], where the Multilevel Stochastic Collocation method has also been included. The hierarchy between the methods is in agreement with the case $d=1$, with the Multilevel Stochastic Collocation being comparable or slightly better than single level methods, but worse than the MISC approaches as expected.

Concerning the accuracy of the theoretical estimate: since for this test $\widetilde{{{r_{}}}}_1 = \widetilde{{{r_{}}}}_2 = \widetilde{{{r_{}}}}_2 = 2$ and $\widetilde{\gamma}_1 = \widetilde{\gamma}_2 = \widetilde{\gamma}_3=1$, $\zeta = 2$ still holds, while this time $\mathfrak{z}=3$; hence, the predicted convergence rate is $W^{-\zeta}(\log W)^{(\zeta +1)(\mathfrak{z}-1)} = W^{-2} (\log W)^6$. The plots suggest that the theoretical estimates might be slightly too optimistic when $N$ increases but it is important to recall that Theorem \[thm:misc\_complexity\] gives only an asymptotic result, and the plot could be negatively influenced by pre-asymptotic effects. Observe also that in this case there are a few data points where a-posteriori MISC is not better than a-priori MISC; this observation can be ascribed to the fact that a-posteriori MISC is optimal only with respect to the upper bound in . In other words, a-posteriori MISC selects the contributions according to the absolute value of the contributions but then the MISC estimator is computed by summing signed contributions. Hence, cancellations between contributions with similar sizes and opposite signs will occur.

Finally, we remark that, in our calculations, MLSC and SGSC were not able to achieve very small errors, unlike MISC. This is due to a limitation in the linear solver we are using that allows systems with only up to $2^{17}$ degrees of freedom (around 1GB of memory) to be solved. These “single-level” methods hit that limit sooner than MISC since they entail solving a very large system that comes from isotropically discretizing all three spatial dimensions.

![Results for test $D=3$, case $N=1$ (top), $N=5$ (center) and $N=10$ (bottom). []{data-label="fig:res_D=3"}](COOK_d3_N1.pdf "fig:"){width="80.00000%"}\
![Results for test $D=3$, case $N=1$ (top), $N=5$ (center) and $N=10$ (bottom). []{data-label="fig:res_D=3"}](COOK_d3_N5.pdf "fig:"){width="80.00000%"}\
![Results for test $D=3$, case $N=1$ (top), $N=5$ (center) and $N=10$ (bottom). []{data-label="fig:res_D=3"}](COOK_d3_N10.pdf "fig:"){width="80.00000%"}

Conclusions {#s:conclusions}
===========

In this work, we have proposed [MISC, a combination technique method]{} to solve UQ problems, optimizing both the deterministic and stochastic resolution levels simultaneously to minimize the computational cost. [A distinctive feature of MISC is that its construction is based on the notion of profit of the mixed differences composing it, rather than balancing the total error contributions arising from the deterministic and stochastic components.]{} We have detailed a complexity analysis and derived a convergence theorem showing that in certain cases the convergence of the method is essentially dictated by the convergence properties of the deterministic solver. We have then verified the effectiveness of the method proposed on a couple of numerical test cases, comparing its performance with other methods available in the literature. The results obtained are encouraging, as they suggest that the proposed methodology is more effective than the other methods considered here. The theoretical results have been also found to be consistent with the numerical results to a satisfactory extent. As a final remark, we observe that the methodology presented here is not limited to the spatial or temporal discretization parameters of the deterministic problem, but could also be applied to other discretization parameters, such as smoothing parameters or artificial viscosities.

#### Acknowledgement

F. Nobile and L. Tamellini received support from the Center for ADvanced MOdeling Science (CADMOS) and partial support by the Swiss National Science Foundation under the Project No. 140574 “Efficient numerical methods for flow and transport phenomena in heterogeneous random porous media”. R. Tempone is a member of the KAUST Strategic Research Initiative, Center for Uncertainty Quantification in Computational Sciences and Engineering.

Proof of Theorem \[thm:misc\_complexity\]
=========================================

The following technical lemmas are needed in the convergence proof.

\[lemma:bound-sum-with-int\] For ${{\ensuremath{{\boldsymbol x}}}} \in (1,\infty)^D$, define $\lfloor {{\ensuremath{{\boldsymbol x}}}} \rfloor = \left(\lfloor x_i\rfloor\right)_{i=1}^D$. For any $f: (1,\infty)^D \to {\mathbb{R}}$ and $g: (1,\infty)^D \to {\mathbb{R}}_+$, $$\sum_{\{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^D  \::\: f({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}) \leq 0\}} g({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}) =
    \int_{\{{{\ensuremath{{\boldsymbol x}}}} \in (1,\infty)^D  \::\: f(\lfloor {{\ensuremath{{\boldsymbol x}}}} \rfloor) \leq 0\}} g(\lfloor {{\ensuremath{{\boldsymbol x}}}} \rfloor) {~\text{d}{{\ensuremath{{\boldsymbol x}}}}}$$ holds. Moreover, if $g$ and $f$ are increasing, then $$\sum_{\{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^D  \::\: f({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}) \leq 0\}} g({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}})
    \leq
    \int_{\{{{\ensuremath{{\boldsymbol x}}}} \in (1,\infty)^D  \::\: f({{\ensuremath{{\boldsymbol x}}}} -1) \leq
      0\}} g({{\ensuremath{{\boldsymbol x}}}}) {~\text{d}{{\ensuremath{{\boldsymbol x}}}}},$$ and if $g$ and $f$ are decreasing, then $$\sum_{\{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^D  \::\: f({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}) \leq 0\}} g({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}})
    \leq
    \int_{\{{{\ensuremath{{\boldsymbol x}}}} \in (1,\infty)^D \::\: f({{\ensuremath{{\boldsymbol x}}}}) \leq
      0\}} g({{\ensuremath{{\boldsymbol x}}}}-1) {~\text{d}{{\ensuremath{{\boldsymbol x}}}}}.$$

We have $$\begin{aligned}
    \sum_{\{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^D  \::\: f({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}) \leq 0\}} g({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}})
   &= \sum_{\{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^D  \::\: f({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}) \leq 0\}}
    g({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}) \int_{{{\ensuremath{{\boldsymbol x}}}} \in [0,1]^{D}} {~\text{d}{{\ensuremath{{\boldsymbol x}}}}}\\
   &= \sum_{\{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^D  \::\: f({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}) \leq 0\}}
     \int_{{{\ensuremath{{\boldsymbol x}}}} \in [0,1]^{D}}
    g(\lfloor {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}+ {{\ensuremath{{\boldsymbol x}}}} \rfloor) {~\text{d}{{\ensuremath{{\boldsymbol x}}}}}\\
   &=
\int_{\{{{\ensuremath{{\boldsymbol x}}}} \in (1,\infty)^D  \::\: f(\lfloor {{\ensuremath{{\boldsymbol x}}}} \rfloor) \leq
      0\}} g(\lfloor {{\ensuremath{{\boldsymbol x}}}} \rfloor) {~\text{d}{{\ensuremath{{\boldsymbol x}}}}}.
  \end{aligned}$$ Combining these inequalities with $x-1 \leq \lfloor x \rfloor \leq x$ finishes the proof.

\[lemma:bound-inf-sum-exp-lin\] Assume ${{\ensuremath{{\boldsymbol a}}}} \in {\mathbb{R}}_+^D$, ${{\ensuremath{{\boldsymbol b}}}} \in {\mathbb{R}}_+^D$ and $L>|{{\ensuremath{{\boldsymbol a}}}}|$. Then, $$\sum_{\{{{\ensuremath{{\boldsymbol x}}}} \in {\mathbb{N}}_+^D \: :\: \sum_{i=1}^D a_i e^{b_i x_i} +
    b_i x_i > L \}} \exp\left(-\sum_{i=1}^D a_i
  e^{b_i x_i}\right) \leq
\left(\prod_{i=1}^D\frac{\exp(2a_i)}{a_i^2}\right)
\exp\left(-L\right) (L+1)^{2D+1}.$$

Define the set $$\mathcal P = \left\{ \left(e^{b_i x_i}\right)_{i=1}^D \::\: {{\ensuremath{{\boldsymbol x}}}}
    \in {\mathbb{N}}_+^D\right\},$$ and define $\lfloor {{\ensuremath{{\boldsymbol y}}}} \rfloor = \left(\lfloor y_i \rfloor\right)_{i=1}^D$. Then $$\aligned \sum_{\{{{\ensuremath{{\boldsymbol x}}}} \in {\mathbb{N}}_+^D \: :\: \sum_{i=1}^D a_i e^{b_ix} +
    b_i x_i > L \}} \exp\left(-\sum_{i=1}^D a_i
  e^{b_i x_i}\right)  &=
  \sum_{\{{{\ensuremath{{\boldsymbol y}}}} \in \mathcal P \: :\: {{\ensuremath{{\boldsymbol a}}}} \cdot {{\ensuremath{{\boldsymbol y}}}} +
    |\log({{\ensuremath{{\boldsymbol y}}}})| > L \}} \exp(-{{\ensuremath{{\boldsymbol a}}}} \cdot {{\ensuremath{{\boldsymbol y}}}}) \\
  &\leq \sum_{\{{{\ensuremath{{\boldsymbol y}}}} \in \mathcal P \: :\: {{\ensuremath{{\boldsymbol a}}}} \cdot (\lfloor
    {{\ensuremath{{\boldsymbol y}}}} \rfloor + {{\ensuremath{{\boldsymbol 1}}}}) +
    |\log(\lfloor {{\ensuremath{{\boldsymbol y}}}}\rfloor+{{\ensuremath{{\boldsymbol 1}}}})| > L \}} \exp(-{{\ensuremath{{\boldsymbol a}}}} \lfloor {{\ensuremath{{\boldsymbol y}}}} \rfloor) \\
&\leq \sum_{\{{{\ensuremath{{\boldsymbol y}}}} \in {\mathbb{N}}_+^D \: :\: {{\ensuremath{{\boldsymbol a}}}} \cdot  {{\ensuremath{{\boldsymbol y}}}} +
    |\log( {{\ensuremath{{\boldsymbol y}}}}+{{\ensuremath{{\boldsymbol 1}}}})| > L-|{{\ensuremath{{\boldsymbol a}}}}| \}} \exp(-{{\ensuremath{{\boldsymbol a}}}} \cdot {{\ensuremath{{\boldsymbol y}}}}) \\
&\leq \int_{\{{{\ensuremath{{\boldsymbol y}}}} \in (1, \infty)^D \: :\: {{\ensuremath{{\boldsymbol a}}}} \cdot {{\ensuremath{{\boldsymbol y}}}} +
    |\log( {{\ensuremath{{\boldsymbol y}}}}+{{\ensuremath{{\boldsymbol 1}}}})| > L-|{{\ensuremath{{\boldsymbol a}}}}| \}} \exp(-{{\ensuremath{{\boldsymbol a}}}} \cdot ({{\ensuremath{{\boldsymbol y}}}} - {{\ensuremath{{\boldsymbol 1}}}})) {~\text{d}{{\ensuremath{{\boldsymbol y}}}}}.
\endaligned$$ Letting $z_i = a_i y_i + \log(y_i+1)$ and $p(z_i) = y_i \leq
\frac{z_i}{a_i}$, then $$\begin{aligned}
  & \sum_{\{{{\ensuremath{{\boldsymbol x}}}} \in {\mathbb{N}}_+^D \: :\: \sum_{i=1}^D a_i
    e^{b_ix} + b_i x_i > L \}} \exp{\left(-\sum_{i=1}^D a_i e^{b_i
        x_i}\right)} \\
  & = \exp(|{{\ensuremath{{\boldsymbol a}}}}|)\int_{\{{{\ensuremath{{\boldsymbol y}}}} \in (1, \infty)^D \: :\: {{\ensuremath{{\boldsymbol a}}}}
    \cdot {{\ensuremath{{\boldsymbol y}}}} + |\log({{\ensuremath{{\boldsymbol y}}}}+{{\ensuremath{{\boldsymbol 1}}}})| > L-|{{\ensuremath{{\boldsymbol a}}}}| \}} \exp(-{{\ensuremath{{\boldsymbol a}}}}
  \cdot {{\ensuremath{{\boldsymbol y}}}} - |\log({{\ensuremath{{\boldsymbol y}}}}+ {{\ensuremath{{\boldsymbol 1}}}})| + |\log({{\ensuremath{{\boldsymbol y}}}}+{{\ensuremath{{\boldsymbol 1}}}})|)
  {~\text{d}{{\ensuremath{{\boldsymbol y}}}}} \\
  &= \exp(|{{\ensuremath{{\boldsymbol a}}}}|) \int_{\{{{\ensuremath{{\boldsymbol z}}}} \in \otimes_{i=1}^D(a_i + \log(2),
    \infty) \: :\: |{{\ensuremath{{\boldsymbol z}}}}| > L-|{{\ensuremath{{\boldsymbol a}}}}| \}} \exp(- |{{\ensuremath{{\boldsymbol z}}}}|)
  \prod_{i=1}^D \left( \frac{p(z_i) + 1}{a_i +
      \frac{1}{p(z_i)+1}}\right){~\text{d}{{\ensuremath{{\boldsymbol z}}}}} \\
  &\leq \exp(|{{\ensuremath{{\boldsymbol a}}}}|) \int_{\{{{\ensuremath{{\boldsymbol z}}}} \in \otimes_{i=1}^D(a_i +
    \log(2), \infty) \: :\: |{{\ensuremath{{\boldsymbol z}}}}| > L-|{{\ensuremath{{\boldsymbol a}}}}| \}} \exp(-|{{\ensuremath{{\boldsymbol z}}}}|)
  \prod_{i=1}^D \left( \frac{z_i +
      a_i}{a_i^2}\right){~\text{d}{{\ensuremath{{\boldsymbol z}}}}} \\
  &\leq \left(\prod_{i=1}^D\frac{\exp(a_i)}{a_i^2}\right) \int_{\{{{\ensuremath{{\boldsymbol z}}}} \in \otimes_{i=1}^D(a_i + \log(2), \infty) \: :\: |{{\ensuremath{{\boldsymbol z}}}}| >
    L-|{{\ensuremath{{\boldsymbol a}}}}| \}}
  \exp(-|{{\ensuremath{{\boldsymbol z}}}}| + |\log({{\ensuremath{{\boldsymbol z}}}} + {{\ensuremath{{\boldsymbol a}}}})|){~\text{d}{{\ensuremath{{\boldsymbol z}}}}} \\
  &= \left(\prod_{i=1}^D\frac{\exp(2 a_i)}{a_i^2}\right) \int_{\{{{\ensuremath{{\boldsymbol x}}}} \in \otimes_{i=1}^D(\log(2), \infty) \: :\: |{{\ensuremath{{\boldsymbol x}}}}| > L \}}
  \exp(-|{{\ensuremath{{\boldsymbol x}}}}| + |\log({{\ensuremath{{\boldsymbol x}}}})|){~\text{d}{{\ensuremath{{\boldsymbol x}}}}} \\
  &\leq \left(\prod_{i=1}^D\frac{\exp(2 a_i)}{a_i^2}\right)
  \int_{\{{{\ensuremath{{\boldsymbol z}}}} \in \otimes_{i=1}^D(0, \infty) \: :\: |{{\ensuremath{{\boldsymbol z}}}}| > L
    \}} \exp(-|{{\ensuremath{{\boldsymbol z}}}}| + |\log({{\ensuremath{{\boldsymbol z}}}})|){~\text{d}{{\ensuremath{{\boldsymbol z}}}}} .
  \end{aligned}$$ Now let us prove, by induction on $D$, that we have $$\int_{\{{{\ensuremath{{\boldsymbol z}}}} \in {\mathbb{R}}_+^D \: :\:
  |{{\ensuremath{{\boldsymbol z}}}}| > L \}}
 \exp(-|{{\ensuremath{{\boldsymbol z}}}}| + |\log({{\ensuremath{{\boldsymbol z}}}})|){~\text{d}{{\ensuremath{{\boldsymbol z}}}}} \leq \exp(-L)(L+1)^{2D-1}\,\,.$$ For $D=1$, the inequality is a trivial equality that can be obtained with integration by parts. Assume the inequality is true for $D$ and let us prove it for $D+1$: $$\begin{aligned}
\int_{\{{{\ensuremath{{\boldsymbol z}}}} \in {\mathbb{R}}_+^{D+1} \: :\: |{{\ensuremath{{\boldsymbol z}}}}| > L
    \}} \exp(-|{{\ensuremath{{\boldsymbol z}}}}| + \log({{\ensuremath{{\boldsymbol z}}}} )){~\text{d}{{\ensuremath{{\boldsymbol z}}}}} &=
  \int_{L}^\infty y \exp(-y) \int_{\{{{\ensuremath{{\boldsymbol x}}}} \in
    {\mathbb{R}}_+^{D} \}}
  \exp(-|{{\ensuremath{{\boldsymbol x}}}}| + \log({{\ensuremath{{\boldsymbol x}}}} )){~\text{d}{{\ensuremath{{\boldsymbol x}}}}} {~\text{d}y} \\
  &+\int_{0}^L y\exp(-y) \int_{\{{{\ensuremath{{\boldsymbol x}}}} \in
    {\mathbb{R}}_+^{D} \::\: |{{\ensuremath{{\boldsymbol x}}}}| > L-y \}}
  \exp(-|{{\ensuremath{{\boldsymbol x}}}}| + \log({{\ensuremath{{\boldsymbol x}}}})){~\text{d}{{\ensuremath{{\boldsymbol x}}}}} {~\text{d}y} \\
  &\leq \exp(-L) (L+1) \\
  &+ \int_{0}^L y\exp(-y) \exp(-L+y)(L-y+1)^{2D-1} {~\text{d}y}
  \\
  &\leq \exp(-L)\left( L+1+ L^2 (L+1)^{2D-1}\right)\\
&\leq \exp(-L)\left( L+1\right)\left(1+ L (L+1)^{2D-1}\right)\\
&\leq \exp(-L) (L+1)^{2(D+1)-1}\,\,.
  \end{aligned}$$ Finally, substituting back, we get the result.

Given ${{\ensuremath{{\boldsymbol a}}}} \in {\mathbb{R}}_+^D$ and $A>0$, let $\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, A)$ denote the number of occurrences of $A$ in ${{\ensuremath{{\boldsymbol a}}}}$, $$\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, A) = \#\{i = 1, \ldots, d \::\: a_i = A\}.$$

Assume $k \in {\mathbb{N}}$, ${{\ensuremath{{\boldsymbol a}}}} \in {\mathbb{R}}_+^D, L>|{{\ensuremath{{\boldsymbol a}}}}|$. Then, the following bounds hold true: $$\begin{aligned}
  \int_{\{{{\ensuremath{{\boldsymbol x}}}} \in {\mathbb{R}}_+^D \::\: |{{\ensuremath{{\boldsymbol x}}}}| > L\}} \exp(-{{\ensuremath{{\boldsymbol a}}}} \cdot{{\ensuremath{{\boldsymbol x}}}}){~\text{d}{{\ensuremath{{\boldsymbol x}}}}}
    \leq \mathfrak{B}_D({{\ensuremath{{\boldsymbol a}}}}) \exp(-\min({{\ensuremath{{\boldsymbol a}}}}) L) L^{\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, \min({{\ensuremath{{\boldsymbol a}}}}))-1}\end{aligned}$$ where $\mathfrak{B}_D({{\ensuremath{{\boldsymbol a}}}})$ is a positive constant independent of $L$.

See [@abdullatif.etal:MultiIndexMC Lemma B.3] for a proof of the inequality and the value of $\mathfrak{B}_D({{\ensuremath{{\boldsymbol a}}}})$. Moreover, a proof of a consistent equality for the case ${{\ensuremath{{\boldsymbol a}}}}={{\ensuremath{{\boldsymbol 1}}}}$ can be found in [@griebel:tensor Proposition 2.3].

Assume $k \in {\mathbb{N}}$, ${{\ensuremath{{\boldsymbol a}}}} \in {\mathbb{R}}_+^D, L>|{{\ensuremath{{\boldsymbol a}}}}|$. Then, the following bound holds: $$\begin{aligned}
  \int_{\{{{\ensuremath{{\boldsymbol x}}}} \in {\mathbb{R}}_+^D \: :\: |{{\ensuremath{{\boldsymbol x}}}}| \leq L \}}
  \exp\left({{\ensuremath{{\boldsymbol a}}}} \cdot {{\ensuremath{{\boldsymbol x}}}} \right) \left(L - |{{\ensuremath{{\boldsymbol x}}}}| \right)^k {~\text{d}{{\ensuremath{{\boldsymbol x}}}}}
  \leq \mathfrak A_{D}({{\ensuremath{{\boldsymbol a}}}}, k) \exp(\max({{\ensuremath{{\boldsymbol a}}}}) L) L^{\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, \max({{\ensuremath{{\boldsymbol a}}}}))-1},\end{aligned}$$ where $$\begin{aligned}
& \mathfrak A_{D}({{\ensuremath{{\boldsymbol a}}}}, k) = \frac{k!}{(\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, \max({{\ensuremath{{\boldsymbol a}}}}))-1)! \max({{\ensuremath{{\boldsymbol a}}}})^{k+1}}
  \left(\prod_{a_i < \max({{\ensuremath{{\boldsymbol a}}}})} \frac{1}{ \max({{\ensuremath{{\boldsymbol a}}}}) -
  a_i}\right).
$$

Without loss of generality, assume that $a_{i} \geq a_{i+1}$ for all $i=1 \ldots D$, such that $a_1 = \max({{\ensuremath{{\boldsymbol a}}}})$. We prove the result by induction on $D$. For $D=1$, we have $$\aligned \int_0^L \exp\left(a x \right) \left(L - x \right)^k {~\text{d}x}
  &= \frac{k!}{a^{k+1}}\left(\exp(aL) - \sum_{i=0}^{k} \frac{a^i L^{i}}{i!}\right) \\
&\leq \frac{k! \exp(aL)}{a^{k+1}}.
\endaligned$$

Next, assume that the result is valid for a given $D>1$ and ${{\ensuremath{{\boldsymbol a}}}}\in{\mathbb{R}}_+^D$ where $a_{i} \geq a_{i+1}$ for all $i=1 \ldots D$, such that $a_1 = \max({{\ensuremath{{\boldsymbol a}}}})$. Let $b \leq a_1$ and define a new vector $\widetilde{{{\ensuremath{{\boldsymbol a}}}}} = ({{\ensuremath{{\boldsymbol a}}}}, b) \in
{\mathbb{R}}_+^{D+1}$. We have $$\begin{aligned}
& \int_{\{({{\ensuremath{{\boldsymbol x}}}}, y) \in {\mathbb{R}}_+^{D+1} \: :\: y +
    |{{\ensuremath{{\boldsymbol x}}}}| \leq L \}} \exp\left(b y + {{\ensuremath{{\boldsymbol a}}}} \cdot {{\ensuremath{{\boldsymbol x}}}}
  \right)
  \left(L - y - |{{\ensuremath{{\boldsymbol x}}}}| \right)^k {~\text{d}y} {~\text{d}{{\ensuremath{{\boldsymbol x}}}}} \\
  &=\int_0^L {\exp\left(b y \right)} \int_{\{{{\ensuremath{{\boldsymbol x}}}} \in
    {\mathbb{R}}_+^{D} \: :\: |{{\ensuremath{{\boldsymbol x}}}}| \leq L-y \}} \exp\left({{\ensuremath{{\boldsymbol a}}}} \cdot
    {{\ensuremath{{\boldsymbol x}}}} \right) \left(L - y -|{{\ensuremath{{\boldsymbol x}}}}| \right)^k {~\text{d}{{\ensuremath{{\boldsymbol x}}}}}
  {~\text{d}y} \\
  &\leq \mathfrak{A}_D({{\ensuremath{{\boldsymbol a}}}}, k) {\exp(a_1 L)} \int_0^L {\exp\left((b-a_1) y \right)} (L-y)^{\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, a_1)-1}
  {~\text{d}y}.
  \end{aligned}$$ We distinguish between two cases:

1.  $b < a_1$ then $\mathfrak{n}(\widetilde{{{\ensuremath{{\boldsymbol a}}}}}, a_1) = \mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, a_1)$ and $$\aligned \int_0^L {\exp\left(-(a_1-b)
          y \right)} (L-y)^{\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, a_1)-1} {~\text{d}y}
      &\leq L^{\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, a_1)-1}
      \int_0^\infty {\exp\left(-(a_1-b) y \right)}  {~\text{d}y} \\
    &\leq L^{\mathfrak{n}(\widetilde{{{\ensuremath{{\boldsymbol a}}}}}, a_1)-1} \frac{1}{a_1-b},
    \endaligned$$ and in this case $$\aligned \mathfrak{A}_D({{\ensuremath{{\boldsymbol a}}}}, k) \frac{1}{a_1-b} &=
    \frac{k!}{(\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, a_1)-1)! a_1^{k+1}}
      \left(\prod_{a_i < a_1} \frac{1}{ a_1 - a_i}\right) \left(\frac{1}{a_1-b}\right)\\
    &= \mathfrak{A}_{D+1}(\widetilde{{{\ensuremath{{\boldsymbol a}}}}}, k).
    \endaligned$$

2.  $b = a_1$ then $\mathfrak{n}(\widetilde{{{\ensuremath{{\boldsymbol a}}}}}, a_1) = \mathfrak{n}({{\ensuremath{{\boldsymbol a}}}},
      a_1)+1$ and $$\aligned
     \int_0^L (L-y)^{\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, a_1)-1}
      {~\text{d}y} &= \frac{L^{\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, a_1)}}{\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, a_1)}
    = \frac{L^{\mathfrak{n}(\widetilde{{{\ensuremath{{\boldsymbol a}}}}}, a_1)-1}}{\mathfrak{n}(\widetilde{{{\ensuremath{{\boldsymbol a}}}}}, a_1)-1},
    \endaligned$$ and again $$\aligned \frac{1}{\mathfrak{n}(\widetilde{{{\ensuremath{{\boldsymbol a}}}}}, a_1)-1}  \mathfrak{A}_D({{\ensuremath{{\boldsymbol a}}}}, k)  &=
    \frac{1}{\mathfrak{n}(\widetilde{{{\ensuremath{{\boldsymbol a}}}}}, a_1)-1}  \cdot \frac{k!}{(\mathfrak{n}({{\ensuremath{{\boldsymbol a}}}}, a_1)-1)! a_1^{k+1}}
      \left(\prod_{a_i < a_1} \frac{1}{ a_1 - a_i}\right)\\
    &= \mathfrak{A}_{D+1}(\widetilde{{{\ensuremath{{\boldsymbol a}}}}}, k).
    \endaligned$$

The bounds and can be obtained by elementary algebraic operations combining Assumptions \[assump:growth\_of\_dof\] and \[assump:dW\_dE\_model\]; for instance, $$\begin{aligned}
\Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}^\textnormal{det} &
\leq C_{\textnormal{work}}^{\textnormal{det}} \prod_{i=1}^D h_i^{-\widetilde{\gamma}_i}
= C_{\textnormal{work}}^{\textnormal{det}} \prod_{i=1}^D (h_0 2^{-\alpha_i})^{-\widetilde{\gamma}_i}
= C_{\textnormal{work}}^{\textnormal{det}} h_0^{-|\widetilde{{{\ensuremath{{\boldsymbol \gamma}}}}}|} \prod_{i=1}^D 2^{\widetilde{\gamma}_i \alpha_i}
= C_{\textnormal{work}}^{\textnormal{det}} h_0^{-|\widetilde{{{\ensuremath{{\boldsymbol \gamma}}}}}|} \prod_{i=1}^D e^{\widetilde{\gamma}_i \alpha_i \log 2}, \\
\Delta W_{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}^\textnormal{stoc}  &
\leq C_{\textnormal{work}}^{\textnormal{stoc}} \prod_{n=1}^N 2^{\beta_n}
= C_{\textnormal{work}}^{\textnormal{stoc}} \prod_{n=1}^N e^{\beta_n \log 2},\end{aligned}$$ from which follows by setting ${\ensuremath{C_{\text{\textnormal{work}}}}}= C_{\textnormal{work}}^{\textnormal{det}}
h_0^{-|\widetilde{{{\ensuremath{{\boldsymbol \gamma}}}}}|}
C_{\textnormal{work}}^{\textnormal{stoc}}$. The proof is then divided into two steps.

#### Step 1: Work Estimate

Observe that $\Xi_i = \frac{\gamma_{i}}{\gamma_i+{{r_{i}}}}  < 1$ for all $i=1,\ldots D$ and that $\mathfrak{z} = \mathfrak{n}({{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}, \chi)$. Thanks to equations and , and using Lemma \[lemma:bound-sum-with-int\], the total work satisfies $$\begin{aligned}
  &{{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathcal{I}^*(L)\right]}} =
\sum_{({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}) \in \mathcal{I}^*(L)} \Delta
  W_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \\
  \leq & {{\ensuremath{C_{\text{\textnormal{work}}}}}}
         \sum_{\left\{
         ({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}) \in {\mathbb{N}}_+^{D+N} \::\:
         \sum_{i=1}^D ({{r_{i}}} + \gamma_i)\alpha_i
         + \sum_{j=1}^N \delta \beta_j + g_j e^{\delta \beta_j} \leq
          L
         \right\}}
         \exp\left(\sum_{i=1}^D \gamma_i \alpha_i + \delta |{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}|
         \right) \\
  \leq & {{\ensuremath{C_{\text{\textnormal{work}}}}}}
         \int_{\left\{({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}) \in (1,\infty)^{D+N} \::\:
         \sum_{i=1}^D ({{r_{i}}} + \gamma_i)(\alpha_i-1)
         + \sum_{j=1}^N \delta (\beta_j-1) + g_j e^{\delta (\beta_j-1)} \leq
          L
         \right\}}
         \exp\left(\sum_{i=1}^D \gamma_i \alpha_i + \delta |{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}|
         \right) {~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}{~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}.\end{aligned}$$ Next, let $\overline{\beta_j} = g_j e^{\delta(\beta_j-1)}$ and $\overline{\alpha_i} = ({{r_{i}}} + \gamma_i) (\alpha_i-1)$. We have $$\begin{aligned}
  {{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathcal{I}^*(L)\right]}} \leq
   & {{\ensuremath{C_{\text{\textnormal{work}}}}}} \left(\prod_{j=1}^N \frac{2}{g_j \delta}\right)
      \left(\prod_{i=1}^D \frac{\exp(\gamma_i)}{{{r_{i}}} + \gamma_i}
      \right) \\
      &\qquad \int_{\left\{ (\overline{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}, \overline{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}) \in
          {\mathbb{R}}_+^D \times (\otimes_{j=1}^N (g_j, \infty))
          \::\:
          |\overline{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}|
          + |\overline{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| + |\log \overline{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| \leq
          L + |\log {{\ensuremath{{\boldsymbol g}}}}| \right\}}
      \exp\left({{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}\cdot \overline {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\right) {~\text{d}\overline{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}}{~\text{d}\overline{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}.
    \end{aligned}$$ Dropping the over-line notation and defining $\widetilde L = L + |\log
{{\ensuremath{{\boldsymbol g}}}}|$ and $\mathscr{C}_{\text{W},1}$ to be the constant factor, we obtain $$\begin{aligned}
  &{{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathcal{I}^*(L)\right]}}
    \leq \mathscr{C}_{\text{W},1}
    \int_{\left\{ ({{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}, {{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}) \in
{\mathbb{R}}_+^D \times (\otimes_{j=1}^N (g_j, \infty)) \::\:
    |{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}|
    + |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| + |\log {{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| \leq  \widetilde L \right\}}
    \exp\left({{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}\cdot  {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\right) {~\text{d}{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}}{~\text{d}{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}}\\
  &= \mathscr{C}_{\text{W},1}
    \int_{\left\{ {{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \in \otimes_{j=1}^N (g_j, \infty) \::\:
    |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| + |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| \leq  \widetilde L \right\}}  \int_{\left\{
    {{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}} \in {\mathbb{R}}_+^D \::\:
    |{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}| \leq  \widetilde L - |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| - |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| \right\}}
    \exp\left({{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}\cdot  {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\right) {~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}} {~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \\
  &\leq \mathscr{C}_{\text{W},1}
    \mathfrak{A}_D\left( {{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}, 0\right)
    \int_{\left\{  {{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \in \otimes_{j=1}^N (g_j, \infty) \::\:
    |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| + |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| \leq  \widetilde L \right\}} \exp\Big( \chi \left(\widetilde L - |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| -
    |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}|\right)\Big) \left(\widetilde L - |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| -
    |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}|\right)^{\mathfrak{z}-1} {~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}\end{aligned}$$ Define $\mathscr{C}_{\text{W},2} = \mathscr{C}_{\text{W},1}
    \mathfrak{A}_D\left( {{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}, 0\right) \exp(\chi \widetilde L)$, then $$\begin{aligned}
{{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathcal{I}^*(L)\right]}}
  &\leq \mathscr{C}_{\text{W},2}
    \int_{\left\{  {{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \in \otimes_{j=1}^N (g_j, \infty) \::\:
    |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| + |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| \leq  \widetilde L \right\}} \exp\left(- \chi \left(|{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| +
    |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}|\right)\right) \left(\widetilde L - |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| -
    |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}|\right)^{\mathfrak{z}-1} {~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}\\
  &\leq \mathscr{C}_{\text{W},2}
     \left({\widetilde L} - |{{\ensuremath{{\boldsymbol g}}}}| -
    |\log {{\ensuremath{{\boldsymbol g}}}} |\right)^{\mathfrak{z}-1}
    \int_{\left\{  {{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \in \otimes_{j=1}^N (g_j, \infty) \::\:
    |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| + |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| \leq  \widetilde L \right\}} \exp\left(-\chi \left( |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| +
    |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}|\right)\right)  {~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}.\end{aligned}$$ Since $\chi > 0$, the previous integral is bounded for all $\widetilde L$ and we have $$\begin{aligned}
        {{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathcal{I}^*(L)\right]}} &\leq \mathscr{C}_{\text{W}} \exp(\chi
                            L) \left( L - |{{\ensuremath{{\boldsymbol g}}}}|\right)^{\mathfrak{z}-1}\leq \mathscr{C}_{\text{W}} \exp(\chi
                            L) L^{\mathfrak{z}-1},\end{aligned}$$ where $$\mathscr{C}_{\text{W}} = {{\ensuremath{C_{\text{\textnormal{work}}}}}} \left(\prod_{j=1}^N \frac{2
    g_i^{\chi}}{g_j \log 2}\right)
      \left(\prod_{i=1}^D \frac{\exp(\gamma_i)}{{{r_{i}}} + \gamma_i} \right)
\mathfrak{A}_D\left( {{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}, 0\right)
\int_{\left\{{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}} \in \otimes_{j=1}^N (g_j, \infty) \right\}} \exp\left(-\chi \left( |{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}| +
  |\log{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}|\right)\right)  {~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}.$$ Substituting yields $$\begin{aligned}
  {{\ensuremath{\mathrm{Work}}\mspace{-2mu}\left[\mathcal{I}^*(L)\right]}}
    &\leq W_{\max} \left(1 - \frac{(\mathfrak{z}-1)
\log \left( \frac{\log\left(\frac{W_{\max}}{\mathscr{C}_{\text{W}}}\right)}{\chi}
      \right)
      }
{\log\left(\frac{W_{\max}} {\mathscr{C}_{\text{W}}}\right)}
\right)^{\mathfrak{z}-1}.
  \end{aligned}$$ From here it is easy to see that if is satisfied, then follows.

#### Step 2: Error Estimate

Thanks to equations and , the total error satisfies $$\begin{aligned}
    &{{\ensuremath{\mathrm{Error}}\mspace{-2mu}\left[\mathcal{I}^*(L)\right]}}
        \leq \sum_{({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}) \notin \mathcal I^*} \Delta E_{{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}}\\
    \leq & {\ensuremath{C_{\text{\textnormal{error}}}}}\sum_{\left\{ ({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}, {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}) \in {\mathbb{N}}_+^{D+N} \::\:
           \sum_{i=1}^D ({{r_{i}}} + \gamma_i)\alpha_i
           + \sum_{j=1}^N \delta \beta_j + g_j e^{\delta{\beta_j}} > L \right\}}
           \exp\left(
           -\sum_{i=1}^D {{r_{i}}} \alpha_i - \sum_{j=1}^N g_j
           e^{\delta \beta_j} \right) \\
    = & {\ensuremath{C_{\text{\textnormal{error}}}}}\sum_{\left\{ ({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}) \in {\mathbb{N}}_+^{D+N}
        \::\:
        \sum_{i=1}^D ({{r_{i}}} + \gamma_i)\alpha_i > L \right\}}
        \exp\left(
        -\sum_{i=1}^D {{r_{i}}} \alpha_i - \sum_{j=1}^N g_j
        e^{\delta \beta_j} \right) \\
         &+ {\ensuremath{C_{\text{\textnormal{error}}}}}\sum_{\left\{  {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^{D} \::\:
           \sum_{i=1}^D ({{r_{i}}} + \gamma_i)\alpha_i  \leq L
           \right\}}
                      \sum_{\left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}\in {\mathbb{N}}_+^{N} \::\:
           \sum_{j=1}^N \delta \beta_j + g_j e^{\delta{\beta_j}} > L
           - \sum_{i=1}^D ({{r_{i}}} + \gamma_i)\alpha_i
           \right\}}
           \exp\left(-\sum_{i=1}^D {{r_{i}}} \alpha_i -\sum_{j=1}^N g_j e^{\delta \beta_j} \right).
   \end{aligned}$$ Looking at the first term, let $\eta_i = \frac{{{r_{i}}}}{\gamma_i+{{r_{i}}}} < 1$ and ${{\ensuremath{{{\ensuremath{{\boldsymbol \eta}}}}}}}=
\left(\eta_i\right)_{i=1}^D$ and note that $\mathfrak{z} = \#\left\{ i = 1 \ldots D:  \eta_i =
   \min({{\ensuremath{{{\ensuremath{{\boldsymbol \eta}}}}}}})\right\}$. Then $$\begin{aligned}
  &\sum_{\left\{ ({{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}},{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}) \in {\mathbb{N}}_+^{D+N}
    \::\:
    \sum_{i=1}^D ({{r_{i}}} + \gamma_i)\alpha_i > L \right\}}
    \exp\left(
    -\sum_{i=1}^D {{r_{i}}} \alpha_i - \sum_{j=1}^N g_j
    e^{\delta \beta_j} \right)
  \\
  &= \left(\sum_{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}\in {\mathbb{N}}_+^{N}}
    \exp\left( - \sum_{j=1}^N g_j e^{\delta \beta_j} \right)\right)
    \left(\sum_{
    \left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^{D}
    \::\:
    \sum_{i=1}^D ({{r_{i}}} + \gamma_i)\alpha_i > L \right\}}
    \exp\left(-\sum_{i=1}^D {{r_{i}}} \alpha_i \right)\right)\\
  &\leq \mathscr{C}_{\text{E},1}
\int_{
    \left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in (1,\infty)^{D}\::\:
    \sum_{i=1}^D ({{r_{i}}} + \gamma_i)\alpha_i > L \right\}}
    \exp\left(-\sum_{i=1}^D {{r_{i}}} \left(\alpha_i-1\right) \right) {~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}\\
  &= \mathscr{C}_{\text{E},1} \left(\prod_{i=1}^D
\frac{\exp({{r_{i}}})}{{{r_{i}}} + \gamma_i} \right)
\int_{
    \left\{ {{\ensuremath{{\boldsymbol x}}}} \in \otimes_{i=1}^D ({{r_{i}}}+\gamma_i,\infty) \::\:
    |{{\ensuremath{{\boldsymbol x}}}}| > L \right\}}
    \exp\left(-\sum_{i=1}^D \frac{{{r_{i}}}}{{{r_{i}}}+\gamma_i} x_i \right)
    {~\text{d}{{\ensuremath{{\boldsymbol x}}}}} \\
&\leq \mathscr{C}_{\text{E},2} \exp\left(-\min({{\ensuremath{{{\ensuremath{{\boldsymbol \eta}}}}}}}) L\right) L^{\mathfrak{z}-1},\end{aligned}$$ where $$\mathscr{C}_{\text{E},2} =  \mathfrak{B}_D({{\ensuremath{{{\ensuremath{{\boldsymbol \eta}}}}}}})
  \left(\prod_{i=1}^D \frac{\exp({{r_{i}}})}{{{r_{i}}} + \gamma_i}
  \right) \sum_{{{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}\in {\mathbb{N}}_+^{N}}
    \exp\left( - \sum_{j=1}^N g_j e^{\delta \beta_j} \right).$$ For the second term, letting $H = L - \sum_{i=1}^D ({{r_{i}}}+\gamma_i)\alpha_i$, we can bound the sum using Lemma \[lemma:bound-inf-sum-exp-lin\]: $$\sum_{\left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}\in {\mathbb{N}}_+^{N} \::\:
           \sum_{j=1}^N \delta \beta_j + g_j e^{\delta{\beta_j}} > H \right\}}
           \exp\left(-\sum_{j=1}^N g_j e^{\delta \beta_j} \right)
\leq \left(\prod_{j=1}^N \frac{\exp(2g_j)}{g_j^2} \right) \exp(-H)
(H+1)^{2N-1}.$$ Defining $\mathscr C_{E,3} = \prod_{j=1}^N {\exp(2g_j)}{g_j^{-2}}$ and substituting back $$\begin{aligned}
 &\sum_{\left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^{D} \::\: \sum_{i=1}^D
      ({{r_{i}}} + \gamma_i)\alpha_i \leq L \right\}}
  \exp\left(-\sum_{i=1}^D {{r_{i}}} \alpha_i \right)
  \sum_{\left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \beta}}}}}}}\in {\mathbb{N}}_+^{N} \::\: \sum_{j=1}^N \delta
      \beta_j + g_j e^{\delta{\beta_j}} > L - \sum_{i=1}^D ({{r_{i}}}
      + \gamma_i)\alpha_i \right\}}
  \exp\left(-\sum_{j=1}^N g_j e^{\delta \beta_j} \right) \\
  &\leq \mathscr C_{E,3}\sum_{\left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{N}}_+^{D} \::\: \sum_{i=1}^D
      ({{r_{i}}} + \gamma_i)\alpha_i \leq L \right\}}
  \exp\left(-L + \sum_{i=1}^D \gamma_i \alpha_i \right)
  \left(L +1 - \sum_{i=1}^D ({{r_{i}}}+\gamma_i)\alpha_i\right)^{2N-1}\\
  &= \mathscr C_{E,3}  \int_{\left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in (1,\infty)^{D} \::\: \sum_{i=1}^D
      ({{r_{i}}} + \gamma_i)\lfloor \alpha_i \rfloor \leq L \right\}}
  \exp\left(-L + \sum_{i=1}^D \gamma_i \lfloor \alpha_i \rfloor \right)
  \left(L +1- \sum_{i=1}^D ({{r_{i}}}+\gamma_i)\lfloor \alpha_i
    \rfloor\right)^{2N-1}{~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}\\
&\leq \mathscr C_{E,3}    \int_{\left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in (1,\infty)^{D} \::\: \sum_{i=1}^D
      ({{r_{i}}} + \gamma_i)(\alpha_i-1) \leq L \right\}}
  \exp\left(-L + \sum_{i=1}^D \gamma_i  \alpha_i \right)
\left(L+1 - \sum_{i=1}^D ({{r_{i}}}+\gamma_i)( \alpha_i-1)\right)^{2N-1}
{~\text{d}{{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}}\\
&= \mathscr C_{E,3}\left(\prod_{i=1}^D \frac{\exp(\gamma_i)}{\gamma_i + {{r_{i}}}} \right)
  {\exp ( - L )}
  \int_{\left\{ {{\ensuremath{{{\ensuremath{{\boldsymbol \alpha}}}}}}}\in {\mathbb{R}}_+^{D} \::\: |{{\ensuremath{{\boldsymbol x}}}}| \leq L \right\}}
  \exp\left({{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}\cdot {{\ensuremath{{\boldsymbol x}}}} \right)
\left(L+1 - |{{\ensuremath{{\boldsymbol x}}}}|\right)^{2N-1}
{~\text{d}{{\ensuremath{{\boldsymbol x}}}}}\\
  &\leq \mathscr{C}_{\text{E}, 4} {\exp ( (\chi-1) L )}
  L^{\mathfrak{z}-1},\end{aligned}$$ where $$\mathscr{C}_{\text{E}, 4} = \left(\prod_{j=1}^N
    \frac{\exp(2g_j)}{g_j^2} \right) \left(\prod_{i=1}^D \frac{\exp(\gamma_i)}{\gamma_i + {{r_{i}}}} \right) \mathfrak{A}_D({{\ensuremath{{{\ensuremath{{\boldsymbol \Xi}}}}}}}, 2N-1).$$ Finally, noting that $$\chi-1 = -\min({{\ensuremath{{{\ensuremath{{\boldsymbol \eta}}}}}}}),$$ we have the error estimate $$\begin{aligned}
          {{\ensuremath{\mathrm{Error}}\mspace{-2mu}\left[\mathcal{I}^*(L)\right]}} \leq
  {\ensuremath{C_{\text{\textnormal{error}}}}}\left(\mathscr{C}_{\text{E}, 2}
  +\mathscr{C}_{\text{E}, 4} \right) \exp(-\min({{\ensuremath{{{\ensuremath{{\boldsymbol \eta}}}}}}}) L) L^{\mathfrak{z}-1}.\end{aligned}$$ Then, substituting $L$ from and evaluating the limit gives .

[10]{} url \#1[`#1`]{}urlprefixhref \#1\#2[\#2]{} \#1[\#1]{}

R. G. Ghanem, P. D. Spanos, Stochastic finite elements: a spectral approach, Springer-Verlag, New York, 1991.

O. P. Le Ma[î]{}tre, O. M. Knio, Spectral methods for uncertainty quantification, Scientific Computation, Springer, New York, 2010, with applications to computational fluid dynamics. [](http://dx.doi.org/10.1007/978-90-481-3520-2).

H. G. Matthies, A. Keese, Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations, Computer Methods in Applied Mechanics and Engineering 194 (12-16) (2005) 1295–1331.

R. A. Todor, C. Schwab, [Convergence rates for sparse chaos approximations of elliptic problems with stochastic coefficients]{}, IMA J Numer Anal 27 (2) (2007) 232–261. [](http://dx.doi.org/10.1093/imanum/drl025).

D. Xiu, G. Karniadakis, The [W]{}iener-[A]{}skey polynomial chaos for stochastic differential equations, [SIAM]{} Journal on Scientific Computing 24 (2) (2002) 619–644.

I. Babuška, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, [SIAM]{} Review 52 (2) (2010) 317–355.

H. J. Bungartz, M. Griebel, Sparse grids, Acta Numerica 13 (2004) 147–269.

F. Nobile, R. Tempone, C. Webster, An anisotropic sparse grid stochastic collocation method for partial differential equations with random input data, [SIAM]{} Journal on Numerical Analysis 46 (5) (2008) 2411–2442.

D. Xiu, J. Hesthaven, High-order collocation methods for differential equations with random inputs, [SIAM]{} Journal on Scientific Computing 27 (3) (2005) 1118–1139.

B. Khoromskij, C. Schwab, Tensor-structured galerkin approximation of parametric and stochastic elliptic pdes, [SIAM]{} Journal on Scientific Computing 33 (1) (2011) 364–385.

B. Khoromskij, I. Oseledets, Quantics-tt collocation approximation of parameter-dependent and stochastic elliptic pdes, Computational Methods in Applied Mathematics 10 (4) (2010) 376–394.

A. Nouy, Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms, Computer Methods in Applied Mechanics and Engineering 197 (51) (2008) 4718–4736.

J. Ballani, L. Grasedyck, [H]{}ierarchical [T]{}ensor [A]{}pproximation of [O]{}utput [Q]{}uantities of [P]{}arameter-[D]{}ependent [PDE]{}s, [SIAM/ASA]{} Journal on Uncertainty Quantification 3 (1) (2015) 852–872. [](http://dx.doi.org/10.1137/140960980).

S. Boyaval, C. Le Bris, T. Lelièvre, Y. Maday, N. Nguyen, A. Patera, Reduced basis techniques for stochastic problems, Archives of Computational Methods in Engineering 17 (4) (2010) 435–454.

P. Chen, A. Quarteroni, G. Rozza, Comparison between reduced basis and stochastic collocation methods for elliptic problems, Journal of Scientific Computing 59 (1) (2014) 187–216. [](http://dx.doi.org/10.1007/s10915-013-9764-2).

A. Cohen, R. Devore, C. Schwab, Analytic regularity and polynomial approximation of parametric and stochastic elliptic [PDE’S]{}, Analysis and Applications 9 (1) (2011) 11–47.

S. Heinrich, Multilevel [M]{}onte [C]{}arlo methods, in: Large-Scale Scientific Computing, Vol. 2179 of Lecture Notes in Computer Science, Springer Berlin Heidelberg, 2001, pp. 58–67.

M. B. Giles, Multilevel [M]{}onte [C]{}arlo path simulation, Operations Research 56 (3) (2008) 607–617.

A. Barth, C. Schwab, N. Zollinger, Multi-level [M]{}onte [C]{}arlo finite element method for elliptic [PDE]{}s with stochastic coefficients, Numerische Mathematik 119 (1) (2011) 123–161.

A. Barth, A. Lang, C. Schwab, Multilevel [M]{}onte [C]{}arlo method for parabolic stochastic partial differential equations, [BIT]{} Numerical Mathematics 53 (1) (2013) 3–27.

J. Charrier, R. Scheichl, A. Teckentrup, Finite element error analysis of elliptic pdes with random coefficients and its application to multilevel [M]{}onte [C]{}arlo methods, [SIAM]{} Journal on Numerical Analysis 51 (1) (2013) 322–352.

K. Cliffe, M. Giles, R. Scheichl, A. Teckentrup, Multilevel [M]{}onte [C]{}arlo methods and applications to elliptic pdes with random coefficients, Computing and Visualization in Science 14 (1) (2011) 3–15. [](http://dx.doi.org/10.1007/s00791-011-0160-x).

S. Mishra, C. Schwab, J. Sukys, Multi-level [M]{}onte [C]{}arlo finite volume methods for nonlinear systems of conservation laws in multi-dimensions, Journal of Computational Physics 231 (8) (2012) 3365–3388.

A. Teckentrup, P. Jantsch, C. G. Webster, M. Gunzburger, A [M]{}ultilevel [S]{}tochastic [C]{}ollocation [M]{}ethod for [P]{}artial [D]{}ifferential [E]{}quations with [R]{}andom [I]{}nput [D]{}ata, [SIAM/ASA]{} Journal on Uncertainty Quantification 3 (1) (2015) 1046–1074. [](http://dx.doi.org/10.1137/140969002).

H. W. van Wyk, Multilevel sparse grid methods for elliptic partial differential equations with random coefficients, arXiv preprint arXiv:1404.0963 (2014).

H. Harbrecht, M. Peters, M. Siebenmorgen, On multilevel quadrature for elliptic stochastic partial differential equations, in: Sparse Grids and Applications, Vol. 88 of Lecture Notes in Computational Science and Engineering, Springer, 2013, pp. 161–179.

F. Y. Kuo, C. Schwab, I. Sloan, Multi-level [Q]{}uasi-[M]{}onte [C]{}arlo [F]{}inite [E]{}lement [M]{}ethods for a [C]{}lass of [E]{}lliptic [PDE]{}s with [R]{}andom [C]{}oefficients, Foundations of Computational Mathematics 15 (2) (2015) 411–449.

F. Nobile, F. Tesei, A [M]{}ulti [L]{}evel [M]{}onte [C]{}arlo method with control variate for elliptic [PDE]{}s with log-normal coefficients, Stochastic Partial Differential Equations: Analysis and Computations 3 (3) (2015) 398–444. [](http://dx.doi.org/10.1007/s40072-015-0055-9).

A.-L. Haji-Ali, F. Nobile, R. Tempone, [M]{}ulti-[i]{}ndex [M]{}onte [C]{}arlo: when sparsity meets sampling, Numerische Mathematik 132 (2015) 767–806. [](http://dx.doi.org/10.1007/s00211-015-0734-5).

H. J. Bungartz, M. Griebel, D. Röschke, C. Zenger, Pointwise convergence of the combination technique for the [L]{}aplace equation, East-West Journal of Numerical Mathematics 2 (1994) 21–45.

M. Griebel, M. Schneider, C. Zenger, A combination technique for the solution of sparse grid problems, in: P. de Groen, R. Beauwens (Eds.), [Iterative Methods in Linear Algebra]{}, IMACS, Elsevier, North Holland, 1992, pp. 263–281.

M. Hegland, J. Garcke, V. Challis, The combination technique and some generalisations, Linear Algebra and its Applications 420 (2–3) (2007) 249–275. [](http://dx.doi.org/10.1016/j.laa.2006.07.014).

M. Griebel, H. Harbrecht, On the convergence of the combination technique, in: J. Garcke, D. Pflüger (Eds.), Sparse Grids and Applications - Munich 2012, Vol. 97 of Lecture Notes in Computational Science and Engineering, Springer International Publishing, 2014, pp. 55–74. [](http://dx.doi.org/10.1007/978-3-319-04537-5_3).

F. Nobile, L. Tamellini, R. Tempone, Convergence of quasi-optimal sparse-grid approximation of [H]{}ilbert-space-valued functions: application to random elliptic [PDE]{}s, Numerische [M]{}athematik[](http://dx.doi.org/10.1007/s00211-015-0773-y).

J. Beck, F. Nobile, L. Tamellini, R. Tempone, On the optimal polynomial approximation of stochastic [PDE]{}s by [G]{}alerkin and collocation methods, Mathematical Models and Methods in Applied Sciences 22 (09) (2012) 1250023.

J. Beck, F. Nobile, L. Tamellini, R. Tempone, A [Q]{}uasi-optimal [S]{}parse [G]{}rids [P]{}rocedure for [G]{}roundwater [F]{}lows, in: Spectral and High Order Methods for Partial Differential Equations - ICOSAHOM 2012, Vol. 95 of Lecture Notes in Computational Science and Engineering, Springer, 2014, pp. 1–16.

M. Griebel, S. Knapek, Optimized general sparse grid approximation spaces for operator equations, Mathematics of Computation 78 (268) (2009) 2223–2257. [](http://dx.doi.org/10.1090/S0025-5718-09-02248-0).

M. Bieri, A sparse composite collocation finite element method for elliptic [SPDE]{}s., [SIAM]{} Journal on Numerical Analysis 49 (6) (2011) 2277–2301.

T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: [CAD]{}, finite elements, nurbs, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (39–41) (2005) 4135–4195. [](http://dx.doi.org/10.1016/j.cma.2004.10.008).

W. J. Gordon, C. A. Hall, Construction of curvilinear co-ordinate systems and applications to mesh generation, International Journal for Numerical Methods in Engineering 7 (4) (1973) 461–477. [](http://dx.doi.org/10.1002/nme.1620070405).

A. Quarteroni, A. Valli, Domain Decomposition Methods for Partial Differential Equations, Numerical mathematics and scientific computation, Clarendon Press, 1999.

L. Trefethen, Approximation Theory and Approximation Practice, Society for Industrial and Applied Mathematics, 2013.

L. N. Trefethen, Is [G]{}auss quadrature better than [C]{}lenshaw-[C]{}urtis?, [SIAM]{} Review 50 (1) (2008) 67–87.

A. Chkifa, On the lebesgue constant of leja sequences for the complex unit disk and of their real projection, Journal of Approximation Theory 166 (0) (2013) 176–200.

F. Nobile, L. Tamellini, R. Tempone, Comparison of [C]{}lenshaw-[C]{}urtis and [L]{}eja [Q]{}uasi-[O]{}ptimal [S]{}parse [G]{}rids for the [A]{}pproximation of [R]{}andom [PDE]{}s, in: R. M. Kirby, M. Berzins, J. S. Hesthaven (Eds.), Spectral and High Order Methods for Partial Differential Equations - ICOSAHOM ’14, Vol. 106 of Lecture Notes in Computational Science and Engineering, Springer International Publishing, 2015, pp. 475–482. [](http://dx.doi.org/10.1007/978-3-319-19800-2_44).

A. Narayan, J. D. Jakeman, Adaptive [L]{}eja [S]{}parse [G]{}rid [C]{}onstructions for [S]{}tochastic [C]{}ollocation and [H]{}igh-[D]{}imensional [A]{}pproximation, [SIAM]{} Journal on Scientific Computing 36 (6) (2014) A2952–A2983.

A. Genz, B. D. Keister, Fully symmetric interpolatory rules for multiple integrals over infinite regions with [G]{}aussian weight, Journal of Computational and Applied Mathematics 71 (2) (1996) 299–309.

G. W. Wasilkowski, H. Wozniakowski, Explicit cost bounds of algorithms for multivariate tensor product problems, Journal of Complexity 11 (1) (1995) 1–56.

S. Martello, P. Toth, Knapsack problems: algorithms and computer implementations, Wiley-Interscience series in discrete mathematics and optimization, J. Wiley & Sons, 1990.

D. D[ũ]{}ng, M. Griebel, Hyperbolic cross approximation in infinite dimensions, Journal of Complexity[](http://dx.doi.org/10.1016/j.jco.2015.09.006).

M. Griebel, J. Oettershagen, On tensor product approximation of analytic functions, Journal of Approximation Theory (2016) In print[](http://dx.doi.org/10.1016/j.jat.2016.02.006).

J. Bäck, F. Nobile, L. Tamellini, R. Tempone, Stochastic spectral [G]{}alerkin and collocation methods for [PDE]{}s with random coefficients: a numerical comparison, in: Spectral and High Order Methods for Partial Differential Equations, Vol. 76 of Lecture Notes in Computational Science and Engineering, Springer, 2011, pp. 43–62.
