Gradient-Free Training of Spiking Neural Networks via
Low-Rank Evolution Strategies


Abstract

Spiking Neural Networks (SNNs) offer compelling energy efficiency on neuromorphic hardware, yet their training remains challenging because the discrete spike threshold is non-differentiable. Surrogate-gradient methods sidestep this by approximating the derivative, but they impose backpropagation infrastructure that is incompatible with on-chip learning. Evolution Strategies (ES) are a natural gradient-free alternative, yet their computational cost scales with the number of parameters, making them impractical for large weight matrices.

We present a method for training SNNs using Eggroll, a low-rank factorisation of ES perturbations that reduces per-generation memory from \(\mathcal{O}(mn)\) to \(\mathcal{O}(r(m{+}n))\). Combining Eggroll with a Leaky Integrate-and-Fire SNN on N-MNIST, we demonstrate that gradient-free training achieves % test accuracy while reducing per-generation wall-clock time by \(\times\) relative to full-rank ES. Our results demonstrate EGGROLL is viable for SNN training, with a clear accuracy-speed tradeoff, compatible with training on neuromorphic hardware without surrogate gradients.

1 Introduction↩︎

Energy Efficiency is a huge focus of modern AI research. A lot of research is focused on reducing the energy usage of AI architecture. Intel’s Loihi  [1] showed us that SNNs could be a potential solution on real hardware. The dominant paradaigm in current AI is the transformer architecture  [2]. Transformers are trained on GPUs and use hundreds of watts in their training process. While the artificial neural networks with transformers show remarkable performance, they use a significant amount of energy. The human brain, in comparison, is extremely efficient, only using  20 watts. SNNs attempt to capture this efficiency: neurons integrate inputs over time and fire only when a membrane potential threshold is crossed, enabling sparse, event-driven computation on dedicated neuromorphic chips such as Intel Loihi [1] and IBM TrueNorth [3].

1.0.0.1 The training problem.

Despite their efficiency during inference, SNNs are still quite difficult to train. The spiking generation function(Heavside step) has a zero gradient almost everywhere. This makes it impossible to do standard backpropagation, since it requires calculating the gradient. The common solution to this is using surrogate gradients  [4]. Surrogate gradients work by replacing the true gradient of a function with an approximation, which can be differentiated. This allows backpropagation to proceed through layers where it would otherwise be impossible. However, surrogate gradients still require autograd infrastructure, which is not compatible with on chip learning on neuromorphic hardware.

1.0.0.2 Evolution Strategies.

ES [5] is a group of loosely bio inspired training methods that are alternatives to backpropagation that do not require gradients to be calculated. In ES, weight vectors are randomly perturbed and then a fitness evaluation is done on these perturbations. However these methods are very computationally expensive, incurring \(O(P mn)\) memory for each generation, for a weight matrix of dimensions \(m × n\) and population \(P\), making it very challenging to scale on larger networks.

1.0.0.3 Our contribution.

We integrate Eggroll [6] — which replaces each full-rank perturbation with a low-rank product \(\mathbf{AB}^\top\), \(\mathbf{A} \in \mathbb{R}^{m \times r}\), \(\mathbf{B} \in \mathbb{R}^{n \times r}\) — into the SNN training loop. This reduces per-generation cost to \(\mathcal{O}(r(m{+}n))\) while preserving the gradient-free property.

Concretely, our contributions are:

  1. we characterize EGGROLL’s behavior on non-differentiable spike functions Eggroll with SNNs, enabling gradient-free training without surrogate approximations (Section 3).

  2. A rank ablation showing that accuracy remains relatively stable as \(r\) decreases, with a clear efficiency–accuracy Pareto frontier (Section 4).

  3. Comparison against vanilla ES and surrogate-gradient BPTT on N-MNIST. N-MNIST is a native neuromorphic dataset, with no rate-coding approximation used, making this a more principled evaluation than the one done on a static MNIST dataset. (Section 4).

2 Background↩︎

2.1 N-MNIST Dataset↩︎

