Sparse Covariance Neural Networks

## Abstract

Abstract Covariance Neural Networks (VNNs) perform graph convolutions on the covariance matrix of tabular data and achieve success in a variety of applications. However, the empirical covariance matrix on which the VNNs operate may contain many spurious correlations, making VNNs’ performance inconsistent due to these noisy estimates
and decreasing their computational efficiency.
To tackle this issue, we put forth Sparse coVariance Neural Networks (S-VNNs), a framework that applies sparsification techniques on the sample covariance matrix before convolution. When the true covariance matrix is sparse, we propose hard and soft thresholding to improve covariance estimation and reduce computational cost. Instead, when the true covariance is dense, we propose stochastic sparsification where data correlations are dropped in probability according to principled strategies. We show that S-VNNs are more stable than nominal VNNs as well as sparse principal component analysis. By analyzing the impact of sparsification on their behavior, we provide novel connections between S-VNN stability and data distribution.
We support our theoretical findings with experimental results on various application scenarios, ranging from brain data to human action recognition, and
show an improved task performance, stability, and computational efficiency of S-VNNs compared with nominal VNNs.

## 1 Introduction

Covariance-based data processing is key to machine learning pipelines due to its ability to whiten data distributions, identify principal directions, and estimate the interdependencies among features. Such advantages have been shown in several applications including brain connectivity estimation (Bessadok et al., ; Qiao et al., ), financial data (de Miranda Cardoso et al., ; Wang & Aste, ) and human action recognition through motion sensors measurements (Liao et al., ; Wang et al., ).
One prominent approach is Principal Component Analysis (PCA), which maximizes the variance of data points by projecting them onto the eigenvectors of their covariance matrix (Jolliffe & Cadima, ).
However, PCA-based approaches are unstable to errors in the covariance estimation, i.e., a poor estimate of the covariance matrix and its eigenvectors may lead to unpredictably bad results, especially in high-dimensional and low-data regimes or when the covariance matrix eigenvalues are close to each other (Jolliffe, ; Paul, ; Baik et al., ).
To overcome this issue, coVariance Neural Networks (VNNs) have been proposed (Sihag et al., ). VNNs apply graph convolutions to the covariance matrix, an operation that, similarly to PCA, learns the importance of the principal directions of the data. Due to their graph convolution operation, they inherit the stability property of GNNs (Gama et al., ) and extend it to errors in covariance estimation, which guarantees more reliable performance even in the presence of poor covariance estimates since their output difference when operating on the sample and true covariance is bounded (Sihag et al., ; ; ).
Nevertheless, VNNs operate on the sample covariance matrix, which results in two limitations. First, VNNs are sensitive to settings in which covariance estimation is problematic, especially in high-dimensional and low-data regimes (Paul, ; Baik et al., ; Féral & Péché, ). Second, VNNs are computationally expensive and memory inefficient as the sample covariance matrix is typically dense no matter whether the true covariance matrix is dense or sparse due to sample estimation errors,
which is especially problematic on high-dimensional datasets.

To overcome these limitations, we propose sparsification-based covariance regularizers for VNNs to improve the covariance matrix estimation.
We perform stability analysis for sparse VNNs (S-VNNs) to investigate the effects of covariance sparsification techniques on the performance of VNNs.
When true covariance matrices are sparse, we propose hard and soft thresholding strategies and show that S-VNNs are more robust to both the finite-data estimation error and the covariance sparsification thresholding.
For generic scenarios where true covariance matrices are not necessarily sparse, we propose a stochastic sparsification approach and prove the stability of S-VNNs in this setting. Our findings quantify a trade-off between stability and sparsification (i.e., efficiency). That is, dropping covariances according to their absolute values leads to higher stability but, in some cases, low sparsification, whereas dropping a fixed percentage of covariance values improves sparsification at the expense of stability. We corroborate our findings with experiments on synthetic and real datasets highlighting the benefits of sparsity on covariance-based neural processing both in terms of performance and computation. Our specific contributions are as follows.

- (C1)
Sparse covariance neural networks. We propose covariance sparsification techniques for VNNs to improve performance and reduce computation. When the true covariance is sparse, we perform hard and soft covariance thresholding. When the true covariance is dense, we propose a stochastic sparsification framework that allows controlling the sparsification level, and design specific strategies in this framework based on application requirements.
- (C2)
Stability to covariance sparsification. We characterize the effects of covariance sparsification techniques on the VNN performance through stability analysis. We show that S-VNNs can not only reduce the computation but also improve the stability under appropriate settings, compared to the nominal VNN and the sparse PCA.
- (C3)
Empirical validation. We validate our theoretical findings with experiments on real and synthetic datasets, demonstrating the importance of sparsification in low-data settings and the improved efficiency of S-VNNs on brain data and human action recognition use cases.

## 2 Problem Formulation

Consider $t$ samples $\{\mathbf{x}_{i}\}_{i=1}^{t}$ of a random vector $\mathbf{x}\in\mathbb{R}^{N}$ with mean $\bm{\mu}=\mathbb{E}[\mathbf{x}]$, covariance $\mathbf{C}=\mathbb{E}[(\mathbf{x}-\bm{\mu})(\mathbf{x}-\bm{\mu})^{\mathsf{T}}]$ and respective estimates $\bm{\hat{\mu}}=1/t\sum_{i=1}^{t}\mathbf{x}_{i}$, $\mathbf{\hat{C}}=1/t\sum_{i=1}^{t}(\mathbf{x}_{i}-\bm{\hat{\mu}})(\mathbf{x}_{
i}-\bm{\hat{\mu}})^{\mathsf{T}}$. The estimated covariance $\mathbf{\hat{C}}$ captures the data interdependencies and can be used
for principal component analysis (PCA) via its eigendecomposition $\mathbf{\hat{C}}=\mathbf{\hat{V}}\mathbf{\hat{\Lambda}}\mathbf{\hat{V}}^{
\mathsf{T}}$ with eigenvectors $\mathbf{\hat{V}}=[\mathbf{\hat{v}}_{0},\dots,\mathbf{\hat{v}}_{N-1}]\in\mathbb
{R}^{N\times N}$ and eigenvalues $\mathbf{\hat{\Lambda}}=\text{diag}(\hat{\lambda}_{0},\dots,\hat{\lambda}_{N-1})$, to maximize the variance of the data points by projecting them on its eigenvectors, i.e., $\mathbf{\tilde{x}}=\mathbf{\hat{V}}^{\mathsf{T}}\mathbf{x}$.
However, this finite-sample PCA is unstable as the estimated eigenvectors $\mathbf{\hat{V}}$ may differ significantly from the true directions $\mathbf{V}$ especially when the true eigenvalues in $\mathbf{\Lambda}$ are close to each other (Jolliffe, ).

To overcome this issue, Sihag et al. () proposed graph convolutional covariance filters as an extension of PCA that is stable to finite-sample estimation uncertainties even with close eigenvalues.
Each data point is seen as a node of a graph and the data sample $\mathbf{x}$ as a graph signal, where the covariance $\mathbf{C}$ represents the graph structure. Covariance filters build on the standard graph convolutional filters, which shift a signal $\mathbf{x}$ over a graph represented by the shift operator $\mathbf{S}$ (e.g., adjacency matrix or Laplacian) according to the equation $\sum_{k=0}^{K}h_{k}\mathbf{S}^{k}\mathbf{x}$, where $K$ is the order of the filter and $h_{k}$ are learnable parameters (Gama et al., ; Ortega et al., ; Isufi et al., ). Several popular GNNs, such as GCN (Kipf & Welling, ), are special instances of the generic graph convolutional filter.
Building on this, the covariance filter of order $K$ and the VNN layer, which assembles a filter bank of $F_{\text{in}}\times F_{\text{out}}$ covariance filters and a non-linear activation function $\sigma(\cdot)$, are

$$ $\displaystyle\mathbf{u}=\sum_{k=0}^{K}h_{k}\mathbf{\hat{C}}^{k}\mathbf{x}= \mathbf{H}(\mathbf{\hat{C}})\mathbf{x}~{}~{}\textnormal{and}~{}~{}\mathbf{u}^{ l}_{f}=\sigma\left(\sum_{g=1}^{F_{\text{in}}}\mathbf{H}^{l}(\mathbf{\hat{C}}) \mathbf{u}^{l-1}_{g}\right)~{}f=1,\ldots,F_{\text{out}},~{}l=1,\ldots,L.$ (1) $$

This operation is related to PCA as it processes the eigenvectors of the covariance matrix (see Appendix A for details).
We denote the VNN architecture as $\mathsf{\Phi}(\mathbf{x},\mathbf{\hat{C}},\mathcal{H})$, where $\mathcal{H}=\{h_{klgf}\}$ contains all parameters.
The output of the last layer $\mathbf{u}^{L}=\Phi(\mathbf{x},\mathbf{\hat{C}},\mathcal{H})$ contains the final representations generated by VNN and can be directly used for a downstream task (e.g., classification or regression) or further processed by a readout layer.
The model parameters $\mathcal{H}$ are optimized to minimize a task-specific loss (e.g., cross-entropy for classification or mean squared error for regression tasks) over a training set.
Inheriting the stability to perturbations of neural network structures (Gama et al., ; Bruna & Mallat, ), VNNs have been shown stable to errors in covariance estimation, as stated in the following theorem.

###### Theorem 1 (Stability of VNNs (Sihag et al., 2022 ) ) .

Consider a VNN $\mathsf{\Phi}(\mathbf{x},\mathbf{C},\mathcal{H})$ of $L$ layers with $F$ covariance filters $\mathbf{H}(\mathbf{C})$ per layer and nonlinearities $\sigma(\cdot)$ such that $|\sigma(a)-\sigma(b)|\leq|a-b|$. Let each filter’s frequency response be Lipschitz with constant $P$(^1^11For a definition of Lipschitz filter see Def. 1.). Consider a generic data sample $\mathbf{x}\sim\mathcal{N}(\mathbf{0},\mathbf{C})$ such that $\|\mathbf{x}\|\leq 1$. Then, with probability $1-\mathcal{O}(1/t)$ it holds that

$$ $\begin{gathered}\|\mathbf{H}(\mathbf{\hat{C}})\!-\!\mathbf{H}(\mathbf{C})\|\! \leq\!\frac{P}{\sqrt{t}}\mathcal{O}\left(\!\!\sqrt{N}\!+\!\frac{\|\mathbf{C}\| \sqrt{\log(Nt)}}{\nu t}\!\right)\!\!=\!\beta\!~{}\textnormal{and}~{}\|\mathsf{ \Phi}(\mathbf{x},\mathbf{C},\mathcal{H})\!-\!\mathsf{\Phi}(\mathbf{x},\mathbf{ \hat{C}},\mathcal{H})\|\!\leq\!LF^{L\!-\!1}\!\beta\end{gathered}$ $$

for some $\nu>0$ capturing the data distribution, and $\|\cdot\|$ denoting the $2-$ or the spectral norm.

This notion of stability follows an extensive line of research on stability of GNNs to generic graph perturbations (Gama et al., ; Cerviño et al., ; Arghal et al., ; Parada-Mayorga & Ribeiro, ; Keriven et al., ; Ruiz et al., ), but it specifically considers covariance estimation errors, which is fundamental for VNNs. Indeed, since in practical applications the true covariance is not available and VNNs operate on an estimate of it, this notion of stability acts as a certificate of performance guarantee w.r.t. the ideal scenario when working with the true covariance matrix and characterizes the impact of data and model characteristics on the performance.
More in detail,
the output difference of VNN caused by the sample estimation error is bounded proportionally by the square root of the number of samples $t$. When $t$ increases, the bound approaches zero, ultimately, showing convergence to the true data distribution.
This result is specific to covariance estimation errors, whereas generic GNN stability bounds do not decrease with the number of samples.
This improves the stability w.r.t. PCA, which we provide in the following lemma.

###### Lemma 1 .

Consider a true covariance matrix $\mathbf{C}$ and a sample covariance estimate $\mathbf{\hat{C}}$ with respective eigendecompositions $\mathbf{C}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{\mathsf{T}}$ and $\mathbf{\hat{C}}=\mathbf{\hat{V}}\mathbf{\hat{\Lambda}}\mathbf{\hat{V}}^{
\mathsf{T}}$.
Then, for any signal $\mathbf{x}$ with $\|\mathbf{x}\|\leq 1$, it holds with high probability that
$\|\mathbf{V}^{\mathsf{T}}\mathbf{x}-\mathbf{\hat{V}}^{\mathsf{T}}\mathbf{x}\|
\leq\mathcal{O}(t^{-1/2}(\min_{i,j\neq i}|\lambda_{i}-\lambda_{j}|)^{-1})$.

Lemma 1 shows that PCA stability is inversely proportional to the smallest gap between covariance eigenvalues, leading to unstable behaviors when the eigenvalues are close. VNNs do not suffer from this as the covariance filter can exhibit a stable response to close eigenvalues at the expense of lower discriminability (Sihag et al., ).
Because of this advantage in finite-data settings, VNNs have been shown effective in covariance-based learning tasks, both on static and temporal data (Sihag et al., ; ; Cavallo et al., ; ). However, they have two major limitations. First, when the true covariance matrix is sparse, the finite-sample estimate contains spurious correlations (Baik et al., ; Paul, ), affecting the performance of the VNN significantly. Second, when the true covariance is dense, a VNN on the finite-sample estimate is limited by its high quadratic computational complexity in the data dimension, restricting its applicability to low-dimensional settings only.

Paper objective. Our goal is to overcome the above limitations by sparsifying the sample covariance matrix for the VNN. While PCA-based processing sparsification has been studied (Bickel & Levina, ; Deshpande & Montanari, ), its extension to the VNN is challenging because it is unclear: (i) how to perform sparsification depending on whether the true covariance matrix is sparse or dense; (ii) what are the effects of sparsification on VNN stability. When the true covariance matrix is sparse, we apply
thresholding-based sparsification strategies for the sample covariance matrix, and prove that they result in stability improvements
(Thms. 2-3). When the true covariance matrix is dense, we put forth a stochastic sparsification framework in a form akin to dropout (Def. 4), characterize its impact on the VNN stability (Thm. 4), and propose principled sparsification strategies based on these findings (Sec. 4.3). All proofs are collected in the appendix.

###### Theorem 1 (Stability of VNNs (Sihag et al., 2022 ) ) .

###### Lemma 1 .

## 3 Sparse True Covariance

In applications involving brain and spectroscopic imaging, weather forecasting, and finance, the sample covariance matrix is notoriously prone to spurious correlations due to limited sample size (Bickel & Levina, ; Jobson & Korkie, ; Ledoit & Wolf, ).
Sparsity-based regularizers for the covariance estimation have been proposed as a backup (Bickel & Levina, ). In this section, we consider the scenario where true convariance matrices are sparse, propose two sparsification strategies, hard and soft thresholding, and analyze their effects on the VNN embeddings.
The choice of hard and soft thresholding instead of other regularized covariance estimations (Ledoit & Wolf, ; Bickel & Levina, ; Bien & Tibshirani, ; Friedman et al., ) lies in the fact that thresholding yields a fast approach and allows theoretical tractability.

Instrumental tool for our analysis will be the frequency response of a covariance filter. By computing the graph/covariance Fourier transform of input and output signal of the covariance filter (cf. equation 1) we get
$\mathbf{\tilde{u}}=\mathbf{\hat{V}}^{\mathsf{T}}\mathbf{u}=\mathbf{\hat{V}}^{
\mathsf{T}}\sum_{k=0}^{K}h_{k}[\mathbf{\hat{V}}\mathbf{\hat{\Lambda}}\mathbf{
\hat{V}}^{\mathsf{T}}]^{k}\mathbf{x}=\sum_{k=0}^{K}h_{k}\mathbf{\hat{\Lambda}}
^{k}\mathbf{\hat{V}}^{\mathsf{T}}\mathbf{x}$, which, for the $i$-th entry, leads to
$\mathbf{\tilde{u}}_{i}=\sum_{k=0}^{K}h_{k}\hat{\lambda}_{i}^{k}\mathbf{\tilde{
x}}_{i}=h(\hat{\lambda}_{i})\mathbf{\tilde{x}}_{i}$.
This indicates that the frequency response of the covariance filter is a polynomial $h(\lambda)=\sum_{k=0}^{K}h_{k}\lambda^{k}$ with frequency variable $\lambda$ specified on the eigenvalues $\hat{\lambda}_{i}$.

###### Definition 1 (Lipschitz covariance filter) .

The covariance filter is Lipschitz with constant $P$ if, for every pair of eigenvalues $\lambda_{i},\lambda_{j}$, $i\neq j$, its frequency response satisfies: $|h(\lambda_{i})-h(\lambda_{j})|\leq P|\lambda_{i}-\lambda_{j}|$.

