# Adaptive stratified sampling for non-smooth problems

Per Petterssona, Sebastian Krumscheidb
a NORCE Norwegian Research Centre, N-5838 Bergen, Norway
b RWTH Aachen University, 52072 Aachen, Germany

###### Abstract

Science and engineering problems subject to uncertainty are frequently both
computationally expensive and feature non-smooth parameter dependence, making
standard Monte Carlo too slow, and excluding efficient use of accelerated
uncertainty quantification methods relying on strict smoothness assumptions.
To remedy these challenges, we propose an adaptive stratification method
suitable for non-smooth problems and with significantly reduced variance
compared to Monte Carlo sampling. The stratification is iteratively refined
and samples are added sequentially to satisfy an allocation criterion
combining the benefits of proportional and optimal sampling. Theoretical
estimates are provided for the expected performance and probability of failure
to correctly estimate essential statistics. We devise a practical adaptive
stratification method with strata of the same kind of geometrical shapes,
cost-effective refinement satisfying a greedy variance reduction criterion.
Numerical experiments corroborate the theoretical findings and exhibit
speedups of up to three orders of magnitude compared to standard Monte Carlo
sampling.

##  1 Introduction

Many complex engineering problems are non-smooth functions of physical
parameters with unknown values that can be modeled as random variables.
Problems where the solution or its derivative is non-smooth or even
discontinuous in uncertain physical parameters include high-speed flow in
computational fluid dynamics [16], porous media flows [3] and transport [35],
weather and climate predictions [26], and geomechanics [2]. Due to the non-
smooth parameter dependence, quantifying the uncertainty in the solution of
these problems commonly involves repeatedly sampling from the parameter
domain. As a single evaluation of the underlying model is in general very
computationally demanding, efficient sampling is essential to keep the total
number of samples at a minimum to obtain some prescribed statistical tolerance
accuracy. For classic Monte Carlo based sampling methods, for example, the
error decays as $CN^{-1/2}$, where $C$ is typically well approximated by a
constant depending on the problem at hand, and $N$ is the number of samples.
Increasing the number of samples by two orders of magnitude to get a single
order of magnitude error reduction is often numerically intractable, in
particular when the model at hand describes a complex physical problem.
Variance reduction techniques aim at reducing the overall computational cost
by providing an estimator with reduced variance compared to the Monte Carlo
variance $CN^{-1/2}$ for fixed $N$; see, e.g., [18, Ch. 9] or [1, Ch. V] and
the references therein. Improving upon the canonical Monte Carlo rate of
$-1/2$, for example by using quasi-Monte Carlo sampling [24, 9] or spectral
methods [40], typically requires a regular (i.e., differentiable) dependence
of the quantity of interest with respect to the uncertain parameters.
Conversely, general-purpose variance reduction techniques that are effective
when there is a non-smooth parameter dependence aim at decreasing the constant
$C$. Variance reduction via generalized (approximate) control variates, such
as the Multi-Level and Multi-Index Monte Carlo methods, has become a very
popular approach for a wide range of applications due to its computational
efficiency; see, e.g., [4, 12, 15] for a general account and [23, 19, 13, 20]
for variants tailored to, and applications for, particular non-smooth
quantities of interest. These methods are in particular suitable when a
hierarchy of correlated models of different degree of fidelity can be
established, e.g., multiple physical grids of different resolution. For
problems such as those modeling fractured porous media for example, it is
however often prohibitive to introduce conforming grids of different degrees
of refinement.

An alternative means to obtain variance reduction that does not rely on the
concept of a hierarchy of different models (for instance when only a single
unstructured grid is available) is offered by stratified sampling methods [1,
18]. The idea is straightforward: the stochastic domain is partitioned into
disjoint subsets, so-called strata, and a suitable number of samples are drawn
from each stratum. A quantity of interest can be computed as a function of
local mean values of each stratum. The number of samples can be chosen
differently from standard Monte Carlo sampling, which offers the potential to
achieve an estimator with significantly reduced standard error. The optimal
number of samples per stratum is a function of the local variance and the size
of the stratum. Collectively, this leads to a variance reduction compared to
Monte Carlo sampling. Thus, the number of times an often expensive numerical
simulator needs to be solved can be reduced with significant computational
speedup as a result. The existing literature is mainly concerned with either
the proportional or the optimal allocation of samples given a fixed
stratification of the stochastic domain. In this work, we introduce a _hybrid_
sample allocation rule as a linear combination of proportional and optimal
allocation. Moreover, the optimal stratification of the stochastic domain is
not assumed known _apriori_. We therefore devise an original method where we
start from a single stratum and adaptively subdivide the stochastic domain
into finer strata, while assigning new samples to asymptotically satisfy a
sampling distribution defined by some prescribed allocation rule.

Stratified sampling is usually restricted to low dimensions due to exponential
growth in the number of strata with increasing dimension. This problem can be
overcome by using, e.g., latin hypercube sampling where only the marginal
distributions are stratified [18]. Here we use an adaptive strategy to
attenuate the dimensional restriction. It is noteworthy, however, that while
adaptivity enables to use stratified sampling for moderate dimensions, its
efficiency will nonetheless degenerate in higher dimensions. In [11] the
authors introduce an asymptotically optimal stratified sampling estimator for
a fixed stratification using the standard empirical variance estimators to
update the stratum parameters. Adaptivity of the stratification was suggested
in [30], where recursive adaptive stratification was used to bisect a
parallelepiped domain based on the maximum of the squares of the differences
between minimum and maximum sample values in the tentative stratification. An
adaptive stratification algorithm using a sequentially updated stratification
matrix determining stratum boundaries was proposed in [10]. Refined stratified
sampling where single-sample strata are sequentially bisected, and the
stratification enriched by a single-sample stratum at the time was introduced
in [34]. This framework was later combined with hierarchical Latin hypercube
sampling to target very high-dimensional problems [33].

Adaptive stratification is a special case of adaptive decomposition of the
stochastic domain. In [8] the authors describe an adaptive Quasi Monte Carlo
method on a hypercube stratification with bisections of strata maximizing an
error indicator. The work [39] introduces a simplex tessellation of stochastic
parameter space to discretize irregular domains. Targeting discontinuous
functions, the authors of [17] suggested a stochastic domain decomposition
based on solution discontinuities identified via polynomial annihilation. A
Voronoi tessellation of the random space based on random samples, together
with localized surrogate modeling to capture solution discontinuities, has
been proposed in [31]. Also based on polynomial annihilation for discontinuity
detection is the work [14] using support vector machines for building
stochastic surrogate models.

The paper is structured as follows. Section 2 describes the classical
stratified sampling estimator. Moreover, it introduces the novel hybrid
stratified sampling estimator. Section 3 presents the new adaptive
stratification algorithm together with a variance minimizing stratum splitting
strategy as well as a probabilistic robustness analysis due to exact stratum
variances being replaced by sample estimates. In Section 4, we discuss
practical aspects of the proposed adaptive sampling algorithm, including
choices of stratum shapes and efficient strata splitting strategies. Numerical
results are presented in Section 5, where the method is tested on a hierarchy
of different problems, ranging from synthetic test cases to those describing
simplified geomechanics during CO2 injection into an aquifer. The paper ends
with a discussion and outlook in Section 6.

##  2 Stratified Sampling Estimators

We consider a quantity of interest $Q$ that is given as a scalar-valued
measurable function $f\colon\mathbb{R}^{n}\to\mathbb{R}$ of an $n$-dimensional
random vector
$\mathbf{Y}=(Y_{1},\dots,Y_{n})\in\mathfrak{U}\subset\mathbb{R}^{n}$, that is,
$Q=f(\mathbf{Y})$ constitutes a random variable on some probability space
$(\Omega,\mathcal{F},\mathbb{P})$. Specifically, throughout this work we
assume that $Q\in L^{2}(\Omega,\mathcal{F},\mathbb{P})$. Furthermore, we
suppose that the components $Y_{k}$ ($k=1,\ldots,n$) are mutually independent
random variables with finite variance and known, but not necessarily identical
probability distributions. If $F_{k}$ denotes the cumulative distribution
function of $Y_{k}$, $k=1,\ldots,n$, then the random variables $Y_{k}$ and
$F^{-1}_{Y_{k}}(U_{k})$ for a uniformly distributed random variable $U_{k}$ on
$[0,1]$ have the same distribution in view of the inverse probability integral
transform. One may thus write $Y_{k}=F^{-1}_{Y_{k}}(U_{k})$ for $k=1,\ldots,n$
with $U_{1},\ldots,U_{n}$ mutually independent and identically distributed
(i.i.d.), $U_{1}\sim U[0,1]$. Without loss of generality, we may therefore
assume that the stochastic domain $\mathfrak{U}$ is the hypercube
$\mathfrak{U}=[0,1]^{n}$ and that $Y_{k}\sim U[0,1]$ for $k=1,\dots,n$,
viewing the particular choices of $F_{k}$ as part of the “model” $f$, which,
throughout this work, is assumed to be “complicated” and computationally
expensive to evaluate.

###  2.1 Basic definition and properties

Stratified sampling is an estimation technique to approximate
$\mathbb{E}\left(Q\right)=\mathbb{E}\left(f(\mathbf{Y})\right)$ with reduced
variance compared to classic Monte Carlo sampling [1, 18]. The variance
reduction is achieved by dividing the stochastic domain $\mathfrak{U}$ into
multiple disjoint regions, so-called strata, aiming at reducing the variation
in each stratum. Specifically, let $\mathcal{S}$ be a stratification of the
domain $\mathfrak{U}$, in the sense that
$\mathfrak{U}=\cup_{S\in\mathcal{S}}S$ and $S\cap T=\emptyset$ for
$S,T\in\mathcal{S}$ with $S\not=T$. For each stratum $S\in\mathcal{S}$ we
define $Q_{S}\colon\Omega\to\mathbb{R}$ as the random variable with
distribution of $Q=f(\mathbf{Y})$ conditioned upon $\mathbf{Y}\in S$, that is,
the distribution $\mathbb{P}_{Q_{S}}$ of $Q_{S}$ is given by

$A\mapsto\mathbb{P}_{Q_{S}}(A)\equiv\mathbb{P}(Q_{S}\in
A)=\frac{\mathbb{P}(Q\in A,\mathbf{Y}\in S)}{\mathbb{P}(\mathbf{Y}\in
S)}=\frac{\mathbb{P}\bigl{(}f(\mathbf{Y})\in A,\mathbf{Y}\in
S\bigr{)}}{\mathbb{P}(\mathbf{Y}\in S)}\;,$

for any $A\in\mathcal{F}$. Let $p_{S}:=\mathbb{P}(\mathbf{Y}\in S)$ denote the
measure (or “size”) of $S$. Then

$\mathbb{E}\left(Q\right)=\mathbb{E}\left(f(\mathbf{Y})\right)=\sum_{S\in\mathcal{S}}p_{S}\mathbb{E}\left(Q_{S}\right)=\sum_{S\in\mathcal{S}}p_{S}\mathbb{E}\left(f(\mathbf{Y})|\mathbf{Y}\in
S\right)$

by the law of total probability. The stratified sampling estimator $\hat{Q}$
of $\mathbb{E}\left(Q\right)$ is obtained by estimating the expected value
$\mathbb{E}\left(Q_{S}\right)$ in each stratum $S$ by a Monte Carlo average
based on $N_{S}\in\mathbb{N}$ i.i.d. samples of $Q_{S}$:

$\hat{Q}:=\sum_{S\in\mathcal{S}}\frac{p_{S}}{N_{S}}\sum_{j=1}^{N_{S}}Q_{S}^{(j)}\;,\quad
Q_{S}^{(j)}\sim\mathbb{P}_{Q_{S}}\quad 1\leq j\leq N_{S}\;.$ (1)

Consequently, $\hat{Q}$ relies on a total of $N:=\sum_{S\in\mathcal{S}}N_{S}$
number of samples and constitutes an unbiased estimator of
$\mathbb{E}\left(Q\right)$. For (1) to provide a practical estimator, the
measure $p_{S}$ has to be known for every stratum $S\in\mathcal{S}$ and,
furthermore, it needs to be possible to sample from the distribution of
${Q_{S}}$. Throughout this work, we will assume that this is the case; see
Sect. 4 for possibilities how to realize this in practice. The estimator’s
variance is then given by

$\mathbb{V}\left(\hat{Q}\right)=\sum_{S\in\mathcal{S}}\frac{p_{S}^{2}\sigma_{S}^{2}}{N_{S}}=\frac{1}{N}\sum_{S\in\mathcal{S}}\frac{p_{S}^{2}\sigma_{S}^{2}}{N_{S}/N},$
(2)

where
$\sigma_{S}^{2}:=\mathbb{V}\left(Q_{S}\right)=\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in
S\right)$. If the numbers of samples $N_{S}$ are selected such that
$\lim_{N\to\infty}\frac{N}{N_{S}}<\infty$ for all $S\in\mathcal{S}$, then
$\lim_{N\to\infty}N\,\mathbb{V}\left(\hat{Q}\right)<\infty$ and the stratified
sampling estimator $\hat{Q}$ satisfies the central limit theorem, in the sense
that

$\frac{\hat{Q}-\mathbb{E}\left(Q\right)}{\sqrt{\mathbb{V}\left(\hat{Q}\right)}}\underset{N\to\infty}{\Rightarrow}\mathcal{N}(0,1)\;,$
(3)

see, e.g., [1] and the references therein. In practice, when $N$ is
sufficiently large, this asymptotic normality can be used to report also an
approximate confidence region for the point estimate $\hat{Q}$. In that case,
the natural variance estimator for (2) is $\hat{V}$, which is obtained by
replacing the unknown variances $\sigma_{S}^{2}$ in each stratum by the
empirical variances $\hat{\sigma}_{S}^{2}$. For example, let
$z_{(1+p)/2}:=\Phi^{-1}(\frac{1+p}{2})$ be the $(1+p)/2$ quantile of the
standard normal distribution, then

$p\approx\mathbb{P}\left(\left|\hat{Q}-\mathbb{E}\left(Q\right)\right|\leq
z_{\frac{1+p}{2}}{\mathbb{V}\left(\hat{Q}\right)}^{1/2}\right)\;,$ (4)

in view of the asymptotic normality of $\hat{Q}$.

###  2.2 Proportional and optimal sample allocation rules

To make the stratified sampling estimator practical, one has to select a rule
for the number of samples $N_{S}$ in each stratum $S\in\mathcal{S}$. There are
two popular choices, namely proportional allocation and optimal sample
allocation, that is the number of samples are chosen according to

$N^{\text{prop}}_{S}:=p_{S}N\quad\text{and}\quad
N^{\text{opt}}_{S}:=\frac{p_{S}\sigma_{S}}{\sum_{S\in\mathcal{S}}p_{S}\sigma_{S}}N\;,$
(5)

