AD+},{\mathbf{v}}_{i}\\!>\\!|$ for $i\in\\{1,\dots,50\\}$, where
${\mathbf{v}}_{i}$ is the $i$-th eigenvector of covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ associated with its $i$-largest eigenvalue.
The bars are evaluated from the mean of $|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{i}\\!>\\!|$ obtained for individuals in the AD+ group
(results for eigenvectors associated with coefficient of variation of
$|\\!<\\!\bar{\mathbf{r}}_{\sf AD+},{\mathbf{v}}_{i}\\!>\\!|$ larger than
$30\%$ across the AD+ group excluded). For every individual in AD+ group, the
associations of its regional residuals with eigenvectors of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ were evaluated as the mean of those over
$100$ nominal VNN models (trained on the OASIS-3 dataset). The eigenvectors
associated with top three largest values for $|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{i}\\!>\\!|$ are plotted on the brain template in Panel b.
Subcallosal region in the right hemisphere was associated with the element
with the largest magnitude in ${\mathbf{v}}_{3}$ and ${\mathbf{v}}_{4}$ and is
highlighted with a red circle in the corresponding plots. Panel c displays the
boxplots for the distributions of $|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{i}\\!>\\!|$ for ${\mathbf{v}}_{1}$, ${\mathbf{v}}_{3}$, and
${\mathbf{v}}_{4}$ across the AD+ group.

####  5.4.2 Regional residuals derived from VNNs transferred from FTDC100 to
OASIS-3 are correlated with principal components of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$

For this set of experiments, we evaluated the regional residuals for all
individuals in the OASIS-3 dataset using VNNs that were trained on FTDC100
dataset and transferred to process OASIS-3. Figure 9 displays the bar plot for
the mean of the inner product between the regional residuals for AD+ group and
first $30$ eigenvectors of $\hat{\mathbf{C}}_{148}^{\sf AD+}$. Similar to Fig.
8a, the regional residuals derived from VNNs that were transferred from
FTDC100 to OASIS-3 had the largest correlations with the first, third, and
fourth eigenvectors of the anatomical covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ (${|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{4}\\!>\\!|=0.56\pm 0.013}$, ${|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{3}\\!>\\!|=0.362\pm 0.014}$,
${|\\!<\\!\bar{\mathbf{r}}_{\sf AD+},{\mathbf{v}}_{1}\\!>\\!|=0.29\pm 0.004}$
). These observations implied that the VNN models trained on FTDC100 dataset
were able to leverage the relevant eigenvectors for brain age prediction in
OASIS-3 dataset in the formation of regional residuals. From the findings in
Fig. 8a and Fig. 9, we can conclude that the VNN models trained to predict
chronological age in both OASIS-3 and FTDC100 were able to leverage the
principal components of $\hat{\mathbf{C}}_{148}^{\sf AD+}$ that were
associated with the regional profile of elevated $\Delta$-Age in Fig. 5.3.1.
The first, third, and fourth eigenvectors of $\hat{\mathbf{C}}_{148}^{\sf
AD+}$ were also the three most significantly associated with the regional
residuals derived from VNNs that were transferred from FTDC300 or FTDC500
datasets to OASIS-3 dataset. The results in this context are included in Fig.
19 in Appendix J.

Figure 9: Mean inner product between the normalized vector of regional
residuals (norm = 1) of VNN outputs (transferred from FTDC100 to OASIS-3)
obtained from AD+ group and the eigenvectors of $\hat{\mathbf{C}}_{148}^{\sf
AD+}$ (covariance matrix of combined HC and AD+ group). This figure
illustrates a bar plot for $|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{i}\\!>\\!|$ for $i\in\\{1,\dots,30\\}$, where
${\mathbf{v}}_{i}$ is the $i$-th eigenvector of covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ associated with its $i$-largest eigenvalue.
The bars are evaluated from the mean of $|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{i}\\!>\\!|$ obtained for individuals in the AD+ group
(results for eigenvectors associated with coefficient of variation of
$|\\!<\\!\bar{\mathbf{r}}_{\sf AD+},{\mathbf{v}}_{i}\\!>\\!|$ larger than
$30\%$ excluded). For every individual in AD+ group, the association of its
regional residuals with eigenvectors of $\hat{\mathbf{C}}_{148}^{\sf AD+}$
were evaluated over $100$ nominal VNN models (trained on the FTDC100 dataset).
Figure 10: Regional profiles derived from the robustness of regional residuals
being elevated for AD+ group with respect to HC group in OASIS-3 using
different sets of VNNs. Panel a projects the robustness of the regional
residuals being elevated for AD+ with respect to HC group across the 100 VNN
models on a brain template, where the VNNs were trained on the OASIS-3
dataset. The number of VNN models for which a brain region was deemed to have
a significantly higher regional residual for AD+ group with respect to HC
group quantifies the robustness of that brain region as a contributor to the
observed elevated $\Delta$-Age in AD+ with respect to HC. Panel b displays the
regional profiles obtained from the analysis of regional residuals for AD+ and
HC groups, where the regional residuals were derived from $100$ VNNs trained
to predict chronological age for HC group in FTDC100. Panel c displays the
corresponding results derived from $100$ VNNs that were randomly initialized.

####  5.4.3 Regional profiles for brain regions with elevated regional
residuals in AD+ group are qualitatively preserved when VNNs are transferred
from FTDC100 to OASIS-3

Next, we checked whether the consistency in correlations with principal
components in Fig. 8a and Fig. 9 resulted in qualitative consistency in the
regional profiles derived from the VNNs trained on FTDC100 dataset for regions
associated with elevated regional residuals in AD+ group with respect to HC
group in OASIS-3. Figure 10b illustrates the projection of the robustness of
the observed significantly higher regional residuals for AD+ group with
respect to HC group on a brain template, where the robustness was evaluated
using $100$ nominal VNN models in a similar fashion as for the results in Fig.
5.3.1b. For clarity, the regional profile for elevated $\Delta$-Age obtained
from the VNNs that were trained on OASIS-3 dataset have been included in Fig.
10a.

The regional profiles obtained using VNNs trained on FTDC100 in Fig. 10b were
largely consistent with that in Fig. 10a. This observation implied that the
VNNs trained on FTDC100 dataset possessed the ability to extract regional
profiles characteristic of accelerated aging in OASIS-3 dataset. This
observation could be attributed to the fact that the VNNs trained on FTDC100
were able to exploit similar principal components of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ as the VNNs that were trained on OASIS-3
dataset.

By training a VNN to predict chronological age in healthy individuals, we aim
to gain the information about healthy aging, which is expected to be
instrumental in detecting patterns of accelerated aging. Our results thus far
have shown that the markers of accelerated aging identified by VNNs trained in
this manner can be linked to the principal components of the anatomical
covariance matrix irrespective of quantitative transferability. Thus, we
conjecture that training the VNNs to predict chronological age helps fine tune
their ability to manipulate the input data using principal components that
were relevant to the elevated $\Delta$-Age effect in AD+ group. To further
support this conjecture, we evaluated the regional profiles for brain regions
that exhibited elevated regional residuals for AD+ group with respect to HC
group in OASIS-3 using 100 VNNs that were randomly initialized (i.e., not
trained whatsoever) and had the same architecture as the VNNs that were
trained on OASIS-3. Figure 10c demonstrates the regional profiles derived from
randomly initialized VNNs on OASIS-3 dataset. Figure 10c shows the robustness
of the elevated regional residuals for AD+ group with respect to HC group was
severely diminished with respect to the parallel results in Fig. 10a and Fig.
10b.

Our experiments here have demonstrated that the ability of VNNs to extract
information about accelerated aging from cortical thickness data relies on
their ability to exploit the eigenvectors or principal components of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ and not necessarily the performance on the
task of predicting chronological age of healthy individuals. Thus, training
the VNNs to predict chronological age in HC group helped fine tune their
parameters to exploit the eigenvectors or principal components of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ relevant to $\Delta$-Age in AD+ group. The
observations above facilitate the decoupling of the task of brain age
prediction from the objective of achieving near perfect performance on the
task of chronological age prediction in the HC group.

Additional Experiments: The associations between regional residuals and
eigenvectors of $\hat{\mathbf{C}}_{148}^{\sf AD+}$ and the regional profiles
for brain regions with elevated regional residuals in AD+ group when VNNs were
transferred from FTDC300 or FTDC500 are included in Fig. 19 in Appendix J. The
results in Appendix J further demonstrate that the most significant
correlations observed in Fig. 6a for VNNs trained on OASIS-3 were largely
preserved when VNNs were transferred from FTDC300 or FTDC500 datasets to
OASIS-3. Additional experiments in Appendix L help establish the stability of
the regional profiles in Fig. 5.3.1b to perturbations in the covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ and demonstrate that the VNN models used
here were not overfit on the composition of the population used to estimate
the anatomical covariance matrix. In Appendix M, we report that improving the
performance of VNNs on the chronological age prediction task using an adaptive
readout may penalize the regional interpretability offered by VNNs in the
context of brain age (besides sacrificing the scale-free characteristic and
the transferability guarantees for VNNs).

##  6 Discussion

Graph convolution operator on covariance matrix, termed as a coVariance
filter, forms the backbone of VNN architecture. The coefficients of the
coVariance filter characterize its ability to manipulate the data according to
the eigenspectrum of the covariance matrix to achieve a learning objective.
Thus, statistical inference using VNNs draws similarities with PCA-driven
statistical approaches. However, PCA conventionally operates within the
feature space of a given dataset and hence, does not provide any notion of
similarity between principal components of datasets with different number of
features. In this paper, we have studied the key property of transferability
of VNN models, which allows VNNs to be transferable between datasets with
similar characteristics but different number of features. The notion of
similarity between datasets consisting of different number of features is
borrowed from the existing theory of graphons that studies limits of dense
graphs [81]. Specifically, our theoretical results show that if there exists a
sequence of covariance matrices that converges to a continuous limit object in
the asymptote of infinite number of features, the VNNs can be transferred
between any two covariance matrices of such a sequence for statistical
inference. The underlying theoretical results rely on the convergence of the
eigenspectrum of a continuous approximation of covariance matrices, which
result in convergence of the coVariance filter outputs for covariance matrices
belonging to a converging sequence, and subsequently, the convergence of VNN
outputs. Our experiments pertain to dense anatomical covariance matrices and
therefore, graphon model-based analyses were certainly appropriate to study
transferability of VNNs. We also note that the parallel results on
transferability of GNNs in [79] are restricted to graphs with no self-loops or
unlabeled graphs, and these results can be recovered from Theorem 3 when we
have $\zeta=1$ (i.e., all features in the dataset have equal variance).
Furthermore, we remark that sparse covariance matrices are also of practical
interest as they can help manage computational complexity [104]. Therefore,
studying VNN transferability over sparse covariance matrices is a future
direction of interest.

We also observed that functions of VNN outputs were correlated with several
eigenvectors of the underlying covariance matrix (Fig. 8), thus, validating
our hypothesis based on Lemma 1 that VNNs manipulate the data according to the
eigenspectrum of the covariance matrix. Hence, information processing with
VNNs draws similarities with PCA based analysis. Unlike PCA, VNNs are
guaranteed to be stable to perturbations in the underlying covariance matrix
(Theorem 2). Previously, we have empirically validated the stability of
regression performance of VNNs to perturbations in the covariance matrix in
[1]. Figure M demonstrates that the stability of VNNs to perturbations in the
covariance matrix also extended to the observations relevant to gauging the
interpretability of VNN outputs in the context of brain age.

From a data analytics perspective, we empirically showed that VNNs were able
to extract information about chronological age from cortical thickness
features of healthy controls and this information was transferable with
minimal difference across cortical thickness datasets curated according to
different scales of Schaefer’s atlas. These observations were predicted by our
theoretical results as the covariance matrices derived from cortical thickness
features at different scales of Schaefer’s atlas formed a converging sequence.
Although we used Schaefer’s atlas to demonstrate the transferability of VNNs
within the scope of our theoretical results, our findings are not restricted
to Schaefer’s atlas and can readily be extended to any other brain atlases
that accommodate multi-scale parcellations. We also note that this
transferability of performance seemingly came at the cost of performance.
Specifically, in Appendix M, we report that the regression performance of VNNs
could be improved significantly using adaptive readouts. Use of adaptive
readouts sacrifices end-to-end transferability and scale-free aspect of VNNs
if the number of learnable features in the readout function is dependent on
the dimensionality of the data. Therefore, further study of VNNs with adaptive
readouts may be necessary for performance-oriented applications.

Despite growing interest in similar multi-scale datasets in neuroscience [51,
52, 53], the analyses strategies that optimally identify or leverage the
redundancy in such datasets across different scales are currently lacking. In
this context, VNNs also provide a novel analysis framework that can readily
transfer across multiple scales. Analyses on datasets curated according to
different brain atlases also revealed the limits of transferability of VNNs
(Table 1 and Table 2). For instance, the VNNs trained on cortical thickness
features from FTDC100 datasets (Schaefer’s atlas) did not retain regression
performance on OASIS-3 dataset. These observations imply that VNNs may need to
be augmented with mappings between brain atlases to optimally generalize their
performance on chronological age prediction beyond the multi-scale cortical
thickness datasets.

In our experiments, the VNNs were trained on datasets extracted from a
population of older adults (age $>$ 50 years). However, the changes in the
brain are not uniform across the anatomical regions for different age groups -
some regions grow more in development [105, 106] and atrophy more in aging
than other regions [107, 108]. Thus, factors related to the impacts of
development and neurodegeneration on the brain may render the VNNs used in
this paper inadequate for capturing aging related information from cortical
thickness data for a population with a broader or a different age profile than
the datasets studied here. From a technical perspective, such biological
factors may lead to violation of the assumptions necessary for the theoretical
guarantees on the transferability of VNNs when assessing transference between
neuroimaging datasets consisting of healthy controls from different age
groups. Moreover, similar considerations may be necessary for longitudinal
analyses where the VNNs trained on the data from a specific age group may not
be equipped to perform accurate inference at all stages of the lifespan of an
individual.

While our focus in this paper has been strictly on cortical thickness datasets
that were derived from structural MRI, there exist multiple neuroimaging
modalities, such as diffusion MRI, functional MRI, and electroencephalogram
(EEG) that provide distinct insights into brain activity and structure. These
modalities typically have smaller signal-to-noise ratio (SNR) as compared to
structural MRI [109, 110]. Hence, it is imperative to investigate the
transferability of VNNs over different modalities of neuroimaging in order to
gauge widespread applicability of their transferability across multi-scale
datasets.

Brain age prediction task provides a unique statistical challenge as the brain
age has no ground truth and the machine learning models used must be able to
capture the changes driven by age-related neurodegeneration that lead to
accelerated aging. In this context, our study of VNNs provides a foundational
contribution to the methodology of brain age prediction. Our results showed
that besides facilitating transferability, the non-adaptive readout function
was instrumental to characterizing the regional interpretability of VNN models
in the brain age prediction task. Specifically, the final regression output
could be written as a mean of entities of an $m$-dimensional vector (as in
(32)) and therefore, allowed us to evaluate the contribution of each brain
region to the final regression output. Thus, VNN models could also provide a
regional profile if an elevated brain age was observed in morbidity.
Furthermore, the regional profiles extracted by VNN models were correlated
with certain principal components of the anatomical covariance matrix.
Importantly, training the VNNs to predict chronological age helped fine tune
their parameters to exploit the relevant principal components of the
anatomical covariance matrix. Thus, the role of the age-bias correction step
was restricted to projecting the VNN outputs onto a space where one could
observe biological aging with respect to the chronological age from a layman’s
perspective.

Interestingly, the VNNs transferred from FTDC datasets to OASIS-3 datasets
were able to exploit the same brain age-related principal components of the
anatomical covariance matrix of OASIS-3 dataset as the VNNs trained on OASIS-3
dataset, which resulted in consistent regional profiles for the brain regions
with elevated regional residuals in the AD+ group for the two sets of VNNs.
This observation further implied that the regional profiles identified by VNNs
were robust to various factors characterizing the heterogeneity across FTDC
and OASIS-3 datasets (such as distinct quality of neuroimaging data, contrasts
across scanners in MRI acquisition, and different cohort compositions). Thus,
our findings suggest that the convolution operation modeled by coVariance
filters of VNNs provides a useful analytic tool to derive interpretable,
spatially robust, and reproducible information from cortical thickness
datasets that is relevant for brain age prediction.

The results in Appendix M showed that improving the performance on the
chronological age prediction task may not necessarily improve brain age
prediction, either in terms of higher $\Delta$-Age or better characterization
of brain regions perceived to be the contributors to elevated $\Delta$-Age.
Brain age is a coarse metric that is expected to be elevated as compared to
chronological age in various neurodegenerative conditions but may not have
enough discriminability to discriminate between them. Therefore, by
associating $\Delta$-Age with a regional profile, VNNs also provide a feasible
tool to distinguish pathologies if the distributions of $\Delta$-Age for them
are overlapping.

