# Transferability of coVariance Neural Networks and Application to
Interpretable Brain Age Prediction using Anatomical Features

Saurabh Sihag  University of Pennsylvania, Philadelphia, PA. Gonzalo Mateos
University of Rochester, Rochester, NY.  Corey T. McMillan  University of
Pennsylvania, Philadelphia, PA. Alejandro Ribeiro  University of Pennsylvania,
Philadelphia, PA.

###### Abstract

Graph convolutional networks (GCN) leverage topology-driven graph
convolutional operations to combine information across the graph for inference
tasks. In our recent work, we have studied GCNs with covariance matrices as
graphs in the form of coVariance neural networks (VNNs) that draw similarities
with traditional PCA-driven data analysis approaches while offering
significant advantages over them. In this paper, we first focus on
theoretically characterizing the transferability of VNNs. The notion of
transferability is motivated from the intuitive expectation that learning
models could generalize to “compatible” datasets (possibly of different
dimensionalities) with minimal effort. VNNs inherit the scale-free data
processing architecture from GCNs and here, we show that VNNs exhibit
transferability of performance (without re-training) over datasets whose
covariance matrices converge to a limit object. Multi-scale neuroimaging
datasets enable the study of the brain at multiple scales and hence, provide
an ideal scenario to validate the theoretical results on the transferability
of VNNs. To gauge the advantages offered by VNNs in neuroimaging data
analysis, we focus on the task of “brain age” prediction using cortical
thickness features. In clinical neuroscience, there has been an increased
interest in machine learning algorithms derived from MRI cortical thickness
features which provide estimates of “brain age” that deviate from
chronological age. Importantly, discordance between brain age and
chronological age (“brain age gap”) can reflect increased vulnerability or
resilience toward neurological disease or cognitive impairments. We leverage
the architecture of VNNs to extend beyond the coarse metric of brain age gap
in Alzheimer’s disease (AD) and make two important observations: (i) VNNs can
assign anatomical interpretability to elevated brain age gap in AD by
identifying contributing brain regions, and (ii) the interpretability offered
by VNNs is contingent on their ability to exploit specific principal
components of the anatomical covariance matrix. We further leverage the
transferability of VNNs to cross validate the aforementioned observations
across datasets of different dimensionalities.

Under review. Portions of this manuscript have appeared in [1] and [2]. See
Section 7 for information regarding data and code availability.

##  1 Introduction

In various modern applications, the number of features (denoted by $m$) in a
dataset is a fundamental component of data acquisition that is typically a
characteristic of the desired application and logistics involved [3, 4]. Most
machine learning algorithms and statistical inference approaches are designed
over a pre-defined feature set and hence, their computational and sample
complexities inherently depend on the dimensionality $m$ [5, 6]. In this
paper, we study a deep learning framework called coVariance neural networks
(VNN) [1] that is based on graph neural networks operating on sample
covariance matrix from a given dataset and is scale-free, i.e., the number of
learnable parameters in VNN is independent of the dimensionality $m$ of the
dataset (Fig. 1). The scale-free aspect of VNNs makes it feasible for them to
be _transferable_ between datasets of different dimensionalities without any
changes to their architecture, i.e., VNNs trained on a dataset with
dimensionality $m=m_{1}$ can process another dataset with dimensionality
$m=m_{2}$ with the same set of learned parameters.

While a larger number of features $m$ in a dataset may imply higher resolution
or quality of data collected, too many features can lead to challenges related
to storage, computational complexity, and interpretability of statistical
models for effective data analysis. On the other hand, a dataset with too few
features may be devoid of enough relevant information for accuracy and
inference. Nevertheless, one can intuitively expect some correspondence
between a dataset consisting of $m=m_{1}$ features and another dataset
consisting of $m=m_{2}$ features if both sets of features describe a similar
phenomenon at different scales or resolutions. Conventional data analysis
approaches (e.g., principal component analysis) and machine learning
algorithms are unable to exploit or accommodate this aspect of similarity
between datasets when the number of learnable parameters for inference is
determined by the dimensionality $m$ of the dataset. Hence, such statistical
models need to be re-designed from scratch if the data dimensionality changes.
In this context, we provide the theoretical conditions on the covariance
matrices under which the performance achieved by a VNN model for an inference
task over one dataset can be transferred to that over another dataset with
different dimensionality without re-training or changes to the VNN model (Fig.
3). Besides the methodological gains over traditional data analysis
approaches, VNN frameworks also offer advantages in managing computational
complexity. Indeed, under appropriate conditions, a model trained on a lower
feature count dataset can be directly applied for inference from a higher
feature count dataset.

Neuroimaging is an example of a modern application in which the number of
features is highly variable across datasets, but different datasets contain
similar information [7, 8]. Specifically, MRI data can be represented in many
scales ranging from single voxels (typically $\sim 1\text{ mm}^{3}$) to
regions-of-interest (ROIs) derived from multi-scale brain atlases that range
from dozens to thousands of parcellations (e.g., from 100 to 1000 number of
parcellations in a multi-scale brain atlas [9, 10]). Multi-scale brain atlases
provide mappings that divide the brain cortex into different number of
parcellations and therefore, associated neuroimaging datasets describe similar
information over the brain cortex at different scales. Hence, we empirically
validate the theoretical results for transferability of VNNs on a set of
neuroimaging datasets curated according to different scales of a commonly used
brain atlas.

The simplicity of VNN framework allows us to analyze the contributions of each
feature to the final statistical outcome of a VNN model. This aspect is
significant for applications that use neuroimaging datasets, as the features
therein are associated with distinct brain regions. Hence, VNNs make it
feasible to assign anatomical or regional _interpretability_ to the learning
outcomes when trained on datasets consisting of anatomical features. Motivated
by these observations, we pursue the task of “brain age” prediction from
cortical thickness features derived from magnetic resonance imaging (MRI) as
an application of VNNs (Fig. 4). Critically, individuals may experience age-
related effects at different rates, captured by so-called “biological aging”.
Hence, accelerated aging (e.g., when biological age is elevated as compared to
chronological age) may predict age-related vulnerabilities like risk for
cognitive decline or neurological conditions like Alzheimer’s disease and
related dementias (ADRD) [11]. In this domain, the metric of interest is
_brain age gap_ , i.e., the difference between the biological age and the
chronological age. We use the notation $\Delta$-Age to refer to the brain age
gap. Since brain age has no ground truth, $\Delta$-Age is essentially a
qualitative metric that is expected to be elevated in individuals with
underlying neurodegenerative condition as compared to the healthy population.
Numerous existing studies based on a large spectrum of machine learning
approaches report elevated $\Delta$-Age in neurodegenerative conditions,
including Alzheimer’s disease [12] and schizophrenia [13]. However, several
criticisms for brain age evaluation approaches using machine learning have
also been identified. Major criticisms include the coarseness of $\Delta$-Age
that results in lack of specificity of brain regions contributing to elevated
brain age; and unexplained reliance on the prediction accuracy for
chronological age in the design of these machine learning models [14]. The
architecture of VNN enables us to propose a principled framework for brain age
prediction that accommodates interpretability by isolating the contributing
brain regions to elevated $\Delta$-Age in neurodegeneration. These
contributing brain regions are identified by analyzing the group differences
between the group with the neurodegenerative condition and healthy controls
with respect to the contributions of different features to the corresponding
predictions made by VNNs (Fig. 4). A layman overview of the advantage offered
by VNNs in terms of adding interpretability to brain age prediction is
included in Appendix A.

Together, the interpretability and transferability aspects of VNNs provide
novel insights into the role of training the model to predict chronological
age in the brain age prediction framework. A significant portion of existing
literature that studies brain age using deep learning models considers the
ability of their models to accurately predict chronological age (time since
birth) for healthy controls [15, 16, 17] as a relevant metric for assessing
quality of their methodological approach. Simultaneously, deep learning models
that have a relatively moderate fit on the chronological age of healthy
controls can also provide better insights into brain age than the ones with a
tighter fit [16]. Thus, there is a lack of conceptual clarity in the role of
the quality of fit against chronological age of healthy controls in predicting
a meaningful brain age [18], and this issue is unlikely to be addressed by any
statistical approach for evaluating brain age if it lacks transparency
(Appendix A).

It is reasonable to postulate that an accurate or a near-perfect prediction of
chronological age in healthy controls does not necessarily equip the deep
learning models to accommodate neurodegeneration-driven changes in the brain.
The study of interpretability and transferability aspects of VNNs demonstrates
that the contributing brain regions behind elevated $\Delta$-Age in
neurodegeneration may be qualitatively recovered after transferring VNNs from
one dataset to another. Moreover, we also find that the ability of VNNs to
exploit specific eigenvectors of the anatomical covariance matrix is the
underlying factor behind the anatomical interpretability offered by VNNs in
the context of $\Delta$-Age. Hence, VNNs can facilitate a principled
decoupling of the objective of inferring brain age from the aim of achieving a
near-perfect performance on the task of predicting chronological age for
healthy controls.

Contributions: Our contributions in this paper are summarized as follows:

  * -

Transferability of VNNs: We theoretically characterize the transferability of
VNNs between datasets of different dimensionalities. For a dataset with
$m_{1}$ features and covariance matrix ${\mathbf{C}}_{m_{1}}$ and another
dataset with $m_{2}$ features and covariance matrix ${\mathbf{C}}_{m_{2}}$, we
demonstrate that the outputs of a VNN when initialized on
${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$ are close in some sense
under appropriate conditions on covariance matrices ${\mathbf{C}}_{m_{1}}$ and
${\mathbf{C}}_{m_{2}}$ (see Theorem 3). The theoretical results on
transferability of VNNs were validated on a regression task based on a set of
cortical thickness datasets curated according to different scales of a
commonly used multi-scale brain atlas (Fig. 5 and Tables 1 and 2).

  * -

Brain age prediction using VNNs: We deployed VNNs for the task of brain age
prediction and compared the $\Delta$-Age between healthy controls and
individuals with AD diagnosis. The insights gained in this set of experiments
are summarized below:

    * a)

VNNs provide anatomically interpretable $\Delta$-Age in neurodegeneration:
$\Delta$-Age in individuals with AD diagnosis was elevated as compared to
healthy controls. The simplicity of the VNN architecture allowed us to
characterize the regional contributors to the elevated $\Delta$-Age, thus,
adding anatomical interpretability to $\Delta$-Age (Fig. 5.3.1). On a multi-
scale dataset, the transferability of VNNs also helped cross-validate the
spatial robustness of the observed regional profiles for $\Delta$-Age in
neurodegeneration (Fig. 20 in Appendix K).

    * b)

Anatomical interpretability was correlated with eigenvectors of the anatomical
covariance matrix: Our experiments demonstrated that there was a correlation
between specific eigenvectors of the anatomical covariance matrix and the
features that facilitated anatomical interpretability for $\Delta$-Age (Fig.
8). Thus, $\Delta$-Age was linked to the ability of VNNs to exploit specific
eigenvectors of the anatomical covariance matrix.

    * c)

Clarity in the role of training VNNs to predict chronological age: Our
experiments conclusively showed that training the VNNs to predict
chronological age enabled them to exploit the eigenvectors of the anatomical
covariance matrix associated with elevated $\Delta$-Age in neurodegeneration
(Fig. 9, Fig. 10, and Fig. 19). The observations made in this context
facilitated decoupling of the brain age prediction task from the objective of
achieving a near perfect performance on the prediction of chronological age in
healthy controls.

Next, we provide a literature review pertinent to our contributions in this
paper.

Related Literature. Graph neural networks (GNNs) are a widely popular
adaptation of convolutional neural networks to graph-structured data [19, 20,
21]. Graphs are natural descriptors of complex, spatially-distributed
phenomena and therefore, graph-structured datasets are prevalent in a variety
of application domains [22], including physical infrastructure [23], social
network analysis [24], biology [25], network neuroscience [26] and natural
sciences [27]. Processing of graph-structured data faces various practical
challenges (like generalizability, reproducibility, scalability) [22] and
therefore, a number of variants of GNNs have been proposed to address them. We
refer the reader to recent survey articles in [21] and [28] that categorize
GNN architectures according to diverse criteria, including mathematical
formulations, algorithms, and hardware/software implementations. Convolutional
GNNs [29], graph autoencoders [30], recurrent GNNs, and gated GNNs [31] are
among a few prominently studied and applied categories of GNNs.

The taxonomy pertinent to this paper is that of graph convolutional networks
[29]. GCNs typically rely on an information aggregation procedure (referred to
as graph convolutions) over a graph structure for data processing. Several
implementation strategies for graph convolution operations have been proposed
in the literature, including spectral convolutions [32], Chebyshev polynomials
[33], ordinary polynomials [34], and diffusion based representations [35].
GCNs admit the properties of stability to topological perturbations and
transferability across graphs of different sizes in various settings [36, 37,
38, 39], which makes them a well-motivated data analysis tool for graph-
structured data.

In our recent work in [1], we studied coVariance neural networks (VNN), which
are GCNs with sample covariance matrices as graph and polynomial graph filters
as convolution operation. Covariance matrices and principal component analysis
(PCA) form the two cornerstones of non-parametric analyses in real world
applications that have spatially distributed, multi-variate data acquisition
protocols, including neuroimaging [40], computer vision [41, 42], weather
modeling [43], traffic flow analysis [44], and cloud computing [45]. Our
results in [1] established the following significant observations i) there
exist similarities between the spectral analysis of graph convolution on
covariance matrix and the standard PCA transformation; and ii) VNNs are robust
to the number of samples used to estimate the sample covariance matrix, thus,
overcoming a potential source of instability and irreproducibility of PCA
based statistical inference [46, 47].

The transferability of GNNs from training graphs to some compatible family of
test graphs has been previously studied from different theoretical
perspectives [48, 39, 49, 50]. The notion of transferability of GNNs broadly
encapsulates the intuition that GNNs may be able to retain their performance
for some inference task when applied over test graphs (irrespective of the
size) that describe the same phenomenon as the training graphs. In this
context, the study in [48] considers transferability of GNNs over graphs that
represent the discretization of underlying topological spaces. Several studies
also consider GNN transferability over graphs that belong to a converging
sequence that approaches a limiting object in the asymptote of a large number
of nodes [39, 50]. In [49], the similarity between the ego graph distributions
(derived from graph topology) formed the workhorse for assessing
transferability of GNNs. Transferability of GNNs also provides advantages in
terms of computational complexity, which for GNNs scales as ${\cal O}(m^{2})$
for dense graphs with $m$ nodes. In this paper, we extend the notion of
transferability to VNNs and establish the transference over covariance
matrices of different sizes that converge in some sense. In this context,
transferability is not feasible for traditional PCA-driven statistical models,
as the principal components are restricted within the feature space of the
original dataset and need to be re-evaluated if the number of features change.
Specifically, PCA does not provide any notion of similarity between the
principal components extracted from a covariance matrix of size $m_{1}\times
m_{1}$ and that from another covariance matrix of size $m_{2}\times m_{2}$.
Thus, transferability of VNNs is broadly relevant to the domain of
multivariate statistics.

