Learning Spatiotemporal Dynamical Systems from Point Process Observations

## Abstract

Abstract Spatiotemporal processes are a fundamental tool for modeling dynamics across various domains, from heat propagation in materials to oceanic and atmospheric flows. However, currently available neural network-based modeling approaches fall short when faced with data collected randomly over time and space, as is often the case with sensor networks in real-world applications like crowdsourced earthquake detection or pollution monitoring. In response, we developed a new spatiotemporal method that effectively handles such randomly sampled data. Our model integrates techniques from amortized variational inference, neural differential equations, neural point processes, and implicit neural representations to predict both the dynamics of the system and the probabilistic locations and timings of future observations. It outperforms existing methods on challenging spatiotemporal datasets by offering substantial improvements in predictive accuracy and computational efficiency, making it a useful tool for modeling and understanding complex dynamical systems observed under realistic, unconstrained conditions.

## 1 Introduction

Figure: Figure 1: Top: Previous models assume dense observations at every time point and fixed spatiotemporal grids. Bottom: Our model works with extremely sparse observations and predicts where the next observations will happen.
Refer to caption: extracted/5636476/figures/grid_many_modes_vs_ours_one_node.png

In this work, we are dealing with the modeling of spatiotemporal dynamical systems whose dynamics are driven by partial differential equations. Such systems are ubiquitous and range from heat propagation in microstructures to the dynamics of oceanic currents. We use the data-driven modeling approach, where we observe a system at various time points and spatial locations and use the collected data to learn a model. In practice, the data is often collected by sensor networks that make measurements at random time points and random spatial locations. Such a measurement approach is used, for example, in crowdsourced earthquake monitoring where smartphones are used as measurement devices (Minson et al., ; Kong et al., ), in oceanographic monitoring where measurements are made by floating buoys (Albaladejo et al., ; Xu et al., ; Marin-Perianu et al., ), and for air pollution monitoring with vehicle-mounted sensors (Ma et al., ; Ghanem et al., ). This approach offers several advantages as it needs no sensor synchronization and allows the sensors to move freely. However, its random nature makes modeling more challenging as it requires us to capture both the system dynamics and the random observation process.

To the best of our knowledge, existing neural network-based methods cannot model randomly observed spatiotemporal dynamical systems, mainly for two reasons. First, they do not model observation times and locations, hence they cannot predict where and when the next observation will be made. Second, they assume the sensors form a fixed spatial grid and record data simultaneously, which is not the case in our setup where the data could come from as little as a single sensor at each time point (Fig. 1). For example, some earlier methods assume the observations are made on a fixed and regular spatiotemporal grid (Long et al., ; Geneva and Zabaras, ), other methods can handle irregular, but still fixed, observation locations (Iakovlev et al., ; Lienen and Günnemann, ), and other works go further and allow the observation locations to change over time (Pfaff et al., ; Yin et al., ) but fix the observation times and assume dense observations. Whereas, another line of research has proposed methods that model only the spatiotemporal observation process (Chen et al., ; Zhu et al., ; Zhou et al., ; Zhou and Yu, ; Du et al., ).

We fill the existing gap and propose a model for randomly observed spatiotemporal dynamical systems. Our model uses ideas from amortized variational inference (Kingma and Welling, ), neural differential equations (Chen et al., ; Rackauckas et al., ), neural point processes (Mei and Eisner, ; Chen et al., ), and implicit neural representations (Chen et al., ; Yin et al., ) to efficiently learn the underlying system dynamics as well as to model the random observation process. Our model maps initial observations to the latent initial state via a Transformer encoder (Vaswani et al., ). Then, it computes the latent trajectory using neural ODEs (Chen et al., ). Finally, it uses implicit neural representations and points along the latent trajectory to parameterize the point process and observation distribution.
Our model shows strong empirical results outperforming other models from the literature on challenging spatiotemporal datasets.

## 2 Background

### 2.1 Spatiotemporal Point Processes

Spatiotemporal point processes (STPP) model sequences of events occurring in space and time. Each event has an associated event time $t_{i}\in\mathbb{R}_{\geq 0}$ and event location ${\bm{x}}_{i}\in\mathbb{R}^{d}$. Given an event history $\mathcal{H}_{t}\triangleq\{(t_{i},{\bm{x}}_{i})\ |\ t_{i}<t\}$ with all events up to time $t$, we can characterize an STPP by its conditional intensity function

$$ $\displaystyle\lambda^{*}(t,{\bm{x}})\triangleq\lim_{\delta t\downarrow 0,\ \delta r\downarrow 0}\frac{\mathbb{P}(t_{i}\in[t,t+\delta t],{\bm{x}}_{i}\in B _{\delta r}({\bm{x}})|\mathcal{H}_{t})}{\delta t|B_{\delta r}({\bm{x}})|},$ (1) $$

where $\delta t$ denotes an infinitesimal time interval, and $B_{\delta r}({\bm{x}})$ denotes a $\delta r$-ball centered at ${\bm{x}}$. Given a history $\mathcal{H}_{t}$ with $i-1$ events, $\lambda^{*}(t,{\bm{x}})$ describes the instantaneous probability of the next, $i$th, event occurring at time $t$ and location ${\bm{x}}$. Given a sequence of $N$ events $\{(t_{i},{\bm{x}}_{i})\}_{i=1}^{N}$ on a bounded domain $A\subset[0,T]\times\mathbb{R}^{d}$, the log-likelihood for the STPP is evaluated as (Daley et al., )

$$ $\displaystyle\log p(\{(t_{i},{\bm{x}}_{i})\}_{i=1}^{N})=\sum_{i=1}^{N}{\log \lambda^{*}(t_{i},{\bm{x}}_{i})}-\int_{A}{\lambda^{*}(t,{\bm{x}})d{\bm{x}}dt}.$ (2) $$

### 2.2 Ordinary and Partial Differential Equations

Given a deterministic continuous-time dynamic system with state ${\bm{z}}(t)$, we can describe the evolution of its state in terms of an ordinary differential equation (ODE)

$$ $\displaystyle\frac{d{\bm{z}}(t)}{dt}=f(t,{\bm{z}}(t)).$ (3) $$

For an initial state ${\bm{z}}_{1}$ at time $t_{1}$ we can solve the ODE to obtain the system state ${\bm{z}}(t)$ at later times $t>t_{1}$. The solution exists and is unique if $f$ is continuous in time and Lipschitz continuous in state (Coddington et al., ), and can be obtained either analytically or using numerical ODE solvers (Heirer et al., ). In this work we solve ODEs numerically using ODE solvers from torchdiffeq (Chen, ) package. Similarly, the dynamics of spatiotemporal systems with state ${\bm{z}}(t,{\bm{x}})$ defined over both space and time is described in terms of a partial differential equation (PDE)

$$ $\displaystyle\frac{\partial{\bm{z}}(t,{\bm{x}})}{\partial t}=F(t,{\bm{x}},{\bm {z}}(t,{\bm{x}}),\nabla{\bm{z}}(t,{\bm{x}}))$ (4) $$

which incorporates both temporal and spatial derivatives, indicated by $\frac{\partial}{\partial t}$ and $\nabla$, respectively.

## 3 Problem Setup

In this work, we model spatiotemporal systems from data. The data consists of multiple trajectories collected by observing the system over a period of time. A trajectory consists of $N$ triplets $\{(t_{i},{\bm{x}}_{i},{\bm{y}}_{i})\}_{i=1}^{N}$, where ${\bm{y}}_{i}\in\mathbb{R}^{d_{\bm{y}}}$ is a system observation, with $t_{i}\in\mathbb{R}_{\geq 0}$ and ${\bm{x}}_{i}\in\mathbb{R}^{d_{\bm{x}}}$ being the corresponding observation time and location. The number of observations $N$ can change across different trajectories in the dataset. Due to randomness of the observation process we assume the observation times and locations do not overlap (neither within a trajectory nor across trajectories), resulting in a single observation ${\bm{y}}_{i}$ per time point and location. For brevity, we describe our method for a single observed trajectory, but extension to multiple trajectories is straightforward.

We assume the data generating process (DGP) involves a latent state ${\bm{u}}(t,{\bm{x}})\in\mathbb{R}^{d_{{\bm{u}}}}$ that is defined for $t\in\mathbb{R}_{\geq 0}$ and ${\bm{x}}\in\mathbb{R}^{d_{\bm{x}}}$. We assume the following DGP with some initial time $t_{1}$:

