# Lasso Monte Carlo, a Variation on Multi Fidelity Methods for High
Dimensional Uncertainty Quantification

Arnau Albà  Paul Scherrer Institut, Forschungstrasse 111, 5232 Villigen,
Switzerland  ETH Zürich, Rämistrasse 101, 8092 Zurich Switzerland  Romana
Boiger  Paul Scherrer Institut, Forschungstrasse 111, 5232 Villigen,
Switzerland  Dimitri Rochman  Paul Scherrer Institut, Forschungstrasse 111,
5232 Villigen, Switzerland  Andreas Adelmann<EMAIL_ADDRESS>Paul
Scherrer Institut, Forschungstrasse 111, 5232 Villigen, Switzerland

###### Abstract

Uncertainty quantification (UQ) is an active area of research, and an
essential technique used in all fields of science and engineering. The most
common methods for UQ are Monte Carlo and surrogate-modelling. The former
method is dimensionality independent but has slow convergence, while the
latter method has been shown to yield large computational speedups with
respect to Monte Carlo. However, surrogate models suffer from the so-called
curse of dimensionality, and become costly to train for high-dimensional
problems, where UQ might become computationally prohibitive. In this paper we
present a new technique, Lasso Monte Carlo (LMC), which combines a Lasso
surrogate model with the multifidelity Monte Carlo technique, in order to
perform UQ in high-dimensional settings, at a reduced computational cost. We
provide mathematical guarantees for the unbiasedness of the method, and show
that LMC can be more accurate than simple Monte Carlo. The theory is
numerically tested with benchmarks on toy problems, as well as on a real
example of UQ from the field of nuclear engineering. In all presented examples
LMC is more accurate than simple Monte Carlo and other multifidelity methods.
Thanks to LMC, computational costs are reduced by more than a factor of 5 with
respect to simple MC, in relevant cases.

##  1 Introduction

Uncertainty Quantification (UQ) aims to calculate the effect of unknown or
uncertain system parameters on the outcome of an experiment or computation. It
is an active area of research, and an essential tool to test the robustness
and accuracy of methods used in many domains of science and engineering, such
as risk assessment in civil engineering [1], design and optimisation of
particle accelerators [2, 3], weather prediction [4], medical physics [5], and
nuclear engineering [6, 7, 8].

The UQ process can be described as follows: let $f(\boldsymbol{x})$ be a
deterministic function which represents the numerical experiment

$f\colon\begin{array}[t]{ >{\displaystyle}r >{{}}c<{{}} >{\displaystyle}l
}\mathbb{R}^{d}&\to&\mathbb{R}^{m}\\\
\boldsymbol{x}&\mapsto&f(\boldsymbol{x})\,,\end{array}$

with input and output dimensions of size $d$ and $m$, respectively. Let
$\boldsymbol{x}_{0}=(x_{1},x_{2},...,x_{d})$ be an input vector, with an
associated uncertainty. The uncertainty can be modelled by letting the input
be random variable $X$, centred around $\boldsymbol{x}_{0}$, such that
$\operatorname*{\mathbb{E}}[X]=\boldsymbol{x}_{0}$. A common approach is to
model the input with a multivariate normal distribution
$X\sim\mathcal{N}(\boldsymbol{x}_{0},\Sigma)$, where $\Sigma$ is a known
covariance matrix. The aim of UQ, and more specifically response variability
methods, is to estimate the mean $\mu$ and variance $\sigma^{2}$ of the output
distribution $f(X)$ which is then written as
$f(\boldsymbol{x}_{0})=\mu\pm\sigma$. Without loss of generality, in the rest
of the paper it is assumed that the output is one-dimensional, i.e. $m=1$.

When $f$ is a black-box function one has to rely on non-intrusive UQ methods,
such as Monte Carlo (henceforth referred to as simple MC) [9, 10] or surrogate
modelling [11, 12].

With simple MC, $N$ independent and identically distributed (i.i.d).
multivariate random variables $X_{1},X_{2},...,X_{N}$ are sampled to obtain a
set of input vectors
$\boldsymbol{x}_{1},\boldsymbol{x}_{2},...,\boldsymbol{x}_{N}$. Then $f$ is
evaluated at each input to obtain a set of outputs
$f(\boldsymbol{x}_{1}),f(\boldsymbol{x}_{2}),...,f(\boldsymbol{x}_{N})$, from
which the sample mean and sample variance are calculated. Monte Carlo methods
are known to converge as $\operatorname*{\mathcal{O}}(N^{-{\frac{1}{2}}})$
which, if $f$ is expensive to evaluate, can make such methods computationally
expensive or even prohibitive.

The slow convergence of the errors with simple MC can in some cases be
bypassed by using surrogate models [11]. With this approach, the $N$ input-
output samples of $f$ are used to train a surrogate model $f^{(N)}$. Then
$f^{(N)}$ is evaluated $M$ times, and the outputs are used to estimate the
mean and variance. The advantage of this method is that $f^{(N)}$ is
computationally cheaper than $f$, and evaluating it $M$ times with $M\gg N$
has negligible runtime. The bottleneck of this method is in obtaining the $N$
samples for the training set. If $f$ is low-dimensional such that $N\gg d$,
the surrogate model is likely to have a small bias, however in high-
dimensional cases with $N<d$ the surrogate will be biased (see the curse of
dimensionality [13, 14]), making also the mean and variance estimations
biased. The specific case where $N<d$, and where increasing $N$ is not
possible or computationally expensive, is the main focus of this paper. The
aim is to find a method more accurate than simple MC for a fixed computational
budget $N$.

In this regard, multifidelity Monte Carlo (MFMC) [15, 16, 17, 18] offers a
promising approach. Multifidelity methods combine low and high fidelity models
to accelerate the solution of an outer-loop application, like uncertainty
quantification, sensitivity analysis or optimisation. The low fidelity models,
that are either simplified models, projection-based models or data-fit
surrogates, are used for speed up, while the high fidelity models are kept in
the loop for accuracy and/or convergence. Thus MFMC offers a framework to
combine a biased surrogate model $f^{(n)}$ with the high fidelity model $f$,
in such a manner that unbiased estimates of the mean and variance are
computed. The crucial point with MFMC is then, for a given number of high
fidelity evaluations $N$, optimising the trade-off between how many of them
are used in the multifidelity estimators and how many for training the
surrogate. This optimisation problem has been recently addressed in [19].
Nevertheless, there is no guarantee that, for a given $N$, the MFMC estimates
will be more accurate than simple MC, especially in high-dimensional cases
with $N<d$, where $f^{(n)}$ is likely to have a large bias.

To address the challenges presented by high-dimensional UQ in existing
approaches, we introduce the Lasso Monte Carlo (LMC) method, a variation on
MFMC. With LMC we propose a new data management strategy in which the high
fidelity samples are reused several times both to train multiple surrogates
and in the multifidelity estimators. The resulting algorithm is such that the
estimations are guaranteed to be equally or more accurate than simple MC and
MFMC, under certain assumptions. This new approach can be viewed as a variance
reduction technique, based on a two-level version of MFMC, and Lasso
regression (least absolute shrinkage and selection operator), simultaneously
performing regression analysis, variable selection and regularization [20].

It is worth mentioning the relation between multifidelity, multilevel, and
multi-indexing methods: Multilevel Monte Carlo (MLMC) [21] is a special case
within the broader framework of multifidelity methods, in which typically a
hierarchy of low fidelity models is derived by varying a parameter (e.g. mesh
width). A generalization of MLMC is introduced in [22], so called multi-index
Monte Carlo methods, which use multidimensional levels and high-order mixed
differences to reduce the computational cost and the variance of the resulting
estimator.

The remainder of the paper is organised as follows: we start with a review of
MFMC for the estimation of central moments, and in particular we derive the
expressions for the two-level estimators. This is followed by a discussion on
the trade-off between accuracy and computational costs of MFMC. The new method
LMC is then introduced, and we prove that it is equally or more accurate than
simple MC. We then review the theory behind the Lasso regression method, and
show how it can be used in the LMC algorithm. Finally, in section 3 LMC is
benchmarked on a variety of examples. Proofs for all the theorems and lemmas
are provided in appendix A.

###  1.1 Notation & Assumptions

Throughout the paper, bold letters represent vectors
$\boldsymbol{x}=(x_{1},x_{2},...,x_{d})$, with $d$ the dimension. A lower case
letter $x$ is a realisation of a random variable $X$, and in the case where
$X$ is a multivariate random variable of dimension $d$, a realisation of it
will be a vector $(x_{1},x_{2},...,x_{d})$. For a random variable $X$ with
probability density function (PDF) $\phi(x)$, we calculate the expectation
value or mean with
$\operatorname*{\mathbb{E}}[X]=\int_{\mathbb{R}}x\,\phi(x)\,dx\,$. We also use
$\operatorname*{\mathbb{E}}[f(X)]$ or $\operatorname*{\mathbb{E}}[f]$ for the
mean of the function $f$, whose input follows the distribution of $X$. The
variance is defined as
$\operatorname*{\mathrm{Var}}[X]=\operatorname*{\mathbb{E}}\Big{[}\big{(}X-\operatorname*{\mathbb{E}}[X]\big{)}^{2}\Big{]}$,
the covariance
$\operatorname*{\mathrm{Cov}}[X,Y]=\operatorname*{\mathbb{E}}\Big{[}\big{(}X-\operatorname*{\mathbb{E}}[X]\big{)}\big{(}Y-\operatorname*{\mathbb{E}}[Y]\big{)}\Big{]}$,
the fourth central moment
$m_{4}[X]=\operatorname*{\mathbb{E}}\Big{[}\big{(}X-\operatorname*{\mathbb{E}}[X]\big{)}^{4}\Big{]}$,
and finally the multivariate second moment
$m_{2,2}[X,Y]=\operatorname*{\mathbb{E}}\Big{[}\big{(}X-\operatorname*{\mathbb{E}}[X]\big{)}^{2}\big{(}Y-\operatorname*{\mathbb{E}}[Y]\big{)}^{2}\Big{]}$.

The function space
$L^{p}\big{(}\mathbb{R}^{d},\phi(\boldsymbol{x})d\boldsymbol{x}\big{)}$ for
$1\leq p<\infty$ is defined as the space of functions satisfying

$\left(\int_{\mathbb{R}^{d}}|f(\boldsymbol{x})|^{p}\phi(\boldsymbol{x})d\boldsymbol{x}\right)^{1/p}<\infty\,$.

The function $f$ is the ground truth, i.e. the expensive model, while
$f^{(n)}$ is a cheap-to-evaluate surrogate model that was fitted to a training
set of $n$ samples of $f$. Similarly, $\mu_{N}$ and $\mu_{N}^{(n)}$ are the
sample estimators for the mean of a sample set of size $N$, calculated with
the true and surrogate model respectively. Also
$\operatorname*{{\sigma^{2}}}2_{N}$ and
$\operatorname*{{\sigma^{2}}}2_{N}^{(n)}$ are the sample estimators of the
variance of a sample set of size $N$, computed with the true and surrogate
models respectively. We write the two-level estimators as $\mu_{N,M}^{(n)}$
and $\operatorname*{{\sigma^{2}}}2_{N,M}^{(n)}$, and the LMC estimators as
$\mathcal{M}_{N,M}$ and $\Sigma^{2}_{N,M}$. It is assumed that the
computational cost of training and evaluating $f^{(n)}$ is negligible compared
to the cost of evaluating $f$. Therefore the cost of an estimator is given by
$N$, the number of evaluations of $f$.

A normal distribution with mean $\boldsymbol{\mu}$ and covariance matrix
$\Sigma$ is written as $\mathcal{N}(\boldsymbol{\mu},\Sigma)$, and $U[a,b]$ is
a uniform distribution between $a$ and $b$ with $a<b$. With mean squared error
(MSE) is defined as
$\operatorname*{\text{MSE}}\left(\mu_{N},\operatorname*{\mathbb{E}}[f]\right)=\operatorname*{\mathbb{E}}\left[\left(\mu_{N}-\operatorname*{\mathbb{E}}[f]\right)^{2}\right]$.

Without loss of generality, we assume that all the sets of data used have been
centred around zero, such that
$\frac{1}{N}\sum_{i=1}^{N}f(\boldsymbol{x}_{i})=0\,,\quad\text{and}\quad\frac{1}{N}\sum_{i=1}^{N}x_{ik}=0\,,\quad
k=1,2,...,d\,$.

The following assumption is made for any surrogate model $f^{(n)}$, and
section 2.4 discusses its validity.

###### Assumption 1.1

Let $f^{(n)}$ be a surrogate model that has been fitted on a training set of
$n$ inputs $\boldsymbol{x}_{1},\boldsymbol{x}_{2},...,\boldsymbol{x}_{n}$
sampled from i.i.d. random variables distributed as $X$, and $n$ outputs
$f(\boldsymbol{x}_{1}),f(\boldsymbol{x}_{2}),...,f(\boldsymbol{x}_{n})$. Then
the following inequalities are satisfied

$\displaystyle\operatorname*{\mathrm{Var}}\left[f-f^{(n)}\right]\leq\leavevmode\nobreak\
\operatorname*{\mathrm{Var}}\left[f\right]\,,$ (1a)
$\displaystyle\begin{split}&m_{2,2}\left[f+f^{(n)},f-f^{(n)}\right]+\frac{1}{N-1}\operatorname*{\mathrm{Var}}\left[f+f^{(n)}\right]\operatorname*{\mathrm{Var}}\left[f-f^{(n)}\right]\\\
&\quad-\frac{N-2}{N-1}\left(\operatorname*{\mathrm{Var}}\left[f\right]-\operatorname*{\mathrm{Var}}\left[f^{(n)}\right]\right)^{2}\leq\leavevmode\nobreak\
m_{4}[f]-\frac{N-3}{N-1}{\operatorname*{\mathrm{Var}}}^{2}[f]\,,\end{split}$
(1b)

for any integers $n>0$ and $N>3$.

##  2 Theory

###  2.1 Two-Level Monte Carlo

The proposed method, LMC, is based on multilevel and multifidelity Monte Carlo
methods [21, 16, 15]. As with any MC method, one wants to estimate
$\operatorname*{\mathbb{E}}[f]$ for some function $f$. Several models
$f_{1},f_{2},...,f_{L}$ are available that approximate $f$, which have
increasing cost and increasing accuracy, i.e. $f_{L}$ is the most accurate and
expensive model, while $f_{1}$ is computationally cheap but inaccurate. The
difference in naming between multilevel and multifidelity methods comes from
how these levels of accuracy are defined: in multilevel Monte Carlo the levels
are obtained by coarsening or refining a grid, or changing the step-size of
the integrator, whereas in multifidelity Monte Carlo the different functions
are given by any general lower-order model such as data-fit surrogates or
projection-based models. In both cases, however, the goal of the method is to
reduce the overall computational cost of computing
$\operatorname*{\mathbb{E}}[f]$ with respect to traditional MC, by optimally
balancing the amount of evaluations at each level of accuracy.

Originally MLMC was developed for estimating only the mean of a distribution,
but in more recent years it has been extended [23, 24, 25] to estimate higher
order moments of the distribution, which are necessary for UQ.

