# Efficient Online Variational Estimation via Monte Carlo Sampling

Mathis Chagneux, Mathias Müller, Pierre Gloaguen, Sylvain Le Corff, Jimmy Olsson

► To cite this version:

Mathis Chagneux, Mathias Müller, Pierre Gloaguen, Sylvain Le Corff, Jimmy Olsson. Efficient Online Variational Estimation via Monte Carlo Sampling. 2026. hal-05492599

HAL Id: hal-05492599

https://hal.science/hal-05492599v1

Preprint submitted on 6 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Efficient Online Variational Estimation via Monte Carlo Sampling

Mathis Chagneux†, Mathias Muller‡, Pierre Gloaguen*, Sylvain Le Corff⊤, and Jimmy Olsson‡

†Télécom Paris, Institut Polytechnique de Paris, Palaiseau.

‡KTH Royal Institute of Technology.

*Université Bretagne Sud, UMR CNRS 6205, LMBA, F-56000 Vannes, France.

⊤Sorbonne Université, Université Paris Cité, CNRS, Laboratoire de Probabilités, Statistique et Modélisation, LPSM, F-75005 Paris, France.

## Abstract

This article addresses online variational estimation in parametric state-space models. We propose a new procedure for efficiently computing the evidence lower bound and its gradient in a streaming-data setting, where observations arrive sequentially. The algorithm allows for the simultaneous training of the model parameters and the distribution of the latent states given the observations. It is based on i.i.d. Monte Carlo sampling, coupled with a well-chosen deep architecture, enabling both computational efficiency and flexibility. The performance of the method is illustrated on both synthetic data and real-world air-quality data. The proposed approach is theoretically motivated by the existence of an asymptotic contrast function and the ergodicity of the underlying Markov chain, and applies more generally to the computation of additive expectations under posterior distributions in state-space models.

## 1 Introduction

This work considers state-space models (SSMs) where the law of observations $(Y_t)_{t \in \mathbb{N}}$ is governed by an unobserved, or ‘hidden’, Markov chain $(X_t)_{t \in \mathbb{N}}$, and the finite dimensional distributions of $(X_t, Y_t)_{t \in \mathbb{N}}$ are given by parametric distributions indexed by $\theta \in \Theta$. Learning the parameter $\theta$ in this context is a complex task, as it usually requires access to the joint-smoothing distributions, defined at time $t$ as the conditional distribution of the latent states $X_{0:t}^1$ given the corresponding observations $Y_{0:t}$. This paper addresses the challenging problem of online estimation in SSMs, which consists of sequentially learning both $\theta$ and the smoothing distributions as data stream in real time.

A first possible approach to this task is based on sequential Monte Carlo (SMC) smoothing (see [Chopin et al., 2020] and the references therein). Such SMC-based methods for online learning come with strong theoretical guarantees [Le Corff and Fort, 2013, Olsson and Westerborn, 2017, Tadić and Doucet, 2020, Gao et al., 2025], but typically suffer from the curse of dimensionality in the state dimension. This recently motivated the use of variational inference (VI) [Blei et al., 2017] as an alternative, as VI has demonstrated strong empirical performance for online inference in large-scale settings, particularly through stochastic variational frameworks [Hoffman et al., 2013, Broderick et al., 2013].

Taking a variational approach, the smoothing distributions are approximated by simpler distributions, the variational distributions, depending on some unknown parameter $\phi \in \Phi$. Both $\theta$ and $\phi$ are then learned by maximizing a proxy of the log likelihood, the evidence lower bound (ELBO), which is generally done using gradient ascent, thus requiring the computation of the ELBO’s gradient. A key feature of these approaches is the choice of the variational distribution, which must be suited for online learning. [Campbell et al., 2021] rely on structured variational distributions that both mimic the Markovian form of the smoothing distributions (see Section 3) and are well suited to online learning. The authors explore online variational additive smoothing for the recursive computation of the ELBO and its gradients, using a Bellman-type recursion inspired by reinforcement learning and recursive maximum likelihood methods. A major limitation of this algorithm is that, at each iteration $t$, it requires solving an inner optimization problem to learn a regression function serving as a proxy for a conditional expectation. This step can be computationally expensive and depends critically on the appropriate choice of the regression class. More broadly, a key drawback of existing online variational learning methods is the lack of theoretical guarantees for the proposed algorithms.

<footnote n="1">¹$\mathbf{a}_{u:v}$ is a short-hand notation for $(a_u, \dots, a_v)$.</footnote>
In this paper, we develop a theoretically grounded framework for online variational learning. Our approach is rooted in a stochastic approximation perspective, and, following [Mastrototaro et al., 2025], we show that for a structured variational distributions parameterized by flexible function classes (such as deep neural networks), online variational learning amounts to maximizing a lower bound (COLBO) on the limiting time-normalized asymptotic log likelihood, also known as the asymptotic contrast function. Building on this theoretical framework, we first formulate an ideal—but generally intractable—stochastic approximation algorithm that maximizes the COLBO, and then propose a Monte Carlo version of then same, using an efficient importance sampling approach that avoids any regression task and outperforms the algorithm of [Campbell et al., 2021] in terms of computation time. The main contributions of this paper can be summarized as follows:

* We provide a theoretically grounded framework for online variational learning, showing that it can be viewed as Robbins-Monro algorithm.

* We propose a computationally efficient online estimator of the COLBO its gradient in the context of SSMs. In contrast to computationally intensive SMC or Markov chain Monte Carlo (MCMC) methods, our algorithm, which we refer to as *Recursive Monte Carlo Variational Inference* (RMCVI) to emphasize its iterative structure, relies on simple i.i.d. samples from the marginal variational distributions.

* The proposed algorithm is not limited to the online computation and optimization of the COLBO and can be directly adapted to compute more general expectations, training losses, or gradients under distributions that admit a Markovian structure.

* Experimentally, we demonstrate the performance of our estimator both on synthetic and real world data.

## 2 Related work

Our methodology draws on recent advances in smoothing methods for SSMs by (i) proposing a Monte Carlo approach for approximating conditional expectations and (ii) relying on a structured variational family.

*SMC for online learning in SSMs.* The original approaches to online smoothing is based on SMC methods. We refer the reader to [Douc et al., 2014, Section 11] for a presentation of the general concepts underlying SMC algorithms, and to [Olsson and Westerborn, 2017, Gloaguen et al., 2022, Dau and Chopin, 2022] for more recent developments and applications of these methods to online smoothing. Theoretical guarantees for online optimization of the log-likelihood using SMC can be found in [Le Corff and Fort, 2013, Olsson and Westerborn, 2017, Tadić and Doucet, 2020, Gao et al., 2025].

*Offline VI for SSMs.* Early works on VI for SSMs primarily focus on offline estimation [Johnson et al., 2016, Krishnan et al., 2017, Lin et al., 2018], i.e., they require prior access to the entire observation sequence $Y_{0:t}$ in order to compute the gradients of the ELBO. These approaches rely on a *forward* factorization of the variational distribution, which is incompatible with online learning.

*Online VI for SSMs.* [Marino et al., 2018, Zhao and Park, 2020, Dowling et al., 2023] opts to trade smoothing for filtering by targeting the marginal distributions at each timestep with variational distributions that depend only on observations up to $t$. By designing a new variational family, [Campbell et al., 2021] show that the ELBO and its gradient can be recursively expressed via conditional expectations, providing a natural framework for online learning. The authors approximate these conditional expectations at each time step using functional approximations.

*Online Variational SMC.* A recent approach to improving SMC for online learning is to adapt the particle proposal dynamics by optimizing a variational objective [Zhao et al., 2022, Mastrototaro and Olsson, 2024, Mastrototaro et al., 2025]. In this framework, the variational method aims to approximate the locally optimal proposal while simultaneously learning the model parameters, rather than targeting the full joint smoothing distribution.

*Theoretical guarantees for VI.* In the offline setting, [Chagneux et al., 2024] established the first theoretical results on error control in VI for SSM, building on the variational family of [Campbell et al., 2021]. Moreover, [Mastrototaro et al., 2025] provided an online variational SMC framework that provably maximizes a well-defined asymptotic contrast function via stochastic approximation, which serves as the conceptual inspiration for our work.

## 3 Model and background

Consider an SSM $(X_t, Y_t)_{t \in \mathbb{N}}$, where $(X_t)_{t \in \mathbb{N}}$ is a discrete-time Markov chain on $\mathcal{X} := \mathbb{R}^{d_x}$. The distribution of $X_0$ has density $\chi$ with respect to the Lebesgue measure $\mu$ and for all $t \in \mathbb{N}$, the conditional distribution of $X_{t+1}$ given $X_{0:t}$ depends only on $X_t$ and has transition density $m_t^\theta(X_t, \cdot)$. In SSMs, it is assumed that the states of the Markov
chain are only partially observed through an observation process $(Y_t)_{t \in \mathbb{N}}$ taking on values in $Y := \mathbb{R}^{d_y}$. For every $t \in \mathbb{N}$, the observations $Y_{0:t}$ are assumed to be conditionally independent given $X_{0:t}$ and such that the conditional distribution of each $Y_s$, $s \in [0, t]$, given $X_{0:t}$ depends only on the corresponding $X_t$ and has density $g_t^\theta(X_t, \cdot)$ with respect to the Lebesgue measure. The model is then defined, for all $t \in \mathbb{N}$, by the joint distributions

$$p_{0:t}^{\theta}(x_{0:t}, y_{0:t}) := \chi(x_0)g_0^{\theta}(x_0, y_0) \prod_{s=1}^{t} m_s^{\theta}(x_{s-1}, x_s)g_s^{\theta}(x_s, y_s) \quad (1)$$

of the hidden states and the observations.

A classical learning task in SSMs is state inference, which consists of estimating the joint-smoothing distribution, i.e. the conditional distribution of $X_{0:t}$ given $Y_{0:t}$, given by

$$\pi_{0:t}^{\theta}(x_{0:t}) := \frac{p_{0:t}^{\theta}(x_{0:t}, y_{0:t})}{p_{0:t}^{\theta}(y_{0:t})},$$

where $p_{0:t}^{\theta}(y_{0:t}) := \int p_{0:t}^{\theta}(x_{0:t}, y_{0:t}) dx_{0:t}$ is the observed-data likelihood. The marginal of this joint distribution with respect to the state $x_t$ at time $t$ is known as the filtering distribution at time $t$, and its density with respect to the Lebesgue measure is denoted by $\pi_t^\theta$. It is straightforward to show that the density of the joint-smoothing distribution satisfies the so-called backward decomposition

$$\pi_{0:t}^{\theta}(x_{0:t}) = \pi_t^{\theta}(x_t) \prod_{s=1}^{t} b_{s-1|s}^{\theta}(x_s, x_{s-1}), \quad (2)$$

where each backward kernel

$$b_{s-1|s}^{\theta}(x_s, x_{s-1}) \propto m_s^{\theta}(x_{s-1}, x_s)\pi_{s-1}^{\theta}(x_{s-1}) \quad (3)$$

is the conditional probability density function of $X_{s-1}$ given $(X_s, Y_{0:s-1})$. The backward decomposition stems from the fact that the hidden process is still Markov when evolving conditionally on the observations, with time-inhomogeneous transition densities (3). However, since the filtering distributions are intractable, the backward kernels generally lack closed-form expressions.

In variational approaches, the smoothing distribution $\pi_t^\theta$ is approximated by selecting a candidate from a parametric family $\{q_{0:t}^\phi\}_{\phi \in \Phi}$, known as the variational family, where $\Phi$ is a parameter space. This parameter is learned jointly with $\theta$ by maximizing the ELBO

$$\mathcal{L}_t^{\theta, \phi} = \mathbb{E}_{q_{0:t}^{\phi}} \left[ \log \frac{p_{0:t}^{\theta}(X_{0:t}, Y_{0:t})}{q_{0:t}^{\phi}(X_{0:t})} \right], \quad (4)$$

where $\mathbb{E}_{q_{0:t}^{\phi}}$ denotes expectation under $q_{0:t}^\phi$. Needless to say, the form of the variational family is crucial in this approach. Motivated by the backward decomposition, some works impose structure on the variational family through a factorization of $q_{0:t}^\phi$. A variational counterpart of (2), introduced by [Campbell et al., 2021], is given by

$$q_{0:t}^{\phi}(x_{0:t}) = q_t^{\phi}(x_t) \prod_{s=1}^{t} q_{s-1|s}^{\phi}(x_s, x_{s-1}), \quad (5)$$

where $q_t^\phi$ (resp. $q_{s-1|s}^\phi(x_s, \cdot)$) are user-designed probability density functions, the variational kernels, whose parameters are learned from data. A decisive advantage of this factorization is that it respects the true dependencies in (2). Moreover, [Chagneux et al., 2024] established an upper bound on the error when expectations with respect to the smoothing distribution are approximated by expectations with respect to variational distributions that satisfy this backward factorization.

**Defining the variational distributions.** The variational inference framework requires the definition of the variational distributions involved in (5), i.e. the set of distributions $(q_t^\phi, q_{t|t+1}^\phi)_{t \in \mathbb{N}}$. We here define a new variational family in a recursive manner, with shared parameters over time. This new design based on backward factorization is efficient in terms of online parameter learning (as the number of parameters does not grow with $t$) and creates a link between variational kernels that will ensure an efficient importance sampling procedure in Section 5.

Our online learning challenge requires that variational distributions (i) can be recursively defined using streaming data $(y_t)_{t \in \mathbb{N}}$, (ii) are related to each other to mimic the relation given by (3) between the backward kernel and the
filtering distribution, (iii) are easy to sample from to perform Monte Carlo approximations. For this purpose, each $q_t^\phi$ is chosen as a parametric distribution belonging to the exponential family (in our experiments, the Gaussian family in $\mathcal{X}$), defined by some parameter $\eta_t$ belonging to a parameter space $\mathcal{E}$. More precisely, we define intermediate quantities $(a_t)_{t \in \mathbb{N}}$ belonging to some user-defined space $\mathcal{A}$, initialized at some arbitrary value $a_0 \in \mathcal{A}$ and governed by a deterministic recursion $a_t = \mathcal{A}^\phi(a_{t-1}, y_t)$. Based on these quantities, we let, for each $t$, $\eta_t = f^\phi(a_t)$. Here the mappings $\mathcal{A}^\phi$ and $f^\phi$ are used-defined. This framework creates a link between variational filtering distributions, in the spirit of the filtering recursions in SSMs. The variational backward kernels are then defined from on the basis of this flow of distributions by setting, for all $t \in \mathbb{N}_{>0}$,

$$q_{t-1|t}^{\phi}(x_t, x_{t-1}) \propto q_{t-1}^{\phi}(x_{t-1}) \psi_t^{\phi}(x_{t-1}, x_t), \quad (6)$$

where $(\psi_t^\phi)_{t \ge 0}$ are potential functions on $\mathcal{X}^2$ of form $\psi_t^\phi(x_{t-1}, x_t) = \exp(\langle \tilde{\eta}_t^\phi(x_t), T(x_{t-1}) \rangle)$, with $\tilde{\eta}_t^\phi(x_t) = \tilde{f}^\phi(x_t)$ and $T(x_{t-1})$ being a natural parameter and a sufficient statistic, respectively, for the chosen exponential family. Eqn. (6) ensures that $q_{t-1|t}^\phi(x_t, \cdot)$ will be a probability density function with natural parameter $\eta_{t-1|t}^\phi = \eta_{t-1}^\phi + \tilde{\eta}_t^\phi$. In this convenient setting, the backward kernels $q_{t-1|t}^\phi$ can have arbitrarily complex dependencies on $x_t$, while their densities are derived analytically from the potentials. This enables straightforward Monte Carlo sampling procedures and direct computations of normalizing constants (which are required in our proposed algorithm, e.g., in (16) below), while at the same time avoiding the reduction of our variational kernels to mere transformations or linearizations (e.g., linear Gaussian kernels). In is important to note that the parameters of functions $\mathcal{A}^\phi, f^\phi, \tilde{f}^\phi$ are shared across time, leading to an amortized framework. In our experiments, these functions are neural networks, and $\phi$ are their weights.

# 4 Online variational learning

**The asymptotic contrast function and the COLBO.** In the context of maximum likelihood estimation, the online learning of an unknown model parameter $\theta$ is known as recursive maximum likelihood (RML) [Le Gland and Mevel, 1997]. RML focuses on maximizing the asymptotic contrast function $\lambda(\theta) := \lim_{t \to \infty} t^{-1} \log p_\theta(Y_{0:t})$ (a.s.), which serves as a foundational objective in this setting. If the data are generated by an SSM belonging to the parametric family of interest, characterized by a 'true' parameter $\theta^*$, then, under suitable identifiability conditions, the asymptotic contrast is maximised at $\theta^*$. Consequently, the maximum likelihood estimator (MLE) is strongly consistent in the sense that it converges almost surely to $\theta^*$ as $t$ tends to infinity.

Since the asymptotic contrast $\lambda(\theta)$ is intractable, we use a similar approach to that of [Mastrototaro et al., 2025] and instead aim to maximise online, with respect to $(\theta, \phi)$, the contrast lower bound (COLBO) given by

$$\ell(\theta, \phi) := \lim_{t \to \infty} \frac{1}{t} \mathcal{L}_t^{\theta, \phi} \le \lambda(\theta) \quad (\text{a.s.}). \quad (7)$$

**Stochastic approximation viewpoint.** Following standard RML ideas, online variational learning seek to maximize $\ell(\theta, \phi)$ by updating $(\theta, \phi)$ in the direction of its gradient. Because $\ell(\theta, \phi)$ is defined as a long-run time average, it is natural to pursue a stochastic approximation approach with the goal of solving $\nabla_{\theta,\phi} \ell(\theta, \phi) = 0$. Indeed, defining $\mathcal{G}_t^{\theta, \phi} := \nabla_{\theta,\phi} \mathcal{L}_t^{\theta, \phi} - \nabla_{\theta,\phi} \mathcal{L}_{t-1}^{\theta, \phi}$, we may write

$$\frac{1}{t} \nabla_{\theta, \phi} \mathcal{L}_t^{\theta, \phi} = \frac{1}{t} \sum_{s=1}^{t} \mathcal{G}_s^{\theta, \phi} + \frac{1}{t} \nabla_{\theta, \phi} \mathcal{L}_0^{\theta, \phi}. \quad (8)$$

