Amortized Control of Continuous State Space Feynman-Kac Model for Irregular Time Series

## Abstract

Abstract Many real-world datasets, such as healthcare, climate, and economics, are often collected as irregular time series, which poses challenges for accurate modeling. In this paper, we propose the Amortized Control of continuous State Space Model (ACSSM) for continuous dynamical modeling of time series for irregular and discrete observations. We first present a multi-marginal Doob’s h ℎ h italic_h -transform to construct a continuous dynamical system conditioned on these irregular observations. Following this, we introduce a variational inference algorithm with a tight evidence lower bound (ELBO), leveraging stochastic optimal control (SOC) theory to approximate the intractable Doob’s h ℎ h italic_h -transform and simulate the conditioned dynamics. To improve efficiency and scalability during both training and inference, ACSSM employs amortized inference to decouple representation learning from the latent dynamics. Additionally, it incorporates a simulation-free latent dynamics framework and a transformer-based data assimilation scheme, facilitating parallel inference of the latent states and ELBO computation. Through empirical evaluations across a variety of real-world datasets, ACSSM demonstrates superior performance in tasks such as classification, regression, interpolation, and extrapolation, while maintaining computational efficiency.

## 1 Introduction

State space models (SSMs) are widely used for modeling time series data (Fraccaro et al., ; Rangapuram et al., ; de Bézenac et al., ) but typically assume evenly spaced observations. However, many real-world datasets, such as healthcare (Silva et al., ), climate (Menne et al., ), are often collected as irregular time series, which poses challenges for accurate modeling. It provokes the neural differential equations (Chen et al., ; Rubanova et al., ; Li et al., ; Kidger et al., ; Zeng et al., ) for a continuous-time dynamical modeling of a time-series. Moreover, by incorporating stochastic dynamical models with SSM models, the application of continuous-dynamical state space models (CD-SSM) (Jazwinski, ) with neural networks have emerged which allows for more flexible modeling (Schirmer et al., ; Ansari et al., ). However, they typically rely on the Bayesian recursion to approximate the desired latent posterior distribution (filtering/smoothing distribution) involves sequential prediction and updating, resulting in computational costs that scale with the length of observations (Särkkä & García-Fernández, ).

Instead of relying on Bayesian recursion, we propose an Amortized Control of continuous State Space Model (ACSSM), a variational inference (VI) algorithm that directly approximates the posterior distribution using a measure-theoretic formulation, enable a simulation of continuous trajectories from such distributions. To achieve this, we introduce a Feynman-Kac model (Del Moral, ; Chopin et al., ), which facilitates a thorough sequential analysis. Within this formulation, we propose a multi-marginal Doob’s $h$-transform to provide conditional latent dynamics described by stochastic differential equations (SDEs) which conditioned on a irregularly sampled, set of discrete observations. It extends the previous Doob’s $h$-transform (Chopin et al., ) by incorporating multi-marginal constraints instead of a single terminal constraint.

Considering that the estimation of the Doob’s $h$-transform is infeasible in general, we utilize the theory of stochastic optimal control (SOC) (Fleming & Soner, ; Carmona, ) to simulate the conditioned SDEs by approximating the Doob’s $h$-transform. It allow us to propose a tight evidence lower bound (ELBO) for the aforementioned VI algorithm by establishing a fundamental connection between the partial differential equations (PDEs) associated with Doob’s $h$-transform and SOC. In the Sequential Monte Carlo (SMC) literatures, Doob’s $h$-transform is analogous to the twist-function (Guarniero et al., ). In this context, Heng et al. () proposed an algorithm focused on approximating the twisted transition kernel instead of the underlying dynamics. Recently, a concurrent study Lu & Wang () presented a continuous-time expansion. However, both studies primarily concentrate on approximation methods rather than practical real-world applications.

In practical situations, the computation of ELBO for a VI algorithm might impractical due to the instability and high memory demands associated with gradient computation of the approximated stochastic dynamics over the entire sequence interval (Liu et al., ; Park et al., ). To address this issue, we propose two efficient modeling approaches:
1) We establish amortized inference by introducing an auxiliary variable to the latent space, generated by a neural network encoder-decoder. It separates the learning of representations in high-dimensional time-series from generative modeling in the latent dynamics (Fraccaro et al., ), offering a more flexible parameterization. Moreover, it allows inference of the posterior distribution for a novel time-series sequence without relying on Bayesian recursion by incorporating the learned control function.
2) We use the simulation-free property, which allows for parallel computation of latent marginal distributions over a given interval. Moreover, we explore a more flexible linear approximation of the drift function in controlled SDEs to enhance the efficiency of the proposed controlled dynamics.

We evaluated ACSSM on several time-series tasks across various real-world datasets. Our experiments show that ACSSM consistently outperforms existing baseline models in each tasks, demonstrating its effectiveness in capturing the underlying dynamics of irregular time-series. Additionally, ACSSM achieves significant computational efficiency, enabling faster training times compared to dynamics-based models that rely on numerical simulations.

We summarize our contributions as follows:

- •
We extend the theory of Doob’s $h$-transform, previsoly studied with a single marginal constraint, to a multi-marginal cases. This indicates the existence of a class of conditioned SDEs that depend on future observations, where the solutions of these SDEs lead to the true posterior path measure within the framework of CD-SSM.
- •
We reformulate the simulation of conditioned SDEs as a SOC problem to approximating an impractical Doob’s $h$-transform. By leveraging the connection between SOC theory and Doob’s $h$-transform, we propose a variational inference algorithm with a tight ELBO.
- •
For practical real-world applications, we introduce an efficient and scalable modeling approach that enables parallelization of latent dynamic simulation and ELBO computation.
- •
We demonstrate its superior performance across various real-world irregularly sampled time-series tasks, including per-time point classification, regression, and sequence interpolation and extrapolation, all with computational efficiency.

#### Notation

Throughout this paper, we denote path measure by $\mathbb{P}^{(\cdot)}$, defined on the space of continuous functions $\Omega=C([0,T],\mathbb{R}^{d})$. We sometimes denote with $\mathbb{P}$ the expectation as $\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}[\cdot]=\mathbb{E}_{\mathbb{P}}[\cdot|
\mathbf{X}_{t}=\mathbf{x}]$, where the stochastic processes corresponding to $\mathbb{P}^{(\cdot)}$ are represented as $\mathbf{X}^{(\cdot)}$ and their time-marginal distribution at time $t\in[0,T]$ is given by the push-forward measure $\mu^{(\cdot)}_{t}:=(\mathbf{X}^{(\cdot)}_{t})_{\#}\mathbb{P}^{(\cdot)}$ with marginal density $\mathbf{p}^{(\cdot)}_{t}$. This marginal density represents the Radon-Nikodym derivative $d\mu_{t}^{(\cdot)}(\mathbf{x})=\mathbf{p}^{(\cdot)}_{t}(\mathbf{x})d\mathcal{L
}(\mathbf{x})$, where $\mathcal{L}$ denotes the Lebesgue measure. Additionally, for a function $\mathcal{V}:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$, we define the first and second derivatives with respect to $\mathbf{x}\in\mathbb{R}^{d}$ as $\nabla_{\mathbf{x}}\mathcal{V}$ and $\nabla_{\mathbf{xx}}\mathcal{V}$, respectively, and the derivative with respect to time $t\in[0,T]$ as $\partial_{t}\mathcal{V}$. For a sequence of functions $\{\mathcal{V}_{i}\}_{i\in[1:k]}$, we will denote $\mathcal{V}_{i}(t,\mathbf{x}):=\mathcal{V}_{i,t}$ and $[1:k]=\{1,\cdots,k\}$. Finally, the Kullback-Leibler (KL) divergence between two probability measures $\mu$ and $\nu$ is defined as $D_{\text{KL}}(\mu|\nu)=\int_{\mathbb{R}^{d}}\log\frac{d\mu}{d\nu}(\mathbf{x})d
\mu(\mathbf{x})$ when $\mu$ is absolutely continuous with respect to $\nu$, and $D_{\text{KL}}(\mu|\nu)=+\infty$ otherwise.

## 2 Preliminaries

#### Continuous-Discrete State Space Model

Consider for a set of time-stamps (regular or irregular) $\{t_{i}\}_{i=0}^{k}$ over a interval $\mathcal{T}=[0,T]$, $\textit{i}.\textit{e}.,0=t_{0}\leq\cdots\leq t_{k}=T$. The CD-SSM assumes a continuous-time Markov state trajectory $\mathbf{X}_{0:T}$ in latent space $\mathbb{R}^{d}$ is given as a solution of the SDE:

$$ $\textit{(Prior State)}\quad d\mathbf{X}_{t}=b(t,\mathbf{X}_{t})dt+d\mathbf{W}_ {t},$ (1) $$

where $\mathbf{X}_{0}\sim\mu_{0}$ and $\{\mathbf{W}_{t}\}_{t\in[0,T]}$ is a $\mathbb{R}^{d}$-valued Wiener process that is independent of the $\mu_{0}$. Since $\mathbf{X}_{t}$ is Markov process, the time-evolution of marginal distribution $\mu_{t}$ is governed by a transition density, which is the solution to the Fokker-Planck equation assocaited with $\mathbf{X}_{t}$. This allows us to define a path measure $\mathbb{P}$ that represent the weak solutions of the SDE in (1) over an interval $[0,T]$(^1^11See details on Appendix A.2.).

For a measurement model $g_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}_{t_{i}})$, we consider the case that we have only access to the realization of the (latent) observation process at each discrete-time stamps $\{t_{i}\}_{i\in[1:k]}$, $\textit{i}.\textit{e}.,\mathbf{y}_{t_{i}}\sim g_{i}(\mathbf{y}_{t_{i}}|\mathbf
{X}_{t_{i}}),\forall i\in[1:k]$. In this paper, our goal is to infer the classes of SDEs which inducing the filtering/smoothing path measure $\mathbb{P}^{\star}:=\mathbb{P}^{\star}(\cdot|\mathcal{H}_{t_{k}})$, the conditional distribution over the interval $[0,T]$ for a given $\mathbb{P}$ and a set of observations up to time $t_{k}$, $\mathcal{H}_{t_{k}}=\{\mathbf{y}_{t_{i}}|i\leq k\}$:

$$ $\displaystyle\textit{(Posterior Dist.)}\quad d\mathbb{P}^{\star}(\mathbf{X}_{0 :T}|\mathcal{H}_{t_{k}})$ $\displaystyle=\frac{1}{\mathbf{Z}(\mathcal{H}_{t_{k}})}\prod_{i=1}^{K}g_{i}( \mathbf{y}_{t_{i}}|\mathbf{X}_{t_{i}})d\mathbb{P}(\mathbf{X}_{0:T})$ (2) $$

where the normalizing constant $\mathbf{Z}(\mathcal{H}_{t_{k}})=\mathbb{E}_{\mathbb{P}}\left[\prod_{i=1}^{K}g_
{i}(\mathbf{y}_{t_{i}}|\mathbf{X}_{t_{i}})\right]$ serve as a observations likelihood. The path measure formulation of the posterior distribution described in (2) referred to as Feynman-Kac models. See (Del Moral, ; Chopin et al., ) for a more comprehensive understanding.

## 3 Controlled Continuous-Discrete State Space Model

In this section, we introduce our proposed model, ACSSM. First, we present the Multi-marginal Doob’s $h$-transform, outlining the continuous dynamics for $\mathbb{P}^{\star}$ in Section 3.1. Then, in Section 3.2, we frame the VI for approximating $\mathbb{P}^{\star}$ using SOC. To support scalable real-world applications, we discuss efficient modeling and amortized inference in Section 3.3(^2^22Proofs and detailed derivations are provided in Appendix A.).

### 3.1 Multi Marginal Doob’s h ℎ h italic_h -transform

Before applying VI to approximate the posterior distribution $\mathbb{P}^{\star}$ in SOC, we first show that a class of SDEs exists whose solutions induce a path measure equivalent to $\mathbb{P}^{\star}$ in (2). This formulation provides a valuable insight for defining an appropriate objective function for the SOC problem in the next section. To do so, we first define a sequence of normalized potential functions $\{f_{i}\}_{i\in[1:k]}$, where each $f_{i}:\mathbb{R}^{d}\to\mathbb{R}_{+}$, for all $i\in[1:k]$,

$$ $f_{i}(\mathbf{y}_{t_{i}}|\mathbf{x}_{t_{i}})=\frac{g_{i}(\mathbf{y}_{t_{i}}| \mathbf{x}_{t_{i}})}{\mathbf{L}_{i}(g_{i})},$ (3) $$

where $\mathbf{L}_{i}(g_{i})=\int_{\mathbb{R}^{d}}g_{t_{i}}(\mathbf{y}_{t_{i}}|
\mathbf{x}_{t_{i}})d\mathbb{P}(\mathbf{x}_{0:T})$, for all $i\in[1:k]$ is the normalization constant. Then, we can observe that the potential functions $\{f_{i}\}_{i\in[1:k]}$ defined in (3) satisfying the normalizing property $\textit{i}.\textit{e}.,\mathbb{E}_{\mathbb{P}}[\prod_{i=1}^{k}f_{i}(\mathbf{y}
_{t_{i}}|\mathbf{x}_{t_{i}})]=1$ and $d\mathbb{P}^{\star}(\mathbf{x}_{0:T}|\mathcal{H}_{t_{k}})=\prod_{i=1}^{k}f_{i}
(\mathbf{y}_{t_{i}}|\mathbf{x}_{t_{i}})d\mathbb{P}(\mathbf{x}_{0:T})$ from (2). Now, with the choice of reference measure $\mathbb{P}$ induced by Markov process in (1), we can define the conditional SDEs conditioned on $\mathcal{H}_{t_{k}}$ which inducing the desired path measure $\mathbb{P}^{\star}$. Note that this is an extension of the original Doob’s $h$-transform (Doob, ), incorporating multiple marginal constraints. Below, we summarize the relevant result.

###### Theorem 3.1 (Multi-marginal Doob’s $h$ -transform) .

Let us define a sequence of functions $\{h_{i}\}_{i\in[1:k]}$, where each $h_{i}:[t_{i-1},t_{i})\times\mathbb{R}^{d}\to\mathbb{R}_{+}$, for all $i\in[1:k]$, is a conditional expectation $h_{i}(t,\mathbf{x}_{t}):=\mathbb{E}_{\mathbb{P}}\left[\prod_{j\geq i}^{k}f_{j}
(\mathbf{y}_{t_{j}}|\mathbf{X}_{t_{j}})|\mathbf{X}_{t}=\mathbf{x}_{t}\right]$,
where $\{f_{i}\}_{i\in[1:k]}$ is defined in (3). Now, we define a function $h:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}_{+}$ by integrating the functions $\{h_{i}\}_{i\in[1:k]}$,

$$ $h(t,\mathbf{x}):=\sum_{i=1}^{k}h_{i}(t,\mathbf{x})\mathbf{1}_{[t_{i-1},t_{i})} (t).$ (4) $$

Then, with the initial condition $\mu^{\star}_{0}(d\mathbf{x}_{0})=h_{1}(t_{0},\mathbf{x}_{0})\mu_{0}(d\mathbf{x
}_{0})$,
the solution of the following conditional SDE inducing
the posterior path measure $\mathbb{P}^{\star}$ in (2):

$$ $\text{(Conditioned State)}\quad d\mathbf{X}^{\star}_{t}=\left[b(t,\mathbf{X}^{ \star}_{t})+\nabla_{\mathbf{x}}\log h(t,\mathbf{X}^{\star}_{t})\right]dt+d \mathbf{W}_{t}$ (5) $$

Theorem 3.1 demonstrates that we can obtain sample trajectories from $\mathbb{P}^{\star}$ in (2) by simulating the dynamics in (5). However, estimating the functions $\{h_{i}\}_{i\in[1:k]}$ requires both the estimation of the sequence of normalization constants $\{\mathbf{L}_{i}\}_{i\in[1:k]}$ and the computation of conditional expectations, which is infeasible in general. For these reasons, we propose a VI algorithm to approximate the functions $\{h_{i}\}_{i\in[1:k]}$. We derive the variational bound for the VI by exploiting the theory of SOC.

###### Theorem 3.1 (Multi-marginal Doob’s h ℎ h italic_h -transform) .

### 3.2 Stochastic Optimal Control

The SOC (Fleming & Soner, ; Carmona, ) is a mathematical framework that deals with the problem of finding control policies in order to achieve certain object. In this paper, we define following control-affine SDE, adjusting the prior dynamics in (1) with a Markov control $\alpha:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$:

$$ $\textit{(Controlled State)}\quad d\mathbf{X}^{\alpha}_{t}=\left[b(t,\mathbf{X} ^{\alpha}_{t})+\alpha(t,\mathbf{X}^{\alpha}_{t})\right]dt+d\tilde{\mathbf{W}}_ {t},$ (6) $$

where $\mathbf{X}^{\alpha}_{0}\sim\mu_{0}$. We refer to the SDE in (6) as controlled SDE. We can expect that for a well-defined function set $\alpha\in\mathbb{A}$, the class of controlled SDE (6) encompass the SDE in (5). This implies that the desired path measure $\mathbb{P}^{\star}$ can be achieved through the SOC formulation. In general, the goal of SOC is to find the optimal control policy $\alpha^{\star}$ that minimizes a given arbitrary cost function $\mathcal{J}(t,\mathbf{x}_{t},\alpha)$ $\textit{i}.\textit{e}.,\alpha^{\star}(t,\mathbf{x}_{t})=\operatorname*{arg\,
min}_{\alpha\in\mathbb{A}}\mathcal{J}(t,\mathbf{x}_{t},\alpha)$ and determine the value function $\mathcal{V}(t,\mathbf{x}_{t})=\min_{\alpha\in\mathbb{A}}\mathcal{J}(t,\mathbf{
x}_{t},\alpha)$,
where $\mathcal{V}(t,\mathbf{x}_{t}):=\mathcal{J}(t,\mathbf{x}_{t},\alpha^{\star})
\leq\mathcal{J}(t,\mathbf{x}_{t},\alpha)$ holds for any $\alpha\in\mathbb{A}$. Below, we demonstrate how, with a carefully chosen cost function, the theory of SOC establishes a connection between two classes of SDEs (5) and (6). This connection enables the development of a variational inference algorithm with a tight evidence lower bound (ELBO). To this end, we consider the following cost function:

$$ $\mathcal{J}(t,\mathbf{x}_{t},\alpha)=\mathbb{E}_{\mathbb{P}^{\alpha}}\left[ \int_{t}^{T}\frac{1}{2}\left\lVert\alpha(s,\mathbf{X}^{\alpha}_{s})\right \rVert^{2}ds-\sum_{i:\{t\leq t_{i}\}}\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^ {\alpha}_{t_{i}})|\mathbf{X}^{\alpha}_{t}=\mathbf{x}_{t}\right].$ (7) $$

Then, the value function $\mathcal{V}$ for (7) satisfies the dynamic programming principle (Carmona, ):

###### Theorem 3.2 (Dynamic Programming Principle) .

Let us consider a sequence of left continuous functions $\{\mathcal{V}_{i}\}_{i\in[1:k+1]}$, where each $\mathcal{V}_{i}\in C^{1,2}([t_{i-1},t_{i})\times\mathbb{R}^{d})$

$$ $\mathcal{V}_{i}(t,\mathbf{x}_{t}):=\min_{\alpha\in\mathbb{A}}\mathbb{E}_{ \mathbb{P}^{\alpha}}\left[\int_{t_{i-1}}^{t_{i}}\frac{1}{2}\left\lVert\alpha_{ s}\right\rVert^{2}ds-\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}} )+\mathcal{V}_{i+1}(t_{i},\mathbf{X}^{\alpha}_{t_{i}})|\mathbf{X}_{t}=\mathbf{ x}_{t}\right],$ (8) $$

for all $i\in[1:k]$ and $\mathcal{V}_{k+1}=0$. Then, for any $0\leq t\leq u\leq T$, the value function $\mathcal{V}$ for the cost function in (7) satisfying the recursion defined as follows:

$$ $\mathcal{V}(t,\mathbf{x}_{t})=\min_{\alpha\in\mathbb{A}}\mathbb{E}_{\mathbb{P} ^{\alpha}}\left[\int_{t}^{t_{I(u)}}\frac{1}{2}\left\lVert\alpha_{s}\right \rVert^{2}ds-\sum_{i:\{t\leq t_{i}\leq t_{I(u)}\}}\log f_{i}+\mathcal{V}_{I(u) +1}(t_{I(u)},\mathbf{X}^{\alpha}_{t_{I(u)}})|\mathbf{X}^{\alpha}_{t}=\mathbf{x }_{t}\right],$ (9) $$

with the indexing function $I(u)=\max\{i\in[1:k]|t_{i}\leq u\}$ and $f_{i}=f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})$

The value functions presented in Theorem 3.2 suggest that the objective of the optimal control policy $\alpha$ for the interval $[t_{i-1},t_{i})$ is not just to minimize the negative log-potential $-\log f_{i}(\mathbf{y}_{t_{i}}|\cdot)$ for the immediate observation $\mathbf{y}_{t_{i}}$. Instead, it also involves considering future costs $\{\mathcal{V}_{j}\}_{j\in[i+1:k]}$ and the corresponding future observations $\{\mathbf{y}_{t_{j}}\}_{j\in[i+1:k]}$, since $\mathcal{V}_{i}$ follows a recursive structure. Since our goal is to approximate $\mathbb{P}^{\star}$  in (5), it is natural that the optimal control policy $\alpha^{\star}$ should reflect the future observations $\{\mathbf{y}_{t_{j}}\}_{j\in[i+1:k]}$, as the $h$-function inherently does. Next, we will derive the explicit form of the optimal control policy.

###### Theorem 3.3 (Verification Theorem) .

Suppose there exist a sequence of left continuous functions $\mathcal{V}_{i}(t,\mathbf{x})\in C^{1,2}([t_{i-1},t_{i}),\mathbb{R}^{d})$, for all $i\in[1:k]$, satisfying the following Hamiltonian-Jacobi-Bellman (HJB) equation:

$$ $\displaystyle\partial_{t}\mathcal{V}_{i,t}+\mathcal{A}_{t}\mathcal{V}_{i,t}+ \min_{\alpha\in\mathbb{A}}\left[(\nabla_{\mathbf{x}}\mathcal{V}_{i,t})^{\top} \alpha_{i,t}+\frac{1}{2}\left\lVert\alpha_{i,t}\right\rVert^{2}\right]=0,\quad t _{i-1}\leq t<t_{i}$ (10) $\displaystyle\mathcal{V}_{i}(t_{i},\mathbf{x})=-\log f_{i}(\mathbf{y}_{t_{i}}| \mathbf{x})+\mathcal{V}_{i+1}(t_{i},\mathbf{x}),\quad t=t_{i},\quad\forall i \in[1:k],$ (11) $$