Currently, there is no clear benchmark for an acceptable performance of a
machine learning model on the chronological age prediction task in order for
it to predict brain age as most existing approaches in this domain lack
interpretability. For instance, the study in [18] recognizes brain age
prediction using ‘loosely fitted’ models on chronological age as a potential
pitfall in the analyses while another reports better brain age prediction
using a ‘moderately’ fitted deep learning model [16]. The study in [15]
associates voxel-wise interpretability to brain age but the link between the
interpretability patterns and accuracy of chronological age prediction by
their approach is not discussed. One motivation to perform brain age
prediction with models that have high quality fit on the chronological age is
presumably to mitigate the impact of age-bias correction step on the final
result [18]. This criticism appears to be justified when brain age estimation
approaches cannot isolate the abnormalities that lead them to predict elevated
$\Delta$-Age. However, the simplicity of VNN models allows us to analyze the
deviations in the intermediate steps of brain age estimation before age-bias
correction is applied. In this context, we note that VNNs with unweighted mean
as a readout function could not achieve perfect prediction of chronological
age for healthy controls and yet, their brain age predictions in individuals
with AD were associated with robust, transferable regional profiles. Thus,
VNNs seem to be methodologically adept at finding neurodegeneration-driven
factors that contribute to elevated $\Delta$-Age. The insights provided above
in the context of VNNs and brain age are simply infeasible for brain age
prediction approaches if they rely on complex and non-transparent deep
learning models despite having millions of learnable parameters.

Based on the discussion thus far, we conclude that learning to predict
chronological age in healthy controls is instrumental for VNNs to provide
interepretability to elevated $\Delta$-Age. However, a near-prefect
chronological age prediction for healthy controls by itself may not be a
determinant of the quality of brain age prediction in neurodegeneration. From
a broader perspective, brain age prediction even for healthy controls is a
complex task due to various factors that can contribute to accelerated aging
in the absence of an adverse health condition [111, 112, 113]. While we do not
claim that the VNNs provide the ‘best’ brain age prediction on any metric, our
experiments have convincingly demonstrated that VNNs are able to extract the
sufficient, robust information from cortical thickness datasets for
anatomically interpretable and justified brain age prediction in
neurodegeneration. It is possible that VNN based brain age predictions and
associated regional interpretations could further be optimized or improved
upon in some manner. For instance, incorporating the metrics of aging from DNA
methylation aging [114] in the training of VNNs is a promising future
direction that can help expand our understanding of aging.

Existing studies, including this paper, fall short at concretely defining the
notion of optimal brain age prediction. However, we note that fine-tuning the
fit of machine learning models on chronological age in order to gain the
desired $\Delta$-Age in neurodegeneration or correlations of $\Delta$-Age with
auxiliary measures potentially makes such models overfit on the $\Delta$-Age
itself. Hence, a larger focus is needed on principled statistical approaches
for brain age prediction that can capture the factors that lead to accelerated
aging. Locally interpretable and theoretically grounded deep learning models
such as VNNs can provide a feasible, promising future direction to build
statistically and conceptually legitimate brain age prediction models in
broader contexts.

##  7 Data and Code Availability

MRI and clinical data for individuals in the FTDC datasets may be requested
through https://www.pennbindlab.com/data-sharing and upon review by the
University of Pennsylvania Neurodegenerative Data Sharing Committee, access
will be granted upon reasonable request. OASIS-3 dataset is publicly available
and hosted on central.xnat.org. Code for demonstrating transferability of VNNs
and brain age evaluation is available at
https://github.com/pennbindlab/VNN_Brain_Age. Requests for details regarding
IDs of individuals in OASIS-3 and source data for all figures may be sent to
<EMAIL_ADDRESS>
##  8 Acknowledgements

The MRI data for FTDC datasets were provided by the Penn Frontotemporal
Degeneration Center (NIH AG066597) and Penn Institute on Aging. Cortical
thickness data were made available by Penn Image Computing and Science Lab at
University of Pennsylvania. OASIS-3 data were provided by Longitudinal
Multimodal Neuroimaging: Principal Investigators: T. Benzinger, D. Marcus, J.
Morris; NIH P30 AG066444, P50 AG00561, P30 NS09857781, P01 AG026276, P01
AG003991, R01 AG043434, UL1 TR000448, R01 EB009352.

##  Appendix A An Abstract Overview of VNN-based Brain Age Prediction

Figure 11 provides an abstract overview of the general procedure of evaluating
brain age using machine learning (ML) models. From Fig. 11, we note that if
the ML model is a black box, it may be infeasible to capture the contributors
to elevated age-gap in Step 3. Furthermore, in this context, it is also
unclear whether age-bias correction step influences final $\Delta$-Age
prediction through some statistical artifact [18]. Hence, it can be desirable
to minimize the role of age-bias correction in $\Delta$-Age evaluation by
selecting an ML model that achieves a near perfect fit on chronological age of
healthy controls in Step 1. However, there is no guarantee that achieving an
‘perfect fit’ on true age of healthy controls will enable the ML model to
capture the impact of neurodegeneration in individuals with neurodegeneration.

Figure 11: A general overview of brain age evaluation using machine learning
algorithms in the existing literature. Step 1 consists of training a machine
learning (ML) model to predict chronological age (true age) for healthy
controls. If the correlation between predicted age and true age is smaller
than $1$, an age-bias exists in ML model outputs as the age for older
individuals tends to be under-estimated and that for younger individuals tends
to be over-estimated. To correct for this bias, a linear regression based
model is applied on the ML model outputs in Step 2. Under the hypothesis that
ML model can capture accelerated aging in age-related neurodegeneration, it is
expected that $\Delta$-Age for individuals with neurodegeneration will be
significantly higher than those of healthy controls (Step 3).

VNNs allow us to analyze the contribution of each feature (brain region) to
the final output. Hence, by analyzing the elevations in contributions of
different brain regions via studying group differences in regional residuals,
we are able to characterize the brain regions that contribute to accelerated
aging (or larger $\Delta$-Age) in individuals with age-related
neurodegeneration (Fig. 12). Thus, we can verify that VNNs captured
neurodegeneration-driven effects that eventually led to elevated $\Delta$-Age
for an individual. Our experiments show that VNNs do not obtain a perfect fit
on chronological age of healthy individuals. Hence, age-bias correction is
important to appropriately project the VNN model outputs via a linear model
into an appropriate space such that a clinician can observe an elevated
$\Delta$-Age effect in individuals with neurodegenerative condition (AD in
this paper). Based on these observations, we remark that VNNs provide an
interpretable framework for brain age prediction.

Figure 12: Interpretability offered by VNNs in brain age prediction. By
analyzing the final layer outputs of VNNs, we can isolate brain regions that
have larger regional residuals for individuals with AD with respect to healthy
controls. Furthermore, the elevated regional residuals in these brain regions
eventually contribute to elevated $\Delta$-Age after age-bias correction.

##  Appendix B Proof of Theorem 2

We reprise the proof from [1]. To start with, we note that the coVariance
filters with respect to $\hat{\mathbf{C}}$ and ${\mathbf{C}}$ are given by

$\displaystyle{\mathbf{H}}(\hat{\mathbf{C}})=\sum\limits_{k=0}^{m}h_{k}\hat{\mathbf{C}}^{k}\quad\text{and}\quad{\mathbf{H}}({\mathbf{C}})=\sum\limits_{k=0}^{m}h_{k}{\mathbf{C}}^{k}\;.$
(37)

We start by characterizing the perturbation of sample covariance matrix
$\hat{\mathbf{C}}$ with respect to ${\mathbf{C}}$ in Lemma 1. To this end, we
define

$\displaystyle{\mathbf{E}}\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\hat{\mathbf{C}}-{\mathbf{C}}\;,$
(38)

and ${\mathbf{I}}_{m}$ as an $m\times m$ identity matrix. Also, the eigenvalue
decomposition of sample covariance matrix $\hat{\mathbf{C}}$ is given by

$\displaystyle\hat{\mathbf{C}}=\hat{\mathbf{V}}\hat{\Lambda}\hat{\mathbf{V}}^{\sf
T}\;,$ (39)

where $\hat{\mathbf{V}}=[\hat{\mathbf{v}}_{1},\dots,\hat{\mathbf{v}}_{m}]$ is
the matrix constituted by orthonormal eigenvectors of $\hat{\mathbf{C}}$ and
$\hat{\Lambda}={\sf diag}(\hat{\lambda}_{1},\dots,\hat{\lambda}_{m})$ is the
diagonal matrix of eigenvalues of $\hat{\mathbf{C}}$, such that,
$\hat{\lambda}_{1}\geq\hat{\lambda}_{2}\dots\geq\hat{\lambda}_{m}$. Clearly,
the eigenvalues and eigenvectors of $\hat{\mathbf{C}}$ are estimates of the
eigenvalues and eigenvectors of the true covariance matrix ${\mathbf{C}}$.

######  Lemma 2.

Consider an ensemble covariance matrix ${\mathbf{C}}$ with the
eigendecomposition in (3) and a sample covariance matrix $\hat{\mathbf{C}}$
with the eigendecomposition in (39). For any eigenvalue $\lambda_{i}>0$ of
${\mathbf{C}}$, the perturbation ${\mathbf{E}}$ satisfies

$\displaystyle{\mathbf{E}}{\mathbf{v}}_{i}=\beta_{i}\delta{\mathbf{v}}_{i}+\delta\lambda_{i}{\mathbf{v}}_{i}+(\delta\lambda_{i}{\mathbf{I}}_{m}-{\mathbf{E}})\delta{\mathbf{v}}_{i}$
(40)

where

$\displaystyle\beta_{i}\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}(\lambda_{i}{\mathbf{I}}_{m}-{\mathbf{C}}),\quad\delta{\mathbf{v}}_{i}\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\hat{\mathbf{v}}_{i}-{\mathbf{v}}_{i},\quad\delta\lambda_{i}\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\hat{\lambda}_{i}-\lambda_{i}\;.$
(41)

###### Proof.

Note that from the definition of eigenvectors and eigenvalues, we have

$\displaystyle\hat{\mathbf{C}}\hat{\mathbf{v}}_{i}=\hat{\lambda}_{i}\hat{\mathbf{v}}_{i}\;.$
(42)

We can rewrite (42) in terms of perturbations with respect to the ensemble
covariance matrix ${\mathbf{C}}$ and the outputs of its eigendecomposition as
follows:

$\displaystyle(\hat{\mathbf{C}}-{\mathbf{C}})({\mathbf{v}}_{i}+\delta{\mathbf{v}}_{i})+{\mathbf{C}}({\mathbf{v}}_{i}+\delta{\mathbf{v}}_{i})=(\lambda_{i}+\delta\lambda_{i})({\mathbf{v}}_{i}+\delta{\mathbf{v}}_{i})\;,$
(43)

where we have used $\hat{\lambda}_{i}=\lambda_{i}+\delta\lambda_{i}$ and
$\hat{\mathbf{v}}_{i}={\mathbf{v}}_{i}+\delta{\mathbf{v}}_{i}$. Using the fact
that ${\mathbf{C}}{\mathbf{v}}_{i}=\lambda_{i}{\mathbf{v}}_{i}$ and
rearranging the terms in (43), we have

$\displaystyle(\hat{\mathbf{C}}-{\mathbf{C}}){\mathbf{v}}_{i}=(\lambda_{i}{\mathbf{I}}_{m}-{\mathbf{C}})\delta{\mathbf{v}}_{i}+\delta\lambda_{i}({\mathbf{v}}_{i}+\delta{\mathbf{v}}_{i})-(\hat{\mathbf{C}}-{\mathbf{C}})\delta{\mathbf{v}}_{i}\;.$
(44)

By setting ${\mathbf{E}}=\hat{\mathbf{C}}-{\mathbf{C}}$ and
$\beta_{i}=\lambda_{i}{\mathbf{I}}_{m}-{\mathbf{C}}$, we can rewrite (44) as

$\displaystyle{\mathbf{E}}{\mathbf{v}}_{i}=\beta_{i}\delta{\mathbf{v}}_{i}+\delta\lambda_{i}{\mathbf{v}}_{i}+(\delta\lambda_{i}{\mathbf{I}}_{m}-{\mathbf{E}})\delta{\mathbf{v}}_{i}\;.$
(45)

∎

We next establish the first order approximation for $\hat{\mathbf{C}}^{k}$ in
terms of ${\mathbf{C}}$ and ${\mathbf{E}}$. The first order approximation of
$\hat{\mathbf{C}}^{k}$ is given by

$\displaystyle({\mathbf{C}}+{\mathbf{E}})^{k}={\mathbf{C}}^{k}+\sum\limits_{r=0}^{m}{\mathbf{C}}^{r}{\mathbf{E}}{\mathbf{C}}^{k-r-1}+\tilde{\mathbf{E}}\;,$
(46)

where $\tilde{\mathbf{E}}$ satisfies
$\|\tilde{\mathbf{E}}\|\leq\sum\limits_{r=2}^{k}{k\choose
r}\|{\mathbf{E}}\|^{r}\|{\mathbf{C}}\|^{k-r}$. Using (46), we have

$\displaystyle{\mathbf{H}}(\hat{\mathbf{C}})-{\mathbf{H}}({\mathbf{C}})$
$\displaystyle=\sum\limits_{k=0}^{m}h_{k}[({\mathbf{C}}+{\mathbf{E}})^{k}-{\mathbf{C}}^{k}]\;,$
(47)
$\displaystyle=\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}{\mathbf{E}}{\mathbf{C}}^{k-r-1}+\tilde{\mathbf{E}}\;,$
(48)

where $\tilde{\mathbf{E}}$ satisfies $\|\tilde{\mathbf{E}}\|^{2}={\cal
O}(\|{\mathbf{E}}\|^{2})$ [38]. The focus of our subsequent analysis will be
the first term in (48). For a random data sample
${\mathbf{x}}=[x_{1},\dots,x_{m}]^{\sf T}$, such that,
$\|{\mathbf{x}}\|_{2}\leq 1$ and ${\mathbf{x}}\in{\mathbb{R}}^{m\times 1}$,
its Fourier transform with respect to ${\mathbf{C}}$ is given by
$\tilde{\mathbf{x}}={\mathbf{V}}^{\sf T}{\mathbf{x}}$, where
$\tilde{\mathbf{x}}=[\tilde{x}_{1},\dots,\tilde{x}_{m}]^{\sf T}$ [1]. The
relationship $\tilde{\mathbf{x}}$ and ${\mathbf{x}}$ can be expressed as

$\displaystyle{\mathbf{x}}=\sum\limits_{i=1}^{m}\tilde{x}_{i}{\mathbf{v}}_{i}\;.$
(49)

Multiplying both sides in (48) by ${\mathbf{x}}$ and by leveraging (49), we
get

$\displaystyle[{\mathbf{H}}(\hat{\mathbf{C}})-{\mathbf{H}}({\mathbf{C}})]{\mathbf{x}}$
$\displaystyle=\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}{\mathbf{E}}{\mathbf{C}}^{k-r-1}{\mathbf{x}}+\tilde{\mathbf{E}}{\mathbf{x}}\;,$
(50)
$\displaystyle=\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k=1}{\mathbf{C}}^{r}{\mathbf{E}}{\mathbf{C}}^{k-r-1}{\mathbf{v}}_{i}+\tilde{\mathbf{E}}{\mathbf{x}}\;,$
(51)
$\displaystyle=\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}{\mathbf{E}}{\mathbf{v}}_{i}+\tilde{\mathbf{E}}{\mathbf{x}}\;,$
(52)

where we have used ${\mathbf{C}}{\mathbf{v}}_{i}=\lambda_{i}{\mathbf{v}}_{i}$
in the transition from (51) to (52). We focus only on the first term in (52)
and leverage the result from Lemma 1 that expands
${\mathbf{E}}{\mathbf{v}}_{i}$ to get

$\displaystyle\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}{\mathbf{E}}{\mathbf{v}}_{i}$
$\displaystyle=\underbrace{\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}\beta_{i}\delta{\mathbf{v}}_{i}}_{\text{Term
1}}$
$\displaystyle\enskip+\underbrace{\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}\delta\lambda_{i}{\mathbf{v}}_{i}}_{\text{Term
2}}$
$\displaystyle\enskip+\underbrace{\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}(\delta\lambda_{i}{\mathbf{I}}_{m}-{\mathbf{E}})\delta{\mathbf{v}}_{i}}_{\text{Term
3}}\;.$ (53)

Next, we analyze term 1, term 2, and term 3 in (B) separately.

Analysis of Term 1 in (B). In the analysis of term 1, we start by noting that

$\displaystyle\beta_{i}$
$\displaystyle=\lambda_{i}{\mathbf{I}}_{m}-{\mathbf{C}}\;,$ (54)
$\displaystyle=\sum\limits_{j=1}^{m}(\lambda_{i}-\lambda_{j}){\mathbf{v}}_{j}{\mathbf{v}}_{j}^{\sf
T}\;,$ (55)
$\displaystyle={\mathbf{V}}(\lambda_{i}{\mathbf{I}}_{m}-\Lambda){\mathbf{V}}^{\sf
T}\;.$ (56)

Using (56) and $\delta{\mathbf{v}}_{i}={\mathbf{u}}_{i}-{\mathbf{v}}_{i}$ in
term 1 in (B), we have

$\displaystyle\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}{\mathbf{V}}(\lambda_{i}{\mathbf{I}}_{m}-\Lambda){\mathbf{V}}^{\sf
T}(\hat{\mathbf{v}}_{i}-{\mathbf{v}}_{i})\;.$ (57)

Using ${\mathbf{C}}^{r}={\mathbf{V}}\Lambda^{r}{\mathbf{V}}^{\sf T}$ in (57),
term 1 in (B) is equivalent to

$\displaystyle\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}\lambda_{i}^{k-r-1}{\mathbf{V}}\Lambda^{r}(\lambda_{i}{\mathbf{I}}_{m}-\Lambda){\mathbf{V}}^{\sf
T}(\hat{\mathbf{v}}_{i}-{\mathbf{v}}_{i})\;,$ (58)
$\displaystyle=\sum\limits_{i=1}^{m}\tilde{x}_{i}{\mathbf{V}}{\mathbf{L}}_{i}{\mathbf{V}}^{\sf
T}(\hat{\mathbf{v}}_{i}-{\mathbf{v}}_{i})\;,$ (59)

