Spiking Fourier Graph Operators (SpikF-GO)

• Spiking Fourier Graph Operators (SpikF-GO) is a spiking neural architecture for multivariate time series forecasting that transforms each scalar observation into a graph node and applies Fourier-domain mixing.
• The architecture employs a hypervariate graph formulation with spike-based FFT, Hard Concrete frequency gating, and Complex LIF units to capture temporal and cross-variable dependencies.
• Empirical results on eight benchmarks demonstrate that SpikF-GO achieves state-of-the-art forecasting performance with improved energy efficiency compared to traditional ANN models.

Spiking Fourier Graph Operators (SpikF-GO) are a spiking neural architecture for multivariate time series forecasting that combines a hypervariate graph formulation—where every scalar observation becomes a graph node—with spike-driven spectral processing in the Fourier domain. The model introduces a Hard Concrete frequency gate for learnable sparse frequency selection and a Complex LIF gate that applies independent spiking neurons to real and imaginary Fourier components, with a variant that injects Central Pattern Generator-based positional encodings prior to spectral mixing. In the reported evaluation on eight benchmarks under a unified experimental protocol, SpikF-GO achieves the best average rank among all SNN methods and outperforms its ANN counterpart, FourierGNN, at reduced energy cost, while remaining competitive at substantially smaller embedding dimensions (Bakhshaliyev et al., 11 Jun 2026).

1. Historical positioning and problem formulation

SpikF-GO was introduced in the context of a specific limitation in spiking time-series forecasting: most existing SNN-based TSF models process variables independently, so a multivariate input matrix $X\in\mathbb{R}^{L\times N}$ is handled effectively as $N$ independent univariate series. The reported consequence is that cross-variable structure—spatial correlations, co-movement, and time-varying interactions—is left implicit rather than being modeled by an explicit graph mechanism. On the ANN side, graph-based models such as STGCN, AGCRN, StemGNN, MTGNN, and FourierGNN explicitly encode inter-variable dependencies via adjacency matrices and graph operations, with FourierGNN being particularly relevant because it treats every scalar value $x_{t,n}$ in the input window as a graph node and performs graph mixing by multiplying in Fourier space over the flattened node axis (Bakhshaliyev et al., 11 Jun 2026).

Within that landscape, SpikF-GO combines FourierGNN’s hypervariate graph perspective with a spike-based computation pipeline. The model adopts the view that each scalar observation at each time step is a node, performs graph mixing purely by spectral operators in the Fourier domain, and replaces continuous ANN computation by spike-driven computation. The inputs are encoded as spike trains; graph mixing is implemented with a Spiking FFT, complex linear operators, and Complex LIF gates in the Fourier domain; and frequency sparsity is enforced by a Hard Concrete gate. The paper describes this as a graph-based SNN that explicitly models intra-series temporal dependencies, inter-series correlations, and time-varying cross-variable interactions in a unified structure, and as a Fourier-based SNN that preserves binary, event-driven computation for real and imaginary Fourier components while pruning frequencies (Bakhshaliyev et al., 11 Jun 2026).

The work also positions itself against earlier spiking forecasting architectures. SpikeRNN, SpikeTCN, Spikformer, TS-GRU, TS-TCN, TS-Former, TS-LIF-based models, and SpikF are treated as important baselines, but not as graph-based multivariate models. The authors state that SpikF-GO is the first to bring graph-based multivariate modeling into the spiking domain for TSF and the first to provide a unified comparison across SNN forecasting architectures under a common experimental protocol (Bakhshaliyev et al., 11 Jun 2026).

2. Hypervariate graph construction and spike encoding

The model starts from a batch input

$$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$

where $B$ is batch size, $L$ is the look-back window, and $N$ is the number of variables. SpikF-GO then applies a hypervariate graph encoder that flattens the $(L,N)$ structure into

$M = L N$

nodes by a map

$$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$

For each batch element,

$N$0

Each flattened index $N$1 corresponds to a specific pair $N$2, so the node set spans both time and variable axes (Bakhshaliyev et al., 11 Jun 2026).

Embedding is performed by broadcasting a shared learnable embedding vector $N$3 over nodes: $N$4 This is followed by an affine transform $N$5 and RMSNorm,

$N$6

