# Nonasymptotic Convergence Rate of Quasi-Monte Carlo: Applications to Linear
Elliptic PDEs with Lognormal Coefficients and Importance Samplings

Yang Liu111KAUST SRI Center for Uncertainty Quantification in Computational
Science and Engineering 0000-0003-0778-3872<EMAIL_ADDRESS>Raúl
Tempone22footnotemark: 2333Alexander von Humboldt Professor in Mathematics for
Uncertainty Quantification, RWTH Aachen University, 52062 Aachen, Germany
0000-0003-1967-4446 Computer, Electrical and Mathematical Sciences and
Engineering,
4700 King Abdullah University of Science and Technology (KAUST),
Thuwal 23955-6900, Kingdom of Saudi Arabia

###### Abstract

This study analyzes the nonasymptotic convergence behavior of the quasi-Monte
Carlo (QMC) method with applications to linear elliptic partial differential
equations (PDEs) with lognormal coefficients. Building upon the error analysis
presented in (Owen, 2006), we derive a nonasymptotic convergence estimate
depending on the specific integrands, the input dimensionality, and the finite
number of samples used in the QMC quadrature. We discuss the effects of the
variance and dimensionality of the input random variable. Then, we apply the
QMC method with importance sampling (IS) to approximate deterministic, real-
valued, bounded linear functionals that depend on the solution of a linear
elliptic PDE with a lognormal diffusivity coefficient in bounded domains of
$\mathbb{R}^{d}$, where the random coefficient is modeled as a stationary
Gaussian random field parameterized by the trigonometric and wavelet-type
basis. We propose two types of IS distributions, analyze their effects on the
QMC convergence rate, and observe the improvements.

###### keywords:

Quasi-Monte Carlo , Importance sampling , Finite elements , Partial
differential equations with random data , Lognormal diffusion , Wavelets

###### MSC:

[2020] 65C05 , 65N50 , 65N22 , 35R60

##  1 Introduction

The quasi-Monte Carlo (QMC) method computes the expectation of a random
variable using deterministic low-discrepancy sequences. The QMC method offers
a convergence rate of $\mathcal{O}(n^{-1+\epsilon})$ for $\epsilon>0$,
surpassing the Monte Carlo convergence rate of $\mathcal{O}(n^{-1/2})$.
However, the effectiveness of the QMC method depends on the regularity of the
integral and the integration dimension. The integrand variation and
integration dimension dictate the classical upper bound for the QMC method,
known as the Koksma–Hlawka inequality [34].

Concerns about the quality of point sets arise as the integration dimension
increases [32, 49, 28, 12]. For instance, the Halton sequence may exhibit
pathological behavior in certain dimensions; thus, researchers have proposed
remedies [11, 39].

Despite these challenges, the QMC method has been successful in relatively
high dimensions due to the “low effective dimension” [47, 48]. This notion
first arose from the analysis of variance (ANOVA) decomposition in [9].
Several studies have aimed to apply various approaches to minimize the
effective dimension, particularly in financial applications [7, 6, 46].
Nevertheless, an equivalence theorem in [50], an extension of the worst-case
integration error analysis in [41], stated that no decomposition method
(includes Brownian bridge, principal component analysis, etc.) is consistently
superior to other methods in terms of different payoff functions. Nonetheless,
considering the exact forms of the payoff function, the authors of [23]
introduced the linear transformation method, exploiting the nonuniqueness of
the covariance matrix decomposition to find the most important dimensions. The
design of the transformation matrix’s columns optimizes the variance
contribution of each random variable.

In addition to the effective dimension, the regularity of the integrand also
affects the QMC integration. In finance applications, option pricing problems
can sometimes involve discontinuous functions. The orthogonal transformation
was proposed in [51] to transform the discontinuity, involving a linear
combination of the random variables, to become parallel to the axes. The
authors of [24] discussed the connection between orthogonal and linear
transformation methods. The authors in [19] characterized the “QMC-friendly”
axis-parallel discontinuity and explicitly derived the convergence rate for
integrands with such discontinuities, supported by numerical experiments.
Although the optimal QMC convergence rates are not recovered, the
superiorities to the MC methods are shown.

Conditioning is a well-known approach to reduce variance, but it can also
improve the smoothness [18, 52, 7, 4]. In [17], the authors demonstrated that,
under certain conditions, all terms from the ANOVA decomposition except the
one with the highest order have infinite smoothness. In addition, in [13], the
conditions under which preintegration works were provided. In [30], the
authors proposed preintegrating in a subspace comprising a linear combination
of random input variables. Moreover, previous studies [25, 5] exploited the
regularity with Fourier transformation when the original function exhibits
non-smoothness.

A partial differential equation (PDE) with random coefficients is another
application in which the QMC method accelerates the convergence of computing
the statistical properties of a quantity of interest (QoI), often taking the
form of a functional of the PDE solution. In [44], the authors proved that the
QMC worst-case error is independent of the dimension if the integrand belongs
to a class of weighted Sobolev spaces. Another study [33] considered the
weights of the form “product and order dependent” to minimize the worst-case
error and provide a fast method to design the lattices according to the
weights. The dimension-independent worst-case QMC errors of elliptic PDEs with
affine uniform and lognormal coefficients were analyzed in [29] and [16]. In
the latter case, the space is equipped with additional weight functions.
Moreover, in [21, 26], the authors considered the “product weights” in the
weighted space and lattice rule design. More recent studies [15] and [40] have
considered the median-of-the-mean QMC estimator using lattice rule (without
specific design) and digital sequence, respectively.

Despite the careful analysis of the asymptotic convergence rates in the
literature, the ideal QMC convergence rate $\mathcal{O}(n^{-1})$ is not always
observed in practice. In this work, we aim to study the nonasymptotic
convergence behavior of the QMC method and explain the often observed
suboptimal convergence rates.

Specifically, certain integration problems involve integrand domains different
from $[0,1]$, such as the integration with respect to (w.r.t) the Gaussian
measure on $\mathbb{R}$. The inverse cumulative distribution function (CDF)
transformation is applied for compatibility with QMC methods. However, this
transformation often introduces singularities at the boundaries. The work [37]
characterized the boundary singularity with the boundary growth rate and made
connections with the asymptotic QMC convergence rate. Particularly, the author
provided examples of integrands which blow up at boundaries but still lead to
optimal QMC convergence rate. Building on this work, we aim to analyze the
nonasymptotic QMC convergence rate for some examples involving lognormal
random variables to explain the observed suboptimal rates. Last, the
importance sampling (IS) is a well-known method for variance reduction. We
also aim to discover the extent to which IS with certain proposal
distributions improves the convergence rate.

Some recent publications have delved into the realm of QMC methods for
potentially unbounded integrands. One study [14] combined the robust mean
estimator with QMC sampling. In parallel, the work [20] and [35] studied the
QMC method with IS. While these studies provided valuable insights into the
asymptotic convergence rate, they did not address the nonasymptotic behaviors
frequently observed in real-world scenarios. In contrast, our work contributes
to developing a convergence rate model for QMC methods with a finite number of
samples. This innovative approach offers potential pathways to enhance the
practically-observed convergence rate.

The paper is organized as follows. Section 2 discusses the nonasymptotic
convergence rate. Next, Section 3 presents the convergence rate analysis for
two examples: the expectation of a lognormal random variable and elliptic PDEs
with lognormal coefficients. Then, Section 4 evaluates the effects of two
kinds of IS distributions. Section 5 details the numerical results. Finally,
Section 6 presents the conclusions.

##  2 Nonasymptotic convergence rate for the randomized QMC method

This section follows the proofs in [37] and accordingly modifies them to
establish the nonasymptotic results. First, we introduce the notation. We are
interested in the following integration problem:

$\displaystyle I(g)=\int_{[0,1]^{s}}g(\mathbf{t})d\mathbf{t},$ (1)

where $g:[0,1]^{s}\to\mathbb{R}$. We let
$\mathcal{P}=\\{\mathbf{t}_{1},\mathbf{t}_{2},\dotsc,\mathbf{t}_{n}\\}$ be a
point set in $[0,1]^{s}$.

The QMC estimator for the integrand $g$ is given by

$\displaystyle\hat{I}_{n}(g)=\frac{1}{n}\sum_{i=1}^{n}g(\mathbf{t}_{i}),$ (2)

where the $\\{\mathbf{t_{i}}\\}$ is a predesigned deterministic low-
discrepancy sequence [34, 10]. A notion to describe how well the points in
$\mathcal{P}$ are uniformly distributed is the star discrepancy
$\Delta_{\mathcal{P}}^{*}$, given by

$\displaystyle\Delta_{\mathcal{P}}^{*}:=\sup_{\mathbf{x}\in[0,1]^{s}}\left\lvert\Delta_{n}(\mathbf{x};\mathcal{P})\right\rvert,$
(3)

where
$\Delta_{n}(\mathbf{x};\mathcal{P})=\frac{1}{n}\sum_{j=1}^{n}\mathbbm{1}_{t_{j}\in[0,\mathbf{x})}-\prod_{k=1}^{s}x_{k}$,
which is the difference of the measure of $[0,\mathbf{x})$ and the proportion
of points that belong to this set, where $[0,\mathbf{x})$ denotes the tensor
product of each dimension, i.e. $[0,\mathbf{x})=\prod_{k=1}^{s}[0,x_{k})$. We
expect a small difference for an evenly distributed point set. For some low-
discrepancy sequences with fixed length $n$, we have

$\displaystyle\Delta_{\mathcal{P}}^{*}=\mathcal{O}(n^{-1}(\log n)^{s}),$ (4)

where the exponent of the log term becomes $s-1$ for an infinite sequence [34,
10].

In this work we consider functions with a continuous first-order mixed
derivative. The variation in the Hardy–Krause sense can be computed as
follows:

$\displaystyle
V_{HK}(g)=\sum_{\emptyset\neq\mathfrak{u}\subseteq\\{1,2,\dotsc,s\\}}\int_{[0,1]^{\lvert\mathfrak{u}\rvert}}\left\lvert{\partial^{\mathfrak{u}}g(\mathbf{t^{\mathfrak{u}}}\colon\mathbf{1^{-\mathfrak{u}}})}\right\rvert
d\mathbf{t^{\mathfrak{u}}},$ (5)

where $\lvert\mathfrak{u}\rvert$ is the cardinality of the set $\mathfrak{u}$,
$\mathbf{y}={\mathbf{t^{\mathfrak{u}}}\colon\mathbf{1^{-\mathfrak{u}}}}\in[0,1]^{s}$
denotes a point in $[0,1]^{s}$ with ${y^{j}}={t^{j}}$ for $j\in\mathfrak{u}$,
and $y^{j}=1$ otherwise. For the set
$\mathfrak{u}=\\{\mathfrak{u}_{1},\dotsc,\mathfrak{u}_{\lvert\mathfrak{u}\rvert}\\}\subseteq\\{1,\dotsc,s\\}$,
the mixed derivative $\partial^{\mathfrak{u}}g$ is explicitly given by

$\partial^{\mathfrak{u}}g(\mathbf{t})=\frac{\partial^{\lvert\mathfrak{u}\rvert}g}{\partial\mathbf{{t}}_{\mathfrak{u}}}(\mathbf{{t}})=\frac{\partial}{\partial
t_{\mathfrak{u}_{\lvert\mathfrak{u}\rvert}}}\cdots\frac{\partial}{\partial
t_{\mathfrak{u}_{1}}}g(\mathbf{{t}})$

where the continuity ensures that the order of differentiation can be switched
while maintaining derivative invariance.

The Koksma–Hlawka inequality provides an error estimate for the QMC method,
given by

$\displaystyle\left\lvert I(g)-\hat{I}_{n}(g)\right\rvert\leq
V_{HK}(g)\Delta_{\mathcal{P}}^{*},$ (6)

where $V_{HK}(g)$ is the variation of $g$ in the Hardy–Krause sense.

However, the use of a deterministic point set yields biased results.
Randomization techniques have been introduced to address this problem, giving
rise to the randomized QMC (RQMC) unbiased estimator [31]:

$\displaystyle\hat{I}_{R,n}(g)=\frac{1}{R}\sum_{r=1}^{R}\frac{1}{n}\sum_{i=1}^{n}g(\mathbf{t}_{i}\oplus\bm{\Delta}_{r})=\frac{1}{R}\sum_{r=1}^{R}\hat{I}_{n}^{(r)}(g),$
(7)

where $\mathbf{t}_{i}$ is the $i$th deterministic QMC quadrature point,
$\bm{\Delta}_{r}$ represents the $r$th randomization, and $\oplus$ denotes the
randomization operation. An example of such an operation is the random shift,
where $\bm{\Delta}_{r}\sim U[0,1]^{s}$ and
$\mathbf{a}\oplus\mathbf{b}=(\mathbf{a}+\mathbf{b})\textrm{ mod }1$, with the
modulo taken component-wise. This random shift approach is easy to implement
and provides an unbiased integral estimator. Apart from the unbiasedness, the
randomization also provides an convenient access to the estimator variance.

To prepare for the analysis of the RQMC method, we study the behavior of the
samples from the uniform distribution in $[0,1]^{s}$. More specifically, we
consider the problem of evaluating the integral

$\displaystyle\begin{split}I(g)&=\int_{[0,1]^{s}}g(\mathbf{t})d\mathbf{t}\\\
&=\int_{\mathbb{R}^{s}}g\circ\Phi(\mathbf{y})\rho(\mathbf{y})d\mathbf{y}\end{split}$
(8)

where $\circ$ denotes the function composition,
$\mathbf{y}=\Phi^{-1}(\mathbf{t})$, $\Phi^{-1}:[0,1]^{s}\to\mathbb{R}^{s}$ the
inverse CDF of $s$-dimensional standard normal distribution, and $\rho$
represents its probability density function. Section 3 presents two concrete
examples of $g$. Before we analyze such an integrand, we study the behavior of
uniform random variables.