where a minimum is attained by $\alpha^{\star}_{i}(t,\mathbf{x})=\nabla_{\mathbf{x}}\mathcal{V}_{i}(t,\mathbf{
x})$. Now, define a function $\alpha:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ by integrating the optimal controls $\{\alpha_{i}\}_{i\in[1:k]}$,

$$ $\alpha^{\star}(t,\mathbf{x}):=\sum_{i=1}^{k}\alpha^{\star}_{i}(t,\mathbf{x}) \mathbf{1}_{[t_{i-1},t_{i})}(t)$ (12) $$

Then, $\mathcal{V}(t,\mathbf{x}_{t})=\mathcal{J}(t,\mathbf{x}_{t},\alpha^{\star})\leq
\mathcal{J}(t,\mathbf{x}_{t},\alpha)$ holds for any $(t,\mathbf{x}_{t})\in[0,T]\times\mathbb{R}^{d}$ and $\alpha\in\mathbb{A}$.

Note that the optimal control (12) for the cost function (7) share a similar structure as in (4). Since the theory of PDEs establishes a fundamental link between various classes of PDEs and SDEs (Richter & Berner, ; Berner et al., ), it allows us to reveal the inherent connection between (12) and (4), thereby enable us to simulate the conditional SDE (5) in an alternative way.

###### Lemma 3.4 (Hopf-Cole Transformation) .

The $h$ function satisfying the following linear PDE:

$$ $\displaystyle\partial_{t}h_{i,t}+\mathcal{A}_{t}h_{i,t}=0,\quad t_{i-1}\leq t< t_{i}$ (13) $\displaystyle h_{i}(t_{i},\mathbf{x})=f_{i}(\mathbf{y}_{t_{i}}|\mathbf{x})h_{i +1}(t_{i},\mathbf{x}),\quad t=t_{i},\quad\forall i\in[1:k].$ (14) $$

Moreover, for a logarithm transformation $\mathcal{V}=-\log h$, $\mathcal{V}$ satisfying the HJB equation in (10-11).

According to Lemma 3.4, the solution of linear PDE in (13-14) is negative exponential to the solution of the HJB equation in (10-11). Therefore, it leads to the following corollary:

###### Corollary 3.5 (Optimal Control) .

For optimal control $\alpha^{\star}$ induced by the cost function (7) with dynamics (6), it satisfies $\alpha^{\star}=\nabla_{\mathbf{x}}\log h$. In other words, we can simulate the conditional SDEs in (5) by simulating the controlled SDE (6) with optimal control $\alpha^{\star}$.

Corollary 3.5 states that the Markov process induced by $\alpha^{\star}$ and $\nabla_{\mathbf{x}}\log h$ is equivalent under same initial condition. However, comparing the $\mathbb{P}^{\alpha}$ induced by the controlled dynamics in (6) with an initial condition $\mu_{0}$, the conditioned dynamics (5) has an intial condition $\mu_{0}^{\star}$ to inducing the desired path measure $\mathbb{P}^{\star}$. In other words, although we find the optimal control $\alpha^{\star}$, the constant discrepancy between $\mu_{0}$ and $\mu^{\star}_{0}$ still remain, thereby keeping $\mathbb{P}^{\star}$ and $\mathbb{P}^{\alpha^{\star}}$ misaligned. Fortunately, a surrogate cost function can be derived from the variational representation under the KL-divergence, allowing us to find the optimal control while minimizing the discrepancy.

###### Theorem 3.6 (Tight Variational Bound) .

Let us assume that the path measure $\mathbb{P}^{\alpha}$ induced by (6) for any $\alpha\in\mathbb{A}$ satisfies $D_{\text{KL}}(\mathbb{P}^{\alpha}|\mathbb{P}^{\star})<\infty$. Then, for a cost function $\mathcal{J}$ in (7) and $\mu^{\star}_{0}$ in (5), it holds:

$$ $\displaystyle D_{\text{KL}}(\mathbb{P}^{\alpha}|\mathbb{P}^{\star})=D_{\text{ KL}}(\mu_{0}|\mu^{\star}_{0})+\mathbb{E}_{\mathbf{x}_{0}\sim\mu_{0}}\left[ \mathcal{J}(0,\mathbf{x}_{0},\alpha)\right]=\mathcal{L}(\alpha)+\log\mathbf{Z} (\mathcal{H}_{t_{k}})\geq 0,$ (15) $$

where the objective function $\mathcal{L}(\alpha)$ (a negative ELBO) is given by:

$$ $\mathcal{L}(\alpha)=\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\int_{0}^{T}\frac{1} {2}\left\lVert\alpha(s,\mathbf{X}^{\alpha}_{s})\right\rVert^{2}ds-\sum_{i}^{k} \log g_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})\right]\geq-\log \mathbf{Z}(\mathcal{H}_{t_{k}}).$ (16) $$

Moreover, assume that $\mathcal{L}(\alpha)$ has a global minimum $\alpha^{\star}=\operatorname*{arg\,min}_{\alpha\in\mathbb{A}}\mathcal{L}(\alpha)$. Then the equality holds in (16) $\textit{i}.\textit{e}.,\mathcal{L}(\alpha^{\star})=-\log\mathbf{Z}(\mathcal{H}
_{t_{k}})$ and $\mu^{\star}_{0}=\mu_{0}$ almost everywhere with respect to $\mu_{0}$.

Theorem 3.6 suggests that the optimal control $\alpha^{\star}=\operatorname*{arg\,min}_{\alpha}\mathcal{L}(\alpha)$ provides the tight variational bound for the likelihood functions $\{g_{i}\}_{i\in[1:k]}$ and the prior path measure $\mathbb{P}$ induced by (1). Furthermore, $\alpha^{\star}$ ensures that $\mu_{0}^{\star}=\mu_{0}$, indicating that simulating the controlled SDE in (6) with $\alpha^{\star}$ and initial condition $\mu_{0}$ generates trajectories from the posterior path measure $\mathbb{P}^{\star}$ in (2).

In practice, the optimal control $\alpha^{\star}$ can be approximated using a highly flexible neural network, which serves as a function approximator ($\textit{i}.\textit{e}.,\alpha:=\alpha(\cdot;\theta)$), optimized through gradient descent-based optimization (Li et al., ; Zhang & Chen, ; Vargas et al., ).

However, applying gradient descent-based optimization necessitates computing gradients through the simulated diffusion process over the interval $[0,T]$ to estimate the objective function (16) for neural network training. This approach can become slow, unstable, and memory-intensive as the time horizon or dimension of latent space increases (Iakovlev et al., ; Park et al., ). It contrasts with the philosophy of many recent generative models (Ho et al., ; Song et al., ), which aim to decompose the generative process and solve the sub-problems jointly. Additionally, for inference, it requires numerical simulations such as Euler-Maruyama solvers (Kloeden & Platen, ), which can also be time-consuming for a long time series. It motivated us to propose an efficient and scalable modeling approach for real-world applications described in the next section.

Figure: Figure 1: Conceptual illustration. Given the observed time stamps $\mathcal{T}=\{t_{i}\}_{i\in[1:4]}$ and the unseen time stamps ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathcal{T}_ {u}}$($\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\times}$ in figure), the encoder maps the input time series $\{\mathbf{o}_{t}\}_{t\in\mathcal{T}}$ into auxiliary variables $\{\mathbf{y}_{t}\}_{t\in\mathcal{T}}$. These auxiliary variables are then utilized to compute the control policies $\{\alpha\}_{i\in[1:5]}$ through a masked attention mechanism that relies on two different assimilation schemes. The computed policies $\{\alpha\}_{i\in[1:5]}$ control the prior dynamics $\mathbb{P}$ over the interval $[0,T]$ to approximate the posterior $\mathbb{P}^{\star}$ in the latent space. Finally, the sample path ${\color[rgb]{0.1796875,0.45703125,0.71484375}\definecolor[named]{ pgfstrokecolor}{rgb}{0.1796875,0.45703125,0.71484375}\mathbf{X}^{\alpha}_{0:T} }\sim\mathbb{P}^{\alpha}$ are decoded to generate predictions across the complete time stamps ${\color[rgb]{0.484375,0.625,0.3984375}\definecolor[named]{pgfstrokecolor}{rgb} {0.484375,0.625,0.3984375}\mathcal{T}^{\prime}}=\mathcal{T}\cup\mathcal{T}_{u}$(${\color[rgb]{0.484375,0.625,0.3984375}\definecolor[named]{pgfstrokecolor}{rgb} {0.484375,0.625,0.3984375}\circ}$ in figure), over the entire interval $[0,T]$.
Refer to caption: x1.png

###### Theorem 3.2 (Dynamic Programming Principle) .

###### Theorem 3.3 (Verification Theorem) .

###### Lemma 3.4 (Hopf-Cole Transformation) .

###### Corollary 3.5 (Optimal Control) .

###### Theorem 3.6 (Tight Variational Bound) .

### 3.3 Efficient Modeling of the Latent System

The linear approximation of the general drift function provides a simulation-free property for the dynamics and significantly enhances scalability while ensuring high-quality performance (Deng et al., ). Motivate by this property, we investigate a class of linear SDEs to improve the efficiency of the proposed controlled dynamics, ensuring superior performance compared to other baselines. We introduce the following affine linear SDEs:

$$ $d\mathbf{X}_{t}=\left[-\mathbf{A}\mathbf{X}_{t}+\alpha\right]dt+d\mathbf{W}_{t },\quad\text{where}\quad\mu_{0}\sim\mathcal{N}(\mathbf{m}_{0},\mathbf{\Sigma}_ {0}),$ (17) $$

and a matrix $\mathbf{A}\in\mathbb{R}^{d\times d}$ and a vector $\alpha\in\mathbb{R}^{d}$. The solutions for $\mathbf{X}_{t}$ in (17) has a closed-form Gaussian distribution for any $t\in[0,T]$, where the mean $\mathbf{m}_{t}$ and covariance $\mathbf{\Sigma}_{t}$ can be explicitly computed by solving the ODEs (Särkkä & Solin, ):

$$ $\displaystyle\mathbf{m}_{t}=e^{-\mathbf{A}t}\mathbf{m}_{0}-\mathbf{A}^{-1}(e^{ -\mathbf{A}t}-\mathbf{I})\alpha$ (18) $\displaystyle\mathbf{\Sigma}_{t}=e^{-\mathbf{A}t}\mathbf{\Sigma}_{0}e^{- \mathbf{A}^{\top}t}+\int_{0}^{t}e^{-\mathbf{A}(t-s)}e^{\mathbf{A}^{\top}(t-s)}ds.$ (19) $$

However, calculating the moments $\mathbf{m}_{t}$ and $\mathbf{\Sigma}_{t}$ in (18) involves computing matrix exponentials, inversions, and performing numerical integration. These operations can be computationally intensive, especially for large matrices or when high precision is required. These computations can be simplifies by restricting the matrix $\mathbf{A}$ to be a diagonal or semi-positive definite (SPD).

###### Remark 3.7 (Diagonalization) .

Since SPD matrix $\mathbf{A}$ admits the eigen-decomposition $\mathbf{A}=\mathbf{E}\mathbf{D}\mathbf{E}^{\top}$ with $\mathbf{E}\in\mathbb{R}^{d\times d}$ and $\mathbf{D}\in\text{diag}(\mathbb{R}^{d})$, the process $\mathbf{X}_{t}$ expressed in a standard basis can be transformed to a $\hat{\mathbf{X}}_{t}$ which have diagonalized drift function. In the space spanned by the eigen-basis $\mathbf{E}$, the dynamics in (1) can be rewritten into:

$$ $d\hat{\mathbf{X}}_{t}=\left[-\mathbf{D}\hat{\mathbf{X}}_{t}+\hat{\alpha}\right ]dt+d\hat{\mathbf{W}}_{t},\quad\text{where}\quad\hat{\mu}_{0}\sim\mathcal{N}( \hat{\mathbf{m}_{0}},\hat{\mathbf{\Sigma}_{0}}),$ (20) $$

$\hat{\mathbf{X}}_{t}=\mathbf{E}^{\top}\mathbf{X}_{t}$, $\hat{\alpha}=\mathbf{E}^{\top}\alpha$, $\hat{\mathbf{W}}_{t}=\mathbf{E}^{\top}\mathbf{W}_{t}$, $\hat{\mathbf{m}}_{t}=\mathbf{E}^{\top}\mathbf{m}_{t}$ and $\hat{\mathbf{\Sigma}}_{t}=\mathbf{E}^{\top}\mathbf{\Sigma}_{t}\mathbf{E}$. Note that $\hat{\mathbf{W}}_{t}\stackrel{{\scriptstyle d}}{{=}}\mathbf{E}^{\top}\mathbf{W
}_{t}$ for any $t\in[0,T]$ due to the orthonormality of $\mathbf{E}$, so $\hat{\mathbf{W}}_{t}$ can be regarded as a standard Wiener process. Because of $\mathbf{D}=\text{diag}(\bm{\lambda})$, where $\bm{\lambda}=\{\lambda_{1},\cdots,\lambda_{d}\}$ and each $\lambda_{i}\geq 0$ for all $i\in[1:d]$, we can obtain the state distributions of $\mathbf{X}_{t}$ for any $t\in[0,T]$ by solving ODEs in (18-19) analytically, without the need for numerical computation. The results are then transforming back to the standard basis $\textit{i}.\textit{e}.,\mathbf{m}_{t}=\mathbf{E}\hat{\mathbf{m}}_{t},\mathbf{
\Sigma}_{t}=\mathbf{E}\hat{\mathbf{\Sigma}}_{t}\mathbf{E}^{T}$.

###### Remark 3.7 (Diagonalization) .

#### Locally Linear Approximation

To leverage the advantages of linear SDEs in (17) which offer simulation-free property, we aim to linearize the drift function in (6). However, the naïve formulation described in (17) may limit the expressiveness of the latent dynamics for real-world applications. Hence, we introduce a parameterization strategy inspired by (Becker et al., ; Klushyn et al., ), which leverage neural networks to enhance the flexibility by fully incorporating an attentive structure with a given observations $\mathbf{y}_{0:T}$ while maintaining a linear formulation:

$$ $\displaystyle d\mathbf{X}^{\alpha}_{t}=\left[-\mathbf{A}_{t}\mathbf{X}_{t}+ \alpha_{t}\right]dt+d\tilde{\mathbf{W}}_{t},$ (21) $$

where the approximated drift function is constructed as affine formulation with following components:

$$ $\mathbf{A}_{t}=\sum_{l=1}^{L}w^{(l)}_{\theta}(\mathbf{z}_{t})\mathbf{A}^{(l)}, \quad\alpha_{t}=\mathbf{B}_{\theta}\mathbf{z}_{t}.$ (22) $$

Figure: Figure 2: Two types of information assimilation schemes.
Refer to caption: x2.png

The matrix $\mathbf{A}_{t}$ is given by a convex combination of $L$ trainable base matrices $\{\mathbf{A}^{(l)}\}_{l\in[1:L]}$, where the weights $\mathbf{w}_{\theta}=\text{softmax}(\mathbf{f}_{\theta}(\mathbf{z}_{t}))$ are produced by the neural network $\mathbf{f}_{\theta}$. Additionally, $\mathbf{B}_{\theta}\in\mathbb{R}^{d\times d}$ is a trainable matrix. The latent variable $\mathbf{z}_{t}$ is produced by the transformer $\mathbf{T}_{\theta}$, which encodes the given observations $\mathbf{y}_{0:T}$, depending on the task-specific information assimilation scheme.

Figure 2 illustrates two assimilation schemes using marked attention mechanism(^3^33See Appendix C.2 for more details.): the history assimilation scheme, the transformer $\mathbf{T}_{\theta}$ encodes information up to the current time $t$ and outputs $\mathbf{z}_{t}$, $\textit{i}.\textit{e}.,\mathbf{z}_{t}=\mathbf{T}_{\theta}(\mathcal{H}_{t})$, and the full assimilation scheme, the transformer $\mathbf{T}_{\theta}$ encodes information over the entire interval $[0,T]$ and outputs $\mathbf{z}_{t}$, $\textit{i}.\textit{e}.,\mathbf{z}_{t}=\mathbf{T}_{\theta}(\mathcal{H}_{T})$. Note that this general formulation brought from the control formulation enables more flexible use of information encoded from observations. In contrast, previous Kalman-filtering based CD-SSM method (Schirmer et al., ) relies on recurrent updates, which limits them to typically using historical information only.

Furthermore, since observations are updated only at discrete time steps $\{t_{i}\}_{i\in[1:k]}$, the latent variables $\mathbf{z}_{t}$ remain constant within any interval $t\in[t_{i-1},t_{i})$ for all $i\in[1:k]$, making $\mathbf{A}_{i}$ and $\alpha_{i}$ constant as well. As a result, the dynamics (21) remain linear over local intervals. This structure enables us to derive a closed-form solution for the intermediate latent states.

###### Theorem 3.8 (Simulation-free estimation) .

Let us consider sequences of SPD matrices $\{\mathbf{A}_{i}\}_{i\in[1:k]}$ and vectors $\{\alpha_{i}\}_{i\in[1:k]}$ and the following control-affine SDEs for all $i\in[1:k]$:

$$ $d\mathbf{X}_{t}=\left[-\mathbf{A}_{i}\mathbf{X}_{t}+\alpha_{i}\right]dt+\sigma d \mathbf{W}_{t},\quad t\in[t_{i-1},t_{i}).$ (23) $$

Then, the solution of (23) condition to the initial distribution $\mu_{0}=\mathcal{N}(\mathbf{m}_{0},\mathbf{\Sigma}_{0})$ is a Gaussian process $\mathcal{N}(\mathbf{m}_{t},\mathbf{\Sigma}_{t})$ for all $t\in[0,T]$, where the first two moments is given by

$$ $\displaystyle\mathbf{m}_{t}=e^{-(t-t_{0})\mathbf{A}_{i}}\mathbf{m}_{0}-\sum_{j =1}^{i-1}\left(e^{-(t_{i-1}-t_{j})\mathbf{A}_{j}}\mathbf{A}_{j}^{-1}(e^{-(t_{j }-t_{j-1})\mathbf{A}_{j}}-\mathbf{I})\alpha_{j}\right)-\mathbf{A}_{i}^{-1}(e^{ -(t-t_{i-1})\mathbf{A}_{i}}-\mathbf{I})\alpha_{i},$ $\displaystyle\mathbf{\Sigma}_{t}=e^{-2(t-t_{0})\mathbf{A}_{i}}\mathbf{\Sigma}_ {0}-\sum_{j=1}^{i-1}\left(e^{-2(t_{i-1}-t_{j})\mathbf{A}_{j}}\mathbf{A}_{j}^{- 1}\frac{(e^{-2(t_{j}-t_{j-1})\mathbf{A}_{j}}-\mathbf{I})}{2}\right)-\mathbf{A} _{i}^{-1}\frac{(e^{-2(t-t_{i-1})\mathbf{A}_{i}}-\mathbf{I})}{2}.$ $$

###### Theorem 3.8 (Simulation-free estimation) .

#### Parallel Computation

Given an associative operator $\otimes$ and a sequence of elements $[s_{t_{1}},\cdots s_{t_{K}}]$, the parallel scan algorithm computes the all-prefix-sum which returns the sequence

$$ $[s_{t_{1}},(s_{t_{1}}\otimes s_{t_{2}}),\cdots,(s_{t_{1}}\otimes s_{t_{2}} \otimes\cdots\otimes s_{t_{K}})]$ (24) $$

in $\mathcal{O}(\log K)$ time. Leveraging the linear formulation described in Theorem 3.8 and the inherent parallel nature of the transformer architecture for sequential structure, our method can be integrated with the parallel scan algorithm (Blelloch, ) resulting efficient computation of the marginal Gaussian distributions by computing both moments $\{\mathbf{m}_{t}\}_{t\in[0,T]}$ and $\{\mathbf{\Sigma}_{t}\}_{t\in[0,T]}$ in a parallel sense(^4^44See Appendix B for more details.).

###### Remark 3.9 (Non-Markov Control) .

Note that (21-22) involves approximating the Markov control by a non-Markov control $\alpha(\mathcal{H}_{t}):=\alpha^{\theta}$, parameterized by neural network $\theta$. However, Theorem 3.3 establishes that the optimal control should be Markov, as it is verified by the HJB equation (Van Handel, ). In our case, we expect that with a high-capacity neural network $\theta$, the local minimum $\theta^{M}$, obtained after $M$ gradient descent steps $\theta^{m+1}=\theta^{m}-\nabla_{\theta}\mathcal{L}(\alpha^{\theta^{m}})$ yields $\mathcal{L}(\alpha^{\star})\approx\mathcal{L}(\alpha^{\theta^{m}\to\theta^{M}})$.

###### Remark 3.9 (Non-Markov Control) .

#### Amortization

Moreover, to enhance flexibility, we treat $\mathbf{y}_{0:T}$ as an auxiliary variable in the latent space which is produced by a neural network encoder $q_{\phi}$ applied to the given time series data $\mathbf{o}_{0:T}$ $\textit{i}.\textit{e}.,\mathbf{y}_{0:T}\sim q_{\phi}(\mathbf{y}_{0:T}|\mathbf{
o}_{0:T})$, where it is factorized as

$$ $q_{\phi}(\mathbf{y}_{0:T}|\mathbf{o}_{0:T})=\prod_{i=1}^{k}q_{\phi}(\mathbf{y} _{t_{i}}|\mathbf{o}_{t_{i}})=\prod_{i=1}^{k}\mathcal{N}(\mathbf{y}_{t_{i}}| \mathbf{q}_{\phi}(\mathbf{o}_{t_{i}}),\mathbf{\Sigma}_{q})$ (25) $$

with a fixed variance $\mathbf{\Sigma}_{q}$. Additionally, it enables the modeling of nonlinear emission distributions through a neural network decoder $p_{\psi}(\mathbf{o}_{0:T}|\mathbf{y}_{0:T})$, where it is factorized as

$$ $p_{\psi}(\mathbf{o}_{0:T}\mid\mathbf{y}_{0:T})=\prod_{i=1}^{k}p_{\psi}(\mathbf {o}_{t_{i}}\mid\mathbf{y}_{t_{i}}),$ (26) $$

with the likelihood function $p_{\psi}$ depending on the task at hand. This formulation separates the representation learning of $\mathbf{y}_{0:T}$ from the dynamics of $\mathbf{X}^{\alpha}_{0:T}$, allowing for a more flexible parameterization. Typically, the information capturing the underlying physical dynamics resides in a much lower-dimensional space compared to the original sequence (Fraccaro et al., ). Therefore, separating generative modeling through lower-dimensional latent dynamics from representation modeling in the original space (e.g., pixel values in an image sequence) offers greater flexibility in parameterization.