A defining architectural point is that no explicit adjacency matrix or Laplacian is constructed. Instead, the Fourier-domain linear maps across the node axis implicitly define the graph operator, in a manner described as conceptually analogous to a graph convolution with a data-driven adjacency (Bakhshaliyev et al., 11 Jun 2026).

The continuous TS time step $N$7 is aligned with the SNN time step $N$8 through

$N$9

so each time-series step is subdivided into $x_{t,n}$0 spiking steps. For each SNN step $x_{t,n}$1,

$x_{t,n}$2

with scalar parameters $x_{t,n}$3. A LIF encoder then generates the spike tensor

$x_{t,n}$4

This tensor constitutes the input to the spectral block (Bakhshaliyev et al., 11 Jun 2026).

A common misconception is that the “graph” in SpikF-GO refers to a predefined spatial graph in the usual GNN sense. In fact, the formulation is minimalistic and nonparametric: every scalar observation becomes a node, and graph mixing is induced implicitly by Fourier-domain operators over the flattened node axis rather than by an explicitly specified adjacency matrix (Bakhshaliyev et al., 11 Jun 2026).

3. Spectral graph mixing and Complex LIF dynamics

The central operator block is the Spiking Fourier Graph Operator, or S-FGO. The spike tensor $x_{t,n}$5 is transformed along the node dimension $x_{t,n}$6 by a Spiking FFT: $x_{t,n}$7 At fixed SNN step, batch index, and channel, the underlying transform is the Fourier transform over the hypervariate node axis,

$x_{t,n}$8

Because $x_{t,n}$9 is applied along the node axis, any elementwise complex linear transform in the Fourier domain corresponds to a linear operator in the node domain. The subsequent complex linear layers therefore act as Fourier Graph Operators over the hypervariate graph, with log-linear complexity in $$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$0 due to FFT usage (Bakhshaliyev et al., 11 Jun 2026).

The neuron model is a discrete-time LIF unit with charging, thresholding, and reset-leak dynamics: $$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$1

$$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$2

$$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$3

Here $$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$4 is input current, $$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$5 is the decay factor, $$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$6 is the threshold, and $$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$7 is the reset voltage. Training uses an arctangent surrogate gradient,

$$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$8

with slope parameter $$\mathbf{X}_{\text{in}} \in \mathbb{R}^{B \times L \times N},$$9 (Bakhshaliyev et al., 11 Jun 2026).

A distinctive contribution is the Complex LIF gate, which acts directly on complex Fourier features. For a complex tensor

$B$0

the gate is defined by

$B$1

The real and imaginary parts are fed to independent LIF neuron populations. The resulting binary mask is the logical OR of the two spike signals, so a complex feature is retained if either real or imaginary neuron spikes and is zeroed otherwise. This preserves event-driven computation within the spectral domain itself (Bakhshaliyev et al., 11 Jun 2026).

Inside each spectral layer, the model applies complex affine normalization, complex linear transformation, and Complex LIF gating in a stacked form: $B$2 Residual connections with learned scaling are also included. The paper states the presence of residuals without writing a full explicit residual equation; a plausible implication is that the spectral block is organized as a complex residual stack rather than as an isolated feedforward transform (Bakhshaliyev et al., 11 Jun 2026).

After spectral processing, the model applies a spiking inverse FFT, reshapes the output back from $B$3 nodes to $B$4, performs linear temporal compression from $B$5 to $B$6, and then applies LIF, linear projection, time averaging over SNN steps, GELU, and a final linear layer to produce forecasts

$B$7

The decoder therefore remains mostly spike-based but is not purely spiking, since it includes linear projections and GELU (Bakhshaliyev et al., 11 Jun 2026).

4. Frequency gating, sparsity control, and CPG augmentation

SpikF-GO introduces a Hard Concrete frequency gate because the number of Fourier bins,

$B$8

can be large when $B$9. The gate assigns one parameterized stochastic binary relaxation to each frequency bin $L$0, with learnable log-odds $L$1. During training, the continuous relaxed gate is obtained by sampling $L$2, applying a binary-concrete transform with temperature $L$3, stretching to $L$4, and clipping to $L$5: $L$6

$L$7

The spectral tensor is then masked as

$L$8

broadcast over $L$9. At inference, deterministic binarization is used: $N$0 This yields a fixed sparse set of active frequencies during deployment (Bakhshaliyev et al., 11 Jun 2026).

The sparsity objective is an $N$1-style penalty on the expected number of active gates,

$N$2

