A. M. Stuart

K. J. H. Law, A. M. Stuart

Data assimilation leads naturally to a Bayesian formulation in which the posterior probability distribution of the system state, given the observations, plays a central conceptual role. The aim of this paper is to use this Bayesian posterior probability distribution as a gold standard against which to evaluate various commonly used data assimilation algorithms. A key aspect of geophysical data assimilation is the high dimensionality and low predictability of the computational model. With this in mind, yet with the goal of allowing an explicit and accurate computation of the posterior distribution, we study the 2D Navier-Stokes equations in a periodic geometry. We compute the posterior probability distribution by state-of-the-art statistical sampling techniques. The commonly used algorithms that we evaluate against this accurate gold standard, as quantified by comparing the relative error in reproducing its moments, are 4DVAR and a variety of sequential filtering approximations based on 3DVAR and on extended and ensemble Kalman filters. The primary conclusions are that: (i) with appropriate parameter choices, approximate filters can perform well in reproducing the mean of the desired probability distribution; (ii) however they typically perform poorly when attempting to reproduce the covariance; (iii) this poor performance is compounded by the need to modify the covariance, in order to induce stability. Thus, whilst filters can be a useful tool in predicting mean behavior, they should be viewed with caution as predictors of uncertainty. These conclusions are intrinsic to the algorithms and will not change if the model complexity is increased, for example by employing a smaller viscosity, or by using a detailed NWP model.

R. Scheichl, A. M. Stuart, A. L. Teckentrup

We are interested in computing the expectation of a functional of a PDE solution under a Bayesian posterior distribution. Using Bayes' rule, we reduce the problem to estimating the ratio of two related prior expectations. For a model elliptic problem, we provide a full convergence and complexity analysis of the ratio estimator in the case where Monte Carlo, quasi-Monte Carlo or multilevel Monte Carlo methods are used as estimators for the two prior expectations. We show that the computational complexity of the ratio estimator to achieve a given accuracy is the same as the corresponding complexity of the individual estimators for the numerator and the denominator. We {also include numerical simulations, in the context of the model elliptic problem, which demonstrate the effectiveness of the approach.

G. A. Pavliotis, A. M. Stuart, K. C. Zygalakis

In this paper we study the problem of the numerical calculation (by Monte Carlo Methods) of the effective diffusivity for a particle moving in a periodic divergent-free velocity filed, in the limit of vanishing molecular diffusion. In this limit traditional numerical methods typically fail, since they do not represent accurately the geometry of the underlying deterministic dynamics. We propose a stochastic splitting method that takes into account the volume preserving property of the equations motion in the absence of noise, and when inertial effects can be neglected. An extension of the method is then proposed for the cases where the noise has a non trivial time-correlation structure and when inertial effects cannot be neglected. Modified equations are used to perform backward error analysis. The new stochastic geometric integrators are shown to outperform standard Euler-based integrators. Various asymptotic limits of physical interest are investigated by means of numerical experiments, using the new integrators.

S. L. Cotter, M. Dashti, A. M. Stuart

Inverse problems are often ill-posed, with solutions that depend sensitively on data. In any numerical approach to the solution of such problems, regularization of some form is needed to counteract the resulting instability. This paper is based on an approach to regularization, employing a Bayesian formulation of the problem, which leads to a notion of well-posedness for inverse problems, at the level of probability measures. The stability which results from this well-posedness may be used as the basis for quantifying the approximation, in finite dimensional spaces, of inverse problems for functions. This paper contains a theory which utilizes the stability to estimate the distance between the true and approximate posterior distributions, in the Hellinger metric, in terms of error estimates for approximation of the underlying forward problem. This is potentially useful as it allows for the transfer of estimates from the numerical analysis of forward problems into estimates for the solution of the related inverse problem. In particular controlling differences in the Hellinger metric leads to control on the differences between expected values of polynomially bounded functions and operators, including the mean and covariance operator. The ideas are illustrated with the classical inverse problem for the heat equation, and then applied to some more complicated non-Gaussian inverse problems arising in data assimilation, involving determination of the initial condition for the Stokes or Navier-Stokes equation from Lagrangian and Eulerian observations respectively.

H. Lamba, J. C. Mattingly, A. M. Stuart

The understanding of adaptive algorithms for SDEs is an open area where many issues related to both convergence and stability (long time behaviour) of algorithms are unresolved. This paper considers a very simple adaptive algorithm, based on controlling only the drift component of a time-step. Both convergence and stability are studied. The primary issue in the convergence analysis is that the adaptive method does not necessarily drive the time-steps to zero with the user-input tolerance. This possibility must be quantified and shown to have low probability. The primary issue in the stability analysis is ergodicity. It is assumed that the noise is non-degenerate, so that the diffusion process is elliptic, and the drift is assumed to satisfy a coercivity condition. The SDE is then geometrically ergodic (converges to statistical equilibrium exponentially quickly). If the drift is not linearly bounded then explicit fixed time-step approximations, such as the Euler-Maruyama scheme, may fail to be ergodic. In this work, it is shown that the simple adaptive time-stepping strategy cures this problem. In addition to proving ergodicity, an exponential moment bound is also proved, generalizing a result known to hold for the SDE itself.