A New Smoothing Technique based on the Parallel Concatenation of Forward/Backward Bayesian Filters: Turbo Smoothing

TL;DR

This paper introduces turbo smoothing, a new method combining forward and backward Bayesian filters via parallel concatenation, improving complexity-accuracy tradeoff in smoothing for linear Gaussian systems.

Contribution

It extends turbo filtering concepts to smoothing, proposing algorithms that outperform recent methods in complexity and accuracy for certain systems.

Findings

Achieves better complexity-accuracy tradeoff

Demonstrates improved performance over recent smoothing techniques

Validates algorithms with numerical results on linear Gaussian systems

Abstract

Recently, a novel method for developing filtering algorithms, based on the parallel concatenation of Bayesian filters and called turbo filtering, has been proposed. In this manuscript we show how the same conceptual approach can be exploited to devise a new smoothing method, called turbo smoothing. A turbo smoother combines a turbo filter, employed in its forward pass, with the parallel concatenation of two backward information filters used in its backward pass. As a specific application of our general theory, a detailed derivation of two turbo smoothing algorithms for conditionally linear Gaussian systems is illustrated. Numerical results for a specific dynamic system evidence that these algorithms can achieve a better complexity-accuracy tradeoff than other smoothing techniques recently appeared in the literature.

Figures

Equations

$${\mathbf{x}}_{l+1}={\mathbf{f}}_l\left({\mathbf{x}}_l\right)+{\mathbf{w}}_l$$

$${\mathbf{y}}_l$$

$${\mathbf{x}}_{l+1}={\mathbf{F}}_l{\mathbf{x}}_l+{\mathbf{u}}_l+{\mathbf{w}}_l$$

$${\mathbf{y}}_l={\mathbf{H}}_l^T{\mathbf{x}}_l+{\mathbf{v}}_l+{\mathbf{e}}_l,$$

$${\mathbf{x}}_{l+1}^{(Z)}={\mathbf{A}}_l^{(Z)}\left({\mathbf{x}}_l^{(N)}\right){\mathbf{x}}_l^{(L)}+{\mathbf{f}}_l^{(Z)}\left({\mathbf{x}}_l^{(N)}\right)+{\mathbf{w}}_l^{(Z)}$$

$${\mathbf{y}}_l={\mathbf{g}}_l\left({\mathbf{x}}_l^{(N)}\right)+{\mathbf{B}}_l\left({\mathbf{x}}_l^{(N)}\right){\mathbf{x}}_l^{(L)}+{\mathbf{e}}_l$$

$${\vec{m}}_{fp}\left({\mathbf{x}}_l\right)\triangleq f\left({\mathbf{x}}_l,{\mathbf{y}}_{1:(l-1)}\right),$$

$${\overleftarrow{m}}_{be}\left({\mathbf{x}}_{l+1}\right)\triangleq f\left({\mathbf{y}}_{(l+1):T}|{\mathbf{x}}_{l+1}\right),$$

$$f\left({\mathbf{x}}_l,{\mathbf{y}}_{1:T}\right)={\vec{m}}_{fp}\left({\mathbf{x}}_l\right){\overleftarrow{m}}_{be}\left({\mathbf{x}}_l\right)$$

$$f\left({\mathbf{x}}_l,{\mathbf{y}}_{1:T}\right)={\vec{m}}_{fe}\left({\mathbf{x}}_l\right){\overleftarrow{m}}_{bp}\left({\mathbf{x}}_l\right),$$

$$f\left({\mathbf{x}}_l,{\mathbf{y}}_{1:T}\right)$$

$$=$$

$${\mathbf{z}}_l^{(N)}\triangleq {\mathbf{x}}_{l+1}^{(L)}-{\mathbf{A}}_l^{(L)}\left({\mathbf{x}}_l^{(N)}\right){\mathbf{x}}_l^{(L)}\text{,}$$

$${\mathbf{z}}_l^{(N)}={\mathbf{f}}_l^{(L)}\left({\mathbf{x}}_l^{(N)}\right)+{\mathbf{w}}_l^{(L)}$$

$${\mathbf{z}}_l^{(L)}\triangleq {\mathbf{x}}_{l+1}^{(N)}-{\mathbf{f}}_l^{(N)}\left({\mathbf{x}}_l^{(N)}\right);$$

$${\mathbf{z}}_l^{(L)}={\mathbf{A}}_l^{(N)}\left({\mathbf{x}}_l^{(N)}\right){\mathbf{x}}_l^{(L)}+{\mathbf{w}}_l^{(N)}$$

$${\vec{m}}_{fp}\left({\mathbf{x}}_l\right)\triangleq \mathcal{N}\left({\mathbf{x}}_l;{\eta }_{fp,l},{\mathbf{C}}_{fp,l}\right).$$

$${\vec{m}}_{fe1}\left({\mathbf{x}}_l\right)\triangleq \mathcal{N}\left({\mathbf{x}}_l;{\eta }_{fe1,l},{\mathbf{C}}_{fe1,l}\right),$$

$${\mathbf{W}}_{fe1,l}={\mathbf{H}}_l{\mathbf{W}}_e{\mathbf{H}}_l^T+{\mathbf{W}}_{fp,l}$$

$${\mathbf{w}}_{fe1,l}={\mathbf{H}}_l{\mathbf{W}}_e\left({\mathbf{y}}_l-{\mathbf{v}}_l\right)+{\mathbf{w}}_{fp,l},$$

$${\vec{m}}_{fp,j}\left({\mathbf{x}}_l^{(N)}\right)\triangleq \delta \left({\mathbf{x}}_l^{(N)}-{\mathbf{x}}_{fp,l,j}^{(N)}\right),$$

$${\vec{m}}_{fe1,j}\left({\mathbf{x}}_l^{(N)}\right)\triangleq w_{fe,l,j}\delta \left({\mathbf{x}}_l^{(N)}-{\mathbf{x}}_{fp,l,j}^{(N)}\right).$$

$${\overleftarrow{m}}_{be}\left({\mathbf{x}}_{l+1}\right)\triangleq \mathcal{N}\left({\mathbf{x}}_{l+1};{\eta }_{be,l+1},{\mathbf{C}}_{be,l+1}\right)$$

$${\overleftarrow{m}}_{be}\left({\mathbf{x}}_{l+1}^{(N)}\right)\triangleq \delta \left({\mathbf{x}}_{l+1}^{(N)}-{\mathbf{x}}_{be,l+1}^{(N)}\right),$$

$${\overleftarrow{m}}_{bp}\left({\mathbf{x}}_l\right)\triangleq \mathcal{N}\left({\mathbf{x}}_l;{\eta }_{bp,l},{\mathbf{C}}_{bp,l}\right)$$

$${\overleftarrow{m}}_{bp,j}\left({\mathbf{x}}_l^{(N)}\right)\triangleq w_{bp,l,j}$$

$${\overleftarrow{m}}_{be1}\left({\mathbf{x}}_l\right)\triangleq \mathcal{N}\left({\mathbf{x}}_l;{\eta }_{be1,l},{\mathbf{C}}_{be1,l}\right)$$

$${\overleftarrow{m}}_{be2}\left({\mathbf{x}}_l\right)\triangleq \mathcal{N}\left({\mathbf{x}}_l;{\eta }_{be2,l},{\mathbf{C}}_{be2,l}\right)={\overleftarrow{m}}_{be}\left({\mathbf{x}}_l\right),$$

$${\overleftarrow{m}}_{be1,j}\left({\mathbf{x}}_l^{(N)}\right)\triangleq w_{be1,l,j}$$

$${\overleftarrow{m}}_{be}\left({\mathbf{x}}_l^{(N)}\right)=\delta \left({\mathbf{x}}_l^{(N)}-{\mathbf{x}}_{be,l}^{(N)}\right)\text{,}$$

Full text

Giorgio M. Vitetta, Pasquale Di Viesti and Emilio Sirignano

University of Modena and Reggio Emilia

Department of Engineering ”Enzo Ferrari”

Via P. Vivarelli 10/1, 41125 Modena - Italy

email: giorgio.vitetta@unimore.it, pasquale.diviesti2@unibo.it, emilio.sirignano@unimore.it

Keywords: Hidden Markov Model, Smoothing, Factor Graph, Particle Filter, Kalman Filter, Parallel Concatenation, Sum-Product Algorithm, Turbo Processing.

1 Introduction

The problem of Bayesian smoothing for a state space model (SSM) concerns the development of recursive algorithms able to estimate the probability density function (pdf) of the model state on a given observation interval, given a batch of noisy measurements acquired over it [1]; the estimated pdf is known as a smoothed or smoothing pdf. A general strategy for solving this problem is based on the so called two-filter smoothing formula [2]-[3]; in fact, this formula allows to compute the required smoothing density by merging the statistical information generated in the forward pass of a Bayesian filtering method with those evaluated in the backward pass of a different filtering method, paired with the first one and known as backward information filtering (BIF). Unluckily, closed form solutions for this strategy can be derived for linear Gaussian and linear Gaussian mixture models only [1], [4]. For this reason, all the existing smoothing algorithms based on the above mentioned formula and applicable to general nonlinear models are approximate and are based on sequential Monte Carlo techniques (e.g., see [2], [5], [6] and references therein). Unluckily, the adoption of these algorithms, known as particle smoothers, may be hindered by their complexity, which becomes unmanageable when the dimension of the sample space for the considered SSM is large.

Recently, a factor graph approach has been exploited to devise a new filtering method, based on the parallel concatenation of two (constituent) Bayesian filters and called *turbo filtering *(TF) [7]. In this manuscript, a new smoothing technique that employs TF in its forward pass and a new BIF scheme, based on the parallel concatenation of two backward information filters, is developed. Our derivation of the new BIF method, called backward information turbo filtering (BITF), is based on a general graphical model; this allows us to: a) represent any BITF algorithm as the interconnection of *two soft-in soft-out *(SISO) processing modules; b) represent the iterative processing accomplished by these modules as a message passing technique; c) derive the expressions of the passed messages by applying the sum-product algorithm (SPA) [8], [9], together with a specific scheduling procedure, to the graphical model itself; d) show how the statistical information generated by a BITF algorithm in the backward pass can be merged with those produced by the paired TF technique in the forward pass in order to evaluate the required smoothed pdfs. To exemplify the usefulness of this smoothing method, called turbo-smoothing (TS) in the following, we take into consideration the TF algorithms proposed in [7] for the class of conditionally linear Gaussian (CLG) SSMs and derive a BITF algorithm paired with them. This approach leads to the development of two new TS algorithms, one generating an estimate of the joint smoothing density over the whole observation interval, the other one an estimate of the marginal smoothing densities over the same interval. Our computer simulations for a specific CLG SSM evidence that, in the considered case, the derived TS algorithms perform very closely to the Rao-Blackwellized particle smoothing (RBPS) technique proposed in [10] and to the particle smoothers devised in [11].

The remaining part of this manuscript is organized as follows. A description of the considered SSMs is illustrated in Section 2. In Section 3, a general graphical model on which the processing accomplished in BITF and TS is based is illustrated; then, a specific instance of it, referring to a CLG SSM, is developed and the messages passed over it in BITF are defined. In Section 4, the scheduling and the computation of such messages are described, specific TS algorithms are developed, and the differences and similarities between these algorithms and other smoothing techniques are briefly analysed. A comparison, in terms of accuracy and execution time, between the proposed techniques and three smoothers recently appeared in the literature is provided in Section 5 for a specific CLG SSM. Finally, some conclusions are offered in Section 6.

Notations: The same notations as refs. [11], [7] and [12] are adopted.

2 Model Description

