Skip to content

Tutorial 74: Auto-Critical Reservoir Computing

Zero-hyperparameter Liquid State Machine with mean-field auto-criticality. The reservoir self-tunes to edge-of-chaos dynamics for maximum computational capacity. Only the readout layer needs training.

AutoCriticalReservoir

import numpy as np
from sc_neurocore.reservoir import AutoCriticalReservoir

# Create reservoir: 64 inputs, 1000 recurrent neurons, 10 outputs
reservoir = AutoCriticalReservoir(
    n_inputs=64,
    n_neurons=1000,
    n_outputs=10,
)

# Train readout on data (reservoir weights are fixed)
train_x = np.random.randn(500, 64)
train_y = np.eye(10)[np.random.randint(0, 10, 500)]
test_x = np.random.randn(100, 64)

predictions = reservoir.train_and_predict(train_x, train_y, test_x)

# Criticality metrics
metrics = reservoir.metrics(test_x)
print(metrics.summary())

Why Auto-Criticality

Standard reservoirs require manual tuning of spectral radius, input scaling, and leak rate. The auto-critical reservoir adjusts these automatically using mean-field theory to maintain the edge of chaos — the regime where computational capacity is maximized.

API Reference

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}"
        )