$\boldsymbol{E}_{G}$, $\boldsymbol{F}_{G}$ and $\boldsymbol{M}_{G}$. By [32,
Theorem 56], Item 3 is equivalent to the existence of a doubly-stochastic map
$X\colon\mathbb{C}\mathcal{T}^{\circ}_{G}\to\mathbb{C}\mathcal{T}^{\circ}_{H}$
such that $X\boldsymbol{B}_{G}=\boldsymbol{B}_{H}X$ for all
$\boldsymbol{B}\in\\{\boldsymbol{I},\boldsymbol{C},\boldsymbol{G},\boldsymbol{E},\boldsymbol{M}\\}$.
This is the content of Item 4. ∎

##  Appendix C A Comonadic Strategy for Homomorphism Indistinguishability

We propose to view our Theorem 1.5 as the instantiation of a novel strategy
for proving characterisations of homomorphism indistinguishability in terms of
matrix equations for graph classes which so far have resisted the strategy
outlined in Section 1.2. This novel strategy is inspired by the comonadic
framework developed in [2, 4, 15, 3], etc. It is in large parts generic, in
the sense that it does not require prior knowledge of the graph class in
questions. This is achieved by leveraging the properties of a comonad
characterising the graph class.

###  C.1 Preliminaries

We first sketch the framework of comonads on the category of graphs before
commenting on how it can be used to derive matrix equations. See [15, 38] for
further details. Let $\mathsf{Graph}$ denote the category of finite graphs,
i.e. with objects being finite graphs and morphisms being graph homomorphisms.
A _comonad_ on $\mathsf{Graph}$ is a tuple $(\mathfrak{C},\varepsilon,\delta)$
where $\mathfrak{C}\colon\mathsf{Graph}\to\mathsf{Graph}$ is a functor and
$\varepsilon\colon\mathfrak{C}\to 1_{\mathsf{Graph}}$ and
$\delta\colon\mathfrak{C}\to\mathfrak{C}\mathfrak{C}$ are natural
transformations satisfying the following diagrams for all objects $A$ of
$\mathsf{Graph}$:

$Y_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{2}$${\mathfrak{C}A}$${\mathfrak{C}\mathfrak{C}A}$${\mathfrak{C}\mathfrak{C}A}$${\mathfrak{C}\mathfrak{C}\mathfrak{C}A}$$\scriptstyle{\delta_{A}}$$\scriptstyle{\delta_{A}}$$\scriptstyle{\mathfrak{C}\delta_{A}}$$\scriptstyle{\delta_{\mathfrak{C}A}}$
$Y_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{2}$${\mathfrak{C}A}$${\mathfrak{C}\mathfrak{C}A}$${\mathfrak{C}\mathfrak{C}A}$${\mathfrak{C}A}$$\scriptstyle{\delta_{A}}$$\scriptstyle{\delta_{A}}$$\scriptstyle{\operatorname{id}_{\mathfrak{C}A}}$$\scriptstyle{\mathfrak{C}\varepsilon_{A}}$$\scriptstyle{\varepsilon_{\mathfrak{C}A}}$

######  Example C.1 ([2, 15]).

The _pebbling comonad_ maps a graph $A$ to the graph $\mathfrak{P}_{k,d}A$
with vertex set being the set of sequences
$[(p_{1},a_{1}),\dots,(p_{\ell},a_{\ell})]$ for $p_{i}\in[k]$, $a_{i}\in
V(A)$, and $i\in[\ell]$, $\ell\leq d$. There is an edge between
$[(p_{1},a_{1}),\dots,(p_{\ell},a_{\ell})]$ and
$[(p^{\prime}_{1},a^{\prime}_{1}),\dots,(p^{\prime}_{\ell^{\prime}},a^{\prime}_{\ell^{\prime}})]$
if and only if one sequence is and initial segment of the other, and if wlog
the first sequence is contained in the second then
$p_{\ell}\not\in\\{p^{\prime}_{\ell+1},\dots,p^{\prime}_{\ell^{\prime}}\\}$,
and $a_{\ell}a^{\prime}_{\ell^{\prime}}\in E(A)$. The homomorphisms
$\varepsilon_{A}\colon\mathfrak{P}_{k,d}A\to A$ sends
$[(p_{1},a_{1}),\dots,(p_{\ell},a_{\ell})]$ to $a_{\ell}$. It is readily
verified that this is a morphisms in $\mathsf{Graph}$. The homomorphisms
$\delta_{A}\colon\mathfrak{P}_{k,d}A\to\mathfrak{P}_{k,d}\mathfrak{P}_{k,d}A$
maps

$[(p_{1},a_{1}),\dots,(p_{\ell},a_{\ell})]\mapsto[(p_{1},s_{1}),\dots,(p_{\ell},s_{\ell})]$

where $s_{j}=[(p_{1},a_{1}),\dots,(p_{j},a_{j})]$ for all $j\leq\ell$. The
triple $(\mathfrak{P}_{k,d},\varepsilon,\delta)$ specifies the pebbling
comonad.

A _coalgebra_ for a comonad $(\mathfrak{C},\varepsilon,\delta)$ is a pair
$(A,\alpha)$ where $A\in\operatorname{obj}\mathsf{Graph}$ and $\alpha\colon
A\to\mathfrak{C}A$ is a morphism such that the following diagrams commute:

$Y_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{2}$${A}$${\mathfrak{C}A}$${\mathfrak{C}A}$${\mathfrak{C}\mathfrak{C}A}$$\scriptstyle{\alpha}$$\scriptstyle{\alpha}$$\scriptstyle{\mathfrak{C}\alpha}$$\scriptstyle{\delta_{A}}$
$Y_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{2}$${A}$${\mathfrak{C}A}$${A}$$\scriptstyle{\alpha}$$\scriptstyle{\operatorname{id}_{A}}$$\scriptstyle{\varepsilon_{A}}$

The significance of coalgebras for this paper stems from the following result:

######  Theorem C.2 ([15, Section IV.B]).

Let $A\in\operatorname{obj}\mathsf{Graph}$. Then the following are equivalent:

  1. 1.

$A$ admits a $k$-pebble forest cover of depth $\leq d$,

  2. 2.

$A$ admits a $\mathfrak{P}_{k,d}$-coalgebra $A\to\mathfrak{P}_{k,d}A$.

In fact, there is a bijective correspondence between these objects. We
therefore think of coalgebras as graph decompositions.

The _Eilenberg–Moore category_ $\operatorname{EM}(\mathfrak{C})$ of a comonad
$\mathfrak{C}$ is the category whose objects are $\mathfrak{C}$-coalgebras and
whose morphisms are $\mathfrak{C}$-coalgebra morphisms. A _morphism
$h\colon(A,\alpha)\to(B,\beta)$ of $\mathfrak{C}$-coalgebras_ is a morphisms
$h\colon A\to B$ such that the following diagram commutes:

$Y_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{2}$${A}$${\mathfrak{C}A}$${B}$${\mathfrak{C}B}$$\scriptstyle{h}$$\scriptstyle{\alpha}$$\scriptstyle{\mathfrak{C}h}$$\scriptstyle{\beta}$

Informally, the Eilenberg–Moore category is the category of decomposed graphs.
The _forgetful functor_
$U^{\mathfrak{C}}\colon\operatorname{EM}(\mathfrak{C})\to\mathsf{Graph}$ which
sends $(A,\alpha)$ to $A$ associates a decomposed graph to its underlying
graph. Thus, by Theorem C.2, a graph is in the image of
$U^{\mathfrak{P}_{k,d}}$ if and only if it is admits a $k$-pebble forest cover
of depth $\leq d$.

######  Fact C.3 (cf. e.g. [38]).

Let $\mathfrak{C}$ be a comonad over a category $\mathcal{A}$. Then
$\operatorname{EM}(\mathfrak{C})$ has all colimits that $\mathcal{A}$ has. The
forgetful functor
$U^{\mathfrak{C}}\colon\operatorname{EM}(\mathfrak{C})\to\mathcal{A}$ creates
them.

From the perspective of homomorphism counting, the characterisation of the
adjunction underlying the comonad $\mathfrak{C}$ in terms of homomorphism
functor is of utmost importance. Let
$F^{\mathfrak{C}}\colon\mathsf{Graph}\to\operatorname{EM}(\mathfrak{C})$
denote the _free functor_ mapping $A$ to $(\mathfrak{C}A,\delta_{A})$, i.e.
$F^{\mathfrak{C}}\vdash U^{\mathfrak{C}}$.

######  Fact C.4 ([49, Corollary 5.4.23] or e.g. [15]).

Let $(A,\alpha)\in\operatorname{obj}\operatorname{EM}(\mathfrak{C})$ and
$G\in\operatorname{obj}\mathsf{Graph}$. Then

$\hom_{\mathsf{Graph}}(U^{\mathfrak{C}}(A,\alpha),G)\cong\hom_{\operatorname{EM}(\mathfrak{C})}((A,\alpha),F^{\mathfrak{C}}G).$

Denote by $\widehat{\ell}\colon(A,\alpha)\to F^{\mathfrak{C}}G$ the morphism
corresponding to $\ell\colon U^{\mathfrak{C}}(A,\alpha)\to G$ under this
bijection.

In other words, counting homomorphisms in $\mathsf{Graph}$ from $A$ to $G$
where $A$ admits a $\mathfrak{C}$-coalgebra amounts to counting homomorphisms
in $\operatorname{EM}(\mathfrak{C})$ from a $\mathfrak{C}$-coalgebra of $A$ to
the free $\mathfrak{C}$-coalgebra $(\mathfrak{C}G,\delta_{G})$ of $G$.

###  C.2 Bilabelled Objects and their Coalgebraic and Augmented Homomorphism
Representation

We develop a categorical language for bilabelled objects. The following
definitions resemble those in [43]. Fix a category $\mathcal{C}$, which has
all finite pushouts. An example for such a category is $\mathsf{Graph}$ but
also $\operatorname{EM}(\mathfrak{C})$ for a comonad $\mathfrak{C}$ on
$\mathsf{Graph}$, cf. Fact C.3.

######  Definition C.5 (Categories of Bilabelled Objects).

For $L\in\operatorname{obj}\mathcal{C}$, define
$\operatorname{Labelled}_{\mathcal{C}}(L,K)$ to be the category with objects
being pairs of $\mathcal{C}$-morphisms $L\to A\leftarrow K$ and morphisms
being $\mathcal{C}$-morphisms $h$ such that the following diagram commutes:

$Y_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{2}$${L}$${A}$${K}$${B}$$\scriptstyle{h}$

For a family of objects $\mathcal{L}\subseteq\operatorname{obj}\mathcal{L}$,
write $\operatorname{Labelled}_{\mathcal{C}}(\mathcal{L})$ for the union of
all $\operatorname{Labelled}_{\mathcal{C}}(L,K)$ for $L,K\in\mathcal{L}$. We
use boldface letters $\boldsymbol{A}=(L\to A\leftarrow K)$ to indicate
bilabelled objects.

An instance of a category of bilabelled objects can be constructed from a
comonad $\mathfrak{C}$ on $\mathsf{Graph}$ and its Eilenberg–Moore category
$\operatorname{EM}(\mathfrak{C})$ as follows. Fix a graph
$L\in\operatorname{obj}\mathsf{Graph}$. Write $\mathcal{L}$ for the set of all
$\mathfrak{C}$-coalgebras $L\to\mathfrak{C}L$ of $L$. Let
$\Lambda(\mathfrak{C},L)\coloneqq\operatorname{obj}\operatorname{Labelled}_{\operatorname{EM}(\mathfrak{C})}(\mathcal{L})\sqcup\\{\bot\\}$
where $\bot$ is a fresh symbol. In other words, $\Lambda(\mathfrak{C},L)$
comprises all bilabelled objects of the form
$(L,\lambda)\to(A,\alpha)\leftarrow(L,\lambda^{\prime})$ and $\bot$. The set
$\Lambda(\mathfrak{C},L)$ is the family of bilabelled graphs of our interest.

