Ensemble updating of binary state vectors by maximising the expected number of unchanged components

TL;DR

This paper introduces a novel ensemble updating method for binary state vectors that maximizes the expected number of unchanged components, improving filtering accuracy over naive approaches in binary state estimation tasks.

Contribution

The authors adapt ensemble filtering techniques to binary variables using Markov chains, proposing a method that maximizes component stability during updates, which is a novel approach in this context.

Findings

Reduces the Frobenius norm of estimation error by half compared to naive methods.

Demonstrates improved accuracy in a petroleum reservoir simulation example.

Method can be generalized to more complex state spaces.

Abstract

In recent years, several ensemble-based filtering methods have been proposed and studied. The main challenge in such procedures is the updating of a prior ensemble to a posterior ensemble at every step of the filtering recursions. In the famous ensemble Kalman filter, the assumption of a linear-Gaussian state space model is introduced in order to overcome this issue, and the prior ensemble is updated with a linear shift closely related to the traditional Kalman filter equations. In the current article, we consider how the ideas underlying the ensemble Kalman filter can be applied when the components of the state vectors are binary variables. While the ensemble Kalman filter relies on Gaussian approximations of the forecast and filtering distributions, we instead use first order Markov chains. To update the prior ensemble, we simulate samples from a distribution constructed such that the expected number of equal components in a prior and posterior state vector is maximised. We demonstrate the performance of our approach in a simulation example inspired by the movement of oil and water in a petroleum reservoir, where also a more na\"{i}ve updating approach is applied for comparison. Here, we observe that the Frobenius norm of the difference between the estimated and the true marginal filtering probabilities is reduced to the half with our method compared to the na\"{i}ve approach, indicating that our method is superior. Finally, we discuss how our methodology can be generalised from the binary setting to more complicated situations.

Figures

Tables

Table 1. (a)

$k$	1
$t_k^{\ast }$	0.400000
$q_k^{\ast 0}\left(t_k^{\ast }\right)$	1.000000
$q_k^{\ast 1}\left(t_k^{\ast }\right)$	0.211299

$k$	1
$t_k^{\ast }$	0.400000
$q_k^{\ast 0}\left(t_k^{\ast }\right)$	1.000000
$q_k^{\ast 1}\left(t_k^{\ast }\right)$	0.211299

Table 3. (a)

$k$	2	3	4
$t_k^{\ast }$	0.305356	0.308676	0.281108
$q_k^{\ast 00}\left(t_k^{\ast }\right)$	1.000000	1.000000	0.853968
$q_k^{\ast 01}\left(t_k^{\ast }\right)$	0.481489	0.212926	0.000000
$q_k^{\ast 10}\left(t_k^{\ast }\right)$	1.000000	0.860986	0.546043
$q_k^{\ast 11}\left(t_k^{\ast }\right)$	0.097118	0.000000	0.000000

Table 4. Table 2: Probabilities defining the true model p ​ ( x t | x t − 1 ) 𝑝 conditional superscript 𝑥 𝑡 superscript 𝑥 𝑡 1 p(x^{t}|x^{t-1}) used to simulate a true chain { x t } t = 1 T superscript subscript superscript 𝑥 𝑡 𝑡 1 𝑇 \{x^{t}\}_{t=1}^{T} .

$x_{i-1}^{t-1}$	$x_i^{t-1}$	$x_{i+1}^{t-1}$	$P(x_i^t=1|x_{i-1}^t=1,x_{i-1:i+1}^{t-1})$	$P(x_i^t=1|x_{i-1}^t=0,x_{i-1:i+1}^{t-1})$
0	0	0	0.0100	0.0050
1	0	0	0.0400	0.0100
0	1	0	0.9999	0.9800
1	1	0	0.9999	0.9900
0	0	1	0.0400	0.0400
1	0	1	0.9800	0.0400
0	1	1	0.9999	0.9800
1	1	1	0.9999	0.9800

Equations

$$x^1\sim p(x^1),$$

$$x^t|x^{t-1}\sim p(x^t|x^{t-1}),t>1,$$

$$y^t|x^t\sim p(y^t|x^t),t\ge 1.$$

$$\text{i)}p(x^t|y^{1:t-1})={\int }_{{\Omega }_X}p(x^t|x^{t-1})p(x^{t-1}|y^{1:t-1})\text{d}{x}^{t-1},$$

$$\text{ii)}p(x^t|y^{1:t})=\frac{p(x^t|y^{1:t-1})p(y^t|x^t)}{{\int }_{{\Omega }_X}p(x^t|y^{1:t-1})p(y^t|x^t)\text{d}{x}^t}.$$

$$x^1\sim \mathcal{N}(x^1|{\mu }^1,Q^1),$$

$$x^t|x^{t-1}=A^t{x}^{t-1}+{ϵ}^t,{ϵ}^t\sim {\mathcal{N}}_n(ϵ|0,Q^t),$$

$$y^t|x^t=H^t{x}^t+{\omega }^t,{\omega }^t\sim {\mathcal{N}}_k(\omega |0,R^t),$$

$$y|x=Hx+ϵ,ϵ\sim {\mathcal{N}}_k(ϵ;0,R).$$

$$\tilde{\mu }=\mu +K(y-H\mu ),$$

$$\tilde{\Sigma }=(I-KH)\Sigma ,$$

$$\tilde{x}=x+K(y-Hx+ϵ),$$

$${\tilde{x}}^{(i)}=x^{(i)}+\hat{K}(y-H{x}^{(i)}+{ϵ}^{(i)}),i=1,\dots ,M.$$

$${\tilde{x}}^{(i)}=\hat{\mu }+\hat{K}(y-H\hat{\mu })+B(x^{(i)}-\hat{\mu }),$$

$$B\hat{\Sigma }{B}^{\prime }=(I-\hat{K}H)\hat{\Sigma }.$$

$$f(\tilde{x}|y)={\int }_{{\Omega }_X}f(x)q(\tilde{x}|x,y)\text{d}x.$$

$${\tilde{x}}^{(i)}\sim q(\tilde{x}|x^{(i)},y),$$

$$q(\tilde{x}|x,y)=\mathcal{N}\left(\tilde{x};x+\hat{K}(y-Hx),\hat{K}R{\hat{K}}^{\prime }\right).$$

$$q(\tilde{x}|x,y)=\delta \left(\tilde{x};\hat{\mu }+\hat{K}(y-H\hat{\mu })+B(x-\hat{\mu })\right).$$

$$f(x)=f(x_1)f(x_2|x_1)\cdots f(x_n|x_{n-1}).$$

$$p(y|x)=p(y_1|x_1)\cdots p(y_n|x_n).$$

$$f(\tilde{x}|y)=\sum _{x\in {\Omega }_X}f(x)q(\tilde{x}|x,y).$$

$$q^{\ast }(\tilde{x}|x,y)=\underset{q(\tilde{x}|x,y)}{\mathrm{argmax}}{\text{E}}_{\pi }\left[\sum _{i=1}^{n}1(x_i={\tilde{x}}_i)\right],$$

$$\pi (\tilde{x},x|y)=f(x)q(\tilde{x}|x,y),$$

$$q(\tilde{x}|x,y)=q({\tilde{x}}_1|x_1,y)q({\tilde{x}}_2|{\tilde{x}}_1,x_2,y)q({\tilde{x}}_3|{\tilde{x}}_2,x_3,y)\cdots q({\tilde{x}}_n|{\tilde{x}}_{n-1},x_n,y).$$

$$q^{\ast }(\tilde{x}|x,y)=q^{\ast }({\tilde{x}}_1|x_1,y)q^{\ast }({\tilde{x}}_2|{\tilde{x}}_1,x_2,y)q^{\ast }({\tilde{x}}_3|{\tilde{x}}_2,x_3,y)\cdots q^{\ast }({\tilde{x}}_n|{\tilde{x}}_{n-1},x_n,y).$$

$${\pi }_k$$

$$q_k$$

$$\underset{q_{k:n}}{\max }{\text{E}}_{\pi }\left[\sum _{i=k}^{n}1(x_i={\tilde{x}}_i)\right]$$

$$\underset{q_{(k-1):n}}{\max }{\text{E}}_{\pi }\left[\sum _{i=k-1}^{n}1(x_i={\tilde{x}}_i)\right]$$

$${\pi }_{t_k}({\tilde{x}}_{k-1}=0,x_k=0|y)$$

$${\pi }_{t_k}({\tilde{x}}_{k-1}=0,x_k=1|y)$$

$${\pi }_{t_k}({\tilde{x}}_{k-1}=1,x_k=0|y)$$

Full text

1 Introduction

A state-space model consists of a latent $\{x^{t}\}_{t=1}^{\infty}$ process and an observed $\{y^{t}\}_{t=1}^{\infty}$ process, where $y^{t}$ is a partial observation of $x^{t}$. More specifically, the $y^{t}$’s are assumed to be conditionally independent given the $x^{t}$ process and $y^{t}$ only depends on $x^{t}$. Estimation of the latent variable at time $t$, $x^{t}$, given all observations up to this time, $y^{1:t}=(y^{1},\ldots,y^{t})$, is known as the filtering or data assimilation problem. In the linear Gaussian situation an easy to compute and exact solution is available by the famous Kalman filter. In most non-linear or non-Gaussian situations no computationally feasible exact solution exists and ensemble methods are therefore frequently adopted. The distribution $p(x^{t}|y^{1:t})$ is then not analytically available, but is represented by a set of realisations $\tilde{x}^{t(1)},\ldots,\tilde{x}^{t(M)}$ from this filtering distribution. Assuming such an ensemble of realisations to be available for time $t-1$, the filtering problem is solved for time $t$ in two steps. First, based on the Markov chain model for the $x^{t}$ process, each $\tilde{x}^{t-1(i)}$ is used to simulate a corresponding forecast realisation $x^{t(i)}$, which marginally are independent samples from $p(x^{t}|y^{1:t-1})$. This is known as the forecast or prediction step. Second, an update step is performed, where each $x^{t(i)}$ is updated to take into account the new observation $y^{t}$ and the result is an updated ensemble $\tilde{x}^{t(1)},\ldots,\tilde{x}^{t(M)}$ which represents the filtering distribution at time $t$, $p(x^{t}|y^{1:t})$. The updating step is the difficult one and the different strategies that have been proposed can be classified into two classes, particle filters and ensemble Kalman filters.

In particle filters (Doucet et al., 2001) each filtering realisation $\tilde{x}^{t(i)}$ comes with an associated weight $\tilde{w}^{t(i)}$, and the pair $(\tilde{w}^{t(i)},\tilde{x}^{t(i)})$ is called a particle. In the forecast step a forecast particle $(w^{t(i)},x^{t(i)})$ is generated from each filtering particle $(\tilde{w}^{t-1(i)},\tilde{x}^{t-1(i)})$ by generating $x^{t(i)}$ from $\tilde{x}^{t-1(i)}$ as discussed above and by keeping the weight unchanged, i.e. $w^{t(i)}=\tilde{w}^{t-1(i)}$. The updating step consists of two parts. First the weights are updated by multiplying each forecast weight $w^{t(i)}$ by the associated likelihood value $p(y^{t}|x^{t(i)})$, keeping the $x^{t}$ component of the particles unchanged. Thereafter a re-sampling may be performed, where $(\tilde{w}^{t(i)},\tilde{x}^{t(i)}),i=1,\ldots,M$ are generated by sampling the $\tilde{x}^{t(i)}$’s independently from $x^{t(i)},i=1,\ldots,M$ with probabilities proportional to the updated weights, and thereafter setting all the new filtering weights $\tilde{w}^{t(i)}$ equal to one. Different criteria can be used to decide whether or not the re-sampling should be done. The particle filter is very general in that it can be formulated for any Markov $x^{t}$ process and any observation distribution $p(y^{t}|x^{t})$. However, when running the particle filter one quite often ends up with particle depletion, meaning that a significant fraction of the particles ends up with negligible weights, which in practice requires the number of particles to grow exponentially with the dimension of the state vector $x^{t}$. To cope with the particle depletion problem various modifications of the basic particle filter approach described here have been proposed, see for example Carpenter et al. (1999), Gilks and Berzuini (2001) and Doucet and Johansen (2011) and references therein.

The ensemble Kalman filter (Evensen, 1994; Burgers et al., 1998) uses approximations in the update step, and thereby produces only an approximate solution to the filtering problem. In the update step it starts by using the forecast samples $x^{t(i)},i=1,\ldots,M$ to estimate a Gaussian approximation to the forecast distribution $p(x^{t}|y^{1:t-1})$. This is combined with an assumed Gaussian observation distribution $p(y^{t}|x^{t})$ to obtain a Gaussian approximation to the filtering distribution $p(x^{t}|y^{1:t})$. Based on this Gaussian approximation the filtering ensemble is generated by sampling $\tilde{x}^{t(i)}$ independently from Gaussian distributions, where the mean of $\tilde{x}^{t(i)}$ equals $x^{t(i)}$ plus a shift which depends on the approximate Gaussian filtering distribution. The associated variance is chosen so that the marginal distribution of the generated filtering sample $\tilde{x}^{t(i)}$ is equal to the Gaussian approximation to $p(x^{t}|x^{1:t})$ when the forecast sample $x^{t(i)}$ is assumed to be distributed according to the Gaussian approximation to $p(x^{t}|y^{1:t-1})$. The basic ensemble Kalman filter described here is known to have a tendency to under-estimate the variance in the filtering distribution and various remedies have been proposed to correct for this, see for example the discussions in Anderson (2007a, b) and Sætrom and Omre (2013). The square root filter (Whitaker and Hamill, 2002; Tippett et al., 2003) is a special variant of the ensemble Kalman filter where the update step is deterministic. The filtering ensemble is then generated from the forecast ensemble only by adding a shift to each ensemble element. Here the size of the shift is chosen to get that the marginal distribution of the filtering realisations is equal to the approximated Gaussian filtering distribution.

