Parallel Phase Oscillator SSM
Implementation & benchmark of arXiv:2608.07754 — Phase State Space Models: Parallel, Surrogate-Free Training of Spiking Networks
Source paper: Phase State Space Models: Parallel, Surrogate-Free Training of Spiking Networks arXiv:2608.07754 ⓘ · analyzed Sep 1, 2026
AI-generated research hypothesis, automatically tested. Not peer-reviewed.
Idea description
Replace real diagonal state-space channels with complex damped oscillators whose hidden states encode both amplitude and phase. Train with parallel causal convolution and deploy with the equivalent one-step recurrence, allowing the same layer to support efficient batched training and low-memory streaming inference.
Formulas
Mathematical statement
The paper defines resonate-and-fire dynamics as dU/dt = k U(t) + w I(t), where U(t) is a complex hidden potential, I(t) is a real input current, w is a complex input weight, and k is the complex state coefficient. It sets k = lambda + i omega, where lambda < 0 is damping, omega is angular frequency, and i is the imaginary unit. For a timestep Delta and piecewise-constant input I_t, exact integration gives U_{t+1} = a U_t + b I_t, with a = exp(k Delta) and b = w (exp(k Delta) - 1) / k when k is nonzero. Since |a| = exp(lambda Delta) < 1, the recurrence is stable. Unrolling gives U_t = a^t U_0 + sum from j = 0 to t-1 of a^(t-1-j) b I_j, which is a causal convolution with kernel h_r = a^r b. The oscillatory factor exp(i omega Delta) stores phase, while exp(lambda Delta) controls memory decay.
Implementation notes
(1) Integration point: add the layer after an input projection in a sequence model, replacing or augmenting a real diagonal SSM. For input X with shape batch by time by dimension, project X to real currents I with M oscillator channels. Maintain a complex state U for each channel, then map the concatenated real and imaginary parts back to the model dimension with a linear projection.
(2) Pseudocode:
rho, omega, w = learned_parameters[M]
lambda = -softplus(rho)
k = lambda + 1j * omega
a = exp(k * dt)
b = where(abs(k) > eps, w * (a - 1) / k, w * dt)
I = input_projection(X)
if training:
h[r] = a**r * b for r in range(T)
U = causal_fft_convolution(I, h) + a**arange(T) * U0
else:
U = U0
for t in range(T):
U = a * U + b * I[:, t]
save(U)
Y = output_projection(concat(real(U), imag(U)))
(3) Compute a and b analytically from rho, omega, w, and dt. Compute the training convolution using FFT with zero-padding to at least twice the sequence length and retain only causal positions. Monitor max absolute a during training; the parameterization lambda = -softplus(rho) enforces stability. For a real-valued implementation, either concatenate real and imaginary parts or use conjugate oscillator pairs with frequencies omega and -omega.
(4) First experiment: train a 4-layer, 128-channel model on Sequential MNIST or Speech Commands and compare it with a same-parameter real diagonal SSM and a GRU. Use identical optimizer, batch size, and hidden-state width. Measure validation accuracy, training throughput, peak memory, and streaming tokens or timesteps per second. Success means better accuracy on long periodic dependencies than the real SSM at equal parameters, faster training than recurrent unrolling, and stable loss without gradient explosions. Ablate learned versus fixed omega and compare unconstrained frequencies with the Nyquist constraint abs(omega) times dt <= pi.
Verification
This idea has not been verified yet.
Verification happens in two stages: Stage 1 — a mechanism check on a toy system confirms the claimed mathematical phenomenon reproduces; Stage 2 — a benchmark implements the idea on a real (small) neural network task and compares it against a tuned baseline over 8 paired seeds with a permutation test.
Artifacts
Artifacts unavailable.