Rao-Blackwellized Particle Smoothing as Message Passing

Rao-Blackwellized Particle Smoothing as Message Passing

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.

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 evaluated at point is denoted ; represents the pdf of a Gaussian random vector characterized by the mean and covariance matrix evaluated at point ; the precision (or weight) matrix associated with the covariance matrix is denoted , whereas the transformed mean vector is denoted .

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 -th interval is represented by the -dimensional real vector ; this is partitioned in a) its -dimensional linear component and b) its -dimensional nonlinear component (with and ). The update equations of the linear and nonlinear components are given by

 x(L)l+1=A(L)l(x(N)l)x(L)l+f(L)l(x(N)l)+w(L)l, (1)

and

 x(N)l+1=f(N)l(x(N)l)+A(N)l(x(N)l)x(L)l+w(N)l, (2)

respectively; here, () is a time-varying -dimensional (-dimensional) real function, () is a time-varying () real matrix and () is the -th element of the process noise sequence (), which consists of - dimensional (-dimensional) independent and identically distributed (iid) noise vectors (statistical independence between and is also assumed for simplicity). Moreover, in the -th interval some noisy observations, collected in the measurement vector

 yl≜[y0,l,y1,l,...,yP−1,l]T=hl(x(N)l)+Bl(x(N)l)x(L)l+el, (3)

are available about ; here, is a time-varying real matrix, is a time-varying -dimensional real function and the -th element of the measurement noise sequence consisting of -dimensional iid noise vectors and independent of both and . 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 (where represents the length of the observation interval), given a) the initial pdf and b) the -dimensional measurement vector .

3 A FG-Based Representation of the Smoothing Problem

In this Section we formulate the computation of the marginal smoothed density (with ) 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 and the observation model ; 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 in place of the posterior pdf . 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])

 f(xl,y1:T) =f(yl:T∣∣xl,y1:(l−1))f(xl,y1:(l−1)) =f(yl:T|xl)f(xl,y1:(l−1)) (4)

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 and the conditional pdf , respectively, for any .

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 for its initialization. In the MU step of its -th recursion (with ) the joint pdf

 f(xl,y1:l)=f(xl,y1:(l−1))f(yl|xl) (5)

is computed on the basis of pdf , and the new measurement vector . In the TU step, instead, the pdf (5) is exploited to compute the pdf

 (6)

representing a prediction about the future state .

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

 f(y(l+1):T|xl)=∫f(y(l+1):T|xl+1)f(xl+1|xl)dxl (7)

and

 f(yl:T|xl)=f(y(l+1):T|xl)f(yl|xl), (8)

respectively. Note that this procedure requires the knowledge of the pdf 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.

 →mfp(xl)≜f(xl,y1:(l−1)), (9)

the forward messages emerging from the equality node and that passed along the edge associated with are given by and , 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

 ←mbe(xl+1)≜f(y(l+1):T|xl+1) (10)

enters the FG along the half edge associated with (the flow of backward messages is indicated by blue arrows in Fig. 1), the backward message  emerging from the node associated with the pdf is given by (see (7))

 ←mbp(xl) =∫←mbe(xl)f(xl+1|xl)dxl =∫f(y(l+1):T|xl+1)f(xl+1|xl)dxl =f(y(l+1):T|xl). (11)

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

 f(yl|xl)←mbp(xl)=f(yl|xl)f(y(l+1):T|xl) =f(yl:T|xl)=←mbe(xl) (12)

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 can be evaluated as (see (4), (9) and (12))

 f(xl,y1:T)=→mfp(xl)←mbe(xl), (13)

with (note that and )

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 from that for , i.e. the evaluation of   from that of . In practice, from a graphical viewpoint, two sub-graphs, one referring to smoothing for , the other one to smoothing for , are developed first; then, they are merged by adding five distinct equality nodes, associated with the variables (namely, , , , and ) 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 can benefit not only from the measurement , but also from the so called pseudo-measurement (see (2))

 (14)

which, from a statistical viewpoint, is characterized by the pdf . Similarly, the pseudo-measurement (see (1))

 (15)