respectively. In practice these rules are, of course, only used up to integer
rounding. For the discussion in this section, we will, however, consider the
number of samples in strata as continuous variables for simplicity.

The latter allocation rule in (5) is optimal in the sense that it provides the
stratification estimator with the smallest variance for a given stratification
$\mathcal{S}$ with a total of $N$ samples. Specifically, the estimator’s
variance ${V}_{\text{opt}}$ using optimal sample allocation
$N_{S}=N_{S}^{\text{opt}}$ and the estimator’s variance ${V}_{\text{prop}}$
using proportional sample allocation $N_{S}=N_{S}^{\text{prop}}$ satisfy,

${V}_{\text{opt}}=\frac{1}{N}{\Bigl{(}\sum_{S\in\mathcal{S}}p_{S}\sigma_{S}\Bigr{)}}^{2}\leq\frac{1}{N}\sum_{S\in\mathcal{S}}p_{S}\sigma_{S}^{2}={V}_{\text{prop}}\leq\frac{\mathbb{V}\left(Q\right)}{N}\;,$
(6)

where the first inequality is a consequence of Jensen’s inequality and the
second inequality follows from the law of total variance. Notice that the
preceding display indicates that stratification with either sample allocation
rule never increases the variance compared to classic Monte Carlo sampling
with $N$ samples.

While the optimal sample allocation rule provides an estimator with minimal
variance, its practical implementation faces the difficulty that the strata
standard deviations $\sigma_{S}$ are typically unknown. Possible remedies
include estimating the standard deviations using a pilot run, see, e.g., [18],
or via an adaptive procedure by sequentially allocating samples, such as the
algorithm introduced in [11] which asymptotically achieves the optimal
variance ${V}_{\text{opt}}$. We reiterate that in this work we are concerned
with the case where $f$ represents a complex model that may be computationally
expensive to evaluate. Pilot runs to estimate the variances in each stratum
are therefore not affordable.

####  2.2.1 Probability of misestimating the estimator’s variance

In addition to potential practical difficulties, the choice of the sample
allocation rule may also affect the accuracy of the stratification estimator’s
variance estimate, which is, for example, used in (4) to offer approximate
confidence intervals. In practice, the natural variance estimators for the
stratification estimator using proportional and optimal sample allocation are

$\hat{V}_{\text{prop}}:=\frac{1}{N}\sum_{S\in\mathcal{S}}p_{S}\hat{\sigma}_{S}^{2}\quad\text{and}\quad\hat{V}_{\text{opt}}:=\frac{1}{N}{\Bigl{(}\sum_{S\in\mathcal{S}}p_{S}\hat{\sigma}_{S}\Bigr{)}}^{2}\;,$
(7)

respectively, where $\hat{\sigma}^{2}_{S}$ denotes the empirical variance

$\hat{\sigma}^{2}_{S}:=\frac{1}{N_{S}-1}\sum_{j=1}^{N_{S}}{\left(Q_{S}^{(j)}-\mkern
1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern
1.5mu_{S}\right)}^{2}\;,\quad\mkern
1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern
1.5mu_{S}:=\frac{1}{N_{S}}\sum_{j=1}^{N_{S}}Q_{S}^{(j)}\;,$

in stratum $S$ based on $N_{S}$ i.i.d. samples
$Q_{S}^{(1)},\dots,Q_{S}^{(N_{S})}$ and
$\hat{\sigma}_{S}\equiv\sqrt{\hat{\sigma}^{2}_{S}}$ is the corresponding
empirical standard deviation. The variance estimates in (7), in particular
$\hat{V}_{\text{opt}}$, may have practical limitations for strata with small
sample sizes $N_{S}$, because $\hat{\sigma}_{S}$ is not an unbiased estimator
for the standard deviation $\sigma_{S}$. In fact, $\hat{\sigma}_{S}$ is only
an asymptotically unbiased estimator in contrast to $\hat{\sigma}^{2}_{S}$,
which is an unbiased estimator for the variance $\sigma_{S}^{2}$. The Lemma
below states concentration inequalities for quantifying the deviations of both
of the stratification estimator’s variance estimates from the theorized
values. Their proof is based on arguments similar to the ones used for classic
Chernoff bounds; see, e.g., [37, Chap. 2].

######  Lemma 1.

For a given stratification $\mathcal{S}$, let $\hat{V}_{\text{prop}}$ and
$\hat{V}_{\text{opt}}$ be as in (7). Suppose that for every stratum
$S\in\mathcal{S}$ there exist finite constants $0\leq M_{S}$, such that the
empirical standard deviations satisfy $\hat{\sigma}_{S}\leq M_{S}$ almost
surely. Then the variance estimator for proportional sample allocation
satisfies

$\mathbb{P}\Bigl{(}\bigl{|}\hat{V}_{\text{prop}}-{V}_{\text{prop}}\bigr{|}\geq\vartheta\Bigr{)}\leq
2\exp{\left(-\frac{2\vartheta^{2}N^{2}}{\sum_{S\in\mathcal{S}}{p_{S}}^{2}{M_{S}}^{4}}\right)}\;,$

for any $\vartheta>0$. Conversely, the variance estimator for optimal sample
allocation satisfies

$\mathbb{P}\Bigl{(}\bigl{|}\hat{V}_{\text{opt}}-{V}_{\text{opt}}\bigr{|}\geq\vartheta\Bigr{)}\leq
2\exp{\left(-\frac{2{\left(\lvert B\rvert-\vartheta
N\right)}^{2}}{{\left(\sum_{S\in\mathcal{S}}p_{S}M_{S}\right)}^{4}}\right)}\;,$

for any $\vartheta>\frac{\lvert B\rvert}{N}$, where
$B:=\sum_{S,T\in\mathcal{S},S\not=T}p_{S}p_{T}\bigl{(}b_{S}b_{T}+b_{S}\sigma_{T}+b_{T}\sigma_{S}\bigr{)}$
denotes the estimator’s bias with
$b_{S}:=\mathbb{E}(\hat{\sigma}_{S})-\sigma_{S}$ for $S\in\mathcal{S}$.

###### Proof.

Let $\vartheta\geq 0$. It follows from Markov’s inequality that

$\displaystyle\mathbb{P}\Bigl{(}\bigl{|}\hat{V}-\mathcal{V}\bigr{|}\geq\vartheta\Bigr{)}$
$\displaystyle=\mathbb{P}\bigl{(}N(\hat{V}-\mathcal{V})\geq\vartheta
N\bigr{)}+\mathbb{P}\bigl{(}N(\hat{V}-\mathcal{V})\leq-\vartheta N\bigr{)}$
(8) $\displaystyle=\mathbb{P}\bigl{(}e^{tN(\hat{V}-\mathcal{V})}\geq
e^{t\vartheta N}\bigr{)}+\mathbb{P}\bigl{(}e^{-tN(\hat{V}-\mathcal{V})}\geq
e^{t\vartheta N}\bigr{)}$ $\displaystyle\leq e^{-t\vartheta
N}\Bigl{[}\mathbb{E}\left(e^{tN(\hat{V}-\mathcal{V})}\right)+\mathbb{E}\left(e^{-tN(\hat{V}-\mathcal{V})}\right)\Bigr{]}\;,$

for any $t>0$, where $V$ and $\hat{V}$ are placeholders for any stratified
sampling estimator’s true variance and its estimated variance, respectively.

We begin with the case of proportional sample allocation. That is, we use (8)
with $\hat{V}\equiv\hat{V}_{\text{prop}}$ and ${V}\equiv{V}_{\text{prop}}$. As
the variance estimators $\hat{\sigma}_{S}^{2}$ and $\hat{\sigma}_{T}^{2}$ are
independent for strata if $S\not=T$, it follows from Hoeffding’s lemma and the
hypotheses that

$\mathbb{E}\left(e^{\pm
tN(\hat{V}-\mathcal{V})}\right)=\prod_{S\in\mathcal{S}}\mathbb{E}\left(e^{\pm
t({\hat{\sigma}}_{S}^{2}-{\sigma_{S}^{2}})p_{S}}\right)\leq\prod_{S\in\mathcal{S}}e^{\frac{t^{2}p_{S}^{2}}{8}M_{S}^{4}}=e^{\frac{t^{2}}{8}\sum_{S\in\mathcal{S}}p_{S}^{2}\,M_{S}^{4}}\;,$

since $\mathbb{E}\left(\hat{\sigma}_{S}^{2}\right)=\sigma_{S}^{2}$ for all
$S\in\mathcal{S}$. Consequently, we obtain the bound

$\mathbb{P}\Bigl{(}\bigl{|}\hat{V}-\mathcal{V}\bigr{|}\geq\vartheta\Bigr{)}\leq
2e^{\left(\frac{t^{2}}{8}\sum_{S\in\mathcal{S}}p_{S}^{2}\,M_{S}^{4}-t\vartheta
N\right)}\;.$

for any $t>0$. Minimizing the right-hand side over $t>0$ yields the claim.

For the case of optimal sample allocation, we proceed similarly. In fact,
using $\hat{V}\equiv\hat{V}_{\text{opt}}$ and ${V}\equiv{V}_{\text{opt}}$,
Hoeffding’s lemma yields

$\mathbb{E}\left(e^{\pm tN(\hat{V}-\mathcal{V})}\right)=\mathbb{E}\left(e^{\pm
t\sum_{S,T\in\mathcal{S}}p_{S}p_{T}(\hat{\sigma}_{S}\hat{\sigma}_{T}-\sigma_{S}\sigma_{T})}\right)\leq
e^{\pm tB+\frac{t^{2}}{8}{\left(\sum_{S\in\mathcal{S}}p_{S}M_{S}\right)}^{4}}$

in this case, where we have used that

$-{\Bigl{(}\sum_{S\in\mathcal{S}}p_{S}\sigma_{S}\Bigr{)}}^{2}\leq
N\bigl{(}\hat{V}-\mathcal{V}\bigr{)}\leq{\Bigl{(}\sum_{S\in\mathcal{S}}p_{S}M_{S}\Bigr{)}}^{2}-{\Bigl{(}\sum_{S\in\mathcal{S}}p_{S}\sigma_{S}\Bigr{)}}^{2}\;.$

Here, the bias term $B$ is given as

$B:=N\mathbb{E}\left(\hat{V}-V\right)=\sum_{S,T\in\mathcal{S},S\not=T}p_{S}p_{T}\bigl{(}b_{S}b_{T}+b_{S}\sigma_{T}+b_{T}\sigma_{S}\bigr{)}\;,$

with $b_{S}:=\mathbb{E}(\hat{\sigma}_{S})-\sigma_{S}$. Combining these bounds
with (8) eventually yields

$\mathbb{P}\Bigl{(}\bigl{|}\hat{V}-{V}\bigr{|}\geq\vartheta\Bigr{)}\leq\min_{t>0}2e^{t\left(\lvert
B\rvert-\vartheta
N\right)+\frac{t^{2}}{8}{\left(\sum_{S\in\mathcal{S}}p_{S}M_{S}\right)}^{4}}=2\exp{\left(-\frac{2{\left(\lvert
B\rvert-\vartheta
N\right)}^{2}}{{\left(\sum_{S\in\mathcal{S}}p_{S}M_{S}\right)}^{4}}\right)}\;,$

provided that $\vartheta>\lvert B\rvert/N$, thus completing the proof. ∎

The Lemma above shows that when using optimal sample allocation, the natural
stratification estimator’s variance estimator is affected by a bias term for
small samples sizes $N$, which originates from the biased strata standard
deviation estimators. In contrast, proportional allocation is not affected by
such a bias. Moreover, a sufficient condition for the hypothesis that the
empirical standard deviations satisfy $\hat{\sigma}_{S}\leq M_{S}<\infty$
almost surely in every stratum $S\in\mathcal{S}$ is that the function $f$ is
bounded on every stratum. The concentration inequalities stated in Lemma 1
indicate that the right-hand sides can be decreased for a stratification
$\mathcal{S}$ that isolates highly varying regions of $f$ in strata of small
measure. The adaptive procedure introduced in this work, which will be
detailed in Sect. 3, takes advantage of this observation. Finally, it is
noteworthy that Lemma 1 may be strengthened by using sharper bounds on the
empirical standard deviations. For the purpose of this work and to highlight
the effects of biased strata standard deviation estimators, the version
presented here is sufficient, however.

###  2.3 Hybrid sample allocation rules

Instead of using either proportional or optimal allocation for the
stratification procedure presented in this work, it will be useful to consider
sample allocation rules that are “between” both versions. Specifically, let
$\alpha\in[0,1]$ and consider the hybrid sample allocation rule

$N_{S}^{\alpha}:=(1-\alpha)N_{S}^{\text{prop}}+\alpha
N_{S}^{\text{opt}}=p_{S}N\left(1+\alpha(\overline{\sigma}_{S}-1)\right)\;,$
(9)

where

$\overline{\sigma}_{S}:=\frac{\sigma_{S}}{\sum_{T\in\mathcal{S}}p_{T}\sigma_{T}}\;,$

so that that number of total samples $N=\sum_{S\in\mathcal{S}}N_{S}^{\alpha}$
is as before. Notice that the hybrid allocation rule contains, in particular,
both proportional allocation ($\alpha=0$) and optimal allocation ($\alpha=1$)
as special cases. For the parameterized sample allocation rule
$N_{S}^{\alpha}$, $\alpha\in[0,1]$, we define the family of hybrid stratified
sampling estimators in the natural way by

$\hat{Q}_{\alpha}:=\sum_{S\in\mathcal{S}}\frac{p_{S}}{N_{S}^{\alpha}}\sum_{j=1}^{N_{S}^{\alpha}}Q_{S}^{(j)}\;,\quad
Q_{S}^{(j)}\sim\mathbb{P}_{Q_{S}}\;,$ (10)

for all $1\leq j\leq N_{S}^{\alpha}$ and $S\in\mathcal{S}$. The variance of
the estimator $\hat{Q}_{\alpha}$ is thus given by

$V_{\alpha}:=\mathbb{V}\left(\hat{Q}_{\alpha}\right)=\frac{1}{N}\sum_{S\in\mathcal{S}}\frac{p_{S}\sigma_{S}^{2}}{1+\alpha(\overline{\sigma}_{S}-1)}\;.$
(11)

An immediate consequence is that the variance of the hybrid stratified
sampling estimator satisfies the bound

$V_{\alpha}\leq\min\left\\{\frac{V_{0}}{1-\alpha},\frac{V_{1}}{\alpha}\right\\}$
(12)

