Search

We analyse and implement a quasi-Monte Carlo (QMC) finite element method (FEM) for the forward problem of uncertainty quantification (UQ) for the Helmholtz equation with random coefficients, both in the second-order and zero-order terms of the equation, thus modelling wave scattering in random media. The problem is formulated on the infinite propagation domain, after scattering by the heterogeneity, and also (possibly) a bounded impenetrable scatterer. The spatial discretization scheme includes truncation to a bounded domain via a perfectly matched layer (PML) technique and then FEM approximation. A special case is the problem of an incident plane wave being scattered by a bounded sound-soft impenetrable obstacle surrounded by a random heterogeneous medium, or more simply, just scattering by the random medium. The random coefficients are assumed to be affine separable expansions with infinitely many independent uniformly distributed and bounded random parameters. As quantities of interest for the UQ, we consider the expectation of general linear functionals of the solution, with a special case being the far-field pattern of the scattered field. The numerical method consists of (a) dimension truncation in parameter space, (b) application of an adapted QMC method to compute expected values, and (c) computation of samples of the PDE solution via PML truncation and FEM approximation. Our error estimates are explicit in $s$ (the dimension truncation parameter), $N$ (the number of QMC points), $h$ (the FEM grid size) and (most importantly), $k$ (the Helmholtz wavenumber). The method is also exponentially accurate with respect to the PML truncation radius. Illustrative numerical experiments are given.

In this paper we consider Deep Neural Networks (DNNs) with a smooth activation function as surrogates for high-dimensional functions that are somewhat smooth but costly to evaluate. We consider the standard (non-periodic) DNNs as well as propose a new model of periodic DNNs which are especially suited for a class of periodic target functions when Quasi-Monte Carlo lattice points are used as training points. We study the regularity of DNNs by obtaining explicit bounds on all mixed derivatives with respect to the input parameters. The bounds depend on the neural network parameters as well as the choice of activation function. By imposing restrictions on the network parameters to match the regularity features of the target functions, we prove that DNNs with $N$ tailor-constructed lattice training points can achieve the generalization error (or $L_2$ approximation error) bound ${\tt tol} + \mathcal{O}(N^{-r/2})$, where ${\tt tol}\in (0,1)$ is the tolerance achieved by the training error in practice, and $r = 1/p^*$, with $p^*$ being the ``summability exponent'' of a sequence that characterises the decay of the input variables in the target functions, and with the implied constant independent of the dimensionality of the input data. We apply our analysis to popular models of parametric elliptic PDEs in uncertainty quantification. In our numerical experiments, we restrict the network parameters during training by adding tailored regularization terms, and we show that for an algebraic equation mimicking the parametric PDE problems the DNNs trained with tailored regularization perform significantly better.

In this paper, we apply quasi-Monte Carlo (QMC) methods with an initial preintegration step to estimate cumulative distribution functions and probability density functions in uncertainty quantification (UQ). The distribution and density functions correspond to a quantity of interest involving the solution to an elliptic partial differential equation (PDE) with a lognormally distributed coefficient and a normally distributed source term. There is extensive previous work on using QMC to compute expected values in UQ, which have proven very successful in tackling a range of different PDE problems. However, the use of QMC for density estimation applied to UQ problems will be explored here for the first time. Density estimation presents a more difficult challenge compared to computing the expected value due to discontinuities present in the integral formulations of both the distribution and density. Our strategy is to use preintegration to eliminate the discontinuity by integrating out a carefully selected random parameter, so that QMC can be used to approximate the remaining integral. First, we establish regularity results for the PDE quantity of interest that are required for smoothing by preintegration to be effective. We then show that an $N$-point lattice rule can be constructed for the integrands corresponding to the distribution and density, such that after preintegration the QMC error is of order $\mathcal{O}(N^{-1+\epsilon})$ for arbitrarily small $\epsilon>0$. This is the same rate achieved for computing the expected value of the quantity of interest. Numerical results are presented to reaffirm our theory.

We show that a very simple randomised algorithm for numerical integration can produce a near optimal rate of convergence for integrals of functions in the $d$-dimensional weighted Korobov space. This algorithm uses a lattice rule with a fixed generating vector and the only random element is the choice of the number of function evaluations. For a given computational budget $n$ of a maximum allowed number of function evaluations, we uniformly pick a prime $p$ in the range $n/2 < p \le n$. We show error bounds for the randomised error, which is defined as the worst case expected error, of the form $O(n^{-\alpha - 1/2 + \delta})$, with $\delta > 0$, for a Korobov space with smoothness $\alpha > 1/2$ and general weights. The implied constant in the bound is dimension-independent given the usual conditions on the weights. We present an algorithm that can construct suitable generating vectors \emph{offline} ahead of time at cost $O(d n^4 / \ln n)$ when the weight parameters defining the Korobov spaces are so-called product weights. For this case, numerical experiments confirm our theory that the new randomised algorithm achieves the near optimal rate of the randomised error., Comment: 29 pages, 2 figures

