# Exploiting Heterogeneity in Timescales for Sparse Recurrent Spiking Neural
Networks for Energy-Efficient Edge Computing

Biswadeep Chakraborty  Department of Electrical and Computer Engineering
Georgia Institute of Technology
Atlanta, USA
<EMAIL_ADDRESS>Saibal Mukhopadhyay  Department of Electrical and
Computer Engineering
Georgia Institute of Technology
Atlanta, USA
<EMAIL_ADDRESS>
###### Abstract

Spiking Neural Networks (SNNs) represent the forefront of neuromorphic
computing, promising energy-efficient and biologically plausible models for
complex tasks. This paper weaves together three groundbreaking studies that
revolutionize SNN performance through the introduction of heterogeneity in
neuron and synapse dynamics. We explore the transformative impact of
Heterogeneous Recurrent Spiking Neural Networks (HRSNNs), supported by
rigorous analytical frameworks and novel pruning methods like Lyapunov Noise
Pruning (LNP). Our findings reveal how heterogeneity not only enhances
classification performance but also reduces spiking activity, leading to more
efficient and robust networks. By bridging theoretical insights with practical
applications, this comprehensive summary highlights the potential of SNNs to
outperform traditional neural networks while maintaining lower computational
costs. Join us on a journey through the cutting-edge advancements that pave
the way for the future of intelligent, energy-efficient neural computing.

###### Index Terms:

component, formatting, style, styling, insert

##  I Introduction

Spiking Neural Networks (SNNs), often referred to as the third generation of
neural networks, have garnered significant attention due to their potential
for lower operating power when mapped to hardware. SNNs, particularly those
using leaky integrate-and-fire (LIF) neurons, have demonstrated classification
performance comparable to deep neural networks (DNNs). However, the majority
of these models rely on supervised training algorithms like backpropagation-
through-time (BPTT) [1, 2, 3], which are highly data-dependent and struggle
with limited training data and generalization [4, 5]. Moreover, BPTT-trained
models require complex architectures with a large number of neurons to achieve
good performance. While unsupervised learning methods such as Spike-Timing
Dependent Plasticity (STDP) have been introduced, they typically underperform
compared to their supervised counterparts due to the high complexity of
training [6].

To address these challenges, this paper introduces a series of studies focused
on enhancing SNN performance through heterogeneity in network dynamics. The
first paper presents a Heterogeneous Recurrent Spiking Neural Network (HRSNN)
with variability in both LIF neuron parameters and STDP dynamics [7]. Previous
works have suggested that heterogeneity in neuron time constants can improve
classification performance [8, 9, 10, 11], but lacked theoretical explanations
for these improvements. Our study not only explores how heterogeneity in both
neuronal and synaptic parameters can enhance performance with less training
data and fewer connections but also leverages a novel Bayesian Optimization
(BO) method for hyperparameter tuning. This approach scales well for larger,
more complex tasks, such as action recognition, that were previously
infeasible.

The second paper builds on these insights by providing rigorous analytical
results that explain the effects of heterogeneity in LIF and STDP parameters
[12]. We demonstrate through mathematical modeling and empirical validation
that heterogeneity improves the linear separation properties of unsupervised
SNN models. This paper establishes a robust theoretical framework, linking the
observed performance enhancements to fundamental properties of heterogeneous
networks. We show that optimizing neuronal and synaptic heterogeneity can
reduce spiking activity while improving memory capacity, defined as the
network’s ability to learn and retain information [13, 14].

Despite the benefits of heterogeneity, the increased computational cost
remains a significant concern. The third paper addresses this by introducing
Lyapunov Noise Pruning (LNP), a novel method that exploits the heterogeneity
to prune the network efficiently [15]. By identifying and removing redundant
or noisy connections, LNP maintains high performance while significantly
reducing computational demands. Grounded in Lyapunov stability theory, this
task-agnostic pruning strategy ensures the pruned network remains stable and
effective across various applications without the need for extensive
retraining.

The key contributions of these three papers are as follows:

  * •

Heterogeneous Recurrent Spiking Neural Networks (HRSNNs): Introduction of
heterogeneity in LIF and STDP dynamics, demonstrating enhanced robustness and
adaptability of SNNs [7].

  * •

Analytical Insights into HRSNNs: Detailed mathematical modeling and empirical
verification of the effects of heterogeneity, providing guidelines for
optimizing SNNs [12].

  * •

Lyapunov Noise Pruning (LNP): Development of a novel, task-agnostic pruning
method leveraging heterogeneity to reduce computational costs while
maintaining network performance [15].

The findings from these studies underscore the transformative potential of
incorporating heterogeneity in SNNs. By enhancing robustness, improving
learning dynamics, and mitigating computational costs, these advancements pave
the way for more practical and powerful neural network models. This
consolidated summary aims to provide researchers and practitioners with a
comprehensive overview of these key developments, fostering further innovation
and application in the field of SNNs.

##  II Related Works

The field of Spiking Neural Networks (SNNs) is enriched by various
foundational and contemporary studies. SNNs leverage the temporal dynamics of
spiking neurons, which has been a significant focus in neuroscience and
machine learning. Early work on SNNs by [16] established the fundamental
concepts of spiking neurons and their bio-inspired learning mechanisms, such
as Spike-Timing Dependent Plasticity (STDP) [17][18, 19].

Empirical results have shown the effectiveness of SNNs in tasks like
spatiotemporal data classification [20, 21], sequence-to-sequence mapping
[22], object detection [23, 24], and universal function approximation [25,
26]. The energy efficiency of SNNs, attributed to their sparse firing
activity, is a critical advantage highlighted in studies by [27], [28], and
[29]. However, optimizing spiking activity while maintaining performance is a
key challenge addressed by recent empirical and analytical studies [30, 31].

The concept of heterogeneity in SNNs has been explored to improve performance
and efficiency. Studies by [32], [33], and [34] have shown that biological
systems naturally employ heterogeneous dynamics to enhance robustness and
adaptability. In the context of SNNs, [8] and [10] demonstrated that
heterogeneity in neuronal time constants could improve classification
performance, though they lacked a theoretical foundation. Our work builds on
these insights, providing both empirical evidence and theoretical explanations
for the benefits of heterogeneity.

Recent advancements in optimizing SNNs have been marked by the development of
novel pruning methods. The Lyapunov Noise Pruning (LNP) method introduced by
[15] represents a significant step forward. LNP leverages heterogeneity to
identify and prune redundant connections, thereby reducing computational
complexity while maintaining network stability and performance. This approach
contrasts with traditional activity-based pruning methods and highlights the
importance of theoretical grounding in developing efficient neural network
models.

In addition to these foundational works, studies on the practical applications
of SNNs in various domains underscore their versatility. [35] and [36]
explored the use of SNNs in brain-computer interfacing and EEG electrode
optimization, respectively, demonstrating the potential for SNNs in biomedical
applications. [23] and [37] further highlighted the application of SNNs in
object detection and edge intelligence.

The integration of heterogeneous dynamics in SNNs for improved performance and
energy efficiency has also been explored by [7] and [38], providing a robust
framework for future research in this area. The intersection of theoretical
insights and practical applications positions SNNs as a promising technology
for next-generation neural network models.

##  III Methods

###  III-A Recurrent Spiking Neural Network

SNN consists of spiking neurons connected with synapses. The spiking LIF is
defined by the following equations:

$\tau_{m}\frac{dv}{dt}=a+R_{m}I-v;v=v_{\text{reset }},\text{ if
}v>v_{\text{threshold }}$ (1)

where $R_{m}$ is membrane resistance, $\tau_{m}=R_{m}C_{m}$ is time constant
and $C_{m}$ is membrane capacitance. $a$ is the resting potential. $I$ is the
sum of current from all input synapses connected to the neuron. A spike is
generated when membrane potential $v$ crosses the threshold, and the neuron
enters refractory period $r$, during which the neuron maintains its membrane
potential at $v_{\text{reset }}.$We construct the HRSNN from the baseline
recurrent spiking network (RSNN) consisting of three layers: (1) an input
encoding layer ($\mathcal{I}$), (2) a recurrent spiking layer ($\mathcal{R}$),
and (3) an output decoding layer ($\mathcal{O}$). The recurrent layer consists
of excitatory and inhibitory neurons, distributed in a ratio of
$N_{E}:N_{I}=4:1$. The PSPs of post-synaptic neurons produced by the
excitatory neurons are positive, while those produced by the inhibitory
neurons are negative. We used a biologically plausible LIF neuron model and
trained the model using STDP rules.

From here on, we refer to connections between $\mathcal{I}$ and $\mathcal{R}$
neurons as $\mathcal{S}_{\mathcal{IR}}$ connections, inter-recurrent layer
connections as $\mathcal{S}_{\mathcal{RR}}$, and $\mathcal{R}$ to
$\mathcal{O}$ as $\mathcal{S}_{\mathcal{RO}}$. We created
$\mathcal{S}_{\mathcal{RR}}$ connections using probabilities based on
Euclidean distance, $D(i,j)$, between any two neurons $i,j$:

$P(i,j)=C\cdot\exp\left(-\left(\frac{D(i,j)}{\lambda}\right)^{2}\right)$ (2)

with closer neurons having higher connection probability. Parameters $C$ and
$\lambda$ set the amplitude and horizontal shift, respectively, of the
probability distribution. $\mathcal{I}$ contains excitatory encoding neurons,
which convert input data into spike trains. $\mathcal{S}_{IR}$ only randomly
chooses $30\%$ of the excitatory and inhibitory neurons in $\mathcal{R}$ as
the post-synaptic neuron. The connection probability between the encoding
neurons and neurons in the $\mathcal{R}$ is defined by a uniform probability
$\mathcal{P}_{\mathcal{IR}}$, which, together with $\lambda$, will be used to
encode the architecture of the HRSNN and optimized using BO. In this work,
each neuron received projections from some randomly selected neurons in
$\mathcal{R}$.

We used unsupervised, local learning to the spiking recurrent model by letting
STDP change each $\mathcal{S}_{\mathcal{RR}}$ and $\mathcal{S}_{\mathcal{IR}}$
connection, modeled as:

$\frac{dW}{dt}=A_{+}T_{pre}\sum_{o}\delta\left(t-t_{\text{post
}}^{o}\right)-A_{-}T_{\text{post }}\sum_{i}\delta\left(t-t_{\text{pre
}}^{i}\right)$ (3)

where $A_{+},A_{-}$ are the potentiation/depression learning rates and
$T_{\text{pre }}/T_{\text{post }}$ are the pre/post-synaptic trace variables,
modeled as,

$\displaystyle\tau_{+}^{*}\frac{dT_{\text{pre }}}{dt}$
$\displaystyle=-T_{\text{pre }}+a_{+}\sum_{i}\delta\left(t-t_{\text{pre
}}^{i}\right)$ (4) $\displaystyle\tau_{-}^{*}\frac{dT_{\text{post }}}{dt}$
$\displaystyle=-T_{\text{post }}+a_{-}\sum_{o}\delta\left(t-t_{\text{post
}}^{o}\right)$ (5)

where $a_{+},a_{-}$ are the discrete contributions of each spike to the trace
variable, $\tau_{+}^{*},\tau_{-}^{*}$ are the decay time constants,
$t_{\text{pre }}^{i}$ and $t_{\text{post }}^{o}$ are the times of the pre-
synaptic and post-synaptic spikes, respectively.

Figure 1: Concept of HRSNN with variable Neuronal and Synaptic Dynamics

Heterogeneous LIF Neurons The use of multiple timescales in spiking neural
networks has several underlying benefits, like increasing the memory capacity
of the network. In this paper, we propose the usage of heterogeneous LIF
neurons with different membrane time constants and threshold voltages, thereby
giving rise to multiple timescales. Due to differential effects of excitatory
and inhibitory heterogeneity on the gain and asynchronous state of sparse
cortical networks [39], [40], we use different gamma distributions for both
the excitatory and inhibitory LIF neurons. This is also inspired by the
brain’s biological observations, where the time constants for excitatory
neurons are larger than the time constants for the inhibitory neurons. Thus,
we incorporate the heterogeneity in our Recurrent Spiking Neural Network by
using different membrane time constants $\tau$ for each LIF neuron in
$\mathcal{R}$. This gives rise to a distribution for the time constants of the
LIF neurons in $\mathcal{R}$.

Heterogeneous STDP Experiments on different brain regions and diverse neuronal
types have revealed a wide variety of STDP forms that vary in plasticity
direction, temporal dependence, and the involvement of signaling pathways
[41], [42], [43]. As described by Pool et al. [44], one of the most striking
aspects of this plasticity mechanism in synaptic efficacy is that the STDP
windows display a great variety of forms in different parts of the nervous
system. However, most STDP models used in Spiking Neural Networks are
homogeneous with uniform timescale distribution. Thus, we explore the
advantages of using heterogeneities in several hyperparameters discussed
above. This paper considers heterogeneity in the scaling function constants
($A_{+},A_{-}$) and the decay time constants ($\tau_{+},\tau_{-}$).

###  III-B Classification Property of HRSNN

We theoretically compare the performance of the heterogeneous spiking
recurrent model with its homogeneous counterpart using a binary classification
problem. The ability of HRSNN to distinguish between many inputs is studied
through the lens of the edge-of-chaos dynamics of the spiking recurrent neural
network, similar to the case in spiking reservoirs shown by Legenstein et al.
[45]. Also, $\mathcal{R}$ possesses a fading memory due to its short-term
synaptic plasticity and recurrent connectivity. For each stimulus, the final
state of the $\mathcal{R}$, i.e., the state at the end of each stimulus,
carries the most information. The authors showed that the rank of the final
state matrix $F$ reflects the separation property of a kernel:
$F=\left[S(1)\quad S(2)\quad\cdots\quad S(N)\right]^{T}$ where $S(i)$ is the
final state vector of $\mathcal{R}$ for the stimulus $i$. Each element of $F$
represents one neuron’s response to all the $N$ stimuli. A higher rank in $F$
indicates better kernel separation if all $N$ inputs are from $N$ distinct
classes.

The effective rank is calculated using Singular Value Decomposition (SVD) on
$F$, and then taking the number of singular values that contain $99\%$ of the
sum in the diagonal matrix as the rank. i.e. $F=U\Sigma V^{T}$ where $U$ and
$V$ are unitary matrices, and $\Sigma$ is a diagonal matrix
$\operatorname{diag}\left(\lambda_{1},\lambda_{2},\lambda_{3},\ldots,\lambda_{N}\right)$
that contains non-negative singular values such that
$(\lambda_{1}\geq\lambda_{2}\cdots\geq\lambda_{N})$.