###  2.1 Some properties of the uniform distribution

Lemma 1 introduces a useful property of the uniform distribution.

######  Lemma 1 (A bound for an $s$-dimensional uniform random variable).

Let each $\mathbf{{t}}_{i}$ in the sequence
$\\{\mathbf{{t}}_{i}\\}_{i=1}^{\infty}$ be uniformly distributed over
$[0,1]^{s}$. Define $E_{n}$ as the event $\\{\prod_{j=1}^{s}t_{n}^{j}\leq
C\cdot n^{-r}\\}$ for $r>1$ and $C>0$. Then, we have

$\displaystyle\textnormal{Pr}\left(E_{n}\enspace i.o.\right)=0,$ (9)

where $i.o.$ stands for “infinitely often.”

The sequence $\\{{t}_{i}\\}$ in Lemma 1 is not necessarily independent, as is
the case for the RQMC method, where the points are desired to exhibit a
negative correlation [53]. Lemma 1 is a slight modification of Lemma 4.1 in
[37], where the minimum condition is removed. The proof follows [37], except
we corrected an upper bound.

###### Proof.

Let $\mathbf{t}_{n}\sim U[0,1]^{s}$. We have

$\begin{split}\textnormal{Pr}\left(\prod_{j=1}^{s}t_{n}^{j}\leq C\cdot
n^{-r}\right)&=\textnormal{Pr}\left(-2\log\left(\prod_{j=1}^{s}t_{n}^{j}\right)\geq
2r\log(n)-2\log(C)\right)\\\ &=\textnormal{Pr}\left(\chi^{2}_{(2s)}\geq
2r\log(n)-2\log(C)\right)\\\
&=\int_{2r\log(n)-2\log(C)}^{\infty}\frac{z^{s-1}e^{-z/2}}{2^{s}\Gamma(s)}dz\\\
&=\int_{r\log(n)-\log(C)}^{\infty}\frac{y^{s-1}e^{-y}}{\Gamma(s)}dy,\end{split}$
(10)

where $\chi^{2}_{(2s)}$ denotes a chi-squared distribution with $2s$ degrees
of freedom and $\Gamma(s)$ denotes the gamma function evaluated at $s$. Now,
we use the upper bound derived in [43] for the incomplete gamma function:

$\displaystyle\Gamma(s,x)=\int_{x}^{\infty}y^{s-1}e^{-y}dy\leq G_{s}(x)$ (11)

where

$\displaystyle G_{s}(x)=\frac{(x+b_{s})^{s}-x^{s}}{sb_{s}}e^{-x}$ (12)

with $b_{s}=\Gamma(s+1)^{\frac{1}{s-1}}$. Then, we can bound (10) by

$\displaystyle\int_{r\log(n)-\log(C)}^{\infty}\frac{y^{s-1}e^{-y}}{\Gamma(s)}dy$
$\displaystyle\leq\frac{G_{s}(\log(\frac{n^{r}}{C}))}{\Gamma(s)}$
$\displaystyle=\frac{\left(\log\left(\frac{n^{r}}{C}\right)+\Gamma(s+1)^{\frac{1}{s-1}}\right)^{s}-\left(\log(\frac{n^{r}}{C})\right)^{s}}{\Gamma(s)s\Gamma(s+1)^{\frac{1}{s-1}}}\cdot\frac{C}{n^{r}}.$
(13)

The last expression in (13) is summable, for $n=1,\dotsc,+\infty$, yielding
$\sum_{n=1}^{\infty}\textrm{Pr}(E_{n})<\infty$. Thus, by the Borel–Cantelli
lemma,

$\displaystyle\textnormal{Pr}\left(E_{n}\enspace i.o.\right)$
$\displaystyle=\textnormal{Pr}\left(\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}E_{k}\right)$
$\displaystyle=\lim_{n\to\infty}\textnormal{Pr}\left(\bigcup_{k=n}^{\infty}E_{k}\right)$
$\displaystyle\leq\lim_{n\to\infty}\sum_{k=n}^{\infty}\textnormal{Pr}\left(E_{k}\right)$
$\displaystyle=0.$

∎

######  Remark 1 (Results in earlier literature).

Lemma 4.1 in [37] states a stronger result, namely that

$\displaystyle\textnormal{Pr}\left(\min_{1\leq i\leq
n}\prod_{j=1}^{s}t_{i}^{j}\leq C\cdot n^{-r}\enspace i.o.\right)=0.$ (14)

However, the additional condition

$\displaystyle\min_{1\leq i\leq n}$

leads to a conclusion that is different from the statement (14). Specifically,
if we assume $\\{t_{i}\\}_{i=1}^{n}$ are independently and identically
distributed (i.i.d.) $U[0,1]^{s}$ and define the event $F_{n}=\\{\min_{1\leq
i\leq n}\prod_{j=1}^{s}t_{i}^{j}\leq n^{-r}\\}$. We choose $C=1$ for
simplicity. Then,

$\displaystyle\textnormal{Pr}\left(F_{n}\enspace i.o.\right)$
$\displaystyle=\textnormal{Pr}\left(\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}F_{k}\right)$
$\displaystyle=\lim_{n\to\infty}\textnormal{Pr}\left(\bigcup_{k=n}^{\infty}F_{k}\right)$
$\displaystyle\geq\lim_{n\to\infty}\textnormal{Pr}\left(F_{n}\right)$
$\displaystyle=\lim_{n\to\infty}\left(1-\prod_{i=1}^{n}\textnormal{Pr}\left(\prod_{j=1}^{s}t_{i}^{j}>n^{-r}\right)\right).$

For example, when $s=1$,

$\displaystyle\lim_{n\to\infty}\left(1-\prod_{i=1}^{n}\textnormal{Pr}\left(t_{i}>n^{-r}\right)\right)$
$\displaystyle=\lim_{n\to\infty}\left(1-\left(1-n^{-r}\right)^{n}\right)$
$\displaystyle=1-\exp\left(-\frac{1}{r}\right).$

######  Remark 2 (The nonunique constant $C$).

The constant $C$ in Lemma 1 is not determined uniquely. Nevertheless, we
retain this constant $C$ in the estimation model.

The event $E_{n}$, as defined in the above proof, only occurs finitely many
times almost surely. The origin is bounded away in this sense. Through
symmetry arguments, we can similarly show that all the $2^{s}$ corners are
bounded away from the following hyperbolic set:

$\displaystyle K_{n,s}=\left\\{\mathbf{t}\in[0,1]^{s}\mid\prod_{1\leq j\leq
s}\min(t_{j},1-t_{j})\geq Cn^{-1}\right\\},$ (15)

where we use the boundary case $n^{-1}$ instead of $n^{-r}$ from Lemma 1. We
also set $Cn^{-1}\leq 1$ to exclude trivial cases.

Figure 1: Illustration of the hyperbolic set $K_{n,s}$ introduced in (15)
using the linear scale (left) and log scale (right) for various $n$, where
$C=1$.

From Lemma 1, we see that the uniform distribution samples “diffuse” to the
corner at a certain rate.

Figure 2 plots the $n$ samples from the uniform distribution $U[0,1]^{2}$ and
the corresponding reference boundaries $t_{1}t_{2}=n^{-1}$. The samples
approach the corner when size increases and a small portion of samples lie
outside their corresponding reference boundaries.

Figure 2: Samples from a two-dimensional uniform distribution $U[0,1]^{2}$
with various sample sizes $n$, and the corresponding reference boundaries
$t_{1}t_{2}=n^{-1}$.

###  2.2 Integrand with infinite variation

We are often interested in integrands with infinite variation, for which the
Koksma–Hlawka inequality (6) is not informative. However, cases remain where
the QMC methods work. The study [19] considered discontinuities inside the
integration domain and proved that the convergence rate is still superior to
the Monte Carlo method as long as some axes are parallel to the discontinuity
interface. Another study [37] explored integrands that blow up at the boundary
of $[0,1]^{s}$ and derived convergence rates based on a boundary growth
condition. We are interested in the latter case. Following the work [37], the
set $K_{n,s}$ was introduced to split the integration domain and the Sobol’
low-variation extension [45] $\tilde{g}$ extends $g$ from $K_{n,s}$ to
$[0,1]^{s}$. The extension $\tilde{g}$ depends on $n$ and $s$ through the
hyperbolic set $K_{n,s}$, but we suppress these dependences for the sake of
simpler notations. Though the original name in the literatrue states “low-
variation”, we assure the readers that the extension $\tilde{g}$ indeed has
finite variation. For completeness, we provide the details here.

To this end, some notations need to be introduced. For an index set
$\mathfrak{u}\subseteq\\{1,2,\dotsc,s\\}$, $-\mathfrak{u}$ denotes its
complement $\\{1,2,\dotsc,s\\}\setminus\mathfrak{u}$. The notation
$\mathbf{z^{\mathfrak{u}}}:\mathbf{c^{-\mathfrak{u}}}$ is used to denote the
point $\mathbf{y}\in[0,1]^{s}$, where $y^{j}=z^{j}$ for $j\in\mathfrak{u}$ and
$y^{j}=c^{j}$ for $j\notin\mathfrak{u}$, “concatenating” the vectors
$\mathbf{z^{\mathfrak{u}}}$ and $\mathbf{c^{\mathfrak{-u}}}$. For simplicity,
we assume the derivative of the integrand $\partial^{\mathfrak{u}}g$ exists in
$(0,1)^{s}$, for $\mathfrak{u}\subseteq\\{1,2,\dotsc,s\\}$. Given an anchor
point $\mathbf{c}\in K_{n,s}$ and using the fundamental theorem of calculus,
we obtain

$\displaystyle
g(\mathbf{t})=g(\mathbf{c})+\sum_{\mathfrak{u}\neq\emptyset}\int_{[\mathbf{c^{\mathfrak{u}}},\mathbf{t^{\mathfrak{u}}}]}\partial^{\mathfrak{u}}g(\mathbf{z^{\mathfrak{u}}}:\mathbf{c^{-\mathfrak{u}}})d\mathbf{z^{\mathfrak{u}}}.$
(16)

The Sobol’ low-variation extension $\tilde{g}:\mathbb{R}^{s}\to\mathbb{R}$ is
given by

$\displaystyle\tilde{g}(\mathbf{t})=g(\mathbf{c})+\sum_{\mathfrak{u}\neq\emptyset}\int_{[\mathbf{c^{\mathfrak{u}}},\mathbf{t^{\mathfrak{u}}}]}\mathbbm{1}_{\\{\mathbf{z^{\mathfrak{u}}}:\mathbf{c^{-\mathfrak{u}}}\in
K_{n,s}\\}}~{}\partial^{\mathfrak{u}}g(\mathbf{z^{\mathfrak{u}}}:\mathbf{c^{-\mathfrak{u}}})d\mathbf{z^{\mathfrak{u}}}.$
(17)

Following [37, 45], we apply a three-epsilon argument to bound the integration
error. Recalling the notations on the exact integration (1) and QMC estimation
(2), we bound:

$\displaystyle\begin{split}\lvert I(g)-\hat{I}_{n}(g)\rvert&\leq\lvert
I(g)-I(\tilde{g})\rvert+\lvert
I(\tilde{g})-\hat{I}_{n}(\tilde{g})\rvert+\lvert\hat{I}_{n}(\tilde{g})-\hat{I}_{n}(g)\rvert\\\
&\leq\int_{[0,1]^{s}-K_{n,s}}\lvert
g-\tilde{g}\rvert+\Delta^{*}(\mathbf{t}_{1},\dotsc,\mathbf{t}_{n})V_{HK}(\tilde{g})+\frac{1}{n}\sum_{i=1}^{n}\lvert\tilde{g}(\mathbf{t}_{i})-g(\mathbf{t}_{i})\rvert.\end{split}$
(18)

Observe that $g$ and $\tilde{g}$ coincide in $K_{n,s}$. Moreover, we have

$\displaystyle\mathbb{E}\left[\frac{1}{n}\sum_{i=1}^{n}\left\lvert\tilde{g}(\mathbf{t}_{i})-g(\mathbf{t}_{i})\right\rvert\right]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[\left\lvert\tilde{g}(\mathbf{t}_{i})-g(\mathbf{t}_{i})\right\rvert\right]$
$\displaystyle=\int_{[0,1]^{s}-K_{n,s}}\lvert g-\tilde{g}\rvert.$ (19)

Using the above inequalities, we obtain the following finite upper bound for
the RQMC method:

$\displaystyle\begin{split}\mathbb{E}[\lvert I(g)-\hat{I}_{n}(g)\rvert]&\leq
2\int_{[0,1]^{s}-K_{n,s}}\lvert
g-\tilde{g}\rvert+\mathbb{E}\left[\Delta^{*}(\mathbf{t}_{1},\dotsc,\mathbf{t}_{n})\right]V_{HK}(\tilde{g}).\end{split}$
(20)

Notice that the above bound (20) depends on the value of $C$, which is used to
construct the set $K_{n,s}$ in Equation (15). When $C=0$, this becomes the
classical Koksma–Hlawka inequality, and when $C=\frac{n}{2^{s}}$, the
hyperbolic set $K_{n,s}$ vanishes and the right hand side becomes
$2\int_{[0,1]^{s}}\lvert g\rvert$. For $0<C<\frac{n}{2^{s}}$, the asymptotic
upper bound also reduces to the Koksma–Hlawka inequality, as $n\to\infty$.

To examine the behavior of the difference $g-\tilde{g}$ and the variation
$V_{HK}(\tilde{g})$ as the set $K_{n,s}$ grows when $n$ increases, we assume
the following boundary growth condition (in the nonasymptotic case) outside
the set $K_{n,s}$.