The kernel interpolant in a reproducing kernel Hilbert space is optimal in the worst-case sense among all approximations of a function using the same set of function values. In this paper, we compare two search criteria to construct lattice point sets for use in lattice-based kernel approximation. The first candidate, $\calP_n^*$, is based on the power function that appears in machine learning literature. The second, $\calS_n^*$, is a search criterion used for generating lattices for approximation using truncated Fourier series. We find that the empirical difference in error between the lattices constructed using $\calP_n^*$ and $\calS_n^*$ is marginal. The criterion $\calS_n^*$ is preferred as it is computationally more efficient and has a proven error bound., Comment: 16 pages, 3 figures

We describe a fast method for solving elliptic partial differential equations (PDEs) with uncertain coefficients using kernel interpolation at a lattice point set. By representing the input random field of the system using the model proposed by Kaarnioja, Kuo, and Sloan (SIAM J.~Numer.~Anal.~2020), in which a countable number of independent random variables enter the random field as periodic functions, it was shown by Kaarnioja, Kazashi, Kuo, Nobile, and Sloan (Numer.~Math.~2022) that the lattice-based kernel interpolant can be constructed for the PDE solution as a function of the stochastic variables in a highly efficient manner using fast Fourier transform (FFT). In this work, we discuss the connection between our model and the popular ``affine and uniform model'' studied widely in the literature of uncertainty quantification for PDEs with uncertain coefficients. We also propose a new class of weights entering the construction of the kernel interpolant -- \emph{serendipitous weights} -- which dramatically improve the computational performance of the kernel interpolant for PDE problems with uncertain coefficients, and allow us to tackle function approximation problems up to very high dimensionalities. Numerical experiments are presented to showcase the performance of the serendipitous weights.

In this paper we propose and analyse a method for estimating three quantities related to an Asian option: the fair price, the cumulative distribution function, and the probability density. The method involves preintegration with respect to one well chosen integration variable to obtain a smooth function of the remaining variables, followed by the application of a tailored lattice Quasi-Monte Carlo rule to integrate over the remaining variables.

We consider uncertainty quantification for the Poisson problem subject to domain uncertainty. For the stochastic parameterization of the random domain, we use the model recently introduced by Kaarnioja, Kuo, and Sloan (SIAM J. Numer. Anal., 2020) in which a countably infinite number of independent random variables enter the random field as periodic functions. We develop lattice quasi-Monte Carlo (QMC) cubature rules for computing the expected value of the solution to the Poisson problem subject to domain uncertainty. These QMC rules can be shown to exhibit higher order cubature convergence rates permitted by the periodic setting independently of the stochastic dimension of the problem. In addition, we present a complete error analysis for the problem by taking into account the approximation errors incurred by truncating the input random field to a finite number of terms and discretizing the spatial domain using finite elements. The paper concludes with numerical experiments demonstrating the theoretical error estimates., Comment: 38 pages, 3 figures

We approximate $d$-variate periodic functions in weighted Korobov spaces with general weight parameters using $n$ function values at lattice points. We do not limit $n$ to be a prime number, as in currently available literature, but allow any number of points, including powers of $2$, thus providing the fundamental theory for construction of embedded lattice sequences. Our results are constructive in that we provide a component-by-component algorithm which constructs a suitable generating vector for a given number of points or even a range of numbers of points. It does so without needing to construct the index set on which the functions will be represented. The resulting generating vector can then be used to approximate functions in the underlying weighted Korobov space. We analyse the approximation error in the worst-case setting under both the $L_2$ and $L_{\infty}$ norms. Our component-by-component construction under the $L_2$ norm achieves the best possible rate of convergence for lattice-based algorithms, and the theory can be applied to lattice-based kernel methods and splines. Depending on the value of the smoothness parameter $\alpha$, we propose two variants of the search criterion in the construction under the $L_{\infty}$ norm, extending previous results which hold only for product-type weight parameters and prime $n$. We also provide a theoretical upper bound showing that embedded lattice sequences are essentially as good as lattice rules with a fixed value of $n$. Under some standard assumptions on the weight parameters, the worst-case error bound is independent of $d$., Comment: 27 pages, 4 figures

We study the application of a tailored quasi-Monte Carlo (QMC) method to a class of optimal control problems subject to parabolic partial differential equation (PDE) constraints under uncertainty: the state in our setting is the solution of a parabolic PDE with a random thermal diffusion coefficient, steered by a control function. To account for the presence of uncertainty in the optimal control problem, the objective function is composed with a risk measure. We focus on two risk measures, both involving high-dimensional integrals over the stochastic variables: the expected value and the (nonlinear) entropic risk measure. The high-dimensional integrals are computed numerically using specially designed QMC methods and, under moderate assumptions on the input random field, the error rate is shown to be essentially linear, independently of the stochastic dimension of the problem -- and thereby superior to ordinary Monte Carlo methods. Numerical results demonstrate the effectiveness of our method.

