Rao-Blackwellized Particle Smoothing as Message Passing

TL;DR

This paper introduces a novel Rao-Blackwellized particle smoother for conditionally linear Gaussian state-space models, using a factor graph approach to improve fixed-lag smoothing accuracy and efficiency.

Contribution

It develops a new Rao-Blackwellized particle smoothing method tailored for conditionally linear Gaussian models, enhancing existing smoothing techniques with a factor graph framework.

Findings

Improved smoothing accuracy demonstrated on a specific model.

Reduced computational complexity compared to traditional methods.

Effective point mass approximation of the joint smoothing distribution.

Abstract

In this manuscript the fixed-lag smoothing problem for conditionally linear Gaussian state-space models is investigated from a factor graph perspective. More specifically, after formulating Bayesian smoothing for an arbitrary state-space model as forward-backward message passing over a factor graph, we focus on the above mentioned class of models and derive a novel Rao-Blackwellized particle smoother for it. Then, we show how our technique can be modified to estimate a point mass approximation of the so called joint smoothing distribution. Finally, the estimation accuracy and the computational requirements of our smoothing algorithms are analysed for a specific state-space model.

Tables

Table 1. Table 1: Mathematical rules for the evaluation of the message m o ​ u ​ t ​ ( 𝐱 ) subscript 𝑚 𝑜 𝑢 𝑡 𝐱 m_{out}(\mathbf{x}) , emerging from an equality node fed by the input messages m i ​ n , 1 ​ ( 𝐱 ) subscript 𝑚 𝑖 𝑛 1 𝐱 m_{in,1}(\mathbf{x}) and m i ​ n , 2 ​ ( 𝐱 ) subscript 𝑚 𝑖 𝑛 2 𝐱 m_{in,2}(\mathbf{x}) .

Formula no.	$m_{in,1}\left(𝐱\right)$	$m_{in,2}\left(𝐱\right)$	$m_{out}\left(𝐱\right)$
1	$\delta \left(𝐱-𝐚\right)$	$f\left(𝐱\right)$	$f\left(𝐚\right)\delta \left(𝐱-𝐚\right)$
2	$𝒩(𝐱;{\eta }_1,{𝐂}_1)$	$𝒩(𝐱;{\eta }_2,{𝐂}_2)$	$𝒩(𝐱;\eta ,𝐂),$ $𝐰={𝐰}_1+{𝐰}_2$, $𝐖={𝐖}_1+{𝐖}_2$
3	$𝒩(𝐱;{\eta }_1,{𝐂}_1)$	$𝒩(𝐜;\mathrm{𝐀𝐱}+𝐛,{𝐂}_2)$	$𝒩(𝐱;\eta ,𝐂),$ $𝐰={𝐰}_1+{𝐀}^T{𝐖}_2𝐜$, $𝐖={𝐖}_1+{𝐀}^T{𝐖}_2𝐀$

Table 2. Table 2: Mathematical rules for the evaluation of the message m o ​ u ​ t ​ ( 𝐱 2 ) subscript 𝑚 𝑜 𝑢 𝑡 subscript 𝐱 2 m_{out}(\mathbf{x}_{2}) , emerging from a function node f ​ ( 𝐱 1 , 𝐱 2 ) 𝑓 subscript 𝐱 1 subscript 𝐱 2 f(\mathbf{x}_{1},\mathbf{x}_{2}) on the basis of the input message m i ​ n , 1 ​ ( 𝐱 1 ) subscript 𝑚 𝑖 𝑛 1 subscript 𝐱 1 m_{in,1}(\mathbf{x}_{1}) ; note that in formula no. 4 N 𝑁 N denotes the size of the vector 𝐱 1 subscript 𝐱 1 \mathbf{x}_{1} , and that both m o ​ u ​ t ​ ( 𝐱 2 ) subscript 𝑚 𝑜 𝑢 𝑡 subscript 𝐱 2 m_{out}(\mathbf{x}_{2}) and f ​ ( 𝐱 1 , 𝐱 2 ) 𝑓 subscript 𝐱 1 subscript 𝐱 2 f(\mathbf{x}_{1},\mathbf{x}_{2}) are independent of 𝐱 2 subscript 𝐱 2 \mathbf{x}_{2} .

Formula no.	$m_{in}\left({𝐱}_1\right)$	$f({𝐱}_1,{𝐱}_2)$	$m_{out}\left({𝐱}_2\right)$
1	$𝒩({𝐱}_1;{\eta }_1,{𝐂}_1)$	$𝒩({𝐱}_2;{\mathrm{𝐀𝐱}}_1+𝐠,{𝐂}_2)$	$𝒩({𝐱}_2;𝐀{\eta }_1+𝐠,{\mathrm{𝐀𝐂}}_1{𝐀}_l^T+{𝐂}_2)$
2	$\delta \left({𝐱}_1-𝐚\right)$	$𝒩({𝐱}_2;{\mathrm{𝐀𝐱}}_1+𝐠,{𝐂}_2)$	$𝒩({𝐱}_2;\mathrm{𝐀𝐚}+𝐠,{𝐂}_2)$
3	$\delta \left({𝐱}_1-𝐚\right)$	$𝒩({𝐱}_1;{\mathrm{𝐀𝐱}}_2,{𝐂}_2)$	$𝒩(𝐚;{\mathrm{𝐀𝐱}}_2,{𝐂}_2)$
4	$𝒩({𝐱}_1;{\eta }_1,{𝐂}_1)$	$𝒩({𝐱}_1;{\eta }_2,{𝐂}_2)$	$\begin{array}{c}K\exp \left\{\frac{1}{2}\left[{\eta }^T𝐖\eta -{\eta }_1^T{𝐖}_1{\eta }_1-{\eta }_2^T{𝐖}_2{\eta }_2\right]\right\}\\ 𝐰={𝐖}_1{\eta }_1+{𝐖}_2{\eta }_2,𝐖={𝐖}_1+{𝐖}_2,\\ K={\left(det\left({𝐂}_1+{𝐂}_2\right)\right)}^{-N/2}\end{array}$
5	$𝒩({𝐱}_1;{\eta }_1,{𝐂}_1)$	$𝒩({𝐱}_1;𝐠+𝐀{𝐱}_2,{𝐂}_2)$	$\begin{array}{c}𝒩({𝐱}_2;\eta ,𝐂)\\ 𝐰={𝐀}^T{𝐖}_2\left[{𝐂}_3{𝐖}_1{\eta }_1-\left[𝐈-{𝐂}_3{𝐖}_2\right]𝐠\right]\\ 𝐖={𝐀}^T{𝐖}_2\left[𝐈-{𝐂}_3{𝐖}_2\right]𝐀,{𝐂}_3\triangleq {\left[{𝐖}_1+{𝐖}_2\right]}^{-1}\end{array}$

Equations

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

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

$${\mathbf{y}}_l\triangleq [y_{0,l},y_{1,l},...,y_{P-1,l}{]}^T={\mathbf{h}}_l\left({\mathbf{x}}_l^{(N)}\right)+{\mathbf{B}}_l\left({\mathbf{x}}_l^{(N)}\right){\mathbf{x}}_l^{(L)}+{\mathbf{e}}_l,$$

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

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

$$f\left({\mathbf{x}}_l,{\mathbf{y}}_{1:l}\right)=f\left({\mathbf{x}}_l,{\mathbf{y}}_{1:(l-1)}\right)f\left({\mathbf{y}}_l|{\mathbf{x}}_l\right)$$

$$f\left({\mathbf{x}}_{l+1},{\mathbf{y}}_{1:l}\right)=\int f\left({\mathbf{x}}_{l+1}|{\mathbf{x}}_l\right)f\left({\mathbf{x}}_l,{\mathbf{y}}_{1:l}\right)d{\mathbf{x}}_l,$$

$$f\left({\mathbf{y}}_{(l+1):T}|{\mathbf{x}}_l\right)=\int f\left({\mathbf{y}}_{(l+1):T}|{\mathbf{x}}_{l+1}\right)f\left({\mathbf{x}}_{l+1}|{\mathbf{x}}_l\right)d{\mathbf{x}}_l$$

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

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

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

$${\overleftarrow{m}}_{bp}\left({\mathbf{x}}_l\right)$$

$$=\int f({\mathbf{y}}_{(l+1):T}|{\mathbf{x}}_{l+1})f\left({\mathbf{x}}_{l+1}|{\mathbf{x}}_l\right)d{\mathbf{x}}_l$$

$$=f\left({\mathbf{y}}_{(l+1):T}|{\mathbf{x}}_l\right)\text{.}$$

$$f\left({\mathbf{y}}_l|{\mathbf{x}}_l\right){\overleftarrow{m}}_{bp}\left({\mathbf{x}}_l\right)=f\left({\mathbf{y}}_l|{\mathbf{x}}_l\right)f\left({\mathbf{y}}_{(l+1):T}|{\mathbf{x}}_l\right)$$

$$=f\left({\mathbf{y}}_{l:T}|{\mathbf{x}}_l\right)={\overleftarrow{m}}_{be}\left({\mathbf{x}}_l\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)\text{,}$$

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

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

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

$${\vec{m}}_{fp,j}\left({\mathbf{x}}_l^{(L)}\right)\triangleq \mathcal{N}\left({\mathbf{x}}_l^{(L)};{\eta }_{fp,l,j}^{(L)},{\mathbf{C}}_{fp,l,j}^{(L)}\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}}_{be}\left({\mathbf{x}}_{l+1}^{(L)}\right)\triangleq \mathcal{N}\left({\mathbf{x}}_{l+1}^{(L)};{\eta }_{be,l+1}^{(L)},{\mathbf{C}}_{be,l+1}^{(L)}\right),$$

$$m_{out}\left({\mathbf{x}}_2\right)=\int m_{in}\left({\mathbf{x}}_1\right)f\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)d{\mathbf{x}}_1$$

$${\overleftarrow{m}}_{1,j}\left({\mathbf{x}}_l^{(L)}\right)$$

