Bernstein–von Mises theorems for
time evolution equations

Denote by $L^{2}(\Omega)$ the Hilbert space of square-integrable vector fields over a bounded domain $\Omega\subset\mathbb{R}^{d}$. For initial condition $\theta=u(0)$, consider states $u_{\theta}(t)$ at times $t>0$ of a dynamical system evolving in $L^{2}(\Omega)$. Examples we have in mind include the Navier-Stokes, reaction-diffusion or McKean-Vlasov systems [43, 39, 9, 41], where $u_{\theta}$ models the velocity field of an incompressible fluid, the concentration of a chemical substance, or the distribution of interacting particles, respectively. The infinitesimal dynamics of such systems are described by a non-linear partial differential equation (PDE)

where $\Delta$ is the Laplacian and $F$ is a functional modelling the underlying non-linearity.

In scientific applications, the initial condition $\theta$ is uncertain and moreover, physical observations are necessarily discrete and include statistical error. For instance, following the paradigm of probabilistic numerics [13, 10, 40, 3], an experimenter may sample the state $u_{\theta}$ of the system at positions $(t_{i},\omega_{i})_{i=1}^{N}$ drawn uniformly and independently at random from the time-space cylinder $[0,T]\times\Omega$, with $T>0$ a fixed time horizon. One thus observes the vector $Z^{(N)}=(Y_{i},t_{i},\omega_{i})_{i=1}^{N}$ from the regression equation

where the $\varepsilon_{i}\sim\mathcal{N}(0,I)$ are independent and identically distributed (iid) Gaussian random noise. The infinite product probability measure describing the law of the sequence $Z^{(\infty)}$ arising from initial condition $\theta$ will be denoted by $P_{\theta}^{\mathbb{N}}$.

The uncertainty about the state $\theta$ of the system before measurements are taken can be updated to a posterior distribution on the $u_{\theta}$’s given $Z^{(N)}$, by applying the rules of conditional probability. This follows seminal ideas of Laplace (Chapter VI in [24]) and nowadays is called Bayesian inference [19] or sometimes more specifically – when considering dynamical systems in applied sciences – ‘data assimilation’ or ‘filtering’ [16, 38, 25, 26, 23]. In infinite-dimensional settings where $\theta$ is a function, a typical approach is to model the initial condition by a Gaussian random field $(\theta(x):x\in\Omega)$ whose mean function is based on past knowledge (or in absence of it, equals zero) and whose covariance is an inverse power of the Laplace operator on $\Omega$. If we denote by $\Pi$ the induced ‘prior’ probability distribution on $\theta$, the updating step (i.e., Bayes’ theorem) generates the posterior Gibbs’ probability distribution for $\theta|Z^{(N)}$ of the form

We also obtain an update for the whole dynamical system $u_{\theta}(t)$ whose marginal time distributions on $L^{2}(\Omega)$ are given by the image measures

For the last observation times $t_{(N)}=\max_{i\leq N}t_{i}\in(0,T],$ the preceding laws are sometimes called filtering distributions; for $t>t_{(N)}$ they can be used to forecast future states of the system, while for $t<t_{(N)}$ the are referred to as ‘smoothing distributions’, see, e.g., p.xii in [25] or p.173 in [38] for this terminology from the data assimilation literature. Various iterative algorithms exist that aim to compute these posterior measures by numerical (e.g., MCMC, or particle filter) methods [40, 11, 22, 38, 25, 4, 21]. The resulting inferences are widely used in scientific application areas, but statistical theory for the performance of such methods in non-linear settings remains elusive.

