# Fast Temporal Wavelet Graph Neural Networks

Duc Thien Nguyen  Manh Duc Tuan Nguyen  Truong Son Hy  Risi Kondor

###### Abstract

Spatio-temporal signals forecasting plays an important role in numerous
domains, especially in neuroscience and transportation. The task is
challenging due to the highly intricate spatial structure, as well as the non-
linear temporal dynamics of the network. To facilitate reliable and timely
forecast for the human brain and traffic networks, we propose the Fast
Temporal Wavelet Graph Neural Networks (FTWGNN) that is both time- and memory-
efficient for learning tasks on timeseries data with the underlying graph
structure, thanks to the theories of multiresolution analysis and wavelet
theory on discrete spaces. We employ Multiresolution Matrix Factorization
(MMF) (Kondor et al., 2014) to factorize the highly dense graph structure and
compute the corresponding sparse wavelet basis that allows us to construct
fast wavelet convolution as the backbone of our novel architecture.
Experimental results on real-world PEMS-BAY, METR-LA traffic datasets and
AJILE12 ECoG dataset show that FTWGNN is competitive with the state-of-the-
arts while maintaining a low computational footprint. Our PyTorch
implementation is publicly available at https://github.com/HySonLab/TWGNN.

Multiresolution matrix factorization, multiresolution analysis, wavelet
theory, graph wavelets, wavelet convolution, spatio-temporal architecture.

##  1 Introduction

Time series modeling has been a quest in a wide range of academic fields and
industrial applications, including neuroscience (Pourahmadi & Noorbaloochi,
2016) and traffic modeling (Li et al., 2018). Traditionally, model-based
approaches such as autoregressive (AR) and Support Vector Regression (Smola &
Schölkopf, 2004) require domain-knowledge as well as stationary assumption,
which are often violated by the complex and non-linear structure of neural and
traffic data.

Recently, there has been intensive research with promising results on the
traffic forecasting problem using deep learning such as Recurrent Neural
Network (RNN) (Qin et al., 2017), LSTM (Koprinska et al., 2018), and graph-
learning using Tranformer (Xu et al., 2020). On the other hand, forcasting in
neuroscience has been focusing mainly on long-term evolution of brain network
structure based on fMRI data, such as predicting brain connectivities of an
Alzheimer’s disease after several months (Bessadok et al., 2022), where
existing methods are GCN-based (Göktaş et al., 2020) or GAN-based graph
autoencoder (Gürler et al., 2020). Meanwhile, research on instantaneous time
series forecasting of electroencephalogram (EEG) or electrocorticography
(ECoG) remains untouched, even though EEG and ECoG are often cheaper and
quicker to obtain than fMRI, while short-term forecasting may be beneficial
for patients with strokes or epilepsy (Shoeibi et al., 2022).

In graph representation learning, a dense adjacency matrix expressing a
densely connected graph can be a waste of computational resources, while
physically, it may fail to capture the local “smoothness” of the network. To
tackle such problems, a mathematical framework called Multiresolution Matrix
Factorization (MMF) (Kondor et al., 2014) has been adopted to “sparsify” the
adjacency and graph Laplacian matrices of highly dense graphs. MMF is unusual
amongst fast matrix factorization algorithms in that it does not make a low
rank assumption. Multiresolution matrix factorization (MMF) is an alternative
paradigm that is designed to capture structure at multiple different scales.
This makes MMF especially well suited to modeling certain types of graphs with
complex multiscale or hierarchical strucutre (Hy & Kondor, 2022), compressing
hierarchical matrices (e.g., kernel/gram matrices) (Teneva et al., 2016; Ding
et al., 2017), and other applications in computer vision (Ithapu et al.,
2017). One important aspect of MMF is its ability to construct wavelets on
graphs and matrices during the factorization process (Kondor et al., 2014; Hy
& Kondor, 2022). The wavelet basis inferred by MMF tends to be highly sparse,
that allows the corresponding wavelet transform to be executed efficiently via
sparse matrix multiplication. (Hy & Kondor, 2022) exploited this property to
construct fast wavelet convolution and consequentially wavelet neural networks
learning on graphs for graph classification and node classification tasks. In
this work, we propose the incorporation of fast wavelet convolution based on
MMF to build a time- and memory-efficient temporal architecture learning on
timeseries data with the underlying graph structure.

From the aforementioned arguments, we propose the Fast Temporal Wavelet Graph
Neural Network (FTWGNN) for graph time series forecasting, in which the MMF
theory is utilized to describe the local smoothness of the network as well as
to accelerate the calculations. Experiments on real-world traffic and ECoG
datasets show competitive performance along with remarkably smaller
computational footprint of FTWGNN. In summary:

  * •

We model the spatial domain of the graph time series as a diffusion process,
in which the theories of multiresolution analysis and wavelet theory are
adopted. We employ Multiresolution Matrix Factorization (MMF) to factorize the
underlying graph structure and derive its sparse wavelet basis.

  * •

We propose the Fast Temporal Wavelet Graph Neural Network (FTWGNN), an end-to-
end model capable of modeling spatiotemporal structures.

  * •

We tested on two real-world traffic datasets and an ECoG dataset and achieved
competitive results to state-of-the-art methods with remarkable reduction in
computational time.

##  2 Related work

A spatial-temporal forecasting task utilizes spatial-temporal data information
gathered from various sensors to predict their future states. Traditional
approaches, such as the autoregressive integrated moving average (ARIMA),
k-nearest neighbors algorithm (kNN), and support vector machine (SVM), can
only take into account temporal information without considering spatial
features (Van Lint & Van Hinsbergen, 2012; Jeong et al., 2013). Aside from
traditional approaches, deep neural networks are proposed to model much more
complex spatial-temporal relationships. Specifically, by using an extended
fully-connected LSTM with embedded convolutional layers, FC-LSTM (Sutskever et
al., 2014) specifically combines CNN and LSTM to model spatial and temporal
relations. When predicting traffic, ST-ResNet (Zhang et al., 2017) uses a deep
residual CNN network, revealing the powerful capabilities of the residual
network. Despite the impressive results obtained, traffic forecastin scenarios
with graph-structured data is incompatible with all of the aforementioned
methods because they are built for grid data. For learning tasks on graphs,
node representations in GNNs (Kipf & Welling, 2016) uses a neighborhood
aggregation scheme, which involves sampling and aggregating the features of
nearby nodes. Since temporal-spatial data such as traffic data or brain
network is a well-known type of non-Euclidean structured graph data, great
efforts have been made to use graph convolution methods in traffic
forecasting. As an illustration, DCRNN (Li et al., 2018) models traffic flow
as a diffusion process and uses directed graph bidirectional random walks to
model spatial dependency.

In the field of image and signal processing, processing is more efficient and
simpler in a sparse representation where fewer coefficients reveal the
information that we are searching for. Based on this motivation,
Multiresolution Analysis (MRA) has been proposed by (Mallat, 1989) as a design
for multiscale signal approximation in which the sparse representations can be
constructed by decomposing signals over elementary waveforms chosen in a
family called wavelets. Besides Fourier transforms, the discovery of wavelet
orthogonal bases such as Haar (Haar, 1910) and Daubechies (Daubechies, 1988)
has opened the door to new transforms such as continuous and discrete wavelet
transforms and the fast wavelet transform algorithm that have become crucial
for several computer applications (Mallat, 2008).

(Kondor et al., 2014) and (Hy & Kondor, 2022) have introduced Multiresolution
Matrix Factorization (MMF) as a novel method for constructing sparse wavelet
transforms of functions defined on the nodes of an arbitrary graph while
giving a multiresolution approximation of hierarchical matrices. MMF is
closely related to other works on constructing wavelet bases on discrete
spaces, including wavelets defined based on diagonalizing the diffusion
operator or the normalized graph Laplacian (Coifman & Maggioni, 2006) (Hammond
et al., 2011) and multiresolution on trees (Gavish et al., 2010) (Bickel &
Ritov, 2008).

##  3 Background

###  3.1 Multiresolution Matrix Factorization

Most commonly used matrix factorization algorithms, such as principal
component analysis (PCA), singular value decomposition (SVD), or non-negative
matrix factorization (NMF) are inherently single-level algorithms. Saying that
a symmetric matrix ${\bm{A}}\in\mathbb{R}^{n\times n}$ is of rank $r\ll n$
means that it can be expressed in terms of a dictionary of $r$ mutually
orthogonal unit vectors $\\{u_{1},u_{2},\dots,u_{r}\\}$ in the form

${\bm{A}}=\sum_{i=1}^{r}\lambda_{i}u_{i}u_{i}^{T},$