$$\cdot {\overleftarrow{m}}_{be}\left({\mathbf{x}}_{l+1}^{(L)}\right){\vec{m}}_{fp,j}\left({\mathbf{x}}_l^{(N)}\right)d{\mathbf{x}}_{l+1}^{(L)}d{\mathbf{x}}_l^{(N)}$$

$$=\mathcal{N}\left({\mathbf{x}}_l^{(L)};{\eta }_{1,l,j}^{(L)},{\mathbf{C}}_{1,l,j}^{(L)}\right),$$

$${\mathbf{w}}_{1,l,j}^{(L)}$$

$$\cdot \left[{\bar{\mathbf{C}}}_{l+1}{\mathbf{w}}_{be,l+1}^{(L)}-{\mathbf{P}}_l^{(L)}{\mathbf{f}}_{l,j}^{(L)}\right],$$

$${\mathbf{W}}_{1,l,j}^{(L)}\triangleq {\left({\mathbf{C}}_{1,l,j}^{(L)}\right)}^{-1}={\left({\mathbf{A}}_{l,j}^{(L)}\right)}^T{\mathbf{W}}_w^{(L)}{\mathbf{P}}_l^{(L)}{\mathbf{A}}_{l,j}^{(L)},$$

$${\vec{m}}_j\left({\mathbf{z}}_l^{(L)}\right)=f\left({\mathbf{z}}_l^{(L)}|{\mathbf{x}}_{l/(l-1),j}^{(N)},{\tilde{\mathbf{x}}}_{l+1}^{(N)}\right)=\delta \left({\mathbf{z}}_l^{(L)}-{\mathbf{z}}_{l,j}^{(L)}\right),$$

$${\overleftarrow{m}}_{2,j}\left({\mathbf{x}}_l^{(L)}\right)$$

$$\cdot {\overleftarrow{m}}_j\left({\mathbf{z}}_l^{(L)}\right){\vec{m}}_{fp,j}\left({\mathbf{x}}_l^{(N)}\right)d{\mathbf{x}}_l^{(N)}d{\mathbf{z}}_l^{(L)}$$

$$=\mathcal{N}\left({\mathbf{z}}_{l,j}^{(L)};{\mathbf{A}}_{l,j}^{(N)}{\mathbf{x}}_l^{(L)},{\mathbf{C}}_w^{(N)}\right),$$

$${\overleftarrow{m}}_{3,j}\left({\mathbf{x}}_l^{(L)}\right)$$

Full text

Giorgio M. Vitetta, Emilio Sirignano and Francesco Montorsi

University of Modena and Reggio Emilia

Department of Engineering ”Enzo Ferrari”

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

email: giorgio.vitetta@unimore.it, emilio.sirignano@unimore.it, francesco.montorsi@gmail.com

Keywords: State Space Representation, Hidden Markov Model, Filtering, Smoothing, Marginalized Particle Filter, Belief Propagation.

1 Introduction

Bayesian filtering and Bayesian smoothing for state space models (SSMs) are two interrelated problems that have received significant attention for a number of years [1]. Bayesian filtering allows to recursively estimate, through a prediction/update mechanism, the probability density function (pdf) of the current state of any SSM, given the history of some observed data up to the current time. Unluckily, the general formulas describing the Bayesian filtering recursion (e.g., see [2, eqs. (4)-(5)]) admit closed form solutions for linear Gaussian and linear Gaussian mixture SSMs [1] only. On the contrary, approximate solutions are available for general nonlinear models; these are based on sequential Monte Carlo (SMC) techniques (also known as particle filtering methods) which represent a powerful tool for numerical approximations [3]-[5].

Bayesian smoothing, instead, exploits an entire batch of measurements to generate a significantly better estimate of the pdf (i.e., a smoothed or smoothing pdf) of SSM state over a given observation interval. Two general methods are available in the literature for recursively calculating smoothing densities, namely the forward filtering-backward smoothing recursion [4], [7] and the method based on the two-filter smoothing formula [8]-[10]. In both cases the computation of smoothing densities requires combining the predicted and/or filtered densities generated by a standard Bayesian filtering method with those produced by a recursive backward technique (known as backward information filtering, BIF, in the case of two-filter smoothing). Similarly as filtering, closed form solutions for Bayesian smoothing are available for linear Gaussian and linear Gaussian mixture models [1], [11]. This has motivated the development of various SMC approximations (also known as particle smoothers) for the above mentioned two methods in the case of nonlinear SSMs (e.g., see [4], [6], [8], [9], [12]-[15] and references therein).

While SMC methods can be directly applied to an arbitrary nonlinear SSM for both filtering and smoothing, it has been recognized that their estimation accuracy can be improved in the case of conditionally linear Gaussian (CLG) SSMs. In fact, the linear substructure of such models can be marginalised, so reducing the dimension of their SMC space [16], [17]. This idea has led to the development of important SMC techniques for filtering and smoothing, known as Rao-Blackwellized particle filtering (also dubbed marginalized particle filtering, MPF) [17] and Rao-Blackwellized particle smoothing (RBPS) [13], [14], [19], respectively.

Recently, the filtering problem for CLG SSMs has been investigated from a factor graph (FG) perspective in [20], where a novel interpretation of MPF as a forward only message passing algorithm over a specific FG has been provided and a novel extension of it, dubbed turbo filtering (TF), has been derived. In this manuscript, the same conceptual approach is employed to provide new insights in the fixed-interval smoothing problem [13] and to develop a novel solution for it. The proposed solution is represented by a novel RBPS method (dubbed Rao-Blackwellized serial smoothing, RBSS) having the following relevant features: a) it can be derived applying the well known sum-product algorithm (SPA) [22], [23], together with a specific scheduling procedure, to the same FG developed in [20] for a CLG SSM; b) unlike the RBPS methods devised in [13] and [14], it can be employed for a SSM in which both the linear and nonlinear state components influence each another; c) its computational complexity is appreciably smaller than that required by the other RBPS techniques; d) it benefits, unlike all the other RBPS techniques, from the exploitation of all the available pseudo-measurements and the ex novo computation of the weights for the particles generated in its forward recursion; e) it can be easily modified to compute the joint smoothing distribution over the entire observation interval (the resulting algorithm is called extended RBSS, ERBSS, in the following). Our simulation results evidence that, for the considered SSM, RBSS achieves a good accuracy-complexity tradeoff and that, in particular, it is slightly outperformed by ERBSS in state estimation accuracy, which, however, at the price, however, of a substantially higher computational cost.

It is worth mentioning that the application of FG methods to Bayesian smoothing is not new. However, as far as we know, the few results available in the technical literature about this topic refer to the case of linear Gaussian SSMs only [22], [24], [25], whereas we exclusively focus on the case in which the mathematical laws expressing state dynamics and/or available observations are nonlinear.

The remaining part of this manuscript is organized as follows. The model of the considered CLG SSM is briefly illustrated in Section 2. A representation of the smoothing problem through Forney-style FGs for both an arbitrary SSM and a CLG SSM is provided in Section 3. In Section 4 the RBSS technique is developed applying the SPA and proper message scheduling strategies to the FG derived for a CLG SSM; moreover, it is shown how it can be modified to estimate a point mass approximation of the joint smoothing distribution. Our FG-based smoothing algorithms are compared, in terms of accuracy and computational effort, in Section 5. Finally, some conclusions are offered in Section 6.

Notations: The probability density function (pdf) of a random vector $\mathbf{R}$ evaluated at point $\mathbf{r}$ is denoted $f(\mathbf{r})$; $\mathcal{N}\left(\mathbf{r};\mathbf{\eta_{r}},\mathbf{C_{r}}\right)$ represents the pdf of a Gaussian random vector $\mathbf{R}$ characterized by the mean $\mathbf{\eta_{r}}$ and covariance matrix $\mathbf{\mathbf{C_{r}}}$ evaluated at point $\mathbf{r}$; the *precision *(or weight) matrix associated with the covariance matrix $\mathbf{\mathbf{C_{r}}}$ is denoted $\mathbf{\mathbf{W_{r}}}$, whereas the transformed mean vector $\mathbf{\mathbf{W_{r}}\eta_{r}}$ is denoted $\mathbf{\mathbf{w_{r}}}$.

2 System Model

In the following we focus on the discrete-time CLG SSM described in [20], [21]. In brief, the SSM hidden state in the $l$-th interval is represented by the $D$-dimensional real vector $\mathbf{x}_{l}\triangleq[x_{0,l},x_{1,l},...,$ $x_{D-1,l}]^{T}$; this is partitioned in a) its $D_{L}$-dimensional *linear component *$\mathbf{x}_{l}^{(L)}\triangleq[x_{0,l}^{(L)},x_{1,l}^{(L)},...,x_{D_{L}-1,l}^{(L)}]^{T}$ and b) its $D_{N}$-dimensional *nonlinear component *$\mathbf{x}_{l}^{(N)}\triangleq[x_{0,l}^{(N)},x_{1,l}^{(N)},...,x_{D_{N}-1,l}^{(L)}]^{T}$ (with $D_{L}<D$ and $D_{N}=D-D_{L}$). The update equations of the linear and nonlinear components are given by

[TABLE]

and

[TABLE]

respectively; here, $\mathbf{f}_{l}^{(L)}\left(\mathbf{x}\right)$ ($\mathbf{f}_{l}^{(N)}\left(\mathbf{x}\right)$) is a time-varying $D_{L}$-dimensional ($D_{N}$-dimensional) real function, $\mathbf{A}_{l}^{(L)}(\mathbf{x}_{l}^{(N)})$ ($\mathbf{A}_{l}^{(N)}(\mathbf{x}_{l}^{(N)})$) is a time-varying $D_{L}\times D_{L}$ ($D_{N}\times D_{L}$) real matrix and $\mathbf{w}_{l}^{(L)}$ ($\mathbf{w}_{l}^{(N)}$) is the $l$-th element of the process noise sequence $\{\mathbf{w}_{k}^{(L)}\}$ ($\{\mathbf{w}_{k}^{(N)}\}$), which consists of $D_{L}$- dimensional ($D_{N}$-dimensional) independent and identically distributed (iid) noise* *vectors (statistical independence between $\{\mathbf{w}_{k}^{(L)}\}$ and $\{\mathbf{w}_{k}^{(N)}\}$ is also assumed for simplicity). Moreover, in the $l$-th interval some noisy observations, collected in the measurement vector

