Skip to content

Reservoir Computing — Auto-Critical LSM

Liquid State Machine with mean-field auto-criticality tuning. Zero hyperparameter tuning: the critical weight W_c = θ / (2βNp) is computed analytically from threshold (θ), leak (β), neuron count (N), and connectivity (p).

Theory

At the critical regime, exactly half the neurons fire in each refractory period, maximizing the reservoir's computational capacity (kernel quality). The reservoir projects inputs into a high-dimensional nonlinear state space; only the readout layer (ridge regression) is trained.

The critical weight formula derives from mean-field analysis: the expected number of post-synaptic spikes per pre-synaptic spike equals 1.0 at criticality, producing a branching ratio of exactly 1.

Components

  • AutoCriticalReservoir — Full LSM pipeline with auto-criticality tuning.
Parameter Default Meaning
n_inputs (required) Input dimension
n_neurons 1000 Reservoir size
n_outputs 10 Readout dimension
threshold 1.0 LIF spike threshold
leak 0.1 Membrane leak factor (0-1)
connectivity 0.1 Synapse existence probability
seed 42 RNG seed

Key methods:

  • step(x) — Process one timestep, return spike vector
  • run(inputs) — Run full input sequence (T, n_inputs) → state matrix (T, n_neurons)
  • fit_readout(states, targets, ridge) — Train readout via ridge regression
  • predict(states) — Predict from reservoir states
  • train_and_predict(train_in, train_tgt, test_in) — Full pipeline

  • metrics(inputs) → ReservoirMetrics — firing fraction, criticality error, kernel quality, spectral radius

  • ReservoirMetrics — Quality metrics dataclass with summary() method.

Usage

from sc_neurocore.reservoir import AutoCriticalReservoir

# Zero-config reservoir
res = AutoCriticalReservoir(n_inputs=2, n_neurons=500)
print(f"Spectral radius: {res.spectral_radius:.3f}")
print(f"Critical weight: {res.w_critical:.6f}")

# Run and train
import numpy as np
train_in = np.random.randn(200, 2)
train_tgt = np.sin(np.arange(200)).reshape(-1, 1)
test_in = np.random.randn(50, 2)

preds = res.train_and_predict(train_in, train_tgt, test_in)

# Quality assessment
metrics = res.metrics(train_in)
print(metrics.summary())

Reference: Scientific Reports 2025 — mean-field analytical framework for configuring spiking reservoirs at the critical regime.

See Tutorial 74: Reservoir Computing.

sc_neurocore.reservoir

Liquid State Machine with mean-field auto-criticality tuning.

AutoCriticalReservoir

Spiking Liquid State Machine with automatic criticality tuning.

Parameters

n_inputs : int n_neurons : int Reservoir size. n_outputs : int Readout dimension. threshold : float LIF spike threshold. leak : float Membrane leak factor (0-1). Higher = faster decay. connectivity : float Fraction of possible synapses that exist (sparsity). seed : int