### 3.4 Training and Inference

#### Training Objective Function

We jointly train, in an end-to-end manner, the amortization parameters $\{\phi,\psi\}$ for the encoder-decoder pair, along with the parameters of the latent dynamics $\theta=\{\mathbf{f}_{\theta},\mathbf{B}_{\theta},\mathbf{T}_{\theta},\mathbf{m
}_{0},\mathbf{\Sigma}_{0},\{\mathbf{A}^{(l)}\}_{l\in[1:L]}\}$, which include the parameters required for controlled latent dynamics. The training is achieved by maximizing the evidence lower bound (ELBO) of the observation log-likelihood for a given the time series $\mathbf{o}_{0:T}$:

$$ $\displaystyle\log p_{\psi}(\mathbf{o}_{0:T})$ $\displaystyle\geq\mathbb{E}_{\mathcal{H}_{T}\sim q_{\phi}(\mathbf{y}_{0:T}| \mathbf{o}_{0:T})}\left[\log\frac{\prod_{i=0}^{K}p_{\psi}(\mathbf{o}_{t_{i}}| \mathbf{y}_{t_{i}})g(\mathbf{y}_{0:T})}{\prod_{i=0}^{K}q_{\phi}(\mathbf{y}_{t_ {i}}|\mathbf{o}_{t_{i}})}\right]$ (27) $\displaystyle\gtrapprox\mathbb{E}_{\mathcal{H}_{T}\sim q_{\phi}(\mathbf{y}_{0: T}|\mathbf{o}_{0:T})}\sum_{i=1}^{K}\left[\log p_{\psi}(\mathbf{o}_{t_{i}}| \mathbf{y}_{t_{i}})-\mathcal{L}(\theta)\right]=\text{ELBO}(\psi,\phi,\theta)$ (28) $$

Since $\mathbf{Z}(\mathcal{H}_{t_{k}})=g(\mathbf{y}_{0:T})$, the prior over auxiliary variable $g(\mathbf{y}_{0:T})$ can be computed using the ELBO $\mathcal{L}(\theta)$ proposed in (7) as part of our variational inference procedure for the latent posterior $\mathbb{P}^{\star}$ in proposed Sec 2, including all latent parameters $\theta$ containing the control $\alpha^{\theta}$. Note that our modeling is computationally favorable for both training and inference, since estimating marginal distributions in a latent space can be parallelized. This allows for the generation of a latent trajectory over the entire interval without the need for numerical simulations. The overall training and inference processes are summarized in the Algorithm 1 and Algorithm 2 in Appendix, respectively.

## 4 Experiment

In this section, we present empirical results demonstrating the effectiveness of ACSSM in modeling real-world irregular time-series data. The primary objective was to evaluate its capability to capture the underlying dynamics across various datasets. To demonstrate the applicability of the ACSSM, we conducted experiments on four tasks: per-time regression/classification and sequence interpolation/extrapolation, using four datasets: Human Activity, USHCN (Menne et al., ), and Physionet (Silva et al., ). We compare our approach against various baselines including RNN architecture (RKN-$\Delta_{t}$ (Becker et al., ), GRU-$\Delta_{t}$ (Chung et al., ), GRU-D (Che et al., )) as well as dynamics-based models (Latent-ODE (Chen et al., ; Rubanova et al., ), ODE-RNN (Rubanova et al., ), GRU-ODE-B (De Brouwer et al., ), CRU (Schirmer et al., ), Latent-SDE${}_{\text{H}}$ (Zeng et al., )) and attention-based models (mTAND (Shukla & Marlin, )), which have been developed for modeling irregular time series data. We reported the averaged results over five runs with different seed. The best results are highlighted in bold, while the second-best results are shown in blue. Additional experimental details can be found in Appendix C.

### 4.1 Per time point classification & \& & regression

**Table 1: Test Accuracy ($\%$).**
| Model | Acc |
| --- | --- |
| Latent-ODE^† | $87.0\pm 2.8$ |
| Latent-SDE${}_{\text{H}}$^‡ | $90.6\pm 0.4$ |
| mTAND^† | $91.1\pm 0.2$ |
| ACSSM (Ours) | $\mathbf{91.4\pm 0.4}$ |

#### Human Activity Classification

We investigated the classification performance of our proposed model. For this purpose, we trained the model on the Human Activity dataset, which contains time-series data from five individuals performing various activities such as walking, sitting, lying, and standing etc. The dataset includes 12 features in total, representing 3D positions captured by four sensors attached to the belt, chest, and both ankles. Following the pre-processing approach proposed by (Rubanova et al., ), the dataset comprises 6,554 sequences, each with 211 time points. The task is to classify each time point into one of seven activities.

Table 1 reports the test accuracy, showing that ACSSM outperforms all baseline models. We employed the full assimilation scheme to maintain consistency with the other baselines, which infer the latent state using full observation. It is important to note that the two dynamical models, Latent-ODE and Latent-SDE${}_{\text{H}}$, incorporate a parameterized vector field in their latent dynamics, thereby rely on numerical solvers to infer intermediate states. Besides, mTAND is an attention-based method that does not depend on dynamical or state-space models, thus avoiding numerical simulation. We believe the significant performance improvement of our model comes from its integration of an attention mechanism into dynamical models. Simulation-free dynamics avoid numerical simulations while maintaining temporal structure, leading to more stable learning.

**Table 2: Test MSE ($\times 10^{-3}$).**
| Model | MSE |
| --- | --- |
| Latent-ODE^† | $15.70\pm 0.29$ |
| CRU^† | $4.63\pm 1.07$ |
| Latent-SDE${}_{\text{H}}$^‡ | $3.84\pm 0.35$ |
| S5^⋆ | $3.41\pm 0.27$ |
| mTAND^‡ | $3.20\pm 0.60$ |
| ACSSM (Ours) | $\mathbf{2.98\pm 0.30}$ |

#### Pendulum Regression

Next, We explored the problem of sequence regression using pendulum experiment (Becker et al., ), where the goal is to infer the sine and cosine of the pendulum angle from irregularly observed, noisy pendulum images (Schirmer et al., ). To assess our performance, we compared it with previous dynamics-based models, reporting the regression MSE on a held-out test set as shown in Table 2. We employed the full assimilation scheme. The experimental results demonstrated that our proposed method outperformed existing models, delivering superior performance. These findings highlight that, even when linearizing the drift function, the amortization and proposed neural network based locally linear dynamics in Sec 3.3 preserves the expressivity of our approach, enabling more accurate inference of non-linear systems.

Moreover, particularly in comparison to CRU, we believe that the significant performance improvement stems from the fundamental differences in how information is leveraged. To infer intermediate angular values, utilizing not only past information but also future positions of the pendulum can enhance the accuracy of these predictions. In this regard, while CRU relies solely on the past positions of the pendulum for its predictions, our model is capable of utilizing both past and future positions.

### 4.2 Sequence Interpolation & \& & Extrapolation

#### Datasets

We benchmark the models on two real-world datasets, USHCN and Physionet. The USHCN dataset (Menne et al., ) contains 1,218 daily measurements from weather stations across the US with five variables over a four-year period. The Physionet dataset (Silva et al., ) contains 8,000 multivariate clinical time-series of 41 features recorded over 48 hours.(^5^55See Appendix C.1 for details on the datasets.)

#### Interpolation

We begin by evaluating the effectiveness of ACSSM on the interpolation task. Following the approach of Schirmer et al. () and Rubanova et al. (), each model is required to infer all time points $t\in\mathcal{T^{\prime}}$ based on a subset of observations $\mathbf{o}_{t\in\mathcal{T}}$ where $\mathcal{T}\subseteq\mathcal{T}^{\prime}$. For the interpolation task, the encoded observations $\mathbf{y}_{t\in\mathcal{T}}$ were assimilated by using the full assimilation scheme for the construction of the accurate smoothing distribution. The interpolation results presented in Table 3, where we report the test MSE evaluated on the entire time points $\mathcal{T}^{\prime}$. For all datasets, ACSSM outperforms other baselines in terms of test MSE. It clearly indicate the expressiveness of the ACSSM, trajectories $\mathbf{X}^{\alpha}_{t\in\mathcal{T}^{\prime}}$ sampled from approximated path measure over the entire interval $\mathcal{T^{\prime}}$ are contain sufficient information for generating accurate predictions.

#### Extrapolation

We evaluated ACSSM’s performance on the extrapolation task following the experimental setup of Schirmer et al. (). Each model infer values for all time stamps $t\in\mathcal{T}^{\prime}$, where $\mathcal{T}^{\prime}$ denotes the union of observed time stamps $\mathcal{T}=\{t_{i}\}_{i\in[1:k]}$ and unseen time stamps $\mathcal{T}_{u}=\{t_{i}\}_{i\in[k+1:N]}$, $\textit{i}.\textit{e}.,$ $\mathcal{T}^{\prime}=\mathcal{T}\cup\mathcal{T}_{u}$. For the Physionet dataset, input time stamps $\mathcal{T}$ covered the first 24 hours, while target time stamps $\mathcal{T}^{\prime}$ spanned the rest hours. In the USHCN dataset, the timeline was split evenly, with $t_{k}=\frac{N}{2}$. We report the test MSE for unseen time stamps $\mathcal{T}_{u}=\mathcal{T}^{\prime}-\mathcal{T}$ based on the observations on time stamps $\mathcal{T}$. For modeling an accurate filtering distribution, we employed the history assimilation scheme for this task. As illustrated in Table 3, ACSSM consistently outperformed all baseline models in terms of MSE on the USHCN dataset, achieving a significant performance gain over the second-best model. For the Physionet dataset, ACSSM exhibited comparable performance.

#### Computational Efficiency

To evaluate the training costs in comparison to dynamics-based models that depend on numerical simulations, we re-ran the CRU model on the same hardware used for training our model, indicated by ^∗ in Table 3. Specifically, we utilized a single NVIDIA RTX A6000 GPU. As illustrated in Table 3, ACSSM significantly lowers training costs compared to dynamics-based models. Notable, ACSSM demonstrated a runtime that was more than 16$-$25$\times$ faster than CRU, even though both models aim to approximate $\mathbb{P}^{\star}$ as well. It highlight that the latent modelling approach discussed in Sec 3.3 improves efficiency while achieving an accurate approximation of $\mathbb{P}^{\star}$.

**Table 3: Test MSE ($\times 10^{-2}$) for inter/extra-polation on USHCN and Physionet.**
| Model | Interpolation | Extrapolation | Runtime (sec./epoch) |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| USHCN | Physionet | USHCN | Physionet | USHCN | Physionet |  |
| mTAND^† | 1.766 ± 0.009 | 0.208 ± 0.025 | 2.360 ± 0.038 | 0.340 ± 0.020 | 7 | 10 |
| RKN-$\Delta_{t}^{{\dagger}}$ | 0.009 ± 0.002 | 0.186 ± 0.030 | 1.491 ± 0.272 | 0.703 ± 0.050 | 94 | 39 |
| GRU-$\Delta_{t}^{{\dagger}}$ | 0.090 ± 0.059 | 0.271 ± 0.057 | 2.081 ± 0.054 | 0.870 ± 0.077 | 3 | 5 |
| GRU-D^† | 0.944 ± 0.011 | 0.338 ± 0.027 | 1.718 ± 0.015 | 0.873 ± 0.071 | 292 | 5736 |
| Latent-ODE^† | 1.798 ± 0.009 | 0.212 ± 0.027 | 2.034 ± 0.005 | 0.725 ± 0.072 | 110 | 791 |
| ODE-RNN^† | 0.831 ± 0.008 | 0.236 ± 0.009 | 1.955 ± 0.466 | 0.467 ± 0.006 | 81 | 299 |
| GRU-ODE-B^† | 0.841 ± 0.142 | 0.521 ± 0.038 | 5.437 ± 1.020 | 0.798 ± 0.071 | 389 | 90 |
| CRU^† | 0.016 ± 0.006 | 0.182 ± 0.091 | 1.273 ± 0.066 | 0.629 ± 0.093 | 122 (57.8)^∗ | 114 (63.5)^∗ |
| ACSSM (Ours) | 0.006 ± 0.001 | 0.116 ± 0.011 | 0.941± 0.014 | 0.627 ± 0.019 | 2.3 | 3.8 |

## 5 Conclusion and limitation

In this work, we proposed the method for modeling time series with irregular and discrete observations, which we called ACSSM. By using a multi-marginal Doob’s $h$-transform and a variational inference algorithm by exploiting the theory of SOC, ACSSM efficiently simulates conditioned dynamics. It leverages amortized inference, a simulation-free latent dynamics framework, and a transformer-based data assimilation scheme for scalable and parallel inference. Empirical results show that ACSSM outperforms in various tasks such as classification, regression, and extrapolation while maintaining computational efficiency across real-world datasets.

Although we present the theoretical basis of our method, a thorough analysis in followings remain an open challenge. The variational gap due to the linear approximation may lead to cumulative errors over time, which requires further examination specially for a long time-series such as LLM. Unlike SMC variants (Heng et al., ; Lu & Wang, ) that use particle-based importance weighting to mitigate approximation errors, ACSSM depends on high-capacity neural networks to accurately approximate the optimal control. Incorporating multiple controls, akin to a multi-agent dynamics approach (Han & Hu, ), may alleviate these challenges by enhancing flexibility and robustness.

## Reproducibility Statement

On the theoretical part, all proofs and assumptions are left to Appendix A due the space constraint.
The training and inference algorithms are detailed in Algorithm 1 and Algorithm 2, respectively. Additional implementation details, such as data preprocessing are included in Appendix C. We believe these details are sufficient for interested readers to reproduce the results.

## Ethics Statement

In this work, we proposed a method for modeling irregular time series for practical applications, suggesting that ACSSM does not directly influence ethical or societal issues in a positive or negative way. However, because ACSSM can be applied to health-care datasets, we believe it has the potential to benefit society by improving health and well-being of people.

## Appendix A Proofs and Derivations

In this section, we present the proofs and derivations for all relevant theorems, lemmas, and corollaries We first restate the core concepts of stochastic calculus, which will be used without further explanation.

Assumptions. Throughout the paper, we work with a probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\in[0,T]},\mathbb{P})$, where the filtration $\{\mathcal{F}_{t}\}_{t\in[0,T]}$ supports an $\mathbb{R}^{d}$-valued $\mathcal{F}_{t}$-adapted Wiener process $\mathbf{W}_{t}$ for all $t\in[0,T]$. It is important to note that the $\mathbb{P}$-null set is included in $\mathcal{F}_{0}$, indicating that any event with probability zero at time 0 0 is measurable in the initial $\sigma$-algebra.

We assume that $b$ and $\alpha$ satisfy following conditions:

- •
(Lipschitz condition): For any $t\in[0,T]$, $w\in\Omega$, and $\mathbf{x},\mathbf{x}^{\prime}\in\mathbb{R}^{d}$, where $c_{0}>0$ is a Lipschitz constant, the functions satisfy the inequality $|b(t,w,\mathbf{x})-b(t,w,\mathbf{x}^{\prime})|\leq c_{0}|\mathbf{x}-\mathbf{x}
^{\prime}|$.
- •
(Linear growth condition): For every $\mathbf{x}\in\mathbb{R}^{d}$, the $\mathcal{F}_{t}$-progressively measurable processes $b(t,\mathbf{x})_{t\in[0,T]}$ satisfy $\mathbb{E}\left[\int_{0}^{T}|b_{s}|^{2}ds\right]<\infty$ and $|b(t,\mathbf{x})|\leq c_{1}(1+|\mathbf{x}|)$ for $t\in[0,T]$ and $c_{1}>0$.
- •
(Control function): For any $t\in[0,T]$, $w\in\Omega$, $\mathbf{x}\in\mathbb{R}^{d}$, and $\theta,\theta^{\prime}\in\Theta$, the control function $\alpha$ is $L$-Lipschitz function, $|\alpha(t,\mathbf{x},\theta)-\alpha(t,\mathbf{x},\theta^{\prime})|\leq L|
\theta-\theta^{\prime}|$. Moreover it satisfy $\mathbb{E}\left[\int_{0}^{T}|\alpha_{s}^{2}|ds\right]<\infty$.

###### Definition A.1 (Infinitesimal Generator) .

Let us consider an Itô diffusion process of the form:

$$ $d\mathbf{X}_{t}=b(t,\mathbf{X}_{t})dt+\sigma(t)^{\top}d\mathbf{W}_{t},$ (29) $$

Then, an infinitesimal generator of the above diffusion process is given by:

$$ $\displaystyle\mathcal{A}_{t}f=\lim_{t\downarrow 0^{+}}\frac{\mathbb{E}\left[f( \mathbf{X}_{t})\right]-f(\mathbf{x})}{t}$ $\displaystyle=\nabla_{\mathbf{x}}f^{\top}b+\frac{1}{2}\text{Trace}\left[\sigma \sigma^{\top}\nabla_{\mathbf{xx}}f\right].$ (30) $$

###### Definition A.2 (Path Measure) .