combined with MSE in the total loss

$N$3

The paper reports that replacing this gate with a fixed Top-$N$4 frequency selection degrades performance, which is presented as evidence that data-driven frequency sparsity is more effective than heuristic selection (Bakhshaliyev et al., 11 Jun 2026).

The variant SpikF-GO w/ CPG adds Central Pattern Generator-based positional encodings before spectral mixing. The paper states that it uses the exact CPG module from Lv et al. as implemented in the original work. Although the CPG equations are not reproduced, the role is explicit: it injects positional signal into the spiking representation so as to strengthen long-range temporal modeling. Empirically, the CPG-enhanced variant consistently outperforms plain SpikF-GO on most datasets and achieves the best average rank. This suggests that the spectral block alone does not eliminate the need for positional structure in spiking forecasting, especially at larger horizons (Bakhshaliyev et al., 11 Jun 2026).

5. Optimization protocol, efficiency profile, and empirical performance

The forecasting objective is the standard mean squared error over batch, horizon, and variables: $N$5 Evaluation uses

$N$6

Gradient propagation proceeds through LIF neurons via the arctangent surrogate, through the Complex LIF gate via a Straight-Through Estimator treatment of the logical OR and masking, and through the Hard Concrete gate by reparameterization in $N$7 (Bakhshaliyev et al., 11 Jun 2026).

The implementation is reported in PyTorch 2.5.1 on NVIDIA RTX 4090. All models use Reversible Instance Normalization, and the study corrects a data leakage bug in the original FourierGNN code by using training-set statistics only. All SNN models use $N$8 steps except the SpikF baseline, which uses $N$9 as in the original model. SpikF-GO uses $(L,N)$0 spectral layers and embedding size $(L,N)$1 by default, with the CPG variant using the same structural hyperparameters (Bakhshaliyev et al., 11 Jun 2026).

The energy argument rests on three factors: spike-based operations instead of full-precision MACs, sparse frequency selection by the Hard Concrete gate, and the possibility of drastically shrinking the embedding dimension. Using the assumptions from Yao et al. on 45nm technology, the study adopts 4.6 pJ per FLOP and 0.9 pJ per SOP. On Solar with $(L,N)$2 and $(L,N)$3, the reported theoretical energy estimates are:

• FourierGNN: $(L,N)$4,
• SpikF-GO ($(L,N)$5): $(L,N)$6, i.e. $(L,N)$7 less,
• SpikF baseline: $(L,N)$8, i.e. $(L,N)$9 less than FourierGNN,
• SpikF-GO with $M = L N$0: $M = L N$1, i.e. $M = L N$2 less than FourierGNN (Bakhshaliyev et al., 11 Jun 2026).

The eight benchmarks are ECG, COVID-19, Solar, Electricity, Traffic, METR-LA, PEMS-BAY, and Wiki. The main tables focus on input length 12 and forecast length 12. On average rank across all eight datasets, SpikF-GO w/ CPG is reported as best on both metrics, with $M = L N$3 rank 2.4 and MAE rank 2.3. SpikF-GO is second-best in average $M = L N$4 with rank 2.8 and has MAE rank 3.8. FourierGNN has average rank 4.8 for $M = L N$5 and 3.1 for MAE. The SpikF baseline is reported at approximately 4.4 for $M = L N$6 and 4.3 for MAE among SNNs (Bakhshaliyev et al., 11 Jun 2026).

Several dataset-specific results are highlighted. On COVID-19, FourierGNN achieves $M = L N$7 and MAE $M = L N$8, SpikF reaches $M = L N$9 and MAE $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$0, and SpikF-GO improves to $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$1 and MAE $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$2. On Traffic, FourierGNN attains $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$3 and MAE $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$4, while SpikF-GO w/ CPG reaches $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$5 and MAE $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$6. On PEMS-BAY, FourierGNN reports $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$7 and MAE $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$8, whereas SpikF-GO w/ CPG reaches $$\mathcal{H}:\mathbb{R}^{B\times L\times N}\to\mathbb{R}^{B\times M}.$$9 and ties MAE at $N$00 (Bakhshaliyev et al., 11 Jun 2026).