Lipschitz covariance filters control the variability of the frequency response through the constant $P$. Higher $P$ allows the frequency response to generate different outputs at close eigenvalues and improves the filter’s discriminability, but may degrade its stability as shown in the following.

###### Definition 1 (Lipschitz covariance filter) .

### 3.1 Hard thresholding

Hard thresholding removes the covariance entries below a given value and it has been shown to improve covariance estimation in high-dimensional low-data settings (Bickel & Levina, ).

###### Definition 2 (Hard thresholding) .

Given the sample covariance matrix $\mathbf{\hat{C}}$ and a coefficient $\tau>0$, the hard thresholding function is
$\eta(\mathbf{\hat{C}})_{ij}=\hat{c}_{ij}$ if $|\hat{c}_{ij}|\geq\tau/\sqrt{t}$, 0 otherwise.

The hard threshold is inversely dependent on the number of samples as $1/\sqrt{t}$. This follows the intuition that the non-thresholded estimator approaches the true covariance as $t$ increases, hence, the sparsification is less needed and disappears in the limit of $t\rightarrow\infty$.
Hard thresholding provides a more reliable covariance estimate because it removes small spurious finite-sample errors and, thus, improves the performance of VNNs when the true covariance is sparse. We analyze the impact of hard thresholding on the VNN stability in the following theorem.

###### Theorem 2 .

Let the true covariance $\mathbf{C}$ belong to the sparse class $\mathcal{C}=\{\mathbf{C}:c_{ii}\leq M,\sum_{j=1}^{N}\mathbbm{1}[c_{ij}\neq 0]
\leq c_{0},\forall i\}$ where $M>0$ is a constant, $\mathbbm{1}(\cdot)$ is the indicator function and $c_{0}$ is the maximum number of non-zero elements in each row of $\mathbf{C}$. Consider a hard-thresholded sample covariance matrix $\mathbf{\bar{C}}$ following Def. 2 with $\tau=M^{\prime}\sqrt{\log N}$ and $M^{\prime}$ large enough. With a high probability, it holds that

$$ $\displaystyle\|\mathbf{H}(\mathbf{\bar{C}})\mathbf{x}-\mathbf{H}(\mathbf{C}) \mathbf{x}\|\leq t^{-1/2}Pc_{0}\sqrt{N\log{N}}(1+\sqrt{2N})+\mathcal{O}\left(t ^{-1}\right).$ $$

Thm. 2 shows that hard-thresholded covariance filters –and consequently VNN from the results in Thm. 1– are stable to covariance estimation errors and that they
converge to the respective filter and VNN operating with the true covariance as $1/\sqrt{t}$ in the number of data samples. Importantly, hard thresholding improves the stability with a tighter bound than dense VNNs (Sihag et al., ) and sparse PCA (Bickel & Levina, ).

Comparison with dense VNN. Contrasting the results of Thm. 2 and Thm. 1, we see that both bounds decrease with the same order $\mathcal{O}(t^{-1/2})$. However, the hard-thresholded VNN bound depends linearly on the number of non-zero elements in a row $c_{0}$, whereas that of the dense VNN depends on the spectral norm $\|\mathbf{C}\|$. For sparse covariance, $c_{0}\ll\|\mathbf{C}\|$ and the hard thresholding improves stability.

Comparison with sparse PCA.
To provide further insights on the stability of sparse
covariance filters, we contrast it also with the stability of sparse PCA.

###### Proposition 1 .

Consider a true covariance matrix $\mathbf{C}$ and the thresholded sample covariance estimate $\mathbf{\bar{C}}$ with respective eigendecompositions $\mathbf{C}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{\mathsf{T}}$ and $\mathbf{\bar{C}}=\mathbf{\bar{V}}\mathbf{\bar{\Lambda}}\mathbf{\bar{V}}^{
\mathsf{T}}$. Then, for any signal $\mathbf{x}$ with $\|\mathbf{x}\|\leq 1$, it holds with high probability that

$$ $\displaystyle\|\mathbf{V}^{\mathsf{T}}\mathbf{x}-\mathbf{\bar{V}}^{\mathsf{T}} \mathbf{x}\|\leq t^{-1/2}(\min_{i}|\lambda_{i}-\lambda_{i+1}|)^{-1}c_{0}N\sqrt {2\log{N}}$ (2) $$

That is, the stability of sparse PCA is inversely proportional to the minimum eigenvalue gap. This term is not present in the stability of S-VNN, i.e., Thm. 2, as it is absorbed by the filter Lipschitz constant $P$: indeed, the covariance filter can exhibit a stable frequency behavior for close eigenvalues at the expense of lower discriminability, but the latter is compensated in the subsequent VNN layers that use cascades of filterbanks and non-linearities. Thus, hard-thresholded VNNs attain improved stability w.r.t. sparse PCA in sparse covariance setting.

###### Definition 2 (Hard thresholding) .

###### Theorem 2 .

###### Proposition 1 .

### 3.2 Soft thresholding

When the data follows the spiked covariance model, such as in electrocardiogram or brain image data (Johnstone & Lu, ), soft thresholding has been studied as a better alternative to achieve more reliable covariance estimates (Deshpande & Montanari, ).
In this model, data points follow $\mathbf{x}_{i}=\sum_{q=1}^{r}\sqrt{\beta_{q}}u_{q,i}\mathbf{v}_{q}+\mathbf{z}_
{i}$,
where $\mathbf{v}_{1},\dots,\mathbf{v}_{r}\in\mathbb{R}^{N}$ are orthonormal vectors with exactly $c_{0}$ non-zero entries of magnitudes lower-bounded by $\theta/\sqrt{c_{0}}$ for some constant $\theta>0$, $u_{q,i}\sim\mathcal{N}(0,1)$ and $\mathbf{z}_{i}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ are i.i.d., and $\beta_{q}\in\mathbb{R}_{+}$ is a measure of the signal-to-noise ratio.
In this case, we sparsify the covariance estimate with soft thresholding.

###### Definition 3 (Soft thresholding) .

Given the sample covariance matrix $\mathbf{\hat{C}}$ and a coefficient $\tau>0$, we define the soft thresholding function as
$\eta(\mathbf{\hat{C}})_{ij}=\hat{c}_{ij}-\textnormal{sign}(\hat{c}_{ij})\tau/
\sqrt{t}$ if $|\hat{c}_{ij}|>\tau/\sqrt{t}$, 0 otherwise.

Unlike hard thresholding which only removes small noisy values, soft thresholding subtracts a value from all entries to remove noise also from non-zeroed coefficients.
Again, we set the threshold to decrease with the number of samples analogously to (Deshpande & Montanari, ) as more accurate covariance estimates reduce the need for sparsification.
The following theorem links soft thresholding to the VNN stability.

###### Theorem 3 .

Consider a soft-thresholded estimate of the covariance matrix $\mathbf{\bar{C}}$ as per Def. 3 with $\tau=M^{\prime}\sqrt{\log(N/c_{0}^{2})}$ and $M^{\prime},C$ two large enough constants.
Let the eigenvalues of the true covariance $\{\lambda_{i}\}_{i=0}^{N-1}$ be all distinct. Then, the following holds with high probability:

$$ $\displaystyle\|\mathbf{H}(\mathbf{\bar{C}})\mathbf{x}\!-\!\mathbf{H}(\mathbf{C })\mathbf{x}\|\!\leq\!t^{-1/2}P\sqrt{N}Cc_{0}\max(1,\!\lambda_{\textnormal{max }})\sqrt{\max\left(\log(N/c_{0}^{2}),1\right)}(1+\sqrt{2N})\!+\!\mathcal{O} \left(t^{-1}\!\right)\!.$ $$

Thm. 3 shows that VNNs are stable also when the covariance estimate is soft thresholded and the stability bound decreases with the number of available samples as $t^{-1/2}$.
The main takeaways for hard thresholding w.r.t. dense VNNs and sparse PCA hold also for soft-thresholded VNNs.
Furthermore, soft thresholding provides better stability than hard thresholding in low-data settings under spiked covariance, as the bound in Thm. 3 contains $\sqrt{\log(N/c_{0}^{2})}$ which is smaller than $\sqrt{\log N}$ in Thm. 2.

Computational complexity.
A sparse VNN layer has a computational complexity of order $\mathcal{O}(\|\mathbf{\bar{C}}\|_{0}KF_{\text{in}}F_{\text{out}})$, where $\|\mathbf{\bar{C}}\|_{0}$ is the number of non-zero values of the thresholded sample covariance. This is significantly better than a VNN with dense covariance matrix, which is of order $\mathcal{O}(N^{2}KF_{\text{in}}F_{\text{out}})$, since $\|\mathbf{\bar{C}}\|_{0}$ is generally much smaller than $N^{2}$.

###### Definition 3 (Soft thresholding) .

###### Theorem 3 .

## 4 Generic True Covariance

Often, the true covariance matrix and its estimate are both dense which makes the VNN computationally heavy. While we can use hard or soft thresholding or any graph sparsification approach, the resulting sparse VNN may be significantly different from the VNN operating on the true covariance, ultimately, leading to divergent outputs and degraded performance.
To reduce the computational cost in a tractable manner, we propose a stochastic sparsification framework in a form akin to dropout. Such an approach is general and does not have any structure assumption on the covariance matrix.

###### Definition 4 (Stochastic sparsification) .

Let $\mathbf{\Delta}$ be a matrix of the same support as the sample covariance $\mathbf{\hat{C}}$ with entries $\delta_{ij}=\delta_{ji}=1$ with probability $p_{ij}$ and 0 otherwise (i.e., $\textnormal{Bernoulli}(p_{ij})$) and $\delta_{ii}=1$.
A sparsified covariance matrix is $\mathbf{\tilde{C}}=\mathbf{\Delta}\odot\mathbf{\hat{C}}$, where $\odot$ is the elementwise product.

Stochastic sparsification generates random sparsified matrices $\mathbf{\tilde{C}}$ of the sample covariance matrix via element-wise independent sampling in a way that $\mathbf{\tilde{C}}$ preserves: (i) the symmetric property in consistency with the covariance principle; and (ii) the variance of the data points on the main diagonal. The sparsification can be controlled by the probabilities $p_{ij}$ depending on application requirements.

Figure: Figure 1: Two-dimensional generalized frequency response $h(\bm{\lambda})$ of stochastic sparsified covariance filter. $h(\bm{\lambda})$ is determined by the coefficients $\{h_{k}\}_{k=0}^{K}$ and is independent of the covariance realizations. For specific covariance realizations, $h(\bm{\lambda})$ is instantiated on the corresponding multivariate frequencies (e.g., $\lambda_{11},\lambda_{12}$ and $\lambda_{21},\lambda_{22}$ for two different realizations $\mathbf{\tilde{C}}_{1}\mathbf{\tilde{C}}_{2}$ and $\mathbf{\tilde{C}}^{\prime}_{1}\mathbf{\tilde{C}}^{\prime}_{2}$).
Refer to caption: x1.png

###### Definition 4 (Stochastic sparsification) .

### 4.1 VNNs with stochastic covariance sparsification

We now investigate the effects of stochastic sparsification on VNNs’ stability by framing the stochastic sparsification as a stochastic perturbation problem following (Gao et al., ) and defining stochastic sparsified covariance filters.

###### Definition 5 (Stochastic covariance filter) .

Given a sequence of i.i.d. realizations $\mathbf{\tilde{C}}_{k},\dots,\mathbf{\tilde{C}}_{1}$ of the stochastic sparsified covariance $\mathbf{\tilde{C}}$ independent on data distribution, a covariance filter $\mathbf{H}(\mathbf{\tilde{C}})$ performs convolution of a generic signal $\mathbf{x}$ as $\mathbf{\tilde{u}}=\mathbf{H}(\mathbf{\tilde{C}})\mathbf{x}=\sum_{k=0}^{K}h_{k
}\mathbf{\tilde{C}}_{k}\dots\mathbf{\tilde{C}}_{1}\mathbf{\tilde{C}}_{0}
\mathbf{x}$ with $\mathbf{\tilde{C}}_{0}=\mathbf{I}$.

Stochastic filters shift the signal $\mathbf{x}$ over $K$ different random realizations $\{\mathbf{\tilde{C}}_{k}\}_{k=1}^{K}$ of $\mathbf{\tilde{C}}$ rather than shifting over a fixed $\mathbf{\hat{C}}$.
This leads to a generalized frequency interpretation.

###### Definition 6 (Generalized covariance filter frequency response) .

The generalized frequency response of the stochastic covariance filter is a multivariate function of the form
$h(\bm{\lambda})=\sum_{k=0}^{K}h_{k}\prod_{\kappa=0}^{k}\lambda_{\kappa}$
where $\bm{\lambda}=[\lambda_{1},\dots,\lambda_{K}]^{\mathsf{T}}$ with $\lambda_{0}=1$ is a generic vector variable and each frequency variable $\lambda_{k}$ corresponds to the covariance realization $\mathbf{\tilde{C}}_{k}$ for $k=1,\ldots,K$.

The derivations of such a frequency response are reported in Appendix E. Fig. 1 illustrates an example of $K=2$. This extends the deterministic case, in which all realizations and corresponding eigenvalues are equal with $\prod_{\kappa=0}^{k}\lambda_{\kappa}=\lambda^{k}$.
Given the generalized frequency response, we generalize the integral Lipschitz property for stability analysis.

###### Definition 7 (Lipschitz gradient) .

Consider the analytic generalized frequency response $h(\bm{\lambda})$ and two instantiations $\bm{\lambda}_{1}=[\lambda_{11},\dots,\lambda_{1K}]^{\mathsf{T}}$ and $\bm{\lambda}_{2}=[\lambda_{21},\dots,\lambda_{2K}]^{\mathsf{T}}$ of vector variable $\bm{\lambda}$. Consider also the auxiliary vector that concatenates the first $k$ entries of $\bm{\lambda}_{2}$ and the last $K-k$ entries of $\bm{\lambda}_{1}$, i.e., $\bm{\lambda}^{(k)}\!=\![\lambda_{21},\dots,\lambda_{2k},\lambda_{1(k+1)},\dots
,\lambda_{1K}]^{\mathsf{T}}$. The Lipschitz gradient of $h(\bm{\lambda})$ between $\bm{\lambda}_{1}$ and $\bm{\lambda}_{2}$ is

$$ $\displaystyle\nabla_{L}h(\bm{\lambda}_{1},\bm{\lambda}_{2})=\left[\partial h( \bm{\lambda}^{(1)})/\partial\lambda_{1},\dots,\partial h(\bm{\lambda}^{(K)})/ \partial\lambda_{K}\right]^{\mathsf{T}}$ $$

where $\partial h(\bm{\lambda}^{(k)})/\partial\lambda_{k}$ is the partial derivative w.r.t. $\lambda_{k}$ at $\bm{\lambda}^{(k)}$.

###### Definition 8 (Generalized integral Lipschitz filter) .

A covariance filter is generalized integral Lipschitz if there exists a constant $P$ s.t. $\|\nabla_{L}h(\bm{\lambda}_{1},\bm{\lambda}_{2})\|_{2}\leq P$ and $\|\bm{\lambda}_{1}\odot\nabla_{L}h(\bm{\lambda}_{1},\bm{\lambda}_{2})\|_{2}
\leq P.$

The Lipschitz gradient characterizes the variability of $h(\bm{\lambda})$ since, for two multivariate frequency vectors $\bm{\lambda}_{1},\bm{\lambda}_{2}$, we have $h(\bm{\lambda}_{2})-h(\bm{\lambda}_{1})=\nabla_{L}^{\mathsf{T}}h(\bm{\lambda}_
{1},\bm{\lambda}_{2})(\bm{\lambda}_{2}-\bm{\lambda}_{1})$.
The generalized integral Lipschitz filter limits this variability to be at most linear in the multidimensional space and to decrease as the frequency $\bm{\lambda}$ is specified at large values, which extends the standard integral Lipschitz property to the multivariate frequency domain.

###### Definition 5 (Stochastic covariance filter) .

###### Definition 6 (Generalized covariance filter frequency response) .

###### Definition 7 (Lipschitz gradient) .

###### Definition 8 (Generalized integral Lipschitz filter) .

### 4.2 Stability of VNNs with stochastic covariance sparsification

These preliminaries allow us to analyze the stability of VNNs to stochastic covariance sparsification.

###### Theorem 4 .

Consider a randomly sparsified covariance filter $\mathbf{H}(\mathbf{\tilde{C}})$ [cf. Def. 5] that is generalized integral Lipschitz with constant P [cf. Def. 8]. Let also $\mathbf{H}(\mathbf{C})$ denote the same filter (fixed parameters) operating on the true covariance matrix. Then, for a generic signal $\|\mathbf{x}\|\leq 1$, the expected squared difference between the two filters can be upper-bounded as

