Frances Y. Kuo

Josef Dick, Frances Y. Kuo, Quoc T. Le Gia et al.

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 $α= \lfloor 1/p \rfloor+1$ in $s$ dimensions with $N$ points can be constructed using a fast component-by-component algorithm, in $\mathcal{O}(α\,s\, N\log N + α^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.

Ivan G. Graham, Frances Y. Kuo, Dirk Nuyens et al.

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érn 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 + ν^{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 $ν$ the Matérn smoothness parameter.) We show that the sampling cube can become smaller as the correlation length decreases when $h_0$ and $ν$ 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ève eigenvalues of the covariance operator.

Frances Y. Kuo, Dirk Nuyens

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.

Ivan G. Graham, Frances Y. Kuo, Dirk Nuyens et al.

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.

Andreas Griewank, Frances Y. Kuo, Hernan Leövey et al.

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.

Frances Y. Kuo, Christoph Schwab, Ian H. Sloan

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+δ})$ for arbitrary $δ>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.

Josef Dick, Frances Y. Kuo, Quoc T. Le Gia et al.

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.

Michael B. Giles, Frances Y. Kuo, Ian H. Sloan

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.

Ronald Cools, Frances Y. Kuo, Dirk Nuyens et al.

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.

Peter Kritzer, Frances Y. Kuo, Dirk Nuyens et al.

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 $α\ge 0$ and product weights $1\geγ_1\geγ_2\ge\cdots>0$, where the functions are continuous and periodic when $α>1/2$. The algorithm is based on rank-$1$ lattice rules with a random number of points~$n$. For the case $α>1/2$, we prove that the algorithm achieves almost the optimal order of convergence of $\mathcal{O}(n^{-α-1/2})$, where the implied constant is independent of the dimension~$d$ if the weights satisfy $\sum_{j=1}^\infty γ_j^{1/α}<\infty$. The same rate of convergence holds for the more general case $α>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/(α+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.

Frances Y. Kuo, Dirk Nuyens

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.

Yoshihito Kazashi, Frances Y. Kuo, Ian H. Sloan

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.

Frances Y. Kuo, Dirk Nuyens

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.

Yoshihito Kazashi, Frances Y. Kuo, Ian H. Sloan

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.

Alexander Keller, Frances Y. Kuo, Dirk Nuyens et al.

This survey article is concerned with the application of lattice rules to Deep Neural Networks (DNNs), lattice rules being a family of quasi-Monte Carlo methods. They have demonstrated effectiveness in various contexts for high-dimensional integration and function approximation. They are extremely easy to implement thanks to their very simple formulation -- all that is required is a good integer generating vector of length matching the dimensionality of the problem. In recent years there has been a burst of research activities on the application and theory of DNNs. We review our recent article on using lattice rules as training points for DNNs with a smooth activation function, where we obtained explicit regularity bounds of the DNNs. By imposing restrictions on the network parameters to match the regularity features of the target function, we prove that DNNs with tailored lattice training points can achieve good theoretical generalization error bounds, with implied constants independent of the input dimension. We also demonstrate numerically that DNNs trained with our tailored regularization perform significantly better than with standard $\ell_2$ regularization.

Alexander D. Gilbert, Ivan G. Graham, Frances Y. Kuo et al.

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+δ})$, where $δ>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.

Alexander D. Gilbert, Frances Y. Kuo, Ian H. Sloan

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źniakowski. 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.

Frances Y. Kuo, Dirk Nuyens, Leszek Plaskota et al.

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 $γ_{\mathfrak{u}}$, appearing in the norm of the infinite-dimensional function space, are somehow known. Often such weights $γ_{\mathfrak{u}}$ were determined by minimizing an error bound depending on the $B_{\mathfrak{u}}$, the $γ_{\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.

Frances Y. Kuo, Robert Scheichl, Christoph Schwab et al.

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^{-θ})$ with $θ< 2$, and in practice even $θ\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.