Iterative State Estimation in Non-linear Dynamical Systems Using Approximate Expectation Propagation
Authors: Sanket Kamthe
So Takao
Shakir Mohamed
Marc Peter Deisenroth
Venue: Accepted by TMLR ()

Abstract: Bayesian inference in non-linear dynamical systems seeks to find good posterior approximations of a latent state given a sequence of observations. Gaussian filters and smoothers, including the (extended/unscented) Kalman filter/smoother, which are commonly used in engineering applications, yield Gaussian posteriors on the latent state. While they are computationally efficient, they are often criticised for their crude approximation of the posterior state distribution. In this paper, we address this criticism by proposing a message passing scheme for iterative state estimation in non-linear dynamical systems, which yields more informative (Gaussian) posteriors on the latent states.  Our message passing scheme is based on expectation propagation (EP). We prove that classical Rauch--Tung--Striebel (RTS) smoothers, such as the extended Kalman smoother (EKS) or the unscented Kalman smoother (UKS), are special cases of our message passing scheme. Running the message passing scheme more than once can lead to significant improvements of the classical RTS smoothers, so that more informative state estimates can be obtained. We address potential convergence issues of EP by generalising our state estimation framework to damped updates and the consideration of general $\alpha$-divergences.

Published in Transactions on Machine Learning Research (06/2022)
Iterative State Estimation in Non-linear Dynamical Systems Using Approximate Expectation Propagation
Sanket. Kamthe"
sanketkomthe@@gaallcoo
Deportment of Computing, Imperial College London Chief Technology Office, Morgan Chase
So Takao"
80. so.takao &wc.l.ac.cc.c..m uk
UCL Centre for Artificial Intelligence, University College London
Shakir Mohamed
Department of Computer Science, University College London
Mare Peter Deiscnroth
UCL Centre for Artificial Intelligence, University College London
Reviewed on OpenReview: https: // openreview. por/forun? id- myyt4ffdo.
Abstract
Bayesian inference in non-linear dynamical systems seeks find good posterior approxi- mations of latent state given sequence of observations. Gaussian filters and smoothers, including the poxtended/unscentee engineering applications, yield Gaussian posteriors on the latent state. While they are computationally efficient, they are often criticised for their crude approximation of the pos- terior state distribution. this paper, we address this criticism by proposing a message passing scheme for iterative state estimation in non-linear dynamical systems, which yields more informative (Gaussian) posteriors on the latent states. Our message passing scheme based on expectation propagation (EP). We prove that classical Mauch-Tuune Striebel (RTS) smoothers, such as the extended Kalman smoother (EKS) the unscented Kalman smoother (UKS), are special cases of our message passing scheme. Running the message passing scheme more than once can lead to significant improvements of the classical RIS smoothers, so that more informative state estimates can be obtained. We address potential convergence issues of EP by generalising our state estimation framework to damped updates and the consideration of general a-divergences.
Kalman filter/mmoohher which commonly used in
1 Introduction
Dynamical systems (state space models) powerful tools for describing (latent) time-series with noisy observations. Dynamical systems are commonly used in signal processing, control theory, and climate science. Given sequence of noisy observations, estimating the (posterior) distribution of the latent states in a dynamical system plays an important role, e.g., in optimal control theory (Aström, 2006; Bertsekas, 2005), data assimilation (Lakshmivarahan & Stensrud, 2009), or robotics (Thrun et al., 2005; Barfoot, 2017).
For linear systems with Gaussian noise, an optimal estimate of posterior distributions of the latent states conditioned on the observations up to the current time step is given in closed form by the Kalman fil- ter (Kalman, 1960). The estimate conditioned all (including future) measurements of finite time series is given by the Rauch Tung Striebel (RTS) (Rauch et al., 1965) smoother. The efficiency of the Kalman
"Equal contribution. Also with DeepMind.

---

Published in Transactions on Machine Learning Research (06/2022)
filter and the RTS smoother stems from the closure property of Gaussian densities under linear/affin trans- formations, conditioning, and marginalisation. These closure properties allow us to derive recursive filtering and smoothing equations in closed form (Anderson & Moore, 2005).
In non-linear systems, these computationally convenient properties no longer hold; in particular, a non-linear transformation of a Gaussian density yields a non-Gaussian density. Typically, we tackle this challenge either by approximating a non-Gaussian density by a Gaussian density (Maybeck, 1979; Julier et al., 2000; Arasaratnam & Haykin, 2009) by using sampling-ber bbatn methods, c.g., the particle filter (Doucet et al., 2000), to represent non-Gaussian densities by means of particles. In this article, we focus on the former, i.e. finding better Gaussian approximations for non-Gaussian densities.
By making a Gaussian approximation to the distribution of the latent states, we are in fact either implicilly or explicitly linearising the non-linear transformation. The quality of these approximations depends on the mechanism used for the lincarisation. A crude approximation may introduce large errors, which we can alleviate by iteratively correcting previous estimmates (Lefebvre et al., 2004). For example, in the iterative extended Kalman filter (IEKF) (Corporation, 1974), measurement updates are carried out iteratively until convergence, where the point around which one linearises the measurement function is refined iteratively, This is shown to be equivalent to Gauß-Newton method for MAP estimation (Bell & Cathey, 1993). An iterative version of the extended Kalman smoother has also been proposed (Bell, 1994), which, again is equivalent to a GauB Newton method for MAP estimation. Recent progress (Raitoharju et al., 2018; Gultekin & Paisley, 2017) shows that the Bayesian framework and iterative updates can give significant improvement over standard filtering and smoothing.
All Gaussian filtering systems assume that we can project the true postcrior onto a Gaussian (Thrun et al., 2005) with a reasonable approximation error. More generally, this can be interpreted as approximate assumed density filter (ADF) with a Gaussian assumption on densities. In the ADF, the approximate posterior is obtained by minimising the Kulbback-Leibler (KI.) divergence between the true posterior and the assumed density. For a Gaussian assumed density, this amounts to computing the first two moments of the true posterior and approximating it with a Gaussian that possesses the same first two moments. Tn the Assumed Density Smoother, as well in the iterative Kalman smoothers, there no mechanism to revisit previously filtered estimates and correct them based on all the observations.
This shortcoming was addressed by expectation propagation (EP) proposed by Minka (2001a;b); Opper & Winther (2001). EP extends the ADF by adding a mechanism to iteratively update previous computations without double-counting observations. For linear Gaussian systems, one full iteration of EP is identical to the Kalman filter-smoother (Minka, 2001a; Heskes & Zoeter, 2002; Qi & Minka, 2003). Hence, EP is optimal for linear systems.
One can extend EP to non-linear systems in a manner that is similar to the generalised Gaussian filters. This extension is motivated by results in (Dehaene & Barthelmé, 2016), who show that EP has an order of magnitude lower asymptotic error bound compared to a canonical Gaussian approximation for non- linear systems. Furthermore, applying EP to non-linear systems allows us to perform multiple sweeps of RIS smoothing, iterating over measurement and transition non-linearities conditioned on all available observations. We note that the IEKF systems can only iterate over measurement non-linearities and the IEKS (Bell, 1994), being a MAP estimator, does not properly account for posterior uncortainties. EP can be implemented a message passing algorithm. Hence, it can be used as a distributed (and asynchronous) state estimator and is not limited to the sequential nature of iterative Kalman filters.
One challenge with implementing EP is that its efficacy relies on the case of computing moments with respect to the true posterior; for non-linear systems, the exact posterior is often intractable. In this work, we propose to approximate the true posterior by linearisation methods, such as the unscented transform (Julier et al., 2000). Hence, we only implement approrimale EP when we apply EP to non-linear systems. We prove that this approximate EP exactly recovers the respective RIS smoothers, e.g., with the unscented transform approximation we obtain the unscented Kalman smoother (Särkkä, 2008), and if we explicitly linearise functions, we obtain the extended Kalman smoother. This allows us to extend EP to all non- linear Gaussian fitter//mmooterrs, and generalise them to a setting where posterior state distributions can be iteratively refined to yield better Gaussian approximations.
2

---

