Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Authors: Adrien Cortes
Remi Rehm
Victor Letzelter
Venue: ICML 2025 poster (2025)
Keywords: Time-series forecasting
Quantization
Probabilistic methods
Conditional Distribution Estimation
Winner-takes-all
Diversity
Multiple Choice Learning
Primary area: deep_learning->sequential_models_time_series
TL;DR: We introduce TimeMCL, leveraging the Multiple Choice Learning (MCL) paradigm to forecast multiple plausible time series futures.

Abstract: We introduce $\texttt{TimeMCL}$, a method leveraging the Multiple Choice Learning (MCL) paradigm to forecast multiple plausible time series futures. Our approach employs a neural network with multiple heads and utilizes the Winner-Takes-All (WTA) loss to promote diversity among predictions. MCL has recently gained attention due to its simplicity and ability to address ill-posed and ambiguous tasks. We propose an adaptation of this framework for time-series forecasting, presenting it as an efficient method to predict diverse futures, which we relate to its implicit *quantization* objective. We provide insights into our approach using synthetic data and evaluate it on real-world time series, demonstrating its promising performance at a light computational cost.

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Adrien Cortés
Rémi Rehm
Victor Letzelter
Abstract
place an explicit distribution on the model's output and per- form maximum likelihood estimation (Salinas al., 2020; Alexandrov et al.. 2020). While such parametric methods can be computationally efficient, they may depend heavily the choice of distribution family, reducing their flexibility in capturing complex uncertainties (Gneiting & Katzfuss, 2014).
We introduce T TimeMCL, a method leveraging the Multiple Choice Learning (MCL) paradigm to forecast multiple plausible time series futures. Our approach employs a neural network with multiple heads and utilizes the Winner-Takes-Al (WTA) loss promote diversity among predic- tions. MCL has recently gained attention due its simplicity and ability to address ill ill-posed and ambiguous tasks. We propose adaptation of this framework for time-series forecasting, pre- senting it as an efficient method to predict diverse futures, which we relate to its implicit quanti- zation objective. We provide insights into our approach using synthetic data and evaluate on real-world time series, demonstrating its promis- ing performance at a light computational cost.
In parallel, the success of general-purpont generative models, especially conditional diffusion models (Ho et al., 2020) such as TimeGrad (Rasul al., 2021a), has led to strong empirical performance in high-dimensional time series forecasting. However, these models tend to be computationally expensive at inference, particularly when multiple what-if scenarios, which we will refer to as hypotheses, need to be generated in real-time. Moreover, there is often no explicit mechanism to guarantee sufficiently diverse hypotheses within a single model pass.
To address these limitations, propose TimeMCL- novel approach based on Multiple Choice Learning (MCL) techniques- that produces diverse and plausible predictions via a single forward pass.
1. Introduction
Predicting the weather of the upcoming weekend or the stock prices of next month with perfect accuracy would un- doubtedly be useful. Unfortunately, time-series forecasting a highly ill ill-posed problem. In many cases, the available input information insufficient reduce our uncertainty about the estimation of the underlying stochastic process and the data itself may contain noise. Consequently, the best forecaster can do is estimate plausible future trajectories, ideally along with the probability of each outcome.
Contributions.
We introduce TimeMCL, new approach for time series forecasting that adapts the -inner-Takes-AIl (WTA) loss to generate multiple plausible futures.
We show that TimeMCL can be viewed as a functional quanfizer, and we illustrate its theoretical properties on synthetic data.
Because temporal data are highly structured and typically come with imput-outpul pairs that require no additional man- ual annotation, autoregressive neural networks have become the de facto standard for forecasting high-dimensional time series from historical data and exogenous covariates (Ran- gapuram et al., 2018; Salinas et al., 2019; Benidis et al., 2020). To capture predictive uncertainty, practitioners often
We evaluate our method on real-world benchmarks, demonstrating that TimeMCI, efficiently produces a diverse of forecasts with just few samples in single forward pass.'
'Equal contribution 'LTCI, Télécom Paris, Institut Polytech- nique de Paris, France VValeo.ai. ai, Paris, France. Correspondence to: Adrien Cortés <ad. cortes43 gmail.com>, Victor Letzelter <letzeltee.victorr frr.
2. Related Work
Ambiguity and need for diversity time series forecast- ing. In recent years, scquence-to-sequeee neural networks (Hochreiter & Schmidhuber, 1997; Sutskever al., 2014;
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
"Code available http:///ghhaaccoom///eettttttttt

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Chung et al., 2014; Torres et al., 2021) have become increas- ingly effective in time series forecasting, often surpassing classical techniques (Hyndman et al., 2008; Hyndman & Khandakar, 2008). Yet capturing the inherent ambiguity in future outcomes remains a critical challenge, especially in high-dimensional settings (Ashok et al., 2024). Salinas et al. (2020) proposed a probabilistic autoregressive global model capable of fitting and forecasting high-dimensional time series while highlighting the need to quantify uncer- tainty in predictions. Building on this line of research, Rasul et al. (2021a) introduced conditional diffusion model that summarizes the past values of the time series into a hidden state, then performs a diffusion process-conneee conditioned on this state-to generate forecasts. Rasul et al. (2021b) lb) re- tained the conditional architecture but replaced the diffusion mechanism with a normalizing-flov generator. While these methods are capable of modeling uncertainty, computational efficiency remains crucial factor, particularly in real-time scenarios (Chen & Boccelli, 2018) where multiple plausi- ble futures must be generated. To address these challenges, general-purpose we introduce a new family of general purpose autoregres-
omitted in the following for conciseness. Our focus is on scenarios where the conditional distribution may exhibit multiple modes (mulii-modalityy). motivating richer repre- sentation than single-mean regressor.
To address this issue, the goal of probabilistic time series forecasting is to capture conditional distributions over future time series given past values with model po. with parameters 0, whose likelihood can be expressed as
T II
po 11:t0-11
pe (ate T1a-1). t=to
Once trained, ancestral sampling methods can be used in- fer sequence-level predictions. Let us illustrate this scheme using for instance networks (RNNs) from Graves (2013).
based recurrent neural
When considering hidden-variables- based models, the ba- sie building block of sequence-to- to-seqeence architectures (Sutskever et al., 2014), one often implicitly parametrize the model with, up to a (log) normalization constant,
sive time series forecasters based on the Alll (WTA) principle, leveraging its quantization properties to produce diverse and realistic forecasts in single forward.
log po( Fto:t-1, 11:0-111 = -((fo(he-)), xx). (3)
where ((:-) be an appropriate loss, e.g., the mean squared error f(a, 3) & | and all the context and previ- ous states is encapsulated into hidden state h4-1 E H. Assuming the vanilla form for recurrent networks, the hid- den state propagation is represented by model so with h4-1 = so (2 h2- and fo is the final projection.
Optimal vector quantization & Multiple choice learn- ing for conditional distribution estimation. Quantization is concerned with finding the best fimitely supported ap- proximation of a probability measure (Bennett, 1948; Du al., 1999; Pagès, 2015; Chevallier, 2018). In the context of time series forecasting, this translates quantizing the conditional probability distribution of the target time series. Multiple Choice Learning (MCL) with Winner-Takes- (WTA) loss Gouzman-Rivers et al., 2012; Lee et al., 2016) provides a practical framework for such conditional quan- tization through multi-head networks, which act as a fixed set of codevectors (also called hypotheses) (Rupprecht et al., 2017; Letzelter et al., 2024; Perera et al., 2024). While MCL has thus far been explored in various applications, notably computer vision tasks (Lee et al., 2017; Rupprecht et al., 2017: Tian et al., 2019), we adapt it here to predict a quantized representation of the conditional probability distribution of future time series values, using a training scheme specifically tailored for this setting.
Once trained, the predictions can then be performed by first encoding the past sequence F1:6--1 with a hidden state by applying recursively sg: heg se (2 (,g.1,,,.,,,((,,,,,. (a(2,,,s(22) s0 ho) .......,,,, )), where ho is ar- bitrary initial hidden state. Then the recurrent model can be unrolled, i.e., turned into autoregressive mode by decod- ing the predicted sequence by applying recursively 89, this time over its own predicted outputs, leading the prediction ftoT ~ po (2 Too:T-1, h4,). In the following, we will denote Fo F1 + #to:T F0(11:to-1) 22' the unrolled network.
4. TimeMCL method
TimeMCL leverages the Winner-Takes-All principle (Guzman-Rivera et al., 2012; Lee et al., 2016), which was originally introduced to address ambiguous tasks. WTA naturally extends scenarios in which future time-series trajectories exhibit multimodality (e.g., seasonality, regime switches, and abrupt events). We build on the Multiple Choice Learning framework, which produces K distinct "hypotheses" through multiple heads. Not only does our esti- mator allow us to produce plausible hypotheses; T TimeMCL provably quantizes the target distribution of futures, and
3. Problem setup and notations
Let (IEr) E xT++1 represent multivariate time series on X' = RDD over time indexes [tt, TI. where i E [1, T1. We aim to learn the conditional law
p( (1:to-1,
(1)
of future values T (over the interval [to, TD) given past observations - and covariates C1:T, the latter being
2

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
is therefore expected to infer the K most representative predictions of the target distribution.
where the binary cross entropy BCE(p, g) - - plog(()
(1 p) log(1 q), aligns the predicted and target winner assignation probabilities. The final training objective is compound loss L = L + BLS, where B > 0 is the confidence loss weight. See Figure for an illustration of the components of TimeMCL.
4.1. Training scheme
A key insight behind WTA is that learning K separate hy- potheses _t, ....., _x E to with objective that ef- fectively induces tesselation (Du et al., 1999) of the target space into K cells (one for each hypothesis) aims to capture the best possible information of the target distribution with K points.
R
TimeMCL works as an alternative the vanilla maximum likelihood estimation (MLE) of (2). Let pb. Pg.-...P7" r be K models with parameters (01, 1.1,.--3 0x) ), for which one can associate heads .a.....ff SKK H X' using the hidden- state representation as in Section 3. In our implementation, the K models have shared So for hidden-state propagation, and differ only by their final heads f and one may define the complete models with f 8g.
Figure I. Components of TimeMCL. Prediction heads are in Lightcoral color, and Score heads are Blue. Rectangles contain functions, and circles contain features.
The
consists of the following training scheme for each data point (21::0-1,
4.2. Inference and Sampling
Once trained, T TimeMCL provides K plausible predictions a ~ P!..... &K ~ P5 by ancestral sampling on each of the models. When using recurrent neural networks, the unrolling procedure descrihed in Section 3 can be applied by first encoding the input sequence 8q to obtain a hidden state he, and unrolling the autoregressive model. As in Section 3, we encapsulate these operation with unrolled networks 3 Fg...., FKK The scores are computed in the same way, using the score heads 7 instead of the prediction heads fă. and we denote r their unrolled networks. To get a single score associated with each predicted trajectory, we averaged the predicted scores over the sequence.
1. We compute the negaaive-log-likeliood of each model
Σ t=
log (2, x1a-11)
where log PV (2ze x14-11) = - e($ (),), for each head k E [1, KT.
2. We pick the "winner" k* = arg ming Cb. and we back- propagate only through that winning head (I**).
This two-step optimization allows optimize the loss in an alternating fashion, despite the non-differentiabbiity of the min operator. Note that the latter can be computed batch- wise on the historical data, computing the Winner head of each batch index, with a loss that writes as:
In cases where the ambiguity reduced, such that for short- horizon forecasts when only one prediction is required, one might pick the best head according to the predicted score or sample from them in proportion to some confidence measure. In a longer autoregressive forecast, we can consider the K outputs at each time step -thus producing a set of possible futures from the single forward pass.
CWTA (01..... ... 0x) A Exyrl I min ( (tgrr) klI....K
(4)
where the expectation is taken over (Iisto-1 Fto:T) p(Z1:)) While the WTA Loss trains several models with the aim of producing several trajectories, heads as in Letzelter et al. (2023) 7/....,,% % H + 1] to learn to predict the probability of head being the winner and avoid overconfident heads. The latter are trained with
4.3. Taking advantage of Winner-takes-all variants
While the WTA Loss has been proven effective for handling ambiguous tasks (Lee et al., 2016; Seo et al., 2020; Garcia et al., 2021), some heads may theoretically be under-trained (Rupprecht et al., 2017). This may occur when a single mode (or a few modes) dominates the target distribution, or due to suboptimal initialization, similar to what can happen in Lloyd's 's algorithm for K K-Means clustering (Lloyd, 1982; Arthur & Vassilvitskii, 2007). In this case, the scoring loss L ensures setting a low probability to those concerning hypotheses so that the latter can be ignored at inference.
also use score
L A Emar
Σ kEl...KKK
BCE [1 1",,%% 1..%%(me (he-1))
3

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
is possible to mitigate this issue, and therefore improve the performance of the estimator by using relaxation techniques of the min operator. In this case, the best head loss in (4) can be substituted with weighted loss from the different heads:
5.1. TimeMCL is stationary conditional functional quantizer
For simplicity and without loss of generality, let tem- porarily assume that 3 = 0, i.e., only the WTA Loss LWTA is optimized. When predicting the future of time series given its context -1 we are effectively observing, dur- ing supervision, a (partial) path realization Ttg:T of an un- derlying stochastic process, Following standard functional data analysis (Bosq. 2000; Ramsay & Silverman, 2005), we assume that each time-series trajectory Tto:T belongs to Banach space (E, |I- For concreteness, one may consider E = (ito. TI) endowed with the usual I norm, though any separable Banach space is admissible for the theoretical arguments to hold. In this space, the distance d(...) induced by || || enables us to define the Voronoi tessellation over future paths Itg:T as described below.
EWTA (0),...,.z) A Bxyr Σ qu Lo (11:1t-11 (to:T)
(6)
where the coefficients qu 2 O sum to one, and allow to assign some weight to non-winner hypotheses.
This idea was originally suggested through Relaxed Winner- takes-all (R-WTA) loss proposed by Rupprecht et al. (2017), which suggested back propagattg not only on the winning head k-, but also on the non-winner. In this case Ok- = 1
and O = K-1 for k + k* (See (21) in Appendix C.3). More recently, Perera et al. (2024) proposed an annealed method inspired from Deterministic annealing (Rose al.,
The 1,2 quantization problem aims at finding the best ap- proximation of random vector X, using K points in E. The quality of the approximation using at most K points [fi...., fK }, is generally measured with the distortion, de- fined as: D2 = Ex fe (f(....fx d(X, f) which is finite over X admits a second order moment. Note that our learning setup involves quantization in conditional form, i.e. the random variable of interest depends context.
1990) using a softmin operator:
1 Tto:
(T) Z ( 11 1, ::o::: exp(
T
with
We state hereafter our main theoretical result, showing that T TimeMCL can provably perform functional quantization of the target space of plausible trajectories. It can be seen as an adaptation of Proposition 5.2 in Letzelter et al. (2024) to the case of functional quantization of time series.
K 2
Lp(11 -1, Fea.T)
exp
T
where the temperature T is annealed during training, egg. considering at training step f, T(t) T(0)p² with p < At higher temperatures, the target assignment is effectively softened, making the carly stages of training casier.
Proposition 5.1. See Proposition A.4 in the Appendix. Un- der the assumptions that:
1. The batch size is big enough so that the difference between the WTA risk and its empirical version can be neglected (Assumption A.
We implemented these two variants, within our T TimeMCL model in the Gl uonts framework (Alexandrov et al., 2020). Based on our experience, these methods are mean- ingful variations of WTA that were worth exploring, as they demonstrated improvements over vanilla WTA in certain configurations.
2. The neural network we are considering is expressive enough that minimizing the risk is equivalent to minimizing the input input-dependent risk for each context x14.-1 (Assumption A.2).
5. Theoretical analysis and interpretation
The training has converged and CTA has reached a local minima (Assumption A.3).
In this section, based on the notion of functional quanti- zation, we provide insights into the interpretation of our approach. particular, we show that, under squared error, the K heads form a Voronoi tessellation of future trajectories and act as conditional codebook. This viewpoint explains how WTA theoretically captures the best possible way the conditional law over stochastic processes, given a sampling budget of K predictions. Our claims are then illustrated through a synthetic data example, specifically focusing on certain Gaussian Processes.
Then, TimeMCL is conditional stationary quantizer for each sampled window (11: that is, for each k E [1, K]:
F (0110-11 = E[zwer Ftoe E X (x(2110-1)), (8)
where
X (2 = fooOT cb < LG Vr + +b)

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
We denoted by abuse of notations Ci Ct (1 Too:)) for simplicity and the same for Lig This makes T TimeMCL akin to conditional and gradient-based version of K- Means over the set of possible future trajectories.
5.2. Smoothness of TimeMCL predictions for (most) time series
As we have shown, under certain assumptions, if the model reaches a stationary point.
F (#1: = E [a0. Fio:T E [x((
Sketch of proof. See Proposition A,4 for the full proof. The demonstration of this result is made by first leveraging As- sumption A.1 1 to re-write the WTA Loss in the form
and interpret the prediction as mean of different hypothetical trajectories. Since most stochastic time series contain centered noise-caused by various factors such as measurement errors or random events-the averaging pro- cess tends to eliminate this noise, resulting in a smooth appearance. We observe this phenomenon consistently in the real examples we visualized. We also noticed that the appearance of smoothness increases training progresses. This property of mean predictions reinforces our conviction that the model is providing representative trajectories rather than just random samples.
K Σ
Cb dp(atgrr m1 (9)
where each FroT (1 picks the head k it clos- est to. Under the expressivity assumption (Assumption A.2), (9) comes down to optimizing functional (14) of the hypotheses position for each fixed context 11:ta-1- We assume that, during training, our predictor generates tra- jectories solely from the context (i.e., independently of the observed values), effectively as though F) F K were used directly for training.
5.3. Synthetic data example
We evaluate TimeMCL on three synthetic processes: Brow- nian motion, a Brownian bridge, and an AR(5) process. While the Brownian motion serves as a simple example with minimal context dependence, the Brownian bridge introduces a time-conditioned structure, and the AR(5) pro- cess tests the model's ability handle stronger context dependencies. Training is performed randomly sampled trajectories with appropriate conditioning, and quantization is assessed using theoretical references: Karhunen-Loèwe based quantization for Brownian motion and the Brownian bridge, and Lloyd's ' algorithm for the AR(5) process (Ap- pendix B).
From here, we leverage the fact that TimeMCL is a two-step training procedure, and then the alternating optimization argument from Rupprecht et al. (2017) (Theorem 1) applies to obtain the optimal centroids. This is performed using L square loss for E, from the vanishing gradient condition on the optimized functional (Assumption A.3). We also say that in this case, Voronoi tesselation on the trajectory space induced by the hypothesis is centroidal (Du et al., 1999).
Now if we assume 3 > 0, we can show the following proposition.
For this toy experiment, we used three-layer MLP that predicts the entire sequence at once from its context (see Appendix B for details). To keep the implementation as straightforward as possible, we kept the model's 's parameteri- zation lightweight and omitted the score heads.
Proposition 5.2. Under similar assumptions as Proposi- tion 5.1, one can show that a necessary optimality condition for the score heads is that
As shown in Figure 2, TimeMCL consistently produces smooth and near- rear-optimal quantizations across all settings. For Brownian motion, the predicted trajectories closely align with the theoretical optimal quantization, demonstrating the model's ability to learn conditional distributions from min- imal context. In the Brownian bridge setting, where the conditioning depends on both time and value, the model successfully captures the structural constraints, leading to coherent and well-quantized trajectories. The AR(5) pro- cess presents a more complex challenge due to its stronger temporal dependencies, yet TimeMCL effectively utilizes the past observations to produce long-horizon predictions that remain consistent with Lloyd's quantization. The re- sults highlight the model's ability to condition on past ob- servations, effectively adapting to different processes and maintaining predictive stability over extended time horizons.
1b (1 = (xtg.T E Xu (1) (1) (10)
Proof. Full proof in Appendix A.6.
TimeMCL can thus be viewed as a conditional vector quan- tization scheme (Gersho & Gray, 1992), where each head k a code vector (in functional form).
By conditioning on the past data 21:6, TimeMCL effec- tively learns a family of partitions in the time-series trajec- tory space. If the number of hypotheses is large, increasing K under the above assumptions typically reduces recon- struction error at a rate akin to K 22dd in line with classical quantization theory (Zador, 1982).
5

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
trajectory
Target quantization
Prodictions
Brownian Motion
Brownian Bridge
AR(p)
1:25
0.8
1.11
1.00
0.6
Amplitude
0.75
000
0.0
0.22
-6.5
-0.2
0.00
-0.4
-1.0
00.25
0.6
0.0
0.2
0.44
0.6
0.8
1.0
0.4
0.6
0.8
1.0
0.0 0.2
0.4
0.6
Time
Time
Time
Figure 2. Conditional Quantization of Stochastic Processes with TimeMCI. The Figure shows the predictions of Tí meMCI. three synthetic datasets as described Section 5.3 and Appendix B. Predictions TimeMCL are shown in lightcoral color, and target quantization each those toy experiments. Brownian Motion, Brownian Bridge, and AR(p) are increasingly complex terms of conditioning dependencies in the context window (See Section 5.3). We see that in those three cases, TimeMCL predictions nicely approximate the shape of target quantization justiinn its interpretation conditional functional quantizer.
shown brown. We used 10 hypotheses here, three-layer MLP backbone with score-heads disabled in
6. Empirical evaluation
Metrics. We evaluate our approach with six different met- rics, each of them being detailed in Appendix C.2. First, we considered the Distortion, which is computed as
In this section, we empirically validate our method, with experiments on real-world time series. Compared to the experiments in Section 5.3, the underlying law of the data- generating process is not known, making the task more realistic. The goal is to compare TimeMCL with state-of- the-art probabilistic time series forecasters, emphasizing its balance between quantization, predictive performance, and computational efficiency.
1
Σ min d (2 D2 N Il.....
consisting of computing the Euclidean distance d between the target series and closest hypothesis, averaged over the test set with N samples. This allows to compare fairly with other baselines when the sample size K is fixed. To assess computational efficiency during inference, port both the inference Floating Point Operations (FLOPs) and run-time for each baseline in Table 2. As a means to validate the theoretical claim presented in Section 5.2, we also compute the Total Variation (TV), which quanti- fies the smoothness of the predicted trajectories Table 3. Finally, for comprehensive comparison, erage Root Mean Square Error (RMSE) and the Continuous Ranked Probability Score (CRPS) (summed over all time series dimensions). These results are reported in Tables 7 and 8 in the Appendix.
6.1. Experimental setup
re-
Datasets. Our approach evaluated on six well-established benchmark datasets taken from Gluonts library (Alexan- drov al., 2020), preprocessed exactly in Salinas al. (2019); Rasul al. (2021a). Each dataset consists of strictly positive and bounded real-valued time series. The char- acteristics of these datasets are summarized in Table in the Appendix. SOLAR contains hourly solar power outputs from 137 sites with strong daily seasonality. ElEcTRICITY records hourly power consumption for 370 clients, exhibit- ing daily and weekly periodicities. EXCHANGE tracks daily exchange rates for eight currencies, showing less seasonality and being influenced by macroeconomic factors. TRAFFIC provides hourly occupancy rates from 963 road sensors, cap- turing rush-hour peaks and weekly patterns. TAXI comprises time series of NYC taxi rides recorded at 1,214 locations. WIKIPEDIA includes daily views of 2,000 Wikipedia pages (Gasthaus et al., 2019). See Lai et al. (2018) and Appendix C.I for an extensive description.
include the Av-
Baselines. We considered the following baselines: ETS (Hyndman et al., 2008), Tactis2 (Ashok et al., 2024). DeepAR (Salinas al., 2020), TempFlov (Rasul al., 2021b), and TimeGrad (Rasul et al., 2021a). These were compared against T meMCI. with two relaxation techniques: Relaxed-WTH (Rupprecht et al., 2017) and aMCL (Perera et al., 2024). Note that both of these Multiple Choice Learn- ing variants use score heads as in Letzelter et al. (2023) (with 0.5).
6

---

Wannerltakes-all for Multivariate Probabilistic Time Series Forecasting
TeegiLow
Tactis2
TimeGrad
TempFlow
TimeNCL (annealed)
TimeMCL (relaxed)
Figure 3. Qualitative results. Visualization of the predictions on the SOLAR dataset, comparing Tactis2 (purple), TimeGrad (red), TempFlow (orange), and TimeMCL trained with two relaxation techniques (blue). Each model predicts sixteen hypotheses, dashed lines represent the predicted mean the conditional distribution which was computed as the weighted sum the predicted hypotheses by the scores for TimeMCL. Each row represents different dimension (D 137 here). The light yellow and the score intensities are displayed as shaded blue lines, as per the scale on the right. DeepAR is not included in this visualization as it does perform competitively. In this example. TimeGrad, TempF Low and Tactis2 generate meaningful hypotheses. However, Ti meMC1.'s 's predictions are noticeably smoother and more diverse.
highlights the prediction window,
Architectures. We compare DeepAR, TempFlow, and TimeGrad with T imeMCL, using the same neural net- work backbone: an RNN with LSTM cells, as in the orig- inal implementations (Hochreiter & Schmidhuber, 1997). This ensures fair comparison conditions. In these experi- ments, each hypothesis head in TimeMCL and the projec- tion layer of DeepAR consists of a single linear layer. Mean- while, the noise prediction in TimeGrad is implemented with dilated ConvNet featuring residual connections (Van Den Oord et al., 2016; Kong et al., 2021; Rasul et al., 20) Additionally, we include comparisons with methods based on transformer backbones, such as the transformer-baset version of TempF sempFlow, 1ow, (named Trf. TempF Fow (Rasul et al., 2021b)) and Tactis2 (Drouin al., 2022; Ashok et al. 2024), which leverages copulas for modeling depen- dencies. Note that ETS does not use neural network.
Evaluation. The evaluation dataset is divided into mul- tiple non-overlapping subsets, each containing sufficient points for both context and prediction lengths, allowing comprehensive assessment across different temporal seg- ments. To compute TimeMCL metrics while accounting for predicted hypothesis probabilities, we resample with replacement from the K hypotheses obtained in single for- ward pass, weighting them by their predicted scores before computing the metries.
6.2. Results
Tables I and 2 show Distortion performance, FLOPs, and run-time, comparing T TimeMCL with two tested WTT vari- ants with 16 hypotheses - WTA-Relaxed (Rupprecht et al., 2017) and aMCL (Perera et al., 2024) - against the base- lines. Table 3 displays Total Variation comparison, as a means to quantify smoothness.
Training. Training is conducted using the Adam optimizer with an initial learning rate of 10-3 for 200 training epochs. During each epoch, 30 batches of size 200 are sampled from the historical data, considering random windows with context set equal to the prediction length. We used validation split of size equal to 10 times the prediction length. Except for Tact i s2, which uses Z-Score normalization, the data are scaled by computing the mean value dimension by dimension over the context and dividing the target by this mean. This scaling follows the Ti TimeGrad experimental setup (Rasul et al., 2021a), ensuring consistency. The model is trained the scaled data, and the inverse transformation is applied later for prediction.
Distortion and Computation Cost. Table I demonstrates that T TimeMCL, particularly when trained with its Relaxed variant, achieves competitive performance compared to the other models when the number of hypotheses is fixed (with K = 16 here). This is especially promising, as Tactis s2 and TimeGrad, which are the strongest competitors in terms of distortion, ineur significantly higher FLOPs and run-time (see Table 2). A similar trend is observed in Table 6 in the Appendix, which shows results for 8 hypotheses. This behavior is expected since TimeMCL explicitly op- timizes for distortion. It's worth noting that, among the neural-based methods, TimeMCL is the second most com-

---

Winner-takes-all for Multivariate Probabilistic Time Scries Forecasting
Table Distortion Risk with 16 Hypotheses. TimeMCL (R.) ) and TimeMCL (A. ) correspond to the relaxed and annealed variants. ETS, Trf. TempFlow and Tactis2, columns are in gray because they don't share the same backbone as the other baselines. The test scores are averaged over four training seeds for each model. Best scores are in bold, and second- best are underlined.
ETS
Trf TempFlow
Tactis2
TimeGrad
DeepAR
TempFlow
TimeNCL (R.)
TimeMCL (A.)
ELEC.
14041 + 877
13519 = 2789
11616 = 2160
9872 i 668
133107 = 2870
14836 +: 396
11604 + 1042
11611 = 1206
EXCH.
0.051
0.062 = 10.15
0.030 + 0.002
0.038 0.005
0.061 + 0.000
0.051 0.007
0.035 0.005
0.044 0.005
SOLAR
641.32 3.57
374,69 1. 20.97
358.01 1. 38.47
360.6 34.79
748.68 18.92
371.14 1 18.2
280.03 15.22
305.5 4,05
TRAFFIC
2.64 3.64:20.001
1.26 = 0.02
0.84 + 0.04
0.78 = 0.05
2.12 + 0.08
1.21 = 0.04
0.68 1: 0.01
0.72 = 0.05
TAXI
583.52 0.41
278.22 1: 13.8
243.63 1: 9.1
209.64 1. 3.14
407 1. 14.09
268.65 1: 10.69
187.81 1: 6.17
229, 269.2661 1: 22.06
WIKI
715150 :: 4742
515285
16000
254675 3870
261361 + 1028
368087 = 8159
382554 + 4149
261950 9504
267624 8433
Table 2. Computational cost at inference of neural-based methods for K 16 hypotheses EXCHANGE. FLOPs computed with a single batch of size 1. Run-times are averaged over 15 random seeds.
Trf. TempFlow
Tactis2
TimeGrad
DeepAR
TempFlow
TimeMCL
FLOPs. (1)
2.04 10a
1.85 10
3.05 x 10P
4.65 10
9:29 10"
8.83 10
RUN TIME, (L)
2.17 + 0.23
8.69 + 0.36
211.572200
0.00
1.39 + 0.03
1.12 0.01
ure 3 demonstrates the ability of the hypotheses to predict multiple futures with the aim of capturing different modes of the target data distribution. Tactis2, TimeGrad, TempFlow struggle produce predictions that deviate significantly from the mean, with most of their predictions being sampled from the same mode of the distribution. In- terestingly, T TimeMCL demonstrates the ability to generate predictions far from the mean with a non-negligible prob- ability, indicating that the model successfully captures dif- ferent modes and estimates the probability of each mode, including possibly rare events.
Table 3. Total Variation (L) comparison for K = 16 hypotheses.
See Table 11 for the full scores.
Tactlo2
TimeDrad
VongPlow
TinoMOL
7AmeMML (A.
Eusc.
154252
312445 1004
454597
128/75
245908
Excu.
0.237
0.606
1.104 -0100
0.631
0.042
SolaR
5383
3053
3341
TRAFFI
11,063
12:901
18.314 200117
5 1.662243
5.714 -0:3
TAXI
6452
4330 16
4147.2 1192
709 86 ex
701A61
2634547
2458503
12410009
putationally efficient model, just behind DeepAR, while achieving significantly better distortion scores. On that account, TimeMCL strikes a promising trade-off between computational cost and performance. For more details, refer to Appendices C.2.4 and C.2.5.
6.3. Effeet of the number of hypotheses
Table 4 presents a comparison of performance as a function of the number of hypotheses on the test split of the SOLAR dataset. See also Table 10 for an evaluation of the run-time as function of the number of hypotheses. As expected, the methods generally show improved performance as the number of hypotheses increases, although this also leads to longer run times. However, the performance improvement is not always strictly monotonic with TimeMCL. This sug- gests that our method still has room for refinement, as some hypotheses may remain slightly underutilized. We suspect that this may be partially due to the choice of scaling by the mean, which could be suboptimal for initializing i meMCI. Indeed, as mentioned above, we scale the data by dividing by the mean estimated dimension by dimension, which gen- erally results in a target of constant sign when time series values lie far from the origin, which may be suboptimal when hypotheses are randomly initialized around the origin. To address this, we plan investigate the impact of different scalers and data pre- prepprccsssing techniques, such as Z-score normalization or reversible instance normalization (Kim et al., 2021), which could promote better use of hypotheses.
Smoothness. Table 3 displays the Total Variation defined as TV = > ET || - 2/1122 which quantifies the average smoothness of the predicted trajectories (lower is more smooth). We see that TimeMCI, when trained either with the annealed or relaxed variant, provides significantly smoother trajectories compared to the baselines, further confirming the claim of Section 5.2 as a consequence of Proposition 5.1.
Comparing TimeMCL with the baselines on standard metrics. Table 7 and 8 provides performance on CRPS and RMSE, respectively. We see that, except for CRPS and RMSE on EXCHANGE (for which Tactis is better), TimeGrad outperforms the other baselines. We also ob- serve that TimeMCL is competitive despite optimizing a completely different training objective, at fraction of TimeGrad's and Tactis's computational cost.
Qualitative comparison. We qualitatively compare the predictions of TimeMCL with those of the baselines on the SOLAR, TRAFFIC, and ELECTRICITY datasets, as shown in Figure and in Figures 7 and the Appendix. Fig-
8

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Observations
TimeMCL mecm
TineNCL hypothess
TimoMCL forecasts
600
100000
0,903
25.0
100000
0.000
25.9
B31-USD
BTC USD
29000
808-1150 1S-110
STC-3SD
2631
0,000
0.77
600
30390
0.830
87.
500
9000
20.0
0.00
1230
12.99
1200
den
12m
12:39
8r0
12.0
dmg
12:0
mon
Figure 4. Cryptocurrency price forecasting for January 4, 2025 with K hypotheses. Setup described in Scction 6.4, Legend
the same as in Figure 3. We see that TimeMCL produces at least one hypothesis that closely matches the realized trajectory. We expect each hypothesis to specialize in distinct market scenario, for example, a trend reversal or sudden collapse.
both the smoothness of the TimeMCL trajectories and their capacity to capture rare events, Implementation details and further analysis can be found in Appendix C.4.
Table 4. Efffec of the number of hypotheses on the Distortion Risk for SOLAR. Results are averaged over four training seeds. TimeMCL (R.) and TimeMCL (A.) ) corresponds to the relaxed and annealed variants. See Table 9 for the full results.
7. Conclusion
K' Tactis2 TimeGrad TempFlow
TimaMCL (R.) )
TimeMCL (A.)
433.76 su 398.96 ...
422.24 mm
422.75 cm
422.75 aur
2
410.23
380.99
204.89
384,28
418.27
We introduced TimeMCL, a new model for time series fore- casting, designed to predict multiple plausible scenarios that can be viewed as an optimal quantization of the future distribution of the process we aim to predict. This model effectively captures different modes of distribution while providing smooth predictions. Moreover, unlike traditional models that rely primarily on the mean and variance of their predictions for interpretation, TimeMCL provides directly interpretable predictions by capturing multiple plausible future scenarios.
3
377.71 am 377.1 ax
398.96
356.42
358.48 s
4
376.91 1 375.25 1:221
386.28 am
326.62 2a0
347.45 aus
van 374.96 1:2
385.54
376.25
323.2
8
367.54 1:54 371.97
379.93
305.63
347.82
16
358.01
sue 360.6
371.14 2:
280.03 201
305.5 a.m
6.4. Financial time series
Ti MeM.. was evaluated a corpus of years of hourly cryptocurrency prices (15 correlated tickers collected from YahooF inance). Table 5 lists the assets, along with their pairwise correlations and price scales. Due to the wide variance in price magnitudes across assets, we apply Z- score normalization rather than mean scaling. All models are trained with K- = hypotheses. TimeMCL employs the annealed winner-takes-all loss and is compared with Tactis2 is,, TimeGrad, and TempFlow. As reported in Table 12, TimeMCL demonstrates strong performance in Distortion, RMSE, and CRPS, while also providing smooth trajectories at a reasonable computation cost (measured in FLOPs). Among the baselines, Tactiss s2 and TempF Low (with a transformer backbone) achieved competitive results, whereas TimeGrad, DeepAR, and the RNN version of TempFlow were less effective. Figure 4 shows represen- tative trajectories of our approach. At least one T MeMCC hypothesis aligns closely with the realized path, illustrating that individual hypotheses can specialize in distinct market regimes (e.g., trend reversals, sudden collapses). Figure 6 contrasts these predictions with the baselines, highlighting
immMCCL could be useful when combined with other back- bones, it complements state-of-the-an approaches with its objective function. In this paper, implemented the model using recurrent neural network to model temporal dependencies. Exploring the same approach with different architectures and datasets, such as transformer based back- bones when scaling the dataset size, would be a valuable direction for future research.
Limitations. We found that the choice of scaler can be limitation in TimeMCI. An inadequate scaler may bias the model toward certain hypotheses. While relaxation tech- niques help by softening hard winner assignments, we plan to explore advanced normalization methods to further im- prove T TimeMCL's 's vanilla setup. Another possible limita- tion is that TimeMCL requires the number of predictions to be predefined beforehand. Exploring dynamic rearrange- ments of hypotheses when adding new ones without full retraining is left for future work.
9

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Acknowledgements
Chung. J., Gulcehre, C., Cho, K., and Bengio, Y. Empirical evaluation of gated recurrent neural networks sequence modeling. arXiv preprint urXiv:1412.3555, 2014.
This work was funded by the French Association for Tech- nological Research (ANRT CIFRE contract 2022-1854) and the LISTEN Laboratory of Télécom Paris. also benefited from access to the HPCC resources of IDRIS (allocation 2024- AD01 1014345) by GENCI. We are grateful to the reviewers for their insightful comments.
Corlay, S. The nystr)," om method for functional quantiza- tion with an application to the fractional brownian motion. arXiv preprint 1009. 2010.
Dinh, L., Sohl-Dickstein, J., and Bengio, Density estima- tion using real nvp. In ICLR, 2017.
Impact Statement
Drouin, A., Marcotte, É., and Chapados, N. Tactis: Transfomeer-ttennionn copulas for time series. In ICML, 2022.
This paper presents work whose goal advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Du, Q., Faber, V., and Gunzburger, M. Centroidal voronoi tessellations: Applications and algorithms. SIAM Review, 1999.
References
Alexandrov, A., Benidis, K., Bohlke-Schneider, M.,
Garcia, N. C., Bargal, S. A., Ablavsky, V., Morerio, P., Murino, V., and Sclaroff, S. Distillation multiple choice learning for multimodal action recognition. In WACV, 2021.
Flunkert, V., Gasthaus, J., Januschowski, T., Maddix, D. C., Rangapuram, S., Salinas, D., Schulz, J., Stella, I.., Türkmen, A. C., and Wang. Gluonts: Probabilistic and neural time series modeling in python. JMLR, 2020.
Gasthaus, J., Benidis, K., Wang. Y., Rangapuram, S. S., Salinas, D., Flunkert, V., and Januschowski, T. Proba- bilistic forecasting with spline quantile function rnns. In AISTATS, 2019.
Arcozzi, N., Campanino, M., and Giambartolomei, G. The Karhunen-Locd Theorem. PhD thesis, Master's thesis, University of Bologna, 2015.
Arthur, D. and Vassilvitskii. S. k-means++: the advantages of careful seeding. In SODA. Society for Industrial and Applied Mathematics, 2007.
Gersho, A. and Gray, R. M. Vector Quantization and Sig- nal Compression. Springer Science & Business Media, Boston, MA, 1992.
Ashok, A., Marcotte, É.,, Zantedeschi, V., Chapados, N., and Drouin, A. TACTiS-2: Better, Faster, Simpler Attentional Copulas for Multivariate Time Series. In ICLR, 2024.
Gneiting, T. and Katzfuss, M. Probabilistic forecasting. Annual Review of Statistics and Its Application, 2014.
Godahewa, R. W., Bergmeir, C., Webb, G. I., Hyndman, R., and Montero-Manso, P. Monash time series forecasting archive. In NeurlPS, 2021.
Asuncion, A., Newman, D., et al. Uci machine learning repository, 2007.
Benidis, K., Rangapuram, S. S., Flunkert, V., Gasthaus, J., Wang, Y., Bohlke-Schneider, M., Salinas, D., Stella, L.. Callot, L., and Januschowski, T. Neural forecast- ing: Introduction and literature overview. arXiv preprint arXiv:2004. 10240, 2020.
Graves, A. Generating sequences with recurrent neural networks. arXiv preprint arXiv:1308.088. 2013.
Guzman-Rivera A., Batra, D., and Kohli, P. Multiple choice learning: Learning to produce multiple structured outputs. In NeurlPS, 2012.
Bennett, W. R. Spectra of quantized signals. The Bell System Technical Journal, 1948.
Hendrycks, D. and Gimpel, K. Gaussian error linear units (gelus). arXiv preprint arXiv: (tl::6666.0... 2016.
Bosq, D. Linear Processes Function Spaces: Theory and Applications. Springer, 2000.
Ho, J., Jain, A., and Abbeel, P. Denoising diffusion proba- bilistic models. In NeurIPS, 2020.
Chen, and Boccelli, D. L. Real-time forecasting and visualization toolkit for melti-seasonal time series. Envi- ronmental Modelling & Software, 2018.
Hochreiter, S. and Schmidhuber, J. Long short-term memory. Neural computation, 1997.
Chevallier, J. Uniform decomposition of probability mea- sures: quantization, clustering and rate of convergence. Journal of Applied Probability, 2018.
Huang, C.-W., Krueger, D.,, Lacoste, A., and Courville, A.
Neural autoregressive flows, In ICML, 2018
10

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Hyndman, R., Koehler, A. B., Ord, J. K., and Snyder, R. D. Forecasting with exponential smoothing: the state space approach. Springer Science & Business Media, 2008.
Pages, G. and Printems, J. Optimal quantization for finance: from random vectors to stochastic processes. In Hand- book of numerical analysis. Elsevier, 2009.
Hyndman, R. J. and Khandakar, Y. Automatic time se- ries forecasting: the forecast package for Journal of statistical software, 2008.
Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., and Lakshminarayanan, B. Normalizing flows for probabilistic modeling and inference. JMLR, 2021.
Perera, D., Letzelter, V., Mariotte, T., Cortés, A., Chen, M.,
Karatzas, I., Shreve, S. E., Karatzas, I., and Shreve, S. E. Brownian motion. Brownian motion and stochustic cal- culus, 1998.
Essid, S., and Richard, G. Annealed multiple choice learn- ing: Overcoming limitations of winner-takes-all with nealing. In NeurIPS, 2024.
Kim, T., Kim, Tae, Y., Park, C.,, Choi, J.-H,, and Choo, Reversible instance normalization for accurate time-series forecasting against distribution shift. In ICLR, 2021.
Ramsay, O. and Silverman, B. W. Functional Data Anal- ysis. Springer, 2005.
Rangapuram, S. S., Seeger, M. W., Gasthaus, J., Stella, L., Wang, Y., and Januschowski, T. Deep state space models for time series forecasting. In NeurIPS, 2018.
Kong, 7 Ping, W., Huang, J., Zhao, K., and Catanzaro, B. Diffwave: A versatile diffusion model for audio synthesis. In ICLR, 2021.
Rasul, K., Seward, C.. Schuster, I., and Vollgraf, R. Au- toregressive denoising diffusion models for multivariate probabilistic time series forecasting. In ICML, 2021a.
Lai, G., Chang, W.-C., Yang, Y., and Liu, H. Modeling long- and short-term temporal patterns with deep neural networks. SIGIR, 2018.
Rasul, K., Sheikh, A.-S., Schuster, I., Bergmann, U. M., and Vollgraf, R. Multivariate probabilistic time series forecasting via conditioned normalizing flows. In ICLR, 2021b.
Lee, S., Purushwalkam Shiva Prakash, S., Cogswell, M., Ranjan, V., Crandall, D., and Batra, D. Stochastic multi- ple choice learning for training diverse deep ensembles. In NeurlPS, 2016.
Rose, K., Gurewitz, E., and Fox, G. deterministic anneal- ing approach to clustering. Pattern Recognition Letters, 1990.
Lee, S., Purushwalkam, S., Centeno, M., Zaidan, A., and Batra, D. Confident multiple choice learning. In ICML, 2017.
Rupprecht, C., Laina, I., DiPietro, R., Baust, M., Tombari, F., Navab, N., and Hager, G. D. Learning in an uncertain world: Representing ambiguity through multiple hypothe- ses. In ICCV, 2017.
Letzelter, V., Fontaine, M., Chen, M., Pérez, P., Essid, S., and Richard, G. Resilient multiple choice learning: A learned scoring scheme with application audio scene analysis. In NeurlPS, 2023.
Salinas, D., BohlBB-hhhhheeeeeee Schnke-Schheider, M., Callot, L.,, Medico, R., and Gasthaus, High-dimensional multivariate forecast- ing with low-rank gaussian copula processes. In NeurIPS, 2019.
Letzelter, V., Perera, D., Rommel, C., Fontaine, M., Essid, S., Richard, G., and Perez, P. Winner-takess learners are geae conditional density estimators. In ICML, 2024.
Salinas, D.,, Flunkert, V., Gasthaus, J., and Januschowski, T. Deepar: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 2020.
Lloyd, S. Least squares quantization pem. IEEE Trans- actions on Information Theory, 1982.
Loubes, J.-M. and Pelletier, B. Prediction by quantization of conditional distribution. Electronic Journal of Statistics, 2017.
Seo, Y., Lee, K., Clavera Gilaberte, I., Kurutach, T., Shin, J., and Abbeel, Trajectory-wise multiple choice learning for dynamics generalization in reinforcement learning. In NeurlPS, 2020.
Matheson, J. E. and Winkler, R. L. Scoring rules for con- tinuous probability distributions. Management science, 1976.
Sutskever, I., Vinyals, O., and Le, Q. V. Sequence to se- quence learning with neural networks. In NeurIPS, 2014.
Pagès, G. Introduction to vector quantization and its appli- cations for numerics. ESAIM: Proceedings and Surveys, 2015.
Tian, K., Xu, Y., Zhou, S., and Guan, J. Versatile multiple choice learning and its application to vision computing. In CVPR, 2019.

---

Winner-takes-all for Multivariate Probabilistic Time Scries Forecasting
Torres, J. F., Hadjout, D., Sebaa, A., Martinez-Alvarez, F., and Troncoso, A. Deep learning for time series forecast- ing: survey, Big Data, 2021.
Tsay, R. S. Analysis of financial time series. John Eiley and Sons, 2005.
Van Den Oord, A., Dieleman, S., Zen, H., Simonyan, K., Vinyals, O., Graves, A., Kalchbrenner, N., Senior, A., Kavukeuoglu, K., et al. Wavenet: A generative model for raw audio. arXiv preprint 1609.03399 2016.
Zador, P. L. Asymptotic quantization error of continuous sig- nals and the quantization dimension. IEEE Transactions on Information Theory, 1982.
12

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Organization of the Appendix
The Appendix is organized as follows. Appendix A contains the proofs of the theoretical results, establishing that TimeMCL can be interpreted functional quantizer. Appendix B describes the synthetic data used in the toy example. Appendix details the empirical evaluation on real data, covering the datasets (Appendix C.I), the experimental procedure and evaluation metrics (Appendix C.2), and the baseline models (Appendix C.3). Finally, Appendix C.4 presents additional results for the financial time series.
A. Theoretical results
Following the main paper notations, we assume each time series lives in X' R.P. We refer the context and target sequences 71:to-1 x1:to-1 and Tta:T E X respectively.
A.1. Proof of Proposition 5.1
Let us consider the following assumptions. Assumption A.1 (True risk minimization). The batch size is big enough so that the difference between the CWTA risk and its empirical version can be neglected.
In the Assumption A.2 that follows, we considered the forecaster (or unrolled) network Fg function X1:6--1 + 2 that is directly optimized during training, ie., teacher foreing is disabled. This implies that the same forward computation applies in training and inference. In particular, the model predictions on the target window xT training time depend only on the past context -1: We adopt this simplified setup in the toy experiments presented in Section 5.3 and in Appendix B.
Assumption A.2 (Expressiveness). The family of neural networks considered is expressive enough so that minimizing the expected training risk reduces to minimizing the input-dependent risk for each individual input'. Formally, this means
inf Ex., (( F.,T1 = -S 1Ee
inf Beo [ [e,,, [6,r)) dp( Ap(#1110-1) 1
(11)
where è is the set of possible neural network parameters, and l may be the mean square error loss & (Zto:T, Ztoer) A
EL 1 zzll2.
Assumption A.3 (Optimality). The training has converged and CWTA has reached a local minima.
Proposition A.4. Under the Assumptions A. A.2 and A.3, TimeMCL is a conditional stationary quantizer for each sampled window 11 Tt. that is
F (1) E[ FtoiT (x(
(12)
where
Xo(IIto-1) = Tto:T E XtoiT (F < (T Vr k].
1, lto:T) - | (2 zell². This makes TimeMCL akin where accordingly with Assumption A.2; (x 1,
to conditional and gradient- bassed version of K-Means over the set of possible future trajectories.
futuur
Proof. Assuming we are performing true risk minimization (Assumption A.1), the WTA Loss LTA can be re-written as
K
CWTA
/
c) (2140-1, Rog:) dp(x6, F1: -ddo-1)pp((16--.))
(13)
Beger E. (1
due Chasles relation because the Voronoi tesselation forms partition of the output space. Under the expressiveness Assumption A.2, (13) comes down to optimizing
K
Farag-1 (z)
( 11ito-1),
(14)
Xe(miao 1)
k=1
"See also Assumption in Perera et al. (2024) and Section Loubes & Pelletier (2017).
13

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
for each fixed context -1 where each z E X and
T (( xtt)
1=
To prove that each model (11:t-11) converges to the conditional mean of the future trajectories that fall into its corresponding Voronoi cell, follow the decoupling strategy of Rupprecht et al. (2017) and Letzelter et al. (2024), Define two sets of variables:
Generators {g ( which induce a partition of the output space Tt&T via the Voronoi diagram:
Xx(g) (toor X | 9 (2110o < |
Vr + k).
Centroids {= (2112-11) which are the points used to compute the intra-cell T.? loss.
Decoupling generators vs. centroids through alternating optimization. In T TimeMCL. the update is done in a two-step fashion: first, the winners are computed for each point in the batch, and then the latter updated. This two-step optimization allows us bypass the non-differentiabiiity of the min operator in (4).
Let us decouple the two variables and define the following functional to optimize:
K Σ
Flg, 2;
|[ 2 (1116-1)11 p( 11:10-1) dz
Nx(a)
R=
We now consider minimizing F(g, 2; 11:tt-1) with respect g and in alternating manner:
(a) Fixing the partition (g) and optimizing the centroids (z). In this case, the gradient update of 2 has direction:
OF
((=11)-))Vol(((()))) - -Sa Tta:TP((Ito.T [ (ddee))))))))))
O 2z (g,
xx((s)
where (o)( - Sxa(0e p( 11 Thus each z (1) is updated in the direction of the cell conditional mean, similarly to the Lloyd algorithm (Lloyd, 1982).
(b) Fixing the centroids (2) and optimizing the partition (g). Conversely, if is fixed, one can reduce F 2; P1:to -1) by ensuring that X (g) is indeed the Voronoi cell generated by ~ (x Indeed, for any deviation from strict Voronoi partition, there would exist subsets of Tt. incorrectly assigned to some cell, and reassigning them the nearest 2 (1) would lower the overall loss (See also Rupprecht et al. (2017), Theorem 1).
By alternating between these two steps, any stationary point of F(g, 20 21:1001)) must satisfy
2 (I1tto-) -=((0100-1))
Such configuration is centroidal Voronoi tessellation (see, e.g., Du et al., 1999).
Reverting the original notation of the WTA network, identify g (x 1) 2 =* (x11to))) optimality, thus each prediction head coincides with the mean of the distribution restricted to its Voronoi cell:
F (T1:t0-1) E Wooel
Fto: (x(11tt-11]-
Consequently, at the global minimum the WTA objective, we obtain the desired conditional centroidal Voronoi tessellation.
14

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
A.2. Proof of Proposition 5.2 (Scoring Heads Unbiased Estimators)
Once trained. predictions from score heads are performed through the unrolled version of the score heads 1 which can be viewed as function 1 x [0. 1] when averaging the predicted scores the sequence. In the following, assume that the r is directly optimized during training, similarly to the approach used for the prediction heads in Section A. We prove that when the training of T' TimeMCL with / > 0 has converged globally, the (unrolled) scoring heads r (2 match the conditional probability mass of their respective Voronoi region X (2)::o1)). This statement generalizes the arguments from Letzelter et al. (2024) to our time-series setup.
Assumption A.5 (Global optimality for both centroids and scores), In addition having converged to global minimum of the centroid objective (13), the TimeMCL model also reaches global minimum for its scoring objective
(0)
BCE [[ (2 r (Z1100 p( -10 -1,7)) drig d
where BCE(.,)) denotes the binary cross-entropy, and F1:t-11 H 11) is the optimal Voronoi tessellation corresponding to the converged heads {
Proposition A.6 (Unbiased Estimator of Voronoi Cell Mass). We assume, as in Assumption A.I, that the batch size big
enough so that the difference between the L and empirical version can be neglected. Under Assumption A.5, and assuming perfect expressiveness of the score heads in Assumption A.2, for any context Ti:to-1 and index (1,.... K}, the optimal scoring head satisfies
D (3 PP (Ntg:T Agg (1)1:)-1) (1:t6-
Proof. We adapt the derivation of the unbiased property our time-series setup. Recall that the optimal scoring objective (see Assumption A.5) may be written
min
E/ X.(21n)-))
BCE [ Ng (2 1,, r (2 (2150-1) p p( 21 p( dr dæ 1sto -1
(15)
Under the true risk minimization and expressiveness assumptions, and since the input-dependent risk decomposes over each E {I....,Kt K}. we can focus on a fixed T1:t-11 and single index the sum. For that fixed #1s0-1 and kt the part of (15) to optimize is:
BCE
1 E Xe (e I)l, 1 (F1tt-1)) p(acor [1:to dx
Writing out the binary cross cross-entropy explicitly, we get:
- log (1 p( xisto- :o f Ix Vi (miin 1)
log Log(1-11)) (Zi ( p( x1sto-1) do
(x 1)
Let us denote mi ( = P E (x( x1) Then the above expression becomes
m4 (m log( (T) (T - (1 my (x(10-111) log (1 1% (x1:
We recognize this the binary cross-entropy between the probabilities 12ys (1 and my (I The unique global minimum of that scalar binary cross-entropy is attained at
(1 (= P( E Xg (#1 F1
This equality must hold for each context 01:to-1 and each index k, completing the proof.
B. Experimental details on the synthetic datasets
All toy examples were conducted using neural network with a three-layer fully connected backbone (200 hidden units, ReLU activation) and K prediction heads, each generating a trajectory of given length. Training uses a Winner-Takes-AN loss combined with squared error minimization, optimized with Adam using a learning rate of 10 33 a batch size 4096, and 500 iterations. For each of the setups, we used the Relaxed version of TimeMCL with & 0.05.
15

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Brownian Motion. In order to conduct interpretable experiments and generate synthetic data, we employed Brownian motion. A Brownian motion {W: 1220 is a Gaussian process characterized by mean E E[W.] 0 and covariance function Co, Wa) min(t, s) for t, 2 0. This process possesses several properties of interest (Karatzas et al., 1998). In particular, the Markov property, which ensures that the process {W - W, J220 is itself a Brownian motion starting at O, independent of the process (W, 0 I & S 8). This implies that to condition a Brownian motion on its trajectory up to time I, suffices to condition on its value Wt- The process beyond time t behaves as if were a new Brownian motion originating at W, In conditional quantization, we aim to quantize Brownian motion at a given prediction horizon, conditioned on its value time f. To compute the target quantization, we followed the methodology detailed in (Pages & Printems, 2009), leveraging the Karhunen-Loève decomposition (Arcozzi et al., 2015) of the Brownian motion on & E [0, 1]:
We Σ V AW (t),
(16)
where En ~ 10, are independent standard Gaussian random variables. Here, AW 1
=²(n-})² are the eigenvalues
of the covariance operator of Brownian motion, and en (t) v2sin (n - (1) are the corresponding eigenfunctions (Corlay, 2010), To quantize Brownian motion, we first truncate its Karhunen-Loève decomposition (16) to terms. Then, the coefficients En are quantized into K, optimal levels {ai {a(Kk)
Kx, (K=)), using the quantile function of the normal distribution. All possible combinations of quantized coefficients are generated via the Cartesian product. In this case, for each multi-index (i1:-..; im) E IIm [1, 1, it is possible to define quantizer
(K1 Kmn 33
(1)
A
(Kn) W Oin en (1)
which takes [I=-1 Kr, possible values. Following Pages & Printems (2009), product quantization with exactly K codevectors requires optimizing two parameters: the truncation parameter m and the quantization level allocation {K, J==1, with IT.= Kn K. Pages & Printems (2009) provides a list of optimal tuples for these parameters. our numerical example, we select 2 and (K1, K2) (5 resulting in quantizer with x 2 10 codevectors.
TimeMCL was trained 250-step sequences randomly drawn from Brownian motion on [0, 1]. discretized with 500 steps. At each iteration, a new sample was generated, with the model conditioned on the last observed value and tasked with predicting the next 249 steps. The generated trajectories are smooth and closely match the theoretical optimal ones. The model demonstrates consistent long-term predictions and effectively learns conditional quantization.
Brownian Bridge. Brownian motion is unrepresentative of most natural processes, as its future evolution is independent of the conditioning time. In contrast, many real-world processes depend on when they are observed. A Brownian bridge, viewed as a Brownian motion conditioned to reach a fixed value at T = 1, better captures such dependencies. A continuous process (Ejoo) a Brownian bridge on the interval [0, T] if and only if it has the same distribution as {W, / \eeoor where [Wereejo,r; is a standard Brownian motion. In our case, T' 1, and the process B is referred to standard Brownian bridge. As need to conditionally quantize the Brownian bridge, note that if we pick random time [o, 1], the Brownian bridge conditioned on its value at this point remains Brownian bridge. To obtain the optimal quantization, the same process as for Brownian motion is used, with adapted eigenvalues and eigenvectors in the Karhunen-Loove decomposition (See Corlay (2010), page 3),
2 sin TTT)
T A
(1)
V
T2 n n21.
In our numerical example, we select the same values for m and (K1, K2) as those chosen for the Brownian motion. TimeMCL, like for Brownian motion, is trained on random samples from a Brownian bridge starting at Bg = 1 and ending at B, = 1. With 500 discretization points, predicts 250 steps, conditioned on the prediction start time and value.
Autoregressive model AR(p). The two first experiments involved Brownian motion, where only the last value matters for prediction, and the Brownian bridge, which depends on the last observed time and value. While these processes assess meMCL's quantization ability, they lack strong context dependence. Real-world applications require modeling context to
16

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
predict and quantize future trajectories effectively. Consider an autoregressive process
X, o +
(17)
where (ez), is a white noise with variance o². We generate an autoregressive sequence where the first five values follow a normal distribution with scale a. The process runs for an initial warm-up period of 100 steps. After this, we predict and quantize the next 250 steps using the last 100 observations, i.e., {Xto -999, -9.....; Xt,) ). The five values preceding in form crucial context- -iithout it, correctly quantizing the future distribution is impossible. TimeMCL is expected to exploit this information for optimal quantization. To the best of our knowledge, no simple analytical solution exists for quantizing an autoregressive process. To approximate an optimal quantization, the process is simulated multiple times from to. conditioned on its context, to construct a conditional distribution. The Lloyd algorithm, run with a very large number simulations (105 sampled trajectories), produces a set quantized trajectories that serve as references for TimeMCL.
TimeMCL is trained on 250-step sequences randomly sampled from 500-step trajectories of the AR(S) process, with a new batch generated at each iteration. To ensure consistency with Brownian bridge and motion experiments, time is normalized by dividing all steps by the total number of steps. We took Φ1 = 0.4, d2 d3 0.2, 04 Φ5 = 0.1, and ø 0.06. The trajectories predicted by TimeMCL closely align with the optimal ones found by Lloyd's algorithm, exhibiting smooth behavior. Various values and orders of the autoregressive process, both stationary and non-stationary, were tested, consistently yielding results aligned with Lloyd's algorithm.
C. Experiments with real data
C.1. Datasets
We evaluate approach six well-established benehmark datasets taken from Gluonts (Alexandrov et al., 2020), each containing strictly positive and bounded real-valued series. See Lai et al. (2018) for an extensive description.
The SOLAR dataset (Lai al., 2018) consists of hourly aggregated solar power outputs from 137 photovoltaic sites, spanning 7009 time steps. Strong daily seasonality is typically observed due to the day-night cycle.
The ELECTRICITY dataset (Asuncion et al., 2007) contains hourly power consumption data from 370 clients across 5833 time steps. Demand patterns often exhibit both daily and weekly periodicities. driven by regular human activity and business operations.
The EXCHANGE dataset (Lai et al., 2018) features daily exchange rates for eight different currencies, with 6071 observations per series. Unlike energy or traffic data, exchange rates often lack clear seasonal patterns and are influenced by broader macroeconomic factors.
The TRAFFIC dataset (Asuncion et al., 2007) comprises occupancy rates (ranging from 0 to 1) measured hourly by 963 road sensors over 4001 time steps. These series generally display recurrent hour peaks as well as differences between weekdays and weekends.
The TAXI dataset consists of traffic time-series data of New York City taxi rides, recorded at 1214 locations every 30 minutes in January 2015 (training set) and January 2016 (test set). We used the preprocessed version from (Salinas et al., 2019).
The WIKIPEDIA dataset (Gasthaus et al., 2019) contains daily views of 2000 Wikipedia pages.
Table 5 provides an overview of the main characteristics of these six datasets. Note for each dataset, we used the official train/test split. We dedicate 10 times the number of prediction steps for validation, at the end of the training data.
C.2. Metries
In all the following, each model produces K hypotheses {&t T Additionally Time-MCL predicts K scores (one per trajectory), that we denote as ......, .......A. 2K E [0, 1]. with E 2 = 1, omitting the input T1:t--1 for conciseness. The observed trajectory denoted as : E Dx We used the Gluonts integrated watour to compute the
17

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Table 5. Datasets characteristics.
DATASET
DIMENSION D DOMAIN X
FREQ.
TIME STEPS
PRED STEPS
t0 +1
SOLAR
137
IR+
HOUR
7,009
24
ELECTRICITY
370
R
HOUR
5,833
24
EXCHANGE
IR+
DAY
6,071
30
TRAFFIC
963
(0, 1)
HOUR
4,001
24
TAXI
1,214
N
30-MIN
1,488
24
WIKIPEDIA
2,000
N
DAY
792
30
metrics, and we customized it to include the distortion. For evaluation over the whole test set, we have N input-targer pairs (11 r). E N... We denote by s. the k-th hypothesis generated for input #1 at time L.
In the following, we define the summation over dimension operator as A: E ERP. PP= Σd Id R. where x denotes the d-th dimension of the vector x.
C.2.1. RMSE (RoOT MEAN SQUARE ERROR)
Let us denote by & (x the conditional mean estimator of p( (1:t--1 given the probabilistic model evaluating. For TimeMCL, used z Σ n f..r and for the other methods that do not score heads,
F K Ek= T. RMSE is defined as:
RMSE
1 T to
(A(zt) - A(
(18)
In Gluonts, we can access this metric under the name rmse. The prefix m.sum indicates in (18) that the RMSE is computed after aggregating the target and the prediction dimension by dimension.
C.2.2. CRPS (CONTINUOUS RANKED PROBABILITY SCORE)
In the following, we denote by Qx (q) the quantile of order a real random variable X, with E [0, 1] (when it exists).
Let us define the random variable X which takes the value A(alii) with probability 2n. or uniformly with probability K when using no score heads, respectively. Let us introduce:
Lylzio.r (b)
Σ A(x:) Qx
xx (xx(a) Qx;(a)) -0). (
A(x)),
where L, and T are referred as the quant leLoss and abs.target.sum in Gluonts.
CRPS (Matheson & Winkler, 1976) is computed as
L, ( (ho.))
CRPS A
i=
(19)
O
I (
GEQ
where we used @ {0.05, 0.1,... 0.95J. Equation (19) is referred as sum mean wQuantileloss in the python library, CRPS-Sum (Salinas et al., 2019; Ashok et al., 2024) in the related literature.
C.2.3. DISTORSION
Distortion, also referred to as the Oracle or Quantization error in the literature (Pages & Printems, 2009; Lee et al., 2016; Perera et al., 2024) was defined as
D
T
D2
min Ky
Σ
j.d
T to
+
d=1
This metric was integrated into the Eval uator class from Gluonts.
18

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Table Distortion comparison with Hypotheses. Results are averaged over four training seeds. ETS, Trf. TempFlow and Tactis2, columns are in gray because they don't share the same backbone as the other baselines. Best scores are in bold, and second best are underlined.
ETS
TrE. Temmpllo Pow
Factis2
TimeGrad
DeepAR
TempPllow
TimeMCL (R.)
TimeMCL (A.
ELEC.
15132 =: 1662
14160 = 3370
12557 :: 2405
10172 1: 680
144832
3081
15822 :: 531
12155 414
11509 :: 703
EXCH.
0,057 1: 0.002
0,066 0.012
0.031 0.002
0.006
0.062 0.001
0.051 0.007
0.039 0.001
0.042 0.005
SOLAR
654,48 + 3.45
387.15 18:57
367.54 35.4
371.97
33.31
770.48 + 19.91
379.93 22.38
305.63 +: 28.24
347.82 28.82
TRAFFIC
2.65 +: 0.0
1.28 + 0.02
0.85 + 0.04
0.79 + 0.05
2.14 + 0.06
1.23 + 0.04
0.69 + 0.02
0.71 + 0.01
TAXI
589.15 = 0.75
284.25 = 13.94
245.09 + 8.79
211.71 = 4.02
429.01 = 14.88
272.84 = 12.14
187.04 2.45
256. 16 + 24.95
WIKI
722166 = 5248
528546
11657
255328
3386
261463 + 1047
371035
6578
388511 + 2472
271878 = 17358
258468 + 3108
C.2.4. FLOPS (FLOATING POINT OPERATIONS)
FLOPs serve as an approximate measure of the model's computational load. For each baseline, we only considered FLOPs evaluation at inference, with a single batch of size on the EXCHANGE dataset. We used the F FlopCountAna lysis function from the fvcore library measure it.
C.2.5. INFERENCE RUN-TIME
Inference time was computed on single NVIDIA GeForee RTX 2080 Ti, while making sure it is the only process that runs on the node. It was measured over the whole test set of EXCHANGE (make_evaluse ion.prediiceco ct ions function in Gl uontss, using a batch size of 64. For fair and accurate time comparison with respect to the number of hypotheses, we disable parallel sampling for each baseline (by setting num.para wwwraaleeelsaappl to 1).
C.2.6. TOTAL VARIATION
We quantify the average smoothness of the sampled trajectories using the Total Variation (TV). Given K sampled trajectories. we defined it as
N N
K
24.1 1/.111)
TV
F241
(20)
K
k=1
For Ti Teme,,L, we sampled the hypotheses in proportion to their scores for Total Variation computation instead of uniformly as in (20) for the other baselines. The smoother the predictions of probabilistic estimator, the smaller its total variation.
C.3. Experimental details
C.3.1. ARCHITECTURES, TRAINING, AND INFERENCE
Recurrent Neural Network Backbone. Time-MCI TempFlow, TimeGrad, and DeepAR share the same RNN backbone, which takes as input a concatenation of feature types. Let be the batch size, T the time series length, and Durget the target dimension. The input includes (i) lagged target values, normalized and shaped as B T (Tiags X Durget) and (ii) Fourier features of shape B T (2 x Ny)- The combined input the RNN thus has shape B x T x (TTigs X Durgat) + (2 x N)). Our RNN implementation follows Rasul et al. (2021a) and consists of L 2 layers of Long Short-Term Memory (LSTM) cells, each containing C = 40 hidden units. The term Ny. which is associated with Fourier transforms, will be further detailed in the following sections.
Sampling and dataloader The time series is first segmented into past and future windows based on predefined lengths and lead time sampling slicing points randomly during training and using the last timestamp during predic- After transformation, the dataloader assembles inputs used across all models. It separates past and future targets and adds several structured features: (i) Target Dimension Indicator (B x Dearee) uniquely labels each (fi) an Observed Values Indicator replaces missing target entries with dummies and flags them as observed (1) or missing (0); (iii) Padding Mask marks synthetic padding to ensure uniform sequence lengths, preventing them from influencing training. Additionally, temporal features of shape B x Tped x 2Ny encode periodic patterns using sine and cosine projections, where N₂ denotes the number of unique values for a given frequency (e.g.. N₁ = 12 for months).
tion (Inst anttttaceesseee
target dimension;
Transformer Backbone, For the transformer-based version of TempFlow (Trf. TempFllow and for Tact is, we
19

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Table 7. RMSE (L) comparison for hypotheses.
Tvf. TempFlow
Tactis2
TimeGrad
DeepAR
TempFlow
TimeMCL (R.)
PimeMCL. (A.)
ELEC.
25771 +: 1117
20559 4: 11011
19342 2 5196
11597 1: 931
145977 = 3561
23573 2: 2614
18571 +: 3063
19151 1: 2657
ExCH.
0.089 0.01
0,147 0.043
0.062 0.007
0.077 D.22
0,086 0.002
0.021
0:084 +: 0.002
0,103 0.021
SOLAR
3661.86
4: 28.87
1121.74::100 213.06
3531,02 439.59
3244.92 1: 247.16
7199.09 +: 107.71
3801.72 4: 122.09
3639.23 606.34
3896.18 4: 168.71
TRANFIC
15:24 0.09
27.57
11.54
4.44 1.07
4: 1,48
26.77 1.29
5,48 +: 0.71
5.73 4
TAXI
9779.81
4: 13:42
4223.88 =: 561.11
2662.52 2: 58.79
2357,85
1: 133.05
7323.32 :: 215.35
4136.39
1: 581.06
6855
1: 3456.01
8519.15 +: 2691.98
+: 39186
1152828
171868
564176
115598
505869
60432
864607
105223
1210500 412274
1264339
166347
Table 8. CRPS-Sum (L) comparison for K = 16 hypotheses. Results averaged over four training seeds.
ETS
Tef. TempFlow
Tactis2
TimeGrad
DeepAR
TempFlow
TimeMCL. (R.) )
TimeMCL (A.)
ELBC.
0.0601 +: 0.0045
0.0434 :: 0.0214
0,0409 + 0.0105
0.0235 = 0.0018
0.4234 0.0101
0.0503 + 0.0038
0,0488 = 0.0074
0.0481 0.0086
EXCH.
0.0079 1: 0.001
0.0138 1 0.0039
0.0066 1: 0.0014
0.0074 0.0029
0.0081 0.0002
0.0091 1 0.0019
0.0103 1: 0.0002
0.0131 0.0031
SOLAR
0.509 = 0.0041
0.4528 : 0.0271
0,3919 = 0.0378
0.3789 = 0.0356
1.028 # 0.0225
0.4352 = 0.0215
0,3831 + 0.0658
0,3966 = 0.0079
TRAFFIC
0.2007 0.0014
0.394 0.0092
0.1339 +: 0.0166
0.0561
0.0137
0.5257
0.0205
0.3833 + 0.0271
0.067 0.0102
0.0651 + 0.009
TAXI
0.8948 1. 0.0017
0.288 0.0467
0.1756 1 0.0056
0.1449
0.0118
0.5223
0.0239
0.3106 1 0.0512
0,4438 1. 0.1965
0.5803 0,2206
WIKI
0.0789 1 0.0041
0,1309 + 0.0316
0.0649 = 0.019
0.0567 0.0084
0.6588
+ 0.2225
0.0905 + 0.0051
0.142 + 0.0624
0,1479 = 0.0392
employed the same architectures as in their respective original implementations (Rasul et al., 2021b; Ashok et al., 2024).
Optimization Training is conducted using the Adam optimizer with a learning rate of 10-3, and following an 'LR on plateau' scheduler (See Pytorch documentation), a weight decay of 10 and 200 training epochs. Additionally, a separate validation split is used with temporal size equal to 10 times the prediction length. When training TimeMCL, the WTA loss (and its variants) was divided by the prediction length T lo + 1. All models are trained using gradient norm clipping with threshold of 10. Training and evaluation are performed in single-precision (float32) with Py PyTorch.
Inference. We used the official experimental protocol for evaluation in this benchmark (e.g.. (Rasul et al., 2021a)). The official test dataset is divided into multiple non-overlapping subsets, each containing sufficient points for both context and prediction lengths, allowing comprehensive assessment across different temporal segments. To compute the TimeMCL metrics while respecting the probabilities given by the scores, we applied a simple trick: after a single forward pass, we obtain K hypotheses along with their associated probabilities (derived from the scores). To ensure that the hypotheses contribute to the metrics in proportion to their assigned probabilities, we performed resampling with replacement from these K hypotheses based on their respective probabilities before computing metrics.
C.3.2. BASELINES
DeepAR. (Salinas et al., 2020). The output distribution as a multivariate normal distribution with mean vector Ja E IR&,, a low-rank covariance factor L E Dxr and a diagonal covariance vector a E RD, where the covariance matrix Σ is defined as Σ = = + diag(o), (o). thereby ensuring efficient parameterization and capturing primary dependencies among the target dimensions. In all the experiments, we set the rank to
TimeGrad. (Rasul et al., 2021a). For the diffusion process, as in Rasul al. (2021a) we employ 100 diffusion steps with linear variance schedule, transitioning from 3 = 1 10-4 to 8100 = 0.1. The loss function utilized is the L loss. The residual architecture is integrated within the diffusion component to facilitate effective training and capture complex temporal dependencies. This residual network consists of 8 residual layers, each with 8 residual channels. The dilation cycle length is set to 2, which controls the expansion of the receptive field without increasing the number of parameters.
TempFlow Rasul et al. (2021b). TempFlow applies a conditional normalizing flow (Papamakarios et al., 2021) that is invertible and has tractable Jacobian, yielding an exact likelihood. Batch normalization after each coupling layer and per-series mean scaling further stabilizes optimization. We adopt the original TempFlor configuration from Rasul et al. (2021b): a two-layer LSTM conditioner with 40 cells and dropout 0.1, fed with a 100-step context, and a Real-valuee Non- VoluVolume Preserving (RealNVP) (Dinh et al., 2017) head composed of 3 coupling blocks whose scale/shift networks have 2 hidden layers with 100 hidden units. The transformers version of TempFlow Trf. . TempF1 1 low uses 8 heads with 3 encoder and 3 decoding layers, with an embedding dimension of 32, and a feedforward dimension of 128 with a dropout rate of 0.1 and Gelu (Hendrycks & Gimpel, 2016) activation layers.
20

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Table Distorsion (4) Comparison for Solar Dataset. Results averaged over four training seeds, this table, the distortion computed with a variable number of hypotheses K for each baseline, as in Table 4 of the main paper.
ETS
Trf. Tempp oow
Tactis2
TimeGrad
DeepAR
TempFlow
TimeMCI (R.)
TimeMCI (A.
1
685.77 +: 23.08
468.84 32.54
433.76 +: 19.02
398,96 + 32.38
838.14 + 26.13
422.24 i: 10.64
422.75 + 84.28
422.75 +: 84.28
2
678.49 + 13.29
424.7 + 14.47
410.23 + 12,98
380.99 + 31.06
817.75 15.28
404.89
20.53
384,28 + 19.92
418.27 + 95.95
3
678.94 10.6
409.97 20.37
377.71 +: 34.55
377.1 30.87
794.45
19.39
356.42 + 44.11
358.48 + 20.29
671.97 +: 16.52
402.49 20.51
376.91 +: 34.72
375.25 4: 32.23
789.66 18.97
18.5 386.28
326.62 +: 20.89
347.45 + 14.59
666.44
400.63 +: 18.4
374.76 32.31
374.96 32.36
776.94
17.93
385.54
19.14
376.25 +: 43.23
323.2 18.3
654.48
387.15 18.57
367.54 35.4
371.97
33.31
770.48
19.91
379.93
22.38
305.63 28.24
347.82
28.82
16
641.32 3:57
374.69 20.97
358.01 +: 38.47
360.6 34.79
748.68 18.92
371.14 = 18.2
280.03 + 15.22
305.5 = 4.05
Tactis2 (Ashok et al., 2024). Unlike the RNN-based TimeGrad, TempFlow and Time- Tact i encodes each time step as token that includes the value, binary missingness flag, static covariates, and positional code. Two Transformer encoders process these tokens produce separate embeddings: one for the marginal distributions and another for the dependence structure. A hyper-network transforms the marginal embedding into the parameters of Deep-Sigmoide Flow, which maps the observation to its probability integral transform. An attentional copula conditioned on the dependence embedding captures cross-series interactions. The training follows a two-stage curriculum: first, learning the marginals associated with each dimension; then freezing them and fitting the copula. We use the default hyperparameters reported by Ashok et al. (2024) for fred-md (Godahewa et al., 2021), detailed in this notebook from the official code. Our setup uses 5-dimensional time-series embeddings. Both the Marginal CDF Encoder and the Attentional Copulas Encoder have two layers with one attention head (dimension 16, no dropout). The Transformer Decoder has one attention layer (dimension 8) with 3 heads, a 2-layer MI.P (hidden dimension 48), and 20 histogram bins. The Deep Sigmoidal Flow marginals (Huang et al., 2018) consist of 2 flow layers of dimension 8 and 2-layer MLP (hidden dimension 48). The training comprises 20 epochs in phase and 180 epochs in phase 2. Each time series is standardized by subtracting its training-set mean and dividing by its training-set standard deviation, ensuring zero mean and unit variance before model ingestion. Optimization parameters are kept the same as in the other baselines.
ETS. Exponential smoothing (Hyndman et al.. 2008) is a non-neural model. Holt- Winters Exponential Smoothing was applied with additive trend and additive seasonality, assuming seasonal period of 24. Model parameters were estimated through automatic optimization using the statsmodels library.
TimeMCL (Relaxed) The Relaxed- WTA Loss (Rupprecht et al., 2017) was computed for each pair (I (tag)) as
K Σ 8=1,
(1 s) C3 (1 + K
C3(
(21)
with 0.1. We found that this provides good trade-off for handling under-trainee hypotheses without deteriorating the distortion performance. However, we did not specifically tune this value for each dataset.
TimeMCL (Annealed) The annealed MCL (aMCL) Loss (Perera et al., 2024) was computed with an exponential temperature scheduler; T(t) = Top', where t is the number of the training epoch, To is the initial temperature, and p is the decay factor. each temperature T, the aMCL loss computed for each pair (I Tto:T) as
K
qu (T) (2100111111
(22)
where the coefficients
(T) 21: exp (- C -1 (.001))).))))))
(L -1, Rtg:T)
Σ cxp
T
T
(23)
are detached from the computational graph. In our experiments with aMCL, we set = 0.95 and To = 10, and we set as limit temperature Tim 5 x 10-4 before switching back to the vanilla WTA mode.
21

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Table 10. Inference time (in seconds) for EXCHANGE datasct. Results averaged over 15 random seeds. For accurate time comparison with respect to K, we disable parallel computation of the samples for each baseline. ETS, which doesn't t require neural networks is also included for completeness (in gray). See Appendix C.2.5 for details.
K
ETS
Trf. TempFlow
Tactis2
TimeGrad
DeepAR
TempFlow
TimeMCL
0.06 = 0.01
0.19 + 0.02
0,56 + 0.04
14,92 = 0.23
0.09 = 0.01
0.13 0.01
0.08 + 0.01
2
0.06 = 0.00
0.34 0.01
1,10 +10++000 0.05
30.04 = 0.49
0.13 + 0.01
0.01
0.11 + 0.01
3
0.06 + 0.00
0.48 + 0.02
+
45.10 = 0.50
0.18 + 0.01
0.00+0000
0.14 + 0.01
4
0.06 + 0.00
0.63
2.18 0.08
60.07 = 0.97
0.22 + 0.01
0.01
0.18 + 0.01
5
0.07 + 0.00
0.79 0.04
2.75 +
75.69 + 0.86
0.26 0.01
0.49
0.23 + 0.02
8
0.06 = 0.00
1.22 0.05
4.32 0.11
119.67 = 0.97
0.38 = 0.02
0.73 0.03
0.39 + 0.01
16
0.06 = 0.00
2.47 0.23
8.69 + 0.36
241.57 = 2.24
0.70 = 0.04
1.39 = 0.03
1.12 + 0.04
Table 11. Total Variation (L) comparison for K = 16 hypotheses. ETS, Trf. TempFlow and Tactis2 columns in gray because they don't share the same backbone as the other baseline. Best methods are bold, second best are underlined. Results averaged over four training seeds. Total Variation quantifies the average smoothness of the predicted trajectories (lower is smooth).
ETS
TrÉ. TempFflow
Tactie2
TimeGrad
DeepAR
TempFlow
TimeNCL (R.)
TineMCL (A.)
BLBC.
315611 1: 1858
451750 1 61031
284232
23269
372443
18499
45H430X 1: 97542
4340V7 :: 34623
220375
33961
245MIK it 30505
EXCH.
0.426
0.0052
0.633 - 0.0764
0,237 + 0.0355
0.606 - 0.0208
2,075 0.034
1,104 + 0.1552
0.831
0.0075
0.842 + 0.041
10163
3653
3391
TRAFFIC
18.782 0.016
20.56 0.92
11,063 + 0.84
12.991 + 1.661
69.603 + 2.61
18:114 0.512
5.766 + 0.19
5,714 0.2X3
TAXI
4385.41 + 11.46
5331.4= 517.54
5452.77 324.07
4232
713.9
7932.26 711.87
4147.2 362.87
099.86-1 552.07
701.61
199.38
WIKI
19692997 d: 269644
18843511 d: 1553022
2634547 J: 597323
2458503 #: 151593
9983566
434670
12410909
319517
10690 9532
271611
359320
Table 12. Results neural networks-based methods on cryptocurrency dataset. Here, K = 4, and the results were averaged over three random seeds. Here. T TimeMCL follows the same experimental setup in the previous benchmark, except that used Z-Score normalization (instead mean scaling) during training.
Trf. TempFlow
Tactis2
TimeGrad
DeepAR
TempFlov
TimeMCL (A.
Distortion
1334.441
245.746
1898.896
255.785
1834.202 1 188. 188.147
2437.798 :: 94.987
1870.915 1: 174.434
1400.396 144.5
Total Variation
8895.81
1198.262
4209.638
1047.498
12797.597
3040.59
32352.73
3054.369
9957.855
2982.198
2174.019
CRPS
0.014 + 0.001
0.018 + 0.001
0.018 0.003
0.019 0.001
0.02 0.001
0.016 0.004
RMSE
2275.986 1: 165.877
2515.755
138.424
2642.089 1: 366 366.955
2743.841 1: 93.297
2756.95 1. 68 08.01
2528.18 1 687.728
FLOPS
3.89 x10
9.81 x0"
9.13 x 10
1.74 x 10
1.94 x10"
9.98 x10"
22

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
C.4. Evaluation on financial data
An application of Time-MCL is in forecasting financial time series, such as asset prices. Rather than working directly with raw prices, it is common practice to use the log returns (Tsay, 2005), defined 2 = log P - log Pt where P represents the price of asset d at time t. These log returns are challenging to forecast. Capturing the extreme tails of their distribution can yield particularly valuable insights for financial applications.
Hourly Return Correlations
Ticker
Min.
Max.
Points.
ADA-USD
ADA
0.24
1.31
14 614
80I-080
BCH
180.45
711.92
14614
BTC-08D
BTC
25001441
108
240.09
14614
DOGE-USD
DOGE
0.06
0.48
14
BOS USD
EOS
0.40
1.50
614
ETC USD
0.6
ETC
14.52
39.67
614
ETB-USD
ETH
528.68
4076.34
14613
LINK-USD
LINK
5.77
30.44
14414
LIC-USD
0.4
LTC
51.16
144.77
607
MATIC-USD
MATIC
0.26
1.28
14 607
SOL-USD
SOL
17.49
287.81
14613
VET-USS
0.2
VET
0.01
0.08
14614
XLM USD
XLM
0.08
0.60
14614
XMR USD
XMR
101.70
241.71
14 614
XRP USD
XRP
0.40
3.39
14 614
ADA-ESD
DCR-US5
BTC-USD
DOGE-USD
BOS USD
ETC ETC USD
RTB-ISD
LINK LSD
MTC-USD
MATIC-E3D
SDI- USD
VET USD
XLMAUSD
XXXR USD
0.0
Figure 5. Cryptocurrency dataset description. (Left) Pair-wise matrix correlations of log returns between the time series, (Right) Yahoo
Finance ticker symbols, number of hourly observations per asset (Points), and the corresponding price scales (Min, Max), thereby summarizing both data volume and cross-asset dependence.
We obtained cryptocurrency data from YahooF inance and was divided into three splits: training (2023-07-01 to 2024- 09-01), validation (2024-09-02 to 2024-12-05), and test (2024-12-08 to 2025-03-01), with an hourly resolution. Figure provides overview of the assets collected for training, including return correlations among these assets, their price ranges, and the number of data points available. Missing data points were handled using forward filling.
Setup. We trained and evaluated the previously introduced neural-based baselines on this dataset, using K = 4 hypotheses. Each model was trained for 100 epochs, with 100 iterations per epoch and a batch size of 64. The prediction length and context length were both set to 24. The experimental setup was the same as in the previous benchmark, except that we applied Z-score normalization instead of mean scaling during training, which better handles the wide range of asset price scales. With the exception of Tactist is2, which already uses Z-score normalization, we observed that this change significantly improved the performance of all baselines, except for DeepAR where retained mean scaling. Note that, because the validation losses were monotonically decreasing, we used the final training checkpoint for testing.
Results. Quantitative results on quality (Distortion, CRPS, RMSE), smoothness (Total Variation), and computational cost (FLOPs) are presented Table 12. TimeMCL produces the smoothest predictions overall. Among non non-trrnsformer-bbe baselines (T meieeGad, DeepAR, and TempFlow) TimeMCL consistently outperforms the others terms of quality. In this setup, it also performs competitively with Tact 22, which was a strong competitor in the previous benchmark. Unlike in the previous benchmark, we found that the Transformer-based variant of TempFlow (Trf. TempF low) now outperforms the original TempF 1 low and slightly surpasses TimeMCL in quality metrics, though it requires nearly 40 times more FLOPs. In the future, we plan to refine this comparison by implementing a Transformer-based version of T meM...
Cryptocurrency price predictions are visualized in Figure 6. Each row represents target dimension, and each column corresponds to method listed in Table 12. Following the notation introduced in Figure, TiMeMCL predictions are shown in shades of blue, with color intensity reflecting the associated score. We observe that TimeMCL produces smoother predictions compared to other methods and effectively captures different modes in the conditional distribution. For clarity, only a subset of the cryptocurrencies is shown in the figure; however, all models were trained to jointly predict all cryptocurrencies.
23

---

Winner-takes-all for Multivariate Probabilistic Time Scries Forecasting
Oiservations
Tactie2 hypotheses
Teapflow hypothese
TimeGrad lypotheess
TineHCL hypotheses
Tactis2 mean
Trf Teapflow urean
TimeGrad menn
TimeMCL meac
Tactis2
Trf. Tempfflow
TimeGrad
TimeMCL
BCH USD
470-
100000
10000
BTC-USD SD
95000
36300
97000
0.075
BOS. USD
0.985
0.800
0.505
ETC USD
28.5
26.0
20.5
15:
1810
1m0
10.0
10.00
1240
150
10.0
18
120
Figure 6. Cryptocurrency Price Prediction on 04/01/2025 with K = 4 Hypotheses.
24

---

Winner-takes-all for Multivariate Probabilistic Time Series Forecasting
Tactli Speehom
Tacktad
Tradbool laye
Trandrad rinat
Teup Tea Tepesflnee
Trap Pae nowa
Tactis2
TimeGrad
TempFlow
TimeMCL (relaxed)
500
0.5
0.4
6.8
6470
0.0
8.0
Figure 7. Qualitative results on the TRAFFIC dataset, with same experimental setup as Figure hypotheses.
Oberrstaee
TactLad apellame
TactLod SCIME
Tranfiad Ixpethone
TiaeGrud CONE
TeupFFoe le pepetlooe
TeepFloe saics
TeNeREL Ayputionn
Tactis2
TimeGrad
TempFlow
TimeMCL (relaxed)
300
Figure Qualitative results on the ELECTRICITY dataset, with same experimental setup as in Figure with 8 hypotheses.
25


Reviews:
[{"review_id": "vz7yjdR2lu", "reviewer": "Reviewer_A4Lw", "summary": "This paper proposes the idea to generate a diverse set of forecast trajectories instead of a single trajectory as typically done. Prior ideas of doing this involved sampling from the output distribution such as in TimeGrad using a diffusion process or sampling from other models generating a distribution. However such methods do not necessarily produce diverse outputs covering the whole space of outputs.\n\nThis paper proposes using the idea of winner takes it all (WTA) to learn a tessalation of the output space (similar to K means) using multiple output heads and use it to generate several representative outputs. Experimental results show superior performance metrics compared to other comparable baselines. While the proposed methods do not perform well in terms of CRPS or RMSE, this is expected as the method does not aim at producing mean or median forecasts which minimize traditional losses.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are grateful to the reviewer for their positive feedback and insightful comments.\n\n### Sensibility to initialization\n\n> K-means is quite sensitive to initialization. Can you comment on the sensitivity of WTA to the initialization of the hypotheses?\n\nAs MCL can be seen as a conditional and gradient-based variant of K-Means, it inherits some of its limitations, in particular the sensibility to initialization. This is also related to the known collapse issue (Rupprecht et al.) in the MCL literature, where some of the hypotheses may never be chosen, leading to suboptimal Distortion performance. This is the reason why we decided to leverage WTA variants (aMCL, Relaxed-WTA), for TimeMCL, which makes the algorithm more robust to these issues. \n\nNote also that previous works (Letzelter et al., 2023; Shekarforoush et al.,2024) have noticed that this collapse issue can, in some settings, be naturally solved by the randomness of the data distribution. As mentioned in the Limitations paragraph of the submission, further work will include enhanced normalization techniques to further improve the quality of the optimum of the vanilla TimeMCL.\n\n### Motivation clarification\n\n> The motivation behind the applicability of the paper to real life forecasting problems is unclear.\n> What are some important applications for this method in real life scenarios? In what settings would someone want to generate forecasts from multiple hypothesis? \n\nIndeed, compared to generative models, we believe time-MCL has the ability to capture rare events or \"modes\" in the conditional distribution. This is for instance illustrated in Figure 2 (middle) on the Solar dataset where one of the hypotheses is capturing a rare event. This can indeed also be useful for Stock Market prediction, for capturing trend reversal. For clarification, we included a use case with financial data in the [rebuttal pdf](https://anonymous.4open.science/r/TimeMCL_ICML-E616/TimeMCL_ICML.pdf) (See e.g., Figure B) for which details are provided in the answers of Reviewer 2HFL and Reviewer uVhD.\n\nLetzelter, Victor, Mathieu Fontaine, Mickaël Chen, Patrick Pérez, Slim Essid, and Gaël Richard. \"Resilient Multiple Choice Learning: A learned scoring scheme with application to audio scene analysis.\" In NeurIPS, 2023\n\nShekarforoush, Shayan, David Lindell, Marcus A. Brubaker, and David J. Fleet. \"CryoSPIN: Improving Ab-Initio Cryo-EM Reconstruction with Semi-Amortized Pose Inference.\" In NeurIPS, 2024", "claims_and_evidence": "All methodological claims are supported by experimental results and analysis showing that they work as intended. Theoretical claims are also evidenced with proofs in the appendix.", "methods_and_evaluation_criteria": "The proposed approach is suitable for the problem studied. WTA has been shown to work for other domains such as vision and this paper shows that it works for time series forecasting as well. The paper uses the standard benchmark datasets for evaluation, and evaluations are extensive and sufficient.", "theoretical_claims": "Theoretical proofs in the supplementary material were not checked.", "experimental_designs_or_analyses": "The experimental setup is sound\n- Using synthetic time series to test the correctness of the approach is reasonable and experimental results show that the method is working as intended.\n- The results are evaluated using distortion - a reasonable metric to test the accuracy of a diverse set of outputs.", "relation_to_broader_scientific_literature": "The idea is significant and highly relevant to the time series forecasting literature. Generating a diverse set of predictions is a relevant problem in machine learning in general (especially wrt to Generative models), such ideas being extended to forecasting is highly relevant to the time series community.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "This paper presents a novel approach to model multi-modal outputs for time series forecasting. As such, i advocate for the acceptance of this paper.\n\nHowever, the motivation behind the applicability of the paper to real life forecasting problems is unclear.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "It seems that the proposed method is akin to the k means algorithm. K-means is quite sensitive to intialization. Can you comment on the sensitivity of WTA to the initialization of the hypotheses?\n\nIt is possibly the case that the loss is such (convex or some similar property) that any initialization would lead to the hypothesis shifting appropriately and aligning with the actual distribution during training.\n\nWhat are some important applications for this method in real life scenarios? In what settings would someone want to generate forecasts from multiple hypothesis? An examples that come to mind is stock market predictions - it would be important to look at a diverse set of forecasts to understand any extreme predictions. Or product demand forecasts where a retail company may want to prepare for different scenarios.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0OS0ZVM2AS", "reviewer": "Reviewer_2HFL", "summary": "This work introduces TimeMCL, a time series forecasting model that looks to project plausible future scenarios and their associated probabilities to better forecast multimodal distributions. The model learns multiple heads, as well as scores associated to each head to estimate the probability of a given head being correct. Training these heads with vanilla Winner-takes-all (WTA) loss can result in under-trained heads for unimodal distributions, so the authors explore relaxations that enable more distributed gradients. The authors provide a theoretical analysis of TimeMCL, demonstrating that their combined architecture and training scheme leads to a voronoi tesselation of future trajectories, which they show on synthetic datasets. They also conduct experiments on six typical TS datasets from GluonTS wrt Distortion, FLOPs, RMSE and CRPS.\n\nThe main claimed contributions are:\n- the TimeMCL approach that takes a backbone and trains multiple heads for it using (two variations of relaxed) Winner-Takes-All loss (with scoring heads)\n- the theoretical analysis of TimeMCL as a functional quantizer\n- the evaluation of TimeMCL with an RNN backbone on synthetic and real-world benchmarks", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their feedback. The [rebuttal pdf](https://anonymous.4open.science/r/TimeMCL_ICML-E616/TimeMCL_ICML.pdf) is attached to the response. \n\n### Clarification of the theory\n\n> TimeMCL is a stationary conditional functional quantizer [...], but it's unclear what value this provides as the underlying clustering problem is NP-hard, and you're dealing with very high-dimensional data.\n\nThe reviewer raises a valid point. However, as shown in Figure 1, our toy example qualitatively demonstrates that our training scheme closely aligns with the target conditional quantizer in practice, even when forecasting across 250 time steps (see also Appendix B for details).\n\n> Maybe [...] this algorithm converges to the optimal voronoi tesselation [...], but it would help to analyze the rate of convergence theoretically and experimentally.\n\nTheorem 2 of Loubes & Pelletier (2017) provides an asymptotic upper bound on the distortion error with respect to the number of training pairs in quantizer learning. Extending this result to neural networks trained with WTA Loss is a promising direction for future work.\n\n### Additional baselines and metrics \n\n> A runtime analysis would be nice [...].\n> It would be good to add in other relevant baselines, e.g. TACTiS and simple baselines\n> You could also assess the models using actual multivariate metrics, e.g. energy score [...]\n> It would be nice to also provide [...] a measure of diversity, [...] and smoothness.\n\nIn response to the reviewer, we added experiments with a simple baseline, ETS exponential smoothing (Hyndman et al., 2008), which does not rely on neural networks. We also included TempFlow (Rasul et al., 2021), a normalizing flow-based using both RNN and Transformer backbones, as well as Tactis-2 (Ashok et al., 2024), a Copula method with a Transformer backbone.\n\nThese methods were evaluated across the same six datasets. In addition to the metrics from the original paper (Distortion, RMSE, and CRPS-Sum), we followed the reviewers' suggestions and included Inference runtime, Smoothness, and Energy Score.\n\n### Results analysis\n\n**Distortion Comparison (Table A,E,F)**.\n\n* TimeMCL outperforms TempFlow, and Tactis proves to be a strong competitor in terms of Distortion, though being slower at inference (see Table H).\n\n* An ablation study on the number of hypotheses (Table E) shows that TimeMCL consistently achieves the best performance (except with $K=1$).\n\n**Inference run-time (Table H)**\n\n* Among neural network-based methods, TimeMCL and DeepAR demonstrate the best trade-off between speed and performance.\n\n* ETS achieves the fastest inference, with weaker performance otherwise.\n\n**Smoothness Analysis (Table B)**\n\n* TimeMCL achieves the best smoothness scores, as measured by Total Variation (averaged over predictions).\n\n* This supports our theoretical claim from Section 5.2, which predicts that TimeMCL generates smoother trajectories.\n\n**Additional metrics (Tables C & D)**\n\n* TimeMCL does not significantly improve RMSE, CRPS, or the Energy Score (Table G), since it does not directly optimize these metrics.\n\nWe did not conduct further experiments with Diversity, as we believe Distortion implicitly captures it.\n\nWhen comparing TimeMCL (using an RNN backbone) with methods using transformer backbones, it's hard to tell if performance improvements are due to the training method or the backbone. To clarify, we plan to implement TimeMCL with a transformer backbone, following Tactis's training details.\n\n### Additional datasets\n\n> It is difficult to know how much multivariate correlation is present in these datasets. Therefore, performance on synthetic multivariate tasks would be insightful\n\nIn response, we conducted experiments on a new dataset of correlated financial cryptocurrencies time series (see Table J), with the correlation matrix in Table K.\n\nTimeMCL was trained with an aMCL loss and compared to previous baselines. Results in Table I show that our method remains competitive, excelling in both Distortion and Smoothness, with strong CRPS performance. Visualizations are in Figure B.\n\n### Additional questions\n\n> Can the model accept covariates?\n\nThe model can accept covariates. In previous implementations these typically serve as additional concatenated input features.\n\n> To compute TimeMCL metrics [...], we resample with replacement from the K hypotheses obtained [...]. Can you explain this is more detail please?\n\nOur implementation extends GluonTS (Alexandrov et al., 2020) with minimal code changes. Instead of modifying evaluation functions, we resample from the K hypotheses, weighted by their probabilities, allowing us to use existing evaluation functions (e.g., CRPS) without rewriting metrics for TimeMCL.\n\n> Did you look at methods that might enable a variable number of heads? \n\nIndeed, the number of predictions must be predefined beforehand. Exploring dynamic \"rearrangements\" of hypotheses when adding new ones, without full retraining, is left for future work.", "claims_and_evidence": "> TimeMCL forecasts diverse possible futures & provides smooth forecasts\n\n  - The is visualized in figures 2-4. It would be nice to also provide some sort of quantitative measure of diversity here, e.g. Frechet or even just dispersion averaged across series and time points. The same goes for smoothness.\n\n> TimeMCL is a stationary conditional functional quantizer\n\n  - This seems to be true under the proposed assumptions, but it's unclear what value this provides as the underlying clustering problem is NP-hard, and you're dealing with very high-dimensional data.\n\n> TimeMCL is compared against SOTA probabilistic forecasters\n\n  - The TACTiS models that you cited are not compared against, which is odd.\n    - https://arxiv.org/abs/2202.03528\n    - https://arxiv.org/abs/2310.01327\n  - Missing simple baselines against which to compare on the real-world datasets, e.g. naive, drift, ARIMA, exp-smoothing, etc.", "methods_and_evaluation_criteria": "The datasets are typical for ML forecasting, although it is difficult to know how much multivariate correlation is present in these datasets. Therefore, performance on synthetic multivariate tasks would be insightful, e.g. correlated brownian motion and/or VAR processes.", "theoretical_claims": "- I reviewed the proof of proposition 5.1 briefly. I am mostly unfamiliar with the clustering literature, but from what I gather, k-means is NP-hard, so I'm not sure what this formulation actually brings given the extremely high dimensionality of multivariate trajectory data. Maybe in the limit of expressivity and time, this algorithm converges to the optimal voronoi tesselation due to the two-step formulation, but it would help to analyze the rate of convergence theoretically and experimentally.\n\n- It's unclear why the architectures for sections 5 and 6 differ.\n\n- I did not review the proof of proposition 5.2", "experimental_designs_or_analyses": "- A training/inference runtime analysis would be nice in addition to flops\n\n- It would be good to add in other relevant baselines, e.g. TACTiS and simple baselines\n\n- You could also assess the models using actual multivariate metrics, e.g. energy score, variogram, etc. (see https://arxiv.org/abs/2407.00650)", "relation_to_broader_scientific_literature": "This paper relates to a recent line of work about time series forecasting using machine learning methods. Specifically, the paper belongs to two niches: multivariate forecasting and probabilistic forecasting. Cited methods, included DeepAR, TimeGRAD and TACTiS, figure among this line of work as well.", "essential_references_not_discussed": "- Missing reference to distortion being standard in quantization\n\n- Other references to multivariate probabilistic forecasting:\n  - https://arxiv.org/abs/2410.02168\n  - https://ojs.aaai.org/index.php/AAAI/article/view/29085", "other_strengths_and_weaknesses": "- Extremely clear and well-written paper\n\n- The main points that can move my score are:\n  - experiments with multivariate synthetic tasks\n  - multivariate metrics (energy score, variogram)\n  - comparison to TACTiS-2", "other_comments_or_suggestions": "Typos: \n\n- paragraph at line 84 has some weird verb conjugation, e.g. \"one have\"\n- missing period line 100-101 col 2 before \"WTA\".\n- the notation, and especially the overload of x, feels a bit clumsy, e.g. on line 112-113 you switch to the superscript for the hypotheses, which took a minute to notice since the subscript is omitted (I imagine to avoid it being too heavy). You can probably just use \\mathbf{x} for a vector to prevent the cognitive break during reading.\n- line 120 col 1: you might want to reiterate that these are the \"final projection\" heads here instead of just \"hidden state representation\"\n- Tables 6 and 7 belong in the main text, as they are key evaluation results with more established metrics than distortion.\n- Line 167 col 2: \"finite over if\" \n- Line 171 col 2: the first clause of this paragraph is incomprehensible.", "questions_for_authors": "> covariates c1:T , the latter being omitted in the following for conciseness \n\ncan the model accept covariates?\n\n---\n\n> To compute TimeMCL metrics while respecting hypothesis probabilities, we resample with replacement from the K hypotheses obtained in a single forward pass, weighting them by their assigned probabilities before computing metrics.\n\nCan you explain this is more detail please? I'm not sure I follow, I figured that the score heads provided the probabilities.\n\n---\n\nDid you look at methods that might enable a variable number of heads? Seems like the optimal tesselation is conditional on the number of heads, so you can't really \"extend\" the number of heads.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "jNisv2SFQI", "reviewer": "Reviewer_uVhD", "summary": "The paper presents a new method, TimeMCL, for time series forecasting. The proposed method uses Multiple Choice Learning using Winner-Takes-All (WTA) loss to forecast multiple plausible time series future. The paper uses synthetic data to show that TimeMCL is a functional quantizer. The proposed TimeMCL is compared with two Diffusion methods on standard datasets.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their relevant remarks. The [rebuttal pdf](https://anonymous.4open.science/r/TimeMCL_ICML-E616/TimeMCL_ICML.pdf) is attached to the response. Figure A of the pdf will be included in the paper.\n\n### Comparison with more SOTA methods \n\n> The results will be more validated if compared with more SOTA methods.\n\nWe conducted this benchmark using DeepAR and TimeGrad to demonstrate that, with consistent settings across baselines (e.g., backbone, data scaler, training details), our approach offers competitive distortion at a low computational cost. Since we used an RNN backbone, we felt comparing it with other architectures, such as transformers, would complicate conclusions. However, we acknowledge the reviewer’s point that additional baselines could strengthen our work.\n\nTo address the reviewer's comment and enhance our evaluation, we added additional models. Specifically, we included Tactis-2 (Ashok et al, 2024), a transformer-based model based on non-parametric copulas, and TempFlow (Rasul et al, 2020), which uses conditioned normalizing flows (with both RNN and transformer backbones). For completeness and as suggested by Reviewer 2HFL, we also included exponential smoothing (ETS) as a simple baseline without neural networks.\n\n### Results analysis\n\n**Distortion Comparison (Table A)**.\n\n* TimeMCL outperforms TempFlow, both when using the same RNN backbone and when TempFlow is based on a Transformer.\n\n* Tactis proves to be a strong competitor in terms of Distortion, though at a significantly higher computational cost (see Table H).\n\n* We conducted an ablation study on the number of hypotheses (Table E). We observed that TimeMCL consistently achieves the best performance (except when using only one hypothesis).\n\n**Inference run-time (Table H)**\n\n* Among neural methods, TimeMCL and DeepAR demonstrate the best trade-off between speed and performance.\n\n* ETS achieves the fastest inference but exhibits weaker performance on other metrics.\n\n**Smoothness Analysis (Table B)**\n\n* TimeMCL achieves the best smoothness scores, as measured by Total Variation.\n\n* This supports our theoretical claim from Section 5.2, which predicts that TimeMCL generates smoother trajectories.\n\n**Additional metrics (Tables C & D)**\n\n* * TimeMCL does not significantly improve RMSE, CRPS, or the Energy Score (Table G), as expected, since it does not directly optimize these metrics.\n\nWe plan to extend this comparison by implementing timeMCL with a transformers-based architecture as the backbone, and adhering to Tactis's training details for more accurate performance comparison.\n\n> How are the graph tranformer methods, e.g. STGNN, compare with the diffusion methods in general and with TimeMCL in particular?\n\nRegarding spatio-temporal graph transformer methods, does the reviewer have a specific method in mind? Most STGNN methods we found are tailored to specific tasks (e.g., traffic prediction in Luo et al., 2023). However, we did identify stemGNN (Cao et al., 2020), which integrates graph and attention mechanisms and is evaluated on similar data.\n\nHowever, comparing our approach with stemGNN is difficult, as it is non-probabilistic and generates only one prediction per input. As future work, exploring a graph transformer-based approach in place of our RNN backbone could be promising.\n\n### Need for datasets from more domains, e.g., Stock markets\n\n> How will the proposed approach perform in other domains e.g. Stock Market?\n\nWe thank the reviewer for their advice. Accordingly, we performed experiments on a dataset of correlated financial time series, consisting of 15 correlated Cryptocurrencies (See Tables J and K of the rebuttal pdf).\n\nOn this dataset, TimeMCL was trained with annealed winner takes all loss, and compared with the previous baselines, and the results, presented in Table I further demonstrate the competitiveness of the method in terms of Distortion and Smoothness, with also good CRPS here.\n\nWe also provide a visualization of the predictions, along with the baselines in Figure B, which is akin to Figure 2 of the main paper.\n\n### Missing references  \n\n> I believe other SOTA works must be discussed including the Graph Deep learning methods e.g. STGNN.\n\nTo make our benchmarks more comprehensive, we’ve included Tactis2, TempFlow (with its transformer-based variant), and non-neural exponential smoothing (ETS) methods (Hyndman et al., 2008). Additionally, STGNN will be referenced in the paper.\n\nLuo, X., Zhu, C., Zhang, D., & Li, Q. (2023). Stg4traffic: A survey and benchmark of spatial-temporal graph neural networks for traffic prediction\n\nCao, Defu, Yujing Wang, Juanyong Duan, Ce Zhang, Xia Zhu, Congrui Huang, Yunhai Tong et al. \"Spectral temporal graph neural network for multivariate time-series forecasting.\" In NeurIPS 2020\n\nRasul, Kashif, Abdul-Saboor Sheikh, Ingmar Schuster, Urs M. Bergmann, and Roland Vollgraf. \"Multivariate Probabilistic Time Series Forecasting via Conditioned Normalizing Flows.\" In ICLR 2021", "claims_and_evidence": "1. The main claim that TimeMCL is a functional quantizer is supported by mathematical proofs and the use of synthetic data\n2. The forecasting ablility of TimeMCL is supported by comparing with 2 diffusion methods", "methods_and_evaluation_criteria": "The proposed method is evaluated on standard datasets. The comparisons are done with 2 SOTA diffusion methods.\n\nThere should be more comparisons with other SOTA methods", "theoretical_claims": "The paper presented proofs and used synthetic data to support TimeMCL as a Functional Quantizer. I did not check the correctness of the mathematics and proofs.", "experimental_designs_or_analyses": "The datasets used are usually the ones that are used in SOTA time series forecasting methods. Two SOTA diffusion methods are used for comparisons. The results will be more validated if compared with more SOTA methods.", "relation_to_broader_scientific_literature": "The paper proposes a WTA approach for time series forecasting and I believe does not relate to broader scientific literature.", "essential_references_not_discussed": "The works referenced in the paper are adequate. Although I believe other SOTA works must be discussed including the Graph Deep learning methods e.g. STGNN.                  .", "other_strengths_and_weaknesses": "Strengths:\n1. Use of synthetic data \n2. Mathematical equations and proofs\n3. The use of multiple performance criteria\nWeaknesses\n1. Need to consider datasets from more domains e.g. stock market\n2. Need to compare with other SOTA time series forecasting methods\n3. In comparing with TimeGrade and DeepAr, there is an assumption that these are the 2 best methods for time series forecasting.", "other_comments_or_suggestions": "1. Please add a figure that shows the network architecture of the proposed method and a figure that visualy describes the proposed technique.\n2. Figures and tables should be closer to the text descriptions e.g. figure 1 is on a different page than the description, same is true of table 2 & 3.", "questions_for_authors": "1. How will the proposed approach perform on other domains e.g. Stock Market?\n2. How are the graph tranformer methods, e.g. STGNN, compare with the diffusion methods in general and with TimeMCL in particular?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {"ethical_review_concerns": "NA"}}, {"review_id": "eBnMdjezfG", "reviewer": "Reviewer_zTPp", "summary": "This paper addresses a time-series forecasting problem, where the model generates multiple forecasting for each timestamp. The authors propose TimeMCL, a method based on the Multiple Choice Learning (MCL) paradigm that can output multiple plausible forecasting with multiple prediction heads and score heads. TimeMCL is trained with Winner-Takes-All (WTA) loss and score head loss to produce diverse forecasting from multiple heads due to the nature of the WTA loss, computing the gradient for the head with minimum loss value. The authors claim that TimeMCL can be viewed as a functional quantizer with the theoretical analysis. Experimental results on multiple datasets demonstrated that the proposed method performs better than the baselines in terms of the Distortion metric while it performs comparably with the baselines in terms of the standard metrics.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their positive feedback on the paper.\n\n> In Eq.5, why is not $x_{t-1}$ inputted to the function gamma?\n\nThis is because, with our notations, $\\gamma^k_{\\theta}$ corresponds only to the head. The full model writes as $\\gamma_\\theta^k \\circ s_\\theta$, where $s_{\\theta}$ is the backbone (See section 3,4 and Figure A of the [rebuttal pdf](https://anonymous.4open.science/r/TimeMCL_ICML-E616/TimeMCL_ICML.pdf)).\n\n> In l.134 of the right-hand side in p.2, why is not $x_{t-1}$ inputted to the function $f^k_\\theta$?\n\nSame as above. This is because, with our notations, $f^k_{\\theta}$ Only the head. The full model writes as $f_\\theta^k \\circ s_\\theta$, where $s_{\\theta}$ is the backbone. This will be made clearer with an illustration as suggested by Reviewer uVhD. \n\n> The authors mentioned that the score head can avoid overconfident heads, but how they avoid them is not described.\n\nOverconfidence is a known issue in Multiple Choice Learning (Rupprecht et al., Lee et al), where, some heads associated with very low probability zones may not be distinguishable from plausible hypotheses at inference time. Score-heads allow to solve this issue by learning the probability of each head to avoid this situation.   \n\n> Can we not replace min in the WTA loss with the gamma and jointly train with the WTA loss, which can be a more straightforward approach? It can be similar to Eq.9.\n\nDoes the reviewer suggest training a loss that looks like a weighted sum of the $L_{\\theta}^{k}(x_{1:t_{0}-1},x_{t_{0}:T})$ by the $K$ predicted scores? The suggestion of the reviewer is very interesting, as it would when the hypotheses are fixed, encourage the score associated with the lowest distance to increase (and the other to decrease). However, this would make the prediction head loss dependent on the values of the scores, which may produce different training dynamics as the current model. In the current version, the score head objective depends on the position of the hypotheses, but not the opposite. This would definitely be a promising try for further work.                \n\n> Proposition 5.2. is analyzed only with the binary cross entropy for the score heads. However, s in TimeMCL are shared in the score heads and heads, and TimeMCL is trained based on the compound loss with the WTA loss, which is not the direct case of Proposition 5.2.\n\nIndeed, Proposition 5.2 is analyzed only with the binary cross entropy for the score heads. In this proposition, we implicitly assume that the prediction heads have already converged, so that the task of learning the probability mass of each cell becomes duable for the prediction heads. \n\nIndeed, in according with Letzelter et al. (2024b) (Section C.1.2), we observed the WTA training scheme leads to a fast convergence of the predictions, while the scoring heads are slightly slower to train because they need the prediction heads to have already converged to do so. \n\nThis will be made clearer in the assumption of this Proposition 5.1.\n\n> Why is the distortion a fair metric ? \n\nThe Distortion is known from the quantization literature (Pagès.G.,2015), as a way to assess the quality of a set of $K$ samples $z_k$, $k = 1,...,K$, for quantizing a target distribution $p$, with:\n\n$$ D_2 := \\int_{\\mathcal{X}} \\min_{k=1, \\ldots, K} \\left\\| z_k - x \\right\\|_{2}^2 \\mathrm{d}p(x). $$\n  \nIn our setup, the distortion we are considering is the generalization of the above where the samples $z_k$ are functions of the context. It writes in the context of time series\n\n$$ \\int_{\\mathcal{X}^{T}} \\min_{k=1, \\ldots, K} \\left\\| z_k\\left(x_{1: t_0-1}\\right) - x_{t_0: T} \\right\\|^2 \\; \\mathrm{d}p(x_{1: t_0-1}, x_{t_0: T}) \\simeq \\frac{1}{N} \\sum_{i} \\min_{k=1, \\ldots, K} \\left\\| z_k\\left(x^i_{1: t_0-1}\\right) - x^i_{t_0: T} \\right\\|^2. $$\n\nwhere $(x^i_{1: t_0-1},x^i_{t_0: T})$ are samples from $p(x_{1: t_0-1},x_{t_0: T})$. \n\nIt implicitly assesses how well the predicted samples cover a target distribution with a given set of samples. The distortion is a fair metric provided that the same number of hypotheses is used for each baseline, as the distortion is expected to improve with $K$, as per the Rate distortion curve (Gray, 1989) in the optimal case.\n\nGray, Robert M. Source coding theory. Vol. 83. Springer Science & Business Media, 1989.", "claims_and_evidence": "* In Eq.5, why is not x_{t-1} inputted to the function gamma?\n* In l.134 of the right-hand side in p.2, why is not x_{t-1} inputted to the function f^k_theta?\n* The authors mentioned that the score head can avoid overconfident heads, but how they avoid them is not described.\n* Can we not replace min in the WTA loss with the gamma and jointly train with the WTA loss, which can be a more straightforward approach? It can be similar to Eq.9.\n* Proposition 5.2. is analyzed only with the binary cross entropy for the score heads. However, s in TimeMCL are shared in the score heads and heads, and TimeMCL is trained based on the compound loss with the WTA loss, which is not the direct case of Proposition 5.2.", "methods_and_evaluation_criteria": "* Why is the Distortion fair metric? How is the Distortion computed for the baselines? Can we use score head to choose the hypotheses and use the standard metrics?\n* The paragraph \"Comparing TimeMCL with the baselines on standard\nmetrics\" describes only Tables 6 and 7 in the Appendix. It cannot be appropriate in terms of the page limit regulation.", "theoretical_claims": "* The descriptions around l.215 of the right-hand side in p.4 are not self-complete. For example, z is not defined.", "experimental_designs_or_analyses": "Please see Methods And Evaluation Criteria.", "relation_to_broader_scientific_literature": "TimeMCL, a method based on the Multiple Choice Learning (MCL) paradigm, can be novel and practically important.", "essential_references_not_discussed": "NA", "other_strengths_and_weaknesses": "Clarity issues:\n* In l.82 of the right-hand side in p.2, the subscripts may be wrong in the equation.\n* In l.87 of the right-hand side in p.2, the right most \"...\" may be unnecesally.\n* The character T is used in multiple meanings in forecasting horizon and annealed temperature.\n* Figures 1 and 2 appear in reversed order.", "other_comments_or_suggestions": "NA", "questions_for_authors": "See above", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper introduces TimeMCL, a novel method for multivariate probabilistic time series forecasting that uses a multiple-choice learning approach with a winner-takes-all loss to generate diverse future predictions. Reviews highlight its ability to produce a range of plausible forecasts and its theoretical grounding as a functional quantizer. While it shows strong performance on a custom distortion metric designed for diverse outputs, it doesn't necessarily excel on traditional metrics like RMSE and CRPS, as it prioritizes diversity over mean accuracy.

Pros include its innovative approach to generating diverse forecasts, supported by theoretical analysis and promising results on a dedicated metric. Reviewers also noted the clarity of the writing. Cons involve weaker performance on standard forecasting metrics and the need for more comparisons with a broader range of state-of-the-art methods, and a clearer motivation for real-world applicability. Some theoretical aspects and experimental details also required clarification.

The authors addressed concerns by providing additional experiments with more baselines and on a new financial dataset, including multivariate evaluation metrics, and clarified theoretical points and the model's architecture. They also discussed potential real-world applications and acknowledged limitations like the fixed number of prediction heads.

Considering the novelty of the approach, the theoretical backing, and the authors' efforts to address the reviewers' concerns with additional experiments and clarifications, the unanimous consensus is for acceptance. The method offers a different perspective on time series forecasting by focusing on capturing a spectrum of possibilities, which could be valuable in specific applications despite not optimizing traditional error metrics.
Decision: Accept (poster)