Interpreted through the lens of ergodic theory, it is tempting to see the long term limit of the right-hand side of (8) as an expectation, allowing $\nabla_{\theta,\phi} \ell(\theta, \phi)$ to be expressed as a mean field (i.e., the deterministic drift) that governs the long-run behavior of Robbins-Monro stochastic updates. In this idealized framework, it is natural to use the observed gradient increments to build a sequence $(\theta_t, \phi_t)_{t \in \mathbb{N}}$ leading to an ideal procedure summarized in Algorithm 1. However, since each term $\mathcal{G}_s^{\theta, \phi}$ depends on the whole historical record $Y_{0:s}$ and $(Y_t)_{t \in \mathbb{N}}$ is not a Markov process, the existence of the limit of (8) as $t$ tends to infinity is non-trivial. Actually, to the best of our knowledge, no theoretical results exist justifying the existence of the COLBO objective (7) and its gradient, both of which are necessary to place this learning procedure on firm theoretical ground. In the coming sections we provide theoretical results motivating this existence for the variational family of Section 3 (Eqn. (5) and (6)). In particular, we show that the COLBO objective and its gradient can be justified via the law of large numbers, applied to a suitably constructed Markov chain. On the basis of this justification, Algorithm 1 can be motivated as a stochastic approximation scheme with state-dependent Markov noise.
**Recursive expression of the ELBO and its gradient.** In the following, we assume that we are given a sequence $(y_t)_{t \in \mathbb{N}}$ of observations, and leave the dependence on these implicit in the notation. Write $\ell_0^{\theta, \phi}(x_{-1}, x_0) := \log(\chi(x_0)g_0^{\theta}(x_0, y_0))$, $q_{-1|0}^{\phi}(x_0, x_{-1}) = 1$ and for $t \in \mathbb{N}_{>0}$,

$$ \ell_t^{\theta, \phi}(x_{t-1}, x_t) := \log \left( \frac{m_t^{\theta}(x_{t-1}, x_t)g_t^{\theta}(x_t, y_t)}{q_{t-1|t}^{\phi}(x_t, x_{t-1})} \right), \quad (9) $$

which allows to rewrite the ELBO (4) as

$$ \mathcal{L}_t^{\theta, \phi} = \mathbb{E}_{q_{0:t}^{\phi}} \left[ \sum_{s=0}^{t} \ell_s^{\theta, \phi}(x_{s-1}, x_s) - \log q_t^{\phi}(x_t) \right]. $$

The following proposition provides recursive formulas for computing both the ELBO and its gradient. For brevity, we let $\mathbb{E}_{t-1|t}^{\phi, x_t}$ denote expectation under $q_{t-1|t}^{\phi}(x_t, \cdot)$.

**Proposition 4.1.** For every $t \in \mathbb{N}$ and $(\theta, \phi) \in \Theta \times \Phi$, the ELBO and its gradient are given by

$$ \begin{aligned} \mathcal{L}_t^{\theta, \phi} &= \mathbb{E}_{q_t^{\phi}}[h_t(X_t)] - \mathbb{E}_{q_t^{\phi}}[\log q_t^{\phi}(X_t)], \\ \nabla_{\phi} \mathcal{L}_t^{\theta, \phi} &= \mathbb{E}_{q_t^{\phi}}[\nabla_{\phi} \log q_t^{\phi}(X_t) h_t(X_t) + u_t(X_t)], \\ \nabla_{\theta} \mathcal{L}_t^{\theta, \phi} &= \mathbb{E}_{q_t^{\phi}}[v_t(X_t)], \end{aligned} $$

where the real-valued function $h_t$ on $\mathcal{X}$ and its gradients $u_t := \nabla_{\phi}h_t$ and $v_t := \nabla_{\theta}h_t$ satisfy the recursions

$$ \begin{aligned} h_t(x_t) &= \mathbb{E}_{t-1|t}^{\phi, x_t} \left[ h_{t-1}(X_{t-1}) + \ell_t^{\theta, \phi}(X_{t-1}, x_t) \right] \\ u_t(x_t) &= \mathbb{E}_{t-1|t}^{\phi, x_t} \left[ u_{t-1}(X_{t-1}) + \nabla_{\phi} \log q_{t-1|t}^{\phi}(x_t, X_{t-1}) \{h_{t-1}(x_{t-1}) + \ell_t^{\theta, \phi}(x_{t-1}, x_t)\} \right] \\ v_t(x_t) &= \mathbb{E}_{t-1|t}^{\phi, x_t} \left[ v_{t-1}(X_{t-1}) + \nabla_{\theta} \ell_t^{\theta, \phi}(X_{t-1}, x_t) \right], \end{aligned} $$

with $h_0(x_0) = \ell_0^{\theta, \phi}(x_{-1}, x_0)$, $u_0(x_0) = 0$, and $v_0(x_0) = \nabla_{\theta}\ell_0^{\theta, \phi}(x_{-1}, x_0)$.

*Proof.* See Section B. $\square$

Proposition 4.1 is of twofold interest. First, it provides a recursive scheme for the online computation of the ELBO, which will lead to a learning algorithm in Section 5. Second, it highlights the natural quantities to consider when studying the existence of the COLBO objective and its gradient.

**Existence of the COLBO.** From now on, we assume that the observed data is generated by some SSM $(X_t, Y_t)_{t \in \mathbb{N}}$, which does not necessarily belong to the parametric family considered in Section 3. Under this assumption, it is easy to see that also the process $(Z_t)_{t \in \mathbb{N}}$, where $Z_t := (X_t, Y_t, h_t, u_t, v_t, a_t)$, with $h_t, u_t$, and $v_t$ being the functions defined recursively in Theorem 4.1 and $a_t$ being the intermediate quantities used in the parameterization of $q_t^{\phi}$, is a Markov chain. The state space and Markov kernel of $(Z_t)_{t \in \mathbb{N}}$ are denoted by $(\mathcal{Z}, \mathcal{Z})$ and $T^{\theta, \phi}$, respectively (see Section C,

**Algorithm 1** Ideal algorithm (exact recursions)

1: For each $t \in \mathbb{N}$, at the arrival of $y_t$, using Prop. 4.1:

2: Compute $h_t, u_t, v_t$.

3: Compute $\nabla_{\theta} \mathcal{L}_t^{\theta_t, \phi_t}$ and $\nabla_{\phi} \mathcal{L}_t^{\theta_t, \phi_t}$.

4: Update:

$$ \begin{aligned} \phi_{t+1} &\leftarrow \phi_t + \gamma_{t+1}^{\phi} \left( \nabla_{\phi} \mathcal{L}_t^{\theta_t, \phi_t} - \nabla_{\phi} \mathcal{L}_{t-1}^{\theta_t, \phi_t} \right), \\ \theta_{t+1} &\leftarrow \theta_t + \gamma_{t+1}^{\theta} \left( \nabla_{\theta} \mathcal{L}_t^{\theta_t, \phi_t} - \nabla_{\theta} \mathcal{L}_{t-1}^{\theta_t, \phi_t} \right), \end{aligned} $$

where $(\gamma_{t+1}^{\theta}, \gamma_{t+1}^{\phi})$ are learning rates satisfying the usual Robbins Monro conditions.
Eqn. (34), for details). The Markov property follows from the assumed SSM dynamics of $(X_t, Y_t)_{t \in \mathbb{N}}$, along with the fact that the updates of Proposition 4.1, as well as the update of $a_t$ from $a_{t-1}$, are performed recursively based on the current observation $Y_t$. Denote also by $S^\phi$ the Markov kernel of the marginal chain $(X_t, Y_t, a_t)_{t \in \mathbb{N}}$. We will establish the exponential forgetting of the extended chain $(Z_t)_{t \in \mathbb{N}}$ under the following assumptions.

**Assumption 4.2.** There exist $\pi \in \mathbf{M}_1(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ and $\alpha \in (0,1)$ such that for every $t \in \mathbb{N}_{>0}$, $\phi \in \Phi$, and $(x, y, a) \in \mathbf{X} \times \mathbf{Y} \times \mathcal{A}$,

$$ \| (S^\phi)^t(x, y, a) - \pi\|_{TV} \le \alpha^t. $$

Theorem 4.2 is discussed in Section C.2. The following assumptions are purely technical.

**Assumption 4.3.** There exists $c \in \mathbb{R}_{>0}$ such that for every $\phi \in \Phi$, $a_t \in \mathcal{A}$, and $(x_{s+1}, y_{s+1}) \in \mathbf{X} \times \mathbf{Y}$,

$$ (i) \quad \mathbb{E}_{s|s+1}^{\phi, x_{s+1}} \left[ |\nabla_{\phi} \log q_{s|s+1}^{\phi}(x_{s+1}, X_s)|^2 \right] \le c^2, $$

$$ (ii) \quad \mathbb{E}_{s|s+1}^{\phi, x_{s+1}} \left[ |\ell_s^{\theta, \phi}(X_s, x_{s+1})|^2 \right] \le c^2, $$

$$ (iii) \quad \mathbb{E}_{s|s+1}^{\phi, x_{s+1}} \left[ |\nabla_{\theta} \ell_s^{\theta, \phi}(X_s, x_{s+1})| \right] \le c. $$

**Assumption 4.4.** There exist constants $0 < \varepsilon^- < \varepsilon^+ < \infty$ such that, for every $t \in \mathbb{N}_{>0}$, $(x_{t-1}, x_t) \in \mathbf{X}^2$, and $\phi \in \Phi$,

$$ \varepsilon^- \le \psi_t^\phi(x_{t-1}, x_t) \le \varepsilon^+. $$

The following theorem establishes the geometric ergodicity of the extended Markov chain, if not for all $f$ in the space $\mathbb{F}_b(\mathcal{Z})$ of bounded measurable functions on $\mathbb{Z}$, so at least for a subclass $\mathcal{L}(\mathcal{Z}) \subset \mathbb{F}_b(\mathcal{Z})$ of Lipschitz functions. More precisely, $f \in \mathcal{L}(\mathcal{Z})$ if there exists $\varphi \in \mathbb{F}_b(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ such that for every $(x, y, a, h, h', u, u', v, v') \in \mathbf{X} \times \mathbf{Y} \times \mathcal{A} \times \mathbb{F}_b(\mathcal{X})^6$, writing $\mathbf{s} = (x, y, a)$,

$$ (i) \quad |f(\mathbf{s}, h, u, v)| \le \varphi(\mathbf{s}), $$

$$ (ii) \quad |f(\mathbf{s}, h, u, v) - f(\mathbf{s}, h', u', v')| \le \varphi(\mathbf{s}) (\mathrm{osc}(h - h') + \mathrm{osc}(u - u') + \mathrm{osc}(v - v')). $$

**Theorem 4.5 (Geometric ergodicity of $(Z_t)_{t \in \mathbb{N}}$).** Assume 4.2–4.4. Then there exist $\rho \in (0,1)$ and a functional $\zeta : \mathbb{F}_b(\mathcal{X})^6 \to \mathbb{R}_{>0}$ such that for every $(\theta, \phi) \in \Theta \times \Phi$, $t \in \mathbb{N}$, $f \in \mathcal{L}(\mathcal{Z})$, $z = (x_0, y_0, a_0, h_0, u_0, v_0) \in \mathbb{Z}$, and $z' = (x'_0, y'_0, a'_0, h'_0, u'_0, v'_0) \in \mathbb{Z}$,

$$ |(T^{\theta, \phi})^t f(z) - (T^{\theta, \phi})^t f(z')| \le \zeta(h_0, h'_0, u_0, u'_0, v_0, v'_0) \|\varphi\|_{\infty} \rho^t. \quad (10) $$

Moreover, there exists a functional $\bar{\zeta} : \mathbb{F}_b(\mathcal{X})^3 \to \mathbb{R}_{>0}$ and a kernel $\Pi^{\theta, \phi}$ such that for every $f \in \mathcal{L}(\mathcal{Z})$, $\Pi^{\theta, \phi} f$ is constant and for every $t \in \mathbb{N}$ and $z = (x_0, y_0, a_0, h_0, u_0, v_0) \in \mathbb{Z}$,

$$ |(T^{\theta, \phi})^t f(z) - \Pi^{\theta, \phi} f| \le \bar{\zeta}(h_0, u_0, v_0) \|\varphi\|_{\infty} \rho^t. \quad (11) $$

*Proof.* See Section C.

Although the contraction (11) does not hold in the total variation norm (due to the restriction to test functions in $\mathcal{L}(\mathcal{Z})$), the quantity $\Pi^{\theta, \phi}$ provided by the same theorem can be regarded as a candidate for the unique stationary distribution of $(Z_t)_{t \in \mathbb{N}}$. Moreover, by Theorem 4.1, each term $\mathcal{G}_s^{\theta, \phi} = \mathcal{G}^{\theta, \phi}(Z_{s-1:s})$ depends explicitly on the consecutive states $Z_{s-1:s}$ of the extended chain. By the law of large numbers for Markov chains, we may expect that (a.s.),

$$ \lim_{t \to \infty} \frac{1}{t} \nabla_{\theta, \phi} \mathcal{L}_t^{\theta, \phi} = \iint \mathcal{G}^{\theta, \phi} \langle z, z' \rangle \Pi^{\theta, \phi}(dz) T^{\theta, \phi}(z, dz'). \quad (12) $$

Letting the limit (12) serve as the mean field of a stochastic approximation scheme with state-dependent Markov noise [Karimi et al., 2019], a recursive Robbins–Monro algorithm finding a stationary point of the COLBO gradient is given by

$$ (\theta_{t+1}, \phi_{t+1}) \leftarrow (\theta_t, \phi_t) + \gamma_{t+1} \mathcal{G}^{\theta_t, \phi_t} (Z_{t-1}, Z_t), $$

and $Z_{t+1} \sim T^{\theta_{t+1}, \phi_{t+1}}(Z_t, \cdot)$, where $(\gamma_t)_{t \in \mathbb{N}_{>0}}$ is a sequence of step sizes satisfying the usual assumptions. This procedure is summarized in Algorithm 1, which uses distinct step-size sequences $(\gamma_t^\theta)_{t \in \mathbb{N}_{>0}}$ and $(\gamma_t^\phi)_{t \in \mathbb{N}_{>0}}$ for updating the model and variational parameters.
# 5 Online Monte Carlo approximation

We now derive a practical version of the ideal Algorithm 1. Theorem 4.1 suggests that it is possible to estimate the ELBO and its gradient recursively. The key feature of our Monte Carlo algorithm is that each conditional expectation in the recursion only needs to be estimated on a finite support, bypassing the regression step required at each time step in [Campbell et al., 2021]. This results in a more efficient procedure, as confirmed empirically in Section 6.2. The algorithm proceeds as follows:

First, sample $\{\xi_0^i\}_{i=1}^N \stackrel{\text{i.i.d.}}{\sim} q_0^\phi$, and set

$$ \hat{h}_0^{\phi,i} = h_0(\xi_0^i), \quad \hat{u}_0^{\phi,i} = u_0(\xi_0^i), \quad \hat{v}_0^{\theta,i} = v_0(\xi_0^i). $$

At time $t \in \mathbb{N}_{>0}$, having access to a Monte Carlo sample $\{\xi_{t-1}^i\}_{i=1}^N$ from $q_{t-1}^\phi$ and approximations $\hat{h}_{t-1}^{\phi,i}$, $\hat{u}_{t-1}^{\phi,i}$, $\hat{v}_{t-1}^{\theta,i}$ of $h_{t-1}(\xi_{t-1}^i)$, $u_{t-1}(\xi_{t-1}^i)$, $v_{t-1}(\xi_{t-1}^i)$, respectively, sample independently $\{\xi_t^i\}_{i=1}^N$ from $q_t^\phi$ and update

$$ \hat{h}_t^{\phi,i} = \sum_{j=1}^{N} \bar{w}_{t-1|t}^{\phi,i,j} (\hat{h}_{t-1}^{\phi,j} + \ell_t^{\theta,\phi}(\xi_{t-1}^j, \xi_t^i)), \quad (13) $$

$$ \hat{u}_t^{\phi,i} = \sum_{i=1}^{N} \bar{w}_{t-1|t}^{\phi,i,j} \left\{ \hat{u}_{t-1}^{\phi,j} + \nabla_{\phi} \log q_{t-1|t}^{\phi}(\xi_t^i, \xi_{t-1}^j) (\hat{h}_{t-1}^{\phi,j} + \ell_t^{\theta,\phi}(\xi_{t-1}^j, \xi_t^i)) \right\}, \quad (14) $$

$$ \hat{v}_t^{\theta,i} = \sum_{i=1}^{N} \bar{w}_{t-1|t}^{\phi,i,j} \left( \hat{v}_{t-1}^{\theta,j} + \nabla_{\theta} \ell_t^{\theta,\phi}(\xi_{t-1}^j, \xi_t^i) \right), \quad (15) $$

where

$$ \bar{w}_{t-1|t}^{\phi,i,j} := \frac{q_{t-1|t}^{\phi}(\xi_t^i, \xi_{t-1}^j)/q_{t-1}^{\phi}(\xi_{t-1}^j)}{\sum_{k=1}^{N} q_{t-1|t}^{\phi}(\xi_t^i, \xi_{t-1}^k)/q_{t-1}^{\phi}(\xi_{t-1}^k)}. \quad (16) $$

Estimators (13–15) are self-normalized importance sampling estimators of the updates of Proposition 4.1, and (16) provides the (shared) importance weights of these estimators. Note that we cannot perform direct Monte Carlo approximation on the basis of samples from $q_{t-1|t}^\phi(\xi_t^i, \cdot)$, as we would not have access to any approximations of the values of the functionals $h_{t-1}$ and $u_{t-1}$ at the sampled points. The use of importance sampling is therefore a prerequisite for updating the approximations. Moreover, note that the design of variational distributions imposed by (6) creates a link between the target $q_{t-1|t}^\phi$ distribution and $q_t^\phi$, making the latter a natural proposal distribution.

Once $\{(\hat{h}_t^{\phi,i}, \hat{u}_t^{\phi,i}, \hat{v}_t^{\theta,i})\}_{i=1}^N$ are computed, approximations of the ELBO and its gradient at time $t$ are obtained by

$$ \hat{\mathcal{L}}_t^{\theta,\phi} = \frac{1}{N} \sum_{i=1}^{N} (\hat{h}_t^{\phi,i} - \log q_t^\phi(\xi_t^i)), \quad (17) $$

$$ \hat{\nabla}_{\phi} \mathcal{L}_t^{\theta,\phi} = \frac{1}{N} \sum_{i=1}^{N} \left( \nabla_{\phi} \log q_t^{\phi}(\xi_t^i) \hat{h}_t^{\phi,i} + \hat{u}_t^{\phi,i} \right), \quad (18) $$

$$ \hat{\nabla}_{\theta} \mathcal{L}_t^{\theta,\phi} = \frac{1}{N} \sum_{i=1}^{N} \hat{v}_t^{\theta,i}. \quad (19) $$

The full procedure, which we refer to as RMCVI (Recursive Monte Carlo Variational Inference), is detailed in Algorithm 2 in Section A. Appendix D provides refinements to increase computational efficiency and reduce the variance of the gradient estimator. It is worth noting that RMCVI computes estimated gradients using both $(\phi_t, \theta_t)$ and $(\phi_{t-1}, \theta_{t-1})$ (Eqns. (24–25), Appendix A), which introduces a deviation from the ideal updates of Algorithm 1. Such approximations—commonly employed in RML settings—are crucial for enabling a feasible practical online implementation.

# 6 Experiments

We now evaluate the proposed algorithm on several streaming data inference tasks. Our goal is to demonstrate that the method can jointly learn both the latent posterior approximation and the model parameters in a fully sequential manner while requiring substantially less computation than existing online approaches (Table 1).
<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Evolution of $\hat{\mathcal{L}}_t^{\theta, \phi}/t$ during the online learning in the linear Gaussian SSM. Vertical lines indicate times at which state estimation is made on a test sequence. (b) State estimation on a test sequence (only one particular dimension is displayed).</figcaption></figure>

In all experiments, the variational filtering distributions $q_t^\phi$ are chosen to be in the Gaussian family. For the non-linear models (Sections 6.2 and 6.3) we implement the deterministic recursion $a_t = \mathcal{A}^\phi(a_{t-1}, y_t)$ as an RNN, where both the update function $\mathcal{A}^\phi$ and the parameter mapping $f^\phi$ are parameterized by MLPs with tanh activation functions. The variational backward kernels are defined via the potentials $\psi_t^\phi$ in (6), which are parameterized by similar neural networks (except for the linear Gaussian case, where we exploit analytical conjugation to derive exact backward kernels). In addition to this architecture, some control-variate tricks are implemented to reduce the variance of the gradient estimator (see details in Appendix E and D).

## 6.1 Linear-Gaussian HMM

We first assess our algorithm on a linear Gaussian SSM. Consider $(X_t, Y_t)_{t \in \mathbb{N}}$ in $\mathbb{R}^{d_x} \times \mathbb{R}^{d_y}$, with $X_0 \sim \mathcal{N}(\mu_0, Q_0)$, $Y_0 = GX_0 + \varepsilon_0$, and, for all $t \in \mathbb{N}_{>0}$,

$$ 
\begin{aligned}
X_t &= F X_{t-1} + \nu_t, & \nu_t &\sim \mathcal{N}(0, I_{d_x}), \\
Y_t &= G X_t + \varepsilon_t, & \varepsilon_t &\sim \mathcal{N}(0, I_{d_y}),
\end{aligned}
$$

where $\theta = (F, G)$ is the parameter to be learned. In this case, the true smoothing distributions are Gaussian and can be computed via closed-form recursions (the Kalman smoother), providing an analytical reference for evaluation. The variational family $q_{0:t}^\phi$ is parameterized using Gaussian conditionals and marginals; see Section E.1 for details.

The first experiment is run with $d_x = d_y = 10$, and a sequence of $T = 50000$ observations. Figure 1 shows the evolution of the ELBO through the learning of $(\theta, \phi)$ as well as the posterior mean of the hidden states on a test sequence of observations never seen by the model. The right panel shows a particular dimension of the hidden state for the test sequence, as well as the learned posterior distribution at different iterations. As the number of observations grows, our estimator learns a mapping that produces the true posterior distribution, corresponding to the oracle Kalman smoother. Figure 2 shows the MAE in parameter estimation of $F$ (a) and $G$ (b). Our method achieves lower error than the regression based method of [Campbell et al., 2021] in roughly one quarter of the runtime when using $N = 1000$ samples for the importance weights. Overall, these results highlight the computational efficiency of the proposed updates without compromising statistical accuracy.

## 6.2 Chaotic recurrent neural network.

We now consider the model used in [Campbell et al., 2021], where $X_0 \sim \mathcal{N}(0, Q)$, $Y_0 = X_0 + \varepsilon_0$, and, for $t \in \mathbb{N}_{>0}$,

$$ 
\begin{aligned}
X_t &= X_{t-1} + \frac{\Delta}{\rho} (\gamma W \tanh(X_{t-1}) - X_{t-1}) + \eta_t, \\
Y_t &= X_t + \varepsilon_t,
\end{aligned}
$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Model-parameter learning in the linear-Gaussian HMM. Mean-absolute errors of the transition (F) and emission (G) matrices for our method and [Campbell et al., 2021].</figcaption></figure>

where $(\eta_t)_{t \in \mathbb{N}_{>0}}$ and $(\varepsilon_t)_{t \in \mathbb{N}}$ are mutually independent sequences of i.i.d. $\mathcal{N}(0, Q)$ and Student-t random variables, respectively. We set $d_x = d_y = 5$, and the use same true parameters as [Campbell et al., 2021] (see Appendix E.2) and $N = 500$ importance samples.

**Filtering and one-step smoothing.** To compare our method with the one of [Campbell et al., 2021] on this model, we reproduce the 1-step smoothing experiment of their work [Campbell et al., 2021, Appendix B.2]. Specifically, we evaluate the ability of both approaches to estimate the conditional laws of $X_{t-1}$ given $Y_{0:t}$ and of $X_t$ given $Y_{0:t}$ by learning $\mathbb{E}_{q_{t-1:t}^{\phi}}[X_{t-1}]$ and $\mathbb{E}_{q_t^{\phi}}[X_t]$. In order to perform the same comparison, we mimic the non-amortized framework of the original paper (details are provided Appendix E.2). Table 1 reports the 1-step smoothing and filtering errors as defined in Eqn. (61). We also report average computation time per gradient steps. With comparable errors, RMCVI is about 5 times faster than the regression approach.

**Online learning.** Moving beyond the fixed-parameter setting, we evaluate our method in a true streaming regime where both the variational parameters and selected generative parameters are learned online. Concretely, $T = 5 \times 10^5$ observations are processed with updates of the parameters $\gamma$ and $\rho$. Figure 3 (a) shows the MAE between parameter estimates and the true value. The central panel shows the dynamics of $\hat{\mathcal{L}}_t^{\theta, \phi}/t$ with checkpoint markers, at which state estimation on held-out sequences is performed (right panel) for a specific state dimension. These results show that the proposed scheme remains stable and accurate while simultaneously learning $(\rho, \gamma)$ in this highly nonlinear regime.

## 6.3 Air-quality Data

We evaluate the framework on the UCI Air-Quality dataset [De Vito, 2008], which consists of hourly averaged responses from a chemical sensor array alongside meteorological data. We process the data into an 8-dimensional observation vector $y_t \in \mathbb{R}^8$ spanning approximately one year ($T \approx 9,300$ steps). This benchmark is characterized by frequent periods of sensor failure, where valid signals are absent for extended durations. Rather than imputing these irregularities offline, we process the stream directly to rigorously test the method's ability to maintain coherent belief states during blackout periods. Visualizations of the data and more details on the signal characteristics and preprocessing are provided in Appendix E.3.

<table><caption>Table 1: Time per gradient step, 1-step smoothing and filtering RMSE ($\times 10^{-2}$) (defined in Eqn. (61)) for the chaotic RNN.</caption><thead><tr><th>Method</th><th>1-Smooth.</th><th>Filt.</th><th>Time</th></tr></thead><tbody><tr><td>RMCVI (ours)</td><td>8.9 (0.2)</td><td>10.3 (0.2)</td><td>1 ms</td></tr><tr><td>[Campbell et al., 2021]</td><td>9.2 (0.2)</td><td>10.3 (0.2)</td><td>4.8 ms</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Parameter MAE. (b) Approximate $\tilde{\mathcal{L}}_t/t$ with checkpoint markers. (c) State estimation on a test sequence for one latent dimension. Colored lines/markers correspond to the same checkpoints.</figcaption></figure>

We model the air quality dynamics using a non-linear Gaussian SSM with residual transitions. Let $X_t \in \mathbb{R}^5$ and $Y_t \in \mathbb{R}^8$ denote the latent state and observations respectively. The generative process is defined as:

$$ 
\begin{aligned}
X_t &= X_{t-1} + f_{\theta}(X_{t-1}) + \nu_t, & \nu_t &\sim \mathcal{N}(0, Q_{\theta}) \\
Y_t &= g_{\theta}(X_t) + \varepsilon_t, & \varepsilon_t &\sim \mathcal{N}(0, R_{\theta})
\end{aligned}
$$

where $f_{\theta}$ and $g_{\theta}$ are neural networks parameterized by $\theta$ (with tanh activations), and the noise terms have diagonal covariance matrices, which are also learned. Here we used $N = 20$ importance samples for the Monte Carlo estimates.

**Online predictive performance.** Our primary focus is the model's performance in an online setting where all parameters must be learned from a cold start. We evaluate the model's ability to learn complex dynamics by measuring the one-step-ahead prediction RMSE on the five primary pollutants (CO, NO$_x$, NO$_2$, C$_6$H$_6$, O$_3$). We compare our RMCVI method against two baselines: a probabilistic online LSTM [Salinas et al., 2020] with a Gaussian output trained sequentially via maximum likelihood, and Online Variational SMC (OVSMC, [Mastrototaro and Olsson, 2024]).

Figure 4 demonstrates the robustness and accuracy of our approach. As shown in Panel (a), we introduce a smoothing experiment with artificial sensor failure. While the purely autoregressive LSTM is limited to filtering and thus tracks the corrupted signal, RMCVI effectively recovers the underlying ground truth. Similarly, panel (b) confirms that RMCVI matches the predictive performance of the LSTM and outperforms OVSMC. Thus, RMCVI combines the forecasting power of autoregressive networks with the advantage of variational smoothing, all while maintaining lower computational costs than particle-based methods.

# 7 Conclusion

We introduced a theoretically grounded online variational learning algorithm for SSM. The performance of our method is assessed with synthetic and real-world datasets, where it is shown to be more efficient than recent alternatives both for smoothing and prediction tasks. An important future work concerns the theoretical analysis of the online Monte Carlo version of the algorithm, which is crucial to obtain quantitative guarantees and optimize hyperparameters.

# References

[Blei et al., 2017] Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017). Variational inference: A review for statisticians. *Journal of the American Statistical Association*, 112(518):859–877.