Published in Transactions on Machine Learning Research (06/2022)
Contribution
1. We prove that any Gaussian filter for non-linear dynamical systems is equivalent to single iteration of appresimale EP. This equivalence will allow us to implement (approximate) EP using Kalman-fllter-syyle updates and without explicitly computing gradients of the log-partition function. Additionally, we can improve non-linear Gaussian smoothers by iteratively refining the marginal posterior distributions using approximate EP.
Because of approximate EP's equivalence with Gaussian smoothing, we can further improve estimation of the marginal posteriors by pre-existing techniques in the EP literature: (1) Damping, which reduces oscillatory behaviour; (2) Power EP, which allows us to interpolate between EP (mode averaging) and variational Bayes (mode seeking) for iterative smoothing in non-linear state-space models.
We use a classic non-linear system benchmark, the uniform nonlinear growth model (UNGM) (Doucel et al., 2000) with a fixed seed as running example to illustrate properties of the iterative smoothing scheme. The article is organised as follows. In Section 2, we establish the notation and review classic Kalman filtering and smoothing, which we then apply to the UNGM benchmark. Section 3, we introduce EP a way for iteratively computing posterior distributions on latent states. Furthermore, of derivative-baset (approximate) EP is equivalent to smoothing with non-lincar Kalman filters, which does not require any derivatives. We demonstrate that EP further improves the standard smoothing results of the UNGM example. Section 4.2 then addresses potential limit-cycle behaviours of EP by using damped updates (Heskes & Zoeter, 2003), and then further improves the results by using a-divergence EP (power EP) (Minka, 2004). In Section 5, we analyse the efficacy of our method and its sensitivity with respect to a, the damping parameter, and the linearisation technique.
prove that the first iteration
2 Gaussian Filtering and Smoothing
Consider the state-space models of the form
x; = f(x(-1) + wis
(1)
y. h(x;)) + Vt.
(2)
where 2e E R is an unobserve//latent state and y, E R is measurement at time step t 1, ...... T. Moreover, f and h are the (non-linear) transition and measurement functions, respectively, and the noise processes w, and v, are ii.d. zero-mean Gaussian with diagonal covariances Q and R, respectively. The initial state distribution is given by p(zo) = N(za Mo E..). The purpose of this section is to provide background details on Gaussian filtering and smoothing algorithmss on the stale-space model (1) (2).
2.1 Optimal Filtering and Smoothing
Bayesian filtering corresponds to finding a posterior distribution p(ær(v1:t) on the unobserved latent state æt, given observations F1r : (y1..... y2) up until the most recent time step It proceeds by sequentially computing time update and measurement update (Anderson & Moore, 2005), as explain below.
Given filtering distribution 1/:::)111::-); the time update
p p(2:///a-1) = p(: dx,-1
(3)
computes predictive distribution of the state æ; given observations up to time step t - 1. Here, p( N((x, (x= is the transition probability, which follows from the dynamical system (1). The measurement update then incorporates the observation ye at time t by computing
p( xz 2
"Cdde p(
(4)
day,
3

---

Published in Transactions on Machine Learning Research (06/2022)
which yields the filtering distribution. Here, p(y, (x)) = (3, hh(z)), R) is the observation likelihood, which follows from (2).
On the other hand, in Bayesian smoothing, the posterior distribution all observations in the time series is computed. Here, we consider the Rauch- -ann--ttree---b Striebel (RTS) smoother, which computes this posterior iteratively by applying the smoothing update
of the latent state æ, given
p(Z - / plas/eee1 1 V1r))((
= / plzes y Via) (aa (l (1r))d
(5)
where we used the Markov property p V1a) p (1.r))
Equations (3) (5) are generally intractable for non-linear systems. However, if we focus on Gaussian approx- imations to the filtering and smoothing distributions, we can derive equations for the time update, measure- ment update and the smoothing update, provided that we can approximate the densities p(æ;-1 3. (1 p( y: 11)=-1) and p(a,, 1++1 3/:e)) respectively by joint Gaussian distributions (Deisenroth, 2010; Särkkä et al., 2015).
To this, let us take for example approximation of p(æt; Ye Vil-1:t1) a joint Gaussian
(
Σ
Eay E
pl y 12l =N
2c y:
(6)
Σ
2c f/c-1
The subscripts t|t 1 indicate that the statistics of random variables (e.g., the mean) at time step t are conditioned on information/ poaaaaaiiooobsevatiinn up to time step / - 1 and the superseripts indicate whether the quantity is observed latent or observed space (for example Σay -j == cov(æ, Vel Vt-1))-
By exploiting the joint Gaussianity in (6), the measurement update equation can be computed by Gaussian conditioning, leading us to the approximate filter distribution M(x://r. x (N 2(jp), where
Míje = pij-11 + Ke 1%)p1)).
(7)
E = 11p111 K 8
(8)
:/2 Here, the matrix K, = 2y E - - corresponds to the Kalman gain in linear systems. Now, to arrive
at the Gaussian approximation (6), we do moment matching, which requires one to compute the means and covariances of the true distribution p(z,, v.lu...). This has the following exact expressions
p /
(9)
- f(x==i)((((==)) p( - 1 Moje-1 (4 Q,
(10)
1 1|1-1 ((a.))((
(11)
f|t-1 ((lz:)(kaax)))ppprrrrl.l.. dx; php-- (10) + R,
(12)
22 - (h(x))))) p(x, dæ, u .- (1t
(13)
Generally, these cannot be computed in closed form when f and h are non-linear. Hence, we resort to methods that approximate these integrals, such as linearising the non-linear terms or applying the unscented transform (Deisenroth & Ohlsson, 2011; Särkkä, 2013).
Similar update equations be obtained for the RIS smoother by approximating the joint distribution p(#; 1 lV by joint Gaussian (Deisenroth et al., 2012), which leads to closed-form equations for the mean and covariance of the smoothing distribution, given by
P pdj + Le - 11+1)))
(14)
Eip E/+1444 Lı [Et+TT Ext11)) LT,
(15)
4

---

Published in Transactions on Machine Learning Research (06/2022)
respectively. Here, the matrix L, 2i.+1je (E(+1j)) is called the smoothing gain, where E 21,+1jp denotes the cross-covariance between 2y and 11+1, which has the true expression
E/o+ = / x (f(x)) (( dæ, H
(16)
Again, for non-linear transition function f, the integral is intractable. Hence, we approximate it by means of, for example, linearisation or the unscented transformation.
For the time update, if we similarly approximate the distribution p( 1, 2: by a joint Gaussian, we can obtain an approximation to the predictive distribution (3) by a simple Gaussian marginalisation.
2.2 Implicit Linearisation in Non-linear Gaussian Filters and Smoothers
Applying moment matching approximate the true distributions p(z 21b3a11)) p(a,, y ValMia-1 p(a;, Fe+1 1y).) as joint Gaussians in Section 2.1 implicitly lincarises the respective transition and mea- surement functions, two random variables are jointly Gaussian if and only if they are related by affine transformation up lo zero-mean Gaussian noise. This is stated more precisely in the following Lemma. Lemma 1 (Implicit lincarisation). random variable (a, y) is jointly Gaussian if and only if there exist matrices M E REXD, P ER and a vector v E R such that y = Mæ e, where e ~ ~N(0, P).
Proof. If (xr y) is jointly Gaussian, we can write
(
~N .
1 11
))
(17)
where E0e (E"M).. On the other hand, y Mæ + u te with ~ Noo, P) were hold, then for x N(px, )*,, we have
x ( N ((M M
uh + ME"
Σ
E'M M MET M + PP)
(18)
We see that the RHS of expressions (17) and (18) agree and only if
M Σ" (""))--
(19)
p = Σ" - M
(20)
v = nDx -
(21)
It is therefore always possible to write y = Mx + v + & for random variables x and that are jointly Gaussian and vice versa. Note that the matrix P in (20) is guaranteed to be positive semidefinite as is exactly the expression for the covariance of the conditional distribution p(yla:).
We illustrate this idea of linearisation that arises as result of forming joint Gaussians (and vice versa) in the following two examples. In the first, look the case when linearise the nonlinear transition function and how that results in joint Gaussian approximation between two random variables. In the second example, look at using particle-based methods to approximate a general joint distribution p(æ,3 y) as joint Gaussians and how that implicitly leads to a linear relationship between æ and y.
the Taylor approximation to explicitly
Example 1 (Taylor approximation). Let ~ N( Σ²) and y|æ ~Ng((x) P), where g is non-linear differentiable function. Since the joint distribution p(æ, y) N(x ",,"", N(ylg()))) ( P) is non-Gaussian, we approximate it by a Gaussian by means of explicitly lincarising around μ².
g(x) R g(x") + M(x - (+))
(22)
where M = V.( This gives us the linear relation y = Matute, NO P), where v = g(u*)~ and the following Gaussian approximation to p(x, y) holds:
-M
approx N
((11E)))
E
MY M + +P))
EN T'
MM
(23)
y
5

---

Published in Transactions on Machine Learning Research (06/2022)
8
EP
EP
UKF
8
UKF
7
UKS
UKS
RMSE
NLL
6.
2
10 20 30
40 50
10
20
30
40
50
EP Iterations
EP Iterations
Figure 1: Root mean square error (RMSE) and negative log-likelihood (NLL) performance of estimation algorithms applied to the UNGM example (single seed). The unscented Kalman filter and the unscented RTS smoother form the basclines (dashed lines). EP starts from the RTS smoothing bascline, as its first iteration is identical to the RTS smoother. We can see that the EP estimate improves on the smoothing baselines with an increasing number of EP iterations. The oscillatory behaviour of EP can be addressed by damping, which we discuss in Section 4.1.
Example 2 (Monte Carlo and unscented transformation). Alternatively, we approximate the joint dis- tribution p(æ, y) in the previous example via. a partiele-based representation, that is, represent p(æ) using collection of particles (which can be determined deterministically as in the unscented transformation or stochastically as the particle filter) and propagate them forward by to compute the moments μr", Σ", Eey empirically. We can then use the computed moments to deduce the implicit linearisation y = Mx+u+e with e~ N0, P), where M, , P, v are computed as (19)-(21) The Monte Carlo and unscented Iransformation considered in our experiments in Section 5 follow this framework.
Commonly used Bayesian filters and smoothers for non-lincar systems are the extended Kalman filter (EKF), the unscented Kalman filter (UKF) (Julier et al., 2000), and the cubature Kalman filter (CKF) (Arasarainam & Haykin, 2009), all of which compute means and covariances of the unobserved state dt- Since all update equations of these non-linear filters and smoothers these methods either explicitly (e.g., EKF) or implicitly (c.g., UKF, CKF) linearise the transition and measurement operators in (1)-(2,, which incurs approximation errors.
be expressed (7), (8), and (14), (15) respectively,
To account for the approximation errors, we may improve the fltering/smoothiin results by iteratively refining the state estimate. For example, in the iterative extended Kalman filter (IEKF) (Senne, 1972), an iterative measurement update used to improve the standard EKF solution. On the other hand, similar results for smoothing are not well explored. However, obtained in filtering by ileratively linearising the transition function well. This is challenging, as, unlike the measurement function that only depends on the measurement at time step t, the smoothing result depends on all observations from time steps to T.
would expect to similar improvements
The approach that we propose here to apply Expectation Propagation (EP), a general approximate inference scheme that iteratively refine approximations of the target distribution (Minka, 2001a). Standard EP only applicable for linear functions, but it exactly recovers the Kalman and RTS smoother for linear systems (Minka, 2001a). Similar to based iterative filters, expect the EP updates to iteratively improve smoothing results in the non-linear case.
2.3 RTS Smoothing Applied Non-linear Benchmark
Before describing EP in detail (sce Section 3), we motivate the usefulness of its iterative refinement of the posterior state distribution instance of the well-known non-linear benchmark:

---

Published in Transactions on Machine Learning Research (06/2022)
Jc (x:)
p(2)+1122)
(a(()())
fa(())
(a
zt
att
Sa (x++1)
fa(œ;)
fa (2((1)))
(a) Factor graph
(b) Pully-factored factor graph. Factors between nodes are split.
Figure Factor-graph representations of a dynamical system. (a) Exact factor graph of a dynamical system. Factors (black squares) correspond conditional probability distributions, e.g., p(æ+41 (æxx) for the factor between D2 and 2++1- Circles are random variables. Observations Yt not displayed deterministic. (b) The fully-factored factor graph illustrates how local messages are being sent to update the marginal posterior at each time step. The splitted factors are obtained by applying belief propagation on (a), that is (xx))) = f p(æ++1
they are
[( [ (#t++) 1) (a ));,-1 and fo (2++1) = p(ae4 [2,)/o (a..) (axs (xx, d.x.
21 11-1
25m
8 ccos(1.2(t 1)) + wr:.
(24)
T4
V/ 20
Ub:
(25)
where we use the same set-up as in (Doucet et al., 2000): we (N,, 2 ~ (V, 10) and p(ro) M(xo00,5), 5). We consider trajectories of 100 time steps. This benchmark is challenging for state estimation as the transition density around x x 0 is bimodal and symmetric. The scaled quadratic measurement function makes these modes indistinguishable.
We compare the results of EP with the unscented Kalman filter and smoother (van der Merwe & Wan (2004)). Results shown in Figure 1, where plotted the root mean square error (RMSE) and the negative log- likelihood (NLL) of its predictions, all measured in latent space (see Appendix A.3 for the precise expressions of these metrics). While the RMSE only mcasures the quality of the posterior mean, the NLL is a measure that allows us to assess the quality of the approximate posterior mean and covariance. Observe that the results from RTS smoothing and the first iteration EP are identical, a result we prove in Section 3.4. We see that EP iterative smoothing generally improves the results obtained by single pass. In the next section, we develop the expectation propagation algorithm used generate the solid red curve in Figure 1.
3 Expectation Propagation
After this brief motivation, we are now ready to describe expectation propagation (FP) more detail.
In the context of dynamical systems, EP is typically dcrived using factor-graph formulation (Qi, 2005; Deisenroth & Mohamed, 2012), where the (marginal) distribution over the latent state p(æ, ( i.e., the smoothing distribution, is represented as product of factors/ merss//ssseess fi. so that p(æx|U1..) = II. fi see Figure 2 for an illustration. There are three types of factors ft in this context: forward, backward, and measurement messages, which we will denote by the symbols D., <1, 4, respectively. The forward factor fo is a mhesage/informaat
that is carried from all the latent states and observations prior to the current latent state zt. The backwards factor fa sends information from future states to Ht; and the measurement factor Sa incorporates information from the current observation. Thus, each (marginal) state posterior smoothing distribution can be written product of factors/ factors/messager
p(Z:lurr) (xx)
fa( fa(xx)
(26)
= (forward) (measurement) (backward).
2 A thorough comparison across different filters and seeds discussed Section 5.

---

Published in Transactions on Machine Learning Research (06/2022)
For the non-linear systems we consider here, the exact factors f((z)) above and therefore the posterior distributions p((a//1r)) are intractable.
EP (Minka, 2001b;a) is a deterministic Bayesian inference algorithm that can be used to approximate the complex and an often intractable smoothing distribution p(=z/3).r) = II. f((x,) by a distribution of the form ((x;)) = [I; (2rz), which a product of simpler and tractable factors qi- The underlying assumption is that if we can find good approximating factors qi, we may get better estimates of the intractable posterior distribution p(æ, (.r)). In practice, restrict the factors 9; to families of certain distributions, such that they are closed under multiplication and division (Seeger, 2005), e.g., Gaussians. We get approximatations in this family by a procedure known as projection.
Definition 1. Let F be a family of probability distributions. The KL projection of p onto F is a distribution in F that minimises KL (pllg) over all distributions g E F. We use the shorthand notation projz [-] to define the operator that projects arbitrary distributions onto the family F:
projz [p] == argmin KL (pllq).
(27)
gEF
In EP, we assume that the approximations ((x)) and the factors q((z)) come from the exponential family, This has the useful property that. the operation projæ [p] simply amounts to finding q E F whose expected sufficient statistics correspond to those of p (Seeger, 2005). For example, if F is the family of normal distributions, then proj= [p] = [M, Σ) where je and E are the mean and covariance of p.
3.1 EP Algorithm
The iterative refinement of each factor qi is carried out as follows in EP (Minka, 2001a). First, we compute the cavity distribution
(x,)) x ((
(28)
which contains all the information about the distribution of m, except for the 7th factor, Lee. the "rest of the distribution" (Minka, 2004).
Second, we form the (unnormalised) tilted distribution f((x))g,((;)) by combining the "true" factor f(())) with the cavity distribution (,,(x)) and determine an updated factor (:(m)), so that q (x)))) approxi- mates fa(2;) (Z,)) as close as possible in the KL sense. We achieve this via projection (Definition 1). This gives us the updated posterior g( and approximate factor gi (,,,)) which reads
q(x,) = = [f((x)))
(29)
(x.,,) new x g(z;)
(30)
(xx))
respectively. Note that the cavity distribution q, provides the "context" within which to find the updated approximate factor qs: the Gaussian approximation 9, = f; need only hold locally, where the posterior defined by the remaining factors d has the most concentration of mass.
For linear systems with Gaussian noise, the true factors fax are also Gaussian, hence the projection step and the implied Gaussian approximation are exact. This means that in this case, the factors 9 are unique and optimal, i.e., 9i = fi (Minka, 2001a). The Kalman smoother obtains these results in a single forward-backwart sweep through the factor graph. For non-lincar systems, however, computing the exact factors fi are intractable. Hence, we use approximations to compute these factors (Qi & Minka, 2003), which requires more than a single iteration of EP to converge. If we iterate multiple times, all factors contain information from observations at all times. When correctly set up, we can at cach itcration improve the quality of the factors and obtain closer approximation to the true smoothing distribution (Vrr).
3.2 EP Updates for Generalised Gaussian Filtering/Smoothing
In the following, we provide explicit computational detail of the EP algorithm when the approximating family F is the family of Gaussians. For simplicity, we drop the subscript F from the notation projz
8

---

Published in Transactions on Machine Learning Research (06/2022)
Algorithm Gaussian EP for Dynamical Systems
Init: Set all factors q to (0 oof); Set g(æ1) = p(æ1) and marginals g(æx>1) 1000 repeat
for t= 1 to T' do
4:
for all factors (xx)), where = D. A. do Compute cavity distribution:
5:
q((())) x a(<;)/q
(31)
6:
Update posterior by projection:
now (x:) proj
(32)
Update factor:
(33)
end for end for
10: until Convergence or maximum number of iterations exceeded
for the projection in this case. The pseudo-code in Alg. 1 describes the main steps of Gaussian EP. This iteratively updates the marginal ((xx)) and the messages qi(x)), E {D, < that form this marginal.
Compute cavity distribution: First, the cavily distribution (28) ( is computed (step 5 in Alg, 1) by removing factor (xx)) from the marginal g(xx). Since g(x)) and qx(x;) are Gaussian, the cavity distribution is also (scalod) Gaussian,
(x,2,) x (<<:)/q((x)) N(x (A) )
(34)
The terms M. Σ) can be computed via the Gaussian division rule (132)-(134) (Appendix A.2),
Σ E (E)) (E))-1)-
(35)
14 / Σ ( p4 (E))
(36)
Update posterior by projection: Next, the Gaussian projection of the tilted distribution (x,)4, is computed (see Def. 1). This requires one to compute the moments of f((x))9;((Z;) which can be deduced from the derivatives of the log-partition function (Minka, 2001a;b; 2008)
log Zi( 1 x = log / ((2))))
(37)
The moments of the updated (marginal) posterior ( can then be computed as
14 new uh +
(38)
new Σ
mI 2V DEWE
(39)
where W/ == d log Z 7
(40)
and Mv = log Z: /z)
(41)
Update factor: Finally, we obtain the updated factor o. (see (30) again via the Gaussian division rule, since both the updated marginal and the cavity distribution (unnormalised) Gaussians.

---

Published in Transactions on Machine Learning Research (06/2022)
In some cases, it is casier to compute the updated factor q.n" first; the updated posterior is then obtained by multiplying the updated factor with the cavity distribution: quew gfth Tr--
Remark 1.
Remark 2. While the focus of the paper is on smoothing, i.e., the offline inference setting, message passing in the online setting via a moving window (fixed-lag smoothing) could also be implemented. Here, we would use the messages of the last & time steps to find an approximate filtering distribution p(ar:ltt-erd kr). Iterations within the considered window would still give improved filtering distributions.
3.3 Implicit Lincarisation Requires Explicit Consideration
We now highlight an important detail in the computation of the derivatives V MAA V xyy, which are required in the EP update formulas (39). For nonlinear systems, these cannot be computed in closed-form, hence we use approximate derivatives practice, which obtained by replacing the partition function Z; by its approximation, denoted Zi, in the computation of (40) (41). As we shall see more concretely in the next section, this approximation in general takes Gaussian form, making it possible to compute derivatives.
As an illustrative example, consider the partition function Z Sa (m:)gq (xx dx., corresponding to the measurement message, where Sa (xx)) = p(ya = Nwy, (((+). R). Approximating the tilted distribution p(y, x ( by a joint Gaussian in æ, and y/., we get the approximation
Z = / Nv
Øt y,
My u u2/r
Eera
)) dx, = NV..I lu E21a),
(42)
Mg1a
which takes the form of Gaussian. Here, Me uiie D//4 M correspond to approximations of the integrals Ey [x(=))) covgta [x()) and covgia [x,, hxx))l respectively, which can be obtained using the same tech- niques as we saw in the examples in Section 2.2, by applying the unscented transform, linearisation etc. The \a terms pig r/ and E are equivalent to p4 and L) the superscripts were simply added to distinguish them from moments in the observed space.
Note that Zo in (42) is Gaussian in the observed space ye while the derivatives (40)-(41) are taken with respect to moments in the latent space at. Using the chain rule, the derivatives thus read
dlog Z,
Ölog Za Oppia
Olog Za ax
(43)
M/a
du. 2/a
Outla
Oppia
aΣyL 0p, r/a
E\A = dlog Z4
Oog ZA 0pt 1/4
log 2, VA
(44)
d 0pg
0Σ,1
which requires knowledge of the derivatives of moments the observed space with respect to moments in the latent space (terms highlighted in red). Crucially, these do not necessarily vanish when æ, and y, have non-zero correlation we expect (n v/ 2)) to change with (u, a/a 1 if they are correlated.
Computing the terms in red can be achieved if we have an explicit relation between the random variables æ, and yr- Since It and y, are jointly Gaussian in the approximation we made to our partition function (42), then by Lemma 1, they are related by an alfine transformation up to zero-mean Gaussian noise, that is, y: Jix, + un + 7, for some matrix J, EXD, vector u, E R. and n, (O R). Thus, have
uihe = J,H." 1/ + u.
(45)
Σ/1A J, JJ. RR.
(46)
We now make an important assumption that the matrix/vectors J., u, do not depend on the mo- ments μ, ais E a) which is equivalent to saying that the approximate conditional distribution ((), = N(y, Jı + u. R) is independent of the moments pl a/A EF/ This is a reasonable assumption since the likelihood p(y,aar (x)) of our original, unapproximated system is independent of the prior p(æx).
10

---

Published in Transactions on Machine Learning Research (06/2022)
With this assumption, using (45)- we get
0upla J, EXD
Oug
v/A EXDXD
O
a "k [Je]as [Tele-
-EEH ODEREXAN EXEXD
REEEDD
Op
0 | Ik
(47)
Thus, from (43)-(44), we obtain the following closod-form formulas for the derivatives in question
Mia (y/
pi ( [ 1 ERRR. IxD
(48)
V D/A 2 a/a MA J, ( E ERDx,
(49)
where we used the identities (126) in Appendix A.1 to compute the derivatives of log- Gaussian
3.4 Derivative-based EP Updates and Kalman Smoothers
Section 2.1 provides generic equations for Bayesian filtering (7) and RTS smoothing (14) (15) in non- lincar systems with Gaussian state distributions; the only requirement to be able compute means and covariances of Lwo joint distributions to apply these updale equations. In Section 3.2, based on derivatives of the log- oog-partiionn function, derived generic update equations (38)-(39) for computing means and covari- ances of posterior distributions of the latent states of a non-linear dynamical system, but now formulated as a pressage-passigg. algorithm, which can potentially be iterated and may lead to more refined posterior distributions than RTS smoothing.
In the following, we will show that these generic update equations are identical for a single iteration of EP when updates are performed in the following specific order: As in Gaussian filtering, we start with the time update "forward factor" (D>) update, followed by a measurement update "up factor" (4) update and repeat until we reach t T. Then move backwards until reach t = 0 by updating the "back factor" (4) similar to the backward pass of RTS smoother. The implication of this statement is that we can interpret the EP computations in the context of classical Gaussian filtering and RTS smoothing for non-lincar dynamical systems. These equivalences include all classical filters and RTS smoothors for lincar and non-linear systems, including the extended/ wttpe////ssneeeee/eoeee cuubature Kalman filters/mmoohhers.
Forward message To compute the mean and covariance for the forward message fo (æt), we obtain it directly instead of using the derivalive-based update formulas (38) (39), as this turns out to be much simpler. Note that the approximate forward message has the expression
(x,) / f p(z (a D /
M(x:If(E(11) (
(50)
which we from the factor graphs in Figure 2. Due to the non-linearity of f, the above integral is intractable. Hence, approximate the integrand by joint Gaussian, denoted (( using techniques such as lincarisation and unscented transform (sce examples in Section 2.2 for details). By Lemma 1, we then obtain the Gaussian approximation
L (xx) = f (;) = -N(x,) M 14 + U1-1, M M Pr-I),
(51)
to the true factor fo (xx) for some rule in the last equality. Now, since Pi-1 is the covariance of the conditional distribution p(z, (see proof of Lemma 1), we have P = Q. This results in the update formulas
Pr-1, where we used (18) and the Gaussian marginalisation
, = M, + 11-1.
(52)
Σ = M.
+
(33)
for the forward message. In accordance with Remark 1, we can then also update the approximate marginal posterior g(æ:) by a simple Gaussian multiplication, i.e., g(ar)) xx qo (m.))(pp))))
11

---

Published in Transactions on Machine Learning Research (06/2022)
Note that the approach that we use here to obtain the forward update formula (i.e., form joint Gaussian approximation to the integrand of (50) and marginalise) is identical to how we obtain a general time update in standard Bayesian filtering (see Section 2.1). The only difference is that in Bayesian filtering, the term corresponding to the cavity distribution (1-(( in (50) is given by the filtering distribution p(z)-11B/ra-1) (see (3) in Section 2.1). However, if we refer to Table for the correspondence between EP messages and filtering distributions, we see that the very first iteration of EP, the cavity distribution (x4(41))) is exactly equal to ((
Hence, we have the following result.
Theorem 1. In the very first iteration of EP, computing the approrimate forward messuge qp (xs) via the moments of Jo (axx) in (50) corresponds to the lime updale in classical Bayesian fillering algorillen, where the moments are given by (9)-
Measurement message For the measurement message, we employ the derivative-based update formulas (38)) We have already explained how to compute the (approximate) derivatives VAya and VIna for this case in Section 3.3, where the resulting expressions are given by (48) (49). Using this, we can then obtain the updated posterior q via (38)-(39) and updated factor qa by Gaussian division; see (30).
Theorem 2. The moments
H4 = 14 x/A
+
E
uia
(54)
Σ, = E Er A
V M4
- 2V
Ex
(55)
uia
of the update posterior ((x)), which are computed via the derivatives of log Za, are identical to the mo-
ments (7) of the filtering distribution (measurement update)
Mdjt 17/1-11 + Ky 14 1/t-11
(56)
Eie Eg 21/t-1 K, 1|t-1
(57)
obtained by any Gaussian filter, where
K²:= : E (2)) -1
(58)
is the "Kalman gain".
Proof. Since we are applying the EP forward and measurement updates in sequence and recalling the fact that we initialised the backward messages (a(x)) with an uninformative prior, we see that
n (a:)) = 40 (ax)) = N(: M D W 0, ol) = -(<x pij--
(59)
where we used the relation in Table to identify the forward message q= (xx) with the predictive distribution p(Z:/V1a-1) R Ai - 1 This gives us ii r/4 = - and E Eij-1 which allows us to
re-express the standard Kalman update (57) in terms of EP messages,
M h a + Ki )
(60)
E/ D.AA K.
(61)
where K, K;= =.=EAvA (M) We start by showing that the expressions for the means in (54) and (60) are identical, which is true if
m Kily, ))
(62)
holds. To see this, we use (19) and (48) to re-express the left- hand side of (62) as
Σ Hia
(48) E Ju | b) (1 - 4 Me (19) x (M ( - Smui
(g), 11 (
(63)
م (و))))))))) E (EE))) ( 14 (/ Eooia (E (y) u. ()))),
(64)
Ki
12

---

Published in Transactions on Machine Learning Research (06/2022)
Filtering/smoothingg distributions
Corresponding EP messages
Predictive distribution (time update) p(a 111--11
t- (a:)
Filtering distribution (measurement update) p(æilv):))
da (2:) (
Smoothing distribution p(z:/111.))
da (2))) = g((x))
Table 1: Glossary of
distributions and corresponding EP messages during the first itera- tion of EP. These equivalences can be verified inductively by comparing how they are both computed. Note that all posteriors are only approximate (Gaussian) in the non- non-linear case.
which verifies identity (62), thus concluding the proof for the mean.
To show the corresponding identity for the covariances in (55) and (61), we need to show that
VV E K, K,D,etA
(65)
MiA
holds.
Using (19) and (49), the hant-hand&e side (65) becomes
4/ م/A - 22 Mnrs
م (49)
1 J! (
(66)
(19)
E (2) 1 - M (D))) - pqpi ( (2)) -
(67)
(68)
=
concluding the proof of Theorem 2.
Theorems 1-2 state that computing the Gaussian approximate posterior state distributions with derivative- based updates of the forward and measurement messages is equivalent to the time and measurement updates in standard Gaussian filtering for non-linear systems.
Backward message
For the backward message, use the partition function
9 = (() x / [a(ax))N(aa)) 1 2)
Ze (14
(69)
where
fa(ax.)= / p(a0+1 (x)) (a (da+1)
(70)
to update the backward message (x(x;) and the posterior distribution p(zx(u117) using the derivativo-based update formulas (38)- (39). The true factor fa(#,) in (70) accounts for the coupling between of and # Substituting the expression for fa)) in (70) into (69) and reordering the integration yields
Za / (4411)
dx
(71)
(m
The inner integral in (71) corresponds to the forward message (50) at time step t+ which can be approxi- mated by dp (#++1) x (( /ff11
DP 1);1), as discussed earlier. Then, using the identity (131) in Appendix A.2, the integral in (71) can be computed analytically, giving us the following Gaussian approximation
Ż - | Di + 2)H.)
(72)
13

---

Published in Transactions on Machine Learning Research (06/2022)
We can now similar trick as before to compute the derivatives N V a/ in the backward message updates. Recall that when obtaining the approximate forward message t> (T++1) & (a we have approximated the distribution (()++1/=))) #++1 = a, + ve + E, for some matrix M,, vector vt, and & ~ (O Q) by Lemma 1. Then, assuming as before that M, we: Q are independent of (H, 14 Σ 9))3 deduce from (52) (53) that
by a joint Gaussian q( 2c+1), giving us the relation
0p M
0p
SE [+1 0,
0 E itlis [M.]a[Malj
0,
(73)
0p.
4
σΣ,3
dp.
0E,
kl
giving us
dlog Za
Ölog Za
0p
Olog Za 021+1
#/ =
dy, 14
0pt41 0p.-
0Σ 0p..
19
(p 10
(Ep + -(+1)-1 Mc.
'Me
(74)
log 7,
Olog Zg 0up+1
i log Za 0E
V
M, + 1 M,
P/
(75)
dΣ, 14
0pßt1 0
0Σ 0Σ
14
2
where again, used the identities (126) (127) in Appendix A.1 to compute the derivatives of log Za.
Theorem 3. The backward message updates via the derivatives of log Za
μ. = uu
(76)
E. = Σ
4/ A 24 E/ 1
2V
(77)
are identical to those obtained by the RTS smoothing update
Mūj = Mūje + LI - H(t1je);
(78)
EjjT Eije Li EEitIT - E1+111) L,,
(79)
where
L = :h:
DE 1,141 (2t+1))=
(80)
To prove this, we use the following auxiliary result.
Lemma 2. For any t, the following identity holds
(>) up)= 21'(m, up),,
(81)
where Me Eg are the moments of the full approximate posterior q(x,) N(xx "E,EE I
1u Ap
)
Proof. Using identities (129)-(130) in Appendix A.2, expand the terms the right-hand side of (81)
as
Σ = (Σ) 111 + (E Σ 1 = (1X)) + (E n²²
(82)
(83)
Subtracting (83) from (82) then gives us
(M, ( pP) (E) -+p)
Σ,
(84)
as required.
0
3This equivalent to saying that the approximate transition probability $(2.+1 =N(#0 M, Q) is independent of n4 14 s) which again reasonable assumption.
14

---

Published in Transactions on Machine Learning Research (06/2022)
Proof of Theorem 3. Using the relations in Table 1, we can re-express the right-hand side of the RTS smooth- ing equations (78) in terms of EP messages as
MIIT - μi \4 + Le(1++1 - p(+1);
(85)
Enp = E + L, [E+1 21+1) L
(86)
Furthermore, noting the implicit linearisation I++1 M, + Vs + E, and using (52), we have
/f [ - - 11+ pla ) ( mn dæ
11
(87)
Tab.1 /| (xt 1 (EL41 p T 2 (m:)4 (31+1) dxidæ
(88)
(52)
|(02= (2e - m1 (M =T g¹ (xx) + Bgii ( ui )ET
(89)
MM,
(90)
which gives
L = Eitt1 ( = M (E(+))
(91)
We now show the equivalence in the expressions for the means (76) and (85). To prove this, it suffices to show that
E A/a L.(M441 -Miz1)
(92)
By identity (136) in Appendix A.2, we have that
(Et + -1 ) = (E(+1) ((EL)) + (2(R))
(D)
(+1)
(93)
E1+
Inserting the expression (74) for V a\a and invoking Lemma 2, the left-hand side of (92) becomes
(1 - m(+1)
(94)
10
(11i: - 1l+1))
(95)
Lem.2 M (
E Exi (#e+1 - -1it1)
(96)
=1
=
2 M (1141 - mP+1) (91) Le(Het1 - -uP
(97)
This proves the equivalence in the expressions for the mean.
Next, we show the equivalence in the expressions for the covariances (77) and (86). To do this, it suffices show that
M/ 2V 2 Σ L. LE -E L.
(98)
First, use the matrix identity (137) in Appendix A.2, to get
(D + x1p,) = (E (E/+1) -1 (Et+1) -1 + (DD -
-1
(99)
=E++1
= - E(+1) (
(100)
15

---

Published in Transactions on Machine Learning Research (06/2022)
In (99), used the fact that the sum of the precision matrices of all messages gives the precision matrix of the marginal. Note that the precision of the cavity distribution is the sum of the precisions of all messages except the forward message. Combining (100) with (75), the left-hand side of (98) yields
- M, Eit + E M,E, 4 (100)
Σ M. (ER Ei+1 2141) D
M.D
(101)
=/
which concludes the proof.
Theorems 3 state that for a single iteration of EP, the Gaussian filtering and smoothing updates and the serivative-based EP updates with implicit lincarisation are identical. This means that the first iteration of EP can be implemented using derivative-free updates, i.e., the classical update equations for Gaussian filters and smoothers. For subsequent EP iterations, we can incorporate updated messages by using the following equations, which very much resemble the general Gaussian sttering//smooliie equations:
Measurement message
14, 14 14 + Ki 1)
(102)
E, Σ/ K2 K,DKaio
(103)
where for x, ~ (x,), we define*
Mi 1/4 E [k(x=)], Mg/4 cov,iaal(z20 + R.
(104)
D. = cov (
K, Ex E (x:/10)
(105)
Backward message
14, p4
14 Le(1441 14+1)
(106)
Σ, Σ + L. E1+1 Db
(107)
L, = cov a [Ke, f f(xx)) (11) -1
(108)
Here, we identify pa) = ( = u (z)+1) and = =
co + De....
Remark 3. As shown in Algorithm 1, we initialise all approximate factors (xx)) = No, ol), i (D, with zero mean and infinite variance (zero precision). This implies that the division and multiplication by these densities corresponds to subtracting/adlue adding zeros to precision matrices, respectively. Therefore, before updating 4 for the very first time, it holds that (,(x)) g(æ;), i.e., the cavity computation step leaves the posterior unchanged as division by vero-precision Gaussian (subtracting zero) has no effect. Therefore, in the first EP iteration, the generic derivative-fre frree update equations in (102)-(103) and (106)-(107) are identical to generic Gaussian filter/smoothr updates in (56)-(57) and (78)-(79), respectively.
To summarize, we showed that a single iteration of EP recovers classical Gaussian filters/ Elters/RTS smoothers special cases. We can use this close relationship to implement EP without the explicit use of derivatives, which is often more numerically stable for non-linear approximations (Seeger, 2005).
Remark (Computational Complexity). The derivative-free updates, which are nearly identical to the generic updates for Gaussian flters/smmoohers, tells us that the computational complexity of EP for state estima- tion is O(FF), where K the number of EP iterations, and F is the computational complexity of the corresponding Gaussian filter/smol mstoorrmmohher. This means that EP has the same computational complexity as a standard Gaussian fimoorrsm for each EP iteration.
16

---

Published in Transactions on Machine Learning Research (06/2022)
8
EP (a= 1.0, y=0.
EP (a 1.0, y-
EP with damping (a - 1-0,
EP with damping (a 1 y=
7
with damping (a 009, 0
Power damping (a= 0.9. y=
RMSF
6-
NLL
5-
4-
3
2
10 20 30
40
50
10
20
30
40
50
EP Iterations
EP Iterations
Figure 3: RMSE and NLL for EP (red), damped EP (green), and damped power EP (grey). We extend the rumning example from Figure 1, but omit the satering/ssmoothinn baseline graphs to focus the effects of damping and power. Both damped EP and damped power EP reduce the oscillations seen in EP substantially. Using power EP, we further improve the performance of EP and damped EP. The power of EP and damping factor are identical, which improves numerical stability as the power and damping cancel out in the computations of the updated marginals, see (116) in Algorithm 2.
4 Improving EP by Damping and a-Divergence
Revisiting Figure 1, we see that although EP improves the results of the RTS smoother, EP does not converge to fixed point but oscillates. In general, it is possible that in the case of non-linear systems, approximate EP may deteriorate and produce results inferior to ones of the RTS smoother. It may even diverge and fail to perform the estimation task. One can employ heuristics to determine termination conditions, but there are no established stopping or convergence criteria. In this section, we introduce damping and minimisation of the a-divergence as possible remodics to this well-known shortcoming of EP.
4.1 Damped EP
Oscillatory behaviour of EP, as observed in Figure 1, is well-known (Heskes & Zoeter, 2003). Standard EP has an iteration scheme without lermination criteria, and is not guaranteed to converge (Minka, 2001a). EP operates on one factor at a time. However, one can establish convergence criterion based on the covariance of all the factors in the algorithm. This is called double-loop EP (Heskes & Zoeter, 2003). In a double-loop EP, the inner loop is standard EP with damped updates for each factor. The outer loop is used to pick the optimal damping rate for the inner EP update based all the factors. However, this outer loop is prohibitively expensive for most practical applications. Moreover, non-linear systems, as discussed in this article, use approcimate EP, and the approximation would lose any convergence guarantees associated with the double-loop version.
Instead of computing an optimal damping for cach iteration, we can use fixed damping rate for the EP updates, which results in the so-called damped EP (Minka, 2004). This damping helps the algorithm to avoid limit cycles (Heskes & Zoeter, 2003). This changes the EP updates of the marginal and the factor qi in (29)-(30) to
of ((xx)) proj [f((m)))
(109)
q( new x (((x)))((xx))
(110)
q (xxx))ew new (((x))aew.
(111)
q (xI)
where o' corresponds to the undamped update in (29). Damping can be interpreted as a log-opinion pool between the current estimate and an updated estimale, where the damping factor E (0, 1] indicates how
4 As usual, the integrals below are only computed approximately, egg., by using linearisation unscented transform.
17

---

Published in Transactions on Machine Learning Research (06/2022)
much weight is given to the current estimate. For 7 1, we obtain the classical EP updates from (39)-(30). whereas for 7 = 0, we would not update al all.
The damping factor 7 is a hyper-parameter, which can be tuned depending on the filtering task. In Figure 3, we plot the effect that damping and a-divergence (this will be explained in the next section) has on EP performance, where we use the Unscented Transform to implicitly linearise the system (see also Figures 10 and 11 in Appendix C for the equivalent plots using the Monte-Carlo Transform and Taylor Transform to linearise the system, respectively). We see that damping factor of 0.9 reduces oscillations and stabilises EP in this particular example. We chose 0.9 to illustrate the effect of damping oscillations; in Section 5, we provide a more thorough sensitivity analysis for %
4.2 a- D-Divergence
In some cases, the damping factor 7 alone may not improve the results of EP. This is particularly true for non-linear systems, where the true posteriors are often multi-modal. EP, as discussed so far in this article, minimizes the KT. divergence KT. (pllg) between the true posterior and its approximation, so that the approximate posterior q has significant probability mass where p has significant probability mass, i.e., it averages between modes; this direction of the KL divergence is not modo-socking (Bishop, 2006). Depending the application, one may want lo encourage the algorithm to cover most of the probability mass (i.e., mode averaging behaviour), in which case standard EP may be desirable. On the other hand, it may be desirable in some applications that the algorithm picks one of the modes rather than averaging them, e.g., in the case of tracking car moving past a junction. To achieve this flexibility, we can modify the optimisation objective to encourage the desired behaviour by using power EP instead (Minka, 2004). Power EP minimises the a-divergence (Amari, 1985) instead of the KL-divergence used by standard EP (Minka, 2004). The a-divergence is a generalisation of the KL divergence (Zhu & Rohwer, 1995). We follow (Minka, 2005; Zhu & Rohwer, 1995) and define the a-divergence as
(xll)))) = ap(z) + (1 a))(z) p(z)" dx a(1 aa)
(112)
where (pll4) = 0 if and only if p=p = q and positive otherwise. For a = 1, we obtain the Kl.-divergence KL (pllg) (Zhu & Rohwer, 1995), for a = 0.5, the objective corresponds to the Hellinger distance and for a= = O, we obtain the KL-divergence KL (4llp) again, except with the arguments reversed. Clearly, the case a= recovers standard EP and the case = 0 recovers variational Bayesian inference. Thus, power EP can be seen a bridge between these two commonly used algorithms (Minka, 2005).
To increase the flexibility of our distributed state-estimation algorithm, propose to use the a-divergence as our cost function to which we find the approximate Gaussian posterior. That is, we wish to find a Gaussian o' such that g' = arg min, Da (pllg) in our projection step, instead of g' = arg min, KL. (pll4) as in EP. Here, the 'truc' distribution p is given by the tilted distribution p(z)= = (2))) (2) as usual. Typically it would be difficult. to compute or even approximale (112). However, we can work around this challenge by using the following property of a-divergences (Minka, 2005):
Proposition 1. Any stationary point q of (pllg) is also stationary point of proj [p()))((=)''-))
Algorithm 2 describes EP (with damping) using the a-divergence to find the approximate posterior, instead of the usual KL divergence minimisation. Setting p(æ) xx f())))((()), one can check by a straightforward calculation that proj (((=))((=====>>>))))))><)))))) the projection step in this modified EP reads
proj [[)())))))(())]] where we defined dara(ar) x (xr). Thus
d' = arg min (plll)
d((x)) = proj (x))))))
(117)
giving us (114) in step 6 of the algorithm.
We see from (117) that an interpretation of power EP is that it allows for the fractional inclusion of a true factor f., which then yields partial update to the corresponding approximate fractional factor of in step of the algorithm. In Appendix B, we detail the explicit computations of the projection step (117) for state-estimation problems, analogous the computations did in Section 3.4 for standard E.P.
18