where $u_{1},\dots,u_{r}$ are the normalized eigenvectors of $A$ and
$\lambda_{1},\dots,\lambda_{r}$ are the corresponding eigenvalues. This is the
decomposition that PCA finds, and it corresponds to factorizing ${\bm{A}}$ in
the form

${\bm{A}}={\bm{U}}^{T}{\bm{H}}{\bm{U}},$ (1)

where ${\bm{U}}$ is an orthogonal matrix and ${\bm{H}}$ is a diagonal matrix
with the eigenvalues of ${\bm{A}}$ on its diagonal. The drawback of PCA is
that eigenvectors are almost always dense, while matrices occuring in learning
problems, especially those related to graphs, often have strong locality
properties, in the sense that they are more closely couple certain clusters of
nearby coordinates than those farther apart with respect to the underlying
topology. In such cases, modeling $A$ in terms of a basis of global
eigenfunctions is both computationally wasteful and conceptually unreasonable:
a localized dictionary would be more appropriate. In contrast to PCA, (Kondor
et al., 2014) proposed Multiresolution Matrix Factorization, or MMF for short,
to construct a sparse hierarchical system of $L$-level dictionaries. The
corresponding matrix factorization is of the form

${\bm{A}}={\bm{U}}_{1}^{T}{\bm{U}}_{2}^{T}\dots\bm{U}_{L}^{T}{\bm{H}}{\bm{U}}_{L}\dots\bm{U}_{2}{\bm{U}}_{1},$

where ${\bm{H}}$ is close to diagonal and ${\bm{U}}_{1},\dots,{\bm{U}}_{L}$
are sparse orthogonal matrices with the following constraints:

  * •

Each ${\bm{U}}_{\ell}$ is $k$-point rotation (i.e. Givens rotation) for some
small $k$, meaning that it only rotates $k$ coordinates at a time. Formally,
Def. 3.1 defines the $k$-point rotation matrix.

  * •

There is a nested sequence of sets
${\mathbb{S}}_{L}\subseteq\cdots\subseteq{\mathbb{S}}_{1}\subseteq{\mathbb{S}}_{0}=[n]$
such that the coordinates rotated by ${\bm{U}}_{\ell}$ are a subset of
${\mathbb{S}}_{\ell}$.

  * •

${\bm{H}}$ is an ${\mathbb{S}}_{L}$-core-diagonal matrix that is formally
defined in Def. 3.2.

We formally define MMF in Defs. 3.3 and 3.4. A special case of MMF is the
Jacobi eigenvalue algorithm (Jacobi, 1846) in which each ${\bm{U}}_{\ell}$ is
a 2-point rotation (i.e. $k=2$).

######  Definition 3.1.

We say that ${\bm{U}}\in\mathbb{R}^{n\times n}$ is an elementary rotation of
order $k$ (also called as a $k$-point rotation) if it is an orthogonal matrix
of the form

${\bm{U}}={\bm{I}}_{n-k}\oplus_{(i_{1},\cdots,i_{k})}{\bm{O}}$

for some ${\mathbb{I}}=\\{i_{1},\cdots,i_{k}\\}\subseteq[n]$ and
${\bm{O}}\in{\mathbb{S}}{\mathbb{O}}(k)$. We denote the set of all such
matrices as ${\mathbb{S}}{\mathbb{O}}_{k}(n)$.

######  Definition 3.2.

Given a set ${\mathbb{S}}\subseteq[n]$, we say that a matrix
${\bm{H}}\in\mathbb{R}^{n\times n}$ is ${\mathbb{S}}$-core-diagonal if
${\bm{H}}_{i,j}=0$ unless $i,j\in{\mathbb{S}}$ or $i=j$. Equivalently,
${\bm{H}}$ is ${\mathbb{S}}$-core-diagonal if it can be written in the form
${\bm{H}}={\bm{D}}\oplus_{{\mathbb{S}}}\overline{{\bm{H}}}$, for some
$\overline{H}\in\mathbb{R}^{\lvert{\mathbb{S}}\rvert\times\lvert{\mathbb{S}}\rvert}$
and ${\bm{D}}$ is diagonal. We denote the set of all ${\mathbb{S}}$-core-
diagonal symmetric matrices of dimension $n$ as
${\mathbb{H}}^{{\mathbb{S}}}_{n}$.

######  Definition 3.3.

Given an appropriate subset ${\mathbb{O}}$ of the group
${\mathbb{S}}{\mathbb{O}}(n)$ of $n$-dimensional rotation matrices, a depth
parameter $L\in\mathbb{N}$, and a sequence of integers $n=d_{0}\geq d_{1}\geq
d_{2}\geq\dots\geq d_{L}\geq 1$, a Multiresolution Matrix Factorization (MMF)
of a symmetric matrix ${\bm{A}}\in\mathbb{R}^{n\times n}$ over ${\mathbb{O}}$
is a factorization of the form

${\bm{A}}={\bm{U}}_{1}^{T}{\bm{U}}_{2}^{T}\dots\bm{U}_{L}^{T}{\bm{H}}{\bm{U}}_{L}\dots\bm{U}_{2}{\bm{U}}_{1},$
(2)

where each ${\bm{U}}_{\ell}\in{\mathbb{O}}$ satisfies
$[{\bm{U}}_{\ell}]_{[n]\setminus{\mathbb{S}}_{\ell-1},[n]\setminus{\mathbb{S}}_{\ell-1}}={\bm{I}}_{n-d_{\ell}}$
for some nested sequence of sets
${\mathbb{S}}_{L}\subseteq\cdots\subseteq{\mathbb{S}}_{1}\subseteq{\mathbb{S}}_{0}=[n]$
with $\lvert{\mathbb{S}}_{\ell}\rvert=d_{\ell}$, and
${\bm{H}}\in{\mathbb{H}}^{{\mathbb{S}}_{L}}_{n}$ is an
${\mathbb{S}}_{L}$-core-diagonal matrix.

######  Definition 3.4.

We say that a symmetric matrix ${\bm{A}}\in\mathbb{R}^{n\times n}$ is fully
multiresolution factorizable over
${\mathbb{O}}\subset{\mathbb{S}}{\mathbb{O}}(n)$ with $(d_{1},\dots,d_{L})$ if
it has a decomposition of the form described in Def. 3.3.

###  3.2 Multiresolution analysis

(Kondor et al., 2014) has shown that MMF mirrors the classical theory of
multiresolution analysis (MRA) on the real line (Mallat, 1989) to discrete
spaces. The functional analytic view of wavelets is provided by MRA, which,
similarly to Fourier analysis, is a way of filtering some function space into
a sequence of subspaces

$\dots\subset{\mathbb{V}}_{-1}\subset{\mathbb{V}}_{0}\subset{\mathbb{V}}_{1}\subset{\mathbb{V}}_{2}\subset\dots$
(3)

However, it is best to conceptualize (3) as an iterative process of splitting
each ${\mathbb{V}}_{\ell}$ into the orthogonal sum
${\mathbb{V}}_{\ell}={\mathbb{V}}_{\ell+1}\oplus{\mathbb{W}}_{\ell+1}$ of a
smoother part ${\mathbb{V}}_{\ell+1}$, called the approximation space; and a
rougher part ${\mathbb{W}}_{\ell+1}$, called the detail space (see Fig. 1).
Each ${\mathbb{V}}_{\ell}$ has an orthonormal basis
$\Phi_{\ell}\triangleq\\{\phi_{m}^{\ell}\\}_{m}$ in which each $\phi$ is
called a father wavelet. Each complementary space ${\mathbb{W}}_{\ell}$ is
also spanned by an orthonormal basis
$\Psi_{\ell}\triangleq\\{\psi_{m}^{\ell}\\}_{m}$ in which each $\psi$ is
called a mother wavelet. In MMF, each individual rotation
${\bm{U}}_{\ell}:{\mathbb{V}}_{\ell-1}\rightarrow{\mathbb{V}}_{\ell}\oplus{\mathbb{W}}_{\ell}$
is a sparse basis transform that expresses $\Phi_{\ell}\cup\Psi_{\ell}$ in the
previous basis $\Phi_{\ell-1}$ such that:

$\phi_{m}^{\ell}=\sum_{i=1}^{\text{dim}({\mathbb{V}}_{\ell-1})}[{\bm{U}}_{\ell}]_{m,i}\phi_{i}^{\ell-1},$
$\psi_{m}^{\ell}=\sum_{i=1}^{\text{dim}({\mathbb{V}}_{\ell-1})}[{\bm{U}}_{\ell}]_{m+\text{dim}({\mathbb{V}}_{\ell-1}),i}\phi_{i}^{\ell-1},$

