CayleyNets: Graph Convolutional Neural Networks with Complex Rational Spectral Filters

TL;DR

CayleyNets introduce a novel spectral graph convolutional architecture using rational complex functions, enabling efficient, localized filtering that improves performance across various graph-based learning tasks.

Contribution

The paper presents Cayley polynomials for spectral filtering, offering a flexible, scalable, and effective method for deep learning on graphs with diverse Laplacian constructions.

Findings

Outperforms existing spectral graph convolutional methods in classification tasks

Filters are localized and scale linearly with data size for sparse graphs

Effective across multiple graph learning applications

Abstract

The rise of graph-structured data such as social networks, regulatory networks, citation graphs, and functional brain networks, in combination with resounding success of deep learning in various applications, has brought the interest in generalizing deep learning models to non-Euclidean domains. In this paper, we introduce a new spectral domain convolutional architecture for deep learning on graphs. The core ingredient of our model is a new class of parametric rational complex functions (Cayley polynomials) allowing to efficiently compute spectral filters on graphs that specialize on frequency bands of interest. Our model generates rich spectral filters that are localized in space, scales linearly with the size of the input data for sparsely-connected graphs, and can handle different constructions of Laplacian operators. Extensive experimental results show the superior performance of our approach, in comparison to other spectral domain convolutional architectures, on spectral image classification, community detection, vertex classification and matrix completion tasks.

Tables

Table 1. Table 1: Test accuracy obtained with different methods on the MNIST dataset.

Model	Order	Accuracy	#Params
LeNet5	-	99.33%	1.66M
ChebNet	25	99.14%	1.66M
CayleyNet	12	99.18%	1.66M

Table 2. Table 2: Test accuracy of different methods on the standard split of the CORA dataset.

Method	Accuracy	#Params
DCNN [1]	72.3 $\pm$ 0.8 %	23K
CayleyNet64 features	81.0 $\pm$ 0.5 %	92K
GCN [23]	81.6 $\pm$ 0.4 %	23K
MoNet [30]	81.7 $\pm$ 0.5 %	23K
CayleyNet16 features	81.9 $\pm$ 0.7 %	23K
GAT [39]	83.0 $\pm$ 0.7 %	92K

Table 3. Table 4: Performance (RMSE) of different matrix completion methods on the MovieLens dataset.

Method	RMSE
MC [7]	0.973
IMC [20, 40]	1.653
GMC [21]	0.996
GRALS [32]	0.945
sRGCNNCheby,r=4 [31]	0.929
sRGCNNCheby,r=8 [31]	0.925
sRGCNNCayley	0.922

Equations

$${\mathbf{f}}_l^{\mathrm{out}}=\xi \left(\sum _{l^{\prime }=1}^p{\mathbf{\Phi }}_k{\hat{\mathbf{G}}}_{l,l^{\prime }}{\mathbf{\Phi }}_k^{\top }{\mathbf{f}}_{l^{\prime }}^{\mathrm{in}}\right).$$

$$g_{\mathbf{\alpha }}(\tilde{\lambda })=\sum _{j=0}^r{\alpha }_j{T}_j(\tilde{\lambda })\text{vspace}-0.5mm$$

$$T_n(\cos (\theta ))=\cos (n\theta ),$$

$$g_{\mathbf{c},h}(\lambda)=c_{0}+2{\rm Re}\Big{\{}\sum_{j=1}^{r}c_{j}(h\lambda-i)^{j}(h\lambda+i)^{-j}\Big{\}}$$

$$\mathbf{Gf}=g_{\mathbf{c},h}(\mathbf{\Delta })\mathbf{f}=c_0\mathbf{f}+2\mathrm{Re}\{\sum _{j=1}^r{c}_j(h\mathbf{\Delta }-i\mathbf{I}{)}^j(h\mathbf{\Delta }+i\mathbf{I}{)}^{-j}\mathbf{f}\},$$

$$\mathbf{G}=c_0\mathbf{I}+\sum _{j=1}^r{c}_j{\mathcal{C}}^j(h\mathbf{\Delta })+\overline{c_j}{\mathcal{C}}^{-j}(h\mathbf{\Delta }).$$

$$\mathcal{C}^{j}(h\boldsymbol{\Delta})=\boldsymbol{\Phi}\mathrm{diag}\big{(}\big{[}\mathcal{C}(h\lambda_{1})\big{]}^{j},\ldots,\big{[}\mathcal{C}(h\lambda_{n})\big{]}^{j}\big{)}\boldsymbol{\Phi}^{\top}.$$

$${\mathbf{y}}_0=\mathbf{f},(h\mathbf{\Delta }+i\mathbf{I}){\mathbf{y}}_j=(h\mathbf{\Delta }-i\mathbf{I}){\mathbf{y}}_{j-1},j=1,\dots ,r.$$

$$\begin{array}{rl}{\tilde{\mathbf{y}}}_j^{(k+1)}=& \mathbf{J}{\tilde{\mathbf{y}}}_j^{(k)}+{\mathbf{b}}_j\\ {\mathbf{b}}_j=& (\mathrm{Diag}(h\mathbf{\Delta }+i\mathbf{I}){)}^{-1}(h\mathbf{\Delta }-i\mathbf{I}){\tilde{\mathbf{y}}}_{j-1},\end{array}$$

$$\mathcal{C}(h\mathbf{\Delta }){\tilde{\mathbf{y}}}_{j-1}-{\tilde{\mathbf{y}}}_j={\mathbf{J}}^{K+1}\mathcal{C}(h\mathbf{\Delta }){\tilde{\mathbf{y}}}_{j-1}$$

$$e_j=\frac{{\Vert {\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}-{\tilde{\mathbf{y}}}_j\Vert }_2}{{\Vert {\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}\Vert }_2}.$$

$$\begin{array}{rl}{e}_j\le & \frac{{\Vert {\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}-\mathcal{C}(h\mathbf{\Delta }){\tilde{\mathbf{y}}}_{j-1}\Vert }_2}{{\Vert {\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}\Vert }_2}+\frac{{\Vert \mathcal{C}(h\mathbf{\Delta }){\tilde{\mathbf{y}}}_{j-1}-{\tilde{\mathbf{y}}}_j\Vert }_2}{{\Vert {\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}\Vert }_2}\\ =& \frac{{\Vert {\mathcal{C}}^{j-1}(h\mathbf{\Delta })\mathbf{f}-{\tilde{\mathbf{y}}}_{j-1}\Vert }_2}{{\Vert {\mathcal{C}}^{j-1}(h\mathbf{\Delta })\mathbf{f}\Vert }_2}+\frac{{\Vert {\mathbf{J}}^{K+1}\mathcal{C}(h\mathbf{\Delta }){\tilde{\mathbf{y}}}_{j-1}\Vert }_2}{{\Vert \mathbf{f}\Vert }_2}\\ \le & e_{j-1}+{\Vert {\mathbf{J}}^{K+1}\Vert }_2\frac{{\Vert \mathcal{C}(h\mathbf{\Delta }){\tilde{\mathbf{y}}}_{j-1}\Vert }_2}{{\Vert \mathbf{f}\Vert }_2}\\ =& e_{j-1}+{\Vert {\mathbf{J}}^{K+1}\Vert }_2\frac{{\Vert {\tilde{\mathbf{y}}}_{j-1}\Vert }_2}{{\Vert \mathbf{f}\Vert }_2}\\ \le & e_{j-1}+{\Vert {\mathbf{J}}^{K+1}\Vert }_2(1+e_{j-1})\end{array}$$

$$\begin{array}{rl}{\Vert {\tilde{\mathbf{y}}}_{j-1}\Vert }_2\le & {\Vert {\mathcal{C}}^{j-1}(h\mathbf{\Delta })\mathbf{f}\Vert }_2+{\Vert {\mathcal{C}}^{j-1}(h\mathbf{\Delta })\mathbf{f}-{\tilde{\mathbf{y}}}_{j-1}\Vert }_2\\ =& {\Vert \mathbf{f}\Vert }_2+{\Vert \mathbf{f}\Vert }_2{e}_{j-1}.\end{array}$$

$$\begin{array}{rl}{e}_j\le & e_{j-1}+\sqrt{n}{\Vert \mathbf{J}\Vert }_{\infty }^{K+1}(1+e_{j-1})\\ =& (1+\sqrt{n}{\kappa }^{K+1})e_{j-1}+\sqrt{n}{\kappa }^{K+1}.\end{array}$$

$$e_j\le (1+\sqrt{n}{\kappa }^{K+1}{)}^j-1=j\sqrt{n}{\kappa }^{K+1}+\mathcal{O}({\kappa }^{2K+2}).$$

$$e_j\le j\sqrt{n}{\kappa }^{K+1}.$$

$$\mathbf{J}=-(hd\mathbf{I}+i\mathbf{I}{)}^{-1}h(\mathbf{\Delta }-d\mathbf{I})=\frac{h}{hd+i}(d\mathbf{I}-\mathbf{\Delta })=\frac{h}{hd+i}\mathbf{W}.$$

