Dec 17, 2013 - State-space models are successfully used in many areas of science, engineer- ing and economics to model time series and dynamical syste...

0 downloads 7 Views 879KB Size

arXiv:1306.2861v2 [stat.ML] 17 Dec 2013

Roger Frigola1 , Fredrik Lindsten2 , Thomas B. Sch¨on2,3 and Carl E. Rasmussen1 1. Dept. of Engineering, University of Cambridge, UK, {rf342,cer54}@cam.ac.uk 2. Div. of Automatic Control, Link¨oping University, Sweden, [email protected] 3. Dept. of Information Technology, Uppsala University, Sweden, [email protected]

Abstract State-space models are successfully used in many areas of science, engineering and economics to model time series and dynamical systems. We present a fully Bayesian approach to inference and learning (i.e. state estimation and system identification) in nonlinear nonparametric state-space models. We place a Gaussian process prior over the state transition dynamics, resulting in a flexible model able to capture complex dynamical phenomena. To enable efficient inference, we marginalize over the transition dynamics function and infer directly the joint smoothing distribution using specially tailored Particle Markov Chain Monte Carlo samplers. Once a sample from the smoothing distribution is computed, the state transition predictive distribution can be formulated analytically. Our approach preserves the full nonparametric expressivity of the model and can make use of sparse Gaussian processes to greatly reduce computational complexity.

1

Introduction

State-space models (SSMs) constitute a popular and general class of models in the context of time series and dynamical systems. Their main feature is the presence of a latent variable, the state xt ∈ X , Rnx , which condenses all aspects of the system that can have an impact on its future. A discrete-time SSM with nonlinear dynamics can be represented as xt+1 = f (xt , ut ) + vt , (1a) yt = g(xt , ut ) + et , (1b) where ut denotes a known external input, yt denotes the measurements, vt and et denote i.i.d. noises acting on the dynamics and the measurements, respectively. The function f encodes the dynamics and g describes the relationship between the observation and the unobserved states. We are primarily concerned with the problem of learning general nonlinear SSMs. The aim is to find a model that can adaptively increase its complexity when more data is available. To this effect, we employ a Bayesian nonparametric model for the dynamics (1a). This provides a flexible model that is not constrained by any limiting assumptions caused by postulating a particular functional form. More specifically, we place a Gaussian process (GP) prior [1] over the unknown function f . The resulting model is a generalization of the standard parametric SSM. The functional form of the observation model g is assumed to be known, possibly parameterized by a finite dimensional parameter. This is often a natural assumption, for instance in engineering applications where g corresponds to a sensor model – we typically know what the sensors are measuring, at least up to some unknown parameters. Furthermore, using too flexible models for both f and g can result in problems of non-identifiability. We adopt a fully Bayesian approach whereby we find a posterior distribution over all the latent entities of interest, namely the state transition function f , the hidden state trajectory x0:T , {xi }Ti=0 1

and any hyper-parameter θ of the model. This is in contrast with existing approaches for using GPs to model SSMs, which tend to model the GP using a finite set of target points, in effect making the model parametric [2]. Inferring the distribution over the state trajectory p(x0:T | y0:T , u0:T ) is an important problem in itself known as smoothing. We use a tailored particle Markov Chain Monte Carlo (PMCMC) algorithm [3] to efficiently sample from the smoothing distribution whilst marginalizing over the state transition function. This contrasts with conventional approaches to smoothing which require a fixed model of the transition dynamics. Once we have obtained an approximation of the smoothing distribution, with the dynamics of the model marginalized out, learning the function f is straightforward since its posterior is available in closed form given the state trajectory. Our only approximation is that of the sampling algorithm. We report very good mixing enabled by the use of recently developed PMCMC samplers [4] and the exact marginalization of the transition dynamics. There is by now a rich literature on GP-based SSMs. For instance, Deisenroth et al. [5, 6] presented refined approximation methods for filtering and smoothing for already learned GP dynamics and measurement functions. In fact, the method proposed in the present paper provides a vital component needed for these inference methods, namely that of learning the GP model in the first place. Turner et al. [2] applied the EM algorithm to obtain a maximum likelihood estimate of parametric models which had the form of GPs where both inputs and outputs were parameters to be optimized. This type of approach can be traced back to [7] where Ghahramani and Roweis applied EM to learn models based on radial basis functions. Wang et al. [8] learn a SSM with GPs by finding a MAP estimate of the latent variables and hyper-parameters. They apply the learning in cases where the dimension of the observation vector is much higher than that of the latent state in what becomes a form of dynamic dimensionality reduction. This procedure would have the risk of overfitting in the common situation where the state-space is high-dimensional and there is significant uncertainty in the smoothing distribution.