[TABLE]

are available about $\mathbf{x}_{l}$; here, $\mathbf{B}_{l}(\mathbf{x}_{l}^{(N)})$ is a time-varying $P\times D_{L}$ real matrix, $\mathbf{h}_{l}(\mathbf{x}_{l}^{(N)})$ is a time-varying $P$-dimensional real function and $\mathbf{e}_{l}$ the $l$-th element of the measurement noise sequence $\left\{\mathbf{e}_{k}\right\}$ consisting of $P$-dimensional iid noise vectors and independent of both $\{\mathbf{w}_{k}^{(N)}\}$ and $\{\mathbf{w}_{k}^{(L)}\}$. In the following Section we mainly focus on the so-called fixed-interval *smoothing problem *[13]; this consists of computing the sequence of posterior densities$\{f(\mathbf{x}_{l}|\mathbf{y}_{1:N}),\,l=1,2,...,T\}$ (where $T$ represents the length of the observation interval), given a) the initial pdf $f(\mathbf{x}_{1})$ and b) the $T\cdot P$-dimensional measurement vector $\mathbf{y}_{1:T}=\left[\mathbf{y}_{1}^{T},\mathbf{y}_{2}^{T},...,\mathbf{y}_{T}^{T}\right]^{T}$.

3 A FG-Based Representation of the Smoothing

Problem

In this Section we formulate the computation of the marginal smoothed density $f(\mathbf{x}_{l}|\mathbf{y}_{1:T})$ (with $l=1,2,...,T$) as a message passing algorithm over a specific FG for the following two cases: C.1) a SSM whose statistical behavior is characterized by the Markov model $f(\mathbf{x}_{l+1}|\mathbf{x}_{l})$ and the observation model $f(\mathbf{y}_{l}|\mathbf{x}_{l})$; C.2) a SSM having the additional property of being CLG (see the previous Section).

In case C.1 we take into consideration the joint pdf $f(\mathbf{x}_{l},\mathbf{y}_{1:T})$ in place of the posterior pdf $f(\mathbf{x}_{l}|\mathbf{y}_{1:T})$. This choice is motivated by the fact that: a) the computation of the former pdf can be easily formulated as a recursive message passing algorithm over a proper FG, since, as shown below, this involves only products and sums of products; b) the former pdf, being proportional to the latter one, is represented by the same FG (this issue is discussed in [22, Sec. II, p. 1297]). Note that the validity of statement a) relies on the following mathematical results: a) the factorization (e.g., see [8, Sec. 3])

[TABLE]

for the pdf of interest; b) the availability of recursive methods, known as Bayesian filtering [2] (and called forward filtering, FF, in the following for clarity) and *backward information filtering *(BIF; e.g., see [8]) for computing the joint pdf $f(\mathbf{x}_{l},\mathbf{y}_{1:(l-1)})$ and the conditional pdf $f(\mathbf{y}_{l:T}|\mathbf{x}_{l})$, respectively, for any $l$.

As far as FF is concerned, the formulation illustrated in [20, Sec. 2] is adopted here; this consists of a measurement update (MU) step followed by a time update (TU) step and assumes the a priori knowledge of the pdf $f(\mathbf{x}_{1})$ for its initialization. In the MU step of its $l$-th recursion (with $l=1,2,...,T$) the joint pdf

[TABLE]

is computed on the basis of pdf $f(\mathbf{x}_{l},\mathbf{y}_{1:(l-1)})$, and the new measurement vector $\mathbf{y}_{l}$. In the TU step, instead, the pdf $f\left(\mathbf{x}_{l},\mathbf{y}_{1:l}\right)$ (5) is exploited to compute the pdf

[TABLE]

representing a prediction about the future state $\mathbf{x}_{l+1}$.

A conceptually similar recursive procedure can be easily developed for the $(T-l)$-th recursion of BIF (with $l=T-1,T-2,...,1$). In fact, this can be formulated as a TU step followed by a MU step; these are expressed by

[TABLE]

and

[TABLE]

respectively. Note that this procedure requires the knowledge of the pdf $f(\mathbf{y}_{T}|\mathbf{x}_{T})$ for its initialization (see (7)).

Eqs. (5)-(8) show that each of the FF (or BIF) recursions involves only products of pdfs and a sum (i.e., an integration) of products. For this reason, based on the general rules about graphical models illustrated in [22, Sect. II], such recursions can be interpreted as specific instances of the SPA111In a Forney-style FG, such a rule can be formulated as follows [22]: the message emerging from a node f along some edge x is formed as the product of f and all the incoming messages along all the edges that enter the node f except x, summed over all the involved variables except x. applied to the cycle free FG of Fig. 1 (where the simplified notation of [22] is employed).

More specifically, it is easy to show that eqs. (5) and (6) can be seen as a SPA-based algorithm for forward message passing over the FG shown in Fig. 1 (the flow of forward messages is indicated by red arrows in the considered figure). In fact, if the FG is fed by the message222In the following the acronyms be, fp and sm are employed in the subscripts of various messages, so that readers can easily understand their meaning; in fact, the messages these acronyms refer to represent a form of backward estimation, forward prediction and smoothing, respectively.

[TABLE]

the forward messages emerging from the equality node and that passed along the edge associated with $\mathbf{x}_{l+1}$ are given by $\vec{m}_{fe}\left(\mathbf{x}_{l}\right)=f\left(\mathbf{x}_{l},\mathbf{y}_{1:l}\right)$ and $f(\mathbf{x}_{l+1},\mathbf{y}_{1:l})=\vec{m}_{fp}\left(\mathbf{x}_{l+1}\right)$, respectively [20], [21]. A similar interpretation can be provided for eqs. (7) and (8), which, however, can be reformulated as a SPA-based algorithm for backward message passing over the considered FG. In fact, if the input message

[TABLE]

enters the FG along the half edge associated with $\mathbf{x}_{l+1}$ (the flow of backward messages is indicated by blue arrows in Fig. 1), the backward message $\overset{\leftarrow}{m}_{bp}\left(\mathbf{x}_{l}\right)$  emerging from the node associated with the pdf $f(\mathbf{x}_{l+1}|\mathbf{x}_{l})$ is given by (see (7))

[TABLE]

Therefore, the message going out of the equality node in the backward direction can be evaluated as (see (8) and (10))

[TABLE]

and this concludes our proof.

These results easily lead to the conclusion that, once the forward and backward message passing algorithms illustrated above have been carried out over the entire observation interval, the smoothed pdf $f\left(\mathbf{x}_{l},\mathbf{y}_{1:T}\right)$ can be evaluated as (see (4), (9) and (12))

[TABLE]

with $l=1,2,...,T$ (note that $\overset{\leftarrow}{m}_{be}\left(\mathbf{x}_{T}\right)=1$ and $\vec{m}_{fp}\left(\mathbf{x}_{1}\right)=f(\mathbf{x}_{1})$)

The FG we develop for case C.2 is based not only on that analysed for case C.1, but also on the idea of representing a mixed linear/nonlinear SSM as the concatenation of two interacting sub-models, one referring to the linear component of system state, the other one to its nonlinear component [20]. This suggests to decouple the smoothing problem for $\mathbf{x}_{l}^{(L)}$ from that for $\mathbf{x}_{l}^{(N)}$, i.e. the evaluation of * *$f(\mathbf{x}_{l}^{(L)}|\mathbf{y}_{1:T})$ from that of $f(\mathbf{x}_{l}^{(N)}|\mathbf{y}_{1:T})$. In practice, from a graphical viewpoint, two sub-graphs, one referring to smoothing for $\mathbf{x}_{l}^{(L)}$, the other one to smoothing for $\mathbf{x}_{l}^{(N)}$, are developed first; then, they are merged by adding five distinct equality nodes, associated with the variables (namely, $\mathbf{y}_{l}$, $\mathbf{x}_{l}^{(L)}$, $\mathbf{x}_{l}^{(N)}$, $\mathbf{x}_{l+1}^{(L)}$ and $\mathbf{x}_{l+1}^{(N)}$) shared by such sub-graphs. This leads to the FG illustrated in Fig. 2, in which the sub-graph referring to the linear (nonlinear) state component is identified by red (blue) lines, whereas the equality nodes added to merge them are identified by black lines. Note that the sub-graph for the linear (nonlinear) component is derived under the assumption that the nonlinear (linear) component is known. Consequently, smoothing for the linear component $\mathbf{x}_{l}^{(L)}$ can benefit not only from the measurement $\mathbf{y}_{l}$, but also from the so called pseudo-measurement (see (2))

[TABLE]

which, from a statistical viewpoint, is characterized by the pdf $f(\mathbf{z}_{l}^{(L)}|\mathbf{x}_{l}^{(L)},\mathbf{x}_{l}^{(N)})$. Similarly, the pseudo-measurement (see (1))

[TABLE]

characterized by the pdf $f(\mathbf{z}_{l}^{(N)}|\mathbf{x}_{l}^{(N)})$, can be exploited in smoothing for the nonlinear component $\mathbf{x}_{l}^{(N)}$.  These considerations explain why the upper (lower) sub-graph shown in Fig. 2 contains an additional node representing the pdf $f(\mathbf{z}_{l}^{(L)}|\mathbf{x}_{l}^{(L)},\mathbf{x}_{l}^{(N)})$ ($f(\mathbf{z}_{l}^{(N)}|\mathbf{x}_{l}^{(N)})$) and a specific node not referring to the above mentioned pdf factorizations, but representing the transformation from the couple $(\mathbf{x}_{l}^{(N)},\mathbf{x}_{l+1}^{(N)})$ to $\mathbf{z}_{l}^{(L)}$ ($(\mathbf{x}_{l}^{(L)},\mathbf{x}_{l+1}^{(L)})$ to $\mathbf{z}_{l}^{(N)}$); the last peculiarity, evidenced by the presence of an arrow on all the edges connected to such a node, has to be carefully kept into account when deriving message passing algorithms.

