Particle Smoothing for Conditionally Linear Gaussian Models as Message Passing Over Factor Graphs

In this paper, 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 a forward–backward message passing over a factor graph, we focus on the above-mentioned class of models and derive two novel particle smoothers for it. Both the proposed techniques are based on the well-known two-filter smoothing approach and employ marginalized particle filtering in their forward pass. However, on the one hand, the first smoothing technique can only be employed to improve the accuracy of state estimates with respect to that achieved by forward filtering. On the other hand, the second method, that belongs to the class of Rao–Blackwellized particle smoothers, also provides a point mass approximation of the so-called joint smoothing distribution. Finally, our smoothing algorithms are compared, in terms of estimation accuracy and computational requirements, with a Rao–Blackwellized particle smoother recently proposed by Lindsten et al.  (“Rao–Blackwellized particle smoothers for conditionally linear Gaussian models,”  IEEE J. Sel. Topics Signal Process., vol. 10, no. 2, pp. 353–365, 2016).


I. 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 a 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 twofilter 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 sample 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], [18] and Rao-Blackwellized particle smoothing (RBPS) [13], [14], [20], respectively.
Recently, the filtering problem for CLG SSMs has been investigated from a factor graph (FG) perspective in [21], 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 fixedinterval smoothing problem [13] and to develop novel solutions for it.The proposed solutions are represented by two novel particle smoothing methods, the first one dubbed serial particle smoothing (SPS), the second one Rao-Blackwellized serial smoothing (RBSS).These methods share the following relevant features: a) they are based on the two-filter smoothing formula and employ MPF in their forward pass; b) they can be derived applying the well known sum-product algorithm (SPA) [23], [24], together with a specific scheduling procedure, to the same FG developed in [21] and [22] for a CLG SSM; c) unlike the RBPS methods devised in [13] and [14], they can be employed for a SSM in which both the linear and nonlinear state components influence each other.Moreover, the computational requirements of the SPS technique are comparable with that of MPF, since a single estimate of the nonlinear state trajectory is generated in its backward pass; consequently, it is substantially simpler than the RBPS methods illustrated in [13], [14] and [20].On the contrary, the computational load of the RBPS method is significantly higher, being close to that of the method developed in [20]; in fact, these methods employ similar procedures to generate a joint smoothing distribution over the entire observation interval.
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 [23], [25], [26], 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 II.A representation of the smoothing problem through Forney-style FGs for both an arbitrary SSM and a CLG SSM is provided in Section III.In Section IV the SPS and the RBSS techniques are developed for a CLG SSM, and are compared with that of other particle smoothing methods.Section V is devoted to comparing, in terms of accuracy and computational effort, our FG-based smoothing algorithms with the RBPS method developed in [20].Finally, some conclusions are offered in Section VI.
Notations: The probability density function (pdf) of a random vector R evaluated at point r is denoted f (r); N (r; η r , C r ) represents the pdf of a Gaussian random vector R characterized by the mean η r and covariance matrix C r evaluated at point r; the precision (or weight) matrix associated with the covariance matrix C r is denoted W r , whereas the transformed mean vector W r η r is denoted w r .

