Multivariate Time Series Forecasting needs Cross Variable Loss

# Multivariate Time Series Forecasting needs Cross Variable Loss
Source: [https://arxiv.org/html/2608.05742](https://arxiv.org/html/2608.05742)
Kuiye Ding1,Yifan Hu2,Hanchen Wang1,Hao Xue3 1University of Technology Sydney,2Tsinghua University, 3The Hong Kong University of Science and Technology \(Guangzhou\) kuiye\.ding@student\.uts\.edu\.au, hanchen\.wang@uts\.edu\.au huyf25@mails\.tsinghua\.edu\.cn, haoxue@hkust\-gz\.edu\.cn

###### Abstract

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 proposeCross\-VariableLoss\(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\.

## 1Introduction

Multivariate time series forecasting \(MTSF\) is a fundamental task across diverse domains, where the core challenge lies in accurately capturing complex synergistic relationships among variables\. To model these temporal dynamics, recent pioneering works\(Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1); Zhao and Shen,[2024](https://arxiv.org/html/2608.05742#bib.bib2); Qiu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib3); Hu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib4); Ding et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib5)\)have predominantly focused their efforts on extracting representations from historical input features, a paradigm known as input\-side channel dependence modeling\. In stark contrast, the intrinsic correlations among future predictions at the output\-side have been largely overlooked\. Current mainstream Direct Forecasting \(DF\) paradigms\(Wang et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib6)\)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\. This asymmetry heavily prioritizes the input while neglecting the output\. Consequently, it not only severs the physical co\-evolution among future variables but also induces severe structural errors\.

On one hand, blindly decoupling prediction targets directly violates the intrinsic physical and operational laws of multivariate systems\. In the real world, time series are rarely generated independently; rather, they are governed by shared system dynamics and exhibit strong co\-evolutionary characteristics\. For instance, in quantitative financial markets\(Hu et al\.,[2025b](https://arxiv.org/html/2608.05742#bib.bib7)\), the returns of correlated assets often display highly synchronous surges, plunges, or hedging effects when subjected to macroeconomic shocks\(Hu et al\.,[2025c](https://arxiv.org/html/2608.05742#bib.bib8)\)\. Similarly, in traffic flow forecasting, the traffic volumes recorded by hundreds of sensors within an urban network are driven by shared meteorological conditions \(e\.g\., sudden rainstorms causing widespread congestion\), and synchronized daily commuting patterns\. This results in a profound co\-evolution of future traffic dynamics, as illustrated in Figure[5](https://arxiv.org/html/2608.05742#A2.F5)\. If a model isolates these variables at the output stage, the predicted curves for individual variables might appear numerically well\-fitted\. However, when combined into a joint multivariate forecast, they frequently produce contradictory results that violate true physical constraints and market rules\.

Furthermore, overlooking these synergistic relationships induces critical structural errors at the optimization level\. Modern MTSF models rely heavily on point\-wise loss functions, such as Mean Squared Error \(MSE\)\(Wang et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib9); Qiu et al\.,[2026](https://arxiv.org/html/2608.05742#bib.bib10)\), for optimization\. From the perspective of probabilistic inference, this point\-wise optimization is 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 an unbridgeable theoretical divide, termed an Objective Gap, between existing standard loss functions and the true joint spatio\-temporal distribution of the data\. Consequently, even if two models achieve remarkably similar point\-wise MSE, their cross\-variable covariance structures might differ drastically, completely failing to capture the joint geometric topology of the original data\.

To overcome these dual physical and mathematical deficiencies, the optimization objective must transcend the singular pursuit of absolute point\-wise accuracy and explicitly constrain the relative structural consistency among variables\. This implies that we must capture not only the concurrent interactions \(synchronous effects\) between variables at the same time step but also the lagged propagation of system perturbations as time evolves \(asynchronous effects\)\. To elegantly unify these two phenomena, we introduce a graph\-structured perspective, as graphs are inherently suited for depicting complex, non\-uniform pair\-wise interactions among multivariate data\. Specifically, we propose a novel graph\-based structural regularization scheme termed Cross\-Variable Loss \(CvLoss\)\. By constructing a cross\-variable graph, CvLoss connects different variables at identical time steps to capture synchronous concurrent interactions, while also linking variables across different steps to model asynchronous lagged effects\. Essentially, CvLoss imposes anℓ1\\ell\_\{1\}\-norm Graph Total Variation \(Graph TV\) regularization on the prediction residual field, compelling the model to strictly adhere to the underlying relative structural consistency while minimizing absolute errors\. Furthermore, as a plug\-and\-play regularizer, CvLoss seamlessly integrates into various mainstream forecasting architectures\. Our introduced edge\-sampling mechanism significantly reduces computational complexity during training, and crucially, the regularizer imposes zero overhead during inference, ensuring testing speed remains completely unaffected while maintaining extremely high predictive accuracy\.

Our main contributions are summarized as follows:

- •We reveal a critical limitation in current MTSF paradigms, demonstrating that standard point\-wise learning objectives neglect inherent inter\-channel error dependencies among future predictions\.
- •We propose a novel Cross\-Variable Loss acting as a structural regularizer, explicitly constraining the inter\-channel relationships of predicted values via a residual consistency formulation\.
- •Extensive experiments demonstrate our proposed objective consistently improves various state\-of\-the\-art models, achieving superior predictive performance with affordable computational overhead\.

## 2Preliminaries and Related Work

### 2\.1Problem definition

In this paper, we focus on the MTSF problem\. Throughout, italic capitals denote matrices \(e\.g\.,XX,YY,ZZ,EE\), bold lowercase letters denote vectors \(e\.g\.,𝐞\\mathbf\{e\}\), calligraphic capitals denote sets \(e\.g\.,ℰ\\mathcal\{E\}\), and lowercase letters denote scalars and indices \(e\.g\.,tt,dd,α\\alpha\); the dimension constantsHH,TT,DD,PP,LLandNNare written as capitals by convention\. Given a multivariate time\-series dataset withDDvariables, the input historical sequence is defined asX∈ℝH×DX\\in\\mathbb\{R\}^\{H\\times D\}, and the ground\-truth label sequence is defined asY∈ℝT×DY\\in\\mathbb\{R\}^\{T\\times D\}, whereHHis the lookback window andTTis the forecast horizon\. Modern time series forecasting models are primarily trained under a multitask learning manner, known as the direct forecasting \(DF\) paradigm\. Therefore, the target is to learn a multi\-output predictorfθ:ℝH×D→ℝT×Df\_\{\\theta\}:\\mathbb\{R\}^\{H\\times D\}\\rightarrow\\mathbb\{R\}^\{T\\times D\}, parameterised by the trainable weightsθ\\thetaof the forecasting backbone, that directly producesTT\-step forecasts for allDDvariables simultaneously, yielding the predictionY^=fθ​\(X\)\\hat\{Y\}=f\_\{\\theta\}\(X\)\. Unlike current paradigms that decouple multivariate outputs into isolated scalar predictions, our formulation explicitly maintains the joint structure to capture cross\-channel dependencies and shared temporal dynamics\.

### 2\.2Learning objectives in time\-series forecasting

Recent learning objectives extend point\-wise MSE to capture temporal dynamics through shape alignment\(Liu et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib11); Le Guen and Thome,[2019](https://arxiv.org/html/2608.05742#bib.bib12); Cuturi and Blondel,[2017](https://arxiv.org/html/2608.05742#bib.bib13)\), likelihood maximization\(Wang et al\.,[2026a](https://arxiv.org/html/2608.05742#bib.bib14),[2025b](https://arxiv.org/html/2608.05742#bib.bib15); Lange et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib16); Wang et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib6)\), distribution balancing\(Pan et al\.,[2026](https://arxiv.org/html/2608.05742#bib.bib17); Wang et al\.,[2026b](https://arxiv.org/html/2608.05742#bib.bib18)\), and intra\-sequence modeling via decomposition\(Qiu et al\.,[2025b](https://arxiv.org/html/2608.05742#bib.bib19)\)or local step alignment\(Xiong et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib20); Kudrat et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib21)\)\. Despite these advancements, existing approaches fundamentally treat multivariate forecasting as a collection of isolated scalar predictions\. By optimizing exclusively for intra\-channel temporal structures, they completely neglect explicit constraints on cross\-variable dependencies\. Consequently, these methods fail to capture the joint geometric structure and correlated disturbances among multivariate outputs\.

### 2\.3Cross\-variable Dependencies in time\-series forecasting

Effectively modeling cross\-variable dependencies is fundamental for MTSF\(Wang et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib9); Liang et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib22)\)\. Existing paradigms typically employ Channel Dependence \(CD\) strategies to explicitly capture inter\-channel correlations, as seen in representative architectures like iTransformer\(Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\), TSMixer\(Chen et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib23)\), Crossformer\(Zhang and Yan,[2023](https://arxiv.org/html/2608.05742#bib.bib24)\), and TimesNet\(Wu et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib25)\)\. However, CD models often overfit to spurious correlations and struggle with noise interference from irrelevant variables\. To address this, recent advancements introduce Channel Partiality \(CP\) to selectively filter dependencies\. For instance, DUET\(Qiu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib3)\)proposes a Channel Clustering Module to filter noisy channels via frequency\-domain soft clustering\. Similarly, TimeFilter\(Hu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib4)\)constructs a patch\-specific spatial\-temporal graph with dynamic routing to retain essential correlations while discarding irrelevant noise\. Other approaches like MTGNN\(Wu et al\.,[2020](https://arxiv.org/html/2608.05742#bib.bib26)\)and MCformer\(Han et al\.,[2024a](https://arxiv.org/html/2608.05742#bib.bib27)\)also explore sparse interactions to balance accuracy and efficiency\. Despite these significant improvements, current models predominantly focus on extracting representations from historical observations\(Zhao et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib28); Wei et al\.,[2026](https://arxiv.org/html/2608.05742#bib.bib29)\)\. The complex structural geometry and co\-evolving cross\-variable dependencies among future predictions remain fundamentally underexplored\.

![Refer to caption](https://arxiv.org/html/2608.05742v1/x1.png)\(a\)Ground\-truth structure\.![Refer to caption](https://arxiv.org/html/2608.05742v1/x2.png)\(b\)Baseline structure\.![Refer to caption](https://arxiv.org/html/2608.05742v1/x3.png)\(c\)Improvement by CvLoss\.![Refer to caption](https://arxiv.org/html/2608.05742v1/x4.png)\(d\)ECL Snapshot\.Figure 1:Motivating example on ECL with input lengthH=96H=96and forecast horizonT=96T=96\. \(a,b\) Cross\-variable correlation structures of the ground\-truth future series and of the series predicted by a plain iTransformer\(Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\)\. \(c\) Entries on which the same backbone trained with CvLoss attains a smaller structural error than the plain baseline, covering84\.83%84\.83\\%of all cross\-variable entries\. \(d\) Forecast snapshot on a representative variable\. The protocol used to compute \(a\)–\(c\) is defined in Appendix[B\.3](https://arxiv.org/html/2608.05742#A2.SS3)\.

## 3Proposed Method

### 3\.1Motivation

Modern time series forecasting models are primarily trained under the multitask learning manner, known as the direct forecasting \(DF\) paradigm\(Wang et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib6)\)\. Given an input historyX∈ℝH×DX\\in\\mathbb\{R\}^\{H\\times D\}, a multi\-output predictorfθ:ℝH×D→ℝT×Df\_\{\\theta\}:\\mathbb\{R\}^\{H\\times D\}\\rightarrow\\mathbb\{R\}^\{T\\times D\}directly producesTT\-step forecasts forDDvariables:Y^=fθ​\(X\)\.\\hat\{Y\}=f\_\{\\theta\}\(X\)\.The standard training objective is the point\-wise MSE:

ℒdf=1T​∑t=1T1D​‖Yt,:−Y^t,:‖22\\mathcal\{L\}\_\{\\mathrm\{df\}\}=\\frac\{1\}\{T\}\\sum\_\{t=1\}^\{T\}\\frac\{1\}\{D\}\\\|Y\_\{t,:\}\-\\hat\{Y\}\_\{t,:\}\\\|\_\{2\}^\{2\}\(1\)
This objective is simple and effective, but it only penalizes prediction errors at individual time\-variable entries\. Although modern forecasting architectures may implicitly capture cross\-variable interactions through shared parameters, such dependencies are not explicitly reflected in the objective in Eq\. \([1](https://arxiv.org/html/2608.05742#S3.E1)\), which does not constrain the*relative error structure*across variables\. As a result, two models with similar point\-wise MSE may still exhibit substantially different cross\-variable correlated patterns\.

Figure[1\(a\)](https://arxiv.org/html/2608.05742#S2.F1.sf1)makes this concrete on ECL\. Panels \(a\) and \(b\) show the cross\-variable correlation matrix of the ground\-truth future window and of the window predicted by a plain iTransformer trained with MSE: the predicted matrix is visibly flatter and misplaces much of the off\-diagonal structure, even though the point\-wise error of that model is competitive\. Panel \(c\) marks the entries on which the identical backbone trained with our objective recovers the correlation more accurately, covering84\.83%84\.83\\%of all cross\-variable entries, and panel \(d\) shows the corresponding forecast snapshot\. The correlation structure of the future window is therefore something a point\-wise objective leaves largely unconstrained, and something that can be improved without changing the architecture\. This observation motivates the analysis below\.

From a broader perspective, multivariate forecasting is not merely a collection of isolated scalar prediction problems, since different variables often evolve with shared dynamics and correlated disturbances\. This view is also consistent with multi\-task learning, where explicit coupling across related targets can improve generalization\(Evgeniou et al\.,[2005](https://arxiv.org/html/2608.05742#bib.bib30)\)\. Consequently, relying purely on point\-wise estimation may fail to fully exploit the shared structure across variables\.

Motivated by these perspectives, we seek to augment the vanilla DF objective with an explicit structural regularizer that captures dependencies among variables\. Rather than replacing point\-wise supervision, this additional term is intended to complement it by encouraging consistency in the residual structure across related variables\. To clarify why such a structural term may be beneficial, we first examine the standard DF objective from a probabilistic perspective\.

Let𝐞=vec​\(Y^−Y\)∈ℝT​D\\mathbf\{e\}=\\mathrm\{vec\}\(\\hat\{Y\}\-Y\)\\in\\mathbb\{R\}^\{TD\}denote the flattened prediction error vector over the entire forecast horizon and all variables\. Suppose the residuals follow a zero\-mean multivariate Gaussian distribution𝐞∼𝒩​\(𝟎,𝚺S​T\)\\mathbf\{e\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\mathbf\{\\Sigma\}\_\{ST\}\), where𝚺S​T∈ℝT​D×T​D\\mathbf\{\\Sigma\}\_\{ST\}\\in\\mathbb\{R\}^\{TD\\times TD\}is the spatiotemporal covariance matrix that captures both cross\-variable and cross\-temporal dependencies\.

###### Theorem 3\.1\(Objective Gap with Spatiotemporal Dependencies\)

Ignoring additive constants independent of the predictions, the scaled negative log\-likelihood \(NLL\) under the true covariance𝚺S​T\\mathbf\{\\Sigma\}\_\{ST\}isℒnll=12​T​D​𝐞⊤​𝚺S​T−1​𝐞\.\\mathcal\{L\}\_\{\\mathrm\{nll\}\}=\\frac\{1\}\{2TD\}\\mathbf\{e\}^\{\\top\}\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}\\mathbf\{e\}\.On the other hand, the standard DF objective in Eq\. \([1](https://arxiv.org/html/2608.05742#S3.E1)\) can be written asℒdf=1T​D​𝐞⊤​𝐞\.\\mathcal\{L\}\_\{\\mathrm\{df\}\}=\\frac\{1\}\{TD\}\\mathbf\{e\}^\{\\top\}\\mathbf\{e\}\.Therefore, compared with the isotropic Gaussian NLL under varianceσ2\\sigma^\{2\}, the discrepancy between the two objectives is

Δ=ℒnll−12​σ2​ℒdf=12​T​D​𝐞⊤​\(𝚺S​T−1−1σ2​𝐈\)​𝐞,\\Delta=\\mathcal\{L\}\_\{\\mathrm\{nll\}\}\-\\frac\{1\}\{2\\sigma^\{2\}\}\\mathcal\{L\}\_\{\\mathrm\{df\}\}=\\frac\{1\}\{2TD\}\\mathbf\{e\}^\{\\top\}\\left\(\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}\-\\frac\{1\}\{\\sigma^\{2\}\}\\mathbf\{I\}\\right\)\\mathbf\{e\},\(2\)where𝐈∈ℝT​D×T​D\\mathbf\{I\}\\in\\mathbb\{R\}^\{TD\\times TD\}is the identity matrix\. In real\-world multivariate time series, the presence of cross\-variable and lagged dependencies ensures that the true precision matrix𝚺S​T−1\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}contains non\-zero off\-diagonal elements\. Consequently,𝚺S​T−1≠1σ2​𝐈\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}\\neq\\frac\{1\}\{\\sigma^\{2\}\}\\mathbf\{I\}, meaning a strictly non\-zero objective gap \(Δ≠0\\Delta\\neq 0\) inherently exists because the standard point\-wise DF objective neglects these underlying structural correlations\.

Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)shows that the standard point\-wise DF objective coincides with the isotropic Gaussian NLL only when the residual covariance is spherical, i\.e\.,𝚺S​T=σ2​𝐈\\mathbf\{\\Sigma\}\_\{ST\}=\\sigma^\{2\}\\mathbf\{I\}\. When substantial cross\-variable dependencies are present, potentially together with lagged temporal interactions, minimizing point\-wise MSE alone does not explicitly account for the underlying residual dependency structure\. However, in real\-world systems, dynamic physical processes often induce lagged effects—such as a disturbance in one variable propagating to another variable after several time steps\. Consequently, off\-diagonal elements in𝚺S​T\\mathbf\{\\Sigma\}\_\{ST\}, which represent cross\-variable and cross\-temporal correlations, are significantly non\-zero\. As a result, point\-wise MSE may fail to reflect the true spatiotemporal error geometry and may not fully exploit correlated residual structure during training\.

##### Scope of the theorem\.

We state precisely what Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)does and does not establish, since the distinction governs the rest of the paper\. It establishes that*some*non\-diagonal term is necessary: whenever the residual precision matrix has non\-zero off\-diagonal entries, the point\-wise objective leaves a strictly non\-zero gapΔ\\Delta\. It does*not*identify which non\-diagonal term should be used, and in particular it does not derive the specific objective we adopt in Section[3\.2](https://arxiv.org/html/2608.05742#S3.SS2)\. Estimating𝚺S​T−1\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}itself is not a viable route at this scale\. On ECL withT=720T=720andD=321D=321that matrix has more than5×10105\\times 10^\{10\}entries, so the estimation problem is ill\-posed, and we do not attempt it\. What the theorem licenses is the weaker and sufficient claim that the objective should be*relation\-aware*: it must couple residuals that belong to different variables instead of scoring them independently\. Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)below makes the resulting family explicit, and the empirical choice of norm within that family is examined in Appendix[C\.4](https://arxiv.org/html/2608.05742#A3.SS4)\.

###### Proposition 3\.2\(MSE with a squared graph penalty is a structured Gaussian NLL\)

Letℰ\\mathcal\{E\}be a cross\-variable edge set over theNNpatch nodes, letA~=A⊗IL\\tilde\{A\}=A\\otimes I\_\{L\}be the lifted incidence matrix of Section[3\.2](https://arxiv.org/html/2608.05742#S3.SS2), and letLap​\(ℰ\)=A~⊤​A~\\mathrm\{Lap\}\(\\mathcal\{E\}\)=\\tilde\{A\}^\{\\top\}\\tilde\{A\}be the corresponding graph Laplacian\. Define the structured precision family

𝚺−1​\(λ\)=1σ2​\(𝐈\+λ​Lap​\(ℰ\)\),λ≥0\.\\mathbf\{\\Sigma\}^\{\-1\}\(\\lambda\)=\\frac\{1\}\{\\sigma^\{2\}\}\\bigl\(\\mathbf\{I\}\+\\lambda\\,\\mathrm\{Lap\}\(\\mathcal\{E\}\)\\bigr\),\\qquad\\lambda\\geq 0\.\(3\)Then, writing𝐫=vec​\(E\)\\mathbf\{r\}=\\mathrm\{vec\}\(E\)for the flattened residual field, the following identity holds for every𝐫\\mathbf\{r\}:

𝐫⊤​𝚺−1​\(λ\)​𝐫=1σ2​\(‖𝐫‖22\+λ​∑\(i,j\)∈ℰ‖𝐞i−𝐞j‖22\)\.\\mathbf\{r\}^\{\\top\}\\mathbf\{\\Sigma\}^\{\-1\}\(\\lambda\)\\,\\mathbf\{r\}=\\frac\{1\}\{\\sigma^\{2\}\}\\Bigl\(\\\|\\mathbf\{r\}\\\|\_\{2\}^\{2\}\+\\lambda\\\!\\\!\\sum\_\{\(i,j\)\\in\\mathcal\{E\}\}\\\!\\\!\\\|\\mathbf\{e\}\_\{i\}\-\\mathbf\{e\}\_\{j\}\\\|\_\{2\}^\{2\}\\Bigr\)\.\(4\)

Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)is proved in Appendix[A\.2](https://arxiv.org/html/2608.05742#A1.SS2)\. It identifies exactly the object that Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)shows MSE to be missing: point\-wise MSE plus a*squared*graph total\-variation penalty overℰ\\mathcal\{E\}is the negative log\-likelihood of a Gaussian whose precision matrix has supportℰ\\mathcal\{E\}off the diagonal, soℰ\\mathcal\{E\}is not an arbitrary analogy but the conditional\-independence structure of the residual field\. Settingλ=0\\lambda=0recovers the spherical case of Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1), in whichΔ=0\\Delta=0\. This also settles the status of the default topology: becauseℰ\\mathcal\{E\}is the support of the off\-diagonal entries, a complete cross\-variable graph imposes*no*zero constraint on any of them and is therefore the weakest assumption in the family, whereas the empty graph, i\.e\. plain MSE, is the member that Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)rules out\.

### 3\.2Predicting Multivariate Time Series with Cross Variable Loss Constraints

#### 3\.2\.1Definition

##### Patch nodes\.

LetY∈ℝT×DY\\in\\mathbb\{R\}^\{T\\times D\}be the target andY^=fθ​\(X\)∈ℝT×D\\hat\{Y\}=f\_\{\\theta\}\(X\)\\in\\mathbb\{R\}^\{T\\times D\}be the prediction\. We split the time axis intoPPnon\-overlapping patches of lengthLL\(Nie et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib31); Ding et al\.,[2026](https://arxiv.org/html/2608.05742#bib.bib32)\), so thatT=P×LT=P\\times L\. Each node is defined as a pairv=\(p,d\)v=\(p,d\), wherep∈\{1,…,P\}p\\in\\\{1,\\dots,P\\\}indexes the patch andd∈\{1,…,D\}d\\in\\\{1,\\dots,D\\\}indexes the variable\. Let𝐳v∈ℝL\\mathbf\{z\}\_\{v\}\\in\\mathbb\{R\}^\{L\}and𝐳^v∈ℝL\\hat\{\\mathbf\{z\}\}\_\{v\}\\in\\mathbb\{R\}^\{L\}denote the ground\-truth and predicted patch vectors at nodevv, respectively\. Stacking all nodes gives

Z∈ℝN×L,Z^∈ℝN×L,N=P×D\.Z\\in\\mathbb\{R\}^\{N\\times L\},\\qquad\\hat\{Z\}\\in\\mathbb\{R\}^\{N\\times L\},\\qquad N=P\\times D\.\(5\)We define the node\-wise residuals as

𝐞v=𝐳^v−𝐳v∈ℝL,\\mathbf\{e\}\_\{v\}=\\hat\{\\mathbf\{z\}\}\_\{v\}\-\\mathbf\{z\}\_\{v\}\\in\\mathbb\{R\}^\{L\},\(6\)and stack them into the residual matrixE∈ℝN×LE\\in\\mathbb\{R\}^\{N\\times L\}, whose vectorisationvec​\(E\)∈ℝN​L\\mathrm\{vec\}\(E\)\\in\\mathbb\{R\}^\{NL\}is a permutation of the flattened error vector𝐞\\mathbf\{e\}of Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1), sinceN​L=T​DNL=TD\. ForL=1L=1, patch nodes recover point nodes\.

##### Cross\-variable edges\.

To capture the spatiotemporal dependencies across different series, we introduce an edge setℰ\\mathcal\{E\}to define the cross\-variable graph topology\. Guided by Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1), rather than explicitly estimating the intractable precision matrix, we utilizeℰ\\mathcal\{E\}as a practical surrogate structure to regularize residual discrepancies\. This formulation is highly flexible and can easily accommodate different structural priors\. Specifically, we design three variants for constructingℰ\\mathcal\{E\}: \(1\)*Synchronous cross\-variable edges*, which strictly capture concurrent interactions at the same time step \(i\.e\., connecting identical patch indices\); \(2\)*Asynchronous cross\-variable edges*, which exclusively model lagged dependencies across different time steps; and \(3\)*Fully connected cross\-variable edges*, which include both synchronous and lagged interactions\. We adopt the fully connected graph as our default setting to maximize the coverage of potential interactions without relying on domain\-specific prior knowledge\. We evaluate these topological designs in the ablation studies in Section[4\.5](https://arxiv.org/html/2608.05742#S4.SS5)\.

##### Difference operator\.

Fix an arbitrary orientation for each undirected edge and letA∈ℝ\|ℰ\|×NA\\in\\mathbb\{R\}^\{\|\\mathcal\{E\}\|\\times N\}be the incidence matrix: for an edgek=\(i,j\)k=\(i,j\), thekk\-th row contains\+1\+1at positionii,−1\-1at positionjj, and0elsewhere\. Since each node is anLL\-dimensional patch, we lift the incidence operator as

A~:=A⊗IL\.\\tilde\{A\}:=A\\otimes I\_\{L\}\.\(7\)Letvec​\(E\)∈ℝN​L\\mathrm\{vec\}\(E\)\\in\\mathbb\{R\}^\{NL\}denote the vectorization ofEE\. ThenA~​vec​\(E\)∈ℝ\|ℰ\|​L\\tilde\{A\}\\,\\mathrm\{vec\}\(E\)\\in\\mathbb\{R\}^\{\|\\mathcal\{E\}\|L\}stacks the residual differences𝐞i−𝐞j∈ℝL\\mathbf\{e\}\_\{i\}\-\\mathbf\{e\}\_\{j\}\\in\\mathbb\{R\}^\{L\}over all edges\.

#### 3\.2\.2Cross\-Variable Loss

Building upon the defined graph topology, we now formalize our objective\. First, by reshaping the predictions into patch matrices, the standard DF loss can be equivalently rewritten as

ℒdf=1N×L​‖Z^−Z‖F2=1T×D​‖Y^−Y‖F2\.\\mathcal\{L\}\_\{\\mathrm\{df\}\}=\\frac\{1\}\{N\\times L\}\\\|\\hat\{Z\}\-Z\\\|\_\{F\}^\{2\}=\\frac\{1\}\{T\\times D\}\\\|\\hat\{Y\}\-Y\\\|\_\{F\}^\{2\}\.\(8\)To capture the relative structural geometry between variables, for each edge\(i,j\)∈ℰ\(i,j\)\\in\\mathcal\{E\}, we define the predicted and true patch differences as

Δi​jpred=𝐳^i−𝐳^j,Δi​jtrue=𝐳i−𝐳j\.\\Delta^\{\\mathrm\{pred\}\}\_\{ij\}=\\hat\{\\mathbf\{z\}\}\_\{i\}\-\\hat\{\\mathbf\{z\}\}\_\{j\},\\qquad\\Delta^\{\\mathrm\{true\}\}\_\{ij\}=\\mathbf\{z\}\_\{i\}\-\\mathbf\{z\}\_\{j\}\.\(9\)We then define CvLoss as the mean absolute discrepancy of these structural differences:

ℒcv:=1\|ℰ\|×L​∑\(i,j\)∈ℰ‖Δi​jpred−Δi​jtrue‖1\.\\mathcal\{L\}\_\{\\mathrm\{cv\}\}:=\\frac\{1\}\{\|\\mathcal\{E\}\|\\times L\}\\sum\_\{\(i,j\)\\in\\mathcal\{E\}\}\\bigl\\\|\\Delta^\{\\mathrm\{pred\}\}\_\{ij\}\-\\Delta^\{\\mathrm\{true\}\}\_\{ij\}\\bigr\\\|\_\{1\}\.\(10\)
Using the node\-wise residual𝐞i=𝐳^i−𝐳i\\mathbf\{e\}\_\{i\}=\\hat\{\\mathbf\{z\}\}\_\{i\}\-\\mathbf\{z\}\_\{i\}, the structural discrepancy on each edge simplifies neatly:

Δi​jpred−Δi​jtrue=\(𝐳^i−𝐳^j\)−\(𝐳i−𝐳j\)=𝐞i−𝐞j\.\\Delta^\{\\mathrm\{pred\}\}\_\{ij\}\-\\Delta^\{\\mathrm\{true\}\}\_\{ij\}=\(\\hat\{\\mathbf\{z\}\}\_\{i\}\-\\hat\{\\mathbf\{z\}\}\_\{j\}\)\-\(\\mathbf\{z\}\_\{i\}\-\\mathbf\{z\}\_\{j\}\)=\\mathbf\{e\}\_\{i\}\-\\mathbf\{e\}\_\{j\}\.\(11\)Therefore, CvLoss admits a compact graph\-theoretic form via the incidence matrixA~\\tilde\{A\}:

ℒcv=1\|ℰ\|×L​∑\(i,j\)∈ℰ‖𝐞i−𝐞j‖1=1\|ℰ\|×L​‖A~​vec​\(E\)‖1\.\\mathcal\{L\}\_\{\\mathrm\{cv\}\}=\\frac\{1\}\{\|\\mathcal\{E\}\|\\times L\}\\sum\_\{\(i,j\)\\in\\mathcal\{E\}\}\\\|\\mathbf\{e\}\_\{i\}\-\\mathbf\{e\}\_\{j\}\\\|\_\{1\}=\\frac\{1\}\{\|\\mathcal\{E\}\|\\times L\}\\\|\\tilde\{A\}\\,\\mathrm\{vec\}\(E\)\\\|\_\{1\}\.\(12\)This computes the total variation of prediction errors across the graph\.

![Refer to caption](https://arxiv.org/html/2608.05742v1/x5.png)Figure 2:Overview of the workflow CvLoss\.
#### 3\.2\.3Interpretation

The role of CvLoss is to explicitly regularize the cross\-variable residual structure that is not directly controlled by the point\-wise DF loss\. In particular,ℒdf\\mathcal\{L\}\_\{\\mathrm\{df\}\}encourages each residual𝐞v\\mathbf\{e\}\_\{v\}to be small in absolute magnitude, whileℒcv\\mathcal\{L\}\_\{\\mathrm\{cv\}\}encourages connected nodes to have compatible residual patterns\. This induces an explicit coupling across variables at the objective level\.

From the graph perspective,ℒcv\\mathcal\{L\}\_\{\\mathrm\{cv\}\}is anℓ1\\ell\_\{1\}\-type graph total variation regularizer on the residual field\. Unlike point\-wise MSE, which evaluates each entry independently, CvLoss constrains the edge\-wise discrepancies of prediction errors across variables\. This encourages the model to exploit shared temporal structure and improves coordination among related variables\.

We use theℓ1\\ell\_\{1\}norm because it is less sensitive to a few large edge\-wise mismatches caused by localized anomalies or variable\-specific shocks, while still promoting graph\-level residual consistency\. Replacing the squared penalty of Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)byℓ1\\ell\_\{1\}yields a pairwise Gibbs field over the same graph rather than a Gaussian one, so theℓ1\\ell\_\{1\}form is a task\-specific empirical choice rather than a consequence of the proposition; Appendix[C\.4](https://arxiv.org/html/2608.05742#A3.SS4)compares the two norms directly under an identical protocol\.

### 3\.3Model Implementation

Our framework augments vanilla DF with the CvLoss constraint, complementing point\-wise supervision with cross\-variable regularization\.

As shown in[Figure 2](https://arxiv.org/html/2608.05742#S3.F2), the normalized historyXXis fed into the model to produceY^=fθ​\(X\)\\hat\{Y\}=f\_\{\\theta\}\(X\)\. We compute the point\-wise errorℒdf\\mathcal\{L\}\_\{\\mathrm\{df\}\}, map forecasts and labels into patch nodes, and compute the structural errorℒcv\\mathcal\{L\}\_\{\\mathrm\{cv\}\}overℰ\\mathcal\{E\}\. Since the incidence operations are linear,ℒcv\\mathcal\{L\}\_\{\\mathrm\{cv\}\}is optimized end\-to\-end\.

Finally, the independent DF error and the structural cross\-variable error are fused\. A weighting parameterα∈\(0,1\)\\alpha\\in\(0,1\)controls the relative contribution of each objective:

ℒα:=\(1−α\)⋅ℒdf\+α⋅ℒcv\.\\mathcal\{L\}\_\{\\alpha\}:=\(1\-\\alpha\)\\cdot\\mathcal\{L\}\_\{\\mathrm\{df\}\}\+\\alpha\\cdot\\mathcal\{L\}\_\{\\mathrm\{cv\}\}\.\(13\)Throughout the paperα\\alphais a*fixed*scalar selected on the validation split, together with the patch lengthLLthe only quantity we tune; it is never learned during training and never varies across batches, variables or edges\. Appendix[B\.2](https://arxiv.org/html/2608.05742#A2.SS2.SSS0.Px1)records the selection protocol, and Appendix[D\.2](https://arxiv.org/html/2608.05742#A4.SS2)reports a separate diagnostic study whose adaptive coefficients are not values ofα\\alpha\.

By sharing statistical strength across the multivariate series, our framework complements the DF paradigm with explicit cross\-variable regularization while preserving its practical advantages, including efficient inference and multi\-task learning capabilities\. It is worth noting that the standardℒdf\\mathcal\{L\}\_\{\\mathrm\{df\}\}term serves as a mathematically necessary absolute anchor\. It prevents uniform drift and trivial translation\-invariant solutions \(where𝐞i=𝐞j=𝐜≠𝟎\\mathbf\{e\}\_\{i\}=\\mathbf\{e\}\_\{j\}=\\mathbf\{c\}\\neq\\mathbf\{0\}\) even when the structural penalty dominates \(i\.e\.,α→1\\alpha\\to 1\)\. Furthermore, our proposed regularization scheme is entirely model\-agnostic, seamlessly compatible with various forecasting architecturesfθf\_\{\\theta\}\(e\.g\., Transformers, CNNs, and MLPs\)\.

## 4Experiments

To demonstrate the efficacy of CvLoss, the following aspects deserve empirical investigation:

- •Performance:Does CvLoss perform well?In Section[4\.2](https://arxiv.org/html/2608.05742#S4.SS2), we benchmark CvLoss against state\-of\-the\-art baselines, and in Section[4\.3](https://arxiv.org/html/2608.05742#S4.SS3), we compare it with alternative learning objectives\.
- •Gain:Why does it work?In Section[4\.5](https://arxiv.org/html/2608.05742#S4.SS5), we perform an ablative study, dissecting the individual components of CvLoss and clarifying their contributions to forecast accuracy\.
- •Generality:Does it support other models?In Section[4\.4](https://arxiv.org/html/2608.05742#S4.SS4), we examine its compatibility with various architectures and datasets, with further results in Appendix[C\.5](https://arxiv.org/html/2608.05742#A3.SS5)\.
- •Sensitivity:Is it sensitive to hyperparameters?In Section[4\.6](https://arxiv.org/html/2608.05742#S4.SS6), we analyze the sensitivity of CvLoss to the hyperparameterα\\alpha, showing stable performance across a broad parameter range\.
- •Efficiency:What is the computational cost of it?In Appendix[C\.10](https://arxiv.org/html/2608.05742#A3.SS10), we evaluate the running cost of CvLoss across different scenarios\.

### 4\.1Setup

##### Datasets\.

Following\(Wu et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib25)\), we evaluate on standard benchmarks: four ETT subsets, ECL, Traffic, Weather, and PEMS\(Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\)\. All datasets are chronologically split into training, validation, and test sets\. Statistics are in Appendix[B\.1](https://arxiv.org/html/2608.05742#A2.SS1)\.

##### Baselines\.

We compare CvLoss to a range of competitive baselines, categorized as: \(1\) Transformer\-based models: PatchTST\(Nie et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib31)\), iTransformer\(Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\), TQNet\(Lin et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib33)\), Crossformer\(Zhang and Yan,[2023](https://arxiv.org/html/2608.05742#bib.bib24)\), FEDformer\(Zhou et al\.,[2022](https://arxiv.org/html/2608.05742#bib.bib34)\)and Leddam\(Yu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib35)\); \(\(2\) Non\-transformer models: DLinear\(Zeng et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib36)\), SOFTS\(Han et al\.,[2024b](https://arxiv.org/html/2608.05742#bib.bib37)\), TimesNet\(Wu et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib25)\)and TimeFilter\(Hu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib4)\)\.

##### Implementation\.

Baseline implementations closely follow the official codebase from TimeFilter\(Hu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib4)\)and TQNet\(Lin et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib33)\)\. To ensure fair comparison, the drop\-last trick is disabled for all models, as recommended in\(Qiu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib38)\)\. All models are trained with the Adam optimizer\(Kingma and Ba,[2015](https://arxiv.org/html/2608.05742#bib.bib39)\)\. When integrating CvLoss into a baseline forecast model, we retain all hyperparameters from the public benchmarks\(Hu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib4); Lin et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib33)\), only tuningα\\alphaand the patch size\. Experiments are run on Intel\(R\) Xeon\(R\) Platinum 8470Q with 8 NVIDIA RTX 5090 GPUs\. Further implementation details are provided in Appendix[B](https://arxiv.org/html/2608.05742#A2)\.

### 4\.2Overall performance

Table 1:Overall forecasting performance\.ModelsCvLoss

TimeFilter

TQNet

iTransformer

Leddam

SOFTS

PatchTST

Crossformer

TimesNet

DLinear

FEDformer

\(Ours\)

\(2025\)

\(2025\)

\(2024\)

\(2024\)

\(2023\)

\(2023\)

\(2023\)

\(2023\)

\(2023\)

\(2022\)

Metrics

MSE

MAE

MSE

MAE

MSE

MAE

MSE

MAE

MSE

MAE

MSE

MAE

MSE

MAE

MSE

MAE

MSE

MAE

MSE

MAE

MSE

MAE

ETTm1

0\.372

0\.381

0\.3770\.3930\.3770\.3930\.407

0\.410

0\.386

0\.397

0\.393

0\.403

0\.387

0\.400

0\.513

0\.496

0\.400

0\.406

0\.404

0\.408

0\.448

0\.452

ETTm2

0\.270

0\.317

0\.2720\.3210\.278

0\.322

0\.288

0\.332

0\.281

0\.325

0\.287

0\.330

0\.281

0\.326

0\.757

0\.610

0\.291

0\.333

0\.354

0\.402

0\.305

0\.349

ETTh1

0\.419

0\.427

0\.4200\.4280\.441

0\.434

0\.454

0\.447

0\.431

0\.429

0\.449

0\.442

0\.469

0\.454

0\.529

0\.522

0\.458

0\.450

0\.461

0\.457

0\.440

0\.460

ETTh2

0\.359

0\.391

0\.3650\.3980\.378

0\.403

0\.383

0\.407

0\.373

0\.399

0\.385

0\.408

0\.387

0\.407

0\.942

0\.684

0\.414

0\.427

0\.563

0\.519

0\.437

0\.449

Weather

0\.236

0\.260

0\.2400\.270

0\.242

0\.2680\.258

0\.279

0\.242

0\.272

0\.255

0\.278

0\.259

0\.281

0\.259

0\.315

0\.259

0\.287

0\.265

0\.315

0\.309

0\.360

ECL

0\.157

0\.252

0\.1590\.2560\.164

0\.258

0\.178

0\.270

0\.169

0\.263

0\.174

0\.264

0\.216

0\.304

0\.244

0\.334

0\.193

0\.295

0\.225

0\.319

0\.214

0\.327

Traffic

0\.407

0\.254

0\.4080\.269

0\.445

0\.276

0\.428

0\.282

0\.467

0\.294

0\.4090\.2670\.555

0\.362

0\.550

0\.304

0\.620

0\.336

0\.625

0\.383

0\.610

0\.376

Solar

0\.2180\.254

0\.228

0\.2620\.197

0\.255

0\.233

0\.262

0\.230

0\.264

0\.229

0\.256

0\.270

0\.307

0\.641

0\.639

0\.301

0\.319

0\.330

0\.401

0\.292

0\.381

PEMS03

0\.079

0\.183

0\.084

0\.191

0\.0800\.1860\.096

0\.204

0\.101

0\.210

0\.087

0\.192

0\.151

0\.265

0\.138

0\.253

0\.119

0\.271

0\.219

0\.295

0\.167

0\.291

PEMS04

0\.078

0\.181

0\.083

0\.186

0\.0810\.1850\.098

0\.207

0\.102

0\.213

0\.091

0\.196

0\.162

0\.273

0\.145

0\.267

0\.109

0\.220

0\.236

0\.350

0\.195

0\.308

PEMS07

0\.063

0\.156

0\.071

0\.170

0\.0650\.1610\.088

0\.190

0\.087

0\.192

0\.075

0\.173

0\.166

0\.270

0\.181

0\.272

0\.106

0\.208

0\.241

0\.343

0\.133

0\.282

PEMS08

0\.080

0\.182

0\.0830\.1860\.105

0\.203

0\.127

0\.212

0\.102

0\.211

0\.114

0\.208

0\.238

0\.289

0\.232

0\.270

0\.150

0\.244

0\.347

0\.421

0\.234

0\.326

- •Note: We fix the input length as 96 for all baselines following\(Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\)\.Boldandunderlineddenote best and second\-best results, respectively\.*Avg*indicates average results over multiple forecast horizons\. CvLoss using TimeFilter\(Hu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib4)\)as backbone\.

Table 2:Backbone\-controlled summary of CvLoss\.Every comparison contrasts a forecasting backbone trained with its standard objective against the*identical*backbone, data split, and training protocol trained with CvLoss, so the loss term is the only variable\. A*cell*is one \(backbone, dataset, metric\) triple evaluated at the dataset\-level average over forecast horizons; W/T/L counts cells on which CvLoss is better, equal, or worse, and the last two columns give the mean relative error reduction over the cells of that block\.- •Note: The seven distinct backbones are TimeFilter, TQNet, PDF, CFPT, iTransformer, DLinear and PatchTST; the twelve distinct datasets are the four ETT subsets, Weather, ECL, Traffic, Solar and PEMS03/04/07/08\. Blocks follow the evaluation protocol of the work each one compares against, so they are overlapping rather than disjoint samples and the total is a summary rather than a pooled statistical test\.*The single loss*is TQNet/ETTh1 MAE in Table[3](https://arxiv.org/html/2608.05742#S4.T3)\(0\.434→0\.4380\.434\\\!\\rightarrow\\\!0\.438\)\.*The two ties*are DLinear/Weather MSE and PatchTST/ECL MSE in Table[19](https://arxiv.org/html/2608.05742#A3.T19), which are equal at the reported three\-decimal precision\. Counting the same evidence at the finer horizon level rather than at the dataset average gives 96/0/0 over the 96 comparisons of Table[10](https://arxiv.org/html/2608.05742#A2.T10)and 31/1/0 over the 32 comparisons of Table[4](https://arxiv.org/html/2608.05742#S4.T4)\.

Table[1](https://arxiv.org/html/2608.05742#S4.T1)compares long\- and short\-term forecasting results\. Adding CvLoss consistently improves accuracy across datasets; on Traffic, for example, MAE decreases by 0\.014\. The gains come from explicit structural regularization over synchronous and asynchronous cross\-variable dependencies\.

##### Backbone\-controlled attribution\.

Table[1](https://arxiv.org/html/2608.05742#S4.T1)answers two questions at once, and they should not be conflated\. The pair of columns CvLoss \(Ours\) and TimeFilter is backbone\-controlled: CvLoss*is*TimeFilter, trained on the same data with the same protocol and differing only in the loss term, and it improves all2424dataset\-level entries \(1212datasets×2\\times\\,2metrics\) with no tie and no loss\. The remaining columns answer a different question, namely whether the backbone we build on is itself competitive, and contribute nothing to the attribution\.

Because a single backbone cannot rule out an interaction between CvLoss and TimeFilter, Table[2](https://arxiv.org/html/2608.05742#S4.T2)consolidates every fixed\-backbone comparison in this paper into one view\. Across114114such comparisons, spanning seven backbones from three architecture families \(Transformer, linear, and hybrid decomposition designs\), twelve datasets and both metrics, CvLoss improves111111, ties22and loses11, with a mean relative error reduction of3\.39%3\.39\\%in MSE and3\.61%3\.61\\%in MAE\. The effect is therefore a property of the objective rather than of any one architecture\. Two qualifications belong with that number\. First, the blocks follow the protocols of the works they compare against and overlap in backbones and datasets, so the total summarises the evidence rather than constituting a pooled significance test; per\-block paired effect sizes with confidence intervals are reported in Appendix[C\.9](https://arxiv.org/html/2608.05742#A3.SS9)\. Second, the reductions are unevenly distributed\. They are large where the data carry strong cross\-variable structure that the backbone leaves unmodelled \(Traffic on iTransformer,−17\.7%\-17\.7\\%MSE; PEMS07 on TimeFilter,−11\.3%\-11\.3\\%\), and small where either condition fails \(ETTh1,−0\.24%\-0\.24\\%MSE on TimeFilter\)\. Section[4\.6](https://arxiv.org/html/2608.05742#S4.SS6)and Appendix[B\.3](https://arxiv.org/html/2608.05742#A2.SS3)examine this dependence directly\.

##### Examples\.

A qualitative comparison between forecasts generated by DF versus CvLoss is presented in Figure[3](https://arxiv.org/html/2608.05742#S4.F3)\. While the model trained with DF captures general periodic trends, it often fails to model accurate amplitudes and intricate temporal dynamics\. For example, on the ETTm2 dataset, it suffers from severe amplitude attenuation, completely missing the prominent ground\-truth peaks\. Similarly, on ECL, it overestimates the high periodic peaks, such as those observed around the 1750th and 250th steps\. In contrast, CvLoss accurately captures these large peak magnitudes and complex fluctuations\.

![Refer to caption](https://arxiv.org/html/2608.05742v1/x6.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x7.png)

\(a\)ETTm2 snapshot\.
![Refer to caption](https://arxiv.org/html/2608.05742v1/x8.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x9.png)

\(b\)ECL snapshot\.

Figure 3:The forecast sequence of DF \(in blue\) and CvLoss \(in red\), with historical lengthH=96\\mathrm\{H\}=96\.

### 4\.3Learning objective comparison

Table 3:Comparable results with other learning objectives for time\-series forecasting\.- •Note:Boldandunderlineddenote best and second\-best results, respectively\. Follow the settings of QDF\(Wang et al\.,[2026a](https://arxiv.org/html/2608.05742#bib.bib14)\)\.

CvLoss is compared with 11 learning objectives from four categories: shape alignment, likelihood maximization, distribution balancing, and decomposition\. Table[3](https://arxiv.org/html/2608.05742#S4.T3)summarizes average performance on four datasets, with full horizon\-level results in Appendix Table[14](https://arxiv.org/html/2608.05742#A3.T14)\.

Shape\-alignment losses \(e\.g\., Soft\-DTW\(Cuturi and Blondel,[2017](https://arxiv.org/html/2608.05742#bib.bib13)\), Koopman\(Lange et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib16)\)\) and likelihood objectives \(e\.g\., Time\-o1\(Wang et al\.,[2025b](https://arxiv.org/html/2608.05742#bib.bib15)\), FreDF\(Wang et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib6)\), QDF\(Wang et al\.,[2026a](https://arxiv.org/html/2608.05742#bib.bib14)\)\) improve over MSE, but still treat multivariate outputs as scalar predictions\. They therefore leave the spatiotemporal objective gap \(Δ≠0\\Delta\\neq 0\) in Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)unaddressed\. By penalizing edge\-wise residual discrepancies, CvLoss ranks first in 20 settings \(Appendix Table[14](https://arxiv.org/html/2608.05742#A3.T14)\) and generalizes across TQNet\(Lin et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib33)\)and PDF\(Dai et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib40)\)\(Appendix Table[13](https://arxiv.org/html/2608.05742#A3.T13)\)\.

### 4\.4Generalization studies

We test CvLoss on TimeFilter, TQNet, CFPT\(Kou et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib41)\), and iTransformer\. As shown in[Figure 4](https://arxiv.org/html/2608.05742#S4.F4), CvLoss improves different architectures and datasets\. On ECL, all four backbones improve, with iTransformer reducing MSE by 8\.6% and MAE by 7\.9%\. On Weather, MSE drops by 1\.3%–2\.3% and MAE by 1\.9%–3\.3%, showing that CvLoss is broadly plug\-and\-play\.

![Refer to caption](https://arxiv.org/html/2608.05742v1/x10.png)\(a\)ECL with MSE
![Refer to caption](https://arxiv.org/html/2608.05742v1/x11.png)\(b\)ECL with MAE
![Refer to caption](https://arxiv.org/html/2608.05742v1/x12.png)\(c\)Weather with MSE
![Refer to caption](https://arxiv.org/html/2608.05742v1/x13.png)\(d\)Weather with MAE

Figure 4:Improvement of CvLoss applied to different forecast models, shown with colored bars for means over forecast lengths \(96, 192, 336, 720\) and error bars for 50% confidence intervals\.Table 4:Ablation study results\.- •Note:Boldandunderlineddenote best and second\-best results, respectively\. “Sync\.” and “Async\.” are abbreviations for synchronous and asynchronous\.

### 4\.5Ablation studies

To investigate the impact of cross\-variable interaction designs within our proposed CVL, we conduct ablation studies on the construction of the edge setℰ\\mathcal\{E\}\. Specifically, we compare the following variants:

- •CvLoss†utilizes a synchronous edge setℰ\\mathcal\{E\}that strictly captures cross\-variable interactions at the same time step\.
- •CvLoss‡utilizes an asynchronous edge setℰ\\mathcal\{E\}that exclusively models lagged cross\-variable interactions across different time steps \(disallowing same\-time interactions\)\.
- •CvLoss is our default setting, whereℰ\\mathcal\{E\}is a fully connected cross\-variable graph that permits both synchronous and lagged interactions\.

As shown in Table[4](https://arxiv.org/html/2608.05742#S4.T4), the default CvLoss generally outperforms both CvLoss†and CvLoss‡in most settings\. This confirms that integrating both synchronous and asynchronous interactions withinℰ\\mathcal\{E\}is crucial for capturing complex dependencies and yields the best overall forecasting performance\.

### 4\.6Hyperparameter sensitivity

Tables[6](https://arxiv.org/html/2608.05742#S4.T6)and[6](https://arxiv.org/html/2608.05742#S4.T6)report the sensitivity of CvLoss to the loss weightα\\alphaon TQNet and TimeFilter\. Increasingα\\alphafrom 0 yields consistent improvements, indicating the benefit of incorporating cross\-variable structural dependencies into the learning objective across different backbones and datasets\. For example, CvLoss reduces the MSE of TQNet from 0\.379 to 0\.365 on ETTh2, and that of TimeFilter from 0\.159 to 0\.157 on ECL\. However, using the structural penalty \(α=1\\alpha=1\) is not always optimal; on ECL, both models perform best aroundα=0\.5\\alpha=0\.5and degrade whenα\\alphareaches 1\. This suggests that CvLoss is most effective when combined with point\-wise MSE, which anchors forecasting accuracy and prevents translation\-invariant solutions under structural regularization\.

Table 5:Varyingα\\alpharesults of TQNet- •Note:Boldandunderlineddenote the best and second\-best results\.

Table 6:Varyingα\\alpharesults of TimeFilter\.- •Note:Boldandunderlineddenote the best and second\-best results\.

## 5Conclusion

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\. Extensive experiments show that CvLoss consistently improves SOTA forecasting backbones, outperforms competitive learning objectives, and stays efficient in practice\. These results suggest that modeling the structural geometry of future variables is a simple but effective direction for improving multivariate forecasting\.

## 6Limitations

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\. Learning adaptive graph structures and patch sizes could further improve efficiency and interpretability\. In addition, this work mainly evaluates deterministic forecasting on regular benchmark datasets; extending CvLoss to probabilistic forecasting, irregular sampling, and missing\-value scenarios is left for future work\.

Two limitations concern the relationship between our analysis and the objective we actually optimise, and we state them explicitly\. First, the implemented loss is not the object that Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)literally implicates\. The theorem shows that a non\-diagonal precision matrix is necessary, and Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)identifies the corresponding Gaussian family, but that family is generated by a*squared*graph penalty; theℓ1\\ell\_\{1\}penalty we deploy defines a pairwise Gibbs field over the same graph that is not Gaussian\. The choice is empirical rather than a theoretical consequence: Appendix[C\.4](https://arxiv.org/html/2608.05742#A3.SS4)shows that both norms help and thatℓ1\\ell\_\{1\}helps more under our protocol\. The theory should therefore be read as motivating a relation\-aware constraint rather than as deriving this one\. Second, the default uniform fully connected graph assigns the same weight to every cross\-variable pair, whereas real dependencies are heterogeneous\. We use it deliberately: by Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)the edge set is the support of the off\-diagonal entries, so a complete graph is the least informative choice available rather than an estimated structure\. It is nonetheless an uninformative prior\. The heterogeneous alternatives already present in the paper, namely top\-KKselection by the measured structural discrepancy \(Appendix[D](https://arxiv.org/html/2608.05742#A4)\) and the synchronous and asynchronous supports of Section[4\.5](https://arxiv.org/html/2608.05742#S4.SS5), are only partial substitutes for a learned weighting\.

## References

- Liu et al\. \(2024\)Yong Liu, Tengge Hu, Haoran Zhang, Haixu Wu, Shiyu Wang, Lintao Ma, and Mingsheng Long\.itransformer: Inverted transformers are effective for time series forecasting\.In*International Conference on Learning Representations*, 2024\.
- Zhao and Shen \(2024\)Lifan Zhao and Yanyan Shen\.Rethinking channel dependence for multivariate time series forecasting: Learning from leading indicators\.In*The Twelfth International Conference on Learning Representations*, 2024\.
- Qiu et al\. \(2025a\)Xiangfei Qiu, Xingjian Wu, Yan Lin, Chenjuan Guo, Jilin Hu, and Bin Yang\.Duet: Dual clustering enhanced multivariate time series forecasting\.In*SIGKDD*, pages 1185–1196, 2025a\.
- Hu et al\. \(2025a\)Yifan Hu, Guibin Zhang, Peiyuan Liu, Disen Lan, Naiqi Li, Dawei Cheng, Tao Dai, Shu\-Tao Xia, and Shirui Pan\.Timefilter: Patch\-specific spatial\-temporal graph filtration for time series forecasting\.In*Forty\-second International Conference on Machine Learning*, 2025a\.
- Ding et al\. \(2025\)Kuiye Ding, Fanda Fan, Yao Wang, Ruijie Jian, Xiaorui Wang, Luqi Gong, Yishan Jiang, Chunjie Luo, and Jianfeng Zhan\.Dualsg: A dual\-stream explicit semantic\-guided multivariate time series forecasting framework\.In*Proceedings of the 33rd ACM International Conference on Multimedia*, pages 508–517, 2025\.
- Wang et al\. \(2025a\)Hao Wang, Licheng Pan, Zhichao Chen, Degui Yang, Sen Zhang, Yifei Yang, Xinggao Liu, Haoxuan Li, and Dacheng Tao\.Fredf: Learning to forecast in the frequency domain\.In*ICLR*, 2025a\.
- Hu et al\. \(2025b\)Yifan Hu, Yuante Li, Peiyuan Liu, Yuxia Zhu, Naiqi Li, Tao Dai, Shu tao Xia, Dawei Cheng, and Changjun Jiang\.Fintsb: A comprehensive and practical benchmark for financial time series forecasting\.*arXiv preprint arXiv:2502\.18834*, 2025b\.
- Hu et al\. \(2025c\)Yifan Hu, Peiyuan Liu, Yuante Li, Dawei Cheng, Naiqi Li, Tao Dai, Jigang Bao, and Xia Shu\-Tao\.Finmamba: Market\-aware graph enhanced multi\-level mamba for stock movement prediction\.*arXiv preprint arXiv:2502\.06707*, 2025c\.
- Wang et al\. \(2024\)Yuxuan Wang, Haixu Wu, Jiaxiang Dong, Yong Liu, Mingsheng Long, and Jianmin Wang\.Deep time series models: A comprehensive survey and benchmark\.2024\.
- Qiu et al\. \(2026\)Xiangfei Qiu, Hanyin Cheng, Xingjian Wu, Junkai Lu, Jilin Hu, Chenjuan Guo, Christian S\. Jensen, and Bin Yang\.A comprehensive survey of deep learning for multivariate time series forecasting: A channel strategy perspective, 2026\.
- Liu et al\. \(2021\)Xiang Liu, Naiqi Li, and Shu\-Tao Xia\.Gdtw: A novel differentiable dtw loss for time series tasks\.In*Proc\. IEEE Int\. Conf\. Acoust\. Speech Signal Process\.*, pages 2860–2864\. IEEE, 2021\.
- Le Guen and Thome \(2019\)Vincent Le Guen and Nicolas Thome\.Shape and time distortion loss for training deep time series forecasting models\.*Proc\. Adv\. Neural Inf\. Process\. Syst\.*, 32, 2019\.
- Cuturi and Blondel \(2017\)Marco Cuturi and Mathieu Blondel\.Soft\-dtw: a differentiable loss function for time\-series\.In*Proc\. Int\. Conf\. Mach\. Learn\.*, pages 894–903\. PMLR, 2017\.
- Wang et al\. \(2026a\)Hao Wang, Licheng Pan, Yuan Lu, Zhichao Chen, Tianqiao Liu, Shuting He, Zhixuan Chu, Qingsong Wen, Haoxuan Li, and Zhouchen Lin\.Quadratic direct forecast for training multi\-step time\-series forecast models\.In*ICLR*, pages 1–9, 2026a\.
- Wang et al\. \(2025b\)Hao Wang, Licheng Pan, Zhichao Chen, Xu Chen, Qingyang Dai, Lei Wang, Haoxuan Li, and Zhouchen Lin\.Time\-o1: Time\-series forecasting needs transformed label alignment\.*Proc\. Adv\. Neural Inf\. Process\. Syst\.*, 2025b\.
- Lange et al\. \(2021\)Henning Lange, Steven L Brunton, and J Nathan Kutz\.From fourier to koopman: Spectral methods for long\-term time series prediction\.*Journal of Machine Learning Research*, 22\(41\):1–38, 2021\.
- Pan et al\. \(2026\)Licheng Pan, Hao Wang, Haocheng Yang, Yuqi Li, Qingsong Wen, Xiaoxi Li, Zhichao Chen, Haoxuan Li, Zhixuan Chu, and Yuan Lu\.Deep time\-series forecasting needs kernelized moment balancing, 2026\.
- Wang et al\. \(2026b\)Hao Wang, Licheng Pan, Yuan Lu, Zhixuan Chu, Xiaoxi Li, Shuting He, Zhichao Chen, Haoxuan Li, Qingsong Wen, and Zhouchen Lin\.Distdf: Time\-series forecasting needs joint\-distribution wasserstein alignment\.In*Proc\. Int\. Conf\. Learn\. Represent\.*, pages 1–9, 2026b\.
- Qiu et al\. \(2025b\)Xiangfei Qiu, Xingjian Wu, Hanyin Cheng, Xvyuan Liu, Chenjuan Guo, Jilin Hu, and Bin Yang\.Dbloss: Decomposition\-based loss function for time series forecasting\.In*NeurIPS*, 2025b\.
- Xiong et al\. \(2025\)Qi Xiong, Kai Tang, Minbo Ma, Ji Zhang, Jie Xu, and Tianrui Li\.Modeling temporal dependencies within the target for long\-term time series forecasting\.*IEEE Transactions on Knowledge and Data Engineering*, 37\(12\):7300–7314, 2025\.doi:10\.1109/TKDE\.2025\.3609415\.
- Kudrat et al\. \(2025\)Dilfira Kudrat, Zongxia Xie, Yanru Sun, Tianyu Jia, and Qinghua Hu\.Patch\-wise structural loss for time series forecasting\.In*Forty\-second International Conference on Machine Learning*, 2025\.
- Liang et al\. \(2024\)Yuxuan Liang, Haomin Wen, Yuqi Nie, Yushan Jiang, Ming Jin, Dongjin Song, Shirui Pan, and Qingsong Wen\.Foundation models for time series analysis: A tutorial and survey\.In*Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining*, 2024\.
- Chen et al\. \(2023\)Si\-An Chen, Chun\-Liang Li, Sercan O Arik, Nathanael Christian Yoder, and Tomas Pfister\.TSMixer: An all\-MLP architecture for time series forecast\-ing\.*Transactions on Machine Learning Research*, 2023\.ISSN 2835\-8856\.
- Zhang and Yan \(2023\)Yunhao Zhang and Junchi Yan\.Crossformer: Transformer utilizing cross\-dimension dependency for multivariate time series forecasting\.In*International Conference on Learning Representations*, 2023\.
- Wu et al\. \(2023\)Haixu Wu, Tengge Hu, Yong Liu, Hang Zhou, Jianmin Wang, and Mingsheng Long\.Timesnet: Temporal 2d\-variation modeling for general time series analysis\.In*International Conference on Learning Representations*, 2023\.
- Wu et al\. \(2020\)Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang\.Connecting the dots: Multivariate time series forecasting with graph neural networks\.In*Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining*, 2020\.
- Han et al\. \(2024a\)Wenyong Han, Tao Zhu, Liming Chen, Huansheng Ning, Yang Luo, and Yaping Wan\.Mcformer: Multivariate time series forecasting with mixed\-channels transformer\.*IEEE Internet of Things Journal*, 11\(17\):28320–28329, 2024a\.
- Zhao et al\. \(2024\)TingYu Zhao, Bo Peng, Yuan Sun, DaiPeng Yang, ZhenGuang Zhang, and Xi Wu\.Rethinking superpixel segmentation from biologically inspired mechanisms\.*Applied Soft Computing*, 156:111467, 2024\.
- Wei et al\. \(2026\)Jiawen Wei, Lan Jiang, Pengbo Wei, Ziwen Ye, Teng Song, Chen Chen, and Guangrui Ma\.TARFVAE: Efficient one\-step generative time series forecasting via TARFLOW based VAE\.In*The Thirty\-ninth Annual Conference on Neural Information Processing Systems*, 2026\.
- Evgeniou et al\. \(2005\)Theodoros Evgeniou, Charles A\. Micchelli, and Massimiliano Pontil\.Learning multiple tasks with kernel methods\.*Journal of Machine Learning Research*, 6\(21\):615–637, 2005\.
- Nie et al\. \(2023\)Yuqi Nie, Nam H Nguyen, Phanwadee Sinthong, and Jayant Kalagnanam\.A time series is worth 64 words: Long\-term forecasting with transformers\.In*International Conference on Learning Representations*, 2023\.
- Ding et al\. \(2026\)Kuiye Ding, Fanda Fan, Chunyi Hou, Zheya Wang, Lei Wang, Zhengxin Yang, and Jianfeng Zhan\.Timemosaic: Temporal heterogeneity guided time series forecasting via adaptive granularity patch and segment\-wise decoding\.*Proceedings of the AAAI Conference on Artificial Intelligence*, 40\(25\):20790–20798, Mar\. 2026\.doi:10\.1609/aaai\.v40i25\.39218\.URL[https://ojs\.aaai\.org/index\.php/AAAI/article/view/39218](https://ojs.aaai.org/index.php/AAAI/article/view/39218)\.
- Lin et al\. \(2025\)Shengsheng Lin, Haojun Chen, Haijie Wu, Chunyun Qiu, and Weiwei Lin\.Temporal query network for efficient multivariate time series forecasting\.In*Forty\-second International Conference on Machine Learning*, 2025\.
- Zhou et al\. \(2022\)Tian Zhou, Ziqing Ma, Qingsong Wen, Xue Wang, Liang Sun, and Rong Jin\.FEDformer: Frequency enhanced decomposed transformer for long\-term series forecasting\.In*International Conference on Machine Learning*, pages 27268–27286\. PMLR, 2022\.
- Yu et al\. \(2024\)Guoqi Yu, Jing Zou, Xiaowei Hu, Angelica I Aviles\-Rivero, Jing Qin, and Shujun Wang\.Revitalizing multivariate time series forecasting: Learnable decomposition with inter\-series dependencies and intra\-series variations modeling\.In*Forty\-first International Conference on Machine Learning*, 2024\.
- Zeng et al\. \(2023\)Ailing Zeng, Muxi Chen, Lei Zhang, and Qiang Xu\.Are transformers effective for time series forecasting?In*Proc\. AAAI Conf\. Artif\. Intell\.*, 2023\.
- Han et al\. \(2024b\)Lu Han, Xu\-Yang Chen, Han\-Jia Ye, and De\-Chuan Zhan\.SOFTS: Efficient multivariate time series forecasting with series\-core fusion\.In*Advances in Neural Information Processing Systems*, 2024b\.
- Qiu et al\. \(2024\)Xiangfei Qiu, Jilin Hu, Lekui Zhou, Xingjian Wu, Junyang Du, Buang Zhang, Chenjuan Guo, Aoying Zhou, Christian S\. Jensen, Zhenli Sheng, and Bin Yang\.Tfb: Towards comprehensive and fair benchmarking of time series forecasting methods\.In*VLDB*, pages 2363–2377, 2024\.
- Kingma and Ba \(2015\)Diederik P\. Kingma and Jimmy Ba\.Adam: A method for stochastic optimization\.In*International Conference on Learning Representations*, pages 1–9, 2015\.
- Dai et al\. \(2024\)Tao Dai, Beiliang Wu, Peiyuan Liu, Naiqi Li, Jigang Bao, Yong Jiang, and Shu\-Tao Xia\.Periodicity decoupling framework for long\-term series forecasting\.*International Conference on Learning Representations*, 2024\.
- Kou et al\. \(2025\)Feifei Kou, Jiahao Wang, Lei Shi, Yuhan Yao, Yawen Li, Suguo Zhu, Zhongbao Zhang, and Junping Du\.Cfpt: Empowering time series forecasting through cross\-frequency interaction and periodic\-aware timestamp modeling\.In*Forty\-second International Conference on Machine Learning*, 2025\.
- Li et al\. \(2021\)Jianxin Li, Xiong Hui, and Wancai Zhang\.Informer: Beyond efficient transformer for long sequence time\-series forecasting\.In*The Thirty\-Fifth AAAI Conference on Artificial Intelligence*, 2021\.
- Wu et al\. \(2021\)Haixu Wu, Jiehui Xu, Jianmin Wang, and Mingsheng Long\.Autoformer: Decomposition transformers with Auto\-Correlation for long\-term series forecasting\.In*Advances in Neural Information Processing Systems*, 2021\.
- Liu et al\. \(2022\)Minhao Liu, Ailing Zeng, Muxi Chen, Zhijian Xu, Qiuxia Lai, Lingna Ma, and Qiang Xu\.Scinet: time series modeling and forecasting with sample convolution and interaction\.In*Advances in Neural Information Processing Systems*, 2022\.
- Liu et al\. \(2025\)Peiyuan Liu, Beiliang Wu, Yifan Hu, Naiqi Li, Tao Dai, Jigang Bao, and Shu\-Tao Xia\.Timebridge: Non\-stationarity matters for long\-term time series forecasting\.*International Conference on Machine Learning*, 2025\.

## Appendix ATheoretical Justification

In this section, we provide the detailed derivation of Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1), which reveals the suboptimality of the standard DF objective under real\-world spatiotemporal correlations, followed by the proof of Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)\.

### A\.1Proof of Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)

Proof\.From a probabilistic perspective, finding the optimal parametersθ\\thetafor the forecasting modelfθf\_\{\\theta\}is equivalent to maximizing the likelihood of the ground\-truth targetYYgiven the input historyXX\. Let

𝐞=vec​\(Y^−Y\)∈ℝT​D\\mathbf\{e\}=\\mathrm\{vec\}\(\\hat\{Y\}\-Y\)\\in\\mathbb\{R\}^\{TD\}be the flattened forecast error vector\.

Assume that the ground\-truth residuals follow a zero\-mean multivariate Gaussian distribution, i\.e\.,

𝐞∼𝒩​\(𝟎,𝚺S​T\),\\mathbf\{e\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\mathbf\{\\Sigma\}\_\{ST\}\),where𝚺S​T∈ℝT​D×T​D\\mathbf\{\\Sigma\}\_\{ST\}\\in\\mathbb\{R\}^\{TD\\times TD\}is the true spatiotemporal covariance matrix\. The probability density function of the residuals is given by:

p​\(𝐞\)=1\(2​π\)T​D​det\(𝚺S​T\)​exp⁡\(−12​𝐞⊤​𝚺S​T−1​𝐞\)p\(\\mathbf\{e\}\)=\\frac\{1\}\{\\sqrt\{\(2\\pi\)^\{TD\}\\det\(\\mathbf\{\\Sigma\}\_\{ST\}\)\}\}\\exp\\left\(\-\\frac\{1\}\{2\}\\mathbf\{e\}^\{\\top\}\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}\\mathbf\{e\}\\right\)
Minimizing the NLL of this distribution yields the optimal objective function\. Taking the negative logarithm, we have:

−log⁡p​\(𝐞\)=T​D2​log⁡\(2​π\)\+12​log​det\(𝚺S​T\)\+12​𝐞⊤​𝚺S​T−1​𝐞\-\\log p\(\\mathbf\{e\}\)=\\frac\{TD\}\{2\}\\log\(2\\pi\)\+\\frac\{1\}\{2\}\\log\\det\(\\mathbf\{\\Sigma\}\_\{ST\}\)\+\\frac\{1\}\{2\}\\mathbf\{e\}^\{\\top\}\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}\\mathbf\{e\}
Assuming the covariance matrix𝚺S​T\\mathbf\{\\Sigma\}\_\{ST\}represents the intrinsic data distribution and is independent of the model parametersθ\\theta, the first two terms are constants with respect to the optimization process\. By dropping these constants and scaling the objective by1T​D\\frac\{1\}\{TD\}to obtain the average loss per prediction entry, we define the practical NLL objective as:

ℒnll=12​T​D​𝐞⊤​𝚺S​T−1​𝐞\\mathcal\{L\}\_\{\\mathrm\{nll\}\}=\\frac\{1\}\{2TD\}\\mathbf\{e\}^\{\\top\}\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}\\mathbf\{e\}
Now, consider the standard point\-wise MSE objective used in the DF paradigm, defined as:

ℒdf=1T​∑t=1T1D​‖Yt,:−Y^t,:‖22=1T​D​𝐞⊤​𝐞\\mathcal\{L\}\_\{\\mathrm\{df\}\}=\\frac\{1\}\{T\}\\sum\_\{t=1\}^\{T\}\\frac\{1\}\{D\}\\\|Y\_\{t,:\}\-\\hat\{Y\}\_\{t,:\}\\\|\_\{2\}^\{2\}=\\frac\{1\}\{TD\}\\mathbf\{e\}^\{\\top\}\\mathbf\{e\}
Notice thatℒdf\\mathcal\{L\}\_\{\\mathrm\{df\}\}can be probabilistically interpreted as the negative log\-likelihood under a strictly naive assumption: the residuals of all variables across all time steps are entirely independent and share a uniform varianceσ2\\sigma^\{2\}\. Under this assumption, the covariance matrix reduces to a spherical structure𝚺S​T=σ2​𝐈\\mathbf\{\\Sigma\}\_\{ST\}=\\sigma^\{2\}\\mathbf\{I\}, where𝐈∈ℝT​D×T​D\\mathbf\{I\}\\in\\mathbb\{R\}^\{TD\\times TD\}is the identity matrix\. In this special case, the NLL naturally degrades to:

ℒnll​\_​naive=12​T​D​𝐞⊤​\(σ2​𝐈\)−1​𝐞=12​σ2​\(1T​D​𝐞⊤​𝐞\)=12​σ2​ℒdf\\mathcal\{L\}\_\{\\mathrm\{nll\\\_naive\}\}=\\frac\{1\}\{2TD\}\\mathbf\{e\}^\{\\top\}\(\\sigma^\{2\}\\mathbf\{I\}\)^\{\-1\}\\mathbf\{e\}=\\frac\{1\}\{2\\sigma^\{2\}\}\\left\(\\frac\{1\}\{TD\}\\mathbf\{e\}^\{\\top\}\\mathbf\{e\}\\right\)=\\frac\{1\}\{2\\sigma^\{2\}\}\\mathcal\{L\}\_\{\\mathrm\{df\}\}
To formally quantify the structural mismatch between the standard DF paradigm and the optimal estimation under true spatiotemporal correlations, we define the structural gapΔ\\Deltaas the difference between the true optimal objectiveℒnll\\mathcal\{L\}\_\{\\mathrm\{nll\}\}and the scaled DF objective12​σ2​ℒdf\\frac\{1\}\{2\\sigma^\{2\}\}\\mathcal\{L\}\_\{\\mathrm\{df\}\}:

Δ=ℒnll−12​σ2​ℒdf=12​T​D​𝐞⊤​𝚺S​T−1​𝐞−12​T​D​σ2​𝐞⊤​𝐈𝐞\\Delta=\\mathcal\{L\}\_\{\\mathrm\{nll\}\}\-\\frac\{1\}\{2\\sigma^\{2\}\}\\mathcal\{L\}\_\{\\mathrm\{df\}\}=\\frac\{1\}\{2TD\}\\mathbf\{e\}^\{\\top\}\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}\\mathbf\{e\}\-\\frac\{1\}\{2TD\\sigma^\{2\}\}\\mathbf\{e\}^\{\\top\}\\mathbf\{I\}\\mathbf\{e\}
Factoring out the common terms, we obtain:

Δ=12​T​D​𝐞⊤​\(𝚺S​T−1−1σ2​𝐈\)​𝐞\\Delta=\\frac\{1\}\{2TD\}\\mathbf\{e\}^\{\\top\}\\left\(\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}\-\\frac\{1\}\{\\sigma^\{2\}\}\\mathbf\{I\}\\right\)\\mathbf\{e\}
In real\-world multivariate time series, dynamic physical processes \(e\.g\., disturbance propagation and lagged correlations\) dictate that𝚺S​T\\mathbf\{\\Sigma\}\_\{ST\}is neither diagonal nor uniform\. Consequently,𝚺S​T−1≠1σ2​𝐈\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}\\neq\\frac\{1\}\{\\sigma^\{2\}\}\\mathbf\{I\}, meaning the structural gapΔ≠0\\Delta\\neq 0\. Therefore, relying exclusively onℒdf\\mathcal\{L\}\_\{\\mathrm\{df\}\}implicitly ignores the off\-diagonal elements of the precision matrix𝚺S​T−1\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}, leading to statistically suboptimal parameter estimation\. This concludes the proof\.

### A\.2Proof of Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)

Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)shows that MSE is the correct objective only when the residual precision matrix is a multiple of the identity, but it does not say which non\-diagonal precision matrix should replace it, and estimating𝚺S​T−1\\mathbf\{\\Sigma\}\_\{ST\}^\{\-1\}from data is not feasible at the scales considered here\. Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)closes this gap from the opposite direction: instead of estimating a precision matrix, we fix a one\-parameter family of precision matrices whose off\-diagonal support is a cross\-variable graph, and show that the negative log\-likelihood of that family is exactly MSE plus a squared graph total\-variation penalty\.

Proof\.Recall from Section[3\.2](https://arxiv.org/html/2608.05742#S3.SS2)thatA∈ℝ\|ℰ\|×NA\\in\\mathbb\{R\}^\{\|\\mathcal\{E\}\|\\times N\}is the incidence matrix of the cross\-variable graph under an arbitrary edge orientation, thatA~=A⊗IL\\tilde\{A\}=A\\otimes I\_\{L\}is its lifted counterpart, and that𝐫=vec​\(E\)∈ℝN​L\\mathbf\{r\}=\\mathrm\{vec\}\(E\)\\in\\mathbb\{R\}^\{NL\}collects the node\-wise residuals𝐞1,…,𝐞N∈ℝL\\mathbf\{e\}\_\{1\},\\dots,\\mathbf\{e\}\_\{N\}\\in\\mathbb\{R\}^\{L\}\. By construction, thekk\-th block ofA~​𝐫\\tilde\{A\}\\mathbf\{r\}associated with the edgek=\(i,j\)k=\(i,j\)equals𝐞i−𝐞j\\mathbf\{e\}\_\{i\}\-\\mathbf\{e\}\_\{j\}, so

‖A~​𝐫‖22=∑\(i,j\)∈ℰ‖𝐞i−𝐞j‖22\.\\\|\\tilde\{A\}\\mathbf\{r\}\\\|\_\{2\}^\{2\}=\\sum\_\{\(i,j\)\\in\\mathcal\{E\}\}\\\|\\mathbf\{e\}\_\{i\}\-\\mathbf\{e\}\_\{j\}\\\|\_\{2\}^\{2\}\.WritingLap​\(ℰ\)=A~⊤​A~\\mathrm\{Lap\}\(\\mathcal\{E\}\)=\\tilde\{A\}^\{\\top\}\\tilde\{A\}for the lifted graph Laplacian, the same quantity is the associated quadratic form,

𝐫⊤​Lap​\(ℰ\)​𝐫=𝐫⊤​A~⊤​A~​𝐫=‖A~​𝐫‖22=∑\(i,j\)∈ℰ‖𝐞i−𝐞j‖22\.\\mathbf\{r\}^\{\\top\}\\mathrm\{Lap\}\(\\mathcal\{E\}\)\\,\\mathbf\{r\}=\\mathbf\{r\}^\{\\top\}\\tilde\{A\}^\{\\top\}\\tilde\{A\}\\,\\mathbf\{r\}=\\\|\\tilde\{A\}\\mathbf\{r\}\\\|\_\{2\}^\{2\}=\\sum\_\{\(i,j\)\\in\\mathcal\{E\}\}\\\|\\mathbf\{e\}\_\{i\}\-\\mathbf\{e\}\_\{j\}\\\|\_\{2\}^\{2\}\.Substituting the structured precision family of Eq\. \([3](https://arxiv.org/html/2608.05742#S3.E3)\) and expanding by linearity gives

𝐫⊤​𝚺−1​\(λ\)​𝐫=1σ2​\(𝐫⊤​𝐫\+λ​𝐫⊤​Lap​\(ℰ\)​𝐫\)=1σ2​\(‖𝐫‖22\+λ​∑\(i,j\)∈ℰ‖𝐞i−𝐞j‖22\),\\mathbf\{r\}^\{\\top\}\\mathbf\{\\Sigma\}^\{\-1\}\(\\lambda\)\\,\\mathbf\{r\}=\\frac\{1\}\{\\sigma^\{2\}\}\\Bigl\(\\mathbf\{r\}^\{\\top\}\\mathbf\{r\}\+\\lambda\\,\\mathbf\{r\}^\{\\top\}\\mathrm\{Lap\}\(\\mathcal\{E\}\)\\,\\mathbf\{r\}\\Bigr\)=\\frac\{1\}\{\\sigma^\{2\}\}\\Bigl\(\\\|\\mathbf\{r\}\\\|\_\{2\}^\{2\}\+\\lambda\\\!\\\!\\sum\_\{\(i,j\)\\in\\mathcal\{E\}\}\\\!\\\!\\\|\\mathbf\{e\}\_\{i\}\-\\mathbf\{e\}\_\{j\}\\\|\_\{2\}^\{2\}\\Bigr\),which is Eq\. \([4](https://arxiv.org/html/2608.05742#S3.E4)\)\. This concludes the proof\.

##### Remarks\.

Four consequences are worth recording, since they determine how the theory should and should not be read\.

- •The family is a valid precision family\.Lap​\(ℰ\)\\mathrm\{Lap\}\(\\mathcal\{E\}\)is positive semi\-definite, so𝐈\+λ​Lap​\(ℰ\)≻0\\mathbf\{I\}\+\\lambda\\,\\mathrm\{Lap\}\(\\mathcal\{E\}\)\\succ 0for everyλ≥0\\lambda\\geq 0and𝚺−1​\(λ\)\\mathbf\{\\Sigma\}^\{\-1\}\(\\lambda\)is a genuine precision matrix\. Its off\-diagonal entries are−λ/σ2\-\\lambda/\\sigma^\{2\}exactly on the node pairs inℰ\\mathcal\{E\}and zero elsewhere, soℰ\\mathcal\{E\}is the conditional\-independence structure of the corresponding Gaussian field rather than an analogy\.
- •It contains the spherical case\.Atλ=0\\lambda=0we recover𝚺−1=σ−2​𝐈\\mathbf\{\\Sigma\}^\{\-1\}=\\sigma^\{\-2\}\\mathbf\{I\}, which is precisely the situation in which the objective gapΔ\\Deltaof Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)vanishes\. The theorem and the loss are therefore two ends of one family rather than unrelated statements\.
- •It matches the deployed weighting\.Dividing the objectiveℒα=\(1−α\)​ℒdf\+α​ℒcv\\mathcal\{L\}\_\{\\alpha\}=\(1\-\\alpha\)\\mathcal\{L\}\_\{\\mathrm\{df\}\}\+\\alpha\\mathcal\{L\}\_\{\\mathrm\{cv\}\}by\(1−α\)\(1\-\\alpha\)leaves the minimiser unchanged and yieldsℒdf\+λ​ℒcv\\mathcal\{L\}\_\{\\mathrm\{df\}\}\+\\lambda\\mathcal\{L\}\_\{\\mathrm\{cv\}\}withλ=α/\(1−α\)\\lambda=\\alpha/\(1\-\\alpha\)\. Values ofα\\alphaclose to11therefore correspond to a largeλ\\lambdawithin the family above, not to a weight outside the stated range\(0,1\)\(0,1\)\.
- •The complete graph is the weakest member\.Sinceℰ\\mathcal\{E\}is the support of the off\-diagonal entries, takingℰ\\mathcal\{E\}to be the complete cross\-variable graph imposes*no*zero constraint on any cross\-variable pair; sparser choices impose strictly more\. The empty graph, i\.e\. plain MSE, is the member excluded by Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)whenever cross\-variable dependencies are present\.

Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)concerns the squared penalty\. The objective we deploy in Eq\. \([12](https://arxiv.org/html/2608.05742#S3.E12)\) uses theℓ1\\ell\_\{1\}norm on the same edge differences, which yields the pairwise field

p1​\(𝐫∣ℰ\)∝exp⁡\[−a​∑v=1N‖𝐞v‖22−b​∑\(i,j\)∈ℰ‖𝐞i−𝐞j‖1\],a=1−αN​L\>0,b=α\|ℰ\|​L\>0\.p\_\{1\}\(\\mathbf\{r\}\\mid\\mathcal\{E\}\)\\;\\propto\\;\\exp\\Bigl\[\-a\\sum\_\{v=1\}^\{N\}\\\|\\mathbf\{e\}\_\{v\}\\\|\_\{2\}^\{2\}\-b\\\!\\\!\\sum\_\{\(i,j\)\\in\\mathcal\{E\}\}\\\!\\\!\\\|\\mathbf\{e\}\_\{i\}\-\\mathbf\{e\}\_\{j\}\\\|\_\{1\}\\Bigr\],\\qquad a=\\frac\{1\-\\alpha\}\{NL\}\>0,\\quad b=\\frac\{\\alpha\}\{\|\\mathcal\{E\}\|L\}\>0\.This factorises into node and edge potentials and is therefore a pairwise Gibbs Markov random field with respect to the same graph, but it is not Gaussian, and we consequently restrict the precision\-matrix reading to theℓ2\\ell\_\{2\}case\. We also do not claim thatℰ\\mathcal\{E\}recovers a sparse precision support: the default topology is complete, with random edge sampling used only for scalability\. Theℓ1\\ell\_\{1\}form is adopted because squared penalties can be dominated by a small number of large edge\-wise discrepancies arising from localized anomalies or variable\-specific shocks\. This is a task\-specific empirical hypothesis about time\-series residuals rather than a universal claim\. A comparable separation between the norm used for analysis and the norm used in practice appears in FreDF\[Wang et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib6)\], whose analysis isℓ2\\ell\_\{2\}/Gaussian while its deployed frequency\-domain loss isℓ1\\ell\_\{1\}\. Appendix[C\.4](https://arxiv.org/html/2608.05742#A3.SS4)tests the two norms directly\.

## Appendix BReproduction Details

### B\.1Dataset descriptions

Our empirical evaluation is conducted on a diverse collection of widely\-used time series forecasting benchmarks\. Each dataset presents distinct characteristics in terms of dimensionality and temporal resolution\. A summary is provided in[Table 7](https://arxiv.org/html/2608.05742#A2.T7)\.

- •ETT\[Li et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib42)\]: This dataset includes seven electricity transformer indicators logged between July 2016 and July 2018\. It is partitioned into four subsets according to sampling rates: ETTh1 and ETTh2 for hourly data, and ETTm1 and ETTm2 for 15\-minute intervals\.
- •Weather\[Wu et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib43)\]: Encompasses 21 meteorological factors collected every 10 minutes throughout the year 2020 from the Max Planck Institute for Biogeochemistry\.
- •ECL\[Wu et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib43)\]: Records the hourly power usage patterns of 321 individual customers\.
- •Traffic\[Wu et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib43)\]: Monitors hourly road occupancy via 862 sensors located on freeways in the San Francisco Bay Area from 2015 to 2016\.
- •
- •PEMS\[Liu et al\.,[2022](https://arxiv.org/html/2608.05742#bib.bib44)\]: Derived from public California highway traffic records, aggregated at 5\-minute steps\. We adopt two standard subsets, PEMS03 and PEMS08, for our study\.

Table 7:Dataset description\.- •Note:Ddenotes the number of variates\.*Frequency*denotes the sampling interval of time points\.*Train, Validation, Test*denotes the number of samples employed in each split\. The taxonomy aligns with\[Wu et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib25)\]\.

##### Why these datasets\.

The twelve benchmarks above are the standard suite of\[Wu et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib25), Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\]and were selected before any results were obtained; we record the reasoning here because the choice interacts with what CvLoss is designed to do\. Three considerations apply\. First, we require coverage of the variable\-count range, since CvLoss regularizes interactions between variables and its cost and its benefit both depend onDD: the suite spansD=7D=7\(the four ETT subsets\) toD=862D=862\(Traffic\), with Weather \(2121\), Solar \(137137\), PEMS08 \(170170\), PEMS04 \(307307\), ECL \(321321\), PEMS03 \(358358\) and PEMS07 \(883883\) in between\. Second, we require coverage of the cross\-variable correlation range, which Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)predicts should govern the size of the effect; the suite spans weakly correlated channels \(ETT, average CCC0\.140\.14–0\.180\.18\) to strongly correlated ones \(Solar,0\.910\.91; PEMS07,0\.780\.78\), as reported in Table[20](https://arxiv.org/html/2608.05742#A4.T20)\. Third, both the long\-horizon regime \(T∈\{96,192,336,720\}T\\in\\\{96,192,336,720\\\}\) and the short\-horizon regime used for PEMS \(T∈\{12,24,36,48\}T\\in\\\{12,24,36,48\\\}\) are represented, so the conclusions are not specific to one horizon scale\. All four ETT subsets are included, at both the hourly \(ETTh1, ETTh2\) and the1515\-minute \(ETTm1, ETTm2\) sampling rates\.

We deliberately exclude the Exchange dataset, following the protocol of\[Hu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib4)\]\. Exchange is close to a driftless random walk, a property documented since\[Liu et al\.,[2022](https://arxiv.org/html/2608.05742#bib.bib44)\], and on such a series the trivial predictor that repeats the last observed value is competitive with deep forecasters\. A benchmark on which the degenerate baseline is already near\-optimal cannot discriminate between learning objectives, which is the comparison this paper makes\. The point is compounded for our specific claim: Exchange consists of eight loosely coupled national exchange rates, so it carries little of the cross\-variable structure that CvLoss is designed to exploit, and including it would neither support nor challenge the mechanism under test\. We note the exclusion explicitly rather than silently, since the omission is a choice about what the evidence can show\.

Following established protocols\[Qiu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib38), Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\], all datasets are segmented chronologically into training, validation, and testing portions\. Regarding the ETT, Weather, ECL, Solar and Traffic datasets, we maintain a constant input window of 96 time steps and assess model performance across four distinct forecasting horizons: 96, 192, 336, and 720\. For the PEMS datasets, while the historical look\-back period remains 96, the evaluation focuses on shorter prediction intervals of 12, 24, 36, and 48 steps\. Throughout our experimental phase, we disable thedropping\-last trickduring the final test set assessment to ensure that every data point in the terminal batch is fully utilized\.

### B\.2Implementation details of model training

To ensure a strictly fair comparison, we build all experiments on top of publicly released benchmark implementations and retain the original model architectures and hyperparameter configurations\. When integrating CvLoss, we do not modify the forecasting backbone, data preprocessing, optimizer, training schedule, or evaluation protocol; only the proposed objective term and its associated weight are added\. We also disable the drop\-last trick for all methods following\[Qiu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib38)\], so that every test sample is evaluated consistently\. To establish a fair comparison, we reproduced all baseline models using their official, publicly available implementations, primarily sourcing from the TQNet\[Lin et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib33)\], TimeFilter\[Hu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib4)\], iTransformer\[Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\], TimeBridge\[Liu et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib45)\]and Time\-Series\-Library\[Wang et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib9)\]repositories\. The reproducibility of these baseline results was verified prior to our experiments\. All models were trained to minimize the MSE loss function using the Adam optimizer\[Kingma and Ba,[2015](https://arxiv.org/html/2608.05742#bib.bib39)\]\. All baselines apply to the publicly available parameters in their library without any changes\. To prevent overfitting, we employed an early stopping mechanism that terminates training if the validation loss fails to improve for three consecutive epochs\.

##### What is tuned, and on which split\.

Two quantities are tuned for CvLoss: the loss weightα\\alphaand the patch lengthLL\. Nothing else is\. Both are selected on the*validation*split and never on the test split, and all methods, baselines included, use early stopping on the validation loss with a patience of three epochs\. Every number reported in the main text and in the appendix uses the fixed convex objectiveℒα=\(1−α\)​ℒdf\+α​ℒcv\\mathcal\{L\}\_\{\\alpha\}=\(1\-\\alpha\)\\mathcal\{L\}\_\{\\mathrm\{df\}\}\+\\alpha\\mathcal\{L\}\_\{\\mathrm\{cv\}\}with a single scalarα∈\(0,1\)\\alpha\\in\(0,1\)chosen this way;α\\alphais never learned during training, and there is no configuration in which it varies across batches, variables or edges\. The sensitivity tables report the full sweep overα\\alpha, whereas the main tables report the validation\-selected setting; the two are consistent by construction, since theα=0\\alpha=0rows of Tables[12](https://arxiv.org/html/2608.05742#A3.T12)and[12](https://arxiv.org/html/2608.05742#A3.T12)reproduce the unregularised baselines exactly\. The separate quantities in Table[20](https://arxiv.org/html/2608.05742#A4.T20)are a diagnostic and are*not*values ofα\\alpha; this is stated with that table\.

##### Baselines are tuned at least as favourably as our method\.

This deserves a precise statement, because a plug\-in objective can appear to win simply by receiving more tuning effort than the models it is added to\. The comparison here is arranged so that the opposite holds\. Baselines are the authors’ own public implementations, run with the hyperparameters published in their libraries and verified to reproduce beforehand; we change none of them\. When CvLoss is added, the backbone, data preprocessing, normalisation, optimizer, learning\-rate schedule, batch size, early\-stopping rule and evaluation protocol are all left exactly as the baseline defines them, and only the additional objective term and its weight are introduced\. Our method therefore has strictly less freedom than the models it is compared against: two scalars against a fully tuned configuration\. The drop\-last trick is disabled for every model, ours included, following\[Qiu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib38)\], so no method benefits from a truncated final test batch\.

##### No ground\-truth future statistics enter the objective\.

In every configuration used for the reported results, the edge setℰ\\mathcal\{E\}is fixed a priori, since it is the complete cross\-variable graph of Section[3\.2](https://arxiv.org/html/2608.05742#S3.SS2), and therefore depends on no property of the future window\. The quantity CvLoss consumes is the residualE=Z^−ZE=\\hat\{Z\}\-Z, which is the same supervision signal that the standard DF objective already consumes in Eq\. \([1](https://arxiv.org/html/2608.05742#S3.E1)\); no additional statistic ofYYis used\. The one variant that inspects the ground\-truth correlation matrix𝐂​\(Y\)\\mathbf\{C\}\(Y\)is the top\-KKedge selection of Appendix[D](https://arxiv.org/html/2608.05742#A4), which is used*during training only*to concentrate computation on the worst\-aligned variable pairs, and which produced none of the numbers in this paper\. In all configurations the CvLoss term leaves the computation graph at inference: it contributes no parameters, no buffers and no forward\-pass computation to the deployed model, which is why the inference latency in Figure[11](https://arxiv.org/html/2608.05742#A3.F11)sits at parity with the unmodified backbone\.

### B\.3Additional Visualization of Cross\-Variable Correlation Structures

![Refer to caption](https://arxiv.org/html/2608.05742v1/x14.png)\(a\)Ground\-truth structure\.![Refer to caption](https://arxiv.org/html/2608.05742v1/x15.png)\(b\)Baseline structure\.![Refer to caption](https://arxiv.org/html/2608.05742v1/x16.png)\(c\)Improvement by CvLoss\.![Refer to caption](https://arxiv.org/html/2608.05742v1/x17.png)\(d\)PEMS03 Snapshot\.
Figure 5:Motivating example on PEMS03 with input lengthH=96H=96and forecast horizonT=48T=48\. \(a,b\) Cross\-variable correlation structures of the ground truth and the plain iTransformer\[Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\]prediction\. \(c\) Entries where CvLoss reduces structural error over the baseline, covering 61\.98% of cross\-variable pairs\. \(d\) Forecast snapshot: orange denotes the baseline prediction, black denotes the ground truth, and blue denotes ours \(CvLoss\)\. See details in Appendix[B\.3](https://arxiv.org/html/2608.05742#A2.SS3)\.To further examine how CvLoss affects the recovery of cross\-variable dependencies, we provide additional visualizations on the ETTh2 and Weather datasets\. These cases follow the same qualitative protocol as Figure[1\(a\)](https://arxiv.org/html/2608.05742#S2.F1.sf1), while using input lengthH=96H=96and forecast horizonT=96T=96\.

LetY∈ℝT×DY\\in\\mathbb\{R\}^\{T\\times D\}denote the ground\-truth future series and letY^base,Y^cv∈ℝT×D\\hat\{Y\}^\{\\text\{base\}\},\\hat\{Y\}^\{\\text\{cv\}\}\\in\\mathbb\{R\}^\{T\\times D\}denote the future series predicted by the baseline iTransformer and by iTransformer trained with CvLoss, respectively\. For any future seriesS∈ℝT×DS\\in\\mathbb\{R\}^\{T\\times D\}, we define its cross\-variable correlation structure as

C​\(S\)=\[\|Corr​\(S:,i,S:,j\)\|\]i,j=1D∈ℝD×D,C\(S\)=\\left\[\\,\\big\|\\mathrm\{Corr\}\(S\_\{:,i\},S\_\{:,j\}\)\\big\|\\,\\right\]\_\{i,j=1\}^\{D\}\\in\\mathbb\{R\}^\{D\\times D\},whereCorr​\(⋅,⋅\)\\mathrm\{Corr\}\(\\cdot,\\cdot\)denotes the Pearson correlation coefficient\. Accordingly, we visualize

Cgt=C​\(Y\),Cbase=C​\(Y^base\)\.C^\{\\mathrm\{gt\}\}=C\(Y\),\\qquad C^\{\\mathrm\{base\}\}=C\(\\hat\{Y\}^\{\\mathrm\{base\}\}\)\.
To quantify structural deviation from the ground\-truth pattern, we define the structural error matrices as

SEbase=\|Cbase−Cgt\|,SEcv=\|C​\(Y^cv\)−Cgt\|,\\mathrm\{SE\}^\{\\mathrm\{base\}\}=\\big\|C^\{\\mathrm\{base\}\}\-C^\{\\mathrm\{gt\}\}\\big\|,\\qquad\\mathrm\{SE\}^\{\\mathrm\{cv\}\}=\\big\|C\(\\hat\{Y\}^\{\\mathrm\{cv\}\}\)\-C^\{\\mathrm\{gt\}\}\\big\|,where the absolute value is taken element\-wise\. The improvement induced by CvLoss is then measured by

Δ=SEbase−SEcv\.\\Delta=\\mathrm\{SE\}^\{\\mathrm\{base\}\}\-\\mathrm\{SE\}^\{\\mathrm\{cv\}\}\.A positive entryΔi​j\>0\\Delta\_\{ij\}\>0means that CvLoss produces a smaller structural error than the baseline on the\(i,j\)\(i,j\)\-th cross\-variable entry, i\.e\., the recovered correlation is closer to the ground\-truth structure\. Therefore, panel \(c\) visualizes the positive part ofΔ\\Delta, highlighting where CvLoss improves over the baseline in recovering cross\-variable dependencies\.

These supplementary cases complement the ECL example of Figure[1\(a\)](https://arxiv.org/html/2608.05742#S2.F1.sf1)in the main text and the PEMS03 example of Figure[5](https://arxiv.org/html/2608.05742#A2.F5), and verify that the structural benefit of CvLoss is not restricted to a single dataset\. As shown in Figures[6\(a\)](https://arxiv.org/html/2608.05742#A2.F6.sf1)and[7\(a\)](https://arxiv.org/html/2608.05742#A2.F7.sf1), CvLoss consistently improves the recovery of cross\-variable correlation structures across different data regimes, including both low\-dimensional and higher\-dimensional multivariate forecasting benchmarks\. This observation is consistent with our main claim that CvLoss acts as a structural regularizer on multivariate outputs and improves the preservation of cross\-variable dependencies beyond point\-wise supervision alone\.

![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/heatmap_vis/ETTh2_96_96_gt.png)\(a\)Ground\-truth structure\.![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/heatmap_vis/ETTh2_96_96_baseline.png)\(b\)Baseline structure\.![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/heatmap_vis/ETTh2_96_96_delta_positive.png)\(c\)Improvement by CvLoss\.Figure 6:The cross\-variable correlation structures of future series givenXXon the ETTh2 dataset, with input lengthH=96H=96and forecast horizonT=96T=96\. The correlation matrices are computed from the ground\-truth future series \(a\) and the future series predicted by iTransformer\[Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\]\(b\), while \(c\) shows the entries on which CvLoss yields smaller structural errors than iTransformer \(62%\)\.
![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/heatmap_vis/weather_96_96_gt.png)\(a\)Ground\-truth structure\.![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/heatmap_vis/weather_96_96_baseline.png)\(b\)Baseline structure\.![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/heatmap_vis/weather_96_96_delta_positive.png)\(c\)Improvement by CvLoss\.Figure 7:The cross\-variable correlation structures of future series givenXXon the Weather dataset, with input lengthH=96H=96and forecast horizonT=96T=96\. The correlation matrices are computed from the ground\-truth future series \(a\) and the future series predicted by iTransformer\[Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\]\(b\), while \(c\) shows the entries on which CvLoss yields smaller structural errors than iTransformer\(44%\)\.
Table 8:The comprehensive results on the long\-term forecasting task\.Note: We fix the input length as 96 following\[Liu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib1)\]\.Boldtypeface highlights the top performance for each metric, whileunderlinedtext denotes the second\-best results\.*Avg*indicates the results averaged over forecasting lengths: T=96, 192, 336 and 720\.

Table 9:Full results of short\-term forecasting\.Note: The input sequence length is set to 96 for all baselines\. All results are averaged across four different forecasting horizons:T∈\{12,24,48\}T\\in\\\{12,24,48\\\}\. The best and second\-best results are highlighted inboldandunderlined, respectively\.

Table 10:Backbone\-controlled comparison on four forecasting backbones: each pair contrasts the unmodified model with the identical model trained with CvLoss\. The*Wins*row counts, for each backbone, the horizon–dataset–metric cells on which CvLoss is better, out of3×4×2=243\\times 4\\times 2=24\. These counts are produced by the evaluation script from the*unrounded*outputs, not by comparing the three\-decimal values printed here: seven cells are equal at the displayed precision but are improvements at the fourth decimal, so the totals of24/2424/24per backbone and96/9696/96overall are exact\. We keep three decimals for consistency with the other tables rather than widening the display\.

## Appendix CMore Experimental Results

### C\.1Overall performance

Additional overall results are reported in Tables[8](https://arxiv.org/html/2608.05742#A2.T8)and[9](https://arxiv.org/html/2608.05742#A2.T9), where we provide performance under different forecast horizons\.

### C\.2Patch length sensitivity

Additional results for the patch length hyperparameter are reported in Figure[8](https://arxiv.org/html/2608.05742#A3.F8), covering the ETTh2 and ETTm2 datasets under the forecast horizonsT=192T=192andT=336T=336\.

![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/etth2_mse_legend_center.png)\(a\)ETTh2 with MSE
![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/etth2_mae_legend_center.png)\(b\)ETTh2 with MAE
![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/ettm2_mse_legend_center.png)\(c\)ETTm2 with MSE
![Refer to caption](https://arxiv.org/html/2608.05742v1/fig/ettm2_mae_legend_center.png)\(d\)ETTm2 with MAE

Figure 8:Sensitivity analysis of the patch length hyperparameter\. The plots illustrate the forecasting performance in terms of MSE and MAE on the ETTh2 \(a, b\) and ETTm2 \(c, d\) datasets\. The x\-axis represents the ratio of the forecast horizon to the patch length \(T / patch length\), while the lines denote different forecast horizons \(T=192T=192andT=336T=336\)\. The nearly horizontal curves demonstrate that the predictive accuracy remains consistent regardless of the specific patch length chosen, indicating that CvLoss is highly robust to this hyperparameter\.Table 11:Varyingα\\alpharesults of iTransformer- •Note:Boldandunderlineddenote the best and second\-best results\.

Table 12:Varyingα\\alpharesults of TimeBridge\.- •Note:Boldandunderlineddenote the best and second\-best results\.

### C\.3Comparison with different learning objectives

##### Baselines\.

Since this paper focuses on the design of learning objectives, we select competitive learning objectives tailored for training forecast models as baselines\. For the comprehensive evaluation in Table[14](https://arxiv.org/html/2608.05742#A3.T14), we compare against 11 objectives, which can be categorized as follows: ❶shape\-alignment objectives:GDTW\[Liu et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib11)\], Dilate\[Le Guen and Thome,[2019](https://arxiv.org/html/2608.05742#bib.bib12)\], and Soft\-DTW\[Cuturi and Blondel,[2017](https://arxiv.org/html/2608.05742#bib.bib13)\]; ❷likelihood maximization objectives:QDF\[Wang et al\.,[2026a](https://arxiv.org/html/2608.05742#bib.bib14)\], Time\-o1\[Wang et al\.,[2025b](https://arxiv.org/html/2608.05742#bib.bib15)\], Koopman\[Lange et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib16)\], FreDF\[Wang et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib6)\], and MSE \(also denoted as DF\); ❸distribution balancing objectives:KMB\-DF\[Pan et al\.,[2026](https://arxiv.org/html/2608.05742#bib.bib17)\], and DistDF\[Wang et al\.,[2026b](https://arxiv.org/html/2608.05742#bib.bib18)\]; and ❹decomposition\-based objectives:DBLoss\[Qiu et al\.,[2025b](https://arxiv.org/html/2608.05742#bib.bib19)\]\. For the cross\-model evaluation in Table[13](https://arxiv.org/html/2608.05742#A3.T13), we select a representative subset of 7 objectives: QDF, Time\-o1, FreDF, Koopman, Soft\-DTW, DF, and DBLoss\. The implementation of baselines follows the official codebase from\[Pan et al\.,[2026](https://arxiv.org/html/2608.05742#bib.bib17)\]\.

##### Implementation\.

To thoroughly assess the effectiveness of our proposed objective, we evaluate it under two distinct experimental settings\. In the first setting \(Table[14](https://arxiv.org/html/2608.05742#A3.T14)\), we employ CFPT\[Kou et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib41)\]as the default forecast model to compare all 10 learning objectives across six standard benchmarks \(ETTh1, ETTh2, ETTm1, ETTm2, Weather, and ECL\)\. In the second setting \(Table[13](https://arxiv.org/html/2608.05742#A3.T13)\), we adopt TQNet and PDF as the backbone models to verify the generalizability of our method across different architectures\. This evaluation is conducted on four datasets \(ETTh1, ETTm1, ECL, and Weather\) against the aforementioned subset of 6 learning objectives\. To ensure fair comparison, the drop\-last trick is disabled for all baselines, as recommended in\[Qiu et al\.,[2024](https://arxiv.org/html/2608.05742#bib.bib38)\]\. All objectives are trained with the Adam optimizer\[Kingma and Ba,[2015](https://arxiv.org/html/2608.05742#bib.bib39)\]\. When integrating CvLoss to train a forecast model, we retain all hyperparameters from the public benchmarks\[Kou et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib41)\], only tuningα\\alphaand the segment patch length\.

##### On the weakest competing objectives\.

Two of the eleven objectives, Koopman\[Lange et al\.,[2021](https://arxiv.org/html/2608.05742#bib.bib16)\]and Soft\-DTW\[Cuturi and Blondel,[2017](https://arxiv.org/html/2608.05742#bib.bib13)\], diverge sharply on some datasets\. Koopman reaches0\.4990\.499MSE on ETTm1 with TQNet, and Soft\-DTW reaches4\.5024\.502in Table[14](https://arxiv.org/html/2608.05742#A3.T14), so a reader may reasonably suspect a reproduction failure that would inflate our win counts\. It is not one\. Our figures reproduce the published numbers of the objective\-comparison protocol we adopt: on the PDF backbone, Table 2 of QDF\[Wang et al\.,[2026a](https://arxiv.org/html/2608.05742#bib.bib14)\]reports Koopman at0\.587/0\.4850\.587/0\.485on ETTm1 and Soft\-DTW at0\.695/0\.5480\.695/0\.548on ECL and1\.296/0\.4521\.296/0\.452on Weather, matching ours\. These older shape\-alignment objectives simply degrade on some multivariate benchmarks, as their own authors’ evaluations show\. None of our conclusions rests on them: neither ranks first in any setting, so removing both leaves the first\-place count of CvLoss in Table[14](https://arxiv.org/html/2608.05742#A3.T14)unchanged\. The objectives that actually compete with CvLoss are QDF, KMB\-DF, DBLoss, Time\-o1 and FreDF, and they differ from it by a few thousandths\.

Table 13:Comparable results with different learning objectives\.Table 14:Full results on the multi\-step forecasting task\. The length of history window is set to 96 for all baselines\.Avgindicates the results averaged over forecasting lengths: T=96, 192, 336 and 720 for ETT, ECL, and Weather\. Follow the settings of QDF\[Pan et al\.,[2026](https://arxiv.org/html/2608.05742#bib.bib17)\]

### C\.4Choice of the edge\-wise norm:ℓ1\\ell\_\{1\}versusℓ2\\ell\_\{2\}

Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)shows that MSE plus a*squared*graph total\-variation penalty is the negative log\-likelihood of a Gaussian with a structured precision matrix, whereas the objective we deploy in Eq\. \([12](https://arxiv.org/html/2608.05742#S3.E12)\) uses theℓ1\\ell\_\{1\}norm on the same edge differences\. The theory therefore backs theℓ2\\ell\_\{2\}form, and theℓ1\\ell\_\{1\}form is an empirical preference\. That asymmetry could hide a substantive weakness\. Ifℓ2\\ell\_\{2\}were as good, the theory\-backed variant would be available at no cost; ifℓ2\\ell\_\{2\}were far worse, the mechanism would not be what the theory describes\. We therefore test the two norms directly under an identical protocol, varying nothing but the norm\.

Table 15:Effect of the edge\-wise norm\. All settings are identical apart from the penalty:*Base*is the unmodified backbone,\+ℓ2\+\\ell\_\{2\}replaces∥⋅∥1\\\|\\cdot\\\|\_\{1\}in Eq\. \([12](https://arxiv.org/html/2608.05742#S3.E12)\) by the squaredℓ2\\ell\_\{2\}norm of Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2), and\+ℓ1\+\\ell\_\{1\}is CvLoss as deployed\. Values are averaged overT∈\{96,192,336,720\}T\\in\\\{96,192,336,720\\\}; the best entry in each row and metric is inbold\.Table[15](https://arxiv.org/html/2608.05742#A3.T15)reports the outcome\. Relative to the unregularised base, theℓ2\\ell\_\{2\}variant improves five of the eight averaged cells, matches the base in two, and slightly degrades Weather/TQNet MSE\. Relation\-aware coupling is therefore beneficial in its own right, which is the prediction Proposition[3\.2](https://arxiv.org/html/2608.05742#S3.Thmtheorem2)actually makes, and the benefit is not an artefact of the particular norm\. Theℓ1\\ell\_\{1\}variant is nonetheless better thanℓ2\\ell\_\{2\}in all eight cells\. We read this as consistent with the motivation given in Section[3\.2](https://arxiv.org/html/2608.05742#S3.SS2): a squared edge penalty is dominated by the few pairs with the largest structural mismatch, which in these series are typically localized anomalies rather than the systematic co\-evolution the term is meant to enforce\. The practical conclusion is that the theory identifies the right family and the norm is chosen within it on evidence, not thatℓ1\\ell\_\{1\}is derived\.

### C\.5Generalization studies

Additional experimental results of varying forecast models are available in Table[10](https://arxiv.org/html/2608.05742#A2.T10), Table[13](https://arxiv.org/html/2608.05742#A3.T13)and Fig[9](https://arxiv.org/html/2608.05742#A3.F9)\.

### C\.6Hyperparameter Sensitivity

We evaluate the stability ofCvLossregarding the loss weightα\\alphaand patch lengthLL\. Loss Weightα\\alpha: Tables[12](https://arxiv.org/html/2608.05742#A3.T12)and[12](https://arxiv.org/html/2608.05742#A3.T12)show accuracy consistently improving asα\\alphaincreases\. Optimal performance typically occurs atα≈0\.5\\alpha\\approx 0\.5or11, where structural regularization complements point\-wise accuracy\. Patch LengthLL: Figure[8](https://arxiv.org/html/2608.05742#A3.F8)demonstrates that forecasting performance remains stable across variousT/LT/Lratios\. The horizontal curves confirm thatCvLossis highly robust to patch size, requiring minimal tuning to achieve significant predictive gains\.

![Refer to caption](https://arxiv.org/html/2608.05742v1/x18.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x19.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x20.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x21.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x22.png)\(a\)ETTm1
![Refer to caption](https://arxiv.org/html/2608.05742v1/x23.png)\(b\)ETTm2
![Refer to caption](https://arxiv.org/html/2608.05742v1/x24.png)\(c\)ETTh1
![Refer to caption](https://arxiv.org/html/2608.05742v1/x25.png)\(d\)ETTh2

Figure 9:Performance of different forecast models with and without CvLossṪhe forecast errors are averaged over forecast lengths and the error bars represent 50% confidence intervals\.
### C\.7Case study with TQNet and TimeFilter of varying historical lengths

Table 16:Varying input sequence length results on the Weather dataset\.Additional experimental results of varying historical lengths are available in[Table 16](https://arxiv.org/html/2608.05742#A3.T16), complementing the fixed length of 96 used in the main text\. The forecast models selected include TimeFilter\[Hu et al\.,[2025a](https://arxiv.org/html/2608.05742#bib.bib4)\], which is the recent state\-of\-the\-art forecast model, and TQNet\[Lin et al\.,[2025](https://arxiv.org/html/2608.05742#bib.bib33)\]\. The results demonstrate that CvLoss consistently improves both forecast models across different historical sequence lengths\.

### C\.8Random Seed Sensitivity

Additional experimental results of random seed sensitivity are available in Table[17](https://arxiv.org/html/2608.05742#A3.T17), where we report the mean and standard deviation of results obtained from experiments conducted with seven different random seeds \(2020–2026\)\. The results indicate minimal sensitivity of the proposed method to random initialization, as all standard deviations remain below 0\.004\.

Table 17:CvLoss ablation on ECL and Weather\. Results are reported as mean±stdover seeds 2020–2026\. Lower is better\.
### C\.9Paired effect sizes

The standard deviations in Table[17](https://arxiv.org/html/2608.05742#A3.T17)are*marginal*: each summarises the spread of one arm across seeds\. Comparing two methods by asking whether their marginal intervals overlap is the wrong instrument here, and it is a conservative one, because the two arms are not independent samples\. Both are trained with the same seed indices, the same data ordering and the same initialisation scheme, so a large part of the seed\-to\-seed variance is shared and cancels in the difference\. This subsection therefore reports the paired quantity directly\.

For each cell and each seeds∈\{2020,…,2026\}s\\in\\\{2020,\\dots,2026\\\}we form the paired difference

ds=ℳbase,s−ℳC​v​L​o​s​s,s,d\_\{s\}=\\mathcal\{M\}\_\{\\mathrm\{base\},s\}\-\\mathcal\{M\}\_\{CvLoss,s\},\(14\)whereℳ\\mathcal\{M\}is the reported error metric, so thatds\>0d\_\{s\}\>0means CvLoss reduces the error on that seed\. Because both arms use the same seed indices, the mean paired difference coincides with the difference of the two reported seven\-seed means,d¯=17​∑sds=ℳ¯base−ℳ¯C​v​L​o​s​s\\bar\{d\}=\\frac\{1\}\{7\}\\sum\_\{s\}d\_\{s\}=\\bar\{\\mathcal\{M\}\}\_\{\\mathrm\{base\}\}\-\\bar\{\\mathcal\{M\}\}\_\{CvLoss\}\. We summarise each cell byd¯\\bar\{d\}, by the relative reductiond¯/ℳ¯base\\bar\{d\}/\\bar\{\\mathcal\{M\}\}\_\{\\mathrm\{base\}\}, and by the paired95%95\\%intervald¯±t0\.975,6​SD​\(ds\)/7\\bar\{d\}\\pm t\_\{0\.975,6\}\\,\\mathrm\{SD\}\(d\_\{s\}\)/\\sqrt\{7\}witht0\.975,6=2\.447t\_\{0\.975,6\}=2\.447, computed from the unrounded per\-seed outputs\.

Of the4040cells covered by Table[17](https://arxiv.org/html/2608.05742#A3.T17),2727are already separated by their marginal intervals\. Table[18](https://arxiv.org/html/2608.05742#A3.T18)reports the paired analysis for the remaining1313, which are exactly the cases the marginal criterion leaves unresolved\.

Table 18:Paired seven\-seed effects for the1313cells whose marginal intervals overlap in Table[17](https://arxiv.org/html/2608.05742#A3.T17)\.d¯\\bar\{d\}is the mean paired reduction in the stated metric,*Rel\.*the corresponding relative reduction, and the interval is the paired95%95\\%confidence interval ford¯\\bar\{d\}\. Positive values favour CvLoss\. No additional training was performed: these are the same seven completed runs re\-analysed in paired form\.All1313paired means favour CvLoss, and1111of the1313paired intervals exclude zero\. The two that do not, Weather/TQNet atT=336T=336andT=720T=720, both MSE, are positive but small, and we report them as inconclusive rather than as wins; they are the two settings in which the effect of CvLoss is genuinely at the resolution of the experiment\. Across all4040cells the median relative reduction is1\.66%1\.66\\%with an interquartile range of1\.131\.13percentage points, which we regard as the honest summary of the size of the effect on these two datasets\.

Two remarks on how these numbers should not be used\. First, we do not report a pooled significance test over the full grid of metric–horizon–dataset cells\. Such a test would treat the cells as independent observations, which they are not: MSE and MAE are computed from the same runs, and horizons within a dataset share the data, the split and the configuration\. A rank\-based test over the same cells would inherit the same dependence, so substituting one statistic for another would not repair the assumption\. Second, none of the paper’s claims requires it\. The evidence is the direction and consistency of the controlled comparisons themselves, summarised in Table[2](https://arxiv.org/html/2608.05742#S4.T2), together with the paired intervals above, and these are descriptive statements about the evaluated settings rather than inferences to a population of tasks\.

### C\.10Complexity

We evaluate the computational characteristics of CvLoss from two perspectives: \(i\) the training and inference overhead introduced by the original formulation, and \(ii\) the scalability improvement brought by random edge sampling\.

We first measure the runtime overhead of CvLoss on ETTh2 and ECL with the batch size fixed to 32\. As shown in Fig\.[10](https://arxiv.org/html/2608.05742#A3.F10), CvLoss introduces additional computation only during training, since it regularizes cross\-variable dependence through an extra loss term\. On the small\-scale ETTh2 benchmark, this overhead is marginal in both the forward and backward passes\. On the larger ECL benchmark, the overhead becomes more visible, especially in backpropagation, which is expected because modeling cross\-variable interactions is more expensive when the number of variables is large\. Importantly, although the prediction lengthTTvaries from 32 to 1024, the runtime remains in a similar range rather than increasing sharply, indicating that the added cost is well controlled across different horizons\.

At inference time, CvLoss is removed from the computation graph and therefore does not participate in prediction\. As shown in Fig\.[11](https://arxiv.org/html/2608.05742#A3.F11), the inference latency of the model trained with CvLoss remains highly comparable to that of the baseline iTransformer on both datasets\. The ratio plots further show that the latency fluctuates around parity, confirming that CvLoss introduces no practical deployment\-time overhead\.

While the above results demonstrate that CvLoss is inexpensive at inference and manageable during training, the original formulation still builds a fully connected cross\-variable graph\. When the number of patches isPPand the variable dimension isDD, the graph containsP×DP\\times Dnodes, and the number of pairwise interactions scales asO​\(P2​D2\)O\(P^\{2\}D^\{2\}\)\. This quadratic growth becomes a scalability bottleneck on high\-dimensional datasets such as ECL\. To address this issue, we further introduce*random edge sampling*, which limits the number of evaluated interactions to at most 1000 per batch\. In practice, this truncates the cost of CvLoss from rapidly increasing quadratic complexity to an effectively constant level\. We emphasise that random edge sampling is a scalability device studied in this subsection alone: it produced none of the accuracy numbers reported elsewhere in the paper, all of which use the complete cross\-variable graph without sampling or selection\.

Fig\.[12](https://arxiv.org/html/2608.05742#A3.F12)compares the training runtime before and after sampling\. On ETTh2, where the variable dimension is small, the full CvLoss remains tractable, but the sampled version still makes the runtime flatter as the number of patches increases\. On ECL, the benefit is much more pronounced: the runtime of the full CvLoss grows rapidly with the number of patches, whereas the sampled version remains almost unchanged\. This confirms that random edge sampling effectively removes the main scalability bottleneck of the original formulation\.

Fig\.[13](https://arxiv.org/html/2608.05742#A3.F13)further shows that random edge sampling substantially improves training efficiency without sacrificing predictive performance\. As the number of patches increases, the speedup grows rapidly, while the relative MSE change remains close to zero\.

Fig\.[12](https://arxiv.org/html/2608.05742#A3.F12)summarizes this observation from an efficiency–accuracy perspective\. The results show that large computational gains can be achieved with negligible forecasting degradation, suggesting that a small number of randomly sampled interactions is sufficient to preserve the regularization effect of CvLoss in practice\.

![Refer to caption](https://arxiv.org/html/2608.05742v1/x26.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x27.png)

\(a\)ETTh2: forward \(left\) and backward \(right\)\.
![Refer to caption](https://arxiv.org/html/2608.05742v1/x28.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x29.png)

\(b\)ECL: forward \(left\) and backward \(right\)\.

Figure 10:Training latency on ETTh2 and ECL\. CvLoss adds little overhead on ETTh2 and a larger, mainly backward\-pass overhead on ECL\.![Refer to caption](https://arxiv.org/html/2608.05742v1/x30.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x31.png)

\(a\)ETTh2: inference latency \(left\) and latency ratio \(right\)\.
![Refer to caption](https://arxiv.org/html/2608.05742v1/x32.png)

![Refer to caption](https://arxiv.org/html/2608.05742v1/x33.png)

\(b\)ECL: inference latency \(left\) and latency ratio \(right\)\.

Figure 11:Inference latency on ETTh2 and ECL\. CvLoss remains comparable to the baseline, with ratios close to parity\.![Refer to caption](https://arxiv.org/html/2608.05742v1/x34.png)\(a\)ETTh2 runtime\.
![Refer to caption](https://arxiv.org/html/2608.05742v1/x35.png)\(b\)ECL runtime\.

Figure 12:Training runtime versus the number of patches\. Dashed lines denote random edge sampling with at most 1000 edges\.![Refer to caption](https://arxiv.org/html/2608.05742v1/x36.png)\(a\)ECL speedup\.
![Refer to caption](https://arxiv.org/html/2608.05742v1/x37.png)\(b\)MSE change\.
![Refer to caption](https://arxiv.org/html/2608.05742v1/x38.png)\(c\)Trade\-off\.

Figure 13:Optimization benefit of random edge sampling\. The speedup grows rapidly with the number of patches, while the relative MSE change remains close to zero\. The trade\-off plot summarizes that large speedups are achieved with negligible performance degradation\.Table 19:Comprehensive long\-term forecasting results grouped by backbone models in Benchmark TS\-Library\[Wu et al\.,[2023](https://arxiv.org/html/2608.05742#bib.bib25)\]\.Note: Average results \(*Avg*\) are calculated overT∈\{96,192,336,720\}T\\in\\\{96,192,336,720\\\}\. Evaluation for thebestandsecond\-bestresults are independently conducted within theiTransformer, theDLinear, and thePatchTST\. "Ours" refers to different enhancements applied to the respective base models\.

## Appendix DComplexity Analysis and Edge Selection via Structural Discrepancy

In the main text, we introduced CvLoss as a lightweight plug\-in loss function that regularizes the structural consistency among variables\. While a fully connected cross\-variable graph effectively captures comprehensive spatiotemporal interactions, calculating the loss over all possible variable pairs yields a quadratic computational complexity of𝒪​\(P2​D2\)\\mathcal\{O\}\(P^\{2\}D^\{2\}\)\. This poses a significant scalability challenge for high\-dimensional multivariate systems\. To address this bottleneck, we introduce a targeted edge selection mechanism based on structural discrepancy to significantly reduce the computational complexity\.

##### Status of this mechanism\.

Three edge constructions appear in this paper, and it matters which produced which result\. The complete cross\-variable graph of Section[3\.2](https://arxiv.org/html/2608.05742#S3.SS2), with no sampling and no selection, produced*every*accuracy number in the main text and in the appendix\. The random edge sampling of Appendix[C\.10](https://arxiv.org/html/2608.05742#A3.SS10)and the top\-KKselection described below are scalability devices for theO​\(P2​D2\)O\(P^\{2\}D^\{2\}\)edge count, and are evaluated only in the efficiency studies of Figures[12](https://arxiv.org/html/2608.05742#A3.F12)and[13](https://arxiv.org/html/2608.05742#A3.F13)\. The synchronous and asynchronous supports compared in Section[4\.5](https://arxiv.org/html/2608.05742#S4.SS5)are ablations of the topology, not the default\. The distinction is important for one further reason: the top\-KKrule below is the only construction in the paper that reads the ground\-truth correlation matrix𝐂​\(Y\)\\mathbf\{C\}\(Y\), and it does so during training only, on the training split\. No configuration uses any statistic of the future window at test time\.

### D\.1Targeting Maximum Structural Discrepancy

As established in our theoretical analysis, the standard MSE fundamentally ignores latent cross\-variable relationships in its design\. Consequently, models trained purely with MSE often produce predictions with significant structural discrepancies\. This phenomenon is visually evident in the cross\-variable correlation heatmaps presented in our motivating examples \(e\.g\., Figure[1\(a\)](https://arxiv.org/html/2608.05742#S2.F1.sf1)\)\.

The heatmaps intuitively represent the structural relationships between variables, where each point corresponds to the correlation between a specific variable pair\. By comparing the predicted correlation heatmap𝐂​\(Y^\)\\mathbf\{C\}\(\\hat\{Y\}\)with the ground\-truth correlation heatmap𝐂​\(Y\)\\mathbf\{C\}\(Y\), we observe clear differences indicating where the model fails to capture the true co\-evolutionary dynamics\.

To optimize computational efficiency, we do not need to compute the first\-order difference loss over all points in the heatmap\. Instead, we can selectively compute the loss exclusively for those variable pairs that exhibit the largest discrepancy between the prediction and the ground truth\. Specifically, we measure the structural error for each variable pair\(i,j\)\(i,j\)asΔi​j=\|𝐂​\(Y^\)i​j−𝐂​\(Y\)i​j\|\\Delta\_\{ij\}=\|\\mathbf\{C\}\(\\hat\{Y\}\)\_\{ij\}\-\\mathbf\{C\}\(Y\)\_\{ij\}\|\. We then select the top\-KKpairs with the maximumΔi​j\\Delta\_\{ij\}to construct our edge setℰ\\mathcal\{E\}\.

By calculating the CvLoss strictly on these most severely misaligned pairs, we force the model to correct its most critical structural failures\. This targeted approach dramatically reduces the computational complexity from𝒪​\(D2\)\\mathcal\{O\}\(D^\{2\}\)to a manageable𝒪​\(K\)\\mathcal\{O\}\(K\), demonstrating that CvLoss can be efficiently scaled to high\-dimensional multivariate time series forecasting without suffering from quadratic computational overhead\.

### D\.2Diagnostic study: does the objective track cross\-variable structure?

This subsection is a diagnostic and does not describe the deployed method\.Everywhere else in this paper, the two loss terms are combined by the fixed convex weightα∈\(0,1\)\\alpha\\in\(0,1\)of Eq\. \([13](https://arxiv.org/html/2608.05742#S3.E13)\), selected on validation as described in Appendix[B\.2](https://arxiv.org/html/2608.05742#A2.SS2.SSS0.Px1)\. We therefore ask a separate question: if the relative importance of the structural term were free to adapt, would it track how much cross\-variable structure a dataset actually contains? To answer it, we run an auxiliary experiment in which the two terms carry*unnormalised*coefficients learned during training\. The quantities in Table[20](https://arxiv.org/html/2608.05742#A4.T20)are those auxiliary coefficients\. They are not values ofα\\alpha, they are not constrained to\(0,1\)\(0,1\), and they produced no result reported anywhere else in the paper; we report them only for the correspondence described below\.

With that caveat, Table[20](https://arxiv.org/html/2608.05742#A4.T20)places the converged coefficients alongside the statistical properties of each dataset, specifically the average Concordance Correlation Coefficient \(CCC\) and Pearson Correlation Coefficient\.

We observe a clear positive correspondence between the channel\-wise correlation \(CCC/Pearson\), the dimensionality \(DD\), and the converged coefficient on the structural term\. For datasets with fewer variables and lower cross\-channel correlation, such as the ETT series \(D=7D=7, CCC≈0\.14\\approx 0\.14–0\.180\.18\), the optimization converges to relatively small and comparable coefficients for both MSE and the structural penalty, indicating a balanced reliance\.

Conversely, for high\-dimensional spatial\-temporal graphs exhibiting strong internal consistency, such as the PEMS traffic datasets \(D≥170D\\geq 170, CCC\>0\.77\>0\.77\) and Solar \(D=137D=137, CCC\>0\.91\>0\.91\), the converged coefficient on the structural term grows substantially\. This is the behaviour Theorem[3\.1](https://arxiv.org/html/2608.05742#S3.Thmtheorem1)predicts: the objective gap it identifies is larger precisely when the residual precision matrix is further from spherical, so a dataset with a rich, highly synchronized cross\-variable geometry has more for the structural term to correct\. The same pattern is visible in the deployed method, where the largest accuracy gains occur on the most correlated datasets \(PEMS07, CCC0\.780\.78,−11\.3%\-11\.3\\%MSE; Solar, CCC0\.910\.91,−4\.4%\-4\.4\\%\) and the smallest on the least correlated \(the ETT series, CCC0\.140\.14–0\.180\.18\)\. The practical implication is that channel correlation can be measured before training to anticipate whether the term will help\.

Table 20:Diagnostic study only — these are not values ofα\\alpha\.Converged coefficients of an auxiliary experiment in which the two loss terms carry*unnormalised*learned weights, reported alongside the channel\-wise correlation of each dataset\. They are unconstrained in magnitude, are not the convex weightα∈\(0,1\)\\alpha\\in\(0,1\)of Eq\. \([13](https://arxiv.org/html/2608.05742#S3.E13)\), and produced none of the accuracy results in this paper; see Appendix[D\.2](https://arxiv.org/html/2608.05742#A4.SS2)\. The datasets are ordered by their number of variables \(DD\)\.

## Appendix EBroader Impacts

CvLoss is a general training objective for multivariate time series forecasting\. Its potential positive impact is to improve forecasting reliability in applications such as energy management, traffic monitoring, weather analysis, and other operational systems where variables evolve jointly\. Better forecasts may support more efficient resource allocation and planning\. We do not identify direct negative societal impacts specific to CvLoss, since the method does not introduce new sensitive data, human\-subject data, generative capabilities, or deployment\-specific decision rules\. Any societal risk would mainly come from downstream misuse or misinterpretation of forecasting models in application domains, rather than from the proposed loss itself\.

## Appendix FStatement on the Use of Large Language Models \(LLMs\)

In accordance with the conference guidelines, we disclose our use of Large Language Models \(LLMs\) in the preparation of this paper as follows:

We used LLMs \(specifically, OpenAI GPT\-5\.2, GPT\-5\.4 and Google Gemini 3\)*solely for checking grammar errors and improving the readability of the manuscript*\. The LLMs*were not involved in research ideation, the development of research contributions, experiment design, data analysis, or interpretation of results*\. All substantive content and scientific claims were created entirely by the authors\. The authors have reviewed all LLM\-assisted text to ensure accuracy and originality, and take full responsibility for the contents of the paper\. The LLMs are not listed as an author\.