Given the FG of Fig. 2, we would like to follow the same line of reasoning as that illustrated for the graphical model of Fig. 1. In particular, given the input backward messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(L)})\triangleq f(\mathbf{y}_{(l+1):T},\mathbf{z}_{(l+1):T}^{(L)},\mathbf{x}_{l+1}^{(L)})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(N)})\triangleq f(\mathbf{y}_{(l+1):T},\mathbf{z}_{(l+1):T}^{(N)},\mathbf{x}_{l+1}^{(N)})$, we would like to derive a BIF algorithm based on this FG (FF has already been investigated in [20] and [21]) and generating the output backward messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(L)})=f(\mathbf{y}_{l:T},\mathbf{z}_{l:T}^{(L)},\mathbf{x}_{l}^{(L)})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(N)})=f(\mathbf{y}_{l:T},\mathbf{z}_{l:T}^{(N)},\mathbf{x}_{l}^{(N)})$ on the basis of the available a priori information and the noisy measurement $\mathbf{y}_{l}$. Unluckily, the new FG, unlike the one represented in Fig. 1, is not cycle-free, so that any application of the SPA to it unavoidably leads to* approximate solutions* [23], whatever message scheduling procedure is adopted. In the following Section we show that the RBSS technique we propose represents one of such solutions.

4 Particle Smoothing as Message Passing

In this Section we first illustrate some assumptions about the statistical properties of the SSM defined in Section 2. Then, we develop the RBSS technique and compare its most relevant features with those of the other RBPS algorithms available in the technical literature. Finally, we show how this technique can be modified to estimate the *joint smoothing density *$f(\mathbf{x}_{1:T}|\mathbf{y}_{1:T})$.

4.1 Statistical properties of the considered SSM

Even if the FG representation shown in Fig. 2 can be employed for any mixed linear/nonlinear system described by eqs. (1)-(3), the methods derived in this Section apply, like MPF [17] and TF [20], to the specific class of GLG* *SSMs. For this reason, following [20], [21] we assume that: a) the process noise $\{\mathbf{w}_{k}^{(L)}\}$ ($\{\mathbf{w}_{k}^{(N)}\}$) is Gaussian and all its elements have zero mean and covariance $\mathbf{C}_{w}^{(L)}$ ($\mathbf{C}_{w}^{(N)}$) for any $l$; b) the measurement noise$\ \{\mathbf{e}_{k}^{(L)}\}$ is Gaussian having zero mean and covariance matrix $\mathbf{C}_{e}$ for any $l$; c) all the above mentioned Gaussian processes are statistically independent. Under these assumptions, the pdfs $f(\mathbf{y}_{l}|\mathbf{x}_{l}^{(L)},\mathbf{x}_{l}^{(N})$, $f(\mathbf{z}_{l}^{(L)}|\mathbf{x}_{l}^{(L)})$ and $f(\mathbf{x}_{l+1}^{(L)}|\mathbf{x}_{l}^{(L)},\mathbf{x}_{l}^{(N)})$ are Gaussian with mean (covariance matrix) $\mathbf{B}_{l}(\mathbf{x}_{l}^{(N)})\mathbf{x}_{l}^{(L)}+\mathbf{h}_{l}(\mathbf{x}_{l}^{(N)})$, $\mathbf{A}_{l}^{(N)}(\mathbf{x}_{l}^{(N)})\,\mathbf{x}_{l}^{(L)}$ and $\mathbf{f}_{l}^{(L)}(\mathbf{x}_{l}^{(N)})+\mathbf{A}_{l}^{(L)}(\mathbf{x}_{l}^{(N)})\,\mathbf{x}_{l}^{(L)}$, respectively ($\mathbf{C}_{e}$, $\mathbf{C}_{w}^{(N)}$ and $\mathbf{C}_{w}^{(L)}$, respectively). Similarly, the pdfs $f(\mathbf{z}_{l}^{(N)}|\mathbf{x}_{l}^{(N)})$ and $f(\mathbf{x}_{l+1}^{(N)}|\mathbf{x}_{l}^{(N)},\mathbf{x}_{l}^{(L)})$ are Gaussian with mean (covariance matrix) $\mathbf{f}_{l}^{(L)}(\mathbf{x}_{l}^{(N)})$ and $\mathbf{f}_{l}^{(N)}(\mathbf{x}_{l}^{(N)})+\mathbf{A}_{l}^{(N)}(\mathbf{x}_{l}^{(N)})\,\mathbf{x}_{l}^{(L)}$, respectively ($\mathbf{C}_{w}^{(L)}$ and $\mathbf{C}_{w}^{(N)}$, respectively).

4.2 Derivation of the Rao-Blacwellized serial

smoother

The FF algorithm employed in the forward pass of the proposed RBSS is represented by MPF333Note that TF can be employed in place of MPF in the forward pass of RBSS. However, our computer simulations have evidenced that, in the presence of strong measurement and/or process noise (like in the scenarios considered in Section 5), this choice doe not provide any performance improvement with respect to MPF.. In its $(l-1)$-th recursion (with $l=2,3,...,T$), the particle set $\{\mathbf{x}_{l/(l-1),j}^{(N)},j=0,1,...,N_{p}-1\}$, consisting of $N_{p}$ distinct particles, is predicted for the nonlinear state component $\mathbf{x}_{l}^{(N)}$ (TU for this component); the weight $w_{l/(l-1),j}$ assigned to the particle $\mathbf{x}_{l/(l-1),j}^{(N)}$ is equal to $1/N_{p}$ for any $j$, since the use of particle resampling in each recursion is assumed. The particle weights are updated in the MU of the following (i.e., $l$-th) recursion on the basis of the new measurement $\mathbf{y}_{l}$ (MU for the nonlinear component): the new weights are denoted $\{w_{l/l,j},j=0,1,...,N_{p}-1\}$ in the following and, generally speaking, are all different. This is followed by particle resampling, that generates the new particle set $\{\mathbf{x}_{l/l,j}^{(N)},j=0,1,...,N_{p}-1\}$ (usually containing multiple copies of the most likely particles of the set $\{\mathbf{x}_{l/(l-1),j}^{(N)}\}$). A conceptually similar procedure is followed for the linear state component, for which a particle-dependent Gaussian representation is adopted. In particular, in the following, the Gaussian model predicted for $\mathbf{x}_{l}^{(L)}$ in the $(l-1)$-th recursion (TU for the linear state component) and associated with $\mathbf{x}_{l/(l-1),j}^{(N)}$ is denoted $\mathcal{N}(\mathbf{x}_{l}^{(L)};\mathbf{\eta}_{fp,l,j}^{(L)},\mathbf{C}_{fp,l,j}^{(L)})$. Note that only a portion of these Gaussian models is usually updated in the MU of the next (i.e., $l$-th) recursion; in fact, this task follows particle resampling, which typically leads to discarding a fraction of the particles collected in the set $\{\mathbf{x}_{l/(l-1),j}^{(N)}\}$.

The recursive algorithm developed for the backward pass of the RBSS technique results from the application of the SPA to the FG shown in Fig. 2, and accomplishes BIF and smoothing (i.e., the merge of statistical information generated by FF and BIF). Each of its recursions consists of two parts, the first concerning the linear state component, the second one the nonlinear state component; moreover, these parts are executed serially. The message scheduling employed in the $(T-l)$-th recursion of BIF and smoothing (with $l=T-1,T-2,...,1$) is summarized in Fig. 3, where the edges involved in the first (second) part are identified by continuous (dashed) lines. Similarly to MPF, most of the processing tasks which both parts consist of can be formulated with reference to a single particle; this explains why the notation adopted for the messages appearing in Fig. 3 includes the subscript $j$, that represents the index of the particle (namely, the particle $\mathbf{x}_{l/(l-1),j}^{(N)}$) representing $\mathbf{x}_{l}^{(N)}$ within the considered recursion.

Before providing a detailed description of the messages passed in the graphical model of Fig. 3, all the messages feeding the considered recursion (i.e., its input messages) and those emerging from it (i.e., its output messages) must be defined. The input messages can be divided in two groups. The first group consists of the messages $\vec{m}_{fp,j}(\mathbf{x}_{l}^{(L)})$ and $\vec{m}_{fp,j}(\mathbf{x}_{l}^{(N)})$, that are predicted the $(l-1)$-th recursion of the forward pass; the second one, instead, is made of the messages $\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l+1}^{(N)})$ and $\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l+1}^{(L)})$, that are generated in the $(T-l-1)$-th recursion of the backward pass. The messages of the first group are defined as

[TABLE]

and

[TABLE]

and can be interpreted as the $j$-th hypothesis about a) the value (namely, $\mathbf{x}_{l/(l-1),j}^{(N)}$) taken on by the (hidden) nonlinear state component $\mathbf{x}_{l}^{(N)}$ and b) the statistical representation of the (hidden) linear state component $\mathbf{x}_{l}^{(L)}$ associated with such a value, respectively. In the $l$-th recursion of FF, the likelihood of this hypothesis is assessed by evaluating the above mentioned weight $w_{l/l,j}$; such a weight, however, is ignored in the backward pass. This choice is motivated by the our belief that, if such a weight is computed ex novo, its accuracy can be improved thanks to the availability of both more refined (i.e., smoothed) statistical information about $\mathbf{x}_{l}^{(L)}$ and additional (backward) information about $\mathbf{x}_{l+1}^{(N)}$ (see (18) and (19) below).

The input messages of the second group are defined as

[TABLE]

and

[TABLE]

