Finite-horizon square-root Kalman filter by rectangular factors and QR

Fix a finite horizon $T\ge 0$, positive state dimensions $d_k$ for $0\le k\le T$, and positive observation dimensions $p_k$ for $1\le k\le T$. Inputs are real finite matrices $F_k\in {\mathbb{R}}^{d_k\times d_{k-1}}$, $H_k\in {\mathbb{R}}^{p_k\times d_k}$, observation vectors $y_k\in {\mathbb{R}}^{p_k}$, and an initial mean $m_0\in {\mathbb{R}}^{d_0}$. For each covariance supply a rectangular real factor $L_0\in {\mathbb{R}}^{d_0\times r_0},L_{Q,k}\in {\mathbb{R}}^{d_k\times r_{Q,k}},L_{R,k}\in {\mathbb{R}}^{p_k\times r_{R,k}},$ where each column count is any finite nonnegative integer. The represented covariances are $P_0=L_0{L}_0^T$, $Q_k=L_{Q,k}{L}_{Q,k}^T$, and $R_k=L_{R,k}{L}_{R,k}^T$. Thus zero columns and rank-deficient factors are allowed. The admissible statistical inputs are those of Linear-Gaussian state-space models with independent noises for which the canonical innovation $S_k=H_k{P}_k^{-}{H}_k^T+R_k$ is positive definite at each step. The The Axiom of Choice assumption is inherited only for the interpretation as a posterior; the finite matrix procedure itself uses no probability-space choice.

Choose a unit roundoff $u$ in the no-underflow operationwise model of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. During an update, whenever an entry $a$ of a supplied matrix, factor, or vector is read as an arithmetic operand or copied into a QR prearray, first convert it to $\mathrm{fl}(a)$ in that format. Return range-failure at that read if conversion fails under the same no-underflow rule. Read operands from left to right; for products and block copies, visit output entries column by column, with increasing row index within each column. Repeat the conversion on each such read; the supplied arrays themselves remain unchanged. The formulas below use these converted entries during execution, while $P_0,Q_k,R_k,F_k,H_k,y_k,m_0$ denote the original exact target data in comparisons with the canonical recursion. In particular, conversion errors must be included in any certified local error bounds. For $T=0$, return the input pair immediately without reading or converting its entries. Every matrix or vector product below is evaluated as left-to-right dot products. A nonzero nonnormal intermediate, nonfinite result, zero innovation diagonal, or failed factorization returns a status and the first failing step; it is not coerced into a covariance estimate. Each arithmetic and factorization failure has a declared return; there is no iterative convergence loop or unbounded refinement search.

The QR operation used below is the following unpivoted reduced procedure on $A\in {\mathbb{R}}^{m\times n}$ with floating-point entries and $m\ge n\ge 1$. Initialize a working array $R^{(0)}=A$ and an empty ordered reflector list. For $j=1,\dots ,n$, take $x=R_{j:m,j}^{(j-1)}$ and $B=R_{j:m,j+1:n}^{(j-1)}$, call Scaled Householder reflector formation and application on $(x,B)$, and, on success, replace the active block by its returned $C_{+}$. If the call fails, return its first status and step. Store the exact orthogonal reflector represented by the returned record, or the identity for its identity branch. After $n$ successful steps, return the computed upper-triangular block $R=R_{1:n,1:n}^{(n)}$ and the ordered reflector list representing the exact product $Q_{\mathrm{f}\mathrm{u}\mathrm{l}\mathrm{l}}$ whose first $n$ columns are the reduced $Q$ factor. Do not materialize a rounded dense $Q$. The conventions are those of Full, reduced, and column-pivoted computational QR factorisations.

Initialize $m_{0\mathrm{\mid }0}=m_0$ and $L_{0\mathrm{\mid }0}=L_0$. At step $k$, perform the following operations in order.

1. Compute the rounded mean prediction ${\widehat{m}}_{k\mathrm{\mid }k-1}=\mathrm{fl}(F_k{\widehat{m}}_{k-1\mathrm{\mid }k-1})$. Form the rounded factor product ${\widehat{X}}_k=\mathrm{fl}(F_k{L}_{k-1\mathrm{\mid }k-1})$. Concatenate ${\widehat{G}}_k=[{\widehat{X}}_k{L}_{Q,k}]$ and append zero columns until its column count is $q_k=\max (d_k,r_{k-1\mathrm{\mid }k-1}+r_{Q,k})$. Apply the above reduced Householder QR procedure to ${\widehat{G}}_k^T\in {\mathbb{R}}^{q_k\times d_k}$. On success let $R_k^p$ be its reduced upper factor and set $L_{k\mathrm{\mid }k-1}=(R_k^p{)}^T$. No diagonal positivity or full-rank condition is imposed on this prediction factor.