The ablation study isolates three design choices. A Temporal-Only variant, which removes cross-variable modeling and processes each variable independently, drops $N$01 from $N$02 on Solar, $N$03 on Traffic, and $N$04 on COVID-19. Replacing the Hard Concrete gate by fixed Top-$N$05 frequency selection lowers performance, for example from Solar $N$06 and COVID $N$07. Replacing RMSNorm by Scale-Shift leaves Solar essentially unchanged but yields slight drops on Traffic and COVID-19. Hyperparameter sensitivity analyses show that performance improves up to 8 or 12 spiking steps and then plateaus, that there is little gain beyond $N$08, and that $N$09 remains stable from embedding size $N$10 to $N$11 (Bakhshaliyev et al., 11 Jun 2026).

A recurrent point of interpretation is that GPU runtime is longer for SNN models because multiple SNN steps must be simulated. The paper treats this as an ANN-hardware artifact rather than as a contradiction of the efficiency claim, arguing that on neuromorphic chips the SOP-versus-FLOP advantage is the relevant quantity. A necessary qualification is that the reported energy analysis is theoretical and derived from operation counts rather than from direct measurement on actual neuromorphic hardware (Bakhshaliyev et al., 11 Jun 2026).

6. Mathematical context, related graph Fourier frameworks, and limitations

SpikF-GO belongs to a broader line of work on graph Fourier operators, but it instantiates that line in a distinctive way. Unlike Laplacian-eigenbasis graph signal processing, the model does not construct an explicit adjacency matrix, Laplacian, or graph shift operator. Its “graph” is the hypervariate flattening of the $N$12 input, and its graph operator is realized by Fourier-domain mixing along that node axis. This differs from projector-based GFT formulations on defective adjacency matrices, where the spectral domain is a direct sum of Jordan subspaces with spectral projectors $N$13, and from unitary graph shift formulations in which a graph discrete Fourier transform is defined by an orthogonal or unitary graph shift operator (Deri et al., 2017, Dees et al., 2019).

At the algorithmic level, SpikF-GO also differs from earlier efforts to obtain fast graph Fourier transforms via sparse factorizations. “Approximate fast graph Fourier transforms via multi-layer sparse approximations” develops a greedy approximate diagonalization of the graph Laplacian using sparse orthogonal factors built from Givens rotations, yielding an FFT-like, layerwise structure for approximate graph Fourier transforms. That construction provides an explicit precedent for interpreting graph spectral operators as multi-layer sparse networks, but SpikF-GO does not approximate a Laplacian eigenbasis in that manner; instead, it performs S-FFT on the hypervariate node axis and learns complex spectral operators directly in the resulting frequency space (Magoarou et al., 2016).

These distinctions matter because the term “Fourier graph operator” can invite incompatible intuitions. In classical graph signal processing, the Fourier transform is usually tied to an eigendecomposition of a graph operator such as a Laplacian, an adjacency matrix, a spectral projector family, or a unitary shift. In SpikF-GO, by contrast, the operator is “graph-based” because the node set is the hypervariate graph over all scalar observations and because Fourier-domain mixing over that node axis induces a learned graph operator in the original node domain. This suggests a hybrid status: the model is graph-structured and spectral, but not a direct implementation of a canonical GFT in the sense of those earlier theories (Bakhshaliyev et al., 11 Jun 2026).

The reported limitations are correspondingly practical rather than foundational. First, the node count is $N$14, so very large $N$15 or $N$16 can make S-FFT and spectral layers heavy even when their asymptotic cost is log-linear in $N$17. Second, the model is global over all nodes and does not exploit locality or multi-scale graph hierarchies. Third, complex-valued spiking requires two LIF populations per complex feature, which makes implementation in standard SNN toolkits more involved. Fourth, SNN training with BPTT and surrogate gradients remains relatively costly and sensitive to hyperparameters, and Hard Concrete gating adds stochasticity that requires careful selection of temperatures and priors. Finally, the energy evaluation is theoretical rather than hardware-measured, so direct validation on neuromorphic platforms remains an open requirement (Bakhshaliyev et al., 11 Jun 2026).

The extension directions identified by the authors include scaling to high-dimensional, long-horizon settings through hierarchical graph decompositions or multi-resolution spectral methods, and validating the architecture on real neuromorphic platforms. A plausible implication is that future versions of SpikF-GO may converge toward more explicit graph operator hierarchies, potentially bringing it closer to the sparse multi-layer spectral constructions studied elsewhere in graph Fourier research (Bakhshaliyev et al., 11 Jun 2026, Magoarou et al., 2016).