Elucidating the Design Choice of Probability Paths in Flow Matching for Forecasting

## Abstract

Abstract Flow matching has recently emerged as a powerful paradigm for generative modeling and has been extended to probabilistic time series forecasting in latent spaces. However, the impact of the specific choice of probability path model on forecasting performance remains under-explored.
In this work, we demonstrate that forecasting spatio-temporal data with flow matching is highly sensitive to the selection of the probability path model.
Motivated by this insight, we propose a novel probability path model designed to improve forecasting performance.
Our empirical results across various dynamical system benchmarks show that our model achieves faster convergence during training and improved predictive performance compared to existing probability path models.
Importantly, our approach is efficient during inference, requiring only a few sampling steps.
This makes our proposed model practical for real-world applications and opens new avenues for probabilistic forecasting.

## 1 Introduction

Generative modeling has achieved remarkable success in recent years, especially for generating high-dimensional objects by learning mappings from simple, easily-sampled reference distributions, $\pi_{0}$, to complex target distributions, $\pi_{1}$.
In particular, diffusion models have pushed the capabilities of generating realistic samples across various data modalities, including images [, , ], videos [, , ], and spatio-temporal scientific data like climate and weather patterns [, ].
Despite their impressive performance, diffusion models often come with high computational costs during training and inference.
Additionally, they typically assume a Gaussian reference distribution, which may not be optimal for all data types and can limit modeling flexibility.

One promising alternative is flow matching, where the mappings are learned via a stochastic process that transforms $\pi_{0}$ into $\pi_{1}$ through a random ordinary differential equation (ODE), approximating its marginal vector flow [, , , , ].
While score-based generative models [, , , ] are specific instances of flow matching, with Gaussian transition densities, the general framework allows for a broader class of interpolating paths.
This flexibility can lead to deterministic sampling schemes that are faster and require fewer steps [].
Recent work has demonstrated the remarkable capabilities of flow matching models for generating high-dimensional images [] and discrete data [].

Building on this, flow matching in latent space has recently been applied to forecasting spatio-temporal data [] — predicting future frames in videos.
This approach leverages latent representations to capture the complex dynamics inherent in temporal data.
However, spatio-temporal forecasting, especially for video and dynamical systems data, presents unique challenges.
A video prediction model capable of generalizing to new, unseen scenarios must implicitly “understand” the scene: detecting and classifying objects, learning how they move and interact, estimating their 3D shapes and positions, and modeling the physical laws governing the environment [].
Similarly, accurate climate and weather forecasting requires capturing intricate physical processes and interactions across multiple scales [, ].

We observe that, in the context of spatio-temporal forecasting, the performance of flow matching is highly sensitive to the choice of the probability path model, an important topic which has not been widely explored within a unified framework.
Different probability paths can significantly impact the accuracy and convergence of forecasting models, particularly when dealing with complex dynamical systems characterized by partial differential equations (PDEs) and chaotic behaviors.
Motivated by this observation, we propose a novel probability path model specifically designed for probabilistic forecasting of dynamical systems.
Our choice of probability path results in flow matching models that can more effectively capture the underlying dynamics of spatio-temporal data, leading to more accurate and reliable forecasts.
Extending previous approaches, we provide a theoretical framework and efficient algorithms tailored for probabilistic forecasting using flow matching in latent space.
By investigating suitable choices of probability paths for modeling dynamical systems data, we provide insights into why our proposed model outperforms existing methods.
In particular, our probability path model leverages the inherent continuity and correlation in the spatio-temporal data to shorten the interpolating path when compared to other common choices of path, resulting in more stable training and faster convergence.

In summary, our proposed probability path model not only outperforms existing models on several forecasting tasks involving PDEs and other dynamical systems, but it also achieves faster convergence during training.
Moreover, we require only a few sampling steps during inference.
This efficiency enhances the practicality of flow matching approaches for real-world applications where computational resources and time are critical constraints.
Our main contributions are the following.

- •
Theoretical Framework and Efficient Algorithms:
We present a theoretical framework and efficient algorithms for applying flow matching in latent space to the probabilistic forecasting of dynamical systems (see Algorithms 1-2), extending the approach of [] and [].
Our approach is specifically tailored for time series data, enabling effective modeling of complex temporal dependencies inherent in dynamical systems.
- •
Novel Probability Path Model:
We propose a new probability path model (see Section 4.2) specifically designed for modeling dynamical systems data.
We provide intuitions to understand why our model leads to smoother training loss curve and faster convergence when compared to other models.
We also provide theoretical insights to show that the variance of the vector field (VF) generating our proposed path can be lower than that of the optimal transport VF proposed by [] for sufficiently correlated spatio-temporal samples (see Theorem 3 in App. C).
- •
Empirical Validation:
We provide extensive empirical results to demonstrate that our proposed probability path model can outperform other flow matching models on several forecasting tasks involving PDEs and other dynamical systems (see Section 6).
Our results demonstrate that the proposed probability path model outperforms existing flow matching models, achieving faster convergence during training and improved predictive performance.