and represent part of the statistical information generated in the previous (i.e., the $(T-l-1)$-th) recursion of the backward pass. In particular, as explained in detail below, the messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(N)})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(L)})$ convey the final estimate $\mathbf{x}_{be,l+1}^{(N)}$ (i.e., a single particle representation) of $\mathbf{x}_{l+1}^{(N)}$ and a simplified statistical representation of $\mathbf{x}_{l+1}^{(L)}$, respectively. This explains why the RBSS, in the $(T-l)$-th recursion of its backward pass, processes the input messages (16)-(19) to compute an estimate, denoted $\mathbf{x}_{be,l}^{(N)}$, of $\mathbf{x}_{l}^{(N)}$ and a simplified statistical model, denoted $\mathcal{N}(\mathbf{x}_{l}^{(L)};\mathbf{\eta}_{be,l}^{(L)},\mathbf{C}_{be,l}^{(L)})$, for $\mathbf{x}_{l}^{(L)}$; these information are conveyed by the output messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(N)})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(L)})$, respectively. The evaluation of these messages is based, as already mentioned above, on the scheduling illustrated in Fig. 3 and on the formulas listed in Tables 1 and 2 (actually, the only formulas missing in these Tables are those employed in the evaluation of the message $\overset{\leftarrow}{m}_{j}(\mathbf{z}_{l}^{(N)})$ (42)  and, in particular, of its parameters $\mathbf{\eta}_{\mathbf{z},l,j}^{(N)}$ (44) and $\mathbf{C}_{\mathbf{z},l,j}^{(N)}$ (45); mathematical details about this can be found in [20, Sec. 6]). Such formulas refer to the computation of the message $m_{out}\left(\mathbf{x}\right)=m_{in,1}\left(\mathbf{x}\right)m_{in,2}\left(\mathbf{x}\right)$ (emerging from an equality node fed by the messages $m_{in,1}\left(\mathbf{x}\right)$ and $m_{in,2}\left(\mathbf{x}\right)$) and

[TABLE]

(emerging from a function node $f\left(\mathbf{x}_{1},\mathbf{x}_{2}\right)$ fed by the message $m_{in,}\left(\mathbf{x}_{1}\right)$), respectively; moreover, they are provided by [22, Table 2, p. 1303] or can be easily derived on the basis of standard mathematical results about Gaussian random variables. For this reason, in the following description of the RBSS backward pass, we provide, for each message, a simple code identifying the specific formula on which its evaluation is based; in particular, the notation TX-Y is employed to identify formula no. Y appearing in Table X. Moreover, to ease the interpretation of the proposed signal processing tasks executed within the RBSS algorithm, the message passing accomplished in the considered recursion is divided in the seven steps described below; steps 1-3 and steps 4-6 refer to the two parts of the message passing shown in Fig. 3, whereas the last step concern the evaluation of: a) the smoothed pdf of $\mathbf{x}_{l}$ and the pdfs of its components; b) the output messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(N)})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(L)})$.

1. Time update for $\mathbf{x}_{l}^{(L)}$ - Compute the message (see T2-5, (16) and (19))

[TABLE]

where

[TABLE]

[TABLE]

$\mathbf{A}_{l,j}^{(L)}\triangleq\mathbf{A}_{l}^{(L)}(\mathbf{x}_{l/(l-1),j}^{(N)})$, $\mathbf{W}_{w}^{(L)}\triangleq(\mathbf{C}_{w}^{(L)})^{-1}$, $\mathbf{P}_{l}^{(L)}\triangleq\mathbf{I}_{D_{L}}-\mathbf{\bar{C}}_{l+1}\mathbf{W}_{w}^{(L)}$, $\mathbf{\bar{C}}_{l+1}\triangleq(\mathbf{W}_{w}^{(L)}+\mathbf{W}_{be,l+1}^{(L)})^{-1}$, $\mathbf{W}_{be,l+1}^{(L)}\triangleq(\mathbf{C}_{be,l+1}^{(L)})^{-1}$, $\mathbf{f}_{l,j}^{(L)}\triangleq\mathbf{f}_{l}^{(L)}(\mathbf{x}_{l/(l-1),j}^{(N)})$ and $\mathbf{w}_{be,l+1}^{(L)}\triangleq\mathbf{W}_{be,l+1}^{(L)}\mathbf{\eta}_{be,l+1}^{(L)}$.

2. Measurement update for $\mathbf{x}_{l}^{(L)}$ - Compute: a) the message

[TABLE]

where $\mathbf{z}_{l,j}^{(L)}\triangleq\mathbf{x}_{be,l+1}^{(N)}-\mathbf{f}_{l,j}^{(N)}$ and $\mathbf{f}_{l,j}^{(N)}\triangleq\mathbf{f}_{l}^{(N)}(\mathbf{x}_{l/(l-1),j}^{(N)})$; b) the messages (see T2-3, T1-3, T2-2 and T1-2, respectively; see also (16), (21) and (24))

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Here,

[TABLE]

[TABLE]

$\mathbf{A}_{l,j}^{(N)}\triangleq\mathbf{A}_{l}^{(N)}(\mathbf{x}_{l/(l-1),j}^{(N)})$, $\mathbf{W}_{w}^{(N)}\triangleq[\mathbf{C}_{w}^{(N)}]^{-1}$,

[TABLE]

[TABLE]

$\mathbf{B}_{l,j}\triangleq\mathbf{B}_{l}(\mathbf{x}_{l/(l-1),j}^{(N)})$, $\mathbf{h}_{l,j}\triangleq\mathbf{h}_{l}(\mathbf{x}_{l/(l-1),j}^{(N)})$, $\mathbf{W}_{e}\triangleq\mathbf{C}_{e}^{-1}$,

[TABLE]

and

[TABLE]

3. Merge of forward and backward messages about $\mathbf{x}_{l}^{(L)}$ - Compute the message (see (13), (17), (29), T1-2 and Fig. 3)

[TABLE]

where

[TABLE]

[TABLE]

$\mathbf{W}_{fp,l,j}^{(L)}\triangleq(\mathbf{C}_{fp,l,j}^{(L)})^{-1}$ and $\mathbf{w}_{fp,l,j}^{(L)}\triangleq\mathbf{W}_{fp,l,j}^{(L)}\mathbf{\eta}_{fp,l,j}^{(L)}$.

4. Time update for $\mathbf{x}_{l}^{(N)}$ - Compute the message (see T2-1, (18) and (36))

[TABLE]

where

[TABLE]

and

[TABLE]

5. Measurement update for $\mathbf{x}_{l}^{(N)}$ - Compute: a) the message

[TABLE]

and the message (see T3-1)

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$K_{l,j}=(\det(\mathbf{C}_{\mathbf{z},l,j}^{(N)}+\mathbf{C}_{w}^{(L)}))^{-D_{L}/2}$, $\mathbf{W}_{\mathbf{z},l,j}^{(N)}\triangleq(\mathbf{C}_{\mathbf{z},l,j}^{(N)})^{-1}$ and $\mathbf{w}_{\mathbf{z},l,j}^{(N)}\triangleq\mathbf{W}_{\mathbf{z},l,j}^{(N)}\mathbf{\eta}_{\mathbf{z},l,j}^{(N)}$; b) the messages (see T1-1 and T2-1, respectively)

[TABLE]

and

[TABLE]

where $\mathbf{\eta}_{4,l,j}^{(N)}=\mathbf{B}_{l,j}\mathbf{\eta}_{sm,l,j}^{(L)}+\mathbf{h}_{l,j}$ and $\mathbf{C}_{4,l,j}^{(N)}=\mathbf{B}_{l,j}\mathbf{C}_{sm,l,j}^{(L)}(\mathbf{B}_{l,j})^{T}+\mathbf{C}_{e}$.

6. Merge of forward and backward messages about $\mathbf{x}_{l}^{(N)}$ - This requires: a) computing the messages (see (48) and (49))

[TABLE]

and (see (16) and T1-1)

[TABLE]

b) normalising the weight set $\{W_{l,j}\,\}$, i.e. generating the new weight

[TABLE]

for $j=0,1,,...,N_{p}-1$; c) setting

[TABLE]

for $j=0,1,,...,N_{p}-1$.

7. *Generation of smoothed pdfs and input messages for the next recursion *- Compute: a) the pdfs

[TABLE]

[TABLE]

and

[TABLE]

that represent approximations of the marginal smoothed pdfs of $\mathbf{x}_{l},$ $\mathbf{x}_{l}^{(N)}$ and $\mathbf{x}_{l}^{(L)}$, respectively; b) the input messages

[TABLE]

and

[TABLE]

for the next recursion; here,

[TABLE]

[TABLE]

and

[TABLE]

After completing step 7, the $(T-l)$-th recursion of the RBSS technique is over. Then, the recursion index $l$ is decreased by one; if it equals zero, the backward pass is over, otherwise a new recursion is started. Note also that the first recursion of the backward pass requires the knowledge of its input messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{T}^{(N)})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{T}^{(L)})$, whose evaluation is based on the statistical information generated in the last recursion of the forward pass. In fact, in our work these messages are defined as

[TABLE]

and

[TABLE]

respectively; here, $\mathbf{x}_{fe,T}^{(N)}\triangleq\sum\limits_{j=0}^{N_{p}-1}w_{T/T,j}\,\mathbf{x}_{T/(T-1),j}^{(N)}$, whereas the parameters $\mathbf{\eta}_{fe,T}^{(L)}$ and $\mathbf{C}_{fe,T}^{(L)}$ of (63) are evaluated on the basis of formulas (60) and (61), but employing, in place of the Gaussian messages $\{\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l}^{(N)})\}$ (see (50)), the messages $\{\mathcal{N}(\mathbf{x}_{T}^{(L)};\mathbf{\eta}_{fe,T,j}^{(L)},\mathbf{C}_{fe,T,j}^{(L)})\}$ generated by the MU for the linear state component in the last (i.e., in the $T$-th) recursion of FF.

The RBSS algorithm illustrated above deserves various comments, that are listed below.

The message flow in the backward pass proceeds in a reverse order with respect to the forward pass (a similar scheduling in the backward pass has been adopted in [14]); in fact, in MPF the evaluation of particle weights and the prediction of new particles for the next recursion (accomplished in the MU and in the TU, respectively, for the nonlinear state component) precedes the MU and the TU for the linear state component. Moreover, unlike TF, a single pass is accomplished over the FG. 2. 2.

