Bullet Trains: Parallelizing Training of Temporally Precise Spiking Neural Networks

Todd Morrill; Christian Pehle; Anthony Zador

cs.NEFeb 27, 2026

Bullet Trains: Parallelizing Training of Temporally Precise Spiking Neural Networks

Authors: Todd Morrill, Christian Pehle, Anthony Zador

Organizations: Department of Computer Science, Columbia University, New York, NY, USA · Department of Neuroscience, Cold Spring Harbor Laboratory, Cold Spring Harbor, NY, USA

Abstract

Continuous-time, event-native spiking neural networks (SNNs) operate strictly on spike events, treating spike timing and ordering as the representation rather than an artifact of time discretization. This viewpoint aligns with biological computation and with the native resolution of event sensors and neuromorphic processors, while enabling compute and memory that scale with the number of events. However, two challenges hinder practical, end-to-end trainable event-based SNN systems: 1) exact charge--fire--reset dynamics impose inherently sequential processing of input spikes, and 2) precise spike times must be solved without time bins. We address both. First, we use parallel associative scans to consume multiple input spikes at once, yielding up to 44x speedups over sequential simulation while retaining exact hard-reset dynamics. Second, we implement differentiable spike time solvers that compute spike times to machine-precision without discrete-time approximations or restrictive analytic assumptions. We demonstrate the viability of training SNNs using our solutions on four event-based datasets on GPUs.

Figures & tables

Figure 1: Our parallel method achieves up to 44x speedups over serial processing of spike events on several event-based datasets while retaining exact hard-reset dynamics. The plot shows results for the Spiking Heidelberg Digits (SHD) dataset.

Figure 2: Our contributions visualized—state trajectory of a LIF neuron consuming input spikes in parallel and producing an output spike with our Newton-Raphson solver. A shows the neuron jumping from the initial state

\mathbf{s}_{0}

at time 0 (composed of membrane potential and synaptic current) directly to all future states

\mathbf{s}_{1}=f_{1}(\mathbf{s}_{0}),\mathbf{s}_{2}=f_{2}(\mathbf{s}_{0}),\ldots

in parallel using associative scans. We use lightweight analytical checks to determine if an output spike will occur in the interval between any consecutive pair of input spikes by computing—in parallel—

t_{V_{\max}}(\mathbf{s})

, the time of maximum voltage starting at the interval’s left endpoint, and

V(t_{V_{\max}})

the voltage at that time. The voltage state is hard reset at the spike time. B shows our iterative Newton-Raphson solver finding the output spike time

t^{\star}

, computed to machine-precision.

Figure 3: Discrete-time binning causes quantization errors and loses spike ordering. Here, two spikes occurring at different times (0.0021s and 0.0027s) both round to 0.003s, making them indistinguishable to downstream neurons.

Figure 4: An example input queue feeding a neuron and producing an output queue. The input chunk size is 4 spikes, which means 4 input spikes are consumed in parallel. After consuming 3 input spikes in chunk 1, the neuron produces an output spike. The work done to consume the 4th input spike in chunk 1 is discarded, and the next chunk resumes processing input spikes immediately after the output spike time.

Figure 5: Spike time as a function of synaptic weight for different time constant configurations. Discrete-time approximation with

\Delta t=1\text{ms}

diverges from the ground truth spike time at nearly all weight values, while our method closely tracks the ground truth spike time.

Figure 6: Training time per batch as we vary chunk size on the SHD dataset with hidden dimensions ranging from 128 to 512 and batch size ranging from 32 to 128. Averages are taken over 100 batches.

Figure 7: Convergence of the Newton-Raphson and Bisection solvers towards the ground truth spike time.

Figure 8: Speedup of our parallel SNN implementation over a sequential event-based SNN implementation with exact hard-reset dynamics on the SHD dataset as we vary the hidden dimension and batch size. The model architecture uses 2 feedforward hidden layers and delays. Averages are taken over 100 batches (forward and backward pass) over 3 runs and chunk size is set to 128.

Dataset	Model	Exact Gradients	Continuous Spike Times	Parallelized	Compute Depth	Memory	Accuracy
MNIST	1F350H, $\tau_{m}=2\tau_{s}$ ( Göltz et al., 2021 )	✓	✓	x	$O(N)$	$O(N)$	97.20 (0.10)
	1F350H ( Wunderlich & Pehle, 2021 )	✓	✓	x	$O(N)$	$O(N)$	97.60 (0.10)
	1F350H (Ours)	✓	✓	✓	$O(C\log{K})$	$O(N)$	98.04 (0.07)
