MDFEM: Multivariate decomposition finite element method for elliptic PDEs with lognormal diffusion coefficients using higher-order QMC and FEM

We introduce the multivariate decomposition finite element method for elliptic PDEs with lognormal diffusion coefficient $a=\exp(Z)$ where $Z$ is a Gaussian random field defined by an infinite series expansion $Z(\boldsymbol{y}) = \sum_{j \ge 1} y_j \, \phi_j$ with $y_j \sim \mathcal{N}(0,1)$ and a given sequence of functions $\{\phi_j\}_{j \ge 1}$. We use the MDFEM to approximate the expected value of a linear functional of the solution of the PDE which is an infinite-dimensional integral over the parameter space. The proposed algorithm uses the multivariate decomposition method to compute the infinite-dimensional integral by a decomposition into finite-dimensional integrals, which we resolve using quasi-Monte Carlo methods, and for which we use the finite element method to solve different instances of the PDE. We develop higher-order quasi-Monte Carlo rules for integration over the finite-dimensional Euclidean space with respect to the Gaussian distribution by use of a truncation strategy. By linear transformations of interlaced polynomial lattice rules from the unit cube to a multivariate box of the Euclidean space we achieve higher-order convergence rates for functions belonging to a class of anchored Gaussian Sobolev spaces while taking into account the truncation error. Under appropriate conditions, the MDFEM achieves higher-order convergence rates in term of error versus cost, i.e., to achieve an accuracy of $O(\epsilon)$ the computational cost is $O(\epsilon^{-1/\lambda-d'/\lambda}) = O(\epsilon^{-(p^* + d'/\tau)/(1-p^*)})$ where $\epsilon^{-1/\lambda}$ and $\epsilon^{-d'/\lambda}$ are respectively the cost of the quasi-Monte Carlo cubature and the finite element approximations, with $d' = d \, (1+\delta')$ for some $\delta' \ge 0$ and $d$ the physical dimension, and $0 < p^* \le (2 + d'/\tau)^{-1}$ is a parameter representing the sparsity of $\{\phi_j\}_{j \ge 1}$.Comment: 48 page

Univariate interpolation by exponential functions and gaussian RBFs for generic sets of nodes

We consider interpolation of univariate functions on arbitrary sets of nodes by Gaussian radial basis functions or by exponential functions. We derive closed-form expressions for the interpolation error based on the Harish-Chandra-Itzykson-Zuber formula. We then prove the exponential convergence of interpolation for functions analytic in a sufficiently large domain. As an application, we prove the global exponential convergence of optimization by expected improvement for such functions.Comment: Some stylistic improvements and added references following feedback from the reviewer

Hot new directions for quasi-Monte Carlo research in step with applications

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

Pareto Smoothed Importance Sampling

Importance weighting is a general way to adjust Monte Carlo integration to account for draws from the wrong distribution, but the resulting estimate can be noisy when the importance ratios have a heavy right tail. This routinely occurs when there are aspects of the target distribution that are not well captured by the approximating distribution, in which case more stable estimates can be obtained by modifying extreme importance ratios. We present a new method for stabilizing importance weights using a generalized Pareto distribution fit to the upper tail of the distribution of the simulated importance ratios. The method, which empirically performs better than existing methods for stabilizing importance sampling estimates, includes stabilized effective sample size estimates, Monte Carlo error estimates and convergence diagnostics.Comment: Major revision: 1) proofs for consistency, finite variance, and asymptotic normality, 2) justification of k<0.7 with theoretical computational complexity analysis, 3) major rewrit