Let us consider a path-sequence of random variables $\{\mathbf{X}_{t_{i}}\}_{i\in[1:N]}$ over an interval $0=t_{i}\leq\cdots\leq t_{N}=T$, where each $\mathbf{X}_{t_{i}}$ taking values in measurable space $(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$ with a Borel $\sigma$-algebra $\mathcal{B}(\mathbb{R}^{d})$ . Then, the set of random variables $\{\mathbf{X}_{t_{i}}\}_{i\in[1:N]}$ is given by the probability measure $\mathbb{P}\in C([0,T],\mathbb{R}^{d})$:

$$ $\mathbb{P}(\mathbf{X}_{t_{0}}\in d\mathbf{x}_{t_{0}},\cdots,\mathbf{X}_{t_{N}} \in d\mathbf{x}_{t_{N}})=\mathbb{P}(d\mathbf{x}_{0})\prod_{i=1}^{N}\mathbf{P}_ {i}(\mathbf{x}_{t_{i-1}},d\mathbf{x}_{t_{i}}),$ (31) $$

where $\{\mathbf{P}_{i}\}_{i=0}^{N}$ is a sequence of probability kernels from $(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$ to $(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$, for any event $A\in\mathcal{B}(\mathbb{R}^{d})$, $\mathbf{P}_{i}(\mathbf{x}_{t_{i-1}},A)=\int_{A}\mathbf{p}_{i}(\mathbf{x}_{t_{i
-1}},\mathbf{x}_{t_{i}})d\mathbf{x}_{t_{i}}$, where $\mathbf{p}_{i}(\mathbf{x}_{t_{i-1}},\mathbf{x}_{t_{i}}):=\mathbf{p}(t_{i},
\mathbf{x}_{t_{i}}|t_{i-1},\mathbf{x}_{t_{i-1}})$ is a transition density obtained by a solution of the Fokker-Placnk equation (Risken & Frank, ):

$$ $\partial_{t}\mathbf{p}_{t}(\mathbf{x}_{t})=\mathcal{A}^{\star}_{t}\mathbf{p}_{ t}=-\nabla_{\mathbf{x}}\cdot(b\mathbf{p}_{t})+\frac{1}{2}\textit{Trace}\left[ \sigma\sigma^{\top}\nabla_{\mathbf{xx}}\mathbf{p}_{t}\right],$ (32) $$

where $\mathcal{A}^{\star}$ is an adjoint operator of the generator in (30) and $\mathbf{p}_{t}$ is the Radon-Nikodym density of $\mu_{t}$ with respect to the Lebesgue measure. By taking $N\rightarrow\infty$, we get the path measure $\mathbb{P}(\mathbf{X}_{0:T}\in d\mathbf{x}_{0:T})$ which describes the weak solutions of the SDE of the form in equation (29).

###### Lemma A.3 (Itô’s formula) .

Let $v(t,x)$ be $C^{1}$ in $t$ and $C^{2}$ in $x$ and let $\mathbf{X}_{t}$ be the Itô diffusion process of the form in equation (29).
Then, the stochastic process $v(t,\mathbf{X}_{t})$ is also an Itô diffusion process satisfying:

$$ $dv(t,\mathbf{X}_{t})=\left[\partial_{t}v(t,\mathbf{X}_{t})+\mathcal{A}_{t}v(t, \mathbf{X}_{t})\right]dt+\nabla_{\mathbf{x}}v(t,\mathbf{X}_{t})^{\top}\sigma(t )d\mathbf{W}_{t}.$ (33) $$

###### Theorem A.4 (Girsanov Theorem) .

Consider the two Itô diffusion processes of form

$$ $\displaystyle d\mathbf{X}_{t}=b(t,\mathbf{X}_{t})dt+\sigma(t,\mathbf{X}_{t})^{ \top}d\mathbf{W}_{t},\quad t\in[0,T],$ (34) $\displaystyle d\mathbf{Y}_{t}=\tilde{b}(t,\mathbf{Y}_{t})dt+\sigma(t,\mathbf{Y }_{t})^{\top}d\mathbf{W}_{t},\quad t\in[0,T],$ (35) $$

where both drift functions $b,\tilde{b}$ and the diffusion function $\sigma$ assumed to be invertible are adapted to $\mathcal{F}_{t}$ and $\mathbf{W}_{[0,T]}$ is $\mathbb{P}$-Wiener process. Moreover, consider $\mathbb{P}$ as the path measures induced by (34). Let us define $H_{t}:=\sigma^{-1}(\tilde{b}-b)$ which is assumed to be satisfying the Novikov’s condition ($\textit{i}.\textit{e}.,\mathbb{E}_{\mathbb{P}}\left[\exp\left(\frac{1}{2}\int_
{0}^{T}\left\lVert H_{s}\right\rVert^{2}ds\right)\right]<\infty$), and the $\mathbb{P}$-martingale process

$$ $\mathbf{M}_{t}:=\exp\left(\int_{0}^{1}H_{s}^{\top}d\mathbf{W}_{s}-\frac{1}{2} \int_{0}^{t}\left\lVert H_{s}\right\rVert^{2}ds\right)$ (36) $$

satisfies $\mathbb{E}_{\mathbb{P}}[\mathbf{M}_{T}]=1$. Then for the path measure $\mathbb{Q}$ given as $d\mathbb{Q}=\mathbf{M}_{T}d\mathbb{P}$,
the process $\tilde{W}_{t}=W_{t}-\int_{0}^{t}H_{s}ds$ is a $\mathbb{Q}$-Wiener process and $\mathbf{Y}_{t}$ can be represented as

$$ $d\mathbf{Y}_{t}=b(t,\mathbf{Y}_{t})dt+\sigma(t,\mathbf{Y}_{t})^{\top}d\tilde{W }_{t},\quad t\in[0,T].$ (37) $$

Therefore $\mathbb{Q}$-law of the process $\mathbf{Y}_{t}$ is same as $\mathbb{P}$-law of the process $\mathbf{X}_{t}$.

###### Definition A.1 (Infinitesimal Generator) .

###### Definition A.2 (Path Measure) .

###### Lemma A.3 (Itô’s formula) .

###### Theorem A.4 (Girsanov Theorem) .

### A.1 Proof of Theorem 3.1

We start the section by showing the normalizing property of $\{f_{i}\}_{i\in[1:k]}$ in (3). By definition, it satisfied that

$$ $\prod_{i=1}^{k}\mathbf{L}_{i}(g_{i})=\prod_{i=1}^{k}\int_{\mathbb{R}^{d}}g_{i} (\mathbf{y}_{t_{i}}|\mathbf{x}_{t_{i}})d\mathbb{P}(\mathbf{x}_{0:T})\stackrel{ {\scriptstyle(i)}}{{=}}\mathbb{E}_{\mathbb{P}}\left[\prod_{i=1}^{k}g_{i}( \mathbf{y}_{t_{i}}|\mathbf{x}_{t_{i}})\right]=\mathbf{Z}(\mathcal{H}_{t_{k}}),$ (38) $$

where $(i)$ follows from the conditional indenpendency of $\mathbf{y}_{t_{i}}$ given $\mathbf{x}_{t_{i}}$ for all $i\in[1:k]$. Hence, we get the normalizing property:

$$ $\mathbb{E}_{\mathbb{P}}\left[\prod_{i=1}^{k}f_{i}(\mathbf{x}_{t_{i}})\right]= \mathbb{E}_{\mathbb{P}}\left[\frac{\prod_{i=1}^{k}g_{i}(\mathbf{y}_{t_{i}}| \mathbf{x}_{t_{i}})}{\prod_{i=1}^{k}\mathbf{L}_{t_{i}}(g_{t_{i}})}\right]= \frac{1}{\mathbf{Z}(\mathcal{H}_{t_{k}})}\mathbb{E}_{\mathbb{P}}\left[\prod_{i =1}^{k}g_{i}(\mathbf{y}_{t_{i}}|\mathbf{x}_{t_{i}})\right]=1.$ (39) $$

###### Theorem 3.1 (Multi-marginal Doob’s $h$ -transform) .

Let us define a sequence of functions $\{h_{i}\}_{i\in[1:k]}$, where each $h_{i}:[t_{i-1},t_{i})\times\mathbb{R}^{d}\to\mathbb{R}_{+}$, for all $i\in[1:k]$, is a conditional expectation $h_{i}(t,\mathbf{x}_{t}):=\mathbb{E}_{\mathbb{P}}\left[\prod_{j\geq i}^{k}f_{j}
(\mathbf{y}_{t_{j}}|\mathbf{X}_{t_{j}})|\mathbf{X}_{t}=\mathbf{x}_{t}\right]$,
where $\{f_{i}\}_{i\in[1:k]}$ is defined in (3). Now, we define a function $h:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}_{+}$ by integrating the functions $\{h_{i}\}_{i\in[1:k]}$,

$$ $h(t,\mathbf{x}):=\sum_{i=1}^{k}h_{i}(t,\mathbf{x})\mathbf{1}_{[t_{i-1},t_{i})} (t).$ (40) $$

Then, with the initial condition $\mu^{\star}_{0}(d\mathbf{x}_{0})=h_{1}(t_{0},\mathbf{x}_{0})\mu_{0}(d\mathbf{x
}_{0})$,
the solution of the following conditional SDE inducing
the posterior path measure $\mathbb{P}^{\star}$ in (2):

$$ $\text{(Conditioned State)}\quad d\mathbf{X}^{\star}_{t}=\left[b(t,\mathbf{X}^{ \star}_{t})+\nabla_{\mathbf{x}}\log h(t,\mathbf{X}^{\star}_{t})\right]dt+d \mathbf{W}_{t}$ (41) $$

###### Proof.

We start with the interval $[t_{i-1},t_{i})$ without loss of generality. For all $t\in[t_{i-1},t_{i})$ and for any $A_{t_{i}}\subset\mathbb{R}^{d}$, the transition kernel of the conditioned process is defined by the transition kernel of the original Markov process $\mathbf{P}_{i}$ and $h$-function defined in (3.1):

$$ $\mathbf{P}^{h_{i}}_{i}(\mathbf{X}_{t_{i}}\in A|\mathbf{X}_{t}=\mathbf{x}):= \mathbf{P}^{h_{i}}_{i}(\mathbf{x}_{t},A)=\frac{h_{i}(t_{i},\mathbf{X}_{t_{i}}) }{h_{i}(t,\mathbf{x}_{t})}\mathbf{P}_{t_{i}}(\mathbf{x}_{t},d\mathbf{x}_{t_{i} }).$ (42) $$

By the definition of $h$, the transition kernel $\mathbf{P}^{h}_{t_{i}}(\mathbf{x}_{t},A)$ is a probability kernel:

$$ $\displaystyle\int_{\mathbb{R}^{d}}\mathbf{P}^{h_{i}}_{i}(\mathbf{x}_{t},d \mathbf{x}_{t_{i}})$ $\displaystyle=\int_{\mathbb{R}^{d}}\frac{h_{i}(t_{i},\mathbf{X}_{t_{i}})}{h_{i }(t,\mathbf{x}_{t})}\mathbf{P}_{i}(\mathbf{x}_{t},d\mathbf{x}_{t_{i}})$ (43) $\displaystyle=\frac{\int_{\mathbb{R}^{d}}h_{i}(t_{i},\mathbf{X}_{t_{i}}) \mathbf{P}_{i}(\mathbf{x}_{t},d\mathbf{x}_{t_{i}})}{h_{i}(t,\mathbf{x}_{t})}$ (44) $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\frac{\int_{\mathbb{R}^{d}}f_{i} (\mathbf{y}_{t_{i}}|\mathbf{X}_{t_{i}})h_{i+1}(t_{i},\mathbf{X}_{t_{i}}) \mathbf{P}_{i}(\mathbf{x}_{t},d\mathbf{x}_{t_{i}})}{{\mathbb{E}_{\mathbb{P}} \left[\prod_{j=i}^{k}f_{j}(\mathbf{y}_{t_{j}}|\mathbf{X}_{t_{j}})|\mathbf{X}_{ t}=\mathbf{x}_{t}\right]}}$ (45) $\displaystyle=\frac{\mathbb{E}_{\mathbb{P}}\left[\prod_{j=i}^{k}f_{j}(\mathbf{ y}_{t_{j}}|\mathbf{X}_{t_{j}})|\mathbf{X}_{t}=\mathbf{x}_{t}\right]}{\mathbb{E }_{\mathbb{P}}\left[\prod_{j=i}^{k}f_{j}(\mathbf{y}_{t_{j}}|\mathbf{X}_{t_{j}} )|\mathbf{X}_{t}=\mathbf{x}_{t}\right]}=1,$ (46) $$

where $(i)$ follows from the recursion established by the definition of $h$ in Theorem 3.1:

$$ $h_{i}(t_{i},\mathbf{x}_{t_{i}})=f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}_{t_{i}})h_ {i+1}(t_{i},\mathbf{x}_{t_{i}}),\quad\forall i\in[1:k-1].$ (47) $$

Now, the infinitesimal generator for $\mathbf{P}^{h}$ can be computed for any $\varphi\in C^{1,2}([0,T]\times\mathbb{R}^{d})$ and for all $\mathbb{P}$-almost $\mathbf{x}\in\mathbb{R}^{d}$:

$$ $\displaystyle\mathcal{A}^{h_{i}}_{t}\varphi_{t}$ $\displaystyle=\lim_{s\downarrow 0}\frac{\mathbb{E}_{\mathbb{P}^{h}}\left[ \varphi(t_{s},\mathbf{X}_{t+s})|\mathbf{X}_{t}=\mathbf{x}\right]-\varphi(t, \mathbf{x})}{s}$ (48) $\displaystyle=\lim_{s\downarrow 0}\frac{\mathbb{E}_{\mathbb{P}}\left[[\varphi( t_{s},\mathbf{X}_{t+s})-\varphi(t,\mathbf{x})]\frac{\mathbf{P}^{h_{i}}_{i}( \mathbf{x},d\mathbf{x}_{t+s})}{\mathbf{P}_{i}(\mathbf{x},d\mathbf{x}_{t+s})}| \mathbf{X}_{t}=\mathbf{x}\right]}{s}$ (49) $\displaystyle=\lim_{s\downarrow 0}\frac{\mathbb{E}_{\mathbb{P}}\left[[\varphi( t_{s},\mathbf{X}_{t+s})-\varphi(t,\mathbf{x})]\frac{h_{i}(t+s,\mathbf{X}_{t+s} )}{h_{i}(t,\mathbf{x})}|\mathbf{X}_{t}=\mathbf{x}\right]}{s}$ (50) $\displaystyle=\lim_{s\downarrow 0}\frac{\mathbb{E}_{\mathbb{P}}\left[[\varphi( t_{s},\mathbf{X}_{t+s})-\varphi(t,\mathbf{x})]\left[\frac{h_{i}(t+s,\mathbf{X} _{t+s})-h_{i}(t,\mathbf{x})}{h_{i}(t,\mathbf{x})}+1\right]|\mathbf{X}_{t}= \mathbf{x}\right]}{s}$ (51) $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\mathcal{A}_{t}\varphi_{t}+\lim_ {s\downarrow 0}\frac{\mathbb{E}_{\mathbb{P}}\left[[\varphi(t_{s},\mathbf{X}_{t +s})-\varphi(t,\mathbf{x})]\left[h_{i}(t+s,\mathbf{X}_{t+s})-h_{i}(t,\mathbf{x })\right]|\mathbf{X}_{t}=\mathbf{x}\right]}{sh_{i}(t,\mathbf{x})},$ (52) $$

where $(i)$ follows from the definition of the infinitesimal generator.
Now, we shall compute the second term of the RHS of the equation (52) adapted from (Léonard, ). By employing basic stochastic calculus, for a stochastic process $\varphi_{t}=\varphi(t,\mathbf{X}_{t})$ and $h_{i,t}=h_{i}(t,\mathbf{X}_{t})$,

$$ $\displaystyle\varphi_{t+s}h_{i,t+s}=\varphi_{0}h_{i,0}+\int_{0}^{t+s}\varphi_{ u}dh_{i,u}+\int_{0}^{t+s}h_{i,u}d\varphi_{u}+[\varphi,h_{i}]_{t+s}$ (53) $\displaystyle\varphi_{t}h_{i,t}=\varphi_{0}h_{i,0}+\int_{0}^{t}\varphi_{u}dh_{ i,u}+\int_{0}^{t}h_{i,u}d\varphi_{u}+[\varphi,h_{i}]_{t},$ (54) $$

where $[\varphi,h_{i}]_{t}=\int_{0}^{t}d\varphi_{t}dh_{i,t}$ is quadratic variation of $\varphi$ and $h_{i}$.
By subtracting equation (54) from equation (53),

$$ $\varphi_{t+s}h_{i,t+s}-\varphi_{t}h_{i,t}=\int_{t}^{t+s}\varphi_{u}dh_{i,u}+ \int_{t}^{t+s}h_{i,u}d\varphi_{u}+[\varphi,h_{i}]_{t+s}-[\varphi,h_{i}]_{t}$ (55) $$

Applying integration by parts leads to the following equation

$$ $\displaystyle(\varphi_{t+s}-\varphi_{t})(h_{i,t+s}-h_{i,t})$ $\displaystyle=\varphi_{t+s}h_{i,t+s}-\varphi_{t}h_{i,t}-\varphi_{t}(h_{i,t+s}- h_{i,t})-h_{i,t}(\varphi_{t+s}-\varphi_{t})$ (56) $\displaystyle=\int_{t}^{t+s}(\varphi_{u}-\varphi_{t})dh_{i,u}+\int_{t}^{t+s}(h _{i,u}-h_{i,t})d\varphi_{u}+[\varphi,h_{i}]_{t+s}-[\varphi,h_{i}]_{t},$ $$

where $\varphi_{t}(h_{i,t+s}-h_{i,t})=\int_{t}^{t+s}\varphi_{t}dh_{i,u}$. Moreover, by applying Itô’s formula, we get

$$ $d\varphi_{t}=\mathcal{A}_{t}\varphi_{t}dt+(\nabla_{\mathbf{x}}\varphi)^{\top}d \mathbf{W}_{t},\quad dh_{i,t}=\mathcal{A}_{t}h_{i,t}dt+(\nabla_{\mathbf{x}}h_{ i,t})^{\top}d\mathbf{W}_{t}.$ (57) $$

Therefore, since $\mathbf{W}_{t}$ is $\mathbb{P}$-martingale,

$$ $\displaystyle\mathbb{E}_{\mathbb{P}}\left[(\varphi_{t+s}-\varphi_{t})(h_{i,t+s }-h_{i,t})|\mathbf{X}_{t}=\mathbf{x}\right]$ (58) $\displaystyle=\mathbb{E}_{\mathbb{P}}\left[\int_{s}^{t+s}(\varphi_{u}-\varphi_ {t})\mathcal{A}_{u}h_{i,u}du+\int_{s}^{t+s}(h_{i,u}-h_{i,t})\mathcal{A}_{u} \varphi_{u}du+[\varphi,h_{i}]_{t+s}-[\varphi,h_{i}]_{t}|\mathbf{X}_{t}=\mathbf {x}\right]$ (59) $\displaystyle=\underbrace{\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[\int_{t} ^{t+s}(\varphi_{u}-\varphi_{t})\mathcal{A}_{u}h_{i,u}du\right]}_{\textbf{(A)}} +\underbrace{\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[\int_{t}^{t+s}(h_{i,u }-h_{i,t})\mathcal{A}_{u}\varphi_{u}du\right]}_{\textbf{(B)}}+\underbrace{ \mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[[\varphi,h_{i}]_{t+s}-[\varphi,h_{ i}]_{t}\right]}_{\textbf{(C)}}.$ $$

For a first term (A), the Hölder’s inequality with $1/p+1/q=1$ and $p,q\geq 1$ yields,

$$ $\displaystyle\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[\int_{t}^{t+s}( \varphi_{u}-\varphi_{t})\mathcal{A}_{u}h_{i,u}du\right]$ $\displaystyle\leq\left(\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[\int_{t}^{t +s}|\varphi_{u}-\varphi_{t}|^{q}du\right]\right)^{1/q}\left(\mathbb{E}^{t, \mathbf{x}}_{\mathbb{P}}\left[\int_{t}^{t+s}|\mathcal{A}_{t}h_{i,u}|^{p}du \right]\right)^{1/p}$ (60) $\displaystyle=\left(\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[\int_{t}^{t+s} |\varphi_{u}-\varphi_{t}|^{q}du\right]\right)^{1/q}\left(\int_{t}^{t+s}\mathbb {E}^{t,\mathbf{x}}_{\mathbb{P}}\left[|\mathcal{A}_{t}h_{i,u}|^{p}\right]du \right)^{1/p}$ (61) $$

Given the bounded and continuous function $\varphi$, the dominated convergence theorem yields $(\lim_{s\downarrow 0}\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[\int_{t}^{t+s
}|\varphi_{u}-\varphi_{t}|^{q}du\right])^{1/q}=(\mathbb{E}^{t,\mathbf{x}}_{
\mathbb{P}}\left[\lim_{s\downarrow 0}\int_{t}^{t+s}|\varphi_{u}-\varphi_{t}|^{
q}du\right])^{1/q}=0$ and since $h_{i}\in C^{1,2}([t_{i-1},t_{i}),\mathbb{R}^{d})$ and boundedness of $b$ in Assumptions A, following inequality holds

$$ $|\mathcal{A}_{u}h_{i,u}|^{p}\leq|\partial_{t}h_{i,t}|^{p}+|(\nabla_{\mathbf{x} }h_{i,t}^{T})b|^{p}+|\frac{1}{2}\textit{Trace}\left[\nabla_{\mathbf{x}\mathbf{ x}}h_{i,t}\right]|^{p}<\infty,$ (62) $$

for any $u\in[t_{i-1},t_{i})$ and $\mathbb{P}$ almost surely. In other words, $\sup_{u\in[t,t+s]}\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[|\mathcal{A}_{u}
h_{i,u}|^{p}\right]<\infty,\forall t\in[0,T],s>0,\text{and}\;p>1$. Therefore we get $\lim_{s\downarrow 0}\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[\int_{t}^{t+s}
(\varphi_{u}-\varphi_{t})\mathcal{A}_{u}h_{i,u}du\right]=0$
and we can get the similar result for the second term (B), $\textit{i}.\textit{e}.,\lim_{s\downarrow 0}\mathbb{E}^{t,\mathbf{x}}_{\mathbb{
P}}\left[\int_{t}^{t+s}(h_{i,u}-h_{i,t})\mathcal{A}_{u}\varphi_{u}du\right]=0$.

For the last term (C), by the definition of the quadratic variation of $\varphi$ and $h_{i}$ yields:

$$ $\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[[\varphi,h_{i}]_{t+s}-[\varphi,h_{ i}]_{t}\right]=\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[\int_{t}^{t+s}d \varphi_{u}dh_{i,u}\right]=\mathbb{E}^{t,\mathbf{x}}_{\mathbb{P}}\left[\int_{t }^{t+s}(\nabla_{\mathbf{x}}\varphi_{u})^{\top}\nabla_{\mathbf{x}}h_{i,u}du\right]$ (63) $$

Subsequently, by taking the limit from (52), we get the following result:

$$ $\displaystyle\lim_{s\downarrow 0}\frac{\mathbb{E}_{\mathbb{P}}\left[[\varphi(t _{s},\mathbf{X}_{t+s})-\varphi(t,\mathbf{x})]\left[h_{i}(t+s,\mathbf{X}_{t+s}) -h_{i}(t,\mathbf{x})\right]|\mathbf{X}_{t}=\mathbf{x}\right]}{sh_{i}(t,\mathbf {x})}$ (64) $\displaystyle=\lim_{s\downarrow 0}\frac{\mathbb{E}_{\mathbb{P}}\left[\int_{t}^ {t+s}(\nabla_{\mathbf{x}}\varphi(u,\mathbf{X}_{u}))^{\top}\nabla_{\mathbf{x}}h _{i}(u,\mathbf{X}_{u})du|\mathbf{X}_{t}=\mathbf{x}\right]}{sh_{i}(t,\mathbf{x})}$ (65) $\displaystyle=(\nabla_{\mathbf{x}}\varphi(t,\mathbf{X}_{t}))^{\top}\nabla_{ \mathbf{x}}\log h_{i}(t,\mathbf{X}_{t}),$ (66) $$

Therefore, the infinitesimal generator for $\mathbf{P}^{h_{i}}_{i}$ is defined by:

$$ $\mathcal{A}^{h_{i}}_{t}\varphi_{t}=\mathcal{A}_{t}\varphi_{t}+(\nabla_{\mathbf {x}}\varphi_{t})^{\top}\nabla_{\mathbf{x}}\log h_{i,t}$ (67) $$

which shows that the conditioned SDE for an interval $[t_{i-1},t_{i})$ is given by

$$ $d\mathbf{X}^{h}_{t}=\left[b(t,\mathbf{X}_{t})+\nabla_{\mathbf{x}}\log h_{i}(t, \mathbf{X}_{t})\right]dt+d\mathbf{W}_{t}$ (68) $$

Hence, integrating the generators over the entire interval yields:

$$ $\displaystyle\mathcal{A}^{h}_{t}\varphi$ $\displaystyle=\mathcal{A}_{t}\varphi_{t}+\sum_{i=1}^{k}\left[(\nabla_{\mathbf{ x}}\varphi_{t})^{\top}\nabla_{\mathbf{x}}\log h_{i,t}\right]\mathbf{1}_{[t_{i- 1},t_{i})}(t)$ (69) $\displaystyle=\mathcal{A}_{t}\varphi_{t}+(\nabla_{\mathbf{x}}\varphi_{t})^{ \top}\nabla_{\mathbf{x}}\log h_{t}.$ (70) $$

It implies that the conditioned SDE for the entire interval $[0,T]$ is given by:

$$ $d\mathbf{X}^{h}_{t}=\left[b(t,\mathbf{X}_{t})+\nabla_{\mathbf{x}}\log h(t, \mathbf{X}_{t})\right]dt+d\mathbf{W}_{t}.$ (71) $$

Now, assume that $\mathbf{X}^{h}_{t}\sim\mu^{\star}_{0}(\mathbf{x})$ where $\mu^{\star}_{0}(\mathbf{x})$ is absolutely continuous with $\mu_{0}(\mathbf{x})$. Then, the path measure induced by $\mathbf{X}^{h}_{t}$ can be computed as:

$$ $\displaystyle d\mathbb{P}^{h}(\mathbf{x}_{0:T})$ $\displaystyle=d\mu^{\star}_{0}(\mathbf{x}_{0})\prod_{i=1}^{k}\left[\prod_{j=1} ^{N}\mathbf{P}_{i(j)}^{h_{i}}(\mathbf{x}_{t_{i(j-1)}},d\mathbf{x}_{t_{i}(j)})\right]$ (72) $\displaystyle=d\mu^{\star}_{0}(\mathbf{x}_{0})\prod_{i=1}^{k}\frac{h_{i}(t_{i} ,\mathbf{x}_{t_{i}})}{h_{i}(t_{i-1},\mathbf{x}_{t_{i-1}})}\left[\prod_{j=1}^{N }\mathbf{P}_{i(j)}(\mathbf{x}_{t_{i(j-1)}},d\mathbf{x}_{t_{i}(j)})\right]$ (73) $\displaystyle=d\mu^{\star}_{0}(\mathbf{x}_{0})\prod_{i=1}^{k}\frac{h_{i+1}(t_{ i},\mathbf{x}_{t_{i}})f_{i}(\mathbf{y}_{t_{i}}|\mathbf{x}_{t_{i}})}{h_{i}(t_{i -1},\mathbf{x}_{t_{i-1}})}\left[\prod_{j=1}^{N}\mathbf{P}_{i(j)}(\mathbf{x}_{t _{i(j-1)}},d\mathbf{x}_{t_{i}(j)})\right]$ (74) $\displaystyle\stackrel{{\scriptstyle N\uparrow\infty}}{{=}}\frac{d\mu^{\star}_ {0}}{d\mu_{0}}(\mathbf{x}_{0})\frac{h_{k+1}(t_{k},\mathbf{x}_{t_{k}})}{h_{1}(t _{0},\mathbf{x}_{0})}\prod_{i=1}^{k}f_{i}(\mathbf{y}_{t_{i}}|\mathbf{x}_{t_{i} })d\mathbb{P}(\mathbf{x}_{0:T})$ (75) $$