In this manuscript we focus on a discrete-time SSM whose $D$-dimensional hidden state in the $l$-th interval is denoted $\mathbf{x}_{l}\triangleq[x_{0,l},x_{1,l},...,$ $x_{D-1,l}]^{T}$, and whose state update and measurement models are expressed by

[TABLE]

and

[TABLE]

respectively. Here, $\mathbf{f}_{l}\left(\mathbf{x}_{l}\right)$ ($\mathbf{h}_{l}\left(\mathbf{x}_{l}\right)$) is a time-varying $D$-dimensional ($P$-dimensional) real function and $\mathbf{w}_{l}$ ($\mathbf{e}_{l}$) the $l$-th element of the process (measurement) noise sequence $\left\{\mathbf{w}_{k}\right\}$ ($\left\{\mathbf{e}_{k}\right\}$); this sequence consists of $D$-dimensional ($P$-dimensional) independent and identically distributed (iid) Gaussian noise vectors, each characterized by a zero mean and a covariance matrix $\mathbf{C}_{w}$ ($\mathbf{C}_{e}$). Moreover, statistical independence between $\left\{\mathbf{e}_{k}\right\}$ and $\{\mathbf{w}_{k}\}$ is assumed.

In the following, two additional mathematical representations for the considered SSM are also exploited. The first one is approximate, being employed by an extended Kalman filter (EKF); in fact, it is based on the linearized versions of eqs. (1) and (2), namely (e.g., see [1, pp. 194-195])

[TABLE]

and

[TABLE]

respectively; here, $\mathbf{F}_{l}\triangleq[\partial\mathbf{f}_{l}\left(\mathbf{x}\right)/\partial\mathbf{x}]_{\mathbf{x=x}_{fe,l}}$, $\mathbf{x}_{fe,l}$ is the (forward) estimate of $\mathbf{x}_{l}$ evaluated by the EKF in its $l$-th recursion, $\mathbf{u}_{l}\triangleq\mathbf{f}_{l}\left(\mathbf{x}_{fe,l}\right)-\mathbf{F}_{l}\mathbf{x}_{fe,l}$, $\mathbf{H}_{l}^{T}\triangleq[\partial\mathbf{h}_{l}\left(\mathbf{x}\right)/\partial\mathbf{x}]_{\mathbf{x=x}_{fp,l}}$, $\mathbf{x}_{fp,l}$ is the (forward) prediction $\mathbf{x}_{l}$ computed by the EKF in its $(l-1)$-th recursion and $\mathbf{v}_{l}\triangleq\mathbf{h}_{l}\left(\mathbf{x}_{fp,l}\right)-\mathbf{H}_{l}^{T}\mathbf{x}_{fp,l}$.

The second representation is based on the additional assumption that the SSM described by eqs. (1)-(2) is CLG [10], [13], so that its state vector in the $l$-th interval can be partitioned as $\mathbf{x}_{l}=[(\mathbf{x}_{l}^{(L)})^{T},(\mathbf{x}_{l}^{(N)})^{T}]^{T}$; here, $\mathbf{x}_{l}^{(L)}\triangleq[x_{0,l}^{(L)}$, $x_{1,l}^{(L)},...,x_{D_{L}-1,l}^{(L)}]^{T}$ ($\mathbf{x}_{l}^{(N)}\triangleq[x_{0,l}^{(N)},x_{1,l}^{(N)},...,x_{D_{N}-1,l}^{(N)}]^{T}$) is the so called *linear *(nonlinear) component of $\mathbf{x}_{l}$, with $D_{L}<D$ ($D_{N}=D-D_{L}$). For this reason, following [11], [12] and [13], the models

[TABLE]

and

[TABLE]

are adopted for the update of the linear ($Z=L$) and nonlinear ($Z=N$) components and for the measurement vector, respectively. In the state update model (5), $\mathbf{f}_{l}^{(Z)}(\mathbf{x}_{l}^{(N)})$ ($\mathbf{A}_{l}^{(Z)}(\mathbf{x}_{l}^{(N)})$) is a time-varying $D_{Z}$-dimensional real function ($D_{Z}\times D_{L}$ real matrix) and $\mathbf{w}_{l}^{(Z)}$ consists of the first $D_{L}$ (last $D_{N}$) elements of $\mathbf{w}_{l}$ if $Z=L$ ($Z=N$); independence between $\{\mathbf{w}_{k}^{(L)}\}$ and $\{\mathbf{w}_{k}^{(N)}\}$ is also assumed for simplicity and the covariance matrix $\mathbf{w}_{k}^{(L)}$ ($\mathbf{w}_{k}^{(N)}$) is denoted $\mathbf{C}_{w}^{(L)}$ ($\mathbf{C}_{w}^{(N)}$). In the measurement model (6), instead, $\mathbf{g}_{l}(\mathbf{x}_{l}^{(N)})$ ($\mathbf{B}_{l}(\mathbf{x}_{l}^{(N)})$) is a time-varying $P$-dimensional real function ($P\times D_{L}$ real matrix).

In the next two Sections we focus on the problem of developing algorithms for the estimation of a) the joint smoothed pdf $f(\mathbf{x}_{1:T}|\mathbf{y}_{1:T})$ (problem P.1) and b) the sequence of marginal *smoothed pdfs *$\{f(\mathbf{x}_{l}|\mathbf{y}_{1:T}),\,l=1,2,...,T\}$ (problem P.2); here, $T$ is the duration of the observation interval and $\mathbf{y}_{1:T}=\left[\mathbf{y}_{1}^{T},\mathbf{y}_{2}^{T},...,\mathbf{y}_{T}^{T}\right]^{T}$ is a $T\cdot P$-dimensional vector. It is important to point out that: a) in solving both problems P.1 and P.2, the prior knowledge of the pdf $f(\mathbf{x}_{1})$ of the initial state is assumed; b) in principle, if an estimate of the joint pdf $f(\mathbf{x}_{1:T}|\mathbf{y}_{1:T})$ is available, estimates of all the posterior $\{f(\mathbf{x}_{l}|\mathbf{y}_{1:T})\}$ can be evaluated by marginalization.

3 Graphical Modelling for the Parallel Concatenation of Bayesian

Information Filters

In this Section, we derive the graphical models on which BITF and TS techniques are based. More specifically, starting from the factor graph representing Bayesian smoothing [11], we first develop a general graphical model for the parallel concatenation of two backward information filters. Then, a specific instance of this model is devised for the case in which the forward filters are an EKF and a particle filter (PF), and the considered SSM is CLG.

3.1 Graphical Model for the Parallel Concatenation of Bayesian

Information Filters and Message Passing over it

The development of our BIF algorithms is based on the graphical approach illustrated in ref. [11, Sec. III]. This approach consists in representing Bayesian filtering and BIF as two recursive algorithms that compute, on the basis of the SPA, a set of probabilistic messages passed along the same (cycle free) factor graph; this graph is illustrated Fig. 1-a) and refers to a SSM characterized by the Markov model $f(\mathbf{x}_{l+1}|\mathbf{x}_{l})$ and the measurement model $f(\mathbf{y}_{l}|\mathbf{x}_{l})$. More specifically, in the $l$-th recursion of Bayesian filtering*, *messages are passed along the considered graph in the forward direction; moreover, the messages $\vec{m}_{fe}\left(\mathbf{x}_{l}\right)=f\left(\mathbf{x}_{l},\mathbf{y}_{1:l}\right)$ and $\vec{m}_{fp}\left(\mathbf{x}_{l+1}\right)=f(\mathbf{x}_{l+1},\mathbf{y}_{1:l})$ (denoted $FE_{l}$ and $FP_{l+1}$, respectively, in Fig. 1 and conveying a forward estimate of $\mathbf{x}_{l}$ and a forward prediction of $\mathbf{x}_{l+1}$, respectively) are computed on the basis of the input message

[TABLE]

for $l=2$, $3$, $...$, $T$ ($\vec{m}_{fp}\left(\mathbf{x}_{1}\right)=f(\mathbf{x}_{1})$). Dually, in the $(T-l)$-th recursion of BIF, messages are passed along the considered graph in the backward direction, and the messages $\overset{\leftarrow}{m}_{bp}\left(\mathbf{x}_{l}\right)=f(\mathbf{y}_{(l+1):T}|\mathbf{x}_{l})$ and $\overset{\leftarrow}{m}_{be}\left(\mathbf{x}_{l}\right)=f(\mathbf{y}_{l:T}|\mathbf{x}_{l})$ (denoted $BP_{l}$ and $BE_{l}$, respectively, in Fig. 1 and conveying a backward prediction of $\mathbf{x}_{l}$ and a backward estimate of $\mathbf{x}_{l+1}$, respectively) are computed on the basis of the input message

[TABLE]

with $l=T-2,T-3,...,1$ ($\overset{\leftarrow}{m}_{be}\left(\mathbf{x}_{T}\right)=f\left(\mathbf{y}_{T}|\mathbf{x}_{T}\right)$). Once the backward pass is over, a solution to problem P.2 becomes available, since the marginal smoothed pdf $f\left(\mathbf{x}_{l},\mathbf{y}_{1:T}\right)$ can be evaluated as111Note that, similarly as [11] and [12], the joint pdf $f(\mathbf{x}_{l},\mathbf{y}_{1:T})$ is considered here in place of the posterior pdf $f(\mathbf{x}_{l}|\mathbf{y}_{1:T})$.

[TABLE]

or, equivalently, as

[TABLE]

with $l=1,2,...,T$. Note that, from a graphical viewpoint, formulas (9) and (10) can be related with the two different partitionings of the graph shown in Fig. 1-a) (where a specific partitioning is identified by a brown dashed vertical line cutting the graph in two parts).

In ref. [7] it has been also shown that the factor graph illustrated in Fig. 1-a) can be employed as a building block in the development of a larger graphical model that represents a turbo filtering scheme, i.e. the parallel concatenation of two (constituent) Bayesian filters (denoted F1 and F2 in the following). In this model, the graphs referring to F1 and F2 are interconnected in order to allow the mutual exchange of statistical information in the form of pseudo-measurements (conveyed by probabilistic messages). From a graphical viewpoint, the exploitation of these additional information in each filter requires:

a) modifying the graph shown in Fig. 1-a) in a way that each constituent filter can benefit from the pseudo-measurements provided by the other filter through an additional measurement update;

b) developing message passing algorithms over a proper graphical model for

1. the conversion of the statistical information generated by each constituent filter into a form useful to the other one and 2) the generation, inside each constituent filter, of the statistical information to be made available to the other filter.

As far as the need expressed at point a) is concerned, the graph of Fig. 1-a) can be easily modified by adding a new equality node and a new edge along which the message $m_{pm}\left(\mathbf{x}_{l}\right)$, conveying pseudo-measurement information, is passed; this results in the factor graph shown in Fig. 1-b). Note that, in the new graphical model, two forward estimates (backward estimates) are computed in the forward (backward) pass. The first estimate, represented by $\vec{m}_{fe1}\left(\mathbf{x}_{l}\right)$ ($\overset{\leftarrow}{m}_{be1}\left(\mathbf{x}_{l}\right)$) is generated by merging $\vec{m}_{fp}\left(\mathbf{x}_{l}\right)$ ($\overset{\leftarrow}{m}_{bp}\left(\mathbf{x}_{l}\right)$) with the message $m_{ms}\left(\mathbf{x}_{l}\right)$ ($m_{pm}\left(\mathbf{x}_{l}\right)$) conveying measurement (pseudo-measurement) information, whereas the second one, represented by $\vec{m}_{fe2}\left(\mathbf{x}_{l}\right)$ ($\overset{\leftarrow}{m}_{be2}\left(\mathbf{x}_{l}\right)=\overset{\leftarrow}{m}_{be}\left(\mathbf{x}_{l}\right)$), is evaluated by merging $\vec{m}_{fe1}\left(\mathbf{x}_{l}\right)$ ($\overset{\leftarrow}{m}_{be1}\left(\mathbf{x}_{l}\right)$) with the message $m_{pm}\left(\mathbf{x}_{l}\right)$ ($m_{ms}\left(\mathbf{x}_{l}\right)$). Moreover, similarly as the previous case, the smoothed pdf $f\left(\mathbf{x}_{l},\mathbf{y}_{1:T}\right)$ can be computed as