2

Gaussian Process State-Space Model

We describe the generative probabilistic model of the Gaussian process SSM (GP-SSM) represented in Figure 1b by f (xt ) ∼ GP mθx (xt ), kθx (xt , x0t ) , (2a) xt+1 | ft ∼ N (xt+1 | ft , Q), yt | xt ∼ p(yt | xt , θ y ),

(2b) (2c)

and x0 ∼ p(x0 ), where we avoid notational clutter by omitting the conditioning on the known inputs ut . In addition, we put a prior p(θ) over the various hyper-parameters θ = {θ x , θ y , Q}. Also, note that the measurement model (2c) and the prior on x0 can take any form since we do not rely on their properties for efficient inference. The GP is fully described by its mean function and its covariance function. An interesting property of the GP-SSM is that any a priori insight into the dynamics of the system can be readily encoded in the mean function. This is useful, since it is often possible to capture the main properties of the dynamics, e.g. by using a simple parametric model or a model based on first principles. Such

(a) Standard GP regression

(b) GP-SSM

Figure 1: Graphical models for standard GP regression and the GP-SSM model. The thick horizontal bars represent sets of fully connected nodes. 2

simple models may be insufficient on their own, but useful together with the GP-SSM, as the GP is flexible enough to model complex departures from the mean function. If no specific prior model is available, the linear mean function m(xt ) = xt is a good generic choice. Interestingly, the prior information encoded in this model will normally be more vague than the prior information encoded in parametric models. The measurement model (2c) implicitly contains the observation function g and the distribution of the i.i.d. measurement noise et .

3

Inference over States and Hyper-parameters

Direct learning of the function f in (2a) from input/output data {u0:T −1 , y0:T } is challenging since the states x0:T are not observed. Most (if not all) previous approaches attack this problem by reverting to a parametric representation of f which is learned alongside the states. We address this problem in a fundamentally different way by marginalizing out f , allowing us to respect the nonparametric nature of the model. A challenge with this approach is that marginalization of f will introduce dependencies across time for the state variables that lead to the loss of the Markovian structure of the state-process. However, recently developed inference methods, combining sequential Monte Carlo (SMC) and Markov chain Monte Carlo (MCMC) allow us to tackle this problem. We discuss marginalization of f in Section 3.1 and present the inference algorithms in Sections 3.2 and 3.3. 3.1

Marginalizing out the State Transition Function

Targeting the joint posterior distribution of the hyper-parameters, the latent states and the latent function f is problematic due to the strong dependencies between x0:T and f . We therefore marginalize the dynamical function from the model, and instead target the distribution p(θ, x0:T | y1:T ) (recall that conditioning on u0:T −1 is implicit). In the MCMC literature, this is referred to as collapsing [9]. Hence, we first need to find an expression for the marginal prior p(θ, x0:T ) = p(x0:T | θ)p(θ). Focusing on p(x0:T | θ) we note that, although this distribution is not Gaussian, it can be represented as a product of Gaussians. Omitting the dependence on θ in the notation, we obtain p(x1:T | θ, x0 ) =

T Y

p(xt | θ, x0:t−1 ) =

T Y

N xt | µt (x0:t−1 ), Σt (x0:t−1 ) ,

(3a)

t=1

t=1

with e −1 (x1:t−1 − m0:t−2 ), µt (x0:t−1 ) = mt−1 + Kt−1,0:t−2 K 0:t−2 −1 e e Σt (x0:t−1 ) = Kt−1 − Kt−1,0:t−2 K K> 0:t−2

t−1,0:t−2

(3b) (3c)