######  Assumption 1 (Nonasymptotic boundary growth condition).

For the integrand $g$, we assume the following boundary growth condition (21)

$\displaystyle\lvert\partial^{\mathfrak{u}}g(\mathbf{t})\rvert\leq
B(\mathbf{t_{c}})\prod_{j=1}^{s}{t_{j}}^{-A_{j}(t_{c}^{j})-1_{j\in\mathfrak{u}}}\quad\mathbf{0}<\mathbf{t}\leq\mathbf{t_{c}},\quad\mathbf{t_{c}}\in\partial
K_{n,s},$ (21)

for a given constant $B(\mathbf{t_{c}})$,
$\forall\mathfrak{u}\subseteq\\{1,2,\dotsc,s\\}$. The inequality
$\mathbf{0}<\mathbf{t}\leq\mathbf{t_{c}}$ holds component-wise.

In Assumption 1, we emphasize the dependence $A_{j}=A_{j}(\mathbf{t_{c}})$,
which is crucial in the following nonasymptotic analysis. For now, we consider
an anisotropic integrand.

######  Example 1 (An anisotropic integrand).

$g(\mathbf{t})=\exp(2\Phi^{-1}(1-t_{1})+\Phi^{-1}(1-t_{2}))\quad\mathbf{t}\in[0,1]^{2}.$
(22)

As is often the case, the function 22 in our example is not isotropic w.r.t.
the coordinates. Figure 3 provides an illustration of the sets $K_{n,s}$ with
$s=2$ and the contour of the function 22 on the log scale. The value of
$\lvert\partial^{\mathfrak{u}}g(\mathbf{t})\rvert$ may not be constant for
$\mathbf{t}\in\partial K_{n,s}$.

Figure 3: Illustration of the set $K_{n,s}$ in linear scale with $C=1$ and a
contour plot of the function
$g(\mathbf{t})=\exp(2\Phi^{-1}(1-t_{1})+\Phi^{-1}(1-t_{2}))$ in log scale.

This finding motivates the determination of an upper bound of the following
term:

$\displaystyle\prod_{j=1}^{s}{t_{c}^{j}}^{-A_{j}(t_{c}^{j})},~{}\textrm{where}~{}\mathbf{t_{c}}\in\partial
K_{n,s}$

in dimension $s>1$, which is equivalent to the following optimization problem:

$\displaystyle\begin{split}\max_{\mathbf{t_{c}}}&~{}\sum_{j=1}^{s}-A_{j}\log
t_{c}^{j}\\\ s.t.&~{}\sum_{j=1}^{s}\log t_{c}^{j}=\log\delta,\end{split}$ (23)

for a given $\delta=Cn^{-1}>0$. We define $\mathbf{t_{c}^{*}}$ as the
optimizer of (23) and

$A_{j}^{*}=A_{j}(\mathbf{t_{c}^{*}}),$ (24)

for $j=1,\dotsc,s$. In one dimension, the point $\mathbf{t_{c}^{*}}$ is
explicitly given by $Cn^{-1}$. For $s\geq 2$, and $\mathbf{t_{c}^{*}}$ is
determined by solving the optimization problem (23).

Following the proof of Lemma 5.1 in [37], for an integrand $g$ satisfying the
boundary growth condition (21), the difference of the integrand $g$ and its
Sobol’ low-variation extension $\tilde{g}$ is bounded by

$\displaystyle\lvert
g(\mathbf{t})-\tilde{g}(\mathbf{t})\rvert\leq\tilde{B}(\mathbf{t_{c}})\prod_{j=1}^{s}{t_{j}}^{-A_{j}(t_{c}^{j})}\quad
0\leq\mathbf{t}\leq\mathbf{t_{c}},\quad\mathbf{t_{c}}\in\partial K_{n,s},$
(25)

where
$\tilde{B}(\mathbf{t_{c}})=B\prod_{j=1}^{s}(1+\frac{1}{A_{j}(t_{c}^{j})})$.
Notice that the power $-A_{j}$ has a dependence on the specific choice of
$\mathbf{t_{c}}\in K_{n,s}$. Consider the optimization problem (23), we
further consider the bound of $\lvert g-\tilde{g}\rvert$ with the only
dependence on $K_{n,s}$ in Lemma 2.

######  Lemma 2 (A bound for the difference of the integrand and its Sobol’
low-variation extension).

Let $\mathbf{t}\notin K_{n,s}$, for integrand $g$ satisfying the boundary
growth condition (21), we have

$\displaystyle\lvert
g(\mathbf{t})-\tilde{g}(\mathbf{t})\rvert\leq\tilde{B}(\mathbf{t_{c}^{*}})\prod_{j=1}^{s}{t_{j}}^{-A_{j}^{*}}$
(26)

where $A_{j}^{*}$ is determined by the optimization problem (23).

###### Proof.

Directly using Equation (23) and (25) leads to the result. ∎

From Lemma 2, following Lemma 5.4 and the proof of Theorem 5.5 in [37], we
have

$\displaystyle\int_{[0,1]^{s}-K_{n,s}}\lvert g-\tilde{g}\rvert\leq
B_{1}n^{\max_{j}A_{j}^{*}-1},$ (27)

and

$\displaystyle V_{HK}(\tilde{g})\leq B_{2}n^{\max_{j}A_{j}^{*}},$ (28)

for finite $B_{1}$ and $B_{2}$. Now, we have obtained the ingredients to
present the nonasymptotic upper bound for the RQMC upper bound.

######  Theorem 3 (Nonasymptotic QMC error estimate).

For a continuous integrand $g$ on $[0,1]^{s}$, whose mixed first-order
derivative exists, the expected integration for the RQMC method with a
quadrature size of $n$ satisfies

$\displaystyle\mathbb{E}[\lvert I(g)-\hat{I}_{n}(g)\rvert]\leq
C_{1}n^{-1+\max_{j}A_{j}^{*}}+C_{2}n^{-1+\epsilon+\max_{j}A_{j}^{*}}$ (29)

where $\epsilon>0$, and $C_{1}$ and $C_{2}$ are finite constants that depend
on $g$ and $s$.

###### Proof.

Substitute Equation (27) and Equation (28) into (20) and notice that, for low-
discrepancy sequences,

$\displaystyle\mathbb{E}\left[\Delta^{*}(\mathbf{t}_{1},\dotsc,\mathbf{t}_{n})\right]=\mathcal{O}(n^{-1}(\log
n)^{s})$ (30)

(see [38], for reference). ∎

Theorem 3 is the contribution of this work to the nonasymptotic QMC error
estimate model. We observe that the convergence rate of the RQMC method
depends on $\max_{j}A_{j}^{*}$, which depends on the QMC quadrature size $n$.

######  Remark 3 (Nonasymptotics rather than singularity).

In some cases, for instances, the two examples considered in the next Section
3, both integrands exhibit singularities at the boundary. However, the value
of $A_{j}^{*}$ from (24) converges to 0 as $n\to\infty$. In this situation,
the convergence rate for the RQMC method becomes the optimal asymptotic rate
$\mathcal{O}(n^{-1+\epsilon})$ for $\epsilon>0$, in the presence of the
singularities and hence infinite variations. However, the optimal convergence
rate may not be observed for a finite sample size $n$, as $A_{j}^{*}$ may not
converge to 0, resulting a suboptimal convergence rate indicated by Theorem 3.
The next section illustrates the nonasymptotics with two concrete examples.

##  3 Two examples of infinite variation integrands

This section delves into two significant integration examples: the expectation
estimate of lognormal random variables and the QoIs relating to the solution
of elliptic PDEs with lognormal coefficients. These examples shed light on the
intricate nature of infinite variation integrands.

###  3.1 Lognormal random variable

We first analyze the integration problem (1), where
$g=\exp(\bm{\sigma}^{T}{\Phi}^{-1})$, $\bm{\sigma}\in\mathbb{R}^{s}_{+}$, and
${\Phi}^{-1}:[0,1]^{s}\to\mathbb{R}^{s}$ applies the inverse standard normal
distribution CDF component-wise. The partial derivative of $g$ w.r.t.
$\mathbf{t}_{\mathfrak{u}}$ is given by

$\displaystyle\frac{\partial g}{\partial\mathbf{t}_{\mathfrak{u}}}$
$\displaystyle=\frac{\partial}{\partial\mathbf{t}_{\mathfrak{u}}}\exp\left(\sum_{j=1}^{s}\sigma_{j}\Phi^{-1}(t_{j})\right)$
$\displaystyle=\exp\left(\sum_{j=1}^{s}\sigma_{j}\Phi^{-1}(t_{j})\right)\cdot\prod_{j\in\mathfrak{u}}\sigma_{j}{\partial^{j}\Phi^{-1}(t_{j})}.$
(31)

Next, we consider the case where $t_{j}\to 0$ for $j=1,\dotsc,s$. We use the
following approximation:

$\displaystyle\begin{split}\Phi^{-1}(t_{j})&=-\sqrt{-2\log(t_{j})}+o(1)\quad
t\to 0\\\ &\approx-\sqrt{-2\log(t_{j})}.\end{split}$ (32)

The case where $t_{j}\to 1$ can be derived similarly due to symmetry:

$\displaystyle\Phi^{-1}(t_{j})$
$\displaystyle\approx\sqrt{-2\log(1-t_{j})}\quad t_{j}\to 1$ (33)

(see [42], Chapter 3.9, for instance).

Using the approximation of $\Phi^{-1}$ (33), we can simplify the partial
derivative of $g$ as

$\displaystyle\frac{\partial g}{\partial\mathbf{t}_{\mathfrak{u}}}$
$\displaystyle=\exp\left(\sum_{j=1}^{s}\sigma_{j}\Phi^{-1}(t_{j})\right)\cdot\prod_{j\in\mathfrak{u}}\sigma_{j}{\partial^{j}\Phi^{-1}({t}_{j})}$
$\displaystyle\approx\exp\left(\sum_{j=1}^{s}\sigma_{j}{\sqrt{-2\log(1-{t}_{j})}}\right)\cdot\prod_{j\in\mathfrak{u}}\sigma_{j}{\partial^{j}\sqrt{-2\log(1-{t}_{j})}}$
$\displaystyle=\exp\left(\sum_{j=1}^{s}\sigma_{j}{\sqrt{-2\log(1-{t}_{j})}}\right)\cdot\prod_{j\in\mathfrak{u}}\sigma_{j}\frac{1}{(1-{t}_{j})\sqrt{-2\log(1-{t}_{j})}}.$
(34)

We define $\mathbf{v}=1-\mathbf{t}$ to simplify the notation, and the term
$\frac{1}{(1-{t})\sqrt{-2\log(1-{t})}}$ behaves like
$\mathcal{O}({v}^{-1+\epsilon})$ for any $\epsilon>0$ as ${v}\to 0$. To check
the boundary growth condition (21), it remains to study the local growth
behavior of the term

$\displaystyle{h(\mathbf{v}):=\exp\left(\sum_{j=1}^{s}\sigma_{j}{\sqrt{-2\log{v}_{j}}}\right)}=\exp\left(\sum_{j=1}^{s}\sigma_{j}{\sqrt{-2\log(1-{t}_{j})}}\right).$
(35)

######  Lemma 4 (First-order Taylor approximation).

The function

$\displaystyle
h(\mathbf{v}):=\exp\left(\sum_{j=1}^{s}\sigma_{j}{\sqrt{-2\log{v}_{j}}}\right)$
(36)

satisfies the following bound:

$\displaystyle{h(\mathbf{v})}\leq
B(\mathbf{v})\prod_{j=1}^{s}({v}_{j})^{-A_{j}}\quad{\mathbf{v}}\notin
K_{n,s},$ (37)

where $A_{j}=\frac{\sigma_{j}}{\sqrt{-2\log v_{c}^{j}}}$,
$\mathbf{v_{c}}\in\partial K_{n,s}$.

###### Proof.

We take the logarithm of both sides of Equation (37):

$\displaystyle\log{h(\mathbf{v})}=\sum_{j=1}^{s}\sigma_{j}\sqrt{-2\log({v_{j}})}\leq-\sum_{j=1}^{s}{A}_{j}\log(v_{j})+\log
B(\mathbf{v}).$ (38)

We further apply the change of variable to simplify the notation. For
$\mathbf{v_{c}}\in\partial K_{n,s}$, let $\mathbf{z}=-\log\mathbf{v}$, and
$\mathbf{z_{c}}=-\log\mathbf{v_{c}}$. By a first-order Taylor approximation of
(38),

$\displaystyle\begin{split}\sum_{j=1}^{s}\sigma_{j}\sqrt{2z_{j}}&\approx\nabla_{\mathbf{z}=\mathbf{z_{c}}}\left(\sum_{j=1}^{s}\sigma_{j}\sqrt{2{z_{j}}}\right)\cdot(\mathbf{z}-\mathbf{z_{c}})+\left(\sum_{j=1}^{s}\sigma_{j}\sqrt{2z_{c}^{j}}\right)\\\
&=\sum_{j=1}^{s}\frac{\sigma_{j}}{\sqrt{2z_{c}^{j}}}{z}^{j}+\sum_{j=1}^{s}\sigma_{j}\sqrt{\frac{z_{c}^{j}}{2}}\end{split}$
(39)

where we use $\lvert\mathbf{z}-\mathbf{z_{c}}\rvert<\epsilon$ for a small
threshold $\epsilon>0$. By comparing Equations (37) and (39), we deduce that

$\displaystyle{A}_{j}$
$\displaystyle=\frac{\sigma_{j}}{\sqrt{2z_{c}^{j}}}=\frac{\sigma_{j}}{\sqrt{-2\log
v_{c}^{j}}}=\frac{\sigma_{j}}{\sqrt{-2\log(1-t_{c}^{j})}}.$ (40)