$$ $\mathbb{E}[\|\mathbf{H}(\mathbf{C})\mathbf{x}-\mathbf{H}(\mathbf{\tilde{C}}) \mathbf{x}\|^{2}]\leq\underbrace{NP^{2}Q+\mathcal{O}((1-p_{1})(1-p_{2}))}_{ \textnormal{sparsification error}}+\underbrace{\frac{P^{2}}{t}\mathcal{O}\left (N+\frac{\|\mathbf{C}\|^{2}\log(Nt)}{\nu^{2}t^{2}}\right)}_{\textnormal{ covariance uncertainty}}$ $$

where $Q=\sum_{i=1}^{N}\sum_{n=1}^{N}\hat{c}_{in}^{2}(1-p_{in})$, $p_{1},p_{2}$ are generic probabilities such that $\mathcal{O}((1-p_{1})(1-p_{2}))$ is dominated by the linear terms in the probability value in $Q$ and $\nu$ is defined in Thm. 1.

Thm. 4 identifies two main factors that affect the stability of stochastic sparsified covariance filters: the covariance uncertainty and the sparsification error.

Covariance uncertainty.
This term is analogous to the dense case in Thm. 1. It decreases with the number of samples as $\mathcal{O}(1/t)$ and is dominated by the sparsification error for sufficiently large
$t$. Note that
a large $t$ does not necessarily correspond to a small covariance estimation error here because the stochastic sparsification may still remove large covariance values; ultimately, making our analysis principally different from (Sihag et al., ) which assumes small covariance perturbations.

Sparsification error.
This term decreases as the sampling probabilities $p_{ij}\to 1$, corresponding to an improved stability but a lower sparsification. This indicates a trade-off between the perturbation effect caused by the covariance sparsification error and the computational cost saved by the stochastic sparsification. The stability constant depends on the filter Lipschitz constant $P$, the data dimension $N$, and the coupling between the covariance values $\hat{c}_{ij}$ and the sampling probabilities $p_{ij}$, i.e., $Q$. A larger $P$ allows for a higher filter discriminability as the frequency response can change more quickly but leads to worse stability. Data with high dimensionality $N$ results in a larger graph and increases the effect of sparsification on the VNN stability.
More importantly, $Q$ represents the interplay between sample covariance values $\hat{c}_{ij}$ and their corresponding sampling probabilities $p_{ij}$, which can be used as a design choice to develop different stochastic sparsification strategies in a principled manner as we shall discuss in the next section. Finally, note that the results in Thm. 4 and its proof generalize and differ substantially from the results in (Gao et al., , Theorem 1) as we consider here edge-specific probabilities $p_{ij}$ rather than an identical probability $p_{ij}=p$ and we identify a connection between VNN stability and the data distribution through their covariances $\hat{c}_{ij}$, which is not present in (Gao et al., ).

###### Theorem 4 .

### 4.3 Stochastic sparsification strategies

Following our theoretical analysis, we propose two strategies to assign the sparsification probabilities to covariance values and assess their impact on the stability of VNN.

Absolute covariance values (ACV). Despite being dense, the covariance matrix may contain some values of lower magnitude that may be reminiscent of spurious correlations. We want therefore to drop these small covariance values with a high probability to save computational cost while keeping large covariance values to maintain the useful dependencies. Thus, we define the probability as $p_{ij}=|\hat{c}_{ij}|/\hat{c}_{\text{max}}$, where $\hat{c}_{\text{max}}=\max_{i,j}|\hat{c}_{ij}|$. Fig. 2 (left) shows that the stability term $q_{in}$ within $Q$ [cf. Thm. 4] is small when the covariance value $\hat{c}_{in}$ approaches zero or the maximal value $\hat{c}_{\max}$. This is because small covariance values have little impact on the VNN stability even if they are more likely to be dropped, and large covariance values (close to $\hat{c}_{\text{max}}$) are less likely to be dropped though their removal would affect stability more.
It peaks around $\hat{c}_{in}=2/3\hat{c}_{\text{max}}$, which identifies an interval in which ACV may have a notable impact on the stability.

Ranked covariance values (RCV). While ACV improves efficiency by preserving stronger correlations, it does not allow any control over the amount of sparsification (e.g., if all covariance values are high, very few are dropped). To overcome this, we define a set of probability values with a desired mean $p$ and assign them to covariance values based on their positions in this absolute ranking. Formally, consider an ordered set of probabilities $\mathcal{P}=\{p^{\prime}_{1},\dots,p^{\prime}_{N^{\prime}}\}$, where $p^{\prime}_{i}\leq p^{\prime}_{i+1}$, and each $p^{\prime}_{i}$ is sampled from $\mathcal{N}(p,\sigma)$ with $\sigma=\min((1-p)/3,p/3))$ (such that the number of values not in $[0,1]$, which we clip to the interval, is negligible) and $N^{\prime}$ is the number of probability values for assignment. We set $p_{ij}=p^{\prime}_{k}$ where $k=|\{c_{lm}:|c_{lm}|<|c_{ij}|;l,m=1,\dots,N\}|$ is the position of $c_{ij}$ in the ordered ranking of absolute covariances.
Consequently, the expected percentage of dropped covariances is $1-p$, which allows controlling the sparsification level at the risk of removing useful covariances.
In practice, the value of $p$ can be treated as a hyperparameter and tuned through cross-validation or on a validation set based on the desired accuracy and stability.

Fig. 2 (right) shows how the term $Q$ changes for ACV and RCV with different means of covariance value distribution. RCV with a smaller $p$ corresponds to lower stability especially when the covariance values are high, because they are dropped regardless of their value. ACV, instead, maintains a consistent level of stability since it balances the covariance magnitudes and the dropping probabilities.

Computational complexity.
The expected computational complexity of a VNN layer with stochastic sparsification is of the order $\mathcal{O}(\mathbb{E}[\|\mathbf{\tilde{C}}\|_{0}]KF_{\text{in}}F_{\text{out}})$, where $\mathbb{E}[\|\mathbf{\tilde{C}}\|_{0}]$ is the expected number of non-zero elements of the sparsified covariance. For ACV, we have that $\mathbb{E}[\|\mathbf{\tilde{C}}\|_{0}]=\hat{c}_{\text{mean}}(\|\mathbf{\hat{C}
}\|_{0}-N)+N$, where $\hat{c}_{\text{mean}}=1/N^{2}\sum_{i,j=1,\dots,N,i\neq j}|\hat{c}_{ij}|$, so the complexity depends on the data. Instead, for the RCV we have $\mathbb{E}[\|\mathbf{\tilde{C}}\|_{0}]=p(\|\mathbf{\hat{C}}\|_{0}-N)+N$, so the complexity can be controlled through $p$.

###### Remark 1 .

The stochastic sparsification with RCV generalizes the hard thresholding approach in Sec. 3, which can be recovered by setting probabilities as $p_{ij}=1$ if $|c_{ij}|>\tau/\sqrt{t}$, 0 0 otherwise.

###### Remark 2 .

Both thresholding and stochastic sparsification provide symmetric covariance estimates but do not necessarily preserve positive semidefiniteness (PSD), which is common in related works (Bickel & Levina, ; Deshpande & Montanari, ; Li et al., ; Liao et al., ) and places our techniques within the scope of the literature. Bickel & Levina () also provide a sufficient condition for PSD under hard-thresholding, which we report in Appendix G.

Figure: Figure 2: Insights into stability for stochastic sparsification. (Left) Stability bound elements $q_{in}=\hat{c}_{in}^{2}(1-|\hat{c}_{in}|/\hat{c}_{\text{max}})$ with $Q=\sum_{i=1}^{N}\sum_{n=1}^{N}q_{in}$ for ACV [cf. Thm. 4]. (Right) Stability term $Q$ [cf. Thm. 4] for covariance values distributions with different means and stochastic sparsifications.
Refer to caption: x2.png

###### Remark 1 .

###### Remark 2 .

## 5 Numerical Results

We corroborate the proposed methods with experiments on real and synthetic datasets, targeting the following objectives:
validate our theoretical results on
(O1) thresholding when the true covariance matrix is sparse and
(O2) stochastic sparsification for a generic true covariance;
(O3) compare S-VNNs and VNNs on real datasets.
We provide details and additional results in Appendix I.

### 5.1 Stability of sparse VNN

Experimental setup – (O1)-(O2). We generate synthetic data with a controlled covariance by sampling data points $\mathbf{x}_{i}\sim\mathcal{N}(\mathbf{0},\mathbf{C})$.
Then, we create regression targets as $y_{i}=\mathbf{w}^{\mathsf{T}}\mathbf{x}_{i}+u$, where $\mathbf{w}$ is a vector with elements $w_{j}\sim\text{Uniform}(0,1)$ and $u\sim\mathcal{N}(0,3)$ is a noise term. We generate 1000 samples of size $N=100$ and divide the data into train/validation/test splits of size 80%/10%/10%.

Sparse true covariance – (O1).
We evaluate the stability of S-VNN on a synthetic dataset with sparse covariance (SparseCov). We train a VNN with the true covariance and test it with different covariance estimates obtained with a varying number of samples. We compare the original sample covariance, hard thresholding (Hthr), and soft thresholding (Sthr). We use PCA-SVM as baseline, i.e., we transform data with PCA followed by an RBF-SVM for regression.
Fig. 3 (left) reports the MAE w.r.t. the number of samples for covariance estimation.
Overall, hard and soft thresholding provide better covariance estimates, ultimately, improving the performance (i.e., lower MAE) and increasing the stability compared to the nominal VNN. VNNs are more stable than PCA-based models, because their performance is less affected by covariance perturbations as analyzed in (Sihag et al., ). These results support our theoretical discussion in Sec. 3, i.e., thresholding allows VNNs to maintain stability in sparse covariance settings.

Generic true covariance – (O2).
We assess the stochastic sparsification using synthetic datasets with large covariance values (LargeCov) and with small covariance values (the synthetic linear regression dataset in (Sihag et al., , Appendix E.4) with tail=0.2, here called SmallCov); see Sec. I.1 for the distributions of the covariance values. We train a VNN with the true covariance matrix and test it with the stochastically sparsified versions.
Fig. 3 (center) and (right) compare the ACV and RCV with different values of $p$. For LargeCov, we see that a high sparsification leads to greater instability (RCV for small $p$) as large covariance values are dropped which hinders performance. For SmallCov, instead, all sparsification approaches affect stability only lightly as removing small covariance values does not lead to performance changes. ACV remains stable because it sparsifies little even when the covariance values are large. In turn, this corroborates our theoretical observations in Sec. 4.1.

Figure: Figure 3: Stability of VNN. (Left) MAE with hard and soft thresholding on SparseCov (standard deviations are of order at most $10^{-2}$ for VNNs and $10^{-1}$ for PCA). (Center) MAE for stochastic sparsification on LargeCov. (Right) MAE for stochastic sparsification on SmallCov.
Refer to caption: x4.png

Figure: Figure 4: Accuracy and time for a forward pass for sparse and dense VNNs. From left to right, the results are for Epilepsy, CNI, MHEALTH and Realdisp. RCV results are for $p=0.75,0.5,0.25$ from top to bottom, and for Hthr and Sthr we report results for 3 different thresholds to achieve analogous sparsification. Legend is shared.
Refer to caption: x7.png

### 5.2 Sparse VNN on real data

Experimental setup – (O3).
We consider four real datasets of brain measurements and human action recognition data.
Epilepsy (Kramer et al., ) consists of an electrocardiography time series equally split before and after an epileptic seizure event. We perform binary classification to predict whether a time sample was recorded before or after the seizure.
CNI (Schirmer et al., ) contains resting-state fMRI time series of patients with attention deficit hyperactivity disorder (ADHD) and neurotypical controls (NC), which we use as labels for binary classification.
MHEALTH (Banos et al., ) and Realdisp (Banos et al., ) contain measurements of wearable devices placed on subjects performing different actions. The goal is to classify the action performed by each subject. Additional details are reported in Sec. I.1.
We compare various sparsification approaches (thresholding, RCV, ACV) against the dense VNN. For completeness, we report additional results for various PCA-based classifiers and GraphSAGE in Sec. I.3, although our objective is not to achieve state-of-the-art performance on these datasets but to reduce the impact of spurious correlations while improving computational efficiency of VNNs.
For Hthr and Sthr, we set the threshold to achieve a level of sparsification comparable to RCV with $p=0.75,0.5,0.25$ to ease the comparison.

Results – (O3). From Fig. 4, we see that the S-VNNs improve substantially on the computation time and increase the accuracy w.r.t. the nominal VNN.
This is particularly emphasized for Epilepsy, CNI, and MHEALTH, where most S-VNNs perform better than dense VNN, likely because the sample covariance matrices in these datasets contain spurious correlations.
On Realdisp, the accuracy improvement is smaller, indicating that in this dataset most correlations carry useful information to solve the task.
Both for thresholding and RCV, changing the value of the threshold or $p$ does not lead to large changes in accuracy, meaning that the majority of covariance values are not relevant for performance, while it affects the time efficiency of the model by increasing sparsity.
Finally, on 3 out of 4 datasets, stochastic sparsification achieves a slightly superior performance than deterministic thresholding, likely due to the regularization effect of different realizations of a sparsified covariance over different runs in a form akin to dropout-like techniques (Papp et al., ; Rong et al., ).

## 6 Related Works

Sparse PCA. The instability of the finite-sample covariance matrix is widely studied for PCA (Jolliffe, ; Paul, ; Baik et al., ) and various regularized estimators have been proposed including covariance shrinkage (Ledoit & Wolf, ; ), lasso penalties (Bien & Tibshirani, ; Han Liu & Zhao, ), and thresholding (Bickel & Levina, ; Deshpande & Montanari, ). We leverage here hard and soft thresholding given their benefits in PCA in low-data sparse-covariance settings and study their impact on VNN stability.

Stability of GNNs-VNNs. VNNs can be seen as graph convolutional neural networks operating on the covariance graph (Sihag et al., ), which draws analogies with PCA, but have been shown to be more stable extending the small perturbation analysis of GNNs (Gama et al., ; Kenlay et al., ; ; Levie et al., ; Maskey et al., ; Gao et al., ) to covariances.
Due to their robustness, they have proved successful in brain data processing (Sihag et al., ), transferability to large-dimension data (Sihag et al., ), temporal settings (Cavallo et al., ) and biased datasets (Cavallo et al., ).
All these results, however, hold for the dense covariance matrix which is suboptimal in low-data regimes and computationally heavy. Here, we propose S-VNNs and study their stability when the true covariance matrix is sparse and dense. This merits a different treatment than the perturbation assumption studied in (Sihag et al., ) and overcomes the limitations of dense VNNs. Our stochastically sparsified VNN draws analogies with the stability of GNNs to random link drops (Gao et al., ; ) but we generalize those findings to the more challenging case where all sparsified probabilities differ rather than being identical. We also propose design strategies for the sparsified probabilities based on our theoretical analysis.

Graph sparsification in GNNs. Graph sparsification techniques have been used for multiple purposes in GNNs but predominantly to improve: (i) scalability in large and dense graphs (Hamilton et al., ; Zhang et al., ; Zeng et al., ; Peng et al., ; Srinivasa et al., ; Ye & Ji, ); (ii) GNN expressiveness in graph classification tasks and alleviate oversmoothing (Papp et al., ; Rong et al., ; Fang et al., ; Zheng et al., ; Morris et al., ); and (iii) enhance interpretability (Rathee et al., ; Li et al., ; Naber et al., ).
Our S-VNNs aim on the one hand to improve scalability and on the other to preserve the tractability and theoretical links with covariance-based data processing.
While our sparsification approaches can be applied also to weighted graphs (see Appendix H), in this work we focus on covariance sparsification to characterize theoretically how this affects stability w.r.t. covariance estimation errors and spurious correlations.

## 7 Conclusion

In this work, we develop sparse covariance neural networks (S-VNNs) and study the effects of sparsification techniques on covariance neural networks.
We show that S-VNNs are more stable to the finite-data sample effect and more computationally efficient. In particular, when the true covariance is sparse, we propose S-VNNs with hard/soft thresholding and show that they
transfer to the nominal VNN with a rate inversely proportional to the square root of the number of data points. When the true covariance is dense, we put forth a stochastic sparsified approach to reduce the computational complexity while maintaining mathematical tractability and stability. In the framework of stochastic sparsification, we propose principled strategies, i.e., ACV and RCV, that align with the theoretical observations to tune sparsification impact and computation efficiency. Experimental results on four real datasets show that the proposed S-VNNs improve substantially on the computation time and achieve competitive or better performance compared with nominal VNNs, thus validating that spurious correlations are present and damaging in real datasets.
Although VNNs can be applied to several other regularized covariance estimators (Ledoit & Wolf, ; Bickel & Levina, ; Bien & Tibshirani, ; Friedman et al., ) to reduce the impact of spurious correlations, we here focus on thresholding since, unlike other estimators, it allows to control the sparsity of the covariance matrix (fundamental for time and memory efficiency), and allows for theoretical tractability.
Future work will investigate other covariance estimators and characterize the regularization effect of stochastic sparsification, which allows the model to observe multiple instances of the same covariance with potential benefits for robustness, expressiveness and reduced oversmoothing.