$$e_j\le e_{j-1}+{\Vert \mathbf{J}\Vert }_2^{K+1}(1+e_{j-1})=e_{j-1}+{\kappa }^{K+1}(1+e_{j-1}).$$

$$\begin{array}{rl}\frac{{\Vert \mathbf{Gf}-\widetilde{\mathbf{Gf}}\Vert }_2}{{\Vert \mathbf{f}\Vert }_2}\le & 2\sum _{j=1}^r\left|c_j\right|\frac{{\Vert {\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}-{\tilde{\mathbf{y}}}_j\Vert }_2}{{\Vert {\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}\Vert }_2}\\ =& 2\sum _{j=1}^r\left|c_j\right|e_j\le 2\sum _{j=1}^r{M}_j\left|c_j\right|{\kappa }^{K+1}.\end{array}$$

$$\begin{gathered} {\Vert \mathbf{G}{\mathbf{\delta }}_m-\mathbf{G}{\mathbf{\delta }}_m{|}_{{\mathcal{N}}_{l,m}}\Vert }_2 \\ \le {\Vert \mathbf{G}{\mathbf{\delta }}_m-{\widetilde{\mathbf{G}\mathbf{\delta }}}_m\Vert }_2+{\Vert {\widetilde{\mathbf{G}\mathbf{\delta }}}_m-\mathbf{G}{\mathbf{\delta }}_m{|}_{{\mathcal{N}}_{l,m}}\Vert }_2 \\ \le {\Vert \mathbf{G}{\mathbf{\delta }}_m-{\widetilde{\mathbf{G}\mathbf{\delta }}}_m\Vert }_2+{\Vert {\widetilde{\mathbf{G}\mathbf{\delta }}}_m-\mathbf{G}{\mathbf{\delta }}_m\Vert }_2 \\ =2{\Vert \mathbf{G}{\mathbf{\delta }}_m-{\widetilde{\mathbf{G}\mathbf{\delta }}}_m\Vert }_2\le 4M{\kappa }^{K+1}{\Vert {\mathbf{\delta }}_m\Vert }_2 \\ =4M({\kappa }^{1/r}{)}^l. \end{gathered}$$

$${\mathbf{G}}_{\mathbf{c},h}\mathbf{f}=c_0\mathbf{f}+\sum _{j=1}^{r}2\mathrm{Re}\{c_j{\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}\}.$$

$$S(\mathbf{c},h)=F({\mathbf{G}}_{\mathbf{c},h}\mathbf{f})$$

$$S(\mathbf{a}+ϵ\mathbf{\eta })-S(\mathbf{a})=ϵ⟨D[S;\mathbf{a}],\mathbf{\eta }⟩+o(ϵ)$$

$$\frac{{\mathbf{G}}_{\mathbf{c}+ϵ\mathbf{\eta },h}\mathbf{f}-{\mathbf{G}}_{\mathbf{c},h}\mathbf{f}}{ϵ}={\mathbf{G}}_{\mathbf{\eta },h}\mathbf{f}$$

$$={\rm Re}\big{\{}\eta_{0}\mathbf{f}+\sum_{j=1}^{r}2\eta_{j}\mathcal{C}^{j}(h\boldsymbol{\Delta})\mathbf{f}\big{\}}.$$

$$D_{\mathbf{c}}[\mathbf{G}_{\mathbf{c},h}\mathbf{f};(h,\mathbf{c})]_{j}=\left\{\begin{array}[]{ccc}\mathbf{f}&,&j=0\\ 2\mathcal{C}^{j}(h\boldsymbol{\Delta})\mathbf{f}&,&{\rm else}\end{array}\right.$$

$$\frac{S(\mathbf{c}+ϵ\mathbf{\eta })-S(\mathbf{c})}{ϵ}\approx \nabla F({\mathbf{G}}_{\mathbf{c},h}\mathbf{f})\cdot \mathrm{Re}\left\{D_{\mathbf{c}}[{\mathbf{G}}_{\mathbf{c},h}\mathbf{f}]\mathbf{\eta }\right\}$$

$$\begin{gathered} \frac{S(\mathbf{c}+ϵ\eta )-S(\mathbf{c})}{ϵ} \\ \approx \mathrm{Re}\left\{\nabla F\cdot \left({\eta }_0\mathbf{f}+\sum _{j=1}^{r}2{\eta }_j{\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}\right)\right\} \\ =\mathrm{Re}\left\{[\nabla F\cdot \mathbf{f}]{\eta }_0+\sum _{j=1}^{r}2[\nabla F\cdot {\mathcal{C}}^j(h\mathbf{\Delta })\mathbf{f}]{\eta }_j\right\} \\ =\mathrm{Re}⟨\nabla F\cdot D_{\mathbf{c}}[{\mathbf{G}}_{\mathbf{c},h}\mathbf{f};(h,\mathbf{c})],\eta ⟩. \end{gathered}$$

$$\nabla S=\nabla F\cdot D[p;(\mathbf{c},f)].$$

$$\frac{\partial }{\partial h}{g}_{\mathbf{c},h}(\lambda )=\sum _{j=1}^r{c}_{j}j{\mathcal{C}}^{j-1}(h\lambda ){\mathcal{C}}^{\prime }(h\lambda )\lambda -\overline{c_j}j{\mathcal{C}}^{-j-1}(h\lambda ){\mathcal{C}}^{\prime }(h\lambda )\lambda$$

Full text

1 Introduction

In many domains, one has to deal with large-scale data with underlying non-Euclidean structure. Prominent examples of such data are social networks, genetic regulatory networks, functional networks of the brain, and 3D shapes represented as discrete manifolds. The recent success of deep neural networks and, in particular, convolutional neural networks (CNNs) [25] have raised the interest in geometric deep learning techniques trying to extend these models to data residing on graphs and manifolds. In this paper we focus on spectral graph CNNs. Geometric deep learning approaches, and specifically spectral graph CNNs, have been successfully applied to computer graphics and vision [28, 2, 4, 3, 30], brain imaging [24], and drug design [12] problems, to mention a few. For a comprehensive presentation of methods and applications of deep learning on graphs and manifolds, we refer the reader to the review paper [5].

1.1 Related work

The earliest neural network formulation on graphs was proposed by [13] and [35], combining random walks with recurrent neural networks (their paper has recently enjoyed renewed interest in [26, 38]). The first CNN-type architecture on graphs was proposed by [6]. One of the key challenges of extending CNNs to graphs is the lack of vector-space structure and shift-invariance making the classical notion of convolution elusive. Bruna *et al. *formulated convolution-like operations in the spectral domain, using the graph Laplacian eigenbasis as an analogy of the Fourier transform [37]. [10] proposed an efficient filtering scheme using recurrent Chebyshev polynomials applied on the Laplacian operator. As opposed to [6], Chebyshev filters are defined as functions $\mathbb{R}\rightarrow\mathbb{R}$ applied on the spectrum, as in [37]. This makes filters learned on one graph generalizable to other graphs. [23] simplified this architecture using filters operating on 1-hop neighborhoods of the graph. [1] proposed a Diffusion CNN architecture based on powers of the degree-normalized transition matrix. [30] (and later, [15]) proposed a spatial-domain generalization of CNNs to graphs using local patch operators represented as Gaussian mixture models, showing a significant advantage of such models in generalizing across different graphs. In [31], spectral graph CNNs were extended to multiple graphs and applied to matrix completion and recommender system problems.

1.2 Main contribution

In this paper, we construct graph CNNs employing an efficient spectral filtering scheme based on the new class of Cayley polynomials that enjoys similar advantages of the Chebyshev filters [10] such as localization and linear complexity in the number of edges. The main advantage of our filters over [10] is their ability to detect narrow frequency bands of importance during training, and to specialize on them while being well-localized on the graph. We demonstrate experimentally that this affords our method greater flexibility, making it perform better than ChebNets on a broad range of graph learning problems.

1.3 Notation

We use $a,\mathbf{a}$, and $\mathbf{A}$ to denote scalars, vectors, and matrices, respectively. $\bar{z}$ denotes the conjugate of a complex number, $\mathrm{Re}\{z\}$ its real part, and $i=\sqrt{-1}$ denotes the imaginary unit. $\mathrm{diag}(a_{1},\ldots,a_{n})$ denotes an $n\times n$ diagonal matrix with diagonal elements $a_{1},\ldots,a_{n}$. $\mathrm{Diag}(\mathbf{A})=\mathrm{diag}(a_{11},\ldots,a_{nn})$ denotes an $n\times n$ diagonal matrix obtained by setting to zero the off-diagonal elements of $\mathbf{A}$. $\mathrm{Off}(\mathbf{A})=\mathbf{A}-\mathrm{Diag}(\mathbf{A})$ denotes the matrix containing only the off-diagonal elements of $\mathbf{A}$. $\mathbf{I}$ is the identity matrix and $\mathbf{A}\circ\mathbf{B}$ denotes the Hadamard (element-wise) product of matrices $\mathbf{A}$ and $\mathbf{B}$. Proofs are given in the appendix.

2 Spectral techniques for deep learning on graphs

2.1 Spectral graph theory

Let $\mathcal{G}=(\{1,\ldots,n\},\mathcal{E},\mathbf{W})$ be an undirected weighted graph, represented by a symmetric adjacency matrix $\mathbf{W}=(w_{ij})$. We define $w_{ij}=0$ if $(i,j)\notin\mathcal{E}$ and $w_{ij}>0$ if $(i,j)\in\mathcal{E}$. We denote by $\mathcal{N}_{k,m}$ the $k$-hop neighborhood of vertex $m$, containing vertices that are at most $k$ edges away from $m$. The unnormalized graph Laplacian is an $n\times n$ symmetric positive-semidefinite matrix $\boldsymbol{\Delta}_{u}=\mathbf{D}-\mathbf{W}$, where $\mathbf{D}=\mathrm{diag}(\sum_{j\neq i}w_{ij})$ is the degree matrix. The normalized graph Laplacian is defined as $\boldsymbol{\Delta}_{\mathrm{n}}=\mathbf{D}^{-1/2}\boldsymbol{\Delta}_{u}\mathbf{D}^{-1/2}=\mathbf{I}-\mathbf{D}^{-1/2}\mathbf{W}\mathbf{D}^{-1/2}$. In the following, we use the generic notation $\boldsymbol{\Delta}$ to refer to some Laplacian.

Since both normalized and unnormalized Laplacian are symmetric and positive semi-definite matrices, they admit an eigendecomposition $\boldsymbol{\Delta}=\boldsymbol{\Phi}\boldsymbol{\Lambda}\boldsymbol{\Phi}^{\top}$, where $\boldsymbol{\Phi}=(\boldsymbol{\phi}_{1},\ldots\boldsymbol{\phi}_{n})$ are the orthonormal eigenvectors and $\boldsymbol{\Lambda}=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{n})$ is the diagonal matrix of corresponding non-negative eigenvalues (spectrum) $0=\lambda_{1}\leq\lambda_{2}\leq\ldots\leq\lambda_{n}$. The eigenvectors play the role of Fourier atoms in classical harmonic analysis and the eigenvalues can be interpreted as (the square of) frequencies. Given a signal $\mathbf{f}=(f_{1},\ldots,f_{n})^{\top}$ on the vertices of graph $\mathcal{G}$, its graph Fourier transform is given by $\hat{\mathbf{f}}=\boldsymbol{\Phi}^{\top}\mathbf{f}$. Given two signals $\mathbf{f},\mathbf{g}$ on the graph, their spectral convolution can be defined as the element-wise product of the Fourier transforms, \mathbf{f}\star\mathbf{g}=\boldsymbol{\Phi}\big{(}(\boldsymbol{\Phi}^{\top}\mathbf{g})\circ(\boldsymbol{\Phi}^{\top}\mathbf{f})\big{)}=\boldsymbol{\Phi}\,\mathrm{diag}(\hat{g}_{1},\ldots,\hat{g}_{n})\hat{\mathbf{f}}, which corresponds to the property referred to as the Convolution Theorem in the Euclidean case.

