arXiv:1605.03017v3 [math.ST] 27 Jul 2016

Marginalized Particle Filtering and Related
Filtering Techniques as Message Passing
July 29, 2016
Abstract
In this manuscript a factor graph approach is employed to investigate the recursive filtering problem for mixed linear/nonlinear state-space
models. Our approach allows us to show that: a) the factor graph characterizing the considered filtering problem is not cycle free; b) in the case
of conditionally linear Gaussian systems, applying the sum-product rule,
together with different scheduling procedures for message passing, to this
graph results in both known and novel filtering techniques. In particular,
it is proved that, on the one hand, adopting a specific message scheduling
for forward only message passing leads to marginalized particle filtering
in a natural fashion; on the other hand, if iterative strategies for message passing are employed, a novel filtering method, dubbed turbo filter
for its conceptual resemblance to the turbo decoding methods devised for
concatenated channel codes, can be developed.

Giorgio M. Vitetta, Emilio Sirignano, Francesco Montorsi and Matteo Sola
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, sola.matteo87@gmail.com
Keywords: State Space Representation, Hidden Markov Model, Particle
Filter, Belief Propagation, Turbo Processing.

1

G. M. Vitetta et al., Marginalized Particle Filtering . . .

1

2

Introduction

The nonlinear filtering problem consists of inferring the posterior distribution of
the hidden state of a nonlinear dynamic system from a set of past and present
measurements [1]. It is well known that, if a nonlinear dynamic system can be
described by a state-space model (SSM), a general sequential procedure, based
on the Bayes’ rule and known as Bayesian filtering, can be easily derived for
recursively computing the posterior distribution of the system current state
[1]. Unluckily, Bayesian filtering is analytically tractable in few cases for the
following two reasons [2]: a) one of the two steps it consists of requires multidimensional integration which, in most cases, does not admit a closed form
solution; b) the functional form of the required posterior distribution may not
be preserved over successive recursions. For this reason, sequential techniques
employed in practice are based on various analytical approximations and, consequently, generate a functional approximation of the desired distribution. Such
techniques are commonly divided into local and global methods on the basis of
the way posterior distributions are approximated [3, 4, 5]. Local methods, like
extended Kalman filtering [6] and unscented filtering [7], are computationally
efficient, but may suffer from the problem of error accumulation over time. On
the contrary, global methods, like sequential Monte Carlo methods [8, 9] (also
known as particle filtering, PF, methods [10, 11, 12]) and point mass filtering
[5, 13] may achieve high accuracy at the price, however, of an unmanageable
complexity and numerical problems in the presence of a large dimension of system state [14]. These considerations have motivated various research activities
focused on the development of novel Bayesian filters able to achieve high accuracy under given computational constraints. Significant results in this research
area concern the use of the new representations for complex distributions, like
belief condensation filtering [3], and the development of novel filtering techniques combining local and global methods, like marginalized particle filtering
(MPF) [15, 16], and other methods originating from it [4, 17, 18]. Note that
the last class of methods applies to mixed nonlinear/nonlinear models [19], that
is to models whose state can be partitioned in a conditionally linear portion
(usually called linear state variable) and in a nonlinear portion (representing
the remaining part of system state and called nonlinear state variable). This
partitioning of system state allows to combine a global method (e.g., particle
filtering) operating on the nonlinear state variable with a local technique (e.g.,
Kalman filtering) involving the linear state variable only.
In this manuscript the factor graph (FG) approach illustrated by Loeliger et

G. M. Vitetta et al., Marginalized Particle Filtering . . .

3

al. in [20] is employed to revisit the problem of recursive Bayesian filtering for
mixed linear/nonlinear models from a perspective substantially different from
that adopted in MPF [15]. This allows us to shed new light on the problem of
filtering for mixed linear/nonlinear models, providing a new interpretation of
MPF and paving the way for the development of new filtering techniques. In
particular, based on this approach, we are able to show that: a) the considered
filtering problem can be formulated as a message massing problem over a specific
FG, which, unluckily, is not cycle free; b) in the case of a conditionally linear
Gaussian (CLG) SSM [19], MPF results from the application of the sum-product
algorithm (SPA) [20, 21], together with a specific scheduling procedure for forward only message passing, to this graph; c) our graphical representation leads,
in a natural fashion, to the development of novel filtering methods simplifying
and/or generalising it. As far as the last point is concerned, in our work specific
attention is paid to the development of a novel iterative filtering technique that
exploits the exchange of probabilistic (i.e., soft) messages to progressively refine
the posteriors of the linear and nonlinear state variables within each recursion
and is dubbed turbo filtering (TF) for its conceptual resemblance to the iterative
(i.e., turbo) decoding of concatenated channel codes.
It is important to point that our approach has been inspired by various
ideas and results already available in the technical literature concerning different
research areas; here, we limit to mention the following relevant facts:
• A mixed linear/nonlinear Markov system can be represented as the concatenation of two interacting subsystems, one governed by linear dynamics, the other one accounting for a nonlinear behavior; conceptually related
(finite state) Markov models can be found in data communications and,
in particular, in concatenated channel coding (e.g., turbo coding [22]) and
in coded transmissions over inter-symbol interference channels for which
turbo decoding methods [22, 23] and turbo equalization techniques [24] have
been developed, respectively1 .
• Factor graphs play an essential role in the derivation and interpretation
of turbo decoding and equalization [20, 26] (for instance, turbo decoding
techniques emerge in a natural fashion from graphical models of codes
[27]).
• Both Kalman filtering and particle methods can be viewed as message
1 Note that these classes of algorithms can be seen as specific applications of the so called
turbo principle [25], [35, Par. 10.5.1]

4

G. M. Vitetta et al., Marginalized Particle Filtering . . .

passing procedures on factor graphs, as shown in [20, 21] and in [28], respectively.
• Various methods to progressively refine distributional approximations through
multiple iterations have been developed in the field of Bayesian inference
on dynamic systems (even if implementations substantially different from
that we devise have been proposed), and, in particular, in expectation
propagation in Bayesian networks [29, 30] and in variational Bayesian
filtering [4]. Consequently, various links to previous work on Bayesian
inference on graphical models and variational Bayes methods [31] can be
also established.
The remaining part of this manuscript is organized as follows. The mathematical model of the considered class of mixed linear/nonlinear systems is illustrated in Section 2, whereas a representation of the filtering problem for these
systems through a proper FG is provided in Section 3. Then, it is shown that
applying the SPA and proper message scheduling strategies to this FG leads
to MPF in Section 4. This approach paves the way, in a natural fashion, for
the development of simplifications and generalizations of MPF, that are devised
in Sections 5 and 6, respectively. The novel filtering methods proposed in this
manuscript are compared, in terms of accuracy and computational effort, with
MPF in Section 7. Finally, some conclusions are offered in Section 8.
Notations: The probability density function (pdf) of a random vector R
evaluated at point r is denoted f (r); N (r; ηr , Cr ) represents the pdf of a Gaussian random vector R characterized by the mean ηr and covariance matrix Cr
evaluated at point r; the precision (or weight ) matrix associated with the covariance matrix Cr is denoted Wr , whereas the transformed mean vector Wr ηr
is denoted wr ; xi denotes the i-th element of the vector x.

2

System Model

In the following we focus on a discrete-time mixed linear/nonlinear SSM [15],
whose hidden state in the l-th interval is represented by a D-dimensional real
vector xl , [x0,l , x1,l , ..., xD−1,l ]T . We assume that this vector can be partitioned as

T 
 T T
(L)
(N )
xl = xl
, xl
,
(1)
(L)

(L)

(L)

(L)

(N )

(N )

(N )

(L)

where xl , [x0,l , x1,l , ..., xDL −1,l ]T (xl , [x0,l , x1,l , ..., xDN −1,l ]T ) is the
so called linear (nonlinear ) component of xl (1), with DL < D (DN = D − DL ).

5

G. M. Vitetta et al., Marginalized Particle Filtering . . .
(L)

This partitioning of xl is accomplished as follows. First, xl is identified as
that portion of xl (1) characterized by the following two properties:
1. Conditionally linear dynamics - This means that its update equation, con(N )
ditioned on xl , is linear, so that




(L)
(L)
(N )
(L)
(L)
(N )
(L)
xl+1 = Al
xl
xl + fl
xl
+ wl ,
(2)
(L)

(L)

(N )

where fl (x) is a time-varying DL -dimensional real function, Al (xl )
(L)
is a time-varying DL × DL real matrix and wl is the l-th element of
(L)
the process noise sequence {wk }, which consists of DL - dimensional
independent and identically distributed (iid) noise vectors.
2. Conditionally linear (or almost linear ) dependence of all the available
measurements on it - In other words, these quantities, conditioned on
(N )
(L)
xl , exhibit a linear dependence on xl (additional details about this
feature are provided below).
(N )

Then, xl
is generated by putting together all the components of xl that
(L)
do not belong to xl . For this reason, generally speaking, this vector is characterized by at least one of the following two properties:
a) Nonlinear dynamics - The update equation




(N )
(N )
(N )
(N )
(N )
(L)
(N )
xl+1 = fl
xl
+ Al
xl
xl + wl
(3)

is assumed in the following for the nonlinear component of system state, where
(N ) (N )
(N )
Al (xl ) is a time-varying DN × DL real matrix, fl (x) is a time-varying
(N )
DN -dimensional real function and wl is the l-th element of the process noise
(N )
sequence {wk }, which consists of DN -dimensional iid noise vectors and is
(L)
statistically independent of {wk }.
b) A nonlinear dependence of all the available measurements on it (further
details are provided below).
In the following Section we focus on the so-called filtering problem, which
concerns the evaluation of the posterior pdf f (xt |y1:t ) at an instant t > 1, given
a) the initial pdf f (x1 ) and b) the t · P -dimensional measurement vector
T

y1:t = y1T , y2T , ..., ytT ,

(4)

where yl , [y0,l , y1,l , ..., yP −1,l ]T denotes the P -dimensional real vector collecting all the noisy measurements available at time l. As already mentioned

G. M. Vitetta et al., Marginalized Particle Filtering . . .

6

above, the measurement vector yl exhibits a linear (nonlinear ) dependence on
(L)
(N )
xl (xl ), so that the model [18]




(N )
(L)
(N )
xl + el
(5)
+ Bl xl
yl = hl xl
(N )

(N )

can be adopted, where Bl (xl ) is a time-varying P × DL real matrix, hl (xl )
is a time-varying P -dimensional real function and el the l-th element of the
measurement noise sequence {ek } consisting of P -dimensional iid noise vectors
(N )
(L)
and independent of both {wk } and {wk }.

3

Representation of the Filtering Problem via
Factor Graphs

Generally speaking, the filtering problem for a SSM described by the Markov
model f (xl+1 |xl ) and the observation model f (yl |xl ) for any l concerns the
computation of the posterior pdf f (xt |y1:t ) for t ≥ 1 by means of a recursive
procedure [1]. It is well known that, if the pdf f (x1 ) is known, a general Bayesian
recursive procedure, consisting of a measurement update step followed by a time
update step, can be employed. In practice, in the first step of the l-th recursion
(with l = 1, 2, ..., t) the conditional pdf

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

1

f yl y1:(l−1)

(6)

is computed on the basis of pdf f (xl |y1:(l−1) ) (evaluated in the last step of the
previous recursion2 ), the present measurement vector yl and the pdf
Z


f yl y1:(l−1) = f (yl |xl ) f xl y1:(l−1) dxl .
(7)
In the second step f (xl |y1:l ) (6) is exploited to compute the pdf
Z
f (xl+1 |y1:l ) = f (xl+1 |xl ) f (xl |y1:l ) dxl ,

(8)

which represents a prediction about the future state xl+1 . It is important to
point out that: 1) the term 1/f (yl |y1:(l−1) ) appearing in the right hand side
(RHS) of (6) represents a normalization factor ; 2) both (7) and (8) require
integration with respect to xl and this may represent a formidable task when
the dimensionality of xl is large and/or the pdfs appearing in the integrands are
2 Note that in the first recursion (i.e., for l = 1) f (x |y
l 1:(l−1) ) = f (x1 |y1:0 ) = f (x1 ) and
f (yl |y1:(l−1) ) = f (y1 |y1:0 ) = f (y1 ), so that f (x1 |y1 ) = f (x1 )f (y1 |x1 )/f (y1 ).

G. M. Vitetta et al., Marginalized Particle Filtering . . .

7

not Gaussian; 3) this recursive procedure lends itself to be efficiently represented
by a message passing algorithm over a proper FG3 [28], in which each factor
of a product of functions is represented by a distinct node (a rectangle in our
diagrams), whereas each variable is associated with a specific (and, usually,
unoriented) edge or half edge. As far as the last point is concerned, we also
note that the derivation of this FG relies on the fact that the a posteriori pdf
f (xt |y1:t ) has the same FG as the joint pdf f (xt , y1:t ) (see [20, Sec. II, p. 1297])
and the last pdf can be computed recursively through a procedure similar to
that illustrated above, but in which the measurement update (6) and the time
update (8) are replaced by

