# Single-shot blind deconvolution with coded aperture

Hideyuki Muneta  Department of Information and Physical Sciences, Graduate
School of Information Science and Technology, Osaka University, 1-5 Yamadaoka,
Suita, Osaka 565-0871, Japan  Ryoichi Horisaki  Department of Information
Physics and Computing, Graduate School of Information Science and Technology,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
Corresponding author<EMAIL_ADDRESS>Yohei Nishizaki  Systems
and Control Laboratory, Environmental Technology Research Division, Osaka
Research Institute of Industrial Science and Technology, 1-6-50, Morinomiya,
Joto-ku, Osaka 536-8553, Japan  Makoto Naruse  Department of Information
Physics and Computing, Graduate School of Information Science and Technology,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Jun
Tanida  Department of Information and Physical Sciences, Graduate School of
Information Science and Technology, Osaka University, 1-5 Yamadaoka, Suita,
Osaka 565-0871, Japan

###### Abstract

In this paper, we present a method for single-shot blind deconvolution
incorporating a coded aperture (CA). In this method, we utilize the CA,
inserted on the pupil plane, as support constraints in blind deconvolution.
Not only an object but also a point spread function of turbulence are
estimated from a single captured image by a reconstruction algorithm with the
CA support. The proposed method is demonstrated by a simulation and an
experiment in which point sources are recovered under severe turbulence.

††journal: ao

##  1 Introduction

Blind deconvolution is a category of deconvolution methods and has been
studied for optical imaging through turbulence for astronomy, biomedicine,
security, and so on [1, 2, 3, 4]. In non-blind deconvolution, an object image
is retrieved from a known point spread function (PSF) and the captured image.
On the other hand, blind deconvolution estimates both an object image and a
PSF from a single or multiple captured images. Blind deconvolution is in high
demand in a wide range of applications, as mentioned above, compared to non-
blind deconvolution because PSFs may be unknown in many practical situations.
However, blind deconvolution is an ill-posed problem, especially in the case
of single-shot observations, and it is difficult to solve stably. Therefore,
conventionally, blind deconvolution requires various assumptions to be valid,
such as a limited size of PSFs (i.e. weak aberrations) and certain prior
knowledge about the imaging conditions.

Coded apertures (CAs) are typically used in computational imaging and have
been applied to both non-blind and blind deconvolutions [5, 6, 7]. In the case
of non-blind deconvolution, CAs have been employed for improving the condition
of inverse defocus problems in imaging systems [8, 9, 10, 11, 12]. Blind
deconvolution based on CAs has also been reported with multi-shot modalities
for imaging through turbulence, where a single CA or multiple CAs are assumed.
Aperture masking interferometry in astronomy is a single-CA-based approach for
compensating for atmospheric turbulence [13, 14, 15, 16, 17, 18]. This method
observes a large number of short-exposure images of a stationary object
through different scattering processes (i.e. dynamically changing turbulence)
by using a single CA for estimating the object’s spectrum or the object
itself. In aperture masking interferometry, the pinholes on the CA must be
aligned as non-redundant configurations. This means that the pinholes on the
CA are very sparsely arranged; hence the light efficiency is low in this
method. The multiple-CAs-based approach also has been used for blind
deconvolution, where the object is observed through a stationary scattering
process with different CAs sequentially [19, 20]. As mentioned so far, a
serious issue of these blind deconvolution methods in the literature is the
limited capability in imaging dynamically changing objects due to the multi-
shot CA modalities.

In this paper, we present a method for single-shot blind deconvolution using a
CA with high light efficiency. We consider the CA as a support of the pupil
function during the blind deconvolution process. By using the CA support, the
number of estimated variables on the pupil plane is reduced by partially
blocking the light from the object, and the condition of the inverse problem
is improved. This method may allow us to solve the single-shot blind
deconvolution problem stably by using not only computations but also optical
hardware. Our study will contribute to imaging through turbulence in various
applications by enhancing imaging speeds, tolerating strong turbulence, and
simplifying optical hardware.

##  2 Method

Figure 1: Schematic diagram of single-shot blind deconvolution with a coded
aperture (CA).

A diagram of the proposed method is shown in Fig. 1. In our method, an object
$f$ illuminated with spatially incoherent light is captured as a single
intensity image $g$ through turbulence and the CA. The CA is located on the
pupil plane in the imaging optics and is described as a binary random pattern
$A$. By assuming shift-invariant aberration of the turbulence, the forward
model is written with a convolution process as

$\displaystyle\begin{split}&g=f\otimes h_{\rm i},\\\ &h_{\rm
i}=|\mathcal{F}^{-1}[P\times A]|^{2},\end{split}$ (1)

where $h_{\rm i}$ is the point spread function (PSF) of incoherent light
through the turbulence and the imaging optics, and $P$ is the transfer
function of the aberrations by the turbulence under coherent light, which we
call the aberration pupil function. These variables are two-dimensional.
Lowercase variables are ones in the spatial domain, and uppercase variables
are ones in the Fourier domain, respectively. Here, $\otimes$ denotes the
convolution operator, $\times$ denotes the element-wise product operator, and
$\mathcal{F}^{-1}$ is the inverse Fourier transform, respectively.

The inverse problem of the single-shot blind deconvolution with the CA in Eq.
(1) is written as

$\displaystyle\widehat{f}=\operatorname*{arg\,min}_{\widehat{f},\widehat{P}(\\{u|A(u)=1\\})}\|g-\widehat{f}\otimes|\mathcal{F}^{-1}[\widehat{P}\times
A]|^{2}\|_{2}^{2},$ (2)

where $u$ is a set of frequency coordinates, $\|\bullet\|_{2}$ is the
$\ell_{2}$-norm, and the accent symbol on $\widehat{\quad}$ denotes estimated
variables. Instead of estimating the PSF $h$ in the conventional blind
deconvolution, we retrieve the aberration pupil function $P$ on the binary CA
support $A$ simultaneously with the reconstruction of the object $f$. In this
case, the estimated variables on the aberration pupil function $P$ are only at
the frequency coordinates where the light passes through the CA ($A(u)=1$),
and variables at the frequency coordinates where the light is blocked by the
CA ($A(u)=0$) do not need to be estimated. Therefore, we can limit the number
of estimated variables by using the support $A$ even when the PSF $h$ is large
due to strong aberrations.

Figure 2: Flow of the reconstruction algorithm. $e_{g}$: the error of the
captured image. $\widehat{F}$: the estimated Fourier spectrum of
$\widehat{f}$. $\widehat{H_{\rm i}}$: the estimated incoherent transfer
function. $\widehat{h_{\rm c}}$: the estimated coherent PSF. $e_{h_{\rm c}}$:
the error of the coherent PSF. IFT: the inverse Fourier transform.

We solve the inverse problem in Eq. (2) based on the ptychographic iterative
engine (PIE), which has been used for coherent diffractive imaging [21, 22].
The flow of the reconstruction is shown in Fig. 2. We initially set the
estimated original image $\widehat{f}$ and the estimated aberration pupil
function $\widehat{P}$, and calculate the simulated captured image
$\widehat{g}$ by using the forward model of Eq. (1). Then, the error $e_{g}$
between the real and simulated captured images is computed as

$\displaystyle e_{g}=g-\widehat{g}.$ (3)

Based on the PIE, this error separately back-propagates to the Fourier
spectrum $\widehat{F}$ of the estimated object image and the estimated
transfer function $\widehat{H_{\rm i}}$ of incoherent light, which is the
Fourier transform of $\widehat{h_{\rm i}}$, as

$\displaystyle\widehat{F}=\widehat{F}+\frac{\widehat{H_{\rm
i}}^{*}\times\mathcal{F}[e_{g}]}{(1-{\alpha}_{F})|\widehat{H_{\rm
i}}|^{2}+{\alpha}_{F}{|\widehat{H_{\rm i}}|^{2}_{\rm max}}},$ (4)
$\displaystyle\widehat{H_{\rm i}}=\widehat{H_{\rm
i}}+\frac{\widehat{F}^{*}\times\mathcal{F}[e_{g}]}{(1-{\alpha}_{H_{\rm
i}})|\widehat{F}|^{2}+{\alpha}_{H_{\rm i}}{|\widehat{F}|^{2}_{\rm max}}},$ (5)

where $\mathcal{F}$ is the Fourier transform, the superscript of ∗ denotes the
complex conjugation operator, and the subscript of $\rm max$ is the maximum
value on an image. The parameters $\alpha_{F}$ and $\alpha_{H_{\rm i}}$ are
constants for the PIE-based feedback to $\widehat{F}$ and $\widehat{H_{\rm
i}}$, respectively. The estimated object image $\widehat{f}$ is the inverse
Fourier transform of $\widehat{F}$ in Eq. (4).

Next, we estimate the aberration pupil function $\widehat{P}$, as shown in the
right half-plane of Fig. 2. The tentative PSF $\widehat{h_{\rm c}}^{\prime}$
of coherent light is described with $\widehat{P}$ and $A$ as