---

Published in Transactions on Machine Learning Research (06/2022)
Algorithm 2 Damped Power EP for Dynamical Systems
Init: Set all factors of to (0, ool); Set g(z1) = p(z1) and marginals ((#1>1) NN (0
repeat
for t= 1 to T do
4:
for all factors L' where i = D., A, do Compute cavity distribution:
5:
(ari(e)) x qf
(113)
6:
Project:
d(x))) = proj (xx) (x)r(
(114)
Update fractional factor (with damping):
g g'(x))
q (at)
of (x+) new of (xt)²
x o8 (2tc)
(115)
((x;))
8:
Update posterior (with damping):
qf (2x) new
(2)))))))
1/a
( (xx()) g(xr))
3/a
g(x1) new x q((x))
= g(xx))
(116)
9, (xx))
end for 10: end for
11: until Convergence or maximum number of iterations exceeded
The advantage of applying damped power EP to our running example is demonstrated in Figure 3. We see that the RMSE/NLL are both lower than that the damped EP. The same choice of damping y and a is convenient mathomatically, since the power and damping cancel each other out, which simplifics the posterior update step 8 in Algorithm 2.
5 Results
In the following, we evaluate the proposed iterative state estimator three benchmark problems. We the root mcan square error (RMSE) and the negative log-likelihood (NLL) evaluate its performance, whose precise expressions can be found in Appendix A.3.
5.1 Non-linear Growth Model
First consider the non-linear benchmark from Section 2.3. The example shown in Figure 3 uses the unscented transform for single seed. We now repeat the same experiment with 10 different random seeds to evaluate the robustness of different algorithms. Our first goal is to evaluate power EP for different values of power a on the non-linear benchmark problem using different transforms.
5.1.1 Effects of Damping, Power and Linearisation Methods
We test a-divergence EP with three different linearisation techniques, the Taylor transform (TT), the (scaled) unscented transform (UT) and the Monte Carlo transform (MCT). For the scaled UT we use parameters (OVT, BUTT, KUT) (1, 2, 3) for the transition function and (out, KUT) = 21, 2, for the measurement function throughout the paper. Similar to the TT (employed by the EKS) and the UT (employed by the UKS), the MCT computes Gaussian approximations; however, the means and covariances of these Gaussian approximations are computed by sampling from the (Gaussian) cavity distribution and then using a Monte Carlo estimate to obtain the desired means and covariances. We refer to this the Monte-Carlo Kalman
19