e 0 . Equation (3) follows from the fact that, once for t ≥ 2 and µ1 (x0 ) = m0 , Σ1 (x0 ) = K conditioned on x0:t−1 , a one-step prediction for the state variable is a standard GP prediction. Here, > we have defined the mean vector m0:t−1 , m(x0 )> . . . m(xt−1 )> and the (nx t) × (nx t) positive definite matrix K0:t−1 with block entries [K0:t−1 ]i,j = k(xi−1 , xj−1 ). We use two sets of e 0:t−1 = indices, as in Kt−1,0:t−2 , to refer to the off-diagonal blocks of K0:t−1 . We also define K K0:t−1 + It ⊗ Q. We can also express (3a) more succinctly as, e −1 (x1:t − m0:t−1 )). (4) e 0:t−1 |− 12 exp(− 1 (x1:t − m0:t−1 )> K p(x1:t | θ, x0 ) = |(2π)nx t K 0:t−1 2 This expression looks very much like a multivariate Gaussian density function. However, we eme 0:t−1 depend (nonlinearly) on the argument phasize that this is not the case since both m0:t−1 and K x1:t . In fact, (4) will typically be very far from Gaussian. 3.2

Sequential Monte Carlo

With the prior (4) in place, we now turn to posterior inference and we start by considering the joint smoothing distribution p(x0:T | θ, y0:T ). The sequential nature of the proposed model suggests the use of SMC. Though most well known for filtering in Markovian SSMs – see [10, 11] for an introduction – SMC is applicable also for non-Markovian latent variable models. We seek to api proximate the sequence of distributions p(x0:t | θ, y0:t ), for t = 0, . . . , T . Let {xi0:t−1 , wt−1 }N i=1 3

be a collection of weighted particles approximating p(x0:t−1 | θ, y0:t−1 ) by the empirical distribuPN i tion, pb(x0:t−1 | θ, y0:t−1 ) , i=1 wt−1 δxi0:t−1 (x0:t−1 ). Here, δz (x) is a point-mass located at z. To propagate this sample to time t, we introduce the auxiliary variables {ait }N i=1 , referred to as ancestor indices. The variable ait is the index of the ancestor particle at time t − 1, of particle xit . j Hence, xit is generated by first sampling ait with P(ait = j) = wt−1 . Then, xit is generated as, ai

t xit ∼ p(xt | θ, x0:t−1 , y0:t ),

(5) ait

for i = 1, . . . , N . The particle trajectories are then augmented according to xi0:t = {x0:t−1 , xit }. Sampling from the one-step predictive density is a simple (and sensible) choice, but we may also consider other proposal distributions. In the above formulation the resampling step is implicit and corresponds to sampling the ancestor indices (cf. the auxiliary particle filter, [12]). Finally, the particles are weighted according to the measurement model, wti ∝ p(yt | θ, xit ) for i = 1, . . . , N , where the weights are normalized to sum to 1. 3.3

Particle Markov Chain Monte Carlo

There are two shortcomings of SMC: (i) it does not handle inference over hyper-parameters; (ii) despite the fact that the sampler targets the joint smoothing distribution, it does in general not provide an accurate approximation of the full joint distribution due to path degeneracy. That is, the successive resampling steps cause the particle diversity to be very low for time points t far from the final time instant T . To address these issues, we propose to use a particle Markov chain Monte Carlo (PMCMC, [3, 13]) sampler. PMCMC relies on SMC to generate samples of the highly correlated state trajectory within an MCMC sampler. We employ a specific PMCMC sampler referred to as particle Gibbs with ancestor sampling (PGAS, [4]), given in Algorithm 1. PGAS uses Gibbs-like steps for the state trajectory x0:T and the hyper-parameters θ, respectively. That is, we sample first x0:T given θ, then θ given x0:T , etc. However, the full conditionals are not explicitly available. Instead, we draw samples from specially tailored Markov kernels, leaving these conditionals invariant. We address these steps in the subsequent sections. Algorithm 1 Particle Gibbs with ancestor sampling (PGAS) 1. Set θ[0] and x1:T [0] arbitrarily. 2. For ` ≥ 1 do (a) Draw θ[`] conditionally on x0:T [` − 1] and y0:T as discussed in Section 3.3.2. (b) Run CPF-AS (see [4]) targeting p(x0:T | θ[`], y0:T ), conditionally on x0:T [` − 1]. (c) Sample k with P(k = i) = wTi and set x1:T [`] = xk1:T . 3. end 3.3.1