$\displaystyle\widehat{h_{\rm c}}^{\prime}=\mathcal{F}^{-1}[\widehat{P}\times
A].$ (6)

Its rectified PSF $\widehat{h_{\rm c}}$ is calculated with an amplitude
constraint with $\widehat{H_{\rm i}}$ in Eq. (5) as

$\displaystyle\widehat{h_{\rm c}}=\sqrt{\mathcal{F}^{-1}[\widehat{H_{\rm
i}}]}{\rm exp}(i\times{\rm arg}(\widehat{h_{\rm c}}^{\prime})),$ (7)

where $i$ is the imaginary number, and ${\rm arg}(\bullet)$ is the argument of
complex variables. The error $e_{h_{\rm c}}$ between $\widehat{h_{\rm c}}$ and
$\widehat{h_{\rm c}}^{\prime}$ is written as

$\displaystyle e_{h_{\rm c}}=\widehat{h_{\rm c}}-\widehat{h_{\rm
c}}^{\prime},$ (8)

and it back-propagates to the tentative aberration pupil function
$\widehat{P}^{\prime}$ as

$\displaystyle\widehat{P}^{\prime}=\widehat{P}+\frac{A^{*}\times\mathcal{F}[e_{h_{\rm
c}}]}{(1-{\alpha}_{P})|A|^{2}+{\alpha}_{P}{|A|^{2}_{\rm max}}},$ (9)

where the parameter $\alpha_{P}$ is a constant for the PIE-based feedback to
$\widehat{P}^{\prime}$. The rectified aberration pupil function $\widehat{P}$
of $\widehat{P}^{\prime}$ is calculated with an amplitude constraint by
assuming turbulence with no light absorption (or spatially invariant light
absorption) as

$\displaystyle\widehat{P}={\rm exp}(i\times{\rm arg}(\widehat{P}^{\prime})).$
(10)

This process is iterated until $\|e_{g}\|_{2}^{2}$ converges.

##  3 Demonstration

Figure 3: Simulation results. 3 The object image and 3 the CA. 3-3 The
aberration phase map, 3-3 the captured image, and 3-3 the reconstructed image.
333 Simulation results with eight Zernike coefficients ($N=8$). 333 Simulation
results with sixteen Zernike coefficients ($N=16$). 333 Simulation results
with twenty-four Zernike coefficients ($N=24$).

We demonstrated the proposed method in simulation and optical experiments. In
the simulation, we performed image reconstruction under aberrations with
different strengths, as shown in Fig. 3. We generated the phase map $w$ of the
aberrations by using Zernike coefficients and calculated the aberration pupil
function $P$ as

$\displaystyle w(r,\theta)=\sum_{n=1}^{N}a_{n}z_{n}(r,\theta),$ (11)
$\displaystyle P={\rm exp}(i\times w),$ (12)

where $r$ is the radial distance, $\theta$ is the azimuthal angle, $a_{n}$ is
the Zernike coefficients, $z_{n}$ denotes the Zernike series, and $N$ is the
polynomial order assumed in the simulation [23]. The order $N$ was set to 8,
16, and 24 to change the strength of the aberration as in the phase maps of
Figs. 3-3, where the coefficients $a_{n}$ were randomly generated between $-1$
and 1. In the simulation, the image size was set to $128\times 128$ pixels.
The object image in Fig. 3 was composed of sparsely arranged point sources
with different sizes, assuming astronomical observations. The CA in Fig. 3 was
a binary random pattern. Its transmittance ratio, which was defined as the
ratio of the number of white pixels to the total pixel count on the CA, was 50
%. The captured images in Figs. 3-3 were calculated by using the forward model
in Eq. (1), where white Gaussian noise with a signal-to-noise ratio (SNR) of
30 dB was added. The images in Figs. 3-3 were reconstructed from these single
captured images with 100,000 iterations of the loop in Fig. 2. The ambiguity
of image translations in blind deconvolution was compensated for based on the
cross-correlation between the object and reconstructed images. The peak SNRs
(PSNRs) between the original and reconstructed images in Figs. 3-3 were 36.6
dB, 42.9 dB, and 36.1 dB, respectively. Therefore, in this simulation, the
proposed method was successfully demonstrated with various strengths of
aberrations.

Figure 4: Experimental setup. SLM: spatial light modulator.

Next, we experimentally demonstrated the proposed method with the optical
setup shown in Fig. 4. The setup consisted of two 4f systems, where the
upstream one was for introducing the aberrations and the downstream one
implemented imaging through the CA. The focal length of lenses for the 4f
systems was 100 mm. The object was composed of multiple holes on a flat
aluminum foil, as shown in Fig. 5, and it was illuminated with spatially
incoherent light from a green light emitting diode (LED: LST1-01F06-GRN1-00
manufactured by New Energy, wavelength: 530 nm). Light from the object was
disturbed through the upstream 4f system, where a lens with a focal length of
200 mm was located on the pupil plane to serve as turbid media in Fig. 4 for
the shift-invariant aberration. The aberrated image was captured through the
downstream 4f system with the CA implemented by a spatial light modulator
(SLM: LC2012 manufactured by HOLOEYE, pixel count: $1024\times 768$, pixel
pitch: 36 µm) in the amplitude mode and a monochrome image sensor (PL-D7512
manufactured by PixeLink, pixel count: $4096\times 3000$, pixel pitch: 3.45
µm).

Figure 5: Experimental results. 5 The direct image of the object. 5-5 The
captured images with different transmittance ratios. 5-5 The reconstructed
images from Figs. 5-5. The transmittance ratios are shown in the upper parts
on each of the images. The length of the scale bar on Fig. 5 is 0.5 mm.

In the experiment, binary random patterns with different transmittance ratios
were used for the CA. As shown in Fig. 5, the transmittance ratio was changed
from 10 % to 100 % with intervals of 10 %. The exposure time was set to 0.13
seconds at all the transmittance ratios, and it was chosen for the
transmittance ratio of 100 % (the brightest case) to prevent light saturation.
The captured images with each transmittance ratio are shown in Figs. 5-5,
where the point sources on the object plane were not visually recognizable.
These captured images were the results obtained after downsampling with a
factor of three and cropping the central area of $200\times 200$ pixels. The
reconstruction results from the captured images are shown in Figs. 5-5. By
comparisons between them and Fig. 5, the object was successfully reconstructed
based on the proposed method with the transmittance ratios between 30 % and 60
%. Possible reasons for the failure cases may be noisy measurements due to
lower light intensities when the transmittance ratios were less than 30 %, and
insufficient support constraints with denser CAs when the transmittance ratios
were larger than 60 %, respectively.

##  4 Conclusion

In summary, we proposed and demonstrated a single-shot imaging method for
blind deconvolution with a CA. Single-shot blind deconvolution is an ill-posed
inverse problem. In the proposed method, the CA located on the pupil plane in
the imaging optics works as the support constraint in the inverse problem. An
object is captured with the CA-based imaging optics through turbulence, and it
is reconstructed from the single captured image by using the PIE-based
algorithm incorporating the CA-based support constraint. The proposed method
was numerically and experimentally demonstrated by reconstructing point
sources from images captured through severe aberrations. Our method realized
single-shot blind deconvolution by using compact and low-cost imaging hardware
with a high light efficiency.

The reconstruction algorithm in this study assumed shift-invariant and
adsorption-free turbulence. These assumptions may be resolved by modifying the
algorithm although they are acceptable in various applications. The
demonstrations here can be directly applied to astronomical observations and
biomedical fluorescence microscopy, where many imaging targets are composed of
incoherent point sources [24, 25, 26, 27, 28]. Therefore, our single-shot
blind deconvolution method will contribute to various fields, including
astronomy and biology.