II. SYSTEM MODEL
In the following we focus on the discrete-time CLG SSM described in [21], [22].In brief, the SSM hidden state in the l-th interval is represented by the D-dimensional The update equations of the linear and nonlinear components are given by respectively; here, f ) is the l-th element of the process noise sequence {w k } is also assumed for simplicity).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 (x l |y 1:T ), l = 1, 2, ..., T } (where T represents the length of the observation interval), given a) the initial pdf f (x 1 ) and b) the where with l = 1, 2, ..., T , is a P -dimensional vector collecting the noisy observations available about x l in the l-th interval.Here, ) is a time-varying P -dimensional real function and e l the l-th element of the measurement noise sequence {e k } consisting of P -dimensional iid noise vectors (all independent of both {w

III. A FG-BASED REPRESENTATION OF SMOOTHING
In this Section we formulate the computation of the marginal smoothed density f (x l |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 (x l+1 |x l ) and the observation model f (y l |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 (x l , y 1:T ) in place of the posterior pdf f (x l |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 [23, Sec. II, p. 1297]).Note that the validity of statement a) relies on the following mathematical results: 1) the factorization (e.g., see [8,Sec. 3]) for the pdf of interest; 2) 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 (x l , y 1:(l−1) ) and the conditional pdf f (y l:T |x l ), respectively, for any l.
As far as FF is concerned, the formulation illustrated in [21, 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 (x 1 ) for its initialization.In the MU step of its l-th recursion (with l = 1, 2, ..., T ) the joint pdf Fig. 1: Graphical representation of the message passing for the evaluation of the joint pdf f (x l+1 , y 1:l ) and of the conditional pdf f (y l:T |x l ) on the basis of eqs.( 5)-( 6) and ( 7)-( 8), respectively (the forward and backward message flows are indicated by red and blue arrows, respectively) is computed on the basis of pdf f (x l , y 1:(l−1) ), and the new measurement vector y l .In the TU step, instead, the pdf f (x l , y 1:l ) (5) is exploited to compute the pdf representing a prediction about the future state 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 respectively.Note that this procedure requires the knowledge of the pdf f (y T |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 [23, Sect.II], such recursions can be interpreted as specific instances of the SPA1 applied to the cycle free FG of Fig. 1 (where the simplified notation of [23] 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 message2 the forward message emerging from the equality node and that passed along the edge associated with x l+1 are given by m f e (x l ) = f (x l , y 1:l ) and f (x l+1 , y 1:l ) = m f p (x l+1 ), respectively [21], [22].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 ← m be (x l+1 ) f y (l+1):T |x l+1 (10) enters the FG along the half edge associated with x l+1 (the flow of backward messages is indicated by blue arrows in Fig. 1), the backward message ← m bp (x l ) (emerging from the node associated with the pdf f (x l+1 |x l )) is given by (see (7)) Therefore, the backward message emerging from the equality node can be evaluated as (see ( 8) and ( 10)) 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 (x l , y 1:T ) can be evaluated as (see ( 4), ( 9) and ( 12)) with l = 1, 2, ..., T (note that ← m be (x T ) = 1 and m f p (x 1 ) = f (x 1 )) or, alternatively, as It is worth noting that, despite the conceptual simplicity of the procedure illustrated above, its implementation can represent a formidable task, mainly because of the multidimensional integration required in ( 6) and ( 7).This has motivated the development of the so called Rao-Blackwellization approach, according to which the state vector x l is partitioned in a nonlinear component x l ; moreover, from a statistical viewpoint, x (N ) l is represented through a set of weighted particles, whereas x (L) l through a set of a particle-dependent Gaussian pdfs (i.e., through a Gaussian mixture, GM).In [21] and [22] it has been shown that Rao-Blackwellized filtering can be seen as an instance of the SPA applied to the FG we develop for case C.2; the new FG 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 [21].For this reason, two distinct subgraphs are drawn, one referring to smoothing for x can benefit not only from the measurement y l , but also from the vector (see ( 2)) ) which, from a statistical viewpoint, is characterized by the pdf f (z ).Similarly, the vector (see (1)) ) characterized by the pdf f (z ), can be exploited in smoothing for the nonlinear component x ).This leads to the overall 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.This graphical model deserves the following important comments: • No approximation is made in deriving it.
• Its upper (lower) sub-graph contains an additional node representing the pseudo-measurement pdf f (z )) and a specific node not referring to a pdf factorization, but representing the transformation from the couple (x ); the last peculiarity is evidenced by the presence of an arrow on all the edges connected to such a node.
• Unlike the FG represented in Fig. 1, it is not cycle-free.
This property is due to the fact that, generally speaking, smoothing for x (N ) l is not decouplable from that for x (L) l ; in other words, uncertainty in each component of system state needs to be accounted for in the distribution of the other state component.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 (marginal) messages l+1 ), we are interested in deriving a BIF algorithm3 based on this FG and generating the output backward (marginal) messages ← m be (x ) on the basis of the available a priori information and the noisy measurement y l .Note that, given the FG of of Fig. 2, the evaluation of the message ← m be (x l ).It is well known that, on the one hand, in cycle-free graphical models exact marginalization can be accomplished by means of the SPA; on the other hand, the application of the SPA to a graphical model containing cycles unavoidably leads to approximate solutions [24], whatever message scheduling is adopted.Despite this, we believe that, like in the related problem of filtering for CLG SSMs [21], [22], this approach can lead to the development of accurate and computationally efficient smoothing algorithms, as shown in the following Section.

