Flow Matching with Gaussian Process Priors for Probabilistic Time Series Forecasting

We thank the reviewer for their comprehensive feedback. In the following, we address their concerns.

Comment: The model is only discussed and validated on univariate time-series.
Response: We have evaluated our model on univariate time series as they are present in many real-world problems, and, more importantly, metrics specifically designed for multivariate forecasting, such as the CRPS-sum, have several limitations, particularly when it comes to evaluating the performance of individual dimensions. As shown in [1], the CRPS-sum can obscure the performance of the model on each dimension due to its aggregation process.

Technically, TSFlow is not limited to univariate time series and can be extended to the multivariate case by adjusting the neural network $u_\theta(t,x_t)$ to handle multiple features, i.e., $u_\theta(t,x_t):[0,1]\times\mathbb{R}^{K\times L}\to\mathbb{R}^{K\times L}$.

For this, one could use a channel-independent architecture such as PatchTST [2] or employ a feature interaction layer (e.g., linear layer or a transformer) across the feature dimension to process the multivariate time series.

Comment: While there are many proposed options in the model, there lacks an ablation or a comprehensive comparison of different choices. Neither is there any theoretical insights in the proposed model.
Response: We have extended our initial ablation in Appendix B.1. We ablate different components of our conditional model, number of neural function evaluations, compare our unconditional model with and without conditional prior sampling, and finally, we compare the asymmetric Laplace distribution against a Gaussian as guidance score. We are happy to include further ablations if requested.

Comment: The presentation of the model needs a bit more clarification. The problem formulation can potentially be expanded a bit more (see questions below). The discussion of different training/inference choices is a bit dense. Maybe some concise pseudocode or flowchart in the main text could be helpful.
Response: We follow the reviewer's suggestions. We have addressed the questions in the following and included a pseudocode of the unconditional generation in the main text. If the reviewer wishes, we will also include the algorithm of the conditional training or conditional prior sampling.

Comment: If my understanding is correct, the overall idea of the model is to push an initial GP to the target time series. Therefore, the sequence length $L$ in the problem formulation section is equivalent to the vector field dimension $d$ in the background section. Can the author(s) kindly confirm if this is true? It would be nice to clarify the setting a little bit. I think the main source of confusion is that there are two time indices in this paper: the time in flow-matching and the time in the time series, and these are orthogonal to each other. I think this should be made clear somewhere in the paper.
Response: The reviewer is correct that the sequence length $L$ corresponds to the vector field dimension $d$ in the background section. To address this, we have clarified the setting in the manuscript. Specifically, we now explicitly state that the flow-matching time $t$ and the time-series index are orthogonal, and we have updated the notation by using $\tau$ to represent the time indices within the time series time. Additionally, we have clarified that $L$ corresponds to $d$ in the background section.

Comment: Missing Experiments on Guided Generation (3.1.3): I could not find experiments on Guided Generation (3.1.3) in the paper. Where can I find them?
Response: The results are shown in Table 3 under the name TSFlow-Uncond. We apologize for the confusion and have clarified this in the text. Additionally, we show the results of the guided generation without conditional prior sampling in Table 8 in Appendix B.1.2 and ablate the guidance score function in Appendix B.1.3.

Comment: LPS Score Clarification: The Linear Predictive Score (LPS) is not sufficiently explained. Please provide more details on how the LPS is calculated, specifically on what the linear model is regressed against and what the output represents.
Response: The Linear Predictive Score (LPS) is computed as follows: First, synthetic samples $\mathbf{y}=(\mathbf{y}^p,\mathbf{y}^f)$ are drawn from the generative model. A linear regression model is then trained to map the past of the time series ($\mathbf{y}^p$) to its future values ($\mathbf{y}^f$), using the true future values for supervision. Since it is a linear regression model, the parameters can be computed in closed form without requiring iterative optimization. Finally, the LPS is then computed as the CRPS of the linear model's forecasts on the real test set. We have provided additional details in Appendix A.5 and added a reference to this explanation in the main text.

Comment: GP Hyperparameter Optimization: Can you perform a small experiment to evaluate the difference in performance with/without GP hyperparameter optimization?
Response: We followed the reviewer's recommendation and conducted the requested experiment. For this evaluation, we used the periodic kernel, parameterized with a learnable parameter $\gamma$: $K_{\mathrm{PE}}(\tau, \tau') = \exp(\gamma\sin^2(d))$

We initialized $\gamma=1.0$ and optimized it by maximizing the likelihood using 10,000 samples from the real datasets. The optimization was performed for 1,000 iterations with Adam. In the following table, we report the $W_2^2$ distances of the generative model trained with the priors under both fixed and fitted values of $\gamma$.