---

Published in Transactions on Machine Learning Research (06/2022)
Power: 1.0, Damping: 1.0
Power: 1.0, Damping: 0.9
Power: 0.9, Damping: 0.9
MCT
RMSE
7.5
7.5
7.5
2.5
2.5-
2.5
10
20
40
50
10
30
40
10
20
50
400
TT
400
TT
TT
200
200
200
50
10
20
40
10
30
50
UT
MCT
UT
MCT
UT
MCT
10
5-
5-
5-
50
20
50
EP Iterations
EP Iterations
EP Iterations
(a) EP, UNCM
(b) Damped EP, UNCM
(c) Dampod power EP, UNGM
Figure 4: RMSE and NLL performance multiple runs of the experimental sct-up corresponding to Figure 3. Other linearisation methods, namely the Monte Carlo transform (MCT) and the first-order Taylor transform (TT) are shown alongside the unscented transform (UT). Solid lines show average performance across 10 independent evaluations and shaded regions show one standard deviation away from the mean. Dashed lines show average performance of the corresponding Gaussian smoother across the same 10 seeds. (a) (7,a)= = 1); (b) () = (1, (c) (7, (0.9, 0.9). MCT is the most stable and performs well for wide range of damping factors and powers. In contrast, the TT performs poorly on the NLL and is plotted separately from the UT and the MCT due to the difference in scales.
smoother (abbreviated as MCKS). We used 10 samples throughout this paper. The posteriors updated according to (102) (103) and (106) (107). The MCT determines the exact moments in the limit of infinitely many samples; that means, in the context of this article, the MCT serves as a apst-approximmtion baseline, which, however, is computationally more demanding than the TT or the UT.
In Figure 4a, we run standard approximate EP (a = 1.0, 7 = 1.0) using all three transforms for approximate inference (TT, UT, MCT). In this setting, the first. iteration is identical to the respective Gaussian smoother. We that both the TT and UT based iterations, in terms of RMSE and NLL, oscillate significantly; for the MCT, the RMSE improves rapidly in the beginning, however, subsequent updates do not improve the results much further; only slight oscillations occur. Note that the overall higher accuracy and stability of the MCT is due to the large number of samples (10*) used for statistical linearisation in this one-dimensional system.
Changing the damping rate to 7 = 0.9 (Figure 4b), as used in the running example, we see that there are fewer oscillations in the RMSE across all transformations (TT, UT, MCT). Further reductions oscillations can be observed by setting 7 vo smaller values, but at the cost of requiring more EP iterations for convergence. Note that damping does not necessarily improve the RMSE (or the variability of the results) significantly, compared with Figure4b 4a. An exception is the UT, where damping does improve the results slightly. However, the main impact of the damping factor on the RMSE is therefore to reduce oscillations. The same could be said for the NL.I in the UT and MCT examples. However, for the TT, damping appears to impact NL.L. negatively, increasing the variability of the results even further across seeds.
In Figure 4c, we use a 0.9 and damping ? = 0.9 as in Figure All transforms now show more stable RMSE curves with visibly more consistent behaviour for TT. The performance, in terms of RMSE, is slightly
20

---

Published in Transactions on Machine Learning Research (06/2022)
improved across all transformations compared to Figure 4b. For the NLL, the same observations hold with the exception of the TT, where the performance is generally worse compared to the previous two cases, with high spread across seeds after around iteration 10. Thus, with the exception of the TT, slightly decreasing the power and damping rate improves the overall results on average, which, for the UT, we found to be consistent across all random seeds used this experiment. For the MCT, this true for nine out of ten soeds that we tested, and when improvement is observed, it more pronounced compared to the UT (sec, e.g., Figure 10 in Appendix C).
In Table 2, we display the RMSE and the NLL results of damped power EP after 50 iterations for different combinations of a and y. This is compared against the results of the corresponding Kalman-type smoother for the given type of linearisation (e.g. the EKS when using the Taylor transform for linearisation) and the IEKS (Bell, 1994), which we take as baseline methods. Note that the IEKS is a MAP estimator (Bell, 1994) and therefore only produces pointwise predictions (s+)?=1- To evaluate the NLI., we supplement these predictions with covariance matrices (P.).-1: obtained by the RTS smoother at every iteration. We see in general that, EP iterations improve predictions in terms of the RMSE and NLL. However, setting the a parameter too low could result in bad NLL values, which we discuss in more details in the next section. The IEKS bascline does not produce good prodictions in this example; damped power EP overall does better, although there are some exceptions, e.g. the NLL values when a is small.
Taylor transform RMSE
Unscented transform
Monte-Carde transform
Baselines
NLL
RMSE
NLL
RMSE
NLL
EKS/UKS/MCKS
8.73 + 1.99
45.8 + 22.9
7.60 + 0.73
9.54 + 3.28
5.17 + 0.69
3.35 + 0.94
IEKS
15.7 = 4.0
98.8 = 22.8
2 = 1.0
RMSE
NLL
RMSE
NLL
RMSE
NLL
Power EP (a 1.0)
7.58 1.37
91.8 + 52.5
5.74 + 0.90
6.63 2.13
4.17 + 0.65
2.84 + 0.80
Power EP (a 0.8)
6.68 + 1.47
166 + 190
5.07 + 0.99
7.58 + 3.45
3.63 + 0.80
2.81 0.97
Power EP (a = 0.6)
6.87 1.49
97.8 = 87.3
5.69 1.11
9.39 4.77
3.17 + 0.84
3.52 1.44
Power EP (a 0.4)
7.26 1.62
200 + 223
5.92 1.38
43.0
32.8
3.04 + 0.84
Power EP (a 0.2)
7.28 4 1.27
136 + 114
6.01 + 1.40
197 + 161
3.81 + 1.0G
189 + 127
= 0.8
RMSE
NLL
RMSE
NLL
RMSE
NLL
Power EP (a = 1.0)
7.70 + 1.75
91.8
85.9
4.90 + 0.95
4.47
4.14 0.74
2.72 + 0.66
Power EP (a 0.8)
7.31 + 1.08
5.10 38.0
3.00 6.59
3.33 0.60
2.41 0.75
Power EP (a = 0.6)
7.65 + 1.06
93.9 63.9
5.26 + 0.98
9.14
2.94 + 0.70
4.00 + 4.66
Power EP (a = 0.4)
7.71 1.23
180 117
4.97 1.03
28.2 + 20.8
2.86 + 0.73
31.0 34.5
Power EP (a 0.2)
8.08 + 1.29
130 + 94
1.80
463 - 474
3.30 0.47
258 + 176
0.6
RMSE
NLL
RMSE
NLL
RMSE
NLL
Power EP (a = 1.0)
8.29 1.68
104 85.2
5.05 1 0.83
5.20
2.25
4.11 0.71
2.68 + 0.65
Power EP (a 0.8)
7.53
109 = 79.1
4.70
5.00 2.11
3.41 0.77
2.58 1.04
Power EP (a 0.6)
7.80 1.13
91.22
5.25 + 1.13
7.70 +
2.96 + 0.70
3.67 4.06
Power EP (a = 0.4)
7.89 + 1.06
+
5.14 0.99
21.6 + 16.0
2.74 0.55
2 35.8
Power EP (a = 0.2)
8.15 + 1.25
157 1
5.56 + 0.81
336 + 272
2.80 + 0.62
247 + 183
Table 2: Comparison of RMSE and NLL performance of itcrated damped power EP for the UNGM exper- iment. We compare the results of damped power EP with (0.2, 0.6, 0.8, 1.0) and 7 {0.6, 0.8, 1.0) against the corresponding Gaussian smoother. We also consider the IEKS as baseline. The mean and standard deviation is displayed across 10 different seeds after 50 EP iterations (we also take 50 iterations for the IEKS), for all transformation types (TT, UT and MCT). We highlight in bold the best results for the RMSE and the NLL. In general, iterating EP improves the RMSE and the NLL, provided a is not too small. An exception is the NLL when using the TT for linearisation, which case iterating EP gives worse results. The IEKS performed poorly in this example and power EP yields better results in general.
5.1.2 Sensitivity Analysis
For a more detailed sensitivity analysis, we perform a full sweep across both power a [0.1, 1.0] and damping 7 E [0.1, 1.0), with a fixed number of 50 EP iterations per experiment. We plot the resulting heat maps of the mean of both RMSE and NLL across ten different seeds and all three transforms in Figures 6 and 7.
21

---

Published in Transactions on Machine Learning Research (06/2022)
15-
Iteration: 1, Damping: 0.5
Iteration: 10, Damping: 0.5
15
Iteration: 50, Damping: 0.5
MCT
π
UT
MCT
MCT
RMSE
10
10
5
5
0.1
0.4
0.7
1.0
0.1
0.4
0.7
1.0
0.4
0.7
1.0
MCT
750-
MCT
200
NL1
500
250
0:
0.1
100
0.1
0.4
0.7
0.1
0.1
Power
Power
Power
(a) After EP iteration
(b) After 10 EP iterations
(c) Afte EP iterations
Figure 5: Sensitivity analysis of the proposed algorithm across different choices of power a. The damping is constant for all the experiments and set to 7 = 0.5. We average across 10 dilferent seeds with shaded areas indicating the standard deviation. (a) For the first iteration, we see that the results are largely insensitive to power, with RMSE slightly preferring lower powers. (b) and (c) show the corresponding 10th and 50th iteration, respectively. For the UT and the MCT, the RMSE slightly better at lower powers, whereas the NLL much worse (and increasing with iterations). The TT less sensitive to power overall, except in the 10th iteration, there sudden jump in the RMSE and NLL at a 1.0.
The RMSE in Figure 6a for the Taylor transform shows optimal performance in the top-right corner, for moderately high values of power and damping a, E [0.8, 1.0). Decreasing both parameters leads to increas- ingly worse performance, especially at very low damping rates of around E [031, The NLL (Figure 7a) is less sensitive to the power, but is similarly affected negatively by decreasing damping rates.
The worst performance, in terms of both RMSE and NLL, mutually occurs under the combination of low damping rates y E (0.1,0.3 and high power 7 E [0.9, 1.0). This due slow convergence at low damping rates and can be checked to converge to better values by iterating sufficiently many times. This is shown in Figure 12 (Appendix C), where we plot the RMSE and NLL values of damped EP after 500 iterations. The damping factor is set to 7 0.1 and the TT was used for linearisation. It is somewhat surprising that the RMSE and NLL eventually converge to reasonable values, despite reaching values several orders of magnitude higher (O(10)) for the RMSE and O)10") for the NLL!) during the itorations.
The RMSE plot (Figure 6a) also shows that, while generally prefers high power and damping rates, fixing a 1.0 can affect the performance negatively. This does not occur for the NLL (Figure 7a), which, on the other hand, sees sub-optimal performance when the damping rate is set to 1.0. This implics that when using the TT for this example, the optimal performance can be achieved by setting both and y close to, but not equal to 1.0.
For the UT (Figures 6b and 7b), the optimal performance in terms of RMSE occurs for a wider range of parameters compared to the TT, at around a, E (0.4, 1.0), and the worst performance occurs at the corner case a = 0.1, 0.1. Note that the NLL. value is still fairly modest this corner case despite performing badly at other values of damping, with fixed a = 0.1. This can be explained by the slow convergence when the damping factor is small. We see that the NLL is overall insensitive to damping. However, it is sensitive to power. A low power (around < 0.3) harms the performance substantially. This is due to the fact that in the low-power region, power EP is mode-sccking and gives over-confident predictions at incorrect modes. This can be seen for example in Figure Sa, where we plot the outputs of power EP with a = 0.1 and 7 0.4. In the figure, we see that in the 10th and 50th iterations, the posterior of the latent states does not cover the ground truth well at several places due to overconfident predictions. Contrast this with Figure 8b, which is the same plot except with a 0.9, where much of the ground truth is covered and therefore gives much better NLL score.
22

---

Published in Transactions on Machine Learning Research (06/2022)
1.0
>0
1.0
>0
1.0
>0
8
Damping
Compleg
0.6
0.6
0.6
6
6
0.4
0.4
0.4
4
0.2
0.2
0.2
0.2
0.4
0.6
0.8
1.0
0.4
0.6
0.8
1.0
Power
Power
Power
(a) Taylor transform, UNGM
(b) Unscented transform, UNGM
(c) Monte Carlo transform, UNGM
Figure 6: RMSE heatmap for the Uniform Non-Linear Growth model experiment. We fixed the number of iterations to 50 and plot the average RMSE over 10 independent runs of the proposed method for varying powers and damping factors. We vary both power and damping between (0.1,1.0) for different. transformation types. Optimal performance is achieved for mid-to-high damping and power factors. Very low damping rates yicld poor performance due to slow convergence.
1.0
>500
1.0
>500
>500
0.8
400
0.8
400
0.8
400
Damping
300
Damping
300
Damping
300
0.6
0.6
200
200
200
100
100
100
0.2
0.2
0.2
0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
Power
Power
Power
(a) Taylor transform, UNGM
(b) Unscented transform, UNGM
(c) Monte Carlo transform, UNGM
Figure 7: Negative Log-Likelihodd (NLL) heatmap for the UNGM experiment. Again, we fixed the number of iterations to and plot the average NLL over independent runs. Similar to the RMSE, the best performance occurs for mid-to-high power and damping factors. However in contrast to the RMSE, worst performance in the UT and the MCT occur at very low powers. The poor performance in the bottom-right corner in the MCT is due to slow convergence.
For the MCT, Figure 7e shows that the NLL results are similar to the ones obtained with the UT, with low power around a < 0.3 being harmful to the performance. In contrast, the result in Figure 6c, while subtle, indicates that the RMSE slightly prefers mid-to-low power at around a < 0.6. This indicated more clcarly in Figure 5, where we plottod the RMSE and NLL against the power different EP iterations. The significant degradation in NLL due to low power however is clearly not worth the slight gain in RMSE. Hence, it is optimal to use mid-to- high power of around > 0.6 in this example. We see that both the RMSE and NLL results in Figures 6c and 7c are less sensitive to damping rates; only the RMSE is affected negatively by low damping rates of around a < 0.2, due to the slower convergence rate in this region.
5.2 Bearings-only Tracking of Turning Target
We now test our proposed method for classic bearings-only tracking problem that is commonly encountered in passive sensor environments. In this problem, we assume that single target is executing maneuvering turn, with four sensors measuring the angles between the target and the sensors, similar to the example
23

---

Published in Transactions on Machine Learning Research (06/2022)
(lteration 1) RMSE=8.11, NLL=3.47
(Iteration 1) RMSE=7.67, NLL=3.53
20-
Prediction
20
Prediction
0
0-
-20
-20
20
60
80
100
20
40
60
80
100
(Iteration 10) RMSE=6.00 08, NLL=38.86
(Iteration 10) RMSE=6.1
NLL=3.74
20
20
0
0-
-20
-20
20
60
80
100
20
60
80
100
(lteration 50) RMSE=6.17, 17, NLL=681.48
(lteration 50) RMSE=4.74, NLL=2.65
20
20
0
-20
20
20
60
80
100
20
60
80
100
Time steps
Time steps
(a) a = 0.1. While the RMSE values generally improve, the NLL values get progressively worse with iterations as a result. of the modo-soeking behaviour of power EP with a = 0.1.
(b) a = 0.9. The RMSE decreases with iterations. However, in contrast to the results (a), we see that the NLI decreases too due to better uncertainty coverage.
Figure 8: Outputs of damped power EP for single run of the UNGM experiment. The ground-truth trajectory is shown in red; the shaded blue areas indicate the standard deviation around the mean (solid 0.4. The unscented transform is
blue) of the posterior state estimate using EP with damping factor used for linearisation. We show the results at EP ilerations 1, 10 and 50 for (a) a 0.1; (b) 0.9.
considered in Särkkä & Hartikainen (2010). The dynamical system given by
2, = (w + w1-1, w.- ~ Q), y, = h(x:)) - Dt.
(118)
v. ~ N(O, R),
where the state is given by (I1, 21, x2, 222 w) E R such that the position of the target is (21, the corresponding velocities are & (#1, respectively, and w is the time-varying turning rate of the maneuver. The transition matrix A (w) takes the form
1
sin(wat)
l-cos(wDe)
0
cos(w
O
sin(W 0
A(w)
0
1-cos(wD1)
sin(war)
0
(119)
0
sin (t)
Cts(wat)
0
0
where At is the timestep between two consecutive measurements.
The measurement model h(x)
(N' (x)) (()) E R4 is given by
22 (02-2)) 8
h h(()) tan
i (1,2,3,4)
(120)
which measures the angle between the target and the i-th sensor s (of. sb) Following the same setup as in the target tracking example considered by Arasaratnam Haykin (2009), we choose the initial state distribution p(æ0) NoH,, Eo), where
H4 = (1000m, 300ms 1000m, Oms -1 -1
(121)
Σg ttaa((10m') 10m2s- 100m², 10m2s 10 1011100ardd rrdd3-2),
(122)
24

---

EKS Published in Transactions on Machine Learning Research (06/2022)
and choose model error covariance to be of the form Q (ilag(i M, 922t), where q1 = 0.1m²s
d2 = 1.75 x 10-1s-,, and
A A
M
(123)
Since the measurements Yr only depend on the position variables, the cross-corrclation in (123) is important for reconstructing the full latent state. The observation errors for the four sensors are assumed be independent and given by R r0), where v10 x 10 rrad for all {1, 2, 4}. We consider trajectories of length 50s and set the sampling frequency to At 1s. The sensors placed the four corners of the square s' (-2k,,, 20km), s (20km, 20km), s = (-20km, 20km), s (20km,-00km) and assumed to produce measurements al every time step.
Position (z)
Velocity (2)
Angular velocity (w)
TT
RMSE
NLL
RMSE
NLL
RMSE (x10-9)
NLL
43.1 + 5.5
11.0 + 1.2
11.1 2.4
11.9
1.44 0.35
-2.80
TEKS
24.8 + 1.8
8.33 + 0.26
7.10 2.33
5.01
0.81
1.26 0.33
-3.03 0.10
PEP (a 1.0)
43.1 5.5
11.0
2..1
2.3
1.44
-.28
PEP (a = 0.8)
3.8
9.82 + 2.14
11.9
1.26 0.31
-3.00
PEP (a = 0.6)
3.9
9.62
9.80 2.13
1.31
1.33
-.19
PEP (a 0.4)
36.7
9.79 2.12
12.1
1.26 0.33
-.19
PEP (a = 0.2)
36.7 5.9
9.61 + 0.82
9.78 + 2.11
12.2 + 3.3
1.26 + 0.33
UT
RMSE
NLL
RMSE
NLL
RMSE (x10 ")
NLL
UKS
26.5 1.9
8.58 0.49
7.64 1.80
5.66 1 0.42
1.26 0.28
3.04 1 0.07
PEP (a 1.0)
1.9
0.49
7.64 1.80
5.66
0.42
1.26 + 0.28
3.04
PEP (a = 0.8)
2.12
8.51
7.50 + 1.87
5.37 0.48
1.26 0.29
-3.04 0.08
PEP (Q 0.6)
2.1
8.50 0.36
7.49 1.86
0.44
1.26 0.29
-3.04
PEP (a 0.4)
26.2
8.50 1.50 0.36
7.48 + 1.85
-.46
1.26 0.29
-3.04
PEP (a = 0.2)
26.2 2.1
8.50 1 0.36
7.48 1.85
5.32 + 0.45
1.26 1 0.28
3.04 0.08
MCT
RMSE
NLL
RMSE
NLL
RMSE (x10-*)
NLL
MCKS
28.6 3.1
8.70
9.85 2.17
6.03 = 0.48
1.66 = 0.35
22.5
0.12
PEP (a = 1.0)
26.4 4 2.2
8.57 + 0.58
7.68 + 1.84
5.69 + 0.55
1.27 0.29
-3.03 + 0.08
PEP (a = 0.8)
26.1
4.440
7.58 1.89
5.40
-3.3 0.08
PEP (a 0.6)
2.2
8.50
7.49
0.55
1.27
-3.03
PEP {a 0.4)
2.2
8.50
1.93
5.34
1.26 0.30
-3.04
PEP (a = 0.2)
26.3 2.1
8.52 1522:00 0.43
15551.88
5.38 + 0.55
1.27 + 0.29
-3.03 0.08
Table 3: Comparison RMSE and NLL performance iterated power EP (abbreviated as PEP) for the bearings-only turning target tracking problem. We compare the results for a {0.2, 0.4, 0.8, 1.0]. The damping factor is fixed at 7 1.0 for all experiments; changing the damping factor did affect the results much. We display the mean and standard deviation across 10 different seeds after 50 E.P iterations, for all transformation types (TT, UT and MCT). As baselines, also display the results of the corresponding Gaussian smoother and the IEKS after 50 ilerations. For each variable, the best results in terms of the RMSE and the NLL are highlighted in bold. We observe that when a = 1.0, i.e., standard approximate EP, iteration does not improve upon the results of the corresponding Kalman-type smoother for the TT and UT (improvementt seen however for the MCT). Choosing smaller values of a can improve the result, which is most signficant the position variable (z) when using the TT. However overall, the best performance is obtained by the LEKS, with power EP performing competitively only in the variable.
We compare the performance of itcrated power EP for various values of a, including a = 1.0, which is equivalent to standard approximate EP. In all the experiments, we take 50 EP iterations, which we verified is enough to ensure convergence, and set the damping factor to 7 1.0 (i.e., no damping). No oscillatory bchaviour was obscrved during the EP iterations (except for the MCT, where oscillations occurred primarily due its stochastic nature; taking larger number of samples helped reduce the oscillations) and therefore damping was necessary to aid convergence. For reference, also compared the results against the
25

---

Published in Transactions on Machine Learning Research (06/2022)
corresponding Kalman-type smoother for the given type of lincarisation and the IEKS (Bell, 1994), which form our baseline methods.
Below are further details on the configurations used in the experiments
The parameters of the scaled UT were set to (aVT, BCT, KUT) = 0,, for both the transition and measurement functions, in accordance with Särkkä & Hartikainen (2010).
10a particles were used in all the MCT experiments.
We used the values E {0.2, 0.4, 0.6, 0.8, 1.0) for the power.
The results of the experiments are displayed in Table 3. For this system, we observe that simply iterating the approximate EP (i.e., power EP with a 1.0) does not improve the results of the corresponding Kalman- type smoother when we use the TT and UT for linearisation. This is true for both the RMSE and NLL performance across all variables (position, velocity, angular velocity). However, by taking smaller values of a, we overall gain in performance. For the TT, significant improvement is obtained in the RMSE for the position variable when taking S 0.8, with the best performance occurring at a {0.2, 0.4) although the results are very close when a E {0.6, 0.8). We also see improvements in the RMSE for the velocity and angular velocity variables by taking a S 0.8, although they are less pronounced than the former. The NLL also soes some improvement by taking smaller values of a, except for the velocity variable, where taking a I 0.6 appears to slightly deteriorate the NLL performance. This likely due overconfident predictions owing to the mode-seeking nature of power EP for small a. While similar improvements in the RMSE and NLL are observed for the UT, the difference is not significant (in fact, difference at all in the angular velocity variable). This is likely due to the UKS already performing very well on this task and therefore has smaller margins for improvements.
On the other hand, when using the MCT for linearisation, we see that most of the improvements in the RMSE and the NLL comes from iterating as opposed to taking smaller a, although laking smaller a can help squeeze out better performance. This is evident from the fact that the difference in values between the result from the Kalman smoother and EP after 50 iterations between the results of power EP with different a. Interestingly, a can help to reduce the oscillations during the EP iterations even with moderate number of particles (in this case, 10*). We show this in Figure 13 (Appendix C), where we plotted the average RMSE and NLL across all variables (position, velocity and angular velocity) against the number of EP iterations for 7) E {(1.0, 1.0), (1.0, 0.8), (0.8, 1.0)). When a = 1.0, we see that damping reduces the oscillations slightly in both the RMSE and NLL, however the oscillations in the NLL is still very large. Taking a 0.8 reduces this oscillation in the NLL dramatically, even when damping is applied (although the RMSE becomes slightly oscillatory again).
on average greater than the difference also observed that smaller values for
While power EP does give better results than the corresponding Gaussian smoother, it does not outperform the IEKS, which overall gives the best performance in this example. Only in the w variable does power EP yield competitive results.
Overall, we conclude that for this system, the IEKS performs the best, although power EP can give good results. We see that simply iterating the EP does not alter the smoother results when the TT and UT are used for linearisation. However, notable improvements can be observed by considering smaller values of a in power EP. When using the MCT, we see benefits from iterating the approximate EP, even when = 1.0. Morcover, further gains in performance can be obtained by taking smaller values of a. The results for the IEKS however are hard to beat even with these improvements, with power EP only slightly outperforming the IEKS in the w variable.
5.3 Lorenz 96 Model
In our final example, we test our method on the Lorenz 96 model (Lorenz, 1996), another commonly used benchmark for nonlinear filters/smoothers (Brajard et al., 2020; Fablet et al., 2021; Ott et al., 2004). The
26

---

Published in Transactions on Machine Learning Research (06/2022)
model is defined by a system of ODEs
dz;
(zi+1 - 1i-2)0;-1 Di F, i ,,..,,d d
di
(124)
where
T-1 = 1:
50 = Fd Ja+1 = F1-
Here, d d>4 4 is the dimension of the system and F is forcing constant, which we set to F 8. For our transition model f. we consider fourth-order Runge-Kutta discrctisation of system (124) with step size At = 0.05. The model error covariance is Q = 0.1I. For the measurement model, we consider the quadratic operator h(xx) = (xf....,,. =²)T E Rdd with error covariance I
Since we are free to choose the state dimension of the system, we investigate how our method performs with increasing dimensions. We evaluate this for systems with dimensions d (5, 10, 20, 40, 100, 100,, Our initial state Mo is obtained by spinning-up the model (124) using the fourth-order Runge-Kutta (RK4) discretisation, starting from the state (0.01, O, o,..., 0) E Rd and taking 99 RK4 timesteps with step size At = 0.05. We set the initial state distribution to p(zo) = =a (o: | I).
The RMSE and NLL results of power EP applicd to the Lorenz 96 model across various values of d is shown in Table 4. We used the values a E (0.1, 1.0} for the power factor and took 10 EP iterations with no damping for all the experiments. Oscillatory behaviour was not observed during the EP iterations and adding damping did not change the results by much. For the linearisation methods, we considered the TT and the UT, with the same UT paramcters (ouT, 9ur, KUT) that we used for the UNGM experiment (Section 5.1). We found the MCT with 10a particles to be numerically unstable at d 2 40 for some seed settings and therefore do not report results for MCT here. 5
We see that overall, iterating EP helps to improve both the NLL and RMSE performance over the corre- sponding Kalman-type smoothers. In contrast, the values of a do not affect the performance by much. For the TT, improvements by EP iteration observed for the RMSE across all values of d, even outperforming the IEKS. At low dimensions (d 10) we see similar improvement in the NLL. However, for d 2 20, EP iterations appear to affect the NLL results negatively for some values of Moreover, the obtained NLL values are worse by several orders of magnitude compared to those obtained for d S 10. To investigate this further, we plotted the results for individual seeds in Figure 9 with d = 20 and (a, 7) = 1.0). We see that in nine out of len cases, EP iterations actually improve the NLL performance. However, there is a single outlier case that starts out significantly worse than the rest in the first iteration and becomes worse with subsequent itcrations for the NLL (for the RMSE, itcration appears to always improve the results, which means that the uncertainty estimates become increasingly overconfident). These outlier scenarios that start out badly are therefore responsible for the significantly worse NLL values and large spread that we see for d > 20, even though in most cases, EP iterations actually improve the NLL.
For the UT, we see improvements in the RMSE and NLL with EP iterations for all experiments that we conducted. Furthermore, we see that the gain in performance due to iterations becomes increasingly pronounced with increasing d, especially for the NLL. For example, at d = 5, we only see an improvement in the NLL by a few decimal points upon iterating. However at d = 200,, we see a performance gain of 0(102) much more than what had anticipated. While less obvious, this trend of increasing performance gain with increasing dimensions also holds for the RMSE. In Figure 14 (Appendix C), we show Hovmöller diagrams illustrating the absolute errors for the predictions made by the Unscented Kalman smoother and the errors made by approximate EP with (a, = (1.0, 1.0) after 10 iterations. The state dimension is set to d d= 200. The improvements in predictions due to EP itcration is clear-the errors made by EP after 10 iterations is visibly smaller compared to the errors made after one iteration (ie.. the UKS).
Overall, we conclude that when using the UT for linearisation, EP iterations have the effect of improving predictions in this system, which becomes more significant at higher dimensions. Thus, our method may prove to be more beneficial when applied to higher-dimensional systems, but with the downside of increased computational cost. While the cost of our method only scales linearly in the number of EP iterations, the cost of rumning a single iteration is already expensive for nigher-dimensiona systems. Therefore, it may not be feasible in practice to run many iterations.
5 EExpermments with > 10a many particles would have required unjustified computational resources and time.
27

---

Published in Transactions on Machine Learning Research (06/2022)
Taylor transform
Unscented transform
== 5
RMSE
NLL
RMSE
NLL
EKS/UKS
0.43 0.04
-2.47 0.49
0.42 0.05
-2.61 0.38
IEKS
0.41 0.04
-2.69
0.40
Power EP (a 1.0)
0.40 0.03
-2.77 0.38
0.40 + 0.03
-2.78 0.36
Power EP (a 0.8)
0.40 0.03
- 0.37
0.40 = 0.03
-2.8 = 0.36 0.36
Power EP (a 0.1)
0.41 + 0.03
-2.77 0.37
0.40 + 0.03
-2.78 0.36
d= 10
RMSE
NLL
RMSE
NLL
EKS/UKS
0.60 + 0.04
-5.19 + 0.45
0.60 + 0.02
5.38 + 0.29
IEKS
0.58 + 0.03
-5.71 = 0.30
Power EP (a = 1.0)
0.56 +: 0.02
-5.86 + 0.31
0.56 + 0.02
- 5.85 + 0.29
Power EP (a = 0.8)
0.56 0.02
5.88 = 0.31
0.56 - 0.02
5.88 = 0.30
Power EP (a 0.1)
0.37 0.03
-5.84 0.35
0.56 + 0.02
-5.86
0.33
20
RMSE
NLL
RMSE
NLL
EKS/UKS
2.49 4.29
1543 + 4626
0.95 0.09
9.33
1.19
IEKS
2.37 4.22
1789 + 5390
Power EP (a 1.0)
2.25 1 4.26
1637 = 4943
0.84 - 0.05
-11.0
Power EP (a 0.8)
2.34 + 4.23
1741 + 5252
0.82 + 0.04
-11.3 0.5
Power EP (a = 0.1)
2.36 + 4.23
1792 + 5400
0.87 + 0.12
-00.5
d = 40
RMSE
NLL
RMSE
NLL
EKS/UKS
3.27 3.92
1462 + 3184
1.41 + 0.09
16.5 = 1.19
IEKS
3.15 3.85
1540 + 3388
Power EP (a = 1.0)
3.11 3.86
1424 + 3146
1.21 + 0.05
21.6 1.22
Power EP (a 0.8)
3.11 3.85
1495 + 3298
1.20 + 0.05
21.9 + 1.3
Power EP (a 0.1)
3.13 3.86
1544 + 3398
1.28 0.20
18.33 9.5
d 100
RMSE
NLL
RMSE
NLL
EKS/UKS
3.04 2.01
741 + 2064
2.54 0.08
22.1 2.31
IEKS
2.90 + 1.92
850 2410
Power EP (a 1.0)
2.85 + 1.94
741 +: 2125
1.98 + 0.08
- 49.2 + 2.35
Power EP (a = 0.8)
2.85 + 1.93
821 + 2342
1.97 0.05
-1.1 + 0.9
Power EP (a = 0.1)
2.88 + 1.93
849 + 2410
1.97 + 0.12
50.6 + 4.2
d = 200
RMSE
NLL
RMSE
NLL
EKS/UKS
5.09 2.93
1399 + 2374
4.47 0.17
8.72 + 4.13
IEKS
4.92 2.91
1753 + 3024
Power EP (a = 1.0)
4.84 2.96
1604 + 2833
2.95 0.05
89.0 1.67
Power EP (a 0.8)
4.83 2.96
1698 + 2968
2.80 0.04
-103
Power EP (a 0.1)
4.87 + 2.93
1752 + 3027 2.61 + 0.03
-11
Table 4: Comparison of RMSE and NLL performance of iterated power EP for the Lorenz 96 experiment across different dimensions d E {5, 10, 20, 40, 100, 200). We compare the results for E {0.1, 0.8, 1.0) and consider only the Taylor transform (TT) and the Unscented transform (TT). The Monte Carlo transform is nunerically unstable at high dimensions. We display the mean and standard deviation across 10 different seeds after 10 EP iterations with = 1.0. As baselines, we also display the results of the corresponding Gaussian smoother and the IEKS after 10 iterations. For each dimension, we highlight bold, the best results in the RMSE and the NLL. We see that in general, iterating EP improves the performance of the smoother, which can even be quile significant at high dimensions for the UT. An exception is the NLL performance at moderate-to-hige! dimensions when using the TT. The poor performance seen here is due to outlier results. If we exclude these, the result remains overall positive (see Figure 9).
6 Conclusion
We proposed an iterative state estimator for non-linear dynamical systems based Expectation Propagation (FP). We show that any Gaussian filter smoother be shown special case of the proposed method.
28

---

Published in Transactions on Machine Learning Research (06/2022)
15.3
17000
15.2
16000
15.1
RMSE
NU
1.5
50
1.0
2
4
6
8 10
2
4
8
10
EP Iterations
EP Iterations
Figure 9: RMSE and NLL performance the Lorenz 96 experiment where we used the TT for linearisation and set (a, 7) = (1.0, 1.0). The state dimension is set to = 20. We display the results for 10 independent seeds. For most seeds, we observe that EP iterations improve predictions. However, there is one outlier that starts off badly after the first iteration. In this case, while subsequent iterations improve the RMSE performance, it deteriorates the NLL performance, leading to the results seen in Table 4, where the NLL performance becomes worse average with iterations.
The standard derivative-based EP method is numerically unstable for non-linear systems (Seeger, 2005). We show that by explicitly accounting for the implicit linearisation typically used in Gaussian smoothers, use derivative-free EP updates, which closely resemble the generic Gaussian smoothing update equations, to build an iterative state estima The proposed method gets the best of both worlds, numerical stability of Gaussian state estimators and iterative improvements of the EP algorithm. We further improve the algorithm by using damping and power EP to improve convergence and performance.
Our experiments show the efficacy of EP for state estimation improve standard (non-iterative) Gaussian filtering and smoothing methods three benchmark problems. The iterative EP method still needs the user to choose the linearisation algorithm, which is identical to choosing the type of Gaussian filter in the standard setting. In the three examples, we showcase that iterating our approximate EP state-estimation algorithms (with the help of damping and power EP) generally lead to improved results in terms of the RMSE and NLL performance compared with their Kalman filter/mmoother counterpart. Further improvements are also observed by (1) applying damping, which aids in convergence when EP is highly oscillatory, as was the case in the UNGM experiment, and (2) replacing standard EP with power EP, which makes the predictions more mode secking. The latter particularly beneficial in the bearings-only tracking experiment, where simply taking EP iterations did not result in improved performance (when taking the TT and UT for linearisation). However, setting a to smaller values lead to improved results. We also saw in the Lorenz 96 experiment that iterating EP may be particularly effective for higher-dimenssonal systems as we saw much more significant improvements at higher dimensions when considering the UT for lincarisation. This may be potentially useful for data assimilation problems weather forecasting (Carrassi et al., 2018; Evensen, 2009), which deals with aatering//smoothin
in very high dimensions.
We have also compared our algorithm with another iterative smoother, namely the IEKS (Bell, 1994), which is an extension of the EKS that iteratively refines the linearisation point after every RTS smoother pass. While this can produce good results, sometimes even better than our EP algorithm as we saw in the bearings- only tracking example, these are merely MAP estimators and strictly speaking, do not predict uncertainlies (although in our experiments, we have supplemented their outputs with covariance matrices generated by the RTS smoother at every iteration). EP on the other hand, outputs uncertainty naturally and morcover refines them iteratively without double counting observations. We have seen that on complex tasks such the UNGM and the L96 experiments, power E.P generally outperforms the IEKS.
29

---

Published in Transactions on Machine Learning Research (06/2022)
Overall, EP-based state estimation generalises classical Gaussian smoothers (we proved that they are special case), but it allows for iterative refinement of the state estimates. Furthermore, the message-pasiie passing formulation lends itself to asynchronous and distributed state estimation; we leave this for future work.
References
Shun-ichi Amari. Cntrrential-Gemmttiiin Methods in Statistics, volume 28 of Lecture Notes in Statistics. Springer, 1985.
Brian D. O. Anderson and John B. Moore. Optimal Fittering. Dover Publications, 2005.
Tenkaran Arasaralnam and Simon Haykin. Cubature Kalman Filters. IEEE Transactions on Aulomatic
Control, 54(6):1254-1268, 2009.
Karl J. Aström. Introduction to Stochastic Control Theory. Dover Publications, 2006.
Timothy D. Barfoot. State Estimation for Robotics. Cambridge University Press, 2017.
Bradley M Bell. The iterated Kalman smoother as a Gauss Newton method. SIAM Journal on Optimization, 1994.
Bradley M. Bell and Frederick W. Cathey. The Iterated Kalman Filter Update Gauss-Newton Method.
IEEE Transactions on Automatic Control, 38(2):294-25) 297, 1993.
Dimitri P. Bertsekas. Dynamic Programming and Optimal Control, volume 1. Athena Scientific, 3rd edition, 2005.
Christopher M. Bishop. Pattern Recográtion and Machine Learning. Springer, 2006.
Julien Brajard, Alberto Carrassi, Mare Bocquet, and Laurent Bertino. Combining Data Assimilation and Machine Learning to Emulate Dynamical Model from Sparse and Noisy Observations: A Case Study with the Lorenz 96 Model. Journal of Computational Science, 44:101171, 2020.
Alberto Carrassi, Mare Bocquet, Laurent Bertino, and Geir Evensen. Data assimilation the geosciences: An overview methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change, 2018.
The Analytic Sciences Corporation. Applied Optimal Estimation. MIT Press, 1974.
Guillaume P. Dehacne and Simon Barthelmé. Bounding Errors of Expectation Propagation. In Advances in Newral Information Processing Systems, 2016.
Mare Deisenroth. Efficient Reinforcement Learning Using Gaussian Processes. KIT Scientifie Publishing, 2010.
Mare P. Deisenroth and Shakir Mohamed. Expectation Propagation in Gaussian Process Dynamical Systems. In Aduances in Neural Information Processing Systems, 2012.
Marc P. Deisenroth and Henrik Ohlsson. A General Perspoctive on Gaussian Filtering and Smoothing: Explaining Current and Deriving New Algorithms. In Proceedings of the American Control Conference, 2011.
Mare P. Deisenroth, Ryan D. Turner, Marco F. Huber, Uwe D. Hanebeck, and Carl E. Rasmussen. Robust Filtering and Smoothing with Gaussian Processes. IEEE Transactions Automatic Control, 57(7):1865 1871, 2012.
Arnaud Doucet, Simon J Godsill, and Christophe Andrieu. On Sequential Monte Carlo Sampling Methods for Bayesian Filtering. Stalisties and Compuling, 10:197-208 208, 2000.
Geir Evensen. Data Assimilation: The Ensemble Kalman Filter. Springer, 2009.
30

---

Published in Transactions on Machine Learning Research (06/2022)
Ronan Fablet, Bertrand Chapron, Lucas Drumetz, Etienne Mémin, Olivier Pannekoucke, and François Rousseau. Learning Variational Data Assimilation Models and Solvers. Journal of Advonces in Mod- eling Earth Systems, 2022S0222222 2021.
San Gultekin and John Paisley. Nonlinear Kalman Filtering With Divergence Minimization. IEEE Trans- actions on Signal Processing, 65(23):1 13, 2017.
Tom Heskes and Onno Zocter. Expectation Propagation for Approximate Inference in Dynamic Bayesian Networks. In Proceedings of the International Conference on Uncertointy in Artificial Intelligence, 2002.
Tom Heskes and Onno Zoeter. Extended Version: Expectation Propagation for Approximate Inference Dynamic Bayesian Networks. Technical report, University of Nijmegen, 2003.
Simon J. Julier, Jeffrey K. Uhlmann, and Hugh F. Durrant- Whyte. A New Method for the Nonlinear Transformation of Means and Covariances in Filters and Estimators. IEEE Transactions on Automatic Control, 457774477448, 2000.
Rudolf E. Kalman. A New Approach to Linear Filtering and Prediction Problems. Journal of Basic Engi- neering, 82(1):35-45) 45, 1960.
S. Lakshmivarahan and David J. Stensrud. Ensemble Kalman Filter: Application to Meteorological Data Assimilation. JEEE Control Systems, 29(3):34-46, 2009.
Tine Lefebvre, Herman Bruyninckx, and Joris De Schutter. Kalman Filters for Non-linear Systems: A Comparison of Performance. International Journal of Control, 77(7):639-66) 2004.
Edward N. Lorenz. Predictability; A problem Partly Solved. In Seminar on Predictability, 1996.
Peter S. Maybeck. Stochastic Models, Estimation, and Control, volume 141 of Mathematics in Science and Engineering. Academic Press, 1979.
Thomas P. Minka. A Family of Algorithmss for Approzimale Bayesion Inference. PhD thesis, Massachusetts Institute of Technology, Cambridge, MA, USA, 2001a.
Thomas Minka. Expectation Propagation for Approximate Bayesian Inference. In Proceedings of the Conference on Uncertainty in Artificial Intelligence, pp. 362 369, 2001b.
Thomas P. Minka. Power EP. Technical Report MSR-TR-2004-144 Microsoft Rescarch, 2004.
Thomas Minka. Divergence Measures and Message Passing. Technical report, Microsoft Research, 2005.
Thomas P. Minka. EP: A Quick Reference. Technical report, Microsoft. Rescarch, 2008.
Manfred Opper and Ole Winther. Tractable Approximations for Probabilistic Models: The Adaptive TAP Mean Field Approach. Physical Review Letlers, 86(17):5, 2001.
Edward Ott, Brian R. Hunt, Istvan Szunyogh, Alckscy Zimin, Eric J. Kostelich, Mattoo Corazza, Eu- genia Kalnay, D. J. Patil, and James A. Yorke. A Local Ensemble Kalman Filter for Atmospheric Data Assimilation. Tellus A: Dynamic Meteorology and Oceanography, 56(5):415 428, 2004.
Kaare Brandt Petersen and Michael Syskind Pedersen. The Matrix Cookbook, 2012.
Yuan Qi. Extending Expectation Propagation for Grophical Models. PhD thesis, Massachusetts Institute of Technology, 2005.
Yuan Qi and Thomas Minka. Expectation Propagation for Signal Detection in Flat-Fading Channels. In Proceedings of the IEEE International Symposium on Information Theory, 2003.
Matti Raitoharju, Lennart Svensson, Ángel García- Feréádeee and Robert Piché. Damped Posterior Linearization Filter. IEEE Signal Processing Letters, 25(4):536-544 540, 2018.
31

---

Published in Transactions on Machine Learning Research (06/2022)
H. E. Rauch, F. Tung, and C. T. Striebel. Maximum Likelihood Estimates of Linear Dynamic Systems. ATAA Journal, 10011445-1450, 1965.
Simo Särkkä. Unscented Rauch- Rauch-Tung-Sa Tung-Striebel Smoother. IEEE Tronsactions on Automatic Control, 53(3): 845-849, 2008.
Simo Särkkä. Bayesian Filtering and Smoothing. Cambridge University Press, 2013. ISBN 9781 139344203.
Simo Särkkä and Jouni Hartikainen. On Gaussian Optimal Smoothing of Non-linear State Space Models. IEEE Transactions on Aulomatic Control, 55:1938-199 1941, 2010.
Simo Särkkä, Jouní Hartikainen, Lennart Svensson, and Fredrik Sandblom. On the Relation Between Gaus- sian Process Quadratures and Sigma-Point Methods. orXiv:1504, 05994, 2015.
Matthias W. Seeger. Expectation Propagation for Exponential Families. Technical report, University of California Berkeley, 2005.
Kenneth D. Senne. Stochastic Processes and Filtering Theory, IEEE Transactions on Automatic Control, 17(5))752-753, 1972.
Sebastian Thrun, Wolfram Burgard, and Dieter Fox. Probabilistic Robotics. MIT Press, 2005.
Rudolph van der Merwe and Eric A. Wan. Sigma-Point Kalman Filters for Integrated Navigation. In Proceedings of the Annual Meeting of The Institute of Navigation, 2004.
Huaiyu Zhu and Richard Rohwer. Information Geometric Measurements of Generalisation. Technical report, Aston University, 1995.
32

---

Published in Transactions on Machine Learning Research (06/2022)
A Appendix: Some Useful Mathematical Identities
A.1 Identities for Gaussians
For random vector E RD, define the Gaussian probability density function
N( nr Σ) = (2x)
Σ exp (- 1 (x: p1) Σ (x 11)
(125)
where la E R. and E E RDD are the mean and covariance function, respectively.
Partial derivatives of a log-Gaussian with respect to its mean and covariance are (Petersen & Pedersen, 2012):
llg pe 23
(x - 14) - Σ-
(126)
0μ
i Nog pl, E)
T Σ Σ- Σ (x u)( μ1)
(127)
0Σ
2
For Gaussian distributions W(x A), N( b, B), their multiplication and division results in another (unnormalised) Gaussian. For multiplication, we obtain
N(x a, A)W b, B) cN(x c)
(128)
C= (A + B 1-1
(129)
c C-ACAA- -2 BB++Bl+ B-bb) c N(a b, A B),
(130)
(131)
where in (131), c not a random variable (but a constant); however, we use the N notation for the functional
form of c. The division of two Gaussians is given by
N( a, A) N( b, B) NN(x
(132)
B-11-- C= (A-1 ----------
(133)
c C(A -² BB-
(134)
Finally, raising Gaussian to the power of a > 0 gives
N(ola,a A).
(135)
A.2 Identities for Matrix Inversion
A useful identity for matrix inversion (as a special case of the Woodbury identity) is
4-++B-ly A(A + B
(136)
A A(A + B A.
(137)
A.3 Metrics
Given the ground truth trajectory (æx) and outputs (m, Σe) of a given state-estimation algorithm for time steps 1,..., T, evaluate its performance using the root mean square error (RMSE) and the negative
33

---

Published in Transactions on Machine Learning Research (06/2022)
log- ligllikli (NLL). For ae RD, they are given by
RMSE
Dt 1ell?
(138)
T
NLL T
Σ t=1
log (2 Ex),
(139)
log logN(a;lu Ex) - log(2T)
2 log E (2x - 14) 11) ", (x. " 14),
(140)
respectively. While the RMSE only tells us how close our predictive mean is the ground truth, the NLL also evaluates the quality of our predictive uncertainties it penalises heavily when the ground truth does not lie within some credible region around the predictive mean, or when the uncertainty is much larger than nccessary even if it captures the ground truth.
B Appendix: Projections for Power EP
Here, we derive the expressions for the projections proj (Zx)) that are necessary in the implemen- tation of power EP. Recall that we denoted by (x); == g(x,)// ( (zx) the cavity distribution with respect to the fractional factor of (æt). In our implementation of power EP, we keep track of the fractional factors (xx) instead of the full factors (xx). To make this explicit, we adopt the notation rx(2;)) := of (xx) below.
Forward Projection Recall that the forward factor at time step is given by
(x = / p(zz x
(141)
In practice, we compute the cavity distribution 41 (x from the fractional factor (=( by taking
( = q(x,,-))/ 1)rr'd (x1-1).
(142)
Following the arguments in Section 3.4, the message (141) can approximated by
(x,) Nax, M,-144 14 + u1- M, M, + +Q) = -N(x; [M )
(143)
where M,-1 and Vt-1 are matrix and vector respectively, obtained by implicit linearisation (Lemma 1). Thus, we find the following approximation to the fractional forward factor (x:):
To (Zx) C (xx) (135) N(x I |4f a ).).
(144)
Mcasurement Projection As in standard EP, we derive the measurement projection for power EP via the partition function, which in this case, reads
Za = f2 (1 = N (h(i)) R)) aa (
(145)
Let (a)a(())) (a me 2as Da DA.) By following the same procedure as in Section 3.4, we obtain the following approximation to Zz:
2.- N J,44 a.aa + u. J. Me J. + a RR) = N 14 ande D. Eo.)
(146)
where as before, we used the implicit linearisation y, = J, + u, + e. The gradients of this approximate log partition function can thus be computed explicitly, which can be checked be
dlog Za
(y, - -M, ( (Ea.sia)-IJ.
(147)
du a/i
V dlog 2a
1
V Vinl VT uis - -J, (T
(148)
d
2
34

---

Published in Transactions on Machine Learning Research (06/2022)
This gives the update formulas (i.e., the moments of proj (Z1) Maax ple El))
mi Mr aizl 2a V
m/s
(149)
MM made Eo J J (3/, - 14 1
(150)
Mi -i + Ki - )
(151)
Σ = Eqmilo Ea.mlo u/o V M/S - 2V E/a 1 Ennnio
(152)
Eamila - D Da.ela (2 J. )
(153)
Eq.ilo K K
(154)
where k = Da.yla (Dailla) the Kalman gain malrix and M Ea.avia = Eg Da.1A J.
Backward Projection Again, we use the derivatives of the log partition function compute the backward projection proj (2)))
(are]- (xre)- The partition function corresponding to this case reads
Z4 = / fa
= - (x1+))dx2+1 ) ga (o)))d
(155)
Using (131), we obtain the following useful identity, which allows us to move the power inside the integral
(ʃ p( (e, (o) (21+1)dxz4))) )
(f f(x)). Q) 0,0 (21+1)
(
(156)
(131)
[N(F(a))) u M4t1 Q 1//B a
21+1)
(157)
(135)
N(f(a) ( aif -Q + -'' [
(158)
(131)
/ p(a441 20) (2)+1)dx211
(159)
Thus, using Fubini's theorem to switch the order of integrals, we get
/ 1 (204 )dze++) ) Cay (ala(()))e
Z4 =
(160)
(159)
/ ( pl p(a1411/0) xx) g (e(+1)dx11 das (a dx
(161)
Fubini
[ plaa. qu (a;) Ri 1d
(162)
Now, using the implicit linearisation #++1 M... v. e, we get the following approximation to Zạ:
2g = / N
Mi al + Uu, M, M. 4 M. +a - ai (#++1)dx0++
(163)
(131)
10 0.6 DnD +a 141)). N( ptt1
(164)
This gives the following expressions for the gradients:
dlog Zo
V p
(n) (.b) ( a- - -
(165)
dui a/
1i4 dlog Ea 1
(Vula
un aid M. (EER + -1 - -MM.)
(166)
2
35

---

Published in Transactions on Machine Learning Research (06/2022)
Now set Va.p((111)) Ma M441 E and define ((x(+1) qa, 90 (2 (I++1), (x)+1)) whose mean and covariance we denote by Pu+1 and Σ++1 respectively. In the case a 1, the distribution can be checked to equivalent to the marginal ((#44)). Since we have M, = 21,41 (2) (2 (follows from (19)), the update formulas for backward messages (i.e., moments of proj (a) (ax)] = Na 14. E() can be checked to be
o.6
H4 = 14
+ Dal9
M, DaD 141 +
+a -=11) -1(1))1(((0)) (1/t1 - (.D))
-1241)
(167)
14 a
+ Li(
a.D)
(168)
Σ. 0010
-Eal M, (2 + - (1+1) ( 1
(169)
2. + ( E(41) L.
(170)
1,14- (2 is the smoother gain. This can be checked in a similar way that we derived the standard RTS smoother update equations from the derivative-based backward EP updates in Section 3.4.
where -
C Appendix: Additional Figures
6.0
EP 1-00y
4.5
EP (a L-
EP with damping (a y
EP with damping (a- 1--1.0,y-...
5.5
Power EP with damping y= 0.9)
4.0
Power EP with damping (a= 0.9, y=
RMSE
5.0
3.5
4.5
3.0
2.5
10
20 EP Iterations
30
40
50
10
30
40
50
EP Iterations
Figure 10: RMSE and NLL for EP (red), damped EP (green), and damped power EP (grey) for the UNGM experiment. The Monte Carlo transform (MCT) is used as the lincarisation method here. Damping and power improves the performance significantly and reduces the oscillation seen the NLL for standard EP.
300
1.0, y-
(a=110, y=
with damping 0.9)
250
EP with damping Lgg-1.,y-0..
11
Power EP with damping 0.9 0.9
Power EP with damping (a - 0.9, y 0.9)
RMSE
200
10-
150
9-
100
10
20
30
40
50
10
20
30
40
50
Iterations
EP Iterations
Figure 11: RMSE and NLL for EP (red), damped EP (green), and damped power EP (grey) for the UNGM experiment. The Taylor transform (TT) is used as the linearisation method here. Both damped EP and damped power EP reduce the oscillations seen in EP for the RMSE, with damped power EP even improving the performance. However, this is not reflected in the NLL, where power appears to degrade performance.
36

---

Published in Transactions on Machine Learning Research (06/2022)
Power: 1.0, Damping: 0.1
1e11
Power: 1.0, Damping: 0.1
RMSE
NLI
0.5
-0.5
-4-
100
200 300
400
500
100
200
300
400
500
EP Iterations
Iterations
Figure 12: RMSE and NLL performance over multiple runs of the UNGM experiment, recorded for 500 iterations. We only consider the corner case (a = 1.0, = 0.1) and used the Taylor transform (TT) for lincarisation. The solid lines show the mean across ten different seeds and the shaded areas indicate the standard deviation. Initially, both the RMSE and NLL increase rapidly for a couple of iterations. However, they peak around 25 30 iterations before starting to decrease. Eventually, they converge to good values after around 200 iterations for the RMSE and 100 iterations for the NLL.
Power: 1.0, Damping: 1.0
Power: 1.0, Damping: 0.8
Power: 0.8, Damping: 1.0
35
MCT
35
MCT
35
MCT
RMSE
30
25
25
25
10
30
10
20
40
20
50
20
MCT
20
MCT
20
MCT
NLL
10
10 20
30 EP Iterations
40
50
10
20
30
40
50
10
30 EP Iterations
20
40
50
EP Iterations
Figure 13: RMSE and NLL performance of the bearings-only turning target tracking problem, averaged over all variables (position, velocity, angular velocity). The MCT is used for lincarisation. Solid line shows average performance across 10 random seeds and shaded regions indicate standard deviation. Dashed line shows the result of the Monte Carlo Kalman smoother. The NLL is highly oscillatory when a 1.0. Damping does little to reduce this oscillation, however is easily controlled by setting 0.8.
37

---

Published in Transactions on Machine Learning Research (06/2022)
0
15
0
2.5
25
25
25
220
10
50
50
Components
75
75
75
-1.5
5
100
100
100
-1.0
125
125
125
150
150
0.5
175
-5
175
175
0.0
0
25
25
0
25
Time
Time
Time
(a) Ground truth (a simulation of the Lorenz 96 model with d 200).
(b) Absolute emrors for the predictions (JÆTT FFeea))- The left panel displays the result from the UKS. The right panel displays the result of EP after 10 iterations.
Figure 14: Hovmöller diagram depicting the time-evolution of each component in the Lorenz 96 system with 200 dimensions. In (a), we display the ground truth trajectory of a single simulation of the Lorenz 96 model. In (b), we display the absolute errors for the predictions made by the Unscented Kalman smoother (left) vs approximate EP after 10 iterations (right). We used the UT for the linearisation method in approximate EP. We overall improvements in prediction quality by iterating EP.
38


Reviews:
[{"review_id": "WwAzrfXE2c", "reviewer": "Reviewer_NB9x", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: \n- The paper is well-written overall (subject to minor improvements) and the mathematical results are correct to the best of my knowledge. \n\nWeakness:\n- I found the paper bit vague at the beginning on its objective: if the objective is filtering p(x_t| y_1,..,y_t) or smoothing p(x_t| y_1,..,y_T). There is a significant difference between the two. filtering is recursive, does not need it store the past observations, and give estimate of the state of the system at current time. Hence, it is useful and important for control. Smoothing requires storing all the observation from time zero to the current time, not scalable with the time horizon T, and give a better estimate of the state in the past which is not necessarily useful for control. So the paper needs to be more clear on its objective, and the presented motivation based on control does not seem reasonable.       \n\n- I am not sure why it was stated that implementing the Kalman smother iteratively is difficult. I found this article on iterative Kalman smoother [1]. This is important as it forms the correct baseline for comparison.  \n\n- Figure 1 should be explained better: how is RMSE calculated? (how many runs?) To which time step x_t does the reported rmse correspond to? \n\n- The paper should present approaches on how to make the algorithm scalable with the time horizon to be implemented in real-time manner. How to extend the approximate Gaussians from T to T+1. Can a time-window procedure be implemented? \n\n- I did not understand how the formulas for the factors in (50) and (70) are obtained. Is there an expression for the factors in (26) in terms of distributions? I would appreciate if the message passing is explained better for readers like me who are unfamiliar with the topic. \n\n- Minor: it will be better to only number equations that are referenced. \n\n[1] Bell, Bradley M. \"The iterated Kalman smoother as a Gauss–Newton method.\" SIAM Journal on Optimization 4.3 (1994): 626-636. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents an algorithm to approximate the smoothing distribution for a nonlinear stochastic dynamical system. The approximation is in terms of a Gaussian distribution. The algorithm approximates the smoothing distribution p(x_t|y_1,...y_T) with product of certain Gaussian factors, and these factors are updated through an iterative message passing procedure going from t=0 to t=T. It is shown that a single iteration of the algorithm is equivalent to the Kalman smoother, and the algorithm can incorporate any linearization procedure. An extension of the algorithm with damping effects is also presented to remove the oscillatory behavior. The algorithm is illustrated with several numerical experiments in comparison with Kalman smoother.     ", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "My suggestions are reflected in the weakness points written above. ", "broader_impact_concerns": "I did not observe any broader impact concerns.  ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "toRJEW0QyV", "reviewer": "Reviewer_JrZp", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: The paper uses implicit linearization ideas to propose a generalization of nonlinear filters, which enables an improvement over standard \"one-pass\" methods. This is a welcome contribution. The derivation could be clarified (see below) but I found the idea important and useful.\n\nWeaknesses: The notation seems slightly confusing and the numerical simulations can be improved.\n\n\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes an expectation propagation method that generalizes standard nonlinear filtering techniques. Interestingly, one iteration of the proposed method recovers some of the well-known nonlinear filters, while running many iterations has the promise of improving the properties of the estimators. It could be an interesting contribution to the filtering techniques.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "(i) While the authors seem to use the standard EP notation, I found it slightly confusing. It could be introduced better if the authors somehow find a way to introduce the proper conditioning to the variational distributions (e.g. what's called \"marginal\" q(x_t) is often conditioned on y_{1:t} or y_{1:T} and this is not clear from the context easily). Also while for the linear case, the correspondences in Table 1 are true, for the nonlinear case, these are approximations and maybe could be made clearer.\n\n(ii) I urge authors to study a nonlinear example where things are more complex and conduct a sensitivity analysis for increasing dimension. This does not necessarily need to be a positive story but it could be important how much improvement is made, if one runs the EKF (i.e. a one iteration of this method) or iterate according to EP updates in a high-dimensional scenario. A Lorenz 96 example often used in the literature could be a good example and can be used to test unscented Kalman filter with a linear observation model. Note that this system allows one to test things with increasing dimension. Any other linear (but not independent) high dimensional system can also be tried like in the bearings-only scenario. Results with varying (increasing) d would be crucial for us to understand the usability of this methodology in a realistic setup.", "broader_impact_concerns": "There is no immediate ethical concerns regarding this work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "3ekDquR2rW", "reviewer": "Reviewer_NPfV", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n * The paper is clearly written.\n * The method is simple to use and can use many of the already existing techniques.\n * More iterations of EP often lead to better results.\n\nWeaknesses:\n\n * There is not that much novelty in the proposed approach, which is based on already known techniques. However, this is not a critical point for this journal.\n * The experimental section is weak since it is based on only two problems that are synthetic.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper describes the use of approximate expectation propagation (EP) to carry out approximate inference in non-linear dynamical systems. The method consists in approximating non-Gaussian distributions using a Gaussian distribution. This is done in an iterative process by using expectation propagation approximate factors. This factors are found by minimizing the more general alpha-divergence, instead of the typical KL divergence. The algorithm, after one iteration, is similar to already existing techniques. Performing more iterations is expected to lead to improvements, similar to those obtained of EP over ADF. This is verified in two experiments involving synthetic systems where the ground truth is known.  The results obtained indicate that damping the EP updates is needed for stability and that minimizing alpha-divergences can sometimes lead to better results.\n", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Eq. (28) should be propto.\n\nIn Fig. 3 damping reduces oscillations, but does not provide convergence. Why is this? Is it due to the approximate EP updates? Does increasing the amount of damping improve more convergence? From my experience, non-convergent EP is very rare and by using enough damping EP always seems to converge in practice.\n\nIt is not clear what is the computational cost of the proposed method.\n\nApproximate EP updates where the exact normalization constant of the factor times the cavity is also proposed and use in this work which could be relevant:\n\nBui, T. D., Hernández-Lobato, D., Hernández-Lobato, J. M., Li, Y., and Turner, R. E. (2016). Deep Gaussian\nprocesses for regression using approximate expectation propagation. In International Conference on\nMachine Learning, pages 1472–1481.\n\nI am curious to know if a more principled approach in which the cavity distributions are all the same, as considered in that work, could be also applied to this setting. That could have the advantage of leaving aside the pesky EP update equations and could allow to use standard optimization methods. May be the authors, can comment on that.\n\nUsing only 50 iterations for EP at most seems like a small number. The authors also indicate that and in the appendix they use up to 500. Is there a reason for not using more than 50 iterations? I wonder if the authors can use something to check convergence. May be looking at changes in the approximate factors from one iteration to another?\n\nI think that overall this is a nice paper that will be interesting for the community. However, the weakest part is the experimental section. I would suggest that the authors try to strengthen it by including at least another problem.\n", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers agreed that the contribution is interesting and relevant to the TMLR community. A common concern/suggestion about the experiments was addressed. After the rebuttal period, all reviewers were unanimous in their decision to accept this submission. 
Decision: Accept (poster)