The Gaussian approximations used in the ensemble Kalman filter update step limits the use of this filter type to continuous variables, whereas the particle filter setup can be used for both continuous and categorical variables. In the literature there exists a few attempts to use the ensemble Kalman filter setup also for categorical variables, see in particular Oliver et al. (2011). The strategy then used for the update step is first to map the categorical variables over to continuous variables, perform the update step as before in the continuous space, and finally map the updated continuous variables back to corresponding categorical variables. Our goal in this article is to study how the basic ensemble Kalman filter idea can be used for categorical variables without having to map the categorical variables over to a continuous space. As discussed above the update step is the difficult one in ensemble filtering methods. The basic ensemble Kalman filter update starts by estimating a Gaussian approximation to the forecast distribution $p(x^{t}|y^{1:t-1})$. More generally one may use another parametric class than the Gaussian. For categorical variables the simplest alternative is to consider a first-order Markov chain, which is what we focus on in this article. Having a computationally feasible approximation for the forecast distribution we can find a corresponding approximate filtering distribution. Given the forecast ensemble the question then is from which distribution to simulate the filtering ensemble to obtain that the filtering realisations marginally is distributed according to the given approximate filtering distribution, corresponding to the property for the standard ensemble Kalman filter. In this article we develop in detail one possible way to do this when the elements of the state vector are binary, the approximate forecast distribution is a first-order Markov chain, and the observation distribution has a specifically simple form.

The article has the following layout. First, in Section 2 we specify the general state-space model and the associated filtering problem and present the ensemble Kalman filter. Next, in Section 3 we formulate the general ensemble updating strategy briefly discussed above. Then, in Section 4, we limit the attention to a situation where the elements of the state vector are binary and discuss how the update step can be performed in this case. After that, we present two numerical experiments with simulated data in Section 5. Finally, in Section 6 we give some final remarks and briefly discuss how our setup can be generalised to a situation with more than two classes and an assumed higher-order Markov chain model for the forecast distribution.

2 Preliminaries

In this section, we review some basic theoretical aspects of ensemble-based filtering methods. The material presented should provide the reader with the necessary background for understanding our proposed approach and it also establishes some of the notations used throughout the article.

2.1 Review of the filtering problem

The filtering problem in statistics can be nicely illustrated with a graphical model; see Figure 1. Here, $\{x^{t}\}_{t=1}^{\infty}$ represents a time series of unobserved states and $\{y^{t}\}_{t=1}^{\infty}$ a corresponding time series of observations. Each state $x^{t}$ is $n$-dimensional and can take on values in a state space $\Omega_{X}$, while each observation $y^{t}$ is $k$-dimensional and can take on values in a state space $\Omega_{Y}$. The series of unobserved states, called the state process, constitutes a first order Markov chain with initial state $x^{1}$, initial distribution $p(x^{1})$, and transition probabilities $p(x^{t}|x^{t-1}),t>1$. For each state $x^{t},t\geq 1$, there is a corresponding observation $y^{t}$. The observations are assumed conditionally independent given the state process, with $y^{t}$ depending on $\{x^{t}\}_{t=1}^{\infty}$ only through $x^{t}$, according to some likelihood model $p(y^{t}|x^{t})$. To summarise, the model is specified by

[TABLE]

The objective of the filtering problem is, for each $t$, to compute the so-called filtering distribution, $p(x^{t}|y^{1:t})$, that is, the distribution of $x^{t}$ given all observations up to this time, $y^{1:t}=(y^{1},\dots,y^{t})$. Because of the assumptions about the state and observation processes, it can be shown (see e.g. Künsch (2000)) that the series of filtering distributions can be computed recursively according to the following equations:

[TABLE]

As one can see, the recursions evolve as a two-step process, each iteration consisting of i) a prediction step, and ii) an update step. In the prediction, or forecast step, one computes the predictive, or forecast, distribution $p(x^{t}|y^{1:t-1})$, while in the update step, one computes the filtering distribution $p(x^{t}|y^{1:t})$ by conditioning the predictive distribution on the incoming observation $y^{t}$ through application of Bayes’ rule. The update step can be formulated as a standard Bayesian inference problem, with $p(x^{t}|y^{1:t-1})$ becoming the prior, $p(y^{t}|x^{t})$ the likelihood, and $p(x^{t}|y^{1:t})$ the posterior.

There are two important special cases where the analytical solutions to the filtering recursions, (2.2a) and (2.2b), can be computed exactly. The first case is the hidden Markov model (HMM). Here, the state space $\Omega_{X}$ consists of a finite number of states, and the integrals in (2.2a) and (2.2b) reduce to finite sums. If the number of states in $\Omega_{X}$ is large, however, the summations become computer-intensive, rendering the filtering recursions computationally intractable. The second case is the linear Gaussian state space model, which can be formulated as follows:

[TABLE]

[TABLE]

[TABLE]

where $A^{t}\in\mathbb{R}^{n\times n}$ and $H^{t}\in\mathbb{R}^{k\times n}$ are non-random linear operators, $Q^{t}\in\mathbb{R}^{n\times n}$ and $R^{t}\in\mathbb{R}^{k\times k}$ are covariance matrices, and $x^{1},\epsilon^{1},\epsilon^{2},\dots,\omega^{1},\omega^{2},\dots$ are all independent. In this case, the sequence of predictive and filtering distributions are Gaussian, and the filtering recursions lead to the famous Kalman filter (Kalman, 1960).

In general, we are unable to evaluate the integrals in (2.2a) and (2.2b), leaving the filtering recursions intractable. Approximate solutions therefore become necessary. The most common approach in this regard is the class of ensemble-based methods, where a set of samples, called an ensemble, is used to empirically represent the sequence of predictive and filtering distributions. Starting from an initial ensemble $\{x^{1(1)},\dots,x^{1(M)}\}$ of $M$ independent realisations from the Markov chain initial model $p(x^{1})$, the idea is to advance this ensemble forward in time according to the model dynamics. As the original filtering recursions, the propagation of the ensemble alternate between an update step and a prediction step. Specifically, suppose at time $t\geq 1$ that an ensemble $\{x^{t(1)},\dots,x^{t(M)}\}$ of independent realisations from the forecast distribution $p(x^{t}|y^{1:t-1})$ is available. We then want to update this forecast ensemble by conditioning on the incoming observation $y^{t}$ in order to obtain an updated, or posterior, ensemble $\{\tilde{x}^{t(1)},\dots,\tilde{x}^{t(M)}\}$ with independent realisations from the filtering distribution $p(x^{t}|y^{1:t})$. If we are able to carry out this updating, we proceed and propagate the updated ensemble $\{\tilde{x}^{t(1)},\dots,\tilde{x}^{t(M)}\}$ one step forward in time by using the Markov chain prior model to simulate $x^{t+1(i)}|\tilde{x}^{t(i)}\sim p(x^{t+1}|\tilde{x}^{t(i)})$ for each $i$. This produces a new forecast ensemble for time $t+1$ with independent realisations from the forecast distribution $p(x^{t+1}|y^{1:t})$. However, while we are often able to simulate from the Markov forward model $p(x^{t}|x^{t-1})$, there is no straightforward way for carrying out the update of the prior ensemble $\{\tilde{x}^{t(1)},\dots,\tilde{x}^{t(M)}\}$ to a posterior ensemble $\{\tilde{x}^{t(1)},\dots,\tilde{x}^{t(M)}\}$. Therefore, ensemble methods require approximations in the update step. Consequently, the assumption we make at the beginning of each time step $t$, i.e. that $x^{t(1)},\dots,x^{t(M)}$ are exact and independent realisations from $p(x^{t}|y^{1:t-1})$, holds only approximately, except in the initial time step.

In the remains of this article, we focus primarily on the challenging updating of a prior ensemble $\{x^{t(1)},\dots,x^{t(M)}\}$ into a posterior ensemble $\{\tilde{x}^{t(1)},\dots,\tilde{x}^{t(M)}\}$ at a specific time step $t$. We refer to this task as the ensemble updating problem. For simplicity, we omit the time superscript $t$ and the $y^{1:t-1}$ from the notations in the remaining sections, as these quantities are fixed throughout. That is, we write $x$ instead of $x^{t}$, $p(x)$ instead of $p(x^{t}|y^{1:t-1})$, $p(x|y)$ instead of $p(x^{t}|y^{1:t})$, and so on.

2.2 The ensemble Kalman filter

The ensemble Kalman filter (EnKF), first introduced in the geophysics literature by Evensen (1994), is an approximate ensemble-based method that relies on Gaussian approximations to overcome the difficult updating of the prior ensemble. The updating is done in terms of a linear shift of each ensemble member, closely related to the traditional Kalman filter equations. The literature on the EnKF is extensive, but some basic references include Evensen (1994) and Burgers et al. (1998). Here, we only provide a brief presentation. For simplicity reasons, we restrict the focus to linear, Gaussian observational models; that is, the relationship between $y$ and $x$ is assumed linear with additive zero-mean Gaussian noise,

[TABLE]

There exist two main classes of EnKFs, stochastic filters and deterministic, or so-called square root filters, differing in whether the updating of the ensemble is carried out stochastically or deterministically. The stochastic EnKF is probably the most common version, and we begin our below presentation of the EnKF by focusing on this method.

Consider first a linear Gaussian state space model as introduced in the previous section. Under this linear Gaussian model, it follows from the Kalman filter recursions that the current predictive, or prior model $p(x)$ is a Gaussian, $\mathcal{N}_{n}(x;\mu,\Sigma)$, with analytically tractable mean $\mu$ and covariance $\Sigma$. Furthermore, the current filtering, or posterior distribution $p(x|y)$ is a Gaussian, $\mathcal{N}_{n}(x;\tilde{\mu},\tilde{\Sigma})$, with mean $\tilde{\mu}$ and covariance $\tilde{\Sigma}$ analytically available from the Kalman filter update formulas as

[TABLE]

[TABLE]

where $K=\Sigma H^{\prime}(H\Sigma H^{\prime}+R)^{-1}$ is the Kalman gain. The stochastic EnKF update is based on the following fact: If $x\sim\mathcal{N}_{n}(x;\mu,\Sigma)$ and $\epsilon\sim\mathcal{N}_{k}(\epsilon;0,R)$ are independent random samples, then

[TABLE]

is a random sample from $\mathcal{N}_{n}(x;\tilde{\mu},\tilde{\Sigma})$. The verification of this result is straightforward. Clearly, under the assumption that the prior ensemble $\{x^{(1)},\dots,x^{(M)}\}$ contains independent samples from the Gaussian $\mathcal{N}_{n}(x;\mu,\Sigma)$, one theoretically valid way to obtain the updated ensemble is to plug each $x^{(i)}$ into (2.4). However, for computational reasons, the prior covariance matrix $\Sigma$ and hence the Kalman gain $K$ can be intractable for high-dimensional systems. To circumvent this issue, stochastic EnKF passes each prior sample $x^{(i)}$ through a linear shift identical to (2.4), but with the true Kalman gain $K$ replaced with an empirical estimate $\hat{K}$ inferred from the prior ensemble,

[TABLE]

In the EnKF literature, each term $Hx^{(i)}-\epsilon^{(i)}$ is typically referred to as a ”perturbed” observation. Under the Gaussian-linear assumptions, the update (2.5) produces approximate samples from the Gaussian posterior model $\mathcal{N}_{n}(x;\tilde{\mu},\tilde{\Sigma})$ that corresponds to the Gaussian prior model $\mathcal{N}_{n}(x;\mu,\Sigma)$. The update is consistent in the sense that as the ensemble size goes to infinity, the distribution of the updated samples approaches the Gaussian $\mathcal{N}_{n}(x;\tilde{\mu},\tilde{\Sigma})$, that is, the solution of the Kalman filter.

Although the EnKF update is based on Gaussian assumptions about the predictive and filtering distributions, it can be applied even if these assumptions are not met. Naturally, bias is in this case introduced, and the updated samples will not converge in distribution to the true posterior $p(x|y)$. However, since the update is a linear combination of the $x^{(i)}$’s, non-Gaussian properties present in the true prior and posterior models can, to some extent, be captured.

Deterministic EnKFs instead use a non-random linear transformation to update the ensemble. In the following, let $\hat{\mu}$ and $\hat{\Sigma}$ denote the sample mean and sample covariance, respectively, of the prior ensemble. Further, let $\hat{\tilde{\mu}}$ and $\hat{\tilde{\Sigma}}$ denote the mean and covariance, respectively, of the Gaussian posterior model $\mathcal{N}_{n}(x;\hat{\tilde{\mu}},\hat{\tilde{\Sigma}})$ corresponding to the Gaussian prior approximation $\mathcal{N}_{n}(x;\hat{\mu},\hat{\Sigma})$. Generally, the update equation of a square root EnKF can be written as

[TABLE]

where $B\in\mathbb{R}^{n\times n}$ is a solution to the quadratic matrix equation

[TABLE]

