Search

We present a multiscale mixed finite element method for solving second order elliptic equations with general $L^{\infty}$-coefficients arising from flow in highly heterogeneous porous media. Our approach is based on a multiscale spectral generalized finite element method (MS-GFEM) and exploits the superior local mass conservation properties of mixed finite elements. Following the MS-GFEM framework, optimal local approximation spaces are built for the velocity field by solving local eigenvalue problems over generalized harmonic spaces. The resulting global velocity space is then enriched suitably to ensure inf-sup stability. We develop the mixed MS-GFEM for both continuous and discrete formulations, with Raviart-Thomas based mixed finite elements underlying the discrete method. Exponential convergence with respect to local degrees of freedom is proven at both the continuous and discrete levels. Numerical results are presented to support the theory and to validate the proposed method.

Uncertainty Quantification (UQ) is vital to safety-critical model-based analyses, but the widespread adoption of sophisticated UQ methods is limited by technical complexity. In this paper, we introduce UM-Bridge (the UQ and Modeling Bridge), a high-level abstraction and software protocol that facilitates universal interoperability of UQ software with simulation codes. It breaks down the technical complexity of advanced UQ applications and enables separation of concerns between experts. UM-Bridge democratizes UQ by allowing effective interdisciplinary collaboration, accelerating the development of advanced UQ methods, and making it easy to perform UQ analyses from prototype to High Performance Computing (HPC) scale. In addition, we present a library of ready-to-run UQ benchmark problems, all easily accessible through UM-Bridge. These benchmarks support UQ methodology research, enabling reproducible performance comparisons. We demonstrate UM-Bridge with several scientific applications, harnessing HPC resources even using UQ codes not designed with HPC support., Comment: Fix broken reference, add R. Akers as author for contribution to tritium benchmark

We present a flexible method for computing Bayesian optimal experimental designs (BOEDs) for inverse problems with intractable posteriors. The approach is applicable to a wide range of BOED problems and can accommodate various optimality criteria, prior distributions and noise models. The key to our approach is the construction of a transport-map-based surrogate to the joint probability law of the design, observational and inference random variables. This order-preserving transport map is constructed using tensor trains and can be used to efficiently sample from (and evaluate approximate densities of) conditional distributions that are required in the evaluation of many commonly-used optimality criteria. The algorithm is also extended to sequential data acquisition problems, where experiments can be performed in sequence to update the state of knowledge about the unknown parameters. The sequential BOED problem is made computationally feasible by preconditioning the approximation of the joint density at the current stage using transport maps constructed at previous stages. The flexibility of our approach in finding optimal designs is illustrated with some numerical examples inspired by disease modeling and the reconstruction of subsurface structures in aquifers.

We introduce a new Projected Rayleigh Quotient Iteration aimed at improving the convergence behaviour of classic Rayleigh Quotient iteration (RQI) by incorporating approximate information about the target eigenvector at each step. While classic RQI exhibits local cubic convergence for Hermitian matrices, its global behaviour can be unpredictable, whereby it may converge to an eigenvalue far away from the target, even when started with accurate initial conditions. This problem is exacerbated when the eigenvalues are closely spaced. The key idea of the new algorithm is at each step to add a complex-valued projection to the original matrix (that depends on the current eigenvector approximation), such that the unwanted eigenvalues are lifted into the complex plane while the target stays close to the real line, thereby increasing the spacing between the target eigenvalue and the rest of the spectrum. Making better use of the eigenvector approximation leads to more robust convergence behaviour and the new method converges reliably to the correct target eigenpair for a significantly wider range of initial vectors than does classic RQI. We prove that the method converges locally cubically and we present several numerical examples demonstrating the improved global convergence behaviour. In particular, we apply it to compute eigenvalues in a band-gap spectrum of a Sturm-Liouville operator used to model photonic crystal fibres, where the target and unwanted eigenvalues are closely spaced. The examples show that the new method converges to the desired eigenpair even when the eigenvalue spacing is very small, often succeeding when classic RQI fails., Comment: 23 pages, 6 figures, 1 table

At the heart of the Met Office climate and weather forecasting capabilities lies a sophisticated numerical model which solves the equations of large-scale atmospheric flow. Since this model uses semi-implicit time-stepping, it requires the repeated solution of a large sparse system of linear equations with hundreds of millions of unknowns. This is one of the computational bottlenecks of operational forecasts and efficient numerical algorithms are crucial to ensure optimal performance. We developed and implemented a bespoke multigrid solver to address this challenge. Our solver reduces the time for solving the linear system by a factor two, compared to the previously used BiCGStab method. This leads to significant improvements of overall model performance: global forecasts can be produced 10%-15% faster. Multigrid also avoids stagnating convergence of the iterative scheme in single precision. By allowing better utilisation of computational resources, our work has led to estimated annual cost savings of GBP 300k for the Met Office., Comment: 6 pages, 2 figures, 1 table; This is a preprint of the following chapter: A. Malcolm, E. H. Mueller, R. Scheichl: "Improving Met Office Weather and Climate Forecasts with Bespoke Multigrid Solvers", to be published in "More UK Success Stories in Industrial Mathematics", edited by Philip J. Aston, 2024, Springer reproduced with permission of Springer Nature Switzerland AG

