# The fine line between dead neurons and sparsity
in binarized spiking neural networks

Jason K. Eshraghian  Wei D. Lu

###### Abstract

Spiking neural networks can compensate for quantization error by encoding
information either in the temporal domain, or by processing discretized
quantities in hidden states of higher precision. In theory, a wide dynamic
range state-space enables multiple binarized inputs to be accumulated
together, thus improving the representational capacity of individual neurons.
This may be achieved by increasing the firing threshold, but make it too high
and sparse spike activity turns into no spike emission. In this paper, we
propose the use of ‘threshold annealing’ as a warm-up method for firing
thresholds. We show it enables the propagation of spikes across multiple
layers where neurons would otherwise cease to fire, and in doing so, achieve
highly competitive results on four diverse datasets, despite using binarized
weights. Source code is available at https://github.com/jeshraghian/snn-tha/.

Machine Learning, ICML

##  1 Introduction

Memory access and data communication remain amongst the most expensive
instructions in deep learning workloads (Hashemi et al., 2018; Kozyrakis et
al., 2010). To complement Feynman111And his American Physical Society lecture
titled ‘There is Plenty of Room at the Bottom’ (Feynman, 1959)., there is
plenty of room for optimization at the top of the stack: learning algorithms
can be tweaked to decrease memory access frequency using activation sparsity,
and to reduce memory usage by weight quantization.

Activation sparsity: spiking neural networks (SNNs) are an example of models
that leverage activation sparsity ${\boldsymbol{z}\in\\{0,1\\}^{n}}$, where an
activation ${z^{i}=0}$ allows input-weight-multiply steps to be bypassed as
the product is guaranteed to also return ‘0’. Memory access to read out
weights can be skipped (Roy et al., 2019; Eshraghian et al., 2021).

Reduced weight precision: parameter quantization directly reduces the number
of memory cells that must be read at any compute cycle. In the extreme,
binarized neural networks (BNNs) constrain each weight to
$w^{ij}\in\\{-1,1\\}$ at the cost of reduced representational capacity and
challenging learning convergence (Hubara et al., 2016; Courbariaux et al.,
2016).

(a)

(b)

(c)

Figure 1: (a) Neuron response to the three possible post-synaptic currents
$x^{j}$: (i) $x^{j}=0$: No response if $w^{ij}=0$, (ii) $x^{j}=-1$: Negative
deflection with short-term memory but no spiking activity, (iii) $x^{j}=+1$:
Spike propagation but no short-term memory due to reset. (b) High dynamic
range (HDR) state precision comes at the expense of suppressed network
activity. (c) Gradual evolution of the dynamic range (DR) balances precision
and spike emission. Threshold annealing follows a similar principle, but the
threshold $\theta$ is increased over training iterations instead of time.
Spiking neuron parameters are thus fixed at run time, reducing computation and
memory costs.

SNNs can amortize the cost of binarization by: 1) distributing information in
the temporal domain rather than the activation amplitude (Voelker et al.,
2020; Kheradpisheh et al., 2021), and 2) accumulating binarized data in the
high-precision, dynamically evolving state-space of spiking neurons (Schaefer
& Joshi, 2020; Laborieux et al., 2021). In both cases, discrete activations
are applied to a continuous domain (either time or neuron state).

The dynamic range of the neuron state can be increased by using a large firing
threshold. This is crucial for binarized SNNs (BSNNs), as it extends the range
over which discrete inputs may be accumulated. However, indefinitely
increasing the threshold impedes neuronal activity (and therefore, learning),
leading to the ‘dead neuron problem’. To offset this issue, we propose
‘threshold annealing’, a warm-up technique that improves classification
accuracy in supervised learning tasks across four datasets (two static, two
dynamic). This is achieved by dynamically increasing the threshold outside of
the training-loop which enables:

  * •

early epochs to achieve sufficient network activity for low-precision learning
to take place (i.e., by using a small threshold to promote adequate neuronal
firing),

  * •

later epochs to utilize the wide dynamic range state-space of spiking neurons
for fine-tuning at increased precision (i.e., large thresholds that enable
accumulation of multiple binarized quantities).

Our findings show state-of-the-art performance on all four datasets when
compared to similarly lightweight networks, and a reduction of dead neurons on
the most challenging dataset by 71% when compared to BSNNs without threshold
annealing, leading to faster training convergence.

##  2 Background

###  2.1 Spiking Neuron Model

SNNs adopt the same topologies modeled by conventional neural networks, but
with the artificial neuron traded for a spiking neuron model. A single-state
leaky integrate-and-fire neuron derived using the forward Euler method
(Appendix A) governed by the following discrete-time dynamics is used
(Lapique, 1907; Gerstner & Kistler, 2002):

$u^{j}_{t+1}=\beta u^{j}_{t}+\sum_{i}w^{ij}z^{i}_{t+1}-z_{t}^{j}\theta,$ (1)

$z^{j}_{t}=\begin{cases}1,&\text{\rm if $u^{j}_{t}\geq\theta$}\\\
0,&\text{otherwise},\\\ \end{cases}$ (2)

where $u^{j}_{t}$ is the hidden state (membrane potential) of neuron $j$ at
time $t$; $\beta$ is the decay rate of membrane potential; $w^{ij}$ is the
weight between neurons $i$ and $j$; and a spike $z_{t}^{j}$ is elicited if the
potential exceeds the threshold $\theta$. The final term of Equation 1 resets
the state by subtracting the threshold $\theta$ each time an output spike is
generated.

###  2.2 Binarized SNNs

In BSNNs, each input spike $z^{i}_{t+1}\in\\{0,1\\}$ is scaled by a binarized
weight $w^{ij}\in\\{-1,1\\}$, resulting in a trinary post-synaptic potential
$x^{j}_{t+1}=w^{ij}z_{t+1}^{i}\in\\{-1,0,1\\}$. This is added to the hidden
state of a downstream neuron $u^{j}_{t}$, which triggers its own spike given
sufficient excitation $u^{j}\geq\theta$. With threshold
normalization,222Normalizing $\theta=1$ is common practice as it both enables
the surrogate gradient term $\partial z^{j}/\partial u^{j}$ to be interpreted
as a firing probability (Shrestha & Orchard, 2018; Zenke & Vogels, 2021), and
is also simple (Booij & tat Nguyen, 2005). the three possible outcomes while
the membrane potential is at rest are depicted in Figure 1(a). In case (i), a
nulled activation means no state change. In case (ii), the neuron is inhibited
and the hidden state follows a continuous trajectory back to its resting
value. This is the only case where memory dynamics persist. In case (iii), the
state immediately reaches the firing threshold $u^{j}=\theta$ upon spike
arrival, and is subsequently forced to reset $u^{j}\leftarrow 0$. This is the
only case where information is propagated to deeper layers. A problem emerges:
short-term memory and spike emission are mutually exclusive.

Using a large threshold is the most obvious solution (Lu & Sengupta, 2020),
and achieves several benefits illustrated in Figure 1(b):

  * •

the dynamic range of the state-space is increased,

  * •

multiple inputs may be accumulated, and

  * •

excitatory inputs can now be sustained as short-term, recurrent memory.

But these benefits do not indefinitely hold - too high a threshold obstructs
the state from ever reaching it, causing spiking sparsity to extinguish into
dead neurons. Large thresholds also increase susceptibility to vanishing
gradients from an effect that is independent of dead neurons; surrogate
gradients are typical activation functions (e.g., sigmoid, tanh, fast-sigmoid)
centered about the firing threshold. Shift the surrogate function too far, and
many neurons will dominantly operate in the subthreshold regime where
$\partial z^{i}/\partial u^{i}\rightarrow 0$.

Another drawback of using large thresholds is that it may interfere with
gradient-based methods that rely on low thresholds to promote more neuronal
activity, which is often the case for spike-time derivative methods and
single-spike training schemes (Bohte et al., 2002; Booij & tat Nguyen, 2005;
Xu et al., 2013). These conflicting trade-offs must somehow be balanced.

###  2.3 Related Work