Reproducibility statement.
The proofs of all theoretical results are collected in the appendix.
For experimental results, we make the code available at the repository [https://github.com/andrea-cavallo-98/Sparse-VNN](https://github.com/andrea-cavallo-98/Sparse-VNN) and we provide references to all the datasets in the following. The preprocessing functions for all datasets are made available in the code.

- •
The synthetic datasets SparseCov, LargeCov and SmallCov are available in the repository.
- •
Realdisp can be downloaded at [https://archive.ics.uci.edu/dataset/305/realdisp+activity+recognition+dataset](https://archive.ics.uci.edu/dataset/305/realdisp+activity+recognition+dataset).
- •
MHEALTH can be downloaded at [https://archive.ics.uci.edu/dataset/319/mhealth+dataset](https://archive.ics.uci.edu/dataset/319/mhealth+dataset).
- •
CNI can be downloaded at [http://www.brainconnectivity.net/challenge.html](http://www.brainconnectivity.net/challenge.html).
- •
Epilepsy can be downloaded at [https://math.bu.edu/people/kolaczyk/datasets.html](https://math.bu.edu/people/kolaczyk/datasets.html) and we provide the preprocessed version in the repository.

We report additional experimental details for reproducibility in Sec. I.2.

Acknowledgements.
Part of this work was supported by the TU Delft AI Labs program, the NWO OTP GraSPA proposal #19497, and the NWO VENI proposal 222.032.

## Appendix A Background on PCA and Connection with VNNs

PCA is a change of basis that aims at minimizing the redundancy of the transformed data while maintaining their statistical properties of interest (Shlens, ; Jolliffe, ).
PCA achieves this by projecting the data onto a space where the covariance matrix of the transformed data is diagonal, i.e., the off-diagonal elements, which characterize the redundancy of the data, are all zeros.
More formally, given a dataset $\mathbf{X}\in\mathbb{R}^{N\times t}$, where each column is a data sample $\mathbf{x}\in\mathbb{R}^{N}$, PCA applies a transformation $\mathbf{\tilde{X}}=\mathbf{T}^{\mathsf{T}}\mathbf{X}$ with $\mathbf{T}\in\mathbb{R}^{N\times N}$ such that the covariance of the transformed data $\mathbf{\tilde{X}}\mathbf{\tilde{X}}^{\mathsf{T}}$ is diagonal.
This goal is achieved by setting $\mathbf{T}=\mathbf{\hat{V}}$, where $\mathbf{\hat{C}}=\mathbf{X}\mathbf{X}^{\mathsf{T}}=\mathbf{\hat{V}}\mathbf{
\hat{\Lambda}}\mathbf{\hat{V}}^{\mathsf{T}}$ is the eigendecomposition of the covariance matrix of the data, since $\mathbf{\tilde{X}}\mathbf{\tilde{X}}^{\mathsf{T}}=\mathbf{\hat{V}}^{\mathsf{T}
}\mathbf{X}\mathbf{X}^{\mathsf{T}}\mathbf{\hat{V}}=\mathbf{\hat{V}}^{\mathsf{T
}}(\mathbf{\hat{V}}\mathbf{\hat{\Lambda}}\mathbf{\hat{V}}^{\mathsf{T}})\mathbf
{\hat{V}}=\mathbf{\hat{\Lambda}}$ due to the orthonormality of the covariance eigenvectors. PCA is also commonly used for dimensionality reduction by selecting only some of the eigenvectors of the covariance matrix to define the new basis, and several techniques exist to select the number of components, such as looking at the eigenvalues of the covariance and selecting the number corresponding to the largest eigenvalue drop.
PCA is equivalent to the graph Fourier transform (GFT) for a graph defined by the covariance $\mathbf{\hat{C}}$ (Gama et al., ; Isufi et al., ). Indeed, the GFT projects data on the space defined by the eigenvectors of the graph shift operator $\mathbf{S}\in\mathbb{R}^{N\times N}$ (a matrix describing the graph structure, e.g., adjacency or Laplacian), i.e., $\mathbf{\tilde{X}}=\mathbf{U}^{\mathsf{T}}\mathbf{X}$ given the eigendecomposition $\mathbf{S}=\mathbf{U}\mathbf{W}\mathbf{U}^{\mathsf{T}}$.
If we consider the graph described by the sample covariance matrix, i.e., $\mathbf{S}=\mathbf{\hat{C}}$, then PCA and GFT perform the same operation.

The connection between PCA and VNNs emerges by looking at the action of the covariance filter in the spectral domain. For a data sample $\mathbf{x}$, taking the GFT of the covariance filter leads to $\mathbf{\tilde{u}}=\mathbf{\hat{V}}^{\mathsf{T}}\mathbf{u}=\mathbf{\hat{V}}^{
\mathsf{T}}\sum_{k=0}^{K}h_{k}[\mathbf{\hat{V}}\mathbf{\hat{\Lambda}}\mathbf{
\hat{V}}^{\mathsf{T}}]^{k}\mathbf{x}=\sum_{k=0}^{K}h_{k}\mathbf{\hat{\Lambda}}
^{k}\mathbf{\hat{V}}^{\mathsf{T}}\mathbf{x}$, equivalent to
$\sum_{k=0}^{K}h_{k}\hat{\lambda}_{i}^{k}\mathbf{\tilde{x}}_{i}=h(\hat{\lambda}
_{i})\mathbf{\tilde{x}}_{i}$ for the $i$-th entry; that is, the filter is a polynomial in the covariance eigenvalues that modulates each component of the GFT (i.e., PCA) of the signal $\mathbf{x}$ separately (cf. Sec. 3). This shows that the covariance filter processes the principal components of the data, and there exists a set of coefficients $h_{k}$ such that a covariance filter performs the same operation as PCA (Sihag et al., , Theorem 1).

## Appendix B Proof of Lemma 1

We provide a bound for the stability of PCA, i.e., $\|\mathbf{V}^{\mathsf{T}}\mathbf{x}-\mathbf{\hat{V}}^{\mathsf{T}}\mathbf{x}\|$. We have that

$$ $\displaystyle\begin{split}\left\|\mathbf{V}^{\mathsf{T}}\mathbf{x}-\mathbf{ \hat{V}}^{\mathsf{T}}\mathbf{x}\right\|&=\left\|\sum_{i=1}^{N}(\mathbf{v}_{i}- \mathbf{\hat{v}}_{i})^{\mathsf{T}}x_{i}\right\|\leq\sum_{i=1}^{N}\|\mathbf{v}_ {i}-\mathbf{\hat{v}}_{i}\||x_{i}|\\ &\leq\sum_{i=1}^{N}\|\mathbf{v}_{i}-\mathbf{\hat{v}}_{i}\|\leq N\max_{i}\left \|\mathbf{v}_{i}-\mathbf{\hat{v}}_{i}\right\|.\end{split}$ (3) $$

where we used triangle inequality, the fact that $\|\mathbf{x}\|\leq 1$ and the fact that $\|(\mathbf{v}_{i}-\mathbf{\hat{v}}_{i})^{\mathsf{T}}\|=\|\mathbf{v}_{i}-
\mathbf{\hat{v}}_{i}\|$.
Using the sin-theta theorem (Stewart & guang Sun, , Theorem V.3.6) (cf. equation 18), we have that

$$ $\displaystyle\max_{i}\left\|\mathbf{v}_{i}-\mathbf{\hat{v}}_{i}\right\|\leq \frac{\sqrt{2}\|\mathbf{E}\|}{\min_{i}|\lambda_{i}-\lambda_{i+1}|}$ (4) $$

where $\|\mathbf{E}\|=\|\mathbf{C}-\mathbf{\hat{C}}\|$.
We now leverage the result from (Vershynin, , Theorem 5.6.1) which shows that $\|\mathbf{E}\|\leq\mathcal{O}(t^{-1/2})$ with high probability, and this concludes the proof.
∎

## Appendix C Proof of Thms. 2 and 3

### C.1 Preliminaries

We begin by providing a Lemma that will be used in the following.

###### Lemma 2 .

Consider the true and the thresholded sample covariance matrices and their respective eigendecompositions $\mathbf{C}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{\mathsf{T}}$ and $\mathbf{\bar{C}}=\mathbf{\bar{V}}\bar{\mathbf{\Lambda}}\mathbf{\bar{V}}^{
\mathsf{T}}$. The following holds:

$$ $\displaystyle\frac{1}{\sqrt{2}}|\mathbf{v}_{j}^{\mathsf{T}}\mathbf{\bar{v}}_{i }|\leq\sqrt{1-|\mathbf{v}_{j}^{\mathsf{T}}\mathbf{\bar{v}}_{j}|}\quad\forall i \neq j.$ (5) $$

###### Proof.

We first observe that, due to the orthonormality of $\mathbf{V}$ and $\mathbf{\bar{V}}$, we have

$$ $\displaystyle\frac{1}{2}\|\mathbf{v}_{j}-\mathbf{\bar{v}}_{j}\|^{2}=\frac{1}{2 }(\mathbf{v}_{j}-\mathbf{\bar{v}}_{j})^{\mathsf{T}}(\mathbf{v}_{j}-\mathbf{ \bar{v}}_{j})=$ (6) $\displaystyle\frac{1}{2}(\mathbf{v}_{j}^{\mathsf{T}}\mathbf{v}_{j}-2\mathbf{v} _{j}^{\mathsf{T}}\mathbf{\bar{v}}_{j}+\mathbf{\bar{v}}_{j}^{\mathsf{T}}\mathbf {\bar{v}}_{j})=1-|\mathbf{v}_{j}^{\mathsf{T}}\mathbf{\bar{v}}_{j}|$ (7) $$

where we assumed $\mathbf{v}_{j}^{\mathsf{T}}\mathbf{\bar{v}}_{j}\geq 0$. Note that this assumption does not lose generality because eigenvectors are invariant to change of sign, i.e., if $\mathbf{v}_{j}$ is an eigenvector of $\mathbf{C}$, $-\mathbf{v}_{j}$ is equivalently an eigenvector.
Now, we write $\mathbf{v}_{j}-\mathbf{\bar{v}}_{j}=\delta\mathbf{v}_{j}$ as the sum of two components, $\delta\mathbf{v}_{j}=\delta\mathbf{v}_{j\parallel}+\delta\mathbf{v}_{j\perp}$, that are respectively parallel and perpendicular to $\mathbf{v}_{j}$. Due to triangle inequality, we have that $\|\delta\mathbf{v}_{j\perp}\|+\|\delta\mathbf{v}_{j\parallel}\|\leq\|\delta
\mathbf{v}_{j\perp}\|\leq\|\delta\mathbf{v}_{j}\|$. Since the columns of $\mathbf{V}$ form an orthonormal basis of $\mathbb{R}^{N}$, the perpendicular component of $\delta\mathbf{v}_{j}$ is the sum of the projections of $\mathbf{\bar{v}}_{j}$ on all columns of $\mathbf{V}$ except $\mathbf{v}_{j}$, i.e., $\delta\mathbf{v}_{j\perp}=\sum_{i=1,i\neq j}^{N}(\mathbf{\bar{v}}_{j}^{\mathsf
{T}}\mathbf{v}_{i})\mathbf{v}_{i}$ and the norm becomes $\|\delta\mathbf{v}_{j\perp}\|=\sqrt{\sum_{i=1,i\neq j}^{N}(\mathbf{\bar{v}}_{j
}^{\mathsf{T}}\mathbf{v}_{i})^{2}}$.
Due to triangle inequality and orthonormality of $\mathbf{V}$, for any $i\neq j$, we have

$$ $\displaystyle\frac{1}{\sqrt{2}}|\mathbf{\bar{v}}_{j}^{\mathsf{T}}\mathbf{v}_{i }|\leq\frac{1}{\sqrt{2}}\|\delta\mathbf{v}_{j\perp}\|\leq\frac{1}{\sqrt{2}}\| \mathbf{v}_{j}-\mathbf{\bar{v}}_{j}\|=\sqrt{1-|\mathbf{v}_{j}^{\mathsf{T}} \mathbf{\bar{v}}_{j}|}.$ (8) $$

∎

We now provide a general result on VNN stability that will be the starting point to prove Thms. 2 and 3.

###### Proposition 2 (Stability of VNNs) .

Consider a Lipschitz covariance filter $\mathbf{H}(\mathbf{C})$ with constant $P$.
Let $\mathbf{\bar{C}}$ be an estimated covariance matrix of $\mathbf{C}$. Then, for any generic signal $\mathbf{x}$ where $\|\mathbf{x}\|\leq 1$ and estimation error $\mathbf{E}=\mathbf{C}-\mathbf{\bar{C}}$ where $\|\mathbf{E}\|\ll 1$, it holds that

$$ $\displaystyle\|\mathbf{H}(\mathbf{\bar{C}})\mathbf{x}-\mathbf{H}(\mathbf{C}) \mathbf{x}\|\leq P\sqrt{N}\|\mathbf{E}\|(1+\sqrt{2N})+\mathcal{O}(\|\mathbf{E} \|^{2}).$ (9) $$

###### Proof.

We will make use of the eigendecompositions $\mathbf{C}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{\mathsf{T}}$ and $\mathbf{\bar{C}}=\mathbf{\bar{V}}\bar{\mathbf{\Lambda}}\mathbf{\bar{V}}^{
\mathsf{T}}$, where $\mathbf{V}$ and $\mathbf{\bar{V}}$ contain in their columns the eigenvectors of $\mathbf{C}$ and $\mathbf{\bar{C}}$, respectively, and $\mathbf{\Lambda}$ and $\bar{\mathbf{\Lambda}}$ contain their eigenvalues on the diagonal. Following (Sihag et al., , eq. (32)-(39)), under the assumption of small perturbation, i.e., $\|\mathbf{E}\|\ll 1$, the stability bound is the sum of three terms:

$$ $\displaystyle\mathbf{H}(\mathbf{\bar{C}})\mathbf{x}-\mathbf{H}(\mathbf{C}) \mathbf{x}\approx\underbrace{\sum_{i=0}^{N}\tilde{x}_{i}\sum_{k=0}^{K}h_{k} \sum_{r=0}^{k-1}\mathbf{C}^{r}\lambda_{i}^{k-r-1}(\bar{\lambda}_{i}-\lambda_{i })\mathbf{v}_{i}}_{\textnormal{term 1}}$ (10) $\displaystyle+\underbrace{\sum_{i=0}^{N}\tilde{x}_{i}\sum_{k=0}^{K}h_{k}\sum_{ r=0}^{k-1}\mathbf{C}^{r}\lambda_{i}^{k-r-1}(\lambda_{i}\mathbf{I}_{N}-\mathbf{ C})(\mathbf{\bar{v}}_{i}-\mathbf{v}_{i})}_{\textnormal{term 2}}$ (11) $\displaystyle+\underbrace{\sum_{i=0}^{N}\tilde{x}_{i}\sum_{k=0}^{K}h_{k}\sum_{ r=0}^{k-1}\mathbf{C}^{r}\lambda_{i}^{k-r-1}((\bar{\lambda}_{i}-\lambda_{i}) \mathbf{I}_{N}-\mathbf{E})(\mathbf{\bar{v}}_{i}-\mathbf{v}_{i})}_{\textnormal{ term 3}}$ (12) $$

where $\tilde{x}_{i}$ is the $i$-th component of the Covariance Fourier Transform of a generic graph signal $\mathbf{x}$, i.e., $\mathbf{\tilde{x}}=\mathbf{V}^{\mathsf{T}}\mathbf{x}$, and we inverted term 1 and 2 compared to (Sihag et al., ) for ease of explanation in the following. We now analyze the three terms separately.

Term 1.
Leveraging (Sihag et al., , eq. (59)-(61)), we have that

$$ $\displaystyle\sum_{i=0}^{N}\tilde{x}_{i}\sum_{k=0}^{K}h_{k}\sum_{r=0}^{k-1} \mathbf{C}^{r}\lambda_{i}^{k-r-1}(\bar{\lambda}_{i}-\lambda_{i})\mathbf{v}_{i} =\sum_{i=0}^{N}\tilde{x}_{i}h^{\prime}(\lambda_{i})(\bar{\lambda}_{i}-\lambda_ {i})\mathbf{v}_{i}$ (13) $$

where $h^{\prime}(\lambda_{i})$ is the derivative of the frequency response of the filter and is bounded in absolute value by $P$.
Using Weyl’s theorem (Golub & van Loan, , Theorem 8.1.6), we have that $|\bar{\lambda}_{i}-\lambda_{i}|\leq\|\mathbf{E}\|$.
Therefore, by taking the norm and using the fact that $\sum_{i=1}^{N}|\tilde{x}_{i}|\leq\sqrt{N}\|\mathbf{\tilde{x}}\|=\sqrt{N}\|
\mathbf{x}\|\leq\sqrt{N}$, we get

$$ $\displaystyle\left\|\sum_{i=0}^{N}\tilde{x}_{i}h^{\prime}(\lambda_{i})(\bar{ \lambda}_{i}-\lambda_{i})\mathbf{v}_{i}\right\|\leq\sum_{i=0}^{N}|\tilde{x}_{i }||h^{\prime}(\lambda_{i})|\|\mathbf{E}\|\|\mathbf{v}_{i}\|\leq$ (14) $\displaystyle\sum_{i=0}^{N}|\tilde{x}_{i}|P\|\mathbf{E}\|\leq$ (15) $\displaystyle P\sqrt{N}\|\mathbf{E}\|.$ (16) $$

Term 2.
Following (Sihag et al., , eq. (40)-(54)), the norm of term 2 is bounded by

$$ $\displaystyle\left\|\sum_{i=0}^{N}\tilde{x}_{i}\sum_{k=0}^{K}h_{k}\sum_{r=0}^{ k-1}\mathbf{C}^{r}\lambda_{i}^{k-r-1}(\lambda_{i}\mathbf{I}_{N}-\mathbf{C})( \mathbf{\bar{v}}_{i}-\mathbf{v}_{i})\right\|\leq\sqrt{N}\sum_{i=1}^{N}|\tilde{ x}_{i}|\max_{j}|h(\lambda_{i})-h(\lambda_{j})||\mathbf{v}_{j}^{\mathsf{T}} \mathbf{\bar{v}}_{i}|$ (17) $$

where $h(\lambda)$ is the frequency response of the covariance filter. Note that in equation 17 we have a term $\sqrt{N}$ that does not appear in (Sihag et al., , eq. (54)) because we consider operator norm instead of uniform norm.
We now leverage the sin-theta theorem (Stewart & guang Sun, , Theorem V.3.6): for $j=1,\dots,N$, we have that

$$ $\frac{1}{2}\|\mathbf{v}_{j}-\mathbf{\bar{v}}_{j}\|^{2}=1-|\mathbf{v}_{j}^{ \mathsf{T}}\mathbf{\bar{v}}_{j}|\leq\frac{\|\mathbf{E}\|^{2}}{\min_{i}|\lambda _{i}-\lambda_{i+1}|^{2}}.$ (18) $$

By using Lemma 2, we obtain:

$$ $\displaystyle|\mathbf{v}_{j}^{\mathsf{T}}\mathbf{\bar{v}}_{i}|\leq\sqrt{2} \sqrt{1-|\mathbf{v}_{j}^{\mathsf{T}}\mathbf{\bar{v}}_{j}|}\leq\frac{\sqrt{2}\| \mathbf{E}\|}{\min_{i}|\lambda_{i}-\lambda_{i+1}|}\quad\forall i\neq j.$ (19) $$

By plugging equation 19 into equation 17, we have

$$ $\displaystyle\sqrt{N}\sum_{i=1}^{N}|\tilde{x}_{i}|\max_{j\neq i}|h(\lambda_{i} )-h(\lambda_{j})||\mathbf{v}_{j}^{\mathsf{T}}\mathbf{\bar{v}}_{i}|\leq$ (20) $\displaystyle\sqrt{2N}\|\mathbf{E}\|\sum_{i=1}^{N}|\tilde{x}_{i}|\max_{j\neq i }\frac{|h(\lambda_{i})-h(\lambda_{j})|}{\min_{i}|\lambda_{i}-\lambda_{i+1}|}\leq$ (21) $\displaystyle\sqrt{2N}\|\mathbf{E}\|\sum_{i=1}^{N}|\tilde{x}_{i}|\max_{j\neq i }\frac{|h(\lambda_{i})-h(\lambda_{j})|}{|\lambda_{i}-\lambda_{j}|}\leq$ (22) $\displaystyle\sqrt{2N}\|\mathbf{E}\|\sum_{i=1}^{N}|\tilde{x}_{i}|P\leq$ (23) $\displaystyle\sqrt{2}N\|\mathbf{E}\|P$ (24) $$

where we used the Lipschitz property of the filter and the fact that $\sum_{i=1}^{N}|\tilde{x}_{i}|\leq\sqrt{N}\|\mathbf{x}\|$.

Term 3.
This term depends on the second-order error $((\bar{\lambda}_{i}-\lambda_{i})\mathbf{I}_{N}-\mathbf{E})(\mathbf{\bar{v}}_{i
}-\mathbf{v}_{i})$.
From Weyl’s theorem (Golub & van Loan, , Theorem 8.1.6), we have that $\|(\bar{\lambda}_{i}-\lambda_{i})\mathbf{I}_{N}-\mathbf{E}\|\leq 2\|\mathbf{E}\|$.
From equation 18, we have that $\|\mathbf{\bar{v}}_{i}-\mathbf{v}_{i}\|=\mathcal{O}(\|\mathbf{E}\|)$. Therefore, the product of the two is $\|(\bar{\lambda}_{i}-\lambda_{i})\mathbf{I}_{N}-\mathbf{E}\|\|\mathbf{\bar{v}}
_{i}-\mathbf{v}_{i}\|\approx\mathcal{O}(\|\mathbf{E}\|^{2})$, which is dominated by the other two terms in equation 10 and equation 11 under small perturbation assumption.

We complete the proof by merging the terms in equation 24 and equation 16.

∎

###### Lemma 2 .

###### Proof.

###### Proposition 2 (Stability of VNNs) .

###### Proof.

### C.2 Proof of Thm. 2

From (Bickel & Levina, , Theorem 1), given a true covariance $\mathbf{C}$ and a hard-thresholded sample covariance $\mathbf{\bar{C}}$, it holds with high probability that

$$ $\displaystyle\|\mathbf{\bar{C}}-\mathbf{C}\|=\|\mathbf{E}\|\leq c_{0}\sqrt{ \frac{\log N}{t}}.$ (25) $$

By replacing $\|\mathbf{E}\|$ in equation 9, the claim follows.
∎

### C.3 Proof of Thm. 3

From (Deshpande & Montanari, , Theorem 1), given a true covariance $\mathbf{C}$ and a soft-thresholded sample covariance $\mathbf{\bar{C}}$, it holds with high probability that

$$ $\displaystyle\|\mathbf{\bar{C}}-\mathbf{C}\|=\|\mathbf{E}\|\leq\frac{1}{\sqrt{ t}}Cc_{0}\max(1,\lambda_{0})\sqrt{\max\left(\log\frac{N}{c_{0}^{2}},1\right)}$ (26) $$

for a generic constant $C$.
By replacing $\|\mathbf{E}\|$ in equation 9, the claim follows.
∎

## Appendix D Proof of Proposition 1

We have that

$$ $\displaystyle\begin{split}\left\|\mathbf{V}^{\mathsf{T}}\mathbf{x}-\mathbf{ \bar{V}}^{\mathsf{T}}\mathbf{x}\right\|&=\left\|\sum_{i=1}^{N}(\mathbf{v}_{i}- \mathbf{\bar{v}}_{i})^{\mathsf{T}}x_{i}\right\|\leq\sum_{i=1}^{N}\|\mathbf{v}_ {i}-\mathbf{\bar{v}}_{i}\||x_{i}|\\ &\leq\sum_{i=1}^{N}\|\mathbf{v}_{i}-\mathbf{\bar{v}}_{i}\|\leq N\max_{i}\left \|\mathbf{v}_{i}-\mathbf{\bar{v}}_{i}\right\|.\end{split}$ (27) $$

where we used triangle inequality, the fact that $\|\mathbf{x}\|\leq 1$ and the fact that $\|(\mathbf{v}_{i}-\mathbf{\bar{v}}_{i})^{\mathsf{T}}\|=\|\mathbf{v}_{i}-
\mathbf{\bar{v}}_{i}\|$.
Using the sin-theta theorem (Stewart & guang Sun, , Theorem V.3.6) (cf. equation 18), we have that $\max_{i}\left\|\mathbf{v}_{i}-\mathbf{\bar{v}}_{i}\right\|\leq\sqrt{2}\|
\mathbf{E}\|/\min_{i}|\lambda_{i}-\lambda_{i+1}|$ where $\|\mathbf{E}\|=\|\mathbf{C}-\mathbf{\bar{C}}\|$. Now, we use the expression for $\|\mathbf{E}\|$ with hard thresholding in (Bickel & Levina, , Theorem 1), i.e.,

$$ $\displaystyle\|\mathbf{E}\|\leq c_{0}\sqrt{\frac{\log N}{t}},$ (28) $$

and the claim follows.
∎

## Appendix E Generalized Covariance Filter Frequency Response

We provide here additional details related to the derivation of the generalized covariance filter frequency response in Def. 6 following (Gao et al., ).
Consider a covariance filter operating on a sequence of $K$ sparsified covariances $\mathbf{\tilde{C}}_{1}\dots\mathbf{\tilde{C}}_{K}$ with eigendecomposition $\mathbf{\tilde{C}}_{k}=\mathbf{\tilde{V}}_{k}\mathbf{\tilde{\mathbf{\Lambda}}}
_{k}\mathbf{\tilde{V}}_{k}^{\mathsf{T}}$ (where $\mathbf{\tilde{V}}_{k}=[\mathbf{\tilde{v}}_{k1},\dots,\mathbf{\tilde{v}}_{kN}]$ and $\mathbf{\tilde{\mathbf{\Lambda}}}_{k}=\text{diag}(\tilde{\lambda}_{k1},\dots,
\tilde{\lambda}_{kN})$) according to Def. 5.
For matrix $\mathbf{\tilde{C}}_{1}$, we can express a signal $\mathbf{x}$ as $\mathbf{x}=\sum_{i_{1}=1}^{N}\hat{x}_{1i_{1}}\mathbf{\tilde{v}}_{1i_{1}}$, where $\mathbf{\hat{x}}_{1}=[\hat{x}_{11},\dots,\hat{x}_{1N}]^{\mathsf{T}}$ is the covariance Fourier transform of $\mathbf{x}$ w.r.t. $\mathbf{\tilde{C}}_{1}$.
By performing a graph signal shift, we obtain

$$ $\displaystyle\mathbf{x}^{(1)}=\mathbf{\tilde{C}}_{1}\mathbf{x}=\mathbf{\tilde{ C}}_{1}\sum_{i_{1}=1}^{N}\hat{x}_{1i_{1}}\mathbf{\tilde{v}}_{1i_{1}}=\sum_{i_{ 1}=1}^{N}\hat{x}_{1i_{1}}\tilde{\lambda}_{1i_{1}}\mathbf{\tilde{v}}_{1i_{1}}.$ (29) $$

When performing a second shift, i.e., $\mathbf{x}^{(2)}=\mathbf{\tilde{C}}_{2}\mathbf{\tilde{C}}_{1}\mathbf{x}$, we decompose each eigenvector $\mathbf{\tilde{v}}_{1i_{1}}$ by taking its covariance Fourier transform w.r.t. $\mathbf{\tilde{C}}_{2}$, i.e., we get $\mathbf{\tilde{v}}_{1i_{1}}=\sum_{i_{2}=1}^{N}\hat{x}_{2i_{1}i_{2}}\mathbf{
\tilde{v}}_{2i_{2}}$ where $\mathbf{\tilde{x}}_{2i_{1}}=[\hat{x}_{2i_{1}1},\dots,\hat{x}_{2i_{1}N}]$ is the covariance Fourier transform of $\mathbf{\hat{v}}_{1i_{1}}$ over $\mathbf{\tilde{C}}_{2}$.
Performing this operation on all eigenvectors $\mathbf{\tilde{v}}_{11},\dots,\mathbf{\tilde{v}}_{1N}$ and using equation 29, we can write

$$ $\displaystyle\mathbf{x}^{(2)}=\mathbf{\tilde{C}}_{2}\sum_{i_{1}=1}^{N}\hat{x}_ {1i_{1}}\tilde{\lambda}_{1i_{1}}\mathbf{\tilde{v}}_{1i_{1}}=\sum_{i_{2}=1}^{N} \sum_{i_{1}=1}^{N}\hat{x}_{2i_{1}i_{2}}\hat{x}_{1i_{1}}\tilde{\lambda}_{2i_{2} }\tilde{\lambda}_{1i_{1}}\mathbf{\tilde{v}}_{1i_{1}}.$ (30) $$

Generalizing this to $k$ shifts, we get

$$ $\displaystyle\mathbf{x}^{(k)}=\sum_{i_{k}=1}^{N}\dots\sum_{i_{1}=1}^{N}\hat{x} _{ki_{k-1}i_{k}}\dots\hat{x}_{2i_{1}i_{2}}\hat{x}_{1i_{1}}\prod_{j=1}^{k} \tilde{\lambda}_{ji_{j}}\mathbf{\tilde{v}}_{ki_{k}}$ (31) $$

and by aggregating the $K+1$ shifted signals according to Def. 5 we obtain

$$ $\displaystyle\mathbf{\tilde{u}}=\sum_{i_{K}=1}^{N}\dots\sum_{i_{1}=1}^{N}\hat{ x}_{Ki_{K-1}i_{K}}\dots\hat{x}_{2i_{1}i_{2}}\hat{x}_{1i_{1}}\sum_{k=1}^{K}h_{k }\prod_{j=1}^{k}\tilde{\lambda}_{ji_{j}}\mathbf{\tilde{v}}_{ki_{k}}.$ (32) $$

From equation 32 we see that this representation involves all eigenvalues $\mathbf{\tilde{\mathbf{\Lambda}}}_{K},\dots,\mathbf{\tilde{\mathbf{\Lambda}}}_
{1}$ and eigenvectors $\mathbf{\tilde{V}}_{K},\dots,\mathbf{\tilde{V}}_{1}$ of the sequence of covariance realizations. Therefore, we can consider the coefficients $\{\hat{x}_{1i_{1}}\}_{i_{1}=1}^{N}$ and $\{\hat{x}_{2i_{j}i_{j+1}}\}_{j=1}^{K-1}$ as the generalized covariance Fourier transform of signal $\mathbf{x}$ on the sequence of sparsified covariances $\mathbf{\tilde{V}}_{K},\dots,\mathbf{\tilde{V}}_{1}$. This supports the definition of a generalized frequency response for a covariance filter over random sparsified covariances in Def. 6.

## Appendix F Proof of Thm. 4

Consider a covariance filter $\mathbf{H}(\mathbf{\tilde{C}})\mathbf{x}$ operating on the matrix $\mathbf{\tilde{C}}=\mathbf{E}+\mathbf{\hat{C}}$, which represents a sparsified random matrix that is a copy of $\mathbf{\hat{C}}$ whose elements $c_{ij}$ are dropped with probability $(1-p_{ij})$, where $\mathbf{\hat{C}}$ is the sample covariance matrix and $\mathbf{C}$ is the true covariance matrix. Then, the stability of the covariance filter depends on two sources of error. Indeed, by adding and subtracting $\mathbf{H}(\mathbf{\hat{C}})\mathbf{x}$ and leveraging the triangle inequality, we obtain

$$ $\displaystyle\|\mathbf{H}(\mathbf{C})\mathbf{x}-\mathbf{H}(\mathbf{\tilde{C}}) \mathbf{x}\|^{2}=\|\mathbf{H}(\mathbf{C})\mathbf{x}-\mathbf{H}(\mathbf{\hat{C} })\mathbf{x}+\mathbf{H}(\mathbf{\hat{C}})\mathbf{x}-\mathbf{H}(\mathbf{\tilde{ C}})\mathbf{x}\|^{2}\leq$ (33) $\displaystyle\|\mathbf{H}(\mathbf{C})\mathbf{x}-\mathbf{H}(\mathbf{\hat{C}}) \mathbf{x}\|^{2}+\|\mathbf{H}(\mathbf{\hat{C}})\mathbf{x}-\mathbf{H}(\mathbf{ \tilde{C}})\mathbf{x}\|^{2}=\alpha+\beta$ (34) $$

where $\alpha$ is an instability error due to the uncertainties in the covariance matrix estimate and $\beta$ depends on the stochastic sparsification of the sample covariance matrix. We analyze and provide an expression for these two terms in the remainder of the proof.

### F.1 Covariance uncertainty error

The effect of covariance estimation errors on the stability of VNN is deterministic and, consequently, can be bounded by the square of the bound in Thm. 1:

$$ $\displaystyle\|\mathbf{H}(\mathbf{C})\mathbf{x}-\mathbf{H}(\mathbf{\hat{C}}) \mathbf{x}\|^{2}\leq\beta^{2}=\frac{P^{2}}{t}\mathcal{O}\left(N+\frac{\| \mathbf{C}\|^{2}\log(Nt)}{\nu^{2}t^{2}}\right)$ (35) $$

where $t$ is the number of data samples used to estimate the covariance matrix.

### F.2 Covariance sparsification error

To analyze the stability of VNN to the stochastic sparsification of the covariance matrix, we leverage and extend previous results on the stability of GNNs to stochastic graph perturbations.
We begin by providing some lemmas that will be used in the main statement.

###### Lemma 3 .

Consider a random covariance matrix $\mathbf{\tilde{C}}_{r}=\mathbf{\hat{C}}+\mathbf{E}_{r}$, where $\mathbf{\tilde{C}}_{r}$ is a copy of the sample covariance matrix $\mathbf{\hat{C}}$ whose elements $\hat{c}_{ij}$ are dropped independently with probability $1-p_{ij}$ and $\mathbf{E}_{r}$ is its distance from $\mathbf{\hat{C}}$. Under the conditions of Def. 4, i.e., edge perturbations are undirected (i.e., $\mathbf{E}_{r}$ is symmetric) and no perturbations occur on the diagonal (i.e., $\mathbf{E}_{r}$ has zeros on the diagonal), it holds that

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r}^{2}])=\sum_{i=1}^{N}\sum _{n=1}^{N}\hat{c}_{in}^{2}(1-p_{in})=Q.$ (36) $$

###### Proof.

Each entry of the random error matrix $\mathbf{E}_{r}$ can be represented as $[\mathbf{E}_{r}]_{ij}=-\delta_{ij}c_{ij}$, where $\delta_{ij}$ is a Bernoulli variable that is one with probability $1-p_{ij}$ and zero with probability $p_{ij}$. By performing the matrix multiplication, we can express each element of $\mathbf{E}_{r}^{2}$ and, consequently, its expectation, as

$$ $\displaystyle[\mathbf{E}_{r}^{2}]_{ij}=\sum_{n=1}^{N}\hat{c}_{in}\hat{c}_{nj} \delta_{in}\delta_{nj},\quad[\mathbb{E}[\mathbf{E}_{r}^{2}]]_{ij}=\sum_{n=1}^{ N}\hat{c}_{in}\hat{c}_{nj}\mathbb{E}[\delta_{in}\delta_{nj}].$ (37) $$

The Bernoulli variables $\delta_{ij}$ are independent except for $\delta_{ij}=\delta_{ji}$ given the symmetry of $\mathbf{E}_{r}$. Therefore, we have

$$ $\displaystyle\mathbb{E}[\delta_{in}\delta_{nj}]=\begin{cases}(1-p_{in})(1-p_{ nj})&\text{if }i\neq j\\ (1-p_{in})&\text{if }i=j\end{cases}$ (38) $$

We can now compute the trace by summing the elements on the diagonal and using the fact that $\mathbf{E}_{r}$ is symmetric, i.e.,

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r}^{2}])=\sum_{i=1}^{N}\sum _{n=1}^{N}\hat{c}_{in}\hat{c}_{ni}\mathbb{E}[\delta_{in}\delta_{ni}]=\hat{c}_{ in}^{2}(1-p_{in}).$ (39) $$

