Efficient Bayesian inference for large chaotic dynamical systems

Estimating parameters of chaotic geophysical models is challenging due to their inherent unpredictability. These models cannot be calibrated with standard least squares or filtering methods if observations are temporally sparse. Obvious remedies, such as averaging over temporal and spatial data to characterize the mean behavior, do not capture the subtleties of the underlying dynamics. We perform Bayesian inference of parameters in high-dimensional and computationally demanding chaotic dynamical systems by combining two approaches: (i) measuring model–data mismatch by comparing chaotic attractors and (ii) mitigating the computational cost of inference by using surrogate models. Specifically, we construct a likelihood function suited to chaotic models by evaluating a distribution over distances between points in the phase space; this distribution defines a summary statistic that depends on the geometry of the attractor, rather than on pointwise matching of trajectories. This statistic is computationally expensive to simulate, compounding the usual challenges of Bayesian computation with physical models. Thus, we develop an inexpensive surrogate for the log likelihood with the local approximation Markov chain Monte Carlo method, which in our simulations reduces the time required for accurate inference by orders of magnitude. We investigate the behavior of the resulting algorithm with two smaller-scale problems and then use a quasi-geostrophic model to demonstrate its large-scale application.

1 Introduction

Time evolution of many geophysical dynamical systems is chaotic. Chaoticity means that state of a system sufficiently far in the future cannot be predicted even if we know the dynamics and the initial conditions very precisely. Commonly used examples of chaotic systems include climate, weather, and the solar system.

A system being chaotic does not mean that it is random: the dynamics of models of chaotic systems are still determined by parameters, which may be either deterministic or random (Gelman et al., 2013). For example, Monte Carlo methods may be used to simulate future climate variability, but the distribution of possible climates will depend on the parameters of the climate model, and using the wrong model parameter distribution will result in potentially biased results with inaccurate uncertainties. For this reason, parameter estimation in chaotic models is an important problem for a range of geophysical applications. This paper focuses on Bayesian approaches to parameter inference in settings where (a) model dynamics are chaotic, and (b) sequential observations of the system are obtained so rarely that the model behavior has become unpredictable.

Parameters of a dynamical system model are most commonly inferred by minimizing a cost function that captures model–observation mismatch (Tarantola, 2005). In the Bayesian setting (Gelman et al., 2013), modeling this mismatch probabilistically yields a likelihood function, which enables maximum likelihood estimation or fully Bayesian inference. In fully Bayesian inference, the problem is further regularized by prescribing a prior distribution on the model state. In practice, Bayesian inference is often realized via Markov chain Monte Carlo (MCMC) methods (Gamerman, 1997; Robert and Casella, 2004).

This straightforward strategy – for instance, using the squared Euclidean distance between model outputs and data to construct a Gaussian likelihood – is, however, inadequate for chaotic models, where small changes in parameters, or even in the tolerances used for numerical solvers, can lead to arbitrarily large differences in model outputs (Rougier, 2013). Furthermore, modeling dynamical systems is often computationally very demanding, which makes sample generation time consuming. Since successful application of MCMC generally requires large numbers of model evaluations, performing Bayesian inference with MCMC is often not possible.

The present works combines two recent methods to tackle these problems due to model chaoticity and computational cost. Chaoticity is tamed by using correlation integral likelihood (CIL) (Haario et al., 2015), which is able to constrain the parameters of chaotic dynamical systems. We couple CIL with local approximation MCMC (LA-MCMC) (Conrad et al., 2016), which is a surrogate modeling technique that makes asymptotically exact posterior characterization feasible for computationally expensive models. We show how combining these methods enables a Bayesian approach to infer the parameters of chaotic high-dimensional models and quantify their uncertainties in situations previously discussed as intractable (Rougier, 2013). Moreover, we introduce several computational improvements to further enhance the applicability of the approach.

The CIL method is based on the concept of fractal dimension from mathematical physics, which broadly speaking characterizes the space-filling properties of the trajectory of a dynamical system. Earlier work (e.g., Cencini et al., 2010) describes a number of different approaches for estimating the fractal dimension. Our previous work extends this concept: instead of computing the fractal dimension of a single trajectory, a similar computation measures the distance between different model trajectories (Haario et al., 2015), based on which a specific summary statistic, called feature vector, is computed. The modification provides a normally distributed statistic of the data, which is sensitive to changes in the underlying attractor from which the data were sampled. Statistics that are sensitive to changes in the attractor yield likelihood functions that can better constrain the model parameters and therefore also result in more meaningful parameter posterior distributions.

The LA-MCMC method (Conrad et al., 2016, 2018) approximates the computationally expensive log-likelihood function using local polynomial regression. In this method, the MCMC sampler directly uses the approximation of the log likelihood to construct proposals and evaluate the Metropolis acceptance probability. Infrequently but regularly adding “full” likelihood evaluations to the point set used to construct the local polynomial regression continually improves the approximation, however. Expensive full likelihood evaluations are thus used only to construct the approximation or “surrogate” model. Conrad et al. (2016) show that, given an appropriate mechanism for triggering likelihood evaluations, the resulting Markov chain converges to the true posterior distribution while reducing the number of expensive likelihood evaluations (and hence forward model simulations) by orders of magnitude. Davis et al. (2020) show that LA-MCMC converges with approximately the expected $\mathrm{1}/\sqrt{T}$ error decay rate after a finite number of steps T, and Davis (2018) introduces a numerical parameter that ensures convergence even if only noisy estimates of the target density are available. This modification is useful for the chaotic systems studied here.

The rest of this paper is organized as follows. Section 2 reviews some additional background literature and related work and Sect. 3 describes the methodologies used in this work, including the CIL, the stochastic LA-MCMC algorithm, and the merging of these two approaches. Section 4 is dedicated to numerical experiments, where the CIL/LA-MCMC approach is applied to several, progressively more demanding examples. These examples are followed by a concluding discussion in Sect. 5.

2 Background and related work

Traditional parameter estimation methods, which directly utilize the model–observation mismatch, constrain the modeling to limited time intervals when the model is chaotic. This avoids the eventual divergence (chaotic behavior) of orbits that are initially close to each other. A classical example is variational data assimilation for weather prediction, where the initial states of the model are estimated using observational data and algorithms such as 4D-Var, after which a short-time forecast can be simulated (Asch et al., 2016).

Sequential data assimilation methods, such as the Kalman filter (KF) (Law et al., 2015), allow parameter estimation by recursively updating both the model state and the model parameters by conditioning them on observational data obtained over sufficiently short timescales. With methods such as state augmentation, model parameters can be updated as part of the filtering problem (Liu and West, 2001). Alternatively, the state values can be integrated out to obtain the marginal likelihood over the model parameters (Durbin and Koopman, 2012). Hakkarainen et al. (2012) use this filter likelihood approach to estimate parameters of chaotic systems. For models with strongly non-linear dynamics, ensemble filtering methods provide a useful alternative to the extended Kalman filter or variational methods; see Houtekamer and Zhang (2016) for a recent review of various ensemble variants.

