Numerically Robust Fixed-Point Smoothing Without State Augmentation
Authors: Nicholas Krämer
Venue: Accepted by TMLR ()

Abstract: Practical implementations of Gaussian smoothing algorithms have received a great deal of attention in the last 60 years. However, almost all work focuses on estimating complete time series (``fixed-interval smoothing'', $\mathcal{O}(K)$ memory) through variations of the Rauch--Tung--Striebel smoother, rarely on estimating the initial states (``fixed-point smoothing'', $\mathcal{O}(1)$ memory). Since fixed-point smoothing is a crucial component of algorithms for dynamical systems with unknown initial conditions, we close this gap by introducing a new formulation of a Gaussian fixed-point smoother. In contrast to prior approaches, our perspective admits a numerically robust Cholesky-based form (without downdates) and avoids state augmentation, which would needlessly inflate the state-space model and reduce the numerical practicality of any fixed-point smoother code. The experiments demonstrate how a JAX implementation of our algorithm matches the runtime of the fastest methods and the robustness of the most robust techniques while existing implementations must always sacrifice one for the other.

Published in Transactions on Machine Learning Research (01/2025)
Numerically Robust Fixed-Point Smoothing Without State Augmentation
Nicholas Krämer
pekro Oduu. dk
Technical University of Denmark
Kongens Lyngby, Denmark
Reviewed on OpenReview: https: opennevii. net/ forun? Ld-LVEEE
Abstract
Practical implementations of Gaussian smoothing algorithms have received a great deal of attention in the last 60 years. However, almost all work focuses on estimating complete time series ("fixed-inderoal smoothing", O(K) memory) through variations of the Rauch-Tung- Striebel smoother, rarely estimating the initial states ("fixed-point smoothing", O(I) memory). Since fixed-point smoothing is a crucial component of algorithms for dynamical systems with unknown initial conditions, we close this gap by introducing a new formulation of a Gaussian fixed-point smoother. In contrast to prior approaches, our perspective admits a numerically robust. Cholesky-based form (without downdates) and avoids state augmentation, which would needlessly inflate the state-space model and reduce the numerical practicality of any fixed-point smoother code. The experiments demonstrate how JAX implementation of our algorithm matches the runtime of the fastest methods and the robustness of the most robust techniques while existing implementations must always sacrifice one for the other.
Code:
1 Introduction
Linear Gaussian state-space models and Bayesian filtering and smoothing enjoy numerous applications in areas like tracking, navigation, control (Grewal and Andrews, 2014; Särkkä and Svensson, 2023). They also serve as the computational backbone of contemporary machine learning methods that revolve around time-series, streaming data, or sequence modelling: for example, message passing, Gaussian processes, and probabilistic numerics (Grewal and Andrews, 2014; Särkkä and Solin, 2019; Murphy, 2023; Hennig et al., 2022). Filtering and smoothing have also been used for constructing and training neural networks (Singhal and Wu, 1988; Gu and Dao, 2023; Chang et al., 2023) or continual learning (Sliwa al., 2024). All of these applications require fast and robust algorithms for state-space modcls; we propose new in this paper.
Recent developments in Bayesian smoothing focus subjects like linearisation (e.g. et al., 2016), temporal parallelisation (e.g. Särkkä and implementations (e.g. Yaghoobi et al., 2022). However, they all exclusively focus on fixed-interval smoothing, which targets the full time-series p(zo:K y1:k) in O(K) memory, never on fixed- smoothing, which only yields the initial value p(ro y1:k) but does so in O(1) memory. Even though fixed-point smoothing a crucial component of, for example, estimating past locations of a spacecraft (Meditch, 1969), it has yet to receive much attention in the literature state estimation in dynamical systems. The experiments in Section 4 demonstrate fixed-point smoothing applications in probabilistic numerics and in tracking problem. We anticipate that probabilistic mumerical solvers for differential equations will especially benefit from improving the practicality of fixed-point smoothers. The reason is that these algorithms closely connect to filtering and smoothing algorithms (Schober et al., 2019), especially the numerically robust kind (Krämer and Hennig, 2024), and that unknown initial conditions are typical for paramcter estimation problems in differential equations and scientific machine learning (Rackauckas et al., 2020). Section 4.2 revisits probabilistic numerics prime application for fixed-point smoothing.
2020), or numerically robust

---