∎

###### Lemma 4 .

Consider a covariance filter $\mathbf{H}(\mathbf{C})$ with coefficients $\{h_{k}\}_{k=0}^{K}$ and generalized integral Lipschitz frequency response with constant $P$. Given some realizations of a random matrix $\mathbf{\tilde{C}}_{r}=\mathbf{E}_{r}+\mathbf{\hat{C}}$, for any signal $\mathbf{x}$, it holds that

$$ $\displaystyle\mathbb{E}\left[\sum_{r=1}^{K}\mathsf{trace}\left(\sum_{k=r}^{K} \sum_{l=r}^{K}h_{k}h_{l}\mathbf{E}_{r}\mathbf{\hat{C}}^{k+l-2r}\mathbf{E}_{r} \mathbf{\hat{C}}^{r-1}\mathbf{x}\mathbf{x}^{\mathsf{T}}\mathbf{\hat{C}}^{r-1} \right)\right]\leq NP^{2}\|\mathbf{x}\|^{2}Q$ (40) $$

with $Q$ defined in equation 36.

###### Proof.

Following (Gao et al., , eq. (B.13)-(B.15)), we rewrite the term in equation 40 as

$$ $\displaystyle\sum_{i=1}^{N}\hat{x}_{i}^{2}\sum_{r=1}^{K}\mathsf{trace}\left( \sum_{k=r}^{K}\sum_{l=r}^{K}h_{k}h_{l}\hat{\lambda}_{i}^{2r-2}\mathbf{\hat{C}} ^{k+l-2r}\mathbb{E}[\mathbf{E}_{r}^{2}]\right)$ (41) $$