[Broderick et al., 2013] Broderick, T., Boyd, N., Wibisono, A., Wilson, A. C., and Jordan, M. I. (2013). Streaming variational bayes. *Advances in neural information processing systems*, 26.

[Campbell et al., 2021] Campbell, A., Shi, Y., Rainforth, T., and Doucet, A. (2021). Online variational filtering and parameter learning. *Advances in Neural Information Processing Systems*, 34.
[Chagneux et al., 2024] Chagneux, M., Gassiat, É., Gloaguen, P., and Le Corff, S. (2024). Additive smoothing error in backward variational inference for general state-space models. *Journal of Machine Learning Research*.

[Chopin et al., 2020] Chopin, N., Papaspiliopoulos, O., et al. (2020). *An introduction to sequential Monte Carlo*, volume 4. Springer.

[Dau and Chopin, 2022] Dau, H.-D. and Chopin, N. (2022). On the complexity of backward smoothing algorithms. *arXiv preprint arXiv:2207.00976*.

[De Vito, 2008] De Vito, S. (2008). *Air Quality*. UCI Machine Learning Repository. DOI: https://doi.org/10.24432/C59K5F.

[Douc et al., 2014] Douc, R., Moulines, E., and Stoffer, D. (2014). *Nonlinear time series: theory, methods and applications with R examples*. CRC Press.

[Dowling et al., 2023] Dowling, M., Zhao, Y., and Park, I. M. (2023). Real-time variational method for learning neural trajectory and its dynamics. In *The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023*. OpenReview.net.

[Gao et al., 2025] Gao, Y., Li, W., and Chen, R. (2025). Parameter estimation in state space models using particle importance sampling. *Proceedings of Machine Learning Research*, 258:1252–1260.

[Gloaguen et al., 2022] Gloaguen, P., Le Corff, S., and Olsson, J. (2022). A pseudo-marginal sequential Monte Carlo online smoothing algorithm. *Bernoulli*, 28(4):2606 – 2633.

[Hoffman et al., 2013] Hoffman, M. D., Blei, D. M., Wang, C., and Paisley, J. (2013). Stochastic variational inference. *the Journal of machine Learning research*, 14(1):1303–1347.

[Johnson et al., 2016] Johnson, M. J., Duvenaud, D. K., Wiltschko, A., Adams, R. P., and Datta, S. R. (2016). Composing graphical models with neural networks for structured representations and fast inference. *Advances in neural information processing systems (NeurIPS)*, 29.

[Karimi et al., 2019] Karimi, B., Miasojedow, B., Moulines, E., and Wai, H.-T. (2019). Non-asymptotic analysis of biased stochastic approximation scheme. *Proceedings of Machine Learning Research*, 99:1–31.

[Kingma and Ba, 2015] Kingma, D. P. and Ba, J. (2015). Adam: A method for stochastic optimization. In *3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings*.

[Krishnan et al., 2017] Krishnan, R., Shalit, U., and Sontag, D. (2017). Structured inference networks for nonlinear state space models. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 30.

[Le Corff and Fort, 2013] Le Corff, S. and Fort, G. (2013). Convergence of a particle-based approximation of the block online expectation maximization algorithm. *ACM Transactions on Modeling and Computer Simulation (TOMACS)*, 23(1):1–22.

[Le Gland and Mevel, 1997] Le Gland, F. and Mevel, L. (1997). Recursive estimation in HMMs. In *Proc. IEEE Conf. Decis. Control*, pages 3468–3473.

[Lin et al., 2018] Lin, W., Khan, M. E., and Hubacher, N. (2018). Variational message passing with structured inference networks. In *International Conference on Learning Representations*.

[Marino et al., 2018] Marino, J., Cvitkovic, M., and Yue, Y. (2018). A general method for amortizing variational filtering. In *Advances in neural information processing systems (NeurIPS)*, volume 31.

[Mastrototaro et al., 2025] Mastrototaro, A., Müller, M., and Olsson, J. (2025). Recursive learning of asymptotic variational objectives. *International Conference on Artificial Intelligence and Statistics (AISTATS)*.

[Mastrototaro and Olsson, 2024] Mastrototaro, A. and Olsson, J. (2024). Online variational sequential monte carlo. *Proceedings of the 41st International Conference on Machine Learning (ICML)*.

[Mohamed et al., 2020] Mohamed, S., Rosca, M., Figurnov, M., and Mnih, A. (2020). Monte carlo gradient estimation in machine learning. *The Journal of Machine Learning Research*, 21(1):5183–5244.
[Olsson and Westerborn, 2017] Olsson, J. and Westerborn, J. (2017). Efficient particle-based online smoothing in general hidden Markov models: the PaRIS algorithm. *Bernoulli*, 23(3):1951–1996.

[Salinas et al., 2020] Salinas, D., Flunkert, V., Gasthaus, J., and Januschowski, T. (2020). Deepar: Probabilistic forecasting with autoregressive recurrent networks. *International Journal of Forecasting*, 36(3):1181–1191.

[Sun et al., 2020] Sun, C., Song, M., Hong, S., and Li, H. (2020). A review of designs and applications of echo state networks. *arXiv preprint arXiv:2012.02974*.

[Tadić and Doucet, 2005] Tadić, V. B. and Doucet, A. (2005). Exponential forgetting and geometric ergodicity for optimal filtering in general state-space models. *Stochastic Processes and their Applications*, 115.

[Tadić and Doucet, 2020] Tadić, V. B. and Doucet, A. (2020). Asymptotic properties of recursive particle maximum likelihood estimation. *IEEE Transactions on Information Theory*, 67(3):1825–1848.

[Yoshida and Miyato, 2017] Yoshida, Y. and Miyato, T. (2017). Spectral norm regularization for improving the generalizability of deep learning. *arXiv preprint arXiv:1705.10941*.

[Zhao et al., 2022] Zhao, Y., Nassar, J., Jordan, I., Bugallo, M., and Park, I. M. (2022). Streaming variational monte carlo. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 45(1):1150–1161.

[Zhao and Park, 2020] Zhao, Y. and Park, I. M. (2020). Variational online learning of neural dynamics. *Frontiers in computational neuroscience*, 14:71.
# A Full online gradient estimator

<caption>**Algorithm 2** Online gradient estimator</caption>

**Input:** initial parameter estimate $(\theta_0, \phi_0)$, step sizes $\{\gamma_t^\theta, \gamma_t^\phi\}_{t \ge 1}$.

Sample $\{\xi_0^i\}_{i=1}^N$ independently from $q_0^{\phi_0}$

For $1 \le i \le n$, set $\hat{h}_0^{\phi_0, i} \leftarrow h_0(\xi_0^i)$, $\hat{u}_0^{\phi_0, i} \leftarrow 0$, $\hat{v}_0^{\phi_0, i} \leftarrow v_0(\xi_0^i)$

set $\hat{\nabla}_\phi \mathcal{L}_0^{\theta_0, \phi_0} \leftarrow N^{-1} \sum_{i=1}^N \nabla \log q_0^{\phi_0}(\xi_0^i) \hat{h}_0^{\phi_0, i}$ and $\hat{\nabla}_\theta \mathcal{L}_0^{\theta_0, \phi_0} \leftarrow N^{-1} \sum_{i=1}^N \hat{v}_0^{\phi_0, i}$

update $\phi_1 \leftarrow \phi_0 + \gamma_1^\phi \hat{\nabla}_\phi \mathcal{L}_0^{\theta_0, \phi_0}$ and $\theta_1 \leftarrow \theta_0 + \gamma_1^\theta \hat{\nabla}_\theta \mathcal{L}_0^{\theta_0, \phi_0}$

**for** $t \leftarrow 1$ **to** $n-1$ **do**

sample $\{\xi_t^i\}_{i=1}^N$ independently from $q_t^{\phi_t}$

set

$$ \bar{w}_{t-1|t}^{\phi_t, i, j} \leftarrow \frac{q_{t-1|t}^{\phi_t}(\xi_t^i, \xi_{t-1}^j)/q_{t-1|t}^{\phi_{t-1}}(\xi_{t-1}^j)}{\sum_{k=1}^{N} q_{t-1|t}^{\phi_t}(\xi_t^i, \xi_{t-1}^k)/q_{t-1|t}^{\phi_{t-1}}(\xi_{t-1}^k)}; \quad (20) $$

update

$$ \hat{h}_{t}^{\phi_{t}, i} \leftarrow \sum_{j=1}^{N} \bar{w}_{t-1|t}^{\phi_{t}, i, j} \left( \hat{h}_{t-1}^{\phi_{t-1}, j} + \ell_{t}^{\theta, \phi_{t}}(\xi_{t-1}^{j}, \xi_{t}^{i}) \right), \quad (21) $$

$$ \hat{u}_{t}^{\phi_{t}, i} \leftarrow \sum_{i=1}^{N} \bar{w}_{t-1|t}^{\phi_{t}, i, j} \left\{ \hat{u}_{t-1}^{\phi_{t-1}, j} + \nabla_{\phi} \log q_{t-1|t}^{\phi_{t}}(\xi_{t}^{i}, \xi_{t-1}^{j}) \left( \hat{h}_{t-1}^{\phi_{t-1}, j} + \ell_{t}^{\theta, \phi_{t}}(\xi_{t-1}^{j}, \xi_{t}^{i}) \right) \right\}, \quad (22) $$

$$ \hat{v}_{t}^{\theta_{t}, i} \leftarrow \sum_{i=1}^{N} \bar{w}_{t-1|t}^{\phi_{t}, i, j} \left\{ \hat{v}_{t-1}^{\theta_{t-1}, j} + \nabla_{\theta} \ell_{t}^{\theta, \phi_{t}}(\xi_{t-1}^{j}, \xi_{t}^{i}) \right\}; \quad (23) $$