Research code is available at [https://github.com/erichson/minFlowMatching/tree/main](https://github.com/erichson/minFlowMatching/tree/main).

## 2 Related Work

Modeling and forecasting dynamical systems. Traditionally dynamical systems arising in scientific applications have been studied through mathematical models based on physical laws, such as ODEs and PDEs []. These methods often require significant domain knowledge and strong assumptions, limiting their flexibility in data-driven scenarios where governing equations are unknown. Statistical and machine learning (ML) methods have emerged as powerful alternatives to study these systems. Classical approaches include autoregressive integrated moving average (ARIMA) and state-space models such as Kalman filters [], which have long been used for time series forecasting but struggle with highly nonlinear dynamics. ML methods such as recurrent neural networks (RNNs) [] and variants for spatio-temporal data [, ] are capable of learning complex temporal dependencies but they are challenging to train. Neural ODEs [, ] offer a structured approach to continuous-time modeling by integrating neural networks with differential equation solvers. However, these models still face challenges with generalizing to new systems, and learning continuous physical dynamics [].

Generative modeling for time series forecasting. Generative models, particularly score-based diffusion models and flow-based approaches, have gained significant attention for learning complex data distributions in high-dimensional spaces such as images [, , ], videos [, ], and dynamical systems [, , ]. In the context of time series forecasting, diffusion models like TimeGrad [, ] incorporate conditioning on previous time steps into both the forecasting process and the loss function. Building upon TimeGrad, CSDI [] enhances performance in imputation tasks by replacing the RNN encoder with a transformer. ScoreGrad [] adapts this framework to a score-based diffusion model for multivariate probabilistic forecasting. Methods such as DSPD and CSPD [], based on DDPM and SDE respectively, model dynamics as continuous functions and introduce time-correlated noise functions. Another relevant work, SSSD${}^{\text{S4}}$ [], uses state-space models (S4) to encode time series and performs diffusion over the temporal domain instead of across multivariate components. While these models demonstrate strong performance, they often require many sampling steps during generation. Related to our approach, [] proposes an SDE framework utilizing stochastic interpolants [] for probabilistic forecasting, but their sampler also requires many steps during generation. Another recent work is [], which introduces a flow matching framework for time series modeling within the data space; however, it concentrates on small ODE datasets and does not address forecasting tasks.

Applications of diffusion models to spatio-temporal data have primarily focused on video generation []. For spatio-temporal forecasting on sensor grids, graph neural networks (GNNs) are widely employed in tasks such as traffic prediction [, , ] and air quality forecasting []. Another recent work is DYffusion [], which utilizes a non-Markovian diffusion process to propagate images temporally, similar to DDIM. This method iteratively refines initial predictions at larger time steps by conditioning them on predictions at smaller time steps.

## 3 Flow Matching for Probabilistic Forecasting

In this section, we first introduce the objective of probabilistic forecasting; and then we discuss how flow matching can be used for learning conditional distributions in latent space.

Probabilistic forecasting framework.
Suppose that we are given a training set of $n$ trajectories, with each trajectory of length $m$, $S_{n}=\{(\mathbf{x}^{1:m})^{(i)}\}_{i=1,\dots,n}$, where $(\mathbf{x}^{1:m})^{(i)}=((x^{1})^{(i)},\dots,(x^{m})^{(i)})$, (with the $(x^{l})^{(i)}\in\mathbb{R}^{d}$), coming from an underlying continuous-time dynamical system.
For simplicity, we denote the trajectories as $\mathbf{x}^{1:m}=(x^{1},\dots,x^{m})$ unless there is a need to specify the corresponding $n$.
The trajectories are observed at arbitrary time points $\mathbf{t}_{1:m}=(t_{1},\dots,t_{m})$ such that $x^{i}:=x(t_{i})\in\mathbb{R}^{d}$ and $(x(t))_{t\in[t_{1},t_{m}]}$ are the observed states of the ground truth system.
In practice, we may have access to only few trajectories, i.e., $n$ is small or even $n=1$, and the trajectories themselves may be observed at different time stamps.

The goal of probabilistic forecasting is to predict the distribution of the upcoming $l$ elements given the first $k$ elements, where $m=l+k$:

$$ $q(x^{k+1},\dots,x^{k+l}|x^{1},\dots,x^{k})=\prod_{i=1}^{l}q(x^{k+i}|x^{1}, \dots,x^{k+i-1}).$ (1) $$

We propose to model each one-step predictive conditional distribution in Eq. (1) via a probability density path.
Instead of using score-based diffusion models to specify the path, we choose latent flow matching, a simpler method to train generative models.
With flow matching, we directly work with probability paths, and we can simply avoid reasoning about diffusion processes altogether.

Flow matching in latent space.
Let $z^{\tau}=\mathcal{E}(x^{\tau})$ for $\tau=1,\dots,m$, where $\mathcal{E}$ denotes a pre-trained encoder that maps from the data space to a lower dimensional latent space. Working in the latent space, our goal is to approximate the ground truth distribution $q(z^{\tau}|x^{1},\dots,x^{\tau-1})$ by the parametric distribution $p(z^{\tau}|z^{\tau-1})$, which can then be decoded as $x^{\tau}=\mathcal{D}(z^{\tau})$. The latent dynamics can be modeled by an ODE:

$$ $\dot{Z}_{t}=u_{t}(Z_{t}),$ (2) $$

where $u_{t}$ is the vector field describing the instantenous rate of change of the state at time $t$. Learning the dynamics of the system is equivalent to approximating the vector field $u_{t}$.
This can be done by regressing a neural network using the mean squared error (MSE) loss.

Following the idea of flow matching, we infer the dynamics of the system generating ${\bf z}$ from the collection of latent observables by learning a time-dependent vector field $v_{t}:[0,1]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$, $t\in[0,1]$, such that the ODE

$$ $\dot{\phi}_{t}(Z)=v_{t}(\phi_{t}(Z)),\ \ \phi_{0}(Z)=Z,$ (3) $$

defines a time-dependent diffeomorphic map (called a flow), $\phi_{t}(Z):[0,1]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$, that pushes a reference distribution $p_{0}(Z)$ towards the distribution $p_{1}(Z)\approx q(Z)$ along some probability density path $p_{t}(Z)$ and the corresponding vector field $u_{t}(Z)$.
In other words, $p_{t}=[\phi_{t}]_{*}p_{0}$, where $[\cdot]_{*}$ denotes the push-forward operation.
Here, $q$ is the ground truth distribution, $p$ denotes a probability density path, i.e., $p:[0,1]\times\mathbb{R}^{d}\to\mathbb{R}_{>0}$, and $\int p_{t}(Z)dZ=1$.
We also write $Z_{t}=\phi(Z)$; and thus the ODE can be written as $\dot{Z}_{t}=v_{t}(Z_{t}),Z_{0}=Z$.
Typically the reference distribution $p_{0}$ is chosen to be the standard Gaussian [, ].

In other words, the main goal of flow matching is to learn a deterministic coupling between $p$ and $q$ by learning a vector field $v_{t}$ such that the solution to the ODE (3) satisfies $Z_{0}\sim p$ and $Z_{1}\sim q$. When ${\bf Z}=(Z_{t})_{t\in[0,1]}$ solves Eq. (3) for a given function $v_{t}$, we say that ${\bf Z}$ is a flow with the vector field $v_{t}$.
If we have such a vector field, then $(Z_{0},Z_{1})$ is a coupling of $(p,q)$.
If we can sample from $p$, then we can generate approximate samples from the coupling by sampling $Z_{0}\sim p$ and numerically integrating Eq. (3). This can be viewed as a continuous normalizing flow [].

If one were given a pre-defined probability path $p_{t}(Z)$ and the corresponding vector field $u_{t}(Z)$ that generates the path, then one could parametrize $v_{t}(Z)$ with a neural network $v^{\theta}_{t}(Z)$, with $\theta$ the learnable parameter, and solve the least square regression by minimizing the flow matching loss:

$$ $\min_{\theta}\mathcal{L}_{fm}(\theta):=\mathbb{E}_{t,p_{t}(Z)}\omega(t)\|v^{ \theta}_{t}(Z)-u_{t}(Z)\|^{2},$ (4) $$

where $t\in\mathcal{U}[0,1]$, $Z\sim p_{t}(Z)$ and $\omega(t)>0$ is a weighting function. We take $\omega(t)=1$ unless specified otherwise.
However, we do not have prior knowledge for choosing $p_{t}$ and $u_{t}$, and there are many choices that can satisfy $p_{1}\approx q$. Moreover, we do not have access to a closed form $u_{t}$ that generates the desired $p_{t}$. We shall follow the approach of [] and construct a target probability path by mixing simpler conditional probability paths. This probability path is the marginal probability path:

$$ $p_{t}(Z)=\int p_{t}(Z|\tilde{Z})q(\tilde{Z})d\tilde{Z},$ (5) $$

obtained by marginalizing the conditional probability density paths $p_{t}(Z|\tilde{Z})$ over observed latent trajectories $\tilde{Z}$, with $p_{0}(Z|\tilde{Z})=p(Z)$ and $p_{1}(Z|\tilde{Z})=\mathcal{N}(Z|\tilde{Z},\epsilon^{2}I)$ for a small $\epsilon>0$. Doing so gives us a marginal probability $p_{1}$ which is a mixture distribution that closely approximates $q$. Then, assuming that $p_{t}(Z)>0$ for all $Z$ and $t$, we can also define a marginal vector field as:

$$ $u_{t}(Z)=\int u_{t}(Z|\tilde{Z})\frac{p_{t}(Z|\tilde{Z})q(\tilde{Z})}{p_{t}(Z) }d\tilde{Z},$ (6) $$

where $u_{t}(Z|\tilde{Z})$ is a conditional vector field (conditioned on the latent trajectory $\tilde{Z}$).
It turns out that this way of mixing the conditional vector fields leads to the correct vector field for generating the marginal probability path (5). We can then break down the intractable marginal VF into simpler conditional VFs which depends on a single sample.

To deal with the intractable integrals in Eq. (5)-(6) which complicates computation of an unbiased estimator of $\mathcal{L}_{fm}$, we shall minimize the conditional loss proposed by []:

$$ $\min_{\theta}\mathcal{L}_{cfm}(\theta):=\mathbb{E}_{t,p_{t}(Z|\tilde{Z}),q( \tilde{Z})}\omega(t)\|v^{\theta}_{t}(Z)-u_{t}(Z|\tilde{Z})\|^{2},$ (7) $$

where $t\in\mathcal{U}[0,1]$, $\tilde{Z}\sim q(\tilde{Z})$, $Z\sim p_{t}(Z|\tilde{Z})$ and $u_{t}(Z|\tilde{Z})$ is the vector field defined per sample $\tilde{Z}$ that generates the conditional probability path $p_{t}(Z|\tilde{Z})$. Importantly, one can show that the solution of (7) is guaranteed to converge to the same result in (4); see Theorem 1 in App. C. Therefore, the conditional flow matching loss can match the pre-defined target probability path, constructing the flow that pushes $p_{0}$ towards $p_{1}$. Since both the probability path and VF are defined per sample, we can sample unbiased estimates of the conditional loss efficiently, particularly so with suitable choices of conditional probability paths and VFs.

## 4 Probability Path Models for Probabilistic Forecasting

In this section, we describe the family of probability paths that we consider for flow matching, and we propose an improved model for probabilistic forecasting of spatio-temporal data.

### 4.1 Common Probability Path Models

The family of Gaussian conditional probability paths gives us tractable choices to work with since the relevant quantities in Eq. (7) and thus the conditional flow can be defined explicitly. Therefore, we will work with Gaussian probability paths. Moreover, we are going to solve (7) over the dataset of all transition pairs $\mathcal{D}=\{(z^{\tau-1},z^{\tau})\}_{\tau=2,\dots,m}$, and use a pair of points for $\tilde{Z}$, setting $\tilde{Z}=(Z_{0},Z_{1})\in\mathcal{D}$.
In particular, we consider the following class of models for the probability path:

$$ $p_{t}(Z|\tilde{Z}:=(Z_{0},Z_{1}))=\mathcal{N}(Z|a_{t}Z_{0}+b_{t}Z_{1},c^{2}_{t }I),$ (8) $$

where $a_{t}$, $b_{t}$ and $c_{t}$ are differentiable time-dependent functions on $[0,1]$, and $I$ denotes the identity.

Table 1 provides five different choices of probability paths, including our proposed choice, that we study here.
The optimal transport (OT) VF model described in Table 1 was initially proposed by [], and setting $\epsilon_{min}=0$ gives us the rectified flow model of [], which proposed connecting data and noise on a straight line.
The stochastic interpolant model in Table 1 is the one considered by [].
The VE and VP-diffusion conditional VFs (derived with Theorem 2) coincide with the vector fields governing the Probability Flow ODE for the VE and VP diffusion paths proposed in [].
It has been shown that combining diffusion conditional VFs with the flow matching objective leads to a training alternative that is more stable and robust when compared to existing score matching approaches [].

**Table 1: Choices of probability density paths that we study in this paper.**
| Model | $a_{t}$ | $b_{t}$ | $c_{t}^{2}$ |
| --- | --- | --- | --- |
| VE-diffusion | 1 | 0 | $\sigma_{1-t}^{2}$, where $\sigma_{t}$ is increasing in $t$, $\sigma_{0}=0$ |
| VP-diffusion | $e^{-\frac{1}{2}T(1-t)}$ | 0 | $1-e^{-T(1-t)},T(t)=\int_{0}^{t}\beta(s)ds$, $\beta=$ noise scale |
| OT-VF / rectified flow | $t$ | 0 | $(1-(1-\epsilon_{min})t)^{2}$, $\epsilon_{min}\geq 0$ |
| Stochastic interpolant | $1-t$ | $t$ or $t^{2}$ | $\epsilon^{2}t(1-t)^{2}$, $\epsilon>0$ |
| Ours | $1-t$ | $t$ | $\sigma_{min}^{2}+\sigma^{2}t(1-t)$, $\sigma_{min},\sigma\geq 0$ |

As remarked in [], there are many choices of vector fields that generate a given probability path. We shall use the simplest vector field that generates flow whose map is affine linear. Let $p_{t}(Z|\tilde{Z})$ be the Gaussian probability path defined in Eq. (8) and consider the flow map $\psi_{t}$ defined as $\psi_{t}(Z)=a_{t}Z_{0}+b_{t}Z_{1}+c_{t}Z$ with $c_{t}>0$. Then the unique vector field that defines $\psi_{t}$ is (see Theorem 2 and the proof in App. C):

$$ $u_{t}(Z|\tilde{Z})=\frac{c_{t}^{\prime}}{c_{t}}(Z-(a_{t}Z_{0}+b_{t}Z_{1}))+a_{ t}^{\prime}Z_{0}+b_{t}^{\prime}Z_{1},$ (9) $$

where prime denotes derivative with respect to $t$, and $u_{t}(Z|\tilde{Z})$ generates the Gaussian path $p_{t}(Z|\tilde{Z})$.

In view of this, minimizing the conditional loss becomes:

$$ $\min_{\theta}\mathcal{L}_{cfm}(\theta):=\mathbb{E}_{t,z^{\tau},z^{\tau-1},p_{t }(Z|z^{\tau},z^{\tau-1})}\omega(t)\bigg{\|}v^{\theta}_{t}(Z)-\frac{c_{t}^{ \prime}}{c_{t}}(Z-(a_{t}z^{\tau}+b_{t}z^{\tau-1}))-a_{t}^{\prime}z^{\tau}-b_{t }^{\prime}z^{\tau-1}\bigg{\|}^{2},$ $$

where $t\sim\mathcal{U}[0,1]$, $Z\sim p_{t}(Z|z^{\tau},z^{\tau-1})$ and ${\bf z}\sim q({\bf z})$.
We refer to this as the Flow Matching loss parametrization and work with this parametrization. There are other parametrizations: most popular ones are the Score Matching loss, Score Flow loss and DDPM loss. See App. B for a comparison of different loss parametrizations and App. A for connections to SDE based generative models.

### 4.2 A Novel Probability Path Model

We propose to choose $a_{t}=1-t$, $b_{t}=t$, $c^{2}_{t}=\sigma_{min}^{2}+\sigma^{2}t(1-t)$, in which case we have the probability path described by:

$$ $p_{t}(Z|\tilde{Z})=\mathcal{N}(Z|tZ_{1}+(1-t)Z_{0},(\sigma_{min}^{2}+\sigma^{2 }t(1-t))I),$ (10) $$

which transports a Gaussian distribution centered around $Z_{0}$ with variance $\sigma_{min}^{2}$ at $t=0$ to a Gaussian distribution centered around $Z_{1}$ at $t=1$ with variance $\sigma_{min}^{2}$. Here $\sigma_{min},\sigma\geq 0$ are tunable parameters.
In the case when $\sigma_{min}=0$, it describes a Brownian bridge that interpolates between $Z_{0}$ and $Z_{1}$ []. To ensure numerical stability when sampling $t\sim\mathcal{U}[0,1]$, it is beneficial to use a small $\sigma_{min}>0$.
Note that $\sigma^{2}$ is a scale factor determining the magnitude of fluctuations around the path interpolating between $Z_{0}$ and $Z_{1}$. The variance $c_{t}^{2}$ is minimum with the value of $\sigma_{min}^{2}$ at the endpoints $t=0$ and $t=1$, and the maximum variance is $\sigma_{min}^{2}+\sigma^{2}/4$ which occurs in the middle of the path at $t=1/2$.

The corresponding vector field that defines the flow is then given by (applying Theorem 2):

$$ $u_{t}(Z|\tilde{Z})=Z_{1}-Z_{0}+\frac{\sigma^{2}}{2}\frac{1-2t}{\sigma_{min}^{2 }+\sigma^{2}t(1-t)}(Z-(tZ_{1}+(1-t)Z_{0})).$ (11) $$

We remark that the conditional vector field $u_{t}(Z|\tilde{Z})$ is linear in $Z$ and $\tilde{Z}$.
This choice of the vector field corresponds to the choice of a linear model for the dynamics of the latent variable $Z$, and thus is reminiscent of the Koopman formalism [].

We expect that our proposed probability path model improves upon the other considered models, as it takes advantage of the inherent continuity and correlation in the spatio-temporal data.
Intuitively, for time series samples whose underlying dynamics are continuous and obey a physical law, a Gaussian sample is typically further from the time series samples, so the distance between a Gaussian sample and a time series sample should generally be larger than the distance between consecutive time series samples (which can be highly correlated).
Therefore, using a probability path that connects consecutive time series samples could lead to faster convergence and more stable training, when compared to using a path that simply connects a time series sample to a Gaussian sample, since the resulting flow model in the former case can better capture the true system dynamics with less effort.

Moreover, if the consecutive samples are sufficiently correlated, then the variance of the vector field corresponding to our proposed probability path model can be lower than the variance of the vector field corresponding to the other choices of probability paths. We refer to Theorem 3 in App. C.3 for such comparison result for our proposed model and the optimal transport VF model of [] (see also the discussions in Section C.3).

## 5 An Efficient Probabilistic Forecasting Algorithm

In this section, we present efficient algorithms for training and inferencing the flow matching model.

Recently [] proposed an efficient algorithm for latent flow matching for the task of video prediction, using the probability path generated by the optimal transport VF of []. To enable efficient training, we shall follow [] and leverage the
iterative nature of sampling from the learned flow and use a single random conditioning element from the past at each iteration. However, our method differs from [] as we shall use different probability paths and target vector fields.

Training.
We set $Z_{1}$ to be the target element and $Z_{0}$ to be the reference element chosen to be the previous element before the target element, i.e., if $Z_{1}=z^{\tau}$, then $Z_{0}=z^{\tau-1}$, for $\tau=2,\dots,m$. Note that this differs from [], where $\tilde{Z}=z_{1}$ (i.e., they do not use a reference element to define their probability path, whereas we use a pair of elements $(Z_{0},Z_{1})$). In this way, our probability path model maps a distribution centered around a previous state to the distribution of the current states, which is more natural from the point of view of probabilistic forecasting whose goal is to obtain an ensemble of forecasts. Algorithm 1 summarizes the training procedure of our method.

Both the autoencoder and the VF neural network can also be jointly trained in an end-to-end manner, but our results show that separating the training can lead to improved performance. Moreover, doing so allows us to better assess the impact of using different probability paths.

Figure: Algorithm 1 Flow matching for spatio-temporal data

Figure: Algorithm 2 One-step ahead forecasting with forward Euler

Inference.
We use an ODE sampler during inference to generate forecasts.
The ODE sampler is described as follows.
Let $(Y_{i}^{\tau})_{i=0,\dots,N}$ denote the generation process, where $N$ is the number of integration steps and the superscript $\tau$ denotes the time index for which the generation/forecast is intended for. Given the previous elements $(x^{1},\dots,x^{T-1})$ of a time series sample, in order to generate the next element (i.e., the $T$-th element), we start with sampling the initial condition $Y_{0}^{T}$ from $\mathcal{N}(z^{T-1},\sigma_{sam}^{2}I)$ for some small $\sigma_{sam}\geq 0$, where $z^{T-1}=\mathcal{E}(x^{T-1})$. This is in contrast to the procedure of [], which simply uses a mean-zero Gaussian sample instead. We then use an ODE solver to integrate the learned vector field along the time interval $[0,1]$ to obtain $Y_{N}^{T}$. We use $\mathcal{D}(Y_{N}^{T})$ as an estimate of $x^{T}$, and forecasting is done autoregressively.
Algorithm 2 summarizes this procedure when the sampling is done using the forward Euler scheme. Note that we can also use computationally more expensive numerical schemes such as the Runge-Kutta (RK) schemes.

## 6 Empirical Results

In this section, we present our main empirical results to elucidate the design choice of probability paths.
We focus on PDE dynamics forecasting tasks here (additional results can be found in App. D).
We test the performance of our probability path model, i.e., Eq. (10) with $a_{t}=1-t$, $b_{t}=t$ and $c_{t}=\sqrt{\sigma_{min}^{2}+\sigma^{2}t(1-t)}$ on these tasks. We pick $\sigma_{min}=0.001$, and treat $\sigma$ and $\sigma_{sam}$ as tunable parameters.
We compare our proposed model with four other models of probability paths:

- •
RIVER []: RIVER uses the OT-VF model in Table 1, i.e., $a_{t}=0$, $b_{t}=t$, $c_{t}=1-(1-\epsilon_{min})t$, choosing $\epsilon_{min}=10^{-7}$.
- •
VE-diffusion in Table 1: We use $\sigma_{t}=\sigma_{min}\sqrt{\left(\frac{\sigma_{max}}{\sigma_{min}}\right)^{2
t}-1}$ with $\sigma_{min}=0.01$, $\sigma_{max}=0.1$, and sample $t$ uniformly from $[0,1-\epsilon]$ with $\epsilon=10^{-5}$ following [].
- •
VP-diffusion in Table 1: We use $\beta(s)=\beta_{min}+s(\beta_{max}-\beta_{min})$ where $\beta_{min}=0.1,\beta_{max}=20$ and $t$ is sampled from $\mathcal{U}[0,1-\epsilon]$ with $\epsilon=10^{-5}$, following []. Thus, $T(s)=s\beta_{min}+\frac{1}{2}s^{2}(\beta_{max}-\beta_{min})$.
- •
The stochastic interpolant path in Table 1: We consider the path proposed by [] and use the suggested choice of $a_{t}=1-t$, $b_{t}=t^{2}$ and $c_{t}=\epsilon(1-t)\sqrt{t}$ (see Eq. (2) in [] and note that $Var((1-t)W_{t})=(1-t)^{2}t$ for the standard Wiener process $W_{t}$).
We choose $\epsilon=0.01$. This is a path that is similar to ours, but with the variance $c_{t}^{2}$ chosen such that the maximum occurs at $t=1/\sqrt{3}$ instead of at the middle of the path at $t=1/2$.
Later, we will see that different forms of variance can lead to vastly different performance in the considered tasks.

Evaluation metrics. We evaluate the models using the following metrics.
First, we use the standard mean squared error (MSE) and the relative Frobenius norm error (RFNE) to measure the difference between predicted and true snapshots.
Second, we compute metrics such as the peak signal-to-noise ratio (PSNR), and the structural similarity index measure (SSIM) to further quantify the quality and similarity of the generated snapshots [].
Third, we use the Pearson correlation coefficient to assess the correlation between predicted and true snapshots.

Training details.
We use an autoencoder (AE) to embed the training data into a low-dimensional latent space, which enables the model to capture the most relevant features of the data while reducing dimensionality [].
We then train a flow matching model in this latent space.
Training generative models in latent space has also been considered by [] for score matching models and by [] for flow matching models.
To train the AE, we minimize a loss function that consists of reconstruction error, in terms of MSE, between the input data and its reconstructed version from the latent space.
The choice of the AE architecture is tailored to the complexity of the dataset (see App. E for details).

### 6.1 Probabilistic Forecasting of Dynamical Systems

We evaluate the performance of our proposed probability path model on challenging dynamical systems to demonstrate its effectiveness in forecasting complex continuous dynamics. Specifically, we consider the following tasks (for details see App. E.1):

- •
Fluid Flow Past a Cylinder: This task involves forecasting the vorticity of a fluid flowing past a cylinder. The model conditions on the first 5 frames and predicts the subsequent 20 frames at a resolution of $64\times 64$ with 1 channel representing vorticity.
- •
Shallow-Water Equation: This dataset models the dynamics of shallow-water equations [], capturing essential aspects of geophysical fluid flows. We use the first 5 frames for conditioning and predict the next 15 frames at a resolution of $128\times 128$ with 1 channel representing horizontal flow velocity.
- •
Diffusion-Reaction Equation: This dataset models the dynamics of a 2D diffusion-reaction equation [].
We use the first 5 frames for conditioning and predict 15 future frames at a 128 $\times$ 128 with 2 channels representing velocity in the $x$ and $y$ directions.
- •
Incompressible Navier-Stokes Equation: As a more challenging benchmark, we consider forecasting the dynamics of a 2D incompressible Navier-Stokes equation [].
We use the first 5 frames for conditioning and predict the next 20 frames at a resolution of $512\times 512$ with 2 channels representing velocity in the $x$ and $y$ directions.

Table 2 summarizes the performance of our model compared to other probability path models across all tasks. It can be seen that our probability path model achieves the lowest test MSE and RFNE across all tasks, indicating more accurate forecasts. Moreover, the higher PSNR and SSIM scores indicate that our model better preserves spatial structures in the predictions.
Despite the similarity of our proposed model with the stochastic interpolant of [], in that both models use consecutive samples to define the path, our model outperforms the stochastic interpolant model, suggesting that choosing the maximum variance to occur at the middle of the path is a better choice.
Importantly, our model is highly efficient during inference time since it requires only 10 sampling steps; this is significantly fewer than the 50 steps needed by other models.

Figure 1 shows the Pearson correlation coefficients of the predicted snapshots over time for all models. Our model’s predictions shows a slower decay of correlation coefficients compared to other models, indicating better temporal consistency and long-term predictive capability. Correlation coefficients about $95\%$ indicate performance on par with physics-based numerical simulators.

**Table 2: Results for forecasting dynamical systems using different probability path models for flow matching. Results are averaged over 5 generations obtained with 10 sampling steps using RK4.**
| Task | Model | Test MSE ($\downarrow$) | Test RFNE ($\downarrow$) | PSNR ($\uparrow$) | SSIM ($\uparrow$) |  |
| --- | --- | --- | --- | --- | --- | --- |
| Flow past Cylinder | RIVER | 3.05e-03 | 5.70e-02 | 42.66 | 0.98 |  |
| VE-diffusion | 2.75e-01 | 5.21e-01 | 26.81 | 0.52 |  |  |
| VP-diffusion | 3.20e-03 | 5.20e-01 | 42.60 | 0.97 |  |  |
| Stochastic interpolant | 3.39e-03 | 6.09e-02 | 41.82 | 0.97 |  |  |
| Ours | 3.80e-04 | 2.30e-02 | 48.89 | 0.99 |  |  |
| Shallow-Water | RIVER | 9.29e-04 | 1.5e-01 | 34.90 | 0.91 |  |
| VE-diffusion | 1.23e-02 | 5.5e-01 | 28.76 | 0.57 |  |  |
| VP-diffusion | 1.31e-03 | 1.8e-01 | 34.33 | 0.88 |  |  |
| Stochastic interpolant | 1.06e-03 | 1.5e-01 | 35.57 | 0.88 |  |  |
| Ours | 6.90e-04 | 1.3e-01 | 36.10 | 0.93 |  |  |
| Diffusion-Reaction | RIVER | 2.37e-03 | 2.07e-01 | 38.47 | 0.84 |  |
| VE-diffusion | 9.74e-02 | 1.7 | 32.83 | 0.35 |  |  |
| VP-diffusion | 1.72e-02 | 6.48e-01 | 34.65 | 0.53 |  |  |
| Stochastic interpolant | 6.17e-02 | 8.62e-01 | 45.64 | 0.76 |  |  |
| Ours | 3.56e-04 | 1.16e-01 | 34.34 | 0.89 |  |  |
| Navier-Stokes | RIVER | 1.98e-04 | 2.5e-01 | 38.71 | 0.85 |  |
| VE-diffusion | 1.21e-04 | 2.4 | 29.53 | 0.27 |  |  |
| VP-diffusion | 7.62e-02 | 1.85e-01 | 40.73 | 0.93 |  |  |
| Stochastic interpolant | 7.26e-05 | 1.8e-01 | 37.81 | 0.88 |  |  |
| Ours | 1.90e-05 | 9.8e-02 | 41.52 | 0.96 |  |  |

Figure: Figure 1: Pearson correlation coefficient to assess the correlation between predicted and true snapshots at various prediction steps for different probability path models. Our probability path model shows the best performance on all three tasks. Note that the first 5 snapshots are the conditioning snapshots.

Figure 2 compares the training loss curves of our model with others trained on the fluid flow past a cylinder and the shallow water equation task. Our method leads to faster convergence during training and smoother loss curves.
This suggests that our model requires fewer iterations to generate high-quality samples when compared to other flow matching models.

### 6.2 Ablation Study

To further assess our model, we conducted an ablation study focusing on the impact of various hyperparameters. Specifically, we study the impact of the values of $\sigma$, the choice of sampler, and the number of sampling steps during inference. For small $\sigma_{sam}$, we find that this parameter has negligible impact on test performance, so we fixed $\sigma_{\text{sam}}=0$ for all experiments in this section.

Impact of $\sigma$ on training stability. Figure 3 illustrates the effect of different $\sigma$ values on the training loss curve for our method on the fluid flow past a cylinder task. We observed that larger values of $\sigma$ (e.g., $\sigma=0.1$) resulted in smoother loss curves and more stable convergence during training.

Effect of $\sigma$, sampler choice, and sampling steps on test performance. Table 3 investigates how different values of $\sigma$, the choice of sampler (Euler or RK4), and the number of sampling steps affect test performance. It can be seen, that even with as few as 5 sampling steps using the Euler scheme, our model perform reasonably well. However, increasing the number of sampling steps or employing the more computationally intensive RK4 sampler can help to lead to better results.

Figure: Figure 2: Training loss for different models of probability path for the fluid flow past a cylinder task. Our model leads to fastest convergence and smoothest loss curve among all models.

Figure: Figure 3: Training loss for different values of $\sigma$ using our probability path model. The left subplot shows results for the Euler sampler, and the right subplot for the RK4 sampler. We see that the loss curve is sensitive to the choice of $\sigma$, with larger values of $\sigma$ giving smoother loss curves.

**Table 3: Ablation study for the fluid flow past a cylinder task. Results are averaged over 5 generations.**
| $\sigma$ | sampler | sampling steps | Test MSE ($\downarrow$) | Test RFNE ($\downarrow$) | PSNR ($\uparrow$) | SSIM ($\uparrow$) |
| --- | --- | --- | --- | --- | --- | --- |
| 0.0 | Euler | 5 | 1.78e-03 | 4.50e-02 | 44.12 | 0.98 |
| 0.01 | Euler | 5 | 1.79e-03 | 4.50e-02 | 44.13 | 0.98 |
| 0.1 | Euler | 5 | 1.94e-03 | 4.50e-02 | 44.26 | 0.98 |
| 0.0 | Euler | 10 | 7.14e-04 | 2.90e-02 | 47.07 | 0.99 |
| 0.01 | Euler | 10 | 3.81e-04 | 2.30e-02 | 48.73 | 0.99 |
| 0.1 | Euler | 10 | 2.83e-03 | 5.10e-02 | 44.14 | 0.98 |
| 0.0 | RK4 | 10 | 3.89e-04 | 2.20e-02 | 49.22 | 0.99 |
| 0.01 | RK4 | 10 | 3.80e-04 | 2.30e-02 | 48.89 | 0.99 |
| 0.1 | RK4 | 10 | 6.49e-03 | 7.80e-02 | 40.70 | 0.97 |
| 0.0 | Euler | 20 | 5.51e-04 | 2.70e-02 | 47.44 | 0.99 |
| 0.01 | Euler | 20 | 7.32e-04 | 3.10e-02 | 46.85 | 0.99 |
| 0.1 | Euler | 20 | 7.93e-04 | 3.10e-02 | 46.80 | 0.99 |
| 0.0 | RK4 | 20 | 9.31e-04 | 3.30e-02 | 46.42 | 0.99 |
| 0.01 | RK4 | 20 | 6.62e-04 | 2.80e-02 | 47.62 | 0.99 |
| 0.1 | RK4 | 20 | 6.58e-04 | 2.90e-02 | 47.14 | 0.99 |

### 6.3 Discussion

Our empirical results showed that our proposed model consistently outperformed other models across different forecasting tasks involving different types of dynamical systems. Our model shows improved training efficiency, with faster convergence reducing the computational resources and time required for model training. Moreover, our model is efficient during inference time since it only requires a few sampling steps, making it practical for real-world applications where computational efficiency is crucial. Additionally, the model maintained better temporal consistency, as indicated by a slower decay of Pearson correlation coefficients over longer prediction horizons. These findings validate the effectiveness of our approach in modeling complex dynamical systems.

Our ablation study further validated the advantage of our proposed probability path model. We found that larger $\sigma$ values not only contributed to smoother training loss curves but also enhanced the overall stability and efficiency of the model. Moreover, we saw that $\sigma_{\text{sam}}=0$ can be fixed without compromising accuracy. The fact that our model achieved improved performance even with the simplest sampler (Euler scheme) and a minimal number of sampling steps (as few as five) validates its practical applicability, especially in scenarios where computational resources and time are limited.

## 7 Conclusion

In this work, we investigated the use of flow matching in latent space for probabilistic forecasting of spatio-temporal dynamics, providing a theoretical framework and an efficient algorithm. We demystified the critical role of the probability path design in this setting and proposed an improved probability path model. Our model leverages the inherent continuity and correlation in the spatio-temporal data, leading to more stable training and faster convergence. Our empirical evaluations on several dynamical systems forecasting tasks demonstrated that our model performs favorably when compared to existing models. These findings highlight that while flow matching techniques hold great promise for probabilistic forecasting of spatio-temporal dynamics, it is important to select an appropriate probability path model to achieve optimal test performance.

#### Acknowledgments

The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725 (NAISS 2024/5-269).

SHL would like to acknowledge support from the Wallenberg Initiative on Networks and Quantum Information (WINQ) and the Swedish Research Council (VR/2021-03648).
AY was supported by the SciAI Center, funded by the Office of Naval Research under Grant Number N00014-23-1-2729. NBE would like to acknowledge LBL’s LDRD initiative for providing partial support.

## Appendix

This appendix is organized as follows.
In App. A, we provide some remarks on the connection of flow matching models to other generative models.
In App. B, we provide commonly used loss parametrizations and compare them with our flow matching loss.
In App. C, we provide theorems and proofs to justify the discussions in Section 3 (see Theorem 1-3).
In App. D, we provide additional empirical results.
In App. E, we provide the missing experimental details.

## Appendix A Connection to SDE Based Generative Model Frameworks

In this section, we provide some remarks on the connection between flow matching models and SDE based generative models [].

Consider the following continuous-time Gaussian latent variable model []:
$Z_{t}=\mathcal{E}(X_{t})$,

$$ $\displaystyle Z_{t}$ $\displaystyle=A_{t}Z_{0}+B_{t}Z_{1}+L_{t}\epsilon,\ \ t\in[0,1],$ (13) $\displaystyle X_{t}$ $\displaystyle=\mathcal{D}(Z_{t}),$ (14) $$

where $t$ is the continuous variable, $X_{0}\in\mathbb{R}^{d}$ represent data samples, $Z\in\mathbb{R}^{p}$ is the latent variable, $\epsilon\sim\mathcal{N}(0,I)$ is independent of the random variables $X_{0},Z_{0},Z_{1}$. Here $A_{t},B_{t}$ and $L_{t}\geq 0$ are pre-specified coefficients which are possibly matrix-valued and time-dependent, $\mathcal{D}$ and $\mathcal{E}$ denote the decoder and encoder map respectively, and $\mathcal{D}\circ\mathcal{E}=I$. Note that $Z_{0}$
and $Z_{1}$ are initial and terminal point of the path $(Z_{t})_{t\in[0,1]}$ in the latent space.

The above latent variable model can be identified (up to equivalence in law for each $t$) with the linear SDE of the form:

$$ $d\hat{Z}_{t}=F_{t}\hat{Z}_{t}dt+H_{t}Z_{1}dt+G_{t}dW_{t},\ \ \hat{Z}_{0}=Z_{0} ,\ \ t\in[0,1),$ (15) $$

where $(W_{t})_{t\in[0,1]}$ is the standard Wiener process.
By matching the moments, we obtain

$$ $\displaystyle F_{t}$ $\displaystyle=\dot{A}_{t}A_{t}^{-1},$ (16) $\displaystyle H_{t}$ $\displaystyle=\dot{B}_{t}-\dot{A}_{t}A_{t}^{-1}B_{t},$ (17) $\displaystyle G_{t}G_{t}^{T}$ $\displaystyle=\dot{L}_{t}L_{t}^{T}+L_{t}\dot{L}_{t}^{T}-\dot{A}_{t}A_{t}^{-1}L _{t}L_{t}^{T}-L_{t}L_{t}^{T}A_{t}^{-T}(\dot{A}_{t})^{T},$ (18) $$

where the overdot denotes derivative with respect to $t$ and $A^{T}$ denotes the transpose of $A$.

Under the above formulation, various existing generative models such as DDPM [], VP-SDE and VE-SDE of
[, ], the critically damped SDE of [], the flow matching models in
[, , ] and the stochastic interpolants of [] can be recovered, and new models can be derived.

The following proposition establishes
the connection between flow matching using our proposed probability path model, the Gaussian latent variable model (13) and the linear SDE model (15).

###### Proposition 1 .

For every $t\in[0,1]$, the $Z_{t}$ defined in Eq. (10) can be identified, up to equivalence in law, with the $Z_{t}$ generated by the latent variable model (13) with $A_{t}=(1-t)I$, $B_{t}=tI$, $L_{t}=\sqrt{\sigma_{min}^{2}+\sigma^{2}t(1-t)}I$. For $t\in[0,1)$, it can also be identified with the solution $\hat{Z}_{t}$ of the linear SDE (15) with $F_{t}=-I/(1-t)$, $H_{t}=(1+\frac{t}{1-t})I$ and $G_{t}=\sqrt{\sigma^{2}+\frac{2\sigma_{min}^{2}}{1-t}}I$. Moreover, $\lim_{t\to 1}\hat{Z}_{t}=^{d}Z_{1}+\sigma_{min}\epsilon$, where $\epsilon\sim\mathcal{N}(0,I)$ and $=^{d}$ denotes equivalence in distribution.

###### Proof.

The identification follows from matching the moments of $Z_{t}$ and $\hat{Z}_{t}$, i.e., applying Eq. (16)-(18).

To prove that $\lim_{t\to 1}\hat{Z}_{t}=^{d}Z_{1}+\sigma_{min}\epsilon$,
we use the explicit solution of the SDE:

$$ $\hat{Z}_{t}=\Phi(t,0)Z_{0}+\int_{0}^{t}\Phi(t,s)H_{s}Z_{1}ds+\int_{0}^{t}\Phi( t,s)G_{s}dW_{s},$ $$

where $\Phi(t,s)$ is the fundamental solution of the homogeneous equation $d\Phi(t,s)=F_{t}\Phi(t,s)dt$ with $\Phi(s,s)=I$.
For our $F_{t}=-I/(1-t)$, we have $\Phi(t,s)=\exp(-\int_{s}^{t}\frac{1}{1-u}du)I=(\frac{1-t}{1-s})I$.
Substituting this and the formula for $H_{t}$ into the solution, we obtain
$\hat{Z}_{t}=(1-t)Z_{0}+\int_{0}^{t}(\frac{1-t}{1-s})(1+\frac{s}{1-s})Z_{1}ds+
\int_{0}^{t}(\frac{1-t}{1-s})\sqrt{\sigma^{2}+\frac{2\sigma_{min}^{2}}{1-s}}
IdW_{s}$.

Now, let us examine each term as $t\to 1$.
First,
$(1-t)Z_{0}\to 0$ as $t\to 1$ and
$\int_{0}^{t}(\frac{1-t}{1-s})(1+\frac{s}{1-s})Z_{1}ds=tZ_{1}\to Z_{1}$ as $t\to 1$.
It remains to deal with the stochastic integral term
$M_{t}:=\int_{0}^{t}(\frac{1-t}{1-s})\sqrt{\sigma^{2}+\frac{2\sigma_{min}^{2}}{
1-s}}dW_{s}$. Note that $M_{t}$ is an Itô integral that has zero mean, i.e. $\mathbb{E}M_{t}=0$, and using Itô’s formula [],

$$ $\displaystyle\mathbb{E}M_{t}^{2}$ $\displaystyle=\int_{0}^{t}((1-t)/(1-s))^{2}\left(\sigma^{2}+\frac{2\sigma_{min }^{2}}{1-s}\right)ds$ (19) $\displaystyle=(1-t)^{2}\left[\sigma^{2}\int_{0}^{t}\frac{1}{(1-s)^{2}}ds+2 \sigma_{min}^{2}\int_{0}^{t}\frac{1}{(1-s)^{3}}ds\right]$ (20) $\displaystyle=(1-t)^{2}\left[\sigma^{2}\left(\frac{t}{1-t}\right)+\sigma_{min} ^{2}\left(\frac{1}{(1-t)^{2}}-1\right)\right]$ (21) $\displaystyle=(1-t)t\sigma^{2}+\sigma_{min}^{2}t(2-t),$ (22) $$

which tends to $\sigma_{min}^{2}$ as $t\to 1$. Combining the above results, $\lim_{t\to 1}\hat{Z}_{t}\sim\mathcal{N}(Z_{1},\sigma_{min}^{2}I)$.

∎

###### Remark 1 .

The linear SDE in Proposition 1 is not the unique SDE that Eq. (10) can be identified with at each $t$. For instance, $Z_{t}$ in Eq. (10) can also be identified with the solution $\hat{Z}_{t}$ of the SDE:

$$ $d\hat{Z}_{t}=(Z_{1}-Z_{0})dt+\frac{\sigma^{2}}{2}\frac{1-2t}{\sqrt{\sigma_{min }^{2}+\sigma^{2}t(1-t)}}dW_{t},\ \ t\in[0,1].$ (23) $$

It is straightfoward to check that $\hat{Z}_{t}$ and $Z_{t}$ are both Gaussian with the same mean and variance for all $t\in[0,1]$. Note that the probability path defined in Eq. (10) does not specify the covariance. However, since during training and inference, $Z_{t}$ and $Z_{s}$ are sampled using independent Gaussian samples, $Cov(Z_{t},Z_{s})=0$ for $t\neq s$. On the other hand, the SDEs that we have constructed above specify the covariance. Thus, while $Z_{t}=^{d}\hat{Z}_{t}$ for all $t\in[0,1]$, the paths $(Z_{t})_{t\in[0,1]}$ and $(\hat{Z}_{t})_{t\in[0,1]}$ may not be equivalent in law as they may have different covariance.
It is natural to ask how the choice of covariance structure affects the model performance. Note that generating the SDE path $(\hat{Z}_{t})_{t\in[0,1]}$ requires simulation of the SDE. This beats the original purpose of flow matching, which champions a simulation-free approach that directly works with a pre-specified probability path.

###### Proposition 1 .

###### Proof.

###### Remark 1 .

## Appendix B On Different Loss Parametrizations

In this section, we list popular choices of loss parametrization considered in the literature and connect them to our flow matching loss. We refer to [] for a more comprehensive discussion. Recall that the Gaussian path that we consider is: $Z_{t}=a_{t}Z_{0}+b_{t}Z_{1}+c_{t}\xi$, where $\xi\sim\mathcal{N}(0,I)$. In general, these loss parametrizations take the form of:

$$ $\mathcal{L}(\theta):=\mathbb{E}_{t,p_{t}(Z|\tilde{Z}),q(\tilde{Z})}\omega(t)\| m^{\theta}_{t}(Z)-m_{t}(Z,\tilde{Z})\|^{2},$ (24) $$

where $\omega(t)>0$ is a weighting function, $m_{t}(Z,\tilde{Z})$ is the object (conditioned on $\tilde{Z}$) to be learnt and $m_{t}^{\theta}$ is a neural network model used to learn the object of interest. Depending on which object one would like to learn/match, we have different loss parametrizations.

Flow matching. The flow matching loss that we focus in this paper is:

$$ $\mathcal{L}_{cfm}(\theta):=\mathbb{E}_{t,p_{t}(Z|\tilde{Z}),q(\tilde{Z})} \omega(t)\|v^{\theta}_{t}(Z)-u_{t}(Z|\tilde{Z})\|^{2},$ (25) $$

where one aims to learn the flow generating vector field:

$$ $u_{t}(Z|\tilde{Z})=\frac{c_{t}^{\prime}}{c_{t}}(Z-(a_{t}Z_{0}+b_{t}Z_{1}))+a_{ t}^{\prime}Z_{0}+b_{t}^{\prime}Z_{1},$ (26) $$

Score matching. The score matching loss is:

$$ $\mathcal{L}_{sm}(\theta):=\mathbb{E}_{t,p_{t}(Z|\tilde{Z}),q(\tilde{Z})} \lambda(t)\|s^{\theta}_{t}(Z)-\nabla\log p_{t}(Z|\tilde{Z})\|^{2},$ (27) $$

where $\lambda(t)>0$ is a weighting function and one aims to learn the score function:

$$ $\nabla\log p_{t}(Z|\tilde{Z})=\frac{a_{t}Z_{0}+b_{t}Z_{1}-Z}{c_{t}^{2}}.$ (28) $$

If $\lambda(t)=c_{t}^{2}$, then this reduces to the original score matching loss [], whereas if $\lambda(t)=\beta(1-t)$, this becomes the score flow loss [].

Noise matching.
The noise matching loss is:

$$ $\mathcal{L}_{nm}(\theta):=\mathbb{E}_{t,p_{t}(Z|\tilde{Z}),q(\tilde{Z})}\| \epsilon^{\theta}_{t}(Z)-\epsilon_{t}(Z|\tilde{Z})\|^{2},$ (29) $$

where one aims to learn the noise:

$$ $\epsilon_{t}(Z|\tilde{Z})=\frac{Z-(a_{t}Z_{0}+b_{t}Z_{1})}{c_{t}}.$ (30) $$

## Appendix C Theoretical Results and Proofs

In this section, we provide theorems and proofs to justify the discussions in Section 3 and Section 4.

### C.1 Connecting Flow Matching with Conditional Flow Matching

The following theorem justifies the claim that minimizing $\mathcal{L}_{fm}$ is equivalent to minimizing $\mathcal{L}_{cfm}$.

###### Theorem 1 .

If the conditional vector field $u_{t}(Z|\tilde{Z})$ generates the conditional probability path $p_{t}(Z|\tilde{Z})$, then the marginal vector field $u_{t}$ in Eq. (6) generates the marginal probability path $p_{t}$ in Eq. (5). Moreover, if $p_{t}(Z)>0$ for all $t,Z$, then $\mathcal{L}_{fm}$ and $\mathcal{L}_{cfm}$ are equal up to a constant independent of $\theta$.

###### Proof.

The proof is a straightforward extension of the proofs of Theorem 1-2 in [] from conditioning on data samples to conditioning on latent samples and allowing an arbitrary weighting function $\omega(t)$.

Suppose that the conditional vector field $u_{t}(Z|\tilde{Z})$ generates the conditional probability path $p_{t}(Z|\tilde{Z})$, we would like to show that the marginal vector field $u_{t}$ in Eq. (6) generates the marginal probability path $p_{t}$ in Eq. (5). To show this, it suffices to verify that $p_{t}$ and $u_{t}$ satisfy the continuity equation:

$$ $\frac{d}{dt}p_{t}(Z)+div(p_{t}(Z)u_{t}(Z))=0,$ (31) $$

where the divergence operator, $div$, is defined with respect to the latent variable $Z=(Z^{1},\dots,Z^{d})$, i.e., $div=\sum_{i=1}^{d}\frac{\partial}{\partial Z^{i}}$.

We begin by taking derivative of $p_{t}(Z)$ with respect to time:

$$ $\displaystyle\frac{d}{dt}p_{t}(Z)$ $\displaystyle=\frac{d}{dt}\int p_{t}(Z|\tilde{Z})q(\tilde{Z})d\tilde{Z}$ (32) $\displaystyle=\int\frac{d}{dt}p_{t}(Z|\tilde{Z})q(\tilde{Z})d\tilde{Z}$ (33) $\displaystyle=-\int div(u_{t}(Z|\tilde{Z})p_{t}(Z|\tilde{Z}))q(\tilde{Z})d \tilde{Z}$ (34) $\displaystyle=-div\left(\int u_{t}(Z|\tilde{Z})p_{t}(Z|\tilde{Z})q(\tilde{Z})d \tilde{Z}\right)$ (35) $\displaystyle=-div(u_{t}(Z)p_{t}(Z)).$ (36) $$

In the third line, we use the fact that $u_{t}(\cdot|\tilde{Z})$ generates $p_{t}(\cdot|\tilde{Z})$. In the last line, we use Eq. (6). In the second and forth line above, the exchange of integration and differentiation can be justified by assuming that the integrands satisfy the regularity conditions of the Leibniz rule.

Next, we would like to show that if $p_{t}(Z)>0$ for all $t,Z$, then $\mathcal{L}_{fm}$ and $\mathcal{L}_{cfm}$ are equal up to a constant independent of $\theta$. We follow [] and assume that $q(Z)$ and $p_{t}(Z|\tilde{Z})$ are decreasing to zero sufficiently fast as $\|Z\|\to 0$, that $u_{t},v_{t},\nabla_{\theta}v_{t}$ are bounded, so that all integrals exist and exchange of integration order is justified via Fubini’s theorem.

Using the bilinearity of the 2-norm, we decompose the squared losses as:

$$ $\displaystyle\|v_{t}^{\theta}(Z)-u_{t}(Z)\|^{2}$ $\displaystyle=\|v_{t}^{\theta}(Z)\|^{2}-2\langle v^{\theta}_{t}(Z),u_{t}(Z) \rangle+\|u_{t}(Z)\|^{2},$ (37) $\displaystyle\|v^{\theta}_{t}(Z)-u_{t}(Z|\tilde{Z})\|^{2}$ $\displaystyle=\|v^{\theta}_{t}(Z)\|^{2}-2\langle v^{\theta}_{t}(Z),u_{t}(Z| \tilde{Z})\rangle+\|u_{t}(Z|\tilde{Z})\|^{2}.$ (38) $$

Now,

$$ $\displaystyle\mathbb{E}_{p_{t}(Z)}\|v^{\theta}_{t}(Z)\|^{2}$ $\displaystyle=\int\|v^{\theta}_{t}(Z)\|^{2}p_{t}(Z)dZ$ (39) $\displaystyle=\int\int\|v^{\theta}_{t}(Z)\|^{2}p_{t}(Z|\tilde{Z})q(\tilde{Z})d \tilde{Z}dZ$ (40) $\displaystyle=\mathbb{E}_{q(\tilde{Z}),p_{t}(Z|\tilde{Z})}\|v^{\theta}_{t}(Z) \|^{2},$ (41) $$

where we use Eq. (5) in the second equality above and exchange the order of integration in the third equality.

Next, we compute:

$$ $\displaystyle\mathbb{E}_{p_{t}(Z)}\langle v^{\theta}_{t}(Z),u_{t}(Z)\rangle$ $\displaystyle=\int\bigg{\langle}v^{\theta}_{t}(Z),\frac{\int u_{t}(Z|\tilde{Z} )p_{t}(Z|\tilde{Z})q(\tilde{Z})d\tilde{Z}}{p_{t}(Z)}\bigg{\rangle}p_{t}(Z)dZ$ (42) $\displaystyle=\int\bigg{\langle}v^{\theta}_{t}(Z),\int u_{t}(Z|\tilde{Z})p_{t} (Z|\tilde{Z})q(\tilde{Z})d\tilde{Z}\bigg{\rangle}dZ$ (43) $\displaystyle=\int\int\langle v^{\theta}_{t}(Z),u_{t}(Z|\tilde{Z})\rangle p_{t }(Z|\tilde{Z})q(\tilde{Z})d\tilde{Z}dZ$ (44) $\displaystyle=\mathbb{E}_{q(\tilde{Z}),p_{t}(Z|\tilde{Z})}\langle v^{\theta}_{ t}(Z),u_{t}(Z|\tilde{Z}),$ (45) $$

where we first plug in Eq. (6) and then exchange the order the integration in order to arrive at the last equality.

Finally, noting that $u_{t}$ are $\omega(t)$ independent of $\theta$ (and are thus irrelevant for computing the loss gradients), we have proved the desired result.
∎

###### Theorem 1 .

###### Proof.

### C.2 Identifying the Vector Field that Generates the Gaussian Paths

Similar to Theorem 3 in [], we have the following result, which identifies the unique vector field that generates the Gaussian probability path.

###### Theorem 2 .

Let $p_{t}(Z|\tilde{Z})$ be the Gaussian probability path defined in Eq. (8) and consider the flow map $\psi_{t}$ defined as $\psi_{t}(Z)=a_{t}Z_{0}+b_{t}Z_{1}+c_{t}Z$ with $c_{t}>0$. Then the unique vector field that defines $\psi_{t}$ is:

$$ $u_{t}(Z|\tilde{Z})=\frac{c_{t}^{\prime}}{c_{t}}(Z-(a_{t}Z_{0}+b_{t}Z_{1}))+a_{ t}^{\prime}Z_{0}+b_{t}^{\prime}Z_{1},$ (46) $$

where prime denotes derivative with respect to $t$, and $u_{t}(Z|\tilde{Z})$ generates the Gaussian path $p_{t}(Z|\tilde{Z})$.

###### Proof.

Let $p_{t}(Z|\tilde{Z})$ be the Gaussian probability path defined in Eq. (8) and consider the flow map $\psi_{t}$ defined as $\psi_{t}(Z)=a_{t}Z_{0}+b_{t}Z_{1}+c_{t}Z$. We would like to show that the unique vector field that defines $\psi_{t}$ is:

$$ $u_{t}(Z|\tilde{Z})=\frac{c_{t}^{\prime}}{c_{t}}(Z-(a_{t}Z_{0}+b_{t}Z_{1}))+a_{ t}^{\prime}Z_{0}+b_{t}^{\prime}Z_{1},$ (47) $$

and $u_{t}(Z|\tilde{Z})$ generates the Gaussian path $p_{t}(Z|\tilde{Z})$.

We denote $w_{t}=u_{t}(Z|\tilde{Z})$ for notational simplicity. Then,

$$ $\frac{d}{dt}\psi_{t}(Z)=w_{t}(\psi_{t}(Z)).$ (48) $$

Since $\psi_{t}$ is invertible (as $c_{t}>0$), we let $Z=\psi^{-1}(Y)$ and obtain

$$ $\psi_{t}^{\prime}(\psi^{-1}(Y))=w_{t}(Y),$ (49) $$

where the prime denotes derivative with respect to $t$ and we have used the apostrosphe notation for the derivative to indicate that $\psi_{t}^{\prime}$ is evaluated at $\psi^{-1}(Y)$.

Inverting $\psi_{t}(Z)$ gives:

$$ $\psi_{t}^{-1}(Y)=\frac{Y-\mu_{t}(\tilde{Z})}{c_{t}},$ (50) $$

where $\mu_{t}(\tilde{Z}):=a_{t}Z_{0}+b_{t}Z_{1}$.

Differentiating $\psi_{t}$ with respect to $t$ gives $\psi_{t}^{\prime}(Z)=c_{t}^{\prime}Z+\mu_{t}^{\prime}(\tilde{Z})$.

Plugging the last two equations into Eq. (49), we obtain:

$$ $w_{t}(Y)=\frac{c_{t}^{\prime}}{c_{t}}(Y-\mu_{t}(\tilde{Z}))+\mu_{t}^{\prime}( \tilde{Z})$ (51) $$

which is the result that we wanted to show.

∎

###### Theorem 2 .

###### Proof.

### C.3 Comparing the Variance of Different Vector Field Models

We show that under reasonable assumptions, the variance of the vector field corresponding to our proposed probability path model is lower than the variance of the vector field corresponding to the other choices of probability paths. Here the variance is taken with respect to the randomness in the samples $z^{\tau}$ and the Gaussian samples drawn during gradient descent updates. To simplify our analysis and to facilitate discussion, we only compare our probability path model to the model generated by the optimal transport VF of [], setting $\epsilon_{min}:=0$ (rectified flow) without loss of generality.

We start with the following lemma.

###### Lemma 1 .

Let $A,B,C,D$ be random vectors where $C$ and $D$ are independent, both $A$ and $B$ are independent of $C$ and $D$ (but $A$ and $B$ could be dependent). If $Cov(A,B)\geq(Var(C)-Var(D)+Var(B))/2$, then

$$ $Var(A+D)\geq Var(A-B+C)=Var(B-A+C),$ (52) $$

where $A\geq B$ means that $A-B$ is positive semidefinite.

###### Proof.

We compute:

$$ $\displaystyle Var(A+D)$ $\displaystyle=Var(A-B+C+B+D-C)$ (53) $\displaystyle=Var(A-B+C)+Var(B+D)+Var(C)+2Cov(A-B+C,B+D)$ $\displaystyle\ \ \ \ -2Cov(A-B+C,C)-2Cov(B+D,C)$ (54) $\displaystyle=Var(A-B+C)+Var(B+D)-Var(C)+2Cov(A-B+C,B)$ $\displaystyle\ \ \ \ +2Cov(A-B+C,D)$ (55) $\displaystyle=Var(A-B+C)+Var(B+D)-Var(C)+2Cov(A,B)-2Var(B)$ (56) $\displaystyle=Var(A-B+C)-Var(C)+Var(D)+2Cov(A,B)-Var(B),$ (57) $$

where we have simply rearranged the terms in the first equality, used the formula $Var(A+B+C)=Var(A)+Var(B)+Var(C)+2Cov(A,B)+2Cov(A,C)+2Cov(B,C)$, bilinearity of covariance, the facts that $Cov(A,A)=Var(A)$ and $Var(cA)=c^{2}A$ for a scalar $c$, as well as the assumption that both $A$ and $B$ are independent of $C$, $D$ in the last four equalities.

Therefore, if $-Var(C)+Var(D)+2Cov(A,B)-Var(B)\geq 0$, then we have $Var(A+D)\geq Var(A-B+C)$.

∎

Now, we consider the vector fields that correspond to the above two models of probability path, denoting the vector field that corresponds to our model as $u_{t}$ and the vector field that corresponds to the rectified flow model as $\tilde{u}_{t}$. They generate the probability path $Z_{t}=tz^{\tau}+(1-t)z^{\tau-1}+c_{t}\xi$ and $\tilde{Z}_{t}=tz^{\tau-1}+(1-t)\eta$ respectively, where $\xi,\eta\sim\mathcal{N}(0,I)$, $c_{t}=\sqrt{\sigma_{min}^{2}+\sigma^{2}t(1-t)}$ and $t\in[0,1]$.

Applying Eq. (9), we have:

$$ $\displaystyle u_{t}(Z_{t}|z^{\tau-1},z^{\tau})$ $\displaystyle=z^{\tau}-z^{\tau-1}+c_{t}^{\prime}\xi,$ (58) $\displaystyle\tilde{u}_{t}(\tilde{Z}_{t}|z^{\tau-1})$ $\displaystyle=z^{\tau-1}-\eta.$ (59) $$

###### Theorem 3 .

Suppose that $(z^{\tau})_{\tau=1,\dots,m}$, with the $z^{\tau}\in\mathbb{R}^{d}$, is a discrete-time stochastic process with nonzero correlation in time and let $t\in[0,1]$ be given.

If $Cov(z^{\tau-1},z^{\tau})\geq\frac{1}{2}\left(\left(\frac{\sigma^{4}}{4\sigma_{
min}^{2}}-1\right)I+Var(z^{\tau})\right)$, then
$Var(\tilde{u}_{t}(\tilde{Z}_{t}|z^{\tau-1}))\geq Var(u_{t}(Z_{t}|z^{\tau-1},z^
{\tau})).$

###### Proof.

Note that $c_{t}^{\prime}=\frac{\sigma^{2}(1-2t)}{2\sqrt{\sigma_{min}^{2}+\sigma^{2}t(1-t
)}}$, $Var(c_{t}^{\prime}\xi)=(c_{t}^{\prime})^{2}I$ and $Var(-\eta)=I$. Therefore, using these and applying Lemma 1 with $A:=z^{\tau-1}$, $B:=z^{\tau}$, $C:=c_{t}^{\prime}\xi$ and $D:=-\eta$, allow us to establish the claim that $Var(\tilde{u}_{t}(Z_{t}|z^{\tau-1}))\geq Var(u_{t}(Z_{t}|z^{\tau-1},z^{\tau}))$ if

$$ $Cov(z^{\tau-1},z^{\tau})\geq\frac{1}{2}\left(\frac{\sigma^{4}(1-2t)^{2}I}{4( \sigma_{min}^{2}+\sigma^{2}t(1-t))}+Var(z^{\tau})-I\right).$ (60) $$

Since the function $f(t):=\frac{\sigma^{4}(1-2t)^{2}}{4(\sigma_{min}^{2}+\sigma^{2}t(1-t))}$ is maximized at the endpoints $t=0,1$ with the maximum value of $\sigma^{4}/4\sigma_{min}^{2}$, the desired result stated in the theorem follows.
∎

Theorem 3 implies that if the consecutive latent variables $z^{\tau}$, $z^{\tau-1}$ are sufficiently correlated and $\sigma$ is chosen small enough relative to $\sigma_{min}$, then the variance of the vector field that corresponds to our probability path model is lower than that corresponds to the rectified flow model.

Discussions.
Let us denote $v_{t}^{\theta}(Z):=v^{\theta}_{t}(z^{\tau},z^{\tau-1},\xi)$ to show the explicit dependence of the vector field neural net on the random samples $t,z^{\tau}$, $z^{\tau-1}$ and $\xi\sim\mathcal{N}(0,I)$ drawn during each update of gradient descent during training. During each gradient descent update, our model involves computation of

$$ $\displaystyle\nabla_{\theta}\mathcal{L}_{cfm}(\theta;t,\xi,z^{\tau},z^{\tau-1})$ $\displaystyle=2\nabla_{\theta}v_{t}^{\theta}(z^{\tau},z^{\tau-1},\xi)^{T}\cdot (v_{t}^{\theta}(z^{\tau},z^{\tau-1},\xi)-u_{t}(Z_{t}|z^{\tau},z^{\tau-1}))$ (61) $\displaystyle=:2\nabla_{\theta}v_{t}^{\theta}(z^{\tau},z^{\tau-1},\xi)^{T} \cdot\Delta_{t}^{\theta}(z^{\tau-1},z^{\tau},\xi),$ (62) $$

with $t\sim\mathcal{U}[0,1]$, $\xi\sim\mathcal{N}(0,I)$ and the latent samples $z^{\tau},z^{\tau-1}$ drawn randomly.

Similarly, for the rectified flow model let us denote $\tilde{v}_{t}^{\theta}(Z):=v^{\theta}_{t}(z^{\tau-1},\eta)$ to show the explicit dependence of the vector field neural net on the random samples $t$, $z^{\tau-1}$ and $\eta\sim\mathcal{N}(0,I)$ drawn during each update of gradient descent during training.
Each update of gradient descent using the rectified flow model involves computation of

$$ $\displaystyle\nabla_{\theta}\tilde{\mathcal{L}}_{cfm}(\theta;t,\eta,z^{\tau-1})$ $\displaystyle=2\nabla_{\theta}\tilde{v}_{t}^{\theta}(z^{\tau-1},\eta)^{T}\cdot (\tilde{v}_{t}^{\theta}(z^{\tau-1},\eta)-\tilde{u}_{t}(\tilde{Z}_{t}|z^{\tau-1 }))$ (63) $\displaystyle=:2\nabla_{\theta}\tilde{v}_{t}^{\theta}(z^{\tau-1},\eta)^{T} \cdot\tilde{\Delta}_{t}^{\theta}(z^{\tau-1},\eta),$ (64) $$

with $t\sim\mathcal{U}[0,1]$, $\eta\sim\mathcal{N}(0,I)$ and the latent sample $z^{\tau-1}$ drawn randomly.

Lower gradient variance results in smoother training loss curve and potentially faster convergence, so it is useful to compare the variances of the loss gradient for the two models.
However, the variances are highly dependent on $\nabla_{\theta}v_{t}^{\theta}$, $\nabla_{\theta}\tilde{v}_{t}^{\theta}$ and their covariance with the other random vectors appearing in Eq. (61) and Eq. (63), making such comparison challenging without strong assumptions. Heuristically, the difference in the variances of the loss gradient during each update for the two models is primarily determined by the difference between $Var(\Delta_{t}^{\theta}(z^{\tau-1},z^{\tau},\xi))$ and $Var(\tilde{\Delta}_{t}^{\theta}(z^{\tau-1},\eta))$ if $\nabla_{\theta}v_{t}^{\theta}$ and $\nabla_{\theta}\tilde{v}_{t}^{\theta}$ are relatively stable. In this case, we have $Var(\Delta_{t}^{\theta}(z^{\tau-1},z^{\tau},\xi))\leq Var(\tilde{\Delta}_{t}^{
\theta}(z^{\tau-1},\eta))$ if we suppose the assumptions in Theorem 3, $Var(\tilde{v}_{t}^{\theta})\geq Var(v_{t}^{\theta})$ and $Cov(v_{t}^{\theta},u_{t})\geq Cov(\tilde{v}_{t}^{\theta},\tilde{u}_{t})$.

The implications of Theorem 3 together with the heuristics above could partially explain why our probability path model leads to smoother loss curve and faster convergence (see Figure 2) compared to the RIVER method of [].
On the other hand, the dependence of the lower bound in the theorem on $\sigma$ and $\sigma_{min}$ suggests that using values of $\sigma$ that is relatively large enough might not keep the variance of the vector field low, which could partially explain the phenomenon displayed in Figure 2, where using $\sigma=0.1$ and $\sigma_{min}=0.001$ leads to large loss fluctuations.

###### Lemma 1 .

###### Proof.

###### Theorem 3 .

###### Proof.

## Appendix D Additional Empirical Results

In this section, we provide additional experimental results.

Figure 4 provides visual results of the predicted snapshots by our model for each task. The visualizations highlight our model’s ability to capture complex flow patterns and dynamics.

Figure: Figure 4: Visualization of predicted frames using our model of probability path for the considered tasks. From top to bottom: fluid flow past a cylinder, shallow-water equation, diffusion-reaction equation, and incompressible Navier-Stokes equation. In each case, GT indicates the ground truth frames and pred. indicates the predicted frames.

## Appendix E Experimental Details

In this section, we provide the experimental details for the tasks considered in Section 6.

### E.1 Details on the Datasets

Fluid flow past a cylinder.
We use the fluid flow past a stationary cylinder at a Reynolds number of 100 as a simple test problem. This fluid flow is a canonical problem in fluid dynamics characterized by a periodically shedding wake structure [, ]. The flow dynamics are governed by the two-dimensional incompressible Navier–Stokes equations:

$$ $\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}=-\frac {1}{\rho}\nabla p+\nu\nabla^{2}\mathbf{u},$ $$