2.2 Spectral CNNs

[6] used the spectral definition of convolution to generalize CNNs on graphs, with a spectral convolutional layer of the form

[TABLE]

Here the $n\times p$ and $n\times q$ matrices $\mathbf{F}^{\mathrm{in}}=(\mathbf{f}^{\mathrm{in}}_{1},\ldots,\mathbf{f}^{\mathrm{in}}_{p})$ and $\mathbf{F}^{\mathrm{out}}=(\mathbf{f}^{\mathrm{out}}_{1},\ldots,\mathbf{f}^{\mathrm{out}}_{q})$ represent respectively the $p$- and $q$-dimensional input and output signals on the vertices of the graph, $\boldsymbol{\Phi}_{k}=(\boldsymbol{\phi}_{1},\ldots,\boldsymbol{\phi}_{k})$ is an $n\times k$ matrix of the first eigenvectors, $\hat{\mathbf{G}}_{l,l^{\prime}}=\mathrm{diag}(\hat{g}_{l,l^{\prime},1},\ldots,\hat{g}_{l,l^{\prime},k})$ is a $k\times k$ diagonal matrix of spectral multipliers representing a learnable filter in the frequency domain, and $\xi$ is a nonlinearity (e.g., ReLU) applied on the vertex-wise function values. Pooling is performed by means of graph coarsening, which, given a graph with $n$ vertices, produces a graph with $n^{\prime}<n$ vertices and transfers signals from the vertices of the fine graph to those of the coarse one. Assuming $k=\mathcal{O}(n)$ Laplacian eigenvectors are used, a spectral convolutional layer requires $\mathcal{O}(pqk)=\mathcal{O}(n)$ parameters to train.

In addition, an informal approach for keeping the filters localized in the spectral domain was proposed. The idea is to learn just a few spectral coefficients of the filter, and obtain the rest using interpolation. This also keeps the number of filter parameters $\mathcal{O}(1)$. The spatial locality property simulates local receptive fields [8], and is important for the interpretability of convolutions as filters. Moreover, for spatial implementation of (1), such as ChebNet and CayleyNet, small receptive fields typically indicate sparse implementations.

This framework has several major drawbacks. First, the computation of the forward and inverse graph Fourier transforms incur expensive $\mathcal{O}(n^{2})$ multiplication by the matrices $\boldsymbol{\Phi},\boldsymbol{\Phi}^{\top}$, as there is no FFT-like algorithms on general graphs.

Second, the spectral filter coefficients are basis dependent, and consequently, a spectral CNN model learned on one graph cannot be transferred to another graph. This is as opposed to [37], where the frequency responses of the filters coefficients are represented as $\hat{g}_{i}=g(\lambda_{i})$, where $g(\lambda)$ is a smooth transfer function of frequency $\lambda$. Applying such filter to signal $\mathbf{f}$ can be expressed as $\mathbf{G}\mathbf{f}=g(\boldsymbol{\Delta})\mathbf{f}=\boldsymbol{\Phi}g(\boldsymbol{\Lambda})\boldsymbol{\Phi}^{\top}\mathbf{f}=\boldsymbol{\Phi}\,\mathrm{diag}(g(\lambda_{1}),\ldots,g(\lambda_{n}))\boldsymbol{\Phi}^{\top}\mathbf{f}$, where applying a function to a matrix is understood in the operator functional calculus sense (applying the function to the matrix eigenvalues).

It is noteworthy to mention alternative functional calculus driven approaches to define convolution. In [33] filters are defined as functions of the adjacency matrix $g(\mathbf{W})$, and in [34] the problem of ordering the eigenvalues of $\mathbf{W}$ according to a natural notion of frequency was addressed.

2.3 ChebNet

[10] used polynomial filters represented in the Chebyshev basis

[TABLE]

applied to rescaled frequency $\tilde{\lambda}\in[-1,1]$; here, $\boldsymbol{\alpha}$ is the $(r+1)$-dimensional vector of polynomial coefficients parametrizing the filter and optimized for during the training, and $T_{j}(\lambda)=2\lambda T_{j-1}(\lambda)-T_{j-2}(\lambda)$ denotes the Chebyshev polynomial of degree $j$ defined in a recursive manner with $T_{1}(\lambda)=\lambda$ and $T_{0}(\lambda)=1$. Chebyshev polynomials form an orthogonal basis for the space of polynomials of order $r$ on $[-1,1]$. Applying the filter is performed by $g_{\boldsymbol{\alpha}}(\tilde{\boldsymbol{\Delta}})\mathbf{f}$, where $\tilde{\boldsymbol{\Delta}}=2\lambda_{n}^{-1}\boldsymbol{\Delta}-\mathbf{I}$ is the rescaled Laplacian such that its eigenvalues $\tilde{\boldsymbol{\Lambda}}=2\lambda_{n}^{-1}\boldsymbol{\Lambda}-\mathbf{I}$ are in the interval $[-1,1]$.

Such an approach has several important advantages. First, since $g_{\boldsymbol{\alpha}}(\tilde{\boldsymbol{\Delta}})=\sum_{j=0}^{r}\alpha_{j}T_{j}(\tilde{\boldsymbol{\Delta}})$ contains only matrix powers, additions, and multiplications by scalar, it can be computed avoiding the explicit expensive $\mathcal{O}(n^{3})$ computation of the Laplacian eigenvectors. Furthermore, due to the recursive definition of the Chebyshev polynomials, the computation of the filter $g_{\boldsymbol{\alpha}}(\boldsymbol{\Delta})\mathbf{f}$ entails applying the Laplacian $r$ times, resulting in $\mathcal{O}(rn)$ operations assuming that the Laplacian is a sparse matrix with $\mathcal{O}(1)$ non-zero elements in each row (a valid hypothesis for most real-world graphs that are sparsely connected). From another point of view, in each multiplication by the Laplacian, neighbors in the graph exchange data, and there are overall $r$ such neighbor exchanges. Second, the number of parameters is $\mathcal{O}(1)$ as $r$ is independent of the graph size $n$. Third, since the Laplacian is a local operator affecting only 1-hop neighbors of a vertex and a polynomial of degree $r$ of the Laplacian affects only $r$-hops, the resulting filters have guaranteed spatial localization.