The distributions $\hat{\Pi}_{t,N}$ are random probability measures in function space that involve the non-linearity of the underlying evolution operator both in the Gibbs’ measure (3) and in the push-forward (4), and have no simple closed form representation. A first sanity check for this Bayesian approach would be to establish its ‘posterior consistency’ [19]; namely that $\hat{\Pi}_{t,N}$ converges in an appropriate sense to Dirac measure $\delta_{u_{\theta_{0}}}$ at the state $u_{\theta_{0}}(t)$ arising from ‘ground truth’ initial condition $\theta_{0}$, at least as the signal to noise ratio increases (i.e., as sample size $N\to\infty$), and with high $P_{\theta_{0}}^{\mathbb{N}}$-probability. For the $2d$-Navier-Stokes model such a result was proved recently in [36] – see Remark 2 for more discussion.

In this article we aim to further advance our understanding of the statistical behaviour of the random measures $\hat{\Pi}_{t,N}$ for large sample sizes $N\to\infty$. Specifically, following another idea of Laplace, we investigate whether the non-Gaussian posterior measures $\hat{\Pi}_{t,N}$ are perhaps approximately Gaussian in a suitable sense, a phenomenon that goes by the name of Bernstein-von Mises (BvM) theorem in the mathematical statistics literature. Such results hold in ‘regular’ finite-dimensional statistical models [44] but are more intricate in high- or infinite-dimensional settings [18, 5, 6, 19, 33]. We will devise a template tailored to dynamical systems of parabolic type (1) that allows to show that in a suitable sense, uniformly in fixed time windows, and with high $P_{\theta_{0}}^{\mathbb{N}}$-probability,

where $\tilde{\theta}_{N}$ is the posterior mean vector and where the limit law $\mathcal{N}(0,\mathcal{C}_{t})$ can be characterised by the solution at time $t$ of a linear stochastic PDE with certain ‘rough’ Gaussian initial conditions. We demonstrate how the theory works for a concrete nonlinear reaction-diffusion system and discuss structurally similar non-linear parabolic PDEs (1) in Remark 2.

These Gaussian approximations will be shown to hold in the uniform norm over $\Omega$, and at the convergence rate $1/\sqrt{N}$. This is possible by first proving a Gaussian approximation in weak ‘Schwartz’-type topologies for the initial conditions $\theta$ (partly following ideas from [5, 6, 30] in different models), and then using strong smoothing properties of the semigroup underlying (1) at positive times $t>0$. Some discussion of the last point and its relationship to non-existence results [18] for infinite-dimensional Bernstein-von Mises theorems can be found in Remark 3.

We let $\Omega=[0,1]^{d}$ denote the $d$-dimensional torus and consider periodic solutions to the PDE

where $\Delta$ is the Laplacian and $f$ is a given non-linear ‘reaction’ term that is smooth, compactly supported and time-independent, acting on $u$ by composition $f(u)=f\circ u$ point-wise. The unique solutions of this PDE induce a dynamical system evolving in $L^{2}(\Omega)$, see [43, 39]. In applications the initial condition $\theta$ may represent the distribution of some substance and (6) describes its evolution over time where the non-linear reaction term models the amount of the substance that is created or destroyed in dependence $f(u(t))$ of the current state $u(t)$. For the proofs we will assume the dimension of the state space $\Omega=[0,1]^{d}$ to be $d\leq 3$ and we also only consider scalar $u$ taking values in $\mathbb{R}$ rather than general vector fields $u:\Omega\to\mathbb{R}^{d}$. Extensions to $d\geq 4$ are mostly of a technical nature (involving Schauder theory, as in Sec. 5.2 of [27]) while coupled systems of equations could be dealt with by mostly notational changes as long as $f=\nabla F$ is a gradient vector field of smooth $F:\Omega\to\mathbb{R}$.

We model the uncertainty about the initial condition $\theta$ of the system by a $\gamma$-regular Gaussian random field $(\theta(x):x\in\Omega)$ whose prior $Law(\theta)$ on the space $L^{2}_{0}(\Omega)=L^{2}(\Omega)\cap\{h:\int_{\Omega}h=0\}$ arises from