Published in Transactions on Machine Learning Research (01/2025)
Contributions We introduce implementation of numerically robust fixed-point smoothing. Our approach avoids the typical construction of fixed-point smoothers via siate-augmentation (Biswas and Mahalanabis, 1972; Smith and Roberts, 1982), which increases the dimension of the state model and thus makes estimation needlessly expensive. And unlike previous work on fixed-point smoothing (Moditch, 1967a; Särkkä and Hartikainen, 2010; Rauch, 1963; Meditch, 1967b; 1976; Meditch and Hostetter, 1973; Nishimura, 1969), our perspective is not tied to any particular paramctrisation of Gaussian variables. Instead, our proposed algorithm enjoys numerical robustness and compatibility with data streams, while maintaining minimal complexity. There exist tools for each of those desiderata, but our algorithm is the first to deliver them all at once. We demonstrate the algorithm's 's efficiency and robustness on sequence of test problems, including a probabilistic mumerical method for boundary value problems (Krämer and Hennig, 2021), and show how to use fixed-point smoothing for parameter estimation in Gaussian state-space models.
Notation Enumerated sets are abbreviated with subscripts, for example To:K = (f.....,KKk) TK) and Vi-K == (..,...... .... .x For sequential conditioning of Gaussian variables (like in the Kalman filter equations), indicate the most. recently visited data points with subscripts in the parameter vectors, for example, p(a* (1*k-1) (N(mejs-1, Cz (m CK k) p(æk 91.k). Ln is the identity matrix with n rows and columns. All covariance matrices shall be symmetric and positive semidefinite. Like existing work on numerically robust state-space model algorithms (e.g. Grewal and Andrews, 2014; Krämer, 2024), we define the generalised Cholesky factor by of a covariance matrix Σ as any matrix that satisfies = [xx( This definition includes the "true" Cholesky factor Σ is positive definite but applies to semidefinite Σ; for instance, the zero-matrix is the generalised Cholesky factor of the zero-matrix (which does not admit a Cholesky decompasition). There numerous ways of parametrising multivariate Gaussian distributions, but exclusively focus two kinds. Like Yaghoobi et al. (2022), we distinguish:
Covariance-baste parametrisations: Paramotrisc multivariate Gaussian distributions with means and covariance matrices. Covariance--usee parametrisattons are the slandard approach.
Cholesky-basee based parametrisations: Parametrise multivariate Gaussian distributions with means and generalised Cholesky factors of covariance matrices instead of covariance matrices. Manipulating Gaussian variables in Cholesky-based paranctrisations replaces the addition and subtraction of covariance matrices with QR decompositions, which improves the numerical robustness at the cost of a slightly increased runtime; it leads to methods like the square-root Kalman filter (Bennett, 1965).
Other forms, such as the information canonical form of multivariate Gaussian distribution (Murphy, 2022), are not directly relevant to this work.
2 Problem statement: Fixed-point smoothing
2.1 Background filtering and smoothing
Linear Gaussian state-space models This work only discusses linear, discrete state-space models with additive Gaussian noise because such models are the starting point for Bayesian filtering and smoothing. Nonlincar extensions are future work. For integers d and D, let Ak, E DxD and HI..., Hx E d be linear operators and Coo, B,,.., Bx E DxD and R,,.., Rx E [dddd be covariance matrices. Introduce a vector mojo E RDD Assume that observations yıkk are available according to the state-space model (and potentially as a data stream)
To =, 0, Ik A ba yh = HATE + Tk, k= 1,.,..,
(1)
with pairwise independent Gaussian random variables
0 N((moo Cojo) by N/ Bk), TK N7 RR), k ==
(2)
In order to simplify the notation, Equations and 2 assume that the Gaussian variables b and Fı:K have zero mean and that there is no observation of the initial state vo 0. The experiments in Section 4
2

---

Published in Transactions on Machine Learning Research (01/2025)
Table 1: Estimation in Gaussion state-space models. Other estimation tasks exist, for example, fixed-lag smoothing. However, this article focuses on fixed-point smoothing and its relationship with filtering and fixed-interval smoothing. "RTS smoother": "Rauch Tung Striebel smoother".
Task
Target
Filtering (e.g. Kalman filter)
Real-time updates {p(pa I Via))
Fixed-interval smoothing (e.g. RTS smoother)
Full time series p(ro:K y1.k) Initial state p(2o y1:k)
Fixod-point smoothing
demonstrate successful fixed-point smoothing even when violating those two assumptions. Estimating Fı:K
from observations Vi-K is a standard setup for filtering and smoothing (Särkkä and Svensson, 2023).
Kalman filtering Different estimators target different conditional distributions (Table 1). For example, the Kalman filter (Kalman, 1960) computes {p(a* ()=k)) by initialising p( (1:0) = p(xr)) and alternating
Prediction:
p(r=-1 Vi:k-1) - p(rk V1-k-1) p(ak B1 - p(ax yisk)
k = 1,.... K
(3a)
Update:
k = 1, 1=....,K K.
(3b)
Since all operations are linear and all distributions Gaussian, the recursions are available closed form. Detailed iterations are in Appendix A. An essential extension of the Kalman filter is the square-root Kalman filler (Bennett, 1965; Andrews, 1968) (see also Grewal and Andrews (2014)), which yields the identical distributions the Kalman filter but manipulates Cholesky factors of covariance matrices instead of covariance matrices. The advantage of the Cholesky-basee implementation of the Kalman filter over the covariance-based version is that covariance matrices are guaranteed to remain symmetric and positive semidefinite, whereas accumulated round-off errors can make the covariance- based Kalman filter break down. In this work, we introduce Cholesky-basee implementations for fixed-point smoothing, among other things. Both the Kalman filter and its Cholesky-based extension cost (3) runtime and O(D²) memory, assuming D 2 d (otherwise, exchange D for d). The Cholesky-based filter is slightly more expensive because it uses QR decompositions instead of matrix multiplication; Appendix A contrasts Cholesky- based and covariance- based Kalman filtering.
Fixed- interval smoothing While the filter computes the terminal state p(zK }1:k), the fized-interval smoother targets the entire trajectory p( (1.k). For linear Gaussian state-space models, fixed-interval smoothing is implemented by the (fixedinntevall) Rauch-Tung Taihh-TungStiii -trng-Se smoother (Rauch et al., 1965): Rolying on independent noise in Equation 2, the Rauch Tung Striebel smoother factorises the conditional distribution backwards in time according (interpret y1.0 0 for notational brevity),
K
p(Z0:k (1:k) p(ak (..k) II p(aR-1 Fk, B1:k-1). k=1
(4)
The first term, p(Ix (1:k), is the filtering distribution at the final state and usually computed with Kalman filter. The remaining terms, p(za-1 Tk, V1-A-1), 1), can be assembled with the prediction step in the filtering pass (Equation 3a), which preserves the Kalman filter's O(KD*) runtime complexity but increases the memory consumption from O(D?) to O(KD?). Like the Kalman filter, the Rauch-Tung Stricbel covariance- and Cholesky-based parametrisations (Park and Kailath, 1995) in roughly the same complexity. In practice, Cholesky- arithmetic costs slightly more than bovarianco-based arithmetic but enjoys an increase in numerical robustness. Appendix B contrasts both implementations.
Tung--ttieee smoother allows
Towards fixed-point smoothing Finally, we turn to the fized-point smoothing problem: computing the conditional distribution of the initial state p(zo y1:k) conditioned on all obscrvations (Figure 1). The central difficulty for fixed-point smoothing is that we ask for O(K) algorithms even though all observations Vi-k are in the "future" of wo, which lacks an immediate sequential formulation. That said, one algorithm for scquentially solving the fixed-point smoothing problem involves Rauch-Tung-SS Stricbel smoothing: The fixed-point smoothing problem could be solved by assembling p(Zo:K y1-k) with the Rauch Rauch-Tung- Striebel smoother, which yields parametrisation of p(ao yı:k) in closed form. Unfortunately, this solution is not.
3

---

Published in Transactions on Machine Learning Research (01/2025)
1/1
OK
U1
yK
ID
I
OK
30
I
JK
QoI: Fixed-point
QoI: Fixed-interval
Figure 1: Fixed-point (left) versus fixed-interval smoothing problem (right) factor graphs. The shaded
variables are observed. "Qol": "Quanlity of interest".
very efficient: the Rauch Tung Striebel smoother stores parametrisations of all conditional distributions {p(*x-1 ZR. (11a-1))k= which is problematic for long time-scries because it consumes O(KD2) memory. In contrast, state-of-the-ann fixed-point smoothers require only O(D²) memory.
2.2 Existing approaches to fixed-point smoothing
The stale of the in fixed-point smoothing does not fixed-interval smoothers to compute p(a0 [xi:)) There two more efficient approaches: Either solve the fixed-point smoothing problem by computing the solution smoothing via state-augmented filtering" below), we derive recursions for how p(zo 91:k) evolves with = 1,...., K (Meditch, 1967b;a; 1969, "fixed-point recursions" " below). Both solutions have advantages and disadvantages. We discuss both before explaining how our algorithm fixes their shortcomings.
specific high-dimensional filtering problem (Biswas and Mahalanabis, 1972, "fixed-point
Fixed- Pointpoint smoothing via augmentd filtering Fixed-point smoothing can be implemented with a Kalman filter, more precisely, by running filter on the state-space model (Biswas and Mahalanabis, 1972)
It Tn
A&
0
xxx-
+
(40)
LD
bx: yk (Hs
0)
Ikk TD
+ rk: k k= L,.,,., K
(5)
0
Tp
with initial condition
( (ojo))
Cooo Cojo
Coju
)
[Mojo LCo. Lc
0)
LCoja LCoi
2) 9))
0
(6)
20
~
Cojo
=
(mopo
The covariance of the initial condition rank-deficient. However, the generalised Cholesky factor remains well-defined according to Equation 6. The difference between the state-space models in Equation and Equation 1 is that Equation 5 tracks the augmented state (Ir, instead of 24. As result, the Kalman filter applied to Equations and yiclds paramctrisations of the conditional distributions {p(# zo ) K in ((K(2))) runtime and (((2))² memory. The factor "2D" stems from doubling the state-space dimension. The initial condition p( y1-K) emerges from the augmented state p(rK, zo 91:K) ) in closed form. Relying on the Kalman filter formulation inherits all computational benefits Kalman filters, most importantly, the Kalman filter's Cholesky-based form with its attractive numerical robustness. Nonetheless, doubling of the state-space dimension unfortunate bocause it incrcases the runtime roughly by factor 8 (because (2D)² 8D*) and the memory by factor 4 (because (2D)² 4D²). However, long as one commits to covariance-based arithmetic, this increase can be avoided:
Fixed-pole point recursions Suppose we temporarily ignore Cholesky- forms and commit to covariance- based parametrisations of Gaussian distributions. In that case, we can improve the efficiency of the fixed-point smoother. turns out that during the forward pass, the mean and covariance of p(ro y1sk) follow a recursion that can be run parallel to the Kalman filter pass (Meditch, 1969; Särkkä and Hartikainen, 2010; Meditch, 1967b). Särkkä and Hartikainen (2010) explain how implement this recursion: start by initialising mojo and Cojo according to Equation and set Gojo T. Then, iterate from (moj--1> Cojk-1 Go to
4

---

Published in Transactions on Machine Learning Research (01/2025)
Table 2: Existing approaches to the fixed-point smoothing problem. The "fixed-point recursion" represents Meditch (1969; 1967b) through Equation 7. (1(1) memory" stands in contrast to the O(K) memory of fixed-interval smoother, and the "low-dimensional state" refers to avoiding aggeeauugennttion,, which affects both runtime and memory because it doubles the dimension of the state-space model.
Method
( memory
Low-dimensional state Cholesky-based
Via Rauch Rauch-Tung Striebel smoothing Via filtering with augmented state Fixed-point recursion (Equation 7)
x
x
/
This work (Algorithms to 3)
/
(mojk: Capes Gojk) using the recursion (Särkkä and Hartikainen, 2010; Särkkä, 2013)
Gojs Go k- Gk-1
(7a)
mojk = majk-1 + Gojz ( 11](=11)
(7b)
Coje Co Goj CCIIk-)((GO)))
(7c)
Next to (Moj--1, Cojk-1 (ojR--, Gojk-1), these formulas only depend on the output of the prediction step as well as the smoothing gains Gk-1jk (smoothing gains: Appendix B). For k the fixed-point gain intentionally equals the smoothing gain, which Equation 7 expresses by initialising Gojo = I. The recursion in Equation 7 can be implemented to run simultancously with the forward filtering pass. Unfortunately, even though this implementation avoids doubling the size of the state-space model and enjoys O(D²) memory and O(KD*) runtime, a Cholesky-based formulation has been unknown (until now). Thus, Equation cannot applied to problems where numerical robustness is critical, like the probabilistic numerical simulation of differential equations (Krämer, 2024). The main contribution of this paper to generalisc Equation 7 enable Cholesky-based parametrisations:
Contribution. Compute the solution of the fized-point smoothing problem p(zo yı:k) in O(KD*) runtime, O(D²) memory, using any parametrisation of Gaussian variables, and without state augmentation; Table 2.
3
The method: Numerically robust fixed-point smoothing
Our approach to numerically robust fixed-point smoothing involves two steps: First, we derive a recursion for the conditional distribution p(zp Fx Visk) instead of one for the joint distribution p( 2x V1.)), I..,, K. Second, we implement the recursion in Cholesky-basee based parametrisations without losing closed-form, constant- memory updates. As byproduct of this derivation, other parametrisations of Gaussian distributions (like the information form) become possible, too. However, focus Cholesky-based implementations for their numerical robustness and leave other choices to future work.
3.1
A new recursion for fixed- point smoothing
A derivation of new fixed-point smoothing recursion follows. Similar to the state-augmented filtering perspective, the target distribution p(zo yı:k) can be written as the marginal of the augmented state
p(ao Vl:k) = ploo, oK yi.k) drk-
(8)
Computing p(zo Wiik becomes a byproduct of computing p( Φ (ikk), and the latter admits sequential formulation. This perspective was essential to implementing fixed-point smoothing with state-angmented filtering. The augmented state factorises into the conditional
pk B1.K) p(2o æ VI.K)((KK v..K).
(9)
5

---

Published in Transactions on Machine Learning Research (01/2025)
Since p[#K U1:K) is already computed in the forward pass, it suffices to derive recursion for the fixed-poii point conditional p(zo æx. V..K). The backward factorisation of the smoothing distribution in Equation 4 implies
K p( p Ik, Vi:k dr1-----
p(zo εK. Vr:k) = /
(10)
Marginalising over all intermediate states 21:K-1 like in Equation 10 turns Rauch Tung Striebel smoother into fixed-point smoother. An O(K) runtime implementation now emerges from observing how the integral in Equation 10 rearranges to a nested sequence of single- variable integrals,
p IK. B1.K) = / 1I p(ax-1 æus 11) dr1:K-1
K
(11a)
= f V plan xi)p(a 22, vaa) dai plaz I
vas 1a. (1:2) dxz)) dx.
(11b)
f (f p(zo 22, (2.81)((22 aa M1:2) dx2 dæ
(11c)
This rearranging of the integrals is essential to the derivation because it implies forward-in iel---iiittmme recursion: Algorithm 1 (Fixed-point smoother). To compute the solution to the fized-point smoothing problem, assemble p(zk V1:k) with a Kalman filter and evaluate p(ao IK y1:k) ) as follows. (To simplify the index-related notation in this algorithm, read y1:-1 = V1.0 0.)
Initialise p( To. 11) = (Go + Pojo- Poo) with Gojo = ID, Pojo = 0, and Pojo =
For l 1,... K, iterate from is to k;
p(ao æks 311x / p p( 2k 11:k- 2)p( I 11.k-1) dr
(12)
The conditionals p(xx-1 Ib. b1:k-1) are from the Rauch Mrach-Tung-aae Tlun-SStribbel smoother (Equation 4). 3. Marginalise p( V1.K) from p(ak Viik) and p(zo IK. -1.K-1) via Equation 8.
Equation 12 turns a Bauch-Tung-Striees smoother into fixed-point smoother: The recursion requires access to p(r=-1 Øk, (1:k-1) computed by the forward-pass of the Rauch-Tunn Teng-Striebee Striebel smoother. Therefore, Equation 12 runs concurrently to the forward filtering pass.
3.2 Cholesky-based implementation
The missing link for Algorithm 1 is the merging of two affinely related conditionals in Equation 12. Since both conditionals in Equation 12 are affine and Gaussian, their combination p(zo T&i V1:k-1) is Gaussian and available closed form for all 1..... K. Its parametrisation depends on the parametrisation of each input distribution. For covariance- and Cholesky-based arithmetic, it looks follows:
Algorithm 2 (Covariance-based implementation of Equation 12). Recall the initialisation of Algoritla and the convention y1:-1 = y1:0 = 0. For any k = 1...,, K, if the fized point and smoothing conditionals (see Equation 4 Appendix B for the latter) are porometrised by
p(zo (1.k2) -V((o)A- + Pojk- 101--1 Pojk-1)
(13a)
p(#R-1 Ik, 11) -N(Gk-1JRT& + 1|k Pk-1k))
(13b)
for given Gilz E RDDx Pil RDD and Paj E DxD, then, the next flzed-point conditional is
p( Ty 11 = -GooO + Pujks Polk),
(14a)
Goj G k- Gx-
(14b)
Pojk Go k- IPk-1 1 & Pojk-1:
(14c)
Pojk = Gu (inbu---P--I(GGGR--)) Pojk-11
(14d)
Implementing this operation costs ODD3) mouting-poina operations.

---

Published in Transactions on Machine Learning Research (01/2025)
Algorithm 3 (Cholesky-basse based implementation of Equation 12). For any k I..,K iJ the previous fized-point conditional and smoothing transition (see Equation 4 & Appendir B) are
K (again 91:-1 V1.0 0),
p 26-1,31k-2) (G + Pojk-1, LPj (L
(15a)
p( Tke ( (V + Pk-1 Li (Lee-
(15b)
for given Gilz E DDD, Pijj RD, and Lp RDDD, then, the next fixed- conditional is
p(
Tki 11:k11 = =N(Goju + Paja, Pojj (L
(16a)
Gojs = G
(16b)
Polk = 10jk--1P--114 + Pujk-1,
(16c)
LPol =
(16d)
where R. the upper triangular matrix returned by the QR-decomposition
(L. (GGOKK)))) (L .)
QR
(17)
Implementing this operation costs O(D*) floating-point operations.
Algorithm 2 follows from the rules for manipulating affinely related Gaussian distributions (e.g. Särkkä, 2013, Appendix A.1). Algorithm 3 is indeed the Cholesky-based version of Algorithu 2: The expression in Equation 16d is the generalised Cholesky factor the expression in Equation 14d because
(L ( (GoR--))" (L =PP-IIR PR-1IA (18)
II (LP = AT 2 = (Gojk- LPo
holds. The combination of Algorithms and a novel, Cholesky-basee implementation of a fixed-point smoother. Similarly, the combination of Algorithms and is a novel, covariance-based implementation of a fixed-point smoother. With a small modification Algorithm 1, combining Algorithm with the covariance-based formulas in Algorithm recovers Equation 7:
Proposition 1. If the combination of Algorithms and computes the marginals
p(zo v1:k) = f p(zn IE. Di.k 1)p(E vıa) drk.
(19)
at every l: = 1,... K instend of only the final time-step, the recursion reduces Equation 7.
Proof. Induction; Appendix C.
Both Algorithms and cost O(D*) Noating-point operations per step. The kth iteration in Algorithm needs access to the same quantities as the kth iteration of the forward pass of Rauch Tung Striebel smoother, which can be implemented in O(D²) memory. Unlike the Rauch Tung Striebel smoother, Algorithm 1 does not store intermediate results, which is why it consumes O(D?) memory in total (instead of (²))²)) Both Algorithm 2 and Equation 7 require three matrix-matrix: and one matrix-vector multiplication per step, so they are equally efficient (Algorithm 2 uses fewer matrix additions. However, the addition of matrices is fast). Therefore, it will be no loss of significance that we always implement covariance- based fxxed-point smoothing via Algorithm instead of Equation 7 in our experiments.
Computational complexity

---

Published in Transactions on Machine Learning Research (01/2025)
Covariance-based
Cholesky-base
Via filter
Known
Known
Via Rauch-Tung Striebel Stribbtr
Known
Known
Fixed-poing recursion
Known, but Prop.
Our contribution
Experiment II: Robustness
Experiment Efficiency
Figure 2: Outline of the memory-, runtime,, and robustuess-related demonstratious.
4
Experiments
The experiments serve two purposes. To start with, they investigate whether the proposed Cholesky-based implementation (Algorithm 3) of the fixed-point smoother recursion (Algorithm 1) holds its promises about memory, runtime, and numerical robustness. According to the theory in Section 3, we should observe:
Memory: Slightly lower memory demands than a state-augmented Kalman filter; drastically lower memory demands than a Rauch -ung Striebel smoother.
Wall-time: Faster than state-augmented, Cholesky- Kalman filter; roughly fast Rauch-Tung-Strien Sttriee smoother.
Combining Algorithm with Cholesky-based parametrisations (Algorithm 3) significantly more robust than combining it with covarianco-basee parametrisations (Algorithm 2); comparable to Cholesky-based implementation of Kalman filter.
Numerical robustness:
These three phenomena will be studied in two experiments, one for runtime/ muntie/mmmorry memory and one for numerical robustness (outline: Figure 2). The problem set for these two experiments includes a toy problem for the former and application in probabilistic numerics for the latter. Afterwards, case study in tracking shows how to use the fixed-point smoother for estimating the initial parameter in state-space model.
Hardware and code All experiments run on the CPU of a consumer-grade grade laptop and finish within few minutes. Our JAX implementation (Bradbury et al., 2018) Kalman filters, Brrch-TungSSiiee smoothers, and fixed-point smoothers al.
https://github..om wtt.a//.ccooo/aaoooooo
We implement the existing fixed-point smoother recursions Equation 7 by combining Algorithm with covariance-based parametrisations, recall Proposition 1.
4.1 Experiment I: How efficient is the fixed-point recursion?
Sections 2.2 and 3 mention three approaches to fixed-point smoothing: a detour via Rauch Tung Striebel smoothing, state-angmented Kalman filtering, and our recursion in Algorithm 1. In theory, Algorithm should be the most efficient: it does not inflate the state-space model like state-angmented Kalman filter does and requires ((1) instead of O(K) memory, unlike the Rauch Tung Striebel smoother. This first experiment demonstrates that these effects are visible in practice.
Motivation
Problem setup In this first example, we only measure the execution time and memory requirement of three approaches to fixed-point smoothing. Both memory and runtime depend only on the size of a state-space model and not on the difficulty of the estimation task. Therefore, we consider a state-space model where all system matrices and vectors are populated with random values.
Implementation We choose K 1, 000, vary d, set the size of the hidden state D 2d, and use Cholesky-based arithmetic for all estimators. We take Equation I and introduce nonzero bias in all noises

---

Published in Transactions on Machine Learning Research (01/2025)
Table Runtime in seconds (wall-time, best of three runs). Lower better. The column-wise lowest are bold & shaded. Randomly populated model. 1, 000 steps. All mothods use Cholesky-based paramctrisations.
d=2 2
-=5
d= 10
d 20
d 50
d= 100
Via Ranch Tung Striebel
5.8 x 10-3
1.8 x 10-2
6.5 x 10-2
4.5 x 10-1
2.9 x 10
1.2 x 10"
Via filter
5.0 10
2.1 x 10-2
8.0 10-2
7.1 x 10-1
4.6 x 10
1.6 10
Algorithm
6.4 10"a
1.8 x
6.4 10
4.4 x 10
2.9 10
1.2 x 10
Table 4: Memory in bytes (we use 32-bit arithmetic). Lower is better. The column-wise lowest entries are bold and shaded. Randomly populated model. K = 1, 000 steps. All methods use Cholesky-based parametrisations.
d = 2
= 5
d = 10
d= 20
d = 50
= 100
Via Rauch Mauch-TTung-Ti Strribel
2.2 x 10a
1,2 x 10t
4.9 10
1.9 x 107
1,2 x 10
4.8 10
Via filter
2.8 x 10a
1.6 x 10
6.5 10
2.5 x 10*
1.6 x 109
6.4 10
Algorithm 1
2.2 x 10
1.2 x 10*
4.9 x 10a
1.9 x 10
1.2 x 10
4.8 x 10
in Equation 2. Then, we randomly populate all system matrices in the state-space model with independent. samples from (N, 1/K²). Afterwards, we sample yı:K from this state-space model to generate toy data. A K-dependent covariance controls that the samples Vı-K remain finite in 32-bit floating-poin arithmetic. However, the precise values of the model and data do not matter for this experiment only their size does. Finally, we vary d and measure the runtime and memory requirements of each of the three methods. To measure runtime, we use wall time in seconds. We display the best of three runs instead of the average because all codes are deterministic (no data-dependence, differences is background machine noise. Selecting the best of three runs minimises this noise much possible. To mcasure memory consumption, we count the number of floating-poiis values the estimators carry from step to step, multiply this by 32 (because we use 32-bit arithmetic), and translate bits into byles.
randomness). thus the main driver for runtime
Evaluation The runtime results are in Table and the memory results in Table 4. The runtime data shows that except for d 2, the Ranch-Tung-Stricen Sttricbll smoother and Algorithm 1 are faster than the stale-augmented filter. Algorithm though it computes strictly more at every iteration. However, the Rauch- Tung Striebel smoother executes two loops, one forwards and one backwards, whereas Algorithm 1 only executes forward loop. This discrepancy might lead to the performance improvement from Rauch-Tung- Striibel smoothing to Algorithm 1. The memory data shows that the Rauch- Tung Striebel smoothing solution and the Algorithm have identical Tung-Striebel memory requirements per step. Both consume less storage than the state-augmented filter (per step). As expected, the Rauch Tung Striebel smoother requires exactly K-times the memory of Algorithm 1, which would be infeasible for long time series over nigh-dimensional states. In general, this experiment underlines the claimed efficiency of our new fixed-point smoother recursion in Cholesky-based parametrisation. The following experiment will answer the question of whether Cholesky-based arithmetic is necessary.
marginally faster than the Rauch Tung Striebel smoother code, even
4.2 Experiment II: How much more robust is the Cholesky-based code?
Motivation Cholesky-based implementations replace matrix addition and matrix multiplication with QR decompositions. They are more expensive than covariance-based implementations because the decomposition is more expensive than matrix multiplication. However, there are numerous situations where the gains in numerical robustness are worth the increase in runtime. One example is the probabilistic numerical simulation of differential equations, which use integrated Wiener processes (Schober et al., 2014; 2019) with high order and small timo-stops. For these applications, Cholesky-based implementations are the only feasible approach (Krämer and Hennig, 2024). This second experiment shows that Algorithm 1 unlocks point smoothing for probabilistic numerical simulation through Algorithm 3.
9

---

Published in Transactions on Machine Learning Research (01/2025)
Problem setup
We solve boundary value problem based on ordinary differential equation.
More specifically, we solve the 15th in the collection of test problems by Mazzia and Cash (2015) (Figure 3),
& di a - tu((), --1)==(1)) 11
10
(20)
Output
We follow the procedure for solving boundary value problems with a probabilistic numerical method by Krämer and Hennig (2021) for the most part. However, we skip the iterated Kalman smoother because the differential equation is linear and skip mesh refinement and expectation maximisation to keep this demonstration simple. The focus lies numerical robustness. We choose twice-inte vrice-integrated Wiener process prior and discretise it on K equispaced points in [-1, I], which yields the latent transitions (let At == 1/K)
Input
Figure 3: 15th Boundary value
problem (Maziia and Cash, 2015).
(A1)
(A
(Ar)
(At)*
A ==
At
1
Ål
2 B₁
(At)
20
(Ar)
(At)²
h =
L,,,.., K.
(21)
8
2
1
(At)a
(A
A
The dynamics are constant because the grid points are evenly spaced. The means of the process noise are zero, like in Equation 2. The hidden state
d d F& = u((x)) u(tx),
dt dez ultaa) k= L,,,.., K
(22)
tracks the differential equation solution u and its derivatives, including ar u, which means that the residual 10-3 d² u- - tu(t) a linear function of rg. We introduce the model for the constraints (3x shall be the nonzero mean of the observation noise Th in Equations and 2)
H& = ( (-tx tk 0 10 10-
Bx = O,
Rg = 0,
k= == K
(23a)
Hx = (1 0 0).
8K = -1,
Rx = O.
(23b)
Note how the constraint at the final time-point encodes the right-hand side boundary condition, not the differential equation. We choose the initial mean maio 0, and initial covariance Cojo diag(0, 1) to represent the left-hand side boundary condition. Estimating Fo.k from yı.k == 0 solves the boundary value problem ((Krämer and Hennig, 2021); Figure 3 displays the mean of the fixed-interval smoothing solution). Krämer and Hcnnig (2021) explain how ostimating the initial condition p(an (1:k) is important for parameter estimation problems. In other words, fixed-point smoothers are relevant for boundary value problem simulation. We are the first to use them for this task.
Implementation For this demonstration, consider the Cholesky-basee implementation of the Kalman filter the gold standard for numerical robustness because it has been used successfully in numerous similar problems (Krämer and Hennig, 2024; Bosch et al., 2021; Krämer and Hennig, 2021; Krämer, 2024). The Cholesky- state-augmented Kalman filter provides reference solution of the fixed-point equations. We run Cholesky-based (Algorithm 3) and covariance- based code (Algorithm 2) for the fixed-point smoother recursions (Algorithm 1) and measure how much the estimated initial means deviate from the Kalman-filter reference terms of the absolute root-mean-square error. In infinite precision, the results would be identical. In finite precision, all deviations should be due to a loss of stability. We vary K to investigate how the step-size At affects the results, with the intuition that smaller steps lead to worse conditioning in B& and that this effect makes the estimation task more difficult (Krämer and Hennig, 2024).
Evaluation The results are in Table 5. They indicate how the Cholesky-based code is significantly more robust than the oovariance-bssed code. The covariance-based implementation delivers only three meaningful digits for K 10 points, diverges for more grid points, and breaks for K 500. In contrast, the Cholesky- based code delivers meaningful approximations across all grids. This consistency underlines how much more robust our Cholesky-based code is and how it unlocks fixed-point smoothing for probabilistic numerics.
10

---

Published in Transactions on Machine Learning Research (01/2025)
Table 5: Deviation from the Cholesky-based, state-augmented Kalman filter on the boundary value problem.
Lower is better and close to machino-procisson desirable. The column-wise lowest values are bold and coloured. " "K": number of grid points. Double precision (64-bit arithmetic).
--based:
10
K 50
K = 100
200
K 500
K = 1, 000
Covariance
2.9 x 10 -3
3.4 x 10-1
1.1 x 10"
1.0 x 10a
2.1 x 10¹
NaN
NaN
Cholesky
2.0 x 10
5.0 10
4.2 10"
7.9 10
1.3 10"
6.1 x 10
3.4 10
4.3
Case study: Estimating parameters of a state-space model
The previous two experiments have emphasized the practicality of our proposed method. We conclude by applying fixed-point smoothing to parameter estimation in tracking model. This study aims emulate applications of fixed-point smoothing in navigation tasks (Meditch, 1969) in setup that is easy to reproduce. However, the demonstration also links to the previous boundary value problem experiment through expectation maximisation (Krämer and Hennig, 2021).
Movitation
Problem setup The task is to estimate the mean parameter of an unknown initial condition in a car tracking example. Define the Wiener velocity model on K 10 equispaced points (At = 1/10),
A = 1/2
12
At-12 12 12
(Ar)
3 l
(Ae) 12 HA (I2 0), Rs = 0.12 I2-
= B₁
(24)
(At)
2 12
At.
All biases are zero, by = 0, rh = 0. Textbooks on Bayesian filtering and smoothing (Särkkä, 2013; Särkkä and Svensson, 2023) use this Wiener velocity model to estimate the trajectory of car. We populate m and Le with samples from a standard normal distribution and sample an initial condition
To (01
02
01 02) (mmojo Lago (Looj))
(25)
We sample artificial observations yı.k. From here on, the initial mean maja will be treated as unknown.
Implementation We combine fixed-point smoothing with expectation maximisation (Dempster et al., 1977) to calibrate mojo using the data y1.k. The expectation maximisation update for the initial mean in a linear Gaussian state-space model is mnw == mojk (Särkkä and Svensson, 2023), and this update repeats until convergence. Here, Majk is the mean of p(zo k). We implement the fixed-point smoother recursion in Cholesky- arithmetic and run expectation maximisation for three iterations. We initialise the mean guess by sampling all entrics independently from a centred normal distribution with a variance of 100. We track the data's marginal likelihood ("evidence"), computing it online during the forward filtering pass.
The results are in Figure 4. They show how the combination of expectation maximisation with fixed-point smoothing recovers the initial mean alrcady after three iterations. The evidence increases at every iteration, but that is normal for expectation maximisation (Wu, 1983). In conclusion, fixed-point smoothing is a viable parameter estimation technique in a state-space model.
Evaluation
5 Discussion
Limitations and future work Our new approach to fixed-point smoothing strictly generalises existing techniques, but inherits some of their limitations: Even though the new recursion in Algorithm is independent of the type of state-space model, implementing the fixed-point smoother in closed form assumes a linear Gaussian setup. Future work should explore robust fixed-point smoothing in nonlinear state-space models, for example, through posterior lincarisation (Carcía-Feráándoz et al., 2016). Like all Gaussian smoothing algorithms, our methods have cubic complexity in the state-space dimension because of the matrix-malrix arithmetic or QR decompositions, respectively. Finally, other parametrisations Gaussian variables could be used with Algorithm 1 instead of covariance or Cholesky- paramotrisations ( (Algorithms 2 and 3); for example, the information form or ensembles (Murphy, 2022; Houtekamer and Mitchell, 2005). Combining these alternatives with numerically robust fixed-point smoothing could be interesting avenues for future work.
11

---

Published in Transactions on Machine Learning Research (01/2025)
Evidence: -17.92
i =
Evidence: 0.11
i =
Evidence: 0.13
PDF (8x
20
Ttairaio
Target
Target Nolay data
Target
Notay data
Notsy data
0..4
-0.3
-0.2
-0.1
0.0
0.4
0.3
-022
-0.1
0.0
- 0..4
-
0.3
-022
-0.1
0.0
Tterrddo Target Noisy data
Thwrate
Iterate
PDF (8s)
Target
Target
20
Nolay
Nolay data
0.88
-..7
-0.6
00.5
0.4
0.8
U.7
0.6
-0.5
0.4
0.8
U.7
0.6
00.5
-0.4
Realisation
Realisation
Realisation
Figure 4: Initial distributions p(01,02 Oz v1:k) of the car tracking model after rumning the fixed-point smoother (recall x = (0,, 02, 0,, 02)). Left to right: After three iterations, the combination of expectation maximisation with the fixed-point smoother finds the correct initial mean 0 = (0), (2) the state-space model. Top: First coordinate p(0 f1.k). Bottom: second coordinate p(02 V1.k). "PDF": "Probability density function".
Conclusion This paper presented a new recursion for fixed-point smoothing in Gaussian stale-space models: Algorithm 1. It has lower memory consumption than existing approaches because it avoids state-augmentation and foregocs storing all intermediate results of Rauch-Tung -ure-Striebb parametrisations of Gaussian distributions, and we use this perspective to derive a Cholesky-based form of fixed-point smoothers. As result, our method matches the speed of the fastest and the robustness of the most robust methods, as has been demonstrated in three simulations of varying difficulty. Through these successes, our contribution hopefully revives fixed-point smoothing as part of the literature on Gaussian filters and smoothers. We anticipate notable performance improvements for algorithms the interface smoothing and dynamical systems due to the newfound ability to reliably use non-standard smoothing algorithms.
Stricbel smoother. It also allows arbitrary
Acknowledgements
This work was supported by a rescarch grant (42062) from VILLUM FONDEN. The work was partly fundod the Novo Nordisk Foundation through the Center for Basic Machine Learning Research in Life Science (NNF200C0066606), This project received funding from the European Research Council (ERC) under the European Union's Horizon programme (grant agreement 101125993).
References
Mohinder S Grewal and Angus P Andrews. Kalman filtering: Theory and Practice with MATLAB. John Wiley & Sons, 2014.
Simo Särkkä and Lennart Svensson. Bayesian Filtering and Smoothing. Cambridge University Press, 2023.
Simo Särkkä and Arno Solin. Applied Stochastic Differential Equations, volume 10. Cambridge University Press, 2019.
Kevin P Murphy. Probabilistic Machine Learning: Advanced Topics. MIT Press, 2023.
Philipp Hennig, Michael A Osborne, and Hans P Kersting. Probabilistic Numerics: Compulation as Machine Learning. Cambridge University Press, 2022.
Sharad Singhal and Lance Wu. Training multilayer perceptrons with the extended kalman algorithm. Advances in neural information processing systems, 1, 1988.
Albert Gu and Tri Dao. Mamba: Linear-time sequence modeling with selective state spaces. arXiv preprint 2Vo:2312.00752 2023.
12

---

Published in Transactions on Machine Learning Research (01/2025)
Peter G Chang, Gerardo Durán-Martin, Alex Shestopaloff, Matt Jones, and Kevin Patrick Murphy. Low-rank extended kalman filtering for online learning of neural networks from streaming data. In Conference on Lifelong Learning Agents, pages 1025 1071. PMLR, 2023.
Joanna Sliwa, Frank Schneider, Nathanacl Bosch, Agustinus Kristiadi, and Philipp Hennig. Efficient weight-space Laplace Gaussian filtering and smoothing for sequential deep learning. arXio preprint arXiv:2410.06800, 2024.
Angel F Garnii--errnndee, Lennart. Svensson, and Simo Särkkä. Iterated posterior linearization smoother.
IEEE Transactions of Automatic Control, 62(4):2056 2063, 2016.
Simo Särkkä and Angel F Farríi-Fernández. Temporal parallelization of Bayesian smoothers. Transactions of Automatic Control, 66(1):299 306, 2020.
IEEE
Fatemeh Yaghoobi, Adrien Corenflos, Sakira Hassan, and Simo Särkkä. Parallel square-root statistical linear regression for inference in nonlincar state space models. Preprint on (rXi::2277.04200 2022.
James S Meditch. Stochastic Optimal Linear Estimation and Control. McGraw-Hill, 1969.
Michael Schober, Simo Särkkä, and Philipp Hennig. probabilistic model for the numerical solution of initial value problems. Slalislics and Compuling, 29(1):99-18))) 122, 122, 2019.
Nicholas Krämer and Philipp Hennig. Stable implementation of probabilistic ODE solvers. Journal of Machine Learning Research, 25(111):1-29, 2024.
Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Domínic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scientific machine learning. arXiv preprint ar xrx::2001.04388 2020.
KK Biswas and AK Mahalanabis. An approach to fixed-point smoothing problems. IEEE Transactions on Aerospace and Electroníc Systems, (5):676-682) 1972.
MWA Smith and AP Roberts. The robustness the fixed point smoothing algorithm. Journal of Information and Optimization Sciences, 3(2):87-106 104, 1982.
JS Meditch. Optimal fixed-point continuous linear smoothing. In Joint Automatic Control Conference, number 5, pages 249 257, 1967a.
Simo Särkkä and Jouni Hartikainen. On Gaussian optimal smoothing of non-linear state space models. 1EEE Transactions of Aulomalic Control, 55(8):1938 1941, 2010.
H Rauch. Solutions to the linear smoothing problem. IEEE Transactions of Automatic Control, 8(4):311 372, 1963.
James S Meditch. Orthogonal projection and discrete optimal linear smoothing. SIAM Journal on Control, 5 (1):74-88 19, 1967b.
James Meditch. Initial- and lagging-state observers. Information Sciences, 11(1):55-00 67, 1976.
JS Meditch and GH Hostelter. Estimation of initial conditions in linear systems. In Proceedings of the 1 11th Annual Atierton Conference on Circuit and System Theory, pages 226-235, 1973.
T Nishimura. new approach to estimation of initial conditions and smoothing problems. JEEF Transactions on Aerospace and Electronic Systems, (5):828 836, 1969.
Nicholas Krämer and Philipp Hennig. Lincar-time probabilistic solution of boundary value problems. Advances in Neural Information Processing Systems, 34:11160 11171, 2021.
Peter Nicholas Krämer. Implementing probabilistic numerical solvers for differential equations. PhD thesis, Universität Tübingen, 2024.
13

---

Published in Transactions on Machine Learning Research (01/2025)
John M Bennett. Triangular factors of modified matrices. Numerische Mathematik, 7:217-221, 1965.
Kevin P Murphy. Probabilistic Machine Learwing: an Inlroduction. MIT press, 2022.
Rudolph E Kalman. A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1):35-45 45, 1960.
Angus Andrews. A square root formulation of the kalman covariance equations. AIAA Journal, 6(6): 1165-1166, 1968.
Herbert E Rauch, F TTng, and Charlotte T Striebel. Maximum likelihood estimates of linear dynamic systems. ATAA journal, 3(8):1445-1450, 1965.
Poogycon Park and Thomas Kailath. Square-root RTS smoothing algorithms. International Journal of Control, 62(5):1049 1060, 1995.
Simo Särkkä. Boyesion Filtering and Smoothing. Cambridge University Press, 2013. (Note: This is the first edition. The second edition (Särkkä and Svensson, 2023) does not discuss fixed-point smoothing.).
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang, JAX: composable transformations of Python+NumPy programs, 2018. URL http://github. com/ wtt.cccm//oo..eee google/ jangle/ja...
Michael Schober, David K Duvenaud, and Philipp Hennig, Probabilistic ODE solvers with Runge- Kutta means. Advances in Neural Information Processing Systems, 27, 2014.
Francesca Mazzia and Jeff R Cash. A Fortran test sot for boundary value problem solvers. In AIP Conference Proceedings, volume 1648. AIP Publishing, 2015.
Nathanael Bosch, Philipp Hennig, and Filip Tronarp. Calibrated adaptive probabilistic ODE solvers. Tn International Conference Artificial Intelligence and Statistics, pages 3466 3474. PMLR, 2021.
Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodologieal), 39(1):1-22, 1977.
CF Jeff Wu. On the convergence properties of the EM algorithm. The Annals of Statistics, pages 95 95-103, 1983.
Peter L Houtckamer and Herschel L Mitchell. Ensemble Kalman filtering. Quarterly Journal of the Royal Meteorological Sociely: Journal of the Almospheric Sciences, Applied Meteorology and Physical Oceanography, 131(613):3269 3289, 2005.
Stuart Gibson and Brett Ninness. Robust maximum-likelihood estimation of multivariable dynamic systems. Automatica, 1166 2005.
Matthias Seeger. Low rank updates for the Cholesky decomposition. 2004.
A Kalman filter recursions
The Kalman filter (Kalman, 1960) computes p(rx B1.K) by initialising p(æo V1.0)) = p(zo) and alternating
Prodiction:
p(ak-1 V1:k-1) p(akk V1:k-1) p(ak Vi - p(ax Vick)
k 1,..., K
(26a)
Update:
k = 1, K.
(26b)
Since all operations are linear and all distributions are Gaussian, the recursions are available in closed form. Their parametrisation depends the parametrisation of Gaussian variables follows.
14

---

Published in Transactions on Machine Learning Research (01/2025)
A.1 Covariance-bass based parametrisation
Let -.....,, K. The prediction step computes the predicted distribution p(ux 11 (mmaj--1 CA(R-1) from the previous filtering distribution p(a=-1 Visk-1) (mm==1J Ck- by
m 1Ij=-1 CaJE-1 = C (At)) + BB.
(27)
The update step takes the predicted distribution, forms the joint distribution
p yk 11*x-1) p(u #R)P(I* V1-8-1) = - RR)NNmaas CHRK-1)
(28)
and computes the update as p(æ y1:k) N((MjR CMRR);
S HymAJK-1.
(29a)
S HHCCHR-I (HR)) RSs
(29b)
Z% = CRJK- (H)x) (Sxkk-
(29c)
malk maKk-1 Zz (y/: 3k(1k-1
(29d)
Cajk = CHJR-1 - - (Zk)
(29c)
Iterating these two steps from ho 1, K yields p(zx I V1:k).
A.2 Cholesky-based parametrisation
The Cholesky-based parametrisation predicts the mean like the covariance-based parametrisation. The covariance prediction is replaced by the QR decomposition
DR (1 (LO
(L (T
(30)
followed by setting T-CHA-I Re (Krämer and Hennig, 2024; Grewal and Andrews, 2014). The logic mirrors that in Algorithm 3. The update step in Cholesky-based arithmetic amounts to QR-decomposition of (Gibson and Ninness, 2005)
Ri R2 2
(LRe) ((LC)
(HK)) (LCmaa))T)
(31)
R3
followed by setting
LS ( Zu (())) (19))-1 LocA = -PRa)
(32)
Then, use Z& to evaluate majk and iterate - k + 1. If p(yk V1.k-1) is needed, evalaute Suj=1 like in the covarianco-based implementation. Choosing the QR dccomposition by Gibson and Ninness (2005) avoids implementing Gaussian updates via Cholesky downdates (e.g. Yaghoobi et al., 2022), which are sometimes numerically unstable (Seeger, 2004). We refer to Krämer (2024, Chapter 4) for why the above assignments yield the correct distributions.
B Rauch-Tung-SSte Striebel smoother recursions
The Rauch Tung Striebel smoother proceeds similarly to the Kalman filter.
B.1 Covariance-based parametrisation
The prediction step the smoother involves computing p(rk y1:k 1) like in the filter. It further evaluates
p(xx-1 1 Ik. (1:k-1) =::))==(GGGG) (N 1 + Pk- 1|k, PR-1IR),
(33a)
Gb-1lk = C2-11k-1 (Co)k-1)
(33b)
Pk-1 = m G
(33c)
Pk-1ik = PR- G C (
(33d)
This conditional distribution stored after every step. The rest of the step proceeds like in the Kalman filter.
15

---

Published in Transactions on Machine Learning Research (01/2025)
B.2 Cholesky-based parametrisation
The Cholesky-based paramelrisation compules the same conditional distribution the covariance-based parametrisation but replaces matrix multiplication and matrix addition with another QR decomposition; specifically, the QR dccomposition (Gibson and Ninness, 2005)
Ri
H2 Rz)
(L ((LC)-II )T (Ax)T (Lo))1))7)
(34)
This QR dccomposition is followed by sctting
Le = (9t)", GR-1|K= = ((
Lp--11 (Rs)
(35)
This step computes the Cholesky factors of Ckjt-1 and 1k- well the smoothing gain G--1jK 1|k in a single sweep. It replaces the prediction step of the Cholesky-based Kalman filter. Compute majk-1 and Pa- l/k like in the covariance-based parametrisation. Store these quantities and proceed with the updale step of the Cholesky-based Kalman filter. Again, refer to Krämer (2024, Chapter 4) for why these assignments yield the correct distributions.
C Proof of Proposition
Recall the notation from Algorithm 2, most importantly, Gojo 1, Pojo O, and Pojo = O. Algorithm compules the Iransition
p( Ix V1---1) = (GGo + Poja: Poja)-
(36)
A forward pass with a Kalman filter gives
p(ax V1R) = N(maj Cuj)).
(37)
Then, due to the rules of manipulating Gaussian distributions,
p V1sk) = [Goojmmej + Pojk (onne)) + Polk)
(38)
follows. To show that this matches the result of Equation 7, it sulfices to show
Poja moja-1 Go m
(39a)
Poj Cojk-1 - Goja GoJoC))--((ooss)
(39b)
because the remaining terms already coincide. We induction show Equation 39.
Initialisation For k = 1, we get
Poj1 GojoPoj1 Pojo = Poj1 = mojo moj
(40)
which shows Equation 39a (compare the smoothing recursion in Appendix B for the last step), as well as
Poja Gojo Pojo Pojo Poj1
(41)
which gives Equation 39b.
Step
Now, let Equation 39 hold for some Then, (abbreviate "RTSS": "Rauch-Tune Striebel smoother")
Poja+1 = -GojPPEj- + Pojk
(Algorithm 2)
GoJoPEj-L + mojk-1 GojamBja-1
(Equation 39a holds for k)
Guja(mank - G4R4+1784411) + mojk-1 Gujemaj--
(Expand PkJ++1 as RTSS)
Gojamajk - Go GG
Mojk-1 - Majk-
(Resolve parentheses)
Moja-1 Goo.(mujj- mmJN) Go M+1
(Goj& GojR+1)
majk Gojk-1m4411k
(Equation 7)
16

---

Published in Transactions on Machine Learning Research (01/2025)
Equation 39a is complete. Similarly,
Poje+1 GojePARtti (Goje)) + Poje
(Algorithm 2)
G + COJR-1 CRIoCcIb-1GGoa))
(Equation 39b holds for k))
GujK(C -uJa(CaJ+=-aaHHH
(Gook) + - -R
(Expand PR|K+1 as in RTSS)
(GOJCR Go
+ Cojk-1 Goja (Go
(Expand parentheses, Goj&GE GOJR+1)
Cojk-1 - Gojo k-1 (Go)) Go
(Rearrange) (Equation 7)
Cujk - (onR+1CC+1N(Gooaaa..
the covariance recursion must hold. Equation 39 holds for all 2 1, and the statement is complete.
17


Reviews:
[{"review_id": "HkIgop7dzW", "reviewer": "Reviewer_z28E", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "In the previous section, I outlined the advantages of this work over prior approaches. However, I believe the importance of addressing the fixed-point smoothing problem has not been sufficiently emphasized. Additionally, since this work is closely related to the fixed-point recursions method, particularly equations (7a) to (7c), further explanation is needed as to why the authors did not derive the Cholesky formulations directly from these equations. More explanations will be given in the next section.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work addresses the fixed-point smoothing problem, where the goal is to determine the posterior distribution of the initial state given all observations, specifically $p(x_0 \\mid y_{1:K})$. This problem is particularly important for dynamical systems with unknown initial conditions. The authors focus on a linear, discrete-time state-space model with additive Gaussian noise and present a numerically robust solution based on Cholesky factorization.\n\nA key advantage of this approach over the previous method, \"fixed-point smoothing via space-augmented filtering,\" is that while both rely on Cholesky factorization, the new method eliminates the need for state augmentation, which reduces both runtime and memory usage.\n\n\nThis work is closely related to the \"fixed-point recursions\" method by Särkkä and Hartikainen (2010), which introduced covariance-based recursions for the fixed-point smoothing problem. A key benefit of Särkkä and Hartikainen's method is that it avoids doubling the size of the state-space model, resulting in lower memory usage and faster runtime. Building on this, the current work proposes a Cholesky-based formulation of the \"fixed-point recursions\" method, maintaining the same efficiency in terms of memory usage and runtime.\nAll derivations for the proposed method (Algorithms 1 through 3) have been verified and found to be correct.\n\nThe numerical examples implemented using JAX effectively demonstrate and support the theoretical claims.\n\nReference:\n\nSimo Särkkä and Jouni Hartikainen. On Gaussian optimal smoothing of non-linear state space models. IEEE\nTransactions of Automatic Control, 55(8):1938–1941, 2010.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Major issues:\n1. After Equation 7, it is noted that these formulas do not have a Cholesky-based formulation. Does this imply that it is not possible to derive Cholesky-based formulations directly from these equations, rather than relying on Algorithms 1 and 3?\n\n2. Since this work focuses on linear discrete-time state-space models and claims that it can be further generalized to nonlinear models through techniques such as posterior linearization, it is essential to derive closed-form Cholesky-based formulations for all steps involved in finding $p(x_0 \\mid y_{1\n})$ using Algorithms 1 and 3, based on the state-space model in Equation 1 and the Cholesky form of the covariances in Equation 2. For instance, what is the formulation for $L_{P_{k-1 \\mid k}}$ in terms of the model parameters?\n\nMinor issues:\n1. In Equation 2, $b_k \\sim \\mathcal{N}(0, B_k)$ is given. What is the difference between B_k and Q_k which introduced above Equation 1?\n2. Above Equation 3; Kalman filter computes $p(x_k \\mid y_{1:k})$ for $k = 1, \\ldots, K$ not only $p(x_K \\mid y_{1:K})$.\n3. Feel free to use the correct name \"Särkkä\" instead of \"Sarkka\" for the reference Simo Särkkä and Jouni Hartikainen (2010).\n4. Regarding Equation 11, it is worth noting that $p(x_0 \\mid x_K, y_{1:K}) = p(x_0 \\mid x_K, y_{1:K -1} ) $, that is, given $x_K$, $x_0$ is independent from $y_k$, which we also see through rearranging the integrals.\n5. Step 2 of Algorithm 1: the iteration is from $k = 1$ to $k = K$. What does iterate from $k-1$ to $k$ mean, it seems to be extra.\n6. In Step 3 of Algorithm 1, there is a missing parenthesis in  $p(x_0 \\mid y_{1:K})$.\n7. Section 3.2, in the second line, $p(x_0 \\mid x_k , y_{1:k})$ should be $p(x_0 \\mid x_k , y_{1:k-1})$.\n8. In Algorithm 2 and 3, the notation $N$ is used instead of $K$.\n9. The subscript in Equations (16c) and (16d) should be $0|k$.\n10. If I understand correctly, all the methods in Tables 3 and 4 are Cholesky-based. If so, please indicate this in the tables.\n11. In the first experiment, what is the justification for using only three runs? Is that sufficient? Additionally, why not use the average runtime instead of selecting the best one?", "broader_impact_concerns": "There are no substantial concerns that necessitate a broader impact discussion.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "mtsV8ccnL1", "reviewer": "Reviewer_TN7Q", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This work is well written and tackels an important problem. The contributions listed in the manuscript seem at first glance correct and meaningful. I am by no means an expert in this literature, so I may not be the best person to assess their relevance within it. This is my main concern, that the derivations in this work, although correct, may be either existing in the literature, or very derivative. That being said, the paper makes a good point in convincing the reader their results are correct. One of my main concerns is the presentation of related and historical work. The introduction is extremely short and devles into technical details almost immediately, this is an issue in my opinion because it does not give the reader any bird's eye view of the field.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work introduces a new approach for fixed point smoothing that generalizes existing techniques such as the Racuh-Tung-Striebel smoother, filtering via an augmented state, and fixed point recursions. This is achieved by producing a new recursive derivation for fixed point smoothing, and carefully deriving update rules that can avoid an eugemnted state space computation. This work concludes with a variety of experiments designed to show the computational advantages of this new methodology.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Following up the above comments I would like to see the introduction improved and contain a richer set of related works and background for the reader.", "broader_impact_concerns": "No broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "gTyGRKhp5d", "reviewer": "Reviewer_UyAV", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n\n* simple method that solves a well-defined problem neatly\n* beautiful explanation of the filter recursions, including square root form, which is surprisingly rare in the literature\n* simple but challenging example problems\n\n# weaknesses\n\n* occasionally unclear text (see requested changes)\n* the authors miss one easy \"win\": as they point out, their innovation allows arbitrary Gaussian parameterisation of the fixed point recursion. If their goal to to minmise cost *and* stability, they should probably characterise the efficiency of the information parameterisation also, since that parameterisation can lead to efficient results some times.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors, by joggling the order of marginalisation and factorisation in filtering equations, devise an algorithm which appears to shift the pareto frontier of efficiency and stability for the problem of identifying starting conditions in a filter problem", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Minor clarity issue: sec 2.2, the intro paragraph refers to two things which don't match exactly the following two headings:\n\n>The state of the art in fixed-point smoothing does not use fixed-interval smoothers to compute p(x0 | x1:K ).\n>There are two more efficient approaches: Either we solve the fixed-point smoothing problem by computing\n>the solution to a specific high-dimensional filtering problem (Biswas and Mahalanabis, 1972), or we derive\n>recursions for how p(x0 | y1:k) evolves with k = 1, ..., K (Meditch, 1967b;a; 1969). \n\nWould that be better rewritten as the following?\n\n>The state of the art in fixed-point smoothing does not use fixed-interval smoothers to compute p(x0 | x1:K ).\n>There are two more efficient approaches: Either we solve the fixed-point smoothing problem by computing\n>the solution to a specific high-dimensional filtering problem (_state-augmented filtering_, Biswas and Mahalanabis, 1972), or we derive\n>_fixed-point recursions_ for how p(x0 | y1:k) evolves with k = 1, ..., K (Meditch, 1967b;a; 1969).\n\nOr have I misunderstood?\n\nThere is another weird problem with flow in section 2.2.\n\n>**Fixed-point recursions** Suppose we temporarily ignore Cholesky-based forms and commit to covariance-\n>based parametrisations of Gaussian distributions. In that case, we can improve the efficiency of the fixed-point\n>smoother\n\nI think you can delete the bits about CHolesky because in this point in the paper, while w have described THAT the cholesky decompositon can be used, we haven't actually seen it in action; everything so far has been in terms of covariances.\n\n>**Fixed-point recursions** We can improve the efficiency of the fixed-point\n>smoother in terms of covariances\n\n7a: can you clarify? set k=1. Then we have\n\n$$\n\\begin{aligned}\nG_{0|1}&=G_{0|0}G_{0|1}\\\\\\\\\n&=G_{0|1}\\\\\\\\\n\\end{aligned}\n$$\n\nThis looks under-specified; Anything times the identity is itself. What is going wrong here? Have I introduced an off-by-one error? Should $G$ be coupled with $C$? Something else?", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers all agreed that the claims were well-supported and there was a sufficient audience for the work. During the discussion period, the paper was revised by the authors both for clarity and to motivate the need for this work. As there are no major outstanding concerns from the reviewers, I recommend acceptance as is.
Decision: Accept (poster)