There has been a growing interest in multi-scale datasets in neuroscience [51,
52, 53, 8, 54]. These datasets rely on brain atlases or templates that allow a
multi-scale parcellation of the brain surface (for instance, Schaefer’s atlas
[9] and Lausanne atlas [10]). A multi-scale brain atlas partitions the brain
cortex into a variable number of regions at different scales. However,
statistically sound approaches that optimally leverage the redundancy of
information in datasets consisting of features at multiple scales are
currently lacking. In this paper, we leverage the cortical thickness datasets
curated according to multi-scale brain atlases to validate the transferability
of VNNs. In scenarios where the theoretical guarantees for VNN transferability
do not hold between datasets curated according to different brain atlases,
appropriate mappings that account for the differences between different brain
atlases [7, 8] may be necessary to scale VNN outputs for good quantitative
performance on the chronological age prediction task when VNNs are transferred
from one dataset to another. However, these mappings are unlikely to possess
any information regarding neurodegeneration or brain age and hence, are not
studied here.

The human aging process is characterized by progressive anatomical and
functional changes in the brain [55]. The $\Delta$-Age for a pathology can be
seen as a scalar representation of longitudinal, pathology-driven atypical
changes in the brain [12]. Data from neuroimaging modalities, including
structural magnetic resonance imaging (MRI), functional MRI, and positron
emission tomography, capture the changes of the brain due to neurodegeneration
and healthy aging [56, 57]. Thus, inferring brain age from different
neuroimaging modalities has been an active area of research [58, 59, 60, 61,
62, 63]. In this paper, we leverage the datasets of cortical thickness
measures derived from structural MRI images in OASIS-3 dataset to study brain
age [64]. Cortical thickness evolves with normal aging [65] and is impacted
due to neurodegeneration in AD [66, 67]. Thus, the age-related and disease
severity related variations also appear in anatomical covariance matrices
evaluated from the correlation among the cortical thickness measures across a
population [40].

Due to inherent interpretability offered by VNNs to their statistical
outcomes, VNNs provide an explainable regional profile to the elevated
$\Delta$-Age in neurodegeneration. Limited focus has been on comparable
studies in this regard that associate brain age gaps with regional profiles
[15, 68]. The study in [15] adopts a convolutional neural network approach to
infer brain age from MRI images directly and assigns importance to brain
regions in evaluating the brain age. In principle, the interpretability
offered by VNNs in the context of brain age is similar, as we infer a regional
profile for $\Delta$-Age by isolating the brain regions that are contributors
to the elevated $\Delta$-Age in neurodegeneration. In addition, the regional
profile identified by VNNs is correlated with specific eigenvectors or the
principal components of the anatomical covariance matrix. Hence, the
$\Delta$-Age inferred by our framework is driven by the ability of a VNN to
manipulate the input data according to certain principal components of the
anatomical covariance matrix. Also, VNNs are significantly less complex deep
learning models as compared to those studied in [15]. Our results demonstrate
that the VNNs trained with less than $2000$ learnable parameters exhibit
(spatially robust) regional interpretability in the context of brain age in
AD. In general, the regional expressivity offered by VNNs is in stark contrast
to a multitude of existing relevant studies that rely on less transparent
statistical approaches and further use post-hoc analyses (such as ablation
analysis [69, 70, 71] or exploring correlations with region-specific markers
[58] and psychiatric symptoms [72, 73]) to assign interpretability to a
scalar, elevated $\Delta$-Age effect.

##  2 coVariance Neural Networks

VNNs operate on covariance matrices and have similar architecture as GNNs
111GCNs and GNNs are used interchangeably in the rest of the paper.. We start
by providing preliminary definitions pertaining to the architecture and
discuss the theoretical properties associated with VNNs later.

Consider an $m-$dimensional random vector
${\mathbf{x}}\in{\mathbb{R}}^{m\times 1}$ whose ensemble covariance matrix is
defined as

$\displaystyle{\mathbf{C}}\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}{\mathbb{E}}[({\mathbf{x}}-{\mathbb{E}}[{\mathbf{x}}])({\mathbf{x}}-{\mathbb{E}}[{\mathbf{x}}])^{\sf
T}]\;,$ (1)

where $\cdot^{\sf T}$ is the transpose operator and ${\mathbb{E}}[\cdot]$ is
the expectation with respect to the probability distribution of
${\mathbf{x}}$. In practice, we usually have access to a dataset that provides
us with the statistical information about ${\mathbf{x}}$. Therefore, we also
consider a dataset consisting of $n$ random, independent and identically
distributed (i.i.d) samples of ${\mathbf{x}}$, given by
${\mathbf{x}}_{i}\in{\mathbb{R}}^{m\times 1},\forall i\in\\{1,\dots,n\\}$,
where the dataset can be represented in the matrix form as
${\mathbf{X}}_{n}=[{\mathbf{x}}_{1},\dots,{\mathbf{x}}_{n}]$. Using
${\mathbf{X}}_{n}$, we estimate the ensemble covariance matrix, conventionally
referred to as the sample covariance matrix as follows

$\displaystyle\hat{\mathbf{C}}\triangleq\frac{1}{n-1}\sum\limits_{i=1}^{n}({\mathbf{x}}_{i}-\bar{\mathbf{x}})({\mathbf{x}}_{i}-\bar{\mathbf{x}})^{\sf
T}\;,$ (2)

where $\bar{\mathbf{x}}$ is the sample mean of samples in ${\mathbf{X}}_{n}$.
Next, we discuss the motivation behind studying VNNs separately from GNNs.

###  2.1 Motivation

Covariance matrices are ubiquitous in real world applications that have
spatially distributed, multi-variate data acquisition protocols [40, 43, 44,
45]. The eigenvectors of covariance matrices are termed as principal
components of the dataset and constitute the well-known PCA transformation
[74]. Our discussion in this section shows that the spectral domain
representation of the graph convolution operator instantiated on the
covariance matrix as graph yields the PCA transformation. This observation
suggests that learning with a GNN with covariance matrix as a graph is
achieved (at least in part) by manipulation of the data according to the
eigenbasis of the covariance matrix. Therefore, this paper provides a
conceptual contribution towards the study of GNNs instantiated on covariance
matrices in the form of VNNs. We also provide theoretical guarantees that
result in significant advantages over models that perform statistical
inference using PCA.

The covariance matrix ${\mathbf{C}}$ can be viewed as the adjacency matrix of
a graph representing the stochastic structure of the vector ${\mathbf{x}}$,
where the $m$ dimensions of ${\mathbf{x}}$ can be thought of as the nodes of
an $m$-node, undirected graph and its edges represent the pairwise covariance
between elements in ${\mathbf{x}}$. Furthermore, the eigenvalues of
${\mathbf{C}}$ encode the variability of the dataset along different
directions in an orthogonal space determined by the associated eigenvectors or
principal components.

In graph signal processing, a vector defined on the nodes of the graph is
viewed as the graph signal and the projection of a graph signal on the
eigenbasis of the graph yields the graph Fourier transform [75]. The graph
Fourier transform provides a systematic mathematical tool to analyze
convolutional filters over graphs [29, 76]. Interestingly, the classical
Fourier transform and graph Fourier transform converge over a discrete,
periodic time series represented on a directed, cyclic graph [77]. Similarly
to the graph Fourier transform, we can define the coVariance Fourier transform
as the projection of a random instance ${\mathbf{x}}$222For ease of notation,
we will subsequently use the notation ${\mathbf{x}}$ to refer to a random
instance of the random vector whose covariance matrix is ${\mathbf{C}}$. on
the eigenvectors of the covariance matrix ${\mathbf{C}}$ [1, Definition 1].
The definition of coVariance Fourier transform from [1] is stated next. For
this purpose, we leverage the eigendecomposition of ${\mathbf{C}}$ given by

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

where ${\mathbf{V}}=[{\mathbf{v}}_{1},\dots,{\mathbf{v}}_{m}]$ is a matrix of
size $m\times m$ with its columns as the eigenvectors and
${\mathbf{\Lambda}}={\sf diag}(\lambda_{1},\dots,\lambda_{m})$ is a diagonal
matrix with its diagonal elements representing the eigenvalues of
${\mathbf{C}}$.

######  Definition 1 (coVariance Fourier Transform).

The coVariance Fourier transform (VFT) of a random sample ${\mathbf{x}}$ is
defined as its projection on the eigenspace of ${\mathbf{C}}$ and is given by

$\displaystyle\tilde{\mathbf{x}}\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}{\mathbf{V}}^{\sf
T}{\mathbf{x}}\;.$ (4)

The $i$-th entry of $\tilde{\mathbf{x}}$, i.e., $[\tilde{\mathbf{x}}]_{i}$
represents the projection of ${\mathbf{x}}$ on eigenvector ${\mathbf{v}}_{i}$
and hence, it is associated with the eigenvalue $\lambda_{i}$. Thus, the
similarity between PCA transformation and VFT in (4) implies that eigenvalue
$\lambda_{i}$ encodes the variability of the dataset ${\mathbf{X}}_{n}$ in the
direction of the principal component ${\mathbf{v}}_{i}$. In this context, the
eigenvalues of the covariance matrix are the mathematical equivalent of the
notion of graph frequencies in graph signal processing [75].

GNNs with convolutional filters that rely on a _linear shift-and-sum_
operation fundamentally exhibit the stability to changes in graph topology
[38]. Since the sample covariance matrix $\hat{\mathbf{C}}$ is likely to be
perturbed with respect to ${\mathbf{C}}$ [78], stability is desirable to
mitigate the impact of number of samples on statistical inference. Motivated
by this observation, we define the notion of coVariance filters (VF) that are
polynomials in the covariance matrix and characterize the convolution
operation in VNNs.

######  Definition 2 (coVariance Filters).

Given a set of real valued, scalar parameters $\\{h_{k}\\}_{k=0}^{K}$, the
coVariance filter on a covariance matrix ${\mathbf{C}}$ is defined as

$\displaystyle{\mathbf{H}}({\mathbf{C}})\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\sum\limits_{k=0}^{K}h_{k}{\mathbf{C}}^{k}\;.$
(5)

Furthermore, the output of the covariance filter ${\mathbf{H}}({\mathbf{C}})$
for an input ${\mathbf{x}}$ is given by

$\displaystyle{\mathbf{z}}={\mathbf{H}}({\mathbf{C}}){\mathbf{x}}\;.$ (6)

The application of coVariance filter ${\mathbf{H}}({\mathbf{C}})$ on an input
${\mathbf{x}}$ translates to combining information across different sized
neighborhoods. To elucidate this observation, consider a coVariance filter
with $K=1$ and $h_{0}=0$. In this scenario, the $i$-th element of
${\mathbf{z}}$ is evaluated as

$\displaystyle[{\mathbf{z}}]_{i}=h_{1}\sum\limits_{j=1}^{m}[{\mathbf{C}}]_{ij}[{\mathbf{x}}]_{j}\;.$
(7)

Thus, $[{\mathbf{z}}]_{i}$ represents the linear combination of all elements
according to the $i$-th row of ${\mathbf{C}}$. This observation implies that
the neighborhood of $i$-th dimension in ${\mathbf{x}}$ (derived from the graph
representation ${\mathbf{C}}$) determines the outcome $[{\mathbf{z}}]_{i}$ of
the convolution operation in (7). For $K>1$, the convolution operation
combines information across multi-hop neighborhoods (up to $K$-hop) according
to the weights $h_{k}$.

The spectral analysis of the covariance filtering operation in Definition 2
via VFT of the filter output ${\mathbf{z}}$ yields the frequency response of
the covariance filter and reveals the similarities between covariance
filtering and PCA. After taking the VFT of ${\mathbf{z}}$, we have

$\displaystyle\tilde{\mathbf{z}}$ $\displaystyle={\mathbf{V}}^{\sf
T}{\mathbf{H}}({\mathbf{C}}){\mathbf{x}}\;,$ (8)
$\displaystyle=\sum\limits_{k=0}^{K}h_{k}{\mathbf{\Lambda}}^{k}{\mathbf{V}}^{\sf
T}{\mathbf{x}}=\sum\limits_{k=0}^{K}h_{k}{\mathbf{\Lambda}}^{k}\tilde{\mathbf{x}}\;,$
(9)

where $\tilde{\mathbf{x}}={\mathbf{V}}^{\sf T}{\mathbf{x}}$ is the covariance
Fourier transform of ${\mathbf{x}}$ and (9) follows from (8) from the
orthonormality of eigenvectors of ${\mathbf{C}}$. The frequency response of
the coVariance filter depends on the filter taps $\\{h_{k}\\}$ and the
eigenvalues of ${\mathbf{C}}$ and is given by

$\displaystyle h(\lambda)=\sum\limits_{k=0}^{K}h_{k}\lambda^{k}\;.$ (10)

Furthermore, since $\tilde{\mathbf{x}}$ is a projection of ${\mathbf{x}}$ on
the eigenvector space ${\mathbf{V}}$ and
$[\tilde{\mathbf{z}}]_{i}=h(\lambda_{i})[{\mathbf{V}}^{\sf
T}{\mathbf{x}}]_{i}$, the $i$-th element of $\tilde{\mathbf{z}}$ yields
similarities with the standard PCA transformation. This observation is
formalized in Lemma 1.

######  Lemma 1 (Spectrum of coVariance Filter and PCA).

Given a covariance matrix ${\mathbf{C}}$ with eigendecomposition in (3), if
the PCA transformation of input ${\mathbf{x}}$ is given by
${\mathbf{y}}={\mathbf{V}}^{\sf T}{\mathbf{x}}$, there exists a filter bank of
coVariance filters $\\{{\mathbf{H}}_{i}({\mathbf{C}}):i\in\\{1,\dots,m\\}\\}$,
such that, the score of the projection of input ${\mathbf{x}}$ on eigenvector
${\mathbf{v}}_{i}$ can be recovered by the application of a coVariance filter
${\mathbf{H}}_{i}({\mathbf{C}})$ as:

$\displaystyle[{\mathbf{y}}]_{i}={\mathbf{v}}_{i}^{\sf
T}{\mathbf{H}}_{i}({\mathbf{C}}){\mathbf{x}}\;,$ (11)

where the frequency response $h_{i}(\lambda)$ of the filter
${\mathbf{H}}_{i}({\mathbf{C}})$ is given by

$h_{i}(\lambda)=\begin{cases}\omega_{i},\quad\text{if}\quad\lambda=\lambda_{i}\;,\\\
0,\quad\text{otherwise}\end{cases}\;.$ (12)

Lemma 1 establishes equivalence between processing data samples with PCA and
processing data samples with a specific polynomial on the covariance matrix.

Our previous work in [1] showed that in contrast to PCA involving eigenvectors
of the sample covariance matrix, information processing with polynomials of
the sample covariance matrix can be stable to finite sample induced
perturbations. Indeed, in practice, we may only have access to the sample
covariance matrix $\hat{\mathbf{C}}$ which is an estimate of ${\mathbf{C}}$.
Since $\hat{\mathbf{C}}$ is a consistent estimator of ${\mathbf{C}}$, the
eigenvalues and eigenvectors of $\hat{\mathbf{C}}$ approach those of
${\mathbf{C}}$ in the limit of infinite number of samples, i.e.,
$n\rightarrow\infty$. However, for a finite number of samples $n$, the
eigenvectors and eigenvalues of $\hat{\mathbf{C}}$ are perturbed with respect
to those of ${\mathbf{C}}$ [78]. In principle, statistical inference using PCA
can be prone to instability due to eigenvectors corresponding to eigenvalues
of the ensemble covariance matrix that are close [46] and, thus, lead to
irreproducible statistical conclusions [47]. In this context, we informally
state Theorem 2 from [1].

