Dirk Nuyens

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.

Josef Dick, Dirk Nuyens, Friedrich Pillichshammer

The aim of this paper is to show that one can achieve convergence rates of $N^{-α+ δ}$ for $α> 1/2$ (and for $δ> 0$ arbitrarily small) for nonperiodic $α$-smooth cosine series using lattice rules without random shifting. The smoothness of the functions can be measured by the decay rate of the cosine coefficients. For a specific choice of the parameters the cosine series space coincides with the unanchored Sobolev space of smoothness 1. We study the embeddings of various reproducing kernel Hilbert spaces and numerical integration in the cosine series function space and show that by applying the so-called tent transformation to a lattice rule one can achieve the (almost) optimal rate of convergence of the integration error. The same holds true for symmetrized lattice rules for the tensor product of the direct sum of the Korobov space and cosine series space, but with a stronger dependence on the dimension in this case.

Jan Baldeaux, Josef Dick, Gunther Leobacher et al.

We show how to obtain a fast component-by-component construction algorithm for higher order polynomial lattice rules. Such rules are useful for multivariate quadrature of high-dimensional smooth functions over the unit cube as they achieve the near optimal order of convergence. The main problem addressed in this paper is to find an efficient way of computing the worst-case error. A general algorithm is presented and explicit expressions for base~2 are given. To obtain an efficient component-by-component construction algorithm we exploit the structure of the underlying cyclic group. We compare our new higher order multivariate quadrature rules to existing quadrature rules based on higher order digital nets by computing their worst-case error. These numerical results show that the higher order polynomial lattice rules improve upon the known constructions of quasi-Monte Carlo rules based on higher order digital nets.

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.

Pieterjan Robbe, Dirk Nuyens, Stefan Vandewalle

The Multilevel Monte Carlo method is an efficient variance reduction technique. It uses a sequence of coarse approximations to reduce the computational cost in uncertainty quantification applications. The method is nowadays often considered to be the method of choice for solving PDEs with random coefficients when many uncertainties are involved. When using Full Multigrid to solve the deterministic problem, coarse solutions obtained by the solver can be recycled as samples in the Multilevel Monte Carlo method, as was pointed out by Kumar, Oosterlee and Dwight [Int. J. Uncertain. Quantif., 7 (2017), pp. 57--81]. In this article, an alternative approach is considered, using Quasi-Monte Carlo points, to speed up convergence. Additionally, our method comes with an improved variance estimate which is also valid in case of the Monte Carlo based approach. The new method is illustrated on the example of an elliptic PDE with lognormal diffusion coefficient. Numerical results for a variety of random fields with different smoothness parameters in the Matérn covariance function show that sample recycling is more efficient when the input random field is nonsmooth.

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.

Nico Achtsis, Ronald Cools, Dirk Nuyens

We develop a conditional sampling scheme for pricing knock-out barrier options under the Linear Transformations (LT) algorithm from Imai and Tan (2006). We compare our new method to an existing conditional Monte Carlo scheme from Glasserman and Staum (2001), and show that a substantial variance reduction is achieved. We extend the method to allow pricing knock-in barrier options and introduce a root-finding method to obtain a further variance reduction. The effectiveness of the new method is supported by numerical results.

Nico Achtsis, Ronald Cools, Dirk Nuyens

We propose a quasi-Monte Carlo algorithm for pricing knock-out and knock-in barrier options under the Heston (1993) stochastic volatility model. This is done by modifying the LT method from Imai and Tan (2006) for the Heston model such that the first uniform variable does not influence the stochastic volatility path and then conditionally modifying its marginals to fulfill the barrier condition(s). We show this method is unbiased and never does worse than the unconditional algorithm. Additionally the conditioning is combined with a root finding method to also force positive payouts. The effectiveness of this method is shown by extensive numerical results.

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.

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.

Dirk Nuyens, Gowri Suryanarayana, Markus Weimar

We study multivariate integration of functions that are invariant under permutations (of subsets) of their arguments. We find an upper bound for the $n$th minimal worst case error and show that under certain conditions, it can be bounded independent of the number of dimensions. In particular, we study the application of unshifted and randomly shifted rank-$1$ lattice rules in such a problem setting. We derive conditions under which multivariate integration is polynomially or strongly polynomially tractable with the Monte Carlo rate of convergence $O(n^{-1/2})$. Furthermore, we prove that those tractability results can be achieved with shifted lattice rules and that the shifts are indeed necessary. Finally, we show the existence of rank-$1$ lattice rules whose worst case error on the permutation- and shift-invariant spaces converge with (almost) optimal rate. That is, we derive error bounds of the form $O(n^{-λ/2})$ for all $1 \leq λ< 2 α$, where $α$ denotes the smoothness of the spaces. Keywords: Numerical integration, Quadrature, Cubature, Quasi-Monte Carlo methods, Rank-1 lattice rules.

