\section{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
\cite{ghanem.spanos:book,lemaitre:book,matthies.keese:galerkin,todor.schwab:convergence,xiu.karniadakis:wiener}
or collocation strategies based on sparse grids
(see e.g. \cite{babuska.nobile.eal:stochastic2,b.griebel:acta,nobile.tempone.eal:aniso,xiu.hesthaven:high}),
low-rank techniques \cite{khoromskij.schwab:tensor,khoromskij.oseledets:tensor,nouy:2008a,ballani.gras:tensor}
and reduced basis methods (see e.g. \cite{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., \cite{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
\emph{Multilevel Monte Carlo} (MLMC). It was first
proposed in \cite{Heinrich:MLMC} for applications in parametric
integration and extended to weak approximation of stochastic
differential equations in \cite{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$, $\E{F}$.
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,
\begin{equation}\label{eq:MLMC-telescopic}
  \E{F} \approx \E{F_L} =
\E{F_0} + \sum_{\ell=1}^L\E{F_\ell - F_{\ell-1}},
\end{equation}
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 $\E{F_\ell-F_{\ell-1}}$ 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
\begin{equation}\label{eq:MLMC-estim}
  \E{F} \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,
\end{equation}
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. \cite{bsz11,bls13,scheichl.charrier:MLMC,scheichl.giles:MLMC,mss12}.
Recent works \cite{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 \cite{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
\cite{abdullatif.etal:MultiIndexMC}, that differs from the Multilevel
Monte Carlo method in that the telescoping idea presented in equations
\eqref{eq:MLMC-telescopic}-\eqref{eq:MLMC-estim} 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
\cite{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
\cite{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
\cite{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
\cite{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 \cite{griebel.knapek:optimized}).
A somewhat analogous approach was proposed in \cite{bieri:sparse.tensor.coll}, where
the number of solves per discretization level is prescribed \emph{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
\ref{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 \ref{s:method}, and our
main theorem detailing the complexity of MISC for a particular choice
of an index set is presented in Section \ref{s:complexity}.  Finally,
Section \ref{s:numerics} presents some numerical tests, while Section
\ref{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:
\begin{itemize}
  \item $\nset$ denotes the set of integer numbers including zero;
  \item $\nset_+$ denotes the set of positive integer numbers, i.e. excluding zero;
  \item $\rset_+$ denotes the set of positive real numbers, $\rset_+ = \{ r \in \rset : r > 0\}$;
  \item $\bm{1}$ denotes a vector whose components are always equal to one;
  \item $\vec e^\kappa_\ell$ denotes the $\ell$-th canonical vector in $\rset^\kappa$,
    i.e., $(\vec 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 $\vec v \in \rset^N$, we will write
    $\vec v - \vec e_1$ instead of $\vec v - \vec e_1^N$;
  \item given $\vec v \in \rset^N$, $|\vec v| = \sum_{n=1}^N v_n$,
    $\max(\vec v) = \max_{n=1,\ldots N}v_n$ and $\min(\vec v) = \min_{n=1,\ldots N}v_n$;
  \item given $\vec v \in \rset^N$ and $f:\rset \to \rset$,
    $f(\vec v)$ denotes the vector obtained by applying $f$ to each component of $\vec v$,
    $f(\vec v) = [f(v_1),f(v_2),\cdots,f(v_N)] \in \rset^N$;
  \item given $\vec v, \vec w \in \rset^N$, the inequality $\vec v > \vec w$ holds true if and only if $v_n > w_n$ $\forall n=1,\ldots,N$.
  \item given $\vec v \in \rset^D$ and $\vec w \in \rset^N$, $[\vec v, \vec w] = (v_1,\ldots,v_D,w_1,\ldots, w_N) \in \rset^{D+N}$.
\end{itemize}

 %
\section{Problem setting}\label{s:problem-setting}

Let $\mathscr{B} \subset \rset^{d}$, $d=1,2,3$, be an open
hyper-rectangular domain (referred to hereafter as the ``physical domain'')
and let $\vec 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 \rset$ 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(\vec y)d\vec 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(\vec y)d\vec y)$ is a complete probability space.
Consider the following generic PDE, together with the assumption stated next:
\begin{problem}\label{pb:strong_form_yy}
Find $u:\mathscr{B} \times \Gamma \rightarrow \rset$ such that for
$\rho$-almost every $\vec y \in \Gamma$
\[%
  \begin{cases}
    \mathcal{L}(u;\vec x,\vec y) = \mathscr{F}(\vec x)  & \vec x\in \mathscr{B},\\
    u(\vec x,\vec y)=0 %
 & \vec x\in \partial \mathscr{B}.
  \end{cases}
\]
\end{problem}
\begin{assumption}[Well posedness]\label{assump:wellposed}
  Problem \ref{pb:strong_form_yy} is well posed in some Hilbert
  space $V$ for $\rho$-almost every $\vec y \in \Gamma$.
\end{assumption}
The solution of Problem \ref{pb:strong_form_yy} can be
seen as an $N$-variate Hilbert-space valued function $u(\vec 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\`eve or Fourier expansion, possibly truncated after a finite number of terms,
see, e.g., \cite{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(\vec y) \|_{V}^2 \rho(\vec y)d\vec y < \infty\right\}.
$
Finally, given some functional of the solution $u$, $\Theta: V \rightarrow \rset$,
we denote by $F : \Gamma \rightarrow \rset$
the $N$-variate real-valued function assigning to each realization $\vec y \in \Gamma$ the corresponding
value of $\Theta[u]$ (quantity of interest), i.e., $F(\vec y) = \Theta[u(\cdot,\vec y)]$,
and we aim at estimating its expected value,
\[\E{F} = \int_\Gamma F(\vec y) \rho(\vec y) d\vec y.\]
\begin{example}\label{example:ell}
As a motivating example, consider the following elliptic problem:
find $u:\mathscr{B} \times \Gamma \rightarrow \rset$
such that for $\rho$-almost every $\vec y \in \Gamma$
\begin{equation}\label{eq:ell_PDE_strong_yy_example}
  \begin{cases}
    -\div( a(\vec x,\vec y)\nabla u(\vec x,\vec y) ) =\mathscr{F}(\vec x)  & \vec x\in \mathscr{B},\\
    u(\vec x,\vec y)=h(\vec x) & \vec x\in\partial \mathscr{B},
  \end{cases}
\end{equation}
holds, where $\div$ and $\nabla$ denote differentiation with respect to the physical
variables, $\vec x$, only, and the function $a: \mathscr{B} \times \Gamma \rightarrow \rset$
is bounded away from $0$ and $\infty$, i.e., there exist two constants, $a_{min}, a_{max}$, such that
\begin{equation}\label{eq:bounding_a}
0< a_{min} \leq a(\vec x,\vec y) \leq a_{max} <\infty,
\quad \forall \vec x \in \mathscr{B} \mbox{ and for } \rho\mbox{-almost every } \vec y \in \Gamma.
\end{equation}
This boundedness condition guarantees that Assumption \ref{assump:wellposed} is satisfied,
i.e. the equation is well posed for $\rho$-almost every $\vec 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., \cite{babuska.nobile.eal:stochastic2}.
This is the example we will focus on in Section \ref{s:numerics}, where we will test
numerically the performance of the Multi-Index Stochastic Collocation method that we will
detail in Section \ref{s:method}.
\end{example}

\begin{remark}
  The method that we present in the following sections can be also applied
  to more general problems than Problem \ref{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.
\end{remark}
\begin{remark}
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 \cite{hughes:IGA}
or Transfinite Interpolation \cite{gordon.hall:transfinite.interp})
or by a Domain Decomposition approach \cite{quarteroni1999domain}
if the domain can be obtained as a union of hyper-rectangles.
\end{remark}

\subsection{Approximation along the deterministic and stochastic dimensions} \label{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 \ref{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 $\valpha \in \nset^D_+$, we denote by
$u^\valpha(\vec x,\vec y)$ the approximation of $u$ obtained from setting $h_i = h_{i,\alpha_i}$,
with the implicit assumption that $u^\valpha(\vec x,\vec y) \to u(\vec x,\vec y)$
as $\min_{1\leq i \leq D} \alpha_i \to \infty$ for $\rho$-almost every $\vec y \in \Gamma$;
similarly, we also write $F^\valpha(\vec y) = \Theta[u^\valpha(\cdot,\vec y)]$.
For instance, we could solve the problem stated in Example \ref{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^\valpha$ 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^\valpha$ is at least a continuous function over $\Gamma$
(the smoother $F^\valpha$ is, the more effective the Lagrangian approximation will be);
for instance, for the problem stated in Example \ref{example:ell}, it can be shown under
moderate assumptions on $a(\vec x,\vec y)$ that $F$ and $F^\valpha$ are $\vec y$-analytic, see, e.g.,
\cite{back.nobile.eal:optimal,cohen.devore.schwab:nterm2}; we will
return to this point in Section \ref{s:numerics}.

To derive a generic tensor Lagrangian interpolation of $F^\valpha$, we first introduce the set
$\cset^0(\Gamma_n)$ of real-valued continuous functions over $\Gamma_n$,
and the subspace of polynomials of degree at most $q$ over $\Gamma_n$,
$\pset^q(\Gamma_n) \subset \cset^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 \nset_+}$, 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 : \cset^0(\Gamma_n) \rightarrow \pset^{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 \cset^0(\Gamma_n)$.

The $N$-variate Lagrangian interpolant can then be built by a tensorization
of univariate interpolants:
denote by $\cset^0(\Gamma)$ the space of real-valued $N$-variate continuous functions over $\Gamma$
and by $\pset^\vec q(\Gamma) = \bigotimes_{n=1}^N \pset^{q_n}(\Gamma_n)$ the subspace of
polynomials of degree at most $q_n$ over $\Gamma_n$, with $\vec q = (q_1,\ldots,q_N) \in \nset^N$,
and consider a multi-index $\vbeta \in \nset_+^N$ assigning the interpolation level in
each direction, $y_n$; the multivariate interpolant can then be written as
\begin{align*}
\mathscr{U}^{m(\vbeta)} : \cset^0(\Gamma) \rightarrow \pset^{m(\vbeta) -\bm{1}}(\Gamma), \qquad
\mathscr{U}^{m(\vbeta)}[F^\valpha] (\vec y)
&= \left( \mathscr{U}^{m(\beta_1)}_1 \otimes \cdots \otimes \mathscr{U}^{m(\beta_N)}_N \right) [F^\valpha](\vec y).%
\end{align*}

The set of collocation points needed to build the
tensor interpolant $\mathscr{U}^{m(\vbeta)}[u] (\vec y)$
is the tensor grid $\mathscr{T}^{m(\vbeta)} = \times_{n=1}^N \mathscr{H}_n^{m(\beta_n)}$
with cardinality $ \#\mathscr{T}^{m(\vbeta)} = \prod_{n=1}^N m(\beta_n)$.
Observe that the Lagrangian interpolant immediately induces an $N$-variate quadrature formula,
$\mathscr{Q}^{m(\vbeta)} : \cset^0(\Gamma) \rightarrow \rset$,
\begin{align*}
& \mathscr{Q}^{m(\vbeta)}[F^\valpha]
= \E{{\mathscr{U}}^{m(\vbeta)}[F^\valpha](\vec y)}
= \sum_{j=1}^{\# \mathscr{T}^{m(\vbeta)} } F^\valpha(\widehat{\vec y}_j) \varpi_j,
\end{align*}
where $\widehat{\vec y}_j \in \mathscr{T}^{m(\vbeta)}$ and the quadrature weights $\varpi_j$
are the expected values of the Lagrangian polynomials centered
in $\widehat{\vec 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(\vbeta)}$ and $\mathscr{Q}^{m(\vbeta)}$.
Common choices are Gaussian quadrature points
like Gauss-Legendre for uniform measures or Gauss-Hermite for Gaussian measures,
cf. e.g., \cite{trefethen:book2013}, which are however
\emph{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 \emph{nested} collocation points, and specifically with
Clenshaw-Curtis points \cite{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 \cite{nobile.eal:optimal-sparse-grids,Chkifa:leja},
whose performance is somehow equivalent to that of Clenshaw-Curtis
for quadrature purposes, see \cite{nobile.etal:leja,narayan:Leja}.
Clenshaw-Curtis points are defined as
\begin{equation}\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),
\end{equation}
together with the following level-to-nodes relation, $m(i_n)$, that ensures their nestedness:
\begin{equation}\label{eq:m_CC}
m(0)=0,\,\, m(1)=1,\,\, m(i_n)=2^{i_n-1}+1.
\end{equation}
We conclude this section by introducing the following operator norm,
which acts as a ``Lebesgue constant''
from $\cset^0(\Gamma)$ to $L^2_\rho(\Gamma)$:
\begin{equation}\label{eq:leb_def}
\mathbb{M}^{m(\vbeta)} = \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)}.
\end{equation}
In particular, for the Clenshaw-Curtis points, it is possible to bound $\mathbb{M}_n^{m(\beta_n)}$ as:
\begin{equation}\label{eq:CC_leb_est}
  \mathbb{M}^{m(\vbeta)} \leq \mathbb{M}_{est}^{m(\vbeta)} = \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}
\end{equation}
See \cite{nobile.eal:optimal-sparse-grids} and references therein.

\begin{remark}
  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
  \cite{genz.keister:kpnquad,back.nobile.eal:lognormal}; we mention
  also the recent work \cite{narayan:Leja} on generalized Leja points
  that can be used for arbitrary measures on unbounded domains.
\end{remark}

\section{Multi-Index Stochastic Collocation}\label{s:method}

It is easy to see that an accurate approximation of $\E{F}$ by a
direct tensor technique as the one just introduced,
$\E{F} \approx \mathscr{Q}^{m(\vbeta)} [F^\valpha]$, 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
\cite{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
\eqref{eq:MLMC-telescopic} and \eqref{eq:MLMC-estim}.  Denoting
$\mathscr{Q}^{m(\vbeta)} [F^\valpha] = \descqoi{\valpha}{\vbeta}$,
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{align}\label{eq:misc_delta_operations}
& \Delta_{i}^\textnormal{det}[\descqoi{\valpha}{\vbeta}] =
\begin{cases}
  \descqoi{\valpha}{\vbeta} - \descqoi{\valpha-\vec e_i}{\vbeta}, & \text{if }\alpha_i > 1,\\
  \descqoi{\valpha}{\vbeta} & \text{if } \alpha_i=1,
  \end{cases} \\% \label{eq:delta_det_def} \\
& \Delta_{j}^\textnormal{stoc}[\descqoi{\valpha}{\vbeta}]    =
\begin{cases}
  \descqoi{\valpha}{\vbeta} - \descqoi{\valpha}{\vbeta-\vec e_j}, & \text{if }\beta_j > 1,\\
  \descqoi{\valpha}{\vbeta} & \text{if } \beta_j=1.
  \end{cases} %
\end{align}
We then define the first-order tensor difference operators,
\begin{align}
& \vec \Delta^\textnormal{det}[\descqoi{\valpha}{\vbeta}]
= \bigotimes_{i=1}^D \Delta_i^\textnormal{det}[\descqoi{\valpha}{\vbeta}]
= \Delta_1^\textnormal{det} \left[ \, \Delta_2^\textnormal{det} \left[ \, \cdots
  \Delta_D^\textnormal{det} \left[ \descqoi{\valpha}{\vbeta} \right] \, \right] \, \right] %
= \sum_{\vec j \in \{0,1\}^D} (-1)^{|\vec j|} \descqoi{\valpha-\vec j}{\vbeta},
\label{eq:delta_det_def} \\
& \vec \Delta^\textnormal{stoc}[\descqoi{\valpha}{\vbeta}]
= \bigotimes_{j=1}^N \Delta_j^\textnormal{stoc}[\descqoi{\valpha}{\vbeta}]
= \sum_{\vec j \in \{0,1\}^N} (-1)^{|\vec j|} \descqoi{\valpha}{\vbeta-\vec j} \, .
\label{eq:delta_stoc_def}
\end{align}
with the convention that $F_{\vec v, \vec w}=0$ whenever a component
of $\vec v$ or $\vec w$ is zero. Observe that computing
$\vec \Delta^\textnormal{det}[\descqoi{\valpha}{\vbeta}]$ actually
requires up to $2^D$ solver calls, and analogously applying
$\vec \Delta^\textnormal{stoc}[\descqoi{\valpha}{\vbeta}]$ requires
interpolating $F^\valpha$ on up to $2^N$ tensor grids; for
instance, if $D=N=2$ and $\valpha,\vbeta > \bm{1}$, we have
\begin{align*}
  &\vec \Delta^\textnormal{det} [\descqoi{\valpha}{\vbeta}]
  = \Delta^\textnormal{det}_2 \left[ \, \Delta^\textnormal{det}_1 \left[ \, \descqoi{\valpha}{\vbeta} \, \right] \, \right]
   = \Delta^\textnormal{det}_2 [ \descqoi{\valpha}{\vbeta} - \descqoi{\valpha-\vec e_1}{\vbeta}]
   = \descqoi{\valpha}{\vbeta} - \descqoi{\valpha-\vec e_1}{\vbeta}
   - \descqoi{\valpha-\vec e_2}{\vbeta} + \descqoi{\valpha-\bm{1}}{\vbeta}, \\[2pt]%
  &\vec \Delta^\textnormal{stoc} [\descqoi{\valpha}{\vbeta}]
  = \descqoi{\valpha}{\vbeta} - \descqoi{\valpha}{\vbeta-\vec e_1}
   - \descqoi{\valpha}{\vbeta-\vec e_2} + \descqoi{\valpha}{\vbeta-\bm{1}}. \nonumber
\end{align*}
Finally, letting
$\vec \Delta[\descqoi{\valpha}{\vbeta}] =
\vec \Delta^\textnormal{stoc}[ \vec \Delta^\textnormal{det}[\descqoi{\valpha}{\vbeta}]]$,
we define the Multi-Index Stochastic Collocation (MISC) estimator of $\E{F}$ as
\begin{equation}\label{eq:misc_estimator}
  %
  \mathscr{M}_{\mathcal{I}}[F]
  = \sum_{[\valpha, \vbeta] \in \mathcal I} \vec \Delta[\descqoi{\valpha}{\vbeta}]
  = \sum_{ [\valpha, \vbeta] \in \mathcal I}  c_{\valpha,\vbeta} \descqoi{\valpha}{\vbeta},
\end{equation}
where $\mathcal{I} \subset \nset_+^{D+N}$ and
$c_{\valpha,\vbeta} \in \zset$. %
Observe that many of the coefficients in \eqref{eq:misc_estimator}, $c_{\valpha,\vbeta}$,
may be zero: in particular, $c_{\valpha,\vbeta}$ is zero whenever
$[\valpha,\vbeta]+\vec{j} \in \mathcal{I}$
$\forall \vec{j}\in\{0,1\}^{N+D}$.  Similarly to the analogous sparse grid construction
\cite{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{align*}
  \forall \,[\valpha, \vbeta] \in \mathcal{I}, \quad
  \begin{cases}
    \valpha - \vec{e}_i \in \mathcal{I} \mbox{ for } 1 \leq i \leq D \mbox{ and } \alpha_i > 1,\\
    \vbeta - \vec{e}_j \in \mathcal{I} \mbox{ for } 1 \leq j \leq N \mbox{ and } \beta_j > 1.
  \end{cases}
\end{align*}

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

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

The efficiency of the MISC method in equation \eqref{eq:misc_estimator}
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 \ref{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
\cite{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_{\valpha, \vbeta}$, and ``error contribution'',
$\Delta E_{\valpha, \vbeta}$, for each hierarchical surplus operator,
$\vec \Delta[\descqoi\valpha\vbeta]$. %
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 $\vec \Delta[\descqoi\valpha\vbeta]$ 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 \eqref{eq:delta_det_def} and \eqref{eq:delta_stoc_def};
the error contribution measures instead how much the error $|\E{F} - \mathscr{M}_\mathcal{I}[F]| $
would decrease once the operator $\vec \Delta[\descqoi\valpha\vbeta]$
has been added to $\mathscr{M}_\mathcal{I}[F]$. %
In formulas, we define
\[\Delta W_{\valpha, \vbeta}
= \work{\mathscr{M}_{\mathcal{I} \cup \{[\valpha,\vbeta]\}}[F]} - \work{\mathscr{M}_{\mathcal{I}}[F]}
= \work{\vec \Delta[\descqoi\valpha\vbeta]},\]
so that
\begin{equation}\label{eq:work_decomp}
  \work{\mathscr{M}_\mathcal{I}[F]} = \sum_{[\valpha, \vbeta] \in \mathcal{I}} \Delta W_{\valpha, \vbeta},
\end{equation}
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 $(\valpha, \vbeta)$ 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 \eqref{eq:misc_estimator}, $c_{\valpha,\vbeta}$, 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 \{[\valpha,\vbeta]\}$, hence it would be possible not to compute the
corresponding approximations $\descqoi{\valpha}{\vbeta}$.
This approach is discussed in Section~\ref{ss:impl}.

Similarly, we define
\[\Delta E_{\valpha, \vbeta}
 = \Big| {\mathscr{M}_\mathcal{I \cup \{[\valpha,\vbeta]\}}[F]} - {\mathscr{M}_\mathcal{I}[F]} \Big|
 = \left| {\vec \Delta[\descqoi\valpha\vbeta]} \right|.
\]
Thus, by construction, the error of the MISC estimator \eqref{eq:misc_estimator}
can be bounded as the sum of the error contributions not included in
the estimator $\mathscr{M}_\mathcal{I}[F]$, %
\begin{align}\label{eq:error-decomp}
\error{\mathscr{M}_\mathcal{I}[F]} =
|\E{F} - \mathscr{M}_\mathcal{I}[F]|
& = \left| {\sum_{[\valpha, \vbeta] \notin \mathcal I} \vec \Delta[\descqoi\valpha\vbeta]} \right| \nonumber \\[2mm]
 & \leq \sum_{[\valpha, \vbeta] \notin \mathcal I} \left| { \vec \Delta[\descqoi\valpha\vbeta] } \right|
  = \sum_{[\valpha, \vbeta] \notin \mathcal I} \Delta E_{\valpha, \vbeta}.
\end{align}
Consequently, a quasi-optimal set $\mathcal{I}$ can be computed by solving the following
``binary knapsack problem'' \cite{martello:knapsack}:
\begin{align}
  \mbox{maximize} &\sum_{[\valpha, \vbeta] \in \nset_+^{D+N}}\Delta E_{\valpha, \vbeta} x_{\valpha, \vbeta} \nonumber \\
  \mbox{such that }& \sum_{[\valpha, \vbeta] \in \nset_+^{D+N}}\Delta W_{\valpha, \vbeta} x_{\valpha, \vbeta} \leq W_{max}, \label{eq:knapsack} \\
  & \qquad x_{\valpha, \vbeta} \in \{0,1\}, \nonumber
\end{align}
and setting
$\mathcal{I}=\{ [\valpha, \vbeta] \in \nset_+^{D+N}: x_{\valpha,
  \vbeta} = 1\}$.
Observe that such a set is only ``quasi'' optimal since the error
decomposition \eqref{eq:error-decomp} 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_{\valpha, \vbeta}$ 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 \cite{martello:knapsack}:
\begin{enumerate}
\item compute the ``profit'' of each hierarchical surplus, i.e., the quantity
  \[%
    P_{\valpha, \vbeta} = \frac{\Delta E_{\valpha, \vbeta}}{\Delta W_{\valpha, \vbeta}};
  \]
\item sort the hierarchical surpluses by decreasing profit;
\item add the hierarchical surpluses to $\mathscr{M}_\mathcal{I}[F]$ according to such order
  until the constraint on the maximum work is fulfilled.
\end{enumerate}
Note that by construction $x_{\valpha, \vbeta}=1$ for all the multi-indices included
in the selection except for the last one, for which $x_{\valpha, \vbeta}<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_{\valpha, \vbeta}$ have integer values (see also \cite{b.griebel:acta});
observe that this integer
solution is also the solution of the original binary knapsack problem \eqref{eq:knapsack}
with the new value of $W_{max}$ in the work constraint.
Thus, if the quantities $\Delta E_{\valpha,\vbeta}$ and $\Delta W_{\valpha,\vbeta}$ were available,
the quasi-optimal index set for the MISC estimator could be computed as
\begin{equation}\label{eq:opt_set}
 \mathcal{I} = \mathcal{I}(\epsilon) =\left\{ [\valpha, \vbeta] \in \nset_+^{D+N}\::\:
   \frac{\Delta E_{\valpha,\vbeta}}{\Delta W_{\valpha,\vbeta}} \geq \epsilon \right\},
\end{equation}
for a suitable $\epsilon>0$.%
                            %

\begin{remark}\label{remark:webster_and_bieri}
  The MISC setting could in principle include the Multilevel Stochastic Collocation method
  proposed in \cite{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 \cite{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
  \cite{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
  \cite{bieri:sparse.tensor.coll}, by considering again only one
  deterministic discretization parameter, i.e., $D=1$, and then
  setting
  \begin{equation}\label{eq:iso-MISC}
    \mathcal{I} = \left\{ [\alpha,\vbeta] \in \nset_+^{1+N} : \alpha + \sum_{n=1}^N \beta_n \leq w \right\},
  \end{equation}
  with $w \in \nset_+$.  In other words, the approach in \cite{bieri:sparse.tensor.coll}
  is based neither on profit nor on error balancing.
\end{remark}

 %
\section{Complexity analysis of the MISC method}\label{s:complexity}

In this section, we assume suitable models for the error and work
contributions, $\Delta E_{\valpha,\vbeta}$ and
$\Delta W_{\valpha,\vbeta}$ (which are numerically verified in
Section \ref{s:numerics} for the problem in Example \ref{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.

\begin{assumption}\label{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}.\]
\end{assumption}
\begin{assumption}\label{assump:dW_dE_factor}
  The error and work contributions, $\Delta E_{\valpha,\vbeta}$ and $\Delta W_{\valpha,\vbeta}$,
  can be bounded as products of two terms,
  \[%
    \Delta E_{\valpha,\vbeta} \leq \Delta E_\valpha^\textnormal{det}  \Delta E_\vbeta^\textnormal{stoc}  \qquad \text{and}\qquad
    \Delta W_{\valpha,\vbeta} \leq \Delta W_\valpha^\textnormal{det}  \Delta
    W_\vbeta^\textnormal{stoc},  %
\]
  where $\Delta W_\valpha^\textnormal{det}$ and $\Delta E_\valpha^\textnormal{det}$ denote the cost and
  the error contributions due to the deterministic difference operator,
  $\Delta^\textnormal{det}[\descqoi{\valpha}{\vbeta}]$,
  and similarly $\Delta W_\vbeta^\textnormal{stoc}$ and $\Delta E_\vbeta^\textnormal{stoc}$ denote the
  cost and the error contribution
  due to the stochastic difference operator, $\Delta^\textnormal{stoc}[\descqoi{\valpha}{\vbeta}]$,
  cf. equations \eqref{eq:delta_det_def}-\eqref{eq:delta_stoc_def}.
\end{assumption}
\begin{assumption}\label{assump:dW_dE_model}
  The following bounds hold true for the factors appearing in Assumption \ref{assump:dW_dE_factor}:
  \begin{align}
    & \Delta W_\valpha^\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_\valpha^\textnormal{det} \leq C_{\mathrm{error}}^{\textnormal{det}} \prod_{i=1}^D (h_{i,\alpha_i})^{\widetilde{\dRate{}_i} } , \label{eq:dE_alpha_est}\\
    & \Delta W_\vbeta^\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_\vbeta^\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{align}
  for some rates $\widetilde{\gamma}_i, \widetilde{\dRate{}}_i, \widetilde{g}_i > 0$.
\end{assumption}

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
\cite{griebel:infdim, griebel:tensor} for similar arguments for
different choices of finite and infinite dimensional index sets.
\begin{restatable}[MISC computational complexity]{theorem}{misccomplexity}
\label{thm:misc_complexity}
Under Assumptions \ref{assump:growth_of_dof} to \ref{assump:dW_dE_model}, the bounds
for the factors appearing in Assumption \ref{assump:dW_dE_factor} can be equivalently rewritten as
  \begin{subequations}
  \begin{align}
    \label{eq:deltaW_ass}
    & \Delta W_{\valpha,\vbeta} \leq \CWork e^{\sum_{i=1}^D \gamma_i \alpha_i} e^{ \delta |\vbeta|}, \\
    \label{eq:deltaE_ass}
    & \Delta E_{\valpha,\vbeta} \leq \CError e^{ - \sum_{i=1}^D \dRate{i}\alpha_i} e^{- \sum_{j=1}^N g_j \exp(\delta {\beta_j})},
  \end{align}
  \end{subequations}
  with $\gamma_i = \widetilde{\gamma}_i\log 2$, $\dRate{i} =\widetilde{\dRate{}}_i\log 2$,
  $\delta = \log 2$ and $g_i = \widetilde{g}_i C_{m,low}$. Define the following set
  \begin{equation}\label{eq:opt_set_detailed_logform}
    \mathcal{I}^*(L) =
    \left\{ [\valpha, \vbeta] \in \nset_+^{D+N}\::\:
      \sum_{i=1}^D (\dRate{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 \rset_+.
  \end{equation}
  Then there exists a constant $\mathscr C_{\text{W}}$ such that, for any $W_{\max}$ satisfying
  \begin{equation}
    \label{eq:maxW_ass}
    W_{\max} \geq \mathscr C_{\text{W}} \exp\left(\chi\right),
      %
  \end{equation}
  and choosing $L$ as
  \begin{equation}\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),
  \end{equation}
  with $\vXi = \left( \frac{\gamma_1}{\gamma_1 + r_1},\ldots,\frac{\gamma_D}{\gamma_D + r_D} \right)$,
  $\chi = \max(\vXi)$,  $\zeta=\min_{i=1,\ldots,D} \frac{\dRate{i}}{\gamma_i}$ and
  $\mathfrak{z} = \#\{i=1,\ldots D \::\: \frac{\dRate{i}}{\gamma_i} =
  \zeta\}$, the MISC estimator $\mathscr M_{\mathcal I^*(L(W_{\max}))}$ satisfies
  \begin{subequations}
    \label{eq:misc_complexity}
    \begin{align}
      &\work{\mathscr M_{\mathcal I^*(L(W_{\max}))}} \leq W_{\max}, \label{eq:misc_complexity_work}
      \\[2mm]
      &\limsup_{W_{\max} \uparrow \infty} \frac
{\error{\mathscr M_{\mathcal I^*(L(W_{\max}))}}}{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{align}
  \end{subequations}
\end{restatable}

\begin{remark}\label{rem:this-is-the-set}
  The set $\mathcal{I}^*$ proposed in Theorem \ref{thm:misc_complexity} can be obtained
  by assuming that the bounds in equations \eqref{eq:deltaW_ass} and \eqref{eq:deltaE_ass}
  are actually equalities and by using the definition of the quasi-optimal set \eqref{eq:iso-MISC}:
  \begin{align*}
    \mathcal{I}^*
    & = \left\{ [\valpha,\vbeta] \in \nset_+^{D+N} : \frac{\Delta E_{\valpha,\vbeta}}{\Delta W_{\valpha,\vbeta}} \geq \epsilon \right \} \\
    & = \left\{ [\valpha,\vbeta] \in \nset_+^{D+N} :
      \frac{e^{ - \sum_{i=1}^D \dRate{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 |\vbeta|}} \geq \epsilon \right \} \\
    & = \left\{ [\valpha, \vbeta] \in \nset_+^{D+N}\::\:
      \sum_{i=1}^D (\dRate{i}+\gamma_i) \alpha_i + \sum_{i=1}^N (\delta \beta_i +  g_i e^{\delta \beta_i} ) \leq L \right\},
  \end{align*}
  where the last equality holds with $L = - \log \epsilon$.
\end{remark}

\begin{remark}\label{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
  \eqref{eq:deltaE_ass}, 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 \ref{thm:misc_complexity}
  only depends on $\gamma_i$ and $\dRate{i}$, i.e.,
  the cost and error rates of the deterministic solver, respectively.
  This observation coincides with that in \cite[page 2299]{bieri:sparse.tensor.coll}:
  ``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 \cite{bieri:sparse.tensor.coll},
  MISC takes greater advantage of this fact, since it is based
  on an optimization procedure, cf. equation \eqref{eq:opt_set};
  this performance improvement is well documented by the comparison
  between the two methods shown in the next section.
  %
  %
  %
  %
  Figure \ref{fig:MISC_set} shows the multi-indices included in
  $\mathcal{I}$ according to \eqref{eq:opt_set_detailed_logform}
  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.
\begin{figure}[h]
  \centering
  \includegraphics[width=0.5\textwidth]{opt_MISC_set.pdf}
  \caption{Index sets $\mathcal{I}(L)$ for $D=N=1$, according to
    equation \eqref{eq:opt_set_detailed_logform}.}
  \label{fig:MISC_set}
\end{figure}
\end{remark}

\begin{remark}
  Theorem~\ref{thm:misc_complexity} is only valid in case
  Assumptions~\ref{assump:wellposed}--~\ref{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.
\end{remark}

 %
\section{Example and numerical evidence}\label{s:numerics}

In this section, we test the effectiveness of the MISC approximation on some instances
of the general elliptic equation \eqref{eq:ell_PDE_strong_yy_example} in Example \ref{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}(\vec x)=1$.
As for the random diffusion coefficient, we set %
\begin{equation}\label{eq:log-uniform-a}
      a(\vec x,\vec y)= e^{\gamma_N(\vec x,\vec y)}, \quad
      \gamma_N(\vec x,\vec y) = \sum_{n=1}^N \lambda_n \psi_n(\vec x) y_n,
\end{equation}
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(\vec 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~\ref{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 \eqref{eq:bounding_a}, $a_{min}$ and $a_{max}$,
that in turn assures the well posedness of the problem.
Finally, the quantity of interest is defined as
\begin{equation}\label{eq:qoi_def}
  F(\vec y) = \int_\mathscr{B} u(\vec x,\vec y) Q(\vec x) d\vec x,
  \quad Q(\vec x) = \frac{1}{(\sigma\sqrt{2\pi})^d}\exp \left( -\frac{\|\vec x-\vec x_0\|_2^2}{2\sigma^2} \right)
\end{equation}
with $\sigma=0.16$ and locations $\vec x_0=0.3$ for $d=1$ and
$\vec x_0=[0.3, 0.2, 0.6]$ for $d=3$. We also make the choice $h_0=1/3$
in Assumption~\ref{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.

\begin{table}[h]
\centering
\begin{tabular}{l|rrrrrrrrrr}
$n$ & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\
\hline
$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 \\
\hline
\end{tabular}
\caption{Included functions for $d=3$ in \eqref{eq:log-uniform-a}. Here
  $\psi_n(\vec x) = \phi_{i(n)}(x_1)\phi_{j(n)}(x_2)\phi_{k(n)}(x_3)$.}\label{tab:modes}
\end{table}

\subsection{Verifying bounds on work and error contributions}
In this subsection we discuss the validity of Assumptions \ref{assump:growth_of_dof}
to \ref{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.

\paragraph{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 \ref{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:
\[\work{F_\valpha} \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 \eqref{eq:deltaW_alpha_est}
with the following argument: since computing
$\vec \Delta^\textnormal{det}[F^\valpha]$ requires up to $2^D$
solver calls, each on a different grid (cf. equation
\eqref{eq:delta_det_def}), we have
\begin{align*}
\Delta W_\valpha^\textnormal{det}
& = \work{\vec \Delta^\textnormal{det}[ F_\valpha ]}
= \sum_{\vec j \in \{0,1\}^D} \work{F_{\valpha-\vec j}} \\
& \leq C_{\mathrm{GMRES}}\sum_{\vec 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_{\vec 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{align*}
i.e., bound \eqref{eq:deltaW_alpha_est} 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_\valpha^\textnormal{det}$, we observe numerically
that bound \eqref{eq:dE_alpha_est} holds true in practice with $\widetilde{\dRate{}}_i=2$,
$i=1,\ldots,D$, due to the fact that $a \in C^{\infty}(\mathscr{B})$
for $\rho$-almost every $\vec y \in \Gamma$,
$\mathscr{F} \in C^{\infty}(\mathscr{B})$
and the function $Q$ appearing in the quantity of interest \eqref{eq:qoi_def}
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 \ref{assump:dW_dE_factor} is valid
(we will numerically verify it later in this section),
in Figure \ref{fig:check_rates_det} we show the value of
$\Delta E_{\valpha,\vbeta} = \Delta E_\valpha^\textnormal{det} \Delta E_\vbeta^\textnormal{stoc}$ for fixed $\vbeta=\bm{1}$
and variable $\valpha = j \bar{\valpha} + \bm{1}, j=1,2,\ldots$,
as well as the corresponding value of the bound \eqref{eq:dE_alpha_est} for $\Delta E_\valpha^\textnormal{det}$.
The line obtained by choosing $\bar{\valpha} = [1,\, 0,\, 0]$ confirms that the size of $\Delta E_\valpha^\textnormal{det}$
indeed decreases exponentially fast with respect to $\alpha_1$, and by fitting the
computed values of $\Delta E_\valpha^\textnormal{det}$
we obtain that the convergence rate is $\widetilde{\dRate{}}_j=2$, as previously mentioned;
analogous conclusions can be obtained for $\alpha_2$ and $\alpha_3$
by setting $\bar{\valpha} = [0,\, 1,\, 0]$ (shown in Figure \ref{fig:check_rates_det})
and $\bar{\valpha} = [0,\, 0,\, 1]$ (not shown). Most importantly,
confirmation of the product structure of $\Delta E_\valpha^\textnormal{det}$ can be obtained by observing, e.g.,
the decay of $\Delta E_\valpha^\textnormal{det}$ for $\bar{\valpha} = [1,\, 1,\, 0]$ and $\bar{\valpha} = [1,\, 1,\,1]$.

\begin{figure}
  \subfigure[For fixed ${\vbeta=\bm{1}}$ and variable
  ${\valpha = k\bar{\valpha}+\bm{1}}$.\label{fig:check_rates_det}]
  {\includegraphics[page=1,width=0.48\textwidth]{COOK_d3_N5_rates.pdf}}
  \subfigure[For fixed ${\valpha=\bm{1}}$ and variable
  ${\vbeta = k\bar{\vbeta}+\bm{1}}$.
  \label{fig:check_rates_stoc}]
  {\includegraphics[page=2,width=0.48\textwidth]{COOK_d3_N5_rates.pdf}}
  \caption{Verifying the validity of the bound \eqref{eq:deltaE_ass}
    for the value of $|\Delta E_{\valpha,\vbeta}|$ for the test case
    with $D=3$ and $N=5$.  The {\em dashed} lines are based on the model
    in \eqref{eq:deltaE_ass} with ${\tilde{\dRate{}}_i=2}$ for all
    $i=1,2,3$ and $g_j$ as in Table~\ref{tab:g_values} for
    $j=1,\ldots, 5$. The {\em solid} lines are based on computed
    values. }
\end{figure}

\paragraph{Stochastic discretization} The interpolation over the parameter space is based
on the tensorized Lagrangian interpolation technique with Clenshaw-Curtis
points explained in Section \ref{s:interpolation}, cf. eqs. \eqref{eq:CC_def} and \eqref{eq:m_CC}.
In particular, due to the nestedness of the Clenshaw-Curtis points,
adding the operator $\Delta^\textnormal{stoc}[\descqoi{\valpha}{\vbeta}]$ to the MISC estimator will require
$\Delta W_\vbeta^\textnormal{stoc} = \prod_{j=1}^N \big( m(\beta_j) - m(\beta_j-1) \big)$ new collocation points,
which in view of equation \eqref{eq:m_CC}  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 \ref{assump:growth_of_dof}
and bound \eqref{eq:deltaW_beta_est} in Assumption \ref{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_\vbeta^\textnormal{stoc}$ is not uniquely defined, i.e., it
depends on the set $\mathcal{I}$ to which  $\Delta^\textnormal{stoc}[\descqoi{\valpha}{\vbeta}]$ is being added,
see, e.g., \cite[Example 1 in Section 3]{nobile.eal:optimal-sparse-grids}.

Finally, to discuss the validity of bound \eqref{eq:deltaE_beta_est} for $\Delta E_\vbeta^\textnormal{stoc}$,
we rely on the theory developed in our previous works
\cite{back.nobile.eal:optimal,nobile.eal:optimal-sparse-grids}.
We begin by introducing the Chebyshev polynomials of the first kind
$\ortp_{q}(y)$ for $y \in [-1,1]$, which are defined by the relation
\[\ortp_{q}(\cos(\theta))
= \cos(q \theta), \quad 0 \leq \theta \leq \pi, \quad q \in \nset.\]
Then, for any multi-index $\vec q \in \nset^N$, we consider the
$N$-variate Chebyshev polynomials $\oortp_\vec q(\vec y) = \prod_{n=1}^N \ortp_{q_n}(y_n)$
and introduce the spectral expansion of $f:[-1,1]^N\rightarrow \rset$
over  $\{\oortp_\vec q\}_{\vec q \in \nset^N}$,
\[%
f(\vec y) = \sum_{\vec q \in \nset^N} f_\vec q \oortp_\vec q(\vec y),
\quad
f_\vec q = \int_\Gamma f(\vec y) \oortp_\vec q(\vec y) \prod_{n=1}^N
\frac{1}{\sqrt{1-y_n^2}} d\vec y,
\]
Next, given any $\ellrad_{1},\ellrad_{2},\ldots,\ellrad_{N}>1$
we introduce the Bernstein polyellipse
$\mathcal{E}_{\ellrad_1,\ldots,\ellrad_N} =\prod_{n=1}^N \mathcal{E}_{n,\ellrad_n}$,
where $\mathcal{E}_{n,\ellrad_n}$ denotes the ellipses in the complex plane defined as
\begin{align*}%
  \mathcal{E}_{n,\ellrad_n} = \bigg\{ z_n \in \mathbb{C} : \,\,
  \real{z} \leq \frac{\ellrad_{n} + \ellrad_{n}^{-1}}{2}\cos{\phi}, \,\,\,
  \im{z} \leq \frac{\ellrad_{n} - \ellrad_{n}^{-1}}{2}\sin{\phi},
  \,\, \phi \in [0, 2\pi) \bigg\},
\end{align*}
and recall the following lemma (see \cite[Lemma 2]{nobile.eal:optimal-sparse-grids} for a proof).

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

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

\begin{lemma}\label{lemma:qoi_polyellipse}
The quantity of interest $F(\vec y)=\Theta[u(\cdot,\vec y)]$ is analytic in a Bernstein polyellipse
with parameters $\ellrad_n = \tau_n + \sqrt{\tau_n^2+1}$, for any $\tau_n < \frac{\pi}{2 N \lambda_n}$.
\end{lemma}
\begin{proof}
Equation \eqref{eq:ell_PDE_strong_yy_example} can be extended in the complex domain
by replacing $\vec y$ with $\vec z \in \mathbb{C}^N$, and is analytic in the set
$\Sigma = \left\{ \vec z \in \mathbb{C}^N : \Re e{[a(\vec x,\vec z)]} >0 \right\}$,
see, e.g., \cite{babuska.nobile.eal:stochastic2}. By writing $z_n = b_n + i c_n$, we have
\begin{align*}
a(\vec x,\vec z)
& = \exp \left( \sum_n z_n \lambda_n \psi_n(\vec x) \right)
  = \exp \left( \sum_n b_n \lambda_n \psi_n(\vec x) \right) \exp \left( \sum_n i c_n \lambda_n \psi_n(\vec x) \right) \\
& = \exp \left( \sum_n b_n \lambda_n \psi_n(\vec x) \right)
	\left[ \cos \left(\sum_n c_n \lambda_n \psi_n(\vec x) \right) + i \sin \left( \sum_n c_n \lambda_n \psi_n(\vec x) \right) \right]
\end{align*}
so that the region $\Sigma$ can be rewritten as
\[
\Sigma =
\left\{ \vec z = \vec b + i \vec c \in \mathbb{C}^N : \cos \left( \sum_n c_n \lambda_n \psi_n(\vec x) \right) >0 , \forall \vec x \in \mathscr{B} \right\}.
\]
Such a region includes the smaller region
\[\Sigma_2 = \left\{ \vec z = \vec b + i \vec 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\{ \vec z = \vec b + i \vec 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 \eqref{eq:log-uniform-a}.
Next we let $\tau_n = \frac{\pi}{2 N \lambda_n}$
and define the following subregion of $\Sigma_3$:
\[\Sigma_4 = \bigg\{ \vec z = \vec b + i \vec 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 $\ellrad_n$ such that
\[\frac{\ellrad_n - \ellrad_n^{-1}}{2} = \tau_n
\Rightarrow \ellrad_n^2 - 1 - 2\tau_n\ellrad_n = 0
\Rightarrow \ellrad_n = \tau_n + \sqrt{\tau_n^2+1}
\]
in which $u(\vec x,\vec 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$.
\end{proof}
\begin{remark}
  Incidentally, we remark that the choice of $\tau_n$ considered in Lemma \ref{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$.
\end{remark}

\begin{lemma}
  \begin{equation}\label{eq:deltaE_beta_gen}
    \Delta E_\vbeta^\textnormal{stoc} \leq C_E e^{-\sum_{n=1}^N g^*_n m(\beta_n-1)} \mathbb{M}^{m(\vbeta)}
  \end{equation}
  holds, where $C_E= 4^N B \prod_{n=1}^N \frac{1}{1-e^{-g^*_n}}$, $B$ as in Lemma \ref{lemma:inc_cheb_conv},
  $g^*_n = \log \ellrad_n $ with $\ellrad_n$ as in Lemma \ref{lemma:qoi_polyellipse},
  and $\mathbb{M}^{m(\vbeta)}$ has been defined in equation \eqref{eq:leb_def}.
  %
  %
  %
  %
  %
\end{lemma}
\begin{proof}
  Combining Lemmas \ref{lemma:inc_cheb_conv} and \ref{lemma:qoi_polyellipse},
  we obtain that the Chebyshev coefficients of $F$ can be bounded as
  \[| F_\vec q | \leq C_\textnormal{Cheb}(\vec q) \prod_{n=1}^N e^{-g^*_n q_n},\]
  with $g^*_n = \log \ellrad_n = \log (\tau_n + \sqrt{\tau_n^2+1})$ and $\tau_n$
  as in Lemma \ref{lemma:qoi_polyellipse}. Then, the result can
  be obtained following the same argument of
  \cite[Lemma 5]{nobile.eal:optimal-sparse-grids}.
\end{proof}

To conclude, we first observe that $\mathbb{M}^{m(\vbeta)}$
grows logarithmically with respect to $m(\vbeta)$, see eq. \eqref{eq:CC_leb_est},
so it is asymptotically negligible in the estimate above, i.e. we can write
  \[\Delta E_\vbeta^\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 \eqref{eq:m_CC} implies that $m(i-1) \geq \frac{m(i)-1}{2}$.
We can finally write
\begin{align*}
  \Delta E_\vbeta^\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)}, %
\end{align*}
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 \eqref{eq:deltaE_beta_est} in Assumption \ref{assump:dW_dE_model} is valid for the test we are considering.
Finally, we point out that in practice we work with the expression \eqref{eq:deltaE_ass},
whose rates $g_n$ are actually better estimated numerically,
using the same procedure used to obtain the deterministic rates
$\tilde {\dRate{}}_j=2$:
we choose a sufficiently fine spatial resolution level $\valpha$, consider a variable $\vbeta = j\bar{\vbeta} + \bm{1}$
and fit the (simplified) model $\Delta E_\vbeta^\textnormal{stoc} \leq C \prod_{n=1}^N e^{-g_n 2^{\beta_n}}$.
The values obtained are reported in Table \ref{tab:g_values}, and they are found
to be equal for the case $d=1$ and $d=3$
(see also \cite{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 $\valpha$ 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 \ref{rem:only-space-matters}.
Figure \ref{fig:check_rates_stoc} then shows the validity of the bound
$\Delta E_\vbeta^\textnormal{stoc} \leq C \prod_{n=1}^N e^{-g_n 2^{\beta_n}}$
comparing for fixed $\valpha=\bm{1}$ and $\vbeta = j\bar{\vbeta} + \bm{1}$
the value of $\Delta E_{\valpha}^\textnormal{det} \Delta E_\vbeta^\textnormal{stoc}$ and the corresponding estimate.

\begin{table}[tbp]
  \centering
  \begin{tabular}{ccccccccccc}
    \hline
	& $g_1$	& $g_2$	& $g_3$	& $g_4$	& $g_5$	& $g_6$	& $g_7$	& $g_8$	& $g_9$	& $g_{10}$	\\
    \hline
	& 2.4855& 2.8174& 4.5044&4.1938 &4.7459 &6.8444	& 7.1513& 7.8622& 8.6584& 9.4545	\\
    \hline
  \end{tabular}
  \caption{Values of rates $g$ for the test cases considered.}
  \label{tab:g_values}
\end{table}

\paragraph{Stochastic-deterministic product structure}
We conclude this section by verifying Assumption \ref{assump:dW_dE_factor}, i.e., the fact that the error
contribution can be factorized as $\Delta E_{\valpha, \vbeta} = \Delta E_\valpha^\textnormal{det} \Delta E_\vbeta^\textnormal{stoc}$ and
that an analogous decomposition holds for $\Delta W_{\valpha,\vbeta}$. While the latter is trivial,
to verify the former we employ the same strategy used to verify the models
for $\Delta E_\valpha^\textnormal{det}$ and  $\Delta E_\vbeta^\textnormal{stoc}$, this time letting both $\valpha$
and $\vbeta$ change for every point,
i.e., $\valpha = j \bar{\valpha} + \bm{1}$ and  $\vbeta = j \vbeta_0 + \bm{1}$.
Figure \ref{fig:check_mixed_rates} shows the comparison between the computed value of
$\Delta E_{\valpha,\vbeta}$ and their estimated counterpart and confirms the validity
of the product structure assumption.

\begin{figure}
  \centering
  \includegraphics[page=3,width=0.5\textwidth]{COOK_d3_N5_rates.pdf}
  \caption{Comparison of $|\Delta E_{\valpha,\vbeta}|$ for
    $\vbeta=j\bar{\vbeta}+\bm{1}$ and $\valpha = j\bar{\valpha}+\bm{1}$
    for the test case with $D=3$ and $N=5$.
    The {\em dashed} lines are based on the model
    in \eqref{eq:deltaE_ass} with ${\tilde{\dRate{}}_i=2}$ for all
    $i=1,2,3$ and $g_j$ as in Table~\ref{tab:g_values} for
    $j=1,\ldots, 5$. The {\em solid} lines are based on computed
    values.}
  \label{fig:check_mixed_rates}
\end{figure}

\subsection{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 $\E{F}$
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 \eqref{eq:work_decomp} and \eqref{eq:deltaW_ass}. The names used here for the methods
are also used in the legends of the figures showing the convergence plots.

\begin{description}

\item[``a-priori'' MISC] refers to the MISC method with index
  set $\mathcal I$ defined by \eqref{eq:opt_set}, where
  $\Delta W_{\valpha, \vbeta}$ and $\Delta E_{\valpha, \vbeta}$ are
  taken to equal their upper bounds in \eqref{eq:deltaW_ass} and
  \eqref{eq:deltaE_ass}, respectively. The resulting set is explicitly
  written in~\eqref{eq:opt_set_detailed_logform}.
  The convergence rate of this set is predicted by Theorem \ref{thm:misc_complexity},
  cf. Remark \ref{rem:this-is-the-set}. Note that we do not need to determine
  the value of the constants $\CWork$ and $\CError$ since they can be
  absorbed in the parameter $\epsilon$ in \eqref{eq:opt_set}.
  %

\item[``a-posteriori'' MISC] refers to the MISC method with
  index set $\mathcal I$ defined by \eqref{eq:opt_set}, where
  $\Delta W_{\valpha, \vbeta}$ is taken to equal its upper bound in
  \eqref{eq:deltaW_ass}, and $\Delta E_{\valpha, \vbeta}$ is instead
  computed explicitly as $\left| \Delta[\descqoi{\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 \eqref{eq:deltaE_ass}.

\item[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 \ref{remark:webster_and_bieri};
  we recall this is not exactly the MLSC method that was implemented in
  \cite{teckentrup.etal:MLSC}, see again Remark \ref{remark:webster_and_bieri}.
  Just as with MISC, we consider both the ``a-priori'' and ``a-posteriori'' version of MLSC,
  where $\Delta E_{\valpha, \vbeta}$ is taken to be equal to its upper \eqref{eq:deltaE_ass}
  in the former case and assessed by direct numerical evaluation in the latter case.

\item[SCC] refers to the ``Sparse Composite Collocation method''
  in Remark \ref{remark:webster_and_bieri}, see equation \eqref{eq:iso-MISC}.

\item[MIMC] refers to the Multi-Index Monte Carlo method as detailed in
  \cite{abdullatif.etal:MultiIndexMC}, for which the complexity
  $\Order{W_{\max}^{-0.5}}$ can be estimated
  %
   for the test case at hand and as long as $d<4$.

\item[SGSC] refers to the quasi-optimal Sparse Grids Stochastic Collocation (SGSC)
  with fixed spatial discretization as proposed in
  \cite{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 \ref{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.
\end{description}

\begin{figure}[htp]
  \centering
  \includegraphics[page=1,width=0.7\textwidth]{COOK_d3_N10.pdf}
  \caption{Envelope of SGSC convergence curves for the test case with $d=3$ and $N=10$.}
  \label{fig:envelope}
\end{figure}

\subsection{Implementation details}
\label{ss:impl}
To implement MISC, we need two components:
\begin{enumerate}
\item Given a profit level parameter, $\epsilon$, we build the quasi-optimal set
  $\mathcal I$ based on~\eqref{eq:opt_set}, \eqref{eq:deltaW_ass}
  and~\eqref{eq:deltaE_ass}. One method to achieve this is
  to exploit the fact that this set is downward closed and use the
  following recursive algorithm.
  \begin{verbatim}
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
  \end{verbatim}
\item Given the set, $\mathcal I(L)$, we
  evaluate \eqref{eq:misc_estimator}. Here we have two choices:
\begin{itemize}
\item Evaluate the individual terms
  $\vec \Delta[\descqoi{\valpha}{\vbeta}]$ for every
  $\valpha,\vbeta \in \mathcal I$. To do so, we use the operator
  defined in~\eqref{eq:misc_delta_operations} 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 $\descqoi{\valpha}{\vbeta}$ 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.
\item Use the combination form of~\eqref{eq:misc_estimator} and
  only compute the terms that have $ c_{\valpha,\vbeta} \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.
\end{itemize}
\end{enumerate}

\subsection{Test with $D=1$}

Here we consider three different numbers of stochastic variables,
namely $N=1,5,10$. Results are shown in Figure \ref{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 \ref{thm:misc_complexity}: in this case,
$\zeta = \frac{\dRate{}_1}{\gamma_1} = \frac{\widetilde{\dRate{}}_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.

\begin{figure}[htbp]
  \centering
  \includegraphics[page=2, width=0.8\textwidth]{COOK_d1_N1.pdf}\\
  \includegraphics[page=2, width=0.8\textwidth]{COOK_d1_N5.pdf}\\
  \includegraphics[page=2, width=0.8\textwidth]{COOK_d1_N10.pdf}\\
  \caption{Results for test $D=1$, case $N=1$ (top), $N=5$ (center) and $N=10$ (bottom).}
  \label{fig:res_D=1}
\end{figure}

\subsection{Test with $D=3$}
In this case, we obtain the convergence curves shown in
Figure \ref{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{\dRate{}}_1 = \widetilde{\dRate{}}_2 = \widetilde{\dRate{}}_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 \ref{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 \eqref{eq:error-decomp}.  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.

\begin{figure}[htbp]
  \centering
  \includegraphics[page=2,width=0.8\textwidth]{COOK_d3_N1.pdf}\\
  \includegraphics[page=2,width=0.8\textwidth]{COOK_d3_N5.pdf}\\
  \includegraphics[page=2,width=0.8\textwidth]{COOK_d3_N10.pdf}
  \caption{Results for test $D=3$, case $N=1$ (top), $N=5$ (center)
    and $N=10$ (bottom). }
  \label{fig:res_D=3}
\end{figure}

 %
\section{Conclusions}\label{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.

\paragraph{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.