Fast Temporal Wavelet Graph Neural Networks
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
Duc Thien Nguyen, Manh Duc Tuan Nguyen, Truong Son Hy, Risi Kondor
# Fast Temporal Wavelet Graph Neural Networks

###### 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](https://github.com/HySonLab/TWGNN).

Machine Learning, ICML

## 1 Introduction

Time series modeling has been a quest in a wide range of academic fields and industrial applications, including neuroscience (Pourahmadi and Noorbalooloch, 2016) and traffic modeling (Li et al., 2018). Traditionally, model-based approaches such as autoregressive (AR) and Support Vector Regression (Smola and Scholkopf, 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 (Goktas et al., 2020) or GAN-based graph autoencoder (Gurter 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 and 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 and 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 and 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

### 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 \(\mathbf{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},\ldots,u_{r}\}\) in the form

\[\mathbf{A}=\sum_{i=1}^{r}\lambda_{i}u_{i}u_{i}^{T}\,,\]

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

\[\mathbf{A}=\mathbf{U}^{T}\mathbf{H}\mathbf{U},\]
where \(\mathbf{U}\) is an orthogonal matrix and \(\mathbf{H}\) is a diagonal matrix with the eigenvalues of \(\mathbf{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

\[\mathbf{A}=\mathbf{U}_{1}^{T}\mathbf{U}_{2}^{T}\ldots\mathbf{U}_{L}^{T}\mathbf{HU}_{L}\ldots\mathbf{U} _{2}\mathbf{U}_{1},\]

where \(\mathbf{H}\) is close to diagonal and \(\mathbf{U}_{1},\ldots,\mathbf{U}_{L}\) are sparse orthogonal matrices with the following constraints:

* Each \(\mathbf{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 \(\mathbf{U}_{\ell}\) are a subset of \(\mathbb{S}_{\ell}\).
* \(\mathbf{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 \(\mathbf{U}_{\ell}\) is a 2-point rotation (i.e. \(k=2\)).

**Definition 3.1**.: We say that \(\mathbf{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

\[\mathbf{U}=\mathbf{I}_{n-k}\oplus_{(i_{1},\cdots,i_{k})}\mathbf{O}\]

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

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

**Definition 3.3**.: Given an appropriate subset \(\mathbb{O}\) of the group \(\mathbb{SO}(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\cdots\geq d_{L}\geq 1\), a **Multiresolution Matrix Factorization (MMF)** of a symmetric matrix \(\mathbf{A}\in\mathbb{R}^{n\times n}\) over \(\mathbb{O}\) is a factorization of the form

\[\mathbf{A}=\mathbf{U}_{1}^{T}\mathbf{U}_{2}^{T}\ldots\mathbf{U}_{L}^{T}\mathbf{HU}_{L}\ldots\mathbf{U} _{2}\mathbf{U}_{1}, \tag{2}\]

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

**Definition 3.4**.: We say that a symmetric matrix \(\mathbf{A}\in\mathbb{R}^{n\times n}\) is **fully multiresolution factorizable** over \(\mathbb{O}\subset\mathbb{SO}(n)\) with \((d_{1},\ldots,d_{L})\) if it has a decomposition of the form described in Def. 3.3.

### 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

\[\cdots\subset\mathbb{V}_{-1}\subset\mathbb{V}_{0}\subset\mathbb{V}_{1}\subset \mathbb{V}_{2}\subset\ldots \tag{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 \(\mathbf{U}_{\ell}\colon\mathbb{V}_{\ell-1}\to\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}^{\dim(\mathbb{V}_{\ell-1})}[\mathbf{U}_{\ell}]_{m,i}\phi _{i}^{\ell-1},\]

\[\psi_{m}^{\ell}=\sum_{i=1}^{\dim(\mathbb{V}_{\ell-1})}[\mathbf{U}_{\ell}]_{m+\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 \(\dim(\mathbb{V}_{\ell})=d_{\ell}=|\mathbb{S}_{\ell}|\). In the \(\Phi_{1}\cup\Psi_{1}\) basis, \(\mathbf{A}\) compresses into \(\mathbf{A}_{1}=\mathbf{U}_{1}\mathbf{A}\mathbf{U}_{1}^{T}\). In the \(\Phi_{2}\cup\Psi_{2}\cup\Psi_{1}\) basis, it becomes \(\mathbf{A}_{2}=\mathbf{U}_{2}\mathbf{U}_{1}\mathbf{A}\mathbf{U}_{1}^{T}\mathbf{U}_{2}^{T}\), and so on. Finally, in the \(\Phi_{L}\cup\Psi_{L}\cup\cdots\cup\Psi_{1}\) basis, it takes on the form \(\mathbf{A}_{L}=\mathbf{H}=\mathbf{U}_{L}\ldots\mathbf{U}_{2}\mathbf{U}_{1}\mathbf{A}\mathbf{U}_{1}^{T}\mathbf{U} _{2}^{T}\ldots\mathbf{U}_{L}^{T}\) that consists of four distinct blocks (supposingly that we permute the rows/columns accordingly):

\[\mathbf{H}=\begin{pmatrix}\mathbf{H}_{\Phi,\Phi}&\mathbf{H}_{\Psi,\Psi}\\ \mathbf{H}_{\Psi,\Phi}&\mathbf{H}_{\Psi,\Psi}\end{pmatrix},\]

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

\[\mathbf{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_{L}}c_{m}^{\ell}\psi_{m}^{\ell}{\psi_{m}^{\ell}}^{T},\]
where \(h_{i,j}\) coefficients are the entries of the \(\mathbf{H}_{\Phi,\Phi}\) block, and \(c^{\ell}_{m}=\langle\psi^{\ell}_{m},\mathbf{A}\psi^{\ell}_{m}\rangle\) wavelet frequencies are the diagonal elements of the \(\mathbf{H}_{\Psi,\Psi}\) block.

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

\[\mathbf{A}\to\mathbf{U}_{1}\mathbf{A}\mathbf{U}_{1}^{T}\to\cdots\to\mathbf{U}_{L}\ldots\mathbf{U}_{1} \mathbf{A}\mathbf{U}_{1}^{T}\ldots\mathbf{U}_{L}^{T}\doteq\mathbf{H}, \tag{4}\]

so the corresponding matrix factorization must be a multi-level factorization

\[\mathbf{A}\approx\mathbf{U}_{1}^{T}\mathbf{U}_{2}^{T}\ldots\mathbf{U}_{\ell}^{T}\mathbf{H}\mathbf{U}_ {\ell}\ldots\mathbf{U}_{2}\mathbf{U}_{1}. \tag{5}\]

### MMF optimization problem

Finding the best MMF factorization to a symmetric matrix \(\mathbf{A}\) involves solving

\[\min_{\begin{subarray}{c}\mathbb{S}_{L}\in\mathbb{S}\subseteq\mathbb{S}_{1} \subseteq\mathbb{S}_{0}=[n]\\ \mathbf{H}\in\mathbb{H}_{n}^{\mathbb{S}_{L}};\,\mathbf{U}_{1},\ldots,\mathbf{U}_{L}\in \mathbb{O}\end{subarray}}\|\mathbf{A}-\mathbf{U}_{1}^{T}\ldots\mathbf{U}_{L}^{T}\mathbf{H}\mathbf{ U}_{L}\ldots\mathbf{U}_{1}\|. \tag{6}\]

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

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

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

\[\min_{\begin{subarray}{c}\mathbb{S}_{L}\subseteq\cdots\subseteq\mathbb{S}_{1} \subseteq\mathbb{S}_{0}=[n]\\ \end{subarray}}\min_{\mathbf{U}_{1},\ldots,\mathbf{U}_{L}\in\mathbb{O}}\|\mathbf{U}_{L} \ldots\mathbf{U}_{1}\mathbf{A}\mathbf{U}_{1}^{T}\ldots\mathbf{U}_{L}^{T}\|_{\text{resi}}^{2}. \tag{8}\]

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

* For the inner optimization, the variables (i.e. Givens rotations \(\mathbf{U}_{1},\ldots,\mathbf{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.

### Stiefel manifold optimization

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

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

where \(\mathbf{I}_{p}\) is the identity matrix and \(\mathcal{F}(\mathbf{X}):\mathbb{R}^{n\times p}\to\mathbb{R}\) is a differentiable function. The feasible set \(\mathcal{V}_{p}(\mathbb{R}^{n})=\{\mathbf{X}\in\mathbb{R}^{n\times p}:\mathbf{X}^{T} \mathbf{X}=\mathbf{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_{\mathbf{X}_{1}\in\mathcal{V}_{p_{1}}(\mathbb{R}^{n_{1}}),\ldots,\mathbf{X}_{q} \in\mathcal{V}_{p_{q}}(\mathbb{R}^{n_{q}})}\mathcal{F}(\mathbf{X}_{1},\ldots,\mathbf{ X}_{q}) \tag{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 \(\mathbf{X}_{\ell}=\mathbf{O}_{\ell}\in\mathbb{R}^{k\times k}\), where \(\mathbf{U}_{\ell}=\mathbf{I}_{n-k}\oplus_{\mathbb{I}_{\ell}}\mathbf{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}(\mathbf{O}_{1},\ldots,\mathbf{O}_{L})=\|\mathbf{U}_{L}\ldots\mathbf{U}_{1}\mathbf{A} \mathbf{U}_{1}^{T}\ldots\mathbf{U}_{L}^{T}\|_{\text{resi}}^{2}. \tag{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

### 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 \(\mathbf{A}\in\mathbb{R}^{n\times n}\). Without the loss of generality, we assume that \(\mathbf{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 \(\mathbf{A}_{i,j}\). Given a graph signal \(\mathbf{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),\quad\text{for each}\;\;v\in V,\]

where \(\alpha_{m}^{\ell}=(f,\psi_{m}^{\ell})\) and \(\beta_{m}=(f,\phi_{m}^{L})\) 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},\ldots,\psi^{L}\}\), each corresponds to a resolution (see Figure 6 for visualizations).
* The rows of \(\mathbf{H}=\mathbf{A}_{L}\) make exactly \(n-L\) father wavelets \(\overline{\phi}=\{\phi_{m}^{L}=\mathbf{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 \(\mathbf{f}^{(k-1)}\) of size \(|V|\times F_{k-1}\) into an output \(\mathbf{f}^{(k)}\) of size \(|V|\times F_{k}\) as

\[\mathbf{f}^{(k)}_{:,j}=\sigma\Bigg{(}\mathbf{W}\sum_{i=1}^{F_{k-1}}\mathbf{g}^{(k)}_{i,j} \mathbf{W}^{T}\mathbf{f}^{(k-1)}_{:,i}\Bigg{)}\quad\text{for}\;\;j=1,\ldots,F_{k}, \tag{12}\]

where \(\mathbf{W}\) is our wavelet basis matrix as we concatenate \(\overline{\phi}\) and \(\overline{\psi}\) column-by-column, \(\mathbf{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 \(\mathbf{f}\) into the spectral domain (i.e. \(\hat{\mathbf{f}}=\mathbf{W}^{T}\mathbf{f}\) is the forward transform and \(\mathbf{f}=\mathbf{W}\hat{\mathbf{f}}\) is the inverse transform), then a learnable filter \(\mathbf{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 \(\mathbf{W}\) is _sparse_, both the wavelet transform and its inverse transform can be implemented efficiently via sparse matrix multiplication.

### 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=(\mathbf{X},\mathbf{A})\) represent an undirected graph, where \(\mathbf{X}=[\mathbf{x}_{1}^{T},\ldots,\mathbf{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{\mathbf{A}}\in\mathbb{R}^{N\times N}\) in which \(\sum_{j}\tilde{\mathbf{A}}_{ij}=1\forall i\). In the simplest case, when \(\tilde{\mathbf{L}}=\mathbf{I}-\tilde{\mathbf{A}}\) is the normalized random walk matrix, the diffusion process on graph is governed by the following

Figure 3: Architecture for the Wavelet Convolutional Gated Recurrent Unit. **WC**: graph wavelet convolution given MMF’s wavelet basis

Figure 2: Architecture of Fast Temporal Wavelet Neural Network. **WC:** graph wavelet convolution given MMF’s wavelet basis.
equation (Chamberlain et al., 2021):

\[\frac{\mathrm{d}\mathbf{X}(t)}{\mathrm{d}t}=(\tilde{\mathbf{A}}-\mathbf{I})\mathbf{X}(t) \tag{13}\]

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

\[\mathbf{X}(k) =\mathbf{X}(k-1)+(\tilde{\mathbf{A}}-\mathbf{I})\mathbf{X}(k-1)\] \[=\mathbf{X}(k-1)-\tilde{\mathbf{L}}\mathbf{X}(k-1)\] \[=\tilde{\mathbf{A}}\mathbf{X}(k-1)\] \[=\tilde{\mathbf{A}}^{k}\mathbf{X}(0) \tag{14}\]

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

\[\mathbf{X}_{:,d}\star_{\mathcal{G}}f_{\mathbf{\theta}}=\sum_{k=0}^{K-1}\theta_{k} \tilde{\mathbf{A}}^{k}\mathbf{X}_{:,d}\quad\forall d\in\{1,\dots,D\} \tag{15}\]

where \(\mathbf{\Theta}\in\mathbb{R}^{K\times 2}\) are the parameters for the filter.

Temporal Dependency ModelThe 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).

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

where \(\mathbf{X}(t),\mathbf{H}(t)\) denote the input and output of at time \(t\), while \(\mathbf{r}^{(t)},\mathbf{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 \(\mathbf{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 \(\mathbf{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 \(\mathbf{W}\), that can be computed as product \(\mathbf{U}_{1}\mathbf{U}_{2}\cdots\mathbf{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 and 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 hierarchical 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 \(\{\mathbf{U}_{\ell}\}_{\ell=1}^{L}\) and the wavelet basis \(\mathbf{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](https://github.com/HySonLab/TWGNN). Implementation of Multiresolution Matrix Factorization and graph wavelet computation (Hy and Kondor, 2022) is publicly available at [https://github.com/risilab/Learnable_MMF](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 METRLA 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. 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.

Adjacency matrixAccording to DCRNN (Li et al., 2018), the traffic sensor network is expressed by an adja
cency 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\begin{cases}\exp\!\left(-\frac{\text{dist}(v_{i},v_{j})}{ \sigma^{2}}\right)\!,&\text{dist}(v_{i},v_{j})\leq k\\ 0,&\text{otherwise}\end{cases}\,, \tag{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 and Roweis, 2003). In particular, for the matrix data \(\mathbf{X}=[\mathbf{x}_{1},\dots,\mathbf{x}_{N}]\in\mathbb{R}^{T\times N}\) where \(\mathbf{x}_{i}\) denotes the time series data of node \(i\) for \(i\in\{1,\dots,N\}\), an adjacency matrix \(\mathbf{A}\) is identified to gather all the coefficients of the affine dependencies among \(\{\mathbf{x}_{i}\}_{i=1}^{N}\) by solving the following optimization problem

\[\mathbf{A}\coloneqq\arg\min_{\mathbf{\hat{A}}\in\mathbb{R}^{N\times N}} \left\|\mathbf{X}-\mathbf{X}\hat{\mathbf{A}}^{T}\right\|_{\text{F}}^{2}+\lambda_{A} \big{\|}\hat{\mathbf{A}}\big{\|}_{1}\] \[\text{s.to}\quad\mathbf{1}_{N}^{T}\hat{\mathbf{A}}=\mathbf{1}_{N}^{T}\,,\quad \text{diag}(\hat{\mathbf{A}})=\mathbf{0}\,, \tag{17}\]

where the constraint \(\mathbf{1}_{N}^{T}\hat{\mathbf{A}}=\mathbf{1}_{N}^{T}\) realizes the affine compositions, while \(\text{diag}(\hat{\mathbf{A}})=\mathbf{0}\) excludes the trivial solution \(\hat{\mathbf{A}}=\mathbf{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}},\dots,x_{j_{k}}\}\), thus \(\hat{A}\) is regularized by the \(l_{1}\)-norm loss \(\big{\|}\hat{\mathbf{A}}\big{\|}_{1}\) to be sparse. Task (17) is a composite convex minimization problem with affine constraints, hence can be solved by (Slavakis and Yamada, 2018).

### 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 (CalorTrans) Performance Measurement System (PeMS).

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{\mathbf{A}}\coloneqq\frac{1}{2}(\mathbf{A}+\mathbf{A}^{T})\) is taken to compute the wavelet basis matrix \(\mathbf{W}\) following Sec. 3.2.

\begin{table}
\begin{tabular}{c||c|c|c c c c c c c} \hline \hline Dataset & \(T\) & Metric & HA & ARIMA\({}_{kal}\) & VAR & SVR & FNN & FC-LSTM & DCRNN & FTWGNN \\ \hline \multirow{6}{*}{METR-LA} & \multirow{3}{*}{15 min} & MAE & 4.16 & 3.99 & 4.42 & 3.99 & 3.99 & 3.44 & 2.77 & **2.70** \\  & & RMSE & 7.80 & 8.21 & 7.89 & 8.45 & 7.94 & 6.30 & 5.38 & **5.15** \\  & & MAPE & 13.00\% & 9.6\% & 10.2\% & 9.3\% & 9.9\% & 9.6\% & 7.3\% & **6.81\%** \\ \cline{2-10}  & \multirow{3}{*}{30 min} & MAE & 4.16 & 5.15 & 5.41 & 5.05 & 4.23 & 3.77 & 3.15 & **3.0** \\  & & 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\%** \\ \cline{2-10}  & \multirow{3}{*}{60 min} & MAE & 4.16 & 6.90 & 6.52 & 6.72 & 4.49 & 4.37 & 3.60 & **3.42** \\  & & 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\%** \\ \hline \hline \multirow{6}{*}{PEMS-BAY} & \multirow{3}{*}{15 min} & MAE & 2.88 & 1.62 & 1.74 & 1.85 & 2.20 & 2.05 & 1.38 & **1.14** \\  & & 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\%** \\ \cline{2-10}  & \multirow{3}{*}{30 min} & MAE & 2.88 & 2.33 & 2.32 & 2.48 & 2.30 & 2.20 & 1.74 & **1.5** \\  & & MAPE & 6.8\% & 5.4\% & 5.0\% & 5.5\% & 5.43\% & 5.2\% & 3.9\% & **3.15\%** \\ \cline{2-10}  & \multirow{3}{*}{60 min} & MAE & 2.88 & 3.38 & 2.93 & 3.28 & 2.46 & 2.37 & 2.07 & **1.79** \\  & & RMSE & 5.59 & 6.5 & 5.44 & 7.08 & 4.98 & 4.96 & 4.74 & **3.99** \\ \cline{2-10}  & \multirow{3}{*}{60 min} & MAE & 2.88 & 3.37 & 6.5\% & 8.0\% & 5.89\% & 5.7\% & 4.9\% & **4.14\%** \\ \cline{1-1} \cline{2-10}  & \multirow{3}{*}{30 min} & MAE & 2.88 & 2.33 & 2.32 & 2.48 & 2.30 & 2.20 & 1.74 & **1.5** \\ \cline{1-1}  & & MAPE & 6.8\% & 5.4\% & 5.0\% & 5.5\% & 5.43\% & 5.2\% & 3.9\% & **3.15\%** \\ \cline{1-1} \cline{2-10}  & \multirow{3}{*}{60 min} & MAE & 2.88 & 3.38 & 2.93 & 3.28 & 2.46 & 2.37 & 2.07 & **1.79** \\ \cline{1-1}  & & RMSE & 5.59 & 6.5 & 5.44 & 7.08 & 4.98 & 4.96 & 4.74 & **3.99** \\ \cline{1-1} \cline{2-10}  & \multirow{3}{*}{60 min} & MAE & 6.8\% & 8.3\% & 6.5\% & 8.0\% & 5.89\% & 5.7\% & 4.9\% & **4.14\%** \\ \hline \hline \end{tabular}
\end{table}
Table 2: Training time/epoch between DCRNN and FTWGNN.

\begin{table}
\begin{tabular}{c||c|c c|c} \hline \hline Dataset & \(T\) & DCRNN & FTWGNN & Speedup \\ \hline \multirow{3}{*}{METR-LA} & 15 min & 350s & **217s** & 1.61x \\ \cline{2-5}  & 30 min & 620s & **163s** & 3.80x \\ \cline{2-5}  & 60 min & 1800s & **136s** & 13.2x \\ \cline{2-5}  & 15 min & 427s & **150s** & 2.84x \\ \cline{2-5}  & 30 min & 900s & **173s** & 5.20x \\ \cline{2-5}  & 60 min & 1800s
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.

### Brain networks

_Annotated Joints in Long-term Electrocorticography for 12 human participants_ (AJILE12), publicly available at (Peterson et al., 2022), records intracranial neural activity via the invasive ECoG, which involves implanting electrodes directly under the skull (Peterson et al., 2022). 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 3(a)). The adjacency matrix \(\mathbf{A}\) is obtained by solving task (17), then the wavelet basis matrix \(\mathbf{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 3(b) 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 \(\mathbf{W}\) for each dataset. It is important to note that \(\mathbf{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 Net

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

\begin{table}
\begin{tabular}{c||c|c|c c c c c c c} \hline Dataset & \(T\) & Metric & HA & VAR & LR & SVR & LSTM & DCRNN & FTWGNN \\ \hline \multirow{8}{*}{AJILE12} & \multirow{4}{*}{1 sec} & MAE & 0.88 & 0.16 & 0.27 & 0.27 & 0.07 & 0.05 & **0.03** \\  & & 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\%** \\ \cline{2-10}  & \multirow{4}{*}{5 sec} & MAE & 0.88 & 0.66 & 0.69 & 0.69 & 0.39 & 0.16 & **0.11** \\  & & RMSE & 1.23 & 0.96 & 0.92 & 0.93 & 0.52 & 0.24 & **0.15** \\  & & MAPE & 320\% & 221\% & 376\% & 339\% & 147\% & 64\% & **57\%** \\ \cline{2-10}  & \multirow{4}{*}{1 sec} & MAE & 0.88 & 0.82 & 0.86 & 0.86 & 0.87 & 0.78 & **0.70** \\  & & RMSE & 1.23 & 1.15 & 1.13 & 1.13 & 1.14 & 1.01 & **0.93** \\ \cline{1-1}  & & MAPE & 320\% & 320\% & 448\% & 479\% & 330\% & 294\% & **254\%** \\ \hline \end{tabular}
\end{table}
Table 3: Performance comparison on ECoG signals forecast.
work_ (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.