Inspired by mechanisms that promote stability in the primary sensory cortex,
threshold adaptation has been incorporated into computational models of
synaptic plasticity (Bienenstock et al., 1982; Zhang & Linden, 2003).
Adaptation has been used as a homeostatic mechanism to regulate network
activity in deep learning-inspired SNNs trained via local learning rules
(Diehl & Cook, 2015) and gradient-based optimization (Bellec et al., 2018;
Shaban et al., 2021). In all previous cases, threshold adaptation evolves over
the course of simulated time steps (Figure 1(c)). It has shown to be a
potential approach to retaining long-term memory dynamics in SNNs using
biologically plausible means, with the cost of each neuron updating its
threshold trajectory at each time step. While this scales memory complexity
linearly with the number of neurons $\mathcal{O}(n)$ during the forward-pass,
the threshold must be recorded at every time step during training to implement
the backpropagation (BPTT) algorithm, which additionally scales memory usage
with time as well $\mathcal{O}(nT)$ (Rumelhart et al., 1986; Werbos, 1990).
This is overhead that we aim to eliminate in a compute and memory bound
system, such as a BSNN, by i) moving the threshold evolution outside of the
weight update loop, and ii) using a shared threshold for each layer.

By eliminating the temporal variation of the threshold (i.e., it is only
varied between training iterations), the purposes of membrane potential and
threshold can be better decoupled. The membrane potential can be thought of as
the mechanism for temporal data retention, while threshold warm-up leads to an
evolving dynamic range state-space which enables inter-layer spike
propagation.

##  3 Threshold Annealing

Algorithm 1 Threshold Annealing

0: $\theta_{0}\in\mathbb{R}_{>0}$: Initial threshold

0: $\theta_{\infty}\in\mathbb{R}_{>\theta_{0}}$: Final threshold

0: $\alpha_{\theta}\in\mathbb{R}_{>0}:$ Inverse threshold time constant

0: $w:$ Initial parameters

0: $\eta:$ Learning rate

$\theta_{\gamma}\leftarrow\theta_{0}$ (Initialize threshold)

$s\leftarrow\partial\tilde{z}(\theta_{0})/\partial u$ (Initialize surrogate
gradient using $\theta_{0}$)

$\gamma\leftarrow 0$ (Initialize training step)

while $w$ not converged do

$\gamma\leftarrow\gamma+1$

$g_{\gamma}\leftarrow\nabla_{w}f_{\gamma}(s,w_{\gamma-1})$ (Get gradients
w.r.t. surrogate gradient $s$ using initial threshold $\theta_{0}$)

$w_{\gamma}\leftarrow w_{\gamma-1}-\eta\cdot g_{\gamma}$ (Update parameters
e.g. using SGD or Adam)

$\theta_{\gamma}\leftarrow\theta_{\gamma-1}+\alpha(\theta_{\infty}-\theta_{\gamma-1})$
(Update threshold without updating $s$)

end while

return: $w_{\gamma}$ (Resulting parameters)

Out-of-the-loop threshold updates reduces complexity: See Algorithm 1 for
pseudo-code of our proposed threshold annealing method. Appendix B provides
more detailed pseudo-code that highlights how temporal evolution is isolated
from spiking dynamics. A neuron’s firing threshold exponentially relaxes from
an initial value $\theta_{0}$ to a larger final value $\theta_{\infty}$ with
time constant $\tau_{\theta}$:

$\tau_{\theta}\frac{d\theta_{\gamma}}{d\gamma}=-\theta_{\gamma}+\theta_{\infty},\quad\theta_{\infty}>\theta_{\gamma}$
(3)

where all neurons in each layer share the following explicit update rule at
each training iteration $\gamma$:

$\theta_{\gamma+1}=\theta_{\gamma}+\alpha(\theta_{\infty}-\theta_{\gamma})$
(4)

and $\alpha=\text{exp}(-1/\tau_{\theta})$ is the inverse time constant of the
threshold, noting that time here refers to the span of training iterations
rather than the temporal dynamics of a spiking neuron. By only updating the
threshold outside of the temporal evolution of each forward-pass (i.e.,
$\theta_{\gamma}$ is fixed during inference), memory no longer scales with
time, and it also does not scale with iterations $\gamma$ as the BPTT
algorithm truncates the computational graph between iterations. For shared
thresholds per layer, memory scales with $\mathcal{O}(L)<<\mathcal{O}(nT)$,
where $L$ is the number of layers in the network.

Fast to slow evolution balances weight updates with threshold updates: The
threshold dynamics are designed to follow the evolution of the effective
learning rate. For example, when using the Adam optimizer, the estimation of
moments (and therefore, weight updates) exponentially decays over iterations
(Kingma & Ba, 2014). Aligning the rate of threshold and weight updates is
desirable as it avoids either of the two having an overpowering effect on the
network: excessively large threshold updates will clip too many spikes at the
start of training (similarly to Figure 1(b)), whereas exceedingly small
threshold updates fails to sufficiently extend the dynamic range of the state-
space (the problems depicted in Figure 1(a)).

To ensure that a threshold update does not lead to the elimination of any
spikes whilst still expanding the state-space of a neuron, the threshold
change $\Delta\theta_{\gamma}=\theta_{\gamma+1}-\theta_{\gamma}$ should be
bounded by the distance between the new membrane potential at iteration
$\gamma+1$ from the original threshold at $\gamma$:

$\Delta\theta_{\gamma}<u_{\gamma+1}-\theta_{\gamma}$ (5)

Figure 2: Membrane potential $u^{j}$ changes across iterations as a function
of weight updates. If the threshold update across training iterations
$\Delta\theta_{\gamma}$ is greater than the change of membrane potential, then
a spike from iteration $\gamma$ will no longer occur at $\gamma+1$. Imposing
the constraint from Equation 5 prevents this. Note: the reset mechanism is not
depicted for visual clarity.

This is visually depicted in Figure 2, where neglecting the condition from
Equation 5 would suppress spikes. If this condition is neglected for all of
time $t$, the neuron would cease to contribute to the loss function, and lack
a gradient signal to contribute to learning, i.e., a dead neuron. From the
derivation provided in Appendix C, the following two constraints can be
imposed on $\theta_{\gamma}$ and $\theta_{\infty}$ to prevent threshold
annealing from contributing to dead neurons:

$(1-\alpha)\theta_{\gamma}+\alpha\theta_{\infty}<\mathbf{z}\cdot~{}\Big{(}\mathbf{w}_{\gamma}^{~{}:,j}+\eta\frac{\partial\mathcal{L}}{\partial\mathbf{w}^{~{}:,j}_{\gamma}}\Big{)}$
(6)

$\theta_{\infty}<\mathbf{z}\cdot~{}\mathbf{w}_{\infty}^{~{}:,j}$ (7)

where $\mathbf{z}\in\mathbb{R}^{i}$ is the input vector, and
$\mathbf{w}_{\gamma}^{~{}:,j}\in\mathbb{R}^{i}$ is the weight vector connected
to the $j^{th}$ neuron of the first layer, and $\mathbf{w}_{\infty}^{~{}:,j}$
is the final weight vector at the end of training.

The weight update term on the right-hand side of Equation 6 assumes stochastic
gradient descent (SGD) is used, where for adaptive optimizers the derivative
term often decays exponentially with network converges. As $\theta_{\infty}$
is a constant term, the evolution of $\theta_{\gamma}$ should also
exponentially diminish. This is adhered to in Equation 4, and balances the two
objectives of i) using a large threshold, while ii) not inducing dead neurons
due to threshold updates.

More intuitively, during early training iterations, the weight updates are
large enough to ‘catch up’ to the effect of a large threshold step. This is
because the effect of weight decay is minimal during early training
iterations. As the network is trained for more epochs, threshold updates are
reduced to match that of weight decay.

In the experiments that follow, we adhere to exponential decay of threshold
update magnitude as per Equation 4 and choose $\theta_{\infty}$,
$\theta_{\gamma}$ and $\alpha$ based on a hyperparameter search. For sparsely
populated input data, i.e., where the mean of the input features $\mathbf{z}$
over time is quite small (e.g., the Spiking Heidelberg Digits dataset (Cramer
et al., 2020)), the optimal value of $\alpha$ reduces by several orders of
magnitude to slow down the rate at which $\theta_{\gamma}$ converges to
$\theta_{\infty}$ (see Section D.3). This result is independently predicted by
Equation 6.

Clamped surrogate gradients boost the weight update signal for subthreshold
neurons:  The conventional method to account for the non-differentiability of
Equation 2 is to use a surrogate gradient which substitutes the Heaviside
operator with a threshold-shifted, differentiable alternative e.g., the
sigmoid function, in the backward pass (see Figure 3). If this approach was
strictly adhered to, then the surrogate function would shift along with the
threshold over training iterations.