In the following paragraphs, a two-level version of MFMC is derived, in which
two levels of accuracy are considered: the expensive, unbiased, true model
$f$, and a computationally cheap, biased, surrogate model $f^{(n)}$.

####  2.1.1 Mean Estimator

Let $f(\boldsymbol{x})\in
L^{2}\big{(}\mathbb{R}^{d},\phi(\boldsymbol{x})d\boldsymbol{x}\big{)}$ be a
function, whose input is distributed according to the multivariate
$d$-dimensional random variable $X$, with probability density function (PDF)
$\phi(\boldsymbol{x})$. Let
$\boldsymbol{x}_{1},\boldsymbol{x}_{2},...,\boldsymbol{x}_{N}$ and
$\boldsymbol{z}_{1},\boldsymbol{z}_{2},...,\boldsymbol{z}_{M}$ be two sets of
input samples of size $N$ and $M$, drawn from the i.i.d. random variables
distributed as $X$. The aim is to estimate the mean
$\operatorname*{\mathbb{E}}[f]$, with the minimum number of evaluations of the
function $f$. The simple MC estimator is

$\mu_{N}=\frac{1}{N}\sum_{i=1}^{N}f(\boldsymbol{x}_{i})\,.$ (2)

The mean squared error of the estimator is

$\operatorname*{\text{MSE}}\Big{(}\mu_{N},\operatorname*{\mathbb{E}}[f]\Big{)}=\frac{\operatorname*{\mathrm{Var}}[f]}{N}\,.$
(3)

Therefore it is an unbiased estimator since
$\lim_{N\rightarrow\infty}\mu_{N}=\operatorname*{\mathbb{E}}[f]$. With this
method we require
$N=\frac{\operatorname*{\mathrm{Var}}[f]}{\operatorname*{\text{MSE}}}$ samples
to obtain an estimation with a mean squared error of
$\operatorname*{\text{MSE}}$.

Now let $f^{(N)}\in
L^{2}\big{(}\mathbb{R}^{d},\phi(\boldsymbol{x})d\boldsymbol{x}\big{)}$ be a
surrogate model that was trained on $N$ evaluations of $f$, and is much
cheaper to evaluate. Using this surrogate model to compute the sample mean,
the estimator is

$\mu_{M}^{(N)}=\frac{1}{M}\sum_{i=1}^{M}f^{(N)}(\boldsymbol{z}_{i})\,.$ (4)

Note that this estimator has the same cost as (2) since the number of
evaluations of $f$ is the same. The error is

$\operatorname*{\text{MSE}}\Big{(}\mu^{(N)}_{M},\operatorname*{\mathbb{E}}[f]\Big{)}=\left(\operatorname*{\mathbb{E}}[f^{(N)}]-\operatorname*{\mathbb{E}}[f]\right)^{2}+\frac{\operatorname*{\mathrm{Var}}\left[f^{(N)}\right]}{M}\,.$
(5)

The first term is the bias, and the second term is the variance. The variance
term quickly vanishes if we assume that $f^{(N)}$ has a negligible runtime,
and that thus $M\rightarrow\infty$. However, the bias term can only be reduced
by increasing the training set size $N$ (and hence improving the accuracy of
the surrogate model), which might be impossible or computationally demanding.
This is especially problematic in high-dimensional cases, since the volume of
the input space to be sampled increases exponentially with $d$ [14], and
unless the training set was large $N\gg d$ the surrogate will be heavily
biased. Additionally, even if the training set were large and
$N\rightarrow\infty$, the bias would in general not decay to zero due to model
bias.

Let us now combine the surrogate model with $f$ into a two-level estimator.
For this, assume that the set of $N$ samples is split into a training subset
$\boldsymbol{x}_{1},\boldsymbol{x}_{2},...,\boldsymbol{x}_{n}$ of size $n$
with $1\leq n\leq N$, and an evaluation subset
$\boldsymbol{x}_{n+1},\boldsymbol{x}_{n+2},...,\boldsymbol{x}_{N}$ of size
$N-n$. The training samples are used to fit a surrogate $f^{(n)}$. Then the
two-level estimator reads

$\mu_{N-n,M}^{(n)}=\frac{1}{M}\sum_{i=1}^{M}f^{(n)}(\boldsymbol{z}_{i})+\frac{1}{N-n}\sum_{i=n+1}^{N}f(\boldsymbol{x}_{i})-f^{(n)}(\boldsymbol{x}_{i})=\mu^{(n)}_{M}+\mu_{N-n}-\mu_{N-n}^{(n)}\,.$
(6)

As with the previous estimator, the cost is $N$, but in this case it is an
unbiased estimate since the error

$\operatorname*{\text{MSE}}\Big{(}\mu_{N-n,M}^{(n)},\operatorname*{\mathbb{E}}[f]\Big{)}=\frac{\operatorname*{\mathrm{Var}}\left[f^{(n)}\right]}{M}+\frac{\operatorname*{\mathrm{Var}}\left[f-f^{(n)}\right]}{N-n}$
(7)

has variance terms, but no bias term. Here again we assume that $f^{(n)}$ has
negligible runtime and that $M\rightarrow\infty$, and thus the first term
vanishes. Then, the following statement can be made regarding the MSE of
simple MC and of the two-level estimator:

$\lim_{M\rightarrow\infty}\operatorname*{\text{MSE}}\Big{(}\mu_{N-n,M}^{(n)},\operatorname*{\mathbb{E}}[f]\Big{)}\leq\operatorname*{\text{MSE}}\Big{(}\mu_{N},\operatorname*{\mathbb{E}}[f]\Big{)}\quad\iff\quad\frac{N}{N-n}\leq\frac{\operatorname*{\mathrm{Var}}[f]}{\operatorname*{\mathrm{Var}}\left[f-f^{(n)}\right]}\,.$
(8)

Therefore, for a given computational budget $N$, the two-level estimator is
equally or more accurate than simple MC if and only if $n$ and $f^{(n)}$ are
such that (8) is satisfied. Note that the fraction on the far right of (8) is
guaranteed to be larger than $1$ by assumption 1a. Also in this fraction note
that the denominator could in principle be zero, however we assume that this
will never happen due to the model bias of $f^{(n)}$.

####  2.1.2 Variance Estimation

Let $f(\boldsymbol{x})\in
L^{4}\big{(}\mathbb{R}^{d},\phi(\boldsymbol{x})d\boldsymbol{x}\big{)}$. Let
there be two sets of input samples of size $N$ and $M$, as in section 2.1.1.
Then the simple MC estimator for the variance is

$\sigma^{2}_{N}=\frac{1}{N-1}\sum_{i=1}^{N}\left(f(\boldsymbol{x}_{i})-\sum_{j=1}^{N}\frac{f(\boldsymbol{x}_{j})}{N}\right)^{2}\,,$
(9)

which is unbiased and has an error

$\operatorname*{\text{MSE}}\Big{(}\sigma^{2}_{N},\operatorname*{\mathrm{Var}}[f]\Big{)}=\frac{1}{N}\left(m_{4}[f]-\frac{N-3}{N-1}{\operatorname*{\mathrm{Var}}}^{2}[f]\right)\,.$
(10)

Using a surrogate model $f^{(N)}\in
L^{4}\big{(}\mathbb{R}^{d},\phi(\boldsymbol{x})d\boldsymbol{x}\big{)}$ trained
on $N$ evaluations of $f$, the variance estimator is

${\operatorname*{{\sigma^{2}}}2}_{M}^{(N)}=\frac{1}{M-1}\sum_{i=1}^{M}\left(f^{(N)}(\boldsymbol{x}_{i})-\sum_{j=1}^{M}\frac{f^{(N)}(\boldsymbol{x}_{j})}{M}\right)^{2}\,,$
(11)

which has an error

$\operatorname*{\text{MSE}}\Big{(}{\operatorname*{{\sigma^{2}}}2}_{M}^{(N)},\operatorname*{\mathrm{Var}}[f]\Big{)}=\Big{(}\operatorname*{\mathrm{Var}}[f^{(N)}]-\operatorname*{\mathrm{Var}}[f]\Big{)}^{2}+\frac{1}{M}\left(m_{4}[f^{(N)}]-\frac{M-3}{M-1}{\operatorname*{\mathrm{Var}}}^{2}[f^{(N)}]\right)\,.$

The surrogate estimation of the variance is biased. The second term vanishes
if $M\rightarrow\infty$, while the first term is affected by model bias, and
decays slowly due to the curse of dimensionality.

Now assume that the $N$ input samples are split into a subset of $n$ training
samples which are used to fit a surrogate $f^{(n)}$, and a subset of $N-n$
evaluation samples. Then the two-level estimator for the variance is

${\operatorname*{{\sigma^{2}}}2}^{(n)}_{N-n,M}={\operatorname*{{\sigma^{2}}}2}_{M}^{(n)}+\sigma^{2}_{N-n}-{\operatorname*{{\sigma^{2}}}2}_{N-n}^{(n)}\,,$
(12)

This estimator has a computational cost of $N$ and is unbiased. The error is

$\begin{split}\operatorname*{\text{MSE}}&\left({\operatorname*{{\sigma^{2}}}2}^{(n)}_{N-n,M},\operatorname*{\mathrm{Var}}[f]\right)\\\
&=\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}^{(n)}_{M}\right]+\operatorname*{\mathrm{Var}}\left[\sigma^{2}_{N-n}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N-n}\right]\\\
&=\frac{1}{M}\left(m_{4}[f^{(n)}]-\frac{M-3}{M-1}{\operatorname*{\mathrm{Var}}}^{2}[f^{(n)}]\right)\\\
&\quad+\frac{1}{N-n}\Big{(}m_{2,2}\left[f+f^{(n)},f-f^{(n)}\right]+\frac{1}{N-n-1}\operatorname*{\mathrm{Var}}[f+f^{(n)}]\operatorname*{\mathrm{Var}}[f-f^{(n)}]\\\
&\hskip
56.9055pt-\frac{N-n-2}{N-n-1}\left(\operatorname*{\mathrm{Var}}[f]-\operatorname*{\mathrm{Var}}[f^{(n)}]\right)^{2}\Big{)}\,.\end{split}$
(13)

Similarly to the mean estimator, the error (13) has two variance terms. We
assume that the first term vanishes as $M\rightarrow\infty$. Then, the
following statement can be made regarding the MSE of simple MC and of the two-
level estimator:

$\begin{split}&\lim_{M\rightarrow\infty}\operatorname*{\text{MSE}}\left({\operatorname*{{\sigma^{2}}}2}^{(n)}_{N-n,M},\operatorname*{\mathrm{Var}}[f]\right)\leq\operatorname*{\text{MSE}}\Big{(}\sigma^{2}_{N},\operatorname*{\mathrm{Var}}[f]\Big{)}\quad\iff\\\
&\frac{N}{N-n}\leq\Big{(}m_{4}[f]-\frac{N-3}{N-1}{\operatorname*{\mathrm{Var}}}^{2}[f]\Big{)}\Big{(}m_{2,2}\left[f+f^{(n)},f-f^{(n)}\right]\\\
&\hskip
34.14322pt+\frac{1}{N-n-1}\operatorname*{\mathrm{Var}}[f+f^{(n)}]\operatorname*{\mathrm{Var}}[f-f^{(n)}]-\frac{N-n-2}{N-n-1}\left(\operatorname*{\mathrm{Var}}[f]-\operatorname*{\mathrm{Var}}[f^{(n)}]\right)^{2}\Big{)}^{-1}\,.\end{split}$
(14)

Therefore, for a given $N$, the two-level estimator is equally or more
accurate than simple MC if and only if $n$ and $f^{(n)}$ are such that (14) is
satisfied.

###  2.2 Lasso Monte Carlo (LMC)

Consider a fixed computational budget of $N$ evaluations of $f$, and a fixed
choice of surrogate model. Then, the mean squared errors of the two-level
estimators (7, 13) are controlled by $n$. There is a trade-off in the choice
of this number, a larger $n$ increases the training set size and potentially
reduces the bias of $f^{(n)}$, while it also increases the variance of the
estimator since the number of evaluation samples $N-n$ is decreased. In the
default implementation of MFMC (henceforth referred to as static MFMC), the
choice of $n$ and the training of the surrogate are done beforehand,
independently of the MSE. In recent years adaptive MFMC [19] (sometimes also
referred to as context-aware surrogate models [26]) were introduced. With this
approach, an upper bound of the MSE is estimated, and used to calculate the
optimal $n$ that minimises the errors. Both of these approaches converge
faster than simple MC [15, 19], however there is no guarantee that the MFMC
estimates are more accurate than those of simple MC. Indeed, for a given $N$,
it is possible that there exists no choice of $n$ in $\\{1,...,N\\}$ such that
conditions (8, 14) are satisfied.

A possible way to guarantee that the two-level estimators are more accurate
than simple MC, is by considering a larger computational budget of $n+N$. Let
$\boldsymbol{x}_{1},\boldsymbol{x}_{2},...,\boldsymbol{x}_{n}$ be the training
set and $\boldsymbol{x}_{n+1},\boldsymbol{x}_{n+2},...,\boldsymbol{x}_{n+N}$
be the evaluation set. Then the estimators (6,12) become $\mu^{(n)}_{N,M}$ and
${\operatorname*{{\sigma^{2}}}2}^{(n)}_{N,M}$, which would have the following
errors:

$\begin{split}\operatorname*{\text{MSE}}&\Big{(}\mu_{N,M}^{(n)},\operatorname*{\mathbb{E}}[f]\Big{)}=\frac{\operatorname*{\mathrm{Var}}\left[f-f^{(n)}\right]}{N}+\frac{\operatorname*{\mathrm{Var}}\left[f^{(n)}\right]}{M}\\\
\operatorname*{\text{MSE}}&\left({\operatorname*{{\sigma^{2}}}2}^{(n)}_{N,M},\operatorname*{\mathrm{Var}}[f]\right)=\frac{1}{N}\Big{(}m_{2,2}\left[f+f^{(n)},f-f^{(n)}\right]\\\
&\hskip
56.9055pt+\frac{1}{N-1}\operatorname*{\mathrm{Var}}[f+f^{(n)}]\operatorname*{\mathrm{Var}}[f-f^{(n)}]\\\
&\hskip
56.9055pt-\frac{N-2}{N-1}\left(\operatorname*{\mathrm{Var}}[f]-\operatorname*{\mathrm{Var}}[f^{(n)}]\right)^{2}\Big{)}\\\
&\hskip
56.9055pt+\frac{1}{M}\left(m_{4}[f^{(n)}]-\frac{M-3}{M-1}{\operatorname*{\mathrm{Var}}}^{2}[f^{(n)}]\right)\,.\end{split}$
(15)

These estimators are guaranteed to be more accurate than simple MC due to
assumptions 1:

$\begin{split}&\lim_{M\rightarrow\infty}\operatorname*{\text{MSE}}\Big{(}\mu_{N,M}^{(n)},\operatorname*{\mathbb{E}}[f]\Big{)}\leq\operatorname*{\text{MSE}}\Big{(}\mu_{N},\operatorname*{\mathbb{E}}[f]\Big{)}\quad\iff\quad\text{Assumption
}\ref{eq:assumption_a}\,,\\\
&\lim_{M\rightarrow\infty}\operatorname*{\text{MSE}}\left({\operatorname*{{\sigma^{2}}}2}^{(n)}_{N,M},\operatorname*{\mathrm{Var}}[f]\right)\leq\operatorname*{\text{MSE}}\Big{(}\sigma^{2}_{N},\operatorname*{\mathrm{Var}}[f]\Big{)}\quad\iff\quad\text{Assumption
}\ref{eq:assumption_b}\,.\end{split}$ (16)