and


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

Z

f (xl+1 |xl ) f (xl , y1:l ) dxl ,

(9)

(10)

respectively, so that the evaluation of the above mentioned normalization factor
is no more required. In fact, eqs. (9) and (10) involve only products of pdfs
and a sum (i.e., integration) of products, so that they can be represented by
means of the FG shown in Fig. 1 (where, following [20], a simplified notation
is used for the involved pdfs and the equality constraint node, which represents
an equality constraint “function”, that is a Dirac delta function). Since this
FG is cycle free, the pdf f (xl , y1:l ) can be evaluated applying the well known
sum-product rule 4 (i.e., the SPA) to it, i.e. developing a proper mechanism for
passing probabilistic messages along this FG (the flow of messages is indicated
by red arrows in Fig. 1). In fact, if the input message m
~ in (xl ) = f (xl , y1:(l−1) )
enters this FG, the message going out of the equality node is given by
m
~ e (xl ) = m
~ in (xl ) f (yl |xl ) ,

(11)

so that m
~ e (xl ) = f (xl , y1:l ) (see (9)); then, the message emerging from the
function node referring to the pdf f (xl+1 |xl ) is expressed by
Z
m
~ out (xl+1 ) = f (xl+1 |xl ) m
~ e (xl ) dxl ,
(12)
so that m
~ out (xl+1 ) = f (xl+1 , y1:l ) = m
~ in (xl+1 ) (see (10)). From this result it
can be easily inferred that the pdf f (xt , y1:t ) (and, up to a scale factor, the pdf
3 Forney-style

factor graphs are always considered in the following [20].
a Forney-style FG, such a rule can be formulated as follows [20]: 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.
4 In

G. M. Vitetta et al., Marginalized Particle Filtering . . .

8

xl
min ( xl )

=

f y l / xl

f xl +1 / xl

yl

me ( xl )
x l +1

mout ( xl +1 ) = min ( xl +1 )

Figure 1: Factor graph representing (9) and (10). The SPA message flow characterizing the l-th recursion of Bayesian filtering is indicated by red arrows.
f (xt |y1:t )) results from the application of the SPA to the overall FG originating
from the ordered concatenation of multiple subgraphs, each structured like the
one shown in Fig. 1 and associated with l = 1, 2, ..., t. In this graph the flow of
messages produced by the SPA proceeds from left to right, i.e. the pdf f (xt , y1:t )
is generated by a forward only message passing. Note also that, in principle, the
desired pdf f (xt , y1:t ) is computed as the product between two messages, one
for each direction, reaching the rightmost half edge of the overall FG, but one of
− (x ) = 1.
the two incoming messages for that edge is the constant function ←
m
he t
Unluckily, as the size D of xl (1) gets large, the computational burden associated with (6)-(8) (or, equivalently, (9) and (10)) becomes unmanageable. In
principle, a substantial complexity reduction can be achieved decoupling 5 the
(L)
(N )
(L)
filtering problem for xl from that for xl , i.e. the evaluation of f (xl |y1:l )
(N )
from that of f (xl |y1:l ). In fact, this approach potentially provides the following two benefits: a) a given filtering problem is turned into a couple of filtering
problems of smaller dimensionality and b) some form of computationally efficient
standard filtering (e.g., Kalman or extended Kalman filtering) can be hopefully
(L)
exploited for the linear portion xl of the state vector xl (1). As a matter of
fact, these principles have been exploited in devising MPF [10, 15] and, as it will
become clearer in the following, they must be always kept into account in the
derivation of our FG representation. Before illustrating this derivation, however, the measurement and state models on which such a representation relies
need to be clearly defined; for this reason, these models are analysed in detail
in the following part of this Section. To begin, let us concentrate on the models
(L)
(L)
involved in the filtering problem for xl , i.e. on the evaluation of f (xl |y1:l ),
5 Note

that, generally speaking, in the considered problem the coupling of the filtering
(L)
(N)
problem for xl with that for xl
is due not only to the structure of the update equations
(2) and (3), but also to the measurement vector yl (5), since this exhibits a mixed dependence
on the two components of the state vector xl (1).

9

G. M. Vitetta et al., Marginalized Particle Filtering . . .
(N )

under the assumption that the nonlinear portion xl
of the system state is
known for any l. In this case, the evaluation of this pdf can benefit not only
from the knowledge of yl (5), but also from that of the quantity (see (3))




(L)
(N )
(N )
(N )
(N )
(N )
(L)
(N )
zl , xl+1 − fl
xl
= Al
xl
xl + wl ,
(13)

which can be interpreted as a pseudo-measurement [15], since it does not originate from real measurements, but from the constraints expressed by the state
equation (3). This leads to considering the overall observation model


(L)
(L)
(N )
f yl , zl xl , xl

 

(14)
(L)
(N )
(L)
(L)
(N )
= f yl xl , xl
f zl xl , xl
(L)

for xl , where


(L)
(N )
f yl xl , xl
= f (el )|e =y −B x(N ) x(L) −h x(N )  ,
l

l

l

l

l

l

(15)

l

and




(L)
(L)
(N )
(N )
f zl xl , xl
= f wl

(N )

wl

(L)

=zl

(N )

−Al



(L)
(N )
xl
xl

.

(16)

If the observation model (14) and the state model (see (2))


(L)
(L)
(N )
f xl+1 xl , xl






(L)
(L)
(N )
(L)
(N )
(L)
xl
− Al
xl
xl
= fw(L) xl+1 − fl

(17)

(L)

are adopted for xl , the graph identified by the blue lines and rectangles appearing in Fig. 1 can be drawn. Then, in principle, if the sum-product rule is
(N )
(N )
applied to it under the assumption that the couple (xl , xl+1 ) is known for any
l, the expressions of the messages flowing in the overall graph for the evaluation
(L)
(L)
of f (xt , y1:t , z1:t ) can be easily derived. It is important to point out that:
• This graph contains a node which does not refer to the above mentioned
density factorizations6 , but represents the transformation from the couple
(N )
(N )
(L)
(xl , xl+1 ) to zl (see (13)); this feature of the graph has to be carefully
kept into account when deriving message passing algorithms.
(L)

(L)

(N )

• Generally speaking, the evaluation of the conditional pdf f (zl |xl , xl )
(N )
(N )
requires the knowledge of the joint pdf of xl+1 and xl
conditioned on
(L)
xl (see (13)).
6 This peculiarity is also evidenced by the presence of an arrow on all the edges connected
to such a node.

10

G. M. Vitetta et al., Marginalized Particle Filtering . . .

The same line of reasoning can be followed for the filtering problem concern(N )
(L)
ing xl . Consequently, in this case the linear portion xl of the system state
is assumed to be known for any l and the pseudo-measurement (see (2))




(N )
(L)
(L)
(N )
(L)
(L)
(N )
(L)
zl , xl+1 − Al
xl
xl = fl
xl
+ wl
(18)
is defined. This leads to the overall observation model


(N )
(N )
(L)
f yl , zl
xl , xl
=

 

(N )
(L)
(N )
(N )
f yl xl , xl
f zl
xl
(N )

(N )

(N )

(19)

(L)

(L)

(N )

for xl , where f (zl |xl ) can be expressed similarly as f (zl |xl , xl )
(see (16)). Then, if the observation model (19) and the state model


(N )
(N )
(L)
f xl+1 xl , xl






(20)
(N )
(N )
(N )
(N )
(N )
(L)
xl
− Al
xl
xl
= fwN xl+1 − fl
(N )

are adopted for xl , the red graph of Fig. 2 can be drawn and exploited in
(N )
(N )
a similar way as the blue graph for the evaluation of f (xt , y1:t , z1:t ), under
(L)
(L)
the assumption that the couple (xl , xl+1 ) is known for any l. Note also
that, similarly to what has been mentioned earlier about the conditional pdf
(L) (L)
(N )
(N ) (N )
f (zl |xl , xl ), the evaluation of f (zl |xl ) requires the knowledge of the
(L)
(L)
(N )
joint pdf of xl+1 and xl conditioned on xl .
Finally, merging the blue graph with the red one (i.e., adding four equality
(L)
(N )
(L)
(N )
constraint nodes for the variables xl , xl , xl+1 and xl+1 shared by the red
graph and the blue one) produces the overall FG illustrated in Fig. 2. Given this
FG, we would like to follow the same line of reasoning as that adopted for the FG
(L)
(L)
of Fig. 1. In other words, given the input messages m
~ in (xl ) = f (xl , y1:(l−1) )
(N )
(N )
and m
~ in (xl ) = f (xl , y1:(l−1) ) (entering the FG along the half edges asso(L)
(N )
ciated with xl and xl , respectively), we would like to derive a forward
only message passing algorithm based on this FG and generating the output
(L)
(L)
(N )
(N )
messages m
~ out (xl+1 ) = f (xl+1 , y1:l ) and m
~ out (xl+1 ) = f (xl+1 , y1:l ) (emerging
(L)
(N )
from the FG along the half edges associated with xl+1 and xl+1 , respectively)
on the basis of the available a priori information and the noisy measurement
yl . Unluckily, the new FG, unlike that shown in Fig. 1, is not cycle-free, so
that any application of the SPA to it unavoidably leads to approximate results
[21], whatever message scheduling procedure [21, 27] is adopted. This consideration must be carefully kept into account in both the derivation of MPF as
a message passing algorithm and in the development of possible modifications
and generalizations of this technique, as it will become clearer in Sections 4-6.

11

G. M. Vitetta et al., Marginalized Particle Filtering . . .


mout ( xl(+L1) )

xl( L )


min ( xl( L ) )

( L)
xl+
1

f x( L ) / x( L ) , x ( N )
l +1

=

=

l

l

=

=
f z ( L ) / x ( L ) ,x ( N )
l

z

l

l

( L)
l

xl( N ) , xl(+N1) → z l( L )

=

f y / x ( L ) ,x ( N )
l

l

l

yl

=

xl( L ) , xl(+L1) → z l( N )
z l( N )
f z ( N ) / x( N )
l

x

(N )
l

=

min ( xl( N ) )

=

l

=

=
f x( N ) / x( N ) , x( L )
l +1
l
l

mout ( xl(+N1) )

(N )
xl+
1

Figure 2: Overall factor graph resulting from the merge of two subgraphs, one
(L)
(N )
referring to filtering for xl (in blue), the other one to that for xl (in red).
The equality constraint nodes introduced to connect these subgraphs are identi(L)
(N )
fied by black lines. The flow of the messages along the half edges xl and xl
(L)
(N )
(input) and that of the messages along the half edges xl+1 and xl+1 (output)
are indicated by green arrows.

G. M. Vitetta et al., Marginalized Particle Filtering . . .

4

12

Message Passing in Marginalized Particle Filtering

In the following Section we show how the equations describing the l-th recursion of MPF result from the application of the SPA to the FG shown in Fig. 2.
However, before illustrating the detailed derivation of such equations, it is important to discuss the following relevant issues. First of all, the MPF technique
has been developed for the specific class of GLG SSMs [15, 19], to which we always refer in the following discussion. In particular, in the following we assume
(L)
(N )
that: a) {wk } ({wk }) is a Gaussian random process and all its elements
(L)
(N )
(L)
have zero mean and covariance Cw (Cw ) for any l; b) {ek } is a Gaussian
random process having zero mean and covariance matrix Ce for any l; c) all the
above mentioned Gaussian processes are statistically independent. Under these
(L)
(N
(L) (L)
(L) (L)
(N )
assumptions, the pdfs f (yl |xl , xl ), f (zl |xl ) and f (xl+1 |xl , xl ) (see
(N ) (L)
(N )
(15)-(17)) are Gaussian with mean (covariance matrix) Bl (xl )xl +hl (xl ),
(N ) (N )
(L)
(L) (N )
(L) (N )
(L)
(N )
Al (xl ) xl and fl (xl ) + Al (xl ) xl , respectively (Ce , Cw and
(N ) (N )
(N ) (N )
(L)
(L)
Cw , respectively). Similarly, the pdfs f (zl |xl ) and f (xl+1 |xl , xl )
(L) (N )
(see (19) and (20)) are Gaussian with mean (covariance matrix) fl (xl ) and
(N )
(N ) (N )
(N ) (N )
(L)
(L)
fl (xl ) + Al (xl ) xl , respectively (Cw and Cw , respectively).
Secondly, as explained below in detail, the MPF can be interpreted as a
forward only message passing algorithm operating over the FG shown in Fig.
2. The scheduling procedure adopted for MPF unavoidably leads to ignoring
(N )
the evaluation of the pseudo-measurement zl (18). For this reason, in the following we refer to the simplified FG shown in Fig. 3, which has been obtained
from that illustrated in Fig. 2 removing the block representing the transfor(L)
(L)
(N )
mation from (xl , xl+1 ) to zl
and the edges referring to the evaluation of
the last vector. Note that: a) in the new graph the block referring to the pdf
(L)
(N
f (yl |xl , xl ) appears twice, since this pdf is involved in the two subgraphs
shown in Fig. 2; b) some brown edges and equality nodes have been added to
(N )
(L)
feed such blocks with min (xl ) and min (xl ), since these represents the only a
(N )
(L)
priori information available about xl and xl , respectively, at the beginning
of the l-th recursion.
Thirdly, in MPF a particle-based model and a Gaussian model are adopted
(N )
(L)
for the input and the output messages referring to xl and xl , respectively,
and the functional structure of the generated messages is preserved in each
recursion. More specifically, on the one hand, the a priori information avail(N )
able about xl
at the beginning of the l-th recursion is represented by a set
(N )
(N )
of Np particles Sl/(l−1) = {xl/(l−1),j , j = 0, 1, ..., Np − 1} and their weights