######  Theorem 1 (Stability of coVariance filter).

Consider a dataset with sample covariance matrix $\hat{\mathbf{C}}$ formed by
$n$ samples and the counterpart ensemble covariance matrix ${\mathbf{C}}$.
Under Assumption 1 in [1], the following holds with high probability:

$\displaystyle\left\lVert{\mathbf{H}}(\hat{\mathbf{C}})-{\mathbf{H}}({\mathbf{C}})\right\rVert={\cal
O}\left(\frac{\nu}{n^{\frac{1}{2}-\varepsilon}}\right)\;,$ (13)

for some $\nu>0$ and $\varepsilon\in(0,1/2)$.

Theorem 1 establishes that information processing using a polynomial of the
covariance matrix offers stability with respect to the perturbations between
the sample covariance matrix $\hat{\mathbf{C}}$ and ${\mathbf{C}}$ [1]. Also,
as a corollary to Theorem 1, we can state that the difference between outputs
of covariance filters instantiated on distinct sample covariance matrices are
bounded. These observations imply that statistical inference based on
covariance filters are characterized by robustness to the effects of finite
sample size and, thus, result in consistent statistical outcomes with high
confidence. No such guarantees are offered by PCA. Next, we discuss the
architecture of VNNs that is based on covariance filters, which results in
VNNs inheriting the stability offered by coVariance filters.

###  2.2 Architecture

We begin with the description of a coVariance perceptron that forms one layer
of the VNN architecture. For this purpose, we leverage the definition of a
pointwise, nonlinear activation function $\sigma(\cdot)$, such that, for
${\mathbf{x}}=[x_{1},\dots,x_{m}]$, we have
$\sigma({\mathbf{x}})=[\sigma(x_{1}),\dots,\sigma(x_{m})]$. Examples of point-
wise, nonlinear activation functions are ${\sf ReLU}$ and $\sf tanh$.

######  Definition 3 (coVariance Perceptron).

For a given pointwise nonlinear activation function $\sigma(\cdot)$, input
${\mathbf{x}}$, a coVariance filter ${\mathbf{H}}({\mathbf{C}})$ and its
corresponding set of filter taps ${\mathcal{H}}$, the coVariance perceptron is
defined as

$\displaystyle\Phi({\mathbf{x}};{\mathbf{C}},{\mathcal{H}})\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\sigma({\mathbf{H}}({\mathbf{C}}){\mathbf{x}})\;.$
(14)

A VNN can be constructed by stacking perceptrons to form multi-layer
information processing architecture. This observation is formalized next.

######  Remark 1 (Multi-layer VNN).

Consider an $L$-layer architecture formed by stacking $L$ coVariance
perceptrons defined in Definition 3. In this scenario, we denote the
coVariance filter in layer $\ell$ of a VNN by
${\mathbf{H}}_{\ell}({\mathbf{C}})$ and its corresponding set of filter taps
are given by ${\mathcal{H}}_{\ell}$. For a given pointwise nonlinear
activation function $\sigma(\cdot)$, the relationship between the input
${\mathbf{x}}_{\ell-1}$ and the output ${\mathbf{x}}_{\ell}$ for the
coVariance perceptron in the $\ell$-th layer is given by

$\displaystyle{\mathbf{x}}_{\ell}=\sigma({\mathbf{H}}_{\ell}({\mathbf{C}}){\mathbf{x}}_{\ell-1})\;\quad\text{for
}\quad\ell\in\\{1,\dots,L\\},$ (15)

where ${\mathbf{x}}_{0}$ is the input ${\mathbf{x}}$. We refer to this
$L$-layer architecture as an $L$-layer VNN.

Figure 1: Basics of VNN architecture. Panel a illustrates that the coVariance
filter ${\mathbf{H}}({\mathbf{C}})$ is a polynomial in ${\mathbf{C}}$ and its
application on input ${\mathbf{x}}$. Panel b shows the construction of a
coVariance perceptron based on coVariance filter ${\mathbf{H}}({\mathbf{C}})$
and pointwise nonlinearity $\sigma$. coVariance perceptron specifies one layer
of VNN. Panel c shows a basic multi-layer VNN architecture formed by stacking
three coVariance perceptrons.

A pictorial illustration of a multi-layer VNN is included in Fig. 1.
Furthermore, similar to other deep learning models, sufficient expressive
power can be facilitated in the VNN architecture by incorporating multiple
input multiple output (MIMO) processing at every layer. Formally, consider a
VNN layer $\ell$ that can process $F_{\ell-1}$ number of $m$-dimensional
inputs and outputs $F_{\ell}$ number of $m$-dimensional outputs via
$F_{\ell-1}\times F_{\ell}$ number of filter banks [38]. In this scenario, the
input is specified as ${\mathbf{X}}_{\sf in}=[{\mathbf{x}}_{\sf
in}[1],\dots,{\mathbf{x}}_{\sf in}[F_{\sf in}]]$, and the output is specified
as ${\mathbf{X}}_{\sf out}=[{\mathbf{x}}_{\sf out}[1],\dots,{\mathbf{x}}_{\sf
out}[F_{\sf out}]]$. The relationship between the $f$-th output
${\mathbf{x}}_{\sf out}[f]$ and the input ${\mathbf{x}}_{\sf in}$ is given by

$\displaystyle{\mathbf{x}}_{\sf out}[f]$
$\displaystyle=\sigma\left(\sum\limits_{g=1}^{F_{\sf
in}}{\mathbf{H}}_{fg}({\mathbf{C}}){\mathbf{x}}_{\sf in}[g]\right)\;,$ (16)

where ${\mathbf{H}}_{fg}({\mathbf{C}})$ is the coVariance filter that
processes ${\mathbf{x}}_{\sf in}[g]$. Without loss of generality, we focus the
subsequent discussion on the scenario when we have
$F_{\ell}=F,\forall\ell\in\\{1,\dots,L\\}$. In this case, the set of all
filter taps is given by ${\cal H}=\\{{\mathcal{H}}_{fg}^{\ell}\\},\forall
f,g\in\\{1,\dots,F\\},\ell\in\\{1,\dots,L\\}$, where
${\mathcal{H}}_{fg}=\\{h_{fg}^{\ell}[k]\\}_{k=0}^{K}$ and $h_{fg}^{\ell}[k]$
is the $k$-th filter tap for filter ${\mathbf{H}}_{fg}({\mathbf{C}})$. Thus,
we can compactly represent a multi-layer VNN architecture capable of MIMO
processing via the notation $\Phi({\mathbf{x}};{\mathbf{C}},{\cal H})$ as the
set of filter taps ${\mathcal{H}}$ captures the full span of its architecture.
We also use the notation $\Phi({\mathbf{x}};{\mathbf{C}},{\cal H})$ to denote
the output of the VNN. For a VNN with $F$ number of $m$-dimensional outputs in
the final layer, the size of the VNN output
$\Phi({\mathbf{x}};{\mathbf{C}},{\cal H})$ is $m\times F$.

The output $\Phi({\mathbf{x}};{\mathbf{C}},{\cal H})$ is succeeded by a
readout function that maps $\Phi({\mathbf{x}};{\mathbf{C}},{\cal H})$ to the
desired output. In this paper, we assume non-adaptive or non-learnable readout
function (e.g., mean, max or min functions) which preserves the property of
permutation invariance for VNN model. Furthermore, a non-adaptive readout
function is essential for the transferability property of VNNs (discussed in
Section 2.3).

It is imperative to study the robustness of VNN outputs to the number of
samples $n$ in order to guarantee reproducibility of VNN statistical outcomes.
Specifically, it is desirable that the change in VNN outputs is controlled or
bounded when the architecture is trained using sample covariance matrices
estimated from $n_{1}$ or $n_{2}$ samples when $n_{1}\neq n_{2}$. In Theorem
2, we informally state the result on the stability of VNNs by analyzing
$\|\Phi({\mathbf{x}};\hat{\mathbf{C}},{\mathcal{H}})-\Phi({\mathbf{x}};{\mathbf{C}},{\mathcal{H}})\|$,
i.e., the difference between the VNN outputs for the sample covariance matrix
$\hat{\mathbf{C}}$ and the ensemble covariance matrix ${\mathbf{C}}$. This
Theorem was also previously established in [1].

######  Theorem 2 (VNN Stability).

Consider an ensemble covariance matrix ${\mathbf{C}}$ and its estimate
$\hat{\mathbf{C}}$ formed from $n$ samples. Given a bank of coVariance filters
with filter taps
${\mathcal{H}}=\\{{\mathcal{H}}_{fg}^{\ell}:f\in\\{1,\dots,F\\},\ell\in\\{1,\dots,L\\}\\}$,
the coVariance filters are stable and satisfy

$\displaystyle\|{\mathbf{H}}_{fg}^{\ell}(\hat{\mathbf{C}})-{\mathbf{H}}_{fg}^{\ell}({\mathbf{C}})\|\leq\alpha_{n}\;,$
(17)

for some $\alpha_{n}>0$ with high probability (Theorem 1). Also, for a
pointwise nonlinearity function $\sigma(\cdot)$, such that,
$|\sigma(a)-\sigma(b)|\leq|a-b|$, we have

$\displaystyle\|\Phi({\mathbf{x}};\hat{\mathbf{C}},{\mathcal{H}})-\Phi({\mathbf{x}};{\mathbf{C}},{\mathcal{H}})\|\leq
LF^{L}\alpha_{n}\;.$ (18)

###### Proof.

See Appendix B in the Supplementary Material. ∎

The parameter $\alpha_{n}$ in (17) represents the finite sample effect on the
perturbations in $\hat{\mathbf{C}}$ with respect to ${\mathbf{C}}$ and is
borrowed from Theorem 1. By leveraging the perturbation theory of covariance
matrices to analyze the stability of coVariance filters, we also show in the
proof of Theorem 2 that $\alpha_{n}$ scales as $1/n^{\frac{1}{2}-\varepsilon}$
for some $\varepsilon\in(0,1/2)$ with respect to the number of samples $n$. We
note that the bound in (18) becomes looser with increase in $F$ or $L$ which
is consistent with the parallel result for GNNs [38]. However, without the
analysis of the lower bounds on
$\|\Phi({\mathbf{x}};\hat{\mathbf{C}},{\mathcal{H}})-\Phi({\mathbf{x}};{\mathbf{C}},{\mathcal{H}})\|$,
we cannot claim that the robustness of VNNs indeed worsens with increase in
$F$ or $L$. Moreover, we remark that VNNs sacrifice discriminability between
eigenvectors associated with close eigenvalues to achieve stability [1]. As a
corollary, we also state that Theorem 2 can readily be extended to
characterize the difference between VNN outputs corresponding to sample
covariance matrices estimated from $n_{1}$ and $n_{2}$ samples via (18) and
application of triangle inequality.

###  2.3 Transferability of VNNs

The notion of transferability of VNNs is made feasible by the properties of
coVariance filters (Definition 2) that can be instantiated on covariance
matrices of any dimension and therefore, the VNNs can readily be transferred
to process a dataset of a different dimensionality. From the perspective of
implementation, transferability of VNNs to a dataset of different number of
features amounts to replacing the covariance matrix ${\mathbf{C}}$ in a VNN
model $\Phi(\cdot;{\mathbf{C}},{\mathcal{H}})$ with a covariance matrix of
another size, while keeping the parameters ${\mathcal{H}}$ fixed. Since we
consider covariance matrices of different dimensionalities, we denote a
covariance matrix ${\mathbf{C}}$ of size $m\times m$ by ${\mathbf{C}}_{m}$.
Informally, we can state our objective for assessing transferability as
follows.

Informal Problem Statement for VNN Transferability. Given a data point
${\mathbf{x}}_{m_{1}}$ from a dataset with $m_{1}$ features and associated
covariance matrix ${\mathbf{C}}_{m_{1}}$, and another data point
${\mathbf{x}}_{m_{2}}$ from a dataset with $m_{2}$ features and associated
covariance matrix ${\mathbf{C}}_{m_{2}}$, we aim to characterize the
conditions under which the VNN outputs
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ and
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ converge. When
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ and
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ converge, we
can conclude that the parameters ${\mathcal{H}}$ are transferable between two
datasets consisting of $m_{1}$ and $m_{2}$ features.