Preintegration is a technique for high-dimensional integration over $d$-dimensional Euclidean space, which is designed to reduce an integral whose integrand contains kinks or jumps to a $(d-1)$-dimensional integral of a smooth function. The resulting smoothness allows efficient evaluation of the $(d-1)$-dimensional integral by a Quasi-Monte Carlo or Sparse Grid method. The technique is similar to conditional sampling in statistical contexts, but the intention is different: in conditional sampling the aim is to reduce the variance, rather than to achieve smoothness. Preintegration involves an initial integration with respect to one well chosen real-valued variable. Griebel, Kuo, Sloan [Math. Comp. 82 (2013), 383--400] and Griewank, Kuo, Le\"ovey, Sloan [J. Comput. Appl. Maths. 344 (2018), 259--274] showed that the resulting $(d-1)$-dimensional integrand is indeed smooth under appropriate conditions, including a key assumption -- the integrand of the smooth function underlying the kink or jump is strictly monotone with respect to the chosen special variable when all other variables are held fixed. The question addressed in this paper is whether this monotonicity property with respect to one well chosen variable is necessary. We show here that the answer is essentially yes, in the sense that without this property the resulting $(d-1)$-dimensional integrand is generally not smooth, having square-root or other singularities.

In this paper, we analyse a method for approximating the distribution function and density of a random variable that depends in a non-trivial way on a possibly high number of independent random variables, each with support on the whole real line. Starting with the integral formulations of the distribution and density, the method involves smoothing the original integrand by preintegration with respect to one suitably chosen variable, and then applying a suitable quasi-Monte Carlo (QMC) method to compute the integral of the resulting smoother function. Interpolation is then used to reconstruct the distribution or density on an interval. The preintegration technique is a special case of conditional sampling, a method that has previously been applied to a wide range of problems in statistics and computational finance. In particular, the pointwise approximation studied in this work is a specific case of the conditional density estimator previously considered in L'Ecuyer et al., arXiv:1906.04607. Our theory provides a rigorous regularity analysis of the preintegrated function, which is then used to show that the errors of the pointwise and interpolated estimators can both achieve nearly first-order convergence. Numerical results support the theory.

We investigate the application of efficient recursive numerical integration strategies to models in lattice gauge theory from quantum field theory. Given the coupling structure of the physics problems and the group structure within lattice cubature rules for numerical integration, we show how to approach these problems efficiently by means of Fast Fourier Transform techniques. In particular, we consider applications to the quantum mechanical rotor and compact U(1) lattice gauge theory, where the physical dimensions are two and three. This proceedings article reviews our results presented in J. Comput. Phys 443 (2021) 110527., Comment: arXiv admin note: text overlap with arXiv:2011.05451

We prove that a variant of the classical Sobolev space of first-order dominating mixed smoothness is equivalent (under a certain condition) to the unanchored ANOVA space on $\mathbb{R}^d$, for $d \geq 1$. Both spaces are Hilbert spaces involving weight functions, which determine the behaviour as different variables tend to $\pm \infty$, and weight parameters, which represent the influence of different subsets of variables. The unanchored ANOVA space on $\mathbb{R}^d$ was initially introduced by Nichols & Kuo in 2014 to analyse the error of quasi-Monte Carlo (QMC) approximations for integrals on unbounded domains; whereas the classical Sobolev space of dominating mixed smoothness was used as the setting in a series of papers by Griebel, Kuo & Sloan on the smoothing effect of integration, in an effort to develop a rigorous theory on why QMC methods work so well for certain non-smooth integrands with kinks or jumps coming from option pricing problems. In this same setting, Griewank, Kuo, Le\"ovey & Sloan in 2018 subsequently extended these ideas by developing a practical smoothing by preintegration technique to approximate integrals of such functions with kinks or jumps. We first prove the equivalence in one dimension (itself a non-trivial task), before following a similar, but more complicated, strategy to prove the equivalence for general dimensions. As a consequence of this equivalence, we analyse applying QMC combined with a preintegration step to approximate the fair price of an Asian option, and prove that the error of such an approximation using $N$ points converges at a rate close to $1/N$.

High dimensional integrals are abundant in many fields of research including quantum physics. The aim of this paper is to develop efficient recursive strategies to tackle a class of high dimensional integrals having a special product structure with low order couplings, motivated by models in lattice gauge theory from quantum field theory. A novel element of this work is the potential benefit in using lattice cubature rules. The group structure within lattice rules combined with the special structure in the physics integrands may allow efficient computations based on Fast Fourier Transforms. Applications to the quantum mechanical rotor and compact $U(1)$ lattice gauge theory in two and three dimensions are considered.