Source code in src/sc_neurocore/reservoir/auto_reservoir.py 
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
class AutoCriticalReservoir:
    """Spiking Liquid State Machine with automatic criticality tuning.

    Parameters
    ----------
    n_inputs : int
    n_neurons : int
        Reservoir size.
    n_outputs : int
        Readout dimension.
    threshold : float
        LIF spike threshold.
    leak : float
        Membrane leak factor (0-1). Higher = faster decay.
    connectivity : float
        Fraction of possible synapses that exist (sparsity).
    seed : int
    """

    def __init__(
        self,
        n_inputs: int,
        n_neurons: int = 1000,
        n_outputs: int = 10,
        threshold: float = 1.0,
        leak: float = 0.1,
        connectivity: float = 0.1,
        seed: int = 42,
    ):
        self.n_inputs = n_inputs
        self.n_neurons = n_neurons
        self.n_outputs = n_outputs
        self.threshold = threshold
        self.leak = leak
        self.connectivity = connectivity

        rng = np.random.RandomState(seed)

        # Mean-field critical weight: W_c = theta / (2 * beta * N * p)
        effective_n = n_neurons * connectivity
        self.w_critical = threshold / max(2.0 * leak * effective_n, 1e-8)

        # Reservoir weights: sparse, scaled to W_critical
        mask = rng.random((n_neurons, n_neurons)) < connectivity
        np.fill_diagonal(mask, False)
        self.W_res = rng.randn(n_neurons, n_neurons) * self.w_critical
        self.W_res *= mask

        # Input weights
        self.W_in = rng.randn(n_neurons, n_inputs) * np.sqrt(2.0 / n_inputs)

        # Readout weights (trained by ridge regression)
        self.W_out = np.zeros((n_outputs, n_neurons))

        # State
        self._v = np.zeros(n_neurons)
        self._spikes = np.zeros(n_neurons)

    @property
    def spectral_radius(self) -> float:
        eigvals = np.abs(np.linalg.eigvals(self.W_res))
        return float(eigvals.max()) if len(eigvals) > 0 else 0.0

    def reset(self):
        self._v = np.zeros(self.n_neurons)
        self._spikes = np.zeros(self.n_neurons)

    def step(self, x: np.ndarray) -> np.ndarray:
        """Process one timestep, return reservoir state (spikes)."""
        current = self.W_in @ x + self.W_res @ self._spikes
        self._v = (1 - self.leak) * self._v + self.leak * current
        self._spikes = (self._v >= self.threshold).astype(np.float64)
        self._v -= self._spikes * self.threshold
        return self._spikes.copy()

    def run(self, inputs: np.ndarray) -> np.ndarray:
        """Run input sequence through reservoir, return state matrix.

        Parameters
        ----------
        inputs : ndarray of shape (T, n_inputs)

        Returns
        -------
        ndarray of shape (T, n_neurons)
        """
        self.reset()
        T = inputs.shape[0]
        states = np.zeros((T, self.n_neurons))
        for t in range(T):
            states[t] = self.step(inputs[t])
        return states

    def fit_readout(self, states: np.ndarray, targets: np.ndarray, ridge: float = 1e-4):
        """Train readout via ridge regression.

        Parameters
        ----------
        states : ndarray of shape (T, n_neurons)
        targets : ndarray of shape (T, n_outputs)
        ridge : float
            Regularization strength.
        """
        # W_out = targets^T @ states @ (states^T @ states + ridge*I)^{-1}
        S = states
        reg = ridge * np.eye(self.n_neurons)
        self.W_out = np.linalg.solve(S.T @ S + reg, S.T @ targets).T

    def predict(self, states: np.ndarray) -> np.ndarray:
        """Predict from reservoir states."""
        return states @ self.W_out.T

    def train_and_predict(
        self, train_inputs: np.ndarray, train_targets: np.ndarray, test_inputs: np.ndarray
    ) -> np.ndarray:
        """Full pipeline: run train, fit readout, run test, predict."""
        train_states = self.run(train_inputs)
        self.fit_readout(train_states, train_targets)
        test_states = self.run(test_inputs)
        return self.predict(test_states)

    def metrics(self, inputs: np.ndarray) -> ReservoirMetrics:
        """Compute reservoir quality metrics."""
        states = self.run(inputs)
        firing_fraction = float(states.mean())
        criticality_error = abs(firing_fraction - 0.5)

        # Kernel quality: rank of state matrix normalized by timesteps
        rank = np.linalg.matrix_rank(states)
        kernel_quality = rank / max(states.shape[0], 1)

        return ReservoirMetrics(
            firing_fraction=firing_fraction,
            criticality_error=criticality_error,
            kernel_quality=kernel_quality,
            spectral_radius=self.spectral_radius,
        )

step(x)

Process one timestep, return reservoir state (spikes).

Source code in src/sc_neurocore/reservoir/auto_reservoir.py 
112
113
114
115
116
117
118
def step(self, x: np.ndarray) -> np.ndarray:
    """Process one timestep, return reservoir state (spikes)."""
    current = self.W_in @ x + self.W_res @ self._spikes
    self._v = (1 - self.leak) * self._v + self.leak * current
    self._spikes = (self._v >= self.threshold).astype(np.float64)
    self._v -= self._spikes * self.threshold
    return self._spikes.copy()

