PRF: Parallel Resonate and Fire Neuron for Long Sequence Learning in Spiking Neural Networks

## Abstract

Abstract Recently, there is growing demand for effective and efficient long sequence modeling, with State Space Models (SSMs) proving to be effective for long sequence tasks. To further reduce energy consumption, SSMs can be adapted to Spiking Neural Networks (SNNs) using spiking functions. However, current spiking-formalized SSMs approaches still rely on float-point matrix-vector multiplication during inference, undermining SNNs’ energy advantage. In this work, we address the efficiency and performance challenges of long sequence learning in SNNs simultaneously. First, we propose a decoupled reset method for parallel spiking neuron training, reducing the typical Leaky Integrate-and-Fire (LIF) model’s training time from O ⁢ ( L 2 ) 𝑂 superscript 𝐿 2 O(L^{2}) italic_O ( italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) to O ⁢ ( L ⁢ log ⁡ L ) 𝑂 𝐿 𝐿 O(L\log L) italic_O ( italic_L roman_log italic_L ) , effectively speeding up the training by 6.57 × 6.57\times 6.57 × to 16.50 × 16.50\times 16.50 × on sequence lengths 1 , 024 1 024 1,024 1 , 024 to 32 , 768 32 768 32,768 32 , 768 . To our best knowledge, this is the first time that parallel computation with a reset mechanism is implemented achieving equivalence to its sequential counterpart. Secondly, to capture long-range dependencies, we propose a Parallel Resonate and Fire (PRF) neuron, which leverages an oscillating membrane potential driven by a resonate mechanism from a differentiable reset function in the complex domain. The PRF enables efficient long sequence learning while maintaining parallel training. Finally, we demonstrate that the proposed spike-driven architecture using PRF achieves performance comparable to Structured SSMs (S4), with two orders of magnitude reduction in energy consumption, outperforming Transformer on Long Range Arena tasks. 1 1 1 The GitHub repository will be released after paper accepted.

## 1 Introduction

Long Sequence Modeling. Long sequence modeling is a fundamental problem in machine learning. This significant advancement has yielded wide-ranging impact across various fields (such as Transformer Vaswani et al. () and Mamba Gu & Dao ()), including reinforcement learning (e.g., robotics and autonomous driving) Chen et al. (), autoregressive task (e.g., large language models) Zhao et al. (), generative tasks (e.g., diffusion model) Peebles & Xie () and etc. The Transformer architecture Vaswani et al. (), which combines token mixing with self-attention and channel mixing with dense matrices, has been successfully applied in these fields, since sequence data (with length $L$) is well modeled by the self-attention mechanism.

State space models (SSMs) Gu et al. () effectively address the limitations of self-attention machanism in processing long sequences by reducing computational complexity from $O(L^{2})$ to $O(L)$ during inference Feng et al. (). The recurrent form of SSMs allows themselves to scale to longer sequence lengths more efficiently. Furthermore, the (Structured) SSMs utilize orthogonal polynomial bases to initialize recurrent weights Gu et al. () for token mixing, and to compress input history Gu et al. () into hidden state that enables long sequence learning yielding superior performance than Transformer models on long sequence tasks Gu et al. ().
As a result, variants of SSMs (e.g. Mamba Gu & Dao ()) are successfully applied across various fields. However, SSMs still requires large number of float-point matrix-vector multiplications for both token mixing and channel mixing to effectively extract information. These dense matrix multiplication computations scale quadratically with model size $(D)$, consuming significant amount of energy during inference.

Spiking Neural Networks. SNNs benefit from the ability to convert the Float-Point (FP) Multiply Accumulate (MAC) computations into sparse Accumulate (AC) computations through their spike-driven mechanism Hu et al. (); Yao et al. (). Prior arts incorporated the spiking function at the input of token and channel mixing in SSMs to reduce energy consumption in channel mixing Stan & Rhodes (); Bal & Sengupta (); Shen et al. (). However, the token mixing operation within the SSMs blocks has yet been optimized, resulting in extensive FP matrix-vector multiplications and thus compromising the overall energy efficiency Stan & Rhodes (). To avoid FP matrix-vector multiplication computations in both channel mixing and token mixing, it is also crucial to reduce the complexity of token mixing. Therefore, we must revisit and design novel the spiking neurons in SNNs to enabl effective long sequence learning capabilities, while maintaining computation efficiency.

The Challenges for SNNs. There are two challenges in implementing long sequence learning in SNNs. (i) First, the Backpropagation Through Time (BPTT) is the commonly used training method for SNNs Wu et al. (). However, BPTT leads to a quadratic growth in training time with respect to sequence length, scaling as $O(L^{2})$ Kag & Saligrama (). While there are some existing SNN efficient training methods that can improve training efficiency Bellec et al. (); Xiao et al. (); Yin et al. (), they are still outperformed by BPTT for longer sequences Meng et al. (). (ii) Secondly, commonly used spiking neuron models struggle to capture long-range dependencies. Specifically, the widely used LIF neuron has difficulty in capturing and distinguishing dependencies in membrane potential over long intervals, which limits its performance on long-range tasks.Previous work attempted to address the training efficiency and performance improvement separately by improving neuron models Fang et al. (); Spieler et al. (). Fang et al. (); Spieler et al. (). However, efforts made to tackle these challenges on long sequence tasks at the same time remains unveiled. In this work, we aim for addressing these two challenges simultaneously. Our contributions are summarized as follows:

- •
To accelerate the training process, we propose a novel decoupled reset method to implement parallel training that is equivalent to sequential training. This approach accelerates the back propagation by three orders of magnitude and can be applied to any types of spiking neurons.
- •
To effectively extract long range dependencies, we propose the Parallel Resonate and Fire (PRF) neuron with oscillating membrane potential in complex domain leveraging an adaptive and differentiable reset mechanism.
- •
To minimize inference energy, we further incorporate PRF into the design of Spike Driven Temporal and Channel Mixer (SD-TCM) module. This module achieves performance comparable to S4 in long range arena tasks while reducing the inference energy by two orders of magnitude.

## 2 Related Work

##### State Space Models

A general discrete form of SSMs is given by the equation: $\bm{u}_{t}=\bm{A}\bm{u}_{t-1}+\bm{B}\bm{x}_{t}$, $\bm{y}_{t}=\bm{C}\bm{u}_{t}$, where $\bm{A}\in\mathbb{R}^{H\times H}$, $\bm{B}\in\mathbb{R}^{D\times H}$, $\bm{C}\in\mathbb{R}^{H\times D}$ matrices is shape with the model size $D$ and hidden size $H$. The success of the Structured SSMs (S4) arises from the fact that the coefficients of the orthogonal bases are solved to fit an arbitrary sequence curve Gu et al. (). By leveraging orthogonal polynomial bases for initializing structured matrices $\bm{A}$ and $\bm{B}$, S4 effectively compresses input history and outperforms transformers in long-range sequence tasks Gu et al. (). Subsequent variants of SSMs have also achieved great success on long sequence task Gu et al. (); Goel et al. (); Gu et al. (; ); Orvieto et al. (). However, these SSMs-based model still exist power consumption issues that scale quadratically with model size $D$. This is largely due to the use of the dense matrices for channel mixing after $\bm{y}_{t}$, such as in the GLU block Dauphin et al. (), which requires $2D^{2}$ computation in matrices, leading to a significant number of FP-MAC operations.

##### Spikinglized SSMs

The spike mechanism can alleviate the energy problem in dense matrices by converting the FP-MAC as sparse FP-AC computation. Some spiking-formalized (spikinglized) approaches integrate the spike function into SSMs after token mixing to reduce the FP-MAC computation of channel mixing. For instance, Oliver et al.Stan & Rhodes () make intersection of SNNs with S4D model Gu et al. () for long-range sequence modelling, by adding Heaviside function to token mixing output, $\bm{y}_{t}$, at each SSMs layer. Similarly, Abhronil et al.Bal & Sengupta () combine stochastic spiking function at the output of $\bm{y}_{t}$. These spikinglized method can harvest the long sequence learning capability of SSMs models, but also retain nonlinear activation computation and FP matrix-vector multiplications, as they retain the $\bm{A}$, $\bm{B}$ and $\bm{C}$ matrices during recurrent inference. However, this retention limits the energy efficiency advantages of SNNs and poses challenges for deploying the model on neuromorphic chips. Therefore, we aim to further optimize the inference process by reducing matrices $\bm{A}$ and $\bm{B}$ to vectors and eliminating $\bm{C}$. This requires rethinking the role of spiking neurons in long sequence learning.

## 3 Problem Formulation

In this section, we first introduce the commonly used Leaky Integrate-and-Fire (LIF) neuron, then we describe the two main challenge for long range learning ability with spiking neuron.

### 3.1 The Leaky Integrate-and-Fire (LIF) Neuron

The LIF model is a widely used spiking neuron model. It simulates neurons by integrating input signals and firing spikes when the membrane potential exceeds a threshold. The dynamic of membrane potential $u(t)$ is followed by:

$$ $\frac{du(t)}{dt}=-\frac{1}{\tau}\left(u(t)-u_{\mathrm{reset}}\right)+\frac{R}{ \tau}c(t),$ (1) $$

where the $c(t)$ is the input current. The constants $\tau$, $u_{\mathrm{reset}}$, and $R$ denote the membrane time constant, reset potential, and resistance. We use $R=\tau$ as in previous work. When $u(t)$ reaches the threshold $V_{\mathrm{th}}$, the neuron fires a spike, and then $u(t)$ resets to $u_{\mathrm{reset}}$. The discrete LIF model is expressed as:

$$ $\displaystyle u_{t}$ $\displaystyle=\beta\cdot(u_{t-1}-V_{\mathrm{th}}s_{t-1})+c_{t},$ (2) $\displaystyle s_{t}$ $\displaystyle=\mathcal{H}(u_{t}-V_{\mathrm{th}})=\begin{cases}1,&\text{if }u_{ t}\geq V_{\mathrm{th}}\\ 0,&\text{otherwise}\end{cases},$ (3) $$

where $\beta\triangleq 1-\frac{1}{\tau}\in(0,1)$, and the discrete timesteps $t=1,2,\ldots,T$, the initial situation $s_{0}=u_{0}=0$. After firing, the membrane potential resets according to previous spike $s_{t-1}$.
This spiking neuron face two primary challenges for long sequence tasks: (i) The coupled reset prevents parallel training along timesteps, causing training time significantly for long sequences. (ii) The commonly used LIF neuron model struggles to capture long-range dependencies, hindering the performance on long sequence tasks. The overview as shown in Figure 1.

Figure: Figure 1: (a) The reset mechanism prevents parallelization as the state update relies on previous spike output, causing $O(L^{2})$ timing cost. (b) The proposed decoupled reset mechanism enables parallel computation. (c) Fast decay causes long-range dependencies to vanish in the membrane potential. (d) Slow decay could generate dependency over long sequence, but causes dependency ambiguity. (e) The resonate mechanism with membrane potential helps distinguish relevant inputs.
Refer to caption: x1.png

### 3.2 Problem 1: Coupled Reset Prevents Parallelism along Timesteps

The first problem is that the training cost of LIF-based sequential computation increases exponentially as the number of timesteps extends, with a complexity of $O(L^{2})$. This occurs because BPTT method requires the computational graph to expand along the time dimension Kag & Saligrama (). While the sequential computation of linear combination can be parallelized to reduce the time cost complexity from $O(L^{2})$ to $O(LlogL)$ during training, as done in SSMs. However, the $u_{t}$ dependency on previous spikes output $s_{t-1}$ with nonlinear Heaviside function. By recursively expanding $u_{t}$ and simplifying in Equation 4:

$$ $\displaystyle u_{t}=c_{t}+\beta c_{t-1}+...+\beta^{t-1}c_{1}-V_{\mathrm{th}}( \beta s_{t-1}+\beta^{2}s_{t-2}+...+\beta^{t-1}s_{1}).$ (4) $$

The reset mechanism causes $u_{t}$ to depend on all previous spike outputs. This coupling hinders the neuron’s ability to perform parallel computations, further limiting effective parallel computation for both $u_{t}$ and $s_{t}$, as illustrated in Figure 1 (a). This forces the spiking neuron to compute sequentially. Consequently, as the number of timesteps increases during training, the time cost grows quadratically. To address this issue, we propose the decoupled reset method, which facilitates parallel computation of spiking neurons with the reset mechanism (Illustrated in Figure 1(b) and discussed in Method 4.1).

### 3.3 Problem 2: Commonly Used LIF Struggle with Long-Range Dependencies

The second problem is that the commonly used LIF model struggles to capture and distinguish the dependencies in the membrane potential over long interval. This limitation arises due to the dilemma of decay factor $\beta$ in the membrane potential dynamics. When $\beta$ is small, the membrane potential decreases rapidly. This Fast Decay cause the neuron to quickly forget past inputs (Figure 1 (c)). For example, in Case $1$, if there are two inputs, $c_{t_{0}}$ and $c_{T}$, separated by $T$ time steps, result in $u_{T}$ becoming independent of $c_{t_{0}}$. Similarly, in Case $2$, with zero input $c_{t_{0}}$, the result is identical with Case $1$. This failure to retain long-term information leads to the vanishing of dependencies. Conversely, when $\beta$ is large, the membrane potential retains its value over an extended period $T$ (Figure 1 (d)). This Slow Decay can solve the dependency vanishing problem, but causing ambiguity between closely spaced inputs. For instance, input $c_{t_{0}}=c_{t_{0}^{\prime}}$ with $t_{0}^{\prime}=t_{0}+\delta t$ result in similar membrane potentials, there is small difference between the membrane potential $u_{T}=\beta^{(T)}\cdot c_{t_{0}}+c_{T}$ in Case $1$ and the $u_{T}=\beta^{(T-\delta t)}\cdot c_{t_{0}^{\prime}}+c_{T}$ in Case $2$, as shown in Figure 1 (d). To address this challenge, we propose enhancing the neuron’s dynamics by incorporating resonate mechanism. This allows neurons to remain sensitive to input, even after a long interval of $T$. (Illustrated in Figure 1 (e) and discussed in detail in Method 4.2).