This paper deals with the kernel-based approximation of a multivariate periodic function by interpolation at the points of an integration lattice -- a setting that, as pointed out by Zeng, Leung, Hickernell (MCQMC2004, 2006) and Zeng, Kritzer, Hickernell (Constr. Approx., 2009), allows fast evaluation by fast Fourier transform, so avoiding the need for a linear solver. The main contribution of the paper is the application to the approximation problem for uncertainty quantification of elliptic partial differential equations, with the diffusion coefficient given by a random field that is periodic in the stochastic variables, in the model proposed recently by Kaarnioja, Kuo, Sloan (SIAM J. Numer. Anal., 2020). The paper gives a full error analysis, and full details of the construction of lattices needed to ensure a good (but inevitably not optimal) rate of convergence and an error bound independent of dimension. Numerical experiments support the theory., Comment: 37 pages, 5 figures

We propose and analyze a quasi-Monte Carlo (QMC) algorithm for efficient simulation of wave propagation modeled by the Helmholtz equation in a bounded region in which the refractive index is random and spatially heterogenous. Our focus is on the case in which the region can contain multiple wavelengths. We bypass the usual sign-indefiniteness of the Helmholtz problem by switching to an alternative sign-definite formulation recently developed by Ganesh and Morgenstern (Numerical Algorithms, 83, 1441-1487, 2020). The price to pay is that the regularity analysis required for QMC methods becomes much more technical. Nevertheless we obtain a complete analysis with error comprising stochastic dimension truncation error, finite element error and cubature error, with results comparable to those obtained for the diffusion problem.

We study an optimal control problem under uncertainty, where the target function is the solution of an elliptic partial differential equation with random coefficients, steered by a control function. The robust formulation of the optimization problem is stated as a high-dimensional integration problem over the stochastic variables. It is well known that carrying out a high-dimensional numerical integration of this kind using a Monte Carlo method has a notoriously slow convergence rate; meanwhile, a faster rate of convergence can potentially be obtained by using sparse grid quadratures, but these lead to discretized systems that are non-convex due to the involvement of negative quadrature weights. In this paper, we analyze instead the application of a quasi-Monte Carlo method, which retains the desirable convexity structure of the system and has a faster convergence rate compared to ordinary Monte Carlo methods. In particular, we show that under moderate assumptions on the decay of the input random field, the error rate obtained by using a specially designed, randomly shifted rank-1 lattice quadrature rule is essentially inversely proportional to the number of quadrature nodes. The overall discretization error of the problem, consisting of the dimension truncation error, finite element discretization error and quasi-Monte Carlo quadrature error, is derived in detail. We assess the theoretical findings in numerical experiments.

In a recent paper by the same authors, we provided a theoretical foundation for the component-by-component (CBC) construction of lattice algorithms for multivariate $L_2$ approximation in the worst case setting, for functions in a periodic space with general weight parameters. The construction led to an error bound that achieves the best possible rate of convergence for lattice algorithms. Previously available literature covered only weights of a simple form commonly known as product weights. In this paper we address the computational aspect of the construction. We develop fast CBC construction of lattice algorithms for special forms of weight parameters, including the so-called POD weights and SPOD weights which arise from PDE applications, making the lattice algorithms truly applicable in practice. With $d$ denoting the dimension and $n$ the number of lattice points, we show that the construction cost is $\mathcal{O}(d\,n\log(n) + d^2\log(d)\,n)$ for POD weights, and $\mathcal{O}(d\,n\log(n) + d^3\sigma^2\,n)$ for SPOD weights of degree $\sigma\ge 2$. The resulting lattice generating vectors can be used in other lattice-based approximation algorithms, including kernel methods or splines.

This paper provides the theoretical foundation for the construction of lattice algorithms for multivariate $L_2$ approximation in the worst case setting, for functions in a periodic space with general weight parameters. Our construction leads to an error bound that achieves the best possible rate of convergence for lattice algorithms. This work is motivated by PDE applications in which bounds on the norm of the functions to be approximated require special forms of weight parameters (so-called POD weights or SPOD weights), as opposed to the simple product weights covered by the existing literature. Our result can be applied to other lattice-based approximation algorithms, including kernel methods or splines.

We consider rank-1 lattices for integration and reconstruction of functions with series expansion supported on a finite index set. We explore the connection between the periodic Fourier space and the non-periodic cosine space and Chebyshev space, via tent transform and then cosine transform, to transfer known results from the periodic setting into new insights for the non-periodic settings. Fast discrete cosine transform can be applied for the reconstruction phase. To reduce the size of the auxiliary index set in the associated component-by-component (CBC) construction for the lattice generating vectors, we work with a bi-orthonormal set of basis functions, leading to three methods for function reconstruction in the non-periodic settings. We provide new theory and efficient algorithmic strategies for the CBC construction. We also interpret our results in the context of general function approximation and discrete least-squares approximation.

