# COS-GNN: Continuous Spiking Graph Neural Networks

Nan Yin  Mengzhu Wang  Li Shen  Hitesh Laxmichand Patel  Baopu Li  Bin Gu
Huan Xiong

###### Abstract

Continuous graph neural networks (CGNNs) have garnered significant attention
due to their ability to generalize existing discrete graph neural networks
(GNNs) by introducing continuous dynamics. They typically draw inspiration
from diffusion-based methods to introduce a novel propagation scheme, which is
analyzed using ordinary differential equations (ODE). However, the
implementation of CGNNs requires significant computational power, making them
challenging to deploy on battery-powered devices. Inspired by recent spiking
neural networks (SNNs), which emulate a biological inference process and
provide an energy-efficient neural architecture, we incorporate the SNNs with
CGNNs in a unified framework, named Continuous Spiking Graph Neural Networks
(COS-GNN). We employ SNNs for graph node representation at each time step,
which are further integrated into the ODE process along with time. To enhance
information preservation and mitigate information loss in SNNs, we introduce
the high-order structure of COS-GNN, which utilizes the second-order ODE for
spiking representation and continuous propagation. Moreover, we provide the
theoretical proof that COS-GNN effectively mitigates the issues of exploding
and vanishing gradients, enabling us to capture long-range dependencies
between nodes. Experimental results on graph-based learning tasks demonstrate
the effectiveness of the proposed COS-GNN over competitive baselines.

Machine Learning, ICML

##  1 Introduction

Continuous graph neural networks (CGNNs) play a vital role in shaping modern
societies and have been extensively studied across multiple fields in science
and engineering, such as social network prediction (Zhang et al., 2022b;
Hafiene et al., 2020; Liao et al., 2021), COVID-19 forecasting (Luo et al.,
2023), and interacting dynamic system learning (Huang et al., 2021).
Typically, CGNNs extend neural ordinary differential equation (ODE) methods to
graph neural networks (GNNs) to model the continuous change of node
representations, efficiently addressing dynamic system learning problems
(Rusch et al., 2022; Xhonneux et al., 2020). This makes modeling and
understanding the intricate dynamics of complex relational systems a major
area of research with various applications (Hsieh et al., 2021; Yin et al.,
2024b, 2022b, 2023a; Pang et al., 2023).

Recently proposed CGNNs (Battaglia et al., 2016; Kipf et al., 2018) approaches
typically utilize GNNs to obtain node representations at each timestamp, which
are subsequently employed for trend prediction. However, these methods are not
able to solve the irregularly sampled issue (Huang et al., 2020, 2021). In
contrast, recent ODE works (Poli et al., 2019; Gupta et al., 2022) have proven
to be effective in modeling system dynamics when dealing with missing data
(Chen et al., 2018). They extend ODE to model interacting dynamical systems,
which essentially integrate GNNs with ODE to capture spatio-temporal
relationships within dynamical systems. However, for modeling long-term node
dependencies, CGNNs tend to consume a substantial amount of energy, posing
challenges for their deployment on battery-powered devices.

Inspired by recent work of spiking neural networks (SNNs) (Brette et al.,
2007), which transform continuous input into discrete spikes, providing a more
intuitive and streamlined approach to inference and model training compared
with conventional networks (Zhang et al., 2022a; Maass, 1997), we aim to
integrate the energy-efficient attributes of SNNs into CGNNs to confront
foundational challenges. However, the combination of SNNs and CGNNs is
difficult due to the following challenges: (1) How to incorporate the SNNs and
CGNNs into a unified framework? SNNs and CGNNs operate by integrating
information along two distinct time dimensions: the latency dimension in SNNs
and the time dimension in graph-based models. An important challenge is how to
amalgamate these dual processes into a unified framework, preserving the
energy-efficient characteristics of SNNs while also harnessing the dynamic
learning capabilities of CGNNs. (2) How to alleviate the problem of
information loss of SNNs? SNNs fulfill the need for low-power consumption by
discretizing continuous features. However, this binarization process often
results in a substantial loss of fine-grained information, leading to
diminished performance. Consequently, mitigating the issue of information loss
in SNNs represents another significant challenge. (3) How to guarantee the
stability of the proposed framework? Traditional graph-based techniques
encounter the challenge of exploding and vanishing gradient problems when
modeling long-rang dependencies of GNNs. Hence, devising a stable model for
graph learning with theoretical guarantee constitutes the third prominent
challenge.

To tackle these challenges, we propose a framework named Continuous Spiking
Graph Neural Networks (COS-GNN) for graph-based learning tasks. Specifically,
to address the first challenge, we approach it by considering SNNs as a type
of ODE and integrating it with CGNNs into the partial differential equations
(PDE) framework, denoted as COS-GNN-1st. The PDE structure models the initial
representation using SNNs by considering the inner time latency at each time
step and evaluates the dynamics over time with CGNNs. Furthermore, to address
the information loss problem caused by SNNs, we derive a high-order spike
representation and introduce the second-order PDE structure, referred to as
COS-GNN-2nd. Moreover, we provide empirical evidence that COS-GNN successfully
mitigates issues related to exploding and vanishing gradients, enabling the
capture of long-range dependencies among nodes. We perform comprehensive
evaluations of COS-GNN against state-of-the-art methods on various graph-based
datasets, showcasing the efficacy and versatility of proposed approach.

The contributions can be summarized as follows:

  * •

Novel Architecture. We incorporate SNNs and CGNNs into a PDE framework, which
reserves energy-efficient properties of SNNs while preserving the ability to
capture continuous changes in CGNNs. Additionally, we introduce the high-order
PDE structure (COS-GNN-2nd), which addresses the information loss challenge
encountered in COS-GNN-1st.

  * •

High-order Spike Representation. We are the first to derive the second-order
spike representation and study the backpropagation of second-order spike ODE
to mitigate the information loss problem. This innovation is subsequently
applied in the COS-GNN-2nd.

  * •

Theoretical Analysis. We provide a theoretical proof demonstrating that COS-
GNN effectively mitigates the issue of exploding and vanishing gradients,
ensuring the stability of our proposed method.

  * •

Extensive Experiments. We evaluate the proposed COS-GNN on extensive graph-
based learning datasets, which evaluate that our proposed COS-GNN outperforms
the variety of state-of-the-art methods.

##  2 Related Work

Continuous Graph Neural Networks. Recently, numerous methods based on CGNNs
have emerged for modeling dynamic interaction systems (Battaglia et al., 2016;
Kipf et al., 2018; Chen et al., 2018; Ju et al., 2024). These methods commonly
employ GNNs to initialize node representations at discrete timestamps, which
are then utilized for predicting node behaviors. Nevertheless, these discrete
methodologies often necessitate the presence of all nodes at each timestamp
(Huang et al., 2020, 2021; Yin et al., 2023b, 2022a), which is challenging to
achieve in real-world scenarios. In contrast, ODE has proven to be effective
in modeling system dynamics when dealing with missing data (Chen et al.,
2018). Recent works (Poli et al., 2019; Gupta et al., 2022) involve
initializing state representations with GNNs, followed by the establishment of
a neural ODE model for both nodes and edges, guiding the evolution of the
dynamical system. However, CGNNs exhibit substantial energy consumption when
modeling long-term node dependencies. In our work, we integrate energy-
efficient SNNs into CGNNs, thereby retaining the low energy characteristics of
SNNs while harnessing the dynamic learning capabilities of CGNNs.

Spiking Neural Networks. SNNs (Brette et al., 2007; Maass, 1997) have emerged
as a promising solution for addressing graph machine learning challenges.
Currently, there are two main directions in the field of SNNs. The first
direction focuses on establishing a connection between SNNs and Artificial
Neural Networks (ANNs), allowing for the conversion from ANNs to SNNs
(Rueckauer et al., 2017; Rathi et al., 2020). However, these methods often
necessitate a relatively extensive time latency to achieve comparable
performance to ANNs, resulting in higher latency and typically increased
energy consumption. The second involves the direct training of SNNs using
backpropagation techniques (Bohte et al., 2000; Esser et al., 2015; Bellec et
al., 2018). These approaches follow the backpropagation through time (BPTT)
framework, offering reduced latency but demanding substantial training memory.
Nevertheless, there is a scarcity of research concentrating on dynamic spiking
graphs, and the solitary existing work (Li et al., 2023; Yin et al., 2024a)
merely integrates SNNs into dynamic graphs, falling short in effectively
capturing the continuous changes and subtle dynamics inherent in time series
data. Our approach ingeniously combines SNNs with CGNNs to effectively capture
the continuous changes and subtle dynamics within graph time series data.

##  3 Preliminaries

###  3.1 Continuous Graph Neural Networks

Given a graph $G=(\mathcal{V},\mathcal{E})$ with the node set $\mathcal{V}$
and the edge set $\mathcal{E}$. $\mathbf{X}\in\mathbb{R}^{|\mathcal{V}|\times
d}$ is the node feature matrix, $d$ is the feature dimension. The binary
adjacency matrix denoted as
$\mathbf{A}\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}$, where $a_{ij}=1$
denotes there exists a connection between node $i$ and $j$, and vice versa.
Our goal is to learn a node representation $\mathbf{H}$ for downstairs tasks.

First-order CGNNs. The first graph ODE-based CGNN method is proposed by
(Xhonneux et al., 2020). Considering the Simple GNN (Wu et al., 2019) with
$\mathbf{H}_{n+1}=\mathbf{AH}_{n}+\mathbf{H}_{0}$, the solution of the ODE is
given by:

$\displaystyle\frac{d\mathbf{H}(t)}{dt}=ln\mathbf{AH}(t)+\mathbf{E},$ (1)
$\displaystyle\mathbf{H}(t)=(\mathbf{A}$
$\displaystyle-\mathbf{I})^{-1}(e^{(\mathbf{A}-\mathbf{I})t}-\mathbf{I})\mathbf{E}+e^{(\mathbf{A}-\mathbf{I})t}\mathbf{E},$

where $\mathbf{E}=\varepsilon\mathbf{(}X)$ is the output of the encoder
$\varepsilon$ and the initial value
$\mathbf{H}(0)=(ln\mathbf{A})^{-1}(\mathbf{A}-\mathbf{I})\mathbf{E}$.

Second-order CGNNs. To model high-order correlations in long-term temporal
trends, (Rusch et al., 2022) first propose the second-order graph ODE, which
is represented as:

$\mathbf{X}^{{}^{\prime\prime}}=\sigma(\mathbf{F}_{\theta}(\mathbf{X},t))-\gamma\mathbf{X}-\alpha\mathbf{X}^{{}^{\prime}},$
(2)

where
$\left(\mathbf{F}_{\theta}(\mathbf{X},t)\right)_{i}=\mathbf{F}_{\theta}\left(\mathbf{X}_{i}(t),\mathbf{X}_{j}(t),t\right)$
is a learnable coupling function with parameters $\theta$. Due to the
unavailability of an analytical solution for Eqn. 2, GraphCON (Rusch et al.,
2022) addresses it through an iterative numerical solver employing a suitable
time discretization method. GraphCON utilizes the IMEX (implicit-explicit)
time-stepping scheme, an extension of the symplectic Euler method (Hairer et
al., 1993) that accommodates systems with an additional damping term.

$\displaystyle\mathbf{Y}^{n}=$ $\displaystyle\mathbf{Y}^{n-1}+\Delta
t\left[\sigma(\mathbf{F}_{\theta}(\mathbf{X}^{n-1},t^{n-1}))\right.$ (3)
$\displaystyle\left.-\gamma\mathbf{X}^{n-1}-\alpha\mathbf{Y}^{n-1}\right],$
$\displaystyle\mathbf{X}^{n}=$ $\displaystyle\mathbf{X}^{n-1}+\Delta
t\mathbf{Y}^{n},\,n=1,\cdots,N,$

where $\Delta t>0$ is a fixed time-step and $\mathbf{Y}^{n}$, $\mathbf{X}^{n}$
denote the hidden node features at time $t^{n}=n\Delta t$.

