MATHICSE Technical Report : A continuation multilevel Monte Carlo algorithm

We propose a novel Continuation Multi Level Monte Carlo (CMLMC) algorithm for weak approximation of stochastic models that are described in terms of differential equations either driven by random measures or with random coefficients. The CMLMC algorithm solves the given approximation problem for a sequence of decreasing tolerances, ending with the desired one. CMLMC assumes discretization hierarchies that are defined a priori for each level and are geometrically refined across levels. The actual choice of computational work across levels is based on parametric models for the average cost per sample and the corresponding weak and strong errors. These parameters are calibrated using Bayesian estimation, taking particular notice of the deepest levels of the discretization hierarchy, where only few realizations are available to produce the estimates. The resulting CMLMC estimator exhibits a nontrivial splitting between bias and statistical contributions. We also show the asymptotic normality of the statistical error in the MLMC estimator and justify in this way our error estimate that allows prescribing both required accuracy and confidence in the final result. Numerical examples substantiate the above results and illustrate the corresponding computational savings.

A continuation multilevel Monte Carlo algorithm

Fabio Nobile, Erik Gustaf Bogislaw Von Schwerin

We propose a novel Continuation Multi Level Monte Carlo (CMLMC) algorithm for weak approximation of stochastic models. The CMLMC algorithm solves the given approximation problem for a sequence of decreasing tolerances, ending when the required error tolerance is satisfied. CMLMC assumes discretization hierarchies that are defined a priori for each level and are geometrically refined across levels. The actual choice of computational work across levels is based on parametric models for the average cost per sample and the corresponding variance and weak error. These parameters are calibrated using Bayesian estimation, taking particular notice of the deepest levels of the discretization hierarchy, where only few realizations are available to produce the estimates. The resulting CMLMC estimator exhibits a non-trivial splitting between bias and statistical contributions. We also show the asymptotic normality of the statistical error in the MLMC estimator and justify in this way our error estimate that allows prescribing both required accuracy and confidence in the final result. Numerical results substantiate the above results and illustrate the corresponding computational savings in examples that are described in terms of differential equations either driven by random measures or with random coefficients.

2015

Probabilistic and Bayesian methods for uncertainty quantification of deterministic and stochastic differential equations

Giacomo Garegnani

In this thesis we explore uncertainty quantification of forward and inverse problems involving differential equations. Differential equations are widely employed for modeling natural and social phenomena, with applications in engineering, chemistry, meteorology, and economics. Mathematical models of complex systems in these fields require numerical methods, which introduce uncertainties in the outcome. Moreover, there has recently been a steep rise in the availability of data, which also come with an uncertainty. Therefore, blending mathematical models, data, and their respective uncertainties is nowadays of the utmost importance. The first part of this thesis is dedicated to two novel methods for multiscale inverse problems. We first consider an elliptic partial differential equation (PDE), with a diffusion tensor oscillating at a small scale. Given noisy observations of the solution, we consider the problem of inferring a slow-scale parametrization of the multiscale tensor. For this purpose, we combine numerical homogenization, which yields a single-scale surrogate of the full model, and the ensemble Kalman filter. The scheme we propose is accurate in the homogenized limit, and outperforms existing methods in terms of computational cost. We then study the error due to the mismatch between the full and the homogenized models, and show how to combine statistical techniques for model misspecification and our scheme. We then move to multiscale diffusion processes, and consider the problem of inferring effective dynamics from multiscale observations. A homogenized single-scale equation reproducing the full model exists also in this case. The effective model, though, is the subject of the inference procedure, and not only a computational tool. The resulting issue of model misspecification is usually bypassed by subsampling at an appropriate rate, which is non-trivial to choose, and which may give misleading results. We avoid subsampling by designing a novel technique based on filtered data, and show how to modify classical estimators and obtain an effective equation consistently with homogenization. Our technique is robust and can be employed as a black-box tool for inferring effective surrogates of complex stochastic models.In the second part we present two novel schemes belonging to the field of probabilistic numerics, whose purpose is to provide a statistical description of the uncertainty due to numerical discretization. We first consider ordinary differential equations (ODEs), and introduce a probabilistic integrator based on random time steps and Runge-Kutta methods (RTS-RK). Tuning the distribution of the time steps, we generate a probability measure on the solution which allows for a consistent uncertainty quantification of numerical errors. Unlike previous probabilistic methods in literature, our scheme inherits the geometric properties of the underlying deterministic integrators. In particular, we show long-time energy conservation when the RTS-RK is applied to Hamiltonian ODEs. We employ the idea of randomizing the discretization to propose a random mesh finite element method (RM-FEM) for elliptic PDEs. We prove that the measure induced by the RM-FEM on the solution can be employed to derive a posteriori error estimators. Hence, the RM-FEM provides a consistent statistical characterization of numerical errors. For both our novel schemes, we demonstrate the usefulness of the probabilistic approach in Bayesian inverse problems.

EPFL2021

Propagation Models for Biochemical Reaction Networks

Maria-Emanuela-Canini Mateescu

In this thesis we investigate different ways of approximating the solution of the chemical master equation (CME). The CME is a system of differential equations that models the stochastic transient behaviour of biochemical reaction networks. It does so by describing the time evolution of probability distribution over the states of a Markov chain that represents a biological network, and thus its stochasticity is only implicit. The transient solution of a CME is the vector of probabilities over the states of the corresponding Markov chain at a certain time point t, and it has traditionally been obtained by applying methods that are general to continuous-time Markov chains: uniformization, Krylov subspace methods, and general ordinary differential equation (ODE) solvers such as the fourth order Runge-Kutta method. Even though biochemical reaction networks are the main application of our work, some of our results are presented in the more general framework of propagation models (PM), a computational formalism that we introduce in the first part of this thesis. Each propagation model N has two associated propagation processes, one in discrete-time and a second one in continuous-time. These propagation processes propagate a generic mass through a discrete state space. For example, in order to model a CME, N propagates probability mass. In the discrete-time case the propagation is done step-wise, while in the continuous-time case it is done in a continuous flow defined by a differential equation. Again, in the case of the chemical master equation, this differential equation is the equivalent of the chemical master equation itself where probability mass is propagated through a discrete state space. Discrete-time propagation processes can encode methods such as the uniformization method and the fourth order Runge-Kutta integration method that we have mentioned above, and thus by optimizing propagation algorithms we optimize both of these methods simultaneously. In the second part of our thesis, we define stochastic hybrid models that approximate the stochastic behaviour of biochemical reaction networks by treating some variables of the system deterministically. This deterministic approximation is done for species with large populations, for which stochasticity does not play an important role. We propose three such hybrid models, which we introduce from the coarsest to the most refined one: (i) the first one replaces some variables of the system with their overall expectations, (ii) the second one replaces some variables of the system with their expectations conditioned on the values of the stochastic variables, (iii) and finally, the third one, splits each variable into a stochastic part (for low valuations) and a deterministic part (for high valuations), while tracking the conditional expectation of the deterministic part. For each of these algorithms we give the corresponding propagation models that propagate not only probabilities but also the respective continuous approximations for the deterministic variables.