The Kalman Filter
Max Welling
California Institute of Technology 136-93 Pasadena, CA 91125 welling@vision.caltech.edu
Until now we have only looked at systems without any dynamics or structure over time. The models had spatial structure but were de ned at one moment in time, i.e. they were \static". In this lecture we will analyse a powerfull dynamic model, which could be characterized as factor analysis over time, although the number of observed features is not necessarily larger than the number of factors. The idea is again to compute only the mean and the covariance statistics, i.e. to characterize the probabilities by Gaussians. This has the advantage of being completely tractible (for strongly nonlinear systems the Gaussian assumption can no longer hold). The power of the Kalman Filter (KF), is that it operates on-line. This implies, that to compute the best estimate of the state and its uncertainty, we can update the previous estimates by the new measurement. This implies that we don't have to consider all the previous data again, to compute the optimal estimates; we only need to consider the estimates from the previous time step and the new measurement. So what are KFs usually used for, or what do they model? It is not hard to motivate the KF, because in practice it can be used for almost \everthing that moves". Popular applications include, navigation, guidance, radar tracking, sonar ranging, satellite orbit computation, stock prize prediction, etc. These applications can summarized as denoising, tracking and control. It is used in all sorts of elds, like engineering, seismology, bioengineering, econometrics etc. For instance, when the Eagle landed on the moon, it did so with a KF. Also, gyroscopes in airplanes use KFs. And the list goes on and on. The main idea is that we like to estimate a state of some sort (location and velocity of airplane) and its uncertainty. However, we do not directly observe these states. We only observe some measurements from an array of sensors, which are noisy. As an additional complication, the states evolve in time, also with its own noise or uncertainties. The question now becomes, how can we still optimally use our measurements to estimate the unobserved (hidden) states and their uncertainties. In the following we will rst describe the Kalman Filter and derive the KF equations. We assume that the parameters of the system are xed (known). Then we derive the Kalman Smoother equations, which allow us to use measurements forward in time to help predict the state at the current time better. Because these estimates are usually less noisy than the if we used measurements up till the current time only, we say that we smooth the state estimates. Next, we will show how to employ the Kalman Filter and smoother equations to e ciently estimate the parameters of the model from training data. This will then be, not surprisingly, another instance of the EM algorithm. In the appendix you will nd some usefull lemma's and matrix equalities, together with the derivation of the \lag-one smoother", which is needed for learning. Let me brie y describe an example. Consider a ship at sea which has lost its bearings. They need to estimate their current position using the positions of the stars. In the meantime the ships moves on on a wavy sea. The question becomes, how can we incorporate the measurement optimally, to estimate the ships location at sea. Let's assume that the model for the ship's location over time is given by, yt+1 = yt + c + wt (1)

1 Introduction

2 Introductory Example

^ Let's also assume we have some estimate yt of the location at some time t and some uncertainty t2 . How does this change as the ships sails for one second. Of course, it will drift in the direction of its velocity by a distance c. However, the uncertainty grows, due to the waves. This is expressed as, ^ ^ yt+1 = yt + c (3)
1

This contains a drift term (constant velocity c), and a noise term1 . A noisy measurement is described by, (2) xt = yt + vt

The constant drift term will not be used in the main body of these classnotes

1

The uncertainty is given by,

(4) If, we do not do any measurements, the uncertainty in the position will keep growing, until we have no clue anymore as to where we are. If we add the information of a measurement, the nal estimate is weighted average between the observed position and the previous estimate, t 2 = 2+ 2 +1 t w

^ yt0 +1 =