## References

  * [1] G. R. Ayers and J. C. Dainty, “Iterative blind deconvolution method and its applications,” Opt. Lett. 13, 547–549 (1988).
  * [2] T. J. Schulz, “Multiframe blind deconvolution of astronomical images,” J. Opt. Soc. Am. A 10, 1064–1073 (1993).
  * [3] D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Processing Magazine 13, 43–64 (1996).
  * [4] S. Chaudhuri, R. Velmurugan, and R. Rameshan, _Blind Image Deconvolution_ (Springer, 2014).
  * [5] J. G. Ables, “Fourier transform photography: A new method for X-ray astronomy,” Publications of the Astronomical Society of Australia 1, 172–173 (1968).
  * [6] R. H. Dicke, “Scatter-hole cameras for X-rays and gamma rays,” Publications of the Astronomical Society of Australia 153, L101–106 (1968).
  * [7] E. E. Fenimore and T. M. Cannon, “Coded aperture imaging with uniformly redundant arrays,” Appl. Opt. 17, 337–347 (1978).
  * [8] A. Veeraraghavan, R. Raskar, A. Agrawal, A. Mohan, and J. Tumblin, “Dappled photography: Mask enhanced cameras for heterodyned light fields and coded aperture refocusing,” ACM Trans. Graph. 26, 69 (2007).
  * [9] A. Levin, R. Fergus, F. Durand, and W. T. Freeman, “Image and depth from a conventional camera with a coded aperture,” ACM Trans. Graph. 26, 70 (2007).
  * [10] C. Zhou, S. Lin, and S. K. Nayar, “Coded aperture pairs for depth from defocus and defocus deblurring,” International Journal of Computer Vision 93, 53–72 (2011).
  * [11] M. S. Asif, A. Ayremlou, A. Sankaranarayanan, A. Veeraraghavan, and R. G. Baraniuk, “Flatcam: Thin, lensless cameras using coded aperture and computation,” IEEE Transactions on Computational Imaging 3, 384–397 (2017).
  * [12] R. Horisaki, Y. Okamoto, and J. Tanida, “Deeply coded aperture for lensless imaging,” Opt. Lett. 45, 3131–3134 (2020).
  * [13] A. Labeyrie, “Attainment of diffraction limited resolution in large telescopes by Fourier analysing speckle patterns in star images,” Astronomy and Astrophysics 6, 85 (1970).
  * [14] J. E. Baldwin, C. A. Haniff, C. D. Mackay, and P. J. Warner, “Closure phase in high-resolution optical imaging,” Nature 320, 595–597 (1986).
  * [15] C. A. Haniff, C. D. Mackay, D. J. Titterington, D. Sivia, J. E. Baldwin, and P. J. Warner, “The first images from optical aperture synthesis,” Nature 328, 694–696 (1987).
  * [16] C. A. Haniff and D. F. Buscher, “Diffraction-limited imaging with partially redundant masks. i. infrared imaging of bright objects,” J. Opt. Soc. Am. A 9, 203–218 (1992).
  * [17] X. Peng, G. J. Ruane, M. B. Quadrelli, and G. A. Swartzlander, “Randomized apertures: high resolution imaging in far field,” Opt. Express 25, 18296–18313 (2017).
  * [18] S. Sallum and J. Eisner, “Data reduction and image reconstruction techniques for non-redundant masking,” The Astrophysical Journal Supplement Series 233, 9 (2017).
  * [19] R. Horstmeyer, X. Ou, J. Chung, G. Zheng, and C. Yang, “Overlapped Fourier coding for optical aberration removal,” Opt. Express 22, 24062–24080 (2014).
  * [20] J. Chung, G. W. Martinez, K. C. Lencioni, S. R. Sadda, and C. Yang, “Computational aberration compensation by coded-aperture-based correction of aberration obtained from optical fourier coding and blur estimation,” Optica 6, 647–661 (2019).
  * [21] A. M. Maiden, M. J. Humphry, and J. M. Rodenburg, “Ptychographic transmission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A 29, 1606–1614 (2012).
  * [22] R. Horisaki, K. Fujii, and J. Tanida, “Diffusion-based single-shot diffraction tomography,” Opt. Lett. 44, 1964–1967 (2019).
  * [23] J. D. Schmidt, _Numerical simulation of optical wave propagation with examples in MATLAB_ (SPIE, 2010).
  * [24] J. Bobin, J.-L. Starck, and R. Ottensamer, “Compressed sensing in astronomy,” IEEE Journal of Selected Topics in Signal Processing 2, 718–726 (2008).
  * [25] M. Bertero, P. Boccacci, G. Desiderà, and G. Vicidomini, “Image deblurring with poisson data: from cells to galaxies,” Inverse Problems 25, 123006 (2009).
  * [26] E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M. W. Davidson, J. Lippincott-Schwartz, and H. F. Hess, “Imaging intracellular fluorescent proteins at nanometer resolution,” Science 313, 1642–1645 (2006).
  * [27] M. J. Rust, M. Bates, and X. Zhuang, “Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM),” Nature Methods 3, 793–796 (2006).
  * [28] K. F. Tehrani, J. Xu, Y. Zhang, P. Shen, and P. Kner, “Adaptive optics stochastic optical reconstruction microscopy (AO-STORM) using a genetic algorithm,” Opt. Express 23, 13677–13692 (2015).


# Co-learning synaptic delays, weights and adaptation in spiking neural
networks

Lucas Deckers, Laurens Van Damme, Ing Jyh Tsang, Werner Van Leekwijck & Steven
Latré
IDLab, University of Antwerp - imec
Sint-Pietersvliet 7
Antwerp, Belgium
<EMAIL_ADDRESS>
###### Abstract

Spiking neural networks (SNN) distinguish themselves from artificial neural
networks (ANN) because of their inherent temporal processing and spike-based
computations, enabling a power-efficient implementation in neuromorphic
hardware. In this paper, we demonstrate that data processing with spiking
neurons can be enhanced by co-learning the connection weights with two other
biologically inspired neuronal features: 1) a set of parameters describing
neuronal adaptation processes and 2) synaptic propagation delays. The former
allows the spiking neuron to learn how to specifically react to incoming
spikes based on its past. The trained adaptation parameters result in neuronal
heterogeneity, which is found in the brain and also leads to a greater variety
in available spike patterns. The latter enables to learn to explicitly
correlate patterns that are temporally distanced. Synaptic delays reflect the
time an action potential requires to travel from one neuron to another. We
show that each of the co-learned features separately leads to an improvement
over the baseline SNN and that the combination of both leads to state-of-the-
art SNN results on all speech recognition datasets investigated with a simple
2-hidden layer feed-forward network. Our SNN outperforms the ANN on the
neuromorpic datasets (Spiking Heidelberg Digits and Spiking Speech Commands),
even with fewer trainable parameters. On the 35-class Google Speech Commands
dataset, our SNN also outperforms a GRU of similar size. Our work presents
brain-inspired improvements to SNN that enable them to excel over an
equivalent ANN of similar size on tasks with rich temporal dynamics.

##  1 Introduction

Spiking neural networks (SNN), seen as the third generation of neural network
models (Maass, 1997), have recently attracted growing attention as a low-power
alternative for artificial neural networks (ANN). Unlike ANN implementations
(García-Martín et al., 2019), SNN can enable power-efficient processing on
specific neuromorphic hardware such as SENeCa (Yousefzadeh et al., 2022), the
Intel Loihi2 (Orchard et al., 2021) or IBM TrueNorth (DeBole et al., 2019).
The main power gains can be attributed to SNN inherent event-based
computations (by means of spikes) and sparsity, reducing the number of
multiplications that are required. A recent study (Stöckl & Maass, 2021)
illustrated these potential SNN energy savings.

Recent years have shown great progress in learning algorithms for deep SNN.
Especially training SNN based on backpropagation-through-time with surrogate
gradients (Neftci et al., 2019), helped overcome the problem of the non-
differentiability, introduced by thresholding mechanism in a spiking neuron.
These advances enabled SNN to move to deeper and more complicated model
architectures using attention mechanisms (Yao et al., 2023) or transformers
(Zhou et al., 2023) and (Zhu et al., 2023). The main issue related to SNN
however persists: frequently the SNN model does not perform as well as the
equivalent ANN.

In biology, researchers have found that the brain is equipped with a plethora
of powers to process spike trains adequately. One of those is the axonal
delay. The transmission speed of an action potential is known to depend on the
myelination of the axon and thus determine how a spike is delayed (Purves et
al., 2001). Furthermore, these delays are crucial in sensory processing
(Orchard & Etienne-Cummings, 2014) and known to adapt during the learning
process (Lin & Faber, 2002). A delay-enabled spiking network was also found to
be able to compute a richer class of functions than a threshold circuit with
adjustable weights (Maass & Schmitt, 1999). Another differentiator between
classical ANN and the brain is the widespread heterogeneity and neuronal
adaptation processes taking place. Neurons of all forms and shapes are found,
each enabling a wide array processing functions (Gerstner & Kistler, 2002).
Moreover, biological neurons exhibit slow dynamical processes that act at
longer time scales, enabling processing of events that are temporally
distanced in an implicit manner. These adaptive processes often limit the
number of spikes produced.

Combining the ability to optimize the weights, delays and training neuronal
parameters can lead to a more diverse and possibly improved internal
representation. Figure 1 shows the responses of A) a typical leaky integrate-
and-fire (LIF) neuron, B) a neuron with delayed input spike trains and C) a
neuron with delayed input spike trains and neuronal adaptation, for a spike
train of four equidistant input spikes and all equal connection weights. It
can clearly be seen that the output spikes patterns are vastly different. The
LIF neuron (A) will always respond the same while the others provide a wider
range of possible results because of the trainable delay (B) and non-linear
trainable adaptation processes (C). Both extensions show to be complementary
in providing further memory. The axonal delay allows the neuron to explicitly
correlate incoming spikes at longer timescales whereas the adaptation
implicitly alters a neurons behavior based on its past regime.