∎

Notice that ${A}_{j}\to 0$ as $t_{j}\to 1$ for $j=1,\dotsc,s$. Theorem 3
indicates that the QMC method achieves the asymptotic convergence rate
$\mathcal{O}(n^{-1+\epsilon})$.

Let us now consider the optimization problem (23):

$\displaystyle\begin{split}\max_{\mathbf{v_{c}}}&~{}\sum_{j=1}^{s}-A_{j}\log
v_{c}^{j}\\\ s.t.&~{}-\sum_{j=1}^{s}\log
v_{c}^{j}=\log\delta^{-1},\end{split}$ (41)

for a given $\delta=Cn^{-1}$. The Lagrangian is given by

$\displaystyle\begin{split}\mathcal{L}&=\sum_{j=1}^{s}-A_{j}\log
v_{c}^{j}+\lambda\left(\sum_{j=1}^{s}-\log v_{c}^{j}-\log\delta^{-1}\right)\\\
&=\sum_{j=1}^{s}\sigma_{j}\sqrt{\frac{-\log
v_{c}^{j}}{2}}+\lambda\left(\sum_{j=1}^{s}-\log
v_{c}^{j}-\log\delta^{-1}\right).\end{split}$ (42)

The optimizer $-\log{v_{c}^{j}}^{*}$ for Equation (41) is given by

$\displaystyle-\log{v_{c}^{j}}^{*}=\frac{\sigma_{j}^{2}}{\sum_{j=1}^{s}\sigma_{j}^{2}}\log\delta^{-1}.$
(43)

Moreover, the optimal value of $A_{j}$, denoted by $A_{j}^{*}$, is given by

$\displaystyle
A_{j}^{*}=\sqrt{\frac{\sum_{j=1}^{s}\sigma_{j}^{2}}{2\log\delta^{-1}}}=\sqrt{\frac{\sum_{j=1}^{s}\sigma_{j}^{2}}{2(\log
n-\log C)}}.$ (44)

Hence, $\max_{j}A_{j}^{*}=\sqrt{\frac{\sum_{j=1}^{s}\sigma_{j}^{2}}{2(\log
n-\log C)}}$. For a finite sample size $n$, considering the convergence model
(29) in Theorem 3, the convergence rate improves as $n$ increases, or the
variance $\sigma_{j}^{2}$ decreases for each dimension $j$.

###  3.2 Elliptic partial differential equations with lognormal coefficients

We consider an elliptic PDE on a Lipschitz domain
$\mathcal{D}\in\mathbb{R}^{d}$, in the following form:

$\displaystyle-\nabla_{\mathbf{x}}\cdot\left(a(\mathbf{x};\omega)\nabla_{\mathbf{x}}u(\mathbf{x};\omega)\right)$
$\displaystyle=f(\mathbf{x};\omega)$ for $\mathbf{x}\in\mathcal{D}$, (45a)
$\displaystyle u(\mathbf{x};\omega)$ $\displaystyle=0$ $\displaystyle\text{for
$\mathbf{x}\in\partial\mathcal{D}$},$ (45b)

with almost all $\omega\in\Omega$, where $\Omega$ corresponds to a complete
probability space tuple $(\Omega,\mathcal{F},\mathbb{P})$. The differential
operators $\nabla\cdot$ and $\nabla$ are taken w.r.t. the spatial variable
$\mathbf{x}$, and $\partial_{n}$ is the outward normal derivative operator. We
explore the equations in (45) in the weak form. For a suitable functional
space $V$ (e.g., $V=H^{1}_{0}(\mathcal{D})$), we seek $u(\mathbf{x};\omega)\in
V$ such that

$\begin{split}\left(a(\mathbf{x};\omega)\nabla_{\mathbf{x}}u(\mathbf{x};\omega),\nabla_{\mathbf{x}}v(\mathbf{x};\omega)\right)&=\left(f(\mathbf{x};\omega),v(\mathbf{x};\omega)\right)\quad\\\
&\text{for all }v\in V\text{ and almost all }\omega\in\Omega.\end{split}$ (46)

The QoI $Q(\omega)$ is given in the following general form:

$\displaystyle Q(\mathbf{y}(\omega))=\mathcal{G}(u(\cdot;\omega)),$ (47)

where $\mathcal{G}\in V^{\prime}$, $V^{\prime}$ denotes the dual space of $V$.
We are interested in computing $\mathbb{E}[Q]$. The integrand in (1) is given
by $g=Q\circ{\Phi}^{-1}$.

We consider $a(\mathbf{x};\omega)$, which takes the following form:

$\displaystyle
a(\mathbf{x};\omega)=\exp\left(\sum_{j=1}^{s}\mathrm{y}_{j}\psi_{j}\right),$
(48)

with $s\in\mathbb{N}^{+}$ and $y_{j}$ are i.i.d. samples from
$\mathcal{N}(0,1)$, $j=1,\dotsc,s$. Moreover, we define
$b_{j}=\lVert\psi_{j}\rVert_{L^{\infty}(\mathcal{D})}<+\infty$, for
$j=1,\dotsc,s$.

From [16], an upper bound for the derivative of $Q$ w.r.t. the random variable
$\mathbf{y}$ is given by

$\displaystyle\left\lvert{\partial^{\mathfrak{u}}}Q({\mathbf{y}})\right\rvert$
$\displaystyle\leq\lVert\mathcal{G}\rVert_{V^{\prime}}\lVert\partial^{\mathfrak{u}}u(\cdot,\mathbf{y})\rVert_{V}$
$\displaystyle\leq\frac{\lvert\mathfrak{u}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}\rvert}}\left(\prod_{j\in\mathfrak{u}}b_{j}\right)\underbrace{\lVert
f\rVert_{V^{\prime}}\lVert\mathcal{G}\rVert_{V^{\prime}}}_{K^{*}}\prod_{j=1}^{s}\exp(b_{j}\lvert{\mathrm{y}}_{j}\rvert)$
(49)

The derivative of the integrand $g$ is as follows:

$\displaystyle\begin{split}\left\lvert{\partial}^{{\mathfrak{u}}}g(\bm{t})\right\rvert&=\left\lvert\partial^{\mathfrak{u}}Q(\Phi^{-1}(\mathbf{t}))\right\rvert=\lvert\partial^{\mathfrak{u}}Q\rvert\cdot\prod_{j\in\mathfrak{u}}\left\lvert\partial^{j}\Phi^{-1}(t_{j})\right\rvert\\\
&\leq K^{*}\frac{\lvert\mathfrak{u}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}\rvert}}\left(\prod_{j\in\mathfrak{u}}b_{j}\right)\cdot\prod_{j=1}^{s}\exp(b_{j}\lvert\mathrm{y}_{j}\rvert)\left\lvert\partial^{j}\Phi^{-1}({t}_{j})\right\rvert\\\
&\approx K^{*}\frac{\lvert\mathfrak{u}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}\rvert}}\left(\prod_{j\in\mathfrak{u}}b_{j}\right)\prod_{j=1}^{s}b_{j}\exp\left(b_{j}{\sqrt{-2\log(t_{j})}}\right)\cdot\prod_{j=1}^{s}\frac{1}{t_{j}\sqrt{-2\log(t_{j})}}\end{split}$
(50)

where we apply the approximation of $\Phi^{-1}$ (32) as $t_{j}$ approaches 0,
where $j=1,\dotsc,s$. Unlike the example in Section 3.1, the considered
integrand $g$ has singularities when $t_{j}\to 0$ and $t_{j}\to 1$. Due to the
symmetry in each dimension, we only need to consider the singularity at 0. To
investigate the local growth behavior of the derivative, we define

$\displaystyle h(\mathbf{t})=K^{*}\frac{\lvert\mathfrak{u}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}\rvert}}\left(\prod_{j\in\mathfrak{u}}b_{j}\right)\prod_{j=1}^{s}b_{j}\exp\left(b_{j}{\sqrt{-2\log(t_{j})}}\right).$
(51)

Similar to the analysis in Section. 3.1, we apply a first-order Taylor
approximation. We take the logarithm of both sides of (51):

$\displaystyle\log\left\lvert
h(\mathbf{t})\right\rvert=\sum_{j=1}^{s}b_{j}\sqrt{-2\log({t_{j}})}+\log\left(K^{*}\frac{\lvert\mathfrak{u}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}\rvert}}\left(\prod_{j\in\mathfrak{u}}b_{j}\right)\prod_{j=1}^{s}b_{j}\right).$
(52)

Thus, for

$\displaystyle{A}_{j}$ $\displaystyle=\frac{b_{j}}{\sqrt{-2\log t_{c}^{j}}},$
(53)

we have

$\displaystyle\left\lvert h(\mathbf{t})\right\rvert\leq
B(\mathbf{t})\prod_{j=1}^{s}({t}_{j})^{-A_{j}}\quad{\mathbf{t}}\notin K_{n,s}$
(54)

where

$\displaystyle B(\mathbf{t})=\left(K^{*}\frac{\lvert\mathfrak{u}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}\rvert}}\left(\prod_{j\in\mathfrak{u}}b_{j}\right)\prod_{j=1}^{s}b_{j}\right)\cdot\exp\left(\sum_{j=1}^{s}b_{j}\sqrt{\frac{-\log
t_{c}^{j}}{2}}\right).$

Similarly, we can determine $A_{j}^{*}$ as follows:

$\displaystyle A_{j}^{*}=\sqrt{\frac{\sum_{j=1}^{s}b_{j}^{2}}{2(\log n-\log
C)}}$ (55)

and $\max_{j}A_{j}^{*}=\sqrt{\frac{\sum_{j=1}^{s}b_{j}^{2}}{2(\log n-\log
C)}}$.

In this section, we have considered two concrete examples with unbounded
variation and applied the theory developed in Section 2 to analyze their
convergence rates. Consistent with the theory, the convergence rate can be
improved by increasing the sample size $n$. However, whether we can further
enhance the convergence rate by modifying the sample distribution and placing
more points near the singular corners is unclear. In the next section, we
discuss the effects of importance sampling.

##  4 Importance sampling

This section provides two kinds of IS proposal distributions and analyzes
their effects on the convergence rate.

###  4.1 First proposal distribution

This section proposes a Gaussian distribution with scaled variance to
distribute more samples closer to the corners. This kind of IS was studied in
[27] for multivariate functions that belong to certain Sobolev spaces, where
the worst-case integration error was analyzed and optimized.

Let
$g\in\\{\exp\circ(\bm{\sigma}^{T}\Phi^{-1}),Q\circ\Phi^{-1}\\}:[0,1]^{s}\to\mathbb{R}$,
which are the two integrands considered in Section 3. We introduce the
component-wise multiplication notation $\odot$, such that
$\bm{\alpha}\odot\mathbf{y}=\\{\alpha_{1}\mathrm{y}_{1},\alpha_{2}\mathrm{y}_{2},\dotsc,a_{s}\mathrm{y}_{s}\\}=\mathrm{diag}(\bm{\alpha})\mathbf{y}$
when $\bm{\alpha},\mathbf{y}\in\mathbb{R}^{s}$. The integral of $g$ over
$[0,1]^{s}$ is

$\displaystyle\begin{split}I(g)&=\int_{[0,1]^{s}}g(\mathbf{t})d\mathbf{t}\\\
&=\int_{\mathbb{R}^{s}}\nu(\mathbf{y})\rho(\mathbf{y})d\mathbf{y}\\\
&=\prod_{j=1}^{s}{\alpha_{j}}\int_{\mathbb{R}^{s}}\nu(\bm{\alpha}\odot\mathbf{y})\rho({\bm{\alpha}}\odot\mathbf{y})d\mathbf{y}\\\
&=\prod_{j=1}^{s}{\alpha_{j}}\int_{[0,1]^{s}}\nu(\bm{\alpha}\odot\Phi^{-1}(\bm{t}))\cdot\frac{\rho(\bm{\alpha}\odot\Phi^{-1}(\bm{t}))}{\rho(\Phi^{-1}(\bm{t}))}d\bm{t},\end{split}$
(56)

where $\nu=g\circ\Phi:\mathbb{R}^{s}\to\mathbb{R}$,
$\bm{\alpha}\in\mathbb{R}^{s}_{+}$, and $\rho$ is the $s$-dimensional standard
normal distribution density. The integrand with IS $g_{\textrm{IS}}$ is given
by

$\displaystyle\begin{split}g_{\textrm{IS}}(\bm{t})&=\left(\prod_{j=1}^{s}{\alpha_{j}}\right)\nu(\bm{\alpha}\odot\Phi^{-1}(\bm{t}))\cdot\frac{\rho(\bm{\alpha}\odot\Phi^{-1}(\bm{t}))}{\rho(\Phi^{-1}(\bm{t}))}\\\
&=\left(\prod_{j=1}^{s}{\alpha_{j}}\right)\nu(\bar{\mathbf{y}})\cdot\prod_{j=1}^{s}\exp\left(-\frac{\mathrm{y}_{j}^{2}}{2}({\alpha}_{j}^{2}-1)\right)\end{split}$
(57)

where we define $\mathbf{y}=\Phi^{-1}(\bm{t})$ and
$\bar{\mathbf{y}}=\bm{\alpha}\odot\mathbf{y}=\bm{\alpha}\odot\Phi^{-1}(\bm{t})$.
The function $g_{\textrm{IS}}$ no longer has the singularity at the boundaries
for ${\alpha}_{j}>1$, $j=1,2,\dotsc,s$. Figure 4 illustrates the one-
dimensional integrand with and without IS in Example 1.

Figure 4: One-dimensional Integrands $g$ and $g_{\mathrm{IS}}$ in Example 1,
where $s=1,\sigma=1.0$, $\alpha>1$. The IS eliminates the singularity at the
boundary. The legend IS_A stands for the IS of the first type.