Note that $B$ is not unique except in the univariate case. This gives rise to a variety of square root algorithms; see Tippett et al. (2003) for further explanation. As such, several square root formulations have been proposed in the literature, including, but not limited to, Whitaker and Hamill (2002), Anderson (2001), and Bishop et al. (2001). The non-random update of square root EnKF ensures that the sample mean and sample covariance of the posterior ensemble equal $\hat{\tilde{\mu}}$ and $\hat{\tilde{\Sigma}}$ exactly. This is different from stochastic EnKFs where, under linear-Gaussian assumptions, the sample mean and sample covariance of the posterior ensemble only equal $\hat{\tilde{\mu}}$ and $\hat{\tilde{\Sigma}}$ in expectation.

3 A general ensemble updating framework

In this section, we present a general framework for ensemble updating. Both the EnKF and the update strategy for binary vectors proposed in this article can be viewed as special cases within the framework.

3.1 The framework

For convenience, let us first quickly review the ensemble updating problem. Starting out, we have a prior ensemble, $\{x^{(1)},\dots,x^{(M)}\}$, which we assume contain independent realisations from a prior model $p(x)$. The prior model $p(x)$ is typically intractable in this context, either computationally or analytically, or both. Given an incoming observation $y$ and a corresponding likelihood model $p(y|x)$ the goal is to update the prior ensemble according to Bayes’ rule in order to obtain a posterior ensemble, $\{\tilde{x}^{(1)},\dots,\tilde{x}^{(M)}\}$, with independent realisations from the posterior model $p(x|y)$. However, carrying out this update exactly is generally unfeasible and approximate strategies are required.

Conceptually, the proposed framework is quite simple. It involves three main steps as follows. First, we replace the intractable model, $p(x|y)\propto p(x)p(y|x)$, with a simpler model, $f(x|y)\propto f(x)p(y|x)$, using information from the prior ensemble to construct the approximate prior model $f(x)$. Throughout, we will refer to the model $f(x|y)\propto f(x)p(y|x)$ as the assumed model. Notice that the likelihood model $p(y|x)$ has not been replaced; for simplicity, we assume that this model is already tractable. Second, we put forward a distribution conditional on $x$ and $y$, denoted $q(\tilde{x}|x,y)$, obeying the following property:

[TABLE]

Third, we update the prior ensemble by generating samples from this conditional distribution,

[TABLE]

To understand the framework, one must note that under the assumption that the assumed model is the correct one, the prior samples have distribution $f(x)$ and the updated samples should have distribution $f(x|y)$. If one is able to compute and sample from $f(x|y)$, one straightforward way to obtain the updated samples is to sample directly from $f(x|y)$. However, since the assumed model is not really the correct one, that is probably not the best way to proceed. The prior ensemble contains valuable information about the true model that may not have been captured by the assumed model, and by straightforward simulation from $f(x|y)$ this information is lost. To capture more information from the prior ensemble, it is probably advantageous to simulate conditionally on the prior samples. This is why we introduce the conditional distribution $q(\tilde{x}|x,y)$. The criterion in (3.1) ensures that the marginal distribution of the updated samples $\tilde{x}^{(1)},\dots,\tilde{x}^{(M)}$ produced by $q(\tilde{x}|x,y)$ still would be $f(x|y)$ given that the assumed model is the correct one. However, since the assumed model is not the correct one, the marginal distribution of the updated samples is not $f(x|y)$, but some other distribution, hopefully one closer to the true posterior model $p(x|y)$.

There are two especially important things about the proposed framework that must be taken care of in a practical application. First, we need to select an assumed prior $f(x)$ which, combined with the likelihood model $p(y|x)$, returns a tractable posterior $f(x|y)$. Second, we need to construct the updating distribution $q(\tilde{x}|x,y)$. Clearly, there can be many, or infinitely many different $q(\tilde{x}|x,y)$ which all fulfil the constraint in (3.1). Below, we consider the proposed framework in two specific situations. The first case corresponds to the EnKF where $f(x)$, $p(y|x)$ and $q(\tilde{x}|x,y)$ are all Gaussians. The second case is the developed method of this article where $f(x)$ and $p(y|x)$ constitute a hidden Markov model with binary states and $q(\tilde{x}|x,y)$ represents a transition matrix.

3.2 The EnKF as a special case

In the EnKF, the assumed prior model $f(x)$ is a Gaussian distribution. Combined with a Gaussian-linear likelihood $p(y|x)$, this leads to a Gaussian assumed posterior $f(x|y)$. The conditional distribution $q(\tilde{x}|x,y)$, which arises from the EnKF linear update, takes a different form depending on whether the filter is stochastic or deterministic. In stochastic EnKF, the linear update (2.5) yields a Gaussian $q(\tilde{x}|x,y)$ with mean equal to $x+\hat{K}(y-Hx)$ and covariance equal to $\hat{K}R\hat{K}^{\prime}$, that is

[TABLE]

In square root EnKF, the case is a bit different. Because the update now is deterministic, $q(\tilde{x}|x,y)$ has zero covariance and becomes a degenerate Gaussian distribution, or a delta function located at the value to which $x$ is moved, that is

[TABLE]

3.3 The proposed method as a special case

Suppose $x=(x_{1},\dots,x_{n})$ is a vector of $n$ binary variables, $x_{i}\in\{0,1\}$, and that $x$ is spatially arranged along a line. A possible assumed prior model for $x$ is then a first order Markov chain,

[TABLE]

Furthermore, suppose that for each variable $x_{i}$ there is a corresponding observation, $y_{i}$, so that $y=(y_{1},\dots,y_{n})$, and suppose that the $y_{i}$’s are conditionally independent given $x$, with $y_{i}$ depending on $x$ only through $x_{i}$,

[TABLE]

This combination of $f(x)$ and $p(y|x)$ constitutes a hidden Markov model, as introduced in Section 2. It follows that the corresponding assumed posterior $f(x|y)$ is also a first order Markov chain for which all quantities of interest are possible to compute.

Now, since $\Omega_{X}=\{0,1\}^{n}$ is a discrete sample space, we rewrite the constraint in (3.1) as a sum,

[TABLE]

As in the general case, there can be many, or infinitely many, valid solutions of $q(\tilde{x}|x,y)$. Because of the discrete context, $q(\tilde{x}|x,y)$ now represents a transition matrix. Brute force, this transition matrix involves $2^{n}(2^{n}-1)$ parameters, and the constraint (3.3) leads to a system of $2^{n}-1$ linear equations on these parameters. Even for moderate $n$, solving such a problem is too complicated. In the next section, we propose a simplifying strategy where we impose Markovian properties on $q(\tilde{x}|x,y)$, formulate an optimality criterion for $q(\tilde{x}|x,y)$, and use dynamic programming to construct the optimal solution.

4 Ensemble updating of binary state vectors

This section continues on the situation introduced in Section 3.3. In particular, we focus on the construction of $q(\tilde{x}|x,y)$. We start out in Section 4.1 proposing an optimality criterion and enforce Markovian properties on $q(\tilde{x}|x,y)$. Thereafter, in Section 4.2, we present a dynamic programming (DP) algorithm that constructs the optimal solution in a backward-forward recursive manner. Finally, in Section 4.3 we take a closer look at some more technical aspects of the backward recursion of this DP algorithm.

4.1 Choice of optimality criterion

As mentioned in the previous section, there are infinitely many valid solutions of $q(\tilde{x}|x,y)$. For us, however, it is sufficient with one solution, preferably an optimal solution, $q^{*}(\tilde{x}|x,y)$, with respect to some criterion. In this regard, an optimality criterion needs to be specified. To specify this criterion, we argue that in order for $q(\tilde{x}|x,y)$ to retain information from the prior ensemble and capture important properties of the true prior and posterior model, it should not make unnecessary changes to the prior samples. That is, as we update each prior sample $x^{(i)}$, we should take new information from the incoming observation $y$ into account and, to a certain extent, push $x^{(i)}$ towards $y$, but the adjustment we make should be minimal. With this in mind, we propose to define the optimal solution $q^{*}(\tilde{x}|x,y)$ as the one that maximises the expected number of variables, or components, of $x$ that remain unchanged after the update to $\tilde{x}$. Mathematically,

[TABLE]

where the subscript $\pi$ indicates that the expectation is taken over the following joint distribution,

[TABLE]

that is, the joint distribution of $x$ and $\tilde{x}$ given that the assumed models $f(x)$ and $f(x|y)$ are correct.

Brute force, the construction of $q(\tilde{x}|x,y)$ involves the specification of $2^{n}(2^{n}-1)=\mathcal{O}(4^{n})$ parameters. Clearly, this is intolerable even for moderate $n$. To reduce the number of parameters, we propose to enforce Markovian properties on $\pi(\tilde{x},x|y)$ as illustrated graphically in Figure 2. Then, $q(\tilde{x}|x,y)$ can be factorised as

[TABLE]

Consequently, the number of parameters reduces to $2+4(n-1)=\mathcal{O}(n)$, namely two parameters for the first factor $q(\tilde{x}_{1}|x_{1},y)$, and four parameters for each $q(\tilde{x}_{k}|\tilde{x}_{k-1},x_{k},y)$, $k=2,\dots,n$. Another, and just as important consequence of the Markovian structure, is that the optimal solution $q^{*}(\tilde{x}|x,y)$ can be efficiently computed using dynamic programming. Following (4.3), the optimal solution can also be factorised as

[TABLE]

The next section presents a DP algorithm where the $n$ factors in (4.4) are constructed recursively.

4.2 Dynamic programming

Here, we describe a DP algorithm for constructing the optimal solution $q^{*}(\tilde{x}|x,y)$ introduced in the previous section. The algorithm involves a backward recursion followed by a forward recursion. The main challenge is the backward recursion and the details therein are a bit technical. For simplicity, this section provides an overall description of the algorithm, while the more technical aspects of the backward recursion are considered separately in Section 4.3. Throughout, we use the notation $\pi(\tilde{x}_{i:j},x_{k:l}|y)$, $1\leq i\leq j\leq n$, $1\leq k\leq l\leq n$, to denote the joint distribution of $\tilde{x}_{i:j}=(\tilde{x}_{i},\dots,\tilde{x}_{j})$ and $x_{k:l}=(x_{k},\dots,x_{l})$ assuming $x$ is distributed according to $f(x)$ and $\tilde{x}$ is simulated using $q(\tilde{x}|x,y)$. Furthermore, we introduce the following simplifying notations:

[TABLE]

The backward recursion of the DP algorithm involves recursive computation of the quantities

[TABLE]

for $k=n,n-1,\dots,1$. In words, (4.5) represents the largest possible contribution of the partial expectation $\text{E}_{\pi}\left[\sum_{i=k}^{n}1(x_{i}=\tilde{x}_{i})\right]$ to the full expectation $\text{E}_{\pi}\left[\sum_{i=1}^{n}1(x_{i}=\tilde{x}_{i})\right]$ that can be obtained for a fixed $\pi(\tilde{x}_{1:k-1},x_{1:k}|y)$. The recursion uses the fact that, for $k\geq 2$, the Markovian structure of $\pi(x,\tilde{x}|y)$ yields

[TABLE]

suggesting that the full maximum value in (4.1) can be computed recursively by recursive maximisation over $q_{n}$, $q_{n-1}$, $\dots$, $q_{1}$.

An essential aspect of the backward recursion are the distributions $\pi_{1},\dots,\pi_{n}$. At each step $k$, we compute (4.5) as a function of $\pi_{k}$. Essentially, each $\pi_{k}$, $k\geq 2$, consists of four numbers, or parameters, one for each possible configuration of the pair $(\tilde{x}_{k-1},x_{k})$. However, one parameter is lost since $\pi(\tilde{x}_{k-1},x_{k}|y)$ is a probability distribution so that the four numbers necessarily sum to one. Another two parameters are lost due to the constraint (3.1) which here entails that $\pi(\tilde{x}_{k-1}|y)=f(\tilde{x}_{k-1}|y)$ and $\pi(x_{k}|y)=f(x_{k})$. Thereby, only one parameter, call it $t_{k}$, remains. This parameter $t_{k}$ is free to vary within an interval $[t_{k}^{\min},t_{k}^{\max}]$, where the bounds $t_{k}^{\min}$ and $t_{k}^{\max}$ are determined by the probabilistic nature of $\pi_{k}$. An example parametrisation is to set $t_{k}=\pi(\tilde{x}_{k-1}=0,x_{k}=0|y)$, which is the approach taken in this work. Below, the notation $\pi_{t_{k}}(\tilde{x}_{k-1},x_{k}|y)$ will, when appropriate, be used instead of $\pi(\tilde{x}_{k-1},x_{k}|y)$, in order to express the dependence on $t_{k}$ more explicitly. The chosen parameter $t_{k}$ leads to a parametrisation of $\pi_{k}$ as follows,

[TABLE]

and the bounds of the interval $[t_{k}^{\min},t_{k}^{\max}]$ are given as

[TABLE]

[TABLE]

For $k=1$, the situation is a bit different, since there is only one variable, $x_{1}$, involved in $\pi_{1}$. In fact, due to the criterion (3.1), we have $\pi(x_{1}|y)=f(x_{1})$. Consequently, $t_{1}$ is not a parameter free to vary within a certain range, but a fixed number, for example either $f(x_{1}=0)$ or $f(x_{1}=1)$. Here, we take $t_{1}=f(x_{1}=0)$.

