Exact linear-Gaussian Kalman filter

Fix a finite horizon $T\ge 0$, positive state dimensions $d_0,\dots ,d_T$, and positive observation dimensions $p_1,\dots ,p_T$. The input is a linear-Gaussian model from Linear-Gaussian state-space models with independent noises, with matrices $F_k\in {\mathbb{R}}^{d_k\times d_{k-1}}$ and $H_k\in {\mathbb{R}}^{p_k\times d_k}$, initial mean $m_0$, and symmetric positive-semidefinite covariances $P_0,Q_k,R_k$. Its state and observation equations are $X_k=F_k{X}_{k-1}+q_k,Y_k=H_k{X}_k+r_k,$ where $q_k\sim N(0,Q_k)$ and $r_k\sim N(0,R_k)$ are centered and mutually independent across all times and of $X_0$.

The admissible data are vectors $y_k\in {\mathbb{R}}^{p_k}$ for which every innovation covariance encountered below is positive definite. Positive definite means $u^T{S}_{k}u>0$ for every nonzero $u\in {\mathbb{R}}^{p_k}$. The algorithm returns predicted and updated means and covariances and the innovation and gain at each time. Initialize $m_{0\mathrm{\mid }0}=m_0,P_{0\mathrm{\mid }0}=P_0\mathrm{.}$ For $k=1,\dots ,T$, perform:

1. Predict: $m_{k\mathrm{\mid }k-1}=F_k{m}_{k-1\mathrm{\mid }k-1},P_{k\mathrm{\mid }k-1}=F_k{P}_{k-1\mathrm{\mid }k-1}{F}_k^T+Q_k\mathrm{.}$
2. Form the innovation and its covariance: $v_k=y_k-H_k{m}_{k\mathrm{\mid }k-1},S_k=H_k{P}_{k\mathrm{\mid }k-1}{H}_k^T+R_k\mathrm{.}$
3. Factor $S_k=L_k{L}_k^T$ with lower-triangular $L_k$ having positive diagonal, using A matrix admits a Cholesky factorisation with positive diagonal exactly when it is Hermitian positive definite, and that factor is unique. Solve $L_k{U}_k=H_k{P}_{k\mathrm{\mid }k-1}$ and $L_k^T{K}_k^T=U_k$ by triangular solves. This defines $K_k^T=S_k^{-1}{H}_k{P}_{k\mathrm{\mid }k-1}$, hence $K_k=P_{k\mathrm{\mid }k-1}{H}_k^T{S}_k^{-1}$, without explicitly forming $S_k^{-1}$.
4. Update: $m_{k\mathrm{\mid }k}=m_{k\mathrm{\mid }k-1}+K_k{v}_k,P_{k\mathrm{\mid }k}=P_{k\mathrm{\mid }k-1}-K_k{S}_k{K}_k^T\mathrm{.}$

For $T=0$ the loop is empty and the returned pair is $(m_0,P_0)$. No inverse of $P_0$, $Q_k$, or $R_k$ is used. This item defines the procedure and makes no complexity, rounding-error, or numerical-stability claim.