Filtering-based approaches generally introduce additional tuning parameters, such as the length of the assimilation time window, the model error covariance matrix, and covariance inflation parameters. These choices have an impact on model parameter estimation and may introduce bias. Indeed, as discussed in Hakkarainen et al. (2013), changing the filtering method requires updating the parameters of the dynamical model. Alternatives to KF-based parameter estimation methods that do not require ensemble filtering include operational ensemble prediction systems (EPSs); for example, ensemble parameter calibration methods by Jarvinen et al. (2011) and Laine et al. (2011) have been applied to the Integrated Forecast System (IFS) weather models (Ollinaho et al., 2012, 2013, 2014) at the European Centre for Medium-Range Weather Forecasts (ECMWF). However, these approaches are heuristic and again limited to relatively short predictive windows.

Climate model parameters have in previous studies (e.g., Roeckner et al., 2003; ECMWF, 2013; Stevens et al., 2013) been calibrated by matching summary statistics of quantities of interest, such as top-of-atmosphere radiation, with the corresponding summary statistics from reanalysis data or output from competing models. The vast majority of these approaches produce only point estimates. A fully Bayesian parameter inversion was performed by Järvinen et al. (2010), who inferred closure parameters of a large-scale computationally intensive climate model, ECHAM5, using MCMC and several different summary statistics.

Computational limitations make applying algorithms such as MCMC challenging for weather and climate models. Generating even very short MCMC chains may require methods such as parallel sampling and early rejection for tractability (Solonen et al., 2012). Moreover, even if these computational challenges can be overcome, finding statistics that actually constrain the parameters is difficult, and inference results can be thus be inconclusive. The failure of the summary statistic approach in Järvinen et al. (2010) can be explained intuitively: the chosen statistics average out too much information and therefore fail to characterize the geometry of the underlying chaotic attractor in a meaningful way.

Several Monte Carlo methods have been presented to tackle expensive or intractable likelihoods; see, e.g., Luego et al. (2020) for a recent comprehensive literature review. Two notable such methods are approximate Bayesian computation (ABC) (Beaumont et al., 2002) and pseudo-marginal sampling (Andrieu and Roberts, 2009). The approach that is most closely related to the one presented in this paper is Bayesian inference using synthetic likelihoods, which was proposed as an alternative to ABC (Wood, 2010; Price et al., 2018). Recent work by Morzfeld et al. (2018) describes another feature vector approach for data assimilation. For more details and comparisons among these approaches, see the discussion below in Sect. 3.1.

In this work, we employ a different summary statistic, where the observations are considered as samples from the underlying attractor. Due to the nature of the summary statistic used in CIL, the observation time stamps are not explicitly used. This allows arbitrarily sparse observation time series, and consecutive observations may be farther than any window where the system remains predictable. To the best of our knowledge, parameter estimation in this setting has not been discussed in the literature. Another difference with the synthetic likelihood approach is that it involves regenerating data for computing the likelihood at every new model parameter value, which would be computationally unfeasible in our setting.

3 Methods

3.1 Correlation integral likelihood

We first construct a likelihood function that models the observations by comparing certain summary statistics of the observations to the corresponding statistics of a trajectory simulated from the chaotic model. As a source of statistics, we will choose the correlation integral, which depends on the fractal dimension of the chaotic attractor. Unlike other statistics – such as the ergodic averages of a trajectory – the correlation integral is able to constrain the parameters of a chaotic model (Haario et al., 2015).

Let us denote by

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}}=f(\mathit{u}\mathbf{,}\mathit{\theta }),\mathbf{u}(t=\mathrm{0})={\mathbf{u}}_{\mathrm{0}},\end{array}$$

a dynamical system with state u(t)∈ℝⁿ, initial state u₀∈ℝⁿ, and parameters θ∈ℝ^q. The time-discretized system, with time steps $t_i\in \mathit{\{}{t}_{\mathrm{1}},\mathrm{\dots },t_{\mathit{\tau }}\mathit{\}}$ denoting selected observation points, can be written as

$$\begin{array}{}\text{(2)}& {\mathbf{u}}_i\equiv \mathbf{u}\left(t_i\right)=F(t_i;{\mathbf{u}}_{\mathrm{0}},\mathit{\theta }).\end{array}$$

Either the full state u_i∈ℝⁿ or a subset ${\mathit{s}}_i\in {\mathbb{R}}^{d\le n}$ of the state components are observed. We will use $\mathbf{S}=\mathit{\{}{\mathit{s}}_{\mathrm{1}},\mathrm{\dots },{\mathit{s}}_{\mathit{\tau }}\mathit{\}}$ to denote a collection of these observables at successive times.

Using the model–observation mismatch at a collection of times t_i to constrain the value of the parameters θ is not suitable when the system (1) has chaotic dynamics, since the state vector values s_i are unpredictable after a finite time interval. Though long-time trajectories s(t) of chaotic systems are not predictable in the time domain, they do, however, represent samples from an underlying attractor in the phase space. The states are generated deterministically, but the model's chaotic nature allows us to interpret the states as samples from a particular θ-dependent distribution. Yet obvious choices for summary statistics T that depend on the observed states S, such as ergodic averages, ignore important aspects of the dynamics and are thus unable to constrain the model parameters. For example, the statistic $T\left(\mathbf{S}\right)=\frac{\mathrm{1}}{\mathit{\tau }}{\sum }_{i=\mathrm{1}}^{\mathit{\tau }}{\mathit{s}}_i$ is easy to compute and is normally distributed in the limit τ→∞ (under appropriate conditions), but this ergodic mean says very little about the shape of the chaotic attractor.

Instead, we need a summary statistic that retains information relevant for parameter estimation but still defines a computationally tractable likelihood. To this end, Haario et al. (2015) devised the CIL, which retains enough information about the attractor to constrain the model parameters. We first review the CIL and then discuss how to make evaluation of the likelihood tractable.

We will use the CIL to evaluate the “difference” between two chaotic attractors. For this purpose, we will first describe how to statistically characterize the geometry of a given attractor, given suitable observations S. In particular, constructing the CIL likelihood will require three steps: (i) computing distances between observables sampled from a given attractor; (ii) evaluating the empirical cumulative distribution function (ECDF) of these distances and deriving certain summary statistics T from the ECDF; and (iii) estimating the mean and covariance of T by repeating steps (i) and (ii).

