arXiv:1207.2940v5 [stat.ML] 17 Aug 2016

Expectation Propagation in
Gaussian Process Dynamical Systems:
Extended Version
Marc Peter Deisenroth∗
Department of Computing
Imperial College London
UK
m.deisenroth@imperial.ac.uk

Shakir Mohamed∗
Google Deepmind
London
UK
shakir@google.com

Abstract
Rich and complex time-series data, such as those generated from engineering systems,
financial markets, videos or neural recordings, are now a common feature of modern data
analysis. Explaining the phenomena underlying these diverse data sets requires flexible
and accurate models. In this paper, we promote Gaussian process dynamical systems
(GPDS) as a rich model class that is appropriate for such analysis. In particular, we
present a message passing algorithm for approximate inference in GPDSs based on expectation propagation. By posing inference as a general message passing problem, we iterate
forward-backward smoothing. Thus, we obtain more accurate posterior distributions over
latent structures, resulting in improved predictive performance compared to state-of-theart GPDS smoothers, which are special cases of our general message passing algorithm.
Hence, we provide a unifying approach within which to contextualize message passing in
GPDSs.

1

Introduction

The Kalman filter and its extensions [1], such as the extended and unscented Kalman filters [8], are principled statistical models that have been widely used for some of the most
challenging and mission-critical applications in automatic control, robotics, machine learning, and economics. Indeed, wherever complex time-series are found, Kalman filters have
been successfully applied for Bayesian state estimation. However, in practice, time series often have an unknown dynamical structure, and they are high dimensional and noisy, violating
many of the assumptions made in established approaches for state estimation. In this paper,
we look beyond traditional linear dynamical systems and advance the state-of the-art in state
estimation by developing novel inference algorithms for the class of nonlinear Gaussian process
dynamical systems (GPDS).
GPDSs are non-parametric generalizations of state-space models that allow for inference
in time series, using Gaussian process (GP) probability distributions over nonlinear transition
and measurement dynamics. GPDSs are thus able to capture complex dynamical structure
with few assumptions, making them of broad interest. This interest has sparked the development of general approaches for filtering and smoothing in GPDSs, such as [9, 4, 6]. In
this paper, we further develop inference algorithms for GPDSs and make the following contributions: (1) We develop an iterative local message passing framework for GPDSs based
*Authors contributed equally. Appeared in Advances in Neural Information Processing Systems 25, pp. 2609–2617,
2012 [5].

1

on Expectation Propagation (EP) [12, 11], which allows for refinement of the posterior distribution and, hence, improved inference. (2) We show that the general message-passing
framework recovers the EP updates for existing dynamical systems as a special case and expose the implicit modeling assumptions made in these models. We show that EP in GPDSs
encapsulates all GPDS forward-backward smoothers [6] as a special case and transforms them
into iterative algorithms yielding more accurate inference.

2

Gaussian Process Dynamical Systems

Gaussian process dynamical systems are a general class of discrete-time state-space models
with
xt = h(xt−1 ) + wt ,
z t = g(xt ) + v t ,

wt ∼ N (0, Q) ,

v t ∼ N (0, R) ,

h ∼ GP h ,

g ∼ GP g ,

(1)
(2)