Treating uncertainties in models is essential in many fields of science and engineering. Uncertainty quantification (UQ) on complex and computationally costly numerical models necessitates a combination of efficient model solvers, advanced UQ methods and HPC-scale resources. The resulting technical complexities as well as lack of separation of concerns between UQ and model experts is holding back many interesting UQ applications. The aim of this paper is to close the gap between advanced UQ methods and advanced models by removing the hurdle of complex software stack integration, which in turn will offer a straightforward way to scale even prototype-grade UQ applications to high-performance resources. We achieve this goal by introducing a parallel software architecture based on UM-Bridge, a universal interface for linking UQ and models. We present three realistic applications from different areas of science and engineering, scaling from single machines to large clusters on the Google Cloud Platform.

We develop new multilevel Monte Carlo (MLMC) methods to estimate the expectation of the smallest eigenvalue of a stochastic convection-diffusion operator with random coefficients. The MLMC method is based on a sequence of finite element (FE) discretizations of the eigenvalue problem on a hierarchy of increasingly finer meshes. For the discretized, algebraic eigenproblems we use both the Rayleigh quotient (RQ) iteration and implicitly restarted Arnoldi (IRA), providing an analysis of the cost in each case. By studying the variance on each level and adapting classical FE error bounds to the stochastic setting, we are able to bound the total error of our MLMC estimator and provide a complexity analysis. As expected, the complexity bound for our MLMC estimator is superior to plain Monte Carlo. To improve the efficiency of the MLMC further, we exploit the hierarchy of meshes and use coarser approximations as starting values for the eigensolvers on finer ones. To improve the stability of the MLMC method for convection-dominated problems, we employ two additional strategies. First, we consider the streamline upwind Petrov--Galerkin formulation of the discrete eigenvalue problem, which allows us to start the MLMC method on coarser meshes than is possible with standard FEs. Second, we apply a homotopy method to add stability to the eigensolver for each sample. Finally, we present a multilevel quasi-Monte Carlo method that replaces Monte Carlo with a quasi-Monte Carlo (QMC) rule on each level. Due to the faster convergence of QMC, this improves the overall complexity. We provide detailed numerical results comparing our different strategies to demonstrate the practical feasibility of the MLMC method in different use cases. The results support our complexity analysis and further demonstrate the superiority over plain Monte Carlo in all cases.

In this paper, the first large-scale application of multiscale-spectral generalized finite element methods (MS-GFEM) to composite aero-structures is presented. The crucial novelty lies in the introduction of A-harmonicity in the local approximation spaces, which in contrast to [Babuska, Lipton, Multiscale Model. Simul. 9, 2011] is enforced more efficiently via a constraint in the local eigenproblems. This significant modification leads to excellent approximation properties, which turn out to be essential to capture accurately material strains and stresses with a low dimensional approximation space, hence maximising model order reduction. The implementation of the framework in the DUNE software package, as well as a detailed description of all components of the method are presented and exemplified on a composite laminated beam under compressive loading. The excellent parallel scalability of the method, as well as its superior performance compared to the related, previously introduced GenEO method are demonstrated on two realistic application cases, including a C-shaped wing spar with complex geometry. Further, by allowing low-cost approximate solves for closely related models or geometries this efficient, novel technology provides the basis for future applications in optimisation or uncertainty quantification on challenging problems in composite aero-structures.

We propose a deep importance sampling method that is suitable for estimating rare event probabilities in high-dimensional problems. We approximate the optimal importance distribution in a general importance sampling problem as the pushforward of a reference distribution under a composition of order-preserving transformations, in which each transformation is formed by a squared tensor-train decomposition. The squared tensor-train decomposition provides a scalable ansatz for building order-preserving high-dimensional transformations via density approximations. The use of composition of maps moving along a sequence of bridging densities alleviates the difficulty of directly approximating concentrated density functions. To compute expectations over unnormalized probability distributions, we design a ratio estimator that estimates the normalizing constant using a separate importance distribution, again constructed via a composition of transformations in tensor-train format. This offers better theoretical variance reduction compared with self-normalized importance sampling, and thus opens the door to efficient computation of rare event probabilities in Bayesian inference problems. Numerical experiments on problems constrained by differential equations show little to no increase in the computational complexity with the event probability going to zero, and allow to compute hitherto unattainable estimates of rare event probabilities for complex, high-dimensional posterior densities.

In this work we propose an adaptive multilevel version of subset simulation to estimate the probability of rare events for complex physical systems. Given a sequence of nested failure domains of increasing size, the rare event probability is expressed as a product of conditional probabilities. The proposed new estimator uses different model resolutions and varying numbers of samples across the hierarchy of nested failure sets. In order to dramatically reduce the computational cost, we construct the intermediate failure sets such that only a small number of expensive high-resolution model evaluations are needed, whilst the majority of samples can be taken from inexpensive low-resolution simulations. A key idea in our new estimator is the use of a posteriori error estimators combined with a selective mesh refinement strategy to guarantee the critical subset property that may be violated when changing model resolution from one failure set to the next. The efficiency gains and the statistical properties of the estimator are investigated both theoretically via shaking transformations, as well as numerically. On a model problem from subsurface flow, the new multilevel estimator achieves gains of more than a factor 60 over standard subset simulation for a practically relevant relative error of 25%.