######  Example C.6.

With $L$ the edge-less graph with $d$ vertices,
$\Lambda(\mathfrak{P}_{k+1,d},L)$ resembles $\mathcal{WL}_{k}^{d}$ as
introduced in Definition 4.1. Differences arise since Definition 4.1 includes
only certain pebble forest covers arising in Theorem 2.2.

We proceed to define a suitable homomorphism representation. First, we
generalise the homomorphism representation of [44, 33] which does not encode
any information about the labels. To avoid confusion with homomorphism counts,
write $\operatorname{mor}(L,G)$ for the set of morphisms $L\to G$ in
$\mathsf{Graph}$. Note that since we consider finite graphs only, this set is
finite. Recall the notation $\widehat{\ell}$, $\widehat{k}$, etc. from Fact
C.4.

######  Definition C.7.

Let $G$ be a graph. The _comonadic homomorphism representation_ is the map
$\Lambda(\mathfrak{C},L)\to\mathbb{C}^{\operatorname{mor}(L,G)\times\operatorname{mor}(K,G)}$
mapping
$\boldsymbol{A}=\left((L,\lambda)\to(A,\alpha)\leftarrow(K,\kappa)\right)$ to
the matrix
$\boldsymbol{A}_{G}\in\mathbb{C}^{\operatorname{mor}(L,G)\times\operatorname{mor}(K,G)}$
with entry $(\ell,k)$ being the number of homomorphisms $A\to G$ such that the
diagram

$Y_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{2}$${(L,\lambda)}$${(A,\alpha)}$${(K,\kappa)}$${F^{\mathfrak{C}}G}$$\scriptstyle{\hat{\ell}}$$\scriptstyle{\hat{k}}$

commutes. Furthermore, $\bot$ is mapped to the zero matrix.

It remains to define a representation which also incorporates information
about the structure of the labels, i.e. their coalgebras.

######  Definition C.8.

Let $G$ be a graph. The _augmented homomorphism representation_ associates
with every element of $\Lambda(\mathfrak{C},L)$ an endomorphism of the
$\mathbb{C}$-vector space

$\mathbb{C}^{\operatorname{mor}(L,G)}\otimes\bigoplus_{\begin{subarray}{c}\mathfrak{C}\text{-coalgebras}\\\
\lambda\colon L\to\mathfrak{C}L\end{subarray}}\mathbb{C}.$

It maps
$\boldsymbol{A}=\left((L,\lambda)\to(A,\alpha)\leftarrow(L,\kappa)\right)$ to
$\widehat{\boldsymbol{A}}_{G}\coloneqq\boldsymbol{A}_{G}\otimes
e_{\lambda}e_{\kappa}^{T}$ where $\boldsymbol{A}_{G}$ is defined in Definition
C.7 and $e_{\lambda},e_{\kappa}\in\bigoplus_{\lambda\colon
L\to\mathfrak{C}L}\mathbb{C}$ are the $\lambda$-th and $\kappa$-th standard
basis vectors respectively. Furthermore, $\bot$ is mapped to the zero matrix.

Note that we do not succeed to parallel the notation of Section 4 entirely. In
Section 4, $\boldsymbol{A}$ and $\widehat{\boldsymbol{A}}$ indicated different
objects with different representations $\boldsymbol{A}_{G}$ and
$\widehat{\boldsymbol{A}}_{G}$. Here, the underlying objects are the same,
only the representations are different. We write
$\boldsymbol{A}=((L,\lambda)\to(A,\alpha)\leftarrow(K,\kappa))$ to indicate
bilabelled coalgebras, $\boldsymbol{A}_{G}$ to indicate the comonadic
homomorphism representation, and $\widehat{\boldsymbol{A}}_{G}$ to indicate
the augmented homomorphism representation.

###  C.3 Operations

The familiar operations on bilabelled objects can also be defined in
categorical terms.

######  Definition C.9 (Operations on Bilabelled Objects).

  1. 1.

Define the _series composition_ of $L\to A\leftarrow K$ and $K\to B\leftarrow
M$, in symbols $(L\to A\leftarrow K)\cdot(K\to B\leftarrow M)$, as the unique
morphisms $L\to C\leftarrow M$ induced by the pushout diagram in Equation 3
where $C$ is the pushout.

$Y_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{2}$${L}$${K}$${M}$${A}$${B}$${C}$
(3)

  2. 2.

Define the _reverse_ of $L\to A\leftarrow K$, in symbols $(L\to A\leftarrow
K)^{*}$, as $K\to A\leftarrow L$.

The operations in Definition C.9 are compatible with the comonadic
homomorphism representation.

######  Lemma C.10.

Let $\boldsymbol{A}=\left((L,\lambda)\to(A,\alpha)\leftarrow(K,\kappa)\right)$
and $\boldsymbol{B}=\left((K,\kappa)\to(B,\beta)\leftarrow(M,\mu)\right)$. Let
$G$ be a graph. Then
$(\boldsymbol{A}\cdot\boldsymbol{B})_{G}=\boldsymbol{A}_{G}\cdot\boldsymbol{B}_{G}$.

###### Proof.

Write $(C,\gamma)$ for the pushout of $a_{2}\colon(K,\kappa)\to(A,\alpha)$ and
$b_{1}\colon(K,\kappa)\to(B,\beta)$. Let $\ell\in\operatorname{mor}(L,G)$ and
$m\in\operatorname{mor}(M,G)$. Consider the following diagram.

$Y_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{1}$$X_{1}$$X_{2}$$Y_{2}$${{(L,\lambda)}}$${{(K,\kappa)}}$${{(M,\mu)}}$${{(A,\alpha)}}$${{(B,\beta)}}$${{(C,\gamma)}}$${F^{\mathfrak{C}}G}$$\scriptstyle{a_{1}}$$\scriptstyle{\widehat{\ell}}$$\scriptstyle{a_{2}}$$\scriptstyle{b_{1}}$$\scriptstyle{b_{2}}$$\scriptstyle{\widehat{m}}$

By the definition of the pushout and Fact C.4, the set of morphisms
$(C,\gamma)\to F^{\mathfrak{C}}G$ making the above diagram commute is in
bijection with the set of triples $(k,a,b)$ where
$k\in\operatorname{mor}(K,G)$, $a\colon(A,\alpha)\to F^{\mathfrak{C}}G$,
$b\colon(B,\beta)\to F^{\mathfrak{C}}G$, such that
$aa_{2}=\widehat{k}=bb_{1}$, $\widehat{\ell}=aa_{1}$, and
$\widehat{m}=bb_{2}$. This implies the claim. ∎

Lemma C.10 has immediate consequences for the augmented homomorphism
representation.

######  Corollary C.11.

Let $\boldsymbol{A}=\left((L,\lambda)\to(A,\alpha)\leftarrow(K,\kappa)\right)$
and $\boldsymbol{B}=\left((K,\kappa)\to(B,\beta)\leftarrow(M,\mu)\right)$. Let
$G$ be a graph. Then
$(\widehat{\boldsymbol{A}\cdot\boldsymbol{B}})_{G}=\widehat{\boldsymbol{A}}_{G}\cdot\widehat{\boldsymbol{B}}_{G}$.

With series composition and reversal, $\Lambda(\mathfrak{C},L)$ has the
structure of an _involution semigroup_ , i.e. series composition is an
associative binary operation and reversal induces an involution map
${}^{*}\colon\Lambda(\mathfrak{C},L)\to\Lambda(\mathfrak{C},L)$ such that
$(g\cdot h)^{*}=h^{*}\cdot g^{*}$ and $(g^{*})^{*}=g$ for all
$g,h\in\Lambda(\mathfrak{C},L)$. The series composition of
$(L,\lambda_{1})\to(A,\alpha)\leftarrow(L,\lambda_{2})$ and
$(L,\lambda_{3})\to(A,\alpha)\leftarrow(L,\lambda_{4})$ is defined via
Definition C.9 if $\lambda_{2}=\lambda_{3}$. Otherwise, the result is defined
to be $\bot$. Also the series composition of any element with $\bot$ on either
side is $\bot$. The involution operation is defined as the reversal of
Definition C.9. Again, $(\bot)^{*}\coloneqq\bot$.

As before, the unlabelling operation permits us to pass from a bilabelled
coalgebra to the underlying graph.

######  Definition C.12.

The _unlabelling_ of
$\boldsymbol{A}=\left((L,\lambda)\to(A,\alpha)\leftarrow(K,\kappa)\right)$, in
symbols $\operatorname{soe}\boldsymbol{A}$, is $A$.

Unlabelling corresponds to sum-of-entries under the comonadic and the
augmented homomorphism representation.

######  Lemma C.13.

For every graph $G$ and every
$\boldsymbol{A}=\left((L,\lambda)\to(A,\alpha)\leftarrow(K,\kappa)\right)$,

$\hom_{\mathsf{Graph}}(\operatorname{soe}\boldsymbol{A},G)=\operatorname{soe}\boldsymbol{A}_{G}=\operatorname{soe}\widehat{\boldsymbol{A}}_{G}.$

###### Proof.

By Fact C.4,
$\hom_{\mathsf{Graph}}(A,G)=\hom_{\operatorname{EM}(\mathfrak{C})}((A,\alpha),F^{\mathfrak{C}}G)$.
By Definition C.7, this is the sum-of-entries of $\boldsymbol{A}_{G}$. Recall
from Definition C.8, that only one entry in the second factor of the Kronecker
product $\widehat{\boldsymbol{A}}_{G}$ is non-zero. This entry equals one and
hence
$\operatorname{soe}\boldsymbol{A}_{G}=\operatorname{soe}\widehat{\boldsymbol{A}}_{G}$.
∎

###  C.4 Finite Generation

For the case of the pebbling comonad, finite generation was proven in
Proposition 4.7 after devising a suitable set of generators. Proving finite
generation for the involution semigroup of bilabelled objects over arbitrary
comonads on $\mathsf{Graph}$, remains a problem for further investigations.

######  Definition C.14.

A subset $\Gamma\subseteq\Lambda(\mathfrak{C},L)$ is a _set of generators_ for
$\Lambda(\mathfrak{C},L)$ if for every
$\widehat{\boldsymbol{A}}\in\Lambda(\mathfrak{C},L)$ there exist finitely many
$\widehat{\boldsymbol{B}}^{1},\dots,\widehat{\boldsymbol{B}}^{r}\in\Gamma$
such that
$\widehat{\boldsymbol{A}}=\widehat{\boldsymbol{B}}^{1}\cdot\dots\cdot\widehat{\boldsymbol{B}}^{r}$.

###  C.5 Matrix Equations

It remains to deduce matrix equations characterising homomorphism
indistinguishability over the set of graphs underlying
$\Lambda(\mathfrak{C},L)$. Note that finite generation is not among the
assumptions of the following theorem. It is only needed to obtain a finite
system of matrix equations.

######  Theorem C.15.

Let $\mathfrak{C}$ be a comonad on $\mathsf{Graph}$. Let
$L\in\operatorname{im}U^{\mathfrak{C}}$ be a graph such that
$\Lambda(\mathfrak{C},L)$ admits a set of generators $\Gamma$, which is closed
under reversal. Let $G$ and $H$ be graphs. Then the following are equivalent:

  1. 1.