G. M. Vitetta et al., Marginalized Particle Filtering . . .

13

{wl/(l−1),j , j = 0, 1, ..., Np − 1)}; following [15], we assume that such weights
are uniform (in other words, wl/(l−1),j = 1/Np for j = 0, 1, ..., Np − 1), so
that they can be ignored in the following derivation. On the other hand,
(L)
the a priori information available about xl is represented by a set of Gaussian pdfs, each associated with a specific particle; in particular, the Gaussian
(L)
(L) (L)
model N (xl ; ηl/(l.−1),j , Cl/(l−1),j ) is associated with the j-th particle (with
j = 0, 1, ..., Np − 1) at the beginning of the same recursion. From the last point
it can be inferred that, in developing a message passing algorithm that represents the MPF technique, we can focus on: a) a single particle contained in the
(N )
(N )
input message min (xl ) and, in particular, on the j-th particle xl/(l−1),j ; b)
(L)

(L)

(L)

on the Gaussian model N (xl ; ηl/(l.−1),j , Cl/(l−1),j ) associated with that particle. For this reason, we assume that, at the beginning of the l-th recursion,
(L)
(N )
our knowledge about xl and xl is condensed in the message




(L)
(L) (L)
(L)
(21)
= N xl ; ηl/(l.−1),j , Cl/(l−1),j
m
~ in,j xl

and in the message





(N )
(N )
(N )
= δ xl − xl/(l−1),j ,
m
~ in,j xl

(22)

respectively, with j = 0, 1, ..., Np − 1; these are processed to generate the corre(L)
(N )
sponding output messages m
~ out,j (xl+1 ) and m
~ out,j (xl+1 ), which are required to
(L)
(N )
have the same functional form as m
~ in,j (xl ) (21) and m
~ in,j (xl ) (22), respec(N )
tively. For this reason, the algorithm for computing m
~ out,j (xl+1 ) is expected
(N )
to generate a new particle x(l+1)/l,j with a (uniform) weight w(l+1)/l,j = 1/Np ;
(L)

similarly, that for evaluating m
~ out,j (xl+1 ) is expected to produce a new Gaus(N )
(L)
(L)
(L)
sian pdf N (xl+1 ; η(l+1)/l,j , C(l+1)/l,j ) associated with the particle x(l+1)/l,j (note
that, in deriving this pdf, possible scale factors are unrelevant and, consequently,
can be dropped).
(L)
(N )
Given m
~ in,j (xl ) (21) and m
~ in,j (xl ) (22), if the SPA is applied to the
considered graph and the message scheduling illustrated in Fig. 3 (and, as a
matter of fact, adopted in MPF) is employed, the steps described below are
(L)
(N )
carried out to evaluate m
~ out,j (xl+1 ) and m
~ out,j (xl+1 ) in the l-th recursion of
MPF7 .
(N )
1. Measurement update for xl
- This step aims at updating the weight
(N )
of the j-th particle xl/(l−1),j on the basis of the new measurements yl (this
7 In the following derivations some mathematical results about Gassian random variables
(e.g., see [32, Par. 2.3.3]) and Gaussian message passing in linear models (e.g., see [20, Table
2, p. 1303]) are exploited. As far as the MPF formulation is concerned, we always refer to
that given by algorithm 1 in [15, Sec. II].

14

G. M. Vitetta et al., Marginalized Particle Filtering . . .


min , j ( x (l L ) )

x (l L )

yl


min , j ( x l( L ) )


m2, j ( x (l L ) )


m2, j ( x l( L ) )


m4, j ( x l( L ) )

=

f x ( L ) / x ( L ) ,x ( N )
l +1

=
fy


m1, j ( x l( L ) )
l

,/ x l( L ) ,x l( N )

=

m3, j ( x (l L ) )

l


m j ( z l( L ) )

z

l

x (l +L1)


m5, j ( x l(+L1) )

= mout , j ( x (l +L1) )

f z( L ) / x ( L ) ,x ( N )
l


min , j ( x l( N ) )

l

=
l


m2, j ( x l( N ) )

( L)
l

x (l N ) , x l(+N1) → z l( L )

m5, j ( x (l +N1) )

m2, j ( x l( N ) )

=


m2, j ( x (l L ) )


m2, j ( x (l L ) )


min , j ( x l( L ) )

fy

l

=
x (l N )


m2, j ( x l( N ) )


m1, j ( x l( N ) )

,/ x l( L ) ,x l( N )

=

min , j ( x l( N ) )


m2, j ( x (l N ) )

yl


m2, j ( x (l N ) )


m5, j ( x (l +N1) )

=

=
f x ( N ) / x ( N ) ,x ( L )
l +1

l

l


mout , j ( x l(+N1) )

x (l +N1)

Figure 3: Overall factor graph for the representation of MPF processing; this
graph is obtained from the one shown in Fig. 3 removing the part referring to
(N )
the evaluation of zl and inserting two new (brown) equality constraints and
(L)
(N )
some (brown) edges referring to xl and xl . The message flow characterizing
MPF and referring to the j-th particle is also shown.

15

G. M. Vitetta et al., Marginalized Particle Filtering . . .

corresponds to step 2) of algorithm 1 in [15, Sec. II]). It involves the computation
of the messages


 Z 


(N )
(L)
(N )
(L)
(L)
m
~ 1,j xl
m
~ in,j xl
= f yl xl , xl
dxl
(23)
and






(N )
(N )
(N )
,
m
~ 1,j xl
=m
~ in,j xl
m
~ 2,j xl

(24)







(N )
(N )
(L)
(N )
(N )
ηl/(l.−1),j + hl xl
η1,l,j xl
, Bl xl

(26)

which provides the new importance weight for the considered particle (see Fig.
(N )
(L)
3). Substituting the expression of f (yl xl , xl ) (see (15)) and (21) in (23)
produces, after some manipulation (see the Appendix)







(N )
(N )
(N )
(N )
(N )
= N yl ; η1,l,j xl
, C1,l,j xl
,
(25)
m
~ 1,j xl
where

and

T





(N )
(N )
(L)
(N )
(N )
+ Ce .
Cl/(l−1),j Bl xl
C1,l,j xl
, Bl xl
(N )

(27)

(N )

Then, substituting m
~ in,j (xl ) (22) and m
~ 1,j (xl ) (25) in (24) yields




(N )
(N )
(N )
= wl,j δ xl − xl/(l−1),j ,
m
~ 2,j xl

(28)

where8



(N )
(N )
wl,j , N yl ; η1,l,j , C1,l,j

(29)
(N )

is the new particle weight combining the a priori information about xl with
the information provided by the new measurements; here (see (26) and (27))


(L)
(N )
(N )
(N )
(30)
η1,l,j , η1,l,j xl/(l−1),j = Bl,j ηl/(l.−1),j + hl,j

and



(L)
(N )
(N )
(N )
C1,l,j , C1,l,j xl/(l−1),j = Bl,j Cl/(l−1),j BTl,j + Ce ,

(31)

(N )

(N )

with hl,j , hl (xl/(l−1),j ) and Bl,j , Bl (xl/(l−1),j ). In MPF, after normal(w)

ization9 of the particle weights {wl,j } (i.e., after dividing them by Pl
(N)

,

8 In evaluating this weight, the factor [det(C
−P /2 appearing in the expression of the
1,l,j )]
involved Gaussian pdf is usually neglected, since this entails a negligible loss in estimation
accuracy.
9 Note that normalization requires the knowledge of all the weights {w } (29) and that,
l,j
unlike it, all the previous and following tasks can be carried out in parallel (i.e., on a particleby-particle basis).

16

G. M. Vitetta et al., Marginalized Particle Filtering . . .
Np −1

X

wl,j ), particle resampling with replacement is accomplished (this corre-

j=0

sponds to step 3) of algorithm 1 in [15, Sec. II]). Note that, even if resampling
does not emerge from the application of SPA to the considered graph, its use,
as it will become clearer at the end of this Section, plays an important role in
(N )
the generation of the new particles for xl+1 . Moreover, it can be easily incorporated in our message passing; in fact, resampling simply entails that Np parti(N )
cles {xl/(l−1),j } and their associated weights {wl,j } (29) are replaced by the new
(N )

(N )

particles {xl/l,j } (forming the new set Sl/l ) and their weights {wl/l,j = 1/Np },
(N )

respectively. Consequently, m
~ 2 (xl ) (28) is replaced by




(N )
(N )
(N )
= δ xl − xl/l,j ,
m
~ 2,j xl

(32)

since the particle weight does not depend on the index j. It is also worth
(L)
mentioning that, after resampling, the set of Gaussian messages {m
~ in,j (xl )}
(21) needs to be properly reordered and that the messages associated with all
the discarded particles are not propagated to the next steps.
(L)
2. Measurement update for xl - This step aims at updating our statistical
(L)
knowlege about xl on the basis of the new measurement yl (and corresponds
to step 3-a) of algorithm 1 in [15, Sec. II]). It involves the computation of the
messages



 Z 

(N )
(N )
(L)
(N )
(L)
dxl
(33)
m
~ in,j xl
= f yl xl , xl
m
~ 1,j xl
and







(L)
(L)
(L)
,
m
~ 1,j xl
=m
~ in,j xl
m
~ 2,j xl

(34)
(L)

which represents the output of the measurement update for xl (see Fig. 3).
Substituting (15) and (22) in (33) produces (see [32, Par. 2.3.3, eq. (2.115)])




(L)
(L)
(35)
= N yl ; Bl,j xl + hl,j , Ce
m
~ 1,j xl
which, after some manipulation (in which unrelevant scale factors are dropped),
can be put in the Gaussian form




(L)
(L) (L)
(L)
= N xl ; η1,l,j , C1,l,j ,
(36)
m
~ 1,j xl

with

(L)

(L)

(L)

w1,l,j , W1,l,j η1,l,j = BTl,j We (yl − hl,j ) ,

−1
(L)
(L)
W1,l,j , C1,l,j
= BTl,j We Bl,j

(37)
(38)

17

G. M. Vitetta et al., Marginalized Particle Filtering . . .
and We , C−1
e . Then, substituting (21) and (36) in (34) yields




(L)
(L)
(L) (L)
= N xl ; ηl/(l−1),j , Cl/(l−1),j
m
~ 2,j xl


(L) (L)
(L)
·N xl ; η1,l,j , C1,l,j ,

which can be reformulated as




(L)
(L) (L)
(L)
= N xl ; η2,l,j , C2,l,j ,
m
~ 2,j xl

(39)

(40)

if scale factors are ignored; here,
(L)

(L)

(L)

(L)

(41)


−1
(L)
(L)
(L)
(L)
W2,l,j , C2,l,j
= Wl/(l−1),j + W1,l,j ,

(42)

(L)

(L)

(L)

(L)

(L)

(L)

w2,l,j , W2,l,j η2,l,j = wl/(l−1),j + w1,l,j ,

Wl/(l−1),j , (Cl/(l−1),j )−1 and wl/(l−1),j , Wl/(l−1),j ηl/(l−1),j .
(N )

3. Time update for xl - This step aims at generating the j-th particle for
and its associated weight (this corresponds to step 3-b) of algorithm 1 in
(N )
[15, Sec. II]); these information are conveyed by the message m
~ 5,j (xl+1 ), which
can be expressed as (see Fig. 3)

 RR 

(N )
(N )
(L)
(N )
f xl+1 xl , xl
m
~ 5,j xl+1 =




(43)
(N )
(L)
(L)
(N )
·m
~ 2,j xl
m
~ 2,j xl
dxl dxl .

(N )
xl+1

The double integral appearing in the RHS of the last equation can be evaluated
(N )
as follows. First of all, substituting m
~ 2,j (xl ) (32) in (43) yields




 Z 
(L)
(L)
(N )
(N )
(L)
(N )
dxl .
(44)
~ 2,j xl
m
~ 5,j xl+1 = f xl+1 xl , xl/l,j m
(N )

(N )

(L)

(L)

Then, substituting the expression of f (xl+1 |xl , xl ) (see (20)) and m
~ 2,j (xl )
(40) in (44) yields, after some manipulation, the Gaussian message (see [32, Par.
2.3.3, eq. (2.115)])