run(inputs)

Run input sequence through reservoir, return state matrix.

Parameters

inputs : ndarray of shape (T, n_inputs)

Returns

ndarray of shape (T, n_neurons)

Source code in src/sc_neurocore/reservoir/auto_reservoir.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def run(self, inputs: np.ndarray) -> np.ndarray:
    """Run input sequence through reservoir, return state matrix.

    Parameters
    ----------
    inputs : ndarray of shape (T, n_inputs)

    Returns
    -------
    ndarray of shape (T, n_neurons)
    """
    self.reset()
    T = inputs.shape[0]
    states = np.zeros((T, self.n_neurons))
    for t in range(T):
        states[t] = self.step(inputs[t])
    return states

fit_readout(states, targets, ridge=0.0001)

Train readout via ridge regression.

Parameters

states : ndarray of shape (T, n_neurons) targets : ndarray of shape (T, n_outputs) ridge : float Regularization strength.

Source code in src/sc_neurocore/reservoir/auto_reservoir.py 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
def fit_readout(self, states: np.ndarray, targets: np.ndarray, ridge: float = 1e-4):
    """Train readout via ridge regression.

    Parameters
    ----------
    states : ndarray of shape (T, n_neurons)
    targets : ndarray of shape (T, n_outputs)
    ridge : float
        Regularization strength.
    """
    # W_out = targets^T @ states @ (states^T @ states + ridge*I)^{-1}
    S = states
    reg = ridge * np.eye(self.n_neurons)
    self.W_out = np.linalg.solve(S.T @ S + reg, S.T @ targets).T

predict(states)

Predict from reservoir states.

Source code in src/sc_neurocore/reservoir/auto_reservoir.py 
153
154
155
def predict(self, states: np.ndarray) -> np.ndarray:
    """Predict from reservoir states."""
    return states @ self.W_out.T

train_and_predict(train_inputs, train_targets, test_inputs)

Full pipeline: run train, fit readout, run test, predict.

Source code in src/sc_neurocore/reservoir/auto_reservoir.py 
157
158
159
160
161
162
163
164
def train_and_predict(
    self, train_inputs: np.ndarray, train_targets: np.ndarray, test_inputs: np.ndarray
) -> np.ndarray:
    """Full pipeline: run train, fit readout, run test, predict."""
    train_states = self.run(train_inputs)
    self.fit_readout(train_states, train_targets)
    test_states = self.run(test_inputs)
    return self.predict(test_states)

metrics(inputs)

Compute reservoir quality metrics.

Source code in src/sc_neurocore/reservoir/auto_reservoir.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
def metrics(self, inputs: np.ndarray) -> ReservoirMetrics:
    """Compute reservoir quality metrics."""
    states = self.run(inputs)
    firing_fraction = float(states.mean())
    criticality_error = abs(firing_fraction - 0.5)

    # Kernel quality: rank of state matrix normalized by timesteps
    rank = np.linalg.matrix_rank(states)
    kernel_quality = rank / max(states.shape[0], 1)

    return ReservoirMetrics(
        firing_fraction=firing_fraction,
        criticality_error=criticality_error,
        kernel_quality=kernel_quality,
        spectral_radius=self.spectral_radius,
    )

ReservoirMetrics dataclass

Reservoir quality metrics.

Source code in src/sc_neurocore/reservoir/auto_reservoir.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
@dataclass
class ReservoirMetrics:
    """Reservoir quality metrics."""

    firing_fraction: float  # fraction of neurons active
    criticality_error: float  # |firing_fraction - 0.5|
    kernel_quality: float  # linear separability of reservoir states
    spectral_radius: float

    def summary(self) -> str:
        return (
            f"Reservoir: firing={self.firing_fraction:.3f}, "
            f"criticality_err={self.criticality_error:.4f}, "
            f"kernel_q={self.kernel_quality:.3f}, "
            f"spectral_r={self.spectral_radius:.3f}"
        )