Note that the VNN outputs $\Phi(\cdot;{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ and
$\Phi(\cdot;{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ have distinct
dimensionalities if $m_{1}\neq m_{2}$ and therefore, a direct comparison
between them is not natural. Fundamentally, it is imperative to provide a
mathematical framework to compare vectors and covariance matrices of different
sizes in order to be able to analyze the transferability of VNNs. To
facilitate such a comparison between vectors of different sizes, we consider a
simple mapping that represents the vector on a continuous interval $[0,1]$.
Specifically, given an $m$-dimensional vector
${\mathbf{x}}=[x_{1},\dots,x_{m}]$, we can define a continuous representation
of ${\mathbf{x}}$ as a function $y_{{\mathbf{x}}}:[0,1]\mapsto{\mathbb{R}}$,
such that, $y_{{\mathbf{x}}}(u)=x_{i}$ for $u\in{\mathcal{U}}_{i}$, where
${\mathcal{U}}_{i}$ is a pre-defined interval associated with the $i$-th
element of ${\mathbf{x}}$. Similarly, we can map a covariance matrix
${\mathbf{C}}_{m}$ to a compact space defined by $[0,1]^{2}$ using the mapping
${\mathbf{W}}_{{\mathbf{C}}_{m}}:[0,1]^{2}\mapsto{\mathbb{R}}$ where we have
${\mathbf{W}}_{{\mathbf{C}}_{m}}(u,v)=[{\mathbf{C}}_{m}]_{ij}$ for
$u\in{\mathcal{U}}_{i}$ and $v\in{\mathcal{U}}_{j}$. A pictorial illustration
of $y_{{\mathbf{x}}}$ and ${\mathbf{W}}_{{\mathbf{C}}}$ for covariance matrix
${\mathbf{C}}$ is included in Fig 2. Therein, the intervals
${\mathcal{U}}_{i}$ are parameterized by variables $\rho_{i}$, which will be
discussed subsequently in (20). Note that we can recover ${\mathbf{x}}$ from
$y_{{\mathbf{x}}}$ and vice-versa (similarly for ${\mathbf{C}}_{m}$ and
${\mathbf{W}}_{{\mathbf{C}}_{m}}$). Hence, for data points
${\mathbf{x}}_{m_{1}}$ and ${\mathbf{x}}_{m_{2}}$ consisting of $m_{1}$ and
$m_{2}$ elements, respectively, the closeness of continuous representations
$y_{{\mathbf{x}}_{m_{1}}}$ and $y_{{\mathbf{x}}_{m_{2}}}$ can be used as a
metric to assess the similarity between data points in multi-scale datasets.
This observation also extends to the comparison between covariance matrices
${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$.

Figure 2: Representations of $m$-dimensional vector ${\mathbf{x}}$ and
associated covariance matrix ${\mathbf{C}}$ in the continuous domain. a.
Representation for ${\mathbf{x}}$ is obtained by discretizing the space
$[0,1]$. b. Representation ${\mathbf{W}}_{{\mathbf{C}}}$ for ${\mathbf{C}}$ is
obtained by discretizing the space $[0,1]^{2}$ according to (20). Thus,
${\mathbf{W}}_{{\mathbf{C}}}$ retains the symmetry of ${\mathbf{C}}$. The area
spanned by the diagonal elements of ${\mathbf{C}}$ is marked in blue. The size
of the square area allotted to a diagonal element is proportional to its
value. Other parts of the grid accommodate the off-diagonal elements of
${\mathbf{C}}$.

For a VNN architecture with $F$ number of outputs in the final layer, the
dimensionality of VNN outputs
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ is
$m_{1}\times F$ and that for
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ is
$m_{2}\times F$. Thus, we can compare
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ and
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ via the
comparison between the continuous representations of every column in outputs
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ and
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$, where the
continuous representations are defined in the same fashion as
$y_{{\mathbf{x}}}$ above. For VNN
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$, we use the
notation $y_{m_{1}}[f]$ to denote the continuous representation of $f$-th
output in $\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$,
i.e.,
$y_{[\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})]_{f}}$.
Similar to the relationship between $y_{{\mathbf{x}}}$ and ${\mathbf{x}}$, the
$f$-th VNN output
$[\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})]_{f}$ and its
continuous representation $y_{m_{1}}[f]$ are operationally interchangeable
(see Appendix C for details). Using the continuous representations above, we
can now describe the assessment of transferability of VNNs more concretely.

Formal Problem Statement for VNN Transferability: Consider two VNNs
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ and
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ instantiated
on data with $m_{1}$ and $m_{2}$ features, respectively. If we have the
following conditions: (a) the continuous approximations of inputs
${\mathbf{x}}_{m_{1}}$ and ${\mathbf{x}}_{m_{2}}$ are close, i.e.,
$\|y_{{\mathbf{x}}_{m_{1}}}-y_{{\mathbf{x}}_{m_{2}}}\|_{2}$ is bounded; and
(b) the continuous approximations of covariance matrices
${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$ are close, i.e.,
$\|{\mathbf{W}}_{{\mathbf{C}}_{m_{1}}}-{\mathbf{W}}_{{\mathbf{C}}_{m_{2}}}\|_{2}$
is bounded; we aim to show that the continuous representations of VNN outputs
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ and
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ are close,
i.e., $\|y_{m_{1}}[f]-y_{m_{2}}[f]\|_{2}$ is bounded for all
$f\in\\{1,\dots,F\\}$.

Next, we informally state the main result of this section that establishes the
transferability between VNNs processing datasets consisting of $m_{1}$ and
$m_{2}$ features.

######  Theorem 3 (Transferability of VNNs).

Consider two VNNs
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ and
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ consisting of
$L$ layers and $F$ outputs per layer. If the continuous representations of
inputs and covariance matrices are close, i.e.,
$\|y_{{\mathbf{x}}_{m_{1}}}-y_{{\mathbf{x}}_{m_{2}}}\|_{2}$ and
$\|{\mathbf{W}}_{{\mathbf{C}}_{m_{1}}}-{\mathbf{W}}_{{\mathbf{C}}_{m_{2}}}\|_{2}$
are bounded, and assumptions A1 and A2 are satisfied (described in (20) and
(24)), we have

$\displaystyle\|y_{m_{1}}[f]-y_{m_{2}}[f]\|_{2}\leq
LF^{L}\beta\Big{(}\frac{1}{m_{1}^{3\zeta/2-1}}+\frac{1}{m_{2}^{3\zeta/2-1}}\Big{)}\;,\forall
f\in\\{1,\dots,F\\}\;,$ (19)

for some $\beta>0$ and $\zeta\in(2/3,1]$.

###### Proof.

See Appendix E. ∎

Theorem 3 implies that continuous representations of all $F$ outputs of the
respective final layers 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}})$ converge with
increase in $m_{1}$ and $m_{2}$. Since the continuous representation
$y_{m_{1}}[f]$ and VNN output
$[\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})]_{f}$ are
operationally interchangeable, we expect the measures of central tendency
(e.g., mean, median) of outputs
$[\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})]_{f}$ and
$[\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})]_{f}$ to
converge as well. By extension, we also expect the measures of central
tendency for $\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$
and $\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ to
converge if Theorem 3 holds. In this context, if the readout function for VNN
is unweighted mean, we expect the statistical outcomes 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}})$ to be close
and this convergence to be stronger for large $m_{1}$ and $m_{2}$. The impact
of Theorem 3 is broad, as we have shown that the parameters ${\mathcal{H}}$
can be “scale-free” while preserving the performance over an inference task.
Specifically, a VNN can be instantiated on a dataset of different
dimensionality than the training dataset and the VNN recovers close
statistical outcomes for the same parameters ${\mathcal{H}}$ for both
datasets, provided the data samples and covariance matrices of the training
dataset and the new dataset are close in terms of their continuous
representations. Thus, VNNs also offer a significant advantage over PCA-based
analysis approaches as the principal components are restricted within the
feature space of a dataset and do not provide any mathematical insight into
the structure of another dataset of different dimensionality even when the
datasets may be related. An overview of the transferability of VNNs is
illustrated in Fig. 3. While the rigorous details behind the proof of Theorem
3 are relegated to Appendix E, we briefly discuss the major mathematical
concepts and assumptions that have enabled us to establish (19) in Theorem 3
in the highlighted text on page 19.

Thus far, we have discussed the transferability of VNNs in terms of achieving
statistical outcomes (e.g., outcome of a regression model) that are close
using datasets of different dimensionalities, given that the datasets are
aligned in some way. As discussed previously, a multi-scale neuroimaging
dataset provides an ideal setting for validating the theoretical guarantees of
Theorem 3 in this context. However, note that the VNNs also provide
expressivity at the feature-level at the final layer. For instance, if VNN is
deployed for a regression task and the readout layer is a simple mean
function, the VNN final layer outputs could be leveraged to characterize the
contributions of each feature in the dataset to the final regression outcome.
In applications based on neuroimaging, this observation can be of great
interest, as each feature in a neuroimaging dataset is typically associated
with a distinct brain region. Thus, VNNs naturally provide a feasible way to
_interpret_ the final statistical outcomes. Moreover, if such an
interpretability is achieved using VNNs, we can further intertwine our
analysis with the transferability of VNNs to test its spatial robustness
across datasets of different dimensionalities. The observations made here
motivated us to pursue investigating the utility of VNNs in characterizing the
contributing brain regions to elevated $\Delta$-Age due to age-related
neurodegeneration.

