control//state estimation//Kalman filter//covariance propagation//square-root Kalman filter

A square-root Kalman filter is an implementation of the Kalman filter that never stores the covariance \(\ca{P}\) itself but a triangular factor \(S\) with \(\ca{P}=SS^{\mathsf T}\), propagating and updating that factor directly, and it exists to keep the filter numerically sound on hardware with limited precision: a 24-state navigation filter in float32 on a microcontroller, running for hours. Mathematically it computes the same estimate as the ordinary filter; numerically it survives where the ordinary one degrades.

A square-root Kalman filter is an implementation of the Kalman filter that never stores the covariance P\ca{P}P itself but a triangular factor SSS with P=SST\ca{P}=SS^{\mathsf T}P=SST, propagating and updating that factor directly, and it exists to keep the filter numerically sound on hardware with limited precision: a 24-state navigation filter in float32 on a microcontroller, running for hours. Mathematically it computes the same estimate as the ordinary filter; numerically it survives where the ordinary one degrades.

The problem it solves is one of arithmetic. The textbook update P=(I−KH)P−\ca{P}=(I-KH)\cb{P^-}P=(I−KH)P− subtracts two nearly equal matrices whenever a measurement is precise, and rounding errors pile up step after step (covariance propagation). After thousands of steps P\ca{P}P can stop being symmetric or acquire a negative variance on its diagonal, an uncertainty that is impossible, and from then on the gain is nonsense and the filter diverges. Because a factor SSS has the square root of P\ca{P}P's condition number, working with it roughly doubles the useful precision: what needs double precision as P\ca{P}P fits in single precision as SSS. And SSTSS^{\mathsf T}SST is symmetric and non-negative by construction, whatever the rounding does to SSS.

There is a ladder of fixes, from a line of code to a different filter.

Forcing symmetry after each step, P←(P+PT)/2\ca{P}\leftarrow(\ca{P}+\ca{P}^{\mathsf T})/2P←(P+PT)/2, and a floor under the diagonal are patches. The Joseph form writes the update as a sum of two terms that are symmetric and positive by construction, at the cost of a few more multiplications, and stays correct even with a suboptimal gain:

P=(I−KH) P− (I−KH)T+KRKT.\ca{P}=(I-KH)\,\cb{P^-}\,(I-KH)^{\mathsf T}+K\cc{R}K^{\mathsf T}.P=(I−KH)P−(I−KH)T+KRKT.

The square-root filter is the top rung, for long runs on low-precision hardware.

The factor is kept with standard tools of matrix factorization: a Cholesky decomposition to start, and QR or rank-one updates at each step. The UD filter stores P=UDUT\ca{P}=UDU^{\mathsf T}P=UDUT with UUU unit triangular and DDD diagonal, which avoids square roots altogether and is a long-standing favourite of embedded navigation.

It is old and proven. A square-root formulation (Potter's) was used in the Apollo navigation software, whose computer had nothing like modern floating point, and the unscented Kalman filter has square-root versions for the same reason, since its sigma points need a factor of P\ca{P}P anyway.

Its costs are code complexity and debugging effort, more than computation: the step count is similar, but the equations are harder to read and to check against a reference implementation. On a desktop in double precision with a well-scaled state, the Joseph form or plain symmetrization is usually enough (variable scaling).