A statistical and computational theory for robust and sparse Kalman smoothing

Kalman smoothers reconstruct the state of a dynamical system starting from noisy output samples. While the classical estimator relies on quadratic penalization of process deviations and measurement errors, extensions that exploit Piecewise Linear Quadratic (PLQ) penalties have been recently proposed in the literature. These new formulations include smoothers robust with respect to outliers in the data, and smoothers that keep better track of fast system dynamics, e.g. jumps in the state values. In addition to L2, well known examples of PLQ penalties include the L1, Huber and Vapnik losses. In this paper, we use a dual representation for PLQ penalties to build a statistical modeling framework and a computational theory for Kalman smoothing. We develop a statistical framework by establishing conditions required to interpret PLQ penalties as negative logs of true probability densities. Then, we present a computational framework, based on interior-point methods, that solves the Kalman smoothing problem with PLQ penalties and maintains the linear complexity in the size of the time series, just as in the L2 case. The framework presented extends the computational efficiency of the Mayne-Fraser and Rauch-Tung-Striebel algorithms to a much broader non-smooth setting, and includes many known robust and sparse smoothers as special cases.

💡 Research Summary

The paper presents a comprehensive theory for robust and sparse Kalman smoothing by extending the classical quadratic‑penalty framework to a broad class of Piecewise Linear‑Quadratic (PLQ) penalties. The authors begin by defining PLQ functions through a polyhedral set U and a positive‑semidefinite matrix M, yielding a convex function θ_{U,M}(w)=sup_{u∈U}{⟨u,w⟩−½⟨u,Mu⟩}. By allowing an affine transformation b+By, they obtain a family of functions ρ(y)=θ_{U,M}(b+By) that encompass many familiar loss functions such as the L₂ (quadratic), L₁ (absolute), Huber, and Vapnik ε‑insensitive losses.

A central contribution is the statistical interpretation of these PLQ penalties as negative log‑densities of proper probability distributions. The authors prove two key theorems: (1) coercivity of ρ (i.e., ρ(y)→∞ as ‖y‖→∞) guarantees integrability of exp(−ρ), and (2) coercivity is equivalent to the condition Bᵀ cone(U)°={0}. This condition is easy to verify for the canonical penalties, establishing that each corresponds to a legitimate density on an affine subspace.

With this statistical foundation, the paper formulates the Kalman smoothing problem for a linear state‑space model x_k=G_k x_{k−1}+w_k, z_k=H_k x_k+v_k where the process noise w_k and measurement noise v_k follow PLQ densities. The MAP estimator becomes the minimization of a sum of PLQ penalties applied to the process and measurement residuals, expressed compactly in equation (3.3). Because the PLQ sets U_w and U_v are polyhedral, the problem is an Extended Linear‑Quadratic Program (ELQP).

Instead of solving the ELQP directly, the authors derive the Karush‑Kuhn‑Tucker (KKT) conditions (equation (3.4)), which retain a block‑tridiagonal structure thanks to the bidiagonal nature of the state transition matrix G and the block‑diagonal nature of H, Q, R, and the transformation matrices B_w, B_v. This structure is crucial for computational efficiency.

The algorithmic core is an interior‑point (IP) method applied to the relaxed KKT system. Complementarity conditions (s·q=0) are replaced by μ‑scaled equations (s q=μ 1), and a damped Newton iteration solves the resulting nonlinear system. Because each Newton step involves solving a linear system with the same block‑tridiagonal sparsity, it can be performed in O(N n³+N m) time, where N is the number of time steps, n the state dimension, and m the measurement dimension. The number of IP iterations is typically bounded (10–20), so the overall complexity scales linearly with the horizon length, matching the classic Mayne‑Fraser and Rauch‑Tung‑Striebel smoothers.

The authors verify that the four canonical penalties satisfy the non‑degeneracy condition Null(M)∩U^∞={0}, ensuring that the IP algorithm applies without modification. Consequently, the proposed framework unifies robust (outlier‑resistant) smoothing via Huber or Vapnik losses and sparse smoothing via L₁ or L₀‑type approximations under a single computational scheme.

In the concluding section, the paper emphasizes that the PLQ‑based Kalman smoother offers a statistically sound, computationally efficient, and highly flexible tool for a wide range of applications where Gaussian assumptions are violated—such as fault‑tolerant navigation, financial time‑series with jumps, and compressed‑sensing‑style state estimation. By preserving linear‑time complexity while handling non‑smooth penalties, the work bridges a critical gap between robust statistical modeling and real‑time signal processing.