Ian H. Sloan

Johann S. Brauchart, Edward B. Saff, Ian H. Sloan et al.

We study equal weight numerical integration, or Quasi Monte Carlo (QMC) rules, for functions in a Sobolev space $H^s(S^d)$ with smoothness parameter $s>d/2$ defined over the unit sphere $S^d$ in $R^{d+1}$. Focusing on $N$-point sets that achieve optimal order QMC error bounds (as is the case for efficient spherical designs), we are led to introduce the concept of QMC designs: these are sequences of $N$-point node sets $X_N$ on $S^d$ such that the worst-case error of the corresponding QMC rules satisfy a bound of order $O(N^{-s/d})$ as $N\to\infty$ with an implied constant that depends on the $H^s(S^d)$-norm. We provide methods for generation and numerical testing of QMC designs. As a consequence of a recent result of Bondarenko et al. on the existence of spherical designs with appropriate number of points, we show that minimizers of the $N$-point energy for the reproducing kernel for $H^s(S^d)$, $s>d/2$, form a sequence of QMC designs for $H^s(S^d)$. Furthermore, without appealing to the Bondarenko et al. result, we prove that point sets that maximize the sum of suitable powers of the Euclidean distance between pairs of points form a sequence of QMC designs for $H^s(S^d)$ with $s\in(d/2,d/2+1)$. Numerical experiments suggest that many familiar sequences of point sets on the sphere (equal area, spiral, minimal [Coulomb or log.] energy, and Fekete points) are QMC designs for appropriate values of $s$. For comparison purposes we show that sets of random points that are independently and uniformly distributed on the sphere do not constitute QMC designs for any $s>d/2$. If $(X_N)$ is a sequence of QMC designs for $H^s(S^d)$, we prove that it is also a sequence of QMC designs for $\mathbb{H}^{s'}(S^d)$ for all $s'\in(d/2,s)$. This leads to the question of determining the supremum of such $s$, for which we provide estimates based on computations for the aforementioned sequences.

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.

Yu Guang Wang, Quoc T. Le Gia, Ian H. Sloan et al.

Spherical needlets are highly localized radial polynomials on the sphere $\mathbb{S}^{d}\subset \mathbb{R}^{d+1}$, $d\ge 2$, with centers at the nodes of a suitable cubature rule. The original semidiscrete spherical needlet approximation of Narcowich, Petrushev and Ward is not computable, in that the needlet coefficients depend on inner product integrals. In this work we approximate these integrals by a second quadrature rule with an appropriate degree of precision, to construct a fully discrete needlet approximation. We prove that the resulting approximation is equivalent to filtered hyperinterpolation, that is to a filtered Fourier-Laplace series partial sum with inner products replaced by appropriate cubature sums. It follows that the $\mathbb{L}_{p}$-error of discrete needlet approximation of order $J$ for $1 \le p \le \infty$ and $s > d/p$ has for a function $f$ in the Sobolev space $\mathbb{W}_{p}^{s}(\mathbb{S}^{d})$ the optimal rate of convergence in the sense of optimal recovery, namely $\mathcal{O}\bigl(2^{-J s}\bigr)$. Moreover, this is achieved with a filter function that is of smoothness class $C^{\lfloor \frac{d+3}{2}\rfloor}$, in contrast to the usually assumed $C^{\infty}$. A numerical experiment for a class of functions in known Sobolev smoothness classes gives $\mathbb{L}_2$ errors for the fully discrete needlet approximation that are almost identical to those for the original semidiscrete needlet approximation. Another experiment uses needlets over the whole sphere for the lower levels together with high-level needlets with centers restricted to a local region. The resulting errors are reduced in the local region away from the boundary, indicating that local refinement in special regions is a promising strategy.

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.

Michael Feischl, Frances Kuo, Ian H. Sloan

We use the $H$-matrix technology to compute the approximate square root of a covariance matrix in linear cost. This allows us to generate normal and log-normal random fields on general point sets with optimal cost. We derive rigorous error estimates which show convergence of the method. Our approach requires only mild assumptions on the covariance function and on the point set. Therefore, it might be also a nice alternative to the circulant embedding approach which applies only to regular grids and stationary covariance functions.

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.

Quoc T. Le Gia, Ian H. Sloan, Yu Guang Wang et al.

In this paper we establish a multiscale approximation for random fields on the sphere using spherical needlets --- a class of spherical wavelets. We prove that the semidiscrete needlet decomposition converges in mean and pointwise senses for weakly isotropic random fields on $\mathbb{S}^{d}$, $d\ge2$. For numerical implementation, we construct a fully discrete needlet approximation of a smooth $2$-weakly isotropic random field on $\mathbb{S}^{d}$ and prove that the approximation error for fully discrete needlets has the same convergence order as that for semidiscrete needlets. Numerical examples are carried out for fully discrete needlet approximations of Gaussian random fields and compared to a discrete version of the truncated Fourier expansion.

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.

Q. T. Le Gia, Ian H. Sloan, Andrew J. Wathen

This paper proposes a new preconditioning scheme for a linear system with a saddle-point structure arising from a hybrid approximation scheme on the sphere, an approximation scheme that combines (local) spherical radial basis functions and (global) spherical polynomials. Making use of a recently derived inf-sup condition [13] and the Brezzi stability and convergence theorem for this approximation scheme, we show that the linear system can be optimally preconditioned with a suitable block-diagonal preconditioner. Numerical experiments with a non-uniform distribution of data points support the theoretical conclusions.

Quoc T. Le Gia, Ian H. Sloan, Robert S. Womersley et al.

This paper discusses sparse isotropic regularization for a random field on the unit sphere $\mathbb{S}^2$ in $\mathbb{R}^{3}$, where the field is expanded in terms of a spherical harmonic basis. A key feature is that the norm used in the regularization term, a hybrid of the $\ell_{1}$ and $\ell_2$-norms, is chosen so that the regularization preserves isotropy, in the sense that if the observed random field is strongly isotropic then so too is the regularized field. The Pareto efficient frontier is used to display the trade-off between the sparsity-inducing norm and the data discrepancy term, in order to help in the choice of a suitable regularization parameter. A numerical example using Cosmic Microwave Background (CMB) data is considered in detail. In particular, the numerical results explore the trade-off between regularization and discrepancy, and show that substantial sparsity can be achieved along with small $L_{2}$ error.

Hui Cao, Sergei V. Pereverzyev, Ian H. Sloan et al.

In this paper, a two-step regularization method is used to solve an ill-posed spherical pseudo-differential equation in the presence of noisy data. For the first step of regularization we approximate the data by means of a spherical polynomial that minimizes a functional with a penalty term consisting of the squared norm in a Sobolev space. The second step is a regularized collocation method. An error bound is obtained in the uniform norm, which is potentially smaller than that for either the noise reduction alone or the regularized collocation alone. We discuss an a posteriori parameter choice, and present some numerical experiments, which support the claimed superiority of the two-step method.

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.

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.

Johann S. Brauchart, Josef Dick, Edward B. Saff et al.

We prove that the covering radius of an $N$-point subset $X_N$ of the unit sphere $S^d \subset R^{d+1}$ is bounded above by a power of the worst-case error for equal weight cubature $\frac{1}{N}\sum_{\mathbf{x} \in X_N}f(\mathbf{x}) \approx \int_{S^d} f \, \mathrm{d} σ_d$ for functions in the Sobolev space $\mathbb{W}_p^s(S^d)$, where $σ_d$ denotes normalized area measure on $S^d.$ These bounds are close to optimal when $s$ is close to $d/p$. Our study of the worst-case error along with results of Brandolini et al. motivate the definition of Quasi-Monte Carlo (QMC) design sequences for $\mathbb{W}_p^s(S^d)$, which have previously been introduced only in the Hilbert space setting $p=2$. We say that a sequence $(X_N)$ of $N$-point configurations is a QMC-design sequence for $\mathbb{W}_p^s(S^d)$ with $s > d/p$ provided the worst-case equal weight cubature error for $X_N$ has order $N^{-s/d}$ as $N \to \infty$, a property that holds, in particular, for a sequence of spherical $t$-designs in which each design has order $t^d$ points. For the case $p = 1$, we deduce that any QMC-design sequence $(X_N)$ for $\mathbb{W}_1^s(S^d)$ with $s > d$ has the optimal covering property; i.e., the covering radius of $X_N$ has order $N^{-1/d}$ as $N \to \infty$. A significant portion of our effort is devoted to the formulation of the worst-case error in terms of a Bessel kernel, and showing that this kernel satisfies a Bernstein type inequality involving the mesh ratio of $X_N$. As a consequence we prove that any QMC-design sequence for $\mathbb{W}_p^s(S^d)$ is also a QMC-design sequence for $\mathbb{W}_{p^\prime}^s(S^d)$ for all $1 \leq p < p^\prime \leq \infty$ and, furthermore, if $(X_N)$ is a quasi-uniform QMC-design sequence for $\mathbb{W}_p^s(S^d)$, then it is also a QMC-design sequence for $\mathbb{W}_p^{s^\prime}(S^d)$ for all $s > s^\prime > d/p$.

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.