$$ $\displaystyle{\bm{u}}(t_{1},\cdot)$ $\displaystyle\sim p({\bm{u}}),\quad\frac{\partial{\bm{u}}(t,{\bm{x}})}{ \partial t}=F({\bm{u}}(t,{\bm{x}}),\nabla{\bm{u}}(t,{\bm{x}}))),$ (5) $\displaystyle t_{i},{\bm{x}}_{i}$ $\displaystyle\sim\mathrm{nhpp}(\lambda({{\bm{u}}(t,{\bm{x}})})),\quad$ $\displaystyle i=1,\dots,N,$ (6) $\displaystyle{\bm{y}}_{i}$ $\displaystyle\sim p({\bm{y}}_{i}|g({\bm{u}}(t_{i},{\bm{x}}_{i}))),\quad$ $\displaystyle i=1,\dots,N.$ (7) $$

In that process, we first sample the latent field ${\bm{u}}(t_{1},\cdot)$ from $p({\bm{u}})$ and then specify its dynamics by a PDE (Eq. 5). Next, we sample the observation times $t_{i}$ and locations ${\bm{x}}_{i}$ from a non-homogeneous Poisson process (nhpp) with intensity $\lambda$ that is a function of the latent state (Eq. 6). Finally, we sample ${\bm{y}}_{i}$ from the observation distribution parameterized by a function $g$ (Eq. 7). All datasets in this work were generated according to this process (see Appendix A for details). We assume the data generating process is fully unknown, and our goal is to construct and learn its model from the data.

##### Limitations.

Using Poisson process for sampling event times and locations can be seen as a limitation in some applications because the assumption of non-overlapping time points may not hold in some real-world scenarios where effectively simultaneous observations might occur. Similarly, the lack of interaction between the event occurrence and system dynamics may not hold if observations affect the system, for example by triggering human actions or affecting sensor availability.

## 4 Methods

In this section we describe our proposed model, its components, and the parameter inference method.

### 4.1 Model

Our goal is to model data generated by the true space-time continuous DGP (Eqs. 5-7). Our generative model is defined as:

$$ $\displaystyle{\bm{z}}(t_{1})$ $\displaystyle\sim p({\bm{z}}(t_{1})),\quad\frac{d{\bm{z}}(\tau)}{dt}=f({\bm{z} }(\tau)),$ (8) $\displaystyle{\bm{u}}(t,x)$ $\displaystyle=\phi({\bm{z}}(t),{\bm{x}}),$ (9) $\displaystyle t_{i},{\bm{x}}_{i}$ $\displaystyle\sim\mathrm{nhpp}(\lambda({{\bm{u}}(t,{\bm{x}})})),\quad$ $\displaystyle i=1,\dots,N,$ (10) $\displaystyle{\bm{y}}_{i}$ $\displaystyle\sim p({\bm{y}}_{i}|g({\bm{u}}(t_{i},{\bm{x}}_{i}))),\quad$ $\displaystyle i=1,\dots,N.$ (11) $$

We describe each of the model components (Eqs. 8-11) in detail below. A diagram of our model is shown in Figure 2.

##### Latent dynamics (Eq. 8 ).

Our goal is to model the latent PDE dynamics of the DGP (Eq. 5). To that end, we introduce a low-dimensional temporal latent state ${\bm{z}}(\tau)\in\mathbb{R}^{d_{\bm{z}}}$ with $\tau\in\mathbb{R}_{\geq 0}$ that encodes the spatiotemporal field ${\bm{u}}(\tau,\cdot)$, and introduce a neural network parameterized ODE governing its dynamics (Eq. 8). We use a low-dimensional latent state because it allows to simulate the dynamics considerably quicker than a full-grid spatiotemporal discretization (Wu et al., ; Yin et al., ). We observed that the simulation times to numerically solve the ODE were long and increasing with the number of observations. We found that this happens because of a bottleneck in existing ODE solvers caused by extremely dense time grids (such as in Fig. 3). Such time grids force the ODE solver to choose small step sizes which prevents it from adaptively selecting the optimal step size thus reducing its efficiency. To improve the computational efficiency, we move away from the original dense time grid $t_{1},\dots,t_{N}$ to an auxiliary sparse grid $\tau_{1},\dots,\tau_{n}$, with $n\ll N$, where the first and last time points coincide with the original grid: $\tau_{1}=t_{1}$, $\tau_{n}=t_{N}$ (see Fig. 2). Then, we numerically solve for the latent state ${\bm{z}}(t)$ only at the sparse grid, which allows the ODE solver to choose the optimal step size and results in approximately an order of magnitude (up to 9x) faster simulations (see Sec. 5). Finally, we use the latent state evaluations ${\bm{z}}(\tau_{1}),\dots,{\bm{z}}(\tau_{n})$ to approximate the original latent state ${\bm{z}}(t)$ at any time point via interpolation as $\tilde{{\bm{z}}}(t)=\mathrm{interpolate}(t;{\bm{z}}(\tau_{1}),\dots,{\bm{z}}(
\tau_{n}))$.

##### Latent state decoding (Eq. 9 ).

Next, we need to recover the latent spatiotemporal state ${\bm{u}}(t,{\bm{x}})$ from the low-dimensional representation $\tilde{{\bm{z}}}(t)$. We define the latent spatiotemporal state ${\bm{u}}(t,{\bm{x}})$ as a function $\phi(\tilde{{\bm{z}}}(t),{\bm{x}})$ of the latent state $\tilde{{\bm{z}}}(t)$ and evaluation location ${\bm{x}}$ (Eq. 9). In particular, we parameterize $\phi(\tilde{{\bm{z}}}(t),{\bm{x}})$ by an MLP which takes $\tilde{{\bm{z}}}(t)+\mathrm{proj}({\bm{x}})$ as the input and maps it directly to ${\bm{u}}(t,{\bm{x}})$, where $\mathrm{proj}({\bm{x}})\in\mathbb{R}^{d_{\bm{z}}}$ is a trainable linear mapping of ${\bm{x}}$ to $\mathbb{R}^{d_{\bm{z}}}$. This formulation results in space-time continuous ${\bm{u}}(t,{\bm{x}})$ which can be evaluated at any time point and spatial location, which is required to define the intensity function $\lambda({\bm{u}}(t,{\bm{x}}))$ and observation function $g({\bm{u}}(t,{\bm{x}}))$.

##### Intensity function (Eq. 10 ).

We parameterize the intensity function $\lambda({\bm{u}}(t,{\bm{x}}))$ by an MLP which takes ${\bm{u}}(t,{\bm{x}})$ as the input and maps it to intensity of the Poisson process (Eq. 10). To ensure non-negative values of the intensity function and improve its numerical stability we further exponentiate the output of the MLP and add a small constant to it.

##### Observation function (Eq. 11 ).

The observation function $g({\bm{u}}(t,{\bm{x}}))$ maps the latent spatiotemporal state ${\bm{u}}(t,{\bm{x}})$ to parameters of the observation model (Eq. 11). Therefore, its exact structure depends on the observation model. In this work the observation model is a normal distribution, where the variance is fixed and the observation function $g({\bm{u}}(t,{\bm{x}}))$ returns the mean.

Details of the model definition and parameterization can be found in Appendix B and C.

Figure: Figure 2: Diagram of our model. Black crosses indicate observation times, while white crosses indicate observation locations. The initial observations (context) are mapped by the encoder to the distribution of the latent initial state $q({\bm{z}}_{1})$. The initial state ${\bm{z}}_{1}$ is then sampled from that distribution and evolved through time using neural ODEs. For computational efficiency, the latent trajectory is evaluated only at the sparse time grid $\tau_{1},\dots,\tau_{n}$ and the latent state at other time points is evaluated via interpolation $\tilde{{\bm{z}}}(t)$. The latent state $\tilde{{\bm{z}}}(t)$ is then mapped by $\phi$ to the spatiotemporal state ${\bm{u}}(t,{\bm{x}})$, which is used to parameterize the point process and observation distribution (via mappings $\lambda$ and $g$, respectively) for predicting subsequent observation times and locations.
Refer to caption: extracted/5636476/figures/main_figure.png

### 4.2 Parameter and Latent State Inference