SHD	2F256HD ( Hammouamri et al., 2024 )	x	x	x	$O(T)$	$O(T)$	95.07 (0.24)
	2F512HD ( Mészáros et al., 2025 )	✓	x	x	$O(T)$	$O(N)$	93.10
	2F256HD (Ours)	✓	✓	✓	$O(C\log{K})$	$O(N)$	94.77 (0.46)

Table 1: Classification accuracy comparison on event-based datasets. Model architectures are defined by number of layers and connectivity pattern (F for feedforward, R for recurrent), number of hidden units (H), and use of delays (D) (e.g., 2F512HD means 2 layers, feedforward connectivity, 512 hidden units, and use of delays).

N

is the number of spikes processed by a single neuron and

T

is the number of time steps in the simulation. Our metrics are averaged over 5 random seeds and standard deviation is reported in parentheses.

Figure 9: Impact of temporal discretization on Yin-Yang classification. We trained our method under different discrete-time constraints by quantizing spike times to multiples of

\Delta t

during training and testing. Accuracy degrades significantly as

\Delta t

increases, dropping to near-chance (

33\%

) at

\Delta t\geq 2

ms. Continuous spike times (

\Delta t\to 0

, our standard method) maintain the highest accuracy, demonstrating that discrete-time methods would require

\Delta t\ll 1

ms—incurring substantial computational costs—to preserve temporal information.

Appendix figures & tables8 assets

Supplementary material from the paper’s appendix.

Appendix

Parameter	Yin-Yang Value	MNIST Value	SHD Value	SSC Value
Optimizer	Adam	Adam	Adam	Adam
Learning Rate Schedule	Cosine Decay	Exponential	Cosine Decay	Cosine Decay
Cosine Decay Steps	6000	—	10000	88500
Learning Rate	0.02	0.02	0.001	0.001
Learning Rate End Value	0.0001	0.001	0.0001	0.0001
Learning Rate Decay Factor	—	0.95	—	—

Appendix

Table 2: System hyperparameters across the 3 datasets used.

Figure 10: Speedup of our parallel SNN implementation over a sequential event-based SNN implementation with exact hard-reset dynamics on the SSC dataset as we vary the hidden dimension and batch size. The model architecture uses 2 feedforward hidden layers and delays. Averages are taken over 100 batches (forward and backward pass) over 3 runs and chunk size is set to 128.

Hidden Layer Size	Batch Size	Computation Mode	Avg. Runtime (std. dev.)
128	32	Parallel	0.193 (0.003)
	32	Serial	8.075 (0.212)
	64	Parallel	0.284 (0.003)
	64	Serial	8.195 (0.033)
	128	Parallel	0.457 (0.004)
	128	Serial	8.510 (0.067)

Appendix

Table 3: Average training time (forward and backward pass) per batch in seconds on the SSC dataset for our parallel SNN implementation versus a sequential event-based SNN implementation with exact hard-reset dynamics. The model architecture uses 2 feedforward hidden layers and delays. Averages are taken over 100 batches (forward and backward pass) over 3 runs and chunk size is set to 128.

Hidden Layer Size	Batch Size	Computation Mode	Avg. Runtime (std. dev.)
128	32	Parallel	0.189 (0.001)
	32	Serial	8.352 (0.022)
	64	Parallel	0.277 (0.002)
	64	Serial	8.167 (0.044)
	128	Parallel	0.449 (0.003)
	128	Serial	8.614 (0.034)

Appendix

Table 4: Average training time (forward and backward pass) per batch in seconds on the SHD dataset for our parallel SNN implementation versus a sequential event-based SNN implementation with exact hard-reset dynamics. The model architecture uses 2 feedforward hidden layers and delays. Averages are taken over 100 batches (forward and backward pass) over 3 runs and chunk size is set to 128.

Method	Epoch Training Time (s) $\downarrow$	Max Memory (MiB) $\downarrow$
256 ( Hammouamri et al., 2024 )	157.89	8,978
512 ( Hammouamri et al., 2024 )	338.94	16,366
256 ( Mészáros et al., 2025 )	38.93	4,824
512 ( Mészáros et al., 2025 )	141.35	8,632
256 (Ours)	44.66	7,967
512 (Ours)	111.55	26,184

Appendix

Table 5: Absolute speed and memory comparison of our method to prior work on the SHD dataset with batch size 256. All methods use 2 hidden layers with delays and 256 or 512 hidden units per layer. * denotes metrics reported in Mészáros et al. (2025) on an RTX A5000; all other numbers were measured on an NVIDIA H100 (80GB). See text for details.

Figure 11: Spike metrics for 5 runs on the SSC dataset with 512 hidden units and chunk size 128. On average, layer one and two neurons consume approximately 100% of input spikes. We also note that both layer one and layer two retain a very high percentage of work done (e.g., 80-90%+).