We observe that if we have in nite con dence in the measurement v ! 0, then the new location estimate is simply equal to the measurement. Also, if we have in nite con dence in the previous estimate, the measurement is ignored. For the new uncertainty we nd, This is also easy to interpret, since it says that if one of the uncertainties dissappears, the total uncertainty disappears, since the other estimate can simply be ignored in that case. Notice that the uncertainty always decreases or stays equal, by adding the measurement. The estimates for the location and uncertainty, incorporatating the measurement, can be rewritten as follows, ^ ^ ^ yt0 +1 = yt+1 + Kt+1(xt+1 yt+1) (7) 02 2 (8) Kt+1) t+1 t+1 = (1
2 2 0 2 = t+1 v t+1 2 2 t+1 + v

2 2 t+1 v 2+1 + v yt+1 + t+1 + v xt+1 2 2^ 2 t

(5)

(6)

Kt+1 =

From this we can see that the state estimate is corrected by the measurement error, multiplied by a gain factor (the Kalman gain). If the gain is zero, no attention is paid to the measurement, if its one, we simply use the measurement as our new state estimate. Similarly for the uncertainties, if the measurement is in nitely accurate, the gain is one, which implies that there is no uncertainty left. On the other hand, if the measurement is worthless, the gain is zero, which therefore does not decrease the overall uncertainty. The above one dimensional example will be generalized to higher dimensions in the following.

2 t+1 2+1 + v 2 t

(9)

3 The Model

Let us rst introduce the state of the KF yt at time t. The state is a vector of dimension d and remains unobserved. At every time t we alo have a k dimenional vector of observations xt , which depend on the state and some additive Gaussian noise. We will assume the following dynamical model for the KF: yt+1 = Ayt + wt (10) xt = Byt + vt (11) Note that the dynamics is governed by a Markov process, i.e. the state at yt+1 is independent of all other states, given yt . The evolution noise and the measurement noise are assumed white and Gaussian, i.e. distributed according to, w Gw 0; Q] (12) v Gv 0; R] (13) The noise vectors vt and wt are also assumed to be uncorrelated with the the state yt . From this we simply derive, E yt ; vk ] = 0 8 t; k (14) E yt; wk ] = 0 t k (15) E xt ; vk ] = 0 t k 1 (16) E xt; wk ] = 0 t k (17) (18) E vt; wk ] = 0 8 t; k E vt ; vk ] = 0 t 6= k; = R t = k (19) (20) E wt; wk ] = 0 t 6= k; = Q t = k 2