In step 1 a one-step ahead prediction is evaluated for $\mathbf{x}_{l}^{(L)}$ on the basis of the pdf of $\mathbf{x}_{l+1}^{(L)}$ (provided by the particle-independent message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(L)})$ (19)). A conceptually similar task is carried out for $\mathbf{x}_{l}^{(N)}$ in step 4. However, in the last case, pdf prediction does not involve the generation of new particles (like in the TU step of MPF), but only the computation of new weights for the particles originating from the forward pass. For this reason, the support of the pdf $\hat{f}(\mathbf{x}_{l}^{(N)}|\mathbf{y}_{1:N})$ (55) estimated for $\mathbf{x}_{l}^{(N)}$ in the backward pass remains exactly the same as that of the corresponding filtered pdf computed in the forward pass. 3. 3.

In step 2 the pdf $\overset{\leftarrow}{m}_{1,j}(\mathbf{x}_{l}^{(L)})$ (21) emerging from step 1 is refined on the basis of a) the measurement $\mathbf{y}_{l}$ and b) the pseudo-measurement $\mathbf{z}_{l,j}^{(L)}$, which depends on the particle index $j$ through $\mathbf{x}_{l,j}^{(N)}$ only (since a single particle is available for $\mathbf{x}_{l+1}^{(N)}$). Even if this entails a loss of diversity in the pseudo-measurement set $\{\mathbf{z}_{l,j}^{(L)}\}$ with respect to the corresponding set generated by MPF in the forward pass, the use of these quantities in state estimation is still beneficial. Incidentally, we note that no attention to the exploitation of pseudo-measurements $\mathbf{z}_{l}^{(L)}$ and $\mathbf{z}_{l}^{(N)}$ is paid in the development of the other RBPS methods available in the literature, even if these quantities are known to play an important role in state estimation [17], [20], [21]. 4. 4.

In step 3 the merge of the forward message $\vec{m}_{fp,j}(\mathbf{x}_{l}^{(L)})$ with the backward message $\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l}^{(L)})$ results in the ‘smoothed’ message $m_{sm,j}(\mathbf{x}_{l}^{(L)})$ (36), which is expected to provide a more refined statistical representation of $\mathbf{x}_{l}^{(L)}$ than $\vec{m}_{fp,j}(\mathbf{x}_{l}^{(L)})$ or $\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l}^{(L)})$ alone (under the assumption that $\mathbf{x}_{l}^{(N)}=\mathbf{x}_{l/(l-1),j}^{(L)}$) and, consequently, to improve the accuracy of the particle weights evaluated in steps 4 and 5; note also that $m_{sm,j}(\mathbf{x}_{1}^{(L)})=\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{1}^{(L)})$ and $m_{sm,j}(\mathbf{x}_{T}^{(L)})=\vec{m}_{fp,j}(\mathbf{x}_{T}^{(L)})$ should be assumed, since at the instant $l=1$ ($l=T$) only a backward estimate (a forward prediction) is available for $\mathbf{x}_{l}^{(L)}$. 5. 5.

In step 3 the equivalence between the expressions (27) and (28) is motivated by the fact that they differ by a scale factor and that scale factors can be always neglected in passing Gaussian messages [22]. 6. 6.

In step 5 the factors $w_{1,l,j}$, $w_{2,l,j}$ and $w_{4,l,j}$ of the overall weight $W_{l,j}$ (50) are related to the state transition $\mathbf{x}_{l+1}^{(N)}\rightarrow\mathbf{x}_{l}^{(N)}$, to the statistical representation of $\mathbf{z}_{l}^{(N)}$ (conveyed by the Gaussian message $\overset{\leftarrow}{m}_{j}(\mathbf{z}_{l}^{(N)})$ (42)) and to the measurement $\mathbf{y}_{l}$, respectively. Note also that: a) the weight $w_{1,l,j}$ depends on the (particle-independent) estimate $\mathbf{x}_{be,l+1}^{(L)}$, which can be interpreted as an additional pseudo-measurement originating from our knowledge of the future (and, consequently, unavailable in the forward pass); b) the weight $w_{2,l,j}$ (43) cannot be computed in the forward pass because of the scheduling adopted in MPF (the TU for the nonlinear state component represents the last step accomplished in each recursion of MPF); c) the weight $w_{4,l,j}$ corresponds to the weight $w_{l/l,j}$ computed by MPF in the forward pass but, as already mentioned at point 4), is expected to be more accurate thanks to the availability of more refined statistical information about $\mathbf{x}_{l}^{(L)}$ (conveyed by the message $m_{sm,j}(\mathbf{x}_{l}^{(L)})$ (36) in place of $m_{fp,j}(\mathbf{x}_{l}^{(L)})$ (17)). 7. 7.

Steps 1-6 need to be repeated $N_{p}$ times, once for each particle of the set $\{\mathbf{x}_{l/(l-1),j}^{(N)}\}$; in practice, this task can be parallelized, since the processing executed for any particle within these steps is not influenced from that carried out for all the other particles. 8. 8.

The expressions of the weights $w_{1,l,j}$, $w_{2,l,j}$ and $w_{4,l,j}$ have similar mathematical structure (see (39), (43) and (49), respectively) in the sense that they are given by the product of an exponential with a particle-dependent factor. An approximate evaluation of these weights can be obtained neglecting the contribution of a such a factor in each of their expressions. As a matter of fact, our computer simulations have evidenced that, at least for the considered SSM, this simplification does not entail a visible loss in RBSS accuracy. However, if used, it requires the adoption of weight normalization for each of the three weight sets; consequently, the overall weight $W_{l,j}$ (see (50)) is computed as

[TABLE]

where $\tilde{w}_{k,l,j}\,\triangleq w_{k,l,j}/\sum_{j=0}^{N_{p}-1}w_{k,l,j}$ for $k=1,2$ and $4$. 9. 9.

The final particle weights $\{W_{sm,l,j}\}$ (see (52)) are employed to generate both the final estimate $\mathbf{x}_{be,l}^{(N)}$ (59) of $\mathbf{x}_{l}^{(N)}$ and the $N_{p}$-component Gaussian mixture (GM) $\hat{f}(\mathbf{x}_{l}^{(L)}|\mathbf{y}_{1:N})$ (56), expressing our final estimate of the pdf of $\mathbf{x}_{l}^{(L)}$. This GM, however, is not passed to the next recursion as it is, since this would be make the complexity of our message passing algorithm unmanageable. This is the reason why this pdf is condensed in the Gaussian message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(L)})$ (58) by means of a standard transformation, expressed by formulas (60) and (61), and preserving both the mean and the covariance matrix of the GM itself (e.g., see [27, Sect. 4]).

Our final comment concerns the smoothing of the linear state component and has been inspired by the considerations illustrated in [19, Par. IV-D], where it is stressed that in Rao-Blackwellized methods the statistics for the linear state component need to be computed conditionally on the considered nonlinear state trajectories. As a matter of fact, our RBSS algorithm generates a single estimate of nonlinear state trajectory in its backward pass (the $l$-th point of this trajectory is represented by $\mathbf{x}_{be,l}^{(N)}$ with $l=1,2,...,T-1$ and by $\mathbf{x}_{fe,T}^{(N)}$ for $l=T$); however, the statistical models for the linear state components associated with this trajectory (see $\hat{f}(\mathbf{x}_{l}^{(L)}|\mathbf{y}_{1:N})$ (56) or its condensed representation $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(L)})$ (58)) do not satisfy the above mentioned condition, since they do not actually refer to a specific nonlinear state trajectory. This suggests that, once the RBSS algorithm has been carried out, more refined statistics for the linear state component could be computed by:

Carrying out, first of all, a new forward pass under the assumption that the nonlinear state component is known and, in particular, $\mathbf{x}_{l}^{(N)}=\mathbf{x}_{be,l}^{(N)}$ for $l=1,2,..,T-1$ and $\mathbf{x}_{T}^{(N)}=\mathbf{x}_{fe,T}^{(N)}$; this produces a single message $\vec{m}_{fp}(\mathbf{x}_{l}^{(L)})\triangleq\mathcal{N}(\mathbf{x}_{l}^{(L)};\mathbf{\eta}_{fp,l}^{(L)},\mathbf{C}_{fp,l}^{(L)})$ in place of the $N_{p}$ messages $\{\vec{m}_{fp,j}(\mathbf{x}_{l}^{(L)})\}$ (see (17)) for $l=2,..,T$. 2. 2.

Then, accomplishing a new backward pass under the same assumption as the previous point; this generates a single Gaussian message $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(L)})\triangleq\mathcal{N}(\mathbf{x}_{l}^{(L)};\mathbf{\eta}_{be,l}^{(L)},\mathbf{C}_{be,l}^{(L)})$ in place of the $N_{p}$ messages $\{\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l}^{(L)})\}$ (see (29)) for $l=T-1,T-2,..,1$ (note that $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{T}^{(L)})$ is still given by (63)). 3. 3.

Finally, merging $\vec{m}_{fp}(\mathbf{x}_{l}^{(L)})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(L)})$ in the message $m_{sm}(\mathbf{x}_{l}^{(L)})=\mathcal{N}(\mathbf{x}_{l}^{(L)};\mathbf{\eta}_{sm,l}^{(L)},\mathbf{C}_{sm,l}^{(L)})$, with $l=2,1,..,T-1$ ($m_{sm}(\mathbf{x}_{1}^{(L)})=\overset{\leftarrow}{m}_{be}(\mathbf{x}_{1}^{(L)})$ and $m_{sm}(\mathbf{x}_{T}^{(L)})=\vec{m}_{fp}(\mathbf{x}_{T}^{(L)})$ are assumed) on the basis of (36)-(38), so that a new final estimate $\mathbf{\eta}_{sm,l}^{(L)}$ is available for $\mathbf{x}_{l}^{(L)}$.

We believe that, even if this procedure is conceptually appealing, the improvement it may provide in the estimation accuracy for the linear state component is influenced by a) the number of modes of the density of $\mathbf{x}_{l}^{(L)}$ (since the adopted unimodal model for this state component might provide a poor statistical representation of it) and b) the presence of large errors, at specific instants, in the estimated nonlinear state trajectory.