Apart from the parametrisation of $\pi_{k}$, an essential feature of each $\pi_{k},$ for $k\geq 2,$ is its dependence on $\pi_{k-1}$ and $q_{k-1}$. This connection is due to the particular structure of $\pi(\tilde{x},x|y)$ cf. Figure 2. Generally, for $k\geq 3$, we know that $\pi_{k}$, or $\pi(\tilde{x}_{k-1},x_{k}|y)$, can be deduced from the joint distribution $\pi(\tilde{x}_{k-2},\tilde{x}_{k-1},x_{k-1},x_{k}|y)$ by summing out the variables $\tilde{x}_{k-2}$ and $x_{k-1}$,

[TABLE]

where $\pi(\tilde{x}_{k-2},\tilde{x}_{k-1},x_{k-1},x_{k}|y)$ can be written in the particular form

[TABLE]

Similarly, for the special case $k=2$, we can compute $\pi(\tilde{x}_{1},x_{2}|y)$ by summing out $x_{2}$ from $\pi(\tilde{x}_{1},x_{1},x_{2}|y)$,

[TABLE]

where $\pi(\tilde{x}_{1},x_{1},x_{2}|y)$ can be factorised as

[TABLE]

Plugging in $\tilde{x}_{k-1}=0$ and $x_{k}=0$ in (4.10), and using that $\pi_{k-1}$ is parametrised by $t_{k-1}$, we obtain a formula for $t_{k}$ in terms of $t_{k-1}$ and $q_{k-1}$, $k\geq 3$. Likewise, plugging in $\tilde{x}_{1}=0$ and $x_{2}=0$ in (4.11), and using that $f(x_{1}=0)=t_{1}$, we obtain a formula for $t_{2}$ in terms of $t_{1}$ and $q_{1}$. To express the dependence of $t_{k}$ on $t_{k-1}$ and $q_{k-1}$, $k\geq 2$, we will use the notation

[TABLE]

In some of the following equations, it will be necessary to explicitly express that (4.5) is a function of $t_{k}$. For this reason, we now define

[TABLE]

Similarly, we need a notation for the argument of the maximum in (4.6) as a function of $t_{k}$:

[TABLE]

[TABLE]

If $q^{*}_{t_{k}}(\tilde{x}_{k}|\tilde{x}_{k-1},x_{k},y)$ and $q^{*}_{t_{1}}(\tilde{x}_{1}|x_{1},y)$ are discussed in a context where the specific values of the involved variables are not important, simpler notations are preferable. In this regard, we also introduce

[TABLE]

Also, we need a notation for $\text{E}_{\pi}[1(x_{k}=\tilde{x}_{k})]$ indicating that this is a function of both $t_{k}$ and $q_{k}$,

[TABLE]

The backward recursion computes $\text{E}_{k:n}^{*}(t_{k})$ recursively for $k=n,n-1,...,1$, each step performing a maximisation over $q_{k}$ as a function of the parameter $t_{k}$. The recursion is initialised by

[TABLE]

[TABLE]

Then, for $k=n-1,n-2,\dots,1$, the recursion proceeds according to

[TABLE]

[TABLE]

Note that at the final step of the backward recursion, where $k=1$, we compute $\text{E}^{*}_{1:n}(t_{1})$ and $q^{*}_{1}(t_{1})$. Now, since we have one specific value for $t_{1}$, we also obtain one specific value for $\text{E}^{*}_{1:n}(t_{1})$ and corresponding specific values for $q^{*}_{1}(t_{1})$. This completes the backward recursion.

After the backward recursion, the forward recursion can proceed. Here, we recursively compute the specific values for $t_{2},t_{3},\dots,t_{n}$ corresponding to $\text{E}^{*}_{1:n}(t_{1})$ and $q^{*}_{1}(t_{1})$. Hence, we recursively obtain the optimal values $q^{*}(\tilde{x}_{2}|\tilde{x}_{1},x_{2},y)$, $q^{*}(\tilde{x}_{3}|\tilde{x}_{2},x_{3},y)$, $\dots$, $q^{*}(\tilde{x}_{n}|\tilde{x}_{n-1},x_{n},y)$ in (4.4). The forward recursion is initialised by

[TABLE]

[TABLE]

Then, for $k=2,3,\dots,n$, the recursion proceeds according to

[TABLE]

[TABLE]

When the forward recursion terminates, the optimal solution $q^{*}(\tilde{x}|x,y)$ is readily available.

4.3 Parametric, piecewise linear programming

In this section, we look further into the backward recursion of the DP algorithm described in the previous section. As we shall see, each step of this recursion involves the set-up of an optimisation problem that we refer to as a ’parametric, piecewise linear program’; namely, an optimisation problem with a piecewise linear objective function subject to a set of linear constraints, which at each step $k$ is solved as a function of the parameter $t_{k}$. For simplicity of writing in following equations, we introduce the notations:

[TABLE]

for $i,j\in\{0,1\}$ and $k\geq 2$.

Reconsider the initial step of the backward recursion. Here, the aim is to compute $E^{*}_{n}(t_{n})$ in (4.13) and $q_{n}^{*}(t_{n})$ in (4.14). The objective function at this step, $\text{E}_{n}(t_{n},q_{n})$, can be computed as

[TABLE]

Since $\pi^{01}_{n}(t_{n})+\pi^{11}_{n}(t_{n})=f(x_{n}=1)$, we can, after rearranging the terms, rewrite this function as

[TABLE]

As a function of the parameter $t_{n}\in[t_{n}^{\min},t_{n}^{\max}]$, we are interested in computing the solution of $q_{n}$ which maximises (4.19). In this regard, one needs to take the constraint (3.1) into account. Specifically, the constraint entails at this step that

[TABLE]

for all $\tilde{x}_{n-1},\tilde{x}_{n}\in\{0,1\}$. Hence, using that $\pi(\tilde{x}_{n-1},\tilde{x}_{n},x_{n}|y)$ $=\pi(\tilde{x}_{n-1},x_{n}|y)q(\tilde{x}_{n}|\tilde{x}_{n-1},x_{n},y)$, and that $\pi(\tilde{x}_{n-1},\tilde{x}_{n}|y)$ follows by summing out $x_{n}$ from $\pi(\tilde{x}_{n-1},\tilde{x}_{n},x_{n}|y)$, we see that $q_{n}$ must fulfil

[TABLE]

This requirement leads to four linear equations of which two are linearly independent, one where we set $\tilde{x}_{n-1}=0$ and one where we set $\tilde{x}_{n-1}=1$. Using the notations in (4.18a)-(LABEL:eq:pi_k_notation), the two linearly independent equations can be written as

[TABLE]

Additionally, we know that $q_{n}^{00},q_{n}^{01},q_{n}^{10},$ and $q_{n}^{11}$ can only take values within the interval $[0,1]$,

[TABLE]

To summarise, we want, as a function of the parameter $t_{n}\in[t_{n}^{\min},t_{n}^{\max}]$, to compute the solutions of $q_{n}^{00},q_{n}^{01},q_{n}^{10}$, and $q_{n}^{11}$ which maximise the function (4.19) subject to the constraints (4.20)-(4.22). For any fixed $t_{n}$, this is a maximisation problem where both the objective function and all the constraints are linear in $q_{n}^{00},q_{n}^{01},q_{n}^{10}$, and $q_{n}^{11}$. As such, the maximisation problem can, for a given value of $t_{n}$, be formulated as a linear program and solved accordingly. In Appendix A, we show that the optimal solutions $q_{n}^{*00}(t_{n})$, $q_{n}^{*01}(t_{n})$, $q_{n}^{*10}(t_{n})$, and $q_{n}^{*11}(t_{n})$ are piecewise-defined functions of $t_{n}$ and easy to compute analytically. Furthermore, we show that the corresponding function $E^{*}_{n}(t_{n})$, obtained by inserting $q_{n}^{*00}(t_{n})$, $q_{n}^{*01}(t_{n})$, $q_{n}^{*10}(t_{n})$, and $q_{n}^{*11}(t_{n})$ into (4.19), is a continuous piecewise linear (CPL) function of $t_{n}$.

Next, consider the intermediate steps of the backward recursion, that is $k=n-1,n-2,\dots,2.$ At each such step, the aim is to compute $E^{*}_{k:n}(t_{k})$ in (4.15) and $q_{{k}}^{*}(t_{k})$ in (4.16). The objective function at each step reads

[TABLE]

which is to be maximised with respect to $q_{k}$. Here, the first term, $E_{k}(t_{k},q_{k})$, can be computed as

[TABLE]

The second term, $E^{*}_{(k+1):n}(t_{k+1}(t_{k},q_{k}))$, is a CPL function of $t_{k+1}$. For $k=n-1$, this result is immediate, since we know from the first iteration that $E_{n}^{*}(t_{n})$ is CPL. For $k<n-1$, the result is explained in Appendix A. Since $t_{k+1}(t_{k},q_{k})$ is linear in $q_{k}$, it follows that $E^{*}_{k+1}(t_{k+1}(t_{k},q_{k}))$ is CPL in $q_{k}$ for any given $t_{k}\in[t_{k}^{\min},t_{k}^{\max}]$. Hence, the objective function (4.23) is also CPL in $q_{k}$ for any $t_{k}\in[t_{k}^{\min},t_{k}^{\max}]$. As in the first backward step, we have the following equality and inequality constraints for $q_{k}$:

[TABLE]

and

[TABLE]

Additionally, we need to incorporate constraints ensuring that $q_{k}$ and $t_{k}$ return a value $t_{k+1}$ within the interval $[t_{k+1}^{\min},t_{k+1}^{\max}]$, with $t_{k+1}^{\min}$ and $t_{k+1}^{\max}$ given by (4.8) and (4.9), respectively. That is, we require

[TABLE]

From (4.10), the formula $t_{k+1}(t_{k},q_{k})$ follows as

[TABLE]

Clearly, for any fixed $t_{k}\in[t_{k}^{\min},t_{k}^{\max}]$, all the constraints (4.25)-(4.27) are linear in $q_{k}$. However, the objective function in (4.23) is only piecewise linear. As such, we are not faced with a standard linear program, but a ’piecewise linear program’. Piecewise linear programs are a well-studied field of linear optimisation and several techniques for solving such problems have been proposed and studied, see for instance Fourer (1985, 1988, 1992). The most straightforward approach is to solve the standard linear program corresponding to each line segment of the objective function separately, and afterwards compare the solutions and store the overall optimum. This technique can be inefficient and is not recommended if the number of pieces of the objective function is relatively large. However, in our case, the objective functions normally consist of only a few pieces. For example, in the simulation experiment of Section 5.2, where a model $q(\tilde{x}|x,y)$ was constructed as much as 1,000 times, the largest number of intervals observed was 10 and the average number of intervals was 4.35. Therefore, we consider the straightforward approach as a convenient method for solving our piecewise linear programs, although we are aware that more elegant strategies exist and might have their advantages. Further details of our solution are presented below.

First, some new notations needs to be introduced. For each $2\leq k\leq n$, we let $M_{k}$ denote the number of pieces, or intervals, of $E^{*}_{k:n}(t_{k})$, and we let $t_{k}^{B(j)},$ $j=1,\dots,M_{k}+1$, denote the corresponding breakpoints. Note that for the first and last breakpoints, we have $t_{k}^{B(1)}=t_{k}^{\min}$ and $t_{k}^{B(M_{k}+1)}=t_{k}^{\max}$. Further, we let ${I}_{k}^{(j)}=[t_{k}^{B(j)},t_{k}^{B(j+1)}]\subseteq[t_{k}^{\min},t_{k}^{\max}]$ denote interval no. $j$, and $\mathcal{S}_{k}=\{1,2,\dots,M_{k}\}$ the set of interval indices. For each $j\in\mathcal{S}_{k}$, $E^{*}_{k:n}(t_{k})$ is defined by a linear function, which we denote by $E^{*(j)}_{k}(t_{k})$, whose intercept and slope we denote by $a_{k}^{(j)}$ and $b_{k}^{(j)}$, respectively.

Each linear piece, $E^{*(j)}_{k+1}(t_{k+1})$, of the piecewise linear function $E^{*}_{(k+1):n}(t_{k+1})$ leads to a standard parametric linear program. Specifically, if $E_{(k+1):n}^{*}(t_{k+1}(t_{k},q_{k}))$ in (4.23) is replaced by $E_{(k+1):n}^{*(j)}(t_{k+1}(t_{k},q_{k}))$, we obtain an objective function

[TABLE]

which is linear, not piecewise linear, as a function of $q_{k}$. The corresponding constraints for $q_{k}$ are given in (4.25) and (4.26), but instead of (4.27), we require that $t_{k}$ and $q_{k}$ return a value $t_{k+1}(t_{k},q_{k})$ within the interval $I_{k+1}^{(j)}$,

[TABLE]

Using (4.24), (4.28), and that $E^{*(j)}_{k+1}(t_{k+1})=a_{k+1}^{(j)}+b_{k+1}^{(j)}t_{k+1}$, we can for each $j\in\mathcal{S}_{k+1}$ rewrite (4.29) as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

To summarize, we obtain for each $j\in\mathcal{S}_{k+1}$ a standard parametric linear program, with the objective function given in (4.31) and the constraints given in (4.25), (4.26), and (4.30). Solving the parametric linear program corresponding to each $j\in\mathcal{S}_{k+1}$, yields the following quantities:

[TABLE]

[TABLE]

The overall maximum value $E_{k:n}^{*}(t_{k})$ and corresponding argument $q_{k}^{*}(t_{k})$ are then available as

[TABLE]

and

[TABLE]

where

[TABLE]