Figure 3: Overview of transferability of VNNs. $y_{m_{1}}[f]$ is the
continuous representation of $f$-th output of VNN
$\Phi({\mathbf{x}}_{m_{1}};{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ that is
instantiated on data and covariance matrix with $m_{1}$ features. Similarly,
$y_{m_{2}}[f]$ represents the $f$-th output of VNN
$\Phi({\mathbf{x}}_{m_{2}};{\mathbf{C}}_{m_{2}},{\mathcal{H}})$ that is
instantiated on dataset with $m_{2}$ features. If the continuous counterparts
of covariance matrices ${\mathbf{C}}_{m_{1}}$ and ${\mathbf{C}}_{m_{2}}$,
i.e., ${\mathbf{W}}_{{\mathbf{C}}_{m_{1}}}$ and
${\mathbf{W}}_{{\mathbf{C}}_{m_{2}}}$, belong to a sequence that converges to
a graphon ${\mathbf{W}}$ (Definition 4) and the continuous representations of
inputs ${\mathbf{x}}_{m_{1}}$ and ${\mathbf{x}}_{m_{1}}$ are close, the
convergence between $y_{m_{1}}[f]$ and $y_{m_{2}}[f]$ is characterized in
terms of $m_{1}$ and $m_{2}$ in Theorem 3.

Mathematical Foundations of Transferability

The continuous representations of graph signals and graphs have previously
been leveraged to study transferability of GNNs under the domain of graphon
information processing [79]. Specifically, GNNs can be transferable between
graphs belonging to a converging sequence if the graphs in this sequence
converge to a limit object called _graphon_ as the number of nodes approaches
infinity [80]. In a similar fashion, we leverage the theory of graphons [80]
and graphon signal processing [79] to establish Theorem 3. Graphons are the
limits of _dense_ graphs (i.e., graphs with number of edges of the order
$\Theta(m^{2})$) [81] and hence, appropriate to study limits of covariance
matrices that are typically dense. The definition of a graphon is provided in
Definition 4.

######  Definition 4 (Graphon).

A graphon is a bounded, symmetric, measurable function
${{\mathbf{W}}:[0,1]^{2}\mapsto[-1,1]}$. Under the setting where a covariance
matrix is viewed as a weighted graph and normalized such that its largest
eigenvalue is $1$, we claim that a sequence of covariance matrices
$\\{{\mathbf{C}}_{m}\\}$ being convergent implies that the sequence of their
counterpart continuous representations, i.e.,
$\\{{\mathbf{W}}_{{\mathbf{C}}_{m}}\\}$, converges to some graphon
${\mathbf{W}}$ if ${\mathbf{W}}_{{\mathbf{C}}_{m}}$ is appropriately
constructed from ${\mathbf{C}}_{m}$. This claim is based on generalizing [80,
Corollary 3.9] to our setting, and the formal statement in this regard is
included in Remark 2 in Appendix C. Theorem 3 holds for any pair of covariance
matrices in the converging sequence $\\{{\mathbf{C}}_{m}\\}$ (under certain
assumptions that will be discussed shortly) and thus, parameters
${\mathcal{H}}$ can be transferred between any two VNNs instantiated on
distinct covariance matrices in this sequence. The construction of
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ relies on appropriately defining the
intervals ${\mathcal{U}}_{i}$ and is described in the following steps. a.
Partition the interval $[0,1]$ into $m$ disjoint intervals
$[{\mathcal{U}}_{1},\dots,{\mathcal{U}}_{m}]$, such that,
$\displaystyle{\mathcal{U}}_{i}=\begin{cases}[0,\rho_{1}]\quad\text{if}\quad
i=1\\\ (\rho_{i-1},\rho_{i}]\quad\text{if}\quad
i\in\\{2,\dots,m\\}\end{cases}\;,$ (20) where
$\displaystyle\rho_{i}\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}\frac{1}{{\sf
tr}({\mathbf{C}}_{m})}\sum\limits_{j=1}^{i}[{\mathbf{C}}_{m}]_{jj}\;,$ (21)
and ${\sf tr}({\mathbf{C}}_{m})$ is the trace of ${\mathbf{C}}_{m}$. Clearly,
$\rho_{m}=1$. b. The relationship between feature $i$ and feature $j$ is given
by ${\mathbf{W}}_{{\bf C}_{m}}(u_{i},u_{j})=[{\mathbf{C}}_{m}]_{ij}$ for
$u_{i}\in{\mathcal{U}}_{i},u_{j}\in{\mathcal{U}}_{j}$. If the continuous
representation ${\mathbf{W}}_{{\mathbf{C}}_{m}}$ of ${\mathbf{C}}_{m}$ is
constructed according to the above steps, we refer to
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ as the graphon approximation of
${\mathbf{C}}_{m}$. Thus, the graphon limit ${\mathbf{W}}$ forms the schema
for which the covariance matrix ${\mathbf{C}}_{m}$ represents the covariance
realization at resolution $m$. Next, we note that the result (19) is
contingent upon two main assumptions related to the covariance matrix
${\mathbf{C}}_{m}$ and the graphon limit ${\mathbf{W}}$. These assumptions are
mentioned below. A1. ($(\Omega,\zeta)$-dominant property of covariance
matrices) For the sequence $\\{{\mathbf{C}}_{m}\\}$, there exist positive
constants $\Omega$ and $\zeta$, such that, we have $\displaystyle\frac{1}{{\sf
tr}({\mathbf{C}}_{m})}\max_{j\in\\{1,\dots,m\\}}[{\mathbf{C}}_{m}]_{jj}\leq\frac{\Omega}{m^{\zeta}}\;,$
(22) for all finite $m$. Our analysis in Appendix E shows that Theorem 3 holds
for any $\Omega>0$ and $\zeta\in(2/3,1]$. The property in 22 implies that
$\displaystyle\frac{1}{{\sf
tr}({\mathbf{C}}_{m})}\max_{j\in\\{1,\dots,m\\}}[{\mathbf{C}}_{m}]_{jj}\rightarrow
0$ (23) as $m\rightarrow\infty$. We refer to the covariance matrix
${\mathbf{C}}_{m}$ satisfying the property in (22) as being
$(\Omega,\zeta)$-dominant. We also note that (23) suggests that the variance
profile of individual features in the dataset, characterized by their
corresponding diagonal elements in the covariance matrix, must not be
concentrated within a small subset of features. A2. (Lipschitz continuity of
Graphon.) If ${\mathbf{W}}$ is the limit of the sequence
$\\{{\mathbf{W}}_{{\mathbf{C}}_{m}}\\}$ as $m\rightarrow\infty$, then for
Theorem 3 to hold, ${\mathbf{W}}$ must satisfy
$\displaystyle|{\mathbf{W}}(u_{1},v_{1})-{\mathbf{W}}(u_{2},v_{2})|\leq\alpha(|u_{1}-u_{2}|+|v_{1}-v_{2}|)\;,$
(24) for any $u_{1},v_{1},u_{2},v_{2}\in[0,1]$. Any graphon satisfying (24) is
termed as an $\alpha$-Lipschitz graphon. The Lipschitz continuity of graphon
${\mathbf{W}}$ determines the smoothness of the information present between
any two coordinates $(u_{1},v_{1})$ and $(u_{2},v_{2})$. Therefore, it is
intuitively expected that for a given $m$, a graphon ${\mathbf{W}}$ with
smaller Lipschitz constant $\alpha$ will be better approximated by
${\mathbf{W}}_{{\mathbf{C}}_{m}}$ as compared to that with higher Lipschitz
constant.

Figure 4: Overview of brain age prediction framework using VNNs. The VNN is
trained as a regression model by optimizing filter taps ${\mathcal{H}}$ to
predict chronological age using cortical thickness features for HC group. The
VNN prediction is formed via a readout function that averages the outputs at
the final layer of VNN. For brain age prediction, we replace the covariance
matrix in the trained VNN model with the anatomical covariance matrix from the
combined population of HC group and individuals with AD. For each individual,
VNNs first form the scalar output by combining outputs at the final layer of
VNN via the readout function, and then bias-correction is applied to finally
evaluate the brain age. Statistical analysis of the final layer outputs of the
VNN is performed across the combined population to identify locally elevated
residual effects in neurodegeneration. These locally elevated residual effects
eventually contribute to the elevated $\Delta$-Age effect and hence, provide a
regional profile to elevated $\Delta$-Age due to age-related
neurodegeneration. If the transferability of performance on the chronological
age prediction task holds for different datasets (according to Theorem 3), we
expect to see similar $\Delta$-Age distributions across them without re-
training. However, the regional profiles do not hinge on the performance of
VNN in predicting chronological age of HC group and are derived before age-
bias correction is applied to calculate $\Delta$-Age for an individual. Hence,
if the regional profiles are characteristic of accelerated aging, we expect to
observe similar regional profiles after transferring VNNs to different
datasets even when Theorem 3 may not hold (for instance when transferring the
VNN model from FTDC datasets to OASIS-3 dataset).

##  3 Application: Brain Age Prediction

$\Delta$-Age is a known biomarker of cognitive decline and neurodegeneration
[82, 14]. However, in the absence of a ground truth, the notion of
$\Delta$-Age is abstract and has a limited clinical utility without
identification of the main contributors to the elevated brain age due to
neurodegeneration. In this paper, we focus on Alzheimer’s disease (AD) as an
example of pathological age-related neurodegeneration. Age is a major risk
factor for AD and hence, AD is characterized by biological traits that signify
accelerated aging [83]. We leverage the architecture of VNNs to provide a
regional perspective to brain age prediction and our results demonstrate that
the elevated $\Delta$-Age in AD is accompanied by abnormalities in various
regions of interest characteristic of AD. The description of datasets is
included in subsection 4.1 in the Methods section. Next, we discuss how the
analysis of outputs at the final layer of VNNs may capture the impact of age-
related neurodegeneration on various regions of the brain.

###  3.1 Identification of Regions Associated with Neurodegeneration using VNN
Architecture

To start with, we used VNN as a regression model to fit multivariate cortical
thickness data against chronological age for healthy controls (HC). Details on
training and selected architecture of VNN models used for the brain age
analysis are provided in the subsection 4.2 in the Methods section. Thus, in
our experiments, a VNN processed the cortical thickness data through a multi-
layer architecture consisting of the anatomical covariance matrix derived from
cortical thickness of the HC group. The final prediction by the VNN model was
determined by unweighted aggregation of the final layer outputs which can be
conceptualized as an unweighted mean of age predictions at individual brain
regions. Therefore, VNN architecture allowed us to compute “regional
residuals” (scalar output at a given region derived from VNN final layer
output - aggregated VNN output or age estimate formed by VNN) at each brain
region to assess their contribution to the final output of VNN. We refer the
reader to subsection 4.4.1 in the Methods section for concrete details
regarding calculations of regional residuals.

The VNN model trained on the data from HC group according to the procedure in
subsection 4.2 captured the information about healthy aging from cortical
thickness and associated anatomical covariance matrix. Note that this
information may not be enough to predict chronological age accurately, i.e.,
we may not achieve a perfect fit on chronological age of HC group. However,
our theoretical results regarding stability of VNNs in Theorem 2 and
transferability of VNNs in Theorem 3 dictate that the regression performance
achieved by VNNs in this context is expected to be stable to perturbations in
the anatomical covariance matrix due to change in sample size and transferable
across datasets of different dimensionalities if their respective covariance
matrices satisfy certain theoretical assumptions.

Next, we tested these trained VNN models on the combined groups of HC and AD
by replacing the anatomical covariance matrix from the HC group with the
anatomical covariance matrix from the combined group of HC and individuals
with AD. Since capturing the accelerated aging is of interest, we hypothesized
that the brain regions characteristic of AD (for instance, regions with higher
cortical atrophy) with respect to HC group were likely to have elevated
regional residuals for individuals with AD. Also, since outputs of VNNs are
obtained via manipulation of the data according to principal components or
eigenvectors of the covariance matrix (Section 2.1), we hypothesized that the
regional residuals may be correlated with the principal components of the
anatomical covariance matrix of the combined group. Furthermore, the stability
of VNNs established in Theorem 2 implied that any group differences observed
in regional residuals were expected to be stable to perturbations in the
covariance matrix and thus, stable to the composition of the combined
population of HC and individuals with AD used for estimating the anatomical
covariance matrix.

The scale-free architecture of VNNs allowed us to gauge the spatial and
statistical robustness of any observed elevated regional residuals in age-
related neurodegeneration in scenarios when the transferability of performance
was guaranteed (Theorem 3 and Fig. 3) and when it was not guaranteed. Our
discussion in Appendix F demonstrates that the cortical thickness datasets
curated according to different parcellations of Schaefer’s atlas (referred to
as FTDC datasets and described in Section 4.1) laid within the scope of
theoretical guarantees on transferability as described in Fig. 3. By
extension, the performance of VNN models trained for chronological age
prediction of HC group was transferable across datasets organized according to
different resolutions of Schaefer’s atlas. Hence, assessments of the regions
pertaining to elevated regional residuals across different scales of
Schaefer’s atlas could provide a proof of concept for the transferability
property holding in brain age prediction. Experiments pertaining to this
aspect are included in Appendix K.

We focus on the setting where the theoretical guarantees of transferability of
performance (Theorem 3) may not hold to assess whether transferring VNNs
across datasets curated according to distinct brain atlases resulted in
retaining the regional profiles pertinent to AD-related neurodegeneration. In
this context, we studied the regional profiles when VNNs were transferred from
a dataset organized according to Schaefer’s atlas to a dataset organized
according to Destrieux (DKT) [84] atlas. The datasets in this context were
derived from distinct populations. The dataset organized according to DKT
atlas laid outside the purview of theoretical guarantees on transferability of
VNNs trained on datasets organized according to different scales of Schaefer’s
atlas (see Fig. 13 in Appendix F and associated discussion). Lack of
theoretical guarantees on the VNN transferability between datasets organized
according to DKT atlas and those according to Schaefer’s atlas imply that VNNs
may need to be augmented by a mapping that accounts for the differences
between DKT atlas and Schaefer’s atlas to achieve comparable performance on
the chronological age prediction task. Hence, while the performance on the
task of chronological age prediction was not guaranteed to be trivially
transferable between datasets curated according to Schaefer’s atlas and DKT
atlas, we hypothesized that the observed effects of elevated regional
residuals may be observed even after transferring VNNs if such effects were
driven by cortical atrophy or neurodegeneration-related features in the
anatomical covariance matrix. Therefore, evaluating the regional residuals
after transferring VNNs between datasets curated according to distinct brain
atlases could allow decoupling of the potential contributors of elevated brain
age from the objective of achieving near perfect performance over
chronological age prediction in the HC group 333We remark that the discussion
here is not restricted to Schaefer’s atlas and DKT atlas and could potentially
be extended to any pair of datasets curated according to distinct brain
atlases, such that, VNNs may not be trivially transferable between them..

We further expanded the experiments above by investigating the correlations
between regional residuals and clinical dementia rating (CDR) metrics. CDR sum
of boxes scores are commonly used in clinical and research settings to stage
dementia severity [85]. A higher CDR score is associated with more severe
cognitive and functional status. If our hypothesis that the elevated regional
residuals were driven by age-related neurodegeneration was valid, we expected
to observe a significant alignment between the span of brain regions with
elevated regional residuals and the span of brain regions whose regional
residuals were correlated with CDR metrics.

###  3.2 Individual-level Brain Age Prediction

Next, we focused on predicting the $\Delta$-Age for an individual from VNN
outputs using a procedure that is consistent with other studies in the
literature. The residuals evaluated via the difference between the
chronological age and predicted age for individuals in the HC group by VNNs
are prone to bias. Specifically, age for younger individuals tends to be over-
estimated and that for older individuals tends to be under-estimated [86, 87].
This phenomenon may appear, for example, when the correlation between the
chronological age and predicted age is significantly smaller than 1. In order
to correct for this age-related bias in the residuals, we used a linear
regression based approach [88] and evaluated $\Delta$-Age. We hypothesized
$\Delta$-Age to be significantly elevated in the group of individuals with AD
with respect to that for those in the HC group. The methodological details on
age-bias correction and the procedure to evaluate $\Delta$-Age are included in
subsection 4.4.2 in the Methods section. Also, by extension of the
transferability property of VNNs, we have previously reported similar
distributions of $\Delta$-Age across datasets curated according to different
scales of a multi-scale brain atlas in [2].

##  4 Methods

###  4.1 Data

We consider datasets from two independent populations in this paper, as
described below.

Multi-scale FTDC Datasets. These datasets consist of the cortical thickness
data extracted at different resolutions from healthy controls (HC; $n=105$,
age = $62.6\pm 7.62$ years, 57 females). For each individual, the cortical
thickness data was curated according to multi-resolution Schaefer atlas [9],
at 100 parcel, 300 parcel, and 500 parcel resolutions with finer resolution
cortical thickness estimates with increasing number of parcellations. The ANTs
cortical thickness pipeline [89, 90] was used to derive mean cortical
thickness within each atlas parcel using 3T T1-weighted MRIs ( 1mm isotropic
resolution). We report results on three datasets: FTDC100, FTDC300 and
FTDC500, that constitute the cortical thickness datasets corresponding to 100,
300 and 500 cortical thickness features, respectively. Also, the FTDC100,
FTDC300, and FTDC500 datasets are jointly referred to as FTDC datasets.

OASIS-3 Dataset. This dataset was derived from publicly available freesurfer
estimates of cortical thickness (hosted on central.xnat.org), as previously
reported [64], and comprised of cognitively normal individuals (HC; $n=652$,
age = $67.76\pm 7.88$ years, $382$ females) and individuals with AD dementia
diagnosis and at various stages of cognitive decline ($n=209$, age = $74.61\pm
7.13$ years, $102$ females). The cortical thickness features were curated
according to the Destrieux (DKT) atlas [84](consisting of $148$ cortical
regions). In the context of transferability, OASIS-3 provided a dataset
curated according to a distinct brain atlas than the multi-scale Schaefer’s
atlas for FTDC datasets and hence, allowed us to investigate the VNN
transferability beyond the setting with multi-resolution datasets. For clarity
of exposition of the brain age prediction method, any dementia staging to
subdivide the group of individuals with AD dementia diagnosis into mild
cognitive impairment (MCI) or AD was not performed and we use the label AD+ to
refer to this group. The individuals in AD+ group were significantly older
than those in HC group (t-test: $p$-value = $4.15\times 10^{-27}$). The
boxplots for the distributions of chronological age for HC and AD+ groups are
included in Fig. 16. For $206$ individuals in the AD+ group, the CDR sum of
boxes scores evaluated within one year (365 days) from the MRI scan were
available (CDR sum of boxes = $3.38\pm 1.73$). The CDR sum of boxes scores for
this population were evaluated according to [91].

The data from individuals in HC group for FTDC and OASIS-3 were leveraged to
investigate transferability of VNNs as governed by the theoretical framework
provided in Section 2. Furthermore, we primarily used the OASIS-3 dataset to
derive conclusions about $\Delta$-Age as a marker of accelerated aging in the
AD+ group. The observations made after transferring VNNs between FTDC and
OASIS-3 datasets were used to gain methodological insights regarding the role
of training VNNs to predict chronological age in the brain age prediction
framework. Besides the datasets described above, we also extended the
transferability aspect of VNNs to the brain age prediction framework on a
small, multi-scale dataset consisting of individuals with AD diagnosis
(independent of OASIS-3) that was curated according to different scales of the
Schaefer’s atlas (Appendix K) in order to demonstrate the spatial robustness
of the regional profiles obtained by our brain age prediction framework across
multiple scales and validation of the findings on OASIS-3 dataset.

###  4.2 VNN Learning

We use the VNN model for a regression task where a multi-variate feature set
is regressed to a scalar quantity. Note that the VNN output of the
architecture represented by
$\Phi({\mathbf{x}};\hat{\mathbf{C}},{\mathcal{H}})$444We use the notation
$\hat{\mathbf{C}}$ for covariance matrix in the Methods section as the VNN
architecture is modeled on the sample covariance matrix in practical
implementation. for one $m$-dimensional input is of dimension $m\times F$ if
the VNN architecture has $F$ $m$-dimensional outputs in the final layer. The
regression output is determined by a readout layer which evaluates an
unweighted mean of all outputs of the final layer of VNN. Therefore, the
regression output for an input ${\mathbf{x}}$ is given by

$\displaystyle\hat{y}=\frac{1}{Fm}\sum\limits_{j=1}^{m}\sum\limits_{f=1}^{F}[\Phi({\mathbf{x}};\hat{\mathbf{C}},{\mathcal{H}})]_{jf}\;.$
(25)

Prediction using unweighted mean at the output implies that the VNN model
exhibits permutation-invariance (i.e., the final output is independent of the
permutation of the input features and covariance matrix) and transferability.
Although the performance of the VNN model could potentially be improved by
adding a learnable or an adaptive readout function (e.g., weighted mean or a
single layer perceptron) that maps the final layer outputs of VNN to scalar
$\hat{y}$ via a learnable mapping, our subsequent experiments implied a
negative impact on the interpretability of VNN model in the context of brain
age prediction. For a regression task, the training dataset
$\\{{\mathbf{x}}_{i},y_{i}\\}_{i=1}^{n}$ is leveraged to learn the filter taps
in ${\mathcal{H}}$ for the VNN $\Phi(\cdot;\hat{\mathbf{C}},{\mathcal{H}})$
such that they minimize the training loss, i.e.,

$\displaystyle{\mathcal{H}}_{\sf
opt}=\min_{{\mathcal{H}}}\frac{1}{n}\sum\limits_{i=1}^{n}\ell(\hat{y}_{i},y_{i})\;,$
(26)

where

$\displaystyle\hat{y}_{i}=\frac{1}{Fm}\sum\limits_{j=1}^{m}\sum\limits_{f=1}^{F}[\Phi({\mathbf{x}}_{i};\hat{\mathbf{C}},{\mathcal{H}})]_{jf}\;,$
(27)

and $\ell(\cdot)$ is the mean-squared error (MSE) loss function.

In our experiments, we trained four sets of VNN models; one each for the HC
group in FTDC100, FTDC300, FTDC500, and OASIS-3 datasets. The training process
was similar for all VNNs. We randomly split the dataset into an approximately
$90/10$ train/test split. Thus, the test sets for FTDC datasets consisted of
$10$ individuals and that in OASIS-3 dataset consisted of $65$ individuals.
The sample covariance matrix was evaluated using all samples in the training
set ($n=95$ for FTDC datasets and $n=587$ for OASIS-3 dataset) and we had the
sample covariance matrix $\hat{\mathbf{C}}$ of size $m\times m$ (where $m=100$
for FTDC100, $m=300$ for FTDC300, $m=500$ for FTDC500, and $m=148$ for
OASIS-3). Furthermore, for all datasets, $\hat{\mathbf{C}}$ was normalized
such that its maximum eigenvalue was $1$. Next, the training set was randomly
split internally, such that, the VNN was trained with respect to the mean
squared error loss between the predicted age and the true age in $n=84$
samples for FTDC datasets and $n=513$ samples for OASIS-3 dataset. The loss
was optimized using batch stochastic gradient descent with Adam optimizer
available in PyTorch library [92] for up to $100$ epochs. The batch size was
$34$ for FTDC100 dataset, $8$ for FTDC300 dataset, $12$ for FTDC500 dataset,
and $78$ for OASIS-3 dataset. The VNN model with the best minimum mean squared
error performance on the remaining samples in the training set (which acted as
a validation set) was included in the set of nominal models for this
permutation of the training set. For each dataset, we trained and validated
the VNN models over $100$ permutations of the complete training set of $n=95$
samples for each of the FTDC datasets and $n=587$ samples for OASIS-3 dataset,
thus, leading to $100$ trained VNN models (also referred to as nominal models)
per dataset.

The hyperparameters for the VNN architecture and learning rate of the
optimizer were chosen according to a hyperparameter search procedure using the
package Optuna [93]. For FTDC100, the VNN had a $L=2$-layer architecture, with
a filter bank such that we had $F=26$ and $2$ filter taps in each layer. The
learning rate for the optimizer was $0.058$. The number of learnable
parameters for this VNN was $1456$. For FTDC300, the VNN had a $L=2$-layer
architecture, with a filter bank such that we had $F=39$ and $3$ filter taps
in the first layer and $2$ filter taps in the second layer. The learning rate
for the optimizer was $0.0241$. The number of learnable parameters for this
VNN was $3237$. For FTDC500, the VNN model had a $L=2$-layers with a filter
bank such that we had $F=27$ and $4$ filter taps in the first layer and $2$
filter taps in the second layer. The number of learnable parameters for this
VNN was $1620$. The learning rate for the Adam optimizer was set to $0.0631$.
For OASIS-3, the VNN model had a $L=2$-layers with a filter bank such that we
had $F=5$ and $6$ filter taps in the first layer and $10$ filter taps in the
second layer. The learning rate for the Adam optimizer was set to $0.059$. The
number of learnable parameters for this VNN was $290$.

Note that the parameters of all VNNs were learned on the HC group in their
respective datasets and no subsequent training was performed for brain age
prediction. The above procedure for training the VNNs by splitting the
datasets into training/validation/test sets was performed to ensure that the
VNNs were not overfit on the training set. The performance of VNNs on test
sets for different datasets are reported in Appendix G. However, our results
primarily focus on the settings of transferability and brain age prediction,
both of which are reported on the complete datasets in the results section.

###  4.3 Transferability of VNNs

In our experiments, we empirically studied the transferability of the VNN
models across different datasets described in Section 4.1. In general,
transferring a VNN model from dataset A to dataset B implies that the VNN was
trained for an inference task on dataset A and is being tested for the same
inference task on dataset B. The scale-free aspect of VNN architecture allows
transferring of VNNs between datasets of different dimensionalities.

We first describe the procedure for evaluating VNN transferability on multi-
scale cortical thickness datasets derived from the same population. Let
$\hat{\mathbf{C}}_{m_{1}}$ denote the sample covariance matrix of size
$m_{1}\times m_{1}$ from the dataset on which the VNN was trained. Here,
$\hat{\mathbf{C}}_{m_{1}}$ was normalized such that its maximum eigenvalue was
$1$ in order to reconcile with the definition of graphon in Definition 4.
Consider an individual whose chronological age is $y$ and cortical thickness
data available as ${\mathbf{x}}_{m_{1}}$ of size $m_{1}\times 1$ and
${\mathbf{x}}_{m_{2}}$ of size $m_{2}\times 1$. Therefore, the VNN model
prediction for the individual formed by the trained VNN model is given by

$\displaystyle\hat{y}_{1}=\frac{1}{Fm_{1}}\sum\limits_{j=1}^{m_{1}}\sum\limits_{f=1}^{F}[\Phi({\mathbf{x}}_{m_{1}};\hat{\mathbf{C}}_{m_{1}},{\mathcal{H}})]_{jf}\;,$
(28)

where the set of filter taps ${\mathcal{H}}$ are determined by training on the
cortical thickness dataset with $m_{1}$ features and we use the notation
$\Phi({\mathbf{x}}_{m_{1}};\hat{\mathbf{C}}_{m_{1}},{\mathcal{H}})$ to denote
the output at the final layer of the VNN. To empirically validate the
theoretical results on the transferability of VNNs, we aim to evaluate the
change in performance when the covariance matrix $\hat{\mathbf{C}}_{m_{1}}$ is
replaced with $\hat{\mathbf{C}}_{m_{2}}$, where $\hat{\mathbf{C}}_{m_{2}}$ is
generated from dataset with $m_{2}$ cortical thickness features. Therefore,
the predicted age based on the cortical thickness features
${\mathbf{x}}_{m_{2}}$ for the same individual is given by

$\displaystyle\hat{y}_{2}=\frac{1}{Fm_{2}}\sum\limits_{j=1}^{m_{2}}\sum\limits_{f=1}^{F}[\Phi({\mathbf{x}}_{m_{2}};\hat{\mathbf{C}}_{m_{2}},{\mathcal{H}})]_{jf}\;,$
(29)

where the filter taps ${\mathcal{H}}$ are the same as those in (28). If the
VNN model was transferable from a dataset with $m_{1}$ cortical thickness
features to that with $m_{2}$ features, we expected the predictions
$\hat{y}_{1}$ and $\hat{y}_{2}$ to be close.

The procedure to investigate VNN transferability between datasets of different
dimensionalities and distinct populations is similar as above with the
following extensions. In principle, to evaluate the transferability from
dataset B to dataset A, we compared the performance of VNNs that were trained
on dataset A and the VNNs that were transferred from dataset B to dataset A.
Clearly, if the VNNs that were transferred from dataset B to dataset A
achieved comparable performance as that of the VNNs that were trained on
dataset A, we could conclude that VNNs exhibited transference from dataset B
to dataset A.

###  4.4 Brain Age Prediction

The VNN models trained as regression models for predicting chronological age
using the cortical thickness data from healthy controls were expected to
capture the effect of healthy aging in cortical thickness features. Next, we
describe our strategy to characterize brain age with regional interpretability
in the context of AD. We focus on the brain age results determined for
OASIS-3. Additional observations made on multi-scale datasets in the context
of brain age are included in Appendix K.

####  4.4.1 Regional Analyses of VNN Residuals in Age-Related
Neurodegeneration

The VNN architecture allows us to associate a scalar output with each
dimension among the $m$ dimensions in the final layer. Specifically, we have

$\displaystyle{\mathbf{p}}_{i}=\frac{1}{F}\sum\limits_{f=1}^{F}[\Phi({\mathbf{x}}_{i};\hat{\mathbf{C}}_{m},{\mathcal{H}})]_{f}\;,$
(30)

where ${\mathbf{p}}_{i}$ is the vector denoting the mean of all final layer
outputs associated with filters in the filter bank at the final layer. Note
that the mean of all elements of ${\mathbf{p}}_{i}$ is the prediction
$\hat{y}_{i}$ formed in (28). In the context of cortical thickness datasets,
we can associate each element of ${\mathbf{p}}_{i}$ with a distinct brain
region. Therefore, the vector ${\mathbf{p}}_{i}$ is a vector of “regional
contributions” to the output $\hat{y}_{i}$ by the VNN. The parameters
${\mathcal{H}}$ were learnt over the HC group as described previously and kept
unchanged in the subsequent analyses. The subsequent details in this section
reflect that we primarily focused on the OASIS-3 dataset to study brain age.
For OASIS-3 dataset, we use the notation $\hat{\mathbf{C}}_{148}$ for the
covariance matrix formed by the cortical thickness features from HC group.

Next, we leveraged (30) to study the effect of neurodegeneration on brain
regions. For this purpose, in the OASIS-3 dataset, we evaluated the covariance
matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$ from the combined cortical thickness
data of HC and AD+ groups. As a consequence of the stability of Theorem 3, we
expect the inference drawn from VNNs to be stable to be stable to
perturbations in the covariance matrix. Therefore, subsequent evaluations
using VNNs are expected to be stable to the composition of combined HC and AD+
groups used to estimate the anatomical covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$.

For every individual in the combined dataset of HC and AD+ groups, we
processed their cortical thickness data ${\mathbf{x}}$ through the model
$\Phi({\mathbf{x}};\hat{\mathbf{C}}_{148}^{\sf AD+},{\mathcal{H}})$ where
parameters ${\mathcal{H}}$ were learnt in the regression task on the data from
HC group as described previously. Hence, the vector of mean of all final layer
outputs for cortical thickness input ${\mathbf{x}}$ is given by

$\displaystyle{\mathbf{p}}=\frac{1}{F}\sum\limits_{f=1}^{F}[\Phi({\mathbf{x}};\hat{\mathbf{C}}_{148}^{\sf
AD+},{\mathcal{H}})]_{f}\;,$ (31)

and VNN output is given by

$\displaystyle\hat{y}=\frac{1}{148}\sum\limits_{j=1}^{148}[{\mathbf{p}}]_{j}\;.$
(32)

Furthermore, we define the residual for feature $a$ (or brain region
represented by feature $a$ in this case) as

$\displaystyle[{\mathbf{r}}]_{a}\stackrel{{\scriptstyle\scriptscriptstyle\triangle}}{{=}}[{\mathbf{p}}]_{a}-\hat{y}\;.$
(33)

Thus, equation 33 allows us to characterize the residuals with respect to the
VNN output $\hat{y}$ at the regional level for individual brain regions for an
individual with cortical thickness data ${\mathbf{x}}$. Henceforth, we refer
to the residuals evaluated according to (33) as “regional residuals”.
Therefore, when an individual was predicted to have a brain age higher than
their chronological age due to neurodegenerative condition, we hypothesized
such an observation to be an aggregated effect of contributions from certain
biologically plausible brain regions. The brain regions contributing to the
observed higher brain age could be characterized at a regional level by the
analysis of regional residuals as defined in (33). Thus, the elements of the
residual vector ${\mathbf{r}}$ can potentially act as a biomarker that can
enable the isolation of brain regions affected due to age-related
neurodegeneration.

In our experiments, for a given VNN model, we evaluated the residual vector
${\mathbf{r}}$ for every individual in the dataset. Also, for experiments on
OASIS-3 dataset, we denote the population of residual vectors for healthy
controls as ${\mathbf{r}}_{\sf HC}$, individuals in AD+ group as
${\mathbf{r}}_{\sf AD+}$. The length of the residual vectors is the same as
the number of cortical thickness features in the dataset. Since each element
of the residual vector is associated with a distinct brain region, we
performed ANOVA to test for group differences between individuals in HC and
AD+ groups. Also, since elevation in $\Delta$-Age is the biomarker of interest
in this analysis, we hypothesized that the brain regions that exhibited higher
means for regional residuals for AD+ group than HC group would be the most
relevant to capturing accelerated aging. Hence, we report our results only for
brain regions that showed elevated regional residual distribution in AD+ group
with respect to HC group. Further, the group difference between AD+ and HC
groups in the residual vector element for a brain region was deemed
significant if it met the following criteria: i) the corrected $p$-value
(Bonferroni correction) for the clinical diagnosis label in the ANOVA model
was smaller than $0.05$, and ii) the uncorrected $p$-value for clinical
diagnosis label in ANCOVA model with age and sex as covariates was smaller
than $0.05$. An example for this regional analysis is included in Appendix H.