We infer the model parameters and posterior of the latent initial state ${\bm{z}}_{1}\triangleq{\bm{z}}(t_{1})$ using amortized variational inference (Kingma and Welling, ). To approximate the posterior of ${\bm{z}}_{1}$ we define an approximate posterior distribution $q({\bm{z}}_{1};\psi)$ with local variational parameters $\psi$. Since there are multiple trajectories in the dataset and each trajectory has its own initial state, we need to define and optimize $\psi$ separately for each trajectory. To avoid optimizing $\psi$ for each trajectory we use amortization and define $\psi$ as the output of an encoder

$$ $\displaystyle\psi=\mathrm{Encoder}(\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{ N_{\mathrm{ctx}}}),$ (12) $$

which maps initial observations $\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}}$ (or the latest $N_{\mathrm{ctx}}$ observations, which we call the context) to the local variational parameters. The context size is defined by $N_{\mathrm{ctx}}$ and could be different for each trajectory (see Appendix A for details). We discuss the structure of the encoder in the next section.

In variational inference we minimize the Kullback-Leibler divergence between the approximate posterior $q({\bm{z}}_{1};\psi)$ and the true posterior $p({\bm{z}}_{1}|\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}})$ (Blei et al., ):

$$ $\displaystyle\mathrm{KL}\left[q({\bm{z}}_{1};\psi)\lVert p({\bm{z}}_{1}|\{t_{i },{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}})\right],$ (13) $$

which is equivalent to maximizing the evidence lower bound (ELBO), which for our model is

$$ $\displaystyle\mathcal{L}$ $\displaystyle=\underbrace{\sum_{i=1}^{N}{\mathbb{E}_{q({\bm{z}}_{1};\psi)} \left[\ln p({\bm{y}}_{i}|t_{i},{\bm{x}}_{i},{\bm{z}}_{1})\right]}}_{\textit{(i )}\text{ Expected observation log-lik.}}+\underbrace{\mathbb{E}_{q({\bm{z}}_{1 };\psi)}\left[\sum_{i=1}^{N}\ln\lambda({\bm{u}}(t_{i},{\bm{x}}_{i}))-\int{ \lambda({\bm{u}}(t,{\bm{x}}))d{\bm{x}}dt}\right]}_{\textit{(ii)}\text{ Expected STPP log-lik.}}$ (14) $\displaystyle\quad-\underbrace{\mathrm{KL}[q({\bm{z}}_{1};\psi)\lVert p({\bm{z }}_{1})].}_{\textit{(iii)}\text{ KL between prior and posterior}}$ (15) $$

The ELBO is maximized wrt. the model and encoder parameters. Appendix B contains detailed derivation of the ELBO, and fully specifies the model and the approximate posterior. While the term (iii) can be computed analytically, computation of terms (i) and (ii) involves approximations: Monte Carlo integration for the expectations and intensity integral, and numerical ODE solvers for the solution of the initial value problems. Appendix B details the computation of ELBO.

### 4.3 Encoder

Our encoder maps the context $\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}}$ to the local variational parameters $\psi$ of $q({\bm{z}}_{1};\psi)$ using Eq. 12. The encoder works by first mapping elements of the context into a high-dimensional feature space, then processing the resulting sequence using a stacks of transformer encoder layers (Vaswani et al., ), and finally mapping the transformer’s output to variational parameters $\psi$. Next, we discuss this process in more details.

First, we convert the context sequence $\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}}$ to a sequence of vectors $\{{\bm{b}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}}$, where

$$ $\displaystyle{\bm{b}}_{i}=\mathrm{proj}(t_{i})+\mathrm{proj}({\bm{x}}_{i})+ \mathrm{proj}({\bm{y}}_{i}),$ (16) $$

with $\mathrm{proj}$ being separate trainable linear projections. We do not add positional encodings since $t_{i}$ and ${\bm{x}}_{i}$ already provide all required positional information. We further add an aggregation token [AGG] with a learnable representation and obtain the following input sequence: $\{\{{\bm{b}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}},\text{[AGG]}\}$. Next, a stack of transformer encoder layers processes that input sequence into the output sequence $\{\{\overline{{\bm{b}}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}},\overline{\text{[AGG]}}\}$. Finally, we map the output aggregation token $\overline{\text{[AGG]}}$ to the variational parameters $\psi$. This mapping depends on the form of the approximate posterior (see Appendix B for detailed specification of the mapping).

Figure: Figure 3: Examples of trajectories from our datasets. White crosses indicate the observation locations ${\bm{x}}$, marks on the horizontal time grid indicate the corresponding observation times $t$, and color-coded fields denote the system state with observations ${\bm{y}}$ at white crosses $(t,{\bm{x}})$. Top: Burger’s (1D) dataset. The whole trajectory is plotted, with 1-D measurement location ${\bm{x}}$ along the vertical and time $t$ along the horizontal directions. Middle and Bottom: Shallow Water (2D) and Navier-Stokes (2D) datasets, respectively. Due to the 2D nature of the systems, only snapshots of the system state are plotted.
Refer to caption: extracted/5636476/figures/dataset_traj_vis.png

## 5 Experiments

In this section we demonstrate properties of our method and compare it against other methods from the literature. Our datasets are generated by three commonly-used PDE systems: Burgers’ (models nonlinear 1D wave propagation), Shallow Water (models 2D wave propagation under the gravity), and Navier-Stokes with transport (models the spread of a pollutant in a liquid over a 2D domain). The observation times and locations are sampled from a non-homogeneous Poisson process with state-dependent intensity function using the thinning algorithm (Lewis and Shedler, ). See Appendix A for details about the dataset generation. Figure 3 shows examples of trajectories from the datasets. For all datasets our model has at most 3 million parameters, and training takes at most 1.5 hours on a single GeForce RTX 3080 GPU. The training is done for 25k iterations with learning rate 3e-4 and batch size 32. We use the adaptive ODE solver (dopri5) from torchdiffeq package with relative and absolute tolerance set to 1e-5. As performance metrics we use mean absolute error (MAE) for the observations ${\bm{y}}_{i}$, and event-averaged log-likelihood of the point process for $t_{i}$ and ${\bm{x}}_{i}$, both evaluated on the test set. See Appendix C for our training setup and for details about the model architecture and hyperparameters.

Figure: Figure 4: Test MAE and log-likelihood vs context size.
Refer to caption: extracted/5636476/figures/exp_7_results.png

##### Context Size.

Our encoder maps the initial observations (context) $\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}}$ to parameters $\psi$ of the approximate posterior $q({\bm{z}}_{1};\psi)$. Here we look at how accuracy of the state predictions and process likelihood is affected by the context size. We train and test our model with different context sizes: full ($N_{\mathrm{ctx}}$ to 1), half of the context ($N_{\mathrm{ctx}}/2$ to 1), etc., and show the results in Figure 4. We see that both MAE and process likelihood quickly improve as we increase the context size, but the improvements saturate for larger contexts. This indicates that the encoder mostly uses observations that are close to the time point at which the latent state should be inferred, and does not utilize observations that are too far away from it.

Figure: Figure 5: Test MAE and process log-likelihood for different temporal grid resolutions ($n$) and interpolation methods (nn = nearest neighbor, lin = linear).
Refer to caption: extracted/5636476/figures/exp_4_results.png

##### Latent State Interpolation.

As discussed in Section 4.1, we do not evaluate the latent state at the full time grid $t_{1},\dots,t_{N}$ using the ODE solver (we call it the sequential method). Instead, we use the ODE solver to evaluate the latent state only at a sparser time grid $\tau_{1},\dots,\tau_{n}$, and then interpolate the evaluations to the full grid. Here we look at how this approach affects training times, and investigate the effects of the grid resolution $n$ and different interpolation methods. We vary the resolution from coarse $n=2$ to fine $n=50$, and test nearest neighbor and linear interpolation methods. Figure 6 shows that both MAE and process log-likelihood improve with the grid resolution, but only when we use nearest neighbor interpolation. With linear interpolation, the model achieves its best performance with just $n=2$, meaning that interpolation is done between only two points in the latent space, and further increasing the resolution does not improve the results. Both interpolation methods achieve the same optimal performance when $n\geq 20$. Furthermore, Table 6 shows that our method results in up to 9x faster latent state evaluation than the sequential method, which in our case translates in up to 4x faster training (computed with $n=20$).

##### Interaction Between Observation and Process Models.