in which $\Phi_{0}$ is the standard basis, i.e. $\phi_{m}^{0}=e_{m}$; and
$\text{dim}({\mathbb{V}}_{\ell})=d_{\ell}=\lvert{\mathbb{S}}_{\ell}\rvert$. In
the $\Phi_{1}\cup\Psi_{1}$ basis, ${\bm{A}}$ compresses into
${\bm{A}}_{1}={\bm{U}}_{1}{\bm{A}}{\bm{U}}_{1}^{T}$. In the
$\Phi_{2}\cup\Psi_{2}\cup\Psi_{1}$ basis, it becomes
${\bm{A}}_{2}={\bm{U}}_{2}{\bm{U}}_{1}{\bm{A}}{\bm{U}}_{1}^{T}{\bm{U}}_{2}^{T}$,
and so on. Finally, in the $\Phi_{L}\cup\Psi_{L}\cup\dots\cup\Psi_{1}$ basis,
it takes on the form
${\bm{A}}_{L}={\bm{H}}={\bm{U}}_{L}\dots\bm{U}_{2}{\bm{U}}_{1}{\bm{A}}{\bm{U}}_{1}^{T}{\bm{U}}_{2}^{T}\dots\bm{U}_{L}^{T}$
that consists of four distinct blocks (supposingly that we permute the
rows/columns accordingly):

${\bm{H}}=\begin{pmatrix}{\bm{H}}_{\Phi,\Phi}&{\bm{H}}_{\Phi,\Psi}\\\
{\bm{H}}_{\Psi,\Phi}&{\bm{H}}_{\Psi,\Psi}\end{pmatrix},$

where
${\bm{H}}_{\Phi,\Phi}\in\mathbb{R}^{\text{dim}({\mathbb{V}}_{L})\times\text{dim}({\mathbb{V}}_{L})}$
is effectively ${\bm{A}}$ compressed to ${\mathbb{V}}_{L}$,
${\bm{H}}_{\Phi,\Psi}={\bm{H}}_{\Psi,\Phi}^{T}=0$ and ${\bm{H}}_{\Psi,\Psi}$
is diagonal. MMF approximates ${\bm{A}}$ in the form

${\bm{A}}\approx\sum_{i,j=1}^{d_{L}}h_{i,j}\phi_{i}^{L}{\phi_{j}^{L}}^{T}+\sum_{\ell=1}^{L}\sum_{m=1}^{d_{\ell}}c_{m}^{\ell}\psi_{m}^{\ell}{\psi_{m}^{\ell}}^{T},$

where $h_{i,j}$ coefficients are the entries of the ${\bm{H}}_{\Phi,\Phi}$
block, and
$c_{m}^{\ell}=\langle\psi_{m}^{\ell},{\bm{A}}\psi_{m}^{\ell}\rangle$ wavelet
frequencies are the diagonal elements of the ${\bm{H}}_{\Psi,\Psi}$ block.

In particular, the dictionary vectors corresponding to certain rows of
${\bm{U}}_{1}$ are interpreted as level one wavelets, the dictionary vectors
corresponding to certain rows of ${\bm{U}}_{2}{\bm{U}}_{1}$ are interpreted as
level two wavelets, and so on. One thing that is immediately clear is that
whereas Eq. (1) diagonalizes ${\bm{A}}$ in a single step, multiresolution
analysis will involve a sequence of basis transforms
${\bm{U}}_{1},{\bm{U}}_{2},\dots,{\bm{U}}_{L}$, transforming ${\bm{A}}$ step
by step as

${\bm{A}}\rightarrow{\bm{U}}_{1}{\bm{A}}{\bm{U}}_{1}^{T}\rightarrow\dots\rightarrow{\bm{U}}_{L}\dots\bm{U}_{1}{\bm{A}}{\bm{U}}_{1}^{T}\dots\bm{U}_{L}^{T}\triangleq{\bm{H}},$
(4)

so the corresponding matrix factorization must be a multilevel factorization

${\bm{A}}\approx{\bm{U}}_{1}^{T}{\bm{U}}_{2}^{T}\dots\bm{U}_{\ell}^{T}{\bm{H}}{\bm{U}}_{\ell}\dots\bm{U}_{2}{\bm{U}}_{1}.$
(5)