for any $\alpha\in[0,1]$, since
$1+\alpha(\overline{\sigma}_{S}-1)\geq\alpha\overline{\sigma}_{S}$ as well as
$1+\alpha(\overline{\sigma}_{S}-1)\geq 1-\alpha$. Moreover, direct
calculations show that

$V_{0}-V_{\alpha}=\frac{\alpha}{N}\sum_{S\in\mathcal{S}}\frac{p_{S}\sigma_{S}^{2}(\overline{\sigma}_{S}-1)}{1+\alpha(\overline{\sigma}_{S}-1)}\;,\quad\text{and}\quad
V_{\alpha}-V_{1}=\frac{1-\alpha}{N}\sum_{S\in\mathcal{S}}\frac{p_{S}\sigma_{S}^{2}(\overline{\sigma}_{S}-1)}{\overline{\sigma}_{S}\left(1+\alpha(\overline{\sigma}_{S}-1)\right)}\;,$

so that one expects that $V_{\alpha}\approx V_{0}$ and $V_{\alpha}\approx
V_{1}$ for $\alpha\approx 0$ and $\alpha\approx 1$, respectively.

For notational convenience, we introduce the vector notations

$\bm{\sigma}:={(\sigma_{S})}_{S\in\mathcal{S}}\in\mathbb{R}^{\lvert\mathcal{S}\rvert}\quad\text{and}\quad\bm{p}:={(p_{S})}_{S\in\mathcal{S}}\in\mathbb{R}^{\lvert\mathcal{S}\rvert}\;,$

so that
$\overline{\sigma}_{S}=\frac{\sigma_{S}}{\langle\bm{p},\bm{\sigma}\rangle}$,
using the standard Euclidean inner product in
$\mathbb{R}^{\lvert\mathcal{S}\rvert}$. For any $\alpha\in[0,1]$, we define
the variance constant
$C_{\alpha}\colon\mathbb{R}^{\lvert\mathcal{S}\rvert}\to\mathbb{R}$ as a
function of $\bm{\sigma}$ via

$C_{\alpha}(\bm{\sigma}):=\langle\bm{p},\bm{\sigma}\rangle\sum_{S\in\mathcal{S}}\frac{p_{S}\sigma_{S}^{2}}{\alpha\sigma_{S}+(1-\alpha)\langle\bm{p},\bm{\sigma}\rangle}\;,$

so that the hybrid stratification estimator’s variance in (11) can be written
as

$V_{\alpha}=\frac{C_{\alpha}(\bm{\sigma})}{N}\;.$

Moreover, the hybrid stratified sampling estimator $\hat{Q}_{\alpha}$
satisfies the central limit theorem, cf. (3),

$\sqrt{N}\left(\hat{Q}_{\alpha}-\mathbb{E}\left(Q\right)\right)\underset{N\to\infty}{\Rightarrow}\mathcal{N}\left(0,C_{\alpha}(\bm{\sigma})\right)\;,$

for any $\alpha\in[0,1]$, provided that

$\alpha\sigma_{S}+(1-\alpha)\langle\bm{p},\bm{\sigma}\rangle\not=0\quad\forall
S\in\mathcal{S}\;,$ (13)

which is derived from the condition $\lim_{N\to\infty}\frac{N}{N_{S}}<\infty$,
underlying the central limit theorem (3) for general stratified sampling
estimators. Notice that this condition rules out the degenerate extreme case
of a “perfect” stratification $\mathcal{S}$ with $\bm{\sigma}=0$, but also the
case of optimal sample allocation ($\alpha=1$) whenever $\sigma_{S}=0$ for
some $S\in\mathcal{S}$.

####  2.3.1 Asymptotic distribution of the empirical variance constant

As before, to make the central limit theorem above practical, e.g., to provide
confidence estimates, the variance constant $C_{\alpha}(\bm{\sigma})$ needs to
be estimated. The natural estimator based on the empirical standard deviations
$\hat{\bm{\sigma}}:={(\hat{\sigma}_{S})}_{S\in\mathcal{S}}$ thus is
$C_{\alpha}(\hat{\bm{\sigma}})$. For a finite sample size $N$, the hybrid
stratified sampling estimator will, of course, satisfy a similar concentration
inequality to the ones discussed in Sect. 2.2.1 for $\alpha\in\\{0,1\\}$.
Indeed, for $\alpha>0$ there will be a bias due to
$\mathbb{E}(\hat{\sigma}_{S})-\sigma_{S}\not=0$ for $S\in\mathcal{S}$, which
will vanish asymptotically as $N\to\infty$. Complementary to the finite sample
size concentration inequalities presented in Lemma 1, here we discuss the
asymptotic distribution of the empirical variance constant
$C_{\alpha}(\hat{\bm{\sigma}})$ as $N\to\infty$. Specifically, the result
below uses the delta method, see, e.g., [32, 1], which relies on the gradient
$\nabla C_{\alpha}(\bm{\sigma})$ of the variance constant
$C_{\alpha}(\bm{\sigma})$ with respect to $\bm{\sigma}$. We report an explicit
expression of the gradient $\nabla C_{\alpha}(\bm{\sigma})$ in Appendix A for
the reader’s convenience.

######  Lemma 2.

Suppose that $\mathcal{S}$ is such that $Q_{S}\in
L^{4}(\Omega,\mathcal{F},\mathbb{P})$ for all $S\in\mathcal{S}$. Suppose
further that $\alpha\in[0,1]$ is such that condition (13) holds. Then

$\sqrt{N}\left(C_{\alpha}(\hat{\bm{\sigma}})-C_{\alpha}(\bm{\sigma})\right)\underset{N\to\infty}{\Rightarrow}\mathcal{N}\left(0,\varsigma_{\alpha}^{2}\right)\;,\quad\varsigma_{\alpha}^{2}:=\langle\nabla
C_{\alpha}(\bm{\sigma}),\Sigma_{\alpha}\nabla
C_{\alpha}(\bm{\sigma})\rangle\;,$

where $\Sigma_{\alpha}:=\mathrm{diag}\left(\varsigma_{S,\alpha}^{2}\colon
S\in\mathcal{S}\right)$ with

$\varsigma_{S,\alpha}^{2}:=\begin{cases}\frac{\sigma_{S}^{2}\left(\kappa_{S}-1\right)\langle\bm{p},\bm{\sigma}\rangle}{4p_{S}\left((1-\alpha)\langle\bm{p},\bm{\sigma}\rangle+\alpha\sigma_{S}\right)}\;,&\text{if}\quad\sigma_{S}\not=0\;,\\\
0\;,&\text{if}\quad\sigma_{S}=0\;,\end{cases}$

and $\kappa_{S}:=\mathbb{E}\left({\lvert
Q_{S}-\mathbb{E}\left(Q_{S}\right)\rvert}^{4}\right)/\sigma_{S}^{4}$ denoting
the kurtosis of $Q_{S}$ in stratum $S\in\mathcal{S}$.

###### Proof.

The claim is a consequence of the multivariate delta method, e.g., [32, Ch.
3.3] and the fact that the empirical standard deviation in each stratum
satisfies a central limit theorem [7, Ex. 3.6], which is also proved using the
delta method. Specifically, we have

$\sqrt{N_{S}^{\alpha}}\left(\hat{\sigma_{S}}-\sigma_{S}\right)\Rightarrow\mathcal{N}\left(0,\frac{\sigma_{S}^{2}\left(\kappa_{S}-1\right)}{4}\right)$

as $N_{S}^{\alpha}\to\infty$, which is ensured under the hypothesis. Observing
that the different estimators $\hat{\sigma_{S}}$ are independent across strata
completes the proof. ∎

####  2.3.2 Variance reduction for fixed $N$ and uniform Cartesian
stratification

The previous discussion focused on the effects of the hybrid allocation rule
and the sample size $N$ for a given stratification. In this section, we
address the effects of the hybrid allocation rule on the variance reduction
obtained with respect to the regularity of the function $f$ that gives
$Q=f(\mathbf{Y})$. We recall that this work is in particular motivated by
problems for which we expect $f$ to be non-smooth. Below, we study the
variance reduction with respect to the size of a uniform Cartesian
stratification depending on the regularity of $f$, which is an extension of
the work in [1, Chap. V.7]. In particular, we will focus on two special cases
of $f$ that will guide the adaptive stratification procedure in the following.

######  Lemma 3.

Consider a uniform Cartesian stratification $\mathcal{S}$ of
$\mathfrak{U}\subset\mathbb{R}^{n}$. If $f\in C^{1}(\mathfrak{U})$, then the
variance of the hybrid stratification estimator $\hat{Q}_{\alpha}$ satisfies

$V_{\alpha}\leq\frac{nC}{3\,N}\,{\lvert\mathcal{S}\rvert}^{-2/n}\min\left(\frac{1}{\alpha},\frac{1}{1-\alpha}\right)\leq\frac{2nC}{3\,N}\,{\lvert\mathcal{S}\rvert}^{-2/n}\;,$

for any $\alpha\in[0,1]$, where $C=\sup_{\mathbf{u}\in\mathfrak{U}}\|\nabla
f(\mathbf{u})\|_{2}^{2}\leq\|f\|_{C^{1}(\mathfrak{U})}^{2}$.

If $f\colon\mathfrak{U}\to\mathbb{R}$ is a piecewise constant function with a
jump discontinuity of size $\delta>0$ across a curve $\Gamma$ in
$\mathfrak{U}$, then the stratification estimator’s variance is bounded by

$V_{\alpha}\leq\frac{\delta^{2}}{4\,N}\begin{cases}\frac{\lvert\mathcal{T}\rvert}{\lvert\mathcal{S}\rvert}\;,&\alpha=0\;,\\\
\frac{\lvert\mathcal{T}\rvert}{\lvert\mathcal{S}\rvert}\,\min\left\\{\frac{1}{1-\alpha},\frac{1}{\alpha}\,\frac{\lvert\mathcal{T}\rvert}{\lvert\mathcal{S}\rvert}\right\\}\;,&0<\alpha<1\;,\\\
{\left(\frac{\lvert\mathcal{T}\rvert}{\lvert\mathcal{S}\rvert}\right)}^{2}\;,&\alpha=1\;,\end{cases}$

where $\mathcal{T}\subset\mathcal{S}$ denotes the set of all strata that
contain $\Gamma$.

###### Proof.

For a uniform Cartesian stratification, we have
$p_{S}=1/\lvert\mathcal{S}\rvert$ for all $S\in\mathcal{S}$, so that the
stratification estimator’s variance can be written as

$V_{\alpha}=\frac{\|\bm{\sigma}\|_{1}}{N\,\lvert\mathcal{S}\rvert}\sum_{S\in\mathcal{S}}\frac{\sigma_{S}^{2}}{(1-\alpha)\|\bm{\sigma}\|_{1}+\alpha\sigma_{S}\lvert\mathcal{S}\rvert}\;,$

where $\|\bm{\sigma}\|_{1}=\sum_{S\in\mathcal{S}}\sigma_{S}$. We proceed by
bounding the local variances $\sigma_{S}^{2}$ in each stratum
$S\in\mathcal{S}$. To do so, we distinguish the two regularity cases of the
function $f\colon\mathfrak{U}\to\mathbb{R}$.

We begin with the case $f\in C^{1}(\mathfrak{U})$. It follows from a Taylor
expansion of $f$ and the Cauchy–Schwarz inequality that

$\sigma_{S}^{2}=\mathbb{V}\left(Q_{S}\right)=\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in
S\right)\leq\sup_{\mathbf{u}\in\mathfrak{U}}\|\nabla
f(\mathbf{u})\|_{2}^{2}\,\mathbb{E}\left(\|\mathbf{Y}\|_{2}^{2}\right){\lvert\mathcal{S}\rvert}^{-2/n}\;.$

Noting that $\mathbb{E}\left(\|\mathbf{Y}\|_{2}^{2}\right)=n/3$, we obtain the
bound $\sigma_{S}^{2}\leq n\,C\,3^{-1}{\lvert\mathcal{S}\rvert}^{-2/n}$ for
all $S\in\mathcal{S}$, where $C=\sup_{\mathbf{u}\in\mathfrak{U}}\|\nabla
f(\mathbf{u})\|_{2}^{2}$. We thus have
$\|\bm{\sigma}\|_{1}^{2}\leq\frac{n\,C}{3}\,{\lvert\mathcal{S}\rvert}^{2-2/n}$.
Consequently, for $\alpha\in\\{0,1\\}$ the estimator’s variances satisfy

$V_{0}=\frac{1}{N\,\lvert\mathcal{S}\rvert}\sum_{S\in\mathcal{S}}\sigma_{S}^{2}\leq\frac{n\,C}{\,3\,N}{\lvert\mathcal{S}\rvert}^{-2/n}\quad\text{and}\quad
V_{1}=\frac{\|\bm{\sigma}\|_{1}^{2}}{N\,{\lvert\mathcal{S}\rvert}^{2}}\leq\frac{n\,C}{3\,N}\,{\lvert\mathcal{S}\rvert}^{-2/n}\;,$

respectively. For the intermediate values, $0<\alpha<1$ we use inequality (12)
together with the bounds above to complete the first part of the proof.

Next, we consider the case of $f\colon\mathfrak{U}\to\mathbb{R}$ being a
piecewise constant function. In that case, we have $\sigma_{S}^{2}=0$ for all
$S\in\mathcal{S}\setminus\mathcal{T}$. Conversely, for $S\in\mathcal{T}$ the
quantity $Q_{S}$ is a random variable taking two distinct values with
probabilities proportional to the sizes of the stratum subdivision by
$\Gamma$. Consequently, $\sigma_{S}^{2}\leq\delta^{2}/4$, where $\delta>0$ is
the height of the jump discontinuity. For $\alpha\in\\{0,1\\}$ the estimator’s
variances therefore satisfy

$V_{0}=\frac{1}{N\,\lvert\mathcal{S}\rvert}\sum_{S\in\mathcal{S}}\sigma_{S}^{2}\leq\frac{\delta^{2}{\lvert\mathcal{T}\rvert}}{4\,N\,{\lvert\mathcal{S}\rvert}}\quad\text{and}\quad
V_{1}=\frac{\|\bm{\sigma}\|_{1}^{2}}{N\,{\lvert\mathcal{S}\rvert}^{2}}\leq\frac{\delta^{2}{\lvert\mathcal{T}\rvert}^{2}}{4\,N\,{\lvert\mathcal{S}\rvert}^{2}}\;,$

respectively, and from (12) we conclude that

$V_{\alpha}\leq\frac{\delta^{2}}{4\,N}\,\frac{{\lvert\mathcal{T}\rvert}}{\lvert\mathcal{S}\rvert}\min\left\\{\frac{1}{1-\alpha},\frac{1}{\alpha}\,\frac{\lvert\mathcal{T}\rvert}{\lvert\mathcal{S}\rvert}\right\\}$