$G$ and $H$ are homomorphism indistinguishable over the graphs
$\operatorname{soe}\boldsymbol{A}$ for
$\boldsymbol{A}\in\Lambda(\mathfrak{C},L)$.

  2. 2.

There exists a pseudo-stochastic matrix $X$ such that
$X\widehat{\boldsymbol{A}}_{G}=\widehat{\boldsymbol{A}}_{H}X$ for all
$\widehat{\boldsymbol{A}}\in\Gamma$.

###### Proof.

The backward direction is easy to see, cf. proof of Theorem 1.5. The forward
direction follows from Theorem 2.4 observing that
$\operatorname{soe}\widehat{\boldsymbol{A}}_{G}=\operatorname{soe}\boldsymbol{A}_{G}$
for all $\boldsymbol{A}$ and $G$. ∎

###  C.6 The Strategy

We conclude by summarising the comonadic strategy to derive homomorphism
indistinguishability characterisations in terms of matrix equations. This
approach is exemplified with less notational overhead in Section 4 for the
case of the pebbling comonad $\mathfrak{P}_{k,d}$.

  1. 1.

Define a family of labelled graphs and their homomorphism matrices

Fix a graph $L$ and define $\Lambda(\mathfrak{C},L)$ as in Section C.2 as the
family of bilabelled graphs of interest. They are of the form
$(L,\lambda)\to(A,\alpha)\leftarrow(L,\lambda^{\prime})$ where $\lambda$ and
$\lambda^{\prime}$ are coalgebras of $L$ and $A$ is a graph of interest. The
representation in terms of matrices is defined in Definition C.8. It encodes
the homomorphism counts and information about the labelling, i.e. the
coalgebras $\lambda$ and $\lambda^{\prime}$.

For the pebbling comonad, Section 4.1 features the set $\mathcal{WL}_{k}^{d}$.
It comprises tuples $(\boldsymbol{F},p_{\textup{in}},p_{\textup{out}})$ where
$\boldsymbol{F}$ is a bilabelled graph possessing a certain pebble forest
cover. This resembles the above definition. The functions $p_{\textup{in}}$
and $p_{\textup{out}}$ play the role of $\lambda$ and $\lambda^{\prime}$. The
object $L$ is the graph with $d$ isolated vertices placed on a path in its
forest cover.

  2. 2.

Define suitable combinatorial and algebraic operations

Series composition and reversal are the relevant combinatorial operations, cf.
Section C.3. They correspond under the augmented homomorphism representation
to matrix product and transposition. In the coalgebraic case, closure under
these operations follows essentially from basic properties of the
Eilenberg–Moore category $\operatorname{EM}(\mathfrak{C})$. Remembering the
coalgebras $\lambda$ and $\lambda^{\prime}$ also on the algebraic side ensures
that no homomorphism matrices of graphs outside the class are produced.

For the pebbling comonad, the relevant lemmata were established in Section
4.2.

  3. 3.

Prove finite generation

Finite generation as defined in Section C.4 is necessary to obtain a finite
system of matrix equations in the subsequent step. This has been proven in
Section 4.3 for the pebbling comonad. Proving finite generation for arbitrary
comonads on $\mathsf{Graph}$ remains an open problem.

  4. 4.

Recover a matrix equation using algebraic and representation-theoretic
techniques.

The desired matrix equations can be obtained using the representation-
theoretic framework of [33], cf. Theorem 2.4. This resulted in Theorem C.15 in
the general case and Theorem 1.5 in the case of the pebbling comonad.

###  C.7 Outlook

We propose to apply this strategy to other comonads and hence to other graph
classes in the future. In [3], it was recently shown that for every graph
class $\mathcal{F}$ which is closed under isomorphism, finite coproducts, and
summands, cf. ibidem, there exists a comonad $\mathfrak{C}$ such that
$F\in\mathcal{F}$ if and only if $F$ admits a $\mathfrak{C}$-coalgebra. Thus,
it may be argued that all graph classes of interest can be characterised using
a comonad. However, the comonads constructed in [3] are only suitable for
counting homomorphisms and fail to represent the graph decompositions, which
were a conceptual ingredient for proving finite generation for Theorem 1.5.

For the future, we propose to eradicate issues with finiteness in Theorem C.15
to derive a general theorem characterising homomorphism indistinguishability
over graph classes characterised by a comonad in terms of matrix equations.


[1, 2]Sanskar Amgain These authors contributed equally to this work.

These authors contributed equally to this work.

1]NAAMII, Lalitpur, Nepal 2]Base Gene Therapeutics Limited, UK 3]University of
Aberdeen, UK 4]WEISS, and University College London, UK

# Investigation of Federated Learning Algorithms for Retinal Optical Coherence
Tomography Image Classification with Statistical Heterogeneity

