# Learnable Spectral Wavelets on Dynamic Graphs to Capture Global Interactions

Anson Bastos,1 Abhishek Nadgeri,2 Kuldeep Singh,3 Toyotaro Suzumura,4 Manish
Singh1

###### Abstract

Learning on evolving(dynamic) graphs has caught the attention of researchers
as static methods exhibit limited performance in this setting. The existing
methods for dynamic graphs learn spatial features by local neighborhood
aggregation, which essentially only captures the low pass signals and local
interactions. In this work, we go beyond current approaches to incorporate
global features for effectively learning representations of a dynamically
evolving graph. We propose to do so by capturing the spectrum of the dynamic
graph. Since static methods to learn the graph spectrum would not consider the
history of the evolution of the spectrum as the graph evolves with time, we
propose a novel approach to learn the graph wavelets to capture this evolving
spectra. Further, we propose a framework that integrates the dynamically
captured spectra in the form of these learnable wavelets into spatial features
for incorporating local and global interactions. Experiments on eight standard
datasets show that our method significantly outperforms related methods on
various tasks for dynamic graphs.

## Introduction

Recently there has been tremendous progress in the domain of Graph
Representation Learning (Khoshraftar and An 2022). The aim here is to develop
novel methods to learn the features of graphs in a vector space. Such
approaches have successful applications in the domain of image recognition
(Han et al. 2022), computational chemistry (Ying et al. 2021), drug discovery
(Stärk et al. 2022) and Natural Language processing (Wu et al. 2021a).
Although effective, the underlying graphs are static in nature.

In many real-world scenarios, graphs are dynamic, for example, social
networks, citation graphs, bank transactions, etc. For such cases, various
approaches have been developed (see survey by Kazemi et al. (2020)). Broadly,
these methods aim to learn the evolving nature of graphs through spatial
features relying on local neighborhood aggregation (local dependencies)
(Pareja et al. 2020; Goyal, Chhetri, and Canedo 2020). For example,
researchers (Pareja et al. 2020; Shi et al. 2021) have resorted to using GNNs
along with RNNs to capture the dynamic evolving nature of graphs. With initial
successes, these methods are inherently limited to capturing local
interactions, missing out on other important information. For example, in the
case of a dynamic money transaction graph, having a skewed representation of
fraud users among genuine users, there would exist links between these fraud
and genuine users, thus, giving rise to high-frequency signals. In the case of
local neighborhood aggregation (attending to low-frequency signals), the
majority node(genuine) pattern will cause the fraudulent pattern to diminish.
Thus, capturing global properties becomes necessary. Here, the global
information would help identify the fraud pattern eventually assisting in
identifying the criminal. Similarly, in the citation graph, the local
properties will help to understand the category of paper. In contrast, global
properties will help to understand the amount of interdisciplinary research
across research domains. Hence, learning the global dependencies is crucial
for dynamic graphs which is a relatively unexplored area in its scientific
literature.

Figure 1: For Brain dataset (Xu et al. 2019b) at a given timestep, the figure
(a) shows the learned filter functions at different scales by our proposed
framework resembling band reject, low pass and all pass filters. The three
diagrams to the right show the corresponding wavelets for a node (colored in
yellow in figures (b,c,d) at lowest tip of the Brain) in the brain stem
region. The graph nodes indicate the Regions of Interest(ROI) in brain. As
scales change, the ROIs that the wavelet focuses on change from
concentrated(local interactions) to diffused(global interactions). Moreover,
the wavelets respect the brain structure and focus on the concerned
region(brain stem in this case) thus mitigating noise due to interactions from
unrelated regions.

In this paper we aim to encompass global dependencies in dynamic graphs,
beyond local neighborhood aggregation, and view it through the lens of
spectral graph theory (Hammond, Vandergheynst, and Gribonval 2019). For the
same, we introduce a novel concept of learnable spectral graph wavelets to
capture global interactions of dynamic graph. The concept of learnable wavelet
has following motivation: (i) Computing the spectra using full eigen
decomposition is computationally expensive and requires $\mathcal{O}(N^{3})$
time. Moreover, this gets even more computationally costly by a factor of the
number of timesteps($T$) considered for dynamic graphs that evolve over every
timestep. Spectral graph wavelets can be computed efficiently in
$\mathcal{O}(N+|\mathcal{E}|)$ as we shall see in the following sections. (ii)
Wavelets are sparser as compared to the graph Fourier transform adding to the
computational benefits (Tremblay and Borgnat 2014). (iii) Wavelets give a
sense of localization in the vertex domain of the graphs thus enabling
interpretability of the convolutions while also being able to capture the
global properties of the graph by changing the scale parameter (c.f., Figure
1). iv) Focusing on dynamic learnable wavelets helps in capturing global
properties in that evolving spectra as the graph changes with time, where
static and non-learnable wavelet methods (Xu et al. 2019a) show empirical
limitations.

Furthermore, build on the recent success of capturing neighborhood features
for evolving graphs, we propose to learn homogeneous representation of spatial
and spectral features to apprehend both local and global dependencies. Our
approach is very intuitive as it keeps the proven local information intact
whilst adding the global properties through learnable graph wavelet approach.
Similar approaches have been proven in computer vision (Srinivas et al. 2021)
and NLP (Prasad and Kan 2019) where such restrictive inductive bias work well
on learning local properties but miss out on global interactions. Also, the
local aggregation leads to problems such as over smoothing and may not work
well on heterophilic graphs (Wu et al. 2021b). Thus, using our methods will
help to alleviate the drawbacks of the existing popular methods on dynamic
graphs.

Our key contributions are two-fold: 1) a novel approach to learning spectral
wavelets on dynamic graphs for capturing global dependencies (with its
theoretical foundations), 2) a novel framework named DEFT that combines
spectral features obtained using learnable wavelets into spatial feature of
the evolving graphs. For effective use in downstream tasks, DEFT integrates
the spatial and spectral features into homogeneous representations which
allows capturing shift invariance (Oppenheim and Schafer 1975) among the node
features that could arise from the temporal nature.

## Related Works

There has been considerable work on static graphs from the spatial and
spectral domain perspective. Some works such as (Kipf and Welling 2017;
Hamilton, Ying, and Leskovec 2017; Veličković et al. 2018) focus on
effectively learning spatial properties. Similarly, efforts such as (Levie et
al. 2018; Balcilar et al. 2020; Bastos et al. 2022) have looked at graphs from
a spectral perspective. GWNN (Xu et al. 2019a) has proposed to use spectral
graph wavelets on static graphs. However, it obtains static wavelets using
heat kernel-based filters. Unlike GWNN, we learn spectral graph wavelets for
dynamically evolving graphs. We observe in the experiment section that
learnable wavelets perform significantly better for dynamic graphs compared to
static wavelet methods.

One straight-forward way to use methods developed for static graphs on
evolving dynamic graphs is to use RNN modules such as GRU and LSTM in addition
to GNN modules to capture the evolving graph dynamism. This idea has been
explored in works such as (Seo et al. 2016; Narayan and Roe 2018; Manessi,
Rozza, and Manzo 2020). However, these models suffer a performance drop if new
nodes are introduced in the graph; the GNNs may not be able to adapt to the
evolving graph structure. Thus, EvolveGCN (Pareja et al. 2020) has introduced
to use RNNs to learn the parameters of the evolving GNN. The GNN, a GCN (Kipf
and Welling 2017) in this case, is thus used ahead of the RNN module that
captures the graph dynamism and offers promising results on dynamic graphs.
EvolveGCN has limitation that it generalizes to unseen nodes in future
timesteps, which is not always the case in real-world scenarios. Xu et al.
(2019c) proposed similar approach by learning parameters of a GAT with an RNN
that focus on graph topology discrepancies.

Furthermore, autoencoder-based methods have been introduced, such as in (Goyal
et al. 2017; Goyal, Chhetri, and Canedo 2020; Xu, Singh, and Karniadakis 2022)
that focus on reconstructing the graphs in future timesteps as an objective.
Xiang, Huang, and Wang (2022) proposed a meta-learning framework in which the
objective is to predict the graph at future timestep. However, these works
focus only on learning spatial properties to capture local dependencies and
ignores the global dependencies that may emerge due to dynamic nature of the
graph.

## Preliminaries

Consider a graph with vertices and edges as $(\mathcal{V},\mathcal{E})$ and
adjacency matrix $A$. The laplacian($L=D-A$) can be decomposed into its
orthogonal basis, namely the eigenvectors($U$) and eigenvalues($\Lambda$)
as:$L=U\Lambda U^{*}$. Let $X\in R^{N\times d}$ be the signal on the nodes of
the graph. The Fourier Transform $\hat{X}$ of $X$ is then given as:
$\hat{X}=U^{*}X$.

Spectral graph wavelet transforms (Hammond, Vandergheynst, and Gribonval 2011)
are obtained by functions of the laplacian $L$. Consider a linear self-adjoint
operator($g(L)$) that acts on each component in the spectral space. We define
a parameter $s$ for the scaling in the spectral space. The spectral graph
wavelet at any given vertex $n$ is defined as the impulse($\delta_{n}$)
response of the wavelet operator at vertex $n$:

$\psi_{s,n}(m)=\sum_{k=1}^{N}g(s\lambda_{k})U_{k}^{*}(n)U_{k}(m)$

The n-th wavelet coefficients at scale $s$ represented as $W_{f}(s,n)$ can be
obtained by taking the inner product of the function $f$ in the vertex domain
with these wavelets as

$W_{f}(s,n)=\left<\psi_{s,n},f\right>=\sum_{k=1}^{N}g(s\lambda_{k})\hat{f}(k)U_{k}(m)$
(1)

In our work, we propose to learn the wavelet coefficients for dynamic graphs,
where the exact form of the scaling function $g(s\lambda_{k})$ is
parameterized.

Figure 2: Figure (a) shows the DEFT architecture. The yellow outer box shows
the modules that learn the weight parameters of the GNNs in an evolving
manner. These parameters are then given to the spectral and spatial modules to
learn the corresponding features. The learned features then propagates to the
homogeneous representation module (HRM) followed by the aggregation module(AM)
which together forms the integration module (red outer box) of our framework.
Figure (b) explains the components of the spectral module. Similar to image
(a), the yellow outer box learns the filter functions for the wavelet
approximators.

## Method

Figure 2 illustrates the proposed DEFT (Dynamic wavElets For global
inTeractions) framework that comprises the following modules: 1) Spectral
component: focuses on global features of the graph in the form of learnable
spectral graph wavelets. 2) Spatial component: the necessity of this component
is to mitigate the unsmooth spectrum phenomenon i.e. the node features gets
correlated with eigenvectors corresponding to high frequencies causing
information loss (Yang et al. 2022). To resolve this, spatial module focus on
the local neighborhood of the graph in addition to that captured by the
wavelets. 3) Integration module: Finally, for a homogeneous representation of
the global properties captured by the spectral component and the local
properties learned by the spatial component, we propose a module that learns
invariant representations and use these in an aggregation module for
downstream tasks.

### Spectral Module

We aim to capture global features without requiring the full eigen
decomposition of the laplacian. Here, we propose to approximate the wavelet
operator using some function. We use the Chebyshev polynomials to be
consistent with the literature (Hammond, Vandergheynst, and Gribonval 2011).
It is known that for any function $h$ with domain $y\in[-1,1]$ there exists a
convergent Chebyshev series:

$h(y)=\sum_{k=0}^{\infty}c_{k}T_{k}(x)$

with the Chebyshev polynomials satisfying the recurrence relation
$T_{k}(y)=2yT_{k-1}(y)-T_{k-2}(y),T_{0}=1,T_{1}=y$ and the coefficients
$c_{k}$ are given by:

$c_{k}=\frac{2}{\pi}\int_{-1}^{1}\frac{T_{k}(y)h(y)}{\sqrt{1-y^{2}}}dy=\frac{2}{\pi}\int_{0}^{\pi}cos(k\theta)h(cos(\theta))d\theta$

In order to approximate the function $g(s_{j}x)$ for each scale $j$, we need
to bring the domain $x$ in $[-1,1]$. Noting that $0\leq x\leq\lambda_{max}$
for the laplacian $L$, we perform the transform
$y=\frac{x-a}{a},a=\frac{\lambda_{max}}{2}$. We now define
$\overline{T_{k}}(x)=T_{k}(y)=T_{k}(\frac{x-a}{a})$ and the approximation for
$g$ looks as below

$g(s_{j}x)=\sum_{k=0}^{\infty}c_{j,k}\overline{T_{k}}(x)$ (2)

with the coefficients given by,

$c_{j,k}=\frac{2}{\pi}\int_{0}^{\pi}cos(k\theta)g(s_{j}(a(cos(\theta)+1)))d\theta$

We truncate the polynomial to $M$ terms, which is the filter order. The
coefficients $c_{j,k}$ which are analytical coefficients of the filter
function as desired are approximated using functions parameterized by GNNs and
MLPs, as we shall see next. $f_{c}^{s}$ is the parameterized form of it,
obtained in spectral module at scale s. A GNN is used to perform message
passing over the input graph along with the node features $v$ at layer $l$ for
neighborhood $\mathcal{N}$.

$v_{im}^{l}=A_{f}({v_{j}^{l}|v_{j}\in\mathcal{N}(v_{i})}),\ \ \ \
v_{i}^{l+1}=U_{f}(v_{i}^{l},v_{im}^{l})$

