Covariance Density Neural Networks
Authors: Om Roy
Yashar Moshfeghi
Keith M Smith
Venue: Accepted by TMLR ()

Abstract: Graph neural networks have re-defined how we model and predict on network data but
there lacks a consensus on choosing the correct underlying graph structure on which to
model signals. CoVariance Neural Networks (VNN) address this issue by using the sample
covariance matrix as a Graph Shift Operator (GSO). Here, we improve on the performance
of VNNs by constructing a Density Matrix where we consider the sample Covariance matrix
as a quasi-Hamiltonian of the system in the space of random variables. Crucially, using this
density matrix as the GSO allows components of the data to be extracted at different scales,
allowing enhanced discriminability and performance. We show that this approach allows
explicit control of the stability-discriminability trade-off of the network, provides enhanced
robustness to noise compared to VNNs, and outperforms them in useful real-life applications
where the underlying covariance matrix is informative. In particular, we show that our
model can achieve strong performance in subject-independent Brain Computer Interface
EEG motor imagery classification, outperforming EEGnet while being faster. This shows
how covariance density neural networks provide a basis for the notoriously difficult task of
transferability of BCIs when evaluated on unseen individuals.

Covariance Density Neural Networks
Om Roy!." Yashar Moshfeghi', Keith Smith
Computer and Information Sciences, University of Strathclyde, Glasgow, Scotland, United Kingdom, G1 1XQ
o o.roy 20220umi strath.ac.uk.
Abstract
Graph neural networks have re-defined how we model and predict on network data but there lacks a consensus on choosing the correct underlying graph structure on which to model signals. VoVaiiance Neural Networks (VNN) address this issue by using the sample covariance matrix Graph Shift Operator (GSO). Here, we improve the performance of VNNs by constructing Density Matrix where we consider the sample Covariance matrix as a quasi- Manillnannlnniaan of the system in the space of random variables. Crucially, using this density matrix as the GSO allows components of the data to be extracted at different scales, allowing enhanced discriminability and performance. We show that this approach allows explicit control of the -abblityydiscrmmmnaabbl trade-off of the network, provides enhanced robustness to noise compared to VNNs, and outperforms them in useful real-life applications where the underlying covariance matrix informative. In particular, we show that our model can achieve strong performance in EEG motor imagery classification, outperforming F.EGnet while being faster. This shows how covariance density neural networks provide a basis for the notoriously difficult task of transferability of BCIs when evaluated on unseen individuals, while providing a principled, tuneable control over the stability-disrrimeee discriminability trade-off via the inverse temperature parameter B.
Brain Computer Interface
1
Introduction
Graph Neural Networks (GNN) have served as an essential platform for the modelling of Network data
(Gama et al.
2020 Scarselli
2009
Veličković
2018 Ruiz ct al. 2021 Defferrard et al. 2016 Isufi
et al.
2020 Kipf & Welling 2017
Vignac el al. 2020 Much of the recent. developments in GNNs have
stemmed from Graph Signal Processing (GSP) (Oreea et al. 2018
Sandryhaila & Moura, 2013 and other
2013 Lambiotte arcas of Network Science (Peixoto 2015 Guimerà & Pardo 2009; Domenico et al.
ct al. 2014 Arenas et al. 2006). A key component of GSP is the Graph Shift Operator (GSO), whose
oogen-lcoompsitioon underlies spectral filtering operations analogous to the Discrete Fourier Transform.
However, in many application domains, including neuroscience, finance, and genomics, we observe multivariate signals without known underlying graph topology. Standard GNNs require a pre-specified graph, and learning one from data can introduce instability, overfitting, and noise sensitivity, especially under limited samples. A fundamental limitation of GSP is that the GSO is typically unrelated to the graph signal data itself (mith et al. 2019).
CoVariance Neural Networks (VNNs) (Sihag et al. 2022) address this by using the sample covariance matrix, computed directly from the data, as the GSO. This approach has intuitive links to Principal Component Analysis (PCA) and inherits desirable transferability properties of GNNs to data of different. dimensions (Shag et al. 2022). However, VNNs face two key limitations:
No control over the dialllttydissriiinnbbllyl tradc-off: VNNs inherently discriminate signal differences lying in the low-voriande eigenspace, while high-variance principal components remain indistinguishable. There is mechanism to tune which spectral components the network can resolve.

---

2. Noise sensitivity: Since VNNs discriminate in the low-variance subspace, noise concentrated in these directions common in low signal-to-noise ratio (SNR) settings (such as EEG data) directly degrades performance.
In this work, we propose Covariance Density Neural Networks (CDNNs), which resolve both limitations. We construct a density matrix p(C) = e BC /Tr(e- from the sample covariance matrix C, treating C as a quasi-Hamittonia in the space of random variables. Using this density matrix as the GSO introduces qualitatively new capabilities.
Density matrices have long been used in Quantum Mechanies and information theory (Borst & Theunissen Wilde 2013) to describe physical systems Quantum States (Domenico & Biamonte, 2016 Acín
1999
et al.
2007 Braunstein & Caves 1994. Zdeborová & Krzakala, 2016 Hübner 1992). A density matrix is a
positive semi-definite, self-adjoint matrix with a trace of one acting on the Hilbert Space of the system (Acin eel al. 2007). This gives a probabilistic interpretation to the eigenvalues of the matrix where the eigenmodes represent different probabilistic states while also allowing an entropic interpretation. The Hamiltonian of the system determines the evolution of the system defined by a density matrix. Naturally, the graph Laplacian has appeared as a candidate to represent this Hamiltonian operator. However, the question remains, what. if the true underlying graph structure is unknown?
In this vein, we show that the sample covariance matrix, under reasonable assumptions, can act as spectral surrogate to the Graph Laplacian and thus play the role of a Quasi-Hamillonian We propose a novel density matrix constructed from the sample covariance matrix and define convolutions on this operator as the basis of our Covariance Density Neural Networks (CDNN). As a consequence of this formulation we present novel information theoretic tools such as the multi-scale Von Neumann entropy for Covariance matrices, a measure of entropy that can be applied to singular, low-rank covariance matrices. We also create multi-scale filter banks the sample covariance matrix which show improves performance. Further, we show empirically and theoretically the importance of the 3 (inverse temperature) parameter and how it allows us to control the discriminability and stability of the network.
Empirically, we show that CDNNs outperform VNNs in financial forccasting where the underlying covariance matrix may be informative (Wu et al., 2020). CDNNs also show euhanced performance in the analysis of neurological signals, in particular, we demoustrate strong performance classifying unseen individuals' EEG brain signals in Motor Imagery (MI) tasks (Tangermann et al. 2012 Xu et al. 2022). This leverages the transferability Property of VNNs and combines it with the increased stability and discriminability of CDNNs to classify 4-class motor imagery signals when evaluation is done on a test individual not seen in training "Tangermann et al. 2012). We show that our approach outperforms benchmarks in the field of EEG signal classification such as EEGNet GEGNet/Keutayee Keutayeva et al. 2024) while being faster. This is a step towards real-time
Brain Computer Interfaces (BCI) (Xu et al. 2022 Yang et al., 2021: Zhang & Liu 2018 Zhang et al. 2019
2023
where high accuracy on unseen individuals coupled with low training times are highly desirable.
2
Background and Motivation
Overview
Graph Signal Processing in brief. In classical signal processing, time series is a sequence of values indexed by time, and the shift operator (delay by one sample) underpins the Fourier transform and filtering. Graph Signal Processing (GSP) generalises this idea: graph signal assigns a value to each node of a graph, and the Graph Shift Operator (GSO), typically the adjacency or Laplacian matrix, plays the role of the shift. Multiplying a signal by the GSO "propagates" values along edges, and repeated application (powers of the GSO) defines graph convolution filters that mix information at increasing neighborhoods. The cigen-decomposition of the GSO provides graph Fourier basis analogous to the DFT, enabling spectral analysis and filtering of signals on graphs.
Graph Neural Networks as learnable graph filters. Graph Neural Networks (GNNs) are essentially graph filters with learnable coefficients. typical GNN layer can be written as a (Ek hx,5,)), where S is the GSO, by are trainable filter coefficients, and a(-) is a point-wise non-linearity (e.g., ReI.U). The coefficients
2

---

{h&} determine the spectral response of the filter (how strongly different graph frequencies are amplified attenuated) and are learned from data via backpropagation. Thus, GNNs combine the representational power of polynomial graph filters with the flexibility of learned parameterization and non-linear activation functions.
Covariance as a graph operator. When no graph topology is available, the sample covariance matrix Cn offers a natural data-driven alternative: its (i, entry measures the linear association between variables and j, providing implicit "edge weights." Using Cn as the GSO means that filtering a signal propagates information according to the statistical dependencies in the data, and the cigenvectors of C, which are the principal components, serve the graph Fourier basis.
Why a density matrix? A standard covariance-based based filter can only discriminate signal differences in the low-variance eigenspace (see Theorem 1). The density matrix transformation p(C) ( -TT((e=TTT inverts this: for 3 > 0, large covariance eigenvalues map to small density eigenvalues, shifting high-variance components into the discriminable region. The paramcter thus provides explicit, tuncable control over which spectral components the network vesolves and provides noise-disrinninabiity tradeofT. By combining multiple 3 values in filter bank, CDNNs can discriminate across all spectral scales simultaneously. The normalization by the trace induces scale invariance and stability to the operator while allowing interpretability measurcs.
and transfer of information-theom
We now introduce formal definitions that link these frameworks together.
Definition 1 (Graph Convolutional Filter). Let h = //o....,,x/! T be filter coefficients. Let be the GSO. A graph convolutional filter of order K is the linear map
K
H(x)
hx S H(S)x, x.
k=
where H(S) = >k Ex S*
Definition& 2 (Graph Fourier Transform (GFT) Isufi et al. 2024 For diogonalizable GSO S VA V-1 with eigenvectors V and eigenvalues A, the GFT of a graph signal æ is == V x, and the inverse GFT is
x v
Definition (Covariance and Sample Covariance). Let X E Rdd with mean μr. The covariance matrir is given by
C= B x) 1.
(1)
From n samples x(1) ...., x the sample covoriance is defined as
Cn
X X)
(2)
where
X
(3)
As n + 00, the sample covoriance converges to the true covariance: C₂ + C.
Definition 4 (coVariance Fourier Transform (VFT) (Sihag et al. 2022)). Let C = UEU be the eigendecomposiiionn of the sample coveriance. Then the VFT of a random sample æ is U n.
Quantum mechanical density matrices establish probability structure for system states (Rao 1945 Domenico & Biamonte, 2016). They must be Hermitian (p = p!), positive semi-definite (p 2 0), and have unit trace (Tr(p) =). Importantly, the following density matrix has recently been Proposed Domenico & Biamonte 2016
-L
Z Tre
Z
where L is the graph Taplacian, and B > 0 is a parameter controlling the diffusion.
3

---

This allows the calculation of the spectral entropy of these matrices where eigenmodes define probabilistic routes of diffusion and B allows the calculation of the entropy at different scales. Importantly, the sub- additivity of entropy is also maintained, a feat not achieved previously (Araki & Lieb 1970 Anand et al.) 2011). Taking inspiration from this we observe that there is a distinct relationship between graph Laplacian's and inverse covariance (precision) matrices and that graph learning is linked to sparse covariance estiimation (Dittrich & Matz, 2020).
It is established that the estimation of the combinatorial graph Laplacian (CGL) matrix from observed signals is usually performed via the CGL estimator (CGLE) Pavez 2022.. Specifically, the estimated Laplacian matrix L is obtained as
È arg min LEC,(E)
(-log log det (L) tr(LC.)),
(4)
where Ep) denotes the set of valid combinatorial graph Taplacian matrices consistent with the edge set E,
and Cn is the sample covariance matrix the observed graph signals.
Note that the estimation problem in equation is grounded in tight spectral connection between the graph Laplacian and the covariance of graph-stationas signals. Under the Gaussian-Marker assumption C, L the covariance can be expressed as spectral function of L and therefore shares the same eigenvectors (Dong ct al. 2019: Pavez & Ortega 2016). This means that, when the number of observations greatly exceeds the graph order N, the sample covariance C, forms a consistent surrogate for the Taplacian eigenspace. The CGLE in equation g can thus be interpreted as selecting edge weights so that L aligns with the empirical second-order statistics while retaining the common eigenbasis. It is thus intuitive to interpret the sample covariance matrix as signed or dual data-driven Laplacian. We refer the reader to Appendix A.6 for more details.
It has also been shown that the eigenvectors of the Laplacian optimally decorrelate signals in a Gaussian Markov random field (GMRF) model where the precision matrix (inverse covariance matrix) is defined as the Laplacian (Zhang & Florencio 2013). This Property makes precision matrix eigenvectors valuable for tasks like signal compression. However, the covariance matrix C shares the same eigenvectors as the precision matrix (up to eigenvalue scaling), retaining the decorrelation Properties without requiring explicit, inversion of the covariance matrix. Moreover, C directly encodes global pairwise dependencies and serves more practical graph shift operator in graph neural networks (GNNs) (Keriven & Peyré 2019 Zügner & Günnemann 2019: Isuli et al. 2024 Bruna et al. 2014), as it avoids the sparsily constraints of the precision matrix and the numerical instability associated with its estimation (lou
2017
Mestre 2008
A real signal æ; on a graph has Dirichlet energy
x Lxx = Σ
wisla; - x;)² Σ day +
(-
Σ Wij a.T;);
(iE))E
(1.))E
self
interaction
where d, = E, Wiy is the vertex degree. For centered random variables x; with covariance matrix C,
x TCxx= Σ Var(z.) of
Σ (o)
the diagonal variances Var(æ;) supply node potentials while off-dingonal covariances provide pairwise couplings. Thus C may be viewed as Hamiittnmmtoononn its spectrum governs the same collective modes that minimise Dirichlet energy for L.
In certain use cases when the covariance matrix is used as a "graph" constructed from the data itself (i.e brain networks) the conceptualization of the covariance matrix as a "qumitonaamllnniann" makes intuitive sense
We emphasisc that the term "quasi-Hamiltoniie is used spectral analogy highlighting the rolc of C in shaping the system's spectral modes and governing sigual diffusion,rather than implying strict quantun-mechanical interpretation. The analogy is grounded in the shared mathematical structure: just as the Hamiltonian determines the Gibbs state in statistical mechanics, C determines our density operator.
4

---

and may have practical benefits (the Hamiltonian is usually considered the energy operator of the system in quantum mechanics). Furthermore, when there is no predefined underlying graph structure the covariance matrix can act as a spectral surrogate to the graph Laplacian. This is the motivation of Covariance Density Neural Networks.
Takeaway
When no graph topology is available, we demonstrate that the sample covariance matrix can serve as spectral proxy for the true graph Laplacian, effectively functioning as a quasi-Hamiltonn or energy operator our density-matrix formulation.
3 Covariance Density Neural Networks
In this section we introduce the theoretical foundations of Covariance Density Neural Networks. We first introduce the Covariance Density Matrix and Filter.
Definition (Covariance Density Operator and Filter). Define the density operator associated with a covoriance mairia C' as:
BBC
pl CC Tr(e-Do)-
For an input vector x and real-valued parameters {ha] the output of the Covariance Density filter is defined as:
K z= H(H(p(C))=== Σ hp lo 6)
Analogously, can define the Covariance Density Perceptron and its filter bank representation follows: Definition 6 (Covariance Density Perceptron). Consider Covariance Density filter Hp( C)) and given nonlinear activation function o-)). For an input x, the Covariance Density Perceptron is defined as:
(ax; C, H o
Definition 7 (Multi-Scale Covariance Density Filter Bank/ Bank/Layeri Layer). For covariance density filter bank consisting of Fout scales with Fin m-dimensional inputs for each 8m the m-th scale is defined as:
Fin
Doue[m]] =0 Σ H (p(C)) (sns(g) (0]
and the fonal malli-scale oulpul of the layer is
Ffsnal A A(Fowe[1), Dous [2], .... Fous Fous| Fout)),
where A an aggregation operator (e.g., concatenation, summation, or mean).
order to show that graph neural networks are stable, the assumptions of integral Lipschitz continuity of the network's GSO frequency response is traditionally assumed. For Covariance Density Networks however, explicitly derive the Composite Lipschitz constant with respect to the density transformation is instructive for explaining the role of the inverse temperature 3 in the network.
2While they exist, graph learning pruning approaches can introduce instability, over-fitting, and noise sensitivity, especially when the true connectivity is uncertain or the number of samples is limited. By contrast, using the sample covariance (and thus the density matrix) as fixed surrogate (similar to the Laplacian) provides stable, interpretable operator without the risk of learning spurious edges.
5