The above model could be considered as a factor analysis model over time, i.e. at every instant we have a FA model, where the factors now depend on the factors of a previous time step. The initial state y1 is distributed according to y1 G ; ]: (21) It is easy to generalize this model to include a `drift' and external inputs. The drift is a constant change expressed by adding to the dynamical equation. yt+1 = Ayt + + wt 0 = A0 yt + wt (22)
0 where we incorporated the constant again through the rede nitions A0 = A; ] and yt = yt ; 1]. For the external inputs, we let the evolution depend on a l dimensional vector of inputs, ut , as follows, yt+1 = Ayt + Cut + wt (23)

This model is used when we want to control the system. We could even let the parameters fA; ; C; g depend on time. However, in the rest of this chapter we will assume the simplest case, i.e. a linear evolution without drift or inputs, and a linear measurement equation, with white uncorrelated noise. Since the initial state is Gaussian and the evolution equation is linear, this implies that the state at later times will remain Gaussian. We want to be able to estimate the state and the covariance of the state at any time t, given a set of observations x = fx1 ; :::; x g. If is equal to the current time t we say that we lter the state. If is smaller than t we say that we predict the state and nally, if is larger the t we say that we smooth the state. The main probability of interest is, p(yt jx ); (24) since it conveys all the information about the state yt at time t, given all the observations up to time . Since we are assuming that this probability is Gaussian, we only need to calculate its mean and covariance, denoted by ^ yt = E ytjx ] (25) ~ ~ Pt = E yt yt jx ] (26) ^ ~ where we de ned yt = yt yt , i.e. the state prediction error. Notice that these quantities still depend on the random variables x and are therefore random variables themselves. We will now prove however that the covariance P does actually not depend on x . P may be considered as a parameter therefore in the ~ following. To proof the above claim we simply show that the correlation between the random variables yt and x vanishes. For normally distributed random variables this implies that they are independent. ~ ^ Lemma 3 The random variables yt = yt yt and x = fx1 ; :::; x g are independent. ^ E yt x ] E yt (x ) x ] = Z dyt dx p(yt ; x ) yt x dyt dx p(yt ; x
Z Z

4 General Properties

proof

) yt x dyt dx p(x ) p(yt jx ) yt x = 0 since p(yt ; x ) = p(x ) p(yt jx ) From this we derive the following corollary, ^ ^ ~ ~ Pt1 ;t2 = E yt1 yt2 ] = E yt1 yt2 ] E yt1 yt2 ]

Z

dx p(x

)x

Z

dyt p(yt jx

) yt

(27)

^t Another usefull result we will need is the fact that the predicted measurement error "t = xt Byt 1 is independent of the measurements xt 1 . The predicted measurement error is also called the innovation, since it represents that part of the new measurement xt that can not be predicted using knowledge of the xt 1 measurements, since they are independent. 3

^ Lemma 4 The random variables "t = xt Bytt 1 and xt 1 = fx1; :::; xt 1g are independent. Proof The proof is simple, and proceeds again by proving that they are uncorrelated, ^ E "t; xt 1] = E xtxt 1 ] BE ytt 1 xt 1 ] t 1 ^t = BE yt x ] + BE vt xt 1 ] BE yt 1 xt 1 ] ~t = BE yt 1 xt 1 ] + BE vt xt 1 ] = 0; ~t where we used the result that yt 1 is independent of xt 1 proved directly above and (16). Notice, that since the innovation "t 1 is a function of xt 1 , this also implies the following corollary, for = 1; :::; t 1 (28) ^ Before we proceed we would like to remark we have not well motivated yt as the preferred estimate of the state yt given data x . Other possible choices could be, the most likely state given the data, or the most likely sequence of states y1 ; :::; yt given the data. It turns out that for Gaussian distributed random ^ variables, these objectives are equivalent. On top of that, it turns out that the yt is also the minimal variance estimator, i.e. it minimizes Pt .

E "t" ] = 0

5 Kalman Filter Equations
We will now proceed to derive the Kalman lter equations (i.e. = t). First write, p(yt jxt ) = p(xt jyt ) p(yt jxt 1 ) ; p(xt jxt 1 )

(29) (30) (31) (32) (33)

where

p(yt jx

t

1) =

Z

dyt 1 p(yt jyt 1 ) p(yt 1 jxt 1 )

The denominator in (29) is an unimportant normalization factor. The remaining densities are given by, = Gx Byt ; R] p(yt jyt 1 ) = Gy Ayt 1 ; Q] p(xt jyt ) t t

The equations (29) and (30) have interesting interpretations as a reactive reenforcement due to an observation and a di usion equation for the Gaussian probability between observations. Equation (29) is basically applying Bayes law in the same way as we did for Bayesian learning, i.e. we are updating the probability distribution of an unknown random variable by including evidence. This has the e ect of making the distribution peakier, i.e. less uncertain. The second equation (30) evolves the hidden state from one instant to the next, without considering more evidence. This has the e ect of making the distribution less peakier, i.e. introducing more uncertainty. Together, these equations express p(yt jxt ) in terms of p(yt 1 jxt 1 ) and may be used recursively. It is now easy to verify that in case of a prediction only the equation (30) remains, i.e. Z p(yt jx ) = dyt 1 p(yt jyt 1 ) p(yt 1 jx ) t, i.e. smoothing, is more di cult and will be dealt with later. Let us return to the ^t ltering case. We like to calculate the mean yt 1 and covariance Pt 1 of the pdf p(yt jxt 1 ), expressed in t t 1 t 1 ^ terms of the mean yt 1 and covariance Pt 1 of the density p(yt 1 jxt 1 ). Notice that these two estimators determine the densities completely, since they are Gaussian. For the mean we nd, ^ ytt
1 =

E ytjxt 1 ] = AE yt 1 jxt 1 ] + E wt 1 jxt 1 ] ^t 1 = Ayt 1
4

(35)

where we have used (17) and the fact that w has zero mean (12). For the covariance we write we rst write, ~ ytt
1 = (yt

~t 1 and notice that wt 1 is independent of yt 1 since this is a function of yt 1 and xt 1 , while wt 1 is independent of both (using 15 and 17). Thus we can write,

^ ytt 1 ) = Ayt 1 + wt ~t 1 = Ayt 1 + wt

1 1

^ Aytt

1 1

Ptt

1 =

^t Next we wish to calculate yt and Ptt interms of the above calculated quantities, using equation (29). We write, ^ G yt 1 ; Ptt 1 ] Gx Byt ; R] : p(yt jxt ) = y t (37) p(x jxt 1 ) t t

~ y E ytt 1 (~tt 1 )T ] ~t 1 ~t 1 = E (Ayt 1 + wt 1 )(Ayt 1 + wt 1 )T ] T ~ t 1 yt 1 = AE yt 1 (~t 1 )T ]AT + E wt 1 wt 1 ] = APt 1 AT + Q t 1

(36)

We will now use Lemma 1, applying (74) to the second term in the numerator of (37), and then applying (75) to the result of that we nd, p(yt jxt ) = k2 (xt ) Gyt

t

^t ( Pt 1 ] 1 + BT R 1 B) 1 ( Pt 1 ] 1 yt 1 + BT R 1 xt ); ( Pt 1 ] 1 + BT R 1B) 1 ]; (38) t t t Because we know that the multplication of two Gaussians is again a Gaussian and moreover that (38) must be normalized with respect to yt , we deduce the the factor k2 (xt ) = 1. Finally we must use Lemma 2 to show that p(yt jxt ) is a Gaussian distribution with the following mean and covariance, ^ ^ ^ ytt = ytt 1 + Kt(xt Bytt 1 ) t 1 t Pt = (I KtB)Pt Kt = Ptt 1 BT (R + BPtt 1 BT )
1

(39) (40) (41)

^0 where Kt is called the Kalman gain. These equations are initialized by y1 = and P0 = . These equations, 1 together with (35) and (36) constitute the celebrated Kalman Filter equations and allow one to estimate the state of the system on-line, i.e. every new observation can be used recursively, given the information that was already received before. Notice that the gain factor Kt grows if the measurement covariance R becomes smaller, thus putting more weight on the measurement residual (di erence between predicted an actual measurement). Also if the noise covariance Pt 1 becomes smaller, less emphasis is put on the measurement t residual. It is also instructive to notice that the evolution of the state noise Pt (and therefore the Kalman t gain Kt ), is independent of the measurements and states and may be precomputed. From a numerical point of view, equation (40) is not preferable, due to the fact that it is a di erence of two postive de nite matrices, which is not guaranteed to result in a positive de nite matrix, and may lead to numerical instabilities. This however, is easily xed by noting that from (41) we can derive,

KtRKT = (I KtB)Ptt 1BT KT ; t which can be used to rewrite (40) as,

(42) (43)

Ptt = (I KtB)Ptt 1 (I Kt B)T + KtRKT ; t which is a sum of two postive de nite matrices!

6 Kalman Smoother Equations x Next we want to solve the smoothing problem. This implies that we are going to include later measurements ; > t, to improve our estimates of the states yt . The resultant estimates will be smoother (less noisy). 5

First concentrate on the mean E yt 1 jx ], for > t. We will now invoke the corollary of Lemma 1, ^ ^ identifying y = yt 1 , x = yt , y = yt 1 , x = yt and Gz ; ] = p(yt 1 ; yt jx ). Using these identi cations we nd, ^ ^ yt 1 = E yt 1 jx ] = E yt 1 jyt = yt ; x ]: (44) Next we write, p(yt 1 ; yt ; xt 1 ; xt ; :::; x ) p(yt 1 ; yt jx ) = p(x ) t 1 t 1 t 1 t 1 j = p(xt ; :::; x jyt 1 ; yt ; x ) p(ypt(yt )1 ; x ) p(yt 1 jx ) p(x ) x t 1 y ) (y j = p(xt ; :::; x jp(tx p:::;txytjx1 ) 1p)(yt 1 jx ) t t; = k1 (yt ; x ) p(yt jyt 1 ) p(yt 1 jxt 1 ) ^t 1 t 1 = k1 (yt ; x ) Gy Ayt 1 ; Q] Gy 1 yt 1 ; Pt 1 ] (45) In the same spirit as the derivation for the kalman lter (see derivation around 37), we will now use Lemma 1, applying (74) to the second term in (45), and then applying (75) to the result of that, p(yt 1 ; yt jx ) p(yt 1 jyt ; x ) = p(yt jx ) 1 ]; ^t 1 = k2 (yt ; x ) Gy 1 ( Pt 1 ] 1 + AT Q 1 A) 1 ( Pt 1 ] 1 yt 1 + AT Q 1 yt ); ( Pt 1 ] 1 + AT Q 1 A) (46) t 1 t 1 t 1 Notice the similarity between (38) and (46). Because we know that the multplication of two Gaussians is again a Gaussian and moreover that (46) must be normalized with respect to yt 1 , we deduce the the factor k4 (yt ; x ) = 1. Finally we invoke (44) and Lemma 2 to nd, ^ 1 ^ yt 1 = ytt 1 + Jt 1 (^t ytt 1 ) y ^ (47) t 1 T t 1 1 Jt 1 = Pt 1A Pt ] (48) ^ where we used (36) in the last line. This is initialized with y , computed from the Kalman Filter equations. For the covariance we rst observe, ^ ~ 1 ~ ^ 1 yt 1 + Jt 1 yt = ytt 1 + Jt 1 Aytt 1 ; (49) where we used (35) and (47). Multiplying both sides with their respective transpose from the right, taking ~ ^ expectations, and using the fact that yt 1 is independent of yt (since the latter is a function of x and t 1 t 1 ~ ^ Lemma 3) and similarly for yt 1 and yt , gives us, ^ ^ t ^ 1^ 1 Pt 1 + Jt 1E yt yt ]JT 1 = Ptt 1 + Jt 1AE ytt 1 ytt 1 ]AT JtT 1: (50) 1 Then we use, ^ ^ E yt yt ] = E yt yt] Pt (51) T = E (Ayt 1 + wt 1 )(Ayt 1 + wt 1 ) ] Pt (52) = AE yt 1 yt 1 ]AT + Q Pt (53) where (15) and (20) and the corollary following Lemma 3 was used. Analoguesly, ^ 1^ 1 (54) E ytt 1 ytt 1 ] = E yt 1 yt 1] Ptt 1 1 Putting these together and using (36), we nd Pt 1 = Ptt 1 + Jt 1(Pt Ptt 1)JT 1 (55) t 1 which is initialized by P , computed from the Kalman Filter equations. Equations (47), (55) and (48) are the so called Kalman smoother equations. If we want to estimate a state yt given data x with > t, then we rst apply the Kalman lter equations recursively until we have reached the state y . While moving ^t ^t forward we store the values for yt , yt 1 , Pt and Pt 1 , t = 1::: . Then we move backward by applying the t t smoother equations, until we have reached the state t we would like to estimate. Because we include more observations in the estimation of the state, the result will be less noisy as compared to the Kalman lter result, hence the name smoother. t t t

6

We will now proceed to estimate the parameters f ; ; B; R; A; Qg of the Kalman lter model using EM. We consider the states yt as hidden variables, while x are the observations. We assume we have observed N sequences of length . The joint probability of the complete data is given by, p(y ; x

7 Parameter Estimation for the Kalman Model
) = p(y1 )
Y
t

For EM, we are interested in the expectation of the joint pdf over the posterior density,
Q

=2

p(yt jyt 1 )

Y t =1

p(xt jyt )

(56)

= =

N XZ

N 1 X Z dy p(y jx ) (d + k) log(2 ) n 2 n=1 + log det + ( 1) log det Q + log det R + (y1 )T 1 (y1 )

n

=1

dy p(y

jxn ) log p(y ; xn )]

+ +

X t X t =2

(yt Ayt 1 )T Q 1(yt Ayt 1 ) (57)

Inspection of this objective function reveals that the only su cient statistics that need to be calculated in the E-step are, ^ (58) E ytjxn ] = yt;n t = 1; :::; ^ ^ E ytyt jxn ] = Pt + yt;nyt;n Mn t = 1; :::; (59) t n ^ ^ E ytyt 1 jxn ] = Pt;t 1 + yt;nyt 1;n Mt;t 1 t = 2; :::; (60) Fortunately, except for the last one, these are precisely the quantities that can be calculated through the Kalman Filter and Smoother recursions. The last quantity, which is called the lag-one covariance smoother is computed in the appendix C. It is given by the following recursion, Pt 1;t 2 = Ptt 1JT 2 + Jt 1 (Pt;t 1 APtt 1)JT 2 : (61) 1 t 1 t which is initialized by, P ; 1 = (I K B)AP 1: (62) 1 We now concentrate on the M-step, which maximizes (57) with respect to the parameters of the joint density only, i.e. the parameters present in the posterior are held xed. Taking derivatives with respect to and equating to zero gives,
X @ Q= @ n=1
N

=1

(xt;n Byt )T R 1 (xt;n Byt )]

1 (^ y1;n =1

)=0) (63) (^1;n )T + y
N X n

1 new = N

N X

For

this implies,
@ @

n

^ y1;n ))

1Q = N

1 2

new

1 =N

N X n

N 1 X(Mn y ^ 2 n=1 1 1;n

T

T

=1

Mn 1

new new

T

1 = P1 + N 7

=1

(^1;n y

new )(^ 1;n

y

new )

T

(64)

where we used that Pt is independent of xn . Taking derivatives with respect to A gives.
XX @ Q= (Q 1 Mn 1 t;t @A n=1 t=2
N

Q 1AMn 1) ) t
=1 t=2

Anew =
For Q we nd,
@ Q= @Q

"

N XX n

=1 t=2

M

#"

n t;t

1

N XX n

M

# 1

n t

1

(65)

1( 2

Qnew = N ( 1 1) (Mn Anew Mn 1;t); t t n=1 t=2 where we used (65), and Mn 1;t = (Mn 1 )T . For B we have, t t;t
N XX

N 1 X X(Mn AMn 1)N Q 2 Mt;t 1AT + AMn 1AT ) ) t t 1;t t n=1 t=2

(66)

XX @ ^ Q= (R 1 xt;n yt;n @B n=1 t=1

N

R 1BMn) ) t
=1 t=1

Bnew =

"

N XX n

And nally we have for R, N @ 1 1 X X(x x ^ T Q = NQ @R 2 2 n=1 t=1 t;n t;n xt;n yt;n B N XX ^ R = 1 (x x B y x ); new =1 t=1

^ xt;nyt;n

#"

N XX n

M

# 1

n t

(67)

^ Byt;nxt;n + BMnBT ) ) t (68)

N

n

where we used (67). Alternating E-steps and M-steps will thus converge to the maximum likelihood estimates of these parameters. Notice that the recursion equations ful l a double role. They may be used to e ciently compute the E-step in the learning problem, and, once the parameters are xed, to estimate the optimal state and state-covariance of the dynamical system, possibly on line.

=1 t=1

t;n t;n

new t;n t;n

8 Computation of Likelihood
L =
N X n

To monitor the total log-likelihood of the system we may calculate,
=1

log p(xn )] =
Z

N XX n

The mean and covariance of the Gaussians p(xt;n jxt 1;n ) can be computed as follows (omitting n for notational convenience), ^ x 1 = t t

=1 t=2

X log p(xt;n jxt 1 )] + log p(x1;n )]: n n N

=1

(69)

= = =

Z Z Z

dxt p(xt jxt 1 ) xt dxt dxt dyt
Z Z

dyt p(xt jyt ) p(yt jxt 1 ) xt dyt t ^ Gx Byt ; R] Gy ytt 1 ; Ptt 1 ] xt t t

^t = Byt 1 :

^ Gy ytt 1 ; Ptt 1 ] Byt (70) 8

For p(x1 ) the above calculation gives, Similarly, for the covariance we nd,

^1 x0 = B ^ ^ xtt 1 xtt 1) t (71)

Htt

1 =

Z

= =

Z Z

dxt p(xt jxt 1 ) (xt xt dxt dyt
Z

dyt t ^ ^ ^ Gx Byt ; R] Gy ytt 1 ; Ptt 1 ] (xt xt xtt 1 xtt 1 ) t The covariance for p(x1 ) is then given by,

^ ^ ^ Gy ytt 1 ; Ptt 1 ] (R + Byt yt BT xtt 1 xtt 1 ) ^t ^t ^t ^t = R + B(Pt 1 + yt 1 yt 1 )BT Byt 1 yt 1 BT t t 1 T = R + BPt B :

(72) (73)

H0 = R + B BT 1

A Lemma's identities, Lemma 1 Let Gy

;

] denote a normal density with mean

and covariance . We have the following (74) (75)

Gx Ay; ] = k1 (x) Gy (AT 1 A) 1 AT 1 x; (AT 1 A) 1 ]; Gy a; A] Gy b; B] = Gy (A 1 + B 1 ) 1 (A 1 a + B 1 b); (A 1 + B 1 ) 1 ] Ga b; A + B]:
Also, if we write z = y; x], = y; x

] and y; y

=

h

y y jx

Z

yy xy

yx xx

i

, then we have, (76)

Gz ; ] Gx x ; xx] = Gy

dx Gz ;

] = Gy

yy yx

] ( xx ) 1(

x

x);

yy

yx

(

xx

) 1

xy

]

(77)

As a corollary we notice that

Lemma 2 Consider a d aT Pa > 08 a (i.e. positive eigenvalues). The following equalities hold, (P 1 + BT R 1 B) 1 = P PBT (BPBT + R) 1 BP (P 1 + BT R 1 B) 1 BT R 1 = PBT (BPBT + R) 1

E y jx = x ] = E y ] d matrix P > 0, a k k matrix R > 0 and a k

(78) d matrix B, where P > 0 implies (79) (80)

B Matrix Identities

In the derivations to follow, the following idetities are useful, aT Ab = tr AbaT ] tr AB] = tr BA] @ tr AB] = BT @A
@ tr AT B] @A

(81) (82) (83) (84) (85) (86)

= B log det A] = log det A 1 ] @ log det A] = (AT ) 1 @A 9

In the E-step of the EM algorithm that estimated the parameters of the KF, we needed the lter and smoother equations. However, one quantity remains undetermined, which is the so called lag-one-covariance-smoother, Pt;t 1. This quantity is only needed backwards. First we will derive the initial value the backward recursions,

C Lag-One Covariance Smoother
Ptt;t
1 =

~t 1 ~t where we used (11), (39) and (47) with = t. Next, we write this out and use that both yt 1 and yt 1 are independent of vt , which is proved using (14) and (16). This will give t (90) Pt;t 1 = Ptt;t 1 1 Ptt 1BT KT JT 1 KtBPtt;t 1 1 + Kt(BPtt 1 BT + R)KT JT 1 t t t t t 1 = (I Kt B)APt 1 ; (91) where in the last ~t independent of yt t line we used (41) and Pt 1 1 = APt 1 , which is proved using the fact that wt 1 is t;t 1 1 (15, 17). If we set t = we nd the initial condition, 1

~~ E ytt ytt 1 ] ~t ^t ~t = E fyt 1 Kt (xt Byt 1 )gfyt ~t ~t ~t = E fyt 1 Kt (Byt 1 + vt )gfyt

1 1 1 1

Jt Jt

^ 1 Kt (xt By 1 )g] t 1 ~ 1 Kt (Byt + vt )g] t t

(87) (88) (89)

P;

1 = (I

K B)AP

1: 1

(92)

The derivation for the backward recursion is somewhat elaborate. First we write from (47) and (35), ^ ~ 1 ~ ^ 1 yt 1 + Jt 1yt = ytt 1 + Jt 1 Aytt 1 (93) t 2 t 2 ^ ~ ~ ^ yt 2 + Jt 2 yt 1 = yt 2 + Jt 2 Ayt 2 (94) Next, we equate, ^ ^ y E (^t 1 + Jt 1 yt )(^t 2 + Jt 2yt 2 )T ] = y t 1 t 1 t 2 ^ ^ 2 E (^t 1 + Jt 1Ayt 1 )(^t 2 + Jt 2 Aytt 2 )T ] ) y y (95) (96) where lemma 3 was used several times to get rid of cross-terms. We will now rewrite some of the terms appearing above, ^ ^ E yt yt 1] = E ytyt 1 ] Pt;t 1 = E (Ayt 1 + wt 1)(Ayt 2 + wt 2)T ] Pt;t 1 = AE yt 1yt 2 ]AT + AE yt 1wt 2 ] Pt;t 1 = AE yt 1yt 2 ]AT + AE (Ayt 2 + wt 2)wt 2 ] Pt;t 1 = AE yt 1yt 2 ]AT + AQ Pt;t 1: (97) Next, ~ 1~ 2 ~ 2 ^ 2 ~ 2 E ytt 1 ytt 2 ] = E fytt 1 Kt 1(xt 1 Bytt 1 )gytt 2 ] = ~ 2 Ptt 2;t 2 Kt 1E (Bytt 1 + vt 1)~tt 2 ] = y 2 1 Ptt 2;t 2 Kt 1BPtt 2;t 2 ; (98) 1 1 where (39) and (11) was used. Next, ^ 1~ 2 E ytt 1 ytt 2 ] = ^ 2 ~ 2 ~ 2 E fytt 1 + Kt 1(Bytt 1 + vt 1)gytt 2 ] = Kt 1BPtt 2;t 2 : (99) 1 10 ^ ^ Pt 1;t 2 + Jt 1E yt yt 1 ]JT 2 = t t 1 t 2 t 1 t 2 ~ ~ ^ ~ ^ 1^ 2 E yt 1 yt 2 ] + Jt 1 AE yt 1 yt 2 ] + Jt 1 AE ytt 1 ytt 2 ]AT JT 2 ; t

Finally we have, ^ 1^ 2 E ytt 1 ytt 2 ] = E (^tt 1 + Kt 1"t 1)^tt 2 ] = y 2 y 2 ^ 2^ 2 E ytt 1 ytt 2 ] = E yt 1yt 2 ] Ptt 2;t 2; 1 where the most important ingredient was lemma 4. Putting this together, we have Pt 1;t 2 = Ptt 2;t 2 Kt 1BPtt 2;t 2 + 1 1

(100)

Jt 1 AKt 1BPtt Jt 1 Pt;t 1JT 2 t

2 1;t 2

Jt 1APtt

2 AT JT t 2 1;t 2

Jt 1 AQJT 2 + t

(101)

The second line can be rewriten as follows,

Ptt 1JT 2 1 t

(I Kt 1 B)Pt t t (I Kt 1 B)Pt

2 1;t 2 = 2 Pt 2 ] 1 APt 2 = 1 t 1 t 2

(102)

where we used (40) and (48). The third line can be rewritten as,

Jt 1A(Kt 1 BPtt 2;t 2 Ptt 2;t 2 AT JT 2 QJT 2) = t t 1 1 t 2 t 2 T T Jt 1A(Kt 1 BPt 1 Pt 1;t 2A Q)Jt 2 = Jt 1A(Kt 1 BPtt 2 Ptt 2)JT 2 = 1 1 t Jt 1A(I Kt 1B)Ptt 2 JT 2 = 1 t Jt 1APtt 1JT 2 ; 1 t where again (40) and (48) was used, and the fact that,

(103) (104)

Ptt 2 = Ptt 1

2 AT + Q 1;t 2

which can be derived analoguesly to (36). Finally, putting this together we derive the lag-one covariance smoother, Pt 1;t 2 = Ptt 1JT 2 + Jt 1 (Pt;t 1 APtt 1)JT 2 : (105) 1 t 1 t

11