Here, $A_{f},U_{f}$ are the aggregation and update functions, respectively.
The update function could contain a non-linearity such as leaky ReLU for
better expressivity of the learned function. Since we intend to learn the
filter coefficients $f_{c}\in R^{M}$ for the concerned graph($G$), we apply a
pooling layer to get an intermediate vector representation($v_{G}\in
R^{d_{1}}$) from the output of the GNN. The pooling layer converts a set of
vectors(one for each node of the graph) to a single vector representation. For
the final filter coefficients $f_{c}$, we apply a two-layer MLP with
activation($\sigma$) to $v_{g}$

$f_{c}=W_{2}\ \ \sigma\left(W_{1}v_{g}\right)\\\ $ (3)

where $W_{1}\in R^{d_{2}\times d_{1}},W_{2}\in R^{d_{2}\times M}$ are
learnable weights. Since the two-layer MLP is a universal approximator (Hornik
1991) we can be assured of the existence of a function in this space that
learns the desired mapping to the filter coefficients. In principle, any
message passing GNN can be used to perform the update and aggregation steps.
This process can be repeated with multiple GNNs for learning multiple filter
functions. As we consider dynamic graphs, we would like to evolve the
parameters of the GNN with time (Dynamic Parameter learning Module of Figure
2, common for both spectral and spatial components). Inspired from (Pareja et
al. 2020), we use an RNN module for generating the parameters for the GNN in
layer $l$ at time $t$:

$W_{t}^{l}=RNN(H_{t}^{l},W_{t-1}^{l})$ (4)

where $W_{t}^{l}$ and $H_{t}^{l}$ are the hidden state and input at layer $l$
and time $t$ of the RNN. In the below pseudo code, we outline our method to
evolve the Spectral Module($ESpectral$) for dynamically learning filter
coefficients per timestep

1:function $f_{c_{t}}=ESpectral$($A_{t},H_{t}^{(l)},W_{t-1}^{(l)}$)

2: $W_{t}^{(l)}=RNN(W_{t-1}^{(l)})$

3: $H_{t}^{(l+1)}=GNN(A_{t},H_{t}^{(l)},W_{t}^{(l)})$

4: $v_{g_{t}}=Pool(A_{t},H_{t}^{(l+1)})$

5: $f_{c_{t}}=W_{2}\ \ \sigma\left(W_{1}v_{g_{t}}\right)$

6:end function

After learning GNN parameters, we need to learn filter coefficients for
evolving graph. Learning the filter coefficients $f_{c}^{s}$ at a given
timestep, we can obtain the wavelet operator $g(L)$ at scale $s=1$ using
equation 2. For operators at a given scale $s>0$, we could learn different
parameterizations of the filter function at every scale. Note that
approximating the functions at different scales in this manner would incur a
storage cost of $\mathcal{O}(JN)$ for the filter coefficients. Along with this
the storage and computation complexity would increase $J$ times for
dynamically computing the filter coefficients from the GNNs. Thus we propose a
“rescaling trick” wherein to obtain the operator at any scale $s>0$, we
perform the appropriate change of variables in equation 2 to get $g(sL)$,
keeping the coefficients $c_{j,k}$ fixed. That is for a scale $s$, instead of
mapping $\lambda\xrightarrow[]{}g(\lambda)$ it would be mapped to
$\lambda\xrightarrow[]{}g(s\lambda)$. It helps maintain parameter
efficiency(by a factor of $J$) as the GNN weights(and also filter
coefficients) are tied across all the scales. Note here the exact filter
learned would vary since the maximum frequency is the same. Hence, a bandpass
at scale ”one” may get converted to a highpass at scale ”two”. Once we
approximate $g(\lambda)$, we use it in learning the wavelet coefficients
(output of spectral module) $W_{f}(s,n)$ as described in equation 1.

Now, we give an approximation bound between the learned and desired filter
function under the given framework for dynamic graphs with $N$ nodes and
varying edges and signals. Please see appendix for all proofs.

######  Lemma 1.

Consider $G^{t}(\lambda)$ to be the filter function at time $t$. Assume the
Markov property to hold in the dynamic setting where the desired filter
function($G^{t+1}(\Lambda)$) at time $t+1$ depends on the past state at time
$t$($G^{t}(\Lambda)$). Consider this mapping between the past state and the
current states to be captured by an arbitrary functional $f$ such that
$G^{t+1}(\lambda)=f(G^{t}(\lambda_{1}),G^{t}(\lambda_{2}),\dots
G^{t}(\lambda_{N}),\lambda)$ and we assume $f$ to be $L$ lipschitz continuous.
Further, let $C_{t}=U_{t}G^{t}(\lambda)U_{t}^{T}\in R^{N\times N}$ represents
the convolution support of the desired spectral filter and $C_{t}^{a}$ be the
learnt convolution support at time $t$. Then, we have:

(i) $\left\lVert C_{t+1}^{a}-C_{t+1}\right\rVert_{F}\leq
LN^{2}\sqrt{\left\lVert
C_{t}^{a}-C_{t}\right\rVert_{F}^{2}+\epsilon_{ca}^{2}}+\epsilon_{fa}$

(ii) $\left\lVert C_{t+1}^{a}-C_{t}^{a}\right\rVert_{F}\leq\left\lVert
C_{t+1}-C_{t}\right\rVert_{F}+2\sqrt{N}\epsilon_{ca}$

where $\epsilon_{ca}$ and $\epsilon_{fa}$ are the filter polynomial(Chebyshev)
and function approximation errors and depends on the order of the polynomial,
number of training samples, model size etc.

Above result gives us a relation between the error at times $t$ and $t+1$ and
has a factor of $N^{2}L$. Thus it requires the filter function to be
smooth($L<\frac{1}{N^{2}}$) for convergence because under the given Markov
assumptions, the past errors could accumulate in future timesteps. This is
precisely why we need a gated model like GRU/LSTM in the Dynamic Parameter
Generation module that can decouple the approximation of the filter function
at a given timestep from the error in previous timesteps.

### Spatial Module

Yang et al. (2022) concluded that high-frequency components of the signal on
the graph get strengthened with a high correlation with each other, and the
smooth signals become weak (unsmooth spectrum) i.e., the cosine similarity
between the transformed signal and the low eigenvector reduces with the
layers. In our setting, we illustrate that the factor by which the signal
corresponding to the low-frequency component gets weakened is directly
proportional to the magnitude of the frequency response at that frequency.

######  Lemma 2.

