Dirk Nuyens

Alexander D. Gilbert, Frances Y. Kuo, Dirk Nuyens et al.

We study the application of a quasi-Monte Carlo (QMC) method to a class of semi-linear parabolic reaction-diffusion partial differential equations used to model tumor growth. Mathematical models of tumor growth are largely phenomenological in nature, capturing infiltration of the tumor into surrounding healthy tissue, proliferation of the existing tumor, and patient response to therapies, such as chemotherapy and radiotherapy. Considerable inter-patient variability, inherent heterogeneity of the disease, sparse and noisy data collection, and model inadequacy all contribute to significant uncertainty in the model parameters. It is crucial that these uncertainties can be efficiently propagated through the model to compute quantities of interest (QoIs), which in turn may be used to inform clinical decisions. We show that QMC methods can be successful in computing expectations of meaningful QoIs. Well-posedness results are developed for the model and used to show a theoretical error bound for the case of uniform random fields. The theoretical linear error rate, which is superior to that of standard Monte Carlo, is verified numerically. Encouraging computational results are also provided for lognormal random fields, prompting further theoretical development.

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.

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.

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.

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.

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.

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.

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.

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.