where, for all $i\in[1:k]$, we define a increasing sequence $\{i(j)\}_{j\in[0:N]}$ with $i(0)=i-1$, $i(1)=i-1+\frac{1}{N}$ and $i(N)=i$. Hence, for a $d\mu_{0}^{\star}(\mathbf{x}_{0})=h_{1}(t_{0},\mathbf{x}_{0})d\mu_{0}(\mathbf{x
}_{0})$ and $h_{k+1}=1$ yields

$$ $\displaystyle d\mathbb{P}^{h}(\mathbf{x}_{0:T})$ $\displaystyle=\prod_{i=1}^{k}f_{i}(\mathbf{y}_{t_{i}}|\mathbf{x}_{t_{i}})d \mathbb{P}(\mathbf{x}_{0:T})$ (76) $\displaystyle=\frac{1}{\mathbf{Z}(\mathcal{H}_{t_{k}})}\prod_{i=1}^{k}g_{i}( \mathbf{y}_{t_{i}}|\mathbf{x}_{t_{i}})d\mathbb{P}(\mathbf{x}_{0:T})$ (77) $\displaystyle=d\mathbb{P}^{\star}(\mathbf{x}_{0:T}).$ (78) $$

It concludes the proof.
∎

###### Theorem 3.1 (Multi-marginal Doob’s h ℎ h italic_h -transform) .

###### Proof.

### A.2 Proof of Theorem 3.2

###### Theorem 3.2 (Dynamic Programming Principle) .

Let us consider a sequence of left continuous functions $\{\mathcal{V}_{i}\}_{i\in[1:k+1]}$, where each $\mathcal{V}_{i}\in C^{1,2}([t_{i-1},t_{i})\times\mathbb{R}^{d})$

$$ $\mathcal{V}_{i}(t,\mathbf{x}_{t}):=\min_{\alpha\in\mathbb{A}}\mathbb{E}_{ \mathbb{P}^{\alpha}}\left[\int_{t_{i-1}}^{t_{i}}\frac{1}{2}\left\lVert\alpha_{ s}\right\rVert^{2}ds-\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}} )+\mathcal{V}_{i+1}(t_{i},\mathbf{X}^{\alpha}_{t_{i}})|\mathbf{X}_{t}=\mathbf{ x}_{t}\right],$ (79) $$

for all $i\in[1:k]$ and $\mathcal{V}_{k+1}=0$. Then, for any $0\leq t\leq u\leq T$, the value function $\mathcal{V}$ for the cost function in (7) satisfying the recursion defined as follows:

$$ $\mathcal{V}(t,\mathbf{x}_{t})=\min_{\alpha\in\mathbb{A}}\mathbb{E}_{\mathbb{P} ^{\alpha}}\left[\int_{t}^{t_{I(u)}}\frac{1}{2}\left\lVert\alpha_{s}\right \rVert^{2}ds-\sum_{i:\{t\leq t_{i}\leq t_{I(u)}\}}\log f_{i}(\mathbf{y}_{t_{i} }|\mathbf{X}^{\alpha}_{t_{i}})+\mathcal{V}_{I(u)+1}(t_{I(u)},\mathbf{X}^{ \alpha}_{t_{I(u)}})|\mathbf{X}^{\alpha}_{t}=\mathbf{x}_{t}\right],$ (80) $$

with the indexing function $I(u)=\max\{i\in[1:k]|t_{i}\leq u\}.$

###### Proof.

Following the approach used in the proof of the standard dynamic programming principle with the flow property induced by Markov control (Van Handel, ), we can apply similar methods to our cost function. We start the proof by establishing the recursion of $\mathcal{J}$. Let us define the sequence of left continuous cost functions $\{\mathcal{J}_{i}\}_{i\in[1:k+1]}$:

$$ $\mathcal{J}_{i}(t,\mathbf{x}_{t},\alpha):=\mathbb{E}^{t,\mathbf{x}_{t}}_{ \mathbb{P}^{\alpha}}\left[\int_{t_{i-1}}^{t_{i}}\frac{1}{2}\left\lVert\alpha_{ s}\right\rVert^{2}ds-\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}} )+\mathcal{J}_{i+1}(t_{i},\mathbf{X}^{\alpha}_{t_{i}},\alpha)\right],$ (81) $$

where we denote $\mathbb{E}^{\mathbf{t},\mathbf{x}}_{\mathbb{P}}[\cdot]=\mathbb{E}_{\mathbb{P}}
[\cdot|\mathbf{X}_{t}=\mathbf{x}]$, for all $i\in[1:k]$ and $\mathcal{J}_{k+1}=0$. Since $\mathbb{P}^{\alpha}$ is Markov process, it satisfying following recursion, for any $0\leq t\leq u\leq T$,

$$ $\mathcal{J}(t,\mathbf{x}_{t},\alpha)=\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P} ^{\alpha}}\left[\int_{t}^{t_{I(u)}}\frac{1}{2}\left\lVert\alpha_{s}\right \rVert^{2}ds-\sum_{i:\{t\leq t_{i}\leq t_{I(u)}\}}\log f_{i}(\mathbf{y}_{t_{i} }|\mathbf{X}^{\alpha}_{t_{i}})+\mathcal{J}_{I(u)+1}(t_{I(u)},\mathbf{X}^{ \alpha}_{t_{I(u)}},\alpha)\right],$ (82) $$

with the indexing function $I(u)=\max\{i\in[1:k]|t_{i}\leq u\}$.

For any $\epsilon>0$, there exists a control $\alpha^{\prime}\in\mathbb{A}[t,T]$ such that

$$ $\displaystyle\mathcal{V}(t,\mathbf{x})+\epsilon\geq\mathcal{J}(t,\mathbf{x}, \alpha^{\prime})$ (83) $\displaystyle=\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha^{\prime}}} \left[\int_{t}^{t_{I(u)}}\frac{1}{2}\left\lVert\alpha^{\prime}_{s}\right\rVert ^{2}ds-\sum_{i:\{t\leq t_{i}\leq t_{I(u)}\}}\log f_{i}(\mathbf{y}_{t_{i}}| \mathbf{X}^{\alpha^{\prime}}_{t_{i}})+\mathcal{J}_{I(u)+1}(t_{I(u)},\mathbf{X} ^{\alpha^{\prime}}_{t_{I(u)}},\alpha^{\prime})\right]$ (84) $\displaystyle\geq\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha^{\prime}}} \left[\int_{t}^{t_{I(u)}}\frac{1}{2}\left\lVert\alpha^{\prime}_{s}\right\rVert ^{2}ds-\sum_{i:\{t\leq t_{i}\leq t_{I(u)}\}}\log f_{i}(\mathbf{y}_{t_{i}}| \mathbf{X}^{\alpha^{\prime}}_{t_{i}})+\mathcal{V}_{I(u)+1}(t_{I(u)},\mathbf{X} ^{\alpha^{\prime}}_{t_{I(u)}})\right]$ (85) $\displaystyle\geq\inf_{\alpha^{\prime}\in\mathbb{A}[t,T]}\mathbb{E}^{t,\mathbf {x}_{t}}_{\mathbb{P}^{\alpha^{\prime}}}\left[\int_{t}^{t_{I(u)}}\frac{1}{2} \left\lVert\alpha^{\prime}_{s}\right\rVert^{2}ds-\sum_{i:\{t\leq t_{i}\leq t_{ I(u)}\}}\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha^{\prime}}_{t_{i}})+ \mathcal{V}_{I(u)+1}(t_{I(u)},\mathbf{X}^{\alpha^{\prime}}_{t_{I(u)}})\right]$ (86) $$

Since $\epsilon$ was arbitrary, limiting $\epsilon\to 0$ we get:

$$ $\mathcal{V}(t,\mathbf{x})\geq\inf_{\alpha^{\prime}\in\mathbb{A}[t,T]}\mathbb{E }^{t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha^{\prime}}}\left[\int_{t}^{t_{I(u)}} \frac{1}{2}\left\lVert\alpha^{\prime}_{s}\right\rVert^{2}ds-\sum_{i:\{t\leq t_ {i}\leq t_{I(u)}\}}\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha^{\prime}}_ {t_{i}})+\mathcal{V}_{I(u)+1}(t_{I(u)},\mathbf{X}^{\alpha^{\prime}}_{t_{I(u)}} )\right].$ (87) $$

For the reverse direction, consider the control $\tilde{\alpha}\in\mathbb{A}[t,T]$ obtained from integrating:

$$ $\tilde{\alpha}_{s}:=\begin{cases}\alpha^{1}_{s},\quad s\in[t,t_{I(u)})\\ \alpha^{2}_{s}\quad s\in[t_{I(u)},T].\end{cases}$ (88) $$

Then, by following the definition of the value function

$$ $\displaystyle\mathcal{J}(t,\mathbf{x},\tilde{\alpha})\geq\inf_{\alpha^{2}\in \mathbb{A}[t_{I(u)},T]}\mathcal{J}(t,\mathbf{x},\tilde{\alpha})$ (89) $\displaystyle=\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha^{1}}}\left[ \int_{t}^{t_{I(u)}}\frac{1}{2}\left\lVert\alpha^{1}_{s}\right\rVert^{2}ds-\sum _{i:\{t\leq t_{i}\leq t_{I(u)}\}}\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{ \alpha^{1}}_{t_{i}})+\mathcal{V}_{I(u)+1}(t_{I(u)},\mathbf{X}^{\alpha^{1}}_{t_ {I(u)}})\right]$ (90) $\displaystyle\geq\min_{\alpha^{1}\in\mathbb{A}[t,t_{I(u)})}\mathbb{E}^{t, \mathbf{x}_{t}}_{\mathbb{P}^{\alpha^{1}}}\left[\int_{t}^{t_{I(u)}}\frac{1}{2} \left\lVert\alpha^{1}_{s}\right\rVert^{2}ds-\sum_{i:\{t\leq t_{i}\leq t_{I(u)} \}}\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha^{1}}_{t_{i}})+\mathcal{V}_ {I(u)+1}(t_{I(u)},\mathbf{X}^{\alpha^{1}}_{t_{I(u)}})\right]$ (91) $\displaystyle=\inf_{\tilde{\alpha}\in\mathbb{A}[t,T]}\mathbb{E}^{t,\mathbf{x}_ {t}}_{\mathbb{P}^{\tilde{\alpha}}}\left[\int_{t}^{t_{I(u)}}\frac{1}{2}\left \lVert\tilde{\alpha}_{s}\right\rVert^{2}ds-\sum_{i:\{t\leq t_{i}\leq t_{I(u)} \}}\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\tilde{\alpha}}_{t_{i}})+\mathcal {V}_{I(u)+1}(t_{I(u)},\mathbf{X}^{\tilde{\alpha}}_{t_{I(u)}})\right]$ (92) $\displaystyle\geq\mathcal{V}(t,\mathbf{x}).$ (93) $$

Combining both inequalities in (87, 92-93), we arrive at the desired result:

$$ $\mathcal{V}(t,\mathbf{x})=\min_{\alpha\in\mathbb{A}}\mathbb{E}_{\mathbb{P}^{ \alpha}}\left[\int_{t}^{t_{I(u)}}\frac{1}{2}\left\lVert\alpha_{s}\right\rVert^ {2}ds-\sum_{i:\{t\leq t_{i}\leq t_{I(u)}\}}\log f_{i}(\mathbf{y}_{t_{i}}| \mathbf{X}^{\alpha}_{t_{i}})+\mathcal{V}_{I(u)+1}(t_{I(u)},\mathbf{X}^{\alpha} _{t_{I(u)}})|\mathbf{X}^{\alpha}_{t}=\mathbf{x}_{t}\right].$ (94) $$

This concludes the proof.
∎

###### Theorem 3.2 (Dynamic Programming Principle) .

###### Proof.

### A.3 Proof of Theorem 3.3

###### Theorem 3.3 (Verification Theorem) .

Suppose there exist a sequence of left continuous functions $\mathcal{V}_{i}(t,\mathbf{x})\in C^{1,2}([t_{i-1},t_{i}),\mathbb{R}^{d})$, for all $i\in[1:k]$, satisfying the following Hamiltonian-Jacobi-Bellman (HJB) equation:

$$ $\displaystyle\partial_{t}\mathcal{V}_{i,t}+\mathcal{A}_{t}\mathcal{V}_{i,t}+ \min_{\alpha\in\mathbb{A}}\left[(\nabla_{\mathbf{x}}\mathcal{V}_{i,t})^{\top} \alpha_{i,t}+\frac{1}{2}\left\lVert\alpha_{i,t}\right\rVert^{2}\right]=0,\quad t _{i-1}\leq t<t_{i}$ (95) $\displaystyle\mathcal{V}_{i}(t_{i},\mathbf{x})=-\log f_{i}(\mathbf{y}_{t_{i}}| \mathbf{x})+\mathcal{V}_{i+1}(t_{i},\mathbf{x}),\quad t=t_{i},\quad\forall i \in[1:k],$ (96) $$

where a minimum is attained by $\alpha^{\star}_{i}(t,\mathbf{x})=\nabla_{\mathbf{x}}\mathcal{V}_{i}(t,\mathbf{
x})$. Now, define a function $\alpha:[0,t_{k}]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ by integrating the optimal controls $\{\alpha_{i}\}_{i\in\{1,\cdots,k\}}$,

$$ $\alpha^{\star}(t,\mathbf{x}):=\sum_{i=1}^{k}\alpha^{\star}_{i}(t,\mathbf{x}) \mathbf{1}_{[t_{i-1},t_{i})}(t)$ (97) $$

Then, $\mathcal{J}(t,\mathbf{x}_{t},\alpha^{\star})\leq\mathcal{J}(t,\mathbf{x}_{t},\alpha)$ holds for any $(t,\mathbf{x}_{t})\in[0,T]\times\mathbb{R}^{d}$ and $\alpha\in\mathbb{A}$. In other words, $\alpha^{\star}$ is optimal control for $\mathcal{V}$ in (9).

###### Proof.

Without loss of generality, consider $t\in[t_{i-1},t_{i})$. By applying the Itô’s formula to the value function $\mathcal{V}$ and taking expectation with respect to $\mathbb{P}^{\alpha}$, we obtain

$$ $\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha}}\left[\mathcal{V}_{i}(t_{i} ,\mathbf{X}^{\alpha}_{t_{i}})\right]=\mathcal{V}_{i}(t,\mathbf{x})+\mathbb{E}^ {t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha}}\left[\int_{t}^{t_{i}}\left(\partial_{ t}\mathcal{V}_{i,s}+\mathcal{A}_{t}\mathcal{V}_{i,s}+(\nabla_{\mathbf{x}} \mathcal{V}_{i,s})^{\top}\alpha_{i,s}\right)ds\right],$ (98) $$

where we denote $\mathbb{E}^{\mathbf{t},\mathbf{x}}_{\mathbb{P}}[\cdot]=\mathbb{E}_{\mathbb{P}}
[\cdot|\mathbf{X}_{t}=\mathbf{x}]$. By adding the Lagrangian term $\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha}}\left[\int_{t}^{t_{i}}\frac
{1}{2}\left\lVert\alpha_{i,s}\right\rVert^{2}ds\right]$ to both sides of (98), we get for the LHS of (98)

$$ LHS $\displaystyle=\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha}}\left[ \mathcal{V}_{i}(t_{i},\mathbf{X}^{\alpha}_{t_{i}})\right]+\mathbb{E}^{t, \mathbf{x}_{t}}_{\mathbb{P}^{\alpha}}\left[\int_{t}^{t_{i}}\frac{1}{2}\left \lVert\alpha_{i,s}\right\rVert^{2}ds\right]$ (99) $\displaystyle=\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha}}\left[ \mathcal{V}_{i}(t_{i},\mathbf{X}^{\alpha}_{t_{i}})+\int_{t}^{t_{i}}\frac{1}{2} \left\lVert\alpha_{i,s}\right\rVert^{2}ds\right]$ (100) $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\mathbb{E}^{t,\mathbf{x}_{t}}_{ \mathbb{P}^{\alpha}}\left[\int_{t}^{t_{i}}\frac{1}{2}\left\lVert\alpha_{i,s} \right\rVert^{2}ds-\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})+ \mathcal{V}_{i+1}(t_{i},\mathbf{X}^{\alpha}_{t_{i}})\right]$ (101) $\displaystyle=\mathcal{J}_{i}(t,\mathbf{x},\alpha),$ (102) $$

where $(i)$ follows from the definition of HJB equation in (11). Now for the RHS of (98), we have:

$$ $\displaystyle\text{RHS}=\mathcal{V}_{i}(t,\mathbf{x})+\mathbb{E}^{t,\mathbf{x} _{t}}_{\mathbb{P}^{\alpha}}\left[\int_{t}^{t_{i}}\left(\partial_{t}\mathcal{V} _{i,s}+\mathcal{A}_{t}\mathcal{V}_{i,s}+(\nabla_{\mathbf{x}}\mathcal{V}_{i,s}) ^{\top}\alpha_{i,s}\right)ds\right]+\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^ {\alpha}}\left[\int_{t}^{t_{i}}\frac{1}{2}\left\lVert\alpha_{i,s}\right\rVert^ {2}ds\right]$ (103) $\displaystyle=\mathcal{V}_{i}(t,\mathbf{x})+\mathbb{E}^{t,\mathbf{x}_{t}}_{ \mathbb{P}^{\alpha}}\left[\int_{t}^{t_{i}}\left(\partial_{t}\mathcal{V}_{i,s}+ \mathcal{A}_{t}\mathcal{V}_{i,s}+\left[(\nabla_{\mathbf{x}}\mathcal{V}_{i,s})^ {\top}\alpha_{i,s}+\frac{1}{2}\left\lVert\alpha_{i,s}\right\rVert^{2}\right] \right)ds\right]$ (104) $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\mathcal{V}_{i}(t,\mathbf{x})+ \mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^{\alpha}}\left[\int_{t}^{t_{i}}\left (\left[(\nabla_{\mathbf{x}}\mathcal{V}_{i,s})^{\top}\alpha_{i,s}+\frac{1}{2} \left\lVert\alpha_{i,s}\right\rVert^{2}\right]-\min_{\alpha\in\mathbb{A}}\left [(\nabla_{\mathbf{x}}\mathcal{V}_{i,s})^{\top}\alpha_{i,s}+\frac{1}{2}\left \lVert\alpha_{i,s}\right\rVert^{2}\right]\right)ds\right],$ (105) $$

where $(i)$ follows from the definition of HJB equation in (10).
Therefore, we get the following result:

$$ $\displaystyle\mathcal{J}_{i}(t,\mathbf{x},\alpha)=\mathcal{V}_{i}(t,\mathbf{x})$ $\displaystyle\hskip 25.60747pt+\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb{P}^{ \alpha}}\left[\int_{t}^{t_{i}}\left(\left[(\nabla_{\mathbf{x}}\mathcal{V}_{i,s })^{\top}\alpha_{i,s}+\frac{1}{2}\left\lVert\alpha_{i,s}\right\rVert^{2}\right ]-\min_{\alpha\in\mathbb{A}}\left[(\nabla_{\mathbf{x}}\mathcal{V}_{i,s})^{\top }\alpha_{i,s}+\frac{1}{2}\left\lVert\alpha_{i,s}\right\rVert^{2}\right]\right) ds\right].$ (106) $$

Due to the fact that

$$ $\int_{t}^{t_{i}}\left(\left[(\nabla_{\mathbf{x}}\mathcal{V}_{i,s})^{\top} \alpha_{i,s}+\frac{1}{2}\left\lVert\alpha_{i,s}\right\rVert^{2}\right]-\min_{ \alpha\in\mathbb{A}}\left[(\nabla_{\mathbf{x}}\mathcal{V}_{i,s})^{\top}\alpha_ {i,s}+\frac{1}{2}\left\lVert\alpha_{i,s}\right\rVert^{2}\right]\right)ds\geq 0,$ $$

for all $t\in[t_{i-1},t_{i})$ and $\mathbb{P}^{\alpha}$ almost $\mathbf{x}$, we conclude that $\mathcal{J}_{i}(t,\mathbf{x},\alpha)\geq\mathcal{V}_{i}(t,\mathbf{x})$, where the equality holds for $\alpha^{\star}_{i,t}=\min_{\alpha\in\mathbb{A}}\left[(\nabla_{\mathbf{x}}
\mathcal{V}_{i,t})^{\top}\alpha_{i,t}+\frac{1}{2}\left\lVert\alpha_{i,t}\right
\rVert^{2}\right]=-\nabla_{\mathbf{x}}\mathcal{V}_{i,t}$. Additionaly, it implies that $\mathcal{J}_{i}(t,\mathbf{x},\alpha^{\star})=\mathcal{V}_{i}(t,\mathbf{x})$. Subsequently, for any $t\in[t_{i-1},t_{i})$, the recursion in (9) yields:

$$ $\displaystyle\mathcal{V}(t,\mathbf{x}_{t})$ $\displaystyle=\min_{\alpha\in\mathbb{A}}\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb {P}^{\alpha}}\left[\int_{t}^{t_{I(u)}}\frac{1}{2}\left\lVert\alpha_{s}\right \rVert^{2}ds-\sum_{i:\{t\leq t_{i}\leq t_{I(u)}\}}\log f_{i}(\mathbf{y}_{t_{i} }|\mathbf{X}^{\alpha}_{t_{i}})+\mathcal{V}_{I(u)+1}(t_{I(u)},\mathbf{X}^{ \alpha}_{t_{I(u)}})\right]$ (107) $\displaystyle=\min_{\alpha\in\mathbb{A}}\mathbb{E}^{t,\mathbf{x}_{t}}_{\mathbb {P}^{\alpha}}\left[\int_{t}^{t_{i}}\frac{1}{2}\left\lVert\alpha_{i,s}\right \rVert^{2}ds-\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})+ \mathcal{V}_{i+1}(t_{i},\mathbf{X}^{\alpha}_{t_{i}})\right]$ (108) $\displaystyle=\mathcal{V}_{i}(t,\mathbf{x}_{t}).$ (109) $$

This implies that the optimal control for the value function $\mathcal{V}$ in (9) over the interval $t\in[t_{i-1},t_{i})$ is $\alpha^{\star}_{i}$.
Finally, $\mathcal{V}$ in (9) can be represented as the integrated form :

$$ $\mathcal{V}(t,\mathbf{x}_{t})=\sum_{i=1}^{k}\mathcal{V}_{i}(t,\mathbf{x}_{t}) \mathbf{1}_{[t_{i-1},t_{i})}.$ (110) $$