Intuitively, the CIL thus interprets observations of a chaotic trajectory as samples from a fixed distribution over phase space. It allows the time between observations to be arbitrarily large – importantly, much longer than the system's non-chaotic prediction window.

Now we describe the CIL construction in detail. Suppose that we have collected a data set S comprising observations of the dynamical system of interest. Let S be split into n_epo different subsets called epochs. The epochs can, in principle, be any subsets of length N from the reference data set S. In this paper, however, we restrict the epochs to be time-consecutive intervals of N evenly spaced observations. Let ${\mathit{s}}^k=\mathit{\{}{\mathit{s}}_i^k{\mathit{\}}}_{i=\mathrm{1}}^N$ and ${\mathbf{s}}^l=\mathit{\{}{\mathit{s}}_j^l{\mathit{\}}}_{j=\mathrm{1}}^N$, with $\mathrm{1}\le k,l\le n_{\text{epo}}$ and k≠l, be two such disjoint epochs. The individual observable vectors ${\mathit{s}}_i^k\in {\mathbb{R}}^d$ and ${\mathit{s}}_j^l\in {\mathbb{R}}^d$ comprising each epoch come from time intervals $[t_{kN+\mathrm{1}},t_{(k+\mathrm{1})N}]$ and $[t_{lN+\mathrm{1}},t_{(l+\mathrm{1})N}]$, respectively. In other words, superscripts refer to different epochs and subscripts refer to the time points within those epochs. Haario et al. (2015) then define the modified correlation integral sum $C(R,N,{\mathit{s}}^k,{\mathit{s}}^l)$ by counting all pairs of observations that are less than a distance R>0 from each other:

$$\begin{array}{}\text{(3)}& {\displaystyle }C(R,N,{\mathit{s}}^k,{\mathit{s}}^l)={\displaystyle \frac{\mathrm{1}}{N^{\mathrm{2}}}}\sum _{i,j\le N}{\mathbf{1}}_{[\mathrm{0},R]}\left(∥{\mathit{s}}_i^k-{\mathit{s}}_j^l∥\right),\end{array}$$

where 1 denotes the indicator function and $\parallel \cdot \parallel$ is the Euclidean norm on ℝ^d. In the physics literature, evaluating Eq. (3) in the limit R→0, with k=l and i≠j, numerically approximates the fractal dimension of the attractor that produced s^k=s^l (Grassberger and Procaccia, 1983a, b). Here, we instead use Eq. (3) to characterize the distribution of distances between s^k and s^l at all relevant scales. We assume that the state space is bounded; therefore, an R₀ covering all pairwise distances in Eq. (3) exists. For a prescribed set of radii $R_m=R_{\mathrm{0}}{b}^{-m}$, with b>1 and $m=\mathrm{0},\mathrm{\dots },M$, Eq. (3) defines a discretization of the ECDF of the distances $\parallel {\mathit{s}}_i^k-{\mathit{s}}_j^l\parallel$, with discretization boundaries given by the numbers R_m.

Now we define $y_m^{k,l}=C(R_m,N,{\mathit{s}}^k,{\mathit{s}}^l)$ as components of a statistic $T({\mathit{s}}^k,{\mathit{s}}^l)={\mathit{y}}^{k,l}:=(y_{\mathrm{0}}^{k,l},\mathrm{\dots },y_M^{k,l})$. This statistic is also called the feature vector. According to Borovkova et al. (2001) and Neumeyer (2004), the vectors y^k,l are normally distributed, and the estimates of the mean and covariance converge at the rate $\sqrt{n_{\mathrm{epo}}}$ to their limit points. This is a generalization of the classical result of Donsker (1951), which applies to Independent and identically distributed samples from a scalar-valued distribution. We characterize this normal distribution by subsampling the full data set S. Specifically, we approximate the mean μ and covariance Σ of T by the sample mean and sample covariance of the set $\mathit{\{}{\mathit{y}}^{k,l}:\mathrm{1}\le k,l\le n_{\text{epo}},k\ne l\mathit{\}}$, evaluated for all $\frac{\mathrm{1}}{\mathrm{2}}{n}_{\mathrm{epo}}(n_{\mathrm{epo}}-\mathrm{1})$ pairs of epochs (s^k,s^l) using fixed values of R₀, b, M, and N.

The Gaussian distribution of T effectively characterizes the geometry of the attractor represented in the data set S. Now we wish to use this distribution to infer the parameters θ. Given a candidate parameter value $\stackrel{\mathrm{̃}}{\mathit{\theta }}$, we use the model to generate states ${\mathit{s}}^{*}\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)=\mathit{\left\{}{\mathit{s}}_i^{*}\right(\stackrel{\mathrm{̃}}{\mathit{\theta }}){\mathit{\}}}_{i=\mathrm{1}}^N$ for the length of a single epoch. We then evaluate the statistics $y_m^{k,*}=C(R_m,N,{\mathit{s}}^k,{\mathit{s}}^{*}(\stackrel{\mathrm{̃}}{\mathit{\theta }}\left)\right)$ as in Eq. (3), by computing the distances between elements of ${\mathit{s}}^{*}\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)$ and the states of an epoch s^k selected from the data S. Combining these statistics into a feature vector ${\mathit{y}}^{k,*}\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)=(y_m^{k,*}{)}_{m=\mathrm{0}}^M$, we can write a noisy estimate of the log-likelihood function:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\log p\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right|{\mathit{s}}^k)& =-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathit{y}}^{k,*}\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)-\mathit{\mu }\right)}^{\top }{\mathrm{\Sigma }}^{-\mathrm{1}}\left({\mathit{y}}^{k,*}\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)-\mathit{\mu }\right)\\ & +\text{constant}.\end{array}\end{array}$$

Comparing ${\mathit{s}}^{*}\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)$ with other epochs drawn from the data set S, however, will produce different realizations of the feature vector. We thus average the resulting log likelihoods over all epochs:

$$\begin{array}{}\text{(5)}& \log p\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\mathrm{|}\mathbf{S}\right)={\displaystyle \frac{\mathrm{1}}{n_{\text{epo}}}}\sum _{k=\mathrm{1}}^{n_{\text{epo}}}\log p\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\mathrm{|}{\mathit{s}}^k\right).\end{array}$$

This averaging, which involves evaluating Eq. (4) n_epo times, involves only new distance computations and is thus relatively cheap relative to time integration of the dynamical model.

Because the feature vectors ${\mathit{y}}^{k,*}$ are random for any finite N, and because the number of epochs n_epo is also finite, the log likelihood in Eq. (5) is necessarily random. It is then useful to view Eq. (5) as estimate of an underlying true log likelihood. We are therefore in a setting where cannot evaluate the unnormalized posterior density exactly; we only have access to a noisy approximation of it. Previous work (Springer et al., 2019) has demonstrated that derivative-free optimizers such as the differential evolution (DE) algorithm can successfully identify the posterior mode in this setting, yielding a point estimate of θ. In the fully Bayesian setting, one could characterize the posterior p(θ|S) using pseudo-marginal MCMC methods (Andrieu and Roberts, 2009) but at significant computational expense. Below, we will use a surrogate model constructed adaptively during MCMC sampling to reduce this computational burden.