for any $0<\alpha<1$, which completes the proof. ∎

The preceding lemma indicates that for a differentiable function $f$, optimal
allocation is not superior to proportional allocation, nor to any other hybrid
allocation, asymptotically as the stratification is refined. Conversely, for
piecewise constant functions, optimal allocation is superior to other
allocation rules. However, hybrid allocation with $\alpha\in(0,1]$
sufficiently large offers _quasi-optimal_ asymptotic variance reduction as the
stratification is refined, that is, as
$\gamma=\frac{\lvert\mathcal{T}\rvert}{\lvert\mathcal{S}\rvert}<1$ decreases.
Indeed, if $\alpha\in(0,1]$ satisfies
$\alpha\geq\frac{\gamma_{0}}{1+\gamma_{0}}$, for example $\alpha=\gamma_{0}$,
then $V_{\alpha}\lesssim\gamma^{2}$ for all $\gamma\leq\gamma_{0}$. Finally,
it is also noteworthy that Lemma 3 shows that stratified sampling estimators
face the “curse of dimensionality”, in the sense that the variance reduction
effects deteriorate drastically as the dimension $n$ increases. To mitigate
this effect to some extent, we will use an adaptive procedure for dynamically
creating the stratification tailored to function $f$ that exhiit localized
variability, e.g., discontinuities. This is described next.

##  3 Adaptive Stratification

The optimal allocation of samples in each stratum is not known a priori, since
the stratum standard deviations are not known. Nor is the optimal
stratification of the stochastic domain known in advance, which would minimize
the stratified sampling estimator’s variance; cf. Lemmas 1 and 3. The first
problem will be alleviated similarly to the method proposed in [11], where the
samples are allocated iteratively to fixed strata where the local standard
deviations are computed to satisfy the rules (9) in the large sample size
limit. To remedy the second problem, we propose a method that uses local
variance estimates and a greedy variance reduction approach to determine how
to adaptively split the domain into new strata.

###  3.1 Stratification algorithm

For clarity, we provide a relatively high-level summary of the proposed
adaptive stratification procedure in Algorithm 1, followed by more technical
descriptions of the individual parts in subsequent subsections. We start from
a stratification, which may be based on a priori knowledge of the dependence
of the quantity of interest $Q=f(\mathbf{Y})$. If such information does not
exist, the initial stratification may consist of a few pre-determined strata
or just be a single stratum, i.e., the entire stochastic domain. We
consecutively distribute $N_{\text{new}}$ samples to the current
stratification consisting of $N_{\text{strata}}$ strata, chosen so that the
resulting sample distribution will be as close as possible to the hybrid
allocation rule (9) for a given value of $\alpha$. It is worth noting that
depending on previous sample allocation and due to integer rounding, it may
not be possible to exactly satisfy the theoretical sampling rates.
Furthermore, the current stratification may be quite different from the final
stratification, and the allocation of samples is at best optimal with respect
to the current situation. For each stratum, local properties are updated
(e.g., local means, standard deviations, and number of samples). Next, we use
a greedy approach to refine the stratification by splitting the stratum (or
strata) along some hyperplane that results in the largest variance reduction
within a finite candidate set of refined stratifications. As an additional
splitting criterion, we only consider strata that contain a minimum number of
samples for splitting to make sure that the splitting decision is an informed
one. These two steps are repeated until an upper limit $N_{\text{max}}$ for
the total number of samples has been reached.

One may note that refinement cannot increase the variance of the estimator,
cf. Appendix B, but a poor choice of splitting may require many subsequent
splittings to better adapt to the variability of the quantity of interest
$Q=f(\mathbf{Y})$, and each refinement is only performed after new samples
have been distributed. Hence, for practical reasons of the proposed method,
splitting should be performed with some care. The samples of the split stratum
are distributed to the new strata via (linear) sorting. If the hybrid
parameter $\alpha\in[0,1[$ is set to be dynamic (see Sect. 4.4), then it is
updated before the process of allocating new samples and splitting is repeated
until the limit of the computational budget is reached. Algorithm 1 describes
the workflow in pseudocode and contains references to the subsections where
the method’s components are described in more detail.

Algorithm 1 Adaptive Stratified Sampling.

1:Evaluate initial sample set $\\{Q^{(j)}\\}_{j=1}^{N_{\textup{new}}}$ of size
$N_{\textup{new}}$ using prop. allocation ($\alpha=0$).

2:$N\leftarrow N_{\textup{new}}$.

3:Initialize coarse stratification into $N_{\textup{strata}}$ strata. (_Sect.
4.1_)

4:Determine variance reduction effect from tentative splits.

5:while $N\leq N_{\textup{max}}$ do

6: % _Refine stratification by splitting strata._

7: Find candidate stratum $S_{\textup{split}}$ and hyperplane to split.
(_Sect. 4.2_)

8: if Splitting criteria satisfied then

9: Split stratum $S_{\textup{split}}$ into two (equal size) strata:

10: Redistribute samples in $S_{\textup{split}}$ to the new strata.

11: Compute stratum parameters for the new strata, remove the parent stratum.

12: $N_{\textup{strata}}\leftarrow N_{\textup{strata}}+1$.

13: end

14: % _Add new samples according to current stratification._

15: Determine number of new samples $N_{\textup{new}}$ to allocate so that
$N+N_{\text{new}}\leq N_{\text{max}}$.

16: for $S\in\mathcal{S}$ do

17: Add $N_{\textup{new},S}$ new random samples to stratum $S$. (_Sect. 4.3_)

18: Update stratum parameters. (_Sect. 4.3.1_)

19: Determine variance reduction effect from tentative splits. (_Sect. 3.2_)

20: end

21: Update hybrid allocation parameter $\alpha$. (_Sect. 4.4_)

22: $N\leftarrow N+N_{\textup{new}}$

23:end

24:Output: estimator $\hat{Q}_{\alpha}$ (Eq. (10)).

###  3.2 Variance minimizing splitting

Inaccurate estimates of the stratum standard deviations may lead to lack of
variance reduction of the standard stratification estimator [5]. Likewise, by
symmetry of $\sigma_{S}$ and $p_{S}$ in the case of optimal allocation, c.f.
(5) and (6), an inaccurate estimate of an unknown measures $p_{s}$ may also
result in poor performance. As a remedy, it is natural to choose
stratifications $\mathcal{S}$ with well-defined and easy to compute measures
$p_{S}$ for all $S\in\mathcal{S}$. An additional attractive and thus desirable
feature of an adaptive stratification is that each new stratum that is created
is a subset of only one existing stratum. Otherwise, a new stratum might have
regions that have been sampled with different rates as determined by different
estimates $\hat{\sigma}_{S}$ which would require special treatment in
subsequent allocation of samples as a compensation for non-uniform sampling.
This would be a problem for all hybrid allocation rules with $\alpha>0$. In
this work, we employ stratifications where the stratum measures are easy to
compute exactly, and new strata are created by splitting an existing stratum
into one or more subsets, each of the same geometrical shape as the parent
stratum. In particular, in Sect. 4 we will describe in detail stratifications
based on hyperrectangles and simplicies.

We use a greedy strategy for stratum splitting. The rationale for greedy
strategies is that while there is no way of predicting the final
stratification, the best we can do is to efficiently reduce the variance based
on the information given by the current stratification. The basis for
implementing this strategy is the expression (11) for the variance of the
estimator for a fixed stratification with hybrid sample allocation. We seek
the splitting that leads to the maximum variance reduction when one (or more)
existing stratum is bisected, resulting in a refined stratification. Recalling
the notation introduced in Sect. 2.3, the stratification estimator’s variance
for a given stratification $\mathcal{S}$ and a fixed total number of samples
$N$ can be written as

$\frac{1}{N}\sum_{S\in\mathcal{S}}\frac{p_{S}\sigma_{S}^{2}}{1+\alpha\left(\bar{\sigma}_{S}(\mathcal{S})-1\right)}=:V(\mathcal{S})\;,$

where we have suppressed the dependence on the hybrid parameter $\alpha$,
which is assumed to be fixed in this subsection. Suppose now that stratum
$T\in\mathcal{S}$ is split into two equi-probable, disjoint strata $T^{+}$ and
$T^{-}$ so that $T={T^{-}}\cup{T^{+}}$ and $p_{T}=2p_{T^{\pm}}$. Denote by
$\mathcal{S}_{[T]}$ the refined stratification that is obtained by this
splitting, that is $\mathcal{S}_{[T]}:=(\mathcal{S}\setminus T)\cup T^{+}\cup
T^{-}$. The variance of the corresponding stratification estimator based on
the same total number of samples
$\sum_{S\in\mathcal{S}}N_{S}=N=\sum_{S\in\mathcal{S_{[T]}}}N_{S}$ can thus be
written as

$\displaystyle V(\mathcal{S}_{[T]})$
$\displaystyle=\frac{1}{N}\sum_{S\in\mathcal{S}_{[T]}}\frac{p_{S}\sigma_{S}^{2}}{1+\alpha\left(\bar{\sigma}_{S}(\mathcal{S}_{[T]})-1\right)}$
$\displaystyle=\frac{1}{N}\left(\sum_{S\in\mathcal{S}\setminus
T}\frac{p_{S}\sigma_{S}^{2}}{1+\alpha\left(\bar{\sigma}_{S}(\mathcal{S}_{[T]})-1\right)}+\frac{p_{T}}{2}\sum_{S\in\left\\{T^{+},T^{-}\right\\}}\frac{\sigma_{S}^{2}}{1+\alpha\left(\bar{\sigma}_{S}(\mathcal{S}_{[T]})-1\right)}\right)\;.$

Consequently, the variance reduction obtained by the splitting
$T={T^{-}}\cup{T^{+}}$ compared to not splitting is

$\displaystyle N\left(V(\mathcal{S})-V(\mathcal{S}_{[T]})\right)$
$\displaystyle=p_{T}\left(\frac{\sigma_{T}^{2}}{1+\alpha\left(\bar{\sigma}_{T}(\mathcal{S})-1\right)}-\frac{1}{2}\sum_{S\in\left\\{T^{+},T^{-}\right\\}}\frac{\sigma_{S}^{2}}{1+\alpha\left(\bar{\sigma}_{S}(\mathcal{S}_{[T]})-1\right)}\right)$
$\displaystyle\qquad+\alpha\sum_{S\in\mathcal{S}\setminus
T}\frac{p_{S}\sigma_{S}^{2}\left(\bar{\sigma}_{S}(\mathcal{S}_{[T]})-\bar{\sigma}_{S}(\mathcal{S})\right)}{\left(1+\alpha\left(\bar{\sigma}_{S}(\mathcal{S})-1\right)\right)\left(1+\alpha\left(\bar{\sigma}_{S}(\mathcal{S}_{[T]})-1\right)\right)}\;.$

In the special case of proportional allocation ($\alpha=0$), the expression
simplifies to

$N\left(V(\mathcal{S})-V(\mathcal{S}_{[T]})\right)=p_{T}\left(\sigma_{T}^{2}-\frac{\sigma_{T^{-}}^{2}+\sigma_{T^{+}}^{2}}{2}\right)\;.$

and for optimal allocation,

$N\left(V(\mathcal{S})-V(\mathcal{S}_{[T]})\right)=-p_{T}\Delta_{T}\left(p_{T}\sigma_{T}+2\sum_{S\in\mathcal{S}\setminus
T}p_{S}\sigma_{S}\right)-p_{T}^{2}\Delta_{T}\frac{\sigma_{{T}^{-}}+\sigma_{{T}^{+}}}{2}\;,$

respectively, where the notation
$\Delta_{T}:=\frac{\sigma_{T^{-}}+\sigma_{T^{+}}}{2}-\sigma_{T}$ has been
used.

The _greedy_ splitting strategy then entails finding the stratum
$T\in\mathcal{S}$ that when split across a hyperplane with index $j$ (that
corresponds to some ordering of a set of permissible splitting planes)
provides the largest variance reduction, that is

$(S,j)_{\text{split}}\in\operatorname*{arg\,max}_{T,j}N\left(V(\mathcal{S})-V(\mathcal{S}_{{[T]}_{j}})\right)\;,$
(14)

where $\mathcal{S}_{{[T]}_{j}}$ denotes the refined stratification when the
stratum $T$ is split across the hyperplane $j$ as determined by the shape of
strata. For instance, if the stratification is defined by a Cartesian grid,
the hyperplane will be perpendicular to the coordinate axis $j$.

###  3.3 Probability of failure to identify a stratum that should be split

For any stratum $S\in\mathcal{S}$, recall that $\hat{\sigma}_{S}^{2}$ denotes
the usual empirical variance of $Q_{S}$ based on the $N_{S}\equiv
N_{S}^{\alpha}$ available samples, which is an unbiased estimator of
$\sigma_{S}^{2}\equiv\mathbb{V}\left(Q_{S}\right)$. Notice that
$\sigma_{S}^{2}=0$ implies $\hat{\sigma}_{S}^{2}=0$ and, conversely, that
$\hat{\sigma}_{S}^{2}>0$ implies $\sigma_{S}^{2}>0$. The proposed adaptive
splitting procedure’s efficiency may deteriorate if the empirical variance of
a stratum drastically underestimates the true variance of that stratum, so
that identification and refinement of a high-variability stratum fails. This
can happen, for example, if a stratum does not contain enough samples, so that
a discontinuity is not observed given these samples. To quantify this
“failure” of drastically underestimating the local variance in a stratum $S$,
suppose therefore that $\sigma_{S}^{2}>0$. As this implies that
$\hat{\sigma}_{S}^{2}\geq 0$, it follows from the Paley–Zygmund inequality
combined with the Cauchy–Schwarz inequality that the probability of
underestimating the local variance of $Q_{S}$ is bounded by

$\mathbb{P}(\hat{\sigma}_{S}^{2}\leq\theta\sigma_{S}^{2})\leq
1-\frac{(1-\theta)^{2}}{(1-\theta)^{2}+N_{S}^{-1}\left(\kappa_{S}-\frac{N_{S}-3}{N_{S}-1}\right)}=\frac{\kappa_{S}-\frac{N_{S}-3}{N_{S}-1}}{N_{S}(1-\theta)^{2}+\left(\kappa_{S}-\frac{N_{S}-3}{N_{S}-1}\right)}\;,$
(15)

for any $\theta\in[0,1]$, where $\kappa_{S}$ denotes the kurtosis of $Q_{S}$.
This inequality shows that the probability of underestimating the variance of
$Q_{S}$ by a factor $\theta$ is of order $\mathcal{O}(1/N_{S})$. In fact, if
the kurtosis $\kappa_{S}$ of $Q_{S}$ was known, one could use the upper bound
to determine the required sample size $N_{S}$ to guarantee that
$\mathbb{P}(\hat{\sigma}_{S}^{2}\leq\theta\sigma_{S}^{2})\leq p_{\text{crit}}$
for some prescribed tolerance $p_{\text{crit}}\in(0,1)$ and $\theta\in[0,1]$.
However, the upper bound also reveals that it is affected by the kurtosis
$\kappa_{S}$ in stratum $S$, in the sense that a higher kurtosis leads to a
higher probability of underestimating the variance.

An exemplary borderline case that is insightful in the context of the proposed
stratified sampling procedure for heterogeneous functions $f$ constitutes the
situation where the stratum $S$ contains a region $R\subset S$, $S\setminus
R\not=\emptyset$, on which $f$ is constant, and all available samples are
contained in $R$. That is, $\mathbf{Y}^{(1)},\dots,\mathbf{Y}^{(N)}\in
R\subset S$ so that $Q^{(i)}_{S}=f(\mathbf{Y}^{(i)})\equiv f_{R}$ for all
$i=1,\dots,N_{S}$. It follows that $\hat{\sigma}_{S}^{2}=0$, even though
$\sigma_{S}^{2}>0$, since not the entire stratum has been sampled from and the
region $S\setminus R$ has been missed. Consequently, this can lead to a poor
performance of the adaptive procedure due to missing to identify a high
variance stratum which should have been subdivided. As a matter of fact,
minimizing the risk of falsely estimating a stratum variance to be zero is
particularly important when using Algorithm 1 with optimal sample allocation
($\alpha=1$), since such a drastic underestimation would not be recovered and
thus remain a performance bottleneck.

Although the kurtosis is known for some distributions, the upper bound in (15)
is not of immediate use in practice where the kurtosis is unknown. In fact,
even estimating the kurtosis empirically based on the available samples will
be undefined in that case. A possible remedy is to use a rule-of-thumb, such
as assuming that $Q_{S}$ follows a normal distribution ($\kappa=3$), a uniform
distribution on some interval ($\kappa=9/5$), or some other distribution with
known kurtosis in stratum $S$; see Fig. 1 for examples of the upper bound on
the failure probability in (15) for two different distributions.

Figure 1: Upper bound of the failure probability for a Gaussian distribution
(left) and a uniform distribution (right) based on $N_{S}$ samples.

####  3.3.1 Kernel-based estimation

Alternatively, a sample-based approach with moment approximations by means of
a kernel density estimator (KDE) could be useful. In fact, even in the special
case when obtaining only samples from the region $R$ where $f$ is constant,
this KDE approach is possible and yields, interestingly, the Gaussian kurtosis
rule-of-thumb approximation as we will see below.

In general, the KDE-based moments are given by

$m_{k,\scriptscriptstyle{\textup{KDE}}}(Q_{S})=\frac{1}{\delta
N_{S}}\sum_{i=1}^{N_{S}}\int
q^{k}K\left(\frac{q-Q^{(i)}_{S}}{\delta}\right)\,\textup{d}q\;,$

where $K$ is a symmetric probability density function on $\mathbb{R}$, e.g.,
that of a standard Normal distribution, and $\delta>0$ denotes the bandwidth,
which controls the smoothing [18, Chap. 8.5]. For example, this KDE-moment
approximation with a Gaussian kernel yields

$\displaystyle m_{1,\scriptscriptstyle{\textup{KDE}}}(Q_{S})$
$\displaystyle=\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}Q_{S}^{(i)}\;,$ (16)
$\displaystyle m_{2,\scriptscriptstyle{\textup{KDE}}}(Q_{S})$
$\displaystyle=\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}(Q_{S}^{(i)})^{2}+\delta^{2},$
$\displaystyle m_{3,\scriptscriptstyle{\textup{KDE}}}(Q_{S})$
$\displaystyle=\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}(Q_{S}^{(i)})^{3}+3\delta^{2}\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}Q_{S}^{(i)},$
$\displaystyle m_{4,\scriptscriptstyle{\textup{KDE}}}(Q_{S})$
$\displaystyle=\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}(Q_{S}^{(i)})^{4}+6\delta^{2}\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}(Q_{S}^{(i)})^{2}+3\delta^{4}\;,$