where ${\mathbf{L}}_{i}$ is a diagonal matrix whose $j$-th element is given by

$\displaystyle[{\mathbf{L}}_{i}]_{j}$
$\displaystyle=\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}(\lambda_{i}-\lambda_{j})\lambda_{i}^{k-r-1}\lambda_{j}^{r}\;,$
(60)
$\displaystyle=\sum\limits_{k=0}^{m}h_{k}(\lambda_{i}-\lambda_{j})\frac{\lambda_{i}^{k}-\lambda_{j}^{k}}{\lambda_{i}-\lambda_{j}}\;,$
(61)
$\displaystyle=\sum\limits_{k=0}^{m}h_{k}\lambda_{i}^{k}-\sum\limits_{k=0}^{m}h_{k}\lambda_{j}^{k}\;,$
(62) $\displaystyle=h(\lambda_{i})-h(\lambda_{j})\;,$ (63)

where $h(\lambda_{i})$ is the frequency response of the coVariance filter and
is defined in (10). Therefore, we have ${\mathbf{L}}_{i}={\sf
diag}([h(\lambda_{i})-h(\lambda_{j})]_{j})$. Next, in (59), we note that

$\displaystyle{\mathbf{V}}^{\sf
T}(\hat{\mathbf{v}}_{i}-{\mathbf{v}}_{i})=[{\mathbf{v}}_{1}^{\sf
T}(\hat{\mathbf{v}}_{i}-{\mathbf{v}}_{i}),\cdots,{\mathbf{v}}_{m}^{\sf
T}(\hat{\mathbf{v}}_{i}-{\mathbf{v}}_{i})]^{\sf T}\;.$ (64)

Using (64) and (63) in (59) and ${\mathbf{v}}_{j}^{\sf
T}{\mathbf{v}}_{i}=0,\forall j\neq i$, we deduce that the term 1 in (B) is
equivalent to

$\displaystyle\sum\limits_{i=1}^{m}\tilde{x}_{i}{\mathbf{V}}{\mathbf{L}}_{i}{\mathbf{V}}^{\sf
T}(\hat{\mathbf{v}}_{i}-{\mathbf{v}}_{i})=\sum\limits_{i=1}^{m}\tilde{x}_{i}{\mathbf{V}}{\mathbf{J}}_{i}\;,$
(65)

where the $j$-th element of ${\mathbf{J}}_{i}$ is given by

$\displaystyle[{\mathbf{J}}_{i}]_{j}=\begin{cases}&0\;,\quad\text{if }\enskip
j=i\;,\\\ &(h(\lambda_{i})-h(\lambda_{j})){\mathbf{v}}_{j}^{\sf
T}\hat{\mathbf{v}}_{i}\;,\enskip\text{otherwise}\end{cases}\;.$ (66)

For the stability analysis, we are interested in the norm of term 1.
Therefore, by noting the equivalence between the term 1 in (B) and (65), after
taking the norm, we have

$\displaystyle\left\lVert\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}\beta_{i}\delta{\mathbf{v}}_{i}\right\rVert_{2}$
$\displaystyle=\left\lVert\sum\limits_{i=1}^{m}\tilde{x}_{i}{\mathbf{V}}{\mathbf{J}}_{i}\right\rVert_{2}\;,$
(67)
$\displaystyle\leq\sqrt{m}\sum\limits_{i=1}^{m}|\tilde{x}_{i}|\max_{j,i\neq
j}|h(\lambda_{i})-h(\lambda_{j})||{\mathbf{v}}_{j}^{\sf
T}\hat{\mathbf{v}}_{i}|\;.$ (68)

Note that ${\mathbf{v}}_{j}^{\sf T}\hat{\mathbf{v}}_{i}$ is the inner product
between the eigenvector ${\mathbf{v}}_{j}$ of the ensemble covariance matrix
${\mathbf{C}}$ and the eigenvector $\hat{\mathbf{v}}_{i}$ of the sample
covariance matrix $\hat{\mathbf{C}}$. The bounds on ${\mathbf{v}}_{j}^{\sf
T}\hat{\mathbf{v}}_{i}$ in terms of the number of data samples $n$ have been
studied in the existing literature. Here, we leverage the result from [78,
Theorem 4.1] to conclude that if ${\sf
sgn}(\lambda_{j}-\lambda_{i})2\hat{\lambda}_{j}>{\sf
sgn}(\lambda_{j}-\lambda_{i})(\lambda_{j}-\lambda_{i})$ for
$\lambda_{i}\neq\lambda_{j}$, the condition

$\displaystyle\left\lVert\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}\beta_{i}\delta{\mathbf{v}}_{i}\right\rVert_{2}\leq\sqrt{m}\sum\limits_{i=1}^{m}|\tilde{x}_{i}|\max_{j,i\neq
j}|h(\lambda_{i})-h(\lambda_{j})|\frac{2k_{i}}{n^{1/2-\varepsilon}|\lambda_{i}-\lambda_{j}|}\;,$
(69)

is true with probability at least $\left(1-\frac{1}{n^{2\varepsilon}}\right)$
for some $\varepsilon\in(0,1/2]$, where
$k_{i}=\Big{(}{\mathbb{E}}[\|{\mathbf{X}}{\mathbf{X}}^{\sf
T}{\mathbf{v}}_{i}\|_{2}^{2}]-\lambda_{i}^{2}\Big{)}^{\frac{1}{2}}$.
Furthermore, we note that the condition ${\sf
sgn}(\lambda_{j}-\lambda_{i})2\hat{\lambda}_{j}>{\sf
sgn}(\lambda_{j}-\lambda_{i})(\lambda_{j}-\lambda_{i})$ is satisfied with
probability at least $1-\frac{2k_{i}^{2}}{|\lambda_{i}-\lambda_{j}|}$ [78,
Corollary 4.2], which via a union bound and first order approximation from
Taylor series implies that (69) is true with probability at least
$1-\frac{1}{n^{2\varepsilon}}-\frac{2\kappa m}{n}$ for $\kappa$ defined as

$\displaystyle\kappa\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\max_{i,j:\lambda_{i}\neq\lambda_{j}}\frac{k_{i}^{2}}{|\lambda_{i}-\lambda_{j}|}\quad\text{where}\quad
k_{\sf
min}\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\min_{i\in\\{1,\dots,m\\},\lambda_{i}>0}k_{i}\;.$
(70)

Therefore, for a coVariance filter with the property

$\displaystyle\frac{|h(\lambda_{i})-h(\lambda_{j})|}{|\lambda_{i}-\lambda_{j}|}\leq\frac{M}{k_{i}}\;,$
(71)

for some real constant $M>0$, the condition in (69) is equivalent to

$\displaystyle\left\lVert\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}\beta_{i}\delta{\mathbf{v}}_{i}\right\rVert_{2}\leq\frac{2\sqrt{m}M}{n^{\frac{1}{2}-\varepsilon}}\sum\limits_{i=1}^{m}|\tilde{x}_{i}|\;,$
(72)

which holds with probability at least
$1-\frac{1}{n^{2\varepsilon}}-\frac{2\kappa m}{n}$. Furthermore, note that
$\sum\limits_{i=1}^{m}|\tilde{x}_{i}|\leq\sqrt{m}\|{\mathbf{x}}\|_{2}$. When
the random sample ${\mathbf{x}}$ satisfies $\|{\mathbf{x}}\|_{2}\leq 1$, we
have

$\displaystyle{\mathbb{P}}\left(\left\lVert\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}\beta_{i}\delta{\mathbf{v}}_{i}\right\rVert_{2}\leq\frac{2}{n^{\frac{1}{2}-\varepsilon}}mM\right)\geq
1-\frac{1}{n^{2\varepsilon}}-\frac{2\kappa m}{n}\;,$ (73)

for any $\varepsilon\in(0,1/2]$.

Analysis of Term 2 in (B). Using
${\mathbf{C}}{\mathbf{v}}_{i}=\lambda_{i}{\mathbf{v}}_{i}$, we note that term
2 in (B) is equivalent to

$\displaystyle\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}\delta\lambda_{i}{\mathbf{v}}_{i}$
$\displaystyle=\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}\lambda_{i}^{k-1}\delta\lambda_{i}{\mathbf{v}}_{i}\;,$
(74)
$\displaystyle=\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}kh_{k}\lambda_{i}^{k-1}\delta\lambda_{i}{\mathbf{v}}_{i}\;,$
(75)
$\displaystyle=\sum\limits_{i=1}^{m}\tilde{x}_{i}h^{\prime}(\lambda_{i})\delta\lambda_{i}{\mathbf{v}}_{i}\;.$
(76)

Next, using Weyl’s theorem [115, Theorem 8.1.6], we note that
$\|{\mathbf{E}}\|\leq\alpha$ implies that $|\delta\lambda_{i}|\leq\alpha$ for
any $\alpha>0$. For a random instance ${\mathbf{x}}$ of a random vector
${\mathbf{X}}$ whose probability distribution is supported within a ball of
radius $1$ w.l.o.g, such that, $\|{\mathbf{x}}\|_{2}\leq 1$, we have

$\displaystyle{\mathbb{P}}\left({\mathbf{E}}\leq
B\left(\frac{\|{\mathbf{C}}\|\sqrt{\log
m+u}}{\sqrt{n}}+\frac{(1+\|{\mathbf{C}}\|)(\log m+u)}{n}\right)\right)\geq
1-2^{-u}\;,$ (77)

for some constant $B>0$ and $u>0$. The result in (77) follows directly from
[116, Theorem 5.6.1]. Therefore, using (76), we have

$\displaystyle\left\lVert\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}\delta\lambda_{i}{\mathbf{v}}_{i}\right\rVert_{2}$
$\displaystyle\leq\sum\limits_{i=1}^{m}|\tilde{x}_{i}||h^{\prime}(\lambda_{i})||\delta\lambda_{i}|\|{\mathbf{v}}_{i}\|_{2}\;.$
(78)

Using (77), $|h^{\prime}(\lambda_{i})|\leq M/k_{\sf min}$ (where $k_{\sf
min}=\min_{i\in\\{1,\dots,m\\},\lambda_{i}>0}k_{i}$) from (71), and
$\|{\mathbf{v}}_{i}\|_{2}=1$, we have

$\displaystyle{\mathbb{P}}\left(\left\lVert\sum\limits_{i=1}^{m}\tilde{x}_{i}\sum\limits_{k=0}^{m}h_{k}\sum\limits_{r=0}^{k-1}{\mathbf{C}}^{r}\lambda_{i}^{k-r-1}\delta\lambda_{i}{\mathbf{v}}_{i}\right\rVert_{2}\right.$
$\displaystyle\quad\left.\leq\frac{A}{k_{\sf
min}}\sqrt{m}M\left(\frac{\|{\mathbf{C}}\|\sqrt{\log
m+u}}{\sqrt{n}}+\frac{(1+\|{\mathbf{C}}\|)(\log m+u)}{n}\right)\right)\geq
1-2^{-u}\;,$ (79)

for some constant $A>0$ and $u>0$.

Analysis of Term 3 in (B). We remark that the term 3 in (B) consists of second
order error terms that diminish faster with the number of samples $n$ as
compared to term 1 and term 2. This can be concluded from the following
observations. Firstly, from (77) and Weyl’s theorem, we note that
$\|\delta\lambda_{i}{\mathbf{I}}_{m}-{\mathbf{E}}\|\leq 2\|{\mathbf{E}}\|$ and
$\|{\mathbf{E}}\|\simeq{\cal O}(1/\sqrt{n})$ with high probability. Secondly,
the upper bound on the term $\delta{\mathbf{v}}_{i}$ also has a similar
scaling behavior as ${\mathbf{E}}$ [117].

Therefore, the second order error term
$(\delta\lambda_{i}{\mathbf{I}}_{m}-{\mathbf{E}})\delta{\mathbf{v}}_{i}$
diminishes at a rate faster than ${\cal O}(1/\sqrt{n})$, which is faster as
compared to terms 1 and 2, that individually scale as ${\cal
O}(1/n^{1/2-\varepsilon})$ for $\varepsilon\in(0,1/2]$ and ${\cal
O}(1/\sqrt{n})$, respectively. Finally, by noting that the condition on
$\|[{\mathbf{H}}(\hat{\mathbf{C}})-{\mathbf{H}}({\mathbf{C}})]{\mathbf{x}}\|_{2}$
reduces to the condition on operator norm
$\|[{\mathbf{H}}(\hat{\mathbf{C}})-{\mathbf{H}}({\mathbf{C}})]\|$ over the
search space $\max_{{\mathbf{x}}\in{\mathbb{R}}^{m\times
1},\|{\mathbf{x}}\|_{2}\leq
1}\\{\|[{\mathbf{H}}(\hat{\mathbf{C}})-{\mathbf{H}}({\mathbf{C}})]{\mathbf{x}}\|_{2}\\}$
and that the terms scaling at $1/\sqrt{n}$ or slower in (73) and (B) dominate
the scaling behavior of the upper bound on
$\|{\mathbf{H}}(\hat{\mathbf{C}})-{\mathbf{H}}({\mathbf{C}})\|$, we arrive at
the following theorem.

######  Theorem 4 (Stability of coVariance Filter).

Consider a random vector ${\mathbf{X}}\in{\mathbb{R}}^{m\times 1}$ , such
that, its corresponding covariance matrix is given by
${\mathbf{C}}={\mathbb{E}}[({\mathbf{X}}-{\mathbb{E}}[{\mathbf{X}}])({\mathbf{X}}-{\mathbb{E}}[{\mathbf{X}}])^{\sf
T}]$. For a sample covariance matrix $\hat{\mathbf{C}}$ formed using $n$ i.i.d
instances of ${\mathbf{X}}$ and a random instance ${\mathbf{x}}$ of
${\mathbf{X}}$, such that, $\|{\mathbf{x}}\|_{2}\leq 1$ and under assumption
(71), the following holds with probability at least
$(1-{n^{-2\varepsilon}}-2\kappa m/n)(1-1/n)$ for any $\varepsilon\in(0,1/2]$:

$\displaystyle\left\lVert{\mathbf{H}}(\hat{\mathbf{C}})-{\mathbf{H}}({\mathbf{C}})\right\rVert=\frac{M}{n^{\frac{1}{2}-\varepsilon}}\cdot{\cal
O}\left(m+\frac{\sqrt{m}\|{\mathbf{C}}\|\sqrt{\log mn}}{k_{\sf
min}n^{\varepsilon}}\right)\;.$ (80)

Thus, the right-hand side of (80) provides the structure of $\alpha_{n}$ in
Theorem 2. From Theorem 4, we note that
$\|{\mathbf{H}}(\hat{\mathbf{C}})-{\mathbf{H}}({\mathbf{C}})\|$ decays with
the number of samples $n$ at least at the rate of
$1/n^{\frac{1}{2}-\varepsilon}$. Thus, we conclude that the stability of the
coVariance filter improves as the number of samples $n$ increases. This
observation is along the expected lines as the estimate $\hat{\mathbf{C}}$
becomes closer to the ensemble covariance matrix ${\mathbf{C}}$ by the virtue
of the law of large numbers. We refer the reader to [1] for further discussion
on the implication of the assumption in (71) on filter design. The rest of the
proof follows directly from [1, Theorem 3].

##  Appendix C Graphon Information Processing

The theory of graphons has previously been leveraged to study the
transferability of GNNs between graphs in the same graphon family [39]. The
proof of Theorem 3 relies on establishing the transferability of VNNs between
datasets in the setting where their corresponding covariance matrices belong
to a converging sequence characterized by a graphon. Our main objective in
this section is to show that data processing over coVariance filter can
equivalently be represented in the continuous domain using its graphon
approximation. Establishing this property will ultimately allow us to compare
VNNs instantiated on covariance matrices derived from datasets with different
numbers of features. We begin with some preliminaries for graphons.

###  C.1 Preliminaries.

The definition of a graphon is re-stated below.

######  Definition 5 (Graphon).

A graphon is a bounded, symmetric, measurable function
${{\mathbf{W}}:[0,1]^{2}\mapsto[-1,1]}$.

Using the theory of convergence of graphons and interpreting the covariance
matrix as a weighted graph representation of data, a graphon ${\mathbf{W}}$
exists as a limiting object for the sequence of graphon approximations
$\\{{\mathbf{W}}_{{\mathbf{C}}_{m}}\\}$ if the sequence of covariance matrices
$\\{{\mathbf{C}}_{m}\\}$ converges in the _cut distance_ [80]. In Remark 2, we
formalize the sufficient condition for existence of a limit object for a given
sequence of covariance matrices. The statement in Remark 2 is an extension of
[80, Corollary 3.9] to our setting where covariance matrices are viewed as
weighted graphs.

######  Remark 2 (Graphon as limit object [80]).