Therefore, the optimal control for $\mathcal{V}$ in (9) becomes $\alpha(t,\mathbf{x})=\sum_{i=1}^{k}\alpha^{\star}_{i}(t,\mathbf{x})\mathbf{1}_
{[t_{i-1},t_{i})}(t)$.
∎

###### Theorem 3.3 (Verification Theorem) .

###### Proof.

### A.4 Proof of Lemma 3.4

We first restate the Feynman-Kac formula (Oksendal, ; Baldi, ), which gives a probabilistic representation of the solution to certain PDEs using expectations of stochastic processes. It relies on the fact that the conditional expectation is martingale process.

###### Lemma A.5 (The Feynman-Kac formula) .

Let us define $f\in C^{2}(\mathbb{R}^{d})$ and $g\in C(\mathbb{R}^{d})$. Then, a function $h(t,\mathbf{x}_{t})=\mathbb{E}_{\mathbb{P}}\left[e^{-\int_{t}^{T}f(s,\mathbf{X
}_{s})ds}g(\mathbf{X}_{T})|\mathbf{X}_{t}=\mathbf{x}_{t}\right]$ is a solution of the following linear PDE:

$$ $\displaystyle\partial_{t}h_{t}+\mathcal{A}_{t}h_{t}-fh_{t}=0,\quad 0\leq t<T,$ (111) $\displaystyle h(t,\mathbf{x})=g(\mathbf{X}_{T}),\quad t=T.$ (112) $$

###### Proof.

Define the process $\mathbf{Y}_{t}=e^{-\int_{t}^{T}f(s,\mathbf{X}_{s})ds}h(t,\mathbf{X}_{t})$. Since $h$ is a conditional expectation with respect to $\mathbb{P}$, implies that $\mathbf{Y}_{t}$ is martingale process. By applinyg Itô formula, we have:

$$ $d\mathbf{Y}_{t}=-f(t,\mathbf{X}_{t})e^{-\int_{t}^{T}f(s,\mathbf{X}_{s})ds}h(t, \mathbf{X}_{t})dt+e^{-\int_{t}^{T}f(s,\mathbf{X}_{s})ds}dh(t,\mathbf{X}_{t}).$ (113) $$

Next, we apply Itô formula to $h(t,\mathbf{X}_{t})$:

$$ $dh(t,\mathbf{X}_{t})=\left(\frac{\partial h}{\partial t}+\mathcal{A}_{t}h \right)dt+\nabla_{\mathbf{x}}h(t,\mathbf{X}_{t})^{\top}d\mathbf{W}_{t},$ (114) $$

Thus, by substituting equation (114) into equation (113), we get

$$ $d\mathbf{Y}_{t}=-f(t,\mathbf{X}_{t})Y_{t}dt+e^{-\int_{t}^{T}f(s,\mathbf{X}_{s} )ds}\left(\frac{\partial h}{\partial t}+\mathcal{A}_{t}h\right)dt+e^{-\int_{t} ^{T}f(s,\mathbf{X}_{s})ds}\nabla_{\mathbf{x}}h(t,\mathbf{X}_{t})^{\top}d \mathbf{W}_{t}.$ (115) $$

For $\mathbf{Y}_{t}$ to be a martingale process, it will have zero drift. Therefore, it implies that

$$ $\frac{\partial h}{\partial t}(t,\mathbf{X}_{t})+\mathcal{A}_{t}h(t,\mathbf{X}_ {t})-f(t,\mathbf{X}_{t})h(t,\mathbf{X}_{t})=0,$ (116) $$

where $h(T,\mathbf{X}_{T})=g(\mathbf{X}_{T})$ by definition. It concludes the proof.
∎

Now let us provide the proof of the Lemma 3.4.

###### Lemma 3.4 (Hopf-Cole Transformation) .

The $h$ function satisfying the following linear PDE:

$$ $\displaystyle\partial_{t}h_{i,t}+\mathcal{A}_{t}h_{i,t}=0,\quad t_{i-1}\leq t< t_{i}$ (117) $\displaystyle h_{i}(t_{i},\mathbf{x})=f_{i}(\mathbf{y}_{t_{i}}|\mathbf{x})h_{i +1}(t_{i},\mathbf{x}),\quad t=t_{i},\quad\forall i\in[1:k].$ (118) $$

Moreover, for a logarithm transformation $\mathcal{V}=-\log h$, $\mathcal{V}$ satisfying the HJB equation in (10-11).

###### Proof.

The linaer PDE presented in (117-118) can be directly derived from the function $h_{i}(t,\mathbf{x}_{t})=\mathbb{E}_{\mathbb{P}}\left[\prod_{j\geq i-1}^{k}f_{j
}(\mathbf{X}_{t_{j}})|\mathbf{X}_{t}=\mathbf{x}_{t}\right]$. This is done by applying the Feynman-Kac formula in Lemma A.5 with $f:=0$ and $g=\prod_{j=i-1}^{k}f_{j}(\mathbf{X}_{t_{j}})$. Now, let us consider the function $\mathcal{V}_{i}(t,\mathbf{x})=-\log h_{i}(t,\mathbf{x})$ (or $h_{i}(t,\mathbf{x})=e^{-\mathcal{V}_{i}(t,\mathbf{x})}$ for all $i\in[1:k]$ and compute

$$ $\partial_{t}h_{i,t}=-h_{i,t}\partial_{t}\mathcal{V}_{i,t},\;\;\nabla_{\mathbf{ x}}h_{i,t}=-h_{i,t}\nabla_{\mathbf{x}}\mathcal{V}_{i,t},\;\;\nabla_{\mathbf{xx }}h_{i,t}=h_{i,t}(\left\lVert\nabla_{\mathbf{x}}\mathcal{V}_{i,t}\right\rVert^ {2}-\nabla_{\mathbf{xx}}\mathcal{V}_{i,t}).$ (119) $$

Then, it is straightforward to compute:

$$ $\displaystyle h_{i,t}\partial_{t}\mathcal{V}_{i,t}$ $\displaystyle=-\partial_{t}h_{i,t}\stackrel{{\scriptstyle(i)}}{{=}}\mathcal{A} _{t}h_{i,t}$ (120) $\displaystyle=(\nabla_{\mathbf{x}}h_{i,t})^{\top}b_{t}+\frac{1}{2}\textit{ Trace}\left[\nabla_{\mathbf{xx}}h_{i,t}\right]$ (121) $\displaystyle=(-h_{i,t}\nabla_{\mathbf{x}}\mathcal{V}_{i,t})^{\top}b_{t}+\frac {1}{2}\textit{Trace}\left[h_{i,t}(\left\lVert\nabla_{\mathbf{x}}\mathcal{V}_{i ,t}\right\rVert^{2}-h_{i,t}\nabla_{\mathbf{xx}}\mathcal{V}_{i,t})\right]$ (122) $\displaystyle=(-h_{i,t}\nabla_{\mathbf{x}}\mathcal{V}_{i,t})^{\top}b_{t}+\frac {1}{2}\textit{Trace}\left[h_{i,t}\left\lVert\nabla_{\mathbf{x}}\mathcal{V}_{i, t}\right\rVert^{2}\right]-\frac{1}{2}\textit{Trace}\left[h_{i,t}\nabla_{ \mathbf{xx}}\mathcal{V}_{i,t}\right],$ (123) $$

where $(i)$ follows from (117). Now, we can simplify (123) by dividing both sides with $h>0$:

$$ $\displaystyle\partial_{t}\mathcal{V}_{i,t}$ $\displaystyle=(-\nabla_{\mathbf{x}}\mathcal{V}_{i,t})^{\top}b_{t}+\frac{1}{2} \textit{Trace}\left[\left\lVert\nabla_{\mathbf{x}}\mathcal{V}_{i,t}\right \rVert^{2}\right]-\frac{1}{2}\textit{Trace}\left[\nabla_{\mathbf{xx}}\mathcal{ V}_{i,t}\right]$ (124) $\displaystyle=-\mathcal{A}_{t}\mathcal{V}_{i,t}+\frac{1}{2}\left\lVert\nabla_{ \mathbf{x}}\mathcal{V}_{i,t}\right\rVert^{2}$ (125) $$

Therefore, combining the above results, we have

$$ $\partial_{t}\mathcal{V}_{i,t}+\mathcal{A}_{t}\mathcal{V}_{i,t}-\frac{1}{2} \left\lVert\nabla_{\mathbf{x}}\mathcal{V}_{i,t}\right\rVert^{2}=0,\quad \mathcal{V}_{i}(t_{i},\mathbf{x})=-\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{x})+ \mathcal{V}_{i+1}(t_{i},\mathbf{x}).$ (126) $$

Since $\min_{\alpha\in\mathbb{A}}\left[(\nabla_{\mathbf{x}}\mathcal{V}_{i,t})^{\top}
\alpha_{i,t}+\frac{1}{2}\left\lVert\alpha_{i,t}\right\rVert^{2}\right]=-\frac{
1}{2}\left\lVert\nabla_{\mathbf{x}}\mathcal{V}_{i,t}\right\rVert^{2}$, this concludes the proof.
∎

###### Lemma A.5 (The Feynman-Kac formula) .

###### Proof.

###### Lemma 3.4 (Hopf-Cole Transformation) .

###### Proof.

### A.5 Proof of Corollary 3.5

###### Corollary 3.5 (Optimal Control) .

For optimal control $\alpha^{\star}$ induced by the cost function (7) with dynamics (6), it satisfies $\alpha^{\star}=\nabla_{\mathbf{x}}\log h$. In other words, we can simulate the conditional SDEs in (5) by simulating the controlled SDE (6) with optimal control $\alpha^{\star}$.

###### Proof.

Lemma 3.4 implies that we can obtain the relation $-\mathcal{V}_{i}(t,\mathbf{x})=\log h_{i}(t,\mathbf{x})$. Moreover, by following the definition of the optimal control $\alpha^{\star}$ in (12) and the value function $\mathcal{V}$ in (110), it suggest that $\alpha^{\star}(t,\mathbf{x})=-\nabla_{\mathbf{x}}\mathcal{V}(t,\mathbf{x})$ for all $t\in[0,T]$. Finally, combining the definition of the $h$ function in (4) and the results from Lemma 3.4, we can conclude that $\alpha^{\star}(t,\mathbf{x})=-\nabla_{\mathbf{x}}\mathcal{V}(t,\mathbf{x})=
\nabla_{\mathbf{x}}\log h(t,\mathbf{x})$ for all $t\in[0,T]$.
∎

###### Corollary 3.5 (Optimal Control) .

###### Proof.

### A.6 Proof of Theorem 3.6

The Donsker–Varadhan variational principle provides a variational formula for the large deviations of functionals of Brownian motion, often related to free-energy minimization problems (Boué & Dupuis, ). Moreover, through Girsanov’s theorem, this principle extends to a wide range of Markov processes, including Itô diffusion processes (Hartmann et al., ; Tzen & Raginsky, ).

###### Lemma A.6 (Donsker–Varadhan Variational Principle) .

For a bounded and measurable functions $\mathcal{W}:C([0,T],\mathbb{R}^{d})\to\mathbb{R}$, following relation holds:

$$ $-\log\mathbb{E}_{\mathbf{X}\sim\mathbb{P}}\left[e^{-\mathcal{W}(\mathbf{X}_{0: T})}\right]=\min_{\mathbb{Q}\ll\mathbb{P}}\left[\mathbb{E}_{\mathbf{Y}\sim \mathbb{Q}}\left[\mathcal{W}(\mathbf{Y}_{[0,T]})\right]+D_{\text{KL}}(\mathbb{ Q}|\mathbb{P})\right]$ (127) $$

###### Proof.

The proof relies on the change of measure and the Jensen’s inequality:

$$ $\displaystyle-\log\mathbb{E}_{\mathbf{X}\sim\mathbb{P}}\left[e^{-\mathcal{W}( \mathbf{X}_{0:T})}\right]$ $\displaystyle=-\log\mathbb{E}_{\mathbf{Y}\sim\mathbb{Q}}\left[e^{-\mathcal{W}( \mathbf{Y}_{0:T})}\frac{d\mathbb{P}}{d\mathbb{Q}}(\mathbf{Y}_{0:T})\right]$ (128) $\displaystyle\leq\mathbb{E}_{\mathbf{Y}\sim\mathbb{Q}}\left[\mathcal{W}( \mathbf{Y}_{0:T})-\log\frac{d\mathbb{P}}{d\mathbb{Q}}(\mathbf{Y}_{0:T})\right]$ (129) $\displaystyle=\mathbb{E}_{\mathbf{Y}\sim\mathbb{Q}}\left[\mathcal{W}(\mathbf{Y }_{0:T})\right]+D_{\text{KL}}(\mathbb{Q}|\mathbb{P}),$ (130) $$

where the equality holds if and only if

$$ $\frac{d\mathbb{Q}}{d\mathbb{P}}(\mathbf{Y}_{0:T})=e^{-\log\mathbb{E}_{\mathbf{ X}\sim\mathbb{P}}\left[e^{-\mathcal{W}(\mathbf{X}_{0:T})}\right]-\mathcal{W}( \mathbf{Y}_{0:T})}.$ (131) $$

It concludes the proof.
∎

###### Theorem 3.6 (Tight Variational Bound) .

Let us assume that the path measure $\mathbb{P}^{\alpha}$ induced by (6) for any $\alpha\in\mathbb{A}$ satisfies $D_{\text{KL}}(\mathbb{P}^{\alpha}|\mathbb{P}^{\star})<\infty$. Then, for a cost function $\mathcal{J}$ in (7) and $\mu^{\star}_{0}$ in (5), it holds:

$$ $\displaystyle D_{\text{KL}}(\mathbb{P}^{\alpha}|\mathbb{P}^{\star})=D_{\text{ KL}}(\mu_{0}|\mu^{\star}_{0})+\mathbb{E}_{\mathbf{x}_{0}\sim\mu_{0}}\left[ \mathcal{J}(0,\mathbf{x}_{0},\alpha)\right]=\mathcal{L}(\alpha)+\log g( \mathcal{H}_{t_{k}})\geq 0,$ (132) $$

where the objective function $\mathcal{L}(\alpha)$ (negative ELBO) is given by:

$$ $\mathcal{L}(\alpha)=\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\int_{0}^{T}\frac{1} {2}\left\lVert\alpha(s,\mathbf{X}^{\alpha}_{s})\right\rVert^{2}ds-\sum_{i}^{k} \log g_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})\right]\geq-\log g( \mathcal{H}_{t_{k}}).$ (133) $$

Moreover, assume that $\mathcal{L}(\alpha)$ has a global minimum $\alpha^{\star}=\operatorname*{arg\,min}_{\alpha\in\mathbb{A}}\mathcal{L}(\alpha)$. Then the equality holds in (16) $\textit{i}.\textit{e}.,\mathcal{L}(\alpha^{\star})=-\log g(\mathcal{H}_{t_{k}})$ and $\mu^{\star}_{0}=\mu_{0}$ almost everywhere with respect to $\mu_{0}$.

###### Proof.

We begin by deriving the KL-divergence between $\mathbb{P}^{\alpha}$ and $\mathbb{P}^{\star}$:

$$ $\displaystyle D_{\text{KL}}(\mathbb{P}^{\alpha}|\mathbb{P}^{\star})$ $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}D_{\text{KL}}(\mu_{0}|\mu^{\star }_{0})+\mathbb{E}_{\mathbf{x}_{0}\sim\mu_{0}}\left[D_{\text{KL}}(\mathbb{P}^{ \alpha}(\cdot|\mathbf{X}^{\alpha}_{0})|\mathbb{P}^{\star}(\cdot|\mathbf{X}^{ \alpha}_{0}))|\mathbf{X}^{\star}_{0}=\mathbf{x}_{0}\right]$ (134) $\displaystyle=$ $\displaystyle D_{\text{KL}}(\mu_{0}|\mu^{\star}_{0})+\mathbb{E}_{\mathbf{x}_{0 }\sim\mu_{0}}\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\log\frac{d\mathbb{P}^{ \alpha}}{d\mathbb{P}^{\star}}(\mathbf{X}^{\alpha}_{0:T})|\mathbf{X}^{\alpha}_{ 0}=\mathbf{x}_{0}\right]$ (135) $\displaystyle=$ $\displaystyle D_{\text{KL}}(\mu_{0}|\mu^{\star}_{0})+\mathbb{E}_{\mathbf{x}_{0 }\sim\mu_{0}}\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\log\frac{d\mathbb{P}^{ \alpha}}{d\mathbb{P}}(\mathbf{X}^{\alpha}_{0:T})+\log\frac{d\mathbb{P}}{d \mathbb{P}^{\star}}(\mathbf{X}^{\alpha}_{0:T})|\mathbf{X}^{\alpha}_{0}=\mathbf {x}_{0}\right],$ (136) $$

where $(i)$ follows from the disintegration theorem (Léonard, ). For the second term of the RHS of (136), applying the Girsanov’s theorem in A.4, we have:

$$ $\displaystyle\mathbb{E}^{0,\mathbf{x}_{0}}_{\mathbb{P}^{\alpha}}\left[\log \frac{d\mathbb{P}^{\alpha}}{d\mathbb{P}}(\mathbf{X}^{\alpha}_{0:T})\right]= \mathbb{E}^{0,\mathbf{x}_{0}}_{\mathbb{P}^{\alpha}}\left[\int_{0}^{T}\alpha_{t }d\tilde{\mathbf{W}}_{t}+\int_{0}^{T}\frac{1}{2}\left\lVert\alpha_{t}\right \rVert^{2}dt\right]=\mathbb{E}^{0,\mathbf{x}_{0}}_{\mathbb{P}^{\alpha}}\left[ \int_{0}^{T}\frac{1}{2}\left\lVert\alpha_{t}\right\rVert^{2}dt\right].$ (137) $$

Moreover, by the definition of $\mathbb{P}^{\star}$, we get:

$$ $\mathbb{E}^{0,\mathbf{x}_{0}}_{\mathbb{P}^{\alpha}}\left[\log\frac{d\mathbb{P} }{d\mathbb{P}^{\star}}(\mathbf{X}^{\alpha}_{0:T})\right]=\mathbb{E}^{0,\mathbf {x}_{0}}_{\mathbb{P}^{\alpha}}\left[-\sum_{i=1}^{k}\log f_{i}(\mathbf{y}_{t_{i }}|\mathbf{X}^{\alpha}_{t_{i}})\right].$ (138) $$

Combining the results from (137) and (138), we obtain:

$$ $\displaystyle\mathbb{E}^{0,\mathbf{x}_{0}}_{\mathbb{P}^{\alpha}}\left[\log \frac{d\mathbb{P}^{\alpha}}{d\mathbb{P}}(\mathbf{X}^{\alpha}_{0:T})+\log\frac{ d\mathbb{P}}{d\mathbb{P}^{\star}}(\mathbf{X}^{\alpha}_{0:T})\right]$ $\displaystyle=\mathbb{E}^{0,\mathbf{x}_{0}}_{\mathbb{P}^{\alpha}}\left[\int_{0 }^{T}\frac{1}{2}\left\lVert\alpha(t,\mathbf{X}^{\alpha}_{t})\right\rVert^{2}dt -\sum_{i=1}^{k}\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})\right]$ (139) $\displaystyle=\mathbb{E}^{0,\mathbf{x}_{0}}_{\mathbb{P}^{\alpha}}\left[\int_{0 }^{T}\frac{1}{2}\left\lVert\alpha(t,\mathbf{X}^{\alpha}_{t})\right\rVert^{2}dt -\sum_{i=1}^{k}\log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})\right]$ (140) $\displaystyle=\mathcal{J}(0,\mathbf{x}_{0},\alpha).$ (141) $$

Moreover, by the definition of $\{f_{i}\}_{i\in[1:k]}$, we get

$$ $\displaystyle\mathcal{J}(0,\mathbf{x}_{0},\alpha)$ $\displaystyle=\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\int_{0}^{T}\frac{1}{2} \left\lVert\alpha(t,\mathbf{X}^{\alpha}_{t})\right\rVert^{2}dt-\sum_{i=1}^{k} \log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}_{t_{i}})|\mathbf{X}_{0}^{\alpha}= \mathbf{x}_{0}\right]$ (142) $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\mathbb{E}_{\mathbb{P}^{\alpha}} \left[\int_{0}^{T}\frac{1}{2}\left\lVert\alpha(t,\mathbf{X}^{\alpha}_{t}) \right\rVert^{2}dt-\sum_{i=1}^{k}\log g_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}_{t_{ i}})|\mathbf{X}_{0}^{\alpha}=\mathbf{x}_{0}\right]+\log\mathbf{Z}(\mathcal{H}_ {t_{k}}),$ (143) $$

where $(i)$ follows from normalizing property in (39). Hence, it result the that:

$$ $\displaystyle D_{\text{KL}}(\mathbb{P}^{\alpha}|\mathbb{P}^{\star})$ $\displaystyle=D_{\text{KL}}(\mu_{0}|\mu^{\star}_{0})+\mathbb{E}_{\mathbf{x}_{0 }\sim\mu_{0}}\left[\mathcal{J}(0,\mathbf{x}_{0},\alpha)\right]$ (144) $\displaystyle=D_{\text{KL}}(\mu_{0}|\mu^{\star}_{0})+\mathcal{L}(\alpha)+\log \mathbf{Z}(\mathcal{H}_{t_{k}}).$ (145) $$

For a KL-divergence term $D_{\text{KL}}(\mu_{0}|\mu^{\star}_{0})$, since $d\mu^{\star}_{0}(\mathbf{x}_{0})=h_{1}(t_{0},\mathbf{x}_{0})d\mu_{0}(\mathbf{x
}_{0})$, we obtain

$$ $\displaystyle D_{\text{KL}}(\mu_{0}|\mu^{\star}_{0})$ $\displaystyle=\mathbb{E}_{\mathbf{x}_{0}\sim\mu_{0}}\left[\log\frac{d\mu_{0}}{ d\mu^{\star}_{0}}(\mathbf{x}_{0})\right]=\mathbb{E}_{\mathbf{x}_{0}\sim\mu_{0} }\left[-\log h_{1}(0,\mathbf{x}_{0})\right]$ (146) $\displaystyle=\mathbb{E}_{\mathbf{x}_{0}\sim\mu_{0}}\left[\mathcal{V}(0, \mathbf{x}_{0})\right]$ (147) $\displaystyle=\mathbb{E}_{\mathbf{x}_{0}\sim\mu_{0}}\left[\mathcal{J}(0, \mathbf{x}_{0},\alpha^{\star})\right]$ (148) $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\mathbb{E}_{\mathbf{x}_{0}\sim \mu_{0}}\left[\tilde{\mathcal{J}}(0,\mathbf{x}_{0},\alpha^{\star})\right]+\log \mathbf{Z}(\mathcal{H}_{t_{k}})$ (149) $\displaystyle\stackrel{{\scriptstyle(ii)}}{{=}}\mathcal{L}(\alpha^{\star})+ \log\mathbf{Z}(\mathcal{H}_{t_{k}}),$ (150) $$