###  3.2 Spiking Neural Networks

First-order SNNs. In contrast to traditional artificial neural networks, SNNs
convert input data into binary spikes over time, with each neuron in the SNNs
maintaining a membrane potential that accumulates input spikes. A spike is
produced as an output when the membrane potential exceeds a threshold. And the
first-order SNNs is formulated as:

$\displaystyle u_{\tau+1,i}=\lambda(u_{\tau,i}-$ $\displaystyle
V_{th}s_{\tau,i})+\sum_{j}w_{ij}s_{\tau,j}+b,$ (4) $\displaystyle
s_{\tau+1,i}=$ $\displaystyle\mathbb{H}(u_{\tau+1,i}-V_{th}),$

where $\mathbb{H}(x)$ is the Heaviside step function, which is the non-
differentiable spiking function. $s_{\tau,i}$ is the binary spike train of
neuron $i$, and $\lambda$ is the constant. $w_{ij}$ and $b$ are the weights
and bias of each neuron.

Second-order SNNs. The first-order neuron models assume that an input voltage
spike causes an immediate change in synaptic current, affecting the membrane
potential. However, in reality, a spike leads to the gradual release of
neurotransmitters from the pre-synaptic neuron to the post-synaptic neuron. To
capture the temporal dynamics, we utilize the synaptic conductance-based LIF
model, which considers the gradual changes in input current over time. To
solve this, (Eshraghian et al., 2023) propose the second-order SNN, which is
formulated as:

$\left\\{\begin{array}[]{lr}I_{\tau+1}=\alpha I_{\tau}+WX_{\tau+1},&\\\
u_{\tau+1,i}=\beta u_{\tau,i}+I_{\tau+1,i}-R,&\\\
s_{\tau,i}=\mathbb{H}(u_{\tau+1,i}-V_{th}),&\end{array}\right.$ (5)

where $\alpha=exp(-\Delta t/\tau_{syn})$ and $\beta=exp(-\Delta
t/\tau_{mem})$, $\tau_{syn}$ models the time constant of the synaptic current
in an analogous way to how $\tau_{mem}$ models the time constant of the
membrane potential.

##  4 Methodology

In this part, we present the proposed COS-GNN for modeling the continuous
spiking graph neural networks. COS-GNN combines SNNs with CGNNs in a framework
of PDE, which preserves the advantage of SNNs and CGNNs simultaneously. To
mitigate the problem of information loss attributed to SNNs, we involve the
derivation of second-order spike representation and differentiation for
second-order SNNs, and then introduce the high-order PDE, referred to as COS-
GNN-2nd. Finally, we present a theoretical proof to ensure that our proposed
COS-GNN effectively mitigates the challenges associated with gradient
exploding and vanishing. The framework of COS-GNN is shown in Figure 1.

###  4.1 First-order COS-GNN

Specifically, SNNs propagate information within the time latency $\tau$, and
the CGNNs evaluate the feature evolution along with time $t$. We propose the
first-order COS-GNN to incorporate SNNs with CGNNs in the form of PDE with
their respective timelines. Intuitively, information is interactively
propagated in dual time dimensions by both the SNNs and the CGNNs, which is
formulated as:

######  Proposition 4.1.

Define the first-order SNNs as
$\frac{du_{t}^{\tau}}{d\tau}=g(u_{t}^{\tau},\tau)$, and first-order CGNNs as
$\frac{du_{t}^{\tau}}{dt}=f(u^{\tau}_{t},t)$, then the first-order COS-GNN can
be formulated as:

$\displaystyle u_{T}^{N}=$ $\displaystyle
2\int_{0}^{T-1}f\left(\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy$ (6)
$\displaystyle+\int_{T-1}^{T}f\left(\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy,$

or:

$\displaystyle u_{T}^{N}=$ $\displaystyle
2\int_{0}^{N}g\left(\int_{0}^{T-1}f(u_{y}^{x},x)dy\right)dx$ (7)
$\displaystyle+\int_{T-1}^{T}f\left(\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy.$

where $N$ is the total latency of SNNs, and $T$ is the time step length of
CGNNs, $u_{y}^{x}$ denotes the neuron membrane on latency $x\in[0,N]$ and time
step $y\in[0,T]$. Details of derivation are shown in Appendix A.

Figure 1: Overview of the proposed COS-GNN. The proposed COS-GNN takes a graph
with node features as input, which are initially encoded using the SNNs
(first-order or second-order). Subsequently, a first-order or second-order
continuous graph process is employed to evolve the representation over
different time step. Finally, the representation is projected for downstream
tasks.

To provide a comprehensive understanding of Proposition 4.1, we approach it
from two distinct perspectives. First, in Eqn. 6, within the temporal
dimension where $0\leq t\leq T$, the SNNs initially compute spike
representations at each time step $t$. These representations are then employed
in the evolution of the CGNNs process. Second, as outlined in Eqn. 7, we focus
to the latency dimension $\tau$ with $0\leq\tau\leq N$. Then, the CGNNs model
the evolution of node embeddings along a time series $t$ and subsequently
integrate them across the latency dimension.

In our implementation of the first-order COS-GNN, we employ Eqn. 6 by
initially calculating spike representations at each time step, followed by
modeling the evolution of node embeddings. As described in (Meng et al.,
2022), the first-order spike representation at time step $t$ is denoted as:
$\mathbf{H}(t)=\frac{\sum_{\tau=1}^{N}\lambda^{N-\tau}s_{\tau}^{t}}{\sum_{\tau=1}^{N}\lambda^{N-\tau}}$.
By combining Eqn. 1, we have:

$\frac{d\mathbf{H}(t)}{dt}=ln\mathbf{AH}(t)+\frac{\sum_{\tau=1}^{N}\lambda^{N-\tau}s_{\tau}^{0}}{\sum_{\tau=1}^{N}\lambda^{N-\tau}},$
(8)

where $s_{\tau}^{t}$ is the binary spike on latency $\tau$ and time step $t$,
and $\lambda=exp(-\frac{\Delta t}{\kappa})$ with $\Delta t\ll\kappa$, $\kappa$
is the time constant.

###  4.2 Second-order Spiking Neural Networks

The proposed first-order COS-GNN effectively addresses the challenge of
combining SNNs with CGNNs to achieve energy-efficient modeling for continuous
graph learning. However, the first-order SNNs typically suffer from the
information loss issue (Yin et al., 2023c). Motivated by recent advancements
in high-order structure modeling (Luo et al., 2023), which addresses high-
order correlations of nodes, we introduce the second-order COS-GNN to tackle
the aforementioned issue. In this part, we begin by deriving the second-order
spike representation and investigating the backpropagation of second-order
SNNs.

Second-order SNNs Forward Propagation. We set the forward propagation layer of
SNNs to $L$. According to Eq. 5, the spiking dynamics for SNNs can be
formulated as:

$\displaystyle u^{i}(n)=$
$\displaystyle\beta^{i}u^{i}(n-1)+(1-\beta^{i})\frac{V_{th}^{i-1}}{\Delta
t}\left(\alpha^{i}I^{i-1}(n-1)\right.$
$\displaystyle\left.+\mathbf{W}^{i}s^{i-1}(n)\right)-V^{i}_{th}s^{i}(n),$

where $i=1,\cdots,L$ denotes the $i$-th layer, $s^{0}$ and $s^{i}$ denote the
input and output of SNNs. $I^{i}$ is the input of the $i$-th layer,
$n=1,\cdots,N$ is the time step and $N$ is the latency.
$\alpha^{i}=exp(-\Delta\tau/\tau_{syn}^{i})$,
$\beta^{i}=exp(-\Delta\tau/\tau_{mem}^{i})$ and
$0<\Delta\tau\ll\\{\tau_{syn}^{i},\tau_{mem}^{i}\\}$.

Second-order Spike Representation. Considering the second-order SNNs model
defined by Eqn. 5, we first define the weighted average input current as
$\hat{I}(N)=\frac{1}{(\beta-\alpha)^{2}}\frac{\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})I_{in}(n)}{\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})}$,
and the scaled weighted firing rate as
$\hat{a}(N)=\frac{1}{\beta^{2}}\frac{V_{th}\sum_{n=0}^{N-1}\beta^{N-n}s(n)}{\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})\Delta\tau}$.
Then, we use $\hat{a}(N)$ denote the spike representation of
$\\{s(n)\\}_{n=1}^{N}$. Similarly to the first-order spike representation
(Meng et al., 2022), we directly determine the relationship between
$\hat{I}(N)$ and $\hat{a}(N)$ using a differentiable mapping. Specifically, by
combing Eqn. 5, we have:

$\displaystyle u(\tau+1)=$ $\displaystyle\beta u(\tau)+\alpha
I_{syn}(\tau)+I_{input}(\tau)-V_{th}s(\tau)$ (9) $\displaystyle=$
$\displaystyle\beta^{2}u(\tau-k+1)+\alpha\sum_{i=0}^{k-1}\beta^{i}I_{syn}(\tau-i)$
$\displaystyle+\sum_{i=0}^{k-1}\beta^{i}(I_{input}(\tau-i)-V_{th}s(\tau-i)).$

By summing Eq. 9 over $n=1$ to $N$, we have:

$\displaystyle u(N)=$
$\displaystyle\frac{1}{\beta-\alpha}\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})I_{in}(n)$
(10) $\displaystyle-\frac{1}{\beta}\sum_{n=0}^{N-1}\beta^{N-n}V_{th}s(n).$

Dividing Eq. 10 by $\Delta\tau\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})$:

$\displaystyle\hat{a}(N)=$
$\displaystyle\frac{\beta-\alpha}{\beta}\frac{\hat{I}(N)}{\Delta\tau}-\frac{u(N)}{\Delta\tau\beta\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})}$
$\displaystyle\approx$
$\displaystyle\frac{\tau_{syn}\tau_{mem}}{\tau_{mem}-\tau_{syn}}\hat{I}(N)$
$\displaystyle-\frac{u(N)}{\Delta\tau\beta\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})},$

since $\lim\limits_{\Delta\tau\to
0}\frac{1-\alpha/\beta}{\Delta\tau}=\frac{\tau_{syn}\tau_{mem}}{\tau_{mem}-\tau_{syn}}$
and $\Delta\tau\ll\frac{1}{\tau_{syn}}-\frac{1}{\tau_{mem}}$, we can
approximate $\frac{\beta-\alpha}{\beta\Delta\tau}$ by
$\frac{\tau_{syn}\tau_{mem}}{\tau_{mem}-\tau_{syn}}$. Following (Meng et al.,
2022), and take $\hat{a}(N)\in[0,\frac{V_{th}}{\Delta\tau}]$ into
consideration and assume $V_{th}$ is small, we ignore the term
$\frac{u(N)}{\Delta\tau\beta\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})}$, and
approximate $\hat{a}(N)$ with:

$\lim\limits_{N\to\infty}\hat{a}(N)\approx
clamp\left(\frac{\tau_{syn}\tau_{mem}}{\tau_{mem}-\tau_{syn}}\hat{I}(N),0,\frac{V_{th}}{\Delta\tau}\right),$
(11)

where $clamp(x,a,b)=max(a,min(x,b))$. During the training of the second-order
SNNs, we have Proposition 4.2, which is similar to (Meng et al., 2022), and
the detailed derivation is shown in the Appendix B.

######  Proposition 4.2.

Define
$\hat{a}^{0}(N)=\frac{\sum_{n=0}^{N-1}\beta^{N-n-2}_{i}s^{0}(n)}{\sum_{n=0}^{N-1}(\beta^{N-n}_{i}-\alpha^{N-n}_{i})\Delta\tau}$
and
$\hat{a}^{i}(N)=\frac{V_{th}^{i}\sum_{n=0}^{N-1}\beta^{N-n-2}_{i}s^{i}(n)}{\sum_{n=0}^{N-1}(\beta^{N-n}_{i}-\alpha^{N-n}_{i})\Delta\tau}$,
$\forall i=1,\cdots,L$, where $\alpha^{i}=exp(-\Delta\tau/\tau_{syn}^{i})$ and
$\beta^{i}=exp(-\Delta\tau/\tau_{mem}^{i})$. Further, define the
differentiable mappings