As previously mentioned, and as shown in Appendix A, $E_{k:n}^{*}(t_{k})$ is a CPL function of $t_{k}$. As such, $E_{k:n}^{*}(t_{k})$ is fully specified by its breakpoints and corresponding function values. The breakpoints of $E_{k:n}^{*}(t_{k})$ can be computed beforehand. Thereby, we can obtain $E_{k:n}^{*}(t_{k})$ for all values of $t_{k}$ quite efficiently, as we only need to solve our parametric, piecewise linear program at the breakpoints of $E_{k:n}^{*}(t_{k})$.

Finally, consider the last step of the backward recursion, $k=1$. Here, the aim is to compute $q^{*}_{t_{1}}(\tilde{x}_{1}|x_{1},y)$ and $E_{1:n}^{*}(t_{1})$. Essentially, this step proceeds in the same fashion as the intermediate steps, but some technicalities are a bit different since there are only two variables involved in $q_{1}$, namely $q_{1}^{0}=q(\tilde{x}_{1}=0|x_{1}=0,y)$ and $q_{1}^{1}=q(\tilde{x}_{1}=0|x_{1}=1,y)$. Also, $t_{1}$ is not a parameter free to vary within a certain range, but a fixed number, namely $t_{1}=f(x_{1}=0)$, meaning that we obtain specific values for $q^{*}_{t_{1}}(\tilde{x}_{1}|x_{1},y)$ and $E_{1:n}^{*}(t_{1})$. The function we want to maximise at this final backward step, with respect to $q_{1}$, is

[TABLE]

where now, recalling that $\pi(x_{1}|y)=f(x_{1})$, the first term, $E_{1}(t_{1},q_{1})$, can be written as

[TABLE]

Again, as in the intermediate steps, we have a piecewise linear, not a linear, objective function. To determine the constraints for $q_{1}$, we note that the requirement (3.1) for $q(\tilde{x}|x,y)$ entails that

[TABLE]

Thereby, since $t_{1}=f(x_{1}=0)$ and using that $f(\tilde{x}_{1}|y)=\sum_{x_{1}}\pi(\tilde{x}_{1},x_{1}|y)$ and $\pi(\tilde{x}_{1},x_{1}|y)=f(x_{1})q(\tilde{x}_{1}|x_{1},y)$, we see that the following requirement must be met by $q(\tilde{x}_{1}|x_{1},y)$:

[TABLE]

Additionally, we have the inequality constraints

[TABLE]

So, we are faced with a piecewise linear program, with the piecewise linear objective function (4.34) and the linear constraints (4.36) and (4.37). Again, we proceed by iterating through each linear piece of $E^{*}_{2:n}(t_{2}(t_{1},q_{1}))$, solving the standard linear program corresponding to each piece separately. That is, for each $j\in\mathcal{S}_{2}$, we replace $E^{*}_{2:n}(t_{2}(t_{1},q_{1}))$ in (4.34) by $E^{*(j)}_{2:n}(t_{2}(t_{1},q_{1}))$ and consider instead the objective function

[TABLE]

which is linear, not piecewise linear, as a function of $q_{1}$. As we did for each subproblem $j\in\mathcal{S}_{k+1}$ in every intermediate backward iteration, we must for each subproblem $j\in\mathcal{S}_{2}$ incorporate the inequality constraints

[TABLE]

where now the formula $t_{2}(t_{1},q_{1})$ follows from (4.11) and (4.12) as

[TABLE]

Using (4.35), (4.40), and that $E^{*(j)}_{2:n}(t_{2})=a_{2}^{(j)}+b_{2}^{(j)}t_{2}$, we can rewrite the function in (4.38) as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

To summarize, we obtain for each $j\in\mathcal{S}_{2}$ a standard linear program, where the aim is to maximise the objective function (4.41) with respect to $q_{1}$ subject to the constraints (4.36), (4.37) and (4.39). This program is solved for $t_{1}=f(x_{1}=0)$. Analogously to (4.32) and (4.33), let

[TABLE]

[TABLE]

Ultimately, we obtain

[TABLE]

and

[TABLE]

where

[TABLE]

5 Numerical experiments

In this section, we present two empirical studies with simulated data. In Section 5.1, we present a toy example where the assumed prior $f(x)$ is a given stationary Markov chain of length $n=4$. Here, we focus on the construction of $q(\tilde{x}|x,y)$ for this assumed prior model, not on the application of it in an ensemble-based context. In Section 5.2, we consider a higher-dimensional and ensemble-based example, inspired by the movement, or flow, of water and oil in a petroleum reservoir.

5.1 Toy example

Suppose the assumed prior $f(x)$ is a Markov chain of length $n=4$ with homogenous transition probabilities $f(x_{k}=0|x_{k-1}=0)=0.7$ and $f(x_{k}=1|x_{k-1}=1)=0.8$ for $k\geq 2$, and initial distribution $f(x_{1})$ equal to the associated limiting distribution. Thus, the Markov chain $f(x)$ is stationary with marginal probabilities $f(x_{k}=0)=0.40$ for each $k=1,2,3,4$. Further, suppose every factor $f(y_{i}|x_{i})$ of the likelihood model $f(y|x)$ is a Gaussian distribution with mean $x_{i}$ and standard deviation $\sigma=2$, and consider the observation vector $y=(-0.681,-1.585,0.007,3.103)$. The corresponding posterior Markov chain model $f(x|y)$ then have the transition probabilities

[TABLE]

and marginal distributions

[TABLE]

Based on the prior model $f(x)$ and the posterior model $f(x|y)$, we can construct $q^{*}(\tilde{x}|x,y)$ as discussed above. For this simple example, this involves computing fourteen quantities, namely $q_{1}^{*0}(t_{1}^{*})=q^{*}(\tilde{x}_{1}=0|x_{1}=0,y)$, $q_{1}^{*1}(t_{1}^{*})=q^{*}(\tilde{x}_{1}=0|x_{1}=1,y)$, $q_{k}^{*ij}(t_{k}^{*})=q^{*}(\tilde{x}_{k}=0|\tilde{x}_{k-1}=i,x_{k}=j,y)$, for $k=2,3,4$, and $i,j=0,1$. As described in the previous section, the construction of $q(\tilde{x}|x,y)$ involves a backward recursion and a forward recursion. In the backward recursion, we compute $E_{k:n}^{*}(t_{k})$ and $q_{k}^{*00}(t_{k})$, for $k=4,3,2$. The results for these quantities are presented in Figure 3.

In the forward recursion, we start out computing the optimal solution of the first factor, $q^{*}(\tilde{x}_{1}|x_{1},y)$, and then compute the remaining optimal parameter values $t_{2}^{*}$, $t_{3}^{*}$ and $t_{4}^{*}$ and corresponding optimal solutions $q_{k}^{*ij}(t_{k}^{*})$, $k=2,3,4$, $i,j=0,1$. The results from the forward recursion are given in Table 1.

Taking a closer look at the results for the optimal solution $q^{*}(\tilde{x}|x,y)$, we see that many of the probabilities $q_{k}^{*ij}(t_{k}^{*})$ are either zero or one. This feature can be formally explained mathematically (see Appendix A), but is also quite an intuitive result, and has to do with how the probabilities of the prior model $f(x)$ differ from the corresponding probabilities of the posterior model $f(x|y)$. Often, if $f(x_{k}=0)<f(x_{k}=0|y)$, we obtain $q_{k}^{*00}(t_{k}^{*})=1$ and $q_{k}^{*10}(t_{k}^{*})=1$, while $q_{k}^{*01}(t_{k}^{*})$ and $q_{k}^{*11}(t_{k}^{*})$ take values somewhere between zero and one. Thus, if we have a prior sample $x$ with $x_{k}=0$, the update of $x$ to $\tilde{x}$ is always such that $\tilde{x}_{k}=0$. Specifically, in our toy example, this is the case for $k=2$; that is, we have $f(x_{2}=0)<f(x_{2}=0|y)$, and obtained $q_{2}^{*00}(t_{2}^{*})=1$ and $q_{2}^{*10}(t_{2}^{*})=1$. Likewise, if $f(x_{k}=0)>f(x_{k}=0|y)$, we often obtain $q_{k}^{*01}(t_{k}^{*})=0$ and $q_{k}^{*11}(t_{k}^{*})=0$, while $q_{k}^{*00}(t_{k}^{*})$ and $q_{k}^{*10}(t_{k}^{*})$ take values somewhere between zero and one. Thus, if we have a prior sample $x$ with $x_{k}=1$, the update of $x$ to $\tilde{x}$ is always such that $\tilde{x}_{k}=1$. In our toy example, this is the case for $k=4$; that is, we have $f(x_{4}=0)>f(x_{4}=0|y)$, and obtained $q_{4}^{*01}(t_{4}^{*})=0$ and $q_{4}^{*10}(t_{4}^{*})=0$. However, the model $q(\tilde{x}|x,y)$ is not only constructed so that the marginal probabilities in (5.2) are fulfilled, but also so that the posterior transition probabilities in (5.1) are reproduced. In our toy example, we see for example that for $k=3$ we obtained $q_{3}^{*10}(t_{3}^{*})<1$ even if $f(x_{3}=0)<f(x_{3}=0|y)$. Instead, we observe another deterministic term, namely $q_{3}^{*11}(t_{3}^{*})=0$.

5.2 Ensemble-based, higher-dimensional example with simulated data

Until now, we have focused on the ensemble updating problem at a specific time step of the filtering recursions. However, in a practical application, one is interested in the filtering problem as a whole and needs to cope with the ensemble updating problem sequentially for $t=1,2,\dots,T$. We now address this issue and investigate the application of our proposed approach in this context. That is, we reconsider the situation with an unobserved Markov process, $\{x^{t}\}_{t=1}^{T}$, and a corresponding time series of observations, $\{y^{t}\}_{t=1}^{T}$, and at every time step $t=1,\dots,T$, we construct a distribution $q(\tilde{x}^{t}|x^{t},y^{1:t})$ in order to update the prior ensemble $\mathcal{X}^{t}=\{x^{t(1)},x^{t(2)},\dots,x^{t(M)}\}$ to a posterior ensemble $\widetilde{\mathcal{X}}^{t}=\{\tilde{x}^{t(1)},\tilde{x}^{t(2)},\dots,\tilde{x}^{t(M)}\}$. Below, we first present the set-up of our simulation example in Section 5.2.1, and thereafter study the performance of our approach in Sections 5.2.2 and 5.2.3.

5.2.1 Specification of simulation example

To construct a simulation example we must first define the $\{x^{t}\}_{t=1}^{T}$ Markov chain. We set $T=100$ and let $x^{t}=(x_{1}^{t},\ldots,x_{n}^{t})$ be an $n=400$ dimensional vector of binary variables $x_{i}^{t}\in\{0,1\}$ for each $t=1,\ldots,T$. To simplify the specification of the transition probabilities $p(x^{t}|x^{t-1})$ we make two Markov assumptions. First, conditioned on $x^{t-1}$ we assume the elements in $x^{t}$ to be a Markov chain so that

[TABLE]

The second Markov assumption we make is that

[TABLE]

for $i=2,\ldots,n-1$, i.e. the value in element $i$ at time $t$ only depends on the values in elements $i-1$, $i$ and $i+1$ at the previous time step. For $i=1$ and $i=n$ we make the corresponding Markov assumptions

[TABLE]

To specify the $x^{t}$ Markov process we thereby need to specify $p(x_{i}^{t}|x_{i-1}^{t},x_{i-1}^{t-1},x_{i}^{t-1},x_{i+1}^{t-1})$ for $t=2,\ldots,T$ and $i=2,\ldots,n$ and the corresponding probabilities for $t=1$ and for $i=1$ and $n$.

To get a reasonable test for how our proposed ensemble updating procedure works we want an $\{x^{t}\}_{t=1}^{T}$ process with a quite strong dependence between $x^{t-1}$ and $x^{t}$, also when conditioning on observed data. Moreover, conditioned on $y^{1:t}$, the elements in $x^{t}$ should not be first order Markov so that the true model differ from the assumed Markov model defined in Sections 3.1 and 3.3. In the following we first discuss the choice of $p(x_{i}^{t}|x_{i-1}^{t},x_{i-1}^{t-1},x_{i}^{t-1},x_{i+1}^{t-1})$ for $t=2,\ldots,T$ and $i=2,\ldots,n$ and thereafter specify how these are modified for $t=1$ and for $i=1$ and $n$. When specifying the probabilities we are inspired by the process of how water comes through to an oil producing well in a petroleum reservoir, but without claiming our model to be a very realistic model for this situation. Thereby $t$ represents time and $i$ is location in the well. We let $x_{i}^{t}=0$ represent the presence of oil at location or node $i$ at time $t$ and correspondingly $x_{i}^{t}=1$ represents the presence of water. In the start we assume oil is present in the whole well, but as time goes by more and more water is present and at time $t=T$ water has become the dominating fluid in the well. Whenever $x_{i}^{t-1}=1$ we therefore want $x_{i}^{t}=1$ with very high probability, especially if also $x_{i-1}^{t}=1$. If $x_{i}^{t-1}=0$ we correspondingly want a high probability for $x_{i}^{t}=0$ unless $x_{i-1}^{t}=1$ and $x_{i-1}^{t-1}=x_{i+1}^{t-1}=1$. Trying different sets of parameter values according to these rules we found that the values specified in Table 2 gave realisations consistent with the requirements discussed above.

One realisation from this model is shown in Figure 4(a),

where black and white represent [math] (oil) and $1$ (water), respectively. The corresponding probabilities when $t=1$ and for $i=1$ and $n$ we simply defined from the values in Table 2 by defining all values lying outside the $\{(i,t):i=1,\ldots,n;t=1,\ldots,T\}$ lattice to be zero. In particular this implies that at time $t=0$, which is outside the lattice, oil is present in the whole well. In the following we consider the realisation shown in Figure 4(a) to be the (unknown) true $x^{t}$ process.