where $(i)$ follows from:

$$ $\displaystyle\mathcal{J}(0,\mathbf{x}_{0},\alpha)$ $\displaystyle=\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\int_{0}^{T}\frac{1}{2} \left\lVert\alpha(s,\mathbf{X}^{\alpha}_{s})\right\rVert^{2}ds-\sum_{i}^{k} \log f_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})|\mathbf{X}^{\alpha} _{0}=\mathbf{x}_{0}\right]$ (151) $\displaystyle=\underbrace{\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\int_{0}^{T} \frac{1}{2}\left\lVert\alpha(s,\mathbf{X}^{\alpha}_{s})\right\rVert^{2}ds-\sum _{i}^{k}\log g_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})|\mathbf{X}^ {\alpha}_{0}=\mathbf{x}_{0}\right]}_{\tilde{\mathcal{J}}(0,\mathbf{x}_{0}, \alpha)}+\log\mathbf{Z}(\mathcal{H}_{t_{k}}).$ (152) $$

Thus, we obtain $\alpha^{\star}=\operatorname*{arg\,min}_{\alpha\in\mathbb{A}}\mathcal{J}(0,
\mathbf{x}_{0},\alpha)=\operatorname*{arg\,min}_{\alpha\in\mathbb{A}}\tilde{
\mathcal{J}}(0,\mathbf{x}_{0},\alpha)$ as $\log\mathbf{Z}(\mathcal{H}_{t_{k}})$ is constant. It result that $\mathcal{J}(0,\mathbf{x}_{0},\alpha^{\star})=\tilde{\mathcal{J}}(0,\mathbf{x}_
{0},\alpha^{\star})+\log\mathbf{Z}(\mathcal{H}_{t_{k}})$. Additionally, $(ii)$ follows from

$$ $\mathbb{E}_{\mathbf{x}_{0}\sim\mu_{0}}\left[\min_{\alpha\in\mathbb{A}}\tilde{ \mathcal{J}}(0,\mathbf{x}_{0},\alpha)\right]=\min_{\alpha\in\mathbb{A}}\mathbb {E}_{\mathbf{x}_{0}\sim\mu_{0}}\left[\tilde{\mathcal{J}}(0,\mathbf{x}_{0}, \alpha)\right]=\min_{\alpha\in\mathbb{A}}\mathcal{L}(\alpha)=\mathcal{L}( \alpha^{\star})$ (153) $$

since the minimization is independent of the initial condition, as implied by the disintegration theorem (Léonard, ). Hence, since the Donsker–Varadhan variational principle in Lemma A.6 with $\mathcal{W}(\mathbf{X}^{\alpha}_{0:T})=-\sum_{i=1}^{k}\log g_{i}(\mathbf{y}_{t
_{i}}|\mathbf{X}^{\alpha}_{t_{i}})$ yields

$$ $\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\int_{0}^{T}\frac{1}{2}\left\lVert\alpha (t,\mathbf{X}^{\alpha}_{t})\right\rVert^{2}dt-\sum_{i=1}^{k}\log g_{i}(\mathbf {y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}})\right]+\log\mathbb{E}_{\mathbb{P}} \left[\prod_{i=1}^{k}g_{i}(\mathbf{y}_{t_{i}}|\mathbf{X}^{\alpha}_{t_{i}}) \right]\geq 0,$ (154) $$

it implies that $D_{\text{KL}}(\mu_{0}|\mu^{\star}_{0})=\mathcal{L}(\alpha^{\star})+\log\mathbf
{Z}(\mathcal{H}_{t_{k}})=0$ for the optimal control $\alpha^{\star}$. In other words, $\mu_{0}^{\star}=\mu_{0}$ almost everywhere with respect to $\mu_{0}$ and from the variational bound

$$ $\displaystyle D_{\text{KL}}(\mathbb{P}^{\alpha}|\mathbb{P}^{\star})=\mathcal{L }(\alpha)+\log\mathbf{Z}(\mathcal{H}_{t_{k}})\geq 0,$ (155) $$

the equality holds if and only if $\alpha\to\alpha^{\star}$. It concludes the proof.
∎

###### Lemma A.6 (Donsker–Varadhan Variational Principle) .

###### Proof.

###### Theorem 3.6 (Tight Variational Bound) .

###### Proof.

### A.7 Derivation of Amortized ELBO in ( 28 ).

Let $\mathbf{o}_{0:T}$ is given time-series data. Then, for an auxiliary variable $\mathbf{y}_{0:T}\sim q_{\phi}(\mathbf{y}_{0:T}|\mathbf{o}_{0:T})$, the ELBO is given as

$$ $\displaystyle\log p_{\psi}(\mathbf{o}_{0:T})\geq\mathbb{E}_{q_{\phi}(\mathbf{y }_{0:T}|\mathbf{o}_{0:T})}\left[\log\frac{\prod_{i=1}^{K}p_{\psi}(\mathbf{o}_{ t_{i}}|\mathbf{y}_{t_{i}})g(\mathbf{y}_{0:T})}{\prod_{i=1}^{K}q_{\phi}(\mathbf {y}_{t_{i}}|\mathbf{o}_{t_{i}})}\right]$ (156) $\displaystyle=\mathbb{E}_{q_{\phi}(\mathbf{y}_{0:T}|\mathbf{o}_{0:T})}\left[ \sum_{i=1}^{K}\log p_{\psi}(\mathbf{o}_{t_{i}}|\mathbf{y}_{t_{i}})+\log g( \mathbf{y}_{0:T})-\log\prod_{i=1}^{K}q_{\phi}(\mathbf{y}_{t_{i}}|\mathbf{o}_{t _{i}})\right]$ (157) $\displaystyle\geq\mathbb{E}_{q_{\phi}(\mathbf{y}_{0:T}|\mathbf{o}_{0:T})}\left [\sum_{i=1}^{K}\log p_{\psi}(\mathbf{o}_{t_{i}}|\mathbf{y}_{t_{i}})-\mathcal{L }(\alpha)-\log\prod_{i=1}^{K}q_{\phi}(\mathbf{y}_{t_{i}}|\mathbf{o}_{t_{i}})\right]$ (158) $\displaystyle=\mathbb{E}_{q_{\phi}(\mathbf{y}_{0:T}|\mathbf{o}_{0:T})}\left[ \sum_{i=1}^{K}\log p_{\psi}(\mathbf{o}_{t_{i}}|\mathbf{y}_{t_{i}})-\mathcal{L} (\alpha)-\sum_{i=1}^{K}\log q_{\phi}(\mathbf{y}_{t_{i}}|\mathbf{o}_{t_{i}})\right]$ (159) $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\mathbb{E}_{q_{\phi}(\mathbf{y}_ {0:T}|\mathbf{o}_{0:T})}\left[\sum_{i=1}^{K}\log p_{\psi}(\mathbf{o}_{t_{i}}| \mathbf{y}_{t_{i}})-\mathcal{L}(\alpha)\right],$ (160) $$

where $(i)$ follows from $\mathbb{E}_{q_{\phi}(\mathbf{y}_{0:T}|\mathbf{o}_{0:T})}\left[-\sum_{i=1}^{K}
\log q_{\phi}(\mathbf{y}_{t_{i}}|\mathbf{o}_{t_{i}})\right]=C$ since $q_{\phi}$ is Gaussian distribution with constant covariance matrix.

### A.8 Proof of Theorem 3.8

###### Theorem 3.8 (Simulation-free estimation) .

Let us consider sequences of SPD matrices $\{\mathbf{A}_{i}\}_{i\in[1:k]}$ and vectors $\{\alpha_{i}\}_{i\in[1:k]}$ and the following control-affine SDEs for all $i\in[1:k]$:

$$ $d\mathbf{X}_{t}=\left[-\mathbf{A}_{i}\mathbf{X}_{t}+\alpha_{i}\right]dt+\sigma d \mathbf{W}_{t},\quad t\in[t_{i-1},t_{i}).$ (161) $$

Then, the solution of (23) condition to the initial distribution $\mu_{0}=\mathcal{N}(\mathbf{m}_{0},\mathbf{\Sigma}_{0})$ is a Gaussian process $\mathcal{N}(\mathbf{m}_{t},\mathbf{\Sigma}_{t})$ for all $t\in[0,T]$, where the first two moments is given by

$$ $\displaystyle\mathbf{m}_{t}=e^{-(t-t_{0})\mathbf{A}_{i}}\mathbf{m}_{0}-\sum_{j =1}^{i-1}\left(e^{-(t-t_{k})\mathbf{A}_{j}}\mathbf{A}_{j}^{-1}(e^{-(t_{j}-t_{k -1})\mathbf{A}_{j}}-\mathbf{I})\alpha_{j}\right)-\mathbf{A}_{i}^{-1}(-e^{(t-t_ {i-1})\mathbf{A}_{i}}-\mathbf{I})\alpha_{i},$ $\displaystyle\mathbf{\Sigma}_{t}=e^{-2(t-t_{0})\mathbf{A}_{i}}\mathbf{\Sigma}_ {0}-\sum_{j=1}^{i-1}\left(e^{-2(t-t_{j})\mathbf{A}_{j}}\mathbf{A}_{j}^{-1} \frac{(e^{-2(t_{j}-t_{j-1})\mathbf{A}_{j}}-\mathbf{I})}{2}\right)-\mathbf{A}_{ i}^{-1}\frac{(e^{-2(t-t_{i-1})\mathbf{A}_{i}}-\mathbf{I})}{2}.$ $$

###### Proof.

Note that since $\{\mathbf{A}_{i}\}_{i\in[1:k]}$ are SPD matrices, we can apply the transformation outlined in Remark 3.7. For simplicity, we will treat $\{\mathbf{A}_{i}\}_{i\in[1:k]}$ as diagonal matrices, and express the dynamics in the projected form using the eigenbasis as it can be effortlessly transforming back to standard basis. Now, for any $t\in[t_{i},t_{i+1})$, the solution to (23) at time $t$ is given as

$$ $\mathbf{X}_{t}=e^{-\Delta_{i}(t)\mathbf{A}_{i}}\left(\mathbf{X}_{t_{i}}+\int_{ t_{i}}^{t}e^{\Delta_{i}(s)\mathbf{A}_{i}}\alpha_{i}ds+\int_{t_{i}}^{t}e^{ \Delta_{i}(s)\mathbf{A}_{i}}d\mathbf{W}_{s}\right),\Delta_{i}(t)=\begin{cases} t-t_{i},\;\text{for }t>t_{i}\\ 0,\;\text{for }t\leq t_{i}.\end{cases}$ (162) $$

Given that we have defined $\mathbf{X}_{0}\sim\mathcal{N}(\mathbf{m}_{0},\mathbf{\Sigma}_{0})$, $\mathbf{X}_{t_{i}}$ is a Gaussian process for any $i\in[1:k]$.
The first two moments of Gaussian process can be computed from (162). First, given $\mathbf{A}_{i}$ is diagonal, the integral can be computed as $\int_{t_{i}}^{t}e^{\Delta_{i}(s)\mathbf{A}_{i}}ds=-\mathbf{A}_{i}^{-1}(\mathbf
{I}-e^{\Delta_{i}(t)\mathbf{A}_{i}})$ and $\mathbf{M}_{i}(t):=\int_{t_{i}}^{t}e^{\Delta_{i}(s)\mathbf{A}_{i}}d\mathbf{W}_
{s}$ is a martingale process with respect to $\mathbb{P}^{\alpha}$ $\textit{i}.\textit{e}.,\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\mathbf{M}_{i}(t)
\right]=0$. Hence, since $\alpha_{i}$ is time-invariant vector, the mean $\mathbb{E}_{\mathbb{P}}[\mathbf{X}_{t}]=\mathbf{m}_{t}$ can be computed as

$$ $\mathbf{m}_{t}=e^{-\Delta_{i}(t)\mathbf{A}_{i}}\mathbf{m}_{t_{i}}-\mathbf{A}_{ i}^{-1}(e^{-\Delta_{i}(t)\mathbf{A}_{i}}-\mathbf{I})\alpha_{i}.$ (163) $$

Secondly, for a covariance $\mathbb{E}_{\mathbb{P}}[(\mathbf{X}_{t}-\mathbf{m}_{t})(\mathbf{X}_{t}-\mathbf
{m}_{t})^{\top}]=\mathbf{\Sigma}_{t}$, we can compute

$$ $\displaystyle\mathbf{\Sigma}_{t}=\mathbb{E}_{\mathbb{P}^{\alpha}}\left[e^{-2 \Delta_{i}(t)\mathbf{A}_{i}}\left(\mathbf{X}_{t_{i}}-\mathbf{m}_{t_{i}}+ \mathbf{M}_{i}(t)\right)\left(\mathbf{X}_{t_{i}}-\mathbf{m}_{t_{i}}+\mathbf{M} _{i}(t)\right)^{\top}\right]$ (164) $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}e^{-2\Delta_{i}(t)\mathbf{A}_{i} }\mathbb{E}_{\mathbb{P}^{\alpha}}\left[(\mathbf{X}_{t_{i}}-\mathbf{m}_{t_{i}}) (\mathbf{X}_{t_{i}}-\mathbf{m}_{t_{i}})^{\top}+\left\lVert\mathbf{M}_{i}(t) \right\rVert_{2}^{2}\right]$ (165) $\displaystyle\stackrel{{\scriptstyle(ii)}}{{=}}e^{-2\Delta_{i}(t)\mathbf{A}_{i }}\mathbf{\Sigma}_{t_{i}}-\frac{1}{2}\mathbf{A}_{i}^{-1}(e^{-2\Delta_{i}(t) \mathbf{A}_{i}}-\mathbf{I}),$ (166) $$

where $(i)$ follows from the fact that $\mathbf{M}_{i}(t)$ is a martingale and we use Itô isometry in $(ii)$:

$$ $\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\left\lVert\mathbf{M}_{i}(t)\right\rVert _{2}^{2}\right]=\mathbb{E}_{\mathbb{P}^{\alpha}}\left[\int_{t_{i}}^{t}\left \lVert e^{\Delta_{i}(s)\mathbf{A}_{i}}\right\rVert_{2}^{2}ds\right]=-\frac{1}{ 2}\mathbf{A}_{i}^{-1}(e^{2\Delta_{i}(t)\mathbf{A}_{i}}-\mathbf{I}).$ (167) $$

Hence, we get the Gaussian law of $\mathbf{X}_{t}$ at time $t\in[t_{i},t_{i+1})$, which is given by

$$ $\mathcal{N}\left(e^{-\Delta_{i}(t)\mathbf{A}_{i}}\mathbf{m}_{t_{i}}-\mathbf{A} ^{-1}(e^{-\Delta_{i}(t)\mathbf{A}_{i}}-\mathbf{I})\alpha_{i},e^{-2\Delta_{i}(t )\mathbf{A}_{i}}\mathbf{\Sigma}_{t_{i}}-\frac{1}{2}\mathbf{A}_{i}^{-1}(e^{-2 \Delta_{i}(t)\mathbf{A}_{i}}-\mathbf{I})\right).$ (168) $$

Furthermore, given recurrence forms of mean (163) and covariance (166), the first two moments of Gaussian distribution for each time steps $t_{i}$ can be computed sequentially. For simplicity, assume that $\mathbf{A}=\mathbf{A}_{i}=\mathbf{A}_{j}$ for all $i,j\in[1:k]$, as this can be easily extended to the case where $\mathbf{A}_{i}\neq\mathbf{A}_{j}$ and we assume that $\alpha_{j\geq k}=\alpha_{k}$. Then, for a mean $\mathbf{m}_{t_{i}}$ we have,

$$ $\displaystyle\mathbf{m}_{t_{1}}$ $\displaystyle=e^{-\Delta_{0}(t_{1})\mathbf{A}}\mathbf{m}_{t_{0}}-\mathbf{A}^{- 1}(e^{-\Delta_{0}(t_{1})\mathbf{A}}-\mathbf{I})\alpha_{1}$ (169) $\displaystyle\mathbf{m}_{t_{2}}$ $\displaystyle=e^{-\Delta_{0}(t_{2})\mathbf{A}}\mathbf{m}_{t_{0}}-e^{-\Delta_{1 }(t_{2})\mathbf{A}}\mathbf{A}^{-1}(e^{\Delta_{0}(t_{1})\mathbf{A}}-\mathbf{I}) \alpha_{1}-\mathbf{A}^{-1}(e^{-\Delta_{1}(t_{2})\mathbf{A}}-\mathbf{I})\alpha_ {2}$ (170) $\displaystyle\vdots$ (171) $\displaystyle\mathbf{m}_{t_{i}}$ $\displaystyle=e^{-\Delta_{0}(t_{i})\mathbf{A}}\mathbf{m}_{t_{0}}-\sum_{m=1}^{i }\left(e^{-\Delta_{m}(t_{i})\mathbf{A}}\mathbf{A}^{-1}(e^{\Delta_{m-1}(t_{m}) \mathbf{A}}-\mathbf{I})\alpha_{m}\right)$ (172) $$

Moreover, for a covariance $\mathbf{\Sigma}_{t_{i}}$, similar calculation yields

$$ $\displaystyle\mathbf{\Sigma}_{t_{1}}$ $\displaystyle=e^{-2\Delta_{0}(t_{1})\mathbf{A}}\mathbf{\Sigma}_{t_{0}}-\frac{1 }{2}\mathbf{A}^{-1}(e^{-2\Delta_{0}(t_{1})\mathbf{A}}-\mathbf{I})$ (173) $\displaystyle\mathbf{\Sigma}_{t_{2}}$ $\displaystyle=e^{-2\Delta_{0}(t_{2})\mathbf{A}}\mathbf{\Sigma}_{t_{0}}-\frac{1 }{2}e^{-2\Delta_{1}(t_{2})\mathbf{A}}\mathbf{A}^{-1}(e^{-2\Delta_{0}(t_{1}) \mathbf{A}}-\mathbf{I})-\frac{1}{2}\mathbf{A}^{-1}(e^{-2\Delta_{1}(t_{2}) \mathbf{A}}-\mathbf{I})$ (174) $\displaystyle\vdots$ (175) $\displaystyle\mathbf{\Sigma}_{t_{i}}$ $\displaystyle=e^{-2\Delta_{0}(t_{i})\mathbf{A}}\mathbf{\Sigma}_{t_{0}}-\frac{1 }{2}\sum_{m=1}^{i}\left(e^{-2\Delta_{m}(t_{i})\mathbf{A}}\mathbf{A}^{-1}(e^{-2 \Delta_{m-1}(t_{m})\mathbf{A}}-\mathbf{I})\right)$ (176) $$

Therefore, given a quadruplet ($\mathbf{m}_{t_{0}},\mathbf{\Sigma}_{t_{0}},\mathbf{A},\{\alpha_{i}\}_{i\in[1:k]}$), the mean and covariance for any time $t\in[0,T]$ can be computed. This concludes proof.
∎

###### Theorem 3.8 (Simulation-free estimation) .

###### Proof.

## Appendix B Parallel Scan

The parallelization of scan operation so called all-prefix-sums algorithms (Blelloch, ) have been well studied. Given an associative operator $\otimes$ and a sequence of elements $[s_{t_{1}},\cdots s_{t_{K}}]$, the parallel scan algorithm computes the all-prefix-sum which returns the sequence $[s_{t_{1}},(s_{t_{1}}\otimes s_{t_{2}}),\cdots,(s_{t_{1}}\otimes s_{t_{2}}
\otimes\cdots\otimes s_{t_{K}})]$ in $\mathcal{O}(\log K)$ time. We already have verified that the first two moments $\{\mathbf{m}_{t_{i}},\mathbf{\Sigma}_{t_{i}}\}_{i\in[1:k]}$ of the controlled distributions $\{\mu^{\alpha}_{t_{i}}\}_{i\in[1:k]}$ can be estimated by the linear recurrences in (163, 166),

$$ $\displaystyle\mathbf{m}_{t_{i}}=\hat{\mathbf{A}}_{i-1}\mathbf{m}_{t_{i-1}}+ \hat{\mathbf{B}}_{i-1}\alpha_{i-1}$ (177) $\displaystyle\mathbf{\Sigma}_{t_{i}}=\hat{\mathbf{A}}_{i-1}^{2}\mathbf{\Sigma} _{t_{-1}}+\bar{\mathbf{B}}_{i-1}\mathbf{1}_{d},$ (178) $$

where, for brevity, we define $\hat{\mathbf{A}}_{i-1}=e^{-\Delta_{i-1}(t_{i})\mathbf{A}_{i-1}}$, $\hat{\mathbf{B}}_{i-1}=-\mathbf{A}_{i-1}^{-1}(e^{-\Delta_{i-1}(t_{i})\mathbf{A
}_{i-1}}-\mathbf{I})$, $\bar{\mathbf{B}}_{i-1}=-\frac{1}{2}\mathbf{A}_{i-1}^{-1}(e^{-2\Delta_{i-1}(t_{
i})\mathbf{A}_{i-1}}-\mathbf{I})$ and $\mathbf{1}_{d}=(1,\cdots,1)\in\mathbb{R}^{d}$. For a parallel scan, we will define the sequence of tuple $\{\mathbf{M}_{i}\}_{i\in[1:k]}$, such that each element is $\mathbf{M}_{i}=(\hat{\mathbf{A}}_{i-1},\hat{\mathbf{B}}_{i-1}\alpha_{i-1})$ and $\{\mathbf{S}_{i}\}_{i\in[1:k]}$, such that each element is $\mathbf{S}_{i}=(\hat{\mathbf{A}}^{2}_{i-1},\bar{\mathbf{B}}_{i-1}\mathbf{1}_{d})$ for $\{\mathbf{m}_{t_{i}}\}_{i\in[1:k]}$ and $\{\mathbf{\Sigma}_{t_{i}}\}_{i\in[1:k]}$, respectively.

Now, let us define a binary operator $\otimes$:

$$ $\displaystyle\mathbf{M}_{i}\otimes\mathbf{M}_{i+1}=(\hat{\mathbf{A}}_{i}\circ \hat{\mathbf{A}}_{i-1},\hat{\mathbf{A}}_{i}\circ\hat{\mathbf{B}}_{i-1}\alpha_{ i-1}+\hat{\mathbf{B}}_{i}\alpha_{i})$ (179) $\displaystyle\mathbf{S}_{i}\otimes\mathbf{S}_{i+1}=(\hat{\mathbf{A}}^{2}_{i} \circ\hat{\mathbf{A}}^{2}_{i-1},\hat{\mathbf{A}}^{2}_{i}\circ\bar{\mathbf{B}}_ {i-1}\mathbf{1}_{d}+\bar{\mathbf{B}}_{i}\mathbf{1}_{d})$ (180) $$