(N )
(N ) (N )
(N )
(45)
m
~ 5,j xl+1 = N xl+1 ; η5,l,j , C5,l,j ,
where

(N )

(N )

C5,l,j
(N )

(N )

Al,j , Al

(N )

(N ) (L)

(N )

η5,l,j , Al,j η2,l,j + fl,j ,

T
(N ) (L)
(N )
)
, C(N
+
A
C
A
,
w
l,j
2,l,j
l,j

(N )

(N )

(46)
(47)

(N )

(xl/l,j ) and fl,j , fl (xl/l,j ). Note that, in principle,




(N )
(N )
~ 5,j xl+1 ,
m
~ out,j xl+1 = m

(48)

18

G. M. Vitetta et al., Marginalized Particle Filtering . . .

as it can be easily inferred from Fig. 3. However, as already mentioned above,
(j)
(N )
in MPF the output message m
~ out (xl+1 ) is required to have the same functional
(j) (N )
form as m
~ in (xl ) (22). This result can be achieved a) sampling the Gaussian
(N ) (N )
(N )
(N )
function N (xl+1 ; η5,l,j , C5,l,j ) (see (45)), that is drawing a sample x(l+1)/l,j from
(N )

it and b) assigning to the sample x(l+1)/l,j a probability w(l+1)/l,j equal to the
weight wl/l,j = 1/Np (originating from resampling). It is worth pointing out
(N )
that: 1) the particles {x(l+1)/l,j } form the new set S(l+1)/l ; 2) in accomplishing
step a) of this procedure, it may be useful to introduce artificial noise (this can
be simply done adding the same positive quantity to the diagonal elements of
(N )
the matrix Cw appearing in the RHS of (47)) in order to mitigate the so called
(N )
degeneracy problem [1, 34]. If this approach is adopted, the message m
~ 5,j (xl+1 )
(45) is replaced by




(N )
(N )
(N )
(49)
m
~ 5,j xl+1 = δ xl+1 − x(l+1)/l,j ,
(N )

which emerges from the graph as m
~ out,j (xl+1 ). This message is also used in the
(L)
time update for xl , as illustrated in the next step.
(L)
4. Time update for xl - This step aims at generating a new Gaussian
(L)
(N )
~ 5,j (xl+1 ) =
pdf associated with the j-th particle x(l+1)/l,j and conveyed by m
(L)

m
~ out,j (xl+1 ) (this corresponds to step 3-c) of algorithm 1 in [15, Sec. II]).
However, before doing that, a further measurement update is accomplished on
(L)
the basis of the pseudo-measurement zl (13). This involves the evaluation of
(L)
the messages m
~ j (zl ),

 

 Z

(L)
(L)
(L)
(N )
(L)
(L)
f zl xl , xl
dzl
(50)
= m
~ j zl
m
~ 3,j xl
and






(L)
(L)
(L)
m
~ 4,j xl
=m
~ 2,j xl
m
~ 3,j xl
,

(51)

(L)

as shown in Fig. 3. Generally speaking, the message m
~ j (zl ) can be expressed
as

 RR 

(L)
(L)
(N )
(N )
=
f zl xl , xl+1
m
~ j zl

 

(52)
(N )
(N )
(N )
(N )
(N )
·f xl+1 xl
f xl
dxl dxl+1 .
(N )

(N )

However, since in this case f (xl

) = δ(xl

(N )
x(l+1)/l,j )

(N )
m
~ 2,j (xl )

(N )

(N )

(L)

(N +1)

−xl/l,j ), f (xl+1 |xl ) = δ(xl
(N )
m
~ 5,j (xl+1 )

−

can be assumed (see
(32) and
(49), respectively), eq. (52) easily leads to the expression






(L)
(L)
(N )
(L)
(L)
(N )
(53)
= f zl xl/l,j , x(l+1)/l,j = δ zl − zl,j ,
m
~ j zl

19

G. M. Vitetta et al., Marginalized Particle Filtering . . .

where
(L)

(N )

(N )

zl,j , x(l+1)/l,j − fl,j .

(54)
(L)

(L)

(N )

Then, substituting (53) and the expression of f (zl |xl , xl
(50) yields




(L)
(L)
(N ) (L)
)
.
= N zl,j ; Al,j xl , C(N
m
~ 3,j xl
w

) (see (16)) in

Finally, substituting the last expression and (40) in (51) produces




(L)
(L) (L)
(L)
= N xl ; η2,l,j , C2,l,j
m
~ 4,j xl


(L)
(N ) (L)
(N )
·N zl,j ; Al,j xl , Cw
,

(55)

(56)

which, after some manipulation (in which unrelevant scale factors are dropped),
can be rewritten as




(L)
(L) (L)
(L)
= N xl ; η4,l,j , C4,l,j ,
(57)
m
~ 4,j xl

where


T
(L)
(L) (L)
(L)
(N )
(N ) (L)
w4,l,j , W4,l,j η4,l,j = w2,l,j + Al,j
Ww
zl,j ,

(N )


T

−1
(L)
(N )
(L)
(L)
(N ) (N )
Ww
Al,j
= W2,l,j + Al,j
W4,l,j , C4,l,j

(58)
(59)

(N )

and Ww , [Cw ]−1 .
(L)
The last part of the time update step for xl requires the evaluation of the
output message

 RR 

(L)
(L)
(L)
(N )
f xl+1 xl , xl
m
~ 5,j xl+1 =




(60)
(N )
(L)
(N )
(L)
dxl dxl ,
m
~ 2,j xl
·m
~ 4,j xl
(N )

which, similarly as m
~ 5,j (xl+1 ) (43), requires double integration. Substituting
(N )
m
~ 2,j (xl ) (32) in the RHS of the last expression yields




 Z

(L)
(L)
(L)
(L)
(L)
(N )
dxl .
(61)
~ 4,j xl
f xl+1 xl , xl/l,j m
m
~ 5,j xl+1 =
(L)

(L)

(N )

Then, substituting the expression of f (xl+1 |xl , xl/(l−1),j ) (see (17)) and (57)
in the last equation gives (see [32, Par. 2.3.3, eq. (2.115)])




(L)
(L)
(L)
(L)
m
~ 5,j xl+1 = N xl+1 ; η5,l,j , C5,l,j


(62)
(L)
=m
~ out,j xl+1 ,

where

(L)

(L) (L)

(L)

(L)

η5,l,j , Al,j η4,l,j + fl,j = η(l+1)/l,j ,

(63)

20

G. M. Vitetta et al., Marginalized Particle Filtering . . .

(L)

(L)


T
(L)
(L) (L)
(L)
(L)
C5,l,j , C(L)
+
A
C
A
= C(l+1)/l,j ,
w
l,j
4,l,j
l,j

(N )

(L)

(L)

(64)
(L)

(N )

~ out,j (xl+1 ) (62)
fl,j , fl (xl/l,j ) and Al,j , Al (xl/l,j ). The evaluation of m
concludes the MPF message passing procedure, which needs to be carried out
for each of the Np particles available at the beginning of the l-th recursion.
Note that this procedure needs a proper inizialization (this corresponds to step
1) of algorithm 1 in [15, Sec. II]). In practice, before starting the first recursion
(N )
(N )
(corresponding to l = 1), the set S1/0 = {x1/0,j , j = 0, 1, ..., Np − 1}, consisting
(N )

of Np particles, is generated for x1

sampling the pdf
 Z

(L)
(N )
= f (x1 ) dx1 ,
f x1

(65)
(L)

(L)

(L)

and the same weight w1/0 = 1/Np and Gaussian model N (x1 ; η1/0 , C1/0 ) for
(L)

are assigned to each of them.
Finally, it is worth pointing out that: 1) the processing accomplished in the
(L)
measurement and time update for xl can be interpreted as a form of Kalman
filtering, in which both the real measurement yl and the pseudo-measurement
(L)
(N )
(L)
zl are processed [15]; 2) in the l-th recursion estimates of xl and xl can
P
Np −1
(w)
(L)
(N )
(N )
(see (28)) and as x̂l =
be evaluated as x̂l
=
j=0 wl,j xl/(l−1),j /Pl
PNp −1 (L)
j=0 η4,l,j /Np (see (58)), respectively; 3) the result expressed by eq. (45)
shows that, generally speaking, the statistical representation generated by the
(N )
SPA for the state xl+1 is a Gaussian mixture (GM), whose Np components have
the same weight (equal to 1/Np ) because resampling is always used in step 1.
The last point leads to the conclusion that, if resampling was not accomplished
in the l-th recursion, the weight of the j-th component of this GM would be
proportional by wl,j (29); this would unavoidably raise the problem of sampling
(N )
a GM with unequally weighted components in generating the particle set Sl/(l+1)

x1

(L)

and that of properly handling the resulting pseudo-measurements {zl,j }. These
considerations motivate the use of resampling in each recursion, indipendently
(N )
of the effective sample size [1] characterizing the particle set Sl/(l−1) .

5

Simplifying Marginalized Particle Filtering

The MPF derivation illustrated in the last two Sections unveils the real nature
of MPF and its limitations, and shows the inner structure of the processing
accomplished within each step. For these reasons, it paves the way for the development of new filtering methods related to MPF. In this Section we exploit
our FG-based representation of Bayesian filtering to develop reduced complexity

G. M. Vitetta et al., Marginalized Particle Filtering . . .

21

alternatives to MPF by simplifying the message passing derived in the previous
Section. It is worth mentioning that some methods for reducing MPF computational complexity [16] have been already proposed in the technical literature
[4], [17], [18]. In particular, the method proposed in [4] and [17] is based on rep(N )
resenting the particle set for xl as a single particle (that corresponds to the
center of mass of the set itself), so that a single Kalman filter is employed in up(L)
dating the statistics of the linear component xl ; consequently, the statistical
(L)
knowledge about xl is condensed in the single message




(L)
(L)
(L) (L)
(66)
= N xl ; ηl/(l.−1) , Cl/(l−1) ,
m
~ in xl
(L)

instead of the Np messages {m
~ in,j (xl )} (21). Unluckily, this simplified MPF
(N )
algorithm works well only if the posterior distribution of xl is unimodal. Its
(N )
generalization to the case in which the posterior distribution of xl
is multimodal has been illustrated later in [18]. In the proposed technique the particles
are partitioned into Kl groups or clusters (the parameter Kl is required to equal
(N )
the number of modes of the posterior density of xl ) and each group is represented by a single particle that corresponds to its center of mass; this allows
to reduce the overall number of Kalman filters from Np to Kl . The implementation of this technique requires, however, solving the following two specific
problems: a) identifying the number of modes of the posterior distribution of
(N )
xl ; b) partitioning the particles into clusters according to a grouping method
in each recursion. Unluckily, practical solutions for suche problems have not
been proposed in [18].
Our derivation of simplified algorithms has been only partly inspired by the
manuscripts cited above. In fact, first of all, it relies on the following specific
methods: a) a set of Np equal weight particles {xj ; j = 0, 1, ..., Np − 1} is
represented through its center of mass
x̄ ,

Np −1
1 X
xj ,
Np j=0

(67)

as already suggested in [17] and [18], when the computation of a message refer(L)
(N )
ring to xl involves the particle-based representation of xl ; b) a set of Np
Gaussian messages {N (x; ηj , Cj ); j = 0, 1, ..., Np − 1}, that refer to a set of Np
equal weight particles, is represented as the Np −component GM
Np −1
1 X
N (x; ηj , Cj ) ,
fGM (x) ,
Np j=0

(68)

22

G. M. Vitetta et al., Marginalized Particle Filtering . . .

and this GM is approximated through its projection onto the Gaussian pdf
fG (x) = N (x; ηG , CG ), where ηG and CG are selected as
ηG ,

Np −1
1 X
ηj
Np j=0

(69)

and
PNp −1
CG = (1/Np ) j=0
Cj
PNp −1
T
T
−ηG (ηG ) + (1/Np ) j=0 ηj (ηj )

(70)

respectively, so that the mean and covariance of fGM (x) (68) are preserved
(L)
(e.g., see [33, Sec. IV]). Secondly, as far as the messages {m
~ in,j (xl )} (21)
are concerned, we do not adopt the approximations proposed in [4], [17] and
[18]. In fact, we focus on the following two cases: case #1 - the messages
(L)
{m
~ in,j (xl )} are all different but, when needed in message passing, are con(L)
(L)
densed in the single message (66), where ηl/(l.−1) and Cl/(l−1) are given by (69)
(L)

(L)

and (70), respectively, with ηj = ηl/(l.−1),j and Cj = Cl/(l−1),j ; b) case #2 (L)