Figure 2: An illustrative example showing the Heterogeneous Recurrent Spiking
Neural Network structure. First, we show the temporal encoding method based on
the sensory receptors receiving the difference between two time-adjacent data.
Next, the input sequences are encoded by the encoding neurons that inject the
spike train into $30\%$ neurons in $\mathcal{R}.\mathcal{R}$ contains a $4:1$
ratio of excitatory (green nodes) and inhibitory (orange nodes), where the
neuron parameters are heterogeneous. The synapses are trained using the
heterogeneous STDP method.

Definition: Linear separation property of a neuronal circuit $\mathcal{C}$ for
$m$ different inputs $u_{1},\ldots,u_{m}(t)$ is defined as the rank of the
$n\times m$ matrix $M$ whose columns are the final circuit states
$\mathbf{x}_{u_{i}}\left(t_{0}\right)$ obtained at time $t_{0}$ for the
preceding input stream $u_{i}$.

Following from the definition introduced by Legenstein et al. [45], if the
rank of the matrix $M=m$, then for the inputs $u_{i}$, any given assignment of
target outputs $y_{i}\in\mathbb{R}$ at time $t_{0}$ can be implemented by
$\mathcal{C}$.

We use the rank of the matrix as a measure for the linear separation of a
circuit $C$ for distinct inputs. This leverages the complexity and diversity
of nonlinear operations carried out by $C$ on its input to boost the
classification performance of a subsequent linear decision-hyperplane.

Theorem 1:  Assuming $\mathcal{S}_{u}$ is finite and contains $s$ inputs, let
$r_{\text{Hom}},r_{\text{Het}}$ are the ranks of the $n\times s$ matrices
consisting of the $s$ vectors $\mathbf{x}_{u}\left(t_{0}\right)$ for all
inputs $u$ in $\mathcal{S}_{u}$ for each of Homogeneous and Heterogeneous
RSNNs respectively. Then $r_{\text{Hom}}\leq r_{\text{Het}}$.

Short Proof: Let us fix some inputs $u_{1},\ldots,u_{r}$ in $\mathcal{S}_{u}$
so that the resulting $r$ circuit states
$\mathbf{x}_{u_{i}}\left(t_{0}\right)$ are linearly independent. Using the
Eckart-Young-Mirsky theorem for low-rank approximation, we show that the
number of linearly independent vectors for HeNHeS is greater than or equal to
the number of linearly independent vectors for HoNHoS. The detailed proof is
given in the Supplementary.

Definition : Given $K_{\rho}$ is the modified Bessel function of the second
kind, and $\sigma^{2},\kappa,\rho$ are the variance, length scale, and
smoothness parameters respectively, we define the modified Matern kernel on
the Wasserstein metric space $\mathcal{W}$ between two distributions
$\mathcal{X},\mathcal{X}^{\prime}$ given as

$k\left(\mathcal{X},\mathcal{X}^{\prime}\right)=\sigma^{2}\frac{2^{1-\rho}}{\Gamma(\rho)}\left(\sqrt{2\rho}\frac{\mathcal{W}(\mathcal{X},\mathcal{X}^{\prime})}{\kappa}\right)^{\rho}H_{\rho}\left(\sqrt{2\rho}\frac{\mathcal{}(\mathcal{X},\mathcal{X}^{\prime})}{\kappa}\right)$
(6)

where $\Gamma(.),H(.)$ is the Gamma and Bessel function, respectively.

Theorem 2: The modified Matern function on the Wasserstein metric space
$\mathcal{W}$ is a valid kernel function

Short Proof:  To show that the above function is a kernel function, we need to
prove that Mercer’s theorem holds. i.e., (i) the function is symmetric and
(ii) in finite input space, the Gram matrix of the kernel function is positive
semi-definite. The detailed proof is given in the Supplementary.

###  III-C Optimal Hyperparameter Selection using Bayesian Optimization

While BO is used in various settings, successful applications are often
limited to low-dimensional problems, with fewer than twenty dimensions [46].
Thus, using BO for high-dimensional problems remains a significant challenge.
In our case of optimizing HRSNN model parameters for 2000, we need to optimize
a huge number of parameters, which is extremely difficult for BO. As discussed
by Eriksson et al. [47], suitable function priors are especially important for
good performance. Thus, we used a biologically inspired initialization of the
hyperparameters derived from the human brain (see Supplementary for details).

This paper uses a modified BO to estimate parameter distributions for the LIF
neurons and the STDP dynamics. To learn the probability distribution of the
data, we modify the surrogate model and the acquisition function of the BO to
treat the parameter distributions instead of individual variables. This makes
our modified BO highly scalable over all the variables (dimensions) used. The
loss for the surrogate model’s update is calculated using the Wasserstein
distance between the parameter distributions.

BO uses a Gaussian process to model the distribution of an objective function
and an acquisition function to decide points to evaluate. For data points in a
target dataset $x\in X$ and the corresponding label $y\in Y$, an SNN with
network structure $\mathcal{V}$ and neuron parameters $\mathcal{W}$ acts as a
function $f_{\mathcal{V},\mathcal{W}}(x)$ that maps input data $x$ to
predicted label $\tilde{y}$. The optimization problem in this work is defined
as

$\min_{\mathcal{V},\mathcal{W}}\sum_{x\in X,y\in
Y}\mathcal{L}\left(y,f_{\mathcal{V},\mathcal{W}}(x)\right)$ (7)

where $\mathcal{V}$ is the set of hyperparameters of the neurons in
$\mathcal{R}$ (Details of hyperparameters given in the Supplementary) and
$\mathcal{W}$ is the multi-variate distribution constituting the distributions
of (i) the membrane time constants $\tau_{m-E},\tau_{m-I}$ of the LIF neurons,
(ii) the scaling function constants $(A_{+},A_{-})$ and (iii) the decay time
constants $\tau_{+},\tau_{-}$ for the STDP learning rule in
$\mathcal{S}_{\mathcal{RR}}$.

Again, BO needs a prior distribution of the objective function $f(\vec{x})$ on
the given data
$\mathcal{D}_{1:k}=\left\\{\vec{x}_{1:k},f\left(\vec{x}_{1:k}\right)\right\\}.$
In GP-based BO, it is assumed that the prior distribution of
$f\left(\vec{x}_{1:k}\right)$ follows the multivariate Gaussian distribution,
which follows a Gaussian Process with mean $\vec{\mu}_{\mathcal{D}_{1:k}}$ and
covariance $\vec{\Sigma}_{\mathcal{D}_{1:k}}$. We estimate
$\vec{\Sigma}_{\mathcal{D}_{1:k}}$ using the modified Matern kernel function,
which is given in Eq. 6. In this paper, we use $d(x,x^{\prime})$ as the
Wasserstein distance between the multivariate distributions of the different
parameters. It is to be noted here that for higher-dimensional metric spaces,
we use the Sinkhorn distance as a regularized version of the Wasserstein
distance to approximate the Wasserstein distance [48].

$\mathcal{D}_{1:k}$ are the points that have been evaluated by the objective
function, and the GP will estimate the mean $\vec{\mu}_{\mathcal{D}_{k:n}}$
and variance $\vec{\sigma}_{\mathcal{D}_{k:n}}$ for the rest unevaluated data
$\mathcal{D}_{k:n}$. The acquisition function used in this work is the
expected improvement (EI) of the prediction fitness as:

$EI\left(\vec{x}_{k:n}\right)=\left(\vec{\mu}_{\mathcal{D}_{k:n}}-f\left(x_{\text{best
}}\right)\right)\Phi(\vec{Z})+\vec{\sigma}_{\mathcal{D}_{k:n}}\phi(\vec{Z})$
(8)

where $\Phi(\cdot)$ and $\phi(\cdot)$ denote the probability distribution
function and the cumulative distribution function of the prior distributions,
respectively. $f\left(x_{\text{best }}\right)=\max
f\left(\vec{x}_{1:k}\right)$ is the maximum value that has been evaluated by
the original function $f$ in all evaluated data $\mathcal{D}_{1:k}$ and
$\vec{Z}=\frac{\vec{\mu}_{\mathcal{D}_{k:n}}-f\left(x_{\text{best
}}\right)}{\vec{\sigma}_{\mathcal{D}_{k:n}}}$. BO will choose the data
$x_{j}=$
$\operatorname{argmax}\left\\{EI\left(\vec{x}_{k:n}\right);x_{j}\subseteq\vec{x}_{k:n}\right\\}$
as the next point to be evaluated using the original objective function.

##  IV Analytical Results

###  IV-A Preliminaries and Definitions

Heterogeneity:

We define heterogeneity as a measure of the variability of the hyperparameters
in an RSNN that gives rise to an ensemble of neuronal dynamics.

Entropy is used to measure population diversity. Assuming that the random
variable for the hyperparameters $X$ follows a multivariate Gaussian
Distribution ($X\sim\mathcal{N}(\mu,\Sigma)$), then the differential entropy
of $x$ on the multivariate Gaussian distribution, is
$\displaystyle\mathcal{H}(x)=\frac{n}{2}\ln(2\pi)+\frac{1}{2}\ln|\Sigma|+\frac{n}{2}$.
Now, if we take any density function $q(\mathrm{x})$ that satisfies
$\displaystyle\int q(\mathrm{x})x_{i}x_{j}d\mathrm{x}=\Sigma_{ij}$ and
$p=\mathcal{N}(0,\Sigma),$ then $\mathcal{H}(q)\leq\mathcal{H}(p).$ The
Gaussian distribution maximizes the entropy for a given covariance. Hence, the
log-determinant of the covariance matrix bounds entropy. Thus, for the rest of
the paper, we use the determinant of the covariance matrix to measure the
heterogeneity of the network.

Memory Capacity: Given an input signal x(t), the memory capacity $\mathcal{C}$
of a trained RSNN model is defined as a measure for the ability of the model
to store and recall previous inputs fed into the network [49, 50].
In this paper, we use $\mathcal{C}$ as a measure of the performance of the
model, which is based on the network’s ability to retrieve past information
(for various delays) from the reservoir using the linear combinations of
reservoir unit activations observed at the output. Intuitively, HRSNN can be
interpreted as a set of coupled filters that extract features from the input
signal. The final readout selects the right combination of those features for
classification or prediction. First, the $\tau$-delay $\mathcal{C}$ measures
the performance of the $\mathrm{RC}$ for the task of reconstructing the
delayed version of model input $x(t)$ at delay $\tau$ (i.e., $x(t-\tau)$ ) and
is defined as the squared correlation coefficient between the desired output (
$\tau$-time-step delayed input signal, $x(t-\tau))$ and the observed network
output $y_{\tau}(t)$, given as:

$\mathcal{C}=\lim_{\tau_{\max}\rightarrow\infty}\sum_{\tau=1}^{\tau_{\max}}\mathcal{C}(\tau)=\lim_{\tau_{\max}\rightarrow\infty}\sum_{\tau=1}^{\tau_{\max}}\frac{\operatorname{Cov}^{2}\left(x(t-\tau),y_{\tau}(t)\right)}{\operatorname{Var}(x(t))\operatorname{Var}\left(y_{\tau}(t)\right)},\tau\in\mathbb{N},$
(9)

where $\operatorname{Cov}(\cdot)$ and $\operatorname{Var}(\cdot)$ denote the
covariance function and variance function, respectively. The $y_{\tau}(t)$ is
the model output in this reconstruction task. $\mathcal{C}$ measures the
ability of RSNN to reconstruct precisely the past information of the model
input. Thus, increasing $\mathcal{C}$ indicates the network is capable of
learning a greater number of past input patterns, which in turn, helps in
increasing the performance of the model. For the simulations, we use
$\tau_{\max}=100$.

Spike-Efficiency:  Given an input signal x(t), the spike-efficiency
($\mathcal{E}$) of a trained RSNN model is defined as the ratio of the memory
capacity $\mathcal{C}$ to the average total spike count per neuron
$\tilde{S}$.

$\mathcal{E}$ is an analytical measure used to compare how $\mathcal{C}$ and
hence the model’s performance is improved with per unit spike activity in the
model. Ideally, we want to design a system with high $\mathcal{C}$ using fewer
spikes. Hence we define $\mathcal{E}$ as the ratio of the memory capacity
using $N_{\mathcal{R}}$ neurons $\mathcal{C}(N_{\mathcal{R}})$ to the average
number of spike activations per neuron ($\tilde{S}$) and is given as:

$\mathcal{E}=\frac{\mathcal{C}(N_{\mathcal{R}})}{\frac{\sum_{i=1}^{N_{\mathcal{R}}}S_{i}}{N_{\mathcal{R}}}},\quad\quad
S_{i}=\int_{0}^{T}s_{i}(t)dt\approx N_{\text{post
}}\frac{T}{\int_{t_{ref}}^{\infty}t\Phi_{i}dt}$ (10)

where $N_{\text{post }}$ is the number of postsynaptic neurons, $\Phi_{i}$ is
the inter-spike interval spike frequency for neuron $i$, and $T$ is the total
time. It is to be noted here that the total spike count $S$ is obtained by
counting the total number of spikes in all the neurons in the recurrent layer
until the emission of the first spike at the readout layer.

We present three main analytical findings. Firstly, neuronal dynamic
heterogeneity increases memory capacity by capturing more principal components
from the input space, leading to better performance and improved
$\mathcal{C}$. Secondly, STDP dynamic heterogeneity decreases spike activation
without affecting $\mathcal{C}$, providing better orthogonalization among the
recurrent network states and a more efficient representation of the input
space, lowering higher-order correlation in spike trains. This makes the model
more spike-efficient since the higher-order correlation progressively
decreases the information available through neural population [51, 52].
Finally, incorporating heterogeneity in both neuron and STDP dynamics boosts
the $\mathcal{C}$ to spike activity ratio, i.e., $\mathcal{E}$, which enhances
performance while reducing spike counts.

Memory Capacity: The performance of an RSNN depends on its ability to retain
the memory of previous inputs. To quantify the relationship between the
recurrent layer dynamics and $\mathcal{C}$, we note that extracting
information from the recurrent layer is made using a combination of the
neuronal states. Hence, more linearly independent neurons would offer more
variable states and, thus, more extended memory.