---

Theorem 1 (Composite Lipschitz Conditions for covariance density Filters). The composite frequency response of Covoriance Densily Filler is given by
K Σ
M((())
Z
where hgg are finite coefficients, is tunable parameter, Zk = 6- the partátion function,
and ( are the eigenvalues of the covariance matrix C. Assuming and hy are fovile constants, the coveriance density filter salisfies:
The Lipschitz condition:
[h(p(22) - h(p(Aa)) Sal22 A
The Integrol Lipschitz condition:
(h(o(A2)) -(((A)) 0 A1
0
a
sup (Aitt
Por some with Σk= |k11133
Max Pairwise Difference vs B
Importantly, the theorem holds for both positive and negative 8. We assumee composite Lipschitz bounds on the density-transforned eigenvalues, treating filter coefficients as fixed, even though learned filters may violate integral Tipschitz. Viewing p(C') as GSO, linear filters cannot remain stable and distinguish signals beyond their Lipschitz cutoff mutoff-iee. high-variance principal com- ponents stay indistinguishable. Pointwise nonlinearities can help (Ruiz et al. 2020 but some high-frequency content.
0.025
Max (h(p())) h(p(A)))
Max (h(p(A,) - -(((A)))
atA2 A1l
0.020
0.015
0.010
0.005
0.000
-1.00
-0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
B
Figure 1: Composite Lipschitz Condition Empirical Validation
still eludes discrimination [Pfrommer et al. 2021
Observe, in the covariance donsity context, that the cigenvalues of the covariance matrix C determine high-variance vs low-variance directions in the eigenbasis. we denole these eigenvalues by (A.) ascending order (i.e. Å₁ S < AN), then the density operator has eigenvalues {pr} that are decreasing functions of X when 3 0 and increasing when <.
For the Covariance matrix itself larger X; corresponds to higher variance directions. Consequently, for the covariance density filter (p(C)), the composite integral Lipschitz condition is equivalent to
(A1 + Azl
|p: (M(p(A))) < [3] Σ ha [kksspp
3 "Alhhoggh is rather un-intuitive consider negative values B, allowing its use allows greater exploration the spectral profile, nevertheless such interpretations are not missing in the literature, consider the communicability, matrix introduced by Estrada et al. (Estrada et al. 2012) or even the Softmax function (LeCun et al. 2015 Goodfellow et al. 2016).

---

this forces the derivative of the filter to shrink as Ps grows. Consider the case when B 0, larger x leads to smaller Pi, so one sees a role reversal, i.e. large X are mapped to small density-eigenvalue Px- and small A are mapped to large density-eigenvalue
Hence, imposing an integral Lipschitz bound in the covariance density setting tends to flatten the filter for the large values of Priee.
Lee. for the low-variance eigenvalues A₃-
This leads to a different notion "cut-off" frequency: instead of failing to distinguish signals whose dif- ference lie in the highest-variance subspace, here the filter has trouble distinguishing components ly- ing in the lowest-voriance subspace. These non- discriminable subspaces are spanned by the eigen- vectors whose variance in the covariance matrix is lower. In a PCA interpretation, the high- varianae principal components become the "low-p; modes' and the covariance-density filter is able to discriminate signals that differ primarily in these varianee directions.
Decision boundaries in entropy space
1.-
1.0-
0.8
The trade-off between stability and discriminability thus flips, precisely because X and Ps are inversely related under the inverse exponential map.
0.6
LDA on
LDA on
For positive 3 the exponential map ensures a 'com- pression' of high variance components into the dis- criminable region of the integral Lipschitz filter while for 3 < we return to the original covariance sctting where signal differences that lie the low variance components eigenbasis can be discriminated but this time with larger eigenvalue gap for higher variance components. i.c. the filter behaves similarly to covariance filter. Importantly we can now construct multi-scale filter banks (and thus networks) that are discriminative in both high and low frequencies while being simultaneously more stable, although this in- creased discriminability could introduce noise into the system. The set of discrèminable signals for a
0.4
20
0.60
0.70
0.75
0.80
0.85
0.90
0.95
Spaiee
Figure Discrimination of two singular Gaus- sians using CVNE. We sampled 5= triples (X,Y, 2), estimated covariances on sliding windows (W=128) and computed Snaive, S3 (8 2). Back- ground: 2D LDA decision regions. Dashed: best 1D splits Naive fails (AUC 0.5), von Neumann succeeds (AUC 1.0).
CDNN is thus Pfrrommer et al. 2021 ::
((x,)) E RN ((x:((C)) - (y:: E (How-2)},
8> 0,
D+
((x,))R)) ((x: C) (y: C) E (Vnigh->)) B << 0.
where (Vic) corresponds the eigenspace associated to those As these two spaces are orthogonal complements of each other, CDNNs constructed using multi-scale filter banks consisting of positive and negative B values can discriminate signal differences in both high and low frequencies.
Ps above the relevant cut-off.
In this sense, it is also intuitive to allow the network to learn 8 values from the data itself. In this vein, we show that there is no significant drop in performance when this parameter learned from the data and not hand-picked using domain knowledge (Appendix A.3.1
As shown in Theorem the Lipschitz constant of the filter frequency response is explicitly controlled by B and the filter coefficients. This is particularly useful as discriminability of the eigenvalues of the covariance matrix and the filters stability by tuning 3.
directly manage the trade-off between the
7

---

Another benefit of the density matrix formulation is that it enables the use of information-theoretiee methods.
Given that the eigenvalues of the density matrix now form probability distribution (i.e. summing to 1) we can now define multi-scale Von Neumann Entropy for Covariance matrices (CVNE).
This allows us to determine the 'regularity' of covariance matrix. we show several desirable properties of this formula- tion, such as its ability to provide entropy values for low-rank, singular matrices, thus the CVNE and consequently CDNN's behave as inherent regularizers, stabilizing ill-conditioned co- variance structures. We also prove the sub-additivity of CVNE (Appendix A.7) and show cases where it can be discriminalive compared to naively normalizing the covariance matrix by its trace. For example, such an entropy is less effective when near global variance fluctuations occur (Appendix A.8 As the limiting extreme case of pure global scaling Figure 2 shows that it is in fact scale-blind.
Multiscale von Neumann Entropy for Co- variance Matrices
Let C RNXN be positive semidefinite co- variance matrix. The Multiscale von Neuman Entropy of C is
Sp(C) = TT(((((C) log log)(
For completeness Theorem 2 shows that Covariance Density Neural networks are permutation equivariant, another typical requirement of suitable GNN.
Theorem 2 (Permutation Equivariance of the Covariance Density Filter). For any permutation matrix T E IR define the permauled Covariance matric C T CT and the permuled signal is == T Then the filter H(p) is permutation equivariant, i.e.,
H((p(O)== T H((p(
While the sensitivity of the ensemble covariance matrix to perturbations well-studied we must derive new bounds in this new covariance density reprosentation.
For this, we begin by relating the operator norm the covariance matrix with the covariance density matrix. Lemma (Error Bound for Covariance Density Matrix). Lel C be covariance malria and 5C ils perturbation. Further, let E be the density 'error matriz', then
|&] || F(3, F(3, (1 (1+me (1teoool1111))))
E
R
where R > 0, m = dCm, 1/200 is the indicator of 3 < 0, and
F(8,C,8C) { (100000 -1)
B2 0,
|8| [I|C + 8C|| - [[]] <<0,
with limg-o F(3, C, 1.
Theorem 3 (Stability of Covariance Density Network). Consider an L-layer, F-channel mulli-scale CDNN o(x: P N. H) with Lipschitz constants {a,} and inverse-temperahures {B,). Let Q > O be constant. Assume [h (p) < is 1-Lipschitz, and let k be as in Sihag et al. 2022). For any E (0, 2 with probability
2 1 2e 2km/n,
penalty for B<0
F (x:: - (, < L (II
a4 Q( MC R
F(B,,C,0C
1=1
deneity estimation ervor penalty for 3<0
X m (1+ m
10 (8.IIIC)
m2
O(n 1n-
4-1
a2
density estimation error
dimensionality impact
From the stability bound, the terms F(3, C, 8C) and exp (ι<о) [[][|[]] quantify penalties associated with negative B in the uncertainty in the covariance density estimate, amplifying the covariance uncertainty. These
8

---

terms vanish for positive beta indicating the network is more stable for positive values of beta. As expected, the uncertainty of the covariance density estimate is controlled by the Lipschitz constant which, recall, is controlled completely by the filter coefficients and beta. Thus as beta tends to 0 the network becomes increasingly stable. However each filter in the filter bank of size F can have a different 3, thus including larger or negalive values of will reduce the stability of the network. Also note, that even if there are large errors the estimation of the true covariance matrix (indicated by &C) this can be offset by using smaller values of 3, confirming the de-noising effect. The dimensionality impact term ni reflects the relationship between the sample size n, the dimensionality and the convergence rate.
Takcaway
We propose multi-scale Covariance Density Neural Networks (CDNNs), show that B controls the slabilly-disccininablly tradeoff of the network, and formulate the Multiscale Von Neumann Entropy for Covariance Matrices.
4 Experimental Results
4.1 Covariance Matrix Construction
Financial data. Each asset/featur constitutes a random variable (uode), and each time step within a sliding window constitutes one observation. The sample covariance matrix is computed once over the full training set and captures cross-asset correlations.
EEG data, Each EEG channel is random variable (node), and each time step in trial constitutes one observation. The covariance matrix is computed once over all training subjects' trials timesteps, repeated for all participants and averaged over, producing a global cross-channel covariance that captures shared spatial patterns across individuals.
4.2 Financial Forecasting
Covariance matrices have been deemed useful in predicting financial asset behaviour, offering insights into asset interactions and enabling better risk management. strategies. Building on this, covariance-basel neural networks show significant promise, this promise however, largely affected by the low signal to noise ratio (SNR) present in financial data and thus the estimates of it's underlying covariance structure.
The S&PP 500 dataset (Jander 2023 spans 10 years of multivariate time series data, including features like stock volatility and daily option volumes. The Exchange Rate dataset (Wu et al. 2020) consists of daily exchange rates for different currencies across 7588 open days. Covering the period from January 2, 2020, to February 2, 2024, the US Stock dataset (Kumar 2024) includes various variables across major stocks, indices, commoditics, and cryptocurrencies, offering insights into pricing and trading activity.
Table 1 compares traditional graph-basee models (e.g., Graph Attention Networks (GAT), Graph Isomorpshism Network (GLN), Graph Convolutional Networks (GCN))) with covariance- based models (VNN, CDNN). We also include a hybrid model that first applics the covariance density filter at multiple scales before applying an attention layer.
We attribute the improved performance of CDNN over a standard VNN to the versatility of the multi-scale filter bank, capturing information different seales. Over all datascts, either the CDNN or hybrid model performed the best. This indicates that the added covariance information indeed informative and that enhancing attention based models with this information improves performance.
The robustness of covariance models likely stems from their ability to leverage stable cross-asset correlations driven by macroeconomic or sector factors. This inductive bias makes them inherently risk-aware, ensuring better diversification and consideration of correlated risks. In addition, the noise-robustness of CDNNs shows improved performance over all datasets.

---

Table 1: Mean Absolute Error (MAE) + std dev over seeds for horizons 1, 3, 5. Best; second best.
Data
Mthd
H1
H3
H5
VNN
0.1336+00.0075
0.1786+0.0077 0.1677=0.0212
CDNN
0.1102+00.008
0. 0.1291 0065
0. 0.159-+0.0137
GCN
0.2415-10.0100
0.2334 0.0034
0.237110.0177
Exchange
GTN
0.1939++0.0088.
0.1950-0.0031
0.2073++0.015
GAT
0.1819:10.028
0.1852-0.0230
0.2075.10. 0228
Hybrid
0.1057=00.073
0.1171=0.0120
0.126710.0077
VNN
0.6198+0.0158
0.6392-0.0016
0.6553+0.0008
CDNN
0.5848+000033
0.5983=000049
0.6090=000006
GCN
0.603110.0193
0.639
0.6719:100009
S&P 500
GIN
0.6682++0.0115
0.7014=0.0123
0.7103+0.000
GAT
0.5930=t0.0245
0.6309=0.015
0.503++00019
Hybrid
0.5772-10.0039
O. 6.6000 0075
0.611140. 0027
VNN
0.424410.0159
18822=000088
0. 0. 0063
CDNN
0.3412++0.011
0.4375=0.0129
0.5290=
GCN
0.5572++0.0366
0.5934-0.0255
0.6592+0.0015
US Stock
GIN 0.6332++0, 1248
0.6627-0.0935
0.7113+0.090
GAT 0.43466+0..555
0.5141-0.0542
0.5832++0.02
Hybrid
0.354240.000
0.450610.0015
0.554710.0211
4.3
Transferability of covariance density Networks in Brain Computer Interfaces
Brain-Computer Interfaces (BCIs) have long faced challenges in transferability (Zhang et al.) 2019 Zhang & Liu 2018 Yang et al., 2021 i.e., ensuring model trained on multiple individuals genera lizes to unseen individuals' brain signals. The sample covariance matrix computed purely over training individuals offers a promising solution by capturing stable, global correlations common across individuals. Using this as a shift operator for unseen individuals allows individual differences to be "filtered" by these shared patterns (Roy et al. 2024 2023), this can be best related to the neuro-scientific concept known as functional alignment Bazeille et al. 2021).
Our model is bench-marked on two BCI tasks: the 4-class BCI-2A motor imagery dataset (9 subjects) and the PhysioNet 2-class dataset (105 subjects). For BCI-2A, we train on data from all but one held-out subject, rotating through each as the test individual. On PhysioNet, we evaluate across 10 non-overlapping folds, ensuring each fold's subjects are unseen during training. In both cases, we compute a global covariance matrix from the training subjects and use it as a graph-shift operator to classify the held-out test data.
Table 2: Average Classification accuracies (%) + std dev over subjects/fole /foldsss//ssssssssosss and training speed (s/epoch) on
BCI 2A (4 Class). Best; Second best.
Method
2A (2)
2A (1)
PhysioNet
Average
Speed
GIN
51.1 + 3.33
32,6 + 2.78
78.9 + 4.03
54.20
0.91
GAT
62.67 + 9.36
44.1 + 6.77
81.0 + 5.02
64.10
2.20
VNN
59.3 + 1.86
30.1 + 2.72
79.4 + 3.01
56.30
0.90
CDNN
73.4 + 8.53
50.4 + 7.13
80.02 4 2.74
67.90
1.10
EEGNet
61.4 6.09
45.2 5.47
75.1 + 0.29
60.60
4.40
CSP+L.DA
56.1 + 4.21
33.37
+ 8.99
52.08
+ 3.86
47.20
N/A
Table 22 presents our results, comparing CDNNs with a VNN, more traditional Graph Based approaches such GIN and GAT, CSP approach and the lightweight EEGNet benchmark. For each graph model, we take the convolved signals one time series per channel and concalenate every channel's full sequence of time points into single long feature vector in order to retain the high temporal resolution of EEG signals. That flattened vector is then fed into a simple fully-connected network (MLP), and we apply a Tanh nonlinearity followed by the classification layer.
The classical CSP + LDA approach performed the worst, EEGnet in fact underperformed on the PhysioNet dataset compared simpler graph based models highlighting that, for iedeecc-inddeendent BCIs, less complex models are preferable to avoid over-fitting. While GAT tops the PhysioNet leader board, CDNNs are the strongest model overall, outdoing VNN both in accuracy and efficiency (1.1 s/epoch vs. E.GGeet's 4.5 s), making them a suitable choice for rapid retraining and real-world BCT deployment.
10

---

Density Matrix (ß (B=0.1)
Density Matrix (B=5.0) Entropy: 2.908
Density Matrix 15.1)
Global Covariance
Entropy: 3.071
Entropy: 2.812
Motrix Value
Donsity Matrix (B=0.1)
Density Matrix (B=5.0) Entropy: 2.938
Density Matrix (B=15.1)
Unseen Covariance
Entropy:
Entropy: 2.853
Figure 3: Contour Plot showing distribution of Covariance and Covariance Density Matrix on the Scalp at different scales. The top row indicates the averaged training covariance matrix (Excluding the Unseen Individual) and the bottom row indicates the unseen Individuals test. covariance matrix. We can see that the training covariance matrix relatively similar to the Unseen Individual and that different scales unveil further potential similarities thus increasing classification performance
The large gains the BCI-2A 2-class and 4-class datasets explained by the multi-scale filter bank's ability to capture patterns at different spectral scales of the cross-subject covariance matrix (Figure 3) while also being robust to noise at positive values of B (Theorem 3).
4.4 Computational cost.
The primary additional cost of CDNNs over VNNs is the matrix exponential exp(-9C,.). Since the covariance matrix is fixed, this is computed once via eigen-decomposition and backward passes have identical per-step cost to VNNs: O(Km²) per order-K filter, where m is the number of channels/ baaannels/vrrablles variables. Memory overhead is a single additional m x m matrix per 3 scale. Empirically, with caching enabled, CDNN runtime matches VNN even for m > 1000 (see Appendix A.3, Figure 7). Compared non-graph basclines, CDNN trains at s.1s/epoch on BCI-2A versus EFGNetts 4.4s/epoch- 4x speedup (Table 2).
O(m") cost and cached. Subsequent forward
5 Conclusions and Limitations
We introduce Covariance Density Neural Networks (CDNNs), which treat the covariance matrix as a Hamiltonian-like object build density matrix that acts as graph shift operator a Graph Neural Network. We showed that the sample Covariance Matrix act spectral surrogate to the true graph Laplacian and thus constructed novel density matrix, allowing us to use lools from information theory to define a new entropy measure for covariance matrices with notable benefits. We've also established rigorous stability bounds for the Covariance Density Filter and the network itself, showing how 3 explicitly controls the stability of the network.
In real-world applications where covariance structures matter (o.g., financial or neurological networks), CDNNs show strong performance against traditional benchmarks and deliver significant gains over VNNs. In subject-independeen motor imagery classification for Brain-Computer Interfaces (BCIs), they outperform FEGNet while being significantly faster. Although we relied on simple low-parameter architectures here, we believe that folding CDNNs into more complex neural network architectures can yield further benefits.
11

---

There are some clear limitations in the use of the Covariance matrix as it only encodes second-order statistics, it misses higher-order non-linear dependencies. Future work could integrate mutual-information matrices or Gaussian-process kernels into the density matrix formulation. Furthermore, time-varying or non-stationary signals (e.g. drift in EEG) need dynamic covariance models or online adaptation and exploration beyond static matrices should be explored and can be particularly beneficial in the study of Dynamic Functional Connectivity (DFC) the brain. We have also consider that the computational cost scales with the size of the covariance matrix and matrix exponentials thus may become harder to compute. Since we only need the product of the matrix-vector multiplication approaches such as Krylov Subspace Approximation methods may be useful.
6 Impact Statement
This work may enhance Brain-Computer Interfaces (BCIs) by enabling more reliable brain signal decoding for assistive technologies. However, ensuring secure, anonymized handling of neurological data and addressing fairness across diverse users remain important challenges to prevent biases in medical and neuro-technological applications.
References
How close are the eigenvectors of the sample and actual covariance matrices? In International Conference on Machine Learning, pp. 2228-2237. PMLR, 2017.
Antonio Acín, J. Ignacio Cirac, and Macicj Lewenstein. Entanglement percolation in quantum networks
Nature Physics, 3:256, 2007.
Kartik Anand, Ginestra Bianconi, and Simone Severini. Shannon and von neumann entropy of random networks with heterogeneous expected degree. Physical Review E, 83:036109, 2011.
Huzihiro Araki and Elliott H. Lieb. Entropy inequalities. Communications in Mathematical Physics, 18:160, 1970.
Alex Arenas, Albert Díaz-Guilera, and Conrad J. Vicente Synchronization reveals topological scales in complex networks. Physical Review Letters, 96:114102, 2006.
Bruno Aristimunha, Igor Carrara, Pierre Guetschel, Sara Sedlar, Pedro Rodrigues, Jan Sosulski, Divyesh Narayanan, Erik Bjarcholt, Quentin Barthelemy, Robin T. Schirrmeister, Emmanuel Kalunga, Ludovic Darmet, Callan Gregoire, Ali Abdul Hussain, Ramiro Gatti, Vladislav Goncharenko, Jordy Thielen, Thomas Moreau, Yannick Roy, Vinay Jayaram, Alexandre Barachant, and Sylvain Chevallier. Mother of all BCI benchmarks (MOABB), 2023.
Thomas Bazeille, Elizabeth DuPre, Hugo Richard, Jean-Baptiste Poline, and Bertrand Thirion. An empirical evaluation of functional alignment using inter-subject decoding. Neurolmage, 245:118683, 2021. doi: 2...0066..euroimmge2221...11.
Alexander Borst and Frederic Thounissen. Information theory and neural coding. Nature Neuroscience, 2: 947, 1999.
Samuel L. Braunstein and Carlton M. Caves. Statistical distance and the geometry of quantum states. Physical Review Letters, 72:3439, 1994.
Leo Breiman. Bagging predictors. Machine Learning, 24(2):123-140, 1996.
Joan Bruna, Wojeiech Zaremba, Arthur Szlam, and Yann LeCun. Spectral networks and deep locally connected networks graphs. In 2nd International Conference Learning Representations, pp. 1-14, Banff, AB, April 2014.
12

---

Andrea Cavallo, Maosheng Sabbaqi, and Elvin Isufi. Spatiotemporal covariance neural networks. In Joint European Conference on Machine Learning and Knowledge Discovery in Dalabases, pp. 18 18-34, Cham, August 2024. Springer Nature Switzerland.
Michaöl Defferrard, Xavier Bresson, and Pierre Vandergheynst. Convolutional neural networks graphs with fast localized spectral filtering. In 30th Conference on Neural Information Processing Systems, pp- 3844 3858, Barcelona, Spain, December 2016. Neural Information Processing Foundation.
Thomas Dittrich and Gerald Matz. Signal processing on signed graphs: Fundamentals and potentials. IEEE Signal Processing Magazine, 37(6):86 98, 2020.
Manlio De Domenico and Jacob Biamonte. Spectral entropies as information--heoreliic lools for complex network comparison. Physical Review X, 6(4):041062, 2016. doi: 0..108//hysReee...011000...
Manlio De Domenico, Albert Solè-Ribalta, Emanuele Cozzo, Mikko Kivelä, Yamir Moreno, Mason A. Porter, Sergio Gòmez, and Alex Arenas, Mathematical formulation of multilayer networks. Physical Review X, 3: 041022, 2013.
Xiaowen Dong, Dorina Thanou, Michael Rabbat, and Pascal Frossard. Learning graphs from data: A signal representation perspective. IEEE Signal Processing Magazine, 36(3):44 63, May 2019. doi: 10.1109/MSS MSP. 2018 28.724.
Ernesto Estrada, Naomichi Hatano, and Michele Benzi. The physics of communicability in complex networks. Physics Reports, 514:89, 2012.
Augusto Fontan and Claudio Altafini. On the properties of Laplacian pseudoinverses. arXiv:2109.14587 2021.
Fernando Gama, Elvin Isufi, Geert Leus, and Alejandro Ribeiro. Graphs, convolutions, and neural networks:
From graph filters to graph neural networks. IEEE Signal Processing Magazine, 37(6):128-14) 138, 2020.
Ary L. Goldberger, Luis A. Amaral, Leon Glass, Jeffrey M. Hausdorff, Plamen Ch. Ivanov, Roger G. Mark, Joseph E. Mietus, George B. Moody, Chung-Kang Peng, and H. Eugene Stanley. PhysioBank, Toylioo and PhysioNet: components of new research resource for complex physiologic signals. Circulation, 101 (23):e215-2222 e220, 2000.
Gene H. Golub and Charles F. Van Loan, Matriz computations. JHU Press, 2013.
Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. The Adaptive Computation and Machine Learning Series. The MIT Press, Cambridge, MA, 2016.
Roger Guimerà and Marta networks. Proceedings of the National Academy of Sciences, 106:22073, 2009.
Missing and spurious interactions and the reconstruction of complex
Manfred Hübner. Explicit computation of the bures distance for density matrices. Physics Letters A, 163:239, 1992.
Elvin Tsufi, Fernando Gama, and Alejandro Ribeiro. EdgeNets: Edge varying graph neural networks. arXiv:2001.06620v2 [cs.LG], March 2020.
Elvin Isufi, Fernando Gama, David 1. Shuman, and Santiago Segarra. Graph filters for signal processing and machine learning on graphs. IEEE Transactions on Signal Processing, pp. 1-32, 2024. doi: 10.1109/TSP.
2024.3349788.
Mathis Jander. S&P500 volatility prediction time series data. Kaggle Dataset, 2023.
Tosio Kato. Perturbation Theory for Linear Operators. Springer-Verlat Verlag, Berlin, Heidelberg, New York, 1966.
Nicolas Keriven and Gabriel Peyré. Universal invariant and equivariant graph neural networks. In 33rd Conference on Neural Information Processing Systems, pp. 7092-7101, Vancouver, BC, December 2019. Neural Information Processing Systems Foundation.
13

---

A. Keutayeva, N. Fakhrutdinov, and B. Abibullaev. Compact convolutional transformer for subject- independent motor imagery EEG-bbased BCTs. Scientific Reports, 14:25775, 2024. doi: 10.1038/ 041590-004-73755-4
Diederik P. Kingma and Jimmy Lei Ba. ADAM: A method for stochastic optimization. In 3rd International Conference on Learning Representations, pp. 1-15, San Diego, CA, May 2015.
Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In 5th International Conference on Learning Representations, pp. 1-14, Toulon, France, April 2017.
Saket Kumar. 2019-2024 US stock market data. Kaggle Dataset, 2024.
Renaud Lambiotte, Chan-Charles Delvenne, and Mauricio Barahona. Random walks, markov processes and the multiscale modular organization of complex networks. IEEE Transactions Network Science ond Engineering, 1 1:76, 2014.
Vernon Lawhern, Amelia Solon, Nicholas R. Waylowich, Stephen M. Gordon, Chou Po Hung, and Brent J. Lance. EEGNet: compact convolutional neural network for EEG-based brain brain-computer interfaces. Journal of Neural Engineering, 15, 2018.
Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. Nature, 521(7553)::8 117, 2015.
Xavier Mestre. On the asymptotic behavior of the sample estimates of eigenvalues and eigenvectors of covariance matrices. JEEE Transactions on Signal Processing, (6(11):53333-5))
2008.
Antonio Ortega, Pascal Frossard, Jelena Kovačevié, José M. F. Moura, and Pierre Vandergheynst. Graph signal processing: Overview, challenges, and applications. Proceedings of the IEEE, 106(5):808-88)) 828, May 2018. doi: 10. 2000://.PPP...22.22222222222
Eduardo Pavez. Laplacian constrained precision matrix estimation: Existence and high dimensional consistency. In Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, 2022.
Eduardo Pavez and Antomio Ortega. Generalized Laplacian precision matrix estimation for graph signal processing. In Proceedings of IEEE ICASSP, pp. 6350-6354 Shanghai, 2016. doi: 10.1109/ICASSP:200 7472899.
Tiago P. Pcixoto. Model selection and hypothesis testing groups. Physical Review X, 5:011033, 2015.
for large-scale network models with overlapping
Stefan Pfrommer, Alejandro Ribeiro, and Fernando Gama. Discriminability of single-layer graph neural networks. In ICASSP 2021- 2021 IEEE International Conference Acoustics, Speech and Signal Processing Toronto, ON, Canada, 2021. doi: 10. TCA10///CASPP377222.02 202 202.9941588 .9144583.
(ICASSP), pp.
Calyampudi Radhakrishna Rao. Information and accuracy attainable in the estimation statistical parameters. Bulletin of the Calcutta Mathematical Society, 37:81, 1945.
Om Roy, Yashar Moshfeghi, Agustin Tbanez, Francisco Lopera, Mario A. Parra, and Keith M. Smith. Robust, high temporal-rrsolttion EEG functional connectivity detects increased connectivity coinciding with P300 visual short-term memory binding in both familial and sporadic prodromal Alhheime's disease. In Complez Networks 2023: The 12th International Conference on Complex Networks and Their Applications, pp. 679-682, 2023.
Om Roy, Yashar Moshfeghi, Agustin Ibanez, Francisco Lopera, Mario A. Parra, and Keith M. Smith. FAST functional connectivity implicates P300 connectivity working memory deficits Alzheimer's 's disease. Network Neuroscience, 8(4):1467-149) 2024.
Luana Ruiz, Fernando Gama, Antonio G. Marques, and Alejandro Ribeiro. Invariance preserving localized activation functions for graph neural networks. IEEE Transactions on Signal Processing, 68(1):127-1L) 141, January 2020.
14

---

Luana Ruiz, Fernando Gama, and Alejandro Ribeiro. Graph neural networks: architectures, stability, and transferability. Proceedings of the IEEE, 109(5):660-68)) 2021.
Aliaksei Sandryhaila and José M. F. Moura. Discrete signal processing on graphs. IEEE Transactions on Signal Processing, 61 (7):1644-1165, April 2013.
Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. The graph neural network model. IEEE Transactions on Newral Networks, 20(1):61-8)) 80, 2009.
Saurabh Sihag, Gonzalo Mateos, Corey McMillan, and Alejandro Ribeiro. CoVariance neural networks. In Proceedings of the 36th International Conference on Newral Information Processing Systems (NIPS '22), pp. 17003 17016, Red Hook, NY, USA, 2022. Curran Associates Inc.
Keith Smith, Loukas Spyrou, and Javier Escudero. Graph-variate signal analysis. IEEE Transactions Signal Processing, (722:29333306 2019. doi: 10.1 1011199TTPPP20
/T/.. 201.2008..200 2881658.
Jascha Sohl-Dickstein, Peter Battaglino, and Michael R. DeWeese. Minimum probability flow learning. In Proceedings of the 28th International Conference on Machine Learning, pp. 905 912, June 2011.
M. Tangermann, K. R. Müller, A. Aertsen, et al. Review of the BCT competition IV. Fronliers in Neuroscience, pp. 55, 2012.
Petar Veličkovié, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks. In 6th International Conference on Learning Representations, pp. 1-12, Vancouver, BC, April May 2018.
Roman Vershynin. High-dimensiones probability: An introduction with applications in data science, volume 47. Cambridge University Press, 2018.
Clément Vignac, Guillerno Ortiz-Jiménez, and Pascal Frossard. On the choice of graph neural network architectures. In 45th IEEE International Conference on Acousties, Speech and Signal Processing, pp. Barcclona, Spain, May 2020. IEEE.
Mark M. Wilde. Quantum Information Theory. Cambridge University Press, Cambridge, 2013.
Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chenggi Zhang, Connecting the dots: Multivariate time scries forecasting with graph neural networks. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 753-763, 2020.
F. Xu, G. Dong, J. Li, Q Yang, L. Wang, Y. Zhao, et al. Deep convolution generative ad advessaiial network-based electroencephalogram data augmentation for post-stroke rehabilitation with motor imagery. International Journal of Neural Systems, 32:2250039, 2022. doi: 10.11422/011206652255000
Lie Yang, Yonghao Song, Ke Ma, and Longhan Xie. Motor imagery E.G decoding method based on discriminative feature learning strategy. IEEE Transactions on Newral Systems and Rehabilitation Engineering, 29:368-377 2021. doi: 200100/
Lenka Zdeborová and Florent Krzakala. Statistical physics of inference: Thresholds and algorithms. Advances in Physics, 65:453, 2016.
Cha Zhang and Dinci Florencio. Analyzing the optimality of predictive transform coding using graph-based models. IEEE Signal Processing Letlers, 20(1):106-108 109, January 2013.
Dalin Zhang, Lina Yao, Kaixuan Chen, and Jessica Monaghan. A convolutional recurrent attention ion model for subject-indepenneent EEG signal analysis. IEEE Signal Processing Letters, 26:715-1 719, 2019. doi: N0.1109/LSP.20199 2906824.
H. Zhang, Н. Ji, Yu, J. Li, L. Jin, L. Liu, Z. Bai, and C. Ye. indeeentinde EEG classification based hybrid neural network. Frontiers in Neuroscience, 17:1124089, 2023. doi: 12389/ffnin..22 1124099.
15

---

Qiqi Zhang and Ying Liu. Improving brain computer interface performance by data augmentation with conditional deep convolutional generative adversarial networks. arXiv Preprint 1rXiv:1806.07100, 2018.
Zhi Zhang, Sheng hua Zhong, and Yan Liu. TorchEEGEMO: A deep learning toolbox towards EEG-based emotion recognition. Expert Systems with Applications, pp- 123550, 2024.
Daniel Zügner and Stephan Günnemann. Certifiable robustness and robust training for graph convolutional networks. In 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp- 246-256, Anchorage, AK, August 2019. Association for Computing Machinery.
16

---

A
Appendix
A.1 Hardware Specifications
All experiments were performed on an NVIDIA A100 Tensor Core GPU. The A100 features 6 912 CUDA cores, 3 456 FP64 CUDA cores, and 432 third-generation Tensor Cores, delivering up to 19.5 TFLOPS of FP32 throughput. It is equipped with 40 GB of HBM2e memory and offers a peak memory bandwidth of approximatcly 1.6 TB/s
A.2 Hyperparamter Selection and Model Details
Table 3: Hyperparameters for Financial Datasets
Model
Learning Rate
Epochs
Batch Size
Hidden Dim
Num Layers
Activation
Dropout
Betas
CAT
0.001
500
64
128
ELU
0.2
N/A
0.001
500
128
022
CDNN
0.001
500
128
ELU
0.2
Learned fiterbank, init: 100.01..00.00
GIN
0.001
500
61
128
1
ELU
0.2
N/A
Table 4: Hyperparameters for BC1 2A and PhysioNeet Datasets
Midil
Leurning Rate
Epocha
Barch Size
Hidden Dim
Num Lavurs
Activaton Functions
Dropout
Betax /rdde
0.00
500
Loueter al 1016
Landhem al 2016
Lanture
2015
Lachern 201
Lawbern
2015
N/A
VNN
50
R/A
CDNN
Tunb
I.I. 16
A.3 Data Availability
Zhang et al. 2024 with the Mother of all BCI Bench- marks (MOABB) (aris))Arstmmnhn et al., 2023 wrapper for the Motor Imagery data. The BCT 2A data is also available freely at https:/ www.bbbii.de/comppp.
We used the TorchEEG Python package
The Physionet dataset is available at nttps://physiooee..rrg sorrg/abbut/aae about/daaaaass//. The exchange rate dataset is from (Wu et al. 2020) and was obtained from the Github repository for (Cavallo et al. 2024h. The S&P 500 and US Stock market and
commoditics data are available freely from (Jander 2023] and ( KKmar 2024). For the US Stock datasct only
the price columns for each commodity was kept.
PhysioNet Dataset
This dataset comprises EEG recordings from 109 healthy volunteers. In each trial, visual target appears either the left or right side of a display. The participant then mentally rehearses opening and closing the corresponding hand until the target disappears, and subsequently relaxes. Data were acquired using the BCT2000 system with 64 EEG channels, sampled at 160 H6, and each trial has a duration of 3.1 s. Recordings from subjects 88, 89, 92, and 100 were excluded due to technical issues and excessive rest-state activity, resulting a final cohort of 105 participanis, each completing approximately 45 trials. We ran experiments on a 10 fold split of all individuals, where the test fold participants are unseen during training. All experiments were performed on the same randomization seed for all models for reproducible results. Raw EEG measurements without any preprocessing were used on this data set excpet what was already done in (Goldberger et al.) 2000).
'All Graph Based models treats temporal data as node features the underlying graph. For CAT and GIN the underlying graph fully-connected network. For CDNN and VNN the underlying graph is the sample covariance matrix and is computed over the complete training set and fixed during training. This matrix is then untouched and used for both training and testing. For VNN's the matrix is trace normaliser lized for fair comparison with CDNN's. For the BCT 2A and Physionet aatasets, graph convolution, concatenate every channel's full sequence of time points into single long feature vector which is fed into simple fully-comnected network (MLP), and apply Tanh nonlinearity followed by the classification layer. This is done for all graph models. For CDNN in all Datasets ignore the first un-fittredd term (k 0) reduce noise. architecture except
ELU activation for all models The Adam [Kingma Ba 2015) optimizer was used in all cases. For the financial data we use a 60/20/20 Train/Val/Ta Testt split and for all datasets evaluate on data using the model with lowest. mildation loss during training.
The financial datasets z-score normalized and involve the
17

---

BCI Competition IV 2a Dataset
This dataset contains EEG data from nine healthy participants performing a four-class motor imagery task (tongue, feet, right hand, left hand). Signals were recorded from 22 electrodes at 250 Hz. Each participant completed two sessions (training and testing) on separate days, with each session comprising 288 trials of 4s cach. Model evaluation employed leave-one-subjectoot vros-vvllidation for cach fold, all trials (both sessions) of one subject were held out as the test set, while the trials of the remaining eight subjects formed the training set.
The raw EEG recordings are first reduced to include only the channels of interest: all non-EEG channels (such as MEG and stimulus lines) are removed, leaving only the EEG signals. Next, the remaining signals are scaled by factor of 10 convert from volts to microvolts, ensuring numerical stability and consistency subsequent processing. A zero-phase Butterworth band-pass filter between 0.01 Hz and 38 Hz is then applied to eliminate slow drifts and frequerreeency artifacts. Following filtering, each channel is standardized over time using exponential moving estimate of its mean and variance, with an adaptation rate an initialization block of 1000 samples; this continuous normalization mitigates changes in signal amplitude and variance across the recording. Finally, fixed-length epochs are extracted around each event marker by applying a trial start offset of s before cue onset (with no post-cue offset), producing uniform windows ready for feature extraction and classification.
x 10-3 and
A.3.1 Learnable Filter Bank
We evaluate a filter bank of three learnable 3 values and observe no degradation in accuracy compared to hand picking 3 values. We note that choosing values beforchand based on domain-knowledgo can reduce training cost.
Table Accuracy (%) with learnable filter bank
Dataset
Accuracy
BCI 2A
49.48
PhysioNet
80.09
Avg
64.78
A.4
Ablation Analysis
A.4.1 EEG
We perform our ablation study on the most challenging dataset, the BCI-2A MI dataset. We repeat the exact. analysis over each independent participant using CDNNs with different values of 3 and finally with filter bank of I's. We also compare this to the same model without the covariance-density transform (i.c. VNN). To correct for chance agreement in this multi-class task we report Cohen's kappa score.
18

---

6.5
6.5
Tertkcipnte
Figure 4: Ablation over different configurations for each subject.
Kappa Scores Comparison
0.35
0.30
0.25
0.15
0.10
CDNN (0.1. 5.0. 15.11
0.00
VNN
CONN 10.11
C5NN15.0)
CONN (15.11
Figure 5: Mean Cohen's kappa for each ablation configuration.
We can clearly see that the covariance-density transform gives us a substantial improvement. Lower 3 values, in general, correspond to higher kappa (noise robustness), yet 8 = 5 outperforms 3 0.1 due to increased discriminability. Very large 3 remains noise-sensitive. Filter banks inherit both benefits, yielding the best
performance.
Financial Forecasting
We perform an ablation study on the US Stock dataset (horizon h = 3) to analyze the role of the diffusion parameter B.
19

---

MAE vs B (h=3)
Learnable B trajectories (h=3)
0.075
0.050
0.409
0.025
MAE (test)
0.408
B value
0.000
0.407
-0.025
-0.050
0.406
-0.075
best 8=5 MAE=0.405
00.00
2
3
100
200
300
400
500
B
Epoch
80
P1
82
B
Figure 6: Ablation on the US Stock Exchange Dataset
In the left panel of Fig. 6 we sweep over fixed values of 3 in a single scale CDNN and report the mean absolute error (MAF) on the US Stock dataset. Performance steadily improves as 3 increases, with the best score obtained at B This indicates that stronger diffusion improves performance in this dataset.
In the right panel, we initialize multi scale CDNN with four distinct learnable values close zero and plot their trajectories during training. Interestingly, the B values initially drift slightly negative before steadily increasing towards small positive values (zx 0.05 0.075 by epoch 500).
A.5 Scalability Analysis
Scalability (mean I std over seeds)
VNN
(cached)
20
CDNN (exp each step)
CDNN (cached exp)
ms step
10
0
500
1000
1500
2000
2500
3000
3500
4000
channels C
Figure 7: Runtime scalability of CDNN vs. VNN. Mean + standard deviation of per-step runtime (ms) as function of the number of channels (C). CDNNs scale linearly up to about C'=1024, after which the cost of repeatedly computing the matrix exponential cxp(-BL) grows rapidly. When caching the exponential, CDNN performance matches the near-linear scaling of VNNs, confirming that the additional cost can be effectively handled in practice.
We evaluated the computational scalability CDNNs compared to standard VNNs by measuring the average per-step runtime (forward + backward) across varying numbers of channels (C). Four configurations were tested: (1) VNN, (2) VNN with cached Laplacian, (3) CDNN computing exp(-3L) at each step, and
20

---

(4) CDNN using precomputed exponential (cached exp). Results are averaged over multiple seeds with standard deviation shown error bars.
Overall, CDNNs maintain efficient scaling behaviour for moderate graph sizes, and caching exp(-3L) restores runtime similarity with VNNs even for large-scale settings (C' 1000). Given that in all our experiments we pre-compute the covariance matrix exponential only once this shows that CDNN's are indeed fairly scalable.
A.6 Relationship between Covariance and Laplacian Eigenbasis
Sample Covatianee tiganeeuan (mm)
Tuee Leplacion Bigarestooss
Sumple Covatience digerveciors (omb03)
02
Blenvcctor
-06
Sample Comeriance Eiper vacies Ioe L0000)
Sample Covariance Egsnivecioss (0-200000)
Sample Covariance Eigevactirs (omSco0000)
Figure 8: Using signals generated from Erdös-Rényi Laplacian convolution show convergence of the eigenvectors of the sample covariance matrix to those of the underlying graph Laplacian
Let E R N a (combinatorial or normalized) graph Laplacian with eigen-decomposition
L L= U=UAU,
A= dia , AN).
(5)
where the columns of U are the Laplacian eigenvectors and 0 = A₁ < A2 5... S Ax are the Laplacian eigenvalues. Throughout, E|-] denotes statistical expectation and I is the identity matrix.
21

---

Following the graph-stationary framework of Dong et al. (Dong et al produced by filtering white noise ~ (O, I) through an weder-K polynomial graph fülter
2019 assumee the graph signal x is
au L w (() g(z)
Σ k=
ok 2
(6)
0=0
Because L and gL) commute, they are simultaneously diagonalizable by U:
C=x[xx]] [[L] = oLL) hL),,
(7)
where h(z) g())2. Equation equation 7 implies the covariance matrix C and Laplacian L share the same eigenvectors U: their eigenvalues are merely rescaled. (Dong et al.) 2019
If the graph arises as the precision matrix of GMRE, x ~ No L'): then the covariance is the (Moore Penrose) pseudo-inverss L Therefore
L L'=UA' At = 0, 0, ......"!))
so the eigenbasis is exactly preserved and the eigenvalues are reciprocals (up to the zero mode). Estimating L under Laplacian constraints is therefore equivalent to learning sparse inverse-covartance matrix. (Pavez & Orlega 2016) Conversely, if one has well-conditioned sample covariance, ils inverse supplies a valid Laplacian candidate.
More generally, any graph operator S that is a spectral function of L, i.c. S = f[L) for some analytic f, commutes with and shares ils eigenvectors:
SL LS = = S-U
Hence, whenever the true covariance can be expressed (or approximated) as f(L), the sample covariance obtained from sufficiently many realisations inherits the Laplacian eigenbasis. Results quantifying the perturbation of cigenvectors under finite-sample noise (e.g. matrix-Bernstein or Davis-Kahan bounds) then show that C converges Altafini 2021).
not only in Frobenius norm but also in the subspace spanned by U Foottan &
Propasition 4 (Shared eigenbasis of covariance and Laplacian). Let L be symmetric and C g(L) for some function with g(Ax)
U, and rank = > 2 0. Then C and L are simulianeously diagonalisable by the same orthonormal matris > 0H|. Morcover, if is monotone decreasing, the ordering of variances in C is the reverse of the Laplacian frequency ordering.
Proof. Since C = =[)) a polynomial (or analytic) function of L, commutes with L. Any two real symmetric matrices that commute are simultaneously diagonalisable by an orthogonal matrix; hence they
share U. The rank statement follows because 9(;;) 0 if and only if the corresponding cigenvalue of C
vanishes. If is decreasing, then A < >j = g(>x) > s(A)), reversing the order.
When the number of available graph-signal realisations n is much larger than the graph order N, the sample covariance C₂ itself becomes a spectral surrogate for the Laplacian L. Hence, one can substitute its leading eigenvectors for those of L and thereby circumvent the combinatorial optimisation required by sparse-Taplacian estimation procedures (Pavez & Ortega 2016). Moreover, because graph filters implemented in the covariance domain remain diagonal the shared cigenbasis U, frequency responses may be designed interchangeably on either operator, allowing practitioners work with whichever spectrum is numerically better conditioned casier estimate in given application. Finally, even when the inverse covariance C-1 fails satisfy the strict degree-constraintt. of Laplacian and is merely posilive-semiieelniiee its pseudo-inverse L still admits meaningful node-degree constraints.
interpretation that preserves the common eigenbasis while relaxing
Under broad conditions, stationary graph processes, GMRFs with Laplacian precision, or any model where the covariance is a spectral function of the eigenvectors of the sample covariance converge to those of the graph Laplacian. Consequently, algorithms that exploit the covariance matrix inherit graph-spectral interpretation, while graph filters may be implemented directly in the covariance eigenbasis.
22

---

Proof of Theorem 1
We aim to prove that the density filter frequency response h(pi) satisfies the Lipschitz condition:
(k(p(22) - h(p(A2)) S a0/2- - A1l
for constant a that depends on the filter parameters. The frequency response is defined as:
h&c -BC
Mp(A)))
Σ
Z
where:
-BA
e
i=1
Our focus is to bound the difference (h(p(A2)) h(p(Ai)) and derive meaningful constant. The difference between h(p()2)) and h(p(Ai) is:
hge Bhgk
hye -Ayk
[k(((A2) A(P(A)))
Zk
Zh
Factoring out the common terms, this becomes:
hx
(M((A22) - A(p(A)))
Σ
Z (e
-Br2k
--BAkk)
Using the Mean Value Theorem, know that for each k, there exists E (A1, x2) such that:
e -agk - e -ayk -Bk(A2 A1)e-
Substituting this into the expression, we have:
h& Z
1h(p(A2)) - A(p(Ai))
Σ
(-Bk(A>2 - A1)E -1)e--EA .BEK)
=
Using the Triangle inequality:
K
K
A=
==0
we can bound the summation:
h
- M((()))
<
(-Bk(>>2 - -2)e
z
km 0
Now, factor out LAg A11:
h BEk
-((())))-M((A)))) 1112 -111
7k
By definition of the partition function, z is at least as large e _ck ensuring that zt bound, we simplify:
-
< 1. Using this
Bke AEk
< 1he [|.3|]
Zh
23

---

Substituting this into the inequality, we get:
[k(((()2)) n(o(Ai))) A22 A1
hx |3k]]
Define the constant a as:
a |kx113k]
Thus, we have:
(h(p(A2)) -(((A))) a|/2- A1|-
Thus the covariance density filter Lipchitz Continuous with respect the denisty transformation, assuming constant filter coefficients.
While stability bounds be obtained for relative perturbations using Lipschitz filters more illuminating bounds can be obtained using integral Lipschitz filters. It can be shown that under certain conditions the covariance density filter satisfies being integral Lipschitz continuous.
The Integral Lipschitz condition imposes the additional requirement:
(h(A2) h(21)l 50. 112 A1l
2
To satisfy this condition, the constant a must be bounded by:
1
A21 2
sup
where the supremum is taken over all pairs of eigenvalues A₁ and Ag. Le the following condition must hold:
[Ai 2
sup
If this condition is satisfied, the filter satisfics the Composite Intogral Lipschitz condition with a constant of
1. Further observe The Integral Lipschitz condition can be generalized allow any real positive constant.
0 O. In this case, the condition becomes:
0
a
A
A 2
sup
This means that the filter satisfies the Integral Lipschitz condition with constant as long satisfies the above bound.
It follows that as long as is finite, the filter will satisfy the Integral Lipschitz condition for some choice of theta. The constant 0 can be tuned by adjusting the parameters [/y] (the filter coefficients) and B.
Proof of Theorem 2
We aim to prove that the filter
K
BSk
H(((S)) = hu (Tr(e 8 y
Σ
k=
24

---

is permutation equivariant in the vertex domain. Specifically, for permutation matrix T, want to show:
(
where $ = T ST is the permuted graph shift operator and = T" x is the permuted input signal.
This demonstrates that applying a permutation to the graph and the graph signal results in permutation of the filter output.
For the permuted graph, the shift operator is = T ST. Substituting this into the filter definition, have:
SSk
(((())
Σ
hy (Tr(e
km0
Using the Property of malrix exponentials under similarity transformations:
SS - - c
The trace of a matrix is invariant under cyclic permutations:
(Tr(e (i((e=$))* (Tr(TT AST))* (Ir(e ))..
Substituting these Properties into the filter definition:
T e-SKT H H(p(S)) Σ hp (Tr(e-85)) P5))*
Factor T and T outside the summation:
H(S) T Σ h (Tr(e-08))+
T T= H(ST.
Let the input signal and its permuted version & satisfy x = T x. Applying H($) to x:
2 = = H'
Using the orthogonality of P, where TTT I, this simplifies to:
2 = H())x
Thus, the filter output for the permuted graph and signal is the permuted version of the original filter output:
2 = T
where H(p(S))x.
Proof of Lemma 1
Proof. We decompose the error
3CC
-BC
-CC
-00 (1
-
[1-1
E
+e
2'
2.
25

---

Taking the operator norm and applying the triangle inequality yields
AC
SC
+ lle SC
Z'
Z'
Next, afler an application of Duhamel's formula (Kato 1966)
have
-- -- f (CC+CC) (8-5)) ds.
S(C+6C) - -
Changing the variable via s = 3t (so that ds B dl, t [0, 1]) yields
- -B S( -B(CC+)C) -BL-)( dt.
Taking the operator norm and using submultiplicatiiity gives
< |3 |8C|| //'' ||e (c -10 1|0-000-0))111000 dt. le 10
Define
1: 1/' lle Sc(C IIlle -lle- -(111))
We now consider the two cases:
(a) If 1 2 0: Since and + 8C are positive semidefinite,
l|e ( $ 1, a|l -(1-1) 1 At [0,1].
Hence,
I IS / ltt ldt=
(b) If 3 < 0: Write 8 = -13]. Then
||e -CC+++00 = = ||1 ( < [[
|e 11.11-1 S | |||||
and
Thus,
15 Jo
1311 |C+50|| (9((116) ||l|| di.
Since
B C + SC|| + |B |||| | + |8]t ||C -
- [[cI).
we factor out the constant term:
I < el 1l /' - [cl)) dt.
sC
Setting
d= (IC+
we have
f at
dt e
26

---

Thus, for B < 0
S |
[ 8CI -1|C[]]]
We now define the unified factor
3 0,
C,8)
elaliial
|01001++601111100 -1
3<
|8| || + -
To find the limit 3 - 0:
lim F(3, C, 8-0
For the Right-Hand Hand Limit(8 + 0+):
lim F(B, C, = 8
For the Left-Hand Limit (7))7), Let 8 where + 0f:
((-x,C,OC) - =
x[I 8C| |C|]
Let a + +C|| |||||
F(-x,C,60) e CII e
T0
Taking the limit as + 0+ a standard result gives us:
lim
C||
11
x-0+
za
lim
0-01
30
Therefore:
lim --0-
= =111 =x111
Both the right-hand and left-hand limits B + 0 are equal to 1. Hence,
lim F(B, C,8C) = 8
Continuing,
||e- (CC+CC) _ S0 < 3 F(3
-50
In particular, setting C = C 8C.
|e -BC BC -B0 | |CC| F(3, C,
27

---

Hence, the first term in error decomposition is bounded by
||e
SC
|3]
||| C,
Z'
Z'
Next, we bound the partition function difference. We have
1
1
Z - Z'
Z'
Z
Z
Since
Z' = "=TTe-AiCttoo = 1 TT -(C+&C) - -E-00).
and applying Tr(A) & m ||||| yiekds
((C+50)) - -000)
- S m (ec+-())))) - - s m 3 |C F(3,
Thus,
6 - 1 S |60|
We assume Z'= = R with R 1 and (this holds in high dimensional settings where the covariance matrix is not of full rank and in any this condition by regularizing the covariance matrix such that its smallest eigenvalue is 0), i.e the perturbed partition function is some multiple of the true partition function and thus it follows that
1
m |3
| (33, C,
m | |C F(8, C,OC)
Z'
Z
Z Z'
R
Thus, combining the two parts gives
3
||C]|| C,8)
+ lle 101 m 13 | F(3, C, 8C)
E
2'
Rt
Since Z' = R Z and Z 2 1 implics 1/7' S 1/RR we obtain
| ||0|| F(3, C, (1 (1+m ||e (|e
E|
<
R
Thus VB
| S m exp ( (1820) |8| [[c|]]
Finally, the error bound is given as
E < |||| (,,CC
(1+ m exp 1 (020) |8] (c1])
R
Lemma 2 (Perturbation Theory for the Density Matrix). Consider an ensemble density matriz and a sample density matria P. For any eigenvalue Ps > of P, the perturbation E, satisfies:
=v: = n 8v; + 8p. v; + (pp;1m - dv:.
where
75 = (pit - -P), 8v; = u, Vss 8p, w/ - Pi-
28

---

Proof. Note this proof follows analogously to the standard result in the perturbations of cigenvalues of the covariance matrix from (Sihag et al. 2022). As the density matrix shares the same eigenbasis as the covariance matrix the proof is straightforward with the only differences being in the transformed eigenvalues.
From the definition cigenvectors and eigenvalues, we have
P u; wy U4
We rewrite this terms of perturbations with respect to the ensemble density matrix and the outputs of eigendecomppoition as:
( - P)(v, + ov.) + P(v. + dv.) = (pr + (v: ov.),
where W/ Ps + 6ps and u, = V; + ov;
Using the fact that P v; = P: v, and rearranging terms, we have:
(P Plv: = (p.Im - Pjvv, + op.(v. ov.) - -(P (P- AP--jJv.
By setting = P-PP and 7/- = p.lm P, this simplifies to:
Ev; = n ov, + 6p: v; + (6p.I ov,
Proof of Theorem 3
We note that the density filters with respect to and P are given by
m h pk
m
H(P)
Σ K=
and H(P)
h
(1)
We aim to study the stability of the density filters by analyzing the difference between H(P) and H(P). For this purpose, next establish the first-order approximation for p in terms of P and E. Using P = P + E, the first-order approximation of p is given by
(P+E)*
p
P'EP
E,
(2)
-=0
where E satisfies
E
< Σ 2
(+)
E
P k-T Using (2), we have
H(P) - H(P) = Σ h (P + E)* - -P*]
(3)
k=0
m
h&
P'E D.
(4)
where D satisfies |ID||2 El|22) The focus of our subsequent analysis will be the first term in (4). For a random data sample x = [x1...., _m]" such that ||x]| < R, for some R > and x e gx its VFT with
respect given by VT x, where & [ār...., Fm T. The relationship between and x can be expressed
m
X
J
(5)
29

---

Multiplying both sides of (4) by x and leveraging (5), get
H(()
H(P)
Σ hk
Σ
P'EP x x DD.
(6)
k-1
31
11
P
Ev; + DN,
(7)
1=1
k=0
7=0
m
771
77
m
m
111
P Pi
7/
dv
P
Pi
0p:
V
31
11k
P Pj
(óp
E) dv.
(8)
Where we leveraged the result from Lemma that expands Ev,
Term 1:
This follows very similarly from in (39), we have
Sihag et al 2022
Using 7t = pilm and ōv, = u, - -v in the first term
m
m
k
p
(pilm - P) v.).
(9)
Finally, using P VA in (9), the first term in (8) is equivalent to:
hk
(A - (u;
(10)
1=1
k=0
7=0
i VL vi).
(11)
1=1
where L, diagonal matrix whose j-th element is given by:
m
k-
[. =
hi
Σ (p: - p,)p - r- pJ.
(12)
k=0
r=0
m
hv - pT) / 11
Σ k=
(13)
Pj
m
m
Σ hpp Σ hppp
(14)
k=
h(pp) h(p)).
(15)
After applying the operator norm:
hv
P
P
k
74 ōv:
[ max (h(p)) h(p;)) [[v] u,ll
0=
-=0
Here, vfu, u, the inner product between the eigenvector v, of the ensemble density matrix P and the eigenvector u, of the sample density matrix P. Using the result from
(Vershynin 2018 Theorem 4.1 and lou 2017 we conclude that if
sgu(A; 2x) > sgu(A} (x) Ad),
30

---

for X # Xj, the following condition holds:
m
m
m
1h(p:) h(1p)l
2hg
prr
pi
76 ov:
S v
Tr max J.FI
1> - Xy1 N1/2-8²
k=0
0=0
with probability at least for some E (0, 1/2], where
N
h; = V E XXT v. [2 P
From Theorem 1 we know that
(h(p2) h(11)l S sal22 A21.
Σ
1by 13k1
Note now that we can always choose (i.e by adjusting 3) such that:
(h(p)) h(s)
A Ayl-
hig
Substituting this bound into the above condition and applying a union bound, we have:
k-1 P
2 k k=u |
K
7/2 8v,
Σ |.
11/2-8
i=1
r=0
with probability at least
2km
N
N
where & is as defined in (Vershynin 2018 Furthermore, note that y m
Corollary 4.2 and Sihag el al. 2022
S m
|x|/2
If the random sample x satisfies ||x||2 S Q, then
k
2 2Σ k= I lk213kk|mQ
2km
P
pr
76 dv:
5
N1/2-8
N2c
N
i=1
k=
p=0
0
Term 2: We start with:
Σ
h4
P
P
8p: V,
=
hykp, 18p: V..
A=
Since +-1 < 1 for all Pr (as eigenvalues of P lie within the unit range of the covariance density matrix), we simplify:
172
hakp 8p. v
hulk&pp V..
Using the triangle inequality, relate km |hal|k the constant from Theorem 1, where:
1b ||||
<
k=0
31

---

Substituting this bound, we get:
m
m
m
hgkōpe V; 5
T Bellopelllvel
3
Since | 1, this simplifics further to:
3 h hgkppp
S
To Ts118ps|-
3
==
1=1
From Lemma 2, we have the bound
|
|0C||
C,6)) 1+ m exp (3200) |8] [C]]))
E
By Weyl's theorem, |ōpd
<>
||B|| (Gllb & Loan 2013), SO:
18p;1
[ || (,8)
[1 + ( 18 |[||))]
R
Substituting this bound into our inequality:
haköp, v;
|
(3,,C,
(1+ (1+m) exp{ [ 18| |8
R
Canceling the | terms, obtain
a||6C|| F(3, C,8C) (1+ m exp ((020) (8)[
he k 5p: =
Σ
R
k=0
1=1
Finally, noting that Em |2.1| I /m|xx/2 and if ||x|2 I Q, we have:
Σ
m
972
hu kōp V
a||8C|| F(3, C,8C)
(1 m exp ( 18]|C11)) vmQ.
R
k=U
Term 3: With the same argument regarding the invariance shifts in eigenbasis follows from [3] that:
I E I Ell
||óv:l|
o
with high probability.
Furthermore using the fact that for a random instance of random vector X whose probability distribution is supported within a bounded region w.l.o.g, such that ||x|| < 1, for some constant B > 0 and u > O, we have
log m + u
log m + a
P
|
ICIly
+ (1+ [[]]])
2 1 2
N
N
We expand out Lemma 1 as:
32

---

|31 <
[E|
1+ mexpl 40c0s |8| [||I]] || F(3, C, (,)
|3|
+ +1++11++1+ |C|I) log m + N
log m + u
F F(B,C,0C) R mxxpfnlnoac 18| [[[[]
B B|C
R
N
[3]
log m
u C(3,C,8C)
(1 11 [(<0)] |8| [[cI])
|
R
N
[8]
log m +
F(3, C,8)) 1 expp [sc0] |8| [[Ci]]
+
(1+
[C|]
R
N
Thus:
E =
with high probability.
VN
Which in turns implios:
I6 - EE) = N
Note that for positive 8 this always holds and for negative since from Lemma 1 F(3, C, tends to B lends to 0 we can always pick small to ignore the F term. Thus term 3 diminishes faster with N compared to terms and terms 2 and thus terms and 2 dominate the scaling behaviour of the overall upper bound.
The overall proof completed by noting that the condition on [|( - H[]]]x|| simplifies to the condition on the operator norm (P) H (P)] for any ||x]| < 1.
By unrolling this bound through L layers, and noting that each layer can at most amplify the perturbation by factor of F, obtain the overall network stability bound. The factor LFL-1 appears naturally from composing L layers each with at most F-fold channel combination where for each i E F there are potentially different a and 3.. The filter-level bound then carries through all layers.
Since the probability bounds and constants are inherited from the filter-level analysis, the final network-level bound follows directly, concluding the proof.
A.7 Sub-Additivity for Multi-Scale Von Neumann Entropy for Covariance Matrices
For this definition of entropy for covariance matrices to hold we would like to achieve the desirable sub- additivity Property of a valid entropy measure. The following theorem shows that this Property does indeed hold.
Theorem 5. Let C....., C E NxN be Individual Covariance Matrices. Define
Σ, Σ C
For each matria X, let
Zx = Tr(e- 9x], S(X) 3Tr[X + ln(Zx))
PX
Tr e 9x [
33

---

Then
S
S(Zn
En
Σ ==1
$ Σ S(C).. j=1
Lee. the Von-Neumann Entropy for Covariance Matrices satisfies the sub-addilicity Property.
Proof. The main chunk this proof the same as the one provided by Domenico et al. (Domenico Biamonte 2016) for the Laplacian- basod density matrix, however we use a regularization trick to ensure that the partition function of the covariance density matrix is always 1, ensuring non-negativity.
We will prove this by induction:
Base case (n = 2). We first prove the result for two covariance matrices C' and C?. Let us set ==:=²* C² We want to show
s(C' + C²) < S(CC)+ + S(C²)
We first regularize each matrix that min À 0. For 1,2, define
m
min À | A is an eigenvalue of C² (2 0 since C is p.s.d.).
Let
C = C m I.
Then mi) =
Observe:
S(C)) = S(C*),
because shifting by m² only multiplies e 00' matrix. Concretely,
e by factor et which cancels in the normalized density
e -
e+ -B0
e - = Bm e -00 = Pe. Tr(e-A(Ci1)
e+ Tr(e-9or = Pci S(Ò') = S(C)).
We define
Σn = C2
We rename:
A Č B:= C2 A==A+B.
All three A, C now have min (A) 0. Also
S(A) = S(C)), S(B) = S(CP),
S[A+B)- S(A + B) S(C C²).
Define
PA = c -AA
eBB
PA+B e -BA-BB)
/B
Z
7B
ZA+B
where Zx = T 9X and S(X) = T px| + In(2x).
Because A, B each has zero eigenvalue (or rank
deficiency), Za, Zn, ZA+B 2 1.
We now consider two KT. divergences:
D( ll PA) = (I
0,
D(pA+Bll 20.
Expanding each yields:
34

---

Hence
In = -3(A B) In ZA+B. InpA -3A - ZA.
So
I InAA - +B) A] In In + In = In ZA+B + Zln
Thus
D(PAABlPA)= -HAAAAAAABAALLBAAABB------ -T AA+B + +nZAa
Because =1,
D(pA+BllPA) = - -BTr(BAARB) - In ZA+B + + ZA.
Meanwhile S(A + B) (Tr((A - BA+BB] ZA+B. Observe
- pA -ln AA+B 8TH((A+B)) 3 AAA -In -(A+B)+ PA+B):
Hence
D(PA+BlPP) = + +) PA AIZA.
Since D(pA+BllPA) 2 we get
- S(A + S B + IZZ 2 0.
(1)
Likewise,
- B /ZB D(pA+BlPB) = ((MPA++--MP)))) 2 0.
-
One finds (by the same step):
D(PA+BIIP)) = B)A+B)+8THH) /A+B + +I 2 20.
(2)
We now have:
(L): S(A B) + BT. PA+B] + In AA 2 0,
(2): - S(A + +B) + Br PA+B] + In LB 2 0.
Further consider the fact that covariance matrices and their covariance density counterparts are always positive semi-definite. Therefore, the terms
(3): Tr PA 2
(4): Tr B PB] 2 0
are always non-negative.
Since at least eigenvalue of the covariance matrix is 0 due to the regularization trick (and thus their sum), the term
(5): In Zc
is also always non-negative because e 1, ensuring that the trace term is always least resulting in In(1) 0.
We now add the non-negative terms (1), (2), (3), (4), and observe that the incquality:
D(PA+BBIPB) + D(PA+BllP) + In + 8Tr[AApl + BT 20
We then exploit the fact that for Gibbs-like state the Von Neumann entropy is given by S(A + B) 3 8Tri(A- + B)PA+B] + In ZA-B- After expanding out the KL divergences, recalling that C = A + B, and re-arranging, we get:
22A B) + (A B)pA+B] In ZA+B + S(A) S(B) 0.
35

---

This, after some basic algebraic manipulation, allows us to conclude
B) I S(A) + S(B).
Thus S(A + BB S S(A) + S(B)
Finally, recall A C and B C2, so
S(Č* (2) <
But each Ć has the same entropy as C₁ and Č + C has the same entropy as C' + C². So
s(C² + C²) -=)=s((² 0 v
= S(C))+
This completes the base case 2 in all detail.
Inductive Step. Suppose for some 2 2, S(E)= C²) S S(C)). We show it for k + 1:
k+1
Σ C
Σ C*
c+++
j=1
j=1
By the n 2 sub-additivity (applying it to Σラ= C and C*+l),
k
Σ C +
c < S Eo + S(C)+)).
j=1
1=1
But by induction hypothesis, S SE k أ
C²)
$ Σ أー=1 S(C'). Therefore
k+1
k+1
s
<
S(C') S(C+t)(((+)))) = Σ S(C))
1=1
==1
==1
Hence by induction, for any n,
n S Σ S(C).
j=l
j=l
This completes the proof.
36

---

Original Covartance Matrix Cautsian Distribution
Denaity Hatric (Mi Ercropy)
Dursity Matrix (Max Entropy)
Entrope: 2.0380 B-15.00
Entrooy: 2.9830 (1
Entrocy vs. -Gussiin Dotriburter
044
Exponenba Distr bution
Entropy: 1.0541 et pt0000.00
Entropy: 2.9044 p-1.00
Entrosy vs D- Expccentall Distribution
26-
Origiaal Covaraanc Habrix
Density Matrix (Min Entropv)
Dersity Natrix (Max Entropy)
Distribution
15
9.0
11
Figure 9: We sample covariance matrices from Gaussian, Axuusiia,Exponenttia and Gamma distributions and observe the behaviour for different values of beta. Higher values of 8 relate to increased diffusson across the covariance matrix and thus reduced regularity and thus entropy, smaller values have higher entropy. We can thus conclude that CDNN's in higher entropy state are more stable, corresponding to traditional thermodynamical principles.
A.8 Multi-scale Von-Neumant Entropy for Singular Gaussians
Let X = (X,Y,Z)T N) with rank(E)) <=== 3. The classical differential entropy
h(X) = 1 + 10g(T))
+ log det Σ
(10)
2
2
is defined w.r.t. Lebesgue measure on R If det Σ the right term diverges to -o,, so every closed-form Gaussian entropy estimator explodes. We show that the multi-scale Von Neumann entropy (VNE) overcomes this obstruction and can distinguish singular covariances that differ only by global scale.
Proposition 6 (Density matrix "lifts" singular covariances). Let Σ 0 be x covariance matrir with rank > r < (hence det Σ = 0). For any finite 3 # 0 define the density matrix
Po(2) = exp(-3E)
tr trexp(-022
Then
ps(2)) > and det p2( > 0.
Proof. Diagonalise Σ = Q diag(( Ai.... dr, .). with Q orthogonal and X for S r. Ma-
trix exponential in the same basis is
= Q diag(e -881
ad e 82.
1QQ
Pr= and 1 (repeated p - times). All them are strictly positive, hence whose eigen-values are -A,, ,...., -Ar exp(-BE) > 0 and
detfexp(-B35)
-8 Er d 0.
37

---

Dividing by the positive scalar tr erexp(-932) Ei=1 aA preserves positive-defniteness and scales every eigen-value by the same constant c_t. Therefore
ps()) e -.
Ps(E) > 0,
> 0. ( eE,e=
Proposition guarantees that the Von-Neumann entropy Sg(E) is finite even when is rank- deficient, circumventing the divergence of the classical log det E term.
Take two different rank-2 covariances (both det 0):
2
0
0
0 0
Σ1
0
0
0
Σg
=
0
0
0
0
0
0
For their Von Neumann entropies are
(E = 1.28 bits,
Si(2) = 1.41 bits.
Both matrices are singular, hence the log det term diverges, but the VNE still produces finite, different values- for the "collapsed" Σ,, higher for the more evenly spread Σz. Thus VNE distinguishes degrees of randomness even among rank-deficient covariances where classical formulas fail.
Now consider commonly used naive estimator of the entropy of the covariance matrix.
p
X;
Suaive (Σ)
Σ
Ta loga T..
(11)
Σ
We begin by noting the following observation
Proposition 7 (Scale blindness of the trace normalised surrogate). Let E₁ 2 0 have rank < and Σg = with a > 0. Then Suaive (Σ2) Suaive (E₁).
Proof. Scaling each non-zero eigen-value by multiplies Tr by the same factor, leaving every ratio कं = ArTr unchanged; hence equation 11 is invariant.
38

---

Near-global scale (1.3, 1.2, 1.1)
1.4
1.2
1.0
S
0.8
0.6
0.4
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
Snaive
Figure 10: Scatter of naïve entropy Snaive versus von Neumann entropy S,322 for two covariance regimes Ej (blue) and Σg (orange) at near- glaball scale (eigenvalues multiplied by (1.3, 1.1,
At the extreme global-scale limit (e.g. uniform multiplication of all eigenvalues by a large factor), the naïve entropy remains exactly constant and is therefore blind to any change in overall variance magnitude. Even when the scale change is mild,as in Figure 10 where the spectrum shifts by element wise multiplication of (1.3, 1.2, bothh regimes overlap almost perfectly along the Snaive axis, indicating that normalized-spectrumm methods cannot distinguish them.
By contrast, von Neumann entropy at 3 2, remains sensitive absolute eigenvalue differences even when ratios are nearly equal. In the same near global-scale setting, S cleanly separates the two clusters into distinet vertical bands, capturing both the slight change in global variance and the subtle reshaping of the spectrum.
Quantitatively, a simple LDA classifier thresholded on Snaive yields an AUC of only 0.4992, essentially chance, whereas using S3=2 achieves an AUC of 0.9869, demonstrating that von Neumann entropy vastly outperforms the naïve measure al detecting even mild global-scale changes. In cases such as spike detection in neurological signals, the naive entropy estimator may miss out crucial near global changes that are common neural responses.
Of course, there may be cases where invariance lo scale may be desirable, thus CVNE acts complement to existing entropic measures, rather than replacement.
39

---

A.9 Reconstructing the True Covariance Matrix
Reconstructing the original covariance malrix involves formulating and solving optimization problem that connects the spectral Properties of the covariance matrix to a density-like representation.
Begin with a symmetric positive-definite covariance matrix that can be decomposed as
C UAU",
where U orthogonal and A diag(A1. ( An) contains the eigenvalues. The target distribution is obtained by normalising the eigenvalues,
A
Pj
i=
At the same time we define the covariance density distribution parameterised by an inverse temperature 3 I,
exp(-82.)
98 = E, exppp3B3)
Definition 8 (Moment objective). Let
n
f(3) = 3
> pid: i=1
log Σ e ( j=1
obtained by dropping the constant term E, Pi log from the Kullback Leibler divergence (p Proposition 8 (Strict convexity). The function f: IR + R of Deftnition 8 is twice continuously differentiable and strictly convex. Explicitly,
Ej= A BA, Σラ= e-
((0) =
Σ Pid:
(" Vargo [2] 0,
unless all eigenvalues are identical (a degenerate case).
Z3(B) 2(3)
Proof. Let Z(3) Eye-Ds, Then - (10)=E,ph,)
Because --0--,,>> -00-
E, ABa
(((3) = -m>
md
Σ pdd R [N].
e- R2g
Differentiating once more and applying the quotient rule gives f"( Eq E
[X12= = Var [A]. Because the cigenvalues are not all equal, this variance is strictly positive, so f is strictly convex and C².
Theorem 9 (Existence and uniqueness of the global minimiser). There exists unique value 3 E IR salisfying
E/=1 A-e-
f'(0)=
Σ Pid
Σn -8 Eqg- [A].
This 3 is the global minimiser of f(3) and therefore also mánimises DKL qo).
Proof. By Theorem 8 f is strictly convex, hence possesses at most one stationary point. Recall that ('(3) Xp - E.g. [A], where E, Because the p; are strictly positive and sum to one, X is a strict. convex combination of the eigenvalues; when these are not all equal this gives Amin < Ap < Amax- As 3 - -co the measure q8 concentrates on Amax: so lim f((3 = Xp Amax < O; as 3 - +oc concentrates on Amin: so lim f'(B) = X, Amin > Since f' is continuous and changes sign, the intermediate value theorem guarantecs the existence of I* with f(()*) 0. Strict convexity then makes this point unique and global.
40

---

With 8 obtained (e.g. numerically via Gradient Descent); the optimal covariance-density distribution is
2xp(-3*\,)
a Σ, (xp((33>,))
and the associated density matriz is reconstructed as
M exp(-3*0
tr(exp(-B*CC)
We generate covariance matrix from standard Gaussian matrix and BFGS routine to recover S*. Figures 1 and 12 display the original matrix, its reconstruction, and the alignment of the eigenvalue distribution. Consequently, Covarionce Density Network equipeed with the multi-scale filter bank { is generally capable of matching (and often surpassing) the representational power of classic Covariance Neural Network.
Hcatmaps: Trace-Nore Nore-Nommmm. Cov Reconstructed Frobenius difference 0.006347
Normalized Covariance Matrix
Reconstructed
0.30
0.00
0.06
0.06
9.04
0.04
0.00
Figure 11: Original and reconstructed covariance matrices.
(beta = -0.070)
Original distribution Reconstructed
0.15
0.14
0.13
Probability mass
0.12
0.11
0.10
0.08
Eigenvalue index
Figure 12: Eigenvalues of reconstructed and original covariance matrix.
41

---

Coverance Matto
Receitnnteed Darely Outneal
Tradl
60
200
Di gral Coverarce Marts
Saco Lowals
Oatinall Bitss 3.440
KDI Meles Level
Vaily dutt contann
Naley Covarance Narto
Rectrituced Darelty Prasic
322
Hox
The ceatriiionn
340
geal Owvrrare Marty
kree Level
Dutveel Bell 3555
801 by Livel
Nilly Croar Maari
Aecantturee Drmecly tratie
4.38
330
A20
Coverlarce Mark:
533
-ost
Volly cuttraatnn
4B%
Git
382
Covaranca Natro Lecel
Recenitaceed Damuity Mratis Ozoimal Sate 2.5005
Origial Covertecce Marto
to
AMM
0.5
Figure 13: Synthetic example: Gaussian data generate Ctrue; additive noise yiclds Cnoisy- Matching the noisy and true eigenspectra via 3 produces a density matrix that closely approximates Ctrue- See the connection lo Minimum Probability Flow Learning
et al. 2011).
We can also match realistic distributions by optimizing the beta parameter. For example, we can construct a Covariance Density Matrix with a noisy estimate of the true Covariance and optimize beta to match the noisy estimate the true one (See Figure 10).
While Figure 13 shows the when recreate covariance matrix from standard normal distrbution. We can also consider more challenging cases ie covariance matrices from auto-regressiv processes with non-trivial off-diagonal elements. A more intuitive approach could also be to reconstruct the precision matrix as this would lie in the positive 3 domain.
Figure shows the eigenvalue distribution reconstruction of the true precision and covariance matrix from highly correlated auto-regressive data. We note that while the actual entry by entry reconstruction of the covariance matrix was not ideal but the eigenvalue distribution was well recovered. We encourage future work in this domain.
Empirical Stability Analysis
A.9.1 Stability Analysis of the Operator Norm
The Operator norm difference between original and perturbed matrices is used to evaluate stability. particular we measure:
|||| = || Corigiaal Cperturbeill
and similarly for the density matrices:
|PP|
Poriginal
P perturbed
|]-
Where P perturbed is a perturbation only in the covariance matrix used to compute P.
42

---

Cov KDF.
D.1
Correlated KDF, Noise Opt B: 0.6409
Correlated Prac
0.1
Opt
Opt
40
60
30
40
20
20
10
10
0.00
0.01
0.02
0.03
0.04
-0.2
0.0
0.00
0.03
0.04
1.0
Std Normal Cov KDE, Noise 0.5
Std Normal Prec KDC, Noise 0.5
Correlated Cov KDE, Noise 0.5
Correlated Prec KDC, Noise 0.5
Opt B: -0.1111
Opt B: 0.1199
Ogt B: -0.2196
Opt P: 0.2417
50
30
40
20
10
10
0
0.00
0.06
-0.5
Std Normal KDE, Noise
Stdd
Correlated KDE, Noise Opt -0.1199
Correlated 0.3219
B: -0,0577
Opt 0.0639
40
10
40
30
30
20
2
0.00
0.02
0.04
0.06
0.08
-1
0
0.00
0.02
0.04
0.05
0.08
-0.5
0.0
0.5
1.0
1.5
Stdd Normal Cov KDE, Noise 1.5
Std Normal Prec KDE, Noise 1.5
Correlated Cov KDE, Noise 1.5
Correlated Proc KDE, Noise 1.5
Opt B: 0.0391
Opt B: 0. 0.046
Opt B: 0.0794
Opt P: 0.3259
50
30
15
20
10
30
0
-0.5
Figure 14: Reconstructing the eigenvalue distribution of a (Standard Normal vs highly correlated) Covariance and Precision matrix from a noise perturbed estimate using our moment matching routine.
43

---

tune Operator Norm Differences vs. Perturbation Strength
M w/Ww
opernear Roo uffeeere
0.5
3) Carcciance Mattia
Paturbation Shooggh
Figure 15: Stability of Covariance Density matrices to perturbations in the sample covariance matrix under the Operator norm.
We generate covariance matrices from Gaussian data and add noise perturbations at different levels. Observe that we pay large penaltics for negative üi however, this tends to decrease as 3 tends to Positive values of 3 show enhanced stability compared to the trace-normalized covariance matrix, although as increases, the stability decreases. This decrease is much less pronounced than in the negative scenario. Figure 10 thus corresponds well with Lemma 1.
A.9.2 Stability of CDNN in Synthetic Regression problem
After a Taylor Expansion P can be expressed as:
C
C
3C
P
21
31
(12)
2
As we purely want to observe the effect of C and its transformations the stability to the sample size of the covariance matrix, update as:
X Z
Xabifted = Px
(13)
where P applies the density matrix P to the input
As in Sihag et al. we compare the stability of CDNNs relative to perturbations in the sample covariance matrix, Le. we vary the number of samples used for the construction the covariance matrix. We replicate the exact conditions components, and PCA with Radial Basis Function (RBF) Kernel on random lincar regression problems using the routine sklearn. datasets make titee_rerresiirr in Python, which lets us specify various parameters. We generate cases, with no external noise and one with noise level of 5 (a parameter directly in the Python datasct generation).
Sihag al. and compare our approach with VNNs, Linear regression with PCA
Figure 6 shows the regression performance under no noise. We see that at smaller values of ß we maintain almost perfect stability; however, the MAE performance is weaker than VNNs. We conjecture that in the Friodmann regression problem (Breiman 1996), eigenmodes corresponding lower-variance principal components are the most discriminating, and since VNNs inherently discriminate these components better, they achieve strong performance. At. small values of / these low-variance components not discriminable
44

---

and thus the performance suffers, but 3 increases to larger values clearly see that more low- variance components are shifted to the discriminable eigenspace the performance improves, and at B = 15 we see improved performance compared to VNNs.
Comparison of Stability
Covariance Matrix Shift Operator
Covariance Density Matrix Shift Operator (beta-0.01)
Covariance Density Matrix as Shift Operator (beta=5.0)
Covariance Density Matrix as Shift Operator (beta=15.0)
PCA-LR
PCA-RBF
4.0
3.5
Mean Absolute
3.0
2.0
200
400
800
Samples
Computation
Figure 16: Regression performance under noise-free conditions.
Under noisy conditions, regression models based on density matrices (p) showed superior stability compared to covariance neural networks.
45

---

Comparison of Stability
Covariance Matrix as Shift Operator
Covariance Density Matrix as Shift Operator (beta=0,01)
Covariance Density Matrix Shift Operator (beta=5.0)
Density Shift Operator (beta=15.0)
PCA-LR
PCA-RBF
Absolute
20
200
400
800
Number Samples for Covariance Matrix Computation
Figure 17: Regression performance under noisy conditions.
Figure 7 repeats the experiment with the same conditions but the noise level is increased to 5. We can see that VNNs now suffer a significant drop in performance. This can be attributed to the noise being mostly concentrated in the low-variance components (i.e. the eigenspace that VNNs are best able to discriminate), significantly reducing performance. CDNNs discriminate best in the high-variance cigenspace, so the noise is less likely to affect performance. Thus, even at small B we see stable and stronger performance than VNNs. PCA also operates in this high-variance space and thus does not suffer from noise as much. This suggests that CDNNs inherently exhibit greater robustness to external noise, regardless of the value of positive B.
46


Reviews:
[{"review_id": "B0N3ontlru", "reviewer": "Reviewer_Z2dP", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "**Summary**\nThis paper introduces CDNNs. They use sample covariance matrix of the data as a central operator and transform the covariance matrix into a density matrix and use it to define filtering operations. A parameter controls how much the model trades stability for discriminability. The method is evaluated on financial forecasting and EEG classification tasks. The results show improvements over covariance neural networks and several graph-based baselines.\n\nOverall, the paper is technically strong and the experimental results are promising. However, I found the presentation difficult to follow without a strong background in graph signal processing.\n\n**Strengths**\n\n- Novel and original idea by combining covariance-based models with density-matrix constructions.\n\n- Deep theoretical analysis that shows care in studying stability and robustness.\n\n- Experimental results on finance and EEG tasks are strong and suggest practical relevance.\n\n- The method improves over prior covariance neural networks and shows good robustness to noise.\n\n**Weaknesses and Suggestions**\n\n- Introduction and Background (Graph Signal Processing). The paper heavily relies on concepts from graph signal processing (GSP), such as graph shift operators, spectral filtering, and transferability. However, the current introduction assumes substantial prior knowledge. For readers not specialized in GSP, it is hard to build intuition. I strongly recommend adding a more pedagogical introduction or an intuitive subsection that explains the main ideas of GSP and covariance-based filtering in simpler terms, possibly with a small running example.\n\n- In the experimental section, it is not always clear: what exactly is the input to each model, its dimensionality,  what is the output of the CDNN layers, how temporal data is transformed before classification. Also, how the covariance matrix is computed and reused at training and test time. A simple diagram/example could be used to improve readability.\n\n- Computational Complexity and Practical Cost While runtime is reported, I miss a clearer discussion of: computational complexity with respect to number of nodes / channels, memory cost of the density matrix and matrix exponential, how these costs compare to non-graph baselines.\n\n- Baselines Beyond Graph-Based Models Most comparisons are against graph-based or covariance-based models.  I recommend including additional strong non-graph baselines, such as: TABPFN or similar tabular foundation models (where applicable),  transformer-style baselines with comparable parameter counts, other non-graph time-series models.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The main claims of the paper are supported by theoretical analysis and experimental results", "audience_yn": "Yes", "audience_explanation": "Yes. Researchers familiar with graph signal processing and machine learning for graphs would likely be interested in this work.", "requested_changes": "**Weaknesses and Suggestions**\n\n- Introduction and Background (Graph Signal Processing). The paper heavily relies on concepts from graph signal processing (GSP), such as graph shift operators, spectral filtering, and transferability. However, the current introduction assumes substantial prior knowledge. For readers not specialized in GSP, it is hard to build intuition. I strongly recommend adding a more pedagogical introduction or an intuitive subsection that explains the main ideas of GSP and covariance-based filtering in simpler terms, possibly with a small running example.\n\n- In the experimental section, it is not always clear: what exactly is the input to each model, its dimensionality,  what is the output of the CDNN layers, how temporal data is transformed before classification. Also, how the covariance matrix is computed and reused at training and test time. A simple diagram/example could be used to improve readability.\n\n- Computational Complexity and Practical Cost While runtime is reported, I miss a clearer discussion of: computational complexity with respect to number of nodes / channels, memory cost of the density matrix and matrix exponential, how these costs compare to non-graph baselines.\n\n- Baselines Beyond Graph-Based Models Most comparisons are against graph-based or covariance-based models.  I recommend including additional strong non-graph baselines, such as: TABPFN or similar tabular foundation models (where applicable),  transformer-style baselines with comparable parameter counts, other non-graph time-series models.", "broader_impact_concerns": "Nothing at this point", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "mH4a6BOHim", "reviewer": "Reviewer_wk2E", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Summary and Overall Assessment\n\nThe paper introduces Covariance Density Neural Networks (CDNNs), proposing to construct a density matrix from a covariance matrix and use it as a graph operator in a GNN-like architecture. The authors present some theoretical analyses and empirical studies in financial forecasting and EEG classification.\n\nMy evaluation is largely negative, which could be a result of my unfamiliarity of the field. That being said, the paper lacks a clear motivation and fails to articulate what problem it actually aims to solve. For a reader not already familiar with the graph-signal-processing lineage of GNNs, it is extremely unclear what CDNNs are designed for, why this formulation is needed, and what concrete advantage it offers over prior covariance-based or diffusion-kernel methods. The exposition assumes substantial background knowledge without providing a proper introduction or positioning within the existing literature.", "claims_and_evidence_yn": "No", "claims_explanation": "Main Concerns\n\n1. Poorly defined problem and insufficient motivation.\nThe introduction does not identify a specific gap or limitation in existing GNN or covariance-based frameworks that CDNN meaningfully resolves. The description of the model reads as a reparameterization of the CoVariance Neural Network (VNN) with a temperature-like parameter β, but the reader is left guessing what benefit this confers in practice.\n\n2. Lack of clarity and accessibility.\nThe background material is underdeveloped. Without familiarity with prior work on VNNs or graph diffusion kernels, the technical construction appears unmoored—there is no intuitive explanation of when one would want to use this operator or what kind of data it is intended for. A substantially expanded literature review and clearer motivation are needed to make the contribution intelligible to a general ML audience.\n\n3. Triviality of theoretical analysis.\nTheoretical results (equivariance, Lipschitz and stability bounds) appear to be routine adaptations of standard GNN analyses rather than novel insights. They mainly confirm expected behavior of smooth spectral filters and do not establish new principles.\n\n4. Empirical results difficult to interpret.\nThe experiments show some improvements, but the paper does not convincingly link those gains to the proposed method. As I am not a specialist in GNNs, I cannot judge how compelling the empirical evidence is within that subfield. I am willing to defend my assessment of weak motivation and conceptual triviality, but I acknowledge that I might be missing the paper’s intended main point.\n\nRecommendation\n\nThe work would benefit from a much clearer articulation of the target problem, a grounded comparison to existing covariance- and diffusion-based graph methods, and a deeper explanation of why the proposed construction is necessary. In its current form, the contribution and novelty remain ambiguous, and the theoretical analysis adds little substance beyond standard results.\n\nAgain, although I'm willing to defend my assesment, I have to admit the paper falls beyond my scope of expertise.", "audience_yn": "Yes", "audience_explanation": "The paper is relevant to some audience obviously. Still, it should clarify the background and relevant literatures substantially to attract people who was not familiar with graph-signal-processing and GNNs.", "requested_changes": "Please see my comments above.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "JLOLopzE4n", "reviewer": "Reviewer_HrxG", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Standard GNNs require a pre-specified graph to define connections between nodes; in many cases, this graph is unknown and difficult to design without strong prior knowledge. This paper proposes a covariance-driven alternative by constructing the graph operator directly from data. Its core contribution is to replace the covariance shift operator with a covariance density operator and use (optionally multi-scale) graph filtering to mix node features. The authors provide grounded theoretical development and show improvements over baselines on financial forecasting and EEG classification tasks.", "claims_and_evidence_yn": "Yes", "claims_explanation": "Claims are well supported by theoretical development and empirical results. The authors are also honest about the limitations and explicitly acknowledged that current work is limited to linear dependency and stationary data.", "audience_yn": "Yes", "audience_explanation": "This works provides a plug-and-play module for constructing graph operators directly from data,  which is a timely and useful contribution for the community.", "requested_changes": "**Q1** I feel confused about how to calculate the covariance matrix $C$ for data from different domains (e.g time-histories data has a time dimension along with the feature dimension, while other type of data may only have feature dimension) . What constitutes one \"observation\"? (timepoints or trials for EEG data? )   \n\n**Q2** The term “quasi-Hamiltonian” is best understood as a spectral analogy rather than a literal physical construct. I recommend adding an explicit clarifying sentence in Section 2 stating that the terminology is metaphorical and is intended to highlight the role of (C) in shaping the spectrum, not to claim a strict physical interpretation.\n\n**Q3**  typos  page 2. . “how it alllows us to control” ->  allows \n\npage 8 “ added covariance information is indeed informative and that that that”, remove extra “that” \n\nProof of Theorem 1 : “respect to the denisty transformation” -> density  \n\n**Q4** citation style, please use standard TMLR citation format (authors, year)", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)