A key disadvantage of Chebyshev filters is the fact that using polynomials makes it hard to produce narrow-band filters, as such filters require very high order $r$, and produce unwanted non-local filters. This deficiency is especially pronounced when the Laplacian has clusters of eigenvalues concentrated around a few frequencies with large spectral gap (Figure 3, second to the last). Indeed, in ChebNets the Laplacian eigenvalues are contracted to the band $[-1,1]$, and the clusters of eigenvalues become very concentrated. Now, the frequency response of the Chebyshev filter is a polynomial in $[-1,1]$, which is unable to separate the individual eigenvalues in the clusters due to an uncertainty principle. Such a behavior is characteristic of graphs with community structures, which is very common in many real-world graphs, for instance, social networks.

Let us explain the above phenomenon more accurately. Recall that Chebyshev polynomials are given by

[TABLE]

and form an orthonormal basis of the weighted Lebesgue space L_{2}\big{(}[-1,1],\frac{1}{\sqrt{1-x^{2}}}dx\big{)}. When the variable is changed via $x=\cos(\theta)$, the Chebyshev basis maps to the cosine basis in $[0,\pi]$, and the space L_{2}\big{(}[-1,1],\frac{1}{\sqrt{1-x^{2}}}dx\big{)} maps to the space $L_{2}(0,\pi)$ with the standard Lebesgue measure. Now, suppose that we want to represent a band-pass filter on the narrow band $[\cos(b),\cos(a)]\subset[-1,1]$, with small $\cos(b-a)$. Under the change of variable, this band maps to $[a,b]$ with small $b-a=\epsilon$. In this case, since the characteristic function of $[a,b]$ is the shrinking dilation of the characteristic function of $[-\frac{1}{2},\frac{1}{2}]$ (up to translation), it’s Fourier coefficients are samples from the stretching dilation of the Fourier transform of the characteristic function of $[-\frac{1}{2},\frac{1}{2}]$. As a result, the number of coefficients required to approximate a band pass filter up to some fixed tolerance is inverse proportional to the size of the band. More generally, the number of Chebyshev coefficients required for approximating a filter having features in a given scale, is inverse proportional to the scale. When the Laplacian has a cluster of eigenvalues concentrated around one frequency, a filter that separates these eigenvalues must have features in scale proportional to the radius of the eigenvalue cluster. Therefore, the number of Chebyshev coefficients must be inverse proportional to the cluster size.

To overcome this major drawback, we need a new class of filters, that both entail $\mathcal{O}(r)$ neighbor exchanges, and are able to specialize in narrow bands in frequency.

3 Cayley filters

A key construction of this paper is a family of complex filters that enjoy the advantages of Chebyshev filters while avoiding some of their drawbacks. We define a Cayley polynomial of order $r$ to be a real-valued function with complex coefficients,

[TABLE]

where $\mathbf{c}=(c_{0},\ldots,c_{r})$ is a vector of one real coefficient and $r$ complex coefficients and $h>0$ is the spectral zoom parameter, that will be discussed later. A Cayley filter $\mathbf{G}$ is a spectral filter defined on real signals $\mathbf{f}$ by

[TABLE]

where the parameters $\mathbf{c}$ and $h$ are optimized during training. Similarly to the Chebyshev filters, Cayley filters involve basic matrix operations such as powers, additions, multiplications by scalars, and also inversions. This implies that application of the filter $\mathbf{G}\mathbf{f}$ can be performed without explicit expensive eigendecomposition of the Laplacian operator. Cayley filters are special cases of filters based on general rational functions of the Laplacian, namely ARMA filters [18][19]. For a general rational functions of the Laplacian, calculating the denominator requires a matrix inversion. When the filter is based on arbitrary coefficients, there is no guarantee that the matrix inversions are calculated stably. Guaranteeing stable inversions for arbitrary filter coefficients is important, since the coefficients follow an unknown path during training. For general ARMA filters, the filter can acquire poles arbitrarily close to the spectrum of $\boldsymbol{\Delta}$ during training, so there is no uniform analysis of convergence of the approximate inversions in the filter computation. The motivation to use the subclass of Cayley filters over general rational functions is to guarantee that inversion is uniformly stable. Namely, the number of iterations required for a given approximation error is independent of the filter coefficients, as long as the coefficients are bounded. Moreover, the lower number of parameters in Cayley polynomials in comparison to general rational functions may be beneficial for avoiding overfitting.

In the following, we show that Cayley filters are analytically well behaved; in particular, any smooth spectral filter can be represented as a Cayley polynomial, and low-order filters are localized in the spatial domain. We also discuss numerical implementation and compare Cayley and Chebyshev filters. We show that Cayley filters defined on sparse Laplacians with $\mathcal{O}(1)$ non-zero elements takes $\mathcal{O}(n)$ operations, similarly to Chebyshev filters.

3.1 Analytic properties

Cayley filters are best understood through the Cayley transform, from which their name derives. Denote by $e^{i\mathbb{R}}=\{e^{i\theta}:\theta\in\mathbb{R}\}$ the unit complex circle. The Cayley transform $\mathcal{C}(x)=\frac{x-i}{x+i}$ is a smooth bijection between $\mathbb{R}$ and $e^{i\mathbb{R}}\setminus\{1\}$. The complex matrix $\mathcal{C}(h\boldsymbol{\Delta})=(h\boldsymbol{\Delta}-i\mathbf{I})(h\boldsymbol{\Delta}+i\mathbf{I})^{-1}$ obtained by applying the Cayley transform to the scaled Laplacian $h\boldsymbol{\Delta}$ has its spectrum in $e^{i\mathbb{R}}$ and is thus unitary. Since $z^{-1}=\overline{z}$ for $z\in e^{i\mathbb{R}}$, we can write $\overline{c_{j}\mathcal{C}^{j}(h\boldsymbol{\Delta})}=\overline{c_{j}}\mathcal{C}^{-j}(h\boldsymbol{\Delta})$. Therefore, using $2{\rm Re}\{z\}=z+\overline{z}$, any Cayley filter (4) can be written as a conjugate-even Laurent polynomial w.r.t. $\mathcal{C}(h\boldsymbol{\Delta})$,

[TABLE]

Since the spectrum of $\mathcal{C}(h\boldsymbol{\Delta})$ is in $e^{i\mathbb{R}}$, the operator $\mathcal{C}^{j}(h\boldsymbol{\Delta})$ can be thought of as a multiplication by a pure harmonic in the frequency domain $e^{i\mathbb{R}}$ for any integer power $j$,

[TABLE]

A Cayley filter can be thus seen as a multiplication by a finite Fourier expansions in the frequency domain $e^{i\mathbb{R}}$. Since (5) is conjugate-even, it is a (real-valued) trigonometric polynomial.

Note that any spectral filter can be formulated as a Cayley filter. Indeed, spectral filters $g(\boldsymbol{\Delta})$ are specified by the finite sequence of values $g(\lambda_{1}),\ldots,g(\lambda_{n})$, which can be interpolated by a trigonometric polynomial. Moreover, since trigonometric polynomials are smooth, we expect low order Cayley filters to be well localized in some sense on the graph, as discussed later.

Finally, in definition (4) we use complex coefficients. If $c_{j}\in\mathbb{R}$ then (5) is an even cosine polynomial, and if $c_{j}\in i\mathbb{R}$ then (5) is an odd sine polynomial. Since the spectrum of $h\boldsymbol{\Delta}$ is in $\mathbb{R}_{+}$, it is mapped to the lower half-circle by $\mathcal{C}$, on which both cosine and sine polynomials are complete and can represent any spectral filter. However, it is beneficial to use general complex coefficients, since complex Fourier expansions are overcomplete in the lower half-circle, thus describing a larger variety of spectral filters of the same order without increasing the computational complexity of the filter.

3.2 Spectral zoom

To understand the essential role of the parameter $h$ in the Cayley filter, consider $\mathcal{C}(h\boldsymbol{\Delta})$. Multiplying $\boldsymbol{\Delta}$ by $h$ dilates its spectrum, and applying $\mathcal{C}$ on the result maps the non-negative spectrum to the complex half-circle. The greater $h$ is, the more the spectrum of $h\boldsymbol{\Delta}$ is spread apart in $\mathbb{R}_{+}$, resulting in better spacing of the smaller eigenvalues of $\mathcal{C}(h\boldsymbol{\Delta})$. On the other hand, the smaller $h$ is, the further away the high frequencies of $h\boldsymbol{\Delta}$ are from $\infty$, the better spread apart are the high frequencies of $\mathcal{C}(h\boldsymbol{\Delta})$ in $e^{i\mathbb{R}}$ (see Figure 1). Tuning the parameter $h$ allows thus to ‘zoom’ in to different parts of the spectrum, resulting in filters specialized in different frequency bands.