Lemma 3.1.1: The state of the neuron can be written as follows: $\displaystyle
r_{i}(t)=\sum_{k=0}^{N_{R}}\sum_{n=1}^{N_{R}}\lambda_{n}^{k}\left\langle
v_{n}^{-1},\mathbf{w}^{\mathrm{in}}\right\rangle\left(v_{n}\right)_{i}x(t-k)$,
where $\mathbf{v}_{n},\mathbf{v}_{n}^{-1}\in\mathbf{V}$ are, respectively, the
left and right eigenvectors of $\mathbf{W}$, $\mathbf{w}^{\text{in}}$ are the
input weights, and $\lambda_{n}^{k}\in\mathbf{\lambda}$ belongs to the
diagonal matrix containing the eigenvalues of $\mathbf{W}$;
$\mathbf{a}_{i}=\left[a_{i,0},a_{i,1},\ldots\right]$ represents the
coefficients that the previous inputs $\mathbf{x}_{t}=[x(t),x(t-1),\ldots]$
have on $r_{i}(t)$.

Short Proof: As discussed by [13], the state of the neuron can be represented
as $\mathbf{r}(t)=\mathbf{W}\mathbf{r}(t-1)+\mathbf{w}^{\text{in }}x(t)$,
where $\mathbf{w}^{\text{in }}$ are the input weights. We can simplify this
using the coefficients of the previous inputs and plug this term into the
covariance between two neurons. Hence, writing the input coefficients
$\mathbf{a}$ as a function of the eigenvalues of $\mathbf{W}$,

$\displaystyle\mathbf{r}(t)$
$\displaystyle=\sum_{k=0}^{\infty}\mathbf{W}^{k}\mathbf{w}^{\mathrm{in}}x(t-k)=\sum_{k=0}^{\infty}\left(\mathbf{V}\mathbf{\Lambda}^{k}\mathbf{V}^{-1}\right)\mathbf{w}^{\mathrm{in}}x(t-k)$
$\displaystyle\Rightarrow
r_{i}(t)=\sum_{k=0}^{N_{R}}\sum_{n=1}^{N_{R}}\lambda_{n}^{k}\left\langle
v_{n}^{-1},\mathbf{w}^{\mathrm{in}}\right\rangle\left(v_{n}\right)_{i}x(t-k)\blacksquare$

Theorem 1:  If the memory capacity of the HRSNN and MRSNN networks are denoted
by $\mathcal{C}_{H}$ and $\mathcal{C}_{M}$ respectively, then,
$\displaystyle\mathcal{C}_{H}\geq\mathcal{C}_{M}$, where the heterogeneity in
the neuronal parameters $\mathcal{H}$ varies inversely to the correlation
among the neuronal states measured as
$\sum_{n=1}^{N_{\mathcal{R}}}\sum_{m=1}^{N_{\mathcal{R}}}\operatorname{Cov}^{2}\left(x_{n}(t),x_{m}(t)\right)$
which in turn varies inversely with $\mathcal{C}$.

Intuitive Proof:  [13] showed that the $\mathcal{C}$ increases when the
variance along the projections of the input into the recurrent layer are
uniformly distributed. We show that this can be achieved efficiently by using
heterogeneity in the LIF dynamics. More formally, let us express the
projection in terms of the state space of the recurrent layer. We show that
the raw variance in the neuronal states $\mathcal{J}$ can be written as
$\displaystyle\mathcal{J}=\frac{\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}^{2}(\mathbf{\mathbf{\Sigma}})}{\left(\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}(\mathbf{\mathbf{\Sigma}})\right)^{2}}$
where $\lambda_{n}(\mathbf{\Sigma})$ is the $n$th eigenvalue of
$\mathbf{\Sigma}$. We further show that with higher $\mathcal{H}$, the
magnitude of the eigenvalues of $\mathbf{W}$ decreases and hence leads to a
higher $\mathcal{J}$. Now, we project the inputs into orthogonal directions of
the network state space and model the system as
$\displaystyle\mathbf{r}(t)=\sum_{\tau=1}^{\infty}\mathbf{a}_{\tau}x(t-\tau)+\varepsilon_{r}(t)$
where the vectors $\mathbf{a}_{\tau}\in\mathbb{R}^{N}$ are correspond to the
linearly extractable effect of $x(t-\tau)$ onto $\mathbf{r}(t)$ and
$\varepsilon_{r}(t)$ is the nonlinear contribution of all the inputs onto the
state of $\mathbf{r}(t)$. First, we show that $\mathcal{C}$ increases when the
variance along the projections of the input into the recurrent layer is more
uniform. Intuitively, the variances at directions $\mathbf{a}_{\tau}$ must fit
into the variances of the state space, and since the projections are
orthogonal, the variances must be along orthogonal directions. Hence, we show
that increasing the correlation among the neuronal states increases the
variance of the eigenvalues, which would decrease our memory bound
$\mathcal{C}^{*}$. We show that heterogeneity is inversely proportional to
$\displaystyle\sum_{n=1}^{N_{\mathcal{R}}}\operatorname{Cov}^{2}\left(x_{n}(t),x_{m}(t)\right)$.
We see that increasing the correlations between neuronal states decreases the
heterogeneity of the eigenvalues, which reduces $\mathcal{C}$. We show that
the variance in the neuronal states is bounded by the determinant of the
covariance between the states; hence, covariance increases when the neurons
become correlated. As $\mathcal{H}$ increases, neuronal correlation decreases.
[13] proved that the neuronal state correlation is inversely related to
$\mathcal{C}$. Hence, for HRSNN, with $\mathcal{H}>0$,
$\displaystyle\mathcal{C}_{H}\geq\mathcal{C}_{M}$. $\blacksquare$

Spiking Efficiency We analytically prove that the average firing rate of HRSNN
is lesser than the average firing rate of the MRSNN model by considering a
subnetwork of the HRSNN network and modeling the pre-and post-synaptic spike
trains using a nonlinear interactive Hawkes process with inhibition, as
discussed by [53].

Lemma 3.2.1:  If the neuronal firing rate of the HRSNN network with only
heterogeneity in LTP/LTD dynamics of STDP is represented as $\Phi_{R}$ and
that of MRSNN represented as $\Phi_{M}$, then the HRSNN model promotes
sparsity in the neural firing which can be represented as
$\displaystyle\Phi_{R}<\Phi_{M}$.

Short Proof:  We show that the average firing rate of the model with
heterogeneous STDP (LTP/LTD) dynamics (averaged over the population of
neurons) is lesser than the corresponding average neuronal activation rate for
a model with homogeneous STDP dynamics. We prove this by taking a sub-network
of the HRSNN model. Now, we model the input spike trains of the pre-synaptic
neurons using a multivariate interactive, nonlinear Hawkes process with
multiplicative inhibition. Let us consider a population of neurons of size $N$
that is divided into population $A$ (excitatory) and population $B$
(inhibitory). We use a particular instance of the model given in terms of a
family of counting processes $\left(Z_{t}^{1},\ldots,Z_{t}^{N_{A}}\right.$ )
(population $\left.A\right)$ and
$\left(Z_{t}^{N_{A}+1},\ldots,Z_{t}^{N}\right)$ (population $B$ ) with coupled
conditional stochastic intensities given respectively by $\lambda^{A}$ and
$\lambda^{B}$ as follows:

$\displaystyle\lambda_{t}^{A,N}:=$
$\displaystyle\Phi_{A}\left(\frac{1}{N}\sum_{j\in
A}\int_{0}^{t^{-}}h_{1}(t-u)\mathrm{d}Z_{u}^{j}\right)+$
$\displaystyle\Phi_{B\rightarrow A}\left(\frac{1}{N}\sum_{j\in
B}\int_{0}^{t^{-}}h_{2}(t-u)\mathrm{d}Z_{u}^{j}\right)$
$\displaystyle\lambda_{t}^{B,N}:=$
$\displaystyle\Phi_{B}\left(\frac{1}{N}\sum_{j\in
B}\int_{0}^{t^{-}}h_{3}(t-u)\mathrm{d}Z_{u}^{j}\right)+$
$\displaystyle\Phi_{A\rightarrow B}\left(\frac{1}{N}\sum_{j\in
A}\int_{0}^{t^{-}}h_{4}(t-u)\mathrm{d}Z_{u}^{j}\right)$ (11)

where $A,B$ are the populations of the excitatory and inhibitory neurons,
respectively, $\lambda_{t}^{i}$ is the intensity of neuron $i,\Phi_{i}$ a
positive function denoting the firing rate, and $h_{j\rightarrow i}(t)$ is the
synaptic kernel associated with the synapse between neurons $j$ and $i$.
Hence, we show that the heterogeneous STDP dynamics increase the synaptic
noise due to the heavy tail behavior of the system. This increased synaptic
noise leads to a reduction in the number of spikes of the post-synaptic
neuron. Intuitively, a heterogeneous STDP leads to a non-uniform scaling of
correlated spike-trains leading to de-correlation. Hence, we can say that
heterogeneous STDP models have learned a better-orthogonalized subspace
representation, leading to better encoding of the input space with fewer
spikes. $\blacksquare$

Theorem 2: For a given number of neurons $N_{\mathcal{R}}$, the spike
efficiency of the model
$\mathcal{E}=\frac{\mathcal{C}(N_{\mathcal{R}})}{\tilde{S}}$ for HRSNN
($\mathcal{E}_{R}$) is greater than MRSNN ($\mathcal{E}_{M}$) i.e.,
$\mathcal{E}_{R}\geq\mathcal{E}_{M}$

Short Proof:  First, using Lemma 3.2.1, we show that the number of spikes
decreases when we use heterogeneity in the LTP/LTD Dynamics. Hence, we compare
the efficiencies of HRSNN with that of MRSNN as follows:

$\displaystyle\frac{\mathcal{E}_{R}}{\mathcal{E}_{M}}$
$\displaystyle=\frac{\mathcal{C}_{R}(N_{\mathcal{R}})\times\tilde{S}_{M}}{\tilde{S}_{R}\times\mathcal{C}_{M}(N_{\mathcal{R}})}$
$\displaystyle=\frac{\sum_{\tau=1}^{N_{\mathcal{R}}}\frac{\operatorname{Cov}^{2}\left(x(t-\tau),\mathbf{a}^{R}_{\tau}\mathbf{r}_{R}(t)\right)}{\operatorname{Var}\left(\mathbf{a}^{R}_{\tau}\mathbf{r}_{R}(t)\right)}\times\int\limits_{t_{ref}}^{\infty}t\Phi_{R}dt}{\sum_{\tau=1}^{N_{\mathcal{R}}}\frac{\operatorname{Cov}^{2}\left(x(t-\tau),\mathbf{a}^{M}_{\tau}\mathbf{r}_{M}(t)\right)}{\operatorname{Var}\left(\mathbf{a}^{M}_{\tau}\mathbf{r}_{M}(t)\right)}\times\int\limits_{t_{ref}}^{\infty}t\Phi_{M}dt}$
(12)

Since $S_{R}\leq S_{M}$ and also,the covariance increases when the neurons
become correlated, and as neuronal correlation decreases, $\mathcal{H}_{X}$
increases (Theorem 1), we see that $\mathcal{E}_{R}/\mathcal{E}_{M}\geq
1\Rightarrow\mathcal{E}_{R}\geq\mathcal{E}_{M}\quad\blacksquare$

Optimal Heterogeneity using Bayesian Optimization for Distributions

To get optimal heterogeneity in the neuron and STDP dynamics, we use a
modified Bayesian Optimization (BO) technique. However, using BO for high-
dimensional problems remains a significant challenge. In our case, optimizing
HRSNN model parameters for 5000 neurons requires the optimization of two
parameters per neuron and four parameters per STDP synapse, where standard BO
fails to converge to an optimal solution. However, the parameters to be
optimized are correlated and can be drawn from a probability distribution as
shown by [8]. Thus, we design a modified BO to estimate parameter
distributions instead of individual parameters for the LIF neurons and the
STDP synapses, for which we modify the BO’s surrogate model and acquisition
function. This makes our modified BO highly scalable over all the variables
(dimensions) used. The loss for the surrogate model’s update is calculated
using the Wasserstein distance between the parameter distributions. We use the
modified Matern function on the Wasserstein metric space as a kernel function
for the BO problem. BO uses a Gaussian process to model the distribution of an
objective function and an acquisition function to decide points to evaluate.
For data points $x\in X$ and the corresponding output $y\in Y$, an SNN with
network structure $\mathcal{V}$ and neuron parameters $\mathcal{W}$ acts as a
function $f_{\mathcal{V},\mathcal{W}}(x)$ that maps input data $x$ to $y$. The
optimization problem can be defined as:
$\min_{\mathcal{V},\mathcal{W}}\sum_{x\in X,y\in
Y}\mathcal{L}\left(y,f_{\mathcal{V},\mathcal{W}}(x)\right)$ where
$\mathcal{V}$ is the set of hyperparameters of the neurons in $\mathcal{R}$
and $\mathcal{W}$ is the multi-variate distribution constituting the
distributions of: (i) the membrane time constants $\tau_{m-E},\tau_{m-I}$ of
LIF neurons, (ii) the scaling function constants $(A_{+},A_{-})$ and (iii) the
decay time constants $\tau_{+},\tau_{-}$ for the STDP learning rule in
$\mathcal{S}_{\mathcal{RR}}$.

##  V Lyapunov Noise Pruning (LNP) Method

Algorithm 1 Lyapunov Noise Pruning (LNP) Method

1: Step 1: Synapse Pruning using Spectral Graph Pruning

2: for each $e_{ij}$ in $A$ connecting $n_{i},n_{j}$ do

3: Find $\mathcal{N}(e_{ij})$

4: Compute $\lambda_{k}$ for $k$ in $\mathcal{N}(e_{ij})$

5: Define $W$ using $\text{harmonic mean}(\lambda_{k})$ for $k$ in
$\mathcal{N}(e_{ij})$ (Eq 13)

6: Use $\boldsymbol{b}(t)$ at each node

7: Compute $C$ of firing rates

8: Preserve each $e_{ij}$ with $p_{ij}$ yielding $A^{\text{sparse}}$ for
$i\neq j$

9: end for

10: Step 2: Node Pruning using Betweenness Centrality