**Table 1: Test MAE and log-likelihood of the full model, and after removing the point process (-STPP) or the observation model (-Obs.).**
| Dataset | MAE | MAE<br>(-STPP) | Log-lik. | Log-lik<br>(-Obs.) |
| --- | --- | --- | --- | --- |
| Burgers’ | 0.042 | 0.043 | 3.77 | 3.70 |
| Shallow Water | 0.012 | 0.012 | 3.67 | 3.56 |
| Navier-Stokes | 0.071 | 0.071 | 5.96 | 5.92 |

In our model the observation times and locations are modeled by a point process model (Eq. 10), while observations are modeled by an observation model (Eq. 11). Intuitively, one could expect some interaction between them and ask the question of how important is one model for the performance of the other. In this experiment we look at the effect of removing the point process model on the performance of the observation model and vice versa. We remove the point process model by removing Eq. 10 from our model, so that we do not model the observation times and locations. Similarly, we remove the observation model by removing Eq. 11 from our model and also removing observations ${\bm{y}}_{i}$ from our data, so that there is no information about the observation values. In Table 1 we show that removing the point process model seems to have practically no effect on MAE, while removing the observation model decreases the process likelihood. This shows that if we know the system states we can model the observation times and locations more accurately.

##### Comparison to Other Methods.

We compare our method against two groups of methods from the literature. The first group consists of neural spatiotemporal dynamic models that can model only the observations ${\bm{y}}_{i}$ but not the time points and spatial locations. This group includes FEN (Lienen and Günnemann, ) (graph neural network-based model, simulates the system dynamics at every point on the spatial grid), NSPDE (Salvi et al., ) (space- and time-continuous model simulating the dynamics in the spectral domain), and CNN-ODE (our simple baseline that uses a CNN encoder/decoder to map observations to/from the latent space, and uses neural ODEs (Chen et al., ) to model the latent space dynamics). Note that due to their assumption about multiple observations per time point, dynamic models require time binning (see Appendix C.4). The second group consists of neural spatiotemporal point processes that model only the observation times and locations. This group consists of DSTPP (Yuan et al., ), NSTPP (Chen et al., ), and AutoSTPP (Zhou and Yu, ). We also provide simple baselines for both groups: median predictor for dynamic models, and constant (learnable) intensity for the point process group. Hyperparameters of each method were tuned for the best performance. See Appendix C.4 for details about the models and hyperparameters used.

Table 2 shows the comparison results. We see that most methods from the first group perform rather poorly and fail to beat even the MAE of the median predictor, with the simple CNN-ODE model showing the best results. However, CNN-ODE still performs considerably worse than our model likely because of the loss of information caused by time binning. In the second group most methods showed poor results and failed to beat the simple constant intensity baseline, however AutoSTPP showed very strong performance, although still not outperforming our method.

These experiments show the challenging nature of modeling randomly observed dynamical systems, and demonstrate our method’s properties and strong performance.

**Table 2: Model comparisons. MAE and Log-lik (per event) on test data. The first group of methods (FEN, NSPDE, CNN-ODE) contains neural spatiotemporal dynamical models. The second group (DSTPP, NSTPP, AutoSTPP) contains neural spatiotemporal point process models. Error bars represent one standard error calculated over 5 realizations of a random seed controlling the model initialization.**
| Model | Burgers’ | Shallow Water | Navier-Stokes |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| MAE ($\downarrow$) | Log-lik. ($\uparrow$) | MAE ($\downarrow$) | Log-lik. ($\uparrow$) | MAE ($\downarrow$) | Log-lik. ($\uparrow$) |  |
| Median Pred. | 0.161 | - | 0.148 | - | 0.155 | - |
| FEN | N/A | - | $0.336\pm 0.008$ | - | $0.158\pm 0.006$ | - |
| NSPDE | $0.166\pm 0.002$ | - | $0.339\pm 0.002$ | - | $0.098\pm 0.001$ | - |
| CNN-ODE | $0.076\pm 0.000$ | - | $0.032\pm 0.001$ | - | $0.077\pm 0.002$ | - |
| Const. Intensity | - | 3.59 | - | 2.02 | - | 5.86 |
| DSTPP | - | $0.67\pm 0.13$ | - | $0.08\pm 0.04$ | - | $0.41\pm 0.04$ |
| NSTPP | - | $2.35\pm 0.06$ | - | $1.53\pm 0.13$ | - | $3.83\pm 0.14$ |
| AutoSTPP | - | N/A | - | $3.37\pm 0.01$ | - | $5.91\pm 0.01$ |
| Ours | $\bm{0.043}\pm\bm{0.002}$ | $\bm{3.77}\pm\bm{0.02}$ | $\bm{0.012}\pm\bm{0.001}$ | $\bm{3.62}\pm\bm{0.05}$ | $\bm{0.071}\pm\bm{0.001}$ | $\bm{5.96}\pm\bm{0.01}$ |

## 6 Related Work

##### Neural Point Processes.

Traditionally, intensity functions are simple parametric functions incorporating domain knowledge about the system being modeled. While simple and interpretable, this approach might require strong domain knowledge and might have limited expressivity due to overly simplistic form of the intensity function. These limitations lead to the development of neural point processes which parameterize the intensity function by a neural network, leading to improved flexibility and expressivity. For example Mei and Eisner (); Omi et al. (); Jia and Benson (); Zuo et al. () use neural networks to model time-dependent intensity functions in a flexible and efficient manner, with neural architectures ranging from recurrent neural networks to transformers and neural ODEs. These techniques consistently outperform classical intensity parameterizations and do not require strong domain expertise to choose the right form of the intensity function. Other works extend this idea to marked point processes where intensity is a function of time and also of a discrete or continuous mark. For example, Chen et al. (); Zhu et al. (); Zhou et al. (); Zhou and Yu (); Yuan et al. () use spatial coordinates as a mark and model the intensity using a wide range of methods aimed at improving flexibility and efficiency. Works such as Du et al. (); Xiao et al. (); Boyd et al. () use discrete observations as marks, and Sharma et al. () use continuous marks with underlying dynamic model, which makes their work most similar to ours, with major differences related to the model, training process, and them working only with temporal processes, whereas we are dealing with more general spatiotemporal processes.

##### Neural Spatiotemporal Dynamic Models.

To model the spatiotemporal dynamics our method uses the “encode-process-decode” approach employed in many other works. The main idea of the approach consists of taking a set of initial observations and using it to estimate the latent initial state. A dynamics function is then used to map the latent state forward in time to other time points either discretely (Long et al., ; Han et al., ; Wu et al., ) or continuously (Yildiz et al., ; Salvi et al., ). Finally, the latent state is mapped to the observations using discrete (Yildiz et al., ; Wu et al., ) or continuous (Yin et al., ; Chen et al., ) decoder. Such methods assume multiple observations per time point and do not model the observation times and locations.

##### Implicit Neural Representations.

Our model represents the continuous latent spatiotemporal state ${\bm{u}}(t,{\bm{x}})$ in terms of a function $\phi(\tilde{{\bm{z}}}(t),{\bm{x}})$. This approach to representing continuous field is called implicit neural representations (Park et al., ; Chen and Zhang, ; Mescheder et al., ; Chibane et al., ), and it has found applications in many other works related to modeling of spatiotemporal systems (Chen et al., ; Yin et al., ).

## 7 Conclusion

In this work, we developed a method for modeling of randomly observed spatiotemporal dynamical systems. Along with our method, we proposed an interpolation-based technique to greatly speed up its training times. In the experiments we showed that our method effectively utilizes the context size of various length to improve accuracy of the initial state inference, and that it uses the system state observations to improve the modeling of observation times and locations. We further demonstrated that our method achieves strong performance on challenging spatiotemporal datasets and outperforms other methods from the literature.

## Appendix A Data

As described in Section 3, we assume the following data generating process

$$ $\displaystyle{\bm{u}}(t_{1},\cdot)\sim p({\bm{u}}),\quad\frac{\partial{\bm{u}} (t,{\bm{x}})}{\partial t}=F({\bm{u}}(t,{\bm{x}}),\nabla{\bm{u}}(t,{\bm{x}}))),$ (17) $\displaystyle t_{i},{\bm{x}}_{i}\sim\mathrm{nhpp}(\lambda({{\bm{u}}(t,{\bm{x}} )})),\quad$ $\displaystyle i=1,\dots,N,$ (18) $\displaystyle{\bm{y}}_{i}\sim p({\bm{y}}_{i}|g({\bm{u}}(t_{i},{\bm{x}}_{i}))),\quad$ $\displaystyle i=1,\dots,N,$ (19) $$