Figure 1: Illustration of the variety in responses for different neuron
models: Two input spike trains are processed by three different neurons with
equal weights. Input neuron 1 spikes at $(t_{2}$, $t_{4})$ and input neuron 2
spikes at $(t_{1}$, $t_{3})$. For all neurons we show the evolution of the
membrane potential over time in response to these input spike trains: (a) a
typical leaky-integrate-and-fire (LIF) neuron, spikes at the timestep of the
last incoming spike, $t_{4}$. (b) a LIF neuron with delayed input spikes from
input neuron 2, produces spikes at timesteps $t_{2}$ and $t_{4}$. (c) an
adaptive neuron processes the same delayed input spikes from input neuron 2
and only produces a spike at timestep $t_{3}+d_{23}$, the latency of the
spikes coming from neuron 2.

In this study, we present a SNN with adaptive neurons and synaptic delays that
are co-optimized. Whereas normally, as in ANN, just the synaptic weights are
trained, we show that co-learning the delays and adaptation parameters,
individually enhance the performance of the SNN model, and that the
combination of them also leads to state-of-the-art SNN results. The
contribution of this study is as follows:

  * •

We present a novel SNN model in which both synaptic weights and delays are co-
optimized in collaboration with the neuronal adaptation parameters.

  * •

We analyze the biologically plausible neuronal adaptation parameters and which
effects parameter boundaries have on the neurons working regime and on the SNN
performance on three speech recognition datasets.

  * •

We show that inclusion of these more complex neurons through adaptation and
the addition of trainable delays for every synapse specifically leads to
state-of-the-art results for spiking neural networks. The proposed SNN even
outperforms its non-spiking counterpart with equivalent model size on the
speech recognition problems.

##  2 Related work

Learning algorithms for spiking neural networks. Recently, there has been a
rapid evolution in the development and progress of SNN learning paradigms. In
general, either a trained ANN is converted into a rate-based SNN, aiming at
minimal performance losses due to the ANN-SNN conversion (Deng & Gu, 2021),
(Bu et al., 2022) or the SNN is trained directly as a spiking neural network.
A directly trained SNN is typically trained with backpropagation-through-time
with surrogate gradients (Neftci et al., 2019), (Zenke & Vogels, 2021). These
surrogates are used to approach the derivative of non-differential Heaviside
function, introduced by the spiking mechanism. A recent study went one step
further in using differentiable spikes (Li et al., 2021) for temporal credit
assignment. Moreover, the spike-element-wise (SEW) ResNet (Fang et al., 2021a)
was proposed for training SNN without the vanishing/exploding problem
introduced by surrogate gradients, paving the way for deeper SNN models.

Learning of neuronal parameters and neuronal heterogeneity. Another trend in
SNN research is to learn the optimal distribution of the neuronal (leakage)
parameters (Yin et al., 2021), (Fang et al., 2021b) and (Rathi & Roy, 2023).
Moreover, the Gated LIF neuron (Yao et al., 2022) was proposed to control the
fusion of learnable membrane-related parameters. Additionally, in many works
the heterogeneity of neurons in SNN proved to be beneficial for improved
recognition (Deckers et al., 2022), (Perez-Nieves et al., 2021) (Chakraborty &
Mukhopadhyay, 2023).

Different methods for including neuronal adaptation processes have been
proposed. A first class contains neurons with an adaptive threshold, which is
increased after every spike and exponentially decays over time (Salaj et al.,
2021), (Yin et al., 2021). Others (Falez et al., 2019) proposed a method for
tuning the thresholds with specific target spike timestamps as an objective.
In other work (Gast et al., 2020), the adaptive thresholds were combined with
a synaptic depression model, which showed to replicate both bursting and
steady-state behavior. Another adaptation method is based on adaptation
currents, coupling a secondary variable to the subthreshold membrane potential
and its spike activity (Brunel et al., 2003). This method was successfully
formalized into the AdLIF spiking neuron model (Bittar & Garner, 2022) for SNN
and shown to outperform adaptive threshold based models.

Delays in SNN. Many methods have been proposed for adapting propagation
delays, inspired by spike timing dependent plasticity (Wang et al., 2013) or
based on the ReSuMe learning rule (Zhang et al., 2020). A method for training
per neuron axonal delays based on the SLAYER learning paradigm was proposed
(Shrestha & Orchard, 2018) and extended (Sun et al., 2023) with trainable
delay caps. Recently, the effects of axonal synaptic delay learning were
studied by pruning multiple delay synapses (Patiño-Saucedo et al., 2023),
modeling a one-layer multinomial logistic regression with synaptic delays
(Grimaldi & Perrinet, 2023) and learning delays represented trough 1D
convolutions with learnable spacings (Hammouamri et al., 2023). Similarly, in
order to train synaptic delays, spike trains were transformed into continuous
analog, differentiable signals (Wang et al., 2019). Surprisingly, only
learning the delays (Grappolini & Subramoney, 2023) showed to achieve
comparable performance as only learning the weights.

##  3 Synaptic delay-learning in adaptive spiking neural networks

A spiking neural network (SNN) is a biologically inspired type of neural
network, in which spikes, i.e., binary events are used to communicate between
layers of spiking neurons. In this section we elaborate on the fundamental
properties of spiking neurons and the methods used in this work to train the
synaptic weights, synaptic delays and neuronal parameters in multi-layer
networks of spiking neurons.

###  3.1 Spiking neurons

Spiking neurons differ from classical neurons in ANN because of their inherent
time-dependent processing of input data. Incoming spikes are multiplied by the
synaptic weights and accumulated over time for every neuron. When this
neuronal state, i.e., the membrane potential crosses the spiking threshold, a
neuron emits a spike to a subsequent layer. The membrane potential is
maintained over time and thus creates an internal memory for every individual
neuron.

In this paper, we use the adaptive leaky integrate-and-fire neuron (AdLIF)
model (Bittar & Garner, 2022) with updated parameter boundaries. To highlight
this modification, the neuron model is named AdLIF+. In this model, two
internal states are kept: the membrane potential and the adaptation current,
which provides the neuronal adaptation. Formally, $u[t]$, $w[t]$, and $s[t]$
respectively represent the membrane potential, the adaptation current and the
presence of a spike, at time step t. In the AdLIF neuron model, there are four
trainable neuronal parameters: $\alpha$ and $\beta$ denote the leak of u[t]
and w[t] respectively, while a and b describe the characteristics of the
adaptation current. The adaptation current is coupled with the subthreshold
membrane potential by $a$ while $b$ represents the spike-triggered adaptation.
A neuron generates a spike at timestep t when the membrane potential crosses
the firing threshold, $\theta$, at t. Mathematically, the threshold $\theta$
is represented as a Heaviside-function. The full discrete-time model, with
time step size equal to 1ms, is described in Equation (1). The trainable
neuronal parameters are highlighted in bold.

$\begin{split}&u[t]=\boldsymbol{\alpha}u[t-1]+(1-\boldsymbol{\alpha})(I[t]-w[t-1])-\theta
s[t-1]\\\
&w[t]=\boldsymbol{\beta}w[t-1]+(1-\boldsymbol{\beta})\boldsymbol{\textit{a}}u[t-1]+\boldsymbol{\textit{b}}s[t-1]\\\
&s[t]=u[t]\geq\theta\end{split}$ (1)

The adaptation, implemented by means of this adaptation current $w[t]$ is
affected by both the neurons instantaneous membrane potential and a spike-
triggered fraction. This contrasts with the classic adaptive neuronal
threshold adaptation, in which only the spiking activity is taken into account
(Yin et al., 2021). This model was chosen because of its proven superior
performance in comparison with a leaky integrate-and-fire (LIF) model with an
adaptive neuronal threshold (Bittar & Garner, 2022). The computational graph
of the neuron model, rolled out over time, is shown in Figure 2. The yellow
box, $w[t]$, represents the addition of the adaptation variable to the classic
LIF neuron model and the blue arrows connecting the the internal variables
represent the corresponding trainable neuron parameters.

###  3.2 Training multi-layer SNN

####  3.2.1 Model architecture

Similarly to earlier studies, the SNN in this work consists of a simple feed-
forward network with two hidden layers. Figure 2 shows the network
architecture, rolled out over time. In general, neuron $i$ in hidden layer $l$
receives $I_{i}^{l}$, the pre-synaptic current, which consists of two
elements, as shown in Equation (2). The feed-forward synapses from the
previous layer $l-1$ have associated weights $F_{ji}^{l-1}$ and carry spikes
from the same time step $t$ to neuron j and the neuronal bias, $b_{i}^{l}$.
These input spike trains are summed over all pre-synaptic neurons
$j=1,...,N_{l-1}$. Spike trains are represented by $s[t]\in[0,1]$. Typically,
in this type of architecture all feed-forward are connected in an all-to-all
fashion.

$I_{l}^{i}=\sum_{j=1}^{N_{l-1}}F_{ji}^{l-1}s_{j}^{l-1}[t]+b_{i}^{l}$ (2)

The readout mechanism, which is used to derive the outputs of the SNN,
consists of a single layer. This layer consists of neurons with infinite
threshold. These neurons have no memory and hence the membrane potential is
equal to the the weighted inputs. The output of the SNN model is the sum of
the membrane potentials of the output neurons over time, passed though a
softmax layer. The outputs are shown in Equation (3).