$$ $\nabla\cdot\mathbf{u}=0,$ $$

where $\mathbf{u}=(u,v)$ is the velocity field, $p$ is the pressure, $\rho$ is the fluid density, and $\nu$ is the kinematic viscosity. The vorticity field $\omega$ is obtained from the velocity field via:

$$ $\omega=\nabla\times\mathbf{u},$ $$

providing insights into the rotational characteristics of the flow.

For simulating the data, the Immersed Boundary Projection Method (IBPM) has been used []. The flow tensor has dimensions $199\times 449\times 151$, representing 151 temporal snapshots on a $449\times 199$ spatial grid. We crop and spatially subsample the data which results in a $64\times 64$ spatial field.

Shallow-water equation.
The shallow-water equations, derived from the compressible Navier-Stokes equations, can be used for modeling free-surface flow problems. We consider the 2D equation used in [], which is the following system of hyperbolic PDEs:

$$ $\partial_{t}h+\nabla h{\bf u}=0,\ \ \partial_{t}h{\bf u}+\nabla\left({\bf u}^{ 2}h+\frac{1}{2}g_{r}h^{2}\right)=-g_{r}h\nabla b,$ (65) $$

where ${\bf u}=u,v$ being the velocities in the horizontal and vertical direction respectively, $h$ describes the water depth, and $b$ describes a spatially varying bathymetry. $h{\bf u}$ can be interpreted as the directional
momentum components and $g_{r}$ describes the gravitational acceleration. The mass and momentum conservation properties can hold across shocks in the solution and thus challenging datasets can be
generated. This equation finds application in modeling tsunamis and flooding events.