Earth's temperature variability can be partitioned into internal and externally-forced components. Yet, underlying mechanisms and their relative contributions remain insufficiently understood, especially on decadal to centennial timescales. Important reasons for this are difficulties in isolating internal and externally-forced variability. Here, we provide a physically-motivated emulation of global mean surface temperature (GMST) variability, which allows for the separation of internal and external variations. To this end, we introduce the ``ClimBayes'' software package, which infers climate parameters from a stochastic energy balance model (EBM) with a Bayesian approach. We apply our method to GMST data from temperature observations and 20 last millennium simulations from climate models of intermediate to high complexity. This yields the best estimates of the EBM's forced and forced + internal response, which we refer to as emulated variability. The timescale-dependent variance is obtained from spectral analysis. In particular, we contrast the emulated forced and forced + internal variance on interannual to centennial timescales with that of the GMST target. Our findings show that a stochastic EBM closely approximates the power spectrum and timescale-dependent variance of GMST as simulated by modern climate models. Small deviations at interannual timescales can be attributed to the simplified representation of internal variability and, in particular, the absence of (pseudo-)oscillatory modes in the stochastic EBM. Altogether, we demonstrate the potential of combining Bayesian inference with conceptual climate models to emulate statistics of climate variables across timescales., Comment: The following article has been published in Chaos: An Interdisciplinary Journal of Nonlinear Science and can be found at https://aip.scitation.org/doi/abs/10.1063/5.0106123

We present a multilevel Monte Carlo simulation method for analysing multi-scale physical systems via a hierarchy of coarse-grained representations, to obtain numerically-exact results, at the most detailed level. We apply the method to a mixture of size-asymmetric hard spheres, in the grand canonical ensemble. A three-level version of the method is compared with a previously-studied two-level version. The extra level interpolates between the full mixture and a coarse-grained description where only the large particles are present -- this is achieved by restricting the small particles to regions close to the large ones. The three-level method improves the performance of the estimator, at fixed computational cost. We analyse the asymptotic variance of the estimator, and discuss the mechanisms for the improved performance., Comment: 24 pages, authors' accepted version

Pattern formation in biological tissues plays an important role in the development of living organisms. Since the classical work of Alan Turing, a pre-eminent way of modelling has been through reaction-diffusion mechanisms. More recently, alternative models have been proposed, that link dynamics of diffusing molecular signals with tissue mechanics. In order to distinguish among different models, they should be compared to experimental observations. However, in many experimental situations only the limiting, stationary regime of the pattern formation process is observable, without knowledge of the transient behaviour or the initial state. The unstable nature of the underlying dynamics in all alternative models seriously complicates model and parameter identification, since small changes in the initial condition lead to distinct stationary patterns. To overcome this problem the initial state of the model can be randomised. In the latter case, fixed values of the model parameters correspond to a family of patterns rather than a fixed stationary solution, and standard approaches to compare pattern data directly with model outputs, e.g., in the least squares sense, are not suitable. Instead, statistical characteristics of the patterns should be compared, which is difficult given the typically limited amount of available data in practical applications. To deal with this problem, we extend a recently developed statistical approach for parameter identification using pattern data, the so-called Correlation Integral Likelihood (CIL) method. We suggest modifications that allow increasing the accuracy of the identification process without resizing the data set. The proposed approach is tested using different classes of pattern formation models. For all considered equations, parallel GPU-based implementations of the numerical solvers with efficient time stepping schemes are provided., Comment: Minor errors corrected, additional clarification added, results unchanged

We develop a novel Markov chain Monte Carlo (MCMC) method that exploits a hierarchy of models of increasing complexity to efficiently generate samples from an unnormalized target distribution. Broadly, the method rewrites the Multilevel MCMC approach of Dodwell et al. (2015) in terms of the Delayed Acceptance (DA) MCMC of Christen & Fox (2005). In particular, DA is extended to use a hierarchy of models of arbitrary depth, and allow subchains of arbitrary length. We show that the algorithm satisfies detailed balance, hence is ergodic for the target distribution. Furthermore, multilevel variance reduction is derived that exploits the multiple levels and subchains, and an adaptive multilevel correction to coarse-level biases is developed. Three numerical examples of Bayesian inverse problems are presented that demonstrate the advantages of these novel methods. The software and examples are available in PyMC3., Comment: 29 pages, 12 figures

In this paper, a generalized finite element method (GFEM) with optimal local approximation spaces for solving high-frequency heterogeneous Helmholtz problems is systematically studied. The local spaces are built from selected eigenvectors of carefully designed local eigenvalue problems defined on generalized harmonic spaces. At both continuous and discrete levels, $(i)$ wavenumber explicit and nearly exponential decay rates for local and global approximation errors are obtained without any assumption on the size of subdomains; $(ii)$ a quasi-optimal convergence of the method is established by assuming that the size of subdomains is $O(1/k)$ ($k$ is the wavenumber). A novel resonance effect between the wavenumber and the dimension of local spaces on the decay of error with respect to the oversampling size is implied by the analysis. Furthermore, for fixed dimensions of local spaces, the discrete local errors are proved to converge as $h\rightarrow 0$ ($h$ denoting the mesh size) towards the continuous local errors. The method at the continuous level extends the plane wave partition of unity method [I. Babuska and J. M. Melenk, Int.\;J.\;Numer.\;Methods Eng., 40 (1997), pp.~727--758] to the heterogeneous-coefficients case, and at the discrete level, it delivers an efficient non-iterative domain decomposition method for solving discrete Helmholtz problems resulting from standard FE discretizations. Numerical results are provided to confirm the theoretical analysis and to validate the proposed method.