3.3 Numerical properties

The numerical core of the Cayley filter is the computation of $\mathcal{C}^{j}(h\boldsymbol{\Delta})\mathbf{f}$ for $j=1,\ldots,r$, performed in a sequential manner. Let $\mathbf{y}_{0},\ldots,\mathbf{y}_{r}$ denote the solutions of the following linear recursive system,

[TABLE]

Note that sequentially approximating $\mathbf{y}_{j}$ in (6) using the approximation of $\mathbf{y}_{j-1}$ in the right hand side is stable, since $\mathcal{C}(h\boldsymbol{\Delta})$ is unitary and thus has condition number $1$.

Equations (6) can be solved with matrix inversion exactly, but it costs $\mathcal{O}(n^{3})$. An alternative is to use the Jacobi method,111We remind that the Jacobi method for solving $\mathbf{A}\mathbf{x}=\mathbf{b}$ consists in decomposing $\mathbf{A}=\mathrm{Diag}(\mathbf{A})+\mathrm{Off}(\mathbf{A})$ and obtaining the solution iteratively as $\mathbf{x}^{(k+1)}=-(\mathrm{Diag}(\mathbf{A}))^{-1}\mathrm{Off}(\mathbf{A})\mathbf{x}^{(k)}+(\mathrm{Diag}(\mathbf{A}))^{-1}\mathbf{b}$.

which provides approximate solutions $\tilde{\mathbf{y}}_{j}\approx\mathbf{y}_{j}$. Let $\mathbf{J}=-(\mathrm{Diag}(h\boldsymbol{\Delta}+i\mathbf{I}))^{-1}\mathrm{Off}(h\boldsymbol{\Delta}+i\mathbf{I})$ be the Jacobi iteration matrix associated with equation (6). For the unnormalized Laplacian, $\mathbf{J}=(h\mathbf{D}+i\mathbf{I})^{-1}h\mathbf{W}$. Jacobi iterations for approximating (6) for a given $j$ have the form

[TABLE]

initialized with $\tilde{\mathbf{y}}_{j}^{(0)}=\mathbf{b}_{j}$ and terminated after $K$ iterations, yielding $\tilde{\mathbf{y}}_{j}=\tilde{\mathbf{y}}_{j}^{(K)}$. We denote $\tilde{\mathbf{y}}_{0}={\mathbf{y}}_{0}$. The application of the approximate Cayley filter is given by $\widetilde{\mathbf{G}\mathbf{f}}={c_{0}\tilde{\mathbf{y}}_{0}+2{\rm Re}}\sum_{j=1}^{r}c_{j}\tilde{\mathbf{y}}_{j}\approx\mathbf{G}\mathbf{f}$, and takes $\mathcal{O}(rKn)$ operations under the previous assumption of a sparse Laplacian. The method can be improved by normalizing $\left\|\tilde{\mathbf{y}}_{j}\right\|_{2}=\left\|\mathbf{f}\right\|_{2}$.

Next, we give an error bound for the approximate filter. For the unnormalized Laplacian, let $d=\max_{j}\{d_{j,j}\}$ and $\kappa=\left\|\mathbf{J}\right\|_{\infty}=\frac{hd}{\sqrt{h^{2}d^{2}+1}}<1$. For the normalized Laplacian, we assume that $(h\boldsymbol{\Delta}_{n}+i\mathbf{I})$ is dominant diagonal, which gives $\kappa=\left\|\mathbf{J}\right\|_{\infty}<1$.

Proposition 1**.**

*Under the above assumptions,

$\frac{\|\mathbf{G}\mathbf{f}-\widetilde{\mathbf{G}\mathbf{f}}\|_{2}}{\left\|\mathbf{f}\right\|_{2}}\leq 2M\kappa^{K}$, where $M=\sqrt{n}\sum_{j=1}^{r}j\left|c_{j}\right|$ in the general case, and $M=\sum_{j=1}^{r}j\left|c_{j}\right|$ if the graph is regular.*

Proposition 1 is pessimistic in the general case, while requires strong assumptions in the regular case. We find that in most real life situations the behavior is closer to the regular case. It also follows from Proposition 1 that smaller values of the spectral zoom $h$ result in faster convergence, giving this parameter an additional numerical role of accelerating convergence.

Last, note that a Cayley filter with Jacobi approximation, is based on powers of the Jacobi matrix ${\mathbf{J}}$, in addition to $(h\boldsymbol{\Delta}-i\mathbf{I})$. The Jacobi matrix can be viewed as a general representation matrix of the graph, replacing the standard Laplacian with a matrix that respects the connectivity of the graph (general representation matrices were considered in [18]). This is also true for $(h\boldsymbol{\Delta}-i\mathbf{I})$. In this point of view, learning $h$ is interpreted as learning a general representation matrix. Learning $h$ can be also viewed as learning a normalization of the weights of the graph. The problem of learning the topology of the graph was studied e.g. in [16].

3.4 Complexity

In practice, an accurate inversion of $(h\boldsymbol{\Delta}+i\mathbf{I})$ is not required, since the approximate inverse is combined with learned coefficients, which “compensate”, as necessary, for the inversion inaccuracy. Such behavior is well-documented in the literature in other contexts of model compression and accelerated convergence of iterative algorithms. For example, in [14], sparse signal coding are learned by unrolling iterative shrinkage algorithms (FISTA) into a neural network, where each layer emulates an iteration of the original algorithm but has extra learnable parameters. It is shown that FISTA networks with just a few layers outperform hundreds or thousands of iterations of the original algorithm thanks to the learnable parameters. The above phenomenon is also a common observation for solving sparse linear equations for compressed sensing tasks (see e.g [41]).

In a CayleyNet for a fixed graph, we fix the number of Jacobi iterations. Since the convergence rate depends on $\kappa$, that depends on the graph, different graphs may need different numbers of iterations. The convergence rate also depends on $h$. Since there is a trade-off between the spectral zoom amount $h$, and the accuracy of the approximate inversion, and since $h$ is a learnable parameter, the training finds the right balance between the spectral zoom amount and the inversion accuracy.

To formulate computational complexity results, we consider the case where the number of vertices $n$ is “big”. To formalize this, we consider a sequence of graphs indexed by $n$, and study the asymptotics as $n\rightarrow\infty$. When the graph is sampled from a continuous entity, like a manifold, this asymptotic analysis has a precise meaning. Otherwise, the asymptotic analysis is just a formal way of saying “big $n$”. For every constant of a graph, e.g $d,\kappa$, we add the subscript $n$, indicating the number of vertices of the graph. We assume that there is a global constant $C$, such that the number of edges is bounded by $Cn$. For the unnormalized Laplacian, we assume that $d_{n}$ and $h_{n}$ are bounded, which gives $\kappa_{n}<a<1$ for some $a$ independent of $n$. For the normalized Laplacian, we assume that $\kappa_{n}<a<1$. These assumptions pose regularity on the sequence of graphs. By Proposition 1, fixing the number of Jacobi iterations $K$ and the order of the filter $r$, independently of $n$, keeps the Jacobi error controlled. As a result, the number of parameters of the Cayley filters can be kept $\mathcal{O}(1)$, and for a Laplacian modeled as a sparse matrix, applying a Cayley filter on a signal takes $\mathcal{O}(n)$ operations. Indeed, the Jacobi matrix $\mathbf{J}$ has the same connectivity as the graph, including edges connecting each vertex to itself. In each Jacobi iteration, $\mathbf{J}$ is applied $K$ times, which means that vertices exchange information with their neighbors $K$ times, and one time in the initialization due to one multiplication by $(h\boldsymbol{\Delta}-iI)$. The Jacobi approximation of ${\cal C}(h{\boldsymbol{\Delta}})$ is computed $r$ times, for the $r$ coefficients of the filter, and thus overall there are $(K+1)r$ neighbor exchanges in the method. Thus, for a sparsely connected graph with $\mathcal{O}(n)$ edges, applying a Cayley filter on a signal takes $\mathcal{O}((K+1)rn)$ operations.

3.5 Localization

Unlike Chebyshev filters that have the small $r$-hop support, Cayley filters are rational functions supported on the whole graph. However, it is still true that Cayley filters are well localized on the graph. Let $\mathbf{G}$ be a Cayley filter and $\boldsymbol{\delta}_{m}$ denote a delta-function on the graph, defined as one at vertex $m$ and zero elsewhere. We show that $\mathbf{G}\boldsymbol{\delta}_{m}$ decays fast, in the following sense:

Definition 2** (Exponential decay on graphs).**