4.3 Comparison of the RBSS algorithm with other RBPS

methods

Despite their substantially different structures, the other RBPS methods available in the technical literature [13], [14], [19] share the following relevant features: 1) the computation of an estimate of the joint smoothing density $f(\mathbf{x}_{1:T}|\mathbf{y}_{1:T})$; 2) the reuse of FF particles and weights; 3) the use of resampling in the generation of backward trajectories; 4) the exploitation of Kalman techniques for the linear state component. In the following we provide some details about these features, so that some important differences between such techniques and the RBSS algorithm can be easily understood.

The first feature refers to the fact that these techniques aim at generating realizations from the complete joint smoothing pdf $f(\mathbf{x}_{1:T}|\mathbf{y}_{1:T})$. Each realization consists of a) a trajectory (i.e., a set of $T$ particles, one for each observation instant) for the nonlinear state component and a set of $T$ Gaussian pdfs (one for each observation instant) [13], [19] or b) a trajectory for the entire state [14] (in this case a particle-based representation is adopted for the linear state component too). This approach provides the following relevant advantage: any marginal smoothing density (like those we are interested in) can be easily obtained from the joint density by marginalization (i.e., by discarding the particle sets and the associated Gaussian densities that refer to the instants we are not interested in). This benefit, however, is obtained at the price of a substantial computational complexity in all cases. In fact, the algorithms proposed in [14, p. 443] and [19, p. 357] require to be re-run $M$ times, if $M$ realizations of $f(\mathbf{x}_{1:T}|\mathbf{y}_{1:T})$ are needed; luckily, the processing accomplished in each run reuses all the particles and the weights computed in the forward pass. On the contrary, a single backward pass is accomplished in the algorithm derived in [13, p. 75]; this entails, however, the generation of a new set of weighted particles and Gaussian densities (representing the nonlinear state component and the linear state component, respectively); moreover, the evaluation of marginal smoothing densities is computationally intensive, since it requires merging all the information (particles, weights and Gaussian densities) emerging from both passes (see [13, Par. 4.1.2, p. 80]).

The second feature concerns the fact that the particles and the associated weights generated in the forward pass are reused in the backward pass, even if in different ways. More specifically, in the backward pass of the RPBS techniques of [13] and [19], particles are re-weighted; moreover, each new weight is evaluated as the product of the weight computed in the forward pass for the considered particle with a new weight generated on the basis of backward statistics (see, in particular, step 3)-b)-ii) of Algorithm 1 in [19, p. 357] and the particle smoothing task of Algorithm 4 in [14, p. 443]). On the one hand, the reuse, in the backward pass, of the particles generated in the forward pass greatly simplifies BIF. On the other hand, it places a strong constraint on the support of each of the pdfs computed for nonlinear state component; in fact, such a support is restricted to that identified for the predicted/filtered pdfs in the forward pass. This is the reason why the RBPS technique developed in [13] includes an algorithm for generating, in the backward pass, new particles, which are independent of those computed in the forward pass. The price to be paid for this, however, is represented by the additional computational load due to 1) particle generation in the backward pass and b) the complexity of the method employed for merging forward and backward particles (and their associated weights) to compute the required smoothed densities (see, in particular, [13, Par. 4.1.2, p. 80]).

As far as the third feature is concerned, it is worth mentioning that the use of resampling in [14], [19] is substantially different from that of [13]. In fact, in the first case, resampling is applied to the particle set generated in the TU of each recursion of the forward pass when evaluating a new trajectory in a backward pass; this is motivated by the fact that the mechanism of particle selection can benefit from more refined statistical information, since the new weights generated in the backward pass for the available particle sets are expected to be more reliable than those computed in the forward pass. On the contrary, in the second case, resampling is applied to the new particle set generated in each recursion of the backward pass, exactly like in the forward pass.

Finally, the fourth feature concerns the exploitation of Kalman techniques and, in particular, of Kalman smoothing for the linear state component in the considered RBPS algorithms. Note, however, that a different use of these standard tools is made in the considered manuscripts. In fact, on the one hand, in the RBPS techniques proposed in [13, p. 76] and [14, p. 443] smoothing for linear state component is accomplished within the backward pass and exploits the statistical information about the linear state component generated by Rao-Blackwellized filtering in the forward pass. On the other hand, in [19] the backward pass aims at generating a trajectory for the nonlinear state component only; such a trajectory is based on a) the information generated in the forward pass about this component and b) those generated about the linear state component in the backward pass only. For this reason, in this case, an additional forward pass for the linear state component only is accomplished, under the assumption that the nonlinear state trajectory is known, after that the backward pass has been completed; finally, Kalman smoothing is carried out to merge forward and backward information, as illustrated at the end of the previous Paragraph.

From the considerations illustrated above, it can be easily inferred that, on the one hand, the RBSS algorithm shares feature 4) and part of feature 2) with the other RBPS techniques (in fact, it reuses the FF particles, but not their weights). On the other hand, the RBSS algorithm does not share features 1) and 3); this makes it much faster, since both resampling and the generation of multiple trajectories are time consuming tasks. The other significant differences between the RBSS algorithm and the other methods can be summarized as follows. The algorithms developed in [13] and [14] apply to a mixed linear/nonlinear SSM whose state equation for the nonlinear component (see (1)) does contain the nonlinear term $\mathbf{f}_{l}^{(L)}(\mathbf{x}_{l}^{(N)})$ (see, in particular, [13, eq. (50), p. 75] and [14, eq. (10a), p. 441]); consequently, the only alternative method applicable to the SSM expressed by (1)-(3) in its complete form is represented by the technique devised in [19]. Moreover, as mentioned in the previous Paragraph, the RBSS algorithm, unlike all the other RBPS methods, fully exploits the available pseudo-measurements.

4.4 A message passing algorithm for estimating the joint smoothing

density

Even if backward processing in the RBSS algorithm has been explicitly devised for estimating the marginal smoothing densities $\{f(\mathbf{x}_{l}|\mathbf{y}_{1:T})\}$, the message passing procedure each of its recursion consists of can be easily modified to generate, like the RBPS method proposed in [19], $M$ (equally likely) nonlinear state trajectories providing a point mass approximation of the joint smoothing pdf $f(\mathbf{x}_{1:T}^{(N)}|\mathbf{y}_{1:T})$ (e.g., see [19, eq. 9]). In practice, this requires: a) accomplishing a single forward pass (MPF) followed by $M$ distinct backward passes; b) modifying part of the backward processing devised for RBSS. As far as the last point is concerned, let us focus, like in Paragraph  4.2, on the $(T-l)$-th recursion of the backward pass (with $l=T-1,T-2,...,1$) of the new particle smoother (called enhanced RBSS, ERBSS, in the following). The modifications made within the considered recursion originate from the fact that the nonlinear state trajectory $\{\mathbf{x}_{be,l}^{(N)},l=1,2,...,T\}$ constructed in the ERBSS backward pass consists entirely of particles generated in the forward pass (and not of a linear combination of them, like in RBSS; see (59)). For this reason, we set $\mathbf{x}_{be,l+1}^{(N)}=$ $\mathbf{x}_{(l+1)/l,j_{l+1}}^{(N)}$ and $(\mathbf{\eta}_{be,l+1}^{(L)},\mathbf{C}_{be,l+1}^{(L)})=(\mathbf{\eta}_{sm,l+1,j_{l+1}}^{(L)}$,$\mathbf{C}_{sm,l+1,j_{l+1}}^{(L)})$ in the input messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(N)})$ (18) and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l+1}^{(L)})$ (19), respectively, if the specific particle $\mathbf{x}_{(l+1)/l,j_{l+1}}^{(N)}$ has been selected within the particle set $\{\mathbf{x}_{(l+1)/l,j}^{(N)},j=0,1,...,N_{p}-1\}$ in the previous (i.e., in the $(T-l-1)$-th) recursion; the other two input messages $\vec{m}_{fp,j}(\mathbf{x}_{l}^{(N)})$ (16) and $\vec{m}_{fp,j}(\mathbf{x}_{l}^{(L)})$ (17), however, remain unchanged. ERBSS backward processing can be organized according to seven steps, exactly like RBSS. The first six steps coincide with steps 1-6 of the RBSS algorithm, whereas the remaining one is described below.

7. Sample $\mathbf{x}_{l}^{(N)}$ and  generate *input messages for the next recursion *- This requires: a) drawing a sample (denoted $\mathbf{x}_{l/(l-1),j_{l}}^{(N)}$) from the particle set $\{\mathbf{x}_{l/(l-1),j}^{(N)}\}$, whose elements are characterized by the probabilities $\{\Pr\{\mathbf{x}_{l/(l-1),j}^{(N)}\}=W_{sm,l,j}\,\}$; b) setting $\mathbf{x}_{be,l}^{(N)}=\mathbf{x}_{l/(l-1),j_{l}}^{(N)}$ and $(\mathbf{\eta}_{be,l}^{(L)},\mathbf{C}_{be,l}^{(L)})=(\mathbf{\eta}_{sm,l,,j_{l}}^{(L)}$,$\mathbf{C}_{sm,l,,j_{l}}^{(L)})$, so that the nonlinear backward trajectory is extended by one step, and the input messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(N)})$ (18) and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{l}^{(L)})$ (19) are ready for the next recursion.

The initialization of the ERBSS algorithm requires the knowledge of its input messages $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{T}^{(N)})$ and $\overset{\leftarrow}{m}_{be}(\mathbf{x}_{T}^{(L)})$, that are defined as

[TABLE]

and

[TABLE]

respectively; here, $\mathbf{x}_{T/(T-1),j_{T}}^{(N)}$ denotes the particle selected by sampling the particle set $\{\mathbf{x}_{T/(T-1),j}^{(N)}\}$; the probabilities of its particles are proportional to their weights $\{w_{T/T,j}\}$ generated by the MPF MU for the nonlinear state component in its final recursion.