Consider a neuron with a subthreshold membrane potential at point A, depicted
in Figure 3. Shifting the surrogate function along with the threshold causes
$s^{j}_{t}\rightarrow 0$, i.e., the vanishing gradient problem (contrasted
from the dead neuron problem, where neurons cease to fire).

Figure 3: Top row: output spiking activity $z^{j}_{t}$ is represented with the
Heaviside operator centered about the firing threshold. Over iterations,
threshold annealing pushes the Heaviside operator along the x-axis. Middle
row: Surrogate of the Heaviside operator to address non-differentiability
$\tilde{z}^{j}_{t}$. Our approach does not shift the surrogate function which
reduces the impact of vanishing gradients on subthrehsold neurons (e.g., at
point A). Bottom row: the surrogate gradient $s^{j}_{t}$. The intersection
between a neuron at point A and the non-shifted gradient is larger than that
of the shifted gradient, thus providing latent/quiet neurons with the
opportunity to continue learning, and therefore contribute to network
activity.

(a)

(b)

(c)

(d)

Figure 4: An intuitive demonstration of threshold annealing on a spike timing
task. (a) Full precision weights with threshold $\theta=1$. (b) Binarized
weights with threshold $\theta=1$ leads to unstable neuronal activity. (c)
Binarized weights with a large threshold $\theta=50$ suppresses all firing
activity; the increase in hidden state $u$ is due to an increasing bias. (d)
Binarized weights with threshold annealing from $\theta:5\rightarrow 50$
successfully learns to fire at the desired spike-time at the 400th iteration.
The performance better matches that of the full precision network. Animated
visualizations are provided at the following link: github.com/jeshraghian/snn-
tha/.

To address this issue, our implementation clamps the surrogate gradient to the
initial threshold value $\theta_{0}$, preventing the function from shifting
during the training process. This is highlighted in Algorithm 1, where the
threshold $\theta_{\gamma}$ is updated without updating the surrogate gradient
$s$. Our experiments use the fast sigmoid function:

$s^{j}_{t}=\frac{\partial z^{j}_{t}}{\partial
u^{j}_{t}}\leftarrow\frac{1}{(1+k|\theta_{0}-u^{j}_{t}|)^{2}}$ (8)

where the left arrow denotes substitution, and $k$ is a slope that modulates
the linearity of the surrogate function. While the learning process is
generally robust to the choice of the surrogate gradient (Zenke & Vogels,
2021), we choose to fix the gradient based on the initial threshold
$\theta_{0}$ as it removes the need to re-compute a new surrogate function at
each iteration.

The primary benefit is that the gradient signal is boosted for neurons closer
to their resting potential. Quantitatively, for a neuron at rest
$u_{t}^{j}\approx 0$ and $\theta>>1$, and assuming $k=1$, then
$s_{t}^{j}\approx 1/\theta^{2}$ (based on Equation 8). In practice,
$\theta_{\infty}$ can be close to one order of magnitude larger than
$\theta_{0}$ (see Section D.1 and Section D.3 for examples). If
$\theta_{\gamma}$ is used in Equation 8 instead, this would reduce $s_{t}^{j}$
(and by extension, $\nabla_{w}\mathcal{L}$) by the square of this factor as
$\theta_{\gamma}\rightarrow\theta_{\infty}$.

This benefit does not come for free, as active neurons in later training
iterations (where $u^{j}_{t}\approx\theta_{\infty}$) face the same problem of
a vanishing $s_{t}^{j}$. We claim this is tolerable for two empirical reasons:
i) upon spike emission, the membrane potential of an active neuron is reset
where a low $u^{j}$ provides an opportunity for the neuron to learn at another
time step, and ii) the neuron must have already been trained in earlier
iterations for its membrane potential to reach a large value $\theta_{\infty}$
(see Figure 4(d) for an experiment where $u^{j}$ starts at rest and is trained
to approach $\theta_{\infty}=50$). Quiet neurons are therefore provided with a
better opportunity to contribute to network activity via gradient clamping.

##  4 Experimental Results

We test a full precision SNN and several BSNNs under the following cases i)
with a normalized threshold $\theta=1$, ii) with a fixed threshold determined
by a hyperparameter search, and iii) using threshold annealing, for each
experiment on four datasets ranging from simple to increasingly difficult:
MNIST (LeCun et al., 1998), FashionMNIST (Xiao et al., 2017), DVS128 Gesture
(Amir et al., 2017), and Spiking Heidelberg Digits (SHD) (Cramer et al.,
2020), without dataset augmentation. All experiments were performed on a 16GB
NVIDIA V100 GPU, conducted across 5 trials each. Hyperparameters were selected
based on a sweep over 500 trials for each individual experiments (see Appendix
D for further details). We used snnTorch 0.4.11 to construct spiking neuron
models (Eshraghian et al., 2021) which uses a PyTorch 1.10.1 backend (Paszke
et al., 2019) in Python 3.7.9.

###  4.1 Temporal Coding

Before moving to more complex data-driven tasks, we first provide an
interpretable demonstration of how threshold annealing can improve BSNN
performance on a straightforward spike-time encoding task. A single output
neuron is trained to fire at a target time step by encouraging the membrane
potential to linearly increase to the threshold at $t=75$, and subsequently
remain quiet upon spike emission. A 3-layer BSNN is used with a fully-
connected architecture of 100-1000-1 neurons. The input to the network is a
Poisson spike train simulated across 100 time steps. The mean square error
between the target and actual membrane potential are summed at each time step.
Precise implementation details are provided in Section D.1.

Figure 4 depicts the membrane potential trace of the output neuron over time
at several training iterations $\gamma$. The full precision network in (a)
accomplishes the task with ease, while the BSNN with $\theta=1$ in (b)
highlights the challenges described in Section 2.2, i.e., the neuron is unable
to both emit spikes as well as retain memory. Increasing the threshold of all
neurons to $\theta=50$ completely suppresses spiking activity in the hidden
and output layers. The membrane potential increases due to the bias, and
settles at approximately the midpoint between steady-state and threshold,
where the loss of the triangular membrane potential target is minimized. To
fix the memory leakage problem in (b) and the dead neuron problem in (c),
threshold annealing is applied by gradually warming up $\theta_{0}=5$ to
$\theta_{\infty}=50$ with an inverse threshold time constant $\alpha=5\times
10^{-3}$ (Equation 4). The low threshold at the first iteration $\gamma=0$
enables spike propagation throughout layers, which sustains sufficient
activity to keep the network learning over epochs. At the 400th iteration, the
BSNN with threshold annealing successfully emits a spike at the target time
step. The output neuron continued to remain active at the end of the
simulation, where $\theta_{\gamma}$ was within 0.1% of the $\theta_{\infty}$;
a definitive improvement over use of a fixed threshold in Figure 4(c).
Animations in the source code offer further visual intuition on how threshold
annealing ‘pulls up’ neuronal activity.

###  4.2 Static Datasets

The MNIST (LeCun et al., 1998) and the slightly more challenging FashionMNIST
datasets (Xiao et al., 2017) are used to assess the performance of threshold
annealing on temporally static data. Both datasets consist of 60,000
28$\times$28 greyscale images in the training set, and 10,000 images in the
test set, with 10 output classes of handwritten digits (MNIST) and clothing
items/accessories (FashionMNIST). The convolutional BSNN architecture is
provided in the LABEL:tab:static footnote: 16Conv5 refers to a $5\times 5$
convolutional kernel with $16$ filters; $AP2$ is a $2\times 2$ average-pooling
layer applied to membrane potential, and 1024Dense10 is a dense layer with
1024 inputs and 10 outputs.

The raw datasets are passed to the BSNN for 100 time steps of simulation
without input encoding. The predicted class is set to the neuron with the
highest spike count at the end of the 100 time steps. Hyperparameters and the
optimization process for each experiment are provided in Section D.2, and the
final results are shown in LABEL:tab:static, where our proposed threshold
annealing technique improves the performance of BSNNs by an average of 0.3%
for MNIST, and 0.97% for FashionMNIST over the next best BSNN. The worst-case
scenario of a BSNN with a normalized threshold still performs reasonably well,
on account of static datasets not requiring any memory retention. Performance
degrades further when temporal dynamics are included.

Table 1: Test set accuracies for static datasets.