We selected three commonly used PDE systems: Burgers’, Shallow Water, and Navier-Stokes with transport. Next we discuss the data generating process for each dataset.

### A.1 Burgers’

Figure: Figure 7: Trajectory lengths for Burgers’ dataset.
Refer to caption: extracted/5636476/figures/traj_len_burg.png

The Burgers’ system in 1D is characterized by a 1-dimensional state $u(t,x)$ and the following PDE dynamics:

$$ $\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu \frac{\partial^{2}u}{\partial t^{2}},$ (20) $$

where $\nu$ is the diffusion coefficient. We obtain data for this system from the PDEBench dataset (Takamoto et al., ). The system is simulated on a time interval $[0,2]$ and on a spatial domain $[0,1]$ with periodic boundary conditions. The dataset contains 10000 trajectories, each trajectory is evaluated on a uniform temporal grid with 101 time points, and uniform spatial grid with 256 nodes.

Since the data is discrete, we define continuous field $u(t,x)$ by linearly interpolating the data across space and time. Next, we use the interpolant to define the intensity function

$$ $\displaystyle\lambda(u(t,x))=1600|u(t,x)|,$ (21) $$

and use it to simulate the non-homogeneous Poisson process using the thinning algorithm (Lewis and Shedler, ). As the result, for each the 10000 trajectories, we obtain a sequence of time points and spatial locations $\{t_{i},x_{i}\}_{i=1}^{N}$, where $N$ can be different for each trajectory. Finally, we compute the observations as

$$ $\displaystyle{\bm{y}}_{i}=u(t_{i},x_{i}).$ (22) $$

We use the first 0.5 seconds as the context for the initial state inference, and further filter the dataset to ensure that each trajectory has at least 10 points in the context, resulting in approximately 7000 trajectories. We use 80%/10%/10% split for training/validation/testing.

### A.2 Shallow Water

Figure: Figure 8: Trajectory lengths for Shallow Water dataset.
Refer to caption: extracted/5636476/figures/traj_len_sw.png

The Shallow Water system in 2D is characterized by a 3-dimensional state

$$ $\displaystyle{\bm{u}}(t,x)=[h(t,{\bm{x}}),u(t,{\bm{x}}),v(t,{\bm{x}})]$ (23) $$

and the following PDE dynamics:

$$ $\displaystyle\frac{\partial h}{\partial t}+\frac{\partial(hu)}{\partial x}+ \frac{\partial(hv)}{\partial y}=0,$ (24) $\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+g\frac{\partial h}{\partial x}=0,$ (25) $\displaystyle\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+g\frac{\partial h}{\partial y}=0,$ (26) $$

where $g$ is the gravitational acceleration, $h(t,{\bm{x}})$ is wave height, and $u(t,{\bm{x}})$ and $v(t,{\bm{x}})$ are horizontal and vertical velocities, respectively. We obtain data for this system from the PDEBench dataset (Takamoto et al., ). The system is simulated on a time interval $[0,1]$ and on a spatial domain $[-2.5,2.5]^{2}$. The dataset contains 1000 trajectories, each trajectory is evaluated on a uniform temporal grid with 101 time points, and uniform 128-by-128 spatial grid.

Since the data is discrete, we define continuous field ${\bm{u}}(t,x)$ by linearly interpolating the data across space and time. Next, we use the interpolant to define the intensity function

$$ $\displaystyle\lambda({\bm{u}}(t,{\bm{x}}))=500|h(t,{\bm{x}})-1|,$ (27) $$

indicating that measurements are made only when the wave height is below or above the baseline of 1. We use the intensity function to simulate the non-homogeneous Poisson process using the thinning algorithm (Lewis and Shedler, ). As the result, for each of the 1000 trajectories, we obtain a sequence of time points and spatial locations $\{t_{i},x_{i}\}_{i=1}^{N}$, where $N$ can be different for each trajectory. Finally, we compute the observations as

$$ $\displaystyle{\bm{y}}_{i}=h(t_{i},{\bm{x}}_{i}),$ (28) $$

observing only the wave height.

We use the first 0.1 seconds as the context for the initial state inference, and further filter the dataset to ensure that each trajectory has at least 10 points in the context, resulting in approximately 800 trajectories. We use 80%/10%/10% split for training/validation/testing.

### A.3 Navier-Stokes

Figure: Figure 9: Trajectory lengths for Navier-Stokes dataset.
Refer to caption: extracted/5636476/figures/traj_len_ns.png

The Navier-Stokes system with a transport equation is characterized by a 4-dimensional state

$$ $\displaystyle{\bm{u}}(t,{\bm{x}})=[c(t,{\bm{x}}),u(t,{\bm{x}}),v(t,{\bm{x}}),p (t,{\bm{x}})]$ (29) $$

and the following PDE dynamics:

$$ $\displaystyle\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$ (30) $\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}=-\frac{\partial p}{\partial x}+\nu\left(\frac{ \partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right),$ (31) $\displaystyle\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}=-\frac{\partial p}{\partial y}+\nu\left(\frac{ \partial^{2}v}{\partial x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}\right)+c,$ (32) $\displaystyle\frac{\partial c}{\partial t}=-u\frac{\partial c}{\partial x}-v \frac{\partial c}{\partial y}+\nu\left(\frac{\partial^{2}c}{\partial x^{2}}+ \frac{\partial^{2}c}{\partial y^{2}}\right),$ (33) $$

where the diffusion constant $\nu=0.002$, and $c(t,{\bm{x}})$ is concentration of the transported species, $p(t,{\bm{x}})$ is pressure, and $u(t,{\bm{x}})$ and $v(t,{\bm{x}})$ are horizontal and vertical velocities, respectively. For each trajectory, we start with zero initial velocities and pressure, and the initial scalar field $c_{0}(x,y)$ is generated as:

$$ $\displaystyle\tilde{c}_{0}(x,y)=\sum_{k,l=-N}^{N}{\lambda_{kl}\cos(2\pi(kx+ly) )+\gamma_{kl}\sin(2\pi(kx+ly))},$ (34) $\displaystyle c_{0}(x,y)=\frac{\tilde{c}_{0}(x,y)-\mathrm{min}(\tilde{c}_{0})} {\mathrm{max}(\tilde{c}_{0})-\mathrm{min}(\tilde{c}_{0})},$ (35) $$

where $N=2$ and $\lambda_{kl},\gamma_{kl}\sim\mathcal{N}(0,1)$.

We use PhiFlow (Holl et al., ) to solve the PDEs. The system is simulated on a time interval $[0,1]$ and on a spatial domain $[0,1]^{2}$ with periodic boundary conditions. The solution is evaluated at randomly selected spatial locations and time points. We use $1089$ spatial locations and $25$ time points. The spatial and temporal grids are the same for all trajectories. The dataset contains 1000 trajectories.

Since the data is discrete, we define continuous field ${\bm{u}}(t,x)$ by linearly interpolating the data across space and time. Next, we use the interpolant to define the intensity function

$$ $\displaystyle\lambda({\bm{u}}(t,{\bm{x}}))=2000|c(t,{\bm{x}})|,$ (36) $$

indicating that measurement intensity is proportional to the species concentration. We use the intensity function to simulate the non-homogeneous Poisson process using the thinning algorithm (Lewis and Shedler, ). As the result, for each of the trajectories, we obtain a sequence of time points and spatial locations $\{t_{i},x_{i}\}_{i=1}^{N}$, where $N$ can be different for each trajectory. Finally, we compute the observations as

$$ $\displaystyle{\bm{y}}_{i}=c(t_{i},{\bm{x}}_{i}),$ (37) $$

observing only the species concentrations.

We use the first 0.5 seconds as the context for the initial state inference, and further filter the dataset to ensure that each trajectory has at least 10 points in the context, resulting in all 1000 trajectories satisfying this condition. We use 80%/10%/10% split for training/validation/testing.

## Appendix B Model, Posterior, and ELBO

### B.1 Model

As described in Section 4.1 our model is defined as