the messages {m
~ in,j (xl )} have different means, but their covariance matrices
(L)
(L)
{Cl/(l−1),j } are all equal (their common value is denoted C̃l/(l−1) in the following). In both cases our simplifications aim at minimizing the overall number of
a) Cholesky decompositions for the generation of the new particle set S(l+1)/l
(L)
(such decompositions are required for the Np matrices {C5,l,j } (47)) and b)
matrix inversions; such inversions are needed to compute: a) the Np matrices
(N )
(N )
{W1,l,j , (C1,l,j )−1 } (required in the evaluation of the weights {wl,j } on the
(L)

basis of (29)); b) the Np matrices {Wl/(l−1),j } (required to evaluate the vec(L)

(L)

(L)

tors {w2,l,j } (41) and the matrices {W2,l,j } (42)); c) the Np matrices {C2,l,j }
(L)

(employed in (47)); c) the Np matrices {C4,l,j } (employed in (64)).
Based on the methods illustrated above, our simplified versions of MPF are
(N )
derived as follows. First of all, in the measurement update for xl , we use a
single covariance matrix in the Gaussian pdf appearing in the RHS of (29); in
other words, the j-th weigth wl,j is computed as


(N )
(N )
wl,j , N yl ; η1,l,j , C1,l ,
(71)
(N )

(N )

where C1,l is evaluated on the basis of (70) (see also (69)), setting ηj = η1,l,j
(N )

(30) and Cj = C1,l,j (31) for any j.
(L)

(N )

Second, in the measurement update for xl , the particle set Sl/(l−1) is con(N )

densed in its center of mass x̄l/(l−1) (see (67)). Consequently, the message

23

G. M. Vitetta et al., Marginalized Particle Filtering . . .
(L)

m
~ 1,j (xl ) (35) is replaced by its particle-independent form




(L)
(L)
= N yl ; B̄l xl + h̄l , Ce
m
~ 1 xl

(72)

(N )

(N )

where B̄l , Bl (x̄l/(l−1) ) and h̄l , hl (x̄l/(l−1) ). This message, similarly as (35),
can be put in the Gaussian form (see (36)-(38))




(L)
(L) (L)
(L)
= N xl ; η1,l , C1,l ,
(73)
m
~ 1 xl
where

(L)

(L) (L)

w1,l , W1,l η1,l = B̄Tl We yl − h̄l
and



(74)


−1
(L)
(L)
W1,l , C1,l
= B̄Tl We B̄.

(75)

(L)

This allows us to replace the message m
~ 2,j (xl ) (40) with its particle-independent
counterpart




(L)
(L) (L)
(L)
= N xl ; η2,l , C2,l ,
(76)
m
~ 2 xl
(L)

(L) (L)

(L)

(L)

(L)

(L)

where w2,l , W2,l η2,l and W2,l , (C2,l )−1 are easily obtained from (41) and
(L)

(L)

(42) replacing a) w1,l,j and W1,l,j with w1,l (74) and W1,l (75), respectively;
(L)

(L)

(L)

(L)

(L)

(L)

b) wl/(l−1),j and Wl/(l−1),j with wl/(l−1) , Wl/(l−1) ηl/(l−1) and Wl/(l−1) ,
(L)

(L)

(L)

(Cl/(l−1) )−1 (ηl/(l.−1) and Cl/(l−1) are given by (69) and (70), respectively, with
(L)

(L)

(L)

ηj = ηl/(l.−1),j and Cj = Cl/(l−1),j ). Note that, since the precision matrix W2,l
(L)

is particle-independent, a single matrix C2,l has to be computed for the next
step.
(N )
(L)
Thirdly, in the time update for xl , the message m
~ 2 (xl ) (76) can be used
(L)
(N )
in place of m
~ 2,j (xl ) (40) in the evaluation of m
~ 5,j (xl+1 ) (see (44) and Fig. 3).
This leads to the particle-dependent message




(N )
(N ) (N )
(N )
m
~ 5,j xl+1 = N xl+1 ; η5,l,j , C5,l,j ,
(77)
(N )

(N )

where η5,l,j and C5,l,j are obtained from (46) and (47), respectively, replac(L)

(L)

(L)

(L)

ing η2,l,j and C2,l,j with η2,l and C2,l , respectively. Then, to simplify the
(N )

(N )

generation of the new particle set S(l+1)/l , the covariance matrices {C5,l,j } are
(N )

(N )

condensed in a single matrix C5,l using (70) (see also (69)) with ηj = η5,l,j =
(N ) (L)

(N )

Al,j η2,l + fl,j

(N )

(N )

and Cj = C5,l,j = Cw

(N )

(N )

(L)

(N )

+ Al,j C2,l (Al,j )T ; consequently,

the particle generation mechanism for xl+1 requires computing the Cholesky de(N )
composition of a single matrix (namely, C5,l ), since it is based on the modified
message




(N )
(N ) (N )
(N )
(78)
m
~ 5 xl+1 = N xl+1 ; η5,l,j , C5,l ,

24

G. M. Vitetta et al., Marginalized Particle Filtering . . .
(N )

which depends on the particle index j through the mean vector η5,l,j only.
(L)

Finally, in the time update for xl , the pdf




(L)
(L)
(L)
(N ) (L)
)
f zl xl
= N zl ; Āl xl , C(N
w
(L)

(79)
(N )

is employed in the evaluation of m
~ 3,j (xl ) through (50), where Āl

(N )

= Al

(N )

(x̄l/l )

(N )

(N )

and x̄l/l denotes the center of mass of the particle set Sl/l (see (67)). Then,
(L)

(L)

the messages m
~ 3,j (xl ) (55) and m
~ 4,j (xl ) (57) can be replaced by




(L)
(L)
(N ) (L)
)
= N zl,j ; Āl xl , C(N
m
~ 3,j xl
w

(80)

and





(L)
(L) (L)
(L)
= N xl ; η4,l,j , C4,l ,
m
~ 4,j xl

(81)

respectively, with (see (58) and (59))

and


T
(L)
(L) (L)
(L)
(N )
(N ) (L)
w4,l,j , W4,l η4,l,j = w2,l + Āl
Ww
zl,j

(82)


−1

T
(L)
(L)
(L)
(N )
(N ) (N )
W4,l , C4,l
= W2,l + Āl
Ww
Āl .

(83)

(L)

(L)

(N )

Substituting now m
~ 4,j (xl ) (81) (in place of m
~ 4,j (xl ) (57)) and m
~ 2,j (xl )
(32) in (60), produces, after some manipulation




(L)
(L)
(L)
(L)
(84)
m
~ 5,j xl+1 = N xl+1 ; η5,l,j , C5,l,j ,
where

(L)

(L) (L)

(L)

η5,l,j , Al,j η4,l,j + fl,j
and
(L)

(L)

(L)

C5,l,j , C(L)
w + Al,j C4,l

(L)



(L)

Al,j

(85)
T

.

(86)

Note that: a) since the precision matrix W4,l (83) is particle-independent, a
(L)

single matrix inversion is needed to evaluate the matrix C4,l appearing in (86);
(L)

(L)

b) in case # 2 the matrix set {C5,l,j } is condensed in a single matrix C5,l (this
(L)

(L)

represents the common value C̃(l+1)/l taken on by all the matrices {C(l+1)/l,j }
(L)

processed in the next recursion); c) the matrix C5,l is evaluated on the basis of
(L)

(L)

(69) and (70), setting ηj = η5,l,j and Cj = C5,l,j for any j.
The new filtering techniques, based on MPF and on the simplified messages
derived above, are called simplified MPF (SMPF) in the following; in particular
the acronyms SMPF1 and SMPF2 are used to refer to case #1 and case #2,
respectively.

25

G. M. Vitetta et al., Marginalized Particle Filtering . . .

6

Message Passing in Iterative Filtering Techniques Inspired by Marginalized Particle Filtering

As already mentioned above, the suboptimality of MPF can related to the fact
that the FG underlying the considered filtering problem is not cycle free. It is
well known that the SPA can also be applied to a factor graph with cycles simply
by following the same message propagation rules; however, generally speaking,
this leads to an “iterative” algorithm with no natural termination (and known
as loopy belief propagation), since its messages are passed multiple times on a
given edge [20], [21]. Despite this, some of the most relevant applications of
the SPA have been developed for systems in which the underlying FG does
have cycles, like the one shown in Fig. 1. In the following we show how a
novel iterative technique can be developed for our filtering problem following
this approach. To begin, we note that our interest in iterative methods is also
(N )
motivated by the possibility of exploiting the pseudo-measurement zl (18); in
(N )
fact, the message mj (zl ) referring to this random vector cannot be computed
in MPF because of the adopted scheduling, but can certainly provide additional
information to refine our statistical knowledge about the nonlinear component
(N )
(L)
xl . This message, similarly as the one referring to10 zl , can be put in a
Gaussian form, that is




(N )
(N ) (N )
(N )
= N zl ; ηz,l,j , Cz,l,j ,
(87)
m
~ j zl
(L)

(L)

(N )

since xl+1 and xl , conditioned on xl , are modelled as jointly Gaussian
random vectors. Let us show now how this message can be computed in an
iterative filtering procedure generalising MPF and how it can exploited in such
(L)
a procedure. First of all, we assume that the message m
~ 5,j (xl+1 ) (62), repre(L)
(N )
(N )
senting the pdf of xl+1 conditioned on xl = xl/l,j , is already available when
(L)

the time update for xl (i.e., the first step of MPF) is accomplished. Then,
(L)
(L)
(N )
given m
~ 2,j (xl ) (40) and m
~ 5,j (xl+1 ) (62), the mean and covariance of zl can
be evaluated as (see (18))
(N )

(L)

(L) (L)

ηz,l,j = η5,l,j − Al,j η2,l,j
10 If

(N)

xl

(N)

(L)

= xl/l,j , the random vector zl

(N)

(N)

(88)
(N)

(13) becomes xl+1 − fl

(N)

(xl/l,j ); consequently,
(L)

adopting the Gaussian model (45) for xl+1 results in a Gaussian model for zl

too.

26

G. M. Vitetta et al., Marginalized Particle Filtering . . .

and


T
(N )
(L)
(L) (L)
(L)
Cz,l,j = C5,l,j + Al,j C2,l,j Al,j

T 
T
(L) (L)
(L)
(L)
,
−Al,j Cx,l,j − Cx,l,j
Al,j

(89)

(L)

respectively, where Cx,l,j denotes the cross covariance matrix for the vectors
(L)

xl

(L)

(N )

and xl+1 (conditioned on xl
(L)

(L)

(L)

(N )

~ 2,j (xl ) (40) and the
= xl/l,j ). Given m
(L)

(N )

(L)

(L) (L)

(L)

conditional pdf f (xl+1 |xl , xl/l,j ) = N (xl+1 ; fl,j + Al,j xl , Cw ) (see (17)),
(L)

(L)

(L)

it is easy to show that Cx,l,j = C2,l,j (Al,j )T (e.g., see [32, Par. 2.3.3, eq.
(2.104)]); consequently, eq. (89) can be rewritten as

T
(N )
(L)
(L) (L)
(L)
Cz,l,j = C5,l,j − Al,j C2,l,j Al,j
.

(90)

Equations (88) and (90) represent the desired result, since they provide a com(N )
plete statistical characterization of the message m
~ j (zl ) (87). In principle,
(L)
this message could be exploited in a similar way as that adopted for m
~ j (zl )
(N )
(53); this approach would lead to draw a set of Np samples {zl,j } from the
Gaussian function appearing in the RHS of (87) and to process the resulting
pseudo-measurements to generate a new weight for each particle of the set Sl/l .
However, our computer simulations have shown that this approach is outperformed by more refined method illustrated in the following. In practice, in the
(N )
iterative filtering method we propose the message m
~ j (zl ) (87), similarly as
(L)
m
~ j (zl ) (53), is employed to evaluate the new message11


 

 Z
(N )
(N )
(N )
(N )
(N )
(91)
f zl
xl/l,j dzl ,
= m
~ j zl
m
~ 3,j xl
(N )

(L)

which represents for xl the counterpart of the message m
~ 3,j (xl ) (50). Sub(N ) (N )
(N )
(N )
stituting (87) and the expression of f (zl |xl ) (given xl
= xl/l,j ) in the
RHS of the last expression gives12 (see the Appendix)


(N )
m
~ 3,j xl
 
T
(N )
(N )
(N ) (N )
= D3,l,j · exp 21
η3,l,j W3,l,j η3,l,j
(92)


T

T
(N )
(N ) (N )
(L)
(L) (L)
− ηz,l,j Wz,l,j ηz,l,j − fl,j
Ww fl,j
, pl,j ,
11 Note

(N)

that the following message represents the correlation between the pdf m
~ j (zl

evaluated on the basis of the definition of

(N)
zl

)

(18) and the pdf originating from the fact that
(L)

(N)

(L)

this quantity is expected to equal the random vector fl (xl/l,j ) + wl . For this reason, it
expresses the degree of similarity between these two functions.
12 In our computer simulations the factor D (N) appearing in this weight has been always
3,l,j
neglected, since it negligibly influences estimation accuracy.

27

G. M. Vitetta et al., Marginalized Particle Filtering . . .

where

−1
(N )
(N )
(N )
(L)
W3,l,j , C3,l,j
= Wz,l,j + Ww
,
(N )