where $\hat{x}_{i}$ is the $i$-th component of the Covariance Fourier Transform of the signal $\mathbf{\hat{x}}=\mathbf{\hat{V}}^{\mathsf{T}}\mathbf{x}$ and $\mathbf{\hat{C}}=\mathbf{\hat{V}}\mathbf{\hat{\Lambda}}\mathbf{\hat{V}}^{
\mathsf{T}}$ is the eigendecomposition of $\mathbf{\hat{C}}$ with eigenvectors $\mathbf{\hat{V}}=[\mathbf{\hat{v}}_{1},\dots,\mathbf{\hat{v}}_{N}]^{\mathsf{T}}$ and eigenvalues $\mathbf{\hat{\Lambda}}=\text{diag}(\hat{\lambda}_{1},\dots,\hat{\lambda}_{N})$.
Now, we leverage the property

$$ $\displaystyle\mathsf{trace}(\mathbf{A}\mathbf{B})\leq\frac{\|\mathbf{A}+ \mathbf{A}^{\mathsf{T}}\|}{2}\mathsf{trace}(\mathbf{B})\leq\|\mathbf{A}\| \mathsf{trace}(\mathbf{B})$ (42) $$

which holds for any square matrix $\mathbf{A}$ and positive semi-definite matrix $\mathbf{B}$ (Wang et al., ).
We note that $\mathbb{E}[\mathbf{E}_{r}^{2}]$ is positive semi-definite since $\mathbf{E}_{r}^{2}$ is the square of a symmetric matrix. Therefore, we use equation 42 to write

$$ $\displaystyle\sum_{i=1}^{N}\hat{x}_{i}^{2}\sum_{r=1}^{K}\mathsf{trace}\left( \sum_{k=r}^{K}\sum_{l=r}^{K}h_{k}h_{l}\hat{\lambda}_{i}^{2r-2}\mathbf{\hat{C}} ^{k+l-2r}\mathbb{E}[\mathbf{E}_{r}^{2}]\right)\leq$ (43) $\displaystyle\sum_{i=1}^{N}\hat{x}_{i}^{2}\left\|\sum_{r=1}^{K}\sum_{k=r}^{K} \sum_{l=r}^{K}h_{k}h_{l}\hat{\lambda}_{i}^{2r-2}\mathbf{\hat{C}}^{k+l-2r} \right\|\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r}^{2}])$ (44) $$

where we also used the linearity of the trace operator to move the summation term.

From (Gao et al., , eq. (B.24)), we have that the term within the norm is bounded by the generalized Lipschitz coefficient of the filter, i.e.,

$$ $\displaystyle\left\|\sum_{r=1}^{K}\sum_{k=r}^{K}\sum_{l=r}^{K}h_{k}h_{l}\hat{ \lambda}_{i}^{2r-2}\mathbf{\hat{C}}^{k+l-2r}\right\|\leq P^{2}.$ (45) $$

Then, using Lemma 3, we have that

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r}^{2}])=\sum_{i=1}^{N}\sum _{n=1}^{N}\hat{c}_{in}^{2}(1-p_{in})=Q.$ (46) $$

Substituting these two terms in equation 44 leads to

$$ $\displaystyle\sum_{i=1}^{N}\hat{x}_{i}^{2}P^{2}Q\leq NP^{2}Q\sum_{i=1}^{N}\hat {x}_{i}^{2}=NP^{2}Q\|\mathbf{x}\|^{2},$ (47) $$

where we used the fact that $\sum_{i=1}^{N}\hat{x}_{i}^{2}=\|\mathbf{\hat{x}}\|^{2}=\|\mathbf{V}^{\mathsf{T
}}\mathbf{x}\|^{2}=\|\mathbf{x}\|^{2}$, i.e., the covariance Fourier transform does not modify the norm of the signal.
∎

###### Lemma 5 .

Consider two distinct random covariance realizations $\mathbf{\tilde{C}}_{r_{1}}=\mathbf{\hat{C}}+\mathbf{E}_{r_{1}}$ and $\mathbf{\tilde{C}}_{r_{2}}=\mathbf{\hat{C}}+\mathbf{E}_{r_{2}}$. Then, the following holds:

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r_{1}}\mathbf{E}_{r_{2}}])= \sum_{i=1}^{N}\sum_{n=1}^{N}\hat{c}_{in}^{2}(1-p_{in})^{2}$ (48) $$

###### Proof.

Similarly to Lemma 3, we can express each element of the expected value of the matrix product $\mathbf{E}_{r_{1}}\mathbf{E}_{r_{2}}$ as

$$ $\displaystyle[\mathbb{E}[\mathbf{E}_{r_{1}}\mathbf{E}_{r_{2}}]]_{ij}=\sum_{n=1 }^{N}\hat{c}_{in}\hat{c}_{nj}\mathbb{E}[\delta_{r_{1},in}\delta_{r_{2},nj}]$ (49) $$

where $\delta_{r,ij}$ is a Bernoulli variable relative to realization $\mathbf{E}_{r}$ that is one with probability $1-p_{ij}$ and zero with probability $p_{ij}$. Since $\mathbf{E}_{r_{1}}$ and $\mathbf{E}_{r_{2}}$ are two different realizations, all variables $\delta_{r,ij}$ are independent. Therefore, $\mathbb{E}[\delta_{r_{1},in}\delta_{r_{2},nj}]=(1-p_{in})(1-p_{nj})$ and the trace becomes

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r_{1}}\mathbf{E}_{r_{2}}])= \sum_{i=1}^{N}\sum_{n=1}^{N}\hat{c}_{in}^{2}(1-p_{in})^{2}.$ (50) $$

∎

###### Lemma 6 .

Consider two distinct random covariance realizations $\mathbf{\tilde{C}}_{r_{1}}=\mathbf{\hat{C}}+\mathbf{E}_{r_{1}}$ and $\mathbf{\tilde{C}}_{r_{2}}=\mathbf{\hat{C}}+\mathbf{E}_{r_{2}}$. Then, the following holds:

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r_{1}}^{2}\mathbf{E}_{r_{2} }])=\mathcal{O}((1-p_{1})(1-p_{2})).$ (51) $$

where $p_{1},p_{2}$ are two generic probabilities.

###### Proof.

From Lemma 3 we have that

$$ $\displaystyle[\mathbb{E}[\mathbf{E}_{r_{1}}^{2}]]_{ij}=\sum_{n=1}^{N}\hat{c}_{ in}\hat{c}_{nj}\mathbb{E}[\delta_{r_{1},in}\delta_{r_{1},nj}].$ (52) $$

Therefore, the elements of $\mathbf{E}_{r_{1}}^{2}\mathbf{E}_{r_{2}}$ are

$$ $\displaystyle[\mathbb{E}[\mathbf{E}_{r_{1}}^{2}\mathbf{E}_{r_{2}}]]_{ij}=\sum_ {m=1}^{N}\mathbb{E}[[\mathbf{E}_{r_{1}}^{2}]_{im}\hat{c}_{mj}\delta_{r_{2},mj} ]=\sum_{m=1}^{N}\sum_{n=1}^{N}\hat{c}_{in}\hat{c}_{nm}\hat{c}_{mj}\mathbb{E}[ \delta_{r_{1},in}\delta_{r_{1},nm}\delta_{r_{2},mj}].$ (53) $$

The three variables $\delta_{r_{1},in},\delta_{r_{1},nm},\delta_{r_{2},mj}$ are all independent with the exception of $\delta_{r_{1},in}=\delta_{r_{1},ni}$ due to matrix symmetry.
Therefore,

$$ $\displaystyle\mathbb{E}[\delta_{r_{1},in}\delta_{r_{1},nm}\delta_{r_{2},mj}]= \begin{cases}(1-p_{in})(1-p_{nm})(1-p_{mj})&\text{if }i\neq m\\ (1-p_{in})(1-p_{mj})&\text{if }i=m\end{cases}$ (54) $$

We can now compute the trace:

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r_{1}}^{2}\mathbf{E}_{r_{2} }])=\sum_{i=1}^{N}\sum_{m=1}^{N}\sum_{n=1}^{N}\hat{c}_{in}\hat{c}_{nm}\hat{c}_ {mi}\mathbb{E}[\delta_{r_{1},in}\delta_{r_{1},nm}\delta_{r_{2},mi}].$ (55) $$

Therefore, the trace is a sum of terms of quadratic and cubic order in the probability value, i.e., $\mathcal{O}((1-p_{1})(1-p_{2}))$ or $\mathcal{O}((1-p_{1})(1-p_{2})(1-p_{3}))$, where $p_{1},p_{2},p_{3}$ are generic probabilities values and the quadratic terms dominate the behavior since $1-p<1$. So,

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r_{1}}^{2}\mathbf{E}_{r_{2} }])=\mathcal{O}((1-p_{1})(1-p_{2})).$ (56) $$

∎

###### Lemma 7 .

Consider the sample covariance matrix $\mathbf{\hat{C}}$ and a random sparsified covariance $\mathbf{\tilde{C}}_{r}=\mathbf{\hat{C}}+\mathbf{E}_{r}$ as in Def. 4. We have that

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r}^{3}])=\mathcal{O}((1-p_{ 1})(1-p_{2})).$ (57) $$

where $p_{1},p_{2}$ are two generic probabilities.

###### Proof.

Similarly to Lemma 3, we can write the expected value of each element of $\mathbf{E}^{3}_{r}$ as

$$ $\displaystyle[\mathbb{E}[\mathbf{E}^{3}_{r}]]_{ij}=\sum_{m=1}^{N}\mathbb{E}[[ \mathbf{E}_{r}^{2}]_{im}\hat{c}_{mj}\delta_{mj}]=\sum_{m=1}^{N}\sum_{n=1}^{N} \hat{c}_{in}\hat{c}_{nm}\hat{c}_{mj}\mathbb{E}[\delta_{in}\delta_{nm}\delta_{ mj}]$ (58) $$

and, consequently, the trace

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r}^{3}])=\sum_{i=1}^{N}\sum _{m=1}^{N}\sum_{n=1}^{N}\hat{c}_{in}\hat{c}_{nm}\hat{c}_{mi}\mathbb{E}[\delta_ {in}\delta_{nm}\delta_{mi}]$ (59) $$

where the Bernoulli variables $\delta_{in},\delta_{nm},\delta_{mi}$ are all independent with the exception of $\delta_{in}=\delta_{ni}$ due to matrix symmetry and the terms where $i=n=m=j$. Specifically, we have that

$$ $\displaystyle\mathbb{E}[\delta_{in}\delta_{nm}\delta_{mi}]=\begin{cases}(1-p_{ in})(1-p_{nm})(1-p_{mi})&\text{if }i\neq m\land i\neq n\\ (1-p_{in})(1-p_{mi})&\text{if }i=m\oplus i=n\\ (1-p_{in})&\text{if }i=m=n\end{cases}$ (60) $$

where $\oplus$ denotes the xor operator.
Now we note that, according to Def. 4, the matrix $\mathbf{E}_{r}$ has zeros on the diagonal or, equivalently, $1-p_{ii}=0\quad\forall i$. As a consequence, the terms where $i=m=n$ in equation 59 are all zeros and the trace only contains terms of quadratic and cubic order in the probability value, i.e., $\mathcal{O}((1-p_{1})(1-p_{2}))$ or $\mathcal{O}((1-p_{1})(1-p_{2})(1-p_{3}))$, where $p_{1},p_{2},p_{3}$ are generic probabilities values and the quadratic terms dominate the behavior since $1-p<1$. Therefore,

$$ $\displaystyle\mathsf{trace}(\mathbb{E}[\mathbf{E}_{r}^{3}])=\mathcal{O}((1-p_{ 1})(1-p_{2})).$ (61) $$

∎

