Multivariate Time Series Forecasting needs Cross Variable Loss

The paper addresses a fundamental limitation in multivariate time series forecasting (MTSF). The authors note that while existing research has focused on input-side channel dependence modeling—extracting representations from historical observations—the intrinsic correlations among future predictions at the output-side have been largely overlooked. Modern forecasting models follow the Direct Forecasting (DF) paradigm, which typically decouple the multivariate outputs forcibly, projecting the unified historical embeddings independently into future time steps, thereby degrading MTSF into a collection of isolated scalar prediction problems.

The authors identify dual deficiencies in this approach: On one hand, blindly decoupling prediction targets directly violates the intrinsic physical and operational laws of multivariate systems (since real-world variables co-evolve under shared system dynamics), and Furthermore, overlooking these synergistic relationships induces critical structural errors at the optimization level.

At the optimization level, the paper demonstrates that point-wise loss functions (like MSE) are "fundamentally predicated on a highly restrictive assumption: the prediction residuals across all variables and time steps follow completely independent and isotropic Gaussian distributions (i.e., possessing a diagonal covariance matrix). However, given the aforementioned co-evolutionary phenomena, the true spatio-temporal covariance matrix of the residuals inevitably contains numerous non-zero off-diagonal elements. This creates what the authors term an Objective Gap" between standard loss functions and the true joint spatio-temporal distribution of the data.

The paper formalizes this gap in Theorem 3.1:

"Let e = vec(Ŷ − Y) ∈ RT D denote the flattened prediction error vector... Suppose the residuals follow a zero-mean multivariate Gaussian distribution e ∼ N (0, ΣST)... the scaled negative log-likelihood (NLL) under the true covariance ΣST is Lnll = 2T1D e⊤ Σ−1 ST e. On the other hand, the standard DF objective in Eq. (1) can be written as Ldf = T1D e⊤ e. Therefore, compared with the isotropic Gaussian NLL under variance σ 2, the discrepancy between the two objectives is ∆ = Lnll − 2σ1 2 Ldf = 2T1D e⊤ (Σ−1 ST − σ12 I)e."

The theorem establishes that the standard point-wise DF objective coincides with the isotropic Gaussian NLL only when the residual covariance is spherical, i.e., ΣST = σ 2 I, meaning a strictly non-zero objective gap (∆ ̸= 0) inherently exists because the standard point-wise DF objective neglects these underlying structural correlations.

The authors carefully scope their theorem: It establishes that some non-diagonal term is necessary... It does not identify which non-diagonal term should be used. They clarify that estimating Σ−1 ST directly is intractable: On ECL with T = 720 and D = 321 that matrix has more than 5 × 10 10 entries, so the estimation problem is ill-posed.

Proposition 3.2 then bridges the gap by showing that MSE plus a squared graph total-variation penalty over E is the negative log-likelihood of a Gaussian whose precision matrix has support E off the diagonal. Specifically, for the structured precision family Σ−1(λ) = σ12 I + λ Lap(E), the identity holds:

r⊤ Σ−1 (λ) r = σ1 squared ∥r∥22 + λ P (i,j)∈E ∥ei − ej ∥22

This identifies exactly the object that Theorem 3.1 shows MSE to be missing, confirming that E is not an arbitrary analogy but the conditional-independence structure of the residual field.

The method splits the time axis into P non-overlapping patches of length L (T = P × L), defining each node as v = (p, d) where p indexes patches and d indexes variables. The residuals are computed as ev = ẑv − zv ∈ RL per node.

Three graph topologies are designed: (1) Synchronous edges which strictly capture concurrent interactions at the same time step (i.e., connecting identical patch indices); (2) Asynchronous edges which exclusively model lagged dependencies across different time steps; and (3) Fully connected edges which include both synchronous and lagged interactions. The fully connected graph is the default to maximize the coverage of potential interactions without relying on domain-specific prior knowledge.

Using the incidence matrix à = A ⊗ IL, CvLoss is defined as:

Lcv = E×L 1 P (i,j)∈E ∥ei − ej ∥1 = E×L 1 ∥Ã vec(E)∥1

This computes the total variation of prediction errors across the graph. The final objective combines both losses:

Lα:= (1 − α) · Ldf + α · Lcv

The paper explains that Ldf encourages each residual ev to be small in absolute magnitude, while Lcv encourages connected nodes to have compatible residual patterns. This induces an explicit coupling across variables at the objective level.

Regarding the norm choice: We use the l1 norm because it is less sensitive to a few large edge-wise mismatches caused by localized anomalies or variable-specific shocks. The authors explicitly state the l1 form is a task-specific empirical choice rather than a consequence of the proposition; Appendix C.4 compares the two norms directly under an identical protocol.

The method is evaluated on twelve standard benchmarks (four ETT subsets, Weather, ECL, Traffic, Solar, and four PEMS datasets). The primary results use TimeFilter as the backbone: CvLoss is TimeFilter, trained on the same data with the same protocol and differing only in the loss term, and it improves all 24 dataset-level entries (12 datasets × 2 metrics) with no tie and no loss.

The paper emphasizes backbone-controlled attribution: "Across 114 such comparisons, spanning seven backbones from three architecture families (Transformer, linear, and hybrid decomposition designs), twelve datasets and both metrics, CvLoss improves 111, ties 2 and loses 1, with a mean relative error reduction of 3.39% in MSE and 3.61% in MAE. The effect is therefore a property of the objective rather than of any one architecture."

Notable gains include Traffic on iTransformer, −17.7% MSE; PEMS07 on TimeFilter, −11.3% while effects are small where either condition fails (ETTh1, −0.24% MSE on TimeFilter).

CvLoss is compared against 11 learning objectives across four categories (shape alignment, likelihood maximization, distribution balancing, and decomposition). The results show CvLoss ranks first in 20 settings and that alternative objectives improve over MSE, but still treat multivariate outputs as scalar predictions. They therefore leave the spatiotemporal objective gap (∆ ̸= 0) in Theorem 3.1 unaddressed.

The ablation study compares synchronous-only (CvLoss†), asynchronous-only (CvLoss‡), and fully connected graphs: the default CvLoss generally outperforms both CvLoss† and CvLoss‡ in most settings. This confirms that integrating both synchronous and asynchronous interactions within E is crucial for capturing complex dependencies.

Sensitivity analysis on α shows Increasing α from 0 yields consistent improvements, but using the structural penalty (α = 1) is not always optimal; on ECL, both models perform best around α = 0.5 and degrade when α reaches 1, suggesting CvLoss is most effective when combined with point-wise MSE, which anchors forecasting accuracy and prevents translation-invariant solutions. Patch length sensitivity shows the predictive accuracy remains consistent regardless of the specific patch length chosen, indicating that CvLoss is highly robust to this hyperparameter.

CvLoss improves multiple architectures (TimeFilter, TQNet, CFPT, iTransformer) across datasets. The method imposes zero overhead during inference, ensuring testing speed remains completely unaffected while maintaining extremely high predictive accuracy. A random edge sampling mechanism reduces computational complexity during training from quadratic O(P2D2) to an effectively constant level, achieving large speedups with negligible performance degradation (e.g., 4397× speedup on ECL).

The paper concludes: "This paper studies the overlooked output-side dependencies in multivariate time series forecasting. We show that the standard point-wise DF objective is mismatched with real-world spatiotemporal correlations, and introduce CvLoss to explicitly regularize cross-variable residual consistency. CvLoss is implemented as a graph total variation penalty over forecast patches, allowing both synchronous and lagged interactions to be modeled without changing the inference procedure."

The authors state their results suggest that modeling the structural geometry of future variables is a simple but effective direction for improving multivariate forecasting.