This approach, despite guaranteeing higher accuracy than simple MC, has a
larger computational cost as it requires $N+n$ evaluations of $f$. One could
naively reduce the costs by reusing the same $N$ samples for training
$f^{(N)}$ and evaluating the estimators, thus removing the need for the
additional $n$ training samples. This, however, would introduce bias in the
estimator since the evaluation and training samples must come from i.i.d.
random variables. The LMC algorithm introduced in this work presents an
alternative method, with computational cost $N$, and MSE similar to (15). To
avoid introducing bias, LMC splits the training set into $S$ subsets, and
trains $S$ surrogate models. The full algorithm is detailed in (Algorithm 1).
Note that the LMC algorithm has Lasso in its name, because in this paper the
surrogate model $f^{(n)}$ used in the algorithm was exclusively a Lasso model.
However, algorithm 1 does not make any reference to Lasso and is in principle
agnostic to the method used to train $f^{(n)}$.

Algorithm 1 The Lasso Monte Carlo Algorithm

1:Input sets
$V=\\{\boldsymbol{x}_{1},\boldsymbol{x}_{2},...,\boldsymbol{x}_{N}\\}$ and
$W=\\{\boldsymbol{z}_{1},\boldsymbol{z}_{2},...,\boldsymbol{z}_{M}\\}$, and an
integer $S$ that divides $N$.

2:$N\ll M$

3:Compute the outputs of set $V$:
$f(\boldsymbol{x}_{1}),f(\boldsymbol{x}_{2}),...,f(\boldsymbol{x}_{N})$.

4:Compute $\mu_{N},\sigma^{2}_{N}$ with the simple MC estimators (2, 9).

5:for $\operatorname*{\mathbcal{s}}=1\ldots S$ do

6: Define an evaluation subset of size $\frac{N}{S}$ as
$V_{eval}=\\{\boldsymbol{x}_{N(\operatorname*{\mathbcal{s}}-1)/S+1},\boldsymbol{x}_{N(\operatorname*{\mathbcal{s}}-1)/S+2},...,\boldsymbol{x}_{N\operatorname*{\mathbcal{s}}/S}\\}\subset
V$.

7: The remaining samples define the training subset $V_{train}=V\setminus
V_{eval}$ of size $n=N\frac{S-1}{S}$.

8: Fit $f^{(n)}_{\operatorname*{\mathbcal{s}}}$ to $V_{train}$.

9: With $f^{(n)}_{\operatorname*{\mathbcal{s}}}$, compute the outputs of $W$
and $V_{eval}$.

10: Use $W$ and $V_{eval}$ to compute
$\left(\mu^{(n)}_{N/S,M}\right)_{\operatorname*{\mathbcal{s}}}$ and
$\left({\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S,M}\right)_{\operatorname*{\mathbcal{s}}}$
with the two-level estimators (6,12).

11:end for

12:Compute the LMC mean and variance, by averaging out the two-level
estimations:
$\mathcal{M}_{N,M}=\frac{1}{S}\sum_{\operatorname*{\mathbcal{s}}=1}^{S}\left(\mu^{(n)}_{N/S,M}\right)_{\operatorname*{\mathbcal{s}}}\,,\quad\text{and}\quad\Sigma^{2}_{N,M}=\frac{1}{S}\sum_{\operatorname*{\mathbcal{s}}=1}^{S}\left({\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S,M}\right)_{\operatorname*{\mathbcal{s}}}\,.$
(17)

This algorithm only requires $N$ evaluations of $f$, and therefore has the
same cost as simple MC. Note that the LMC estimators (17) differ from the two-
level estimators (6, 12), since LMC computes the average of several smaller
two-level estimators. It can be shown that, under certain assumptions, the LMC
estimators have similar accuracy to (15), in which case they are guaranteed to
be equally or more accurate than simple MC (see eq. 16). Furthermore, in the
benchmarks section (sec. 3) LMC is found to be more accurate than simple MC,
static MFMC, and adaptive MFMC for all $N$.

######  Theorem 2.1 (LMC Accuracy)

Let $f^{(n)}_{1}=f^{(n)}_{2}=...=f^{(n)}_{S}$. Then the LMC estimators
$\mathcal{M}_{N,M}$ and $\Sigma^{2}_{N,M}$ defined by (17), have the errors
$\lim_{M\rightarrow\infty}\operatorname*{\text{MSE}}\left(\mu^{(n)}_{N,M},\operatorname*{\mathbb{E}}[f]\right)=\lim_{M\rightarrow\infty}\operatorname*{\text{MSE}}\left(\mathcal{M}_{N,M},\operatorname*{\mathbb{E}}[f]\right)$
and
$\lim_{M\rightarrow\infty}\operatorname*{\text{MSE}}\left({\operatorname*{{\sigma^{2}}}2}^{(n)}_{N,M},\operatorname*{\mathrm{Var}}[f]\right)=\lim_{M\rightarrow\infty}\operatorname*{\text{MSE}}\left(\Sigma^{2}_{N,M},\operatorname*{\mathrm{Var}}[f]\right)+\mathcal{O}\left(N^{-2}\right)$.

###### Remark 2.1

The condition $f^{(n)}_{1}=f^{(n)}_{2}=...=f^{(n)}_{S}$ occurs if a surrogate
model $f^{(n)}$ is independent of the elements of the training set, and only
depends on the size $n$ of the training set. This occurs when the training set
is large, or $f^{(n)}$ is a model with small variance.

As a more intuitive interpretation of LMC, we point out that the definition of
$V_{eval}$ and $V_{train}$ in the for-loop of algorithm 1, is equivalent to
performing an S-fold split on the set $V$, as described in figure 1.

Figure 1: Illustration of the LMC algorithm 1, with $S=5$ and $N=500$.
Accordingly, the evaluation sets $V_{eval}$ have size $100$, and the training
sets $V_{train}$ have size $n=400$.

###  2.3 Lasso Regression

In remark 2.1 it was pointed out that a surrogate model with a small variance
is required for the LMC method. For this reason, in this work it was decided
to use a Lasso regression model [20], as it has a regularisation parameter
$\lambda$ for straight-forward tuning of the model variance. This section
introduces Lasso.

The Lasso regression method [20] estimates the weights
$\boldsymbol{\beta}\in\mathbb{R}^{d}$ of a linear model of the type

$f^{(n)}(\boldsymbol{x})=\boldsymbol{\beta}\cdot\boldsymbol{x},$ (18)

such that $f^{(n)}$ approximates $f$. This method estimates
$\boldsymbol{\beta}$ by minimising a loss function on a training set
$\boldsymbol{x}_{1},f(\boldsymbol{x}_{1})$,
$\boldsymbol{x}_{2},f(\boldsymbol{x}_{2})$, …,
$\boldsymbol{x}_{n},f(\boldsymbol{x}_{n})$ of size $n$. The Lasso loss
function is

$\mathcal{L}(\boldsymbol{\beta})=\frac{1}{2}\sum_{i=1}^{n}\left(f(\boldsymbol{x}_{i})-\boldsymbol{\beta}\cdot\boldsymbol{x}_{i}\right)^{2}+\lambda||\boldsymbol{\beta}||_{1}\,,$
(19)

where $\lambda>0$ is a chosen regularisation constant, and
$||\boldsymbol{\beta}||_{1}$ is the $L^{1}$-norm

$||\boldsymbol{\beta}||_{1}=\sum_{k=1}^{d}|\beta_{k}|\,.$

The first term of the Lasso loss function (19) is the usual ordinary least
squares (OLS) loss function. In an overdetermined system, i.e. with $n>d$, the
OLS loss function is strictly convex and has a global minimum at

$\boldsymbol{\beta}^{\text{OLS}}=\left(A^{T}A\right)^{-1}A^{T}\boldsymbol{f}\,,$
(20)

where $A\in\mathbb{R}^{n\times d}$ is a nonsingular matrix containing the
input training data, and
$\boldsymbol{f}=\left(f(\boldsymbol{x}_{1}),f(\boldsymbol{x}_{2}),...,f(\boldsymbol{x}_{n})\right)$
is the output training data, or in general the quantities of interest.

In an underdetermined system, i.e. $n<d$, the OLS loss is weakly convex, and
therefore the problem is ill-posed and has multiple minima. Moreover, if $A$
is close to collinear the OLS solution (20) has large variance. The Lasso
approach simultaneously solves the problems of underdetermination and large
variance by adding a second term to the loss (19): the regularisation term,
also referred to as the constraint or shrinkage term. This term adds a penalty
on the complexity of the model. More specifically, the Lasso regularisation
has the effect of shrinking the model weights $\boldsymbol{\beta}$ towards
zero, and forcing some of them to be exactly zero, thus yielding a sparse
model. In this sense Lasso can be seen as a variable selection strategy [20],
in which the least important inputs are suppressed, and the dimensionality of
the problem is reduced. The regularisation parameter $\lambda$, chosen by the
user, controls the amount of shrinkage: for large $\lambda$ the weight vector
$\boldsymbol{\beta}$ is small in magnitude and highly sparse, and if $\lambda$
is large enough the weight vector only contains zeros (see lemma 2.1).
Oppositely, if $\lambda$ is small the weight vector approaches the OLS
solution, i.e. $\lambda\rightarrow 0^{+}$ leads to
$\boldsymbol{\beta}\rightarrow\boldsymbol{\beta}^{\text{OLS}}$.

######  Lemma 2.1 (Lasso with Null Weights)

Let $\\{x_{ik}\\}_{i=1,2,...,n}^{k=1,2,...,d}$ be the input data, drawn from
some continuous distribution, and
$f(\boldsymbol{x}_{1}),f(\boldsymbol{x}_{2}),...,f(\boldsymbol{x}_{n})$ the
outputs. Let $\boldsymbol{\beta}$ be the weight vector that minimises the
Lasso loss function (19), with a regularisation parameter $\lambda$. Then,

$\boldsymbol{\beta}=(0,0,...,0)\quad\iff\quad\lambda\geq\lambda_{max}\quad\text{with}\quad\lambda_{max}=\max_{k=1,2,...,d}\left|\sum_{i=1}^{n}x_{ik}f(\boldsymbol{x_{i}})\right|\,.$
(21)

###### Remark 2.2

The minimum of the Lasso loss (19) does not have a closed-form expression, and
therefore an optimisation algorithm is required for finding a solution. In
this work we use the scikit-learn [27] implementation of Lasso, which uses the
coordinate descent algorithm [28, 29].

###  2.4 Lasso in LMC

As per theorem 2.1 and conditions (16), the LMC estimators (17) are equally or
more accurate than simple MC if the surrogate model $f^{(n)}$ satisfies
assumptions (1). Up until now these assumptions have been taken for granted,
and in this section we discuss their validity when the surrogate model is in
fact a Lasso model, as described in section 2.3.

####  2.4.1 Lasso in the LMC Mean Estimator

The following theorem shows that a Lasso surrogate model satisfies assumption
(1a) under mild conditions, and therefore the LMC estimator for the mean is
almost always guaranteed to be equally or more accurate than simple MC.

######  Theorem 2.2 (Lasso in Assumption (1a))

Let $f(\boldsymbol{x})\in
L^{2}\left(\mathbb{R}^{d},\phi(\boldsymbol{x})d\boldsymbol{x}\right)$ be a
function, whose input is a multivariate random variable
$X=(X_{1},X_{2},...,X_{d})$, with a PDF $\phi(\boldsymbol{x})$ and nonsingular
covariance matrix $\Sigma$. Let
$\boldsymbol{x}_{1},f(\boldsymbol{x}_{1}),\boldsymbol{x}_{2},f(\boldsymbol{x}_{2}),...,\boldsymbol{x_{n}},f(\boldsymbol{x}_{n})$
be a random sample of input-outputs, called the training set. Let $f^{(n)}$ be
a surrogate model that was trained on this training set by minimising the
Lasso loss function (19), with a regularisation parameter $\lambda$. Then

  1. _a_)

in the asymptotic limit $n\rightarrow\infty$, $f$ and $f^{(n)}$ satisfy
assumption (1a) for any regularisation parameter $|\lambda|<\infty$,

  2. _b_)

in the non-asymptotic limit $n<\infty$, there exists a regularisation
parameter $0<\lambda<\infty$ such that $f$ and $f^{(n)}$ satisfy assumption
(1a)

###### Remark 2.3

Despite the existence of a $\lambda$ that satisfies (1a) in the non-asymptotic
limit, estimating such a $\lambda$ is a nontrivial problem when $n<\infty$. In
general, estimating $\lambda$ with theoretical guarantees requires prior
knowledge of $f$. In this work we use 5-fold cross validation (CV) [30]
(except in benchmark 3.1, see explanation therein), as suggested in the
original Lasso paper [20], despite CV not having any theoretical guarantees in
the non-asymptotic case. We note that there exist approaches that, under
appropriate conditions of $\\{x_{ik}\\}$ and $f$, can compute the Lasso
estimator with certain guarantees [31, 32].

####  2.4.2 Lasso in the LMC Variance estimator

If assumption (1b) is satisfied then the LMC estimator for the variance is
guaranteed to be equally or more accurate than simple MC. However, in lemma
2.2 we show that, a Lasso model sometimes fails to satisfy this assumption.

Nevertheless, as can be seen in section 3, LMC is found to be more accurate
than simple MC and MFMC in all presented cases. The success of LMC can be
explained with the following heuristic argument: in UQ problems the input of
$f(\boldsymbol{x})$ is the perturbation of some nominal value, i.e.
$\boldsymbol{x}=\boldsymbol{x}_{0}+\boldsymbol{\xi}$, with $\boldsymbol{\xi}$
a random value that represents the uncertainty of the input. Then, in the
vicinity of $\boldsymbol{x}_{0}$, the problem can be Taylor expanded into a
linear function with some additional smaller terms
$f(\boldsymbol{x})=f(\boldsymbol{x}_{0})+\boldsymbol{\xi}\cdot\boldsymbol{\nabla}f(\boldsymbol{x}_{0})+\mathcal{O}\left(||\boldsymbol{\xi}||^{2}\right)\,$.
In theorem 2.3 we show that if $f$ is a linear function with added random
noise, then Lasso does satisfy assumption 1b, in which case LMC is guaranteed
to have an equally or smaller MSE than simple MC.

######  Lemma 2.2 (Lasso in Assumption (1b))

Let $f(\boldsymbol{x})\in
L^{4}\left(\mathbb{R}^{d},\phi(\boldsymbol{x})d\boldsymbol{x}\right)$ be a
function, whose input is a multivariate random variable
$X=(X_{1},X_{2},...,X_{d})$. Let $\boldsymbol{x}_{1},f(\boldsymbol{x}_{1})$
$,\boldsymbol{x}_{2},f(\boldsymbol{x}_{2}),$ $...,$
$\boldsymbol{x_{n}},f(\boldsymbol{x}_{n})$ be a random sample of input-
outputs, called the training set. Let $f^{(n)}$ be a surrogate model that was
trained on this training set by minimising the Lasso loss function (19). Then
$f$ and $f^{(n)}$ in general do not satisfy assumption (1b).