[TABLE]

note also that each of these factorisations can be associated with one of the three distinct vertical cuts drawn in Fig. 1-b).

As far as point b) is concerned, in ref. [7] it is shown that, in any TF scheme, all the processing tasks related to the conversion (generation) of the statistical information emerging from (feeding) each constituent filter can be easily incorporated in a single* module*, called soft-in soft-out (SISO) module and whose overall processing can be represented as message passing over a graphical model including the factor graph shown in Fig. 1-b). For this reason, any TF scheme can be devised by linking (i.e., by concatenating) two SISO modules, each incorporating a specific filtering algorithm and exchanging probabilistic information in an iterative fashion. It is also important to point out that the two constituent filters are not required to estimate the whole system state. For this reason, in the following, we assume that: a) the filter Fi estimates the portion $\mathbf{x}_{l}^{(i)}$ (with $i=1$ and $2$) of the state vector $\mathbf{x}_{l}$ (the size of $\mathbf{x}_{l}^{(i)}$ is denoted $D_{i}$, with $D_{i}\leq D$); b) the portion of $\mathbf{x}_{l}$ not included in $\mathbf{x}_{l}^{(i)}$ is denoted $\mathbf{\bar{x}}_{l}^{(i)}$, so that the equalities $\mathbf{x}_{l}=[(\mathbf{x}_{l}^{(1)})^{T},(\mathbf{\bar{x}}_{l}^{(1)})^{T}]^{T}$ or $\mathbf{x}_{l}=[(\mathbf{\bar{x}}_{l}^{(2)})^{T},(\mathbf{x}_{l}^{(2)})^{T}]^{T}$ hold. However, the vector $\mathbf{\bar{x}}_{l}^{(1)}$ ($\mathbf{\bar{x}}_{l}^{(2)}$) is required to be part of (or, at most, to coincide with) $\mathbf{x}_{l}^{(2)}$ ($\mathbf{x}_{l}^{(1)}$), so that an overall estimate of the system state $\mathbf{x}_{l}$ can be always generated on the basis of the posterior pdfs of $\mathbf{x}_{l}^{(1)}$ and $\mathbf{x}_{l}^{(2)}$ evaluated by F1 and F2, respectively. In fact, this constraint on $\mathbf{\bar{x}}_{l}^{(1)}$ and $\mathbf{\bar{x}}_{l}^{(2)}$ leads to the conclusion that, generally speaking, the portion $\mathbf{x}_{l}^{(12)}=[x_{D-D_{2},l},x_{D-D_{2}+1,l},...,x_{D_{1}-1,l}]^{T}$ of $\mathbf{x}_{l}$, collecting $N_{d}\triangleq D_{1}+D_{2}-D$ state variables, is estimated by both F1 and F2, being shared by $\mathbf{x}_{l}^{(1)}$ and $\mathbf{x}_{l}^{(2)}$.

A similar conceptual approach is followed in the remaining part of this Paragraph to derive the general representation of the BIF technique paired with a given TF scheme, that is, briefly, a backward information turbo filtering (BITF) technique. This means that:

1. The general architecture we propose for BITF is based on the parallel concatenation of two constituent Bayesian information filters, that are denoted BIF1 and BIF2 in the following.

2. The processing accomplished by BIF1 (BIF2) is represented as a message passing algorithm over the same graphical model as F1 (F2).

3. BITF processing can be represented as the iterative exchange of probabilistic information between two distinct SISO modules.

4. The $i$-th SISO module (with $i=1$ and $2$) incorporates a specific BIF algorithm, that can be represented as a message passing over a factor graph similar to that shown in Fig. 1-b) and that estimates the portion $\mathbf{x}_{l}^{(i)}$ of $\mathbf{x}_{l}$.

The graphical model developed for the SISO module based on BIF1 is shown in Fig. 2. In this Figure, to ease the interpretation of message passing, three rectangles, labeled as BIF1-IN, BIF1 and BIF1-OUT, have been drawn; this allow us to easily identify the portions of the graphical model involved in a) the conversion of the statistical information provided from BIF2 into a form useful to BIF1, b) BIF1 processing and c) the generation of the statistical information made available by BIF1 to BIF2, respectively.  A detailed description of the signal processing tasks accomplished within each portion is provided below.

BIF1-IN - The statistical information provided by BIF2 to the considered SISO module is condensed in the messages $m_{sm}(\mathbf{x}_{l}^{(2)})$ and $m_{pm}(\mathbf{\bar{x}}_{l}^{(2)})$; these convey a smoothed estimate of $\mathbf{x}_{l}^{(2)}$ and pseudo-measurement information about $\mathbf{\bar{x}}_{l}^{(2)}$, respectively. The first message is processed in two different ways. In fact, on the one hand, it is marginalised in the block labelled by the letter M (see Fig. 2) in order to generate the pdf $m_{sm}(\mathbf{\bar{x}}_{l}^{(1)})$ (do not forget that the state vector $\mathbf{\bar{x}}_{l}^{(1)}$ is included in $\mathbf{x}_{l}^{(2)}$); on the other hand, $m_{sm}(\mathbf{x}_{l}^{(2)})$ is processed jointly with $m_{pm}(\mathbf{\bar{x}}_{l}^{(2)})$ in order to generate the message $m_{pm}(\mathbf{x}_{l}^{(1)})$ conveying pseudo-measurement information about $\mathbf{x}_{l}^{(1)}$ (this is accomplished in the block called PM* conversion*, PMC; see Fig. 2). Then, the messages $m_{sm}(\mathbf{\bar{x}}_{l}^{(1)})$ and $m_{pm}(\mathbf{x}_{l}^{(1)})$ are passed to BIF1.

BIF1 - The message passing accomplished in this part refers to the BIF algorithm paired with F1. The graphical model developed for it and the message passing accomplished over it are based on Fig. 1-b). Note, however, that: a) the message passing aims at computing the (backward) predicted density $\overset{\leftarrow}{m}_{bp}(\mathbf{x}_{l}^{(1)})$ and the (backward) filtered density $\overset{\leftarrow}{m}_{be2}(\mathbf{x}_{l}^{(1)})=\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(1)})$ and on the basis of the backward estimate $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(1)})$ originating from the previous recursion, and of the messages $m_{sm}(\mathbf{\bar{x}}_{l}^{(1)})$ and $m_{pm}(\mathbf{x}_{l}^{(1)})$ provided by BIF1-IN; b) an approximate model of the considered SSM could be adopted in the evaluation of these densities. For this reason, generally speaking, we can assume that the BIF1 algorithm is based on the Markov model $\tilde{f}(\mathbf{x}_{l+1}^{(1)}|\mathbf{x}_{l}^{(1)},\mathbf{\bar{x}}_{l}^{(1)})$ and on the observation model $\tilde{f}(\mathbf{y}_{l}|\mathbf{x}_{l}^{(1)},\mathbf{\bar{x}}_{l}^{(1)})$, representing the exact models $f(\mathbf{x}_{l+1}^{(1)}|\mathbf{x}_{l}^{(1)},\mathbf{\bar{x}}_{l}^{(1)})$ and $f(\mathbf{y}_{l}|\mathbf{x}_{l}^{(1)},\mathbf{\bar{x}}_{l}^{(1)})$, respectively, or approximations of one or both of them. Note also that, in both the second measurement update and the time update accomplished by this algorithm, marginalization with respect to the unknown state component $\mathbf{\bar{x}}_{l}^{(1)}$ is made possible by the availability of the message $m_{sm}(\mathbf{\bar{x}}_{l}^{(1)})$.

BIF1-OUT - This part is fed by the backward estimate $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(1)})$ of $\mathbf{x}_{l+1}^{(1)}$ and by the smoothed estimate $m_{sm}(\mathbf{x}_{l}^{(1)})$ of $\mathbf{x}_{l}^{(1)}$ (available after that the first measurement update has been accomplished by F1). The second message follows two different paths, since a) it is passed to the other SISO module as it is and b) it is jointly processed with $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(1)})$ in order to generate the pseudo-measurement message $m_{pm}(\mathbf{\bar{x}}_{l}^{(1)})$ feeding the other SISO module; the last task is accomplished in the pseudo-measurement generation (PMG) block.

A graphical model structurally identical to the one shown in Fig. 2 can be easily drawn for the SISO module based on BIF2 by interchanging $\mathbf{x}_{l}^{(1)}$ ($\mathbf{\bar{x}}_{l}^{(1)}$) with $\mathbf{x}_{l}^{(2)}$ ($\mathbf{\bar{x}}_{l}^{(2)}$). Merging the graphical model shown in Fig. 2 with its counterpart referring to BIF2 results in the parallel concatenation architecture illustrated in Fig. 3 (details about the underlying graphical model are omitted for simplicity) and on which TS is based. It is important to point out that:

1. The overall graphical model derived for TS, unlike the one illustrated in Fig. 1, is not cycle free; therefore, the application of the SPA to it requires defining a proper message scheduling and, generally speaking, results in iterative algorithms.

2. At the end of the $l$-th recursion of a TS algorithm, two smoothed densities, namely $m_{sm}(\mathbf{x}_{l}^{(1)})$ and $m_{sm}(\mathbf{x}_{l}^{(2)})$, are available. This raises the problem of how these statistical information can be fused in order to get a single pdf for (and, in particular, a single smoothed estimate of) the $N_{d}$-dimensional portion $\mathbf{x}_{l}^{(12)}$ of $\mathbf{x}_{l}$ estimated by both F1/BIF1 and F2/BIF2. Unluckily, this remains an open issue. In our computer simulations, a simple selection strategy has been adopted in state estimation, since one of the two smoothed estimates of $\mathbf{x}_{l}^{(12)}$ has been systematically discarded.

3.2 A Graphical Model for the Parallel Concatenation of the Bayesian

Information Filters Paired with an Extended Kalman Filter and a Particle Filter

In the remaining part of this manuscript we focus on a specific instance of the proposed TS architecture, since we make the same specific choices as [7] for both the SSM and the filters employed in the forward pass. In particular, we focus on the CLG SSM described in Section 2 and assume that:

1. BIF1 is the backward filter associated with an EKF operating over the whole system state (so that $\mathbf{x}_{l}^{(1)}=\mathbf{x}_{l}$ and $\mathbf{\bar{x}}_{l}^{(1)}$ is empty). In other words, BIF1 is a backward Kalman filter based on a linearised model of the considered SSM.

2. BIF2 is a backward filter associated with a PF (in particular, a sequential importance resampling filter [14]) operating on the nonlinear state component only (so that $\mathbf{x}_{l}^{(2)}=\mathbf{x}_{l}^{(N)}$ and $\mathbf{\bar{x}}_{l}^{(2)}=\mathbf{x}_{l}^{(L)}$) and representing it through a set of $N_{p}$ particles (note that $N_{d}=D_{N}$ elements of the system state are shared by the two BIF algorithms). This means that BIF2 is employed to compute new weights for all the elements of the particle set generated by the PF in the forward pass.

