BP on Gaussian Hidden Markov Models: Kalman Filtering and Smoothing

Massachusetts Institute of Technology Department of Electrical Engineering and Computer Science 6.438 Algorithms for Inference Fall 2011 13 BP on Gaussian Hidden Markov Models: Kalman Filtering and Smoothing As we have seen, when our variables of interest are jointly Gaussian, the sum-product algorithm for inference on trees takes a special form in which the messages themselves are Gaussian. As we also saw earlier, a simple but very important class of trees are the Hidden Markov Models (HMMs). In this section of the notes we specialize our Gaussian BP algorithm to the class of Gaussian HMMs, which are equivalently characterized as linear dynamical systems (LDSs). As we will see, what results is a version of the forward-backward algorithm for such models that is referred to as Kalman filtering and smoothing. 13.1 Gaussian HMM (State Space Model) Consider the Gaussian HMM depicted in Fig. 1, where to simplify our development we restrict our attention to zero-mean variables and homogeneous models. In this model, the states xt and observations yt are d- and d' -dimensional vectors, respectively. The Gaussian HMM can be expressed as an LDS by exploiting the representation of the variables in innovation form.1 In particular, first, the states evolve according to the linear dynamics xt+1 = Axt + vt , vt ∼ N(0, Q), x0 ∼ N(0, Λ0 ), (1a) where A is a constant matrix, Λ0 and Q are covariance (i.e., positive semidefinite) matrices, and where {vt } is a sequence of independent random vectors that are also independent of x0 . Second, the observations depend on the state according to the linear measurements yt = Cxt + wt , wt ∼ N(0, R), (1b) where C is a constant matrix and R is a covariance matrix, and where {wt } is a sequence of independent random variables that is also independent of x0 and {vt }. Hence, in terms of our notation, we have the conditional distributions x0 ∼ N(0, Λ0 ) xt+1 |xt ∼ N(Axt , Q) yt |xt ∼ N(Cxt , R). 1 Recall that any pair of jointly Gaussian random vectors u, v can be expressed in innovation form: there exists a matrix G and a Gaussian random vector w (the innovation) that is independent of v such that u = Gv + w. x0 x1 x2 xt y0 y1 y2 yt Figure 1: A Gaussian HMM as a directed tree. In this model, the variables are jointly Gaussian. Such Gaussian HMMs are used to model an extraordinarily broad range of phe­ nomena in practice. 13.2 Inference with Gaussian HMMs In our development of Gaussian BP for trees, we expressed the algorithm in terms of the information form representation of the Gaussian variables involved. In this section, we further specialize our results to the HMM, and express the quantities involved directly in terms of the parameters of the LDS (1). To obtain our results, we start by observing that the joint distribution factors as p(xt0 , y0t ) = p(x0 ) p(y0 |x0 ) p(x1 |x0 ) p(y1 |x1 ) · · · p(yt |xt ), where we have used xt0 as a shorthand to denote the entire sequence x0 , . . . , xt . Substituting the Gaussian form of the constituent conditional distributions and matching corresponding terms, we then obtain � � � � 1 T 1 t t T −1 p(x0 , y0 ) ∝ exp − x0 Λ0 x0 exp − (y0 − Cx0 ) R (y0 − Cx0 ) 2 2 � � 1 T −1 · exp − (x1 − Ax0 ) Q (x1 − Ax0 ) · · · 2 � � 1 T 1 T −1 1 T T −1 T T −1 ∝ exp − x0 Λ0 x0 − y0 R y0 − x0 C R Cx0 + x0 C R y0 2 2 2 � � 1 T −1 1 T T −1 T −1 · exp − x1 Q x1 − x0 A Q Ax0 + x1 Q Ax0 + · · · 2 2 � � 1 T 1 T −1 1 T −1 T −1 T −1 ∝ exp − x0 (Λ0 + C R C + A Q A)x0 − y0 R y0 − x1 Q x1 2 2 2 2 = t t ( ) T −1 T −1 · exp xT 0 C R y0 + x1 Q Ax0 + · · · φk (xk ) k=0 t t ψk−1,k (xk−1 , xk ) k=1 t t ϕk (yk ) k=0 t t ϕ̃k (xk , yk ), k=t where ⎧ 1 T −1 ⎪ − xT + AT Q−1 A)x0 ⎪ 0 (Λ0 + C R C ⎪ M ⎪ 2 ⎪ ⎪ ⎪ �J0 ⎪ ⎪ ⎨ 1 T −1 T −1 T −1 log φk (xk ) = − 2 xk (Q + C R MC + A Q A)xk ⎪ ⎪ �Jk ⎪ ⎪ ⎪ 1 ⎪ T −1 T −1 ⎪ − x (Q + C R C)xk ⎪ ⎪ M ⎩ 2 k k=0 1≤k ≤t−1 k=t �Jt −1 log ψk−1,k (xk−1 , xk ) = xT k Q M A xk−1 , �Lk 1 log ϕk (yk ) = − y0T R−1 y0 , 2 T −1 log ϕ̃k (xk , yk ) = xk CT R M yk . �Mk Next note that y0t are observations, so we condition on these variables. After conditioning on y0t , the joint distribution continues to be Gaussian, and simplifies to the following form pxt0 |y0t (x0t |y0t ) ∝ px0t ,y0t (xt0 , y0t ) ∝ t t k=0 t t φk (xk ) t t ϕ̃k (xk , yk ) k=0 t t ψk−1,k (xk−1 , xk ) k=1 t �t � 1 ( ) T T exp −xT ∝ exp − xk Jk xk + xk Mk yk (−Lk )xk−1 . k M 2 k=0 k=1 �hk Thus the potential for node k is a Gaussian with information parameters (hk , Jk ), and the potential for the edge from node k − 1 to node k is a Gaussian with information parameters (0, Lk ); see Fig. 2. Having now expressed the joint distribution as a product of potentials in Gaussian information form, we can easily read off the information matrices and potential vectors required to apply Gaussian BP. In particular, the initial messages are J0→1 = −L1 J0−1 LT 1 −1 h0→1 = L1 J0 h0 . 3 exp xT 1 ( L1 )x0 x0 exp x1 1 T 2 x0 J0 x0 x2 xt + hT 0 x0 Figure 2: The equivalent undirected tree for Gaussian HMMs, with only the first node and edge potentials shown. Then, the “forward” pass is given by the recursion Ji→i+1 = −Li+1 (Ji + Ji−1→i )−1 LT i+1 hi→i+1 = Li+1 (Ji + Ji−1→i )−1 (hi + hi−1→i ) , i = 0, 1, . . . , t − 1. In turn, the “backward” pass is given by the recursion Ji+1→i = −Li+1 (Ji+1 + Ji+2→i+1 )−1 LT i+1 , hi+1→i = Li+1 (Ji+1 + Ji+2→i+1 )−1 (hi+1 + hi+2→i+1 ). i = t − 1, t − 2, . . . , 0. Finally, the parameters for the marginals for each i are obtained via Ĵi = Ji + Ji−1→i + Ji+1→i ĥi = hi + hi−1→i + hi+1→i , from which the mean vectors and covariance matricses associated with the t marginals are given by ˆ and Σi = Ĵ−1 . µi = Ĵ−1 i hi i As a reminder, the complexity of this algorithm scales as O(td3 ) compared to O((td)3 ) for a naive implementation of the marginalization. The savings is particularly signif­ icant because t » d in typical applications. 13.3 Kalman Filtering and Smoothing The implementation of BP for Gaussian HMMs developed in the previous section has a convenient two-pass structure, consisting of a forward and a backward pass. Note that the data yk enters into the computation through the potential vector hk , and thus in the serial version of BP each piece of the data is used three times: once 4 during the forward pass, once during the backward pass, and once during the marginal computation. There are a variety of ways to rearrange such forward-backward computation for efficiently producing marginals. One such variant corresponds to what is referred to as the Rauch-Tung-Striebel (RTS) algorithm. In this algorithm, the forward and back­ ward passes take the form of what are referred to as Kalman filtering and smoothing, respectively. The Kalman filter was introduced by R. Kalman in his extraordinarily influential work in 1959.2 In contrast to the forward pass of the Gaussian BP procedure de­ veloped in the previous section, the Kalman filter directly generates a particular set of marginals. Specifically, it generates the marginals p(xi |y0i ) for i = 0, 1, . . . , t, i.e., a marginal at node i based on the data only through index i; this is what is meant by the term “filtering.” Each step in the forward pass is typically implemented in two substeps: a prediction substep that generates p(xi |y0i−1 ), followed by an update substep that generates p(xi |y0i ). The backward pass of the RTS algorithm, which implements Kalman smoothing, directly generates the desired marginals p(xi |y0t ), i.e., a marginal at node i based on all the data. Moreover, it does so in a manner that requires only access to the marginals obtained in the forward pass in the computation—once the forward pass is complete, the data is no longer needed. From the above perspectives, the RTS algorithm can be viewed as the Gaussian version of the so-called α, γ variant of the forward-backward algorithm for discrete HMMs. Indeed, suitably normalized, the α messages are the “filtered” marginals, and the γ messages are the “smoothed” marginals. The α, γ forward-backward algorithm and its relation to the α, β forward-backward algorithm we introduced as BP for discrete HMMs are developed in Chapter 12 of Jordan’s notes. There are several versions of the Kalman filter, both in covariance and information forms. As an illustration, we summarize a standard version of the former. As notation, we use µi|j and Σi|j to refer to the mean vector and covariance matrix, respectively, of the distribution p(xi |y0j ). We omit the derivation, which involves substituting the Gaussian forms into the α, γ version of the forward-backward algorithm, with summations replaced by integrals, and evaluating the integrals, the same manner that Gaussian BP was obtained from our original BP equations. More specifically, the forward (filtering) recursion is � α(xi+1 ) = α(xi ) p(xi+1 |xi ) p(yi+1 |xi+1 ) dxi , 2 You are strongly encouraged to read at least the first two pages to appreciate the broader context as well as significance of this work. 5 and the backward (smoothing) recursion is � � � Z α(xi ) p(xi+1 |xi ) dxi+1 , γ(xi ) = γ(xi+1 ) � α(x'i ) p(xi+1 |x'i ) dxi' where, for reference, recall that the γ marginals are expressed in terms of our original β messages via α(xi ) β(xi ) γ(xi ) = . p(y0t ) 13.3.1 Filtering The filtering pass produces the mean and covariance parameters µi|i and Σi|i , re­ spectively, in sequence for i = 0, 1, . . . , t. Each step consists of the following two substeps: Prediction Substep: In this substep, we predict the next state xi+1 from the current observations y0i , whose impact is summarized in the filtered marginals p(xi |y0i ) and thus does not require re-accessing the data: p(xi |y0i ) → p(xi+1 |y0i ), i = 0, 1, . . . , t − 1. In terms of the model parameters (1), the prediction recursion for the associated marginal parameters can be derived to be µi+1|i = Aµi|i Σi+1|i = AΣi|i AT + Q. i = 0, 1, . . . , t − 1, where the initialization of the recursion is µ0|−1 = 0 (2) Σ0|−1 = Λ0 . (3) Update Substep: In this substep, we update the prediction at step i by incorpo­ rating the new data yi+1 , i.e., p(xi+1 |y0i ) → p(xi+1 |y0i+1 ), i = 0, 1, . . . , t − 1. In terms of the model parameters (1), the update recursion for the associated marginal parameters can be derived to be µi+1|i+1 = µi+1|i + Gi+1 (yi+1 − Cµi+1|i ) Σi+1|i+1 = Σi+1|i − Gi+1 CΣi+1|i i = 0, 1, . . . , t − 1, where Gi+1 = Σi+1|i CT (CΣi+1|i CT + R)−1 is a precomputable quantity referred to as the Kalman gain. 6 13.3.2 Smoothing The filtering pass produces the mean and covariance parameters µi|t and Σi|t , respec­ tively, in sequence for i = t, t − 1, . . . , 0. Each step implements p(xi+1 |y0t ) → p(xi |y0t ), i = t, t − 1, . . . , 0. In terms of the model parameters (1), the smoothing recursion for the associated marginal parameters can be derived to be µi|t = µi|i + Fi (µi+1|t − µi+1|i ) Σi|t = Fi (Σi+1|t − Σi+1|i )FT i + Σi|i where Fi = Σi|i AT Σ−1 i+1|i . 7 i = t, t − 1, . . . , 0.