11: for each $n_{i}$ in Network do

12: Compute $B(n_{i})$

13: if $B(n_{i})<\text{threshold}$ then

14: Remove $n_{i}$

15: end if

16: end for

17: Step 3: Delocalization of the Eigenvectors

18: for each $A^{\text{pruned}}$ do

19: Add edges to preserve eigenvectors and maintain stability

20: end for

21: Step 4: Neuronal Timescale Optimization

22: for each pruned model $m$ do

23: Use Lyapunov spectrum $L$ to optimize neuronal timescales $\tau_{i}\forall
i\in\mathcal{R}$ using BO.

24: end for

25: return $A^{\text{pruned}}$

The proposed research presents a pruning algorithm employing spectral graph
pruning and Lyapunov exponents in an unsupervised model. We calculate the
Lyapunov matrix, optimizing for ratios and rates to handle extreme data values
and incorporate all observations. After pruning, nodes with the lowest
betweenness centrality are removed to improve network efficiency, and select
new edges are added during the delocalization phase to maintain stability and
integrity. This method balances structural integrity with computational
efficacy, contributing to advancements in network optimization. Algorithm 1
shows a high-level algorithm for the entire process.

Figure 3: (a) Concept of HRSNN with variable Neuronal Dynamics (b) Figure
showing the task-agnostic pruning and training of the CHRSNN/HRSNN networks
using LNP in comparison to the current approach

Step I: Noise-Pruning of Synapses:  First, we define the Lyapunov matrix of
the network. To formalize the concept of the Lyapunov Matrix, let us consider
a network represented by a graph $G(V,E)$ where $V$ is the set of nodes and
$E$ is the set of edges. For each edge $e_{ij}$ connecting nodes $z$, let
$N(z)$ be the set of neighbors of nodes $z,z=\\{i,j\\}$. Thus, the Lyapunov
exponents corresponding to these neighbors are represented as $\Lambda(N(z))$.
The element $L_{ij}$ of the Lyapunov Matrix ($L$) is then calculated using the
harmonic mean of the Lyapunov exponents of the neighbors of nodes $i,j$ as:

$L_{ij}=\frac{n\cdot|\Lambda(N(i))\cup\Lambda(N(j))|}{\sum_{\lambda\in\Lambda(N(i))\cup\Lambda(N(j))}\frac{1}{\lambda}}$
(13)

where $\lambda$ denotes individual Lyapunov exponents from the set of all such
exponents of nodes $i$ and $j$’s neighbors. $L$, encapsulates the impact of
neighboring nodes on each edge regarding their Lyapunov exponents, which helps
us evaluate the network’s stability and dynamical behavior. Through the
harmonic mean, the matrix accommodates the influence of all neighbors,
including those with extreme Lyapunov exponents, for a balanced depiction of
local dynamics around each edge. Thus, the linearized network around
criticality is represented as:

$\dot{\boldsymbol{x}}=-D\boldsymbol{x}+L\boldsymbol{x}+\boldsymbol{b}(t)=A\boldsymbol{x}+\boldsymbol{b}(t)$
(14)

Here $\boldsymbol{x}$ represents the firing rate of $N$ neurons, with $x_{i}$
specifying the firing rate of neuron $i$. $\boldsymbol{b}(t)$ denotes external
input, including biases, and $L$ is the previously defined Lyapunov matrix
between neurons. $D$ is a diagonal matrix indicating neurons’ intrinsic leak
or excitability, and $A$ is defined as $A=-D+L$. The intrinsic
leak/excitability, $D_{ii}$, quantifies how the firing rate, $x_{i}$, of
neuron $i$ alters without external input or interaction, impacting the neural
network’s overall dynamics along with $\boldsymbol{x}$, $L$, and external
inputs $\boldsymbol{b}(t)$. Positive $D_{ii}$ suggests increased excitability
and firing rate, while negative values indicate reduced neuron activity over
time. We aim to create a sparse network ($A^{\text{sparse}}$) with fewer edges
while maintaining dynamics similar to the original network. The sparse network
is thus represented as:

$\displaystyle\dot{\boldsymbol{x}}$
$\displaystyle=A^{\text{sparse}}\boldsymbol{x}+\boldsymbol{b}(t)$ s.t.
$\displaystyle\quad\left|\boldsymbol{x}^{T}\left(A^{\text{sparse}}-A\right)\boldsymbol{x}\right|\leq\epsilon\left|\boldsymbol{x}^{T}A\boldsymbol{x}\right|\quad\forall\boldsymbol{x}\in\mathbb{R}^{N}$
(15)

for some small $\epsilon>0$. When the network in Eq. 14 is driven by
independent noise at each node, we define
$\boldsymbol{b}(t)=\boldsymbol{b}+\sigma\boldsymbol{\xi}(t)$, where
$\boldsymbol{b}$ is a constant input vector, $\boldsymbol{\xi}$ is a vector of
IID Gaussian white noise, and $\sigma$ is the noise standard deviation. Let
$\Sigma$ be the covariance matrix of the firing rates in response to this
input. The probability $p_{ij}$ for the synapse from neuron $j$ to neuron $i$
with the Lyapunov exponent $l_{ij}$ is defined as:

$p_{ij}=\begin{cases}\rho l_{ij}(\Sigma{ii}+\Sigma{jj}-2\Sigma{ij})&\text{for
}w_{ij}>0\text{ (excitatory)}\\\
\rho|l_{ij}|(\Sigma{ii}+\Sigma{jj}+2\Sigma{ij})&\text{for }w_{ij}<0\text{
(inhibitory)}\end{cases}$ (16)

Here, $\rho$ determines the density of the pruned network. The pruning process
independently preserves each edge with probability $p_{ij}$, yielding
$A^{\text{sparse}}$, where $\displaystyle
A_{ij}^{\text{sparse}}=A_{ij}/p_{ij}$, with probability $p_{ij}$ and $0$
otherwise. For the diagonal elements, denoted as $A_{ii}^{\text{sparse}}$,
representing leak/excitability, we either retain the original diagonal,
setting $A_{ii}^{\text{sparse}}=A_{ii}$, or we introduce a perturbation,
$\Delta_{i}$, defined as the difference in total input to neuron $i$, and
adjust the diagonal as $A_{ii}^{\text{sparse}}=A_{ii}-\Delta_{i}$.
Specifically, $\Delta_{i}=\sum_{j\neq i}|A_{ij}^{\text{sparse}}|-\sum_{j\neq
i}|A_{ij}|$. This perturbation, $\Delta_{i}$, is typically minimal with a zero
mean and is interpreted biologically as a modification in the excitability of
neuron $i$ due to alterations in total input, aligning with the known
homeostatic regulation of excitability.

Step II: Node Pruning based on Betweenness Centrality: In addressing network
optimization, we propose an algorithm specifically designed to prune nodes
with the lowest betweenness centrality ($C_{b}$) in a given graph, thereby
refining the graph to its most influential components. $C_{b}$ quantifies the
influence a node has on information flow within the network. Thus, we
calculate $C_{b}$ for each node and prune the nodes with the least values
below a given threshold, ensuring the retention of nodes integral to the
network’s structural and functional integrity.

Step III: Delocalizing Eigenvectors: To preserve eigenvector delocalization
and enhance long-term prediction performance post-pruning, we introduce a
predetermined number of additional edges to counteract eigenvalue localization
due to network disconnection. Let $G=(V,E)$ and $G^{\prime}=(V,E^{\prime})$
represent the original and pruned graphs, respectively, where
$E^{\prime}\subset E$. We introduce additional edges, $E^{\prime\prime}$, to
maximize degree heterogeneity, $H$, defined as the variance of the degree
distribution $\displaystyle H=\frac{1}{|V|}\sum_{v\in V}(d(v)-\bar{d})^{2}$,
subject to $|E^{\prime}\cup E^{\prime\prime}|\leq L$, where $L$ is a
predetermined limit. This is formalized as an optimization problem to find the
optimal set of additional edges, $E^{\prime\prime}$, enhancing eigenvector
delocalization and improving the pruned network’s predictive performance.
Given graph $G=(V,E)$, with adjacency matrix $A$ and Laplacian matrix $L$, we
analyze the eigenvalues and eigenvectors of $L$ to study eigenvector
localization. To counteract localization due to pruning, we introduce a fixed
number, $m$, of additional edges to maximize the variance of the degree
distribution, $\text{Var}(D)$, within local neighborhoods, formalized as:

$\displaystyle\underset{E^{\prime\prime}}{\max}$
$\displaystyle\quad\text{Var}(D)$ s.t.
$\displaystyle\quad|E^{\prime\prime}|=m,\quad E^{\prime}\cap
E^{\prime\prime}=\emptyset,\quad E^{\prime\prime}\subseteq\bigcup_{v_{i}\in
V}N(v_{i})\times\\{v_{i}\\}$ (17)

This approach ensures eigenvector delocalization while preserving structural
integrity and specified sparsity, optimizing the model’s long-term predictive
performance.

The goal is to maximize the variance of the degree distribution,
$\text{Var}(D)$, by selecting the best set of additional edges
$E^{\prime\prime}$ such that: (i) The number of additional edges is $m$. (ii)
The additional edges are not part of the original edge set $E^{\prime}$. (iii)
The additional edges are selected from the neighborhoods of the vertices.

Step IV: Neuronal Timescale Optimization In optimizing RSNN, known for their
complex dynamical nature, we employ the Lyapunov spectrum to refine neuronal
timescales using Bayesian optimization. During optimization, RSNNs are
inherently unstable due to their variable parameters and pruning processes,
affecting training dynamics’ stability and model learnability. Lyapunov
exponents, as outlined in [54], are crucial for understanding system stability
and are linked to the operational efficiency of RNNs. This paper utilizes the
Lyapunov spectrum as a criterion for optimizing neuronal timescales in pruned
networks, aiming to maintain stability and functionality while minimizing
instability risks inherent in pruning.

##  VI Results

###  VI-A Ablation Studies

Heterogeneity Parameter Importance:  We use SAGE (Shapley Additive Global
importancE) [55], a game-theoretic approach to understand black-box models to
calculate the significance of adding heterogeneity to each parameter for
improving $\mathcal{C}$ and $\tilde{S}$. SAGE summarizes the importance of
each feature based on the predictive power it contributes and considers
complex feature interactions using the principles of Shapley value, with a
higher SAGE value signifying a more important feature. We tested the HRSNN
model using SAGE on the Lorenz96 and the SHD datasets. The results are shown
in Fig. 4. We see that $\tau_{m}$ has the greatest SAGE values for
$\mathcal{C}$, signifying that it has the greatest impact on improving
$\mathcal{C}$ when heterogeneity is added. Conversely, we see that
heterogeneous STDP parameters (viz., $\tau_{\pm},\eta_{\pm}$) play a more
critical role in determining the average neuronal spike activation. Hence, we
confirm the notions proved in Sec. 3 that heterogeneity in neuronal dynamics
improves the $\mathcal{C}$ while heterogeneity in STDP dynamics improves the
spike count. Thus, we need to optimize the heterogeneity of both to achieve
maximum $\mathcal{E}$.

Figure 4: Bar chart showing the global importance of different heterogeneous
parameters using HRSNN on the dataset. The experiments were repeated five
times with different parameters from the same distribution (a) Classification
(b) Prediction

We compare the performance of the HRSNN model with heterogeneity in the LIF
and STDP dynamics (HeNHeS) to the ablation baseline recurrent spiking neural
network models described above. We run five iterations for all the baseline
cases and show the mean and standard deviation of the prediction accuracy of
the network using 2000 neurons. The results are shown in Table I. We see that
the heterogeneity in the LIF neurons and the LTP/LTD dynamics significantly
improve the model’s accuracy and error.

TABLE I: Table comparing the performance of RSNN with Homogeneous and
Heterogeneous LIF neurons using different learning methods with 2000 neurons

Datasets | KTH | DVS128
---|---|---
| Neuron
---
Type
| Homogeneous
---
STDP
| Heterogeneous
---
STDP
BackPropagation |  | Homogeneous
---
STDP
| Heterogeneous
---
STDP
BackPropagation
| Homogeneous
---
LIF
$86.33\pm 4.05$ | $91.37\pm 3.15$ | $94.87\pm 2.03$ | 90.33 $\pm$ 3.41 | 93.37 $\pm$ 3.05 | 97.06 $\pm$ 2.29
| Heterogeneous
---
LIF
92.16 $\pm$ 3.17 | 94.32 $\pm$ 1.71 | 96.84 $\pm$ 1.96 | 92.16 $\pm$ 2.97 | 96.54 $\pm$ 1.82 | 98.12 $\pm$ 1.97

Figure 5: Comparison of performance of HRSNN models for the (a) KTH dataset
and (b) DVS128 dataset for varying number of neurons. The bar graph (left
Y-axis) shows the difference between the accuracies between HeNHeS and HoNHoS
models. The line graphs (right Y-axis) shows the Accuracies (%) for the four
ablation networks (HoNHoS, HeNHoS, HoNHeS, HeNHeS)

###  VI-B Number of neurons

In deep learning, it is an important task to design models with a lesser
number of neurons without undergoing degradation in performance. We
empirically show that heterogeneity plays a critical role in designing spiking
neuron models of smaller sizes. We compare models’ performance and convergence
rates with fewer neurons in $\mathcal{R}$.

Performance Analysis:  We analyze the network performance and error when the
number of neurons is decreased from 2000 to just 100. We report the results
obtained using the HoNHoS and HeNHeS models for the KTH and DVS-Gesture
datasets. The experiments are repeated five times, and the observed mean and
standard deviation of the accuracies are shown in Figs. 5. The graphs show
that as the number of neurons decreases, the difference in accuracy scores
between the homogeneous and the heterogeneous networks increases rapidly.

Convergence Analysis with lesser neurons: Since the complexity of BO increases
exponentially on increasing the search space, optimizing the HRSNN becomes
increasingly difficult as the number of neurons increases. Thus, we compare
the convergence behavior of the HoNHoS and HeNHeS models with 100 and 2000
neurons each. The results are plotted in Fig. 6(a), (b). Despite the huge
number of additional parameters, the convergence behavior of HeNHeS is similar
to that of HoNHoS. Also, it must be noted that once converged, the standard
deviation of the accuracies for HeNHeS is much lesser than that of HoNHoS,
indicating a much more stable model.