set

$$ \hat{\nabla}_{\phi} \mathcal{L}_{t}^{\theta_{t}, \phi_{t}} = \frac{1}{N} \sum_{i=1}^{N} \left\{ \nabla_{\phi} \log q_{t}^{\phi_{t}}(\xi_{t}^{i}) \hat{h}_{t}^{\phi_{t}, i} + \hat{u}_{t}^{\phi_{t}, i} \right\}, \\ \hat{\nabla}_{\theta} \mathcal{L}_{t}^{\theta_{t}, \phi_{t}} = \frac{1}{N} \sum_{i=1}^{N} \hat{v}_{t}^{\theta_{t}, i}. $$

Update

$$ \phi_{t+1} \leftarrow \phi_t + \gamma_{t+1}^{\phi} \left( \hat{\nabla}_{\phi} \mathcal{L}_{t}^{\theta_t, \phi_t} - \hat{\nabla}_{\phi} \mathcal{L}_{t-1}^{\theta_{t-1}, \phi_{t-1}} \right) \quad (24) $$

$$ \theta_{t+1} \leftarrow \theta_t + \gamma_{t+1}^{\theta} \left( \hat{\nabla}_{\theta} \mathcal{L}_{t}^{\theta_t, \phi_t} - \hat{\nabla}_{\theta} \mathcal{L}_{t-1}^{\theta_{t-1}, \phi_{t-1}} \right); \quad (25) $$

end for

# B Proof of proposition 4.1

By definition,

$$ \mathcal{L}_{t}^{\theta, \phi} = \mathbb{E}_{q_{0:t}^{\phi}} \left[ \log \frac{p_{0:t}^{\theta}(X_{0:t}, Y_{0:t})}{q_{0:t}^{\phi}(X_{0:t})} \right] = \mathbb{E}_{q_{0:t}^{\phi}} \left[ \log \frac{\chi(X_0) g_0^{\theta}(X_0, Y_0) \prod_{s=1}^{t} m_s^{\theta}(X_{s-1}, X_s) g_s^{\theta}(X_s, Y_s)}{q_t^{\phi}(X_t) \prod_{s=1}^{t} q_{s-1|s}^{\phi}(X_s, X_{s-1})} \right], $$

which, using (9), yields

$$ \mathcal{L}_{t}^{\theta, \phi} = \mathbb{E}_{q_{0:t}^{\phi}} \left[ \sum_{s=0}^{t} \ell_{s}^{\theta, \phi}(X_{s-1}, X_s) - \log q_{t}^{\phi}(X_t) \right]. $$

Therefore,

$$ \mathcal{L}_{t}^{\theta, \phi} = \mathbb{E}_{q_{t}^{\phi}}[h_{t}(X_{t})] - \mathbb{E}_{q_{t}^{\phi}}[\log q_{t}^{\phi}(X_{t})], \quad (26) $$
where

$$h_t(x_t) = \int \left( \sum_{s=0}^{t} \ell_s^{\theta, \phi}(x_{s-1}, x_s) \right) \prod_{s=1}^{t} q_{s-1|s}^{\phi}(x_s, x_{s-1}) \mathrm{d}x_{0:t-1} = \mathbb{E}_{q_{0:t-1|t}^{\phi}(x_{t,\cdot})} \left[ \ell_{0:t}^{\theta, \phi}(X_{0:t-1}, x_t) \right],$$

with $q_{0:t-1|t}^{\phi}(x_t, x_{0:t-1}) = \prod_{s=1}^{t} q_{s-1|s}^{\phi}(x_s, x_{s-1})$, and we introduced the notation $\ell_{0:t}^{\theta, \phi}(x_{0:t}) = \sum_{s=0}^{t} \ell_s^{\theta, \phi}(x_{s-1}, x_s)$. Then, note that

$$h_t(x_t) = \int (h_{t-1}(x_{t-1}) + \ell_t^{\theta, \phi}(x_{t-1}, x_t)) q_{t-1|t}^{\phi}(x_t, x_{t-1}) \mathrm{d}x_{t-1} \\ = \mathbb{E}_{q_{t-1|t}^{\phi}(x_{t,\cdot})} [h_{t-1}(X_{t-1}) + \ell_t^{\theta, \phi}(X_{t-1}, x_t)].$$

This establishes the recursive expression of the ELBO.

Now, we consider the gradient of the ELBO with respect to $\phi$. Write

$$\begin{aligned} \nabla_{\phi} \mathcal{L}_{t}^{\theta, \phi} &= \nabla_{\phi} \mathbb{E}_{q_{t}^{\phi}} [h_t(X_t)] - \nabla_{\phi} \mathbb{E}_{q_{t}^{\phi}} [\log q_{t}^{\phi}(X_t)] \\ &= \nabla_{\phi} \int (h_t(x_t) - \log q_{t}^{\phi}(x_t)) q_{t}^{\phi}(x_t) \mathrm{d}x_t \\ &= \int (\nabla_{\phi} h_t(x_t) - \nabla_{\phi} \log q_{t}^{\phi}(x_t)) q_{t}^{\phi}(x_t) \mathrm{d}x_t + \int (h_t(x_t) - \log q_{t}^{\phi}(x_t)) \nabla_{\phi} q_{t}^{\phi}(x_t) \mathrm{d}x_t \\ &= \mathbb{E}_{q_{t}^{\phi}} [\nabla_{\phi} h_t(X_t) + (h_t(X_t) - \log q_{t}^{\phi}(X_t)) \nabla_{\phi} \log q_{t}^{\phi}(X_t)], \end{aligned}$$

where we used that $\mathbb{E}_{q_{t}^{\phi}}[\nabla \log q_{t}^{\phi}(X_t)] = 0$. Then, writing $u_t(x_t) = \nabla_{\phi} h_t(x_t)$,

$$\begin{aligned} u_t(x_t) &= \nabla_{\phi} \mathbb{E}_{q_{0:t-1|t}^{\phi}(x_{t,\cdot})} [\ell_{0:t}^{\theta,\phi}(X_{0:t-1}, x_t)] \\ &= \mathbb{E}_{q_{0:t-1|t}^{\phi}(x_{t,\cdot})} \left[ \left( \nabla_{\phi} \log q_{0:t-1|t}^{\phi} \times \ell_{0:t}^{\theta,\phi} \right) (X_{0:t-1}, x_t) \right] + \mathbb{E}_{q_{0:t-1|t}^{\phi}(x_{t,\cdot})} \left[ \nabla_{\phi} \ell_{0:t}^{\theta,\phi} (X_{0:t-1}, x_t) \right]. \end{aligned}$$

Remembering that

$$\ell_{0:t}^{\theta,\phi} = \log(\chi(x_0)g_0^{\theta}(x_0, y_0)) + \sum_{s=1}^{t} \log(m_s^{\theta}(x_{s-1}, x_s)g_s^{\theta}(x_s, y_s)) - \sum_{s=1}^{t} \log q_{s-1|s}^{\phi}(x_s, x_{s-1}),$$

we obtain that $\nabla_{\phi} \ell_{0:t}^{\theta,\phi}(X_{0:t-1}, x_t) = -\nabla_{\phi} \log q_{0:(t-1)|t}^{\phi}(X_{0:t-1}, x_t)$, which has zero expectation under $q_{0:t-1|t}^{\phi}(x_t, \cdot)$. Thus,

$$u_t(x_t) = \mathbb{E}_{q_{0:t-1|t}^{\phi}(x_{t,\cdot})} \left[ \left( \nabla_{\phi} \log q_{0:(t-1)|t}^{\phi} \times \ell_{0:t}^{\theta,\phi} \right) (X_{0:t-1}, x_t) \right].$$

To establish the recursion for $u_t(x_t)$, write

$$u_t(x_t) = \mathbb{E}_{q_{0:t-1|t}^{\phi}(x_{t,\cdot})} \left[ \left( \nabla_{\phi} \log q_{0:(t-2)|t-1}^{\phi}(X_{0:t-1}) + \nabla_{\phi} \log q_{t-1|t}^{\phi}(X_{t-1}, x_t) \right) \times \left( \ell_{0:t-1}^{\theta,\phi}(X_{0:t-1}) + \ell_t^{\theta,\phi}(X_{t-1}, x_t) \right) \right] \quad (27)$$

$$= \mathbb{E}_{q_{t-1|t}^{\phi}(x_{t,\cdot})} [u_{t-1}(X_{t-1})] \quad (28)$$

$$+ \mathbb{E}_{q_{t-1|t}^{\phi}(x_{t,\cdot})} \left[ \nabla_{\phi} \log q_{t-1|t}^{\phi}(X_{t-1}, x_t) \left( \mathbb{E}_{q_{0:(t-2)|t-1}^{\phi}} \left[ \ell_{0:t-1}^{\theta,\phi}(X_{0:t-1}) \right] + \ell_t^{\theta,\phi}(X_{t-1}, x_t) \right) \right] \\ + \mathbb{E}_{q_{t-1|t}^{\phi}(x_{t,\cdot})} \left[ \ell_t^{\theta,\phi}(X_{t-1}, x_t) \times \mathbb{E}_{q_{0:(t-2)|t-1}^{\phi}} \left[ \nabla_{\phi} \log q_{0:(t-2)|t-1}^{\phi}(X_{0:t-1}) \right] \right],$$

which yields

$$u_t(x_t) = \mathbb{E}_{q_{t-1|t}^{\phi}} \left[ u_{t-1}(X_{t-1}) + \nabla_{\phi} \log q_{t-1|t}^{\phi}(X_{t-1}, x_t) \times (h_{t-1}(X_{t-1}) + \ell_t^{\theta,\phi}(X_{t-1}, x_t)) \right],$$

which was to be established.
Finally, let us consider the gradient w.r.t. $\theta$. Using (26), we have that $\nabla_{\theta}\mathcal{L}_{t}^{\theta,\phi} = \mathbb{E}_{q_{t}^{\phi}}[\nabla_{\theta}h_{t}(X_{t})]$. Writing $v_{t} = \nabla_{\theta}h_{t}$, we obtain

$$ 
\begin{aligned}
v_t(x_t) &= \nabla_{\theta} \mathbb{E}_{q_{0:t-1|t}^{\phi}(x_{t,\cdot})} \left[ \ell_{0:t}^{\theta,\phi}(X_{0:t-1}, x_t) \right] \\
&= \mathbb{E}_{q_{0:t-1|t}^{\phi}(x_{t,\cdot})} \left[ \nabla_{\theta} \ell_{0:t-1}^{\theta,\phi}(X_{0:t-2}, X_{t-1}) + \nabla_{\theta} \ell_{t-1}^{\theta,\phi}(X_{t-1}, x_t) \right] \\
&= \mathbb{E}_{q_{0:t-1|t}^{\phi}(x_{t,\cdot})} \left[ \nabla_{\theta} \ell_{0:t-1}^{\theta,\phi}(X_{0:t-2}, X_{t-1}) \right] + \mathbb{E}_{q_{t-1|t}^{\phi}(x_{t,\cdot})} \left[ \nabla_{\theta} \ell_{t-1}^{\theta,\phi}(X_{t-1}, x_t) \right] \\
&= \mathbb{E}_{q_{t-1|t}^{\phi}(x_{t,\cdot})} [v_{t-1}(X_{t-1})] + \mathbb{E}_{q_{t-1|t}^{\phi}(x_{t,\cdot})} \left[ \nabla_{\theta} \ell_{t-1}^{\theta,\phi}(X_{t-1}, x_t) \right],
\end{aligned}
$$

which concludes the proof.

# C Proof of Theorem 4.5

## C.1 Geometric ergodicity of the extended chain

First, some notation. Let $(\mathcal{E}, \mathcal{E})$ be an arbitrary state space. Then a kernel $K$ on $\mathcal{E} \times \mathcal{E}$ induces two endomorphisms, the first acting on the space $\mathbb{F}_b(\mathcal{E})$ of bounded measurable functions on $(\mathcal{E}, \mathcal{E})$ according to

$$ \mathbb{F}_b(\mathcal{E}) \ni f \mapsto Kf(\cdot) := \int f(x) K(\cdot, dx) \in \mathbb{F}_b(\mathcal{E}) $$

and the second acting on the space $\mathcal{M}(\mathcal{E})$ of measures on $(\mathcal{E}, \mathcal{E})$ (we let $\mathcal{M}_1(\mathcal{E}) \subset \mathcal{M}(\mathcal{E})$ denote the subspace of probability measures) according to

$$ \mathcal{M}(\mathcal{E}) \ni \mu \mapsto \mu K(x, \cdot) := \int \mu(dx) K(x, \cdot) \in \mathcal{M}(\mathcal{E}). $$

In addition, the product of two kernels $K$ and $L$ on $(\mathcal{E}, \mathcal{E})$ is defined as the kernel

$$ KL: \mathcal{E} \times \mathcal{E} \ni (x, A) \mapsto \int K(x, dx') L(x', A) $$

on $(\mathcal{E}, \mathcal{E})$. Using this notation, we may define, for any $t \in \mathbb{N}_{>0}$, the power $K^t$ of a kernel $K$ by multiplying $K$ by itself $t-1$ times.

In the following we denote, for every $t \in \mathbb{N}$, by

$$ Q_{t|t+1}^{\phi} f(x_{t+1}) := \int f(x_t) q_{t|t+1}^{\phi}(x_t, x_{t+1}) dx_t, \quad (x_{t+1}, f) \times \mathbb{F}_b(\mathcal{X}), $$

the Markov kernel induced by the transition density $q_{t|t+1}^{\phi}$.

In Section 4 it is assumed that the data generating process $(X_t, Y_t)_{t \in \mathbb{N}}$ is an SSM, and we denote by $R$ its Markov transition kernel of this process. Since the auxiliary states $(a_t)_{t \in \mathbb{N}}$ are generated deterministically from the observations via the mapping $\mathcal{A}^{\phi}$, also the augmented process $(X_t, Y_t, a_t)_{t \in \mathbb{N}}$ is Markov with transition kernel

$$ S^{\phi} f(x_t, y_t, a_t) := \int f(x_{t+1}, y_{t+1}, \mathcal{A}^{\phi}(a_t, y_{t+1})) R(x_t, y_t, d(x_{t+1}, y_{t+1})), \\ (x_t, y_t, a_t, f) \in \mathcal{X} \times \mathcal{Y} \times \mathcal{A} \times \mathbb{F}_b(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A}). $$

For clarity, we now briefly recall the principal assumptions of this work.

**Assumption C.1** (Uniform ergodicity of $(X_t, Y_t, a_t)_{t \in \mathbb{N}}$). There exist $\pi \in \mathcal{M}_1(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ and $\alpha \in (0, 1)$ such that for every $t \in \mathbb{N}_{>0}$, $\phi \in \Phi$, and $(x, y, a) \in \mathcal{X} \times \mathcal{Y} \times \mathcal{A}$,

$$ \| (S^{\phi})^t(x, y, a) - \pi \|_{TV} \le \alpha^t. $$

Theorem C.1 is discussed in Section C.2 below.
**Assumption C.2.** There exist constants $0 < \varepsilon^- < \varepsilon^+ < \infty$ such that for every $\phi \in \Phi$, $t \in \mathbb{N}_{>0}$, and $(x, x') \in \mathcal{X}^2$, $\varepsilon^- \le \psi_t^\phi(x, x') \le \varepsilon^+$.

In addition, we let $\varepsilon := \varepsilon^- / \varepsilon^+$.

**Assumption C.3.** There exists $c \in \mathbb{R}_{>0}$ such that for every $\phi \in \Phi$, $a_t \in \mathbf{A}$, and $(x_{s+1}, y_{s+1}) \in \mathcal{X}^2 \times \mathcal{Y}$,

$$ (i) \quad \int |\nabla_{\phi} \log q_{s|s+1}^{\phi}(x_{s+1}, x_s)|^2 Q_{s|s+1}^{\phi}(x_{s+1}, dx_s) \le c^2, $$

$$ (ii) \quad \int |\ell_{s}^{\theta, \phi}(x_s, x_{s+1})|^2 Q_{s|s+1}^{\phi}(x_{s+1}, dx_s) \le c^2, $$

$$ (iii) \quad \int |\nabla_{\theta} \ell_{s}^{\theta, \phi}(x_s, x_{s+1})| Q_{s|s+1}^{\phi}(x_{s+1}, dx_s) \le c. $$

From now on we let $\mathcal{Z} := \mathbf{X} \times \mathbf{Y} \times \mathbf{A} \times \mathbf{F}_b(\mathcal{X})^3$ denote the state space of the extended chain $(Z_t)_{t \in \mathbb{N}}$ and let $\mathcal{Z} := \mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A} \otimes \mathcal{F}_b(\mathcal{X})^{\otimes 3}$ be the associated $\sigma$-field.