where t = 1, . . . , T . Here, x ∈ RD is a latent state that evolves over time, and z ∈ RE , E ≥ D,
are measurements. We assume i.i.d. additive Gaussian system noise w and measurement
noise v. The central feature of this model class is that both the measurement function g and
the transition function h are not explicitly known or parametrically specified, but instead
described by probability distributions over these functions. The function distributions are
non-parametric Gaussian processes (GPs), and we write h ∼ GP h and g ∼ GP g , respectively.
A GP is a probability distribution p(f ) over functions f that is specified by a mean function
µf and a covariance function kf [16]. Consider a set of training inputs X = [x1 , . . . , xn ]> and
2 ). The posterior
corresponding training targets y = [y1 , . . . yn ]> , yi = f (xi ) + w, w ∼ N (0, σw
predictive distribution at a test input x∗ is Gaussian distributed N (y∗ | µf (x∗ ), σf2 (x∗ )) with
−1
−1
k∗ , where k∗ = kf (X, x∗ ),
mean µf (x∗ ) = k>
y and variance σf2 (x∗ ) = k∗∗ − k>
∗K
∗K
k∗∗ = kf (x∗ , x∗ ), and K is the kernel matrix.
Since the GP is a non-parametric model, its use in GPDSs is desirable since it results in
fewer restrictive model assumptions, compared to dynamical systems based on parametric
function approximators for the transition and measurement functions (1)–(2). In this paper,
we assume that the GP models are trained, i.e., the training inputs and corresponding targets
as well as the GP hyperparameters are known. For both GP h and GP g in the GPDS, we used
zero prior mean functions. As covariance functions kh and kg we use squared- exponential
covariance functions with automatic relevance determination plus a noise covariance function
to account for the noise in (1)–(2).
Existing work for learning GPDSs includes the Gaussian process dynamical model (GPDM)
[21], which tackles the challenging task of analyzing human motion in (high-dimensional)
video sequences. More recently, variational [3] and EM-based [20] approaches for learning
GPDS were proposed. Exact Bayesian inference, i.e., filtering and smoothing, in GPDSs is analytically intractable because of the dependency of the states and measurements on previous
states through the nonlinearity of the GP. We thus make use of approximations to infer the
posterior distributions p(xt |Z) over latent states xt , t = 1, . . . , T , given a set of observations
Z = z 1:T . Existing approximate inference approaches for filtering and forward-backward
smoothing are based on either linearization, particle representations, or moment matching as
approximation strategies [9, 4, 6].
A principled incorporation of the posterior GP model uncertainty into inference in GPDSs
is necessary, but introduces additional uncertainty. In tracking problems where the location
of an object is not directly observed, this additional source of uncertainty can eventually
lead to losing track of the latent state. In this paper, we address this problem and propose
approximate message passing based on EP for more accurate inference. We will show that
forward-backward smoothing in GPDSs [6] benefits from the iterative refinement scheme of
2

qB(xt)

xt

qM(xt)

p(xt+1|xt)

xt+1

qB(xt)

qC(xt+1)

xt

qC(xt) qB(xt+1)

qM(xt)

qM(xt+1)

xt+1

qC(xt+1)

qM(xt+1)

Figure 1: Factor graph (left) and fully factored graph (right) of a general dynamical system.
EP, leading to more accurate posterior distributions over the latent state and, hence, to more
informative predictions and improved decision making.

3

Expectation Propagation in GPDS

Expectation Propagation [11, 12] is a widely-used deterministic algorithm for approximate
Bayesian inference that has been shown to be highly accurate in many problems, including sparse regression models [18], GP classification [10], and inference in dynamical systems [14, 7, 19]. EP is derived using a factor-graph, in which the distribution
Q over the latent
state p(xt |Z) is represented as the product of factors fi (xt ), i.e., p(xt |Z) = i fi (xt ). EP then
specifies an iterative
Q message passing algorithm in which p(xt |Z) is approximated by a distribution q(xt ) = i qi (xt ), using approximate messages qi (xt ). In EP, q and the messages qi
are members of the exponential family, and q is determined such that the the KL-divergence
KL(p||q) is minimized. EP is provably robust for log-concave messages [18] and invariant
under invertible variable transformations [17]. In practice, EP has been shown to be more
accurate than competing approximate inference methods [10, 18].
In the context of the dynamical system (1)–(2), we consider factor graphs of the form
of Fig. 1 with three types of messages: forward, backward, and measurement messages,
denoted by the symbols B, C, M, respectively. For EP inference, we assume a fully-factored
graph, using which we compute the marginal posterior distributions p(x1 |Z), . . . , p(xT |Z),
rather than the full joint distribution p(X|Z) = p(x1 , . . . , xT |Z). Both the states xt and
measurements z t are continuous variables and the messages qi are unnormalized Gaussians,
i.e., qi (xt ) = si N (xt | µi , Σi )

3.1

Implicit Linearizations Require Explicit Consideration

Alg. 1 describes the main steps of Gaussian EP for dynamical systems. For each node xt in the
fully-factored factor graph in Fig. 1, EP computes three messages: a forward, backward, and
measurement message, denoted by qB (xt ), qC (xt ), and qM (xt ), respectively. The EP algorithm
updates the marginal q(xt ) and the messages qi (xt ) in three steps. First, the cavity distribution
q \i (xt ) is computed (step 5 in Alg. 1) by removing qi (xt ) from the marginal q(xt ). Second, in
the projection step, the moments of fi (xt )q \i (xt ) are computed (step 6), where fi is the true
factor. In the exponential family, the required moments can be computed using the derivatives
of the log-partition function (normalizing constant) log Zi of fi (xt )q \i (xt ) [11, 12, 13]. Third,
the moments of the marginal q(xt ) are set to the moments of fi (xt )q \i (xt ), and the message
qi (xt ) is updated (step 7). We apply this procedure repeatedly to all latent states xt , t =
1, . . . , T , until convergence.
EP does not directly fit a Gaussian approximation qi to the non-Gaussian factor fi . Instead,
EP determines the moments of qi in the context of the cavity distribution such that qi =
proj[fi q \i ]/q \i , where proj[·] is the projection operator, returning the moments of its argument.
To update the posterior q(xt ) and the messages qi (xt ), EP computes the log-partition
function log Zi in (4) to complete the projection step. However, for nonlinear transition and

3

Algorithm 1 Gaussian EP for Dynamical Systems
1:
2:
3:
4:
5:

Init: Set all factors qi to N (0, ∞I); Set q(x1 ) = p(x1 ) and marginals q(xt6=1 ) = N (0, 1010 I)
repeat
for t = 1 to T do
for all factors qi (xt ), where i = B, M, C do
Compute cavity distribution q \i (xt ) = q(xt )/qi (xt ) = N (xt | µ\i , Σ\i ) with
−1 −1
Σ\i = (Σ−1
,
t − Σi )

6:

−1
µ\i = Σ\i (Σ−1
t µt − Σi µi )

Determine moments of fi (xt )q \i (xt ), e.g., via the derivatives of
log Zi (µ\i , Σ\i ) = log ∫ fi (xt )q \i (xt )dxt

7:

(3)

(4)

Update the posterior q(xt ) ∝ N (xt | µt , Σt ) and the approximate factor qi (xt ):
µt = µ\i + Σ\i ∇>
m,
∇m := d log Zi /dµ\i ,

\i
Σt = Σ\i − Σ\i (∇>
m ∇m − 2∇s )Σ

∇s := d log Zi /dΣ\i

\i

(5)
(6)
(7)

qi (xt ) = q(xt )/q (xt )
end for
end for
10: until Convergence or maximum number of iterations exceeded
8:
9:

measurement models in (1)–(2), computing Zi involves solving integrals of the form
Z
Z
p(a) = p(a|xt )p(xt )dxt = N (a | m(xt ), S(xt ))N (xt | b, B)dxt ,

(8)

where a = z t for the measurement message, or a = xt+1 for the forward and backward
messages. In nonlinear dynamical systems m(xt ) is a nonlinear measurement or transition
function. In GPDSs, m(xt ) and S(xt ) are the corresponding predictive GP means and covariances, respectively, which are nonlinearly related to xt . Because of the nonlinear dependencies between a and xt , solving (8) is analytically intractable. We propose to approximate
p(a) by a Gaussian distribution N (a | µ̃, Σ̃). This Gaussian approximation is only correct for
a linear relationship a = J xt , where J is independent of xt . Hence, the Gaussian approximation is an implicit linearization of the functional relationship between a and xt , effectively
linearizing either the transition or the measurement models.
When computing EP updates using the derivatives ∇m and ∇s according to (5) it is
crucial to explicitly account for the implicit linearization assumption in the derivatives—
otherwise, the EP updates are inconsistent. For example, in the measurement and the backward message, we directly approximate the partition functions Zi , i ∈ {M, C} by Gaussians
i
Z̃i (a) = N (µ̃i , Σ̃ ). The consistent derivatives d(log Z̃i )/dµ\i and d(log Z̃i )/dΣ\i of Z̃i with
respect to the mean and covariance of the cavity distribution q are obtained by applying the
chain rule, such that
∇m =

d log Z̃i
dµ\i

=

∂ log Z̃i ∂ µ̃i
∂ µ̃i ∂µ\i

∇s =

d log Z̃i
dΣ\i

=

∂ log Z̃i ∂ Σ̃
i
\i
∂ Σ̃ ∂Σ

∂ µ̃i
∂µ\i

i

= J > ∈ RE×D ,

i

∂ Σ̃
∂Σ\i

i

= (a − µ̃i )> (Σ̃ )−1 J > ∈ R1×D ,


i −1 ∂ Σ̃i
Z̃i ∂ log Z̃i
= 12 ∂∂(log
−
(
Σ̃
)
∈ RD×D ,
µ̃i )> ∂ µ̃i
∂Σ\i
= J I4 J > ∈ RE×E×D×D ,

(9)
(10)
(11)

where I4 ∈ RD×D×D×D is an identity tensor. Note that with the implicit linear model a = J xt ,
i
the derivatives ∂ µ̃i /∂Σ\i and ∂ Σ̃ /∂µ\i vanish. Although we approximate Zi by a Gaussian
4

Z̃i , we are still free to choose a method of computing its mean µ̃i and covariance matrix
i
i
Σ̃ , which also influences the computation of J = ∂(µ̃i )/∂µ\i . However, even if µ̃i and Σ̃
i
are general functions of µ\i and Σ\i , the derivatives ∂ µ̃i /∂µ\i and ∂ Σ̃ /∂Σ\i must equal
i
the corresponding partial derivatives in (11), and ∂ µ̃i /∂Σ\i and ∂ Σ̃ /∂µ\i must be set to
0. Hence, the implicit linearization expressed by the Gaussian approximation Z̃i must be
explicitly taken into account in the derivatives to guarantee consistent EP updates.

3.2

Messages in Gaussian Process Dynamical Systems

We now describe each the messages needed for inference in GPDSs, and outline the approximations required to compute the partition function in (4). Updating a message requires a projection to compute the moments of the new posterior marginal q(xt ), followed by a Gaussian
division to update the message itself. For the projection step, we compute approximate parti\i
\i
tion functions Z̃i , where i ∈ {M, B, C}. Using the derivatives d log Z̃i /dµt and d log Z̃i /dΣt ,
we update the marginal q(xt ), see (5).
Measurement Message For the measurement message in a GPDS, the partition function is
Z
Z
\M
\M
\M
\M
(12)
ZM (µt , Σt ) = fM (xt )q\M (xt )dxt ∝ fM (xt )N (xt | µt , Σt )dxt ,
fM (xt ) = p(z t |xt ) = N (z t | µg (xt ), Σg (xt )),

(13)

where fM is the true measurement factor, and µg (xt ) and Σg (xt ) are the predictive mean and
covariance of the measurement GP GP g . In (12), we made it explicit that ZM depends on the
\M
\M
moments µt and Σt of the cavity distribution q\M (xt ). The integral in (12) is of the form
(8), but is intractable since solving it corresponds to a GP prediction with uncertain inputs
[16] which is no longer Gaussian. However, the mean and covariance of a Gaussian approximation Z̃M to ZM can be computed analytically: either using exact moment matching [15, 4],
or approximately by expected linearization of the posterior GP [9]; details are given in the
\M
\M
Appendix. The moments of Z̃M are also functions of the mean µt and variance Σt of the
cavity distribution. By taking the linearization assumption of the Gaussian approximation
into account explicitly (here, we implicitly linearize GP g ) when computing the derivatives,
the EP updates remain consistent, see Sec. 3.1.
Backward Message To update the backward message qC (xt ), we require the partition function
Z
Z
\C
\C
\C
\C
ZC (µt , Σt ) = fC (xt )q\C (xt )dxt ∝ fC (xt )N (xt | µt , Σt )dxt ,
(14)
Z
Z
fC (xt ) = p(xt+1 |xt )q\B (xt+1 )dxt+1 = N (xt+1 | µh (xt ), Σh (xt ))q\B (xt+1 )dxt+1 .
(15)
Here, the true factor fC (xt ) in (15) takes into account the coupling between xt and xt+1 , which
was lost in assuming the full factorization in Fig. 1. The predictive mean and covariance of
GP h are denoted µh (xt ) and Σh (xt ), respectively. Using (15) in (14) and reordering the
integration yields
Z
Z
\C
\C
ZC (µt , Σt ) ∝ q\B (xt+1 ) p(xt+1 |xt )q\C (xt )dxt dxt+1 .
(16)
\C

We approximate the inner integral in (16), which is of the form (8), by N (xt+1 | µ̃\C , Σ̃
\C

by moment matching [15], for instance. Note that µ̃
5

\C

and Σ̃

are functions of

\C
µt

)

and

\C

Σt . This Gaussian approximation implicitly linearizes GP h . Now, (16) can be computed
\C

\B

analytically, and we obtain a Gaussian approximation Z̃C = N (µt+1 | µ̃\C , Σ̃
that allows us to update the moments of q(xt ) and the message qC (xt ).

\B

+ Σt+1 ) of ZC

Forward Message Similarly, for the forward message, the projection step involves computing the partition function
Z
Z
\B
\B
\B
\B
ZB (µt , Σt ) = fB (xt )q\B (xt )dxt = fB (xt )N (xt | µt , Σt )dxt ,
(17)
Z
Z
fB (xt ) = p(xt |xt−1 )q\C (xt−1 )dxt−1 = N (xt | µh (xt−1 ), Σh (xt−1 ))q\C (xt−1 )dxt−1 ,
where the true factor fB (xt ) takes into account the coupling between xt−1 and xt , see Fig. 1.
Here, the true factor fB (xt ) is of the form (8). We propose to approximate fB (xt ) directly by
B
a Gaussian qB (xt ) ∝ N (µ̃B , Σ̃ ). This approximation implicitly linearizes GP h . We obtain
the updated posterior q(xt ) by Gaussian multiplication, i.e., q(xt ) ∝ qB (xt )q\B (xt ). With this
approximation we do not update the forward message in context, i.e., the true factor fB (xt ) is
directly approximated instead of the product fB (xt )q\B (xt ), which can result in suboptimal
approximation.

3.3

EP Updates for General Gaussian Smoothers

We can interpret the EP computations in the context of classical Gaussian filtering and
smoothing [1]. During the forward sweep, the marginal q(xt ) = q\C (xt ) corresponds to the
filter distribution p(xt |z 1:t ). Moreover, the cavity distribution q\M (xt ) corresponds to the time
update p(xt |z 1:t−1 ). In the backward sweep, the marginal q(xt ) is the smoothing distribution
p(xt |Z), incorporating the measurements of the entire time series. The mean and covariance
of Z̃C can be interpreted as the mean and covariance of the time update p(xt+1 |z 1:t ).
Updating the moments of the posterior q(xt ) via the derivatives of the log-partition function recovers exactly the standard Gaussian EP updates in dynamical systems described by
Qi and Minka [14]. For example, when incorporating an updated measurement message, the
\M
\M
\M
zx\M
,
moments in (5) can also be written as µt = µt + K(z t − µz ) and Σt = Σt − KΣt
xz\M
\M \M
xz\M
\M −1
\M
respectively, where Σt
= cov[xt , z t ] and K = Σt
(Σz ) . Here, µz = E[g(xt )] and
\M
Σz = cov[g(xt )] + R, where xt ∼ q\M (xt ). Similarly, the updated moments of q(xt ) with
\C

a new backward message via (5) correspond to the updates [14] µt = µt
\C

\C

\C

\C

\C

\C

+ L(µt+1 − µt+1 )

and Σt = Σt + L(Σt+1 − Σt+1 )L> , where L = cov[xt , xt+1 ](Σt+1 )−1 . Here, we defined
\C

\C

µt+1 = E[h(xt )] and Σt+1 = cov[h(xt )] + Q, where xt ∼ q\C (xt ).
The iterative message-passing algorithm in Alg. 1 provides an EP-based generalization
and a unifying view of existing approaches for smoothing in dynamical systems, e.g., (Extended/Unscented/Cubature) Kalman smoothing and the corresponding GPDS smoothers [6].
Computing the messages via the derivatives of the approximate log-partition functions log Z̃i
recovers not only standard EP updates in dynamical systems [14], but also the standard
Kalman smoothing updates [1].
Using any prediction method (e.g., unscented transformation, linearization), we can compute Gaussian approximations of (8). This influences the computation of log Z̃i and its derivatives with respect to the moments of the cavity distribution, see (9)–(10). Hence, our messagepassing formulation is also general as it includes all conceivable Gaussian filters/smoothers
in (GP)DSs, solely depending on the prediction technique used.

6

Table 1: Performance comparison on the synthetic data set. Lower values are better.
NLLx
MAEx
NLLz

4

EKS
−2.04 ± 0.07
0.03 ± 2.0 × 10−3
−0.69 ± 0.11

EP-EKS
−2.17 ± 0.04
0.03 ± 2.0 × 10−3
−0.73 ± 0.11

GPEKS
−1.67 ± 0.22
0.04 ± 4.6 × 10−2
−0.75 ± 0.08

EP-GPEKS
−1.87 ± 0.14
0.04 ± 4.6 × 10−2
−0.81 ± 0.07

GPADS
+ 1.67 ± 0.37
1.79 ± 0.21
1.93 ± 0.28

EP-GPADS
−1.91 ± 0.10
0.04 ± 4 × 10−3
−0.77 ± 0.07

Experimental Results

We evaluated our proposed EP-based message passing algorithm on three data sets: a synthetic data set, a low-dimensional simulated mechanical system with control inputs, and a
high-dimensional motion-capture data set. We compared to existing state-of-the-art forwardbackward smoothers in GPDSs, specifically the GPEKS [9], which is based on the expected
linearization of the GP models, and the GPADS [6], which uses moment-matching. We refer
to our EP generalizations of these methods as EP-GPEKS and EP-GPADS.
In all our experiments, we evaluated the inference methods using test sequences of measurements Z = [z 1 , . . . , z T ]. We report the negative log-likelihood of predicted measurements
using the observed test sequence (NLLz ). Whenever available, we also compared the inferred
posterior distribution q(X) ≈ p(X|Z) of the latent states with the underlying ground truth
using the average negative log-likelihood (NLLx ) and Mean Absolute Errors (MAEx ). We
terminated EP after 100 iterations or when the average norms of the differences of the means
and covariances of q(X) in two subsequent EP iterations were smaller than 10−6 .

4.1

Synthetic Data

We considered the nonlinear dynamical system
xt+1 = 4 sin(xt ) + w ,

w ∼ N (0, 0.12 ) ,

zt = 4 sin(xt ) + v ,

v ∼ N (0, 0.12 ) .

We used p(x1 ) = N (0, 1) as a prior on the initial latent state. We assumed access to the latent
state and trained the dynamics and measurement GPs using 30 randomly generated points,
resulting in a model with a substantial amount of posterior model uncertainty. The length of
the test trajectory used was T = 20 time steps.
Tab. 1 reports the quality of the inferred posterior distributions of the latent state trajectories using the average NLLx , MAEx , and NLLz (with standard errors), averaged over 10
independent scenarios. For this dataset, we also compared to the Extended Kalman Smoother
(EKS) and an EP-iterated EKS (EP-EKS), as models which make use of the known dynamics.
Iterated forward-backward smoothing with EP (EP-EKS, EP-GPEKS, EP-GPADS) improved
the smoothing posteriors using a single sweep only (EKS, GPEKS, GPADS). The GPADS
had poor performance across all our evaluation criteria for two reasons: First, the GPs were
trained using few data points, resulting in posterior distributions with a high degree of uncertainty. Second, predictive variances using moment-matching are generally conservative
and increased the uncertainty even further. This uncertainty caused the GPADS to quickly
lose track of the period of the state, as shown in Fig. 2(a). By iterating forward-backward
smoothing using EP (EP-GPADS), the posteriors p(xt |Z) were iteratively refined, and the latent state could be followed closely as indicated by both the small blue error bars in Fig. 2(a)
and all performance measures in Tab. 1. EP smoothing typically required a small number of
iterations for the inferred posterior distribution to closely track the true state, Fig. 2(b). On
average, EP required fewer than 10 iterations to converge to a good solution in which the
mean of the latent-state posterior closely matched the ground truth.

7

Average NLL per data point

5
Latent State

2

True state
Posterior state distribution (EP−GPADS)
Posterior state distribution (GPADS)

0

−5
2

4

6

8

10
12
Time step

14

16

18

20

1
GPADS
EP−GPADS

0
−1
−2
5

10

15
EP iteration

20

25

30

(a) Example trajectory distributions with 95% confi- (b) Average NLLx as a function of the EP iteration with
dence bounds.
twice the standard error.

Figure 2: (a) Posterior latent state distributions using EP-GPADS (blue) and the GPADS (gray).
The ground truth is shown in red (dashed). The GPADS quickly loses track of the period of
the state revealed by the large posterior uncertainty. EP with moment matching (EP-GPADS)
in the GPDS iteratively refines the GPADS posterior and can closely follow the true latent
state trajectory. (b) Average NLLx per data point in latent space with standard errors of the
posterior state distributions computed by the GPADS and the EP-GPADS as a function of EP
iterations.

4.2

Pendulum Tracking

We considered a pendulum tracking problem to demonstrate GPDS inference in multidimensional settings, as well as the ability to handle control inputs. The state x of the system is
given by the angle φ measured from being upright and the angular velocity φ̇. The pendulum used has a mass of 1 kg and a length of 1 m, and random torques u ∈ [−2, 2] Nm
were applied for a duration 200 ms (zero-order-hold control). The system noise covariance
was set to Σw = diag(0.32 , 0.12 ). The state was measured indirectly by two bearings sensors with coordinates (x1 , y1 ) = (−2, 0) and (x
respectively, according to
 2 , y2 ) = (−0.5, −0.5),
sin φ−yi 
z = [z1 , z2 ]> + v , v ∼ N 0, diag(0.12 , 0.052 ) with zi = arctan cos
φ−xi , i = 1, 2. We trained
the GP models using 4 randomly generated trajectories of length T = 20 time steps, starting
from an initial state distribution p(x1 ) = N (0, diag(π 2 /162 , 0.52 )) around the upright position.
For testing, we generated 12 random trajectories starting from p(x1 ).
Tab. 2 summarizes the perfor- Table 2: Performance comparison on the pendulummance of the various inference swing data. Lower values are better.
methods.
Generally, the (EPNLLx
MAEx
NLLz
)GPADS performed better than the
GPEKS
−0.35 ± 0.39 0.30 ± 0.02 −2.41 ± 0.047
(EP-)GPEKS across all performance
−0.33 ± 0.44 0.31 ± 0.02 −2.39 ± 0.038
EP-GPEKS
measures. This indicates that the
GPADS
−0.80 ± 0.06 0.30 ± 0.02 −2.37 ± 0.042
(EP-)GPEKS suffered from overEP-GPADS −0.85 ± 0.05 0.29 ± 0.02 −2.40 ± 0.037
confident posteriors compared to
(EP-)GPADS, which is especially pronounced in the degrading NLLx values with increasing EP iterations and the relatively high standard errors. In about 20% of the test cases, the
inference methods based on explicit linearization of the posterior mean function (GPEKS and
EP-GPEKS) ran into numerical problems typical of linearizations [6], i.e., overconfident posterior distributions that caused numerical problems. We excluded these runs from the results
in Tab. 2. The inference algorithms based on moment matching (GPADS and EP-GPADS)
were numerically stable as their predictions are typically more coherent due to conservative
approximations of moment matching.

8

Figure 3: Latent space posterior distribution (95% confidence ellipsoids) of a test trajectory
of the golf-swing motion capture data. The further the ellipsoids are separated the faster the
movement.

4.3

Motion Capture Data

We considered motion capture data (from http://mocap.cs.cmu.edu/, subject 64) containing
10 trials of golf swings recorded at 120 Hz, which we subsampled to 20 Hz. After removing
observation dimensions with no variability we were left with observations z t ∈ R56 , which
were then whitened as a pre-processing step. For trials 1–7 (403 data points), we used the
GPDM [21] to learn MAP estimates of the latent states xt ∈ R3 . These estimated latent states
and their corresponding observations are used to train the GP models GP f and GP g . Trials
8–10 were used as test data without ground truth labels. The GPDM [21] focuses on learning
a GPDS; we are interested in good approximate inference in these models.
Fig. 3 shows the latent-state posterior distribution of a single test sequence (trial 10) obtained from the EP-GPADS. The most significant prediction errors in observed space occurred
in the region corresponding to the yellow/red ellipsoids, which is a low-dimensional embedding of the motion when the golf player hits the ball, i.e., the periods of high acceleration
(poses 3–5).
Tab. 3 summarizes the results of inference on the golf data set in all test trials: Iterating
forward-backward smoothing by means of EP improved the inferred posterior distributions
over the latent states. The posterior distributions in latent space inferred by the EP-GPEKS
were tighter than the ones inferred by the EP-GPADS. The NLLz -values suffered a bit from
this overconfidence, but the predictive performance of the EP-GPADS and EP-GPEKS were
similar. Generally, inference was more difficult in areas with fast movements (poses 3–5 in
Fig. 3) where training data were sparse.
The computational demand the Table 3: Average inference performance (NLLz , motion
two inference methods for GPDSs capture data set). Lower values are better.
we presented is vastly different.
Test trial GPEKS EP-GPEKS GPADS EP-GPADS
High-dimensional approximate inTrial 8
14.20
13.82
14.28
14.09
ference in the motion capture exTrial
9
15.63
14.71
15.19
14.84
ample using moment matching
Trial 10
26.68
25.73
25.64
25.42
(EP-GPADS) was about two orders
of magnitude slower than approximate inference based on linearization of the posterior GP
mean (EP-GPEKS): For updating the posterior and the messages for a single time slice, the
EP-GPEKS required less than 0.5 s, the EP-GPADS took about 20 s. Hence, numerical stability
and more coherent posterior inference with the EP-GPADS trade off against computational
demands.

9

5

Conclusion

We have presented an approximate message passing algorithm based on EP for improved inference and Bayesian state estimation in GP dynamical systems. Our message-passing formulation generalizes current inference methods in GPDSs to iterative forward-backward smoothing. This generalization allows for improved predictions and comprises existing methods for
inference in the wider theory for dynamical systems as a special case. Our new inference approach makes the full power of the GPDS model available for the study of complex time-series
data. Future work includes investigating alternatives to linearization and moment matching
when computing messages, and the more general problem of learning in Gaussian process
dynamical systems.

Acknowledgements
The research leading to these results has received funding from the European Community’s
Seventh Framework Programme (FP7/2007–2013) under grant agreement #270327 (CompLACS) and from the Canadian Institute for Advanced Research (CIFAR). We thank Zhikun
Wang for his help with the motion capture data set.

GP Predictions from Test Input Distributions

A

We will now review two approximations to the predictive distribution
ZZ
p(xt ) =
p(f (xt−1 )|xt−1 )p(xt−1 )df dxt−1 ,

(18)

where f ∼ GP and xt−1 ∼ N (µt−1 , Σt−1 ).

A.1

Moment Matching

In the moment-matching approach, we analytically compute the mean µt and the covariance
Σt of p(xt ). Using the law of iterated expectations, we obtain


µt = Ext−1 Ef [f (xt−1 )|xt−1 ] = Ext−1 [mf (xt−1 )] ,
(19)
where mf is the posterior mean function of the dynamics GP. For target dimension a =
1, . . . , D, we obtain
µat = q >
a βa ,

qai = √

σf2
|Σt−1 Λ−1
a +I|


−1
exp − 12 ν >
i (Σt−1 + Λa ) ν i ,

ν i := (xi − µt−1 )

(20)

for i = 1, . . . , n, where β a = K −1
a ya.
Using the law of iterated variances, the entries of Σt for target dimensions a, b = 1, . . . , D
are


2
σaa
= Ext−1 varf [∆a |xt−1 ] +Ef,xt−1 [∆2a ]−(µat )2 ,
(21)
2
σab
= Ef,xt−1 [∆a ∆b ]−µat µbt ,

a 6= b ,

(22)

2 do not contain an
respectively, where µat is known from (20). The off-diagonal terms σab
additional term Ext−1 [covf [∆a , ∆b |xt−1 ]] because of the conditional independence assumption
used for GP training: Target dimensions do not covary for a given xt−1 .

10

2 and σ 2 , we obtain
For the term common to both σaa
ab

Ef,xt−1 [∆a ∆b ] = β >
a Qβ b ,

Qij =

ka (xi ,µt−1 )kb (xj ,µt−1 )

√

|R|

exp

1 > −1
2 z ij R Σt−1 z ij



(23)

−1
−1
−1
with R := Σt−1 (Λ−1
a + Λb ) + I and z ij := Λa ν i + Λb ν j with ν i taken from (20). Hence,
2
the off-diagonal entries σab of Σt are fully determined by (20) and (22).
2 of Σ contain an additional term
From (21), we see that the diagonal entries σaa
t



2
(24)
Ext−1 varf [∆a |xt−1 ] = σf2a − Tr K −1
a Q + σwa

with Q given in (23). This concludes the computation of Σt .
The moment-matching approximation minimizes the KL divergence KL(p||q) between the
true distribution p and an approximate Gaussian distribution q. This is generally a conservative approximation, i.e., q has probability mass where p has mass [2].

A.2

Linearizing the GP Mean Function

An alternative way of approximating the predictive GP distribution for uncertain test inputs is
to linearize the posterior GP mean function [9]. This is equivalent to computing the expected
linearization of the GP distribution over functions. Given this linearized function, we apply
standard results for mapping Gaussian distributions through linear models. Linearizing the
posterior GP mean function yields to a predicted mean that corresponds to the posterior GP
mean function evaluated at the mean of the input distribution, i.e.,

rai = σf2a exp − 12 (xi − µt−1 )> Λ−1
µat = Ef [fa (µt−1 )] = r >
(25)
a βa ,
a (xi − µt−1 )
for i = 1, . . . , n and target dimensions a = 1, . . . , D, where β a = K −1
a y a . The covariance
matrix Σt of the GP prediction is
Σt = V Σt−1 V > + Σw ,

V =

∂µt
∂µt−1

∂r a
= β>
,
a ∂µ
t−1

(26)

where r a is given in (25) and V is the Jacobian evaluated at µt−1 . In (26), Σw is a diagonal
matrix whose entries are the model uncertainty plus the noise variance evaluated at µt−1 . This
means “model uncertainty” no longer depends on the density of the data points. Instead it is
assumed constant.
Using linearization, the approximation optimality in the KL sense of the moment matching is lost. However, especially in high dimensions, linearization is computationally more
beneficial. This speedup is largely due to the simplified treatment of model uncertainty.

References
[1] B. D. O. Anderson and J. B. Moore. Optimal Filtering. Dover Publications, Mineola, NY,
USA, 2005.
[2] C. M. Bishop. Pattern Recognition and Machine Learning. Information Science and Statistics.
Springer-Verlag, 2006.
[3] A. Damianou, M. K. Titsias, and N. D. Lawrence. Variational Gaussian Process Dynamical Systems. In Advances in Neural Information Processing Systems, 2011.
[4] M. P. Deisenroth, M. F. Huber, and U. D. Hanebeck. Analytic Moment-based Gaussian
Process Filtering. In L. Bouttou and M. L. Littman, editors, Proceedings of the 26th International Conference on Machine Learning, pages 225–232, Montreal, QC, Canada, June 2009.
Omnipress.
11

[5] M. P. Deisenroth and S. Mohamed. Expectation Propagation in Gaussian Process Dynamical Systems. In Advances in Neural Information Processing Systems, 2012.
[6] M. P. Deisenroth, R. Turner, M. Huber, U. D. Hanebeck, and C. E. Rasmussen. Robust
Filtering and Smoothing with Gaussian Processes. IEEE Transactions on Automatic Control,
57(7):1865–1871, 2012. doi:10.1109/TAC.2011.2179426.
[7] T. Heskes and O. Zoeter. Expectation Propagation for Approximate Inference in Dynamic Bayesian Networks. In A. Darwiche and N. Friedman, editors, Proceedings of the
International Conference on Uncertainty in Artificial Intelligence, pages 216–233, 2002.
[8] S. J. Julier and J. K. Uhlmann. Unscented Filtering and Nonlinear Estimation. Proceedings
of the IEEE, 92(3):401–422, March 2004.
[9] J. Ko and D. Fox. GP-BayesFilters: Bayesian Filtering using Gaussian Process Prediction
and Observation Models. Autonomous Robots, 27(1):75–90, July 2009.
[10] M. Kuss and C. E. Rasmussen. Assessing Approximate Inference for Binary Gaussian
Process Classification. Journal of Machine Learning Research, 6:1679–1704, December 2005.
[11] T. P. Minka. Expectation Propagation for Approximate Bayesian Inference. In J. S. Breese
and D. Koller, editors, Proceedings of the 17th Conference on Uncertainty in Artificial Intelligence, pages 362–369, Seattle, WA, USA, August 2001. Morgan Kaufman Publishers.
[12] T. P. Minka. A Family of Algorithms for Approximate Bayesian Inference. PhD thesis, Massachusetts Institute of Technology, Cambridge, MA, USA, January 2001.
[13] T. P. Minka. EP: A Quick Reference. 2008.
[14] Y. Qi and T. Minka. Expectation Propagation for Signal Detection in Flat-Fading Channels. In Proceedings of the IEEE International Symposium on Information Theory, 2003.
[15] J. Quiñonero-Candela, A. Girard, J. Larsen, and C. E. Rasmussen. Propagation of Uncertainty in Bayesian Kernel Models—Application to Multiple-Step Ahead Forecasting. In
IEEE International Conference on Acoustics, Speech and Signal Processing, volume 2, pages
701–704, April 2003.
[16] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. Adaptive
Computation and Machine Learning. The MIT Press, Cambridge, MA, USA, 2006.
[17] M. W. Seeger. Expectation Propagation for Exponential Families. Technical report, University of California Berkeley, 2005.
[18] M. W. Seeger. Bayesian Inference and Optimal Design for the Sparse Linear Model.
Journal of Machine Learning Research, 9:759–813, 2008.
[19] M. Toussaint and C. Goerick. From Motor Learning to Interaction Learning in Robotics,
chapter A Bayesian View on Motor Control and Planning, pages 227–252. SpringerVerlag, 2010.
[20] R. Turner, M. P. Deisenroth, and C. E. Rasmussen. State-Space Inference and Learning
with Gaussian Processes. In Y. W. Teh and M. Titterington, editors, Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, volume JMLR: W&CP
9, pages 868–875, May 2010.
[21] J. M. Wang, D. J. Fleet, and A. Hertzmann. Gaussian Process Dynamical Models for
Human Motion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 30(2):283–
298, 2008.

12