Note that the CIL approach described above already reduces the computational cost of inference by only requiring simulation of the (potentially expensive) chaotic model for a single epoch. We compare each epoch of the data to the same single-epoch model output. Each of these comparisons results in an estimate of the log likelihood, which we then average over data epochs. A larger data set S can reduce the variance of this average, but does not require additional simulations of the dynamical model. Also, we do not require any knowledge about the initial conditions of the model; we omit an initial time interval before extracting ${\mathit{s}}^{*}\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)$ to ensure that the observed trajectory is on the chaotic attractor.

Moreover, the initial values are randomized for all simulations and sampling is started only after the model has integrated beyond the initial, predictable, time window. The independence of the sampled parameter posteriors from the initial values was verified both here and in earlier works by repeated experiments.

Our approach is broadly similar to the synthetic likelihood method (e.g., Wood, 2010; Price et al., 2018) but differs in two key respects: (i) we use a novel summary statistic that is able to characterize chaotic attractors, and (ii) we only need to evaluate the forward model for a single epoch. Comparatively, synthetic likelihoods typically use summary statistics such as auto-covariances at a given lag or regression coefficients. These methods also require long-time integration of the forward model for each candidate parameter value θ, rather than integration for only one epoch. Morzfeld et al. (2018) also discuss several ways of using feature vectors for inference in geophysics. A distinction of the present work is that we use an ECDF-based summary statistic that is provably Gaussian, and we perform extensive Bayesian analysis of the parameter posteriors via novel MCMC methods. These methods are described next.

3.2 Local approximation MCMC

Even with the developments described above, estimating the CIL at each candidate parameter value $\stackrel{\mathrm{̃}}{\mathit{\theta }}$ is computationally intensive. We thus use local approximation MCMC (LA-MCMC) (Conrad et al., 2016, 2018; Davis et al., 2020) – a surrogate modeling method that replaces many of these CIL evaluations with an inexpensive approximation. Replacing expensive density evaluations with a surrogate was first introduced by Sacks et al. (1989) and Kennedy and O'Hagan (2001). LA-MCMC extends these ideas by continually refining the surrogate during sampling, which guarantees convergence.

First introduced in Conrad et al. (2016), LA-MCMC builds local surrogate models for the log likelihood while simultaneously sampling the posterior. The surrogate is incrementally and infinitely refined during sampling and thus tailored to the problem – i.e., made more accurate in regions of high posterior probability. Specifically, the surrogate model is a local polynomial computed by fitting nearby evaluations of the “true” log likelihood. We emphasize that the approximation itself is not locally supported. At each point $\stackrel{\mathrm{̃}}{\mathit{\theta }}$, we locally construct a polynomial approximation, which globally defines a piecewise polynomial surrogate model. This is an important distinction because the piecewise polynomial approximation is not necessarily a probability density function. In fact, the surrogate function may not even be integrable. Despite this challenge, Davis et al. (2020) devise a refinement strategy that ensures convergence and bounds the error after a finite number of samples. In particular, Davis et al. (2020) shows that the error in the approximate Markov chain computed with the local surrogate model decays at approximately the expected $\mathrm{1}/\sqrt{T}$ rate, where T is the number of MCMC steps. Davis (2018) demonstrated that noisy estimates of the likelihood are sufficient to construct the surrogate model and still retain asymptotic convergence. Empirical studies (Conrad et al., 2016, 2018; Davis et al., 2020) on problems of moderate parameter dimension showed that the number of expensive likelihood evaluations per MCMC step can be reduced by orders of magnitude, with no discernable loss of accuracy in posterior expectations.

Here, we briefly summarize one step of the LA-MCMC construction and refer to Davis et al. (2020) for details. Each LA-MCMC step consists of four stages: (i) possibly refine the local polynomial approximation of the log likelihood, (ii) propose a new candidate MCMC state, (iii) compute the acceptance probability, and (iv) accept or reject the proposed state. The major distinction between this algorithm and standard Metropolis–Hastings MCMC is that the acceptance probability in stage (iii) is computed only using the approximation or surrogate model of the log likelihood, at both the current and proposed states. This introduces an error, relative to computation of the acceptance probability with exact likelihood evaluations, but stage (i) of the algorithm is designed to control and incrementally reduce this error at the appropriate rate.

“Refinement” in stage (i) consists of adding a computationally intensive log-likelihood evaluation at some parameter value θ_i, denoted by ℒ(θ_i), to the evaluated set $\mathit{\left\{}\right({\mathit{\theta }}_i,\mathcal{L}\left({\mathit{\theta }}_i\right)){\mathit{\}}}_{i=\mathrm{1}}^K$. These K pairs are used to construct the local approximation via a kernel-weighted local polynomial regression (Kohler, 2002). The values $\mathit{\{}{\mathit{\theta }}_i{\mathit{\}}}_{i=\mathrm{1}}^K$ are called “support points” in this paper. Details on the regression formulation are in Davis et al. (2020); Conrad et al. (2016). As the support points cover the regions of high posterior probability more densely, the accuracy of the local polynomial surrogate will increase. This error is well understood (Kohler, 2002; Conn et al., 2009) and, crucially, takes advantage of smoothness in the underlying true log-likelihood function. This smoothness ultimately allows the cardinality of the evaluated set to be much smaller than the number of MCMC steps.

Intuitively, if the surrogate converges to the true log likelihood, then the samples generated with LA-MCMC will (asymptotically) be drawn from the true posterior distribution. After any finite number of steps, however, the surrogate error introduces a bias into the sampling algorithm. The refinement strategy must therefore ensure that this bias is not the dominant source of error. At the same time, refinements must occur infrequently to ensure that LA-MCMC is computationally cheaper than using the true log likelihood. Davis et al. (2020) analyzes the trade-off between surrogate-induced bias and MCMC variance and proposes a rate-optimal refinement strategy. We use a similar algorithm, only adding an isotropic ℓ² penalty on the polynomial coefficients. More specifically, we rescale the variables so the max/min values are ±1 – often called coded units – then regularize the regression by a penalty parameter (in our cases, the value α=1 was found sufficient). This penalty term modifies the ordinary least squares problem into local ridge regression, which improves performance with noisy likelihoods.