Let $G(\lambda)$ be the frequency response at frequency $\lambda$. Let
$\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{n}$ be the eigenvalues in
descending order and $p_{1},p_{2},\dots p_{n}$ be the corresponding
eigenvectors of the laplacian of the graph. Define $\lambda_{max}$ to be the
eigenvalue at which $G(\lambda)$ is maximum. Let $C^{l}$ represent the
convolution support($UG(\Lambda)U^{T}$) of the spectral filter at layer $l$.
Then the factor by which the cosine similarity between consecutive layers
dampens is
$\underset{l\xrightarrow[]{}\infty}{lim}\frac{|cos(\left<C^{l+1}h,p_{n}\right>)|}{|cos(\left<C^{l}h,p_{n}\right>|}=\frac{G(\lambda_{n})}{G(\lambda_{max})}$.

If $\lambda_{max}$ belongs to one of the higher regions of the spectrum and
$G(\lambda_{n})$ is lower, then as the layers increase, the signals will lose
the low-frequency information. Since the Spectral Module aims to capture the
high-frequency components that may lead to an unsmooth spectrum, to resolve
this, we explicitly strengthen the low-frequency components by using local
neighborhood aggregation. For this, we inherit the message passing GNNs, the
parameters of which are generated using the RNN module inspired from (Pareja
et al. 2020). The below pseudo code outlines the process to evolve the spatial
module($ESpatial$).

1:function
$[H_{t}^{(l+1)},W_{t}^{(l)}]=ESpatial$($A_{t},H_{t}^{(l)},W_{t-1}^{(l)}$)

2: $W_{t}^{(l)}=RNN(W_{t-1}^{(l)})$

3: $H_{t}^{(l+1)}=GNN(A_{t},H_{t}^{(l)},W_{t}^{(l)})$

4:end function

### Integration Module

Homogeneous representation Module(HRM): aims to achieve homogeneous
representations from the spatial and spectral properties along with time
features, which is essential for its usage in downstream tasks. A
straightforward way is by concatenating two representations. However, due to
the dynamic nature of the graph that evolves with time, we propose a learnable
module that provides a notion of distance between the representations and
helps induce position/structure information. Due to graph dynamism, it is
important for the features to satisfy the shift invariance property:

######  Definition 1 (Shift Invariance (Li et al. 2021)).

Any two vectors $v_{i}=f_{1}(v_{i}^{{}^{\prime}})$ and
$v_{j}=f_{1}(v_{j}^{{}^{\prime}})$ satisfy the shift invariance property if
the inner product $\left<v_{i},v_{j}\right>$ is a function of
$v_{i}^{{}^{\prime}}-v_{j}^{{}^{\prime}}$ i.e.
$\left<v_{i},v_{j}\right>=f_{2}(v_{i}^{{}^{\prime}}-v_{j}^{{}^{\prime}})$.

$f_{1}$ is an arbitrary function and $f_{2}$ is a linear transformation of
$f_{1}$. Above property ensures that the relative distance between two nodes
is maintained in the inner product space even if their absolute positions
change(for example with addition of new nodes in the graph with time).
Inspired from Fourier features (Rahimi and Recht 2007, 2008), for a node $i$,
the spectral($v_{gi}\in R^{d_{g}}$) and spatial($v_{li}\in R^{d_{l}}$)
embeddings are concatenated($\|$) along with the timestamp($t\in R^{d_{t}}$)
information if available. Then, it is passed to an MLP.

$v_{gli}=MLP(v_{gi}\|v_{li}\|t)$

In order to obtain Fourier features($v_{ffi}$) from above intermediate
representation, we take the element-wise sine,cosine and concatenate the two
as:

$\displaystyle v_{ffi}=sin(v_{gli})\|cos(v_{gli})$

######  Property 1.

The vector $v$ obtained by concatenation of the element wise sine and cosine
of another vector $v^{{}^{\prime}}$ i.e.
$v=(sin(v^{{}^{\prime}})\|cos(v^{{}^{\prime}}))$, satisfies the shift
invariance property.

The property 1 can be readily noted by observing that
$cos(a-b)=cos(a)cos(b)+sin(a)sin(b)$. Taking the inner
product($\left<.\right>$) of the above features for two nodes($i,j$) gives

$\left<v_{ffi},v_{ffj}\right>=\sum
cos(W_{r}(v_{gli}-v_{glj}))=ff_{W_{r}}(v_{gli}-v_{glj})$

Property 1 is beneficial if these representations are used in an attention-
based model such as (Vaswani et al. 2017; Veličković et al. 2018) as we get a
notion of closeness(similarity) in the embedding space. The final homogeneous
representation for node $i$ ($v_{hri}$) is obtained as:

$\displaystyle v_{hri}=W_{hr2}\ \ \sigma(W_{hr1}(sin(v_{gli})\|cos(v_{gli})))$

where $W_{hr2}\in R^{d\times d},W_{hr1}\in R^{d\times(d_{g}+d_{l}+d_{t})}$ are
learnable weights and $\sigma$ is the activation function.

Aggregation(AM): Once we achieve homogeneous representations for the features,
we can use these in the downstream task by applying a layer of MLP. However,
we also perform another level of aggregation to learn effective
representations. While, in principle, we could use any of the existing message
passing frameworks for this aggregation, we adopt a sparse variant of the
attention mechanism inspired by (Vaswani et al. 2017) for computational
benefits. Specifically, consider $X\in R^{N\times d_{f}}$ to be the node
feature learned from the spatial and spectral modules. Now we define for the
$l$-th layer, $W^{l}_{Q},W^{l}_{K}\in R^{d_{out}\times d_{f}}$ to be the
learnable weight matrices for the query and key of the self attention
respectively. We apply self attention on the transformed features followed by
softmax to get the aggregation weights
$w_{ij}^{l}=softmax(\sum_{d_{k}}\hat{w_{ij}^{l}})$ between nodes $i,j$, where
$\hat{w_{ij}^{l}}=\frac{W^{l}_{Q}X[i]^{T}\odot W^{l}_{K}X[j]^{T}}{d_{out}}$ if
nodes $i$ and $j$ are connected in the graph and 0 otherwise.

### Overall complexity

The complexity of spectral module is
$\mathcal{O}(|E|+N\sum_{j=0}^{J}M_{j})$(cf., appendix) where $M_{j}$ is the
order of the polynomial of the $j$-th filter head. The spatial modules can
compute the features in a $\mathcal{O}(|E|+N)$ complexity. The integration
module further has two components: HRM and AM. The HRM module would have a
computational complexity of $\mathcal{O}(N)$ whereas for the AM, it depends on
the underlying aggregator. In our choices, it would be an $\mathcal{O}(|E|+N)$
complexity. Thus the overall computational complexity comes to
$\mathcal{O}(3|E|+N(3+\sum_{j=0}^{J}M_{j}))=\mathcal{O}(|E|+N)$, for bounded
degree graphs further reduces to $\mathcal{O}(N)$ for a given snapshot at time
$t$.

## Experiments

Now, we present comprehensive experiments to evaluate our proposed framework.
We borrow datasets and its preprocessing/splitting settings used in previous
best baselines (Xiang, Huang, and Wang 2022; Pareja et al. 2020) with an
average of five runs for final values. Code will be released.

Datasets: Table 1 summarizes eight datasets for link prediction, edge
classification, and node classification. Each dataset contains a sequence of
time-ordered graphs. SBM is a synthetic dataset to simulate evolving community
structures. BC-OTC is a who-trusts-whom transaction network where the node
represents users, and the edges represent the ratings that a member gives
others in a range of -10(maximum distrust) to +10(maximum trust). BC-Alpha is
similar to BC-OTC, albeit the transactions are on a different bitcoin network.
UCI dataset is a student community network where nodes represent the students,
and the edges represent the messages exchanged between them. AS dataset
summarizes a temporal communication network indicating traffic flow between
routers. In the Reddit dataset, nodes are the source and target subreddits,
and the edges represent the sentiment’s polarity between users. The Elliptic
dataset represents licit vs illicit transactions on the elliptic network of
bitcoin transactions. The nodes represent transactions and edges represent
payment flows. Finally, for Brain dataset, nodes represent tiny regions/cubes
in the brain, and the edges are their connectivity.

Table 1: Data set statistics and details. In the task column, LP refers to Link Prediction, EC refers to Edge Classification and NC is Node Classification. | # Nodes | # Edges | # Time Steps | Task
---|---|---|---|---
|  |  | (Train / Val / Test) |
BC-OTC | 5,881 | 35,588 | 95 / 14 / 28 | EC
BC-Alpha | 3,777 | 24,173 | 95 / 13 / 28 | EC
Reddit | 55,863 | 858,490 | 122 / 18 / 34 | EC
SBM | 1,000 | 4,870,863 | 35 / 5 / 10 | LP
UCI | 1,899 | 59,835 | 62 / 9 / 17 | LP
AS | 6,474 | 13,895 | 70 / 10 / 20 | LP
Elliptic | 203,769 | 234,355 | 31 / 5 / 13 | NC
Brain | 5,000 | 1,955,488 | 10 / 1 / 1 | NC

Baselines: In order to show the efficacy of our framework, we compare with the
following competitive baselines: (1) Static graph methods: GCN (Kipf and
Welling 2017), GAT (Veličković et al. 2018). (2) Dynamic graph methods such as
GCN-GRU (Pareja et al. 2020), DynGEM (Goyal et al. 2017), GAEN (Shi et al.
2021), the variants of dyngraph2vec (Goyal, Chhetri, and Canedo 2020):
dyngraph2vecAE(v1) and dyngraph2vecAERNN (v2), the best performing two
variants of EvolveGCN (Pareja et al. 2020): EvolveGCN-O and EvolveGCN-H, and
the two variants of LEDG (Xiang, Huang, and Wang 2022): LEDGCN and LEDGAT.

Task and Experiment Settings: For the link prediction task, we report mean
reciprocal rank(MRR) and mean average precision(MAP) as in the baselines. For
MRR, the existing links shall be ranked higher than non existing links. For
each of the given nodes $u,v$, our model gives embeddings in the final layer
$h_{u}^{t},h_{v}^{t}$ respectively. We predict the link between these nodes by
concatenating the vectors for the corresponding nodes. Following the standard
practices, we perform negative sampling and optimize the cross entropy loss
function. For edge classification, our model classifies the edges between two
nodes($u,v$) by concatenating the vectors for the corresponding
nodes($h_{u}^{t},h_{v}^{t}$) at time $t$ and we report the micro F1 score. For
node classification similar to the baselines, our model reports micro F1
(Brain dataset) or minority class F1 (Elliptic dataset) by classifying given
node nodes($u$) as belonging to a certain class at time $t$.

We perform a search for the best model from a set of hyperparameters the range
of which are as follows: Number of layers $\in\\{1,2\\}$, Hidden dimension
$\in\\{32,64,128\\}$, Number of heads $\in\\{4,8,16\\}$, Filter order
$\in\\{4,8,16\\}$, Wavelet scales $\in\left[0.1,10\right]$. The rest of the
parameters and settings are borrowed from previous works (Pareja et al. 2020;
Xiang, Huang, and Wang 2022). All our experiments run on a single Tesla P100
GPU. The code will be public post-acceptance.

Variants of our Framework: In the main results we show three variants of our
method with different aggregators. We use the MLP(DEFT-MLP), GAT(DEFT-GAT) and
sparse Transformer(DEFT-T) as aggregators. The rationale behind selecting GAT
and MLP in addition to the proposed transformer aggregator is to contrast it
with a less expressive static attention mechanism in GAT and a control in MLP.

Table 2: Link prediction results where mean average precision (MAP) and mean
reciprocal rank (MRR) are displayed. Best values are bold, second bests are
underlined.

Datasets | SBM | UCI | AS
---|---|---|---
Metrics | MAP | MRR | MAP | MRR | MAP | MRR
GCN | 0.1894 | 0.0136 | 0.0001 | 0.0468 | 0.0019 | 0.1814
GAT | 0.1751 | 0.0128 | 0.0001 | 0.0468 | 0.0200 | 0.1390
DynGEM | 0.1680 | 0.0139 | 0.0209 | 0.1055 | 0.0529 | 0.1028
GCN-GRU | 0.1898 | 0.0119 | 0.0114 | 0.0985 | 0.0713 | 0.3388
dyngraph2vec V1 | 0.0983 | 0.0079 | 0.0044 | 0.0540 | 0.0331 | 0.0698
dyngraph2vec V2 | 0.1593 | 0.0120 | 0.0205 | 0.0713 | 0.0711 | 0.0493
GAEN | 0.1828 | 0.0078 | 0.0001 | 0.0490 | 0.1304 | 0.0507
EvolveGCN-H | 0.1947 | 0.0141 | 0.0126 | 0.0899 | 0.1534 | 0.3632
EvolveGCN-O | 0.1989 | 0.0138 | 0.0270 | 0.1379 | 0.1139 | 0.2746
LEDG-GCN | 0.1960 | 0.0147 | 0.0324 | 0.1626 | 0.1932 | 0.4694
LEDG-GAT | 0.1822 | 0.0123 | 0.0261 | 0.1492 | 0.2329 | 0.3835
DEFT-MLP(ours) | 0.1658 | 0.0124 | 0.0543 | 0.1715 | 0.1737 | 0.3569
DEFT-GAT(ours) | 0.0966 | 0.0081 | 0.0502 | 0.1702 | 0.0308 | 0.0928
DEFT-T(ours) | 0.2421 | 0.0220 | 0.0501 | 0.2007 | 0.5879 | 0.6471

Figure 3: Performance of edge classification and node classification, with F1
scores on the y-axis. For BC-OTC, BC-Alpha, and Reddit (edge classification)
and Brain(node classification) the F1 score is the micro average. For Elliptic
dataset(node classification), the F1 is the minority F1 because of the class
imbalance and interest in the minority class (illicit transactions).

## Results

### Results on Link Prediction

Table 2 summarizes link prediction results. The key results are: 1) across
datasets, our model (DEFT-T) significantly outperform all baselines which
focus on learning local dependencies (in static and dynamic graph settings).
It illustrates our framework’s effectiveness in learning homogeneous
representations of local and global dependencies in evolving graphs. 2)
Interestingly, a simple MLP aggregator (DEFT-MLP) on the learned spectral and
spatial properties already achieves better results than most of the baselines
on UCI and AS datasets (including ones for dynamic graphs). 3) The transformer
variant (DEFT-T) gives consistently better results than the GAT and MLP
variants indicating the aggregation using the sparse self-attention of
transformers is essential for the framework on these datasets.

### Results on Edge Classification

The results of edge classification on the three datasets (BC-OTC, BC-Alpha,
Reddit) are given in Figure 3 and show the micro averaged F1 score. We observe
from the results that at least one variant of our method outperforms the best
baselines variants: EvolveGCN and LEDG. Our MLP-based configuration achieves
better results than baselines on BC-OTC and BC-Alpha, whereas DEFT-T
consistently outperforms the baselines and other DEFT variants. This further
verifies the efficacy of the model developed under the proposed framework.

### Results on Node Classification

Figure 3 shows the results of node classification on Brain(homophily ratio
(Zhu et al. 2020): 26%) and Elliptic(homophily ratio: 96%) datasets. DEFT
would consider the brain structure on the Brain dataset while capturing global
interactions. DEFT has a substantial performance gain on the Brain dataset,
which has a high heterophily ratio requiring capturing global interactions in
graphs compared to existing models that only perform local aggregation. On the
homophilic Elliptic dataset, DEFT also performs good by capturing low-
frequency signals. It confirms that when the spatial and spectral properties
are combined using an optimal aggregator to learn the local and global
dependencies, the model results in robust performance across tasks.

### Ablation Studies

Table 3: Table shows the link prediction results for the static wavelet and
transformer architectures along with ablations of individual components of our
model.

Datasets | SBM | UCI
---|---|---
Metrics | MAP | MRR | MAP | MRR
GWNN | 0.1789 | 0.0121 | 0.0076 | 0.0820
Transformer | 0.2052 | 0.0174 | 0.0308 | 0.1441
GT | 0.2166 | 0.01805 | 0.0310 | 0.1414
SAN | 0.2143 | 0.0180 | 0.0388 | 0.1822
DEFT-woSpectral | 0.2123 | 0.0171 | 0.0460 | 0.1764
DEFT-woSpatial | 0.2083 | 0.0171 | 0.0421 | 0.1806
DEFT-woHRM | 0.2410 | 0.0217 | 0.0429 | 0.1588
DEFT-staticSpectral | 0.2297 | 0.0198 | 0.0419 | 0.1506
DEFT-T(best configuration) | 0.2421 | 0.0220 | 0.0501 | 0.2007

Model Ablation Study: To study the effect of each module of DEFT, we
systematically remove modules and create several sub-configuration of our best
performing model (DEFT-T): (i) DEFT-woSpectral does not have the spectral
module (ii) DEFT-woSpatial does not have the spatial module (iii) DEFT-woHRM
with the Homogeneous Representation Module removed(the spatial and spectral
features are simply added), (iv) DEFT-staticSpectral has the spectral wavelets
that are learnable but do not evolve with time and in a graph-specific manner.
We also compare our best model with recent transformer-based models used for
static graphs such as transformer (Vaswani et al. 2017), SAN (Kreuzer et al.
2021), GT (Dwivedi and Bresson 2020), and static wavelet baseline GWNN (Xu et
al. 2019a).

The ablation study (Table 3) provides several key insights to understand
proposed DEFT framework. Firstly, the lesser values reported by static (non-
learnable) wavelet GNN (i.e., GWNN) and the static (but learnable) variant of
our model (i.e., DEFT-staticSpectral) conclude that learning dynamic wavelets
benefits in the case of evolving graphs. Secondly, the significant drop in
performance when the spatial or spectral module is removed provides an
essential finding that for dynamic graphs, these modules learn useful and
orthogonal representations. When these modules are combined using an
effectively learned homogeneous representation, the model (DEFT-T) provides
significantly superior performance than its variants. Lastly, the static
transformer models consistently report lesser values than our model.

Figure 4: Effect of the parameters wavelet scale and filter order on the
performance of the model.

Filter parameter selection: In spectral module, there are two parameters to be
tuned: filter(polynomial) order and the wavelet scale. In this study, we try
to understand how selection of parameters affects model performance. From
Figure 4 illustrating results on Brain and Elliptic datasets, we make the
following observations: (1) Wavelet scales: In Brain dataset, the result dips
with higher scale and then increases. For Elliptic, a general increasing trend
is observed. In the case of the learnt filters for the given graphs, it seems
that as the scales decrease the aggregation range of the spectral GNN
increases(as in figure 1). Thus, brain dataset being heterophilic seems to
benefit from scales in the lower region. Larger scales may work well for the
homophilic Elliptic dataset as the aggregation focuses on a relatively smaller
neighborhood, while also capturing the heterophily signals in the minority
illicit transactions. (2) Filter order: On Brain dataset, there is a
decreasing trend with increasing filter order, whereas, for Elliptic, an
optimal result is attained in the middle region. This could be because the
Brain dataset being small needs simpler filters(such as single mode band
reject etc.); with higher filter order, it may have overfit. The larger
Elliptic dataset, on the other hand, benefits from a larger filter order.
Despite being homophilic for the majority nodes, it may need complex filters
to handle the minority nodes due to the peculiar and global nature of illicit
transactions. Thus, it benefits from relatively higher filter orders. These
observations conclude that the parameters of the order of polynomial and the
wavelet scales are dataset dependent.

## Conclusion

Our work proposes a concept of learnable spectral graph wavelets to capture
global dependencies by omitting the need for full eigen decomposition of
dynamic graphs. Furthermore, we implement it in the DEFT framework that
integrates the graph’s spectral and spatial properties for dynamic graphs.
From the results on a wide range of tasks for dynamic graphs, we infer that
the proposed method is able to capture the local(short range) and
global(extended range) properties effectively. Also, it can capture a broad
range of interactions respecting the graph structure, thus avoiding noise
while being sparse(computationally efficient). Future works could build upon
frameworks for joint time and graph Fourier transform for dynamic graphs to
generate wavelets and study if these could be used for graph sampling to
provide computational efficiency on extremely large graphs.

