Grunwald-Letnikov Fractional LIF

The Grunwald-Letnikov (GL) Fractional Leaky Integrate-and-Fire model is the primary fractional neuron implementation in SPIRES. It provides an efficient and numerically robust discretization of fractional-order dynamics using precomputed binomial coefficients and a circular buffer for voltage history.

The GL Fractional Derivative

The Grunwald-Letnikov fractional derivative of order $\alpha$ is defined as a limit of a generalized difference quotient:

where the generalized binomial coefficient is:

For $\alpha = 1$, this reduces to the standard backward difference: $D^1 x(t) \approx (x(t) - x(t - \Delta t)) / \Delta t$.

For $0 \lt \alpha \lt 1$, the sum extends over the entire history of $x$, with weights that decay as a power law. The GL derivative is equivalent to the Caputo derivative for sufficiently smooth functions with consistent initial conditions.

GL Discretization of the FLIF

Applying the GL derivative to the fractional LIF equation:

yields the discrete update rule. Isolating $V_n$ on the left-hand side:

where the GL coefficients are:

These coefficients satisfy the recurrence relation:

Properties of the GL Coefficients

The GL coefficients have important properties that govern the memory behavior of the neuron:

• $c_0(\alpha) = 1$ for all $\alpha$
• $c_1(\alpha) = -\alpha$
• For $0 \lt \alpha \lt 1$, the coefficients alternate in sign and decay as $|c_k| \sim k^{-\alpha-1}$ for large $k$
• The partial sums $\sum_{k=0}^{L} c_k(\alpha)$ converge, which ensures that the truncated history approximation is well-behaved

The power-law decay of the coefficients is the discrete manifestation of the long memory property: contributions from the distant past never fully vanish but become progressively smaller.

Implementation Details

Precomputed Coefficients

At reservoir creation time, SPIRES precomputes the array $\{c_k(\alpha)\}_{k=1}^{L}$ using the recurrence relation. This costs $O(L)$ operations and is done once. The coefficients are stored in a contiguous array for cache-efficient access during the per-step update.

Circular Buffer

Each neuron maintains a circular buffer of length $L$ to store its recent voltage history $\{V_{n-1}, V_{n-2}, \ldots, V_{n-L}\}$. The circular buffer avoids the cost of shifting elements at every time step:

• Write: Store $V_n$ at position $n \bmod L$.
• Read: Access $V_{n-k}$ at position $(n - k) \bmod L$.

This gives $O(1)$ bookkeeping per step (beyond the $O(L)$ history sum).

Per-Step Update Algorithm

At each time step $n$ for neuron $i$:

1. Compute the synaptic input current $I_n^{(i)}$ from incoming spikes.
2. Compute the history sum: $S = \sum_{k=1}^{\min(n,L)} c_k(\alpha) \cdot V_{n-k}^{(i)}$
3. Update: $V_n^{(i)} = \Delta t^\alpha \left(-\frac{V_{n-1}^{(i)} - V_{\text{rest}}}{\tau_m} + I_n^{(i)} + b\right) - S$
4. If $V_n^{(i)} \geq V_{\text{th}}$: emit spike, reset $V_n^{(i)} = V_{\text{reset}}$
5. Write $V_n^{(i)}$ into the circular buffer.

History Length Trade-Off

The truncation length $L$ controls the trade-off between accuracy and computational cost:

The truncation error from using $L$ terms instead of the full history decays as $O(L^{-\alpha})$. For $\alpha = 0.5$, doubling $L$ reduces the error by a factor of $2^{0.5} \approx 1.41$. For $\alpha = 0.9$, the decay is faster ($2^{0.9} \approx 1.87$), meaning shorter histories suffice.

Parameters

SPIRES API

cfg.neuron_type = SPIRES_NEURON_FLIF_GL;

Example Configuration

spires_flif_params params = {
    .alpha   = 0.5,
    .tau_m   = 20.0,
    .v_rest  = -65.0,
    .v_th    = -50.0,
    .v_reset = -65.0,
    .bias    = 0.0,
    .L       = 200,
};
 
spires_reservoir_config cfg = {
    .num_neurons       = 500,
    .num_inputs        = 1,
    .num_outputs       = 1,
    .spectral_radius   = 0.95,
    .ei_ratio          = 0.8,
    .input_strength    = 0.1,
    .connectivity      = 0.1,
    .dt                = 1.0,
    .connectivity_type = SPIRES_CONN_RANDOM,
    .neuron_type       = SPIRES_NEURON_FLIF_GL,
    .neuron_params     = &params,
};

Sweeping the Fractional Order

A common experimental workflow is to sweep $\alpha$ to find the optimal value for a given task:

double alphas[] = {0.3, 0.5, 0.7, 0.9, 1.0};
int n_alphas = sizeof(alphas) / sizeof(alphas[0]);
 
for (int i = 0; i < n_alphas; i++) {
    spires_flif_params params = {
        .alpha   = alphas[i],
        .tau_m   = 20.0,
        .v_rest  = -65.0,
        .v_th    = -50.0,
        .v_reset = -65.0,
        .bias    = 0.0,
        .L       = 200,
    };
    cfg.neuron_params = &params;
 
    spires_reservoir *r = NULL;
    spires_status s = spires_reservoir_create(&cfg, &r);
    if (s != SPIRES_OK) continue;
 
    /* Train, evaluate, record performance for alphas[i] */
 
    spires_reservoir_destroy(r);
}

Complexity Analysis

For a reservoir with $N$ neurons and history length $L$:

• Memory: $O(N \cdot L)$ for voltage histories, plus $O(L)$ for the shared coefficient array
• Time per step: $O(N \cdot L)$ for history sums, parallelized across neurons with OpenMP
• Initialization: $O(L)$ for coefficient precomputation (done once)

The $O(N \cdot L)$ per-step cost is the dominant expense for fractional reservoirs. The SPIRES implementation parallelizes the outer loop over neurons, so with $P$ threads the effective cost is $O((N/P) \cdot L)$.

Why GL is the Primary Fractional Method

The GL discretization is the default choice in SPIRES for several reasons:

1. Direct discretization: No integral approximation is needed. The GL formula naturally produces a discrete update rule.
2. Efficient coefficients: The recurrence relation avoids computing gamma functions at every step.
3. Circular buffer: The implementation is cache-friendly and avoids memory allocation during simulation.
4. Equivalence: For well-behaved initial conditions, the GL and Caputo discretizations produce identical results.
5. Parallelism: The per-neuron history sums are independent and parallelize trivially.

← Caputo Fractional LIF | Diffusive Fractional LIF →