######  Theorem 2.3 (Lasso in Assumption (1a), with a noisy linear function)

Let $f(\boldsymbol{x})\in
L^{4}\left(\mathbb{R}^{d},\phi(\boldsymbol{x})d\boldsymbol{x}\right)$ be a
function with the following properties:

  1. i)

The input vector contains $d$ uncorrelated random variables
$\boldsymbol{X}=(X_{1},X_{2},...,X_{d})$, with a joint PDF
$\phi(\boldsymbol{x})$, and a covariance matrix $\Sigma=\mathbb{1}_{d}$

  2. ii)

$f$ is a linear function with uncorrelated noise, i.e.

$f(\boldsymbol{x})=\boldsymbol{\alpha}\cdot\boldsymbol{x}+\varepsilon\,,$ (22)

with $\boldsymbol{\alpha}$ the true weights of the function, and $\varepsilon$
a realisation of the random variable $\operatorname*{\mathcal{E}}$. The random
variable $\operatorname*{\mathcal{E}}$ has zero mean and finite variance, and
is independent of $X$.

Let
$\boldsymbol{x}_{1},f(\boldsymbol{x}_{1}),\boldsymbol{x}_{2},f(\boldsymbol{x}_{2}),...,\boldsymbol{x_{n}},f(\boldsymbol{x}_{n})$
be a random sample of input-outputs, called the training set. Let $f^{(n)}$ be
a surrogate model that was trained on the training set by minimising the Lasso
loss function (19), with a regularisation parameter $\lambda$. Then

  1. _a_)

in the asymptotic limit $n\rightarrow\infty$, $f$ and $f^{(n)}$ satisfy
assumption 1b for any regularisation parameter with $|\lambda|<\infty$,

  2. _b_)

in the non-asymptotic limit $n<\infty$, there exists a regularisation
parameter $0<\lambda<\infty$ such that $f$ and $f^{(n)}$ satisfy assumption
(1b).

##  3 Numerical Examples

In the following examples the LMC algorithm (1) was used, with $S=5$. All the
calculations were carried out on single CPU processors of type Intel Xeon Gold
6152 with 2.10 GHz and 16 GB of RAM.

The following plots compare LMC to other methods, by comparing the relative
error of the estimators for a given computational budget $N$ (recall that $N$
is the number of evaluations of $f$, and the remaining costs are assumed to be
negligible). To compute the errors, each method is evaluated $20$ times (or
more), where each time a different seed is used to randomly generate the
datasets. The relative errors are calculated with

$\begin{split}\text{Relative error on
mean}&=\frac{1}{20}\sum_{r=1}^{20}\left|\frac{\mu(N,r)-\operatorname*{\mathbb{E}}[f]}{\operatorname*{\mathbb{E}}[f]}\right|,\\\
\text{Relative error on standard
deviation}&=\frac{1}{20}\sum_{r=1}^{20}\left|\frac{\sqrt{\sigma^{2}(N,r)}-\sqrt{\operatorname*{\mathrm{Var}}[f]}}{\sqrt{\operatorname*{\mathrm{Var}}[f]}}\right|\,,\\\
\end{split}$ (23)

where $\mu(N,r)$ and $\sigma^{2}(N,r)$ are the estimators of a given method,
using $N$ samples and the $r$-th initial seed. The error bars in the plots are
given by the standard deviation of the relative errors (23). In the cases
where the analytical expressions for $\operatorname*{\mathbb{E}}[f]$ and
$\operatorname*{\mathrm{Var}}[f]$ are not known, they are approximated with
the sample estimators (2, 9) and a large sample set.

###  3.1 Linear Function

Let $f$ be a linear function with a large input dimension $d=400$:

$\begin{cases}&f(\boldsymbol{x})=\boldsymbol{\alpha}\cdot\boldsymbol{x},\\\
&\text{with
}\boldsymbol{\alpha}=\left(1,\frac{1}{2},\frac{1}{5},\frac{1}{10},\frac{1}{20},\frac{1}{50},\frac{1}{100},\frac{1}{100},...,\frac{1}{100}\right),\end{cases}$

with $\text{dim}(\boldsymbol{\alpha})=d$ and with a normally distributed input
$X\sim\mathcal{N}(\boldsymbol{0},\mathbb{1}_{d})$.

To perform a thorough benchmark, LMC is compared to all other methods
discussed in this work, with an increasing computational budget $N$. Simple MC
is applied as in eqs. (2, 9), and the Lasso approach consists of a simple
surrogate model approach as in eqs. (4, 11). The LMC is applied as in
algorithm 1 with $S=5$. The Static MFMC method is implemented with the two-
level estimators (6, 12) with $n=0.8\cdot N$, and Adaptive MFMC is applied in
the same way, but the optimal $n$ is chosen by testing several values
$n\in[0.1\cdot N,0.2\cdot N,...,0.9\cdot N]$ and choosing the one with the
smallest estimated MSE. Finally, Biased MFMC uses the full set of $N$ samples
for fitting $f^{(N)}$ and then again for evaluating the two-level estimators.
In all of the methods $M=10^{5}$, this ensures that the leading error in the
two-level estimators comes from the $\operatorname*{\mathcal{O}}(N^{-1})$
term. For this benchmark, the $\lambda$ of the Lasso models is chosen such
that the number of non-zero elements of $\boldsymbol{\beta}$ is $\lfloor
0.95\cdot N\rfloor$, where $\lfloor\,\rfloor$ is the floor function. Each
model was evaluated 30 times with different seeds for generating the random
input points.

Figure 2 shows the comparison between simple MC, LMC, and Biased MFMC, being
used to estimate the mean and standard deviation of $f(X)$. For the estimation
of the mean, the LMC and Biased MFMC models perform similarly, indicating that
the Lasso model has a small bias with respect to the mean. They both have a
smaller MSE than simple MC. For the convergence of the standard deviation we
see three different behaviours:

  * •

the MC approach is unbiased and converges towards the true value, albeit with
a large MSE,

  * •

the Biased MFMC approach is clearly biased. This is expected given that the
training and evaluations sets are the same,

  * •

the LMC approach is unbiased and has a smaller MSE than MC.

Figure 2: Estimation of the mean and standard deviation of $f(X)$, with
respect to the computational budget $N$. The three models were run 30 times
each with different initial seeds. The lines are the means of the 30 runs, and
the coloured-in regions are bounded by the standard deviation of the 30 runs.

Figure 3 shows the mean squared error of all methods. It can be clearly seen
that LMC is the most accurate method for the estimation of the mean and
variance, and that it has an MSE smaller or equal than simple MC for all $N$,
as one would expect from the theory (see thm. 2.1).

Figure 3: Mean squared error on the estimation of the mean and variance of
$f(X)$, with respect to the computational budget $N$. Each line is the mean of
30 runs of a given estimator with randomly generated data. To improve
readability, the $95\%$ confidence intervals have been plotted only for the MC
and LMC methods.

###  3.2 Sobol Function

Let $f$ be the nonlinear function known as the Sobol function, with input
dimension $d=400$:

$\begin{cases}&f(\boldsymbol{x})=\prod_{i=1}^{d}\frac{|4x_{i}-2|+c_{i}}{1+c_{i}},\\\
&\text{with
}\boldsymbol{c}=\left(1,2,5,10,20,50,100,500,500,...,500\right),\end{cases}$
(24)

where the input is sampled from $X_{1},X_{2},...,X_{d}$, i.i.d. uniform random
variables, i.e. $X_{k}\sim U[0,1]\quad\forall k=1,2,...,d$. The Sobol function
is a commonly used benchmark for UQ and sensitivity analysis (e.g. [33, 34,
35]), with unit mean, and variance analytically computed with
$\sigma^{2}=\prod_{i=1}^{d}\Bigg{(}\frac{1}{3(1+c_{i})^{2}}+1\Bigg{)}-1$.

Since the Sobol function (24) is symmetric around $\frac{1}{2}$ (see fig. 4),
fitting a Lasso (or other linear) model to it results in a constant model with
$\boldsymbol{\beta}=(0,0,...,0)$. In such a case the two-level estimators (6,
12) become exactly equal to the simple MC estimators, and no advantage would
be gained from using LMC.

Figure 4: Output of the Sobol function with respect to the first 5 input
dimensions, out of 400 dimensions in total. For these plots 200 sample points
were used.

Nevertheless this problem can be circumvented by adding a nonlinearity to the
Lasso model. This is achieved with a nonlinear transformation $\zeta$ on the
input, prior to fitting the surrogate model, which is equivalent to fitting
the model $f^{(n)}=\boldsymbol{\beta}\cdot\boldsymbol{z}$, with
$\boldsymbol{z}=\zeta(\boldsymbol{x})$. If $\zeta$ is adequately chosen, such
that $\operatorname*{\mathrm{Cov}}[f,X_{k}]\neq 0$ for some $k$, the LMC
estimator will no longer be equal to simple MC.

In practice, finding a good nonlinear transformation requires having some
knowledge of $f$. In this case, the analytical expression of the Sobol
function (24) is known, thus one can intuitively choose the nonlinear
transformation $\zeta(\boldsymbol{x})=|\boldsymbol{x}-0.5|$ where the
subtraction and absolute value operations are applied element-wise to
$\boldsymbol{x}$. If $f$ is not known, one might come up with a similar
transformation by plotting a few samples against the input dimensions, as in
figure 4.

A further remark regarding the Sobol function (24) is that it is not a noisy
linear function of the type (22) discussed in thm. 2.3. Therefore, it is not
guaranteed that the LMC variance estimation will be more accurate than simple
MC, since $f$ and $f^{(n)}$ might not satisfy (1b), as explained in lemma 2.2.
In a real scenario, with unknown $f$, LMC and simple MC are simultaneously
evaluated, the MSE is estimated, and the method with the smallest MSE is used.

Figure 5 compares simple MC to LMC with and without the nonlinear
transformation. In this case LMC was implemented with $M=10^{4}$, and the
regularisation parameter $\lambda$ was chosen with 5-fold cross-validation. As
expected, simple MC and the plain LMC implementation have equal errors. By
adding a nonlinear transformation to Lasso, the errors are reduced and LMC
becomes more accurate than simple MC.

Figure 5: Convergence and relative errors on the estimation of the mean and
standard deviation of the Sobol function with input dimension $d=400$, with
respect to the computational budget $N$. The markers are the mean of 20
evaluations with each method, and the error bars are the standard deviation of
the 20 evaluations. Only a few error bars have been plotted to improve the
readability of the plot.

####  3.2.1 Comparison to PCE

It is of interest to compare LMC to other common UQ methodologies, such as
Polynomial Chaos Expansion (PCE) [33]. PCE is a well-established technique for
UQ and sensitivity analysis, and it is often also benchmarked against the
Sobol function (24). A short description of the PCE method is provided in
appendix B.

The PCE implementation used here is the Chaospy library [36]. To fit the PCE
parameters $\beta_{1},\beta_{2},...,\beta_{P}$ we use Lasso regression (see
section 2.3), and the training set is the same set that is used for the LMC
method, a multivariate uniform distribution. Note that there exist other
methods for fitting the PCE parameters (see [33] and [3] for an overview of
common methods). Here the Lasso method was chosen for computing the
parameters, since the training set is smaller than the number of polynomials,
i.e. $N<P$, and as explained in section 2.3, Lasso can deal with this problem
by finding sparse solutions (see sparse PCE [37, 38]).

Figure 6 shows the comparison between MC, LMC with the nonlinear
transformation as discussed above, and PCE with orders $p=3,4$. The $p=5$
calculations are not plotted since they looked identical to the $p=4$ runs. In
this benchmark the Sobol function of dimension $d=8$ was used, as the Chaospy
implementation of PCE became too computationally expensive for larger input
dimensions. Similar to section 3.1, three distinct behaviours can be observed:

  * •

simple MC is unbiased but has a large errors,

  * •

LMC is unbiased and has smaller errors than any of the other methods,

  * •

PCE, a surrogate-based method, is heavily biased, especially in the estimation
of the variance. Furthermore, we note the different behaviour depending on the
order of the PCE. By increasing from order $p=3$ to $p=4$, the number of
parameters to be fitted increases from $P=165$ to $P=495$ (see table 1). As a
consequence, the higher order model has a smaller bias when $N>165$.

Figure 6: Convergence and relative errors on the estimation of the mean and
standard deviation of the Sobol function with input dimension $d=8$, with
respect to the size of the training set $N$. The results are the mean of 20
evaluations with each method, and the error bars and coloured-in regions are
the standard deviation of the 20 evaluations.

###  3.3 FPUT Lattice Problem

The FPUT lattice problem, proposed by Fermi, Pasta, Ulam, and Tsingou in 1955
[39], models a long chain of $P$ nonlinear oscillators. The FPUT equations of
motion, with boundary conditions $x_{0}(t)=0,\,x_{P+1}(t)=1\,\forall
t\in[0,T]$ and initial conditions
$\dot{\boldsymbol{x}}(t=0)=\boldsymbol{v}_{0},\,\boldsymbol{x}(t=0)=\boldsymbol{x}_{0}$,
can be stated as:

$\ddot{x}_{j}(t)=\frac{k}{m_{j}}\left(x_{j+1}+x_{j-1}-2x_{j}\right)\left[1+\alpha(x_{j+1}-x_{j-1})\right],\quad\forall
j=1,2,\ldots,P.$ (25)

With $\boldsymbol{x}(t)=(x_{1},x_{2},...,x_{P})(t)$ and $x_{j}(t)$ we denote
the position of the $j$-th mass in the chain, and $T$ the final time. By
defining $\mathbcal{l}_{j}=x_{j}-x_{j-1}$ as in figure 7, and redefining the
coupling constant $k^{\prime}_{j}\coloneqq\frac{k}{m_{j}}$, the equations of
motion can be rewritten as

$\ddot{x}_{j}=k^{\prime}_{j}\left(\mathbcal{l}_{j+1}-\mathbcal{l}_{j}\right)+\alpha
k^{\prime}_{j}(\mathbcal{l}^{2}_{j+1}-\mathbcal{l}^{2}_{j})\,,\quad\forall
j=1,2,...,P\,,$

with the boundary and initial conditions from (25).

Figure 7: Schematics of the FPUT model.

For the numerical experiments, we fix $P=40$, set arbitrary initial conditions

$\begin{cases}x_{0,j}&=j\frac{1}{P+1}\,,\\\ v_{0,j}&=\frac{1}{5}\sin\left(3\pi
x_{j}\right)\,,\\\ \end{cases}\quad\forall j=1,2,...,P\,,$

and let the final time be $T=500$ . We now define the UQ problem as estimating
the uncertainty of the final kinetic energy $E_{K}(T)$, given an uncertainty
in the spring constants and the nonlinear parameter $\alpha$. This allows for
a nonlinear problem of arbitrary input dimension $d=P+1$. The UQ problem can
be defined as estimating the mean and variance of the following function

$f\colon\begin{array}[t]{ >{\displaystyle}r >{{}}c<{{}} >{\displaystyle}l
}\mathbb{R}^{P+1}&\to&\mathbb{R}\\\
\left(k^{\prime}_{1},k^{\prime}_{2},...,k^{\prime}_{P},\alpha\right)&\mapsto&f(k^{\prime}_{1},k^{\prime}_{2},...,k^{\prime}_{P},\alpha)=E_{K}(T)\,,\end{array}$
(26)