(N )

(N )

(N )

(93)

(L)

(L)
w3,l,j , W3,l,j η3,l,j = wz,l,j + Ww
fl,j ,
(N )

(N )

(N )

(N )

(N )

(L)

(94)

(L)

Wz,l,j , (Cz,l,j )−1 , wz,l,j , Wz,l,j ηz,l,j , Ww , [Cw ]−1 ,

(N )

(N )

(L)

h

i−DL /2
(N )
(N )
D3,l,j , det C̃l,j

(95)
(N )

and C̃l,j , Cz,l,j + Cw . Then, the new message m
~ 3,j (xl
similarly as

) (92) is exploited,

(L)
m
~ 3,j (xl )

(N )

(55), to generate the message (see (51))






(N )
(N )
(N )
,
m
~ 3,j xl
=m
~ 2,j xl
m
~ 4,j xl

(96)

where m
~ 2,j (xl ) is expressed by (28) (i.e., it is the message emerging from the
(N )
measurement update for xl in the absence of resampling). This produces (see
(92))




(N )
(N )
(N )
(97)
= Wl,j δ xl − xl/(l−1),j ,
m
~ 4,j xl
where

Wl,j , wl,j · pl,j

(98)

represents the new weight for the j-th particle; such a weight accounts for both
(N )
the (real) measurement yl and the pseudo-measurement zl,j . Resampling with
(N )

(N )

replacement can now be accomplished for the set Sl/(l−1) , {xl/(l−1),j } on
the basis of the more refined weights {Wl,j } (98); if this is done, the message
(N )
(N )
m
~ 4,j (xl ) takes on the same form as m
~ 2,j (xl ) (32), i.e. it can be expressed
as




(N )
(N )
(N )
(99)
= δ xl − x̃l/l,j ,
m
~ 4,j xl
(N )

with j = 0, 1, ..., Np − 1. Note that, generally speaking, the particle set S̃l/l ,
(N )

(N )

(N )

{x̃l/l,j } produced by resampling in this case is different from Sl/l , {xl/l,j }
(i.e., from the one obtained with MPF), even if both of them originate from the
(N )
same set Sl/(l−1) ; this is due to the fact that the weights {Wl,j } (98) may be
substantially different from the MPF weights {wl,j } because of the factor pl,j .
(N )
(N )
Finally, the new message m
~ 4,j (xl ) (99) is used in place of m
~ 2,j (xl ) in the
(L)
RHS of (43) for the evaluation of the message m
~ 5,j (xl+1 ).
As already state above, our previous derivations rely on the assumption that
(L)
(L)
the message m
~ 5,j (xl+1 ) (62) is available when the time update for xl is carried
out; unluckily, this message becomes available only in the last step of MPF.
However, if MPF is generalised in a way that Nit > 1 iterations (i.e., message

28

G. M. Vitetta et al., Marginalized Particle Filtering . . .

passes) are carried out within the same recursion, in the k-th iteration (with k =
(N )
2, 3, ..., Nit ) the message mj (zl ) (87) can be really evaluated exploiting the
(L)
message m
~ 5,j (xl+1 ) computed in the previous iteration. These considerations
lead, in a natural fashion, to the development of the message passing illustrated
in Fig. 4, which describes the message flow occurring in the k-th iteration
of a new filtering technique, generalising MPF and called turbo filtering (TF)
in the following; note that the superscripts (k) or (k − 1) have been added
to all the messages flowing in the considered graph to identify the iteration
in which they are generated and that the grey circle appearing in the figure
represents a unit delay cell. The processing tasks accomplished by TF can be
summarized as follows. The first part of this technique can be considered as a
(N )
(L)
form of initialization, in which the messages {m
~ 2,j (xl )} and {m
~ 2,j (xl )} are
computed; however, unlike MPF, resampling is not accomplished in the last part
(N )
(N )
of the time update for xl
(so that the messages {m
~ 2,j (xl )} are expressed
(N )
by (28) instead of (32) and refer to the particle set Sl/(l−1) ). Moreover, as
(N )

(L)

shown in Fig. 4, the messages {m
~ 2,j (xl )} and {m
~ 2,j (xl )} emerging from
the first part of TF remain unchanged within the l-th recursion, since they
(N )
(L)
represent the a priori information available about xl and xl , respectively;
consequently, like in any turbo processing method, these information are made
available to all the iterations carried out within each recursion. In the second
part of TF Nit -iterations are accomplished with the aim of progressively refining
(N )
the set S(l+1)/l (the version of this set generated in the k-th iteration is denoted
(k)

(N )

(L)

~ 5,j (xl+1 )}. To achieve
S(l+1)/l [k]) and the associated Gaussian messages {m
these result, in the k-th iteration (with k = 1, 2, ..., Nit ) the ordered computation
(k) (N )
(k) (N )
~ 3,j (xl ) (92)
of the following messages is accomplished: m
~ j (zl ) (87), m
(k)

(N )

(conveying the weight pl,j [k]), m
~ 4,j (xl
(k)

(N )

(L)

(k)

(L)

) (97) (conveying the weight Wl,j [k]),
(k)

(L)

(k)

(L)

~ 5,j (xl+1 )
~ 4,j (xl ) (57) and m
~ j (zl ) (53), m
~ 3,j (xl ) (55), m
m
~ 5,j (xl+1 ) (45), m
13
(62). Moreover, in the k-th iteration resampling is accomplished on the basis of
the particle weights {Wl,j [k]}; generally speaking, this results in a new particle
(N )
(N )
set denoted Sl/l [k] = {xl/(l−1),j [k], j = 0, 1, ..., Np − 1} (which is always a
(N )

subset of Sl/(l−1) ). It is also important to point out that:
(0)

(N )

• In the first iteration (i.e, for k = 1) the messages {m
~ 5,j (xl+1 )} are un(L)

13 Note that, after carrying out resampling, the set of messages {m
~ 2,j (xl )} needs to be
properly reordered, since the messages associated with the discarded particles are not preserved. This modifies the set of particles available in the next iteration and, consequently,
the set of weights {wl,j } associated with them (these weights need to be renormalized after
any change). In the following the notation {wl,j [k]} is adopted to denote the set of weights
employed in the k-th iteration for the evaluation of the overall weights Wl,j [k] according to
(98).

29

G. M. Vitetta et al., Marginalized Particle Filtering . . .
(1)

(N )

(1)

(N )

~ 3,j (xl ) =
defined, so that {m
~ j (zl ) = 1} (and, consequently, {m
1}) must be assumed; moreover, resampling is not accomplished (i.e.,
(N )
(N )
(N )
Sl/l [1] = Sl/l = Sl/(l−1) ), since the weights pl,j appearing in the overall
weight Wl,j (98) become available in the following iterations.
− (x(L) ) = ←
− (x(N ) ) = 1 is assumed for
• In each iteration the equality ←
m
m
he l+1
he l+1
(L)
the two messages entering the FG along the half edges associated with xl+1
(N )
and xl+1 (see Fig. 4), since no information comes from the next recursion.
For this reason, at the end of the last iteration (i.e., for k = Nit ), the
output messages (i.e., the input messages feeding the (l + 1)-th recursion)
are evaluated as (see (48)-(49) and (62))






(N )
(N )
(Nit )
(N )
(N )
− (x(N ) ) = m
x
(100)
m
~
~ 5,jit xl+1 ←
m
~ out,j xl+1 = m
he l+1
5,j
l+1
and







(L)
(L)
(N )
(L)
(N )
− (x(L) ) = m
m
~ 5,jit xl+1 ,
~ 5,jit xl+1 ←
m
~ out,j xl+1 = m
he l+1

(101)

with j = 0, 1, ..., Np − 1.

Another relevant issue concerns the interpretation of the processing tasks
accomplished in the TF technique. In fact, our derivations show that, at the
end of the k-th iteration, the a posteriori statistical information about the j-th
(k) (N )
(N )
(N )
~ 4,j (xl ) (97), which
particle xl/l,j [k] of Sl/l [k] is provided by the message m
conveys the weight (see (98))
(a)

Wl,j [k] = wl,j [k] · pl,j [k] · wl,j ,

(102)

(a)

where wl,j denotes the a priori information available at the beginning of the l(a)

th recursion for the j-th particle ( in our derivation wl,j = 1 has been assumed,
(a)

in place of wl,j = 1/Np , to simplify the notation; see (22)), pl,j [k] is the weight
(N )

(k)

(N )

originating from m
~ j (zl ) (87) and conveyed by m
~ 3,j (xl ) (92), and wl,j [k] is
the weight computed on the basis of the available measurement yl . Taking the
natural logarithm of both sides of (102) produces
(a)

(y)

(z)

Ll,j [k] = Ll,j + Ll,j [k] + Ll,j [k],
(y)

(z)

(103)
(a)

where Ll,j [k] , ln(Wl,j [k]), Ll,j [k] , ln(wl,j [k]), Ll,j , ln(pl,j [k]) and Ll,j ,
(a)

ln(wl,j ). The last equation has exactly the same structure as the well known
formula (see [35, Sec. 10.5, p. 450, eq. (19.15)] or [36, Par. II.C, p. 432, eq.
(20)])
L (uj |y) = L (uj ) + Lc (yj ) + Le (uj )
(104)

30

G. M. Vitetta et al., Marginalized Particle Filtering . . .


min , j ( x (l L ) )

x (l L )

yl


m1, j ( x (l L ) )

f y ,/ x( L ) ,x( N )

l +1

=

=
l




m2, j ( x l( L ) )

=

l

m4,( k j) ( x (l L ) )

mhe ( xl(+L1) ) = 1
f x ( L ) / x ( L ) ,x ( N )


m2,( k j) ( x (l L ) )


min , j ( x (l L ) )

=

l


m2, j ( x l( L ) )

m5,( kj) ( x (l +L1) )


f z( L ) / x ( L )


m (jk ) ( z l( L ) )


mout , j ( x l(+L1) )

l

=


m3,( kj) ( x (l L ) )

l


min , j ( x (l N ) )

l

x (l +L1)

l


m2, j ( x l( N ) )

z (l L )

x (l N ) , x l( +N1) → z l( L )

m5,( kj) ( x (l +N1) )


m2, j ( x (l N ) )

m2, j ( x l( N ) )

=


m2, j ( x l( L ) )

=

m2, j ( x


(N )
l

)


m (jk ) ( z l( N ) )

f y ,/ x( L ) ,x( N )
l


m3,( kj) ( x l( N ) )

l

=
=

min , j ( x (l N ) )

=

 ( k −1) ( L )
z (l N ) m
( xl +1 )
5, j

f z( N ) / x ( N )
l

l

yl


min , j ( x (l N ) )


m2, j ( x l( N ) )


m2,( k j) ( x l( L ) )

m5,( kj) ( x l(+N1) )

=

=


m2, j ( x l( N ) )
x (l N )


m5,( kj) ( x (l +L1) )

x (l L ) , x (l +L1) → z (l N )


min , j ( x (l L ) )

l


m2, j ( x l( L ) )


m4,( k j) ( x l( N ) )

f x ( N ) / x ( N ) ,x ( L )
l +1
l
l

mhe ( x l(+N1) ) = 1


mout , j ( x (l +N1) )
x

(N )
l +1

Figure 4: Message passing over the FG of Fig. 1 for the proposed TF technique.
All the quantities appearing in this figure refer to the k-th iteration of the l-th
recursion; the messages available at the beginning of the considered iteration are
indicated by red arrows, those entering the graph by blue arrows, those leaving
it at the end of the last iteration by black arrows and those computed within
the considered iteration by green arrows.

31

G. M. Vitetta et al., Marginalized Particle Filtering . . .

expressing of the log-likelihood ratio (LLR) available for the j-th information
bit uj at the output of a soft-input soft-output channel decoder operating over
an additive white Gaussian noise (AWGN) channel and fed by: a) the channel
output vector y (whose j-th element yj is generated by the communications
channel in response to a channel symbol conveying uj and is processed to produce the so-called channel LLR Lc (yj )); b) the a priori LLR L (uj ) about uj ;
c) the extrinsic LLR Le (uj ), i.e. a form of soft information available about
uj , but intrinsically not influenced by such a bit (in turbo decoding of concatenated channel codes extrinsic infomation is generated by another channel
decoder with which soft information is exchanged with the aim of progressively
refining data estimates). This correspondence is not only formal, since in eqs.
(103) and (104) terms playing similar roles can be easily identified. For instance,
(y)
(a)
the term Ll,j [k] (Ll,j ) in (103) provides the same kind of information as Lc (yj )
(L (uj )), since these are both related to the noisy data (a priori information)
available about the quantities to be estimated (the system state in one case, an
(z)
information bit the in the other one). What about the term Ll,j [k] appearing
in the RHS of (103)? The link we have established between (103) and (104)
unavoidably leads to the conclusion that such a term should represent the counterpart of the quantity Le (uj ) appearing in (104), i.e. the so called extrinsic
(N )
information (in other words, that part of the information available about xl
(N )
and not intrinsically influenced by xl itself). This interpretation is confirmed
(z)
by the fact that Ll,j [k] is computed on the basis of the statistical knowledge
(L)