IV. PARTICLE SMOOTHING AS MESSAGE PASSING
In this Section we first illustrate some assumptions about the statistical properties of the SSM described in Section II.Then, we develop the SPS and RBSS techniques.Finally, we compare the most relevant features of these techniques with those of the other RBPS algorithms available in the technical literature.

A. Statistical properties of the considered SSM
Even if the graphical model 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 [21], to the specific class of CLG SSMs.For this reason, following [21], [22] we assume that: a) the process noise {w k } is Gaussian having zero mean and covariance matrix C e for any l; c) all the above mentioned Gaussian processes are statistically independent.Under these assumptions, the pdfs f (y l |x w , respectively).

B. Derivation of the particle serial smoother
In developing our first particle smoothing algorithm the following fundamental requirements have been set to limit its computational complexity as much as possible: 1) MPF is employed in its (single) forward pass 4 ; 2) the statistical information generated by the BIF technique adopted in its (single) backward pass can be easily combined with that provided by MPF to generate the marginal smoothed densities of interest; c) the estimates of both the linear component and the nonlinear component are available at the end of the backward pass and each of them is represented by a single trajectory.
As far as the use of MPF is concerned, the following notation is adopted here: • The j-th particle predicted for x (N ) l in the (l − 1)th recursion (with l = 2, 3, ..., T ) of this algorithm is denoted x (N ) l/(l−1),j (with j = 0, 1, ..., N p − 1, where N p is the overall number of particles); moreover, the weight assigned to this particle is denoted w l/(l−1),j (this weight is equal to 1/N p for any j, since the use of particle resampling in each recursion is assumed).
• The weight computed for the particle x (N ) l/(l−1),j in the following (i.e., in the l-th) recursion on the basis of the new measurement y l is denoted w l/l,j (with j = 0, 1, ..., N p − 1).
• The j-th particle available after particle resampling based on the weights {w l/l,j } is denoted x l/l,j (note that the set {x (N ) l/l,j } usually contains multiple copies of the most likely particles of the set {x (N ) l/(l−1),j }).
• The Gaussian model predicted for x (L) l in the (l − 1)th recursion and associated with x f p,l,j ); note, however, that only a portion of these Gaussian models is usually updated in the next (i.e., l-th) recursion on the basis of y l ; in fact, this task follows particle resampling, which typically leads to discarding a fraction of the particles collected in the set {x In the forward pass of the proposed SPS method, the statistical information generated by MPF and saved for further processing is conveyed by the forward messages and with j = 0, 1, ..., N p − 1.Let us now focus on the backward pass of the SPS algorithm and, in particular, tackle the problems of a) developing a recursive BIF algorithm based on the FG of Fig. 2 and b) merging the statistical information it generates with those computed in the forward pass (smoothing).The recursive algorithm we propose results from the application of the SPA to the FG shown in Fig. 2 and is based on the message scheduling illustrated in Figs.
3-(a) and 3-(b), both referring to the (T − l)-th recursion (with l = T − 1, T − 2, ..., 1) and the j-th particle5 x (N ) l/(l−1),j (with j = 0, 1, ..., N p − 1).The processing accomplished by the SPS algorithm within the considered recursion can be divided in three parts, that are executed serially (and this motivates the presence of the adjective 'serial' in the name of the devised smoothing algorithms).The first part, that involves the computation of the messages illustrated in Fig. 3 Note that the last message is expected to provide a more refined statistical representation of x l ) alone; consequently, it can be profitably exploited in the evaluation of new particle weights.This task is accomplished in the second part (see Fig. 3-(b)), which mainly concerns the nonlinear state component; in fact, in this part the message 6 m be,j (x ), conveying a new weight for the particle x (N ) l/(l−1),j , is computed on the basis of the available measurements and pseudo-measurements, and of the message m sm,j (x (L) l ).After running the first two parts for all the available particles, the sets of smoothed messages {m sm,j (x (L) l )} and {m sm,j (x (N ) l )} are available; in the third (and last) part of the SPS algorithm, these messages are fused to generate an estimate of the marginal smoothed pdf of x l .From Figs. 3-(a) and 3-(b) it can be also inferred that the considered recursion is fed not only by the forward messages m f p,j (x (N ) l ) (17) and m f p,j (x (L) l ) (18), but also by the particle-independent backward messages and These messages have been generated in the previous (i.e., the (T − l − 1)-th) recursion and provide statistical information about the backward estimates of x within the considered recursion. 6The backward and the smoothed estimates coincide in this case since all the message { m f p,j (x )} convey a unit weight (see ( 17)).
respectively.Consequently, in the (T − l)-th recursion of the SPS algorithm, the messages ← m be (x ) and ← m be (x (L) l ) must be also computed (see Fig. 3-(b)), so that they become available to the next backward recursion; in practice, they are generated merging the statistical information conveyed by the sets {m be,j (x (L) l )} and {m be,j (x (N ) l )}, respectively.Let us analyse now the three processing steps accomplished in the first part of the SPS technique (see Fig. 3-(a)).The first step aims at evaluating the Gaussian message representing a one-step backward prediction of the pdf of x  3