<EMAIL_ADDRESS>Prashant Shrestha  Sophia Bano  Ignacio del Valle
Torres  Michael Cunniffe  Victor Hernandez  Phil Beales  Binod Bhattarai  [  [
[  [

###### Abstract

Purpose: We apply federated learning to train an OCT image classifier
simulating a realistic scenario with multiple clients and statistical
heterogeneous data distribution where data in the clients lack samples of some
categories entirely.

Methods: We investigate the effectiveness of FedAvg and FedProx to train an
OCT image classification model in a decentralized fashion, addressing privacy
concerns associated with centralizing data. We partitioned a publicly
available OCT dataset across multiple clients under IID and Non-IID settings
and conducted local training on the subsets for each client. We evaluated two
federated learning methods, FedAvg and FedProx for these settings.

Results: Our experiments on the dataset suggest that under IID settings, both
methods perform on par with training on a central data pool. However, the
performance of both algorithms declines as we increase the statistical
heterogeneity across the client data, while FedProx consistently performs
better than FedAvg in the increased heterogeneity settings.

Conclusion: Despite the effectiveness of federated learning in the utilization
of private data across multiple medical institutions, the large number of
clients and heterogeneous distribution of labels deteriorate the performance
of both algorithms. Notably, FedProx appears to be more robust to the
increased heterogeneity.

##  1 Introduction

Optical Coherence Tomography (OCT) is a widely used imaging technique for
diagnosing retinal abnormalities. An OCT image captures the different cross-
sectional layers of the retina, the thickness of which is important for
performing diagnosis [1]. A large number of works have explored the use of
artificial intelligence-based methods for assistive and automated diagnosis
from OCT images. The performance and generalizability of these models are
directly correlated with the amount of data available for training. However,
due to the privacy concerns that arise with sharing raw medical data outside
an institution, multi-institution collaboration for the curation of a large
central database remains a challenging problem.

Few works have recently explored the use of federated learning for OCT image
classification [2, 3] bypassing the need for a central data curation. For
instance, Gholami et al. [2] used federated learning to classify OCT images
into normal and age-related macular degeneration (AMD) categories while Ran et
al. [3] uses FedAvg to detect glaucoma samples. However, these works _do not
consider the statistical heterogeneous setting_ where the local data in
certain clients lack instances of some labels entirely. Furthermore, these
works experiment with a _small number of clients (3-4)_ sidestepping the
potential issues that arise when dealing with many clients such as the
increased communication overhead. Our work studies two prominent federated
learning methods on OCT images simulating such label heterogeneity with a
larger number of clients.

##  2 Methods

Figure 1: Illustration of FedAvg and FedProx training process. In each
communication round, a subset of clients trained locally on their respective
data splits with objective $L_{method}$ are selected for aggregation and used
to update the global model. The updated global model parameters, $w_{g}$ are
copied back to all the clients.

We divided the OCT dataset into disjoint subsets, each associated with a
separate client; simulating separate private data associated with distinct
medical institutions. As shown in Fig. 1, in each communication round, clients
initialize their local parameters with the current global model parameters
$w_{g}$. Each client, $i$ is trained on their respective data split with
objective $L_{i,method}$, the clients send back updated model parameters
$\\{w_{1},w_{2},\ldots w_{k}\\}$ to the server, which averages the parameters
across clients to obtain the updated global model. This process is repeated
until convergence. Since the sharing of information during the training
process occurs only at the level of model parameters, and not raw medical
data, federated learning achieves the aggregation of information from multiple
clients/institutions in a privacy-preserving manner.

To reduce communication overhead, a random subset of the clients is selected
for local training and aggregation. We experiment with two federated learning
methods: FedAvg [4] and FedProx [5]. FedAvg employs standard cross-entropy
loss, $F$ as the local training objective $L$, while FedProx additionally
penalizes the deviation of local parameters from the global model. Our study
utilized the OCT images from Kermany et al. [6] comprising 84,484 OCT images
categorized as CNV, DRUSEN, DME, or NORMAL (Fig. 2(b)).

##  3 Experiments

((a))

((b))

Figure 2: (a) Label distribution of 4 randomly selected clients for IID and 2
SPC Non-IID settings. For the IID setting, the distribution of labels across
clients is similar, while for 2 SPC, every client has multiple missing labels
and distinct distributions (b) Image samples from the dataset

We evaluated the effectiveness of two prominent federated learning methods,
FedAvg and FedProx for the OCT image classification task. With a total of 10
clients, we randomly selected 50% for aggregation to minimize communication
overhead. We experimented with two data configurations: IID and Non-IID
splits. For the IID setting, we randomly partitioned the data into 10 disjoint
subsets, each associated with a separate client. The data splits are thus less
likely to have missing labels and exhibit less heterogeneity in label
distribution across clients. For the Non-IID setting, we grouped the images by
labels and divided them into $10*k$ shards, each containing images belonging
to a single category. Each client is then assigned $k$ shards as the local
data. When $k<n$, $n$ being the total number of categories, this sharding
guarantees that no client possesses samples from all the categories, mirroring
real-life statistical heterogeneity where _not all medical centers have the
same distribution of patient data, and access to all the categories_. Fig 2(a)
illustrates the label distribution for a random subset of clients in both IID
and Non-IID settings.

For our classification model, we initialized the encoder with ResNet18
parameters pretrained on ImageNet. The backbone of the network is frozen and
only the classification head is fine-tuned on the OCT images. We use SGD
optimizer with a learning rate of 0.01 and batch size of 64 to train each
client model for 2 local epochs per communication round. The training is
performed for 100 communication rounds. We use average accuracy to evaluate
the performance of the models under IID and $k\in\\{2,3\\}$ Non-IID settings.

##  4 Results

As shown in Table 1, both FedAvg and FedProx perform well under IID settings
scoring on par with 98.933$\pm$0.27 obtained with centralized training.
However as the statistical heterogeneity across the clients is increased,
under the 3 shards per client (2SPC) and 2 shards per client (2SPC) Non-IID
settings, the performance of both approaches begins to deteriorate. Notably,
FedProx consistently outperforms FedAvg in our Non-IID settings. This is
attributed to the proximal loss term of FedProx that constrains the local
updates to be closer to the global model. This constraint mitigates the effect
of statistical heterogeneity leading to better performance.

Method | IID | 2SPC | 3SPC
---|---|---|---
FedAvg |  98.93 $\pm$ 0.05 |  87.22 $\pm$ 5.97 |  88.32 $\pm$ 2.30
FedProx ($\mu=0.2$) |  98.80 $\pm$ 0.10 |  90.57 $\pm$ 3.09 |  91.01 $\pm$ 2.51
Table 1: Accuracy on test dataset under IID, 2SPC, 3SPC data splits. Results
consist of mean and standard deviation over 3 random seeds.

##  5 Conclusion

Federated learning provides an effective solution to utilize private data
across multiple medical institutions, achieving high performance without
exposing raw medical data. Our study investigated the effectiveness of two
federated learning methods, FedAvg and FedProx under different heterogeneity
settings for OCT image classification. Both methods consistently performed
well across diverse scenarios. However, as the statistical heterogeneity
across clients increases, the performance of both approaches starts to
deteriorate. Notably, FedProx demonstrated greater robustness to the increased
heterogeneity than FedAvg. This suggests that in highly heterogeneous
settings, the application of federated learning for OCT image classification
necessitates more sophisticated approaches. Beyond statistical heterogenity,
future work can investigate the effectiveness of federated learning methods
with separate domain specific medical data in each client, along with dynamic
and evolving medical datasets.

## References

  *   * Al-Mujaini et al. [2013] Al-Mujaini, A., Wali, U.K., Azeem, S.: Optical coherence tomography: clinical applications in medical practice. Oman medical journal 28(2), 86 (2013)
  * Gholami et al. [2023] Gholami, S., Lim, J.I., Leng, T., Ong, S.S.Y., Thompson, A.C., Alam, M.N.: Federated learning for diagnosis of age-related macular degeneration. Frontiers in Medicine 10 (2023)
  * Ran et al. [2023] Ran, A.R., Wang, X., Chan, P.P., Wong, M.O., Yuen, H., Lam, N.M., Chan, N.C., Yip, W.W., Young, A.L., Yung, H.-W., et al.: Developing a privacy-preserving deep learning model for glaucoma detection: a multicentre study with federated learning. British Journal of Ophthalmology (2023)
  * McMahan et al. [2017] McMahan, B., Moore, E., Ramage, D., Hampson, S., Arcas, B.A.: Communication-efficient learning of deep networks from decentralized data. In: Artificial Intelligence and Statistics, pp. 1273–1282 (2017). PMLR
  * Li et al. [2020] Li, T., Sahu, A.K., Zaheer, M., Sanjabi, M., Talwalkar, A., Smith, V.: Federated optimization in heterogeneous networks. MLSys 2, 429–450 (2020)
  * Kermany et al. [2018] Kermany, D.S., Goldbaum, M., Cai, W., Valentim, C.C., Liang, H., Baxter, S.L., McKeown, A., Yang, G., Wu, X., Yan, F., et al.: Identifying medical diagnoses and treatable diseases by image-based deep learning. cell 172(5), 1122–1131 (2018)


# Parallel Spiking Neurons with High Efficiency and Ability to Learn Long-term
Dependencies

Wei Fang1,2, Zhaofei Yu4, Zhaokun Zhou3,2, Ding Chen5,2, Yanqi Chen1,2,
Zhengyu Ma2∗, Timothée Masquelier6, Yonghong Tian1,2,3

1School of Computer Science, Peking University, China
2Peng Cheng Laboratory, China
3School of Electronic and Computer Engineering, Shenzhen Graduate School,
Peking University, China
4School of Artificial Intelligence, Peking University, China
5Department of Computer Science and Engineering, Shanghai Jiao Tong
University, China
6Centre de Recherche Cerveau et Cognition (CERCO), UMR5549 CNRS - Univ.
Toulouse 3, France
Corresponding author

###### Abstract

Vanilla spiking neurons in Spiking Neural Networks (SNNs) use charge-fire-
reset neuronal dynamics, which can only be simulated serially and can hardly
learn long-time dependencies. We find that when removing reset, the neuronal
dynamics can be reformulated in a non-iterative form and parallelized. By
rewriting neuronal dynamics without reset to a general formulation, we propose
the Parallel Spiking Neuron (PSN), which generates hidden states that are
independent of their predecessors, resulting in parallelizable neuronal
dynamics and extremely high simulation speed. The weights of inputs in the PSN
are fully connected, which maximizes the utilization of temporal information.
To avoid the use of future inputs for step-by-step inference, the weights of
the PSN can be masked, resulting in the masked PSN. By sharing weights across
time-steps based on the masked PSN, the sliding PSN is proposed to handle
sequences of varying lengths. We evaluate the PSN family on simulation speed
and temporal/static data classification, and the results show the overwhelming
advantage of the PSN family in efficiency and accuracy. To the best of our
knowledge, this is the first study about parallelizing spiking neurons and can
be a cornerstone for the spiking deep learning research. Our codes are
available at https://github.com/fangwei123456/Parallel-Spiking-Neuron.

##  1 Introduction

Spiking Neural Networks (SNNs) are the next generation [1] of Artificial
Neural Networks (ANNs) using a lower abstraction of the biological neural
system. The spiking neurons are the key components of SNNs, which process
input currents with complex neuronal dynamics and fire spikes as outputs when
their membrane potential reaches the threshold. The SNNs use discrete spikes
to communicate between layers, which enable the event-driven computational
paradigm and have extremely high power efficiency in neuromorphic chips such
as True North[2], Loihi [3] and Tianjic [4]. Due to their high biological
plausibility, SNNs have been regarded by neuroscientists as useful tools for
analyzing, simulating, and learning the biological system [5, 6]. With the
introduction of deep learning methods [7, 8, 9, 10, 11,
kheradpisheh2020temporal], the performance of SNNs on real-world tasks has
been greatly improved and the applications of SNNs are expanded [12, 13]. As
the bridge between neuroscience and computational science, SNNs have attracted
more and more research interest in recent years.

The binary characteristic of spikes causes a lot of information losses, which
is the main factor leading to the lower accuracy of SNNs than ANNs. Among
efforts to bridge the performance gap, the improvement of spiking neurons is a
universal approach and has been widely explored [14, 15, 16, 17, 18, 19, 20].
However, previous works are limited to the existing serial charge-fire-reset
computing paradigm shown in Fig.1(a), which suffers from slow simulation speed
and can hardly learn long-term dependencies.

In this paper, we propose the Parallel Spiking Neuron (PSN) and its variants,
the masked PSN as well as the sliding PSN (SPSN). Fig.1(b) shows the
computational graph of the PSN, which uses a non-predecessor dependency method
to generate hidden states, leading to parallelizable neuronal dynamics and
extremely high simulation speed. Since the connections from inputs to $H[t]$
are direct, the PSN has a better learning ability for long-term dependencies
than other vanilla spiking neurons. To generate $H[t]$ immediately when $X[t]$
arrives, connections from $X[i],i\geq t+1$ can be masked, as the dotted line
shown in Fig.1(b), and then the neuron becomes the masked PSN. The parameters
of the masked PSN can also be set as time-invariant, which derives the sliding
PSN and is more flexible for input sequences with variable lengths. The main
contributions of this paper are as follows:

1) We analyze the influence of removing reset from the widely used charge-
fire-reset neuronal dynamics, showing that the spiking neuron can be
parallelized without reset.

2) By rewriting the neuronal dynamics without reset to a general formula, we
obtain the PSN with fully parallelizable neuronal dynamics. For step-by-step
serial forward and variable length sequences processing, the masked PSN and
the sliding PSN are also derived.

3) We compare the simulation speed of the PSN family and vanilla spiking
neurons, showing an overwhelming performance advantage of the parallel
neuronal dynamics. We evaluate the PSN family on sequential, static, and
neuromorphic data classification tasks and achieve higher accuracy than
previous spiking neurons.

Figure 1: The computational graphs of the vanilla spiking neuron and the
parallel spiking neuron. Figure (a) is cited from [14]. The dotted lines in
Figure (b) are weights that can be masked for step-by-step computation.
$X[t],S[t]$ are the input current and the output spike, and $H[t],V[t]$ are
hidden states at time-step $t$.

##  2 Related Work

###  2.1 Deep Learning for Spiking Neural Networks

The ANN-to-SNN conversion (ANN2SNN) [21, 22, 23, 24, 25] and the surrogate
training [26, 27] are the two main deep learning methods for SNNs. The ANN2SNN
method uses firing rates in SNNs to approximate the activations in ANNs by
converting weights from ANNs to SNNs with extra normalizations to reduce
errors. It achieves close performance to that of ANNs in the ImageNet dataset
[28], but requires many time-steps to estimate accurate firing rates. The
surrogate training method re-defines the gradient of the Heaviside function
used in the spike generating process by a smooth and differentiable surrogate
function. With the surrogate gradients, the SNNs can be trained by
backpropagation through time (BPTT) and gradient descent. The surrogate
training method is not based on rate coding as ANN2SNN, resulting in fewer
time-steps. However, the use of BPTT during training requires more training
time and memory usage than ANN2SNN.

###  2.2 Improvement of Spiking Neurons

Improving the spiking neurons inspired by biological mechanisms or techniques
from ANNs is a practical method. The Parametric Leaky Integrated-and-Fire
(PLIF) spiking neuron [14] incorporates learnable membrane time constants and
is able to combine the learning of both synapse weights and neuronal dynamics.
The gated LIF (GLIF) spiking neuron [15] adds learnable gates to fuse
different integration, decay, and reset methods, which enlarges the
representation space and increases the heterogeneity as well as adaptivity of
spiking neurons. The $k$-based Leaky Integrate-and-Fire (KLIF) neuron [17]
adds a learnable scaling factor to adjust the surrogate gradient and applies a
ReLU activation to avoid the negative membrane potential. The Multi-Level
Firing (MLF) unit [18] contains LIF neurons with thresholds in different
levels, which relieves the gradient vanishing problem and improves expression
ability. Ponghiran et al. [19] improve the spiking neurons by adding input
gates, forget gates, recurrent connections, and multi-bit outputs, which
reduce the mismatch between actual and surrogate gradients.

###  2.3 Improvement of Training Methods and Network Structure

Beyond the spiking neurons, modifying training techniques and network
structure are also effective methods to promote the performance of SNNs.

Normalization The threshold-dependent batch normalization (TDBN) [29]
normalizes features in both temporal and spatial dimensions and eliminates the
variance introduced by thresholds, which is a more suitable normalization
method than the batch normalization [30] for SNNs. Batch Normalization Through
Time (BNTT) [31] decouples BN along time-steps and uses time-wise parameters,
which estimates the distribution of temporal-variant inputs more precisely.
Temporal Effective Batch Normalization (TEBN) [32] adds learnable affine
transformation on each time-step and has a better ability to capture temporal
distributions.

Training Techniques Deep Continuous Local Learning (DECOLLE) [33] uses local
losses to avoid BPTT and implements online learning for SNNs. Online training
through time (OTTT) [34] approximates iterative gradients with temporal
dependencies by eligibility traces updating with time-steps, which also
achieves online training and only requires constant training memory agnostic
to time-steps. Sparse spiking gradient descent [35] only uses neurons whose
membrane potentials are inside a given range for backpropagation, which
achieves a speedup and memory reduction in the backward. The attention
mechanism on membrane potentials [36] eliminates redundant information,
decreases firing rates, and reduces overfitting. Dspike [37] uses the finite
difference to estimate gradients and achieves higher accuracy than plain
surrogate gradient methods. The temporal efficient training (TET) [38]
calculates loss at each time-step and averages them rather than using averaged
outputs to calculate loss, which improves the temporal scalability of SNNs.

Network Structure Spike-Element-Wise ResNet [39] solves the gradient
vanishing/exploding problems of the plain Spiking ResNet caused by sigmoid-
like surrogate functions, which successfully trained the first deep SNN with
more than 150 layers. Spikformer [40] modifies the softmax-based self-
attention in Transformer [41] to the spike-based formulation, which is
appropriate for SNNs and achieves state-of-the-art accuracy on the ImageNet
dataset.

###  2.4 Acceleration for Sequence Processing

Massively parallel computing with graphics processing units (GPUs) [42] is one
of the key factors that drives deep learning research. However, the serial
computing characteristic of the recurrent structure such as SNNs and Recurrent
Neural Networks (RNNs) can not fully exploit the parallel computing power of
GPUs. To accelerate the sequence processing, convolution-based methods [43,
44] discard the recurrent structure and employ the convolutional layers, which
are fully parallelized during training on GPUs. The gated impulse linear
recurrent network [45] is another solution that uses the Parallel Prefix Sum
(Scan) algorithm [46] to parallelize linear recurrences. The Transformer [41]
replaces recurrences with self-attention and positional encoding, which
achieves faster training speed than recurrent and convolutional encoder-
decoder structures.

##  3 Methods

In this section, we introduce the motivation of removing neuronal reset, the
idea of generalizing spiking neurons without reset as the PSN, and the
derivation of the masked PSN and the sliding PSN. Note that we use regular
letters such as $X$ to represent scalars, and bold letters such as $\bm{X}$ to
represent tensors.

###  3.1 Vanilla Spiking Neurons W/WO Reset

Spiking neurons in SNNs have rich neuronal dynamics, which endow SNNs with
temporal information processing ability. In most cases, the behaviors of
spiking neurons can be described by three discrete-time equations [14]:

$\displaystyle H[t]$ $\displaystyle=f(V[t-1],X[t]),$ (1) $\displaystyle S[t]$
$\displaystyle=\Theta(H[t]-V_{th}),$ (2) $\displaystyle V[t]$
$\displaystyle=\begin{cases}H[t]\cdot(1-S[t])+V_{reset}\cdot
S[t],\mathrm{~{}~{}hard~{}reset}\\\ H[t]-V_{th}\cdot
S[t],\mathrm{~{}~{}soft~{}reset}\\\ \end{cases},$ (3)

where $X[t]$ is the input current, $H[t]$ is the membrane potential after
charging but before firing, $V[t]$ is the membrane potential after firing, and
$S[t]$ is the output spike at time-step $t$. $V_{th}$ in Eq.(2) is the
threshold potential, and $V_{reset}$ in Eq.(3) is the reset potential.
$\Theta(x)$ is the Heaviside step function and $\Theta(x)=1$ for all $x\geq
0$, otherwise $\Theta(x)=0$. Eq.(1) is the neuronal charging equation, and $f$
is specific to different spiking neurons. After charging, $H[t]$ will be
compared to $V_{th}$ and determine whether to fire spikes, which is described
by Eq.(2). After firing, the membrane potential will be reset, as Eq.(3)
shows. Note that there are two kinds of resets, which are the hard reset and
the soft reset. If the neuron fires, the hard reset will set $V[t]$ to
$V_{reset}$, while the soft reset will decrease $V[t]$ by $V_{th}$. The hard
reset is widely used in surrogate training for better performance observed in
experiments [47], while the soft reset is preferred in ANN2SNN for lower
conversion errors. When simulating SNNs, the iterative process following
Eqs.(1)-(3) over time-steps is employed, which has a time complexity of
$\mathcal{O}(T)$, where $T$ is the number of time-steps.

For frequently-used spiking neurons, Eq.(1) is linear and can be reformulated
to a non-iterative equation if we can ignore Eq.(3) in some cases, e.g., in
the subthreshold regime $H[t]<V_{th}$ for all $t$ during a period. Then
$V[t]=H[t]$ at all $t$ and we will only use $H[t]$ in this case. More
specifically, we take the Integrate-and-Fire (IF) neuron and the Leaky
Integrate-and-Fire (LIF) neuron as examples. For simplicity, suppose $H[-1]=0$
for both types of neurons, and the resting potential for the LIF neuron is
$0$. The neuronal charge equation, or Eq.(1), for the IF neuron is

$\displaystyle H[t]=H[t-1]+X[t].$ (4)

And Eq.(4) can be easily reformulated as

$\displaystyle H[t]=\sum_{i=0}^{t}X[i].$ (5)

The neuronal charge equation for the LIF neuron is

$\displaystyle H[t]=(1-\frac{1}{\tau_{m}})\cdot H[t-1]+\frac{1}{\tau_{m}}\cdot
X[t],$ (6)

where $\tau_{m}$ is the membrane time constant. Eq.(6) can also be
reformulated as

$\displaystyle
H[t]=\frac{1}{\tau_{m}}\cdot\sum_{i=0}^{t}(1-\frac{1}{\tau_{m}})^{t-i}\cdot
X[i].$ (7)

With the non-iterative equations, $H[t]$ at all time-steps, or
$\bm{H}=\\{H[0],H[1],...,H[T-1]\\}$, can be calculated in parallel. For each
$H[t]$ obtained by the cumulative sum operation, the time complexity can be as
low as $\mathcal{O}(\mathrm{log}(t))\leq\mathcal{O}(\mathrm{log}(T))$ when
using the Parallel Prefix Sum (Scan) algorithm. When $\bm{H}$ is given, the
element-wise operation Eq.(2) can also be applied on the whole $\bm{H}$ and
the time complexity is $\mathcal{O}(1)$ in devices that support parallel
computing such as the CUDA devices. Thus, the whole time complexity for the
neuronal charging and firing reduces to $\mathcal{O}(\mathrm{log}(T))$ when we
can ignore the neuronal resetting.

So, how can we ignore neuronal resetting at all time-steps? A possible method
is setting $V_{th}=+\infty$ to make Eq.(3) become an identity function.
However, the neuron can only output $0$ in this case, which is meaningless.
Another method is to remove the neuronal resetting directly from the neuronal
dynamics, which means only using Eq.(1) and Eq.2. This crude and simple method
may raise further concerns that the absence of neuronal resetting will cause
the neuron to fire uninterruptedly. Note that the firing rate is determined by
inputs, charging, thresholds, and resetting. Thus, the effect of only removing
reset on firing rates is not deterministic. In addition, the experiment
results in the next section will show that the uninterrupted firing issue will
not occur.

###  3.2 Parallel Spiking Neuron

After removing neuronal resetting, the charging equation of the spiking neuron
can be formulated into a non-iterative equation, and the solution for $H[t]$
becomes a cumulative sum operation. It can be found that Eq.(5) and Eq.(7) are
two specific cases of a linear combination of $X[i]$. More specifically, we
formulate a general linear combination as

$\displaystyle H[t]=\sum_{i=0}^{T-1}W_{t,i}\cdot X[i],$ (8)

where $W_{t,i}$ is the weight between $X[i]$ and $H[t]$. Then,
$W_{t,i}=\Theta(t-i)$ for the IF neuron without reset, and
$W_{t,i}=\frac{1}{\tau_{m}}(1-\frac{1}{\tau_{m}})^{t-i}\cdot\Theta(t-i)$ for
the LIF neuron without reset. Based on the above analysis, we propose the
Parallel Spiking Neuron (PSN), whose neuronal dynamics are as follows:

$\displaystyle\bm{H}$
$\displaystyle=\bm{W}\bm{X},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\bm{W}\in\mathbb{R}^{T\times
T},\bm{X}\in\mathbb{R}^{T\times N}$ (9) $\displaystyle\bm{S}$
$\displaystyle=\Theta(\bm{H}-\bm{B}),~{}~{}~{}~{}~{}\bm{B}\in\mathbb{R}^{T},\bm{S}\in\\{0,1\\}^{T\times
N}$ (10)

where $\bm{X}$ is the input sequence, $\bm{W}$ is the learnable weight matrix,
$\bm{H}$ is the hidden state sequence, $\bm{B}$ is the learnable threshold,
and $\bm{S}$ is the binary output spike sequence. $N$ is the number of batch
size, and $T$ is the number of time-steps. Note that $\bm{W}$ is a $T\times T$
matrix, indicating that $H[t]$ can integrate the information from all time-
steps directly, rather than that of the last one or two time-steps. Thus, the
PSN is a $T$-order neuron. No iterative equation is used, and the computation
of PSN is fully parallel. The time complexity of the PSN is determined by
Eq.(9), which is a matrix-matrix multiplication. The simulation of PSN is much
faster than that of vanilla spiking neurons because the matrix-matrix
multiplication is highly optimized in linear algebra libraries such as Intel
MKL and cuBLAS. BPTT is still used during training, but also in a parallel,
rather than a serial formulation.

###  3.3 $k$-Order Masked Parallel Spiking Neuron

The high-order characteristic of the PSN is also a double-edged sword. The
output sequence $\bm{S}$ can only be generated when all $X[t]$ have arrived,
indicating that each spiking neuron layer will increase latency $T$ to the
whole SNN. Such a latency problem is also reported on some time-to-first-spike
SNNs [48, 49, 50]. Note that although the latency is unavoidable, the
throughput can still approach the plain SNN with the multi-stage pipeline
[48]. To solve the latency problem of the PSN, we add a mask $\bm{M}_{k}$ to
multiply $\bm{W}$ element-wisely before generating $\bm{H}$ and propose the
$k$-order masked PSN, whose neuronal charging equation is

$\displaystyle\bm{H}$
$\displaystyle=(\bm{W}\cdot\bm{M}_{k})\bm{X},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\bm{W}\in\mathbb{R}^{T\times
T},\bm{M}_{k}\in\mathbb{R}^{T\times T},\bm{X}\in\mathbb{R}^{T\times N}$ (11)

where $\bm{M}_{k}$ is defined as

$\displaystyle\bm{M}_{k}[i][j]$ $\displaystyle=\begin{cases}1,~{}~{}j\leq
i\leq j+k-1\\\ 0,\mathrm{otherwise}\end{cases}.$ (12)

With the mask $\bm{M}_{k}$, $H[t]$ depends only on the latest $k$ inputs
$\\{X[t-k+1],...,X[t-1],X[t]\\}$ (we assume that $X[i]=0$ for all $i<0$) and
$S[t]$ can be computed and sent to the next layer once $X[t]$ is received.

When training the $k$-order masked PSN, the progressive masking method is
employed, which uses an all-ones matrix $\bm{1}$ as the mask first, and
gradually replaces $\bm{1}$ with $\bm{M}_{k}$. More specifically, we use

$\displaystyle\bm{M}_{k}(\lambda)=\lambda\cdot\bm{M}_{k}+(1-\lambda)\cdot\bm{1}$
(13)

as the mask, and $\lambda$ increases from $0$ to $1$ during training. In the
early stages of training, $\bm{M_{k}}(\lambda)\approx\bm{1}$, indicating that
the future information can be exploited to provide appropriate primary
parameters and help the network converge.

###  3.4 $k$-Order Sliding Parallel Spiking Neuron

The parameters of the PSN and the masked PSN are time-wise, which require
extra operations for processing sequences with variable lengths. To solve this
issue, we modify the parameters to be shared across time-steps and propose the
$k$-order sliding PSN, whose neuronal dynamics are

$\displaystyle H[t]$ $\displaystyle=\sum_{i=0}^{k-1}W_{i}\cdot X[t-k+1+i],$
(14) $\displaystyle S[t]$ $\displaystyle=\Theta(H[t]-V_{th}),$ (15)

where $\bm{W}=[W_{0},W_{1},...,W_{k-1}]\in\mathbb{R}^{k}$ is the learnable
weight, $X[j]=0$ for all $j<0$, and $V_{th}$ is the learnable threshold.
Similar to the masked PSN, $H[t]$ of the sliding PSN depends only on the
latest $k$ inputs as well. More specifically, the weight slides over inputs
and generates hidden states, which is a standard 1D convolution operation.
Additionally, it can also be implemented by matrix-matrix multiplication. When
the input sequence $\bm{X}\in\mathbb{R}^{T\times N}$ arrives and its length
$T$ is known, the matrix $\bm{A}\in\mathbb{R}^{T\times T}$ can be generated as

$\displaystyle\bm{A}[i][j]$
$\displaystyle=\begin{cases}W_{k-1-i+j},~{}~{}i+1-k\leq j\leq i\\\
0,\mathrm{otherwise}\end{cases},$ (16)

then $\bm{H}=\bm{A}\bm{X}$. According to the results of our experiments, using
matrix-matrix multiplication is faster than using convolution.

###  3.5 Summary of the PSN Family

As a summary, Fig.2 shows the comparison of the parameter number $n_{param}$
and the generation of hidden states of the PSN family. Remarkably, using the
PSN in deep SNNs will add a negligible number of parameters because $T$ is
small in directly trained SNNs. For example, using the PSN with $T=4$ will add
340 and 200 parameters in spiking ResNet-18 and VGG-11, which only increase
$0.00291\%$ and $0.00015\%$ parameters of the original SNNs. Additionally, the
PSN family has only one hidden state $\bm{H}$, making them consume less memory
than the vanilla spiking neurons with at least two hidden states
$\bm{H},\bm{V}$ during training, which are analyzed and verified by
experiments in the supplementary materials.

Figure 2: Comparison of the PSN family on the parameter number $n_{param}$ and
the generation of hidden states.

##  4 Experiments

In this section, we report the experiment results of evaluating the PSN family
in aspects of simulation performance and temporal/static data classifying
accuracy. All experiments are based on the SpikingJelly [51] framework. The
details of the training are provided in the supplementary materials.

###  4.1 Simulation Speed Benchmark

We evaluated the simulation speed benchmark of the PSN and the LIF neuron,
which is widely used in deep SNNs and acts as the baseline of the vanilla
spiking neurons. Note that the implementations and simulation speed of the
masked PSN and the sliding PSN are almost identical to the PSN. Thus, we take
the PSN as the benchmark for the PSN family. The benchmark is running in two
modes, which are inference and training. In inference, we use the PyTorch
just-in-time (JIT) compiling to fuse Eq.(1)-(3) of the LIF neuron across all
time-steps into a single function to accelerate by avoiding the calling
overhead of too many tiny CUDA kernels. However, PyTorch JIT does not support
modifying the backward for a fused forward function, so the surrogate method
used in the backward of Eq.(2) breaks the fusion of all three equations in
training. Thus, we have to use three smaller JIT functions at each time-step
for the LIF neuron in training. For the PSN, there is no difference between
the implementations for inference and training, which both involve a matrix-
matrix multiplication and an element-wise operation. We tested the neuron
numbers $N=2^{8},2^{12},2^{16},2^{20}$ and the time-step numbers
$T=2,4,8,16,32,64$, which are typical options for deep SNNs. Denote the
executing time of two neurons as $t_{PSN}$ and $t_{LIF}$ respectively, then
the ratio $\frac{t_{LIF}}{t_{PSN}}$ is shown in Fig.3. It can be found that in
most cases, the simulation of the PSN is much faster than that of the LIF
neuron.

(a) Forward (inference) (b) Forward and backward (training)

Figure 3: The executing time ratio of the LIF neuron and the PSN for one
forward in inference and one forward and backward in training.

###  4.2 Learning Long-term Dependencies

To verify the learning ability for long-term dependencies of different spiking
neurons, we evaluated their performance on classifying sequential CIFAR10 and
CIFAR100. In these tasks, the image will be sent to the SNN column by column,
which is similar to how humans read from left to right. The sequential image
classification task is widely used to evaluate the learning ability of a
network for long-term dependencies. The network structure is modified from
[14] by replacing 2D convolutional/max pooling layers with 1D convolutional
and average pooling layers, and removing the voting layers. Due to the column-
by-column input, $T$ is $32$, which is the width of the image. The compared
neurons include PSN, $32$-order masked PSN, $32$-order sliding PSN, GLIF with
channel-wise parameters, KLIF, PLIF, and the LIF neuron w/wo reset. We use the
optimal reset options determined by ablation experiments for different
neurons, which are using soft reset and detach reset [52] for the PLIF and LIF
neurons, and using soft reset for the KLIF neuron. The GLIF neuron uses
learnable gates to control all options, which do not need to be set manually.
For the masked PSN, $\lambda=\mathrm{min}(1,8\cdot epoch/(epochs-1))$, where
$epoch$ denotes the current training epoch, and $epochs$ denotes the total
number of epochs.

Dataset Neuron | PSN | masked PSN | sliding PSN | GLIF[15] | KLIF[17] | PLIF[14] | LIF | LIF wo reset
---|---|---|---|---|---|---|---|---
Sequential CIFAR10 | 88.45 | 85.81 | 86.70 | 83.66 | 83.26 | 83.49 | 81.50 | 79.50
Sequential CIFAR100 | 62.21 | 60.69 | 62.11 | 58.92 | 57.37 | 57.55 | 55.45 | 53.33

Table 1: The test accuracy (%) of different spiking neurons on sequential
CIFAR.

Tab.1 shows the accuracy of all neurons on sequential CIFAR. It can be found
that the rank of accuracy is PSN $>$ sliding PSN $>$ masked PSN $>$ GLIF $>$
PLIF $>$ KLIF $>$ LIF $>$ LIF wo reset, which is as expected. The PSN exploits
information from all time-steps and has the highest learning ability, but it
may be criticized for unfair comparison because it can read all columns at
once. Therefore, the accuracy of the masked PSN and the sliding PSN is more
convincing because they can still run in a step-by-step mode. Due to the
introduction of high-order information, they work better than any other
traditional spiking neurons. Meanwhile, the GLIF neuron is the best among
traditional spiking neurons, indicating that fusing different kinds of
charging, decaying, and resetting increases learning ability. The LIF neuron
without reset has lower accuracy than the LIF neuron, indicating that the
direct removal of reset drops the learning ability of the vanilla spiking
neurons. Thus, using high-order information as the PSN family is necessary to
fix this decline.

(a) Sequential CIFAR10

(b) Sequential CIFAR100

Figure 4: The accuracy curves and firing rates of different neurons for
classifying sequential CIFAR.

To make clear comparisons, we plotted the accuracy curves during training and
firing rates of each spiking neuron layer of trained SNNs with different
spiking neurons, as shown in Fig.4. These accuracy curves are the moving
averages of 8 epochs, and the light-colored curves show the original accuracy.
The curves of the PSN family are almost always above other curves, which
clearly show their superior performance on convergence speed and final
accuracy. Note that there is a sharp drop in the curves of the masked PSN
around epoch 32, which is caused by the mask becoming completely binary. The
firing rates indicate that uninterrupted firing does not occur with the
removal of neuronal resetting. In general, the PSN and the masked PSN are more
easily to be activated and show higher firing rates than others. Meanwhile,
the firing rates of the sliding PSN and the LIF neuron without reset are
slightly higher than the vanilla spiking neurons. Thus, we conclude that
removal of reset does increase firing rates, but the increment degree depends
on the neuron’s structure. On the other hand, considering the huge improvement
in simulation efficiency and long-term dependencies learning ability, the
slightly higher firing rate is acceptable.

Figure 5: The accuracy-order curve on the sequential CIFAR100.

Fig.5 shows the accuracy-order curve on the sequential CIFAR100 dataset for
the masked PSN and sliding PSN, with the highest accuracy marked by a red ★.
It can be seen that when the order increases from 1 to 4, the accuracy
improves quickly. When we continue to increase the order, the improvement is
marginal and the accuracy even drops. The highest accuracy is obtained by a
large, but not the largest order, which is 20 and 31 for the masked PSN and
sliding PSN, respectively. But in most cases, high accuracy can be guaranteed
with $k=T$. Remarkably, the curve of sliding PSN is almost always above the
curve of masked PSN, and excels that of the PSN when $k=21$ and $k=31$,
indicating that the time-invariant parameters have better generalization
ability to process temporal information.

###  4.3 Static and Neuromorphic Data Classification

Static and neuromorphic data classification tasks are also common benchmarks
for SNNs. We evaluated the PSN family on the static CIFAR10, ImageNet
datasets, and the neuromorphic CIFAR10-DVS [53] dataset. The number of time-
steps for static datasets is small, and the latency is not the main issue.
Thus, we use the PSN in CIFAR10 and ImageNet. According to the results on
sequential CIFAR, the sliding PSN uses fewer parameters but gets higher
accuracy, and we will use it for CIFAR10-DVS. The results are listed in Tab.2.

Dataset | Method | Spiking Network | Time-steps | Accuracy(%)
---|---|---|---|---
CIFAR10 | Dspike[37] | Modified ResNet-18 | 6 | 94.25
TET[38] | ResNet-19 | 6 | 94.50
TDBN[29] | ResNet-19 | 6 | 93.16
TEBN[32] | ResNet-19 | 6 | 94.71
PLIF[14] | PLIF Net | 8 | 93.50
KLIF[17] | Modified PLIF Net | 10 | 92.52
GLIF[15] | ResNet-19 | 6 | 95.03
PSN | Modified PLIF Net | 4 | 95.32
ImageNet | Dspike[37] | ResNet-34 | 6 | 68.19
VGG-16 | 5 | 71.24
TET[38] | SEW ResNet-34 | 4 | 68.00
TDBN[29] | ResNet-34 with double channels | 6 | 67.05
TEBN[32] | SEW ResNet-34 | 4 | 68.28
GLIF[15] | ResNet-34 | 4 | 67.52
SEW ResNet[39] | SEW ResNet-18 | 4 | 63.18
SEW ResNet-34 | 4 | 67.04
PSN | SEW ResNet-18 | 4 | 67.63
| SEW ResNet-34 | 4 | 70.54
CIFAR10-DVS | Dspike[37] | ResNet-18 | 10 | 75.40
TET[38] | VGG | 10 | 83.17
TDBN[29] | ResNet-19 | 10 | 67.80
TEBN[32] | VGG | 10 | 84.90
PLIF[14] | PLIF Net | 20 | 74.80
KLIF[17] | Modified PLIF Net | 15 | 70.90
GLIF[15] | Wide 7B Net | 16 | 78.10
SEW ResNet[39] | Wide 7B Net | 16 | 74.40
sliding PSN ($k=2$) | VGG | 4, 8, 10 | 82.30, 85.30, 85.90

Table 2: Comparison of the PSN family and other methods.

CIFAR10 We modified the network structure of [14], or the PLIF Net, by using
average pooling and removing voting layers. Our PSN uses the shortest time-
steps $T=4$ and achieves the highest accuracy of 95.32%.

ImageNet We used the SEW ResNet [39] for classifying the ImageNet dataset. We
verified the performance of the SPN on SEW ResNet-18/34 and got a stable 3+%
accuracy improvement over the original SEW ResNet with the IF neuron. Our
method achieves 70.54% accuracy and is only second to the Dspike method using
one more time-step and VGG-16, which has a number of parameters 6.3 times as
many as ours.

CIFAR10-DVS We used the VGG structure from [38]. Considering the fact that the
DVS data contains temporal features and $T$ is larger than the static
datasets, we use the $2$-order sliding PSN. We achieve 82.3 % accuracy with
$T=4$, which is the first work to get 80+% accuracy in such few time-steps.
When increasing time-steps to $T=8$, we achieve 85.30% accuracy, which is 2+%
higher than that of the previous SOTA method [32] and we use two fewer time-
steps. With $T=10$, the accuracy is further improved to 85.90%.

##  5 Conclusion

In this paper, we remove reset from neuronal dynamics of vanilla spiking
neurons and propose the PSN family, including the PSN, the masked PSN, and the
sliding PSN. The PSN family can be simulated in parallel, which greatly
accelerates the training of deep SNNs. Weights of inputs in the PSN family can
be fully connected, or masked/shared with a custom order. Experimental results
of simulation speed and temporal/static data classification verify the
efficiency and accuracy of the PSN family when compared with previous spiking
neurons. Our work may be a milestone as well as a cornerstone for modern
spiking neurons.

## Acknowledgments and Disclosure of Funding

This work is supported by grants from the National Natural Science Foundation
of China (62027804, 61825101, 62088102, and 62176003), the major key project
of the Peng Cheng Laboratory (PCL2021A13), and the Agence Nationale de la
Recherche under Grant ANR-20-CE23-0004-04 DeepSee.

## References

  * [1] Wolfgang Maass.  Networks of spiking neurons: the third generation of neural network models.  Neural Networks, 10(9):1659–1671, 1997.
  * [2] Paul A Merolla, John V Arthur, Rodrigo Alvarez-Icaza, Andrew S Cassidy, Jun Sawada, Filipp Akopyan, Bryan L Jackson, Nabil Imam, Chen Guo, Yutaka Nakamura, et al.  A million spiking-neuron integrated circuit with a scalable communication network and interface.  Science, 345(6197):668–673, 2014.
  * [3] Mike Davies, Narayan Srinivasa, Tsung-Han Lin, Gautham Chinya, Yongqiang Cao, Sri Harsha Choday, Georgios Dimou, Prasad Joshi, Nabil Imam, Shweta Jain, Yuyun Liao, Chit-Kwan Lin, Andrew Lines, Ruokun Liu, Deepak Mathaikutty, Steven McCoy, Arnab Paul, Jonathan Tse, Guruguhanathan Venkataramanan, Yi-Hsin Weng, Andreas Wild, Yoonseok Yang, and Hong Wang.  Loihi: A neuromorphic manycore processor with on-chip learning.  IEEE Micro, 38(1):82–99, 2018.
  * [4] Jing Pei, Lei Deng, Sen Song, Mingguo Zhao, Youhui Zhang, Shuang Wu, Guanrui Wang, Zhe Zou, Zhenzhi Wu, Wei He, et al.  Towards artificial general intelligence with hybrid tianjic chip architecture.  Nature, 572(7767):106–111, 2019.
  * [5] Chris Eliasmith, Terrence C. Stewart, Xuan Choo, Trevor Bekolay, Travis DeWolf, Yichuan Tang, and Daniel Rasmussen.  A large-scale model of the functioning brain.  Science, 338(6111):1202–1205, 2012.
  * [6] Marcel Stimberg, Romain Brette, and Dan FM Goodman.  Brian 2, an intuitive and efficient neural simulator.  eLife, 8:e47314, 2019.
  * [7] Emre O Neftci, Hesham Mostafa, and Friedemann Zenke.  Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks.  IEEE Signal Processing Magazine, 36(6):51–63, 2019.
  * [8] Amirhossein Tavanaei, Masoud Ghodrati, Saeed Reza Kheradpisheh, Timothée Masquelier, and Anthony Maida.  Deep learning in spiking neural networks.  Neural Networks, 111:47–63, 2019.
  * [9] Yufei Guo, Xuhui Huang, and Zhe Ma.  Direct learning-based deep spiking neural networks: a review.  Frontiers in Neuroscience, 17:1209795, 2023.
  * [10] Yongqiang Cao, Yang Chen, and Deepak Khosla.  Spiking deep convolutional neural networks for energy-efficient object recognition.  International Journal of Computer Vision, 113(1):54–66, 2015.
  * [11] Bodo Rueckauer, Iulia-Alexandra Lungu, Yuhuang Hu, Michael Pfeiffer, and Shih-Chii Liu.  Conversion of continuous-valued deep networks to efficient event-driven networks for image classification.  Frontiers in Neuroscience, 11:682, 2017.
  * [12] Kaushik Roy, Akhilesh Jaiswal, and Priyadarshini Panda.  Towards spike-based machine intelligence with neuromorphic computing.  Nature, 575(7784):607–617, 2019.
  * [13] Rui Yuan, Qingxi Duan, Pek Jun Tiw, Ge Li, Zhuojian Xiao, Zhaokun Jing, Ke Yang, Chang Liu, Chen Ge, Ru Huang, et al.  A calibratable sensory neuron based on epitaxial vo2 for spike-based neuromorphic multisensory system.  Nature Communications, 13(1):1–12, 2022.
  * [14] Wei Fang, Zhaofei Yu, Yanqi Chen, Timothée Masquelier, Tiejun Huang, and Yonghong Tian.  Incorporating learnable membrane time constant to enhance learning of spiking neural networks.  In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pages 2661–2671, 2021.
  * [15] Xingting Yao, Fanrong Li, Zitao Mo, and Jian Cheng.  GLIF: A unified gated leaky integrate-and-fire neuron for spiking neural networks.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems (NeurIPS), 2022\.
  * [16] Ilyass Hammouamri, Timothée Masquelier, and Dennis George Wilson.  Mitigating catastrophic forgetting in spiking neural networks through threshold modulation.  Transactions on Machine Learning Research, 2022.
  * [17] Chunming Jiang and Yilei Zhang.  Klif: An optimized spiking neuron unit for tuning surrogate gradient slope and membrane potential, 2023.
  * [18] Lang Feng, Qianhui Liu, Huajin Tang, De Ma, and Gang Pan.  Multi-level firing with spiking ds-resnet: Enabling better and deeper directly-trained spiking neural networks.  In Lud De Raedt, editor, International Joint Conferences on Artificial Intelligence Organization (IJCAI), pages 2471–2477, 7 2022.  Main Track.
  * [19] Wachirawit Ponghiran and Kaushik Roy.  Spiking neural networks with improved inherent recurrence dynamics for sequential learning.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), volume 36, pages 8001–8008, 2022.
  * [20] Yufei Guo, Xinyi Tong, Yuanpei Chen, Liwen Zhang, Xiaode Liu, Zhe Ma, and Xuhui Huang.  Recdis-snn: Rectifying membrane potential distribution for directly training spiking neural networks.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 326–335, June 2022.
  * [21] Bing Han, Gopalakrishnan Srinivasan, and Kaushik Roy.  Rmp-snn: Residual membrane potential neuron for enabling deeper high-accuracy and low-latency spiking neural network.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 13558–13567, 2020.
  * [22] Jianhao Ding, Zhaofei Yu, Yonghong Tian, and Tiejun Huang.  Optimal ann-snn conversion for fast and accurate inference in deep spiking neural networks.  In Zhi-Hua Zhou, editor, Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21, pages 2328–2336. International Joint Conferences on Artificial Intelligence Organization, 8 2021.  Main Track.
  * [23] Tong Bu, Wei Fang, Jianhao Ding, PengLin Dai, Zhaofei Yu, and Tiejun Huang.  Optimal ann-snn conversion for high-accuracy and ultra-low-latency spiking neural networks.  In International Conference on Learning Representations (ICLR), 2021\.
  * [24] Shikuang Deng and Shi Gu.  Optimal conversion of conventional artificial neural networks to spiking neural networks.  In International Conference on Learning Representations (ICLR), 2021\.
  * [25] Zecheng Hao, Tong Bu, Jianhao Ding, Tiejun Huang, and Zhaofei Yu.  Reducing ann-snn conversion error through residual membrane potential.  Proceedings of the AAAI Conference on Artificial Intelligence, 37(1):11–21, Jun. 2023.
  * [26] Yujie Wu, Lei Deng, Guoqi Li, Jun Zhu, and Luping Shi.  Spatio-temporal backpropagation for training high-performance spiking neural networks.  Frontiers in Neuroscience, 12, 2018.
  * [27] Sumit Bam Shrestha and Garrick Orchard.  Slayer: Spike layer error reassignment in time.  In Advances in Neural Information Processing Systems (NeurIPS), pages 1419–1428, 2018.
  * [28] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, et al.  Imagenet large scale visual recognition challenge.  International Journal of Computer Vision, 115(3):211–252, 2015\.
  * [29] Hanle Zheng, Yujie Wu, Lei Deng, Yifan Hu, and Guoqi Li.  Going deeper with directly-trained larger spiking neural networks.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), volume 35, pages 11062–11070, 2021.
  * [30] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In International Conference on Machine Learning (ICML), pages 448–456. PMLR, 2015.
  * [31] Youngeun Kim and Priyadarshini Panda.  Revisiting batch normalization for training low-latency deep spiking neural networks from scratch.  Frontiers in Neuroscience, 15, 2021.
  * [32] Chaoteng Duan, Jianhao Ding, Shiyan Chen, Zhaofei Yu, and Tiejun Huang.  Temporal effective batch normalization in spiking neural networks.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems (NeurIPS), 2022\.
  * [33] Jacques Kaiser, Hesham Mostafa, and Emre Neftci.  Synaptic plasticity dynamics for deep continuous local learning (decolle).  Frontiers in Neuroscience, 14, 2020.
  * [34] Mingqing Xiao, Qingyan Meng, Zongpeng Zhang, Di He, and Zhouchen Lin.  Online training through time for spiking neural networks.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems (NeurIPS), 2022\.
  * [35] Nicolas Perez-Nieves and Dan F. M. Goodman.  Sparse spiking gradient descent.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems (NeurIPS), 2021\.
  * [36] Man Yao, Guangshe Zhao, Hengyu Zhang, Yifan Hu, Lei Deng, Yonghong Tian, Bo Xu, and Guoqi Li.  Attention spiking neural networks.  IEEE Transactions on Pattern Analysis and Machine Intelligence, pages 1–18, 2023.
  * [37] Yuhang Li, Yufei Guo, Shanghang Zhang, Shikuang Deng, Yongqing Hai, and Shi Gu.  Differentiable spike: Rethinking gradient-descent for training spiking neural networks.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems (NeurIPS), 2021\.
  * [38] Shikuang Deng, Yuhang Li, Shanghang Zhang, and Shi Gu.  Temporal efficient training of spiking neural network via gradient re-weighting.  In International Conference on Learning Representations (ICLR), 2022\.
  * [39] Wei Fang, Zhaofei Yu, Yanqi Chen, Tiejun Huang, Timothée Masquelier, and Yonghong Tian.  Deep residual learning in spiking neural networks.  Advances in Neural Information Processing Systems (NeurIPS), 34, 2021.
  * [40] Zhaokun Zhou, Yuesheng Zhu, Chao He, Yaowei Wang, Shuicheng YAN, Yonghong Tian, and Li Yuan.  Spikformer: When spiking neural network meets transformer.  In International Conference on Learning Representations (ICLR), 2023\.
  * [41] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in Neural Information Processing Systems (NeurIPS), 30, 2017.
  * [42] Rajat Raina, Anand Madhavan, and Andrew Y Ng.  Large-scale deep unsupervised learning using graphics processors.  In International Conference on Machine Learning (ICML), pages 873–880, 2009.
  * [43] Nal Kalchbrenner, Lasse Espeholt, Karen Simonyan, Aaron van den Oord, Alex Graves, and Koray Kavukcuoglu.  Neural machine translation in linear time.  arXiv preprint arXiv:1610.10099, 2016.
  * [44] Jonas Gehring, Michael Auli, David Grangier, Denis Yarats, and Yann N Dauphin.  Convolutional sequence to sequence learning.  In International Conference on Machine Learning (ICML), pages 1243–1252. PMLR, 2017.
  * [45] Eric Martin and Chris Cundy.  Parallelizing linear recurrent neural nets over sequence length.  In International Conference on Learning Representations (ICLR), 2018\.
  * [46] Mark Harris, Shubhabrata Sengupta, and John D Owens.  Parallel prefix sum (scan) with cuda.  GPU gems, 3(39):851–876, 2007.
  * [47] Eimantas Ledinauskas, Julius Ruseckas, Alfonsas Juršėnas, and Giedrius Buračas.  Training Deep Spiking Neural Networks.  arXiv preprint arXiv:2006.04436, 2020.
  * [48] Seongsik Park, Seijoon Kim, Byunggook Na, and Sungroh Yoon.  T2fsnn: deep spiking neural networks with time-to-first-spike coding.  In 2020 57th ACM/IEEE Design Automation Conference, pages 1–6. IEEE, 2020.
  * [49] Dongwoo Lew, Kyungchul Lee, and Jongsun Park.  A time-to-first-spike coding and conversion aware training for energy-efficient deep spiking neural network processor design.  In Proceedings of the 59th ACM/IEEE Design Automation Conference, pages 265–270, 2022.
  * [50] Lina Bonilla, Jacques Gautrais, Simon Thorpe, and Timothée Masquelier.  Analyzing time-to-first-spike coding schemes.  Frontiers in Neuroscience, 16, 2022.
  * [51] Wei Fang, Yanqi Chen, Jianhao Ding, Zhaofei Yu, Timothée Masquelier, Ding Chen, Liwei Huang, Huihui Zhou, Guoqi Li, and Yonghong Tian.  Spikingjelly: An open-source machine learning infrastructure platform for spike-based intelligence.  Science Advances, 9(40):eadi1480, 2023.
  * [52] Friedemann Zenke and Tim P Vogels.  The remarkable robustness of surrogate gradient learning for instilling complex function in spiking neural networks.  BioRxiv, 2020.
  * [53] Hongmin Li, Hanchao Liu, Xiangyang Ji, Guoqi Li, and Luping Shi.  Cifar10-dvs: An event-stream dataset for object classification.  Frontiers in Neuroscience, 11, 2017.
  * [54] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proceedings of the IEEE international conference on computer vision, pages 1026–1034, 2015.
  * [55] Zhuang Liu, Hanzi Mao, Chao-Yuan Wu, Christoph Feichtenhofer, Trevor Darrell, and Saining Xie.  A convnet for the 2020s.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 11976–11986, 2022.
  * [56] Ilya Loshchilov and Frank Hutter.  Sgdr: Stochastic gradient descent with warm restarts.  arXiv preprint arXiv:1608.03983, 2016.
  * [57] Hongyi Zhang, Moustapha Cisse, Yann N Dauphin, and David Lopez-Paz.  mixup: Beyond empirical risk minimization.  arXiv preprint arXiv:1710.09412, 2017.
  * [58] Sangdoo Yun, Dongyoon Han, Seong Joon Oh, Sanghyuk Chun, Junsuk Choe, and Youngjoon Yoo.  Cutmix: Regularization strategy to train strong classifiers with localizable features.  In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pages 6023–6032, 2019.
  * [59] Samuel G Müller and Frank Hutter.  Trivialaugment: Tuning-free yet state-of-the-art data augmentation.  In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pages 774–782, 2021.
  * [60] Zhun Zhong, Liang Zheng, Guoliang Kang, Shaozi Li, and Yi Yang.  Random erasing data augmentation.  In Proceedings of the AAAI conference on artificial intelligence (AAAI), volume 34, pages 13001–13008, 2020.
  * [61] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna.  Rethinking the inception architecture for computer vision.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition (CVPR), pages 2818–2826, 2016.
  * [62] Qingyan Meng, Mingqing Xiao, Shen Yan, Yisen Wang, Zhouchen Lin, and Zhi-Quan Luo.  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 (CVPR), pages 12444–12453, 2022.
  * [63] Nitin Rathi and Kaushik Roy.  Diet-snn: A low-latency spiking neural network with direct input encoding and leakage and threshold optimization.  IEEE Transactions on Neural Networks and Learning Systems, 2021\.
  * [64] Nitin Rathi, Gopalakrishnan Srinivasan, Priyadarshini Panda, and Kaushik Roy.  Enabling deep spiking neural networks with hybrid conversion and spike timing dependent backpropagation.  In International Conference on Learning Representations (ICLR), 2020\.