$\mathbf{z}^{i}\\!=\\!clamp\\!\left(\\!\frac{\tau_{syn}^{i}\tau_{mem}^{i}}{\tau_{mem}^{i}\\!-\\!\tau_{syn}^{i}}\mathbf{W}^{i}\mathbf{z}^{i-1},0,\frac{V_{th}^{i}}{\Delta\tau}\\!\right)\\!,i=1,\cdots,L.$

If $\lim\limits_{N\to\infty}\hat{a}^{i}(N)=\mathbf{z}^{i}$ for
$i=0,1,\cdots,L-1$, then $\hat{a}^{i+1}(N)\approx\mathbf{z}^{i+1}$ when
$N\to\infty$.

Differentiation on Second-order Spike Representation. In this part, we use the
spike representation to drive the backpropagation training algorithm for
second-order SNNs. With the forward propagation of the $i$-th layers, we get
the output of SNN with $s^{i}=\\{s^{i}(1),\cdots,s^{i}(N)\\}$. We define the
spike representation operator
$r(s)=\frac{1}{\beta^{2}}\frac{V_{th}\sum_{n=0}^{N-1}\beta^{N-n}s(n)}{\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})\Delta\tau}$,
and get the final output $\mathbf{o}^{L}=r(s^{L})$. For the simple second-
order SNN, assuming the loss function as $\mathcal{L}$, we calculate the
gradient $\frac{\partial\mathcal{L}}{\partial\mathbf{W}^{i}}$ as:

$\displaystyle\frac{\partial\mathcal{L}}{\partial\mathbf{W}^{i}}=\frac{\partial\mathcal{L}}{\partial\mathbf{o}^{i}}\frac{\partial\mathbf{o}^{i}}{\partial\mathbf{W}^{i}}=\frac{\partial\mathcal{L}}{\partial\mathbf{o}^{i+1}}\frac{\partial\mathbf{o}^{i+1}}{\partial\mathbf{o}^{i}}\frac{\partial\mathbf{o}^{i}}{\partial\mathbf{W}^{i}},$
(12) $\displaystyle\mathbf{o}^{i}=r(s^{i})\approx
clamp\left(\mathbf{W}^{i}r(s^{i-1}),0,\frac{V_{th}^{i}}{\Delta\tau}\right).$

Therefore, we can compute the gradient of second-order SNNs by calculating
$\frac{\partial\mathbf{o}^{i+1}}{\partial\mathbf{o}^{i}}$ and
$\frac{\partial\mathbf{o}^{i}}{\partial\mathbf{W}^{i}}$ based on Eqn. 12.

###  4.3 Second-order COS-GNN

Having obtained the second-order spike representation for SNNs, we introduce
the second-order COS-GNN. While obtaining an analytical solution for the
second-order COS-GNN may not be feasible, we can derive a conclusion similar
to Proposition 4.1. The specifics are presented as follows.

######  Proposition 4.3.

Define the second-order SNNs as
$\frac{d^{2}u_{t}^{\tau}}{d\tau^{2}}+\delta\frac{du_{t}^{\tau}}{d\tau}=e(u^{\tau}_{t},\tau)$,
and second-order CGNNs as
$\frac{d^{2}u_{t}^{\tau}}{dt^{2}}+\gamma\frac{du_{t}^{\tau}}{dt}=h(u^{\tau}_{t},t)$,
then the second-order COS-GNN follows:

$\displaystyle
u_{t}^{\tau}=\\!\\!\int_{0}^{T}\\!\\!h\left(\int_{0}^{N}e(u_{t}^{\tau})d\tau\right)dt=\\!\\!\int_{0}^{N}\\!\\!e\left(\int_{0}^{T}h(u_{t}^{\tau})dt\right)d\tau,$
$\frac{\partial^{2}u_{t}^{\tau}}{\partial\tau^{2}}+\delta\frac{\partial
u_{t}^{\tau}}{\partial\tau}=g(u_{t}^{\tau}),\quad\frac{\partial^{2}u_{t}^{\tau}}{\partial
t^{2}}+\gamma\frac{\partial u_{t}^{\tau}}{\partial t}=f(u_{t}^{\tau}),$

where
$e(u_{t}^{\tau})=\int_{0}^{N}g(u_{t}^{\tau})d\tau-\delta(u_{t}^{N}-u_{t}^{0})$,
$h(u_{t}^{\tau})=\int_{0}^{T}f(u_{t}^{\tau})dt-\gamma(u_{T}^{\tau}-u_{0}^{\tau})$,
$\frac{\partial e(u_{t}^{\tau})}{\partial\tau}=g(u_{t}^{\tau})$ and
$\frac{\partial h(u_{t}^{\tau})}{\partial t}=f(u_{t}^{\tau})$. $\delta$ and
$\gamma$ are the hyperparameters of second-order SNNs and CGNNs.

The details are derived in Appendix C. Similarly to the first-order COS-GNN,
we implement the second-order COS-GNN by calculating the spike representation
on each time step and then model the node embedding with Eqn. 3. Furthermore,
to optimize the second-order COS-GNN, we analyze the differentiation of the
second-order COS-GNN. Denote the loss function as
$\mathcal{L}=\sum\limits_{i\in\mathcal{V}}\left|\mathbf{X}_{i}^{N}-\mathbf{\bar{X}}_{i}\right|^{2}$,
and $\mathbf{\bar{X}}_{i}$ is the label of node $i$. With the chain rule, we
have:
$\frac{\partial\mathcal{L}}{\partial\mathbf{W}^{l}}=\frac{\partial\mathcal{L}}{\partial\mathbf{o}^{N}_{T}}\frac{\partial\mathbf{o}^{N}_{T}}{\partial\mathbf{o}^{l}_{T}}\frac{\partial\mathbf{o}^{l}_{T}}{\partial\mathbf{W}^{l}}$.
From (Rusch et al., 2022), we have the conclusion that the traditional GNN
model has the problem of gradient exponentially or vanishing, thus, we study
the upper bound of the proposed COS-GNN.

######  Proposition 4.4.

Let $\mathbf{X}^{n}$ and $\mathbf{Y}^{n}$ be the node features, generated by
Eqn. 3, and $\Delta t\ll 1$. The gradient of the second-order CGNNs
$\mathbf{W}_{l}$ is bounded as Eqn. 13, and the gradient of the second-order
SNNs $\mathbf{W}^{k}$ is bounded as Eqn. 14:

$\displaystyle\left|\frac{\partial\mathcal{L}}{\partial\mathbf{W}_{l}}\right|$
$\displaystyle\leq\frac{\beta^{{}^{\prime}}\hat{\mathbf{D}}\Delta t(1+\Gamma
T\Delta t)}{v}\left(\max\limits_{1\leq i\leq
v}(|\mathbf{X}_{i}^{0}|+|\mathbf{Y}_{i}^{0}|)\right)$ (13)
$\displaystyle+\frac{\beta^{{}^{\prime}}\hat{\mathbf{D}}\Delta t(1+\Gamma
T\Delta t)}{v}\left(\max\limits_{1\leq i\leq
v}|\bar{\mathbf{X}_{i}}|+\beta\sqrt{T\Delta t}\right)^{2},$
$\displaystyle\left|\frac{\partial\mathcal{L}}{\partial\mathbf{W}^{k}}\right|\leq$
$\displaystyle\frac{(1+T\Gamma\Delta
t)(1+N\Theta\Delta\tau)V_{th}}{v\beta^{2}\Delta\tau}\left(\max\limits_{1\leq
i\leq v}|\mathbf{X}_{i}^{N}|\right.$ (14)
$\displaystyle\left.+\max\limits_{1\leq i\leq
v}|\bar{\mathbf{X}}_{i}|\right).$

where $\beta=\max\limits_{x}|\sigma(x)|$,
$\beta^{{}^{\prime}}=\max\limits_{x}|\sigma^{{}^{\prime}}(x)|$,
$\hat{D}=\max\limits_{i,j\in\mathcal{V}}\frac{1}{\sqrt{d_{i}d_{j}}}$, and
$\Gamma:=6+4\beta^{{}^{\prime}}\hat{D}\max\limits_{1\leq n\leq
T}||\mathbf{W}^{n}||_{1}$,
$\Theta:=6+4\beta^{{}^{\prime}}\hat{D}\max\limits_{1\leq n\leq
N}||\mathbf{W}^{n}||_{1}.$ $d_{i}$ is the degree of node $i$,
$\mathbf{\bar{X}}_{i}$ is the label of node $i$. Eqn. 13 can be obtained from
(Rusch et al., 2022) directly, and the derivation of the Eqn. 14 is presented
in Appendix D.

The upper bound in Proposition 4.4 demonstrates that the total gradient
remains globally bounded, regardless of the number of CGNNs layers $T$ and
SNNs layers $N$, as long as $\Delta t\sim T^{-1}$ and $\Delta\tau\sim N^{-1}$.
This effectively addresses the issues of exploding and vanishing gradients.

##  5 Experiments

To evaluate the effectiveness of our proposed COS-GNN, we conduct extensive
experiments with COS-GNN across various graph learning tasks, including node
classification and graph classification. We assess the method in two settings:
COS-GNN-1st, using first-order SNNs and second-order CGNNs, and COS-GNN-2nd,
employing second-order SNNs and second-order CGNNs.

###  5.1 Experimental Settings

