Raul Tempone

Joakim Beck, Ben Mansour Dia, Luis FR Espath et al.

In calculating expected information gain in optimal Bayesian experimental design, the computation of the inner loop in the classical double-loop Monte Carlo requires a large number of samples and suffers from underflow if the number of samples is small. These drawbacks can be avoided by using an importance sampling approach. We present a computationally efficient method for optimal Bayesian experimental design that introduces importance sampling based on the Laplace method to the inner loop. We derive the optimal values for the method parameters in which the average computational cost is minimized according to the desired error tolerance. We use three numerical examples to demonstrate the computational efficiency of our method compared with the classical double-loop Monte Carlo, and a more recent single-loop Monte Carlo method that uses the Laplace method as an approximation of the return value of the inner loop. The first example is a scalar problem that is linear in the uncertain parameter. The second example is a nonlinear scalar problem. The third example deals with the optimal sensor placement for an electrical impedance tomography experiment to recover the fiber orientation in laminate composites.

Håkon Hoel, Kody J. H. Law, Raul Tempone

This work embeds a multilevel Monte Carlo sampling strategy into the Monte Carlo step of the ensemble Kalman filter (EnKF) in the setting of finite dimensional signal evolution and noisy discrete-time observations. The signal dynamics is assumed to be governed by a stochastic differential equation (SDE), and a hierarchy of time grids is introduced for multilevel numerical integration of that SDE. The resulting multilevel EnKF is proved to asymptotically outperform EnKF in terms of computational cost versus approximation accuracy. The theoretical results are illustrated numerically.

Kody J. H. Law, Hamidou Tembine, Raul Tempone

The proof of convergence of the standard ensemble Kalman filter (EnKF) from Legland etal. (2011) is extended to non-Gaussian state space models. A density-based deterministic approximation of the mean-field limit EnKF (DMFEnKF) is proposed, consisting of a PDE solver and a quadrature rule. Given a certain minimal order of convergence $κ$ between the two, this extends to the deterministic filter approximation, which is therefore asymptotically superior to standard EnKF when the dimension $d<2κ$. The fidelity of approximation of the true distribution is also established using an extension of total variation metric to random measures. This is limited by a Gaussian bias term arising from non-linearity/non-Gaussianity of the model, which exists for both DMFEnKF and standard EnKF. Numerical results support and extend the theory.

Quan Long, Mohammad Motamed, Raul Tempone

We develop a fast method for optimally designing experiments in the context of statistical seismic source inversion. In particular, we efficiently compute the optimal number and locations of the receivers or seismographs. The seismic source is modeled by a point moment tensor multiplied by a time-dependent function. The parameters include the source location, moment tensor components, and start time and frequency in the time function. The forward problem is modeled by elastodynamic wave equations. We show that the Hessian of the cost functional, which is usually defined as the square of the weighted L2 norm of the difference between the experimental data and the simulated data, is proportional to the measurement time and the number of receivers. Consequently, the posterior distribution of the parameters, in a Bayesian setting, concentrates around the "true" parameters, and we can employ Laplace approximation and speed up the estimation of the expected Kullback-Leibler divergence (expected information gain), the optimality criterion in the experimental design procedure. Since the source parameters span several magnitudes, we use a scaling matrix for efficient control of the condition number of the original Hessian matrix. We use a second-order accurate finite difference method to compute the Hessian matrix and either sparse quadrature or Monte Carlo sampling to carry out numerical integration. We demonstrate the efficiency, accuracy, and applicability of our method on a two-dimensional seismic source inversion problem.

Abdul-Lateef Haji-Ali, Raul Tempone

We address the approximation of functionals depending on a system of particles, described by stochastic differential equations (SDEs), in the mean-field limit when the number of particles approaches infinity. This problem is equivalent to estimating the weak solution of the limiting McKean-Vlasov SDE. To that end, our approach uses systems with finite numbers of particles and a time-stepping scheme. In this case, there are two discretization parameters: the number of time steps and the number of particles. Based on these two parameters, we consider different variants of the Monte Carlo and Multilevel Monte Carlo (MLMC) methods and show that, in the best case, the optimal work complexity of MLMC, to estimate the functional in one typical setting with an error tolerance of $\mathrm{TOL}$, is $\mathcal O\left({\mathrm{TOL}^{-3}}\right)$ when using the partitioning estimator and the Milstein time-stepping scheme. We also consider a method that uses the recent Multi-index Monte Carlo method and show an improved work complexity in the same typical setting of $\mathcal O\left(\mathrm{TOL}^{-2}\log(\mathrm{TOL}^{-1})^2\right)$. Our numerical experiments are carried out on the so-called Kuramoto model, a system of coupled oscillators.

