The Rank-Reduced Kalman Filter: Approximate Dynamical-Low-Rank Filtering In High Dimensions
Authors: Jonathan Schmidt
Philipp Hennig
Jörg Nick
Filip Tronarp
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian
filtering
smoothing
bayesian
state-space models
dynamic-low-rank
high-dimensional
spatio-temporal
Gaussian processes
regression
low rank
state estimation
TL;DR: An approximate Kalman filtering (and smoothing) algorithm is presented that deterministically tracks low-rank representations of the covariance matrices. This makes computations in high-dimensional state-space models tractable.

Abstract: Inference and simulation in the context of high-dimensional dynamical systems remain computationally challenging problems.
Some form of dimensionality reduction is required to make the problem tractable in general.
In this paper, we propose a novel approximate Gaussian filtering and smoothing method
which propagates low-rank approximations of the covariance matrices.
This is accomplished by projecting the Lyapunov equations associated with the prediction step to a manifold of low-rank matrices,
which are then solved by a recently developed, numerically stable, dynamical low-rank integrator.
Meanwhile, the update steps are made tractable by noting that the covariance update only transforms the column space of the covariance matrix, which is low-rank by construction.
The algorithm differentiates itself from existing ensemble-based approaches in that
the low-rank approximations of the covariance matrices are deterministic, rather than stochastic.
Crucially, this enables the method to reproduce the exact Kalman filter as the low-rank dimension approaches the true dimensionality of the problem.
Our method reduces computational complexity from cubic (for the Kalman filter) to quadratic in the state-space size in the worst-case, and can achieve linear complexity if the state-space model satisfies certain criteria.
Through a set of experiments in classical data-assimilation and spatio-temporal regression, we show that the proposed method consistently outperforms the ensemble-based methods in terms of error in the mean and covariance with respect to the exact Kalman filter. This comes at no additional cost in terms of asymptotic computational complexity.

The Rank-Reduced Kalman Filter: Approximate Dynamical-Low-Ram Filtering In High Dimensions
Jonathan Schmidt
Philipp Hennig
Tübingen Al Center, University of Tübingen jonathan schmidt@uni-tiuee tuebingen..nnnn .ee
Tübingen Al Center, University of Tübingen
philipp.. -tuebiigeetueueennnnn de
Jörg Nick
Filip Tronarp
University of Tübingen nickêna, uni meii-tuebingen.de
Lund University wrw... tww.arrnnttttttt.t...... lu. lu.se
Abstract
Inference and simulation in the context of nigh-dimensional dynamical systems remain computationally challenging problems. Some form of dimensionality re- duction is required to make the problem tractable in general. In this paper, we propose novel approximate Gaussian filtering and smoothing method which prop- agates low-rank approximations of the covariance matrices. This is accomplished by projecting the Lyapunov equations associated with the prediction step to a manifold of low-rank matrices, which are then solved by recently developed, numerically stable, dynamical low-rank integrator. Meanwhile, the update steps are made tractable by noting that the covariance update only transforms the column space of the covariance matrix, which low-rank by construction. The algorithm differentiates itself from existing ensemble-based approaches in that the low-rank approximations of the covariance matrices are deterministic, rather than stochastic. Crucially, this enables the method reproduce the exact Kalman filter the low- rank dimension approaches the true dimensionality of the problem. Our method reduces computational complexity from cubic (for the Kalman filter) quadratic in the state-space size in the worst-case, and can achieve linear complexity if the state-space model satisfies certain criteria. Through cal data-assimilation and spatio-temporal regression, we show that the proposed method consistently outperforms the ensemble-based methods in terms of error in the mean and covariance with respect to the exact Kalman filter. This comes at no additional cost in terms of asymptotic computational complexity.
of experiments in classi-
1 Introduction
Spatio-temporal dynamical systems have always played an important role in the applied sciences, such as climate science, numerical weather prediction, geophysics. In precisely these settings, one is often faced with massive amounts of data to be processed. At the same time the interactions are of such complex nature that it imperative include a notion of uncertainty in the model and in its outputs. Both these requirements have a reputation as being computationally expensive. Hence, sensible approximations are indispensable certain scale.
Concretely, we consider state-space models of the form
dx(t) Ax(t) dt ddw(t), y(1)) x(1)) ~ ()),,)), 1= 1,,,..,,
(la) (1b)
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Kalman filter
ours: T=3
11
Figure 1: Rank-r approximations to the true KF covariance (left) for increasing r. The considered problem's 's true rank r* by construction, the state-dimension is n = 1000. The respective low-rank factors are shown above and left of their outer product. For 2 r", the KF estimate is recovered. For r > r", the excess columns the low-rank factor collapse (rightmost plot).
where x E 1 is the latent state, y Rn is the measurement process, and standard Wiener process in R9, Equation (Ia) defines the prior dynamics model via a linear time-invariant (LTI) stochastic differential equation (SDE) with drift matrix E R and dispersion matrix B E Rxx The linear Gaussian observation model Eq. (1b) is specified via the measurement matrix C E mR and measurement-noise covariance E mx where possible dime-dependency has been omitted from the notation. While the probabilistic state-estimation problem may in principle be solved in closed form using the Kalman filter (KF) and smoother [25, 34), we will consider the case where the state dimension n and measurement dimension are very large, making the covariance recursion computationally prohibitive in practice. Such settings include spatio-temporal Gaussian process regression [35, 37, 45, 7, 20], and data assimilation applications such as metcorology and oceanography [17, 11], numerical weather forecasting [21], geoscience [6], inverse problems [9, 41], and brain imaging [16]. While the Gaussian process (GP) literature offers many ways of tractable, approximate inference, approaches based on approximations of the underlying kernel function are not amenable the present problem setting. Without a kernel to approximate, the problem is often solved by ensemble methods, which stochastically propagate low-rank approximations to the covariances in the filtering recursion. This stochasticity introduces unfavorable properties in these methods. To avert this shortcoming, we propose a novel low-rank filtering recursion, which fully deterministic.
Contributions This text develops method for efficient. approximate Gaussian filtering, including smoothing and marginal likelihood computation. 'he method employs the standard predict/comrect formulation of the filtering problem, and can principle be thought of a Kalman filter, in that it represents all covariance matrices as their left square-root factors [19]. However, instead of square matrix square roots, the covariances are approximated by (n x r)- dimensional low-rank factors, where << We propose approximate low-rank prediction- correction loop, together with the computation of the backwards Markov process representing the smoothing posterior, all while preserving this low-rank structure throughout time.
square-root implementation
Crucial to this work are (i) the use of dynamic-low-rann approximation (DLRA) [27, 29, 8] solve Lyapunov equations on the manifold of rank-r matrices in the prediction step and (ii) an efficient correction scheme that updates only the low-dimensional column space of the predicted covariance factor, which introduces no additional approximation error. The computational complexity of the proposed filter is asymptotically equal to those of the widely-used ensemble Kalman filters.
The method offers a deterministic alternative to existing ensemble methods. In particular, it recovers the exact Kalman filter estimate both in the full-rank limit
and if the problem is in fact of rank r < n. The latter is demonstrated in Fig. 1, in which a problem with n 1000 state dimensions and a true rank of 7 approximated with r 3, 7, 11. As soon as exceeds the true rank of the problem, the residual subspace dimensions collapse. The method will be detailed in Section 3 and evaluated on the basis of various experiments in Section All source code is publicly available on GitHub."
-j1
hhttps://githu:

---