Recall that 100 distinct VNN models were trained as regression models on
different permutations of the training set of cortical thickness features from
HC group. We leveraged these trained models to establish the robustness of
observed group differences in the distributions of regional residuals. This
procedure is described next.

Robustness of findings from regional analyses. We performed the regional
analysis described above corresponding to each trained VNN model and tabulated
the number of VNN models for which a brain region was deemed significant in
the regional analysis described above. A higher number of VNN models isolating
a brain region as significant suggested higher robustness of the effect
observed for that brain region. For instance, a brain region with that
exhibited elevated regional residual in AD+ group with respect to HC group
across nearly all the 100 trained VNN models was likely to be a highly robust
contributor to elevated $\Delta$-Age in AD+ group. We used the fsbrain package
in R to project the robustness of significantly elevated regional residuals
for a brain region on the brain template [94].

The scale-free architecture of VNNs is facilitated by the coVariance filter
(Definition 2). Also, the non-adaptive readout function allows us to transfer
VNNs to process datasets of different dimensionalities without any changes to
the architecture. We used these properties of our proposed framework to
investigate “transferability of interpretability” across different datasets
irrespective of whether they laid within the purview of theoretical guarantees
on the transferability of performance of VNNs. This aspect is described next.

Transferability of interpretability. Note that the analysis of regional
residuals is independent of the performance of VNN model in the regression
task. Hence, we leveraged the scale-free aspect of VNNs to evaluate whether
the statistical effects that facilitated certain brain regions to be deemed
significant in the above analyses were preserved after transferring to a
dataset of different dimensionality. For this purpose, we transferred a VNN
model trained on FTDC datasets to OASIS-3 dataset and plotted the brain
regions deemed significant from the analyses of regional residuals of AD+ and
HC groups in OASIS-3 on a brain template. If accelerated aging characteristic
of AD was the driving factor behind the observed statistical differences in
regional residuals of AD+ and HC groups, we expected to observe qualitatively
consistent findings on OASIS-3 dataset for VNNs trained on OASIS-3 dataset and
those transferred from FTDC datasets to OASIS-3 dataset. This observation
would provide evidence for “transferability of interpretability” for VNNs in
brain age across different datasets as VNN models trained on FTDC datasets
could transfer their ability to construct interpretable regional residuals
from FTDC to OASIS-3. Here, we clarify that qualitatively consistency refers
to observing significant group differences between AD+ and HC groups in the
same direction (but not necessarily the same effect size) for regional
residuals obtained from VNNs trained on OASIS-3 and FTDC datasets. By
investigating transferability of interpretability, we could also potentially
decouple the contributing factors behind the observed accelerated aging effect
in AD+ group from the expectation of achieving a near-perfect performance on
chronological age prediction for HC group in existing literature [16, 15]. We
hypothesized that the transferability of interpretability could be driven by a
combination of cortical atrophy and the ability of VNNs to exploit
eigenvectors (principal components) of the covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ that were most relevant to brain age
prediction.

####  4.4.2 Individual-level Brain Age Prediction

Existing studies primarily focus on the gap between estimated brain age and
the chronological age as a biomarker of pathology [12], where the brain age is
a scalar quantity. VNNs provide an interpretable and systemic approach to
brain age via regional analyses as described in Section 4.4.1. However, using
VNNs, we can also obtain a scalar estimate for the brain age through a
procedure consistent with the existing studies in this domain. To evaluate the
brain age from VNN regression output, we first addressed the systemic bias in
the gap between $\hat{y}$ and $y$, where the age may be underestimated for
older individuals and overestimated for younger individuals [86]. Such bias
can confound the interpretations of brain age. Therefore, to correct for this
age-driven bias, we adopted a linear regression model based approach to
correct any age bias in the VNN age estimates [88, 86]. Under this approach,
we followed the following bias correction steps on the VNN estimated age
$\hat{y}$ to obtain the brain age $\hat{y}_{\sf B}$ for an individual with
chronological age $y$ and cortical thickness data ${\mathbf{x}}$:

Step 1. Fit a linear regression model on the complete dataset to determine
coefficients $\alpha$ and $\beta$ in the following linear model:

$\displaystyle\hat{y}-y=\alpha y+\beta\;.$ (34)

Step 2. Evaluate brain age as follows:

$\displaystyle\hat{y}_{\sf B}=\hat{y}-(\alpha y+\beta)\;.$ (35)