## 4 Method

In this section, we present our approach to address two primary challenges in SNNs: parallel training and long-range dependency learning. For parallel training, we decouple the reset mechanism from the integrate computation to implement parallel training while maintaining equivalence with sequential computation. For long-range dependency learning, we propose the Parallel Resonate and Fire (PRF) Neuron by incorporating the reset process as an imaginary part into the time constant, enabling the neuron to achieve long-range learning ability.

### 4.1 decoupled reset for Parallel Computation

To enable parallel computation, we need to decouple the causal relationship with previous spikes, by separating the linear combination part from the nonlinear causal dependency part. To achieve this, we first substitute the Equation 4 into the Equation 3 by expanding the $u_{t}$ in the Heaviside function. Then we merge the reset part into the threshold $V_{\mathrm{th}}$. As such, Equation 3 is rewritten as Equation 5:

$$ $\displaystyle s_{t}=\mathcal{H}(\underbrace{c_{t}+\beta c_{t-1}+...+\beta^{t-1 }c_{1}}_{\text{Leaky and Integrate}\quad\triangleq\quad u_{t}^{\prime}}- \underbrace{V_{\mathrm{th}}\cdot\left((\beta s_{t-1}+\beta^{2}s_{t-2}+...+ \beta^{t-1}s_{1})+1\right)}_{\text{decoupled reset}\quad\triangleq\quad d_{t}}),$ (5) $$

where the linear combination of leaky and integrate part is defined as $u_{t}^{\prime}$, the second part is defined as the decoupled reset $d_{t}$. Now the spike output is:

$$ $\displaystyle s_{t}$ $\displaystyle=\mathcal{H}\left(u^{\prime}_{t}-d_{t}\right)=\begin{cases}1,& \text{if }u^{\prime}_{t}\geq d_{t}\\ 0,&\text{otherwise}\end{cases},$ (6) $$

where this spike output $s_{t}$ depends on $u^{\prime}_{t}$ and $d_{t}$. Note that the first part is linear combination of $u^{\prime}_{t}=\beta u^{\prime}_{t-1}+c_{t}$, we define the sequence of $\mathbf{U}_{T}^{\prime}$ as the ordered set $\{u^{\prime}_{1},u^{\prime}_{2},...,u^{\prime}_{T}\}$ over all timesteps (detailed notation is described in Appendix A). This linear combination of $u^{\prime}_{t}$ can be computed using convolution and further accelerated by converting the convolution into multiplication after applying the Fast Fourier Transform: $\mathbf{U}_{T}^{\prime}=\mathbf{C}_{T}*\mathbf{K}_{T}=\mathcal{F}^{-1}\left(
\mathcal{F}(\mathbf{C}_{T})\odot\mathcal{F}(\mathbf{K}_{T})\right)$, with the $O(LlogL)$ computation complexity, avoiding the $O(L^{2})$ training time in BPTT.

After efficient parallel computation of $u_{t}^{\prime}$, the $d_{t}$ still remains the dependency relationship with previous spikes $s_{t}$. To further decouple this dependency from the spike output $s_{t}$, the key idea is converting the recursive form to the iterative form. We first define the dependency part as $A_{t}$, as shown in Equation 7. Now, we only need to convert the $A_{t}$ from its recursive form into an iterative form. By separating the last spike from all previous spikes, we obtain Equation 8. The left part of this equation corresponds to Equation 6, while the right part refers to the recursive dependency itself, leading to Equation 9 and Equation 10:

$$ $\displaystyle A_{t}\triangleq$ $\displaystyle\beta s_{t-1}+\beta^{2}s_{t-2}+...+\beta^{t-1}s_{1}$ (7) $\displaystyle=$ $\displaystyle\beta s_{t-1}+\beta\left(\beta s_{t-2}+...+\beta^{t-2}s_{1}\right)$ (8) $\displaystyle=$ $\displaystyle\beta\mathcal{H}\left(u^{\prime}_{t-1}-d_{t-1}\right)+\beta A_{t-1}$ (9) $\displaystyle=$ $\displaystyle\begin{cases}\beta\cdot(1+A_{t-1}),&\text{if }u_{t-1}^{\prime} \geq d_{t-1}\\ \beta\cdot A_{t-1},&\text{if }u_{t-1}^{\prime}<d_{t-1}\end{cases},$ (10) $$

where the sequence of $d_{t}$ is calculated according to the sequence of $u^{\prime}_{t}$, by dynamically updating $A_{t}$:

$$ $d_{t}=V_{\mathrm{th}}\cdot(A_{t}+1),\quad A_{0}=u^{\prime}_{0}=0.$ (11) $$

At this point, all recursive forms have been converted into dynamic equations. The formation of $d_{t}$ is completely independent of $s_{t}$. The decoupled reset function $\mathbf{D}_{T}=f_{D}(\mathbf{U}_{T}^{\prime})$ is deduced as shown in Equation 10 and 11, where $d_{t}$ can be dynamically scanned from all $u_{t}^{\prime}$ with $O(L)$ complexity.

In summary, we convert the all calculation of $s_{t}$ with $O(L^{2})$ complexity into a combination of $u^{\prime}_{t}$ with $O(L\log L)$ complexity and subsequently $d_{t}$ with $O(L)$ complexity, achieving a training speed-up of approximately $L/(\log L+1)$. To summarize, this approach facilitates parallel computation:
$\displaystyle\mathbf{U}_{T}^{\prime}$
$\displaystyle=\mathbf{C}_{T}*\mathbf{K}_{T}$
(12a)
$\displaystyle=\mathcal{F}^{-1}\left(\mathcal{F}(\mathbf{C}_{T})\odot\mathcal{F
}(\mathbf{K}_{T})\right)$
(12b)
$\displaystyle\mathbf{D}_{T}$
$\displaystyle=f_{D}(\mathbf{U}_{T}^{\prime})$
(12c)
$\displaystyle\mathbf{S}_{T}$
$\displaystyle=\mathcal{H}(\mathbf{U}_{T}^{\prime}-\mathbf{D}_{T})$
(12d)
$\displaystyle\mathbf{U}_{T}^{\prime}$
$\displaystyle=(u_{1}^{\prime},u_{2}^{\prime},\ldots,u_{T}^{\prime})$
(13a)
$\displaystyle\mathbf{C}_{T}$
$\displaystyle=(c_{1},c_{2},\ldots,c_{T})$
(13b)
$\displaystyle\mathbf{K}_{T}$
$\displaystyle=(\beta^{0},\beta^{1},\ldots,\beta^{T-1})$
(13c)
$\displaystyle\mathbf{D}_{T}$
$\displaystyle=(d_{1},d_{2},\ldots,d_{T})$
(13d)

Where the outputs $\mathbf{S}_{T}=(s_{1},s_{2},\ldots,s_{T})$ from Equation 12d is equivalent with the sequential generated from Equation 3. The kernel vector $\mathbf{K}_{T}$ with $\beta=1-\frac{1}{\tau}$.By combining the above equations, we obtain the parallel computation process, as shown in Algorithm.1 and 2 in Appendix.B.Although parallelized LIF solves the training problem for long sequences (Experiment 5.1), it still does not perform well on long sequences (Experiment 5.2), due to the common issue of long-range dependencies in LIF models (Problem 3.3).

### 4.2 Resonate for Long-Range Learning Ability

To address the challenge of capturing long-range dependencies in spiking neural networks, we introduce the Parallel Resonate-and-Fire (PRF) neuron. This neuron model extends the standard LIF neuron by incorporating a resonance mechanism into its dynamics, enabling it to retain information over longer periods while maintaining computational efficiency.

Firstly, recall the commonly used LIF model from Equation 1, the dynamics are given by: $\frac{du(t)}{dt}=\gamma u(t)-\gamma u_{\mathrm{reset}}+c(t)$, where $\gamma\triangleq-1/\tau$. To introduce more dynamic behavior, we use $r(t)$ and $\theta$ to replace $u_{\mathrm{reset}}$ and $\gamma$ in the reset part, in respectively. Then the dynamic are rewritren as:

$$ $\frac{du(t)}{dt}=\gamma u(t)-\theta r(t)+c(t).$ (14) $$

To introduce membrane potential oscillations, we define a complex decay constant $\tilde{\gamma}=\gamma+i\theta$, where $i=\sqrt{-1}$. In the vanilla LIF model, the reset $r(t)$ is a non-continuous conditional function with $\theta=\gamma$ in Equation 14, and the reset result is instantaneous process:

$$ $\displaystyle\frac{dr(t)}{dt}=\begin{cases}V_{\mathrm{th}}\delta(t),&\text{if }u(t)\geq V_{\mathrm{th}}\\ 0,&\text{if }u(t)<V_{\mathrm{th}}\end{cases}.$ (15) $$

While this reset has proven effective, it poses challenges for efficient computation. We introduce a reset function that is both effective and computationally simple for both forward and backward propagation. We define the reset function to satisfy the condition: $\tilde{u}(t)\triangleq u(t)+ir(t)$. Substituting this definition into Equation 14, we obtain the PRF model:

$$ $\displaystyle\frac{d\tilde{u}(t)}{dt}$ $\displaystyle=\tilde{\gamma}\tilde{u}(t)+c(t),$ (16) $$

where the real part $\Re\{\tilde{u}(t)\}$ corresponds to the membrane potential dynamics. This formulation incorporates the reset mechanism into the time constant via the imaginary component $\theta$, introducing resonance into the neuron’s behavior. When $\theta=0$, the model reduces to the standard LIF neuron without reset. If we extend the Equation 16, we can be expressed in matrix form:

$$ $\begin{pmatrix}\dot{u}(t)\\ \dot{r}(t)\end{pmatrix}=\begin{pmatrix}\gamma&-\theta\\ \theta&\gamma\end{pmatrix}\begin{pmatrix}u(t)\\ r(t)\end{pmatrix}+\begin{pmatrix}c(t)\\ 0\end{pmatrix},$ (17) $$

which give us the insight that this reset process is continuous function:

$$ $\displaystyle\frac{dr(t)}{dt}=\gamma u(t)+\theta r(t),$ (18) $$

where this formulation treats the reset as a continuous function controlled by $\theta$, allowing for a more gradual and reset process.

We then discretize the model (Equation 16) (details are provided in Appendix C), yielding the sequential and parallel formulations of the PRF neuron, as shown in Equation 19 and 20 in respectively:

$$ $\displaystyle\tilde{u}_{t}$ $\displaystyle=\exp(\Delta\tilde{\gamma})\tilde{u}_{t-1}+\Delta c_{t}$ (19a) $\displaystyle s_{t}$ $\displaystyle=\mathcal{H}\left(\Re\{\tilde{u}_{t}\}-V_{\mathrm{th}}\right)$ (19b) $\displaystyle\tilde{\gamma}$ $\displaystyle=\gamma+i\theta$ (19c) $$

$$ $\displaystyle\mathbf{\tilde{U}}_{T}$ $\displaystyle=\mathcal{F}^{-1}\left(\mathcal{F}\left(\mathbf{C}_{T}\right) \odot\mathcal{F}(\mathbf{\tilde{K}}_{T})\right)$ (20a) $\displaystyle\mathbf{S}_{T}$ $\displaystyle=\mathcal{H}\left(\Re\{\mathbf{\tilde{U}}_{T}\}-V_{\mathrm{th}}\right)$ (20b) $$

Now the membrane potential $\widetilde{u}_{t}\in\mathbb{C}$ ,where $\Delta$ is the time step size. The details of parallel computation is described in Algorithm.3. Where $\mathbf{C}_{T}=(c_{1},c_{2},\ldots,c_{T})$ is the input current sequence, and the kernel $\mathbf{\tilde{K}}_{T}$ is calculated by:

$$ $\mathbf{\tilde{K}}_{T}=\left(\Delta A^{(0)},\Delta A^{(1)},\ldots,\Delta A^{(T -1)}\right),$ (21) $$

with $A=\exp\left(\Delta\tilde{\gamma}\right)$. By incorporating the imaginary component $\theta$, the neuron exhibits oscillatory behavior, allowing it to resonate at specific input frequencies (see Appendix D for details). This oscillation phenomenon can be found in biological neurons Izhikevich (). Our model differs from variants of resonant models, like BHRF Higuchi et al. (), which combine adaptive thresholds and refractory mechanisms. Instead, we incorporate the reset directly into the time constant through the imaginary component while enabling efficient parallel computation for training. Furthermore, the PRF neuron can also be deployed on neuromorphic chips for inference after parallel training, requiring only two additional multiplications and one more addition than LIF (see details in Appendix E).

### 4.3 Theory Analysis

(1) Sequential Perspective on Dynamic. Decoupling the reset can be viewed as transforming the reduction in membrane potential caused by increased threshold value. This means that the soft reset mechanism could be equivalent to the adaptive threshold mechanism Bellec et al. () as deduced in the Theorem 1. Furthermore, the LIF model without a reset can be considered a specific instance of a PRF, as demonstrated in Theorem 2. Moreover, the $u_{t}$ in PRF will converge as shown in Theorem 3. This indicates that the membrane potential of PRF is stable and bounded.