Figure 6: Plots showing comparison of the convergence of the BO with
increasing functional evaluations for the (a) KTH and (b)DVSGesture dataset
for varying number of neurons

###  VI-C Sparse Connections

Figure 7: Change in the number of active neurons with network sparsity and
weight variance, for (a)HoNHoS and (b)HeNHeS. The plot is obtained by
interpolating 81 points, and each point is calculated by averaging the results
from 5 randomly initialized HRSNNs.

$\mathcal{S}_{\mathcal{RR}}$ is generated using a probability dependent on the
Euclidean distance between the two neurons, as described by Eq. 2, where
$\lambda$ controls the density of the connection, and $C$ is a constant
depending on the type of the synapses [56].

We performed various simulations using a range of values for the connection
parameter $\lambda$ and synaptic weight scale $W_{\text{scale}}$. Increasing
$\lambda$ will increase the number of synapses. Second, the $W_{\text{scale}}$
parameter determines the mean synaptic strength. Now, a greater
$W_{\text{scale}}$ produces larger weight variance. For a single input video,
the number of active neurons was calculated and plotted against the parameter
values for synaptic weight $W_{\text{scale}}$ and network connectivity
$\lambda$. Active neurons are those that fire at least one spike over the
entire test data set. The results for the HoNHoS and HeNHeS are shown in Figs.
7a, 7b respectively. Each plot in the figure is obtained by interpolating 81
points, and each point is calculated by averaging the results from five
randomly initialized with the parameters specified by the point. The
horizontal axis showing the increase in $\lambda$ is plotted on a linear
scale, while the vertical axis showing the variation in $W_{\text{scale}}$ is
in a log scale. The figure shows the neurons that have responded to the inputs
and reflect the network’s activity level. $W_{\text{scale}}$ is a factor that
enhances the signal transmission within $\mathcal{R}$. As discussed by Markram
et al. [57], the synaptic response that is generated by any action potential
(AP) in a train is given as $EPSP_{n}=W_{\text{scale}}\times\rho_{n}\times
u_{n}$, where $\rho_{n}$ is the fraction of the synaptic efficacy for the
$n$-th AP, and $u_{n}$ is its utilization of synaptic efficacy. Hence, it is
expected that when the $W_{\text{scale}}$ is large, more neurons will fire. As
$\lambda$ increases, more synaptic connections are created, which opens up
more communication channels between the different neurons. As the number of
synapses increases, the rank of the final state matrix used to calculate
separation property also increases. The rank reaches an optimum for
intermediate synapse density, and the number of synapses created increases
steadily as $\lambda$ increases. As $\lambda$ increases, a larger number of
connections creates more dependencies between neurons and decreases the
effective separation ranks when the number of connections becomes too large.
The results for the variation of the effective ranks with $\lambda$ and
$W_{\text{scale}}$ are shown in the Supplementary.

We compare the model’s change in performance with varying sparsity levels in
the connections and plotted in Figs. 7a,b for the HoNHoS and the HeNHeS
models. From the figures, we see that for larger values of $\lambda$, the
performance of both the RSNNs was suboptimal and could not be improved by
tuning the parameter $W_{\text{scale}}$. For a small number of synapses, a
larger $W_{\text{scale}}$ was required to obtain satisfactory performance for
HoNHoS compared to the HeNHeS model. Hence, the large variance of the weights
leads to better performance. Hence, we see that the best testing accuracy for
HeNHeS is achieved with fewer synapses than HoNHoS. It also explains why the
highest testing accuracy for the heterogeneous network (Fig. 7b.) is better
than the homogeneous network (Fig. LABEL:fig:3a), because the red region in
Fig. 7b corresponds to higher $W_{\text{scale}}$ values and thus larger weight
variance than Fig. 7a.

###  VI-D Limited Training Data

In this section, we compare the performance of the HeNHeS to HoNHoS and HeNB-
based spiking recurrent neural networks that are trained with limited training
data. The evaluations performed on the KTH dataset are shown in Fig. 8 as a
stacked bar graph for the differential increments of training data sizes. The
figure shows that using $10\%$ training data, HeNHeS models outperform both
HoNHoS and HeNB for all the cases. The difference between the HeNHeS and HeNB
increases as the number of neurons in the recurrent layer $N_{\mathcal{R}}$
decreases. Also, we see that adding heterogeneity improves the model’s
performance in homogeneous cases. Even when using 2000 neurons, HeNHeS trained
with $10\%$ training data exhibit similar performance to HeNB trained with
$25\%$ of training data. It is to be noted here that for the performance
evaluations of the cases with $10\%$ training data, the same training was
repeated until each model converged.

Figure 8:  Bar graph showing difference in performance for the different
models with increasing training data for the KTH dataset. A similar trend can
be observed for the DVS dataset (shown in Supplementary).

###  VI-E Comparison with Prior Work

In this section, we compare our proposed HRSNN model with other baseline
architectures. We divide this comparison in two parts as discussed below:

  * •

DNN-based Models: We compare the performance and the model complexities of
current state-of-the-art DNN-based models [58, 59, 60, 61, 62] with our
proposed HRSNN models.

  * •

Backpropagation-based SNN Models: We compare the performance of
backpropagation-based SNN models with HoNB and HeNB-based RSNN models. We
observe that backpropagated HRSNN models (HeNB) can achieve similar
performance to DNN models but with much lesser model complexity (measured
using the number of parameters).

    1. 1.

State-of-the-art BP Homogeneous SNN: We compare the performance of current
state-of-the-art backpropagation-based SNN models [63, 64, 65, 66].

    2. 2.

State-of-the-art BP Heterogeneous SNN: We compare the performances of the
current state-of-the-art SNN models, which uses neuronal heterogeneity [8, 67,
68]. We compare the performances and the model complexities of these models.

    3. 3.

Proposed Heterogeneous Backpropagation Models: We introduce two new
backpropagation-based RSNN models. These models are the Homogeneous Neurons
with Backpropagation (HoNB) and the Heterogeneous Neurons with Backpropagation
(HeNB). We use our novel Bayesian Optimization to search for the parameters
for both of these models.

  * •

Unsupervised SNN Models: We also compare the results for some state-of-the-art
unsupervised SNN models with our proposed HRSNN models.

    1. 1.

Homogeneous SNN Models:  We compare the performances of some of the state-of-
the-art unsupervised SNN models which uses homogeneous neuronal parameters
[69, 70, 56].

    2. 2.

HRSNN Models:  We compare the above models with respect to our proposed HRSNN
models using heterogeneity in both neuronal and synaptic parameters. We
compare the model’s performance and the model’s complexity.

We also compare the average neuronal activation of the homogeneous and the
heterogeneous recurrent SNNs for the same given image input for a recurrent
spiking network with 2000 neurons. If we consider the neuronal activation of
neuron $i$ at time $t$ to be $\nu_{i}(t)$, the average neuronal activation
$\bar{\nu}$ for $T$ timesteps is defined as
$\displaystyle\bar{\nu}=\frac{\Sigma_{i=0}^{N_{\mathcal{R}-1}}\Sigma_{t=0}^{T}{\nu_{i}(t)}}{N_{\mathcal{R}}}$.

The results obtained are shown in Table II. The Table shows that the
heterogeneous HRSNN model has a much lesser average neuronal activation than
the homogeneous RSNN and the other unsupervised SNN models. Thus, we conclude
that HeNHeS induces sparse activation and sparse coding of information. Again,
comparing state-of-the-art unsupervised learning models for action recognition
with our proposed HRSNN models, we see that using heterogeneity in the
unsupervised learning models can substantially improve the model’s performance
while having much lesser model complexity.

TABLE II: Table showing the comparison of the Performance and the Model
Complexities for DNN and Supervised and Unsupervised SNN models

| Model | MACs/ACs | RGB Datasets | Event Dataset
---|---|---|---|---
KTH | UCF11 | UCF101 |  | DVS
---
Gesture-128
Supervised Learning Method
DNN |  | PointNet
---
[58]
MAC: $152\times 10^{9}$ | - | - | - | 95.3
| RG-CNN
---
[59]
MAC: $53\times 10^{9}$ | - | - | - | 97.2
| I3D
---
[60]
MAC: $188\times 10^{9}$ | - | 90.9 | - | 96.5
| 3D-ResNet-34
---
[61]
MAC: $78.43\times 10^{9}$ | 94.78 | 83.72 | - | -
| 3D-ResNet-50
---
[61]
MAC: $62.09\times 10^{9}$ | 92.31 | 81.44 | - | -
| TDN
---
[62]
MAC: $69.67\times 10^{9}$ | 99.15 | 98.03 | 97.4 | -
SNN- Supervised (Homogeneous) |  | STBP-tdBN
---
[63]
AC: $15.13\times 10^{7}$ | - | - | - | 96.87
| Shen et al.
---
[64]
AC: $12.14\times 10^{7}$ | - | - | - | 98.26
| Liu et al.
---
[65]
AC: $27.59\times 10^{7}$ | 90.16 | - | - | 92.7
| Panda et al.
---
[66]
AC: $40.4\times 10^{7}$ | - | - | 81.3 | -
| HoNB
---
(2000 Neurons)
AC: $9.54\times 10^{7}$ | 94.87 | 82.89 | 80.25 | 97.06
SNN- Supervised (Heterogeneous) |  | Perez et al.
---
[8]
AC: $8.94\times 10^{7}$ | - | - | - | 82.9
| Fang et al.
---
[68]
AC: $15.32\times 10^{7}$ | - | - | - | 97.22
| She et al. (BPTT)
---
[67]
AC: $13.25\times 10^{7}$ | - | - | - | 98.0
| HeNB
---
(2000 Neurons)
AC: $9.18\times 10^{7}$ | 96.84 | 88.36 | 84.32 | 98.12
| Model |  Number of Neurons |  MACs/ACs/ Avg. Neuron Activation ($\bar{\nu}$) | RGB Datasets | Event Daset
KTH | UCF11 | UCF101 |  | DVS
---
Gesture 128
Unsupervised Learning Method
DNN - Unsupervised |  | MetaUVFS
---
[71]
- | MAC: $58.39\times 10^{9}$ | 90.14 | 80.79 | 76.38 | -
| Soomro et al.
---
[72]
- | MAC: $63\times 10^{9}$ | 84.49 | 73.38 | 61.2 | -
SNN- Unsupervised (Homogeneous) |  | GRN-BCM
---
[69]
1536 | $\bar{\nu}=3.56\times 10^{3}$ | 74.4 | - | - | 77.19
| LSM STDP
---
[70]
135 | $\bar{\nu}=10.12\times 10^{3}$ | 66.7 | - | - | 67.41
| GP-Assisted CMA-ES
---
[56]
500 | $\bar{\nu}=9.23\times 10^{3}$ | 87.64 | - | - | 89.25
RSNN-STDP Unsupervised (Ours) | HoNHoS | 2000 | $\bar{\nu}=3.85\times 10^{3}$ | 86.33 | 75.23 | 74.45 | 90.33
HeNHeS | 500 | $\bar{\nu}=2.93\times 10^{3}$ | 88.04 | 71.42 | 70.16 | 90.15
HeNHeS | 2000 | $\bar{\nu}=2.74\times 10^{3}$ | 94.32 | 79.58 | 77.33 | 96.54

Figure 9: Figure showing the results of ablation studies of BO for the
following three cases: (a) for prediction problem - the radius indicates the
normalized NRMSE loss with a smaller radius indicating lower (better) NRMSE
(b) for classification problem - the radius indicates the normalized accuracy
with larger radius indicating higher (better) accuracy. The numbers represent
the number of neurons used in each model, and the line joins the three
corresponding models with the same model size.

Datasets:  Classification:  We use the Spoken Heidelberg Digits (SHD) spiking
dataset to benchmark the HRSNN model with other standard spiking neural
networks [73].
Prediction:  We use a multiscale Lorenz 96 system [74] which is a set of
coupled nonlinear ODEs and an extension of Lorenz’s original model for
multiscale chaotic variability of weather and climate systems which we use as
a testbed for the prediction capabilities of the HRSNN model [75].

Bayesian Optimization Ablation Studies: First, we perform an ablation study of
BO for the following three cases: (i) Using Memory Capacity $\mathcal{C}$ as
the objective function (ii) Using Average Spike Count $\tilde{S}$ as the
objective function and (iii) Using $\mathcal{E}$ as the objective function. We
optimize both LIF neuron parameter distribution and STDP dynamics
distributions for each. We plot $\mathcal{C}$, $\tilde{S}$, the empirical
spike efficiency $\hat{\mathcal{E}}$, and the observed RMSE of the model
obtained from BO with different numbers of neurons. The results for
classification and prediction problems are shown in Fig. 9(a) and (b),
respectively. Ideally, we want to design networks with high $\mathcal{C}$ and
low spike count, i.e., models in the upper right corner of the graph. The
observed results show that BO using $\mathcal{E}$ as the objective gives the
best accuracy with the fewest spikes. Thus, we can say that this model has
learned a better-orthogonalized subspace representation, leading to better
encoding of the input space with fewer spikes. Hence, for the remainder of
this paper, we focus on this BO model, keeping the $\mathcal{E}$ as the
objective function. This Bayesian Optimization process to search for the
optimal hyperparameters of the model is performed before training and
inference using the model and is generally equivalent to the network
architecture search process used in deep learning. Once we have these optimal
hyper-parameters, we freeze these hyperparameters, learn (unsupervised) the
network parameters (i.e., synaptic weights) of the HRSNN while using the
frozen hyperparameters, and generate the final HRSNN model for inference. In
other words, the hyperparameters, like the distribution of membrane time
constants or the distribution of synaptic time constants for STDP, are fixed
during the learning and inference.

