Spyker: High-performance Library for Spiking Deep Neural Networks
Spiking neural networks (SNNs) have been recently brought to light due to
their promising capabilities. SNNs simulate the brain with higher biological
plausibility compared to previous generations of neural networks. Learning with
fewer samples and consuming less power are among the key features of these
networks. However, the theoretical advantages of SNNs have not been seen in
practice due to the slowness of simulation tools and the impracticality of the
proposed network structures. In this work, we implement a high-performance
library named Spyker using C++/CUDA from scratch that outperforms its
predecessor. Several SNNs are implemented in this work with different learning
rules (spike-timing-dependent plasticity and reinforcement learning) using
Spyker that achieve significantly better runtimes, to prove the practicality of
the library in the simulation of large-scale networks. To our knowledge, no
such tools have been developed to simulate large-scale spiking neural networks
with high performance using a modular structure. Furthermore, a comparison of
the represented stimuli extracted from Spyker to recorded electrophysiology
data is performed to demonstrate the applicability of SNNs in describing the
underlying neural mechanisms of the brain functions. The aim of this library is
to take a significant step toward uncovering the true potential of the brain
computations using SNNs.
Shahriar Rezghi Shirsavar, Mohammad-Reza A. Dehaqani
# Spyker: High-performance Library for Spiking Deep Neural Networks

###### Abstract

Spiking neural networks (SNNs) have been recently brought to light due to their promising capabilities. SNNs simulate the brain with higher biological plausibility compared to previous generations of neural networks. Learning with fewer samples and consuming less power are among the key features of these networks. However, the theoretical advantages of SNNs have not been seen in practice due to the slowness of simulation tools and the impracticality of the proposed network structures. In this work, we implement a high-performance library named Spyker using C++/CUDA from scratch that outperforms its predecessor. Several SNNs are implemented in this work with different learning rules (spike-timing-dependent plasticity and reinforcement learning) using Spyker that achieve significantly better runtimes, to prove the practicality of the library in the simulation of large-scale networks. To our knowledge, no such tools have been developed to simulate large-scale spiking neural networks with high performance using a modular structure. Furthermore, a comparison of the represented stimuli extracted from Spyker to recorded electrophysiology data is performed to demonstrate the applicability of SNNs in describing the underlying neural mechanisms of the brain functions. The aim of this library is to take a significant step toward uncovering the true potential of the brain computations using SNNs.

 Spiking Neural Network, Learning Rules, C++/CUDA, Modular Structure, Biological Plausibility

## I Introduction

The human brain can operate with amazing robustness and energy efficiency. Artificial neural networks (ANNs) aim at modeling the brain, and three generations of these networks have been developed. Each generation of ANNs improves the quality of the modeling of the brain compared to the last. The first generation of ANNs makes use of the McCulloch-Pitts neurons [1]. Although these neurons are inspired by biological neurons, time dynamics are not considered in this model, and the learning rules proposed for them lack power and biological plausibility. These neurons were used in multi-layer perceptron (MLPs) [2] and Hopfield [3] networks.

The second generation of ANNs uses a continuous activation function (ReLU [4] and sigmoid [5], for example) instead of thresholding, which makes them suitable for processing analog signals. They have attracted the attention of researchers in recent years and were able to reach high accuracies [6, 7] (even surpassing humans) and win different challenges [8]. Despite the success of DNNs, there are structural differences between these networks and the human brain. Lack of temporal dynamics, using analog signals for network propagation and activation functions, learning rules without biological roots, and the need for large amounts of data [9] and energy [10] to achieve acceptable results are among these differences.

The third generation of neural networks is spiking neural networks (SNNs). The neural models used in these networks simulate biological neurons more accurately, and the coding mechanisms used in these networks are found in neural communications. Furthermore, the learning rules used in these networks have been discovered in the brain [11, 12, 13]. Having lower energy consumption, learning with fewer samples, and solving more complicated tasks due to time dynamics (several electrophysiological studies emphasize the role of temporal dynamics in neural coding [14, 15]) are some of the advantages of SNNs compared to the second generation of ANNs. SNNs can be used to solve machine learning tasks, study and explore brain functionality, and run on specialized hardware with low power consumption. The research being done on these networks aims to address the disadvantages of DNNs with more realistic modeling of the brain functionality.

Several high-performance well-established frameworks like PyTorch [16], TensorFlow [17], and MXNet [18] have been developed for DNNs in recent years. These libraries have enabled DNNs to achieve new highs in solving machine learning tasks. SNNs are not yet comparable to DNNs due to the lack of fast simulation tools. There have been some attempts, like SpykeTorch [19] and BindsNet [20]. SpykeTorch, written on top of the PyTorch framework, is a simulator for large-scale spiking neural networks (SDNNs). However, it has a slow runtime, and training even simple networks can take up to days to complete. To our knowledge, Spyker is the first toolbox to simulate large-scale networks with high performance, is easy to use, has the flexibility to be used in multiple languages, and has the compatibility to integrate with other commonly used tools. In order to fill this need, we have developed Spyker. Spyker is a C++/CUDA library written from scratch with both C++ and Python interfaces and support for dense and sparse structures. Although Spyker is a stand-alone library, it has a highly flexible API and can work with PyTorch tensors and Numpy arrays. Figure 1 shows an overview of the library. In order to increase performance, small-sized integers are used alongside floating-point numbers. It also uses highly-optimized low-level back-end libraries such as OneDNN and cuDNN to speed up heavy computations such as convolutions and matrix multiplications. Spyker can be compiled on various CPUs to be
optimized locally and take advantage of native CPU-specific instructions.