Let us define
$\varphi(\bm{t})=\prod_{j=1}^{s}{\alpha_{j}}\exp\left(-\frac{(\Phi^{-1}({t_{j}}))^{2}}{2}({\alpha}_{j}^{2}-1)\right)$
for the simplicity of notations. We will analyze the mixed first-order
derivative ${\partial^{\mathfrak{u}}}g_{\mathrm{IS}}$. Specifically, we have

$\displaystyle\begin{split}\left\lvert{\partial^{\mathfrak{u}}}g_{\mathrm{IS}}(\bm{t})\right\rvert&=\left\lvert\sum_{\mathfrak{z}\subseteq\mathfrak{u}}\frac{\partial}{\partial
t_{\mathfrak{u}-\mathfrak{z}}}\nu(\bar{\mathbf{y}})\cdot{\partial^{\mathfrak{z}}}\varphi(\bm{t})\right\rvert.\end{split}$
(58)

We have

$\displaystyle\frac{\partial}{\partial
t_{\mathfrak{u}-\mathfrak{z}}}\nu(\bar{\mathbf{y}})=\frac{\partial}{\partial\bar{\mathrm{y}}_{\mathfrak{u}-\mathfrak{z}}}\nu(\bar{\mathbf{y}})\cdot\prod_{j\in\mathfrak{u}-\mathfrak{z}}\alpha_{j}\frac{\partial\mathrm{y}_{j}}{\partial
t_{j}}.$ (59)

and

$\displaystyle\begin{split}{\partial^{\mathfrak{z}}}\varphi(\bm{t})&=\prod_{j=1}^{s}{\alpha_{j}}\exp\left(-\frac{(\Phi^{-1}({t_{j}}))^{2}}{2}({\alpha}_{j}^{2}-1)\right)\cdot\prod_{k\in\mathfrak{z}}-\frac{\alpha_{k}^{2}-1}{2}\cdot\frac{\partial\left(\Phi^{-1}(t_{k})\right)^{2}}{\partial
t_{k}}\end{split}$ (60)

Notice that $\binom{\mathfrak{u}}{\mathfrak{z}}=1$, for all
$\mathfrak{z}\preceq\mathfrak{u}$. In Example 1 (i.e., when
$\nu=\exp\circ\bm{\sigma}^{T}$), we have

$\displaystyle\partial^{\mathfrak{u}-\mathfrak{z}}\nu(\bar{\mathrm{y}})=\exp\left(\sum_{j=1}^{s}\sigma_{j}\bar{\mathrm{y}}_{j}\right)\cdot\prod_{k\in\mathfrak{u}-\mathfrak{z}}\sigma_{k}.$
(61)

By combining equations (58) to (61), we have

$\displaystyle\begin{split}\left\lvert{\partial^{\mathfrak{u}}}g_{\mathrm{IS}}(\bm{t})\right\rvert&=\left\lvert\sum_{\mathfrak{z}\subseteq\mathfrak{u}}\frac{\partial}{\partial
t_{\mathfrak{u}-\mathfrak{z}}}\nu(\bar{\mathbf{y}})\cdot{\partial^{\mathfrak{u}}}\varphi(\bm{t})\right\rvert\\\
&\leq\prod_{j=1}^{s}{\alpha_{j}}\exp\left(-\frac{(\Phi^{-1}({t_{j}}))^{2}}{2}({\alpha}_{j}^{2}-1)\right)\cdot\prod_{j=1}^{s}\exp(\sigma_{j}\lvert\alpha_{j}\Phi^{-1}(t_{j})\rvert)\\\
&\cdot\sum_{\mathfrak{z}\preceq\mathfrak{u}}\prod_{j\in\mathfrak{u}-\mathfrak{z}}\sigma_{j}\alpha_{j}\frac{\partial\mathrm{y}_{j}}{\partial
t_{j}}\cdot\prod_{k\in\mathfrak{z}}-\frac{\alpha_{k}^{2}-1}{2}\cdot{\frac{\partial\left(\Phi^{-1}(t_{k})\right)^{2}}{\partial
t_{k}}}\\\
&\approx\prod_{j=1}^{s}\sigma_{j}{\alpha_{j}}t_{j}^{{\alpha}_{j}^{2}-1}\cdot\exp\left(\alpha_{j}\sigma_{j}{\sqrt{-2\log(t_{j})}}\right)\\\
&\cdot\sum_{\mathfrak{z}\preceq\mathfrak{u}}\prod_{j\in\mathfrak{s}-\mathfrak{z}}\sigma_{j}\alpha_{j}\frac{1}{t_{j}\sqrt{-2\log(t_{j})}}\prod_{k\in\mathfrak{z}}\frac{\alpha_{k}^{2}-1}{t_{k}},\end{split}$
(62)

where we have substituted $1-t_{j}$ by $t_{j}$ to focus on the boundary growth
condition at the singularity. Before further analysis, we now analyze the
derivative in Example 2.

In Example 2, $\nu=Q$, we apply the following upper bound of
$\left\lvert\partial^{{\mathfrak{u}-\mathfrak{z}}}Q(\bar{\mathbf{y}})\right\rvert$:

$\displaystyle\left\lvert\partial^{{\mathfrak{u}-\mathfrak{z}}}Q(\bar{\mathbf{y}})\right\rvert\leq\frac{\lvert\mathfrak{u}-\mathfrak{z}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}-\mathfrak{z}\rvert}}\left(\prod_{j\in\mathfrak{u}-\mathfrak{z}}b_{j}\right)K^{*}\prod_{j=1}^{s}\exp(b_{j}\lvert\bar{\mathrm{y}}_{j}\rvert)$
(63)

Combining equations (58) to (60) and (63), we obtain

$\displaystyle\begin{split}\left\lvert{\partial^{\mathfrak{u}}}g_{\mathrm{IS}}(\bm{t})\right\rvert&=\left\lvert\sum_{\mathfrak{z}\subseteq\mathfrak{u}}{\partial^{\mathfrak{u}-\mathfrak{z}}}Q(\bm{t})\cdot{\partial^{\mathfrak{u}}}\varphi(\bm{t})\right\rvert\\\
&\leq
K^{*}\prod_{j=1}^{s}{\alpha_{j}}\exp\left(-\frac{(\Phi^{-1}({t_{j}}))^{2}}{2}({\alpha}_{j}^{2}-1)\right)\cdot\prod_{j=1}^{s}\exp(b_{j}\lvert\alpha_{j}\Phi^{-1}(t_{j})\rvert)\\\
&\cdot\sum_{\mathfrak{z}\preceq\mathfrak{u}}\frac{\lvert\mathfrak{u}-\mathfrak{z}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}-\mathfrak{z}\rvert}}\left(\prod_{j\in\mathfrak{s}-\mathfrak{z}}b_{j}\right)\prod_{j\in\mathfrak{u}-\mathfrak{z}}\alpha_{j}\frac{\partial\mathrm{y}_{j}}{\partial
t_{j}}\cdot\prod_{k\in\mathfrak{z}}-\frac{\alpha_{k}^{2}-1}{2}\cdot{\frac{\partial\left(\Phi^{-1}(t_{k})\right)^{2}}{\partial
t_{k}}}\\\ &\approx
K^{*}\prod_{j=1}^{s}{\alpha_{j}}t_{j}^{{\alpha}_{j}^{2}-1}\cdot\exp\left(\alpha_{j}b_{j}{\sqrt{-2\log(t_{j})}}\right)\\\
&\cdot\sum_{\mathfrak{z}\preceq\mathfrak{u}}\frac{\lvert\mathfrak{u}-\mathfrak{z}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}-\mathfrak{z}\rvert}}\prod_{j\in\mathfrak{s}-\mathfrak{z}}b_{j}\alpha_{j}\frac{1}{t_{j}\sqrt{-2\log(t_{j})}}\prod_{k\in\mathfrak{z}}\frac{\alpha_{k}^{2}-1}{t_{k}},\end{split}$
(64)

where we use the approximation $\Phi^{-1}(t)\approx-\sqrt{-2\log(t)}$ and omit
the term $1/\sqrt{-2\log t_{j}}=o(t_{j}^{-1})$, as $t\to 0$. We observe that
the upper bound (64) shares similar structures with (62). We will analyze the
upper bound (64) in detail and show later that the conclusions will apply to
(62) of Example 1 as well.

Similar to the analysis in Section 3.1, we have

$\displaystyle\exp\left(\alpha_{j}b_{j}{\sqrt{-2\log(t_{j})}}\right)\leq
C(\mathbf{t})\cdot(t_{j})^{-B_{j}}$ (65)

where $B_{j}=\frac{\alpha_{j}b_{j}}{\sqrt{-2\log t_{j}}}$. Notice
$\frac{1}{t\sqrt{-2\log(t)}}=o(t^{-1})$ as $t\to 0$. Thus,

$\displaystyle\left\lvert{\partial^{\mathfrak{u}}}g_{\textrm{IS}}(\mathbf{t})\right\rvert\leq
C\prod_{j=1}^{s}(t_{j})^{-1+\alpha_{j}^{2}-1-B_{j}}.$ (66)

The term $\alpha_{j}^{2}-1-B_{j}$ being greater than or equal to $0$ ensures
the term $\max_{j}A_{j}^{*}\leq 0$ in the convergence model (29), which
happens when

$\displaystyle\alpha_{j}^{2}-\alpha_{j}\frac{b_{j}}{\sqrt{-2\log t_{j}}}\geq
1,$ (67)

which is guaranteed by $\alpha_{j}>1$ and a finite $t_{j}^{-1}$ such that
$t_{j}^{-1}\geq\exp\left(\frac{1}{2}\left(\frac{\alpha_{j}b_{j}}{\alpha_{j}^{2}-1}\right)^{2}\right)$.
For the optimality condition in Example 1 we can simply replace $b_{j}$ in
(67) with $\sigma_{j}$. Finite $t_{j}$ for $j=1,\dotsc,s$ imply finite $n$
(via the hyperbolic set (15)), thus the rate $\mathcal{O}(n^{-1})$ is
achievable with a finite $n$, in contrast to the asymptotic optimal rates in
(44) and (55).

Nevertheless, modifying $\alpha$ changes the variance of the integrand
$g_{\mathrm{IS}}$. Thus, we aim to determine the optimal $\mathbf{\alpha}^{*}$
by minimizing the variance of $g_{\mathrm{IS}}$:

$\displaystyle\textrm{Var}(g_{\mathrm{IS}})=\mathbb{E}[g_{\mathrm{IS}}^{2}]-(\mathbb{E}[g_{\mathrm{IS}}])^{2},$
(68)

which is equivalent to find the minimizer for the second-order moment of $g$,
as $\mathbb{E}[g_{\mathrm{IS}}]=\mathbb{E}[g]$:

$\displaystyle\begin{split}\alpha^{*}&=\operatorname*{arg\,min}_{\alpha}\mathbb{E}[g_{\mathrm{IS}}^{2}(\mathbf{t})]\\\
&=\operatorname*{arg\,min}_{\alpha}\int_{[0,1]^{s}}\left(\prod_{j=1}^{s}{\alpha_{j}}\right)^{2}\nu^{2}(\bm{\alpha}\odot\Phi^{-1}(\bm{t}))\cdot\frac{\rho^{2}(\bm{\alpha}\odot\Phi^{-1}(\bm{t}))}{\rho^{2}(\Phi^{-1}(\bm{t}))}d\mathbf{t}\\\
&=\operatorname*{arg\,min}_{\alpha}{\int_{\mathbb{R}^{s}}\left(\prod_{j=1}^{s}{\alpha_{j}}\right)^{2}\nu^{2}(\bm{\alpha}\odot\mathbf{y})\cdot\frac{\rho^{2}(\bm{\alpha}\odot\mathbf{y})}{\rho^{2}(\mathbf{y})}\rho(\mathbf{y})d\mathbf{y}}\\\
&=\operatorname*{arg\,min}_{\alpha}{\int_{\mathbb{R}^{s}}\left(\prod_{j=1}^{s}{\alpha_{j}}\right)\nu^{2}(\mathbf{y})\cdot\frac{\rho^{2}(\mathbf{y})}{\rho^{2}(\bm{\frac{1}{\alpha}}\odot\mathbf{y})}\rho(\bm{\frac{1}{\alpha}}\odot\mathbf{y})d\mathbf{y}}\\\
&=\operatorname*{arg\,min}_{\alpha}{\int_{[0,1]^{s}}\left(\prod_{j=1}^{s}{\alpha_{j}}\right)\nu^{2}(\mathbf{y})\cdot\frac{\rho(\mathbf{y})}{\rho(\bm{\frac{1}{\alpha}}\odot\mathbf{y})}d\bm{t}}.\end{split}$
(69)

When the analytic solution is unavailable, we aim to find an approximate
optimizer $\bar{\alpha}$ using $n$ samples. Specifically, we minimize the
following objective function

$\displaystyle\bar{\alpha}=\operatorname*{arg\,min}_{\alpha}\left(\prod_{j=1}^{s}{\alpha_{j}}\right)\frac{1}{n}\sum_{i=1}^{n}Q^{2}(\mathbf{y}_{i})\cdot\frac{\rho(\mathbf{y}_{i})}{\rho(\bm{\frac{1}{\alpha}}\odot\mathbf{y}_{i})}.$
(70)

In practice, the number $n$ in the pilot run does not exceed the number of
simulations samples.

###  4.2 Second proposal distribution

This section considers another IS proposal distribution. We will only analyze
the Example 2 for brevity, as the case for Example 1 can be derived similarly.
Inspired by the beta distribution, we propose the following distribution
$\rho_{\beta}:[0,1]\to\mathbb{R}_{+}$:

$\displaystyle\rho_{\beta}({w})=\begin{cases}Ct^{\beta-1}&\text{$0\leq
t<\frac{1}{2}$}\\\ C(1-t)^{\beta-1}&\text{$\frac{1}{2}\leq t\leq
1$},\end{cases}$ (71)

with the constant $C=\beta\cdot(\frac{1}{2})^{1-\beta}$. The reason to propose
such a distribution is the ease of computing the CDF
$\Phi_{{\beta}}:[0,1]\to[0,1]$ and its inverse
$\Phi_{{\beta}}^{-1}:[0,1]\to[0,1]$:

$\displaystyle\Phi_{{\beta}}({t})=\begin{cases}\frac{C}{\beta}t^{\beta}&\text{$0\leq
t<\frac{1}{2}$}\\\ 1-\Phi_{{\beta}}(1-{t})&\text{$\frac{1}{2}\leq t\leq
1$},\end{cases}$ (72) $\displaystyle
t=\Phi^{-1}_{\beta}(w)=(2w)^{\frac{1}{\beta}}\cdot\frac{1}{2}\quad\textrm{where
}0\leq w<1/2.$ (73)

Next, we apply IS to compute the integral $I(g)$. Using the density
$\rho_{\beta}$, we can write

$\displaystyle\begin{split}I(g)&=\int_{[0,1]^{s}}g(\mathbf{t})d\mathbf{t}\\\
&=\int_{[0,1]^{s}}g(\mathbf{t})\cdot\frac{\rho_{\beta}(\bm{t})}{\rho_{\beta}(\bm{t})}d\bm{t}\\\
&=\int_{[0,1]^{s}}g(\Phi_{{\beta}}^{-1}(\mathbf{w}))\cdot\frac{1}{\rho_{\beta}(\Phi_{{\beta}}^{-1}(\mathbf{w}))}d\mathbf{w}.\end{split}$
(74)

In Example 2, $g=Q\circ\Phi^{-1}$, and we introduce

$\displaystyle
g_{\mathrm{IS}}(\mathbf{w})=Q\left(\Phi^{-1}(\Phi_{{\beta}}^{-1}(\mathbf{w}))\right)\cdot\frac{1}{\rho_{\beta}(\Phi_{{\beta}}^{-1}(\mathbf{w}))}.$
(75)

Similar to the last section, a possible choice of the parameter $\beta$ can be
determined by minimizing the second-order moment of $g_{\mathrm{IS}}$:

$\displaystyle\begin{split}\beta^{*}&=\operatorname*{arg\,min}_{\beta}\mathbb{E}[g_{\mathrm{IS}}^{2}(\mathbf{w})]\\\
&=\operatorname*{arg\,min}_{\beta}\int_{[0,1]^{s}}\left(\frac{g(\Phi^{-1}(\bm{t}))}{\rho_{\beta}(\mathbf{w})}\right)^{2}\cdot\rho_{\beta}(\mathbf{w})d\mathbf{w}\\\
&=\operatorname*{arg\,min}_{\beta}\int_{[0,1]^{s}}\frac{g^{2}(\Phi^{-1}(\mathbf{w}))}{\rho_{\beta}(\mathbf{w})}d\mathbf{w}.\end{split}$
(76)

Again, we seek an optimizer based on ensembles when the analytical solution is
unavailable.

Next, we study the mixed first-order derivative of $g_{\textrm{IS}}$:

$\displaystyle\begin{split}\left\lvert\frac{\partial}{\partial
w_{\mathfrak{u}}}g_{\mathrm{IS}}(\mathbf{w})\right\rvert&=\left\lvert\sum_{\mathfrak{z}\subseteq\mathfrak{u}}\frac{\partial}{\partial
w_{\mathfrak{u}-\mathfrak{z}}}Q\left(\Phi^{-1}(\Phi_{{\beta}}^{-1}(\mathbf{w}))\right)\cdot\frac{\partial}{\partial
w_{\mathfrak{z}}}\varphi(\mathbf{w})\right\rvert,\\\ \end{split}$ (77)

where
$\varphi(\mathbf{w})=\frac{1}{\rho_{\beta}(\Phi_{{\beta}}^{-1}(\mathbf{w}))}$.
We have

$\displaystyle\frac{\partial}{\partial\mathrm{w}_{\mathfrak{u}-\mathfrak{z}}}Q(\mathbf{w})=\frac{\partial}{\partial\bar{\mathrm{y}}_{\mathfrak{u}-\mathfrak{z}}}Q(\bar{\mathbf{y}})\cdot\prod_{j\in\mathfrak{u}-\mathfrak{z}}\frac{\partial\bar{\mathrm{y}}_{j}}{\partial
t_{j}}\frac{\partial t_{j}}{\partial w_{j}}$ (78)

and

$\displaystyle\begin{split}\frac{\partial}{\partial
u_{\mathfrak{z}}}\varphi(\mathbf{w})&=\frac{\partial}{\partial
u_{\mathfrak{z}}}\prod_{k=1}^{s}\frac{1}{\rho_{\beta_{j}}(\Phi_{{\beta}_{j}}^{-1}(\mathbf{w}))}\\\
&=\prod_{k=1}^{s}\frac{1}{\beta_{k}}(2w_{k})^{\frac{1}{\beta_{k}}-1}\cdot\prod_{j\in\mathfrak{z}}\frac{1-\beta_{j}}{\beta_{j}}w_{j}^{-1},\end{split}$
(79)

where we define $\bar{\mathbf{y}}=\Phi^{-1}(\mathbf{t})$ and
$\mathbf{t}=\Phi^{-1}_{\beta}(\mathbf{w})$. We only consider the case where
$0\leq w_{j}\leq 1/2,\forall j=1,\dotsc,s$ for simplicity, as other cases can
be induced by symmetry. We have

$\displaystyle\frac{\partial t_{j}}{\partial w_{j}}$
$\displaystyle=\frac{1}{\beta_{j}}(2w_{j})^{\frac{1}{\beta_{j}}-1},$ (80)

and

$\displaystyle\begin{split}\frac{\partial\bar{\mathrm{y}}_{j}}{\partial
t_{j}}&=\frac{\partial\Phi^{-1}(t_{j})}{\partial t_{j}}\\\
&\approx-\frac{\partial\sqrt{-2\log t_{j}}}{\partial t_{j}}\\\
&=\frac{1}{t_{j}\sqrt{-2\log t_{j}}}.\end{split}$ (81)

Combining (63) and (78) to (81), we obtain

$\displaystyle\begin{split}\left\lvert\frac{\partial}{\partial
t_{\mathfrak{u}}}g_{\mathrm{IS}}(\mathbf{w})\right\rvert&=\left\lvert\sum_{\mathfrak{z}\subseteq\mathfrak{u}}{\partial^{{\mathfrak{u}-\mathfrak{z}}}}Q\left(\Phi^{-1}(\Phi_{{\beta}}^{-1}(\mathbf{w}))\right)\cdot{\partial^{\mathfrak{z}}}\varphi(\mathbf{w})\right\rvert\\\
&=\left\lvert\sum_{\mathfrak{z}\subseteq\mathfrak{u}}\frac{\partial}{\partial\bar{\mathbf{y}}_{\mathfrak{u}-\mathfrak{z}}}Q(\bar{\mathbf{y}})\cdot\prod_{j\in\mathfrak{u}-\mathfrak{z}}\frac{\partial\bar{\mathrm{y}}_{j}}{\partial
t_{j}}\frac{\partial t_{j}}{\partial w_{j}}\cdot\frac{\partial}{\partial
w_{\mathfrak{z}}}\varphi(\mathbf{w})\right\rvert\\\
&\leq\left\lvert\sum_{\mathfrak{z}\subseteq\mathfrak{u}}\frac{\lvert\mathfrak{u}-\mathfrak{z}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}-\mathfrak{z}\rvert}}\left(\prod_{j\in\mathfrak{u}-\mathfrak{z}}b_{j}\right)K^{*}\prod_{j=1}^{s}\exp\left(b_{j}\left\lvert\Phi^{-1}\left((2w_{j})^{\frac{1}{\beta}}\cdot\frac{1}{2}\right)\right\rvert\right)\right\rvert\\\
&\cdot\prod_{j\in\mathfrak{u}-\mathfrak{z}}\frac{\partial\bar{\mathrm{y}}_{j}}{\partial
t_{j}}\frac{\partial t_{j}}{\partial
w_{j}}\cdot\prod_{j=1}^{s}\frac{1}{\beta_{j}}(2w_{j})^{\frac{1}{\beta_{j}}-1}\cdot\prod_{k\in\mathfrak{z}}\frac{1-\beta_{k}}{\beta_{k}w_{k}}.\end{split}$
(82)

We apply the change of variable
$t_{j}=(2w_{j})^{\frac{1}{\beta_{j}}}\cdot\frac{1}{2}$ for $j=1,\dotsc,s$ and
continue the analysis as follows:

$\displaystyle\begin{split}\left\lvert\frac{\partial}{\partial
t_{\mathfrak{u}}}g_{\mathrm{IS}}(\mathbf{w})\right\rvert&\leq\left\lvert\sum_{\mathfrak{z}\subseteq\mathfrak{u}}\frac{\lvert\mathfrak{u}-\mathfrak{z}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}-\mathfrak{z}\rvert}}\left(\prod_{j\in\mathfrak{u}-\mathfrak{z}}b_{j}\right)K^{*}\prod_{j=1}^{s}\exp\left(b_{j}\left\lvert\Phi^{-1}\left(t_{j}\right)\right\rvert\right)\right\rvert\\\
&\cdot\prod_{j\in\mathfrak{u}-\mathfrak{z}}\frac{\partial\bar{\mathbf{y}}_{j}}{\partial
t_{j}}\frac{\partial t_{j}}{\partial
w_{j}}\cdot\prod_{j=1}^{s}\frac{1}{\beta_{j}}(2w_{j})^{\frac{1}{\beta_{j}}-1}\cdot\prod_{k\in\mathfrak{z}}\frac{1-\beta_{k}}{\beta_{k}w_{k}}\\\
&\approx{\sum_{\mathfrak{z}\subseteq\mathfrak{u}}\frac{\lvert\mathfrak{u}-\mathfrak{z}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}-\mathfrak{z}\rvert}}\left(\prod_{j\in\mathfrak{u}-\mathfrak{z}}b_{j}\right)K^{*}\prod_{j=1}^{s}\exp\left(b_{j}\sqrt{-2\log\left(t_{j}\right)}\right)}\\\
&\cdot\prod_{j\in\mathfrak{u}-\mathfrak{z}}\frac{1}{w_{j}\beta_{j}\sqrt{-2\log\left(t_{j}\right)}}\cdot\prod_{j=1}^{s}\frac{1}{\beta_{j}}(2w_{j})^{\frac{1}{\beta_{j}}-1}\cdot\prod_{k\in\mathfrak{z}}\frac{1-\beta_{k}}{\beta_{k}w_{k}}\\\
&\approx{\sum_{\mathfrak{z}\subseteq\mathfrak{u}}\frac{\lvert\mathfrak{u}-\mathfrak{z}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}-\mathfrak{z}\rvert}}\left(\prod_{j\in\mathfrak{u}-\mathfrak{z}}b_{j}\right)K^{*}\prod_{j=1}^{s}\exp\left(b_{j}\sqrt{-2\log\left(t_{j}\right)}\right)}\\\
&\cdot\prod_{j\in\mathfrak{u}-\mathfrak{z}}\frac{2}{\beta_{j}\sqrt{-2\log\left(t_{j}\right)}}\cdot\prod_{j=1}^{s}\frac{1}{\beta_{j}}(2w_{j})^{\frac{1}{\beta_{j}}-2}\cdot\prod_{k\in\mathfrak{z}}\frac{2(1-\beta_{k})}{\beta_{k}}\\\
&={\sum_{\mathfrak{z}\subseteq\mathfrak{u}}\frac{\lvert\mathfrak{u}-\mathfrak{z}\rvert!}{(\ln
2)^{\lvert\mathfrak{u}-\mathfrak{z}\rvert}}\left(\prod_{j\in\mathfrak{u}-\mathfrak{z}}b_{j}\right)K^{*}\prod_{j=1}^{s}\exp\left(b_{j}\sqrt{-2\log\left(t_{j}\right)}\right)}\\\
&\cdot\prod_{j\in\mathfrak{u}-\mathfrak{z}}\frac{2}{\beta_{j}^{2}\sqrt{-2\log\left(t_{j}\right)}}\cdot\prod_{k\in\mathfrak{z}}\frac{2(1-\beta_{k})}{\beta_{k}^{2}}\cdot\prod_{j=1}^{s}(2t_{j})^{1-2{\beta_{j}}}\\\
&\leq C\prod_{j=1}^{s}\left({t_{j}}\right)^{1-2{\beta_{j}}-B_{j}}\\\
&=C\prod_{j=1}^{s}\left({w_{j}}\right)^{\frac{1}{\beta_{j}}-2-\frac{B_{j}}{\beta_{j}}},\\\
\end{split}$ (83)

where in the last inequality we use,

$\displaystyle\exp\left(b_{j}\sqrt{-2\log\left(t_{j}\right)}\right)\leq
C(t_{j})t_{j}^{-B_{j}}$

with $B_{j}=\frac{b_{j}}{\sqrt{-2\log t_{j}}}$. Compared to the case without
IS ($\bm{\beta}=\mathbf{1}$), an improved regularity of $g_{\mathrm{IS}}$ is
achieved for $0<\beta_{j}<1$, as the exponent becomes larger, i.e.
$\frac{1}{\beta_{j}}-2-\frac{B_{j}}{\beta_{j}}\geq-1-B_{j}$, when $1-B_{j}\geq
0$. Moreover, the condition
$\frac{1}{\beta_{j}}-2-\frac{B_{j}}{\beta_{j}}\geq-1$ for $j=1,\dotsc,s$
ensures the $\mathcal{O}(n^{-1})$ convergence rate, which requires