Results:  We perform an ablation study to evaluate the performance of the
HRSNN model and compare it to standard BP-based spiking models. We study the
performances of both the SHD dataset for classification and the Lorenz system
for prediction. The results are shown in Table III. We compare the Normalized
Root Mean Squared Error (NRMSE) loss (prediction), Accuracy (classification),
Average Spike Count $\tilde{S}$ and the application level empirical spiking
efficiency $\hat{\mathcal{E}}$ calculated as
$\displaystyle\frac{1}{\text{NRMSE}\times\tilde{S}}$ (prediction) and
$\displaystyle\frac{\text{Accuracy}}{\tilde{S}}$ (classification). We perform
the experiments using 5000 neurons in $\mathcal{R}$ on both classification and
prediction datasets. We see that the HRSNN model with heterogeneous LIF and
heterogeneous STDP outperforms other HRSNN and MRSNN models in terms of NRMSE
scores while keeping the $\tilde{S}$ much lower than HRSNN with heterogeneous
LIF and homogeneous STDP. From the experiments, we can conclude that the
heterogeneous LIF neurons have the greatest contribution to improving the
model’s performance. In contrast, heterogeneity in STDP has the most
significant impact on a spike-efficient representation of the data. HRSNN with
heterogeneous LIF and STDP leverages the best of both worlds by achieving the
best RMSE with low spike activations, as seen from Table III. We also compare
the generalizability of the HRSNN vs. MRSNN models, where we empirically show
that the heterogeneity in STDP dynamics helps improve the overall model’s
generalizability. In addition, we discuss how HRSNN reduces the effect of
higher-order correlations, thereby giving rise to a more efficient
representation of the state space.

TABLE III: Table showing the comparison of the Accuracy and NRMSE losses for
the SHD Classification and Lorenz System Prediction tasks, respectively. We
show the average spike rate, calculated as the ratio of the moving average of
the number of spikes in a time interval $T$. For this experiment, we choose
$T=4ms$ and a rolling time span of $2ms$, which is repeated until the first
spike appears in the final layer. Following the works of [76], we show that
the normalized average spike rate is the total number of spikes generated by
all neurons in an RSNN averaged over the time interval $T$. The results marked
with ’*’ denotes we implemented the open-source code for the model and
evaluated the given results.

| Method |  | SHD
---
(Classification)
| Chaotic Lorenz System
---
(Prediction)
| Accuracy
---
$(A)$
| Normalized
---
Avg. Firing
Rate ($\frac{\tilde{S}}{T}$)
| Efficiency
---
($\hat{\mathcal{E}}=\frac{A}{\tilde{S}}$)
NRMSE |  | Normalized
---
Avg. Firing
Rate ($\frac{\tilde{S}}{T}$)
| Efficiency
---
$\hat{\mathcal{E}}=\frac{1}{\text{NRMSE}\times\tilde{S}}$
Unsupervised RSNN |  | MRSNN
---
(Homogeneous LIF,
Homogeneous STDP)
73.58 | -0.508 | $18.44\times 10^{-3}$ | 0.395 | -0.768 | $0.787\times 10^{-3}$
| HRSNN
---
(Heterogeneous LIF,
Homogeneous STDP)
78.87 | 0.277 | $17.19\times 10^{-3}$ | 0.203 | -0.143 | $1.302\times 10^{-3}$
| HRSNN
---
(Homogeneous LIF,
Heterogeneous STDP)
74.03 | -1.292 | $22.47\times 10^{-3}$ | 0.372 | -1.102 | $0.932\times 10^{-3}$
| HRSNN
---
(Heterogeneous LIF,
Heterogeneous STDP)
80.49 | -1.154 | $24.35\times 10^{-3}$ | 0.195 | -1.018 | $1.725\times 10^{-3}$
RSNN with BP |  | MRSNN-BP
---
(Homogeneous LIF,
BP)
81.42 | 0.554 | $16.9\times 10^{-3}$ | 0.182 | 0.857 | $1.16\times 10^{-3}$
| HRSNN-BP
---
(Heterogeneous LIF,
BP)
83.54 | 1.292 | $15.42\times 10^{-3}$ | 0.178 | 1.233 | $1.09\times 10^{-3}$
| Adaptive SRNN
---
[77]
84.46 | $0.831^{*}$ | $17.21^{*}\times 10^{-3}$ | $0.174^{*}$ | $0.941^{*}$ | $1.19^{*}\times 10^{-3}$

Figure 10: Comparative Evaluation of Pruning Methods Across Iterations. Figs.
(a) and (b) show the evolution of the number of synapses and neurons with the
iterations of the LNP and AP algorithms. Fig (c) represents how the RMSE loss
changes when the pruned model after each iteration is trained and tested on
the Lorenz63 dataset

###  VI-F Lyapunov Noise Based Pruning (LNP) Method

The experimental process begins with a randomly initialized HRSNN and CHRSNN.
Pruning algorithms are used to create a sparse network. Each iteration of
pruning results in a sparse network; we experiment with 100 iterations of
pruning. We characterize the neuron and synaptic distributions of the "Sparse
HRSNN" obtained after each pruning iteration to track the reduction of the
complexity of the models with pruning.

We train the sparse HRSNN model obtained after each pruning iteration to
estimate performance on various tasks. Note, the pruning process does not
consider the trained model or its performance during iterations. The sparse
HRSNNs are trained for time-series prediction and image classification tasks.
For the prediction task, the network is trained using 500 timesteps of the
datasets and is subsequently used to predict the following 100 timesteps. For
the classification task, each input image was fed to the input of the network
for Tinput = 100 ms of simulation time in the form of Poisson-distributed
spike trains with firing rates proportional to the intensity of the pixels of
the input images, followed by 100 ms of the empty signal to allow the current
and activity of neurons to decrease. For both tasks, the input data is
converted into spike trains via rate-encoding, forming the high-dimensional
input to the XRSNN. The output spike trains are then processed through a
decoder and readout layer for final predictions or classification results.

Figure 11: Scatter Plot showing Accuracy vs. Avg. SOPs for different pruning
methods on CIFAR10.

Evolution of Complexity of Sparse Models during Pruning: We plot the change in
the number of neurons and synapses for the 100 iteration steps for AP and LNP
Pruning algorithms. The results for the variation of the synapses and neurons
with the iterations of the pruning algorithm are plotted in Figs. 10(a) and
(b), respectively. The LNP methods perform better than the activity-based
pruning method and converge to a model with fewer neurons and synapses. Also,
the inset diagram shows the variance of the distribution as we repeat the
experiment 10 times for each algorithm. We added the results for random
initialization for each step of the iteration of the LNP - initialized 10
different networks, trained them, and showed their performance in 10(c). After
each iteration, we randomly created a network with an equal number of synapses
and edges as found by the LNP method and trained and tested it on the Lorenz63
dataset to get the RMSE loss. We see that the Random initialized network shows
higher variance and shows more jumps signifying the randomly initialized model
is unstable without proper finetuning. In addition, we see the performance
consistently getting worse as the model size keeps decreasing, signifying an
optimal network architecture is more crucial for smaller networks than for
larger networks. We see that the final distribution for both synapses and
neurons of the AP-based models has a higher variance than the LNP algorithm.
This also highlights the stability of the proposed LNP algorithm. In addition
to this, we also plot the final distributions of the timescales observed from
the two methods.

Datasets: We evaluate the performance of the LNP pruning methods for (1) time-
series prediction on chaotic systems (i.e., Lorenz [78] and Rossler systems
[79]) and two real-world datasets - Google Stock Price Prediction and wind
speed prediction datasets [80] and (2) image classification on CIFAR10 and
CIFAR100 datasets [81].

Baselines:  We use the activity-based pruning (AP pruning) method [82] as the
baseline, where we prune the neurons with the lowest activations in each
iteration. This algorithm operates iteratively, pruning the least active
neurons and retraining the model in each iteration. We also compare the LNP
algorithm with current task-dependent state-of-the-art pruning algorithms from
prior works, and the results are shown in Table IV. In addition, we also
introduce the Random initialization method, where after each iteration of the
LNP algorithm, we observe the number of neurons and synapses of the LNP pruned
network and then generate a random Erdos-Renyi graph with the same number of
neurons and synapses. We repeat this process for each step of the iteration.
The results of the Random initialization method are shown in Fig. 10(c).

Evaluation Metric First, we use the standard RMSE loss which is given as
$\displaystyle\operatorname{RMSE}(t)=\sqrt{\frac{1}{D}\sum_{i=1}^{D}\left[\frac{u_{i}^{f}(t)-u_{i}(t)}{\sigma_{i}}\right]^{2}}$
where $D$ is the system dimension, $\sigma$ is the long term standard
deviation of the time series, $\epsilon$ is an arbitrary threshold, and
$u^{f}$ is the forecast. We also use another measure to measure the
performance of prediction called valid prediction time (VPT) [83]. The VPT is
the time $t$ when the accuracy of the forecast exceeds a given threshold. Thus
$\displaystyle VPT(t)=\sum\mathbb{I}(RMSE(t)<\epsilon)$ For these
experiments,we set $\epsilon$ arbitrarily to 0.1. Thus, a higher VPT indicates
a better prediction model.

Energy Efficiency: In assessing the energy consumption of neuromorphic chips,
the central measure is the energy required for a single spike to pass through
a synapse, a notably energy-intensive process [84]. The overall energy
consumption of a Spiking Neural Network (SNN) can be estimated by counting the
synaptic operations (SOPs), analogous to floating-point operations (FLOPs) in
traditional Artificial Neural Networks (ANNs). The energy consumption of an
SNN is calculated as $\displaystyle E=C_{E}\times\text{Total
SOPs}=C_{E}\sum_{i}s_{i}c_{i}$, where $C_{E}$ represents the energy per SOP,
and $\text{Total SOPs}=\sum_{i}s_{i}c_{i}$ is the sum of SOPs. Each
presynaptic neuron $i$ fires $s_{i}$ spikes, connecting to $c_{i}$ synapses,
with every spike contributing to one SOP. Again, for sparse SNNs, the energy
model is reformulated as [85]:
$E=C_{E}\sum_{i}\left(s_{i}\sum_{j}n_{i}^{\text{pre}}\wedge\theta_{ij}\wedge
n_{ij}^{\text{post}}\right)$, where the set
$\left(n_{i}^{\text{pre}},\theta_{ij},n_{ij}^{\text{post}}\right)\in\\{0,1\\}^{3}$
indicates the state of the presynaptic neuron, the synapse, and the
postsynaptic neuron. A value of 1 denotes an active state, while 0 indicates
pruning. The $\wedge$ symbol represents the logical AND operation. Hence, we
calculate the "SOP Ratio" between the unpruned and pruned networks as a metric
for comparison of the energy efficiency of the pruning methods, which
quantifies the energy savings relative to the original, fully connected
(dense) network. This ratio provides a meaningful way to gauge the efficiency
improvements in sparse SNNs compared to their dense counterparts.

TABLE IV: Comparison of Pruning Methods(*=CIFAR10 pruned model trained &
tested on CIFAR100)

|  | CIFAR10 | CIFAR100
---|---|---|---
Method |  | Spiking
---
Model
| Baseline
---
Accuracy
| Accuracy
---
Loss
| Neuron
---
Sparsity
| Synapse
---
Sparsity
| SOP
---
Ratio
| Baseline
---
Accuracy
| Accuracy
---
Loss
| Neuron
---
Sparsity
| Synapse
---
Sparsity
| SOP
---
Ratio
| ADMM
---
[86]
7Conv, 2FC | 89.53 | -3.85 | - | 90 | 2.91 | - | - | - | - | -
| LSNN
---
[87]
6Conv, 2FC | 92.84 | -3.53 | - | 97.96 | 16.59 | - | - | - | - | -
| Grad R
---
[88]
6Conv, 2FC | 92.54 | -0.30 | - | 71.59 | 2.09 | 71.34 | -4.03 | - | 97.65 | 19.45
| IMP
---
[89]
ResNet19 | 93.22 | -0.04 | - | 97.54 | 13.29 | 71.34* | -2.39* | - | 97.54 | 18.69
| STDS
---
[90]
6Conv, 2FC | 92.49 | -0.35 | - | 88.67 | 5.27 | - | - | - | - | -
| ESL-SNN
---
[91]
ResNet19 | 91.09 | -1.7 | - | 95 | 2.11 | 73.48 | -0.99 | - | 95 | 14.22
LNP (ours) | ResNet19 (Random) | 93.29 $\pm$ 0.74 | -2.15 $\pm$ 0.19 | 90.48 | 94.32 | 9.36 $\pm$ 1.28 | 73.32 $\pm$ 0.81 | -3.96 $\pm$ 0.39 | 90.48 | 94.32 | 12.04 $\pm$ 0.41
ResNet19 (converted) |  | -0.04 $\pm$ 0.01 | 93.67 | 98.19 | 22.18 $\pm$ 0.2 |  | -0.11 $\pm$ 0.02 | 94.44 | 98.07 | 31.15 $\pm$ 0.28
HRSNN-STDP | 90.26 $\pm$ 0.88 | -0.76 $\pm$ 0.07 | 94.03 | 95.06 | 39.01 $\pm$ 0.47 | 68.94 $\pm$ 0.73 | -2.16 $\pm$ 0.28 | 92.02 | 94.21 | 48.21 $\pm$ 0.59
HRSNN-BP | 92.37 $\pm$ 0.91 | -0.81 $\pm$ 0.08 | 35.67 $\pm$ 0.41 | 70.12 $\pm$ 0.71 | -2.37 $\pm$ 0.32 | 42.57 $\pm$ 0.63
CHRSNN-STDP | 91.58 $\pm$ 0.83 | -0.62 $\pm$ 0.07 | 95.04 | 96.68 | 50.37 $\pm$ 0.61 | 69.96 $\pm$ 0.68 | -1.65 $\pm$ 0.21 | 93.47 | 97.04 | 57.44 $\pm$ 0.68
CHRSNN-BP | 93.45 $\pm$ 0.87 | -0.74 $\pm$ 0.07 | 45.32 $\pm$ 0.58 | 73.45 $\pm$ 0.66 | -1.11 $\pm$ 0.23 | 50.35 $\pm$ 0.64

Readout Layer Pruning:  The read-out layer is task-dependent and uses
supervised training. In this paper, we do not explicitly prune the read-out
network, but the readout layer is implicitly pruned. The readout layer is a
multi-layer (two or three layers) fully connected network. The first layer
size is equal to the number of neurons sampled from the recurrent layer. We
sample the top 10% of neurons with the greatest betweenness centrality. Thus,
as the number of neurons in the recurrent layer decreases, the size of the
first layer of the read-out network also decreases. The second layer consists
of fixed size with 20 neurons, while the third layer differs between the
classification and the prediction tasks such that for classification, the
number of neurons in the third layer is equal to the number of classes. On the
other hand, the third layer for the prediction task is a single neuron which
gives the prediction output.

TABLE V: Table comparing the performance on the Lorenz 63 and Google datasets.