###### Theorem 1 .

Let $V_{\mathrm{th}}=1$ and $\rho=1$ in Adaptive-LIF model, then the LIF neuron with soft reset model is equivalent to the Adaptive-LIF without reset mechanism. (The proof See Appendix.F.1)

###### Theorem 2 .

Let $\Delta=1$ and $\theta=0$, then the PRF model could degenerate as Valina LIF model without reset mechanism. (The proof see Appendix.F.2)

###### Theorem 3 .

If the inputs $c_{t}\sim\mathcal{N}(0,\sigma^{2})$ follow a normal distribution, then the membrane potential of PRF will converge as a distribution $u_{t}\sim\mathcal{N}(0,\frac{\tau\Delta}{2}\sigma^{2})$, as $t\to+\infty$. (The proof see Appendix.F.3)

(2) Parallel Perspective on Gradients. The issue of membrane potential dependence described in Problem 3.3 is fundamentally a problem of gradients, which are correlated with the kernel and previous spikes. Assuming the input current is $\mathbf{C}^{l}_{T}=\bm{W}^{l}\mathbf{S}^{l-1}_{T}$ across all $T$, the gradients are proportional to:

$$ $\displaystyle\nabla_{\bm{W}^{l}}\mathcal{L}\propto\underbrace{\frac{\partial \mathcal{L}}{\partial u^{l}_{T}}\sum_{t=1}^{T}\frac{\partial u^{l}_{T}}{ \partial u^{l}_{t}}\frac{\partial u^{l}_{t}}{\partial\bm{W}^{l}}}_{\text{ Sequential Perspective}}\propto\underbrace{\sum_{t=1}^{T}\beta^{(t)}s^{l-1}_{t }=\left<\mathbf{K}_{T},\mathbf{S}^{l-1}_{T:1}\right>,}_{\text{Parallel Perspective}}$ (22) $$

here $\left<\cdot,\cdot\right>$ denotes the inner product, $\mathbf{S}_{T}^{l-1}=(s_{1}^{l-1},s_{2}^{l-1},\ldots,s_{T}^{l-1})$ is the sequence of spike outputs from the previous layer $l-1$, and $\mathbf{K}_{T}$ represents the parallelism kernel. The subscript ($T:1$) indicates the sequence reversal. A small $\beta$ causes gradient vanishing due to a narrow receptive field, making it likely for sparse spikes to fall outside this field and yield an inner product close to zero. Conversely, a large $\beta$ results in a long-range, slow-decaying kernel, leading to overly consistent gradient values across spike positions. However, an oscillating kernel with a large $\beta$ can adjust the gradient at different spike positions, smoothing the gradients and improving the neuron’s representational capability (see Appendix G for details). Experiment 5.2 verifies this insight.

###### Theorem 1 .

###### Theorem 2 .

###### Theorem 3 .

### 4.4 Architecture

The Spike-Driven Token and Channel Mixer (SD-TCM) Module (Figure 2, with further details in Appendix H) is inspired by token and channel mixing in Transformer and S4 architectures. For token mixing, we use a PRN Neuron followed by a Linear layer, while for channel mixing, the Neuron is replaced by a Spatial Neuron. The Spatial Neuron, a variant of the LIF Neuron, focuses on instantaneous information at each timestep by setting the time constant $\tau$ close to 1: $\lim_{\tau\to 1^{+}}u_{t}=\left(1-\frac{1}{\tau}\right)u_{t-1}+c_{t}\approx c_
{t}$. This simplifies the output to $s_{t}=\mathcal{H}(c_{t}-V_{\mathrm{th}})$, resembling motor neurons in biology with rapid decay.

During training, a trainable amplitude $\alpha$ acts as a gate, and during inference, $\alpha$ can be merged into the Linear layer: $(\alpha\bm{s}_{t})\times\bm{W}\equiv\bm{s}_{t}\times(\alpha\bm{W})$. Membrane shortcut residual connections are used to maintain event-driven, spike-based communication. As shown in Table 1, the PRF requires only $O(5D)$ computational complexity, while SSMs and Spikinglized SSMs require $O(H^{2}+2DH)$. This makes the architecture rely only on FP-AC and element-wise multiplication, reducing energy consumption and simplifying deployment on neuromorphic chips. As a result, it achieves lower computational complexity and energy consumption (see analysis in Appendix I).

Figure: Figure 2: Diagram of the SD-TCM.
Refer to caption: x2.png

## 5 Experiments

To evaluate the efficiency and effectiveness of the parallel method, as well as the performance improvement on long sequence tasks, we first demonstrate that parallel significantly accelerates training while maintaining equivalence with sequential computation. Next, we explore the PRF neuron’s ability to handle long-range dependencies, showing that kernel oscillations improve both performance and gradient stability. Finally, the SD-TCM module achieves performance comparable to S4 while reducing energy consumption by over 98.57% on Long Range Arena tasks. Detailed experimental setups and training hyperparameters are provided in Appendix J.

### 5.1 The Parallel Process

First, we compare the training runtime across different timesteps in Figure LABEL:fig.seq.par.train.time (left) using three repeated experiments. Beyond 4 timesteps, parallel training consistently outperforms sequential training. For example, at 1,024 timesteps, sequential training takes 4.6 seconds per iteration, while parallel training takes only 0.7 seconds, achieving a speedup of $6.57\times$. This acceleration becomes more significant as the number of timesteps increases, with a speedup of $9.35\times$ at 16,384 timesteps and $16.50\times$ at 32,768 timesteps. The equivalence between sequential and parallel computation is maintained during inference and training, with only a minor accuracy difference during training (details in Appendix K), which we tentatively attribute to numerical error.

The speedup is primarily due to parallel training avoiding the recursive unfolding of the computational graph. The forward and backward passes are accelerated by $1.77\times$ and $132.55\times$, respectively, as shown in Figure LABEL:fig.seq.par.train.time (right). The significant speedup in the backward pass occurs because sequential training requires unfolding the graph at each timestep, which is time-consuming, whereas parallel training computes the graph only once. Further details of the timing for both training and inference can be found in Appendix L.

### 5.2 The Long Range Learning Ability

Although parallelized LIF can speed up training, it still struggles with performance on simple long-sequence tasks, such as sequential-MNIST, as shown in Figure LABEL:fig.long.range.dependency.perspective (left). However, after introducing oscillations in the kernel, the PRF neuron successfully solves the sequential-MNIST problem, achieving both training efficiency and effectiveness.

To better understand how oscillating membrane potential improves performance, we compare the performance of various $\beta$ values, with and without oscillations, on the more challenging permuted-sMNIST dataset. The results after training for 50 epochs are summarized in Figure LABEL:fig.long.range.dependency.perspective (right). (We fit $\Delta=1$ and set $\theta=\frac{\pi}{2}$ while varying the $\beta$ hyperparameter without training.) Without the oscillating term, as $\beta$ increases, accuracy initially improves as expected but then decreases beyond a certain point. In contrast, the introduction of oscillations helps counteract this decline and further enhances performance, indicating that oscillations can improve performance for long sequences.

We further examine the loss landscape contours Li et al. () for each case in Figure LABEL:fig.loss.landscape, corresponding to Figure LABEL:fig.long.range.dependency.perspective (right). Gradients are perpendicular to the contour lines, with sparse or dense contours indicating smaller or larger gradients, respectively. Extremely sparse or dense contours suggest vanishing or exploding gradients. In Figure LABEL:fig.loss.landscape (a), fast decay with a small $\beta$ shows gradient vanishing, with extremely sparse contours. Increasing $\beta$ with moderate decay can alleviate gradient vanishing, but the model still encounters a local optima region (Figure LABEL:fig.loss.landscape (b)). Further increasing $\beta$ with slow decay is expected to fully address the vanishing problem but introduces gradient ambiguity with overlapping contours (Figure LABEL:fig.loss.landscape (c)). In contrast, the introduced oscillation from resonate results in smoother gradients (Figure LABEL:fig.loss.landscape (d)), enabling more effective feature extraction.

As shown in Table 2, the PRF neuron uses only $68.9k$ parameters ($768$ for neurons and $67.2k$ for synapses) to achieve state-of-the-art (SOTA) results while maintaining training efficiency. Previous models required recurrent connections in the linear layers to solve the (p)s-MNIST task. To align the training parameters, we modified the architecture to a feedforward linear layer with a 1-(128)^3-10 structure. Using parallel computation, we completed the entire training process in only $1.60$ hours with a batch size of 256 for 200 epochs.