$u_{out}=\sum_{t}\frac{e^{u_{out}[t]}}{\sum_{j}e^{u_{out,j}[t]}}$ (3)

Figure 2: Left: Computational graph of the AdLIF neuron model, unrolled over
time. The yellow blocks denote the additional adaptation current parameter,
which depends on the membrane potential u[t] and the spike activity, s[t] at
the previous timestep. Right: 2-layer fully connected architecture with
readout of the output neurons membrane potential, unrolled over the time. The
blue connections denote the potential delays in the network. For clarity the
delays starting from time step 1 onwards were omitted.

####  3.2.2 Training procedure

Typically, spiking neural networks (SNN) are trained via backpropagation-
through-time (BPTT) with surrogate gradients (Neftci et al., 2019). In these
methods, the summed membrane potentials of the output neurons, see Equation
(3), are used to constitute the cross-entropy loss of the network, unrolled
over time. The loss with respect to class c for a batch size N is represented
as:

$\mathcal{L}_{c}=\frac{1}{N}\sum^{N}_{n=1}-log(\frac{e^{u_{out,c}}}{\sum_{j}e^{u_{out,j}}})$
(4)

Based on the chain rule in error-backpropagation, the weight update for neuron
$i$ in the penultimate layer $l$ for a sequence of T timesteps is shown in
Equation (5).

$\frac{\delta\mathcal{L}_{c}}{\delta
w^{l}}=\frac{1}{T}\sum^{T}_{t=1}\sum^{t}_{m=0}\frac{\delta\mathcal{L}_{c}[t]}{\delta
u_{out}[m]}\frac{\delta u_{out}[m]}{\delta s_{l}[m]}\frac{\delta
s_{l}[m]}{\delta u_{l}[m]}\frac{\delta u_{l}[m]}{\delta w_{l}}$ (5)

In these methods, surrogate gradients are used to approximate the non-
differentiability, which was introduced by the thresholding mechanism
(Heaviside function), $\frac{\delta s[t]}{\delta u[t]}$, with a differentiable
function in the backward pass of error-backpropagation. For simplicity and
comparability with previous work, the boxcar surrogate gradient function,
shown in Equation (6), was used in this work.

$\frac{\delta s[t]}{\delta u[t]}=\begin{cases}0.5&\text{if $|u[t]-\theta|\leq
0.5$}\\\ 0&\text{otherwise}\end{cases}$ (6)

###  3.3 Introducing synaptic delays

In contrast to ANNs and most other SNN in which connections are defined by a
weight parameter, in our work connections between neurons are characterised by
two synaptic parameters: a weight and a delay. The addition of synaptic delays
leads to an adjusted neuronal processing model. Now, the spike train for
neuron $i$ from layer $l$ is characterised by $I_{l}^{i}$ in Equation (7). The
synaptic delay can be found in $s_{k}^{l}[t-d_{ji}]$, where $d_{ji}$ denotes
the delayed spikes between neuron j and i.

$I_{l}^{i}=\sum_{j=1}^{N_{l-1}}F_{ji}^{l-1}s_{j}^{l-1}[t-d_{ji}]+b_{i}^{l}$
(7)

Training of the synaptic delays is based on the SLAYER learning method
(Shrestha & Orchard, 2018) for training axonal delays, applied to individual
synapses. The delay kernel, $\epsilon_{d}$, is convolved with spike train
s[t], to get the delayed spike kernel, $a^{l}[t]=(\epsilon_{d}*s^{l})[t]$.
Similarly to the case in which delays were equal to 0 (i.e. without delays),
the derivative of the loss with respect to the synaptic weight and delays are
now computed as shown in (8) and (9).

$\frac{\delta\mathcal{L}_{c}}{\delta
w^{l}}=\frac{1}{T}\sum^{T}_{t=1}\sum^{t}_{m=0}\frac{\delta\mathcal{L}_{c}[t]}{\delta
u_{out}[m]}\frac{\delta u_{out}[m]}{\delta a_{l}[m]}\frac{\delta
a_{l}[m]}{\delta u_{l}[m]}\frac{\delta u_{l}[m]}{\delta w_{l}}$ (8)

$\frac{\delta\mathcal{L}_{c}}{\delta
d^{l}}=\frac{1}{T}\sum^{T}_{t=1}\sum_{m=1}^{t}\frac{\delta\mathcal{L}_{c}[t]}{\delta
u_{out}[m]}\frac{\delta u_{out}[m]}{\delta a_{l}[m]}\frac{\delta
a_{l}[m]}{\delta d_{l}}$ (9)

In Equation (9), $\frac{\delta a_{l}}{\delta d_{l}}$ is equal to
$\frac{da}{dt}$ for every unique synapse.

##  4 Experiments

In this section we will elaborate on the experiments conducted for this study
and the corresponding results. First we describe the datasets and the precise
setup for training the SNNs. Following this, we present an analysis of the
AdLIF neuron model and lastly show the full the results.

###  4.1 Datasets

We used three common SNN benchmark speech recognition datasets: the Spiking
Heidelberg Digits (SHD), the Spiking Speech Commands (SSC) (Cramer et al.,
2020) and the Google Speech Commands v0.02 (GSC) dataset (Warden, 2018). The
first two datasets are neuromorphic datasets, in which the original sounds
were converted to spikes, spread out over 700 input channels. The GSC dataset
consists of speech samples. The SHD dataset consists of German and English
spoken digits (0 through 9). The SSC dataset is based on the sounds from the
GSC dataset. Both the SSC and GSC datasets have 35 classes from a large group
of speakers in a non-controlled environment. These larger datasets provide a
more challenging speech recognition task.

Similarly to other works, for the spiking datasets, which were zero-padded and
aligned to 1s, we binned the input data from 700 input channels to 140
channels and 100 timesteps with 10ms bins to ensure uniformity across all
samples. Regarding the GSC dataset, speech data was aligned to 1s by padding
with zeros and thereafter binned in 10 millisecond bins to generate samples of
100 timesteps. Further processing was performed with a Mel filterbank with 40
Mel filters. Since there is no predefined test set for the SHD dataset, we
decided to take both a maximal and an average accuracy on the validation set
across 10 experiments with different random seeds.

###  4.2 Training setup

All neuronal trainable parameters were uniformly initialized between specific
boundaries and subsequently co-learned with the other model trainable
parameters to reflect the neuronal heterogeneity (Perez-Nieves et al., 2021).
These parameters were initialized following a uniform distribution:
$\alpha\in[0.36,0.96]$, $\beta\in[0.96,0.99]$, $a\in[0,1]$ and $b\in[0,2]$.
Different from the AdLIF model (Bittar & Garner, 2022), we extended the
available range of the membrane potential decay parameter $\alpha$ and limited
the dependency of the adaptation current with respect to the membrane
potential, $a$, to remain positive, which showed to stabilize the neuron model
for sparse input data. As indicated earlier, we refer to our variant of this
neuron model as AdLIF+. During the training process, all neuron parameters are
clipped to remain within these boundaries. The spike threshold was fixed at 1
(dimensionless). The weights of all connections were initialised following the
default Xavier uniform distribution and the neuronal biases were set to zero.
For all hidden neurons, the initial membrane potential and adaptation current
were randomly initialized, following a uniform distribution between 0 and 1.

We used the Adam optimizer (Kingma & Ba, 2014) in all experiments with an
initial learning rate of 0.01 for the SHD dataset and 0.001 for the others,
similarly to previous work (Bittar & Garner, 2022). The learning rate for the
delays was always equal to $10*lr_{weigths}$ We used a simple scheduler for
both the weights and delays, which decreased the learning rate by a factor of
0.7 with a patience of 5 epochs. Dropout (Srivastava et al., 2014) was applied
in the hidden layers with rates of 0.5 and 0.25 for SHD and the other datasets
respectively. The available delay values are limited to remain within [0, 25]
time steps for these datasets. The network is initialised without delays i.e.
all delays are equal to zero and the input data was right padded accordingly
for the network to allow delayed spikes. This exact setup was used for all
experiments on all datasets presented in this work and implemented in PyTorch
(Paszke et al., 2019). Our code is based on the Sparch implementation (Bittar
& Garner, 2022).

###  4.3 Analysis of the trainable adaptation parameters

Figure 3: Illustration of different neuron model adaptation parametrizations
and their responses to a fixed train of 12 incoming spikes. We analyse the
number of output spikes and the corresponding evolution of the membrane
potential u[t] over time for a range of $a$ and $b$ neuronal adaptation
parameters. A model without adaptation, the LIF model, is found at where both
$a$ and $b$ are equal to zero. The AdLIF model, which allows both positive and
negative $a$, can possibly result in an unstable spiking regime. In this case
the neuron generates more spikes than the number of incoming spikes, 12 in
this example. We therefore limit the updated neuron model, AdLIF+ to remain
within positive $a$ and $b$ boundaries. For reference, we also show the
parametrizations for 6 and 1 output spikes in purple and green respectively.