Our examples use an adaptive proposal density Haario et al. (2006). This choice deviates slightly from the theory in Davis et al. (2020), which assumes a constant-in-time proposal density. However, this does not necessarily imply that adaptive or gradient-based methods will not converge. In particular, Conrad et al. (2016) show asymptotic convergence using an adaptive proposal density and Conrad et al. (2018) strengthen this result by showing that the Metropolis-adjusted Langevin algorithm, which is a gradient based MCMC method, is asymptotically exact when using a continually refined local polynomial approximation. These results require some additional assumptions about the target density's tail behavior and the stronger rate optimal result from Davis et al. (2020) has not been shown for such algorithms. However, in practice, we see that adaptive methods still work well in our applications. Exploring the theoretical implications of this is interesting and merits further discussion but is beyond the scope of this paper.

The parameters of the algorithm are fixed as given in Davis et al. (2020), for all the examples discussed here: (i) initial error threshold γ₀=1; (ii) error threshold decay rate γ₁=1; (iii) maximum poisedness constant $\overline{\mathrm{\Lambda }}=\mathrm{100}$; (iv) tail-correction parameter η=0 (no tail correction); (v) local polynomial degree p=2. The number of nearest neighbors k used to construct each local polynomial surrogate is chosen to be $k=k_{\mathrm{0}}+(K-k_{\mathrm{0}}{)}^{\mathrm{1}/\mathrm{3}}$ where $k_{\mathrm{0}}=\sqrt{q}D$, q is the dimension of the parameters θ∈ℝ^q, and D is the number of coefficients in the local polynomial approximation of total degree p=2, i.e., $D=(q+\mathrm{2})(q+\mathrm{1})/\mathrm{2}$. If we had k=D, the approximation would be an interpolant. Instead, we oversample by a factor $\sqrt{q}$, as suggested in Conrad et al. (2016), and allow k to grow slowly with the size K of the evaluated set as in Davis (2018). All these details together with example runs can be found in the MATLAB implementation available in the Supplement.

4 Numerical experiments

This section contains numerical experiments to illustrate the methods introduced in the previous sections. As a large-scale example, we characterize the posterior distribution of parameters in the two-layer quasi-geostrophic (QG) model. The computations needed to characterize the posterior distribution with standard MCMC methods in this example would be prohibitive without massive computational resources and are therefore omitted. In contrast, we will show that the LA-MCMC method is able to simulate from the parameter posterior distribution.

Before presenting this example, we first demonstrate that the posteriors produced by LA-MCMC agree with those obtained via exact MCMC sampling methods in cases where the latter are computationally tractable using two examples: the classical Lorenz 63 system and the higher-dimensional Kuramoto–Sivashinsky (KS) model. In both examples, we quantify the computational savings due to LA-MCMC, and in the second we introduce additional ways to speed up computation using parallel (GPU) integration.

Let Δ_t denote the time difference between consecutive observations; one epoch thus contains the times in the interval $[iN{\mathrm{\Delta }}_t,(i+\mathrm{1}\left)N{\mathrm{\Delta }}_t\right)$. The number of data points in one epoch N varies between 1000 and 2000, depending on the example. The training set S consists of a collection of n_epo such intervals. In all the examples, we choose Δ_t to be relatively large, beyond the predictable window. This is more for demonstration purposes than a necessity; the background theory from U statistics allows the subsequent state vectors to be weakly dependent. Numerically, a set of observations that is chosen too densely results in the χ² test failing, and for this reason we recommended to always check for normality before starting the parameter estimation.

For numerical tests, one can either use one long time series or integrate a shorter time interval several times using different initial values to create the training set for the likelihood. For these experiments, the latter method was used with n_epo=64, yielding $\begin{array}{c}{n}_{\mathrm{epo}}\\ \mathrm{2}\end{array}=\mathrm{2016}$ different pairs (s^k,s^l), each of which resulted in an ECDF constructed from N² pairwise distances. According to tests performed while calibrating the algorithm, these values of N and n_epo are sufficient to obtain robust posterior estimates. With less data, the parameter posteriors will be less precise.

The range of the bin radii $R_m,m=\mathrm{0},\mathrm{\dots },M$ is selected by examining the distances within the training set, keeping in mind that a positive variance is needed for every bin to avoid a singular covariance matrix. So the largest radius R₀ can be obtained from

$$\begin{array}{}\text{(6)}& R_{\mathrm{0}}=\underset{k\ne l}{min}\left\{\underset{i,j}{max}∥s_i^k-s_j^l∥\right\}\end{array}$$

over the disjoint subsets of the samples s^k and s^l of length N. The smallest radius is selected by requiring that for all the possible pairs (s^k,s^l), it holds that ${\mathcal{B}}_{R_M}\left({\mathit{s}}_i^k\right)\cap {\mathit{s}}^l\ne \mathrm{\varnothing }$, where ${\mathcal{B}}_{R_M}\left({\mathit{s}}_i^k\right)$ is the ball of radius R_M centered at ${\mathit{s}}_i^k$. That is,

$$\begin{array}{}\text{(7)}& R_M=\underset{k\ne l}{max}\left\{\underset{i,j}{min}∥s_i^k-s_j^l∥\right\}.\end{array}$$

The base value b is obtained by $R_M=R_{\mathrm{0}}{b}^{-M}$, and using this value we fix all the other radii R_m.

As always with histograms, the number of bins M must be selected first. Too small an M loses information, while too large values yield noisy histograms, and this noisiness can be seen also in the ECDFs. However, numerical experiments show that the final results – the parameter posteriors – are not too sensitive to the specific value of M. For instance, for the Lorenz 63 case below, the range of M was varied between 5 and 40, and only a minor decrease of the size of the parameter posteriors was noticed for increasing M. Any slight increase of accuracy comes with a computational cost: higher values of M increase the stochasticity of the likelihood evaluations, which leads to smaller acceptance rates in the MCMC sampling, e.g., from 0.36 to 0.17 to 0.03 for $M=\mathrm{5},\mathrm{15},\mathrm{40}$, respectively, when using the standard adaptive Metropolis sampler. In the examples presented in this section, the length M of the feature vector is fixed to 14 for the Lorenz 63 model and 32 for the higher-dimensional KS and QG models.

To balance the possibly different magnitudes of the components of the state vector, each component is scaled and shifted to the range $[-\mathrm{1},\mathrm{1}]$ before computing the distances. While this scaling could also be performed in other ways, this method worked well in practice for the models considered. The normality of the ensemble of feature vectors is ascertained by comparing the histograms of the quadratic forms in Eq. (4) visually to the appropriate χ² distribution.

In all the three experiments, we create MCMC chains of length 10⁵. However, due to the use of the LA-MCMC approach, the number of full forward model evaluations is much lower, around 1000 or less; we will report these values more specifically below.