(L)

available about xl and xl+1 (see (88) and (90)), which, thanks to (2), does
(N )
provide useful information about xl . The theory of turbo decoding of channel
codes shows that, generally speaking, extrinsic information originates from code
constraints. In our scenario, a similar interpretation can be also provided for
(z)
(L)
Ll,j [k], since xl+1 can be seen as the noisy output of a communication channel,
(L)

(L)

affected by the bias fl,j and the additive noise wl , and over which the code(L) (L)

word Al,j xl

of a rate-1 block code is transmitted in response to the message

(L)
xl .

(z)

These considerations show that, in evaluating Ll,j [k], we are actually
exploiting a sort of ‘code’ constraints, which are mathematically expressed by
(2). The reader can easily verify that a similar interpretation can be provided
(L)
for m
~ j (zl ) (53), which represents the extrinsic information component14 con(L)
(L)
tained in m
~ 4,j (xl ) (57) (conveying our a posteriori information about xl );
(L)
the other two components are represented by the message m
~ 2,j (xl ) (40) (repre14 In

(L)

practice, the mechanism employed to generate zl,j is based on the state update equation (3).

G. M. Vitetta et al., Marginalized Particle Filtering . . .
(L)

32
(L)

senting the measurement information about xl ) and the message m
~ in,j (xl )
(L)
(21) (corresponding to our a priori information about xl ). Consequently, TF
can be seen, in the domain of Bayesian filtering techniques, as the counterpart
of turbo decoding of concatenated codes; this parallelism can be exploited to
provide further insights into iterative filtering techniques. For instance, it is
well known that, in turbo decoding of concatenated channel codes, the extrinsic
information generated by soft decoders become more and more correlated as iterations evolve; this entails that diminishing benefits are provided by additional
iterations. This phenomenon should be observed in TF too for similar reasons
(N )
(N )
and can be motivated by rewriting ηz,l,j (88) and Cz,l,j (90) as
(N )

(L)

(L)

ηz,l,j = fl,j + Al,j
and
(N )

(L)

Cz,l,j = C(L)
w + Al,j

h

h
i
(L)
(L)
η4,l,j − η2,l,j

(L)

(L)

C4,l,j − C2,l,j

i
T
(L)
,
Al,j

(105)
(106)

respectively (thanks to (85) and (86), respectively). In fact, the last two equa(N )
(N )
tions show that the vector ηz,l,j and the matrix Cz,l,j are influenced by the
difference between the statistical information (expressed by a mean vector and
(L)
a covariance matrix) available about xl before processing the pseudomeasure(L)
ment zl and those available after this task has been carried out. In other
(L)
(N )
words, the extrinsic information provided by zl influences zl and viceversa.
Finally, it is important to point out that the proposed analogy between turbo
filtering and turbo decoding suggests the potential limits of the TF technique
(and of any other iterative filtering method relying on the developed FG). In
fact, it is well known that turbo decoding methods do not provide real benefits
below a certain signal-to-noise ratio, i.e. when the quality of the received signal
is so poor that the transmitted coded sequence cannot be recovered. A similar
phenomenon is expected occur with TF too; consequently, this filtering method
could not outperform MPF in the presence of strong measurement noise and/or
fast dynamics affecting the considered SSM.

7

Numerical Results

In this Section MPF and the related filtering methods developed in this manuscript
are compared in terms of accuracy and computational load for a specific CLG
system, characterized by DL = 3, DN = 1 (so that D = 4) and P = 2. The
structure of the considered system has been partly inspired by the example

G. M. Vitetta et al., Marginalized Particle Filtering . . .

proposed by Schön

0.8
(L)
xl+1 =  0
0
and

(L)

in [37] and is characterized by: a) the state models



(N )
cos(xl )
0.2
0


(L)
(L)
)
0.7 −0.2  xl +  − sin(x(N
)  + wl
l
(N )
0.2 0.7
0.5 sin(2xl )



(N )
(N )
(L)
(N )
xl+1 = arctan xl
+ (0.9 0 0) xl + wl
(L)

(N )

33

(107)

(108)

(N )

∼ N (0, (σw )2 I3 ), wl ∼ N (0, (σw )2 ; b) the measurement model



 ! 

(N ) 2
(N )
0 0 0
0.1 xl
· sgn xl
(L)
yl =
+
xl + el
(109)
1 −1 1
0

with wl

with el ∼ N (0, (σe )2 I2 ). Note that the state equation (107), unlike itscounter
(N )
(L)
(N )
part proposed in [37], depends on xl (i.e., it contains a function fl
xl+1 6=
(N )

03 in its RHS), so that TF, which relies on the availability of the vector zl
(18), can be employed for this system.
In our computer simulations the root mean square error (RMSE) has been
evaluated to compare the accuracy of the state estimates generated by different
filtering techniques. More specifically, for each technique two RMSEs have been
computed, one (denoted RM SEL (alg), where ‘alg’ denotes the algorithm this
parameter refers to) representing the square root of the average mean square
(L)
error (MSE) evaluated for the three elements of xl , the other one (denoted
(N )
RM SEN (alg)) referring to the (monodimensional) nonlinear component xl ;
this distinction is important since, as shown by our simulation results, the es(L)
timation accuracy for xl can be quite different from (and is usually smaller
(N )
than) that referring to xl .
As far the assessment of the computational requirements of the investigated
filtering techniques is concerned, MPF (for which an accurate analysis of its
computational complexity is available in [16]) has been taken as a baseline. For
this reason, our comparisons between the considered filtering techniques are
based on the evaluation of a single parameter, denoted ∆c (alg) and representing the percentage variation in computation time of the considered algorithm
(denoted ‘alg’) with respect to MPF (operating with the same parameters and,
in particular, with the same Np as alg).
Moreover, in our computer simulations, the following choices have been
(N )
(L)
made: a) σw = σw = 5 · 10−3 has been selected for the standard devia(L)
(N )
tions of the process noises {wk } and {wk }, unless differently stated; b) the
so called jittering technique [34] has been employed to mitigate the so called
depletion problem in the generation of new particles.

34

G. M. Vitetta et al., Marginalized Particle Filtering . . .

Some results illustrating the dependence of RM SEL and RM SEN on the
number of particles (Np ) for MPF and SMPF#1 are illustrated in Fig. 5 (in this
and in the following figures simulation results are identified by markers, whereas
continuous lines are drawn to ease reading); σe = 10−2 and Np ∈ [50, 300] have
been selected in this case. From these results the following conclusions can be
easily inferred for the considered system:
1. A negligible improvement in the estimation accuracy of both MPF and
SMPF#1 is achieved if the value of Np exceeds 200.
2. A significant gap between RM SEL (MPF) and RM SEN (MPF) exists; this
(L)
is motivated by the fact that, in MPF, the estimation of xl relies on both
(L)
the real measurement yl (109) and the pseudo-measurements {zl,j } (54),
(N )

whereas the estimation of xl

benefits from yl only.

3. The gap between RM SEL (SMPF#1) and RM SEN (SMPF#1) is much
smaller than that observed for MPF, even if the pseudo-measurements
(L)
{zl,j } are also exploited by SMPF#1. This reduction in the RMSE gap
(L)

can be related to the degradation in the estimation accuracy of xl ; for instance, for Np = 200, RM SEL (SMPF#1) is about twice RM SEL (MPF)
(on the contrary, RM SEN (SMPF#1) ∼
= 1.09· RM SEN (MPF)).
Our numerical results have also evidenced that: a) SMPF#1 requires a substantially smaller computational effort than MPF, since ∆c (SMPF#1) approximately ranges in the interval [−62%, −58%] for the considered values of Np ;
b) despite the above mentioned degradation in estimation accuracy (see point
3.), SMPF#1 does not suffer from tracking losses in the considered scenario;
c) SMPF#2 accuracy is almost indentical that of SMPF#1; d) ∆c (SMPF#2)
is approximately lower than ∆c (SMPF#1) by 6% in the considered range for
Np and, consequently, achieves a better complexity-performance tradeoff than
SMPF#1. All this suggests that simplified MPF techniques can be really developed without incurring the serious technical problems that affect the techniques
proposed in [17] (tracking losses and poor RMSE performance) and in [18] (partitioning and update of the particle set into groups within each recursion; see
Section 5).
The dependence of RM SEL and RM SEN on σe (i.e., on the intensity of
the noise affecting the available measurements) has been also analysed for MPF
and SMPF#1. Some results are shown in Fig. 6; Np = 500 and σe ∈ [10−3 , 5 ·
10−2 ] have been selected in this case. These results show that the gap between RM SEL (MPF) and RM SEL (SMPF#1) slightly increases as σe becomes

35

G. M. Vitetta et al., Marginalized Particle Filtering . . .

5

4

SMPF1

RMSE

3

-2

2x10

RMSEL
RMSEN

MPF

0.01
50

100

150

200

250

300

Np
Figure 5: RMSE performance versus Np for the linear component (RM SEL )
and the nonlinear component (RM SEN ) of the state xl for the system described
(L)
by eqs. (107)-(109). MPF and SMPF#1 are considered; in both cases σw =
(N )
σw = 5 · 10−3 and σe = 10−2 have been selected.

36

G. M. Vitetta et al., Marginalized Particle Filtering . . .

6

SMPF1

5
4
3

RMSE

2

MPF
0.01

9
8
7

RMSEL
RMSEN

6
5

10

20

30

40

-3

50x10

σe
Figure 6: RMSE performance versus σe for the linear component (RM SEL ) and
the nonlinear component (RM SEN ) of the state xl for the system described by
(L)
eqs. (107)-(109). MPF and SMPF#1 are considered; in both cases σw =
(N )
σw = 5 · 10−3 and Np = 200 have been selected.
smaller; the opposite occurs for RM SEN (MPF) and RM SEN (SMPF#1). These
(L)
results can be related again to the fact that the pseudo-measurements {zl,j }
(L)

(54) really play a role in the MPF estimation of xl and that the quality of the
information conveyed by these pseudo-measurements indirectly improves as the
real measurements become less noisy.
Some results illustrating the dependence of RM SEL and RM SEN on the
number of particles (Np ) for MPF and TF are illustrated in Fig. 7; σe = 10−2 ,
Np ∈ [1, 300] and Nit = 2 for TF, and Np ∈ [30, 300] for MPF have been chosen
in this case. From these results it is easily inferred that:
1. TF outperforms MPF in tems of both RM SEL and RM SEN ; for instance, RM SEL (MPF) ∼
= 1.71· RM SEL (TF) and RM SEN (MPF) ∼
=
2.86· RM SEN (TF) if Np = 200 is selected.

G. M. Vitetta et al., Marginalized Particle Filtering . . .

37

2. The gap between RM SEL (MPF) and RM SEN (MPF) is substantially
larger than the corresponding gap for TF (in particular, RM SEN (MPF)−
RM SEL (MPF) ∼
= 13.4·[RM SEN (TF)−RM SEL (TF)]; this can be easily
(L)
related to the fact that, unlike MPF, in TF the estimation of both xl
(N )
and xl benefits from the availability of pseudo-measurements.
3. The performance gap between TF and MPF increases as Np gets smaller.
Moreover, TF accuracy starts quickly degrading when Np drops below
11, whereas the same phenomenon starts for Np ∼
= 30 with MPF. Note,
however, that different phenomena occur with MPF and TF when the
particle set is small. In fact, our computer simulations have shown that
MPF suffers from frequent tracking losses for Np < 30 (since it is unable
(N )
to generate a reliable representation of xl
when a limited number of
particles is available); on the contrary, TF does not suffer from the same
problem even if a very small particle set is used.
We believe that last result is really important and can be motivated as
follows. The TF technique, through its feedback mechanism, makes a substantially more efficient use of the available particles than MPF; this results in an
appreciable improvement of both stability and accuracy of state estimation.
Our simulation results have also shown that: a) in the considered scenario a
negligible improvement in estimation accuracy is obtained if Nit > 2 is selected;
b) ∆c (TF) approximately ranges in the interval [70%, 80%] for the considered
values of Np and, consequently, requires a substantially larger computational
effort than MPF if these algorithms operate with the same number of particles. In practice, however, as evidenced by the results shown in Fig. 7, TF
can reliably operate with a very small particle set and, consequently, it outperforms MPF in terms of both performance and complexity if Np is properly
selected. For instance, in the considered scenario TF with Np = 20 particles
achieves a substally better accuracy than MPF with Np = 40 particles, even
if, as evidenced by our computer simulations, they approximately require the
same computation time. It is not difficult to show that similar considerations
hold if SMPF#1 and SMPF#2 are considered in place of MPF. Therefore, our
results suggest that the real key to the complexity reduction of the filtering techniques relying on the FG shown in Fig. 2 is not provided by the approximations
adopted for the MPF processing tasks in Section 5, but by the exploitation of
(N )
the new pseudo-measurement zl (18). We should never forget, however, that
TF cannot

 be adopted for all CLG systems; in fact, it cannot be employed if
(L)
(N )
fl
xl+1 = 0DL in the RHS of (2).

38

G. M. Vitetta et al., Marginalized Particle Filtering . . .

5
4

RMSEL
RMSEN

MPF

3

RMSE

2

0.01
9
8
7

TF

6
5

0

50

100

150
Np

200

250

300

Figure 7: RMSE performance versus Np for the linear component (RM SEL )
and the nonlinear component (RM SEN ) of the state xl for the system described
(L)
(N )
by eqs. (107)-(109). MPF and TF are considered; in both cases σw = σw =
−3
−2
5 · 10 and σe = 10 have been selected.

39

G. M. Vitetta et al., Marginalized Particle Filtering . . .

8
7
6
5

RMSEL
RMSEN

σw = 5 10-3
σw = 10-3

MPF

4
3

RMSE

2

0.01
8
7
6
5

TF

4
3
-2

3

2x10

4

σe
Figure 8: RMSE performance versus σe for the linear component (RM SEL ) and
the nonlinear component (RM SEN ) of the state xl for the system described by
eqs. (107)-(109). MPF and TF are considered; σe = 10−2 and Np = 200 are
(L)
assumed in both cases. As far system noise is concerned, the cases σw =
(N )
(L)
(N )
σw = 5 · 10−3 and σw = σw = 10−1 are taken into consideration.
Finally, the dependence of RM SEL and RM SEN on σe has been assessed for
TF and compared with that characterizing MPF. Some numerical are illustrated
in Fig. 8. In this case, we have selected σe ∈ [1.5 · 10−2 , 5 · 10−2 ] , Np = 200
(L)
(N )
and Nit = 2 for TF, and Np = 200 for MPF; moreover, σw = σw = 5 · 10−3
(L)
(N )
and σw = σw = 10−3 have been considered to analyse the dependence of
the performance gap between TF and MPF on the intensity of process noise
(and, consequently, on system dynamics). These results show that: a) the
performance gap between TF and MPF undergoes small changes if σe varies
in the considered interval; b) on the contrary, a substantial change in this gap
(L)
(N )
is obtained if σw and σw are reduced from 5 · 10−3 to 10−3 . This suggests
that measurement noise and process noise can have different impacts on the
performance gain provided by TF over MPF.