where $E_{K}(T)$ is the kinetic energy of the system at time $T$, calculated
as

$E_{K}(t)=\frac{1}{2}\sum_{j=1}^{N}m_{j}\dot{x}^{2}_{j}(t)\,,$

and the input arguments of (26) are sampled from a normal distribution

$\mathcal{N}\left(\left(1,1,...,1,\frac{1}{2}\right),\sigma^{2}\right)\,,\quad\text{with
}\sigma=10^{-3}\,.$ (27)

The outputs of (26) are obtained by simulating the FPUT model during time
$t\in[0,T]$, and the final kinetic energy $E_{K}(T)$ is computed from the
velocities in the final timestep. To simulate the FPUT model we integrate the
equations of motion (25) using the Scipy implementation of RK45 [40, 41].

Figure 8 shows the convergence of simple MC and LMC, when used on the UQ
problem (26) with input distribution (27). For this benchmark the LMC
algorithm (alg. 1) was used with $M=10^{4}$. We note that the authors also
carried out this UQ problem with different initial conditions and input
distributions, and that in several cases the LMC method did not provide any
advantage over simple MC, as the two methods showed similar accuracy. The
distribution (27) used in figure 8 was selected for being a case where LMC is
clearly advantageous over simple MC.

Figure 8: Convergence and relative errors on the final kinetic energy of the
FPUT simulation. The results are the mean of 20 evaluations with each method,
and the error bars and coloured-in regions are the standard deviation of the
20 evaluations.

###  3.4 Nuclear Assembly Burnup Calculations

Spent nuclear fuel (SNF) can be hazardous for tens of thousands of years after
its removal from a nuclear reactor, due to decay heat, radiation, and possible
criticality excursions. It is therefore essential to accurately calculate
quantities such as the decay heat to reduce the risks and costs of interim
storage and disposal of SNF (a diagram of the fuel’s life cycle is shown in
figure 9). Moreover, SNF calculations must account for uncertainties, either
with a conservative approach, or with the best estimate plus uncertainty
(BEPU) approach. In the latter case, UQ is included in the SNF calculations,
for example by using Monte Carlo sampling of the input covariances.

Figure 9: Diagram of a simplified nuclear fuel cycle with final repository and
no reprocessing.

The decay heat emitted by spent nuclear fuel is a quantity of interest that
needs to be accurately known, to reduce the risks and costs of interim storage
and final disposal of the spent fuel. Decay heat and its uncertainty can be
numerically estimated with nuclear codes that simulate a nuclear fuel assembly
during its irradiation and cooling phases. An example of such a code is CASMO5
[42], which can be roughly described with

$\texttt{CASMO5}:\underbrace{\begin{pmatrix}\text{Fresh fuel parameters}\\\
\text{Irradiation history}\\\ \text{Reactor parameters}\\\ \text{Nuclear
Data}\end{pmatrix}}_{\text{Uncertain
Input}}\rightarrow\begin{pmatrix}\text{Decay Heat}\\\ \text{Isotopic
Content}\\\ \text{etc...}\end{pmatrix}\,.$ (28)

To calculate the uncertainty of the decay heat, a UQ process is carried out to
propagate the uncertainty of the inputs to the outputs of CASMO5. Generally
speaking, nuclear data is a term used to describe the physical quantities
(e.g. cross sections) that appear in the equations that CASMO5 solves. These
quantities are provided in so-called nuclear data libraries (e.g. [43, 44])
which contain the nominal values of nuclear data and their covariance
(uncertainties plus their correlation matrix). Furthermore, in the literature
it is common practice to treat the uncertainty of nuclear data separately from
other sources of uncertainty [45], and, for propagating nuclear data
uncertainty, simple MC is a widely used and accepted method (e.g. [46, 47,
48]). Nuclear data libraries comprise tens of thousands of parameters, with an
associated covariance that propagates to the output of (28). In this example,
the only input and output of interest are nuclear data and decay heat
respectively, and thus the UQ problem can be described with

$\displaystyle f\colon\mathbb{R}^{15\,557}$ $\displaystyle\to$
$\displaystyle\mathbb{R}$ $\displaystyle\left(\textit{nuclear data}\right)$
$\displaystyle\mapsto$ $\displaystyle f\left(\textit{nuclear
data}\right)=\text{Decay Heat}$

with the input nuclear data sampled from a normal distribution with a
covariance matrix $\Sigma$ provided by the nuclear data library.

The output $f\left(\textit{nuclear data}\right)$ is computed with CASMO5,
where the remaining inputs in (28) are kept constant. To generate the input
samples, the in-house code developed at PSI, SHARK-X [49, 50], is used. With
this tool $15\leavevmode\nobreak\ 557$ nuclear data values are perturbed
according to the nuclear covariance matrix of the ENDF/B-VII.1 library [51],
and provided as input to CASMO5. The perturbed values include the scattering
elastic and inelastic cross-sections, (n,2n) reactions, neutron induced
fission, neutron capture, neutron multiplicity, fission yields, and the
fission spectrum, for 19 discrete energy groups. The modelled representative
fuel assembly comes from the Ringhals-2 pressurised water reactor in Sweden
[52]. It consists of 225 pins of UO2 fuel with $3.095\%$ enrichment, and is
irradiated during four reactor cycles up to a total burnup of 35.7 MWd/kgU.
The output quantity corresponds to the decay heat of the representative fuel
assembly after 500 days of cooling.

The UQ process needs to be carried out for every representative spent fuel
assembly that will be disposed of. If the simple MC approach is used,
thousands of CASMO5 simulations are required for each UQ procedure. The
computational time required for each simulation ranges from 10 minutes to
several hours, and in Switzerland over twelve thousand fuel assemblies [53]
are expected for final disposal. Therefore it is of very high interest to
reduce the computational cost of UQ.

Figure 10 shows the estimation of the mean and variance of the decay heat,
expressed in Watts per tonne of heavy metal, where it can be seen that LMC is
clearly more accurate than simple MC. At $N=1000$, the simple MC method has a
mean relative error of $10^{-3}$ and $10^{-2}$ for the mean and standard
deviation respectively, whereas the LMC method requires less than $N<200$ to
obtain the same relative errors. This represents a reduction in computational
costs by more than a factor of 5. For this benchmark the LMC algorithm (alg.
1) was used with $M=6000$.

Figure 10: Convergence and relative errors on the estimation of the mean and
standard deviation of the decay heat at 500 days of cooling. The results are
the mean of 20 evaluations with each method, and the error bars and coloured-
in regions are the standard deviation of the 20 evaluations. Only the final
error bars were included to improve readability of the plot.

##  4 Conclusion & Outlook

High-dimensional UQ problems are ubiquitous in science, and usually approached
with surrogate-based or simple MC methods. These methods become
computationally expensive in high-dimensional settings, and in such problems
the proposed LMC can play an important role. Moreover, given a set of input-
output data that has been previously used for UQ, LMC can be readily applied
without any extra evaluations of $f$, and immediately provide a better
estimate of the UQ.

From the theory point of view, further research is needed to gain a deeper
understanding on the choice of regularisation parameter $\lambda$, and on the
conditions imposed on $f$ to ensure good accuracy. Indeed, despite numerical
examples showing that LMC has smaller errors than other surrogate-based, MC,
and multifidelity methods, the theory is still lacking in the non-asymptotic
case $N<\infty$.

As further possible investigations, it is worth noting that LMC could in
principle be used with other parametric models a part from Lasso. More complex
models have the potential to better fit the data and further reduce the number
of computationally intensive function evaluations. The only requirement is for
the model to use some form of regularisation, such as the $L^{1}$-norm, which
provides sparsity, a form of model selection, and a small model variance.

## Acknowledgements

The authors would like to thank Swissnuclear for partially sponsoring this
work.