The authors acknowledge: "CvLoss uses a predefined cross-variable graph and a fixed patching scheme. Although random edge sampling improves scalability, the graph may still include uninformative interactions, especially for very high-dimensional systems with sparse or domain-specific dependencies. They also note the implemented loss is not the object that Theorem 3.1 literally implicates since the l1 penalty defines a pairwise Gibbs field over the same graph that is not Gaussian, and the default uniform fully connected graph assigns the same weight to every cross-variable pair, whereas real dependencies are heterogeneous. Extending to probabilistic forecasting, irregular sampling, and missing-value scenarios is left for future work."

Based on the paper, I would improve multivariate time-series forecasting AI systems in the following concrete ways.

• Keep the standard MSE/MAE term, but add a graph total-variation penalty over the residual field:

L = (1 - α) · L df + α · L cv

• Define residual nodes as (patch, variable) rather than isolated scalar time steps, and connect nodes with synchronous edges (same time, different variables) and asynchronous edges (lagged cross-variable dependencies).

• Use the l1 total-variation form to make the loss robust to localized anomalies and variable-specific shocks.

• Compute residuals e = ŷ - y over forecast patches.

• Enforce that connected residuals have similar patterns: Σ(i,j)∈E e i - e j 1.

• This directly targets the objective gap identified in the paper, where standard MSE implicitly assumes independent/isotropic residuals.

• Connect every pair of variables across every patch pair to capture both synchronous and lagged interactions without requiring prior knowledge.

• If domain structure is known, restrict the graph to plausible edges to reduce noise and improve interpretability.

• Sample a subset of edges each training iteration to avoid the quadratic cost of full P2D2 edge enumeration.

• This reduces training overhead to an effectively constant level while preserving most of the benefit.

• Apply the new loss only during training.

• The improved system runs at the same inference speed as the original model, with zero extra latency or memory overhead at deployment.

• Use a moderate α (e.g., 0.5) rather than pure structural regularization, because the point-wise MSE anchors accuracy and prevents translation-invariant degenerate solutions.

• Patch length is robust, so it can be set based on seasonality or computation constraints.

• Forecast multiple correlated variables more accurately: It reduces MSE and MAE across diverse multivariate benchmarks, with consistent gains across transformer, linear, and hybrid backbones.

• Model inter-variable dependencies at the output side: Instead of predicting each variable independently, the system explicitly regularizes predictions so that residual patterns across variables and time patches are mutually consistent.

• Capture synchronous and lagged relationships: The graph-based penalty handles both same-time cross-variable correlations and time-delayed dependencies, which are common in real physical, economic, and traffic systems.

• Scale to high-dimensional systems: Random edge sampling allows training on datasets with hundreds or thousands of variables (e.g., electricity, traffic sensor networks) without prohibitive computational cost.

• Remain architecture-agnostic: The loss is a drop-in replacement for standard MSE/MAE, so existing forecasting models can be improved without redesigning their internal architecture.

• Maintain inference efficiency: Because the penalty is only in the training objective, the deployed system produces forecasts at the same speed as before, making it suitable for real-time applications.

• Generalize to other spatiotemporal prediction tasks: The same principle applies to any AI system that outputs a tensor of correlated predictions—e.g., weather field prediction, energy load forecasting, supply chain demand, or multi-agent trajectory prediction—where the loss should reflect the joint structure of the output space.

Multivariate time series forecasting presents unique challenges because future variables often co-evolve under shared system dynamics. While existing studies mainly focus on cross-variable dependencies in historical observations, dependencies among future values are much less explored. Specifically, modern forecasting models largely follow the Direct Forecasting (DF) paradigm, generating multi-step forecasts with point-wise objectives that do not explicitly constrain cross-variable structure. In this work, we show that the DF objective is mismatched in the presence of cross-variable and lagged dependencies, revealing an objective gap. To address this issue, we propose Cross-Variable Loss (CvLoss), a plug-in structural regularizer that constrains forecast residuals on a cross-variable graph. CvLoss penalizes inconsistent edge-wise residual differences over forecast patches, encouraging consistency across both synchronous and asynchronous interactions. Our experiments show that CvLoss consistently improves competitive forecasting models, outperforms representative learning objectives, and is compatible with a variety of forecasting backbones.