arvoelke
commented
6 years ago I spent some time analyzing these different approaches, and arrived at the following argument for replacing nengo.AdaptiveLIF (I'll refer to this here as "ALIF") with the second approach from above that is said to be standard as implemented by Maass et al. (2018) (I'll refer to this here as "ALIFT" for 'ALIF Threshold'). Basically, what it comes down to is that ALIF does not improve the computational capabilities of the network compared to LIF, while ALIFT in fact does.
For the case of ALIF, all of the important bits can be understood from the following code:
def step_math(self, dt, J, output, voltage, ref, adaptation):
"""Implement the AdaptiveLIF nonlinearity."""
n = adaptation
LIF.step_math(self, dt, J - n, output, voltage, ref)
n += (dt / self.tau_n) * (self.inc_n * output - n)From this, we know that the adaptive state n is simply a filtered version of each spike train. Moreover, the filter is a lowpass with time-constant tau_n, scaled by inc_n, and discretized using Euler's method. To be concrete, consider the following LTI system:
alif_synapse = alif.inc_n * nengolib.signal.cont2discrete(
nengolib.Lowpass(alif.tau_n), dt=dt, method='euler') # ALIF code uses EulerBy the way, the filter should instead be using ZOH (the default method in Nengo), an issue I said I would raise in #975, but this difference is subtle.
Regardless, if we apply this filter to the ALIF neurons, then we get exactly the same vector as probing the adaptation state. Then, looking back at the ALIF code, since n is subtracted from the input current J, we can undo the effect of adaptation simply by adding back n to get (J - n) + n = J by linearity.
nengo.Connection(alif.neurons, alif.neurons, transform=1/alif.gain, synapse=alif_synapse)The observation is that an ALIF population is isomorphic to a LIF population with a recurrent diagonal connection. Conversely, any LIF population can be made ALIF by subtracting this same connection. Therefore, the ALIF model offers no more computational power than what's currently made available by LIF. The only "real" benefit of ALIF is that it abstracts away this extra connection, and perhaps makes it a bit faster to simulate.
There is some code at the bottom that empirically validates this equivalence, by adding this connection to convert a population of ALIF into LIF neurons with identical spike trains.
For the case of ALIFT, the adaptation has a divisive effect on the input current, rather than a subtractive effect. This model is implemented at the end of this post. The adaptive state is still exactly the same. The dynamic threshold becomes (1 + n), where n is the same as in ALIF (except my implementation uses ZOH instead of Euler discretization).
alift_synapse = alift.inc_adapt * nengolib.Lowpass(alift.tau_adapt)However, changes to the threshold do not alter the effective input current in a linear manner. Considering the block-diagram J(t) -> [h(t)] -> V(t) where J(t) is the input current, h(t) is the linear filter Lowpass(tau_rc) capturing all aspects of the LIF model apart from reset/threshold mechanics, and V(t) is the membrane voltage, then (1 + n)*J(t) is the input that produces (1 + n)*V(t) since scaling the input to an LTI is equivalent to scaling its output. This isn't true in an exact sense since n is time-varying, but it is pretty darn close (see below). In other words, a dynamic threshold of (1 + n) is approximately the same as dividing the input current by (1 + n). Conversely, we need to multiply J by (1 + n) to undo the adaptation and convert it into a LIF population. This multiplicative interaction is more than just a recurrent connection, since our synapses are linear.
The following demonstrates that, while both ALIF and ALIFT can be converted into a LIF population, the ALIF requires adding n to J with a simple recurrent connection, while the ALIFT requires multiplying J by (1 + n) with a more complicated node.
import numpy as np
import nengo
import nengolib
import matplotlib.pyplot as plt
import seaborn as sns
def install_tuning_curves(ens):
"""Use ens.seed to set ens.{encoders, max_rates, intercepts, gain, bias}."""
rng = np.random.RandomState(ens.seed)
ens.encoders = ens.encoders.sample(ens.n_neurons, ens.dimensions, rng=rng)
ens.max_rates = ens.max_rates.sample(ens.n_neurons, rng=rng)
ens.intercepts = ens.intercepts.sample(ens.n_neurons, rng=rng)
ens.gain, ens.bias = ens.neuron_type.gain_bias(
ens.max_rates, ens.intercepts)
from nengo.neurons import *
from nengo.builder.neurons import *
class AdaptiveLIFT(LIFRate):
probeable = ('spikes', 'voltage', 'refractory_time', 'threshold')
min_voltage = NumberParam('min_voltage', high=0)
tau_adapt = NumberParam('tau_adapt', low=0)
inc_adapt = NumberParam('inc_adapt', low=0)
def __init__(self, tau_rc=0.02, tau_ref=0.002, min_voltage=0,
amplitude=1, tau_adapt=0.1, inc_adapt=0.05):
super(AdaptiveLIFT, self).__init__(
tau_rc=tau_rc, tau_ref=tau_ref, amplitude=amplitude)
self.min_voltage = min_voltage
self.tau_adapt = tau_adapt
self.inc_adapt = inc_adapt
def step_math(self, dt, J, spiked, voltage, refractory_time, threshold):
# reduce all refractory times by dt
refractory_time -= dt
# compute effective dt for each neuron, based on remaining time.
# note that refractory times that have completed midway into this
# timestep will be given a partial timestep, and moreover these will
# be subtracted to zero at the next timestep (or reset by a spike)
delta_t = (dt - refractory_time).clip(0, dt)
# update voltage using discretized lowpass filter
# since v(t) = v(0) + (J - v(0))*(1 - exp(-t/tau)) assuming
# J is constant over the interval [t, t + dt)
voltage -= (J - voltage) * np.expm1(-delta_t / self.tau_rc)
# determine which neurons spiked (set them to 1/dt, else 0)
spiked_mask = voltage > threshold
spiked[:] = spiked_mask * (self.amplitude / dt)
# set v(0) = threshold and solve for t to compute the spike time
# TODO: not sure if this mask is the right way to handle log domain errors
threshold_spiked = threshold[spiked_mask]
m = (voltage[spiked_mask] - threshold_spiked) / (J[spiked_mask] - threshold_spiked)
t_spike = np.zeros_like(m)
t_spike[m < 1] = dt + self.tau_rc * np.log1p(-m[m < 1])
# update threshold using discretized lowpass filter
# applied to the input 1 + spiked * inc_adapt
threshold -= ((1 + self.inc_adapt * spiked - threshold) *
np.expm1(-dt / self.tau_adapt))
# set spiked voltages to zero, refractory times to tau_ref, and
# rectify negative voltages to a floor of min_voltage
voltage[voltage < self.min_voltage] = self.min_voltage
voltage[spiked_mask] = 0
refractory_time[spiked_mask] = self.tau_ref + t_spike
@Builder.register(AdaptiveLIFT)
def build_alift(model, lif, neurons):
model.sig[neurons]['voltage'] = Signal(
np.zeros(neurons.size_in), name="%s.voltage" % neurons)
model.sig[neurons]['refractory_time'] = Signal(
np.zeros(neurons.size_in), name="%s.refractory_time" % neurons)
model.sig[neurons]['threshold'] = Signal(
np.ones(neurons.size_in), name="%s.threshold" % neurons)
model.add_op(SimNeurons(
neurons=lif,
J=model.sig[neurons]['in'],
output=model.sig[neurons]['out'],
states=[model.sig[neurons]['voltage'],
model.sig[neurons]['refractory_time'],
model.sig[neurons]['threshold']]))
ens_kwargs = dict(
n_neurons=200,
dimensions=1,
seed=0
)
t = 1.0
dt = 0.001
tau_probe = 0.005
lif = nengo.LIF()
alif = nengo.AdaptiveLIF(tau_n=0.1, inc_n=0.1)
alift = AdaptiveLIFT(tau_adapt=0.1, inc_adapt=0.1)
with nengo.Network() as model:
u = nengo.Node(output=nengo.processes.WhiteSignal(t, high=10, rms=0.3, y0=0))
x_lif = nengo.Ensemble(neuron_type=lif, **ens_kwargs)
install_tuning_curves(x_lif)
x_alif = nengo.Ensemble(neuron_type=alif, **ens_kwargs)
install_tuning_curves(x_alif)
x_alif_undo = nengo.Ensemble(neuron_type=alif, **ens_kwargs)
install_tuning_curves(x_alif_undo)
x_alift = nengo.Ensemble(neuron_type=alift, **ens_kwargs)
install_tuning_curves(x_alift)
x_alift_undo = nengo.Ensemble(neuron_type=alift, **ens_kwargs)
install_tuning_curves(x_alift_undo)
nengo.Connection(u, x_lif, synapse=None)
nengo.Connection(u, x_alif, synapse=None)
nengo.Connection(u, x_alif_undo, synapse=None)
nengo.Connection(u, x_alift, synapse=None)
nengo.Connection(u, x_alift_undo, synapse=None)
# These two linear filters help explain/undo the two forms of
# adaptation in a Nengo-friendly language
alif_synapse = alif.inc_n * nengolib.signal.cont2discrete(
nengolib.Lowpass(alif.tau_n), dt=dt, method='euler') # ALIF code uses Euler
alift_synapse = alift.inc_adapt * nengolib.Lowpass(alift.tau_adapt)
p_u = nengo.Probe(u, synapse=tau_probe)
p_lif = nengo.Probe(x_lif, synapse=tau_probe)
p_alif = nengo.Probe(x_alif, synapse=tau_probe)
p_alif_undo = nengo.Probe(x_alif_undo, synapse=tau_probe)
p_alift = nengo.Probe(x_alift, synapse=tau_probe)
p_alift_undo = nengo.Probe(x_alift_undo, synapse=tau_probe)
# Extra probes to help explain ALIF and ALIFT respectively
p_adapt_expected = nengo.Probe(x_alif.neurons, 'spikes', synapse=alif_synapse)
p_adapt_actual = nengo.Probe(x_alif.neurons, 'adaptation', synapse=None)
p_threshold_expected = nengo.Probe(x_alift.neurons, 'spikes', synapse=alift_synapse)
p_threshold_actual = nengo.Probe(x_alift.neurons, 'threshold', synapse=None)
# Undo subtractive adaptation (ALIF)
nengo.Connection(x_alif_undo.neurons, x_alif_undo.neurons,
transform=1/x_alif_undo.gain, synapse=alif_synapse)
# Undo divisive adaptation (ALIFT)
J = nengo.Node(size_in=1+x_alift_undo.n_neurons,
output=lambda t, p: (
p[1:]*(p[0, None].dot(x_alift_undo.encoders.T) +
x_alift_undo.bias/x_alift_undo.gain)))
nengo.Connection(u, J[0], synapse=None)
nengo.Connection(x_alift_undo.neurons, J[1:], synapse=alift_synapse)
nengo.Connection(J, x_alift_undo.neurons, synapse=None)
with nengo.Simulator(model, dt=dt, progress_bar=False) as sim:
sim.run(t, progress_bar=False)
fig, ax = plt.subplots(2, 1, sharex=True, figsize=(14, 14))
ax[0].plot(sim.trange(), sim.data[p_alif_undo], alpha=0.8, label="ALIF Undo")
ax[0].plot(sim.trange(), sim.data[p_alif], alpha=0.8, label="ALIF")
ax[0].plot(sim.trange(), sim.data[p_u], lw=3, linestyle='--', zorder=3, label="Reference")
ax[0].legend()
ax[1].plot(sim.trange(), sim.data[p_alift_undo], alpha=0.8, label="ALIFT Undo")
ax[1].plot(sim.trange(), sim.data[p_alift], alpha=0.8, label="ALIFT")
ax[1].plot(sim.trange(), sim.data[p_u], lw=3, linestyle='--', zorder=3, label="Reference")
ax[1].legend()
ax[1].set_xlabel("Time (s)")
plt.show()
# Undoing the ALIF makes it identical to LIF (even at level of spikes)
assert np.allclose(sim.data[p_lif], sim.data[p_alif_undo])
# This is because the adaptive term is equivalent to filtering each spike train
assert np.allclose(sim.data[p_adapt_expected][1:], sim.data[p_adapt_actual][:-1])
# The threshold for ALIFT is also equivalent to filtering each spike train
# (however this has a divisive effect, which cannot be undone simply by addition)
assert np.allclose(1 + sim.data[p_threshold_expected][1:], sim.data[p_threshold_actual][:-1])I have also validated both conversions in the opposite direction (LIF -> ALIF and LIF -> ALIFT using subtraction and division, respectively) although it's not shown in this code.
The bottom three assert statements summarize some of the more important things going on here: (1) a recurrent additive connection converts an ALIF population into an identical LIF population, (2) which is because n is obtained by filtering each neuron with a lowpass, and (3) this also gives the dynamic threshold for ALIFT. However, changes to the threshold have a divisive effect on the input current that cannot be undone by a linear connection. In effect, ALIFT provides the network with new computational capabilities pertaining to multiplicative interactions between spike rates and input currents.
https://github.com/nengo/nengo/blob/d7e348a320e1c903dc5074bea0cb5bc3aca334d1/nengo/neurons.py#L542
Nengo references the
AdaptiveLIFmodel as coming from:However, I checked page 339 of this book, and could not find the approach used by Nengo:
Instead, I found two somewhat close alternatives:
The first approach multiplies the leak term by
(1 + Rg), wheregis the same as thenvariable in Nengo (andRis a constant). But this is not the same as replacing the input current withJ - nsince this misses the multiplicative interaction betweennandV.The second approach is a little imprecise, but I believe is basically saying to substitute the threshold (say 1) with
1 + n. This is the version that is said to be "standard", as used by Maass in his arXiv paper released this week (see equations 3,4 in supplementary). However, this is again not the same as replacing the input current withJ - n.