Year: 2021

Abstract:
We perform approximate inference in state-space models with nonlinear state transitions. Without parameterizing a generative model, we apply Bayesian update formulas using a local linearity approximation parameterized by neural networks. It comes accompanied by a maximum likelihood objective that requires no supervision via uncorrupt observations or ground truth latent states. The optimization backpropagates through a recursion similar to the classical Kalman filter and smoother. Additionally, using an approximate conditional independence, we can perform smoothing without having to parameterize a separate model. In scientific applications, domain knowledge can give a linear approximation of the latent transition maps, which we can easily incorporate into our model. Usage of such domain knowledge is reflected in excellent results (despite our model’s simplicity) on the chaotic Lorenz system compared to fully supervised and variational inference methods. Finally, we show competitive results on an audio denoising experiment.

text: We perform approximate inference in state-space models with nonlinear state transitions. Without parameterizing a generative model, we apply Bayesian update formulas using a local linearity approximation parameterized by neural networks. It comes accompanied by a maximum likelihood objective that requires no supervision via uncorrupt observations or ground truth latent states. The optimization backpropagates through a recursion similar to the classical Kalman filter and smoother. Additionally, using an approximate conditional independence, we can perform smoothing without having to parameterize a separate model. In scientific applications, domain knowledge can give a linear approximation of the latent transition maps, which we can easily incorporate into our model. Usage of such domain knowledge is reflected in excellent results (despite our model’s simplicity) on the chaotic Lorenz system compared to fully supervised and variational inference methods. Finally, we show competitive results on an audio denoising experiment.

heading: 1 INTRODUCTION
text: Many sequential processes in industry and research involve noisy measurements that describe latent dynamics. A state-space model is a type of graphical model that effectively represents such noise-afflicted data (Bishop, 2006). The joint distribution is assumed to factorize according to a directed graph that encodes the dependency between variables using conditional probabilities. One is usually interested in performing inference, meaning to obtain reasonable estimates of the posterior distribution of the latent states or uncorrupt measurements. Approaches involving sampling (Neal et al., 2011), variational inference (Kingma & Welling, 2013), or belief propagation (Koller & Friedman, 2009) have been proposed before. Assuming a hidden Markov process (Koller & Friedman, 2009), the celebrated Kalman filter and smoother (Kalman, 1960; Rauch et al., 1965) are classical approaches to solving the posterior inference problem. However, the Markov assumption, together with linear Gaussian transition and emission probabilities, limit their flexibility. We present filtering and smoothing methods that are related to the classical Kalman filter updates but are augmented with flexible function estimators without using a constrained graphical model. By noting that the filtering and smoothing recursions can be back-propagated through, these estimators can be trained with a principled maximum-likelihood objective reminiscent of the noise2noise objective (Lehtinen et al., 2018; Laine et al., 2019). By using a locally linear transition distribution, the posterior distribution remains tractable despite the use of non-linear function estimators. Further, we show how a linearized smoothing procedure can be applied directly to the filtering distributions, discarding the need to train a separate model for smoothing.
To verify what is claimed, we perform three experiments. (1) A linear dynamics filtering experiment, where we show how our models approximate the optimal solution with sufficient data. We also report that including expert knowledge can yield better estimates of latent states. (2) A more challenging chaotic Lorenz smoothing experiment that shows how our models perform on par with recently proposed supervised models. (3) An audio denoising experiment that uses real-world noise showing practical applicability of the methods.
Our contributions can be summarized as follows.
1. We show that the posterior inference distribution of a state-space model is tractable while parameter estimation is performed by neural networks. This means that we can apply the classical recursive Bayesian updates, akin to the Kalman filter and smoother, with mild assumptions on the generative process.
2. Our proposed method is optimized using maximum likelihood in a self-supervised manner. That is, ground truth values of states and measurements are not assumed to be available for training. Still, despite our model’s simplicity, our experiments show that it performs better or on par with several baselines.
3. We show that the model can be combined with prior knowledge about the transition and emission probabilities, allowing for better applicability in low data regimes and incentivizing the model to provide more interpretable estimates of the latent states.
4. A linearized smoothing approach is presented that does not require explicit additional parameterization and learning of the smoothing distribution.