Spiking neural networks are made of different building blocks (see [21] for more details). The first block is the modeling of the biological neurons. Some examples of this are leaky integrate-and-fire [22], spike-response model [23], and Izhikevich model [24]. Another building block is neural coding, which can be rate coding [25], temporal coding, phase coding and synchrony coding [26], or other coding schemes. The final building block is the learning mechanism. Examples of these mechanisms are STDP [27, 28], R-STDP [29], backpropagation [30], and conversion from ANNs to SNNs [31]. Spyker has a modular implementation of these three blocks that enables its users to build SNNs.

Spyker provides SNN functionality with a high-performance and easy-to-use interface with an open-source and permissive license. It can run on CPU and CUDA devices and has both dense and sparse interfaces. The library introduces new features and fixes most of the shortcomings of its predecessor. The improvements include adding batch processing, strided convolutions, internal padding for convolutions, fully connected layers, and the rate coding mechanism. Compared to its predecessor, the interface of the library is simpler, closer to the current API of deep learning libraries, and more straightforward to use. In this work, several successful network structures are implemented using this library to prove its operability, its runtime speed is compared to SpykeTorch, and the results indicate Spyker can run up to eight times faster. The proposed work is able to reduce the gap between SNNs and DNNs and bring us a step closer to uncovering the true potential of spiking neural networks.

We start with a description of dimensionality of the input arrays and how the spike trains are implemented in the library. Afterward, we provide an explanation of different building blocks of SNNs and how they are implemented in Spyker and modeled in the interface. Then, we implement network structures that have been succesful to prove its operatibility, and we compare the performance of the library to its predecessor on these networks. Furthermore, comparison of the represented stimuli extracted from Spyker to recorded electrophisiology data is performed to demonstrate the applicability of SNNs in describing the underlying neural mechanisms of the brain functions. Finally, we demonstrate an example usage of the library and discuss the impacts of this work and how it can be further improved.

## II Methods

The interface of the Spyker can be better explained when the classes and methods of the interface are grouped by building blocks of SNNs. The categories are feature enhancement, neural coding, neural model, and learning. In this section, the structure of the input to the network is explained. Afterward, the sparse and the dense interfaces are compared. Finally, the building blocks of the library are discussed in detail.

### _Network Input_

Arrays passed through convolutional neural networks that process images are often four-dimensional arrays composed of batch size (B or N), number of channels (C), image height (H), and image width (W). The order can either be BCHW or BHWC (or NCHW or NHWC). SNNs have temporal dynamics, and it is implemented as a dimension that represents time steps in Spyker. The library implements five-dimensional arrays with BTCHW order (T being the time steps). Since DNNs process analog signals, data types used in these networks are (usually four-byte) floating-point numbers. This data type can be computationally expensive compared to a small-sized integer type and take up more space in the memory. Since SNNs process binary signals, Spyker can optionally use eight-bit (or wider) integers alongside floating-point numbers to improve performance further.

### _Dense vs Sparse interface_

The dense interface of Spyker uses the fully allocated memory buffers that are used in neural network computations. However, the sparse interface only needs to hold the indices of the spikes. Conversion between dense and sparse interfaces are provided in the library. The sparse interface has some advantages compared to the dense interface. In the dense interface, the time consumed by each operation is a function of the size of each of the 5 dimensions. However, in the sparse interface, it depends on the number of spikes. This means both memory and time consumed will be greatly reduced when processing sparser signals. Furthermore, since neurons fire at most once when using rank order coding, the increment of the number of time steps will have a smaller effect in the sparse interface compared to the dense interface.

### _Feature Enhancement_

A transformation can be used to enhance features of the input signal (image) before the neural coding process [32, 33, 34]. This results in highlighted features having higher intensities and appearing in earlier time steps, meaning more excitation. Feature enhancement is done through filtering the input here. Various filters are supported in Spyker, and they are introduced in the following subsections.

#### Ii-C1 Difference of Gaussian Filter

The first filter is the Difference of Gaussian (DoG). This filter increases the intensities of edges and other details in the image (see Figure 2 for an example) [35]. It approximates the center-surround properties of the ganglion cells of the retina [36] (see also [37, 38]). This operation is implemented as spyker.DoG(size, filters, pad, device) where size is the size of the width and the height of the filter, filters is a list of DoG filter descriptions (each description takes in two standard deviations), pad is the size of the padding of the image, and device is the device the filter will run on (CPU, GPU or others).