$$ $\displaystyle{\bm{z}}_{1}\sim p({\bm{z}}_{1}),\quad\frac{d{\bm{z}}(\tau)}{dt}= f({\bm{z}}(\tau)),$ (38) $\displaystyle\tilde{{\bm{z}}}(t)=\mathrm{interpolate}(t;{\bm{z}}(\tau_{1}), \dots,{\bm{z}}(\tau_{n})),$ (39) $\displaystyle{\bm{u}}(t,x)=\phi(\tilde{{\bm{z}}}(t),{\bm{x}}),$ (40) $\displaystyle t_{i},{\bm{x}}_{i}\sim\mathrm{nhpp}(\lambda({{\bm{u}}(t,{\bm{x}} )})),\quad$ $\displaystyle i=1,\dots,N,$ (41) $\displaystyle{\bm{y}}_{i}\sim p({\bm{y}}_{i}|g({\bm{u}}(t_{i},{\bm{x}}_{i}))),\quad$ $\displaystyle i=1,\dots,N..$ (42) $$

The corresponding joint distribution is

$$ $\displaystyle p({\bm{z}}_{1},\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{N})=p( {\bm{z}}_{1})p(\{t_{i},{\bm{x}}_{i}\}_{i=1}^{N}|{\bm{z}}_{1})\prod_{i=1}^{N}{p ({\bm{y}}_{i}|{\bm{z}}_{1},t_{i},{\bm{x}}_{i})},$ (43) $$

where

$$ $\displaystyle p({\bm{z}}_{1})=\mathcal{N}(\bm{0},I),$ (44) $\displaystyle p(\{t_{i},{\bm{x}}_{i}\}_{i=1}^{N}|{\bm{z}}_{1})=\prod_{i=1}^{N} {\lambda({\bm{u}}(t_{i},{\bm{x}}_{i}))}\exp{\left(-\int{\lambda({\bm{u}}(t,{ \bm{x}}))d{\bm{x}}dt}\right)},$ (45) $\displaystyle p({\bm{y}}_{i}|{\bm{z}}_{1},t_{i},{\bm{x}}_{i})=\mathcal{N}(g({ \bm{u}}(t_{i},{\bm{x}}_{i})),\sigma_{{\bm{y}}}^{2}I),$ (46) $$

where $\mathcal{N}$ is the normal distribution, $\bm{0}$ is a zero vector, and $I$ is the identity matrix. We set $\sigma_{{\bm{y}}}$ to $10^{-3}$.

### B.2 Posterior

We define the approximate posterior with $\psi=[\psi_{\mu},\psi_{\sigma^{2}}]$ as

$$ $\displaystyle q({\bm{z}}_{1};\psi)=\mathcal{N}({\bm{z}}_{1}|\psi_{\mu},\mathrm {diag}(\psi_{\sigma^{2}})),$ (47) $$

where $\mathcal{N}$ is the normal distribution, and $\mathrm{diag}(\psi_{\sigma^{2}})$ is a matrix with vector $\psi_{\sigma^{2}}$ on the diagonal.

As discussed in Section 4.3, the encoder maps the context $\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}}$ to the local variational parameters $\psi$ (which we break up into $[\psi_{\mu},\psi_{\sigma^{2}}]$). The encoder first converts the context sequence $\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}}$ to a sequence of vectors $\{{\bm{b}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}}$ and then adds an aggregation token, which gives the following input sequence: $\{\{{\bm{b}}_{i}\}_{i=1}^{N_{\mathrm{ctx}}},\text{[AGG]}\}$. We pass it though the transformer encoder and read the value of the aggregation token at the last layer, which we denote by $\overline{\text{[AGG]}}$. Then, we map $\overline{\text{[AGG]}}$ to $[\psi_{\mu},\psi_{\sigma^{2}}]$ as follows:

$$ $\displaystyle\psi_{\mu}=\mathrm{Linear}(\overline{\text{[AGG]}}),$ (48) $\displaystyle\psi_{\sigma^{2}}=\exp\left({\mathrm{Linear}(\overline{\text{[AGG ]}})}\right),$ (49) $$

where $\mathrm{Linear}$ are separate linear mappings.

### B.3 ELBO

Given definitions in the previous sections, we can write the ELBO as

$$ $\displaystyle\mathcal{L}$ $\displaystyle=\int q({\bm{z}}_{1})\ln{\frac{p({\bm{z}}_{1},\{t_{i},{\bm{x}}_{i },{\bm{y}}_{i}\}_{i=1}^{N})}{q({\bm{z}}_{1})}}d{\bm{z}}_{1}$ (50) $\displaystyle=\int q({\bm{z}}_{1})\ln{\frac{\prod_{i=1}^{N}{p({\bm{y}}_{i}|{ \bm{z}}_{1},t_{i},{\bm{x}}_{i})}p(\{t_{i},{\bm{x}}_{i}\}_{i=1}^{N}|{\bm{z}}_{1 })p({\bm{z}}_{1})}{q({\bm{z}}_{1})}}d{\bm{z}}_{1}$ (51) $\displaystyle=\int q({\bm{z}}_{1})\ln{\prod_{i=1}^{N}{p({\bm{y}}_{i}|{\bm{z}}_ {1},t_{i},{\bm{x}}_{i})}}d{\bm{z}}_{1}$ (52) $\displaystyle\quad+\int q({\bm{z}}_{1})\ln{p(\{t_{i},{\bm{x}}_{i}\}_{i=1}^{N}| {\bm{z}}_{1})}d{\bm{z}}_{1}$ (53) $\displaystyle\quad+\int q({\bm{z}}_{1})\ln{\frac{p({\bm{z}}_{1})}{q({\bm{z}}_{ 1})}}d{\bm{z}}_{1}$ (54) $\displaystyle=\underbrace{\sum_{i=1}^{N}{\mathbb{E}_{q({\bm{z}}_{1})}\left[\ln p ({\bm{y}}_{i}|{\bm{z}}_{1},t_{i},{\bm{x}}_{i})\right]}}_{\mathcal{L}_{1}}$ (55) $\displaystyle\quad+\underbrace{\mathbb{E}_{q({\bm{z}}_{1})}\left[\sum_{i=1}^{N }\ln\lambda({\bm{u}}(t_{i},{\bm{x}}_{i}))-\int{\lambda({\bm{u}}(t,{\bm{x}}))d{ \bm{x}}dt}\right]}_{\mathcal{L}_{2}}$ (56) $\displaystyle\quad-\underbrace{\mathrm{KL}[q({\bm{z}}_{1})\lVert p({\bm{z}}_{1 })]}_{\mathcal{L}_{3}}.$ (57) $\displaystyle=\mathcal{L}_{1}+\mathcal{L}_{2}-\mathcal{L}_{3}.$ (58) $$

##### Computing ELBO.

We compute the ELBO as follows:

- 1.
Compute local variational parameter from the context as
$\psi_{\mu},\psi_{\sigma^{2}}=\mathrm{Encoder}(\{t_{i},{\bm{x}}_{i},{\bm{y}}_{i
}\}_{i=1}^{N_{\mathrm{ctx}}})$
- 2.
Sample the latent initial state as ${\bm{z}}_{1}\sim\mathcal{N}({\bm{z}}_{1}|\psi_{\mu},\mathrm{diag}(\psi_{\sigma
^{2}}))$ using reparameterization
- 3.
Evaluate ${\bm{z}}(\tau_{1}),\dots,{\bm{z}}(\tau_{n})$ by solving $\frac{d{\bm{z}}(\tau)}{dt}=f({\bm{z}}(\tau))$ with ${\bm{z}}_{1}$ as the initial condition. We solve the ODE using torchdiffeq package (Chen et al., ).
- 4.
Compute the latent states $\tilde{{\bm{z}}}(t_{1}),\dots,\tilde{{\bm{z}}}(t_{N})$ via interpolation as
$\tilde{{\bm{z}}}(t)=\mathrm{interpolate}(t;{\bm{z}}(\tau_{1}),\dots,{\bm{z}}(
\tau_{n}))$
- 5.
Compute ${\bm{u}}(t_{1},{\bm{x}}_{1}),\dots,{\bm{u}}(t_{N},{\bm{x}}_{N})$ as ${\bm{u}}(t_{i},{\bm{x}}_{i})=\phi(\tilde{{\bm{z}}}(t_{i}),{\bm{x}}_{i})$
- 6.
Compute $\mathcal{L}_{1}$ using Monte Carlo integration
- 7.
Compute $\mathcal{L}_{2}$ using Monte Carlo integration, with the intensity integral computed using torchquad (Gómez et al., ) package.
- 8.
Compute $\mathcal{L}_{3}$ analytically

Sampling is done using reparametrization (Kingma and Welling, ) to allow for exact gradient evaluation via backpropagation. Monte Carlo integration is done using sample size of one.

## Appendix C Experiments Setup

### C.1 Datasets