| MNIST | FashionMNIST
---|---|---
Best | Mean | $\sigma$ | Best | Mean | $\sigma$
flt32 | 99.45 | 99.31 | 0.12 | 91.13 | 91.03 | 0.09
BSNN,$\theta$:1 | 98.46 | 98.44 | 0.02 | 87.11 | 86.70 | 0.36
BSNN | 98.81 | 98.77 | 0.06 | 87.01 | 86.95 | 0.10
Proposed | 99.12 | 99.07 | 0.04 | 88.12 | 87.92 | 0.22

  * 1

16Conv5-AP2-64Conv5-AP2-1024Dense10.

  * 2

$\sigma$: Sample standard deviation.

Table 2: Test set accuracies for temporal datasets.

| DVS128 Gesture | SHD
---|---|---
Best | Mean | $\sigma$ | Best | Mean | $\sigma$
flt32 | 93.05 | 92.87 | 0.32 | 82.51 | 82.27 | 0.27
BSNN,$\theta$:1 | 84.37 | 84.03 | 0.60 | 63.86 | 62.99 | 0.91
BSNN | 89.24 | 88.54 | 0.70 | 62.94 | 62.48 | 0.41
Proposed | 92.36 | 91.32 | 0.56 | 66.35 | 65.68 | 0.66

  * 1

DVS128: 16Conv5-AP2-32Conv5-AP2-8800Dense11.

  * 2

SHD: 700Dense3000Dense20.

  * 3

$\sigma$: Sample standard deviation.

###  4.3 Temporal Datasets

The DVS128 Gesture dataset (Amir et al., 2017) is filmed with an event-based
camera (Patrick et al., 2008) which only processes sufficient changes in
luminance, and consists of 11 different output classes of hand gestures, such
as clapping, arm rotation, and air guitar. To fit the dataset within 100 time
steps, events were integrated over a period of 5 ms such that each sample in
the training and test set are approximately 1 s and 3.5 s in duration,
respectively. Spatial downsampling fits each sample to dimensions of $(2\times
32\times 32)$.

As a more challenging alternative, the SHD dataset was also used which
consists of 20 output classes of audio recordings of the ten digits spoken
aloud in both the English and German languages. The recordings are processed
using a cochlear model for spike encoding (Cramer et al., 2020).

The final test accuracies are shown in Table 2, with experimental and
hyperparameter details in Section D.3. As dataset complexity increases, the
performance of BSNNs degrade further, while threshold annealing offsets the
degradation to a larger degree. For the DVS128 Gesture dataset, when compared
to the full precision network, the baseline BSNNs with and without threshold
normalization drop by 8.8% and 4.4%, respectively. But this performance drop
is reduced significantly to 1.5% when threshold annealing is used during
training. These results indicate that BSNNs have an optimal solution that does
not deviate far from the high precision case; the challenge is finding this
solution. Methods such as threshold annealing make the training process
easier.

Performance drops drastically for the SHD dataset when binarizing weights
across all cases. While threshold annealing improves upon performance by 2.7%
compared to the next best BSNN, it still remains at a large distance from the
high precision performance. While this may partially be accounted for by
complexity of the SHD dataset, it may also be due to the additional
hyperparameters that are introduced which require a larger number of trials
during the hyperparameter search process. For fairness, the number of random
trials was set to 500 for all networks, though threshold annealing would
likely benefit from increasing the search space.

##  5 Discussion

(a)

(b)

Figure 5: Spiking activity across layers with threshold annealing (THA) and
without (BSNN) on temporal datasets. Static dataset distributions are provided
in Section E.1. (a) DVS128 Gesture: Distribution of spikes in the output layer
with threshold annealing is balanced, whereas most neurons are silent without
annealing. (c) SHD: Using a dataset with sparse activity results in a large
number of dead without threshold annealing. Figure 6: Loss curves on dynamic
datasets. Oscillations are a result of using periodic learning rate schedule.

###  5.1 Comparison with lightweight SNNs

There is a general lack of reproducible baselines that implement gradient-
based learning on convolutional BSNNs, likely due to hyperparameter
sensitivity and the associated difficulty in training. The few related
approaches include Lu and Sengupta’s ANN-SNN conversion which is benchmarked
on large-scale static datasets (Lu & Sengupta, 2020). Kheradpisheh et al. use
temporal backpropagation (BS4NN) on dense BSNNs (Kheradpisheh et al., 2021).
We further consider several additional promising lightweight approaches to
training full precision SNNs, including sparse spiking gradient descent (SSGD)
which reduces overhead during gradient calculations (Perez-Nieves & Goodman,
2021), and neural heterogeneity (NH) which compresses the required number of
neurons with the addition of neuron-independent parameters (Perez-Nieves et
al., 2021). Interestingly, BSNNs with threshold annealing can outperform full
precision SNNs in terms of accuracy using considerably less memory for model
parameters. These baselines rely on dense networks only, which is why we
developed our own full precision and BSNN baselines in previous sections.

FashionMNIST: SSGD and NH achieve 86.7% and 87.5%, respectively, using full
precision weights. NH requires 101k full precision weights (3.2 Mbits). BS4NN
achieves 87.3% using 794 Kbits of binarized parameters. Threshold annealing
outperforms in terms of accuracy at 87.9% with far less memory footprint with
12k binarized weights (12 Kbits).

DVS128 Gesture: NH obtains 82.1%, whereas our BSNN approach obtains 91.3% with
less memory overhead. Assuming the input dimensions of the DVS128 Gesture
dataset are downsampled by a factor of 4, the network in NH requires 8.4 Mbits
of full precision weights as against 98 Kbits of binarized weights in our
approach.

SHD: In contrast to the above results, both full precision NH (82.7%) and SSGD
(77.5%) outperform our BSNN (65.7%). This can be attributed to the use of
recurrent connections in NH and SSGD; without recurrent SNNs in SSGD,
performance dropped below 50%. Likewise from Cramer et al.’s original proposal
of the SHD dataset, the test accuracy was 48.6% on a similarly 2-layer dense
network (Cramer et al., 2020), which our BSNN with threshold annealing
outperforms by around 17%. When recurrent connections were added, test
accuracy jumped to 83.2%. Although we use a considerably wider hidden layer
(3000 neurons vs 128 neurons in NH), binarized weights still occupy less
memory (ours: 2.16 Mbits vs NH: 2.95 Mbits).

###  5.2 Spiking Activity and Convergence Rate

Figure 5 shows the distribution of spiking activity across layers both with
and without threshold annealing. This is measured by passing the entire test
set to the network and taking the average spike count for each neuron. For the
fixed threshold network trained on the SHD dataset, 55% and 90% of neurons in
the first and final layer do not contribute to network activity, compared to
approximately 20% with threshold annealing in both layers. This is highly
undesirable at the start of the learning process as gradients cannot
backpropagate through neurons that do not contribute to the final loss, which
causes slower training convergence. This is shown in Figure 6, where the
larger proportion of dead neurons in the fixed threshold network trained on
the SHD dataset leads to worse performance. Averaging across both layers,
threshold annealing reduces the number of dead neurons by 71% in the SHD
network.

The threshold dynamics of the final layer of both networks are also overlayed,
where the rate of increase is slower when trained on SHD than DVS128 Gesture.
This is a result of sparser input activity, and was predicted by Equation 7.
Once the training process is terminated based on early stopping criteria, the
same thresholds of the optimal trial should be used. In all experiments, the
final threshold $\theta_{\gamma}$ was always within 0.1% of the steady-state
threshold $\theta_{\infty}$, so we simply use $\theta_{\infty}$ at test time.

###  5.3 Concluding remarks

In general, our approach to training BSNNs is extremely competitive when
compared to high precision, lightweight networks that consume more memory. As
our implementation is based on a PyTorch backend, the binarized weights are
represented as full precision values and are not explicitly stored with
compressed memory. Rather, the models generated are intended to be ported onto
dedicated neuromorphic hardware that is tailored for single-bit weight
representations (Orchard et al., 2021; Akbarzadeh-Sherbaf et al., 2020; Rahimi
Azghadi et al., 2020).