**Table 2: The neuron for (permute)-sequential MNIST task. The pixel-level image classification captures the hierarchical structure to verify the long-range capture ability. The symbol $*$ denotes including the feedback linear connection, while the symbol of $\uparrow$ and $\downarrow$ indicate that larger and smaller values are better, respectively. The BHRF result is referenced from Higuchi et al. (), while other results are from Zhang et al. ().**
| Task<br>(Length) | Spiking Neuron | Seq.<br>Infer. | Par.<br>Train. | No. Params.<br>(k) $\downarrow$ | Top-1 Test Acc.<br>(%) $\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| sMNIST<br>(784) | LIF | ✓ | ✗ | 85.1 / 155.1* | 72.06 / 89.28 |
| PLIF Fang et al. ()${}^{\text{2021' ICCV}}$ | ✓ | ✗ | 85.1 / 155.1* | 87.92 / 91.79 |  |
| GLIF Yao et al. ()${}^{\text{2022' NeurIPS}}$ | ✓ | ✗ | 87.5 / 157.5* | 95.27 / 96.64 |  |
| TC-LIF Zhang et al. ()${}^{\text{2024' AAAI}}$ | ✓ | ✗ | 85.1 / 155.1* | 97.35 / 99.20 |  |
| ALIF Yin et al. ()${}^{\text{2021' Nat. MI}}$ | ✓ | ✗ | 156.3* | 98.70 |  |
| AHP Rao et al. ()${}^{\text{2022' Nat. MI}}$ | ✓ | ✗ | 68.4* | 96.00 |  |
| BHRF Higuchi et al. ()${}^{\text{2024' ICML}}$ | ✓ | ✗ | 68.9 | 99.1 ± 0.1 |  |
| PRF (Ours) | ✓ | ✓ | 68.9 | 99.18 |  |
| psMNIST<br>(784) | LIF | ✓ | ✗ | 155.1* | 80.36 |
| GLIF Yao et al. ()${}^{\text{2022' NeurIPS}}$ | ✓ | ✗ | 157.5* | 90.47 |  |
| TC-LIF Zhang et al. ()${}^{\text{2024' AAAI}}$ | ✓ | ✗ | 155.1* | 95.36 |  |
| ALIF Yin et al. (; )${}^{\text{2021' Nat. MI}}$ | ✓ | ✗ | 156.3* | 94.30 |  |
| BHRF Higuchi et al. ()${}^{\text{2024' ICML}}$ | ✓ | ✗ | 68.9 | 95.2 |  |
| PRF (Ours) | ✓ | ✓ | 68.9 | 96.87 |  |

### 5.3 Long Range Arena Tasks

To demonstrate the long-range dependency analysis capability of our SD-TCM module, we evaluate it using the Long Range Arena (LRA) benchmark Tay et al. (). This benchmark covers a wide range of classification tasks, including both textual and image domains. For the ListOps, Text, and Retrieval tasks, we use the causal architecture, while S4 employs a bidirectional architecture for all tasks. For the Image and Pathfinder tasks, we use a bidirectional architecture.

**Table 3: Comparison of Accuracy, Parameters and Energy.**
| Metric | Model | ListOps | Text | Retrieval | Image | Pathfinder | Avg. |
| --- | --- | --- | --- | --- | --- | --- | --- |
| En.(mJ) | S4<br>Ours | 5.104<br>0.075 | 3.718<br>0.298 | 24.439<br>0.211 | 19.222<br>0.187 | 6.256<br>0.067 | 11.748<br>0.168 |
| Acc.(%) | S4<br>Ours | 59.60<br>59.20 | 86.82<br>86.33 | 90.90<br>89.88 | 88.65<br>84.77 | 94.20<br>91.76 | 84.03<br>82.39 |
| Par. | S4<br>Ours | 815 k<br>272 k | 843 k<br>830 k | 3.6 M<br>1.1 M | 3.6 M<br>4.1 M | 1.3 M<br>1.3 M | -<br>- |

The SD-TCM module achieves performance comparable to S4 while reducing energy consumption by over 98.57% on average, as shown in Table 4. Specifically, the accuracy on ListOps is 59.60%, compared to S4’s 59.20%. The key advantage is the significant reduction in energy consumption, as shown in Table 4, with ListOps dropping from 5.104 mJ to 0.075 mJ. This energy reduction mainly benefits from the extreme sparsity of spikes, with detailed firing rate statistics provided in Appendix M. Additionally, the model is sensitive to the imaginary part of $\theta$, causing fluctuations with different initialization (see Appendix N for details). Table 4 shows an ablation study on the $\alpha$ parameter in the Text (4096) task. Applying $\alpha$ to both the spatial ($\mathcal{SN}$) and PRF ($\mathcal{TN}$) components slightly lowers accuracy from 85.75% to 85.69%, but applying it only to $\mathcal{SN}$ improves accuracy to 86.33%, demonstrating its benefit for spatial dependencies. Finally, as shown in Table 5, our module achieves performance comparable to the S4 baseline across tasks while preserving spike-driven feature, avoiding nonlinear activation functions and FP MAC operations.

**Table 5: Test Accuracy Comparison on LRA Tasks (%) $(\uparrow)$. ’NL Act. Free’ and ’FP MAC Free’ denote models that do not use nonlinear activation functions or floating-point multiply-accumulate operations in the block, respectively. The underline and bold formatting indicate the SoTA result for Spikinglized SSMs and Improving Neuron methods, respectively.**
| Model | NL Act. | FP MAC | ListOps | Text | Retrieval | Image | Pathfinder | Avg. |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| (Input length) | -Free | -Free | (2,048) | (4,096) | (4,000) | (1,024) | (1,024) |  |
| Random (Lower Bound) | - | - | 10.00 | 50.00 | 50.00 | 10.00 | 50.00 | 34.00 |
| Transformer Vaswani et al. () | ✗ | ✗ | 36.37 | 64.27 | 57.46 | 42.44 | 71.40 | 54.39 |
| S4 (Bidirectional) Gu et al. () | ✗ | ✗ | 59.60 | 86.82 | 90.90 | 88.65 | 94.20 | 84.03 |
| Binary S4D Stan & Rhodes () ${}^{\text{2024' Sci. Rep.}}$ | ✗ | ✗ | 54.80 | 82.50 | 85.03 | 82.00 | 82.60 | 77.39 |
| $\hookrightarrow+$ GSU & GeLU | ✗ | ✗ | 59.60 | 86.50 | 90.22 | 85.00 | 91.30 | 82.52 |
| Stoch. SpikingS4 Bal & Sengupta () ${}^{\text{2024' arXiv}}$ | ✗ | ✗ | 55.70 | 77.62 | 88.48 | 80.10 | 83.41 | 77.06 |
| SpikingSSMs Shen et al. () ${}^{\text{2024' arXiv}}$ | ✗ | ✗ | 60.23 | 80.41 | 88.77 | 88.21 | 93.51 | 82.23 |
| Spiking LMU Liu et al. () ${}^{\text{2024' ICLR}}$ | ✓ | ✗ | 37.30 | 65.80 | 79.76 | 55.65 | 72.68 | 62.23 |
| ELM Neuron Spieler et al. () ${}^{\text{2024' ICLR}}$ | ✗ | ✗ | 44.55 | 75.40 | 84.93 | 49.62 | 71.15 | 69.25 |
| Spike-Driven TCM | ✓ | ✓ | 59.20 | 86.33 | 89.88 | 84.77 | 91.76 | 82.39 |

## 6 Conclusion

This study aims to solve the SNNs problem of parallelization and performance on long sequences. We propose the decoupled reset method, enabling spiking neurons could parallel training. This method can be applied to any type of neuron to speed up. Additionally, we introduce the PRF neuron, incorporating the reset as an imaginary part to formulate oscillations in the membrane potential, which solves the long-range dependency problem. The SD-TCM model, combined with PRF neurons, achieves performance comparable to S4 on the LRA task while reducing energy consumption by two orders of magnitude. However, due to the sensitivity of training to neuron initialization, the PathX problem remains unsolved. This issue could be solved in the future by using better hyperparameters and initialization strategies.

## Appendix A Notation in the Paper

Throughout this paper and in this Appendix, we use the following notations. Matrices are represented by bold italic capital letters, such as $\bm{W}$, while sequences are represented by bold non-italic capital letters, such as $\mathbf{X}_{T}=\{x_{1},x_{2},...,x_{T}\}$. For a function $f(\bm{x}):\mathbb{R}^{d_{1}}\rightarrow\mathbb{R}^{d_{2}}$, we use $\nabla{\bm{x}}f$ instead of $\frac{\partial f}{\partial\bm{x}}$ to denote the first-order derivative of $f$ with respect to $\bm{x}$. The symbols $\odot$ and $\langle\cdot,\cdot\rangle$ represent the element-wise product and the inner product, respectively.

## Appendix B The Algorithm of Pseudo-Code

Algorithms 1 and 2 describe the parallel computation of the LIF model, while Algorithm 3 outlines the computation process for the PRF model.

Figure: Algorithm 1 Parallel Computation of LIF

Figure: Algorithm 2 $f_{D}(U)$ decoupled reset

Figure: Algorithm 3 Parallel Computation of PRF Model

## Appendix C The Discretization of PRF Neurons

Firs, we recall the PRF model:

$$ $\displaystyle\frac{d\tilde{u}(t)}{dt}$ $\displaystyle=\tilde{\gamma}\tilde{u}(t)+c(t),$ (23) $$

where complex membrane potential $\tilde{u}(t)=u(t)+ir(t)$ and a complex decay constant $\tilde{\gamma}=\gamma+i\theta$ with $i=\sqrt{-1}$.

Note that here $c(t)$ is treated as a fixed external current input, which is constant from the point of view of this ODE in $u_{i}(t)$. Let $c_{k}$ denote the average value in each discrete time interval. Assumes the value of a sample of u is held constant for a duration of one sample interval $\delta$.

$$ $c_{t_{k}}=\frac{1}{\Delta t_{k}}\int_{t_{k-1}}^{t_{k}}c(t)dt$ (24) $$

Since we only observe the real part of the membrane potential, the input is a real number. Thus, we can consider the real part $\gamma=-\frac{1}{\tau}$ to describe the model. The model can then be expressed as:

$$ $\displaystyle\frac{du(t)}{dt}$ $\displaystyle=-\frac{1}{\tau}u(t)+\frac{R}{\tau}c(t)$ (25) $\displaystyle e^{\frac{t}{\tau}}\frac{du(t)}{dt}$ $\displaystyle=-\frac{1}{\tau}e^{\frac{t}{\tau}}u(t)+e^{\frac{t}{\tau}}\frac{R} {\tau}c(t)$ $\displaystyle e^{\frac{t}{\tau}}\frac{du(t)}{dt}+\frac{1}{\tau}e^{\frac{t}{ \tau}}u(t)$ $\displaystyle=e^{\frac{t}{\tau}}\frac{R}{\tau}c(t)$ $\displaystyle\frac{d}{dt}\left(e^{\frac{t}{\tau}}u(t)\right)$ $\displaystyle=e^{\frac{t}{\tau}}\frac{R}{\tau}c(t)$ $\displaystyle\int_{t_{k-1}}^{t_{k}}\frac{d}{dt}\left(e^{\frac{t}{\tau}}u(t)\right)$ $\displaystyle=\int_{t_{k-1}}^{t_{k}}e^{\frac{t}{\tau}}\frac{R}{\tau}c(t)dt$ $\displaystyle e^{\frac{t_{k}}{\tau}}u(t_{k})-e^{\frac{t_{k-1}}{\tau}}u(t_{k-1})$ $\displaystyle=\left(e^{\frac{t_{k}}{\tau}}-e^{\frac{t_{k-1}}{\tau}}\right)Rc_{ t_{k}}$ $\displaystyle u(t_{k})$ $\displaystyle=e^{-\frac{\Delta t_{k}}{\tau}}u(t_{k-1})+\left(1-e^{-\frac{ \Delta t_{k}}{\tau}}\right)Rc_{t_{k}}$ $\displaystyle u(t_{k})$ $\displaystyle\approx e^{-\frac{\Delta t_{k}}{\tau}}u(t_{k-1})+\Delta t_{k} \frac{R}{\tau}c_{t_{k}}.$ $$

Rearranging, assuming $R=\tau$ without input decay, we replace $\Delta t_{k}$ with $\Delta$, as done in S4 Gu et al. (). We obtain the discrete form:

$$ $u_{t}=e^{\Delta\gamma}u_{t-1}+\Delta c_{t}.$ (26) $$

Finally, replacing the original part $\tilde{\gamma}=-\gamma+i\theta$ gives the PRF neuron with sequential computation:

$$ $\displaystyle\tilde{u_{t}}$ $\displaystyle=\exp\left(\Delta\left(-\frac{1}{\tau}+i\theta\right)\right) \tilde{u}_{t-1}+\Delta c_{t},$ (27) $\displaystyle s_{t}$ $\displaystyle=\mathcal{H}(\Re\{\tilde{u_{t}}\}-V_{\mathrm{th}}).$ (28) $$

## Appendix D The Frequency Response for PRF neuron

This section mainly discusses the frequency response of the dynamic reset LIF neuron. Recall the membrane potential dynamic equation:

$$ $\displaystyle\tilde{u_{t}}$ $\displaystyle=\exp\left(\Delta\left(-\frac{1}{\tau}+i\theta\right)\right) \tilde{u}_{t-1}+\Delta c_{t},$ (29) $\displaystyle s_{t}$ $\displaystyle=\mathcal{H}(\Re\{\tilde{u_{t}}\}-V_{\mathrm{th}}).$ (30) $$

This dynamic process can be regarded as a damped harmonic oscillator. The real part of the membrane potential, $\Re\{u_{t}\}$, represents the displacement. When the displacement exceeds a certain value, this model will issue a signal. Here, $\Delta$ represents the timestep size, and $c_{t}$ is the input for each timestep, driven by an external force.

First, define $\tilde{\gamma}\triangleq\left(-\frac{1}{\tau}+i\theta\right)$. The model can then be expressed as:

$$ $\displaystyle\tilde{u}_{t}$ $\displaystyle=\exp\left(\Delta\tilde{\gamma}\right)\tilde{u}_{t-1}+\Delta c_{t},$ (31) $\displaystyle\frac{\tilde{u}_{t}-\tilde{u}_{t-1}}{\Delta}$ $\displaystyle=\frac{\exp\left(\Delta\tilde{\gamma}\right)-1}{\Delta}\tilde{u}_ {t-1}+c_{t},$ (32) $$

To explore the numerical effects in the experiment, we obtain the approximate ODE by using $x$ to represent $u$:

$$ $\frac{dx}{d\Delta}-\left(\frac{\exp\left(\Delta\tilde{\gamma}\right)-1}{\Delta }\right)x=c_{t},$ (33) $$

using a first-order Taylor expansion approximation, we obtain:

$$ $\frac{dx}{d\Delta}-\tilde{\gamma}x=c_{t},$ (34) $$

We assume the driven input is an oscillation, $c_{t}=c_{0}\exp(i\omega\Delta)$, with a constant base amplitude $c_{0}$ and a variable angular frequency $\omega$. The position of the membrane potential $x$ will oscillate in resonance as:

$$ $x=x_{\omega}\exp(i\omega\Delta),$ (35) $$

where $x_{\omega}$ is the amplitude as a function of the external excitation frequency. We then have:

$$ $\displaystyle\dot{x}=i\omega x_{\omega}\exp(i\omega\Delta).$ (36) $$

Substituting Equation 36 into Equation 34 gives:

$$ $\displaystyle i\omega x_{\omega}\exp(i\omega\Delta)-\tilde{\gamma}x_{\omega} \exp(i\omega\Delta)$ $\displaystyle=c_{0}\exp(i\omega\Delta)$ (37) $\displaystyle i\omega x_{\omega}-\tilde{\gamma}x_{\omega}$ $\displaystyle=c_{0}$ (38) $\displaystyle\left(i\omega-\tilde{\gamma}\right)x_{\omega}$ $\displaystyle=c_{0}$ (39) $$

Rearranging the equation yields:

$$ $\displaystyle\frac{x_{\omega}}{c_{0}}$ $\displaystyle=\frac{1}{i\omega-\tilde{\gamma}}$ (40) $\displaystyle=\frac{1}{\frac{1}{\tau}+i(\omega-\theta)}.$ (41) $$

Thus, we have:

$$ $\displaystyle\Re\left\{\frac{x_{\omega}}{c_{0}}\right\}=\frac{1/\tau}{(\frac{1 }{\tau})^{2}+(\omega-\theta)^{2}}$ (42) $\displaystyle\Im\left\{\frac{x_{\omega}}{c_{0}}\right\}=\frac{-\omega+\theta}{ (\frac{1}{\tau})^{2}+(\omega-\theta)^{2}}$ (43) $$

The magnitude can be obtained by taking the modulus, which varies with $\omega$ and $\theta$:

$$ $\displaystyle\left\lvert\frac{x_{\omega}}{c_{0}}\right\rvert$ $\displaystyle=\sqrt{\Re\left\{\frac{x_{\omega}}{c_{0}}\right\}^{2}+\Im\left\{ \frac{x_{\omega}}{c_{0}}\right\}^{2}}$ (44) $\displaystyle=\frac{1}{\sqrt{(\frac{1}{\tau})^{2}+(\omega-\theta)^{2}}}$ (45) $$

Assuming $\omega>0$, the value of $\omega$ corresponding to the maximum of the magnitude can be found:

$$ $\displaystyle d\left\lvert\frac{x_{\omega}}{c_{0}}\right\rvert/d\omega$ $\displaystyle=0$ (46) $\displaystyle-\frac{1}{2}\left((\frac{1}{\tau})^{2}+(\omega-\theta)^{2}\right) ^{-\frac{3}{2}}\times 2(\omega-\theta)$ $\displaystyle=0$ (47) $\displaystyle\omega$ $\displaystyle=\theta$ (48) $$

Therefore, the resonant frequency $\omega$ at the point of maximum magnitude coincides with $\theta$, and $\max\left(\left\lvert\frac{x_{\omega}}{c_{0}}\right\rvert\right)=\tau$.

## Appendix E Deployment analysis of PRF Neuron

This model can also be easily deployed on neuromorphic chips. After training, the coefficients can merge together. Recalling the PRF Neuron as described in Equation 19, we explicitly expand the real and imaginary parts as follows:

Figure: Figure 5: The comparison of LIF Neuron (left) and PRF Neuron (right) for inference deployment.
Refer to caption: x3.png

$$ $\displaystyle\tilde{u}_{t}$ $\displaystyle=\exp\left(\Delta\left(-\frac{1}{\tau}+i\theta\right)\right) \tilde{u}_{t-1}+\Delta c_{t}$ (49) $\displaystyle=\left(\exp{\left(-\frac{\Delta}{\tau}\right)}\cos{(\Delta\theta) }+i\exp{\left(-\frac{\Delta}{\tau}\right)}\sin{(\Delta\theta)}\right)\tilde{u} _{t-1}+\Delta c_{t}$ $$

We use the symbols $\Re{\tilde{u}_{t}}\in\mathbb{R}$ and $\Im{\tilde{u}_{t}}\in\mathbb{R}$ to denote the real and imaginary parts of the membrane potential $\tilde{u}_{t}\in\mathbb{C}$, respectively:

$$ $\begin{pmatrix}\Re\{\tilde{u}_{t}\}\\ \Im\{\tilde{u}_{t}\}\end{pmatrix}=\begin{pmatrix}\varphi_{\mathrm{re}}&- \varphi_{\mathrm{im}}\\ \varphi_{\mathrm{im}}&\varphi_{\mathrm{re}}\end{pmatrix}\begin{pmatrix}\Re\{ \tilde{u}_{t-1}\}\\ \Im\{\tilde{u}_{t-1}\}\end{pmatrix}+\begin{pmatrix}\Delta c_{t}\\ 0\end{pmatrix},$ (50) $$

where the coefficients $\varphi_{\mathrm{re}},\varphi_{\mathrm{im}}\in\mathbb{R}$ are the merged parameters:

$$ $\displaystyle\varphi_{\mathrm{re}}$ $\displaystyle=\exp{\left(-\frac{\Delta}{\tau}\right)}\cos{(\Delta\theta)},$ (51) $\displaystyle\varphi_{\mathrm{im}}$ $\displaystyle=\exp{\left(-\frac{\Delta}{\tau}\right)}\sin{(\Delta\theta)}.$ $$

Finally, the spike is output based on the real part $\Re{\tilde{u}_{t}}$. For simplicity, we use $u_{t}$ and $r_{t}$ to denote $\Re{\tilde{u}_{t}}$ and $\Im{\tilde{u}_{t}}$, respectively. Thus, the explicit iteration of the PRF can be written as:

$$ $\displaystyle u_{t}$ $\displaystyle=\varphi_{\mathrm{re}}u_{t-1}-\varphi_{\mathrm{im}}r_{t-1}+\Delta c _{t},$ (52) $\displaystyle r_{t}$ $\displaystyle=\varphi_{\mathrm{im}}u_{t-1}+\varphi_{\mathrm{re}}r_{t-1},$ (53) $\displaystyle s_{t}$ $\displaystyle=\mathcal{H}\left(u_{t}-V_{\mathrm{th}}\right).$ (54) $$

Intuitively, compared to the LIF model, the PRF introduces two additional multiplication operations and one extra addition operation for inference, along with an extra hidden state that needs to be saved, as shown in Figure 5. Furthermore, this neuron model, with its double hidden state, can also be easily deployed on neuromorphic chips, similar to how the AHP model Rao et al. () was deployed on the Loihi chip Davies et al. ().

## Appendix F Proof

### F.1 Proof of Theorem. 1

###### Proof.

Firstly, consider the Adaptive LIF (ALIF) model Bellec et al. (), where the threshold adapts according to recent firing activity. The dynamic threshold is given by Equation 55 to Equation 57:

$$ $\displaystyle A_{t}$ $\displaystyle=V_{\mathrm{th}}+\beta a_{t}$ (55) $\displaystyle z_{t}$ $\displaystyle=\mathcal{H}(u_{t}-A_{t})$ (56) $\displaystyle a_{t+1}$ $\displaystyle=\rho a_{t}+z_{t}$ (57) $$

Here, the decay factor $\rho$ is given by $e^{-\delta t/\tau_{a}}$, where $\tau_{a}$ is the adaptation time constant, as described in ALIF Bellec et al. ().

Intuitively, when $\tau_{a}\gg\delta t$, we have $\rho\rightarrow 1$, simplifying Equation 57 to $a_{t+1}=a_{t}+z_{t}$. The variable $z_{t}$ can be further deduced as follows:

$$ $\displaystyle z_{t}$ $\displaystyle=\mathcal{H}(u_{t}-A_{t})=\begin{cases}0,&\text{if }u_{t}<A_{t}\\ 1,&\text{if }u_{t}\geq A_{t}\end{cases}$ (58) $$

Substituting Equation 58 into Equation 57, we obtain:

$$ $\displaystyle a_{t+1}$ $\displaystyle=\begin{cases}a_{t},&\text{if }u_{t}<A_{t}\\ a_{t}+1,&\text{if }u_{t}\geq A_{t}\end{cases}$ (59) $$

Therefore, Equation 55 can be further expanded using Equation 59:

$$ $\displaystyle A_{t}$ $\displaystyle=V_{\mathrm{th}}+\beta a_{t}$ (60) $\displaystyle a_{t}$ $\displaystyle=\begin{cases}a_{t-1},&\text{if }u_{t}<A_{t}\\ a_{t-1}+1,&\text{if }u_{t}\geq A_{t}\end{cases}$ (61) $$

Finally, we can see that Equation 60 - Equation 61 are equivalent to Equation 10 - Equation 11.
∎

###### Proof.

### F.2 Proof of Theorem. 2

###### Proof.

Firstly, we recall the dynamic iteration of the PRF model:

$$ $\displaystyle\tilde{u}_{t}$ $\displaystyle=\exp\left(\Delta\left(-\frac{1}{\tau}+i\theta\right)\right) \tilde{u}_{t-1}+\Delta c_{t},$ (62) $\displaystyle s_{t}$ $\displaystyle=\mathcal{H}\left(\Re\{\tilde{u}_{t}\}-V_{\mathrm{th}}\right)$ (63) $$

Next, let $\Delta=1$ and $\theta=0$, allowing the model to be rewritten as:

$$ $\displaystyle u_{t}$ $\displaystyle=\exp\left(-\frac{1}{\tau}\right)u_{t-1}+c_{t},$ (64) $\displaystyle s_{t}$ $\displaystyle=\mathcal{H}\left(u_{t}-V_{\mathrm{th}}\right)$ (65) $$

At this point, $u_{t}\in\mathbb{R}$, meaning it only contains a real part, with the decay affecting only the real component. Setting $\theta=0$ can be interpreted as removing the reset process. Furthermore, the exponential term $\exp\left(-\frac{1}{\tau}\right)$ can be approximated using the first-order Taylor expansion:

$$ $\displaystyle\exp\left(-\frac{1}{\tau}\right)$ $\displaystyle\approx 1-\frac{1}{\tau}$ (66) $$

Finally, the dynamic equation can be rewritten as:

$$ $\displaystyle u_{t}$ $\displaystyle=\left(1-\frac{1}{\tau}\right)u_{t-1}+c_{t},$ (67) $\displaystyle s_{t}$ $\displaystyle=\mathcal{H}\left(u_{t}-V_{\mathrm{th}}\right)$ (68) $$

Equation 67 and Equation 68 are equivalent to the standard LIF model without the reset process, as given in Equation 2.
∎

###### Proof.

### F.3 Proof of Theorem. 3

###### Proof.

Firstly, we recall the dynamic iteration of the PRF model:

$$ $\displaystyle\tilde{u}_{t}$ $\displaystyle=\exp\left(\Delta\left(-\frac{1}{\tau}+i\theta\right)\right) \tilde{u}_{t-1}+\Delta c_{t},$ (69) $$

assuming $c_{t}\sim\mathcal{N}(0,\sigma^{2})$ is normally distributed with zero mean and variance $\sigma^{2}$ at time-invariant, and $\Delta$, $\tau$, and $\theta$ are constants. Define the complex constant $A$:

$$ $A=\exp\left(\Delta\left(-\frac{1}{\tau}+i\theta\right)\right)=\exp\left(-\frac {\Delta}{\tau}\right)\exp\left(i\Delta\theta\right).$ (70) $$

Note that the magnitude of $A$ is:

$$ $|A|=e^{-\Delta/\tau}<1,$ (71) $$

since $\tau>\Delta>0$. In further, where the membrane potential $u_{t}=\Re\{\tilde{u}_{t}\}$ is the real part of $\tilde{u}_{t}$.

Next, we expand $\tilde{u}_{t}$ recursively:

$$ $\displaystyle\tilde{u}_{t}$ $\displaystyle=A\tilde{u}_{t-1}+\Delta c_{t}$ (72) $\displaystyle=A(A\tilde{u}_{t-2}+\Delta c_{t-1})+\Delta c_{t}$ (73) $\displaystyle=A^{2}\tilde{u}_{t-2}+A\Delta c_{t-1}+\Delta c_{t}$ (74) $\displaystyle\vdots$ (75) $\displaystyle=A^{t}\tilde{u}_{0}+\Delta\sum_{k=1}^{t}A^{t-k}c_{k}.$ (76) $$

Now, we calculate the expectation and variation of $u_{t}$. since $c_{t}$ are independent and identically distributed with $\mathbb{E}[c_{t}]=0$ and $\operatorname{Var}(c_{t})=\sigma^{2}$, we can compute the expected value of $\tilde{u}_{t}$:

$$ $\displaystyle\mathbb{E}[\tilde{u}_{t}]$ $\displaystyle=\mathbb{E}\left[A^{t}\tilde{u}_{0}+\Delta\sum_{k=1}^{t}A^{t-k}c_ {k}\right]$ (77) $\displaystyle=A^{t}\tilde{u}_{0}+\Delta\sum_{k=1}^{t}A^{t-k}\mathbb{E}[c_{k}]$ (78) $\displaystyle=A^{t}\tilde{u}_{0}.$ (79) $$

As $|A|<1$, it follows that:

$$ $\lim_{t\to\infty}\mathbb{E}[u_{t}]=\Re\{\lim_{t\to\infty}\mathbb{E}[\tilde{u}_ {t}]\}\to 0$ (80) $$

Next, compute the variance of $\tilde{u}_{t}$ :

$$ $\displaystyle\operatorname{Var}(\tilde{u}_{t})$ $\displaystyle=\mathbb{E}\left[|\tilde{u}_{t}|^{2}\right]-|\mathbb{E}[\tilde{u} _{t}]|^{2}$ (81) $\displaystyle=\mathbb{E}\left[\left|A^{t}\tilde{u}_{0}+\Delta\sum_{k=1}^{t}A^{ t-k}c_{k}\right|^{2}\right]-|A^{t}\tilde{u}_{0}|^{2}$ (82) $\displaystyle=2A^{t}\tilde{u}_{0}\Delta\sum_{k=1}^{t}A^{t-k}\mathbb{E}\left[c_ {k}\right]+\Delta^{2}\mathbb{E}\left[\left|\sum_{k=1}^{t}A^{t-k}\mathbb{E}[c_{ k}]\right|^{2}\right]$ (83) $\displaystyle=\Delta^{2}\mathbb{E}\left[\left|\sum_{k=1}^{t}A^{t-k}c_{k}\right |^{2}\right].$ (84) $$

Since $c_{k}$ are independent and have zero mean, we have:

$$ $\displaystyle\mathbb{E}\left[\left|\sum_{k=1}^{t}A^{t-k}c_{k}\right|^{2}\right]$ $\displaystyle=\sum_{k=1}^{t}\sum_{l=1}^{t}A^{t-k}\overline{A^{t-l}}\mathbb{E} \left[c_{k}c_{l}\right]$ (85) $\displaystyle=\sum_{k=1}^{t}|A|^{2(t-k)}\mathbb{E}\left[|c_{k}|^{2}\right], \quad\text{since }\mathbb{E}\left[c_{k}c_{l}\right]=0\text{ for }k\neq l$ (86) $\displaystyle=\sigma^{2}\sum_{k=1}^{t}|A|^{2(t-k)}.$ (87) $$

Therefore, the variance of $u_{t}$ is:

$$ $\operatorname{Var}(u_{t})=\Delta^{2}\sigma^{2}\sum_{k=1}^{t}|A|^{2(t-k)}= \Delta^{2}\sigma^{2}\sum_{m=0}^{t-1}|A|^{2m}.$ (88) $$

Since $|A|=e^{-\Delta/\tau}$ as deduced in Equation 71, we have:

$$ $|A|^{2m}=e^{-2\Delta m/\tau}.$ (89) $$

Thus,

$$ $\operatorname{Var}(u_{t})=\Delta^{2}\sigma^{2}\sum_{m=0}^{t-1}e^{-2\Delta m/ \tau}.$ (90) $$

This is a finite geometric series with first term equals to $1$ and common ratio $r=e^{-2\Delta/\tau}$:

$$ $\sum_{m=0}^{t-1}e^{-2\Delta m/\tau}=\frac{1-r^{t}}{1-r}.$ (91) $$

Therefore,

$$ $\operatorname{Var}(u_{t})=\Delta^{2}\sigma^{2}\frac{1-e^{-2\Delta t/\tau}}{1-e ^{-2\Delta/\tau}}.$ (92) $$

As $t\to\infty$, $e^{-2\Delta t/\tau}\to 0$, so the variance approaches:

$$ $\lim_{t\to\infty}\operatorname{Var}(u_{t})=\frac{1}{1-e^{-2\Delta/\tau}} \approx\frac{\Delta^{2}\sigma^{2}}{2\Delta/\tau}=\frac{\tau\Delta}{2}\sigma^{2},$ (93) $$

this is a finite constant, in summary we could get the distribution after $t\to\infty$:

$$ $u_{t}\sim\mathcal{N}\left(0,\frac{\tau\Delta}{2}\sigma^{2}\right),$ (94) $$

which implies that $u_{t}$ is bounded in probability. This indicates that $u_{t}$ does not diverge but instead stabilizes to a specific distribution related to the input distribution and hyperparameters. It ensures that despite the randomness introduced by the inputs $c_{t}$, the neuron’s response remains predictable in distribution.
∎

###### Proof.

## Appendix G Parallel Perspective on Solving Long-Range Learning Problem

From the subsection of Problem Formulation 2, we gain the insight that the gradient is proportional to the inner product of the kernel and previous layer spikes.

$$ $\displaystyle\nabla_{\bm{W}^{l}}\mathcal{L}\propto\frac{\partial\mathcal{L}}{ \partial u^{l}_{T}}\sum_{t=1}^{T}\frac{\partial u^{l}_{T}}{\partial u^{l}_{t}} \frac{\partial u^{l}_{t}}{\partial\bm{W}^{l}}\propto\underbrace{\sum_{t=1}^{T} \beta^{(t)}s^{l-1}_{t}=\left<\mathbf{K}_{T},\mathbf{S}^{l-1}_{T:1}\right>}_{ \text{Parallel Perspective}}$ (95) $$

To verify this insight, we define three extreme situations (Fast Decay, Slow Decay and Slow Decay with Resonate) to explore this from the parallel kernel perspective, as shown in Figure 6.

Figure: Figure 6: The parallel kernel perspective for learning long-range abilities.
Refer to caption: x5.png

##### Fast Decay

Fast decay may cause the problem of long-range dependency vanishing, as shown in Figure 6(a). The kernel $\mathbf{K}_{T}$ decays rapidly as the time step $T$ increases. In this case, early spikes (represented by $a$) dominate the gradient calculation, while contributions from later spikes (represented by $b$) are almost negligible. As a result, the gradient $\partial\mathcal{L}/\partial\bm{W}^{l}$ is mostly proportional to $a$, leading to the vanishing of long-range dependencies. The network struggles to learn and retain information from distant time steps, causing the gradients to vanish and impairing the learning of long-term dependencies.

##### Slow Decay

Slow Decay may relive the vanishing problem, but may cause gradients ambiguity as shown in Figure.6(b). This situation shows the kernel decays more slowly over time, which balances the contributions from both early and late spikes (denoted as $a_{1},b_{1}$ and $a_{2},b_{2}$). However, this slow decay introduces a new problem—ambiguity. When the contributions from different parts of the sequence are similar (e.g., $a_{1}+b_{1}\approx a_{2}+b_{2}$), the network finds it difficult to distinguish between these cases. This ambiguity can confuse the learning process, as the network may not correctly interpret or differentiate between distinct temporal patterns.

##### Slow Decay with Resonate

Resonate could may relive the ambiguity problem to capture the resonate information under the long range, as shown in Figure.6(c). the kernel $\mathbf{K}_{T}$ oscillates or resonates, effectively capturing contributions from both early and late spikes. This resonance allows the network to amplify and preserve significant information across the entire sequence, represented by $a_{1}$ and $b_{1}$. Such behavior is advantageous for learning long-range dependencies, as it helps maintain the gradient information over time. The network can now better distinguish between different temporal patterns, leading to improved learning and retention of long-term information, which is crucial for tasks involving long sequences.

## Appendix H The Architecture for Long Range Arena Task

This Section introduces the Spike-Driven Temporal and Channel Mixing (SD-TCM) Module, as shown in Figure 7. This design philosophy mainly stems from the token and channel mixing in the transformer and S4 for solving more difficult sequence tasks. Firstly, we introduce the spatial neuron by considering the other side of LIF neurons. Secondly, we present the mixer module, which combines the PRF neuron and spatial neuron. Furthermore, we compare the computation complexity and theoretical energy consumption (Detail in Appendix.I).

Figure: Figure 7: Diagram of the Spike-Driven Temporal and Channel Mixer (SD-TCM) Module.
Refer to caption: x6.png

The design philosophy of this module stems from combining token mixing with channel mixing, a common practice in transformer and S4 modules. The transformer uses a self-attention block for token mixing and a 2-layer MLP for channel mixing Vaswani et al. (). The S4 employs SSMs for token mixing and GLU for channel mixing Gu et al. ().
Firstly, we propose the Spatial neuron utilized for channel mixing. We consider the limitation of the time constant $\tau$ close to the time scale, focusing on the instantaneous information for each timestep:

$$ $\displaystyle\lim_{\tau\to 1^{+}}u_{t}$ $\displaystyle=\left(1-\frac{1}{\tau}\right)u_{t-1}+c_{t}\approx c_{t},$ (96) $$

then the output of Spatial Neuron replaced as $s_{t}=\mathcal{H}(c_{t}-V_{\text{th}})$. This type of neuron can also be observed in motor neurons in biology with extreme huge decay. According to Theorem 2, we gain the insight that the LIF neuron is a subset of the PRF Neuron. The Spatial Neuron is also a subset of the LIF Neuron, aiming to focus on instantaneous spatial information without temporal information. Consequently, we derived a special case of the LIF neuron, which we term the Spatial Neuron. In further, we use the trainable amplitude for Spatial LIF Neuron with the output $\{0,\alpha\}$. Where the amplitude constant $\alpha\in\mathbb{R}^{+}$ is trainable parameters with always initialize as $1$. After training, the amplitude constant $\alpha$ could merge to the following Linear layer during the inference. That means $(\alpha\bm{s}_{t})\times\bm{W}\equiv\bm{s}_{t}\times(\alpha\bm{W})$.

Secondly, we design the SD-TCM module consists of three main components: the Spatial Neuron $\mathcal{SN(\cdot)}$, PRF Neuron $\mathcal{TN(\cdot)}$, fully-connected layer $\mathrm{Linear}(\cdot)$. To keep the spike-driven feature, we use the membrane shortcut residual connect Hu et al. () like spike-driven transformer Yao et al. (). Given a input sequence $I\in\mathbb{R}^{T\times N\times D_{\mathrm{in}}}$, the embedding the sequence of $N$ flattened spike patches with $D$ dimensional channel,

$$ $\displaystyle\mathbf{U}_{T}^{1}$ $\displaystyle=\text{Embedding}(I),$ $\displaystyle I$ $\displaystyle\in\mathbb{R}^{T\times N\times D_{\mathrm{in}}},\quad\mathbf{U}_{ T}\in\mathbb{R}^{T\times N\times D},$ (97) $$

where $T$ denote timestep while aligning with the sequence length. Then the SD-TCM is written as:

$$ $\displaystyle\mathbf{S}_{T}^{l}$ $\displaystyle=\mathcal{TN}(\mathbf{U}_{T}^{l}),$ $\displaystyle\mathbf{S}_{T}^{l}$ $\displaystyle\in\mathbb{B}^{T\times N\times D},\quad l=1,2,\dots,L$ (98) RPE $\displaystyle=\mathbf{U}_{T}^{l}+\text{Linear}(\mathbf{S}_{T}^{l}),$ RPE $\displaystyle\in\mathbb{R}^{T\times N\times D},\quad l=1,2,\dots,L$ (99) $\displaystyle{\mathbf{S}^{\prime}}_{T}^{l}$ $\displaystyle=\mathcal{SN}(\text{RPE}),$ $\displaystyle{\mathbf{S}^{\prime}}_{T}^{l}$ $\displaystyle\in\mathbb{B}^{T\times N\times D},\quad l=1,2,\dots,L$ (100) $\displaystyle\mathbf{U}_{T}^{l+1}$ $\displaystyle=\text{RPE}+\text{Linear}({\mathbf{S}^{\prime}}_{T}^{l}),$ $\displaystyle\mathbf{U}_{T}^{l+1}$ $\displaystyle\in\mathbb{R}^{T\times N\times D},\quad l=1,2,\dots,L$ (101) $$

Where $\mathbb{B}\triangleq\{0,1\}$ is the binary value set. After the $L^{\mathrm{th}}$ layer, the following output as the input of classifier or other head for corresponding task. This block keep the spike-driven with two properties: event-driven and binary spike-based communication. The former means that no computation is triggered when the input is zero. The binary restriction in the latter indicates that there are only additions.

The original S4 layer is unidirectional or causal, which is an unnecessary constraint for the classification tasks appearing in LRA. Goel et al. () propose a bidirectional version of S4 that simply concatenates two S4 convolution kernels back-to-back Gu et al. (). As same the bidirectional model implemented in S4 block. We implement the bidirectional model by replacing the Equation 99 as the following equation:

$$ $\displaystyle\mathbf{S}_{T}^{l}$ $\displaystyle=\text{{Concat}}\left(\mathcal{TN}(\mathbf{U}_{T}^{l}),\text{{Rev }}\left(\mathcal{TN}(\text{{Rev}}(\mathbf{U}_{T}^{l}))\right)\right),$ $\displaystyle\mathbf{S}_{T}^{l}$ $\displaystyle\in\mathbb{B}^{T\times N\times 2D},$ (102) $$

where Concat and Rev means the concatenate and reverse operation along the channel and timestep dimension in respectively. We simply pass the input sequence through an PRF Neuron, and also reverse it and pass it through an independent second PRF Neuron. These spiking outputs are concatenated and passed through a position wise linear layer. Keep the same with S4 and Shashimi Goel et al. (), the following Linear layer ($\bm{W}\in\mathbb{R}^{2D\times D}$) will change the input feature dimension as the double in the bidirectional model.

## Appendix I The Theoretical Analysis of Power consumption

The two tables provide a detailed comparison of inference complexity and energy consumption across various models. Table 6 compares the computational complexity involved in the inference phase for different models. Table 7 evaluates the energy consumption of these models during the token mixing and channel mixing stages.

**Table 6: The comparison for inference complexity. The abstract formulations $x_{t}$, $u_{t}$, and $y_{t}$ represent the input, hidden state, and output, respectively. The symbols $\mathbb{R}$ and $\mathbb{C}$ denote real and complex number sets. Where the symbols $\mathcal{H}$ and $\Re$ denote the Heaviside function and the real part of the complex number.**
|  | Dynamic Equation | Variables | Parameters | Infer. Complexity |
| --- | --- | --- | --- | --- |
| SSMs | $\bm{u}_{t}=\bm{A}\bm{u}_{t-1}+\bm{B}\bm{x}_{t}$<br>$\bm{y}_{t}=\bm{C}\bm{u}_{t}$ | $\bm{x}_{t}\in\mathbb{R}^{D}$<br>$\bm{u}_{t}\in\mathbb{R}^{H}$<br>$\bm{y}_{t}\in\mathbb{R}^{D}$ | $\bm{A}\in\mathbb{R}^{H\times H}$<br>$\bm{B}\in\mathbb{R}^{D\times H}$<br>$\bm{C}\in\mathbb{R}^{H\times D}$ | $O(H^{2}+2DH)$ |
| Spikinglized<br>SSMs | $\bm{u}_{t}=\bm{A}\bm{u}_{t-1}+\bm{B}\bm{x}_{t}$<br>$\bm{y}_{t}=\mathcal{H}(\bm{C}\bm{u}_{t}-V_{\mathrm{th}})$ | $\bm{x}_{t}\in\mathbb{R}^{D}$<br>$\bm{u}_{t}\in\mathbb{R}^{H}$<br>$\bm{y}_{t}\in\{0,1\}^{D}$ | $\bm{A}\in\mathbb{R}^{H\times H}$<br>$\bm{B}\in\mathbb{R}^{D\times H}$<br>$\bm{C}\in\mathbb{R}^{H\times D}$ | $O(H^{2}+2DH)$ |
| Ours | $\bm{u}_{t}=\bm{A}\odot\bm{u}_{t-1}+\bm{B}\odot\bm{x}_{t}$<br>$\bm{y}_{t}=\mathcal{H}(\Re\{\bm{u}_{t}\}-V_{\mathrm{th}})$ | $\bm{x}_{t}\in\mathbb{R}^{D}$<br>$\bm{u}_{t}\in\mathbb{C}^{D}$<br>$\bm{y}_{t}\in\{0,1\}^{D}$ | $\bm{A}\in\mathbb{C}^{D}$<br>$\bm{B}\in\mathbb{R}^{D}$ | $O(5D)$ |

**Table 7: Energy evaluation. $R$ denote the spike firing rates (the proportion of non-zero elements in the neuron output). ($\sigma$: Sigmoid activation function, $\mathcal{H}$: Heaviside function and $\Re$: real part of the complex number, and Ter means the ternary output $\{-1,0,1\}$ with a dynamic threshold).**
|  | Token Mixing | Channel Mixing |  |  |
| --- | --- | --- | --- | --- |
|  | Comput. | Energy | Comput. | Energy |
| S4-LegS | $\bm{u}_{t}=\bm{A}\bm{u}_{t-1}+\bm{B}\bm{x}_{t}$<br>$\bm{y}_{t}=\bm{C}\bm{u}_{t}+\bm{D}\bm{x}_{t}$ | $E_{\mathrm{MAC}}\cdot(H^{2}+DH)$<br>$E_{\mathrm{MAC}}\cdot(HD+D^{2})$ | $\bm{n}_{t}=\textbf{Linear}_{1}(\bm{y}_{t})$<br>$\bm{m}_{t}=\textbf{Linear}_{2}(\bm{y}_{t})$<br>$\bm{o}_{t}=\bm{n}_{t}\odot\sigma(\bm{m}_{t})$ | $E_{\mathrm{MAC}}\cdot D^{2}$<br>$E_{\mathrm{MAC}}\cdot D^{2}$<br>$E_{\mathrm{MAC}}\cdot 2D+E_{\mathrm{M}}\cdot D$ |
| Binary S4D | $\bm{u}_{t}=\bm{A}\bm{u}_{t-1}+\bm{B}\bm{x}_{t}$<br>$\bm{y}_{t}=\mathcal{H}(\bm{C}\bm{u}_{t}+\bm{D}\bm{x}_{t})$ | $E_{\mathrm{MAC}}\cdot(H^{2}+DH)$<br>$E_{\mathrm{MAC}}\cdot(HD+D^{2})$ | $\bm{n}_{t}=\textbf{Linear}_{1}(\bm{y}_{t})$<br>$\bm{m}_{t}=\textbf{Linear}_{2}(\bm{y}_{t})$<br>$\bm{o}_{t}=\bm{n}_{t}\odot\sigma(\bm{m}_{t})$ | $E_{\mathrm{AC}}\cdot R\cdot D^{2}$<br>$E_{\mathrm{AC}}\cdot R\cdot D^{2}$<br>$E_{\mathrm{MAC}}\cdot 2D+E_{\mathrm{M}}\cdot D$ |
| GSU | $\bm{u}_{t}=\bm{A}\bm{u}_{t-1}+\bm{B}\bm{x}_{t}$<br>$\bm{y}_{t}=\bm{C}\bm{u}_{t}+\bm{D}\bm{x}_{t}$ | $E_{\mathrm{MAC}}\cdot(H^{2}+DH)$<br>$E_{\mathrm{MAC}}\cdot(HD+D^{2})$ | $\bm{n}_{t}=\mathrm{Ter}(\bm{W}_{1})\bm{y}_{t}+\bm{b}_{1}$<br>$\bm{m}_{t}=\bm{W}_{2}\mathrm{Ter}(\bm{y}_{t})+\bm{b}_{2}$<br>$\bm{o}_{t}=\textbf{GeLU}(\bm{n}_{t}\odot\bm{m}_{t})$ | $E_{\mathrm{AC}}\cdot R\cdot D^{2}$<br>$E_{\mathrm{AC}}\cdot R\cdot D^{2}$<br>$E_{\mathrm{MAC}}\cdot 2D+E_{\mathrm{M}}\cdot D$ |
| Ours | $\bm{u}_{t}=\bm{A}\odot\bm{u}_{t-1}+\bm{B}\odot\bm{x}_{t}$<br>$\bm{y}_{t}=\mathcal{H}(\Re\{\bm{u}_{t}\}-V_{\mathrm{th}})$<br>$\bm{n}_{t}=\textbf{Linear}_{1}(\bm{y}_{t})+\bm{u}_{t}$ | $E_{\mathrm{M}}\cdot 5D+E_{\mathrm{AC}}\cdot 3D$<br>-<br>$E_{\mathrm{AC}}\cdot(R\cdot D^{2}+D)$ | $\bm{s}_{t}=\mathcal{H}(\bm{n}_{t}-V_{\mathrm{th}})$<br>$\bm{o}_{t}=\textbf{Linear}_{2}(\bm{s}_{t})+\bm{n}_{t}$ | -<br>$E_{\mathrm{AC}}\cdot(R\cdot D^{2}+D)$ |

## Appendix J Description of Datasets and Hyperparameters

All our experiments were conducted on an NVIDIA GeForce RTX 4090 GPU with 24 GB of memory. The specific experimental setup and hyperparameters are detailed in subsection J.1, and the description of the experimental dataset is provided in subsection J.2.

### J.1 Task Specific Hyperparameters

Here we specify any task-specific details, hyperparameter or architectural differences from the defaults outlined above.

#### J.1.1 Sequential MNIST & Permuted Sequential MNIST

For Figure.4, we used a neural network architecture with layers of size 1-64-256-256-10 ($87k$ training parameters) with 256 batch size for training 200 epochs. All experiments were conducted using the same random seeds.

#### J.1.2 Long Range Arena

The total hyperparameters configure is shown in Table 8.

**Table 8: Hyperparameters for LRA Task**
| Task | Depth | Channels | Norm | Pre-norm | Dropout | LR | Neuron LR | B | Epochs | WD | ($\Delta_{min}$, $\Delta_{max}$) |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| ListOps | 8 | 128 | BN | False | 0 | 0.005 | 0.001 | 50 | 40 | 0.05 | (0.001, 0.1) |
| Text | 6 | 256 | BN | True | 0 | 0.005 | 0.001 | 16 | 32 | 0.05 | (0.001, 0.1) |
| Retrieval | 6 | 256 | BN | True | 0 | 0.005 | 0.001 | 32 | 20 | 0.05 | (0.001, 0.1) |
| Image | 6 | 512 | BN | False | 0.1 | 0.005 | 0.001 | 50 | 200 | 0.05 | (0.001, 0.1) |
| Pathfinder | 6 | 256 | BN | True | 0.05 | 0.005 | 0.001 | 64 | 200 | 0.03 | (0.001, 0.1) |

### J.2 Dataset Details

We provide more context and details for (p)s-MNIST and each tasks of the LRA Tay et al. (). Note that we follow the same data pre-processing steps as Gu et al. (), which we also include here for completeness. The following describe mainly refer from Smith et al. ().

- •
Sequential MNIST: (sMNIST) 10-way digit classification from a $28\times 28$ grayscale image of a handwritten digit, where the input image is flattened into a $784$-length scalar sequence.
- •
Permuted Sequential MNIST: (psMNIST) 10-way digit classification from a $28\times 28$ grayscale image of a handwritten digit, where the input image is flattened into a $784$-length scalar sequence. This sequence is then permuted using a fixed order.
- •
ListOps: A lengthened version of the dataset presented by Nangia & Bowman (). Given a nested set of mathematical operations (such as min and max) and integer operands in the range zero to nine, expressed in prefix notation with brackets, compute the integer result of the mathematical expression (e.g. [max 2 6 [ min 9 7 ] 0] $\rightarrow$ 7). Characters are encoded as one-hot vectors, with $17$ unique values possible (opening brackets and operators are grouped into a single token). The sequences are of unequal length, and hence the end of shorter sequences is padded with a fixed indicator value, padded to a maximum length of $2,000$. A reserved end-of-sequence token is appended. There are $10$ different classes, representing the integer result of the expression. There are $96,000$ training sequences, $2,000$ validation sequences, and $2,000$ test sequences. No normalization is applied.
- •
Text: Based off of the iMDB sentiment dataset presented by Maas et al. (). Given a movie review, where characters are encoded as a sequence of integer tokens, classify whether the movie review is positive or negative. Characters are encoded as one-hot vectors, with $129$ unique values possible. Sequences are of unequal length, and are padded to a maximum length of $4,096$. There are two different classes, representing positive and negative sentiment. There are $25,000$ training examples and $25,000$ test examples. No validation set is provided. No normalization is applied.
- •
Retrieval: Based off of the ACL Anthology network corpus presented by Radev et al. (). Given two textual citations, where characters are encoded as a sequence of integer tokens, classify whether the two citations are equivalent. The citations must be compressed separately, before being passed into a final classifier layer. This is to evaluate how effectively the network can represent the text. The decoder head then uses the encoded representation to complete the task. Characters are encoded into a one-hot vector with $97$ unique values. Two paired sequences may be of unequal length, with a maximum sequence length of $4,000$. There are two different classes, representing whether the citations are equivalent or not. There are $147,086$ training pairs, $18,090$ validation pairs, and $17,437$ test pairs. No normalization is applied.
- •
Image: Uses the CIFAR-10 dataset presented by Krizhevsky (). Given a $32\times 32$ grayscale CIFAR-10 image as a one-dimensional raster scan, classify the image into one of ten classes. Sequences are of equal length ($1,024$). There are ten different classes. There are $45,000$ training examples, $5,000$ validation examples, and $10,000$ test examples. RGB pixel values are converted to a grayscale intensities, which are then normalized to have zero mean and unit variance (across the entire dataset).
- •
Pathfinder: Based off of the Pathfinder challenge introduced by Linsley et al. (). A $32\times 32$ grayscale image image shows a start and an end point as a small circle. There are a number of dashed lines on the image. The task is to classify whether there is a dashed line (or path) joining the start and end point. There are two different classes, indicating whether there is a valid path or not. Sequences are all of the same length ($1,024$). There are $160,000$ training examples, $20,000$ validation examples, and $20,000$ test examples. The data is normalized to be in the range $[-1,1]$.

## Appendix K The Equivalence of Sequential and Parallel

First, the parallel computation is equivalent to sequential computation during both inference and training phases. This equivalence is clearly illustrated in Figures 8 and 9. For inference equivalence, Figure 8 shows that when applying random input to an LIF neuron using both sequential and parallel computation, the spiking output remains consistent across both methods.

Figure: Figure 8: Verification of inference equivalence for parallel and sequential computation. With random input, the spiking output from parallel computation is equivalent to that of sequential computation.
Refer to caption: x7.png

Figure: Figure 9: Verification of training equivalence between sequential and parallel training. During inference, all metrics are computed using sequential computation. The red curve represents results from parallel training, while the blue curve corresponds to sequential training for both training and inference. This comparison confirms the equivalence of parallel and sequential training.
Refer to caption: x8.png

For training equivalence, Figure 9 shows the training and testing curves for a CIFAR classification task with 1024-length inputs using a 5-layer MLP, where each linear layer is followed by a neuron. Inference is performed with sequential computation, while training is conducted using either sequential (blue curve) or parallel (red curve) methods. After 64 training epochs, the training and testing curves align closely, indicating that the gradient computations from both sequential and parallel training are nearly identical. The parallel method achieves an $8.58\times$ speedup with a sequence length of $1,024$.

## Appendix L The Parallel Implementation Case Statics Details

The figure presents a comparative analysis of sequential and parallel training processes using the torch.profiler tool for a sample case. The trace recording and periods timeline as shown in Figure.10. The corresponding details data statistics in the Table.10 and 10. We designed a simple case using a fully connected layer (FC: $1\times 10$) and an extremely simple architecture to identify the bottlenecks in sequential and parallel training. We use sequential-MNIST (784 length) with 64 batch size as the input.

Figure: Figure 10: The torch.profiler Paszke et al. () tool was used to visualize runtime during forward and backward propagation. The input length for the sample is 784 sequences with a batch size of 64.
Refer to caption: x9.png

**Table 9: Sequential Training Statistics**
| Functions | Duration ($\mu s$) | Num. of Calls |
| --- | --- | --- |
| Forward |  |  |
| Charge | 146,725 | 784 |
| Reset | 150,120 | 784 |
| Fire | 91,398 | 784 |
| Backprop (Autograd) |  |  |
| GraphBackward | 250,161 | 1564 |
| AtanBackward | 186,086 | 784 |
| SelectBackward | 91,438 | 784 |
| SubBackward | 20,124 | 787 |
| StackBackward | 19,436 | 1 |
| MmBackward | 1,404 | 1 |
| Other |  |  |
| Other | 244,358 | - |
| All | 1,201,250 | - |

In the sequential training timeline, each phase of forward and backward propagation happens one after the other, resulting in a total training time of approximately $1200ms$. These phases include data loading, fully connected forward passes, neuron charge/fire/reset operations, output computation, autograd construction, and error backpropagation. This sequential computation introduces significant delays, particularly during the repetitive neuron charging and resetting in the backward phases, which dominate the computation time.

In contrast, the parallel training timeline reduces the overall training time to around 337 ms by executing multiple forward and backward propagation phases simultaneously. This approach leverages parallel processing to handle operations such as neuron charging and threshold scanning concurrently, thereby reducing redundant delays. The comparison highlights significant efficiency gains in processes like autograd, constructing backward graphs, and error backpropagation.

## Appendix M The Statistics of Fire Rate

The spiking fire rate is derived from the top-1 test accuracy model (Average in Table 11). We extract the fire rate for each layer (12 to 16). On average, the Spatial Neuron ($\mathcal{SN}$) exhibits a higher fire rate than the PRF ($\mathcal{TN}$). The checkpoints and statistical code can be found in the open-source repository.

**Table 11: The average fire rates for each tasks.**
| Tasks | ListOps | Text | Retrieval | Image | Pathfinder |
| --- | --- | --- | --- | --- | --- |
| Avg. Fire Rate (%) | 3.53 | 4.32 | 1.48 | 3.47 | 3.29 |

**Table 12: The fire rate (%) across different layers on ListOps task.**
| Layer | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $\mathcal{TN}$ | 0.0 | 5.17 | 2.50 | 2.83 | 0.80 | 1.17 | 3.02 | 2.22 |
| $\mathcal{SN}$ | 9.60 | 5.29 | 4.51 | 2.63 | 5.58 | 3.57 | 9.57 | 5.07 |

**Table 13: The fire rate (%) across different layers on Text task.**
| Layer | 1 | 2 | 3 | 4 | 5 | 6 |
| --- | --- | --- | --- | --- | --- | --- |
| $\mathcal{TN}$ | 2.11 | 4.30 | 2.79 | 1.66 | 1.02 | 1.48 |
| $\mathcal{SN}$ | 9.20 | 9.90 | 10.11 | 7.18 | 5.55 | 5.25 |

**Table 14: The fire rate (%) across different layers on Retrieval task.**
| Layer | 1 | 2 | 3 | 4 | 5 | 6 |
| --- | --- | --- | --- | --- | --- | --- |
| $\mathcal{TN}$ | 0.66 | 0.42 | 0.42 | 0.49 | 0.48 | 0.87 |
| $\mathcal{SN}$ | 4.79 | 4.24 | 1.54 | 2.46 | 2.50 | 1.91 |

**Table 15: The fire rate (%) across different layers on Image task.**
| Layer | 1 | 2 | 3 | 4 | 5 | 6 |
| --- | --- | --- | --- | --- | --- | --- |
| $\mathcal{TN}$ | 0.22 | 7.10 | 5.25 | 3.67 | 3.90 | 4.24 |
| $\mathcal{SN}$ | 0.44 | 9.90 | 5.71 | 4.35 | 2.77 | 1.07 |

**Table 16: The fire rate (%) across different layers on Pathfinder task.**
| Layer | 1 | 2 | 3 | 4 | 5 | 6 |
| --- | --- | --- | --- | --- | --- | --- |
| $\mathcal{TN}$ | 3.37 | 3.96 | 5.07 | 4.53 | 4.53 | 4.78 |
| $\mathcal{SN}$ | 3.66 | 2.98 | 3.74 | 3.63 | 3.33 | 2.53 |

## Appendix N The Ablation Experiments

We examine the sensitivity of the $\Delta$ and $\theta$ hyper-parameters during initialization. Using a neural network architecture with layers sized 1-64-256-256-10 and a batch size of 256, we fixed $\tau=2$ and set $\Delta$ and $\theta$ as non-trainable scale values. Figure LABEL:fig.init.sensitive illustrates this sensitivity for the sMNIST (left) and psMNIST (right) datasets. The gray frames highlight a shift in sensitive regions from sMNIST to psMNIST, indicating that different data distributions require careful initialization of $\Delta$ and $\theta$. Furthermore, we investigate the impact of varying initialization of $\theta$ values on the performance of the LRA tasks, as shown in Tables 17 - 21. The suitable initialization of $\theta$ is crucial.

**Table 17: ListOps (2048) accuracy for different $\theta$ ranges initialization.**
| $\theta$ | [0, $\pi$/4] | [0, $\pi$/5] | [0, $\pi$/6] | [0, $\pi$/7] | [0, $\pi$/8] |
| --- | --- | --- | --- | --- | --- |
| Acc. (%) | 56.85 | 55.60 | 59.20 | 58.05 | 54.55 |

**Table 18: Text (4096) accuracy for different $\theta$ range initialization.**
| $\theta$ | [0, $\pi$] | [0, $\pi$/2] | [0, $\pi$/4] | [0, $\pi$/8] | [0, $\pi$/16] |
| --- | --- | --- | --- | --- | --- |
| Acc. (%) | 81.27 | 83.28 | 84.94 | 86.33 | 85.32 |

**Table 19: Retrieval (4000) accuracy for different $\theta$ range initialization.**
| $\theta$ | [0, $\pi$/4] | [0, $\pi$/5] | [0, $\pi$/6] | [0, $\pi$/7] | [0, $\pi$/8] |
| --- | --- | --- | --- | --- | --- |
| Acc. (%) | 89.64 | 89.52 | 89.75 | 89.88 | 89.72 |

**Table 20: Image (1024) accuracy for different $\theta$ range initialization.**
| $\theta$ | [0, 2$\pi$/0.15] | [0, 2$\pi$/0.2] | [0, 2$\pi$/0.25] | [0, 2$\pi$/0.3] | [0, 2$\pi$/0.35] |
| --- | --- | --- | --- | --- | --- |
| Acc. (%) | 84.36 | 84.43 | 84.77 | 84.54 | 84.38 |

**Table 21: Pathfinder (1024) accuracy for different $\theta$ range initialization.**
| $\theta$ | [0, 2$\pi$] | [0, 2$\pi$/0.8] | [0, 2$\pi$/0.6] | [0, 2$\pi$/0.4] | [0, 2$\pi$/0.2] |
| --- | --- | --- | --- | --- | --- |
| Acc. (%) | 91.14 | 91.19 | 91.76 | 89.87 | 90.97 |

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: * Any idea why Li et al 2024 does better than you on psMNIST, but not on LRA?

* L244 "At this point, all recursive forms have been converted into dynamic equations." I agree, but still, the dynamic equations should presumably be solved sequentially. So, I don't see how this helps parallelization. More insight is needed here.

* According to Fig 2 the residual connections carry floats, not spikes. This could be a pb for implementation on neuromorphic chips. Could the authors elaborate?

* L505 "This energy reduction mainly benefits from the extreme sparsity of spikes"
Could even higher levels of sparsity could be encouraged via some additional terms in the loss function?
soundness: 3
strengths: The accuracy on the LRA benchmark is trully impressive, almost as good as S4, while it could consume much less energy (at least in theory).
summary: The authors introduce a new spiking neuron model, the parallel resonate and fire (PRF) neuron, that has two advantages over the classic leaky integrate-and-fire (LIF) neuron.

1) it is parallelizable and thus speeds up training and inference times

2) it has a longer-term memory because it can resonate.
weaknesses: * The accuracy on spMNIST is much less convicing. Here two similar recent proposals are much more accurate: Li et al 2024 (cited) and Bal et al 2024 (cited)
These refs should be added to Table 2.
I also suggest adding sequential CIFAR, which is more challenging.

* Other parallel spiking neuron models should be cited and added to Fig 4 left and to the comparison tables when possible, and the differences/similarities should be highlighted:  Fang et al 2023 (cited), Chen et al 2024, Yarga et al 2024

* The resonance seems similar to Baronig et al. 2024, which should be cited, and the differences/similarities should be highlighted.

References:

Baronig M, Ferrand R, Sabathiel S, Legenstein R (2024) Advancing Spatio-Temporal Processing in Spiking Neural Networks through Adaptation. arXiv Available at: http://arxiv.org/abs/2408.07517.

Fang W, Yu Z, Zhou Z, Chen D, Chen Y, Ma Z, Masquelier T, Tian Y (2023) Parallel Spiking Neurons with High Efficiency and Ability to Learn Long-term Dependencies. NeurIPS 1:1–13 Available at: http://arxiv.org/abs/2304.12760.

Zeyu Liu, Gourav Datta, Anni Li, and Peter Anthony Beerel. Lmuformer: Low complexity yet powerful spiking model with legendre memory units. arXiv preprint arXiv:2402.04882, 2024.

Bal M, Sengupta A (2024) P-SpikeSSM: Harnessing Probabilistic Spiking State Space Models for Long-Range Dependency Tasks. arXiv:1–16 Available at: http://arxiv.org/abs/2406.02923.

Chen X, Wu J, Ma C, Yan Y, Wu Y, Tan KC (2024) PMSN: A Parallel Multi-compartment Spiking Neuron for Multi-scale Temporal Processing. arXiv Available at: http://arxiv.org/abs/2408.14917.

Yarga SYA, Wood SUN (2024) Accelerating Spiking Neural Networks with Parallelizable Leaky Integrate-and-Fire Neurons. TechRxiv:1–10.

--

POST REBUTTAL:

All my main concerns have been addressed, so I'm raising my score to 6. I don't go higher because similar parallel neurons have been proposed recently and because I still think the membrane potential (float) residual connections could be an issue for some neuromorphic chips. Still, I think this paper should be accepted.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: - What are the connections between the proposed PRF neuron and SSMs?

- What impact does the reset have on the performance of the PRF model in LRA tasks?

- The existing PSN and PMSN models also aim to address the challenges of parallel training and modeling long-range dependencies. Please discuss the differences between PRF and these models.

[1] Parallel Spiking Neurons with High Efficiency and Ability to Learn Long-term Dependencies.

[2] PMSN: A Parallel Multi-compartment Spiking Neuron for Multi-scale Temporal Processing.
soundness: 3
strengths: - The paper is well-written and easy to read, with excellent presentation. The proposed methods are novel.

- The experiments and analyses are comprehensive. The performance improvements on the challenging LRA tasks in terms of both efficiency and accuracy are significant.

- The PRF neuron shows promise for deployment on neuromorphic chips for inference.
summary: This paper focuses on SNNs for long sequence learning. It first reformulates the reset mechanism of the LIF model into an iterative form to enable parallel training. The authors then propose the PRF neuron to enhance the modeling of long-range dependencies. Both theoretical and experimental results demonstrate the effectiveness and efficiency of the proposed methods.
weaknesses: - The reset mechanism cannot be computed in parallel with the membrane potential integration process, indicating that the overall neuron model is not fully trained in parallel.

- There is a lack of comparison with existing parallel spiking neuron models, such as the PSN family models [1] and PMSN [2] models.


-----
After the discussion, all my concerns have been addressed, and I have raised my score to 8.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: 1. To my knowledge, removing the reset mechanism might be the simplest and most effective approach. Several existing studies have explored parallel implementations, such as PSN[1]. The authors should include an experiment to compare their method with these approaches.
2. The paper validates the proposed method on the S/PS-MNIST dataset and claims to achieve SOTA results. However, the performance on other image sequence datasets is not significant. I suggest that the authors expand the experiments further and compare the proposed method with other advanced sequence models to better demonstrate its advantages.

[1] Fang W, Yu Z, Zhou Z, et al. Parallel spiking neurons with high efficiency and ability to learn long-term dependencies[J]. Advances in Neural Information Processing Systems, 2024, 36.
soundness: 3
strengths: 1. This paper proposes a decoupled reset-based parallel training strategy, reducing the training time complexity of LIF neurons to O(LlogL). The authors claim this is the first application of the reset mechanism in parallel training.

2. The paper introduces a parallel mechanism into Resonate-and-Fire neurons, enhancing the model's temporal processing capability.

3. The paper presents the SD-TCM architecture, which, combined with PRF neurons, was tested on LRA tasks. The authors claim it achieved performance comparable to SSM (S4).
summary: This paper proposes a decoupled reset method for parallel training of spiking neurons, significantly improving the training speed of LIF neurons. Additionally, a parallel mechanism is introduced into the Resonate-and-Fire model, combined with the SD-TCM structure to achieve efficient long-sequence learning. This approach significantly reduces model energy consumption while achieving comparable performance to SSM (S4) and outperforming Transformer in LRA tasks. However, the experimental section is incomplete, particularly the ablation studies. Furthermore, the motivation behind the reset mechanism is not clearly described. Additional experiments are recommended to further validate the proposed approach.
weaknesses: 1. The proposed method's long-sequence learning capability is mainly validated in Section 5.2. However, this validation is limited to the S/PS-MNIST image sequence task, with only marginal improvements compared to current SOTA models.
2. The study lacks effective ablation experiments. In Section 5.3, the model's performance on LRA tasks is evaluated, but the contribution of the SD-TCM structure to the performance was not verified, and the improvement was solely attributed to the temporal processing capability of RF neurons.
3. The paper combines Figures 1(d) and (e) to illustrate the strong temporal processing capability of RF neurons. However, the experimental section lacks relevant visualizations to further explain why RF neurons are more suitable for temporal tasks compared to LIF neurons.
4. The research motivation is unclear. The abstract claims that this is the first implementation of a parallel reset mechanism, yet the experimental section does not evaluate its impact on performance. 
5. The reproducibility of the experiments is limited. The experimental section should include multiple trials to further validate the proposed method's effectiveness.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: see weakness bellow.
soundness: 2
strengths: 1. This paper integrate PRF neurons into the SD-TCM module, which performs comparable to S4 in long-term dependency tasks while reducing inference energy consumption by two orders of magnitude.
2. Improved the previous SNN parallel method to achieve acceleration, and achieved good results on the lra dataset.
summary: This paper proposes a new decoupled reset method to accelerate training, making the effect of parallel training equivalent to sequential training. This method can increase the back propagation speed by three orders of magnitude and is applicable to all types of spiking neurons. In addition, in order to effectively capture long-term dependencies, this paper designs a new type of PRF neuron, which has an oscillating membrane potential in the complex domain and works through an adaptive differentiable reset mechanism.
weaknesses: 1. Since there is no code in the appendix of this article, I have doubts about the authenticity of lra's experimental results.
2. There are too many typos in the article. The citation format of the whole article should be “\citep{}” instead of  “\cite{}”. At the same time, please use the same “spike-driven” format throughout the article.
3. The comparison in Figure 1 is unfair, as SSM is a block and the proposed PRF is just a neuron. This makes the complexity comparison unfair.
4. In addition, I have doubts about the effectiveness of the proposed method. In Table 4, we can see that the improvement gains of different modules are very low. Does this mean that the proposed method relies heavily on its own architecture? Please add additional ablation experiments on other architectures.

Meta-review: This paper proposes a parallel resonate and fire (PRF) neuron model aimed at efficient long-sequence learning. The main claims are: 1) A decoupled reset method that reduces training time from O(L²) to O(LlogL), 2) A PRF neuron design that leverages oscillating membrane potential for capturing long-range dependencies, and 3) A spike-driven architecture using PRF that achieves performance comparable to S4 while reducing energy consumption by two orders of magnitude.

Positive:

- Novel approach to parallel training in SNNs through the decoupled reset

Weaknesses:

- reproducibility concerns due to lack of code availability at submission time
- limited exploration of PRF neuron effectiveness across different architectures
- overlooked the most relevant state-of-the-art methods to compare against: Liquid-S4 (https://arxiv.org/abs/2209.12951) and the entire line of Liquid time-constant networks (https://arxiv.org/abs/2006.04439) work, that shares high degrees of similarities in scope with this model.

I vote for the rejection of the paper.

Decision: Reject