see Condition 1 for details. The choice of $\rho$ provides necessary prior regularisation to prevent the nonlinearities inherent in the posterior Gibbs measure (3) to behave erratically at ‘low temperatures’ $N\to\infty$. A way that achieves this exhibited in [28] takes the form

The posterior distribution $\Pi(\theta|Z^{(N)})$ and its pushforward $\hat{\Pi}_{t,N}$ then arise from data (2) as in (3) and (4), with dynamical system described by the PDE (6) for any fixed $f\in C^{\infty}_{c}(\Omega)$.

Our main Theorem 1 will show that the (non-Gaussian) posterior measures $\hat{\Pi}_{t,N}$ are approximated by a Gaussian law, as announced in (5). This approximation will hold as $N\to\infty$, and in $P_{\theta_{0}}^{\mathbb{N}}$-probability, with ground truth initial condition $\theta_{0}$ that generated the data (2). The most delicate step in the identification of the limit distribution is the construction of a mean zero Gaussian random field $(\vartheta(x):x\in\Omega)$ whose covariance is related to an appropriate inverse of the underlying ‘Fisher information’ operator of the statistical model (see, [44] or Ch.3 in [31] for more on this notion). We will show that this process indeed exists for the reaction-diffusion system (6) and that the law $\mathcal{N}_{\theta_{0}}$ it induces defines a Gaussian Borel measure in the negative Sobolev space $H^{-a}$ whenever $a>1+d/2$. Then let $u_{\theta_{0}}$ be the solution of (6), write $f^{\prime}=df/dx$, and consider the Gaussian process $U$ over $(0,T]\times\Omega$ obtained from the unique weak solution of the following linear time-dependent Schrödinger equation with random initial condition:

Even though the initial condition $\vartheta$ is not point-wise defined as a function (almost surely), our proofs will imply that the weak solutions $U(t,\cdot)$ exist and almost surely define continuous functions on $\Omega$ for $t>0$. This is related to the time-smoothing properties of the parabolic solution operator underlying (1.2) and the fact that the null space of the Schrödinger operator $\Delta+f^{\prime}(\theta_{0})$ governing the ‘explosion’ at time $t=0$ will be seen to be finite-dimensional, see also Remark 3.

We will assume that the true initial condition $\theta_{0}$ (but not the prior model) is smooth to simplify the statement of the following result. But this restriction is not necessary and can be weakened to sufficient Sobolev-regularity of $\theta_{0}$. Fix any time window $0<t_{min}<t_{max}<\infty$ and define the Banach space

of continuous maps on $[t_{\min},t_{\max}]\times\Omega$. Denote by $\mathscr{W}_{1}$ (equal to $W_{1,\mathscr{C}}$ in (20) below) the $1$-Wasserstein distance on the space of probability measures on $\mathscr{C}$.

From convergence in law $\to^{d}$ in the space $\mathscr{C}$ one deduces in particular that $\sqrt{N}(u_{\tilde{\theta}}-u_{\theta_{0}})$ is uniformly tight in $\mathscr{C}$, so the theorem implies that the updated posterior mean estimates $u_{\tilde{\theta}}(t)$ are concentrated near the ‘true’ dynamical system $u_{\theta_{0}}(t)$, and that their uniform $\|\cdot\|_{\mathscr{C}}$ fluctuations at a $\sqrt{N}$ scale are approximately Gaussian.

Theorem 1 has various important applications that we only outline briefly here. One is to uncertainty quantification and the frequentist validity of posterior credible bands which are random subsets $C_{N}\subset\mathscr{C}$ of posterior probability $1-\alpha$ for some fixed $0<\alpha<1$. Theorem 1 and arguments as in the proof of Theorem 7.3.21 in [20] can be combined to show that appropriate $C_{N}$ are of diameter $O_{P_{\theta_{0}}^{\mathbb{N}}}(1/\sqrt{N})$ and satisfy