The next step in specifying the simulation example is to specify an observational process. For this we simply assume one scalar observation $y_{i}^{t}$ for each node $i$ at each time $t$, and assume the elements in $y^{t}=(y_{1}^{t},\ldots,y_{n}^{t})$ to be conditionally independent given $x^{t}$. Furthermore, we let $y_{i}^{t}$ be Gaussian with mean $x_{i}^{t}$ and variance $\sigma^{2}$. As we want the dependence between $x^{t-1}$ and $x^{t}$ to be quite strong also when conditioning on the observations, we need to choose the variance $\sigma^{2}$ reasonably large, so we set $\sigma^{2}=2^{2}$. Given the true $x^{t}$ process shown in Figure 4 we simulate $y_{i}^{t}$ values from the specified Gaussian distribution, and in the following consider these values as observations. An image of these observations is not included, since the variance $\sigma^{2}$ is so high that such an image is not very informative.

Pretending that the $\{x^{t}\}_{t=1}^{T}$ process is unknown and that we only have the observations $\{y^{t}\}_{t=1}^{T}$ available, our aim with this simulation study is to evaluate how well our proposed ensemble based filtering procedure is able to capture the properties of the correct filtering distributions $p(x^{t}|y^{1:t}),t=1,\ldots,T$. To do so we first need to evaluate the properties of the correct filtering distribution. It is possible to get samples from $p(x^{t}|y^{1:t})$ by simulating from $p(x^{1:t}|y^{1:t})$ with a Metropolis–Hastings algorithm, but to a very high computational cost as a separate Metropolis–Hastings run must be performed for each value of $t$. Nevertheless, we do this to get the optimal solution of the filtering problems to which we can compare the results of our proposed ensemble based filtering procedure. In our algorithm for simulating from $p(x^{1:t}|y^{1:t})$ we combine single site Gibbs updates of each element in $x^{1:t}$ with a one-block Metropolis–Hastings update of all elements in $x^{1:t}$. To get a reasonable acceptance rate for the one-block proposals we adopt the approximation procedure introduced in Austad and Tjelmeland (2017) to obtain a partially ordered Markov model (Cressie and Davidson, 1998) approximation to $p(x^{1:t}|y^{1:t})$, propose potential new values for $x^{1:t}$ from this approximate posterior, and accept or reject the proposed values according to the usual Metropolis–Hastings acceptance probability. For each value of $t$ we run the Metropolis–Hastings algorithm for a large number of iterations and discard a burn-in period. From the generated realisations we can then estimate the properties of $p(x^{t}|y^{1:t})$. In particular we can estimate the marginal probabilities $p(x_{i}^{t}=1|y^{1:t})$ as the fraction of realisations in which the simulated $x_{i}^{t}=1$. We denote these estimates of the correct filtering probabilities by $\widehat{p}_{c}(x_{i}^{t}=1|y^{1:t})$. In Figure 4(b) all these estimates are visualised as a grayscale image, where black and white correspond to $\widehat{p}_{c}(x_{i}^{t}=1|y^{1:t})$ equal to zero and one, respectively. It is important to note that Figure 4(b) is not showing the solution of the smoothing problem, but the solution of many filtering problems put together as one image.

Now, we are ready to run our proposed ensemble-based filtering method. The input to our algorithm include the simulated observations $\{y^{t}\}_{t=1}^{T}$, the observational model $p(y^{t}|x^{t})$, the initial model $p(x^{1})$ for the $\{x^{t}\}_{t=1}^{T}$ process, and an ensemble size $M$. Since in real-world situations it is typically necessary to choose $M$ much smaller than $n$ for computational reasons, we chose to set $M=20$. A problem, however, when the ensemble size is this small, is that results may vary a lot from run to run. To quantify this between-run variability, we therefore reran our proposed approach a total of $B=1,000$ times, each time with a new initial ensemble of $M=20$ realisations from the initial model $p(x^{1})$. At each time step $t$ we thus achieved a total of $MB=20,000$ posterior samples of the state vector $x^{t}$, which can be used to construct an estimate, denoted $\hat{p}_{q}(x^{t}|y^{1:t})$, for the true filtering distribution $p(x^{t}|y^{1:t})$.

An important step of our approach is the estimation of a first order Markov chain $f(x^{t}|y^{1:t-1})$ at each time step $t$. Basically, this involves estimating an initial distribution $f(x^{t}_{1}|y^{1:t-1})$ and $n-1$ transition matrices $f(x^{t}_{i+1}|x^{t}_{i},y^{1:t-1})$, $i=1,\dots,n-1$. Since each component $x^{t}_{i}$ is a binary variable, the initial distribution $f(x^{t}_{1}|y^{1:t-1})$ can be represented by one parameter, while the transition matrices $f(x^{t}_{i+1}|x^{t}_{i},y^{1:t-1})$ each require two parameters. In this example, we pursued a Bayesian approach for estimating the parameters. Specifically, if we let $\theta^{t}$ represent a vector containing all the parameters required to specify the model $f(x^{t}|y^{1:t-1})$, we put a prior on $\theta^{t}$, $f(\theta^{t})$, and then set the final estimator for $\theta^{t}$ equal to the mean of the corresponding posterior distribution $f(\theta^{t}|\mathcal{X}^{t})$. In the specification of $f(\theta^{t})$ we assumed that all the parameters in the vector $\theta^{t}$ are independent and that each parameter follows a Beta distribution $\mathcal{B}(\alpha,\beta)$ with parameters $\alpha=2,\beta=2$.

As well as studying the performance of our proposed approach, we are in this experiment interested in studying the results one would get without constructing $q(\tilde{x}^{t}|x^{t},y^{1:t})$ at each time step $t$, instead updating the prior ensemble by simulating independent samples from the assumed Markov chain model $f(x^{t}|y^{1:t})$. Below, we refer to this procedure as the assumed model approach. As with our approach, we reran the assumed model approach $B=1,000$ times, yielding a total of $MB=20,000$ posterior samples of each state vector $x^{t}$, $t=1,\dots,T$. These samples can then be used to construct an estimate, denoted $\hat{p}_{a}(x^{t}|y^{1:t})$, for the true filtering distribution $p(x^{t}|y^{1:t})$. By comparing $\hat{p}_{a}(x^{t}|y^{1:t})$ and $\hat{p}_{q}(x^{t}|y^{1:t})$ with $\hat{p}_{c}(x^{t}|y^{1:t})$, which essentially represents the true model $p(x^{t}|y^{1:t})$, we can get an understanding of how much we gain by doing our approach instead of the much simpler assumed model approach. In the next two sections we investigate how well $\hat{p}_{q}(x^{t}|y^{1:t})$ and $\hat{p}_{a}(x^{t}|y^{1:t})$ capture marginal and joint properties of the true distribution $p(x^{t}|y^{1:t})$, for which the MCMC estimate $\hat{p}_{c}(x^{t}|y^{1:t})$ works as a reference.

5.2.2 Evaluation of marginal distributions

In this section, we are interested in studying how well our approach estimates the marginal filtering distributions $p(x^{t}_{i}|y^{1:t})$. Following the notations introduced above, we let $\hat{p}_{q}(x^{t}_{i}|y^{1:t})$ and $\hat{p}_{a}(x^{t}_{i}|y^{1:t})$ denote estimates of the marginal distribution $p(x^{t}_{i}|y^{1:t})$ obtained with our approach and the assumed model approach, respectively. The values of $\hat{p}_{q}(x^{t}_{i}=1|y^{1:t})$ and $\hat{p}_{a}(x^{t}_{i}=1|y^{1:t})$ are in each case set equal to the mean of the corresponding set of samples of $x_{i}^{t}$. Figures 4(c) and (d) present grayscale images of $\left\{\hat{p}_{q}(x^{t}_{i}=1|y^{1:t})\right\}_{t=1}^{100}$ and $\left\{\hat{p}_{a}(x^{t}_{i}=1|y^{1:t})\right\}_{t=1}^{100}$, respectively. From a visual inspection, the image of $\left\{\hat{p}_{q}(x^{t}_{i}=1|y^{1:t})\right\}_{t=1}^{100}$ is more grey and noisy than that of $\left\{\hat{p}_{c}(x^{t}_{i}=1|y^{1:t})\right\}_{t=1}^{100}$, which contains more tones closer to pure black and white. This is to be expected, since $\left\{\hat{p}_{c}(x^{t}_{i}|y^{1:t})\right\}_{t=1}^{100}$ essentially is the ideal solution, and we do not expect our approach to perform this well. However, the image of $\left\{\hat{p}_{a}(x^{t}_{i}=1|y^{1:t})\right\}_{t=1}^{100}$ is even more grey and noisy than $\left\{\hat{p}_{q}(x^{t}_{i}=1|y^{1:t})\right\}_{t=1}^{100}$, so presumably we do gain something by running our approach instead of the simpler assumed model approach. To investigate this further, we compute the Frobenius norms of the two matrices produced by subtracting the true marginal probabilities $\hat{p}_{c}(x^{t}_{i}=1|y^{1:t})$ from the corresponding estimates $\hat{p}_{q}(x^{t}_{i}=1|y^{1:t})$ and $\hat{p}_{a}(x^{t}_{i}=1|y^{1:t})$. We then obtain the numbers $35.38$ and $63.00$, respectively. That is, the Frobenius norm of the difference between the true and the estimated marginal filtering distributions is reduced to almost the half with our method compared to the assumed model approach, clearly indicating that our method overall provides much better estimates of the marginal distributions $p(x_{i}^{t}|y^{1:t})$.

To look further into the accuracy of the estimates $\hat{p}_{q}(x^{t}_{i}=1|y^{1:t})$ and $\hat{p}_{a}(x^{t}_{i}=1|y^{1:t})$ and to study their variability from run to run, we take a closer look at the results for some specific time steps. For each of these time steps we compute a 90$\%$ quantile interval for each of the estimates $\hat{p}_{q}(x_{i}^{t}=1|y^{1:t})$ and $\hat{p}_{a}(x_{i}^{t}=1|y^{1:t})$, $i=1,\dots,400$. To compute the quantile intervals, recall that our proposed approach and the assumed model approach were both rerun $B=1,000$ times. From run $b=1,\dots,B$ of our approach, we obtained an estimate $\hat{p}_{q}^{(b)}(x_{i}^{t}|y^{1:t})$ of $p(x_{i}^{t}|y^{1:t})$ for each $i$. Likewise, from run $b=1,\dots,B$ of the assumed model approach, we obtained an estimate $\hat{p}_{a}^{(b)}(x_{i}^{t}|y^{1:t})$ of $p(x_{i}^{t}|y^{1:t})$ for each $i$. Hence, for each marginal distribution $p(x_{i}^{t}|y^{1:t})$, we have $B=1,000$ estimates $\{\hat{p}_{q}^{(b)}(x_{i}^{t}|y^{1:t})\}_{b=1}^{B}$ obtained with our approach and $B=1,000$ estimates $\{\hat{p}_{a}^{(b)}(x_{i}^{t}|y^{1:t})\}_{b=1}^{B}$ obtained with the assumed model approach. From these two sets of samples, corresponding quantile intervals for $\hat{p}_{q}(x_{i}^{t}=1|y^{1:t})$ and $\hat{p}_{a}(x_{i}^{t}=1|y^{1:t})$ can be constructed. Figure 5

presents the computed results for time step $t=60$. For simplicity, we have not included similar figures from any of the other time steps that we studied, since they look very much the same as those obtained at time $t=60$. According to Figures 5(a) and (b), it seems that the essentially true value $\hat{p}_{c}(x^{60}_{i}|y^{1:60})$ typically is within the $90\%$ quantile interval corresponding to $\hat{p}_{q}(x_{i}^{60}|y^{1:60})$, but often closer to one of the interval boundaries rather than the estimate $\hat{p}_{q}(x_{i}^{60}|y^{1:60})$ itself. In particular, we note that $\hat{p}_{c}(x^{60}_{i}|y^{1:60})$ often is close to either zero or one, while $\hat{p}_{q}(x_{i}^{60}|y^{1:60})$ is a bit higher than zero or a bit lower than one. This is not unreasonable, since we have used approximations to construct $\hat{p}_{q}(x_{i}^{60}|y^{1:60})$. Thereby, we loose information about the true quantity $\hat{p}_{c}(x^{60}_{i}|y^{1:60})$ and end up with estimated values closer to 0.5. From Figures 5(c) and (d), we observe that this is even more the case for the estimate $\hat{p}_{a}(x_{i}^{60}|y^{1:60})$ whose quantile interval often not even covers $\hat{p}_{c}(x^{60}_{i}|y^{1:60})$.

5.2.3 Evaluation of joint distributions

In this section, we want to evaluate how well our proposed approach manages to capture properties about the joint distribution $p(x^{t}|y^{1:t})$. To do so, we select three specific time steps to study, namely $t=60$, $t=70$, and $t=80$. For each of these steps, we perform two tests on our samples, both concerning a feature we call contact between a pair of nodes of $x^{t}$. So first, we need to explain the concept of contact between a pair of nodes. Consider two components $x_{i}^{t}$ and $x_{j}^{t}$ at a given time step $t$. Given that $x_{i}^{t}$ is equal to one, we say that there is contact between node $i$ and node $j$ in $x^{t}$ if all components of $x^{t}$ between and including node $i$ and node $j$ are equal to one. That is, there is contact between node $i$ and $j$, given that $x_{i}^{t}$ is equal to one, if the function