2. Compute ${\widehat{v}}_k=\mathrm{fl}(y_k-\mathrm{fl}(H_k{\widehat{m}}_{k\mathrm{\mid }k-1}))$. Pad $L_{R,k}$ with zero columns to ${\bar{r}}_k=\max (p_k,r_{R,k})$ columns. Compute ${\widehat{Y}}_k=\mathrm{fl}(H_k{L}_{k\mathrm{\mid }k-1})$ and form the block prearray ${\widehat{Z}}_k=\left[\begin{array}{cc}{L}_{R,k}^{\mathrm{p}\mathrm{a}\mathrm{d}}& {\widehat{Y}}_k\\ 0& L_{k\mathrm{\mid }k-1}\end{array}\right]\mathrm{.}$ Apply the above reduced Householder QR procedure to ${\widehat{Z}}_k^T\in {\mathbb{R}}^{({\bar{r}}_k+d_k)\times (p_k+d_k)}$. Its row count is at least its column count: the prediction step gives $L_{k\mathrm{\mid }k-1}\in {\mathbb{R}}^{d_k\times d_k}$, and ${\bar{r}}_k\ge p_k$. On success transpose its reduced upper factor and partition the resulting lower factor as ${\widehat{L}}_k=\left[\begin{array}{cc}{\widehat{L}}_{11,k}& 0\\ {\widehat{L}}_{21,k}& {\widehat{L}}_{22,k}\end{array}\right],$ with block sizes $p_k$ and $d_k$. If any diagonal entry of ${\widehat{L}}_{11,k}$ is zero, return innovation-factorization-failure. Otherwise, for each $i=1,\dots ,p_k$, multiply column $i$ of the whole factor ${\widehat{L}}_k$ by the sign of its diagonal entry. These exact sign changes make ${\widehat{L}}_{11,k}$ have positive diagonal and preserve ${\widehat{L}}_k{\widehat{L}}_k^T$ and its lower block structure. Use the resulting factor blocks below. No invertibility or positive diagonal is required of ${\widehat{L}}_{22,k}$.

3. Solve ${\widehat{L}}_{11,k}{\widehat{z}}_k={\widehat{v}}_k$ by forward substitution in the order of Forward and backward substitution for triangular linear systems. Set ${\widehat{m}}_{k\mathrm{\mid }k}=\mathrm{fl}({\widehat{m}}_{k\mathrm{\mid }k-1}+{\widehat{L}}_{21,k}{\widehat{z}}_k),L_{k\mathrm{\mid }k}={\widehat{L}}_{22,k}\mathrm{.}$ Do not form a covariance downdate or an explicit matrix inverse. The factor-defined covariance $L_{k\mathrm{\mid }k}{L}_{k\mathrm{\mid }k}^T$ is positive semidefinite as an exact product of the returned entries; no claim is made that a separately rounded materialization of that product remains positive semidefinite.

4. If any arithmetic operation or QR call returns a failure status, return the first such failure and its step; if the leading innovation block has a zero diagonal, return innovation-factorization-failure. Otherwise the finite loop returns the predicted and updated mean/factor arrays with success. The success status records arithmetic and factorization completion only.

For $T=0$ the loop is empty and the output is $(m_0,L_0)$ with no update status. In exact arithmetic, the result implements the covariance and mean map of Exact linear-Gaussian Kalman filter. The prediction factor identity is supplied by QR prediction preserves the canonical covariance factor, and the measurement block identity by The block QR measurement update gives the canonical Kalman update. Under the inherited AC assumption its statistical interpretation is the one in The finite-horizon linear-Gaussian Kalman filter is exact. The nearby covariance supplied by an unstructured QR backward error may introduce cross-correlations between model blocks. Therefore no backward-stability claim within the original independent-noise state-space model, and no uniform-in-time accuracy or stability claim, is part of this procedure.