## References

  * [1] J. B. Nagel, J. Rieckermann, and B. Sudret.  Principal component analysis and sparse polynomial chaos expansions for global sensitivity analysis and model calibration: Application to urban drainage simulation.  Reliability Engineering & System Safety, 195:106737, March 2020\.
  * [2] A. Adelmann.  On Nonintrusive Uncertainty Quantification and Surrogate Model Construction in Particle Accelerator Modeling.  SIAM/ASA Journal on Uncertainty Quantification, 7(2):383–416, January 2019.
  * [3] M. Frey and A. Adelmann.  Global sensitivity analysis on numerical solver parameters of Particle-In-Cell models in particle accelerator systems.  Computer Physics Communications, 258:107577, January 2021.
  * [4] M. R. Allen, P. A. Stott, J. F. B. Mitchell, R. Schnur, and T. L. Delworth.  Quantifying the uncertainty in forecasts of anthropogenic climate change.  Nature, 407(6804):617–620, October 2000.
  * [5] V. Rizzoglio, A. Adelmann, A. Gerbershagen, D. Meer, K. P. Nesteruk, and J. M. Schippers.  Uncertainty quantification analysis and optimization for proton therapy beam lines.  Physica Medica, 75:11–18, July 2020.
  * [6] O. Leray, D. Rochman, P. Grimm, H. Ferroukhi, A. Vasiliev, M. Hursin, G. Perret, and A. Pautz.  Nuclear data uncertainty propagation on spent fuel nuclide compositions.  Annals of Nuclear Energy, 94:603–611, August 2016.
  * [7] D. Rochman, A. Vasiliev, H. Ferroukhi, T. Zhu, S. C. van der Marck, and A. J. Koning.  Nuclear data uncertainty for criticality-safety: Monte Carlo vs. linear perturbation.  Annals of Nuclear Energy, 92:150–160, June 2016.
  * [8] D. A. Rochman, A. Vasiliev, A. Dokhane, and H. Ferroukhi.  Uncertainties for Swiss LWR spent nuclear fuels due to nuclear data.  EPJ Nuclear Sciences & Technologies, 4:6, 2018.
  * [9] J. C. Helton, J. D. Johnson, C. J. Sallaberry, and C. B. Storlie.  Survey of sampling-based methods for uncertainty and sensitivity analysis.  Reliability Engineering & System Safety, 91(10):1175–1209, October 2006.
  * [10] M. Arnst and J. Ponthot.  AN OVERVIEW OF NONINTRUSIVE CHARACTERIZATION, PROPAGATION, AND SENSITIVITY ANALYSIS OF UNCERTAINTIES IN COMPUTATIONAL MECHANICS.  International Journal for Uncertainty Quantification, 4(5):387–421, 2014.
  * [11] B. Sudret, S. Marelli, and J. Wiart.  Surrogate models for uncertainty quantification: An overview.  In 2017 11th European Conference on Antennas and Propagation (EUCAP), pages 793–797, March 2017.
  * [12] N. E. Owen, P. Challenor, P. P. Menon, and S. Bennani.  Comparison of Surrogate-Based Uncertainty Quantification Methods for Computationally Expensive Simulators.  SIAM/ASA Journal on Uncertainty Quantification, April 2017.
  * [13] K. Beyer, J. Goldstein, R. Ramakrishnan, and U. Shaft.  When Is “Nearest Neighbor” Meaningful?  In Catriel Beeri and Peter Buneman, editors, Database Theory — ICDT’99, Lecture Notes in Computer Science, pages 217–235, Berlin, Heidelberg, 1999.
  * [14] M. Koppen.  The curse of dimensionality.  5th online world conference on soft computing in industrial applications (WSC5), 1:4–8, 2000.
  * [15] B. Peherstorfer, K. Willcox, and M. Gunzburger.  Survey of Multifidelity Methods in Uncertainty Propagation, Inference, and Optimization.  SIAM Review, 60(3):550–591, January 2018.
  * [16] B. Peherstorfer, K. Willcox, and M. Gunzburger.  Optimal model management for multifidelity monte carlo estimation.  SIAM Journal on Scientific Computing, 38(5):A3163–A3194, 2016.
  * [17] B. Peherstorfer, M. Gunzburger, and K. Willcox.  Convergence analysis of multifidelity Monte Carlo estimation.  Numerische Mathematik, 139(3):683–707, July 2018.
  * [18] E. Sharifnia and S. H. Tindemans.  Multilevel Monte Carlo with Surrogate Models for Resource Adequacy Assessment.  In 2022 17th International Conference on Probabilistic Methods Applied to Power Systems (PMAPS), pages 1–6, June 2022.
  * [19] B. Peherstorfer.  Multifidelity monte carlo estimation with adaptive low-fidelity models.  SIAM/ASA Journal on Uncertainty Quantification, 7(2):579–603, 2019\.
  * [20] R. Tibshirani.  Regression Shrinkage and Selection via the Lasso.  Journal of the Royal Statistical Society. Series B (Methodological), 58(1):267–288, 1996.
  * [21] M. B. Giles.  Multilevel Monte Carlo Path Simulation.  Operations Research, 56(3):607–617, June 2008.
  * [22] A.-L. Haji-Ali, F. Nobile, and R. Tempone.  Multi-index Monte Carlo: when sparsity meets sampling.  Numerische Mathematik, 132(4):767–806, April 2016.
  * [23] C. Bierig and A. Chernov.  Estimation of arbitrary order central statistical moments by the multilevel Monte Carlo method.  Stochastics and Partial Differential Equations Analysis and Computations, 4(1):3–40, March 2016.
  * [24] S. Krumscheid, F. Nobile, and M. Pisaroni.  Quantifying uncertain system outputs via the multilevel Monte Carlo method — Part I: Central moment estimation.  Journal of Computational Physics, 414:109466, August 2020.
  * [25] E. Qian, B. Peherstorfer, D. O’Malley, V. V. Vesselinov, and K. Willcox.  Multifidelity Monte Carlo Estimation of Variance and Sensitivity Indices.  SIAM/ASA Journal on Uncertainty Quantification, 6(2):683–706, January 2018.
  * [26] T. Alsup and B. Peherstorfer.  Context-aware surrogate modeling for balancing approximation and sampling costs in multifidelity importance sampling and bayesian inverse problems.  SIAM/ASA Journal on Uncertainty Quantification, 11(1):285–319, 2023\.
  * [27] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and É. Duchesnay.  Scikit-learn: Machine Learning in Python.  Journal of Machine Learning Research, 12(85):2825–2830, 2011.
  * [28] W. J. Fu.  Penalized Regressions: The Bridge versus the Lasso.  Journal of Computational and Graphical Statistics, 7(3):397–416, 1998.
  * [29] T. T. Wu and K. Lange.  Coordinate descent algorithms for lasso penalized regression.  The Annals of Applied Statistics, 2(1), March 2008.
  * [30] B. Efron and R. J. Tibshirani.  An Introduction to the Bootstrap.  Springer US, Boston, MA, 1993.
  * [31] A. Belloni, V. Chernozhukov, and L. Wang.  Square-root lasso: pivotal recovery of sparse signals via conic programming.  Biometrika, 98(4):791–806, 2011.
  * [32] M. Chichignoud, J. Lederer, and M. J. Wainwright.  A practical scheme and fast algorithm to tune the lasso with optimality guarantees.  The Journal of Machine Learning Research, 17(1):8162–8181, January 2016.
  * [33] B. Sudret.  Global sensitivity analysis using polynomial chaos expansions.  Reliability Engineering & System Safety, 93(7):964–979, July 2008\.
  * [34] K. Konakli and B. Sudret.  Global sensitivity analysis using low-rank tensor approximations.  Reliability Engineering & System Safety, 156:64–83, December 2016\.
  * [35] L. Christos.  Data-driven uncertainty quantification for high-dimensional engineering problems.  PhD thesis, ETH Zurich, November 2019.
  * [36] J. Feinberg and H. P. Langtangen.  Chaospy: An open source tool for designing methods of uncertainty quantification.  Journal of Computational Science, 11:46–57, November 2015.
  * [37] G. Blatman and B. Sudret.  Adaptive sparse polynomial chaos expansion based on least angle regression.  Journal of Computational Physics, 230(6):2345–2367, March 2011\.
  * [38] N. Lüthen, S. Marelli, and B. Sudret.  Sparse Polynomial Chaos Expansions: Literature Survey and Benchmark.  SIAM/ASA Journal on Uncertainty Quantification, 9(2):593–649, January 2021.
  * [39] T. Dauxois, M. Peyrard, and S. Ruffo.  The Fermi–Pasta–Ulam ‘numerical experiment’: history and pedagogical perspectives.  European Journal of Physics, 26(5):S3–S11, July 2005.
  * [40] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, and P. van Mulbregt.  SciPy 1.0: fundamental algorithms for scientific computing in Python.  Nature Methods, 17(3):261–272, March 2020.
  * [41] J. R. Dormand and P. J. Prince.  A family of embedded Runge-Kutta formulae.  Journal of Computational and Applied Mathematics, 6(1):19–26, March 1980.
  * [42] J. Rhodes, K. Smith, and D. Lee.  CASMO-5 development and applications.  Technical report, American Nuclear Society - ANS; La Grange Park (United States), July 2006.
  * [43] A. J. M. Plompen, O. Cabellos, C. De Saint Jean, M. Fleming, A. Algora, M. Angelone, P. Archier, E. Bauge, O. Bersillon, A. Blokhin, F. Cantargi, A. Chebboubi, C. Diez, H. Duarte, E. Dupont, J. Dyrda, B. Erasmus, L. Fiorito, U. Fischer, D. Flammini, D. Foligno, M. R. Gilbert, J. R. Granada, W. Haeck, F.-J. Hambsch, P. Helgesson, S. Hilaire, I. Hill, M. Hursin, R. Ichou, R. Jacqmin, B. Jansky, C. Jouanne, M. A. Kellett, D. H. Kim, H. I. Kim, I. Kodeli, A. J. Koning, A. Yu. Konobeyev, S. Kopecky, B. Kos, A. Krása, L. C. Leal, N. Leclaire, P. Leconte, Y. O. Lee, H. Leeb, O. Litaize, M. Majerle, J. I. Márquez Damián, F. Michel-Sendis, R. W. Mills, B. Morillon, G. Noguère, M. Pecchia, S. Pelloni, P. Pereslavtsev, R. J. Perry, D. Rochman, A. Röhrmoser, P. Romain, P. Romojaro, D. Roubtsov, P. Sauvan, P. Schillebeeckx, K. H. Schmidt, O. Serot, S. Simakov, I. Sirakov, H. Sjöstrand, A. Stankovskiy, J. C. Sublet, P. Tamagno, A. Trkov, S. van der Marck, F. Álvarez Velarde, R. Villari, T. C. Ware, K. Yokoyama, and G. Žerovnik.  The joint evaluated fission and fusion nuclear data library, JEFF-3.3.  The European Physical Journal A, 56(7):181, July 2020.
  * [44] D. A. Brown, M. B. Chadwick, R. Capote, A. C. Kahler, A. Trkov, M. W. Herman, A. A. Sonzogni, Y. Danon, A. D. Carlson, M. Dunn, D. L. Smith, G. M. Hale, G. Arbanas, R. Arcilla, C. R. Bates, B. Beck, B. Becker, F. Brown, R. J. Casperson, J. Conlin, D. E. Cullen, M. A. Descalle, R. Firestone, T. Gaines, K. H. Guber, A. I. Hawari, J. Holmes, T. D. Johnson, T. Kawano, B. C. Kiedrowski, A. J. Koning, S. Kopecky, L. Leal, J. P. Lestone, C. Lubitz, J. I. Márquez Damián, C. M. Mattoon, E. A. McCutchan, S. Mughabghab, P. Navratil, D. Neudecker, G. P. A. Nobre, G. Noguere, M. Paris, M. T. Pigni, A. J. Plompen, B. Pritychenko, V. G. Pronyaev, D. Roubtsov, D. Rochman, P. Romano, P. Schillebeeckx, S. Simakov, M. Sin, I. Sirakov, B. Sleaford, V. Sobes, E. S. Soukhovitskii, I. Stetcu, P. Talou, I. Thompson, S. van der Marck, L. Welser-Sherrill, D. Wiarda, M. White, J. L. Wormald, R. Q. Wright, M. Zerkle, G. Žerovnik, and Y. Zhu.  ENDF/B-VIII.0: The 8th Major Release of the Nuclear Reaction Data Library with CIELO-project Cross Sections, New Standards and Thermal Scattering Data.  Nuclear Data Sheets, 148:1–142, February 2018.
  * [45] A. Vasiliev, J. Herrero, M. Pecchia, D. Rochman, H. Ferroukhi, and S. Caruso.  Preliminary Assessment of Criticality Safety Constraints for Swiss Spent Nuclear Fuel Loading in Disposal Canisters.  Materials, 12(3):494, January 2019.
  * [46] W. Zwermann, A. Aures, L. Gallner, V. Hannstein, B. Krzykacz-hausmann, K. Velkov, and J. S. Martinez.  Nuclear Data Uncertainty and Sensitivity Analysis with XSUSA for Fuel Assembly Depletion Calculations.  Nuclear Engineering and Technology, 46(3):343–352, June 2014.
  * [47] O. Leray, H. Ferroukhi, M. Hursin, A. Vasiliev, and D. Rochman.  Methodology for core analyses with nuclear data uncertainty quantification and application to Swiss PWR operated cycles.  Annals of Nuclear Energy, 110:547–559, December 2017.
  * [48] D. Rochman, A. Dokhane, A. Vasiliev, H. Ferroukhi, and M. Hursin.  Nuclear data uncertainties for Swiss BWR spent nuclear fuel characteristics.  The European Physical Journal Plus, 135(2):233, February 2020.
  * [49] W. Wieselquist, T. Zhu, A. Vasiliev, and H. Ferroukhi.  PSI Methodologies for Nuclear Data Uncertainty Propagation with CASMO-5M and MCNPX: Results for OECD/NEA UAM Benchmark Phase I.  Science and Technology of Nuclear Installations, 2013:1–15, 2013\.
  * [50] A. Aures, F. Bostelmann, M. Hursin, and O. Leray.  Benchmarking and application of the state-of-the-art uncertainty analysis methods XSUSA and SHARK-X.  Annals of Nuclear Energy, 101:262–269, March 2017.
  * [51] M. B. Chadwick, M. Herman, P. Obložinský, M. E. Dunn, Y. Danon, A. C. Kahler, D. L. Smith, B. Pritychenko, G. Arbanas, R. Arcilla, R. Brewer, D. A. Brown, R. Capote, A. D. Carlson, Y. S. Cho, H. Derrien, K. Guber, G. M. Hale, S. Hoblit, S. Holloway, T. D. Johnson, T. Kawano, B. C. Kiedrowski, H. Kim, S. Kunieda, N. M. Larson, L. Leal, J. P. Lestone, R. C. Little, E. A. McCutchan, R. E. MacFarlane, M. MacInnes, C. M. Mattoon, R. D. McKnight, S. F. Mughabghab, G. P. A. Nobre, G. Palmiotti, A. Palumbo, M. T. Pigni, V. G. Pronyaev, R. O. Sayer, A. A. Sonzogni, N. C. Summers, P. Talou, I. J. Thompson, A. Trkov, R. L. Vogt, S. C. van der Marck, A. Wallner, M. C. White, D. Wiarda, and P. G. Young.  ENDF/B-VII.1 Nuclear Data for Science and Technology: Cross Sections, Covariances, Fission Product Yields and Decay Data.  Nuclear Data Sheets, 112(12):2887–2996, December 2011.
  * [52] F. Sturek, L. Agrenius, and O. Osifo.  Measurements of decay heat in spent nuclear fuel at the Swedish interim storage facility, Clab.  Technical Report R-05-62, Svensk Kärnbränslehantering AB, December 2006\.
  * [53] V. Solans, D. Rochman, H. Ferroukhi, A. Vasiliev, and A. Pautz.  Loading optimization for Swiss used nuclear fuel assemblies into final disposal canisters.  Nuclear Engineering and Design, 370:110897, December 2020.
  * [54] R. J. Tibshirani.  The lasso problem and uniqueness.  Electronic Journal of Statistics, 7:1456–1490, January 2013.
  * [55] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani.  Least angle regression.  The Annals of Statistics, 32(2):407–499, April 2004.
  * [56] S. Rosset and J. Zhu.  Piecewise linear regularized solution paths.  The Annals of Statistics, 35(3):1012–1030, July 2007.
  * [57] J. Mairal and B. Yu.  Complexity analysis of the lasso regularization path.  In Proceedings of the 29th International Coference on International Conference on Machine Learning, ICML’12, pages 1835–1842, Madison, WI, USA, June 2012. Omnipress.
  * [58] A. Saltelli and I. Sobol.  Sensitivity Analysis for Nonlinear Mathematical Models. Numerical Experience Sensitivity Analysis for Nonlinear Mathematical Models. Numerical Experience.  Matematicheskoe Modelirovanie, 7(11):16–28, 1995.
  * [59] R. H. Cameron and W. T. Martin.  The Orthogonal Development of Non-Linear Functionals in Series of Fourier-Hermite Functionals.  Annals of Mathematics, 48(2):385–392, 1947.

##  Appendix A Proofs

######  Proof A.1 (Proof of Theorem (2.1))

$\mathcal{M}_{N,M}=\frac{1}{S}\sum_{{\operatorname*{\mathbcal{s}}}=1}^{S}\left(\mu^{(n)}_{N/S.M}\right)_{\operatorname*{\mathbcal{s}}}=\mu^{(n)}_{M}+\frac{1}{S}\sum_{{\operatorname*{\mathbcal{s}}}=1}^{S}\left(\mu_{N/S}-\mu^{(n)}_{N/S}\right)_{\operatorname*{\mathbcal{s}}}=\mu^{(n)}_{M}+\mu_{N}-\mu^{(n)}_{N}=\mu^{(n)}_{N,M}.$

Therefore $\mathcal{M}_{N,M}=\mu^{(n)}_{N,M}$ and they have the same mean
squared error.

$\displaystyle\operatorname*{\text{MSE}}$
$\displaystyle\left(\Sigma^{2}_{N,M},\operatorname*{\mathrm{Var}}[f]\right)=\operatorname*{\text{MSE}}\left(\frac{1}{S}\sum_{\operatorname*{\mathbcal{s}}=1}^{S}\left({\sigma^{2}}^{(n)}_{N/S,M}\right)_{\operatorname*{\mathbcal{s}}},\operatorname*{\mathrm{Var}}[f]\right)$
$\displaystyle=\operatorname*{\mathbb{E}}\left[\left({\operatorname*{{\sigma^{2}}}2}^{(n)}_{M}-\operatorname*{\mathrm{Var}}[f]+\frac{1}{S}\sum_{\operatorname*{\mathbcal{s}}=1}^{S}\left({\operatorname*{{\sigma^{2}}}2}_{N/S}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S}\right)_{\operatorname*{\mathbcal{s}}}\right)^{2}\right]$
$\displaystyle=\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}^{(n)}_{M}\right]+\operatorname*{\mathrm{Var}}\left[\frac{1}{S}\sum_{\operatorname*{\mathbcal{s}}=1}^{S}\left({\operatorname*{{\sigma^{2}}}2}_{N/S}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S}\right)_{\operatorname*{\mathbcal{s}}}\right]+2\underbrace{\operatorname*{\mathrm{Cov}}\left[{\operatorname*{{\sigma^{2}}}2}^{(n)}_{M},\frac{1}{S}\sum_{\operatorname*{\mathbcal{s}}=1}^{S}\left({\operatorname*{{\sigma^{2}}}2}_{N/S}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S}\right)_{\operatorname*{\mathbcal{s}}}\right]}_{=0}$
$\displaystyle=\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}^{(n)}_{M}\right]+\frac{1}{S^{2}}\sum_{\operatorname*{\mathbcal{s}}=1}^{S}\operatorname*{\mathrm{Var}}\left[\left({\operatorname*{{\sigma^{2}}}2}_{N/S}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S}\right)_{\operatorname*{\mathbcal{s}}}\right]+\frac{1}{S^{2}}\sum_{\begin{subarray}{c}\operatorname*{\mathbcal{s}},t=1\\\
\operatorname*{\mathbcal{s}}\neq
t\end{subarray}}^{S}\underbrace{\operatorname*{\mathrm{Cov}}\left[\left({\operatorname*{{\sigma^{2}}}2}_{N/S}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S}\right)_{\operatorname*{\mathbcal{s}}},\left({\operatorname*{{\sigma^{2}}}2}_{N/S}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S}\right)_{t}\right]}_{=0}$
$\displaystyle=\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}^{(n)}_{M}\right]+\frac{S}{S^{2}}\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}_{N/S}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S}\right]=\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}^{(n)}_{M}\right]+\frac{1}{S}\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}_{N/S}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S}\right]\,.$

The last equation is similar to the MSE of a two-level estimator
${\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S,M}$, which is known from eq. (13).
By replacing the MSE in this equation we have

$\displaystyle\operatorname*{\mathrm{Var}}$
$\displaystyle\left[{\operatorname*{{\sigma^{2}}}2}^{(n)}_{M}\right]+\frac{1}{S}\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}_{N/S}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N/S}\right]=\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}^{(n)}_{M}\right]$
$\displaystyle+\frac{1}{S}\frac{1}{\frac{N}{S}}\left(m_{2,2}\left[f+f^{(n)},f-f^{(n)}\right]+\frac{1}{\frac{N}{S}-1}\operatorname*{\mathrm{Var}}[f+f^{(n)}]\operatorname*{\mathrm{Var}}[f-f^{(n)}]-\frac{\frac{N}{S}-2}{\frac{N}{S}-1}\left(\operatorname*{\mathrm{Var}}[f]-\operatorname*{\mathrm{Var}}[f^{(n)}]\right)^{2}\right)$
$\displaystyle=\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}^{(n)}_{M}\right]+\operatorname*{\mathrm{Var}}\left[{\operatorname*{{\sigma^{2}}}2}_{N}-{\operatorname*{{\sigma^{2}}}2}^{(n)}_{N}\right]+\frac{1}{N}\left(\frac{1}{\frac{N}{S}-1}-\frac{1}{N-1}\right)\operatorname*{\mathrm{Var}}[f+f^{(n)}]\operatorname*{\mathrm{Var}}[f-f^{(n)}]$
$\displaystyle+\frac{1}{N}\left(\frac{1}{\frac{N}{S}-1}-\frac{1}{N-1}\right)\left(\operatorname*{\mathrm{Var}}[f]-\operatorname*{\mathrm{Var}}[f^{(n)}]\right)^{2}$
$\displaystyle=\operatorname*{\text{MSE}}\left({\operatorname*{{\sigma^{2}}}2}^{(n)}_{N,M},\operatorname*{\mathrm{Var}}[f]\right)+\mathcal{O}\left(N^{-2}\right)\,.$