Based on the general models shown in Figs. 2 and 3, the specific graphical model illustrated in Fig. 4 (and referring to the $(T-l)$-th recursion of backward filtering) can be drawn for the considered case. In the following, we provide various details about the adopted notation and the message passing within each constituent filter and from each filter to the other one.

Message passing within BIF1 - BIF1 is based on the approximate statistical models $\tilde{f}(\mathbf{x}_{l+1}|\mathbf{x}_{l})$ and $\tilde{f}(\mathbf{y}_{l}|\mathbf{x}_{l})$; these are derived from the linearised eqs. (3) and (4), respectively. Moreover, the (Gaussian) messages passed over its graph (enclosed within the upper rectangle appearing in Fig. 4) are $\vec{m}_{fp}(\mathbf{x}_{l})$, $m_{ms}(\mathbf{x}_{l})$, $\vec{m}_{fe1}(\mathbf{x}_{l})$, $m_{pm}(\mathbf{x}_{l})$, $\overset{\leftarrow}{m}_{be1}(\mathbf{x}_{l})$, $\overset{\leftarrow}{m}_{be2}(\mathbf{x}_{l})$ ($=\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l})$), $\overset{\leftarrow}{m}_{bp}(\mathbf{x}_{l})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1})$, and are denoted $FP$, $MS$, $FE1$, $PM$, $BE1$, $BE2$ ($BE$), $BP$ and $BE^{{}^{\prime}}$, respectively, to ease reading.

Message passing within BIF2 - BIF2 is based on the exact statistical models $f(\mathbf{x}_{l+1}^{(N)}|\mathbf{x}_{l}^{(N)}$, $\mathbf{x}_{l}^{(L)})$ and $f(\mathbf{y}_{l}|\mathbf{x}_{l}^{(N)},\mathbf{x}_{l}^{(L)})$, that are derived from the eqs. (5) (with $Z=N$) and (6), respectively. Moreover, the messages processed by it and appearing in Fig. 4 refer to the $j$-th particle predicted in the previous (i.e., in the $(l-1)$-th) recursion of forward filtering and denoted $\mathbf{x}_{fp,l,j}^{(N)}$, with $j=0,1,...,N_{p}-1$; such messages are $\vec{m}_{fp,j}(\mathbf{x}_{l}^{(N)})$, $m_{ms,j}(\mathbf{x}_{l}^{(N)})$, $\vec{m}_{fe1,j}(\mathbf{x}_{l}^{(N)})$, $m_{pm,j}(\mathbf{x}_{l}^{(N)})$, $\overset{\leftarrow}{m}_{be1,j}(\mathbf{x}_{l}^{(N)})$, $\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l}^{(N)})$, $\overset{\leftarrow}{m}_{bp,j}(\mathbf{x}_{l}^{(N)})$ and $\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l+1}^{(N)})$, and are denoted $FPN_{j}$, $MSN_{j}$, $FE1N_{j}$, $PMN_{j}$, $BE1N_{j}$, $BEN_{j}$, $BPN_{j}$ and $BEN_{j}^{{}^{\prime}}$, respectively, to ease reading.

Message passing from BIF1 to BIF2 - BIF2 is fed by the message $m_{sm}(\mathbf{x}_{l}^{(L)})$ and the message set $\{m_{pm,j}(\mathbf{x}_{l}^{(N)})\}$ conveying pseudo-measurement information; these messages are computed on the basis of the statistical information made available by BIF1. More specifically, on the one hand, the message $m_{sm}(\mathbf{x}_{l}^{(L)})$ (denoted $SML$) results from the marginalization of $m_{sm}(\mathbf{x}_{l})$ and is employed for marginalising the PF state update and measurement models (i.e., $f(\mathbf{x}_{l+1}^{(N)}|\mathbf{x}_{l}^{(N)}$, $\mathbf{x}_{l}^{(L)})$ and $f(\mathbf{y}_{l}|\mathbf{x}_{l}^{(N)},\mathbf{x}_{l}^{(L)})$, respectively) with respect to $\mathbf{x}_{l}^{(L)}$. On the other hand, the pseudo-measurement message $m_{pm,j}(\mathbf{x}_{l}^{(N)})$ (denoted $PMN_{j}$) is evaluated in the PMG1→2 block by processing the messages $m_{sm}(\mathbf{x}_{l}^{(L)})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(L)})$ (denoted $BEL^{\prime}$ and resulting from the marginalization of $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1})$), under the assumption that $\mathbf{x}_{l}^{(N)}$ is represented by the $j$-th particle (conveyed by the message $m_{sm,j}(\mathbf{x}_{l}^{(N)})$).

As illustrated in the Appendix, the computation of the message $m_{pm,j}(\mathbf{x}_{l}^{(N)})$ involves the evaluation of the pdf of the random vector

[TABLE]

defined on the basis of the state update equation (5) (with $Z=L$) and conditioned on the fact that $\mathbf{x}_{l}^{(N)}=\mathbf{x}_{fp,l,j}^{(N)}$. This pdf, which is computed according to the joint statistical characterization of $\mathbf{x}_{l}^{(L)}$ and $\mathbf{x}_{l+1}^{(L)}$ provided by BIF1, is conveyed by the message $m_{j}(\mathbf{z}_{l}^{(N)})$ (not appearing in Fig. 4). Note also that from eq. (5) (with $Z=L$) the equality

[TABLE]

is easily inferred; the pdf of $\mathbf{z}_{l}^{(N)}$ evaluated on the basis of the right-hand side (RHS) of eq. (15) is denoted $f(\mathbf{z}_{l}^{(N)}|\mathbf{x}_{l}^{(N)})$ in the following.

Message passing from BIF2 to BIF1 - BIF1 is fed by the message $m_{pm}(\mathbf{x}_{l})$ that, unlike the set $\{m_{pm,j}(\mathbf{x}_{l}^{(N)})\}$ passed to BIF2, provides pseudo-measurement information about the whole state $\mathbf{x}_{l}$. This message is generated as follows. The message set $\{m_{sm,j}(\mathbf{x}_{l}^{(N)})\}$ produced by the PF is processed in the PMG2→1 block, that computes the set of $N_{p}$ pseudo-measurement messages $\{m_{pm,j}(\mathbf{x}_{l}^{(L)})\}$ referring to the linear state component only. Then, the two sets $\{m_{pm,j}(\mathbf{x}_{l}^{(L)})\}$ and $\{m_{sm,j}(\mathbf{x}_{l}^{(N)})\}$ are merged in the PMC2→1 block, where the information they convey are converted into the (single) message $m_{pm}(\mathbf{x}_{l})$. Moreover, as illustrated in the Appendix, the message $m_{pm,j}(\mathbf{x}_{l}^{(L)})$ conveys a sample of the random vector [13]

[TABLE]

such a sample is generated under the assumption that $\mathbf{x}_{l}^{(N)}=\mathbf{x}_{fp,l,j}^{(N)}$. The pdf of the random vector $\mathbf{z}_{l}^{(L)}$ is evaluated on the basis of the joint statistical representation of the couple $(\mathbf{x}_{l}^{(N)}$, $\mathbf{x}_{l+1}^{(N)})$ produced by BIF2 and is conveyed by the message $m_{j}(\mathbf{z}_{l}^{(L)})$ (not appearing in Fig. 4); note also that from eq. (5) (with $Z=N$) the equality

[TABLE]

is easily inferred; the pdf of $\mathbf{z}_{l}^{(N)}$ evaluated on the basis of the RHS of eq. (17) is denoted $f(\mathbf{z}_{l}^{(L)}|\mathbf{x}_{l}^{(L)},\mathbf{x}_{l}^{(N)})$ in the following.

The rationale behind the message passing illustrated above can be summarized as follows. The message $m_{pm}(\mathbf{x}_{l})$ is extracted from the statistical information generated by BIF2 and is exploited by BIF1 to refine its backward estimate of the whole state; moreover, merging this estimate with the forward estimate $\vec{m}_{fe1}(\mathbf{x}_{l})$ allows to generate a more accurate statistical representation for $\mathbf{x}_{l}$ and, consequently, for $\mathbf{x}_{l}^{(L)}$ (these are conveyed by $m_{sm}(\mathbf{x}_{l})$ and $m_{sm}(\mathbf{x}_{l}^{(L)})$, respectively); finally, these statistical information are exploited to aid BIF2 in the computation of more refined weights of the particles representing $\mathbf{x}_{l}^{(N)}$.

Given the graphical model shown in Fig. 4 and the messages passed over it, the derivation of a specific BITF algorithm requires: a) defining the mathematical structure of the input messages that feed the $(T-l)$-th recursion of backward filtering and that of the output messages emerging from both backward filtering and smoothing in the same recursion; b) describing message scheduling; c) deriving mathematical expressions for all the computed messages. These issues are analysed in detail in Section 4.

4 Scheduling and Computation of Probabilistic Messages in Turbo

Smoothing Algorithms for CLG Models

In this Section, the specific issues raised at the end of the previous Section and concerning the message passing accomplished over the graphical model shown in Fig. 4 are addressed. For this reason, we first provide various details about a) the messages feeding backward filtering, and b) the messages emerging from it and from the related smoothing. Then, we focus on the scheduling of such messages and on their computation. This allows us to develop two new smoothing techniques, one solving problem P.1, the other one problem P.2. Finally, these techniques are briefly compared with other particle smoothing methods available in the literature.

4.1 Input and Output Messages

The input messages feeding the $(T-l)$-th recursion of backward filtering are generated in the $l$-th recursion of the paired forward filtering and in the previous recursion (i.e., in the $(T-l+1)$-th recursion) of the backward pass. In the following, various details about such messages are provided.

1. Input messages evaluated in the forward pass - A turbo filter, consisting of an EKF (denoted F1) and a PF (denoted F2), is employed in the forward pass of the devised TS algorithms and is run only once. Therefore, the forward predictions/estimates, provided by F1 (F2) and made available to BIF1 (BIF2), are expressed by Gaussian pdfs (sets of weighted particles), each conveyed by a Gaussian message (by a set of particle-dependent messages). The notation adopted in the following for these probabilistic information is summarized below.

Filter F1 - This filter, in its $(l-1)$-th recursion, computes the forward prediction of $\mathbf{x}_{l}$, conveyed by the message222Considerations similar to the ones expressed for $\vec{m}_{fp}(\mathbf{x}_{l})$ (18) and $\vec{m}_{fe1}(\mathbf{x}_{l})$ (19) can be repeated for the messages $\vec{m}_{fp,j}(\mathbf{x}_{l}^{(N)})$ (22) and $\vec{m}_{fe,j}(\mathbf{x}_{l}^{(N)})$ (23), respectively, defined below. (see Fig. 4)

[TABLE]

This message is updated in the $l$-th recursion of F1 on the basis of the measurement $\mathbf{y}_{l}$. This produces the Gaussian message

[TABLE]

representing a forward estimate of $\mathbf{x}_{l}$; the covariance matrix $\mathbf{C}_{fe1,l}$ and the mean vector $\mathbf{\eta}_{fe1,l}$ can be evaluated on the basis of the associated precision matrix (see [11, eqs. (14)-(17)])

[TABLE]

and of the transformed mean vector

[TABLE]

respectively; here, $\mathbf{W}_{e}\triangleq\mathbf{C}_{e}^{-1}$, $\mathbf{W}_{fp,l}\triangleq(\mathbf{C}_{fp,l})^{-1}$ and $\mathbf{w}_{fp,l}\triangleq\mathbf{W}_{fp,l}\mathbf{\eta}_{fp,l}$. The message $\vec{m}_{fp}(\mathbf{x}_{l})$ (18) enters the graphical model developed for BIF1 (see Fig. 4) along the half edge referring to $\mathbf{x}_{l}$.