##  6 Supplementary Material for Parallel Spiking Neurons with High Efficiency
and Ability to Learn Long-term Dependencies

###  6.1 Environments for Simulation Speed Benchmark

The simulation speed benchmark is carried out on a Ubuntu 20.04.6 LTS server
with the Intel Xeon(R) Gold 6226R CPU, 192G memory, and the NVIDIA Quadro RTX
6000 (24G) GPU. The PyTorch vision is 1.12.1, the CUDA version is 11.3.1, and
the GPU driver version is 525.105.17.

###  6.2 Parameters Initialization

For the PSN and the masked PSN, the weights are initialized by the kaiming
uniform [54]. More specifically, the values are initialized from
$\mathcal{U}(-\sqrt{5},\sqrt{5})$. For the sliding PSN, the weight is
initialized by an exponential decay method

$\displaystyle W_{i}=2^{i-k+1},~{}~{}~{}~{}i=0,1,...,k-1$ (17)

which is similar to the weight of the LIF neuron with $\tau_{m}=2$ and
removing reset. The thresholds of the PSN family are all initialized as $1$.

###  6.3 Implementation of the Sliding PSN

As mentioned in the main text, the sliding PSN can be implemented by both the
matrix-matrix multiplication and the 1D convolution. However, the 1D
convolution is much slower than the matrix-matrix multiplication, which is
caused by:

1) The parallelism of convolution is lower than matrix-matrix multiplication,
as we can find that the $1\times 1$ convolution is often implemented by the
fully connected layer in modern deep ANNs, such as ConvNeXt [55].

2) The 1D convolution of the sliding PSN requires the input with a shape of
$(...,1,T)$, while the default data format in SpikingJelly is $(T,N,...)$.
Thus, the transpose operation has to be used to move the time-step dimension,
which copies data and slows down the speed.

###  6.4 Training Hyper-parameters

The main hyper-parameters for different datasets are shown in Tab.3. Unless
otherwise specified, the weight decay (wd) is 0, the momentum is 0.9 for the
SGD optimizer, the learning rate scheduler is the cosine annealing schedule
[56], the automatic mixed precision training is used, and the surrogate
function is the arctan surrogate function
$\sigma(x)=\frac{\alpha}{2(1+(\frac{\pi}{2}\alpha x)^{2})}$ with $\alpha=4$
from [14]. Other training options are listed as follows.

Dataset | Optimizer | Batch Size | Epoch | Learning Rate | GPUs | Loss
---|---|---|---|---|---|---
Sequential CIFAR10/100 |  | SGD
---
AdamW(sliding PSN)
128 | 256 |  | 0.1
---
0.001(sliding PSN)
1 | CE
CIFAR10 | SGD | 128 | 1024 | 0.1 | 1 | CE
ImageNet | SGD |  | 64 (SEW-18)
---
32 (SEW-34)
320 | 0.1 | 8 | TET
CIFAR10-DVS | SGD, wd=5e-4 | 32 | 200 | 0.1 | 2 | TET