heading: 2 RELATED WORK
text: Becker et al. (2019) provide a detailed discussion of recent related work, which we build on here and in table 1. An early method that extends the earlier introduced Kalman filter by allowing nonlinear transitions and emissions is the Extended Kalman filter (Ljung, 1979). It is limited due to the naive approach to locally linearize the transition and emission distributions. Furthermore, the transition and emission mechanisms are usually assumed to be known, or estimated with Expectation Maximization (Moon, 1996). More flexible methods that combine deep learning with variational inference include Black Box Variational Inference (Archer et al., 2015), Structured Inference Networks (Krishnan et al., 2017), Kalman Variational Autoencoder (Fraccaro et al., 2017), Deep Variational Bayes Filters (Karl et al., 2017), Variational Sequential Monte Carlo (Naesseth et al., 2018) and Disentangled Sequential Autoencoder (Yingzhen & Mandt, 2018). However, the lower-bound objective makes the approach less scalable and accurate (see also Becker et al. (2019). Furthermore, all of the above methods explicitly assume a graphical model, imposing a strong but potentially harmful inductive bias. The BackpropKF (Haarnoja et al., 2016) and Recurrent Kalman Network (Becker et al., 2019) move away from variational inference and borrow Bayesian filtering techniques from the Kalman filter. We follow this direction but do not require supervision through ground truth latent states or uncorrupt emissions. Satorras et al. (2019) combine Kalman filters through message passing with graph neural networks to perform hybrid inference. We perform some of their experiments by also incorporating expert knowledge. However, contrary to their approach, we do not need supervision. Finally, concurrently to this work, Revach et al. (2021) develop KalmanNet. It proposes similar techniques but evaluates them in a supervised manner. The authors, however, do suggest that an unsupervised approach can also be feasible. Additionally, we more explicitly state what generative assumptions are required, then target the posterior distribution of interest, and develop the model and objective function from there. Moreover, the current paper includes linearized smoothing (section 6), parameterized smoothing (appendix A), and the recurrent model (appendix C). We also denote theoretical guarantees under the noise2noise objective.
scalable state est. uncertainty noise dir. opt. self-sup.
Ljung (1979) ✓/× ✓ ✓ ✓ × ×
Hochreiter et. al. (1997) ✓ ✓ ✓/× ✓ ✓ × Cho et al. (2014) ✓ ✓ ✓/× ✓ ✓ ×
Wahlström et al. (2015) ✓ ✓ ✓/× × ✓ × Watter et al. (2015) ✓ × ✓ ✓ × ✓
Archer et al. (2015) ✓/× × ✓ ✓ × ✓ Krishnan et al. (2017) ✓ × ✓ ✓ × ✓ Fraccaro et al. (2017) ✓/× × ✓ ✓ × ✓ Karl et al. (2017) ✓ × ✓ ✓ × ✓ Naesseth et al. (2018) ✓ × ✓ ✓ × ✓ Yingzhen et al. (2018) ✓ × ✓ × × ✓
Rangapuram et al. (2018) ✓/× ✓ (1D) ✓ × ✓ ✓ Doerr et al. (2018) × ✓ ✓ ✓ ✓ ✓
Satorras et al. (2019) ✓ ✓ × ✓ ✓ × Haarnoja et al. (2016) ✓ ✓ ✓ ✓ ✓ × Becker et al. (2019) ✓ ✓ ✓ ✓ ✓ ×
Ours ✓/× ✓ ✓ ✓ ✓ ✓
Table 1: We compare whether algorithms are scalable, state estimation can be performed, models provide uncertainty estimates, noisy or missing data can be handled, optimization is performed directly and if supervision is required. “✓/×” means that it depends on the parameterization.

heading: 3 GENERATIVE MODEL ASSUMPTIONS
text: In this section, we explicitly state the model’s generative process assumptions. First, we assume that we can measure (at least) one run of (noise-afflicted) sequential data y0:K := (y0, . . . ,yK), where each yk ∈ RM , k = 0, . . . ,K. We abbreviate: yl:k := (yl, . . . ,yk) and y<k := y0:k−1 and y≤k := y0:k and y−k := (y0:k−1,yk+1:K). We then assume that y0:K is the result of some possibly non-linear probabilistic latent dynamics, i.e., of a distribution p(x0:K), whose variables are given by x0:K := (x0, . . . ,xK) with xk ∈ RN . Each yk is assumed to be drawn from some shared noisy emission probability p(yk | xk). The joint probability is then assumed to factorize as:
p(y0:K ,x0:K) = p(x0:K) K∏ k=0 p(yk | xk). (1)
Further implicit assumptions about the generative model are imposed via inference model choices (see section 7). Note that this factorization encodes several conditional independences like
yk ⊥⊥ (y−k,x−k) | xk. (2) Typical models that follow these assumptions are linear dynamical systems, hidden Markov models, but also nonlinear state-space models with higher-order Markov chains in latent space like presented in fig. 2.
In contrast to other approaches (e.g., Krishnan et al. (2015); Johnson et al. (2016); Krishnan et al. (2017)) where one tries to model the latent dynamics with transition probabilities p(xk | xk−1) and possibly non-linear emission probabilities p(yk | xk), we go the other way around. We assume that all the non-linear dynamics are captured inside the latent distribution p(x0:K), where at this point we make no further assumption about its factorization, and the emission probabilities are (wellapproximated with) a linear Gaussian noise model:
p(yk | xk) = N (yk |Hxk,R) , (3) where the matrix H represents the measurement device and R is the covariance matrix of the independent additive noise. We make a brief argument why this assumption is not too restrictive. First, if one is interested in denoising corrupted measurements, any nonlinear emission can be captured directly inside the latent states xk. To see this, let zk ∈ RN−M denote non-emitted state variables. We then put xk := [yk zk]
⊤, where yk is computed by applying the nonlinear emission to zk. We thus include the measurements in the modeled latent state xk. Then we can put H := [ IM 0M×(N−M) ] . Second, techniques proposed by Laine et al. (2019) allow for nonGaussian noise models, relaxing the need for assumption eq. (3). Third, we can locally linearize the
emission (Ljung, 1979). Finally, industrial or academic applications include cases where emissions are (sparse) Gaussian measurements and the challenging nonlinear dynamics occur in latent space. Examples can be found in MRI imaging (Lustig et al., 2007) and radio astronomy (Thompson et al., 2017).

heading: 4 PARAMETERIZATION
text: In this section, we show how we parameterize the inference model. A lot of the paper’s work relies on established Bayesian filtering machinery. However, for completeness, we like to prove how all the update steps remain valid while using neural networks for function estimation.
Given our noisy measurements y0:K = (y0, . . . ,yK) we want to find good estimates for the latent states x0:K = (x0, . . . ,xK), which generated y0:K . For this, we want to infer the marginal conditional distributions p(xk | y≤k) or p(xk | y<k) (for forecasting), for an online inference approach (filtering); and p(xk | y0:K) or p(xk | y−k), for a full inference approach (smoothing). In the main body of the paper, we only consider filtering. Smoothing can be performed similarly, which is detailed in the supplementary material (appendix A).
We start with the following advantageous parameterization: p(xk | xk−1,y<k) = N ( xk | F̂k|<k xk−1 + êk|<k, Q̂k|<k ) , (4)
where F̂k|<k := F̂k|<k(y<k), êk|<k := êk|<k(y<k) and Q̂k|<k := Q̂k|<k(y<k) are parameterized with neural networks. Next, we have available
p(xk−1 | y≤(k−1)) = N ( xk−1 | x̂(k−1)|≤(k−1), P̂(k−1)|≤(k−1) ) , (5)
i.e., the previous time-step’s conditional marginal distribution of interest. For k = 1, this is some initialization. Otherwise, it is the result of the procedure we are currently describing. We use this distribution to evaluate the marginalization
p(xk | y<k) = ∫ p(xk | xk−1,y<k) p(xk−1 | y<k) dxk−1
= N ( xk | x̂k|<k, P̂k|<k ) ,
(6)
with
x̂k|<k(y<k) = F̂k|<k x̂k−1|≤k−1 + êk|<k, P̂k|<k(y<k) = F̂k|<k P̂k−1|≤k−1 F̂ ⊤ k|<k + Q̂k|<k.
(7)
Note that the distributions under the integral eq. (6) are jointly Gaussian only because of the parameterization eq. (4). Hence, we can evaluate the integral analytically.
Finally, to obtain the conditional p(xk | y≤k) = p(xk | yk,y<k) we use Bayes’ rule:
p(xk | yk,y<k) = p(yk | xk,y<k) · p(xk | y<k)
p(yk | y<k) eq. (2) = p(yk | xk) p(yk | y<k) · p(xk | y<k). (8)
Equation (3) and the result eq. (6) allow us to also get an analytic expression for p(xk | y≤k) = N ( xk | x̂k|≤k, P̂k|≤k ) (9)
with the following abbreviations:
x̂k|≤k(y≤k) := x̂k|<k + K̂k (yk −H x̂k|<k), P̂k|≤k(y≤k) := P̂k|<k − K̂k H P̂k|<k. (10) We introduce the Kalman gain matrix similar to the classical formulas:
K̂k := P̂k|<k H ⊤ ( H P̂k|<k H ⊤ +R )−1 . (11)
Note that taking the matrix inverse at this place in eq. (11) is more efficient than in the standard Gaussian formulas (for reference presented in appendix B) if M ≤ N , which holds for our experiments.
This completes the recursion: we can use eq. (9) for a new time-step k + 1 by plugging it back into eq. (6). We have shown how estimating a local linear transition using neural networks in eq. (4) ensures that all the recursive update steps from the Kalman filter analytically hold without specifying and estimating a generative model.
We note that we could also have parameterized p(xk | y<k) = N ( xk | x̂k|<k(y<k), P̂k|<k(y<k) ) (12)
with x̂k|<k(y<k) and P̂k|<k(y<k) directly estimated by a neural network. This has the advantage that we do not rely on a local linear transition model. However, it also means that we are estimating xk without any form of temporal regularization. Additionally, it is harder to incorporate prior knowledge about the transitions maps into such a model. Nonetheless, we detail this parameterization further in the supplementary material (appendix C) and include its performance in our experiments in section 8.
To conclude the section, we like to discuss some of the limitations of the approach. 1. The Gaussianity assumption of eq. (4) ensures but also restricts eq. (8) and eq. (9) to these forms. That is, we make a direct assumption about the form of the posterior p(xk | y≤k). Defending our case, we like to point out that methods such as variational inference (Krishnan et al., 2017) or posterior regularization (Ganchev et al., 2010) also make assumptions (e.g., mean-field Gaussian) about the posterior. 2. Since we did not explicitly specify a factorization of p(x0:K), we cannot ensure that the distributions we obtain from the above procedure form a posterior to the ground truth generative process. This does not mean, however, that we cannot perform accurate inference. Arguably, not making explicit assumptions about the generative process is preferred to making wrong assumptions and using those for modeling, which can be the case for variational auto-encoders. 3. The local linearity assumption is justifiable if the length between time-steps is sufficiently small. However, note that the model is more flexible than directly parameterizing eq. (12) (see appendix C) since it reduces to that case by putting Fk := 0 for all k.

heading: 5 FITTING
text: We have shown in the previous section how parameterization of a local linear transition model leads to recursive estimation of p(xk | y<k) and p(xk | y≤k) for all k using classical Bayesian filtering formulas. The inference is only effective if the estimates F̂k|<k, êk|<k and Q̂k|<k from eq. (4) are accurate. We can use the parameterization p(xk | y<k) of eq. (6), the emission model p(yk | xk) = N (yk | H xk,R) from eq. (3), and the factorization from eq. (1) to see that an analytic form of the log-likelihood of the data emerges:
p(yk | y<k) = ∫ p(yk | xk) p(xk | y<k) dxk
= N ( yk |H x̂k|<k(y<k),H P̂k|<k(y<k)H⊤ +R ) ,
(13)
log p(y0:K) = K∑ k=0 logN ( yk |H x̂k|<k(y<k),H P̂k|<k(y<k)H⊤ +R ) . (14)
If we put ŷk|<k := H x̂k|<k(y<k) and M̂k|<k := H P̂k|<k(y<k)H⊤ + R, then the maximumlikelihood objective leads to the following loss function, which we can minimize using gradient descent methods w.r.t. all model parameters:
L := K∑
k=0
[ (ŷk|<k − yk)⊤ M̂−1k|<k (ŷk|<k − yk) + log detM̂k|<k ] . (15)
Note that each term in the sum above represents a one-step-ahead self-supervised error term. We thus minimize the prediction residuals ŷk|<k(y<k) − yk in a norm that is inversely scaled with the above covariance matrix, plus a regularizing determinant term, which prevents the covariance matrix from diverging. The arisen loss function is similar to the noise2noise (Lehtinen et al., 2018; Krull et al., 2019; Batson & Royer, 2019; Laine et al., 2019) objective from computer vision literature, combined with a locally linear transition model like Becker et al. (2019). We show in appendix D that this objective will yield correct results (meaning estimating the ground-truth xk) if the noise is independent with E[yk | Hxk] = Hxk. A similar procedure in the causality literature is given by Schölkopf et al. (2016). An algorithmic presentation of performing inference and fitting is presented in appendix E.
Note that after fitting the parameters to the data, eq. (6) can directly be used to do one-step ahead forecasting. Forecasting an arbitrary number of time steps is possible by plugging the new value x̂K+1|K via yK+1 := H x̂K+1|K back into the recurrent model, and so on. This is not a generative model but merely a convenience that we deemed worth mentioning.

heading: 6 LINEARIZED SMOOTHING
text: So far, we have only discussed how to perform filtering. Recall that for smoothing, we are instead interested in the quantity p(xk | y0:K). A smoothing strategy highly similar to the methods described earlier can be obtained by explicitly parameterizing such a model, which we detail in the supplementary material. Here, we introduce a linearized smoothing procedure. The essential advantage is that no additional model has to be trained, which can be costly. Several algorithms stemming from the Kalman filter literature can be applied, such as the RTS smoother (Rauch et al., 1965) and the two-filter smoother (Kitagawa, 1994). To enable this, we need to assume that the conditional mutual information I (xk−1; yk:K | xk,y0:k−1) is small for all k = 1, . . . ,K. In other words, we assume that we approximately have the following conditional independences:
xk−1 ⊥⊥ yk:K | (xk,y0:k−1). (16) To explain the motivation for this requirement, consider the model in fig. 2. If the states of y0:k−1 and xk are known, then the additional information that yk:K has about the latent variable xk−1 would need to be passed along the unblocked deeper paths like yk+1 ← xk+1 ← ek+1 ← ek ← xk−1. Then the assumption of small I(xk−1; yk:K | xk,y0:k−1) can be interpreted as that the deeper paths transport less information than the lower direct paths. If we consider all edges to the xk’s as linear and the edges to the ek’s as non-linear maps, the above could be interpreted as an information-theoretic version of expressing that the non-linear correction terms are small compared to the linear parts in the functional relations between the variables.
We will now show that under the earlier assumptions and eq. (16) we get a Gaussian approximation: p(xk | y0:K) ≈ N ( xk | ẑk, Ĝk ) . We will do backward induction with ẑK := x̂K|≤K and
ĜK := P̂K|≤K . To propagate this to previous time steps k − 1 we use the chain rule: p(xk−1 | y0:K) = ∫ p(xk−1 | xk,y0:K) p(xk | y0:K) dxk, (17)
where the second term is known by backward induction and for the first term we make use of the approximate conditional independence from eq. (16) to get
p(xk−1 | xk,y0:K) = p(xk−1 | yk:K ,xk,y0:k−1) eq. (16) ≈ p(xk−1 | xk,y0:k−1). (18)
The latter was shown to be Gaussian in section 4:
p(xk−1,xk | y0:k−1) = N ([ xk−1 xk ] | [ x̂k−1|≤k−1 x̂k|<k ] , [ P̂k−1|≤k−1 P̂k−1|≤k−1 F̂ ⊤ k|<k
F̂k|<k P̂k−1|≤k−1 P̂k|<k
]) .
(19)
By use of the usual formulas for Gaussians and the reverse Kalman gain matrix Ĵk−1|k we arrive at the following update formulas, k = K, . . . , 1, with ẑK := x̂K|≤K and ĜK := P̂K|≤K :
Ĵk−1|k := P̂k−1|≤k−1 F̂ ⊤ k|<k P̂ −1 k|<k, (20)
Ĝk−1 := P̂k−1|≤k−1 + Ĵk−1|k ( P̂k|≤k − P̂k|<k ) Ĵ⊤k−1|k, (21)
ẑk−1 := x̂k−1|≤k−1 + Ĵk−1|k (ẑk − x̂k|≤k). (22) As such, we can perform inference for all k = 0, . . . ,K with p(xk | y0:K) ≈ N ( xk | ẑk, Ĝk ) .
Algorithmically, the above is presented in appendix E.

heading: 7 RECURRENT NEURAL NETWORK
text: Before going into the experiments section, we briefly explain how we specifically estimate the functions F̂k|<k(y<k), êk|<k(y<k) and Q̂k|<k(y<k) that parameterize the transition probability p(xk | y<k,xk−1) (eq. (4)). The choice of the model here implicitly makes further assumptions about the generative model. If we consider neural networks, the temporal nature of the data suggests recurrent neural networks (Graves et al., 2013), convolutional neural networks (Kalchbrenner et al., 2014), or transformer architectures (Vaswani et al., 2017). Additionally, if the data is image-based, one might further make use of convolutions. For our experiments, we use a Gated Recurrent Unit (GRU) network (Cho et al., 2014), that recursively encodes hidden representations. Therefore, we put
Q̂k|<k := LkL ⊤ k ,
F̂kêk L̂k  := f(h→k ), h→k := GRU(y→k−1,h→k−1), (23) where L̂k is a Cholesky factor and f is a multi-layer perceptron decoder.

heading: 8 EXPERIMENTS
text: We perform three experiments, as motivated in section 1. Technical details on the setup of the experiments can be found in the supplementary material (appendix F). We refer to the model detailed in section 4 as the recursive filter, as it uses the Bayesian update recursion. For smoothing experiments, we use recursive smoother. The model obtained by parameterizing p(xk | y<k) directly (eq. (12)) is referred to as the recurrent filter or recurrent smoother, as it only employs recurrent neural networks (and no Bayesian recursion) to estimate said density directly.

heading: 8.1 LINEAR DYNAMICS
text: In the linear Gaussian case, it is known that the Kalman filter will give the optimal solution. Thus, we can get a lower bound on the test loss. In this toy experiment, we simulate particle tracking under linear dynamics and noisy measurements of the location. We use Newtonian physics equations as prior knowledge. We generate trajectories T = {x0:K ,y0:K} with xk ∈ R6 and yk ∈ R2 according to the differential equations:
ẋ = Ax = [ 0 1 0 0 −c 1 0 −τc 0 ][ p v a ] . (24)
We obtain sparse, noisy measurements yk = Hxk + r with r ∼ N (0,R). H is a selection matrix that returns a two-dimensional position vector. We run this experiment in a filtering setting, i.e.,
we only use past observations. We compare against (1) the raw, noisy measurements which inherently deviate from the clean measurements, (2) the Kalman filter solution where we optimized the transition covariance matrix using clean data (hence supervised), (3) the optimal solution, which is a Kalman filter with ground truth parameters performing exact inference. To estimate the transition maps for the Kalman filter, we use the standard Taylor series of eA·∆t up to the first order. Additionally, we use this expert knowledge as an inductive bias for the recursive filter’s transition maps.
In fig. 4 we depict the test mean squared error (MSE, lower is better) as a function of the number of training samples. Given enough data, the self-supervised models approximate the optimal solution arbitrarily well. Our recursive model significantly outperforms both the Kalman filter and the inference model in the low-data regime by using incorporated expert knowledge. Additionally, we report that the recursive model’s distance to the ground truth latent states is much closer to the optimal solution than both the inference model and Kalman filter. Specifically, we report average mean squared errors of 0.685 for the inference model, 0.241 for the Kalman filter, 0.161 for the recursive model compared to 0.135 for the optimal Kalman filter. Finally, it is worth noting that the recursive model has much less variance as a function of its initialization.

heading: 8.2 LORENZ EQUATIONS
text: We simulate a Lorenz system according to
ẋ = Ax = [ −σ σ 0 ρ− x1 −1 0 x2 0 −β ][ x1 x2 x3 ] . (25)
We have H = I , x ∈ R3 and y ∈ R3. The Lorenz equations model atmospheric convection and form a classic example of chaos. Therefore, performing inference is much more complex than in the linear case. This time, we perform smoothing (see appendix A) and compare against (1) the raw measurements, (2) a supervised Extended Kalman smoother (Ljung, 1979), (3) the variational inference approach of Krishnan et al. (2017), (4) the supervised recursive model of Satorras et al. (2019). Our models include a recurrent smoother, a recursive smoother, and the recursive filter with linearized smoothing (section 6). Transition maps for the Extended Kalman smoother and the recursive models are obtained by taking a second-order Taylor series of eA·∆t. For the supervised extended Kalman filter, we again optimize its covariance estimate using ground truth data.
In fig. 5 we plot the test mean squared error (MSE, lower is better) as a function of the number of examples available for training. It is clear that our methods approach the ground truth states with
more data. This is in contrast to the Extended Kalman smoother, which barely outperforms the noisy measurements. We also see that the recursive models significantly outperform the recurrent model in the low-data regime. The recursive filter with linearized smoothing performs comparably to the other models and even better in low-data regimes. We hypothesize that this is because the required assumption for the linearized smoother holds (section 6) and regularizes the model. The variational method of Krishnan et al. (2017) performs poorly in low-data regimes. Most notably, the supervised method of Satorras et al. (2019) outperforms our models only slightly.

heading: 8.3 AUDIO DENOISING
text: Next, we test the model on non-fabricated data with less ideal noise characteristics. Specifically, we use the SpeechCommands spoken audio dataset (Warden, 2018). Performing inference on spoken audio is challenging, as it arguably requires understanding natural language. To this end, recent progress on synthesizing raw audio has been made (Lakhotia et al., 2021). However, this requires scaling to much larger and more sophisticated neural networks than presented here, which we deem out the current work’s scope. Therefore, we take a subset of the entire dataset, using audio from the classes “tree”, “six”, “eight”, “yes”, and “cat”. We overlay these clean audio sequences {x(1)0:K , . . . ,x (N) 0:K} with various noise classes that the dataset provides. That is, for every noise class C we obtain a set of noisy sequences {y(1)0:K , . . . ,y (N) 0:K}C . We also consider a “combined” class in which we sample from the union of the noise sets. The task is to denoise the audio without having access to clean data. We evaluate the models on non-silent parts of the audio, as performance on those sections is the most interesting. Notably, none of these noise classes is Gaussian distributed.
We show the mean squared error on the test set of all models per noise class in table 2. Our models outperform the Kalman filter, Noise2Noise (Lehtinen et al., 2018), and SIN (Krishnan et al., 2017) unsupervised baselines. We suspect that the relatively poor performance of SIN is due to its generative Markov assumption, regularizing the model too strongly. The poor performance of Noise2Noise is due to the fact that it does not use the current measurement yk to infer xk. Like before, note that the Kalman filter is “supervised” as we optimized its covariance matrix using clean data. The supervised RKN (Becker et al., 2019) outperforms our models on most noise classes, but notably not on white noise. Most of these noise classes have temporal structure, making them predictable from past data. This is confirmed by observing these mean squared error values over the course of training. Initially, the values were better than reported in table 2, but the model increasingly fits the noise over time. Thus, although two of the main assumptions about the model (independent Gaussian noise) are violated, we still can denoise effectively. Since the RKN’s targets are denoised (hence “supervised”), it does not have this problem. In practice, obtaining clean data can be challenging.

heading: 9 CONCLUSION
text: We presented an advantageous parameterization of an inference procedure for nonlinear state-space models with potentially higher-order latent Markov chains. The inference distribution is split into linear and nonlinear parts, allowing for a recursion akin to the Kalman filter and smoother algorithms. Optimization is performed directly using a maximum-likelihood objective that backpropagates through these recursions. Smoothing can be performed similarly, but we additionally proposed linearized smoothing that can directly be applied to the filtering distributions. Our model is simple and builds on established methods from signal processing. Despite this, results showed that it can perform better or on par with fully supervised or variational inference methods.

heading: 10 ETHICS STATEMENT
text: The paper presents a simple method to perform inference using noisy sequential data. Applications can be found throughout society, e.g., tracking particles, denoising images or audio, or estimating system states. While many such examples are for good, there are applications with ethically debatable motivations. Among these could be tracking humans or denoising purposefully corrupted data (e.g., to hide one’s identity).

heading: 11 REPRODUCIBILITY STATEMENT
text: We are in the process of releasing code for the current work. For clarity and reproducibility, the presented methods are available as algorithms in the supplementary material. Furthermore, we made explicit wherever we needed to make an approximation or an assumption.

heading: A PARAMETERIZED SMOOTHING
text: In the main body of the paper, we showed how to parameterize the model for recursive estimation of p(xk | y≤k), k = 0, . . . ,K. Additionally, we provided a linearized smoothing procedure that yields p(xk | y1:K). The disadvantage clearly is the linearization. Here, we show how we can recursively estimate p(xk | y1:K) in a similar sense to the filtering setting. First, we put
p(xk | xk−1,y−k) = N ( xk | F̂k|−k xk−1 + êk|−k, Q̂k|−k ) , (26)
where F̂k|−k(y−k) and êk|−k(y−k) and Q̂k|−k(y−k) are two-sided recurrent neural network outputs similar to section 7. We compute the distribution of interest as follows:
p(xk | y0:K) = p(yk | xk)p(xk | y−k)
p(yk | y−k) (27) ∝ p(yk | xk) ∫ p(xk,xk−1 | y−k)dxk−1 (28)
= p(yk | xk) ∫ p(xk | xk−1,y−k)p(xk−1 | y−k)dxk−1 (29)
xk−1⊥⊥yk|y−k≈ p(yk | xk) ∫ p(xk | xk−1,y−k)p(xk−1 | y0:K)dxk−1 (30)
where we make the approximation to compute eq. (27) efficiently and recursively. It is justified if I(xk−1;yk | y−k) < ϵ for small ϵ. That is, the additional information that yk conveys about xk−1 is marginal if we have all other data. Let
p(xk−1 | y0:K) := N ( xk−1 | x̂k−1|0:K , P̂k−1|0:K ) (31)
be previous time-step’s posterior. Then put
p(xk | xk−1,y−k) := N (x | F̂k|−kxk−1 + êk|−k, Q̂k|−k). (32) where we left out the arguments for the following quantities estimated by an RNN.F̂k|−k(y−k)êk|−k(y−k)
Q̂k|−k(y−k)
 := f(h→k ,h←k ) [h→kh←k ] := [ GRU(h→k−1,yk−1) GRU(h←k+1,yk+1) ] (33)
Applying the integral in eq. (30) we get∫ p(xk | xk−1,y−k)p(xk−1 | y0:K)dxk−1 = N ( xk | x̂k|−k, P̂k|−k ) (34)
where we put
x̂k|−k := F̂k|−kx̂k−1|0:K + êk|−k P̂k|−k := F̂k|−kP̂k−1|0:KF̂ ⊤ k|−k + Q̂k|−k (35)
A data likelihood can be computed as follows p(yk | y−k) = ∫ p(yk | xk)p(xk | y−k)dxk (36)
= ∫ p(yk | xk) ∫ p(xk−1,xk | y−k)dxkxk−1 (37)
= ∫ p(yk | xk) ∫ p(xk | xk−1,y−k)p(xk−1 | y−k)dxkxk−1 (38)
xk−1⊥⊥yk|y−k≈ ∫ p(yk | xk) ∫ p(xk | xk−1,y−k)p(xk−1 | y0:K)dxkxk−1 (39)
where we made the same approximation eq. (30) as before. It evaluates to∫ p(yk | xk)p(xk | y−k)dxk = N ( yk | ŷk|−k,M̂k|−k ) (40)
with
ŷk|−k := Hx̂k|−k M̂k|−k := HP̂k|−kH ⊤ +R (41)
For fitting to the data we now would use a maximum-pseudo-likelihood (Gourieroux et al., 1984) approach by maximizing:
K∑ k=0 log p(yk | y−k) = K∑ k=0 logN ( yk |H x̂k|−k(y−k),H P̂k|−k(y−k)H⊤ +R ) , (42)
leading to minimizing the following self-supervised loss function:
L := K∑
k=0
[ (ŷk|−k − yk)⊤ M̂−1k|−k (ŷk|−k − yk) + log detM̂k|−k ] , (43)
where ŷk|−k := H x̂k|−k(y−k) and M̂k|−k := H P̂k|−k(y−k)H⊤ +R.

heading: B GAUSSIAN CONDITIONING FORMULAS
text: Since many of the calculations used in this work are based on the Gaussian conditioning formulas, we provide them here. If
p(x) = N (x | µ,P ), (44) p(y | x) = N (y |Hx+ b,R), (45)
then
p(y) = N (y |Hµ+ b,R+HPH⊤), (46) p(x | y) = N ( x | Σ [ H⊤R−1 (y − b) + P−1µ ] ,Σ ) , (47)
with
Σ = (P−1 +H⊤R−1H)−1. (48)
C DIRECT PARAMETERIZATION OF p(x | y<k)
We show here how to directly parameterize p(x | y<k) and p(x | y−k). This parameterization is referred to as the recurrent model in our experiments. The procedure is rather straightforward. For filtering, we put p(xk | y<k) = N ( xk | x̂k|−k(y<k), P̂k|<k(y<k) ) . We model:
x̂k|<k := ek, P̂k|<k := LkL ⊤ k ,
[ ek Lk ] := f(h→k ), (49)
where Lk is a cholesky factor and f is a multi-layer perceptron. The argument h→k is recursively given by
h→k := GRU(yk−1,h → k−1), (50)
where we employ a Gated Recurrent Unit (GRU) network (Cho et al., 2014). Note that this parameterization is equivalent to the model described in the main paper with Fk := 0.
For smoothing, p(xk | y−k) = N ( xk | x̂k|−k(y−k), P̂k|−k(y−k) ) .
x̂k|−k := ek, P̂k|−k := LkL ⊤ k ,
[ ek Lk ] := f(h→k ,h ← k ), (51)
where Lk is a cholesky factor and f is a multi-layer perceptron.
h→k := GRU(yk−1,h → k−1), h ← k := GRU(yk+1,h ← k+1). (52)
Once p(x | y<k) (filtering) or p(x | y−k) (smoothing) is obtained, all the procedures for inference and optimization described in the main paper and appendix A remain the same.

heading: D BIAS-VARIANCE-NOISE DECOMPOSITION OF THE SELF-SUPERVISED GENERALIZATION ERROR
text: Any estimate x̂k = x̂k(y−k) for xk that is not dependent on yk will give us a bias-variancenoise decomposition of the generalization error. Note that this setting covers both the filtering and smoothing case. Define “optimal model” ŷ∗k := E [yk | y−k], then under the specified generative model (section 3) we have
E [ ∥ŷk − yk∥22 | y−k ] = E [ ∥ŷk − ŷ∗k + ŷ∗k − yk∥ 2 2 | y−k ] (53)
= E [ ∥ŷk − ŷ∗k∥ 2 2 | y−k ] +Var [yk | y−k] + 2E [ (ŷk − ŷ∗k) ⊤ (ŷ∗k − yk) | y−k ] (54)
= ∥ŷk − ŷ∗k∥ 2 2 +Var [yk | y−k] + 2 (ŷk − ŷ ∗ k) ⊤ E [(ŷ∗k − yk) | y−k]︸ ︷︷ ︸
0
(55)
Var [yk | y−k] = E [ ∥ŷ∗k −Hxk − nk∥ 2 2 | y−k ] (56)
= E [ ∥ŷ∗k −Hxk∥ 2 2 | y−k ] + E [ ∥nk∥22 | y−k ] + 2E [ (Hxk − ŷ∗k) ⊤ nk | y−k ] ︸ ︷︷ ︸
0
(57)
xk,y−k⊥⊥nk = E [ ∥ŷ∗k −Hxk∥ 2 2 | y−k ] + tr(R) (58)
Thus,
E [ ∥ŷk − yk∥22 | y−k ] = ∥ŷk − ŷ∗k∥ 2 2 + E [ ∥ŷ∗k −Hxk∥ 2 2 | y−k ] + tr(R) (59)
Note that ŷ∗k = E [Hxk | y−k] = HE [xk | y−k]. Then, define our model ŷk := Hx̂k. The reducible part of the error becomes
∥ŷ∗k − ŷk∥ 2 2 = ∥H (E[xk | y−k]− x̂k)∥ 2 2 (60)
For this reason, the model output x̂k approaches the optimal model under minimization of the selfsupervised error (perturbed by H).
Additionally, the model ŷk approaches Hxk (the uncorrupted measurement) under this criterion. E [ ∥ŷk − yk∥2 ] = E [ ∥ŷk −Hxk +Hxk − yk∥2 ] (61)
= E [ ∥ŷk −Hxk∥2 ] + E [ ∥Hxk − yk∥2 ] + 2E [ (ŷk −Hxk)⊤ (Hxk − yk) ] (62)
= E [ ∥ŷk −Hxk∥2 ] + E [ ∥Hxk − yk∥2 ] + 2E [(ŷk −Hxk)]⊤ E [Hxk − yk]︸ ︷︷ ︸
0
(63)
= E [ ∥ŷk −Hxk∥2 ] + tr(R) + E [∥Hxk − yk∥]2︸ ︷︷ ︸
0
(64)
If we then take ŷk := Hx̂k then the reducible part of the error approximates the true xk (perturbed by H).
E [ ∥ŷk −Hxk∥2 ] = E [ ∥H (x̂k − xk)∥2 ] (65)

heading: E ALGORITHMS
text: Algorithm 1: Recursive Filter (Inference) input : Data (time-series) y0:K = (y0, . . . ,yK), emission matrices H and R, parameters ϕ. output:
1. For training: Loss value LK and its gradient ∇ϕLK .
2. For inference: x̂k|≤k and P̂k|≤k for k in 0, . . . ,K. Inference is done via: p(xk | y≤k) = N ( xk | x̂k|≤k, P̂k|≤k ) .
3. For linearized smoothing (section 6): F̂k|<k, P̂k|<k, x̂k|≤k and P̂k|≤k, for k in 0, . . . ,K.
4. For forecasting: x̂K+1|<(K+1) and P̂K+1|<(K+1). Forecasting is done via: p(xK+1 | y0:K) = N ( xK+1 | x̂K+1|<(K+1), P̂K+1|<(K+1) ) .
h0 := 0 L−1 := 0 P̂0|<0 := Q̂0 x̂0|<0 := ê0 for k in 0, . . . ,K do
B̂k|<k := ( H P̂k|<k H ⊤ +R )−1
ŷk|<k := H x̂k|<k
Lk := Lk−1 + (yk − ŷk|<k)⊤ B̂k|<k (yk − ŷk|<k)− log det B̂k|<k K̂k := P̂k|<k H ⊤B̂k|<k
P̂k|≤k := P̂k|<k − K̂k H P̂k|<k x̂k|≤k := x̂k|<k + K̂k (yk − ŷk|<k) hk+1 := GRUϕ(hk,yk) êk+1|<(k+1)F̂k+1|<(k+1)
L̂k+1|<(k+1)  := fϕ(hk+1) Q̂k+1|<(k+1) := L̂k+1|<(k+1) L̂ ⊤ k+1|<(k+1)
P̂k+1|<(k+1) := F̂k+1|<(k+1) P̂k|≤k F̂ ⊤ k+1|<(k+1) + Q̂k+1|<(k+1)
x̂k+1|<(k+1) := F̂k+1|<(k+1) x̂k|≤k + êk+1|<(k+1)
end For the training case we use backpropagation through the above loop to compute∇ϕLK .
Algorithm 2: Recursive Filter (Training) input : Data (time-series) y0:K = (y0, . . . ,yK), , emission matrices H and R, initialized parameters ϕ0, number of training rounds I . output: Model parameters ϕ∗ for inference at test-time. for i in 0, . . . , I do
Obtain L(i)K and∇ϕL (i) K from algorithm 1. Run preferred optimizer step to update parameters ϕ with∇ϕL(i)K (and L (i) K ).
end
Algorithm 3: Parameterized Recursive Smoother (Inference) input : Data (time-series) y0:K = (y0, . . . ,yK), emission matrices H and R, initialized parameters ϕ0 output:
1. Loss value LK and its gradient w.r.t. all model parameters∇ϕ(LK) for training.
2. For all k in 0, . . . ,K: x̂k|0:K , P̂k|0:K . These can be used for inference through p(xk | y0:K) = N ( x̂k|0:K , P̂k|0:K ) .
h→0 := 0 h←K+1 := 0 L(i)−1 := 0 P̂0|−0 := Q̂0 x̂0|−0 := ê0 for k in 0, . . . ,K do
B̂k|−k := ( H P̂k|−k H ⊤ +R )−1
ŷk|−k := H x̂k|−k
L(i)k := L (i) k−1 + (yk − ŷk|−k) ⊤ B̂k|−k (yk − ŷk|−k)− log det B̂k|−k K̂k := P̂k|−k H ⊤B̂k|−k
P̂k|−k := P̂k|−k − K̂k H P̂k|−k x̂k|−k := x̂k|−k + K̂k (yk − ŷk|−k) h→k+1 := GRUϕ(h → k ,yk
h←k+1 := GRUϕ(h ← k ,yk+1 êk+1|−(k+1)F̂k+1|−(k+1)
L̂k+1|−(k+1)  := fϕ(h→k+1,h←k+1) Q̂k+1|−(k+1) := L̂k+1|−(k+1) L̂ ⊤ k+1|−(k+1)
P̂k+1|−(k+1) := F̂k+1|−(k+1) P̂k|0:K F̂ ⊤ k+1|−(k+1) + Q̂k+1|−(k+1)
x̂k+1|−(k+1) := F̂k+1|−(k+1) x̂k|0:K + êk+1|−(k+1)
end For training case we also use backpropagation through the above loop to compute∇ϕLK
Algorithm 4: Parameterized Recursive Smoother (Training) input : Data (time-series) y0:K = (y0, . . . ,yK), emission matrices H and R, initialized parameters ϕ0, number of training rounds I output: Model parameters ϕ∗ for inference at test-time. for i in 1, . . . , I do
Obtain ∇ϕL(i)K from algorithm 3. Run preferred optimizer step.
end
Algorithm 5: Linearized Smoother
input : Values of: F̂k|<k, P̂k|<k, P̂k|≤k, x̂k|≤k for all k = 0, . . . ,K, obtained from the recursive filter algorithm. output: Linearly smoothed distributions p(xk | y0:K) = N ( xk | ẑk, Ĝk ) for all
k = 0, . . . ,K. ẑK := x̂K|≤K ĜK := P̂K|≤K for k = K, . . . , 1 do
Ĵk−1|k := P̂k−1|≤k−1 F̂ ⊤ k|<k P̂ −1 k|<k
Ĝk−1 := P̂k−1|≤k−1 + Ĵk−1|k ( P̂k|≤k − P̂k|<k ) Ĵ⊤k−1|k
ẑk−1 := x̂k−1|≤k−1 + Ĵk−1|k (ẑk − x̂k|≤k)
end
Algorithm 6: Recurrent Smoother (Training) input : Training data (time-series) y0:K = (y0, . . . ,yK), emission function N (Hxk,R), initialized parameters ϕ0, number of training iterations n. output: Optimized parameters ϕ∗ for i in 1 to n do
h0 := 0 hK+1 := 0 L(i) := 0 for k in 0 to K do
h→k := GRUϕ(h → k−1,yk−1)
h←k := GRUϕ(h ← k−1,yk+1)[
x̂k|−k L̂k|−k
] := fϕ(h ← k ,h → k )
P̂k|−k := L̂k|−kL̂ ⊤ k|−k
ŷk|−k := Hx̂k|−k B̂k|−k := ( HP̂k|−kH ⊤ +R )−1
L(i)k := L (i) k−1 + (yk − ŷk|−k) ⊤ B̂k|−k (yk − ŷk|−k)− log det B̂k|−k
end Compute ∇ϕL(i)K and apply SGD step with respect to all model parameters, which amounts
to backpropagation through the above calculations. end For the filter variant, the steps that involve the backward direction (←) are left out.
Algorithm 7: Recurrent Smoother (Inference) input : Test data y0:K = (y0, . . . ,yK), trained parameters ϕ∗, emission function N (Hx,R). output: Inferred posteriors p(xk | y0:K) = N (xk | x̂k|0:K , P̂k|0:K) for all k h0 := 0 hK+1 := 0 for k in 0 to K do
h→k := GRUϕ(h → k−1,yk−1)
h←k := GRUϕ(h ← k−1,yk+1)[
x̂k|−k L̂k|−k
] := fϕ(h ← k ,h → k )
P̂k|−k := L̂k|−kL̂ ⊤ k|−k K̂k := P̂k|−kH ⊤ ( HP̂k|−kH ⊤ +R )−1
x̂k|0:K := x̂k|−k + K̂k(yk −Hx̂k|−k) P̂k|0:K := P̂k|−k − K̂kHP̂k|−k
end For the filter variant, the steps that involve the backward direction (←) are left out.

heading: F EXPERIMENTS: DETAILS
text: F.1 LINEAR DYNAMICS
As specified in the main paper, the dynamics are according to
ẋ = Ax = [ 0 1 0 0 −c 1 0 −τc 0 ][ p v a ] . (66)
Since these are linear transitions, we can calculate any transition directly using x(t + ∆t) = eA∆tx(t).
F := [ eA 0 0 eA ] Q := [ Q̄ 0 0 Q̄ ] (67)
We used c = 0.06, τ = 0.17,∆t := 1 and covariance
Q̄ := 0.12 ·  13 0 00 1 0 0 0 3  (68) The matrix exponential is computed using Bader et al. (2019). The parameters for the emission distribution:
H := [ 1 0 0 0 0 0 0 0 0 1 0 0 ] R := 0.52 · [ 1 0 0 1 ] (69)
We simulate a K := 131, 072 trajectory for training, K := 16, 384 trajectory for validation and K := 32, 768 for testing. The F̃ that is used in the (non-optimal) Kalman filter and recursive model is computed as follows:
eA∆t ≈ F̃ := 1∑
n=0
(∆tAn)/n! (70)
F.2 LORENZ EQUATIONS
We simulate a Lorenz system according to
ẋ = Ax = [ −σ σ 0 ρ− x1 −1 0 x2 0 −β ][ x1 x2 x3 ] . (71)
We integrate the system using dt = 0.00001 and sample it uniformly at ∆t = 0.05. We use ρ = 28, σ = 10, β = 8/3. The transition in ∆t arbitrary time-steps is linearly approximated by a Taylor expansion and used in the Kalman smoother and recursive models.
eA|xk∆t ≈ F̃k := 2∑
n=0
(∆tA|xk) n/n! (72)
We simulate K := 131, 072 steps for training, K := 32, 768 for testing and K := 16, 384 for validation. We have H := I and thus x ∈ R3 and y ∈ R3. We use R := 0.52I .

heading: G PARAMETERIZED SMOOTHING: ALTERNATIVE POSTERIOR EVALUATION
text: We would like to point out that the distribution p(xk | y0:K) can be obtained without making assumption eq. (30), which we stretch out here. Initial experiments showed that using these calculations the model did not converge as smoothly as when using the ones stated before. However, it could be of interest to further investigate. Returning to the posterior of interest
p(xk | y0:K) = ∫ p(xk | xk−1,y0:K) p(xk−1 | y0:K) dxk−1 (73)
=
∫ p(yk | xk)
p(yk | xk−1,y−k) p(xk | xk−1,y−k)p(xk−1 | y0:K) dxk−1. (74)
We have
p(xk−1 | y0:K) = N ( xk−1 | x̂k−1|0:K(y0:K), P̂k−1|0:K(y0:K) ) (75)
as the previous time-step’s posterior. p (xk | xk−1,y−k) = N ( xk | F̂k|−kxk−1 + êk|−k, Q̂k|−k ) (76)
where F̂k|−k(y−k), êk|−k(y−k) and Q̂k|−k(y−k) are estimated by a neural network. Combining this with noise model eq. (3) we get:
p(xk | xk−1,y0:K) = N ( xk | F̂k|−k xk−1 + êk|−k +Kk (yk −H (F̂k|−k xk−1 + êk|−k)), Q̂k|−k −Kk H Q̂k|−k ) (77)
= N ( xk | (I −Kn H) ( F̂k|−k xk−1 + êk|−k ) +Kk yk, (I −Kn H) Q̂k|−k ) , (78)
where we directly applied the Woodbury matrix identity to obtain Kalman gain matrix Kk := Q̂k|−k H ⊤ ( H Q̂k|−k H ⊤ +R )−1
(79)
Then, applying the integral we get: p (xk | y0:K) = N ( xk | x̂k|0:K , P̂k|0:K ) , (80)
with:
P̂k|0:K = (I −Kn H) F̂k|−k P̂k−1|0:K F̂⊤k|−k (I −Kn H) ⊤ + (I −Kn H) Q̂k|−k, (81)
and
x̂k|0:K = F̂k|−kx̂k−1|0:K + êk|−k +Kk (yk −H (F̂k|−kx̂k−1|0:K + êk|−k)) (82) = (I −Kn H) ( F̂k|−k x̂k−1|0:K + êk|−k ) +Kk yk. (83)

Review summaries:
A Kalman filter typically requires specifying a state update equation, along with noise covariance matrices for state transition and observation. Determining these can be a challenge in certain scenarios. The manuscript proposes to replace these unknowns with trainable neural networks. The input to the neural networks are the observations, and the output are the unknown parameters (or their factorizations, to ensure positive definiteness). The manuscript validates the approach through numerical experiments.

The authors describe a parametrized inference approach for nonlinear dynamical models with linear observations. They borrow ideas from the well-known Kalman filter, and use Bayesian filtering and smoothing approaches to recursively estimate the states of the nonlinear model. The parameterization makes use of a locally linear transition between subsequent states using neural networks; thus the authors are still able to perform the Kalman-like filtering and smoothing.

The authors proposed a new inference method for state-space model with nonlinear latent states. They showed that their inference model can be incorporated with domain knowledge and optimized in a self-supervised manner. They also demonstrated that their model can achieve competitive results through both simulated and real audio denoising tasks.