, p.1304] (with
and and f Moreover, this message can be put in the equivalent Gaussian form where the precision matrix W and respectively; here, A l/(l−1),j ) and where B l,j B l (x l/(l−1),j ) and h l,j h l (x l/(l−1),j ); this message can be also put in the equivalent Gaussian form where the precision matrix W and w respectively, and can be computed by applying twice eqs.(II.2) and (II.4) of [23, Table 2, p.1303] (referring to the computation of the backward messages emerging from an equality node).This produces the precision matrix The processing accomplished in the second part of the SPS technique can be divided in two steps (see Fig. 3-(b)).In the first step, a new weight is computed for the j-th particle on the basis of the following information: a) the backward estimate x ; c) the measurement y l .For this reason, the new overall weight W l,j for the j-th particle is expressed as a product of three factors (see the lower part of the FG shown in Fig. 3-(b)), i.e. as where the weights w 1,l,j , w 2,l,j and w 4,l,j are related to x (N ) be,l+1 , z (N ) l and y l , respectively, and represent the messages ), respectively.In practice, the weight w 1,l,j is evaluated as where K 1,l,j = (2π det(C x T Wx denotes the square of the norm of the vector x with respect to the positive definite matrix W.
The computation of the weight w 2,l,j requires the knowledge of the Gaussian message that expresses the statistical knowledge acquired about z and respectively.Given the message ← m j (z ) (44), w 2,l,j is computed by correlating it with the pdf f (z l/(l−1),j ) (expressing our prior knowledge about z (N ) l for the considered particle), i.e. as ) (44) in (47) yields, after some manipulation, where 2,l,j = w z,l,j .Finally, the weight w 4,l,j is evaluated as where K 4,l,j (2π det(C B l,j η 4,l,j ) −1 .Substituting the right-hand side (RHS) of ( 41), ( 48) and ( 52) in (40) yields8 where and It worth noting that: a) for a given SSM and under specific operating conditions, the three weights w 1,l,j , w 2,l,j and w 4,l,j may not have the same importance; b) the weight w 4,l,j can be seen as the backward counterpart of the weight w l/l,j evaluated in FF.For these reasons, in some cases, the computational load due to the evaluation of the overall weight W l,j (40) might be substantially reduced by discarding one of the weigths {w p,l,j , p = 1, 2, 4} and/or adopting the approximation w 4,l,j w l/l,j with a limited impact on estimation accuracy.
Once the whole set of weights {W l,j } is available, its normalization must be accomplished; this produces the final (normalised) weight which is conveyed by the messages (see Fig. 3-(b)) for j = 0, 1, , ..., N p − 1.
In the second step the input messages and for the next recursion are generated on the basis of the particle weights {W sm,l,j } (see (58)) and the particle-dependent Gaussian messages { ← m be,j (x )).The vector x (N ) be,l , appearing in the RHS of (60), represents the l-th element of the estimated nonlinear state trajectory and is evaluated as a weighted sum of the particles {x (N ) l/(l−1),j } , i.e. as RBPS algorithms available in the literature, namely with the particle smoothers proposed by Briers et al. in [13, Sec. 4, p. 75], by Fong et al. [14] and by Lindsten et al. [20] (these algorithms are denoted Alg-B, Alg-F and Alg-L, respectively, to ease reading of the remaining part of this manuscript).Despite their substantially different structures, the last three algorithms share the following relevant features: 1) the computation of an estimate of the joint smoothing density f (x 1:T |y 1:T ); 2) the reuse of the 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 SPS/RBSS algorithms can be easily identified.
The first feature refers to the fact that these techniques aim at generating realizations from the complete joint smoothing pdf f (x 1:T |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) in Alg-B and Alg-L, or b) a trajectory for the entire state in Alg-F (since 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, backward filtering and smoothing in Alg-F and Alg-L require to be re-run M times, if M realizations of f (x 1:T |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 Alg-B ; in this pass, however, the generation of a new set of weighted particles and Gaussian densities (representing the nonlinear state component and the linear state component, respectively) is carried out.
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 Alg-F and Alg-L, 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, the particle smoothing task of Algorithm 4 in [14, p. 443] and step 3)-b)-ii) of Algorithm 1 in [20, p. 357]).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 9 .This is the reason why Alg-B includes an algorithm for generating, in its 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 2) the complexity of the processing required for fusing forward and backward particles (and their associated weights) in the computation of 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 Alg-F and Alg-L is substantially different from that of Alg-B.In fact, in the first case, resampling is applied to the particle set produced in the TU of each recursion of the forward pass when generating 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 set 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 KS 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 Alg-B and Alg-F, 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 Alg-L the backward pass aims at generating 9 Note that this occurs in the SPS and in the RBSS techniques too.
a trajectory for the nonlinear state component 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, KS is carried out to merge forward and backward information, as illustrated at the end of the previous Paragraph.
Let us now focus on the similarity of the proposed SPS technique with Alg-B, Alg-F and Alg-L in terms of the four features illustrated above.The SPS 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 (unless the weight w 4,l,j (52) is replaced by the MPF weight w l/l,j , as suggested in Paragraph IV-B).It does not share, however, features 1) and 3); this makes it much faster than Alg-B, Alg-F and Alg-L in the computation of marginal smoothed densities, since both resampling and the generation of multiple trajectories or new particles in the backward pass are time consuming tasks.In fact, the computational complexity associated with the computation of marginal smoothed densities is of order O(N p T M ) for Alg-F and Alg-L, of order O(N 2 p T ) for Alg-B10 ; on the contrary, the computational complexity for the SPS technique is of order O(N p T ) only.
The RBSS technique deserves different comments, since it is structurally very similar to Alg-F and, consequently, its complexity is also of order O(N p T M ).In fact, the main difference between these two algorithms is related to the computation of particle weights in the backward pass.More specifically, in Alg-F, the evaluation of the weight for the j-th particle is not based on (40) (and, consequently, on (55)-(57)), but on the formula (see step b)-ii) of Algorithm 1 in [20, p. 357]) Here, w j l corresponds to w l/l,j (i.e., it represents the j-th particle weight computed in the MPF MU for the nonlinear state component and, consequently, depends on y l ), whereas the remaining factors appearing in the RHS of the last equation are related to the statistical information about x (L) l generated in both the forward pass and the backward pass, and to the pseudo-measurements.More specifically, it can be shown that (analytical details are omitted here for space limitations): • The product det(Λ j t ) −1/2 exp(−η j t /2) appearing in the RHS of (69) quantifies the correlation between the Gaussian pdfs representing x l/(l−1),j , in the forward pass and in the backward pass (η j t is expressed by the sum of three terms, a scalar product and two squared non Euclidean norms; see eq. (21a) of [20, p. 357]).Note, however, that these Gaussian pdfs are not expressed by m f p,j (x (L) l ) (18) and ← m be (x (L) l ) (61).In fact, on the hand, the Gaussian 1) A small improvement in the estimation accuracy of all the considered algorithms is achieved for N p > 200.
2) The MPF and the SPS algorithm closely approach KF and KS accuracy, respectively, as N p increases.3) Even if the RBSS algorithm and Alg-L provide by far richer statistical information than the SPS algorithm, they achieve slightly better accuracy in state estimation for N p < 300, despite their substantially larger computational load.Moreover, for N p > 300 the SPS algorithm performs slightly better than the RBSS algorithm and Alg-L if M = 100 is set; this is due to the fact the overall number of trajectories generated by the last two algorithms is not allowed to increase with N p .In fact, a further improvement in their accuracy can be achieved by increasing M proportionally to N p ; this is evidenced by the trend of the RM SE N curve provided for Alg-L with M = N p .Unluckily, this result is achieved at the price of a substantial increase in the computational load.For instance, the CTB of the Alg-L with M = N p is roughly 2, 3, 4 and 5 times larger than that of the Alg-L with M = 100 for N p = 200, 300, 400 and 500, respectively.Further numerical results (not shown here for space limitations) have also evidenced that: a) as N p increases, the particlebased representation f (p l |y 1:N ) (see (66)) of the pdf of x (N ) l closely approaches the Gaussian pdf computed by KF; b) the RMSE improvement provided by the considered smoothing algorithms over MPF is mainly related to its 'peak shaving' effect in state estimation error (in other words, the amplitude of the spikes appearing in the state estimation error at the end of the forward pass are substantially reduced by smoothing; c) a minor degradation in the estimation accuracy of the SPS and RBSS algorithms is found if w 4,l,j = w l/l,j and w 2,l,j = 1 are set in formula (40) to reduce the computational effort 13required by the evaluation of the particle weights {W l,j } (consequently, only the weight w 1,l,j (41) has to be computed in BIF; note that this requires computing the inverse of the matrix C (N ) 1,l,j (43) only).As far the last point is concerned, it is also worth stressing that the little relevance of the weight w 2,l,j is due to the fact that the correlation appearing in the RHS of (47) exhibits a weak dependence on the selected particle in the presence of strong process noise for the linear state component (i.e., when σ v is large in this case).
All the conclusions illustrated above are also supported by the performance results in Fig. 5, that shows the dependence of RM SE N (blue lines) and RM SE L (red lines) on the number of particles (N p ) for SSM #2.In particular, it is interesting to note that, even in this case, all the RM SE L curves are flat, whereas the RM SE N curves achieve a floor as N p increases.Moreover, the accuracy gap between the SPS algorithm and the RBSS algorithm (or Alg-L) is really small, despite the fact they have substantially different computational requirements.This can be easily inferred from Fig. 6, that shows the dependence of the CT B on the number of particles (N p ) for all the considered algorithms.In particular, these results show that:  1) The accuracy improvement of SPS over MPF is obtained at the price of a limited increase in computational complexity (the SPS and the MPF algorithms have similar computational requirements); 2) The CTB gap between Alg-L (or the RBSS algorithm) and the SPS algorithm is significant (for instance, the computation time of Alg-L is about 50 times larger than that of SPS for N p = 100).
3) The RBSS and its simplified version (SRBSS) have larger computational requirements than Alg-L for the considered SSM, despite the fact that the evaluation of the weights for SRBSS is based on a substantially simpler formula than that employed for Alg-L.This gap can be related to the slightly different processing accomplished for the linear state component (see the last part of Paragraph IV-D).In our work the dependence of RM SE L and RM SE N on the intensity of the measurement noise has been also analysed.Our results (not illustrated here for space limitations) show that the performance gap between MPF and all the considered smoothing gets larger as σ 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, Fig. 6: CTB versus N p for SSM #2; MPF, SPS, RBSS, SRBSS and Alg-L are considered.and this impairs more and more the estimation process in the backward pass.