for any $\delta>0$. That is, except for the first moment, the KDE moments
differ from the empirical moments computed directly from samples by a term
$\mathcal{O}(\delta^{2})$, which indicates the additional smoothing.
Consequently, the KDE-based kurtosis approximation
$\kappa_{\scriptscriptstyle{\textup{KDE}}}(Q_{S})$ of $\kappa_{S}$, that is

$\kappa_{\scriptscriptstyle{\textup{KDE}}}(Q_{S}):=\frac{m_{4,\scriptscriptstyle{\textup{KDE}}}(Q_{S})-4m_{1,\scriptscriptstyle{\textup{KDE}}}(Q_{S})m_{3,\scriptscriptstyle{\textup{KDE}}}(Q_{S})+6m_{1,\scriptscriptstyle{\textup{KDE}}}(Q_{S})^{2}m_{2,\scriptscriptstyle{\textup{KDE}}}(Q_{S})-3m_{1,\scriptscriptstyle{\textup{KDE}}}(Q_{S})^{4}}{{\left(m_{2,\scriptscriptstyle{\textup{KDE}}}(Q_{S})-m_{1,\scriptscriptstyle{\textup{KDE}}}(Q_{S})^{2}\right)}^{2}}\;,$

is always well-defined. In fact, even in the borderline case with
$Q_{S}^{(1)}=\dots=Q_{S}^{(N_{S})}=f_{R}$. Indeed, then we find

$\displaystyle m_{1,\scriptscriptstyle{\textup{KDE}}}(Q_{S})$
$\displaystyle=f_{R}\;,\quad
m_{2,\scriptscriptstyle{\textup{KDE}}}(Q_{S})={f_{R}}^{2}+\delta^{2}\;,$
$\displaystyle m_{3,\scriptscriptstyle{\textup{KDE}}}(Q_{S})$
$\displaystyle={f_{R}}^{3}+3\delta^{2}f_{R}\;,\quad
m_{4,\scriptscriptstyle{\textup{KDE}}}(Q_{S})={f_{R}}^{4}+6\delta^{2}{f_{R}}^{2}+3\delta^{4}\;,$

which yields $\kappa_{\scriptscriptstyle{\textup{KDE}}}(Q_{S})=3$ for any
bandwidth $\delta>0$.

##  4 Computational aspects

In the previous section 3, we have described the generic, conceptual aspects
of the adaptive stratified sampling procure introduced in this work. Here, we
will complement the abstract algorithmic component with implementation-
specific details, including the discussion of two concrete classes of
stratifications.

###  4.1 Stratification geometry

For a practical implementation of Algorithm 1, it is desirable that all strata
have the same geometric shape so that their volumes, i.e., probability
measures $p_{S}$ for $S\in\mathcal{S}$, can be easily and exactly computed.
When a stratification is refined, the new strata should therefore also
maintain the same geometric shape as the parent stratum. Since the sampling
rates vary between the strata, any new strata that inherits existing samples
should be properly contained within an existing stratum, i.e., it should be a
member of a partition of a larger stratum.

####  4.1.1 Hyperrectangular stratification

Adaptively splitting the stochastic domain into hyperrectangles provides a
geometry of strata that is easy to visualize (in lower dimensions) and with
properties that have clear generalizations in multiple dimensions. In fact,
both computing the volumes of hyperrectangles and uniform sampling in
hyperrectangles is straightforward. To maintain the number of strata at
feasible levels, we bisect a stratum into two substrata, i.e., split across a
single dimension rather than across all dimensions simultaneously. In
addition, we split only one stratum at the time. We emphasize that other
splitting strategies may be desirable and are possible within the proposed
framework with very minor modifications. Hyperrectangles are expected to be
very efficient for stratifications of problems where sharp features are nearly
parallel with the coordinate axes in the random domain. The numerical cost of
performing the stratification refinement by greedy maximum variance reduction
described in Sect. 3.2 is dominated by sorting the existing samples according
to them being larger or smaller than a tentative stratum boundary, performed
stratum by stratum and dimension by dimension. Note that the samples do not
need to be fully sorted. Hence, the numerical cost assuming splitting after
distribution of $N_{\text{new}}=\sum_{S\in\mathcal{S}}N_{\text{new},S}$ new
samples is proportional to
$\sum_{S\in\mathcal{S}}nN_{\text{new},S}=nN_{\text{new}}$.

####  4.1.2 Simplex stratification

A stratification based on simplices offers more flexibility compared to a
hyperrectangular stratification where all strata boundaries must be aligned
with the coordinate axes in stochastic space. Drawing uniform samples in an
arbitrary simplex can also be done straightforwardly using standard
methodologies, see, e.g., [18, Ch. 3.3.2]. Since the full stochastic domain of
interest $\mathfrak{U}$ is the unit hypercube in $n$ dimensions, it is clearly
possible to form a (trivial) hyperrectangular stratification with just a
single stratum. The smallest number of simplices required to tessellate the
same hypercube is $n!$. Hence, a full $n$-dimensional simplex stratification
is limited to small $n$, but we remark that an approach where simplex
stratification is only applied to a few random dimensions may remedy this
issue if the remaining dimensions are stratified using hyperrectangles, or not
at all.

To perform the greedy variance reduction stratification refinement, each
sample in the input space is transformed to barycentric coordinates, which
involves solving a linear system of size $n\times n$. Alternatively, for $n$
relatively small compared to the number of samples, the transformation can be
carried out with matrix-vector multiplication for each sample. The total cost
is then the inversion of two matrices of size $n\times n$ after splitting, and
$N_{\text{new}}$ matrix-vector multiplications for transformation of each new
sample to barycentric coordinates. Then, for each of the $n(n+1)/2$ tentative
splitting planes of each simplex, the $N_{\text{new}}$ samples are labeled as
belonging to one of two possible substrata, which amounts to comparing the two
barycentric coordinate entries corresponding to the two vertices that would be
assigned to different strata in case of partitioning.

####  4.1.3 Initialization of the simplex stratification

The coarsest simplex stratification of the hypercube
$\mathfrak{U}\subset\mathbb{R}^{n}$, consisting of $n!$ elements, is not
unique. A suboptimal initial stratification may lead to poor variance
reduction even after subsequent dynamic splitting based on refined estimates
of the stratum variances. In particular, if $n$ is large, one may need to
perform numerous simplex splits to compensate for a suboptimal initial
stratification. As a remedy, the stratification is initialized by choosing the
stratification that minimizes the estimator’s sample variance, given an
initial set of Monte Carlo samples in the hypercube. Whereas the task of
stratified sampling is otherwise to add new samples given an existing
stratification, at this initial stage we instead seek an optimal
stratification given a small set of solution samples. To generate the set of
possible simplex tessellations, we proceed as follows. Kuhn’s decomposition is
used to divide the $n$-hypercube into $n$-simplices [21]. An implementation of
Kuhn’s decomposition, which has also been employed in this work, is described
in [6]. Assuming, as before, that the domain is given by the unit hypercube in
$\mathbb{R}^{n}$, Kuhn’s decomposition has the property that $(0,\dots,0)$ and
$(1,\dots,1)$ are common vertices of all simplex elements. The stratification
is thus characterized by the unique edge crossing all $n$ dimensions. We
generate the set of possible simplex stratifications by rotating the Kuhn
decomposition so that we get all combinations of diagonal edges, i.e., the
$2^{n-1}$ ways to choose two vertices
$\mathbf{v}_{\textup{start}},\mathbf{v}_{\textup{end}}$ with the properties
that $\mathbf{v}_{\textup{start},i}=1-\mathbf{v}_{\textup{end},i}$ for
$i=1,\dots,n$. Figure 2 shows the different initial simplex stratifications
for $n=3$.

Figure 2: Initial simplex stratifications of the hypercube in $n=3$ based on
Kuhn’s decomposition.

For each candidate stratification, the existing initial Monte Carlo samples
are distributed to the simplices through a transformation to barycentric
coordinates, and local statistics are computed. The minimum estimator variance
stratification among the $2^{n}$ candidates is chosen as the initial
stratification of the hypercube.

###  4.2 Identification of bisection planes for splitting

A stratum can theoretically be split along any one of an uncountable number of
planes, but finding a nearly optimal plane may require an unnecessary amount
of work, especially since we expect more information through new samples to
change the shape of the optimal stratification. During stratification
refinement, strata are instead split to obtain the greatest variance reduction
among a finite number of possibilities. To limit the search space of possible
new strata, we consider only all possible bisections into two new members of
the stratum class (here: hyperrectangles or simplices) of all strata and
choose the minimizer as in (14). Next we describe how this is done in the case
of hyperrectangular and simplex stratifications, but note that other options
such as splitting more than one stratum at once could also be implemented
without major changes to the algorithm.

With a hyper-rectangular stratification, splitting must remain aligned with
the coordinate system. Each time we split a stratum across
$n_{\textup{split}}$ dimensions, the number of strata is increased by
$2^{n_{\textup{split}}}-1$. This is not a severe limitation on the choice of
$n_{\textup{split}}$, but we do require at least one sample in each stratum,
so $n_{\textup{split}}$ should not be too large. With $n_{\textup{split}}=1$,
which was also used for the numerical experiments discussed in Sect. 5, there
are $n$ possible ways to split a hyperrectangle by bisection into two new
hyperrectangles, as illustrated for $n=3$ in Fig. 3.

56218734

56218734

56218734

Figure 3:  Splitting of 3D hyperrectangle by bisection parallel to the
coordinate axes.

Simplices can be split in multiple ways to result in a refined simplex
partitioning. Perhaps the simplest way to avoid excessive growth in the number
of strata is to split a candidate simplex in two equal-sized simplices defined
by adding the hyperplane going through the mid-point of one edge and the
remaining vertices that are not the end points of the edge whose midpoint is
used. As a simplex has $n+1$ edges, and each one corresponds to a possible
bisection, there are $(n+1)n/2$ ways to split a simplex using this method.
This is illustrated for $n=3$ in Fig. 4, with six different ways to bisect the
simplex, each producing two new simplices. Clearly, a larger number of
conditional variances for each tentative bisection must be computed, compared
to the case of hyperrectangular dynamic stratification. Bisection leads to
very simple computation of the probability measure of new strata: they are
just 1/2 of the measure of their parent stratum. A straightforward
generalization of the above described bisection method is to split along a
plane where the split point of the edge is not a midpoint. In that case, the
volume can be computed by
$1/n!\left|\det(\mathbf{v}_{2}-\mathbf{v}_{1},\mathbf{v}_{3}-\mathbf{v}_{1},\dots,\mathbf{v}_{n+1}-\mathbf{v}_{1})\right|$,
where $\\{\mathbf{v}_{i}\\}_{i=1}^{n+1}$ is the set of vertices.