######  Proof A.2 (Proof of Lemma (2.1))

The minimum of the loss function is given by the stationary point satisfying

$\frac{\partial}{\partial\beta_{k}}\mathcal{L}(\boldsymbol{\beta})=0\quad\forall\,k=1,2,...,d\,,$

which is a minimum since both terms in (19) are convex. If we write this
equation explicitly, using sub-derivatives since $|\beta_{k}|$ is
nondifferentiable at $\beta_{k}=0$, one finds

$\displaystyle-\sum_{i=1}^{n}x_{ik}\left(f(\boldsymbol{x}_{i})-\boldsymbol{\beta}\cdot\boldsymbol{x}_{i}\right)+\lambda\frac{\partial}{\partial\beta_{k}}|\beta_{k}|=0$
$\displaystyle\iff\sum_{i=1}^{n}f(\boldsymbol{x}_{i})x_{ik}-\sum_{\begin{subarray}{c}l=1\\\
l\neq
k\end{subarray}}^{d}\beta_{l}\sum_{i=1}^{n}x_{il}x_{ik}=\lambda\frac{\partial}{\partial\beta_{k}}|\beta_{k}|+\beta_{k}\sum_{i=1}^{n}x_{ik}^{2}$
$\displaystyle\iff\begin{cases}\frac{1}{\sum_{i=1}^{n}x_{ik}^{2}}\left(\sum_{i=1}^{n}f(\boldsymbol{x}_{i})x_{ik}-\lambda\operatorname*{\text{sgn}}(\beta_{k})-\sum_{\begin{subarray}{c}l=1\\\
l\neq
k\end{subarray}}^{d}\beta_{l}\sum_{i=1}^{n}x_{il}x_{ik}\right)=\beta_{k}&\quad\text{if}\,\beta_{k}\neq
0\,,\\\ &\\\
\sum_{i=1}^{n}f(\boldsymbol{x}_{i})x_{ik}-\sum_{\begin{subarray}{c}l=1\\\
l\neq
k\end{subarray}}^{d}\beta_{l}\sum_{i=1}^{n}x_{il}x_{ik}\in[-\lambda,\lambda]&\quad\text{if}\,\beta_{k}=0\,.\\\
\end{cases}$ (29)

If $\boldsymbol{\beta}=(0,0,...,0)$ is a solution to (29), then
$\lambda\geq\lambda_{max}$, with $\lambda_{max}$ defined in (21). Moreover,
this solution is unique since $\\{x_{ik}\\}$ comes from a continuous
probability distribution [54]. Hence, since
$\boldsymbol{\beta}=\boldsymbol{0}$ is a solution if
$\lambda\geq\lambda_{max}$, it is the only solution and thus
$\lambda\geq\lambda_{max}\Rightarrow\boldsymbol{\beta}=\boldsymbol{0}$.

######  Proof A.3 (Proof of Theorem (2.2) (a))

A weight vector $\boldsymbol{\beta}$ that minimises the Lasso loss (19) must
satisfy (29). If (29) is multiplied on both sides by $\frac{\beta_{k}}{n-1}$
we have

$\beta_{k}\frac{1}{n-1}\sum_{i=1}^{n}f(\boldsymbol{x}_{i})x_{ik}=\beta_{k}\sum_{l=1}^{d}\beta_{l}\frac{1}{n-1}\sum_{i=1}^{n}x_{il}x_{ik}+\frac{\lambda}{n-1}|\beta_{k}|\quad\forall\,k=1,2,...,d\,.$
(30)

In this equation, the left hand side and the first term of the right hand side
contain the unbiased estimator for the covariance, which reads

$\lim_{n\to\infty}\frac{1}{n-1}\sum_{i=1}^{n}\left(z_{i}-\frac{1}{n}\sum_{j=1}^{n}z_{j}\right)\left(y_{i}-\frac{1}{n}\sum_{j=1}^{n}y_{j}\right)=\operatorname*{\mathrm{Cov}}[Z,Y],$

for some random variables $Z$ and $Y$. Then, in the limit
$n\rightarrow\infty$, equation (30) becomes

$\beta_{k}\operatorname*{\mathrm{Cov}}[f,X_{k}]=\sum_{l=1}^{d}\beta_{k}\beta_{l}\operatorname*{\mathrm{Cov}}[X_{l},X_{k}]+\underbrace{\frac{\lambda}{n-1}}_{\rightarrow
0}|\beta_{k}|.$ (31)

Now, by inserting this expression into assumption (1a), the assumption is
shown to be satisfied:

$\begin{split}\operatorname*{\mathrm{Var}}[f-f^{(n)}]&=\operatorname*{\mathrm{Var}}[f]+\operatorname*{\mathrm{Var}}[f^{(n)}]-2\operatorname*{\mathrm{Cov}}(f,f^{(n)})=\operatorname*{\mathrm{Var}}[f]+\operatorname*{\mathrm{Var}}[f^{(n)}]-2\sum_{k=1}^{d}\beta_{k}\operatorname*{\mathrm{Cov}}(f,X_{k})\\\
&=\operatorname*{\mathrm{Var}}[f]+\operatorname*{\mathrm{Var}}[f^{(n)}]-2\sum_{k=1}^{d}\sum_{l=1}^{d}\beta_{k}\beta_{l}\operatorname*{\mathrm{Cov}}(X_{k},X_{l})\\\
&=\operatorname*{\mathrm{Var}}[f]+\operatorname*{\mathrm{Var}}[f^{(n)}]-2\operatorname*{\mathrm{Var}}[f^{(n)}]=\operatorname*{\mathrm{Var}}[f]-\operatorname*{\mathrm{Var}}[f^{(n)}]\leq\operatorname*{\mathrm{Var}}[f].\end{split}$
(32)

######  Proof A.4 (Proof of Theorem (2.2) (b))

For the sake of simplifying notation we define the function

$S(\boldsymbol{\beta})=\operatorname*{\mathrm{Var}}[f-f^{(n)}]-\operatorname*{\mathrm{Var}}[f]\,,$

such that assumption (1a) is satisfied if and only if
$S(\boldsymbol{\beta})\leq 0$. This function can be rewritten as

$S(\boldsymbol{\beta})=\operatorname*{\mathrm{Var}}[f^{(n)}]-2\operatorname*{\mathrm{Cov}}[f,f^{(n)}]=\sum_{k=1}^{d}\sum_{l=1}^{d}\beta_{k}\beta_{l}\operatorname*{\mathrm{Cov}}[X_{k},X_{l}]-2\sum_{k=1}^{d}\beta_{k}\operatorname*{\mathrm{Cov}}[f,X_{k}]\,.$
(33)

Its gradient and Hessian are given by

$\displaystyle\frac{\partial S}{\partial\beta_{m}}(\boldsymbol{\beta})$
$\displaystyle=2\left(\sum_{k=1}^{d}\beta_{k}\operatorname*{\mathrm{Cov}}[X_{k},X_{m}]-\operatorname*{\mathrm{Cov}}[f,X_{m}]\right)\,,$
(34)
$\displaystyle\frac{\partial^{2}S}{\partial\beta_{m}\partial\beta_{n}}(\boldsymbol{\beta})$
$\displaystyle=2\operatorname*{\mathrm{Cov}}[X_{m},X_{n}]=2\,\Sigma_{mn}\,.$
(35)

Since the Hessian (35) is proportional to the covariance matrix $\Sigma\succ
0$, it must be positive definite, and therefore $S(\boldsymbol{\beta})$ is
strictly convex, and thus any local minimum is the global minimum. From (34),
the vector $\boldsymbol{\beta}^{*}$ that minimises $S(\boldsymbol{\beta})$
satisfies

$\operatorname*{\mathrm{Cov}}[f,X_{m}]=\sum_{k=1}^{d}\beta^{*}_{k}\operatorname*{\mathrm{Cov}}[X_{k},X_{m}]\,,$

and it is clear from (32), that $S(\boldsymbol{\beta}^{*})\leq 0$. Therefore
$S$ has a nonpositive global minimum. From the convexity, it follows that
there exists a convex domain $\Omega\subset\mathbb{R}^{d}$, such that

$S(\boldsymbol{\beta})\leq 0\iff\boldsymbol{\beta}\in\Omega\,.$

Additionally we have that the null vector
$\boldsymbol{0}=(0,0,...,0)\in\Omega$, since $S(\boldsymbol{0})=0$.

Now, let $\boldsymbol{\beta}(\lambda)$ be the vector that minimises the loss
function (19), for a given $\lambda>0$. The regularisation path
$\left\\{\boldsymbol{\beta}(\lambda)\,|\,\lambda>0\right\\}$ is continuous and
piecewise linear [55, 56, 57]. Moreover, for $\lambda\rightarrow 0^{+}$ we
have $\boldsymbol{\beta}(\lambda)\rightarrow\boldsymbol{\beta}^{\text{OLS}}$,
and if $\lambda\geq\lambda_{max}$, with $\lambda_{max}$ as in (21), then
$\boldsymbol{\beta}(\lambda)=\boldsymbol{0}$. Therefore
$\boldsymbol{\beta}(\lambda)$ follows a continuous path from
$\boldsymbol{\beta}^{\text{OLS}}$ to $\boldsymbol{0}$ as $\lambda$ increases.
In general, the start of the path is not in $\Omega$, since for a finite
training set with $n<\infty$ it is possible that
$\beta^{\text{OLS}}\notin\Omega$. However, the path ends in $\Omega$ since
$\boldsymbol{0}\in\Omega$. Therefore, there exists at least one point
$\lambda^{*}$ with $0<\lambda^{*}\leq\lambda_{max}$, such that
$\boldsymbol{\beta}(\lambda^{*})\in\Omega$. Thus,
$S(\boldsymbol{\beta}(\lambda^{*}))\leq 0$.

###### Remark A.1

A schematic of a the two-dimensional projection of the domain $\Omega$ is
shown in figure 11. Despite that
$\boldsymbol{\beta}^{\text{OLS}}\notin\Omega$, the Lasso regularisation path
must have some points in $\Omega$, since it is a continuous path that ends at
$\boldsymbol{0}$.

$\beta_{k}$$\beta_{m}$$\boldsymbol{\beta}^{*}$$\boldsymbol{\beta}^{\text{OLS}}$
$\Omega$ Figure 11: Example 2D projection of the domain $\Omega$. The dashed
lines represent three possible regularisation paths.

######  Proof A.5 (Proof of Lemma (2.2))

The left hand side of assumption (1b) can be rewritten as

$\begin{split}&m_{2,2}\left[f+f^{(n)},f-f^{(n)}\right]+\frac{1}{N-1}\operatorname*{\mathrm{Var}}[f+f^{(n)}]\operatorname*{\mathrm{Var}}[f-f^{(n)}]-\frac{N-2}{N-1}\left(\operatorname*{\mathrm{Var}}[f]-\operatorname*{\mathrm{Var}}[f^{(n)}]\right)^{2}\\\
&=m_{4}[f]+m_{4}[f^{(n)}]-2m_{2,2}[f,f^{(n)}]+\frac{1}{N-1}\left({\operatorname*{\mathrm{Var}}}^{2}[f]+{\operatorname*{\mathrm{Var}}}^{2}[f^{(n)}]+2\operatorname*{\mathrm{Var}}[f]\operatorname*{\mathrm{Var}}[f^{(n)}]-4{\operatorname*{\mathrm{Cov}}}^{2}[f,f^{(n)}]\right)\\\
&\quad-\frac{N-2}{N-1}\left({\operatorname*{\mathrm{Var}}}^{2}[f]+{\operatorname*{\mathrm{Var}}}^{2}[f^{(n)}]-2\operatorname*{\mathrm{Var}}[f]\operatorname*{\mathrm{Var}}[f^{(n)}]\right)\\\
&=m_{4}[f]-\frac{N-3}{N-1}{\operatorname*{\mathrm{Var}}}^{2}[f]\\\
&\quad+m_{4}[f^{(n)}]-\frac{N-3}{N-1}{\operatorname*{\mathrm{Var}}}^{2}[f^{(n)}]-2\left(m_{2,2}[f,f^{(n)}]-\operatorname*{\mathrm{Var}}[f]\operatorname*{\mathrm{Var}}[f^{(n)}]\right)-\frac{4}{N-1}{\operatorname*{\mathrm{Cov}}}^{2}[f,f^{(n)}]\,.\\\
\end{split}$

By comparing this expression to the right hand side of (1b), the assumption
can be rewritten as

$\displaystyle m_{4}$
$\displaystyle[f^{(n)}]-\frac{N-3}{N-1}{\operatorname*{\mathrm{Var}}}^{2}[f^{(n)}]-2\left(m_{2,2}[f,f^{(n)}]-\operatorname*{\mathrm{Var}}[f]\operatorname*{\mathrm{Var}}[f^{(n)}]\right)$
$\displaystyle-\frac{4}{N-1}{\operatorname*{\mathrm{Cov}}}^{2}[f,f^{(n)}]\leq
0\quad\forall\,N>2\,,$
$\displaystyle\iff\operatorname*{\mathrm{Var}}[{f^{(n)}}^{2}]-2\operatorname*{\mathrm{Cov}}[f^{2},{f^{(n)}}^{2}]+\frac{2}{N-1}\left({\operatorname*{\mathrm{Var}}}^{2}[f^{(n)}]-2{\operatorname*{\mathrm{Cov}}}^{2}[f,f^{(n)}]\right)\leq
0\quad\forall\,N>2\,,$ (36)

It suffices to show that (1b) is not satisfied in the asymptotic limit
$n\rightarrow\infty$. From eq. (31) we conclude that, in the asymptotic limit,
we have

$\operatorname*{\mathrm{Cov}}[f,f^{(n)}]=\operatorname*{\mathrm{Var}}[f^{(n)}]\,,$

and by injecting this identity into (36) we have the simpler inequality

$\operatorname*{\mathrm{Var}}[{f^{(n)}}^{2}]-2\operatorname*{\mathrm{Cov}}[f^{2},{f^{(n)}}^{2}]-\frac{2}{N-1}{\operatorname*{\mathrm{Var}}}^{2}[f^{(n)}]\leq
0\quad\forall\,N>2\,.$ (37)

Since this condition applies to all $N>2$, we can simply show that it is not
satisfied when $N\rightarrow\infty$. Then assumption (37) becomes

$\operatorname*{\mathrm{Var}}[{f^{(n)}}^{2}]-2\operatorname*{\mathrm{Cov}}[f^{2},{f^{(n)}}^{2}]\leq
0\,.$ (38)

To prove that (38) is in general not satisfied, we provide a one-dimensional
example $f:\mathbb{R}\rightarrow\mathbb{R}$ that does not satisfy (38). The
example is the following: let the input be uniform $X\sim U[-1,1]$, and
$f(x)=\sin(\pi x)$. Then in the limit $n\rightarrow\infty$ we have
$f^{(n)}(x)=\beta x$ with $\beta>0$ and