It is also important to account for how this memory is used. For example,
while the memory overhead of synaptic weights scale with the number of neurons
as $\mathcal{O}(n^{2})$, SNNs in general are power efficient as this memory
can be accessed with low frequency. On the other hand, neuron specific
parameters, such as in neural heterogeneity (Perez-Nieves et al., 2021) and
threshold adaptation (Bellec et al., 2018), have less memory complexity
$\mathcal{O}(n)$, but this memory must be accessed at each time step to update
the hidden state and does not benefit from spike sparsity. By sharing
thresholds for neurons in given layers, threshold annealing not only adds
negligible memory footprint (i.e., one full precision value per layer), but it
does so with good performance. Combining threshold annealing with heterogenous
time constants and recurrent connections could attract the benefits of both
worlds in lightweight BSNN computation, especially where intrinsic circuit-
level variations can be exploited as diverse neuronal behavior (Kang et al.,
2021).

## Acknowledgements

This work was supported, in part, by the Semiconductor Research Corporation
and DARPA, through the Applications Driving Architectures Research Center, and
the National Science Foundation through awards ECCS-191550 and CCF-1900675.
Jason K. Eshraghian is supported by the Forrest Research Fellowship through
the Forrest Research Foundation.

## References

  * Akbarzadeh-Sherbaf et al. (2020) Akbarzadeh-Sherbaf, K., Safari, S., and Vahabie, A.-H.  A digital hardware implementation of spiking neural networks with binary FORCE training.  _Neurocomputing_ , 412:129–142, 2020.
  * Akiba et al. (2019) Akiba, T., Sano, S., Yanase, T., Ohta, T., and Koyama, M.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pp. 2623–2631, 2019.
  * Amir et al. (2017) Amir, A. et al.  A low power, fully event-based gesture recognition system.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 7243–7252, 2017.
  * Azghadi et al. (2020) Azghadi, M. R., Lammie, C., Eshraghian, J. K., Payvand, M., Donati, E., Linares-Barranco, B., and Indiveri, G.  Hardware implementation of deep network accelerators towards healthcare and biomedical applications.  _IEEE Transactions on Biomedical Circuits and Systems_ , 14(6):1138–1159, 2020.
  * Bellec et al. (2018) Bellec, G., Salaj, D., Subramoney, A., Legenstein, R., and Maass, W.  Long short-term memory and learning-to-learn in networks of spiking neurons.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pp. 795–805, 2018.
  * Bergstra et al. (2011) Bergstra, J., Bardenet, R., Bengio, Y., and Kégl, B.  Algorithms for hyper-parameter optimization.  _Advances in Neural Information Processing Systems_ , 24, 2011.
  * Bienenstock et al. (1982) Bienenstock, E. L., Cooper, L. N., and Munro, P. W.  Theory for the development of neuron selectivity: orientation specificity and binocular interaction in visual cortex.  _Journal of Neuroscience_ , 2(1):32–48, 1982\.
  * Bohte et al. (2002) Bohte, S. M., Kok, J. N., and La Poutre, H.  Error-backpropagation in temporally encoded networks of spiking neurons.  _Neurocomputing_ , 48(1-4):17–37, 2002.
  * Booij & tat Nguyen (2005) Booij, O. and tat Nguyen, H.  A gradient descent rule for spiking neurons emitting multiple spikes.  _Information Processing Letters_ , 95(6):552–558, 2005.
  * Courbariaux et al. (2016) Courbariaux, M., Hubara, I., Soudry, D., El-Yaniv, R., and Bengio, Y.  Binarized neural networks: Training deep neural networks with weights and activations constrained to +1 or -1.  _arXiv preprint arXiv:1602.02830_ , 2016.
  * Cramer et al. (2020) Cramer, B., Stradmann, Y., Schemmel, J., and Zenke, F.  The Heidelberg spiking data sets for the systematic evaluation of spiking neural networks.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020\.
  * Diehl & Cook (2015) Diehl, P. U. and Cook, M.  Unsupervised learning of digit recognition using spike-timing-dependent plasticity.  _Frontiers in Computational Neuroscience_ , 9:99, 2015.
  * Eshraghian et al. (2021) Eshraghian, J. K., Ward, M., Neftci, E., Wang, X., Lenz, G., Dwivedi, G., Bennamoun, M., Jeong, D. S., and Lu, W. D.  Training spiking neural networks using lessons from deep learning.  _arXiv preprint arXiv:2109.12894_ , 2021.
  * Feynman (1959) Feynman, R. P.  Plenty of room at the bottom.  In _APS Annual Meeting_ , 1959.
  * Gerstner & Kistler (2002) Gerstner, W. and Kistler, W. M.  _Spiking neuron models: Single neurons, populations, plasticity_.  Cambridge University Press, 2002.
  * Hashemi et al. (2018) Hashemi, M., Swersky, K., Smith, J., Ayers, G., Litz, H., Chang, J., Kozyrakis, C., and Ranganathan, P.  Learning memory access patterns.  In _International Conference on Machine Learning_ , pp. 1919–1928. PMLR, 2018.
  * He et al. (2019) He, T., Zhang, Z., Zhang, H., Zhang, Z., Xie, J., and Li, M.  Bag of tricks for image classification with convolutional neural networks.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 558–567, 2019.
  * Hubara et al. (2016) Hubara, I., Courbariaux, M., Soudry, D., El-Yaniv, R., and Bengio, Y.  Binarized neural networks.  _Advances in Neural Information Processing Systems_ , 29, 2016.
  * Kang et al. (2021) Kang, S. M., Choi, D., Eshraghian, J. K., Zhou, P., Kim, J., Kong, B.-S., Zhu, X., Demirkol, A. S., Ascoli, A., Tetzlaff, R., et al.  How to build a memristive integrate-and-fire model for spiking neuronal signal generation.  _IEEE Transactions on Circuits and Systems I: Regular Papers_ , 68(12):4837–4850, 2021.
  * Kheradpisheh et al. (2021) Kheradpisheh, S. R., Mirsadeghi, M., and Masquelier, T.  BS4NN: Binarized spiking neural networks with temporal coding and learning.  _Neural Processing Letters_ , pp. 1–19, 2021.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kozyrakis et al. (2010) Kozyrakis, C., Kansal, A., Sankar, S., and Vaid, K.  Server engineering insights for large-scale online services.  _IEEE Micro_ , 30(4):8–19, 2010.
  * Laborieux et al. (2021) Laborieux, A., Ernoult, M., Hirtzlin, T., and Querlioz, D.  Synaptic metaplasticity in binarized neural networks.  _Nature Communications_ , 12(1):1–12, 2021.
  * Lapique (1907) Lapique, L.  Recherches quantitatives sur l’excitation electrique des nerfs traitee comme une polarization.  _Journal of Physiology and Pathology_ , 9:620–635, 1907\.
  * LeCun et al. (1998) LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * Loshchilov & Hutter (2017) Loshchilov, I. and Hutter, F.  SGDR: Stochastic gradient descent with warm restarts.  _International Conference on Learning Representations_ , 2017.
  * Lu & Sengupta (2020) Lu, S. and Sengupta, A.  Exploring the connection between binary and spiking neural networks.  _Frontiers in Neuroscience_ , 14:535, 2020.
  * Orchard et al. (2021) Orchard, G., Frady, E. P., Rubin, D. B. D., Sanborn, S., Shrestha, S. B., Sommer, F. T., and Davies, M.  Efficient neuromorphic signal processing with loihi 2.  In _2021 IEEE Workshop on Signal Processing Systems (SiPS)_ , pp. 254–259. IEEE, 2021.
  * Paszke et al. (2019) Paszke, A. et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in Neural Information Processing Systems_ , 32:8026–8037, 2019.
  * Patrick et al. (2008) Patrick, L., Posch, C., and Delbruck, T.  A 128$\times$128 120 db 15$\mu$ s latency asynchronous temporal contrast vision sensor.  _IEEE Journal of Solid-State Circuits_ , 43:566–576, 2008\.
  * Perez-Nieves & Goodman (2021) Perez-Nieves, N. and Goodman, D. F.  Sparse spiking gradient descent.  _arXiv preprint arXiv:2105.08810_ , 2021.
  * Perez-Nieves et al. (2021) Perez-Nieves, N., Leung, V. C., Dragotti, P. L., and Goodman, D. F.  Neural heterogeneity promotes robust learning.  _Nature Communications_ , 12(5791), 2021.
  * Rahimi Azghadi et al. (2020) Rahimi Azghadi, M., Chen, Y.-C., Eshraghian, J. K., Chen, J., Lin, C.-Y., Amirsoleimani, A., Mehonic, A., Kenyon, A. J., Fowler, B., Lee, J. C., et al.  Complementary metal-oxide semiconductor and memristive hardware for neuromorphic computing.  _Advanced Intelligent Systems_ , 2(5):1900189, 2020.
  * Roy et al. (2019) Roy, K., Jaiswal, A., and Panda, P.  Towards spike-based machine intelligence with neuromorphic computing.  _Nature_ , 575(7784):607–617, 2019.
  * Rumelhart et al. (1986) Rumelhart, D. E., Hinton, G. E., and Williams, R. J.  Learning representations by back-propagating errors.  _Nature_ , 323(6088):533–536, 1986.
  * Schaefer & Joshi (2020) Schaefer, C. J. and Joshi, S.  Quantizing spiking neural networks with integers.  In _International Conference on Neuromorphic Systems 2020_ , pp. 1–8, 2020.
  * Shaban et al. (2021) Shaban, A., Bezugam, S. S., and Suri, M.  An adaptive threshold neuron for recurrent spiking neural networks with nanodevice hardware implementation.  _Nature Communications_ , 12(1):1–11, 2021.
  * Shrestha & Orchard (2018) Shrestha, S. B. and Orchard, G.  SLAYER: Spike layer error reassignment in time.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pp. 1419–1428, 2018.
  * Voelker et al. (2020) Voelker, A. R., Rasmussen, D., and Eliasmith, C.  A spike in performance: Training hybrid-spiking neural networks with quantized activation functions.  _arXiv preprint arXiv:2002.03553_ , 2020.
  * Werbos (1990) Werbos, P. J.  Backpropagation through time: What it does and how to do it.  _Proceedings of the IEEE_ , 78(10):1550–1560, 1990.
  * Xiao et al. (2017) Xiao, H., Rasul, K., and Vollgraf, R.  Fashion-MNIST: A novel image dataset for benchmarking machine learning algorithms.  _arXiv preprint arXiv:1708.07747_ , 2017.
  * Xu et al. (2013) Xu, Y., Zeng, X., Han, L., and Yang, J.  A supervised multi-spike learning algorithm based on gradient descent for spiking neural networks.  _Neural Networks_ , 43:99–113, 2013.
  * Zenke & Vogels (2021) Zenke, F. and Vogels, T. P.  The remarkable robustness of surrogate gradient learning for instilling complex function in spiking neural networks.  _Neural Computation_ , 33(4):899–925, 2021.
  * Zhang & Linden (2003) Zhang, W. and Linden, D. J.  The other side of the engram: experience-driven changes in neuronal intrinsic excitability.  _Nature Reviews Neuroscience_ , 4(11):885–900, 2003.