Many studies in uncertainty quantification have been carried out under the assumption of an input random field in which a countable number of independent random variables are each uniformly distributed on an interval, with these random variables entering linearly in the input random field (the so-called affine model). In this paper we consider an alternative model of the random field, in which the random variables have the same uniform distribution on an interval, but the random variables enter the input field as periodic functions. The field is constructed in such a way as to have the same mean and covariance function as the affine random field. Higher moments differ from the affine case, but in general the periodic model seems no less desirable. The new model of the random field is used to compute expected values of a quantity of interest arising from an elliptic PDE with random coefficients. The periodicity is shown to yield a higher order cubature convergence rate of $\mathcal{O}(n^{-1/p})$ independently of the dimension when used in conjunction with rank-1 lattice cubature rules constructed using suitably chosen smoothness-driven product and order dependent weights, where $n$ is the number of lattice points and $p$ is the summability exponent of the fluctuations in the series expansion of the random coefficient. We present numerical examples that assess the performance of our method., Comment: 24 pages, 6 figures

We seek shifted lattice rules that are good for high dimensional integration over the unit cube in the setting of an unanchored weighted Sobolev space of functions with square-integrable mixed first derivatives. Many existing studies rely on random shifting of the lattice, whereas here we work with lattice rules with a deterministic shift. Specifically, we consider "half-shifted" rules, in which each component of the shift is an odd multiple of $1/(2N)$, where $N$ is the number of points in the lattice. We show, by applying the principle that \emph{there is always at least one choice as good as the average}, that for a given generating vector there exists a half-shifted rule whose squared worst-case error differs from the shift-averaged squared worst-case error by a term of order only ${1/N^2}$. Numerical experiments, in which the generating vector is chosen component-by-component (CBC) as for randomly shifted lattices and then the shift by a new "CBC for shift" algorithm, yield encouraging results.

An existence result is presented for the worst-case error of lattice rules for high dimensional integration over the unit cube, in an unanchored weighted space of functions with square-integrable mixed first derivatives. Existing studies rely on random shifting of the lattice to simplify the analysis, whereas in this paper neither shifting nor any other form of randomisation is considered. Given that a certain number-theoretic conjecture holds, it is shown that there exists an $N$-point rank-one lattice rule which gives a worst-case error of order $1/\sqrt{N}$ up to a (dimension-independent) logarithmic factor. Numerical results suggest that the conjecture is plausible.

The component-by-component (CBC) algorithm is a method for constructing good generating vectors for lattice rules for the efficient computation of high-dimensional integrals in the "weighted" function space setting introduced by Sloan and Wo\'zniakowski. The "weights" that define such spaces are needed as inputs into the CBC algorithm, and so a natural question is, for a given problem how does one choose the weights? This paper introduces two new CBC algorithms which, given bounds on the mixed first derivatives of the integrand, produce a randomly shifted lattice rule with a guaranteed bound on the root-mean-square error. This alleviates the need for the user to specify the weights. We deal with "product weights" and "product and order dependent (POD) weights". Numerical tables compare the two algorithms under various assumed bounds on the mixed first derivatives, and provide rigorous upper bounds on the root-mean-square integration error.

We consider the forward problem of uncertainty quantification for the generalised Dirichlet eigenvalue problem for a coercive second order partial differential operator with random coefficients, motivated by problems in structural mechanics, photonic crystals and neutron diffusion. The PDE coefficients are assumed to be uniformly bounded random fields, represented as infinite series parametrised by uniformly distributed i.i.d. random variables. The expectation of the fundamental eigenvalue of this problem is computed by (a) truncating the infinite series which define the coefficients; (b) approximating the resulting truncated problem using lowest order conforming finite elements and a sparse matrix eigenvalue solver; and (c) approximating the resulting finite (but high dimensional) integral by a randomly shifted quasi-Monte Carlo lattice rule, with specially chosen generating vector. We prove error estimates for the combined error, which depend on the truncation dimension $s$, the finite element mesh diameter $h$, and the number of quasi-Monte Carlo samples $N$. Under suitable regularity assumptions, our bounds are of the particular form $\mathcal{O}(h^2+N^{-1+\delta})$, where $\delta>0$ is arbitrary and the hidden constant is independent of the truncation dimension, which needs to grow as $h\to 0$ and $N\to\infty$. Although the eigenvalue problem is nonlinear, which means it is generally considered harder than the analogous source problem, in almost all cases we obtain error bounds that converge at the same rate as the corresponding rate for the source problem. The proof involves a detailed study of the regularity of the fundamental eigenvalue as a function of the random parameters. As a key intermediate result in the analysis, we prove that the spectral gap (between the fundamental and the second eigenvalues) is uniformly positive over all realisations of the random problem.