Datasets. For the node classification, we evaluate COS-GNN on homophilic
(i.e., Cora (McCallum et al., 2000), Citeseer (Sen et al., 2008) and Pubmed
(Namata et al., 2012)) and heterophilic (i.e., Texas, Wisconsin and Cornell
from the
WebKB111http://www.cs.cmu.edu/afs/cs.cmu.edu/project/theo-11/www/wwkb/)
datasets, where high homophily indicates that a node’s features are similar to
those of its neighbors, and heterophily suggests the opposite. The homophily
level is measured according to (Pei et al., 2020), and is reported in Table 1
and 2. In the graph classification task, we utilize the MNIST dataset (LeCun
et al., 1998). To represent the grey-scale images as irregular graphs, we
associate each superpixel (large blob of similar color) with a vertex, and the
spatial adjacency between superpixels with edges. Each graph consists of a
fixed number of 75 superpixels (vertices). To ensure consistent evaluation, we
adopt the standard splitting of 55K-5K-10K for training, validation, and
testing purposes (Rusch et al., 2022).

Table 1: The test accuracy (in %) for node classification on homophilic datasets. The results are calculated by averaging the results of 20 random initializations across 5 random splits. The mean and standard deviation of these results are obtained. Bold numbers means the best performance, and underline numbers indicates the second best performance. | Cora | Citeseer | Pubmed
---|---|---|---
Homophily level | 0.81 | 0.74 | 0.80
GAT-ppr |  81.6$\pm$0.3 |  68.5$\pm$0.2 |  76.7$\pm$0.3
MoNet |  81.3$\pm$1.3 |  71.2$\pm$2.0 |  78.6$\pm$2.3
GraphSage-mean |  79.2$\pm$7.7 |  71.6$\pm$1.9 |  77.4$\pm$2.2
GraphSage-maxpool |  76.6$\pm$1.9 |  67.5$\pm$2.3 |  76.1$\pm$2.3
CGNN |  81.4$\pm$1.6 |  66.9$\pm$1.8 |  66.6$\pm$4.4
GDE |  78.7$\pm$2.2 |  71.8$\pm$1.1 |  73.9$\pm$3.7
GCN |  81.5$\pm$1.3 |  71.9$\pm$1.9 |  77.8$\pm$2.9
GAT |  81.8$\pm$1.3 |  71.4$\pm$1.9 |  78.7$\pm$2.3
SGC |  81.5$\pm$0.4 |  71.7$\pm$0.4 |  79.2$\pm$0.3
GRAND |  83.6$\pm$1.0 |  73.4$\pm$0.5 |  78.8$\pm$1.7
GraphCON-GCN |  81.9$\pm$1.7 |  72.9$\pm$2.1 |  78.8$\pm$2.6
GraphCON-GAT |  83.2$\pm$1.4 |  73.2$\pm$1.8 | 79.5$\pm$1.8
SpikingGCN |  80.7$\pm$0.6 |  72.5$\pm$0.2 |  77.6$\pm$0.5
COS-GNN-1st | 83.3$\pm$2.1 | 73.7$\pm$2.0 |  76.9$\pm$2.7
COS-GNN-2nd | 83.7$\pm$1.3 | 75.2$\pm$2.0 | 79.6$\pm$2.3

Baselines. For the homophilic datasets, we use standard GNN baselines: GCN
(Kipf & Welling, 2017), SGC (Wu et al., 2019), GAT (Velickovic et al., 2017),
MoNet (Monti et al., 2017), GraphSage (Hamilton et al., 2017), CGNN (Xhonneux
et al., 2020), GDE (Poli et al., 2019), GRAND (Chamberlain et al., 2021),
GraphCON (Rusch et al., 2022) and SpikingGCN (Zhu et al., 2022). Due to the
assumption on neighbor feature similarity does not hold in the heterophilic
datasets, we utilize additional specific GNN methods as baselines: GPRGNN
(Chien et al., 2020), H2GCN (Zhu et al., 2020), GCNII (Chen et al., 2020),
Geom-GCN (Pei et al., 2020) and PairNorm (Zhao & Akoglu, 2019). For the graph
classification task, we apply ChebNet (Defferrard et al., 2016), PNCNN (Finzi
et al., 2021), SplineCNN (Fey et al., 2018), GIN (Xu et al., 2019), and
GatedGCN (Bresson & Laurent, 2017) for comparison.

Table 2: The test accuracy (in %) for node classification on heterophilic datasets. All results represent the average performance of the respective model over 10 fixed train/val/test splits. Bold numbers means the best performance, and underline numbers indicates the second best performance. | Texas | Wisconsin | Cornell
---|---|---|---
Homophily level | 0.11 | 0.21 | 0.30
GPRGNN |  78.4$\pm$4.4 |  82.9$\pm$4.2 |  80.3$\pm$8.1
H2GCN |  84.9$\pm$7.2 |  87.7$\pm$5.0 |  82.7$\pm$5.3
GCNII |  77.6$\pm$3.8 |  80.4$\pm$3.4 |  77.9$\pm$3.8
Geom-GCN |  66.8$\pm$2.7 |  64.5$\pm$3.7 |  60.5$\pm$3.7
PairNorm |  60.3$\pm$4.3 |  48.4$\pm$6.1 |  58.9$\pm$3.2
GraphSAGE |  82.4$\pm$6.1 |  81.2$\pm$5.6 |  76.0$\pm$5.0
MLP |  80.8$\pm$4.8 |  85.3$\pm$3.3 |  81.9$\pm$6.4
GCN |  55.1$\pm$5.2 |  51.8$\pm$3.1 |  60.5$\pm$5.3
GAT |  52.2$\pm$6.6 |  49.4$\pm$4.1 |  61.9$\pm$5.1
GraphCON-GCN | 85.4$\pm$4.2 | 87.8$\pm$3.3 | 84.3$\pm$4.8
GraphCON-GAT |  82.2$\pm$4.7 |  85.7$\pm$3.6 |  83.2$\pm$7.0
COS-GNN-1st |  81.6$\pm$6.2 |  84.9$\pm$3.2 |  80.4$\pm$1.9
COS-GNN-2nd | 87.3$\pm$4.2 | 88.8$\pm$2.5 | 83.7$\pm$2.7

Implementation Details. For the homophilic node classification task, we report
the average results of 20 random initialization across 5 random splits. For
the heterophilic node classification task, we present the average performance
of the respective model over 10 fixed train/val/test splits. The results of
baselines are reported in (Rusch et al., 2022). For COS-GNN-1st, we set the
hyperparameter $\lambda$ to 1. As for COS-GNN-2nd, we set the hyperparameters
$\alpha$ and $\beta$ to 1 as default. The time latency $N$ in SNNs are set to
8. For all the methods, we set the hidden size to 64 and the learning rate to
0.001 as default.

Table 3: The test accuracy (in %) for graph classification on MNIST datasets. Bold numbers means the best performance, and underline numbers indicates the second best performance. Model | Test accuracy
---|---
ChebNet (Defferrard et al., 2016) | 75.62
MoNet (Monti et al., 2017) | 91.11
PNCNN (Finzi et al., 2021) | 98.76
SplineCNN (Fey et al., 2018) | 95.22
GIN (Xu et al., 2019) | 97.23
GatedGCN (Bresson & Laurent, 2017) | 97.95
GCN (Kipf & Welling, 2017) | 88.89
GAT (Velickovic et al., 2017) | 96.19
GraphCON-GCN (Rusch et al., 2022) | 98.68
GraphCON-GAT (Rusch et al., 2022) | 98.91
COS-GNN-1st | 98.82
COS-GNN-2nd | 98.92
Table 4: Ablation results. Bold numbers means the best performance. | Cora | Citeseer | Pubmed | Texax | Wisconsin | Cornell | Avg.
---|---|---|---|---|---|---|---
Homophily level | 0.81 | 0.74 | 0.80 | 0.11 | 0.21 | 0.3
COS-GNN-1st-2nd | 83.2$\pm$1.4 | 74.1$\pm$1.4 | 76.3$\pm$2.2 | 81.7$\pm$3.9 | 85.1$\pm$2.8 | 81.0$\pm$1.9 | 80.2
COS-GNN-2nd-1st | 83.5$\pm$1.8 | 73.4$\pm$2.1 | 77.2$\pm$2.3 | 83.1$\pm$3.8 | 84.4$\pm$2.2 | 81.2$\pm$2.7 | 80.5
COS-GNN-1st | 83.3$\pm$2.1 | 73.7$\pm$2.0 | 76.9$\pm$2.7 | 81.6$\pm$6.2 | 84.9$\pm$3.2 | 80.4$\pm$1.9 | 80.1
COS-GNN-2nd | 83.7$\pm$1.3 | 75.2$\pm$2.0 | 79.6$\pm$2.3 | 87.3$\pm$4.2 | 88.8$\pm$2.5 | 83.7$\pm$2.7 | 83.1

###  5.2 Performance Comparision

Homophilic Node Classification. Table 1 shows the results of the proposed COS-
GNN with the comparison of baselines. From the results, we find that: (1)
Compared with the discrete methods (i.e., the baselines excluding GraphCON),
the continuous methods (GraphCON and COS-GNN) achieve the best and second best
performance, indicating that the continuous methods would help to capture the
continuous changes and subtle dynamics from graphs. (2) COS-GNN-1st and COS-
GNN-2nd outperforms other baselines in most cases. We attribute that, even if
SNNs loses some detailed information, COS-GNN can still achieve good
performance on the relatively simple homophilic dataset. Furthermore, the
application of SNNs contributes to improved efficiency in the COS-GNN
framework. (3) COS-GNN-2nd consistently outperforms the COS-GNN-1st. This
highlights the significance of introducing high-order structures to preserve
information and mitigates the information loss issue caused by first-order
SNNs. Although high-order structures suffer higher energy costs compared to
first-order, the performance gains make it worthwhile to deploy them. (4) COS-
GNN-1st and COS-GNN-2nd outperforms the spiking-based method (i.e., Spiking)
in most case. This can be attributed to the utilization of CGNNs, which
efficiently capture the spatio-temporal relationships within dynamical
systems.

Heterophilic Node Classification. Table 2 shows the results of heterophilic
node classification, and we observe that: (1) The traditional message-passing-
based methods (GCN, GAT, GraphSAGE and Geom-GCN) perform worse than the well-
designed methods (GPRGNN, H2GCN, GCNII, GraphCON and COS-GNN) for heterophilic
datasets. This disparity comes from the inaccurate assumption of neighbor
feature similarity, which doesn’t hold in heterophilic datasets. The
propagation of heterophilic information between nodes would degrade the
model’s representation ability, leading to a decline in performance. (2) The
COS-GNN-1st performs less effectively than GraphCON. This is because node
prediction tasks on heterophilic datasets are more influenced by the
characteristics of heterophilic features compared to homophilic datasets.
Consequently, the information loss issue caused by first-order SNNs results in
worse model performance. (3) The COS-GNN-2nd consistently outperforms COS-
GNN-1st, providing further evidence of the effectiveness of high-order
structures in preserving information and mitigating the issue of information
loss.

Table 5: Operations comparison on different datasets. models | Cora | Citeseer | Pubmed
---|---|---|---
GCN | 67.77K | 79.54K | 414.16K
GAT | 308.94K | 349.91K | 1.53M
SGC | 10.03K | 22.22K | 1.50K
COS-GNN-1st | 2.33K | 1.94K | 1.02K
COS-GNN-2nd | 3.97K | 3.42K | 2.78K

Graph Classification. We present the graph classification results of our
proposed COS-GNN alongside comparison baselines in Table 3. From Table 3, we
find that: (1) In the graph classification tasks, CGNNs methods (i.e., COS-GNN
and GraphCON) consistently outperform the baseline methods across all cases.
This underscores the importance of employing a continuous processing approach
when dealing with graph data, enabling the extraction of continuous changes
and subtle dynamics from graphs. (2) The COS-GNN-1st performs worse than the
COS-GNN-2nd, highlighting the significance of incorporating high-order
structures to obtain additional information for prediction, without incurring
significant overhead. (3) The COS-GNN-1st performs worse than GraphCON-GAT and
better than GraphCON-GCN. Compared to GraphCON-GCN, the information loss
caused by SNNs does not critically affect graph representation ability. On the
contrary, the binarization operation of SNNs contributes to reduced energy
consumption. Graph-GAT outperforms COS-GNN-1st, mainly because the GAT method
enhances graph representation. However, Graph-GAT still lags behind COS-
GNN-2nd, indicating that the introduction of high-order structures mitigates
the information loss issue associated with first-order methods.

###  5.3 Ablation Study

We conducted ablation studies to assess the contributions of different
components using two variants, and the results are presented in Table 4.
Specifically, we introduced two model variants: (1) COS-GNN-1st-2nd, which
utilizes the first-order SNNs and second-order CGNNs, and (2) COS-GNN-2nd-1st,
incorporating the second-order SNNs and first-order CGNNs. Table 4 shows that
(1) COS-GNN-2nd consistently outperforms other variations, while COS-
GNN-1st-2nd yields the worst performance. This is because the issue of
information loss is crucial for graph representation, and the incorporation of
high-order SNNs assists in preserving more information, consequently achieving
superior results. (2) In most cases, COS-GNN-2nd-1st outperforms both COS-
GNN-1st and COS-GNN-1st-2nd, suggesting that, compared to the capability of
CGNNs in capturing dynamic node relationships, the ability to mitigate the
issue of information loss is more important.

###  5.4 Energy Efficiency Analysis

To evaluate the energy efficiency of COS-GNN, we apply the metric in (Zhu et
al., 2022), which quantifies the number of operations needed for node
prediction. The results are reported in Table 5. From the results, we find
that the COS-GNN-1st and COS-GNN-2nd have a significant operation reduction
compared with baselines. The conventional operation unit for ANNs on GPUs is
commonly configured as multiply-accumulate (MAC), whereas for SNNs on
neuromorphic chips, it takes the form of synaptic operation (SOP). The SOP is
defined as the change in membrane potential in the LIF nodes. According to (Hu
et al., 2021; Kim et al., 2020), the SOP of SNNs consumes significantly less
energy than the MAC in ANNs, which further underscores the energy efficiency
of COS-GNN.

###  5.5 Sensitivity Analysis

In this part, we examine the sensitivity of the proposed COS-GNN to its
hyperparameters, specifically the time latency parameter ($N$) in SNNs, which
plays a crucial role in the model’s performance. $N$ controls the number of
SNNs propagation steps and is directly related to the training complexity.
Figure 2 shows the results of $N$ across different datasets. We initially vary
the parameter $N$ within the range of $\\{5,6,7,8,9,10,11\\}$ while keeping
other parameters fixed. From the results, we find that, the performance
exhibits a increasing trend initially, followed by stabilization as the value
of $N$ increases. Typically, in SNNs, spiking signals are integrated with
historical information at each time latency. Smaller values of $N$ result in
less information available for graph representation, degrading the
performance. However, large values of $N$ increase model complexity during
training. Striking a balance between model performance and complexity, we set
$N$ to 8 as default.

Figure 2: Sensitivity analysis on time latency $N$ in SNNs across various
datasets. The solid line denotes the results of COS-GNN-1st, and the dotted
line denotes the COS-GNN-2nd.

##  6 Conclusion

In this paper, we address the practical problem of continuous spiking graph
representation learning and propose an effective method named COS-GNN. COS-GNN
integrates SNNs and CGNNs into a unified framework from two distinct time
dimensions, thus retaining the benefits of low-power consumption and fine-
grained feature extraction. Considering that the high-order structure would
help to relieve the problem of information loss, we derive the second-order
spike representation and investigate the backpropagation of second-order SNNs,
and subsequently introduced the second-order COS-GNN. To ensure the stability
of COS-GNN, we further prove that COS-GNN mitigates the gradient exploding and
vanishing problem. Extensive experiments on diverse datasets validate the
efficacy of proposed COS-GNN compared with various competing methods. In
future work, we will explore the higher-order structure for more efficient
continuous graph learning.

##  7 Impact Statements

This work introduces an innovative approach for continuous spiking graph
neural network, with the objective of advancing the machine learning field,
particularly in the domain of graph neural networks. The proposed method has
the potential to substantially enhance the efficiency and scalability of graph
learning tasks. The societal implications of this research are multifaceted.
The introduced method has the capacity to contribute to the development of
more efficient and effective machine learning systems, with potential
applications across various domains, including healthcare, education, and
technology. Such advancements could lead to improved services and products,
ultimately benefiting society as a whole.

## References

  * Battaglia et al. (2016) Battaglia, P., Pascanu, R., Lai, M., Jimenez Rezende, D., et al.  Interaction networks for learning about objects, relations and physics.  _Proceedings of the Conference on Neural Information Processing Systems_ , 2016.
  * Bellec et al. (2018) Bellec, G., Salaj, D., Subramoney, A., Legenstein, R., and Maass, W.  Long short-term memory and learning-to-learn in networks of spiking neurons.  In _Proceedings of the Conference on Neural Information Processing Systems_ , pp. 795–805, 2018.
  * Bohte et al. (2000) Bohte, S., Kok, J., and La Poutre, H.  Spikeprop: Backpropagation for networks of spiking neurons error-backpropagation in a network of spik-ing neurons.  In _ESANN_ , 2000.
  * Bresson & Laurent (2017) Bresson, X. and Laurent, T.  Residual gated graph convnets.  _arXiv preprint arXiv:1711.07553_ , 2017.
  * Brette et al. (2007) Brette, R., Rudolph, M., Carnevale, T., Hines, M., Beeman, D., Bower, J. M., Diesmann, M., Morrison, A., Goodman, P. H., Harris, F. C., et al.  Simulation of networks of spiking neurons: a review of tools and strategies.  _Journal of computational neuroscience_ , 23:349–398, 2007.
  * Chamberlain et al. (2021) Chamberlain, B., Rowbottom, J., Gorinova, M. I., Bronstein, M., Webb, S., and Rossi, E.  Grand: Graph neural diffusion.  In _Proceedings of the International Conference on Machine Learning_ , pp. 1407–1418, 2021.
  * Chen et al. (2020) Chen, M., Wei, Z., Huang, Z., Ding, B., and Li, Y.  Simple and deep graph convolutional networks.  In _Proceedings of the International Conference on Machine Learning_ , pp. 1725–1735, 2020.
  * Chen et al. (2018) Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K.  Neural ordinary differential equations.  In _Proceedings of the Conference on Neural Information Processing Systems_ , 2018.
  * Chien et al. (2020) Chien, E., Peng, J., Li, P., and Milenkovic, O.  Adaptive universal generalized pagerank graph neural network.  In _Proceedings of the International Conference on Learning Representations_ , 2020.
  * Defferrard et al. (2016) Defferrard, M., Bresson, X., and Vandergheynst, P.  Convolutional neural networks on graphs with fast localized spectral filtering.  In _Proceedings of the Conference on Neural Information Processing Systems_ , 2016.
  * Eshraghian et al. (2023) Eshraghian, J. K., Ward, M., Neftci, E. O., Wang, X., Lenz, G., Dwivedi, G., Bennamoun, M., Jeong, D. S., and Lu, W. D.  Training spiking neural networks using lessons from deep learning.  _Proceedings of the IEEE_ , 2023.
  * Esser et al. (2015) Esser, S. K., Appuswamy, R., Merolla, P., Arthur, J. V., and Modha, D. S.  Backpropagation for energy-efficient neuromorphic computing.  In _Proceedings of the Conference on Neural Information Processing Systems_ , pp. 1117–1125, 2015.
  * Fey et al. (2018) Fey, M., Lenssen, J. E., Weichert, F., and Müller, H.  Splinecnn: Fast geometric deep learning with continuous b-spline kernels.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 869–877, 2018.
  * Finzi et al. (2021) Finzi, M. A., Bondesan, R., and Welling, M.  Probabilistic numeric convolutional neural networks.  In _Proceedings of the International Conference on Learning Representations_ , 2021.
  * Gupta et al. (2022) Gupta, J., Vemprala, S., and Kapoor, A.  Learning modular simulations for homogeneous systems.  In _Proceedings of the Conference on Neural Information Processing Systems_ , pp. 14852–14864, 2022.
  * Hafiene et al. (2020) Hafiene, N., Karoui, W., and Romdhane, L. B.  Influential nodes detection in dynamic social networks: A survey.  _Expert Systems with Applications_ , 159:113642, 2020.
  * Hairer et al. (1993) Hairer, E., Norsett, S., and Wanner, G.  _Solving Ordinary Differential Equations I: Nonstiff Problems_ , volume 8.  01 1993.  ISBN 978-3-540-56670-0.  doi: 10.1007/978-3-540-78862-1.
  * Hamilton et al. (2017) Hamilton, W., Ying, Z., and Leskovec, J.  Inductive representation learning on large graphs.  In _Proceedings of the Conference on Neural Information Processing Systems_ , 2017.
  * Hsieh et al. (2021) Hsieh, T.-Y., Wang, S., Sun, Y., and Honavar, V.  Explainable multivariate time series classification: a deep neural network which learns to attend to important variables as well as time intervals.  In _Proceedings of the International ACM Conference on Web Search & Data Mining_, pp. 607–615, 2021.
  * Hu et al. (2021) Hu, Y., Tang, H., and Pan, G.  Spiking deep residual networks.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2021.
  * Huang et al. (2020) Huang, Z., Sun, Y., and Wang, W.  Learning continuous system dynamics from irregularly-sampled partial observations.  In _Proceedings of the Conference on Neural Information Processing Systems_ , volume 33, pp. 16177–16187, 2020.
  * Huang et al. (2021) Huang, Z., Sun, Y., and Wang, W.  Coupled graph ode for learning interacting system dynamics.  In _Proceedings of the International ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, pp. 705–715, 2021.
  * Ju et al. (2024) Ju, W., Yi, S., Wang, Y., Xiao, Z., Mao, Z., Li, H., Gu, Y., Qin, Y., Yin, N., Wang, S., et al.  A survey of graph neural networks in real world: Imbalance, noise, privacy and ood challenges.  _arXiv preprint arXiv:2403.04468_ , 2024.
  * Kim et al. (2020) Kim, S., Park, S., Na, B., and Yoon, S.  Spiking-yolo: spiking neural network for energy-efficient object detection.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , pp. 11270–11277, 2020.
  * Kipf et al. (2018) Kipf, T., Fetaya, E., Wang, K.-C., Welling, M., and Zemel, R.  Neural relational inference for interacting systems.  In _Proceedings of the International Conference on Machine Learning_ , pp. 2688–2697, 2018.
  * Kipf & Welling (2017) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  In _Proceedings of the International Conference on Learning Representations_ , 2017.
  * LeCun et al. (1998) LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * Li et al. (2023) Li, J., Yu, Z., Zhu, Z., Chen, L., Yu, Q., Zheng, Z., Tian, S., Wu, R., and Meng, C.  Scaling up dynamic graph representation learning via spiking neural networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , pp. 8588–8596, 2023.
  * Liao et al. (2021) Liao, S., Liang, S., Meng, Z., and Zhang, Q.  Learning dynamic embeddings for temporal knowledge graphs.  In _Proceedings of the International ACM Conference on Web Search & Data Mining_, pp. 535–543, 2021.
  * Luo et al. (2023) Luo, X., Yuan, J., Huang, Z., Jiang, H., Qin, Y., Ju, W., Zhang, M., and Sun, Y.  Hope: High-order graph ode for modeling interacting dynamics.  In _Proceedings of the International Conference on Machine Learning_ , pp. 23124–23139, 2023.
  * Maass (1997) Maass, W.  Networks of spiking neurons: the third generation of neural network models.  _Neural networks_ , 10(9):1659–1671, 1997.
  * McCallum et al. (2000) McCallum, A. K., Nigam, K., Rennie, J., and Seymore, K.  Automating the construction of internet portals with machine learning.  _Information Retrieval_ , 3:127–163, 2000.
  * Meng et al. (2022) Meng, Q., Xiao, M., Yan, S., Wang, Y., Lin, Z., and Luo, Z.-Q.  Training high-performance low-latency spiking neural networks by differentiation on spike representation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 12444–12453, 2022.
  * Monti et al. (2017) Monti, F., Boscaini, D., Masci, J., Rodola, E., Svoboda, J., and Bronstein, M. M.  Geometric deep learning on graphs and manifolds using mixture model cnns.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 5115–5124, 2017.
  * Namata et al. (2012) Namata, G., London, B., Getoor, L., Huang, B., and Edu, U.  Query-driven active surveying for collective classification.  In _10th International Workshop on Mining and Learning with Graphs_ , volume 8, pp. 1, 2012.
  * Pang et al. (2023) Pang, J., Wang, Z., Tang, J., Xiao, M., and Yin, N.  Sa-gda: Spectral augmentation for graph domain adaptation.  In _Proceedings of the ACM International Conference on Multimedia_ , pp. 309–318, 2023.
  * Pei et al. (2020) Pei, H., Wei, B., Chang, K. C.-C., Lei, Y., and Yang, B.  Geom-gcn: Geometric graph convolutional networks.  _arXiv preprint arXiv:2002.05287_ , 2020.
  * Poli et al. (2019) Poli, M., Massaroli, S., Park, J., Yamashita, A., Asama, H., and Park, J.  Graph neural ordinary differential equations.  _arXiv preprint arXiv:1911.07532_ , 2019.
  * Rathi et al. (2020) Rathi, N., Srinivasan, G., Panda, P., and Roy, K.  Enabling deep spiking neural networks with hybrid conversion and spike timing dependent backpropagation.  In _Proceedings of the International Conference on Machine Learning_ , 2020.
  * Rueckauer et al. (2017) Rueckauer, B., Lungu, I.-A., Hu, Y., Pfeiffer, M., and Liu, S.-C.  Conversion of continuous-valued deep networks to efficient event-driven networks for image classification.  _Frontiers in neuroscience_ , 11:682, 2017.
  * Rusch et al. (2022) Rusch, T. K., Chamberlain, B., Rowbottom, J., Mishra, S., and Bronstein, M.  Graph-coupled oscillator networks.  In _Proceedings of the International Conference on Machine Learning_ , pp. 18888–18909, 2022.
  * Sen et al. (2008) Sen, P., Namata, G., Bilgic, M., Getoor, L., Galligher, B., and Eliassi-Rad, T.  Collective classification in network data.  _AI magazine_ , 29(3):93–93, 2008.
  * Velickovic et al. (2017) Velickovic, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., Bengio, Y., et al.  Graph attention networks.  In _Proceedings of the International Conference on Learning Representations_ , 2017.
  * Wu et al. (2019) Wu, F., Souza, A., Zhang, T., Fifty, C., Yu, T., and Weinberger, K.  Simplifying graph convolutional networks.  In _Proceedings of the International Conference on Machine Learning_ , pp. 6861–6871, 2019.
  * Xhonneux et al. (2020) Xhonneux, L.-P., Qu, M., and Tang, J.  Continuous graph neural networks.  In _Proceedings of the International Conference on Machine Learning_ , 2020.
  * Xu et al. (2019) Xu, K., Hu, W., Leskovec, J., and Jegelka, S.  How powerful are graph neural networks?  In _Proceedings of the International Conference on Learning Representations_ , 2019.
  * Yin et al. (2022a) Yin, N., Feng, F., Luo, Z., Zhang, X., Wang, W., Luo, X., Chen, C., and Hua, X.-S.  Dynamic hypergraph convolutional network.  In _Proceedings of the IEEE Conference on Data Engineering_ , pp. 1621–1634, 2022a.
  * Yin et al. (2022b) Yin, N., Shen, L., Li, B., Wang, M., Luo, X., Chen, C., Luo, Z., and Hua, X.-S.  Deal: An unsupervised domain adaptive framework for graph-level classification.  In _Proceedings of the ACM International Conference on Multimedia_ , pp. 3470–3479, 2022b.
  * Yin et al. (2023a) Yin, N., Shen, L., Wang, M., Lan, L., Ma, Z., Chen, C., Hua, X.-S., and Luo, X.  Coco: A coupled contrastive framework for unsupervised domain adaptive graph classification.  In _Proceedings of the International Conference on Machine Learning_ , pp. 40040–40053, 2023a.
  * Yin et al. (2023b) Yin, N., Shen, L., Xiong, H., Gu, B., Chen, C., Hua, X.-S., Liu, S., and Luo, X.  Messages are never propagated alone: Collaborative hypergraph neural network for time-series forecasting.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2023b.
  * Yin et al. (2023c) Yin, N., Wang, M., Chen, Z., De Masi, G., Gu, B., and Xiong, H.  Dynamic spiking graph neural networks.  _arXiv preprint arXiv:2401.05373_ , 2023c.
  * Yin et al. (2024a) Yin, N., Wang, M., Chen, Z., De Masi, G., Xiong, H., and Gu, B.  Dynamic spiking graph neural networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , pp. 16495–16503, 2024a.
  * Yin et al. (2024b) Yin, N., Wang, M., Chen, Z., Shen, L., Xiong, H., Gu, B., and Luo, X.  DREAM: Dual structured exploration with mixup for open-set graph domain adaption.  In _Proceedings of the International Conference on Learning Representations_ , 2024b.
  * Zhang et al. (2022a) Zhang, D., Zhang, T., Jia, S., Wang, Q., and Xu, B.  Recent advances and new frontiers in spiking neural networks.  In _Proceedings of the International Joint Conference on Artificial Intelligence_ , pp. 5670–5677, 2022a.
  * Zhang et al. (2022b) Zhang, Y., Gao, S., Pei, J., and Huang, H.  Improving social network embedding via new second-order continuous graph neural networks.  In _Proceedings of the International ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, pp. 2515–2523, 2022b.
  * Zhao & Akoglu (2019) Zhao, L. and Akoglu, L.  Pairnorm: Tackling oversmoothing in gnns.  _arXiv preprint arXiv:1909.12223_ , 2019.
  * Zhu et al. (2020) Zhu, J., Yan, Y., Zhao, L., Heimann, M., Akoglu, L., and Koutra, D.  Generalizing graph neural networks beyond homophily.  In _Proceedings of the Conference on Neural Information Processing Systems_ , 2020.
  * Zhu et al. (2022) Zhu, Z., Peng, J., Li, J., Chen, L., Yu, Q., and Luo, S.  Spiking graph convolutional networks.  In _Proceedings of the International Joint Conference on Artificial Intelligence_ , pp. 2434–2440, 2022.

##  Appendix A Proof of Proposition 4.1

Proposition 4.1 Define the first-order SNNs ODE as
$\frac{du_{t}^{\tau}}{d\tau}=g(u_{t}^{\tau},\tau)$, and first-order Graph ODE
as $\frac{du_{t}^{\tau}}{dt}=f(u^{\tau}_{t},t)$, then the first-order graph
PDE network can be formulated as:

$u_{t+1}^{\tau+1}=u_{0}^{0}+\int_{0}^{T}f\left(u^{0}_{y}+\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy+\int_{0}^{N}g\left(u_{0}^{x}+\int_{0}^{T-1}f(u_{y}^{x},y)dy\right)dx.$

Proof.

$\frac{du_{t}^{\tau}}{d\tau}=g(u_{t}^{\tau},\tau),\quad\frac{du_{t}^{\tau}}{dt}=f(u^{\tau}_{t},t),$

$u_{t}^{\tau}$ is a function related to variable $t$ and $\tau$, we have
$\frac{\partial u_{t}^{\tau}}{\partial\tau}=g(u_{t}^{\tau})$ and
$\frac{\partial u_{t}^{\tau}}{\partial t}=f(u_{t}^{\tau})$. Thus,

$u_{t}^{\tau+1}=u_{t}^{\tau}+\int_{\tau}^{\tau+1}g(u_{t}^{x},x)dx,\;\;u_{t+1}^{\tau+1}=u_{t}^{\tau+1}+\int_{t}^{t+1}f(u^{\tau+1}_{y},y)dy,$
(15) $\displaystyle u_{T}^{N}=$ $\displaystyle
u_{T-1}^{N-1}+\int_{T-1}^{T}f(u^{N}_{y},y)dy+\int_{N-1}^{N}g(u_{T-1}^{x},x)dx$
(16) $\displaystyle=$ $\displaystyle
u_{T-2}^{N-2}+\int_{T-2}^{T}f\left(u^{N}_{y},y\right)dy+\int_{N-2}^{N}g(u_{T-1}^{x},x)dx$
$\displaystyle=$ $\displaystyle
u_{0}^{0}+\int_{0}^{T}f\left(u^{N}_{y},y\right)dy+\int_{0}^{N}g(u_{T-1}^{x},x)dx$
$\displaystyle=$ $\displaystyle
u_{0}^{0}+\int_{0}^{T}f\left(u^{N-1}_{y}+\int_{N-1}^{N}g(u_{y}^{x},x)dx\right)dy+\int_{0}^{N}g\left(u_{T-2}^{x}+\int_{T-2}^{T-1}f(u_{y}^{x},y)dy\right)dx$
$\displaystyle=$ $\displaystyle
u_{0}^{0}+\int_{0}^{T}f\left(u^{0}_{y}+\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy+\int_{0}^{N}g\left(u_{0}^{x}+\int_{0}^{T-1}f(u_{y}^{x},y)dy\right)dx.$

By adding the initial state on each time step and latency with $u_{t}^{0}=0$
and $u_{0}^{\tau}=0$, we have:

$\displaystyle u_{T}^{N}=$
$\displaystyle\int_{0}^{T}f\left(\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy+\int_{0}^{N}g\left(\int_{0}^{T-1}f(u_{y}^{x},y)dy\right)dx$
(17) $\displaystyle=$
$\displaystyle\underbrace{\int_{0}^{T-1}f\left(\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy+\int_{0}^{N}g\left(\int_{0}^{T-1}f(u_{y}^{x},y)dy\right)dx}_{first\;\;term}+\underbrace{\int_{T-1}^{T}f\left(\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy}_{second\;\;term}$
$\displaystyle=$ $\displaystyle
2\int_{0}^{T-1}f\left(\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy+\int_{T-1}^{T}f\left(\int_{0}^{N}g(u_{y}^{x},x)dx\right)dy.$

The first term denotes that the SNNs and CGNNs are interactively updated
during the time step $0$ to $T-1$, and the second term denotes that at the
last time step $T$, COS-GNN simply calculates the CGNNs process while ignoring
the SNNs for prediction.

##  Appendix B Proof of Eqn. 11

Proof. From Eq. 9, we have:

$u(\tau+1)=\beta^{2}u(\tau-k+1)+\alpha\sum_{i=0}^{k-1}\beta^{i}I_{syn}(\tau-i)+\sum_{i=0}^{k-1}\beta^{i}(I_{input}(\tau-i)-V_{th}s(\tau-i)),$
(18)
$u(N)=\alpha\sum_{n=0}^{N-1}\beta^{n}I_{syn}(N-n-1)+\sum_{n=0}^{N-1}\beta^{n}(I_{input}(N-n-1)-V_{th}s(N-n-1)).$
(19)

Due to:

$I_{syn}(\tau+1)=\alpha^{k}I_{syn}(\tau-k+1)+\sum_{i=0}^{k}\alpha^{i}I_{input}(\tau-i),$
(20)

we have,

$\displaystyle u(N)=$
$\displaystyle\alpha\sum_{n=0}^{N-1}\beta^{N-n-1}I_{syn}(n)+\sum_{n=0}^{N-1}\beta^{N-n-1}\left(I_{input}(n)-V_{th}s(n)\right)$
$\displaystyle=$
$\displaystyle\alpha\left(\left(\frac{\beta^{N-1}\alpha^{-1}\left(1-(\frac{\alpha}{\beta})^{N}\right)}{{1-\frac{\alpha}{\beta}}}\right)I_{in}(0)+\left(\frac{\beta^{N-2}\alpha^{-1}\left(1-(\frac{\alpha}{\beta})^{N-1}\right)}{{1-\frac{\alpha}{\beta}}}\right)I_{in}(1)+\cdots\right.$
$\displaystyle\left.+\left(\frac{\beta^{N-i}\alpha^{-1}\left(1-(\frac{\alpha}{\beta})^{N-i+1}\right)}{{1-\frac{\alpha}{\beta}}}\right)I_{in}(i-1)+\cdots+(\beta^{2}\alpha^{-1}+\beta+\alpha)I_{in}(N-3)\right.$
$\displaystyle\left.+(\beta\alpha^{-1}+1)I_{in}(N-2)+\alpha^{-1}I_{in}(N-1)\right)-\sum_{n=0}^{N-1}\beta^{N-n-1}V_{th}s(n)$
$\displaystyle=$
$\displaystyle\frac{1}{\beta-\alpha}\left(\left(\beta^{N}\left(1-\left(\frac{\alpha}{\beta}\right)^{N}\right)I_{in}(0)\right)+\dots+\left(\beta^{N-i+1}\left(1-\left(\frac{\alpha}{\beta}\right)^{N-i+1}\right)I_{in}(i-1)\right)\right.$
$\displaystyle\left.+\cdots+(\beta-\alpha)I_{in}(N-1)\right)-\sum_{n=0}^{N-1}\beta^{N-n-1}V_{th}s(n)$
$\displaystyle=$
$\displaystyle\frac{1}{\beta-\alpha}\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})I_{in}(n)-\sum_{n=0}^{N-1}\beta^{N-n-1}V_{th}s(n).$

Define
$\hat{I}(N)=\frac{1}{(\beta-\alpha)^{2}}\frac{\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})I_{in}(n)}{\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})}$,
and
$\hat{a}(N)=\frac{1}{\beta^{2}}\frac{V_{th}\sum_{n=0}^{N-1}\beta^{N-n}s(n)}{\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})}$,
we have:

$\hat{a}(N)=\frac{\beta-\alpha}{\beta}\frac{\hat{I}(N)}{\Delta\tau}-\frac{u(N)}{\Delta\tau\beta\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})}\approx\frac{\tau_{syn}\tau_{mem}}{\tau_{mem}-\tau_{syn}}\hat{I}(N)-\frac{u(N)}{\Delta\tau\beta\sum_{n=0}^{N-1}(\beta^{N-n}-\alpha^{N-n})},$

where $\alpha=exp(-\Delta\tau/\tau_{syn})$,
$\beta=exp(-\Delta\tau/\tau_{mem})$.

##  Appendix C Proof of Proposition 4.3

Proposition 4.3 Define the second-order SNNs as
$\frac{d^{2}u_{t}^{\tau}}{d\tau^{2}}+\delta\frac{du_{t}^{\tau}}{d\tau}=g(u^{\tau}_{t},\tau)$,
and second-order CGNNs as
$\frac{d^{2}u_{t}^{\tau}}{dt^{2}}+\gamma\frac{du_{t}^{\tau}}{dt}=f(u^{\tau}_{t},t)$,
then the second-order COS-GNN is formulated as:

$u_{t}^{\tau}=\int_{0}^{T}h\left(\int_{0}^{N}e(u_{t}^{\tau})d\tau\right)dt=\int_{0}^{N}e\left(\int_{0}^{T}h(u_{t}^{\tau})dt\right)d\tau,$
$\frac{\partial^{2}u_{t}^{\tau}}{\partial\tau^{2}}+\delta\frac{\partial
u_{t}^{\tau}}{\partial\tau}=g(u_{t}^{\tau}),\quad\frac{\partial^{2}u_{t}^{\tau}}{\partial
t^{2}}+\gamma\frac{\partial u_{t}^{\tau}}{\partial t}=f(u_{t}^{\tau}),$