Table 3: Training hyper-parameters for different datasets.

Static/Sequential CIFAR The transforms include random mixup [57] with
$p=1,\alpha=0.2$, random cutmix [58] with $p=1,\alpha=1$, random choice of two
mix methods with $p=0.5$, random horizontal flip with $p=0.5$, trivial augment
[59], normalization, random erasing [60] with $p=0.1$, and label smoothing
[61] with the amount $0.1$. The number of channels is 256 for static CIFAR and
128 for sequential CIFAR.

ImageNet Transforms are identical to [39]. We loaded the pre-trained weights
from the standard ResNet provided by the ANN community as better initialized
parameters. Similar techniques have been widely used in surrogate training
[62, 63, 64].

CIFAR10-DVS Transforms are identical to [32]. We used the almost identical
network structure, training hyper-parameters, and options as [32], except that
we use the plain BN rather than TEBN.

###  6.5 Memory Consumption of the PSN family

The memory complexity of the SNN during training can be approximated as
$\mathcal{O}(\mathcal{W}+T\cdot(\mathcal{X}+\mathcal{H}))$, where
$\mathcal{W}$ is the number of synapses such as convolutional/linear layers
and weights inside neurons such as $\bm{W}$ in the PSN, $T$ is the number of
time-steps, $\mathcal{X}$ is the input/output of all layers at a single time-
step, and $\mathcal{H}$ is the hidden state of all layers at a single time-
step. When comparing the SNN using the vanilla spiking neurons or the PSN
family, the difference in memory consumption mainly depends on $\mathcal{H}$.

