Andrew M. Stuart

Inverse Problems: A Bayesian Perspective

The subject of inverse problems in differential equations is of enormous practical importance, and has also generated substantial mathematical and computational innovation. Typically some form of regularization is required to ameliorate ill-posed behaviour. In this article we review the Bayesian approach to regularization, developing a function space viewpoint on the subject. This approach allows for a full characterization of all possible solutions, and their relative probabilities, whilst simultaneously forcing significant modelling issues to be addressed in a clear and precise fashion. Although expensive to implement, this approach is starting to lie within the range of the available computational resources in many application areas. It also allows for the quantification of uncertainty and risk, something which is increasingly demanded by these applications. Furthermore, the approach is conceptually important for the understanding of simpler, computationally expedient approaches to inverse problems.

Strong Convergence of Euler-Type Methods for Nonlinear Stochastic Differential Equations

Traditional finite-time convergence theory for numerical methods applied to stochastic differential equations (SDEs) requires a global Lipschitz assumption on the drift and diffusion coefficients. In practice, many important SDE models satisfy only a local Lipschitz property and, since Brownian paths can make arbitrarily large excursions, the global Lipschitz-based theory is not directly relevant. In this work we prove strong convergence results under less restrictive conditions. First, we give a convergence result for Euler--Maruyama requiring only that the SDE is locally Lipschitz and that the pth moments of the exact and numerical solution are bounded for some p >2. As an application of this general theory we show that an implicit variant of Euler--Maruyama converges if the diffusion coefficient is globally Lipschitz, but the drift coefficient satisfies only a one-sided Lipschitz condition; this is achieved by showing that the implicit method has bounded moments and may be viewed as an Euler--Maruyama approximation to a perturbed SDE of the same form. Second, we show that the optimal rate of convergence can be recovered if the drift coefficient is also assumed to behave like a polynomial.

Extracting macroscopic dynamics: model problems and algorithms

In many applications, the primary objective of numerical simulation of timeevolving systems is the prediction of coarse-grained, or macroscopic, quantities. The purpose of this review is twofold: first, to describe a number of simple model systems where the coarse-grained or macroscopic behaviour of a system can be explicitly determined from the full, or microscopic, description; and second, to overview some of the emerging algorithmic approaches that have been introduced to extract effective, lower-dimensional, macroscopic dynamics. The model problems we describe may be either stochastic or deterministic in both their microscopic and macroscopic behaviour, leading to four possibilities in the transition from microscopic to macroscopic descriptions. Model problems are given which illustrate all four situations, and mathematical tools for their study are introduced. These model problems are useful in the evaluation of algorithms. We use specific instances of the model problems to illustrate these algorithms. As the subject of algorithm development and analysis is, in many cases, in its infancy, the primary purpose here is to attempt to unify some of the emerging ideas so that individuals new to the field have a structured access to the literature. Furthermore, by discussing the algorithms in the context of the model problems, a platform for understanding existing algorithms and developing new ones is built.

Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise

The ergodic properties of SDEs, and various time discretizations for SDEs, are studied. The ergodicity of SDEs is established by using techniques from the theory of Markov chains on general state spaces, such as that expounded by Meyn-Tweedie. Application of these Markov chain results leads to straightforward proofs of geometric ergodicity for a variety of SDEs, including problems with degenerate noise and for problems with locally Lipschitz vector fields. Applications where this theory can be usefully applied include damped-driven Hamiltonian problems (the Langevin equation), the Lorenz equation with degenerate noise and gradient systems. The same Markov chain theory is then used to study time-discrete approximations of these SDEs. The two primary ingredients for ergodicity are a minorization condition and a Lyapunov condition. It is shown that the minorization condition is robust under approximation. For globally Lipschitz vector fields this is also true of the Lyapunov condition. However in the locally Lipschitz case the Lyapunov condition fails for explicit methods such as Euler-Maruyama; for pathwise approximations it is, in general, only inherited by specially constructed implicit discretizations. Examples of such discretization based on backward Euler methods are given, and approximation of the Langevin equation studied in some detail.

MCMC Methods for Functions: Modifying Old Algorithms to Make Them Faster