Let $\mathbf{f}$ be a signal on the vertices of graph $\mathcal{G}$, $1\leq p\leq\infty$, and $0<\epsilon<1$. Denote by $S\subseteq\{1,\ldots,n\}$ a subset of the vertices and by $S^{c}$ its complement. We say that the $L_{p}$-mass of $\mathbf{f}$ is supported in $S$ up to $\epsilon$ if $\left\|\mathbf{f}|_{S^{c}}\right\|_{p}\leq\epsilon\left\|\mathbf{f}\right\|_{p}$, where $\mathbf{f}|_{S^{c}}=(f_{l})_{l\in S^{c}}$ is the restriction of $\mathbf{f}$ to $S^{c}$. We say that $\mathbf{f}$ has (graph) exponential decay about vertex $m$, if there exists some $\gamma\in(0,1)$ and $c>0$ such that for any $k$, the $L_{p}$-mass of $\mathbf{f}$ is supported in ${\cal N}_{k,m}$ up to $c\gamma^{k}$. Here, ${\cal N}_{k,m}$ is the $k$-hop neighborhood of $m$.

Remark 3**.**

Note that Definition 2 is analogous to classical exponential decay on Euclidean space: $\left|f(x)\right|\leq R\gamma^{-x}$ iff for every ball $B_{\rho}$ of radius $\rho$ about [math], $\|f|_{B_{\rho}^{c}}\|_{\infty}\leq c\gamma^{-\rho}\left\|f\right\|_{\infty}$ with $c=\frac{R}{\left\|f\right\|_{\infty}}$.

Theorem 4**.**

Let $\mathbf{G}$ be a Cayley filter of order $r$. Then, $\mathbf{G}\boldsymbol{\delta}_{m}$ has exponential decay about $m$ in $L_{2}$, with constants $c=4M\frac{1}{\|\mathbf{G}\boldsymbol{\delta}_{m}\|_{2}}$ and $\gamma=\kappa^{1/r}$ (where $M$ and $\kappa$ are from Proposition 1).

3.6 Cayley vs Chebyshev

Below, we compare the two classes of filters:

Chebyshev as a special case of Cayley. For a regular graph with $\mathbf{D}=d\mathbf{I}$, using Jacobi inversion based on zero iteration, we get that any Cayley filter of order $r$ is a polynomial of $\boldsymbol{\Delta}$ in the monomial base \big{(}\frac{h\boldsymbol{\Delta}-i}{hd+i}\big{)}^{j}. In this situation, a Chebyshev filter, which is a real valued polynomial of $\boldsymbol{\Delta}$, is a special case of a Cayley filter.

Spectral zoom and stability. Generally, both Chebyshev polynomials and trigonometric polynomials give stable approximations, optimal for smooth functions. However, this crude statement is over-simplified. One of the drawbacks in Chebyshev filters is the fact that the spectrum of $\boldsymbol{\Delta}$ is always mapped to $[-1,1]$ in a linear manner, making it hard to specialize in small frequency bands. In Cayley filters, this problem is mitigated with the help of the spectral zoom parameter $h$. As an example, consider the community detection problem discussed in the next section. A graph with strong communities has a cluster of small eigenvalues near zero. Ideal filters $g(\boldsymbol{\Delta})$ for extracting the community information should be able to focus on this band of frequencies. Approximating such filters with Cayley polynomials, we zoom in to the band of interest by choosing the right $h$, and then project $g$ onto the space of trigonometric polynomials of order $r$, getting a good and stable approximation (Figure 3, bottom). However, if we project $g$ onto the space of Chebyshev polynomials of order $r$, the interesting part of $g$ concentrated on a small band is smoothed out and lost (Figure 3, second to the last). Thus, projections are not the right way to approximate such filters, and the stability of orthogonal polynomials cannot be invoked. On the other hand, if we want to approximate $g$ on the small band using polynomials, ignoring the behavior away from this band, the approximation will be unstable away from this band; small perturbations in $g$ will result in big perturbations in the Chebyshev filter away from the band. This is due to the fact that any polynomial diverges at infinity, and for an asymptotically small band and polynomials of fixed order, “away from the band behaves like infinity.” For this reason, we say that Cayley filters are more stable than Chebyshev filters.

Regularity. We found that in practice, low-order Cayley filters are able to model both very concentrated impulse-like filters, and wider Gabor-like filters. Cayley filters are able to achieve a wider range of filter supports with less coefficients than Chebyshev filters (Figure 2), making the Cayley class more regular than Chebyshev.

Complexity. Under the assumption of sparse Laplacians, both Cayley and Chebyshev filters incur linear complexity $\mathcal{O}(n)$. Besides, the new filters are equally simple to implement as Chebyshev filters; as seen in Eq.7, they boil down to simple sparse matrix-vector multiplications providing a GPU friendly implementation.

4 Results

4.1 Experimental settings

We test the proposed CayleyNets reproducing the experiments of [10, 23, 30] and using ChebNet [10] as our main baseline method. Pooling and graph coarsening was performed identically to [10]. The hyperparameters are identical to the original experiments, and not optimized. All the methods were implemented in TensorFlow [27]. The experiments were executed on a machine with a 3.5GHz Intel Core i7 CPU, 64GB of RAM, and NVIDIA Titan X GPU with 12GB of RAM. SGD+Momentum and Adam [22] optimization methods were used to train the models in MNIST and the rest of the experiments, respectively. Training and testing were always done on disjoint sets.

4.2 Community detection

We start with an experiment on a synthetic graph consisting of 15 communities with strong connectivity within each community and sparse connectivity across communities (Figure 3, top). Though rather simple, such a dataset allows to study the behavior of different algorithms in controlled settings. On this graph, we generate noisy step signals, defined as $f_{i}=1+\sigma_{i}$ if $i$ belongs to the community, and $f_{i}=\sigma_{i}$ otherwise, where $\sigma_{i}\sim\mathcal{N}(0,0.3)$ is Gaussian i.i.d. noise. The goal is to classify each such signal according to the community it belongs to. The neural network architecture used for this task consisted of a spectral convolutional layer (based on Chebyshev or Cayley filters) with 32 output features, a mean pooling layer, and a softmax classifier for producing the final classification into one of the 15 classes. No regularization has been exploited in this setting. The classification accuracy is shown in Figure 3 (second to the top) along with examples of learned filters (bottom two). We observe that CayleyNet significantly outperforms ChebNet for smaller filter orders, with an improvement as large as 80%. Studying the filter responses, we note that due to the capability to learn the spectral zoom parameter, CayleyNet allows to generate band-pass filters in the low-frequency band that discriminate well the communities (Figure 3 bottom).

4.3 Complexity

We experimentally validated the computational complexity of our model applying filters of different order $r$ to synthetic 15-community graphs of different size $n$ using exact matrix inversion and approximation with different number of Jacobi iterations (Figure 6 in Appendix C). All times have been computed running 30 times the considered models and averaging the final results. As expected, approximate inversion guarantees $\mathcal{O}(n)$ complexity. We further conclude that typically very few Jacobi iterations are required (Figure 4 shows that our model with just one Jacobi iteration outperforms ChebNet for low-order filters on the community detection problem).

4.4 MNIST

Following [10, 30], for a toy example, we approached the classical MNIST digits classification as a learning problem on graphs. Each pixel of an image is a vertex of a graph (regular grid with 8-neighbor connectivity), and pixel color is a signal on the graph. We used a graph CNN architecture with two spectral convolutional layers based on Chebyshev and Cayley filters (producing 32 and 64 output features, respectively), interleaved with pooling layers performing 4-times graph coarsening using the Graclus algorithm [11], and finally a fully-connected layer (this architecture replicates the classical LeNet5, [25], architecture, which is shown for comparison). SGD+Momentum with learning rate equal to 0.02, momentum $m=0.9$, dropout probability $p=0.5$ and weight decay coefficient $\gamma=5\cdot 10^{-4}$ have been applied as described in [10]. MNIST classification results are reported in Table 1. CayleyNet (11 Jacobi iterations) achieves the same (near perfect) accuracy as ChebNet with filters of lower order ($r=12$ vs $25$). Examples of filters learned by ChebNet and CayleyNet are shown in Figure 2. 0.1776 +/- 0.06079 sec and 0.0268 +/- 0.00841 sec are respectively required by CayleyNet and ChebNet for analyzing a batch of 100 images at test time.

4.5 Citation network

Next, we address the problem of vertex classification on graphs using the popular CORA citation graph [36]. Each of the 2,708 vertices of the CORA graph represents a scientific paper, and an undirected unweighted edge represents a citation (5,429 edges in total). For each vertex, a 1,433-dimensional binary feature vector representing the content of the paper is given. The task is to classify each vertex into one of the 7 groundtruth classes, of labels. In the semi-supervised problem (transductive learning), the features of all vertices are known, but labels are given just for a subset of the nodes. The task is to learn a mapping, that takes the features at the nodes as inputs, and gives the labels as outputs. The mapping is trained by minimizing the label error at the nodes with known labels. After training, the the mapping is tested over the nodes in which the labels were unknown during training.