VI. CONCLUSIONS
In this manuscript the smoothing problem for SSMs has been analysed from a FG perspective.This has allowed us to devise two new particle smoothing algorithms for CLG SSMs.The first algorithm, dubbed SPS, generates estimates of the marginal smoothing densities at a limited computational cost.Moreover, as evidenced by our computer simulations for specific SSMs, its accuracy is similar to that achieved by the second algorithm, dubbed RBSS.The RBSS technique, however, generates an estimate of the joint smoothing pdf for the entire observation interval; in addition, its accuracy in state estimation and computational costs are comparable with those of the RBPS algorithm recently proposed by Lindsten et al. in [20].
Our future work concerns the application of FG methods to the problems of filtering and smoothing for other classes of SSMs.
assumption that x (N ) l is known), the other one to smoothing for x (N ) l (under the dual assumption that x (L) l is known).The first (second) sub-graph can be easily obtained from the FG shown in Fig. 1 by including the contribution of a pseudomeasurement, denoted z l ).In fact, smoothing for the linear component x (L) l the two sub-graphs are merged by adding five distinct equality nodes, associated with the shared variables (namely, y l , x

Fig. 2 :
Fig. 2: Factor graph for case C.2.The sub-graph referring to the linear (nonlinear) state component is identified by red (blue) lines, whereas the equality nodes introduced to merge the two sub-graphs by black lines.The direction of the messages passed over the half edges x (L) l and x (N ) l (inputs) and over the half edges x (L) l+1 and x (N ) l+1 (outputs) is indicated by green arrows.
k }) is Gaussian and all its elements have zero mean and covariance C (L) w (C (N ) w ) for any l; b) the measurement noise {e (L) -(a), concerns the linear state component, since it generates the particle-dependent statistical model ← m be,j (x (L) l ) (conveying a particle-dependent backward estimate of x (L) l ) and combines it with the forward model m f p,j (x (L) l ) on the basis of (13); this results in the message m sm,j (x (L) l ), conveying a particledependent statistical model for the smoothed pdf of x (L) l .

l
l−1),j , the random variable x (L) l+1 is generated by adding the output A l−1),j ).Consequently, the precision matrix W (L) 1,l,j and the transformed mean vector w (L) 1,l,j associated with the covariance matrix C (L) 1,l,j and the mean vector η (L) 1,l,j , respectively, can be easily computed by applying eqs.(IV.6)-(IV.8) of [23, Table 4, p.1304] in their backward form (with A → I D L , X → A ) and, then, eqs.(III.5)-(III.6) of [23, Table . The second step can be seen as a MU for the linear state component.In fact, it aims at generating the backward estimate ← m be,j (x that convey the available statistical information about the measurement y l and the pseudo-measurement z ,j = (B l,j ) T W e B l,j

Fig. 3 :
Fig. 3: Representation of the message scheduling employed in the first part (a) and in the second part (b) of the (T − l)-th recursion of SPS backward processing.Blue, green and red arrows are employed to identify the input forward messages, the input/output backward messages and the remaining messages, respectively.

Fig. 5 :
Fig. 5: RMSE performance versus N p for the linear state component (RM SE L ) and the nonlinear state component (RM SE N ) for SSM #2; MPF, SPS, RBSS and Alg-L are considered.