The gap between $\hat{y}_{\sf B}$ and $y$ is typically the biomarker of
interest and is defined below. For an individual with cortical thickness
${\mathbf{x}}$ and chronological age $y$, the brain age gap $\Delta$-Age is
defined as

$\displaystyle\Delta\text{-Age}\triangleq\hat{y}_{\sf B}-y\;,$ (36)

where $\hat{y}_{\sf B}$ is determined from the VNN age estimate $\hat{y}$ from
$\Phi({\mathbf{x}};{\mathbf{C}},{\mathcal{H}})$ and $y$ according to steps in
(34) and (35). The age-bias correction in (34) and (35) was performed for only
HC group to account for bias in the VNN estimates due to healthy aging and
then applied to the AD+ group. Further, the distributions of $\Delta$-Age were
obtained for all individuals in HC and AD+ groups. $\Delta$-Age for AD+ group
was expected to be elevated as compared to HC group as a consequence of
elevated regional residuals derived from the VNN model. To elucidate this, we
consider a toy example where we have two individuals of the same chronological
age $y$ with one belonging to the AD+ group and another to the HC group.
Equation (35) suggests that their corresponding VNN outputs (denoted by
$\hat{y}_{\sf AD+}$ for individual in the AD+ group and $\hat{y}_{\sf HC}$ for
individual in the HC group) are corrected for age-bias using the same term
$\alpha y+\beta$. Hence, $\Delta$-Age for the individual in the AD+ group will
be elevated with respect to that from the HC group only if the VNN prediction
$\hat{y}_{\sf AD+}$ is elevated with respect to $\hat{y}_{\sf HC}$. Since the
VNN predictions $\hat{y}_{\sf AD+}$ and $\hat{y}_{\sf HC}$ can be perceived as
unweighted aggregations of the estimates at the regional level (see (32)),
higher $\hat{y}_{\sf AD+}$ with respect to $\hat{y}_{\sf HC}$ is a direct
consequence of the regional residuals (see (33)) being elevated in AD+ group
with respect to HC group. When the individuals in this example have different
chronological age, the age-bias correction is expected to remove any variance
due to chronological age in $\Delta$-Age. Since the age distributions of AD+
and HC group are different, we also verified that the differences in
$\Delta$-Age for AD+ and HC group were not driven by age or gender differences
via ANCOVA with age and sex as covariates.

##  5 Results

###  5.1 Transferability of VNNs for regression task

We evaluated the transferability of VNN models trained as regression models on
the cortical thickness data of the HC group from FTDC datasets across
different resolutions of Schaefer’s atlas. To begin with, we remark that the
series of covariance matrices formed by cortical thickness features extracted
according to $100-500$ parcellation for HC group in FTDC datasets was
converging (Appendix F). This assessment was pertinent as our theoretical
results in Theorem 3 hold for a converging sequence of covariance matrices.
Furthermore, we also evaluated the distance of covariance matrices from FTDC
datasets from the covariance matrices derived from the HC group of OASIS-3
dataset ($\hat{\mathbf{C}}_{148}$ for DKT atlas). Here, the distances between
the covariance matrix derived from OASIS-3 dataset and those from FTDC
datasets were significantly greater than the pairwise distances for the
covariance matrices associated with different resolutions of Schaefer’s atlas
(Fig. 13 in Appendix F). This observation could potentially be explained by
the differences in constructions of Schaefer’s atlas and DKT atlas.
Specifically, DKT atlas is an anatomic atlas that maps both gyral and sulcal
regions [84]. In contrast, Schaefer’s atlas is derived from a gradient
weighted Markov random field model using functional MRI data [9] that does not
have anatomic boundaries per se. Due to these significant differences in
construction, we did not expect the anatomical covariance matrices from the
data curated according to DKT atlas to be a part of the converging sequence of
covariance matrices from FTDC datasets. Thus, based on the discussion thus
far, we expected the transferability property of VNNs to hold for statistical
inference tasks on FTDC datasets. We also performed exploratory analysis for
transferability from FTDC to OASIS-3 dataset for VNNs trained on FTDC datasets
and vice-versa to assess the degradation in regression performance.

In order to investigate transferability of VNNs, we trained the VNN models for
a regression task between cortical thickness and chronological age for
individuals in HC group according to the procedure described in Section 4.2
for FTDC100, FTDC300, FTDC500, and OASIS-3 datasets. This resulted in $100$
nominal VNN models for each dataset (each model trained on a different
permutation of the training set). The readout layer in the VNNs was non-
adaptive and it evaluated the unweighted mean of the outputs of the final VNN
layer to form an estimate for chronological age. Therefore, the trained VNN
could readily process a dataset with different number of features without any
retraining or alteration in the architecture. The performance outcomes were
quantified in terms of mean absolute error (MAE) and Pearson’s correlation
between the VNN output and ground truth.

Table 1: Transferability across datasets (MAE for VNN regression outputs with respect to the ground truth). Training Testing | FTDC100 (HC) | FTDC300 (HC) | FTDC500 (HC) | OASIS-3 (HC)
---|---|---|---|---
FTDC100 (HC) | $5.39\pm 0.084$ | $5.5\pm 0.101$ | $5.61\pm 0.132$ | $7.55\pm 0.356$
FTDC300 (HC) | $5.39\pm 0.193$ | $5.41\pm 0.167$ | $5.47\pm 0.169$ | $7.22\pm 0.4$
FTDC500 (HC) | $5.43\pm 0.2$ | $5.38\pm 0.15$ | $5.4\pm 0.169$ | $7.14\pm 0.6$
OASIS-3 (HC) | $6.82\pm 0.365$ | $6.98\pm 0.42$ | $6.47\pm 0.34$ | $5.72\pm 0.076$
Table 2: Transferability across datasets (Pearson’s correlation between VNN outputs and ground truth). Training Testing | FTDC100 (HC) | FTDC300 (HC) | FTDC500 (HC) | OASIS-3 (HC)
---|---|---|---|---
FTDC100 (HC) | $0.49\pm 0.017$ | $0.47\pm 0.018$ | $0.468\pm 0.018$ | $0.387\pm 0.021$
FTDC300 (HC) | $0.498\pm 0.05$ | $0.49\pm 0.042$ | $0.486\pm 0.04$ | $0.374\pm 0.06$
FTDC500 (HC) | $0.51\pm 0.021$ | $0.509\pm 0.02$ | $0.51\pm 0.021$ | $0.374\pm 0.031$
OASIS-3 (HC) | $0.454\pm 0.006$ | $0.437\pm 0.004$ | $0.432\pm 0.003$ | $0.43\pm 0.012$
Figure 5: VNN regression outputs versus chronological age after transferring
VNNs across FTDC datasets. Panel a displays the results obtained by averaging
the outputs of all 100 VNN models that were trained on FTDC100 dataset and
tested on FTDC datasets curated according to different resolutions of
Schaefer’s atlas. Mean of regression outputs for all 100 VNN outputs versus
the ground truth (chronological age) are displayed when the VNNs process
FTDC100 dataset ($m=100$; marked by an asterisk since VNNs were trained on
this dataset), FTDC300 dataset ($m=300$) and FTDC500 dataset ($m=500$). Panel
b and Panel c displays similar results as in Panel a for the VNNs that were
trained on FTDC300 and FTDC500 datasets, respectively.

We tabulate MAE in Table 1 and Pearson’s correlation between ground truth and
VNN output in Table 2. For both tables, the row ID provides the dataset on
which VNN models were trained and the column ID indicates the dataset for
which the VNN performance is reported (after transferring the VNNs if training
and testing datasets are different). For instance, the element with row ID
“FTDC100” and column ID “FTDC300” in Table 1 represents the mean and standard
deviation of MAE evaluated on FTDC300 dataset ($m=300$) for the 100 nominal
VNN models trained on FTDC100 dataset ($m=100$). The elements with same row ID
and column ID in Table 1 and Table 2 provide the baseline performance to gauge
performance after transferring VNNs.

The results in Table 1 and Table 2 show that the performance of VNNs in terms
of MAE and correlation between VNN output and ground truth was preserved after
transferring VNNs across FTDC datasets that were curated according to
different resolutions of Schaefer’s atlas. The transferability of VNNs across
FTDC datasets was corroborated by Fig. 5, which illustrates the plots of the
means of the outputs of $100$ VNN models that were trained on FTDC100 dataset
(or FTDC300 or FTDC500 datasets) and used to process FTDC100, FTDC300, and
FTDC500 datasets. We also remark that this experiment is not feasible for PCA-
regression models as the principal components and the regression model from
one dataset cannot be naively transferred to process another dataset that has
a different number of features. Hence, the observations on FTDC datasets in
Table 1, Table 2, and Fig. 5 validate our theoretical results regarding
transferability of VNNs. The plot for VNN output versus chronological age for
the HC group in OASIS-3 dataset is included in Fig. 17 in Appendix I.

Moreover, Table 1 and Table 2 also report the performance of VNN models
trained on FTDC datasets after transferring to OASIS-3 dataset and vice-versa.
As expected from the observations in Fig. 13b, the VNN performances degraded
significantly when processing the datasets curated according to DKT atlas.
However, the VNN models trained on FTDC100, FTDC300 or FTDC500 datasets
retained significant correlation between predicted age and true age for
OASIS-3 dataset (Table 2) and achieved an MAE of $<8$ years (Table 1). Similar
observations were true for VNN models trained on OASIS-3 dataset when
transferred to FTDC datasets. These observations indicated that the outputs of
VNN models trained on FTDC datasets may not be scaled appropriately to
transfer to OASIS-3 dataset or vice-versa and thus, a mapping or a
transformation between Schaefer’s atlas and DKT atlas is necessary to transfer
regression models optimally. Nevertheless, our experiments imply that the VNN
models trained on FTDC datasets (or OASIS-3 dataset) possessed at least some
information about chronological age that could be transferred when processing
the OASIS-3 datasets (or FTDC datasets). This aspect is explored further in
our results in Section 5.2. Although not the focus of our paper, the
observations above provide some context into why fine-tuning pre-trained deep
learning models for a specific application may be useful for statistical
inference from neuroimaging data [7, 95, 96].

###  5.2 VNN regression model outputs for HC group in OASIS-3 are correlated
with the first eigenvector of $\hat{\mathbf{C}}_{148}$

Our results in Lemma 1 suggested that VNN based statistical inference draws
conceptual similarities with PCA-driven analysis. Hence, we further
investigated whether the regression performance by VNNs in Table 1 could be
characterized by contributions of the principal components of the anatomical
covariance matrix. We focus our discussion here on the VNNs trained on OASIS-3
dataset and VNNs transferred from FTDC100 dataset to OASIS-3 dataset (among
the FTDC datasets, FTDC100 dataset had the largest ratio between the number of
samples and number of features). However, the observations made by experiments
on VNNs trained on FTDC300 or FTDC500 datasets were similar to that for those
trained on FTDC100 and are included in Appendix J.

Recall that the final regression output by VNNs is formed by an unweighted
average function as a readout function. Thus, we can equivalently represent
the functionality of the readout as a simple aggregation of the contributions
of different features or brain regions to the final estimate formed by the VNN
(see (30)). Hence, for every individual, we evaluated the mean of the inner
products (also equivalently referred to as dot product between vectors)
between the vectors of contributions of every brain region with the principal
component of the covariance matrix $\hat{\mathbf{C}}_{148}$ for all $100$ VNN
models. Note that a vector of regional contributions was of the same length as
the number of cortical thickness features (i.e., 148 for OASIS-3) and
therefore, each element of this vector was associated with a distinct brain
region. We use the notation ${\mathbf{p}}_{\sf HC}$ to represent the
population of vectors obtained from the HC group in OASIS-3. To evaluate the
inner product, we used $\bar{\mathbf{p}}_{\sf HC}$, which was obtained from
${\mathbf{p}}_{\sf HC}$ after normalization (norm = 1). We denote the
population of inner products across the HC group in OASIS-3 by
$|\\!<\\!\bar{\mathbf{p}}_{\sf HC},{\mathbf{v}}_{i}\\!>\\!|$ for an
eigenvector ${\mathbf{v}}_{i}$ of $\hat{\mathbf{C}}_{148}$. Note that since
all vectors in $\bar{\mathbf{p}}_{\sf HC}$ were normalized to have norm $1$
and the eigenvectors of $\hat{\mathbf{C}}_{148}$ were of length $1$ by
default, $|\\!<\\!\bar{\mathbf{p}}_{\sf HC},{\mathbf{v}}_{i}\\!>\\!|$
represented the population of the cosine of the angles between the vectors in
${\mathbf{p}}_{\sf HC}$ and eigenvectors ${\mathbf{v}}_{i}$ of
$\hat{\mathbf{C}}_{148}$ across the HC group in OASIS-3.

Figure 6a plots the mean of the inner products observed across the HC group
for the first $30$ eigenvectors of $\hat{\mathbf{C}}_{148}$ for VNNs that were
trained on OASIS-3 dataset. The alignment between the first eigenvector of
$\hat{\mathbf{C}}_{148}$ and the vectors of regional contributions to the VNN
output was significantly stronger as compared to other eigenvectors ($0.991\pm
0.0003$ across the HC group) with relatively smaller associations for
eigenvectors ${\mathbf{v}}_{2}$ ($0.04\pm 0.0029$), ${\mathbf{v}}_{3}$
($0.0746\pm 0.0023$), and ${\mathbf{v}}_{4}$ ($0.0716\pm 0.0026$). Thus, it
could be concluded that the VNNs primarily leveraged the information along the
first principal component of $\hat{\mathbf{C}}_{148}$ (i.e.,
${\mathbf{v}}_{1}$) to achieve the performance in Table 1 and Table 2. Figure
6b illustrates the projection of ${\mathbf{v}}_{1}$ on a brain template. In
Fig. 6b, ${\mathbf{v}}_{1}$ predominantly included bilateral anatomic brain
regions in the parahippocampal gyrus, precuneus, inferior medial temporal
gyrus, and precentral gyrus.

Next, we discuss the correlations between the eigenvectors of
$\hat{\mathbf{C}}_{148}$ and the regional contributions derived from VNNs that
were transferred from FTDC100 to OASIS-3. Recall from Table 1 and Table 2 that
the VNNs that were transferred from FTDC100 to OASIS-3 had a diminished
transference of performance in the task of chronological age prediction on the
OASIS-3 dataset. However, similar to Fig. 6a, the first principal component
${\mathbf{v}}_{1}$ had the largest correlation with the regional contributions
derived from the VNNs in this setting ($0.985\pm 0.0015$) and relatively
smaller associations were observed for ${\mathbf{v}}_{3}$ ($0.072\pm 0.004$)
and ${\mathbf{v}}_{4}$ ($0.111\pm 0.005$). Figure 6c displays the mean of the
inner products observed across the HC group for the first $30$ eigenvectors of
$\hat{\mathbf{C}}_{148}$ obtained by the VNNs that were transferred from
FTDC100 to OASIS-3 dataset. Thus, despite diminished quantitative
transferability between FTDC and OASIS-3 datasets (Table 1 and Table 2), the
VNNs that were transferred from the FTDC100 dataset to OASIS-3 dataset
leveraged the first eigenvector of $\hat{\mathbf{C}}_{148}$ to form the
chronological age predictions for HC group in OASIS-3 dataset. This
observation indicated that the VNNs trained to predict chronological age for
one dataset could gain the ability to leverage the relevant eigenvectors of
the anatomical covariance matrix when transferred to process another dataset
irrespective of the transferability of performance over the task of
chronological age prediction.