To present a deep comparison with recent state-of-the-art architectures, we analyze the performance of our model in two different settings: the classic semi-supervised problem presented in [23, 30, 39] with 140 training samples, 500 validation samples and 1,000 test samples and a relaxed version of this that exploits 1,708 vertices for training, 500 for validation and 500 for testing. We opted for a larger amount of training samples in our second experiment, in order to provide an estimate of the quality of CayleyNet in a situation which is less prone to overfitting. This provides a better overview of the goodness of the considered construction, since richer filters are less likely to produce lower performance, as opposed to the typical behavior when the available data is scarce. Cayley operators with matrix inversion have been considered in both settings for our solution. DCNN[1], GCN[23], MoNet[30] and GAT[39] have been used as terms of comparison.

On the standard split, we train CayleyNet by realizing filters as linear combinations of neighborhood descriptors obtained by applying Cayley filters on the signal of features. Two versions of CayleyNet have been implemented for this setting: a lightweight one exploiting two convolutional layers with 16 and 7 output features and a heavier one requiring 64 and 7 output features. This provides a valuable term of comparison with both the solutions presented in [23, 30] and [39] as same number of parameters is respectively required by our implementations. Normalized Laplacian has been used as reference. Adam with learning rate equal to $5\cdot 10^{-3}$, dropout probability $p=0.6$ and weight decay coefficient $\gamma=5\cdot 10^{-4}$ have been used for training. Table 3 presents the results we obtained with our solution, average performance over 50 runs are reported to guarantee accurate estimates. Performance of GCN, MoNet and GAT have been obtained from the respective papers, DCNN has been trained with 1 diffusion layer and 1 hop to guarantee same number of parameters. Our lighter version of CayleyNet outperforms DCNN, GCN and MoNet, while being defeated only by the recent GAT (which however exploits an attention mechanism for better discriminating relevant neighbors). Our heavier CayleyNet shows a significant drop in performance, likely because of overfitting on the small training set.

On our extended split, we analyze the behavior of CayleyNet and ChebNet for a variety of different polynomial orders. Two spectral convolutional layers with 16 and 7 outputs features have been used for implementing the two architectures. Adam with learning rate equal to $10^{-3}$, dropout probability $p=0.5$ and weight decay with coefficient $\gamma=5\cdot 10^{-4}$ have been used for training. Figure 5 presents the results of our analysis. Since ChebNet requires Laplacians with spectra bounded in $[-1,1]$, we consider both the normalized Laplacian (the two left figures), and the scaled unnormalized Laplacian $(2\boldsymbol{\Delta}/{\lambda_{max}}-\mathbf{I})$, where $\boldsymbol{\Delta}$ is the unnormalized Laplacian and $\lambda_{max}$ is its largest eigenvalue (the two right figures). For fair comparison, we fix the order of the filters (top figures) and the overall number of network parameters (bottom figures). In the bottom figures, the Cayley filters are restricted to even cosine polynomials by considering only real filter coefficients. The best CayleyNets consistently outperform the best ChebNets requiring at the same time less parameters (CayleyNet with order $r$ and complex coefficients requires a number of parameters equal to ChebNet with order $2r$). To further complete our analysis, we present the performance obtained by DCNN, GCN, MoNet and GAT on our extended split (Table 3). Two convolutional layers with order $r=1$, 1 head / 1 gaussian kernel, 16 and 7 outputs features have been used for GAT and MoNet222Filters have been realized as: $\mathbf{X}\mathbf{W_{0}}+\tilde{\mathbf{A}}\mathbf{X}\mathbf{W_{1}}$; with $\tilde{\mathbf{A}}$ the corresponding learned adjacency matrix. Attention has been computed for GAT only at second layer to ensure same number of parameters.; 3 convolutional layers with 32 and 16 hidden features have been used for GCN; 2 diffusion layers with 10 hidden features and 2 diffusion hops for DCNN. GCN, MoNet and GAT have been trained with mean cross-entropy, dropout probability $p=0.5$ and weight decay with coefficient $\gamma=5\cdot 10^{-4}$, DCNN has been trained with hinge loss and no regularization (as reported in [1]). CayleyNet appears as the third best approach for solving the considered semi-supervised classification task, and outperforms other spectral CNN methods.

4.6 Recommender system

In our final experiment, we applied CayleyNet to recommendation system, formulated as matrix completion problem on user and item graphs, [31]. The task is, given a sparsely sampled matrix of scores assigned by users (columns) to items (rows), to fill in the missing scores. The similarities between users and items are given in the form of column and row graphs, respectively. [31] approached this problem as learning with a Recurrent Graph CNN (RGCNN) architecture, using an extension of ChebNets to matrices defined on multiple graphs in order to extract spatial features from the score matrix; these features are then fed into an RNN producing a sequential estimation of the missing scores. Here, we repeated verbatim their experiment on the MovieLens dataset ([29]), replacing Chebyshev filters with Cayley filters. Following [31], to train our model we uniformly split the available training scores in two sets of equal dimension, 50% of the provided scores (data scores) are used to initialize the input matrix while the remaining 50% are used as training labels. SRGCNN is trained to reconstruct the missing labels from the few given data scores. At test time we initialize the input matrix only with the considered data scores to provide the network the same conditions it observed at training time. We used separable RGCNN architecture with two CayleyNets of order $r=4$ employing 15 Jacobi iterations. Adam with learning rate equal to $10^{-3}$ and regularization coefficient $\gamma=10^{-10}$ have been used for training333Values obtained from MGCNN implementation available at https://github.com/fmonti/mgcnn/. The results are reported in Table 4. To present a complete comparison we further extended the experiments reported in [31] by training sRGCNN with ChebNets of order 8, this provides an architecture with same number of parameters as the exploited CayleyNet (23K coefficients). Our version of sRGCNN outperforms all the competing methods, including the previous result with Chebyshev filters reported in [31]. sRGCNNs with Chebyshev polynomials of order 4 and 8 respectively require 0.0698 +/- 0.00275 sec and 0.0877 +/- 0.00362 sec at test time, sRGCNN with Cayley polynomials of order 4 and 15 jacobi iterations requires 0.165 +/- 0.00332 sec.

5 Conclusion

In this paper, we introduced a new efficient spectral graph CNN architecture that scales linearly with the dimension of the input data. Our architecture is based on a new class of complex rational Cayley filters that are localized in space, can represent any smooth spectral transfer function, and are highly regular. The key property of our model is its ability to specialize in narrow frequency bands with a small number of filter parameters, while still preserving locality in the spatial domain.

Experimental results on the MNIST, CORA and MovieLens datasets show the good performance of our construction wrt a variety of other approaches, and the superior performance with respect to other spectral CNN methods.

Appendix A Proof of Proposition 1

First note the following classical result for the approximation of $\mathbf{A}\mathbf{x}=\mathbf{b}$ using the Jacobi method: if the initial condition is $\mathbf{x}^{(0)}=\mathbf{0}$, then $(\mathbf{x}-\mathbf{x}^{(k)})=\mathbf{J}^{k}\mathbf{x}$. In our case, note that if we start with initial condition $\tilde{\mathbf{y}}^{(0)}_{j}=0$, the next iteration gives $\tilde{\mathbf{y}}^{(0)}_{j}=\mathbf{b}_{j}$, which is the initial condition from our construction. Therefore, since we are approximating $\mathbf{y}_{j}=\mathcal{C}(h\boldsymbol{\Delta})\tilde{\mathbf{y}}_{j-1}$ by $\tilde{\mathbf{y}}_{j}=\tilde{\mathbf{y}}_{j}^{(K)}$, we have

[TABLE]

Define the approximation error in $\mathcal{C}(h\boldsymbol{\Delta})^{j}\mathbf{f}$ by

[TABLE]

By the triangle inequality, by the fact that $\mathcal{C}^{j}(h\boldsymbol{\Delta})$ is unitary, and by (8)

[TABLE]

where the last inequality is due to

[TABLE]

Now, using standard norm bounds, in the general case we have $\left\|\mathbf{J}^{K+1}\right\|_{2}\leq\sqrt{n}\left\|\mathbf{J}^{K+1}\right\|_{\infty}$. Thus, by $\kappa=\left\|\mathbf{J}\right\|_{\infty}$ we have

[TABLE]

The solution of this recurrent sequence is

[TABLE]

If we use the version of the algorithm, in which each $\tilde{\mathbf{y}}_{j}$ is normalized, we get by (11) $e_{j}\leq e_{j-1}+\sqrt{n}\kappa^{K+1}$. The solution of this recurrent sequence is

[TABLE]

We denote in this case $M_{j}=j\sqrt{n}$

In case the graph is regular, we have $\mathbf{D}=d\mathbf{I}$. In the non-normalized Laplacian case,

[TABLE]