In this section we analyse the effects of the different parameter ranges for
the main adaptation parameters $a$ and $b$ from the AdLIF neuron model, as
presented in Equation 1. Our AdLIF+ model differs from the AdLIF model in two
ways: we extended the available decay window for the membrane potential, the
$\alpha$ parameter, for the model to be able to forget more quickly if needed.
More importantly, we limited the available range of the $a$ parameter, which
allows the current membrane potential u[t] to influence the adaptation current
w[t]. Similarly to adaptation by means of an adaptive threshold, this $a$
parameter mainly limits spikes if the membrane potential was high in the past
and thus provides a homeostatic mechanism, given that this $a$ remains
positive.

In Figure 3, we illustrate that for a possible negative $a$, the behavior of
the adaptation current can lead to a non-controlled, chaotic spiking regime.
In this figure we count the number of spikes that are produced by a neuron,
given a single fixed input spike train of 12 spikes, for different adaptation
parametrizations of $a$ and $b$. Whenever the produced number of spikes is
higher than the number of input spikes, the model enters a chaotic regime in
which a positive feedback loop could be activated and the neuron remains
spiking, even without the presence of input spikes. Figure 3, shows that is
mainly the case for a negative $a$, where apart from a region with a very high
parameter $b$, the spike-triggered fraction of the adaptation current, the
number of output spikes generated is very high. We therefore, unlike other
works, preferred to remain in the top right quadrant with both positive $a$
and $b$, and hence the AdLIf+ name for our neuron model. In this figure, the
base LIF model can be found at the point where $a$ and $b$ are equal to 0 and
no adaptation current w[t] is produced.

###  4.4 Results

In the first experiment we investigate the effects of adding the trainable
parameters to a basic LIF neuron model with similar initialization, adding the
trainable synaptic delays and the interplay of co-learning both as well as the
effect of updated adaptation parametrization boundaries. The results in terms
of classification accuracy and highlighting the individual contributions are
shown in Table 1, as tested on the SHD dataset. Due to its small size, the SHD
dataset is more suitable for conducting exploratory experimemts. When
comparing the LIF model with the AdLIF and the AdLIF+ model, independent of
the inclusion of learning the delay value, we see that the average validation
accuracy as well as its maximal value over 10 experiments is significantly
increased by up to 7.2% and 9,6% respectively. The number of additional
trainable parameters is very limited as it is only increased by about 2.1%.
These results clearly show the added value of trainable adaptation parameters
for SNN models. The AdLIF+ model performs about 2.5% better than the AdLIF
model on average, which is in accordance with the empirical of adaptation
parameters $a$ and $b$ in the analysis from section 4.3. These results show
the importance of adequate trainable neuron parameter boundaries.

Table 1: Effect of the proposed modifications to a basic LIF-based 2 hidden
layer feedforward SNN model with a hidden layer size of 128. Given the limited
size of the SHD validation set, the (average, max) accuracy over 10 runs are
shown.

| Basic LIF | AdLIF |  AdLIF+ | Synaptic delay LIF |  Synaptic delay AdLIF+
---|---|---|---|---|---
# Parameters | 37.9k | 38.7k | 38.7k | 74.8k | 75.8k
Accuracy (%) | (84.49, 86.22) | (91.67, 93.85) | (94.19, 95.14) | (92.18, 93.02) | (95.02, 96,26)

The addition of trainable synaptic delays further increases the accuracy for
both the LIF and AdLIF+ neuron models. This shows that the extra provided
capacity to utilize memory is beneficial for the speech recognition task,
complementing the computational complexity provided by the trained adaptation.
Logically, the number of trainable parameters is doubled as every synapse is
now characterised by both a delay and a weight value.

An overview of the results of our experiments on the full ADLIF+ model with
synaptic delays across all speech recognition datasets is presented in Table
2. We compared our SNN model to state-of-the-art SNN solutions from literature
with a 2-hidden layer architecture and a corresponding state-of-the-art ANN
model for each dataset, shown below the dotted line.

Table 2: Test Accuracy on SHD, SSC and GSC datasets and comparison with state-
of-the-art in SNN for 2 hidden layer feedforward models.

Dataset | Model | Hidden size | # Parameters | Accuracy
---|---|---|---|---
| Adaptive RSNN (Yin et al., 2021) | 128 | / | 90.4%
SHD | RadLIF (Bittar & Garner, 2022) | 1024 | 3.9M | 94.62%
| Axonal delays (Sun et al., 2022) | 128 | 0.1M | 92.36%
| Synaptic delays (Hammouamri et al., 2023) | 256 | 0.2M | 95.07±0.24%
| AdLIF+, synaptic delays (ours) | 128 | 0.076M | 96.26%
| CNN (Cramer et al., 2020) | / | / | 92.4%
| Adaptive RSNN (Yin et al., 2021) | 400 | / | 74.2%
SSC | RadLIF (Bittar & Garner, 2022) | 1024 | 3.9M | 77.4%
| 1D-Conv synaptic delays (Hammouamri et al., 2023) | 512 | 0.7M | 79.77 ± 0.09%
| AdLIF+, synaptic delays (ours) | 512 | 0.71M | 79.81%
| GRU (Bittar & Garner, 2022) | 512 | / | 79.05%
| RSNN, LIF (Zenke & Vogels, 2021) | 256 | / | 85.3%
GSC | RSNN with SFA (Salaj et al., 2021) | 2048* | / | 88.5%
| RadLIF (Bittar & Garner, 2022) | 512 | 0.83M | 94.51%
| Synaptic delays (Hammouamri et al., 2023) | 512 | 0.7M | 94.91 ±0.09%
| AdLIF+, synaptic delays (ours) | 512 | 0.61M | 95.38%
| GRU (Bittar & Garner, 2022) | 512 | / | 94.32%
| Transformer (Gong et al., 2021) | / | / | 98.11%

*Recurrent SNN with a single hidden layer.

For the SHD dataset, the AdLIF+ model with trained synaptic delays matches the
state-of-the-art results of a recently proposed alternative method for
training synaptic delays at just a fraction (less than half) of its number of
trainable parameters and outperforms all other SNN and ANN methods on this
dataset. On the harder SSC and GSC datasets, the AdLIF+ model with trainable
synaptic delays shows on par or better performance than the current state-of-
the-art models in SNN with a similar number of trainable parameters or less.
Furthermore, one additional SNN model with 3 hidden layers was proposed
Hammouamri et al. (2023). This model slightly outperforms ours on the SSC
dataset 80.29 ± 0.06% with an increase of about 45% additional trainable
parameters, but only achieved 95.29 ± 0.11% on the GSC dataset, which is
similar to our SNN with just 2 hidden layers. Given the budget of trainable
parameters, the proposed feed-forward SNN model even outperforms a non-SNN
recurrent model (GRU) with the same preprocessing and moves closer to the
performance of a large ANN model with transformer architecture.

##  5 Conclusion

In this paper, we presented a novel SNN model, the AdLIF+ with trainable
synaptic delays. To the best of our knowledge, this is the first SNN model in
which the synaptic delays are directly learned in coordination with the
neuronal adaptation. Furthermore, we showed that 1) it is possible to co-learn
synaptic weights, delays and neuronal adaptation parameters at the same time,
and 2) co-learning these parameters proved to mutually benefit the
optimization of all learned parameters as shown for three speech recognition
datasets. We highlighted the importance of co-learning the synaptic delays and
neuronal adaptation parameters with the synaptic weights and showed that they
mutually benefit from this co-optimization, reaching state-of-the-art
performance in SNN on all investigated datasets. The superior performance can
be attributed to two additional features: 1) Training the synaptic delays
enables a neuron in the SNN to explicitly correlate temporally distanced
features 2) the trained neuronal adaptation allows a greater variety in spike
patterns, widening the feature space to be explored.

We showed that for a very simple architecture, a 2 hidden layer fully
connected feedforward network, we are able to compete against or outperform
larger ANN models, with a limited number of additional trainable SNN
parameters. The performance of the presented SNN model shows the promise of
SNN research on tasks with rich temporal dynamics, and in particular research
on biologically inspired extensions to existing SNN models.

When comparing to larger ANN models, the performance of the presented SNN
model is lacking. A future step in our research is therefore to investigate
how learning delays and adaptation parameters are influenced by the model
architecture. More advanced architectures such as convolutional spiking neural
networks or experimenting with the training recurrent synapses could further
bridge the gap with ANNs. Another exciting avenue to be explored is the chosen
neuron model. As the proposed AdLIF+ model is still one particular Generalized
integrate-and-fire model version Gerstner & Kistler (2002) and Gerstner et al.
(2014) and many more exist, there are plenty available neuron model extensions
to be investigated. Future research will point out to what extent these will
be useful for application in SNN in the context of 1) their additional
performance in terms of classification accuracy and 2) the additional
complexity for their deployment on dedicated neuromorphic hardware.