$\textstyle{L_{2}({\mathbb{X}})\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\cdots\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{{\mathbb{V}}_{0}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{{\mathbb{V}}_{1}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{{\mathbb{V}}_{2}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\cdots}$$\textstyle{{\mathbb{W}}_{1}}$$\textstyle{{\mathbb{W}}_{2}}$$\textstyle{{\mathbb{W}}_{3}}$
Figure 1:  Multiresolution analysis splits each function space
${\mathbb{V}}_{0},{\mathbb{V}}_{1},\dots$ into the direct sum of a smoother
part ${\mathbb{V}}_{\ell+1}$ and a rougher part ${\mathbb{W}}_{\ell+1}$.

###  3.3 MMF optimization problem

Finding the best MMF factorization to a symmetric matrix ${\bm{A}}$ involves
solving

$\min_{\begin{subarray}{c}{\mathbb{S}}_{L}\subseteq\cdots\subseteq{\mathbb{S}}_{1}\subseteq{\mathbb{S}}_{0}=[n]\\\
{\bm{H}}\in{\mathbb{H}}^{{\mathbb{S}}_{L}}_{n}\mathchar
24635\relax\;{\bm{U}}_{1},\dots,{\bm{U}}_{L}\in{\mathbb{O}}\end{subarray}}\|{\bm{A}}-{\bm{U}}_{1}^{T}\dots\bm{U}_{L}^{T}{\bm{H}}{\bm{U}}_{L}\dots\bm{U}_{1}\|.$
(6)

Assuming that we measure error in the Frobenius norm, (6) is equivalent to

$\min_{\begin{subarray}{c}{\mathbb{S}}_{L}\subseteq\cdots\subseteq{\mathbb{S}}_{1}\subseteq{\mathbb{S}}_{0}=[n]\\\
{\bm{U}}_{1},\dots,{\bm{U}}_{L}\in{\mathbb{O}}\end{subarray}}\|{\bm{U}}_{L}\dots\bm{U}_{1}{\bm{A}}{\bm{U}}_{1}^{T}\dots\bm{U}_{L}^{T}\|^{2}_{\text{resi}},$
(7)

where $\|\cdot\|_{\text{resi}}^{2}$ is the squared residual norm
$\|{\bm{H}}\|_{\text{resi}}^{2}=\sum_{i\neq j\mathchar
24635\relax\;(i,j)\not\in{\mathbb{S}}_{L}\times{\mathbb{S}}_{L}}\lvert{\bm{H}}_{i,j}\rvert^{2}$.
The optimization problem in (6) and (7) is equivalent to the following 2-level
one:

$\min_{{\mathbb{S}}_{L}\subseteq\cdots\subseteq{\mathbb{S}}_{1}\subseteq{\mathbb{S}}_{0}=[n]}\min_{{\bm{U}}_{1},\dots,{\bm{U}}_{L}\in{\mathbb{O}}}\|{\bm{U}}_{L}\dots\bm{U}_{1}{\bm{A}}{\bm{U}}_{1}^{T}\dots\bm{U}_{L}^{T}\|^{2}_{\text{resi}}.$
(8)

There are two fundamental problems in solving this 2-level optimization:

  * •

For the inner optimization, the variables (i.e. Givens rotations
${\bm{U}}_{1},\dots,{\bm{U}}_{L}$) must satisfy the orthogonality constraints.

  * •

For the outer optimization, finding the optimal nested sequence of indices
${\mathbb{S}}_{L}\subseteq\cdots\subseteq{\mathbb{S}}_{1}\subseteq{\mathbb{S}}_{0}=[n]$
is a combinatorics problem, given an exponential search space.

In order to address these above problems, (Hy & Kondor, 2022) proposes a
learning algorithm combining Stiefel manifold optimization and Reinforcement
Learning (RL) for the inner and outer optimization, respectively. In this
paper, we assume that a nested sequence of indices
${\mathbb{S}}_{L}\subseteq\cdots\subseteq{\mathbb{S}}_{1}\subseteq{\mathbb{S}}_{0}=[n]$
is given by a fast heuristics instead of computationally expensive RL. There
are several heuristics to find the nested sequence, for example: clustering
based on similarity between rows (Kondor et al., 2014) (Kondor et al., 2015).
In the next section, we introduce the solution for the inner problem.

###  3.4 Stiefel manifold optimization

In order to solve the inner optimization problem of (8), we consider the
following generic optimization with orthogonality constraints:

$\min_{{\bm{X}}\in\mathbb{R}^{n\times p}}\mathcal{F}({\bm{X}}),\ \
\text{s.t.}\ \ {\bm{X}}^{T}{\bm{X}}={\bm{I}}_{p},$ (9)

where ${\bm{I}}_{p}$ is the identity matrix and
$\mathcal{F}({\bm{X}}):\mathbb{R}^{n\times p}\rightarrow\mathbb{R}$ is a
differentiable function. The feasible set
$\mathcal{V}_{p}(\mathbb{R}^{n})=\\{{\bm{X}}\in\mathbb{R}^{n\times
p}:{\bm{X}}^{T}{\bm{X}}={\bm{I}}_{p}\\}$ is referred to as the Stiefel
manifold of $p$ orthonormal vectors in $\mathbb{R}^{n}$. We will view
$\mathcal{V}_{p}(\mathbb{R}^{n})$ as an embedded submanifold of
$\mathbb{R}^{n\times p}$. In the case there are more than one orthogonal
constraints, (9) is written as

$\min_{{\bm{X}}_{1}\in\mathcal{V}_{p_{1}}(\mathbb{R}^{n_{1}}),\dots,{\bm{X}}_{q}\in\mathcal{V}_{p_{q}}(\mathbb{R}^{n_{q}})}\mathcal{F}({\bm{X}}_{1},\dots,{\bm{X}}_{q})$
(10)

where there are $q$ variables with corresponding $q$ orthogonal constraints.
In the MMF optimization problem (8), suppose we are already given
${\mathbb{S}}_{L}\subseteq\cdots\subseteq{\mathbb{S}}_{1}\subseteq{\mathbb{S}}_{0}=[n]$
meaning that the indices of active rows/columns at each resolution were
already determined, for simplicity. In this case, we have $q=L$ number of
variables such that each variable
${\bm{X}}_{\ell}={\bm{O}}_{\ell}\in\mathbb{R}^{k\times k}$, where
${\bm{U}}_{\ell}={\bm{I}}_{n-k}\oplus_{{\mathbb{I}}_{\ell}}{\bm{O}}_{\ell}\in\mathbb{R}^{n\times
n}$ in which ${\mathbb{I}}_{\ell}$ is a subset of $k$ indices from
${\mathbb{S}}_{\ell}$, must satisfy the orthogonality constraint. The
corresponding objective function is

$\mathcal{F}({\bm{O}}_{1},\dots,{\bm{O}}_{L})=\|{\bm{U}}_{L}\dots\bm{U}_{1}{\bm{A}}{\bm{U}}_{1}^{T}\dots\bm{U}_{L}^{T}\|^{2}_{\text{resi}}.$
(11)

Therefore, we can cast the inner problem of (8) as an optimization problem on
the Stiefel manifold, and solve it by the specialized steepest gradient
descent (Tagare, 2011).

##  4 Method

Figure 2: Architecture of Fast Temporal Wavelet Neural Network. WC: graph
wavelet convolution given MMF’s wavelet basis. Figure 3: Architecture for the
Wavelet Convolutional Gated Recurrent Unit. WC: graph wavelet convolution
given MMF’s wavelet basis

###  4.1 Wavelet basis and convolution on graph

Section 3.2 introduces the theory of multiresolution analysis behind MMF as
well as the construction of a sparse wavelet basis for a symmetric matrix
${\bm{A}}\in\mathbb{R}^{n\times n}$. Without the loss of generality, we assume
that ${\bm{A}}$ is a weight matrix of a weighted undirected graph
$\mathcal{G}=(V,E)$ in which $V=\\{v_{1},..,v_{n}\\}$ is the set of vertices
and $E=\\{(v_{i},v_{j})\\}$ is the set of edges with the weight of edge
$(v_{i},v_{j})$ is given by ${\bm{A}}_{i,j}$. Given a graph signal
${\bm{f}}\in\mathbb{R}^{n}$ that is understood as a function
$f:V\rightarrow\mathbb{R}$ defined on the vertices of the graph, the wavelet
transform (up to level $L$) expresses this graph signal, without loss of
generality $f\in{\mathbb{V}}_{0}$, as:

$f(v)=\sum_{\ell=1}^{L}\sum_{m}\alpha_{m}^{\ell}\psi_{m}^{\ell}(v)+\sum_{m}\beta_{m}\phi_{m}^{L}(v),\
\ \ \ \text{for each}\ \ v\in V,$

where $\alpha_{m}^{\ell}=\langle f,\psi_{m}^{\ell}\rangle$ and
$\beta_{m}=\langle f,\phi_{m}^{L}\rangle$ are the wavelet coefficients. Based
on the wavelet basis construction via MMF detailed in (Hy & Kondor, 2022):

  * •

For $L$ levels of resolution, we get exactly $L$ mother wavelets
$\overline{\psi}=\\{\psi^{1},\psi^{2},\dots,\psi^{L}\\}$, each corresponds to
a resolution (see Figure 6 for visualizations).

  * •

The rows of ${\bm{H}}={\bm{A}}_{L}$ make exactly $n-L$ father wavelets
$\overline{\phi}=\\{\phi^{L}_{m}={\bm{H}}_{m,:}\\}_{m\in{\mathbb{S}}_{L}}$. In
total, a graph of $n$ vertices has exactly $n$ wavelets, both mothers and
fathers.

Analogous to the convolution based on Graph Fourier Transform (Bruna et al.,
2014), each convolution layer $k\in\\{1,..,K\\}$ of wavelet neural network
transforms an input vector ${\bm{f}}^{(k-1)}$ of size $\lvert V\rvert\times
F_{k-1}$ into an output ${\bm{f}}^{(k)}$ of size $\lvert V\rvert\times F_{k}$
as

${\bm{f}}^{(k)}_{:,j}=\sigma\bigg{(}{\bm{W}}\sum_{i=1}^{F_{k-1}}{\bm{g}}^{(k)}_{i,j}{\bm{W}}^{T}{\bm{f}}^{(k-1)}_{:,i}\bigg{)}\
\ \ \ \text{for}\ \ j=1,\dots,F_{k},$ (12)

where ${\bm{W}}$ is our wavelet basis matrix as we concatenate
$\overline{\phi}$ and $\overline{\psi}$ column-by-column,
${\bm{g}}^{(k)}_{i,j}$ is a parameter/filter in the form of a diagonal matrix
learned in spectral domain, and $\sigma$ is an element-wise non-linearity
(e.g., ReLU, sigmoid, etc.). In Eq.(12), first we employ the wavelet transform
of a graph signal ${\bm{f}}$ into the spectral domain (i.e.
$\hat{{\bm{f}}}={\bm{W}}^{T}{\bm{f}}$ is the forward transform and
${\bm{f}}={\bm{W}}\hat{{\bm{f}}}$ is the inverse transform), then a learnable
filter ${\bm{g}}$ to the wavelet coefficients, the inverse transform back to
the spatial domain and let everything through a non-linearity $\sigma$. Since
the wavelet basis matrix ${\bm{W}}$ is sparse, both the wavelet transform and
its inverse transform can be implemented efficiently via sparse matrix
multiplication.

###  4.2 Temporal Wavelet Neural Networks

Capturing spatiotemporal dependencies among time series in various
spatiotemporal forecasting problems demands both spatial and temporal models.
We build our novel Fast Temporal Wavelet Graph Neural Network with the
architectural backbone from Diffusion Convolutional Recurrent Neural Network
(DCRNN) (Li et al., 2018), that combines both spatial and temporal models to
solve these tasks.

Spatial Dependency Model The spatial dynamic in the network is captured by
diffusion process. Let $G=({\bm{X}},{\bm{A}})$ represent an undirected graph,
where
${\bm{X}}=[{\bm{x}}_{1}^{T},\dots,{\bm{x}}_{N}^{T}]^{T}\in\mathbb{R}^{N\times
D}$ denotes signals of $N$ nodes, each has $D$ features. Define further the
right-stochastic edge weights matrix $\tilde{{\bm{A}}}\in\mathbb{R}^{N\times
N}$ in which $\sum_{j}\tilde{{\bm{A}}}_{ij}=1\forall i$. In the simplest case,
when $\tilde{{\bm{L}}}={\bm{I}}-\tilde{{\bm{A}}}$ is the nomalized random walk
matrix, the diffusion process on graph is governed by the following equation
(Chamberlain et al., 2021):

$\frac{{\rm d}{\bm{X}}(t)}{{\rm d}t}=(\tilde{{\bm{A}}}-{\bm{I}}){\bm{X}}(t)$
(13)

where
${\bm{X}}(t)=[{\bm{x}}_{1}^{T}(t),\dots,[{\bm{x}}_{N}^{T}(t)]^{T}\in\mathbb{R}^{N\times
D}$ and ${\bm{X}}(0)={\bm{X}}$. Applying forward Euler discretization with
step size 1, gives:

$\displaystyle{\bm{X}}(k)$
$\displaystyle={\bm{X}}(k-1)+(\tilde{{\bm{A}}}-{\bm{I}}){\bm{X}}(k-1)$
$\displaystyle={\bm{X}}(k-1)-\tilde{{\bm{L}}}{\bm{X}}(k-1)$
$\displaystyle=\tilde{{\bm{A}}}{\bm{X}}(k-1)$
$\displaystyle=\tilde{{\bm{A}}}^{k}{\bm{X}}(0)$ (14)

Eq.14 is similar to the well-established GCN architecture propose in (Kipf &
Welling, 2016). Then, the diffusion convolution operation over a graph signal
$\mathbb{R}^{N\times D}$ and filter $f_{\bm{\theta}}$ is defined as:

${\bm{X}}_{:,d}\star_{\mathcal{G}}f_{\bm{\theta}}=\sum_{k=0}^{K-1}\theta_{k}\tilde{{\bm{A}}}^{k}{\bm{X}}_{:,d}\quad\forall
d\in\\{1,\dots,D\\}$ (15)

where ${\bm{\Theta}}\in\mathbb{R}^{K\times 2}$ are the parameters for the
filter.

Temporal Dependency Model The DCRNN is leveraged from the recurrent neural
networks (RNNs) to model the temporal dependency. In particular, the matrix
multiplications in GRU is replaced with the diffusion convolution, which is
called Diffusion Convolutional Gated Recurrent Unit (DCGRU).

$\displaystyle{\bm{r}}^{(t)}$
$\displaystyle=\sigma({\bm{\Theta}}_{r}\star_{\mathcal{G}}[{\bm{X}}^{(t)},{\bm{H}}^{(t-1)}]+\bm{b})$
$\displaystyle\bm{u}^{(t)}$
$\displaystyle=\sigma({\bm{\Theta}}_{r}\star_{\mathcal{G}}[{\bm{X}}^{(t)},{\bm{H}}^{(t-1)}]+\bm{b_{u}})$
$\displaystyle{\bm{C}}^{(t)}$
$\displaystyle=\tanh({\bm{\Theta}}_{r}\star_{\mathcal{G}}[{\bm{X}}^{(t)},(\bm{r}\odot{\bm{H}}^{(t-1)})]+\bm{b_{c}})$
$\displaystyle{\bm{H}}^{(t)}$
$\displaystyle=\bm{u}^{(t)}\odot{\bm{H}}^{(t-1)}+(1-\bm{u}^{(t)})\odot{\bm{C}}^{(t)}$

where ${\bm{X}}(t),{\bm{H}}(t)$ denote the input and output of at time $t$,
while ${\bm{r}}^{(t)},\bm{u}^{(t)}$ are reset gate and update gate at time
$t$, respectively.

Both the encoder and the decoder are recurrent neural networks with DCGRU
following Sequence-to-Sequence style. To mitigate the distribution differences
between training and testing data, scheduled sampling technique (Bengio et
al., 2015) is used, where the model is fed with either the ground truth with
probability $\epsilon_{i}$ or the prediction by the model with probability
$1-\epsilon_{i}$.

For our novel Fast Temporal Wavelet Graph Neural Network (FTWGNN), the
fundamental difference is that instead of using temporal traffic graph as the
input of DCRNN, we use the sparse wavelet basis matrix ${\bm{W}}$ which is
extracted via MMF (see Section 3.2) and replace the diffusion convolution by
our fast wavelet convolution. Given the sparsity of our wavelet basis, we
significantly reduce the overall computational time and memory usage. Each
Givens rotation matrix ${\bm{U}}_{\ell}$ (see Def. 3.1) is a highly-sparse
orthogonal matrix with a non-zero core of size $K\times K$. The number of non-
zeros in MMF’s wavelet basis ${\bm{W}}$, that can be computed as product
${\bm{U}}_{1}{\bm{U}}_{2}\cdots{\bm{U}}_{L}$, is $O(LK^{2})$ where $L$ is the
number of resolutions (i.e. number of Givens rotation matrices) and $K$ is the
number of columns in a Givens rotation matrix. (Kondor et al., 2014) and (Hy &
Kondor, 2022) have shown in both theory and practice that $L$ only needs to be
in $O(n)$ where $n$ is the number of columns and $K$ small (e.g., 2, 4, 8) to
get a decent approximation/compression for a symmetric hierarchical matrix.
Technically, MMF is able to compress a symmetric hierararchical matrix from
the original quadratic size $n\times n$ to a linear number of non-zero
elements $O(n)$. Practically, all the Givens rotation matrices
$\\{{\bm{U}}_{\ell}\\}_{\ell=1}^{L}$ and the wavelet basis ${\bm{W}}$ can be
stored in Coordinate Format (COO), and the wavelet transform and its inverse
in wavelet convolution (see Eq. 12) can be implemented efficiently by sparse
matrix multiplication in PyTorch’s sparse library (Paszke et al., 2019). The
architecture of our model is shown in Figures 2 and 3.

##  5 Experiments

Our PyTorch implementation is publicly available at
https://github.com/HySonLab/TWGNN. Implementation of Multiresolution Matrix
Factorization and graph wavelet computation (Hy & Kondor, 2022) is publicly
available at https://github.com/risilab/Learnable_MMF.

To showcase the competitive performance and remarkable acceleration of FTWGNN,
we conducted experiments on two well-known traffic forecasting benchmarks
METR-LA and PEMS-BAY, and one challenging ECoG dataset AJILE12. Following (Li
et al., 2018), we compare our model with widely used time series models,
including:  1. HA: Historical Average, which models the traffic flow as a
seasonal process, and uses weighted average of previous seasons as the
prediction;  2. $\text{ARIMA}_{kal}$: Auto-Regressive Integrated Moving
Average model with Kalman filter;  3. VAR: Vector Auto-Regression.  4. SVR:
Support Vector Regression.  5. FNN: Feed forward neural network with two
hidden layers and L2 regularization.  6. FC-LSTM: Recurrent Neural Network
with fully connected LSTM hidden units(Sutskever et al., 2014) Methods are
evaluated on three metrics:  (i) Mean Absolute Error (MAE);  (ii) Mean
Absolute Percentage Error (MAPE); and  (iii) Root Mean Squared Error (RMSE).
FTWGNN and DCRNN are implemented using PyTorch (Paszke et al., 2019) on an
NVIDIA A100-SXM4-80GB GPU. Details about parameters and model structure can be
found in the Appendix A and B.

Dataset | $T$ | Metric | HA | $\text{ARIMA}_{kal}$ | VAR | SVR | FNN | FC-LSTM | DCRNN | FTWGNN
---|---|---|---|---|---|---|---|---|---|---
METR-LA |  | MAE | 4.16 | 3.99 | 4.42 | 3.99 | 3.99 | 3.44 | 2.77 | 2.70
15 min | RMSE | 7.80 | 8.21 | 7.89 | 8.45 | 7.94 | 6.30 | 5.38 | 5.15
| MAPE | 13.0% | 9.6% | 10.2% | 9.3% | 9.9% | 9.6% | 7.3% | 6.81%
| MAE | 4.16 | 5.15 | 5.41 | 5.05 | 4.23 | 3.77 | 3.15 | 3.0
30 min | RMSE | 7.80 | 10.45 | 9.13 | 10.87 | 8.17 | 7.23 | 6.45 | 5.95
| MAPE | 13.0% | 12.7% | 12.7% | 12.1% | 12.9% | 10.9% | 8.80% | 8.01%
| MAE | 4.16 | 6.90 | 6.52 | 6.72 | 4.49 | 4.37 | 3.60 | 3.42
60 min | RMSE | 7.80 | 13.23 | 10.11 | 13.76 | 8.69 | 8.69 | 7.59 | 6.92
| MAPE | 13.0% | 17.4% | 15.8% | 16.7% | 14.0% | 13.2% | 10.5% | 9.83%
PEMS-BAY |  | MAE | 2.88 | 1.62 | 1.74 | 1.85 | 2.20 | 2.05 | 1.38 | 1.14
15 min | RMSE | 5.59 | 3.30 | 3.16 | 3.59 | 4.42 | 4.19 | 2.95 | 2.40
| MAPE | 6.8% | 3.5% | 3.6% | 3.8% | 5.2% | 4.8% | 2.9% | 2.28%
| MAE | 2.88 | 2.33 | 2.32 | 2.48 | 2.30 | 2.20 | 1.74 | 1.5
30 min | RMSE | 5.59 | 4.76 | 4.25 | 5.18 | 4.63 | 4.55 | 3.97 | 3.27
| MAPE | 6.8% | 5.4% | 5.0% | 5.5% | 5.43% | 5.2% | 3.9% | 3.15%
| MAE | 2.88 | 3.38 | 2.93 | 3.28 | 2.46 | 2.37 | 2.07 | 1.79
60 min | RMSE | 5.59 | 6.5 | 5.44 | 7.08 | 4.98 | 4.96 | 4.74 | 3.99
| MAPE | 6.8% | 8.3% | 6.5% | 8.0% | 5.89% | 5.7% | 4.9% | 4.14%

Table 1: Performance comparison of different models for traffic speed forecasting. Dataset | $T$ | DCRNN | FTWGNN | Speedup
---|---|---|---|---
METR-LA | 15 min | 350s | 217s | 1.61x
30 min | 620s | 163s | 3.80x
60 min | 1800s | 136s | 13.23x
PEMS-BAY | 15 min | 427s | 150s | 2.84x
30 min | 900s | 173s | 5.20x
60 min | 1800s | 304s | 5.92x
AJILE12 | 1 sec | 80s | 35s | 2.28x
5 sec | 180s | 80s | 2.25x
15 sec | 350s | 160s | 2.18x
Table 2: Training time/epoch between DCRNN and FTWGNN.

Adjacency matrix According to DCRNN (Li et al., 2018), the traffic sensor
network is expressed by an adjacency matrix which is constructed using
thresholded Gaussian kernel (Shuman et al., 2013). Specifically, for each pair
of sensors $v_{i}$ and $v_{j}$, the edge weight from $v_{i}$ to $v_{j}$,
denoted by $A_{ij}$, is defined as

$A_{ij}\coloneqq\left\\{\begin{aligned}
&\exp(-\frac{\text{dist}(v_{i},v_{j})}{\sigma^{2}}),\quad&\text{dist}(v_{i},v_{j})\leq
k\\\ &0,\quad&\text{otherwise}\end{aligned}\right.\,,$ (16)

where $\text{dist}(v_{i},v_{j})$ denotes the spatial distance from $v_{i}$ to
$v_{j}$, $\sigma$ is the standard deviation of the distances, and $k$ is the
distance threshold.

Nevertheless, such an user-defined adjacency matrix requires expert knowledge,
thus may not work on other domains, e.g., brain networks. In the ECoG time
series forecasting case, the adjacency matrix is computed based on the popular
Local Linear Embedding (LLE) (Saul & Roweis, 2003). In particular, for the
matrix data
${\bm{X}}=[{\bm{x}}_{1},\ldots,{\bm{x}}_{N}]\in{\mathbb{R}}^{T\times N}$ where
${\bm{x}}_{i}$ denotes the time series data of node $i$ for
$i\in\\{1,\ldots,N\\}$, an adjacency matrix ${\bm{A}}$ is identified to gather
all the coefficients of the affine dependencies among
$\\{{\bm{x}}_{i}\\}_{i=1}^{N}$ by solving the following optimization problem

$\displaystyle{\bm{A}}\coloneqq\arg\min_{\hat{{\bm{A}}}\in{\mathbb{R}}^{N\times
N}}{}$
$\displaystyle{}\norm{{\bm{X}}-{\bm{X}}\hat{{\bm{A}}}^{T}}_{\text{F}}^{2}+\lambda_{A}\norm{\hat{{\bm{A}}}}_{1}$
s.to
$\displaystyle{}{\bm{1}}^{T}_{N}\hat{{\bm{A}}}={\bm{1}}^{T}_{N}\,,\quad\text{diag}(\hat{{\bm{A}}})={\bm{0}}\,,$
(17)

where the constraint ${\bm{1}}^{T}_{N}\hat{{\bm{A}}}={\bm{1}}^{T}_{N}$
realizes the affine combitions, while $\text{diag}(\hat{{\bm{A}}})={\bm{0}}$
excludes the trivial solution $\hat{{\bm{A}}}={\bm{I}}_{N}$. Furthermore, to
promote the local smoothness of the graph, each data point $x_{i}$ is assumed
to be approximated by a few neighbors
$\\{x_{j_{1}},x_{j_{2}},\ldots,x_{j_{k}}\\}$, thus $\hat{A}$ is regularized by
the $l_{1}$-norm loss $\norm{\hat{{\bm{A}}}}_{1}$ to be sparse. Task (17) is a
composite convex minimization problem with affine constraints, hence can be
solved by (Slavakis & Yamada, 2018).

###  5.1 Traffic prediction

Two real-world large-scale traffic datasets are considered:

  * •

METR-LA Data of 207 sensors in the highway of Los Angeles County (Jagadish et
al., 2014) over the 4-month period from Mar 1st 2012 to Jun 30th 2012.

  * •

PEMS-BAY Data of 325 sensors in the Bay Area over the 6-month period from Jan
1st 2017 to May 31th 2017 from the California Transportation Agencies
(CalTrans) Performance Measurement System (PeMS).

Figure 4: (a) Augmented ECoG signals by spline interpolation envelopes, (b)
1-second prediction of ECoG signals

Dataset | $T$ | Metric | HA | VAR | LR | SVR | LSTM | DCRNN | FTWGNN
---|---|---|---|---|---|---|---|---|---
AJILE12 |  | MAE | 0.88 | 0.16 | 0.27 | 0.27 | 0.07 | 0.05 | 0.03
1 sec | RMSE | 1.23 | 0.25 | 0.37 | 0.41 | 0.09 | 0.45 | 0.35
| MAPE | 320% | 58% | 136% | 140% | 38% | 7.84% | 5.27%
| MAE | 0.88 | 0.66 | 0.69 | 0.69 | 0.39 | 0.16 | 0.11
5 sec | RMSE | 1.23 | 0.96 | 0.92 | 0.93 | 0.52 | 0.24 | 0.15
| MAPE | 320% | 221% | 376% | 339% | 147% | 64% | 57%
| MAE | 0.88 | 0.82 | 0.86 | 0.86 | 0.87 | 0.78 | 0.70
15 sec | RMSE | 1.23 | 1.15 | 1.13 | 1.13 | 1.14 | 1.01 | 0.93
| MAPE | 320% | 320% | 448% | 479% | 330% | 294% | 254%

Table 3:  Performance comparison on ECoG signals forecast.

The distance function $\text{dist}(v_{i},v_{j})$ in (16) represents the road
network distance from sensor $v_{i}$ to sensor $v_{j}$, producing an
asymmetric adjacency matrix for a directed graph. Therefore, the symmetrized
matrix $\hat{{\bm{A}}}\coloneqq\frac{1}{2}({\bm{A}}+{\bm{A}}^{T})$ is taken to
compute the wavelet basis matrix ${\bm{W}}$ following Sec. 3.2. Parameters can
be found in the Appendix B.

Table 1 shows the evaluation of different approaches on the two traffic
datasets, while Table 2 reports the training time per epoch of FTWGNN and
DCRNN. Overall, although FTWGNN performs better than DCRNN by about only 10%,
it is significantly faster by about 5 times on average.

###  5.2 Brain networks

Annotated Joints in Long-term Electrocorticography for 12 human participants
(AJILE12), publicly available at (Peterson et al., 2022b), records
intracranial neural activity via the invasive ECoG, which involves implanting
electrodes directly under the skull (Peterson et al., 2022a). For each
participant, ECoG recordings are sporadically sampled at 500Hz in $7.4\pm 2.2$
days (mean$\pm$std) from at least 64 electrodes, each of which is encoded with
an unique set of Montreal Neurological Institute (MNI) x, y, z coordinates.

The proposed model is tested on the first one hour of recordings of subject
number 5 with 116 good-quality electrodes. Signals are downsampled to 1Hz,
thus producing a network of 116 nodes, each with $3{,}600$ data points.
Furthermore, the signals are augmented by applying the spline interpolation to
get the upper and lower envelopes along with an average curve (Melia et al.,
2014) (see Figure 4). The adjacency matrix ${\bm{A}}$ is obtained by solving
task (17), then the wavelet basis matrix ${\bm{W}}$ is constructed based on
Sec. 3.2. Parameters can be found in the Appendix B.

Table 3 reports the performance of different methods on the AJILE12 dataset
for 1-, 5-, and 15-second prediction. Generally, errors are much higher than
those in the traffic forecasting problem, since the connections within the
brain network are much more complicated and ambiguous (Breakspear, 2017). High
errors using HA and VAR methods show that the AJILE12 data follows no
particular pattern or periodicity, making long-step prediction extremely
challenging. Despite having a decent performance quantitatively, Figure 4
demonstrates the superior performance of FTWGNN, in which DCRNN fails to
approximate the trend and the magnitude of the signals. Even though FTWGNN
performs well at 1-second prediction, it produces unstable and erroneous
forecast at longer steps of 5 or 15 seconds (see Figure 5). Meanwhile, similar
to traffic prediction case, FTWGNN also sees a remarkable improvement in
computation time by around 2 times on average (see Table 2).

Table 4 shows the sparsity density of our MMF wavelet basis ${\bm{W}}$ for
each dataset. It is important to note that ${\bm{W}}$ is extremely sparse with
no more than $2\%$ of non-zeros, that makes our models run much faster with
significantly low memory usage while achieving competitive results.

##  6 Conclusion

We propose a new class of spatial-temporal graph neural networks based on the
theories of multiresolution analysis and wavelet theory on discrete spaces
with RNN backbone, coined Fast Temporal Wavelet Graph Neural Network (FTWGNN).
Fundamentally, we employ Multiresolution Matrix Factorization to factorize the
underlying graph structure and extract its corresponding sparse wavelet basis
that consequentially allows us to construct efficient wavelet transform and
convolution on graph. Experiments on real-world large-scale datasets show
promising results and computational efficiency of FTGWNN in network time
series modeling including traffic prediction and brain networks. Several
future directions are:  (i) investigating synchronization phenomena in brain
networks (Honda, 2018);  (ii) developing a robust model against
outliers/missing data that appear frequently in practice.

## References

  * Bengio et al. (2015) Bengio, S., Vinyals, O., Jaitly, N., and Shazeer, N.  Scheduled sampling for sequence prediction with recurrent neural networks.  _Advances in neural information processing systems_ , 28, 2015.
  * Bessadok et al. (2022) Bessadok, A., Mahjoub, M. A., and Rekik, I.  Graph neural networks in network neuroscience.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2022.
  * Bickel & Ritov (2008) Bickel, P. J. and Ritov, Y.  Discussion of: Treelets—an adaptive multi-scale basis for sparse unordered data.  _The Annals of Applied Statistics_ , 2(2):474–477, 2008.  ISSN 19326157.  URL http://www.jstor.org/stable/30244209.
  * Breakspear (2017) Breakspear, M.  Dynamic models of large-scale brain activity.  _Nature neuroscience_ , 20(3):340–352, 2017.
  * Bruna et al. (2014) Bruna, J., Zaremba, W., Szlam, A., and Lecun, Y.  Spectral networks and locally connected networks on graphs.  In _International Conference on Learning Representations (ICLR2014), CBLS, April 2014_ , 2014.
  * Chamberlain et al. (2021) Chamberlain, B., Rowbottom, J., Gorinova, M. I., Bronstein, M., Webb, S., and Rossi, E.  Grand: Graph neural diffusion.  In _International Conference on Machine Learning_ , pp. 1407–1418. PMLR, 2021.
  * Coifman & Maggioni (2006) Coifman, R. R. and Maggioni, M.  Diffusion wavelets.  _Applied and Computational Harmonic Analysis_ , 21(1):53–94, 2006.  ISSN 1063-5203.  doi: https://doi.org/10.1016/j.acha.2006.04.004.  URL https://www.sciencedirect.com/science/article/pii/S106352030600056X.  Special Issue: Diffusion Maps and Wavelets.
  * Daubechies (1988) Daubechies, I.  Orthonormal bases of compactly supported wavelets.  _Communications on Pure and Applied Mathematics_ , 41:909–996, 1988.
  * Ding et al. (2017) Ding, Y., Kondor, R., and Eskreis-Winkler, J.  Multiresolution kernel approximation for gaussian process regression.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper/2017/file/850af92f8d9903e7a4e0559a98ecc857-Paper.pdf.
  * Gavish et al. (2010) Gavish, M., Nadler, B., and Coifman, R. R.  Multiscale wavelets on trees, graphs and high dimensional data: Theory and applications to semi supervised learning.  In _Proceedings of the 27th International Conference on International Conference on Machine Learning_ , ICML’10, pp. 367–374, Madison, WI, USA, 2010. Omnipress.  ISBN 9781605589077.
  * Göktaş et al. (2020) Göktaş, A. S., Bessadok, A., and Rekik, I.  Residual embedding similarity-based network selection for predicting brain network evolution trajectory from a single observation.  In _International Workshop on PRedictive Intelligence In MEdicine_ , pp. 12–23. Springer, 2020.
  * Gürler et al. (2020) Gürler, Z., Nebli, A., and Rekik, I.  Foreseeing brain graph evolution over time using deep adversarial network normalizer.  In _International Workshop on Predictive Intelligence In Medicine_ , pp. 111–122. Springer, 2020.
  * Haar (1910) Haar, A.  Zur theorie der orthogonalen funktionensysteme.  _Mathematische Annalen_ , 69:331–371, 1910.
  * Hamilton (2020) Hamilton, J. D.  _Time series analysis_.  Princeton university press, 2020.
  * Hammond et al. (2011) Hammond, D. K., Vandergheynst, P., and Gribonval, R.  Wavelets on graphs via spectral graph theory.  _Applied and Computational Harmonic Analysis_ , 30(2):129–150, 2011.  ISSN 1063-5203.  doi: https://doi.org/10.1016/j.acha.2010.04.005.  URL https://www.sciencedirect.com/science/article/pii/S1063520310000552.
  * Honda (2018) Honda, H.  On mathematical modeling and analysis of brain network.  In _Mathematical Analysis of Continuum Mechanics and Industrial Applications II: Proceedings of the International Conference CoMFoS16 16_ , pp. 169–180. Springer, 2018.
  * Hy & Kondor (2022) Hy, T. S. and Kondor, R.  Multiresolution matrix factorization and wavelet networks on graphs.  In Cloninger, A., Doster, T., Emerson, T., Kaul, M., Ktena, I., Kvinge, H., Miolane, N., Rice, B., Tymochko, S., and Wolf, G. (eds.), _Proceedings of Topological, Algebraic, and Geometric Learning Workshops 2022_ , volume 196 of _Proceedings of Machine Learning Research_ , pp. 172–182. PMLR, 25 Feb–22 Jul 2022.  URL https://proceedings.mlr.press/v196/hy22a.html.
  * Ithapu et al. (2017) Ithapu, V. K., Kondor, R., Johnson, S. C., and Singh, V.  The incremental multiresolution matrix factorization algorithm.  In _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 692–701, 2017.  doi: 10.1109/CVPR.2017.81.
  * Jacobi (1846) Jacobi, C.  Über ein leichtes verfahren die in der theorie der säcularstörungen vorkommenden gleichungen numerisch aufzulösen*).  1846(30):51–94, 1846.  doi: doi:10.1515/crll.1846.30.51.  URL https://doi.org/10.1515/crll.1846.30.51.
  * Jagadish et al. (2014) Jagadish, H. V., Gehrke, J., Labrinidis, A., Papakonstantinou, Y., Patel, J. M., Ramakrishnan, R., and Shahabi, C.  Big data and its technical challenges.  _Commun. ACM_ , 57(7):86–94, 2014.  URL http://dblp.uni-trier.de/db/journals/cacm/cacm57.html#JagadishGLPPRS14.
  * Jeong et al. (2013) Jeong, Y.-S., Byon, Y.-J., Castro-Neto, M. M., and Easa, S. M.  Supervised weighting-online learning algorithm for short-term traffic flow prediction.  _IEEE Transactions on Intelligent Transportation Systems_ , 14(4):1700–1707, 2013.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kipf & Welling (2016) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ , 2016.
  * Kondor et al. (2014) Kondor, R., Teneva, N., and Garg, V.  Multiresolution matrix factorization.  In Xing, E. P. and Jebara, T. (eds.), _Proceedings of the 31st International Conference on Machine Learning_ , volume 32 of _Proceedings of Machine Learning Research_ , pp. 1620–1628, Bejing, China, 22–24 Jun 2014\. PMLR.  URL https://proceedings.mlr.press/v32/kondor14.html.
  * Kondor et al. (2015) Kondor, R., Teneva, N., and Mudrakarta, P. K.  Parallel MMF: a multiresolution approach to matrix computation.  _CoRR_ , abs/1507.04396, 2015.  URL http://arxiv.org/abs/1507.04396.
  * Koprinska et al. (2018) Koprinska, I., Wu, D., and Wang, Z.  Convolutional neural networks for energy time series forecasting.  In _2018 International Joint Conference on Neural Networks (IJCNN)_ , pp. 1–8, 2018.  doi: 10.1109/IJCNN.2018.8489399.
  * Li et al. (2018) Li, Y., Yu, R., Shahabi, C., and Liu, Y.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=SJiHXGWAZ.
  * Mallat (1989) Mallat, S.  A theory for multiresolution signal decomposition: the wavelet representation.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 11(7):674–693, 1989.  doi: 10.1109/34.192463.
  * Mallat (2008) Mallat, S.  _A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way_.  Academic Press, Inc., USA, 3rd edition, 2008.  ISBN 0123743702.
  * Melia et al. (2014) Melia, U., Clariá, F., Vallverdú, M., and Caminal, P.  Filtering and thresholding the analytic signal envelope in order to improve peak and spike noise reduction in eeg signals.  _Medical engineering & physics_, 36(4):547–553, 2014.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Peterson et al. (2022a) Peterson, S. M., Singh, S. H., Dichter, B., Scheid, M., Rao, R. P. N., and Brunton, B. W.  Ajile12: Long-term naturalistic human intracranial neural recordings and pose.  _Scientific Data_ , 9(1):184, Apr 2022a.  ISSN 2052-4463.  doi: 10.1038/s41597-022-01280-y.  URL https://doi.org/10.1038/s41597-022-01280-y.
  * Peterson et al. (2022b) Peterson, S. M., Singh, S. H., Dichter, B., Scheid, M., Rao, R. P. N., and Brunton, B. W.  Ajile12: Long-term naturalistic human intracranial neural recordings and pose (Version 0.220127.0436) [Data set].  _DANDI archive_ , 2022b.  URL https://doi.org/10.48324/dandi.000055/0.220127.0436.
  * Pourahmadi & Noorbaloochi (2016) Pourahmadi, M. and Noorbaloochi, S.  Multivariate time series analysis of neuroscience data: some challenges and opportunities.  _Current Opinion in Neurobiology_ , 37:12–15, 2016.  ISSN 0959-4388.  doi: https://doi.org/10.1016/j.conb.2015.12.006.  URL https://www.sciencedirect.com/science/article/pii/S0959438815001841.  Neurobiology of cognitive behavior.
  * Qin et al. (2017) Qin, Y., Song, D., Cheng, H., Cheng, W., Jiang, G., and Cottrell, G. W.  A dual-stage attention-based recurrent neural network for time series prediction.  In _Proceedings of the 26th International Joint Conference on Artificial Intelligence_ , IJCAI’17, pp. 2627–2633. AAAI Press, 2017.  ISBN 9780999241103.
  * Saul & Roweis (2003) Saul, L. K. and Roweis, S. T.  Think globally, fit locally: unsupervised learning of low dimensional manifolds.  _Journal of machine learning research_ , 4(Jun):119–155, 2003.
  * Shoeibi et al. (2022) Shoeibi, A., Moridian, P., Khodatars, M., Ghassemi, N., Jafari, M., Alizadehsani, R., Kong, Y., Gorriz, J. M., Ramírez, J., Khosravi, A., et al.  An overview of deep learning techniques for epileptic seizures detection and prediction based on neuroimaging modalities: Methods, challenges, and future works.  _Computers in Biology and Medicine_ , pp. 106053, 2022.
  * Shuman et al. (2013) Shuman, D. I., Narang, S. K., Frossard, P., Ortega, A., and Vandergheynst, P.  The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains.  _IEEE Signal Processing Magazine_ , 30(3):83–98, 2013.  doi: 10.1109/MSP.2012.2235192.
  * Slavakis & Yamada (2018) Slavakis, K. and Yamada, I.  Fejér-monotone hybrid steepest descent method for affinely constrained and composite convex minimization tasks.  _Optimization_ , 67(11):1963–2001, 2018.
  * Smola & Schölkopf (2004) Smola, A. J. and Schölkopf, B.  A tutorial on support vector regression.  _Statistics and Computing_ , 14(3):199–222, aug 2004.  ISSN 0960-3174.  doi: 10.1023/B:STCO.0000035301.49549.88.  URL https://doi.org/10.1023/B:STCO.0000035301.49549.88.
  * Sutskever et al. (2014) Sutskever, I., Vinyals, O., and Le, Q. V.  Sequence to sequence learning with neural networks.  _Advances in neural information processing systems_ , 27, 2014.
  * Tagare (2011) Tagare, H. D.  Notes on optimization on stiefel manifolds.  2011\.
  * Teneva et al. (2016) Teneva, N., Mudrakarta, P. K., and Kondor, R.  Multiresolution matrix compression.  In Gretton, A. and Robert, C. C. (eds.), _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics_ , volume 51 of _Proceedings of Machine Learning Research_ , pp. 1441–1449, Cadiz, Spain, 09–11 May 2016. PMLR.  URL https://proceedings.mlr.press/v51/teneva16.html.
  * Van Lint & Van Hinsbergen (2012) Van Lint, J. and Van Hinsbergen, C.  Short-term traffic and travel time prediction models.  _Artificial Intelligence Applications to Critical Transportation Issues_ , 22(1):22–41, 2012.
  * Xu et al. (2020) Xu, M., Dai, W., Liu, C., Gao, X., Lin, W., Qi, G.-J., and Xiong, H.  Spatial-temporal transformer networks for traffic flow forecasting.  _arXiv preprint arXiv:2001.02908_ , 2020.
  * Zhang et al. (2017) Zhang, J., Zheng, Y., and Qi, D.  Deep spatio-temporal residual networks for citywide crowd flows prediction.  In _Thirty-first AAAI conference on artificial intelligence_ , 2017\.

##  Appendix A Baseline settings

HA Historical Average, which models the traffic flow as a seasonal/periodical
process, and uses weighted average of previous seasons/periods as the
prediction.

$\textbf{ARIMA}_{kal}$ Auto-Regressive Integrated Moving Average model with
Kalman filter, implemented by the statsmodel package in Python.

VAR Vector Auto-regressive model (Hamilton, 2020) with orders $(3,0,1)$,
implemented by the Python statsmodel package.

LR Linear Regression with 5 historical observations.

SVR Linear Support Vector Regression (Smola & Schölkopf, 2004) with 5
historical observations.

FNN Feed forward neural network with two hidden layers, each with 256 units.
The initial learning rate is $1e^{-3}$, and the decay rate is $1e^{-1}$ per 20
epochs. In addition, for all hidden layers, dropout with ratio 0.5 and L2
weight decay $1e^{-2}$ is used. The model is trained to minimize the MAE with
batch size 64.

FC-LSTM The encoder-decoder framework using LSTM with peephole (Sutskever et
al., 2014). the encoder and the decoder contain two recurrent layers, each of
which consists of 256 LSTM units, with L1 weight decay rate $2e^{-5}$ and L2
weight decay rate $5e^{-4}$. The initial learning rate is $1e^{-4}$, and the
decay rate is $1e^{-1}$ per 20 epochs.

DCRNN Settings of the Diffusion Convolutional Recurrent Neural Network follow
its original work (Li et al., 2018).

##  Appendix B FTWGNN settings

Data preparation For all datasets, the train/validation/test ratio is
$0.7/0.2/0.1$, divided into batch size of 64.

Adjacency matrix The $k$-neighborhood of the traffic network in Eq. 16 is
thresholded by $k=0.01$, while for the brain network, the parameter
$\lambda_{A}$ is set to $1e^{-5}$ in Task (17).

Wavelet basis For the traffic datasets, 100 mother wavelets are extracted,
i.e., $L=100$, while for the AJILE12 dataset, $L=10$ was used.

Model architecture For the wavelet convolution RNN, both encoder and decoder
contains two recurrent layers, each with 64 units. The initial learning rate
is $1e^{-2}$, decaying by $\frac{1}{10}$ per 20 epochs; the dropout ratio is
$0.1$; and the maximum diffusion step, i.e., $K$, is set to 2. Besides, the
optimizer is the Adam optimizer (Kingma & Ba, 2014).

Dataset | Fourier basis | Wavelet basis
---|---|---
METR-LA | 99.04% | 1.11%
PEMS-BAY | 96.35% | 0.63%
AJILE12 | 100% | 1.81%
Table 4: Sparsity of wavelet basis (i.e. percentage of non-zeros).

##  Appendix C Vizualization of the AJILE12 data

Figure 5: 5-second prediction of ECoG signals.

(a)

(b)

(c)

Figure 6: Visualization of some of the wavelets on the brain network of 116
nodes. The low index wavelets (low $\ell$) are highly localized, whereas the
high index ones are smoother and more dispersed.