characterized by the pdf , can be exploited in smoothing for the nonlinear component .  These considerations explain why the upper (lower) sub-graph shown in Fig. 2 contains an additional node representing the pdf () and a specific node not referring to the above mentioned pdf factorizations, but representing the transformation from the couple to ( to ); 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 and , 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 and on the basis of the available a priori information and the noisy measurement . 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 .

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 () is Gaussian and all its elements have zero mean and covariance () for any ; b) the measurement noise is Gaussian having zero mean and covariance matrix for any ; c) all the above mentioned Gaussian processes are statistically independent. Under these assumptions, the pdfs , and are Gaussian with mean (covariance matrix) , and , respectively (, and , respectively). Similarly, the pdfs and are Gaussian with mean (covariance matrix) and , respectively ( and , 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 -th recursion (with ), the particle set , consisting of distinct particles, is predicted for the nonlinear state component (TU for this component); the weight assigned to the particle is equal to for any , since the use of particle resampling in each recursion is assumed. The particle weights are updated in the MU of the following (i.e., -th) recursion on the basis of the new measurement (MU for the nonlinear component): the new weights are denoted in the following and, generally speaking, are all different. This is followed by particle resampling, that generates the new particle set (usually containing multiple copies of the most likely particles of the set ). 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 in the -th recursion (TU for the linear state component) and associated with is denoted . Note that only a portion of these Gaussian models is usually updated in the MU of the next (i.e., -th) recursion; in fact, this task follows particle resampling, which typically leads to discarding a fraction of the particles collected in the set .

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 -th recursion of BIF and smoothing (with ) 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 , that represents the index of the particle (namely, the particle ) representing 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 and , that are predicted the -th recursion of the forward pass; the second one, instead, is made of the messages and , that are generated in the -th recursion of the backward pass. The messages of the first group are defined as

 →mfp,j(x(N)l)≜δ(x(N)l−x(N)l/(l−1),j) (16)

and

 →mfp,j(x(L)l)≜N(x(L)l;η(L)fp,l,j,C(L)fp,l,j), (17)

and can be interpreted as the -th hypothesis about a) the value (namely, ) taken on by the (hidden) nonlinear state component and b) the statistical representation of the (hidden) linear state component associated with such a value, respectively. In the -th recursion of FF, the likelihood of this hypothesis is assessed by evaluating the above mentioned weight ; 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 and additional (backward) information about (see (18) and (19) below).

The input messages of the second group are defined as

 ←mbe(x(N)l+1)≜δ(x(N)l+1−x(N)be,l+1) (18)

and

 ←mbe(x(L)l+1)≜N(x(L)l+1;η(L)be,l+1,C(L)be,l+1), (19)

and represent part of the statistical information generated in the previous (i.e., the -th) recursion of the backward pass. In particular, as explained in detail below, the messages and convey the final estimate (i.e., a single particle representation) of and a simplified statistical representation of , respectively. This explains why the RBSS, in the -th recursion of its backward pass, processes the input messages (16)-(19) to compute an estimate, denoted , of and a simplified statistical model, denoted , for ; these information are conveyed by the output messages and , 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 (42)  and, in particular, of its parameters (44) and (45); mathematical details about this can be found in [20, Sec. 6]). Such formulas refer to the computation of the message (emerging from an equality node fed by the messages and ) and

 mout(x2)=∫min(x1)f(x1,x2)dx1 (20)

(emerging from a function node fed by the message ), 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 and the pdfs of its components; b) the output messages and .

1. Time update for - Compute the message (see T2-5, (16) and (19))

 ←m1,j(x(L)l) =∫∫f(x(L)l+1∣∣x(L)l,x(N)l) ⋅←mbe(x(L)l+1)→mfp,j(x(N)l)dx(L)l+1dx(N)l =N(x(L)l;η(L)1,l,j,C(L)1,l,j), (21)

where

 w(L)1,l,j ≜W(L)1,l,jη(L)1,l,j=(A(L)l,j)TW(L)w ⋅[¯Cl+1w(L)be,l+1−P(L)lf(L)l,j], (22)
 W(L)1,l,j≜(C(L)1,l,j)−1=(A(L)l,j)TW(L)wP(L)lA(L)l,j, (23)

, , , , , and .

2. Measurement update for - Compute: a) the message

 →mj(z(L)l)=f(z(L)l∣∣x(N)l/(l−1),j,~x(N)l+1)=δ(z(L)l−z(L)l,j), (24)

where and ; b) the messages (see T2-3, T1-3, T2-2 and T1-2, respectively; see also (16), (21) and (24))

 ←m2,j(x(L)l) =∫∫f(z(L)l∣∣x(L)l,x(N)l) =N(z(L)l,j;A(N)l,jx(L)l,C(N)w), (25)
 ←m3,j(x(L)l) =←m1,j(x(L)l)←m2,j(x(L)l) =N(x(L)l;η(L)3,l,j,C(L)3,l,j), (26)
 ←m4,j(x(L)l) =∫f(yl∣∣x(N)l,x(L)l)→mfp,j(x(N)l)dx(N)l =N(yl;Bl,jx(L)l+hl,j,Ce) (27) ≡N(x(L)l;η(L)4,l,j,C(L)4,l,j) (28)

and

 ←mbe,j(x(L)l) =←m3,j(x(L)l)←m4,j(x(L)l) =N(x(L)l;η(L)be,l,j,C(L)be,l,j). (29)

Here,

 w(L)3,l,j≜W(L)3,l,jη(L)3,l,j=w(L)1,l,j+(A(N)l,j)TW(N)wz(L)l,j, (30)
 W(L)3,l,j≜(C(L)3,l,j)−1=W(L)1,l,j+(A(N)l,j)TW(N)wA(N)l,j, (31)

, ,

 w(L)4,l,j≜W(L)4,l,jη(L)4,l,j=(Bl,j)TWe(yl−hl,j)