Filter F2 - This filter, in its $(l-1)$-th recursion, computes the particle set $S_{fp,l}\triangleq\{\mathbf{x}_{fp,l,j}^{(N)},j=0,1,...,N_{p}-1\}$, representing a forward prediction of $\mathbf{x}_{l}^{(N)}$; the weight assigned to the particle $\mathbf{x}_{fp,l,j}^{(N)}$ is equal to $1/N_{p}$ for any $j$, since the use of particle resampling in each recursion is assumed. The statistical information available about $\mathbf{x}_{fp,l,j}^{(N)}$ are conveyed by the message

[TABLE]

with $j=0,1,...,N_{p}-1$. The weight of $\mathbf{x}_{fp,l,j}^{(N)}$ (with $j=0,1,...,N_{p}-1$) is updated by F2 in its $l$-th recursion on the basis of the measurement $\mathbf{y}_{l}$;  the new weight is denoted $w_{fe,l,j}$ and is conveyed by the forward message

[TABLE]

Note that the message set $\{\vec{m}_{fe1,j}(\mathbf{x}_{l}^{(N)})\}$ represents the forward estimate of $\mathbf{x}_{l}^{(N)}$ computed by F2 in its $l$-th recursion and that the message set $\{\vec{m}_{fp,j}(\mathbf{x}_{l}^{(N)})\}$ (see eq. (22)) enters the graphical model developed for BIF2 along the half edge referring to $\mathbf{x}_{l}^{(N)}$ (see Fig. 4).

2. Input messages evaluated in the backward pass - The $(T-l)$-th recursion of backward filtering is fed by the input messages

[TABLE]

and

[TABLE]

that convey the pdf of the backward estimate of $\mathbf{x}_{l+1}$ computed by BIF1 and the backward estimate of $\mathbf{x}_{l+1}^{(N)}$ generated by BIF2, respectively, in the previous recursion.

All the input messages described above are processed to compute: 1) the new backward estimates $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(N)})$, that represent the outputs emerging from the $(T-l)$-th recursion of backward filtering; 2) the required smoothed information (in the form of probabilistic messages) by merging forward and backward messages. In the remaining part of this Paragraph, some essential information about the structure of such messages are provided; details about their computation are given in the next Paragraph.

1. Computation of backward estimates - The computation of the message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l})$ (BIF1) and of the message set $\{\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l}^{(N)})\}$ (BIF2) is accomplished as follows. First, the backward prediction

[TABLE]

of $\mathbf{x}_{l}$ and the message

[TABLE]

conveying a backward weight for the $j$-th particle $\mathbf{x}_{fp,l,j}^{(N)}$ representing $\mathbf{x}_{l}^{(N)}$ (with $j=0,1,...,N_{p}-1$) are computed by BIF1 and BIF2, respectively. Then, in BIF1, the message $\overset{\leftarrow}{m}_{bp}\left(\mathbf{x}_{l}\right)$ (26) is merged with the pseudo-measurement message $m_{pm}(\mathbf{x}_{l})$ and the measurement message $m_{ms}(\mathbf{x}_{l})$ in order to compute

[TABLE]

and (see eq. (24))

[TABLE]

respectively. Similarly, in BIF2, the message $\overset{\leftarrow}{m}_{bp,j}(\mathbf{x}_{l}^{(N)})$ (27) is merged first with the pseudo-measurement message $m_{pm,j}(\mathbf{x}_{l}^{(N)})$ in order to produce the message (see eq. (25))

[TABLE]

conveying a new weight for the $j$-th particle $\mathbf{x}_{fp,l,j}^{(N)}$. Then, the information conveyed by the message set $\{\overset{\leftarrow}{m}_{be1,j}(\mathbf{x}_{l}^{(N)})\}$ is merged with that provided by the measurement-based set $\{m_{ms,j}(\mathbf{x}_{l}^{(N)})\}$ in order to evaluate the message (see eq. (25))

[TABLE]

that conveys a (particle-independent) backward estimate of $\mathbf{x}_{l}^{(N)}$.

2. Computation of smoothed information - In our work, the evaluation of smoothed information is based on the same conceptual approach as [11], [6] and [10]. In fact, the proposed method is based on the following ideas:

a) The joint smoothing pdf $f(\mathbf{x}_{1:T}|\mathbf{y}_{1:T})$ is estimated by providing multiple (say, $M$) realizations of it and a single realization (i.e., a single smoothed state trajectory) is computed in each backward pass; consequently, generating the smoothing output requires running a single forward pass and $M$ distinct backward passes.

b) The factorisation (12) is exploited to evaluate smoothed information, i.e. to merge the statistical information emerging from the forward pass with that computed in any of the $M$ backward passes. In particular, this formula is employed to combine the statistical information made available by F1 (F2) with those generated by BIF1 (BIF2); consequently, the first factor and the second one appearing in the RHS of eq. (12) are expressed by the forward message $\vec{m}_{fe1}(\mathbf{x}_{l})$ (19) and the backward message $\overset{\leftarrow}{m}_{bp}\left(\mathbf{x}_{l}\right)$ (26) (the forward message $\vec{m}_{fe1,j}(\mathbf{x}_{l}^{(N)})$ (23) and the backward message $\overset{\leftarrow}{m}_{bp,j}(\mathbf{x}_{l}^{(N)})$ (27)), respectively, if F1 and BIF1 (F2 and BIF2) are considered.

4.2 Scheduling and Computation of Probabilistic Messages

The message passing algorithm we propose for backward filtering and smoothing is iterative, since, within each recursion of the backward pass, it can accomplish multiple passes over the same edges. Moreover, it results from: a) the adoption of the message scheduling illustrated in Fig. 5, that refers to the $k$-th iteration of the devised algorithm; b) the use of the SPA in the evaluation of all the passed messages. It is also important to mention that the selected scheduling mimics the one employed in [11], which, in turn, has been inspired by [6] and [10]. Based on this scheduling, the computation of the messages passed over the given graphical model can be divided in the three consecutive phases listed below.

I - In this phase, $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1})$ ($BE^{\prime}$) is processed to compute $\overset{\leftarrow}{m}_{bp}(\mathbf{x}_{l})$ ($BP$) and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(L)})$ ($BEL^{\prime}$); moreover, the set $\{m_{pm,j}(\mathbf{x}_{l}^{(L)})\}$ ($PML_{j}$) conveying pseudo-measurement information about $\mathbf{x}_{l}^{(L)}$ is evaluated.

II - In the second phase, an iterative evaluation of the backward estimates of the whole state (BIF1) and of the nonlinear state component (BIF2) is accomplished. More specifically, in the $k$-th iteration of this procedure (with $k=1,2,...,N_{it}$, where $N_{it}$ is the overall number of iterations) the ordered computation of the following messages or sets of $N_{p}$ messages is accomplished in five consecutive steps333Note that the superscript $(k)$ ($(k-1)$) indicates that the associated message is computed in the $k$-th ($(k-1)$-th) iteration of phase II. (see Fig. 5): 1) $\{m_{sm,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$ ($SMN_{j}^{(k)}$), $m_{pm}^{(k)}(\mathbf{x}_{l})$ ($PM^{(k)}$); 2) $\overset{\leftarrow}{m}_{be1}^{(k)}(\mathbf{x}_{l})$ ($BE1^{(k)}$), $m_{sm}^{(k)}(\mathbf{x}_{l})$ ($SM^{(k)}$), $m_{sm}^{(k)}(\mathbf{x}_{l}^{(L)})$ ($SML^{(k)}$); 3) $\{m_{pm,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$ ($PMN_{j}^{(k)}$); 4) $\{\overset{\leftarrow}{m}_{bp,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$ ($BPN_{j}^{(k)}$), $\{\overset{\leftarrow}{m}_{be1,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$ ($BE1N_{j}^{(k)}$); 5) $\{m_{ms,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$ ($MSN_{j}^{(k)}$).

III - In the third phase, the smoothed information $\{m_{sm,j}^{(N_{it}+1)}(\mathbf{x}_{l}^{(N)})\}$ is computed and employed in the evaluation of: a) the output message $m_{be}(\mathbf{x}_{l}^{(N)})$ of BIF1; b) the new pseudo-measurement message $m_{pm}^{(N_{it}+1)}(\mathbf{x}_{l})$. Finally, $m_{pm}^{(N_{it}+1)}(\mathbf{x}_{l})$ is processed to compute $\overset{\leftarrow}{m}_{be1,l}^{(N_{it}+1)}\left(\mathbf{x}_{l}\right)$ and the output message $\overset{\leftarrow}{m}_{be,l}\left(\mathbf{x}_{l}\right)=$ $\overset{\leftarrow}{m}_{be2,l}\left(\mathbf{x}_{l}\right)$ of BIF2.

In the remaining part of this Section, the expressions of all the messages computed in each of the three phases described above are provided; the derivation of these expressions is sketched in the Appendix.

**Phase I **- The message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(L)})$ is computed as

[TABLE]

since it results from the marginalization of $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1})$ (24) with respect to $\mathbf{x}_{l+1}^{(N)}$; in practice, the mean vector $\mathbf{\tilde{\eta}}_{be,l+1}$ and the covariance matrix $\mathbf{\tilde{C}}_{be,l+1}$ are extracted from the parameters $\mathbf{\eta}_{be,l+1}$ and $\mathbf{C}_{be,l+1}$, respectively (since $\mathbf{x}_{l+1}^{(L)}$ consists of the first $D_{L}$ elements of $\mathbf{x}_{l+1}$).

The message $\overset{\leftarrow}{m}_{bp}(\mathbf{x}_{l})$ (26), representing a one-step backward prediction of $\mathbf{x}_{l}$, is computed on the basis of $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1})$ and the pdf $f(\mathbf{x}_{l+1}|\mathbf{x}_{l})$. Its parameters $\mathbf{\eta}_{bp,l}$ and $\mathbf{C}_{bp,l}$ are evaluated on the basis of the precision matrix

[TABLE]

and of the transformed mean vector

[TABLE]

respectively; here, $\mathbf{W}_{be,l+1}\triangleq(\mathbf{C}_{be,l+1})^{-1}$, $\mathbf{P}_{l+1}\triangleq\mathbf{\mathbf{I}}_{D}-\mathbf{W}_{be,l+1}\mathbf{Q}_{l+1}$, $\mathbf{Q}_{l+1}\triangleq(\mathbf{W}_{w}+\mathbf{W}_{be,l+1})^{-1}$, $\mathbf{W}_{w}\triangleq(\mathbf{C}_{w})^{-1}$ and $\mathbf{w}_{be,l+1}\triangleq\mathbf{W}_{be,l+1}\mathbf{\eta}_{be,l+1}$.