Another interesting application our our results concerns guarantees for posterior computation: When Gaussian approximation theorems such as Theorem 1 hold, they furnish a strategy to prove the existence of polynomial-time computational algorithms based on MCMC [2, 22, 37, 31], overcoming potential computational hardness barriers [1]. We note that the hypothesis of a Gaussian approximation is also used in recent work on convergence guarantees for Kalman-type filters in non-linear settings [4].

For real numbers $a,b$, we write $a\lesssim b$ whenever $a\leq Cb$ for some fixed constant $C>0$, and $a\simeq b$ if $a\lesssim b,b\lesssim a$. We write $Z\sim\mu$ when a random variable has law $\mu$, write $\to^{P}$ for convergence in probability, $\to^{d}$ for convergence in distribution, and use the standard $O_{P},o_{P}$ notation for stochastic orders of magnitude, see [44].

For $1\leq p\leq\infty$ and $(\mathscr{Z},\mathcal{A},\mu)$ a measure space, the $L^{p}(\mu)$-spaces are defined in the usual way as all measurable functions $H$ on $\mathscr{Z}$ such that $|H|^{p}$ is $\mu$-integrable, or essentially bounded if $p=\infty$. When $\mathscr{Z}$ is a subset of $\mathbb{R}^{d}$ we take $\mu$ to be Lebesgue measure on the Borel sets, unless specified otherwise. For $X$ a normed linear space, we define the function space $L^{p}([0,T],X)$ of measurable maps $H$ from $[0,T]\to X$ such that $\|H(t)\|_{X}$ lies in $L^{p}([0,T])$. We also define the space $C([0,T],X)$ of continuous maps from $[0,T]\to X$ normed by the supremum norm $\sup_{0<t<T}\|H(t)\|_{X}$ as well as the spaces of weakly differentiable maps

For $\Omega=[0,1]^{d}$ the $d$-dimensional torus (with opposite points identified), $C(\Omega)$ denotes the Banach space of bounded continuous functions defined on $\Omega$, equipped with supremum norm $\|\cdot\|_{\infty}$, and $C^{\gamma}(\Omega)$ denotes the usual spaces of functions defined on $\Omega$ whose partial derivatives up to order $\gamma\in\mathbb{N}$ are bounded functions; for $\gamma\notin\mathbb{N}$ these are the Hölder spaces, and $C^{\infty}(\Omega)=\cap_{\gamma>0}C^{\gamma}(\Omega)$ denotes the set of smooth (infinitely-differentiable) functions defined on $\Omega$. We also introduce the spaces

Further define the Sobolev spaces $H^{\gamma}(\Omega),\gamma\geq 0,$ of functions whose weak partial derivatives up to order $\gamma$ lie in $L^{2}(\Omega)$. The norms $\|\cdot\|_{C^{\gamma}},\|\cdot\|_{H^{\gamma}}$ on $C^{\gamma},H^{\gamma},\gamma\in\mathbb{N},$ are then given by $\|H\|+\sum_{|\alpha|=\gamma}\|D^{\alpha}H\|,$ where $D^{\alpha}$ is the (weak) partial differential operator for multi-index $\alpha$ and the norms $\|\cdot\|$ equal $\|\cdot\|_{\infty}$ or $\|\cdot\|_{L^{2}}$, respectively.

The Sobolev spaces $H^{\gamma}$ carry equivalent sequence space norms

where $\langle\cdot,\cdot\rangle_{L^{2}}$ is the inner product of $L^{2}(\Omega)$, and the $e_{j}$ are the $L^{2}(\Omega)$-orthonormal eigenfunctions of the periodic Laplacian $\Delta$ for (negative) eigenvalues $\lambda_{j}$. Concretely we have $e_{0}=1,\lambda_{0}=0$, and

where $j=0,1,2,\dots$ is an enumeration of the integer lattice $\mathbb{Z}^{d}$ and