See Appendix A for all details about dataset generation. For all datasets we use 80%/10%10% train/validation/test splits. The context size $N_{\mathrm{ctx}}$ is determined for each trajectory separately based on the time interval that we consider to be the context, for Burgers’ and Navier-Stokes it is the first 0.5 seconds (out of 2 seconds), while for Shallow Water it is 0.1 seconds (out of 1 second). Since the number of events occurring within this fixed time interval is different for each trajectory, the context size consequently differs across trajectories. See encoder architecture description below for details on how different context lengths are processed.

### C.2 Training, Validation, and Testing

We use AdamW (Loshchilov and Hutter, ) optimizer with constant learning rate 3e-4 (we use linear warmup for first 250 iterations). We train for 25000 iterations, sampling a random minibatch from the training dataset at every iteration. We do not use any kind of scaling for the ELBO terms.

We compute validation error (MAE on observations) on a single random minibatch from the validation set every 125 iterations. We save the current model weights if the validation error averaged over last 10 validation runs is the smallest so far.

To simulate the model’s dynamics we use differentiable ODE solvers from torchdiffeq package (Chen et al., ). In particular, we use the dopri5 solver with $\text{rtol}=\text{atol}=10^{-5}$ without the adjoint method. To evaluate the intensity integral $\int{\lambda({\bm{u}}(t,{\bm{x}})d{\bm{x}}dt}$ we use Monte Carlo integration from the torchquad (Gómez et al., ) package with 32 randomly sampled points for training, and 256 randomly sampled points for testing, we found these sample sizes to be sufficient for producing small variance of the integral estimation in our case.

### C.3 Encoder and Model Architectures

##### Encoder.

We use padding to efficiently handle input contexts of varying length. The time points, coordinates, and observations are linearly projected to the embedding space (128-dim for Burgers’ and Shallow Water, and 192-dim for Navier-Stokes). Then, the projections are summed and passed through a stack of transformer encoder layers (4 layers for Burgers’ and Shallow Water, and 5 layers for Navier-Stokes, 4 attention heads were used in all cases). We introduce a learnable aggregation token which we append at the end of each input sequence. The final layer’s output is read from the last token of the output sequence (corresponds to the aggregation token) and is mapped to the local variational parameters as discussed in Appendix B.

##### Model.

The dynamics function $f({\bm{z}}(t))$ is an MLP with 3 layers 368 neurons each, and GeLU (Hendrycks and Gimpel, ) nonlinearities. We set resolution of the uniform temporal grid $\tau_{1},\dots,\tau_{n}$ to $n=50$. The dynamics function maps the current latent state ${\bm{z}}(t)$ to its time derivative $\frac{d{\bm{z}}(t)}{dt}$, both are $d_{\bm{z}}$-dimensional vectors.

The function $\phi(\tilde{{\bm{z}}}(t),{\bm{x}})$ is an MLP with 3 hidden layers with width 368 for Burgers’, 256 for Shallow Water, and 512 for Navier-Stokes. We use GeLU (Hendrycks and Gimpel, ) nonlinearities. As the input the MLP takes the sum of $\tilde{{\bm{z}}}(t)\in\mathbb{R}^{d_{{\bm{z}}}}$ with $\mathrm{proj}({\bm{x}})\in\mathbb{R}^{d_{{\bm{z}}}}$, where $\mathrm{proj}$ is a linear projection, and maps the sum to the output ${\bm{u}}(t,{\bm{x}})\in\mathbb{R}^{d_{{\bm{u}}}}$.

The intensity function $\lambda({\bm{u}}(t,{\bm{x}}))$ is an MLP with 3 layers, 256 neurons each, and GeLU (Hendrycks and Gimpel, ) nonlinearities. We further exponentiate the output of the MLP and add a small constant (1e-4) to it to ensure the intensity is positive and to avoid numerical instabilities.

We set the observation function $g({\bm{u}}(t,{\bm{x}}))$ to the identity function, thus $g({\bm{u}}(t,{\bm{x}}))={\bm{u}}(t,{\bm{x}})$.

### C.4 Model Comparison

#### C.4.1 Dynamic Spatiotemporal Models

In this section we discuss neural spatiotemporal models we used for the comparisons. But first, we discuss data preprocessing that we do in order to make these models applicable.

##### Data Preprocessing.

All models in this section assume that dense observations are available at every time point, and some further assume the observations are located on a uniform spatial grid. To make our data satisfy these requirements we use time binning and spatial interpolation. First, we use time binning and divide the original time grid into 10 bins and group all time points into 10 groups based on which bin they belong to. We used 10 time bins as it ensured that there was a sufficient number of time points in each bin while also being sufficient to capture changes of the system state. We denote by $t_{i}$ the time point that represents the bin $i$. Next, we use all observations that belong to bin $i$ and interpolate them to a uniform 32x32 spatial grid. In many cases some points of the spatial grid were outside of the convex hull of the interpolation points, so we set values of such points to -1 (values of the state range from 0 to 1). We used linear interpolation. We denote the resulting interpolated values at time $t_{i}$ as ${\bm{u}}_{i}$ which is the interpolated system state at the 32x32 spatial grid. When we compute the loss between ${\bm{u}}_{i}$ and the corresponding model prediction we use only those spatial locations that had proper interpolation values (i.e. were inside the convex hull of the interpolation points).

##### Median predictor.

This is a simple baseline where we compute the median of the observations ${\bm{y}}_{i}$ on the training data and use it as the prediction on the test data.

##### CNN-ODE.

This is another relatively simple baseline based on latent neural ODEs (Chen et al., ). It operates in some sense similarly to our model as it takes a sequence of initial observations (context) ${\bm{y}}_{1},\dots,{\bm{y}}_{N_{\text{ctx}}}$ and uses a CNN encoder to map the context to the latent initial state ${\bm{z}}_{1}$. Then it maps ${\bm{z}}_{1}$ to ${\bm{z}}_{2},\dots,{\bm{z}}_{N}$ via a latent ODE, and finally decodes ${\bm{z}}_{i}$ to obtain ${\bm{u}}_{i}$ via a CNN decoder. In our case we used the first two observations ${\bm{u}}_{1}$ and ${\bm{u}}_{2}$ as the context.

##### FEN.

This model discretizes the spatial domain into a mesh and uses a graph neural network to model the dynamics of the state at each node of the mesh. The model starts at ${\bm{u}}_{1}$ and maps it directly to ${\bm{u}}_{2},\dots,{\bm{u}}_{N}$ via an ODE solver without using any context. We simulated the dynamics as stationary and autonomous, with the free-form term active. The dynamics function is a 2-layer MLP with the width of 512 neurons and ReLU nonlinearities. All other hyper-parameters were left as defaults as we found that changing them did not improve the results. We used the official code from Lienen and Günnemann ().

##### NSPDE.

This model applies a point-wise transformation to map the system state at every spatial node to a high-dimensional latent representation, then uses the fixed-point method to simulate the dynamics at every node, and finally maps the high-dimensional latent states back to the observation space to obtain the predictions. We used the fixed-point method with maximum number of available modes and a single iteration as we found that larger number of iterations does not improve the results. All other hyper-parameters were left as defaults as we found that changing them did not improve the results. We used the official code from Salvi et al. ().

#### C.4.2 Neural Spatiotemporal Processes

##### Constant Intensity Baseline.

As a simple baseline we use the following intensity function $\lambda(t,x)=c$, where $c$ is a learnable parameter that we fit on the training data.

##### DSTPP.

This model uses a transformer-based encoder to condition a diffusion-based density model on the event history. We use 500 time steps and 500 sampling steps, and train for at most 12 hours. All other hyper-parameters were left as defaults as we found that changing them did not improve the results or made the training too slow. We used the official code from Yuan et al. ().

##### NSTPP.

This model represents the vent history in terms of a latent state governed by a neural ODE, with the latent state being discontinuously updated after each event. The latent state is further used to represent the distribution of spatial locations via a normalizing flow model. We used the attentive CNF version of the model. All other hyper-parameters were left as defaults as we found that changing them did not improve the results. We used the official code from Chen et al. ().

##### AutoSTPP.

This model uses dual network approach for flexible and efficient modeling of spatiotemporal processes. We used 5 product networks with 2 layers 128 neurons each and hyperbolic tangent nonlinearities. The step size was set to 20 and number of steps was set to 101 in each direction. All other hyper-parameters were left as defaults as we found that changing them did not improve the results. We used the official code from Zhou and Yu ().

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: * Can you give me a possible explanation of why it works?
* There is no ablation study on why transformer are used for generating the initial states. Or do you have evidence the initial state part is robust to architecture choice?
* Despite the proposed speedup method, I believe neural-ODE is still untolerably slow and does not scale well. Do you have actual training/ inference time comparison?
soundness: 2
strengths: - The challenge seems well grounded as the sparse data over a large spatial domain is common for many types of problems, e.g., few agent trajectories over a large geographical domain.
- The method looks (possibly) scalable with low-resolution linear interpolation.
- The math formulation is clear and the empirical results are fair. 
- A lot of ablation study, accounting for context size, resolution, removal of components
summary: This is an engineering oriented work that model STPP with intensity driven by a continuous latent states governed by a Neural-ODE, with initial states generated by a transformer encoder. The formulation sounds valid and the proposed benefits are for sparsely distributed data. The main contributions are the new formulation and the interpolation-based speedup technique.
weaknesses: * I don't think the paper really answer the question of why it work on sparse data. There is no theoretical analysis / visualization of how the low-dimensional latent space captures the full dynamics from sparse observations. No discussion of information-theoretic bounds on what can be learned from sparse observations. It is reasonable to expect normalizing-flow based method (like Neural-STPP) not working well because the distribution is too localized, but I don't see why your method have an theoretical advantage over SOTA with kernel-based or closed spatial distribution.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: My main questions relate the points raised in the weaknesses above.
soundness: 3
strengths: -- The problem of learning spatiotemporal point processes is rather important, and any contribution to this problem should be well welcomed by the scientific community.
-- The overall idea of the article is meaningful. 
-- The numerical results are rather good.
summary: This article is devoted to the problem of learning spatiotemporal dynamics from randomly sampled points in space and time. This problem is particularly well suited for the situation where we have sensors that record a system, and we have to predict also the behavior of the sensors during the dynamics (e.g. meteorological sensors that are carried by currents). The method proposed in this article is based on the use of neural ODEs in a learned latent space.
weaknesses: -- Some explanations are not properly given. For instance, I assume that they are using an ODE solver in a latent space because a direct approach would immediately incur into stiffness problems. Why not using a neural PDE solver? Why is it better to learn a latent space and use an ODE solver for a problem that is formulated as a PDE (as in Eqn 5 of the paper)? This is unclear.
-- The latent approach makes the approach less clear, and more out of the control of the user. I suppose the authors have no idea why the encoder creates a certain latent space rather than another. A theoretical approach seems very complicated, in fact the authors limit themselves mostly to empirical results. 
-- It is unclear if a general system can be learned in this way. In a sense, we might think of the encoded latent space as a low-degree approximation of the system, but it might be that certain PDE models stemming from Eqn 5 might not be suitably tackled by such approach. 
-- One of the main claims is that the model is continuous. An interpolation task should be performed in this case to show that they can handle continuity well. They use interpolation in the method, but it is unclear if in an experiment where portions of the trajectories are completely hidden during training, could be recovered during evaluation.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: 1. For what reason is the distribution of the next sensor signal's location predicted? What is the benefit of such a prediction and what computational cost does it impose? If I understand correctly, Table 2 suggests removing the point process model (which simulates the next sensor signal position and time, if I'm correct). At least according to a minimal model when following Occams Razor.
2. The interpolation ablation is very illustrative. Have you tried higher-order interpolations to infer $\hat{z}(t_i)$, i.e., quadratic, cubic? What is the error that incurs from the interpolation compared to modeling the full temporal grid $t_1, \dots, t_N$? Table 1 demonstrates the time improvement when using interpolations; it would be great to see the error associated with the two techniques (Interp. vs Seq.).
3. Have you explored other solvers beyond dopri5, such as Euler, which is much cheaper? Or does the method depend on an adaptive solver to account for potentially different deltas between time steps? Figure 2 somehow suggest that the effectively processed time steps $\tau_m$ are separated by a constant time delta. Is this a requirement of the architecture?
4. How does the latent space dimensionality $d_z$ affect the runtime? Might be interesting to report along with its effect on the parameter count around line 375.
soundness: 4
strengths: _Originality:_ The combination of various different branches from machine learning is original. They are composed in an elegant and versatile way to solve spatiotemporal problems efficiently. The use of latent states enforces abstractions that push for generalizability. Intriguingly, the introduced method does not rely on complex architectures, but mostly emerges out of various MLPs that are well placed and wired.