The evaluation of the set of messages $\{m_{pm,j}(\mathbf{x}_{l}^{(L)})\}$ is based on the message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(N)})$ (25) and on the particle set conveyed by the messages $\{m_{sm,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$ (such a set, being equal to $S_{fp,l}$, is independent of the iteration index $k$; see eq. (40)). In the Appendix it is shown that

[TABLE]

the covariance matrix $\mathbf{\tilde{C}}_{pm,l,j}$ and the mean vector $\mathbf{\tilde{\eta}}_{pm,l,j}$ are computed on the basis of the precision matrix

[TABLE]

and of the transformed mean vector

[TABLE]

respectively; here, $\mathbf{A}_{l,j}^{(N)}\triangleq\mathbf{A}_{l}^{(N)}(\mathbf{x}_{fp,l,j}^{(N)})$,

[TABLE]

is an iteration-independent pseudo-measurement and $\mathbf{f}_{l,j}^{(N)}\triangleq\mathbf{f}_{l}^{(N)}(\mathbf{x}_{fp,l,j}^{(N)})$.

Phase II - A short description of the five steps accomplished in the $k$-th iteration of this phase is provided in the following.

Step 1) Computation of the pseudo-measurements for BIF1- The message $m_{sm,j}^{(k)}(\mathbf{x}_{l}^{(N)})$ is evaluated as444Note that the messages $\overset{\rightarrow}{m}_{fe1,j}^{(k-1)}(\mathbf{x}_{l}^{(N)})\,$ and $\overset{\leftarrow}{m}_{be1,j}^{(k-1)}(\mathbf{x}_{l}^{(N)})$ appearing in the following formula are evaluated in the previous iteration and stored in the delay elements (identified by the letter D in Fig. 5). (see Fig. 5, and eqs. (23) and (30))

[TABLE]

where

[TABLE]

with $w_{fe1,l,j}^{(0)}=w_{fe,l,j}$ (see eq. (23)) and $w_{be1,l,j}^{(0)}=1$ (i.e., $w_{sm,l,j}^{(1)}=w_{fe,l,j}$). Then, the weights $\{w_{sm,l,j}^{(k)}\}$ are normalized; this produces the $j$-th normalised weight

[TABLE]

with $j=0,1,...,N_{p}-1$, where $K_{sm,l}^{(k)}\triangleq 1/\sum\limits_{j=0}^{N_{p}-1}w_{sm,l,j}^{(k)}$. Note that the particles $\{\mathbf{x}_{fp,l,j}^{(N)}\}$ and their new weights $\{W_{sm,l,j}^{(k)}\}$ provide a statistical representation of the smoothed estimate of $\mathbf{x}_{l}^{(N)}$ evaluated in the $k$-th iteration.

Then, the message

[TABLE]

is computed in the block PMC2→1 on the basis of the message sets $\{m_{pm,j}(\mathbf{x}_{l}^{(L)})\}$ (see eq. (35)) and $\{m_{sm,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$; the mean vector $\mathbf{\eta}_{pm,l}^{(k)}$ and the covariance matrix $\mathbf{C}_{pm,l}^{(k)}$ are evaluated as

[TABLE]

and

[TABLE]

respectively, where

[TABLE]

is a $D_{X}$-dimensional mean vector (with $X=L$ and $N)$,

[TABLE]

is a $D_{X}\times D_{Y}$ covariance (or cross-covariance) matrix (with $XY=LL$, $NN$ and $LN)$, $\mathbf{\eta}_{pm,l,j}^{(L)}=\mathbf{\tilde{\eta}}_{pm,l,j}$, $\mathbf{\eta}_{pm,l,j}^{(N)}=\mathbf{x}_{fp,l,j}^{(N)}$, $\mathbf{r}_{pm,l,j}^{(LL)}\triangleq\mathbf{\tilde{C}}_{pm,l,j}+\mathbf{\tilde{\eta}}_{pm,l,j}(\mathbf{\tilde{\eta}}_{pm,l,j})^{T}$, $\mathbf{r}_{pm,l,j}^{(NN)}\triangleq\mathbf{x}_{fp,l,j}^{(N)}(\mathbf{x}_{fp,l,j}^{(N)})^{T}$ and $\mathbf{r}_{pm,l,j}^{(LN)}\triangleq\mathbf{\tilde{\eta}}_{pm,l,j}(\mathbf{x}_{fp,l,j}^{(N)})^{T}$.

Step 2) Computation of the backward and smoothed estimates in BIF1 - The message $\overset{\leftarrow}{m}_{be1}^{(k)}(\mathbf{x}_{l})$ is evaluated as (see Fig. 5)

[TABLE]

where the messages $\overset{\leftarrow}{m}_{bp}(\mathbf{x}_{l})$ and $m_{pm}^{(k)}(\mathbf{x}_{l})$ are given by eq. (26) and eq. (43), respectively. The covariance matrix $\mathbf{C}_{be1,l}^{(k)}$ and the mean vector $\mathbf{\eta}_{be1,l}^{(k)}$ are computed on the basis of the associated precision matrix

[TABLE]

and transformed mean vector

[TABLE]

respectively; here, $\mathbf{W}_{pm,l}^{(k)}\triangleq(\mathbf{C}_{pm,l}^{(k)})^{-1}$, $\mathbf{w}_{pm,l}^{(k)}\triangleq\mathbf{W}_{pm,l}^{(k)}\,\mathbf{\eta}_{pm,l}^{(k)}$, and $\mathbf{W}_{bp,l}$ and $\mathbf{w}_{bp,l}$ are given by eqs. (33) and (34), respectively. From eqs. (50)-(51) the expressions

[TABLE]

and

[TABLE]

can be easily inferred; here, $\mathbf{W}_{l}^{(k)}\triangleq[\mathbf{C}_{pm,l}^{(k)}\mathbf{W}_{bp,l}+\mathbf{I}_{D}]^{-1}$.

Then, the message $m_{sm}^{(k)}(\mathbf{x}_{l})$ is evaluated as (see Fig. 5)

[TABLE]

where the messages $\vec{m}_{fe1}\left(\mathbf{x}_{l}\right)$ and $\overset{\leftarrow}{m}_{be1}^{(k)}\left(\mathbf{x}_{l}\right)$ are given by eqs. (19) and (49), respectively. The covariance matrix $\mathbf{C}_{sm,l}^{(k)}$ and the mean vector $\mathbf{\eta}_{be1,l}^{(k)}$ are computed on the basis of the associated precision matrix

[TABLE]

and transformed mean vector

[TABLE]

respectively. Finally, marginalizing $m_{sm}^{(k)}(\mathbf{x}_{l})$ (55) with respect to $\mathbf{x}_{l}^{(N)}$ results in the message

[TABLE]

where $\mathbf{\tilde{\eta}}_{sm,l}^{(k)}$ and $\mathbf{\tilde{C}}_{sm,l}^{(k)}$ are extracted from the mean $\mathbf{\eta}_{sm,l}^{(k)}$ and the covariance matrix $\mathbf{C}_{sm,l}^{(k)}$ of $m_{sm}^{(k)}(\mathbf{x}_{l})$ (55), respectively (since $\mathbf{x}_{l}^{(L)}$ consists of the first $D_{L}$ elements of $\mathbf{x}_{l}$).

Step 3) Computation of the pseudo-measurements for BIF2** **- The pseudo-measurement information feeding BIF2 is conveyed by the message set $\{m_{pm,j}^{(k)}(\mathbf{x}_{l}^{(N)})\triangleq w_{pm,l,j}^{(k)}\}$, i.e. by a set of new weights for the particles forming the set $S_{fp,l}$. The $j$-th weight is evaluated as

[TABLE]

for any $j$; here,

[TABLE]

$\left\|\mathbf{x}\right\|_{\mathbf{W}}^{2}\triangleq\mathbf{x}^{T}\mathbf{Wx}$ denotes the square of the norm of the vector $\mathbf{x}$ with respect to the positive definite matrix $\mathbf{W}$,

[TABLE]

[TABLE]

$\mathbf{\check{W}}_{z,l,j}^{(k)}\triangleq(\mathbf{\check{C}}_{z,l,j}^{(k)})^{-1}$, $\mathbf{\check{w}}_{z,l,j}^{(k)}\triangleq\mathbf{\check{W}}_{z,l,j}^{(k)}\mathbf{\check{\eta}}_{z,l,j}^{(k)}$, $\mathbf{\check{\eta}}_{z,l,j}^{(k)}$ and $\mathbf{\check{C}}_{z,l,j}^{(k)}$ are expressed by eqs. (97) and (98), respectively, $\mathbf{W}_{w}^{(L)}\triangleq[\mathbf{C}_{w}^{(L)}]^{-1}$, $\mathbf{f}_{l,j}^{(L)}\triangleq\mathbf{f}_{l}^{(L)}(\mathbf{x}_{fp,l,j}^{(N)})$, $D_{pm,l,j}^{(k)}\triangleq[\det(\mathbf{\check{C}}_{l,j}^{(k)})]^{-D_{L}/2}$ and $\mathbf{\check{C}}_{l,j}^{(k)}\triangleq\mathbf{\check{C}}_{z,l,j}^{(k)}+\mathbf{C}_{w}^{(L)}$.

Step 4) Computation of the backward weights in BIF2 - The backward message $\overset{\leftarrow}{m}_{bp,j}^{(k)}(\mathbf{x}_{l}^{(N)})$ (27), i.e. the backward weight (see Fig. 5) is computed as

[TABLE]

where $D_{bp,l,j}^{(k)}=(2\pi\det(\mathbf{C}_{1,l,j}^{(N)}))^{-D_{N}/2}$,

[TABLE]

[TABLE]

$\mathbf{W}_{1,l,j}^{(N)}[k]\triangleq(\mathbf{C}_{1,l,j}^{(N)}[k])^{-1}$ and

[TABLE]

Then, the backward message $\overset{\leftarrow}{m}_{be1,j}^{(k)}(\mathbf{x}_{l}^{(N)})$ is evaluated as (see Fig. 5)

[TABLE]

Based on eqs. (59) and (64), the last formula can be rewritten as

[TABLE]

where

[TABLE]

for any $j$, where $D_{be1,l,j}^{(k)}\triangleq D_{pm,l,j}^{(k)}\,D_{bp,l,j}^{(k)}$ and $Z_{be1,l,j}^{(k)}\triangleq Z_{pm,l,j}^{(k)}+Z_{bp,l,j}^{(k)}$. The messages $\{\overset{\leftarrow}{m}_{be1,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$ (i.e., the weights $\{w_{be1,l,j}^{(k)}\}$) are stored for the next iteration (see step1)).

Step 5) Computation of new measurement-based weights in BIF2- The new measurement-based weight (see Fig. 5)

[TABLE]

is computed on the basis of $m_{sm}^{(k)}(\mathbf{x}_{l}^{(L)})$ (58); here,

[TABLE]

and

[TABLE]

where $\mathbf{B}_{l,j}\triangleq\mathbf{B}_{l}(\mathbf{x}_{fp,l,j}^{(N)})$ and $\mathbf{g}_{l,j}\triangleq\mathbf{g}_{l}(\mathbf{x}_{fp,l,j}^{(N)})$. Then, the $N_{p}$ messages $\{m_{ms,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$ (i.e., the weights $\{w_{fe1,l,j}^{(k)}\,\}$) are stored, since in the next iteration they are employed to generate the message (see Fig. 5, and eqs. (22) and (LABEL:eq:weight_before_resampling))

[TABLE]

and, then, the message $m_{sm,j}^{(k)}(\mathbf{x}_{l}^{(N)})$ (39) (i.e., the smoothed weight $w_{sm,l,j}^{(k)}$ (41)); this concludes the $k$-th iteration. Then, the index $k$ is increased by one, and a new iteration is started by going back to step 1) if $k<N_{it}+1$; otherwise (i.e., if $k=N_{it}+1$, we proceed with the next phase.

Phase III - In this phase, only step 1) and part of step 2) of phase II are carried out in order to compute all the statistical information required for the evaluation of the backward estimates $\overset{\leftarrow}{m}_{be}\left(\mathbf{x}_{l}\right)$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(N)})$, i.e. the outputs generated by BIF1 and BIF2, respectively, in the $l$-th recursion of TS. More specifically, the smoothed information $\{m_{sm,j}^{(N_{it}+1)}(\mathbf{x}_{l}^{(N)})\}$ is computed (as if an additional iteration was started; see eqs. (40)-(41)), the new weights $\{W_{sm,l,j}^{(N_{it}+1)}\}$ are evaluated on the basis of eq. (42) and the set $S_{fp,l}$ is sampled once on the basis of such weights; if the $j_{l}$-th particle (i.e., $\mathbf{x}_{fp,l,j_{l}}^{(N)}$) is selected, we set