#### Ii-C2 Gabor Filter

The following filter is the Gabor filter that determines the presence of specific frequency in content in a specific direction in the image. Research Indicates [39] that the Gabor filter is used in the human visual cortex. The Gabor filter is implemented as spyker.Gabor(size, filters, pad, device). The parameters of this class are the same as the DoG class, but the filters are Gabor filter descriptions, and each description takes in sigma, theta, gamma, lambda, and psi.
#### Ii-C3 Laplacian of Gaussian Filter

The Laplacian of Gaussian (LoG) layer is also implemented in Spyker, and it is approximated using two DoG filters. An LoG filter with standard deviation \(\sigma\) can be approximated using two DoG filters with (\(\sigma\sqrt{2}\), \(\sigma/\sqrt{2}\)) and (\(\sigma/\sqrt{2}\), \(\sigma\sqrt{2}\)) standard deviations. This filter exists in Spyker as spyker.LoG(size, stds, pad, device) where stds are a list of standard deviations needed to describe multiple LoG filters.

#### Ii-C4 Shape of the Filters

The previously explained filters have kernel size \(K_{c}\times K_{h}\times K_{w}\), which are square kernels (\(K_{h}=K_{w}\)). The input can have \(B\times C_{i}\times H_{i}\times W_{i}\) shape which corresponds to batch, channels, height, and width of the input, respectively. The output will have \(B\times C_{o}\times H_{o}\times W_{o}\) shape where:

\[\begin{split} C_{o}&=C_{i}\times K_{c}\\ H_{o}&=H_{i}+2\times P_{h}-K_{h}+1\\ W_{o}&=W_{i}+2\times P_{h}-K_{w}+1\end{split} \tag{1}\]

and \(P_{h}\) and \(P_{w}\) are height and width padding of the filter. The \(K_{c}\) filters are applied to each channel separately.

#### Ii-C5 Zero-phase Component Analysis

Final implemented layer is zero-phase component analysis (ZCA) Whitening. It has been suggested [34] that this transformation can improve the accuracy of SNNs on real-world images. Spyker implements an efficient version of ZCA whitening by taking advantage of routines from highly optimized linear algebra libraries (BLAS and LAPACK) that operate on symmetric matrices. This layer is implemented as spyker.ZCA class which has a fit(array, epsilon) and a call function.

### _Neural Coding_

SNNs process spike trains, but the input consists of analog values (for example, images are made of pixel values). In order to make these inputs suitable for the network, a conversion scheme is needed. The mapping from stimuli to neural responses is called neural coding [40]. Coding schemes implemented in Spyker are explained in the following subsections.

#### Ii-D1 Rate Coding

Out of several coding schemes suggested, rate coding is widely used where the rate of firing of the neurons represents information. In this scheme, the rate of firing is dependent on the intensity of the input value (higher intensity corresponds to faster firing) [25]. The exact time of firing in each neuron is stochastic in nature and may be modeled with a Poisson distribution. A lengthy window of time is required to transmit the information in this coding, and the spikes are not quite sparse.

#### Ii-D2 Temporal Coding

Another popular coding scheme is temporal coding [41]. Recordings in the primary visual cortex show [42] that the response latency decreases with the stimulus contrast. This coding scheme can convey information through the timings of the spikes. Multiple forms of this scheme have been proposed, including rank order coding [43]. Instead of computing the exact timing of each spike, the timings are computed relative to one another in rank order coding. This relative (instead of exact) timing can increase invariance to changes in the input intensity and contrast [43]. It has been suggested [44] that temporal coding might be more efficient in some situations.

#### Ii-D3 Coding in Spyker

Spyker supports rank order and rate coding. The concept of time is implemented with spikes occuring in time steps in this library. Rank order coding maps higher intensities to earlier time steps of a neuron firing. In order to calculate the time step the neuron will fire in, Spyker sorts the intensity values by default. This calculates rank order between spikes, and the spikes will be distributed among time steps evenly. The sorting operation is computationally

Fig. 1: Overview of the Spyker library. Spyker API supports PyTorch tensors and Numpy arrays as well as a built-in data wrapper. The output of Spyker operations have the same container type as the input. The functionality of Spyker can be grouped into subcategories shown in the figure.
expensive (specially on GPUs), and optionally, it can be disabled to have runtime improvements (however, accuracy might be affected). Since processing time steps sequentially is inefficient and time-consuming, Spyker processes all the time steps at once. To this end, when a neuron fires in time step \(t_{i}\), it will also fire at time steps \(t_{i+1}\), \(t_{i+2}\),..., \(t_{n}\) where \(n\) is the number of time steps. An example of this cumulative structure can be seen in Figure 2.

### _Neural Model_

Once the input is filtered and coded, it gets processed by the network. The network is built using fully connected, convolution, integrate-and-fire (IF) activation, pooling, and padding layers. These operations are explained in the following subsections.