Related work The state estimation problem defined by Eq. (1) amounts to linear Gaussian state estimation and can be computed in closed form by Kalman filtering and smoothing. A modern overview can be found in Särkkä and Solin [36]. While the exact filtering and smoothing recursions scale linearly with the number of time points N, they scale cubically in the state dimension n and the measurement dimension m, making their use computationally infeasible in many high-dimensional problems. A widely-used approximate filtering technique used to circumvent this computational burden are ensemble Kalman filters (EnKF) [11-13, 26, 32, 6J. The ensemble is set of randomly sampled state vectors of dimension n, which parametrize the mean and covariance via sample statistics. It is processed successively in centered ensemble can be regarded as a low-rank square-roo factor of the sample covariance matrix. In Section 4 we compare our method to the standard version of the EnKF and the ensemble transform Kalman filter [4], which is categorized as an ensemble-square-rod filter [43].
loop similar to the KF. The mean-
In so far as it aims for computationally tractable Gaussian process regression, the proposed method is related to approaches using the Nyström method [48], inducing point methods [38, 44, 22], random Fourier features [31], and doubly-sparse variational Gaussian processes [1].
The literature offers more contributions related the present work under the broader context of low- rank and dimensionality--educiiii methods. Loiseau et al. [28], Netto and Mili [30], Vijayshankar al. [47] approach low-rank modelling of dynamical systems based on measured data, whereas Farrell and Ioannou [15] are concerned with low-rank approximations of a given model based on classical theory on linear inmeiinvrriant systems. Cressie et al. [10] develop a tractable hicrarchical Bayesian model to approximate nigh-dimensional spatio-temporal dynamics in reduced latent space and condition on measurements by transporting the low-dimensional diffusive process to the full space by a projection that is assumed to be known. In those cases, once a low rank model is obtained for the phenomena under study, state estimation becomes computationally tractable, due to the reduced size of the resulting model. These works build upon the assumption that the state of the system evolves approximately in a low-dimensional manifold, whereas our method works under the assumption that the error in the state estimate evolves approximately in low-dimensional manifold.
2 Background
Throughout this work, for function x(t) and grid (f1....ty) .y) we use the abbreviated notation ar(1)) = xy and (f)) = Ta The task of interest is computing the filtering densities p(a1 yıu) and smoothing densities p(az y1.N) for forl=1 1,..., N. LTT SDEs of the form in Eq. (la) admit an equivalent discrete formulation in terms of linear Gaussian transition densities of the form
x ~ ((hxxi1, Q!).
The transition matrix l(t) E R2x= and process-noise covariance matrix Q(t) E R solve the matrix differential equations (()) = (t) and
((t) AQ(() Q(t)A* + BB",
(3)
in the interval ()) with initial conditions (t:-1) I and Q(ti-1) 0. respectively [36]. Equation (3) is known as a Lyapunov equation, and will be of special importance in this work.
2.1 Square-root filtering
For numerically stable filtering and smoothing, all covariance matrices involved can be represented by matrix square roots (e.g. Cholesky factors). This is known square-root filtering/mmoothing [19]. Let (M1 M2) denote wide block matrix in R4 (d2+ds), where M1 R xdl and M2 Rd xds
are some matrices. Given the filtering covariance at time 4111 as Σn = 3// 2//2 a square-root
factorization of the predicted covariance II, = 11 is given by
a. (DEE/) )
r -1 16 + Qi
and, hence, ( / carried out entirely square-root factors. For more details, see Grewal and Andrews [19].
22 is a matrix square-root of πT. Similarly, the correction step can be
3

---

2.2 Dynamical low-rank approximation
Consider an initial value problem with matrix-valuet flow field F R Rux"" + R
M(1) = F(t,)(())
M(lo) Mo.
(5)
where and potentially very large. Dynamic low-rank approximation (DLRA) methods [27, 8] efficiently compute low-rank factorizations Y(t) = U(t)D(t)V R M(t) of Eq. (5) by solving
(t)) P,(Y(+))
((t),
Y(to) U.D.Vo
(6)
instead. The matrices Ug E RRux, "xxr,,, Do E R' and Vo E R are an initial low-rank factorization of (+o) M(to)) (t)] denotes the projection operator onto the tangent space at Ytt), where (t) lies the manifold of rank-r matrices. the present work, we leverage this technique to obtain a low-rank factor Q of the process-noise covariance at the next prediction location. Thus, we avoid much of the computational cost in the prediction step, otherwise caused by solving the full Lyapunov equation for Q in Eq. (3), for example using matrix- fracii-f decomposition [40, 3]. Concretely, let F(Q) AQ QA + BB'. Then, an approximate low-rank process- noise covariance matrix associated with the prediction step is obtained by integrating
(tt)) P.[Q(t)) F(Q())),
Q(t-1)) UOD =0
from f to t.. At the start the filtering recursion to. an initial low-rank factorization is constructed with random orthogonal matrix Uo and Do O. For all subsequent steps the propagated orthogonal basis U(tz) can be reused. In this work, all mentions of dynamic low-rank integration always refer to the recently developed, numerically stable basis update & Galerkin (BUG) integrator [8], whose error bounds are independent of small singular values. More details are given in Appendix C.
Rank-reduced Kalman filtering
In this section, method for approximate inference in Eq. (1) is developed. The idea is to approximate the full filteringsmoothing covariance matrices of the latent state by an eigendecomposition truncated at the r-th largest eigenvalue. The main challenge is thereby attain linear computational scaling in n and m under appropriate assumptions on the state-spact model. Before proceeding with the filtering and smoothing recursions, it is instructive to examine inference in static low-rank model.
3.1 Efficient inference in static models of low rank
Consider the following latent variable model
x~ (u, (I/"),
y N(C,,,
(8)
where 1 E RRxr of full column rank, r S n, and R E Rmxm Representing the latent state x as
x= = 1222,,
z ~~ Ix.,
reduces the problem to inference the following model
~~ (0,
y v| CIT'/2,, R).
(10)
An efficient inference scheme given by the following proposition, which is proved in Appendix A.I. Proposition 1. Let and y be two random variables governed by Eq. (10). Assume r S m and consider the following singular value decomposition (RT//CITV/3 UDV*, where U,D Rixr and mm Then, defining the whitened residual = R- 2/- Cμ), we get
y ~ NCCH,, R (VD2V* I)R²
(1la)
e, U(T + D²)- x y ~ N(U(( DV"
(11b)
Furthermore, let |-| denote the matrix determinant, The marginal log-likeldhood of given by
ly; CH, R' +)R* = m
log 2 - -log IR
log(1 + D2.))
2
1=
(12)
1
/|o122 + 2 DD
2

---

The below corollary follows from the deterministic relationship between z and x given by Eq. (9). Corollary 1. Let y and x two random variables governed by the model Eq. (8) and r S m. Then
x x1y~ ~ Nu + Ke, E),
(13)
where
k= (/)/UUTTT>---TTT
== UT
(14)
Moreover, a square-root of Σ is readily obtained by
L I D?) 1/2
(15)
3.2 The rank-reduced filtering recursion
In this section, a low-rank the filtering densities and the logarithm of the marginal likelihood.
correction recursion is developed for the purpose of obtaining
The prediction equations Suppose the filtering distribution at time 4-1 is given by
p Y (xxxx 11-1
(16)
where 1
| We begin by detailing how to compute low-rank factor 11, of the predicted covariance. First compute the truncated singular value decomposition (SVD)
1 O²') Ü.D.V
(17)
of the square-root factor (ef. Section 2.1), with Ür Raxr, Bj E R. and V, E Rexr, This can be done in an optimal way by computing the full SVD and truncating it at the r-th largest singular value in (nrr) [18]. The process-noise covariance factor Qj Rexr in Eq. (17) computed using DLRA as described in Section 2.2. The rank-reduced predicted moments at time fy are thus given as
Pi 4:4-1 E R.
(18)
u U,D,
ER
(19)
The update equations
The low-rank update equations follow from Proposition and Corollary 1:
jt 14 + Kjet E TR""
(20)
Σ/- = I1, U. D2) RR
(21)
with the whitened residual es = RR---/// 1/= - CHT) ) and
(R +²C U,D.Vi, K for Up Dr E Rex and V, E mxx The marginal predictive lig-liiel of y, is given by
(22)
logp(y: Y1a-1)= log 2 log R|//1
log (Dr) 11
2
(23)
1e
1 V.
- 2
2
Corollary and by extension, the above filtering recursion-is only valid when r S m, which is exactly the intended usecase for the proposed method. For settings in which < < the correction step follows the square-root correction of a Kalman filter and is detailed in Appendix B.
3.3 Time complexity of the rank-reduced filtering recursion
This section analyzes the computational complexity of the proposed method. In the worst case, the method scales quadratically in the state dimension n and measurement dimension m.. Under favorable conditions, which are met in many real-world applications, the rank-reduced Kalman filter obtains computational complexity of (nr²) as stated by Proposition 2, which proven Appendix A.2. We begin by specifying set of assumptions for Proposition 2.
5

---

Assumption 1.
The maps - Ax, x 4x, and x + BB'x can be evaluated in O
+
Assumption fulfilled naturally in many dynamical systems, in which local operators introduce sparsity into the dynamics, as, e.g., in dinite-difference approximations of spatial differential opera- tors. Another example Kronecker structure in the system matrices arising in spatio-tempone GP regression by assuming the covariance function is separated into product over the respective spatial and temporal components [35, 39, 20, 42].
Assumption 2.
The map x -> Cx can be evaluated in (()).
Assumption 3. The map X + R and the log- log-determinant log [R can be evaluated in (m).
Throughout this work, we refer the situation in which Assumption or 2 do not apply as the "worst case". Assumption 3 is taken for granted as the meesurement-nokee covariance R is often a diagonal matrix, which implies that sensor errors are uncorrelated. This is not only realistic but also commonly imposed in modelling. The above assumptions allow the following proposition.
Proposition 2. Given Assumptions to 3, the proposed method approximates the filtering densities and the marginal likelihood at a cost of O mr² r
The below corollary follows from the book-keeppin in the proof of Proposition 2 in Appendix A.2.
Corollary 2. When Assumption I or Assumption 2 are not satisfied, the worst-case complexity of the proposed low-rank filtering recursion O(nir + nm + m²r).
Overall, the proposed filtering recursions achieve the same asymptotic complexity as the existing ensemble methods under the same set of assumptions. However, our method is entirely deterministic.
3.4 The rank-reduced smoothing recursion
It remains to obtain a recursion for the smoothing densities, which can be shown to be computationally tractable given that the filtering recursion tractable. The smoothing densities are denoted by
p(x yI:N) (x Ar), EN JUN, AN Σn.
(24)
It can be shown that the posterior process has backward Markov representation [5], hence the smoothing marginals may be obtained by
pxx VI.N)
11,+1(xx) x4+1)p(x411 \1:N) dx
(25)
where we call br.4 i the backwards kernel. Consequently, the problem consists of approximating br,411 such that the marginalization Eq. (25) may implemented in a computationally frugal manner.
Approximating the backwards kernel For the linear Gaussian case, the backward kernel [46]
b x1+1) = (xx;; G/81+1 + Pr),
(26)
is parametrized by the smoothing gain G E R=*x,, the shift vector E R², and the covariance of the backwards kernel Py gaxn Let (-) denote the Moore-Penrose pseudoinverse. Then
G = E
Vi - = P. = (L=
G (27)
In the low-rank setting, the backwards kernel is efficiently approximated given the results from above. The approximate smoothing gain is product of a tall, a small quadratic, and a wide matrix:
G Rd 2
Pitz
(u//1)
(1
(28)
:T.ERex*
During filtering. T3 saved alongside the low-rank factors of the prediction and filtering covariances. We proceed to compute a low-rank representation of the tookwards-transideee covariance. Consider the following singular value decomposition, truncated at the r-th largest singular value.
(a G.4 G O,D,VT,
(29)

---

A.
vs. KF menn
vs. KF CoV.
A.
vS. KF mean
B. vs. KF onv.
Frobenius distanie
1.0
Frobenius distance
RRKF
10
RMSE
RMSE
(ours)
0.5
0.5
EnKF
ETKF
0
00000
0.0
00000
0.0
51
100
144
low-rank dimension
low-rank dimension
(a) Linear advection dynamics (true rank = 51).
(b) London air-quality regression.
Figure 2: Performance of the low-rank filters different setlings. In the truly low-rank linear advection problem (a), the RRKF achieves the optimal estimate (+) when r exceeds the true rank of the problem. For the spatio-temporal air-quality regression (h), our method is consistently better.
where Üj E RR Dr e grxr diagonal, V 12*r, and G₁ as given in Eq. (28). Then rank-r factor of P/ is given as
p E IR"
(30)
As for filtering, in the worst case, the cost for low-rank smoothing scales quadratically with the state dimension n. The following proposition, which is proved in Appendix A.3, states that this can be further reduced to linear complexity in n, given the assumptions on the dynamics model are satisfied. Proposition 3. Given Assumption I, computing the smoothing density p(x, I VI:N) costs Onr +
Finally, it is useful mention that realizations of the backwards process posterior samples.
4 Experiments
This section evaluates the proposed method in different experimental settings. The chosen measure of quality is the distance of the approximate low-rank moments to the exact KF. We measure the mean deviations with the cet--mean-sqqaree error (RMSE) and the covariance deviations with the time-averaged relative Frobenius distances. The presented rank-reduced Kalman filter (RRKF) is compared to the ensemble Kalman filter (EnKF) and the ensemble transform Kalman filter (ETKF). All EnKF and ETKF results are given in sample statisties over 20 runs. Section 4.1 shows that for a truly low-rank system, our method recovers the KF estimate up to numerical error as soon as exceeds the true rank of the problem. Section 4.2 tests the method problem with real data. After Section 4.3 evaluates the approximation quality for increasingly low- rank systems in controlled experimental environment, Section 4.4 verifies the stated asymptotic cost of the method. Finally, a large-scale pratio-tempors regression problem is solved in Section 4.5, In all experiments, we compute the stationary mean and low-rank factorization of the stationary covariance matrix of the prior and condition the stationary moments on the first measurement in the respective time-series dataset. More details on the experimental setups are given in Appendix D.
spatio-temporal GP regression
4.1 Linear advection model
The proposed algorithm is evaluated in a standard asstmmssiiil setup in which linear-advcction dynamies with periodic boundary conditions are assimilated to set of simulated data. The setup follows the description in Sakov and Oke [33] and the experimental details are additionally detailed in Appendix D. This problem has rank 51 by construction. Figure 2a shows how the exact KF estimate is recovered by our method for r = 51, while the ensemble methods converge according to Monte-Carlo rate. This is particularly clear for the covariance estimate (Fig. 2a, B.).
4.2 London air-quality regression
This experiment uses hourly data from the London air-quality network [23] between January 2019 and April 2019, which amounts measurements at 72 spatial locations and 2159 points in time. The
7

---

RRKF (ours)
EnKF
ETKF
cumulative spectrum
0.01
t= 0.1
J. 0.25
8 = 1.0
RMSE
1.0
1.0
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
Frobenius
distance
1.0
1.0
1.0
000000
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
111
221
331
441
111
221
331
441
111
221
331
441
111
221
331 441
low-rank dim.
low-rank dim..
low-rank dim..
low-rank dim.
How does the rate of spectral decay influence the low-rank approximation? The larger the spatial length scale of a spatio-temporal Matérn model, the faster the spectrum decays and the faster all methods converge to the KF estimate. The first and second row compare the mean estimates and the covariance estimates, respectively. The RRKF estimate is consistently closer and recovers the KF estimate (0) r n. The cumulative spectrum of the final-step KF covariance is shown in grey.
Figure 3:
A.
Assumptions to
B.
Assumptions 2 &
10"
RRKF
10
(n)
wal time
10
(ours)
wall time
10a
O(n)
((n)
10"
10
O(n))
(nn))
10"1
EnKF
10
O(n")
1
ETKF
10
10a
10
10
10¹
10
10a
10*
state dimension 11
state dimension n
Figure Best-case and worst-case asymptotic complexities of the low-rank filters. When all Assumptions 1 to 3 are satisfied (A.), the elapsed time grows linearly with the state dimension If parts of the assumptions are not met (B.), the cost scales quadratically with the respective dimension.
data used, together with its processing (except the train-test split), is the same as the corresponding experiment by Hamelijnck et al. [20]. The model is spatio-temporal Matérn-32 process with prior hyperparameters that maximize marginal log likelihood. Figure 2b shows that compardd the ensemble methods-ov algorithm is consistently closer the KF estimate and recovers it
4.3 Spatio-temporal Matérn process with varying spatial lengthscale
Consider spatio-temporal Gaussian process x(1) ~ GP ke @ kn) with covariance structure that separable time and space. Let both k; and k2 be of the Matérn family with characteristic spatial and temporal lengthscales l,, f,r and output scales of. o.. Spatio-temporal Matérn processes can be translated to the formulation in Eq. (Ia), as detailed, for instance, by Solin [39]. This setting allows for varying the "low-rankness" by the choice of l,: the larger la, the more information about spatial point is prescribed by its surroundings. As Eg approaches zero, the spatial kernel matrix gets more and more diagonal, encoding spatially independent states, and causing barely decaying spectrum. The spectrum for large &, decays rapidly, making low-rank approximations more accurate for small r.
In the experiment, we consider a spatial domain [0.2] x [0. C IR?, which is subsampled at uniformly-spacee grid with A, = 0.1 (n = trajectory. are drawn from realization of the prior process. Figure 3 demonstrates that, for increasing spatial lengthscales, methods converge faster to true KF estimate. Crucially, however, our method is consistently closer to the optimal estimate and recovers it for
21 = 441). Noisy observations of the full state

---

30 Jan 2023
19 Feb 2023
04 Mar 2023
50
Dota
40
25
20
80
60
Meaan
25
0
Z-score
-0
0
-3
-3
Figure 5: Rainfall Australia. A spatio-temporal GP regression problem with n 911 state dimensions, solved with the RRKF with 1000. The columns are three different days in the time series. The rows show the data, mean estimate, and Z-scores, respectively.
4.4 Runtime
To demonstrate the asymptotic computational complexities given in Section 3.3, investigate the runtimes on two different settings: The first fulfills both Assumptions and 2, whereas the second only fulfills Assumption 2, which means that the cost of prediction is quadratic in the state dimension. Both problems are solved on a temporal grid of size N = 100 with the measurement dimension fixed m 100 and the low-rank dimension fixed
The problem from Section 4.1 serves as the best-case setting. The noise-free Linear-advection dynamics amount multiplication of the state with circulant matrix, which is implemented efficiently using fast Fourier transform. The measurement operator amounts to array indexing. As subplot A. in Fig. 4 confirms, the computational cost grows linearly with the state dimension. As worst-case example we employ a spatio-temporal Matérn process with dense spatial kernel matrix (cf. Section 4.3). In this case, expect the computation time scale quadratically with the state dimension because BB* is dense. This is verified by the right plot in Fig. 4 (B.).
4.5 Large-scale spatio-temporal GP regression on rainfall data
As final experiment, a spatio-temporal GP regression problem is solved on a large gridded analysis data measuring rainfall in Australia. The data is provided by the Australian Water Availability Project (AWAP) [2, 24]. The number spatial data points is n = 3001,, which rules our the cubic-in-n Kalman filter. The smoothing posterior computed at January, 2023 through 5 March, 2023. The total number of observations is thus n N = 236 440. The spatio-tempores prior Matérn model is selected by maximizing the marginal log-likelihood of lower-resolution data set. A set of 6160 spatial points is excluded from the data during filtering and smoothing in order to evaluate spatial interpolation performance. Thus, the measurement dimension is m 6160 = 24751. Figure 5 shows the low-rank approximation to the GP posterior with == 1000 low-rank dimensions. The non-observed spatial locations are framed the horizontal-"S" shape. At measurement locations, the data is described well by the model, while slight smoothing- out effect can observed due to the high-frequency features of the model being truncated. At evaluation points the interpolation aligns with the expectations regarding the simple Matérn model.
40 time points, from 25

---

The Z-score map shows x where & and X denote the smoothing mean and marginal standard deviation and the data point, at time &t. This highlights badly calibrated uncertainty estimates, wherever the model fnder//oversssimaes the data and is divided by small standard deviation, in blue/red. The Z-score distribution receives further attention in Appendix E. All in all, the RRKF achieves high-quality approximate estimates while compressing the state by a factor of 30911 30.
1000
5
Limitations
The main limitation of the proposed approximate probabilistic inference scheme is that the truncation dimensions cuts away covariance information, which is not accounted for. As a result, confidence grows as r shrinks, while estimates should arguably get more uncertain instead with increased compression. In existing low-rank filters, this issue is often counteracted manually by inflating the covariance matrices, on a per-application basis [6, Section 4.4] and calls for a more principled treatment. Finding way to keep track of that residual uncertainty information is beyond the scope of this paper. It remains unclear how to preserve the stated asymptotic complexities while doing so.
Further, all results herein are stated for linear dynamics and observation models. Extensions to non-linear models could include linearization of the according transitions, or cubature methods [34].
Conclusion
By building upon well-established knowledge about optimal compression, we have proposed an algorithm providing principled way to balance approximation accuracy against computational cost in high-dimensioonl state estimation. It combines simple (truncated) singular value decompositions with recent numerical low-rank integrators of large matrix differential equations. We have offered both theoretical and empirical arguments for why it is desirable to use deterministic algorithm when approximating large-scale probabilistic satteesttimttion problems.
Acknowledgments and Disclosure of Funding
The authors gratefully acknowledge financial support by the German Federal Ministry of Education and Research (BMBF) through Project ADIMEM (FKZ 01IS18052B), and financial support by the European Research Council through ERC StG Action 757275 / PANAMA; the DFG Cluster of Excellence "Machine Learning New Perspectives for Science", EXC 2064/1, project number 390727645; the German Federal Ministry of Education and Research (BMBF) through the Tübingen AI Center (FKZ: 011S 10399A;; and funds from the Ministry of Science Research and Arts of the State of Bartt-Wüurteebberg. Jörg Nick is funded by the Deutsche Forschungseemeensccatt (DFG, German Research Foundation) Project-ID 258734477 - SFB 1 1173. Filip Tronarp was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation. The authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Jonathan Schmidt. The authors also thank Nathanael Bosch and Christian Lubich for many valuable discussions and for helpful feedback.
References
[1] Adam, V., Eleftheriadis, S., Artemev, A., Durrande, N., and Hensman, (2020). Doubly sparse variational Gaussian processes. In Chiappa, S. and Calandra, R., editors, Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 2874-2884. PMLR.
[2] Australian Data Archive for Meteorology (ADAM) (2023). Australian Water Availability Project. jsp. [Online; accessed May 2023].
http://www. wwwpwwwwwwwwwwpww.www........... .ow.ggov... www.ggoo..aaaaaaa
[3] Axelsson, P. and Gustafsson, F. (2015). Discrete-time solutions to the continuous differential Lyapunov equation with applications Kalman filtering. IEEE Transactions on Automatic Control, 60(3):632-643.
[4] Bishop, C. H., Etherton, B. J., and Majumdar, S. (2001). Adaptive sampling with the ensemble transform Kalman filter. part i: Theoretical aspects. Monthly weather review, 129(3):420-455
10

---

[5] Cappé, O., Moulines, E., and Ryden, T. (2005). Inference in Hidden Markov Models. Springer Series in Statistics. Springer.
[6] Carrassi, A., Bocquet, M., Bertino, L., and Evensen, G. (2018). Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinarry Reviews: Climate Change, 9(5).
[7] Carron, A., Todescato, M., Carli, R., Schenato, L., and Pillonetto, G. (2016). Machine learning meets Kalman filtering. In 2016 IEEE 55th conference decision and control (CDC), pages IEEE.
[8] Ceruti, G. and Lubich, C. (2022). An unconventional robust integrator for dynamical low-rank approximation. BIT. Numerical Mathematics, 62(1):23-44.
[9] Chada, N. K., Stuart, A. M., and Tong, X. T. (2020). Tikhonov regularization within ensemble Kalman inversion. SIAM Journal on Numerical Analysis, 58(2): 1294.
[10] Cressie, N., Shi, T., and Kang, E. L. (2010). Fixed Rank Filtering for Temporl Data. Journal of Computational and Graphical Statistics, 19(3):24--745.
[1I] Evensen, G. (1994). Sequential data assimilation with a nonlinear quasi-geosttophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research: Oceans, 99(C5):10143-00622
[12] Evensen, G. (2009a). Data assimilation: the ensemble Kalman filter, volume 2. Springer.
[13] Evensen, G. (2009b). The ensemble Kalman filter for combined state and parameter estimation. IEEE Control Systems Magazine, 29(3)833-104
[14] Evensen, G. (2021). Formulating the history matching problem with consistent error statistics. Computational Geosciences, 25(3).
[15] Farrell, B. F. and Ioannou, P. J. (2001). State estimation using reduced-orde kalman filter. Journal of the Almospheric Sciences, 58(23):3666-3680)
[16] Galka, A., Yamashita, O., Ozaki, T., Biscay, R., and Valdés-Sosa Sosa,, P. (2004). A solution to the dynamical inverse problem of EEG generation using spatiotemporal Kalman filtering, Neurolmage, 23(2):435-44
[17] Ghil, M. and Malznotte-Rizzolii P. (1991). Data assimilation in meteorology and oceanography. In Advances in geophysics, volume 33, pages 141-266 Elsevier.
[18] Golub, G. H. and Van Loan, C. F. (2013). Matrix computations. JHU press.
[19] Grewal, M. and Andrews, A. (2011). Kalman Filtering: Theory and Practice Using MATLAB. Wiley.
[20] Hamelijnck, O., Wilkinson, W., Loppi, N., Solin, A., and Damoulas, T. (2021). Spatio-temporal variational Gaussian processes, Advances in Neural Information Processing Systems, 34:23621- 23633.
[21] Hamill, T. M. and Snyder, C. (2000). A hybrid ensemble Kalman filter-3D variational analysis scheme. Monthly Weather Review, 128(8):2905-29)
[22] Hensman, J., Fusi, N., and Lawrence, N. D. (2013). Gaussian processes for big data. In Proceedings of the 29th Conference Uncertainty in Artificial Intelligence (UAI), pages 282-290.
[23] Imperial College London (2020). Londonair London air quality network (LAQN). https:
/ //www. wwww..ooodd londonaii. org uk.
[24] Jones, D., Wang, W., and Fawcett, R. (2009). High-quality spatial climate data-sets for Australia. Australian Meteorological and Oceanographic Journal, (8(00)):33322) -288.
[25] Kalman, F. (1960). new approach to linear filtering and prediction problems. Transactions of the ASME-Journal of Basic Engineering, 82(Series D):35-45.

---

[26] Katzfuss, M., Stroud, J. R., and Wikle, C. K. (2016). Understanding the ensemble Kalman filter. The American Statistician, 70(4):350-3577
[27] Koch, O. and Lubich, C. (2007). Dynamical low-rank approximation. SIAM Journal on Matrix Analysis and Applications, 29(2):434-4544
[28] Loiseau, J.-C., Noack, B. R., and Brunton, S. L. (2018). Sparse reduced-order modelling: sensor-base based dynamics to full-state estimation. Journal of Fluid Mechanics, 844:459-490.
[29] Lubich, C. and Oseledets, I. V. (2014). A projector-spliting integrator for dynamical low-rank approximation. BIT Numerical Mathematics, 54(1):171-188.
[30] Netto, M. and Mili, (2018). A robust data-driven koopman kalman filter for power systems dynamic state estimation. IEEE Transactions on Power Systems, 33(6):722-7227.
[31] Rahimi, and Recht, (2007). Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems (NeurIPS).
[32] Roth, M., Hendeby, G., Fritsche, C., and Gustafsson, F. (2017). The ensemble Kalman filter: a signal processing perspective. EURASIP Journal Advances in Signal Processing, 2017(1):1-16.
[33] Sakov, P. and Oke, P. R. (2008). A deterministic formulation the ensemble Kalman filter: An alternative to ensemble square root filters. Tellus A: Dynamic Meteorology and Oceanography, 60(2):361.
[34] Särkkä, S. (2013). Bayesian Filtering and Smoothing. Cambridge University Press.
[35] Särkkä, S. and Hartikainen, J. (2012). Infinite-dimensional Kalman filtering approach to spatio- temporal Gaussian process regression. In Artificial Intelligence and Statistics, pages 993-1001. PMLR.
[36] Särkkä, S. and Solin, A. (2019). Applied Stochastic Differential Equations. University Press.
Cambridge
[37] Särkkä, S., Solin, A., and Hartikainen, J. (2013). Spatiotemporal learning via infinite- dimensional Bayesian filtering and smoothing: A look at Gaussian process regression through Kalman filtering. JEEE Signal Processing Magazine, 30(4):51-61.
[38] Snelson, E. and Ghahramani, Z. (2005). Sparse Gaussian processes using pseudo-inputs. Advances in Neural Information Processing Systems (NeurIPS), 18:125-112
[39] Solin, A. (2016). Stochastic Differential Equation Methods for Process Regression. Doctoral thesis, Aalto University, School of Science.
Gaussian
[40] Stengel, R. (1994). Optimal Control and Estimation. Dover Publications.
[41] Stuart, A. M. (2010). Inverse problems: Bayesian perspective. Acta numerica, 19:451-559.
[42] Tebbutt, W., Solin, A., and Turner, R. E. (2021). Combining preudo-point and state space approximations for sum-separable Gaussian processes. In de Campos, C. and Maathuis, M. H., editors, Proceedings of the Thirty-Seventl Conference on Uncertainty in Artificial Intelligence, volume 161 of Proceedings of Machine Learning Research, pages 1607-1617. PMLR.
[43] Tippett, M. K., Anderson, J. L., Bishop, C. H., Hamill, T. M., and Whitaker, J. S. (2003). Ensemble square root filters. Monthly weather review, 131(7):145--449)
[44] Titsias, M. (2009). Variational learning of inducing variables in sparse Gaussian processes, In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS). pages 567-574. PMLR.
[45] Todescato, M., Carron, A., Carli, R., Pillonetto, G., and Schenato, L. (2020). Efficient spatio- temporal Gaussian regression via Kalman filtering, Automatica, 118:10
12

---

[46] Tronarp, F.. Bosch, N., and Hennig, P. (2022). Fenrir: Physics-enhanced enhanced regression for initial value problems. In Chaudhuri, K., Jegelka, S., Song, L., Szepesvari, C., Niu, G., and Sabato, S., editors, Proceedings of the 39th International Conference Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 21776-21794. PMLR.
[47] Vijayshankar, S., Nabi, S., Chakrabarty, A., Grover, P., and Benosman, M. (2020). Dynamic mode decomposition and robust estimation: Case study of a 2d turbulent boussinesq flow. In 2020 American Control Conference (ACC), pages 2351-2356. IEEE.
[48] Williams, C. and Seeger, M. (2001). Using the Nyström method to speed up kernel machines.
In Advances in Neural Information Processing Systems (NeurIPS).
13


Reviews:
[{"review_id": "AzvLsq2EaJ", "reviewer": "Reviewer_dWp4", "summary": "Filtering (estimation) and smoothing of large-dimensional state-space models are computationally challenging. Exact Kalman filtering, at least insofar as LTI SDEs are considered, is a mature solution for such tasks but the cubic complexity of such methods makes it computationally intractable. This paper proposes an approximate Gaussian filtering and smoothing method that propagates low-rank approximations of the covariance matrices. To enable this, i) Lyapunov equations are projected into a manifold of lower rank in the prediction step combined with ii) square root filtering, which together offer numerically stable and tractable solution. Four examples are provided, with increasing complexity, and in one case runtime analysis is also carried out to support the claims. Two measures are used RMSE and covariance deviation to assess the method’s performance. Corollary 1 provides a “Kalman-like” result (an important contribution of the paper) while Section 3.2 discusses the algorithm of the method for filtering (3.4 for smoothing) and section 3.3 provides the time complexity of the method.", "questions": "In section 1, Please be more specific about such unfavorable properties (line 45). Overall, the issue of GP is not discussed clearly in this section and elsewhere in the manuscript.\n\nAny insights if such approach can be applied to nonlinear systems? What if the dynamics is unknown (A and B)?\n\nI'd suggest moving section 4.4 as either last section or right after 4.1. Currently the flow of the results section is hard to follow: first three examples are given to demonstrate the accuracy of the model with the given two measures, then a discussion of runtime and finally again another result for rainfall data.\n\nIt would be helpful to include a discussion on the impact of snapshots noise (measurement noise) on the results.", "limitations": "N/A", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The proposed method differs from existing ensemble-based approaches in that the low-rank approximations of the covariance matrices are deterministic, rather than stochastic. This is important because it allows the method to reproduce the exact Kalman filter as the low-rank dimension approaches the true dimensionality of the problem.\n\nThe time complexity analysis is well-written, and authors clearly demonstrate their method reduces computational complexity from cubic (for the Kalman filter) to quadratic in the state-space size at most (and linear in some cases—refer to proposition 3 and assumptions therein).", "weaknesses": "My major issue with this paper is there is a large body of work on estimation of high- and infinite-dimensional dynamical systems and this paper does not position itself in such context. There are numerous works on using sparse identification (Sparse reduced-order modeling: Sensor-based dynamics to full-state estimation), Koopman-based (A Robust Data-Driven Koopman Kalman Filter for Power Systems Dynamic State Estimation), balanced-truncation ROMs and Kalman (State Estimation Using a Reduced-Order Kalman Filter), DMD-based ROM and Kalman (Dynamic mode decomposition and robust estimation: Case study of a 2d turbulent boussinesq flow), etc. In above-mentioned works, usually the high-dimensional governing equations is first projected in some appropriate basis (low-rank), and then an appropriate filter is deigned. Once estimation is done, the reconstruction is projected back into the physical domain. This also enables the application of the method to i) nonlinear dynamical system, ii) unknown systems. Such approaches benefit from lower computational cost, similar to the proposed method in the paper. In a sense authors acknowledge, more or less, such line of work in the last paragraph of Section 1. However, when they write “In contrast to our work, their method assumes that the dynamics unfold entirely in a lower-dimensional space, and conditioning on measurements happens by transporting the low-dimensional diffusive process to the full space by a projection that is also assumed to be known” there is no quantitative or qualitative comparison to support their claim. Let’s say for example 1 considered in this paper, authors first project the PDE into low rank ODE, and then design the standard Kalman filter in such a manifold, and once estimation is propagated in time, the low dimension is converted back into actual space. What are the advantages and disadvantages? Please note I am not simply asking for more comparison (in addition to, say, EnKF) but more a conceptual question of what the real advantage of the proposed method is to state of the art for ROM-based estimation. I’m sure authors are aware of works like “Krylov subspace methods for solving large Lyapunov equations” where large scale Lyapunov equations are solved with efficient algorithms. I hope authors clarify their novelties in the rebuttal.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your careful reading of our paper as well as your thoughtful critique. \nIn our response, we would like to divide the major points of contention into two categories: \"Address literature on reduced order methods\" and \"model order reduction versus our approach\". \nThese will be responded to in order, followed by answers to more specific questions raised.\nIt is our assessment that an elaboration on related work in the camera-ready would be sufficient to address your concerns, \nand it is our hope you come to the same conclusion upon reading our response.\n\nBefore proceeding we shall refer to the papers cited in your review as follows \n\n- [1] Sparse reduced-order modeling: Sensor-based dynamics to full-state estimation\n- [2] A Robust Data-Driven Koopman Kalman Filter for Power Systems Dynamic State Estimation \n- [3] Dynamic mode decomposition and robust estimation: Case study of a 2d turbulent boussinesq flow\n- [4] State Estimation Using a Reduced-Order Kalman Filter \n- [5] Krylov subspace methods for solving large Lyapunov equations\n\nand we will also refer to some of the references of our paper by \n\n- [6] P. Sakov and P. R. Oke, “A deterministic formulation of the ensemble Kalman filter: an alternative to ensemble square root filters,” 2008.\n- [7] J. S. Whitaker and T. M. Hamill, “Ensemble Data Assimilation without Perturbed Observations,” 2002.\n- [8] A. Carrassi, M. Bocquet, L. Bertino, and G. Evensen, “Data assimilation in the geosciences: An overview of methods, issues, and perspectives,” 2018.\n- [9] W. Sacher and P. Bartello, “Sampling Errors in Ensemble Kalman Filtering. Part I: Theory,” 2008.\n\n### Address literature on reduced order methods\n\nThe cited papers [1,2,3] are concerned with low-rank modelling of dynamical systems from measured data.\nThe paper [4] is concerned with low-rank approximations of a given model based on classical theory on linear time invariant systems. \nOnce a low rank model is obtained for the phenomena under study then state estimation becomes computationally trivial, in the sense that the model is not very large. \nThese contributions are of course related to the present paper under the broader context of low-rank/ dimensionality reduction methods. \n\nHowever, they do differ on some key points:  \n\n- [1,2,3] allow themselves to construct a (low-rank) model from measured data. Whereas we assume the model class is fixed and the objective is make inference tractable with minimal violence done to the model. \n- [4], like us, takes the model as given but heavily relies on realization theory of stable linear time-invariant systems to obtain a reduced-order model. \nIn contrast, the dynamic low rank method that we employ can be applied to any matrix _differential_ equation. \nConsequently, we expect our present method to be easier to extend to other model classes (i.e. nonlinear).\n\nFinally, we will touch on the mention of paper [5], which develops methods for obtaining low-rank approximations to the solutions of _algebraic_ Lyapunov equations.\nWhereas the problem we tackle with the dynamic low rank method is to obtain low-rank approximations to the solution of Lyapunov _differential_ equations. \nConsequently, the methods of [5] can not directly be brought to bear on our problem.\n\n### Model order reduction versus our approach \n\nThe review brings up example 1 in the paper to support a discussion of model order reduction (i.e. [1,2,3,4]) versus to the present approach \nof simply tracking a low rank approximation to the covariance matrix while keeping the model intact.\n\nBroadly speaking, the approaches of [1,2,3,4] are expected to work well under the following assumption: \n\n- [A1] The state evolves in a low-dimensional manifold (approximately).\n\nWhereas our approach is expected to work well under the following assumption: \n\n- [A2] the error in the state estimate evolves in a low-dimensional manifold (approximately). \n\nThese two assumptions are quite different indeed and of course come with certain advantages and disadvantages. \nWe can not say how the quality of state estimation would differ between the two approaches under [A1], \nwe do however note that the model reduction approach would be more economical, certainly in storage, but perhaps also in computation. \nOn the other hand, under [A2] it is not clear that the model reduction approach would be appropriate in all cases. \nIn particular it could be the case that the state of the system really does develop in a high-dimensional space whereas the uncertainties are concentrated in a lower dimensional subspace \n-- in this case we expect our proposed approach to be more succesful. \n\n**Q: What is meant by \"unfavorable properties\" of the stochastic nature of the EnKF?**\n\nA: Adding sampled noise to the ensemble introduces sampling error [7,9], which leads to erroneous estimates of the approximate moments and motivates works like [7,6] and more work on deterministic ensemble Kalman filters. It further causes spurious correlations in the covariance estimator, which has to be mitigated, e.g., using localization techniques (e.g. [8]). The main point is that it is preferrable to approximately target the Karhunen-Loève truncation rather than sampling a random subspace.\n\n**Q: Your work only regards linear dynamics systems, whereas existing methods are usually applied to nonlinear dynamics**\n\nA: Indeed our algorithm, as presented in the paper, only applies to linear systems. As you are aware, we explicitly discuss this in the limitations section of our work. Indeed, many state-estimation problems are nonlinear, and this setting will certainly be followed-up upon in future work.\n\n\n**Q: Writing suggestions**\n\nA: Many thanks for taking the time to formulate concrete suggestions in order to help improve our work! We will consider those for the camera-ready version.\n"}, {"review_id": "TVOfSFC1kl", "reviewer": "Reviewer_HH4H", "summary": "The paper proposes a deterministic low rank approximation to the Kalman filter in high dimensional settings where the computational complexity of the full Kalman filter is very expensive. They compare to the ensemble Kalman filter and ensemble transform Kalman filter. ", "questions": "The scenario for this problem is linear dynamic systems? Ensemble style Kalman filters are often used in non-linear dynamic systems and have the greatest advantage there; including computational ones, but that is compared to particle filters etc. . The authors cite using ensemble Kalman filters for computational cheapness, but the vast majority of those papers are from not very good journals and their accuracy is suspect...\n\nIn contrast the real data experiments are for atmospheric problems that generally have non-linear dynamics? \n\nHow are you computing the low rank dimension of the ensemble Kalman filter, which is not a low rank method (figure 2)? Neither the ensemble Kalman filter nor the ensemble transform Kalman filter are explicitly attempting to estimate a low rank approximation. It would be more appropriate to compare against a variant of the Kalman filter that is a low rank approximation (or a sparse one at least). There are tons of variants from just a quick google search. \n", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The paper is well written.", "weaknesses": "See questions. The biggest issue is that the comparisons are not against “fair” competing methods (any low rank Kalman filter). They also do not mention anything about the large swaths of low rank variants of Kalman filters.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your critical assessment of our work!\nWe identified your main points of contention to be (a) that the ensemble Kalman filter is not a low-rank method and (b) we are not comparing our method to a low-rank filter. Both of these points are based on a misconception, which we will explain below, followed by addressing your remaining concerns. We believe that our response demonstrates that the points raised do not warrant a rejection of our work.\n\n**Q: The ensemble Kalman filter is \"not a low-rank method\" and thus not a fair comparison.**\n\nA: This statement is incorrect. The EnKF *is* indeed a low-rank method, where the rank of the covariance approximation is determined by the number of ensemble members. Concretely, $\\Sigma \\approx (\\frac{1}{\\sqrt{r-1}}X)(\\frac{1}{\\sqrt{r-1}}X)^\\top$, where the ensemble $X \\in \\mathbb{R}^{n \\times r}$ is of rank $r$ and $r \\ll n$. When comparing our method to the EnKF in our experiments, the \"low-rank dimension\" of the EnKF thus corresponds to the ensemble size.\nIn view of this, we would agree that our experimental comparisons are not exhaustive. \nHowever, we do maintain that our experimental results are certainly sufficient to demonstrate the merits of our approach. \nFurthermore, on the basis of this review and that of dWp4 we are confident that we can more clearly situate the present contribution within the larger body of work in low-rank methods for estimation and model identification.\n\n**Q: The problem statement contains only linear systems, whereas the EnKF is often used for nonlinear systems.**\n\nA: This paper is the first in what we hope to be a line of work. The method is motivated by being entirely deterministic -- which clearly has benefits over Monte-Carlo-based approaches in some relevant scenarios -- and it is established here on linear systems. Indeed our algorithm, as presented in the paper, only applies to linear systems. As you are aware, we explicitly discuss this in the limitations section of our work. Indeed, many state-estimation problems are nonlinear, and this setting will certainly be followed-up upon in future work. The experiment on rainfall estimation succeeds in showcasing the computational feasibility in a high-dimensional state-estimation problem; the assumption of linear dynamics does not impede this argument.\n\n\n**Q: Cited literature that uses ensemble Kalman filters for computational cheapness is questionable**\n\nA: We fundamentally reject this criticism. The foundational paper by Evensen [1], establishing the first version of the EnKF, presents the method as \"a better alternative than solving the traditional and computationally extremely demanding approximate error covariance equation used in the extended Kalman filter\". Our method solves precisely this problem.\nFurther, in an important foundational piece on the EnKF [2], Houtekamer and Mitchell write:\n\n> \"In particular, to make the ensemble approximation feasible, we have to use a fairly small ensemble with many less members than either the number of model coordinates, or the number of independent observations, or the (unknown) dimension of the dynamical system. To nevertheless obtain good results, we must (i) counter the tendency of the ensemble spread to underestimate the true error, and (ii) localize the ensemble covariances. The localization is severe and leads to imbalance in the initial conditions.\"\n\nThis clearly states the agenda of developing a computationally tractable approximate inference scheme and the need to counteract degeneracies that come with the naïve ensemble-based approximation.\n\n[1] G. Evensen, (1994), Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics, _J. Geophys. Res._,  99(C5),  10143–10162.\n\n[2] P. L. Houtekamer and H. L. Mitchell, “Ensemble Kalman filtering,” Q. J. R. Meteorol. Soc., vol. 131, no. 613, pp. 3269–3289, Oct. 2005.\n\n**Q: Much of the existing low-rank Kalman filters are not mentioned.**\n\nA: We will extend the related work section to position the present contribution in the broader context of dimensionality reduction in high-dimensional state-estimation problems."}, {"review_id": "jIRz2msVVH", "reviewer": "Reviewer_U69V", "summary": "The paper presents a deterministic low-rank approximation of the Kalman filter for estimation in high-dimensional linear Gaussian setting. This setting is important for many problems that arise in meteorology and solved by Ensemble Kalman Filter (EnKF). The paper presents formula for low rank approximation of the covariance for both prediction and correction steps and demonstrates the proposed approach in different scenarios in comparison with EnKF. ", "questions": "1- What is a reference for stochasticity in EnKF introduces unfavorable properties. It will be good to be more concrete here, and explain what are exactly those unfavorable properties, since it is a central motivation for your approach. \n\n2- The computational time plot in figure 4 is a bit surprising. How is that low rank approximation has smaller computational time compared to EnKF?  \n\n3- Can you produce a figure that demonstrates error as a function of raw computational time? \n\n4- What is the reference for O(nr^2) complexity of truncated SVD? ", "limitations": "yes", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- paper is well-written, has clear message\n- it considers an important problem in meteorology\n- the proposed approach is motivated well \n- comprehensive numerical experiments are given", "weaknesses": "- the main weakness is the fact that in most high-dimensional data assimilation problems, one has has only access to a simulator for the forward dynamic update. It is not clear how the proposed approach perform low rank approximation in that case. In EnKF, one simulates each member of the ensemble according to the dynamic model\n- the numerical code for reproducing the results are not given\n- the relevance to NeurIPS audience is weak. \n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your detailed and positive assessment! We hope that we can additionally resolve your main point of criticism in the following and thereby convince you to raise your score.\n\n**Q: It is not clear how to apply this method in case one has to forward the dynamics with a numerical simulator.**\n\nA: This is a very important point, the relevance of which we are very aware of. Since dynamics, which require a numerical simulator for forward propagation, are usually nonlinear, this question will receive systematic treatment in follow-up work, while this paper establishes the foundations of the method on a purely linear setup.\nNote that a simple, linear example of this case is showcased by the linear-advection example (Section 4.1; Fig. 1), which is a common baseline in papers that propose ensemble Kalman filter methods (e.g. [1], [2]). There, the simulation step is a simple matrix-vector product with a circulant (shift) matrix. Instead of individual ensemble members, the mean is simulated forwards and the linear operator is used to forecast the low-rank covariance matrix (as described).\n\n\n**Q: The code for reproducing the results is not given**\n\nA: The code is ready to be made public and will be released upon acceptance of the paper.\n\n\n**Q: The relevance to NeurIPS audience is weak.**\n\nA: We strongly disagree with this statement and hope to convince you of the opposite by addressing it here. The recent momentum in the domain of ML in the sciences, like, for instance, meteorology, geophysics, oceanography, etc., calls for efficient, yet accurate tools for the analysis of high-dimensional time series. Especially in sparse-data regimes, Gaussian-process models and hybrid physics-informed/data-driven ideas from the field of \"data assimilation\" become more and more attractive and often are preferred over purely data-driven approaches. Our method (and planned follow-up work) aims to serve as a well-founded, simple-to-use drop-in replacement for approximate Bayesian filtering/smoothing for exactly such applications and we hope to gain momentum in both the ML community and the sciences.\n\n\n\n**Q: What do you mean by \"unfavorable properties\" that are introduced by the stochasticity of the EnKF?**\n\nA: Adding sampled noise to the ensemble introduces sampling error [3,5], which leads to erroneous estimates of the approximate moments and motivates works like [3,1] and generally all work on deterministic ensemble Kalman filters; and it causes spurious correlations in the covariance estimator, which has then to be mitigated by localization techniques (e.g. [4]). The main point is that it is preferrable to approximately target the Karhunen-Loève truncation rather than sampling a random subspace.\n\n\n**Q: How is the RRKF faster than the EnKF in Fig. 4?**\n\nA: It is always difficult to compare absolute runtimes, due to nuances in the respective implementation. Even though we implemented all the methods ourselves and took great care in optimizing each of them to the same degree, there are likely variations. The point of our runtime analysis lies in the **asymptotic behavior**. This plot is not to show that our algorithm is faster than the EnKF. It is to show that the asymptotic runtimes of the ensemble-based algorithms is preserved by our method, while it has some major advantages over the stochastic counterparts, as described in the paper.\n\n\n**Q: Can you produce a figure that demonstrates error as a function of raw computational time?**\n\nA: Thanks for this suggestion. We included a draft of such figures in the PDF file attached to the general author rebuttal above. We would like to add those to the supplementary material of the camera-ready version.\n\n\n**Q: Why is the truncated SVD complexity $O(nr^2)$?**\n\nA: On tall and wide matrices the rank is upper-bounded by the smaller of both dimensions. It is therefore sufficient to compute only a subset ($r$) of the singular vectors of such a matrix, since the remaining singular values are zero. This is sometimes referred to as \"compact SVD\" (or \"reduced SVD\" (RSVD), or \"economic SVD\").\nThis is detailed, e.g. in Chapter 5 in [6]. We will add this reference to the camera-ready version.\n\n\n[1] P. Sakov and P. R. Oke, “A deterministic formulation of the ensemble Kalman filter: an alternative to ensemble square root filters,” _Tellus A: Dynamic Meteorology and Oceanography_, vol. 60, no. 2, p. 361, Jan. 2008.\n\n[2] G. Evensen, F. C. Vossepoel, and P. J. van Leeuwen, _Data Assimilation Fundamentals: A Unified Formulation of the State and Parameter Estimation Problem_. in Springer Textbooks in Earth Sciences, Geography and Environment, 2022.\n\n[3] J. S. Whitaker and T. M. Hamill, “Ensemble Data Assimilation without Perturbed Observations,” _Mon. Wea. Rev._, vol. 130, no. 7, pp. 1913–1924, Jul. 2002.\n\n[4] A. Carrassi, M. Bocquet, L. Bertino, and G. Evensen, “Data assimilation in the geosciences: An overview of methods, issues, and perspectives,” _WIREs Clim Change_, vol. 9, no. 5, Sep. 2018.\n\n[5] W. Sacher and P. Bartello, “Sampling Errors in Ensemble Kalman Filtering. Part I: Theory,” _Monthly Weather Review_, vol. 136, no. 8, pp. 3035–3049, Aug. 2008.\n\n[6] Golub, Gene H., and Charles F. Van Loan. Matrix computations. JHU press, 2013."}, {"review_id": "ZthYv9FqaI", "reviewer": "Reviewer_2SHM", "summary": "The paper proposes an approximation method for Kalman filtering of intractable, high-dimensional Gauss-Markov settings with linear observations. It adopts existing DLRA methods to compute a discrete-time representation of the sampled process and suggests a recursive filtering and smoothing methods in a low-rank formulation. \n\n", "questions": "1. The recursive methods assumes knowing $\\Sigma_{l-1}$ at each step, but don't mention how to initialize $\\Sigma_{1}$ at $t_1$. Is it just $Q_{t_1}$?\n\n2. The estimations on Rainfall in Australia regression (Fig.5) don't seem very accurate, especially at spatial points where data deviates from the mean value. Authors claim that \"the RRKF achieves high-quality approximate estimates\", but I'm not convinced. Can they support this claim quantitively?\n\nMinor comments:\n\nEqn. (4) and (16) are somewhat unclear. It would be helpful to clearly state matrices dimensions.\n\nIn l.146 it is stated that truncated SVD complexity is $O(nr^2)$, which seems to be right. But can you explain/reference this claim?\n", "limitations": "The method assumes the low rank $r$. This limitation is not discussed.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "To the best of my knowledge, the proposed filtering scheme is novel and is shown to outperform existing EnKFs under the same configurations. As opposed to EnKF approximation methods, this method is deterministic and capable of accurately restoring the optimal filter given the right parameterization. \n\nDerivation of the recursive filter seems correct, and filtering is done in a square-root fashion, which is stable and efficient. The recursion step is shown to be efficient (although under very strong assumptions). The method additionally suggests an efficient calculation for the observation likelihood.", "weaknesses": "The method assumes the low rank $r$, which is a drawback of many low-rank approximation methods. Even if it is beyond the scope of the paper, there should be a short discussion about choosing this parameter or mention this as a limitation.\n\nIn Figs. 2, 3 the RRKF is compared with two other methods, both are Monte-Carlo, and are not expected to restore the exact filter. A righteous comparison should include some deterministic algorithms e.g. [10], [r1] below. \n\n[r1] MLA\tFarrell, Brian F., and Petros J. Ioannou. \"State estimation using a reduced-order Kalman filter.\" Journal of the Atmospheric Sciences 58.23 (2001): 3666-3680\n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your detailed and positive assessment of our work. In the following we will address your remaining concerns, which we hope further improves and solidifies your valuation of the quality of our work.\n\n**Q: The method introduces the low-rank dimension $r$ as a tunable parameter of the method. How would one choose $r$?**\n\nA: Thank you for raising this important point. It is true that the method introduces the low-rank dimension $r$ as a tunable parameter. We will emphasize and discuss this in our revision. In fact, the dynamic low-rank approximation (DLRA) literature offers a potentially quite interesting direction for follow-up work, which might automate the choice of $r$ [2].\n\n\n**Q: The method should be compared to reduced-order Kalman filters ([1]).**\n\nA: [1], like us, takes the model as given but heavily relies on realization theory of stable linear time-invariant systems to obtain a reduced-order model. \nIn contrast, the dynamic low rank method that we employ can be applied to any matrix differential equation. \nConsequently, we expect our present method to be easier to extend to other model classes (i.e. nonlinear). Therefore, we view a comparison to the method presented in [1] as not critical to demonstrate the value of our method. Nevertheless, we will take the opportunity to more clearly situate the present contribution within the larger body of work in low-rank methods for estimation and model identification.\n\n\n**Q: How is the state covariance matrix initialized?**\n\nA: In our experiments, we compute the stationary mean and a low-rank factorization of the stationary covariance matrix of the prior and condition the stationary moments on the first measurement in the respective time-series dataset. We will add this detail to the revised version.\n\n\n**Q: The rainfall predictions in the final experiment could be improved**\n\nA: You are right in noticing that the rainfall model yields less accurate predictions, as opposed to solving a PDE. Essentially, the rainfall data is interpolated approximately with a simple Matérn Gaussian-process model. This, of course, makes the prediction significantly cheaper, as well.\nMost importantly, the focus of this experiment is the high dimensionality of the state-estimation problem, which is solved efficiently and -- given the simple model and the significant rank compression -- accurately. We will clarify the focus of the experiment accordingly in the paper.\n\n\n**Q: The notation for square-root factors of covariances (e.g. Eq. (4) and (16)) is somewhat unclear.**\n\nA: We will clarify the notation in a revised version of the paper. For example, $(\\Phi \\Sigma^{1/2} \\quad Q^{1/2})$ denotes a block matrix that is a square-root factor of the predicted covariance matrix, in that the outer product $(\\Phi \\Sigma^{1/2} \\quad Q^{1/2})(\\Phi \\Sigma^{1/2} \\quad Q^{1/2})^\\top = \\Phi\\Sigma\\Phi^\\top + Q$. If $\\Sigma^{1/2}$ and $Q^{1/2}$ are $n\\times r$ - matrices -- i.e. low-rank square-root factors of $\\Sigma$ and $Q$, respectively -- then the block matrix is in turn a low-rank square-root factor of the predicted covariance matrix.\n\n\n**Q: Why is the truncated SVD complexity $O(nr^2)$?**\n\nA: On tall and wide matrices the rank is upper-bounded by the smaller of both dimensions. It is therefore sufficient to compute only a subset ($r$) of the singular vectors of such a matrix, since the remaining singular values are zero. This is sometimes referred to as \"compact SVD\" (or \"reduced SVD\" (RSVD), or \"economic SVD\").\nThis is detailed, e.g. in Chapter 5 in [3]. We will add this reference to the camera-ready version.\n\n\n[1] B. F. Farrell and P. J. Ioannou, \"State estimation using a reduced order Kalman filter,\" 2001.\n\n[2] G. Ceruti, J. Kusch, and C. Lubich, “A rank-adaptive robust integrator for dynamical low-rank approximation,” 2022.\n\n[3] Golub, Gene H., and Charles F. Van Loan. Matrix computations. JHU press, 2013."}, {"review_id": "CBN9F9XT3T", "reviewer": "Reviewer_e5WD", "summary": "This paper develops a method for performing approximate Kalman filtering and smoothing using a low-rank approximation.  Fully low-rank predict, update and smoothing equations are derived.  Under certain reasonable assumptions, these equations reduce the cubic complexity to linear complexity (and quadratic if the assumptions are violated), while being “approximately exact” compared to Monte Carlo-based ensemble methods.  The method is evaluated on a range of problems and appears to perform favourably.  Sensible limitations and future opportunities are highlighted.\n", "questions": "Also see weaknesses and limitations.\n\nIs operating in the continuous time domain strictly necessary?  My understanding is that a lot of complexity (at least notationally) goes away when you consider discrete time/equally spaced observations.  I think all of your experiments used equally spaced observations as well.  It could be simpler if you present the material in the discrete time domain, and then have a paragraph outlining that this also extends to the continuous time/SDE domain.  \n", "limitations": "A limitation is that the authors do not consider bounding the error introduced by the low-rank approximation.  For instance, the LML in (12) is for the approximate model.  There is no comment on the approximation gap between the approximate model and the full-rank model.  I am not fully sure on this, but I think it is possible to construct a bound based on the truncated parts of the spectrum of the matrix being truncated.  If this isn’t possible, some qualitative discussion would help reinforce this point.  If you were able to construct a bound – any bound! – on the error, then that would put a nice bow on the approximation, as opposed to relying on evaluating the error empirically. \n", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "This paper addresses a critical issue – the super-linear (cubic) complexity of exact Kalman filtering/smoothing (KF/S).  This makes stable and closed-form estimation intractable for anything other than small problems.  The proposed solution is appealing: introduce a low-rank factorization that, through some linear algebra and recent solvers, is linear (or quadratic) in the rank.  While low-rank factorizations are common, I have not seen it applied to KF/S, a domain where it is very naturally appealing.  \n\nThe detail in which the authors walk through the derivation is excellent, and appears to be mathematically sound (although I am not super familiar with this application domain).  The purely theoretical contributions of the paper are not huge; but are well-judged, interesting, apt for the problem at hand (modulo Limitations), and dovetail nicely with recent literature.  \n\nThe experimental evaluation touches on a number of interesting examples that are high dimensional, and for which the proposed method seems to excel.  The final application, although a little off-the-wall, really does a good job in highlighting how high-dimensional this method can go.  Other experiments do a good job of comparing asymptotic to non-asymptotics in a range of domains.  \n\nI think there is definitely an audience at NeurIPS that will enjoy and really benefit from this work.  I also believe this work will engender follow-up work which will further build this direction.  The paper is on the whole well-written (modulo Weakness B) and is nicely visually presented.  Figures are very well prepared, with clearly legible fonts and clear colours.  ", "weaknesses": "The paper, on the whole, is very good, but I do think there are some areas for improvement.\n- A:  The experimental evaluation is, on the whole, good.  Figure 1 is excellent in visually depicting how covariance matrices can be low rank.  WIth that in mind, I wonder if there is a simple, low-dimensional model that the authors can present that can be easily visualised.  I don’t have a great “feel” for any of the examples presented.  The results for this experiment certainly don’t need to be SoTA (since existing methods will probably also perform well).  Some visualisations (PCA projections of state, low-rank factors etc) can then be visualised to highlight the differences between methods and build reader intuition.  This example could be as simple as a synthetic LDS with dimension ~10.  Note:  I wouldn’t expect to see this during the rebuttal period, but I would _strongly_ encourage the authors to consider including it in a camera ready.\n- B:  [I am aware this is a very general and somewhat hard-to-address comment, so please take these more as suggestions]  The paper itself is very notationally dense, to the point of obfuscating the simplicity and applicability of the method.  While the authors are to be commended for their attention to detail, I think parts of the paper can be thinned out, with more technical content being dropped to the supplement.  This will result in a paper that is more immediately digestible, and will have a larger reach and impact.  Some examples:\n  -  I think Sec 3.1 is great for building reader intuition, but this intuition is then lost a little bit by having to wade through long math equations, where the details aren’t actually super important.  \n  - I  wonder if some of the smoothing equations can be cut to the supplement, with just the final forms retained and note directing the interested reader to the supplement.  This newly cleared space can then be used for including higher-level sketch proofs, intuition-building worked examples or text, algorithm blocks, diagrams, experiment visualisations etc. \n  - Including a banner figure (and maybe even an algorithm block) explaining the whole method would be nice, just to visually complement the math.  \n  - Maybe consider including a table outlining the positives and negatives of each approach under particular operating conditions?  Then the benefits are stated clearly and concisely outside of the body text.  \n  - **The key takeaway** is that you have already solved the hard equations and derived an elegant solution!  Therefore, I would encourage you to strip out as much complexity from the main text as possible, and add extra emphasis to the end result – even something as cartoonish as a callout box with the key equations (akin to how Sarkka has a full equation block per model/inference methodology). \n\n## Minor weaknesses:\n- It would be nice if the DLRA/BUG integrator/solver were introduced in the main text.  You highlight it as a key component, but there is not really much discussion of its implementation, assumptions, limitations, complexity etc.  \n- I don't quite understand the $(\\Phi \\Sigma \\quad Q)$ notation in (16) (specifically the white space between the matrices).  Please explicitly define it somewhere. \n\nIf these changes (and the limitations)  were made/commented on, then I would consider upgrading my review. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your strong review and your detailed, positive assessment of our work! We greatly appreciate your appraisal that \"there is definitely an audience at NeurIPS that will enjoy and really benefit from this work\".\nIn the following we will address your remaining concerns, which we hope further improves and solidifies your valuation of the quality of our work.\n\n**Q: Writing suggestions**\n\nA: Many thanks for taking the time to formulate concrete suggestions in order to help improve our work! We will consider them for the camera-ready version.\n\n\n**Q: Is there a simple, low-dimensional model that can be easily visualized, in order to give more intuition into the workings of the method?**\n\nA: We have given this idea some thought. It is not easy to construct an instructive example that is low-dimensional (and thus easily visualizable) and can be further compressed with a low-rank method. The full potential of our contribution unfolds exactly in very large systems in which the state-covariance can be tracked approximately on a low-dimensional manifold. Thus, we went in the opposite direction to focus on this by showcasing the method on a very large dynamical system with the example on rainfall prediction. There, we wanted to convey the intuition that a lot of locations have to be tracked but neighboring points will exhibit similar error estimates, which justifies the low-rank assumption.\n\n\n**Q: Is operating in the continuous time domain strictly necessary?**\n\nA: In many cases, continuous formulations of dynamics are sparse, which does not in general transfer to the discretized model. This would thus not allow for exploiting the sparse properties of the continuous system for computational tractability.\n\n\n**Q: Could a bound to the error of the method be provided?**\n\nA: This brings up an interesting point. However, it is not at present clear how to approach an error analysis in a tractable way. While there are initial ideas for _local_ error bounds, we saw value in postponing their systematic investigation to future work in order to place more focus on reader intuition behind the approach.\n\n\n**Q: It would be nice if the DLRA/BUG integrator/solver were introduced in the main text.**\n\nA: We give an intuition behind the method in Section 2.2 and the concrete algorithm as we use it in Section C of the supplementary material. We further provide all necessary references to the DLRA literature. In our opinion this is a good usage of the limited space for an existing method with a large body of existing literature.\n\n\n**Q: The notation for square-root factors of covariances (e.g. Eq. (16)) is somewhat unclear.**\n\nA: We will clarify the notation in a revised version of the paper. For example, $(\\Phi \\Sigma^{1/2} \\quad Q^{1/2})$ denotes a block matrix that is a square-root factor of the predicted covariance matrix, in that the outer product $(\\Phi \\Sigma^{1/2} \\quad Q^{1/2})(\\Phi \\Sigma^{1/2} \\quad Q^{1/2})^\\top = \\Phi\\Sigma\\Phi^\\top + Q$. If $\\Sigma^{1/2}$ and $Q^{1/2}$ are $n\\times r$ - matrices -- i.e. low-rank square-root factors of $\\Sigma$ and $Q$, respectively -- then the block matrix is in turn a low-rank square-root factor of the predicted covariance matrix."}]

Author rebuttal:
We are very grateful to all reviewers for their high-quality reviews and for providing constructive feedback. We are inspired by the reviewers' confirmation that our work addresses a "critical issue" (e5WD) and an "important problem" (U69V) to which our "proposed solution is appealing" (e5WD).
We would like to thank the reviewers for judging our manuscript as "nicely visually presented" (e5WD), "well-written" (e5WD, U69V, HH4H, dWp4), conveying a "clear message" (U69V); and for highlighting the novelty of the method (e5WD,2SHM) and the "excellent" (e5WD) and "mathematically sound" (e5WD) derivation. We are particularly encouraged by the assessment that "there is definitely an audience at NeurIPS that will enjoy and really benefit from this work" (e5WD).

Further, we are grateful to the reviewers for suggesting improvements, which will greatly benefit a revised version of the manuscript.
The most significant objection is that the relation of the proposed method to other existing approaches to low-rank filtering should be made more explicit. We will take the opportunity to situate the present contribution within the larger body of work in low-rank methods for estimation and model identification in the camera-ready version.

Meta-review: The paper presents a new low-rank approximation method for Gaussian filtering and smoothing in high-dimensional dynamical systems. It significantly departs from previous ensemble-based methods. 

The authors develop a time complexity analysis that is well-detailed and compelling.Addressing the issue of super-linear complexity in exact Kalman filtering and smoothing, the paper proposes a solution using low-rank factorization. This solution reduces complexity to linear or quadratic levels, making it effective for larger problems.

Although the paper doesn't provide major theoretical contributions, it addresses a crucial challenge and aligns with the most recent literature. The paper's strength lies in its practical application demonstrated through a variety of experiments. This validates the proposed method effectively. Feedback suggests further improvements, such as the inclusion of specific baselines, and the authors should do so in the camera-ready version. Based on these observations, the paper is recommended for acceptance.
Decision: Accept (poster)