A sequence of covariance matrices $\\{{\mathbf{C}}_{m}\\}$ is deemed
convergent if they form a Cauchy sequence with respect to the cut distance
$\delta_{\square}$, where the cut distance
$\delta_{\square}({\mathbf{C}}_{m_{1}},{\mathbf{C}}_{m_{2}})$ between
covariance matrices ${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$ is
defined in (93) in Appendix D. Furthermore, for any convergent sequence of
covariance matrices $\\{{\mathbf{C}}_{m}\\}$, the corresponding sequence of
graphon approximations $\\{{\mathbf{W}}_{{\mathbf{C}}_{m}}\\}$ converges to a
graphon.

A distinct feature of the cut distance is that it allows the comparison of
covariance matrices of different sizes. Hence, all covariance matrices whose
graphon approximations converge to a graphon can be considered to be a part of
that graphon family. Moreover, graphon ${\mathbf{W}}$ can be interpreted as
the schema for which the covariance matrix ${\mathbf{C}}_{m}$ represents the
covariance realization at resolution $m$.

###  C.2 Information Processing with Graphons

We next show that a coVariance filter ${\mathbf{H}}({\mathbf{C}}_{m})$ can be
equivalently represented in the continuous domain using convolution operations
over graphon representations ${\mathbf{W}}_{{\mathbf{C}}_{m}}$. Given a
coVariance filter output
${\mathbf{z}}={\mathbf{H}}({\mathbf{C}}_{m}){\mathbf{x}}$, the continuous
representation of ${\mathbf{x}}$ is $y_{{\mathbf{x}}}$ and that of
${\mathbf{C}}_{m}$ is ${\mathbf{W}}_{{\mathbf{C}}_{m}}$. The operation
${\mathbf{C}}{\mathbf{x}}$ is fundamental to the convolution operation in
${\mathbf{H}}({\mathbf{C}}_{m}){\mathbf{x}}$ and therefore, we first provide
its continuous equivalent. For ${\mathbf{s}}={\mathbf{C}}{\mathbf{x}}$, the
$i$-th element of ${\mathbf{s}}$ is given by

$\displaystyle[{\mathbf{s}}]_{i}=\sum\limits_{j=0}^{m}[{\mathbf{C}}_{m}]_{ij}[{\mathbf{x}}]_{j}\;.$
(81)

Thus, $[{\mathbf{s}}]_{i}$ is a linear combination of elements in
${\mathbf{x}}$ according to the $i$-th row of ${\mathbf{C}}_{m}$. In the
continuous space, we can equivalently write (81) as

$\displaystyle
y_{{\mathbf{s}}}(u)=\int_{0}^{1}{\mathbf{W}}_{{\mathbf{C}}_{m}}(u,v)y_{{\mathbf{x}}}(v){\sf
d}v\;,$ (82)

where $y_{{\mathbf{x}}}$ is the continuous representation of ${\mathbf{x}}$
obtained according to the intervals defined in (20). Note that
$y_{{\mathbf{s}}}$ is a continuous representation of ${\mathbf{s}}$, i.e.,
they satisfy $y_{{\mathbf{s}}}(u)=[{\mathbf{s}}]$ for $u\in{\mathcal{U}}_{i}$.
Hence, $y_{{\mathbf{s}}}$ and ${\mathbf{s}}$ can be recovered from each other.
This observation can be extrapolated to define the continuous equivalent of a
coVariance filter. This is feasible because we can write the entity
${\mathbf{C}}_{m}^{k}{\mathbf{x}}$ in ${\mathbf{H}}({\mathbf{C}})$ in a
recursive form. Specifically, if we have
${\mathbf{s}}_{k}={\mathbf{C}}_{m}^{k}{\mathbf{x}}$, then we can rewrite
${\mathbf{s}}_{k}$ as

$\displaystyle{\mathbf{s}}_{k}={\mathbf{C}}_{m}{\mathbf{s}}_{k-1}\;,$ (83)

where ${\mathbf{s}}_{0}={\mathbf{x}}$. Thus, using the same reasoning that
established the equivalence between (81) and (82), we conclude that the
continuous representation $y_{{\mathbf{s}}_{k}}$ of ${\mathbf{s}}_{k}$ can be
recovered via the following operation

$\displaystyle
y_{{\mathbf{s}}_{k}}(u)=\int_{0}^{1}{\mathbf{W}}_{{\mathbf{C}}_{m}}(u,v)y_{{\mathbf{s}}_{k-1}}(v){\sf
d}v\;.$ (84)

Since the coVariance filter output ${\mathbf{z}}$ is a weighted aggregation of
the terms ${\mathbf{s}}_{k}$, we can write its continuous representation
$y_{{\mathbf{z}}}$ as

$\displaystyle
y_{{\mathbf{z}}}(u)=\sum\limits_{k=0}^{K}h_{k}y_{{\mathbf{s}}_{k}}(u)\;.$ (85)

Using the mathematical steps leading up to (85), we have shown that the
continuous representation of the covariance filter output ${\mathbf{z}}$ can
be recovered via the convolution operations over the graphon representation
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ in (82) and (84). Also, ${\mathbf{z}}$ and
$y_{{\mathbf{z}}}$ are operationally interchangeable. Moreover, we can also
extrapolate this correspondence between ${\mathbf{z}}$ and $y_{{\mathbf{z}}}$
to covariance perceptrons and VNNs with multi-layer architecture and MIMO
information processing. The extension of this observation to coVariance
perceptron and a basic VNN is trivial as the coVariance output is evaluated
after application of pointwise non-linearity $\sigma$ on ${\mathbf{z}}$ and a
basic VNN is formed by stacking multiple coVariance perceptrons and number of
inputs and outputs at each layer (i.e., $F$) being set to $1$.

We use the notation ${\mathbf{x}}_{m}$ to denote an input vector with $m$
features. Thus, if VNN output
$\Phi({\mathbf{x}}_{m};{\mathbf{C}}_{m},{\mathcal{H}})$ is of size $m\times 1$
and we have $F=1$ and number of layers $L$, its continuous approximation
$y_{\Phi({\mathbf{x}}_{m};{\mathbf{C}}_{m},{\mathcal{H}})}$ can be recovered
by a convolutional architecture instantiated on
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ with input $y_{{\mathbf{x}}_{m}}$. For a VNN
with MIMO processing, each VNN layer has multiple $m$-dimensional inputs and
multiple $m$-dimensional outputs. Thus, we can equivalently define an
architecture capable of performing MIMO processing that is instantiated on
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ and ${\mathbf{x}}_{m}$ and produces multiple
continuous representations as the output. Such an architecture has previously
been studied in the form of graphon neural networks [79]. In this context, we
define the model
$\tilde{\Phi}(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$
that is modeled via convolution operations over
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ in (84) and has the same architecture as the
VNN $\Phi({\mathbf{x}}_{m};{\mathbf{C}}_{m},{\mathcal{H}})$. Note that the
outputs of
$\tilde{\Phi}(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$
are continuous representations of the outputs of VNN
$\Phi({\mathbf{x}}_{m};{\mathbf{C}}_{m},{\mathcal{H}})$ (see also Fig. 3 for
an illustration). Thus, we can investigate the transferability of parameters
${\mathcal{H}}$ between VNNs instantiated on covariance matrices
${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$ by analyzing the difference
between
$\tilde{\Phi}(y_{{\mathbf{x}}_{m_{1}}};{\mathbf{W}}_{{\mathbf{C}}_{m_{1}}},{\mathcal{H}})$
and
$\tilde{\Phi}(y_{{\mathbf{x}}_{m_{2}}};{\mathbf{W}}_{{\mathbf{C}}_{m_{2}}},{\mathcal{H}})$.

In this context, our analysis hinges on the setting in which the graphon
approximations ${\mathbf{W}}_{{\mathbf{C}}_{m_{1}}}$ and
${\mathbf{W}}_{{\mathbf{C}}_{m_{2}}}$ belong to a sequence of graphon
approximations $\\{{\mathbf{W}}_{{\mathbf{C}}_{m}}\\}$ that converges to a
graphon ${\mathbf{W}}$. Thus, we also consider an information processing
architecture $\tilde{\Phi}(y;{\mathbf{W}},{\mathcal{H}})$ instantiated on
graphon ${\mathbf{W}}$, such that $y$ and continuous representations
$y_{{\mathbf{x}}_{m}}$ always satisfy
$y_{{\mathbf{x}}_{m}}(\rho_{i})=y(\rho_{i}),\forall i\in\\{1,\dots,m\\}$.
Here, we can also interpret $\tilde{\Phi}(y;{\mathbf{W}},{\mathcal{H}})$ as a
generative model with
$\tilde{\Phi}(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$
being an instance of $\tilde{\Phi}(y;{\mathbf{W}},{\mathcal{H}})$ at
resolution $m$. Thus, our analysis of transferability of VNNs also includes
the study of convergence of outputs from
$\tilde{\Phi}(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$
with that from $\tilde{\Phi}(y;{\mathbf{W}},{\mathcal{H}})$.

To this end, we now formally define a convolution filter over a graphon and
characterize its frequency response. We denote the $k$-hop aggregation
(analogous to ${\mathbf{C}}^{k}{\mathbf{x}}$) on
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ and continuous representation
$y_{{\mathbf{x}}_{m}}$ by the operator
$T_{{\mathbf{W}}_{{\mathbf{C}}_{m}}}^{k}y_{{\mathbf{x}}_{m}}$ that is given by

$\displaystyle(T_{{\mathbf{W}}_{{\mathbf{C}}_{m}}}^{k}y_{{\mathbf{x}}_{m}})(u)\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\int_{0}^{1}{\mathbf{W}}_{{\mathbf{C}}_{m}}(u,v)(T_{{\mathbf{W}}_{{\mathbf{C}}_{m}}}^{k-1}y_{{\mathbf{x}}_{m}})(v){\sf
d}v\;,$ (86)

for any $k>1$, where

$\displaystyle(T_{{\mathbf{W}}_{{\mathbf{C}}_{m}}}y_{{\mathbf{x}}_{m}})(u)\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\int_{0}^{1}{\mathbf{W}}_{{\mathbf{C}}_{m}}(u,v)y_{{\mathbf{x}}_{m}}(v){\sf
d}v\;.$ (87)

Thus, based on the discussion above,
$T^{k}_{{\mathbf{W}}_{{\mathbf{C}}_{m}}}y_{{\mathbf{x}}_{m}}$ and
${\mathbf{C}}_{m}^{k}{\mathbf{x}}_{m}$ are operationally interchangeable. We
can also define $k$-hop aggregation over ${\mathbf{W}}$ using the operator
$T_{{\mathbf{W}}}y$ when $y$ is related to $y_{{\mathbf{x}}_{m}}$ by
$y_{{\mathbf{x}}_{m}}(\rho_{i})=y(\rho_{i})$, where $\rho_{i}$ is defined in
(20). Thus, graphon ${\mathbf{W}}$ and the continuous representation $y$ can
be seen as generative models for covariance matrix ${\mathbf{C}}_{m}$ and data
point ${\mathbf{x}}_{m}$. This observation is in parallel to that in the
context of graphs and graphons [79]. We denote the graphon filter for a set of
filter taps ${\mathcal{H}}=\\{h_{k}\\}_{k=0}^{K}$ by
$\Psi(y;{\mathbf{W}},{\mathcal{H}}):[0,1]\rightarrow{\mathbb{R}}$, which is
defined as

$\displaystyle\Psi(y;{\mathbf{W}},{\mathcal{H}})(u)$
$\displaystyle\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\sum\limits_{k=0}^{K}h_{k}(T_{{\mathbf{W}}}^{k}y)(u)\;.$
(88)

Similar to coVariance filter, we can characterize the frequency response of a
graphon filter via using eigendecomposition of ${\mathbf{W}}$ in (88). Because
${\mathbf{W}}$ is bounded and symmetric, we can express the spectral
decomposition of ${\mathbf{W}}$ as

$\displaystyle{\mathbf{W}}(u,v)=\sum\limits_{i\in{\mathbb{Z}}\backslash\\{0\\}}\eta_{i}\Gamma_{i}(u)\Gamma_{i}(v)\;,$
(89)

where $\eta_{i},\forall i\in{\mathbb{Z}}\backslash\\{0\\}$ are eigenvalues and
$\Gamma_{i}$ are the eigensignals of ${\mathbf{W}}$. Therefore, we can re-
write (88) as

$\displaystyle\Psi(y;{\mathbf{W}},{\mathcal{H}})(u)$
$\displaystyle=\sum\limits_{i\in{\mathbb{Z}}\backslash\\{0\\}}\sum\limits_{k=0}^{K}h_{k}\eta_{i}^{k}\Gamma_{i}(u)\int_{0}^{1}\Gamma_{i}(v)y(v){\sf
d}v\;,$ (90)
$\displaystyle=\sum\limits_{i\in{\mathbb{Z}}\backslash\\{0\\}}\tilde{h}(\eta_{i})\Gamma_{i}(u)\int_{0}^{1}\Gamma_{i}(v)y(v){\sf
d}v\;,$ (91)

for $u\in[0,1]$. Note that (90) follows from (88) using (89) and (86), and we
have used the definition
$\tilde{h}(\eta)\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\sum_{k=0}^{K}h_{k}\eta^{k}$
in (91). The term $\tilde{h}(\eta_{i})$ characterizes the frequency response
of a graphon filter and depends on the filter taps $\\{h_{k}\\}$ and the
graphon eigenvalues. The analysis of
$\|\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})-\Psi(y;{\mathbf{W}},{\mathcal{H}})\|_{2}$
in Appendix E reveals that the closeness between the graphon filter
approximation
$\Psi(y_{{\mathbf{x}}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$ and the
graphon filter $\Psi(y;{\mathbf{W}},{\mathcal{H}})$ scales inversely with the
dimension $m$.

##  Appendix D Cut Distance

Here, we borrow the definition of cut distance between two covariance matrices
${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$ with $m_{1}\neq m_{2}$ from
that for weighted graphs in [80]. First, ${\mathbf{C}}_{m_{1}}$ and
${\mathbf{C}}_{m_{2}}$ are normalized such that ${\sf
tr}({\mathbf{C}}_{m_{1}})=1$ and ${\sf tr}({\mathbf{C}}_{m_{2}})=1$. Next, we
define a family of matrices $m_{1}\times m_{2}$ sized matrices ${\cal
K}({\mathbf{C}}_{m_{1}},{\mathbf{C}}_{m_{2}})$, such that, for every element
${\bf K}\in{\cal K}({\mathbf{C}}_{m_{1}},{\mathbf{C}}_{m_{2}})$, we have

$\displaystyle\sum\limits_{j=1}^{m_{2}}[{\bf
K}]_{ij}=[{\mathbf{C}}_{m_{1}}]_{ii}\quad\text{and}\quad\sum\limits_{i=1}^{m_{1}}[{\bf
K}]_{ij}=[{\mathbf{C}}_{m_{2}}]_{jj}\;.$ (92)

The family ${\cal K}({\mathbf{C}}_{m_{1}},{\mathbf{C}}_{m_{2}})$ is referred
to as ‘fractional overlay’ that conceptually defines the mapping between
individual features of ${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$.
Since the covariance matrix provides a weighted graph representation with
features as nodes, the fractional overlay describes the overlap between the
nodes of two weighted graphs with different number of nodes. Next, by
leveraging a member ${\bf K}$ of fractional overlay ${\cal
K}({\mathbf{C}}_{m_{1}},{\mathbf{C}}_{m_{2}})$, we define matrices
${\mathbf{C}}_{m_{1}}[{\bf K}]$ and ${\mathbf{C}}_{m_{2}}[{\bf K}^{\sf T}]$ on
the set of pairs in $[m_{1}]\times[m_{2}]$, where $[m_{1}]=[1,\dots,m_{1}]$.
Hence, the size of matrices ${\mathbf{C}}_{m_{1}}[{\bf K}]$ and
${\mathbf{C}}_{m_{2}}[{\bf K}]$ is $m_{1}m_{2}\times m_{1}m_{2}$. In both
${\mathbf{C}}_{m_{1}}[{\bf K}]$ and ${\mathbf{C}}_{m_{2}}[{\bf K}^{\sf T}]$,
the weight at the diagonal element associated with the element $(i,j)$ in
$[m_{1}]\times[m_{2}]$ is $[{\bf K}]_{ij}$. Furthermore, for $i,j\in[m_{1}]$
and $a,b\in[m_{2}]$, the off-diagonal element associated with the paired
element $((i,a),(j,b))$ is $[{\mathbf{C}}_{m_{1}}]_{ij}$ in
${\mathbf{C}}_{m_{1}}[{\bf K}]$ and $[{\mathbf{C}}_{m_{2}}]_{ab}$ in
${\mathbf{C}}_{m_{2}}[{\bf K}^{\sf T}]$. Since ${\mathbf{C}}_{m_{1}}[{\bf K}]$
and ${\mathbf{C}}_{m_{2}}[{\bf K}^{\sf T}]$ derived from
${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$ are of the same size, the
distance between them can be readily defined. The cut distance between
${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$ is equivalent to the cut
distance between ${\mathbf{C}}_{m_{1}}[{\bf K}]$ and
${\mathbf{C}}_{m_{2}}[{\bf K}^{\sf T}]$ and is defined as

$\displaystyle\delta_{\square}({\mathbf{C}}_{m_{1}},{\mathbf{C}}_{m_{2}})\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\min_{{\bf
K}\in{\cal
K}({\mathbf{C}}_{m_{1}},{\mathbf{C}}_{m_{2}})}d_{\square}({\mathbf{C}}_{m_{1}}[{\bf
K}],{\mathbf{C}}_{m_{2}}[{\bf K}^{\sf T}]])\;,$ (93)

where $d_{\square}$ is a distance metric between matrices of the same size and
defined next. For two matrices ${\mathbf{C}}_{m}$ and ${\mathbf{D}}_{m}$ with
same diagonal elements, the distance
$d_{\square}({\mathbf{C}}_{m},{\mathbf{D}}_{m})$ is defined as

$\displaystyle
d_{\square}({\mathbf{C}}_{m},{\mathbf{D}}_{m})\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\max_{S,T\in[m]}\frac{1}{{\sf
tr}({\mathbf{C}}_{m})^{2}}|e_{{\mathbf{C}}_{m}}(S,T)-e_{{\mathbf{D}}_{m}}(S,T)|\;,$
(94)

where

$\displaystyle
e_{{\mathbf{C}}_{m}}(S,T)\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\sum\limits_{i\in
S,j\in
T}[{\mathbf{C}}_{m}]_{ii}[{\mathbf{C}}_{m}]_{jj}[{\mathbf{C}}_{m}]_{ij}\;.$
(95)

##  Appendix E Proof of Theorem 3

In Theorem 3, we compare the continuous representations of the $f$-th outputs
of VNNs $\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ and
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$. Our
discussion in Appendix C showed that these continuous representations appear
naturally as the outputs of the architectures
$\tilde{\Phi}(y_{{\mathbf{x}}_{m_{1}}};{\mathbf{W}}_{{\mathbf{C}}_{m_{1}}},{\mathcal{H}})$
and
$\tilde{\Phi}(y_{{\mathbf{x}}_{m_{1}}};{\mathbf{W}}_{{\mathbf{C}}_{m_{1}}},{\mathcal{H}})$
instantiated on graphon approximations ${\mathbf{W}}_{{\mathbf{C}}_{m_{1}}}$
and ${\mathbf{W}}_{{\mathbf{C}}_{m_{2}}}$, respectively. Therefore, our
subsequent analysis is focused on the comparisons between their constituent
graphon filters (see Appendix C for definition) that eventually enables us to
establish the convergence between $f$-th outputs of
$\tilde{\Phi}(y_{{\mathbf{x}}_{m_{1}}};{\mathbf{W}}_{{\mathbf{C}}_{m_{1}}},{\mathcal{H}})$
and
$\tilde{\Phi}(y_{{\mathbf{x}}_{m_{1}}};{\mathbf{W}}_{{\mathbf{C}}_{m_{1}}},{\mathcal{H}})$.
We refer the reader to Appendix C for the details on information processing
architecture defined on the graphon approximations and understanding of the
relationship between
$\tilde{\Phi}(y_{{\mathbf{x}}_{m_{1}}};{\mathbf{W}}_{{\mathbf{C}}_{m_{1}}},{\mathcal{H}})$
and corresponding VNN
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$.

We begin by establishing various results pertaining to the comparisons between
${\mathbf{W}}$ and ${\mathbf{W}}_{{\mathbf{C}}_{m}}$, $y$ and
$y_{{\mathbf{x}}_{m}}$, and difference between eigenvalues of two distinct
graphons. We leverage the $(\Omega,\zeta)$-dominant property of sequence of
covariance matrices $\\{{\mathbf{C}}_{m}\\}$ in (22) and the Lipschitz
condition of graphon in (24) to establish the following result.

######  Lemma 3.

Given an $\alpha$-Lipschitz graphon ${\mathbf{W}}$ and
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ as graphon representation of a
$(\Omega,\zeta)$-dominant covariance matrix ${\mathbf{C}}_{m}$, we have

$\displaystyle\|{\mathbf{W}}-{\mathbf{W}}_{{\mathbf{C}}_{m}}\|_{2}\leq{\frac{\alpha\Omega^{3/2}}{m^{3\zeta/2-1}}}\;.$
(96)

###### Proof.

From the construction of ${\mathbf{W}}_{{\mathbf{C}}_{m}}$, we have

$\displaystyle\|{\mathbf{W}}-{\mathbf{W}}_{{\mathbf{C}}_{m}}\|_{2}$
$\displaystyle=\Big{(}\int_{0}^{1}\int_{0}^{1}\|{\mathbf{W}}(u,v)-{\mathbf{W}}_{{\mathbf{C}}_{m}}(u,v)\|^{2}{\sf
d}u{\sf d}v\Big{)}^{\frac{1}{2}}\;,$ (97)
$\displaystyle=\Big{(}\sum\limits_{i,j}\int_{{\mathcal{U}}_{i}}\int_{{\mathcal{U}}_{j}}\|{\mathbf{W}}(u,v)-{\mathbf{W}}_{{\mathbf{C}}_{m}}(u,v)\|^{2}{\sf
d}u{\sf d}v\Big{)}^{\frac{1}{2}}\;.$ (98)

Without loss of generality, we assume that ${\mathcal{U}}_{1}=[0,\rho_{1}]$ is
the largest interval. Using the $\alpha$-Lipschitz continuity of graphon
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ and noting that
${\mathbf{W}}_{{\mathbf{C}}_{m}}(\rho_{i},\rho_{j})={\mathbf{W}}(\rho_{i},\rho_{j})$,
we have

$\displaystyle\|{\mathbf{W}}-{\mathbf{W}}_{{\mathbf{C}}_{m}}\|_{2}$
$\displaystyle\leq\Big{(}\sum\limits_{i,j}\int_{I_{i}}\int_{I_{j}}\alpha^{2}(|u|+|v|)^{2}{\sf
d}u{\sf d}v\Big{)}^{\frac{1}{2}}\;,$ (99)
$\displaystyle\leq\Big{(}m^{2}\int_{0}^{\rho_{1}}\int_{0}^{\rho_{1}}\alpha^{2}(|u|+|v|)^{2}{\sf
d}u{\sf d}v\Big{)}^{\frac{1}{2}}\;,$ (100)
$\displaystyle\leq\Big{(}m^{2}\int_{0}^{\rho_{1}}\int_{0}^{\rho_{1}}\alpha^{2}(|u|+|v|){\sf
d}u{\sf d}v\Big{)}^{\frac{1}{2}}\;,$ (101) $\displaystyle\leq{\alpha
m\rho_{1}^{3/2}}\;.$ (102)

Using the assumption that ${\mathbf{C}}_{m}$ is $(\Omega,\zeta)$-dominant, we
have

$\displaystyle\|{\mathbf{W}}-{\mathbf{W}}_{{\mathbf{C}}_{m}}\|_{2}$
$\displaystyle\leq\frac{\alpha\Omega^{3/2}}{m^{3\zeta/2-1}}\;.$ (103)

∎

Next, we characterize the difference between a graphon signal $y\in
L_{2}([0,1])$ and approximation $y_{{\mathbf{x}}_{m}}$ obtained from a random
sample ${\mathbf{x}}$ in Lemma 4. For this purpose, we have the following
assumption: a graphon signal $y$ satisfies
$|y(a)-y(b)|\leq{\alpha_{2}}|a-b|,\forall a,b\in[0,1]$. We term a graphon
signal satisfying this property as $\alpha_{2}$-Lipschitz graphon signal.

######  Lemma 4.

Given an $\alpha_{2}$-Lipschitz graphon signal $y$ and a graphon signal
approximation $y_{{\mathbf{x}}_{m}}$ obtained from
${\mathbf{x}}_{m}\in{\mathbb{R}}^{m\times 1}$, we have

$\displaystyle\|y-y_{{\mathbf{x}}_{m}}\|_{2}\leq\frac{\alpha_{2}\Omega^{3/2}}{m^{3\zeta/2-1}}\;.$
(104)

###### Proof.

Note that

$\displaystyle\|y-y_{{\mathbf{x}}_{m}}\|_{2}$
$\displaystyle=\sum\limits_{{\mathcal{U}}_{i}}\|y-y_{{\mathbf{x}}_{m}}\|_{L_{2}[I_{i}]}\;,$
(105)
$\displaystyle=\sum\limits_{i=1}^{m}\Big{(}\int_{\rho_{i-1}}^{\rho_{i}}(y(u)-y_{{\mathbf{x}}_{m}}(u))^{2}{\sf
d}u\Big{)}^{\frac{1}{2}}\;,$ (106)

where we have $\rho_{0}=0$. Using the Lipschitz property of graphon signal and
$(\Omega,\zeta)$-property of ${\mathbf{C}}_{m}$, we have

$\displaystyle\|y-y_{{\mathbf{x}}_{m}}\|_{2}\leq
m\Big{(}\alpha_{2}^{2}\int_{0}^{\rho_{1}}u^{2}{\sf d}u\Big{)}^{\frac{1}{2}}$
$\displaystyle\leq\frac{\alpha_{2}\Omega^{3/2}}{m^{3\zeta/2-1}}\;.$ (107)

∎

Next, we state Proposition 4 from [39] that characterizes a bound on the
difference between eigenvalues from two graphons.

######  Lemma 5 (Proposition 4 from [39]).

Consider two graphons ${\mathbf{W}}$ and ${\mathbf{W}}^{\prime}$ with set of
eigenvalues $\\{\eta_{i}\\}_{i=1}^{\infty}$ and
$\\{\beta_{i}\\}_{i=1}^{\infty}$, respectively. Then, for all
$i\in{\mathbb{Z}}^{+}$ , we have

$\displaystyle|\eta_{i}-\beta_{i}|\leq\|T_{{\mathbf{W}}-{\mathbf{W}}^{\prime}}\|_{2}\leq\|{\mathbf{W}}-{\mathbf{W}}^{\prime}\|_{2}\;.$
(108)

Next, we leverage Lemmas 3, 4, and 5 to bound the difference between graphon
convolution $\Psi(y;{\mathbf{W}},{\mathcal{H}})$ and convolution by the
approximation
$\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$
realized from graph filter over ${\mathbf{C}}_{m}$. The subsequent line of
analysis is similar to that in [79]. For this purpose, we have the following
assumptions.

  * A1

The graphon ${\mathbf{W}}$ is $\alpha$-Lipschitz.

  * A2

The graphon signal $y$ is $\alpha_{2}$-Lipschitz.

  * A3

The covariance matrix ${\mathbf{C}}_{m}$ belongs to a convergent
$(\Omega,\zeta)$-dominant sequence of covariance matrices. Also, the
corresponding sequence of graphon approximations belongs to the family of
graphon ${\mathbf{W}}$.

  * A4

The frequency response is band-limited, such that, $|\tilde{h}(\eta)|=0$ for
$\eta\leq\eta_{\sf c}$. Furthermore, we assume that $m_{\sf c}$ largest
eigenvalues of graphon ${\mathbf{W}}$ in terms of magnitude satisfy
$|\eta|>\eta_{\sf c}$ and the set of such eigenvalues is denoted by
${\mathcal{C}}$. Also, the graphon filter is non-expanding and satisfies
$|\tilde{h}(\eta)|\leq 1,\forall\eta$ and
$|\tilde{h}(\eta_{i})-\tilde{h}(\hat{\eta}_{i})|\leq\alpha_{3}|\eta_{i}-\hat{\eta}_{i}|$.

In the following Lemma, we use the notations $\\{\hat{\eta}_{i}\\}$ and
$\\{\hat{\Gamma}_{i}\\}$ for the set of eigenvalues and eigenfunctions,
respectively, of ${\mathbf{W}}_{{\mathbf{C}}_{m}}$.

######  Lemma 6 (Transferability of Graphon Filters).

For a convolution $\Psi(y;{\mathbf{W}},{\mathcal{H}})$ and its approximation
$\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$,
under the assumptions A1-A4 and for $\|y\|_{2}\leq 1$, we have

$\displaystyle\|\Psi(y;{\mathbf{W}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}\leq\frac{\Omega^{3/2}}{m^{3\zeta/2-1}}\left(\alpha_{2}+\alpha\Big{[}\alpha_{3}+\frac{\pi
m_{\sf c}}{2\Delta_{c}}\Big{]}\right)\;,$ (109)

where $\Delta_{c}=\min_{i\neq
j;i,j\in{\mathcal{C}}}\\{|\eta_{i}-\hat{\eta}_{j}|\\}$.

###### Proof.

Note that we can rewrite
$\|\Psi(y;{\mathbf{W}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}$
as

$\displaystyle\|\Psi(y;{\mathbf{W}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}$
$\displaystyle=\|\Psi(y;{\mathbf{W}},{\mathcal{H}})-\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$
$\displaystyle\quad+\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}\;,$
(110)

and using triangle inequality, we have

$\displaystyle\|\Psi(y;{\mathbf{W}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}$
$\displaystyle\leq\underbrace{\|\Psi(y;{\mathbf{W}},{\mathcal{H}})-\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}}_{\text{Term
1}}$
$\displaystyle\quad+\underbrace{\|\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}}_{\text{Term
2}}\;.$ (111)

Next, we analyze Terms 1 and 2 from (E) separately.

Analysis of Term 1. Using the expansion of
$\Psi(y;{\mathbf{W}},{\mathcal{H}})$ and
$\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$, we have

$\displaystyle\|\Psi(y;{\mathbf{W}},{\mathcal{H}})-\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}=\Big{(}\int_{0}^{1}f^{2}(v){\sf
d}v\Big{)}^{1/2}\;,$ (112)

where

$\displaystyle
f(v)=\sum\limits_{i=\in{\mathbb{Z}}\backslash\\{0\\}}\Big{[}\tilde{h}(\eta_{i})\Gamma_{i}(v)\int_{0}^{1}y(u)\Gamma_{i}(u){\sf
d}u-\tilde{h}(\hat{\eta}_{i})\hat{\Gamma}_{i}(v)\int_{0}^{1}y(u)\hat{\Gamma}_{i}(u){\sf
d}u\Big{]}\;.$ (113)

By adding and subtracting
$\tilde{h}(\hat{\eta}_{i})\Gamma_{i}(v)\int_{0}^{1}y(u)\Gamma_{i}(u){\sf d}u$
in (113) and using the triangle inequality, we have

$\displaystyle\|\Psi(y;{\mathbf{W}},{\mathcal{H}})-\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}\leq\Big{(}\int_{0}^{1}f_{1}^{2}(v){\sf
d}v\Big{)}^{1/2}+\Big{(}\int_{0}^{1}f_{2}^{2}(v){\sf
d}v\Big{)}^{1/2}=\|f_{1}\|_{2}+\|f_{2}\|_{2}\;,$ (114)

where

$\displaystyle
f_{1}(v)=\sum\limits_{i\in{\mathbb{Z}}\backslash\\{0\\}}\Big{[}(\tilde{h}(\eta_{i})-\tilde{h}(\hat{\eta}_{i}))\Gamma_{i}(v)\int_{0}^{1}y(u)\Gamma_{i}(u){\sf
d}u\Big{]}\;,$ (115)

and

$\displaystyle
f_{2}(v)=\sum\limits_{i\in{\mathbb{Z}}\backslash\\{0\\}}\Big{[}\tilde{h}(\hat{\eta}_{i})\Gamma_{i}(v)\int_{0}^{1}y(u)(\Gamma_{i}(u)-\hat{\Gamma}_{i}(u)){\sf
d}u\Big{]}\;.$ (116)

Using the Lipschitz property of graphon filter, we have
$|\tilde{h}(\eta_{i})-\tilde{h}(\hat{\eta}_{i})|\leq\alpha_{3}|\eta_{i}-\hat{\eta}_{i}|$.
Therefore, using Lemma 5 and Lemma 3, we have

$\displaystyle\|f_{1}\|_{2}\leq\frac{\alpha_{3}\alpha\Omega^{3/2}}{m^{3\zeta/2-1}}\;.$
(117)

for any $y$ that satisfies $\|y\|_{2}\leq 1$. To analyze $\|f_{2}\|_{2}$, we
note that by Cauchy-Schwarz inequality, we have

$\displaystyle\|f_{2}\|_{2}$
$\displaystyle\leq\sum\limits_{i\in{\mathbb{Z}}\backslash\\{0\\}}|\tilde{h}(\hat{\eta}_{i})|\|\Gamma_{i}\|_{2}\|y(\Gamma_{i}-\hat{\Gamma}_{i})\|_{2}\;,$
(118)
$\displaystyle\leq\sum\limits_{i\in{\mathbb{Z}}\backslash\\{0\\}}|\tilde{h}(\hat{\eta}_{i})|\|\Gamma_{i}-\hat{\Gamma}_{i}\|_{2}\;,$
(119)

where (119) follows from (118), without loss of generality for $\|y\|_{2}=1$,
$\|\Gamma_{i}\|_{2}=1$ and another application of Cauchy-Schwarz inequality.
Next, we note that the integral operator $T_{{\mathbf{W}}}$, such that,
$(T_{{\mathbf{W}}}y)(v)=\int_{0}^{1}{\mathbf{W}}(u,v)y(u){\sf d}u$ is a self-
adjoint Hilbert-Schmidt operator and ${\mathbf{W}}$ admits the spectral
decomposition with $\\{\eta_{i}\\}$ as eigenvalues and $\\{\Gamma_{i}\\}$ as
eigensignals. Therefore, to analyze $\|\Gamma_{i}-\hat{\Gamma}_{i}\|_{2}$, we
note that $\Gamma_{i}$ is projection of operator $T_{{\mathbf{W}}}$ associated
with eigenvalue $\eta_{i}$ and $\hat{\Gamma}_{i}$ is projection of operator
$T_{{\mathbf{W}}_{{\mathbf{C}}_{m}}}$ associated with eigenvalue
$\hat{\eta}_{i}$. By dividing the spectrum of $T_{{\mathbf{W}}}$ as ${\sf
spec}(T_{{\mathbf{W}}})=\\{\eta_{i}\\}\cup\\{\eta_{j}\\}_{j\neq i}$ and that
of $T_{{\mathbf{W}}_{{\mathbf{C}}_{m}}}$ as ${\sf
spec}(T_{{\mathbf{W}}_{{\mathbf{C}}_{m}}})=\\{\hat{\eta}_{i}\\}\cup\\{\hat{\eta}_{j}\\}_{j\neq
i}$, we apply Proposition 2.3 from [118] to have

$\displaystyle\|\Gamma_{i}-\hat{\Gamma}_{i}\|_{2}\leq\frac{\pi}{2}\frac{\|T_{{\mathbf{W}}}-T_{{\mathbf{W}}_{{\mathbf{C}}_{m}}}\|_{2}}{d_{i}}\;,$
(120)

where $d_{i}>0$ is a constant that satisfies $|\eta_{i}-\hat{\eta}_{i+1}|\geq
d_{i}$, $|\eta_{i}-\hat{\eta}_{i-1}|\geq d_{i}$,
$|\eta_{i+1}-\hat{\eta}_{i}|\geq d_{i}$, and $|\eta_{i-1}-\hat{\eta}_{i}|\geq
d_{i}$. Using (120), Lemma 5 and Lemma 3 in (119), we have

$\displaystyle\|f_{2}\|_{2}\leq\frac{\pi\alpha\Omega^{3/2}}{2\Delta_{c}m^{3\zeta/2-1}}\sum\limits_{i\in{\mathbb{Z}}\backslash{0}}|\tilde{h}(\hat{\eta}_{i})|\;,$
(121)

where $\Delta_{c}=\min_{i}d_{i}$. Next, we note that under Assumption A4, the
frequency response of graphon filter is band-limited, i.e., we have
$\sum\limits_{i\in{\mathbb{Z}}\backslash{0}}|\tilde{h}(\hat{\eta}_{i})|\leq
m_{\sf c}$ when the frequency response is limited according to A4. We denote
the set of $m_{\sf c}$ largest eigenvalues (in terms of magnitude) of
${\mathbf{W}}$ by ${\cal C}$. In this scenario, we can rewrite (121) as

$\displaystyle\|f_{2}\|_{2}\leq\frac{\pi\alpha\Omega^{3/2}m_{\sf
c}}{2\Delta_{c}m^{3\zeta/2-1}}\;.$ (122)

Clearly, there is a trade-off between $m_{c}$ and $\zeta$ as we must have
$m_{c}<m^{3\zeta/2-1}$ and $\zeta>2/3$ for (122) to have decreasing behavior
in $m$. Equations 117 and 122 provide the upper bound on Term 1.

Analysis of Term 2. We can expand term 2 as

$\displaystyle\|\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}=\Big{(}\int_{0}^{1}g^{2}(v){\sf
d}v\Big{)}^{1/2}\;,$ (123)

where

$\displaystyle
g(v)=\sum\limits_{i=1}^{\infty}\Big{[}\tilde{h}(\eta_{i})\hat{\Gamma}_{i}(v)\int_{0}^{1}(y(u)-y_{{\mathbf{x}}_{m}}(u))\hat{\Gamma}_{i}(u){\sf
d}u\Big{]}\;.$ (124)

Therefore, using (124), we have

$\displaystyle\|\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}=\|\Psi(y-y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}\;.$
(125)

Note that for a frequency response that satisfies $\tilde{h}(\eta)<1$, the
graphon filter is non-expanding and therefore, we have

$\displaystyle\|\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}\leq\|y-y_{{\mathbf{x}}_{m}}\|_{2}\;.$
(126)

Using Lemma 4, we have

$\displaystyle\|\Psi(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})-\Psi(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})\|_{2}\leq\frac{\alpha_{2}\Omega^{3/2}}{m^{3\zeta/2-1}}\;.$
(127)

Therefore, by combining the upper bounds on Term 1 and Term 2 from (117),
(122), and (127), the proof of Lemma 6 is concluded. ∎

Lemma 6 establishes the transference between the graphon ${\mathbf{W}}$ and
the graphon approximation ${\mathbf{W}}_{{\mathbf{C}}_{m}}$ obtained from the
covariance matrix ${\mathbf{C}}_{m}$. We leverage the result in Lemma 6 to
establishing the transference for graphon neural networks in a similar
setting. We denote the $f$-th output for graphon neural network
$\tilde{\Psi}(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$ with $F$
outputs in the final layer by
$[\tilde{\Psi}(y;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})]_{f}$.

######  Lemma 7 (Transferability of Graphon Neural Networks).

Consider a graphon neural network
$\tilde{\Phi}(\cdot;{\mathbf{W}},{\mathcal{H}})$ with $L$ layers and $F$
outputs per layer and a VNN $\Phi(\cdot;{\mathbf{C}}_{m},{\mathcal{H}})$ with
graphon neural network representation as
$\tilde{\Phi}(\cdot;{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})$. If the
covariance matrix ${\mathbf{C}}_{m}$ belongs to a $(\Omega,\zeta)$-dominant
sequence of covariance matrices and its graphon approximation
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ belongs to a graphon family of
$\alpha$-Lipschitz graphon ${\mathbf{W}}$, then under the assumptions A1-A4,
for $\|y\|_{2}\leq 1$ and $2/3<\beta\leq 1$, we have

$\displaystyle\|[\tilde{\Phi}(y;{\mathbf{W}},{\mathcal{H}})]_{f}-[\tilde{\Phi}(y_{{\mathbf{x}}_{m}};{\mathbf{W}}_{{\mathbf{C}}_{m}},{\mathcal{H}})]_{f}\|_{2}\leq{LF^{L}}\left(\frac{\Omega^{3/2}}{m^{3\zeta/2-1}}\left[\alpha_{2}+\alpha\left[\alpha_{3}+\frac{\pi
m_{\sf c}}{2\Delta_{c}}\right]\right]\right)\;.$ (128)

The proof of Lemma 7 leverages Lemma 6 and accommodates the impact of multi-
layer VNN architecture. We refer the reader to equations (23)-(28) in [39] for
exact analytical steps. Finally, by applying the triangle inequality on (128),
we establish the transference between graphon neural network approximations
$\tilde{\Phi}(\cdot;{\mathbf{W}}_{{\mathbf{C}}_{m_{1}}},{\mathcal{H}})$ and
$\tilde{\Phi}(\cdot;{\mathbf{W}}_{{\mathbf{C}}_{m_{2}}},{\mathcal{H}})$ for
VNNs $\Phi(\cdot;{{\mathbf{C}}_{m_{1}}},{\mathcal{H}})$ and
$\Phi(\cdot;{{\mathbf{C}}_{m_{2}}},{\mathcal{H}})$, respectively, and the
proof of Theorem 3 is concluded.

##  Appendix F Convergence of covariance matrices from FTDC datasets

Cut metric allows the comparison of matrix representations of graphs of
different sizes [81]. Here, an estimated cut norm between a pair of covariance
matrices was evaluated using the cutnorm package in python (implementation
based on [119]). Here, we investigate the convergence of the covariance
matrices formed by datasets curated according to different scales of the
Schaefer’s atlas. For this purpose, we also consider the datasets curated
according to 200 and 400 parcellations of Schaefer’s atlas in addition to the
FTDC datasets. Figure 13a plots the cut norm evaluated for the series
$\\{\hat{\mathbf{C}}_{100},\hat{\mathbf{C}}_{200},\hat{\mathbf{C}}_{300},\hat{\mathbf{C}}_{400},\hat{\mathbf{C}}_{500}\\}$,
where $\hat{\mathbf{C}}_{m}$ is the anatomical covariance matrix constructed
from the cortical thickness features curated according to $m$ parcellations of
Schaefer’s atlas. The cut norm between the consecutive members of this series
reduced with an increase in the number of parcels, thus, consistent with the
properties of a Cauchy sequence in this metric.

Figure 13: Cut norm between covariance matrices. Panel a plots the estimated
cut norm between consecutive elements of the series of covariance matrices
$\\{\hat{\mathbf{C}}_{100},\hat{\mathbf{C}}_{200},\hat{\mathbf{C}}_{300},\hat{\mathbf{C}}_{400},\hat{\mathbf{C}}_{500}\\}$.
Panel b plots the estimated cut norm between the covariance matrices derived
from OASIS-3 dataset according to DK and DKT atlases and the covariance
matrices from FTDC datasets.

We also note that the distance between covariance matrices derived from
OASIS-3 dataset and those from FTDC datasets was significantly greater than
the distances between the covariance matrices associated with different
resolutions of Schaefer’s atlas (Fig. 13b). In Fig. 13b, the covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf OASIS}$ was estimated from the cortical thickness
features curated according to DKT atlas for HC group. The
$\hat{\mathbf{C}}_{68}^{\sf OASIS}$ was estimated from the cortical thickness
features curated according to Desikan-Killiany atlas [120] for the same group.
According to Theorem 3, the VNN performance is transferable across datasets
whose covariance matrices are part of a converging sequence. Hence, we
expected the VNNs to be transferable between the cortical thickness datasets
curated according to different scales of Schaefer’s atlas. However, VNN
transferability was not guaranteed between datasets organized according to
Schaefer’s atlas and those according to DKT atlas.

##  Appendix G Summary of the performance of VNNs on the test sets

Here, we provide the performance of nominal VNNs on the test set for different
sets of VNNs that were trained on FTDC100, FTDC500, and OASIS-3 datasets. The
test set was completely unseen during the training procedure.

OASIS-3: The MAE across $100$ nominal models on the test set was $5.87\pm
0.177$ years and Pearson’s correlation was $0.425\pm 0.006$.

FTDC100: The MAE across $100$ nominal models on the test set was $3.72\pm
0.22$ years and Pearson’s correlation was $0.78\pm 0.01$.

FTDC300: The MAE across $100$ nominal models on the test set was $3.93\pm
0.34$ years and Pearson’s correlation was $0.75\pm 0.008$.

FTDC500: The MAE across $100$ nominal models on the test set was $4.07\pm
0.36$ years and Pearson’s correlation was $0.72\pm 0.007$.

##  Appendix H Illustration of regional residual analysis from VNN model
outputs

In this section, we demonstrate the regional analysis described in Section
4.4.1 for a VNN model that was trained to predict chronological age for HC
group in OASIS-3 dataset. All mathematical notations referred to in this
section are borrowed from Section 4.4.1. Note that no further training was
performed for this VNN model to evaluate brain age or regional residuals.

The covariance matrix in this VNN model was replaced with ${{\mathbf{C}}}^{\sf
AD+}_{148}$ and a combined group of HC and AD+ individuals were processed.
Thus, for each individual $i$ in this combined dataset, we obtained a age
prediction $\hat{y}_{i}$ and a vector of residuals ${\mathbf{r}}_{i}$. The
size of residual vector ${\mathbf{r}}_{i}$ was $148\times 1$ and hence, each
element of ${\mathbf{r}}_{i}$ corresponded to a distinct brain region as
defined by the DKT brain atlas with 148 parcellations. By evaluating the
vector of residuals ${\mathbf{r}}_{i}$ for every individual in the combined
dataset, we were able to form a population of residual vectors from HC group
(referred to as ${\mathbf{r}}_{\sf HC}$) and AD+ group (referred to as
${\mathbf{r}}_{\sf AD+}$). The elements of these residual vectors are referred
to as regional residuals throughout the paper.

Each dimension of these residual vectors was investigated for group
differences between HC and AD+ groups via ANOVA as described in Section 4.4.1.
Thus, for every VNN model, we eventually performed $m=148$ number of ANOVA
tests and evaluated the brain regions for significance in group differences in
their respective residuals. The significance of group differences between the
distributions of regional residuals for HC and AD+ groups corresponding to a
brain region was determined after correcting the $p$-values of ANOVA test for
multiple comparisons via Bonferroni correction (Bonferroni corrected $p$-value
$<0.05$). The group differences were additionally investigated for
significance at an uncorrected level using ANCOVA with age and sex as
covariates.

Figure 14 illustrates the results obtained via ANOVA in this context. The
brain regions are shaded according to the magnitude of the F-values from ANOVA
and all shaded regions were significant according to the criteria provided in
Section 4.4.1. The box plots for various brain regions show that the residuals
were significantly elevated in AD+ group as compared to HC. Although many
brain regions were identified as significant, we note that there was a clear
variability in significance, with some brain regions having more significant
group difference than others. We had 100 trained VNN models for the OASIS-3
dataset and performed similar analyses for each of them. Further, we counted
the number of models for which the above described analysis yielded a brain
region to be significant. A brain region with robust group difference in its
regional residual distribution in HC vs AD+ was expected to be more frequently
labeled as significant by the VNN models. The results of this robustness
analyses on the OASIS-3 dataset are shown in Fig. 5.3.1b.

Figure 15 illustrates the results obtained by analyzing the regional residuals
that were derived by transferring a VNN model from FTDC100 dataset to OASIS-3
dataset. Here, we observed that brain regions similar to those in Fig. 14
exhibited significantly elevated regional residuals in the AD+ group (albeit
with lower significance in terms of $F$-values of the ANOVA test). Thus, Fig.
15 provides the evidence that VNN trained on FTDC100 dataset can transfer the
interpretability from FTDC100 dataset to OASIS-3 dataset despite the lack of
transferability of quantitative performance.

Figure 14: Results depicting the brain regions with significantly elevated
regional residuals for AD+ group with respect to HC group in OASIS-3. The
results here were derived by a VNN model that was trained as a regression
model to predict chronological age from cortical thickness data for HC group
in OASIS-3. The contrasts of the shaded regions are proportional to the
F-values obtained from ANOVA and all shaded regions exhibited elevated
regional residuals in AD+ with respect to HC. Figure 15: Results depicting the
brain regions with significantly elevated regional residuals for AD+ group
with respect to HC group in OASIS-3. The results here were derived by a VNN
model that was transferred from FTDC100 to process OASIS-3 and had been
trained as a regression model to predict chronological age from cortical
thickness data for healthy controls in the FTDC100 dataset. The contrasts of
the shaded regions are proportional to the F-values obtained from ANOVA and
all shaded regions exhibited significant elevated regional residuals in AD+
with respect to HC (Bonferroni-corrected $p$-value smaller than $0.05$).

##  Appendix I Additional details on brain age prediction in OASIS-3

In this section, we provide additional figures and discussions pertaining to
the results for interpretable brain age prediction in Fig. 5.3.1. Figure 16
displays the distributions of chronological age for AD+ and HC groups in
OASIS-3 dataset.

Figure 16: Distribution of chronological age in AD+ and HC groups.

Since VNNs were trained for the regression task, a VNN processed cortical
thickness data and provided an estimate for the chronological age for each
individual. Since we trained 100 VNN models on different permutations of the
training set in OASIS-3, we use the mean of all VNN estimates as the VNN
prediction for an individual. This VNN prediction is further leveraged to form
brain age estimates and $\Delta$-Age metrics. Figure 17a displays the plot for
VNN predictions versus chronological age (ground truth) for the complete HC
group. The Pearson’s correlation between VNN prediction and chronological age
(ground truth) for HC group was $0.43$ which was similar to that reported in
Table 1. This observation implied that there was no degradation in the
performance of VNNs on the task of predicting chronological age for HC group
despite replacing the anatomical covariance matrix only from HC group with
that from the combined HC and AD+ group ($\hat{\mathbf{C}}_{148}^{\sf AD+}$).
However, VNN outputs clearly under-estimated the chronological age for older
individuals and over-estimated the chronological age for individuals on the
younger end of the age distribution for HC group.

Figure 17b displays the plot for VNN predictions versus chronological age
(ground truth) for the complete AD+ group. The Pearson’s correlation between
VNN prediction and chronological age (ground truth) for AD+ group was $0.28$.
We further note that the VNN architecture and our analysis of regional
residuals helped quantify the contribution of each brain region to a data
point in Fig. 17a and Fig. 17b. Hence, the scatter plot in Fig. 17b could be
affected by larger contributions of certain brain regions for AD+ group
relative to the HC group.

Figure 17c illustrates the box plots of residuals evaluated by the difference
between VNN predictions and chronological age for HC and AD+ groups. Figure
17c suggests that the chronological age for AD+ group was underestimated as
compared to that for HC group. This observation was also expected since AD+
group is significantly older than the HC group. However, we expect that the
robust elevated regional residuals from brain regions in Fig. 5.3.1b mitigated
the under-estimation effect due to higher age of AD+ group to some extent.

Figures 17d-f display the results after age-bias correction is applied to the
VNN outputs. As expected, the brain age for HC group in Fig. 17d is largely
concentrated around the line of equality ($x=y$ line). In contrast, the brain
age for AD+ group in Fig. 17e is concentrated above the line of equality.
These effects manifest into the box plots for $\Delta$-Age in Fig. 17f where
we observe the AD+ group to have elevated $\Delta$-Age as compared to HC
group.

VNN architecture facilitated isolation of the effects of accelerated aging
before age-bias correction was applied. Hence, the transformation of VNN
outputs to brain age from Fig. 17a-c to Fig. 17d-f was not surprising.
However, such insights may be infeasible for machine learning approaches that
lack transparency and hence, the impact of deviations due to neurodegeneration
from the healthy control population cannot be interpreted or isolated. In this
context, if the learning model was a black box, Fig. 17a-c may appear to be
counter-intuitive to the goal of detecting accelerated aging in the AD+ group
and the effect of age-bias correction can be unclear, thus, leading to several
criticisms [18].

Figure 17: Supplementary figures to Fig. 5.3.1. Panel a displays the plot of
VNN prediction versus chronological age for HC group. VNN predictions were
obtained as the average of the outputs of $100$ nominal VNNs that were trained
on OASIS-3 and operated on the anatomical covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$. Panel b displays the results similar to
that in panel a for the AD+ group. The solid line in panels a and b is the
least squares line. Panel c includes the boxplots for residuals derived from
the difference between VNN predictions and chronological age for HC and AD+
groups. Panel d and e display the plots for brain age versus chronological age
for HC and AD+ groups, respectively. The solid line in panels d and e is the
identity line. Panel f displays the box plots for $\Delta$-Age in HC and AD+
groups.

##  Appendix J Additional results on associations between regional residuals
derived from VNNs and eigenvectors of $\hat{\mathbf{C}}_{148}^{\sf AD+}$

Figure 18 displays the results for the associations between the regional
contributions to the VNN outputs ($\bar{\mathbf{p}}_{\sf HC}$) and the
eigenvectors of the anatomical covariance matrix $\hat{\mathbf{C}}_{148}$ for
VNNs trained on FTDC300 or FTDC500 datasets and transferred to process the
data from HC group in OASIS-3 dataset. The results in Fig. 18 show that the
first, third, and fourth eigenvectors of $\hat{\mathbf{C}}_{148}$ had the
three highest association with the regional contributions derived from VNNs
trained on FTDC300 or FTDC500 datasets, which was consistent with the results
derived for VNNs trained on OASIS-3 or FTDC100 datasets in Fig. 6. Besides
these three eigenvectors, there were a smaller number of eigenvectors for
which $|\\!<\\!\bar{\mathbf{p}}_{\sf HC},{\mathbf{v}}_{i}\\!>\\!|$ had a
coefficient of variation less than $30\%$ across the HC group for the results
derived from VNNs trained on FTDC300 or FTDC500 datasets in Fig. 18 as
compared to those derived from VNNs trained on OASIS-3 dataset in Fig. 6a.

Figure 18:  Associations between eigenvectors of $\hat{\mathbf{C}}_{148}$ and
regional contributions to the VNN outputs ($\bar{\mathbf{p}}_{\sf HC}$) for
VNNs transferred from FTDC300 or FTDC500 datasets to OASIS-3 dataset. Panel a
illustrates a bar plot for $|\\!<\\!\bar{\mathbf{p}}_{\sf
HC},{\mathbf{v}}_{i}\\!>\\!|$ for $i\in\\{1,\dots,30\\}$, where
${\mathbf{v}}_{i}$ is the $i$-th eigenvector (principal component) of
covariance matrix $\hat{\mathbf{C}}_{148}$ and associated with $i$-largest
eigenvalue in terms of magnitude and the vectors of regional contributions,
$\bar{\mathbf{p}}_{\sf HC}$ were obtained by VNNs that were trained on FTDC300
dataset and transferred to process the data from HC group in OASIS-3 dataset.
The inner product results for eigenvectors with coefficient of variation
greater than $30\%$ across the HC group of OASIS-3 were excluded (and hence,
their contributions set as 0). For every individual in HC group, the
associations between their corresponding vector of regional contributions,
$\bar{\mathbf{p}}_{\sf HC}$ and eigenvectors of $\hat{\mathbf{C}}_{148}$ were
evaluated over $100$ nominal VNN models. Panel b illustrates similar results
as in Panel a for the VNNs that were trained on FTDC500 dataset and
transferred to process the data from HC group in OASIS-3.

Figure 19a displays the associations between the eigenvectors of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ and regional residuals for AD+ group that
were derived from VNNs trained on different datasets. For VNNs trained on
OASIS-3, the three largest associations between regional residuals and
eigenvectors were observed for the first, third, and fourth eigenvectors of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$. Interestingly, this phenomenon was
retained when regional residuals were evaluated using the final layer outputs
of the VNNs that had been transferred from FTDC datasets to OASIS-3 dataset.
However, we also observed that the variance in the inner products was larger
for VNNs trained on FTDC datasets as compared to those trained on OASIS-3
dataset. As a consequence, in Fig. 19a, the number of eigenvectors whose
associations with regional residuals had a coefficient of variation smaller
than $30\%$ diminished considerably for VNNs trained on FTDC300 or FTDC500
datasets.

Figure 19b displays the regional profiles corresponding to the robustness of
the elevated regional residuals in AD+ group with respect to HC group in
OASIS-3 derived from different sets of 100 VNNs. Considering the regional
profile derived from VNNs trained on OASIS-3 dataset as the baseline, the
regional profile derived from the analyses of regional residuals obtained from
VNNs trained on FTDC100 was the most consistent. The regional profiles derived
from the VNNs trained on FTDC300 or FTDC500 datasets retained the most
robustness in parahippocampal and subcallosal regions. These observations were
expected from the results in Fig. 19a as the third and fourth eigenvectors of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ spanned the parahippocampal and subcallosal
regions and were highly correlated with the regional residuals derived from
VNNs trained on FTDC datasets. The results in Fig. 19 further corroborate our
conclusion that the regional profiles characteristic of the contributors to
elevated $\Delta$-Age were dependent on the ability of VNNs to exploit
specific eigenvectors of the anatomical covariance matrix.

Figure 19:  Associations between eigenvectors of $\hat{\mathbf{C}}_{148}^{\sf
AD+}$ and regional residuals for AD+ group; and regional profiles
corresponding to elevated regional residuals in AD+ group. Panel a illustrates
the bar plots for the inner products $|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{i}\\!>\\!|$ between the first $50$ eigenvectors of the
anatomical covariance matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$ and the
normalized regional residuals evaluated across the AD+ group, where the
regional residuals were derived from the final layer outputs of different sets
of VNNs. The results for eigenvectors for which $|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{i}\\!>\\!|$ had a coefficient of variation across the AD+
group greater than $30\%$ were excluded (and hence, have entry $0$ on the bar
plot). Panel b displays the regional profiles for the brain regions associated
with robust and significant elevated regional residuals in AD+ group as
compared to HC group in OASIS-3. The regional profiles obtained from the
analyses of regional residuals that were derived from VNNs trained on OASIS-3
and FTDC datasets are included.

##  Appendix K Transferability of VNNs for evaluating regional profiles and
$\Delta$-Age in multi-scale FTDC datasets

In this section, we provide a proof of concept that transferability between
different scales of Schaefer’s atlas for chronological age prediction lead to
consistent regional profiles for $\Delta$-Age determined by regional residuals
when VNNs were transferred across cortical thickness datasets corresponding to
different resolutions of Schaefer’s atlas. The dataset used for this
experiment is described below.

Dataset: This dataset consisted of $105$ healthy controls from the FTDC
dataset and $67$ individuals with mild cognitive impairment or Alzheimer’s
disease diagnosis (AD+; age = $68.52\pm 9.29$ years, $28$ females). The
cortical thickness data were available at $100$, $300$, and $500$ resolutions
of the Schaefer’s atlas.

The regional residuals were derived for the dataset above using a VNN model
that had been trained to predict chronological age in FTDC100 dataset. The
group differences in the regional residuals for HC and the combined AD+ group
were evaluated using ANOVA for cortical thickness at $100$, $300$, and $500$
resolutions of Schaefer’s atlas. The brain regions with robust, significantly
elevated regional residuals in AD+ with respect to HC are projected on the
brain template for $100$, $300$, and $500$ resolutions in Fig. 20. In Fig. 20,
the isolated brain regions were concentrated in similar regions across all
resolutions of Schaefer’s atlas and were consistent with the brain regions in
Fig. 5.3.1b (with the exception of the precuneus and superior parietal regions
in both hemispheres).

Subsequently, $\Delta$-Age was evaluated for all individuals. The $\Delta$-Age
evaluated from the cortical thickness dataset with $100$ features for AD group
was $3.67\pm 3.73$ years and for HC was $0\pm 2.06$ years. Figure 20b displays
the box plots for $\Delta$-Age in HC and AD+ groups for FTDC100 dataset as
well as those obtained after transferring the VNNs from FTDC100 to FTDC300 and
FTDC500 datasets. Due to the transferability of VNNs across different scales
of Schaefer’s atlas, $\Delta$-Age profiles and brain age versus chronological
age plots are consistent across $100$, $300$, and $500$ resolution data even
when the VNNs were trained only on the FTDC100 dataset.

Figure 20:  Brain age prediction across datasets curated according to multiple
scales of Schaefer’s atlas. Panel a illustrates the regional profiles
consisting of brain regions with robust, elevated regional residuals in the
combined AD+ group with respect to the HC group. The VNNs were trained to
predict chronological age on FTDC100 and the robustness was evaluated over
$100$ nominal VNN models. The regional profiles were obtained for the datasets
with $300$ features and $500$ features after transferring the VNNs from
FTDC100 to FTDC300 and FTDC500. Panel b displays the box plots for
$\Delta$-Age corresponding to the regional profiles in Panel a. Panel c plots
brain age versus chronological age for datasets with $100$, $300$, and $500$
cortical thickness features according to different scales of Schaefer’s atlas.

##  Appendix L Regional profiles corresponding to elevated regional residuals
in AD+ group are stable to the composition of data used to estimate anatomical
covariance matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$

Recall that $\Delta$-Age and associated regional profiles were evaluated using
VNNs that operated upon a composite anatomical covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$. We next checked whether the results
derived from VNNs were stable to the changes in composition of the combined HC
and AD+ group used to estimate the anatomical covariance matrix. For this
purpose, we performed two sets of experiments. In the first set, we included
the whole HC group and gradually varied the number of individuals from the AD+
group to be included to form the estimate $\hat{\mathbf{C}}_{148}^{\sf AD+}$.
Figure Ma includes the results obtained from a randomly selected VNN model
corresponding to anatomical covariance matrix formed by different combinations
of the individuals from HC and AD+ groups. The results in Fig. Ma display the
F-values of ANOVA when the regional residuals from AD+ group are higher than
HC group (Bonferroni corrected $p$-value $<0.05$). The results obtained by the
VNN when it used $\hat{\mathbf{C}}_{148}^{\sf AD+}$ that was estimated from
all $652$ HC individuals and $209$ AD+ individuals displays the most robust
group differences in regions that constitute the regional profile for group
differences in $\Delta$-Age between AD+ and HC in Fig. 5.3.1b. When the
covariance matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$ was perturbed by using
smaller number of individuals from the AD+ group to estimate it, we observed
that the obtained results were preserved till exclusion of about $134$ AD+
individuals. The results obtained via VNN when operating on
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ estimated from $50$ or smaller number of
AD+ individuals and all $652$ HC individuals became noticeably less
significant in terms of $F$-values in general, with the $F$-values in
subcallosal, parahippocampal, temporal pole, and medial temporal lobe regions
in both hemispheres among those affected noticeably. Hence, the removal of AD
related information from the anatomical covariance matrix led to the results
in brain regions characteristic of AD becoming less significant.

Figure Mb illustrates the results obtained for a similar experiment as above,
with the difference that the regional residuals were evaluated for the VNN
when the anatomical covariance matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$ was
perturbed by reducing the number of individuals from the HC group used to
estimate it. Using the result obtained for $\hat{\mathbf{C}}_{148}^{\sf AD+}$
estimated from $652$ HC individuals and $209$ AD+ individuals in Fig. Ma as
the baseline, the results pertaining to ANOVA between regional residuals for
AD+ and HC groups (with AD+ elevated as compared to HC) remained consistent as
long as $100$ or more individuals from HC group were included in forming
$\hat{\mathbf{C}}_{148}^{\sf AD+}$. With less than 100 number of HC
individuals included in $\hat{\mathbf{C}}_{148}^{\sf AD+}$, the results became
noticeably less significant in precuneus and supramarginal regions in the left
hemisphere.

In summary, the group differences observed between the regional residuals for
AD+ and HC groups in OASIS-3 dataset were robust to perturbations in the
covariance matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$ when it was perturbed
from the baseline by using a different combination of HC and AD+ individuals
to estimate it. However, (nearly) complete exclusion of HC or AD+ groups from
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ resulted in loss of significance of the
elevation in regional residuals in AD+ for various regions, including
bilateral parahippocampal and temporal pole regions, and precuneus and
supramarginal regions in the left hemisphere. Thus, both HC and AD+ groups
were relevant to the anatomical covariance matrix $\hat{\mathbf{C}}_{148}^{\sf
AD+}$ that resulted in regional profiles in Fig. 5.3.1 but the regional
profiles in Fig. 5.3.1 were not overfit on the combination of individuals from
HC and AD+ used to estimate $\hat{\mathbf{C}}_{148}^{\sf AD+}$.

##  Appendix M Adaptive readouts may penalize the interpretability of regional
residuals and $\Delta$-Age

Thus far, we have focused on VNNs that operate with a non-adaptive readout
(unweighted average) function. However, it is expected that the performance of
the VNNs on their original task of chronological age prediction could be
improved significantly with the help of an adaptive readout function. Our
experiments showed that this was indeed the case. If a single-layer fully
connected perceptron consisting of $10$ neurons was added to the VNNs with the
same architecture as the ones that were trained on OASIS-3 dataset, we could
improve the performance on the chronological prediction task. For 100 VNNs
with adaptive readout that were trained on random permutations of the training
data, the median MAE for the HC group was $4.64$ years, which was
significantly smaller than the MAE achieved by VNNs with non-adaptive readouts
(Table 1). Among the 100 VNN models with adaptive readouts, we analyzed the
regional residuals for one VNN model with adaptive readout that had the best
performance on chronological age prediction in HC group (test set: MAE =
$4.17$ years, Pearson’s correlation between prediction and ground truth =
$0.73$; complete HC group: MAE = $4.26$ years, Pearson’s correlation between
prediction and ground truth = $0.725$). Our regional residuals revealed no
significant difference between the regional residuals for AD+ group and HC
group. This observation suggested that VNN lost its interpretability due to
the addition of adaptive readout function. Moreover, we also observed a
diminished gap between $\Delta$-Age for AD+ and HC groups determined using
this VNN model ($\Delta$-Age for AD group: $1.58\pm 4.67$ years, $\Delta$-Age
for HC group: $0\pm 3.45$ years, Cohen’s $d$ = 0.384). The findings discussed
here suggest that boosting the performance on chronological age prediction
task by using an adaptive readout function may penalize the interpretability
offered by VNNs with non-adaptive readouts and also diminish the $\Delta$-Age
gap between AD+ and HC groups.

(preceding page): Stability to perturbations in the anatomical covariance
matrix for group differences between AD+ and HC groups observed in regional
residuals. For a VNN model that was trained to predict chronological age for
HC group in OASIS-3 dataset, the regional residuals were first determined
using the anatomical covariance matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$
formed by the cortical thickness data of complete OASIS-3 dataset (i.e., 652
HC individuals and 209 individuals in the AD+ group. The group differences in
regional residuals between AD+ and HC group were investigated according to the
procedure in subsection 4.4.1 in the Methods section. In the procedure
described therein, we evaluated the F-values for the ANOVA test between
regional residuals for AD+ group and HC group. The F-values for the regional
residuals that were elevated in AD+ group with respect to HC group are
projected on the brain template in panel a. The stability of the group
differences to perturbations in $\hat{\mathbf{C}}_{148}^{\sf AD+}$ was further
investigated by varying the composition of cortical thickness data from AD+
and HC groups used to estimate $\hat{\mathbf{C}}_{148}^{\sf AD+}$. Panel a
displays the results obtained via analysis of regional residuals by VNNs that
processes the cortical thickness data from complete OASIS-3 dataset over the
anatomical covariance matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$ estimated from
$652$ HC individuals and varying number of individuals from the AD+group. The
scenarios in which $\hat{\mathbf{C}}_{148}^{\sf AD+}$ was estimated using
cortical thickness data from $652$ HC individuals and no AD+ individuals
formed one extreme of these experiments. Panel b illustrates the results of
similar experiments as in Panel a with the difference that the anatomical
covariance matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$ was estimated using all
$209$ individuals in the AD+ group and varying number of individuals from the
HC group. The results corresponding to $\hat{\mathbf{C}}_{148}^{\sf AD+}$ that
was estimated using $209$ AD+ individuals and $652$ HC individuals formed the
baseline for comparison for all scenarios in Panels a and b.

## References

  * [1] S. Sihag, G. Mateos, C. McMillan, and A. Ribeiro, “coVariance neural networks,” in _Proc. Conference on Neural Information Processing Systems_ , Nov. 2022.
  * [2] ——, “Predicting brain age using transferable coVariance neural networks,” in _Proc. IEEE International Conference on Acoustics, Speech, and Signal Processing_ , Jun. 2023.
  * [3] P. Liu, “A survey of remote-sensing big data,” _Frontiers in Environmental Science_ , vol. 3, p. 45, 2015.
  * [4] B. H. Brinkmann, M. R. Bower, K. A. Stengel, G. A. Worrell, and M. Stead, “Large-scale electrophysiology: acquisition, compression, encryption, and storage of big data,” _Journal of Neuroscience Methods_ , vol. 180, no. 1, pp. 185–192, 2009.
  * [5] J. Cai, J. Luo, S. Wang, and S. Yang, “Feature selection in machine learning: A new perspective,” _Neurocomputing_ , vol. 300, pp. 70–79, 2018.
  * [6] A. Y. Ng, “Feature selection, l 1 vs. l 2 regularization, and rotational invariance,” in _Proceedings of the twenty-first international conference on Machine learning_ , 2004, p. 78.
  * [7] Y. Yang, Y. Zhu, H. Cui, X. Kan, L. He, Y. Guo, and C. Yang, “Data-efficient brain connectome analysis via multi-task meta-learning,” _arXiv preprint arXiv:2206.04486_ , 2022.
  * [8] R. D. Markello, J. Y. Hansen, Z.-Q. Liu, V. Bazinet, G. Shafiei, L. E. Suárez, N. Blostein, J. Seidlitz, S. Baillet, T. D. Satterthwaite _et al._ , “Neuromaps: structural and functional interpretation of brain maps,” _Nature Methods_ , vol. 19, no. 11, pp. 1472–1479, 2022.
  * [9] A. Schaefer, R. Kong, E. M. Gordon, T. O. Laumann, X.-N. Zuo, A. J. Holmes, S. B. Eickhoff, and B. T. Yeo, “Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI,” _Cerebral Cortex_ , vol. 28, no. 9, pp. 3095–3114, 2018.
  * [10] P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C. J. Honey, V. J. Wedeen, and O. Sporns, “Mapping the structural core of human cerebral cortex,” _PLoS Biology_ , vol. 6, no. 7, p. e159, 2008.
  * [11] M. Habes, D. Janowitz, G. Erus, J. Toledo, S. Resnick, J. Doshi, S. Van der Auwera, K. Wittfeld, K. Hegenscheid, N. Hosten _et al._ , “Advanced brain aging: Relationship with epidemiologic and genetic risk factors, and overlap with Alzheimer disease atrophy patterns,” _Translational Psychiatry_ , vol. 6, no. 4, pp. e775–e775, 2016.
  * [12] K. Franke and C. Gaser, “Longitudinal changes in individual brainage in healthy aging, mild cognitive impairment, and Alzheimer’s disease.” _GeroPsych: The Journal of Gerontopsychology and Geriatric Psychiatry_ , vol. 25, no. 4, p. 235, 2012.
  * [13] T. Hajek, K. Franke, M. Kolenic, J. Capkova, M. Matejka, L. Propper, R. Uher, P. Stopkova, T. Novak, T. Paus _et al._ , “Brain age in early stages of bipolar disorders or schizophrenia,” _Schizophrenia Bulletin_ , vol. 45, no. 1, pp. 190–198, 2019.
  * [14] J. H. Cole and K. Franke, “Predicting age using neuroimaging: Innovative brain ageing biomarkers,” _Trends in Neurosciences_ , vol. 40, no. 12, pp. 681–690, 2017.
  * [15] C. Yin, P. Imms, M. Cheng, A. Amgalan, N. F. Chowdhury, R. J. Massett, N. N. Chaudhari, X. Chen, P. M. Thompson, P. Bogdan _et al._ , “Anatomically interpretable deep learning of brain age captures domain-specific cognitive impairment,” _Proceedings of the National Academy of Sciences_ , vol. 120, no. 2, p. e2214634120, 2023.
  * [16] V. M. Bashyam, G. Erus, J. Doshi, M. Habes, I. M. Nasrallah, M. Truelove-Hill, D. Srinivasan, L. Mamourian, R. Pomponio, Y. Fan _et al._ , “MRI signatures of brain age and disease over the lifespan based on a deep brain network and 14 468 individuals worldwide,” _Brain_ , vol. 143, no. 7, pp. 2312–2324, 2020.
  * [17] B. Couvy-Duchesne, J. Faouzi, B. Martin, E. Thibeau-Sutre, A. Wild, M. Ansart, S. Durrleman, D. Dormont, N. Burgos, and O. Colliot, “Ensemble learning of convolutional neural network, support vector machine, and best linear unbiased predictor for brain age prediction: Aramis contribution to the predictive analytics competition 2019 challenge,” _Frontiers in Psychiatry_ , vol. 11, p. 593336, 2020.
  * [18] E. R. Butler, A. Chen, R. Ramadan, T. T. Le, K. Ruparel, T. M. Moore, T. D. Satterthwaite, F. Zhang, H. Shou, R. C. Gur _et al._ , “Pitfalls in brain age analyses,” Wiley Online Library, Tech. Rep., 2021.
  * [19] J. Zhou, G. Cui, S. Hu, Z. Zhang, C. Yang, Z. Liu, L. Wang, C. Li, and M. Sun, “Graph neural networks: A review of methods and applications,” _AI Open_ , vol. 1, pp. 57–81, 2020.
  * [20] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini, “The graph neural network model,” _IEEE Transactions on Neural Networks_ , vol. 20, no. 1, pp. 61–80, 2008.
  * [21] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and P. S Yu, “A comprehensive survey on graph neural networks,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 32, no. 1, pp. 4–24, 2020.
  * [22] S. Sahu, A. Mhedhbi, S. Salihoglu, J. Lin, and M. T. Özsu, “The ubiquity of large graphs and surprising challenges of graph processing,” _Proceedings of the VLDB Endowment_ , vol. 11, no. 4, pp. 420–431, 2017.
  * [23] V. Quintana, A. Simoes-Costa, and A. Mandel, “Power system topological observability using a direct graph-theoretic approach,” _IEEE Transactions on Power Apparatus and Systems_ , no. 3, pp. 617–626, 1982.
  * [24] S. A. Myers, A. Sharma, P. Gupta, and J. Lin, “Information network or social network? The structure of the Twitter follow graph,” in _Proceedings of the 23rd International Conference on World Wide Web_ , 2014, pp. 493–498.
  * [25] A. Fout, J. Byrd, B. Shariat, and A. Ben-Hur, “Protein interface prediction using graph convolutional networks,” _Advances in Neural Information Processing Systems_ , vol. 30, 2017.
  * [26] O. Sporns, “Graph theory methods: applications in brain networks,” _Dialogues in Clinical Neuroscience_ , 2022.
  * [27] A. Sanchez-Gonzalez, N. Heess, J. T. Springenberg, J. Merel, M. Riedmiller, R. Hadsell, and P. Battaglia, “Graph networks as learnable physics engines for inference and control,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 4470–4479.
  * [28] S. Abadal, A. Jain, R. Guirado, J. López-Alonso, and E. Alarcón, “Computing graph neural networks: A survey from algorithms to accelerators,” _ACM Computing Surveys (CSUR)_ , vol. 54, no. 9, pp. 1–38, 2021.
  * [29] S. Zhang, H. Tong, J. Xu, and R. Maciejewski, “Graph convolutional networks: A comprehensive review,” _Computational Social Networks_ , vol. 6, no. 1, pp. 1–23, 2019.
  * [30] T. N. Kipf and M. Welling, “Variational graph auto-encoders,” _arXiv preprint arXiv:1611.07308_ , 2016.
  * [31] Y. Li, R. Zemel, M. Brockschmidt, and D. Tarlow, “Gated graph sequence neural networks,” in _Proceedings of ICLR’16_ , 2016.
  * [32] J. Bruna, W. Zaremba, A. Szlam, and Y. LeCun, “Spectral networks and locally connected networks on graphs,” _arXiv preprint arXiv:1312.6203_ , 2013.
  * [33] D. K. Hammond, P. Vandergheynst, and R. Gribonval, “Wavelets on graphs via spectral graph theory,” _Applied and Computational Harmonic Analysis_ , vol. 30, no. 2, pp. 129–150, 2011.
  * [34] F. Gama, E. Isufi, G. Leus, and A. Ribeiro, “Graphs, convolutions, and neural networks: From graph filters to graph neural networks,” _IEEE Signal Processing Magazine_ , vol. 37, no. 6, pp. 128–138, 2020.
  * [35] J. Atwood and D. Towsley, “Diffusion-convolutional neural networks,” _Advances in neural information processing systems_ , vol. 29, 2016.
  * [36] S. Verma and Z.-L. Zhang, “Stability and generalization of graph convolutional neural networks,” in _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2019, pp. 1539–1548.
  * [37] N. Keriven, A. Bietti, and S. Vaiter, “Convergence and stability of graph convolutional networks on large random graphs,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 21 512–21 523, 2020.
  * [38] F. Gama, J. Bruna, and A. Ribeiro, “Stability properties of graph neural networks,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 5680–5695, 2020.
  * [39] L. Ruiz, L. Chamon, and A. Ribeiro, “Graphon neural networks and the transferability of graph neural networks,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 1702–1712, 2020.
  * [40] A. C. Evans, “Networks of anatomical covariance,” _Neuroimage_ , vol. 80, pp. 489–504, 2013.
  * [41] H. Murase and S. K. Nayar, “Illumination planning for object recognition using parametric eigenspaces,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 16, no. 12, pp. 1219–1227, 1994.
  * [42] F. De la Torre and M. J. Black, “Robust principal component analysis for computer vision,” in _Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001_ , vol. 1. IEEE, 2001, pp. 362–369.
  * [43] D. Stephenson, “Correlation of spatial climate/weather maps and the advantages of using the mahalanobis metric in predictions,” _Tellus A_ , vol. 49, no. 5, pp. 513–527, 1997.
  * [44] H. Shao, W. H. Lam, A. Sumalee, A. Chen, and M. L. Hazelton, “Estimation of mean and covariance of peak hour origin–destination demands from day-to-day traffic counts,” _Transportation Research Part B: Methodological_ , vol. 68, pp. 52–75, 2014.
  * [45] M. N. Ismail, A. Aborujilah, S. Musa, and A. Shahzad, “Detecting flooding based dos attack in cloud computing environment using covariance matrix approach,” in _Proceedings of the 7th international conference on ubiquitous information management and communication_ , 2013, pp. 1–6.
  * [46] I. T. Joliffe and B. Morgan, “Principal component analysis and exploratory factor analysis,” _Statistical Methods in Medical Research_ , vol. 1, no. 1, pp. 69–95, 1992.
  * [47] E. Elhaik, “Principal component analyses (PCA)-based findings in population genetic studies are highly biased and must be reevaluated,” _Scientific Reports_ , vol. 12, no. 1, pp. 1–35, 2022.
  * [48] R. Levie, W. Huang, L. Bucci, M. M. Bronstein, and G. Kutyniok, “Transferability of spectral graph convolutional neural networks.” _J. Mach. Learn. Res._ , vol. 22, pp. 272–1, 2021.
  * [49] Q. Zhu, C. Yang, Y. Xu, H. Wang, C. Zhang, and J. Han, “Transfer learning of graph neural networks with ego-graph information maximization,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 1766–1779, 2021.
  * [50] S. Maskey, R. Levie, and G. Kutyniok, “Transferability of graph neural networks: an extended graphon approach,” _Applied and Computational Harmonic Analysis_ , vol. 63, pp. 48–83, 2023.
  * [51] Y. Zeighami and A. C. Evans, “Association vs. prediction: The impact of cortical surface smoothing and parcellation on brain age,” _Frontiers in big Data_ , vol. 4, p. 637724, 2021.
  * [52] A. I. Luppi and E. A. Stamatakis, “Combining network topology and information theory to construct representative brain networks,” _Network Neuroscience_ , vol. 5, no. 1, pp. 96–124, 2021.
  * [53] J. Royer, R. Rodríguez-Cruces, S. Tavakol, S. Larivière, P. Herholz, Q. Li, R. Vos de Wael, C. Paquola, O. Benkarim, B.-y. Park _et al._ , “An open MRI dataset for multiscale neuroscience,” _Scientific Data_ , vol. 9, no. 1, pp. 1–12, 2022.
  * [54] B. S. Khundrakpam, J. Tohka, A. C. Evans, B. D. C. Group _et al._ , “Prediction of brain maturity based on cortical thickness at different spatial resolutions,” _Neuroimage_ , vol. 111, pp. 350–359, 2015.
  * [55] C. López-Otín _et al._ , “The hallmarks of aging,” _Cell_ , vol. 153, no. 6, pp. 1194–1217, 2013.
  * [56] G. B. Frisoni, N. C. Fox, C. R. Jack, P. Scheltens, and P. M. Thompson, “The clinical use of structural MRI in Alzheimer disease,” _Nature Reviews Neurology_ , vol. 6, no. 2, pp. 67–77, 2010.
  * [57] J. J. Chen, “Functional MRI of brain physiology in aging and neurodegenerative diseases,” _Neuroimage_ , vol. 187, pp. 209–225, 2019\.
  * [58] J. Lee, B. J. Burkett, H.-K. Min, M. L. Senjem, E. S. Lundt, H. Botha, J. Graff-Radford, L. R. Barnard, J. L. Gunter, C. G. Schwarz _et al._ , “Deep learning-based brain age prediction in normal aging and dementia,” _Nature Aging_ , vol. 2, no. 5, pp. 412–424, 2022.
  * [59] W. Yin, L. Li, and F.-X. Wu, “Deep learning for brain disorder diagnosis based on fMRI images,” _Neurocomputing_ , vol. 469, pp. 332–345, 2022.
  * [60] P. R. Millar, B. A. Gordon, P. H. Luckett, T. Benzinger, C. Cruchaga, A. M. Fagan, J. Hassenstab, R. J. Perrin, S. E. Schindler, R. F. Allegri _et al._ , “Multimodal brain age estimates relate to Alzheimer disease biomarkers and cognition in early stages: A cross-sectional observational study,” _bioRxiv_ , 2022.
  * [61] J. Dafflon, W. H. Pinaya, F. Turkheimer, J. H. Cole, R. Leech, M. A. Harris, S. R. Cox, H. C. Whalley, A. M. McIntosh, and P. J. Hellyer, “An automated machine learning approach to predict brain age from cortical anatomical measures,” _Human Brain Mapping_ , vol. 41, no. 13, pp. 3555–3566, 2020\.
  * [62] L. Baecker, J. Dafflon, P. F. Da Costa, R. Garcia-Dias, S. Vieira, C. Scarpazza, V. D. Calhoun, J. R. Sato, A. Mechelli, and W. H. Pinaya, “Brain age prediction: A comparison between machine learning models using region-and voxel-based morphometric data,” _Human Brain Mapping_ , vol. 42, no. 8, pp. 2332–2346, 2021.
  * [63] H. M. Aycheh, J.-K. Seong, J.-H. Shin, D. L. Na, B. Kang, S. W. Seo, and K.-A. Sohn, “Biological brain age prediction using cortical thickness data: A large scale cohort study,” _Frontiers in Aging Neuroscience_ , vol. 10, p. 252, 2018.
  * [64] P. J. LaMontagne, T. L. Benzinger, J. C. Morris, S. Keefe, R. Hornbeck, C. Xiong, E. Grant, J. Hassenstab, K. Moulder, A. G. Vlassenko _et al._ , “OASIS-3: longitudinal neuroimaging, clinical, and cognitive dataset for normal aging and Alzheimer disease,” _MedRxiv_ , 2019.
  * [65] S. M. McGinnis, M. Brickhouse, B. Pascual, and B. C. Dickerson, “Age-related changes in the thickness of cortical zones in humans,” _Brain Topography_ , vol. 24, no. 3, pp. 279–291, 2011.
  * [66] M. Lehmann, S. J. Crutch, G. R. Ridgway, B. H. Ridha, J. Barnes, E. K. Warrington, M. N. Rossor, and N. C. Fox, “Cortical thickness and voxel-based morphometry in posterior cortical atrophy and typical Alzheimer’s disease,” _Neurobiology of aging_ , vol. 32, no. 8, pp. 1466–1476, 2011\.
  * [67] M. R. Sabuncu, R. S. Desikan, J. Sepulcre, B. T. T. Yeo, H. Liu, N. J. Schmansky, M. Reuter, M. W. Weiner, R. L. Buckner, R. A. Sperling _et al._ , “The dynamics of cortical and hippocampal atrophy in Alzheimer disease,” _Archives of neurology_ , vol. 68, no. 8, pp. 1040–1048, 2011.
  * [68] S. G. Popescu, B. Glocker, D. J. Sharp, and J. H. Cole, “Local brain-age: a u-net model,” _Frontiers in Aging Neuroscience_ , vol. 13, p. 761954, 2021\.
  * [69] X. Feng, Z. C. Lipton, J. Yang, S. A. Small, F. A. Provenzano, Alzheimer’s Disease Neuroimaging Initiative, F. L. D. N. Initiative _et al._ , “Estimating brain age based on a uniform healthy population with deep learning and structural magnetic resonance imaging,” _Neurobiology of aging_ , vol. 91, pp. 15–25, 2020.
  * [70] X. Feng, J. Yang, Z. C. Lipton, S. A. Small, F. A. Provenzano, Alzheimer’s Disease Neuroimaging Initiative _et al._ , “Deep learning on MRI affirms the prominence of the hippocampal formation in Alzheimer’s disease classification,” _bioRxiv_ , p. 456277, 2018.
  * [71] A. Essemlali, E. St-Onge, M. Descoteaux, and P.-M. Jodoin, “Understanding Alzheimer disease’s structural connectivity through explainable AI,” in _Medical Imaging with Deep Learning_. PMLR, 2020, pp. 217–229.
  * [72] H. T. Karim, M. Ly, G. Yu, R. Krafty, D. L. Tudorascu, H. J. Aizenstein, and C. Andreescu, “Aging faster: Worry and rumination in late life are associated with greater brain age,” _Neurobiology of Aging_ , vol. 101, pp. 13–21, 2021.
  * [73] F. Liem, G. Varoquaux, J. Kynast, F. Beyer, S. K. Masouleh, J. M. Huntenburg, L. Lampe, M. Rahim, A. Abraham, R. C. Craddock _et al._ , “Predicting brain-age from multimodal imaging data captures cognitive impairment,” _Neuroimage_ , vol. 148, pp. 179–188, 2017.
  * [74] J. Shlens, “A tutorial on principal component analysis,” _arXiv preprint arXiv:1404.1100_ , 2014.
  * [75] A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” _Proceedings of the IEEE_ , vol. 106, no. 5, pp. 808–828, 2018.