For the vanilla spiking neurons, it can be found that the number of hidden
states is at least $2$, which are $\bm{H}$ and $\bm{V}$. The number of hidden
states increases if the spiking neuron integrates more complex neuronal
dynamics such as adaptative thresholds. While for the PSN, only $\bm{H}$ is
used and the number of hidden states is $1$. Thus, it can be found that the
memory consumption of the PSN is $\mathcal{O}(T\cdot N)$ less than the vanilla
spiking neurons, where $N$ is the number of neurons in the SNN.

We design an interesting experiment to verify our analysis. We measure
$\mathcal{M}_{NO}$, which is the memory consumption of training the network
with only synapses but no neurons to estimate
$\mathcal{O}(\mathcal{W}+T\cdot\mathcal{X}$. Then we measure
$\mathcal{M}_{IF}$ and $\mathcal{M}_{PSN}$, which are SNNs with the simplest
vanilla spiking neuron, the IF neuron, and the PSN to estimate
$\mathcal{O}(\mathcal{W}+T\cdot(\mathcal{X}+\mathcal{H}))$. Thus, we can
estimate the memory consumption of hidden states in two SNNs by $\Delta_{IF-
NO}=\mathcal{M}_{IF}-\mathcal{M}_{NO}$ and $\Delta_{PSN-
NO}=\mathcal{M}_{PSN}-\mathcal{M}_{NO}$. We test on the VGG-11 structure, and
the results are shown in Tab.4. It can be found that the $\Delta_{PSN-
NO}<\Delta_{IF-NO}$, and $\frac{\Delta_{IF-NO}}{\Delta_{PSN-NO}}\approx 2$ in
all cases, which matches our analysis. More specifically, the IF neuron uses
$2$ hidden states, and the PSN uses $1$ hidden state, then the ratio of their
extra memory consumption should be $2$. It can also be found that
$\frac{\Delta_{IF-NO}-\Delta_{PSN-NO}}{T\cdot N}$ is approximately a constant,
which also corresponds to our analysis that the memory consumption of the PSN
is $\mathcal{O}(T\cdot N)$ less than the IF neuron.

$T$ | $N$ | $\mathcal{M}_{NO}$ | $\mathcal{M}_{IF}$ | $\mathcal{M}_{PSN}$ | $\Delta_{IF-NO}$ | $\Delta_{PSN-NO}$ | $\frac{\Delta_{IF-NO}}{\Delta_{PSN-NO}}$ | $\frac{\Delta_{IF-NO}-\Delta_{PSN-NO}}{T\cdot N}$
---|---|---|---|---|---|---|---|---
16 | 16 | 35129.9 | 70563.9 | 53675.9 | 35434 | 18546 | 1.9 | 66.0
8 | 16 | 18935.9 | 36863.9 | 28225.9 | 17928 | 9290 | 1.9 | 67.5
16 | 8 | 18933.9 | 36099.9 | 28223.9 | 17166 | 9290 | 1.8 | 61.5

Table 4: Comparison of the memory consumption of training SNNs with the IF
neuron and the PSN.