The N-MNIST dataset  [7] or Neuromorphic MNIST is a spiking version of the popular MNIST dataset that contains static images of hand drawn numbers. It consists of the same 60,000 training and 10,000 test samples as the original MNIST dataset. It was created by mounting the ATIS sensor on a motorized pan tilt unit and moving it while it recorded MNIST images on an LCD display. The sensor outputs asynchronous events with 2 polarities: ON events (higher brightness) and OFF events (lower brightness). This is a better fit for SNNs than the regular MNIST dataset because the data is already a spike train, so no rate coding approximation is required. Due to this, the SNN processes the actual sensor output instead of the surrogate encoding.

2.2 Spiking Neural Networks and the LIF Model↩︎

We model each neuron as a Leaky Integrate-and-Fire unit. At each discrete timestep \(t\), the membrane potential \(V_t\) evolves as \[V_t = \alpha V_{t-1} + \mathbf{w}^\top \mathbf{s}_{t-1}, \label{eq:lif}\tag{1}\] where \(\alpha \in (0,1)\) is the membrane decay constant (leak factor), \(\mathbf{w}\) the synaptic weight vector, and \(\mathbf{s}_{t-1}\) the binary spike vector from the previous layer. A spike is emitted when \(V_t \geq V_\text{th}\), after which \(V_t\) resets to \(V_\text{reset} = 0\). The network output is the spike count (firing rate) over \(T\) timesteps.

2.3 Surrogate Gradient Training↩︎