Sampling the State Trajectories

To sample the state trajectory, PGAS makes use of an SMC-like procedure referred to as a conditional particle filter with ancestor sampling (CPF-AS). This approach is particularly suitable for non-Markovian latent variable models, as it relies only on a forward recursion (see [4]). The difference between a standard particle filter (PF) and the CPF-AS is that, for the latter, one particle at each e0:T = {e eT }. We then samtime step is specified a priori. Let these particles be denoted x x0 , . . . , x et . ple according to (5) only for i = 1, . . . , N − 1. The N th particle is set deterministically: xN t =x To be able to construct the N th particle trajectory, xN t has to be associated with an ancestor particle at time t − 1. This is done by sampling a value for the corresponding ancestor index aN t . Following [4], the ancestor sampling probabilities are computed as i i e t−1|T w ∝ wt−1

et:T }, y0:T ) et:T }) p({xi0:t−1 , x p({xi0:t−1 , x i i ∝ wt−1 = wt−1 p(e xt:T | xi0:t−1 ). (6) i p(x0:t−1 , y0:t−1 ) p(xi0:t−1 )

where the ratio is between the unnormalized target densities up to time T and up to time t − 1, respectively. The second proportionality follows from the mutual conditional independence of the et:T } refers to a path in XT +1 formed by concatenating observations, given the states. Here, {xi0:t−1 , x 4

the two partial trajectories. The above expression can be computed by using the prior over state i e t−1|T trajectories given by (4). The ancestor sampling weights {w }N i=1 are then normalized to sum j N to 1 and the ancestor index aN t is sampled with P(at = j) = wt−1|t .

The conditioning on a prespecified collection of particles implies an invariance property in CPF-AS, e0:T let x e00:T be generated as follows: which is key to our development. More precisely, given x e0:T . 1. Run CPF-AS from time t = 0 to time t = T , conditionally on x e00:T to one of the resulting particle trajectories according to P(e 2. Set x x00:T = xi0:T ) = wTi . e0:T ) on XT +1 , For any N ≥ 2, this procedure defines an ergodic Markov kernel MθN (e x00:T | x leaving the exact smoothing distribution p(x0:T | θ, y0:T ) invariant [4]. Note that this invariance holds for any N ≥ 2, i.e. the number of particles that are used only affect the mixing rate of the kernel MθN . However, it has been experienced in practice that the autocorrelation drops sharply as N increases [4, 14], and for many models a moderate N is enough to obtain a rapidly mixing kernel. 3.3.2

Sampling the Hyper-parameters

Next, we consider sampling the hyper-parameters given a state trajectory and sequence of observations, i.e. from p(θ | x0:T , y0:T ). In the following, we consider the common situation where there are distinct hyper-parameters for the likelihood p(y0:T | x0:T , θ y ) and for the prior over trajectories p(x0:T | θ x ). If the prior over the hyper-parameters factorizes between those two groups we obtain p(θ | x0:T , y0:T ) ∝ p(θ y | x0:T , y0:T ) p(θ x | x0:T ). We can thus proceed to sample the two groups of hyper-parameters independently. Sampling θ y will be straightforward in most cases, in particular if conjugate priors for the likelihood are used. Sampling θ x will, nevertheless, be harder since the covariance function hyper-parameters enter the expression in a non-trivial way. However, we note that once the state trajectory is fixed, we are left with a problem analogous to Gaussian process regression where x0:T −1 are the inputs, x1:T are the outputs and Q is the likelihood covariance matrix. Given that the latent dynamics can be marginalized out analytically, sampling the hyper-parameters with slice sampling is straightforward [15].

4

A Sparse GP-SSM Construction and Implementation Details

A naive implementation of the CPF-AS algorithm will give rise to O(T 4 ) computational complexity, since at each time step t = 1, . . . , T , a matrix of size T × T needs to be factorized. However, it is possible to update and reuse the factors from the previous time step, bringing the total computational complexity down to the familiar O(T 3 ). Furthermore, by introducing a sparse GP model, we can reduce the complexity to O(M 2 T ) where M T . In Section 4.1 we introduce the sparse GP model and in Section 4.2 we provide insight into the efficient implementation of both the vanilla GP and the sparse GP. 4.1

FIC Prior over the State Trajectory

An important alternative to GP-SSM is given by exchanging the vanilla GP prior over f for a sparse counterpart. We do not consider the resulting model to be an approximation to GP-SSM, it is still a GP-SSM, but with a different prior over functions. As a result we expect it to sometimes outperform its non-sparse version in the same way as it happens with their regression siblings [16]. Most sparse GP methods can be formulated in terms of a set of so called inducing variables [17]. These variables live in the space of the latent function and have a set I of corresponding inducing inputs. The assumption is that, conditionally on the inducing variables, the latent function values are mutually independent. Although the inducing variables are marginalized analytically – this is key for the model to remain nonparametric – the inducing inputs have to be chosen in such a way that they, informally speaking, cover the same region of the input space covered by the data. Crucially, in order to achieve computational gains, the number M of inducing variables is selected to be smaller than the original number of data points. In the following, we will use the fully independent conditional (FIC) sparse GP prior as defined in [17] due to its very good empirical performance [16]. As shown in [17], the FIC prior can be obtained by replacing the covariance function k(·, ·) by, k FIC (xi , xj ) = s(xi , xj ) + δij k(xi , xj ) − s(xi , xj ) , (7) 5

where s(xi , xj ) , k(xi , I)k(I, I)−1 k(I, xj ), δij is Kronecker’s delta and we use the convention whereby when k takes a set as one of its arguments it generates a matrix of covariances. Using the Woodbury matrix identity, we can express the one-step predictive density as in (3), with µFIC (x0:t−1 ) = mt−1 + Kt−1,I PKI,0:t−2 Λ−1 t 0:t−2 (x1:t−1 − m0:t−2 ), e t−1 − St−1 + Kt−1,I PKI,t−1 , ΣFIC (x0:t−1 ) = K t

(8a) (8b)

−1 KI,0:t−2 Λ−1 , 0:t−2 K0:t−2,I )

e 0:t−2 − S0:t−2 ] and SA,B , where P , (KI,I + Λ0:t−2 , diag[K −1 KA,I KI,I KI,B . Despite its apparent cumbersomeness, the computational complexity involved in computing the above mean and covariance is O(M 2 t), as opposed to O(t3 ) for (3). The same idea can be used to express (4) in a form which allows for efficient computation. Here diag refers to a block diagonalization if Q is not diagonal. We do not address the problem of choosing the inducing inputs, but note that one option is to use greedy methods (e.g. [18]). The fast forward selection algorithm is appealing due to its very low computational complexity [18]. Moreover, its potential drawback of interference between hyperparameter learning and active set selection is not an issue in our case since hyper-parameters will be fixed for a given run of the particle filter. 4.2

Implementation Details

As pointed out above, it is crucial to reuse computations across time to attain the O(T 3 ) or O(M 2 T ) computational complexity for the vanilla GP and the FIC prior, respectively. We start by discussing the vanilla GP and then briefly comment on the implementation aspects of FIC. There are two costly operations of the CPF-AS algorithm: (i) sampling from the prior (5), requiring the computation of (3b) and (3c) and (ii) evaluating the ancestor sampling probabilities according to (6). Both of these operations can be carried out efficiently by keeping track of a Cholesky faci e et:T −1 }) = Lit Li> torization of the matrix K({x t , for each particle i = 1, . . . , N . Here, 0:t−1 , x i e e 0:T −1 , but where the covariance function et:T −1 }) is a matrix defined analogously to K K({x0:t−1 , x i et:T −1 }. From Lit , it is possible to identify is evaluated for the concatenated state trajectory {x0:t−1 , x sub-matrices corresponding to the Cholesky factors for the covariance matrix Σt (xi0:t−1 ) as well as for the matrices needed to efficiently evaluate the ancestor sampling probabilities (6). It remains to find an efficient update of the Cholesky factor to obtain Lit+1 . As we move from time et will be replaced by xit in the concatenated trajectory. Hence, the t to t + 1 in the algorithm, x i i e e et+1:T −1 }) can be obtained from K({x et:T −1 }) by replacing nx rows and matrix K({x0:t , x 0:t−1 , x columns, corresponding to a rank 2nx update. It follows that we can compute Lit+1 by making nx successive rank one updates and downdates on Lit . In summary, all the operations at a specific time step can be done in O(T 2 ) computations, leading to a total computational complexity of O(T 3 ). For the FIC prior, a naive implementation will give rise to O(M 2 T 2 ) computational complexity. This can be reduced to O(M 2 T ) by keeping track of a factorization for the matrix P. However, to reach the O(M 2 T ) cost all intermediate operations scaling with T has to be avoided, requiring us to reuse not only the matrix factorizations, but also intermediate matrix-vector multiplications.

5

Learning the Dynamics

Algorithm 1 gives us a tool to compute p(x0:T , θ | y1:T ). We now discuss how this can be used to find an explicit model for f . The goal of learning the state transition dynamics is equivalent to that of obtaining a predictive distribution over f ∗ = f (x∗ ), evaluated at an arbitrary test point x∗ , Z p(f ∗ | x∗ , y1:T ) = p(f ∗ | x∗ , x0:T , θ) p(x0:T , θ | y1:T ) dx0:T dθ. (9) Using a sample-based approximation of p(x0:T , θ | y1:T ), this integral can be approximated by p(f ∗ | x∗ , y1:T ) ≈

L

L

`=1

`=1

1X 1X p(f ∗ | x∗ , x0:T [`], θ[`]) = N (f ∗ | µ` (x∗ ), Σ` (x∗ )), L L 6

(10)

where L is the number of samples and µ` (x∗ ) and Σ` (x∗ ) follow the expressions for the predictive distribution in standard GP regression if x0:T −1 [`] are treated as inputs, x1:T [`] are treated as outputs and Q is the likelihood covariance matrix. This mixture of Gaussians is an expressive representation of the predictive density which can, for instance, correctly take into account multimodality arising from ambiguity in the measurements. Although factorized covariance matrices can be pre-computed, the overall computational cost will increase linearly with L.The computational cost can be reduced by thinning the Markov chain using e.g. random sub-sampling or kernel herding [19]. In some situations it could be useful to obtain an approximation from the mixture of Gaussians consisting in a single GP representation. This is the case in applications such as control or real time filtering where the cost of evaluating the mixture of Gaussians can be prohibitive. In those cases one could opt for a pragmatic approach and learn the mapping x∗ 7→ f ∗ from a cloud of points {x0:T [`], f0:T [`]}L `=1 using sparse GP regression. The latent function values f0:T [`] can be easily sampled from the normally distributed p(f0:T [`] | x0:T [`], θ[`]).

6 6.1

Experiments Learning a Nonlinear System Benchmark

Consider a system with dynamics given by xt+1 = axt + bxt /(1 + x2t ) + cut + vt , vt ∼ N (0, q) and observations given by yt = dx2t + et , et ∼ N (0, r), with parameters (a, b, c, d, q, r) = (0.5, 25, 8, 0.05, 10, 1) and a known input ut = cos(1.2(t + 1)). One of the difficulties of this system is that the smoothing density p(x0:T | y0:T ) is multimodal since no information about the sign of xt is available in the observations. The system is simulated for T = 200 time steps, using log-normal priors for the hyper-parameters, and the PGAS sampler is then run for 50 iterations using N = 20 particles. To illustrate the capability of the GP-SSM to make use of a parametric model as baseline, we use a mean function with the same parametric form as the true system, but parameters (a, b, c) = (0.3, 7.5, 0). This function, denoted model B, is manifestly different to the actual state transition (green vs. black surfaces in Figure 2), also demonstrating the flexibility of the GP-SSM. Figure 2 (left) shows the samples of x0:T (red). It is apparent that the distribution covers two alternative state trajectories at particular times (e.g. t = 10). In fact, it is always the case that this bi-modal distribution covers the two states of opposite signs that could have led to the same observation (cyan). In Figure 2 (right) we plot samples from the smoothing distribution, where each circle corresponds to (xt , ut , E[ft ]). Although the parametric model used in the mean function of the GP (green) is clearly not representative of the true dynamics (black), the samples from the smoothing distribution accurately portray the underlying system. The smoothness prior embodied by the GP allows for accurate sampling from the smoothing distribution even when the parametric model of the dynamics fails to capture important features. To measure the predictive capability of the learned transition dynamics, we generate a new dataset consisting of 10 000 time steps and present the RMSE between the predicted value of f (xt , ut ) and the actual one. We compare the results from GP-SSM with the predictions obtained from two parametric models (one with the true model structure and one linear model) and two known models (the ground truth model and model B). We also report results for the sparse GP-SSM using an FIC prior with 40 inducing points. Table 1 summarizes the results, averaged over 10 independent training and test datasets. We also report the RMSE from the joint smoothing sample to the ground truth trajectory. Table 1: RMSE to ground truth values over 10 independent runs. RMSE Ground truth model (known parameters) GP-SSM (proposed, model B mean function) Sparse GP-SSM (proposed, model B mean function) Model B (fixed parameters) Ground truth model, learned parameters Linear model, learned parameters

7

prediction of f ∗ |x∗t , u∗t , data – 1.7 ± 0.2 1.8 ± 0.2 7.1 ± 0.0 0.5 ± 0.2 5.5 ± 0.1

smoothing x0:T |data 2.7 ± 0.5 3.2 ± 0.5 2.7 ± 0.4 13.6 ± 1.1 3.0 ± 0.4 6.0 ± 0.5

20 15

20

Samples Ground truth 1/2 ±(max(yt,0)/d)

15

10 5

f(t)

10

State

0 −5

5

−10

0

−15

−5

−20 1

−10

0.5

−15

0

−20

−0.5

0

10

20

30

40

50

60

Time

−20

−1

u(t)

−15

−10

0

−5

5

10

15

20

x(t)

Figure 2: Left: Smoothing distribution. Right: State transition function (black: actual transition function, green: mean function (model B) and red: smoothing samples). 16

10 2

14

x 12 10 8

θ˙

x˙ 0

2

5

1

0

θ 0

−5

−1

−2

−2

−10 300

350

300

350

300

350

300

350

Figure 3: One step ahead predictive distribution for each of the states of the cart and pole system. Black: ground truth. Colored band: one standard deviation from the mixture of Gaussians predictive. 6.2

Learning a Cart and Pole System

We apply our approach to learn a model of a cart and pole system used in reinforcement learning. The system consists of a cart, with a free-spinning pendulum, rolling on a horizontal track. An external force is applied to the cart. The system’s dynamics can be described by four states and a set of nonlinear ordinary differential equations [20]. We learn a GP-SSM based on 100 observations of the state corrupted with Gaussian noise. Although the training set only explores a small region of the 4-dimensional state space, we can learn a model of the dynamics which can produce one step ahead predictions such the ones in Figure 3. We obtain a predictive distribution in the form of a mixture of Gaussians from which we display the first and second moments. Crucially, the learned model reports different amounts of uncertainty in different regions of the state-space. For instance, note the narrower error-bars on some states between t = 320 and t = 350. This is due to the model being more confident in its predictions in areas that are closer to the training data.

7

Conclusions

We have shown an efficient way to perform fully Bayesian inference and learning in the GP-SSM. A key contribution is that our approach retains the full nonparametric expressivity of the model. This is made possible by marginalizing out the state transition function, which results in a nontrivial inference problem that we solve using a tailored PGAS sampler. A particular characteristic of our approach is that the latent states can be sampled from the smoothing distribution even when the state transition function is unknown. Assumptions about smoothness and parsimony of this function embodied by the GP prior suffice to obtain high-quality smoothing distributions. Once samples from the smoothing distribution are available, they can be used to describe a posterior over the state transition function. This contrasts with the conventional approach to inference in dynamical systems where smoothing is performed conditioned on a model of the state transition dynamics.

8

References [1] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning.

MIT Press, 2006.

[2] R. Turner, M. P. Deisenroth, and C. E. Rasmussen, “State-space inference and learning with Gaussian processes,” in 13th International Conference on Artificial Intelligence and Statistics, ser. W&CP, Y. W. Teh and M. Titterington, Eds., vol. 9, Chia Laguna, Sardinia, Italy, May 13–15 2010, pp. 868–875. [3] C. Andrieu, A. Doucet, and R. Holenstein, “Particle Markov chain Monte Carlo methods,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 72, no. 3, pp. 269–342, 2010. [4] F. Lindsten, M. Jordan, and T. B. Sch¨on, “Ancestor sampling for particle Gibbs,” in Advances in Neural Information Processing Systems 25, P. Bartlett, F. Pereira, C. Burges, L. Bottou, and K. Weinberger, Eds., 2012, pp. 2600–2608. [5] M. Deisenroth, R. Turner, M. Huber, U. Hanebeck, and C. Rasmussen, “Robust filtering and smoothing with Gaussian processes,” IEEE Transactions on Automatic Control, vol. 57, no. 7, pp. 1865 –1871, july 2012. [6] M. Deisenroth and S. Mohamed, “Expectation Propagation in Gaussian process dynamical systems,” in Advances in Neural Information Processing Systems 25, P. Bartlett, F. Pereira, C. Burges, L. Bottou, and K. Weinberger, Eds., 2012, pp. 2618–2626. [7] Z. Ghahramani and S. Roweis, “Learning nonlinear dynamical systems using an EM algorithm,” in Advances in Neural Information Processing Systems 11, M. J. Kearns, S. A. Solla, and D. A. Cohn, Eds. MIT Press, 1999. [8] J. Wang, D. Fleet, and A. Hertzmann, “Gaussian process dynamical models,” in Advances in Neural Information Processing Systems 18, Y. Weiss, B. Sch¨olkopf, and J. Platt, Eds. Cambridge, MA: MIT Press, 2006, pp. 1441–1448. [9] J. S. Liu, Monte Carlo Strategies in Scientific Computing.

Springer, 2001.

[10] A. Doucet and A. Johansen, “A tutorial on particle filtering and smoothing: Fifteen years later,” in The Oxford Handbook of Nonlinear Filtering, D. Crisan and B. Rozovsky, Eds. Oxford University Press, 2011. [11] F. Gustafsson, “Particle filter theory and practice with positioning applications,” IEEE Aerospace and Electronic Systems Magazine, vol. 25, no. 7, pp. 53–82, 2010. [12] M. K. Pitt and N. Shephard, “Filtering via simulation: Auxiliary particle filters,” Journal of the American Statistical Association, vol. 94, no. 446, pp. 590–599, 1999. [13] F. Lindsten and T. B. Sch¨on, “Backward simulation methods for Monte Carlo statistical inference,” Foundations and Trends in Machine Learning, vol. 6, no. 1, pp. 1–143, 2013. [14] F. Lindsten and T. B. Sch¨on, “On the use of backward simulation in the particle Gibbs sampler,” in Proceedings of the 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Kyoto, Japan, Mar. 2012. [15] D. K. Agarwal and A. E. Gelfand, “Slice sampling for simulation based fitting of spatial data models,” Statistics and Computing, vol. 15, no. 1, pp. 61–69, 2005. [16] E. Snelson and Z. Ghahramani, “Sparse Gaussian processes using pseudo-inputs,” in Advances in Neural Information Processing Systems (NIPS), Y. Weiss, B. Sch¨olkopf, and J. Platt, Eds., Cambridge, MA, 2006, pp. 1257–1264. [17] J. Qui˜nonero-Candela and C. E. Rasmussen, “A unifying view of sparse approximate Gaussian process regression,” Journal of Machine Learning Research, vol. 6, pp. 1939–1959, 2005. [18] M. Seeger, C. Williams, and N. Lawrence, “Fast Forward Selection to Speed Up Sparse Gaussian Process Regression,” in Artificial Intelligence and Statistics 9, 2003. [19] Y. Chen, M. Welling, and A. Smola, “Super-samples from kernel herding,” in Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence (UAI 2010), P. Gr¨unwald and P. Spirtes, Eds. AUAI Press, 2010. [20] M. Deisenroth, “Efficient reinforcement learning using Gaussian processes,” Ph.D. dissertation, Karlsruher Institut f¨ur Technologie, 2010.

9