$\beta_{j}\leq 1-\frac{b_{j}}{\sqrt{-2\log t_{j}}}.$ (84)

Similar to the optimality condition (67) in the first proposal distribution,
$0<\beta_{j}<1$ and sufficiently large and finite $t_{j}^{-1}$ for all
$j=1,\dotsc,s$ and thus sufficiently large and finite $n$ are required for the
$\mathcal{O}(n^{-1})$ convergence rate.

##  5 Numerical results

This section tests the convergence rate of the two examples we have
considered, compares them with the analysis, and demonstrates how IS improves
the rate.

###  5.1 Lognormal random variable expectation

This section approximates the expectation of a lognormal random variable, that
is,

${\mathbb{E}\mspace{-2.0mu}\left[g\right]}={\mathbb{E}\mspace{-2.0mu}\left[\exp\left(\sum_{j=1}^{s}\sigma_{j}\xi_{j}\right)\right]},$

where $\xi_{j}$ are i.i.d. $\mathcal{N}(0,1)$. The expectation is analytically
given by

$\prod_{j=1}^{s}\exp\left(\frac{\sigma_{j}^{2}}{2}\right).$

We test the convergence of the root mean squared error (RMSE) of the QMC
method using the Sobol’ sequence. The RMSE for QMC is computed as follows:

$\displaystyle\sqrt{\mathbb{E}\left[\left(\frac{1}{R}\sum_{r=1}^{R}\frac{1}{n}\sum_{i=1}^{n}g(t_{i}\oplus\Delta_{r})-\mathbb{E}[g]\right)^{2}\right]}\approx\sqrt{\frac{1}{R(R-1)}\sum_{r=1}^{R}\left(\hat{I}_{n}^{(r)}-\hat{I}_{R,n}\right)^{2}}$
(85)

where $\hat{I}_{n}^{(r)}=\frac{1}{n}\sum_{i=1}^{n}g(t_{i}\oplus\Delta_{r})$
and
$\hat{I}_{R,n}=\frac{1}{R}\sum_{r=1}^{R}\frac{1}{n}\sum_{i=1}^{n}g(t_{i}\oplus\Delta_{r})$.
In Example 1, we choose the nested uniform scrambling randomization (Owen’s
scrambling). The convergence rate for the RMSE is
$\mathcal{O}(n^{-3/2+\epsilon})$ if the integrand has Lipschitz-continuous
mixed first-order derivatives [36]. The integrand $f$ has a singularity at 0
and thus does not have this property.

Figures 5, 6, and 7 plot the convergence of the RMSE for the original
integrand $g$ and the integrand with importance sampling, $g_{\textrm{IS}}$.
IS improves the integrand regularity in all cases, significantly reducing the
RMSE and improving the convergence rate.

Tables 1, 2, and 3 provide the measured convergence rates $\gamma$ and the
value $(1-\gamma)/\sqrt{\sum_{j=1}^{s}\sigma_{j}^{2}}$. For the same sample
size $n$, the convergence rate $\gamma$ decreases as
$\sqrt{\sum_{j=1}^{s}\sigma_{j}^{2}}$ increases. The value
$(1-\gamma)/\sqrt{\sum_{j=1}^{s}\sigma_{j}^{2}}$ exhibits similarity in
magnitude for all cases in each dimension setting and increases as the
dimension increases. Figure 8 plots the values
$\frac{2(1-\gamma)^{2}}{\sum_{j=1}^{s}\sigma_{j}^{2}}$ (left) and $Cn^{-1}$
with the corresponding 95% confidence intervals. The variation within each
group is also observed to be smaller than the values across groups.

Table 1: Convergence of Example 1 for dimension $s=1$ and standard deviation $\sigma=1.0,2.0,3.0$. The convergence rate is fitted at $n=2^{21}$.  $\sigma$ | Conv. Rate $\gamma$ | Conv. Rate with IS | $(1-\gamma)/\sigma$
---|---|---|---
1.0 | 0.875 | 1.888 | 0.125
2.0 | 0.719 | 1.340 | 0.141
3.0 | 0.582 | 1.088 | 0.139

Figure 5: Example 1, the root mean squared error with $s=1$: $\sigma=1.0$
(upper left), 2.0 (upper right), and 3.0 (bottom). The convergence rate is
fitted at $n=2^{21}$. The legend IS_A stands for the IS of the first type.

Figure 6: Example 1, the root mean squared error with $s=3$:
$\sigma=\\{1.0,1.0,1.0\\}$ (upper left), $\\{2.0,1.0,1.0\\}$ (upper right),
$\\{2.0,1.4,1.0\\}$ (lower left), and $\\{1.0,1.7,1.0\\}$ (lower right). The
convergence rate is fitted centered at $n=2^{21}$. The legend IS_A stands for
the IS of the first type.

Figure 7: Example 1, the root mean squared error with $s=6$: Case I (upper left), Case II (upper right), Case III (lower left), and Case IV (lower right). The convergence rate is fitted at $n=2^{21}$. The legend IS_A stands for the IS of the first type.  Table 2: Convergence of Example 1, for dimension $s=3$ and various standard deviation settings, where the convergence rate is fitted at $n=2^{21}$.  $\sigma$ | Conv. Rate $\gamma$ | Conv. Rate with IS | $(1-\gamma)/\sqrt{\sum_{j=1}^{s}\sigma_{j}^{2}}$
---|---|---|---
{1.0, 1.0, 1.0} | 0.717 | 1.264 | 0.163
{2.0, 1.0, 1.0} | 0.614 | 1.305 | 0.158
{2.0, 1.4, 1.0} | 0.575 | 1.349 | 0.161
{2.0, 1.7, 1.0} | 0.541 | 1.354 | 0.162
Table 3: Convergence of Example 1 for dimension $s=6$ and various standard deviation $\sigma$ settings, where ${\sigma}_{j}={1}$, $j=1,\dotsc,s$ for Case I, and $\sigma_{j}=\frac{\sqrt{6}\xi_{j}}{\sqrt{\sum_{j=1}^{s}\xi_{j}^{2}}}$, where $\xi_{j}\sim\textrm{Lognormal}(0,1)$, $j=1,\dotsc,s$ for Cases II, III, and IV. The convergence rate is fitted at $n=2^{21}$.  $\sigma$ cases | Conv. Rate $\gamma$ | Conv. Rate with IS | $(1-\gamma)/\sqrt{\sum_{j=1}^{s}\sigma_{j}^{2}}$
---|---|---|---
I | 0.577 | 0.962 | 0.173
II | 0.571 | 1.000 | 0.175
III | 0.575 | 1.014 | 0.174
IV | 0.567 | 1.031 | 0.178

Figure 8: Example 1, $\frac{2(1-\gamma)^{2}}{\sum_{j=1}^{s}\sigma_{j}^{2}}$
(left) and $Cn^{-1}$ (right) against the number of samples $n$ across various
dimension and variance settings. The colors represent the dimension $s$,
whereas each scatter point with the same color represents a unique variance
setting. The shaded area is the 95% confidence interval. Both the values
$\frac{2(1-\gamma)^{2}}{\sum_{j=1}^{s}\sigma_{j}^{2}}$(left) and $Cn^{-1}$
decrease as the sample size increases. The values exhibit similarity in
amplitude within each dimension setting.

###  5.2 Elliptic partial differential equations with lognormal coefficient

We first specify the settings of this example. The QoI is the weighted
integration of the solution $u$ over the entire domain $\mathcal{D}$, that is,

$\displaystyle
Q(u)=\int_{\mathcal{D}}g(\mathbf{x})u(\mathbf{x};\cdot)d\mathbf{x},$

where $\mathcal{D}=[-1,1]^{2}$, $g=\rho(;\mu,\Sigma)*1_{D_{0}}$, $*$ denotes
the convolution operator,
$\mathcal{D}_{0}=\left[0.25,0.5\right]\times\left[-0.5,-0.25\right]\subset\mathcal{D}$,
and $\rho(;\mu,\Sigma)$ represents the Gaussian density function in
$\mathbb{R}^{2}$ with mean $\mu=0$ and covariance $\Sigma=\frac{1}{16}\bm{I}$.
The PDE problem (45a) is solved with DEAL.II [1], using the $\mathcal{Q}_{1}$
finite element on a 16$\times$16 mesh. In this example, we employed the Sobol’
sequence with a random linear scramble as the randomization [22].

Recall that the coefficient $a$ involved in the PDE model (45) is given by,

$\displaystyle
a(\mathbf{x};\omega)=\exp\left(\sum_{j=1}^{s}y_{j}\psi_{j}\right).$

where the basis $\psi_{j}$ is connected to the following Matérn covariance
kernel:

$\displaystyle
C(h)=\frac{1}{2^{\nu-1}\Gamma({\nu})}\left(\sqrt{2\nu}\frac{h}{r}\right)^{\nu}K_{\nu}\left(\sqrt{2\nu}\frac{h}{r}\right),$
(86)

where we chose $\nu=4.5$ and $r=1$. In addition, $\Gamma$ is the Gamma
function, and $K_{\nu}$ is the modified Bessel function of the second kind.
Let $\lambda_{j,p}$ and $\theta_{j,p}$ be the Fourier coefficients and
trigonometric functions in the Fourier expansion of the covariance kernel $C$
on an extended domain $\mathcal{D}_{p}=[-\gamma,\gamma]^{d}$, respectively,
with $d=2$. The extension is needed to ensure the positivity of the Fourier
coefficients $\lambda_{j,p}$ [3]. We let

$\displaystyle\begin{split}\psi_{j,p}&=\lambda_{j,p}\theta_{j,p}\\\
\psi_{j}&=\left.\psi_{j,p}\right|_{\mathcal{D}}.\end{split}$ (87)

Apart from the trigonometric-basis expansions, we also apply Meyer wavelet
functions to construct the basis $\psi_{j}$ (for a detailed analysis and
instructions, see Section 4 in [3]).

To validate the proposed convergence model (55), we also considered the
following coefficient:

$\displaystyle
a(\mathbf{x};\omega)=\exp\left(\sum_{j=1}^{s}y_{j}\sigma_{j}\psi_{j}\right),$
(88)

where we generated six samples of $\sigma$ within each dimension setting. In
the first sample, $\sigma=\mathbf{1}$. In the second and third samples,
$\sigma_{j}$ are i.i.d sampled from the uniform distribution $U[1,2]$. In the
last three samples, the $\sigma$ values of the first three cases are
multiplied component-wise by 2. Figure 9 plots the deteriorated rate
$1-\gamma$ against $\sqrt{\sum_{j=1}^{s}\sigma_{j}^{2}b_{j}^{2}}$. In all
considered cases, the rate $1-\gamma$ exhibits an almost linear dependence on
$\sqrt{\sum_{j=1}^{s}\sigma_{j}^{2}b_{j}^{2}}$. Unlike Example 1, the effect
of dimensions on the convergence rates is not obvious in both cases, which was
unexpected because the QMC points quality deteriorates in higher dimensions.
One possible reason for this observation lies in Equation. (49), where we
applied the upper bounds for $Q$ and its derivatives.

Figure 9: Example 2, the value $(1-\gamma)$ against
$\sqrt{\sum_{j=1}^{s}\sigma_{j}^{2}b_{j}^{2}}$ for using the Fourier basis
(left) and wavelet-type basis (right) across various dimension and variance
settings. The colors represent the dimension $s$, whereas each scatter point
with the same color represents a unique variance setting. The dimension-
independent effect is observed in the convergence rate in both cases.

Figure 10 compares the $L^{\infty}$ norm of the basis ${\psi_{j,p}}$ and
${\psi_{j}}$ between trigonometric and wavelet-type bases. The trigonometric
basis is optimal in approximating the random field in the $L^{2}$ sense.
However, the basis $L^{\infty}$ norm, which appears in the convergence model
(55) and (29) is the primary interest. Ideally,
$\lVert\psi_{j}\rVert_{L^{\infty}(\mathcal{D})}$ is desired to converge fast
as $j$ increases to reduce the effective dimension. In the nonasymptotic QMC
convergence model for elliptic PDEs (55), the value $\sum_{j=1}^{s}b_{j}^{2}$
is expected to be small for a fixed dimension $s$. We observed that the values
$b_{j}$ of the wavelet-type basis in the plotted range are smaller than those
of trigonometric basis. The localization properties of the wavelet-type basis
make it even more favorable when the function $\psi_{j,p}$ is restricted on
$\mathcal{D}$ since
$\lVert\psi_{j}\rVert_{L^{\infty}(\mathcal{D})}\leq\lVert\psi_{j,p}\rVert_{L^{\infty}(\mathcal{D}_{p})}$.

Figure 10: Example 2, a comparison of the $L^{2}$ norms
$\lVert\psi_{j,p}\rVert_{L^{2}(\mathcal{D}_{p})}$ and $L^{\infty}$ norms
$\lVert\psi_{j}\rVert_{L^{\infty}(\mathcal{D})}$ between the Fourier basis and
wavelet-type basis.

Figure 11 plots the convergence of the RMSE for the original integrand $g$ and
the integrand with IS $g_{\textrm{IS}}$, where the wavelet-type basis is used
for two variance settings and two dimension settings. In all cases, both types
of IS improve the integrand regularity, reducing the RMSE and improving the
convergence rate. Furthermore, the effectiveness of the IS is more significant
when the dimension is smaller, as IS has more limitations in higher dimensions
to avoid increasing the variance. Moreover, when the original integrand has a
larger $\sigma^{2}$ the IS works better because $\alpha_{j}$ in (66) becomes
larger and $\beta_{j}$ in (83) becomes smaller.