The Lorenz 63 model was integrated with a standard Runge–Kutta solver. The numerical solution of the KS-model is based on our in-house fast Fourier transform (FFT)-based solver, which runs on the GPU side and is built around Nvidia compute unified device architecture (CUDA) toolchain and cuFFT library (which is a part of the CUDA ecosystem). The quasi-geostrophic model employs semi-Lagrangian solver and runs entirely on CPU, but the code has been significantly optimized with performance-critical parts, such as advection operator, compiled using an Intel single program compiler (ISPC) with support of Advanced Vector Extensions 2 (AVX2) vectorization.

4.1 Lorenz 63

We use the classical three-dimensional Lorenz 63 system (Lorenz, 1963) as a simple first example to demonstrate how LA-MCMC can be successfully paired with the CIL and the adaptive Metropolis (AM) algorithm (Haario et al., 2001, 2006) to obtain the posterior distribution for chaotic systems at a greatly reduced computational cost, compared to AM without the local approximation. The time evolution of the state vector $\mathit{s}=(X,Y,Z)$ is given by

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\dot{X}& =\mathit{\sigma }(Y-X),\\ \dot{Y}& =X(\mathit{\rho }-Z)-Y,\\ \dot{Z}& =XY-\mathit{\beta }Z.\end{array}\end{array}$$

This system of equations is often said to describe an extreme simplification of a weather model.

The reference data were generated with parameter values σ=10, ρ=28, and $\mathit{\beta }=\frac{\mathrm{8}}{\mathrm{3}}$ by performing n_epo=64 distinct model simulations, with observations made at 2000 evenly distributed times between [10,20 000]. These observations were perturbed with 5 % multiplicative Gaussian noise. The length of the predictable time window is roughly 7, which is less than the time between consecutive observations. The parameters of the CIL method were obtained as described at the start of Sect. 4, with values M=14, R₀=2.85, and b=1.51.

The set of vectors $\mathit{\left\{}{\mathit{y}}^{k,l}\right|k,l\le n_{\text{epo}}\mathit{\}}$ is shown in Fig. 1 in the log–log scale. The figure shows how the variability of these vectors is quite small. Figure 1 validates the normality assumption for feature vectors.

Figure 1 Left: for all combinations of k and l, the feature vectors for Lorenz 63 and Kuramoto–Sivashinsky, along with concatenated feature vectors for the quasi-geostrophic system are shown. Right: normality check; the χ² density function versus the histograms of the respective negative log-likelihood (NLL) values; see Eq. (4).

Figure 2 Two-dimensional posterior marginal distributions of the parameters of the Lorenz 63 model obtained with LA-MCMC and AM.

Pairwise two-dimensional marginals of the parameter posterior are shown in Fig. 2, both from sampling the posterior with full forward model simulations (AM) and with using the surrogate sampling approach for generating the chain. These two posteriors are almost perfectly superimposed. Indeed, the difference is at the same level as that between repetitions of the standard AM sampling alone.

To get an idea of the computational savings achieved with LA-MCMC, the computation of the MCMC chains of length 10⁵ was repeated 10 times. The cumulative number of full likelihood evaluations is presented in Fig. 3. At the end of the chains, the number of full likelihood evaluations varied between 955 and 1016. Thus, by using LA-MCMC in this setting, remarkable computational savings of up to 2 orders of magnitude are achieved.

Figure 3 Comparison of the cumulative number of full likelihood evaluations while using AM (black line) and LA-MCMC (colored lines). Every colored line correspond to a different chain obtained with LA-MCMC by using the same likelihood.

4.2 The Kuramoto–Sivashinsky model

The second example is the 256-dimensional Kuramoto–Sivashinsky (KS) partial differential equation (PDE) system. The purpose of this example is to introduce ways to improve the computational efficiency by a piecewise parallel integration over the time interval of given data. Also, we demonstrate how decreasing the number of observed components impacts the accuracy of parameter estimation. Even though the posterior evaluation proves to be relatively expensive, direct verification of the results with those obtained by using standard adaptive MCMC is still possible. The Kuramoto–Sivashinsky model is given by the fourth-order PDE:

$$\begin{array}{}\text{(9)}& s_t=-s{s}_x-{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }}}{s}_{xx}-\mathit{\gamma }{s}_{xxxx},\end{array}$$

where $s=s(x,t)$ is a real function of x∈ℝ and $t\in {\mathbb{R}}_{+}$. In addition, it is assumed that s is spatially periodical with period of L, i.e., $s(x+L,t)=s(x,t)$. This experiment uses the parametrization from (Yiorgos Smyrlis, 1996) that maps the spatial domain $[-\frac{L}{\mathrm{2}},\frac{L}{\mathrm{2}}]$ to $[-\mathit{\pi },\mathit{\pi }]$ by setting $\stackrel{\mathrm{̃}}{x}=\frac{\mathrm{2}\mathit{\pi }}{L}x$ and $\stackrel{\mathrm{̃}}{t}={\left(\frac{\mathrm{2}\mathit{\pi }}{L}\right)}^{\mathrm{2}}t$. With L=100, the true value of parameter γ is $(\mathit{\pi }/\mathrm{50}{)}^{\mathrm{2}}\approx \mathrm{0.0039}$, and the true value of η becomes $\frac{\mathrm{1}}{\mathrm{2}}$. These two parameters are the ones that are then estimated with the LA-MCMC method. This system was derived by Kuramoto and Yamada (1976) and Kuramoto (1978) as a model for phase turbulence in reaction–diffusion systems. Sivashinsky (1977) used the same system for describing instabilities of laminar flames.

Assume that the solution for this problem can be represented by a truncation of the Fourier series

$$\begin{array}{}\text{(10)}& s(x,t)=\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\left[A_j\left(t\right)\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L}}jx\right)+B_j\left(t\right)\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L}}jx\right)\right].\end{array}$$

Using this form reduces Eq. (9) to a system of ordinary differential equations for the unknown coefficients A_j(t) and B_j(t),

$$\begin{array}{}\text{(11)}& {\displaystyle }{\dot{A}}_j\left(t\right)={\mathit{\alpha }}_{\mathrm{1}}{j}^{\mathrm{2}}{A}_j\left(t\right)+{\mathit{\alpha }}_{\mathrm{2}}{j}^{\mathrm{4}}{A}_j\left(t\right)+F_{\mathrm{1}}\left(A\right(t\left)\right)\text{(12)}& {\displaystyle }{\dot{B}}_j\left(t\right)={\mathit{\beta }}_{\mathrm{1}}{j}^{\mathrm{2}}{B}_j\left(t\right)+{\mathit{\beta }}_{\mathrm{2}}{j}^{\mathrm{4}}{B}_j\left(t\right)+F_{\mathrm{2}}\left(B\right(t\left)\right),\end{array}$$