Christian Bayer, Markus Siebenmorgen, Raul Tempone

We consider the problem of pricing basket options in a multivariate Black Scholes or Variance Gamma model. From a numerical point of view, pricing such options corresponds to moderate and high dimensional numerical integration problems with non-smooth integrands. Due to this lack of regularity, higher order numerical integration techniques may not be directly available, requiring the use of methods like Monte Carlo specifically designed to work for non-regular problems. We propose to use the inherent smoothing property of the density of the underlying in the above models to mollify the payoff function by means of an exact conditional expectation. The resulting conditional expectation is unbiased and yields a smooth integrand, which is amenable to the efficient use of adaptive sparse grid cubature. Numerical examples indicate that the high-order method may perform orders of magnitude faster compared to Monte Carlo or Quasi Monte Carlo in dimensions up to 35.

Chiheb Ben Hammouda, Alvaro Moraes, Raul Tempone

In biochemically reactive systems with small copy numbers of one or more reactant molecules, the dynamics is dominated by stochastic effects. To approximate those systems, discrete state-space and stochastic simulation approaches have been shown to be more relevant than continuous state-space and deterministic ones. In systems characterized by having simultaneously fast and slow timescales, existing discrete space-state stochastic path simulation methods, such as the stochastic simulation algorithm (SSA) and the explicit tau-leap (Explicit-TL) method, can be very slow. Implicit approximations have been developed to improve numerical stability and provide efficient simulation algorithms for those systems. Here, we propose an efficient Multilevel Monte Carlo (MLMC) method in the spirit of the work by Anderson and Higham (2012). This method uses split-step implicit tau-leap (SSI-TL) at levels where the SSI-TL method is not applicable due to numerical stability issues. We present numerical examples that illustrate the performance of the proposed method.

Jesper Karlsson, Raul Tempone

This work develops novel error expansions with computable leading order terms for the global weak error in the tau-leap discretization of pure jump processes arising in kinetic Monte Carlo models. Accurate computable a posteriori error approximations are the basis for adaptive algorithms; a fundamental tool for numerical simulation of both deterministic and stochastic dynamical systems. These pure jump processes are simulated either by the tau-leap method, or by exact simulation, also referred to as dynamic Monte Carlo, the Gillespie algorithm or the Stochastic simulation algorithm. Two types of estimates are presented: an a priori estimate for the relative error that gives a comparison between the work for the two methods depending on the propensity regime, and an a posteriori estimate with computable leading order term.

Raul Tempone, Soeren Wolfers

We provide a general discussion of Smolyak's algorithm for the acceleration of scientific computations. The algorithm first appeared in Smolyak's work on multidimensional integration and interpolation. Since then, it has been generalized in multiple directions and has been associated with the keywords: sparse grids, hyperbolic cross approximation, combination technique, and multilevel methods. Variants of Smolyak's algorithm have been employed in the computation of high-dimensional integrals in finance, chemistry, and physics, in the numerical solution of partial and stochastic differential equations, and in uncertainty quantification. Motivated by this broad and ever-increasing range of applications, we describe a general framework that summarizes fundamental results and assumptions in a concise application-independent manner.

Fabio Nobile, Raul Tempone, Soeren Wolfers

We provide a framework for the sparse approximation of multilinear problems and show that several problems in uncertainty quantification fit within this framework. In these problems, the value of a multilinear map has to be approximated using approximations of different accuracy and computational work of the arguments of this map. We propose and analyze a generalized version of Smolyak's algorithm, which provides sparse approximation formulas with convergence rates that mitigate the curse of dimension that appears in multilinear approximation problems with a large number of arguments. We apply the general framework to response surface approximation and optimization under uncertainty for parametric partial differential equations using kernel-based approximation. The theoretical results are supplemented by numerical experiments.

Thomas Björk, Anders Szepessy, Raul Tempone et al.

We present Monte Carlo-Euler methods for a weak approximation problem related to the Heath-Jarrow-Morton (HJM) term structure model, based on \Ito stochastic differential equations in infinite dimensional spaces, and prove strong and weak error convergence estimates. The weak error estimates are based on stochastic flows and discrete dual backward problems, and they can be used to identify different error contributions arising from time and maturity discretization as well as the classical statistical error due to finite sampling. Explicit formulas for efficient computation of sharp error approximation are included. Due to the structure of the HJM models considered here, the computational effort devoted to the error estimates is low compared to the work to compute Monte Carlo solutions to the HJM model. Numerical examples with known exact solution are included in order to show the behavior of the estimates.