Adrian Ebert, Hernan Leövey, Dirk Nuyens

The (fast) component-by-component (CBC) algorithm is an efficient tool for the construction of generating vectors for quasi-Monte Carlo rank-1 lattice rules in weighted reproducing kernel Hilbert spaces. We consider product weights, which assigns a weight to each dimension. These weights encode the effect a certain variable (or a group of variables by the product of the individual weights) has. Smaller weights indicate less importance. Kuo (2003) proved that the CBC algorithm achieves the optimal rate of convergence in the respective function spaces, but this does not imply the algorithm will find the generating vector with the smallest worst-case error. In fact it does not. We investigate a generalization of the component-by-component construction that allows for a general successive coordinate search (SCS), based on an initial generating vector, and with the aim of getting closer to the smallest worst-case error. The proposed method admits the same type of worst-case error bounds as the CBC algorithm, independent of the choice of the initial vector. Under the same summability conditions on the weights as in [Kuo,2003] the error bound of the algorithm can be made independent of the dimension $d$ and we achieve the same optimal order of convergence for the function spaces from [Kuo,2003]. Moreover, a fast version of our method, based on the fast CBC algorithm by Nuyens and Cools, is available, reducing the computational cost of the algorithm to $O(d \, n \log(n))$ operations, where $n$ denotes the number of function evaluations. Numerical experiments seeded by a Korobov-type generating vector show that the new SCS algorithm will find better choices than the CBC algorithm and the effect is better when the weights decay slower.

Dirk Nuyens, Gowri Suryanarayana, Markus Weimar

We study multivariate integration of functions that are invariant under the permutation (of a subset) of their arguments. Recently, in Nuyens, Suryanarayana, and Weimar (Adv. Comput. Math. (2016), 42(1):55--84), the authors derived an upper estimate for the $n$th minimal worst case error for such problems, and showed that under certain conditions this upper bound only weakly depends on the dimension. We extend these results by proposing two (semi-) explicit construction schemes. We develop a component-by-component algorithm to find the generating vector for a shifted rank-$1$ lattice rule that obtains a rate of convergence arbitrarily close to $\mathcal{O}(n^{-α})$, where $α>1/2$ denotes the smoothness of our function space and $n$ is the number of cubature nodes. Further, we develop a semi-constructive algorithm that builds on point sets which can be used to approximate the integrands of interest with a small error; the cubature error is then bounded by the error of approximation. Here the same rate of convergence is achieved while the dependence of the error bounds on the dimension $d$ is significantly improved.

Dong T. P. Nguyen, Dirk Nuyens

In this paper we analyze the approximation of multivariate integrals over the Euclidean plane for functions which are analytic. We show explicit upper bounds which attain the exponential rate of convergence. We use an infinite grid with different mesh sizes and lengths in each direction to sample the function, and then truncate it. In our analysis, the mesh sizes and the truncated domain are chosen by optimally balancing the truncation error and the discretization error. This paper derives results in comparable function space settings, extended to $\R^s$, as which were recently obtained in the unit cube by Dick, Larcher, Pillichshammer and Wo{ź}niakowski (2011). They showed that both lattice rules and regular grids, with different mesh sizes in each direction, attain exponential rates, hence motivating us to analyze only cubature formula based on regular meshes. We further also amend the analysis of older publications, e.g., Sloan and Osborn (1987) and Sugihara (1987), using lattice rules on $\R^s$ by taking the truncation error into account and extending them to take the anisotropy of the function space into account.

Adrian Ebert, Peter Kritzer, Dirk Nuyens

In the analysis of using quasi-Monte Carlo (QMC) methods to approximate expectations of a linear functional of the solution of an elliptic PDE with random diffusion coefficient the sensitivity w.r.t. the parameters is often stated in terms of product-and-order-dependent (POD) weights. The (offline) fast component-by-component (CBC) construction of an $N$-point QMC method making use of these POD weights leads to a cost of $\mathcal{O}(s N \log(N) + s^2 N)$ with $s$ the parameter truncation dimension. When $s$ is large this cost is prohibitive. As an alternative Herrmann and Schwab introduced an analysis resulting in product weights to reduce the construction cost to $\mathcal{O}(s N \log(N))$. We here show how the reduced CBC method can be used for POD weights to reduce the cost to $\mathcal{O}(\sum_{j=1}^{\min\{s,s^*\}} (m-w_j+j) \, b^{m-w_j})$, where $N=b^m$ with prime $b$, $w_1 \le \cdots \le w_s$ are nonnegative integers and $s^*$ can be chosen much smaller than $s$ depending on the regularity of the random field expansion as such making it possible to use the POD weights directly. We show a total error estimate for using randomly shifted lattice rules constructed through the reduced CBC construction.

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.