[TABLE]

is equal to one.

Keeping $i$ fixed, we are in our first test interested in studying the probability that there is contact between node $i$ and node $j$ for various values of $j$, given that $x_{i}^{t}$ is equal to one. Mathematically, that means we are interested in

[TABLE]

It is most informative to study (5.3) for a node $i$ whose corresponding component $x_{i}^{t}$ has a high probability of being equal to one. Therefore, we concentrate on estimating (5.3) for three specific choices of $i$, each corresponding to a component $x_{i}^{t}$ with a relatively high probability of being equal to one. According to Figure 4 this appears to be the case for the three nodes $i=115$, $i=210$ and $i=290$ at all three time steps $t=60$, $t=70$ and $t=80$. For each $i$ and $t$, we can then use our three sets of samples of $x^{t}$ to obtain three different estimates of (5.3) for all $j$. Following previous notations, we let $\hat{p}_{c}^{t}(i,j)$ denote the MCMC estimate of $p^{t}(i,j)$, while $\hat{p}_{q}^{t}(i,j)$ and $\hat{p}_{a}^{t}(i,j)$ denote the estimates obtained with our approach and the assumed model approach, respectively. Figure 6 presents the computed results. Comparing the curves representing the estimates $\hat{p}_{c}^{t}(i,j)$, $\hat{p}_{q}^{t}(i,j)$ and $\hat{p}_{a}^{t}(i,j)$, we observe that $\hat{p}_{q}^{t}(i,j)$ and $\hat{p}_{a}^{t}(i,j)$ typically decrease to zero for increasing values of $j$ quicker than $\hat{p}_{c}^{t}(i,j)$. However, we see that $\hat{p}_{a}^{t}(i,j)$ decreases considerably faster than $\hat{p}_{q}^{t}(i,j)$. This makes sense, since the posterior samples used to construct the estimate $\hat{p}_{a}^{t}(i,j)$ are drawn independently from the assumed model $f(x^{t}|y^{1:t})$, not taking the state of the prior samples into account.

In our second test, we focus on the total number of nodes an arbitrary node $i$ with $x_{i}^{t}=1$ is in contact with. We denote this quantity by $L_{i}(x^{t})$. Mathematically, $L_{i}(x^{t})$ can be written

[TABLE]

For each time step $t=60,70$ and 80, we want to study the cumulative distribution of $L_{i}(x^{t})$,

[TABLE]

when randomising over both $i$ and $x^{t}$, with $i\sim\text{unif}\{1,n\}$ and $x^{t}\sim p(x^{t}|y^{1:t})$. Again, we can use our three sets of samples to construct three different estimates of (5.4). That is, we can construct $\hat{F}_{c}(\cdot)$ from the MCMC samples, $\hat{F}_{q}(\cdot)$ from the samples generated with our approach, and $\hat{F}_{a}(\cdot)$ from the samples generated with the assumed model approach. Figure 7 presents the results.

Here, we can see that $\hat{F}_{a}(\cdot)$ is above $\hat{F}_{c}(\cdot)$ at all three time steps $t=60$, 70 and 80, indicating that $L_{i}(x^{t})$ typically is too small and that the assumed model approach seems to underestimate the level of contact between nodes. This makes sense and agrees with the performance of $\hat{p}^{t}_{a}(i,j)$ studied in our first test. According to Figures 7(b) and (c), our estimate $\hat{F}_{q}(\cdot)$ appears to do a better job, since it is relatively close to $\hat{F}_{c}(\cdot)$. However, this is not the case in Figure 7(a). Here, the curve for $\hat{F}_{q}(\cdot)$ is below $\hat{F}_{c}(\cdot)$, suggesting that $L_{i}(x^{t})$ typically is too high. To investigate this further, we also examined corresponding output for other time steps $t$. Then, we observed that for smaller values of $t$, typically smaller than $60$, the curve for $\hat{F}_{q}(\cdot)$ tends to be below $\hat{F}_{c}(\cdot)$, while for larger values of $t$, it tends to be quite close to $\hat{F}_{c}(\cdot)$. This is in fact not so unreasonable, since it is for higher values of $t$ that the value one (i.e. water) is dominant in $x^{t}$. For smaller values of $t$, the value zero (i.e. oil) becomes more and more dominant, and the length of one-valued chains is not supposed to be very high. It appears as if our requirement of ’maximising the expected number of unchanged components’ results in keeping too much information from the prior samples.

6 Closing remarks

In this article, we present an approximate and ensemble-based method for solving the filtering problem. The method is particularly designed for binary state vectors and is based on a generalised view of the well-known ensemble Kalman filter. In the EnKF, a Gaussian approximation $f(x)$ for the true prior is constructed which combined with a linear-Gaussian likelihood yields a Gaussian approximation $f(x|y)$ to the true posterior. Next, the prior ensemble is updated with a linear shift such that the marginal distribution of the updated samples is equal to $f(x|y)$ provided that the distribution of the prior samples is equal to $f(x)$. In our approach, we instead pursue a first order Markov chain for $f(x)$ and combine this with a particular likelihood model such that the corresponding posterior model $f(x|y)$ is also a first order Markov chain. To update the prior samples, we instead of a linear shift construct a distribution $q(\tilde{x}|x,y)$ and simulate the updated samples from this distribution. In the construction of $q(\tilde{x}|x,y)$, we formulate an optimality criterion and, just as in EnKF, require that the distribution of the updated samples is equal to $f(x|y)$ provided that the distribution of the prior samples is equal to $f(x)$. To compute the optimal solution of $q(\tilde{x}|x,y)$ we combine dynamic and linear programming. Based on results from a simulation experiment, the performance of our method seems promising.

The focus of this article is on binary state vectors with a one-dimensional spatial arrangement. Clearly, this is a very simple situation with limited practical interest, considering that most real problems involve at least two spatial dimensions and multiple classes for the state variables. Nevertheless, we consider the work of this article as a first step towards a more applicable method, and in the future we would like to explore possible extensions of our method. Conceptually, most of the material presented in the article can easily be generalised to more complicated situations. Computationally, however, it is more challenging. A generalisation of the material in Sections 3 and 4 to a similar situation with more than two possible classes, involves a growing number of free parameters in the construction of each factor $q(\tilde{x}_{k}|\tilde{x}_{k-1},x_{k},y)$. Specifically, in the case of three classes there will be four parameters involved, while in the case of four classes there will be nine parameters involved. We believe, however, that it is possible to cope with a situation with more than one free parameter via an iterative procedure. Specifically, one can start with some initial values for each of the free parameters and thereafter iteratively optimise with respect to one of the parameters at a time, keeping the other parameters fixed. By iterating until convergence we thereby obtain the optimal solution. How many parameters we are able to deal with using this strategy will depend on how fast convergence is reached and, of course, how much computation time one is willing to use.

Another possible extension of our method is to pursue a higher order Markov chain for the assumed prior model $f(x)$. If this is within reach, a further generalisation to two spatial dimensions may be possible by choosing a Markov mesh model (Abend et al., 1965) for $f(x)$. However, similarly to the case of multiple classes, the computational complexity grows rapidly with the order of the Markov chain. The higher the order, the higher the number of free parameters there will be in the construction of each factor $q(\tilde{x}_{k}|\tilde{x}_{k-1},x_{k},y)$. Computationally we can again imagine to cope with this situation by adopting an iterative optimisation algorithm as discussed above.

An optimality criterion needs to be specified when constructing $q(\tilde{x}|x,y)$. In our work we choose to define the optimal solution as the one that ’maximises the expected number of equal components’. To us this seems like an intuitively reasonable criterion, since we want to retain as much information as possible from the prior samples. However, there may be other criteria that are more suitable and which might improve the performance of our procedure. What optimality criterion that give the best results may even depend on how the true and assumed distributions differ. One may therefore imagine to construct a procedure which at each time $t$ use the prior samples to estimate, or select, the best optimality criterion within a specified class.

Appendix A Appendix

This appendix provides an informal proof of that $E_{k:n}^{*}(t_{k})$, $2\leq k\leq n$, is continuous piecewise linear (CPL). Every iteration of the backward recursion, except the first, relies on this result. The proof is an induction proof and consists of two main steps. First, in Section A.1, we consider the first step of the backward recursion and prove that $E_{n}^{*}(t_{n})$ is CPL. This corresponds to the ’base case’ of our induction proof. Next, in Section A.2, we consider the intermediate steps and prove that $E_{k:n}^{*}(t_{k})$ is also CPL, given that $E_{k+1:n}^{*}(t_{k+1})$ is CPL, $2\leq k<n$. This corresponds to the ’inductive step’ of our induction proof. In Section A.3 of the appendix, we explain how to determine the breakpoints of $E_{k:n}^{*}(t_{k})$, $2\leq k<n$, prior to solving the corresponding parametric, piecewise linear program. This is crucial in order to avoid a numerical computation of $E_{k:n}^{*}(t_{k})$ on a grid of $t_{k}$-values. Throughout the appendix, we assume the reader is familiar with all notations introduced previously in the main parts of the article.

A.1 The first iteration

The parametric linear program of the first backward iteration can easily be computed analytically. Because of the equality constraints, (4.20) and (4.21), we can reformulate the optimisation problem in terms of two variables instead of four. More specifically, we can choose either $q_{n}^{00}$ or $q_{n}^{01}$ from (4.20), together with either $q_{n}^{10}$ or $q_{n}^{11}$ from (4.21), and then reformulate the problem in terms of the two chosen variables. Here, we pursue $q_{n}^{00}$ and $q_{n}^{10}$. By rearranging terms in (4.20) and (4.21), we can write

[TABLE]

Now, if we replace the terms $\pi_{n}^{01}(t_{n})q_{n}^{01}$ and $\pi_{n}^{11}(t_{n})q_{n}^{11}$ in the objective function $E_{n}(t_{n},q_{n})$ cf. (4.19) with the right hand side expressions in (A.1) and (A.2), respectively, we can rewrite $E_{n}(t_{n},q_{n})$ in terms of $q_{n}^{00}$ and $q_{n}^{10}$ as

[TABLE]

where $c_{n}$ is a constant given as

[TABLE]

Furthermore, combining (A.1) and (A.2) with the inequality constraints (4.22) allows us to reformulate the constraints for $q_{n}^{00}$ and $q_{n}^{10}$ as

[TABLE]

[TABLE]

To summarize, we have now obtained a linear program, in which we want to maximise the objective function in (A.3) with respect to the two variables $q_{n}^{00}$ and $q_{n}^{10}$, subject to the constraints (A.4)-(A.5).

If for some fixed $t_{n}\in[t_{n}^{\min},t_{n}^{\max}]$ we consider a coordinate system with $q_{n}^{00}$ along the first axis and $q_{n}^{10}$ along the second, the constraints (A.4)-(A.5) form a rectangular region of feasible solutions, with two edges in the $q_{n}^{00}$-direction and two edges in the $q_{n}^{10}$-direction. The optimal solution lies in a corner point of this region. Since $\pi_{n}^{00}(t_{n})$ and $\pi_{n}^{10}(t_{n})$ are non-negative for any $t_{n}\in[t_{n}^{\min},t_{n}^{\max}]$, it is easily seen from (A.3) that $E_{n}(t_{n},q_{n})$ is maximised with respect to $q_{n}$ when $q_{n}^{00}$ and $q_{n}^{10}$ are as large as possible. Consequently, the optimal solutions of $q_{n}^{00}$ and $q_{n}^{10}$ must equal the upper bounds in (A.4) and (A.5), corresponding to the upper right corner of the rectangular feasible region. That is,

[TABLE]

Clearly, $q^{*00}_{n}(t_{n})$ and $q^{*10}_{n}(t_{n})$ are continuous and piecewise-defined functions of $t_{n}$, since $\pi_{n}^{00}(t_{n})$ and $\pi_{n}^{10}(t_{n})$ are linear functions of $t_{n}$. Specifically, for $t_{n}$-values such that $\pi_{n}^{00}(t_{n})>f_{n}^{00}$, we get $q^{*00}_{n}(t_{n})=f_{n}^{00}/\pi_{n}^{00}(t_{n})$, while for $t_{n}$-values such that $\pi_{n}^{00}(t_{n})\leq f_{n}^{00}$, we get $q^{*00}_{n}(t_{n})=1$. Likewise, for $t_{n}$-values such that $\pi_{n}^{10}(t_{n})>f_{n}^{10}$, we get $q^{*10}_{n}(t_{n})=f_{n}^{10}/\pi_{n}^{10}(t_{n})$, while for $t_{n}$-values such that $\pi_{n}^{10}(t_{n})\leq f_{n}^{10}$, we get $q^{*10}_{n}(t_{n})=1$.

Inserting the optimal solutions $q^{*00}_{n}(t_{n})$ and $q^{*10}_{n}(t_{n})$ into (A.3), returns $E^{*}_{n}(t_{n})$. Doing this, it is easily seen that $E^{*}_{n}(t_{n})$ is a CPL function of $t_{n}$, consisting of maximally three pieces, each piece having a slope equal to either -2, 0 or 2.

A.2 The intermediate iterations