##  Appendix A Spiking Neuron Model

The dynamics of the passive membrane modeled using an RC circuit can be
represented as:

$\tau_{m}\frac{du_{t}}{dt}=-u_{t}+i_{t}r$ (9)

where $u_{t}$ is the hidden state (membrane potential) of the neuron at time
$t$, $\tau_{m}=rc$ is the time constant of the membrane, $r$ and $c$ are the
resistance and capacitance of the membrane, respectively, and $i_{t}$ is the
input current injection at time $t$.

The goal is to convert Equation 9 into a recurrent representation with fewer
hyperparameters to enable training using backprop through time. The forward
Euler method is used to find an approximate solution to Equation 9 without
taking the limit $\Delta t\rightarrow 0$:

$\tau_{m}\frac{u_{t+\Delta t}-u_{t}}{\Delta t}=-u_{t}+i_{t}r$ (10)

For sufficiently small $\Delta t$, this provides a reasonable approximation of
continuous-time integration. Isolating $u_{t+\Delta t}$ gives:

$u_{t+\Delta t}=(1-\frac{\Delta t}{\tau_{m}})u_{t}+\frac{\Delta
t}{\tau_{m}}i_{t}r$ (11)

The change of $u_{t}$ from one time step to the next can be found in absence
of input current:

$\beta=\frac{u_{t+\Delta t}}{u_{t}}=(1-\frac{\Delta t}{\tau_{m}})$ (12)

To reduce the number of free variables in Equation 11, normalize $\Delta t=1$,
and let $r=1\Omega$:

$u_{t+1}=\beta u_{t}+(1-\beta)i_{t}$ (13)

To further simplify the model, let the effect of $(1-\beta)i_{t}$ be
approximated by a learnable weight scaled by the output of the previous layer,
$(1-\beta)i_{t}\leftarrow\sum_{i}w^{ij}z^{i}_{t}$. For computational
efficiency (i.e., to enable inputs to be directly propagated from the input to
the output of the network without delay), a time-shift is applied to
$z_{t}^{i}$:

$u^{j}_{t+1}=\beta u_{t}^{j}+\sum_{i}w^{ij}z_{t+1}^{i}$ (14)

Finally, the term $-z_{t}^{j}\theta$ is added to Equation 14 to account for
the neuronal reset mechanism where the threshold is subtracted from the
membrane potential each time a spike $z_{t}^{j}=1$ occurs:

$u^{j}_{t+1}=\beta u_{t}^{j}+\sum_{i}w^{ij}z_{t+1}^{i}-\theta z_{t}^{j}.$ (15)

##  Appendix B Threshold Annealing Algorithm

Algorithm 2 Threshold Annealing: Same as Algorithm 1, but with temporal
iteration decoupled from minibatch iteration

0: $\theta_{0}\in\mathbb{R}_{>0}$: Initial threshold

0: $\theta_{\infty}\in\mathbb{R}_{>\theta_{0}}$: Final threshold

0: $\alpha_{\theta}\in\mathbb{R}_{>0}:$ Inverse threshold time constant

0: $w:$ Initial parameters

0: $\eta:$ Learning rate

$\theta_{\gamma}\leftarrow\theta_{0}$ (Initialize threshold)

$s\leftarrow\partial\tilde{z}(\theta_{0})/\partial u$ (Initialize surrogate
gradient using $\theta_{0}$)

$\gamma\leftarrow 0$ (Initialize training step)

while $w$ not converged do

$\gamma\leftarrow\gamma+1$

for $t\in\\{1,...,N\\}$ do

$\mathcal{L}\leftarrow f(w_{\gamma-1},\theta_{\gamma})$ (Get loss by iterating
through $N$ time-steps during forward-pass)

$g_{\gamma}\leftarrow\nabla_{w}\mathcal{L}$ (Get gradients w.r.t. surrogate
gradient $s$ using initial threshold $\theta_{0}$)

end for

$w_{\gamma}\leftarrow w_{\gamma-1}-\eta\cdot g_{\gamma}$ (Update parameters
e.g. using SGD or Adam)

$\theta_{\gamma}\leftarrow\theta_{\gamma-1}+\alpha(\theta_{\infty}-\theta_{\gamma-1})$
(Update threshold outside of forward-pass temporal loop, without updating $s$)

end while

return: $w_{\gamma}$ (Resulting parameters)

Algorithm 2 is a more detailed form of Algorithm 1 to highlight how the
temporal-loop during the forward-pass is decoupled from the minibatch training
loop. Homeostasis would place the threshold update $\theta_{\gamma}$ within
the for-loop to stabilize neuronal activity, whereas threshold annealing takes
it out to increase the effective hidden state precision whilst reducing both
memory and computational complexity.

##  Appendix C Upper Bound on the Rate of Threshold Annealing

To minimize the impact threshold updates have on a network’s ability to learn,
we set the constraint that the threshold update cannot exceed the weight
update-induced membrane potential change. Formally, the threshold update
$\Delta\theta_{\gamma}=\theta_{\gamma+1}-\theta_{\gamma}$ is bounded by how
far the new membrane potential at $\gamma+1$ is from the original threshold at
$\gamma$:

$\Delta\theta_{\gamma}<u_{\gamma+1}-\theta_{\gamma}$ (16)

For simplicity, we assume this takes place at the first spike time for an
instantaneous input (which is why notation for time $t$ has been dropped),
though this constraint will be relaxed.
$\Delta\theta_{\gamma}=\theta_{\gamma+1}-\theta_{\gamma}$ can be found from
Equation 4:

$\Delta\theta_{\gamma}=\alpha(\theta_{\infty}-\theta_{\gamma})$ (17)

which is substituted into Equation 16:

$\alpha(\theta_{\infty}-\theta_{\gamma})<u_{\gamma+1}-\theta_{\gamma}$ (18)

