2.2 Dynamical models

Now for concreteness, let us consider in turn each of the four half-updates, referring to the corresponding graphical model, as they apply to the HMM. Since the filter distribution is categorical, the joint distribution represented by Fig. 2.5A is a(nother) mixture model. Marginalizing out the class identities in this model—the time update—therefore amounts to an application of Eq. 2.10. Similarly, the joint distribution of the model in Fig. 2.5B is another mixture model. The measurement update amounts to applying Bayes’ theorem to this mixture model—Eq. 2.6.

The smoother follows the same template. The joint distributions in both Fig. 2.5A (as we have just seen) and Fig. 2.5C are mixture models, and the “future conditioning” and “backward step” correspond, respectively, to “Bayes inversion” (Eq. 2.6) and marginalization (Eq. 2.10) in these two models.

What if we are not actually interested in the filter distribution, that is, if all we want at the end of the day are the smoother distributions, ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{t}\middle|\bm{y}_{1},\ldots,\bm{y}_{{T}}}\right)}$ and ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{t-1},\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{t}\middle|\bm{y}_{1},\ldots,\bm{y}_{{T}}}\right)}$ (at all steps $t$)? Then it can be wasteful and unnecessary to keep renormalizing in the measurement update, Eq. 2.37. We shall see below that it is not wasteful for Gaussian random variables, for which the posterior (cumulants) can be calculated without any explicit calculation of the normalizer. But for discrete random variables, we are nearly doubling the number of operations performed. On the other hand, normalization has agreeable numerical consequences, preventing the numerical underflow that would results from repeated multiplication of numbers less than one. Therefore, we present the following mostly to connect with the historical literature.

Recall again that the goal is to derive the “filtering” equation, ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{t}\middle|\bm{y}_{1},\ldots,\bm{y}_{t}}\right)}$. Here we translate the time update, Eq. 2.36, and the measurement update, Eq. 2.37, into operations on Gaussian random variables. We have lately noted that this required keeping track only of posterior mean vectors and covariance matrices, so to facilitate that process we begin by assigning shorthand symbols to the filter cumulants,

and likewise to the posterior mean and covariance given the observations only up to the previous step:

We hope that these are the only parameters we ever need to keep track of. To show that they are, we proceed by induction.

We expect the covariance, for its part, to shrink when we add new information. The more precisely the emission reflects the hidden state ($\mathbf{K}_{t}$), the larger the fraction of the existing covariance ($\mathbf{{\Upsilon}}_{t|t-1}$) we expect to remove.

We cannot quite say, however, that the uncertainty grows, even though we are not incorporating any new observations and the state transition is noisy: If the underlying dynamical system is strictly stable, i.e. all its eigenvalues are within the unit circle, then the covariance $\mathbf{{A}}\mathbf{{\Upsilon}}_{t-1|t-1}\mathbf{{A}}^{\text{T}}$ will be smaller than $\mathbf{{\Upsilon}}_{t-1|t-1}$; and it could in theory be sufficiently smaller even to account for the additional transition noise ($\mathbf{{\Sigma}}_{\hat{x}}$). Such systems, however, are heading rapidly toward an equilibrium point, and are (perhaps) less common targets for tracking than marginally stable systems, with eigenvalues on the unit circle. For such systems, the uncertainty increases with the time update.

To complete the induction, we need to show that the one-step prediction distribution at step $1$ is a normal distribution over ${\bm{\hat{X}}}_{1}$. And indeed, setting this distribution to be the prior distribution over ${\bm{\hat{X}}}_{1}$, ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{1}}\right)}$, achieves that goal. For compactness, then, we define

so that Eq. 2.47 gives the measurement update for all time steps. Now the time update, Eq. 2.48, and the measurement update, Eq. 2.47, completely define the filter.

Recall the goal: to get, for all $t$, ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{t}\middle|\bm{y}_{1},\ldots,\bm{y}_{{T}}}\right)}$, where we are conditioning on all the observations. Once again we shall show that this distribution is normal, and thus we need only keep track of its first two cumulants:

In deriving our algorithm, it will also be helpful to keep track of the “future-conditioned smoother,” ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{t-1}\middle|\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{t},\bm{y}_{1},\ldots,\bm{y}_{{T}}}\right)}$, so we give names to its cumulants, as well:

Notice that the “future-conditioned” mean is a random variable, because it depends on the random variable ${\bm{\hat{X}}}_{t-1}$. (The future-conditioned covariance, although it also seems to depend on a random variable, will turn out to be non-random.) This raises a minor question of how to represent and store this mean, an issue which did not arise for the HMM. We shall revisit this after the derivation of the smoother.

Since the Kalman filter leaves us with ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{{{T}}}\middle|\bm{y}_{1},\ldots,\bm{y}_{{{T}}}}\right)}$ at the last step, we start there and recurse backwards. This filter distribution is normal, which takes care of the base case of the induction. The induction step, then, assumes that at step $t$ we have ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{t}\middle|\bm{y}_{1},\ldots,\bm{y}_{{T}}}\right)}$, and concludes by showing that we can get ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}_{t-1}\middle|\bm{y}_{1},\ldots,\bm{y}_{{T}}}\right)}$ from it.

Let us try to interpret these equations. Notice that they relate the means and the covariances of four different posterior distributions. Essentially, we are updating the smoother (stepping backwards in time) by updating a smoother-filter disparity:

Both disparities measure how much an estimate changes when we do without the future observations, from $t$ to ${{T}}$. The less these observations matter in estimating ${\bm{\hat{X}}}_{t}$ (i.e., the smaller the right-hand disparity), the less they matter in estimating ${\bm{\hat{X}}}_{t-1}$ (the smaller the left-hand disparity).

We emphasize that both filter estimates, $\bm{\bar{\hat{x}}}_{t-1|t-1}$ and $\bm{\bar{\hat{x}}}_{t|t-1}$, are made on the basis of observations only up to time $t-$1; that is, they differ only by a time update, i.e. by taking into account the state transition. Therefore, if the state transition is particularly noisy, relative to the filter’s current precision, we shouldn’t take $\bm{\bar{\hat{x}}}_{t|t-1}$ so seriously, and accordingly the right-hand disparity should not propagate into much of a left-hand disparity. The “gain,” $\mathbf{J}_{t-1}$, should be small. On the other hand, if the state transition is particularly precise, relative to the filter’s current precision, then the right-hand disparity probably reflects a problem with the filter (at this point in time), and accordingly should propagate into the left-side disparity. The gain should be large.

Apparently, then, the gain $\mathbf{J}_{t-1}$ should encode the precision of the state-transition distribution, relative to the filter distribution at time $t-$1. From its definition in Eq. 2.53 and the original definition of the gain, Eq. 2.15, we see that this is exactly what $\mathbf{J}_{t-1}$ does. It should not be suprising, then, that this gain does not depend on the future observations—indeed, it can be computed in the forward pass (Eq. 2.53)!—even though it is used to update the smoother. The gain doesn’t relate the filter and smoother cumulants; it relates their disparities at two adjacent time steps. And these disparities are related to each other only insofar as time-updating the filter doesn’t much damage its precision.

Similarly, taking into account the future observations shifts, or rather shrinks, the standard deviations of the posterior covariances of ${\bm{\hat{X}}}_{t-1}$ and ${\bm{\hat{X}}}_{t}$ by a proportional amount.

Finally, to completely characterize the posterior distribution, we also require the cross covariance:

Once more, we note that the relevant marginal distributions are normal, so their joint is as well, and we can apply Eq. 2.19:

Let us consider first a detail of the implementation of the smoother. We noted above that the “future-conditioned” mean, ${\bm{\bar{\hat{X}}}}_{t-1|t,{{T}}}$, is a random variable. How do we represent it? We could in theory store the numerical values of the parameters that effect the affine transformations, from Eq. 2.53 (along with, at least implicitly, the transformation itself). But in fact most of the future-conditioning step has been absorbed into the backward step, Eq. 2.54, and in practice the former amounts merely to computation of the gain, $\mathbf{J}_{t-1}$. Breaking the backward step into two pieces was for us merely a conceptual step, a way to exhibit the parallelism of the filter and smoother.

Next, consider the case where the observation space is (much) higher-dimensional than the latent space. This occurs, for example, in neural decoders, where the firing rates of hundreds of neurons are observed, but are thought to represent only some low-dimensional dynamical system. When implementing filter and smoother algorithms for such data, we would have the matrix inversions occur in the latent space……