## Acknowledgments

This work was partly supported by JSPS KAKENHI Grant Number JP21K21280. We
thank anonymous reviewers for constructive feedback.

## References

  * Balcilar et al. (2020) Balcilar, M.; Renton, G.; Héroux, P.; Gaüzère, B.; Adam, S.; and Honeine, P. 2020.  Analyzing the expressive power of graph neural networks in a spectral perspective.  In _International Conference on Learning Representations_.
  * Barron (1991) Barron, A. 1991.  Approximation and Estimation Bounds for Artificial Neural Networks.  _Machine Learning_ , 14: 115–133.
  * Bastos et al. (2022) Bastos, A.; Nadgeri, A.; Singh, K.; Kanezashi, H.; Suzumura, T.; and Mulang’, I. O. 2022.  How Expressive are Transformers in Spectral Domain for Graphs?  _Transactions on Machine Learning Research_.
  * Dwivedi and Bresson (2020) Dwivedi, V. P.; and Bresson, X. 2020.  A Generalization of Transformer Networks to Graphs.  arXiv:2012.09699.
  * Goyal, Chhetri, and Canedo (2020) Goyal, P.; Chhetri, S. R.; and Canedo, A. 2020.  dyngraph2vec: Capturing network dynamics using dynamic graph representation learning.  _Knowl. Based Syst._ , 187.
  * Goyal et al. (2017) Goyal, P.; Kamra, N.; He, X.; and Liu, Y. 2017.  DynGEM: Deep Embedding Method for Dynamic Graphs.  _IJCAI Workshop on Representation Learning for Graphs_.
  * Hamilton, Ying, and Leskovec (2017) Hamilton, W.; Ying, Z.; and Leskovec, J. 2017.  Inductive representation learning on large graphs.  In _Advances in neural information processing systems_ , 1024–1034.
  * Hammond, Vandergheynst, and Gribonval (2011) Hammond, D. K.; Vandergheynst, P.; and Gribonval, R. 2011.  Wavelets on graphs via spectral graph theory.  _Applied and Computational Harmonic Analysis_ , 30(2): 129–150.
  * Hammond, Vandergheynst, and Gribonval (2019) Hammond, D. K.; Vandergheynst, P.; and Gribonval, R. 2019.  The spectral graph wavelet transform: Fundamental theory and fast computation.  In _Vertex-Frequency Analysis of Graph Signals_ , 141–175. Springer.
  * Han et al. (2022) Han, K.; Wang, Y.; Guo, J.; Tang, Y.; and Wu, E. 2022.  Vision GNN: An Image is Worth Graph of Nodes.  _arXiv preprint arXiv:2206.00272_.
  * Hornik (1991) Hornik, K. 1991.  Approximation capabilities of multilayer feedforward networks.  _Neural networks_ , 4(2): 251–257.
  * Kazemi et al. (2020) Kazemi, S. M.; Goel, R.; Jain, K.; Kobyzev, I.; Sethi, A.; Forsyth, P.; and Poupart, P. 2020.  Representation learning for dynamic graphs: A survey.  _J. Mach. Learn. Res._ , 21(70): 1–73.
  * Khoshraftar and An (2022) Khoshraftar, S.; and An, A. 2022.  A Survey on Graph Representation Learning Methods.  _arXiv preprint arXiv:2204.01855_.
  * Kipf and Welling (2017) Kipf, T. N.; and Welling, M. 2017.  Semi-Supervised Classification with Graph Convolutional Networks.  In _5th International Conference on Learning Representations, ICLR 2017_.
  * Kreuzer et al. (2021) Kreuzer, D.; Beaini, D.; Hamilton, W. L.; Létourneau, V.; and Tossou, P. 2021\.  Rethinking Graph Transformers with Spectral Attention.  _NeurlPS 2021_ , abs/2106.03893.
  * Levie et al. (2018) Levie, R.; Monti, F.; Bresson, X.; and Bronstein, M. M. 2018.  Cayleynets: Graph convolutional neural networks with complex rational spectral filters.  _IEEE Transactions on Signal Processing_ , 67(1): 97–109.
  * Li et al. (2021) Li, Y.; Si, S.; Li, G.; Hsieh, C.-J.; and Bengio, S. 2021.  Learnable Fourier Features for Multi-dimensional Spatial Positional Encoding.  In Beygelzimer, A.; Dauphin, Y.; Liang, P.; and Vaughan, J. W., eds., _Advances in Neural Information Processing Systems_.
  * Manessi, Rozza, and Manzo (2020) Manessi, F.; Rozza, A.; and Manzo, M. 2020.  Dynamic graph convolutional networks.  _Pattern Recognit._ , 97.
  * Narayan and Roe (2018) Narayan, A.; and Roe, P. H. 2018.  Learning graph dynamics using deep neural networks.  _IFAC-PapersOnLine_ , 51(2): 433–438.
  * Oppenheim and Schafer (1975) Oppenheim, A. V.; and Schafer, R. W. 1975.  _Digital signal processing_.  Prentice-Hall international editions. Prentice-Hall.  ISBN 0132141078.
  * Pareja et al. (2020) Pareja, A.; Domeniconi, G.; Chen, J.; Ma, T.; Suzumura, T.; Kanezashi, H.; Kaler, T.; Schardl, T. B.; and Leiserson, C. E. 2020.  EvolveGCN: Evolving Graph Convolutional Networks for Dynamic Graphs.  In _The Thirty-Fourth AAAI Conference on Artificial Intelligence, AAAI 2020_ , 5363–5370. AAAI Press.
  * Prasad and Kan (2019) Prasad, A.; and Kan, M.-Y. 2019.  Glocal: Incorporating global information in local convolution for keyphrase extraction.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , 1837–1846.
  * Rahimi and Recht (2007) Rahimi, A.; and Recht, B. 2007.  Random Features for Large-Scale Kernel Machines.  In _Advances in Neural Information Processing Systems 20, Proceedings of the Twenty-First Annual Conference on Neural Information Processing Systems_ , 1177–1184. Curran Associates, Inc.
  * Rahimi and Recht (2008) Rahimi, A.; and Recht, B. 2008.  Weighted Sums of Random Kitchen Sinks: Replacing minimization with randomization in learning.  In _Advances in Neural Information Processing Systems 21, Proceedings of the Twenty-Second Annual Conference on Neural Information Processing Systems_ , 1313–1320. Curran Associates, Inc.
  * Seo et al. (2016) Seo, Y.; Defferrard, M.; Vandergheynst, P.; and Bresson, X. 2016.  Structured Sequence Modeling with Graph Convolutional Recurrent Networks.  _CoRR_ , abs/1612.07659.
  * Shi et al. (2021) Shi, M.; Huang, Y.; Zhu, X.; Tang, Y.; Zhuang, Y.; and Liu, J. 2021.  GAEN: Graph Attention Evolving Networks.  In _IJCAI_ , 1541–1547.
  * Srinivas et al. (2021) Srinivas, A.; Lin, T.-Y.; Parmar, N.; Shlens, J.; Abbeel, P.; and Vaswani, A. 2021\.  Bottleneck transformers for visual recognition.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 16519–16529.
  * Stärk et al. (2022) Stärk, H.; Ganea, O.; Pattanaik, L.; Barzilay, R.; and Jaakkola, T. S. 2022\.  EquiBind: Geometric Deep Learning for Drug Binding Structure Prediction.  _CoRR_ , abs/2202.05146.
  * Tremblay and Borgnat (2014) Tremblay, N.; and Borgnat, P. 2014.  Graph wavelets for multiscale community mining.  _IEEE Transactions on Signal Processing_ , 62(20): 5227–5239.
  * Vaswani et al. (2017) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, Ł.; and Polosukhin, I. 2017.  Attention is all you need.  In _Advances in neural information processing systems_ , 5998–6008.
  * Veličković et al. (2018) Veličković, P.; Cucurull, G.; Casanova, A.; Romero, A.; Liò, P.; and Bengio, Y. 2018.  Graph Attention Networks.  In _International Conference on Learning Representations_.
  * Wu et al. (2021a) Wu, L.; Chen, Y.; Ji, H.; and Liu, B. 2021a.  Deep Learning on Graphs for Natural Language Processing.  In _Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval_ , 2651–2653. Association for Computing Machinery.
  * Wu et al. (2021b) Wu, Z.; Jain, P.; Wright, M.; Mirhoseini, A.; Gonzalez, J. E.; and Stoica, I. 2021b.  Representing long-range context for graph neural networks with global attention.  _Advances in Neural Information Processing Systems_ , 34: 13266–13279.
  * Xiang, Huang, and Wang (2022) Xiang, X.; Huang, T.; and Wang, D. 2022.  Learning to Evolve on Dynamic Graphs (SA).  In _Thirty-Sixth AAAI Conference on Artificial Intelligence, AAAI 2022_. AAAI Press.
  * Xu et al. (2019a) Xu, B.; Shen, H.; Cao, Q.; Qiu, Y.; and Cheng, X. 2019a.  Graph Wavelet Neural Network.  In _International Conference on Learning Representations_.
  * Xu et al. (2019b) Xu, D.; Cheng, W.; Luo, D.; Gu, Y.; Liu, X.; Ni, J.; Zong, B.; Chen, H.; and Zhang, X. 2019b.  Adaptive Neural Network for Node Classification in Dynamic Networks.  In Wang, J.; Shim, K.; and Wu, X., eds., _2019 IEEE International Conference on Data Mining, ICDM 2019, Beijing, China, November 8-11, 2019_ , 1402–1407. IEEE.
  * Xu et al. (2019c) Xu, K.; Hu, W.; Leskovec, J.; and Jegelka, S. 2019c.  How Powerful are Graph Neural Networks?  arXiv:1810.00826.
  * Xu, Singh, and Karniadakis (2022) Xu, M.; Singh, A. V.; and Karniadakis, G. E. 2022.  DynG2G: An Efficient Stochastic Graph Embedding Method for Temporal Graphs.  _IEEE Transactions on Neural Networks and Learning Systems_.
  * Yang et al. (2022) Yang, M.; Shen, Y.; Li, R.; Qi, H.; Zhang, Q.; and Yin, B. 2022.  A New Perspective on the Effects of Spectrum in Graph Neural Networks.  In _International Conference on Machine Learning_ , 25261–25279. PMLR.
  * Ying et al. (2021) Ying, C.; Cai, T.; Luo, S.; Zheng, S.; Ke, G.; He, D.; Shen, Y.; and Liu, T. 2021\.  Do Transformers Really Perform Badly for Graph Representation?  In _Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021_ , 28877–28888.
  * Zhu et al. (2020) Zhu, J.; Yan, Y.; Zhao, L.; Heimann, M.; Akoglu, L.; and Koutra, D. 2020.  Beyond Homophily in Graph Neural Networks: Current Limitations and Effective Designs.  In _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020_.

## Supplementary Material

In this section we outline the proofs for the theoretical results stated in
the main paper. For completion we give the statement along with the proofs.

######  Lemma 3.

Consider $G^{t}(\lambda)$ to be the filter function at time $t$. Assume the
Markov property to hold in the dynamic setting where the desired filter
function($G^{t+1}(\Lambda)$) at time $t+1$ depends on the past state at time
$t$($G^{t}(\Lambda)$). Consider this mapping between the past state and the
current states to be captured by an arbitrary functional $f$ such that
$G^{t+1}(\lambda)=f(G^{t}(\lambda_{1}),G^{t}(\lambda_{2}),\dots
G^{t}(\lambda_{N}),\lambda)$ and we assume $f$ to be $L$ lipschitz continuous.
Further, let $C_{t}=U_{t}G(\lambda)^{t}U_{t}^{T}\in R^{N\times N}$ represent
the convolution support of the desired spectral filter and $C_{t}^{a}$ be the
learnt convolution support at time $t$. Then, we have:

(i) $\left\lVert C_{t+1}^{a}-C_{t+1}\right\rVert_{F}\leq
LN^{2}\sqrt{\left\lVert
C_{t}^{a}-C_{t}\right\rVert_{F}^{2}+\epsilon_{ca}^{2}}+\epsilon_{fa}$

(ii) $\left\lVert C_{t+1}^{a}-C_{t}^{a}\right\rVert_{F}\leq\left\lVert
C_{t+1}-C_{t}\right\rVert_{F}+2\sqrt{N}\epsilon_{ca}$

where $\epsilon_{ca}$ and $\epsilon_{fa}$ are the filter polynomial(Chebyshev)
and function approximation errors and depends on the order of the polynomial,
number of training samples, model size etc.

###### Proof.

(i) The approximate and desired convolution support matrix $C_{t+1}^{a}$ and
$C_{t+1}$ at time $t$ of the spectral filter has the decomposition as:

$\displaystyle C_{t+1}^{a}$
$\displaystyle=U_{t+1}G^{a}_{t+1}(\Lambda)U_{t+1}^{T}$ $\displaystyle C_{t+1}$
$\displaystyle=U_{t+1}G_{t+1}(\Lambda)U_{t+1}^{T}$

where $U_{t+1}\in R^{N\times N}$ is the matrix whose columns contain the
eigenvectors of the graph laplacian at time $t+1$. Now the Frobenius norms of
these two matrices is given by,