Figure 11: Example 2, the root mean squared error with and without two types
of IS for $s=16$ (top), $s=64$ (bottom) with $\sigma^{2}=1$ (left) and
$\sigma^{2}=4$ (right). The fitted convergence rate is marked in the figure.
The random coefficient $a$ is expressed using the wavelet-type basis. Legends
IS_A and IS_B stand for IS of the first and the second type, respectively.

##  6 Conclusion

In this work, we have studied the nonasymptotic QMC convergence rate model for
functions with finite-dimensional inputs. Specifically, we focus on the
expectation of a lognormal random variable and an elliptic PDE problem
characterized by lognormal coefficients. Drawing upon the hyperbolic set
$K_{n,s}$ introduced in the work of Owen [37], we subdivided the integration
domain. This division allowed us to split the QMC integration error into two
distinct contributions, namely: $2\int_{[0,1]^{s}-K_{n,s}}\lvert
g-\tilde{g}\rvert$ and
$\mathbb{E}\left[D_{n}^{*}(t_{1},\dotsc,t_{n})\right]V_{HK}(\tilde{g})$. This
method replaces the right-hand side of the Koksma–Hlawka inequality
$V_{\textrm{HK}}(g)D^{*}_{n}(\mathcal{P})$, which is infinite for the
unbounded functions considered in this work, with two finite terms.

Our primary contribution through this work has been deriving an upper bound
for the integrand derivative via an optimization problem, denoted as (23).
With this in place, we applied it within the convergence model, culminating in
deriving the nonasymptotic QMC convergence rate model. In our quest to verify
the proposed model, we presented numerical examples for the two problems
above. Moreover, the analytical procedures delineated here can find relevance
in varied integration challenges, one notable instance being the option
pricing problem within financial scenarios.

In addition to the above, our work also involved applying two IS
distributions, leading us to evaluate their potential to enhance integrand
regularity. Notably, despite the distinct support of these distributions,
their influence on the integrand appeared broadly analogous.

However, there remain areas that we have not explored in depth. For instance,
we have yet to derive the constant $C$ within the hyperbolic set. Delving into
its dependence on dimensionality and integrand characteristics could be a
promising trajectory for subsequent studies. Moreover, future inquiries could
also encompass the random field’s truncation error, where wavelets demonstrate
a superior edge over the trigonometric basis, as indicated in [2].
Furthermore, examining a multilevel setting with adaptivity, as hinted at in
[8], could be beneficial.

Lastly, it’s worth noting our choice to employ the Sobol’ sequence instead of
exploring weighted spaces or designing lattice rules (e.g., [29], [16]). The
latter could reduce the worst-case integration error for specific function
spaces, potentially leading to a dimension-independent outcome given certain
conditions. This combination of lattice rules and nonasymptotic analysis
remains a compelling avenue for future exploration.

##  7 Acknowledgement

This publication is based on work supported by the Alexander von Humboldt
Foundation and the King Abdullah University of Science and Technology (KAUST)
office of sponsored research (OSR) under Award No. OSR-2019-CRG8-4033. This
work utilized the resources of the Supercomputing Laboratory at King Abdullah
University of Science and Technology (KAUST) in Thuwal, Saudi Arabia. The
authors thank Christian Bayer, Fabio Nobile and Erik von Schwerin for fruitful
discussions. We are also grateful to Michael Samet for the suggestions that
significantly improve the paper’s readibility. We acknowledge the use of the
following open-source software packages: deal.II [1].

## References

  * [1] D. Arndt, W. Bangerth, B. Blais, T. C. Clevenger, M. Fehling, A. V. Grayver, T. Heister, L. Heltai, M. Kronbichler, M. Maier, P. Munch, J.-P. Pelteret, R. Rastak, I. Thomas, B. Turcksin, Z. Wang, and D. Wells, The deal.II library, version 9.2, Journal of Numerical Mathematics, 28 (2020), pp. 131–146.
  * [2] M. Bachmayr, A. Cohen, R. DeVore, and G. Migliorati, Sparse polynomial approximation of parametric elliptic PDEs. Part II: lognormal coefficients, ESAIM: Mathematical Modelling and Numerical Analysis, 51 (2017), pp. 341–363.
  * [3] M. Bachmayr, A. Cohen, and G. Migliorati, Representations of Gaussian random fields and approximation of elliptic PDEs with lognormal coefficients, Journal of Fourier Analysis and Applications, 24 (2018), pp. 621–649.
  * [4] C. Bayer, C. Ben Hammouda, and R. Tempone, Numerical smoothing with hierarchical adaptive sparse grids and quasi-Monte Carlo methods for efficient option pricing, Quantitative Finance, 23 (2023), pp. 209–227.
  * [5] C. Bayer, C. B. Hammouda, A. Papapantoleon, M. Samet, and R. Tempone, Optimal damping with hierarchical adaptive quadrature for efficient Fourier pricing of multi-asset options in Lévy models, ArXiv, abs/2203.08196 (2022).
  * [6] C. Bayer, C. B. Hammouda, and R. Tempone, Numerical smoothing and hierarchical approximations for efficient option pricing and density estimation, arXiv: Computational Finance, (2020).
  * [7] C. Bayer, M. Siebenmorgen, and R. Tempone, Smoothing the payoff for efficient computation of basket option prices, Quantitative Finance, 18 (2016), pp. 491 – 505.
  * [8] J. Beck, Y. Liu, E. von Schwerin, and R. Tempone, Goal-oriented adaptive finite element multilevel Monte Carlo with convergence rates, Computer Methods in Applied Mechanics and Engineering, 402 (2022), p. 115582.
  * [9] R. E. Caflisch, W. J. Morokoff, and A. B. Owen, Valuation of mortgage backed securities using Brownian bridges to reduce effective dimension, vol. 24, Department of Mathematics, University of California, Los Angeles, 1997.
  * [10] J. Dick and F. Pillichshammer, Digital nets and sequences: discrepancy theory and quasi-Monte Carlo integration, Cambridge University Press, 2010.
  * [11] G. Y. Dong and C. Lemieux, Dependence properties of scrambled Halton sequences, Mathematics and Computers in Simulation, 200 (2022), pp. 240–262.
  * [12] M. Gerber and N. Chopin, Sequential quasi Monte Carlo, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77 (2015), pp. 509–579.
  * [13] A. D. Gilbert, F. Y. Kuo, and I. H. Sloan, Preintegration is not smoothing when monotonicity fails, in Advances in Modeling and Simulation: Festschrift for Pierre L’Ecuyer, Springer, 2022, pp. 169–191.
  * [14] E. Gobet, M. Lerasle, and D. Métivier, Mean estimation for randomized quasi Monte Carlo method, (2022).
  * [15] T. Goda and P. L’ecuyer, Construction-free median quasi-Monte Carlo rules for function spaces with unspecified smoothness and general weights, SIAM Journal on Scientific Computing, 44 (2022), pp. A2765–A2788.
  * [16] I. G. Graham, F. Y. Kuo, J. A. Nichols, R. Scheichl, C. Schwab, and I. H. Sloan, Quasi-Monte Carlo finite element methods for elliptic PDEs with lognormal random coefficients, Numerische Mathematik, 131 (2015), pp. 329–368.
  * [17] M. Griebel, F. Y. Kuo, and I. H. Sloan, The smoothing effect of integration in $\mathbb{R}^{d}$ and the ANOVA decomposition, Math. Comput., 82 (2013), pp. 383–400.
  * [18] A. Griewank, F. Y. Kuo, H. Leövey, and I. H. Sloan, High dimensional integration of kinks and jumps—smoothing by preintegration, Journal of Computational and Applied Mathematics, 344 (2018), pp. 259–274.
  * [19] Z. He and X. Wang, On the convergence rate of randomized quasi-Monte Carlo for discontinuous functions, SIAM Journal on Numerical Analysis, 53 (2015), pp. 2488–2503.
  * [20] Z. He, Z. Zheng, and X. Wang, On the error rate of importance sampling with randomized quasi-Monte Carlo, arXiv preprint arXiv:2203.03220, (2022).
  * [21] L. Herrmann and C. Schwab, QMC integration for lognormal-parametric, elliptic PDEs: local supports and product weights, Numerische Mathematik, 141 (2019), pp. 63–102.
  * [22] H. S. Hong and F. J. Hickernell, Algorithm 823: Implementing scrambled digital sequences, ACM Transactions on Mathematical Software (TOMS), 29 (2003), pp. 95–109.
  * [23] J. Imai and K. S. Tan, A general dimension reduction technique for derivative pricing, Journal of Computational Finance, 10 (2006), p. 129.
  * [24] , Pricing derivative securities using integrated quasi-Monte Carlo methods with dimension reduction and discontinuity realignment, SIAM Journal on Scientific Computing, 36 (2014), pp. A2101–A2121.
  * [25] X. Jin and A. X. Zhang, Reclaiming Quasi-Monte Carlo efficiency in portfolio value-at-risk simulation through Fourier transform, Management Science, 52 (2006), pp. 925–938.
  * [26] Y. Kazashi, Quasi-Monte Carlo integration with product weights for elliptic pdes with log-normal coefficients, IMA Journal of Numerical Analysis, 39 (2019), pp. 1563–1593.
  * [27] P. Kritzer, F. Pillichshammer, L. Plaskota, and G. W. Wasilkowski, On efficient weighted integration via a change of variables, Numerische Mathematik, 146 (2018), pp. 545 – 570.
  * [28] F. Y. Kuo and D. Nuyens, A practical guide to quasi-Monte Carlo methods, (2016).
  * [29] F. Y. Kuo, C. Schwab, and I. H. Sloan, Quasi-Monte Carlo finite element methods for a class of elliptic partial differential equations with random coefficients, SIAM J. Numer. Anal., 50 (2012), pp. 3351–3374.
  * [30] S. Liu and A. B. Owen, Preintegration via active subspace, SIAM Journal on Numerical Analysis, 61 (2023), pp. 495–514.
  * [31] P. L’Ecuyer, Randomized quasi-Monte Carlo: An introduction for practitioners, Springer, 2018.
  * [32] W. J. Morokoff and R. E. Caflisch, Quasi-random sequences and their discrepancies, SIAM Journal on Scientific Computing, 15 (1994), pp. 1251–1279.
  * [33] J. A. Nichols and F. Y. Kuo, Fast CBC construction of randomly shifted lattice rules achieving $\mathcal{O}(n^{-1+\delta})$ convergence for unbounded integrands over $\mathbb{R}^{s}$ in weighted spaces with POD weights, Journal of Complexity, 30 (2014), pp. 444–468.
  * [34] H. Niederreiter, Random number generation and quasi-Monte Carlo methods, SIAM, 1992.
  * [35] D. Ouyang, X. Wang, and Z. He, Quasi-Monte Carlo for unbounded integrands with importance sampling, arXiv preprint arXiv:2310.00650, (2023).
  * [36] A. B. Owen, Scrambling Sobol’ and Niederreiter–Xing points, Journal of complexity, 14 (1998), pp. 466–489.
  * [37] A. B. Owen, Halton sequences avoid the origin, SIAM Review, 48 (2006), pp. 487–503.
  * [38] A. B. Owen, Monte Carlo theory, methods and examples, 2013.
  * [39] A. B. Owen, A randomized Halton algorithm in R, arXiv preprint arXiv:1706.02808, (2017).
  * [40] Z. Pan and A. Owen, Super-polynomial accuracy of one dimensional randomized nets using the median of means, Mathematics of Computation, 92 (2023), pp. 805–837.
  * [41] A. Papageorgiou, The Brownian bridge does not offer a consistent advantage in quasi-Monte Carlo integration, Journal of Complexity, 18 (2002), pp. 171–186.
  * [42] J. K. Patel and C. B. Read, Handbook of the normal distribution, vol. 150, CRC Press, 1996.
  * [43] I. Pinelis, Exact lower and upper bounds on the incomplete gamma function, arXiv preprint arXiv:2005.06384, (2020).
  * [44] I. H. Sloan and H. Woźniakowski, When are quasi-Monte Carlo algorithms efficient for high dimensional integrals?, Journal of Complexity, 14 (1998), pp. 1–33.
  * [45] I. Sobol, Calculation of improper integrals using equidistributed sequences, Doklady Akademii Nauk SSSR, 210 (1973), pp. 278–281.
  * [46] X. Wang, Dimension reduction techniques in quasi-Monte Carlo methods for option pricing, INFORMS Journal on Computing, 21 (2009), pp. 488–504.
  * [47] X. Wang and K.-T. Fang, The effective dimension and quasi-Monte Carlo integration, Journal of Complexity, 19 (2003), pp. 101–124.
  * [48] X. Wang and I. H. Sloan, Why are high-dimensional finance problems often of low effective dimension?, SIAM Journal on Scientific Computing, 27 (2005), pp. 159–183.
  * [49] , Low discrepancy sequences in high dimensions: How well are their projections distributed?, Journal of Computational and Applied Mathematics, 213 (2008), pp. 366–386.
  * [50] , Quasi-Monte Carlo methods in financial engineering: An equivalence principle and dimension reduction, Operations Research, 59 (2011), pp. 80–95.
  * [51] X. Wang and K. S. Tan, Pricing and hedging with discontinuous functions: Quasi-Monte Carlo methods and dimension reduction, Management Science, 59 (2013), pp. 376–389.
  * [52] C. Weng, X. Wang, and Z. He, Efficient computation of option prices and greeks by quasi-Monte Carlo method with smoothing and dimension reduction, SIAM Journal on Scientific Computing, 39 (2017), pp. B298–B322.
  * [53] J. Wiart, C. Lemieux, and G. Y. Dong, On the dependence structure and quality of scrambled $(t,m,s)$-nets, Monte Carlo Methods and Applications, 27 (2019), pp. 1 – 26.