Another point to explore is that our work includes a fixed maximal delay,
which needs to be defined before training and requires fine-tuning. In future
work we intend to investigate the effect of co-learning the synaptic weights
and delays on more complex neuron models and model architectures as well as
validate them on non-sound datasets.

## References

  * Bittar & Garner (2022) Alexandre Bittar and Philip N Garner.  A surrogate gradient spiking baseline for speech command recognition.  _Frontiers in Neuroscience_ , 16, 2022.
  * Brunel et al. (2003) Nicolas Brunel, Vincent Hakim, and Magnus JE Richardson.  Firing-rate resonance in a generalized integrate-and-fire neuron with subthreshold resonance.  _Physical Review E_ , 67(5):051916, 2003.
  * Bu et al. (2022) Tong Bu, Wei Fang, Jianhao Ding, PENGLIN DAI, Zhaofei Yu, and Tiejun Huang.  Optimal ANN-SNN conversion for high-accuracy and ultra-low-latency spiking neural networks.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=7B3IJMM1k_M.
  * Chakraborty & Mukhopadhyay (2023) Biswadeep Chakraborty and Saibal Mukhopadhyay.  Heterogeneous recurrent spiking neural network for spatio-temporal classification.  _Frontiers in Neuroscience_ , 17:994517, 2023.
  * Cramer et al. (2020) 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.
  * DeBole et al. (2019) Michael V DeBole, Brian Taba, Arnon Amir, Filipp Akopyan, Alexander Andreopoulos, William P Risk, Jeff Kusnitz, Carlos Ortega Otero, Tapan K Nayak, Rathinakumar Appuswamy, et al.  Truenorth: Accelerating from zero to 64 million neurons in 10 years.  _Computer_ , 52(5):20–29, 2019.
  * Deckers et al. (2022) Lucas Deckers, Ing Jyh Tsang, Werner Van Leekwijck, and Steven Latré.  Extended liquid state machines for speech recognition.  _Frontiers in neuroscience.-Lausanne_ , 16:1–14, 2022.
  * Deng & Gu (2021) Shikuang Deng and Shi Gu.  Optimal conversion of conventional artificial neural networks to spiking neural networks.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=FZ1oTwcXchK.
  * Falez et al. (2019) Pierre Falez, Pierre Tirilly, Ioan Marius Bilasco, Philippe Devienne, and Pierre Boulet.  Multi-layered spiking neural network with target timestamp threshold adaptation and STDP.  In _2019 International Joint Conference on Neural Networks (IJCNN)_ , pp. 1–8. IEEE, 2019.
  * Fang et al. (2021a) Wei Fang, Zhaofei Yu, Yanqi Chen, Tiejun Huang, Timothée Masquelier, and Yonghong Tian.  Deep residual learning in spiking neural networks.  _Advances in Neural Information Processing Systems_ , 34:21056–21069, 2021a.
  * Fang et al. (2021b) 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_ , pp. 2661–2671, 2021b.
  * García-Martín et al. (2019) Eva García-Martín, Crefeda Faviola Rodrigues, Graham Riley, and Håkan Grahn.  Estimation of energy consumption in machine learning.  _Journal of Parallel and Distributed Computing_ , 134:75–88, 2019.
  * Gast et al. (2020) Richard Gast, Helmut Schmidt, and Thomas R Knösche.  A mean-field description of bursting dynamics in spiking neural networks with short-term adaptation.  _Neural computation_ , 32(9):1615–1634, 2020\.
  * Gerstner & Kistler (2002) Wulfram Gerstner and Werner M Kistler.  _Spiking neuron models: Single neurons, populations, plasticity_.  Cambridge university press, 2002.
  * Gerstner et al. (2014) Wulfram Gerstner, Werner M Kistler, Richard Naud, and Liam Paninski.  _Neuronal dynamics: From single neurons to networks and models of cognition_.  Cambridge University Press, 2014.
  * Gong et al. (2021) Yuan Gong, Yu-An Chung, and James Glass.  Ast: Audio spectrogram transformer.  _arXiv preprint arXiv:2104.01778_ , 2021.
  * Grappolini & Subramoney (2023) Edoardo W Grappolini and Anand Subramoney.  Beyond weights: Deep learning in spiking neural networks with pure synaptic-delay training.  _arXiv preprint arXiv:2306.06237_ , 2023.
  * Grimaldi & Perrinet (2023) Antoine Grimaldi and Laurent U Perrinet.  Learning heterogeneous delays in a layer of spiking neurons for fast motion detection.  _arXiv preprint arXiv:2307.14077_ , 2023.
  * Hammouamri et al. (2023) Ilyass Hammouamri, Ismail Khalfaoui-Hassani, and Timothée Masquelier.  Learning delays in spiking neural networks using dilated convolutions with learnable spacings.  _arXiv preprint arXiv:2306.17670_ , 2023.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Li et al. (2021) Yuhang Li, Yufei Guo, Shanghang Zhang, Shikuang Deng, Yongqing Hai, and Shi Gu.  Differentiable spike: Rethinking gradient-descent for training spiking neural networks.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), _Advances in Neural Information Processing Systems_ , volume 34, pp. 23426–23439. Curran Associates, Inc., 2021.  URL https://proceedings.neurips.cc/paper_files/paper/2021/file/c4ca4238a0b923820dcc509a6f75849b-Paper.pdf.
  * Lin & Faber (2002) Jen-Wei Lin and Donald S Faber.  Modulation of synaptic delay during synaptic plasticity.  _Trends in neurosciences_ , 25(9):449–455, 2002\.
  * Maass (1997) Wolfgang Maass.  Networks of spiking neurons: the third generation of neural network models.  _Neural networks_ , 10(9):1659–1671, 1997.
  * Maass & Schmitt (1999) Wolfgang Maass and Michael Schmitt.  On the complexity of learning for spiking neurons with temporal coding.  _Information and Computation_ , 153(1):26–46, 1999.
  * Neftci et al. (2019) Emre O Neftci, Hesham Mostafa, and Friedemann Zenke.  Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks.  _IEEE Signal Processing Magazine_ , 36(6):51–63, 2019.
  * Orchard & Etienne-Cummings (2014) Garrick Orchard and Ralph Etienne-Cummings.  Bioinspired visual motion estimation.  _Proceedings of the IEEE_ , 102(10):1520–1536, 2014.
  * Orchard et al. (2021) Garrick Orchard, E Paxon Frady, Daniel Ben Dayan Rubin, Sophia Sanborn, Sumit Bam Shrestha, Friedrich T Sommer, and Mike Davies.  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) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Patiño-Saucedo et al. (2023) Alberto Patiño-Saucedo, Amirreza Yousefzadeh, Guangzhi Tang, Federico Corradi, Bernabé Linares-Barranco, and Manolis Sifalakis.  Empirical study on the efficiency of spiking neural networks with axonal delays, and algorithm-hardware benchmarking.  In _2023 IEEE International Symposium on Circuits and Systems (ISCAS)_ , pp. 1–5. IEEE, 2023.
  * Perez-Nieves et al. (2021) Nicolas Perez-Nieves, Vincent CH Leung, Pier Luigi Dragotti, and Dan FM Goodman.  Neural heterogeneity promotes robust learning.  _Nature communications_ , 12(1):5791, 2021.
  * Purves et al. (2001) D Purves, GJ Augustine, D Fitzpatrick, LC Katz, AS LaMantia, JO McNamara, SM Williams, et al.  The organization of the nervous system.  _Neuroscience. Purves D, Augustine GJ, Fitzpatrick D, Katz LC, LaMantia AS, McNamara JO, et al, editors. Sunderland, MA: Sinauer Associates_ , 2001.
  * Rathi & Roy (2023) Nitin Rathi and Kaushik Roy.  Diet-snn: A low-latency spiking neural network with direct input encoding and leakage and threshold optimization.  _IEEE Transactions on Neural Networks and Learning Systems_ , 34(6):3174–3182, 2023.  doi: 10.1109/TNNLS.2021.3111897.
  * Salaj et al. (2021) Darjan Salaj, Anand Subramoney, Ceca Kraisnikovic, Guillaume Bellec, Robert Legenstein, and Wolfgang Maass.  Spike frequency adaptation supports network computations on temporally dispersed information.  _Elife_ , 10:e65459, 2021.
  * Shrestha & Orchard (2018) Sumit B Shrestha and Garrick Orchard.  Slayer: Spike layer error reassignment in time.  _Advances in neural information processing systems_ , 31, 2018.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ , 15(1):1929–1958, 2014.
  * Stöckl & Maass (2021) Christoph Stöckl and Wolfgang Maass.  Optimized spiking neurons can classify images with high accuracy through temporal coding with two spikes.  _Nature Machine Intelligence_ , 3(3):230–238, 2021.
  * Sun et al. (2022) Pengfei Sun, Longwei Zhu, and Dick Botteldooren.  Axonal delay as a short-term memory for feed forward deep spiking neural networks.  In _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pp. 8932–8936. IEEE, 2022.
  * Sun et al. (2023) Pengfei Sun, Ehsan Eqlimi, Yansong Chua, Paul Devos, and Dick Botteldooren.  Adaptive axonal delays in feedforward spiking neural networks for accurate spoken word recognition.  _arXiv preprint arXiv:2302.08607_ , 2023.
  * Wang et al. (2013) Runchun Wang, Gregory Cohen, Klaus M Stiefel, Tara Julia Hamilton, Jonathan Tapson, and André van Schaik.  An fpga implementation of a polychronous spiking neural network with delay adaptation.  _Frontiers in neuroscience_ , 7:14, 2013.
  * Wang et al. (2019) Xiangwen Wang, Xianghong Lin, and Xiaochao Dang.  A delay learning algorithm based on spike train kernels for spiking neurons.  _Frontiers in neuroscience_ , 13:252, 2019.
  * Warden (2018) Pete Warden.  Speech commands: A dataset for limited-vocabulary speech recognition.  _arXiv preprint arXiv:1804.03209_ , 2018.
  * Yao et al. (2023) Man Yao, Guangshe Zhao, Hengyu Zhang, Yifan Hu, Lei Deng, Yonghong Tian, Bo Xu, and Guoqi Li.  Attention spiking neural networks.  _IEEE transactions on pattern analysis and machine intelligence_ , 2023.
  * Yao et al. (2022) Xingting Yao, Fanrong Li, Zitao Mo, and Jian Cheng.  Glif: A unified gated leaky integrate-and-fire neuron for spiking neural networks.  _Advances in Neural Information Processing Systems_ , 35:32160–32171, 2022.
  * Yin et al. (2021) Bojian Yin, Federico Corradi, and Sander M Bohté.  Accurate and efficient time-domain classification with adaptive spiking recurrent neural networks.  _Nature Machine Intelligence_ , 3(10):905–913, 2021.
  * Yousefzadeh et al. (2022) Amirreza Yousefzadeh, Gert-Jan Van Schaik, Mohammad Tahghighi, Paul Detterer, Stefano Traferro, Martijn Hijdra, Jan Stuijt, Federico Corradi, Manolis Sifalakis, and Mario Konijnenburg.  Seneca: Scalable energy-efficient neuromorphic computer architecture.  In _2022 IEEE 4th International Conference on Artificial Intelligence Circuits and Systems (AICAS)_ , pp. 371–374. IEEE, 2022.
  * Zenke & Vogels (2021) Friedemann Zenke and Tim P Vogels.  The remarkable robustness of surrogate gradient learning for instilling complex function in spiking neural networks.  _Neural computation_ , 33(4):899–925, 2021.
  * Zhang et al. (2020) Malu Zhang, Jibin Wu, Ammar Belatreche, Zihan Pan, Xiurui Xie, Yansong Chua, Guoqi Li, Hong Qu, and Haizhou Li.  Supervised learning in spiking neural networks with synaptic delay-weight plasticity.  _Neurocomputing_ , 409:103–118, 2020.
  * Zhou et al. (2023) Zhaokun Zhou, Yuesheng Zhu, Chao He, Yaowei Wang, Shuicheng YAN, Yonghong Tian, and Li Yuan.  Spikformer: When spiking neural network meets transformer.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=frE4fUwz_h.
  * Zhu et al. (2023) Rui-Jie Zhu, Qihang Zhao, Guoqi Li, and Jason K. Eshraghian.  SpikeGPT: Generative pre-trained language model with spiking neural networks, 2023.