$\displaystyle\left\lVert C_{t+1}^{a}-C_{t+1}\right\rVert$
$\displaystyle=\left\lVert
U_{t+1}G^{a}_{t+1}(\Lambda)U_{t+1}^{T}-U_{t+1}G_{t+1}(\Lambda)U_{t+1}^{T}\right\rVert$
$\displaystyle=\left\lVert
U_{t+1}(G^{a}_{t+1}(\Lambda)-G_{t+1}(\Lambda))U_{t+1}^{T}\right\rVert$

As the graph evolves at time $t$, let $f_{t}$ represent the function that maps
the desired spectra at time $t$ to the one at time $t+1$ i.e.
$G_{t+1}(\lambda_{i})=f_{t}(\|_{i\in[N]}G^{t}(\lambda_{i}),\lambda_{i}),\forall
i\in[N]$. Let $f_{t}^{a}$ represent a similar mapping from
$G_{t}^{a}(\Lambda)$ to $G_{t+1}^{a}(\Lambda)$. Using the approximation result
of multi-layer neural networks in (Barron 1991), we have
$f_{t}^{a}(\|_{i\in[N]}G^{t}(\lambda_{i}),\lambda_{i})-f_{t}(\|_{i\in[N]}G^{t}(\lambda_{i}),\lambda_{i})=\mathcal{O}(\frac{C_{f}^{2}}{h}+\frac{hd}{N}log(N))=\epsilon_{1}$
where $C_{f}$ is the first absolute moment of the distribution of the Fourier
magnitudes of $f_{t}^{a}$, $h$ the number of hidden units in the network, $d$
is the input dimension and $N_{s}$ the number of training samples, $\|$ is the
concatenation operator. Thus we have

$f_{t}^{a}(\|_{i\in[N]}G^{t}(\lambda_{i}),\lambda)=f_{t}(\|_{i\in[N]}G^{t}(\lambda_{i}),\lambda)+\epsilon_{1}$
(5)

where $\epsilon_{1}=\mathcal{O}(\frac{C_{f}^{2}}{h}+\frac{hd}{N}log(N))$.

Now, from the above equations we upper bound the norm(Frobenius unless
mentioned otherwise) as,

$\displaystyle\left\lVert C_{t+1}^{a}-C_{t+1}\right\rVert=\left\lVert
U_{t+1}(G^{a}_{t+1}(\Lambda)-G_{t+1}(\Lambda))U_{t+1}^{T}\right\rVert$
$\displaystyle=\left\lVert
U_{t+1}(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda))U_{t+1}^{T}\right\rVert$
$\displaystyle=\left\lVert
U_{t+1}(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)+\epsilon_{1}I_{N}-f_{t}(G_{t}(\Lambda),\Lambda))U_{t+1}^{T}\right\rVert$
$\displaystyle=\left\lVert
U_{t+1}(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda))U_{t+1}^{T}+U_{t+1}(\epsilon_{1}I_{N}U_{t+1}^{T})\right\rVert$
$\displaystyle\leq\left\lVert
U_{t+1}(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda))U_{t+1}^{T}\right\rVert+\left\lVert\epsilon_{1}U_{t+1}I_{N}U_{t+1}^{T}\right\rVert$
$\displaystyle=\left\lVert
U_{t+1}(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda))U_{t+1}^{T}\right\rVert+\left\lVert\epsilon_{1}U_{t+1}U_{t+1}^{T}\right\rVert$
$\displaystyle=\left\lVert
U_{t+1}(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda))U_{t+1}^{T}\right\rVert+\epsilon_{1}\left\lVert
I_{N}\right\rVert$ $\displaystyle=\left\lVert
U_{t+1}(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda))U_{t+1}^{T}\right\rVert+\sqrt{N}\epsilon_{1}$
$\displaystyle=\left\lVert
U_{t+1}(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda))U_{t+1}^{T}\right\rVert+\epsilon_{fa}$

where
$\epsilon_{fa}=\sqrt{N}\mathcal{O}(\frac{C_{f}^{2}}{h}+\frac{hd}{N_{s}}log(N_{s}))$
and $f^{a}_{t},f_{t}$ are applied element wise over the domain($\Lambda$ being
discrete). Assuming, $G^{a}_{t}(\Lambda)\approx G_{t}(\Lambda)$(which holds
true asymptotically) we perform a first order Taylor series expansion of the
function $f_{t}(G^{a}_{t}(\Lambda),\Lambda)$ at
$(G^{a}_{t}(\Lambda),\Lambda)=\|_{i\in[N]}G^{t}(\lambda_{i}),\lambda_{i}$ as
below

$\displaystyle f_{t}(G^{a}_{t}(\Lambda),\Lambda)$ $\displaystyle\approx
f_{t}(G_{t}(\Lambda),\Lambda)+\left<(G^{a}_{t}(\Lambda)-G_{t}(\Lambda),0),f^{{}^{\prime}}_{t}(G_{t}(\Lambda),\Lambda)\right>$

where $f^{{}^{\prime}}_{t}(G_{t}(\Lambda),\Lambda)\in R^{N}$ is the vector of
partial derivatives along each basis. Thus in the limit that
$G^{a}_{t}(\Lambda)\xrightarrow[]{}G_{t}(\Lambda)$,

$\displaystyle\underset{G^{a}_{t}(\Lambda)\xrightarrow[]{}G_{t}(\Lambda)}{lim}f_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda)$
$\displaystyle=\left<(G^{a}_{t}(\Lambda)-G_{t}(\Lambda),0),f^{{}^{\prime}}_{t}(G_{t}(\Lambda),\Lambda)\right>$
$\displaystyle\therefore\left\lVert
f_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda)\right\rVert$
$\displaystyle=\left\lVert\left<(G^{a}_{t}(\Lambda)-G_{t}(\Lambda),0),f^{{}^{\prime}}_{t}(G_{t}(\Lambda),\Lambda)\right>\right\rVert$
$\displaystyle\leq\left\lVert\left<(G^{a}_{t}(\Lambda)-G_{t}(\Lambda),0),L\mathbf{1}\right>\right\rVert$
$\displaystyle\leq
L\left\lVert\sum_{i\in[N]}G^{a}_{t}(\lambda_{i})-G_{t}(\lambda_{i})\right\rVert$

Now since $G^{a}_{t}(\lambda_{i})$ is obtained from a truncated
polynomial(Chebyshev) approximation which is parameterized using a (universal)
function approximator, we would have error terms involving the polynomial as
well as function approximation. From (Bastos et al. 2022) which proved such
bounds for the static case we have,

$G^{a}_{t}(\lambda_{i})-G_{t}(\lambda_{i})=\mathcal{O}(\frac{KC_{f}^{2}}{h}+\frac{hK^{2}}{N_{s}}log(N_{s})+K^{-m})$

Here it is assumed that the input dimension is $\mathcal{O}(K)$ which is the
filter order and the desired filter response $G(\lambda)$ has $m+1$ continuous
derivatives on the domain $[-1,1]$. We set
$\mathcal{O}(\frac{KC_{f}^{2}}{h}+\frac{hK^{2}}{N_{s}}log(N_{s})+K^{-m})=\epsilon_{ca}$
for notational simplicity. Now we have that,

$\displaystyle\left\lVert\sum_{i\in[N]}G^{a}_{t}(\lambda_{i})-G_{t}(\lambda_{i})\right\rVert^{2}=\left\lVert
G^{a}_{t}(\Lambda)-G_{t}(\Lambda)\right\rVert^{2}$
$\displaystyle+2\sum_{i\in[N],i\neq
j}(G^{a}_{t}(\lambda_{i})-G_{t}(\lambda_{i}))(G^{a}_{t}(\lambda_{j})-G_{t}(\lambda_{j}))$
$\displaystyle=\left\lVert
G^{a}_{t}(\Lambda)-G_{t}(\Lambda)\right\rVert^{2}+2N^{2}\epsilon_{ca}^{2}$
$\displaystyle\left\lVert\sum_{i\in[N]}G^{a}_{t}(\lambda_{i})-G_{t}(\lambda_{i})\right\rVert=\sqrt{\left\lVert
G^{a}_{t}(\Lambda)-G_{t}(\Lambda)\right\rVert^{2}+2N^{2}\epsilon_{ca}^{2}}$

From the above equation we have,

$\displaystyle\left\lVert
f_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda)\right\rVert\leq
L\left\lVert\sum_{i\in[N]}G^{a}_{t}(\lambda_{i})-G_{t}(\lambda_{i})\right\rVert$
$\displaystyle=L\sqrt{\left\lVert
G^{a}_{t}(\Lambda)-G_{t}(\Lambda)\right\rVert^{2}+2N^{2}\epsilon_{ca}^{2}}$
$\displaystyle=L\sqrt{\left\lVert
U_{t}^{T}U_{t}(G^{a}_{t}(\Lambda)-G_{t}(\Lambda))U_{t}^{T}U_{t}\right\rVert^{2}+2N^{2}\epsilon_{ca}^{2}}$
$\displaystyle=L\sqrt{\left\lVert
U_{t}^{T}(U_{t}G^{a}_{t}(\Lambda)U_{t}^{T}-U_{t}G_{t}(\Lambda)U_{t}^{T})U_{t}\right\rVert^{2}+2N^{2}\epsilon_{ca}^{2}}$
$\displaystyle\leq L\sqrt{\left\lVert
U_{t}^{T}\right\rVert^{2}\left\lVert(U_{t}G^{a}_{t}(\Lambda)U_{t}^{T}-U_{t}G_{t}(\Lambda)U_{t}^{T})\right\rVert^{2}\left\lVert
U_{t}\right\rVert^{2}+2N^{2}\epsilon_{ca}^{2}}$ $\displaystyle\leq
L\sqrt{N^{2}\left\lVert(U_{t}G^{a}_{t}(\Lambda)U_{t}^{T}-U_{t}G_{t}(\Lambda)U_{t}^{T})\right\rVert^{2}+2N^{2}\epsilon_{ca}^{2}}$
$\displaystyle=L\sqrt{N^{2}\left\lVert
C^{a}_{t}-C_{t}\right\rVert^{2}+2N^{2}\epsilon_{ca}^{2}}$
$\displaystyle=LN\sqrt{\left\lVert
C^{a}_{t}-C_{t}\right\rVert^{2}+2\epsilon_{ca}^{2}}$

Using this and the expression for $\left\lVert
C_{t+1}^{a}-C_{t+1}\right\rVert$ we have,

$\displaystyle\left\lVert C_{t+1}^{a}-C_{t+1}\right\rVert\leq\left\lVert
U_{t+1}(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda))U_{t+1}^{T}\right\rVert+\epsilon$
$\displaystyle\leq\left\lVert
U_{t+1}\right\rVert\left\lVert(f^{a}_{t}(G^{a}_{t}(\Lambda),\Lambda)-f_{t}(G_{t}(\Lambda),\Lambda))\right\rVert\left\lVert
U_{t+1}^{T}\right\rVert$ $\displaystyle\leq LN^{2}\sqrt{\left\lVert
C^{a}_{t}-C_{t}\right\rVert^{2}+2\epsilon_{ca}^{2}}+\epsilon_{fa}$

This completes the first part of the proof.

(ii) For the second part we have to upper bound $\left\lVert
C^{a}_{t+1}-C^{a}_{t}\right\rVert$. We have,

$\displaystyle\left\lVert C^{a}_{t+1}-C^{a}_{t}\right\rVert=\left\lVert
U_{t+1}G^{a}_{t+1}(\Lambda)U_{t+1}-U_{t}G^{a}_{t}(\Lambda)U_{t}^{T}\right\rVert$
$\displaystyle=\left\lVert
U_{t+1}(G_{t+1}+\epsilon_{ca}I_{N})(\Lambda)U_{t+1}-U_{t}(G_{t}(\Lambda)+\epsilon_{ca}I_{N})U_{t}^{T}\right\rVert$
$\displaystyle\leq\left\lVert
U_{t+1}G_{t+1}U_{t+1}^{T}-U_{t}G_{t}(\Lambda)U_{t}^{T}\right\rVert+2\epsilon_{ca}\left\lVert
I_{N}\right\rVert$ $\displaystyle=\left\lVert
C_{t+1}-C_{t}\right\rVert+2\sqrt{N}\epsilon_{ca}$

where $\epsilon_{ca}$ in the second step is the same as for the part (i). This
completes the proof.

∎

Part (i) in the above result gives us a relation between the error at times
$t$ and $t+1$. In the expression $LN^{2}\sqrt{\left\lVert
C_{t}^{a}-C_{t}\right\rVert_{F}^{2}+\epsilon_{ca}^{2}}+\epsilon_{fa}$ if the
error of approximating $C_{t}$ if larger than the chebyshev polynomial
approximation error i.e. $\left\lVert
C_{t}^{a}-C_{t}\right\rVert_{F}>>\epsilon_{ca}$, then we have,