Many problems arising in applications result in the need to probe a probability distribution for functions. Examples include Bayesian nonparametric statistics and conditioned diffusion processes. Standard MCMC algorithms typically become arbitrarily slow under the mesh refinement dictated by nonparametric description of the un- known function. We describe an approach to modifying a whole range of MCMC methods, applicable whenever the target measure has density with respect to a Gaussian process or Gaussian random field reference measure, which ensures that their speed of convergence is robust under mesh refinement. Gaussian processes or random fields are fields whose marginal distri- butions, when evaluated at any finite set of N points, are RN-valued Gaussians. The algorithmic approach that we describe is applicable not only when the desired probability measure has density with respect to a Gaussian process or Gaussian random field reference measure, but also to some useful non-Gaussian reference measures constructed through random truncation. In the applications of interest the data is often sparse and the prior specification is an essential part of the over- all modelling strategy. These Gaussian-based reference measures are a very flexible modelling tool, finding wide-ranging application. Examples are shown in density estimation, data assimilation in fluid mechanics, subsurface geophysics and image registration. The key design principle is to formulate the MCMC method so that it is, in principle, applicable for functions; this may be achieved by use of proposals based on carefully chosen time-discretizations of stochas- tic dynamical systems which exactly preserve the Gaussian reference measure. Taking this approach leads to many new algorithms which can be implemented via minor modification of existing algorithms, yet which show enormous speed-up on a wide range of applied problems.

The global dynamics of discrete semilinear parabolic equations

A class of scalar semilinear parabolic equations possessing absorbing sets, a Lyapunov functional, and a global attractor are considered. The gradient structure of the problem implies that, provided all steady states are isolated, solutions approach a steady state as

Space-Time Continuous Analysis of Waveform Relaxation for the Heat Equation

t \to \infty

Optimal tuning of the hybrid Monte Carlo algorithm

. The dynamical properties of various finite difference and finite element schemes for the equations are analysed. The existence of absorbing sets, bounded independently of the mesh size, is proved for the numerical methods. Discrete Lyapunov functions are constructed to show that, under appropriate conditions on the mesh parameters, numerical orbits approach steady state solutions as discrete time increases. However, it is shown that insufficient spatial resolution can introduce deceptively smooth spurious steady solutions and cause the stability properties of the true steady solutions to be incorrectly represented. Furthermore, it is also shown that the explicit Euler scheme introduces spurious solutions with period 2 in the timestep. As a result, the absorbing set is destroyed and there is initial data leading to blow up of the scheme, however small the mesh parameters are taken. To obtain stabilization to a steady state for this scheme, it is necessary to restrict the timestep in terms of the initial data and the space step. Implicit schemes are constructed for which absorbing sets and Lyapunov functions exist under restrictions on the timestep that are independent of initial data and of the space step; both one-step and multistep (BDF) methods are studied.

Runge-Kutta methods for dissipative and gradient dynamical systems

Waveform relaxation algorithms for partial differential equations (PDEs) are traditionally obtained by discretizing the PDE in space and then splitting the discrete operator using matrix splittings. For the semidiscrete heat equation one can show linear convergence on unbounded time intervals and superlinear convergence on bounded time intervals by this approach. However, the bounds depend in general on the mesh parameter and convergence rates deteriorate as one refines the mesh. Motivated by the original development of waveform relaxation in circuit simulation, where the circuits are split in the physical domain into subcircuits, we split the PDE by using overlapping domain decomposition. We prove linear convergence of the algorithm in the continuous case on an infinite time interval, at a rate depending on the size of the overlap. This result remains valid after discretization in space and the convergence rates are robust with respect to mesh refinement. The algorithm is in the class of waveform relaxation algorithms based on overlapping multisplittings. Our analysis quantifies the empirical observation by Jeltsch and Pohl [SIAM J. Sci. Comput., 16 (1995), pp. 40--49] that the convergence rate of a multisplitting algorithm depends on the overlap. Numerical results are presented which support the convergence theory.

Exponential Mean-Square Stability of Numerical Solutions to Stochastic Differential Equations

We investigate the properties of the Hybrid Monte Carlo algorithm (HMC) in high dimensions. HMC develops a Markov chain reversible w.r.t. a given target distribution . by using separable Hamiltonian dynamics with potential -log .. The additional momentum variables are chosen at random from the Boltzmann distribution and the continuous-time Hamiltonian dynamics are then discretised using the leapfrog scheme. The induced bias is removed via a Metropolis- Hastings accept/reject rule. In the simplified scenario of independent, identically distributed components, we prove that, to obtain an O(1) acceptance probability as the dimension d of the state space tends to ., the leapfrog step-size h should be scaled as h=l ×d−1/ 4 . Therefore, in high dimensions, HMC requires O(d1/ 4 ) steps to traverse the state space. We also identify analytically the asymptotically optimal acceptance probability, which turns out to be 0.651 (to three decimal places). This is the choice which optimally balances the cost of generating a proposal, which decreases as l increases (because fewer steps are required to reach the desired final integration time), against the cost related to the average number of proposals required to obtain acceptance, which increases as l increases