Because the spike function \(\Theta(V - V_\text{th})\) is a Heaviside step, its true derivative is zero almost everywhere. Surrogate gradient methods [4] substitute a smooth proxy \(\hat{\Theta}'\) — commonly the fast sigmoid — only during the backward pass, leaving the forward pass unchanged. This recovers gradient flow at the cost of introducing an approximation bias and requiring full backpropagation-through-time (BPTT). Implementing this method requires autograd, which as mentioned earlier, is not compatible with on-chip learning, and makes it unsuitable for neuromorphic hardware.

2.4 OpenAI Evolution Strategies↩︎

[5] propose estimating the gradient of expected fitness \(J(\boldsymbol{\theta})\) as \[\nabla_{\boldsymbol{\theta}} J \approx \frac{1}{2P\sigma} \sum_{i=1}^{P} \bigl[F(\boldsymbol{\theta} + \sigma\boldsymbol{\varepsilon}_i) - F(\boldsymbol{\theta} - \sigma\boldsymbol{\varepsilon}_i)\bigr] \boldsymbol{\varepsilon}_i, \label{eq:es95grad}\tag{2}\] where \(\boldsymbol{\varepsilon}_i \sim \mathcal{N}(\mathbf{0}, \mathbf{I})\), \(P\) is the population size, and \(\sigma\) the perturbation scale. Antithetic (mirrored) sampling reduces variance with no extra sampling budget. The estimate is passed to Adam [8] as if it were a true gradient. The biggest issue with ES is that its memory cost is \(\mathcal{O}(Pmn)\), making it slow for large models.

2.5 EGGROLL: Low-Rank Evolution Strategies↩︎

[6] observe that sampling a full perturbation matrix \(\mathbf{E} \in \mathbb{R}^{m \times n}\) is wasteful: for population \(P\), the perturbation tensor requires \(\mathcal{O}(Pmn)\) memory. Eggroll instead draws \[\mathbf{E}_i = \frac{1}{\sqrt{r}} \mathbf{A}_i \mathbf{B}_i^\top, \quad \mathbf{A}_i \in \mathbb{R}^{m \times r},\; \mathbf{B}_i \in \mathbb{R}^{n \times r}, \label{eq:eggroll}\tag{3}\] normalising by \(\sqrt{r}\) so that each entry has unit variance regardless of rank, without which the variance scales by \(r\), rendering the \(\sigma\) hyperparameter obsolete. Additionally, EGGROLL reconstructs noise on demand using a counter based deterministic random number generator(RNG), so that the perturbations do not have to be stored in memory. This trick helps EGGROLL work on the billion parameter scale. The gradient estimate is then reconstructed via the \((\operatorname{diag}(\mathbf{f})\mathbf{A})^\top \mathbf{B}\) formulation (§4.2 of [6]), which never materialises individual perturbation matrices. Memory per generation drops to \(\mathcal{O}(r(m{+}n))\).

3 Method↩︎

3.1 Network Architecture↩︎

We use a two-layer LIF network: \[2312 \xrightarrow{} 64 \xrightarrow{\text{LIF}} 10 \xrightarrow{\text{LIF}} \text{output}.\] The input is of size 2312 since N-MNIST has a size of \(34 X 34\) pixels, and 2 polarities for each pixel. Each LIF layer shares a common membrane decay \(\beta\) . Initial weights are sampled from \(\mathcal{N}(0,0.3^2)\). The output is the mean spike rate over \(T\) timesteps, converted to class probabilities via softmax for fitness evaluation. Biases are perturbed with independent 1-D Gaussian factors to maintain equivalence with full-rank ES when \(r = \min(m,n)\). Unlike static-image SNNs, the forward pass feeds a different event frame to each timestep not the same image repeated T times. Event counts are clamped to \([0,1]\) to prevent a single pixel dominating membrane potential.

3.2 EGGROLL Integration↩︎

Algorithm 1 details the training loop. Instead of backpropogating through the SNN the usual way, the algorithm starts with network weights and creates many perturbed copies of the network which are subsequently run on a minibatch. The performance of each copy is measured and that score is then converted into a gradient like signal which updates the weights with Adam. The key departure from vanilla ES is the batched forward pass in Line [alg:fwd]: all \(P\) perturbed networks are evaluated simultaneously by expanding the data tensor along a population dimension, avoiding a sequential Python loop over population members. The factors \(\mathbf{A}\), \(\mathbf{B}\) are regenerated from a stored seed rather than cached, keeping GPU memory independent of \(P\). Centered rank normalization is done to make the gradient scale unaffected by absolute fitness values across batches. Antithetic sampling is also used to reduce the variance.

Figure 1: Eggroll training for an SNN

3.3 Fitness Function↩︎

We use negative cross-entropy (log-likelihood) as the fitness signal, rather than raw classification accuracy: \[F(\boldsymbol{\theta}) = -\frac{1}{B}\sum_{b=1}^{B} \log \frac{\exp(\hat{y}_{b,y_b})}{\sum_c \exp(\hat{y}_{b,c})},\] where \(\hat{y}_{b,c}\) is the spike rate for class \(c\) on example \(b\). Log-likelihood provides a smooth landscape for ES to navigate; raw accuracy (a step function over fitness values) gives no gradient signal between threshold crossings.

4 Experiments↩︎

4.1 Setup↩︎

4.1.0.1 Dataset.

N-MNIST [7]: 54,000 training / 6,000 validation / 10,000 test images, generator seed 0. Sensor size \(34 X 34\), 2 polarities. No rate coding is used, the event frames are the spiking input.

4.1.0.2 Baselines.

  • Vanilla ES: full-rank Gaussian perturbations, sequential evaluation, identical hyperparameters.

  • Surrogate-gradient BPTT: fast-sigmoid surrogate, Adam optimiser, same architecture.

4.1.0.3 Compute.

All experiments run on a single NVIDIA RTX 3070 Ti. Wall-clock times are averaged over three seeds.

4.2 Main Results↩︎

Table 1: Test accuracy and per-generation wall-clock time. Mean \(\pm\) std over three seeds.
Method Test Acc.(%) Time / gen (s) Speedup
Surrogate BPTT 94.03\(\pm\)​0.28 31.69 0.41\(\times\)
Vanilla ES 76.17\(\pm\)​6.07 13.00 1.0\(\times\)
Eggroll\(r=1\) 78.75\(\pm\)​0.58 3.41 3.81\(\times\)
Eggroll\(r=2\) 72.27\(\pm\)​12.00 5.96 2.18\(\times\)
Eggroll\(r=4\) 79.21\(\pm\)​0.64 5.82 2.23\(\times\)
Eggroll\(r=8\) 72.80\(\pm\)​6.08 5.88 2.21\(\times\)

The key takeaway from 1 is the significant speedup that EGGROLL provides to ES, for almost no tradeoff in accuracy. Even with a rank of 1, the difference in accuracy is very minute relative to the speedup provided by EGGROLL. Though Surrogate BPTT has a higher accuracy than ES, it is significantly slower.

4.3 Rank Ablation↩︎

Figure 2 plots validation accuracy and per-generation wall-clock time as functions of rank \(r\).

a
b

Figure 2: Rank ablation on N-MNIST.. a — Accuracy vs.rank., b — Wall-clock vs.rank.

The accuracy remains relatively stable across ranks, with the best performance occurring at \(r=1\) and \(r=4\), with larger variance at \(r=2\) and \(r=8\). Wall clock time is lowest at \(r=1\). It increases for \(r=2\), and remains relatively stable across the other 2 ranks.

4.4 Convergence Curves↩︎

Figure  3 plots the convergence of different ES implementations(Naive and with EGGROLL) across generations. Figure  4 plots the average wall clock time for vanilla ES and EGGROLL.

Figure 3: Validation accuracy over generations for all methods (mean \pm std, three seeds). Dashed line: surrogate BPTT final accuracy.
Figure 4: Wall clock time(in seconds) for EGGROLL implementations and vanilla ES

EGGROLL with \(r=4\) converges slightly faster than the other EGGROLL variants, with \(r=1\) slightly behind, but overall the difference in convergence speed is modest. The variance bands (mean \(\pm\) std over three seeds) are generally narrow for \(\mathrm{\small Eggroll}\) with small ranks, indicating more stable convergence than several higher-rank or full-rank baselines. EGGROLL achieves substantially lower wall-clock training time due to its reduced per-generation cost. In particular EGGROLL reaches comparable validation accuracy in fewer seconds than vanilla ES despite similar generation-level convergence behavior.

5 Discussion↩︎

5.0.0.1 Accuracy gap.

Eggroll achieves % accuracy vs.% for surrogate BPTT. We attribute the gap to our limited compute budget as a result of which we could only run 100 generations for the ES training(naive and EGGROLL). Additionally with ES, the gradient estimate can be noisy because of stochastic perturbations and mini-batch variation. Due to our limited compute budget, we used a lightweight two-layer fully connected SNN, which keeps the parameter count modest, but limits our representational capacity relative to deeper architectures. Per-seed results (Appendix 8) show occasional training collapse at \(r=2\), \(r=8\), suggesting sensitivity to initialization that warrants further study.

5.0.0.2 Rank and expressiveness.

The rank ablation (Figure 2) shows that even \(r=1\) retains substantial accuracy, suggesting the useful gradient directions lie in a low-dimensional subspace. This is consistent with the flat minima hypothesis [9] and with findings on intrinsic dimensionality in neural networks [10]. The ablations are empirical evidence for this in SNNs.

5.0.0.3 Limitations.

Our experiments are limited to N-MNIST with a two-layer, fully connected network with no convolutions, and 10 fixed time bins which may discard some temporal resolution. Additionally we have not tested this on real neuromorphic hardware, and have not measured the energy consumption. Scaling to CIFAR10-DVS [11] and N-Caltech101  [7] is the next step, since these are more complex than N-MNIST. Three seeds is insufficient for definitive rank claims; we report per-seed results in Appendix  8 for transparency.

6 Related Work↩︎

6.0.0.1 SNN training methods.

Spike-Timing-Dependent Plasticity (STDP) is a local Hebbian rule that requires no global error signal [12] but struggles to match backprop accuracy. Surrogate gradients [4], [13] are currently the dominant approach for deep SNNs, achieving state-of-the-art accuracy but requiring autograd. Conversion methods [14] train an ANN and convert weights, but incur long simulation times.

6.0.0.2 Evolution Strategies for neural networks.

[5] showed ES competitive with RL on Atari. [15] demonstrated ES on deep networks via compact perturbation seeds. PEPG [16] and CMA-ES [17] are closely related natural-gradient variants. Eggroll [6] is, to our knowledge, the first explicit low-rank factorisation of ES perturbations shown to scale to large models; our work is the first to apply it to SNNs.

6.0.0.3 Neuromorphic hardware constraints.

Loihi [1] and BrainScaleS [18] support on-chip learning but only with local rules or very limited precision. Gradient-free methods are a natural fit for these constraints, motivating this line of work.

6.0.0.4 EGGROLL

 [6] showed that EGGROLL can work on language models and RNNs. As far as we know, our work is first application of EGGROLL on SNNs, that poses the challenge of non differentiable spikes.

7 Conclusion↩︎

Training on SNNs is challenging due to its discrete spike threshold. We have demonstrated that Eggroll— low-rank Evolution Strategies — can train SNNs on N-MNIST without surrogate gradients or backpropagation, achieving % test accuracy at \(\times\) the speed of vanilla ES. The rank \(r=4\) offers the best accuracy-per-second tradeoff in our setup.

7.0.0.1 Future work.

The most immediate extension is scaling to more complex neuromorphic benchmarks, like N-Caltech101  [7] and DVS-CIFAR10  [7] with a convolutional SNN. Beyond accuracy, a key open question is whether Eggroll-trained SNNs produce sparser spike trains than surrogate-gradient methods, which would directly translate to energy savings on neuromorphic hardware. Finally, combining Eggroll with local plasticity rules for the early layers — reserving global ES updates for the readout layer — may offer the best of both worlds.

8 Variance Analysis Across Seeds↩︎

Table 2: Per-seed accuracy and average training time across three seeds.
Method Test Acc.(%) Time / gen (s)
Vanilla ES (seed 1) 80.66 12.41
Vanilla ES (seed 2) 78.59 13.15
Vanilla ES (seed 3) 69.26 13.43
\(\mathrm{\small Eggroll}\) \(r=1\) (seed 1) 79.27 3.04
\(\mathrm{\small Eggroll}\) \(r=1\) (seed 2) 78.13 3.59
\(\mathrm{\small Eggroll}\) \(r=1\) (seed 3) 78.86 3.60
\(\mathrm{\small Eggroll}\) \(r=2\) (seed 1) 79.15 6.06
\(\mathrm{\small Eggroll}\) \(r=2\) (seed 2) 79.25 5.94
\(\mathrm{\small Eggroll}\) \(r=2\) (seed 3) 58.42 5.88
\(\mathrm{\small Eggroll}\) \(r=4\) (seed 1) 79.94 5.90
\(\mathrm{\small Eggroll}\) \(r=4\) (seed 2) 78.98 5.82
\(\mathrm{\small Eggroll}\) \(r=4\) (seed 3) 78.72 5.76
\(\mathrm{\small Eggroll}\) \(r=8\) (seed 1) 79.78 5.75
\(\mathrm{\small Eggroll}\) \(r=8\) (seed 2) 69.94 5.64
\(\mathrm{\small Eggroll}\) \(r=8\) (seed 3) 68.68 6.27
Surrogate BPTT (seed 1) 94.35 32.09
Surrogate BPTT (seed 2) 93.88 31.67
Surrogate BPTT (seed 3) 93.86 31.32

9 EGGROLL Gradient Derivation↩︎

Following [6] §4.2, the gradient estimate for weight matrix \(\mathbf{W}\) under low-rank perturbations is: \[\begin{align} \hat{\nabla}_{\mathbf{W}} J &= \frac{1}{2P\sigma\sqrt{r}} \sum_{i=1}^{P} f_i \mathbf{A}_i \mathbf{B}_i^\top \notag \\ &= \frac{1}{2P\sigma\sqrt{r}} \bigl(\operatorname{diag}(\mathbf{f})\,\mathbf{A}\bigr)^\top \mathbf{B}, \label{eq:eggroll95grad} \end{align}\tag{4}\] where \(\mathbf{A} \in \mathbb{R}^{P \times m \times r}\) and \(\mathbf{B} \in \mathbb{R}^{P \times n \times r}\) are stacked factor tensors, and \(f_i = r_i^+ - r_i^-\) is the antithetic rank difference. Equation 4 requires only two batched matrix multiplications, never materialising the \(P \times m \times n\) perturbation tensor.

References↩︎

[1]
M. Davies, N. Srinivasa, T.-H. Lin, et al., “Loihi: A neuromorphic manycore processor with on-chip learning,” IEEE Micro, vol. 38, no. 1, pp. 82–99, 2018.
[2]
A. Vaswani et al., “Attention is all you need,” CoRR, vol. abs/1706.03762, 2017, [Online]. Available: http://arxiv.org/abs/1706.03762.
[3]
P. A. Merolla, J. V. Arthur, R. Alvarez-Icaza, et al., “A million spiking-neuron integrated circuit with a scalable communication network and interface,” Science, vol. 345, no. 6197, pp. 668–673, 2014.
[4]
E. O. Neftci, H. Mostafa, and F. Zenke, “Surrogate gradient learning in spiking neural networks,” IEEE Signal Processing Magazine, vol. 36, no. 6, pp. 51–63, 2019.
[5]
T. Salimans, J. Ho, X. Chen, S. Sidor, and I. Sutskever, “Evolution strategies as a scalable alternative to reinforcement learning,” arXiv preprint arXiv:1703.03864, 2017.
[6]
S. Sarkar et al., “Evolution strategies at the hyperscale,” arXiv preprint arXiv:2511.16652, 2025.
[7]
G. Orchard, A. Jayawant, G. K. Cohen, and N. Thakor, “Converting static image datasets to spiking neuromorphic datasets using saccades,” Frontiers in Neuroscience, vol. Volume 9 - 2015, 2015, doi: 10.3389/fnins.2015.00437.
[8]
D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
[9]
S. Hochreiter and J. Schmidhuber, “Flat minima,” Neural Computation, vol. 9, no. 1, pp. 1–42, 1997.
[10]
C. Li, H. Farkhoor, R. Liu, and J. Yosinski, “Measuring the intrinsic dimension of objective landscapes,” in International conference on learning representations, 2018.
[11]
H. Li, H. Liu, X. Ji, G. Li, and L. Shi, “CIFAR10-DVS: An event-stream dataset for object classification,” Frontiers in Neuroscience, vol. Volume 11 - 2017, 2017, doi: 10.3389/fnins.2017.00309.
[12]
G. Bi and M. Poo, “Synaptic modifications in cultured hippocampal neurons: Dependence on spike timing, synaptic strength, and postsynaptic cell type,” Journal of Neuroscience, vol. 18, no. 24, pp. 10464–10472, 1998.
[13]
G. Bellec, F. Scherr, A. Subramoney, et al., “A solution to the learning dilemma for recurrent networks of spiking neurons,” Nature Communications, vol. 11, no. 1, p. 3625, 2020.
[14]
Y. Cao, Y. Chen, and D. Khosla, “Spiking deep convolutional neural networks for energy-efficient object recognition,” in International Journal of Computer Vision, 2015, vol. 113, pp. 54–66.
[15]
F. P. Such, V. Madhavan, E. Conti, et al., “Deep neuroevolution: Genetic algorithms are a competitive alternative for training deep neural networks for reinforcement learning,” in arXiv preprint arXiv:1712.06567, 2017.
[16]
F. Sehnke, C. Osendorfer, T. Rückstieß, A. Graves, J. Peters, and J. Schmidhuber, The 18th International Conference on Artificial Neural Networks, ICANN 2008“Parameter-exploring policy gradients,” Neural Networks, vol. 23, no. 4, pp. 551–559, 2010, doi: https://doi.org/10.1016/j.neunet.2009.12.004.
[17]
N. Hansen, “The CMA evolution strategy: A tutorial,” arXiv preprint arXiv:1604.00772, 2016.
[18]
J. Schemmel, D. Brüderle, A. Grübl, et al., “A wafer-scale neuromorphic hardware system for large-scale neural modeling,” in Proceedings of the IEEE international symposium on circuits and systems, 2010, pp. 1947–1950.