9.4 Linear-Gaussian state-space models

All three summands are Gaussians, so we only consider in detail the differentiation of one, the first. Optimizing first with respect to $\mathbf{{A}}$, we find again the familiar normal equations:

where in the last line we intepret the average to be under samples from within, as well as across, sequences.^††margin: implementation note: Since each sequence contributes only ${{T}}-1$ samples, one must remember to subtract $N_{\text{sequences}}$ from the total number of available samples before normalizing. Turning to the covariance matrix,

where again in the penultimate line we changed the interpretation of the average to be within as well as across sequences. The final simplification was carried out with the equation for $\mathbf{{A}}$, exactly the same as with factor analysis (see above).

Derivations precisely analogous lead to formalae for the other cumulants. Here are all of the equations:

In Eq. 9.3, as in the preceding exampels, the bracketed quantities (including the brackets) are the sufficient statistics. They are all averages of either vectors or outer products of vectors, reflecting the quadratic structure inherent in normal distributions. When the states are observed, all these quantities are computed as sample averages under the data distribution. In the context of the EM algorithm, where the states are unobserved, the relevant averaging distributions $\check{p}$ are the RTS smoothing distributions, multiplied by the data distribution:

Once again, computing the sufficient statistics belongs to the E step. So in this case, the E step requires first running the RTS smoother—more specifically, the Kalman filter followed by the RTS smoother—and then computing the sufficient statistics under them. To make this extremely concrete, we note that having the smoother distribution in hand means having a mean (vector) and covariance matrix at every time step, since the distribution is normal. To compute expected outer products, then, we have to combine these together.

For $\mathbf{{A}}$ and $\mathbf{{\Sigma}}_{\hat{x}}$ we need

Notice that all sums start at the second index, and (consequently) are normalized by ${{T}}-1$. That is why the first and third statistics are not identical. If the sequence is long enough, including all ${{T}}$ terms in the first and last statistics will probably have little effect, but both averages can be computed with very little computational cost.

For $\mathbf{{C}}$ and $\mathbf{{\Sigma}}_{\hat{y}|\hat{x}}$, in contrast, we collect statistics for all time. Thus, even though we also need a statistic we have called ${\mathopen{}\mathclose{{}\left\langle{{\bm{\hat{X}}}_{t}{{}{\bm{\hat{X}}}_{t}}^{\text{T}}}}\right\rangle_{{\bm{\hat{X}}}}}$, in this case the average is over all $t$:

We also need average outer products involving the observations, but these can be computed directly rather than via the posterior expectation and covariance:

Finally, the sufficient statistics for the initial cumulants, $\bm{\mu}_{1}$ and $\mathbf{{\Sigma}}_{1}$, require no sums at all, since they rely only on the first time step:

We emphasize here again that, just as in the HMM, the posterior cumulants $\mathbb{E}_{{\bm{\check{X}}}_{1}{}|{\bm{Y}}}{\mathopen{}\mathclose{{}\left[{\bm{\check{X}}}_{1}\middle|\bm{y}_{1},\ldots,\bm{y}_{{T}}{}}\right]}$ and $\text{Cov}_{{\bm{\check{X}}}_{1}{}|{\bm{Y}}}{\mathopen{}\mathclose{{}\left[{\bm{\check{X}}}_{1}\middle|\bm{y}_{1},\ldots,\bm{y}_{{T}}{}}\right]}$ depend on the observations for all time: at least in theory, one must run the filter all the way to the end of the sequence and then the smoother all the way back before computing them.^††margin: implementation note: The second of these can yield a low-rank covariance matrix if only a single trajectory has been observed, and care is sometimes required to avoid inverting a singular matrix in the Kalman filter.