Considering only one spiking neuron, the immediate post-synaptic current
influence at the spike time from Equation 1 for iterations $\gamma$ and
$\gamma+1$ is given below:

$u_{\gamma}=\mathbf{z}\cdot\mathbf{w}_{\gamma}^{~{}:,j}$ (19)

$u_{\gamma+1}=\mathbf{z}\cdot\mathbf{w}_{\gamma+1}^{~{}:,j}$ (20)

where $\mathbf{z}\in\mathbb{R}^{i}$ is the input vector,
$\mathbf{w}_{\gamma}^{~{}:,j}\in\mathbb{R}^{i}$ is the weight vector connected
to the $j^{th}$ neuron of the first layer (and can be generalized to
subsequent layers). For simplicity, assume the weight update is calculated via
stochastic gradient descent (SGD), though noting that the final term below
should be substituted for the preferred optimization update rule:

$\mathbf{w}^{~{}:,j}_{\gamma+1}=\mathbf{w}_{\gamma}^{~{}:,j}+\eta\frac{\partial\mathcal{L}}{\partial\mathbf{w}^{~{}:,j}_{\gamma}}$
(21)

Equation 21 is substituted back into Equation 20:

$u_{\gamma+1}=\mathbf{z}\cdot(\mathbf{w}_{\gamma}^{~{}:,j}+\eta\frac{\partial\mathcal{L}}{\partial\mathbf{w}^{~{}:,j}_{\gamma}})$
(22)

which is then substituted into Equation 18:

$\alpha(\theta_{\infty}-\theta_{\gamma})<\mathbf{z}\cdot(\mathbf{w}_{\gamma}^{~{}:,j}+\eta\frac{\partial\mathcal{L}}{\partial\mathbf{w}^{~{}:,j}_{\gamma}})-\theta_{\gamma}$
(23)

$\implies\underbrace{(1-\alpha)}_{\rm
const.}\theta_{\gamma}+\underbrace{\alpha\theta_{\infty}}_{\rm
const.}<\underbrace{\mathbf{z}}_{\rm
const.}\cdot~{}\Big{(}\mathbf{w}_{\gamma}^{~{}:,j}+\underbrace{\eta\frac{\partial\mathcal{L}}{\partial\mathbf{w}^{~{}:,j}_{\gamma}}}_{\rm
update}\Big{)}$ (24)

All constant terms have been underlined, and the evolution of
$\mathbf{w}^{~{}:,j}_{\gamma}$ is dependent on
$\partial\mathcal{L}/\partial\mathbf{w}^{~{}:,j}_{\gamma}$. When using more
common optimizers, such as Adam or SGD with momentum, the weight update term
must be swapped out. For the Adam optimizer, the update term will follow an
exponential decay as the network approaches training convergence until the
right hand side of Equation 24 stabilizes at
$\mathbf{z}\cdot\mathbf{w}_{\gamma}^{~{}:,j}$. The only reasonable and valid
solution as $\gamma\rightarrow\infty$ is for $\theta_{\gamma}$ to
exponentially relax to $\theta_{\infty}$, where the update term approaches
zero where adaptive/exponential decay is applied, and the term on the left
factors out to:

$\theta_{\infty}<\mathbf{z}\cdot~{}\mathbf{w}_{\infty}^{~{}:,j}$ (25)

$\theta_{\gamma}$ and $\theta_{\infty}$ taking on values much less than the
dot product of the input and weights would be a mathematically valid solution,
but at the cost of each neuron no longer having a sufficiently wide state-
space and the ability to trigger spikes (see Section 2.2 for detail on why
memory and spike emission become mutually exclusive as the threshold is
reduced). Therefore, an exponential relaxation to a steady-state threshold
described in Equation 3 is the optimal solution for exponentially adaptive
weight decays.

Furthermore, $\mathbf{z}$ is a time-varying input (it is only constant across
iterations given it represents the same input data), and so Equation 25 must
be generalized to account for the historical contributions of post-synaptic
current to the membrane state at the spike time (Equation 1).

The above conditions were derived assuming that the threshold update should
not impact spiking activity at all, which is unnecessarily stringent. As long
as each neuron fires at least once, then the immediate weights attached to its
input will contribute to the output loss, have a gradient signal (for dense
layers), and therefore avoid becoming dead neurons. A more reasonable
objective would to be to only introduce the constraint from Equation 25 when a
neuron is emitting at most a single spike across the simulation run time.

##  Appendix D Experimental Results

Parameter Initialization:  For all experiments, parameters were initialized
using the default methods in PyTorch 1.10.1. Dense layers were initialized by
uniformly sampling from U($-\sqrt{a}$, $\sqrt{a}$):

$a=\frac{1}{N_{\rm in}},$ (26)

where $N_{\rm in}$ is the number of input features. Convolutional layer
parameters were also uniformly sampled:

$a=\frac{1}{C_{\rm in}N_{x}N_{y}}$ (27)

where $C_{\rm in}$ is the number of input channels, and $N_{x}$ and $N_{y}$
are the kernel dimensions.

Hyperparameter search:  The following hyperparameters were searched in all
experiments:

  * •

$\beta$: neuron decay rate/membrane potential inverse time constant (Equation
1)

  * •

$\theta$: firing threshold (Equation 2)

  * •

$k$: surrogate gradient slope (Equation 8)

  * •

$\eta$: initial learning rate

  * •

GC: gradient clipping

  * •

WC: weight clipping

  * •

DO: dropout rate (dense layers)

  * •

BN: batch normalization (convolutional layers)

The following hyperparameters were searched in all threshold annealing
experiments:

  * •

$\theta_{0}^{L}$: initial firing threshold of all neurons in layer $L$
(Equation 4)

  * •

$\theta_{\infty}^{L}$: final firing threshold of all neurons in layer $L$
(Equation 4)

  * •

$\alpha^{L}$: inverse time constant of threshold of all neurons in layer $L$
(Equation 4)

Hyperparameters were selected based on 500 separate trials for each experiment
using a tree-structured Parzen Estimator algorithm to randomly sample from the
search space in Optuna (Akiba et al., 2019; Bergstra et al., 2011). Based on
this approach, individual hyperparameters were found for high precision, BSNN,
BSNN with normalized threshold, and BSNN with threshold annealing across all
four datasets. I.e., 16 different sets of hyperparameters were found.

With respect to the threshold:

  * •

High precision and BSNN experiments set the threshold as a free
hyperparameter,

  * •

BSNN with normalized threshold fixes the threshold at $\theta=1$, and

  * •

BSNN with threshold annealing set the initial threshold $\theta_{0}$, final
threshold $\theta_{\infty}$, and inverse time constant $\alpha$ as free
hyperparameters.

Due to the additional hyperparameters in the BSNN with threshold annealing
experiment, the final result would likely benefit from using additional trials
during the hyperparameter search. But for fairness, this was fixed to 500
trials as with all other cases.

Optimization:  Unless otherwise specified, the Adam optimizer was used with
1st and 2nd order moment estimates set to ($\beta_{1},\beta_{2}$) =
($0.9,0.999$). A cosine annealed learning rate with warm restarts is computed
using the following scheduler (Loshchilov & Hutter, 2017; He et al., 2019):

$\eta_{t}=\frac{1}{2}\eta\Big{(}1+\text{cos}\Big{(}\frac{\pi\gamma}{T}\Big{)}\Big{)},$
(28)

where $\eta$ is the initial learning rate, $\gamma$ is the iteration, and $T$
is the period of the schedule. For all cases where cosine annealing is used,
$T$ is set to a period of 10 training epochs.

The fast sigmoid surrogate gradient was used for all spiking neurons (Equation
8).

###  D.1 Temporal Coding

A corresponding notebook to replicate these experiments, in addition to
animated visualizations of the threshold evolution and its impact on membrane
potential, are provided in the source code.

Architecture:  A 3-layer fully connected network was used with 100$-$1000$-$1
neurons.

Hyperparameters:  The hyperparameters of the temporal coding task from Section
4.1 are shown in the table below, and were chosen arbitrarily for the sake of
an intuitive demonstration of threshold annealing.

Table 3: Hyperparameters for Temporal Coding Task