Figure 6: Inner product between the normalized vector of regional
contributions to the VNN outputs ($\bar{\mathbf{p}}_{\sf HC}$) and
eigenvectors of $\hat{\mathbf{C}}_{148}$ (anatomical covariance matrix for HC
group in OASIS-3). 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 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 respective entries 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. The first eigenvector (${\mathbf{v}}_{1}$) had the largest
association with $\bar{\mathbf{p}}_{\sf HC}$ ($0.991\pm 0.0003$ across HC
group). The eigenvector ${\mathbf{v}}_{1}$ is plotted on a brain template in
Panel b. Panel c displays a bar plot for $|\\!<\\!\bar{\mathbf{p}}_{\sf
HC},{\mathbf{v}}_{i}\\!>\\!|$ for the first $30$ eigenvectors of
$\hat{\mathbf{C}}_{148}$ and the vectors of regional contributions,
$\bar{\mathbf{p}}_{\sf HC}$ in the scenario when VNNs were transferred from
FTDC100 to OASIS-3 dataset.

###  5.3 Brain Age Prediction using VNNs in OASIS-3

Next, we performed the regional analyses of residuals obtained from VNN models
for individuals in HC and AD+ groups according to the procedure discussed in
Section 4.4.1. We remark that several existing studies on brain age prediction
have utilized models that achieved better MAE on the HC group than our results
in Table 1 and Table 2 [15, 16, 17, 73]. However, our primary focus here is on
demonstrating that the outputs of the final layers of VNNs trained on data
from HC group enabled us to associate contributions to brain regions that were
characteristic of neurodegeneration in AD. Thus, this section summarizes the
insights provided by VNNs that have neither been explored nor feasible for
most existing brain age evaluation frameworks based on black box
implementation of deep learning models.

To start with, we leveraged the VNN models trained on data from HC group in
OASIS-3 in Section 5.1 to process a combined cortical thickness dataset from
HC and AD+ groups. Also, in this scenario, the covariance matrix was formed
from the combined dataset ($\hat{\mathbf{C}}_{148}^{\sf AD+}$ from cortical
thickness data from HC and AD+ groups). This choice was made in order to
incorporate the neurodegeneration-driven changes in the anatomical covariance
into the model. Next, for every VNN model, we evaluated the distribution of
regional residuals in AD+ and HC groups as described in Section 4.4.1. An
example of the outputs of the regional analysis of a VNN model is included in
Appendix H.

We note that the regional residuals derived from VNNs are independent of the
MAE performance of VNN models in predicting chronological age for the
individuals in HC group. Therefore, besides analyzing brain age in OASIS-3
dataset, we also performed exploratory regional analyses to assess the
interpretability offered by VNN models that were trained on FTDC datasets and
transferred to process OASIS-3 dataset. We aimed to discern the
interpretability offered by VNNs by analyzing the group differences in
regional residuals between AD+ and HC groups and the association of regional
residuals with CDR sum of boxes in the AD+ group.

####  5.3.1 Regional analyses of outputs of VNNs reveal brain regions
characteristic of Alzheimer’s disease

Our results for the OASIS-3 dataset showed that the brain regions
characteristic of AD had significant differences in regional residual
distributions (defined in (33)) for AD+ group as compared to HC group. Note
that the results presented in this section were checked for robustness using
the family of 100 nominal models trained to predict chronological age of
individuals in the HC group in OASIS-3 (described in Section 5.1).

Figure 5.3.1a illustrates the significant group differences in cortical
thickness for AD+ and HC groups (t-test, t-statistic for regions with
Bonferroni corrected $p$-value $<0.05$ are projected on the brain template).
Since AD+ group is characterized by accelerated loss of cortical thickness in
various brain regions, we expected the outcomes of the analysis of regional
residuals to have significant overlap with cortical atrophy patterns in Fig.
5.3.1a. Figure 5.3.1b displays the robustness (determined via analyses of 100
VNN models) for various brain regions in having an elevated regional effect in
their corresponding residual elements for AD+ group with respect to HC group.
The significance of brain regions in each VNN model was determined via
comparisons between HC and AD+ groups and from VNN final layer outputs, where
the parameters of the VNNs were determined via training to predict
chronological age for the HC group in OASIS-3. The most significant regions
with elevated regional residuals in AD+ with respect to HC were concentrated
in bilateral inferior parietal, temporal, entorhinal, parahippocampal,
precuneus, subcallosal, and precentral regions. All these brain regions,
except for precentral and subcallosal, mirrored the cortical atrophy in Fig.
5.3.1a, and these regions are known to be highly interconnected with
hippocampus [97].

Our findings here are also consistent with the existing studies that study
anatomic changes in AD. Inferior parietal and medial temporal regions are well
known regions of atrophy for the typical variant of Alzheimer’s disease [98,
99]. Cortical atrophy in these regions is also accompanied by tau protein
deposition [100]. Precuneus and parahippocampal regions are implicated in the
early stages of AD [101, 102] and precuneus is also a region of interest in
individuals with early-onset AD [103]. The aforementioned brain regions are
also assigned interpretability in the study in [15].

Although the results in Fig. 5.3.1b provided a meaningful regional profile for
AD+ group, we further performed exploratory analysis to check whether the
regional residuals had any clinical significance. To this end, we evaluated
the correlations between CDR sum of boxes and the regional residuals derived
from final layer VNN outputs for the AD+ group for all 100 VNN models. Figure
5.3.1c plots the means of the Pearson’s correlations observed between regional
residuals and CDR sum of boxes across the 100 VNN models for AD+ group on the
brain template. Interestingly, the brain regions with the most significant
correlations aligned with the most robust brain regions in the regional
profile corresponding to elevated regional residuals in the AD+ group in Fig.
5.3.1b. We clarify that our analysis in this context was not restricted to the
highlighted brain regions in Fig. 5.3.1b. Thus, the regional profile extracted
from VNN final layer outputs in Fig. 5.3.1b also captured information about
the dementia severity. Next, we checked if the regional findings in Fig.
5.3.1b and Fig. 5.3.1c translated into elevation in $\Delta$-Age and the
correlation between $\Delta$-Age and CDR sum of boxes for AD+ group.

(preceding page): Interpretable $\Delta$-Age evaluation in OASIS-3 dataset.
Panel a displays the results of group differences in cortical thickness
between AD+ and HC groups. Regions with significant differences (two-sided
t-test, Bonferroni corrected $p$-value $<0.05$) are identified and the
corresponding $t$-statistics are projected on a brain template. Negative
$t$-statistic for a brain region suggests that the AD+ group had significant
cortical atrophy in that region as compared to HC group. Panel b displays the
robustness of the significantly elevated regional residuals for AD+ group with
respect to HC group for different brain regions. For every VNN model in the
set of $100$ nominal VNN models that were trained on HC group in OASIS-3, the
analyses of regional residuals helped isolate brain regions that corresponded
to significantly elevated regional residuals in AD+ group with respect to HC
group (ANOVA: Bonferroni corrected $p$-value $<0.05$, ANCOVA with age and sex
as covariates: $p$-value $<0.05$). After performing this experiment for $100$
VNN models, the robustness of the observed significant effects in a brain
region was evaluated by calculating the number of times a brain region was
identified to have significantly elevated regional residuals in AD+ group with
respect to HC group. The number of times a brain region was linked with
significantly elevated regional residuals in AD+ group with respect to HC
group is projected on the brain template and can be perceived as a marker of
its robustness as a contributor to accelerated aging in AD+ group. Panel c
projects the mean Pearson’s correlation between regional residuals (derived
for each VNN model in the set of 100 nominal VNN models) and CDR sum of boxes
for AD+ group on the brain template. Panel d displays the distribution of
$\Delta$-Age for HC and AD+ groups. The gap between the mean of $\Delta$-Age
for AD+ and HC is $2.36$ years and significant (Cohen’s $d$ = 0.797). The
elevated brain age effect here is characterized by regional profile in Panel
b. Panel e displays the scatter plot for CDR sum of boxes versus $\Delta$-Age
in AD+ group. The correlation between $\Delta$-Age and CDR sum of boxes could
be attributed to the correlations observed between regional residuals and CDR
sum of boxes in Panel c.

####  5.3.2 $\Delta$-Age is elevated in AD+ group and correlated with CDR

We evaluated the $\Delta$-Age for HC and AD+ groups in OASIS-3 according to
the procedure specified in Section 4.4.2. We also investigated the Pearson’s
correlation between $\Delta$-Age and CDR sum of boxes scores in AD+ group.
Figure 5.3.1d illustrates the distributions for $\Delta$-Age for HC and AD+
groups ($\Delta$-Age for HC: $0\pm 2.66$ years, $\Delta$-Age for AD+: $2.36\pm
3.22$ years). The difference in $\Delta$-Age for AD+ and HC groups was
significant (Cohen’s $d$ = 0.797, ANCOVA with age and sex as covariates:
$p$-value = $3.98\times 10^{-22}$, partial $\eta^{2}=0.103$). Also, the
$p$-value associated with age in ANCOVA was $0.61$ and hence, group difference
in $\Delta$-Age was not driven by the difference in the distributions of
chronological age for HC and AD+ groups. Figure 5.3.1e plots the scatter plot
between CDR sum of boxes score and $\Delta$-Age for the AD+ group. The
Pearson’s correlation between CDR sum of boxes score and $\Delta$-Age was
$0.375$ ($p$-value $=2.72\times 10^{-8}$), thus, implying that the
$\Delta$-Age evaluated for AD+ group captured information about dementia
severity. Hence, as expected, the $\Delta$-Age for AD+ was likely to be larger
with an increase in CDR sum of boxes scores. For instance, the mean of
$\Delta$-Age for individuals with CDR sum of boxes greater than 4 was $3.76$
years, and for those with CDR sum of boxes smaller than or equal to 4 was
$1.79$ years.

Given that the age-bias correction procedure is a linear transformation of VNN
outputs, it can readily be concluded that the statistical patterns for
regional residuals in Fig. 5.3.1b and Fig. 5.3.1c lead to elevated
$\Delta$-Age in Fig. 5.3.1d and correlation between $\Delta$-Age and CDR sum
of boxes scores in Fig. 5.3.1e. Thus, a major implication of the findings in
Fig. 5.3.1 is that certain regions contribute abnormally in the context of AD-
related neurodegeneration as compared to HC group to the age estimate formed
by the VNN models and therefore, the analysis of regional residuals derived
from the outputs at the final layer of a VNN model provides a feasible way to
characterize accelerated biological aging in AD+ group with a regional
profile. Additional figures and details pertaining to VNN outputs and brain
age before and after the age-bias correction was applied are included in
Appendix I and they demonstrate the distributions of $\Delta$-Age and brain
age across HC and AD+ groups in OASIS-3.

We have also previously reported in [2] that the gap between $\Delta$-Age for
AD, MCI and HC groups was preserved across different resolutions of Schaefer’s
atlas after transferring VNN models. The results in our previous work in [2]
were a direct consequence of the transferability of VNNs across different
scales of Schaefer’s atlas. For completeness, we further demonstrate that
there was spatial consistency between regional profiles for AD that were
evaluated across different resolutions of Schaefer’s atlas for a distinct
dataset in Appendix K.

###  5.4 Robustness of regional profiles of elevated $\Delta$-Age in AD+ group
hinges on the correlations between regional residuals and principal components
of anatomical covariance matrix $\hat{\mathbf{C}}_{148}^{\sf AD+}$

Thus far, we have independently shown that VNNs trained as regression models
are transferable across multi-scale datasets without any re-training and VNNs
provide an interpretable perspective to accelerated brain age in
neurodegeneration. However, despite lack of quantitative transferability
between FTDC and OASIS-3 datasets, our results in Fig. 6 showed that the VNNs
that were transferred from FTDC100 to OASIS-3 and the VNNs trained on OASIS-3
exploited similar eigenvectors of the anatomical covariance matrix to form
predictions for chronological age for HC group. Here, we explore whether the
regional interpretability provided by VNNs for $\Delta$-Age depends on the
principal components of the anatomical covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$.

####  5.4.1 Regional residuals derived from VNNs trained on OASIS-3 are
correlated with principal components of anatomical covariance matrix

We started by investigating the relationship between regional residuals
derived from VNNs trained on OASIS-3 and the principal components of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$. For this purpose, we evaluated the inner
product of normalized residual vectors (norm = 1) obtained from VNNs trained
on OASIS-3 and the eigenvectors of the covariance matrix
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ for the individuals in AD+ group to
determine whether any specific eigenvectors (principal components) of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ were instrumental to recover the findings
in Fig. 5.3.1b. The normalized residual vector is denoted by
$\bar{\mathbf{r}}_{\sf AD+}$. For every individual, the mean of the absolute
value of the inner product $|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{i}\\!>\\!|$ (where ${\mathbf{v}}_{i}$ is the $i$-th
eigenvector of $\hat{\mathbf{C}}_{148}^{\sf AD+}$) was evaluated for the 100
VNN models that were used to derive the results in Fig. 5.3.1. Figure 8a plots
the mean inner product for eigenvectors associated with $50$ largest
eigenvalues of $\hat{\mathbf{C}}_{148}^{\sf AD+}$. The three largest mean
correlations with the regional residuals in AD+ group were observed for the
fourth principal component of $\hat{\mathbf{C}}_{148}^{\sf AD+}$
(${|\\!<\\!\bar{\mathbf{r}}_{\sf AD+},{\mathbf{v}}_{4}\\!>\\!|=0.6\pm
0.006}$), third principal component ($|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{3}\\!>\\!|=0.384\pm 0.007$), and the first principal
component ($|\\!<\\!\bar{\mathbf{r}}_{\sf
AD+},{\mathbf{v}}_{1}\\!>\\!|=0.294\pm 0.0003$). These principal components
are plotted on a brain template in Fig. 8b. Interestingly, the largest
association was observed for the fourth eigenvector of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$, which spanned the subcallosal and medial
frontal regions in the right hemispheres. The third eigenvector spanned the
subcallosal region in the right hemisphere and temporal pole and
parahippocampal region in both hemispheres. The first principal component of
$\hat{\mathbf{C}}_{148}^{\sf AD+}$ spanned brain regions similar to that by
the first principal component of $\hat{\mathbf{C}}_{148}$ (Fig. 6b). These
observations suggested that the observed regional effects in subcallosal,
parahippocampal and temporal pole regions were the most significant AD-related
contributors to the observed elevated $\Delta$-Age effect in Fig. 5.3.1. These
brain regions were also implicated as contributors to elevated brain age in AD
in an independent dataset in Appendix K.

Figure 8: Mean inner product between the normalized vector of regional
residuals (norm = 1) of VNN outputs (VNNs trained on OASIS-3) obtained from
AD+ group and the eigenvectors of $\hat{\mathbf{C}}_{148}^{\sf AD+}$
(covariance matrix of combined HC and AD+ group), and eigenvector plots. Panel
a illustrates a bar plot for $|\\!<\\!\bar{\mathbf{r}}_{\sf