3210

3210

3210

3210

3210

3210

Figure 4:  Bi-section of 3D simplex along the hyperplane that intersects the
midpoint of an edge and all vertices that do not belong to that edge.

###  4.3 Sequential allocation of samples

The stratification is dynamically updated by bisection (or other means of
splitting) of existing strata. We refer to a cycle of bisection and subsequent
sampling as an _iteration_ of Algorithm 1. At every iteration and for a given
value $\alpha$, $N_{\text{new}}=\sum_{S}N_{\text{new},S}$ samples are newly
distributed and added to the current stratification for an asymptotic (i.e.,
for $N_{\max}$ sufficiently large) sample allocation according to the rates

$q_{S}^{\alpha}:=(1-\alpha)p_{S}+\alpha\bar{\sigma}_{S}p_{S}\;.$

As the hybrid parameter $\alpha\in[0,1]$ is fixed, we will suppress the
explicit dependence and simply write $q_{S}\equiv q_{S}^{\alpha}$. Adding
these $N_{\text{new}}$ samples at the current iteration to the current
stratification then amounts to adding $N_{\text{new},S}$ new samples to
stratum $S$ satisfying

$N_{\text{new},S}=\max(0,\min(\lceil(N_{\text{total}}+N_{\text{new}})q_{S}-N_{S}\rceil,N_{\text{new}})),$

i.e., striving to fulfill $N_{S}=q_{S}(N_{\text{total}}+N_{\text{new}})$ as in
(9), assuming a fixed stratification. Here $\lceil x\rceil$ denotes the
smallest integer greater than or equal to $x$. Even for hybrid sampling rules
where a fraction of samples is proportionally allocated (i.e., $\alpha<1$),
some strata may not be assigned new samples in every iteration. Eventually,
all strata will receive new samples, but for strata where variability has been
underestimated this can require several iterations and may not happen before
the total sampling budget of $N_{\text{max}}$ samples has been reached. To
avoid that a stratum where the standard deviation has mistakenly been assigned
to be zero does not get updated, we may reserve at least one sample for every
stratum, independent of the sampling rates. Hence, we may alternatively use
the allocation rule

$N_{\text{new},S}=1+\max(0,\min(\lceil(N_{\text{total}}+N_{\text{new}}-N_{\text{strata}})q_{S}-N_{S}\rceil,N_{\text{new}}-N_{\text{strata}})).$

In case at least one stratum has been oversampled, e.g., as a result of a
stratum being created from a parent stratum with larger variance, both
sampling rules require a total number of samples exceeding $N_{\text{new}}$.
In these cases, we add samples according to the allocation rule as long as the
total sampling budget is not exceeded.

Since the stratification is dynamic and changes with (almost) every iteration,
there is a trade-off between choosing $N_{\text{new}}$ large enough to
approximately satisfy the rates $q_{S}$, and small enough to save samples for
further refinement of the stratification. A reasonable compromise may be to
always add a constant $c$ times the current number of strata so that every
stratum gets _on average_ $c$ new samples in each iteration,

$N_{\text{new}}=cN_{\text{strata}}.$ (17)

We emphasize that the distribution of samples is not restricted to the choice
(17), and other choices may significantly impact the performance of the
algorithm. At the same time, the results in Sect. 3.3 offer a convenient
probabilistic interpretation of the value for $c$ in (17) as the
(stratification averaged) likelihood of drastically underestimating a local
stratum variance.

####  4.3.1 Updating stratum statistics

Updating the statistics of each stratum in every iteration may become quite
costly for large sample sizes, in particular if only a small number of samples
are added at every iteration so that many iterations are performed. Instead of
recomputing the sample means and standard deviations from scratch using the
standard formulae, we use the following update formulae

$\displaystyle\hat{\mu}_{S}$
$\displaystyle=\frac{N_{\text{old},S}\hat{\mu}_{\text{old},S}+N_{\text{new},S}\hat{\mu}_{\text{new},S}}{N_{\text{old},S}+N_{\text{new},S}},$
(18) $\displaystyle\hat{\sigma}^{2}_{S}$
$\displaystyle=\frac{(N_{\text{old},S}-1)\hat{\sigma}_{\text{old},S}^{2}+(N_{\text{new},S}-1)\hat{\sigma}_{\text{new},S}^{2}+\frac{N_{\text{old},S}N_{\text{new},S}}{N_{\text{old},S}+N_{\text{new},S}}(\hat{\mu}_{\text{old},S}-\hat{\mu}_{\text{new},S})^{2}}{N_{\text{old},S}+N_{\text{new},S}-1}\;,$
(19)

which follows the ideas of Welford’s online algorithm [38]. Here, the
subscripts ’old’ and ’new’ denote two independent sample sets of the same
quantity. These expressions are not only used for the stratum means and
standard deviations, but also to update the corresponding quantities for
potential splits, i.e., possible new strata.

###  4.4 Setting the hybrid sampling parameter $\alpha$

The choice of the hybrid parameter $\alpha\in[0,1]$, which controls the
sampling allocation rule, is affected by two opposing trends: robustness and
variance reduction. In Sect. 3.3 we have discussed that optimal allocation
with $\alpha=1$ can lead to poor performance of the adaptive procedure
described in Algorithm 1 whenever the variances in high-variance strata are
drastically underestimated. This is because the iterative procedure in
Algorithm 1 with $\alpha=1$ may never (e.g., in the case of zero variance
estimates) add new samples to these strata so that the underestimation cannot
be recovered, not even asymptotically. In that sense, the choice $\alpha=1$
lacks robustness. In contrast, proportional allocation (i.e., $\alpha=0$) is
the most robust in that regard as it does not rely on estimates of the strata
standard deviations; cf. also Lemma 1. On the other hand, the results in Sect.
2.3 indicate that optimal allocation with $\alpha=1$ does not only offer
minimal variance for a fixed stratification, but it is also optimal with
respect to the size of a stratification when the function $f$ is
discontinuous. This is in particular relevant for the class of functions of
interest in this work. However, Lemma 3 also reveals that suitable values of
$\alpha$ smaller than $1$ may still offer quasi-optimal variance reduction for
discontinuous functions. The summary of these two opposing effects is
therefore that one should select $\alpha\in[0,1]$ as small as possible yet
large enough to provide quasi-optimal variance reduction. The numerical
examples presented in Sect. 5 contains both results for fixed values of
$\alpha$, and dynamic choice of $\alpha$, to be described next.

####  4.4.1 Dynamic choice of the hybrid parameter

An alternative is to chose $\alpha\in[0,1]$ dynamically across iterations of
the adaptive stratification algorithm described in Algorithm 1. The procedure
that we describe in the following is based on the intuition that at the
beginning of the iterative procedure we may have unreliable estimates of the
strata standard deviations, so that $\alpha$ small is advisable. Eventually,
as the procedure adapts the stratification and increasingly more samples are
being generated, one expects that more reliable estimates of the standard
deviations are available for values of $\alpha<1$, so that $\alpha$ large is
desired for quasi-optimal variance reduction.

To formalize this intuition, let $k\in\mathbb{N}_{0}$ denote the iteration
counter of Algorithm 1 and suppose that the hybrid parameter is initialized
with $\alpha_{0}$ for the first iteration (e.g., $\alpha_{0}=0$ is natural).
Let $\alpha_{k}$ denote the value of the hybrid parameter used during the
$k$-th iteration. At the end of that iteration (line 21 in Algorithm 1), the
value of the hybrid parameter for the iteration $k+1$ should then be selected
in a way that provides a good compromise between variance reduction and
robustness. Here we describe an update based on the stratified sampling
estimators variance $V_{\alpha}=C_{\alpha}(\bm{\sigma})/N$; cf. Sect. 2.3. If
the exact strata standard deviations $\bm{\sigma}$ were known, then the
optimal variance reduction would be achieved for $\alpha=1$. However, the
standard deviations $\bm{\sigma}$ are estimated based on the samples that have
been generated during the iterations so far. Hence, only the empirical
variance constant $C_{\alpha}(\hat{\bm{\sigma}})$ is available. Naively
selecting the hybrid parameter for iteration $k+1$ as
$\alpha_{k+1}=\operatorname*{arg\,min}_{\alpha\in[0,1]}C_{\alpha}(\hat{\bm{\sigma}})$,
where $\hat{\bm{\sigma}}$ contains the estimated strata standard deviations
using the samples available at iteration $k$, may result in spurious effects
due to the lack of robustness, as the empirical variance constant introduces
an error. In fact, in view of the asymptotic distribution characterized in
Lemma 2, the empirical variance constant $C_{\alpha}(\hat{\bm{\sigma}})$ is
random, and it will fluctuate around $C_{\alpha}(\bm{\sigma})$, where those
fluctuations will be asymptotically normally distributed. The size of the
fluctuations can then be approximately quantified via the limit distribution’s
variance, which we write as
$\varsigma_{\alpha}^{2}(\bm{\sigma},\bm{\kappa})=\langle\nabla
C_{\alpha}(\bm{\sigma}),\Sigma_{\alpha}(\bm{\sigma},\bm{\kappa})\nabla
C_{\alpha}(\bm{\sigma})\rangle$ to emphasize the dependence on both the strata
standard deviations $\bm{\sigma}$ and the kurtosis $\bm{\kappa}$ for practical
considerations. Therefore, instead of selecting $\alpha\in[0,1]$ by minimizing
$C_{\alpha}(\hat{\bm{\sigma}})$ at iteration $k$ we incorporate these
fluctuations and consider the minimization of the empirical function
$J_{k}\colon[0,1]\to\mathbb{R}$

$J_{k}(\alpha):=C_{\alpha}(\hat{\bm{\sigma}})+\frac{\varsigma_{\alpha}(\hat{\bm{\sigma}},\hat{\bm{\kappa}})}{\sqrt{N}}\;,$

where the index $k$ indicates that $\hat{\bm{\sigma}}$ and $\hat{\bm{\kappa}}$
contain the estimated strata standard deviations and kurtosis, receptively,
using the $N$ samples available at iteration $k$. That is, for every iteration
$k$, the function $J_{k}$ corresponds to the upper end of the
$68.27\leavevmode\nobreak\ \%$ approximate (asymptotic) confidence interval
for the variance constant $C_{\alpha}(\bm{\sigma})$. Moreover, it satisfies
$\lim_{k\to\infty}J_{k}(\alpha)=C_{\alpha}(\bm{\sigma})$ almost surely, since
$k\to\infty$ implies $N\to\infty$. One therefore expects that the update by
$\alpha_{k+1}=\operatorname*{arg\,min}_{\alpha\in[0,1]}J_{k}(\alpha)$
satisfies $\alpha_{k}\to 1$ asymptotically as $k\to\infty$. To further
robustify the iterative selection of $\alpha$, one could, for example, choose
$\alpha_{k+1}$ as the smallest value in $[0,1]$ that already provides
$\tau\cdot 100\,\%$, $0<\tau\leq 1$ of the optimal (i.e., minimal) upper
confidence band for the variance constant. The complete selection procedure
outlined above is summarized in Algorithm 2.

Algorithm 2 Updating the hybrid sampling parameter $\alpha$.

1:Let $\alpha_{k}$ be the hybrid parameter used for the sample allocation at
the $k$-th iteration of the adaptive stratification procedure.

2:Estimate $\hat{\bm{\sigma}}$ and $\hat{\bm{\kappa}}$ based on the available
$N^{\alpha_{k}}_{S}$ samples in each stratum $S$ of the current stratification
$\mathcal{S}$ at iteration $k$.

3:Determine $J_{k}^{\ast}=\min_{\alpha\in[0,1]}J_{k}(\alpha)$

4:Set $\alpha_{k+1}=\min\bigl{\\{}\alpha\in[0,1]\colon
J_{k}(\alpha)-J_{k}^{\ast}\leq(1-\tau)J_{k}^{\ast}\bigr{\\}}$, or
$\alpha_{k+1}=0$ if the set is empty.

Notice that the choice $\tau=0$ coincides with the case of no additional
robustification. Moreover, step 2 of Algorithm 2 requires estimating both
$\hat{\bm{\sigma}}$ and $\hat{\bm{\kappa}}$ based on the available samples at
iteration $k$. To minimize the risk of drastically misestimating the strata
variances (and thus the kurtosis), cf. inequality (15), we propose to use KDE-
based techniques moment estimates discussed in Sect. 3.3, which introduce
additional smoothing for strata with small observed variability.

##  5 Numerical results

As discussed in Sect. 2, for a given stratification (and if exact sample
allocation rules are available) stratified sampling yields an unbiased
estimator of $\mathbb{E}\left(Q\right)=\mathbb{E}\left(f(\mathbf{Y})\right)$,
whose variance is reduced compared to that of a classic Monte Carlo estimator;
cf. (6). In fact, the ratio of the Monte Carlo estimator’s variance
$\mathbb{V}\left(\hat{Q}_{\textup{MC}}\right)=\mathbb{V}\left(Q\right)/N$ and
a hybrid stratified sampling estimator’s variance
$\mathbb{V}\left(\hat{Q}_{\alpha}\right)$ gives rise to as a measure of the
speedup in view of the central limit theorem (3). Indeed, the speedup

$\textup{speedup}\equiv\frac{\mathbb{V}\left(\hat{Q}_{\textup{MC}}\right)}{\mathbb{V}\left(\hat{Q}_{\alpha}\right)}=\frac{\mathbb{V}\left(Q\right)}{C_{\alpha}(\bm{\sigma})}\;,$
(20)

can be seen as the factor of how many more simulations (i.e., evaluations of
the model $f$) are needed to achieve a certain mean squared error goal when
using a Monte Carlo estimator compared to a hybrid stratified sampling
estimator.

The adaptive stratification procedure introduced in this work aims at
minimizing the stratification estimator’s variance by consecutively refining
high variance strata. That is, the procedure is designed to maximize the
speedup even in the presence of heterogeneous features of $f$. A consequence
of this dynamic sampling-based adaptation is that it is no longer guaranteed
that the estimator at the final iteration will be a statistically unbiased
estimator of $\mathbb{E}\left(Q\right)$ for a finite sample size
$N_{\text{max}}$, due to the dependencies on the previously sampled variables
and strata refinements. However, taking advantage of the hybrid sample
allocation and following the directives introduced in Sect. 3, the adaptive
stratified sampling estimator is still asymptotically unbiased for any
$\alpha\in[0,1)$. Moreover, numerical comparisons (not shown here) of the
empirical adaptive stratified sampling estimators’ bias with the empirical
Monte Carlo bias for the test cases considered in the following indicated that
both terms are of similar size.