In this paper we focus on efficient implementations of the Multivariate Decomposition Method (MDM) for approximating integrals of $\infty$-variate functions. Such $\infty$-variate integrals occur for example as expectations in uncertainty quantification. Starting with the anchored decomposition $f = \sum_{\mathfrak{u}\subset\mathbb{N}} f_\mathfrak{u}$, where the sum is over all finite subsets of $\mathbb{N}$ and each $f_\mathfrak{u}$ depends only on the variables $x_j$ with $j\in\mathfrak{u}$, our MDM algorithm approximates the integral of $f$ by first truncating the sum to some `active set' and then approximating the integral of the remaining functions $f_\mathfrak{u}$ term-by-term using Smolyak or (randomized) quasi-Monte Carlo (QMC) quadratures. The anchored decomposition allows us to compute $f_\mathfrak{u}$ explicitly by function evaluations of $f$. Given the specification of the active set and theoretically derived parameters of the quadrature rules, we exploit structures in both the formula for computing $f_\mathfrak{u}$ and the quadrature rules to develop computationally efficient strategies to implement the MDM in various scenarios. In particular, we avoid repeated function evaluations at the same point. We provide numerical results for a test function to demonstrate the effectiveness of the algorithm.

We show how simple kinks and jumps of otherwise smooth integrands over $\mathbb{R}^d$ can be dealt with by a preliminary integration with respect to a single well chosen variable. It is assumed that this preintegration, or conditional sampling, can be carried out with negligible error, which is the case in particular for option pricing problems. It is proven that under appropriate conditions the preintegrated function of $d-1$ variables belongs to appropriate mixed Sobolev spaces, so potentially allowing high efficiency of Quasi Monte Carlo and Sparse Grid Methods applied to the preintegrated problem. The efficiency of applying Quasi Monte Carlo to the preintegrated function are demonstrated on a digital Asian option using the Principal Component Analysis factorisation of the covariance matrix.

Building on previous research which generalized multilevel Monte Carlo methods using either sparse grids or Quasi-Monte Carlo methods, this paper considers the combination of all these ideas applied to elliptic PDEs with finite-dimensional uncertainty in the coefficients. It shows the potential for the computational cost to achieve an $O(\varepsilon)$ r.m.s. accuracy to be $O(\varepsilon^{-r})$ with $r<2$, independently of the spatial dimension of the PDE.

This article provides a high-level overview of some recent works on the application of quasi-Monte Carlo (QMC) methods to PDEs with random coefficients. It is based on an in-depth survey of a similar title by the same authors, with an accompanying software package which is also briefly discussed here. Embedded in this article is a step-by-step tutorial of the required analysis for the setting known as the uniform case with first order QMC rules. The aim of this article is to provide an easy entry point for QMC experts wanting to start research in this direction and for PDE analysts and practitioners wanting to tap into contemporary QMC theory and methods., Comment: arXiv admin note: text overlap with arXiv:1606.06613

This article provides an overview of some interfaces between the theory of quasi-Monte Carlo (QMC) methods and applications. We summarize three QMC theoretical settings: first order QMC methods in the unit cube $[0,1]^s$ and in $\mathbb{R}^s$, and higher order QMC methods in the unit cube. One important feature is that their error bounds can be independent of the dimension $s$ under appropriate conditions on the function spaces. Another important feature is that good parameters for these QMC methods can be obtained by fast efficient algorithms even when $s$ is large. We outline three different applications and explain how they can tap into the different QMC theory. We also discuss three cost saving strategies that can be combined with QMC in these applications. Many of these recent QMC theory and methods are developed not in isolation, but in close connection with applications.

In a previous paper (J. Comp. Phys. 230 (2011), 3668--3694), the authors proposed a new practical method for computing expected values of functionals of solutions for certain classes of elliptic partial differential equations with random coefficients. This method was based on combining quasi-Monte Carlo (QMC) methods for computing the expected values with circulant embedding methods for sampling the random field on a regular grid. It was found capable of handling fluid flow problems in random heterogeneous media with high stochastic dimension, but a convergence theory was missing. This paper provides a convergence analysis for the method in the case when the QMC method is a specially designed randomly shifted lattice rule. The convergence result depends on the eigenvalues of the underlying nested block circulant matrix and can be independent of the number of stochastic variables under certain assumptions. In fact the QMC analysis applies to general factorisations of the covariance matrix to sample the random field. The error analysis for the underlying fully discrete finite element method allows for locally refined meshes (via interpolation from a regular sampling grid of the random field). Numerical results on a non-regular domain with corner singularities in two spatial dimensions and on a regular domain in three spatial dimensions are included.

In this paper we prove, under mild conditions, that the positive definiteness of the circulant matrix appearing in the circulant embedding method is always guaranteed, provided the enclosing cube is sufficiently large. We examine in detail the case of the Mat\'ern covariance, and prove (for fixed correlation length) that, as $h_0\rightarrow 0$, positive definiteness is guaranteed when the random field is sampled on a cube of size order $(1 + \nu^{1/2} \log h_0^{-1})$ times larger than the size of the physical domain. (Here $h_0$ is the mesh spacing of the regular grid and $\nu$ the Mat\'ern smoothness parameter.) We show that the sampling cube can become smaller as the correlation length decreases when $h_0$ and $\nu$ are fixed. Our results are confirmed by numerical experiments. We prove several results about the decay of the eigenvalues of the circulant matrix. These lead to the conjecture, verified by numerical experiment, that they decay with the same rate as the Karhunen--Lo\`{e}ve eigenvalues of the covariance operator.

We analyze a new random algorithm for numerical integration of $d$-variate functions over $[0,1]^d$ from a weighted Sobolev space with dominating mixed smoothness $\alpha\ge 0$ and product weights $1\ge\gamma_1\ge\gamma_2\ge\cdots>0$, where the functions are continuous and periodic when $\alpha>1/2$. The algorithm is based on rank-$1$ lattice rules with a random number of points~$n$. For the case $\alpha>1/2$, we prove that the algorithm achieves almost the optimal order of convergence of $\mathcal{O}(n^{-\alpha-1/2})$, where the implied constant is independent of the dimension~$d$ if the weights satisfy $\sum_{j=1}^\infty \gamma_j^{1/\alpha}<\infty$. The same rate of convergence holds for the more general case $\alpha>0$ by adding a random shift to the lattice rule with random $n$. This shows, in particular, that the exponent of strong tractability in the randomized setting equals $1/(\alpha+1/2)$, if the weights decay fast enough. We obtain a lower bound to indicate that our results are essentially optimal. This paper is a significant advancement over previous related works with respect to the potential for implementation and the independence of error bounds on the problem dimension. Other known algorithms which achieve the optimal error bounds, such as those based on Frolov's method, are very difficult to implement especially in high dimensions. Here we adapt a lesser-known randomization technique introduced by Bakhvalov in 1961. This algorithm is based on rank-$1$ lattice rules which are very easy to implement given the integer generating vectors. A simple probabilistic approach can be used to obtain suitable generating vectors., Comment: 17 pages

It is known that the Brownian bridge or L\'evy-Ciesielski construction of Brownian paths almost surely converges uniformly to the true Brownian path. In the present article the focus is on the uniform error. In particular, we show constructively that at level $N$, at which there are $d=2^N$ points evaluated on the Brownian path, the uniform error and its square, and the uniform error of geometric Brownian motion, have upper bounds of order $\mathcal{O}(\sqrt{\ln d/d})$, matching the known orders. We apply the results to an option pricing example.

This article provides a survey of recent research efforts on the application of quasi-Monte Carlo (QMC) methods to elliptic partial differential equations (PDEs) with random diffusion coefficients. It considers, and contrasts, the uniform case versus the lognormal case, single-level algorithms versus multi-level algorithms, first order QMC rules versus higher order QMC rules, and deterministic QMC methods versus randomized QMC methods. It gives a summary of the error analysis and proof techniques in a unified view, and provides a practical guide to the software for constructing and generating QMC points tailored to the PDE problems. The analysis for the uniform case can be generalized to cover a range of affine parametric operator equations.

We develop algorithms for multivariate integration and approximation in the weighted half-period cosine space of smooth non-periodic functions. We use specially constructed tent-transformed rank-1 lattice points as cubature nodes for integration and as sampling points for approximation. For both integration and approximation, we study the connection between the worst-case errors of our algorithms in the cosine space and the worst-case errors of some related algorithms in the well-known weighted Korobov space of smooth periodic functions. By exploiting this connection, we are able to obtain constructive worst-case error bounds with good convergence rates for the cosine space., Comment: Journal of Complexity, Available online 26 May 2016

In this paper we present a rigorous cost and error analysis of a multilevel estimator based on randomly shifted Quasi-Monte Carlo (QMC) lattice rules for lognormal diffusion problems. These problems are motivated by uncertainty quantification problems in subsurface flow. We extend the convergence analysis in [Graham et al., Numer. Math. 2014] to multilevel Quasi-Monte Carlo finite element discretizations and give a constructive proof of the dimension-independent convergence of the QMC rules. More precisely, we provide suitable parameters for the construction of such rules that yield the required variance reduction for the multilevel scheme to achieve an $\varepsilon$-error with a cost of $\mathcal{O}(\varepsilon^{-\theta})$ with $\theta < 2$, and in practice even $\theta \approx 1$, for sufficiently fast decaying covariance kernels of the underlying Gaussian random field inputs. This confirms that the computational gains due to the application of multilevel sampling methods and the gains due to the application of QMC methods, both demonstrated in earlier works for the same model problem, are complementary. A series of numerical experiments confirms these gains. The results show that in practice the multilevel QMC method consistently outperforms both the multilevel MC method and the single-level variants even for non-smooth problems., Comment: 32 pages

Quasi-Monte Carlo (QMC) rules $1/N \sum_{n=0}^{N-1} f(\boldsymbol{y}_n A)$ can be used to approximate integrals of the form $\int_{[0,1]^s} f(\boldsymbol{y} A) \,\mathrm{d} \boldsymbol{y}$, where $A$ is a matrix and $\boldsymbol{y}$ is row vector. This type of integral arises for example from the simulation of a normal distribution with a general covariance matrix, from the approximation of the expectation value of solutions of PDEs with random coefficients, or from applications from statistics. In this paper we design QMC quadrature points $\boldsymbol{y}_0, ..., \boldsymbol{y}_{N-1} \in [0,1]^s$ such that for the matrix $Y = (\boldsymbol{y}_{0}^\top, ..., \boldsymbol{y}_{N-1}^\top)^\top$ whose rows are the quadrature points, one can use the fast Fourier transform to compute the matrix-vector product $Y \boldsymbol{a}^\top$, $\boldsymbol{a} \in \mathbb{R}^s$, in $\mathcal{O}(N \log N)$ operations and at most $s-1$ extra additions. The proposed method can be applied to lattice rules, polynomial lattice rules and a certain type of Korobov $p$-set. The approach is illustrated computationally by three numerical experiments. The first test considers the generation of points with normal distribution and general covariance matrix, the second test applies QMC to high-dimensional, affine-parametric, elliptic partial differential equations with uniformly distributed random coefficients, and the third test addresses Finite-Element discretizations of elliptic partial differential equations with high-dimensional, log-normal random input data. All numerical tests show a significant speed-up of the computation times of the fast QMC matrix method compared to a conventional implementation as the dimension becomes large.

We further develop the \emph{Multivariate Decomposition Method} (MDM) for the Lebesgue integration of functions of infinitely many variables $x_1,x_2,x_3,\ldots$ with respect to a corresponding product of a one dimensional probability measure. Although a number of concepts of infinite-dimensional integrals have been used in the literature, questions of uniqueness and compatibility have mostly not been studied. We show that, under appropriate convergence conditions, the Lebesgue integral equals the `anchored' integral, independently of the anchor. The MDM assumes that point values of $f_{\mathfrak{u}}$ are available for important subsets ${\mathfrak{u}}$, at some known cost. In this paper we introduce a new setting, in which it is assumed that each $f_{\mathfrak{u}}$ belongs to a normed space $F_{\mathfrak{u}}$, and that bounds $B_{\mathfrak{u}}$ on $\|f_{\mathfrak{u}}\|_{F_{\mathfrak{u}}}$ are known. This contrasts with the assumption in many papers that weights $\gamma_{\mathfrak{u}}$, appearing in the norm of the infinite-dimensional function space, are somehow known. Often such weights $\gamma_{\mathfrak{u}}$ were determined by minimizing an error bound depending on the $B_{\mathfrak{u}}$, the $\gamma_{\mathfrak{u}}$ \emph{and} the chosen algorithm, resulting in weights that depend on the algorithm. In contrast, in this paper only the bounds $B_{\mathfrak{u}}$ are assumed known. We give two examples in which we specialize the MDM: in the first case $F_{\mathfrak{u}}$ is the $|{\mathfrak{u}}|$-fold tensor product of an anchored reproducing kernel Hilbert space, and in the second case it is a particular non-Hilbert space for integration over an unbounded domain.

