Alexander D. Gilbert

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.

Alexander D. Gilbert, Greg W. Wasilkowski

Constructing active sets is a key part of the Multivariate Decomposition Method. An algorithm for constructing optimal or quasi-optimal active sets is proposed in the paper. By numerical experiments, it is shown that the new method can provide sets that are significantly smaller than the sets constructed by the already existing method. The experiments also show that the superposition dimension could surprisingly be very small, at most 3, when the error demand is not smaller than $10^{-3}$ and the weights decay sufficiently fast.

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.