As already mentioned above, the backward pass described above has to be repeated $M$ times, once for each of the $M$ nonlinear state trajectories; then, smoothing of the linear state component is accomplished for each of them. For this reason, as already explained at the end of Paragraph 4.2, the following tasks are carried out for each nonlinear state trajectory: a) a new forward pass, followed by a new backward pass, is run for the linear state component only (under the assumption that the nonlinear state component is known); b) forward prediction and backward estimation messages are merged.

It is worth stressing that the structure of the proposed ERBSS technique is very similar to that of the Algorithm 2 described in [19, p. 359]; the main differences between these two algorithms can be summarized as follows:

The backward processing developed in [19, p. 359] exploits the knowledge of the particle sets/weights generated in the forward pass, but ignores the associated Gaussian models that represent the forward predictions for the linear state component (actually, the use of such models is limited to the initialization of the backward simulator). Consequently, step 3 of our RBSS algorithm is not accomplished or, equivalently, (37) and (38) are replaced by $\mathbf{W}_{sm,l,j}^{(L)}\triangleq\mathbf{W}_{be,l,j}^{(L)}$ and $\mathbf{w}_{sm,l,j}^{(L)}\triangleq\mathbf{w}_{be,l,j}^{(L)}$ , respectively. From a conceptual viewpoint, two specific motivations can be provided for this specific choice. The first is represented by the fact that, generally speaking, the message $\vec{m}_{fp,j}(\mathbf{x}_{l}^{(L)})$ and the message $\overset{\leftarrow}{m}_{be,j}(\mathbf{x}_{l}^{(L)})$ (see (17) and (29), respectively) refer to a specific forward nonlinear trajectory and to a (unique) backward nonlinear trajectory, respectively, that *do not merge at the considered instant *(i.e., at the instant $t=l$); consequently, fusing these densities may result in poor statistical information and, in particular, may lead to the evaluation of inaccurate weights for the particle set $\{\mathbf{x}_{l/(l-1),j}^{(N)}\}$. The second motivation is represented by the fact that statistical (Gaussian) models generated by backward processing for the linear state component are really conditioned on the selected nonlinear state trajectory; for this reason, once backward processing is over, a new forward pass only has to be carried for each of the $M$ nonlinear trajectories (in other words, unlike the ERBSS technique, an additional backward pass is no more required). 2. 2.

The particle weights evaluated by Algorithm 2 of [19, p. 359] in its backward pass are partly based on the weights $\{w_{l/l,j}\}$ (computed in the forward pass). In particular, the weight $w_{l/l,j}$ replaces $w_{4,l,j}$ in the expression of the overall weights (see $W_{l,j}$ (50) and [19, Algorithm 1, step 3)-b)-ii), p. 357]) for any $j$ and $l$.

Actually, our computer simulations have evidenced that particle smoothing benefits from merging forward and backward information about the linear state component; in fact, this improves both numerical stability of BIF and its estimation accuracy through a more precise evaluation of the overall particle weights $\{W_{l,j}\}$. From a conceptual viewpoint, this choice is motivated by the fact that, as already mentioned at the beginning of Paragraph 4.2, the particle $\mathbf{x}_{l/(l-1),j}^{(N)}$ and its associated Gaussian model $\mathcal{N}(\mathbf{x}_{l}^{(L)};\mathbf{\eta}_{fp,l,j}^{(L)},\mathbf{C}_{fp,l,j}^{(L)})$ should be considered as two parts of the same hypothesis, so that they should be exploited jointly.

5 Numerical Results

In this Section MPF and the smoothing algorithms developed in this manuscript444Our simulations have evidenced that, for the considered SSM, the Algorithm 2 of [19, p. 359] suffers from ill-conditioning and that, even if its square root implementation is adopted, its computational load and accuracy are very close to that of the ERBSS technique. are compared in terms of accuracy and computational load for a specific CLG system, characterized by $D_{L}=3$, $D_{N}=1$ and $P=2$. The structure of the considered system has been inspired by the example proposed in [26] (where it is proposed as a good example for the application of MPF) and is characterized by: a) the state models

[TABLE]

and

[TABLE]

with $\mathbf{w}_{l}^{(L)}\sim\mathcal{N}(0,(\sigma_{w}^{(L)})^{2}\mathbf{I}_{3})$, $w_{l}^{(N)}\sim\mathcal{N}(0,(\sigma_{w}^{(N)})^{2}$; b) the measurement model

[TABLE]

with $\mathbf{e}_{l}\sim\mathcal{N}(0,(\sigma_{e})^{2}\mathbf{I}_{2})$. Note that the state equation (67), unlike its counterpart proposed in [26], depends on $x_{l}^{(N)}$, so that the pseudo-measurement $\mathbf{z}_{l}^{(N)}$ (15) can be evaluated for this system.

In our computer simulations our assessment of state estimation accuracy is based on the evaluation of two root mean square errors (RMSEs), one (denoted $RMSE_{N}($alg$)$, where ‘alg’ denotes the algorithm this parameter refers to) referring to the (monodimensional) nonlinear state component, the other one (denoted $RMSE_{L}($alg$)$) to the (three-dimensional) linear state component; note, however, that the last parameter represents the square root of the average mean square error (MSE) evaluated for the three elements of $\mathbf{x}_{l}^{(L)}$. Our assessment of computational requirements is based, instead, on assessing the average computation time for processing a single block of measurements (this quantity is denoted CTB in the following). Moreover, in our computer simulations, the following choices have been always made: a) $T=200$ has been selected for the length of the observation interval; b) $M=N_{p}$ has been chosen for the EBRSS ($M\lesssim N_{p}$ is recommended in [15]).

Some results illustrating a) the dependence of $RMSE_{L}$ and $RMSE_{N}$ (CTB) on the number of particles ($N_{p}$) for the MPF, the RBSS and ERBSS algorithms are illustrated in Fig. 4 (Fig. 5) 555In these and in the following figures simulation results are identified by markers, whereas continuous lines are drawn to ease reading.; in this case $\sigma_{w}^{(L)}=\sigma_{w}^{(N)}=2\cdot 10^{-1}$ and $\sigma_{e}=3\cdot 10^{-2}$ have been selected. From these results the following conclusions can be easily inferred for the considered scenario:

On the one hand, a negligible improvement in the estimation accuracy of all the considered algorithms is achieved for $N_{p}\geq 100$ (actually, a similar result has been found for other values of $\sigma_{e}$, $\sigma_{w}^{(L)}$ and $\sigma_{w}^{(N)}$); for this reason, $N_{p}$ $=100$ has been selected in all the computer simulations the following results refer to. 2. 2.

The RBSS algorithm outperforms MPF by about $21.12\%$ ($36.5\%$) in terms of $RMSE_{L}$ ($RMSE_{N}$) for $N_{p}=100$. A negligible improvement in RBSS accuracy can be obtained by accomplishing a further smoothing for the linear state component (as explained at the end of Paragraph 4.2); this reason, this possibility is no more considered in the following. Note also that the RBSS improvement is obtained at the price of a limited computational cost, since its CTB is about twice that of MPF. 3. 3.

The ERBBS algorithm provides a by far richer statistical information than the RBSS algorithm, but achieves slightly better accuracy in state estimation and entails a substantially larger computational load, even for small values of $N_{p}$ (for instance, the ERBBS computation time is about 100 times larger than that of RBBS for $N_{p}$ $=100$). Note also that the CTB gap between the EBRSS algorithm and both the RBSS and the MPF techniques becomes larger as $N_{p}$ increases. For this reason, the ERBBS is not taken into consideration anymore in the following simulations. 4. 4.

A relevant gap between $RMSE_{L}($MPF$)$ and $RMSE_{N}($MPF$)$ ($RMSE_{L}($RBSS$)$ and $RMSE_{N}($RBSS$)$) exists; unluckily, the RBSS algorithm is unable to reduce this gap. This can be related to the fact that smoothing accuracy is significantly influenced by that achieved in the forward pass.

A comparison between the MPF and the RBSS state estimation errors has also evidenced that the RMSE improvement provided by the latter algorithm is mainly related to its ‘peak shaving’ effect. In fact, the amplitude of the spikes appearing in the state estimation error at the end of the forward pass are substantially reduced by smoothing. Note, however, that the elements of the system state do not necessarily benefit from this effect in the same way; for instance, for our specific SSM, this effect is stronger for the nonlinear state component than for each of the three elements of the linear state component.

In our work the dependence of $RMSE_{L}$ and $RMSE_{N}$ on the intensity of the process noise and on that of the measurement noise has been also analysed. Some results illustrating the dependence of $RMSE_{L}$ and $RMSE_{N}$ on $\sigma_{e}$ (under the assumption that $\sigma_{w}^{(L)}=$ $\sigma_{w}^{(N)}=2\cdot 10^{-2}$)  are shown in Fig. 6. From these results it is easily inferred that the performance gap between MPF and RBSS shrinks as $\sigma_{e}$ increases; this is due to the fact that a stronger measurement noise results in a poorer quality of the statistical information generated in the forward pass, and this impairs more and more the RBSS estimation process. Other simulation results (not shown here for space limitations) have also evidenced that, for a given intensity of the measurement noise, the gap between $RMSE_{L}($MPF$)$ and $RMSE_{L}($RBSS$)$ (and, similarly, between $RMSE_{N}($MPF$)$ and $RMSE_{N}($RBSS$)$) remains stable as $\sigma_{w}=\sigma_{w}^{(L)}=\sigma_{w}^{(N)}$ changes (in particular, $\sigma_{w}=\in[10^{-2},2\cdot 10^{-1}]$ has been assumed in our simulations).

6 Conclusions

In this manuscript the smoothing problem for SSMs has been analysed from a FG perspective. This has allowed us to devise new RBPS methods for CLG SSMs. Computer simulations for a specific SSM evidence that the RBSS algorithm achieves a good performance-complexity tradeoff. Our future work concerns the application of FG methods to the problems of filtering and smoothing for other classes of SSMs.

Acknowledgment

We would like to thank Dr. Fredrik Lindsten (Uppsala University, Department of Information Technology) for his constructive comments.