#### Iii-E1 Convolution

The integrate-and-fire mechanism is implemented by combining convolution and the IF activation layer. The internal potentials of the neurons are computed using convolution operation, and the IF activation operation produces spikes where neurons have a potential higher than a specified threshold. Multiple layers can be assembled and stacked on top of one another to create deeper structures.

The convolution layer has a kernel with \(C_{o}\times C_{i}\times K_{h}\times K_{w}\) shape. the synaptic weights are initialized randomly with a normal distribution. It performs two-dimensional convolution with support for padding and stride. The input has \(B\times T\times C_{i}\times H_{i}\times W_{i}\) shape which corresponds to batch, time steps, channels, height, and width of the input, respectively. The output has \(B\times T\times C_{o}\times H_{o}\times W_{o}\) shape where:

\[\begin{split} H_{o}&=\lfloor\frac{H_{i}+2\times P _{h}-K_{h}}{S_{h}}\rfloor+1\\ W_{o}&=\lfloor\frac{W_{i}+2\times P_{w}-K_{w}}{S_ {w}}\rfloor+1\end{split} \tag{2}\]

And \(P_{h}\), \(P_{w}\), \(S_{h}\), \(S_{w}\) are the height and width of convolution padding and stride. Padding increases the size of the two-dimensional input before convolution operation by expanding the edges of the input and filling in the new space with a constant value (usually zero). Stride is the number of steps the convolution window takes when it moves on the image. The output of the convolution layers are internal potentials of neurons that need to be passed through an IF activation layer to become output spike trains. This layer is implemented with spyker.Conv(insize, outsize, kernel, stride, pad, mean, std, device) class in Spyker.

#### Iii-E2 Fully Connected

The fully connected layer is combined with the IF activation to model the IF neurons, much similar to what happens in the convolution layers. This layer has a kernel with \(I\times O\) shape. The synaptic weights are initialized

Fig. 2: The figure shows a black and white image being filtered by DoG and Gabor filters. The theta parameter of the Gabor filter is set to -15 degrees. Then the images are coded using rank order coding into four time steps. Spikes are shown with white color on a black background through time steps. Spikes carry on from the previous to the current time step (cumulative structure).
randomly with a normal distribution. The input has \(B\times T\times I\) which corresponds to batch, time steps, and input size, respectively. The output has \(B\times T\times O\) shape. The fully connected layer is represented by spyker.FC(insize, outsize, mean, std, device) in the library.

#### Ii-D3 Pooling

The pooling layer performs two-dimensional max pooling operation with a window size of\(L_{h}\times L_{w}\), a stride of \(S_{h}\times S_{w}\), and a padding of \(P_{h}\), \(P_{w}\). The input has \(B\times T\times C_{i}\times H_{i}\times W_{i}\) shape and the output has \(B\times T\times C_{o}\times H_{o}\times W_{o}\) shape where:

\[\begin{split} H_{o}&=\lfloor\frac{H_{i}+2\times P _{h}-L_{h}}{S_{h}}\rfloor+1\\ W_{o}&=\lfloor\frac{W_{i}+2\times P_{w}-L_{w}}{S_ {w}}\rfloor+1\end{split} \tag{3}\]

The interface of Spyker has the spyker.pool(array, kernel, stride, pad, rates) function to run the pooling operation on the input given the kernel, stride, and padding size. rates argument is the rate of firing of the neurons when rate coding is used. The pooling operation selects neurons that fire earlier when rank order coding is used, and selects neurons that have a higher firing rate when rate coding is used.

### _Learning_

Learning in the brain happens when the strength of connections change between its neurons, and this change in strength is named synaptic plasticity [45]. Learning methods that utilize synaptic plasticity have been developed for SNNs [27, 28, 29].

#### Ii-F1 Spike-timing-dependent Plasticity

One widely recognized synaptic plasticity learning rule is spike-timing-dependent plasticity (STDP) [27, 28]. STDP learning rule operates by adjusting synaptic weights and utilizing the timing of the spikes. A pre-synaptic neuron firing before (after) the postsynaptic neuron results in a strengthened (weakened) connection. STDP allows the neurons to extract and learn frequent features in the input [46]. STDP layer changes synaptic weights with stabilization:

\[\Delta W_{i,j}=\begin{cases}A_{k}^{+}(W_{i,j}-L_{k})(U_{k}-W_{i,j}),&t_{j}\leq t _{i}\\ A_{k}^{-}(W_{i,j}-L_{k})(U_{k}-W_{i,j}),&t_{j}\geq t_{i}\end{cases} \tag{4}\]

where \(A_{k}^{+}\), \(A_{k}^{-}\), \(L_{k}\), \(U_{k}\) are the positive learning rate, negative learning rate, lower bound, and upper bound of the \(k\)th configuration, respectively. If stabilization is not set, then the formula becomes:

\[\Delta W_{i,j}=\begin{cases}A_{k}^{+},&t_{j}\leq t_{i}\\ A_{k}^{-},&t_{j}\geq t_{i}\end{cases} \tag{5}\]

then the weights are computed:

\[W_{i,j}^{+}=max(L_{k},min(U_{k},\Delta W_{i,j})) \tag{6}\]

Input, neurons selected by the winner-take-all mechanism (WTA), and the output are passed to a function belonging to fully connected or convolution layers, and the STDP learning rule is applied. Convolution or fully connected layers in Spyker can have multiple STDP configurations (different learning rules, weight clipping, enabling/disabling stabilizer) implemented as spyker.STDPConfig(positive, negative, stabilize, lower, upper). Each winner neuron can be mapped to an STDP configuration, and that neuron will be updated using the learning rates and such that belongs to the selected configuration. SpykeTorch creates an STDP object for each configuration, and mapping winner neurons to different configurations is done by the user. Compared to SpykeTorch, Spyker provides a more flexible and easy to use API for weight updating and enables batch updating, which improves performance. Samples are processed in mini-batches which increases performance drastically (see the results section), and the batch update rule does not differ from single-sample processing.

#### Ii-F2 Reward-modulated STDP

Another approach is using the reinforcement (RL) learning rule. One method based on RL is reward-modulated STDP [29]. R-STDP adjusts the STDP such that neurons that respond correctly are rewarded, and punished otherwise. It has been suggested [33] that when the input has non-diagnostic frequent features that are less effective in decision-making, R-STDP is able to discard these features and improve the decision-making process. Since convolution and fully connected layers accept STDP configurations as input, R-STDP can be implemented by passing two configurations to a layer (one for rewarding and one for punishing), and mapping each winner neuron to a configuration based on data labels. If one formulates this, \(\Delta W_{i,j}\) will be:

\[\begin{cases}A_{r}^{+}(W_{i,j}-L_{r})(U_{r}-W_{i,j}),&t_{pre}<t_{post}\\ A_{r}^{-}(W_{i,j}-L_{r})(U_{r}-W_{i,j}),&t_{pre}\geq t_{post}\end{cases},\quad \text{if reward}\]

\[\begin{cases}A_{p}^{-}(W_{i,j}-L_{p})(U_{p}-W_{i,j}),&t_{pre}<t_{post}\\ A_{p}^{+}(W_{i,j}-L_{p})(U_{p}-W_{i,j}),&t_{pre}\geq t_{post}\end{cases},\quad \text{if punish} \tag{7}\]

#### Ii-F3 Winner-take-all and Lateral Inhibition

When a neuron fires at a specific location, lateral inhibition [47, 48] operation inhibits other neurons belonging to other neural maps from firing in that location. Lateral inhibition for the convolution operation can be used with spyker.inhibit(array, threshold, inplace) functions. Winner neurons that STDP weight updating will be performed on are selected by the winner-take-all [49, 50] operation. WTA selects neurons that fire earlier, and if the firing time of neurons is the same, then the one that has a higher internal potential will be selected. This operation is implemented with spyker.fcwta(array, radius, count, threshold) for fully connected and spyker.convwta(array, radius, count, threshold) for convolution operations.

## III Results

In this section, we will test the performance of the library against the SpykeTorch library. Afterward, a comparison of the represented stimuli extracted from Spyker to recorded electrophysiology data is conducted to demonstrate the applicability of SNNs in describing the underlying neural mechanisms of brain functions.

### _Library Performance_

In this section, we compare the performance of the library to SpykeTorch on two networks that classify the MNIST dataset.
#### Iv-B1 R-STDP Network

The first network is the Mozfari et al. network [33] which has three convolutional layers. The first layer is trained two times with STDP, the second layer four times with STDP, and the third layer 680 times with R-STDP on the training set while computing the test accuracy at each iteration while training the third layer. We made a small change to the structure of the network (named Alt for alternative). We removed the input padding from the last convolution layer and changed its window size to 4 and the output channels to 400. Results can be seen in Figure 3 and Table I. All the tests are performed on Inte Core i7-9700k with 64G memory and Nvidia Geforce GTX 1080 Ti with 12G memory (Ubuntu 18.04).

In order to compare the results, we test whether the two-sample mean difference confidence interval (99.9%) contains zero. The null hypothesis is having the same means, and the alternative is having different means. The test results indicate that the Spyker Python implementation is faster compared to the SpykeTorch implementation (Confidence intervals are [15477, 15859] and [72607, 80737] for Spyker and SpykeTorch respectively, showing no intersection). Furthermore, the alternative implementation is faster both in the Python implementation with [-3738, -3370] interval and the C++ implementation with [-3828, -3339] interval. As expected, the C++ interface is faster compared to the Python interface with [-1078, -520] interval. The results for the accuracy comparisons show that there are no significant differences ([96.932, 98.169] and [95.996, 97.444] for Python vs SpykeTorch implementations respectively, showing intersection, [-0.89, 0.793] for C++ vs Python, [-0.649, 0.813] for Python alternative vs Python, and [-0.763, 0.971] for C++ alternative vs C++).