$\displaystyle\left\lVert C_{t+1}^{a}-C_{t+1}\right\rVert_{F}$
$\displaystyle\leq LN^{2}\sqrt{\left\lVert
C_{t}^{a}-C_{t}\right\rVert_{F}^{2}+\epsilon_{ca}^{2}}+\epsilon_{fa}$
$\displaystyle\approx LN^{2}\sqrt{\left\lVert
C_{t}^{a}-C_{t}\right\rVert_{F}^{2}}+\epsilon_{fa}$
$\displaystyle=LN^{2}\left\lVert
C_{t}^{a}-C_{t}\right\rVert_{F}+\epsilon_{fa}$
$\displaystyle=(LN^{2})^{t}\left\lVert
C_{0}^{a}-C_{0}\right\rVert_{F}+\frac{\epsilon_{fa}}{1-LN^{2}}$

The last equation follows if $L<\frac{1}{N^{2}}$ in the limit of
$t\xrightarrow[]{}\infty$. Thus convergence is guaranteed if the filter
function is smooth. This is also one of the reasons why we initialise the
filter coeffcients with a one vector(all pass filter). However this may be too
strict a constraint in practice, and in the cases where the initial
approximation error is large and/or the filter function is complex with high
gradients the approximation error of the convolution support may diverge with
time. This issue arises as we have taken the Markov assumption, for
theoretical analysis, for determining the filter function at the next state
i.e. the filter function depends on the current filter function. However in
real world settings the assumption may not hold and empirically learning may
happen from other signals(eg: the current graph input etc.) and the filter
function may not depend on the past state entirely. This is where the LSTM(of
the dynamic parameter module) may help if it’s forget gate decides that the
current input is the useful signal for learning the filter function and not
the previous state, thus decoupling the learning of the filter at current time
step from the errors of the previous ones and empirically attaining a bounded
error with time. Similarly, part (ii) states that under the given assumptions,
the deviation between the approximated convolution supports at consecutive
time steps is directly related to that between the desired supports at the
respective time steps.

######  Lemma 4.