Main statement.
Let $\mathbf{u}=\mathbf{H}(\mathbf{\hat{C}})\mathbf{x}$ and $\mathbf{\tilde{u}}=\mathbf{H}(\mathbf{\tilde{C}})\mathbf{x}$ be the outputs of the deterministic and stochastic filter, respectively. We are interested in the term

$$ $\displaystyle\mathbb{E}[\|\mathbf{u}-\mathbf{\tilde{u}}\|^{2}]=\mathbb{E}[ \mathsf{trace}(\mathbf{u}^{\mathsf{T}}\mathbf{u}+\mathbf{\tilde{u}}^{\mathsf{T }}\mathbf{\tilde{u}}-2\mathbf{u}^{\mathsf{T}}\mathbf{\tilde{u}})]=$ (62) $\displaystyle\mathbb{E}[\mathsf{trace}(\mathbf{\tilde{u}}^{\mathsf{T}}\mathbf{ \tilde{u}}-\mathbf{u}^{\mathsf{T}}\mathbf{u})]+2\mathbb{E}[\mathsf{trace}( \mathbf{u}^{\mathsf{T}}\mathbf{u}-\mathbf{u}^{\mathsf{T}}\mathbf{\tilde{u}})]$ (63) $$

where we added and subtracted $\mathbf{u}^{\mathsf{T}}\mathbf{u}$ and used linearity of expectation and trace.

We represent a random covariance matrix as $\mathbf{\tilde{C}}_{r}=\mathbf{\hat{C}}+\mathbf{E}_{r}$, where $\mathbf{E}_{r}$ is a random matrix that contains the deviation from the true covariance.
Following (Gao et al., , eq. (B.2)-(B.6)), we express the first term in equation 63 as

$$ $\displaystyle\mathbb{E}[\mathsf{trace}(\mathbf{\tilde{u}}^{\mathsf{T}}\mathbf{ \tilde{u}}-\mathbf{u}^{\mathsf{T}}\mathbf{u})]=-2\mathbb{E}[\mathsf{trace}( \mathbf{u}^{\mathsf{T}}\mathbf{u}-\mathbf{u}^{\mathsf{T}}\mathbf{\tilde{u}})]+$ (64) $\displaystyle\sum_{k=1}^{K}\sum_{l=1}^{K}h_{k}h_{l}\mathsf{trace}\left(\mathbb {E}\left[\sum_{r=1}^{\min(k,l)}\mathbf{\hat{C}}^{k-r}\mathbf{E}_{r}\mathbf{ \hat{C}}^{r-1}\mathbf{x}\mathbf{x}^{\mathsf{T}}\mathbf{\hat{C}}^{r-1}\mathbf{E }_{r}\mathbf{\hat{C}}^{l-r}\right]\right)+$ (65) $\displaystyle\sum_{k=0}^{K}\sum_{l=0}^{K}h_{k}h_{l}\mathsf{trace}(\mathbb{E}[ \mathbf{S}_{kl}])$ (66) $$

where the term in equation 64 cancels out with the second term in equation 63,
the term in equation 65 contains cross-products including two error matrices $\mathbf{E}_{r}$ with the same index $r$, and $\mathbf{S}_{kl}$ aggregates the sum of quadratic forms with two error matrices with different indices (i.e., terms of the form
$f_{1}(\mathbf{\hat{C}})\mathbf{E}_{r_{1}}f_{2}(\mathbf{\hat{C}})\mathbf{E}_{r_
{1}}f_{3}(\mathbf{\hat{C}})\mathbf{E}_{r_{2}}f_{4}(\mathbf{\hat{C}})\mathbf{E}
_{r_{2}}f_{5}(\mathbf{\hat{C}})$ or $f_{1}(\mathbf{\hat{C}})\mathbf{E}_{r_{1}}f_{2}(\mathbf{\hat{C}})\mathbf{E}_{r_
{1}}f_{3}(\mathbf{\hat{C}})\mathbf{E}_{r_{2}}f_{4}(\mathbf{\hat{C}})$ for two error matrices $\mathbf{E}_{r_{1}}$ and $\mathbf{E}_{r_{2}}$, $r_{1}\neq r_{2}$). We now proceed to analyze the three terms in equation 64,equation 65,equation 66.

First term.
The term in equation 64 is the opposite of the second term in equation 63, so it cancels out when substituted.

Second term.
Analogously to (Gao et al., , eq. (B.8)), we rewrite the term in equation 65 leveraging the linearity of trace and expectation, the trace cyclic property $\mathsf{trace}(\mathbf{A}\mathbf{B}\mathbf{C})=\mathsf{trace}(\mathbf{C}
\mathbf{A}\mathbf{B})=\mathsf{trace}(\mathbf{B}\mathbf{C}\mathbf{A})$ and rearranging the terms to change the sum limits. Then we exploit Lemma 4 to upper bound it as

$$ $\displaystyle\mathbb{E}\left[\sum_{r=1}^{K}\mathsf{trace}\left(\sum_{k=r}^{K} \sum_{l=r}^{K}h_{k}h_{l}\mathbf{E}_{r}\mathbf{\hat{C}}^{k+l-2r}\mathbf{E}_{r} \mathbf{\hat{C}}^{r-1}\mathbf{x}\mathbf{x}^{\mathsf{T}}\mathbf{\hat{C}}^{r-1} \right)\right]\leq NP^{2}\|\mathbf{x}\|^{2}_{2}Q$ (67) $$

where $Q$ is defined in Lemma 4 and is a sum of terms linear in the probability value.

Third term.
The term in equation 66 can be bounded by an expression similar to equation 44, but with at least two of the terms among $\{\mathbb{E}[\mathbf{E}_{r_{1}}^{2}],\mathbb{E}[\mathbf{E}_{r_{2}}^{2}],
\mathbb{E}[\mathbf{E}_{r_{1}}],\mathbb{E}[\mathbf{E}_{r_{2}}]\}$ within the trace operator. Using Lemmas 5, 6 and 7, we know that these trace terms are of the order $\mathcal{O}((1-p_{1})(1-p_{2}))$ for two generic probability values $p_{1},p_{2}$.
Since the frequency response of the covariance filter $h(\bm{\lambda})$ is bounded, and consequently the coefficients $h_{k}$ are also bounded, we have that

$$ $\displaystyle\sum_{k=0}^{K}\sum_{l=0}^{K}h_{k}h_{l}\mathsf{trace}(\mathbb{E}[ \mathbf{S}_{kl}])=\mathcal{O}((1-p_{1})(1-p_{2})).$ (68) $$

By substituting the three bounds into equation 62, noticing that the terms of quadratic order $\mathcal{O}((1-p_{1})(1-p_{2}))$ are dominated by the linear terms in $Q$ and using the fact that $\|\mathbf{x}\|\leq 1$, we obtain the final bound:

$$ $\displaystyle\mathbb{E}[\|\mathbf{u}-\mathbf{\tilde{u}}\|^{2}]\leq NP^{2}Q+ \mathcal{O}((1-p_{1})(1-p_{2})).$ (69) $$

∎

###### Lemma 3 .

###### Proof.

###### Lemma 4 .

###### Proof.

###### Lemma 5 .

###### Proof.

###### Lemma 6 .

###### Proof.

###### Lemma 7 .

###### Proof.

## Appendix G Sufficient Condition for PSD Thresholded Matrix

Consider a hard-thresholded sample covariance matrix $\mathbf{\bar{C}}$ as per Def. 2. Given the sample covariance matrix $\mathbf{\hat{C}}$, the following is a sufficient condition for $\mathbf{\bar{C}}$ to be positive semidefinite (PSD) (Bickel & Levina, , Section 2):
$\mathbf{\bar{C}}$ is PSD if $\|\mathbf{\bar{C}}-\mathbf{\hat{C}}\|\leq\epsilon$ and $\lambda_{\textnormal{min}}(\mathbf{\hat{C}})>\epsilon$ for an $\epsilon>0$, where $\lambda_{\textnormal{min}}(\cdot)$ computes the smallest eigenvalue. That is, if the true covariance matrix is sparse, then the term $\|\mathbf{\bar{C}}-\mathbf{\hat{C}}\|$ decreases as the number of samples $t$ increases and the sparse estimate is more likely to be PSD.

## Appendix H Extension of Covariance Sparsification Techniques to Graphs

The proposed sparsification techniques can also be applied to weighted graphs.
Hard and soft thresholding of a graph consists of removing all the edges whose absolute weight is below a given threshold, and decreasing the weight of the others in case of soft thresholding.
Stochastic sparsification can also be performed on graphs, where each edge $a_{ij}$ (where $a_{ij}$ is the absolute weight of the edge between nodes $i,j$ or 0 if the edge does not exist) is removed with a probability $1-p_{ij}$. In the case of ACV, $p_{ij}=|a_{ij}|/\max_{i,j}|a_{ij}|$.
For RCV, instead, we can control the level of sparsification through a parameter $p$ by defining an ordered set of probabilities $\mathcal{P}=\{p^{\prime}_{1},\dots,p^{\prime}_{N^{\prime}}\}$, where $p^{\prime}_{i}\leq p^{\prime}_{i+1}$, and each $p^{\prime}_{i}$ is sampled from $\mathcal{N}(p,\sigma)$ with $\sigma=\min((1-p)/3,p/3)$ (such that the number of values not in $[0,1]$, which we clip to the interval, is negligible) and $N^{\prime}$ is the number of probability values for assignment. We then set $p_{ij}=p^{\prime}_{k}$ where $k=|\{a_{lm}:|a_{lm}|<|a_{ij}|;l,m=1,\dots,N\}|$ is the position of $a_{ij}$ in the ordered ranking of absolute edge weights.
This allows to increase the efficiency of GNNs, which is a highly explored research direction (see, for example, Ye & Ji (); Peng et al. (); Rong et al. ()).

However, we here explicitly and purposely focus on VNNs and on covariance sparsification in order to, at the same time, improve VNN time efficiency and preserve the theoretical tractability of their stability with respect to covariance estimation errors, which improves with our techniques under the assumption of a sparse true covariance.

## Appendix I Experimental Setup and Additional Results

**Table 1: Details on datasets.**
| Dataset | Nodes | Node features | Classes | Samples (train/valid/test) |
| --- | --- | --- | --- | --- |
| Epilepsy | 76 | 1 | 2 | 8000 (6400/800/800) |
| CNI | 200 | 122 | 2 | 240 (160/40/40) |
| MHEALTH | 24 | 128 | 12 | 5229 (3162/1003/1064) |
| Realdisp | 117 | 128 | 10 | 5831 (2651/1296/1884) |

### I.1 Datasets

We provide additional details about the four real datasets we use in the experiments. Table 1 summarizes the main information about the datasets and Fig. 7 shows their empirical normalized covariance values distribution.

- •
Epilepsy (Kramer et al., ) consists of an electrocorticography
time series collected during an epilepsy study at the University
of California, San Francisco Epilepsy Center (Kramer et al., ).
The dataset contains 8000 time series samples recorded at 76 locations in the brain, equally split before and after an epileptic seizure event. We use the same data preprocessing technique in (Natali et al., ), and we further z-normalize the samples of each class separately to prevent the classification task from being affected by trivial characteristics such as signal amplitude.
We interpret each time series sample as a signal on a graph with 76 nodes, such that the covariance matrix represents brain functional connectivity. We perform binary classification to predict whether the sample was recorded before or after the seizure. We divide the datasets into train/validation/test splits of size 80%/10%/10%.
- •
CNI (Schirmer et al., ) was released for the Connectomics in NeuroImaging challenge and contains resting-state fMRI time series of patients with attention deficit hyperactivity disorder (ADHD) and neurotypical controls (NC). Each time series consists of 122 graph signals recorded in 200 brain locations (nodes) per patient. The available training set contains 200 patients and the test set 40 patients. We therefore consider a VNN with input feature size 122 over 200 nodes and we perform classification to predict whether the patient is ADHD or NC.
- •
MHEALTH (Banos et al., ) contains measurements of wearable devices placed in the chest, right wrist and left ankle of 10 subjects performing 12 different actions. Since random data splits introduce contaminations due to overlapping windows (Tello et al., ), we use recordings from different subjects for train, validation and test. Specifically, we use subjects 6,10 for validation, 2,9 for test and the rest for training. We do not use the ECG measurements placed in the chest. The goal is to classify the action performed by the subject. The sampling rate is 50 Hz.
- •
Realdisp (Banos et al., ) contains data of 17 subjects performing 33 actions. We use subjects 4, 6, 10 and 11 for validation; subjects 1, 7, 8, 9, 12 and 14 for test and the remaining for training. We use the 10 most common activity labels: walking, jogging, running, cycling, elliptic bike, trunk twist (arms outstretched), rowing, knees (alternatively) bend forward, waist bends forward, trunk twist (elbows bended). We segment the data in MHEALTH and Realdisp creating sliding windows of size 128 (2.56 seconds) for measurements relative to the same activity with 50% overlap.

### I.2 Experimental setup

Hyperparameters.
Table 2 reports the hyperparameters of VNNs used for experiments on different datasets, which we find through a hyperparameter search.
For regression and classification, we average the outputs of VNN over nodes and apply a 2-layer MLP for the final task.
We use Adam optimizer with weight decay 0.001.
For GraphSAGE in Table 3, we use the covariance matrix as a graph and we use the same hyperparameter configuration as VNN except for Realdisp, for which we use 2 layers of size 32 and 100 epochs.

Implementation details.
We repeat all experiments 5 times and, for stochastic experiments, we sample 10 different sparsified covariance matrixes. We report average results and standard deviations. Experiments are run on a 13th Gen Intel Core i7-1365U CPU. We implement models using Pytorch and for GraphSAGE we use Pytorch Geometric (Fey & Lenssen, ).
The sample covariance matrix $\mathbf{\hat{C}}$ is computed once over the complete training set and fixed during training, i.e., it is not affected by batching. For stability experiments, we use different sparsified covariances during test, whereas for the other experiments we keep the same matrix estimated on the training set.

**Table 2: VNN hyperparameters on different datasets.**
| Dataset | $L$ | $F_{1}\{,\dots,F_{L}\}$ | $K$ | Epochs | Learning rate | Batch size |
| --- | --- | --- | --- | --- | --- | --- |
| Epilepsy | 5 | 32,32,32,32,32 | 1 | 200 | 0.001 | 100 |
| CNI | 1 | 64 | 1 | 50 | 0.015 | 50 |
| MHEALTH | 2 | 32,32 | 1 | 200 | 0.015 | 3162 |
| Realdisp | 2 | 128,64 | 1 | 500 | 0.001 | 5831 |
| SparseCov | 2 | 13,13 | 1 | 50 | 0.015 | 800 |
| LargeCov | 2 | 32,32 | 1 | 50 | 0.015 | 800 |
| SmallCov | 2 | 32,32 | 1 | 50 | 0.015 | 800 |

Figure: Figure 5: Stability of VNN and PCA-SVM with dense and sparse covariance on regression task for a synthetic dataset with sparse true covariance. (Left) Regression performance in terms of Mean Absolute Error. (Right) Embedding difference (i.e., $\|\mathsf{\Phi}(\mathbf{x},\mathbf{\bar{C}},\mathcal{H})-\mathsf{\Phi}(\mathbf {x},\mathbf{C},\mathcal{H})\|$) between VNNs/PCA with true and estimated covariance. Standard deviations in our results are of order at most $10^{-2}$ for VNNs and $10^{-1}$ for PCA. Legend is shared.
Refer to caption: x11.png