**Definition C.4.** Let $\mathcal{L}(\mathcal{Z})$ be the set of $f \in \mathbf{F}(\mathcal{Z})$ for which there exists $\varphi \in \mathbf{F}_b(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ such that for every $(x, y, a, h, h', u, u', v, v') \in \mathbf{X} \times \mathbf{Y} \times \mathbf{A} \times \mathbf{F}_b(\mathcal{X})^6$,

$$ (i) \quad |f(x, y, a, h, u, v)| \le \varphi(x, y, a), $$

$$ (ii) \quad |f(x, y, a, h, u, v) - f(x, y, a, h', u', v')| \le \varphi(x, y, a) (\operatorname{osc}(h - h') + \operatorname{osc}(u - u') + \operatorname{osc}(v - v')). $$

The following is a slightly more precise statement of Theorem 4.5, our main result.

**Theorem C.5.** Assume C.1, C.3, and C.2. Then there exist $\rho \in (0, 1)$ and a function $\zeta : \mathbf{F}_b(\mathcal{X})^6 \to \mathbb{R}_{>0}$ such that for every $(\theta, \phi) \in \Theta \times \Phi$, $t \in \mathbb{N}$, $f \in \mathcal{L}(\mathcal{Z})$, $z = (x_0, y_0, a_0, h_0, u_0, v_0) \in \mathcal{Z}$, and $z' = (x'_0, y'_0, a'_0, h'_0, u'_0, v'_0) \in \mathcal{Z}$,

$$ |(T^{\theta, \phi})^t f(z) - (T^{\theta, \phi})^t f(z')| \le \zeta(h_0, h'_0, u_0, u'_0, v_0, v'_0) \|\varphi\|_{\infty} \rho^t. \quad (29) $$

Moreover, there exists a kernel $\Pi^{\theta, \phi}$ on $\mathbf{Z} \times \mathcal{L}(\mathcal{Z})$ such that for every $f \in \mathcal{L}(\mathcal{Z})$, $\Pi^{\theta, \phi} f$ is constant and for every $t \in \mathbb{N}$ and $z = (x_0, y_0, a_0, h_0, u_0, v_0) \in \mathcal{Z}$,

$$ |(T^{\theta, \phi})^t f(z) - \Pi^{\theta, \phi} f| \le \bar{\zeta}(h_0, u_0, v_0) \|\varphi\|_{\infty} \rho^t, \quad (30) $$

where

$$ \bar{\zeta}(h_0, u_0, v_0) := \frac{1}{1-\rho} \int \zeta(h_0, h'_0, u_0, u'_0, v_0, v'_0) T^{\theta, \phi}(z, dz'). \quad (31) $$

We preface the proof of Theorem C.5 with a couple of definitions and lemmas. First, we summarize the recursion for the ELBO using the function-valued mapping

$$ \mathcal{H}^{\theta, \phi}(\cdot, a_t, y_{t+1}) : h_t \mapsto \int (h_t(x_t) + \ell_t^{\theta, \phi}(x_t, x_{t+1})) Q_{t|t+1}^{\phi}(x_{t+1}, dx_t), $$

where the kernel $Q_{t|t+1}^{\phi}$ and the term $\ell_t^{\theta, \phi}$ depend implicitly on $a_t$ and $y_{t+1}$, respectively, implying that $h_{t+1}(x_{t+1}) = \mathcal{H}^{\theta, \phi}(h_t, a_t, y_{t+1})(x_{t+1})$. Based on the latter, we also define, for every $t \in \mathbb{N}_{>0}$ and vector $y_{1:t} \in \mathcal{Y}^t$, the composite versions

$$ \mathcal{H}^{\theta, \phi}(\cdot, a_{0:t-1}, y_{1:t}) : h_0 \mapsto \begin{cases} \mathcal{H}^{\theta, \phi}(h_0, a_0, y_1), & \text{for } t=1, \\ \mathcal{H}^{\theta, \phi}(\mathcal{H}_{t-1}^{\theta, \phi}(h_0, a_{0:t-2}, y_{1:t-1}), a_{t-1}, y_t), & \text{for } t \ge 2. \end{cases} \quad (32) $$

Using the similar notation

$$ \mathcal{V}^{\theta, \phi}(\cdot, v_t, a_t) : v_t \mapsto \int (v_t(x_t) + \nabla_{\theta} \ell_t^{\theta, \phi}(x_t, x_{t+1})) Q_{t|t+1}^{\phi}(x_{t+1}, dx_t), $$

for the recursion of the ELBO gradient with respect to $\theta$, the composite mappings $(\mathcal{V}_t^{\theta, \phi})_{t \in \mathbb{N}_{>0}}$ are defined similarly. Finally, letting

$$ \mathcal{U}^{\theta, \phi}(\cdot, a_t, y_{t+1}) : (h_t, u_t) \mapsto \int (u_t(x_t) + \nabla_{\phi} \log q_{t|t+1}^{\phi}(x_t, x_{t+1}) \{h_t(x_t) + \ell_t^{\theta, \phi}(x_t, x_{t+1})\}) Q_{t|t+1}^{\phi}(x_{t+1}, dx_t) $$
summarize the recursion for the ELBO gradient with respect to $\phi$, so that $u(x_{t+1}) = \mathcal{U}^{\theta,\phi}(h_t, u_t, a_t, y_{t+1})(x_{t+1})$, we also define the compositions

$$ \mathcal{U}_t^{\theta,\phi}(\cdot, a_{0:t-1}, y_{1:t}) : (h_0, u_0) \mapsto \begin{cases} \mathcal{U}^{\theta,\phi}(h_0, u_0, a_0, y_1), & \text{for } t=1, \\ \mathcal{U}^{\theta,\phi}(\mathcal{H}_{t-1}^{\theta,\phi}(h_0, a_{0:t-2}, y_{1:t-1}), \mathcal{U}_{t-1}^{\theta,\phi}(h_0, u_0, a_{0:t-2}, y_{1:t-1}), a_{t-1}, y_t), & \text{for } t \ge 2. \end{cases} \tag{33} $$

Using these definitions, the transition kernel $T^{\theta,\phi}$ of the extended chain $(Z_t)_{t \in \mathbb{N}}$ can be expressed as

$$ T^{\theta,\phi} f(z_t) = \int f(x_{t+1}, y_{t+1}, \mathcal{A}^{\theta,\phi}(a_t, y_{t+1}), \mathcal{H}^{\theta,\phi}(h_t, a_t, y_{t+1}), \mathcal{U}^{\theta,\phi}(h_t, u_t, a_t, y_{t+1}), \mathcal{V}^{\theta,\phi}(v_t, a_t, y_{t+1})) \\ \times R((x_t, y_t), d(x_{t+1}, y_{t+1})), \quad (z_t, f) \in \mathbb{Z} \times \mathbf{F}_b(\mathcal{Z}), \tag{34} $$

where $z_t = (x_t, y_t, a_t, h_t, u_t, v_t)$.

The following lemma establishes the geometric contraction of the function-valued mappings defined above.

**Lemma C.6.** Assume C.2. Then for every $t \in \mathbb{N}$, $y_{1:t} \in \mathcal{Y}^t$, $a_{0:t-1} \in \mathcal{A}^t$, and $(h_0, h'_0, v_0, v'_0) \in \mathbf{F}_b(\mathcal{X})^4$,

$$ (i) \quad \mathrm{osc} \left( \mathcal{H}_t^{\theta,\phi}(h_0, a_{0:t-1}, y_{1:t}) - \mathcal{H}_t^{\theta,\phi}(h'_0, a_{0:t-1}, y_{1:t}) \right) \le (1-\varepsilon)^t \mathrm{osc}(h_0 - h'_0), $$

$$ (ii) \quad \mathrm{osc} \left( \mathcal{Y}_t^{\theta,\phi}(v_0, a_{0:t-1}, y_{1:t}) - \mathcal{Y}_t^{\theta,\phi}(v'_0, a_{0:t-1}, y_{1:t}) \right) \le (1-\varepsilon)^t \mathrm{osc}(v_0 - v'_0), $$

where $\varepsilon \in (0, 1)$ is given in Theorem C.2.

(iii) Assume additionally C.3 (i). Then for every $\varrho \in (1-\varepsilon, 1)$ there exists $d > 0$ such that for every $t \in \mathbb{N}$, $y_{1:t} \in \mathcal{Y}^t$, $a_{0:t-1} \in \mathcal{A}^t$, and $(h_0, h'_0, u_0, u'_0) \in \mathbf{F}_b(\mathcal{X})^4$,

$$ \mathrm{osc} \left( \mathcal{U}_t^{\theta,\phi}(h_0, u_0, a_{0:t-1}, y_{1:t}) - \mathcal{U}_t^{\theta,\phi}(h'_0, u'_0, a_{0:t-1}, y_{1:t}) \right) \le (1-\varepsilon)^t \mathrm{osc}(u_0 - u'_0) + d\varrho^t \mathrm{osc}(h_0 - h'_0). $$

The proof of Lemma C.6 is based on the following lemmas.

**Lemma C.7.** for every $t \in \mathbb{N}_{>0}$, $y_{1:t} \in \mathcal{Y}^t$, $a_{0:t-1} \in \mathcal{A}^t$, $(h_0, u_0) \in \mathbf{F}_b(\mathcal{X})^2$, and $x_t \in \mathcal{X}$,

$$ \mathcal{U}_t^{\theta,\phi}(h_0, u_0, a_{0:t-1}, y_{1:t})(x_t) = \int \dots \int \left( u_0(x_0) + \sum_{s=0}^{t-1} \nabla_{\phi} \log q_{s|s+1}^{\phi}(x_{s+1}, x_s) \{ \mathcal{H}_s^{\theta,\phi}(h_0, a_{0:s-1}, y_{1:s}) + \ell_s^{\theta,\phi}(x_s, x_{s+1}) \} \right) \prod_{s=0}^{t-1} Q_{s|s+1}^{\phi}(x_{s+1}, dx_s) $$

with the convention $\mathcal{H}_0^{\theta,\phi}(h_0, a_{0:-1}, y_{1:0}) := h_0$.

**Proof of Lemma C.7.** We proceed by induction and assume that the claim holds true for $t \in \mathbb{N}_{>0}$. By definition (33),

$$ \begin{aligned} \mathcal{U}_{t+1}^{\theta,\phi}(h_0, u_0, a_{0:t}, y_{1:t+1})(x_{t+1}) &= \mathcal{U}^{\theta,\phi}(\mathcal{H}_t^{\theta,\phi}(h_0, a_{0:t-1}, y_{1:t}), \mathcal{U}_t^{\theta,\phi}(u_0, a_{0:t-1}, y_{1:t}), a_t, y_{t+1}) \\ &= \int \mathcal{U}_t^{\theta,\phi}(u_0, a_{0:t-1}, y_{1:t})(x_t) + \nabla_{\phi} \log q_{t|t+1}^{\phi}(x_{t+1}, x_t) \{ \mathcal{H}_t^{\theta,\phi}(h_0, a_{0:t-1}, y_{1:t})(x_t) + \ell_t^{\theta,\phi}(x_t, x_{t+1}) \} Q_{t|t+1}^{\phi}(x_{t+1}, dx_t). \end{aligned} $$

Now, inserting the induction hypothesis into the right-hand side of the previous expression yields

$$ \begin{aligned} \mathcal{U}_{t+1}^{\theta,\phi}(h_0, u_0, a_{0:t}, y_{1:t+1})(x_{t+1}) &= \int \dots \int \left( u_0(x_0) + \sum_{s=0}^{t} \nabla_{\phi} \log q_{s|s+1}^{\phi}(x_{s+1}, x_s) \{ \mathcal{H}_s^{\theta,\phi}(h_0, a_{0:s-1}, y_{1:s}) + \ell_s^{\theta,\phi}(x_s, x_{s+1}) \} \right) \prod_{s=1}^{t} Q_{s|s+1}^{\phi}(x_{s+1}, dx_s), \end{aligned} $$

which establishes the induction step.

Finally, we note that the base case $t=1$ holds true, since by definition (33)

$$ \begin{aligned} \mathcal{U}_1^{\theta,\phi}(h_0, u_0, a_0, y_1)(x_1) &= \mathcal{U}^{\theta,\phi}(h_0, u_0, a_0, y_1)(x_1) \\ &= \int \left( u_0(x_0) + \nabla_{\phi} \log q_{0|1}^{\phi}(x_1, x_0) \{h_0(x_0) + \ell_0^{\theta,\phi}(x_0, x_{0+1})\} \right) Q_{0|1}^{\phi}(x_1, dx_0). \end{aligned} $$

This completes the proof. $\square$
In the following, let $\beta(M)$ denote the Dobrushin coefficient of a Markov kernel $M$.

**Lemma C.8.** Assume C.2. Then for every $t \in \mathbb{N}$, $\beta(Q_{t|t+1}^{\phi}) \le 1 - \varepsilon$.

*Proof of Theorem C.8.* Pick arbitrarily $(x_{t+1}, f) \in \mathcal{X} \times \mathcal{F}_b(\mathcal{X})$ and write, using definition (6) and Theorem C.2,

$$Q_{t|t+1}^{\phi} f(x_{t+1}) = \frac{\int f(x_t) \psi_t^{\phi}(x_t, x_{t+1}) q_t^{\phi}(dx_t)}{\int \psi_t^{\phi}(x'_t, x_{t+1}) q_t^{\phi}(dx'_t)} \ge \frac{\varepsilon^{-}}{\varepsilon^{+}} q_t^{\phi} f = \varepsilon q_t^{\phi} f,$$

which means that $Q_{t|t+1}^{\phi}$ allows $\mathcal{X}$ as a 1-small set with respect to $(q_t^{\phi}, \varepsilon)$. From this it follows that $\beta(Q_{t|t+1}^{\phi}) \le 1 - \varepsilon$. $\square$

We are now ready to establish Lemma C.6.

*Proof of Lemma C.6.* To establish (i), let $(Q_{s-1|s}^{\phi})_{s=1}^t$ denote the backward transition kernels associated with $(q_{s-1|s}^{\phi})_{s=1}^t$. We may then write, for every $x_t \in \mathcal{X}$,

$$\mathcal{H}_t^{\theta, \phi}(h_0, y_{1:t})(x_t) - \mathcal{H}_t^{\theta, \phi}(h'_0, y_{1:t})(x_t) = Q_{t-1|t}^{\phi} \cdots Q_{0|1}^{\phi}(h_0 - h'_0)(x_t).$$

Now, recall that for every $s$ and $h \in \mathcal{F}_b(\mathcal{X})$,

$$\mathrm{osc}(Q_{s-1|s}^{\phi} h) \le \beta(Q_{s-1|s}^{\phi}) \mathrm{osc}(h), \qquad (35)$$

where $\beta(Q_{s-1|s}^{\phi})$ is the Dobrushin coefficient of $Q_{s-1|s}^{\phi}$, and iterating the bound (35) yields

$$\mathrm{osc}(Q_{t-1|t}^{\phi} \cdots Q_{0|1}^{\phi}(h_0 - h'_0)) \le \left( \prod_{s=1}^{t} \beta(Q_{s-1|s}^{\phi}) \right) \mathrm{osc}(h_0 - h'_0). \qquad (36)$$

From this the claim (i) follows by Theorem C.8.

To prove (ii), note that

$$\begin{aligned} \mathcal{V}_t^{\theta, \phi}(v_0, a_{0:t-1}, y_{1:t})(x_t) - \mathcal{V}_t^{\theta, \phi}(v'_0, a_{0:t-1}, y_{1:t})(x_t) &= Q_{t-1|t}^{\phi} \cdots Q_{0|1}^{\phi}(v_0 - v'_0)(x_t) \\ &= \mathcal{H}_t^{\theta, \phi}(v_0, a_{0:t-1}, y_{1:t})(x_t) - \mathcal{H}_t^{\theta, \phi}(v'_0, a_{0:t-1}, y_{1:t})(x_t). \end{aligned}$$

Thus, (ii) follows immediately from (i).

Finally, to establish (iii), write, using Lemma C.7,

$$\begin{aligned} & \mathcal{U}_t^{\theta, \phi}(h_0, u_0, a_{0:t-1}, y_{1:t})(x_t) - \mathcal{U}_t^{\theta, \phi}(h'_0, u'_0, a_{0:t-1}, y_{1:t})(x_t) \\ &= \int \cdots \int \left( u_0(x_0) - u'_0(x_0) + \sum_{s=0}^{t-1} \nabla_{\phi} \log q_{s|s+1}^{\phi}(x_{s+1}, x_s) \{ \mathcal{H}_s^{\theta, \phi}(h_0, a_{0:s-1}, y_{1:s})(x_s) - \mathcal{H}_s^{\theta, \phi}(h'_0, a_{0:s-1}, y_{1:s})(x_s) \} \right) \\ & \qquad \times \prod_{s=0}^{t-1} Q_{s|s+1}^{\phi}(x_{s+1}, dx_s) \\ &= Q_{t-1|t}^{\phi} \cdots Q_{0|1}^{\phi}(u_0 - u'_0)(x_t) + \sum_{s=0}^{t-1} Q_{t-1|t}^{\phi} \cdots Q_{s+1|s+2}^{\phi} \varphi_s(x_t), \end{aligned} \qquad (37)$$

where we have set

$$\varphi_s(x_{s+1}) := \int \nabla_{\phi} \log q_{s|s+1}^{\phi}(x_{s+1}, x_s) \{ \mathcal{H}_s^{\theta, \phi}(h_0, a_{0:s-1}, y_{1:s})(x_s) - \mathcal{H}_s^{\theta, \phi}(h'_0, a_{0:s-1}, y_{1:s})(x_s) \} Q_{s|s+1}^{\phi}(x_{s+1}, dx_s).$$

Now, note that since

$$\int \nabla_{\phi} \log q_{s|s+1}^{\phi}(x_{s+1}, x_s) Q_{s|s+1}^{\phi}(x_{s+1}, dx_s) = 0,$$

it holds, for every $c \in \mathbb{R}$,

$$\|\varphi_s\|_{\infty} \le \|\mathcal{H}_s^{\theta, \phi}(h_0, a_{0:s-1}, y_{1:s}) - \mathcal{H}_s^{\theta, \phi}(h'_0, a_{0:s-1}, y_{1:s}) - c\|_{\infty} \int |\nabla_{\phi} \log q_{s|s+1}^{\phi}(x_{s+1}, x_s)| Q_{s|s+1}^{\phi}(x_{s+1}, dx_s).$$
Thus, using Theorem C.3(i) and the fact that for all $f \in \mathcal{F}(\mathcal{X})$, $\mathrm{osc}(f) = 2 \inf_{c \in \mathbb{R}} \|f - c\|_{\infty}$, it holds, by (i), that

$$ \|\varphi_s\|_{\infty} \le \frac{1}{2} c \mathrm{osc} \left( \mathcal{H}_s^{\theta, \phi}(h_0, a_{0:s-1}, y_{1:s}) - \mathcal{H}_s^{\theta, \phi}(h'_0, a_{0:s-1}, y_{1:s}) \right) \le \frac{1}{2} c (1 - \varepsilon)^s \mathrm{osc}(h_0 - h'_0). $$

As a consequence, by Theorem C.8,

$$ \begin{aligned} \mathrm{osc}(Q_{t-1|t}^{\phi} \cdots Q_{s+1|s+2}^{\phi} \varphi_s) &\le \left( \prod_{\ell=s+1}^{t-1} \beta(Q_{\ell|\ell+1}^{\phi}) \right) \mathrm{osc}(\varphi_s) \\ &\le c(1-\varepsilon)^{t-s-1} (1-\varepsilon)^s \mathrm{osc}(h_0 - h'_0) \\ &= c(1-\varepsilon)^{t-1} \mathrm{osc}(h_0 - h'_0). \end{aligned} \quad (38) $$

Moreover, since

$$ Q_{t-1|t}^{\phi} \cdots Q_{0|1}^{\phi} (u_0 - u'_0)(x_t) = \mathcal{H}_t^{\theta, \phi}(u_0, a_{0:t-1}, y_{1:t})(x_t) - \mathcal{H}_t^{\theta, \phi}(u'_0, a_{0:t-1}, y_{1:t})(x_t), $$

(i) implies that

$$ \mathrm{osc}(Q_{t-1|t}^{\phi} \cdots Q_{0|1}^{\phi} (u_0 - u'_0)) \le (1 - \varepsilon)^t \mathrm{osc}(u_0 - u'_0). \quad (39) $$

Combining (37), (38), and (39) yields

$$ \mathrm{osc} \left( \mathcal{U}_t^{\theta, \phi}(h_0, u_0, a_{0:t-1}, y_{1:t}) - \mathcal{U}_t^{\theta, \phi}(h'_0, u'_0, a_{0:t-1}, y_{1:t}) \right) \le (1 - \varepsilon)^t \mathrm{osc}(u_0 - u'_0) + ct(1 - \varepsilon)^{t-1} \mathrm{osc}(h_0 - h'_0). $$

Finally, the claim (iii) follows by picking $\varrho \in (1 - \varepsilon, 1)$ and letting $d := ((1 - \varepsilon) \log\{\varrho/(1 - \varepsilon)\})^{-1}$. $\square$

We are now ready to establish Theorem C.5, following the same lines as Theorem C.13.

*Proof of Theorem C.5.* First, denote

$$ f_{s,t}(x_t, y_{s:t}, a_{s-1:t}, h, u, v) := f(x_t, y_t, a_t, \mathcal{H}_{t-s+1}^{\theta, \phi}(h, a_{s-1:t-1}, y_{s:t}), \mathcal{U}_{t-s+1}^{\theta, \phi}(h, u, a_{s-1:t-1}, y_{s:t}), \mathcal{V}_{t-s+1}^{\theta, \phi}(v, a_{s-1:t-1}, y_{s:t})). $$

for $s \in [1, t]$, where we have omitted the dependence on $\theta$ and $\phi$ for brevity. Note that with this notation, for $z = (x_0, y_0, a_0, h_0, u_0, v_0)$,

$$ (T^{\theta, \phi})^t f(z) = \int \cdots \int f_{1,t}(x_t, y_{1:t}, a_{0:t}, h_0, u_0, v_0) \prod_{s=0}^{t-1} K^{\theta, \phi}((x_s, y_s, a_s), d(x_{s+1}, y_{s+1}, a_{s+1})). $$

Now, picking $(\tilde{h}_0, \tilde{u}_0, \tilde{v}_0) \in \mathcal{F}_b(\mathcal{X})^3$ arbitrarily and using the decomposition

$$ \begin{aligned} f_{1,t}(x_t, y_{1:t}, a_{0:t}, h_0, u_0, v_0) = f_{1,t}(x_t, y_{1:t}, a_{0:t}, h_0, u_0, v_0) &- f_{1,t}(x_t, y_{1:t}, a_{0:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) \\ &+ \sum_{s=1}^{t-1} \left( f_{s,t}(x_t, y_{s:t}, a_{s-1:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) - f_{s+1,t}(x_t, y_{s+1:t}, a_{s:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) \right) + f_{t,t}(x_t, y_t, a_{t-1:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0), \end{aligned} $$
which is adopted from [Tadić and Doucet, 2005], we may write,

$$
\begin{align}
& (T^{\theta, \phi})^t f(z) - (T^{\theta, \phi})^t f(z') \nonumber \\
&= \int \cdots \int \left( f_{1,t}(x_t, y_{1:t}, a_{0:t}, h_0, u_0, v_0) - f_{1,t}(x_t, y_{1:t}, a_{0:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) \right) \prod_{s=0}^{t-1} S^{\phi}((x_s, y_s, a_s), d(x_{s+1}, y_{s+1}, a_{s+1})) \nonumber \\
&\quad - \int \cdots \int \left( f_{1,t}(x'_t, y'_{1:t}, a'_{0:t}, h'_0, u'_0, v'_0) - f_{1,t}(x'_t, y'_{1:t}, a'_{0:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) \right) \prod_{s=0}^{t-1} S^{\phi}((x'_s, y'_s, a'_s), d(x'_{s+1}, y'_{s+1}, a'_{s+1})) \nonumber \\
&\quad + \sum_{s=1}^t \int \cdots \int \left( f_{s,t}(x_t, y_{s:t}, a_{s-1:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) - f_{s+1,t}(x_t, y_{s+1:t}, a_{s:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) \right) \nonumber \\
&\qquad \times \Delta S_s^{\phi}((x_0, y_0, a_0), d(x_s, y_s, a_s)) \prod_{\ell=s}^{t-1} S^{\phi}((x_{\ell}, y_{\ell}, a_{\ell}), d(x_{\ell+1}, y_{\ell+1}, a_{\ell+1})) \nonumber \\
&\quad - \sum_{s=1}^t \int \cdots \int \left( f_{s,t}(x_t, y_{s:t}, a_{s-1:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) - f_{s+1,t}(x_t, y_{s+1:t}, a_{s:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) \right) \nonumber \\
&\qquad \times \Delta S_s^{\phi}((x'_0, y'_0, a'_0), d(x_s, y_s, a_s)) \prod_{\ell=s}^{t-1} S^{\phi}((x_{\ell}, y_{\ell}, a_{\ell}), d(x_{\ell+1}, y_{\ell+1}, a_{\ell+1})) \nonumber \\
&\quad + \int f_{t,t}(x_t, y_t, a_{t-1:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) \Delta S_t^{\phi}((x_0, y_0, a_0), d(x_t, y_t, a_t)) \nonumber \\
&\quad - \int f_{t,t}(x_t, y_t, a_{t-1:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0) \Delta S_t^{\phi}((x'_0, y'_0, a'_0), d(x_t, y_t, a_t)), \tag{40}
\end{align}
$$

where we have defined, for $s \in [1, t]$, the signed kernels

$$
\Delta S_s^{\phi} f(x, y, a) := (S^{\phi})^s f(x, y, a) - \pi f, \quad (f, x, y, a) \in \mathbb{F}_b(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A}) \times \mathcal{X} \times \mathcal{Y} \times \mathcal{A}.
$$

Note that by Theorem C.4, for every $s \in [1, t]$,

$$
\begin{align}
& |f_{s,t}(x_t, y_{s:t}, a_{s-1:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0)| \nonumber \\
&\le \varphi(x_t, y_t, a_t) \left( \begin{aligned}[t]
& \operatorname{osc}(\mathcal{H}_{t-s+1}^{\theta,\phi}(\tilde{h}_0, a_{s-1:t-1}, y_{s:t}) - \mathcal{H}_{t-s}^{\theta,\phi}(\tilde{h}_0, a_{s:t-1}, y_{s+1:t})) \\
& + \operatorname{osc}(\mathcal{U}_{t-s+1}^{\theta,\phi}(\tilde{h}_0, \tilde{u}_0, a_{s-1:t-1}, y_{s:t}) - \mathcal{U}_{t-s}^{\theta,\phi}(\tilde{h}_0, \tilde{u}_0, a_{s:t-1}, y_{s+1:t})) \\
& + \operatorname{osc}(\mathcal{V}_{t-s+1}^{\theta,\phi}(\tilde{v}_0, a_{s-1:t-1}, y_{s:t}) - \mathcal{V}_{t-s}^{\theta,\phi}(\tilde{v}_0, a_{s:t-1}, y_{s+1:t}))
\end{aligned} \right). \tag{41}
\end{align}
$$

Using definition (32) and Lemma C.6(i), we conclude that

$$
\begin{align*}
& \operatorname{osc} \left( \mathcal{H}_{t-s+1}^{\theta,\phi}(\tilde{h}_0, a_{s-1:t-1}, y_{s:t}) - \mathcal{H}_{t-s}^{\theta,\phi}(\tilde{h}_0, a_{s:t-1}, y_{s+1:t}) \right) \\
&= \operatorname{osc} \left( \mathcal{H}_{t-s}^{\theta,\phi}(\mathcal{H}^{\theta,\phi}(\tilde{h}_0, y_s), a_{s:t-1}, y_{s+1:t}) - \mathcal{H}_{t-s}^{\theta,\phi}(\tilde{h}_0, a_{s:t-1}, y_{s+1:t}) \right) \\
&\le (1 - \varepsilon)^{t-s} \operatorname{osc} \left( \mathcal{H}^{\theta,\phi}(\tilde{h}_0, y_s) - \tilde{h}_0 \right).
\end{align*}
$$

Since

$$
\begin{align*}
\| \mathcal{H}^{\theta, \phi}(\tilde{h}_0, y_s) - \tilde{h}_0 \|_{\infty} &= \left\| \int (\tilde{h}_0(x) + \ell_{s-1}^{\theta, \phi}(x, \cdot)) Q_{s-1|s}^{\phi}(\cdot, dx) - \tilde{h}_0(\cdot) \right\|_{\infty} \\
&\le 2 \| \tilde{h}_0 \|_{\infty} + c,
\end{align*}
$$

where the constant $c \in \mathbb{R}_{>0}$ is provided by Theorem C.3(ii), it holds that

$$
\operatorname{osc} \left( \mathcal{H}_{t-s+1}^{\theta, \phi} (\tilde{h}_0, a_{s-1:t-1}, y_{s:t}) - \mathcal{H}_{t-s}^{\theta, \phi} (\tilde{h}_0, a_{s:t-1}, y_{s+1:t}) \right) \le 2(2\|\tilde{h}_0\|_{\infty} + c)(1-\varepsilon)^{t-s}. \quad (42)
$$

By the same arguments it is shown that

$$
\operatorname{osc} \left( \mathcal{V}_{t-s+1}^{\theta, \phi} (\tilde{v}_0, a_{s-1:t-1}, y_{s:t}) - \mathcal{V}_{t-s}^{\theta, \phi} (\tilde{v}_0, a_{s:t-1}, y_{s+1:t}) \right) \le 2(2\|\tilde{v}_0\|_{\infty} + c)(1-\varepsilon)^{t-s}. \quad (43)
$$
In addition, similarly, using Lemma C.6(iii),

$$
\begin{align*}
\operatorname{osc} (\mathcal{U}_{t-s+1}^{\theta,\phi}(\tilde{h}_0, \tilde{u}_0, a_{s-1:t-1}, y_{s:t}) - \mathcal{U}_{t-s}^{\theta,\phi}(\tilde{h}_0, \tilde{u}_0, a_{s:t-1}, y_{s+1:t}))
\end{align*}
$$
$$
= \operatorname{osc} (\mathcal{U}_{t-s}^{\theta,\phi}(\mathcal{H}^{\theta,\phi}(\tilde{h}_0, y_s), \mathcal{U}^{\theta,\phi}(\tilde{h}_0, \tilde{u}_0, y_s), a_{s:t-1}, y_{s+1:t}) - \mathcal{U}_{t-s}^{\theta,\phi}(\tilde{h}_0, \tilde{u}_0, a_{s:t-1}, y_{s+1:t}))
$$
$$
\leq (1-\varepsilon)^{t-s} \operatorname{osc} (\mathcal{H}^{\theta,\phi}(\tilde{h}_0, y_s) - \tilde{h}_0) + d\varrho^{t-s} \operatorname{osc} (\mathcal{U}^{\theta,\phi}(\tilde{h}_0, \tilde{u}_0, y_s) - \tilde{u}_0),
$$

and since by Theorem C.3(i–ii) and the Cauchy–Schwarz inequality,

$$
\begin{align*}
\| \mathcal{U}^{\theta, \phi}(\tilde{h}_0, \tilde{u}_0, y_s) - \tilde{u}_0 \|_{\infty} &= \left\| \int (\tilde{u}_0(x) + \nabla_{\phi} \log q_{s-1|s}^{\phi}(\cdot, x) \{\tilde{h}_0(x) + \ell_{s-1}^{\theta,\phi}(x, \cdot)\}) Q_{s-1|s}^{\phi}(\cdot, dx) - \tilde{h}_0(\cdot) \right\|_{\infty} \\
&\le 2 \|\tilde{u}_0\|_{\infty} + c \|\tilde{h}_0\|_{\infty} + c^2,
\end{align*}
$$

we may conclude that

$$
\begin{align}
& \operatorname{osc} \left( \mathcal{U}_{t-s+1}^{\theta,\phi}(\tilde{h}_0, \tilde{u}_0, a_{s-1:t-1}, y_{s:t}) - \mathcal{U}_{t-s}^{\theta,\phi}(\tilde{h}_0, \tilde{u}_0, a_{s:t-1}, y_{s+1:t}) \right) \nonumber \\
& \qquad \le 2(2\|\tilde{h}_0\|_{\infty} + c)(1-\varepsilon)^{t-s} + 2d(2\|\tilde{u}_0\|_{\infty} + c\|\tilde{h}_0\|_{\infty} + c^2)\varrho^{t-s}. \tag{44}
\end{align}
$$

Combining (41), (42), (43), and (44) yields the bound

$$
\begin{align*}
|f_{s,t}(x_t, y_{s:t}, a_{s-1:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0)| - f_{s+1,t}(x_t, y_{s+1:t}, a_{s:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0)|
\end{align*}
$$
$$
\begin{align*}
&\le 2\varphi(x_t, y_t, a_t) \left( 3(2\|\tilde{h}_0\|_{\infty} + c)(1-\varepsilon)^{t-s} + d(2\|\tilde{u}_0\|_{\infty} + c\|\tilde{h}_0\|_{\infty} + c^2)\varrho^{t-s} \right) \\
&\le 2\varphi(x_t, y_t, a_t) \left( (6+cd)\|\tilde{h}_0\|_{\infty} + 2d\|\tilde{u}_0\|_{\infty} + 3c + c^2d \right) \varrho^{t-2}.
\end{align*}
\tag{45}
$$

Similarly,

$$
|f_{1,t}(x_t, y_{1:t}, a_{0:t}, h_0, u_0, v_0) - f_{1,t}(x_t, y_{1:t}, a_{0:t}, \tilde{h}_0, \tilde{u}_0, \tilde{v}_0)| \le \varphi(x_t, y_t, a_t) \zeta''(h_0, u_0, v_0) \varrho^t, \quad (46)
$$

where we have defined the mapping

$$
\zeta''(h_0, u_0, v_0) := 2 \left( (d+1)(\|h_0\| + \|\tilde{h}_0\|_{\infty}) + \|u_0\| + \|\tilde{u}_0\|_{\infty} + \|v_0\| + \|\tilde{v}_0\|_{\infty} \right).
$$

Now, letting

$$
\zeta'(h_0, h'_0, u_0, u'_0, v_0, v'_0) := 2 \left( (6+cd)\|\tilde{h}_0\|_{\infty} + 2d\|\tilde{u}_0\|_{\infty} + 3c + c^2d \right) \lor \zeta''(h_0, u_0, v_0) \lor \zeta''(h'_0, u'_0, v'_0).
$$

Applying the bounds (45) and (46) to the decomposition (40) yields

$$
\begin{align*}
|(T^{\theta, \phi})^t f(z) - (T^{\theta, \phi})^t f(z')| &\le \zeta'(h_0, h'_0, u_0, u'_0, v_0, v'_0) \varrho^t \left( (S^\phi)^t \varphi(x_0, y_0, a_0) + (S^\phi)^t \varphi(x'_0, y'_0, a'_0) \right) \\
&\quad + \zeta'(h_0, h'_0, u_0, u'_0, v_0, v'_0) \sum_{s=1}^t \varrho^{t-s} \left( \Delta S_s^\phi (S^\phi)^{t-s} \varphi(x_0, y_0, a_0) + \Delta S_s^\phi (S^\phi)^{t-s} \varphi(x'_0, y'_0, a'_0) \right) \\
&\quad + \Delta S_t^\phi \varphi(x_0, y_0, a_0) + \Delta S_t^\phi \varphi(x'_0, y'_0, a'_0).
\end{align*}
$$

Now, by applying Theorem C.1 to the right-hand side of (3) we obtain

$$
|(T^{\theta, \phi})^t f(z) - (T^{\theta, \phi})^t f(z')| \le \|\varphi\|_{\infty} (\zeta'(h_0, h'_0, u_0, u'_0, v_0, v'_0) \alpha^t + \zeta'(h_0, h'_0, u_0, u'_0, v_0, v'_0) t(\alpha \lor \varrho)^t + \varrho^t),
$$

from which (29) follows by picking $\rho \in (\alpha \lor \varrho, 1)$ and

$$
\zeta(h_0, h'_0, u_0, u'_0, v_0, v'_0) := 2(\zeta'(h_0, h'_0, u_0, u'_0, v_0, v'_0) + 1)(\alpha \lor \varrho)/\rho + 2\zeta'(h_0, h'_0, u_0, u'_0, v_0, v'_0) (e^{|\log\{(\alpha \lor \varrho)/\rho\}|} - 1).
$$

To prove the second claim, first note that by (29),

$$
\begin{align}
|(T^{\theta, \phi})^{t+1} f(z) - (T^{\theta, \phi})^t f(z)| &\le \int |(T^{\theta, \phi})^t f(z') - (T^{\theta, \phi})^t f(z)| T^{\theta, \phi}(z, dz') \nonumber \\
&\le \|\varphi\|_{\infty} \rho^t \int \zeta(h_0, h'_0, u_0, u'_0, v_0, v'_0) T^{\theta, \phi}(z, dz'). \tag{47}
\end{align}
$$
Now, define the kernel

$$ \Pi^{\theta, \phi} f(z) := f(z) + \sum_{t=0}^{\infty} ((T^{\theta, \phi})^{t+1} f(z) - (T^{\theta, \phi})^t f(z)) $$

on $Z \times \mathcal{L}(\mathcal{Z})$. With this definition, note that by (47),

$$ \begin{aligned} |(T^{\theta, \phi})^t f(z) - \Pi^{\theta, \phi} f(z)| & \le \sum_{s=t}^{\infty} |(T^{\theta, \phi})^{s+1} f(z) - (T^{\theta, \phi})^s f(z)| \\ & \le \bar{\zeta}(h_0, u_0, v_0) \|\varphi\|_{\infty} \rho^t, \end{aligned} $$

where $\bar{\zeta}(h_0, u_0, v_0)$ is defined in (31), which establishes (30). Finally, it remains to prove that the function $\Pi^{\theta, \phi} f$ is constant. For this purpose, pick arbitrarily $(z, z') \in \mathbb{Z}^2$; then, however, by (29) and (30),

$$ |\Pi^{\theta, \phi} f(z) - \Pi^{\theta, \phi} f(z')| \le \inf_{t \in \mathbb{N}} (|(T^{\theta, \phi})^{t+1} f(z) - (T^{\theta, \phi})^t f(z)| + |(T^{\theta, \phi})^t f(z) - \Pi^{\theta, \phi} f(z)| \\ + |(T^{\theta, \phi})^t f(z') - \Pi^{\theta, \phi} f(z')|) = 0, $$

from which the claim follows.

<figure><img src="image_5.png" /></figure>

## C.2 Discussion on Theorem C.1

Recall Theorem C.1, under which it is supposed that for every $\phi \in \Phi$ there exist $\alpha \in (0, 1)$ and $\pi^{\phi} \in \mathbb{M}_1(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ ($\alpha$ being independent of $\phi$) such that for every $t \in \mathbb{N}$, $(x, y, a) \in \mathcal{X} \times \mathcal{Y} \times \mathcal{A}$, and $f \in \mathbb{F}_b(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$,

$$ |(S^{\phi})^t f(x, y, a) - \pi^{\phi} f| \le \|f\|_{\infty} \alpha^t, \quad (48) $$

where

$$ S^{\phi} f(x_t, y_t, a_t) := \int f(x_{t+1}, y_{t+1}, \mathcal{A}^{\phi}(a_t, y_{t+1})) S((x_t, y_t), d(x_{t+1}, y_{t+1})), \\ (x_t, y_t, a_t, f) \in \mathcal{X} \times \mathcal{Y} \times \mathcal{A} \times \mathbb{F}_b(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A}), $$

is the transition kernel of the Markov chain $(X_t, Y_t, a_t)_{t \in \mathbb{N}}$.

In this section, we will present conditions under which (48) applies, if not for all $f \in \mathbb{F}_b(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$, then at least for all $f$ in a certain Lipschitz subclass $\mathcal{L}(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ of $\mathbb{F}(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ to be specified. Recall the mapping $\mathcal{A}^{\phi}$ introduced in Section 3 and define, for every $t \in \mathbb{N}_{>0}$ and vectors $y_{1:t} \in \mathcal{Y}^t$ and $a_{0:t-1} \in \mathcal{A}^t$, the composite versions

$$ \mathcal{A}_t^{\phi}(\cdot, y_{1:t}) := \begin{cases} \mathcal{A}^{\phi}(\cdot, y_1), & \text{for } t=1, \\ \mathcal{A}^{\phi}(\mathcal{A}_{t-1}^{\phi}(\cdot, y_{1:t-1}), y_t), & \text{for } t \ge 2. \end{cases} \quad (49) $$

We will assume that these satisfy the following assumption.

**Assumption C.9.** There exist $c > 0$ and $\kappa \in (0, 1)$ such that for every $\phi \in \Phi$, $t \in \mathbb{N}$, $y_{1:t} \in \mathcal{Y}^t$, and $(a, \tilde{a}) \in \mathcal{A}^2$,

$$ (i) \quad \| \mathcal{A}_t^{\phi}(a, y_{1:t}) - \mathcal{A}_t^{\phi}(\tilde{a}, y_{1:t}) \|_2 \le c \|a - \tilde{a}\|_2 \kappa^t, $$

$$ (ii) \quad \| \mathcal{A}^{\phi}(a, y_t) - a \|_2 \le c. $$

**Remark C.10.** To illustrate Theorem C.9(i), let $\mathcal{A}^{\phi}$ correspond to the following vanilla recurrent neural network architecture, where $a_0 \in \mathbb{R}^d$ is the initial hidden state. For $t \in \mathbb{N}_{>0}$, given a sequence $y_{1:t}$ of observations, the hidden state $a_t$ is given by

$$ a_t = \mathcal{A}^{\phi}(a_{t-1}, y_t) = \sigma(W_{\phi} a_{t-1} + U_{\phi} y_t + b_{\phi}), $$

where

* $W_{\phi} \in \mathbb{R}^{d \times d}$ and $U_{\phi} \in \mathbb{R}^{d \times d_y}$ are weight matrices and $b_{\phi} \in \mathbb{R}^d$ is a bias vector.

* $\sigma: \mathbb{R}^d \to \mathbb{R}^d$ is a Lipschitz mapping with Lipschitz constant $L_{\sigma}$. Typical examples include the ReLU and hyperbolic tangent (tanh) functions, both of which satisfy $L_{\sigma} = 1$.
Consider the spectral norm of $W_\phi$ given by

$$L_{W_\phi} = \|W_\phi\|_2 = \sup_{\|x\|_2=1} \|W_\phi x\|_2,$$

i.e., $L_{W_\phi}$ is the largest singular value of $W_\phi$. Now, let us consider two sequences $(a_t)_{t \in \mathbb{N}}$ and $(\tilde{a}_t)_{t \in \mathbb{N}}$ of hidden states starting with $a_0$ and $\tilde{a}_0$, respectively. For every fixed observation sequence $y_{1:t}$, it holds, by definition (49), that

$$\begin{aligned} \|\mathcal{A}_t^\phi(a_0, y_{1:t}) - \mathcal{A}_t^\phi(\tilde{a}_0, y_{1:t})\|_2 &\le L_\sigma \|W(\mathcal{A}_{t-1}^\phi(a_0, y_{1:t-1}) - \mathcal{A}_{t-1}^\phi(\tilde{a}_0, y_{1:t-1}))\|_2 \\ &\le L_\sigma^t \|W\|_2^t \|a_0 - \tilde{a}_0\|_2. \end{aligned}$$

Typically, we may assume here that $L_\sigma = 1$ (again, this is the case for the ReLU and tanh functions). Theorem C.9(i) then holds if $\|W_\phi\|_2 < 1$. This constraint can be enforced during training via spectral norm regularization [Yoshida and Miyato, 2017]. Alternatively, one may enforce it a priori by fixing the spectral norm to be at most $\rho_{\max} < 1$ and, at each gradient step $k$, projecting the unconstrained weight matrix $\tilde{W}_\phi^{(k)}$ onto the corresponding spectral-norm ball according to

$$W_\phi^{(k)} = \rho_{\max} \frac{\tilde{W}_\phi^{(k)}}{\|\tilde{W}_\phi^{(k)}\|_2}.$$

Such projections are already standard in echo state networks [Sun et al., 2020].

**Assumption C.11** (Uniform ergodicity of $(X_t, Y_t)_{t \in \mathbb{N}}$). There exist $\beta \in (0, 1)$ and $\sigma \in \mathbb{M}_1(\mathcal{X} \otimes \mathcal{Y})$ such that for every $(x, y) \in \mathcal{X} \times \mathcal{Y}$ and $t \in \mathbb{N}_{>0}$,

$$\|R^t((x, y), \cdot) - \sigma\|_{\text{TV}} \le \beta^t.$$

In addition, consider the following Lipschitz class of functions.

**Definition C.12.** Let $\tilde{\mathcal{L}}(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ be the set of $f \in \mathcal{F}(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ for which there exists $\tilde{\varphi} \in \mathcal{F}_b(\mathcal{X} \otimes \mathcal{Y})$ such that for every $(x, y, a, a') \in \mathcal{X} \times \mathcal{Y} \times \mathcal{A}^2$,

$$(i) \quad |f(x, y, a)| \le \tilde{\varphi}(x, y),$$

$$(ii) \quad |f(x, y, a) - f(x, y, a')| \le \tilde{\varphi}(x, y) \|a - a'\|_2.$$

**Proposition C.13.** Assume C.9 and C.11. Then there exist a function $\delta: \mathcal{A}^2 \to \mathbb{R}_{>0}$ and $\alpha \in (0, 1)$ such that for every $\phi \in \Phi$, $t \in \mathbb{N}$, $f \in \tilde{\mathcal{L}}(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$, and $((x, y, a), (x', y', a')) \in (\mathcal{X} \times \mathcal{Y} \times \mathcal{A})^2$,

$$|(S^\phi)^t f(x, y, a) - (S^\phi)^t f(x', y', a')| \le \delta(a, a') \|\tilde{\varphi}\|_\infty \alpha^t. \quad (50)$$

Moreover, there exists a kernel $\pi^\phi$ on $(\mathcal{X} \times \mathcal{Y} \times \mathcal{A}) \times \tilde{\mathcal{L}}(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$ such that for every $f \in \tilde{\mathcal{L}}(\mathcal{X} \otimes \mathcal{Y} \otimes \mathcal{A})$, $\pi^\phi f$ is constant and for every $t \in \mathbb{N}$ and $(x, y, a) \in \mathcal{X} \times \mathcal{Y} \times \mathcal{A}$,

$$|(S^\phi)^t f(x, y, a) - \pi^\phi f| \le \bar{\delta}(a) \|\tilde{\varphi}\|_\alpha^t, \quad (51)$$

where

$$\bar{\delta}(a) := \frac{1}{1-\alpha} \int \delta(a, a') S^\phi((x, y, a), d(x', y', a')). \quad (52)$$

**Proof of Theorem C.13.** We proceed as in the proof of Theorem C.5. To prove the first claim, we introduce the short-hand notation

$$f_{s,t}(x_t, y_{s:t}, a) := f(x_t, y_t, A_{t-s+1}^\phi(a, y_{s:t})) \quad (53)$$

for $s \in [1, t]$, where we have omitted the dependence on $\phi$ for brevity. Note that with this notation,

$$(S^\phi)^t f(z) = \int \cdots \int f_{1,t}(x_t, y_{1:t}, a) R((x, y), d(x_1, y_1)) \prod_{s=1}^{t-1} R((x_s, y_s), d(x_{s+1}, y_{s+1})).$$

Now, picking $\tilde{a} \in \mathcal{A}$ arbitrarily and using the decomposition

$$f_{1,t}(x_t, y_{1:t}, a) = f_{1,t}(x_t, y_{1:t}, a) - f_{1,t}(x_t, y_{1:t}, \tilde{a}) + \sum_{s=1}^{t} (f_{s,t}(x_t, y_{s:t}, \tilde{a}) - f_{s+1,t}(x_t, y_{s+1:t}, \tilde{a})) + f_{t,t}(x_t, y_t, \tilde{a}),$$
we may write,

$$
\begin{align}
& (S^{\phi})^t f(x, y, a) - (S^{\phi})^t f(x', y', a') \nonumber \\
&= \int \cdots \int (f_{1,t}(x_t, y_{1:t}, a) - f_{1,t}(x_t, y_{1:t}, \tilde{a})) R((x, y), d(x_1, y_1)) \prod_{s=1}^{t-1} R((x_s, y_s), d(x_{s+1}, y_{s+1})) \nonumber \\
&\quad - \int \cdots \int (f_{1,t}(x_t, y_{1:t}, a') - f_{1,t}(x_t, y_{1:t}, \tilde{a})) R((x', y'), d(x_1, y_1)) \prod_{s=1}^{t-1} R((x_s, y_s), d(x_{s+1}, y_{s+1})) \nonumber \\
&\quad + \sum_{s=1}^t \int \cdots \int (f_{s,t}(x_t, y_{s:t}, \tilde{a}) - f_{s+1,t}(x_t, y_{s+1:t}, \tilde{a})) \Delta R_s((x, y), d(x_s, y_s)) \prod_{\ell=s}^{t-1} R((x_\ell, y_\ell), d(x_{\ell+1}, y_{\ell+1})) \nonumber \\
&\quad - \sum_{s=1}^t \int \cdots \int (f_{s,t}(x_t, y_{s:t}, \tilde{a}) - f_{s+1,t}(x_t, y_{s+1:t}, \tilde{a})) \Delta R_s((x', y'), d(x_s, y_s)) \prod_{\ell=s}^{t-1} R((x_\ell, y_\ell), d(x_{\ell+1}, y_{\ell+1})) \nonumber \\
&\quad + \int f_{t,t}(x_t, y_t, \tilde{a}) \Delta R_t((x, y), d(x_t, y_t)) - \int f_{t,t}(x_t, y_t, \tilde{a}) \Delta R_t((x', y'), d(x_t, y_t)), \tag{54}
\end{align}
$$

where we have defined, for $s \in [1, t]$, the signed kernels

$$
\Delta R_s f(x, y) := R^s f(x, y) - \sigma f, \quad (f, x, y) \in \mathbb{F}_b(\mathcal{X} \otimes \mathcal{Y}) \times \mathcal{X} \times \mathcal{Y}.
$$

Note that by Theorem C.12 and Theorem C.9, for every $s \in [1, t]$,

$$
\begin{align}
|f_{s,t}(x_t, y_{s:t}, \tilde{a}) - f_{s+1,t}(x_t, y_{s+1:t}, \tilde{a})| &\le \tilde{\varphi}(x_t, y_t) \|\mathcal{A}_{t-s+1}^{\phi}(\tilde{a}, y_{s:t}) - \mathcal{A}_{t-s}^{\phi}(\tilde{a}, y_{s+1:t})\|_2 \nonumber \\
&= \tilde{\varphi}(x_t, y_t) \|\mathcal{A}_{t-s}^{\phi}(\mathcal{A}^{\phi}(\tilde{a}, y_s), y_{s+1:t}) - \mathcal{A}_{t-s}^{\phi}(\tilde{a}, y_{s+1:t})\|_2 \nonumber \\
&\le \tilde{\varphi}(x_t, y_t) c \|\mathcal{A}^{\phi}(\tilde{a}, y_s) - \tilde{a}\|_2 \kappa^{t-s} \nonumber \\
&\le \tilde{\varphi}(x_t, y_t) c^2 \kappa^{t-s}. \tag{55}
\end{align}
$$

Similarly,

$$
\begin{align}
|f_{1,t}(x_t, y_{1:t}, a) - f_{1,t}(x_t, y_{1:t}, \tilde{a})| &\le \tilde{\varphi}(x_t, y_t) c \|a - \tilde{a}\|_2 \nonumber \\
&\le \tilde{\varphi}(x_t, y_t) c (\|a\|_2 + \|\tilde{a}\|) \kappa^t. \tag{56}
\end{align}
$$

Now, let

$$
\delta'(a, a') := c^2 \lor c(\|a\|_2 + \|\tilde{a}\|) \lor c(\|a'\|_2 + \|\tilde{a}\|);
$$

then, applying the bounds (55) and (56) to the decomposition (54) yields

$$
\begin{align}
|(S^{\phi})^t f(x, y, a) - (S^{\phi})^t f(x', y', a')| \nonumber \\
&\le \delta'(a, a') \kappa^t (R^t \tilde{\varphi}(x, y) + R^t \tilde{\varphi}(x', y')) \nonumber \\
&\quad + \delta'(a, a') \sum_{s=1}^t \kappa^{t-s} (\Delta R_s R^{t-s} \tilde{\varphi}(x, y) + \Delta R_s R^{t-s} \tilde{\varphi}(x', y')) \nonumber \\
&\quad + \Delta R_t \tilde{\varphi}(x, y) + \Delta R_t \tilde{\varphi}(x', y'). \tag{57}
\end{align}
$$

Now, by applying Theorem C.11 to the right-hand side of (57) we obtain

$$
|(S^{\phi})^t f(x, y, a) - (S^{\phi})^t f(x', y', a')| \le 2 \|\tilde{\varphi}\|_{\infty} (\delta'(a, a') \kappa^t + \delta'(a, a') t(\kappa \lor \beta)^t + \beta^t),
$$

from which (50) follows by picking $\alpha \in ((\kappa \lor \beta), 1)$ and letting

$$
\delta(a, a') := 2(\delta'(a, a') + 1)(\kappa \lor \beta)/\alpha + 2\delta'(a, a')(e | \log\{(\kappa \lor \beta)/\alpha\}|)^{-1}.
$$

To prove the second claim, first note that by (50),

$$
\begin{align}
& |(S^{\phi})^{t+1} f(x, y, a) - (S^{\phi})^{t} f(x, y, a)| \nonumber \\
&\le \int |(S^{\phi})^{t} f(x', y', a') - (S^{\phi})^{t} f(x, y, a)| S^{\phi}((x, y, a), d(x', y', a')) \nonumber \\
&\le \|\tilde{\varphi}\|_{\infty} \alpha^t \int \delta(a, a') S^{\phi}((x, y, a), d(x', y', a')). \tag{58}
\end{align}
$$
Now, define the kernel

$$ \pi^{\phi} f(x, y, a) := f(x, y, a) + \sum_{t=0}^{\infty} ((S^{\phi})^{t+1} f(x, y, a) - (S^{\phi})^{t} f(x, y, a)) $$

on $X \times Y \times A \times \mathcal{L}(X \otimes Y \otimes A)$. With this definition, note that by (58),

$$ \begin{aligned} |(S^{\phi})^{t} f(x, y, a) - \pi^{\phi} f(x, y, a)| &\le \sum_{s=t}^{\infty} |(S^{\phi})^{s+1} f(x, y, a) - (S^{\phi})^{s} f(x, y, a)| \\ &\le \delta(a) \|\tilde{\phi}\|_{\infty} \alpha^{t}, \end{aligned} $$

where $\delta(a)$ is provided by (52), which establishes (51).

Finally, it remains to prove that the function $\pi^{\phi} f$ is constant. For this purpose, pick arbitrarily $((x, y, a), (x'y', a')) \in (X \times Y \times A)^2$; then, however, by (50) and (51),

$$ \begin{aligned} |\pi^{\phi} f(x, y, a) - \pi^{\phi} f(x', y', a')| &\le \inf_{t \in \mathbb{N}} \left( |(S^{\phi})^{t+1} f(x, y, a) - (S^{\phi})^{t} f(x, y, a)| \right. \\ &\qquad \left. + |(S^{\phi})^{t} f(x, y, a) - \pi^{\phi} f(x, y, a)| \right. \\ &\qquad \left. + |(S^{\phi})^{t} f(x', y', a') - \pi^{\phi} f(x', y', a')| \right) = 0, \end{aligned} $$

from which the claim follows. $\square$

# D Full algorithm using backward sampling and control variate

In this section, we detail Algorithm 3, which presents one iteration of the online gradient ascent algorithm in the amortized scheme, incorporating both backward sampling and control variates for improved efficiency and variance reduction.

**Backward sampling.** Computing the backward weights of (20) has the disadvantage of $O(N^2)$ complexity due to the computation of the normalizing constant, which can be prohibitive when $N$ is large (typically for high-dimensional state spaces). One solution, suggested by [Olsson and Westerborn, 2017] in the context of SMC smoothing, is to use a backward sampling approach. More precisely, at time step $t$, given $\xi_t^i$, one samples independently $M$ indexes $\{j_k\}_{k=1}^M$ from the categorical distribution over $\{1, \dots, N\}$ with weights $\{\bar{w}_{t-1|t}^{\phi, i, j}\}_{j=1}^N$, and replace (21) by $\sum_{k=1}^M (\hat{h}_{t-1}^{\phi, i, j_k} + \ell_t^{\phi, \phi_t}(\xi_{t-1}^{j_k}, \xi_t^i))/M$. [Olsson and Westerborn, 2017] show that even with $M$ much smaller than $N$ (typically, $M=2$), which provides a considerable improvement in complexity, this alternative estimator has only slightly higher variance than the original estimator. Here, noting that $\bar{w}_{t-1|t}^{\phi, i, j} \propto_j \psi_t^{\phi}(\xi_{t-1}^{j}, \xi_t^i)$, backward sampling can, in the case of bounded potential functions, be performed using an accept-reject procedure without having to calculate the normalizing constant of the weights. We refer the reader to [Olsson and Westerborn, 2017, Gloaguen et al., 2022, Dau and Chopin, 2022] for details and alternative backward sampling approaches.

**Variance reduction of the gradient estimator.** Proposition 4.1 involves computing score-function expectations in the form $\mathbb{E}_{q^{\phi}}[\nabla_{\phi} \log q^{\phi}(X) \cdot f(X)]$ for some p.d.f. $q^{\phi}$. As shown in [Mohamed et al., 2020], direct Monte Carlo estimation of the score function leads to high variance and should normally not be used without a suitable variance reduction technique. The most straightforward approach is to design a control variate. Using the fact that $\mathbb{E}_{q^{\phi}}[\nabla_{\phi} \log q^{\phi}(X)] = 0$, the target expectation can be rewritten as $\mathbb{E}_{q^{\phi}}[\nabla_{\phi} \log q^{\phi}(X)\{f(X) - \mathbb{E}_{q^{\phi}}[f(X)]\}]$, which can be estimated with lower variance using a Monte Carlo estimate of $\mathbb{E}_{q^{\phi}}[f(X)]$. In our case, the latter is formed as a by-product of Algorithm 2, and therefore our methodology comes with built-in variance reduction without the need to recompute additional quantities. This accelerated version of Algorithm 2 is described in detail in Algorithm 3 (see Appendix D), which also includes the backward sampling technique described above.

As an alternative to this variance reduction technique, it is natural to consider the reparametrization trick, as it often leads to Monte Carlo estimators with lower variance compared to those obtained using the score function. However, the implementation of the reparametrization trick in this context requires that $\nabla_{\phi} \mathcal{L}_t^{\phi, \phi}$ is expressed as an expectation with respect to a random variable $Z_{0:t}$ that does not depend on $\phi$. Moreover, the recursive expression of this expectation at time $t+1$ must be derivable from its predecessor, which is non-trivial. For example, in the classical case where $q_{0:t}$ is the p.d.f. of a multivariate Gaussian random variable with mean $\mu$ and variance $\Sigma$, and the expectation is taken w.r.t. $Z_{0:t} \sim \mathcal{N}(0, I_{d_x \times (t+1)})$, such a recursion is not feasible as the ELBO is no longer an additive functional when $X_{0:t}$ is replaced by $\mu + \Sigma^{\frac{1}{2}} Z_{0:t}$.
<table><caption>Table 2: RMSE of the predicted marginal means $\mathbb{E}_{q_{0:T}^{\phi}}[X_t]$ w.r.t. the true states $x_t^*$ and average time per gradient step.</caption><thead><tr><th>Gradients</th><th>Δ<sup>φ</sup><sub>T,p</sub> (×10<sup>-2</sup>)</th><th>Avg. time</th></tr></thead><tbody><tr><td>Score-based</td><td>13.5 ± 0.7 (<b>12.2</b>)</td><td>173 ms</td></tr><tr><td>Backward sampling</td><td>11.9 ± 0.4 (<b>11.4</b>)</td><td>17 ms</td></tr></tbody></table>

# E Supplementary details for the numerical experiments in Section 6

## E.1 Appendix for section 6.1 the linear Gaussian SSM.

**Variational Family.** In the linear Gaussian setting, the variational marginals $q_t^\phi$ are parameterized as Gaussian distributions defined by their natural parameters $\eta_t^\phi$. To ensure the backward kernel $q_{t-1|t}^\phi$ remains in the same Gaussian family, the potential is explicitly defined as

$$ \psi_t^\phi(x_{t-1}, x_t) = \exp(\langle \tilde{\eta}_t^\phi(x_t), T(x_{t-1}) \rangle). \qquad (59) $$

This formulation allows the backward density to be derived analytically by simply summing natural parameters, avoiding the need for normalizing constants. The parameters for the backward kernel are updated according to:

$$ \eta_{t-1|t}^{\phi} = \eta_{t-1}^{\phi} + \tilde{\eta}_{t}^{\phi}. \qquad (60) $$

Unlike the general case requiring neural networks, the recursions for these parameters in the linear Gaussian case are analytical, effectively mirroring the smoothing distribution updates of a standard linear Gaussian SSM.

**Parameters for the linear Gaussian SSM.** For the streaming experiment presented in the main text ($T = 50,000$), we learn both model and variational parameters from random initialization. We utilize learning rates of $10^{-3}$ for the variational parameters and $10^{-4}$ for the model parameters. To ensure a rigorous comparison, we replicate the generative settings of [Campbell et al., 2021], using diagonal noise covariance matrices with fixed variances of 0.1 for the transition and 0.25 for the emission.

**Oracle ELBO** As an oracle baseline, we can compute the closed-form ELBO and its associated gradient via the reparameterization trick. Figure 5 displays the evolution of the ELBO in the case of the linear SSM and the offline setting, i.e., when observations are processed through multiple epochs. In this specific experiment, we do not learn the model parameter. For our recursive method, we choose $\Delta = 2$ to truncate the backpropagation, as we observe that $\Delta < 2$ prevents our method from converging altogether, while $\Delta > 2$ only improves convergence speed by a small margin. The experiment is run using 10 different parameters for the generative model, $d_x = d_y = 10$, $T = 500$, and $N = 2$ for the two methods involving Monte Carlo sampling. It shows the convergence of our score-based solution to the correct optimum given by the analytical computations. This is particularly appealing and notably demonstrates that our online gradient-estimation method may perform well using few samples. In practice, we observe that the variance reduction introduced in Section D is crucial in reaching such performance.

## E.2 Appendix for section 6.2 the chaotic RNN

The 1-step smoothing and filtering errors in table 1 are given as

$$ \kappa_T^{(1)} = \frac{1}{T-1} \sum_{t=1}^{T-1} \left( \frac{1}{d_x} \sum_{k=1}^{d_x} \left( \hat{\mathbb{E}}_{q_{t-1:t}^{\phi}} \left[ X_{t-1}^{(k)} \right] - x_{t-1}^{*(k)} \right)^2 \right)^{1/2} \qquad (61) $$

$$ \kappa_T^{(2)} = \frac{1}{T} \sum_{t=1}^{T} \left( \frac{1}{d_x} \sum_{k=1}^{d_x} \left( \hat{\mathbb{E}}_{q_t^{\phi}} \left[ X_t^{(k)} \right] - x_t^{*(k)} \right)^2 \right)^{1/2}, \qquad (62) $$

where $\hat{\mathbb{E}}_q$ denotes the standard Monte Carlo estimate of an expectation w.r.t. $q$. Table 3 reports smoothing and filtering RMSE with respect to the true states at the end of training.
<table><caption>Table 3: Smoothing and filtering RMSE when $\phi$ is learned online together with $(\rho, \gamma)$ in the chaotic RNN. Results are shown for the training stream and for independent sequences from the same generative model.</caption><thead><tr><th>Sequence</th><th>Smoothing RMSE</th><th>Filtering RMSE</th></tr></thead><tbody><tr><td>Training</td><td>0.281</td><td>0.311</td></tr><tr><td>Eval</td><td>0.278 (± 0.01)</td><td>0.305 (± 0.014)</td></tr></tbody></table>

**Parameters for the chaotic RNN** We choose the same generative hyperparameters as [Campbell et al., 2021] with $\Delta = 0.001$, $\tau = 0.025$, $\gamma = 2.5$, 2 degrees of freedom and a scale of 0.1 for the Student-t distribution, and define $Q = \text{diag}(0.01)$. For the joint learning setting, training is performed using the Adam [Kingma and Ba, 2015] optimizer with learning rates of $10^{-3}$ for the variational parameters and $10^{-4}$ for the model parameters. The variational backward kernels follow the parameterization defined in the main text, where the potential functions are parameterized by MLPs with 100 hidden units and tanh activations.

**Implementation settings for the comparison with [Campbell et al., 2021]** In the original paper, the variational distribution is designed in a non-amortized scheme, meaning that the variational parameters are not shared through time. Specifically, we have that $\phi = \{\phi^0, \dots, \phi^t, \dots\}$, and each $\phi^t$ contains the parameter $\eta_t = (\mu_t, \Sigma_t)$ of the distribution $q_t^\phi \sim \mathcal{N}(\mu_t, \Sigma_t)$ and the parameter $\tilde{\eta}_t$ of the function $\psi_t^\phi$. For this latter function, we match the number of parameters of [Campbell et al., 2021] by defining $\psi_t^\phi(x, y) = \exp(\tilde{\eta}_t(y) \cdot T(x))$ with $\tilde{\eta}_t(y) = (\tilde{\eta}_{t,1}(y), \tilde{\eta}_{t,2})$ where $y \mapsto \tilde{\eta}_{t,1}(y)$ is a multi-layer perceptron with 100 neurons from $\mathbb{R}^{d_x}$ to $\mathbb{R}^{d_x}$, and $\tilde{\eta}_{t,2}$ is a negative definite matrix. We follow the optimization schedules of [Campbell et al., 2021] with $K = 500$ gradient steps at each time-step.

**Modifications induced by the non-amortized scheme** In the non-amortized scheme, $\phi = \{\phi^0, \dots, \phi^t\}$ is a set of distinct parameters, each parameter corresponding to a specific time step. In the notations of the article, the estimate $\phi_{t-1}$ of $\phi$ after having processed observations the $y_{0:t-1}$ depends on $\{\phi^0, \dots, \phi^{t-1}\}$. Therefore, the gradient of the ELBO will be w.r.t. $\phi^t$ only. This affects the expression of the statistic $u_t$, and one can see in the expansion of Eqn. (27) that the term (28) will now, when the gradient is taken w.r.t. $\phi^t$, be zero. This means that this term no longer has to be propagated. Indeed, as we set $q_{t-1|t}^\phi(x_t, x_{t-1}) \propto q_{t-1}^{\phi_{t-1}}(x_{t-1})\psi_t^\phi(x_{t-1}, x_t)$, the gradient of the ELBO w.r.t. $\phi^t$ will be

$$ \nabla_{\phi^t} \mathcal{L}_t^{\theta, \phi^t} = \mathbb{E}_{q_t^{\phi^t}} \left[ \mathbb{E}_{q_{t-1|t}^{\phi^t}} \left[ \nabla \log q_t^{\phi}(X_t) \times \ell_t^{\theta, \phi^t}(X_{t-1}, X_t) \right. \right. \\ \left. \left. + \nabla_{\phi^t} \log q_{t-1|t}^{\phi^t}(X_{t-1}, x_t) \times (h_{t-1}(X_{t-1}) + \ell_t^{\theta, \phi^t}(X_{t-1}, X_t)) \right] \right]. $$

This gradient will be estimated using Monte Carlo in the same way as in the algorithm. In [Campbell et al., 2021], the inner conditional expectation is estimated with the regression approach (that we briefly recall below) instead of our importance sampling approach.

**Functional regression approach of [Campbell et al., 2021]** Here, we recall the alternate option used in [Campbell et al., 2021] to propagate approximations of the backward expectations. Denoting $\mathcal{F} = \{g : \mathbb{R}^p \to \mathbb{R}^{d_x}, \mathbb{E}_{q_t^\phi}[\|g(X_t)\|_2] < \infty\}$, $h_t(x)$ satisfies, by definition of conditional expectation,

$$ h_t = \underset{g \in \mathcal{F}}{\operatorname{argmin}} \mathbb{E}_{q_{t-1:t}^{\phi}} \|g(X_t) - [h_{t-1}(X_{t-1}) + \ell_t^{\theta, \phi}(X_{t-1}, X_t)]\|_2, $$

which provides a regressive objective for learning an approximation of $h_t$. In practice, the authors restrict the minimization problem to a subset of $\mathcal{F}$, a parametric family of functions (typically, a neural network) parameterized by $\gamma$ belonging to $\Gamma \subset \mathbb{R}^{d_\gamma}$, and learn this by approximating the expectation using Monte Carlo sampling. More precisely, the authors propose to estimate $h_t$ by $H_{\hat{\gamma}_t}^\phi$, where

$$ \hat{\gamma}_t = \underset{\gamma \in \Gamma}{\operatorname{argmin}} \frac{1}{N} \sum_{k=1}^{N} \| H_{\gamma}^{\phi}(\xi_k^k) - [H_{\hat{\gamma}_{t-1}}^{\phi}(\xi_{t-1}^k) + \tilde{h}_t(\xi_{t-1}^k, \xi_k^k)] \|_2, \quad (63) $$

where $\{(\xi_{t-1}^i, \xi_t^i)\}_{i=1}^N$ is an i.i.d. sample from the variational joint distribution of $(X_{t-1}, X_t)$, which has density $q_{t-1:t}^\phi = q_t^\phi q_{t-1|t}^\phi$. Upon convergence, $H_{\hat{\gamma}_t}^\phi$ is then used in the successive recursions (similar to (21)).
E.3 Appendix for section 6.3

**Air-quality data** The UCI Air Quality dataset contains hourly readings from 5 metal oxide chemical sensors located in a significantly polluted area. Our observation vector $y_t$ comprises Carbon Monoxide (CO), Non-Methane Hydrocarbons (NMHC), Nitrogen Oxides ($NO_x$, $NO_2$), Benzene ($C_6H_6$), Temperature (T), Relative Humidity (RH), and Absolute Humidity (AH). A key challenge of this dataset is the non-uniform data censorship caused by sensor failures. As shown in Figure 6, these dropouts occur at different intervals for different sensors. In the raw data, these are marked by a sentinel value (−200), which we treat as missing values (NaNs) during preprocessing.

**Parameters for the Air Quality Data** The transition function $f_\theta$ and emission function $g_\theta$ are both parameterized by MLPs with 2 hidden layers of 32 units each and tanh activation functions. The noise covariances $Q_\theta$ and $R_\theta$ are learned diagonal matrices. Regarding the variational approximation, the backward potentials are parameterized by MLPs with 32 hidden units and tanh activations. We use $N = 20$ importance samples for the Monte Carlo estimates. Training is performed using the Adam optimizer with learning rates of $10^{-3}$ for the variational parameters and $10^{-4}$ for the model parameters.
<figure><img src="image_6.png" /><figcaption>Figure 4: **Simultaneous Variational Learning and Prediction.** (a) Smoothing reconstruction during sensor failure (black dotted). RMCVI (orange) recovers the Truth (gray) while LSTM (blue) overfits the corruption. (b) Cumulative RMSE for one-step-ahead prediction averaged across all pollution features.</figcaption></figure>
**Algorithm 3** One iteration of the online gradient ascent algorithm (for $t \ge 1$) in the amortized scheme

**Require:**

* Previous statistics $\{(\hat{h}_{t-1}^{\phi_{t-1},i}, \hat{u}_{t-1}^{\phi_{t-1},i}, \hat{v}_{t-1}^{\theta_{t-1},i})\}_{i=1}^N$, and previous samples $\{\xi_{t-1}^i\}_{i=1}^N$;

* Intermediate quantity $a_{t-1}$, parameter estimates $(\theta_t, \phi_t)$, step sizes $(\gamma_t^\theta, \gamma_t^\phi)$;

* New observation $y_t$.

**Ensure:** $\{(\hat{u}_t^{\phi_t, i}, \hat{h}_t^{\phi_t, i})\}_{i=1}^N, \phi_{t+1}, a_t$.

Set $a_t \leftarrow \mathcal{A}^{\phi_t}(a_{t-1}, y_t)$ and $\eta_t^{\phi_t} \leftarrow f^{\phi_t}(a_t)$, the parameters of $q_t^{\phi_t}$

Sample $\{\xi_t^i\}_{i=1}^N$ independently from $q_t^{\phi_t}$

**for** $i \leftarrow 1$ to $N$ **do**

Set $\tilde{\eta}_t^{\phi_t, i} \leftarrow \tilde{f}^{\phi_t}(\xi_t^i)$

**for** $j \leftarrow 1$ to $M$ **do**

// Backward sampling step, $M$ is the number of backward samples

Sample $(j_k)_{k=1}^M$ i.i.d. $\text{Cat}(\{\bar{w}_{t-1|t}^{\phi_t, j}\}_{1 \le j \le N})$ with the weights (20).

**end for**

Set // Recall that each term $\ell_t^{\theta, \phi_t}$ depends on $y_t$.

$$ \hat{h}_t^{\phi_t, i} \leftarrow \frac{1}{M} \sum_{k=1}^{M} \left( \hat{h}_{t-1}^{\phi_{t-1}, j_k} + \ell_t^{\theta, \phi_t} (\xi_{t-1}^{j_k}, \xi_t^i) \right); $$

$$ \hat{u}_t^{\phi_t, i} \leftarrow \frac{1}{M} \sum_{k=1}^{M} \left\{ \hat{u}_{t-1}^{\phi_{t-1}, j_k} + \nabla_{\phi} \log q_{t-1|t}^{\phi_t} (\xi_{t-1}^{j_k}, \xi_t^i) \left( \hat{h}_{t-1}^{\phi_{t-1}, j_k} + \ell_t^{\theta, \phi_t} (\xi_{t-1}^{j_k}, \xi_t^i) - \hat{h}_t^{\phi_t, i} \right) \right\}; $$

$$ \hat{v}_t^{\theta_t, i} \leftarrow \sum_{i=1}^{N} \bar{w}_{t-1|t}^{\phi_t, i, j} \left\{ \hat{v}_{t-1}^{\theta_{t-1}, j} + \nabla_{\theta} \ell_t^{\theta, \phi_t} (\xi_{t-1}^j, \xi_t^i) \right\}; $$

// Note the difference with (22) and the inclusion of control variate $\hat{h}_t^{\phi_t, i}$ for the computation of $\hat{u}_t^{\phi_t, i}$

// $\nabla q_{t-1|t}^{\phi_t}(\xi_{t-1}^{j_k}, \xi_t^i)$ is typically computed with automatic differentiation

**end for**

Set

$$ \hat{\nabla}_{\phi} \mathcal{L}_{t}^{\theta_t, \phi_t} \leftarrow \frac{1}{N} \sum_{i=1}^{N} \left\{ \hat{u}_{t}^{\phi_{t}, i} + \nabla_{\phi} \log q_{t}^{\phi_{t-1}}(\xi_{t}^{i}) \left( \hat{h}_{t}^{\phi_{t}, i} - \frac{1}{N} \sum_{k=1}^{N} \hat{h}_{t}^{\phi_{t}, k} \right) \right\}; $$

$$ \phi_{t+1} \leftarrow \phi_t + \gamma_t^{\phi} \left( \hat{\nabla}_{\phi} \mathcal{L}_{t}^{\theta_t, \phi_t} - \hat{\nabla}_{\phi} \mathcal{L}_{t-1}^{\theta_{t-1}, \phi_{t-1}} \right) $$

$$ \hat{\nabla}_{\theta} \mathcal{L}_{t}^{\theta_t, \phi_t} \leftarrow \frac{1}{N} \sum_{i=1}^{N} \hat{v}_{t}^{\theta_t, i}; $$

$$ \theta_{t+1} \leftarrow \theta_t + \gamma_{t+1}^{\theta} \left( \hat{\nabla}_{\theta} \mathcal{L}_{t}^{\theta_t, \phi_t} - \hat{\nabla}_{\theta} \mathcal{L}_{t-1}^{\theta_{t-1}, \phi_{t-1}} \right). $$

// Note the difference with Algorithm 2 and the inclusion of the control variate $\frac{1}{N} \sum_{i=1}^N \hat{h}_t^{\phi_t, i}$ in the calculation of $\hat{\nabla}_{\phi} \mathcal{L}_{t}^{\theta_t, \phi_t}$

// $\nabla \log q_t^{\phi_{t-1}}(\xi_t^i)$ is typically computed using automatic differentiation
<figure><img src="image_7.png" /><figcaption>Figure 5: Evolution of $\mathcal{L}_T^{\theta, \phi}/T$ computed with three different methods and with three different types of gradients estimates. Full lines: means of the 10 replicates. Shaded lines: standard deviations of the 10 replicates.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Time series of the 8 observed dimensions in the UCI Air Quality dataset. The gaps in the signals correspond to periods of sensor censorship (missing data). The framework must infer the latent state during these blackout periods without direct observation.</figcaption></figure>