Pruning | Model | Training | Avg. SOPs | Lorenz63 | Google Dataset
---|---|---|---|---|---
Method |  | Method | (M) | RMSE | VPT | RMSE | VPT
Unpruned | HRSNN | BP | 815.77 $\pm$ 81.51 | 0.248 $\pm$ 0.058 | 44.17 $\pm$ 6.31 | 0.794 $\pm$ 0.096 | 42.18 $\pm$ 6.22
STDP | 710.76 $\pm$ 79.65 | 0.315 $\pm$ 0.042 | 35.75 $\pm$ 4.65 | $0.905\pm 0.095$ | 32.36 $\pm$ 3.14
CHRSNN | BP | 867.42 $\pm$ 93.12 | 0.235 $\pm$ 0.052 | 47.23 $\pm$ 6.02 | 0.782$\pm$0.091 | 45.28$\pm$5.98
STDP | 744.97 $\pm$ 80.09 | 0.285 $\pm$ 0.021 | 40.17 $\pm$ 5.13 | 1.948 $\pm$ 0.179 | 19.25 $\pm$ 3.54
AP Pruned | HRSNN | BP | 117.52 $\pm$ 14.37 | 1.245 $\pm$ 0.554 | 31.08 $\pm$ 8.23 | 1.457$\pm$ 0.584 | 28.24$\pm$6.98
STDP | 92.68 $\pm$ 10.11 | 1.718 $\pm$ 0.195 | 21.10 $\pm$ 7.22 | 1.948 $\pm$ 0.179 | 19.25 $\pm$ 7.59
CHRSNN | BP | 157.33 $\pm$ 18.87 | 1.114 $\pm$ 0.051 | 33.97 $\pm$ 7.56 | 1.325$\pm$ 0.566 | 26.47$\pm$7.42
STDP | 118.77 $\pm$ 10.59 | 1.596 $\pm$ 0.194 | 29.41 $\pm$ 7.33 | 1.987 $\pm$ 0.191 | 17.68 $\pm$ 7.38
LNP Pruned | HRSNN | BP | 22.87 $\pm$ 2.27 | 0.691 $\pm$ 0.384 | 33.67 $\pm$ 6.88 | 0.855$\pm$0.112 | 33.14$\pm$3.01
STDP | 18.22 $\pm$ 2.03 | 0.705 $\pm$ 0.104 | 32.17 $\pm$ 4.62 | 0.917 $\pm$ 0.124 | 30.25 $\pm$ 3.26
CHRSNN | BP | 19.14 $\pm$ 2.04 | 0.682 $\pm$ 0.312 | 39.15 $\pm$ 6.27 | 0.832$\pm$0.105 | 34.51$\pm$2.87
STDP | 14.79 $\pm$ 1.58 | 0.679 $\pm$ 0.098 | 39.24 $\pm$ 4.15 | 0.901 $\pm$ 0.101 | 32.14$\pm$ 3.05

Performance Comparison in Classification: Table IV shows the comparison of our
model with other state-of-the-art pruning algorithms in current literature. It
must be noted here that these algorithms are task-dependent and use only
synapse pruning, keeping the model architecture fixed. We evaluate the model
on the CIFAR10 & CIFAR100 datasets and observe that our proposed task-agnostic
pruning algorithm performs closely to the current state-of-the-art.

Performance Comparison in Prediction: The pruned model derived by the LNP
method is task-agnostic. As such, we can train the model with either STDP or
gradient-based approaches. In this section, we compare the performance of the
unsupervised STDP-trained model with the supervised surrogate gradient method
to train the pruned HRSNN model [92]. First, we plot the evolution of RMSE
loss with pruning iterations in Fig. 10(c) when trained and evaluated on the
Lorenz 63 dataset. We see that our pruning method shows minimal degradation in
performance compared to the baseline activity pruned method. Further on, Table
V presents the comparative performance of different pruning methods, models,
and training methods on the Lorenz 63 and the Google stock prediction
datasets. It is clear from the results that the LNP-pruned models generally
outperform the Unpruned and AP-pruned counterparts across different models and
training methods. Specifically, LNP pruned models consistently exhibit lower
SOPs, indicating enhanced computational efficiency while maintaining
competitive RMSE and VPT values, which indicate the model’s predictive
accuracy and validity, respectively. This suggests that employing the LNP
pruning method can significantly optimize model performance without
compromising the accuracy of predictions.

Ablation Studies: We conducted an ablation study, where we systematically
examined various combinations of the four sequential steps involved in the LNP
method. This study’s findings are presented graphically, as illustrated in
Fig. 12. At each point (A-E), we train the model and obtain the model’s
accuracy and average. Synaptic Operations (SOPs) to support ablation studies.
The ablation study is done for the HRSNN model, which is trained using STDP
and evaluated on the CIFAR10 dataset.

Figure 12: Plot showing Ablation studies of LNP

In the figure, different line styles and colors represent distinct aspects of
the procedure: the blue line corresponds to Steps 1 and 2 of the LNP process,
the orange line to Step 3, and the green line to Step 4. Solid lines depict
the progression of the original LNP process $(A\rightarrow B\rightarrow
C\rightarrow D\rightarrow E)$, while dotted lines represent the ablation
studies conducted at stages B and C. This visual representation enables a
clear understanding of the individual impact each step exerts on the model’s
performance and efficiency. Additionally, it provides insights into potential
alternative outcomes that might have arisen from employing different
permutations of these steps or omitting certain steps altogether.

##  VII Conclusion

This paper consolidates and summarizes a series of research studies that
collectively advance the field of Spiking Neural Networks (SNNs) through the
introduction of heterogeneity in neuron and synapse dynamics, rigorous
analytical modeling, and innovative pruning methodologies.

The first study introduced the concept of Heterogeneous Recurrent Spiking
Neural Networks (HRSNNs) with variability in LIF neuron parameters and STDP
dynamics, demonstrating significant improvements in performance, robustness,
and efficiency in handling complex spatio-temporal datasets for action
recognition tasks. By incorporating heterogeneity, this research highlighted
the potential of smaller, more efficient models with sparse connections that
require less training data.

Building upon this, the second study provided a comprehensive analytical
framework that elucidates the benefits of heterogeneity in enhancing memory
capacity and reducing spiking activity in SNNs. This work established critical
mathematical properties that align with neurobiological observations,
emphasizing the importance of variability in achieving efficient and robust
neural computations.

The final study addressed the computational challenges associated with
heterogeneity by introducing the Lyapunov Noise Pruning (LNP) method. This
novel pruning strategy leverages the inherent heterogeneity to reduce the
number of neurons and synapses, thereby decreasing computational complexity
while maintaining high performance. The task-agnostic nature of LNP ensures
its applicability across diverse scenarios, making it a versatile tool for
optimizing SNNs.

In conclusion, these studies collectively demonstrate that heterogeneity in
SNNs is not only beneficial but essential for developing high-performance,
energy-efficient neural networks. By bridging theoretical insights with
practical applications, this body of work paves the way for future innovations
in neuromorphic computing and related fields. The advancements presented here
underscore the transformative potential of SNNs, offering new avenues for
research and application in machine learning and artificial intelligence.

## Acknowledgement

This work was supported in part by Semiconductor Research Corporation (SRC).
This work is supported by the Army Research Office and was accomplished under
Grant Number W911NF-19-1-0447. The views and conclusions contained in this
document are those of the authors and should not be interpreted as
representing the official policies, either expressed or implied, of the Army
Research Office or the U.S. Government.