_Quality:_ Claims are well supported with experimental results, which are contrasted against several recent  and competitive ML architectures. Figures are well designed and support the message conveyed in the manuscript.

_Clarity:_ The manuscript is well organized, structured, and written. A rich appendix provides details about the model design, yet a supplementary material to validate the results is missing.

_Significance:_ Results appear significant in terms of how sparse the data is sampled. Three synthetic problems of varying difficulty, as well as a real-world problem demonstrate the applicability of the method. Results are reported in two metrics along with standard deviations, which helps assessing the quality of the forecasts.
summary: A composition of different ML methods is presented to simulate spatiotemporal processes from point observations without access to the whole spatial field. The proposed approach encodes sparse context point observations into a latent state. In this latent state, the dynamic process evolution is integrated over large time steps, while a fine temporal resolution is obtained via interpolation in the latent state. A decoder projects the predicted latent states back to the observation space.
weaknesses: 1. Observation function is constrained to a normal distribution with fixed variance. It would be helpful to add arguments of this design choice, to what extend it limits the expressivity of the model, as well as for what problems this formulation is sufficient.
2. Ablations showing the performance under different spatial and temporal sparsities would be highly informative to understand the quality and limitations of the model at different tasks. Presumably, e.g., Navier-Stokes likely depends on more dense samples compared to Shallow Water. Extending this ablation to the other benchmarked methods would also provide valuable insights about the models' data efficacy.
3. Limitations are not reported. It would be valuable to understand the limits of the method, its computational cost, and the time this architecture needs to train. Also, it is unclear to which extend the method can directly be applied to a task at hand or how much fine tuning is involved.
4. No runtime comparison of the different models provided. If I'm not mistaken, the model must be called for each spatial position of interest in each time step, which amounts to a large number of model calls. Thus, to extend on Table 1, please provide information about the runtime of the entire model when generating a rollout of a spatiotemporal sequence of frames.
5. More details about the differences between the introduced method and AutoSTPP would be valuable, given that these two approaches perform almost equally well. For what reason is your method superior to AutoSTPP?

_Minor Comments_
- Typo in line 306, "withing"
- $N_{\text{ctx}}$ is unclear in Figure 4. What value does the variable take? Would be good to have the actual value. EDIT: C.1 provides this information; I thus suggest to refer to C.1 in the Caption of Figure 4.

---

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: Questions are related to the weaknesses:
Could you address the issue of interpretability and the Poisson process a bit more
soundness: 3
strengths: - The paper is well-written and technically sound. The methodology is clearly presented, and the experimental setup is detailed
- The proposed model is technically sound. It effectively combines techniques from various fields, including neural differential equations, neural point processes and amortized variational inference
- experiments and "ablations studies" are comprehensive, showing the impact of many parameters of the model
summary: This paper presents a novel method for modeling spatiotemporal dynamical systems from point process observations. The model integrates techniques from neural differential equations, neural point processes, implicit neural representations, and amortized variational inference. The authors also introduce a technique to speed training by addressing a computational bottleneck in latent state evaluation. The experimental results demonstrate the effectiveness of the model on challenging spatiotemporal datasets.
weaknesses: - While focusing on predictive capability and computational efficiency, discussing the interpretability of the model would enhance its value. Can something be said about the dynamical system?
- A little more discussion around the limitation of the Poisson process, and potential solution would have been welcome.

Meta-review: The paper proposes a method for learning to model spatio-temporal processes from data that is irregularly sampled in both the spatial and temporal dimensions. It employs an encode-process-decode framework and integrates the following components: point processes (non-homogeneous Poisson processes) for capturing irregular information and dynamics in a latent space, a neural ODE for time-stepping in the latent parameter space, variational inference for learning the parameters of the posterior distribution for the initial latent state, and Implicit neural representations for decoding. The model is evaluated on four datasets and compared against various baselines.

All reviewers acknowledge the novelty of the contribution, which combines existing components in an innovative way. They consider that the authors' claims are well-supported by the experiments and ablation studies. I recommend acceptance.

Decision: Accept (Spotlight)