[TABLE]

so that the message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(N)})$ (25) becomes available at the output of BIF1. On the other hand, the evaluation of the message $\overset{\leftarrow}{m}_{be}\left(\mathbf{x}_{l}\right)$ is accomplished as follows. The messages $m_{pm}^{(N_{it}+1)}(\mathbf{x}_{l})$ and $\overset{\leftarrow}{m}_{be1,l}^{(N_{it}+1)}\left(\mathbf{x}_{l}\right)$ are computed first (see eq. (43) and eqs. (48)-(49), respectively). Then, the BIF2 output message $\overset{\leftarrow}{m}_{be,l}\left(\mathbf{x}_{l}\right)$ is computed as (see Fig. 5)

[TABLE]

where

[TABLE]

is the message conveying the measurement information. Moreover, the covariance matrices $\mathbf{C}_{ms,l}$ and $\mathbf{C}_{be2,l}$, and the mean vectors $\mathbf{\eta}_{ms,l}$ and $\mathbf{\eta}_{be2,l}$ are computed on the basis of the associated precision matrices

[TABLE]

[TABLE]

and of the transformed mean vectors

[TABLE]

[TABLE]

respectively. The $l$-th recursion is now over.

It is important to point out that the first recursion of the backward pass requires the knowledge of the input messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{T})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{T}^{(N)})$. Similarly as any BIF algorithm, the evaluation of these messages in BITF is based on the statistical information generated in the last recursion of the forward pass. In particular, the above mentioned messages are still expressed by eqs. (29) and (31) (with $l=T$ in both formulas), respectively. However, the vector $\mathbf{x}_{be,T}^{(N)}$ is generated by sampling the particle set $S_{fp,T}$ on the basis of the forward weights $\{w_{fe,T,j}\}$, since backward predictions are unavailable at the final instant $l=T$. Therefore, if the $j_{T}$-th particle of $S_{fp,T}$ is selected, we set

[TABLE]

in the message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{T}^{(N)})$ entering the BIF2 in the first recursion (see eq. (25)). As far as BIF1 is concerned, following [11], we choose

[TABLE]

and

[TABLE]

for the message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{T})$.

The general method for BITF and TS developed in this Paragraph is summarized in Algorithm 1.

Algorithm 1 produces all the statistical information required to solve problems P.1 and P.2. Let us now discuss how this can be done in detail. As far as problem** P.1** is concerned, it is useful to point out that Algorithm 1 produces a trajectory $\{\mathbf{x}_{be,l}^{(N)},l=1,2,...,T\}$ for the nonlinear component (see eq. (76)). Another trajectory, representing the time evolution of the linear state component only and denoted $\{\mathbf{x}_{be,l}^{(L)},l=1,2,...,T\}$, can be computed by sampling the message $\vec{m}_{sm}^{(N_{it}+1)}(\mathbf{x}_{l}^{(L)})$ (see eq. (58)) or by simply setting $\mathbf{x}_{be,l}^{(L)}=\mathbf{\tilde{\eta}}_{sm,l}^{(N_{it}+1)}$ (this task can be accomplished in task in step 3-h of Algorithm 1, after sampling the particle set $S_{fp,l}$). The overall algorithm producing this result is called *turbo smoothing algorithm *(TSA) in the following.

The TSA solves problem P.1 and, consequently, problem P.2, since, once it has been run, an approximation of the marginal smoothed pdf at any instant can be simply obtained by marginalization. The last result, however, is achieved at the price of a significant computational cost since $M$ backward passes are required. However, if we are interested in solving problem P.2 only, a simpler particle smoother can be developed following the approach illustrated in [11], so that a single backward pass has to be run. In this pass, the evaluation of the message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(N)})$ (i.e., of the particle $\mathbf{x}_{be,l}^{(N)}$) involves the whole particle set $S_{fp,l}$ and their weights $\{W_{sm,l,j}^{(N_{it}+1)}\}$ (see eq. (42)) evaluated in the last phase of the $l$-th recursion. More specifically, a new smoother is obtained by employing a different method for evaluating $\mathbf{x}_{be,l}^{(N)}$ in step 3-h of Algorithm 1; it consists in computing the smoothed estimate

[TABLE]

of $\mathbf{x}_{l}^{(N)}$ and, then, setting

[TABLE]

The resulting smoother is called simplified turbo smoothing algorithm (STSA) in the following.

Finally, it is important to point out that the computational complexity of the TSA and the STSA can be substantially reduced by reusing the forward weights $\{w_{fe1,l,j}\}$ in all the iterations of phase II, so that step 5) can be skipped; this means that, for any $k$, we set $w_{fe1,l,j}^{(k-1)}=w_{fe1,l,j}$ in the evaluation of the $j$-th particle weight $w_{sm,l,j}^{(k)}$ according to eq. (41) in step 1). Our simulation results have evidenced that, at least for the SSM considered in Section 5, this modification does not have any impact on the estimation accuracy of these algorithms.

4.3 Comparison of the Developed Turbo Smoothing Algorithms with

Related Techniques

The TSA developed in the previous Section is conceptually related to the Rao-Blackwellized particle smoothing (RBPS) techniques proposed by Fong *et al. *[6] and by Lindsten et al. [10] (these algorithms are denoted Alg-B and Alg-L respectively, in the following) and to the RBSS algorithm devised by Vitetta et al. [11]. In fact, all these techniques share with the TSA the following important features: 1) all of them aim at estimating the joint smoothing density over the whole observation interval by generating multiple realizations from it; 2) they accomplish a single forward pass and as many backward passes as the overall number of realizations; 3) they combine Kalman filtering with particle filtering. However, Alg-B, Alg-L and the RBSS algorithm employ, in both their forward and backward passes, as many Kalman filters as the number of particles ($N_{p}$) to generate a particle-dependent estimate of the linear state component only. On the contrary, the TSA employs a single (extended) Kalman filter, that, however, estimates the whole system state. This substantially reduces the memory requirements of particle smooothing and, consequently, the overall number of memory accesses accomplished on the hardware platform on smoothing is run; as evidenced by our numerical results, this feature contributes to making the overall execution time of TSA appreciably shorter than that required by the related algorithms.

On the other hand, the STSA is conceptually related to the SPS algorithm devised by Vitetta et al. [11]. In fact, both algorithms aim at solving problem P.2 only and, consequently, carry out a single backward pass. This property makes them much faster than Alg-B, Alg-L and the RBSS algorithm in the computation of marginal smoothed densities. Finally, note that, similarly as the TS technique, the use of the STSA requires a substantially smaller number of memory accesses than the SPS algorithm.

5 Numerical Results

In this Section we compare, in terms of accuracy and execution time, the TSA and the STSA with Alg-L, the RBSS and the SPS algorithm for a specific CLG SSM. The considered SSM is the same as the SSM#2 defined in [11] and describes the bidimensional motion of an agent. Its state vector in the $l$-th observation interval is defined as $\mathbf{x}_{l}\triangleq[\mathbf{v}_{l}^{T},\mathbf{p}_{l}^{T}]^{T}$, where $\mathbf{v}_{l}\triangleq[v_{x,l},v_{y,l}]^{T}$ and $\mathbf{p}_{l}\triangleq[p_{x,l},p_{y,l}]^{T}$ (corresponding to $\mathbf{x}_{l}^{(L)}$ and $\mathbf{x}_{l}^{(N)}$, respectively) represent the agent velocity and position, respectively (their components are expressed in m/s and in m, respectively). The state update equations are

[TABLE]

and

[TABLE]

where $\rho$ is a forgetting factor (with $0<\rho<1$), $T_{s}$ is the sampling interval, $\mathbf{n}_{v,l}$ is an additive Gaussian noise (AGN) vector characterized by the covariance matrix $\mathbf{I}_{2}$,

[TABLE]

is the acceleration due to a force applied to the agent (and pointing towards the origin of our reference system), $a_{0}$ is a scale factor (expressed in m/s2), $d_{0}$ is a reference distance (expressed in m), and $\mathbf{n}_{p,l}$ is an AGN vector characterized by the covariance matrix $\sigma_{p}^{2}\mathbf{I}_{2}$ and accounting for model inaccuracy. The measurement vector available in the $l$-th interval for state estimation is

[TABLE]

where $\mathbf{e}_{l}\triangleq[\mathbf{e}_{v,l}^{T},\mathbf{e}_{p,l}^{T}]^{T}$ and $\mathbf{e}_{v,l}$ ($\mathbf{e}_{p,l}$) is an AGN vector characterized by the covariance matrix $\sigma_{ev}^{2}\mathbf{I}_{2}$ ($\sigma_{ep}^{2}\mathbf{I}_{2}$).

In our computer simulations, following [11] and [12], the estimation accuracy of the considered smoothing techniques has been assessed by evaluating two root mean square errors (RMSEs), one for the linear state component, the other for the nonlinear one, over an observation interval lasting $T=200$ $T_{s}$; these are denoted $RMSE_{L}($alg$)$ and $RMSE_{N}($alg$)$, respectively, where ‘alg’ is the acronym of the algorithm these parameters refer to. Our assessment of computational requirements is based, instead, on assessing the average computation time required for processing a single block of measurements (this quantity is denoted CTB$($alg$)$ in the following). Moreover, the following values have been selected for the parameters of the considered SSM: $\rho=0.995$, $T_{s}=0.01$ s, $\sigma_{p}$ $=5\cdot 10^{-3}$ m, $\sigma_{e,p}=2\cdot 10^{-2}$ m, $\sigma_{e,v}=2\cdot 10^{-2}$ m/s, $a_{0}=0.5$ m/s2, $d_{0}=5\cdot 10^{-3}$ m and $v_{0}=1$ m/s (the initial position $\mathbf{p}_{0}\triangleq[p_{x,0},p_{y,0}]^{T}$ and the initial velocity $\mathbf{v}_{0}\triangleq[v_{x,0},v_{y,0}]^{T}$ have been set to $[0.01$ m, $0.01$ m$]^{T}$ and $[0.01$ m/s, $0.01$ m/s$]^{T}$, respectively).

Some numerical results showing the dependence of $RMSE_{L}$ and $RMSE_{N}$ on the number of particles ($N_{p}$) for the considered smoothing algorithms are illustrated in Figs. 6 and 7, respectively (simulation results are indicated by markers, whereas continuous lines are drawn to fit them, so facilitating the interpretation of the available data). In this case, $N_{it}=1$ has been selected for both the TSA and the STSA, and the range $[10,150]$ has been considered for $N_{p}$ (since no real improvement is found for $N_{p}\gtrsim 150$). Morever, $RMSE_{L}$ and $RMSE_{N}$ results are also provided for MPF (TF with $N_{it}=1$), since this filtering technique is employed in the forward pass of Alg-L, the RBSS algorithm and the SPS algorithm (the TSA and the STSA); this allows us to assess the improvement in estimation accuracy provided by the backward pass with respect to the forward pass for each smoothing algorithm. These results show that:

1. The TSA, the STSA, Alg-L and the RBSS algorithm achieve similar accuracies in the estimation of both the linear and nonlinear state components.

