Ensemble Kalman methods for high-dimensional hierarchical dynamic space-time models
Matthias Katzfuss, Jonathan R. Stroud, Christopher K. Wikle
TL;DR
This paper introduces a novel class of filtering and smoothing algorithms that combine ensemble Kalman methods with state-space models to enable efficient inference in high-dimensional, nonlinear, non-Gaussian spatio-temporal systems, outperforming existing methods.
Contribution
It develops ensemble Kalman-based algorithms for high-dimensional hierarchical dynamic space-time models, integrating Bayesian inference with dimension reduction techniques.
Findings
Outperforms existing ensemble Kalman, particle filters, and particle MCMC methods.
Effective in high-dimensional, nonlinear, non-Gaussian scenarios.
Demonstrated success with real cloud motion data.
Abstract
We propose a new class of filtering and smoothing methods for inference in high-dimensional, nonlinear, non-Gaussian, spatio-temporal state-space models. The main idea is to combine the ensemble Kalman filter and smoother, developed in the geophysics literature, with state-space algorithms from the statistics literature. Our algorithms address a variety of estimation scenarios, including on-line and off-line state and parameter estimation. We take a Bayesian perspective, for which the goal is to generate samples from the joint posterior distribution of states and parameters. The key benefit of our approach is the use of ensemble Kalman methods for dimension reduction, which allows inference for high-dimensional state vectors. We compare our methods to existing ones, including ensemble Kalman filters, particle filters, and particle MCMC. Using a real data example of cloud motion and data…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13| Filtering/smoothing methods | ||||
| Target distrib. | Dim. of | Cond. Gaussian | Particle | Ensemble |
| , | ( known) | Kalman filter | Particle filter | EnKF (Alg. 2.1.3) |
| ( known) | FFBS | Particle smoother | EnKS (Alg. 2.3) | |
| , | high | RWFFBS | Online MCMC | GEnKF (Alg. 4.1.1) |
| low | RBPF | EARBPF | PEnKF (Alg. 4.1.2) | |
| high | Gibbs with FFBS | Particle Gibbs | GEnKS (Alg. 4.2.1) | |
| low | MH with KF | PMMH | MHEnKS/ | |
| PEnKS (Alg. 4.2.2) | ||||
| Heavy-tailed | Rainfall () | Rainfall ( unknown) | ||||
|---|---|---|---|---|---|---|
| RMSPE | CRPS | RMSPE | CRPS | RMSPE | CRPS | |
| exact | 0.167 | 0.093 | 0.435 | 0.144 | 0.422 | 0.131 |
| GEnKF | 0.185 | 0.103 | 0.452 | 0.152 | 0.939 | 0.220 |
| EnKF | 0.288 | 0.149 | 26.714 | 7.427 | 100 | 100 |
| PF | 0.767 | 0.586 | 2.031 | 1.042 | 3.912 | 1.491 |
| Huber | 0.643 | 0.414 | ||||
| Predictions | EMD to true filtering distributions | |||||||
| MSPE | CRPS | |||||||
| truth | 0.73 | 0.25 | ||||||
| PEnKF | 0.75 | 0.25 | 0.25 | 0.16 | 0.19 | 0.70 | 0.23 | 0.32 |
| EARBPF | 1.26 | 0.33 | 0.27 | 0.25 | 0.24 | 0.71 | 0.32 | 0.37 |
| Parameter | State | |||
|---|---|---|---|---|
| MSPE | CRPS | MSPE | CRPS | |
| GEnKS | 0.002 | 0.024 | 0.710 | 0.478 |
| EnKS+SA | 100 | 100 | 0.914 | 0.540 |
| Particle Gibbs | 0.111 | 0.262 | 12.380 | 2.495 |
| Prior | 0.042 | 0.118 | ||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Ensemble Kalman methods for high-dimensional hierarchical dynamic space-time models
Matthias Katzfuss Texas A&M University ([email protected])
Jonathan R. Stroud Georgetown University
Christopher K. Wikle University of Missouri
Abstract
We propose a new class of filtering and smoothing methods for inference in high-dimensional, nonlinear, non-Gaussian, spatio-temporal state-space models. The main idea is to combine the ensemble Kalman filter and smoother, developed in the geophysics literature, with state-space algorithms from the statistics literature. Our algorithms address a variety of estimation scenarios, including on-line and off-line state and parameter estimation. We take a Bayesian perspective, for which the goal is to generate samples from the joint posterior distribution of states and parameters. The key benefit of our approach is the use of ensemble Kalman methods for dimension reduction, which allows inference for high-dimensional state vectors. We compare our methods to existing ones, including ensemble Kalman filters, particle filters, and particle MCMC. Using a real data example of cloud motion and data simulated under a number of nonlinear and non-Gaussian scenarios, we show that our approaches outperform these existing methods.
Keywords: data assimilation; geoscience applications; Gibbs sampler; particle filter; spatio-temporal statistics; state-space models
1 Introduction
1.1 Inference in state-space models
State-space models (SSMs) are general representations of systems observed over time. A SSM describes the temporal evolution of the system state, and the relationship of the state to observations (e.g., West and Harrison, , 1997; Shumway and Stoffer, , 2006). Here, we analyze spatio-temporal SSMs describing the evolution of one or more discretized spatial fields over time. These models arise, for example, in geoscience and biomedical applications. They are often nonlinear, non-Gaussian, and have very high-dimensional states. Standard statistical approaches such as Markov chain Monte Carlo (MCMC) and particle filters break down for this class of models. Thus, new methods are needed for accurate inference. We propose a number of new filtering and smoothing algorithms for Bayesian inference in these high-dimensional SSMs.
There is an enormous statistics and engineering literature on state and parameter estimation for low-dimensional SSMs. For linear, Gaussian models, the Kalman filter (Kalman, , 1960) provides closed-form posterior distributions for the states conditional on the parameters. For conditionally Gaussian models, inference is often carried out using MCMC methods based on the forward-filter backward-sampling algorithm (Carter and Kohn, , 1994; Frühwirth-Schnatter, , 1994) or the mixture Kalman filter (Chen and Liu, , 2000). MCMC methods have also been developed for more general SSMs, including nonlinear, non-Gaussian models (Carlin et al., , 1992) and non-Gaussian measurement models (Shephard and Pitt, , 1997; Gamerman, , 1998). However, these methods do not scale to high-dimensional states because their computational cost is cubic in the state dimension.
Another popular approach for nonlinear, non-Gaussian SSMs is sequential Monte Carlo (SMC) methods, also known as particle filters (e.g., Doucet et al., , 2001). SMC methods approximate the state distribution at each time by a weighted set of samples or particles. These particles are propagated forward through time according to the state evolution. They are then updated based on new data using a reweighting or resampling step (e.g., Gordon et al., , 1993; Liu and Chen, , 1998). Although these methods are exact with an infinite number of particles, it is well known that they suffer from weight degeneracy (i.e., all weights but one become essentially zero) in high-dimensional problems (e.g., Snyder et al., , 2008). Thus, particle filters do not scale to high dimensions.
In contrast to the statistical approaches above, approximate state-estimation methods in the geoscience literature scale well to high-dimensional problems. In this paper, we focus on a geophysical data assimilation method called the ensemble Kalman filter (EnKF; Evensen, , 1994); see Houtekamer and Zhang, (2016) and Katzfuss et al., (2016) for recent reviews, and Wikle and Berliner, (2007) for a general introduction to data assimilation. The EnKF also uses an ensemble representation that is propagated forward like the particle filter, but instead of reweighting, the EnKF updates the ensemble using a linear “shift” that approximates the best linear update. While the EnKF is based on the assumption of a linear, Gaussian model, it provides accurate inference for many nonlinear models as well. This has allowed the EnKF to be successfully applied to nonlinear spatio-temporal SSMs in millions of dimensions (e.g., Szunyogh et al., , 2008; Whitaker et al., , 2008; Houtekamer et al., , 2014).
The EnKF is also used for parameter estimation and non-Gaussian data assimilation. The most common EnKF method for parameter estimation is state augmentation (Anderson, , 2001), which refers to simply including the parameters in the state vector. While EnKF with state augmentation often works well for certain parameters, such as autoregressive coefficients, the resulting linear update can be highly problematic for several types of parameters, most notably, variance and covariance parameters (Stroud and Bengtsson, , 2007; DelSole and Yang, , 2010). The linear shifting update in the EnKF can also be highly suboptimal for non-Gaussian observations. Alternative non-Gaussian EnKF approaches, including quantile-based filters (Anderson, , 2010) and moment-matching filters (Lei and Bickel, , 2011), still perform poorly for certain classes of measurement distributions. Thus, new EnKF methods are needed for more general state-space models with unknown parameters and non-Gaussian measurements.
1.2 Hierarchical state-space models
Here, we develop new methods to analyze a general class of high-dimensional hierarchical state-space models (HSSMs; e.g., Cressie and Wikle, , 2011). These models consist of four levels, adding two levels to a standard Gaussian SSM: a transformation level, which allows for non-Gaussian observations, and a parameter level, which models the unknown parameters, assumed to be random. This class of HSSMs is very general and can be used to describe most dynamic spatio-temporal systems observed sequentially over time. It allows for unknown time-varying parameters in any part of the model, non-Gaussian and nonlinear observations, mixture models, and even higher-order Markov models by extending the state-space.
More specifically, for discrete time points , our HSSM is given by:
[TABLE]
where is the -dimensional measurement vector at time , is a corresponding -dimensional “latent observation” vector, is the -dimensional unobserved state vector, is the parameter vector, is a known measurement or transformation distribution, is the state evolution operator, and denotes the -variate normal distribution. We assume that and are serially independent given the states and parameters. The model is completed with initial priors and \bm{\theta}_{0}\sim{\color[rgb]{0,0,0}p}(\bm{\theta}_{0}).
The state vector consists of one or more discretized spatial fields, along with any parameters that can be estimated using state augmentation. Conversely, the parameter vector contains any unknown quantities that cannot be handled well by state augmentation. The quantities , , , , and may all depend on the parameter vector . Sometimes we make the dependence explicit by writing and so forth. Further, we view the evolution operator as a “black box” that might be nonlinear and expensive to evaluate or unavailable in closed form. We focus on models where or are very large, usually . Sometimes, we observe so that where is the Dirac delta function.
The HSSM (1)–(4) is very flexible and includes many important SSMs as special cases. If and are known, the model reduces to a standard Gaussian SSM given by (2)–(3) (e.g., Shumway and Stoffer, , 2006). If the evolution is linear with , this becomes a conditionally Gaussian SSM (e.g., Carter and Kohn, , 1994; Chen and Liu, , 2000), where the full conditional distributions for the states are Gaussian. The HSSM also allows for non-Gaussian and nonlinear observations through in (1). This could include exponential family models (e.g., Shephard and Pitt, , 1997; Gamerman, , 1998), nonlinear mean functions (e.g. Carlin et al., , 1992), or threshold models for discrete data or discrete-continuous mixtures (e.g., Sansó and Guenni, , 1999). Finally, the HSSM allows for unknown static or time-varying parameters at any level of the model.
Note that it may be possible to reduce our HSSM to a standard nonlinear, non-Gaussian SSM by integrating out and defining as the state vector. However, we prefer the hierarchical formulation in (1)–(4), because it facilitates model building and computation. Specifically, by conditioning on and , we will use the EnKF for inference on the high-dimensional states.
1.2.1 Examples
We now present some examples of the HSSM in (1)–(4) that will be used for numerical comparisons in Section 5.
Example 1 (Heavy-tailed-data model).
*The distribution of belongs to the class of normal-scale mixtures (Andrews and Mallows, , 1974; West, , 1987) if we assume that is diagonal with independent random parameters on the diagonal. An important member of this class is the -distribution, which is similar to a normal distribution but has heavier tails and is thus more robust to outliers. The following HSSM assumes that the observation noise is independently -distributed with degrees of freedom: *
[TABLE]
where , is known, and follow independent inverse-gamma distributions. The transformation level is a point mass such that .
Next, we consider non-continuous observations, such as binary data and right-skewed data with a point mass at zero (e.g., hourly rainfall amounts). The EnKF and its transformation- or quantile-based extensions (e.g., Anderson, , 2009; Amezcua and Van Leeuwen, , 2014) are particularly poorly suited for assimilation of such observations. In our HSSM framework, this setting is handled using threshold observation models, in which is determined by transformations of the form , , where is a deterministic function containing an indicator function:
Example 2 (Threshold observation model for rainfall).
The observed data are modeled as where is the indicator function. Then follows a rainfall-type distribution, i.e., a mixture of a positive, right-skewed distribution and a point mass at zero, for some that determines the skewness (e.g., Sansó and Guenni, , 1999; Sigrist et al., , 2012). The HSSM is given as:
[TABLE]
where the parameter is either known or assigned a prior distribution.
The threshold framework can also be used for binary data following a Bernoulli distribution with , where is the standard normal CDF. This is called a probit model for . The transformation function is , which corresponds to the rainfall model with parameter .
Example 3 (Dynamic Poisson model).
For spatio-temporal count data, assume and , where accounts for overdispersion (Wikle, , 2002). The state follows a linear evolution, and the model includes six unknown time-varying parameters:
[TABLE]
where the transition matrix is tridiagonal with parameters and is a spatial covariance matrix with variance and range . The unknown time-varying parameters follow a random-walk evolution.
Example 4 (Nonlinear evolution model with unknown static parameter).
The Lorenz-96 model is widely used in atmospheric science as an example of a nonlinear evolution operator with chaotic behavior (Lorenz, , 1996). Lorenz-96 mimics advection at equally-spaced locations along a latitude circle on the globe. The time evolution is determined by
[TABLE]
with cyclic boundary conditions (i.e., , , ). Let denote the evolution operator implied by (5) with forcing parameter and time step , which leads to strong nonlinearities and non-Gaussian forecast distributions. We observe data at each location so that . The HSSM model is:
[TABLE]
where is known, is the nonlinear evolution operator, which is unavailable in closed form, is an unknown parameter, and is a known covariance matrix.
1.3 Extended ensemble Kalman techniques
The goal in this paper is to make Bayesian inference on the states and parameters for HSSMs (1)–(4). Let denote the data up to time , and use similar notation for and . Let denote the conditional distribution of some random variable given . We consider two types of inference. Filtering involves computing the filtering distribution sequentially for each time . Smoothing requires calculating the joint posterior (i.e., smoothing) distribution for fixed time .
We propose a simple general idea for such inference in high-dimensional HSSMs: Take an existing Bayesian inference technique (e.g., a particle filter or a Gibbs sampler), but approximate the part of the algorithm that requires sampling from or integrating out the high-dimensional state vector with an EnKF or with the related ensemble Kalman smoother (EnKS). This leads to a general class of extended ensemble Kalman filters and smoothers, which can be applied to high-dimensional HSSMs. The key is that when conditioning on and , our HSSM reduces to a standard Gaussian SSM with potentially nonlinear evolution, for which the EnKF is often well suited. Table 1 summarizes these algorithms and puts them in the context of existing particle approaches or approaches for linear Gaussian SSMs.
Our proposed filtering and smoothing algorithms can be classified into Ensemble Gibbs and Ensemble Marginalization methods. For Ensemble Gibbs schemes, the EnKF update is used within a Gibbs sampler to draw from the full conditional distribution of given and . For Ensemble Marginalization schemes, the EnKF is used to integrate out and sample from the posterior distribution of the states within a particle filter for the parameters , similar to a mixture KF (Chen and Liu, , 2000). Further connections to existing methods will be provided with each of our proposed algorithms.
Our new extended EnK techniques can be viewed from two perspectives. First, these methods greatly expand the class of problems to which the EnKF and EnKS can be applied. Second, our techniques allow for accurate, approximate inference in higher-dimensional models than standard statistical methods such as particle filters and Gibbs samplers.
In the EnKF literature, several approaches exist for estimating specific tuning parameters, such as inflation and localization parameters (e.g., Wang and Bishop, , 2003; Anderson, 2007a, ). Throughout, we assume that these tuning parameters are known, and we focus instead on inference of model parameters that explicitly appear in the model (1)–(4). Some approaches for state and parameter estimation in highly nonlinear state-space models (Gu and Oliver, , 2007; Chen and Oliver, , 2012; Bocquet and Sakov, , 2013) generate approximate posterior samples from the filtering or smoothing distribution using iterative optimization approaches. Specifically, the ensemble randomized maximum likelihood method (Chen and Oliver, , 2012) recasts posterior simulation as an optimization problem, and uses optimization methods with successive linearization of the model to find the posterior mode. These methods are quite general and can be used for parameter estimation and smoothing inference via state augmentation.
As the EnKF is an approximate technique for inference on the state vector, our extended EnK techniques also provide approximate inference. However, in many complex problems, exact inference is impossible, and approximate techniques such as approximate Bayesian computation (ABC; Sisson et al., , 2007; Marin et al., , 2012), synthetic likelihood (e.g., Wood, , 2010) and integrated nested Laplace approximations (INLA; Rue et al., , 2009; Lindgren et al., , 2011) have generally gained in popularity. While so-called exact approximate approaches (e.g., Andrieu et al., , 2010; Johansen et al., , 2012) result in exact inference given infinite computational resources, we will show here numerically that, in finite computation time, our extended EnK techniques can outperform exact approximate methods even for moderately high dimensions.
1.4 Overview and contributions
This article is organized as follows. In Section 2, we review some existing work on the EnKF and related issues from a statistical perspective. In Section 3, we provide preliminaries to lay the foundation for the development of our extended EnK techniques, including new insights on likelihood approximations in high-dimensional SSMs. In Section 4, we derive and present our extended ensemble Kalman filters and smoothers. In Section 5, we compare our methods to existing approaches using a number of simulated and real data examples. We find that our methods strongly outperform these existing approaches. Conclusions are given in Section 6.
2 Review of ensemble Kalman filter and smoother
In this section, we briefly review the ensemble Kalman filter (EnKF) and its extensions. The EnKF is a sequential Monte Carlo algorithm that generates samples from the state filtering distribution for the SSM (2)-(3), where and are known. Thus, throughout this section, we assume and are known and we often suppress from the conditioning set. Further, let denote an -variate normal density with mean and covariance matrix , evaluated at .
As mentioned in Section 1.2, we assume here that the evolution operator is intractable or computationally expensive, and so all we are able to do with is to “apply” it to a finite and typically rather small number of state vectors. Thus, sequential inference over time must be based on an ensemble of state vectors that represents (i.e., approximates) the relevant distributions at each time point.
More specifically, assume that the filtering ensemble at time is given by , meaning that the filtering distribution at time is approximated as
[TABLE]
where is the Dirac delta (point mass) function. The forecast ensemble is obtained by applying to each member of the filtering ensemble:
[TABLE]
2.1 Updating the state
Based on new data observed at time , we would like to update the forecast ensemble in (7) to obtain a filtering ensemble from the filtering distribution of the state vector, .
2.1.1 Update with a mixture prior
Based on (6), the prior distribution at time would in principle be given by
[TABLE]
where the integral is with respect to the measure induced by the distribution of . However, based on the mixture prior distribution in (8), the filtering distribution updated using the new observation would have the form
[TABLE]
with the filtering or posterior weights
[TABLE]
This reweighting of the th forecast ensemble member from prior weight to posterior weight is essentially what happens in an importance sampler or in the update step of a particle filter that integrates out the innovation error with covariance . This reweighting update is illustrated in a toy example in Figure 1 for the simple case , , and . However, as shown in Snyder et al., (2008) by setting their observation error covariance matrix to be in our notation, the particle weights degenerate, meaning that one mixture component would get all or nearly all of the weight, unless the ensemble size increases exponentially with the effective dimension given by the variance of the particle loglikelihood (see Section 3.2.1 below). Thus, in our high-dimensional setting where only a small ensemble size is computationally affordable, this update based on reweighting is problematic.
2.1.2 Update with a Gaussian prior
In contrast to the particle filter, the EnKF update does not reweight but rather shifts the forecast ensemble. The form of the shift is motivated by the fact that if the evolution operator is linear at each time , then the true prior distribution will be Gaussian. In that case, it is natural to approximate the mixture of Gaussians in the prior distribution (8) by a single Gaussian distribution:
[TABLE]
where and are estimates of the prior mean and covariance matrix (see Section 2.2 below) based on the forecast ensemble obtained in (7). Given this Gaussian prior, the posterior would also be Gaussian:
[TABLE]
where , , and
[TABLE]
is an estimate of the Kalman gain, .
2.1.3 The ensemble Kalman filter update
In the stochastic EnKF update, instead of calculating the posterior distribution in (12) explicitly, a sample from this distribution is obtained via shifting the forecast ensemble members by essentially carrying out conditional simulation (e.g., Katzfuss et al., , 2016). This leads to the following algorithm:
**Algorithm 1: ** **Ensemble Kalman filter (EnKF)
**Start with an initial ensemble . Then, at each time , do the following for each :
1. Forecast Step:
Compute forecast and prior ensemble members, and , respectively, where .
2. Update Step:
Shift each ensemble member as
[TABLE]
where is defined in (13), and is a pseudo-observation with .
The EnKF shifting update in (14) is derived based on assuming is linear and hence that the prior distribution is multivariate Gaussian. The forecast ensemble enters through the estimated Kalman gain , and , , and are evaluated at the parameter value , which is currently still assumed to be known.
The resulting filtering distribution converges to the true filtering distribution as , assuming the estimator converges to the true forecast covariance (see Section 2.2 below). For finite , in contrast to many other SMC methods, the EnKF update does not deteriorate for large , as long as a reasonable estimate of or can be found. It has been shown repeatedly (e.g., Lei et al., , 2010) that this shifting update also works surprisingly well for nonlinear evolution or non-Gaussian priors (see Figure 1 for an illustration), including operational, real-world applications (e.g., Houtekamer and Mitchell, , 2005; Bonavita et al., , 2010; Houtekamer et al., , 2014).
2.2 Estimation of matrices
The EnKF requires estimates of the mean vector and covariance matrix of the prior distribution in (11). The prior covariance matrix is required for the EnKF update through the Kalman gain (13), and both the mean and covariance matrix are needed in the EnKF likelihood in Section 3.2.2 below. The prior mean can simply be estimated by the sample mean of the forecast ensemble; that is, .
Estimation of the prior covariance matrix based on small ensembles is often challenging in applications with large state dimension ,
which can lead to ensemble collapse (e.g., Sætrom and Omre, , 2013) similar to the degeneration of the particle filter. In the EnKF, this collapse is counteracted by variance inflation (e.g., Anderson, 2007a, ), data-dimension-reduction via shrinkage regression (Sætrom and Omre, , 2011), or regularization via tapering. We focus here on tapering (e.g., Furrer et al., , 2006; Anderson, 2007b, ; Furrer and Bengtsson, , 2007), which is a very powerful method to avoid rank deficiency and spurious correlations by assuming that dependence is expected to decrease with increasing spatial distance. More precisely, tapering refers to entrywise multiplication with a sparse positive definite correlation matrix defined based on the distances between the elements. This leads to the expression
with , where
is the sample covariance matrix of , and is the sparse correlation matrix. Other regularization methods for states that are spatially referenced (e.g., Ott et al., , 2004; Hunt et al., , 2007) or not (e.g., Pourahmadi, , 2013) will not be discussed here, but could also be applied in some of our methods described later.
Also note that enters the Kalman gain in (13) through the matrix and the matrix , and so the EnKF can be stable even for massive , as long as is not too large.
2.3 Ensemble Kalman smoother
To approximate the smoothing distribution for some , various extensions of the EnKF to an ensemble Kalman smoother have been proposed (van Leeuwen and Evensen, , 1996; Evensen and van Leeuwen, , 2000; Khare et al., , 2008; Stroud et al., , 2010; Bocquet and Sakov, , 2014). The most widely used smoother is the EnKS algorithm of Evensen and van Leeuwen, (2000). Like the EnKF, the EnKS is a forward-only algorithm that is applicable even for intractable evolution operators . At each time , it essentially carries out an EnKF update on the augmented state that includes the entire history, :
**Algorithm 2: ** **Ensemble Kalman smoother (EnKS)
**Start with an initial ensemble . Then, at each time , do the following for each :
1. Forecast Step:
Compute , where .
2. Update Step:
For compute where is the pseudo-observation with , , and is a regularized version of the sample cross-covariance matrix of and .
The resulting ensemble is then an approximate sample from the smoothing distribution , and more generally, is an approximate sample from .
Note that regularization often results in for for some lag . The computations in Step 2 of Algorithm 2.3 then only have to be carried out for lags .
3 Preliminaries for extended ensemble Kalman techniques
In the previous section, we assumed that the latent and the parameters are all known. To obtain the necessary background and supporting results for the extended EnK techniques to be introduced in Section 4, we now extend the approaches reviewed in Section 2 to the more general HSSMs of the form (1)–(4) with non-Gaussian observations or unknown parameters.
Assume that we have a sample {\color[rgb]{0,0,0}\bm{\theta}_{1:t-1}^{(1)},\ldots,\bm{\theta}_{1:t-1}^{(M)}} from , the filtering distribution of the parameters at time . Except for special cases (see Section 3.1 below), we now need an ensemble of states, corresponding to each parameter sample {\color[rgb]{0,0,0}\bm{\theta}_{1:t-1}^{(i)}}; that is, we need an ensemble of ensembles. Specifically, assume that the joint filtering distribution of parameters and states at time is given by
[TABLE]
the marginal filtering distribution of is \textstyle{\color[rgb]{0,0,0}p(\bm{\theta}_{1:t-1}|\mathbf{z}_{1:t-1})=\textstyle\frac{1}{M}\sum_{i=1}^{M}\delta_{(\bm{\theta}_{1:t-1}^{(i)})}(\bm{\theta}_{1:t-1})}, and the conditional filtering distribution of given is
[TABLE]
Conditional on , the joint prior distribution at time is then given by
[TABLE]
Further, we have
[TABLE]
where we have used the same EnKF approximation as in (11). Note that and are obtained as described in Section 2.2 above using , but only based on the th forecast ensemble , where , .
Based on the approximation in (15), the joint distribution of all variables at time is
[TABLE]
where and may also depend on the parameters .
3.1 Forecast independence
An important special case is given by forecast independence of states and parameters (Frei and Künsch, , 2012), where . Forecast independence leads to enormous simplifications in the algorithms described later, because the prior distribution of the states can be calculated based on the entire ensemble of ensembles from the previous filtering distributions, instead of using only the corresponding ensemble as in (15); that is,
[TABLE]
where and are now based on all ensemble members . From this forecast distribution, we can then obtain the prior distribution by adding to the covariance matrix: .
Because this means that we do not need a separate ensemble for each previous parameter value, we can set and write the forecast ensemble as . Also, because the forecast distribution does not depend on , we only have to propagate the ensemble once from time to in the algorithms described later, which is especially important if is computationally expensive.
Forecast independence can be proved asymptotically if the parameters only appear in (see Frei and Künsch, , 2012). It also arises exactly if the parameters do not appear in and are temporally independent, , because then
[TABLE]
as in (17). Instead of the joint distribution in (16), it can then be shown that
[TABLE]
3.2 Approximation of the likelihood for known
It is also possible to obtain an approximation of the likelihood based on the EnKF. We first study likelihood approximations for the special case of known and observed (i.e., the hierarchical SSM given by (2)–(4)).
The distribution of all data up to time is then given by,
[TABLE]
where we have now made the conditioning on the parameters explicit. However, in the high-dimensional setting with potentially nonlinear evolution considered here, the filtered likelihood at time , \mathcal{L}_{t}(\bm{\theta}_{1:t})\colonequals{\color[rgb]{0,0,0}p(\mathbf{y}_{t}|\bm{\theta}_{1:t},\mathbf{y}_{1:t-1})}, is not available in closed form and cannot be exactly evaluated. Instead, in any parameter-inference procedure, we need to estimate this likelihood for given fixed parameter values based on the forecast ensemble . For the remainder of Section 3.2, we focus on a single time point and assume that is a sample from the true forecast distribution, .
For Bayesian sampling-based inference, the closer the likelihood estimator is to the true likelihood, the better the performance of the parameter-inference procedure. Specifically, unbiased estimators of the true likelihood lead to so-called exact approximate procedures that can be shown to sample from the correct target distribution; but the variability of the estimator is also crucial, in that the variance of the logarithm of the estimator of the likelihood determines the performance of the parameter inference for finite computation time (Andrieu and Roberts, , 2009). Variances smaller than 2 are often recommended for satisfactory performance (e.g., Doucet et al., , 2015).
3.2.1 Particle likelihood
If the mixture prior distribution in (8) were used, the filtered likelihood at time would also be a mixture of normals,
[TABLE]
Except for the fact that we are able to integrate out the innovation error with covariance here due to the linear Gaussian observations, is the likelihood approximation used in several inference procedures based on particle filters (e.g., Andrieu et al., , 2010; Johansen et al., , 2012) and in some approaches in the EnKF literature (e.g., Ueno and Nakamura, , 2014, 2016). Using the notation from (10), we can write this likelihood as . While the and thus also the particle likelihood are unbiased for the true likelihood (Del Moral, , 2004), leading to exact approximate algorithms, the variance of the particle likelihood can be very large in high dimensions, resulting in algorithms with potentially poor performance in finite computation time. We now illustrate this using a simple example.
Example 5.
For a single time point, consider the forecast distribution , with model error covariance , noise covariance , observation matrix , and thus ; that is, . The unknown parameter is .
In the setting of Example 5, Figure 2 shows the distribution of the particle likelihood , which becomes increasingly skewed for increasing . Further, it shows the trace plot for a “pseudo-marginal” Metropolis-Hastings (MH) sampler targeting for using , whose large variances cause the Markov chain to get “stuck” after a particularly large likelihood is obtained. Finally, we explore the ensemble size necessary to keep the variance of the particle loglikelihood, which was averaged over 100 realizations of , to below 2, and this number clearly scales exponentially with the dimension . This result can also be shown analytically:
Proposition 1**.**
In the setting of Example 5, we have as , where the variance is evaluated with respect to , and the are iid from the true forecast distribution.
All proofs are provided in Appendix A.
The results obtained here for a particle or importance-sampling likelihood at a single time point extend directly to a particle filter in a state-space setting. While a particle filter used to approximate a likelihood in a sequential setting will typically include a resampling step at every time point (e.g., Andrieu et al., , 2010), the resulting forecast or prior distribution at the next time point can be at most as good as the iid sampling from the true prior considered here, assuming that is intractable, and so the particle-filter approximation to the likelihood will also degenerate in high dimensions.
3.2.2 EnKF likelihood
Due to the problems of the PF likelihood in high dimensions, we use the EnKF approximation of the prior distribution in (11) to obtain an EnKF likelihood:
[TABLE]
where , , and and may depend on through and .
The EnKF likelihood has been successfully used in several examples and applications (e.g., Mitchell and Houtekamer, , 2000; Stroud et al., , 2010; Frei and Künsch, , 2012; Stroud et al., , 2018). For linear SSMs, the true forecast distribution and likelihood are both exactly Gaussian, and so the EnKF likelihood converges to the true likelihood as , assuming no localization for . For finite , in contrast to the PF likelihood, the EnKF likelihood is a biased estimator of the true likelihood, and thus does not lead to exact approximate inference. Special care should be taken regarding parameters influenced by tapering, as, for example, a spatial-scale parameter in the innovation covariance might be overestimated to compensate for an overly restrictive taper in the forecast covariance (see Section 2.2). The estimated parameter is then useful for running an EnKF with the same taper settings, but it might not be the “true” spatial scale of the innovation covariance in the model (1)–(4).
However, more importantly, the variance of the EnKF (log)likelihood can be relatively small even in high dimensions, as long as a good regularized estimator can be found. This can be viewed as a manifestation of the well-known trade-off between bias and variance that is central to statistical learning with “big” data (e.g., Hastie et al., , 2009). To illustrate this point, we consider the EnKF likelihood in the setting of Example 5 with a diagonal . Figure 2 shows that the EnKF likelihood is less skewed, less variable, and thus closer to the true likelihood than the particle likelihood. This also leads to much better mixing for the MH sampler, and the resulting estimated posterior for (not shown) is nearly indistinguishable from the estimated posterior using the true likelihood. Figure 2(2(d)) indicates that for the EnKF likelihood, only has to increase linearly with the dimension to keep the variance of the loglikelihood bounded. This can also be shown analytically:
Proposition 2**.**
In the setting of Example 5, we have as , where the variance is evaluated with respect to , and the are iid from the true forecast distribution.
To explore whether this result holds more broadly than just in the independent case, we repeated the simulations in Figure 2(2(d)) for a slightly modified version of Example 5 with a tridiagonal forecast covariance matrix . This produced a plot very similar to Figure 2(2(d)).
If the true forecast distribution, and thus also the true likelihood, is non-Gaussian, the EnKF likelihood converges to the Dawid–Sebastiani scoring rule (Dawid and Sebastiani, , 1999) for , and so, from a frequentist perspective, parameter inference based on the EnKF likelihood at a single time point can be interpreted as approximate minimum contrast estimation (Birgé and Massart, , 1993; Pfanzagl, , 1969) under the Dawid–Sebastiani score. Due to the the propriety of this scoring rule, the minimum-contrast estimator can be shown, using the proof of Wald, (1949), to be consistent under regularity conditions (cf. Gneiting and Katzfuss, , 2014, Sect. 3.4.2). From a Bayesian perspective, inference based on the EnKF likelihood essentially amounts to computing the approximate “Gibbs posterior” using the Dawid–Sebastiani score as the empirical risk function (Jiang and Tanner, , 2008; Li et al., , 2013; Katzfuss and Bhattacharya, , 2015). Under certain conditions, most importantly that is consistent and asymptotically normal as , it might be possible to scale and rotate posterior samples of to reflect the true asymptotic uncertainty in the frequentist sense (Shaby, , 2014), but this will be pursued in future work.
3.3 Integrated likelihood
We now extend the EnKF likelihood for known from Section 3.2.2 to the full HSSM with the general transformation layer (1): . This integrated EnKF likelihood of the data at time , with and integrated out, is an important component of some of the algorithms described later. Using (16), we have
[TABLE]
where is defined in (20), and , , and and may depend on .
Evaluating this likelihood requires carrying out an integral over the -dimensional vector . Strategies for doing so depend on the choice of transformation distribution . If the data are Gaussian and the transformation is simply the identity, the integral degenerates and analytical integration is trivial. If it is possible to sample from , the integral can be approximated using these posterior samples (e.g., Chib and Jeliazkov, , 2001; Raftery et al., , 2007). Finally, if is a continuous function in , a Laplace approximation (Tierney and Kadane, , 1986) can be used to approximate the integral. As we are already approximating as a latent Gaussian field in , the Laplace approximation of the integral should often result in negligible additional error, especially if the observations in are independent given and (e.g., Rue et al., , 2009); see Section 5.2 for an example.
4 Extended ensemble Kalman techniques for hierarchical state-space models
For HSSMs of the form (1)–(4) with non-Gaussian observations or unknown parameters, we now propose extended EnK techniques. The general idea is to approximate the inference on the state vector in suitable existing Bayesian inference techniques using the EnKF or the EnKS.
4.1 Extended ensemble Kalman filters
We begin by discussing extended EnKFs that are suitable for filtering inference, which at each time aims to find the filtering distribution given the data collected up to time . Smoothing inference is discussed in Section 4.2 below.
4.1.1 Gibbs ensemble Kalman filter (GEnKF)
First, we consider a combination of an EnKF and a Gibbs sampler (Geman and Geman, , 1984; Gelfand and Smith, , 1990) for filtering inference in HSSMs. For this algorithm we assume
that the parameters do not appear in and are temporally independent (i.e., ). As noted in Section 3.1, this implies forecast independence, which enables inference based on a single ensemble (i.e., ). At each time , our proposed Gibbs ensemble Kalman filter (GEnKF) iteratively samples , , and from their respective full conditional distributions (FCDs), which are proportional to the joint distribution in (3.1). The state is sampled by drawing from the forecast ensemble and using the EnKF update, which allows application of the GEnKF to high-dimensional and nonlinear SSMs.
**Algorithm 3: ** **Gibbs ensemble Kalman filter (GEnKF)
**Starting with an initial ensemble , where (\bm{\theta}_{0}^{(i)},\mathbf{x}_{0|0}^{(i)})\sim{\color[rgb]{0,0,0}p(\bm{\theta}_{0},\mathbf{x}_{0})}=p_{0}(\bm{\theta}_{0})\mathcal{N}_{n}(\mathbf{x}_{0}|\bm{\mu}_{0|0},\bm{\Sigma}_{0|0}), do the following for each :
Obtain the forecast ensemble as , for all . 2. 2.
Find starting values for and for all . 3. 3.
For , iterate between the following steps until convergence:
- (a)
Sample from by sampling uniformly from and computing the EnKF update \mathbf{x}_{t|t}^{(i)}=\mathcal{U}_{t}(\mathbf{\widetilde{x}}_{t|t-1}^{({\color[rgb]{0,0,0}j})}|\mathbf{y}_{t}^{(i)},\mathbf{\widetilde{x}}_{t|t-1}^{(1:M)},\bm{\theta}_{t}^{(i)}) (see (14)). 2. (b)
Sample from {\color[rgb]{0,0,0}p(\mathbf{y}_{t}|\mathbf{x}_{t|t}^{(i)},\bm{\theta}_{t}^{(i)},\mathbf{z}_{1:t})}\propto f_{t}(\mathbf{z}_{t}|\mathbf{y}_{t},\bm{\theta}_{t}^{(i)})\mathcal{N}_{m_{t}}(\mathbf{y}_{t}|\mathbf{H}_{t}(\bm{\theta}_{t}^{(i)})\mathbf{x}_{t|t}^{(i)},\mathbf{R}_{t}(\bm{\theta}_{t}^{(i)})). 3. (c)
Sample from {\color[rgb]{0,0,0}p(\bm{\theta}_{t}|\mathbf{y}_{t}^{(i)},\mathbf{x}_{t|t}^{(i)},\mathbf{z}_{1:t})}\propto f_{t}(\mathbf{z}_{t}|\mathbf{y}_{t}^{(i)},\bm{\theta}_{t})\allowbreak\mathcal{N}_{m_{t}}(\mathbf{y}_{t}^{(i)}|\mathbf{H}_{t}(\bm{\theta}_{t})\mathbf{x}_{t|t}^{(i)},\mathbf{R}_{t}(\bm{\theta}_{t}))\allowbreak{\color[rgb]{0,0,0}\mathcal{N}_{n}(\mathbf{x}_{t}|\widehat{\bm{\mu}}_{t|t-1},\widetilde{\bm{\Sigma}}_{t|t-1}+\mathbf{Q}_{t}(\bm{\theta}_{t}))}\allowbreak{\color[rgb]{0,0,0}p}(\bm{\theta}_{t}).
After convergence, is a joint sample from , and is the filtering ensemble at time .
Note that iterating through Step 3(a)–(c) of Algorithm 4.1.1 can be done completely independently (i.e., in parallel) for . To sample in Step 3(a), we use a systematic sample across , so that each forecast ensemble member is chosen exactly once for each iteration of the Gibbs sampler. It can be shown under mild regularity conditions, that the closer the forecast ensemble is to multivariate normal, the closer the EnKF in Step 3(a) gets to sampling from the true FCD of as increases, and hence the closer the resulting Markov chain will be to the Markov chain that would be obtained using an exact Gibbs sampler, in terms of total variation distance (Alquier et al., , 2016). As the GEnKF requires sampling from the FCDs of and in (b) and (c), respectively, the algorithm is most efficient if the distributions and are such that these FCDs are available in closed form. Otherwise, Metropolis updates are necessary for these variables, which is generally only feasible if the number of parameters in is not too large, and if the data in are independent given .
For computational feasibility of Algorithm 4.1.1, it is crucial to find good starting values for and , in order to minimize the required number of Gibbs iterations. In some applications, it is natural to use the values from the previous time point (i.e., set the initial values as and ). Another possibility is to choose the value of that maximizes the integrated likelihood in (21).
As an example, for the model with heavy-tailed data in Example 1, the full-conditional distribution of is a closed-form inverse-Gamma distribution, and so we can apply the following special case of the GEnKF in Algorithm 4.1.1:
Example 6 (Robust ensemble Kalman filter (REnKF)).
Starting with an initial ensemble , where , do the following for each :
Obtain the forecast ensemble as \mathbf{\widetilde{x}}_{t|t-1}^{(i)}={\color[rgb]{0,0,0}\mathcal{M}_{t}(\mathbf{x}_{t-1|t-1}^{(i)})}, for all . 2. 2.
Set , and for all and . 3. 3.
For , iterate between the following steps until convergence:
- (a)
Sample from by sampling uniformly from and computing the EnKF update (see (14)). 2. (b)
Sample \textstyle\theta_{t,l}^{(i)}|\mathbf{y}_{1:t},\mathbf{x}_{t|t}^{(i)}\stackrel{{\scriptstyle ind}}{{\sim}}\mathcal{IG}\Big{(}\frac{\kappa+1}{2},\frac{\kappa}{2}+\Big{(}\frac{y_{t,l}-(\mathbf{H}_{t}\mathbf{x}_{t|t}^{(i)})_{l}}{{\color[rgb]{0,0,0}\sigma_{t}}}\Big{)}^{2}/2\Big{)}, .
If is unknown with prior , it could be sampled in a collapsed Gibbs sampling scheme prior to updating by sampling from .
Note that the iterative Huber-EnKF (Rao et al., , 2017) seems to be similar in spirit to our robust EnKF, but it uses the Huber norm in place of the loss function implied by a -distribution.
The GEnKF can also be applied to threshold models as in Example 2:
Example 7 (GEnKF for rainfall threshold model).
Starting with an initial ensemble , where , do the following for each :
Obtain the forecast ensemble as , for all . 2. 2.
Find starting values for for all . 3. 3.
For , iterate between the following steps until convergence:
- (a)
Sample from by sampling uniformly from and computing the EnKF update (see (14)). 2. (b)
Sample from {\color[rgb]{0,0,0}p(\theta_{t}|\mathbf{z}_{t},\mathbf{x}_{t})}\propto{\color[rgb]{0,0,0}p(\theta_{t})}\prod_{l:z_{t,l}>0}\mathcal{N}\big{(}z_{t,l}^{1/\theta_{t}}|(\mathbf{H}_{t}\mathbf{x}_{t})_{l},{\color[rgb]{0,0,0}\sigma_{t}^{2}}\big{)}\,z_{t,l}^{(1/\theta_{t})-1}/\theta_{t}. 3. (c)
*For , sample from {\color[rgb]{0,0,0}p(y_{t,l}|z_{t,l},\mathbf{x}_{t},\theta_{t})}=\begin{cases}\delta_{z_{t,l}^{1/\theta_{t}}}(y_{t,l}),&z_{t,l}>0,\\ \mathcal{N}^{-}\big{(}(\mathbf{H}_{t}\mathbf{x}_{t}^{(i)})_{l},\sigma_{t}^{2}\big{)},&z_{t,l}=0,\end{cases} *
where denotes a truncated normal distribution on the negative real line.
A Gibbs ensemble Kalman filter could in principle also be applied if forecast independence does not hold. This is difficult, however, because {\color[rgb]{0,0,0}p(\bm{\theta}_{t},\mathbf{x}_{t}|\mathbf{z}_{1:t})=\int p(\bm{\theta}_{t},\mathbf{x}_{t}|\bm{\theta}_{t-1},\mathbf{z}_{1:t})dP(\bm{\theta}_{t-1}|\mathbf{z}_{1:t}),} where it is straightforward to sample from using a Gibbs sampler similar to above, but it is highly challenging to compute the probability or weight , which will also degenerate in high dimensions. Thus we will not pursue this general case here.
Perhaps the approaches in the literature most closely related to our GEnKF are Myrseth and Omre, (2010) and Tsyrulnikov and Rakitko, (2017), who model the forecast covariance matrix as random with an inverse-Wishart prior. This matrix is not a model parameter, which makes Bayesian inference difficult. Myrseth and Omre, (2010) propose a Gibbs-like algorithm with a single iteration, but they do not take the data into account when calculating the posterior of the forecast covariance. Tsyrulnikov and Rakitko, (2017) seem to ignore the fact that forecast independence does not hold in their model.
4.1.2 Particle ensemble Kalman filter (PEnKF)
For high-dimensional SSMs, particle filters typically degenerate (e.g., Doucet et al., , 2001; Snyder et al., , 2008). While Rao-Blackwellized particle filters (RBPFs; Doucet et al., 2000a, ) can scale to higher dimensions by integrating out part of the state analytically, for example using a Kalman filter (Chen and Liu, , 2000), the integral is generally not available for HSSMs with nonlinear evolution or non-Gaussian data. An exact approximation of RBPFs proposed in Johansen et al., (2012) uses “local” particle filters to carry out the integration, thus employing the particle likelihood from Section 3.2.1, which was shown to scale poorly as well.
As an alternative, we propose here a particle ensemble Kalman filter (PEnKF), which uses the EnKF likelihood from Sections 3.2.2 and 3.3 to approximately integrate out the high-dimensional state. The PEnKF is suitable for HSSMs for which forecast independence does not hold, or for which efficient sampling from the FCDs in Step 3 of Algorithm 4.1.1 is not possible. As a particle filter is applied to , the PEnKF will work best if the number of unknown parameters is not too high. Although this algorithm does not require forecast independence for the parameters as the GEnKF above, forecast independence does result in considerable simplifications (see below).
**Algorithm 4: ** **Particle ensemble Kalman filter (PEnKF)
**Initialize the algorithm with an equally weighted ensemble (\bm{\theta}_{0}^{(i)},\mathbf{x}_{0|0}^{(i,j)})\sim{\color[rgb]{0,0,0}p(\bm{\theta}_{0},\mathbf{x}_{0})}={\color[rgb]{0,0,0}p}(\bm{\theta}_{0})\mathcal{N}_{n}(\mathbf{x}_{0}|\bm{\mu}_{0|0},\bm{\Sigma}_{0|0}) with , . Then, for :
For :
- (a)
Sample a particle from a proposal distribution . 2. (b)
Propagate the corresponding ensemble: , . 3. (c)
Calculate the particle weight:
. 4. (d)
Generate from for (see below). 2. 2.
The filtering distribution is given by
{\color[rgb]{0,0,0}p(\bm{\theta}_{t},\mathbf{x}_{t}|\mathbf{z}_{1:t})}\approx\sum_{i=1}^{M}w_{t}^{(i)}\textstyle\frac{1}{N}\sum_{j=1}^{N}\delta_{(\bm{\theta}_{t}^{(i)},\mathbf{x}_{t|t}^{(i,j)})}(\bm{\theta}_{t},\mathbf{x}_{t}). 3. 3.
If desired, resample to obtain an equally weighted sample (see, e.g., Douc et al., , 2005, for a comparison of resampling schemes).
By noting that , the expression for the weights in Step 1(c) is derived as follows:
[TABLE]
where {\color[rgb]{0,0,0}p(\bm{\theta}_{t},\mathbf{z}_{t}|\bm{\theta}_{1:t-1},\mathbf{z}_{1:t-1})=p(\bm{\theta}_{t}|\bm{\theta}_{1:t-1},\mathbf{z}_{1:t-1})p(\mathbf{z}_{t}|\bm{\theta}_{1:t},\mathbf{z}_{1:t-1})=p(\bm{\theta}_{t}|\bm{\theta}_{t-1})p(\mathbf{z}_{t}|\bm{\theta}_{1:t},\mathbf{z}_{1:t-1})}, and {\color[rgb]{0,0,0}p(\mathbf{z}_{t}|\bm{\theta}_{1:t}^{(i)},\mathbf{z}_{1:t-1})} is approximated by the integrated likelihood in (21).
Thus, the PEnKF is similar to the mixture KF (Chen and Liu, , 2000) and a Rao-Blackwellized particle filter, except that the high-dimensional state is integrated out and sampled using the EnKF approximation. This again allows application to high-dimensional and nonlinear SSMs.
In the numerical examples in Section 5, for simplicity we set the proposal distributions equal to the prior distributions implied by the model: . Of course, it might be possible to adapt existing, more sophisticated strategies for choosing efficient proposal distributions to our setting.
For Step 2 in Algorithm 4.1.2, a sample from needs to be computed for each particle with significantly nonzero weight . This sample can be obtained using a special case of the Gibbs sampler in Algorithm 4.1.1 with held fixed (i.e., skipping Step 3(c)). In settings for which the Laplace approximation works well (see Section 3.3), this Gibbs sampler could consist of only a single iteration, in which is first sampled from the Gaussian approximation to its FCD used in the Laplace approximation, and then is sampled as in Step 3(a) of Algorithm 4.1.1. Here is an example of a PEnKF using this Laplace approach for the Poisson model in Example 3:
Example 8 (PEnKF for dynamic Poisson model).
Initialize and , . Then, for , and each :
Sample from . 2. 2.
Propagate the corresponding ensemble: , . 3. 3.
Using Newton-Raphson, find the mode , and the negative inverse Hessian at the mode, of the function , where . 4. 4.
*Calculate the particle weight , where *
. 5. 5.
Ensemble update: , , where . 6. 6.
Resample with probability proportional to ; then set .
Algorithm 4.1.2 is, in principle, also applicable when is constant over time (i.e., the prior distribution is a point mass at the previous value). However, in this case the proposal distribution also has to be a point mass, which means that over time there will be fewer and fewer unique particles. To avoid this, the particles can be resampled based on a kernel density estimate of the filtering distribution (Liu and West, , 2001; Frei and Künsch, , 2012), but then the state ensemble has to be propagated and updated again based on the newly sampled parameter values.
As described in Section 3.1, in the case of forecast independence, the likelihood can be computed based on the entire filtering ensemble from , and so we can set and carry out Algorithm 4.1.2 with a single ensemble instead of an ensemble of ensembles. This important special case is described in Algorithm 1 of Frei and Künsch, (2012). If does not depend on but is computationally expensive, it can then be advantageous to “oversample” the parameters, so that there are more parameter particles than state ensemble members.
4.2 Smoothing using extended ensemble Kalman techniques
We now discuss extended EnKS methods for HSSMs when smoothing inference is desired; that is, when interest is in for some fixed time period .
4.2.1 Smoothing inference for high-dimensional parameter vectors
If the parameters and are high-dimensional but they can be efficiently sampled from their full-conditional distributions (FCDs), a Gibbs-EnKS can be applied:
**Algorithm 5: ** Gibbs ensemble Kalman smoother (GEnKS)
Initialize and . 2. 2.
Iterate between following steps for until convergence:
- (a)
Obtain samples {\color[rgb]{0,0,0}\mathbf{x}_{1:T|T}^{(i,1:N)}} from {\color[rgb]{0,0,0}p(\mathbf{x}_{1:T}|\bm{\theta}_{1:T}^{(i-1)},\mathbf{y}_{1:T}^{(i-1)})} using the EnKS (Algorithm 2.3), and set {\color[rgb]{0,0,0}\mathbf{x}_{1:T|T}^{(i)}=\mathbf{x}_{1:T|T}^{(i,j)}} for a sampled uniformly at random from . 2. (b)
Draw a sample from {\color[rgb]{0,0,0}p(\mathbf{y}_{1:T}|\mathbf{z}_{1:T},\mathbf{x}_{1:T|T}^{(i)},\bm{\theta}_{1:T}^{(i-1)})}. 3. (c)
Draw a sample from {\color[rgb]{0,0,0}p(\bm{\theta}_{1:T}|\mathbf{x}_{1:T|T}^{(i)},\mathbf{y}_{1:T}^{(i)},\mathbf{z}_{1:T})}.
The GEnKS can be interpreted as an alternative version of the Gibbs-KF of Carter and Kohn, (1994) or the particle-Gibbs algorithm of Andrieu et al., (2010), where the Kalman or particle filter, respectively, is replaced by the EnKF. This allows the GEnKS to be applied to high-dimensional, nonlinear SSMs. In Section 5.3, we conduct a numerical comparison using the Lorenz-96 model from Example 4 between particle Gibbs and the following GEnKS:
Example 9 (GEnKS for Lorenz-96 with static parameter).
Initialize , set , and then iterate between following steps for until convergence:
Obtain samples from using the EnKS (Algorithm 2.3), and set for a sampled uniformly from . 2. 2.
With , draw a sample from
[TABLE]
4.2.2 Smoothing inference for low-dimensional parameters
If cannot be sampled efficiently from its FCD in Step 1(c) of Algorithm 4.2.1 but is low-dimensional, the following smoothing algorithms could be applied.
Static parameters
First, consider the case of a static (i.e., temporally constant) low-dimensional parameter vector. In this case, it is possible to apply a Metropolis-Hastings-EnKS (MHEnKS). This algorithm consists of proposing a new parameter value, say, and then accepting this value with a probability depending on the likelihood
[TABLE]
for , where the terms in the product on the right-hand side are the integrated likelihoods at time in (21), which can be obtained sequentially by carrying out an EnKF and computing the integrated likelihood at each time . Then, we can obtain samples from for each thinned posterior sample using a special case of Algorithm 4.2.1, in which we hold fixed and skip Step 2(c). If the evolution is linear, it is also possible to reuse the filtering ensemble already obtained in the EnKF by applying the EnKS from Stroud et al., (2010).
The MHEnKS can be interpreted as an alternative version of the particle marginal Metropolis-Hastings algorithm of Andrieu et al., (2010), where the particle likelihood in (19) is replaced by the EnKF likelihood in (20). The MHEnKS should hence scale to higher state dimensions.
Time-varying parameters
If varies over time, inference becomes more challenging, in that even if each consists of a small number of elements, the combined vector of unknown parameters at all time points, , can still be rather high-dimensional if is large. Hence, for a low-dimensional time-varying parameter vector, we propose the following particle-EnKS:
**Algorithm 6: ** Particle ensemble Kalman smoother (PEnKS)
Obtain smoothing trajectories and corresponding weights by running the PEnKF in Algorithm 4.1.2 up to time , and saving and potentially resampling the trajectories at each time (Kitagawa, , 1996). 2. 2.
If desired, improve degeneracy by obtaining new state trajectories from for each particle trajectory with significantly nonzero weight , using a special case of Algorithm 4.2.1, in which we hold fixed and skip Step 2(c).
Note that this algorithm might perform poorly when is large, because of degeneracy problems at “early” time points (e.g., Doucet et al., 2000b, ). For Step 1, we can obtain the same simplifications in the case of forecast independence as for the “stand-alone” PEnKF described in Section 4.1.2. For Step 2, we can re-use state ensembles in the case of linear evolution using the EnKS procedure of Stroud et al., (2010).
4.2.3 A more general Markov chain Monte Carlo – ensemble Kalman smoother
An extension of the Gibbs-EnKS in Algorithm 4.2.1, which we will call an MCMC–EnKS, is possible as long as can be split into a set of parameters , where each is itself either low-dimensional or has a closed-form full-conditional distribution, and can be updated as in the Metropolis-Hastings-, Gibbs-, or particle-EnKS. This allows the fitting HSSMs with a variety of unknown parameters.
4.3 Properties of the extended ensemble Kalman algorithms
4.3.1 Convergence
Proposition 3**.**
If the evolution operators are linear and is a consistent estimator of the forecast covariance matrix, the empirical distributions for the states and parameters produced by the PEnKF, GEnKS, MHEnKS, and PEnKS all converge to the true filtering or smoothing distributions under mild regularity conditions, as , , and the number of MCMC iterations all tend to infinity.
Our informal proof is derivative of existing results and provided in Appendix A. Note that is consistent if the tapering range in Section 2.2 tends to infinity as . The convergence in Proposition 3 does not hold exactly in the case of the GEnKF, in that the marginal filtering distribution of integrated over and and the forecast distribution of are not exactly Gaussian anymore, but the GEnKF treats them as such.
For SSMs with nonlinear evolution, the EnKF essentially converges to the best filter among those which update the state linearly based on the data at each time point (e.g., Le Gland et al., , 2011; Law et al., , 2016). This is similar to the notions in so-called linear Bayes estimation that considers linear estimators given only the first and second moments of the prior distribution and likelihood (e.g, Hartigan, , 1969; Goldstein and Wooff, , 2007). In addition, there is a close connection between approximate Bayes computation (ABC) and linear Bayes estimation (e.g., Nott et al., , 2014) and the EnKF can be shown to be equivalent to special cases of ABC algorithms (Nott et al., , 2012).
As argued in Section 3.2.2, the EnKF likelihood converges to the Dawid-Sebastiani score, which is a proper scoring rule. Thus, for nonlinear evolution the extended EnKF and EnKS techniques will not produce exact inference, but the resulting samples of the state and parameters can still provide useful approximations to the desired distributions in many settings.
For the high-dimensional applications of interest here, only small ensembles are computationally feasible, and so the algorithms’ performance for small and is much more important than their asymptotic properties for and tending to infinity. Our algorithms will tend to perform well for HSSMs for which the embedded SSM (2)–(3) is well suited for inference using the EnKF and EnKS. This is the case for many real-world applications (e.g., Houtekamer and Mitchell, , 2005; Bonavita et al., , 2010; Houtekamer et al., , 2014). We have also seen in Section 3.2, that the EnKF approximation to the likelihood can be much less variable than the particle approximation, which should lead to better algorithms. Further numerical comparisons are provided in Section 5 below.
4.3.2 Computational cost
The EnKF requires storing, propagating, and operating on vectors (i.e., ensemble members) of length , and the computational complexity of each EnKF update is for most EnKF variants (e.g., Tippett et al., , 2003). The computations at each time step are highly parallelizable.
For our extended EnKF algorithms, the EnKF generally has to be carried out several times. However, often only a small number of iterations or particles is necessary to obtain a reasonable approximation to the desired distribution (see Section 5 below). Note that in machine learning, often only a few or even a single update (called contrastive divergence) is used to sample from restricted Boltzman machine Markov random fields for training deep neural networks (Hinton, , 2002). Although this adds bias, empirical and some theoretical support suggests that it is an effective and efficient approach in the deep learning context (e.g., Bengio and Delalleau, , 2009). In addition, in the case of forecast independence for the GEnKF and PEnKF, only the update step has to be performed several times, while the forecast step only has to be carried out for a single ensemble. Thus, considering that the most expensive part of the EnKF is often the forward propagation using an expensive evolution model , the increased computational cost of the GEnKF or PEnKF relative to the EnKF can be negligible.
Thus, the EnKF and our proposed extensions are scalable to very high state dimensions and large numbers of observations . Parameters that cannot be included in the state can also be high-dimensional if closed-form full-conditional distributions are available, so that our Gibbs ensemble Kalman techniques can be applied. However, operational, large-scale implementation for millions of dimensions or more is a significant undertaking. While efficient and parallel existing EnKF implementations such as the data assimilation research testbed, or DART (Anderson et al., , 2009), should be suitable for the proposed extensions, the application of EnKF in operational, large-scale environments typically requires serious high-performance-computing expertise and extensive tuning of inflation and localization parameters.
5 Examples and numerical comparisons
The extended EnK methods described in Section 4 are very general and can be applied in many situations. We now present numerical comparisons for the examples in Section 1.2.1.
5.1 Simulation study for non-Gaussian observations
First, we considered three non-Gaussian observation scenarios with temporally independent parameters with closed-form full conditionals, which can be handled using our GEnKF (Algorithm 4.1.1). In particular, we consider heavy-tailed observation noise from a distribution with degrees of freedom from Example 1 (see Figure 3(a)), and the rainfall threshold model from Example 2 with the true assumed known (see Figure 3(b)) or unknown.
We implemented the special cases of the GEnKF presented in Examples 6 and 7 for the heavy-tailed and rainfall data, respectively. In the case of known for the rainfall data, we simply skipped Step 3(b) in Example 7. The GEnKFs were compared to a standard sequential importance resampling (SIR) particle filter, and to the true posterior distribution obtained using an MCMC that was run for 1,000 iterations after 500 burn-in iterations. We also considered a commonly used extension of the EnKF update (Anderson, , 2003): For an observation equation of the form , the pseudo-observation in (14) is then replaced by and the Kalman gain is estimated as , where and are regularized versions of the empirical covariance matrix of and the cross-covariance matrix between and , respectively. For the rainfall data with unknown , we used state augmentation. For the heavy-tailed data, we also implemented the robust Huber-EnKF proposed in Roh et al., (2013). The Huber-EnKF requires a cutoff or clipping value for each observation location, which we selected by sampling 1000 times from the true forecast distribution based on an update consisting of a single observation at the respective location.
Because non-Gaussian observations are not directly related to the evolution model and the forecast step, we focused on the update step at a single time point . We simulated 100 true state vectors of size on a one-dimensional spatial domain with randomly chosen observation locations. The prior of the state was taken to be a normal distribution with known true mean 0.2 (constant over space) and known true covariance matrix based on a powered exponential covariance function with power 1.8 and scale parameter 10, same for all three simulation scenarios. Also, we assumed to be known for all examples.
For the standard EnKF, the Huber-EnKF, and the particle filter (PF), we used the same prior ensemble of size , simulated from the true prior distribution. For our GEnKF we only used a subset of size of this ensemble, and only one Gibbs iteration for rainfall with known , and three iterations for the other two scenarios. For the EnKF and GEnKF, we used the same Wendland taper function with tapering length 20.
The results are shown in Figure 3 and Table 2. We computed both the mean squared prediction error (MSPE), which evaluates the accuracy of the posterior mean as a point prediction, and the continous ranked probability score (CRPS), which evaluates the quality of the entire predictive distribution, including its shape and spread (see, e.g., Gneiting and Katzfuss, , 2014). For the heavy-tailed observations, we compared the state filtering distribution for the different methods to the true state, whereas for rainfall data we considered predictions of the true rainfall, x_{t,l}{\color[rgb]{0,0,0}\mathbbm{1}_{\{x_{t,l}>0\}}}.
The GEnKF outperformed all competitors in all three scenarios, and came relatively close to the true posterior distribution. This is especially remarkable considering the GEnKF consisted of at most three EnKF updates with ensemble size 30, while the competitors each used an ensemble of size 100. Hence, in our simulations, the GEnKF produced the most accurate results at the lowest computational expense.
In Figure 4, we further investigated the computational costs by comparing the computation time and MSPE for the update in the GEnKF and the particle filter for varying ensemble sizes between 5 and 100 and number of particles between 30 and 1 million. The results ignore the potentially significant cost of applying evolution operators, to the benefit of the particle filter. All results are averaged over 20 repetitions. While the GEnKF obtains close to optimal MSPE even for very low computation times, the particle filter produces highly non-optimal MSPEs even for computation times that are several orders of magnitude larger. The differences become even more pronounced for larger datasets.
5.2 Cloud data
We also conducted filtering inference for a cloud-motion dataset shown in Figure 5(a) from Wikle, (2002), representing cloud intensities (i.e., counts) at locations along a spatial transect at time points. The data roughly follow an overdispersed Poisson distribution, and so we assumed the model from Example 3. We randomly set aside 10% of the observations as test data, so that and was a subset of the identity. The observation locations were different for each time period. The evolution was taken to be linear, , where was a tridiagonal matrix with on the main diagonal, and and on the first diagonals above and below, respectively. The evolution error was assumed to exhibit spatial dependence in the form of a Matérn covariance with smoothness 1.5, such that the th element of is given by , where . The unknown parameters were . For , we set , , and . Note that our model addresses the potential extension to time-varying parameters suggested in Wikle, (2002).
We considered two methods for filtering inference on the state and parameters: our PEnKF from Example 8, and an exact approximation of a Rao-Blackwellized particle filter (EARBPF; Johansen et al., , 2012). Forecast independence does not hold for all parameters in this example, and so we used ensembles of size . For tapering of the forecast covariance matrix (see Section 2.2), we chose a Wendland taper with range 8. For the EARBPF, which is similar to the PEnKF except that it uses local particle filters instead of the EnKF and the Laplace approximation to integrate out the state for each parameter particle, we considered one setting with a similar computation time as the PEnKF ( and ), and one setting with a much larger computation time ( and ) that should provide a close approximation to and will henceforth be referred to as the true filtering distribution.
The results, averaged over ten different initial ensembles , are shown in Table 3. In terms of filtering predictions of the test data, the PEnKF performed almost as well as the true filtering distributions, and considerably better than the EARBPF. The PEnKF approximation to the filtering distributions of the six parameters in is also closer than the EARBPF to the true distributions in terms of earth-mover’s distance (computed using the R package by Urbanek and Rubner, , 2012) averaged over time. As an example, Figure 5(b) shows the filtering distributions for one initial ensemble over time of the parameter , which can be viewed as quantifying the amount of movement of the cloud intensities “one pixel to the right” from one time point to the next. This rightward drift, which is clearly present in the data (see Figure 5(a)), is picked up fairly quickly by the true and PEnKF filtering distributions, while the EARBPF approximation even implies a negative between and .
5.3 Smoothing inference for the Lorenz-96 model
Finally, we considered smoothing inference for simulated data based on the popular nonlinear Lorenz-96 model (Lorenz, , 1996)
described in Example 4, with , where is the covariance matrix obtained by a long run of the model. Starting with and , we simulated data from the model for with , assuming with and . The Lorenz-96 equations were solved numerically using an Euler scheme.
We simulated 100 such datasets, with the true drawn for each simulation from the prior. For each dataset, we ran three different methods to obtain approximations to the smoothing distribution using a starting value of . The first method considered was the EnKS in Algorithm 2.3 with 1,000 and state augmentation; that is, was treated as part of the state vector with artificial evolution distribution . We also considered the GEnKS in Example 9. For both EnKS algorithms with used a taper over time and space given by a Wendland correlation function with radius 3 and 8, respectively. Finally, we considered the particle Gibbs sampler (Andrieu et al., , 2010, Sect. 2.4.3), which uses the same distribution as the EnKS in (22) to update , but uses a “conditional” particle filter to sample from the FCD of the state given the data and . For both the GEnKS and the particle Gibbs sampler, we used an ensemble size of and 100 Gibbs iterations, the first 20 of which were considered burn-in.
Due to the nonlinear evolution, it is not possible to obtain the true smoothing distribution of the states and parameter in any reasonable computation time. Hence, we compared the three methods using the same proper scoring rules as in Section 5.1, averaged over the 100 simulated datasets. The results are summarized in Table 4. Our GEnKS algorithm performed best both in terms of inference on the parameter and on the states. The state augmentation method failed, leading to very high (i.e., bad) scores for the posterior on . The particle Gibbs sampler produced worse inference on than simply using the prior distribution (i.e., completely ignoring any information in the data). Some of the results are also summarized in Figure 6, which shows the collapse of the particle Gibbs procedure for inference on and its non-convergence in terms of inference on .
Note that while we considered to be the only unknown parameter here for clearness of exposition, our GEnKS should also work well for unknown parameters in and scalar multiplicative parameters in and , in that closed-form FCDs are available in both of these cases if a normal and inverse-gamma priors, respectively, are assumed.
6 Conclusions
We have introduced a new class of ensemble filtering and smoothing algorithms for Bayesian inference in high-dimensional dynamic spatio-temporal models. The algorithms combine existing approaches from the state-space literature with ensemble Kalman methods from geophysics. The main advantage of our algorithms is that they scale well to high-dimensional models, and they allow for nonlinearity, non-Gaussian observations, and unknown parameters. In contrast, most existing state-space algorithms break down in high dimensions, and most ensemble Kalman approaches are poorly suited for non-Gaussian models with unknown parameters. We illustrated the proposed methods with four examples: two non-Gaussian models, a nonlinear geophysical system, and a real data example of cloud motion. We showed that our methods outperform the main competing approaches, including EnKFs, particle filters, and particle MCMC. Our methods also provide an approach for approximate Bayesian inference for large spatial data, in the special case of only time point.
We have proposed two types of algorithms: Ensemble Gibbs and Ensemble Marginalization methods. For Ensemble Gibbs methods, the EnKF is used for approximate simulation from the full conditional state distribution. In our examples, we show that these approximate Gibbs samplers can provide accurate state and parameter inference and converge within a small number of iterations. For Ensemble Marginalization methods, the EnKF is used to approximately integrate out the states and simulate from their full conditional posterior. In this context, we have studied the theoretical properties of EnKF-based likelihood estimates and shown that an ensemble size that is linear in the effective dimension is sufficient to bound the variance of the estimated loglikelihood. This is in contrast to particle-filter-based likelihood estimates, for which the required number of particles grows exponentially in the effective dimension. We leave as future research further study on the theoretical properties of the proposed algorithms.
Statisticians now acknowledge that fast, approximate methods are necessary to analyze very large datasets. The EnKF, developed in the geophysics community, has been successfully applied for fast, approximate inference on models with millions of state variables. We believe that combining standard statistical approaches with these approximate but scalable EnKF methods provides solutions to problems that might be otherwise unsolvable.
Acknowledgments
Katzfuss’ research was partially supported by National Science Foundation (NSF) Grant DMS–1521676 and NSF CAREER Grant DMS–1654083. Wikle acknowledges the support of NSF grant SES-1132031, funded through the National Science Foundation Census Research Network program. We would like to thank Fuqing Zhang, Amir Nikooienejad, the Associate Editor, and two anonymous referees for helpful comments and suggestions.
Appendix A Proofs
Proof of Proposition 1.
For notational simplicity, we drop the on , and we assume without loss of generality. First, note that follows a noncentral -distribution, and so
[TABLE]
where and . Using properties of the lognormal distribution, we have
[TABLE]
and finally, using the delta-method, we have the asymptotic variance
[TABLE]
∎
Proof of Proposition 2.
We drop the on , and assume without loss of generality. Due to the diagonal tapering, the covariance matrix in the likelihood is a diagonal matrix with th element
[TABLE]
where and . Because , the variance of the mean estimator is negligible relative to that of , and we regard as a constant. Thus, by the delta method,
[TABLE]
Finally, we have , and so . ∎
Proof of Proposition 3.
For linear Gaussian SSMs, the EnKF converges to the Kalman filter as (e.g., Butala et al., , 2008; Mandel et al., , 2011), as long as is a consistent estimator of the forecast covariance matrix. As the EnKS can be viewed as an EnKF on an extended state-space, the same result holds for the EnKS. As stated in Section 3.2.2, the EnKF likelihood also converges to the true likelihood at each time . Combining this with convergence results for MCMC (e.g., Robert and Casella, , 2004) and SMC (e.g., Doucet et al., , 2001) methods under mild regularity conditions, it is clear that the PEnKF, GEnKS, MHEnKS, and PEnKS all converge to the true filtering or smoothing distributions in a HSSM with linear evolution operators , as and or the number of MCMC iterations increase. ∎
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Alquier et al., (2016) Alquier, P., Friel, N., Everitt, R., and Boland, A. (2016). Noisy Monte Carlo: Convergence of Markov chains with approximate transition kernels. Statistics and Computing , 26(1-2):29–47.
- 2Amezcua and Van Leeuwen, (2014) Amezcua, J. and Van Leeuwen, P. J. (2014). Gaussian anamorphosis in the analysis step of the En KF: a joint state-variable/observation approach. Tellus A , 66:1–18.
- 3Anderson et al., (2009) Anderson, J., Hoar, T., Raeder, K., Liu, H., Collins, N., Torn, R., and Avellano, A. (2009). The data assimilation research testbed: A community facility. Bulletin of the American Meteorological Society , 90(9):1283–1296.
- 4Anderson, (2001) Anderson, J. L. (2001). An ensemble adjustment Kalman filter for data assimilation. Monthly Weather Review , 129:2884–2903.
- 5Anderson, (2003) Anderson, J. L. (2003). A local least squares framework for ensemble filtering. Monthly Weather Review , 131:634–642.
- 6(6) Anderson, J. L. (2007 a). An adaptive covariance inflation error correction algorithm for ensemble filters. Tellus A , 59:210–224.
- 7(7) Anderson, J. L. (2007 b). Exploring the need for localization in ensemble data assimilation using an hierarchical ensemble filter. Physica D , 230:99–111.
- 8Anderson, (2009) Anderson, J. L. (2009). Ensemble Kalman filters for large geophysical applications. IEEE Control Systems Magazine , 29(June):66–82.