4 NFEs:

16 NFEs:

Across these experiments, we observe an improvement in 6 out of 10 cases when the kernel parameter $\gamma$ is optimized. Thus, while hyperparameter optimization can lead to better performance in some settings, it is not necessary for the generative model.

Again, we thank the reviewer for their feedback and are happy to address any remaining concerns.

We thank the reviewer for their valuable and comprehensive feedback. In the following, we address their concerns.

Comment: Difficulty in Parsing for Non-Experts. The paper assumes substantial familiarity with flow matching and related generative methods, making it challenging for readers without a deep background in these specific techniques. It took for me considerable time to fully grasp the concepts, which suggests that the paper might also be difficult to read for a broader audience.
Response: We thank the reviewer for the feedback. To make the paper easier to read, we have added more information to our methodology and included a pseudocode of the unconditional training in the main text. We are happy to include further information if requested.

Comment: Lack of Runtime Analysis [...] The paper should clearly state the theoretical runtime complexity of using GP priors and provide empirical runtime comparisons in the experiments. [...].
Response: The theoretical complexity of Gaussian Process priors is $\mathcal{O}(L^3)$, where $L$ denotes the length of the time series. In the unconditional setting, this complexity arises from the Cholesky decomposition $(\mathbf{M}\mathbf{M}^T = \Sigma$, note that we used $\mathbf{M}$ instead of $\mathbf{L}$ to avoid ambiguity) needed to sample from the GP distribution. In the conditional setting, it is due to the inversion of the covariance matrix in Eq. 13. However, both operations only need to be performed once, as the covariances are fixed and do not change throughout the sampling process. Furthermore, $L$ is relatively small in our experiments (at most 360), making these computations minor.

During training, we can sample efficiently using matrix-vector multiplications. In the unconditional setting, this involves computing $\mathbf{x_0} = \mathbf{M} \mathbf{z}$ , while in the conditional setting, we compute $\mu_{f|p} = \Sigma_{fp} \Sigma_{pp}^{-1} \mathbf{y}$ (see Eq. 13). Both operations have a complexity of $\mathcal{O}(L^2)$, which is substantially lower than the initial setup.

We include the runtimes in Appendix A.8 and, for convenience, also present them here. The table shows the training runtimes (seconds per epoch) and evaluation runtimes (seconds for generating all 100 forecasts on 959 time series) for the solar dataset.

As shown in the table, the computational overhead introduced by the non-isotropic priors is minimal, demonstrating that these priors do not substantially impact overall runtime.

Comment: Baseline Comparison A simple baseline for the forecasting tasks could involve using Eq. 6 but with an isotropic Gaussian prior. Please include this method as a comparison partner in Table 3. If it is not a valid approach, it should be explained why such a baseline is excluded.
Response: We include a comparison to an isotropic Gaussian prior in Appendix B.1.1. (see Basemodel in Table 7). If the reviewer wishes, we will include this baseline in Table 3.