## References

  * [1] Yujie Wu, Lei Deng, Guoqi Li, Jun Zhu, and Luping Shi.  Spatio-temporal backpropagation for training high-performance spiking neural networks.  Frontiers in neuroscience, 12:331, 2018.
  * [2] Yingyezhe Jin, Wenrui Zhang, and Peng Li.  Hybrid macro/micro level backpropagation for training deep spiking neural networks.  Advances in neural information processing systems, 31(3):552–561, 2018.
  * [3] Sumit B Shrestha and Garrick Orchard.  Slayer: Spike layer error reassignment in time.  Advances in neural information processing systems, 31:331, 2018.
  * [4] Jesus L Lobo, Javier Del Ser, Albert Bifet, and Nikola Kasabov.  Spiking neural networks and online learning: An overview and perspectives.  Neural Networks, 121:88–100, 2020.
  * [5] Amirhossein Tavanaei, Masoud Ghodrati, Saeed Reza Kheradpisheh, Timothée Masquelier, and Anthony Maida.  Deep learning in spiking neural networks.  Neural Networks, 111:47–63, 2019.
  * [6] Andreea Lazar, Gordon Pipa, and Jochen Triesch.  The combination of stdp and intrinsic plasticity yields complex dynamics in recurrent spiking networks.  In ESANN, volume 1, pages 647–652. Citeseer, 2006.
  * [7] Biswadeep Chakraborty and Saibal Mukhopadhyay.  Heterogeneous recurrent spiking neural network for spatio-temporal classification.  arXiv preprint arXiv:2211.04297, 17(https://doi.org/10.3389/fnins.2023.99451), 2022.
  * [8] Nicolas Perez-Nieves, Vincent CH Leung, Pier Luigi Dragotti, and Dan FM Goodman.  Neural heterogeneity promotes robust learning.  bioRxiv, 12(1):2020–12, 2021.
  * [9] Xueyuan She, Saurabh Dash, and Saibal Mukhopadhyay.  Sequence approximation using feedforward spiking neural network for spatiotemporal learning: Theory and optimization methods.  In International Conference on Learning Representations, volume 32, pages 3028–3039, 2021.
  * [10] Bojian Yin, Federico Corradi, and Sander M Bohte.  Accurate online training of dynamical spiking neural networks through forward propagation through time.  arXiv preprint arXiv:2112.11231, 12(1):1040, 2021.
  * [11] Fleur Zeldenrust, Boris Gutkin, and Sophie Denéve.  Efficient and robust coding in heterogeneous recurrent networks.  PLoS computational biology, 17(4):e1008673, 2021.
  * [12] Biswadeep Chakraborty and Saibal Mukhopadhyay.  Heterogeneous neuronal and synaptic dynamics for spike-efficient unsupervised learning: Theory and design principles.  In Submitted to The Eleventh International Conference on Learning Representations, volume 17, page 994517. Frontiers Media SA, 2023.  Accepted.
  * [13] Pau Vilimelis Aceituno, Gang Yan, and Yang-Yu Liu.  Tailoring echo state networks for optimal learning.  Iscience, 23(9):101440, 2020.
  * [14] Mirko Goldmann, Felix Köster, Kathy Lüdge, and Serhiy Yanchuk.  Deep time-delay reservoir computing: Dynamics and memory capacity.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(9):093124, 2020.
  * [15] Biswadeep Chakraborty, Beomseok Kang, Harshit Kumar, and Saibal Mukhopadhyay.  Sparse spiking neural network: Exploiting heterogeneity in timescales for pruning recurrent SNN.  In The Twelfth International Conference on Learning Representations, 2024.
  * [16] Filip Ponulak and Andrzej Kasinski.  Introduction to spiking neural networks: Information processing, learning and applications.  Acta neurobiologiae experimentalis, 71(4):409–433, 2011.
  * [17] Wulfram Gerstner and Werner M Kistler.  Mathematical formulations of hebbian learning.  Biological cybernetics, 87(5):404–415, 2002.
  * [18] Biswadeep Chakraborty and Saibal Mukhopadhyay.  Characterization of generalizability of spike timing dependent plasticity trained spiking neural networks.  Frontiers in Neuroscience, 15(arXiv preprint arXiv:2105.14677):695357, 2021.
  * [19] Biswadeep Chakraborty and Saibal Mukhopadhyay.  Brain-inspired spiking neural network for online unsupervised time series prediction.  arXiv preprint arXiv:2304.04697, pages 1–8, 2023.
  * [20] Jaehoon Lee, Yasaman Bahri, Roman Novak, Samuel S Schoenholz, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Deep neural networks as gaussian processes.  arXiv preprint arXiv:1711.00165, 2017.
  * [21] Mina A Khoei, Amirreza Yousefzadeh, Arash Pourtaherian, Orlando Moreira, and Jonathan Tapson.  Sparnet: Sparse asynchronous neural network execution for energy efficient inference.  In 2020 2nd IEEE International Conference on Artificial Intelligence Circuits and Systems (AICAS), pages 256–260. IEEE, 2020.
  * [22] Wenrui Zhang and Peng Li.  Temporal spike sequence learning via backpropagation for deep spiking neural networks.  Advances in Neural Information Processing Systems, 33:12022–12033, 2020.
  * [23] Biswadeep Chakraborty, Xueyuan She, and Saibal Mukhopadhyay.  A fully spiking hybrid neural network for energy-efficient object detection.  IEEE Transactions on Image Processing, 30:9014–9029, 2021.
  * [24] Seijoon Kim, Seongsik Park, Byunggook Na, and Sungroh Yoon.  Spiking-yolo: Spiking neural network for energy-efficient object detection.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 11270–11277, 2020.
  * [25] Erol Gelenbe, Zhi-Hong Mao, and Yan-Da Li.  Function approximation with spiked random networks.  IEEE Transactions on Neural Networks, 10(1):3–9, 1999.
  * [26] Nicolangelo Iannella and Andrew D. Back.  A spiking neural network architecture for nonlinear function approximation.  Neural Networks, 14(6):933–939, 2001.
  * [27] Yujie Wu, Lei Deng, Guoqi Li, Jun Zhu, Yuan Xie, and Luping Shi.  Direct training for spiking neural networks: Faster, larger, better.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 1311–1318, 2019.
  * [28] Daehyun Kim, Biswadeep Chakraborty, Xueyuan She, Edward Lee, Beomseok Kang, and Saibal Mukhopadhyay.  Moneta: A processing-in-memory-based hardware platform for the hybrid convolutional spiking neural network with online learning.  Frontiers in Neuroscience, 16:775457, 2022.
  * [29] Gopalakrishnan Srinivasan and Kaushik Roy.  Restocnet: Residual stochastic binary convolutional spiking neural network for memory-efficient neuromorphic computing.  Frontiers in neuroscience, 13:189, 2019.
  * [30] Biswadeep Chakraborty and Saibal Mukhopadhyay.  Topological representations of heterogeneous learning dynamics of recurrent spiking neural networks.  arXiv preprint arXiv:2403.12462, 2024.
  * [31] Biswadeep Chakraborty, Uday Kamal, Xueyuan She, Saurabh Dash, and Saibal Mukhopadhyay.  Brain-inspired spatiotemporal processing algorithms for efficient event-based perception.  In 2023 Design, Automation & Test in Europe Conference & Exhibition (DATE), pages 1–6. IEEE, Publisher Name, 2023.
  * [32] ER De Kloet and JMHM Reul.  Feedback action and tonic influence of corticosteroids on brain function: a concept arising from the heterogeneity of brain receptor systems.  Psychoneuroendocrinology, 12(2):83–105, 1987.
  * [33] Charles Petitpré, Haohao Wu, Anil Sharma, Anna Tokarska, Paula Fontanet, Yiqiao Wang, Françoise Helmbacher, Kevin Yackle, Gilad Silberberg, Saida Hadjab, et al.  Neuronal heterogeneity and stereotyped connectivity in the auditory afferent system.  Nature communications, 9(1):1–13, 2018.
  * [34] Maoz Shamir and Haim Sompolinsky.  Implications of neuronal diversity on population coding.  Neural computation, 18(8):1951–1986, 2006.
  * [35] Biswadeep Chakraborty, Lidia Ghosh, and Amit Konar.  Designing phase-sensitive common spatial pattern filter to improve brain-computer interfacing.  IEEE Transactions on Biomedical Engineering, 67(7):2064–2072, 2019.
  * [36] Biswadeep Chakraborty, Lidia Ghosh, and Amit Konar.  Optimal selection of eeg electrodes using interval type-2 fuzzy-logic-based semiseparating signaling game.  IEEE Transactions on Cybernetics, 51(12):6200–6212, 2020.
  * [37] Minah Lee, Xueyuan She, Biswadeep Chakraborty, Saurabh Dash, Burhan Mudassar, and Saibal Mukhopadhyay.  Reliable edge intelligence in unreliable environment.  In 2021 Design, Automation & Test in Europe Conference & Exhibition (DATE), pages 896–901. IEEE, 2021.
  * [38] Beomseok Kang, Biswadeep Chakraborty, and Saibal Mukhopadhyay.  Unsupervised 3d object learning through neuron activity aware plasticity.  arXiv preprint arXiv:2302.11622, (arXiv preprint arXiv:2302.11622), 2023.
  * [39] Tiago P Carvalho and Dean V Buonomano.  Differential effects of excitatory and inhibitory plasticity on synaptically driven neuronal input-output functions.  Neuron, 61(5):774–785, 2009.
  * [40] Sonja B Hofer, Ho Ko, Bruno Pichler, Joshua Vogelstein, Hana Ros, Hongkui Zeng, Ed Lein, Nicholas A Lesica, and Thomas D Mrsic-Flogel.  Differential connectivity and response dynamics of excitatory and inhibitory neurons in visual cortex.  Nature neuroscience, 14(8):1045–1052, 2011.
  * [41] P Jesper Sjostrom, Ede A Rancz, Arnd Roth, and Michael Hausser.  Dendritic excitability and synaptic plasticity.  Physiological reviews, 88(2):769–840, 2008.
  * [42] Daniel E Feldman.  The spike-timing dependence of plasticity.  Neuron, 75(4):556–571, 2012.
  * [43] Martin Korte and Dietmar Schmitz.  Cellular and system biology of memory: timing, molecules, and beyond.  Physiological reviews, 96(2):647–693, 2016.
  * [44] R Rossi Pool and Germán Mato.  Spike-timing-dependent plasticity and reliability optimization: the role of neuron dynamics.  Neural computation, 23(7):1768–1789, 2011.
  * [45] Robert Legenstein and Wolfgang Maass.  Edge of chaos and prediction of computational performance for neural circuit models.  Neural networks, 20(3):323–334, 2007.
  * [46] Peter I Frazier.  A tutorial on bayesian optimization.  arXiv preprint arXiv:1807.02811, 2018.
  * [47] David Eriksson and Martin Jankowiak.  High-dimensional bayesian optimization with sparse axis-aligned subspaces.  In Uncertainty in Artificial Intelligence, pages 493–503. PMLR, 2021.
  * [48] Jean Feydy, Thibault Séjourné, François-Xavier Vialard, Shun-ichi Amari, Alain Trouvé, and Gabriel Peyré.  Interpolating between optimal transport and mmd using sinkhorn divergences.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2681–2690. PMLR, 2019.
  * [49] Herbert Jaeger.  The “echo state” approach to analysing and training recurrent neural networks-with an erratum note.  Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, 148(34):13, 2001.
  * [50] Herbert Jaeger et al.  Short term memory in echo state networks, volume 5.  GMD-Forschungszentrum Informationstechnik Bremen, Germany, 2001.
  * [51] Fernando Montani, Robin AA Ince, Riccardo Senatore, Ehsan Arabzadeh, Mathew E Diamond, and Stefano Panzeri.  The impact of high-order interactions on the rate of synchronous discharge and information transmission in somatosensory cortex.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367(1901):3297–3310, 2009.
  * [52] Larry F Abbott and Peter Dayan.  The effect of correlated variability on the accuracy of a population code.  Neural computation, 11(1):91–101, 1999.
  * [53] Céline Duval, Eric Luçon, and Christophe Pouzat.  Interacting hawkes processes with multiplicative inhibition.  Stochastic Processes and their Applications, 148:180–226, 2022.
  * [54] Ryan Vogt, Maximilian Puelma Touzel, Eli Shlizerman, and Guillaume Lajoie.  On lyapunov exponents for rnns: Understanding information propagation using dynamical systems tools.  arXiv preprint arXiv:2006.14123, 8:818799, 2020.
  * [55] Ian Covert, Scott M Lundberg, and Su-In Lee.  Understanding global feature contributions with additive importance measures.  Advances in Neural Information Processing Systems, 33:17212–17223, 2020.
  * [56] Yan Zhou, Yaochu Jin, and Jinliang Ding.  Surrogate-assisted evolutionary search of spiking neural architectures in liquid state machines.  Neurocomputing, 406:12–23, 2020.
  * [57] Henry Markram, Joachim Lübke, Michael Frotscher, and Bert Sakmann.  Regulation of synaptic efficacy by coincidence of postsynaptic aps and epsps.  Science, 275(5297):213–215, 1997.
  * [58] Qinyi Wang, Yexin Zhang, Junsong Yuan, and Yilong Lu.  Space-time event clouds for gesture recognition: From rgb cameras to event cameras.  In 2019 IEEE Winter Conference on Applications of Computer Vision (WACV), pages 1826–1835. IEEE, 2019.
  * [59] Yin Bi, Aaron Chadha, Alhabib Abbas, Eirina Bourtsoulatze, and Yiannis Andreopoulos.  Graph-based spatio-temporal feature learning for neuromorphic vision sensing.  IEEE Transactions on Image Processing, 29:9084–9098, 2020.
  * [60] Joao Carreira and Andrew Zisserman.  Quo vadis, action recognition? a new model and the kinetics dataset.  In proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 6299–6308, 2017.
  * [61] Hyunhoon Lee, Young-Seok Kim, Mijung Kim, and Youngjoo Lee.  Low-cost network scheduling of 3d-cnn processing for embedded action recognition.  IEEE Access, 9:83901–83912, 2021.
  * [62] Zijian Wang, Yanting Zhang, Haibo Shi, Lei Cao, Cairong Yan, and Guangwei Xu.  Recurrent spiking neural network with dynamic presynaptic currents based on backpropagation.  International Journal of Intelligent Systems, 2021.
  * [63] Hanle Zheng, Yujie Wu, Lei Deng, Yifan Hu, and Guoqi Li.  Going deeper with directly-trained larger spiking neural networks.  arXiv preprint arXiv:2011.05280, 2020.
  * [64] Guobin Shen, Dongcheng Zhao, and Yi Zeng.  Backpropagation with biologically plausible spatio-temporal adjustment for training deep spiking neural networks.  arXiv preprint arXiv:2110.08858, 2021.
  * [65] Qianhui Liu, Dong Xing, Huajin Tang, De Ma, and Gang Pan.  Event-based action recognition using motion information and spiking neural networks.  In Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21, Z.-H. Zhou, Ed. International Joint Conferences on Artificial Intelligence Organization, volume 8, pages 1743–1749, 2021.
  * [66] Priyadarshini Panda and Narayan Srinivasa.  Learning to recognize actions from limited training examples using a recurrent spiking neural model.  Frontiers in neuroscience, 12:126, 2018.
  * [67] Xueyuan She, Saurabh Dash, Daehyun Kim, and Saibal Mukhopadhyay.  A heterogeneous spiking neural network for unsupervised learning of spatiotemporal patterns.  Frontiers in Neuroscience, 14:1406, 2021.
  * [68] Wei Fang, Zhaofei Yu, Yanqi Chen, Timothée Masquelier, Tiejun Huang, and Yonghong Tian.  Incorporating learnable membrane time constant to enhance learning of spiking neural networks.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 2661–2671, 2021.
  * [69] Yan Meng, Yaochu Jin, and Jun Yin.  Modeling activity-dependent plasticity in bcm spiking neural networks with application to human behavior recognition.  IEEE transactions on neural networks, 22(12):1952–1966, 2011.
  * [70] Vladimir Ivanov and Konstantinos Michmizos.  Increasing liquid state machine performance with edge-of-chaos dynamics organized by astrocyte-modulated plasticity.  Advances in Neural Information Processing Systems, 34, 2021.
  * [71] Jay Patravali, Gaurav Mittal, Ye Yu, Fuxin Li, and Mei Chen.  Unsupervised few-shot action recognition via action-appearance aligned meta-adaptation.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 8484–8494, 2021.
  * [72] Khurram Soomro and Mubarak Shah.  Unsupervised action discovery and localization in videos.  In Proceedings of the IEEE International Conference on Computer Vision, pages 696–705, 2017.
  * [73] Benjamin Cramer, Yannik Stradmann, Johannes Schemmel, and Friedemann Zenke.  The heidelberg spiking data sets for the systematic evaluation of spiking neural networks.  IEEE Transactions on Neural Networks and Learning Systems, 33(7):2744–2757, 2020.
  * [74] Edward N Lorenz.  Predictability: A problem partly solved.  In Proc. Seminar on predictability, volume 1, 1996.
  * [75] Tobias Thornes, Peter Düben, and Tim Palmer.  On the use of scale-dependent precision in earth system modelling.  Quarterly Journal of the Royal Meteorological Society, 143(703):897–908, 2017.
  * [76] Ankita Paul, Stefan Wagner, and Anup Das.  Learning in feedback-driven recurrent spiking neural networks using full-force training.  arXiv preprint arXiv:2205.13585, pages 1–10, 2022.
  * [77] Bojian Yin, Federico Corradi, and Sander M Bohté.  Effective and efficient computation with multiple-timescale spiking recurrent neural networks.  In International Conference on Neuromorphic Systems 2020, pages 1–8, 2020.
  * [78] Meiling Xu, Min Han, CL Philip Chen, and Tie Qiu.  Recurrent broad learning systems for time series prediction.  IEEE transactions on cybernetics, 50(4):1405–1417, 2018.
  * [79] Meiling Xu, Min Han, Tie Qiu, and Hongfei Lin.  Hybrid regularized echo state network for multivariate chaotic time series prediction.  IEEE transactions on cybernetics, 49(6):2305–2315, 2018.
  * [80] Subhrajit Samanta, Mahardhika Pratama, and Suresh Sundaram.  Bayesian neuro-fuzzy inference system for temporal dependence estimation.  IEEE Transactions on Fuzzy Systems, 29(9):2479–2490, 2020.
  * [81] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * [82] Alistair G Rust1, Rod Adams1, Stella George1, and Hamid Bolouri.  Activity-based pruning in developmental artificial neural networks.  In Fourth European Conference on Artificial Life, page 224. MIT Press, 1997.
  * [83] Pantelis R Vlachas, Jaideep Pathak, Brian R Hunt, Themistoklis P Sapsis, Michelle Girvan, Edward Ott, and Petros Koumoutsakos.  Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics.  Neural Networks, 126:191–217, 2020.
  * [84] Steve Furber.  Large-scale neuromorphic computing systems.  Journal of neural engineering, 13(5):051001, 2016.
  * [85] Xinyu Shi, Jianhao Ding, Zecheng Hao, and Zhaofei Yu.  Towards energy efficient spiking neural networks: An unstructured pruning framework.  In The Twelfth International Conference on Learning Representations, 2024.
  * [86] Lei Deng, Yujie Wu, Yifan Hu, Ling Liang, Guoqi Li, Xing Hu, Yufei Ding, Peng Li, and Yuan Xie.  Comprehensive snn compression using admm optimization and activity regularization.  IEEE transactions on neural networks and learning systems, 2021.
  * [87] Guillaume Bellec, Darjan Salaj, Anand Subramoney, Robert Legenstein, and Wolfgang Maass.  Long short-term memory and learning-to-learn in networks of spiking neurons.  Advances in neural information processing systems, 31, 2018.
  * [88] Yanqi Chen, Zhaofei Yu, Wei Fang, Tiejun Huang, and Yonghong Tian.  Pruning of deep spiking neural networks through gradient rewiring.  arXiv preprint arXiv:2105.04916, 2021.
  * [89] Youngeun Kim, Yuhang Li, Hyoungseob Park, Yeshwanth Venkatesha, Ruokai Yin, and Priyadarshini Panda.  Exploring lottery ticket hypothesis in spiking neural networks.  In European Conference on Computer Vision, pages 102–120. Springer, 2022.
  * [90] Yanqi Chen, Zhaofei Yu, Wei Fang, Zhengyu Ma, Tiejun Huang, and Yonghong Tian.  State transition of dendritic spines improves learning of sparse spiking neural networks.  In International Conference on Machine Learning, pages 3701–3715. PMLR, 2022.
  * [91] Jiangrong Shen, Qi Xu, Jian K Liu, Yueming Wang, Gang Pan, and Huajin Tang.  Esl-snns: An evolutionary structure learning strategy for spiking neural networks.  arXiv preprint arXiv:2306.03693, 2023.
  * [92] Emre O Neftci, Hesham Mostafa, and Friedemann Zenke.  Surrogate gradient learning in spiking neural networks.  IEEE Signal Processing Magazine, 36:61–63, 2019.