2. The SPS algorithm is slightly outperformed by the other four smoothing algorithms in terms of $RMSE_{N}$ only; for instance, $RMSE_{N}($SPS$)$ is about $1.11$ times larger than $RMSE_{N}($STSA$)$ for $N_{p}=100$.

3. Even if the RBSS algorithm and the TSA provide by far richer statistical information than their simplified counterparts (i.e., than the SPS algorithm and the STSA, respectively), they do not provide a significant improvement in the accuracy of state estimation; for instance, $RMSE_{N}($SPS$)$ ($RMSE_{N}($STSA$)$) is about $1.12$ ($1.03$) time larger than $RMSE_{N}($RBSS$)$ ($RMSE_{N}($TSA$)$) for $N_{p}=100$.

4. The accuracy improvement in terms of $RMSE_{L}$ ($RMSE_{N}$) provided by all the smoothing algorithms except the SPS (Alg-L, RBSS, TSA and the STSA) is about $24\%$ (roughly $23\%$) with respect to the MPF and TF techniques, for $N_{p}=100$. Moreover, the accuracy improvement in terms of $RMSE_{L}$ ($RMSE_{N}$) achieved by the SPS algorithm is about $24\%$ (about $14\%$) with respect to the MPF technique for $N_{p}=100$.

Note also that, in the considered scenario, TF is slightly outperformed by (perform similarly as) MPF in the estimation of the linear (nonlinear) state component; a similar result is reported in [7] for a different SSM.

Despite their similar accuracies, the considered smoothing algorithms require different computational efforts; this is easily inferred from the numerical results appearing in Fig. 8 and illustrating the dependence of the CTB on $N_{p}$ for all the above mentioned filtering and smoothing algorithms. In fact, these results show that the TSA requires a shorter computation time than Alg-L and the RBSS algorithm; more specifically, CTB$($TSA$)$ is approximately $0.85$ ($0.48$) times smaller than CTB$($Alg-L$)$ (CTB$($RBSS$)$). The same considerations apply to the STSA and the SPS algorithm; in fact, CTB$($STSA$)$ is approximately $0.57$ times smaller than CTB$($SPS$)$. Note also that CTB$($TF$)$ is approximately $0.55$ times smaller than CTB$($MPF$)$ for the same value of $N_{p}$; once again, this result is in agreement with the results shown in [7] for a different SSM.

Finally, all the numerical results illustrated above lead to the conclusion that, in the considered scenario, the TSA and STSA achieve the best accuracy-complexity tradeoff in their categories of smoothing techniques.

6 Conclusions

In this manuscript, factor graph methods have been exploited to formalise the concept of parallel concatenation of Bayesian information filters. This has allowed us to develop a new approximate method for Bayesian smoothing, called turbo smoothing. Two turbo smoothers have been derived for the class of CLG systems and have been compared, in terms of both accuracy and execution time, with other smoothing algorithms for a specific dynamic model. These smoothers have limited requirements in terms of memory; moreover, our simulation results evidence that they perform similarly as their counterparts, but are faster.

Appendix

In this Appendix, the derivation of the expressions of various messages evaluated in each of the three phases the TFA consists of is sketched.

Phase I - Formulas (33) and (34), referring to the message $\overset{\leftarrow}{m}_{bp}(\mathbf{x}_{l})$ (26), can be easily computed by applying eqs. (IV.6)-(IV.8) of  [8, Table 4, p.1304] in their backward form (with $A\rightarrow\mathbf{I}_{D}$, $X\rightarrow\mathbf{F}_{l}\mathbf{x}_{l}$, $Z\rightarrow\mathbf{x}_{l+1}$ and $Y\rightarrow\mathbf{u}_{l}+\mathbf{w}_{l}$) and, then, eqs. (III.5)-(III.6) of [8, Table 3, p.1304] (with $A\rightarrow\mathbf{F}_{l}$, $X\rightarrow\mathbf{x}_{l}$ and $Y\rightarrow\mathbf{F}_{l}\mathbf{x}_{l}$).

The message set $\{m_{pm,j}(\mathbf{x}_{l}^{(L)})\}$ (see eq. (35)) conveys the statistical information provided by the pseudo-measurement $\mathbf{z}_{l}^{(L)}$ (16). The method for computing the message $m_{pm,j}(\mathbf{x}_{l}^{(L)})$ can be represented as a message passing over the graphical model shown in Fig. 9-a). Given $\mathbf{x}_{l}^{(N)}=\mathbf{x}_{fp,l,j}^{(N)}$ (this particle is provided by the message $m_{sm,j}^{(k)}(\mathbf{x}_{l}^{(N)})$ (40)) and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(N)})$ (25), the pseudo-measurement $\mathbf{z}_{l,j}^{(L)}$ (38) associated with the couple $(\mathbf{x}_{fp,l,j}^{(N)}$ , $\mathbf{x}_{be,l+1}^{(N)})$ is computed on the basis of eq. (16); this pseudo-measurement is conveyed by the message (denoted $ZL_{j}$ in Fig. 9-a))

[TABLE]

which is employed in the evaluation of the message (see Fig. 9-(a))

[TABLE]

Then, substituting eq. (93) and $f(\mathbf{z}_{l}^{(L)}|\mathbf{x}_{l}^{(L)},\mathbf{x}_{l}^{(N)})=\mathcal{N(}\mathbf{z}_{l}^{(L)}$; $\mathbf{A}_{l,j}^{(N)}\mathbf{x}_{l}^{(L)},\mathbf{C}_{w}^{(N)})$ (see eq. (17)) in the RHS of eq. (94) yields the message $m_{pm,j}^{(k)}(\mathbf{x}_{l}^{(L)})=\mathcal{\mathcal{N(}}\mathbf{z}_{l,j}^{(L)};\mathbf{A}_{l,j}^{(N)}\mathbf{x}_{l}^{(L)},\mathbf{C}_{w}^{(N)})$ (see [11, App. A, TABLE II, formula no. 3]), that can be easily put in the equivalent Gaussian form (35).

Phase II - Step 1) The procedure we adopt for computing $\overset{\leftarrow}{m}_{pm}^{(k)}(\mathbf{x}_{l})$ (43) on the basis of the sets $\{m_{pm,j}(\mathbf{x}_{l}^{(L)})\}$ and $\{m_{sm,j}^{(k)}(\mathbf{x}_{l}^{(N)})\}$ (see eqs. (35) and (40), respectively) is based on the following considerations. The message $m_{pm,j}(\mathbf{x}_{l}^{(L)})$ is coupled with $m_{sm,j}^{(k)}(\mathbf{x}_{l}^{(N)})$ (for any $j$), since they refer to the same particle set (i.e., $S_{fp,l}$). Moreover, these two messages provide complementary information, because they refer to the two different components of the overall state $\mathbf{x}_{l}$. For these reasons, the statistical information conveyed by the above mentioned sets can be condensed in the joint pdf

[TABLE]

Then, the message $m_{pm}^{(k)}(\mathbf{x}_{l})$ (43) can be computed by projecting this pdf* *onto a single Gaussian pdf; the transformation adopted here to achieve this result and expressed by eqs. (44)-(47) is described in [15, Sec. IV], and ensures that the mean and the covariance of the given pdf are preserved.

Step 2) The expression (49) of $\overset{\leftarrow}{m}_{be1}^{(k)}(\mathbf{x}_{l})$ represents a straightforward application of formula no. 2 of [12, App. A, TABLE I] (with $\mathbf{W}_{1}\rightarrow\mathbf{W}_{bp,l}$, $\mathbf{W}_{2}\rightarrow\mathbf{W}_{pm,l}^{(k)}$, $\mathbf{w}_{1}\rightarrow\mathbf{w}_{bp,l}$ and $\mathbf{w}_{2}\rightarrow\mathbf{w}_{pm,l}^{(k)}$). The same considerations apply to the derivation of the expression (55) of $m_{sm}^{(k)}(\mathbf{x}_{l})$.

Step 3) The algorithm for computing $m_{pm,j}^{(k)}(\mathbf{x}_{l}^{(N)})$ (59)  can be represented as a message passing over the graphical model shown in Fig. 9-b), in which the pseudo-measurement $\mathbf{z}_{l}^{(N)}$ (14) is computed. The expressions of the involved messages can be derived as follows. Given $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(L)})$ (32) and $m_{sm}^{(k)}(\mathbf{x}_{l}^{(L)})$ (58), the message $m_{j}^{(k)}(\mathbf{z}_{l}^{(N)})$ can expressed as (see [7, eqs. (83)-(84)])

[TABLE]

where

[TABLE]

[TABLE]

and $\mathbf{A}_{l,j}^{(L)}=\mathbf{A}_{l}^{(L)}(\mathbf{x}_{fp,l,j}^{(N)})$. Then, $\vec{m}_{j}^{(k)}(\mathbf{z}_{l}^{(N)})$ (96) is exploited to evaluate (see Fig. 9-b))

[TABLE]

Substituting eq. (96) and $f(\mathbf{z}_{l}^{(N)}|\mathbf{x}_{fp,l,j}^{(N)})=\mathcal{N}(\mathbf{z}_{l}^{(N)};\mathbf{f}_{l,j}^{(L)},\mathbf{C}_{w}^{(N)})$ (see eq. (15)) in the RHS of the last expression and evaluating the resulting integral (on the basis of formula no. 4 of [12, App. A, TABLE II]) yields eq. (59).

Step 4) The expression (64) of the weight $w_{bp,l,j}^{(k)}$ is derived as follows. First, we substitute $f(\mathbf{x}_{l+1}^{(N)}/\mathbf{x}_{l}^{(N)},\mathbf{x}_{l}^{(L)})=\mathcal{N(}\mathbf{x}_{l+1}^{(N)}$; $\mathbf{A}_{l}^{(N)}(\mathbf{x}_{l}^{(N)})\mathbf{x}_{l}^{(L)}+\mathbf{f}_{l}^{(N)}(\mathbf{x}_{l}^{(N)}),\mathbf{C}_{w}^{(N)})$ (see eq. (5) with $Z=N$), and the expressions of the messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(N)})$ (25) and $m_{sm}^{(k)}(\mathbf{x}_{l}^{(L)})$ (58) in the RHS of eq. (63). Then, the resulting integral is solved by applying formula no. 1 of [12, App. A, TABLE II] in the integration with respect to $\mathbf{x}_{l}^{(L)}$ and the sifting property of the Dirac delta function in the integration with respect to $\mathbf{x}_{l+1}^{(N)}$.

Step 5) The expression (LABEL:eq:weight_before_resampling) of the weight $w_{fe1,l,j}^{(k)}$ is derived as follows. First, we substitute $f(\mathbf{y}_{l}|\mathbf{x}_{fp,l,j}^{(N)},\,\mathbf{x}_{l}^{(L)})=\mathcal{N(}\mathbf{y}_{l}$; $\mathbf{g}_{l,j}+\mathbf{B}_{l,j}\mathbf{x}_{l}^{(L)},\mathbf{C}_{e})$ (with $\mathbf{B}_{l,j}\triangleq\mathbf{B}_{l}(\mathbf{x}_{fp,l,j}^{(N)})$ and $\mathbf{g}_{l,j}\triangleq\mathbf{g}_{l}(\mathbf{x}_{fp,l,j}^{(N)})$; see eq. (6)), and eq. (58) in eq. (4.2). Then, the resulting integral is solved by applying formula no. 1 of [12, App. A, TABLE II].

Phase III - The expression (78) of $\overset{\leftarrow}{m}_{be2,l}\left(\mathbf{x}_{l}\right)$ results from the application of formula no. 2 of [12, App. A, TABLE I] to eq. (77).