Yiqing Zhou, Karsten Naert, Dirk Nuyens

Accurate knowledge of scrap composition can increase the usage of recycled material to produce steel, reducing the need for raw ore extraction and minimizing environmental impact by conserving natural resources and lowering carbon emissions. First, we introduce two state-space models for the elemental composition of scrap in Electric Arc Furnaces (EAF) and Basic Oxygen Furnaces (BOF): a linear model for elements that transfer entirely into steel, and a non-linear model for elements that partition between steel and slag. The models are fitted with the Kalman filter and the unscented Kalman filter, respectively, using only data already collected in the standard steel production process. Crucially, the resulting scrap composition estimates can in turn be used to predict the elemental composition of future steel production. Second, we analyze how key hyperparameters affect estimation accuracy and stability, and we provide practical guidelines for tuning them from expert knowledge and historical data. Third, we validate the models on real BOF data from ArcelorMittal, using Cu and Cr as representative elements. Both filters outperform windowed non-negative least squares regression, a strong baseline method for scrap composition estimation, yielding reliable real-time estimates of scrap composition.

Yuya Suzuki, Dirk Nuyens

In this paper, we propose a numerical method to approximate the solution of the time-dependent Schrödinger equation with periodic boundary condition in a high-dimensional setting. We discretize space by using the Fourier pseudo-spectral method on rank-$1$ lattice points, and then discretize time by using a higher-order exponential operator splitting method. In this scheme the convergence rate of the time discretization depends on properties of the spatial discretization. We prove that the proposed method, using rank-$1$ lattice points in space, allows to obtain higher-order time convergence, and, additionally, that the necessary condition on the space discretization can be independent of the problem dimension $d$. We illustrate our method by numerical results from 2 to 8 dimensions which show that such higher-order convergence can really be obtained in practice.

Dirk Nuyens, Ronald Cools

In a series of papers, in 1993, 1994 & 1996, Sloan & Niederreiter introduced a modification of lattice rules for non-periodic functions, called "vertex modified lattice rules"', and a particular breed called "optimal vertex modified lattice rules". In the 1994 paper, Niederreiter & Sloan concentrate explicitly on Fibonacci lattice rules, which are a particular good choice of 2-dimensional lattice rules. Error bounds in this series of papers were given related to the star discrepancy. In this paper we pose the problem in terms of the so-called unanchored Sobolev space, which is a reproducing kernel Hilbert space often studied nowadays in which functions have $L_2$-integrable mixed first derivatives. It is known constructively that randomly shifted lattice rules, as well as deterministic tent-transformed lattice rules and deterministic fully symmetrized lattice rules can achieve close to $O(N^{-1})$ convergence in this space, see Sloan, Kuo & Joe (2002) and Dick, Nuyens & Pillichshammer (2014) respectively. We derive a break down of the worst-case error of vertex modified lattice rules in the unanchored Sobolev space in terms of the worst-case error in a Korobov space, a multilinear space and some additional "mixture term". For the 1-dimensional case this worst-case error is obvious and gives an explicit expression for the trapezoidal rule. In the 2-dimensional case this mixture term also takes on an explicit form for which we derive upper and lower bounds. For this case we prove that there exist lattice rules with a nice worst-case error bound with the additional mixture term of the form $N^{-1} \log^2(N)$.

Pieterjan Robbe, Dirk Nuyens, Stefan Vandewalle

We present an adaptive version of the Multi-Index Monte Carlo method, introduced by Haji-Ali, Nobile and Tempone (2016), for simulating PDEs with coefficients that are random fields. A classical technique for sampling from these random fields is the Karhunen-Loève expansion. Our adaptive algorithm is based on the adaptive algorithm used in sparse grid cubature as introduced by Gerstner and Griebel (2003), and automatically chooses the number of terms needed in this expansion, as well as the required spatial discretizations of the PDE model. We apply the method to a simplified model of a heat exchanger with random insulator material, where the stochastic characteristics are modeled as a lognormal random field, and we show consistent computational savings.

Pieterjan Robbe, Dirk Nuyens, Stefan Vandewalle

We present a Multi-Index Quasi-Monte Carlo method for the solution of elliptic partial differential equations with random coefficients. By combining the multi-index sampling idea with randomly shifted rank-1 lattice rules, the algorithm constructs an estimator for the expected value of some functional of the solution. The efficiency of this new method is illustrated on a three-dimensional subsurface flow problem with lognormal diffusion coefficient with underlying Matérn covariance function. This example is particularly challenging because of the small correlation length considered, and thus the large number of uncertainties that must be included. We show numerical evidence that it is possible to achieve a cost inversely proportional to the requested tolerance on the root-mean-square error, for problems with a smoothly varying random field

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.