Comment: Inconsistent Findings on Kernel Choice The periodic kernel minimizes the Wasserstein distance in Figure 2, suggesting it aligns well with the data distribution. However, in Table 1, the periodic kernel does not significantly outperform other kernels [...]. This inconsistency is counterintuitive and warrants further discussion. Moreover, the necessity for different hyperparameter choices across tasks (e.g., generative modeling vs. forecasting) weakens the "one model for all" argument, suggesting more task-specific tuning may be required.
Response: We would like to clarify that to parametrize the Gaussian processes, we chose $\ell=\sqrt{\frac{1}{2}}$, $\ell=1$, and $\ell=\sqrt{2}$, for the SE, OU, and PE kernels, respectively, for all experiments. The reason for this is that they simplify the kernels even further to: $K_{\mathrm{SE}}(\tau, \tau') = \exp(-d^2)$, $K_{\mathrm{OU}}(\tau, \tau') = \exp(-\mid d\mid)$, and $K_{\mathrm{PE}}(\tau, \tau') = \exp(\sin^2(d))$, which are used for all experiments, including generative modeling and forecasting.

Furthermore, the periodic kernel outperforms other methods in downstream tasks, as quantified by the LPS, consistent with the insights presented in Figure 2.

Comment: Could the model be applied to general conditional generation tasks, e.g., imputation. Have the authors tried anything beyond forecasting?
Response: Our methodology is not technically limited to forecasting and can be adapted to general conditional generation tasks, including imputation. However, we chose to focus on forecasting in our experiments due to the broader availability of standardized benchmarks and to ensure fair comparisons with autoregressive models, which do not benefit from bidirectional information passing.

We conducted a forecasting with missing values experiment to explore applications beyond standard forecasting. Following the setup in [4], we evaluated three scenarios: random missing values, blackout missing at the beginning of the context window, and blackout missing at the end of the context window. We masked 50% of the values in each scenario and assessed the forecasting performance using CRPS. Below, we report results for TSFlow (Cond.) and TSDiff. For TSDiff, we report the minimum performance across both versions.

Random Missing:

Blackout Missing Beginning:

Blackout Missing End:

Overall, TSFlow outperforms both versions of TSDiff in 12 out of 18 cases, demonstrating its robustness across these scenarios. If the reviewer wishes, we will include these results in the paper.

We again thank the reviewer for their feedback and are happy to address any remaining concerns.

[1] Yu, Keming, and Rana A. Moyeed. "Bayesian quantile regression." Statistics & Probability Letters 54.4 (2001): 437-447.
[2] Lim, Bryan, et al. "Temporal fusion transformers for interpretable multi-horizon time series forecasting." International Journal of Forecasting 37.4 (2021): 1748-1764.
[3] Gouttes, Adèle, et al. "Probabilistic time series forecasting with implicit quantile networks." arXiv preprint arXiv:2107.03743 (2021).
[4] Kollovieh, Marcel, et al. "Predict, refine, synthesize: Self-guiding diffusion models for probabilistic time series forecasting." Advances in Neural Information Processing Systems 36 (2024).
[5] Tashiro, Yusuke, et al. "Csdi: Conditional score-based diffusion models for probabilistic time series imputation." Advances in Neural Information Processing Systems 34 (2021): 24804-24816.
[6] Alcaraz, Juan Lopez, and Nils Strodthoff. "Diffusion-based Time Series Imputation and Forecasting with Structured State Space Models." Transactions on Machine Learning Research.
[7] Biloš, Marin, et al. "Modeling temporal data as continuous functions with stochastic process diffusion." International Conference on Machine Learning. PMLR, 2023.
[8] Tamir, Ella, et al. "Conditional flow matching for time series modelling." ICML 2024 Workshop on Structured Probabilistic Inference {&} Generative Modeling. 2024.
[9] Kerrigan, Gavin, Giosue Migliorini, and Padhraic Smyth. "Functional flow matching." arXiv preprint arXiv:2305.17209 (2023).

We thank the reviewer for their valuable feedback. In the following, we address their comments and questions.

Comment: There are missing related works [8,9].
Response: We thank the reviewer for pointing out the references. We have included them in our related work section.

While there are multiple differences to [8], the most notable one is that [8] models a dynamical system using CFM, where the time series and ODE share the time dimension. In contrast, TSFlow introduces two orthogonal time dimensions: one along the generative process and another within the time series itself. This allows TSFlow to employ a sequence-to-sequence model across all generation steps.

As for [9], it extends the flow-matching framework to infinite-dimensional spaces and also leverages Gaussian processes. However, their key-focus is different. More specifically, [9] constructs conditional Gaussian measure paths, whereas TSFlow focuses on parameterizing the joint distribution aligning with its sequence modeling framework. This leads to multiple differences, such as our Gaussian process regression and conditional prior sampling or the choice of our architecture.

Comment: The use of equation 9 was a bit unclear to me. Why is this specific form of $q_1$ chosen?
Response: We use this specific form of equation 9 due to its connection with Bayesian quantile regression [1]. Applying the log-likelihood (see Eq. 8) leads to the quantile loss. The quantile loss is a CDF-based score and common practice in probabilistic forecasting [2, 3, 4].

We have added a section in Appendix A.3 discussing this choice. Furthermore, we have added an experiment in Appendix B.2, comparing it to a Gaussian distribution, which yields inferior results.

Comment: In Section 3.1.2, I am guessing that once we sample $x_0 \sim q_0(x_0∣y^p)$, then we use $x_0$ as an initial condition for the flow model to generate new samples $y\mid x_0$. Is this the case? If so, it would help to state this explicitly somewhere in the paper.
Response: The reviewer is correct that conditional prior sampling involves a two-step process. We have clarified this procedure in Section 3.1.2 to make the sampling steps more explicit.

Comment: In Line 303, the authors write "approximating $q_0(x_0 \mid y^p)$ with $q_0(x_0 \mid y^p)$". I think I understand what is meant here, but this seems to be a typo.
Response: The reviewer is correct. The correct phrase should be "factorizing $q_0(\mathbf{x}^p_0 \mid \mathbf{y}^p)$ with $q_0(\mathbf{x}^p_0 \mid \mathbf{y}^p)\,q_0(\mathbf{x}^f_0 \mid \mathbf{y}^p)$", as detailed in the "Gaussian Process Regression" paragraph. We have corrected this typo.

Comment: In lines 154-156 the authors describe a time series as a vector in $\mathbb{R}^L$. Does this mean that the authors only work on time series having a fixed length $L$, or does the method allow for variable-length time series? Similarly, are the authors assuming that the time series all share a fixed discretization (i.e., there are some fixed times $(t_1,\dots,t_L)$ corresponding to $(y_1,\dots,y_L)$, or can the discretization vary across time series? Is this discretization assumed to be uniform, or can it be irregular as well?
Response: Our experiments considered time series of fixed length $L$ with fixed uniform discretizations in hourly and daily intervals (see Tab. 4 in the appendix). This is consistent with previous works [4,5,6]. In theory, TSFlow can be trained on time series of varying lengths $L$. Applying TSFlow to irregular time series would require a slight adaption of the network architecture. More specifically, one would need to include positional embeddings of the time steps as input to the neural network, as done in [7]. We have clarified in the paper that we use regularly sampled time series.

We thank the reviewer for their valuable feedback and comments. In the following, we address their remarks.

Comment: Can the model extend to the multivariate time series problem?
Response: The model can be extended to the multivariate time series problem by adjusting the neural network $u_\theta(t,x_t)$ to multivariate time series, i.e., $u_\theta(t,x_t):[0,1]\times\mathbb{R}^{K\times L}\to\mathbb{R}^{K\times L}$.

To extend the neural network, one could use a channel-independent architecture such as PatchTST [1] or employ a feature interaction layer (e.g. transformer) to process the multivariate time series.

We have decided to perform experiments on univariate data as these present many real-world problems. Furthermore, metrics such as the CRPS-sum specifically designed for multivariate forecasting have several limitations, particularly when evaluating individual dimensions' performance as shown in [2].

Comment: Why does closedness of the prior and data distribution imply easy learning. Do you have any experiments about train efficiency or path efficiency?
Response: Non-isotropic priors introduce patterns that naturally occur in real data, thereby preserving the non-i.i.d. characteristics during the noising/forward process. This introduces an inductive bias throughout the entire generation process. In contrast, the isotropic prior starts with initial samples without interdependencies, requiring the model to learn and introduce these during generation.

We have included Figure 7, which shows the Wasserstein distances of the synthetic samples across different numbers of neural function evaluations (NFEs). The non-isotropic priors demonstrate faster convergence compared to the isotropic prior.

During training, we observed that the non-isotropic priors reached the performance levels of the isotropic prior faster. More specifically, to achieve the best Wasserstein distances (at 16 NFEs) obtained by the isotropic prior, the non-isotropic priors reduce the number of training iterations by nearly 50%. The table below summarizes the percentage of training iterations required by each non-isotropic prior, relative to the isotropic prior, to reach this performance.

Note that we excluded the few cases where the isotropic prior outperformed the non-isotropic priors.

Comment: Can the given prior (Gaussian process) extend to the arbitrary prior?
Response: We thank the reviewer for pointing us to this reference. Yes, the prior in TSFlow can be extended to arbitrary source distributions. Importantly, we do not require access to the PDF or CDF of the distribution; we only need to be able to sample from it. The adaptation only involves modifying the distribution $q_0$, while $q_1$, $p_t$, and $u_t$ remain unchanged.

We included this in the updated manuscript and cited the reference.

Comment: Do you have any theoretical evidence about how the selection of kernel functions effect to the model performance?
Response: Since most datasets consist of a mixture of components, such as periodicity, trends, and noise, it is challenging to assign specific kernels to datasets definitively. However, our results indicate that datasets with a dominant periodic component (e.g., solar, electricity, and traffic) achieve better LPS scores when using a periodic prior (see Table 2). Conversely, on datasets like KDDCup, which are dominated by non-periodic components, the periodic prior struggles to deliver strong performance.

We hope that the questions are answered satisfactorily and are happy to address any upcoming concerns.

[1] Nie, Yuqi, et al. "A Time Series is Worth 64 Words: Long-term Forecasting with Transformers." The Eleventh International Conference on Learning Representations.

[2] Koochali, Alireza, et al. "Random noise vs. state-of-the-art probabilistic forecasting methods: A case study on CRPS-Sum discrimination ability." Applied Sciences 12.10 (2022): 5104.