We can verify that $\otimes$ is associative operator since it satisfying:

$$ $\displaystyle(\mathbf{M}_{s+1}\otimes\mathbf{M}_{t+1})\otimes\mathbf{M}_{u-1}= (\hat{\mathbf{A}}_{t}\circ\hat{\mathbf{A}}_{s},\hat{\mathbf{A}}_{t}\circ\hat{ \mathbf{B}}_{s}\alpha_{s}+\hat{\mathbf{B}}_{t}\alpha_{t})\otimes(\hat{\mathbf{ A}}_{u},\hat{\mathbf{B}}_{u}\alpha_{u})$ (181) $\displaystyle=(\hat{\mathbf{A}}_{u}\circ\hat{\mathbf{A}}_{t}\circ\hat{\mathbf{ A}}_{s},\hat{\mathbf{A}}_{u}\circ\hat{\mathbf{A}}_{t}\circ\hat{\mathbf{B}}_{s} \alpha_{s}+\hat{\mathbf{A}}_{u}\circ\hat{\mathbf{B}}_{t}\alpha_{t}+\hat{ \mathbf{B}}_{u}\alpha_{u})$ (182) $\displaystyle\mathbf{M}_{s+1}\otimes(\mathbf{M}_{t+1}\otimes\mathbf{M}_{u+1})= (\hat{\mathbf{A}}_{s},\hat{\mathbf{B}}_{s}\alpha_{s})\otimes(\hat{\mathbf{A}}_ {u}\circ\hat{\mathbf{A}}_{t},\hat{\mathbf{A}}_{u}\circ\hat{\mathbf{B}}_{t} \alpha_{t}+\hat{\mathbf{B}}_{u}\alpha_{u})$ (183) $\displaystyle=(\hat{\mathbf{A}}_{u}\circ\hat{\mathbf{A}}_{t}\circ\hat{\mathbf{ A}}_{s},\hat{\mathbf{A}}_{u}\circ\hat{\mathbf{A}}_{t}\circ\hat{\mathbf{B}}_{s} \alpha_{s}+\hat{\mathbf{A}}_{u}\circ\hat{\mathbf{B}}_{t}\alpha_{t}+\hat{ \mathbf{B}}_{u}\alpha_{u}).$ (184) $$

Now, the operator $\otimes$ yields similar results for $\mathbf{S}_{i}$, both $\mathbf{M}_{i}$ and $\mathbf{S}_{i}$ can be computed using a parallel scan algorithm. In other words, we can access the marginal distributions for $\mathbb{P}^{\alpha}$ at each observed time stamp in a parallel way. See (Särkkä & García-Fernández, ; Smith et al., ) for a comprehensive example.

## Appendix C Implementation Details

### C.1 Datasets

#### Human Activity

The Human Activity dataset(^6^66[https://doi.org/10.24432/C57G8X](https://doi.org/10.24432/C57G8X) under CC BY 4.0) consists of time series data collected from five individuals performing different activities. Following the preprocessing steps described in (Rubanova et al., ), we obtained $6,554$ sequences, each with $211$ time points and a fixed sequence length of 50 irregularly sampled time stamps. The time range was rescaled to $[0,1]$. The classification task involves assigning each time point to one of seven categories: “walking”, “falling”, “lying”, “sitting”, “standing up”, “on all fours”, or “sitting on the ground”. The dataset was split into 4,194 sequences for training, 1,049 for validation, and 1,311 for testing.

Figure: Figure 3: Example of the pendulum sequence. (Up) The input image sequences $\{\mathbf{o}\}_{t\in\mathcal{T}}$ observed at irregular time stamps. (Down) The angular values of ${\color[rgb]{0.4140625,0.6015625,0.765625}\definecolor[named]{pgfstrokecolor}{ rgb}{0.4140625,0.6015625,0.765625}\text{sin}(\theta_{t})}$ and ${\color[rgb]{1,0.375,0.1484375}\definecolor[named]{pgfstrokecolor}{rgb}{ 1,0.375,0.1484375}\text{cos}(\theta_{t})}$ where $\theta_{t}$ represents the angle of the pendulum at time $t\in[0,100]$, are used as regression targets.
Refer to caption: x3.png

#### Pendulum

The pendulum images were algorithmically generated through numerical simulation as outlined in (Becker et al., ). We followed the setup described in (Schirmer et al., ), where 4,000 image sequences were generated. Each sequence consists of 50 time stamps, irregularly sampled from $T=100$, with each image being a $24\times 24$ pixel representation. The sequences was further corrupted by a correlated noise process, as detailed in (Becker et al., ). For our experiments, we used 2,000 sequences for training and 1,000 sequences for validation and testing.

#### USHCN

The USHCN dataset (Menne et al., )(^7^77[https://data.ess-dive.lbl.gov/view/doi:10.3334/CDIAC/CLI.NDP019](https://data.ess-dive.lbl.gov/view/doi:10.3334/CDIAC/CLI.NDP019) under CC BY 4.0) includes daily measurements from $1,218$ weather stations across the US, covering five variables: precipitation, snowfall, snow depth, and minimum and maximum temperature. We follow the pre-processing steps outlined in (De Brouwer et al., ), but select a subset of $1,168$ stations over a four-year period starting from $1990$, consistent with (Schirmer et al., ). Moreover, we make the time series irregular by subsampling $50\%$ of the time points and randomly removing $20\%$ of the measurements. We normalize the features to lie within the range $[0,1]$ and split into $60\%$ for training, $20\%$ for validation, and $20\%$ for testing.

#### Physionet

The Physionet dataset (Silva et al., )(^8^88[https://physionet.org/content/challenge-2012/1.0.0/](https://physionet.org/content/challenge-2012/1.0.0/) under ODC-BY 1.0) contains $8000$ multivariate clinical time-series obtained from the intensive care unit (ICU). Each time-series includes various clinical features recorded during the first $48$ hours after the patient’s admission to the ICU. We preprocess the data as in (Rubanova et al., ). Although the dataset contains a total of $41$ measurements, we eliminate $4$ static features, $\textit{i}.\textit{e}.,$ age, gender, height, and ICU-type, leaving $37$ time-varying features. We round the time-steps to $6$-minute intervals, following (Schirmer et al., ). We normalize the features to lie within the range $[0,1]$ and split into $60\%$ for training, $20\%$ for validation, and $20\%$ for testing.

**Table 4: Training Hyper-parameters**
| Dataset | Learning Rate | Train Epoch | Time Scale | Batch Size | $\mathbb{R}^{d}$ | $\#$ of base matrices (L) | $\#$ of parameters |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Human Activity | $1\times 10^{-3}$ | 400 | 1/221 | 256 | 288 | 256 | 1.65M |
| Pendulum | $1\times 10^{-3}$ | 500 | 0.1 | 50 | 20 | 15 | 19.6K |
| USHCN | $1\times 10^{-3}$ | 500 | 0.2 | 50 | 20 | 20 | 18.5K |
| Physionet | $1\times 10^{-3}$ | 500 | 0.3 | 100 | 24 | 20 | 28.5K |

### C.2 Masking Scheme

The masked attention mechanism for the assimilation schemes introduced in Sec 3.3, is illustrated in detail in Figure 4.

Figure: Figure 4: Illustraion of Masked Attention. Given the observed time stamps $\mathcal{T}=\{t_{i}\}_{i=1}^{3}$ ($\circ$) and the unseen time stamps ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathcal{T}_ {u}}$ ($\times$), attention scores are computed based on the latent observations $\{\mathbf{y}_{t}\}_{t\in\mathcal{T}}$. (Up) The mask of the History assimilation scheme consists of two components: one masks attention scores related to future time stamps, and the other masks those related to unseen time stamps. (Down) The mask of the Full assimilation scheme only blocks the attention scores corresponding to unseen time stamps. Using these masks, masked attention is calculated. For History assimilation scheme, the latent variables $\mathbf{z}_{t_{i}}$ include information up to time $t_{i}$, while for Full assimilation scheme, $\mathbf{z}_{t_{i}}$ incorporates all available information. Finally, Latent variables corresponding to unseen time stamps are filled with the nearest past latent variable value.
Refer to caption: x4.png

Figure: Algorithm 1 Training ACSSM.

Figure: Algorithm 2 Inference with ACSSM.

### C.3 Training details

#### Training

For all experiments, except for human activity classification, we followed the same experimental setup as CRU (Schirmer et al., )(^9^99[https://github.com/boschresearch/Continuous-Recurrent-Units](https://github.com/boschresearch/Continuous-Recurrent-Units) under AGPL-3.0 license). For the human activity classification task, we used the setup described in mTAND (Shukla & Marlin, )(^10^1010[https://github.com/reml-lab/mTAN](https://github.com/reml-lab/mTAN) under MIT license). For a fair comparison, we kept the number of model parameters similar to mTAND for the Human Activity dataset and CRU for the other datasets. The model was trained using the Adam optimizer (Diederik, ) in all experiments. For the per-point classification and regression tasks, we applied a weight decay of $1\times 10^{-2}$ and applied gradient clipping for classification task, while no weight decay was used for other tasks. Additionally, to prevent overfitting, we limited the training epochs to 200 for the Physionet extrapolation experiments. The remaining training hyper-parameters are detailed in Table 4.

To estimate the objective function (28), for the decoder $p_{\psi}(\cdot|\mathbf{y})=\mathcal{N}(\mathbf{p}_{\psi}(\mathbf{y}),\mathbf{
\Sigma}_{p})$, we used a Gaussian likelihood with a fixed variance of $\mathbf{\Sigma}_{p}=0.01\cdot\mathbf{I}$, following the approach outlined in (Rubanova et al., ) for the Pendulum, USHCN, and Physionet datasets. For the Human Activity dataset, we employed a categorical likelihood. Additionally, for $\mathcal{L}(\alpha)$, the initial conditions $(\mathbf{m}_{0},\mathbf{\Sigma}_{0})$ of the latent state were trainable parameters, initialized randomly. The covariance $\mathbf{\Sigma}_{0}$ was set using an exponential transformation, with a small constant ($\epsilon=10^{-6}$) added to ensure positivity. The latent state transition matrix $\mathbf{A}(\cdot)$ was initialized as $\mathbf{A}=\sum_{l=1}^{L}\mathbf{E}\mathbf{D}_{l}\mathbf{E}^{\top}$, where $\mathbf{E}$ was initialized orthonormally following (Lezcano Casado, ), and the diagonal matrices $\{\mathbf{D}_{l}\}_{l\in[1:L]}$ were initialized randomly and passed through a negative exponential to keep the values negative. For the potentials $\{g_{i}\}_{i\in[1:k]}$, we used a Gaussian likelihood $g_{t}(\cdot|\mathbf{x})=\mathcal{N}(\cdot|r_{t}(\mathbf{m}_{t},\mathcal{H}_{t}
),\mathbf{\Sigma}_{g})$, where $r_{t}(\mathbf{x})=\mathbf{M}_{\theta}\mathbf{m}_{t}+\mathbf{z}_{t}$ with a trainable matrix $\mathbf{M}_{\theta}$, transformer output $\mathbf{z}_{t}=\mathbf{T}_{\theta}(\mathcal{H}_{t})$, and a covariance $\mathbf{\Sigma}_{t}$.

#### Architecture

In all experiments except for the Pendulum dataset, the time series $\mathbf{o}_{0,T}$ was provided with the observation mask concatenated. We was used a dropout rate of 0.2 for a Human Activity, while no dropout rate was used for the other experiments. For our method, the networks used for each dataset are listed in below, where d is the dimension of latent space $\mathbb{R}^{d}$ as described in Table 4.

- •
Human Activity (Input size, I=12)
–
Encoder network $\mathbf{q}_{\phi}$: Input(2$\times$I) $\to$ Linear(d) $\to$ ReLU() $\to$ Linear(d)
–
Transformer $\mathbf{T}_{\theta}$:
Input(d) $\to$ Masked Attn(d) $\to$ FFN(d) $\to$ GRU(d)
–
Weight network $\mathbf{f}_{\theta}$:Input(d) $\to$ Linear(d) $\to$ Softmax()
–
Decoder network $\mathbf{p}_{\psi}$: Input(d) $\to$ Linear(d)
- •
Pendulum (Input size, I=1$\times$24$\times$24)
–
Encoder network $\mathbf{q}_{\phi}$: Input(I) $\to$ Conv2d(1, 12, kernel$_$size=5, stride=4, padding=2) $\to$ ReLU() $\to$ MaxPool2d(kernel$_$size=2, stride=2 $\to$ Conv2d(12,12, kernel$_$size=3, stride=2, padding=1) $\to$ ReLU() $\to$ MaxPool2d(kernel$_$size=2, stride=2) $\to$ Flatten(108) $\to$ Linear(d)
$\to$ ReLU()
–
Transformer $\mathbf{T}_{\theta}$:
$\texttt{Input(d)}\to\texttt{4}\times[\texttt{Masked Attn(d)}\to\texttt{FFN(d)}
]\to\texttt{GRU(d)}$
–
Weight network $\mathbf{f}_{\theta}$:Input(d) $\to$ Linear(d) $\to$ Softmax()
–
Decoder network $\mathbf{p}_{\psi}$: Input(d) $\to$ Linear(d) $\to$ ReLU() $\to$ Linear(d)
- •
USCHCN $\&$ Physionet (Input size, I=5 for USHCN and I=37 for Physionet)
–
Encoder network $\mathbf{q}_{\phi}$: Input(2$\times$I) $\to$ Linear(d) $\to$ ReLU() $\to$ LayerNorm(d) $\to$ Linear(d) $\to$ ReLU() $\to$ LayerNorm(d) $\to$ Linear(d) $\to$ ReLU() $\to$ LayerNorm(d) $\to$ Linear(d)
–
Transformer $\mathbf{T}_{\theta}$:
$\texttt{Input(d)}\to\texttt{4}\times[\texttt{Masked Attn(d)}\to\texttt{FFN(d)}
]\to\texttt{GRU(d)}$
–
Weight network $\mathbf{f}_{\theta}$:Input(d) $\to$ Linear(d) $\to$ Softmax()
–
Decoder network $\mathbf{p}_{\psi}$: Input(d) $\to$ Linear(d) $\to$ ReLU() $\to$ Linear(d) $\to$ ReLU() $\to$ Linear(d)$\to$ ReLU() $\to$ Linear(d)
- •
FFN
–
Input(d) $\to$ LayerNorm(d) $\to$ Linear(d) $\to$ GeLu() $\to$ Linear(d) $\to$ Residual(Input(d))
- •
Masked Attn
–
Input(Q, K, V, mask=M) $\to$ Normalize(Q) $\to$ Linear(Q) $\to$ Linear(K) $\to$ Linear(V) $\to$ Attention(Q, K) $\to$ Masking(M) $\to$ Softmax(d) $\to$ Dropout() $\to$ Matmul(V) $\to$ LayerNorm(d) $\to$ Linear(d) $\to$ Residual(Q)

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 8
presentation: 4
questions: - How many learnable parameters are use for each method? Was each method trained with a similar paramter budget?
- Given several approximations, how does the method perform on larger datasets ?
soundness: 4
strengths: - The paper makes significant theoretical contributions which should be relevant other researchers in this field. 
- The proposed approach has impressive empirical results showing better results with faster training while being theoretically grounded.
summary: The authors propose ACSSM approach for modeling irregular time series which uses continuous-discrete state space models (CD-SSMs).  The authors extend doob's-h transform to multi-marginal case and solve the problem using stochastic optimal control. The formulation leads to a  evidence lower bound (ELBO) and the authors propose a VI based loss function to model the irregular time series. 

Authors make linear approximation for the SDEs, which allows them to be able to perform simulation free estimation and exploit transformers for parallel computations. Overall, the proposed approach is faster and performs better than existing simulation based approaches, neural differential equations and other RNN based approaches tailored to handle irregular time series.
weaknesses: One of the weakness of this work that I see is that it makes several simplifying approximations to make the solution faster/tractable. Authors already acknowledge this but it would have been nice to understand if there are any practical trade-off due to these approximation.

---

confidence: 2
contribution: 4
initial_rating: 8
presentation: 2
questions: - The introduction lists the contributions, but it could benefit from more intuition to guide the reader through the chain of thought. For instance, there is no explanation of what a Feynman-Kac model is or how it facilitates sequential analysis.

- The multi-marginal Doob’s  h-transform is a central component of the approach, but its presentation lacks intuitive guidance, which is reflective of the overall style in the paper. Adding more accessible explanations would enhance understanding.

- Some sentences have syntax errors or missing words. I recommend proofreading the text to improve readability.
soundness: 4
strengths: The diffusion and stochastic optimal control (SOC) perspective for time-series modeling offers a compelling alternative to classical recurrent network methods, as it allows for continuous, flexible representations of latent dynamics that are suited to handle irregular timesteps.

The paper demonstrates superior performance on real-world datasets. Additionally, the efficiency gains in training time make it feasible for high-dimensional data, where classical SDE-based methods might be expensive.
summary: The paper introduces the Amortized Control of Continuous State Space Model (ACSSM) to handle irregularly sampled time-series. ACSSM aims to model the path measure of trajectories in a latent space, conditioned on observations in the data space. By using a latent space, the model captures a flexible and structured representation of the underlying dynamics that generate observed data, which is especially useful for irregularly sampled observations. To construct this conditional path measure, the authors introduce a novel multi-marginal Doob’s h-transform. This extension of the traditional Doob’s transform induces a class of stochastic differential equations (SDEs) that define the desired path measure in the latent space. However, simulating these SDEs directly is computationally infeasible due to the need for expensive normalization constants and conditional expectations. To overcome this challenge, the authors leverage stochastic optimal control to define a variational objective that approximates the optimal control needed to produce the conditioned dynamics. To further enhance computational efficiency, they propose working with affine linear SDEs with known Gaussian perturbation kernels, allowing simulation-free estimation of the latent trajectories and significantly speeding up the inference process.
weaknesses: The absence of a Related Work section limits the reader’s ability to understand how ACSSM compares to existing time-series modeling approaches, particularly those used in irregular time-series contexts, and those compared in the experiments section (e.g., recurrent networks, attention mechanisms, and previous SDE-based models).

The paper leans heavily on measure-theoretic concepts and complex SDE formulations, which could make it difficult for readers not specialized in these areas. More accessible explanations or visual intuitions could enhance understanding.

It remains to be seen how the model scales in practice with very high-dimensional or complex latent spaces, as affine SDE simplifications may reduce the expressiveness of the dynamics in these settings.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 2
questions: -Did you guys try using the affine linear drift function for the latent dynamics? How much better does the learned NN drift function do? 

-It is sometimes unclear where you are using the full and when you are using the history attention mechanism. It is claimed that the authors perform better than CRU because they aren’t just using past information, but also future information. However, this seems like an unfair comparison, am I misunderstanding something?

-How much does the "Parallel Computation" stuff increase the speed of inference? Is this a major contributor of why the method is so fast compared the other methods.
soundness: 3
strengths: -Novel method of solving Continuous State Space model which learns linear dynamics to accommodate irregular time series. 

-Theoretical analysis for the provided algorithm, including the derivation of an ELBO used to solve the dynamics.

-Demonstration of the algorithm on real time series, showing improved performance over other methods, in both Test MSE and also compute time (<5 secs).

-Addresses limitations, such as the errors accumulated from the linear approximation.
summary: The authors introduce a way of amortizing the controller of a state space model (ACSSM) to make it compatible with irregular time series. To do this, they generalize the single-marginal Doob’s $h$-transform to the multi-marginal case. To simulate the resulting continuous dynamics, they use VI to get an ELBO which is then optimized. To make this tractable, they assume that the latent dynamics are locally linear, and use a neural network to get expressive latent dynamics this way. The authors provide a theoretical analysis for this work, showing that the ELBO they obtain is tight, and then offer several experiments where this method shows improved performance over comparable baselines, for both per time classification/regression as well as sequence interpolation/extrapolation.
weaknesses: -Some of the background was a bit hard to follow. As someone who is relatively unfamiliar with the literature in this area, I found that some of the stuff in the methods section were not explained (it’s possible that it was common knowledge). I tried looking at the appendix for a more fleshed out explanation and i still can’t say I’m confident I understand everything going on.

-The paper argues that this method is scalable, and one thing I wish I could understand better is how this would work on environments with more complex dynamics, and possibly even partially observable environments. Did you guys try anything along this route?

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: 1. Can the authors elaborate on the choice of the cost function in Eq(7)? How should one interpret this cost function? Aside from the theoretical benefit, is there any intuition on the choice of this cost function? 
2. Can the author explain the usage of full attention in this scenario? If the full attention is applied to estimate the latent dynamic, does it mean the controls $\alpha_{t_i}$ for the controlled SDE is informed by the future y-observations $o_{t_j}, j > i$? Intuitively the masked attention makes sense to me, but I am unsure about the application of full attention. 
3. It seems like the main experiments are all obtained with the full attention scheme. Can the author provide the result of those experiments using masked attention as an ablation study? When is the masked attention used?
soundness: 3
strengths: 1. The paper is theoretically motivated with novel ideas and methods. 
2. The presentation of the paper is clear. 
3. The paper provides improvements on empirical results. 
4. The paper provides clear instructions on replicating the experiment.
summary: The author models irregular time series through a latent controlled SDE. The author first showed that the problem can be solved by a conditioned SDE, and then established the link between such conditioned SDE to a controlled SDE that optimizes a specific form of cost function. Finally, solving the HJB equation for the cost function can be obtained by simulating the controlled SDE through the Feynman-Kac theorem. The simulation was then simplified by analytically solving for the first and second moments at each timestep. Overall the paper is theoretically well-motivated with good empirical results.
weaknesses: 1. The motivation for specific implementation details is not clear

Meta-review: This paper proposes the Amortized Control of continuous State Space Model for continuous dynamical modeling of time series for irregular and discrete observations. It extends Doob's h-transform to the multi-marginal setting, and defines a variational inference algorithm with a tight ELBO. To speed up training and inference, ACSSM assumes locally linear latent dynamics and employs transformers for parallel computation.

All reviewers praise its novel and significant theoretical contribution to modeling continuous state space with irregular observations. The empirical studies on real-world datasets show superior performance in both computational efficiency and accuracy.

I recommend accepting this paper based on the consensus from reviewers.

Decision: Accept (Oral)