**Table 3: Accuracy (%) on real datasets.**
|  |  | Epilepsy | CNI | MHEALTH | Realdisp |
| --- | --- | --- | --- | --- | --- |
| Baselines | PCA + SVM | 100.0$\pm$0.0 | 50.0$\pm$0.0 | 17.7$\pm$0.0 | 11.9$\pm$0.0 |
| Kernel PCA + SVM | 100.0$\pm$0.0 | 55.0$\pm$0.0 | 17.7$\pm$0.0 | 11.4$\pm$0.0 |  |
| Sparse PCA + SVM | 100.0$\pm$0.0 | 55.0$\pm$0.0 | 17.7$\pm$0.0 | 11.8$\pm$0.0 |  |
| GraphSAGE (Hamilton et al., ) | 57.9$\pm$0.7 | 55.0$\pm$0.7 | 77.0$\pm$1.9 | 67.3$\pm$1.2 |  |
| Dense VNN (Sihag et al., ) | 98.6$\pm$0.4 | 54.0$\pm$3.8 | 87.6$\pm$3.6 | 71.6$\pm$4.0 |  |
| Ours | Hard-thr ($p=0.75$) | 98.9$\pm$0.5 | 53.5$\pm$7.6 | 85.9$\pm$4.3 | 69.0$\pm$3.8 |
| Hard-thr ($p=0.5$) | 98.8$\pm$0.6 | 54.0$\pm$2.9 | 88.6$\pm$1.9 | 70.6$\pm$3.0 |  |
| Hard-thr ($p=0.25$) | 97.6$\pm$0.6 | 54.0$\pm$2.9 | 88.5$\pm$4.2 | 67.2$\pm$8.0 |  |
| Soft-thr ($p=0.75$) | 99.2$\pm$0.3 | 55.0$\pm$4.0 | 89.3$\pm$0.9 | 61.9$\pm$4.7 |  |
| Soft-thr ($p=0.5$) | 99.1$\pm$0.4 | 53.5$\pm$3.8 | 87.6$\pm$3.9 | 73.8$\pm$3.9 |  |
| Soft-thr ($p=0.25$) | 98.2$\pm$0.8 | 55.0$\pm$1.8 | 88.5$\pm$4.0 | 65.8$\pm$7.0 |  |
| RCV ($p=0.75$) | 98.0$\pm$1.2 | 56.0$\pm$6.5 | 90.4$\pm$5.3 | 72.3$\pm$2.8 |  |
| RCV ($p=0.5$) | 99.1$\pm$0.3 | 55.5$\pm$4.8 | 91.0$\pm$2.2 | 69.1$\pm$4.6 |  |
| RCV ($p=0.25$) | 99.4$\pm$0.3 | 55.0$\pm$1.8 | 88.1$\pm$3.1 | 71.6$\pm$1.0 |  |
| ACV | 98.8$\pm$1.0 | 57.5$\pm$1.8 | 92.0$\pm$1.8 | 65.0$\pm$4.3 |  |

**Table 4: Time (sec) for a batch forward pass on real datasets.**
|  |  | Epilepsy | CNI | MHEALTH | Realdisp |
| --- | --- | --- | --- | --- | --- |
| Baselines | PCA + SVM | 0.35$\pm$0.02 | 0.03$\pm$0.02 | 0.08$\pm$0.01 | 0.14$\pm$0.01 |
| Kernel PCA + SVM | 1.90$\pm$0.04 | 0.03$\pm$0.02 | 1.37$\pm$0.14 | 0.78$\pm$0.08 |  |
| Sparse PCA + SVM | 0.57$\pm$0.07 | 0.02$\pm$0.00 | 0.08$\pm$0.00 | 0.14$\pm$0.01 |  |
| GraphSAGE (Hamilton et al., ) | 4.47$\pm$0.18 | 0.65$\pm$0.18 | 0.70$\pm$0.08 | 60.6$\pm$9.7 |  |
| Dense VNN (Sihag et al., ) | 4.49$\pm$0.01 | 0.61$\pm$0.08 | 1.58$\pm$0.07 | 72.6$\pm$1.6 |  |
| Ours | Hard-thr ($p=0.75$) | 3.93$\pm$0.17 | 0.53$\pm$0.01 | 1.24$\pm$0.03 | 55.3$\pm$1.0 |
| Hard-thr ($p=0.5$) | 3.26$\pm$0.31 | 0.39$\pm$0.03 | 0.86$\pm$0.02 | 37.0$\pm$0.8 |  |
| Hard-thr ($p=0.25$) | 2.06$\pm$0.35 | 0.22$\pm$0.00 | 0.47$\pm$0.01 | 19.1$\pm$0.5 |  |
| Soft-thr ($p=0.75$) | 3.91$\pm$0.27 | 0.56$\pm$0.05 | 1.22$\pm$0.05 | 56.0$\pm$1.6 |  |
| Soft-thr ($p=0.5$) | 2.65$\pm$0.08 | 0.38$\pm$0.01 | 0.89$\pm$0.03 | 37.8$\pm$1.0 |  |
| Soft-thr ($p=0.25$) | 1.56$\pm$0.02 | 0.22$\pm$0.00 | 0.52$\pm$0.02 | 19.6$\pm$0.6 |  |
| RCV ($p=0.75$) | 3.90$\pm$0.15 | 0.54$\pm$0.01 | 1.29$\pm$0.06 | 60.1$\pm$0.5 |  |
| RCV ($p=0.5$) | 2.68$\pm$0.04 | 0.40$\pm$0.03 | 1.01$\pm$0.10 | 38.9$\pm$0.8 |  |
| RCV ($p=0.25$) | 1.73$\pm$0.01 | 0.23$\pm$0.01 | 0.64$\pm$0.04 | 21.6$\pm$1.0 |  |
| ACV | 1.27$\pm$0.16 | 0.25$\pm$0.02 | 0.55$\pm$0.04 | 12.8$\pm$0.4 |  |

### I.3 Additional results

Sparse true covariance.
Fig. 5 shows the embedding difference for VNNs and PCA variants with hard and soft thresholding. This confirms the observations in Sec. 5.1 relative to the increased stability of VNNs with thresholding w.r.t. PCA in sparse covariance settings since VNNs’ outputs under sparsification are closer to the embeddings with the true covariance compared to PCA.

Stochastic sparsification.
Fig. 6 shows the stability of VNNs under stochastic sparsification with different strategies on the 4 real-world datasets. On CNI, VNNs remain stable for all sparsification strategies: ACV achieves very close performance to dense VNN on both datasets and RCV performs closely for most values of $p$. This is consistent with the observations in Sec. 5.2, which suggest that most values in these covariances are not useful.
On Epilepsy, MHEALTH and Realdisp, however, the sparsified VNNs appear less stable, despite most covariance values being small (cf. Fig. 7). This can be due to the fact that the VNN on Epilepsy has 5 layers, which lowers stability according to Thm. 1, while the task for HAR datasets is more challenging compared to the task for brain datasets (in this case, there are 10 or 12 classes instead of the 2 for brain datasets), which may lead to lower model stability.

Comparison with PCA and GraphSAGE on real datasets.
We report in Tables 3 and 4 the performance and time efficiency of S-VNNs, dense VNN, three different PCA-based classifiers (standard PCA, kernel PCA (Schölkopf et al., ) and PCA with hard-thresholded covariance matrix, all followed by a kernel SVM classifier) and GraphSAGE (Hamilton et al., ) operating on the covariance matrix as a graph as an example of sparsified GNN.
For PCA-based classifiers, since PCA does not handle multiple features per node, for CNI, MHEALTH and Realdisp, we average the node features before classification.
The results for S-VNNs and VNN are the same as Fig. 4.

PCA-based classifiers are fast, especially those with linear PCA, particularly on CNI, MHEALTH and Realdisp since they operate on the average of all node features. However, their performance is inconsistent. They perform comparably to S-VNNs on Epilepsy and CNI, but very poorly on the HAR datasets, as not considering multiple node features hinders their representation capabilities.
This motivates the need for a more sophisticated and reliable approach such as VNNs.
Furthermore, both sparse and non-linear PCA provably suffer from stability issues, reiterating the advantages of VNNs over them.

GraphSAGE, instead, performs generally worse than S-VNNs (and especially poorly on Epilepsy) and it improves only slightly the time efficiency compared to dense VNN, whereas S-VNNs achieve a more consistent speed-up. This corroborates the effectiveness of our S-VNNs in both improving performance due to the removal of spurious correlation and increasing time efficiency.

###### Remark 3 .

While the non-linear PCA and (sparse) VNNs may, at first sight, seem connected, they are fundamentally different since kernel PCA computes the principal directions of the data projected into a reproducing kernel Hilbert space, whereas VNNs manipulate the eigenvectors of the covariance matrix of the samples in the original space. Therefore, we merely use kernel PCA here as a baseline for additional comparison.

Figure: (a) Epilepsy.
Refer to caption: x13.png

Figure: (a) Epilepsy.
Refer to caption: x17.png

###### Remark 3 .

## Peer Reviews

confidence: 3
contribution: 1
initial_rating: 1
presentation: 1
questions: - How can we know in practice if the covariance of the dataset is sparse?
- How does $\nu$ in thm.1 depends on the data distribution? I only see one data distribution, which depends only on $C$, so $\nu$ is only a function of $C$?
- In def. 1, is it for a fixed matrix? In which set do the eigenvalue pairs $\lambda_i, \lambda_j$ belong? This is unclear.
- Is lemma 1 an original result? If not then a citation is needed here.
- Should $F_{in} = F_{out}$ in (1)? since it seems like $u^l$ is both of size $F_{out}$ and $F_{in}$.
soundness: 2
strengths: - The idea of sparsifying the covariance before feeding it to the network is sound.
summary: This paper builds upon covariance neural networks, which are constructed to process covariance matrices. The authors study the impact of sparsifying the covariance matrix with hard or soft thresholding before feeding it to the network. They demonstrate that this improves the estimation when the true covariance matrix is sparse, and explain how to drop coefficients at random when it is not the case. The authors demonstrate some improvements compared to no thresholding in several experiments.
weaknesses: - This paper is very hard to follow. Several variables are not defined ($V$, $u$, $u_g$, $h_{klfg}$, etc...). Several concepts must be guessed from the text (what are the covariance filters? how does $u$ relate to the $u_g$ in eq. 1? what are the per covariance filters in theorem 1?). Overall, it is hard to understand this paper alone without reading the original paper on covariance neural networks. This paper needs a major rewriting before being ready for publication, which explains my note.
- The novelty is minor. This paper builds upon covariance networks, which are seldom used in practice, and incrementally improves on them by pre-processing the covariances.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: Please address each of the weaknesses above. In particular:
- In Theorem 1 and related results, is it true that $\mathcal{O}$ hides dependence on $\mathcal{H}$, which as in the case of PCA, could be arbitrarily bad?
- What does with high probability mean exactly?
- How are matrices ensured to be PSD?
- What happens in the case of undefined / infinite variance? How should this be reflected in the theorems? Is a Gaussian assumption on x used throughout, or only in the first result?
- Does the theoretical result really show a guaranteed improvement in stability? Or does it only fail to show a decrease in stability?
- Please clarify whether $x$ is drawn from an (unconditional) Gaussian or from a "truncated" Gaussian conditioned on its norm.
soundness: 2
strengths: **Strengths:**
- The proposed architectures come with apparent theoretical guarantees on the closeness of hidden representations or predictions when applied to the approximate versus true covariance matrices. These distances become smaller as the number of samples increase, at a rate of $\mathcal{O}(t^{-1/2})$, ignoring some parameter-dpendent constants.
- Experiments are provided showing the predictive performance and stability of the proposed methods. This experiments agree with the theory, and nicely demonstrate the approach. 
- A very large amount of related literature is cited, and this seems appropriate.
- The theory appears to be mostly sound, roughly based on a Lipschitz assumption. (although there are some isolated places where precision could be improved, see weaknesses below).
summary: **Summary:**
Graph convolutions on the covariance matrix of tabular data are performed by a architecture called Covariance Neural Networks (VNNs). These rely an empirical estimates of the covariance matrix, which is notoriously difficult. The authors propose Sparse VNNs, which sparsifies the covariance matrix before convolution is applied. This is advantageous when the true covariance is sparse, and also when the true covariance is dense, and different techniques are investigated. S-VNNs are compared against VNNs and PCA. Experiments are provided on brain data and human action reconition, showing performance, sensitivity and efficiency benefits over VNNs.
weaknesses: **Weaknesses:**
- Comparing Lemma 1 and Thoerem 1. I am guessing (please correct me if I am wrong!) that somehow hidden inside $\mathcal{O}$ in Theorem 1 are the parameters $\mathcal{H}$. Intuitively, these should learn something similar to PCA, and if the eigenvalues of PCA are close, this constant term in $\mathcal{O}$ will be bad but in terms of $\mathcal{H}$. In Lemma 1, the constant term in terms of the small gap eigenvalues is explicitly given, and it is obvious how this causes instability. Whether my guess about Theorem 1 is correct or not, the authors should state either way about the possibility of dependent of $\mathcal{H}$ on the factors in $\mathcal{O}$. Right now the claim that " VNNs do not suffer from this as the covariance filter can exhibit a stable response to close eigenvalues at the expense of lower discriminability" is not clear, as are similar repeated claims throughout the paper. It seems as though this is also an important consideration for later derived theoretical results (Theorem 2, Proposition 1, ...). **This is my most major important concern.**
- Theorem 1 explicitly writes the probability of the event. The later results only state with high probability. What does with high probability mean?
- In definition 3, how does one ensure that the resulting soft-thresholded matrix is PSD? 
- In definition 4, how does one ensure that the matrix is PSD after dropout? 
- What happens when the true data distribution is heavy tailed, and does not have a (finite) covariance matrix? I guess there should be some condition in the theoretical results, which is currently absent, on the distribution of $x$. In Theorem 1 a Gaussian distribution is used (which has finite variance), and it is not clear in the presentation of later results if a Gausssian is also being used. 
- It is claimed that "For sparse covariance, ... the hard thresholding improves stability." According to my understanding, the authors are able to derive bounds which are tighter for sparse covariances than for dense covariances with hard thresholding. This does not necessarily show that hard thresholding definitively improves stability, only the bound is tighter. If I have understood correctly, perhaps the authors should rephrase their claim (otherwise, am happy to be corrected).



**Minor:**
- I don't understand one sentence in Theorem 1. "Consider a generic data sample $x \sim \mathcal{N} (0, C)$ such that $\Vert x \Vert \leq 1$". Does this mean $x$ is drawn from the conditional distribution obtained by conditioning a Gaussian random vector on the event that the norm is less than or equal to 1? I believe this should be phrased better, in terms of the distribution x is actually drawn from. (I understand this is from another paper, but still it would be nice to improve clarity here).
- Incorrect grammar for paragraph starting line 167/168.
- Limitations of PCA in terms of unstable or poor estimation of eigenvalues is discussed. Do the authors know where does this fit into more robust variants of PCA, exponential family PCA, kernel PCA, etc.? Do such advanced variants of PCA overcome these issues? Either way, it would be nice to mention in the discussion in a sentence around line 39/40.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: 1. See weakness. 

2. Additionally, using t to represent the sample size is somewhat unusual to me, though this is merely a notation issue.
soundness: 2
strengths: The presentation of the paper is clear.
summary: This paper presents a sparse version of covariance neural networks (VNNs), showing that the sparse-VNNs are more stable than nomial VNNs and sparse component analysis.
weaknesses: 1. It is well known that the covariance matrix is not invariant to the scale of the data, making it impractical to set a common threshold for different elements of c_{ij}. Thus, the rationale behind Definition 2 is difficult to justify.

2. From a sparsity perspective, the elements c_{ij}'s should be classified into two categories: zero and nonzero. However, this key point is obscured in Theorem 2, making the results difficult to interpret. In other words, it is challenging to justify that the proposed method effectively denoises the data.

3. The contribution of the paper appears incremental in light of the existing work by Sihag et al. (2022).

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: Although current results can validate the effectiveness of the proposed strategies, the results should be compared with the covariance with sparsity regularizer, not just dense-VNN and robust-PCA.
soundness: 3
strengths: Quality: The paper is well-written, the motivation sounds reasonable, and the proposed sparsification strategies seem good with a theoretical analysis of the stability.
Originality: The paper proposes several sparsification strategies to tackle the spurious correlation and computation cost issues for the VNN. Although the theoretical analysis of the stability of VNN is good, the originality of the proposed solutions is limited. 
Significance: Spurious correlation and sparsity techniques are important.
summary: This paper studies the sparsification strategy upon the spurious correlation and computation efficiency issues for the covariance matrix used in the coVariance Neural Networks (VNN), which is more stable than the Principal component analysis (PCA) methods in the previous studies. This paper proposes hard and soft thresholding strategies if the true covariance matrix is sparse, and two stochastic sparsification techniques, including Absulute covariance values (ACV) and Ranked covariance values (RCV) when the covariance matrix is dense, and theoretically analyze the sparsification error and covariance uncertainty for the stability of VNN. The effectiveness of sparsity techniques is validated through synthetic and real datasets in contrast to the dense-VNN and sparse-PCA on the performance and computational time.
weaknesses: The novelty of the proposed strategies is limited, similar strategies have been used in the study of neural networks like dropout or pruning, and the current results are not enough strong.

Meta-review: The paper proposes Sparse coVariance Neural Networks (S-VNNs), built upon previous work on VNNs, which applies sparsification techniques on the sample covariance matrix before convolution. When the true covariance matrix is sparse, hard and soft thresholdings are proposed to improve covariance estimation and reduce computational cost. When the true covariance is dense, stochastic sparsification is proposed, where data correlations are dropped in probability according to principled strategies. The proposed method is accompanied by numerical experiments.

While reviewers appreciate that the idea of sparsifying the covariance matrix before convolution is sound, they also have concerns about the limited novelty compared to prior work, the weakness of the theoretical results, and the quality of presentation of the current manuscript (see further below).

Decision: Reject