5

G. M. Vitetta et al., Marginalized Particle Filtering . . .

8

40

Conclusions

In this manuscript a FG approach has been employed to analyse the filtering
problem for mixed linear/nolinear models. This has allowed us to: a) prove
that this problem involves a FG which is not cycle free; b) provide a new interpretation of MPF as a forward only message passing algorithm over a specific
FG; c) develop novel filtering algorithms for simplifying or generalising it. In
particular, an important iterative filtering technique, dubbed turbo filtering, has
been devised and its relation with the turbo decoding techniques for concatenated channel codes has been analysed in detail. All the considered filtering
techniques have been compared in terms of both accuracy and computational
requirements for a specific CLG system. The most interesting result emerging from our computer simulations is represented by the clear superiority of
turbo filtering over marginalized particle filtering. In fact, the former technique, through the exploitation of new pseudo-measurements, can achieve a
better accuracy than the latter one at an appreciably smaller computational
load. Our ongoing research activities in this area include the development of
other related filtering techniques and the application of turbo filtering to specific
state estimation problems.

A

Appendix

Given the pdfs f1 (y) , N (y; η1 , C1 ) and f2 (y) , N (y; η1 , C2 ) for the N dimensional vector y, we are interested in evaluating the correlation between
these two functions, i.e. the quantity
Z
c1,2 , f1 (y) · f2 (y) dy
(110)
Substituting the expressions of f1 (y) and f2 (y) in the RHS of the last equation
produces, after some manipulation,



1 T
T
T
(111)
η Wη − η1 W1 η1 − η2 W2 η2
c1,2 = D exp
2
−1
where W1 , C−1
1 , W 2 , C2 ,

W = W1 + W2 ,

(112)

Wη = W1 η1 + W2 η2

(113)

D = (2π det [C1 + C2 ])−N/2 .

(114)

and

G. M. Vitetta et al., Marginalized Particle Filtering . . .

41

References
[1] M. S. Arulampalam, S. Maskell, N. Gordon and T. Clapp, “A Tutorial on
Particle Filters for Online Nonlinear/Non-Gaussian Bayesian Tracking”,
IEEE Trans. Sig. Proc., vol. 50, no. 2, pp. 174-188, Feb. 2002.
[2] F. E. Daum, “Exact Finite-Dimensional Nonlinear Filters”, IEEE Tran.
Aut. Contr., vol. 31, no. 7, pp. 616-622, July 1986.
[3] S. Mazuelas, Y. Shen and M. Z. Win, “Belief Condensation Filtering”,
IEEE Trans. Sig. Proc., vol. 61, no. 18, pp. 4403-4415, Sept. 2013.
[4] V. Smidl and A. Quinn, “Variational Bayesian Filtering”, IEEE Trans. Sig.
Proc., vol. 56, no. 10, pp. 5020-5030, Oct. 2008.
[5] M. S̆imandl, J. Královeca and T. Söderströmc, “Advanced Point-Mass
Method for Nonlinear State Estimation”, Automatica, vol. 42, pp. 11331145, 2006.
[6] B. Anderson and J. Moore, Optimal Filtering, Englewood Cliffs, NJ,
Prentice-Hall, 1979.
[7] S. J. Julier and J. K. Uhlmann, “Unscented Filtering and Nonlinear Estimation”, IEEE Proc., vol. 92, no. 3, pp. 401-422, Mar. 2004.
[8] A. Doucet, J. F. G. de Freitas and N. J. Gordon, “An Introduction to Sequential Monte Carlo methods,” in Sequential Monte Carlo Methods
in Practice, A. Doucet, J. F. G. de Freitas, and N. J. Gordon, Eds. New
York: Springer-Verlag, 2001.
[9] A. Doucet, S. Godsill and C. Andrieu, “On Sequential Monte Carlo Sampling Methods for Bayesian Filtering”, Statist. Comput., vol. 10, no. 3, pp.
197-208, 2000.
[10] F. Gustafsson, “Particle Filter Theory and Practice with Positioning Applications”, IEEE Aerosp. and Electr. Syst. Mag., vol. 25, no. 7, pp. 53-82,
July 2010.
[11] F. Gustafsson, F. Gunnarsson, N. Bergman, U. Forssell, J. Jansson, R.
Karlsson and P. Nordlund, “Particle Filters for Positioning, Navigation,
and Tracking”, IEEE Trans. Sig. Proc., vol. 50, 425-435, 2002.

G. M. Vitetta et al., Marginalized Particle Filtering . . .

42

[12] P. J. Nordlund and F. Gustafsson, “Marginalized Particle Filter for Accurate and Reliable Terrain-Aided Navigation”, IEEE Trans. on Aerosp. and
Elec. Syst., vol. 45, no. 4, pp. 1385-1399, Oct. 2009.
[13] R. Bucy and K. Senne, “Digital Synthesis of Non-Linear Filters”, Automatica, vol. 7, no. 3, pp. 287-298, 1971.
[14] F. Daum and J. Huang, “Curse of Dimensionality and Particle Filters”,
Proc. IEEE Aerosp. Conf., vol. 4, pp. 1979-1993, March 2003.
[15] T. Schön, F. Gustafsson, P.-J. Nordlund, “Marginalized Particle Filters for
Mixed Linear/Nonlinear State-Space Models”, IEEE Trans. Sig. Proc., vol.
53, no. 7, pp. 2279-2289, July 2005.
[16] R. Karlsson, T. Schön, F. Gustafsson, “Complexity Analysis of the
Marginalized Particle Filter”, IEEE Trans. Sig. Proc., vol. 53, no. 11, pp.
4408-4411, Nov. 2005.
[17] F. Mustiere, M. Bolic and M. Bouchard, “A Modified Rao-Blackwellised
Particle Filter”, Proc. of the 2006 IEEE Int. Conf. on Ac., Sp. and Sig.
Proc. (ICASSP 2006), vol. 3, 14-19 May 2006.
[18] T. Lu, M. F. Bugallo and P. M. Djuric, “Simplified Marginalized Particle
Filtering for Tracking Multimodal Posteriors”, Proc. IEEE/SP 14th Workshop on Stat. Sig. Proc. (SSP ’07), pp. 269-273, Madison, WI (USA), 2007.
[19] F. Lindsten, P. Bunch, S. Särkkä, T. B. Schön and S. J. Godsill, “RaoBlackwellized Particle Smoothers for Conditionally Linear Gaussian Models”, IEEE J. Sel. Topics in Sig. Proc., vol. 10, no. 2, pp. 353-365, March
2016.
[20] H.-A. Loeliger, J. Dauwels, Junli Hu, S. Korl, Li Ping, F. R. Kschischang,
“The Factor Graph Approach to Model-Based Signal Processing”, IEEE
Proc., vol. 95, no. 6, pp. 1295-1322, June 2007.
[21] F. R. Kschischang, B. Frey, and H. Loeliger, “Factor Graphs and the SumProduct Algorithm”, IEEE Trans. Inf. Theory, vol. 41, no. 2, pp. 498-519,
Feb. 2001.
[22] C. Berrou and A. Glavieux, “Near Optimum Error Correcting Coding and
Decoding: Turbo-Codes”, IEEE Trans. Commun., vol. 44, no. 10, pp. 1261
- 1271, Oct. 1996.

G. M. Vitetta et al., Marginalized Particle Filtering . . .

43

[23] S. Benedetto, D. Divsalar, G. Montorsi and F. Pollara, “Serial Concatenation of Interleaved Codes: Performance Analysis, Design, and Iterative
Decoding”, IEEE Trans. Inf. Theory, vol. 44, no. 3, pp. 909-926, May 1998.
[24] R. Koetter, A. C. Singer and M. Tüchler, “Turbo Equalization”, IEEE Sig.
Proc. Mag., vol. 21, no. 1, pp. 67-80, Jan. 2004.
[25] J. Hagenauer, “The Turbo Principle: Tutorial Introduction & State of the
Art”, Proc. Int. Symp. Turbo Codes & Related Topics, Brest, France, Sep.
1997, pp. 1-11.
[26] A. P. Worthen and W. E. Stark, “Unified Design of Iterative Receivers using
Factor Graphs”, IEEE Trans. Inf. Theory, vol. 47, no. 2, pp. 843–849, Feb.
2001.
[27] F. R. Kschischang and B. J. Frey, “Iterative Decoding of Compound Codes
by Probability Propagation in Graphical Models”, IEEE J. Sel. Areas Commun., vol. 16, no. 2, pp. 219-230, Feb. 1998.
[28] J. Dauwels, S. Korl and H.-A. Loeliger, “Particle Methods as Message
Passing”, Proc. 2006 IEEE Int. Symp. on Inf. Theory, pp. 2052-2056, 9-14
July 2006.
[29] T. P. Minka, “Expectation Propagation for Approximate Bayesian Inference”, Proc. 17th Annual Conf. Uncertainty in Artif. Intell., pp. 362-369,
Seattle, WA (USA), Aug. 2001.
[30] O. Zoeter and T. Heskes, “Deterministic Approximate Inference Techniques
for Conditionally Gaussian State Space Models”, Statistics and Computing,
vol. 16, no. 3, pp. 279-292, Sep. 2006.
[31] V. Smı́dl and A. Quinn, The Variational Bayes Method in Signal
Processing, Berlin, Germany, Springer, 2005.
[32] C. M. Bishop, Pattern Recognition and Machine Learning, Springer,
2006.
[33] A. R. Runnalls, “Kullback-Leibler Approach to Gaussian Mixture Reduction”, IEEE Trans. on Aer. and Elec. Syst., vol. 43, no. 3, pp. 989-999,
July 2007.
[34] T. Li, M. Bolic, P. Djuric, “Resampling Methods for Particle Filtering:
Classification, Implementation, and Strategies”, IEEE Sig. Proc. Mag., vol.
32, no. 3, pp.70-86, May 2015.

G. M. Vitetta et al., Marginalized Particle Filtering . . .

44

[35] G. M. Vitetta, D. P. Taylor, G. Colavolpe, F. Pancaldi and P. A. Martin,
Wireless Communications: Algorithmic Techniques, John Wiley &
Sons, 2013.
[36] J. Hagenauer, E. Offer and L. Papke, “Iterative decoding of binary block
and convolutional codes”, IEEE Trans. Inf. Theory, vol. 42, no. 2, pp.
429-445, Mar 1996.
[37] T. Schön, “Example Used in Exemplifying the Marginalized
(Rao-Blackwellized) Particle Filter”,
Nov. 2010 (available at
http://users.isy.liu.se/en/rt/schon/Code/RBPF/Document/MPFexample.pdf).