#### Iv-B2 STDP Network

Subsequently, the Kheradpisheh et al. network [32] is used for comparisons. This network is made of two convolutional layers. The first layer is trained 2 times with STDP, and the second layer is trained 20 times with STDP on the training set. The output of the network is classified using the SVM classifier. The elapsed time measured consists of the time needed to train the network on the training set and make predictions for the testing set. The time to utilize SVM is not taken into account because the libraries that simulate the neural network portion are compared here. The results can be seen in in Figure 4 and Table II.

The test results indicate that the Spyker GPU implementation is faster compared to the SpykeTorch implementation (confidence interval [-2728, -2265]). Since the SpykeTorch implementation processes one sample at a time, we also implemented a single sample version on the GPU, and this implementation runs faster compared to the SpykeTorch implementation (confidence interval [-1795, -1338]). There is also an implementation using the sparse interface of the Spyker (that runs on CPU) that is faster than the SpykeTorch implementation on the GPU (confidence interval [-2586, -2120]). These results show that the Spyker implementation is faster while the accuracy is not significantly different ([-0.373, 0.511] for Spyker GPU, [-0.458, 0.603] for single-sample, and [-0.405, 0.549] for sparse implementation, all against the

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline Implementation & Time & Time (S\(\pm\)SD) & Accuracy (\%\(\pm\)SD) & Runs \\ \hline SpykeTorch GPU & 47m30s & 2,850\(\pm\)64 & 98.392\(\pm\)0.093 & 30 \\ \hline Spyker GPU Single & 21m23s & 1,283\(\pm\)66 & 98.465\(\pm\)0.095 & 30 \\ \hline Spyker GPU & 05m53s & 353\(\pm\)9 & 98.461\(\pm\)0.079 & 30 \\ \hline Spyker Sparse & 08m16s & 496\(\pm\)1 & 98.464\(\pm\)0.091 & 30 \\ \hline \end{tabular}
\end{table} TABLE II: Comparisons of the the runtime and accuracy of Spyker against SpykeTorch on the Kheradpisheh et al. network.

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline Implementation & Time & Time (S\(\pm\)SD) & Accuracy (\%\(\pm\)SD) & Runs \\ \hline SpykeTorch & 21h17m & 76,672\(\pm\)916 & 96.720\(\pm\)0.163 & 12 \\ \hline Spyker Python & 04h49m & 15668\(\pm\)52 & 97.550\(\pm\)0.169 & 30 \\ \hline Spyker Python Alt & 03h31m & 12,114\(\pm\)14 & 97.632\(\pm\)0.112 & 30 \\ \hline Spyker C++ & 03h52m & 14,869\(\pm\)50 & 97.502\(\pm\)0.157 & 30 \\ \hline \end{tabular}
\end{table} TABLE I: Comparisons of the the runtime and accuracy of Spyker against SpykeTorch on the Mozafari et al. network.

Fig. 3: Comparison plots of the runtime and accuracy of Spyker against SpykeTorch on the Mozafari et al. network. The plot on the left shows the runtime comparison of Spyker and SpykeTorch implementations. The plot on the right also compares accuracy of the two implementations. Comparisons are between SpykeTorch (ST), implementation using Spyker in Python (SP Py), alternative version using Spyker in Python (SPA Py), and their C++ counterparts (SP C++, SPA C++). The error bars are minimum and maximum values of the samples.
### SpykeTorch implementation).

### _Analyzing the Underlying Structures of the Brain_

In order to demonstrate the use case and the importance of the library in neuroscience research, a similarity analysis is done in this section to compare the biological plausibility of an SNN and a deep CNN model. The neural data needed for the analysis is recorded as spiking activity and LFP signals from Inferior Temporal (IT) cortex using a single electrode (169 sessions from two macaque monkeys, the neural data for the monkeys are pooled together) [51]. The task implemented here is a Rapid Serial Visual Presentation (RSVP). The intervals are 50ms for stimulus and 450ms interstimulus. Eighty-one greyscale images of real-world objects and Gaussian low-pass filtered and high-pass filtered variations of some are shown during the task (total 155 images). The categories of the stimuli are animal face (AF), human face (HF), animal body part(AB), human body part (HB), natual objects (N), and man-made objects (MM).

The SNN used here is structurally similar to the one introduced by Shirasvar et al. [52]. The input of the SNN is resized to 32 and passed through 3 LoG filters with stds of 0.471, 1.099, 2.042. The window sizes of the filters are 7. Then, the output is thresholded and coded into 15 time steps. The first convolution layer has 16 output channels with awindow size of 5 and a padding of 2, and the second convolution layer has 32 output channels with a window size of 3 and a padding of 1. The pooling layers have 2 and 3 window sizes, respectively. The layers are trained 20 times on the images, and the learning rates are doubled after each image until they reach 0.15. Firing times (divided by number of time steps) of the final layer is used as the network output.

The CNN network used here is a ResNet-50 with the classifier layer replaced. The network is not pretrained. The input image is resized to 256 and cropped to 224. The network is trained 15 times on the dataset with Adam optimizer and 0.0001 learning rate. using a linear SVM classifier to classify the 6 categories. the accuracies for the 6 classes are 51.569 \(\pm\) 2.240 (SD), 48.623 \(\pm\) 2.538, and 51.247 \(\pm\) 2.257 for ResNet-50, SNN, and an SVM classifier that is used on the average firing rates of the neural recordings of the monkeys between 150ms and 200m from the onset, respectively. Figure 5 Shows the results of the analysis. The average Kendall's Tau value for the interval between 125ms and 175ms shown in the figure is tested between the SNN and the ResNet. Using a Mann-Whitney U test with the alpha value of 0.001 results in a p-value of 2.028-07, which shows significant difference between the two. This indicates that the SNN has a closer structure to monkey brain.

### _Rate Coding Output_

In this section, we look at the output of an SNN that uses rate coding. The SNN network used here is the Shirasvar et al. [52]. The number of output channels in the convolutional layers are set to 25 and 50. The training is not changed in that 15 time steps are used with rank order coding. However, the inference is done with 300 time steps and rate coding. Afterward, the spike output of 40 neurons are plotted for one testing sample for each class shown in Figure 6. The figure also contains a plot of T-SNE transformed firing rates as output features and the recall score for each class for the average of 30 runs. The accuracy of the 30 runs is 95.635\(\pm\)0.171 on the testing set.

## IV Library Demonstration

In this section, a sample usage of the library is illustrated. The network used here is introduced by Shirasvar et al. [52] to classify the MNIST dataset. The network has two convolutional layers trained with the STDP learning rule. The code shown in this section is only a part of the actual implementation, with the aim of providing a simple example. For the complete implementation, please visit the GitHub repository of Spyker1.

Fig. 4: Comparison plots of the runtime and accuracy of Spyker against SpykeTorch on the Kheradpisheh et al. network. The plot on the left shows shows the runtime comparison of Spyker and SpykeTorch implementations. The plot on the right also compares accuracy of the two implementations. Comparisons are between GPU implementation using SpykeTorch (SP GPU), GPU implementation using Spyker with single-sample instead of batch processing (SP Single), GPU implementation using Spyker (SP GPU), and Sparse CPU implementation using Spyker (SP Sparse). The error bars are minimum and maximum values of the samples.
## IV Conclusion

Fig. 5: Similarity comparison of SNN and ResNet-50 to monkey neural data. The similarity measurement used here is the cosine similarity. The RDM for the monkey is computed for the 50ms interval after the onset. The RDMs are adjusted with histogram equalization. The RSA is calculated with 50ms window size and 5ms stride and 95% confidence interval. Kendall’s Tau measurement is used for the RSA analysis. The RSA is averaged in the interval between 125ms and 175ms and compared in the plot in the top right with 95% confidence interval.

Fig. 6: Raster plot of an SNN network for the MNIST test images. In this figure, 40 neurons are plotted in 300 time steps for 10 samples of the MNIST testing set, each image belonging to one class.
### _Transformation_

The transformation from the input image to the network input consists of feature enhancement and spike coding, shown in Listing 1. Here, a module named Transform is defined that performs the transformation when called. This module applies 3 LoG filters with different standard deviations to the input image with padding to keep the original width and height of the input. The output is stored in 6 channels. Each channel of this output is then coded into fifteen time steps using rank order coding.

```
classTransform: def__init__(self,device): std=[0.471,1.099,2.042] self.fit=spyker.LoG(3,std, pad=3,device=device) def__call__(self,data): data=self.fit(data) spyker.threshold(data,0.01) returnspyker.code(data,15)
```

Listing 1: Implementation of the Transform module

### _Network_

The network has two convolutional layers. Here, a module named Network is defined (shown in Listing 2) to train the neurons and make predictions. Here, the convolution layers are initialized, STDP configurations are set, and the winner selection function is wrapped with a lambda function to keep the hyperparameters in the initialization of the function of the network.

```
classNetwork: def__init__(self,device): self.thresh1,self.thresh2=16,5 self.conv1=spyker.Conv(6,100,5, pad=2,mean=.5,std=.02,device=device) self.conv2=spyker.Conv(100,200,3, pad=1,mean=.5,std=.02,device=device) config1=spyker.STDPConfig(.0004,-.0003) config2=spyker.STDPConfig(.0004,-.0003) self.conv1.stdpconfig=[config1] self.conv2.stdpconfig=[config2] self.wta1=lambdax:spyker.convwa(x,3,5) self.wa2=lambdax:spyker.convwa(x,1,8)
```

Listing 2: Implementation of the Network module

### _Learning_

Training each layer is done in a separate function shown in Listing 3. The training of the layers is done in a sequential order (one layer after another). Training of the first layer is done in the train_layer1 function with the STDP learning rule. Here, the output of the first convolution is computed, and lateral inhibition is performed on it. Then, winner neurons are selected, and STDP weight updating is performed on them. The STDP learning rates in the first layer are multiplied by 1.5 every 2000 samples, and the multiplying process stops once the positive learning rate reaches 0.15. The second layer is trained in a similar way in the train_layer2 function with the STDP learning rule.

```
deftrain_layer1(self,data): output=self.conv1(data) spyker.threshold(output,self.thresh1) spyker.inhibit(output) winners=self.wta1(output) spikes=spyker.fire(output) self.conv1.stdp(data,winners,spikes) deftrain_layer2(self,data): data=self.conv1(data) data=spyker.fire(data,self.thresh1) data=spyker.pool(data,2) output=self.conv2(data) spyker.threshold(output,self.thresh2) spyker.inhibit(output) winners=self.wta2(output) spikes=spyker.fire(output) self.conv2.stdp(data,winners,spikes)
```

Listing 3: The code for training of the network layers

After defining the network module, the process of training and classification is implemented. The training process shown in Listing 4 involves training each layer once with quantization afterward.

```
fordata,targetintrainset: network.train_layer1(transform(data)) spyker.quantize(network.conv1.kernel,0,0.5,1) fordata,targetintrainset: network.train_layer2(transform(data)) spyker.quantize(network.conv2.kernel,0,0.5,1)
```

Listing 4: The training process of the network

### _Inference_

The call operator of the network shown in Listing 5 implements the prediction procedure which processes the input spikes and produces the final network output.

```
def__call__(self,data): data=self.conv1(data) data=spyker.fire(data,self.thresh1) data=spyker.pool(data,2) data=self.conv2(data) data=spyker.fire(data,self.thresh2) data=spyker.pool(data,3) returnspyker.gather(data).flatten(1)
```

Listing 5: Inference function of the network

After training, the output features for every sample in the training set and the testing set are computed (in the gather
function). Then, an SVM classifier is trained on the training set outputs. Finally, predictions are made for the testing set outputs (shown in Listing 6).

```
xtr,ytr=gather(network,transform,train) xte,yte=gather(network,transform,test) svm=LinearSVC(C=2.4).fit(xtr,ytr) pred=svm.predict(xte) accuracy=(pred==testy.numpy()).mean()
```

Listing 6: Implementation of the dimension reduction and classification operations

## V Discussion

Our brain has amazing capabilities. It can learn and perform complicated tasks in a robust manner and with low power consumption. Artificial neural networks have been created to mimic the power of the brain processes. Deep neural networks are ANNs that have had major success in recent years. However, there are structural differences between these networks and the brain, and they encounter problems when it comes to tolerance, energy, and sample efficiency. Spiking neural networks are the next generation of artificial neural networks. SNNs are not a new concept. However, they have been brought to attention recently due to their promising characteristics. The aim of these networks is to build a better model of the brain compared to DNNs.

Several well-established simulation tools exist for DNNs. These tools have allowed DNNs to reach their great success faster and have helped them to computationally scale up. SNNs lack such high-performance simulation tools. There have been some attempts at creating such tools, but they have not been able to live up to expectations. In this work, we introduced Spyker, a high-performance library written from scratch using low-level tools to simulate spiking neural networks on both CPUs and GPUs. Despite being stand-alone, Spyker has great flexibility and the ability to integrate with other tools to create a smooth developing experience. We compared the performance of this library with SpykerTorch, a simulation tool built on the PyTorch framework. We showed that Spyker is multiple times faster compared to this library. Furthermore, to demonstrate the applicability of SNNs in describing the underlying neural mechanisms of the brain functions and the role of Spyker in this field, we compared the similarity of a spiking neural network implemented with this library with the similarity of the ResNet model to a macaque monkey brain. Finally, we illustrated an example implementation to demonstrate the easy and modern interface of the library.

Strong SNN models can be implemented using the Spyker library to solve real-world machine learning problems. Features like fast processing and having a C++ interface alongside the Python interface make this library ready for both research and production. Generalization is an important concept in machine learning and having neural networks that learn and run fast are quite desirable. SNNs have the potential to become state-of-the-art models in machine learning. Other potential use cases of the library is to study and understand how the brain processes information using simulations. In other words, this library enables us to look at neuroscience through the eyes of a brain-inspired neural network.

Although this library has been shown to be perform, there is room for more improvements. Spyker has a sparse interface that runs on the CPU. The sparse interface can be extended to also run on the GPU, and this can improve the performance even further. Furthermore, the support for a larger number of neural models, coding schemes, and learning rules can be added. This helps the library to cover a great range of SNN building blocks. When choosing a model to be deployed on embedded and neuromorphic processors, SNNs are among the top choices due to their energy efficiency. SNNs are often used in neuromorphic computing. Another direction that Spyker can take is in this direction. The computational efficiency of the sparse interface of Spyker can be further improved and made compatible with these types of processors.