where the terms F₁(⋅) and F₂(⋅) are polynomials of the vectors A and B. For details, see Huttunen et al. (2018). The solution can be effectively computed on graphics processors in parallel, and if computational resources allow, several instances of Eq. (9) can be solved in parallel. Even on fast consumer-level laptops, several thousand simulations can be performed in parallel when the discretization of the x dimension contains around 500 points.

A total of 64 epochs of the 256-dimensional KS model are integrated over the time interval [0,150 000], and as in the case of the Lorenz 63 model, the initial predictable time window is discarded and 1024 equidistant measurements from [500,150 000] are selected, with Δ_t≈146. The parameters used for the CIL method were R₀=1801.7, M=32, and b=1.025.

The time needed to integrate the model up to t=150 000 is approximately 103 s with the Nvidia 1070 GPU, implying that generating an MCMC chain with 100 000 samples with standard MCMC algorithms would take almost 4 months. The use of LA-MCMC alone again shortens the time needed by a factor of 100 to around 28 h. However, the calculations can yet be considerably enhanced by parallel computing. In practice, this translates the problem of generating a candidate trajectory of length 150 000 into generating observations from several shorter time intervals. In our example, an efficient division is to perform 128 parallel calculations each of length 4500, with randomized initial values close to the values selected from the training set. Discarding the predictable interval [0,500] and taking eight observations at intervals of 500 yields the same number (1024) of observations as in the initial setting. While the total integration time increases, this reduces the wall-clock time needed for computation of a single candidate simulation from 103 to 2.5 s. The full MCMC chain can be then be generated in 70 h without the surrogate model and in 42 min using LA-MCMC.

Table 1 Parameter values of the four parameter vectors used in the forward KS model simulation examples in Fig. 5. The parameter vectors in the first column labeled 1 are the true parameters, and the second one resides inside the posterior. The last two are outside the posterior. These parameters correspond to points shown in the posterior distribution shown in Fig. 4.

Parameter posterior distributions from the KS system, produced with MCMC both with and without the local approximation surrogate, are shown in Fig. 4. Repeating the calculations several times yielded no meaningful differences in the results. In this experiment, the number of forward model evaluations LA-MCMC needed for generating a chain of length 100 000 was in the range [1131,1221].

Model trajectories from simulations with four different parameter vectors are shown in Fig. 5. These parameter values were (1) the “true” value which was used to generate training data, (2) another parameter from inside the posterior distribution, and (3–4) two other parameters from outside the posterior distribution. These parameters are also shown in Fig. 4. Visually inspecting the outputs, cases 1–3 look similar, while results using parameter vector 4, furthest away from the posterior, are markedly different. Even though the third parameter vector is outside the posterior, the resulting trajectory is not easily distinguishable from cases 1 and 2, indicating that the CIL method differentiates between the trajectories more efficiently.

Figure 4 Posterior distribution of the parameters of the KS system. The parameter values are shown in Table 1, while examples of the respective integrated trajectories are given in Fig. 5.

Figure 5 Example model trajectories from the KS system. Panel (1) shows simulation using the true parameters, the parameters used for (4) are inside the posterior distribution, and (2) and (3) are generated from simulations with parameters outside the posterior distribution, shown in Fig. 4. The values of the parameter vectors 1, 2, 3, and 4 are given in Table 1. The y axis shows the 256-dimensional state vector, and the x axis the time evolution of the system.

Additional experiments were performed to evaluate the stability of the method when not all of the model states were observed. Keeping the setup otherwise fixed, the number of elements of the state vectors observed was reduced from the full 256 step by step to 128, 64, and 32. The resulting MCMC chains are presented in Fig. 6, and as expected, when less is observed, the size of the posterior distribution grows.

Figure 6 Comparison between the KS system's posterior distribution in cases where all or only a part of the states are observed.

4.3 The quasi-geostrophic model

The methodology is here applied to a computationally intensive model, where a brute-force parameter posterior estimation would be too time consuming. We employ the well-known quasi-geostrophic model (Fandry and Leslie, 1984; Pedlosky, 1987) using a dense grid to achieve complex chaotic dynamics in high dimensions. The wall-clock time for one long-time forward model simulation is roughly 10 min, so a naïve calculation of a posterior sample of size 100 000 would take around 2 years. We demonstrate how the application of the methods verified in the two previous examples reduces this time to a few hours.

The QG model approximates the behavior on a latitudinal “stripe” at two given atmospheric heights, projected onto a two-layered cylinder. The model geometry implies periodic boundary conditions, seamlessly stitching together the extreme eastern and western parts of the rectangular spatial domain with coordinates x and y. For the northern and southern edges, user-specified time-independent Dirichlet boundary conditions are used. In addition to these conditions and the topographic constraints, the model parameters include the mean thicknesses of the two interacting atmospheric layers, denoted by H₁ and H₂. The QG model also accounts for the Coriolis force. An example of the two-layer geometry is presented in Fig. 7.

Figure 7 An example of the layer structure of the two-layer quasi-geostrophic model. The terms U₁ and U₂ denote mean zonal flows, respectively, in the top and the bottom layer.

In a non-dimensional form, the QG system can be written as

$$\begin{array}{}\text{(13)}& {\displaystyle }{q}_{\mathrm{1}}& {\displaystyle }=\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{1}}-F_{\mathrm{1}}({\mathit{\psi }}_{\mathrm{1}}-{\mathit{\psi }}_{\mathrm{2}})+\mathit{\beta }y,\text{(14)}& {\displaystyle }{q}_{\mathrm{2}}& {\displaystyle }=\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}-F_{\mathrm{2}}({\mathit{\psi }}_{\mathrm{2}}-{\mathit{\psi }}_{\mathrm{1}})+\mathit{\beta }y+R_s,\end{array}$$

where q_i are potential vorticities, and ψ_i are stream functions with indexes $i=\mathrm{1},\mathrm{2}$ for the upper and the lower layers, respectively. Both the q_i and ψ_i are functions of time t and spatial coordinates x and y. The coefficients $F_i=\frac{f_{\mathrm{0}}^{\mathrm{2}}{L}^{\mathrm{2}}}{\stackrel{\mathrm{´}}{g}{H}_i}$ control how much the model layers interact, and $\mathit{\beta }={\mathit{\beta }}_{\mathrm{0}}L/U$ gives a nondimensional version of β₀, the northward gradient of the Coriolis force that gives rise to faster cyclonic flows closer to the poles. The Coriolis parameter is given by f₀=2Wsin (ℓ), where W is the angular speed of Earth and ℓ is the latitude of interest. L and U give the length and speed scales, respectively, and $\stackrel{\mathrm{´}}{g}$ is a gravity constant. Finally, $R_s(x,y)=\frac{S(x,y)}{\mathit{\eta }{H}_{\mathrm{2}}}$ defines the topography for the lower layer, where $\mathit{\eta }=\frac{U}{f_{\mathrm{0}}L}$ is the Rossby number of the system. For further details, see Fandry and Leslie (1984) and Pedlosky (1987).