$\operatorname*{\mathrm{Cov}}[f^{2},{f^{(n)}}^{2}]=\beta^{2}\operatorname*{\mathrm{Cov}}[f^{2},X^{2}]=-\beta^{2}\frac{1}{4\pi^{2}}<0\,,$

and thus (38) is not satisfied.

######  Proof A.6 (Proof of Theorem (2.3) (a))

Start with the asymptotic expression (37), which assumes $n\rightarrow\infty$.
Its last term is negative, and therefore

$\operatorname*{\mathrm{Var}}[{f^{(n)}}^{2}]-2\operatorname*{\mathrm{Cov}}[f^{2},{f^{(n)}}^{2}]-\underbrace{\frac{2}{N-1}{\operatorname*{\mathrm{Var}}}^{2}[{f^{(n)}}]}_{\geq
0}\leq
0\quad\Leftarrow\quad\operatorname*{\mathrm{Var}}[{f^{(n)}}^{2}]-2\operatorname*{\mathrm{Cov}}[f^{2},{f^{(n)}}^{2}]\leq
0\,.$

We can then simply show that the the right hand inequality is true. Since the
explicit expressions for $f$ and ${f^{(n)}}$ are known, we can write

$\displaystyle\operatorname*{\mathrm{Var}}[{f^{(n)}}^{2}]-2\operatorname*{\mathrm{Cov}}[f^{2},{f^{(n)}}^{2}]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[\sum_{k,l,m,p=1}^{d}\beta_{k}\beta_{l}\beta_{m}\beta_{p}\,X_{k}X_{l}X_{m}X_{p}\right]-\sum_{k,l=1}^{d}\beta_{k}^{2}\beta_{l}^{2}$
$\displaystyle\quad\quad-2\operatorname*{\mathrm{Cov}}\left[(\boldsymbol{\alpha}\cdot\boldsymbol{X}+\operatorname*{\mathcal{E}})^{2},(\boldsymbol{\beta}\cdot\boldsymbol{X})^{2}\right]$
$\displaystyle=\sum_{k=1}^{d}\beta_{k}^{4}\left(\operatorname*{\mathbb{E}}[X_{k}^{4}]-1\right)+\sum_{\begin{subarray}{c}k,l=1\\\
k\neq
l\end{subarray}}^{d}\beta_{k}^{2}\beta_{l}^{2}\left(3\operatorname*{\mathbb{E}}[X_{k}^{2}X_{l}^{2}]-1\right)$
$\displaystyle\quad\quad-2\left(\operatorname*{\mathrm{Cov}}[(\boldsymbol{\alpha}\cdot\boldsymbol{X})^{2},(\boldsymbol{\beta}\cdot\boldsymbol{X})^{2}]+\operatorname*{\mathrm{Cov}}[{\operatorname*{\mathcal{E}}}^{2},(\boldsymbol{\beta}\cdot\boldsymbol{X})^{2}]+2\operatorname*{\mathrm{Cov}}[\operatorname*{\mathcal{E}}(\boldsymbol{\alpha}\cdot\boldsymbol{X}),(\boldsymbol{\alpha}\cdot\boldsymbol{X})^{2}]\right)$
$\displaystyle=\sum_{k=1}^{d}\beta_{k}^{4}\operatorname*{\mathrm{Var}}[X_{k}^{2}]+2\sum_{\begin{subarray}{c}k,l=1\\\
k\neq l\end{subarray}}^{d}\beta^{2}_{k}\beta^{2}_{l}$
$\displaystyle\quad\quad-2\left(\sum_{k,l,m,p=1}^{d}\alpha_{k}\alpha_{l}\beta_{m}\beta_{p}\operatorname*{\mathrm{Cov}}[X_{k}X_{l},X_{m}X_{p}]+\sum_{k,l=1}^{d}\beta_{k}\beta_{l}\underbrace{\operatorname*{\mathrm{Cov}}[{\operatorname*{\mathcal{E}}}^{2},X_{k}X_{l}]}_{=0}+2\sum_{k,l,m=1}^{d}\alpha_{k}\beta_{l}\beta_{m}\underbrace{\operatorname*{\mathrm{Cov}}[\operatorname*{\mathcal{E}}X_{k},X_{l}X_{m}]}_{=0}\right)$
$\displaystyle=\sum_{k=1}^{d}\beta_{k}^{4}\operatorname*{\mathrm{Var}}[X_{k}^{2}]+2\sum_{\begin{subarray}{c}k,l=1\\\
k\neq
l\end{subarray}}^{d}\beta^{2}_{k}\beta^{2}_{l}-2\left(\sum_{k=1}^{d}\alpha^{2}_{k}\beta^{2}_{k}\operatorname*{\mathrm{Var}}[X^{2}_{k}]+\sum_{\begin{subarray}{c}k,l=1\\\
k\neq
l\end{subarray}}^{d}\alpha_{k}^{2}\beta_{l}^{2}\underbrace{\operatorname*{\mathrm{Cov}}[X_{k}^{2},X_{l}^{2}]}_{=0}+2\sum_{\begin{subarray}{c}k,l=1\\\
k\neq
l\end{subarray}}^{d}\alpha_{k}\alpha_{l}\beta_{k}\beta_{l}\underbrace{\operatorname*{\mathrm{Cov}}[X_{k}X_{l},X_{k}X_{l}]}_{=\operatorname*{\mathrm{Var}}[X_{k}]\operatorname*{\mathrm{Var}}[X_{l}]=1}\right)$
$\displaystyle=\sum_{k=1}^{d}\left(\beta_{k}^{4}-2\alpha_{k}^{2}\beta_{k}^{2}\right)\operatorname*{\mathrm{Var}}[X_{k}^{2}]+2\sum_{\begin{subarray}{c}k,l=1\\\
k\neq
l\end{subarray}}^{d}\left(\beta^{2}_{k}\beta^{2}_{l}-2\alpha_{k}\beta_{k}\alpha_{l}\beta_{l}\right)\,.$

From (31), in the asymptotic limit $\boldsymbol{\beta}$ satisfies

$\alpha_{k}=\operatorname*{\mathrm{Cov}}[f,X_{k}]=\sum_{l=1}^{d}\beta_{l}\operatorname*{\mathrm{Cov}}[X_{l},X_{k}]=\underbrace{\operatorname*{\mathrm{Var}}[X_{k}]}_{=1}\beta_{k}=\beta_{k}\,,\quad\forall
k=1,2,...,d\,,$

and hence

$\sum_{k=1}^{d}\left(\beta_{k}^{4}-2\alpha_{k}^{2}\beta_{k}^{2}\right)\operatorname*{\mathrm{Var}}[X_{k}^{2}]+2\sum_{\begin{subarray}{c}k,l=1\\\
k\neq
l\end{subarray}}^{d}\left(\beta^{2}_{k}\beta^{2}_{l}-2\alpha_{k}\beta_{k}\alpha_{l}\beta_{l}\right)=-\sum_{k=1}^{d}\alpha^{4}_{k}\operatorname*{\mathrm{Var}}[X_{k}^{2}]-2\sum_{\begin{subarray}{c}k,l=1\\\
k\neq l\end{subarray}}^{d}\alpha^{2}_{k}\alpha^{2}_{l}\leq 0\,.$

######  Proof A.7 (Proof of Theorem (2.3) (b))

Based on inequality(36), define the function

$S(\boldsymbol{\beta})=\operatorname*{\mathrm{Var}}[{f^{(n)}}^{2}]-2\operatorname*{\mathrm{Cov}}[f^{2},{f^{(n)}}^{2}]+\frac{2}{N-1}\left({\operatorname*{\mathrm{Var}}}^{2}[{f^{(n)}}]-2{\operatorname*{\mathrm{Cov}}}^{2}[f,{f^{(n)}}]\right)\,,$

with $N>2$ an integer, such that (1b) is satisfied if and only if
$S(\boldsymbol{\beta})\leq 0$. This function can be rewritten as

$\displaystyle S(\boldsymbol{\beta})$
$\displaystyle=\sum_{k=1}^{d}\left(\beta_{k}^{4}-2\alpha_{k}^{2}\beta_{k}^{2}\right)\operatorname*{\mathrm{Var}}[X_{k}^{2}]+2\sum_{\begin{subarray}{c}k,l=1\\\
k\neq
l\end{subarray}}^{d}\left(\beta^{2}_{k}\beta^{2}_{l}-2\alpha_{k}\beta_{k}\alpha_{l}\beta_{l}\right)$
$\displaystyle+\frac{2}{N-1}\sum_{k,l=1}^{d}\left(\beta^{2}_{k}\beta^{2}_{l}-2\alpha_{k}\alpha_{l}\beta_{k}\beta_{l}\right)$
$\displaystyle=\sum_{k=1}^{d}\left(\beta_{k}^{4}-2\alpha_{k}^{2}\beta_{k}^{2}\right)\left(\operatorname*{\mathrm{Var}}[X_{k}^{2}]+\frac{2}{N-1}\right)+2\frac{N}{N-1}\sum_{\begin{subarray}{c}k,l=1\\\
k\neq
l\end{subarray}}^{d}\left(\beta^{2}_{k}\beta^{2}_{l}-2\alpha_{k}\alpha_{l}\beta_{k}\beta_{l}\right)\,.$

Define the domain $\Omega$ as

$\Omega=\left\\{\boldsymbol{\beta}\,\big{|}\,0\leq|\beta_{k}|\leq|\alpha_{k}|\,\&\,\operatorname*{\text{sgn}}(\beta_{k})=\operatorname*{\text{sgn}}(\alpha_{k})\quad\forall
k=1,2,...,d\right\\}\,.$

Then

$\boldsymbol{\beta}\in\Omega\Rightarrow S(\boldsymbol{\beta})\leq 0\,.$

Now, let $\boldsymbol{\beta}(\lambda)$ be the vector that minimises the loss
function (19), for a given $\lambda>0$. The regularisation path
$\left\\{\boldsymbol{\beta}(\lambda)\,|\,\lambda>0\right\\}$ is continuous and
piecewise linear [55, 56, 57]. Moreover, for $\lambda\rightarrow 0^{+}$ we
have $\boldsymbol{\beta}(\lambda)\rightarrow\boldsymbol{\beta}^{\text{OLS}}$,
and if $\lambda\geq\lambda_{max}$, with $\lambda_{max}$ as in (21), then
$\boldsymbol{\beta}(\lambda)=\boldsymbol{0}$. Therefore
$\boldsymbol{\beta}(\lambda)$ follows a continuous path from
$\boldsymbol{\beta}^{\text{OLS}}$ to $\boldsymbol{0}$ as $\lambda$ increases.
In general, the start of the path is not in $\Omega$, since for a finite
training set with $N_{tr}<\infty$ it is possible that
$\beta^{\text{OLS}}\notin\Omega$. However, the path ends in $\Omega$ since
$\boldsymbol{0}\in\Omega$. Therefore, there exists at least one point
$\lambda^{*}$ with $0<\lambda^{*}\leq\lambda_{max}$, such that
$\boldsymbol{\beta}(\lambda^{*})\in\Omega$. Thus,
$S(\boldsymbol{\beta}(\lambda^{*}))\leq 0$.

###### Remark A.2

A schematic of a the two-dimensional projection of the domain $\Omega$ is
shown in figure 12. Despite that
$\boldsymbol{\beta}^{\text{OLS}}\notin\Omega$, the Lasso regularisation path
must have some points in $\Omega$, since it is a continuous path that ends at
$\boldsymbol{0}$.

$\beta_{k}$$\beta_{m}$$\boldsymbol{\alpha}$$\boldsymbol{\beta}^{\text{OLS}}$
$\Omega$ Figure 12: Example 2D projection of the domain $\Omega$. The dashed
lines represent three possible regularisation paths.

##  Appendix B Polynomial Chaos Expansion

The PCE method is a surrogate based approach to UQ and sensitivity analysis.
With this method a surrogate model ${f^{(n)}}$ is trained, and can then be
used to compute the moments and Sobol indices [58] of ${f^{(n)}}$. The
surrogate model is of the type

${f^{(n)}}(\boldsymbol{x})=\sum_{\boldsymbol{\alpha}=1}^{P}\beta_{\alpha}\Psi_{\boldsymbol{\alpha}}(\boldsymbol{x}),$
(39)

where $\Psi_{\boldsymbol{1}},\Psi_{\boldsymbol{2}},...,\Psi_{\boldsymbol{P}}$
is a set of multivariate polynomials. These multivariate polynomials consist
in the product of univariate polynomials

$\Psi_{\boldsymbol{\alpha}}(\boldsymbol{x})=\prod_{i=1}^{d}\psi_{\alpha_{i}}(x_{i}),\quad\text{with}\quad\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},...,\alpha_{d}),\quad\text{and}\quad\alpha_{i}>0,$

where $\psi_{\alpha_{i}}(x_{i})$ is a univariate polynomial of degree
$\alpha_{i}$. The univariate polynomials are chosen such that they form an
orthogonal basis with respect to the input distribution, such that we have

$\int_{-\infty}^{\infty}\psi_{m}(x)\psi_{n}(x)\phi(x)\,dx=0\iff n\neq m,$

where $\phi(x)$ is the PDF of the input distribution. As an example, for a
Gaussian input distribution one would choose the Hermite polynomials as the
basis, and for a uniform input distribution we would use the Legendre
polynomials. Given a family of orthogonal polynomials, if
$f(\boldsymbol{x})\in
L^{2}\big{(}\mathbb{R}^{d},\phi(\boldsymbol{x})d\boldsymbol{x}\big{)}$ the
following expansion

$f(\boldsymbol{x})=\sum_{\boldsymbol{\alpha}=1}^{\infty}\beta_{\alpha}\Psi_{\boldsymbol{\alpha}}(\boldsymbol{x})$
(40)

is convergent in the $L2$ sense [59]. However, expression (40) needs to be
truncated to have polynomials of a maximum degree $p$. Given a maximum degree
$p$, and an input dimension $d$, the number of polynomials that we can build
is

$P=\begin{pmatrix}d+p\\\ p\\\ \end{pmatrix}=\frac{(d+p)!}{p!d!}.$ (41)

Hence we end up with expression (39), that needs to have its parameters
$\beta_{1},\beta_{2},...,\beta_{P}$ fitted to a set of given data points. With
this model it is then straightforward to compute the mean and variance:

$\begin{split}\mu_{PCE}&=\beta_{1},\\\
\sigma^{2}_{PCE}&=\sum_{\alpha=2}^{P}\beta^{2}_{i}\operatorname*{\mathbb{E}}\left[\Psi_{i}^{2}\right].\end{split}$
(42)

In section 3.2.1, PCE was used with $d=8$ when comparing it to LMC, whereas
the benchmark with LMC alone used $d=400$. Table 1 shows the the number of
trainable parameters in a PCE model of these dimensions.

P | $d=8$ | $d=400$
---|---|---
$p=2$ | 45 | 80601
$p=3$ | 165 | $\sim 10^{7}$
$p=4$ | 495 | $\sim 10^{9}$
$p=5$ | 1287 | $\sim 10^{11}$
Table 1: Number of multivariate polynomials $P$ that can be generated with a
given orthogonal basis, for an input dimension $d$ and maximum degree $p$.