GenEO (`Generalised Eigenvalue problems on the Overlap') is a method for computing an operator-dependent spectral coarse space to be combined with local solves on subdomains to form a robust parallel domain decomposition preconditioner for elliptic PDEs. It has previously been proved, in the self-adjoint and positive-definite case, that this method, when used as a preconditioner for conjugate gradients, yields iteration numbers which are completely independent of the heterogeneity of the coefficient field of the partial differential operator. We extend this theory to the case of convection-diffusion-reaction problems, which may be non-self-adjoint and indefinite, and whose discretisations are solved with preconditioned GMRES. The GenEO coarse space is defined here using a generalised eigenvalue problem based on a self-adjoint and positive-definite subproblem. We prove estimates on GMRES iteration counts which are independent of the variation of the coefficient of the diffusion term in the operator and depend only very mildly on variations of the other coefficients. These are proved under the assumption that the subdomain diameter is sufficiently small and the eigenvalue tolerance for building the coarse space is sufficiently large. While the iteration number estimates do grow as the non-self-adjointness and indefiniteness of the operator increases, practical tests indicate the deterioration is much milder. Thus we obtain an iterative solver which is efficient in parallel and very effective for a wide range of convection--diffusion--reaction problems.

Numerical models of complex real-world phenomena often necessitate High Performance Computing (HPC). Uncertainties increase problem dimensionality further and pose even greater challenges. We present a parallelization strategy for multilevel Markov chain Monte Carlo, a state-of-the-art, algorithmically scalable Uncertainty Quantification (UQ) algorithm for Bayesian inverse problems, and a new software framework allowing for large-scale parallelism across forward model evaluations and the UQ algorithms themselves. The main scalability challenge presents itself in the form of strong data dependencies introduced by the MLMCMC method, prohibiting trivial parallelization. Our software is released as part of the modular and open-source MIT UQ Library (MUQ), and can easily be coupled with arbitrary user codes. We demonstrate it using the DUNE and the ExaHyPE Engine. The latter provides a realistic, large-scale tsunami model in which identify the source of a tsunami from buoy-elevation data.

This paper is concerned with error estimates of the fully discrete generalized finite element method (GFEM) with optimal local approximation spaces for solving elliptic problems with heterogeneous coefficients. The local approximation spaces are constructed using eigenvectors of local eigenvalue problems solved by the finite element method on some sufficiently fine mesh with mesh size $h$. The error bound of the discrete GFEM approximation is proved to converge as $h\rightarrow 0$ towards that of the continuous GFEM approximation, which was shown to decay nearly exponentially in previous works. Moreover, even for fixed mesh size $h$, a nearly exponential rate of convergence of the local approximation errors with respect to the dimension of the local spaces is established. An efficient and accurate method for solving the discrete eigenvalue problems is proposed by incorporating the discrete $A$-harmonic constraint directly into the eigensolver. Numerical experiments are carried out to confirm the theoretical results and to demonstrate the effectiveness of the method.

We use a two-level simulation method to analyse the critical point associated with demixing of binary hard sphere mixtures. The method exploits an accurate coarse-grained model with two-body and three-body effective interactions. Using this model within the two-level methodology allows computation of properties of the full (fine-grained) mixture. The critical point is located by computing the probability distribution for the number of large particles in the grand canonical ensemble, and matching to the universal form for the $3d$ Ising universality class. The results have a strong and unexpected dependence on the size ratio between large and small particles, which is related to three-body effective interactions, and the geometry of the underlying hard sphere packings., Comment: 11 pages

Highly heterogeneous, anisotropic coefficients, e.g. in the simulation of carbon-fibre composite components, can lead to extremely challenging finite element systems. Direct solvers for the resulting large and sparse linear systems suffer from severe memory requirements and limited parallel scalability, while iterative solvers in general lack robustness. Two-level spectral domain decomposition methods can provide such robustness for symmetric positive definite linear systems, by using coarse spaces based on independent generalized eigenproblems in the subdomains. Rigorous condition number bounds are independent of mesh size, number of subdomains, as well as coefficient contrast. However, their parallel scalability is still limited by the fact that (in order to guarantee robustness) the coarse problem is solved via a direct method. In this paper, we introduce a multilevel variant in the context of subspace correction methods and provide a general convergence theory for its robust convergence for abstract, elliptic variational problems. Assumptions of the theory are verified for conforming, as well as for discontinuous Galerkin methods applied to a scalar diffusion problem. Numerical results illustrate the performance of the method for two- and three-dimensional problems and for various discretization schemes, in the context of scalar diffusion and linear elasticity., Comment: 25 pages

Motivated by recent work on coarse spaces for Helmholtz problems, we provide in this paper a comparative study on the use of spectral coarse spaces of GenEO type for heterogeneous indefinite elliptic problems within an additive overlapping Schwarz method. In particular, we focus here on two different but related formulations of local generalised eigenvalue problems and compare their performance numerically. Even though their behaviour seems to be very similar for several well-known heterogeneous test cases that are mildly indefinite, only one of the coarse spaces has so far been analysed theoretically, while the other one leads to a significantly more robust domain decomposition method when the indefiniteness is increased. We present a summary of upcoming results developing such a theory and describe how the numerical experiments illustrate it., Comment: 26th International Domain Decomposition Conference (DD26) proceedings

In this paper, the generalized finite element method (GFEM) for solving second order elliptic equations with rough coefficients is studied. New optimal local approximation spaces for GFEMs based on local eigenvalue problems involving a partition of unity are presented. These new spaces have advantages over those proposed in [I. Babuska and R. Lipton, Multiscale Model.\;\,Simul., 9 (2011), pp.~373--406]. First, in addition to a nearly exponential decay rate of the local approximation errors with respect to the dimensions of the local spaces, the rate of convergence with respect to the size of the oversampling region is also established. Second, the theoretical results hold for problems with mixed boundary conditions defined on general Lipschitz domains. Finally, an efficient and easy-to-implement technique for generating the discrete $A$-harmonic spaces is proposed which relies on solving an eigenvalue problem associated with the Dirichlet-to-Neumann operator, leading to a substantial reduction in computational cost. Numerical experiments are presented to support the theoretical analysis and to confirm the effectiveness of the new method.

Stochastic PDE eigenvalue problems often arise in the field of uncertainty quantification, whereby one seeks to quantify the uncertainty in an eigenvalue, or its eigenfunction. In this paper we present an efficient multilevel quasi-Monte Carlo (MLQMC) algorithm for computing the expectation of the smallest eigenvalue of an elliptic eigenvalue problem with stochastic coefficients. Each sample evaluation requires the solution of a PDE eigenvalue problem, and so tackling this problem in practice is notoriously computationally difficult. We speed up the approximation of this expectation in four ways: we use a multilevel variance reduction scheme to spread the work over a hierarchy of FE meshes and truncation dimensions; we use QMC methods to efficiently compute the expectations on each level; we exploit the smoothness in parameter space and reuse the eigenvector from a nearby QMC point to reduce the number of iterations of the eigensolver; and we utilise a two-grid discretisation scheme to obtain the eigenvalue on the fine mesh with a single linear solve. The full error analysis of a basic MLQMC algorithm is given in the companion paper [Gilbert and Scheichl, 2022], and so in this paper we focus on how to further improve the efficiency and provide theoretical justification for using nearby QMC points and two-grid methods. Numerical results are presented that show the efficiency of our algorithm, and also show that the four strategies we employ are complementary.

Uncertainty Quantification through Markov Chain Monte Carlo (MCMC) can be prohibitively expensive for target probability densities with expensive likelihood functions, for instance when the evaluation it involves solving a Partial Differential Equation (PDE), as is the case in a wide range of engineering applications. Multilevel Delayed Acceptance (MLDA) with an Adaptive Error Model (AEM) is a novel approach, which alleviates this problem by exploiting a hierarchy of models, with increasing complexity and cost, and correcting the inexpensive models on-the-fly. The method has been integrated within the open-source probabilistic programming package PyMC3 and is available in the latest development version. In this paper, the algorithm is presented along with an illustrative example., Comment: 8 pages, 4 figures, accepted for Machine Learning for Engineering Modeling, Simulation, and Design Workshop at Neural Information Processing Systems 2020

Stochastic PDE eigenvalue problems are useful models for quantifying the uncertainty in several applications from the physical sciences and engineering, e.g., structural vibration analysis, the criticality of a nuclear reactor or photonic crystal structures. In this paper we present a multilevel quasi-Monte Carlo (MLQMC) method for approximating the expectation of the minimal eigenvalue of an elliptic eigenvalue problem with coefficients that are given as a series expansion of countably-many stochastic parameters. The MLQMC algorithm is based on a hierarchy of discretisations of the spatial domain and truncations of the dimension of the stochastic parameter domain. To approximate the expectations, randomly shifted lattice rules are employed. This paper is primarily dedicated to giving a rigorous analysis of the error of this algorithm. A key step in the error analysis requires bounds on the mixed derivatives of the eigenfunction with respect to both the stochastic and spatial variables simultaneously. Under stronger smoothness assumptions on the parametric dependence, our analysis also extends to multilevel higher-order quasi-Monte Carlo rules. An accompanying paper [Gilbert and Scheichl, 2022], focusses on practical extensions of the MLQMC algorithm to improve efficiency, and presents numerical results.

Monte Carlo simulations of quantum field theories on a lattice become increasingly expensive as the continuum limit is approached since the cost per independent sample grows with a high power of the inverse lattice spacing. Simulations on fine lattices suffer from critical slowdown, the rapid growth of autocorrelations in the Markov chain. This causes a strong increase in the number of lattice configurations that have to be generated to obtain statistically significant results. This paper discusses hierarchical sampling methods to tame the growth in autocorrelations. Combined with multilevel variance reduction, this significantly reduces the computational cost of simulations for given tolerances $\epsilon_{\text{disc}}$ on the discretisation error and $\epsilon_{\text{stat}}$ on the statistical error. For observables with lattice errors of order $\alpha$ and integrated autocorrelation times that grow like $\tau_{\mathrm{int}}\propto a^{-z}$, multilevel Monte Carlo (MLMC) reduces the cost from $\mathcal{O}(\epsilon_{\text{stat}}^{-2}\epsilon_{\text{disc}}^{-(1+z)/\alpha})$ to $\mathcal{O}(\epsilon_{\text{stat}}^{-2}\vert\log \epsilon_{\text{disc}} \vert^2+\epsilon_{\text{disc}}^{-1/\alpha})$ or $\mathcal{O}(\epsilon_{\text{stat}}^{-2}+\epsilon_{\text{disc}}^{-1/\alpha})$. Higher gains are expected for simulations of quantum field theories in $D$ dimensions. The efficiency of the approach is demonstrated on two model systems, including a topological oscillator that is badly affected by critical slowdown from topological charge freezing. On fine lattices, the new methods are orders of magnitude faster than standard Hybrid Monte Carlo sampling. For high resolutions, MLMC can be used to accelerate even the cluster algorithm for the topological oscillator. Performance is further improved through perturbative matching which guarantees efficient coupling of theories on the multilevel hierarchy., Comment: 23 pages, 17 figures, 5 tables

Low-rank tensor approximations have shown great potential for uncertainty quantification in high dimensions, for example, to build surrogate models that can be used to speed up large-scale inference problems (Eigel et al., Inverse Problems 34, 2018; Dolgov et al., Statistics & Computing 30, 2020). The feasibility and efficiency of such approaches depends critically on the rank that is necessary to represent or approximate the underlying distribution. In this paper, a-priori rank bounds for approximations in the functional tensor-train representation for the case of Gaussian models are developed. It is shown that under suitable conditions on the precision matrix, the Gaussian density can be approximated to high accuracy without suffering from an exponential growth of complexity as the dimension increases. These results provide a rigorous justification of the suitability and the limitations of low-rank tensor methods in a simple but important model case. Numerical experiments confirm that the rank bounds capture the qualitative behavior of the rank structure when varying the parameters of the precision matrix and the accuracy of the approximation. Finally, the practical relevance of the theoretical results is demonstrated in the context of a Bayesian filtering problem., Comment: 25 pages, 6 figures

We present a non-trivial integration of dimension-independent likelihood-informed (DILI) MCMC (Cui, Law, Marzouk, 2016) and the multilevel MCMC (Dodwell et al., 2015) to explore the hierarchy of posterior distributions. This integration offers several advantages: First, DILI-MCMC employs an intrinsic likelihood-informed subspace (LIS) (Cui et al., 2014) -- which involves a number of forward and adjoint model simulations -- to design accelerated operator-weighted proposals. By exploiting the multilevel structure of the discretised parameters and discretised forward models, we design a Rayleigh-Ritz procedure to significantly reduce the computational effort in building the LIS and operating with DILI proposals. Second, the resulting DILI-MCMC can drastically improve the sampling efficiency of MCMC at each level, and hence reduce the integration error of the multilevel algorithm for fixed CPU time. To be able to fully exploit the power of multilevel MCMC and to reduce the dependencies of samples on different levels for a parallel implementation, we also suggest a new pooling strategy for allocating computational resources across different levels and constructing Markov chains at higher levels conditioned on those simulated on lower levels. Numerical results confirm the improved computational efficiency of the multilevel DILI approach.

We present a method that exploits self-consistent simulation of coarse-grained and fine-grained models, in order to analyse properties of physical systems. The method uses the coarse-grained model to obtain a first estimate of the quantity of interest, before computing a correction by analysing properties of the fine system. We illustrate the method by applying it to the Asakura-Oosawa (AO) model of colloid-polymer mixtures. We show that the liquid-vapour critical point in that system is affected by three-body interactions which are neglected in the corresponding coarse-grained model. We analyse the size of this effect and the nature of the three-body interactions. We also analyse the accuracy of the method, as a function of the associated computational effort., Comment: 17 pages

We present an efficient implementation of the highly robust and scalable GenEO preconditioner in the high-performance PDE framework DUNE. The GenEO coarse space is constructed by combining low energy solutions of a local generalised eigenproblem using a partition of unity. In this paper we demonstrate both weak and strong scaling for the GenEO solver on over 15,000 cores by solving an industrially motivated problem with over 200 million degrees of freedom. Further, we show that for highly complex parameter distributions arising in certain real-world applications, established methods become intractable while GenEO remains fully effective. The purpose of this paper is two-fold: to demonstrate the robustness and high parallel efficiency of the solver and to document the technical details that are crucial to the efficiency of the code., Comment: Fixed error in definition of partition of unity

The periodization of a stationary Gaussian random field on a sufficiently large torus comprising the spatial domain of interest is the basis of various efficient computational methods, such as the classical circulant embedding technique using the fast Fourier transform for generating samples on uniform grids. For the family of Mat\'ern covariances with smoothness index $\nu$ and correlation length $\lambda$, we analyse the nonsmooth periodization (corresponding to classical circulant embedding) and an alternative procedure using a smooth truncation of the covariance function. We solve two open problems: the first concerning the $\nu$-dependent asymptotic decay of eigenvalues of the resulting circulant in the nonsmooth case, the second concerning the required size in terms of $\nu$, $\lambda$ of the torus when using a smooth periodization. In doing this we arrive at a complete characterisation of the performance of these two approaches. Both our theoretical estimates and the numerical tests provided here show substantial advantages of smooth truncation., Comment: 24 pages, 3 figures

Many recent invertible neural architectures are based on coupling block designs where variables are divided in two subsets which serve as inputs of an easily invertible (usually affine) triangular transformation. While such a transformation is invertible, its Jacobian is very sparse and thus may lack expressiveness. This work presents a simple remedy by noting that subdivision and (affine) coupling can be repeated recursively within the resulting subsets, leading to an efficiently invertible block with dense, triangular Jacobian. By formulating our recursive coupling scheme via a hierarchical architecture, HINT allows sampling from a joint distribution p(y,x) and the corresponding posterior p(x|y) using a single invertible network. We evaluate our method on some standard data sets and benchmark its full power for density estimation and Bayesian inference on a novel data set of 2D shapes in Fourier parameterization, which enables consistent visualization of samples for different dimensionalities., Comment: Published at AAAI 2021

We present an analysis of multilevel Monte Carlo techniques for the forward problem of uncertainty quantification for the radiative transport equation, when the coefficients ({\em cross-sections}) are heterogenous random fields. To do this, we first give a new error analysis for the combined spatial and angular discretisation in the deterministic case, with error estimates which are explicit in the coefficients (and allow for very low regularity and jumps). This detailed error analysis is done for the 1D space - 1D angle slab geometry case with classical diamond differencing. Under reasonable assumptions on the statistics of the coefficients, we then prove an error estimate for the random problem in a suitable Bochner space. Because the problem is not self-adjoint, stability can only be proved under a path-dependent mesh resolution condition. This means that, while the Bochner space error estimate is of order $\mathcal{O}(h^\eta)$ for some $\eta$, where $h$ is a (deterministically chosen) mesh diameter, smaller mesh sizes might be needed for some realisations. Under reasonable assumptions we show that the expected cost for computing a typical quantity of interest remains of the same order as for a single sample. This leads to rigorous complexity estimates for Monte Carlo and multilevel Monte Carlo: For particular linear solvers, the multilevel version gives up to two orders of magnitude improvement over Monte Carlo. We provide numerical results supporting the theory.

We propose and analyse a fully adaptive strategy for solving elliptic PDEs with random data in this work. A hierarchical sequence of adaptive mesh refinements for the spatial approximation is combined with adaptive anisotropic sparse Smolyak grids in the stochastic space in such a way as to minimize the computational cost. The novel aspect of our strategy is that the hierarchy of spatial approximations is sample dependent so that the computational effort at each collocation point can be optimised individually. We outline a rigorous analysis for the convergence and computational complexity of the adaptive multilevel algorithm and we provide optimal choices for error tolerances at each level. Two numerical examples demonstrate the reliability of the error control and the significant decrease in the complexity that arises when compared to single level algorithms and multilevel algorithms that employ adaptivity solely in the spatial discretisation or in the collocation procedure., Comment: 26 pages, 7 figures

A key quantity that occurs in the error analysis of several numerical methods for eigenvalue problems is the distance between the eigenvalue of interest and the next nearest eigenvalue. When we are interested in the smallest or fundamental eigenvalue, we call this the spectral or fundamental gap. In a recent manuscript [Gilbert et al., arXiv:1808.02639], the current authors, together with Frances Kuo, studied an elliptic eigenvalue problem with homogeneous Dirichlet boundary conditions, and with coefficients that depend on an infinite number of uniformly distributed stochastic parameters. In this setting, the eigenvalues, and in turn the eigenvalue gap, also depend on the stochastic parameters. Hence, for a robust error analysis one needs to be able to bound the gap over all possible realisations of the parameters, and because the gap depends on infinitely-many random parameters, this is not trivial. This short note presents, in a simplified setting, an important result that was shown in the paper above. Namely, that, under certain decay assumptions on the coefficient, the spectral gap of such a random elliptic eigenvalue problem can be bounded away from 0, uniformly over the entire infinite-dimensional parameter space., Comment: 12 pages, 5 figures. arXiv admin note: text overlap with arXiv:1808.02639

The key innovation in this paper is an open-source, high-performance iterative solver for high contrast, strongly anisotropic elliptic partial differential equations implemented within dune-pdelab. The iterative solver exploits a robust, scalable two-level additive Schwarz preconditioner, GenEO (Spillane et al. 2014). The development of this solver has been motivated by the need to overcome the limitations of commercially available modeling tools for solving structural analysis simulations in aerospace composite applications. Our software toolbox dune-composites encapsulates the mathematical complexities of the underlying packages within an efficient C++ framework, providing an application interface to our new high-performance solver. We illustrate its use on a range of industrially motivated examples, which should enable other scientists to build on and extend dune-composites and the GenEO preconditioner for use in their own applications. We demonstrate the scalability of the solver on more than 15,000 cores of the UK national supercomputer Archer, solving an aerospace composite problem with over 200 million degrees of freedom in a few minutes. This scale of computation brings composites problems that would otherwise be unthinkable into the feasible range. To demonstrate the wider applicability of the new solver, we also confirm the robustness and scalability of the solver on SPE10, a challenging benchmark in subsurface flow/reservoir simulation.

General multivariate distributions are notoriously expensive to sample from, particularly the high-dimensional posterior distributions in PDE-constrained inverse problems. This paper develops a sampler for arbitrary continuous multivariate distributions that is based on low-rank surrogates in the tensor-train format. We construct a tensor-train approximation to the target probability density function using the cross interpolation, which requires a small number of function evaluations. For sufficiently smooth distributions the storage required for the TT approximation is moderate, scaling linearly with dimension. The structure of the tensor-train surrogate allows efficient sampling by the conditional distribution method. Unbiased estimates may be calculated by correcting the transformed random seeds using a Metropolis--Hastings accept/reject step. Moreover, one can use a more efficient quasi-Monte Carlo quadrature that may be corrected either by a control-variate strategy, or by importance weighting. We show that the error in the tensor-train approximation propagates linearly into the Metropolis--Hastings rejection rate and the integrated autocorrelation time of the resulting Markov chain. These methods are demonstrated in three computed examples: fitting failure time of shock absorbers; a PDE-constrained inverse diffusion problem; and sampling from the Rosenbrock distribution. The delayed rejection adaptive Metropolis (DRAM) algorithm is used as a benchmark. We find that the importance-weight corrected quasi-Monte Carlo quadrature performs best in all computed examples, and is orders-of-magnitude more efficient than DRAM across a wide range of approximation accuracies and sample sizes. Indeed, all the methods developed here significantly outperform DRAM in all computed examples., Comment: 32 pages

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+\delta})$, where $\delta>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.

Stein variational gradient descent (SVGD) was recently proposed as a general purpose nonparametric variational inference algorithm [Liu & Wang, NIPS 2016]: it minimizes the Kullback-Leibler divergence between the target distribution and its approximation by implementing a form of functional gradient descent on a reproducing kernel Hilbert space. In this paper, we accelerate and generalize the SVGD algorithm by including second-order information, thereby approximating a Newton-like iteration in function space. We also show how second-order information can lead to more effective choices of kernel. We observe significant computational gains over the original SVGD algorithm in multiple test cases., Comment: 18 pages, 7 figures

We consider the approximate solution of parametric PDEs using the low-rank Tensor Train (TT) decomposition. Such parametric PDEs arise for example in uncertainty quantification problems in engineering applications. We propose an algorithm that is a hybrid of the alternating least squares and the TT cross methods. It computes a TT approximation of the whole solution, which is beneficial when multiple quantities of interest are sought. This might be needed, for example, for the computation of the probability density function (PDF) via the maximum entropy method [Kavehrad and Joseph, IEEE Trans. Comm., 1986]. The new algorithm exploits and preserves the block diagonal structure of the discretized operator in stochastic collocation schemes. This disentangles computations of the spatial and parametric degrees of freedom in the TT representation. In particular, it only requires solving independent PDEs at a few parameter values, thus allowing the use of existing high performance PDE solvers. In our numerical experiments, we apply the new algorithm to the stochastic diffusion equation and compare it with preconditioned steepest descent in the TT format, as well as with (multilevel) quasi-Monte Carlo and dimension-adaptive sparse grids methods. For sufficiently smooth random fields the new approach is orders of magnitude faster., Comment: 29 pages

Finite element (FE) analysis has the potential to offset much of the expensive experimental testing currently required to certify aerospace laminates. However, large numbers of degrees of freedom are necessary to model entire aircraft components whilst accurately resolving micro-scale defects. The new module dune-composites, implemented within DUNE by the authors, provides a tool to efficiently solve large-scale problems using novel iterative solvers. The key innovation is a preconditioner that guarantees a constant number of iterations regardless of the problem size. Its robustness has been shown rigorously in Spillane et al. (Numer. Math. 126, 2014) for isotropic problems. For anisotropic problems in composites it is verified numerically for the first time in this paper. The parallel implementation in DUNE scales almost optimally over thousands of cores. To demonstrate this, we present an original numerical study, varying the shape of a localised wrinkle and the effect this has on the strength of a curved laminate. This requires a high-fidelity mesh containing at least four layers of quadratic elements across each ply and interface layer, underlining the need for dune-composites, which can achieve run times of just over 2 minutes on 2048 cores for realistic composites problems with 173 million degrees of freedom.

We describe modern variants of Monte Carlo methods for Uncertainty Quantification (UQ) of the Neutron Transport Equation, when it is approximated by the discrete ordinates method with diamond differencing. We focus on the mono-energetic 1D slab geometry problem, with isotropic scattering, where the cross-sections are log-normal correlated random fields of possibly low regularity. The paper includes an outline of novel theoretical results on the convergence of the discrete scheme, in the cases of both spatially variable and random cross-sections. We also describe the theory and practice of algorithms for quantifying the uncertainty of a linear functional of the scalar flux, using Monte Carlo and quasi-Monte Carlo methods, and their multilevel variants. A hybrid iterative/direct solver for computing each realisation of the functional is also presented. Numerical experiments show the effectiveness of the hybrid solver and the gains that are possible through quasi-Monte Carlo sampling and multilevel variance reduction. For the multilevel quasi-Monte Carlo method, we observe gains in the computational $\varepsilon$-cost of up to 2 orders of magnitude over the standard Monte Carlo method, and we explain this theoretically. Experiments on problems with up to several thousand stochastic dimensions are included.

A common way to simulate the transport and spread of pollutants in the atmosphere is via stochastic Lagrangian dispersion models. Mathematically, these models describe turbulent transport processes with stochastic differential equations (SDEs). The computational bottleneck is the Monte Carlo algorithm, which simulates the motion of a large number of model particles in a turbulent velocity field; for each particle, a trajectory is calculated with a numerical timestepping method. Choosing an efficient numerical method is particularly important in operational emergency-response applications, such as tracking radioactive clouds from nuclear accidents or predicting the impact of volcanic ash clouds on international aviation, where accurate and timely predictions are essential. In this paper, we investigate the application of the Multilevel Monte Carlo (MLMC) method to simulate the propagation of particles in a representative one-dimensional dispersion scenario in the atmospheric boundary layer. MLMC can be shown to result in asymptotically superior computational complexity and reduced computational cost when compared to the Standard Monte Carlo (StMC) method, which is currently used in atmospheric dispersion modelling. To reduce the absolute cost of the method also in the non-asymptotic regime, it is equally important to choose the best possible numerical timestepping method on each level. To investigate this, we also compare the standard symplectic Euler method, which is used in many operational models, with two improved timestepping algorithms based on SDE splitting methods., Comment: 30 pages, 13 figures, 3 tables, revised version, submitted to Journal of Computational Physics

The computational complexity of naive, sampling-based uncertainty quantification for 3D partial differential equations is extremely high. Multilevel approaches, such as multilevel Monte Carlo (MLMC), can reduce the complexity significantly, but to exploit them fully in a parallel environment, sophisticated scheduling strategies are needed. Often fast algorithms that are executed in parallel are essential to compute fine level samples in 3D, whereas to compute individual coarse level samples only moderate numbers of processors can be employed efficiently. We make use of multiple instances of a parallel multigrid solver combined with advanced load balancing techniques. In particular, we optimize the concurrent execution across the three layers of the MLMC method: parallelization across levels, across samples, and across the spatial grid. The overall efficiency and performance of these methods will be analyzed. Here the scalability window of the multigrid solver is revealed as being essential, i.e., the property that the solution can be computed with a range of process numbers while maintaining good parallel efficiency. We evaluate the new scheduling strategies in a series of numerical tests, and conclude the paper demonstrating large 3D scaling experiments.