Figure 8 An example of the 6050-dimensional state of the quasi-geostrophic model. The contour lines for both the stream function and potential vorticity are shown for both layers. Note the cylindrical boundary conditions.

It is assumed that the motion determined by the model is geostrophic, essentially meaning that potential vorticity of the flow is preserved on both layers:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial q_i}{\partial t}}+u_i{\displaystyle \frac{\partial q_i}{\partial x}}+v_i{\displaystyle \frac{\partial q_i}{\partial y}}=\mathrm{0}.\end{array}$$

Here, u_i and v_i are velocity fields, which are functions of both space and time. They are obtained from the stream functions ψ_i via

$$\begin{array}{}\text{(16)}& u_i=-{\displaystyle \frac{\partial {\mathit{\psi }}_i}{\partial y}},v_i={\displaystyle \frac{\partial {\mathit{\psi }}_i}{\partial x}}.\end{array}$$

Equations (13)–(16) define the spatiotemporal evolution of the quantities q_i,ψ_i, $i=\mathrm{1},\mathrm{2}$.

The numerical integration of this system is carried out using the semi-Lagrangian scheme, where the potential vorticities q_i are computed according to Eq. (15) for given velocities u_i and v_i. With these q_i the stream functions can then be obtained from Eqs. (13) and (14) with a two-stage finite difference scheme.

Finally, the velocity field is updated by Eq. (16) for the next iteration round.

For estimating model parameters from synthetic data, a reference data set is created with 64 epochs each containing N=1000 observations. These data are sampled from the model trajectory with Δ_t=8 (where a time step of length 1 corresponds to 6 h) in the time interval [192,8192], that amounts to a long-range integration of roughly 5–6 years of a climate model. The spatial domain is discretized into a 55×55 grid, which results in consistent chaotic behavior and more complex dynamics than with the often-used 20×40 grid. This is reflected in higher variability in the feature vectors, as seen in the Fig. 1. A snapshot of the 6050-dimensional trajectory of the QG system is displayed in Fig. 8.

The model state is characterized by two distinct fields, the vorticities and stream functions, that naturally are dependent on each other. As shown in (Haario et al., 2015), it is useful to construct separate feature vectors to characterize the dynamics in such situations. For this reason, two separate feature vectors are constructed – one for the potential vorticity on both layers and the other for the stream function.

The Gaussian likelihood of the state is created by stacking these two feature vectors one after another.

The normality of the resulting 2(M+1)-dimensional vector may again be verified as shown in Fig. 1. The number of bins was set to 32, leading to parameter values R₀=55 and b=1.075 for potential vorticity, and R₀=31, and b=1.046 for the stream function.

For parameter estimation, inferring the layer heights from synthetic data is considered. The reference data set with n_epo=64 integrations is produced using the values H₁=5500 and H₂=4500. A single forward model evaluation takes 10 min on a fast laptop, and therefore generating MCMC chains of length 10⁵ with brute force would take around 2 years to run. As previously, using LA-MCMC again reduces the computation time by a factor of 100.

In the experiments performed, the number of forward model evaluations needed was ranging in the interval [682,762], which translates to around 1 week of computing time. As verified with Kuramoto–Sivashinsky example, the forward model integration can be split to segments computed in parallel, which reduced time required to generate data for computing the likelihood further with a factor around 50, corresponding to around 3 h for generating the MCMC chain. The pairwise distances for generating the feature vectors were computed on a GPU, and therefore the required computation time for doing this was negligible compared to the model integration time. The posterior distribution of the two parameters is presented in Fig. 9.

Figure 9 The clearly non-Gaussian posterior distribution of the H₁ and H₂ parameters of the quasi-geostrophic system shows how these parameters anticorrelate with each other.

5 Conclusions and future work

Bayesian parameter estimation with computationally demanding computer models is highly non-trivial. The associated computational challenges often become insurmountable when the model dynamics are chaotic. In this work, we showed it is possible to overcome these challenges by combining the CIL with an MCMC method based on local surrogates of the log-likelihood function (LA-MCMC). The CIL captures changes in the geometry of the underlying attractor of the chaotic system, while local approximation MCMC makes generating long MCMC chains based on this likelihood tractable, with computational savings of roughly 2 orders of magnitude, as shown in Table 2. Our methods were verified by sampling the parameter posteriors of the Lorenz 63 and the Kuramoto–Sivashinsky models, where an (expensive) comparison to exact MCMC with the CIL was still feasible. Then we applied our approach to the quasi-geostrophic model with a deliberately extended grid size. Without CIL, parameter estimation would not have been possible with chaotic models such as these; without LA-MCMC, the generation of long MCMC and sufficiently accurate chains for the higher-resolution QG model parameters would have been computationally intractable. We note that the computational demands of the QG model already get quite close to those of weather models at coarse resolutions. We believe that the approach developed here can provide ways to solve problems such as the climate model closure parameter estimation investigated in Järvinen et al. (2010) or long-time assimilation problems with uncertain model parameters, discussed in Rougier (2013) as unsolved and intractable.

Table 2 Summary of results. This table shows the speed-up due to the CIL/LA-MCMC combination. Since running the quasi-geostrophic model 100 000 times was not possible, the nominal length of the MCMC chain and the speed-up due to LA-MCMC are reported in parentheses in the last column. The numbers of forward model evaluations with LA-MCMC (second row) are rough averages over several MCMC simulations.

There are many potential directions for extension of this work. First, it should be feasible to run parallel LA-MCMC chains that share model evaluations in a single evaluated set; doing so can accelerate the construction of accurate local surrogate models, as demonstrated in Conrad et al. (2018), and is a useful way of harnessing parallel computational resources within surrogate-based MCMC. Extending this approach to higher-dimensional parameters is also of interest. While LA-MCMC has been successfully applied to chains of dimension up to q=12 (Conrad et al., 2018), future work should explore sparsity and other truncations of the local polynomial approximation to improve scaling with dimension. From the CIL perspective, calibrating more complex models, such as weather models, often requires choosing the part of the state vector from which the feature vectors are computed. While computing the likelihood from the full high-dimensional state is computationally feasible, Haario et al. (2015) showed that carefully choosing a subset of the state for the feature vectors performs better. Also, the epochs may need to be chosen sufficiently long to include potential rare events, so that changes in rare event patterns can be identified. This, naturally, will increase the computational cost if one wants to be confident in the inclusion of such events.

While answering these questions will require further work, we believe the research presented in this paper provides a promising and reasonable step towards estimating parameters in the context of expensive operational models.