时空 Kronecker 协方差神经网络

# Spatiotemporal Kronecker Covariance Neural Networks
Source: [https://arxiv.org/html/2609.25326](https://arxiv.org/html/2609.25326)
\{IEEEkeywords\}

Covariance Neural Networks, Spatiotemporal Principal Component Analysis

Student MemberIEEEAthanasios Georgoutsos and Elvin IsufiSenior MemberIEEEAffiliation:Delft University of Technology, Delft, Netherlands

###### Abstract

Multivariate time series contain complex patterns that span across both space and time\. While covariance\-based statistical tools like spatiotemporal Principal Component Analysis \(ST\-PCA\) help identify these patterns, they are limited to linear operations and prone to estimation errors with limited data\. Recent covariance\-based spatiotemporal neural networks offer more stable, non\-linear alternatives, but they ignore correlations across different time steps\. To solve this, we introduce the Kronecker coVariance Neural Network \(KVNN\), a temporal graph neural network that represents the spatiotemporal covariance matrix via a sum of Kronecker products where spatial and temporal dependencies are decoupled\. By implementing filtering operations on spatial and temporal components, KVNNs achieve expressive processing capabilities, admit a rigorous spectral analysis, and are provably stable to finite\-sample estimation errors, ultimately addressing all of ST\-PCA’s limitations\. We show on five real\-world datasets that KVNNs achieve strong forecasting performance, often requiring significantly fewer trainable parameters than competitive methods, and are consistent under estimation noise\.

††corresponding:Corresponding author: Andrea Cavallo \(email: a\.cavallo@tudelft\.nl\)\.††note:This work was supported in part by the TU Delft AI Labs programme, NWO OTP GraSPA proposal \#19497, NWO VENI proposal 222\.032, and by the SURE\-AI Centre grant \#357482, Research Council of Norway\. This work is not related to Athanasios Georgoutsos’ position at Amazon\.com Inc\.## 1INTRODUCTION

Multivariate time series exhibit intricate spatiotemporal correlations whose successful estimation and processing is key for forecasting tasks as well as for modeling and understanding spatiotemporal behaviors\[[21](https://arxiv.org/html/2609.25326#bib.bib14),[14](https://arxiv.org/html/2609.25326#bib.bib15),[17](https://arxiv.org/html/2609.25326#bib.bib17),[11](https://arxiv.org/html/2609.25326#bib.bib16),[13](https://arxiv.org/html/2609.25326#bib.bib18)\]\. ST\-PCA offers a prominent tool to identify spatiotemporal patterns, enhance understanding and tractability of multivariate time series\[[22](https://arxiv.org/html/2609.25326#bib.bib19),[18](https://arxiv.org/html/2609.25326#bib.bib20)\]\. Usually, ST\-PCA relies on the eigendecomposition of a large spatiotemporal covariance matrix, either built directly from an extended space\-time covariance matrix\[[16](https://arxiv.org/html/2609.25326#bib.bib22),[18](https://arxiv.org/html/2609.25326#bib.bib20)\], or on decoupled parametrizations across space and time\. It has found applications in traffic flow forecasting\[[22](https://arxiv.org/html/2609.25326#bib.bib19)\], galaxy dynamics\[[35](https://arxiv.org/html/2609.25326#bib.bib21)\], and online change\-point detection\[[1](https://arxiv.org/html/2609.25326#bib.bib23)\], among others\. Despite its informativeness and widespread use, ST\-PCA suffers from two crucial limitations\. First, it has limited expressivity, since it only offers linear processing of covariance information, and is task\-agnostic\. Second, it is sensitive to finite\-sample errors, which prevents stable performance across repeated trials with different sampling noise\.

A covariance\-based alternative to enhance expressivity and stability are spatiotemporal coVariance Neural Networks \(STVNNs\)\[[6](https://arxiv.org/html/2609.25326#bib.bib3)\]; temporal graph neural networks operating on sample covariance as a proxy of a graph structure of the data\. These are based on the VNN principle\[[29](https://arxiv.org/html/2609.25326#bib.bib4),[28](https://arxiv.org/html/2609.25326#bib.bib39)\], where the covariance matrix of the data is used as an inductive bias in a graph neural network trained end\-to\-end to solve a downstream task, and extend VNNs to include dynamic information via a temporal window of past data\. In this way, STVNNs are highly expressive via non\-linear covariance processing and are provably more stable to finite\-sample estimation errors than PCA\. However, they only consider covariances in individual snapshots, thus ignoring temporally delayed correlations\. This neglects crucial spatiotemporal information and weakens the connection with ST\-PCA\.

To address these limitations, we introduce the Kronecker coVariance Neural Network \(KVNN\), a temporal graph neural network that operates on the spatiotemporal covariance matrix as a graph shift operator, thus accounting for correlations across space and time, and decomposes such covariance as a sum of Kronecker terms, representing spatial and temporal connections separately\. This covariance representation allows decoupled processing of different temporal delays, overcoming STVNN’s limitations\. To further improve efficiency and interpretability, we integrate a group\-sparsity penalty during training, enabling automatic pruning of uninformative spatiotemporal blocks for the task at hand\.

We propose two practical implementations of KVNNs\. The first one relies on fixed temporal matrices and estimated covariances at fixed delays, which admits a rigorous analysis of their spectral behavior and stability to finite\-sample errors\. The second one, instead, uses a low\-rank Kronecker covariance estimate to model long\-term dependencies with limited Kronecker terms, thus increasing parameter efficiency\.

Related works\.The proposed approach could be seen as a deep learning solution to process multivariate time series using the spatiotemporal covariance as an inductive bias\. Alternative deep learning\-based methods for this task also exist, including recurrent models\[[7](https://arxiv.org/html/2609.25326#bib.bib8),[19](https://arxiv.org/html/2609.25326#bib.bib7)\], graph\-time neural networks\[[20](https://arxiv.org/html/2609.25326#bib.bib27),[37](https://arxiv.org/html/2609.25326#bib.bib26),[36](https://arxiv.org/html/2609.25326#bib.bib25),[26](https://arxiv.org/html/2609.25326#bib.bib28),[15](https://arxiv.org/html/2609.25326#bib.bib29)\]and transformers\[[40](https://arxiv.org/html/2609.25326#bib.bib24),[39](https://arxiv.org/html/2609.25326#bib.bib10),[38](https://arxiv.org/html/2609.25326#bib.bib9)\]\. However, these methods either ignore the spatial structure, assume a known graph structure, or they learn it via free parameters or generic attention scores rather than employing an explicit statistical prior, missing a connection to spatiotemporal PCA and its interpretability and stability guarantees\.

The VNN principle that our method relies on was introduced in\[[29](https://arxiv.org/html/2609.25326#bib.bib4),[28](https://arxiv.org/html/2609.25326#bib.bib39)\]and has received increasing attention in methodological extensions\[[4](https://arxiv.org/html/2609.25326#bib.bib6),[5](https://arxiv.org/html/2609.25326#bib.bib34),[25](https://arxiv.org/html/2609.25326#bib.bib35)\]and applications\[[30](https://arxiv.org/html/2609.25326#bib.bib2),[31](https://arxiv.org/html/2609.25326#bib.bib5),[32](https://arxiv.org/html/2609.25326#bib.bib40)\]\. STVNNs\[[6](https://arxiv.org/html/2609.25326#bib.bib3)\]apply this principle to temporal data, but ignore delayed correlations\. We proposed a preliminary version of KVNNs, termed lagged VNNs \(LVNNs\), in\[[10](https://arxiv.org/html/2609.25326#bib.bib12)\]\. These, however, lack expressivity to separately process delayed correlations, and are computationally inefficient as they do not consider Kronecker covariance representations\.

Contributions\.Our contributions are threefold:

\(C1\)We propose the KVNN, a neural network that operates on the Kronecker representation of the spatiotemporal covariance matrix, processing lagged covariances in an expressive and efficient manner\.

\(C2\)We provide a spectral analysis of KVNN, and we prove that its asymptotic stability against finite\-sample estimation errors improves over ST\-PCA\.

\(C3\)We provide extensive numerical experiments on five real\-world datasets, demonstrating that KVNNs achieve highly competitive forecasting performance while benefiting from the parameter efficiency of group\-sparse and low\-rank representations\.

𝐂0,0\{\\lx@inpgf@ignorespaces\\mathbf\{C\}\_\{0,0\}\}𝐂0,1\{\\lx@inpgf@ignorespaces\\mathbf\{C\}\_\{0,1\}\}…\{\\lx@inpgf@ignorespaces\\dots\}𝐂0,T−1\{\\lx@inpgf@ignorespaces\\mathbf\{C\}\_\{0,T\-1\}\}𝐂1,0\{\\lx@inpgf@ignorespaces\\mathbf\{C\}\_\{1,0\}\}𝐂1,1\{\\lx@inpgf@ignorespaces\\mathbf\{C\}\_\{1,1\}\}…\{\\lx@inpgf@ignorespaces\\dots\}𝐂1,T−1\{\\lx@inpgf@ignorespaces\\mathbf\{C\}\_\{1,T\-1\}\}⋮\{\\lx@inpgf@ignorespaces\\vdots\}⋮\{\\lx@inpgf@ignorespaces\\vdots\}⋱\{\\lx@inpgf@ignorespaces\\ddots\}⋮\{\\lx@inpgf@ignorespaces\\vdots\}𝐂T−1,0\{\\lx@inpgf@ignorespaces\\mathbf\{C\}\_\{T\-1,0\}\}𝐂T−1,1\{\\lx@inpgf@ignorespaces\\mathbf\{C\}\_\{T\-1,1\}\}…\{\\lx@inpgf@ignorespaces\\dots\}𝐂T−1,T−1\{\\lx@inpgf@ignorespaces\\mathbf\{C\}\_\{T\-1,T\-1\}\}𝐂~∈ℝN​T×N​T\\tilde\{\\mathbf\{C\}\}\\in\\mathbb\{R\}^\{NT\\times NT\}TTblocks==∑r=1R\\displaystyle\\sum\_\{r=1\}^\{R\}a00\{\\lx@inpgf@ignorespaces a\_\{00\}\}a01\{\\lx@inpgf@ignorespaces a\_\{01\}\}…\{\\lx@inpgf@ignorespaces\\dots\}a10\{\\lx@inpgf@ignorespaces a\_\{10\}\}a11\{\\lx@inpgf@ignorespaces a\_\{11\}\}…\{\\lx@inpgf@ignorespaces\\dots\}⋮\{\\lx@inpgf@ignorespaces\\vdots\}⋮\{\\lx@inpgf@ignorespaces\\vdots\}⋱\{\\lx@inpgf@ignorespaces\\ddots\}𝐀r\\mathbf\{A\}\_\{r\}\(temporal\)T×TT\\times T⊗\\otimes𝐁r\\mathbf\{B\}\_\{r\}\(spatial\)N×NN\\times N→\\toa00​𝐁r\{\\lx@inpgf@ignorespaces a\_\{00\}\\mathbf\{B\}\_\{r\}\}a01​𝐁r\{\\lx@inpgf@ignorespaces a\_\{01\}\\mathbf\{B\}\_\{r\}\}…\{\\lx@inpgf@ignorespaces\\dots\}a10​𝐁r\{\\lx@inpgf@ignorespaces a\_\{10\}\\mathbf\{B\}\_\{r\}\}a11​𝐁r\{\\lx@inpgf@ignorespaces a\_\{11\}\\mathbf\{B\}\_\{r\}\}…\{\\lx@inpgf@ignorespaces\\dots\}⋮\{\\lx@inpgf@ignorespaces\\vdots\}⋮\{\\lx@inpgf@ignorespaces\\vdots\}⋱\{\\lx@inpgf@ignorespaces\\ddots\}block structureTTblocks

t−2t\-2t−1t\-1tt\[𝐱t−2\]1\[\\mathbf\{x\}\_\{t\-2\}\]\_\{1\}\[𝐱t−2\]2\[\\mathbf\{x\}\_\{t\-2\}\]\_\{2\}\[𝐱t−2\]3\[\\mathbf\{x\}\_\{t\-2\}\]\_\{3\}\[𝐱t−1\]1\[\\mathbf\{x\}\_\{t\-1\}\]\_\{1\}\[𝐱t−1\]2\[\\mathbf\{x\}\_\{t\-1\}\]\_\{2\}\[𝐱t−1\]3\[\\mathbf\{x\}\_\{t\-1\}\]\_\{3\}\[𝐱t\]1\[\\mathbf\{x\}\_\{t\}\]\_\{1\}\[𝐱t\]2\[\\mathbf\{x\}\_\{t\}\]\_\{2\}\[𝐱t\]3\[\\mathbf\{x\}\_\{t\}\]\_\{3\}\[𝐂0\]23\[\\mathbf\{C\}\_\{0\}\]\_\{23\}\[𝐂1\]22\[\\mathbf\{C\}\_\{1\}\]\_\{22\}\[𝐂1\]12\[\\mathbf\{C\}\_\{1\}\]\_\{12\}\[𝐂2\]11\[\\mathbf\{C\}\_\{2\}\]\_\{11\}

Figure 1:\(Left\) Kronecker representation of the spatiotemporal covariance𝐂~\\tilde\{\\mathbf\{C\}\}as a sum ofRRdecoupled interactions\. The temporal covariance𝐀r\\mathbf\{A\}\_\{r\}scales the spatial covariance𝐁r\\mathbf\{B\}\_\{r\}to form the block structure\. \(Right\) Spatiotemporal covariance graph on a multivariate time series\. Nodeiiat timettcarries the feature\[𝐱t\]i\[\\mathbf\{x\}\_\{t\}\]\_\{i\}, edges connect nodes within a snapshot \(green\), and across snapshots \(red\)\. Edge weights are the covariances\[𝐂τ\]i​j\[\\mathbf\{C\}\_\{\\tau\}\]\_\{ij\}at lagτ\\tau\.
## 2SPATIOTEMPORAL COVARIANCES

Consider a data matrix𝐗∈ℝN×M′\\mathbf\{X\}\\in\\mathbb\{R\}^\{N\\times M^\{\\prime\}\}where each column𝐱t∈ℝN\\mathbf\{x\}\_\{t\}\\in\\mathbb\{R\}^\{N\}is an observation of a zero\-mean multivariate time series and the observations are in temporal order\. To model the joint spatiotemporal correlations in𝐗\\mathbf\{X\}, we consider a temporal window of sizeTTand we stackTTconsecutive temporal samples into a unique vector𝐱~t=\[𝐱t−T\+1𝖳,𝐱t−T\+2𝖳,…,𝐱t𝖳\]𝖳\\mathbf\{\\tilde\{x\}\}\_\{t\}=\[\\mathbf\{x\}\_\{t\-T\+1\}^\{\\mathsf\{T\}\},\\mathbf\{x\}\_\{t\-T\+2\}^\{\\mathsf\{T\}\},\\dots,\\mathbf\{x\}\_\{t\}^\{\\mathsf\{T\}\}\]^\{\\mathsf\{T\}\}\. By rearranging the data matrix𝐗\\mathbf\{X\}in this way, we get the matrix𝐗~∈ℝN​T×M\\mathbf\{\\tilde\{X\}\}\\in\\mathbb\{R\}^\{NT\\times M\}withM=M′/TM=M^\{\\prime\}/T, which collects observations of the variable𝐱~∈ℝN​T\\mathbf\{\\tilde\{x\}\}\\in\\mathbb\{R\}^\{NT\}with covariance matrix

𝐂~0=𝔼⁡\[𝐱~​𝐱~𝖳\]∈ℝN​T×N​T\.\\mathbf\{\\tilde\{C\}\}\_\{0\}=\\mathbb\{E\}\[\\mathbf\{\\tilde\{x\}\}\\mathbf\{\\tilde\{x\}\}^\{\\mathsf\{T\}\}\]\\in\\mathbb\{R\}^\{NT\\times NT\}\.\(1\)𝐂~0\\mathbf\{\\tilde\{C\}\}\_\{0\}encodes correlations across space and time via a block structure containingT×TT\\times Tblocks𝐂t1,t2=𝔼⁡\[𝐱t1​𝐱t2𝖳\]∈ℝN×N\\mathbf\{C\}\_\{t\_\{1\},t\_\{2\}\}=\\mathbb\{E\}\[\\mathbf\{x\}\_\{t\_\{1\}\}\\mathbf\{x\}\_\{t\_\{2\}\}^\{\\mathsf\{T\}\}\]\\in\\mathbb\{R\}^\{N\\times N\}fort1,t2=0,…​T−1t\_\{1\},t\_\{2\}=0,\\dots T\-1, i\.e\., the covariance among samples at time pointst1,t2t\_\{1\},t\_\{2\}\(cf\. Figure[1](https://arxiv.org/html/2609.25326#S1.F1), left\)\.

𝐂~0\\mathbf\{\\tilde\{C\}\}\_\{0\}can be directly estimated from the samples by stackingTTconsecutive observations and using the sample estimator

𝐂~s=1M​∑t=1M𝐱~t​𝐱~t𝖳\.\\displaystyle\\mathbf\{\\tilde\{C\}\}\_\{s\}=\\frac\{1\}\{M\}\\sum\_\{t=1\}^\{M\}\\mathbf\{\\tilde\{x\}\}\_\{t\}\\mathbf\{\\tilde\{x\}\}\_\{t\}^\{\\mathsf\{T\}\}\.\(2\)However,𝐂~s\\mathbf\{\\tilde\{C\}\}\_\{s\}is high\-dimensional, thus requiring the estimation of a large number of parameters, which may lead to poor results if insufficient data is available, and causes data processing techniques relying on this matrix \(e\.g\., spatiotemporal PCA\) to be computationally inefficient\. Moreover,𝐂~s\\mathbf\{\\tilde\{C\}\}\_\{s\}jointly models spatial and temporal correlations, making it difficult to interpret the resulting principal components\.

Considering its spatiotemporal structure, the matrix𝐂~0\\mathbf\{\\tilde\{C\}\}\_\{0\}can be approximated as a sum ofRRKronecker products

𝐂~=∑r=1R𝐀r⊗𝐁r\\mathbf\{\\tilde\{C\}\}=\\sum\_\{r=1\}^\{R\}\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}\(3\)where matrices𝐀r∈ℝT×T\\mathbf\{A\}\_\{r\}\\in\\mathbb\{R\}^\{T\\times T\}encode temporal relations and matrices𝐁r∈ℝN×N\\mathbf\{B\}\_\{r\}\\in\\mathbb\{R\}^\{N\\times N\}spatial correlations \(cf\. Figure[1](https://arxiv.org/html/2609.25326#S1.F1), left\)\. This representation decouples the spatial and the temporal contributions, allowing for better interpretability and various estimators based on the data assumptions\. In general,𝐂~≠𝐂~0\\mathbf\{\\tilde\{C\}\}\\neq\\mathbf\{\\tilde\{C\}\}\_\{0\}, but for specific choices of Kronecker terms, equality holds\. We present different choices of𝐀r,𝐁r\\mathbf\{A\}\_\{r\},\\mathbf\{B\}\_\{r\}below\.

### 2\.1WIDE\-SENSE STATIONARY DATA

If data come from a stationary distribution, cross\-covariances depend only on the temporal distance between two observations,𝐂t1,t2=𝐂t1−t2=𝐂τ\\mathbf\{C\}\_\{t\_\{1\},t\_\{2\}\}=\\mathbf\{C\}\_\{t\_\{1\}\-t\_\{2\}\}=\\mathbf\{C\}\_\{\\tau\}with𝐂−τ=𝐂τ𝖳\\mathbf\{C\}\_\{\-\\tau\}=\\mathbf\{C\}\_\{\\tau\}^\{\\mathsf\{T\}\}, so all blocks on the same diagonal of𝐂~0\\mathbf\{\\tilde\{C\}\}\_\{0\}in \([1](https://arxiv.org/html/2609.25326#S2.E1)\) coincide\. Since𝐂τ≠𝐂τ𝖳\\mathbf\{C\}\_\{\\tau\}\\neq\\mathbf\{C\}\_\{\\tau\}^\{\\mathsf\{T\}\}in general, the two diagonals at distanceτ\\taucarry different blocks\. To maintain \(skew\-\)symmetric Kronecker terms, we setR=2​T−1R=2T\-1, and we split every lagged covariance𝐂τ\\mathbf\{C\}\_\{\\tau\}into its symmetric and skew\-symmetric parts,𝐂τs=12​\(𝐂τ\+𝐂τ𝖳\)\\mathbf\{C\}\_\{\\tau\}^\{s\}=\\tfrac\{1\}\{2\}\\big\(\\mathbf\{C\}\_\{\\tau\}\+\\mathbf\{C\}\_\{\\tau\}^\{\\mathsf\{T\}\}\\big\), and𝐂τa=12​\(𝐂τ−𝐂τ𝖳\)\\mathbf\{C\}\_\{\\tau\}^\{a\}=\\tfrac\{1\}\{2\}\\big\(\\mathbf\{C\}\_\{\\tau\}\-\\mathbf\{C\}\_\{\\tau\}^\{\\mathsf\{T\}\}\\big\), such that𝐂τ=𝐂τs\+𝐂τa\\mathbf\{C\}\_\{\\tau\}=\\mathbf\{C\}\_\{\\tau\}^\{s\}\+\\mathbf\{C\}\_\{\\tau\}^\{a\}\. The two terms carry different information\. The symmetric part𝐂τs\\mathbf\{C\}\_\{\\tau\}^\{s\}measures how strongly spatial patterns co\-moveτ\\tausteps apart; the skew part𝐂τa\\mathbf\{C\}\_\{\\tau\}^\{a\}measures which of two patterns moves first, i\.e\., the lead–lag structure at delayτ\\tau, and vanishes whenever the lagged covariance is symmetric\.

Then, we use as Kronecker pairs𝐀r,𝐁r\\mathbf\{A\}\_\{r\},\\mathbf\{B\}\_\{r\}the matrixes\(𝐈T,𝐂^0\)\(\\mathbf\{I\}\_\{T\},\\mathbf\{\\hat\{C\}\}\_\{0\}\)and, forτ=1,…,T−1\\tau=1,\\dots,T\-1,

\(𝐃τ\+𝐃τ𝖳,𝐂^τs\)and\(𝐃τ𝖳−𝐃τ,𝐂^τa\),\\big\(\\mathbf\{D\}\_\{\\tau\}\+\\mathbf\{D\}\_\{\\tau\}^\{\\mathsf\{T\}\},\\;\\mathbf\{\\hat\{C\}\}\_\{\\tau\}^\{s\}\\big\)\\quad\\text\{and\}\\quad\\big\(\\mathbf\{D\}\_\{\\tau\}^\{\\mathsf\{T\}\}\-\\mathbf\{D\}\_\{\\tau\},\\;\\mathbf\{\\hat\{C\}\}\_\{\\tau\}^\{a\}\\big\),\(4\)where𝐃τ∈ℝT×T\\mathbf\{D\}\_\{\\tau\}\\in\\mathbb\{R\}^\{T\\times T\}is defined as

\[𝐃τ\]i​j=\{1,if​j−i=τ,0,otherwise,\[\\mathbf\{D\}\_\{\\tau\}\]\_\{ij\}=\\begin\{cases\}1,&\\textnormal\{if \}j\-i=\\tau,\\\\ 0,&\\textnormal\{otherwise,\}\\end\{cases\}\(5\)and𝐂^τ\\mathbf\{\\hat\{C\}\}\_\{\\tau\}is the sample estimate of𝐂τ\\mathbf\{C\}\_\{\\tau\}\. This is an exact rewriting of the two diagonals at distanceτ\\tau, since\(𝐃τ\+𝐃τ𝖳\)⊗𝐂τs\+\(𝐃τ𝖳−𝐃τ\)⊗𝐂τa=𝐃τ⊗𝐂τ𝖳\+𝐃τ𝖳⊗𝐂τ\(\\mathbf\{D\}\_\{\\tau\}\+\\mathbf\{D\}\_\{\\tau\}^\{\\mathsf\{T\}\}\)\\otimes\\mathbf\{C\}\_\{\\tau\}^\{s\}\+\(\\mathbf\{D\}\_\{\\tau\}^\{\\mathsf\{T\}\}\-\\mathbf\{D\}\_\{\\tau\}\)\\otimes\\mathbf\{C\}\_\{\\tau\}^\{a\}=\\mathbf\{D\}\_\{\\tau\}\\otimes\\mathbf\{C\}\_\{\\tau\}^\{\\mathsf\{T\}\}\+\\mathbf\{D\}\_\{\\tau\}^\{\\mathsf\{T\}\}\\otimes\\mathbf\{C\}\_\{\\tau\}, and it requires to estimate only theTTspatial matrices𝐂^0,…,𝐂^T−1\\mathbf\{\\hat\{C\}\}\_\{0\},\\dots,\\mathbf\{\\hat\{C\}\}\_\{T\-1\}\. All temporal matrices in \([4](https://arxiv.org/html/2609.25326#S2.E4)\) are either symmetric or skew\-symmetric, and so are the spatial ones\.

### 2\.2LOW SEPARATION\-RANK COVARIANCE

Another estimator is proposed in\[[33](https://arxiv.org/html/2609.25326#bib.bib11)\], which allows to reduce the number of Kronecker termsRReven further if the covariance is separable in few Kronecker components\. This relies on a truncated SVD of the permuted sample covariance𝐑~=𝒫⁡\(𝐂~s\)∈ℝT2×N2\\mathbf\{\\tilde\{R\}\}=\\mathcal\{P\}\(\\mathbf\{\\tilde\{C\}\}\_\{s\}\)\\in\\mathbb\{R\}^\{T^\{2\}\\times N^\{2\}\}, where𝒫\\mathcal\{P\}is a permutation operator mapping the entries of a matrix of shapeN​T×N​TNT\\times NTto a matrix of shapeT2×N2T^\{2\}\\times N^\{2\}\. Given the SVD𝐑~=𝐔​𝚺​𝐕𝖳\\mathbf\{\\tilde\{R\}\}=\\mathbf\{U\}\\mathbf\{\\Sigma\}\\mathbf\{V\}^\{\\mathsf\{T\}\}, the firstRRleft and right singular vector pairs𝐮r∈ℝT2,𝐯r∈ℝN2\\mathbf\{u\}\_\{r\}\\in\\mathbb\{R\}^\{T^\{2\}\},\\mathbf\{v\}\_\{r\}\\in\\mathbb\{R\}^\{N^\{2\}\},r=1,…,Rr=1,\\dots,Rare reshaped to be the corresponding terms𝐀r,𝐁r\\mathbf\{A\}\_\{r\},\\mathbf\{B\}\_\{r\}, following the property that𝒫⁡\(𝐀r⊗𝐁r\)=𝐚r​𝐛r𝖳\\mathcal\{P\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}\)=\\mathbf\{a\}\_\{r\}\\mathbf\{b\}\_\{r\}^\{\\mathsf\{T\}\}where𝐚r∈ℝT2\\mathbf\{a\}\_\{r\}\\in\\mathbb\{R\}^\{T^\{2\}\},𝐛r∈ℝN2\\mathbf\{b\}\_\{r\}\\in\\mathbb\{R\}^\{N^\{2\}\}\. This allows processing arbitrarily large temporal dependenciesTTwhile still maintaining few terms viaRR\.

## 3KRONECKER COVARIANCE NEURAL NETWORKS

Following the VNN principle\[[29](https://arxiv.org/html/2609.25326#bib.bib4),[28](https://arxiv.org/html/2609.25326#bib.bib39)\], we use spatiotemporal covariance information to build a graph, where the features of samples𝐱t\\mathbf\{x\}\_\{t\}are node signals, and the edges are covariances at different lags \(cf\. Figure[1](https://arxiv.org/html/2609.25326#S1.F1), right\)\. Then, we propose the spatiotemporal Kronecker coVariance Neural Network \(KVNN\), a graph convolutional neural network that operates on such graph, represented via the Kronecker decomposition in \([3](https://arxiv.org/html/2609.25326#S2.E3)\)\. We first define the Kronecker coVariance Filter \(KVF\), which processes an input signal𝐱~\\mathbf\{\\tilde\{x\}\}to produce an output𝐮\\mathbf\{u\}as:

𝐮=∑r=1R𝐇r​\(𝐀r,𝐁r\)​𝐱~=∑r=1R∑k=0Khk​r​𝐀rk⊗𝐁rk​𝐱~,\\mathbf\{u\}=\\sum\_\{r=1\}^\{R\}\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\},\\mathbf\{B\}\_\{r\}\)\\mathbf\{\\tilde\{x\}\}=\\sum\_\{r=1\}^\{R\}\\sum\_\{k=0\}^\{K\}h\_\{kr\}\\mathbf\{A\}\_\{r\}^\{k\}\\otimes\\mathbf\{B\}\_\{r\}^\{k\}\\mathbf\{\\tilde\{x\}\},\(6\)where we used the property\(𝐀⊗𝐁\)k=𝐀k⊗𝐁k\(\\mathbf\{A\}\\otimes\\mathbf\{B\}\)^\{k\}=\\mathbf\{A\}^\{k\}\\otimes\\mathbf\{B\}^\{k\}\. The filter in \([6](https://arxiv.org/html/2609.25326#S3.E6)\) learns different coefficients\{hk​r\}k=0K\\\{h\_\{kr\}\\\}\_\{k=0\}^\{K\}for each Kronecker termrr, thus gaining the flexibility to distinctly process correlation information at different time delays\. To show the inherent flexibility of this approach, we compare it with the potential alternative filter design𝐮=∑k=0Khk​\(∑r=1R𝐀r⊗𝐁r\)k​𝐱~\\mathbf\{u\}=\\sum\_\{k=0\}^\{K\}h\_\{k\}\\left\(\\sum\_\{r=1\}^\{R\}\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}\\right\)^\{k\}\\mathbf\{\\tilde\{x\}\}, adopted by our preliminary work LVNN\[[10](https://arxiv.org/html/2609.25326#bib.bib12)\]\. The latter learns a unique set of coefficientshkh\_\{k\}for all Kronecker terms, losing expressivity, and the covariance powerskkmix the various Kronecker terms𝐀r⊗𝐁r\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}, preventing a separate processing of spatiotemporal information at different time delays\.

A KVNN is a sequence of layers followed by non\-linear activation functionsσ\\sigma, where each layerℓ=1,…,L\\ell=1,\\dots,Lcontains a filterbank ofFin×FoutF\_\{\\textnormal\{in\}\}\\times F\_\{\\textnormal\{out\}\}KVFs processing the input in parallel:

𝐱~f\(ℓ\)=σ\(∑j=1Fin𝐇f​j\(ℓ\)𝐱~j\(ℓ−1\)\),f=1,…,Fout\.\\mathbf\{\\tilde\{x\}\}\_\{f\}^\{\(\\ell\)\}=\\sigma\\left\(\\sum\_\{j=1\}^\{F\_\{\\textnormal\{in\}\}\}\\mathbf\{H\}\_\{fj\}^\{\(\\ell\)\}\\mathbf\{\\tilde\{x\}\}\_\{j\}^\{\(\\ell\-1\)\}\\right\),\\quad f=1,\\dots,F\_\{\\textnormal\{out\}\}\.\(7\)We denote the complete KVNN architecture withΦ⁡\(ℋ,𝐂,𝐱~\)\\Phi\(\\mathcal\{H\},\\mathbf\{C\},\\mathbf\{\\tilde\{x\}\}\)for an input signal𝐱~\\mathbf\{\\tilde\{x\}\}, a covariance representation𝐂\\mathbf\{C\}and withℋ=\{hk​r​f​j\(ℓ\)\}k​r​f​j​ℓ\\mathcal\{H\}=\\\{h\_\{krfj\}^\{\(\\ell\)\}\\\}\_\{krfj\\ell\}the learnable coefficients\.

Implementations\.We define two KVNN implementations which use different Kronecker covariance estimators to achieve different desiderata\.KVNN\-Suses the stationary estimator in Section[2\.1](https://arxiv.org/html/2609.25326#S2.SS1), which allows for theoretical tractability of the frequency response and stability to finite\-sample estimations\.KVNN\-LR, instead, uses the low\-rank estimation in Section[2\.2](https://arxiv.org/html/2609.25326#S2.SS2)\. This decouples the temporal windowTTfrom the number of terms to estimateRR, yielding parameter\-efficient estimators even for large temporal dependencies\. However, it prevents tractability of the spectral and stability analysis\.

Learning sparse weights\.In practice, not all blocks in𝐂~\\mathbf\{\\tilde\{C\}\}might be helpful for the downstream task, while they may introduce computational burden\. Thus, we introduce a group sparsity penalty in the training loss such that the KVNN could ignore covariance terms that are not useful for prediction\. Letℋ=\{hk​r​f​j\}k​r​f​j\\mathcal\{H\}=\\\{h\_\{krfj\}\\\}\_\{krfj\}be the set of all learnable coefficients of a KVNN\. Given a training set\(𝐗tr,𝐘tr\)\(\\mathbf\{X\}\_\{\\textnormal\{tr\}\},\\mathbf\{Y\}\_\{\\textnormal\{tr\}\}\), we train the KVNN to minimize the following objective

minℋ⁡ℒtask​\(𝐘tr,Φ⁡\(ℋ,𝐂,𝐗tr\)\)\+λg​g​\(ℋ\)\\min\_\{\\mathcal\{H\}\}\\mathcal\{L\}\_\{\\textnormal\{task\}\}\(\\mathbf\{Y\}\_\{\\textnormal\{tr\}\},\\Phi\(\\mathcal\{H\},\\mathbf\{C\},\\mathbf\{X\}\_\{\\textnormal\{tr\}\}\)\)\+\\lambda\_\{g\}g\(\\mathcal\{H\}\)\(8\)whereℒtask\\mathcal\{L\}\_\{\\textnormal\{task\}\}is a task\-specific loss \(e\.g\., forecasting MSE\) andg⁡\(ℋ\)g\(\\mathcal\{H\}\)is a group\-sparsity penalty defined asg⁡\(ℋ\)=∑ℓ=1L∑r=1R\(∑k,j,f\(hk​r​f​j\(ℓ\)\)2\)1/2,g\(\\mathcal\{H\}\)=\\sum\_\{\\ell=1\}^\{L\}\\sum\_\{r=1\}^\{R\}\\left\(\\sum\_\{k,j,f\}\(h\_\{krfj\}^\{\(\\ell\)\}\)^\{2\}\\right\)^\{1/2\},which promotes sparsity of entire blocks of coefficientshk​r​f​j\(ℓ\)h\_\{krfj\}^\{\(\\ell\)\}for a given value ofrr, thus ignoring some covariance block terms\. In practice, during training we minimizeg⁡\(ℋ\)g\(\\mathcal\{H\}\)by applying a proximal soft\-thresholding step after each gradient update\. Forλg\>0\\lambda\_\{g\}\>0we additionally prune, before evaluating on validation, a copy of the model in which we set to zero the coefficientshk​r​f​j\(ℓ\)h^\{\(\\ell\)\}\_\{krfj\}whose norm over\(k,f,j\)\(k,f,j\)at a givenrrfalls belowα\\alpha\(set to0\.10\.1\) times the layer’s mean group norm\.

![Refer to caption](https://arxiv.org/html/2609.25326v1/fig_raster.png)Figure 2:Spectrum of the symmetric and asymmetric lagged covariances on the synthetic time series in \(Row 1\)\. \(Row 2\) The largest eigenvector of𝐂τs\\mathbf\{C\}\_\{\\tau\}^\{s\}is non\-zero on cells with persistent patterns, while the last is non\-zero on cells with highly alternating patterns at delayτ\\tau\. \(Row 3\) The eigenvector pair𝐚,𝐛\\mathbf\{a\},\\mathbf\{b\}of𝐂τa\\mathbf\{C\}\_\{\\tau\}^\{a\}is non\-zero for cells that have correlated lead\-follow patterns at delayτ\\tau\.
## 4THEORETICAL ANALYSIS

We present the spectral analysis of KVNN\-S \(along with its filter KVF\-S\), and its stability to finite\-sample errors\. Proofs are collected in the Appendix\.

### 4\.1SPECTRAL ANALYSIS

A unified spectral analysis of the KVF\-S is challenging because the Kronecker terms have different eigenvectors\. However, since the KVF\-S learns different coefficients on distinct Kronecker terms, we analyze each term separately\. This characterizes the filter expressivity and processing of spatiotemporal information and connects it to classical spectral tools\.

For a genericrr, let𝐀r=𝐔r​𝚲r​𝐔r𝖧\\mathbf\{A\}\_\{r\}=\\mathbf\{U\}\_\{r\}\\mathbf\{\\Lambda\}\_\{r\}\\mathbf\{U\}\_\{r\}^\{\\mathsf\{H\}\}and𝐁r=𝐕r​𝐖r​𝐕r𝖧\\mathbf\{B\}\_\{r\}=\\mathbf\{V\}\_\{r\}\\mathbf\{W\}\_\{r\}\\mathbf\{V\}\_\{r\}^\{\\mathsf\{H\}\}be the unitary eigendecompositions of the two factors\. When𝐀r,𝐁r\\mathbf\{A\}\_\{r\},\\mathbf\{B\}\_\{r\}are symmetric,𝐔r,𝐕r\\mathbf\{U\}\_\{r\},\\mathbf\{V\}\_\{r\}are real orthogonal and𝚲r,𝐖r\\mathbf\{\\Lambda\}\_\{r\},\\mathbf\{W\}\_\{r\}are real; when they are skew\-symmetric,𝐔r,𝐕r\\mathbf\{U\}\_\{r\},\\mathbf\{V\}\_\{r\}are complex unitary and𝚲r,𝐖r\\mathbf\{\\Lambda\}\_\{r\},\\mathbf\{W\}\_\{r\}are purely imaginary\. In both cases the productsλr,m​wr,n\\lambda\_\{r,m\}w\_\{r,n\}are real, and the term𝐀r⊗𝐁r\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}is a real symmetric matrix\. We define therr\-thspatiotemporal covariance Fourier transformfor a KVF\-S and a signal𝐱~\\mathbf\{\\tilde\{x\}\}as𝐱^=\(𝐔r⊗𝐕r\)𝖧​𝐱~\\mathbf\{\\hat\{x\}\}=\(\\mathbf\{U\}\_\{r\}\\otimes\\mathbf\{V\}\_\{r\}\)^\{\\mathsf\{H\}\}\\mathbf\{\\tilde\{x\}\}\. For eachrr\(or equivalentlyτ\\tau\), this eigenspectrum has an interpretable meaning\. The temporal factors𝐃τ\+𝐃τ𝖳\\mathbf\{D\}\_\{\\tau\}\+\\mathbf\{D\}\_\{\\tau\}^\{\\mathsf\{T\}\}and𝐃τ𝖳−𝐃τ\\mathbf\{D\}\_\{\\tau\}^\{\\mathsf\{T\}\}\-\\mathbf\{D\}\_\{\\tau\}are theτ\\tau\-step undirected line\-graph shift and its skew counterpart\. Forτ=1\\tau=1their eigenvectors are the discrete sinusoids that graph\-time signal processing adopts as temporal frequencies\[[27](https://arxiv.org/html/2609.25326#bib.bib38)\]\. The spatial factor𝐂τs\\mathbf\{C\}\_\{\\tau\}^\{s\}has as eigenvectors spatial patterns and as eigenvalues their lag\-τ\\tauautocovariance, i\.e\., patterns are ranked by how strongly they carry overτ\\tausteps\. These are the basis of time\-lagged component analysis\[[3](https://arxiv.org/html/2609.25326#bib.bib32),[23](https://arxiv.org/html/2609.25326#bib.bib33)\]\. The spatial factor𝐂τa\\mathbf\{C\}\_\{\\tau\}^\{a\}has eigenvaluesw=±i​μw=\\pm i\\muattached to pairs of orthogonal patterns\(𝐚,𝐛\)\(\\mathbf\{a\},\\mathbf\{b\}\), whereμ\\mumeasures their co\-variation at delayτ\\tau\. Figure[2](https://arxiv.org/html/2609.25326#S3.F2)shows an example\. The first and last eigenvectors of𝐂1s\\mathbf\{C\}\_\{1\}^\{s\}identify, respectively, persistent signals and patterns alternating everyτ=1\\tau=1steps\. The largest eigenpair of𝐂1a\\mathbf\{C\}\_\{1\}^\{a\}, instead, identifies a pattern where one signal leads another one byτ=1\\tau=1step\.

We now analyze the KVF\-S in therr\-th eigenspectrum\. By properties of the Kronecker product, we get

𝐮r\\displaystyle\\mathbf\{u\}\_\{r\}=∑k=0Khk​r​\(𝐔r⊗𝐕r\)​\(𝚲rk⊗𝐖rk\)​\(𝐔r⊗𝐕r\)𝖧​𝐱~\\displaystyle=\\sum\_\{k=0\}^\{K\}h\_\{kr\}\(\\mathbf\{U\}\_\{r\}\\otimes\\mathbf\{V\}\_\{r\}\)\(\\mathbf\{\\Lambda\}\_\{r\}^\{k\}\\otimes\\mathbf\{W\}\_\{r\}^\{k\}\)\(\\mathbf\{U\}\_\{r\}\\otimes\\mathbf\{V\}\_\{r\}\)^\{\\mathsf\{H\}\}\\mathbf\{\\tilde\{x\}\}\(9\)=\(𝐔r⊗𝐕r\)​∑k=0Khk​r​\(𝚲rk⊗𝐖rk\)​𝐱^\.\\displaystyle=\(\\mathbf\{U\}\_\{r\}\\otimes\\mathbf\{V\}\_\{r\}\)\\sum\_\{k=0\}^\{K\}h\_\{kr\}\(\\mathbf\{\\Lambda\}\_\{r\}^\{k\}\\otimes\\mathbf\{W\}\_\{r\}^\{k\}\)\\mathbf\{\\hat\{x\}\}\.\(10\)For theii\-th entry of𝐮\\mathbf\{u\}, indexed by a temporal–spatial pair\(m,n\)\(m,n\), this corresponds to

\[𝐮r\]i=∑k=0Khk​r​λr,mk​wr,nk​\[𝐱^\]i=hr​\(λr,m​wr,n\)​\[𝐱^\]i\.\\displaystyle\[\\mathbf\{u\}\_\{r\}\]\_\{i\}=\\sum\_\{k=0\}^\{K\}h\_\{kr\}\\,\\lambda\_\{r,m\}^\{k\}w\_\{r,n\}^\{k\}\\,\[\\mathbf\{\\hat\{x\}\}\]\_\{i\}=h\_\{r\}\(\\lambda\_\{r,m\}w\_\{r,n\}\)\\,\[\\mathbf\{\\hat\{x\}\}\]\_\{i\}\.\(11\)That is, the KVF\-S learns a polynomial frequency responsehr​\(γ\)h\_\{r\}\(\\gamma\)in the variableγ=λ​w\\gamma=\\lambda w, the product of a temporal and a spatial eigenvalue of therr\-th Kronecker term \(cf\. Figure[3](https://arxiv.org/html/2609.25326#S4.F3)\)\.

![Refer to caption](https://arxiv.org/html/2609.25326v1/fig_response3d.png)Figure 3:Frequency responsehr​\(λ​w\)h\_\{r\}\(\\lambda w\)of the KVF\-S on therr\-th Kronecker term\.
### 4\.2STABILITY ANALYSIS

We analyze the stability of the KVF\-S to finite\-sample estimation errors\. The temporal matrices𝐀r\\mathbf\{A\}\_\{r\}are fixed, so the finite\-sample estimation error lies only in the spatial matrices\. We present the Lipschitz filter assumption, then state our stability result\.

###### Assumption 1

The frequency responsehrh\_\{r\}in \([11](https://arxiv.org/html/2609.25326#S4.E11)\) is Lipschitz with constantPPon an intervalℐr⊂ℝ\\mathcal\{I\}\_\{r\}\\subset\\mathbb\{R\}containing the eigenvalues of both𝐀r⊗𝐁r\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}and𝐀r⊗𝐁^r\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{\\hat\{B\}\}\_\{r\}, i\.e\.,\|hr​\(γ\)−hr​\(γ′\)\|≤P​\|γ−γ′\|\|h\_\{r\}\(\\gamma\)\-h\_\{r\}\(\\gamma^\{\\prime\}\)\|\\leq P\|\\gamma\-\\gamma^\{\\prime\}\|for allγ,γ′∈ℐr\\gamma,\\gamma^\{\\prime\}\\in\\mathcal\{I\}\_\{r\}\.

###### Theorem 1

Consider a KVF satisfying Assumption[1](https://arxiv.org/html/2609.25326#Thmassumption1)and operating on the stationary Kronecker terms \([4](https://arxiv.org/html/2609.25326#S2.E4)\), whose spatial matrices are estimated fromM\>N​TM\>NTindependent, zero\-mean, sub\-Gaussian observations of a wide\-sense stationary process, stacked as in \([2](https://arxiv.org/html/2609.25326#S2.E2)\)\. Let𝐮0\\mathbf\{u\}\_\{0\}be the KVF output on the true lagged covariances and𝐮\\mathbf\{u\}the output on their sample estimates\. Then

‖𝐮0−𝐮‖2=𝒪P​\(T3/2​P​NM\)​‖𝐱~‖2,\\\|\\mathbf\{u\}\_\{0\}\-\\mathbf\{u\}\\\|\_\{2\}=\\mathcal\{O\}\_\{P\}\\\!\\left\(\\frac\{T^\{3/2\}PN\}\{\\sqrt\{M\}\}\\right\)\\\|\\mathbf\{\\tilde\{x\}\}\\\|\_\{2\},\(12\)where𝒪P\\mathcal\{O\}\_\{P\}indicates that the bound holds with probability increasing to11asM→∞M\\to\\infty\.

The stability of the KVF depends on its Lipschitz constantPP: a largerPPallows sharper variations of the response across spatiotemporal modes, which in turn can amplify errors in the estimated spectrum\. The filter converges to the ideal one at rate𝒪\(M−1/2\)\\mathcal\{O\}\(M^\{\-1/2\}\), like standard covariance estimators, and the bound grows withNNandTTas these increase the number of covariance entries to estimate\. This result can be extended to the complete KVNN by applying\[[8](https://arxiv.org/html/2609.25326#bib.bib36), Theorem 4\]\. Notably, the KVNN stability bound does not increase if covariance eigenvalues are close\. This happens, instead, for ST\-PCA, which requires an explicit eigendecomposition \(cf\.\[[6](https://arxiv.org/html/2609.25326#bib.bib3), Proposition 1\]\)\. The advantage of KVNN lies in its spectral processing via polynomial functions bounded by the Lipschitz constant, which ensure that eigenvalue perturbations result in bounded output differences\. This analysis extends analogous observations in the parallel between VNNs and PCA\[[29](https://arxiv.org/html/2609.25326#bib.bib4)\]\. Finally, while independence of the windows is an idealization for stationary data, most of our stability proof does not use it, and it can be replaced by a finite\-memory or spectral\-density condition on the process at the cost of constants \(cf\.\[[2](https://arxiv.org/html/2609.25326#bib.bib37)\]\)\.

Table 1:Test MAE \(×10\\times 10, mean±\\pmstd, lower is better\) by prediction horizonhh\.Bestandsecond bestresults highlighted\.Table 2:Number of learnable parameters for each model in the configuration reported in Table[1](https://arxiv.org/html/2609.25326#S4.T1)\.

## 5NUMERICAL EXPERIMENTS

Our numerical experiments have the following objectives\.\(O1\)Demonstrate that KVNNs are effective at real\-data forecasting and their design choices are motivated by performance improvements w\.r\.t\. relevant baselines\.\(O2\)Show that our stability analysis translates into consistent KVNNs’ performance under finite\-sample covariance estimation errors\.\(O3\)Support the group\-sparsity penalty formulation showing that it identifies meaningful temporal lags\.\(O4\)Analyze the impact of hyperparameters on performance\. Our code is available online111[https://github\.com/andrea\-cavallo\-98/KVNN](https://github.com/andrea-cavallo-98/KVNN)\.

### 5\.1FORECASTING \(O1\)

Datasets\.We use five real\-world multivariate time series\. Three come from the TFB forecasting benchmark\[[24](https://arxiv.org/html/2609.25326#bib.bib13)\]\.NN5containsN=111N=111daily series of cash withdrawals at UK ATMs over791791days, with a weekly seasonality\.CzeLancontainsN=11N=11hourly ecological channels recorded at a Czech forest site over19,93419\{,\}934steps\.AQWancontainsN=11N=11hourly air\-quality channels from a Beijing monitoring station over35,06435\{,\}064steps\. The remaining two are variance forecasting tasks that we construct from publicly available daily price data\. The input collects the daily log returns ofNNassets, and the target is the vector of their within\-day variances on future days, computed with the Garman–Klass estimator\[[9](https://arxiv.org/html/2609.25326#bib.bib30)\]\.Sectorscomprises theN=9N=9original Select Sector SPDR exchange\-traded funds over6,9616\{,\}961trading days \(1999–2026\)\.LargeCapcomprisesN=50N=50large\-capitalization U\.S\. stocks over1,4251\{,\}425days \(2021–2026\)\.

We split every series chronologically into60%60\\%training,20%20\\%validation and20%20\\%test, and standardized with training statistics\. We use CzeLan and AQWan in first differences,𝐱t←𝐱t−𝐱t−1\\mathbf\{x\}\_\{t\}\\leftarrow\\mathbf\{x\}\_\{t\}\-\\mathbf\{x\}\_\{t\-1\}, since in levels they are dominated by persistence, we truncate them and to their first10,00010\{,\}000steps\. We extract windows ofTTconsecutive samples at unit stride\.

Models\.We evaluate three models of the proposed family:KVNN\-S, which uses the stationary Kronecker estimator;KVNN\-LR, which uses the low\-rank estimator; andLVNN\[[10](https://arxiv.org/html/2609.25326#bib.bib12)\], which operates on the full𝐂~s\\mathbf\{\\tilde\{C\}\}\_\{s\}with coefficients shared across lags\. Each baseline isolates one design choice of KVNNs:VNN\[[29](https://arxiv.org/html/2609.25326#bib.bib4)\]ignores temporal information;STVNN\[[6](https://arxiv.org/html/2609.25326#bib.bib3)\]uses a temporal window but neglects delayed correlations;ST\-PCAprojects the windowed input onto the topqqeigenvectors of𝐂~s\\mathbf\{\\tilde\{C\}\}\_\{s\}before a forecasting head, i\.e\., it processes spatiotemporal covariance disjointly from the task; and anLSTM\[[12](https://arxiv.org/html/2609.25326#bib.bib31)\]uses no covariance information\. Finally,Ref\.is the best of naive predictors: the training mean and, where available, persistence and seasonal naive\.

Experimental setup\.We forecast all horizons jointly,h∈\{1,2,3\}h\\in\\\{1,2,3\\\}for NN5, Sectors and LargeCap andh∈\{1,3,6\}h\\in\\\{1,3,6\\\}for CzeLan and AQWan, and report the test MAE on standardized targets as mean±\\pmstandard deviation over33seeds\. We optimize hyperparameters on validation set via a grid search among the following values:TT\(three dataset\-specific values\), depthL∈\{1,2,3\}L\\in\\\{1,2,3\\\}and filter orderK∈\{1,2,3\}K\\in\\\{1,2,3\\\}; in additionλg∈\{0,0\.3,1\}\\lambda\_\{g\}\\in\\\{0,0\.3,1\\\}for KVNN\-S and KVNN\-LR,R∈\{1,3,8\}R\\in\\\{1,3,8\\\}for KVNN\-LR,q∈\{8,32,128\}q\\in\\\{8,32,128\\\}for ST\-PCA, and a width ladder for each baseline\. For KVNN\-S and KVNN\-LR, we fixF=32F=32, while for each baseline we increase the feature size until its largest configuration is at least as large as KVNN’s largest\. All remaining settings are shared by every model and configuration: Adam with learning rate10−210^\{\-2\}, batch size128128, dropout0\.10\.1, a 2\-layer MLP head, a softmax temporal readout, an autoregressive skip initialized to exact persistence, Kronecker blocks normalized to unit spectral norm, and at most600600epochs with early stopping on validation MAE after4040epochs without improvement\. Table[1](https://arxiv.org/html/2609.25326#S4.T1)reports the results, and Table[2](https://arxiv.org/html/2609.25326#S4.T2)the number of trainable parameters of the selected configurations\.

Discussion\.At least one member of the KVNN family \(LVNN, KVNN\-S, KVNN\-LR\) ranks in the top two in all1515dataset–horizon settings; the family occupies both first and second place in1111of them and provides the best model in1212\. The three covariance\-based baselines highlight the value of each design choice: \(i\) the gains over ST\-PCA show that end\-to\-end, task\-aware covariance processing outperforms a disjoint project\-then\-forecast pipeline; \(ii\) the gains over VNN show that temporal information is essential; \(iii\) the gains over STVNN show that delayed cross\-covariances carry information that the snapshot covariance alone misses\.

Among the proposed models, KVNN\-LR is the most consistent: it is the best model in77of1515settings and does so with fewer parameters than KVNN\-S and LVNN on every dataset \(Table[2](https://arxiv.org/html/2609.25326#S4.T2)\), e\.g\.,1\.31\.3k against2222k and1010k on Sectors\. This indicates that a low\-rank Kronecker representation captures the relevant delayed structure more parsimoniously than the explicit lag\-by\-lag parametrization of KVNN\-S\. LVNN is competitive but less consistent, since it cannot weigh different lags differently\.

Every KVNN model improves on Ref\., confirming that the tasks are non\-trivial\. The LSTM, which uses no covariance information, is the least consistent learned model: it barely improves on Ref\. on LargeCap and is among the weakest models on AQWan, yet it is the best model on NN5, where it outperforms KVNN\-S by a small margin, but uses more than three times the parameter count\. This supports the view that the covariance provides a useful inductive bias: it substantially reduces the number of trainable parameters at the cost of, at most, a minor loss in accuracy\.

Figure 4:Stability to finite\-sample covariance estimation errors on the five real\-world datasets\.![Refer to caption](https://arxiv.org/html/2609.25326v1/hyperparams_ss.png)Figure 5:\(Left\) Share of the retained coefficient mass per lag \(λg=2\\lambda\_\{g\}=2\), for a process whose only non\-zero lagged covariance is atτ=q\\tau=q\(rows\)\. The group\-sparsity penalty concentrates the surviving coefficients on the informative lag\. \(Middle\) Impact of hyperparametersT,RT,R\(RRin the legend\)\. \(Right\) Impact of hyperparametersλg,R\\lambda\_\{g\},R\.
### 5\.2STABILITY \(O2\)

We test empirically the finite\-sample stability predicted by Theorem[1](https://arxiv.org/html/2609.25326#Thmtheorem1)\. For each dataset, KVNN\-S, KVNN\-LR and ST\-PCA are trained once with the covariance operator estimated from the full training series and then frozen\. We then re\-estimate the covariance from a contiguous block ofMMconsecutive training samples with a uniformly random start, substitute it into the frozen model, and measure the test MAE\. Since nothing but the covariance operator changes, the resulting degradation isolates the effect of estimation error on the operator from differences in what was learned\. For eachMM, we recompute the covariance2020times with different data subsets and report the mean results±1\\pm 1standard deviation\. All models share one configuration: a single layer withF=32F=32,K=1K=1, an MLP head of sizes\(32,16\)\(32,16\), dropout0\.10\.1, Adam with learning rate10−210^\{\-2\}, batch size128128, and at most300300epochs with early stopping after4040\. We setλg=0\\lambda\_\{g\}=0so that no Kronecker term is pruned, which would alter the operator under study\. KVNN\-S \(R=2​T−1R=2T\-1\) is shown for the three values ofTTin each dataset’s grid; KVNN\-LR \(R=8R=8\) and ST\-PCA \(q=32q=32\) are shown at the intermediate valueTrefT\_\{\\textnormal\{ref\}\}\.

Figure[4](https://arxiv.org/html/2609.25326#S5.F4)shows that KVNNs are markedly more stable than ST\-PCA\. The test MAE of KVNNs stays close to its full\-sample value down to very small fractions of the training data and decreases monotonically asMMgrows, whereas ST\-PCA is uniformly worse, exhibits a larger variance across blocks, and on NN5 degrades sharply whenMMis small\. Within KVNN\-S, a largerTTleads to a larger gap between the small\-MMand full\-sample MAE, consistent with the dependence onTTin Theorem[1](https://arxiv.org/html/2609.25326#Thmtheorem1), although on LargeCap and Sectors the curves are essentially flat and the effect is negligible\. KVNN\-LR is less stable than KVNN\-S due to the explicit SVD required to estimate the Kronecker terms\. Neveretheless, it provides more consistent performance compared to ST\-PCA\.

### 5\.3LAG RECOVERY \(O3\)

We verify that the penaltyg⁡\(ℋ\)g\(\\mathcal\{H\}\)removes the Kronecker terms that carry no correlation in a controlled setting\. We generate a moving\-average process𝐱t=𝜺t\+0\.9​𝐐​𝜺t−q\\mathbf\{x\}\_\{t\}=\\boldsymbol\{\\varepsilon\}\_\{t\}\+0\.9\\,\\mathbf\{Q\}\\,\\boldsymbol\{\\varepsilon\}\_\{t\-q\}with𝜺t∼𝒩⁡\(𝟎,𝐈N\)\\boldsymbol\{\\varepsilon\}\_\{t\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\mathbf\{I\}\_\{N\}\),𝐐\\mathbf\{Q\}a symmetric orthogonal matrix andN=8N=8, so that every𝐂τ\\mathbf\{C\}\_\{\\tau\}is symmetric and each lag is carried by a single Kronecker term\. By construction,𝐂τ=𝟎\\mathbf\{C\}\_\{\\tau\}=\\mathbf\{0\}forτ∉\{0,q\}\\tau\\notin\\\{0,q\\\}, so the corresponding blocks𝐁r\\mathbf\{B\}\_\{r\}, estimated from𝐂^τ\\mathbf\{\\hat\{C\}\}\_\{\\tau\}, contain pure estimation noise\. We forecast𝐱~t\+q\\mathbf\{\\tilde\{x\}\}\_\{t\+q\}from a window ofT=6T=6steps and read out only its last position, so that every cross\-time relation must pass through a Kronecker term\. We train a single\-layer KVNN\-S \(F=16F=16,K=1K=1, blocks normalized to unit spectral norm\) onM=1\.2⋅104M=1\.2\\cdot 10^\{4\}windows for6060epochs with Adam \(η=10−2\\eta=10^\{\-2\}, batch size128128\)\.

Figure[5](https://arxiv.org/html/2609.25326#S5.F5)\(left\) reports, for each true lagqq, the share of the retained coefficient mass carried by each lag atλg=2\\lambda\_\{g\}=2, averaged over three seeds\. The penalty concentrates0\.970\.97,1\.001\.00and1\.001\.00of the mass on the true lag forq=3,4,5q=3,4,5, respectively, leaving at most0\.030\.03on any wrong lag\. Lags 1\-2 are grouped in the same column as they are never the ground truth\.

### 5\.4IMPACT OF HYPERPARAMETERS \(O4\)

Figure[5](https://arxiv.org/html/2609.25326#S5.F5)\(middle and right\) shows the effect ofTT,RRandλg\\lambda\_\{g\}on LargeCap and AQWan\. IncreasingTTimproves all KVNN variants on both datasets, as the model accesses longer temporal dependencies\. A largerTTdoes not require a largerRRin KVNN\-LR\. Often,R=1R=1is insufficient, butR=3R=3already matchesR=8R=8and KVNN\-S at allTT, so a handful of Kronecker terms suffices to represent windows of up to2424steps\. The effect ofλg\\lambda\_\{g\}is model\- and dataset\-dependent\. KVNN\-S benefits from stronger sparsification on LargeCap, KVNN\-LR withR∈\{3,8\}R\\in\\\{3,8\\\}is largely insensitive to it, and KVNN\-LR withR=1R=1slightly degrades on AQWan, as pruning any of its few terms removes useful information\. Overall,λg\\lambda\_\{g\}is not a critical hyperparameter for accuracy, and its main benefit is the reduced number of active Kronecker terms\.

## 6CONCLUSION

We introduced the spatiotemporal Kronecker coVariance Neural Network \(KVNN\), an architecture designed to efficiently and flexibly process the joint spatiotemporal correlations inherent in multivariate time series via a sum\-of\-Kronecker representation\. We provided a spectral analysis and proved that, under stationary assumptions, the KVF maintains strong stability bounds against finite\-sample estimation errors\. Our empirical evaluations on real\-world forecasting tasks validated the forecasting capability of KVNNs on five datasets\.

## APPENDIX

We first provide a useful Lemma, then prove Theorem[1](https://arxiv.org/html/2609.25326#Thmtheorem1)\.

###### Lemma 1

Let𝐌,𝐌^∈ℂN×N\\mathbf\{M\},\\mathbf\{\\hat\{M\}\}\\in\\mathbb\{C\}^\{N\\times N\}be Hermitian with eigenvalues\{αi\}\\\{\\alpha\_\{i\}\\\}and\{βj\}\\\{\\beta\_\{j\}\\\}, and lethhbe a real polynomial with\|h⁡\(αi\)−h⁡\(βj\)\|≤P​\|αi−βj\|\|h\(\\alpha\_\{i\}\)\-h\(\\beta\_\{j\}\)\|\\leq P\|\\alpha\_\{i\}\-\\beta\_\{j\}\|for alli,ji,j\. Then

‖h⁡\(𝐌\)−h⁡\(𝐌^\)‖2\\displaystyle\\\|h\(\\mathbf\{M\}\)\-h\(\\mathbf\{\\hat\{M\}\}\)\\\|\_\{2\}≤P​N​‖𝐌−𝐌^‖2\.\\displaystyle\\leq P\\sqrt\{N\}\\,\\\|\\mathbf\{M\}\-\\mathbf\{\\hat\{M\}\}\\\|\_\{2\}\.\(13\)

Proof :Write𝐌=∑iαi​𝐩i​𝐩i𝖧\\mathbf\{M\}=\\sum\_\{i\}\\alpha\_\{i\}\\mathbf\{p\}\_\{i\}\\mathbf\{p\}\_\{i\}^\{\\mathsf\{H\}\}and𝐌^=∑jβj​𝐪j​𝐪j𝖧\\mathbf\{\\hat\{M\}\}=\\sum\_\{j\}\\beta\_\{j\}\\mathbf\{q\}\_\{j\}\\mathbf\{q\}\_\{j\}^\{\\mathsf\{H\}\}with orthonormal bases\{𝐩i\}\\\{\\mathbf\{p\}\_\{i\}\\\},\{𝐪j\}\\\{\\mathbf\{q\}\_\{j\}\\\}\. Since𝐩i𝖧​h​\(𝐌\)\\mathbf\{p\}\_\{i\}^\{\\mathsf\{H\}\}h\(\\mathbf\{M\}\)andh⁡\(𝐌^\)​𝐪jh\(\\mathbf\{\\hat\{M\}\}\)\\mathbf\{q\}\_\{j\}act through the respective eigenvalues, we get

𝐩i𝖧​\(h⁡\(𝐌\)−h⁡\(𝐌^\)\)​𝐪j\\displaystyle\\mathbf\{p\}\_\{i\}^\{\\mathsf\{H\}\}\\big\(h\(\\mathbf\{M\}\)\-h\(\\mathbf\{\\hat\{M\}\}\)\\big\)\\mathbf\{q\}\_\{j\}=\(h⁡\(αi\)−h⁡\(βj\)\)​𝐩i𝖧​𝐪j\.\\displaystyle=\\big\(h\(\\alpha\_\{i\}\)\-h\(\\beta\_\{j\}\)\\big\)\\,\\mathbf\{p\}\_\{i\}^\{\\mathsf\{H\}\}\\mathbf\{q\}\_\{j\}\.\(14\)By settingh⁡\(x\)=xh\(x\)=x:

𝐩i𝖧​\(𝐌−𝐌^\)​𝐪j\\displaystyle\\mathbf\{p\}\_\{i\}^\{\\mathsf\{H\}\}\\big\(\\mathbf\{M\}\-\\mathbf\{\\hat\{M\}\}\\big\)\\mathbf\{q\}\_\{j\}=\(αi−βj\)​𝐩i𝖧​𝐪j\.\\displaystyle=\\big\(\\alpha\_\{i\}\-\\beta\_\{j\}\\big\)\\,\\mathbf\{p\}\_\{i\}^\{\\mathsf\{H\}\}\\mathbf\{q\}\_\{j\}\.\(15\)
Hence\|𝐩i𝖧​\(h⁡\(𝐌\)−h⁡\(𝐌^\)\)​𝐪j\|≤P​\|𝐩i𝖧​\(𝐌−𝐌^\)​𝐪j\|\|\\mathbf\{p\}\_\{i\}^\{\\mathsf\{H\}\}\(h\(\\mathbf\{M\}\)\-h\(\\mathbf\{\\hat\{M\}\}\)\)\\mathbf\{q\}\_\{j\}\|\\leq P\\,\|\\mathbf\{p\}\_\{i\}^\{\\mathsf\{H\}\}\(\\mathbf\{M\}\-\\mathbf\{\\hat\{M\}\}\)\\mathbf\{q\}\_\{j\}\|for alli,ji,j\(both sides vanish ifαi=βj\\alpha\_\{i\}=\\beta\_\{j\}\)\. For any𝐗\\mathbf\{X\}and unitary𝐏=\[𝐩1,…,𝐩N\]\\mathbf\{P\}=\[\\mathbf\{p\}\_\{1\},\\dots,\\mathbf\{p\}\_\{N\}\],𝐐=\[𝐪1,…,𝐪N\]\\mathbf\{Q\}=\[\\mathbf\{q\}\_\{1\},\\dots,\\mathbf\{q\}\_\{N\}\]we have‖𝐗‖F2=‖𝐏𝖧​𝐗𝐐‖F2=∑i,j\|𝐩i𝖧​𝐗𝐪j\|2\\\|\\mathbf\{X\}\\\|\_\{F\}^\{2\}=\\\|\\mathbf\{P\}^\{\\mathsf\{H\}\}\\mathbf\{X\}\\mathbf\{Q\}\\\|\_\{F\}^\{2\}=\\sum\_\{i,j\}\|\\mathbf\{p\}\_\{i\}^\{\\mathsf\{H\}\}\\mathbf\{X\}\\mathbf\{q\}\_\{j\}\|^\{2\}; summing the squared inequalities overi,ji,jgives‖h⁡\(𝐌\)−h⁡\(𝐌^\)‖F≤P​‖𝐌−𝐌^‖F\\\|h\(\\mathbf\{M\}\)\-h\(\\mathbf\{\\hat\{M\}\}\)\\\|\_\{F\}\\leq P\\\|\\mathbf\{M\}\-\\mathbf\{\\hat\{M\}\}\\\|\_\{F\}\. Finally, the norm inequalities∥⋅∥2≤∥⋅∥F≤N∥⋅∥2\\\|\\cdot\\\|\_\{2\}\\leq\\\|\\cdot\\\|\_\{F\}\\leq\\sqrt\{N\}\\\|\\cdot\\\|\_\{2\}complete the result\.■\\blacksquare

### PROOF OF THEOREM[1](https://arxiv.org/html/2609.25326#Thmtheorem1)

Let𝐇r​\(𝐌\)=∑k=0Khk​r​𝐌k\\mathbf\{H\}\_\{r\}\(\\mathbf\{M\}\)=\\sum\_\{k=0\}^\{K\}h\_\{kr\}\\mathbf\{M\}^\{k\}and𝐄τ=𝐂^τ−𝐂τ\\mathbf\{E\}\_\{\\tau\}=\\mathbf\{\\hat\{C\}\}\_\{\\tau\}\-\\mathbf\{C\}\_\{\\tau\}forτ=0,…,T−1\\tau=0,\\dots,T\-1\. By triangle inequality, we get

∥𝐮0\\displaystyle\\\|\\mathbf\{u\}\_\{0\}−𝐮∥2=∥∑r\[𝐇r\(𝐀r⊗𝐁r\)−𝐇r\(𝐀r⊗𝐁^r\)\]𝐱~∥2\\displaystyle\-\\mathbf\{u\}\\\|\_\{2\}=\\Big\\\|\\sum\_\{r\}\\big\[\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}\)\-\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{\\hat\{B\}\}\_\{r\}\)\\big\]\\mathbf\{\\tilde\{x\}\}\\Big\\\|\_\{2\}≤∑r‖𝐇r​\(𝐀r⊗𝐁r\)−𝐇r​\(𝐀r⊗𝐁^r\)‖2​‖𝐱~‖2\.\\displaystyle\\leq\\sum\_\{r\}\\big\\\|\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}\)\-\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{\\hat\{B\}\}\_\{r\}\)\\big\\\|\_\{2\}\\,\\\|\\mathbf\{\\tilde\{x\}\}\\\|\_\{2\}\.\(16\)Then, for a fixedrr\(we omit the subscriptrrfrom eigenvectors for simplicity\), using the mixed\-product property yields

𝐀r⊗𝐁r=\(𝐔⊗𝐈N\)​\(𝚲⊗𝐁r\)​\(𝐔⊗𝐈N\)𝖧,\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}=\(\\mathbf\{U\}\\otimes\\mathbf\{I\}\_\{N\}\)\(\\mathbf\{\\Lambda\}\\otimes\\mathbf\{B\}\_\{r\}\)\(\\mathbf\{U\}\\otimes\\mathbf\{I\}\_\{N\}\)^\{\\mathsf\{H\}\},\(17\)and since𝐔⊗𝐈N\\mathbf\{U\}\\otimes\\mathbf\{I\}\_\{N\}is unitary,\(𝐀r⊗𝐁r\)k=\(𝐔⊗𝐈N\)​\(𝚲⊗𝐁r\)k​\(𝐔⊗𝐈N\)𝖧\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}\)^\{k\}=\(\\mathbf\{U\}\\otimes\\mathbf\{I\}\_\{N\}\)\(\\mathbf\{\\Lambda\}\\otimes\\mathbf\{B\}\_\{r\}\)^\{k\}\(\\mathbf\{U\}\\otimes\\mathbf\{I\}\_\{N\}\)^\{\\mathsf\{H\}\}for allkk, so that𝐇r​\(𝐀r⊗𝐁r\)=\(𝐔⊗𝐈N\)​𝐇r​\(𝚲⊗𝐁r\)​\(𝐔⊗𝐈N\)𝖧\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}\)=\(\\mathbf\{U\}\\otimes\\mathbf\{I\}\_\{N\}\)\\mathbf\{H\}\_\{r\}\(\\mathbf\{\\Lambda\}\\otimes\\mathbf\{B\}\_\{r\}\)\(\\mathbf\{U\}\\otimes\\mathbf\{I\}\_\{N\}\)^\{\\mathsf\{H\}\}, and likewise for𝐁^r\\mathbf\{\\hat\{B\}\}\_\{r\}\. By unitary invariance of the spectral norm,

‖𝐇r​\(𝐀r⊗𝐁r\)−𝐇r​\(𝐀r⊗𝐁^r\)‖2=\\displaystyle\\big\\\|\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}\)\-\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{\\hat\{B\}\}\_\{r\}\)\\big\\\|\_\{2\}=‖𝐇r​\(𝚲⊗𝐁r\)−𝐇r​\(𝚲⊗𝐁^r\)‖2=\\displaystyle\\big\\\|\\mathbf\{H\}\_\{r\}\(\\mathbf\{\\Lambda\}\\otimes\\mathbf\{B\}\_\{r\}\)\-\\mathbf\{H\}\_\{r\}\(\\mathbf\{\\Lambda\}\\otimes\\mathbf\{\\hat\{B\}\}\_\{r\}\)\\big\\\|\_\{2\}=maxi⁡‖𝐇r​\(λi​𝐁r\)−𝐇r​\(λi​𝐁^r\)‖2,\\displaystyle\\max\_\{i\}\\big\\\|\\mathbf\{H\}\_\{r\}\(\\lambda\_\{i\}\\mathbf\{B\}\_\{r\}\)\-\\mathbf\{H\}\_\{r\}\(\\lambda\_\{i\}\\mathbf\{\\hat\{B\}\}\_\{r\}\)\\big\\\|\_\{2\},\(18\)where the last step uses the fact that the matrix𝚲⊗𝐁r\\mathbf\{\\Lambda\}\\otimes\\mathbf\{B\}\_\{r\}is block\-diagonal withN×NN\\times Nblocksλi​𝐁r\\lambda\_\{i\}\\mathbf\{B\}\_\{r\}, a polynomial acts block\-wise, and the spectral norm of a block\-diagonal matrix is the maximum over its blocks\.

If𝐀r,𝐁r\\mathbf\{A\}\_\{r\},\\mathbf\{B\}\_\{r\}are symmetric,λi\\lambda\_\{i\}is real andλi​𝐁r\\lambda\_\{i\}\\mathbf\{B\}\_\{r\},λi​𝐁^r\\lambda\_\{i\}\\mathbf\{\\hat\{B\}\}\_\{r\}are real symmetric\. If𝐀r,𝐁r\\mathbf\{A\}\_\{r\},\\mathbf\{B\}\_\{r\}are skew\-symmetric,λi=i​θi\\lambda\_\{i\}=i\\theta\_\{i\}withθi∈ℝ\\theta\_\{i\}\\in\\mathbb\{R\}, and\(λi​𝐁r\)𝖧=−i​θi​𝐁r𝖳=i​θi​𝐁r=λi​𝐁r\(\\lambda\_\{i\}\\mathbf\{B\}\_\{r\}\)^\{\\mathsf\{H\}\}=\-i\\theta\_\{i\}\\mathbf\{B\}\_\{r\}^\{\\mathsf\{T\}\}=i\\theta\_\{i\}\\mathbf\{B\}\_\{r\}=\\lambda\_\{i\}\\mathbf\{B\}\_\{r\}, so the blocks are Hermitian; the same holds forλi​𝐁^r\\lambda\_\{i\}\\mathbf\{\\hat\{B\}\}\_\{r\}\. In both cases the eigenvalues ofλi​𝐁r\\lambda\_\{i\}\\mathbf\{B\}\_\{r\}andλi​𝐁^r\\lambda\_\{i\}\\mathbf\{\\hat\{B\}\}\_\{r\}are of the formλi​w\\lambda\_\{i\}wwithwwan eigenvalue of𝐁r\\mathbf\{B\}\_\{r\}or𝐁^r\\mathbf\{\\hat\{B\}\}\_\{r\}, i\.e\., they are eigenvalues of𝐀r⊗𝐁r\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}or𝐀r⊗𝐁^r\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{\\hat\{B\}\}\_\{r\}\. Assumption[1](https://arxiv.org/html/2609.25326#Thmassumption1)and Lemma[1](https://arxiv.org/html/2609.25326#Thmlemma1)give

∥𝐇r\(λi𝐁r\)\\displaystyle\\big\\\|\\mathbf\{H\}\_\{r\}\(\\lambda\_\{i\}\\mathbf\{B\}\_\{r\}\)−𝐇r\(λi𝐁^r\)∥2≤PN∥λi\(𝐁r−𝐁^r\)∥2\\displaystyle\-\\mathbf\{H\}\_\{r\}\(\\lambda\_\{i\}\\mathbf\{\\hat\{B\}\}\_\{r\}\)\\big\\\|\_\{2\}\\leq P\\sqrt\{N\}\\,\\\|\\lambda\_\{i\}\(\\mathbf\{B\}\_\{r\}\-\\mathbf\{\\hat\{B\}\}\_\{r\}\)\\\|\_\{2\}=P​N​\|λi\|​‖𝐁r−𝐁^r‖2\.\\displaystyle=P\\sqrt\{N\}\\,\|\\lambda\_\{i\}\|\\,\\\|\\mathbf\{B\}\_\{r\}\-\\mathbf\{\\hat\{B\}\}\_\{r\}\\\|\_\{2\}\.\(19\)Taking the maximum overiiin \([18](https://arxiv.org/html/2609.25326#Sx1.E18)\) yields

‖𝐇r​\(𝐀r⊗𝐁r\)−𝐇r​\(𝐀r⊗𝐁^r\)‖2≤P​N​‖𝐀r‖2​‖𝐁r−𝐁^r‖2\.\\big\\\|\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{B\}\_\{r\}\)\-\\mathbf\{H\}\_\{r\}\(\\mathbf\{A\}\_\{r\}\\otimes\\mathbf\{\\hat\{B\}\}\_\{r\}\)\\big\\\|\_\{2\}\\leq P\\sqrt\{N\}\\,\\\|\\mathbf\{A\}\_\{r\}\\\|\_\{2\}\\,\\\|\\mathbf\{B\}\_\{r\}\-\\mathbf\{\\hat\{B\}\}\_\{r\}\\\|\_\{2\}\.\(20\)
Forτ=0\\tau=0,‖𝐈T‖2=1\\\|\\mathbf\{I\}\_\{T\}\\\|\_\{2\}=1and𝐁r−𝐁^r=−𝐄0\\mathbf\{B\}\_\{r\}\-\\mathbf\{\\hat\{B\}\}\_\{r\}=\-\\mathbf\{E\}\_\{0\}\. Forτ≥1\\tau\\geq 1, every row and every column of𝐃τ\\mathbf\{D\}\_\{\\tau\}contains at most one non\-zero entry, equal to one, so𝐃τ\\mathbf\{D\}\_\{\\tau\}is a partial permutation matrix and‖𝐃τ‖2=‖𝐃τ⊤‖2=1\\\|\\mathbf\{D\}\_\{\\tau\}\\\|\_\{2\}=\\\|\\mathbf\{D\}\_\{\\tau\}^\{\\top\}\\\|\_\{2\}=1\. By the triangle inequality,‖𝐃τ±𝐃τ⊤‖2≤2\\\|\\mathbf\{D\}\_\{\\tau\}\\pm\\mathbf\{D\}\_\{\\tau\}^\{\\top\}\\\|\_\{2\}\\leq 2\. Moreover,

‖𝐂τs−𝐂^τs‖2=12​‖𝐄τ\+𝐄τ𝖳‖2≤‖𝐄τ‖2,\\displaystyle\\\|\\mathbf\{C\}\_\{\\tau\}^\{s\}\-\\mathbf\{\\hat\{C\}\}\_\{\\tau\}^\{s\}\\\|\_\{2\}=\\tfrac\{1\}\{2\}\\\|\\mathbf\{E\}\_\{\\tau\}\+\\mathbf\{E\}\_\{\\tau\}^\{\\mathsf\{T\}\}\\\|\_\{2\}\\leq\\\|\\mathbf\{E\}\_\{\\tau\}\\\|\_\{2\},\(21\)‖𝐂τa−𝐂^τa‖2=12​‖𝐄τ−𝐄τ𝖳‖2≤‖𝐄τ‖2\.\\displaystyle\\\|\\mathbf\{C\}\_\{\\tau\}^\{a\}\-\\mathbf\{\\hat\{C\}\}\_\{\\tau\}^\{a\}\\\|\_\{2\}=\\tfrac\{1\}\{2\}\\\|\\mathbf\{E\}\_\{\\tau\}\-\\mathbf\{E\}\_\{\\tau\}^\{\\mathsf\{T\}\}\\\|\_\{2\}\\leq\\\|\\mathbf\{E\}\_\{\\tau\}\\\|\_\{2\}\.\(22\)Substituting in \([20](https://arxiv.org/html/2609.25326#Sx1.E20)\) and summing over the2​T−12T\-1terms in \([16](https://arxiv.org/html/2609.25326#Sx1.E16)\),

‖𝐮0−𝐮‖2\\displaystyle\\\|\\mathbf\{u\}\_\{0\}\-\\mathbf\{u\}\\\|\_\{2\}≤P​N​\(‖𝐄0‖2\+2​∑τ=1T−1\(‖𝐄τ‖2\+‖𝐄τ‖2\)\)​‖𝐱~‖2\\displaystyle\\leq P\\sqrt\{N\}\\Big\(\\\|\\mathbf\{E\}\_\{0\}\\\|\_\{2\}\+2\\sum\_\{\\tau=1\}^\{T\-1\}\\big\(\\\|\\mathbf\{E\}\_\{\\tau\}\\\|\_\{2\}\+\\\|\\mathbf\{E\}\_\{\\tau\}\\\|\_\{2\}\\big\)\\Big\)\\\|\\mathbf\{\\tilde\{x\}\}\\\|\_\{2\}≤4​P​N​\(∑τ=0T−1‖𝐄τ‖2\)​‖𝐱~‖2\\displaystyle\\leq 4P\\sqrt\{N\}\\Big\(\\sum\_\{\\tau=0\}^\{T\-1\}\\\|\\mathbf\{E\}\_\{\\tau\}\\\|\_\{2\}\\Big\)\\\|\\mathbf\{\\tilde\{x\}\}\\\|\_\{2\}≤4​P​T​N​maxτ​‖𝐄τ‖2​‖𝐱~‖2\.\\displaystyle\\leq 4PT\\sqrt\{N\}\\,\\max\_\{\\tau\}\\\|\\mathbf\{E\}\_\{\\tau\}\\\|\_\{2\}\\,\\\|\\mathbf\{\\tilde\{x\}\}\\\|\_\{2\}\.\(23\)This deterministic bound holds for any realization of the estimation errors\.

Since𝐂^τ\\mathbf\{\\hat\{C\}\}\_\{\\tau\}is the average of theT−τT\-\\taublocks on theτ\\tau\-th block diagonal of𝐂~s\\mathbf\{\\tilde\{C\}\}\_\{s\}, the error𝐄τ\\mathbf\{E\}\_\{\\tau\}is a convex combination of sub\-matrices of𝐂~s−𝐂~0\\mathbf\{\\tilde\{C\}\}\_\{s\}\-\\mathbf\{\\tilde\{C\}\}\_\{0\}, and hencemaxτ⁡‖𝐄τ‖2≤‖𝐂~s−𝐂~0‖2\\max\_\{\\tau\}\\\|\\mathbf\{E\}\_\{\\tau\}\\\|\_\{2\}\\leq\\\|\\mathbf\{\\tilde\{C\}\}\_\{s\}\-\\mathbf\{\\tilde\{C\}\}\_\{0\}\\\|\_\{2\}\. The windows𝐱~t\\mathbf\{\\tilde\{x\}\}\_\{t\}areMMindependent zero\-mean sub\-Gaussian vectors inℝN​T\\mathbb\{R\}^\{NT\}with covariance𝐂~0\\mathbf\{\\tilde\{C\}\}\_\{0\}, so by\[[34](https://arxiv.org/html/2609.25326#bib.bib1), Theorem 4\.7\.1\]there is a universal constantC\>0C\>0such that for everyu≥0u\\geq 0, with probability at least1−2​e−u1\-2e^\{\-u\},

‖𝐂~s−𝐂~0‖2\\displaystyle\\\|\\mathbf\{\\tilde\{C\}\}\_\{s\}\-\\mathbf\{\\tilde\{C\}\}\_\{0\}\\\|\_\{2\}≤C​κ2​‖𝐂~0‖2​\(N​T\+uM\+N​T\+uM\)\\displaystyle\\leq C\\kappa^\{2\}\\\|\\mathbf\{\\tilde\{C\}\}\_\{0\}\\\|\_\{2\}\\left\(\\sqrt\{\\frac\{NT\+u\}\{M\}\}\+\\frac\{NT\+u\}\{M\}\\right\)≤c​N​T\+uM,\\displaystyle\\leq c\\sqrt\{\\frac\{NT\+u\}\{M\}\},\(24\)whereκ\\kappais the sub\-Gaussian norm of𝐱~\\mathbf\{\\tilde\{x\}\},c=2​C​κ2​‖𝐂~0‖2c=2C\\kappa^\{2\}\\\|\\mathbf\{\\tilde\{C\}\}\_\{0\}\\\|\_\{2\}, and the last step usesM≥N​T\+uM\\geq NT\+u\. Inserting this into \([23](https://arxiv.org/html/2609.25326#Sx1.E23)\) withu=log⁡\(2/δ\)u=\\log\(2/\\delta\)and absorbingδ\\deltaand constants gives \([12](https://arxiv.org/html/2609.25326#S4.E12)\)\.■\\blacksquare

## References

- \[1\]A\. Alanqary, A\. Alomar, and D\. Shah\(2021\)Change point detection via multivariate singular spectrum analysis\.InAdvances in Neural Information Processing Systems,M\. Ranzato, A\. Beygelzimer, Y\. Dauphin, P\.S\. Liang, and J\. W\. Vaughan \(Eds\.\),Vol\.34,pp\. 23218–23230\.External Links:[Link](https://proceedings.neurips.cc/paper_files/paper/2021/file/c348616cd8a86ee661c7c98800678fad-Paper.pdf)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[2\]S\. Basu and G\. Michailidis\(2015\)Regularized estimation in sparse high\-dimensional time series models\.The Annals of Statistics,pp\. 1535–1567\.Cited by:[§4\.2](https://arxiv.org/html/2609.25326#S4.SS2.p2.1)\.
- \[3\]A\. Belouchrani, K\. Abed\-Meraim, J\. Cardoso, and E\. Moulines\(1997\)A blind source separation technique using second\-order statistics\.IEEE Transactions on signal processing45\(2\),pp\. 434–444\.Cited by:[§4\.1](https://arxiv.org/html/2609.25326#S4.SS1.p2.1)\.
- \[4\]A\. Cavallo, Z\. Gao, and E\. Isufi\(2026\)Sparse covariance neural networks\.IEEE Transactions on Signal Processing74\(\),pp\. 3077–3091\.External Links:[Document](https://dx.doi.org/10.1109/TSP.2026.3714283)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p6.1)\.
- \[5\]A\. Cavallo, A\. Raghuvanshi, S\. P\. Chepuri, and E\. Isufi\(2026\)Covariance scattering transforms\.InProceedings of the AAAI Conference on Artificial Intelligence,Vol\.40,pp\. 19925–19933\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p6.1)\.
- \[6\]A\. Cavallo, M\. Sabbaqi, and E\. Isufi\(2024\)Spatiotemporal covariance neural networks\.InJoint European Conf\. on Mach\. Learn\. and Knowledge Discovery in Databases,pp\. 18–34\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p2.1),[§1](https://arxiv.org/html/2609.25326#S1.p6.1),[§4\.2](https://arxiv.org/html/2609.25326#S4.SS2.p2.1),[§5\.1](https://arxiv.org/html/2609.25326#S5.SS1.p3.1)\.
- \[7\]Z\. Che, S\. Purushotham, K\. Cho, D\. Sontag, and Y\. Liu\(2018\)Recurrent neural networks for multivariate time series with missing values\.Scientific reports8\(1\),pp\. 6085\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.
- \[8\]F\. Gama, J\. Bruna, and A\. Ribeiro\(2020\)Stability properties of graph neural networks\.IEEE Transactions on Signal Processing68,pp\. 5680–5695\.Cited by:[§4\.2](https://arxiv.org/html/2609.25326#S4.SS2.p2.1)\.
- \[9\]M\. B\. Garman and M\. J\. Klass\(1980\)On the estimation of security price volatilities from historical data\.Journal of business,pp\. 67–78\.Cited by:[§5\.1](https://arxiv.org/html/2609.25326#S5.SS1.p1.1)\.
- \[10\]A\. Georgoutsos\(2025\)Lagged spatiotemporal covariance neural networks\.Master’s thesis,Delft University of Technology, Faculty of Electrical Engineering, Mathematics and Computer Science\.External Links:[Link](https://resolver.tudelft.nl/uuid:4e9569d0-0361-4fb4-97a8-16f8bef882eb)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p6.1),[§3](https://arxiv.org/html/2609.25326#S3.p1.2),[§5\.1](https://arxiv.org/html/2609.25326#S5.SS1.p3.1)\.
- \[11\]P\. Hewage, A\. Behera, M\. Trovati, E\. Pereira, M\. Ghahremani, F\. Palmieri, and Y\. Liu\(2020\)Temporal convolutional neural \(tcn\) network for an effective weather forecasting using time\-series data from the local weather station\.Soft Comput\.24\(21\),pp\. 16453–16482\.External Links:ISSN 1432\-7643,[Link](https://doi.org/10.1007/s00500-020-04954-0),[Document](https://dx.doi.org/10.1007/s00500-020-04954-0)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[12\]S\. Hochreiter and J\. Schmidhuber\(1997\)Long short\-term memory\.Neural computation9\(8\),pp\. 1735–1780\.Cited by:[§5\.1](https://arxiv.org/html/2609.25326#S5.SS1.p3.1)\.
- \[13\]A\. G\. Ismaeel, K\. Janardhanan, M\. Sankar, Y\. Natarajan, S\. N\. Mahmood, S\. Alani, and A\. H\. Shather\(2023\)Traffic pattern classification in smart cities using deep recurrent neural network\.Sustainability15\(19\),pp\. 14522\.External Links:ISSN 2071\-1050,[Link](http://dx.doi.org/10.3390/su151914522),[Document](https://dx.doi.org/10.3390/su151914522)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[14\]L\. Ji, Y\. Zou, K\. He, and B\. Zhu\(2019\)Carbon futures price forecasting based with arima\-cnn\-lstm model\.Procedia Comput\. Sci\.162\(C\),pp\. 33–38\.External Links:ISSN 1877\-0509,[Link](https://doi.org/10.1016/j.procs.2019.11.254),[Document](https://dx.doi.org/10.1016/j.procs.2019.11.254)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[15\]M\. Jin, H\. Y\. Koh, Q\. Wen, D\. Zambon, C\. Alippi, G\. I\. Webb, I\. King, and S\. Pan\(2024\)A survey on graph neural networks for time series: forecasting, classification, imputation, and anomaly detection\.IEEE transactions on pattern analysis and machine intelligence46\(12\),pp\. 10466–10485\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.
- \[16\]T\.S\. Jombart, S\. Devillard, A\. Dufour, and D\. Pontier\(2008\)Revealing cryptic spatial patterns in genetic variability by a new multivariate method\.Heredity101,pp\. 92–103\.External Links:[Document](https://dx.doi.org/10.1038/hdy.2008.34)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[17\]Z\. Karevan and J\. A\.K\. Suykens\(2020\)Transductive lstm for time\-series prediction: an application to weather forecasting\.Neural Networks125,pp\. 1–9\.External Links:ISSN 0893\-6080,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.neunet.2019.12.030),[Link](https://www.sciencedirect.com/science/article/pii/S0893608020300010)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[18\]M\. Krzyśko, P\. Nijkamp, W\. Ratajczak, W\. Wołyński, and B\. Wenerska\(2024\)Spatio\-temporal principal component analysis\.Spatial Economic Analysis19\(1\),pp\. 8–29\.External Links:[Document](https://dx.doi.org/10.1080/17421772.2023.2237532),[Link](https://doi.org/10.1080/17421772.2023.2237532),https://doi\.org/10\.1080/17421772\.2023\.2237532Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[19\]G\. Lai, W\. Chang, Y\. Yang, and H\. Liu\(2018\)Modeling long\-and short\-term temporal patterns with deep neural networks\.InThe 41st international ACM SIGIR conference on research & development in information retrieval,pp\. 95–104\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.
- \[20\]Y\. Li, R\. Yu, C\. Shahabi, and Y\. Liu\(2017\)Diffusion convolutional recurrent neural network: data\-driven traffic forecasting\.arXiv preprint arXiv:1707\.01926\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.
- \[21\]Z\. Lin\(2018\)Modelling and forecasting the stock market volatility of sse composite index using garch models\.Future Gener\. Comput\. Syst\.79\(P3\),pp\. 960–972\.External Links:ISSN 0167\-739X,[Link](https://doi.org/10.1016/j.future.2017.08.033),[Document](https://dx.doi.org/10.1016/j.future.2017.08.033)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[22\]X\. Mao and P\. Shang\(2019\)Multivariate singular spectrum analysis for traffic time series\.Physica A: Statistical Mechanics and its Applications526,pp\. 121063\.External Links:ISSN 0378\-4371,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.physa.2019.121063),[Link](https://www.sciencedirect.com/science/article/pii/S037843711930648X)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[23\]L\. Molgedey and H\. G\. Schuster\(1994\)Separation of a mixture of independent signals using time delayed correlations\.Physical review letters72\(23\),pp\. 3634\.Cited by:[§4\.1](https://arxiv.org/html/2609.25326#S4.SS1.p2.1)\.
- \[24\]X\. Qiu, J\. Hu, L\. Zhou, X\. Wu, J\. Du, B\. Zhang, C\. Guo, A\. Zhou, C\. S\. Jensen, Z\. Sheng,et al\.\(2024\)Tfb: towards comprehensive and fair benchmarking of time series forecasting methods\.arXiv preprint arXiv:2403\.20150\.Cited by:[§5\.1](https://arxiv.org/html/2609.25326#S5.SS1.p1.1)\.
- \[25\]O\. Roy, Y\. Moshfeghi, and K\. M\. Smith\(2026\)Covariance density neural networks\.Transactions on Machine Learning Research\.Note:External Links:ISSN 2835\-8856,[Link](https://openreview.net/forum?id=TwCkGi5XFB)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p6.1)\.
- \[26\]M\. Sabbaqi and E\. Isufi\(2023\)Graph\-time convolutional neural networks: architecture and theoretical analysis\.IEEE Transactions on Pattern Analysis and Machine Intelligence45\(12\),pp\. 14625–14638\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.
- \[27\]A\. Sandryhaila and J\. M\. Moura\(2014\)Discrete signal processing on graphs: frequency analysis\.IEEE Transactions on signal processing62\(12\),pp\. 3042–3054\.Cited by:[§4\.1](https://arxiv.org/html/2609.25326#S4.SS1.p2.1)\.
- \[28\]S\. Sihag, A\. Cavallo, E\. Isufi, G\. Mateos, and A\. Ribeiro\(2026\)Learning with covariance matrices: principal component analysis meets learning with graphs\.arXiv preprint arXiv:2609\.10490\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p2.1),[§1](https://arxiv.org/html/2609.25326#S1.p6.1),[§3](https://arxiv.org/html/2609.25326#S3.p1.1)\.
- \[29\]S\. Sihag, G\. Mateos, C\. McMillan, and A\. Ribeiro\(2022\)coVariance neural networks\.Advances in Neural Info\. Process\. Syst\.35,pp\. 17003–17016\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p2.1),[§1](https://arxiv.org/html/2609.25326#S1.p6.1),[§3](https://arxiv.org/html/2609.25326#S3.p1.1),[§4\.2](https://arxiv.org/html/2609.25326#S4.SS2.p2.1),[§5\.1](https://arxiv.org/html/2609.25326#S5.SS1.p3.1)\.
- \[30\]S\. Sihag, G\. Mateos, C\. McMillan, and A\. Ribeiro\(2024\)Explainable brain age prediction using covariance neural networks\.Advances in Neural Info\. Process\. Syst\.36\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p6.1)\.
- \[31\]S\. Sihag, G\. Mateos, C\. McMillan, and A\. Ribeiro\(2024\)Transferability of covariance neural networks\.IEEE J\. Sel\. Topics Signal Process\.18\(2\),pp\. 199–215\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p6.1)\.
- \[32\]S\. Sihag, G\. Mateos, and A\. Ribeiro\(2025\)Disentangling neurodegeneration with brain age gap prediction models: a graph signal processing perspective\.IEEE Signal Processing Magazine42\(4\),pp\. 58–77\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p6.1)\.
- \[33\]T\. Tsiligkaridis and A\. O\. Hero\(2013\)Covariance estimation in high dimensions via kronecker product expansions\.IEEE Transactions on Signal Processing61\(21\),pp\. 5347–5360\.Cited by:[§2\.2](https://arxiv.org/html/2609.25326#S2.SS2.p1.1)\.
- \[34\]R\. Vershynin\(2018\)High\-Dimensional Probability: An Introduction with Applications in Data Science\.Cambridge Series in Statistical and Probabilistic Mathematics,Cambridge University Press\.Cited by:[PROOF OF THEOREM](https://arxiv.org/html/2609.25326#Sx1.SSx1.p4.1)\.
- \[35\]M\. D\. Weinberg and M\. S\. Petersen\(2020\)Using multichannel singular spectrum analysis to study galaxy dynamics\.Monthly Notices of the Royal Astronomical Society501\(4\),pp\. 5408–5423\.External Links:ISSN 1365\-2966,[Link](http://dx.doi.org/10.1093/mnras/staa3997),[Document](https://dx.doi.org/10.1093/mnras/staa3997)Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p1.1)\.
- \[36\]Z\. Wu, S\. Pan, F\. Chen, G\. Long, C\. Zhang, and P\. S\. Yu\(2020\)A comprehensive survey on graph neural networks\.IEEE transactions on neural networks and learning systems32\(1\),pp\. 4–24\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.
- \[37\]Z\. Wu, S\. Pan, G\. Long, J\. Jiang, and C\. Zhang\(2019\)Graph wavenet for deep spatial\-temporal graph modeling\.arXiv preprint arXiv:1906\.00121\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.
- \[38\]G\. Zerveas, S\. Jayaraman, D\. Patel, A\. Bhamidipaty, and C\. Eickhoff\(2021\)A transformer\-based framework for multivariate time series representation learning\.InProceedings of the 27th ACM SIGKDD conference on knowledge discovery & data mining,pp\. 2114–2124\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.
- \[39\]Y\. Zhang and J\. Yan\(2023\)Crossformer: transformer utilizing cross\-dimension dependency for multivariate time series forecasting\.InThe eleventh international conference on learning representations,Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.
- \[40\]H\. Zhou, S\. Zhang, J\. Peng, S\. Zhang, J\. Li, H\. Xiong, and W\. Zhang\(2021\)Informer: beyond efficient transformer for long sequence time\-series forecasting\.InProceedings of the AAAI conference on artificial intelligence,Vol\.35,pp\. 11106–11115\.Cited by:[§1](https://arxiv.org/html/2609.25326#S1.p5.1)\.