##  Appendix A Appendix

###  A.1 Training hyperparameters

Table 3 shows an overview of the exact parametetrization of the experiments,
executed for this study. The setups for the SSC and GSC datasets were
identical.

Table 3: Overview of the hyperparameters used in training for all datasets
considered in this study

Dataset | Hidden layer size | Epochs | Batch size | Dropout rate | Lr weights | Lr delays | Initialization | Delay caps
---|---|---|---|---|---|---|---|---
SHD | 128 | 100 | 128 | 0.5 | 0.01 | 0.1 | Xavier uniform | [0,25]
SSC & GSC | 512 | 200 | 32 | 0.25 | 0.001 | 0.01 | Xavier uniform | [0,25]

###  A.2 Number of spikes per neuron per method

One of the properties that affects the efficiency of the SNN model, when
deployed on neuromorphic hardware is the effective number of spikes that are
required when processing a single sample. We therefore analyse the number of
spikes in the hidden layers to assess how the proposed improvements influence
the total number of spikes in the SNN. To account for inter-experimental
variance, we averaged the spike rates over 10 runs on the SHD dataset. The
results are shown in in Table 4.

We see that the number of spikes is not significantly increased by the
addition of any of the considered adaptation methods. The synaptic delays
however result in an increase of the number of spikes in the SNN model.

Model | Basic LIF | AdLIF |  AdLIF+ | Synaptic delay LIF |  Synaptic delay AdLIF+
---|---|---|---|---|---
Number of spikes/neuron | 10.75 | 10.7 | 10.28 | 13.7 | 11.48

Table 4: The average number of spikes per neuron for a single sample, compared
over different SNN models. These results were averaged over 10 runs on the SHD
dataset.

###  A.3 Model architecture

Figure 4 shows the network architecture. The input spikes are passed trough 2
hidden layers. In the readout layer the membrane potential is extracted from
memoryless output neurons. The recurrent connections represent the membrane
potential leakage and the adaptation currents.

Figure 4: Network architecture: a 2-hidden layer feedforward SNN model

###  A.4 Detailed results on SHD dataset

In this section we show the results of the 10 different runs with different
pseudo-random seeds. Figure 5 provides a more general overview than the
average and maximal results, as provided in Table 1. These results again show
the importance of both the neuronal adaptation and the trainable synaptic
delays in the SNN model recognition performance.

Figure 5: Detailed overview of the results on the SHD dataset for all SNN
model setups.

###  A.5 Visualisation: trained neuronal parameters

In this section we visualize the distribution of the trained parameters for
every layer in the fully connected SNN, trained on the SHD dataset. This
particular model achieved an accuracy of 95.45% on the validation set with the
displayed parameter distribution.

Figure 6: Distributions of the trained parameters ($\alpha$, $\beta$, a and b)
in the first hidden layer Figure 7: Distributions of the trained parameters
($\alpha$, $\beta$, a and b) in the second hidden layer

We observe that although all neuronal parameters were initialized uniformly
between their respective boundaries: $a\in$ [0, 1] and $b\in$ [0, 2], the
optimized neuronal parameters follow a different distribution. Mainly
parameter $b$ shows a remarkable results. All other distributions are spread
relatively evenly between their boundaries. We see that apart from some
neurons, most neurons converged to a $b$ value of zero. This suggests that the
parameter linking the adaptation current with the sub-threshold membrane
potential, $a$, is predominantly used for adaptation. This is in accordance
with our empirical results, shown in Figure 3, where mainly a influences the
total number of output spikes for a given input spike train.

###  A.6 Visualisation: trained delays

In our experiments, the delays were initialized at zero, as if there was no
delay in the SNN. In Figure 8, we show the final delay distribution of every
layer in the SNN for the same checkpoint as shown in section (A.5) for the
neuronal parameters (validation accuracy at 95.45%). For the input-Hidden1 and
Hidden1-Hidden2 layers, we see that the optimized delays are centered around 0
and 1 while in the output layer this spread is wider. We hypothesize that in
the first 2 layers, short-term dependencies are prioritized over longer term
ones.

Figure 8: Distributions of the trained delays in the 3 trained layers of the
SNN model

###  A.7 On spiking neuron models

Over the past few years, neuroscientists have proposed many neuron models
Gerstner & Kistler (2002). On one side, the Hodgkin-Huxley model,
characterised by four differential equations, can accurately reproduce
electrophysiological measurements from biological neurons. The biological
accuracy comes at the cost of intrinsic complexity, limiting their use in
simulations of large networks of neurons. On the other side, the Leaky
integrate-and-fire (LIF) neuron, characterised by one differential equation,
was proposed for efficient processing and capturing the core properties of
spiking neurons. Nevertheles, it is known that LIF neurons are not sufficient
to describe the diversity of responses to input patterns found in neurons.

The AdLIF+ neuron, characterised by two differential equations, proposed in
this work, finds the balance between both ends. As shown in Gerstner et al.
(2014) the AdLIF+ model is able to reflect a wide range firing patterns. In
the biophysical interpretation, the sub-threshold adaptation parameter $a$ can
be linked to the slower dynamics of ion channels in relation to the membrane
potential. The spike-triggered adaptation parameter $b$ is linked with ion
channels only opening when the membrane potential is above the spike
threshold.

There are many more biologically inspired additions we can make to the
proposed AdLIF+ model in order to improve its processing power. Non-linear
(exponential) or quadratic (Izhikevich) neurons, multicompartment models,
refractoriness and stochastic neurons are just a few examples of possible
enhancements. Future research should point out to what extent these
modifications are also useful for applications with SNN. Potential gains in
SNN model performance should be weighted with likely, additional complications
in the deployment of complex neuron models on dedicated neuromorphic hardware.