At each intermediate iteration of the backward recursion, we are dealing with a parametric, piecewise linear program, whose analytic solution is, generally, more intricate than that of the parametric linear program of the first iteration. However, proving that the resulting function $E_{k:n}^{*}(t_{k})$ is CPL, provided that $E_{k+1:n}^{*}(t_{k+1})$ is CPL, is not too complicated. Below, we present a proof which can be summarised as follows. First, for each subproblem $j\in\mathcal{S}_{k+1}$ corresponding to the $j$’th linear piece of the previous CPL function $E_{k+1:n}^{*}(t_{k+1})$, we explain that the corners (or possibly edges) of the feasible region that may represent the optimal solution yield a CPL function in $t_{k}$ when inserted into the objective function $E_{k:n}^{(j)}(t_{k},q_{k})$. Second, we argue that since the boundary of the feasible region evolves in a continuous way as a function of $t_{k}$ and since also $E_{k:n}^{(j)}(t_{k},q_{k})$ is continuous in $t_{k}$ and $q_{k}$, any infinitesimal change in $t_{k}$ can only induce an infinitesimal change in the location of the optimal solution. Third, we conclude from these observations that $\tilde{E}_{k:n}^{(j)}(t_{k})$ is CPL for each subproblem $j\in\mathcal{S}_{k+1}$. This means that the final function $E^{*}_{k:n}(t_{k})$ is the maximum of multiple CPL functions. Therefore, $E^{*}_{k:n}(t_{k})$ itself must be piecewise linear. The additional fact that $E^{*}_{k:n}(t_{k})$ is continuous is an immediate consequence of the continuity of the whole optimisation problem and the connection between the subproblems.

As in the first backward step, the equality constraints (4.25) for $q_{k}$ allow us to reformulate the optimisation problem in terms of the two variables $q^{00}_{k}$ and $q^{10}_{k}$. Specifically, for each subproblem $j\in\mathcal{S}_{k+1}$, we can use the equality constraints to write the objective function $E_{k:n}^{(j)}(t_{k},q_{k})$ cf. (4.29) in terms of $q^{00}_{k}$ and $q^{10}_{k}$ as

[TABLE]

where

[TABLE]

and

[TABLE]

The corresponding constraints for $q^{00}_{k}$ and $q^{10}_{k}$ read

[TABLE]

and

[TABLE]

If for some fixed $t_{k}\in[t_{k}^{\min},t_{k}^{\max}]$ we consider a coordinate system with $q_{k}^{00}$ along one axis and $q_{k}^{10}$ along the other, we see that the feasible region formed by the constraints (A.9)-(A.11) is a polygon with maximally six corners. The region is enclosed by two lines in the $q_{k}^{00}$-direction cf. (A.9), two lines in the $q_{k}^{10}$-direction cf. (A.10), and two parallel lines with a negative slope of $-\pi_{k}^{00}(t_{k})/{\pi_{k}^{10}}(t_{k})$ cf. (A.11). Figure 8 illustrates some of the possible shapes that the region can take. Clearly, the optimal solution is located in a corner of the feasible region, possibly along a whole edge.

To understand where along the boundary of the feasible region the optimal solution is located, we note from (A.8) that if $\tilde{\beta}_{k}^{(j)}$ is positive, then $E_{k:n}^{(j)}(t_{k},q_{k})$ is maximised when $q_{k}^{00}$ and $q_{k}^{10}$ are as large as possible, while if $\tilde{\beta}_{k}^{(j)}$ is negative, then $E_{k:n}^{(j)}(t_{k},q_{k})$ is maximised when $q_{k}^{00}$ and $q_{k}^{10}$ are as small as possible. For simplicity, we assume in the following that the feasible region is non-empty. First, consider the case with $\tilde{\beta}_{k}^{(j)}$ positive. Then, we need to check whether or not the upper of the two lines corresponding to the two inequality constraints in (A.11) forms an edge of the feasible region. If this line does not form an edge of the feasible region; see for example the shapes in Figure 8 (a), (c) and (e); we observe that the point $\left(q_{k}^{00(\mathcal{U})}(t_{k}),q_{k}^{10(\mathcal{U})}(t_{k})\right)$, where

[TABLE]

[TABLE]

is a corner. Moreover, this corner represents the optimal solution, since $q_{k}^{00}$ and $q_{k}^{10}$ jointly take their maximal values in this point. Now, if we insert the functions in (A.12) and (A.13) into the objective function $E_{k:n}^{(j)}(t_{k},q_{k})$, we obtain a CPL function in $t_{k}$. Thereby, given that (A.12) and (A.13) represent a corner of the feasible region for all values of $t_{k}$, the resulting function $\tilde{E}_{k:n}^{(j)}(t_{k})$ is CPL in $t_{k}$. If, on the other hand, the upper of the two lines of the constraints (A.11) does represent an edge of the feasible region; see for instance Figures 8(b), (d), (f) and (g); then this whole edge represents the optimal solution. That is, any point along the edge is optimal. This result is due to that the slope of the objective function and the slope of the line for this edge are equal, from which it follows that the objective function takes the same maximal value anywhere along the edge. Now, if we insert $(q_{k}^{00}$, $q_{k}^{10})$-coordinates located on the edge into the objective function $E_{k:n}^{(j)}(t_{k},q_{k})$, we get a function which is constant, and hence CPL, in $t_{k}$. Thereby, given that the edge is part of the feasible region for all values of $t_{k}$, the resulting function $\tilde{E}_{k:n}^{(j)}(t_{k})$ is CPL in $t_{k}$. Next, consider the case with $\tilde{\beta}_{k}^{(j)}$ negative. Then, the situation is equivalent to the case with $\tilde{\beta}_{k}^{(j)}$ positive, but we need to consider the lower part of the feasible region instead of the upper. That is, we need to check whether or not the lower of the two lines corresponding to the constraints in (A.11) forms an edge of the feasible region. If this line does not represent an edge; see for example Figures 8(a), (d) and (f); the optimal solution is found in the lower left corner point, $\left(q_{k}^{00(\mathcal{L})}(t_{k}),q_{k}^{10(\mathcal{L})}(t_{k})\right)$, where

[TABLE]

[TABLE]

Again, if we insert the functions in (A.14) and (A.15) into the objective function $E_{k:n}^{(j)}(t_{k},q_{k})$, we obtain a CPL function in $t_{k}$. Thereby, given that (A.14) and (A.15) represent a corner of the feasible region for all values of $t_{k}$, the resulting function $\tilde{E}_{k:n}^{(j)}(t_{k})$ is CPL in $t_{k}$. If, on the other hand, the lower of the two lines of the constraints (A.11) does represent an edge of the feasible region, then this edge also represents the optimal solution since the objective function takes the same maximal value anywhere along this edge. Now, if we insert $(q_{k}^{00}$, $q_{k}^{10})$-coordinates located on the optimal edge into the objective function $E_{k:n}^{(j)}(t_{k},q_{k})$, we obtain a function which is constant, and hence CPL, in $t_{k}$. Thereby, given that the edge is part of the feasible region for all values of $t_{k}$, the resulting function $\tilde{E}_{k:n}^{(j)}(t_{k})$ is CPL in $t_{k}$.

Because the objective function, $E_{k:n}^{(j)}(t_{k},q_{k})$, as well as all the constraints, (A.9)-(A.11), are continuous in $t_{k}$ and $q_{k}$, it follows that any infinitesimal change $\delta t_{k}$ in $t_{k}$ can only induce corresponding infinitesimal changes in the shape of the feasible region and the value of the objective function. Hence, the optimal solution at any $t_{k}$-value $t_{k}^{\prime}$ must be located in the same corner (or along the same edge) as the optimal solution at the $t_{k}$-value $t_{k}^{\prime}+\delta t_{k}$. We note, however, that it is possible that the infinitesimal change $\delta t_{k}$ may have added or deleted an edge from the region. In this case, it is possible that a single corner represented the optimal solution at $t_{k}^{\prime}$, while a whole edge represents the optimal solution at $t_{k}^{\prime}+\delta t_{k}$, or vice versa. However, this will not cause any discontinuities in the resulting function $\tilde{E}_{k:n}^{(j)}(t_{k})$ because of the continuity of the optimisation problem as a whole. We have already showed that the coordinates describing the evolution of every potentially optimal corner (or edge) as a function of $t_{k}$ return a CPL function in $t_{k}$. Hence, we understand that $\tilde{E}_{k:n}^{(j)}(t_{k})$ must be CPL.

Finally, we obtain the function $E_{k:n}^{*}(t_{k})$ by taking the maximum of the $\tilde{E}_{k:n}^{(j)}(t_{k})$’s. Taking the maximum of a set of continuous piecewise linear functions necessarily produces another piecewise linear, but not necessarily a continuous, function. However, it is obvious without a further proof that $E_{k:n}^{*}(t_{k})$ must be continuous, since all functions in the whole optimisation problem are continuous. Thereby, we can conclude that $E_{k:n}^{*}(t_{k})$ is CPL.

According to numerical experiments, it seems that $q_{k}^{*00}(t_{k})$ and $q_{k}^{*10}(t_{k})$ are analytically given as $q_{k}^{*00}(t_{k})=q_{k}^{00(\mathcal{U})}(t_{k})$ and $q_{k}^{*10}(t_{k})=q_{k}^{10(\mathcal{U})}(t_{k})$, just as in the first backward iteration. However, we have not proved this result, since it is not really important for our application. Yet, we note that if this result can be proved, the computation of $q(\tilde{x}|x,y)$ can be done particularly simple.

A.3 Computing the breakpoints of $E_{k:n}^{*}(t_{k})$

This section concerns computation of the breakpoints of the CPL function $E_{k:n}^{*}(t_{k})$ produced at each intermediate iteration $2\leq k<n$. The breakpoints of $E_{k:n}^{*}(t_{k})$ should be computed prior to solving the corresponding parametric piecewise linear program in order to avoid numerical computation of $E_{k:n}^{*}(t_{k})$ on a grid of $t_{k}$-values. However, it can in some cases be a bit cumbersome and technical to compute the explicit set of $t_{k}$-values representing the breakpoints of $E_{k:n}^{*}(t_{k})$. Fortunately, it is an easier task to compute a slightly larger set of $t_{k}$-values representing potential breakpoints of $E_{k:n}^{*}(t_{k})$, which necessarily includes all of the actual breakpoints. For convenience, we denote in the following the set of actual breakpoints by $A_{k}$ and the larger set of potential breakpoints by $A_{k}^{\prime}\supset A_{k}$. Having computed the set $A_{k}^{\prime}$, we can solve our parametric piecewise linear program for the $t_{k}$-values in this set, and afterwards go through the values of the resulting function $E_{k:n}^{*}(t_{k})$ to check which of the elements in $A_{k}^{\prime}$ that represent actual breakpoints that must be stored in $A_{k}$, and which points that can be omitted.

As explained in Section A.3, the function $E_{n}^{*}(t_{n})$ of the first backward iteration consists of maximally three linear pieces, or equivalently, it has maximally two breakpoints in addition to its two endpoints $t_{n}^{\min}$ and $t_{n}^{\max}$. Since at each intermediate iteration we consider a more complicated parametric piecewise linear program, additional breakpoints can occur in $E_{k:n}^{*}(t_{k})$, with the number of possible breakpoints for $E_{k:n}^{*}(t_{k})$ increasing with the number of breakpoints for $E_{k+1:n}^{*}(t_{k})$ computed at the previous step of the recursion. To compute the set $A_{k}^{\prime}$ of potential breakpoints for $E_{k:n}^{*}(t_{k})$, we need to check for which $t_{k}$-values the corners of the rectangular region formed by the constraints (A.9) and (A.10) intersect with the lines of the constraints (A.11) for each $j\in\mathcal{S}_{k+1}$. Each $t_{k}$-value that causes such an intersection must be included in the set $A_{k}^{\prime}$. To understand why, consider a subproblem $j\in\mathcal{S}_{k+1}$, and assume $\tilde{\beta}^{(j)}_{k}$ is positive. Furthermore, suppose that for all $t_{k}\in[t_{k}^{\min},t_{k}^{\max}]$ the feasible region has a rectangular shape cf. Figure 8(a), meaning that the region is only enclosed by the constraints (A.9) and (A.10), while the extra constraints in (A.11) do not contribute to shaping the region. Then, from Section A.2, we know that the optimal solution lies in the upper right corner given by (A.12) and (A.13) for all $t_{k}$. Moreover, we know that $\tilde{E}_{k:n}^{(j)}(t_{k})$ is CPL with breakpoints corresponding to the breakpoints of (A.12) and (A.13). Now, suppose instead that after some specific value $t_{k}^{\prime}$ the shape of the feasible region changes from a rectangular shape as in Figure 8(a) to a pentagon shape as in Figure 8(f). This means that the upper of the two lines formed by the extra constraints in (A.11) at the $t_{k}$-value $t_{k}^{\prime}$ intersects with the upper right corner point given by (A.12) and (A.13), while for $t_{k}>t_{k}^{\prime}$ the constraints results in that an extra edge is added to the feasible region. From Section A.2, we then know that for $t_{k}>t_{k}^{\prime}$ this extra edge represents the optimal solution and the value of the objective function remains constant as a function of $t_{k}>t_{k}^{\prime}$. Thereby, we understand that a breakpoint may occur in $\tilde{E}_{k:n}^{(j)}(t_{k})$, and hence possibly in $E_{k:n}^{*}(t_{k})$, at the $t_{k}$-value $t_{k}^{\prime}$. If the feasible region were to evolve in a different way than the one considered here, similar arguments can be formulated. In $A_{k}^{\prime}$, we must also include the breakpoints of the functions in (A.12)-(A.15), i.e. the breakpoints of the functions describing the coordinates for the lower left and upper right corner points of the feasible region when the constraints (A.11) do not contribute.