Precision | $\beta$ | $\theta$ | $\theta_{0}$ | $\theta_{\infty}$ | $\alpha$ | $k$ | $\eta$
---|---|---|---|---|---|---|---
flt32 | 0.6 | 2.0 | $-$ | $-$ | $-$ | 5.0 | 1e-3
BSNN; $\theta$=1 | 0.15 | 1.0 | $-$ | $-$ | $-$ | 5.0 | 1e-3
BSNN; $\theta$=50 | 0.15 | 50.0 | $-$ | $-$ | $-$ | 5.0 | 1e-3
Proposed | 0.15 | $-$ | 5.0 | 50.0 | 5e-3 | 5.0 | 1e-3

  * 1

Stochastic Gradient Descent with momentum of 0.9 was used for all temporal
coding experiments. No schedule was used.

Loss:  Let $y_{t}$ be the target membrane potential of the output neuron at
time $t$. The total mean square error is calculated by summing the loss across
all time steps $T$:

$\mathcal{L}_{MSE}=\sum_{t}^{T}(y_{t}-u_{t})^{2}.$ (29)

Peak membrane target:  The peak of the target membrane potential in Figure 4
is set to be 10% above the threshold to increase the probability of spike
emission.

###  D.2 Static Datasets

Hyperparameters:  The hyperparameters used in the static dataset experiments
from Section 4.2 are shown in the table below.

Table 4: Hyperparameters for Static Datasets

Dataset | Precision | $\beta$ | $\theta$ | $k$ | $\eta$ | GC | WC | BN | DO
---|---|---|---|---|---|---|---|---|---
MNIST | flt32 | 0.92 | 2.0 | 6.0 | 1.9e-3 | ✓ | ✓ | ✓ | 0.09
MNIST | BSNN;$\theta$:1 | 0.74 | 1.0 | 4.7 | 6.5e-3 | ✗ | ✗ | ✓ | 0.87
MNIST | BSNN | 0.84 | 6.4 | 2.8 | 2.7e-3 | ✗ | ✗ | ✓ | 0.06
MNIST | Proposed | 0.99 | $-$ | 10.2 | 1.0e-2 | ✗ | ✗ | ✓ | 0.03
FMNIST | flt32 | 0.39 | 1.5 | 7.7 | 2.0e-3 | ✓ | ✓ | ✓ | 0.13
FMNIST | BSNN;$\theta$:1 | 0.82 | 1.0 | 4.3 | 1.2e-3 | ✓ | ✓ | ✓ | 0.86
FMNIST | BSNN | 0.87 | 9.7 | 6.7 | 1.6e-2 | ✓ | ✗ | ✓ | 0.01
FMNIST | Proposed | 0.87 | $-$ | 0.2 | 8.4e-4 | ✗ | ✗ | ✓ | 0.65

  * 1

16Conv5-AP2-64Conv5-AP2-1024Dense10.

  * 2

Batch size of 128 used for all experiments.

  * 3

MNIST flt32 did not use a learning rate schedule.

  * 4

FMNIST Proposed used SGD with momentum set to 0.86.

Table 5: Threshold Annealing Hyperparameters for Static Datasets

Dataset | $\theta_{0}^{1}$ | $\theta_{\infty}^{1}$ | $\alpha^{1}$ | $\theta_{0}^{2}$ | $\theta_{\infty}^{2}$ | $\alpha^{2}$ | $\theta_{0}^{3}$ | $\theta_{\infty}^{3}$ | $\alpha^{3}$
---|---|---|---|---|---|---|---|---|---
MNIST | 11.7 | 16.0 | 2.4e-2 | 14.1 | 30.4 | 0.12 | 0.7 | 4.2 | 1.1e-3
FMNIST | 6.9 | 14.0 | 3.7e-2 | 10.3 | 23.1 | 0.30 | 18.0 | 27.9 | 0.10

  * 1

$\theta_{0}^{L}$ is the initial threshold of all neurons in layer $L$.

  * 2

$\theta_{\infty}^{L}$ is the final threshold of all neurons in layer $L$.

  * 3

$\alpha^{L}$ is the inverse time constant of the threshold for all neurons in
layer $L$.

Loss:  To emulate the limitations of efficient neuromorphic hardware, we
assume the membrane potential is an inaccessible hidden state when calculating
the loss function (Azghadi et al., 2020). Instead, we train our network under
the more challenging constraint of using target spike counts for each class.
The spikes of each output neuron $z^{j}_{t}$ are accumulated over time. The
mean square error loss of the target count $c^{j}$ is calculated and then
summed across the $N$ output classes ($N=10$).

Let $y_{t}$ be the target membrane potential of the output neuron at time $t$.
The total mean square error is calculated by summing the loss across all time
steps $T$. We follow a similar approach to SLAYER (Shrestha & Orchard, 2018)
where the target of the correct class $c$ is to fire 80% of the time, and
incorrect classes are set to fire 20% of the time to avoid excessive
suppression of neuronal activity:

$\mathcal{L}_{MSE}=\sum_{j}^{N}\sum_{t}(c^{j}-z^{j}_{t})^{2}.$ (30)

###  D.3 Temporal Datasets

Hyperparameters:  The hyperparameters used in the temporal dataset experiments
from Section 4.3 are shown in the table below.

Table 6: Hyperparameters for Temporal Datasets

Dataset | Precision | $\beta$ | $\theta$ | $k$ | $\eta$ | GC | WC | BN | DO1 | DO2
---|---|---|---|---|---|---|---|---|---|---
DVS128 | flt32 | 0.72 | 2.5 | 5.4 | 7.5e-3 | ✗ | ✓ | ✓ | 0.29 | $-$
DVS128 | BSNN;$\theta$:1 | 0.1 | 1.0 | 0.1 | 3.6e-3 | ✗ | ✗ | ✗ | 0.64 | $-$
DVS128 | BSNN | 0.78 | 8.7 | 0.3 | 9.6e-4 | ✓ | ✓ | ✗ | 0.06 | $-$
DVS128 | Proposed | 0.93 | $-$ | 0.2 | 1.8e-3 | ✓ | ✗ | ✗ | 0.43 | $-$
SHD | flt32 | 0.98 | 1.7 | 2.0 | 2.7e-4 | ✗ | ✗ | $-$ | 0.05 | 0.04
SHD | BSNN;$\theta$:1 | 0.61 | 1.0 | 0.4 | 1.9e-4 | ✓ | ✓ | $-$ | 0.09 | 0.11
SHD | BSNN | 0.33 | 3.4 | 0.1 | 5.0e-4 | ✓ | ✓ | $-$ | 0.22 | 0.07
SHD | Proposed | 0.95 | $-$ | 0.3 | 6.5e-4 | ✓ | ✓ | $-$ | 0.02 | 0.19

  * 1

DVS128: 16Conv5-AP2-32Conv5-AP2-800Dense11

  * 2

DVS128 flt32 batch size: 16. DVS128 BSNN batch size: 8.

  * 3

DVS128 events integrated over 5 ms.

  * 4

SHD: 700Dense3000Dense20

  * 5

SHD batch size: 32

  * 6

SHD flt32 events integrated over 2 ms; BSNN events integrated over 3 ms.

  * 7

DO1: dropout rate in the first dense layer, DO2: dropout rate in second dense
layer

Table 7: Threshold Annealing Hyperparameters for Temporal Datasets

Dataset | $\theta_{0}^{1}$ | $\theta_{\infty}^{1}$ | $\alpha^{1}$ | $\theta_{0}^{2}$ | $\theta_{\infty}^{2}$ | $\alpha^{2}$ | $\theta_{0}^{3}$ | $\theta_{\infty}^{3}$ | $\alpha^{3}$
---|---|---|---|---|---|---|---|---|---
DVS128 Gesture | 10.4 | 12.2 | 3.3e-3 | 16.6 | 19.1 | 6.1e-3 | 6.8 | 16.5 | 0.17
SHD | 13.5 | 45.2 | 2.8e-5 | 11.2 | 51.1 | 1.36e-5 | $-$ | $-$ | $-$

  * 1

$\theta_{0}^{L}$ is the initial threshold of all neurons in layer $L$.

  * 2

$\theta_{\infty}^{L}$ is the final threshold of all neurons in layer $L$.

  * 3

$\alpha^{L}$ is the inverse time constant of the threshold for all neurons in
layer $L$.

Loss:  For both networks, the same loss function is used as in Equation 30.

##  Appendix E Discussion

###  E.1 Spike Activity on Static Datasets

(a)

(b)

Figure 7: Spiking activity across layers at the start of the training process
of BSNNs with threshold annealing (THA) and without (BSNN) on static datasets.