Recall that in this work, we assume that the total cost of obtaining
$N_{\text{max}}$ samples of the function of interest is vastly higher than the
cost of the dynamic stratification and allocation of samples. In what follows,
the performance of the algorithm presented in previous sections is evaluated
on a range of different test cases. We note that we do not present trivial
cases where e.g., a discontinuity can be exactly tessellated by the kind of
stratification chosen here (hyperrectangles or simplices), but report that for
such cases the speedups compared to standard Monte Carlo become arbitrarily
large. A discontinuity defined by $1/2^{n}$ of a hypersphere is thus a
nontrivial and more insightful test case where the number of stochastic
dimensions $n$ can be varied to investigate how the performance depends on the
dimensionality of the problem (Sect. 5.1). The variance is strictly localized
to strata that contain the boundary of the hypersphere, and zero elsewhere.
That is, the challenging aspect of these test cases for the adaptive
stratification procedure is the identification (and adaptive isolation) of a
discontinuity. Next, an idealized model for the critical pressure in a faulted
reservoir based on analytical physical expressions is investigated (Sect.
5.2). The model exhibits both discontinuities and regions of sharp but
continuous variation. Finally, we perform numerical tests on the Sod test case
from fluid dynamics, and a vertical equilibrium porous medium model describing
subsurface CO2 storage (Sects. 5.3 and 5.4). All these cases exhibit
discontinuous dependence of the quantity of interest (QoI) $Q=f(\mathbf{Y})$
with respect to the stochastic parameters $\mathbf{Y}$, as is depicted in Fig.
5 and described in more details in subsequent subsections.

(a) Fault stress problem, Sect. 5.2. (b) Sod test case, Sect. 5.3.

(c) CO2 storage problem, Sect. 5.4.

Figure 5: Test problems $Q=f(\mathbf{Y})$. Here, $Q$ is the stress threshold
(MPa), density (kg/m3), and normalized CO2 plume height (dimensionless),
respectively.

For all test problems, we present results using hybrid sampling with
$\alpha=0$ (proportional sampling), $\alpha=0.9$, and $\alpha$ determined
dynamically on $[0,\ 0.95]$ using Algorithm 2.

###  5.1 Hyperspherical discontinuity in multiple dimensions

Consider the unit step function with all coordinates within the unit
intervals, and restricted by the hypersphere in $n$ dimensions with center at
the origin and radius $r_{n}$ chosen so that the hypersphere’s volume is equal
to $2^{n-1}$,

$f(\mathbf{Y})=\left\\{\begin{array}[]{ll}1&\mbox{if
}\|\mathbf{Y}\|_{2}^{2}\leq r_{n}^{2}\mbox{ and }0\leq Y_{i}\leq
1,i=1,\dots,n,\\\ 0&\mbox{otherwise}\end{array}\right..$ (21)

An example of the development of the adaptive stratification is shown for the
$n=2$ case in Fig. 6. Starting from 30 uniformly distributed samples in the
unit square, a total of 10,000 samples are distributed adaptively, resulting
in a final stratification consisting of 17 strata. The samples are
progressively allocated to strata that contain the discontinuity, i.e., have
non-zero variances. In addition, 10% of samples every iteration are allocated
proportionally to enable detecting unseen solution features.

(a) 2 strata. (b) 3 strata. (c) 10 strata. (d) 17 strata.

Figure 6: Evolution of simplex tessellation for the quarter of a circle
(hypersperical $n=2$) test case. Hybrid sampling with $\alpha=0.9$ and 30
samples per simplex added at every iteration with a total of 10,000 samples.

Speedups for $n=2,3,4$ are shown in Figs. 7-9, respectively. For $n=2$ and the
largest sample sizes, the speedups are 2-3 orders of magnitude compared to
standard Monte Carlo sampling. For $n=3$, the corresponding speedups are
10-40, and for $n=4$ typically 5-10. Small values of $c$, intuitively
corresponding to more greedy yet more “risky” stratifications in view of Sect.
4.3, tend to slightly increase speedup, but the pattern is not unambiguous.
Dynamically choosing $\alpha$ in most cases resembles sampling with fixed
$\alpha=0.9$ with some notable differences for smaller sample sizes and $n=4$.

(a) Hyperrect., $\alpha=0$. (b) Hyperrect., $\alpha=0.9$. (c) Hyperrect.,
dynamic $\alpha$. (d) Simplices, $\alpha=0$.

(e) Simplices, $\alpha=0.9$. (f) Simplices, dynamic $\alpha$.

Figure 7: Hyperspherical discontinuity problems with $n=2$ random dimensions.
Hyperrectangular and simplex tessellation, sampling with $\alpha$ being 0, 0.9
or set dynamically, and 1000 repetitions, variable $N_{\textup{max}}$ and
sampling constant $c$.

(a) Hyperrect., $\alpha=0$. (b) Hyperrect., $\alpha=0.9$. (c) Hyperrect.,
dynamic $\alpha$. (d) Simplices, $\alpha=0$.

(e) Simplices, $\alpha=0.9$. (f) Simplices, dynamic $\alpha$.

Figure 8: Hyperspherical discontinuity problems with $n=3$ random dimensions.
Hyperrectangular and simplex tessellation, sampling with $\alpha$ being 0, 0.9
or set dynamically, and 1000 repetitions, variable $N_{\textup{max}}$ and
sampling constant $c$.

(a) Hyperrect., $\alpha=0$. (b) Hyperrect., $\alpha=0.9$. (c) Hyperrect.,
dynamic $\alpha$. (d) Simplices, $\alpha=0$.

(e) Simplices, $\alpha=0.9$. (f) Simplices, variable $\alpha$.

Figure 9: Hyperspherical discontinuity problems with $n=4$ random dimensions.
Hyperrectangular and simplex tessellation, sampling with $\alpha$ being 0, 0.9
or set dynamically, and 1000 repetitions, variable $N_{\textup{max}}$ and
sampling constant $c$.

###  5.2 Geomechanics fault surface stress problem

Consider injection through a well into a porous medium with unknown
permeability and unknown location of a fault. We employ a simplified model for
the surface stress at a fault assuming single-phase flow and no poro-elastic
effects of injection, as in [25]. The background stress field only changes at
the injection surface, and hydrostatic pressure is assumed. The QoI is the
stress threshold $S_{\text{thres}}$ that indicates fault stability,

$S_{\text{thres}}=\left\\{\begin{array}[]{ll}\Delta\tau-\Delta\sigma\mu_{F}\Delta
F&\mbox{ if }\Delta\text{CFF}<0,\\\ 0&\mbox{ otherwise}\end{array}\right.$

where $\Delta\text{CFF}=\Delta\sigma\mu_{F}-\Delta\tau$ is a Coulomb fault
failure criterion. The friction drop $\Delta F$ is 0.8, the shear stress
$\Delta\tau$ is 20 MPa, the normal stress $\Delta\sigma$ is the difference
between an original stress assumed to be 50 MPa and the pressure
$p=\rho_{w}g((d+0.5H)+h)$; where $\rho_{w}=1000$ kg/m3 (density water),
$g=9.81$ (gravitational constant), $d=2000$ m (depth of caprock), $H=100$ m
(height of injection formation), and the hydraulic head $h$ is a function of
the well distance and given by Eq. (2.46) in [25]. The stochastic input
parameters are the well distance (in meters) with distribution $U[10,1000]$,
and the friction coefficient $\mu_{F}$ (dimensionless) with lognormal(0.2,
0.7) distribution.

Figure 10 shows the variance ratios for the fault surface problem as a
function of total number of samples, for $\alpha=0,0.9$ and dynamic $\alpha$
using hyperrectangles and simplices for the adaptive stratifications. For this
problem, the best performance is achieved for hyperrectangular stratifications
with proportional allocation ($\alpha=0$) with speedups exceeding three orders
of magnitude compared to Monte Carlo for the largest sample sizes
investigated. This contrasts with the hyperspherical test cases previously
presented, where simplices and nearly optimal allocation yield the best
results.

(a) Hyperrect., $\alpha=0$.

(b) Hyperrect., $\alpha=0.9$.

(c) Hyperrect., dynamic $\alpha$. (d) Simplices, $\alpha=0$.

(e) Simplices, $\alpha=0.9$. (f) Simplices, dynamic $\alpha$.

Figure 10: Fault stress problem. Hyperrectangular and simplex tessellation,
sampling with $\alpha$ being 0, 0.9 or set dynamically, and 1000 repetitions,
variable $N_{\textup{max}}$ and sampling constant $c$.

###  5.3 Sod shock tube problem

The Sod test case describes the time-dependent behavior of two gas phases with
different phase properties, separated by a membrane that is instantly removed
at time 0. The problem can be modeled by the one-dimensional Euler equations
and has been investigated in the uncertainty quantification literature to
illustrate challenges related to discontinuities in stochastic space [29, 36,
28]. We assume the following three uncertain parameters: left state density
$U[0.7,\ 1.3]$, right state density $U[0.05,\ 0.2]$, and initial location of
the membrane with distribution $U[0.45,\ 0.55]$. The quantity of interest is
the density at $(x,t)=(0.7,0.1)$.

The speedups for hyperrectangles and simplices, hybrid sampling with
$\alpha=0$, $\alpha=0.9$, and dynamic $\alpha$, are shown in Fig. 11. The best
performance is obtained with hyperrectangular stratification, despite a curved
hypersurface discontinuity, see Fig. 5(b).

(a) Hyperrect., $\alpha=0$. (b) Hyperrect., $\alpha=0.9$. (c) Hyperrect.,
dynamic $\alpha$. (d) Simplices, $\alpha=0$.

(e) Simplices, $\alpha=0.9$. (f) Simplices, dynamic $\alpha$.

Figure 11: Sod test case described by the Euler equations. Hyperrectangular
and simplex tessellation, sampling with $\alpha$ being 0, 0.9 or set
dynamically, and 1000 repetitions, variable $N_{\textup{max}}$ and sampling
constant $c$.

###  5.4 Vertical equilibrium model of CO2 storage

Subsurface permanent CO2 storage involves injection of large quantities of CO2
into porous reservoirs with largely unknown physical properties. The varying
temporal and physical scales calls for simplified-physics models, achieved,
e.g., by vertical equilibrium assumptions and dimension reduction by
integration over the vertical direction [25]. The result is a nonlinear
hyperbolic conservation law with discontinuous flux function, where the
solution represents the saturation of CO2 as a function of space and time. For
reproducibility, Table 1 contains all values relevant for the setup of the
problem, including distributions for four random parameters: permeability,
brine mobility, CO2 mobility, and background flow rate. A more detailed
description of the problem, including the conservation law and its time-
dependent solution, can be found in [22, 27]. Although not explicitly referred
to in this paper, Table 1 lists notation for various parameters consistent
with [27], also to facilitate reproducibility.

Parameter | Notation | Value/Distribution
---|---|---
Porosity | $\phi$ | 0.15
Residual brine saturation | $S_{\textup{br}}$ | 0.1
Residual CO2 saturation | $S_{\textup{cr}}$ | 0.1
Slope angle | $\theta$ | 0.005
Injection time | $\tau$ | 20 years
Injection rate | $Q_{\textup{inj}}$ |  $1\times 10^{-7}$ m/s
CO2 mobility | $\lambda_{\textup{c}}$ | Uniform $[0.7,\ 1.3]\times 6.25\times 10^{-5}$ ms/kg
Brine mobility | $\lambda_{\textup{b}}$ | Uniform $[0.8,\ 1.2]\times 5\times 10^{-4}$ ms/kg
Background flow | $Q$ | Exponential, mean $1\times 10^{-9}$ m/s
Permeability | $k$ | Lognormal, mean 200 mD, std 50 mD
Table 1:  Parameter setup for the CO2 storage test problem.

The QoI is the height of the CO2 plume 100 m downstream of the injection
location, 600 years after the end of injection. Speedups for the CO2 problem
as a function of the total number of samples are shown in Fig. 12. Speedup is
observed for all sampling schemes, even for the smallest sizes of sample sets.
This problem exhibits zero or small variability away from the discontinuities,
as indicated in Fig. 5(c), making it a particularly suitable candidate for the
proposed method. At the same time, the number of dimensions ($n=3$ due to $Q$
being a function of the ratio of the phase mobilities, rather than their
individual values) already requires a large number of samples to give speedups
greater than an order of magnitude.

(a) Hyperrect., $\alpha=0$. (b) Hyperrect., $\alpha=0.9$. (c) Hyperrect.,
dynamic $\alpha$. (d) Simplices, $\alpha=0$.

(e) Simplices, $\alpha=0.9$. (f) Simplices, dynamic $\alpha$.

Figure 12: CO2 storage problem. Hyperrectangular and simplex tessellation,
sampling with $\alpha$ being 0, 0.9 or set dynamically, and 1000 repetitions,
variable $N_{\textup{max}}$ and sampling constant $c$.

##  6 Conclusions

We have introduced a novel stratified sampling method with adaptive
stratification and sequential, hybridized allocation of samples. For a fixed
stratification, the sample allocation asymptotically approaches a prescribed
linear combination of proportional and optimal allocation. Letting a fraction
of samples be allocated proportionally adds robustness to the adaptive method,
as the stratum standard deviations are not a priori known but estimated on-
the-fly as samples are added. Moreover, a greedy approach is used to split the
stratum that results in the largest reduction of the stratified sampling
estimator’s current variance. Kurtosis dependent estimates are provided that
allow to quantify the probability of drastically underestimating local strata
variances and, thus, of failing to identify strata that should be refined.

To maintain prescribed stratum sampling rates and for ease of implementation,
a new stratum that arises after re-stratification should be contained within a
single parent stratum, and it should retain the same geometrical shape. In
this work, we suggest using hyperrectangles or simplices for the
stratification, as either of these shapes can be bisected and the result is
two new hyperrectangles or simplices. In fact, both classes of tessellations
allow for a flexible partition of the stochastic domain.

The proposed method is anticipated to result in significant speedups for
problems where the variability is localized in random space, e.g., the PDE
solutions describing physical problems with uncertain parameters and
exhibiting steep gradients or discontinuities. In contrast to, e.g., localized
response surface methods based on function approximations, an advantage of the
proposed method is that accurate identification of steep features is not
necessary to obtain good results. Confining sharp features to strata of small
measure is often sufficient to vastly outperform standard Monte Carlo
sampling. This has been verified experimentally through various test cases,
exhibiting speedups of up to three orders of magnitude compared to standard
Monte Carlo.

## Acknowledgements

The authors would like to thank Eirik Keilegavlen at the University of Bergen
for the setup of the fault surface stress problem. The first author was funded
by the Research Council of Norway through the project Quantification of fault-
related leakage risk (FRISK) under project number 294719. A CC BY or
equivalent licence is applied to the author accepted manuscript arising from
this submission, in accordance with the grant’s open access conditions.

##  Appendix A Derivation of the gradient of the variance constant

Introduce the functions

$g_{S}\colon\mathbb{R}^{\lvert\mathcal{S}\rvert}\to\mathbb{R}\;,\quad\bm{\sigma}\mapsto
g_{S}(\bm{\sigma})=\frac{\sigma_{S}}{\langle\bm{p},\bm{\sigma}\rangle}$

for any $S\in\mathcal{S}$. Then

$\bm{\sigma}\mapsto
C_{\alpha}(\bm{\sigma})=\sum_{S\in\mathcal{S}}\frac{p_{S}\sigma_{S}^{2}}{1+\alpha\left(g_{S}(\bm{\sigma})-1\right)}\;,$

so that $V_{\alpha}=C_{\alpha}(\bm{\sigma})/N$. To quantify the robustness of
the variance reduction with respect to perturbation in $\bm{\sigma}$, we
compute the gradient of $C_{\alpha}$. Therefore, we first compute the partial
derivative of $g_{S}$ with respect to the component $\sigma_{U}$ of
$\bm{\sigma}$ as

$\partial_{\sigma_{U}}g_{S}(\bm{\sigma})=\frac{\delta_{S,U}}{\langle\bm{p},\bm{\sigma}\rangle}-\frac{p_{U}\sigma_{S}}{{\lvert\langle\bm{p},\bm{\sigma}\rangle\rvert}^{2}}\;,$

where $\delta_{S,U}$ denotes the Kronecker delta: $\delta_{S,U}=1$ if $S=U$
and $\delta_{S,U}=0$ else. Next we compute the component of the gradient of
$C_{\alpha}$ with respect to $\sigma_{U}$, which, after some algebra, can be
written as:

$\partial_{\sigma_{U}}C_{\alpha}(\bm{\sigma})=\frac{p_{U}\sigma_{U}}{1+\alpha\left(g_{U}(\bm{\sigma})-1\right)}\left(1+\frac{1-\alpha}{1+\alpha\left(g_{U}(\bm{\sigma})-1\right)}\right)+\frac{\alpha
p_{U}}{{\lvert\langle\bm{p},\bm{\sigma}\rangle\rvert}^{2}}\sum_{S\in\mathcal{S}}\frac{p_{S}\sigma_{S}^{3}}{{\left(1+\alpha\left(g_{S}(\bm{\sigma})-1\right)\right)}^{2}}\;.$

Finally, using the fact that

$1+\alpha\left(g_{S}(\bm{\sigma})-1\right)=\frac{\alpha\sigma_{S}+(1-\alpha)\langle\bm{p},\bm{\sigma}\rangle}{\langle\bm{p},\bm{\sigma}\rangle}$

we can eventually write the partial derivative of $C_{\alpha}$ as

$\partial_{\sigma_{U}}C_{\alpha}(\bm{\sigma})=\frac{p_{U}\sigma_{U}\langle\bm{p},\bm{\sigma}\rangle}{\alpha\sigma_{U}+(1-\alpha)\langle\bm{p},\bm{\sigma}\rangle}\left(1+\frac{(1-\alpha)\langle\bm{p},\bm{\sigma}\rangle}{\alpha\sigma_{U}+(1-\alpha)\langle\bm{p},\bm{\sigma}\rangle}\right)+\alpha
p_{U}\sum_{S\in\mathcal{S}}\frac{p_{S}\sigma_{S}^{3}}{{\left(\alpha\sigma_{S}+(1-\alpha)\langle\bm{p},\bm{\sigma}\rangle\right)}^{2}}\;.$

##  Appendix B Further motivation for splitting adaptively

###  B.1 Variance reduction by splitting

Consider an arbitrary stratum $S$ with a partition $S=S_{+}\cup S_{-}$
defining a refined stratification. We are interested in the effect of
splitting on the variance of the estimator. Without loss of generality, assume
that $\mathbb{E}\left(f(\mathbf{Y})|\mathbf{Y}\in S\right)=0$ (we may always
subtract a constant without changing the variance), and let the relative
measure of $S_{+}$ and $S_{-}$ be $q$ and $1-q$ (the measure of $S$ itself
does not matter). Then,

$\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in
S\right)=\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{+}\right)+\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{-}\right)+2\mathbb{C}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{+},f(\mathbf{Y})|\mathbf{Y}\in S_{-}\right)\\\
=\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{+}\right)+\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{-}\right)+\frac{2(1-q)}{q}\left(\mathbb{E}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{-}\right)\right)^{2}\\\ \geq q\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{+}\right)+(1-q)\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in S_{-}\right),$
(22)