We construct quasi-Monte Carlo methods to approximate the expected values of linear functionals of Galerkin discretizations of parametric operator equations which depend on a possibly infinite sequence of parameters. Such problems arise in the numerical solution of differential and integral equations with random field inputs. We analyze the regularity of the solutions with respect to the parameters in terms of the rate of decay of the fluctuations of the input field. If $p\in (0,1]$ denotes the "summability exponent" corresponding to the fluctuations in affine-parametric families of operators, then we prove that deterministic "interlaced polynomial lattice rules" of order $\alpha = \lfloor 1/p \rfloor+1$ in $s$ dimensions with $N$ points can be constructed using a fast component-by-component algorithm, in $\mathcal{O}(\alpha\,s\, N\log N + \alpha^2\,s^2 N)$ operations, to achieve a convergence rate of $\mathcal{O}(N^{-1/p})$, with the implied constant independent of $s$. This dimension-independent convergence rate is superior to the rate $\mathcal{O}(N^{-1/p+1/2})$, for $2/3\leq p\leq 1$ recently established for randomly shifted lattice rules under comparable assumptions. In our analysis we use a non-standard Banach space setting and introduce "smoothness-driven product and order dependent (SPOD)" weights for which we show fast CBC construction., Comment: Equation (3.7) in published version is not correct as stated. This file contains the correct formula

Quasi-Monte Carlo (QMC) methods are applied to multi-level Finite Element (FE) discretizations of elliptic partial differential equations (PDEs) with a random coefficient, to estimate expected values of linear functionals of the solution. The expected value is considered as an infinite-dimensional integral in the parameter space corresponding to the randomness induced by the random coefficient. We use a multi-level algorithm, with the number of QMC points depending on the discretization level, and with a level-dependent dimension truncation strategy. In some scenarios, we show that the overall error is $\mathcal{O}(h^2)$, where $h$ is the finest FE mesh width, or $\mathcal{O}(N^{-1+\delta})$ for arbitrary $\delta>0$, where $N$ is the maximal number of QMC sampling points. For these scenarios, the total work is essentially of the order of one single PDE solve at the finest FE discretization level. The analysis exploits regularity of the parametric solution with respect to both the physical variables (the variables in the physical domain) and the parametric variables (the parameters corresponding to randomness). Families of QMC rules with "POD weights" ("product and order dependent weights") which quantify the relative importance of subsets of the variables are found to be natural for proving convergence rates of QMC errors that are independent of the number of parametric variables.