Figure 12: Spike metrics for 5 runs on the SHD dataset with 256 hidden units and chunk size 128. On average, layer one and two neurons consume approximately 100% of input spikes. We also note that both layer one and layer two retain a very high percentage of work done (e.g., 80-90%+).

Figure 13: Spike metrics for 5 runs on the SHD dataset with 512 hidden units and chunk size 128. On average, layer one and two neurons consume approximately 100% of input spikes. We also note that both layer one and layer two retain a very high percentage of work done (e.g., 80-90%+).

Spiking Neural Networks

Spiking neural networks process information through discrete, time-dependent neural events, aiming for efficient computation and biologically inspired learning.

See all 336 papers

Scalable Training of Continuous-Time Spiking Neural Networks with Differentiable Spike-Time Discretization
Jul 16, 2026Yusuke Sakemi, Tomoya Takeuchi, Takeo Hosomi +1Spiking Neural Networks LIF Neurons
Research Center for Mathematical Engineering, Chiba Institute of Technology, Narashino, Japan · International Research Center for Neurointelligence (WPI-IRCN), The University of Tokyo, Tokyo, Japan · NEC Corporation, Kawasaki, Japan
FidLIF: Parallel Training of High-Fidelity Spiking Neurons with Exact BPTT
Oct 12, 2026Kaixiang He, Yujie WuLIF Neurons Spiking Neural Networks
Department of Computing, The Hong Kong Polytechnic University
Multi-Timescale Conductance Spiking Networks: A Sparse, Gradient-Trainable Framework with Rich Firing Dynamics for Enhanced Temporal Processing
May 12, 2026Alex Fulleda-Garcia, Saray Soldado-Magraner, Josep Maria Margarit-TauléSpiking Neural Networks Neuromorphic Computing
Instituto de Microelectrónica de Barcelona (IMB-CNM), Consejo Superior de Investigaciones Cientíﬁcas (CSIC) Cerdanyola del Vallès, Spain · Department of Neurobiology, University of California Los Angeles (UCLA) Los Angeles, USA · Joint Senior authors

Jul 16, 2026J/K move · Enter open · S save

Scalable Training of Continuous-Time Spiking Neural Networks with Differentiable Spike-Time Discretization

Yusuke Sakemi, Tomoya Takeuchi, Takeo Hosomi, Kazuyuki Aihara

Research Center for Mathematical Engineering, Chiba Institute of Technology, Narashino, Japan · International Research Center for Neurointelligence (WPI-IRCN), The University of Tokyo, Tokyo, Japan · NEC Corporation, Kawasaki, Japan

Spiking Neural Networks LIF Neurons

Open paper AskPDF arXiv

Continuous-time spiking neural networks (SNNs) provide an event-driven framework for temporal computation, computational neuroscience, and neuromorphic hardware. However, training deep continuous-time SNNs is severely constrained by the memory required for exact spike-time computation, which evaluates and retains candidate firing times over intervals determined by presynaptic spike ordering. Here we introduce a memory-efficient training framework based on differentiable spike-time discretization (DSTD) for leaky integrate-and-fire neurons with general membrane and synaptic time constants. DSTD maps irregular presynaptic spikes onto differentiable weighted events at fixed time points, replacing the input-dependent candidate dimension with

M

fixed time intervals while accurately approximating continuous-time membrane-potential dynamics. This reduces candidate-related activation memory from

O(N_{\mathrm{out}}N_{\mathrm{in}})

O(N_{\mathrm{out}}M)

in the case of time-to-first-spike (TTFS) coding, where

N_{\mathrm{in}}

and

N_{\mathrm{out}}

denote the numbers of presynaptic and postsynaptic neurons, respectively. We further introduce synfire-chain-inspired temporal regularization that organizes layer-wise firing windows, mitigates dead-neuron failures, and enables pipeline-like processing. In dense LIF layers, DSTD reduced peak memory consumption by up to approximately 100-fold and training time by up to approximately 20-fold compared with exact spike-time computation. Together, these methods allowed us to train 9-layer convolutional SNNs on CIFAR-10 and 20-layer convolutional SNNs on Fashion-MNIST on a single GPU.

Abstract

Spiking Neural Networks

Explore similar work

Scalable Training of Continuous-Time Spiking Neural Networks with Differentiable Spike-Time Discretization

FidLIF: Parallel Training of High-Fidelity Spiking Neurons with Exact BPTT

Multi-Timescale Conductance Spiking Networks: A Sparse, Gradient-Trainable Framework with Rich Firing Dynamics for Enhanced Temporal Processing