Luis Espath, Pouria Behnoudfar, Raul Tempone

In this work, we further develop the Physics-informed Spectral Learning (PiSL) by Espath et al. \cite{Esp21} based on a discrete $L^2$ projection to solve the discrete Hodge--Helmholtz decomposition from sparse data. Within this physics-informed statistical learning framework, we adaptively build a sparse set of Fourier basis functions with corresponding coefficients by solving a sequence of minimization problems where the set of basis functions is augmented greedily at each optimization problem. Moreover, our PiSL computational framework enjoys spectral (exponential) convergence. We regularize the minimization problems with the seminorm of the fractional Sobolev space in a Tikhonov fashion. In the Fourier setting, the divergence- and curl-free constraints become a finite set of linear algebraic equations. The proposed computational framework combines supervised and unsupervised learning techniques in that we use data concomitantly with the projection onto divergence- and curl-free spaces. We assess the capabilities of our method in various numerical examples including the `Storm of the Century' with satellite data from 1993.

Owen Davis, Mohammad Motamed, Raul Tempone

In this work, we consider the general problem of constructing a neural network surrogate model using multi-fidelity information. Motivated by error-complexity estimates for ReLU neural networks, we formulate the correlation between an inexpensive low-fidelity model and an expensive high-fidelity model as a possibly non-linear residual function. This function defines a mapping between 1) the shared input space of the models along with the low-fidelity model output, and 2) the discrepancy between the outputs of the two models. The computational framework proceeds by training two neural networks to work in concert. The first network learns the residual function on a small set of high- and low-fidelity data. Once trained, this network is used to generate additional synthetic high-fidelity data, which is used in the training of the second network. The trained second network then acts as our surrogate for the high-fidelity quantity of interest. We present four numerical examples to demonstrate the power of the proposed framework, showing that significant savings in computational cost may be achieved when the output predictions are desired to be accurate within small tolerances.

Eike Cramer, Alexander Mitsos, Raul Tempone et al.

Neural networks-based learning of the distribution of non-dispatchable renewable electricity generation from sources such as photovoltaics (PV) and wind as well as load demands has recently gained attention. Normalizing flow density models are particularly well suited for this task due to the training through direct log-likelihood maximization. However, research from the field of image generation has shown that standard normalizing flows can only learn smeared-out versions of manifold distributions. Previous works on normalizing flow-based scenario generation do not address this issue, and the smeared-out distributions result in the sampling of noisy time series. In this paper, we exploit the isometry of the principal component analysis (PCA), which sets up the normalizing flow in a lower-dimensional space while maintaining the direct and computationally efficient likelihood maximization. We train the resulting principal component flow (PCF) on data of PV and wind power generation as well as load demand in Germany in the years 2013 to 2015. The results of this investigation show that the PCF preserves critical features of the original distributions, such as the probability density and frequency behavior of the time series. The application of the PCF is, however, not limited to renewable power generation but rather extends to any data set, time series, or otherwise, which can be efficiently reduced using PCA.

Alexander Litvinenko, Dmitry Logashenko, Raul Tempone et al.

Accurate modeling of contamination in subsurface flow and water aquifers is crucial for agriculture and environmental protection. Here, we demonstrate a parallel method to quantify the propagation of the uncertainty in the dispersal of pollution in subsurface flow. Specifically, we consider the density-driven flow and estimate how uncertainty from permeability and porosity propagates to the solution. We take an Elder-like problem as a numerical benchmark and we use random fields to model the limited knowledge on the porosity and permeability. We construct a low-cost generalized polynomial chaos expansion (gPC) surrogate model, where the gPC coefficients are computed by projection on sparse and full tensor grids. We parallelize both the numerical solver for the deterministic problem based on the multigrid method, and the quadrature over the parametric space

Alexey Chernov, Haakon Hoel, Kody Law et al.

This work embeds a multilevel Monte Carlo (MLMC) sampling strategy into the Monte Carlo step of the ensemble Kalman filter (EnKF), thereby yielding a multilevel ensemble Kalman filter (MLEnKF) which has provably superior asymptotic cost to a given accuracy level. The development of MLEnKF for finite-dimensional state-spaces in the work [20] is here extended to models with infinite-dimensional state- spaces in the form of spatial fields. A concrete example is given to illustrate the results.

Abdul-Lateef Haji-Ali, Fabio Nobile, Lorenzo Tamellini et al.

In this work we introduce the Multi-Index Stochastic Collocation method (MISC) for computing statistics of the solution of a PDE with random data. MISC is a combination technique based on mixed differences of spatial approximations and quadratures over the space of random data. We propose an optimization procedure to select the most effective mixed differences to include in the MISC estimator: such optimization is a crucial step and allows us to build a method that, provided with sufficient solution regularity, is potentially more effective than other multi-level collocation methods already available in literature. We then provide a complexity analysis that assumes decay rates of product type for such mixed differences, showing that in the optimal case the convergence rate of MISC is only dictated by the convergence of the deterministic solver applied to a one dimensional problem. We show the effectiveness of MISC with some computational tests, comparing it with other related methods available in the literature, such as the Multi-Index and Multilevel Monte Carlo, Multilevel Stochastic Collocation, Quasi Optimal Stochastic Collocation and Sparse Composite Collocation methods.

Abdul-Lateef Haji-Ali, Fabio Nobile, Lorenzo Tamellini et al.

We analyze the recent Multi-index Stochastic Collocation (MISC) method for computing statistics of the solution of a partial differential equation (PDEs) with random data, where the random coefficient is parametrized by means of a countable sequence of terms in a suitable expansion. MISC is a combination technique based on mixed differences of spatial approximations and quadratures over the space of random data and, naturally, the error analysis uses the joint regularity of the solution with respect to both the variables in the physical domain and parametric variables. In MISC, the number of problem solutions performed at each discretization level is not determined by balancing the spatial and stochastic components of the error, but rather by suitably extending the knapsack-problem approach employed in the construction of the quasi-optimal sparse-grids and Multi-index Monte Carlo methods. We use a greedy optimization procedure to select the most effective mixed differences to include in the MISC estimator. We apply our theoretical estimates to a linear elliptic PDEs in which the log-diffusion coefficient is modeled as a random field, with a covariance similar to a Matérn model, whose realizations have spatial regularity determined by a scalar parameter. We conduct a complexity analysis based on a summability argument showing algebraic rates of convergence with respect to the overall computational work. The rate of convergence depends on the smoothness parameter, the physical dimensionality and the efficiency of the linear solver. Numerical experiments show the effectiveness of MISC in this infinite-dimensional setting compared with the Multi-index Monte Carlo method and compare the convergence rate against the rates predicted in our theoretical analysis.

Gabriela Malenova, Mohammad Motamed, Olof Runborg et al.

We consider the wave equation with highly oscillatory initial data, where there is uncertainty in the wave speed, initial phase and/or initial amplitude. To estimate quantities of interest related to the solution and their statistics, we combine a high-frequency method based on Gaussian beams with sparse stochastic collocation. Although the wave solution, $u^\varepsilon$, is highly oscillatory in both physical and stochastic spaces, we provide theoretical arguments and numerical evidence that quantities of interest based on local averages of $|u^\varepsilon|^2$ are smooth, with derivatives in the stochastic space uniformly bounded in $\varepsilon$, where $\varepsilon$ denotes the short wavelength. This observable related regularity makes the sparse stochastic collocation approach more efficient than Monte Carlo methods. We present numerical tests that demonstrate this advantage.

Christian Bayer, Alvaro Moraes, Raul Tempone et al.

In this work, we present an extension to the context of Stochastic Reaction Networks (SRNs) of the forward-reverse representation introduced in "Simulation of forward-reverse stochastic representations for conditional diffusions", a 2014 paper by Bayer and Schoenmakers. We apply this stochastic representation in the computation of efficient approximations of expected values of functionals of SNR bridges, i.e., SRNs conditioned to its values in the extremes of given time-intervals. We then employ this SNR bridge-generation technique to the statistical inference problem of approximating the reaction propensities based on discretely observed data. To this end, we introduce a two-phase iterative inference method in which, during phase I, we solve a set of deterministic optimization problems where the SRNs are replaced by their reaction-rate Ordinary Differential Equations (ODEs) approximation; then, during phase II, we apply the Monte Carlo version of the Expectation-Maximization (EM) algorithm starting from the phase I output. By selecting a set of over dispersed seeds as initial points for phase I, the output of parallel runs from our two-phase method is a cluster of approximate maximum likelihood estimates. Our results are illustrated by numerical examples.

Abdul Lateef Haji Ali, Fabio Nobile, Erik von Schwerin et al.

We perform a general optimization of the parameters in the Multilevel Monte Carlo (MLMC) discretization hierarchy based on uniform discretization methods with general approximation orders and computational costs. We optimize hierarchies with geometric and non-geometric sequences of mesh sizes and show that geometric hierarchies, when optimized, are nearly optimal and have the same asymptotic computational complexity as non-geometric optimal hierarchies. We discuss how enforcing constraints on parameters of MLMC hierarchies affects the optimality of these hierarchies. These constraints include an upper and a lower bound on the mesh size or enforcing that the number of samples and the number of discretization elements are integers. We also discuss the optimal tolerance splitting between the bias and the statistical error contributions and its asymptotic behavior. To provide numerical grounds for our theoretical results, we apply these optimized hierarchies together with the Continuation MLMC Algorithm. The first example considers a three-dimensional elliptic partial differential equation with random inputs. Its space discretization is based on continuous piecewise trilinear finite elements and the corresponding linear system is solved by either a direct or an iterative solver. The second example considers a one-dimensional Itô stochastic differential equation discretized by a Milstein scheme.

Abdul-Lateef Haji-Ali, Fabio Nobile, Raul Tempone

We propose and analyze a novel Multi-Index Monte Carlo (MIMC) method for weak approximation of stochastic models that are described in terms of differential equations either driven by random measures or with random coefficients. The MIMC method is both a stochastic version of the combination technique introduced by Zenger, Griebel and collaborators and an extension of the Multilevel Monte Carlo (MLMC) method first described by Heinrich and Giles. Inspired by Giles's seminal work, we use in MIMC high-order mixed differences instead of using first-order differences as in MLMC to reduce the variance of the hierarchical differences dramatically. This in turn yields new and improved complexity results, which are natural generalizations of Giles's MLMC analysis and which increase the domain of the problem parameters for which we achieve the optimal convergence, $\mathcal{O}(\text{TOL}^{-2}).$ Moreover, in MIMC, the rate of increase of required memory with respect to $\text{TOL}$ is independent of the number of directions up to a logarithmic term which allows far more accurate solutions to be calculated for higher dimensions than what is possible when using MLMC. We motivate the setting of MIMC by first focusing on a simple full tensor index set. We then propose a systematic construction of optimal sets of indices for MIMC based on properly defined profits that in turn depend on the average cost per sample and the corresponding weak error and variance. Under standard assumptions on the convergence rates of the weak error, variance and work per sample, the optimal index set turns out to be the total degree (TD) type. In some cases, using optimal index sets, MIMC achieves a better rate for the computational complexity than the corresponding rate when using full tensor index sets...

Alvaro Moraes, Raul Tempone, Pedro Vilanova

In this work, we extend the hybrid Chernoff tau-leap method to the multilevel Monte Carlo (MLMC) setting. Inspired by the work of Anderson and Higham on the tau-leap MLMC method with uniform time steps, we develop a novel algorithm that is able to couple two hybrid Chernoff tau-leap paths at different levels. Using dual-weighted residual expansion techniques, we also develop a new way to estimate the variance of the difference of two consecutive levels and the bias. This is crucial because the computational work required to stabilize the coefficient of variation of the sample estimators of both quantities is often unaffordable for the deepest levels of the MLMC hierarchy. Our method bounds the global computational error to be below a prescribed tolerance, $TOL$, within a given confidence level. This is achieved with nearly optimal computational work. Indeed, the computational complexity of our method is of order $\mathcal{O}\left(TOL^{-2}\right)$, the same as with an exact method, but with a smaller constant. Our numerical examples show substantial gains with respect to the previous single-level approach and the Stochastic Simulation Algorithm.

Hamidou Tembine, Raul Tempone, Pedro Vilanova

In this paper we study iterative procedures for stationary equilibria in games with large number of players. Most of learning algorithms for games with continuous action spaces are limited to strict contraction best reply maps in which the Banach-Picard iteration converges with geometrical convergence rate. When the best reply map is not a contraction, Ishikawa-based learning is proposed. The algorithm is shown to behave well for Lipschitz continuous and pseudo-contractive maps. However, the convergence rate is still unsatisfactory. Several acceleration techniques are presented. We explain how cognitive users can improve the convergence rate based only on few number of measurements. The methodology provides nice properties in mean field games where the payoff function depends only on own-action and the mean of the mean-field (first moment mean-field games). A learning framework that exploits the structure of such games, called, mean-field learning, is proposed. The proposed mean-field learning framework is suitable not only for games but also for non-convex global optimization problems. Then, we introduce mean-field learning without feedback and examine the convergence to equilibria in beauty contest games, which have interesting applications in financial markets. Finally, we provide a fully distributed mean-field learning and its speedup versions for satisfactory solution in wireless networks. We illustrate the convergence rate improvement with numerical examples.