where
$e(u_{t}^{\tau})=\int_{0}^{N}g(u_{t}^{\tau})d\tau-\delta(u_{t}^{N}-u_{t}^{0})$,
and
$h(u_{t}^{\tau})=\int_{0}^{T}f(u_{t}^{\tau})dt-\gamma(u_{T}^{\tau}-u_{0}^{\tau})$.

Proof. Obviously,

$\frac{\partial^{2}u_{t}^{\tau}}{\partial\tau^{2}}+\delta\frac{\partial
u_{t}^{\tau}}{\partial\tau}=g(u_{t}^{\tau}),\quad\frac{\partial^{2}u_{t}^{\tau}}{\partial
t^{2}}+\gamma\frac{\partial u_{t}^{\tau}}{\partial t}=f(u_{t}^{\tau}),$
$so,\quad\frac{\partial
u_{t}^{\tau}}{\partial\tau}+\delta(u_{t}^{N}-u_{t}^{0})=\int_{0}^{N}g(u_{t}^{\tau})d\tau,\quad\frac{\partial
u_{t}^{\tau}}{\partial
t}+\gamma(u_{T}^{\tau}-u_{0}^{\tau})=\int_{0}^{T}f(u_{t}^{\tau})dt.$

Define
$e(u_{t}^{\tau})=\int_{0}^{N}g(u_{t}^{\tau})d\tau-\delta(u_{t}^{N}-u_{t}^{0})$,
and
$h(u_{t}^{\tau})=\int_{0}^{T}f(u_{t}^{\tau})dt-\gamma(u_{T}^{\tau}-u_{0}^{\tau})$,
we have:

$\frac{\partial
u_{t}^{\tau}}{\partial\tau}=e(u_{t}^{\tau}),\quad\frac{\partial
u_{t}^{\tau}}{\partial t}=h(u_{t}^{\tau}),$

thus,

$u_{t}^{\tau}=\int_{0}^{T}h\left(\int_{0}^{N}e(u_{t}^{\tau})d\tau\right)dt=\int_{0}^{N}e\left(\int_{0}^{T}h(u_{t}^{\tau})dt\right)d\tau,$

where $\frac{\partial e(u_{t}^{\tau})}{\partial\tau}=g(u_{t}^{\tau})$ and
$\frac{\partial h(u_{t}^{\tau})}{\partial t}=f(u_{t}^{\tau})$.

##  Appendix D Proof of Proposition 4.4

Proposition 4.4 Let $\mathbf{X}^{n}$ and $\mathbf{Y}^{n}$ be the node
features, generated by Eqn. 3, and $\Delta t\ll 1$. The gradient of the
second-order CGNNs $\mathbf{W}_{l}$ is bounded as Eqn. 21, and the gradient of
the second-order SNNs $\mathbf{W}^{k}$ is bounded as Eqn. 22:

$\displaystyle\left|\frac{\partial\mathcal{L}}{\partial\mathbf{W}_{l}}\right|\leq\frac{\beta^{{}^{\prime}}\hat{\mathbf{D}}\Delta
t(1+\Gamma T\Delta t)}{v}\left(\max\limits_{1\leq i\leq
v}(|\mathbf{X}_{i}^{0}|+|\mathbf{Y}_{i}^{0}|)\right)+\frac{\beta^{{}^{\prime}}\hat{\mathbf{D}}\Delta
t(1+\Gamma T\Delta t)}{v}\left(\max\limits_{1\leq i\leq
v}|\bar{\mathbf{X}_{i}}|+\beta\sqrt{T\Delta t}\right)^{2},$ (21)
$\displaystyle\left|\frac{\partial\mathcal{L}}{\partial\mathbf{W}^{k}}\right|\leq\frac{(1+T\Gamma\Delta
t)(1+N\Theta\Delta\tau)V_{th}}{v\beta^{2}\Delta\tau}\left(\max\limits_{1\leq
i\leq v}|\mathbf{X}_{i}^{N}|+\max\limits_{1\leq i\leq
v}|\bar{\mathbf{X}}_{i}|\right).$ (22)

where $\beta=\max\limits_{x}|\sigma(x)|$,
$\beta^{{}^{\prime}}=\max\limits_{x}|\sigma^{{}^{\prime}}(x)|$,
$\hat{D}=\max\limits_{i,j\in\mathcal{V}}\frac{1}{\sqrt{d_{i}d_{j}}}$, and
$\Gamma:=6+4\beta^{{}^{\prime}}\hat{D}\max\limits_{1\leq n\leq
T}||\mathbf{W}^{n}||_{1}$,
$\Theta:=6+4\beta^{{}^{\prime}}\hat{D}\max\limits_{1\leq n\leq
N}||\mathbf{W}^{n}||_{1}.$ $d_{i}$ is the degree of node $i$,
$\mathbf{\bar{X}}_{i}$ is the label of node $i$.

$\displaystyle\frac{\partial\mathcal{L}}{\partial\bm{W}^{k}}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\bm{Z}_{T}^{N}}\frac{\partial\bm{Z}_{T}^{N}}{\partial\bm{Z}_{l}^{N}}\frac{\partial\bm{Z}_{l}^{N}}{\partial\bm{W}^{k}}=\frac{\partial\mathcal{L}}{\partial\bm{Z}_{T}^{N}}\prod_{n=l+1}^{T}\frac{\partial\bm{Z}_{n}^{N}}{\partial\bm{Z}_{n-1}^{N}}\frac{\partial\bm{Z}_{l}^{N}}{\partial\bm{W}^{k}}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\bm{Z}_{T}^{N}}\prod_{n=l+1}^{T}\frac{\partial\bm{Z}_{n}^{N}}{\partial\bm{Z}_{n-1}^{N}}\frac{\partial\bm{Z}_{l}^{N}}{\partial\bm{Z}_{l}^{k}}\frac{\partial\bm{Z}_{l}^{k}}{\partial\bm{W}^{k}}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\bm{Z}_{T}^{N}}\prod_{n=l+1}^{T}\frac{\partial\bm{Z}_{n}^{N}}{\partial\bm{Z}_{n-1}^{N}}\prod_{i={k+1}}^{N}\frac{\partial\bm{Z}_{l}^{i}}{\partial\bm{Z}_{l}^{i-1}}\frac{\partial\bm{Z}_{l}^{k}}{\partial\bm{W}^{k}},$

From (Rusch et al., 2022), we have:

$\left\lVert\frac{\partial\mathcal{L}}{\partial\bm{Z}_{T}^{N}}\right\rVert_{\infty}\leq\frac{1}{v}\left(\max\limits_{1\leq
i\leq v}|\bm{X}_{i}^{N}|+\max\limits_{1\leq i\leq
v}|\bar{\bm{X}}_{i}|\right),\quad\left\lVert\frac{\partial\bm{Z}_{T}^{N}}{\partial\bm{Z}_{t}^{N}}\right\rVert_{\infty}\leq
1+T\Gamma\Delta t.$ (23)

Due to the second-order SNN has a similar formulation to second-order GNN, we
have a similar conclusion,

$\left\lVert\frac{\partial\bm{Z}_{l}^{N}}{\partial\bm{Z}_{l}^{k}}\right\rVert_{\infty}\leq
1+N\Theta\Delta\tau,$ (24)

with $\beta=\max\limits_{x}|\sigma(x)|$,
$\beta^{{}^{\prime}}=\max\limits_{x}|\sigma^{{}^{\prime}}(x)|$,
$\hat{D}=\max\limits_{i,j\in\mathcal{V}}\frac{1}{\sqrt{d_{i}d_{j}}}$, and
$\Theta:=6+4\beta^{{}^{\prime}}\hat{D}\max\limits_{1\leq n\leq
N}||\bm{W}^{n}||_{1}.$, and with Eq. 12, we have:

$\frac{\partial\bm{Z}_{l}^{k}}{\partial\bm{W}^{k}}\approx
r(\bm{Z}_{l}^{k-1})\leq\frac{V_{th}}{\beta^{2}\Delta\tau}.$ (25)

Multipling 23, 24 and 25, we have the upper bound:

$\frac{\partial\mathcal{L}}{\partial\bm{W}^{k}}\leq\frac{(1+T\Gamma\Delta
t)(1+N\Theta\Delta\tau)V_{th}}{v\beta^{2}\Delta\tau}\left(\max\limits_{1\leq
i\leq v}|\bm{X}_{i}^{N}|+\max\limits_{1\leq i\leq
v}|\bar{\bm{X}}_{i}|\right).$ (26)


# Can the CIB constrain the dark energy?

A. Maniyar1 G. Lagache1 M. Béthermin1 S. Ilić2 1Aix Marseille Univ, CNRS,
CNES, LAM, Marseille, France
2CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 2,
Praha 8 Czech Republic

###### Abstract

Galaxies are often used as tracers of the large scale structure (LSS) to
measure the Integrated Sachs-Wolfe effect (ISW) by cross-correlating the
galaxy survey maps with the Cosmic Microwave Background (CMB) map. We use the
Cosmic Infrared Background (CIB) as a tracer of the LSS to perform a
theoretical CIB-CMB cross-correlation to measure the ISW for different Planck
HFI frequencies. We discuss the detectability of this ISW signal using a
Signal-to-noise ratio analysis and find that the ISW detected this way can
provide us with the highest SNR for a single tracer ranging from 5 to 6.7
(maximum being for 857 GHz) with the CIB and CMB maps extracted over the whole
sky. A Fisher matrix analysis showed that this measurement of the ISW can
improve the constraints on the cosmological parameters; especially the
equation of state of the dark energy $w$ by $\sim 47\%$. Performing a more
realistic analysis including the galactic dust residuals in the CIB maps over
realistic sky fractions shows that the dust power spectra dominate over the
CIB power spectra at $\ell<100$ and ISW can’t be detected with high SNR. We
perform the cross-correlation on the existing CIB-CMB maps over $\sim 11\%$ of
the sky in the southern hemisphere and find that the ISW is not detected with
the existing CIB maps over such small sky fractions.

##  1 Integrated Sachs-Wolfe effect and Cosmic Infrared Background

