State-space RTS smoother is covariance-form despite 'square-root' docstring #668

Category: literature · Severity: medium
Location: gpjax/state_space/inference.py:315 — rts_smoother``

Problem

The docstring claims a 'square-root RTS smoother' but the backward step is covariance-form: it materialises P = L·Lᵀ (inference.py:367) and re-roots via eigh (_psd_sqrt, inference.py:395) each step. This loses the float32 conditioning benefit of square-root filtering and pays an O(d³) eigh per step.

Recommendation

Implement a QR pre-array square-root smoother (Park & Kailath 1996; Särkkä & Solin 2019 §10.7), or fix the docstring to say covariance-form. Prefer the square-root rework for numerical robustness.

_{Filed from the 2026-07-02 GPJax codebase audit (medium tier). Companion spec for the high-severity items: plans/2026-07-02-audit-high-severity-remediation-design.md.}