We use the dataset generated and provided by PDEBench []. The data file ([2D_rdb_NA_NA.h5](2D_rdb_NA_NA.h5)) can be downloaded from [https://github.com/pdebench/PDEBench/tree/main/pdebench/data_download](https://github.com/pdebench/PDEBench/tree/main/pdebench/data_download). The data sample is a series of 101 frames at a 128 $\times$ 128 pixel resolution and come with 1 channel.
The simulation considered in [] is a 2D radial dam break scenario. On a square domain $\Omega=[-2.5,2.5]^{2}$, the
water height is initialized as a circular bump in the center of the domain:

$$ $h(t=0,x,y)=\begin{cases}2,\text{ for }r<\sqrt{x^{2}+y^{2}},\\ 1,\text{ for }r\geq\sqrt{x^{2}+y^{2}},\end{cases}$ $$

with the radius $r$ randomly sampled from $\mathcal{U}(0.3,0.7)$. The dataset is simulated with a finite volume solver using the PyClaw package.
We apply standardization and then normalization to the range of $[-1,1]$ to preprocess the simulated data.

Incompressible Navier-Stokes equation.
The Navier-Stokes equation is the incompressible version of the compressible Navier-Stokes equation, and it can be used to
model hydromechanical systems,
turbulent dynamics and weather. We use the inhomogeenous version of the equation (which includes a vector field forcing term ${\bf u}$) considered by []:

$$ $\nabla\cdot{\bf v}=0,\ \ \ \rho(\partial_{t}{\bf v}+{\bf v}\cdot\nabla{\bf v}) =-\nabla p+\eta\Delta{\bf v}+{\bf u},$ (66) $$

where $\rho$ is the mass density, ${\bf v}$ is the fluid velocity, $p$ is the gas pressure and $\eta$ is shear viscosity.
The initial conditions ${\bf v}_{0}$ and inhomogeneous forcing parameters ${\bf u}$ are each drawn from isotropic Gaussian random fields with truncated power-law decay $\tau$ of the power spectral density and scale $\sigma$, where
$\tau_{v_{0}}=-3$, $\sigma_{v_{0}}=0.15$, $\tau_{u}=-1$, $\sigma_{u}=0.4$. The domain is taken to be the unit square $\Omega=[0,1]^{2}$ and the viscosity $\eta=0.01$.
The equation is numerically simulated using Phiflow. Boundary conditions are taken to be Dirichlet to clamp the field
velocity to zero at the perimeter.

We use the dataset generated and provided by PDEBench []. The data file ([ns_incom_inhom_2d_512-0.h5](ns_incom_inhom_2d_512-0.h5)) can be downloaded from
[https://github.com/pdebench/PDEBench/tree/main/pdebench/data_download](https://github.com/pdebench/PDEBench/tree/main/pdebench/data_download). The data sample is a series of 1000 frames at a 512 $\times$ 512 pixel resolution and come with 2 channels.
We do not apply any data preprocessing procedure here.

Diffusion-reaction equation.
We use the 2D extension of diffusion-reaction equation of [] which describes two non-linearly coupled variables, namely
the activator $u=u(t,x,y)$ and the inhibitor $v=v(t,x,y)$. The equation is given by:

$$ $\displaystyle\partial_{t}u$ $\displaystyle=D_{u}\partial_{xx}u+Du\partial_{yy}u+R_{u},$ (67) $\displaystyle\partial_{t}v$ $\displaystyle=D_{v}\partial_{xx}v+D_{v}\partial_{yy}v+R_{v},$ (68) $$

where $D_{u}$ and $D_{v}$ are the diffusion coefficient for the activator and inhibitor respectively, $R_{u}=R_{u}(u,v)$ and $R_{v}=R_{v}(u,v)$ are the activator and inhibitor reaction function respectively. The
domain of the simulation includes $x\in(-1,1)$, $y\in(-1,1)$, $t\in(0,5]$. This equation can be used for modeling biological pattern formation.

The reaction functions for the activator and inhibitor are defined by the Fitzhugh-Nagumo equation as:
$R_{u}(u,v)=u-u^{3}-k-v$, $R_{v}(u,v)=u-v$,
where $k=5\times 10^{-3}$, and the diffusion coefficients for the activator and inhibitor are $D_{u}=1\times 10^{-3}$ and $D_{v}=5\times 10^{-3}$ respectively. The initial condition is generated as standard Gaussian
noise $u(0,x,y)\sim\mathcal{N}(0,1.0)$ for $x\in(-1,1)$ and $y\in(-1,1)$. We take a no-flow Neumann boundary condition: $D_{u}\partial_{x}u=0$, $D_{v}\partial_{x}v=0$, $D_{u}\partial_{y}u=0$, and $D_{v}\partial_{y}v=0$ for $x,y\in(-1,1)^{2}$.

We use a downsampled version of the dataset provided by PDEBench []. The data file ([2D_diff-react_NA_NA.h5](2D_diff-react_NA_NA.h5)) can be downloaded from [https://github.com/pdebench/PDEBench/tree/main/pdebench/data_download](https://github.com/pdebench/PDEBench/tree/main/pdebench/data_download). The data sample is a series of 101 frames at a 128 $\times$ 128 pixel resolution and come with 2 channels. The sample frames are generated using the finite volume method for spatial discretization, and the time integration is performed using the built-in fourth order Runge-Kutta method in the scipy package.
We do not apply any data preprocessing procedure here.

### E.2 Details on Pre-Training the Autoencoder

We use the same architecture for the encoder and decoder for all the tasks, with the architecture parameters chosen based on the complexity of the task.

The encoder. The encoder first applies a 2D convolution (conv_in) to the input frame, which reduces the number of channels from in_channels to mid_channels, and processes the spatial dimensions. Then, a series of ResidualBlock layers, which progressively process and downsample the feature map, making it smaller in spatial dimensions but more enriched in terms of features, are applied. After the residual blocks, the feature map undergoes an attention process via a multi-head attention layer. This layer helps the encoder focus on important parts of the input, learning relationships between spatial positions in the image. For the post-attention step, the feature map is further processed by residual blocks and normalized, preparing it for the final convolution. The output of the encoder is obtained by applying a final 2D convolution (out_conv), which maps the processed feature map to the desired number of output channels (out_channels).

The decoder. The decoder takes the encoded feature map and transforms it back into an output with similar spatial dimensions as the input. Similar to the encoder, the decoder starts with a convolution that adjusts the number of channels from in_channels to mid_channels. Then, an attention mechanism (similar to the encoder) is applied to focus on important aspects of the encoded features. Next, a series of ResidualBlock layers, combined with UpBlock layers, are used to progressively increase the spatial dimensions of the feature map (upsampling), undoing the compression applied by the encoder. After the upsampling, the output is normalized and passed through a final convolution (out_conv), mapping the internal feature representation to the desired number of output channels (out_channels).

Table 4 summarizes the architecture parameters used for the considered tasks.

**Table 4: Parameters chosen for the encoder (decoder) architecture.**
| Task | Fluid flow | Shallow-water eq. | Navier-Stokes eq. | Diffusion-reaction eq. |
| --- | --- | --- | --- | --- |
| in_channels | 1 (1) | 1 (1) | 2 (2) | 2 (2) |
| out_channels | 1 (1) | 1 (1) | 2 (2) | 2 (2) |
| mid_channels | 64 (128) | 128 (256) | 128 (256) | 128 (256) |

Training details.
We train the autoencoder using AdamW with batch size of 32, no weight decay and $\beta=(0.9,0.999)$. We use the cosine learning rate scheduler with warmup. For the fluid flow past a cylinder task, we train for 2000 epochs and use learning rate of 0.001. For the Navier-Stokes task, we train for 500 epochs and use learning rate of 0.0001. For the other two tasks we train for 5000 epochs and use learning rate of 0.0005.
Our implementation is in PyTorch, and all experiments are run on an NVIDIA A100-SXM4 GPU
with 40 GB VRAM belonging to an internal SLURM cluster.

### E.3 Details on Training the Flow Matching Models

Architecture.
The vector field regressor is a transformer-based model designed to process latent vector fields and predict refined outputs with spatial and temporal dependencies. It uses key parameters like depth and mid_depth, which control the number of transformer encoder layers in the input, middle, and output stages. The state_size and state_res parameters define the number of channels and spatial resolution of the input data, while inner_dim sets the embedding dimension for processing. The model uses learned positional encodings and a time projection to incorporate spatial and temporal context into the input, which can include input_latents, reference_latents, and conditioning_latents. The input is projected into the inner dimension and passed through a series of transformer layers, with intermediate outputs from the input blocks concatenated with the output layers to refine predictions. Finally, the model projects the processed data back to the original spatial resolution and channel size using BatchNorm, producing the final vector field output.

Table 5 summarizes the architecture parameters used for the considered tasks.

**Table 5: Parameters chosen for the vector field neural network.**
| Parameter | Fluid flow | Shallow-water eq. | Navier-Stokes eq. | Diffusion-reaction eq. |
| --- | --- | --- | --- | --- |
| state_size | 4 | 4 | 8 | 4 |
| state_res | [8,8] | [16, 16] | [64, 64] | [16, 16] |
| inner_dim | 512 | 512 | 512 | 512 |
| depth | 4 | 4 | 4 | 4 |
| mid_depth | 5 | 5 | 5 | 5 |

Training details.
For all the considered tasks, we train the regressor using AdamW with batch size of 32, learning rate of 0.0005, no weight decay and $\beta=(0.9,0.999)$. We use the cosine learning rate scheduler with warmup. For the fluid flow past cylinder, we train for 2000 epochs, for the shallow-water equation and diffusion-reaction task we train for 1000 epochs, and for the Navier-Stokes task we train for 100 epochs. Our implementation is in PyTorch, and all experiments are run on an NVIDIA A100-SXM4 GPU
with 40 GB VRAM belonging to an internal SLURM cluster.

### E.4 Details on the Evaluation Metrics

In addition to the standard mean squared error and relative Frobenius norm error (RFNE), we use
the Pearson correlation coefficient to measure the linear relationship between the forecasted frames and the target frames. The range of this coefficient is $[-1,1]$, with zero implying no correlation. Correlations of $-1$ or $+1$ imply an exact linear relationship. Positive correlations imply that as $x$ increases, so does $y$. Negative correlations imply that as $x$ increases, $y$ decreases.
In addition, we use peak signal-to-noise ratio (PSNR) to evaluate the quality of signal representation against corrupting noise, and structural similarity index measure (SSIM) [] to assess perceptual results. The presented results are computed by averaging over batch size and number of sample generations.

### E.5 Standard Deviations for the Presented Results

Table 6-10 provide the standard deviation of the results presented in the main paper.

**Table 6: Standard deviation results for the fluid flow past a cylinder task using different choices of probability paths for flow matching. Results are averaged over 5 generations.**
| Model | Test MSE ($\downarrow$) | Test RFNE ($\downarrow$) | PSNR ($\uparrow$) | SSIM ($\uparrow$) |
| --- | --- | --- | --- | --- |
| RIVER | 1.33e-03 | 3.09e-02 | 1.31 | 5.93e-03 |
| VE-diffusion | 1.75e-02 | 3.60e-01 | 3.25e-01 | 2.31e-02 |
| VP-diffusion | 1.75e-02 | 3.60e-01 | 3.25e-01 | 2.31e-02 |
| Stochastic interpolant | 5.45e-05 | 3.51e-02 | 7.74e-02 | 3.75e-04 |
| Ours ($\sigma=0.01$, $\sigma_{sam}$ = 0, Euler) | 6.66e-06 | 4.49e-03 | 5.34e-02 | 3.48e-05 |

**Table 7: Standard deviation for the results of ablation study for the fluid flow past a cylinder task. Results are averaged over 5 generations.**
| $\sigma$ | sampler | sampling steps | Test MSE ($\downarrow$) | Test RFNE ($\downarrow$) | PSNR ($\uparrow$) | SSIM ($\uparrow$) |
| --- | --- | --- | --- | --- | --- | --- |
| 0.0 | Euler | 5 | 3.57e-05 | 2.42e-02 | 3.98e-02 | 2.51e-04 |
| 0.01 | Euler | 5 | 2.89e-05 | 2.45e-02 | 3.01e-02 | 1.80e-04 |
| 0.1 | Euler | 5 | 1.64e-05 | 2.73e-02 | 2.11e-02 | 7.21e-05 |
| 0.0 | Euler | 10 | 9.08e-06 | 1.25e-02 | 4.08e-02 | 5.73e-05 |
| 0.01 | Euler | 10 | 6.66e-06 | 4.49e-03 | 5.34e-02 | 3.48e-05 |
| 0.1 | Euler | 10 | 2.84e-05 | 4.00e-02 | 3.76e-02 | 2.21e-04 |
| 0.0 | RK4 | 10 | 3.36e-06 | 7.69e-03 | 3.08e-02 | 2.33e-05 |
| 0.01 | RK4 | 10 | 3.14e-06 | 5.43e-03 | 2.56e-02 | 2.22e-05 |
| 0.1 | RK4 | 10 | 5.88e-05 | 5.83e-02 | 2.03e-2 | 1.85e-04 |
| 0.0 | Euler | 20 | 6.35e-06 | 7.36e-03 | 3.90e-02 | 5.90e-05 |
| 0.01 | Euler | 20 | 9.64e-06 | 1.21e-02 | 4.30e-02 | 5.94e-05 |
| 0.1 | Euler | 20 | 5.44e-06 | 1.39e-02 | 1.39e-02 | 1.03e-04 |
| 0.0 | RK4 | 20 | 1.09e-05 | 1.65e-02 | 3.23e-02 | 8.08e-05 |
| 0.01 | RK4 | 20 | 8.35e-06 | 1.32e-02 | 3.29e-02 | 6.54e-05 |
| 0.1 | RK4 | 20 | 3.44e-06 | 1.07e-02 | 1.08e-02 | 4.08e-05 |

**Table 8: Standard deviation results for the shallow-water equation task using different choices of probability paths for flow matching. Results are averaged over 5 generations.**
| Model | Test MSE ($\downarrow$) | Test RFNE ($\downarrow$) | PSNR ($\uparrow$) | SSIM ($\uparrow$) |
| --- | --- | --- | --- | --- |
| RIVER | 2.81e-05 | 8.72e-02 | 1.01e-01 | 1.29e-03 |
| VE-diffusion | 5.01e-04 | 3.16e-01 | 3.50e-01 | 3.25e-03 |
| VP-diffusion | 1.56e-04 | 1.36e-02 | 2.61e-01 | 1.06e-01 |
| Stochastic interpolant | 7.25e-06 | 1.03e-01 | 3.45e-02 | 4.61e-04 |
| Ours ($\sigma$ = 0.01, $\sigma_{sam}$ = 0, RK4) | 1.54e-06 | 7.30e-02 | 3.54e-03 | 2.60e-04 |

**Table 9: Standard deviation results for the diffusion-reaction equation task using different choices of probability paths for flow matching. Results are averaged over 5 generations.**
| Model | Test MSE ($\downarrow$) | Test RFNE ($\downarrow$) | PSNR ($\uparrow$) | SSIM ($\uparrow$) |
| --- | --- | --- | --- | --- |
| RIVER | 4.93e-04 | 2.27e-01 | 7.01e-01 | 1.82e-02 |
| VE-diffusion | 6.65e-03 | 1.04 | 5.59e-01 | 9.11e-03 |
| VP-diffusion | 6.75e-04 | 5.67e-01 | 1.62 | 2.13e-02 |
| Stochastic interpolant | 7.92e-04 | 1.14 | 9.58e-02 | 1.20e-03 |
| Ours ($\sigma$ = 0, $\sigma_{sam}=$ 0, RK4) | 1.02e-06 | 4.76e-02 | 7.96e-03 | 1.11e-04 |

**Table 10: Standard deviation results for the Navier-Stokes equation task using different choices of probability paths for flow matching. Results are averaged over 5 generations.**
| Model | Test MSE ($\downarrow$) | Test RFNE ($\downarrow$) | PSNR ($\uparrow$) | SSIM ($\uparrow$) |
| --- | --- | --- | --- | --- |
| RIVER | 5.79e-05 | 2.46e-01 | 6.14e-01 | 1.13e-02 |
| VE-diffusion | 1.21e-04 | 1.36 | 8.82e-02 | 1.12e-03 |
| VP-diffusion | 2.55e-06 | 1.42e-01 | 8.28e-02 | 6.43e-04 |
| Stochastic interpolant | 2.03e-07 | 1.29e-01 | 9.54e-03 | 1.95e-04 |
| Ours ($\sigma$ = 0.01, $\sigma_{sam}$ = 0, RK4) | 5.76e-08 | 4.88e-02 | 9.32e-03 | 7.42e-05 |

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: If the encoder is trained separately from the flow matching model, what objective function is used to train the encoder?
soundness: 2
strengths: 1. The paper is well-written and accessible.

2. It presents a unified framework for probability path models in the context of flow matching and diffusion models, providing readers with a clearer overview of current generative model research.
summary: This paper addresses spatio-temporal forecasting using latent flow matching. In spatio-temporal forecasting, the performance of flow matching is highly dependent on the choice of the probability path model. To address this, the authors propose a new probability path model that accounts for the inherent continuity and correlation in spatio-temporal data, aiming to shorten the interpolating path.
weaknesses: 1. The main contribution of the paper is a new probability path model for spatio-temporal data. However, the motivation behind this model is presented somewhat vaguely. A more formal and rigorous mathematical illustration would improve clarity.

2. There are numerous existing diffusion model-based methods for spatio-temporal data; a comparison with these methods would strengthen the paper.

3. The baselines all seem to use an encoder. An ablation study demonstrating the benefits of modeling in the latent space would be valuable (Please correct me if I missed this).

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - How are those $s_n$ determined in the sampling algorithm?
- What is the reason for setting the highest variance at the middle and lowest variance at both the start and the end?
- Does the proposed probability path have equivalent form for other variants of diffusion models and bring the same benefits to them, such as the score matching and noise prediction objective?
soundness: 3
strengths: - The paper clearly and thoroughly discusses and compares the various kinds of probabilistic forecasting models, which are informative and insightful.
- The proposed probability path model learns to connect consecutive time series samples, leading to faster convergence and more stable training.
- The experiments are intensive, providing support for the proposed model.
summary: In this paper, the authors discussed the influence of probablity path choice in the context of spatial temporal forecasting. The authors conclude that different probability paths can significantly impact the accuracy and convergence of forecasting models and propose a new probability path model specifically designed for probabilistic forecasting of dynamical systems. Experiments show outperforming performance and faster convergence of the proposed method.
weaknesses: - Some notations are a little bit confusing. For example, do $v_t$ in algorithm 1 and $v_s$ in algorithm 2 mean the same vector field?
- The result for faster convergence and fewer sample steps comes empirically. It can benefit from more therotical derivations or insights for why the proposed probability path yileds better results.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: The authors considered Gaussian probability paths, which are suboptimal fore real world data. Have the authors considered elliptical probability distributions and elliptical mixture models, to fix this issue
The authors used second-order performance metrics, yet the approach is probabilistic and requires probabilistic estimates rather than second-order measures. Which higher-order probabilistic metrics would the authors use in this context.
soundness: 2
strengths: THe ability to perform spatio-temporal forecasting is a major strenght.
A novel theoretical framework for flow-based forecasting of spatio-temporal data. The probabilistic path is parametrized using a neural network, which results in the task to be solved in the form of second-order regression.
By taking into account inherent correlations in spatio-temporal data, the proposed model improves upon the existing models of the kind.
Extensive performance records, and comprehensive ablation study.
summary: The authors examine the impact of the specific choice of the probability path model on the forecasting performance. This is achieved in the context of flow matching, a well know principle which has not been properly examined in forecasting. There, the mappings are learned via stochastic processes, through random differential equations.
weaknesses: The main weakness is the choice of the framework. Stochastic processes and ODE based method are known to underperform for non-Gaussian distributions and non-stationary data.
Especially, the gaussian probabiliy paths cannot deal with real-world data which exhibit fat-tailed distribution.
The conditional probability paths and the  marginalization of distributions are rather standard in GenAI.
The performance metrics are all second-rder (MSE, Frobenius norm, peak signal to noise ratio, etc), which are subobptimal given the  use of probabilistic models. 
Too many concepts are put together: encoders, neural networks (MLP),  regression using MSE, optimal mass transport, diffusion, ODEs, Gaussian models.
Overall, the approach appears more of a "system-building" exercise rather than a deep new framework.

Meta-review: The paper introduces a new probability path model under the framework of latent flow matching to improve forecasting performance for spatio-temporal data. The novelty resides in implementing a new probability path for flow matching within the latent space of a pre-trained autoencoder. Considering the existing latent diffusion models and the close relationship between Gaussian flow matching and Gaussian diffusion, the innovation's significance may appear marginal. The distinction of this new model's approach hinges on demonstrating that the choice of probability path significantly impacts results differently than varying noise scheduling in latent diffusion models. However, the paper lacks a thorough theoretical analysis that explains why this specific probability path is optimally suited for spatio-temporal data. Additionally, there is a notable absence of comparisons with relevant diffusion model baselines. These omissions weaken the case for the proposed model's substantial advancement over existing methods.

Decision: Reject