ISW effect is a potential powerful probe of the dark energy. It’s the low-
redshift counterpart of the Sachs-Wolfe effect happening at the last
scattering surface. The ISW occurs in an Universe not dominated by matter. As
the dark energy starts dominating, it causes the gravitational potential wells
and hills on large scales to decay. As a result, the CMB photons travelling
across them undergo a net gain or loss of energy for potential wells and hills
respectively. It is the dominant contribution to the CMB power spectrum at
large angular scales and has a very small amplitude. Unfortunately, at these
scales, the statistical noise due to cosmic variance is of the same order of
magnitude as the signal making the ISW detection in the CMB alone extremely
challenging. A non-zero spatial correlation between the CMB anisotropies and
tracers of the large scale matter distribution is expected, and can be used to
extract the ISW signal. In this study, we follow the approach of Ilić et
al(2014)[1] and use the CIB as a LSS tracer for this purpose. The CIB is the
weighted integral of the dust heated by the young UV-bright stars within the
galaxies through the cosmic time. Anisotropies in the CIB trace the large
scale distribution of the galaxy density field and hence, the underlying
distribution of the dark matter haloes which host these galaxies, up to a bias
factor. Therefore, the CIB can be used as a tracer of the LSS.

###  1.1 CIB, CIB-CMB lensing modeling

The ISW dominates at the lower multipoles (large angular scales) where the
clustering is dominated by the correlation between the dark matter haloes and
non-linear terms can be ignored. We use the linear CIB model presented in
Maniyar et al (2018)[2]. The linear power spectrum is:

$C_{l}^{\nu\times\nu^{\prime}}=\int\frac{dz}{\chi^{2}}\frac{d\chi}{dz}a^{2}b_{eff}^{2}\bar{j}(\nu,z)\bar{j}(\nu^{\prime},z)P_{lin}(k=l/\chi,z)$
(1)

where $P_{lin}(k,z)$ is the linear theory dark matter power spectrum,
$\chi(z)$ is the comoving distance to redshift $z$ and $a=1/(1+z)$ is the
scale factor. $b_{eff}$ is the mean bias of dark matter haloes hosting dusty
galaxies contributing to the CIB at a given redshift weighted by their
contribution to the emissivities (Planck Collaboration (2014)[3] and Maniyar
et al (2018)[2]) and is called the effective bias. $\bar{j}$ is the mean
comoving emissivity of the CIB galaxies. It can be seen from Eq. 1 that
$b_{eff}$ and $\bar{j}$ are degenerate. To break this degeneracy, we utilise
the CIB-CMB lensing cross-correlation which is given by

$C_{l}^{\nu\phi}=\int
b_{eff}\bar{j}(\nu,z)\frac{3}{l^{2}}\Omega_{m}H_{0}^{2}\bigg{(}\frac{\chi_{*}-\chi}{\chi_{*}\chi}\bigg{)}P_{lin}(k=l/\chi,z)d\chi\,,$
(2)

where $\chi_{*}$ is the comoving distance to the CMB last scattering surface,
$\Omega_{m}$ is the cosmological matter density and $H_{0}$ is the value of
the Hubble’s constant today. From Eq. 2, we see that $C_{l}^{\nu,\phi}$ is
proportional to $b_{eff}$, whereas, $C_{l}^{CIB}$ is proportional to
$b_{eff}^{2}$. Thus, using also the CIB-CMB lensing potential correlation in
the likelihood helps us partially resolve the degeneracy between $b_{eff}$ and
$\bar{j}$ parameters. Details of all the parameters of the model and various
data-sets used to fit the model are given in Maniyar et al (2018)[2] and
Maniyar et al (in prep)[4].

##  2 CIB-ISW cross-correlation and cosmological constraints

The SNR for the ISW obtained through the CIB-CMB cross-correlation is given
as:

$\Big{[}\frac{S}{N}\Big{]}^{2}(\nu)=\sum_{\ell=2}^{\ell_{max}}(2\ell+1)\frac{f_{\mathrm{sky}}\,\big{[}C_{\ell}^{\mathrm{CIBxCMB}}(\nu)\big{]}^{2}}{\big{[}C_{\ell}^{\mathrm{CIBxCMB}}(\nu)\big{]}^{2}+\big{[}C_{\ell}^{\mathrm{CIB}}(\nu)+N_{\ell}^{\mathrm{CIB}}(\nu)\big{]}\,C_{\ell}^{\mathrm{CMB}}}$
(3)

where the N${}_{\ell}^{\mathrm{CIB}}(\nu)$ term contains the galactic dust
residuals left in the CIB maps and fsky represents the fraction of the sky
common to the CMB and the CIB maps.

Figure 1: Left panel shows the predicted cumulative SNR for different Planck
frequencies in the ideal case. The right panel shows the power spectra of the
CIB and the galactic dust at 353 GHz and 3000 GHz with black and red lines
respectively. The dust power spectra have been calculated over 20, and 80% of
the sky.

###  2.1 Ideal case

The ideal case assumes that the CIB and the CMB maps are extracted over the
full sky and are completely dust-free i.e. N${}_{\ell}^{\mathrm{CIB}}(\nu)=0$
and f${}_{\mathrm{sky}}=1$. The only limiting factor is the cosmic variance
and we get the highest SNR in this scenario.
Left panel of Fig. 1 shows the cumulative SNR for the ISW at each frequency.
We observe that the multipoles up to $\ell\leq 50$ contribute the most to the
SNR and it could reach as high as 6.7 for 857 GHz. With a signal observed at
such high SNR, the CIB-ISW cross-correlation can be a competitive probe of the
cosmology and the dark energy. We therefore perform a joint Fisher matrix
analysis with the CMB, the CIB-CMB lensing and CIB-ISW cross-correlation to
compute the improvement in the constraints on the cosmological parameters.
Combining the information from the CIB and CMB alone in a Fisher matrix
provides with constraints on the cosmological parameters. Adding the
information from the CIB-CMB lensing cross-correlation as well as the ISW
obtained through the CIB-CMB cross-correlation further improves these
constraints.

Table 1: The predicted improvement in the constraints on the cosmological
parameters after adding the CIB$\times$Lensing and the ISW signal in the CIB +
CMB Fisher matrix for 100% sky for $\Lambda$CDM model. % improvement on the
CIB+CMB constraints with the ISW over 100% sky
---
Parameters | H0 | $\Omega_{b}h^{2}$ | $\Omega_{c}h^{2}$ | $\tau$ | ns | 109As
[CIB + CMB] + CIB$\times$Lensing | 1.62 | 0.79 | 1.49 | 4.58 | 0.58 | 5.12
[CIB + CMB + CIB$\times$Lensing] + ISW | 6.86 | 2.32 | 8.77 | 16.9 | 4.77 | 16.29
Table 2: The predicted improvement in the constraints on the cosmological
parameters after adding the CIB$\times$Lensing and the ISW signal to the
CIB+CMB+BAO+SNe Ia+H0 prior for 100% sky for the $w$CDM model. % improvement
on the CIB+CMB+Ext constraints with the ISW over 100% sky
---
Parameters | H0 | $w$ | $\Omega_{b}h^{2}$ | $\Omega_{c}h^{2}$ | $\tau$ | ns | 109As
[CIB + CMB + Ext] + CIB$\times$Lensing | 3.83 | 2.60 | 0.33 | 0.52 | 3.67 | 0.24 | 4.20
[CIB + CMB + Ext + CIB$\times$Lensing] + ISW | 31.06 | 49.17 | 5.13 | 26.35 | 22.64 | 12.41 | 21.11

Table 1 shows this improvement in the context of the $\Lambda$CDM model. We
observe an improvement of 6% on the constraints on $\Omega_{m}$ (and hence
$\Omega_{\Lambda}$) and of 11% on $\tau$ (and hence $10^{9}A_{s}$ as they are
positively correlated) after adding the ISW information alone on top of the
CIB + CMB + CIB-CMB lensing cross-correlation information.
In a flat Universe, the ISW directly results from the dark energy and it has a
strong dependence on the equation of state of the dark energy $w=P/\rho$. In
the $\Lambda$CDM model, $w$ is kept constant at -1. In the $w$CDM model, $w$
is a free parameter along with the 6 other cosmological parameters of the
$\Lambda$CDM model. Due to the geometrical degeneracy between $w$ and H0, CMB
data alone can not constrain both these parameters together. Therefore,
external data sets like baryonic acoustic oscillations (BAO) and supernovae
type Ia (SNe Ia), with the low redshift distance measurements are used with
the CMB data to break this degeneracy. As in the case of the $\Lambda$CDM
model, Tab. 2 shows the improvement in the constraints on the cosmological
parameters in the context of the $w$CDM model. We observe an improvement of
47% on the constraints on $w$ and of 18% on $\tau$ and $10^{9}A_{s}$ after
adding the ISW information alone on top of the CIB + CMB + CIB-CMB lensing
cross-correlation + Ext (BAO + SNE Ia + H0 prior) data.

###  2.2 Realistic case

On a large part of the sky, we have strong emission from our own Galaxy. This
emission is much higher than the CIB signal. This prevents us to extract the
CIB signal over the part of the sky where the galactic plane resides and
reduces the available sky fraction for measurement by 30%. Even in the
remaining part of the sky, the CIB maps are contaminated by dust that needs to
be removed. This removal leaves residual noise and has to be accounted for in
the analysis. We calculate the dust power spectra at different frequencies for
different sky fractions. Right panel of Fig. 1 shows the CIB and the galactic
dust power spectra at 353 and 3000 GHz for different sky fractions. The dust
power spectra dominate over the CIB power spectra at lower multipoles where we
expect to extract the ISW signal. As mentioned in Planck Collaboration
(2014)[3], the dust residuals are of the order of 5-10% on the power spectrum.
In our analysis, we take the conservative 10% dust residuals of
C${}_{\ell}^{\mathrm{dust}}(\nu)$ as the noise term. From Eq. 3, we see that
the N${}_{\ell}^{\mathrm{CIB}}(\nu)$ term is quite big compared to the ideal
noiseless ideal case and the SNR drops quite drastically with values from 0.5
for 3000 GHz to 1.5 for 353 GHz for 20% of the sky.

##  3 Conclusions

Using the CIB as a tracer of the LSS, we calculate the theoretical prediction
for the cross-correlation of the CIB with the CMB to extract an ISW signal. We
find that the ISW could be obtained with the highest SNR (reaching up to 6.7
for 857 GHz) for a single tracer with this technique when the CIB and CMB maps
are dust-free and are extracted over the whole sky. We show that the ISW
detected at such high SNR can improve the constraints on the cosmological
parameters quite considerably; especially the constraints on the equation of
state of the dark energy are improved by $\sim 47\%$. However, in the
realistic case, we are limited by the presence of the galactic dust residuals
in the maps and also by the available sky fractions over which the maps are
extracted. When the same analysis is performed with the addition of 10%
galactic dust residual power spectra, they completely dominate over the CIB at
lower multipoles of interest and obtaining an ISW signal with high SNR is not
possible anymore. In the end, we also performed the cross-correlation of the
CIB and the CMB Planck maps over $\sim 11\%$ of the southern sky in a field
named ’GASS’ (Planck Collaboration (2014)[3]). The CIB map in the GASS field
is the cleanest CIB map currently published. We find that with the cross-
correlation performed over such a small sky fraction and with the CIB maps
containing the dust residuals, we don’t detect the ISW signal. In future
though, with the availability of the CIB maps with better galactic dust
removal, this method will be very powerful to detect the ISW with a high SNR
and thereby constraining the dark energy parameters $\Omega_{\Lambda}$ and
$w$.

## Acknowledgements

We acknowledge financial support from the “Programme National de Cosmologie
and Galaxies” (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France, from
the ANR under the contract ANR-15-CE31-0017 and from the OCEVU Labex
(ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the
“Investissements d’Avenir” French government programme managed by the ANR.

## References

## References

  * [1] Ilić, S., Douspis, M., Langer, M., Pénin, A., Lagache, G. 2014, A&A, 571, A30.
  * [2] Maniyar, A., Béthermin, M., & Lagache, G. 2018, ArXiv e-prints, [arXiv:1801.10146].
  * [3] Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A30.
  * [4] Maniyar, A., Lagache, G., Béthermin, M., & Ilić, S. 2018, A&A (to be submitted).