Let $G(\lambda)$ be the frequency response at frequency $\lambda$. Let
$\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{n}$ be the eigenvalues in
descending order and $p_{1},p_{2},\dots p_{n}$ be he corresponding eigenvalues
of the laplacian of the graph. Define $\lambda_{max}$ to be the eigenvalue at
which $G(\lambda)$ is maximum. Let $C^{l}$ represent the convolution
support($UG(\Lambda)U^{T}$) of the spectral filter at layer $l$. Then the
factor by which the cosine similarity between consecutive layers dampens is
$\underset{l\xrightarrow[]{}\infty}{lim}\frac{|cos(\left<C^{l+1}h,p_{n}\right>)|}{|cos(\left<C^{l}h,p_{n}\right>|}=\frac{G(\lambda_{n})}{G(\lambda_{max})}$.

###### Proof.

We begin by noting that the expression
$cos(\left<C^{l}h,p_{n}\right>)=\frac{\alpha_{n}G^{l}(\lambda_{n})}{\sqrt{\sum_{i=1}^{N}\alpha_{i}^{2}G^{2l}(\lambda_{i})}}$
as shown in (Yang et al. 2022), where $\alpha_{i}=\left<h,p_{i}\right>$.

Now, the ratio of these similarities between consecutive layers is,

$\displaystyle\frac{|cos(\left<C^{l+1}h,p_{n}\right>)|}{|cos(\left<C^{l}h,p_{n}\right>|}$
$\displaystyle=\frac{\alpha_{n}G^{l+1}(\lambda_{n})}{\sqrt{\sum_{i=1}^{N}\alpha_{i}^{2}G^{2(l+1)}(\lambda_{i})}}\frac{\sqrt{\sum_{i=1}^{N}\alpha_{i}^{2}G^{2l}(\lambda_{i})}}{\alpha_{n}G^{l}(\lambda_{n})}$
$\displaystyle=G(\lambda_{n})\sqrt{\frac{\sum_{i=1}^{N}\alpha_{i}^{2}G^{2l}(\lambda_{i})}{\sum_{i=1}^{N}\alpha_{i}^{2}G^{2(l+1)}(\lambda_{i})}}$

In the limit of the layers tending to $\infty$,

$\displaystyle\underset{l\xrightarrow[]{}\infty}{lim}\frac{|cos(\left<C^{l+1}h,p_{n}\right>)|}{|cos(\left<C^{l}h,p_{n}\right>|}=\underset{l\xrightarrow[]{}\infty}{lim}G(\lambda_{n})\sqrt{\frac{\sum_{i=1}^{N}\alpha_{i}^{2}G^{2l}(\lambda_{i})}{\sum_{i=1}^{N}\alpha_{i}^{2}G^{2(l+1)}(\lambda_{i})}}$
$\displaystyle=G(\lambda_{n})\sqrt{\frac{G^{2l}(\lambda_{max})\sum_{i=1}^{N}\alpha_{i}^{2}\underset{l\xrightarrow[]{}\infty}{lim}\left(\frac{G^{2l}(\lambda_{i})}{G^{2l}(\lambda_{max})}\right)}{\sum_{i=1}^{N}G^{2(l+1)}(\lambda_{max})\alpha_{i}^{2}\underset{l\xrightarrow[]{}\infty}{lim}\left(\frac{G^{2(l+1)}(\lambda_{i})}{G^{2(l+1)}(\lambda_{max})}\right)}}$
$\displaystyle=G(\lambda_{n})\sqrt{\frac{\alpha_{max}^{2}G^{2l}(\lambda_{max})}{\alpha_{max}^{2}G^{2(l+1)}(\lambda_{max})}}$
$\displaystyle=\frac{G(\lambda_{n})}{G(\lambda_{max})}$

∎

### Discussion on computational complexity of Spectral Module

Learning the filter coefficients $f_{c}^{s}$ at a given timestep, we can
obtain the wavelet operator $g(L)$ at scale $s=1$ using equation
$g(s_{j}x)=\sum_{k=0}^{\infty}c_{j,k}\overline{T_{k}}(x)$. The cost for
computing the terms $T_{k}(L)f$ for a sparse graph would be
$\mathcal{O}(|E|)$. There would be $M$ such terms in the polynomial function.
Also these terms could be reused across different scales. Also the computation
of the weighted sum of these terms by the coefficients at scale $s_{j}$ would
incur a complexity of $\mathcal{O}(N\times M_{j})$. Thus for $J$ scales the
total computational complexity of approximating the wavelet coefficients would
be $\mathcal{O}(|E|+N\sum_{j=0}^{J}M_{j})$. If the graph is of a bounded
degree($d_{G}$) as is the case with most real world graphs this cost would be
reduced to $\mathcal{O}(d_{G}N+N\sum_{j=0}^{J}M_{j})$, at a given snapshot of
the graph at time $t$.


11institutetext: Department of Computing, Imperial College London, London, UK
22institutetext: Data Science Institute, Imperial College London, London, UK
22email<EMAIL_ADDRESS>33institutetext: Department of Brain Sciences,
Imperial College London, London, UK 44institutetext: UK Dementia Research
Institute, Imperial College London, London, UK

# Subject-Specific Lesion Generation and Pseudo-Healthy Synthesis for Multiple
Sclerosis Brain Images

Berke Doga Basaran(✉) 1122 0000-0001-7235-3962
Mengyun Qiao 22 0000-0002-5157-1079 Paul M. Matthews 3344 0000-0002-1619-8328
Wenjia Bai 112233 0000-0003-2943-7698

###### Abstract

Understanding the intensity characteristics of brain lesions is key for
defining image-based biomarkers in neurological studies and for predicting
disease burden and outcome. In this work, we present a novel foreground-based
generative method for modelling the local lesion characteristics that can both
generate synthetic lesions on healthy images and synthesize subject-specific
pseudo-healthy images from pathological images. Furthermore, the proposed
method can be used as a data augmentation module to generate synthetic images
for training brain image segmentation networks. Experiments on multiple
sclerosis (MS) brain images acquired on magnetic resonance imaging (MRI)
demonstrate that the proposed method can generate highly realistic pseudo-
healthy and pseudo-pathological brain images. Data augmentation using the
synthetic images improves the brain image segmentation performance compared to
traditional data augmentation methods as well as a recent lesion-aware data
augmentation technique, CarveMix. The code will be released at
https://github.com/dogabasaran/lesion-synthesis.

###### Keywords:

Lesion generation Image synthesis Generative modelling Data augmentation
Attention mechanisms.

##  1 Introduction

Multiple sclerosis (MS) is a demyelinating neurodegenerative disease of the
nervous system [9]. On brain magnetic resonance (MR) images such as fluid-
attenuated inversion recovery (FLAIR), MS is often manifested as hyperintense
lesions distributed mostly in the white matter [1]. Identification and
understanding the intensity patterns of these lesions thus plays a crucial
role in diagnosing and tracking the progression of the disease. Being able to
model the lesion intensity patterns also enables us to generate synthetic
images for data augmentation purposes and to perform counter-factual inference
to understand disease progression. Although many efforts have been dedicated
in the past to the general image synthesis problem [11, 14, 18], little has
been investigated for subject-specific MS lesion synthesis.

In this work, we propose a novel generative method to model the MS lesion
characteristics. It is able to perform subject-specific lesion synthesis given
a healthy image, as well as synthesise a pseudo-healthy image given a
pathological one. An attention mechanism is designed so as to only modify the
intensities in the lesion region while maintaining the structure in other
brain regions. We demonstrate that the proposed method achieves high realism
in both lesion generation and pseudo-healthy synthesis. In addition, using the
synthetic images for training image segmentation networks, we can improve the
lesion segmentation performance even in a low-data setting.

###  1.1 Related Works

####  1.1.1 Brain Lesion Synthesis

Learning-based brain lesion synthesis has received a lot of attention recently
[22, 3]. Salem et al. uses an encoder-decoder U-Net structure to fuse lesion
masks acquired from MS patients to healthy subjects [20]. Bissoto et al., Jin
et al., and Li et al. utilise generative adversarial networks (GANs) to
synthesize new image samples [2, 13, 15]. Yet, [2] and [15] use semantic maps
from pathological subjects to generate images, which lead to augmented images
which have identical foreground shapes and labels to the samples that they are
drawn from. [13] uses a user-specified tumour mask to generate synthetic
tumours on healthy scans. These methods do not generate diverse data samples
with high foreground variance (lesion/tumour masks) which differ from already
present samples in the training dataset. Reinhold et al. employs a casual
model using variational autoencoders to generate lesions with predefined
lesion load [19]. As for subject-specific pathological to healthy image
synthesis, known as pseudo-healthy synthesis, Xia et al. uses a framework that
consists of a generator, a segmentor and a reconstructor, along with a mask
discriminator which distinguishes segmented masks from pathology masks [24].

####  1.1.2 Data Augmentation

With the assistance of data augmentation methods, deep learning has defined
the state-of-the-art for lesion segmentation [26]. Data augmentation has
improved model generalisability by increasing dataset variance, and has
alleviated some of the issues caused due to scarcity of labelled data.
However, traditional data augmentation methods, which include flipping,
rotating, scaling, intensity alterations, and elastic deformations, do not
significantly alter the original image, e.g. the appearance of lesions. Zhang
et al. proposed CarveMix, derived from CutMix [25], which utilises a lesion-
aware mix-based method to combine carved lesion patches with non-lesion images
[27]. Both of these methods directly employ foreground masks from other
subjects to compose an augmented image with already present lesion shapes and
characteristics. Current data augmentation methods for brain lesions do not
drastically change the local features of lesions, and a technique which can
realistically generate new lesions from a learnt distribution will increase
training sample diversity and segmentation performance.

###  1.2 Contributions

We present a novel brain lesion synthesis method for multiple sclerosis
images. Our contributions are four-fold: (1) Differing from other techniques,
our method is able to preserve subject identity and generate a new lesion
foreground class without directly sampling a lesion mask from another subject.
(2) An attention-based foreground discriminator is introduced to generate
realistic images focusing at the lesion region. (3) The method is able to
generate both synthetic lesions as well as pseudo-healthy images, achieved via
a cyclic structure. This enables generating a diverse set of brain images for
data augmentation purposes. (4) Our method is stochastic, and is able to
generate multiple augmented images with different lesion appearances from a
single input image.

##  2 Methods

The goal is to generate pathological images with lesions, given healthy images
as input, while preserving the subject identity, and vice versa, create
pseudo-healthy images without the lesions in the pathological subjects.
Inspired by CycleGAN [28] and AttentionGAN [23], we adopt a two-branch
architecture to learn both the lesion generation and pseudo-healthy synthesis
tasks. We employ the attention and content decoders from [23], which provide
maps on where and what the important features of the foreground and background
classes of the image are, respectively. However, we encourage the foreground
decoder to learn lesion specific features. In particular, we introduce a
foreground-based discriminator to focus on the generation of the lesion
region. We assume that two unpaired datasets are available for learning: a
healthy dataset, $X_{H}=\\{x_{H}\\}$, and a pathological dataset with lesions
or white matter hyperintensities (WMHs), $X_{P}=\\{x_{P}\\}$. Lower case
$x_{H}$ and $x_{P}$ denote single image samples from these datasets. For the
pathological dataset, a corresponding foreground dataset, $X_{F}=\\{x_{F}\\}$,
is constructed using the lesion masks. Each foreground image is a masked
lesion intensity image, generated by multiplying the original pathological
image from $X_{P}$ with its respective lesion mask.

Figure 1 illustrates the proposed framework, which consists of two branches,
one generator, $G_{P}$, for synthetic lesion generation and the other
generator, $G_{H}$, for pseudo-healthy synthesis. We represent the lesion
generation branch as $H\rightarrow P$, and the pseudo-healthy synthesis branch
as $P\rightarrow H$, where $H$ denotes healthy images, and $P$ denotes
pathological images. The lesion generator, $G_{P}$, takes in a healthy brain
image, $x_{H}$, and generates a pathological image with lesions,
$\hat{x}_{P}$, while also aiming to preserve subject identity. It is composed
of an encoder, $Enc$, which maps the input image into a latent code $z$, a
content decoder, $Dec_{C}$, that decodes $z$ into a foreground content map,
$C^{fore}$, and an attention decoder, $Dec_{A}$, that decodes $z$ into the two
attention maps, namely $A^{fore}$ for foreground and $A^{back}$ for
background. Subsequently, the foreground content, $C^{fore}$, in our case the
lesion-related regions, and attention map, $A^{fore}$, are fused to generate a
foreground output, $O^{fore}$. The input image, $x_{H}$, and background
attention, $A^{back}$, is fused to generate a background output, $O^{back}$.
The summation of $O^{fore}$ and $O^{back}$ produces the generated synthetic
image, $\hat{x}_{P}$. The pseudo-healthy synthesis branch, $G_{H}$, takes an
image with lesions, $x_{P}$, as input and produces an pseudo-healthy image,
$\hat{x}_{H}$. This branch is constructed in a similar way as the lesion
generation branch.

Figure 1: Proposed framework for subject-specific lesion generation and
pseudo-healthy synthesis of brain images, consisting of two generators,
$G_{P}$ and $G_{H}$. See text for detail.

###  2.1 Generators

The two generators, $G_{P}$ and $G_{H}$, are both constructed using the ResNet
architecture with nine blocks. Each generator features an encoder with three
downsampling paths, and two decoders, one for generating content masks, and
one for attention masks. For a given input image, our model generates $n$
attention masks, 1 for background attention, $A^{back}$, and $n-1$ for
foreground attention, $A^{fore}$, and also generates $n-1$ foreground content
masks, $C^{fore}$. Here foreground broadly refers to regions relevant to
lesions and background refers to the rest of the brain image. Producing
multiple foreground attention masks allows lesions of varying shape and
location characteristics to be decoded, thus producing synthetic images with
different lesion distributions. Similar ideas have been explored in [5], where
an input image is disentangled into multiple channels representing different
anatomical structures. The synthetic image $G(x)$ is generated as,

$G(x)=C^{fore}*A^{fore}+x*A^{back}$ (1)

which combines the foreground content, $C^{fore}$, multiplied with their
respective attention masks, $A^{fore}$, and the original input image, $x$,
multiplied with the background attention mask, $A^{back}$. The separation of
foreground and background allows us to explicitly focus on the lesions in the
generating process by developing a lesion-aware foreground discriminator.

###  2.2 Discriminators

To generate realistic brain images, three discriminators are introduced. The
pathological discriminator, $D_{P}$, distinguishes synthetic images with
lesions from real pathological images. The healthy discriminator, $D_{H}$,
distinguishes pseudo-healthy images from real healthy subjects. The third
discriminator, $D_{F}$, which is lesion-aware, encourages the foreground
(synthetic lesions) to look realistic compared to real lesion masks in the
pathological dataset.

###  2.3 Losses

The total generator loss, $\mathcal{L}_{G_{total}}$, consists of three loss
terms: generation loss, $\mathcal{L}_{G}$, cycle-consistency loss,
$\mathcal{L}_{CC}$, and identity loss, $\mathcal{L}_{idt}$. We employ the
least squares GAN loss [17] for the generation loss, also known as
$\mathcal{L}_{2}$ loss. For pseudo-healthy synthesis ($P\rightarrow H$), the
generation loss is formulated as,

$\mathcal{L}_{G_{H}}=\mathbb{E}_{x\sim
X_{P}}\left[(D_{H}(G_{H}(x))-1)^{2}\right]$ (2)

where $D_{H}$ denotes the healthy discriminator for the pseudo-healthy image
$G_{H}(x)$. For lesion generation ($H\rightarrow P$), the generation loss is
formulated as,

$\begin{split}\mathcal{L}_{G_{P}}=\mathbb{E}_{x\sim
X_{H}}\left[\left(\tfrac{1}{2}(D_{P}(G_{P}(x))+D_{F}(O^{fore}_{P}))-1\right)^{2}\right]\end{split}$
(3)

where $D_{P}$ denotes the pathological discriminator for the lesion image
$G_{P}(x)$ and $D_{F}$ denotes the foreground discriminator for the foreground
output $O^{fore}_{P}$, which penalises non-realistic lesion regions.

The cycle-consistency loss, $\mathcal{L}_{CC}$, evaluates a full cycle of an
image through the network, $P\rightarrow H\rightarrow P$ or $H\rightarrow
P\rightarrow H$, against its original input. We encourage the output of these
cycles to be equivalent to the input as we wish to preserve subject identity.
This also ensures that the generated synthetic lesions do not appear at
unrealistic locations in the brain structure. The cycle-consistency loss,
$\mathcal{L}_{CC}$, is formulated as,

$\mathcal{L}_{CC}=\mathbb{E}_{x\sim
X_{H}}\left[\|G_{H}(G_{P}(x))-x\|_{1}\lambda_{H}\right]+\mathbb{E}_{x\sim
X_{P}}\left[\|G_{P}(G_{H}(x))-x\|_{1}\lambda_{P}\right]$ (4)

Finally, we incorporate the identity loss, $\mathcal{L}_{idt}$, with the
objective of getting an identical image out of a generator if its healthiness
is not changed, such that $G_{H}(x_{H})$ would return an image identical to
$x_{H}$. We use the $\mathcal{L}_{1}$ loss for this task.

$\mathcal{L}_{idt}=\mathbb{E}_{x\sim
X_{H}}\left[\|G_{H}(x)-x\|_{1}\lambda_{H}\lambda_{idt}\right]+\mathbb{E}_{x\sim
X_{P}}\left[\|G_{P}(x)-x\|_{1}\lambda_{P}\lambda_{idt}\right]$ (5)

We utilise hyperparameters $\lambda_{H}$, $\lambda_{P}$, and $\lambda_{idt}$
for weighting these losses. The total loss for the generator is summarized as,

$\mathcal{L}_{G_{total}}=\mathcal{L}_{G_{H}}+\mathcal{L}_{G_{P}}+\mathcal{L}_{CC}+\mathcal{L}_{idt}$
(6)

The discriminators aim to minimize the sum of the squared difference between
predicted and expected values for real and synthetic images. We use the least
squares loss for the three discriminators, $D_{H}$, $D_{P}$, $D_{F}$,
formulated as,

$\mathcal{L}_{D_{H}}=\mathbb{E}\left[(D_{H}(x_{H})-1)^{2}+(D_{H}(G_{H}(x_{P})))^{2}\right]$
(7)
$\mathcal{L}_{D_{P}}=\mathbb{E}\left[(D_{P}(x_{P})-1)^{2}+(D_{P}(G_{P}(x_{H})))^{2}\right]$
(8)
$\mathcal{L}_{D_{F}}=\mathbb{E}\left[(D_{F}(x_{F})-1)^{2}+(D_{F}(O^{fore}_{P}))^{2}\right]$
(9)

##  3 Experiments

###  3.1 Evaluation

We evaluate the method performance in data augmentation using the synthetic
data for downstream segmentation task as well as evaluate the image quality of
synthetic data. To assess performance we conduct two experiments, employing
our method as a data augmentation tool and comparing segmentation performance
when trained on only synthetic data. To evaluate quality we ask three raters
to assess healthiness and realism of generated images.

The proposed method is compared with traditional data augmentation (TDA)
techniques which come default with nnU-Net [12], including rotation, scaling,
mirroring, elastic deformation, intensity perturbation, and simulation of low
resolution. We also compare to a state-of-the-art lesion-aware augmentation
method, CarveMix [27]. We utilise offline versions of our method and CarveMix
to double the training dataset for those experiments. When performing data
augmentation using the proposed method, lesion masks are required for the
foreground discriminator. After lesion generation, we perform free-form
deformation registration of the healthy subject image, $x_{H}$, to the
corresponding synthesized pathological image, $\hat{x}_{P}$, threshold the
difference of the two images and generate the lesion masks. We report the mean
and variance of the Dice scores and Hausdorff distances of the segmentations.

###  3.2 Implementation

####  3.2.1 GAN Network

We implement the generative model with a single channel input and single
channel output. We set the generated mask number $n$ to 10, as our experiments
show more variety of realistic lesion characteristics and intensities are
synthesised with a higher $n$. Our generative framework method also features
data augmentation in the form of sagittal mirroring and elastic deformation.

####  3.2.2 Training

The proposed method is developed on PyTorch and trained on one NVIDIA GeForce
RTX 3080. We use the Adam optimizer with an initial learning rate of 0.001 and
0.5 dropout. The network is trained for 400 epochs and with a linear decay of
the learning rate after 200 epochs with a batch size of 1, due to GPU memory
limit. We set the generator loss trade-off hyperparameters $\lambda_{H}$,
$\lambda_{P}$, and $\lambda_{idt}$ to $10$, $10$, and $0.5$, respectively. For
evaluating segmentation performance, nnU-Net is trained for 1,000 epochs using
stochastic gradient descent.

###  3.3 Data

A private dataset is used for training the proposed method, which consists of
FLAIR MRI scans of 120 subjects, including 60 healthy scans and 60 scans with
MS lesions. For evaluation, we train the segmentation networks using the 2016
Multiple Sclerosis Lesion Segmentation dataset (MS2016) [6], and test on the
2008 MICCAI MS Lesion Segmentation dataset (MS2008) [21] and 2015 ISBI
Longitudinal MS Lesion Segmentation dataset (ISBI2015) [4]. We resample all
images into $1\times 1\times 1$ mm3 voxel spacing and extract the centre axial
slice of each image. These images are reshaped into $256\times 256$ dimensions
to pass through the network. The private, MS2008, ISBI2015, and MS2016
datasets provide 120, 20, 21, and 15 2D FLAIR images, respectively.

###  3.4 Results

####  3.4.1 Data Augmentation Performance

We train twelve separate segmentation models using different dataset sizes and
data augmentation methods. We use nnU-Net in 3D full resolution, utilizing
100%, 53.3%, 26.6% and 12.3% of the training dataset, corresponding to 15, 8,
4, and 2 subjects.

Table 1: Mean and standard deviations of the Dice scores ($\%$), at different sizes of training data. Best results are in bold. Asterisks indicate statistical significance (∗: p$\leq$ 0.05, ∗∗: p $\leq$ 0.001) when using a paired Student’s t-test compared to baselines. Size | TDA | TDA+CarveMix[27] | TDA+Proposed
---|---|---|---
MS2008 | ISBI2015 | MS2008 | ISBI2015 | MS2008 | ISBI2015
$100$% | $33.06^{**}_{17.37}$ | $63.20_{16.70}$ | $32.32^{**}_{15.23}$ | $57.82_{16.56}$ | $\textbf{44.77}_{\textbf{15.87}}$ | $\textbf{63.90}_{\textbf{16.45}}$
$53.3$% | $34.03_{17.83}$ | $58.84_{24.54}$ | $35.39_{17.36}$ | $61.25_{19.06}$ | $\textbf{39.70}_{\textbf{18.44}}$ | $\textbf{62.35}_{\textbf{13.21}}$
$26.6$% | $31.55_{18.62}$ | $\textbf{64.43}_{\textbf{12.64}}$ | $\textbf{39.21}_{\textbf{18.86}}$ | $60.28_{21.22}$ | $34.58_{15.47}$ | $60.83_{10.03}$
$13.3$% | $28.96_{23.55}$ | $37.34^{**}_{23.81}$ | $23.57{*}_{23.18}$ | $38.94^{**}_{21.81}$ | $\textbf{32.29}_{\textbf{16.16}}$ | $\textbf{56.64}_{\textbf{13.23}}$
Table 2: Mean and standard deviations of the Hausdorff distances, at different sizes of training data. Best results are in bold. Asterisks indicate statistical significance (∗: p$\leq$ 0.05, ∗∗: p $\leq$ 0.01, ∗∗∗: p $\leq$ 0.001) when using a paired Student’s t-test compared to baselines. Size | TDA | TDA+CarveMix [27] | TDA+Proposed
---|---|---|---
MS2008 | ISBI2015 | MS2008 | ISBI2015 | MS2008 | ISBI2015
$100$% | $3.52^{*}_{1.61}$ | $3.35^{*}_{0.94}$ | $3.24_{1.09}$ | $3.33_{0.89}$ | $\textbf{3.09}_{\textbf{{1.54}}}$ | $\textbf{3.02}_{\textbf{0.80}}$
$53.3$% | $3.16^{**}_{1.10}$ | $3.12_{0.93}$ | $3.30^{***}_{1.26}$ | $3.28_{1.00}$ | $\textbf{2.62}_{\textbf{0.98}}$ | $\textbf{3.07}_{\textbf{1.08}}$
$26.6$% | $3.12_{0.89}$ | $\textbf{3.38}_{\textbf{1.01}}$ | $\textbf{{3.08}}_{\textbf{1.01}}$ | $3.37_{0.94}$ | $3.09_{0.72}$ | $3.86_{1.61}$
$13.3$% | $2.92_{1.13}$ | $3.86^{**}_{1.62}$ | $2.95_{1.03}$ | $3.81^{**}_{1.62}$ | $\textbf{2.90}_{\textbf{0.69}}$ | $\textbf{3.33}_{\textbf{1.17}}$

Table 1 indicates that the average Dice score is increased using the proposed
method in three of the four dataset sizes when compared to TDA and CarveMix,
while also decreasing the standard deviation in most cases. We further report
the mean and variance of the Hausdorff distances in Table 2. Figure 2
visualises the segmentation performance of our method against the benchmark.
Figure 3 demonstrates different stages of the proposed method and augmentation
examples from CarveMix. [27].

Figure 2: Comparison of segmentation performance using different data
augmentation techniques. Lesions are contoured in red. Left to right: Test
subjects with ground truth lesion masks; segmentation using traditional data
augmentation (TDA); segmentation using TDA+CarveMix [27]; segmentation using
TDA+proposed. Figure 3: Different stages of the proposed method. Lesions are
contoured in red. Left to right: Real pathological image, $x_{P}$; pseudo-
healthy image $\hat{x}_{H}$; foreground output for lesion synthesis on pseudo-
healthy image $O^{fore}_{P}$; background output for lesion synthesis
$O^{back}_{P}$; the synthetic pathological image generated from the pseudo-
healthy image $\hat{x}_{P}$; augmented image using CarveMix, unrealistic
lesions circled in yellow.

The generator of the method is stochastic, which employs dropout to generate
multiple diverse augmented images for each given sample, shown in Figure 4.
This increases the diversity of training image samples.

Figure 4: Setting the generative network’s dropout level to 0.5 at training
and inference allows for generating a diverse set of synthetic images (second
to fourth columns) given a single image as input (first column).

####  3.4.2 Synthetic data segmentation performance

We compare segmentation performance when training a model using a fully
synthetic dataset against a model trained on a real dataset. We generate 15
synthetic images from the MS2016 dataset and test on the MS2008 and ISBI2015
datasets. We compare our results to nnU-Net with TDA at 100% dataset size, as
per in Table 1. Results in Table 3 suggest a model trained on our generated
data performs comparable to a model trained on a real dataset.

Table 3: Segmentation performance using fully synthetic data for model training, compared to using real data. Mean and standard deviations of the Dice scores (%) are reported with best results in bold. Test data | Real training data | Synthetic training data
---|---|---
MS2008 | $33.06_{17.37}$ | $\textbf{35.74}_{\textbf{19.85}}$
ISBI2015 | $\textbf{63.20}_{\textbf{16.70}}$ | $52.72_{18.61}$

####  3.4.3 Expert Evaluation

In order to assess quality of generated images, we ask three human raters to
evaluate images in healthiness and realism. We define "healthiness" as a brain
MR image with no visible pathological features, and "realism" as an image
which possesses the same quality and authentic look as one obtained by a
scanner. We collate 50 real healthy, 50 real pathological, 50 pseudo-healthy,
and 50 pseudo-pathological images. Raters are asked to classify each image as
healthy (1) or pathological (0), and real (1) or fake (0). We produce the
results for the real images as a benchmark to compare our generated images.
Table 4 shows that the generated pseudo-healthy and pseudo-pathological images
closely follow real images for the assessed metrics.

Table 4: Rater classification of real and generated images on realism and
healthiness. Asterisk (*) indicates that a lower score is better.

Realism
---
|  Real | Proposed
Healthy images |  85% | 82%
Pathological images |  84% | 82%

Healthiness
---
|  Real | Proposed
Healthy images |  85% | 83%
Pathological images* |  5% | 1%

####  3.4.4 Correlation with Radiological Findings

The spatial heatmap of the synthetic lesions closely resemble those of the
real lesions from clinical data, including the private dataset and the MS2016
dataset. A comparison between the spatial heatmap of synthetic lesions and
real lesions is visualised in Figure 5.

Figure 5: Comparison of spatial heatmaps of real lesions and synthetic
lesions. A: Real lesions from private dataset (200 subjects). B: Real lesions
in MS2016 dataset (60 subjects). C: Synthetic lesions generated using the
proposed method by training on the MS2016 dataset (200 masks).

Furthermore, generated pseudo-pathological images show an increase in the
ventricular volume, while pseudo-healthy images tend to have a decreased
ventricular volume. The difference in ventricular volume is shown in Figure 6.
The change in the ventricular volume is consistent with literature [8, 7, 10,
16].

Figure 6: Apart from generating lesions or pseudo-healthy synthesis, the
proposed method also changes the ventricular volume during image generation,
which is consistent with clinical literature. Top: Pseudo-healthy synthesis
leads to decrease in ventricular volume, delineated in red in the difference
image. Bottom: Lesion generation leads to an increase in ventricular volume,
delineated in red in the difference image.

##  4 Conclusion

We have proposed a novel method for generating lesions and pseudo-healthy
images while preserving subject-specific features, providing a useful tool for
data augmentation in brain image analysis to complement real patient datasets.
It is lesion-aware in the generating process. In quantitative evaluation, the
method improves lesion segmentation performance in downstream tasks. In our
human rater assessment, the generated images achieved a high realism score
close to the real images, and generated lesions follow a spatial distribution
consistent with real clinical datasets.

####  4.0.1 Acknowledgements

This work is supported by the UKRI CDT in AI for Healthcare
http://ai4health.io (Grant No. EP/S023283/1).

## References

  * [1] Bakshi, R., Ariyaratana, S., Benedict, R.H.B., Jacobs, L.: Fluid-Attenuated Inversion Recovery Magnetic Resonance Imaging Detects Cortical and Juxtacortical Multiple Sclerosis Lesions. Archives of Neurology 58(5), 742–748 (2001). https://doi.org/10.1001/archneur.58.5.742
  * [2] Bissoto, A., Perez, F., Valle, E., Avila, S.: Skin Lesion Synthesis with Generative Adversarial Networks. In: OR 2.0 Context-Aware Operating Theaters, Computer Assisted Robotic Endoscopy, Clinical Image-Based Procedures, and Skin Image Analysis. pp. 294–302. Springer International Publishing, Cham (2018). https://doi.org/10.1007/978-3-030-01201-4\\_32
  * [3] Bowles, C., Qin, C., Guerrero, R., et al.: Brain lesion segmentation through image synthesis and outlier detection. NeuroImage: Clinical 16, 643–658 (2017). https://doi.org/10.1016/j.nicl.2017.09.003
  * [4] Carass, A., Roy, S., Jog, A., et al.: Longitudinal multiple sclerosis lesion segmentation: Resource and challenge. NeuroImage 148, 77–102 (2017). https://doi.org/10.1016/j.neuroimage.2016.12.064
  * [5] Chartsias, A., Joyce, T., Papanastasiou, G., et al.: Disentangled representation learning in cardiac image analysis. Medical Image Analysis 58, 101535 (2019). https://doi.org/10.1016/j.media.2019.101535
  * [6] Commowick, O., Istace, A., Kain, M., et al.: Objective evaluation of multiple sclerosis lesion segmentation using a data management and processing infrastructure. Scientific Reports 8(1), 13650 (Sep 2018). https://doi.org/10.1038/s41598-018-31911-7
  * [7] Dalton, C.M., Brex, P.A., Jenkins, R., et al.: Progressive ventricular enlargement in patients with clinically isolated syndromes is associated with the early development of multiple sclerosis. Journal of Neurology, Neurosurgery & Psychiatry 73(2), 141–147 (2002). https://doi.org/10.1136/jnnp.73.2.141
  * [8] Dalton, C.M., Miszkiel, K.A., O‘Connor, P.W., et al.: Ventricular enlargement in MS. Neurology 66(5), 693–698 (2006). https://doi.org/10.1212/01.wnl.0000201183.87175.9f
  * [9] Ghasemi, N., Razavi, S., Nikzad, E.: Multiple Sclerosis: Pathogenesis, Symptoms, Diagnoses and Cell-Based Therapy Citation: Ghasemi N, Razavi Sh, Nikzad E. Multiple sclerosis: pathogenesis, symptoms, diagnoses and cell-based therapy. Cell Journal 19(191), 1–10 (2017). https://doi.org/10.22074/cellj.2016.4867
  * [10] Guptha, S.H., Holroyd, E., Campbell, G.: Progressive lateral ventricular enlargement as a clue to Alzheimer’s disease. The Lancet 359(9322), 2040 (2002). https://doi.org/10.1016/S0140-6736(02)08806-2
  * [11] Huang, H., Yu, P.S., Wang, C.: An introduction to image synthesis with generative adversarial nets (2018). https://doi.org/10.48550/ARXIV.1803.04469
  * [12] Isensee, F., Jaeger, P.F., Kohl, S.A.A., Petersen, J., Maier-Hein, K.H.: nnu-net: a self-configuring method for deep learning-based biomedical image segmentation. Nature Methods 18(2), 203–211 (Feb 2021). https://doi.org/10.1038/s41592-020-01008-z
  * [13] Jin, Q., Cui, H., Sun, C., Meng, Z., Su, R.: Free-form tumor synthesis in computed tomography images via richer generative adversarial network. Knowledge-Based Systems 218, 106753 (2021). https://doi.org/10.1016/j.knosys.2021.106753
  * [14] Karras, T., Laine, S., Aila, T.: A style-based generator architecture for generative adversarial networks. In: 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). pp. 4396–4405 (2019). https://doi.org/10.1109/CVPR.2019.00453
  * [15] Li, Q., Yu, Z., Wang, Y., Zheng, H.: Tumorgan: A multi-modal data augmentation framework for brain tumor segmentation. Sensors 20(15) (2020). https://doi.org/10.3390/s20154203
  * [16] Luxenberg, J.S., Haxby, J.V., Creasey, H., Sundaram, M., Rapoport, S.I.: Rate of ventricular enlargement in dementia of the Alzheimer type correlates with rate of neuropsychological deterioration. Neurology 37(7), 1135 (1987). https://doi.org/10.1212/WNL.37.7.1135
  * [17] Mao, X., Li, Q., Xie, H., et al.: Least squares generative adversarial networks. In: 2017 IEEE International Conference on Computer Vision (ICCV). pp. 2813–2821 (2017). https://doi.org/10.1109/ICCV.2017.304
  * [18] Reed, S., Akata, Z., Yan, X., et al.: Generative adversarial text to image synthesis. In: Proceedings of The 33rd International Conference on Machine Learning. Proceedings of Machine Learning Research, vol. 48, pp. 1060–1069. PMLR, New York, New York, USA (20–22 Jun 2016)
  * [19] Reinhold, J.C., Carass, A., Prince, J.L.: A Structural Causal Model for MR Images of Multiple Sclerosis. In: Medical Image Computing and Computer Assisted Intervention – MICCAI 2021. pp. 782–792. Springer International Publishing, Cham (2021). https://doi.org/10.1007/978-3-030-87240-3\\_75
  * [20] Salem, M., Valverde, S., Cabezas, M., et al.: Multiple Sclerosis Lesion Synthesis in MRI Using an Encoder-Decoder U-NET. IEEE Access 7, 25171–25184 (2019). https://doi.org/10.1109/ACCESS.2019.2900198
  * [21] Styner, M., Lee, J., Chin, B., et al.: 3d segmentation in the clinic: A grand challenge ii: Ms lesion segmentation (11 2008). https://doi.org/10.54294/lmkqvm
  * [22] Sun, L., Wang, J., Huang, Y., et al.: An adversarial learning approach to medical image synthesis for lesion detection. IEEE Journal of Biomedical and Health Informatics 24(8), 2303–2314 (2020). https://doi.org/10.1109/JBHI.2020.2964016
  * [23] Tang, H., Xu, D., Sebe, N., Yan, Y.: Attention-guided generative adversarial networks for unsupervised image-to-image translation. In: 2019 International Joint Conference on Neural Networks (IJCNN). pp. 1–8 (2019). https://doi.org/10.1109/IJCNN.2019.8851881
  * [24] Xia, T., Chartsias, A., Tsaftaris, S.A.: Pseudo-healthy synthesis with pathology disentanglement and adversarial learning. Medical Image Analysis 64, 101719 (2020). https://doi.org/10.1016/j.media.2020.101719
  * [25] Yun, S., Han, D., Chun, S., et al.: Cutmix: Regularization strategy to train strong classifiers with localizable features. In: 2019 IEEE/CVF International Conference on Computer Vision (ICCV). pp. 6022–6031 (2019). https://doi.org/10.1109/ICCV.2019.00612
  * [26] Zeng, C., Gu, L., Liu, Z., Zhao, S.: Review of Deep Learning Approaches for the Segmentation of Multiple Sclerosis Lesions on Brain MRI. Frontiers in Neuroinformatics 14 (2020). https://doi.org/10.3389/fninf.2020.610967
  * [27] Zhang, X., Liu, C., Ou, N., et al.: CarveMix: A Simple Data Augmentation Method for Brain Lesion Segmentation. In: Medical Image Computing and Computer Assisted Intervention – MICCAI 2021. pp. 196–205. Springer International Publishing, Cham (2021)
  * [28] Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired image-to-image translation using cycle-consistent adversarial networks. In: 2017 IEEE International Conference on Computer Vision (ICCV). pp. 2242–2251 (2017). https://doi.org/10.1109/ICCV.2017.244