An equivalent norm is obtained by replacing $(1+\lambda_{j})^{\gamma}$ by $(1+\lambda_{j}^{\gamma})$ and this norm is further equivalent to the graph norm of the image of $L^{2}$ under $(id-\Delta^{\gamma/2})$.

We define topological dual spaces $H^{-\gamma}(\Omega)=(H^{\gamma}(\Omega))^{*}$ with norm

where the supremum may be restricted to $\psi\in C^{\infty}(\Omega)$. If we understand $\langle h,e_{j}\rangle_{L^{2}}$ as the action $T_{h}(e_{j})$ of periodic Schwartz distributions $T_{h}$ on smooth test functions $e_{j}\in C^{\infty}(\Omega)$, then the $h^{-\gamma}$ norms from (13) are equivalent to these dual norms also for negative $-\gamma<0$, and (16) is valid for all $\gamma\in\mathbb{R}$. We define closed subspaces

and note that $H^{0}_{0}(\Omega)=L^{2}_{0}(\Omega)$. On $H^{\gamma}_{0}(\Omega)$ we have equivalent norms as in (13) but using only $\lambda_{j}^{\gamma}$ in place of $(1+\lambda_{j})^{\gamma}$, taking note of the Poincaré inequality $\lambda_{1}\geq 1/2\pi$.

If we set $B^{a}=H^{a}$ for $a>d/2$ and $B^{a}=C^{a}$ for $0\leq a\leq d/2$, then the multiplier inequality for Sobolev norms is

For $a<0$ we can use (16) and (17) to obtain

We also need the following interpolation inequality for Sobolev norms (see (A.16), p.473, in [42] and Lemma 3.27 in [39])

To metrise the distance between two probability measures $\mu,\nu$ on a metric space $(X,d)$, we will use the transportation Wasserstein-$1$ distance,

where $Z_{1}\sim\mu,Z_{2}\sim\nu$ and

See Chapter 6 in [46] for equivalent definitions and further properties.

Let $\Omega=[0,1]^{d},d\in\mathbb{N}$, and let $W$ be a finite-dimensional vector space. If each entry $u_{k}$ of a $W$-valued vector field $u$ lies in a function space $X(\Omega)$ we write $u\in X(\Omega,W)$, and the corresponding norm is denoted by $\|u\|^{2}_{X(\Omega,W)}=\sum_{k}\|u_{k}\|_{X(\Omega)}^{2}$. Often we will just write $X(\Omega)$ for $X(\Omega,W)$. We consider a general forward map

When $W=\mathbb{R}$ this accommodates the solution map of the reaction-diffusion system (6), but further allows dynamical systems of vector fields, such as coupled systems of PDEs. The measurement model (2) is then a special case of the random design regression equation

where the $X_{i}$ are drawn iid from the uniform distribution $\lambda$ on $\mathcal{X}=[0,T]\times\Omega$, independently of the Gaussian noise vectors $\varepsilon_{i}$. Just as after (2) we write $Z^{(N)}:=(Y_{i},X_{i})_{i=1}^{N}$, denote the product measure describing the law of the infinite observation vector $(Y_{i},X_{i})_{i=1}^{\infty}$ by $P_{\theta}^{\mathbb{N}}$, and occasionally write $L^{2}(\mathcal{X})$ for $L^{2}([0,T],L^{2}(\Omega))$. This random design regression setting permits to borrow from the theory developed in [31] as well as the use of tools from concentration of product measures in high dimensions. Modelling explicit dependence structures in the design complicates the development significantly but is in principle possible, for instance as in [32].

We now turn to the prior probability measure for the parameter $\theta\in H^{1}$ – we refer to Ch.2 in [20] for standard background from the theory of Gaussian processes, such as the definition of their reproducing kernel Hilbert spaces (RKHS).

For $e_{j}$ from (15), the prior $\Pi^{\prime}_{\gamma}$ can be represented by the law of a Gaussian random series