The spectral radius of $\boldsymbol{\Delta}$ is bounded by $2d$. This can be shown as follows. a value $\lambda$ is not an eigenvalue of $\boldsymbol{\Delta}$ (namely it is a regular value) if and only if $(\boldsymbol{\Delta}-\lambda\mathbf{I})$ is invertible. Moreover, the matrix $(\boldsymbol{\Delta}-\lambda\mathbf{I})$ is strictly dominant diagonal for any $\left|\lambda\right|>2d$. By Levy–Desplanques theorem ([17] Theorem 6.1.10), any strictly dominant diagonal matrix is invertible, which means that all of the eigenvalues of $\boldsymbol{\Delta}$ are less than $2d$ in their absolute value. As a result, the spectral radius of $(d\mathbf{I}-\boldsymbol{\Delta})$ is realized on the smallest eigenvalue of $\boldsymbol{\Delta}$, namely it is $\left|d-0\right|=d$. This means that the specral radius of $\mathbf{J}$ is $\frac{hd}{\sqrt{h^{2}d^{2}+1}}$. As a result $\left\|\mathbf{J}\right\|_{2}=\frac{hd}{\sqrt{h^{2}d^{2}+1}}=\kappa$. We can now continue from (11) to get

[TABLE]

As before, we get $e_{j}\leq j\kappa^{K+1}+\mathcal{O}(\kappa^{2K+2})$, and $e_{j}\leq j\kappa^{K+1}$ if each $\tilde{\mathbf{y}}_{j}$ is normalized. We denote in this case $M_{j}=j$.

In the case of the normalized Laplacian of a regular graph, the spectral radius of $\boldsymbol{\Delta}_{n}$ is bounded by $2$, and the diagonal entries are all 1. Equation (10) in this case reads

$\mathbf{J}=\frac{h}{h+i}(\mathbf{I}-\boldsymbol{\Delta}_{n})$, and $\mathbf{J}$ has spectral radius $\frac{h}{\sqrt{h^{2}+1}}$. Thus $\left\|\mathbf{J}\right\|_{2}=\frac{h}{\sqrt{h^{2}+1}}=\kappa$ and we continue as before to get

$e_{j}\leq j\kappa^{K+1}$ and $M_{j}=j$.

In all cases, we have by the triangle inequality

[TABLE]

Appendix B Proof of Theorem 4

In this proof we approximate $\mathbf{G}\boldsymbol{\delta}_{m}$ by $\widetilde{\mathbf{G}\boldsymbol{\delta}}_{m}$. Note that the signal $\boldsymbol{\delta}_{m}$ is supported on one vertex, and in the calculation of $\widetilde{\mathbf{G}\boldsymbol{\delta}}_{m}$, each Jacobi iteration increases the support of the signal by 1-hop. Therefore, the support of $\widetilde{\mathbf{G}\boldsymbol{\delta}}_{m}$ is the $r(K+1)$-hop neighborhood ${\cal N}_{r(K+1),m}$ of $m$. Denoting $l=r(K+1)$, and using Proposition 1, we get

[TABLE]

Appendix C Computational Complexity

In Figure 6 we compare the computational complexity of CayleyNey and ChebNet on our community detection problem.

Appendix D Back propagation

In this section we show how to differentiate Cayley filters $\mathbf{G}_{\mathbf{c},h}=g_{c,h}(\boldsymbol{\Delta})$ with respect to the complex coefficient vector $\mathbf{c}=(c_{0},\ldots,c_{r})$ and $h$. Since working with complex parameters is not standard, we explain in detail how this is done. One approach is to simply treat each complex parameter $c_{j}=c_{j}^{R}+ic_{j}^{I}$ as the pair of real parameters $(c_{j}^{R},c_{j}^{I})$, and to explicitly formulate Cayley polynomial with real numbers. This brute force formulation is suitable for automatic back propagation in software like TensorFlow. To justify the calculation from a theoretical standpoint, it is more convenient to consider a general calculus of variation approach to gradient descent.

Our goal is to minimize the loss function with respect to all of the coefficients of all filters in the network. Note that minimization is a set operation in its nature. Namely, a minimal value of a set doesn’t depend on additional structures endowed on the set, such as inner product, topology, Riemannian structure, and so on. This means that we are free to use the vector space structure and the inner product of our choice to define the gradient.

Consider a generic Cayley filter, applied on the real valued signal $\mathbf{f}$

[TABLE]

Let

[TABLE]

denote the dependency of the loss function on $(\mathbf{c},h)$. Our goal is to define an inner product, and calculate a gradient of $S(\mathbf{c},h)$ with respect to the complex coefficient vector $\mathbf{a}=(\mathbf{c},h)$. Given an inner product structure on the space of coefficients, a (variational) gradient at the point $\mathbf{a}$ of a scalar valued function $S$ is a vector $D[S;\mathbf{a}]$ such that for any other coefficient vector $\boldsymbol{\eta}$

[TABLE]

where $\epsilon$ denotes a scalar. If there is no such vector $D[S;\mathbf{a}]$, $S$ is not differentiable at $\mathbf{a}$ by definition. This definition can be extended to vector valued functions.

It can be shown that under the standard complex dot product of the coefficient space, $S$ is not differentiable in general. However, by defining a new inner product in the coefficient space, there is a way to make $S$ differentiable. For intuition, consider the vector space $\mathbb{C}$ with the classical inner product $\left\langle z,w\right\rangle=z\overline{w}$. Here, the space of differentiable functions $S:\mathbb{C}\rightarrow\mathbb{C}$ are the analytic functions. However, if we define the new inner product ${\rm Re}\left\langle z,w\right\rangle=z_{R}w_{R}+z_{I}w_{I}$, $\mathbb{C}$ is isometrically isomorphic to $\mathbb{R}^{2}$ with the standard real inner product. Moreover, if we treat the image space $\mathbb{C}$ of $S$ as $\mathbb{R}^{2}$, the space of differentiable functions $S:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ is the richer space of classical real differentiable functions. This procedure of defining a real vector space from a complex one is called realification ([9] page 117).

Given a general complex Hilbert space $\mathcal{H}$, it’s realification is the real Hilbert space $\mathcal{H}_{\mathbb{R}}$ defined to be $\mathcal{H}$ restricted to multiplication by real scalars. The inner product of two vectors $\mathbf{f},\mathbf{g}$ in $\mathcal{H}_{\mathbb{R}}$ is defined to be the real part of the inner product in $\mathcal{H}$, ${\rm Re}\left\langle\mathbf{f},\mathbf{g}\right\rangle$.

In our case, we treat the coefficient space as $\mathbb{R}^{2r}$. Indeed, $h$ and $c_{0}$ are real, and the rest of the coefficients are complex. Let us differentiate $\mathbf{G}_{\mathbf{c},h}\mathbf{f}$ with respect to $\mathbf{c}$ by definition.

[TABLE]

[TABLE]

Thus, by the definition of the inner product in the (realificated) coefficient space, the differential of $\mathbf{G}_{\mathbf{c},h}\mathbf{f}$ with respect to $\mathbf{c}$ is given by the vector $D_{\mathbf{c}}[\mathbf{G}_{\mathbf{c},h}\mathbf{f};(h,\mathbf{c})]$ with vector valued entries

[TABLE]

Next we show that back propagation works in the realificated space $\mathbb{R}^{2r}$ in the usual way. By (12) and by the chain rule, we can write up to $o(\epsilon)$

[TABLE]

where $\nabla$ denotes the gradient with respect to the signal, and $\cdot$ is the usual real dot product between real valued signals. Let us write in short $\nabla F=\nabla F(\mathbf{G}_{\mathbf{c},h}\mathbf{f})$, and obtain

[TABLE]

This shows, by the definition of the inner product in the (realificated) coefficient space, that

[TABLE]

Namely, the chain rule works in the usual way in the realificated space $\mathbb{R}^{2r}$.

Next we calculate the partial derivative with respect to the spectral zoom $h$. We start by standard Calculus on the function $g_{\mathbf{c},h}(\lambda)$, $\lambda\in\mathbb{R}$, and get

[TABLE]

where $\mathcal{C}^{\prime}(x)=(x+i)^{-1}(1-\mathcal{C}(x))$. We can interpret this as a calculation on the spectrum of $\mathbf{G}_{\mathbf{c},h}=g_{\mathbf{c},h}(\boldsymbol{\Delta})$. Then, by the fact that $\boldsymbol{\Delta}$ is a bounded normal operator, we can carry the calculation to $\mathbf{G}_{\mathbf{c},h}=g_{\mathbf{c},h}(\boldsymbol{\Delta})$ via functional calculus. We thus obtain

[TABLE]

where $\mathcal{C}^{\prime}(h\boldsymbol{\Delta})=(h\boldsymbol{\Delta}+i\mathbf{I})^{-1}(\mathbf{I}-\mathcal{C}(h\boldsymbol{\Delta}))$.

Acknowledgment

MB and FM are partially supported by ERC Consolidator Grant No. 724228 (LEMAN), Google Faculty Research Awards, Amazon AWS ML Research Award, Royal Society Wolfson Research Merit Award, and Rudolf Diesel fellowship at the Institute for Advanced Studies, TU Munich.