where the second equality follows from $S_{+}$ and $S_{-}$ being disjoint and
the relation

$\mathbb{E}\left(f(\mathbf{Y})|\mathbf{Y}\in
S\right)=q\mathbb{E}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{+}\right)+(1-q)\mathbb{E}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{-}\right)=0\;.$

The right-hand side of Eq. (22) equals the variance contribution to the
estimator from the two newly created strata. Thus, we have shown that
splitting of a stratum always reduces the variance of the estimator. that
optimal allocation is maintained after splitting.

###  B.2 Stratification of domains of discontinuous functions

Let $S_{+}$ and $S_{-}$ be two subdomains that forms a partition of
$\mathfrak{U}=[0,\ 1]^{n}$, separated by a hypersurface
$H\subset\mathbb{R}^{n-1}$. Let $f$ be a piecewise constant function on
$\mathbb{R}^{n}$, defined by

$f(\mathbf{Y})=\left\\{\begin{array}[]{ll}c_{+}&\mbox{if }\mathbf{Y}\in
S_{+}\\\ c_{-}&\mbox{if }\mathbf{Y}\in S_{-}\end{array}\right.,$ (23)

where $c_{+}$ and $c_{-}$ are constants, and consider a neighborhood $S_{H}$
around $H$. Denote the measure of $S_{+}$ and $S_{-}$ by
$p_{+}=\mathbb{P}\left(S_{+}\right)$ and $p_{-}=\mathbb{P}\left(S_{-}\right)$,
respectively, and $\tilde{p}_{+}=\mathbb{P}\left(S_{+}\cap S_{H}\right)$ and
$\tilde{p}_{-}=\mathbb{P}\left(S_{-}\cap S_{H}\right)$. Setting
$p=\tilde{p}_{+}/p_{S_{H}}$ we have

$\mathbb{V}\left(f(\mathbf{Y})|\mathbf{Y}\in
S_{H}\right)=\frac{\tilde{p}_{+}}{p_{S_{H}}}\left(1-\frac{\tilde{p}_{+}}{p_{S_{H}}}\right)c_{+}^{2}+\frac{\tilde{p}_{-}}{p_{S_{H}}}\left(1-\frac{\tilde{p}_{-}}{p_{S_{H}}}\right)c_{-}^{2}-2\frac{\tilde{p}_{+}\tilde{p}_{-}}{p_{S_{H}}^{2}}c_{+}c_{-}\\\
=p(1-p)(c_{+}-c_{-})^{2},$ (24)

where we have used that $p_{+}+p_{-}=1$, and
$\tilde{p}_{+}+\tilde{p}_{-}=p_{S_{H}}$. Note that (24) is similar to the
expression for the global variance,

$\mathbb{V}\left(f(\mathbf{Y})\right)=p_{+}(1-p_{+})(c_{+}-c_{-})^{2}.$

To minimize the variance, we see from (24) that the stratum $S_{H}$ should be
small, whereas the local variances in $S_{+}\setminus S_{H}$ and
$S_{-}\setminus S_{H}$ are identically zero (since $f$ is constant).

## References

  * [1] S. Asmussen and P. W. Glynn.  Stochastic simulation: algorithms and analysis, volume 57.  Springer, New York, 2007.
  * [2] L. Cabral Pereira, M. Sánchez, and L. J. do Nascimento Guimarães.  Uncertainty quantification for reservoir geomechanics.  Geomech. Energy Envir., 8:76–84, 2016.  Themed Issue on Selected Papers Symposium of Energy Geotechnics 2015 — Part II.
  * [3] M. Christie, V. Demyanov, and D. Erbas.  Uncertainty quantification for porous media flows.  J. Comput. Phys., 217(1):143–158, 2006.
  * [4] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup.  Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients.  Comput. Visual Sci., 14(1):3–15, 2011.
  * [5] W. G. Cochran.  Sampling Techniques.  John Wiley & Sons, New York, 3rd edition, 1977.
  * [6] François Cuvelier and Gilles Scarella.  Vectorized algorithms for regular tessellations of d-orthotopes and their faces.  working paper or preprint, November 2017.
  * [7] A. DasGupta.  Asymptotic theory of statistics and probability.  Springer, New York, 2008.
  * [8] C. De Luigi and S. Maire.  Adaptive integration and approximation over hyper-rectangular regions with applications to basket options pricing.  Monte Carlo Methods Appl., 16:265–282, 12 2010.
  * [9] J. Dick, F. Y. Kuo, and I. H. Sloan.  High-dimensional integration: The quasi-Monte Carlo way.  Acta Numerica, 22:133–288, 2013.
  * [10] P. Étoré, G. Fort, B. Jourdain, and E. Moulines.  On adaptive stratification.  Ann. Oper. Res., 89(1):127–154, 2011.
  * [11] P. Étoré and B. Jourdain.  Adaptive optimal allocation in stratified sampling methods.  Methodol. Comput. Appl., 12(3):335–360, 2010.
  * [12] M. B. Giles.  Multilevel Monte Carlo methods.  Acta Numer., 24:259–328, 2015.
  * [13] M. B. Giles and A.-L. Haji-Ali.  Multilevel nested simulation for efficient risk estimation.  SIAM/ASA J. Uncertain. Quantif., 7(2):497–525, 2019.
  * [14] A. Gorodetsky and Y. Marzouk.  Efficient localization of discontinuities in complex computational simulations.  SIAM J. Sci. Comput., 36(6):A2584–A2610, 2014.
  * [15] A.-L. Haji-Ali, F. Nobile, and R. Tempone.  Multi Index Monte Carlo: when sparsity meets sampling.  Numer. Math., 132(4):767–806, 2016.
  * [16] S. Hosder, R. Walters, and M. Balch.  Efficient uncertainty quantification applied to the aeroelastic analysis of a transonic wing.  In 46th AIAA Aerospace Sciences Meeting and Exhibit, 2008.
  * [17] J. D. Jakeman, A. Narayan, and D. Xiu.  Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions.  J. Comput. Phys., 242:790– 808, 2013.
  * [18] D. P. Kroese, T. Taimre, and Z. I. Botev.  Handbook of Monte Carlo methods.  Wiley New Jersey, 2011.
  * [19] S. Krumscheid and F. Nobile.  Multilevel Monte Carlo Approximation of Functions.  SIAM/ASA J. Uncertain. Quantif., 6(3):1256–1293, 2018.
  * [20] S. Krumscheid, F. Nobile, and M. Pisaroni.  Quantifying uncertain system outputs via the multilevel Monte Carlo method—Part I: Central moment estimation.  J. Comput. Phys., 414, 2020.
  * [21] H. W. Kuhn.  Some combinatorial lemmas in topology.  IBM J. Res. Dev., 4(5):518–524, November 1960.
  * [22] C. W. MacMinn, M. L. Szulczewski, and R. Juanes.  CO2 migration in saline aquifers. part 1. capillary trapping under slope and groundwater flow.  J. Fluid Mech., 662:329–351, 2010.
  * [23] F. Müller, P. Jenny, and D. W. Meyer.  Multilevel Monte Carlo for two phase flow and Buckley-Leverett transport in random heterogeneous porous media.  J. Comput. Phys., 250:685–702, 2013.
  * [24] H. Niederreiter.  Random number generation and quasi-Monte Carlo methods.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1992.
  * [25] J.M. Nordbotten and M.A. Celia.  Geological Storage of CO2: Modeling Approaches for Large-Scale Simulation.  Wiley, 2011.
  * [26] T. N. Palmer.  Predicting uncertainty in forecasts of weather and climate.  Rep. Prog. Phys., 63(2):71–116, 2000.
  * [27] P. Pettersson.  Stochastic Galerkin formulations for CO2 transport in aquifers: Numerical solutions with uncertain material properties.  Transport Porous Med., 114(2):457–483, 2016.
  * [28] P. Pettersson, G. Iaccarino, and J. Nordström.  A stochastic Galerkin method for the Euler equations with Roe variable transformation.  J. Comput. Phys., 257:481–500, 2014.
  * [29] G. Poëtte, B. Després, and D. Lucor.  Uncertainty quantification for systems of conservation laws.  J. Comput. Phys., 228(7):2443–2467, 2009.
  * [30] W. H. Press and G. R. Farrar.  Recursive stratified sampling for multidimensional Monte Carlo integration.  Comput. Phys., 4(2):190–195, 1990.
  * [31] A. Rushdi, L. P. Swiler, E. T. Phipps, M. D’Elia, and M. S. Ebeida.  Vps: Voronoi piecewise surrogate models for high-dimensional data fitting.  Int. J. Uncertain. Quan., 7(1):1–21, 2017.
  * [32] R. J. Serfling.  Approximation theorems of mathematical statistics.  John Wiley & Sons, Inc., New York, 1980.
  * [33] M. D. Shields.  Refined latinized stratified sampling: A robust sequential sample size extension methodology for high-dimensional latin hypercube and stratified designs.  Int. J. Uncertain. Quan., 6(1):79–97, 2016.
  * [34] M. D. Shields, K. Teferra, A. Hapij, and R. P. Daddazio.  Refined stratified sampling for efficient Monte Carlo based uncertainty quantification.  Reliab. Eng. Syst. Safe., 142:310–325, 2015.
  * [35] D. M. Tartakovsky and S. Broyda.  PDF equations for advective-reactive transport in heterogeneous porous media with uncertain properties.  J. Contam. Hydrol., 120-121:129–140, 2011.
  * [36] J. Tryoen, O. Le Maître, M. Ndjinga, and A. Ern.  Intrusive Galerkin methods with upwinding for uncertain nonlinear hyperbolic systems.  J. Comput. Phys., 229(18):6485–6511, 2010.
  * [37] R. Vershynin.  High-dimensional probability, volume 47 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press, Cambridge, 2018.  An introduction with applications in data science.
  * [38] B. P. Welford.  Note on a method for calculating corrected sums of squares and products.  Technometrics, 4:419–420, 1962.
  * [39] J. A. S. Witteveen and G. Iaccarino.  Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces.  SIAM J. Sci. Comput., 34(2):A814–A838, 2012.
  * [40] D. Xiu.  Numerical methods for stochastic computations: A spectral method approach.  Princeton University Press, Princeton, NJ, 2010.