augmented to independent such copies in each of its coordinates whenever $dim(W)>1$. Shrinking the prior towards zero ensures a degree of a priori regularisation and has been used throughout proofs in the non-linear inverse problems literature since [28]. By the theory of radonifying maps (as in Thm B.1.3 in [31]), $\Pi^{\prime}_{\gamma}$ and $\Pi$ define centred Gaussian Borel probability measures on the separable Hilbert space

which serves as the natural ‘parameter’ space charged almost surely by prior draws. The posterior law of $\theta|Z^{(N)}$ then arises from a dominated model and is given by

see [19] or Sec.1.2.3 in [31].

We now formulate some analytical conditions on $\mathscr{G}$ – these are designed for solution maps $\theta\mapsto\mathscr{G}(\theta)$ of non-linear parabolic PDEs as in (1). We denote the balls in $H^{r}_{0}$ by

and the regularity properties will be required to hold uniformly on such balls with constants that may depend on $B$ and the fixed ground truth $\theta_{0}$. The constants may depend on further parameters $T,W,d,\bar{\gamma},\gamma,\zeta$ which we however suppress in the notation.

We now state a Bernstein-von Mises theorem for the posterior on $\theta$ in the ‘weak’ norm topology of the space $H^{-k-2}$. The posterior mean vector $\tilde{\theta}_{N}=E^{\Pi}[\theta|Z^{(N)}]$ exists as a Bochner integral in $L^{2}_{0}$ and hence also in $H^{-\kappa-2}$, while the limiting Gaussian measure $\mathcal{N}_{\theta_{0}}$ is constructed in Proposition 2 as a tight Borel law on $H^{-k-2}$.

The particular choice of $k$ is not important when combined in conjunction with a parabolic smoothing argument such as the one to be used in the proof of Theorem 1 below. The lower bound on $\gamma$ could be improved slightly by introducing further technicalities (e.g., by Schauder- rather than just energy estimates in Condition 2D).

The general proof strategy follows ideas developed in [6, 7, 30, 29] who derive Bernstein-von Mises theorems in a variety of ‘non-conjugate’ settings. Infinite dimensional (‘functional’) such results have so far only been obtained in ‘direct’ iid models in [6, 8] and for non-linear inverse problems in [30, 35] not covering Gaussian priors. Here we extend the ‘semi-parametric’ proof in [29] for fixed one-dimensional functionals $\langle\theta,\psi\rangle$ to hold uniformly in $\psi$ belonging to a unit ball of $H^{k+2}$, thereby obtaining a Gaussian approximation in the infinite-dimensional $\|\cdot\|_{H^{-k-2}}$ space. We also prove our results in the stronger Wasserstein distance rather than just in a metric for weak convergence, allowing to obtain the asymptotics of the posterior mean and first moments.

For the proof we need to fix a few constants in advance: given $k$ from the theorem and $\zeta>d/2$ from Condition 2D, we take $\kappa,\bar{\kappa}$ satisfying

Further we take $\bar{\gamma}$ satisfying

and assume without loss of generality that Condition 2 holds for this $\bar{\gamma}$. [We note that if Condition 2 holds for some $\bar{\gamma}^{\prime}$, then it also holds for all $\bar{\gamma}\geq\bar{\gamma}^{\prime}$ because $H_{0}^{\bar{\gamma}}\subset H_{0}^{\bar{\gamma}^{\prime}}$.]

Define the span in $L^{2}(\Omega,W)$ of the ($dim(W)$-dimensional tensors of) trigonometric polynomials (15) as

Since the RKHS of the base prior from Condition 1 equals $\mathcal{H}=H^{\gamma}_{0}(\Omega)$ with equivalent norm $\|\cdot\|_{h^{\gamma}}$, and using (14), we have for every $J\in\mathbb{N}$ and some $c>0$ the standard approximation bounds

which we shall use repeatedly in the proofs.