Accompanying Readings: Chapter 2
from IPython.display import YouTubeVideo
YouTubeVideo('KE952yueVLA', width=720, height=400, loop=1, autoplay=0)
from IPython.display import YouTubeVideo
YouTubeVideo('lfNVv0A8QvI', width=720, height=400, loop=1, autoplay=0)
Some sort of mapping between neural activity and a state in the world
Intuitively, we call this "representation"
Critical obvious difference between this and ANNs:
Encoding (nonlinear): $$ a_i = \begin{cases} 1 &\mbox{if } x \mod {2^{i}} > 2^{i-1} \\ 0 &\mbox{otherwise} \end{cases} $$
Decoding (linear): $$ \hat{x} = \sum_i a_i 2^{i-1} $$
Suppose: $x = 13$
Encoding: $a_1 = 1$, $a_2 = 0$, $a_3 = 1$, $a_4 = 1$
Decoding: $\hat{x} = 1*1+0*2+1*4+1*8 = 13$
Write decoder as $\hat{x} = \sum_ia_i d_i$
Linear decoding is nice and simple
The NEF uses linear decoding, but what about the encoding?
$a_i = f_i(x)$
Firing rate goes up as total input current goes up
What is $G_i$?
from IPython.display import YouTubeVideo
YouTubeVideo('hxdPdKbqm_I', width=720, height=400, loop=1, autoplay=0)
# Rectified linear neuron
%pylab inline
import numpy
import nengo
n = nengo.neurons.RectifiedLinear()
J = numpy.linspace(-1,10,100)
plot(J, n.rates(J, gain=30, bias=-25))
xlabel('J (current)')
ylabel('$a$ (Hz)');
Populating the interactive namespace from numpy and matplotlib
$ a = {1 \over {\tau_{ref}-\tau_{RC}ln(1-{1 \over J})}}$
#assume %pylab inline has been run
# Leaky integrate and fire
import numpy
import nengo
n = nengo.neurons.LIFRate(tau_rc=0.02, tau_ref=0.002) #n is a Nengo LIF neuron, these are defaults
J = numpy.linspace(-1,10,100)
plot(J, n.rates(J, gain=1, bias=-2))
xlabel('J (current)')
ylabel('$a$ (Hz)');
Neurons seem to be sensitive to particular values of $x$
What's the mapping between $x$ and $a$?
Sometimes they are fairly straight forward:
Note that the experimenters are graphing $a$, as a function of $x$
$x$ is the volume of the sound in this case
Any ideas?
#assume this has been run
#%pylab inline
import numpy
import nengo
n = nengo.neurons.LIFRate() #n is a Nengo LIF neuron
x = numpy.linspace(-100,0,100)
plot(x, n.rates(x, gain=1, bias=50), 'b') # x*1+50
plot(x, n.rates(x, gain=0.1, bias=10), 'r') # x*0.1+10
plot(x, n.rates(x, gain=0.5, bias=5), 'g') # x*0.05+5
plot(x, n.rates(x, gain=0.1, bias=4), 'c') #x*0.1+4))
xlabel('x')
ylabel('a');
For mapping #1, the NEF uses a linear map: $ J = \alpha x + J^{bias} $
$ J = - \alpha x + J^{bias} $
#assume this has been run
#%pylab inline
import numpy
import nengo
n = nengo.neurons.LIFRate()
e = numpy.array([1.0, 1.0])
e = e/numpy.linalg.norm(e)
a = numpy.linspace(-1,1,50)
b = numpy.linspace(-1,1,50)
X,Y = numpy.meshgrid(a, b)
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig = figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')
p = ax.plot_surface(X, Y, n.rates((X*e[0]+Y*e[1]), gain=1, bias=1.5),
linewidth=0, cstride=1, rstride=1, cmap=pylab.cm.jet)
import nengo
import numpy
n = nengo.neurons.LIFRate()
theta = numpy.linspace(0, 2*numpy.pi, 100)
x = numpy.array([numpy.cos(theta), numpy.sin(theta)])
plot(x[0],x[1])
axis('equal')
e = numpy.array([-1.0, 1.0])
e = e/numpy.linalg.norm(e)
plot([0,e[0]], [0,e[1]],'r')
gain = 1
bias = 2.5
figure()
plot(theta, n.rates(numpy.dot(x.T, e), gain=gain, bias=bias))
plot([numpy.arctan2(e[1],e[0])],0,'rv')
xlabel('angle')
ylabel('firing rate')
xlim(0, 2*numpy.pi);
Encoding
Decoding
The textbook uses $\phi$ for $d$ and $\tilde \phi$ for $e$
But where do we get $d_i$ from?
Find the optimal $d_i$
$ E = \frac{1}{2}\int_{-1}^1 (x-\hat{x})^2 \; dx $
$ \begin{align} E = \frac{1}{2}\int_{-1}^1 \left(x-\sum_i^N a_i d_i \right)^2 \; dx \end{align} $
$ \begin{align} {{\partial E} \over {\partial d_i}} &= {1 \over 2} \int_{-1}^1 2 \left[ x-\sum_j a_j d_j \right] (-a_i) \; dx \\ {{\partial E} \over {\partial d_i}} &= - \int_{-1}^1 a_i x \; dx + \int_{-1}^1 \sum_j a_j d_j a_i \; dx \end{align} $
$ \begin{align} \int_{-1}^1 a_i x \; dx &= \int_{-1}^1 \sum_j(a_j d_j a_i) \; dx \\ \int_{-1}^1 a_i x \; dx &= \sum_j \left(\int_{-1}^1 a_i a_j \; dx\right)d_j \end{align} $
$ \Upsilon = \Gamma d $
where
$ \begin{align} \Upsilon_i &= {1 \over 2} \int_{-1}^1 a_i x \;dx\\ \Gamma_{ij} &= {1 \over 2} \int_{-1}^1 a_i a_j \;dx \end{align} $
$ \begin{align} \sum_x a_i x / S &= \sum_j \left(\sum_x a_i a_j /S \right)d_j \\ \Upsilon &= \Gamma d \end{align} $
where
$ \begin{align} \Upsilon_i &= \sum_x a_i x / S \\ \Gamma_{ij} &= \sum_x a_i a_j / S \end{align} $
So given
$ \Upsilon = \Gamma d $
then
$ d = \Gamma^{-1} \Upsilon $
or, equivalently
$ d_i = \sum_j \Gamma^{-1}_{ij} \Upsilon_j $
import numpy
import nengo
from nengo.utils.ensemble import tuning_curves
from nengo.dists import Uniform
N = 10
model = nengo.Network(label='Neurons')
with model:
neurons = nengo.Ensemble(N, dimensions=1,
max_rates=Uniform(100,200)) #Defaults to LIF neurons,
#with random gains and biases for
#neurons between 100-200hz over -1,1
connection = nengo.Connection(neurons, neurons, #This is just to generate the decoders
solver=nengo.solvers.LstsqL2(reg=0)) #reg=0 means ignore noise
sim = nengo.Simulator(model)
d = sim.data[connection].weights.T
x, A = tuning_curves(neurons, sim)
xhat = numpy.dot(A, d)
pyplot.plot(x, A)
xlabel('x')
ylabel('firing rate (Hz)')
figure()
plot(x, x)
plot(x, xhat)
xlabel('$x$')
ylabel('$\hat{x}$')
ylim(-1, 1)
xlim(-1, 1)
figure()
plot(x, xhat-x)
xlabel('$x$')
ylabel('$\hat{x}-x$')
xlim(-1, 1)
print('RMSE', np.sqrt(np.average((x-xhat)**2)))
RMSE 0.04344424442937684
#Have to run previous python cell first
A_noisy = A + numpy.random.normal(scale=0.2*numpy.max(A), size=A.shape)
xhat = numpy.dot(A_noisy, d)
pyplot.plot(x, A_noisy)
xlabel('x')
ylabel('firing rate (Hz)')
figure()
plot(x, x)
plot(x, xhat)
xlabel('$x$')
ylabel('$\hat{x}$')
ylim(-1, 1)
xlim(-1, 1)
print('RMSE', np.sqrt(np.average((x-xhat)**2)))
RMSE 0.5075755918847665
Include noise while solving for decoders
$ \begin{align} \hat{x} &= \sum_i(a_i+\eta)d_i \\ E &= {1 \over 2} \int_{-1}^1 (x-\hat{x})^2 \;dx d\eta\\ &= {1 \over 2} \int_{-1}^1 \left(x-\sum_i(a_i+\eta)d_i\right)^2 \;dx d\eta\\ &= {1 \over 2} \int_{-1}^1 \left(x-\sum_i a_i d_i - \sum \eta d_i \right)^2 \;dx d\eta \end{align} $
$ \begin{align} E &= {1 \over 2} \int_{-1}^1 \left(x-\sum_i a_i d_i \right)^2 \;dx + \sum_{i,j} d_i d_j <\eta_i \eta_j>_\eta \\ &= {1 \over 2} \int_{-1}^1 \left(x-\sum_i a_i d_i \right)^2 \;dx + \sum_{i} d_i d_i <\eta_i \eta_i>_\eta \end{align} $
$ \begin{align} E = {1 \over 2} \int_{-1}^1 \left(x-\sum_i a_i d_i \right)^2 \;dx + \sigma^2 \sum_i d_i^2 \end{align} $
$ \begin{align} \Gamma_{ij} = \sum_x a_i a_j / S + \sigma^2 \delta_{ij} \end{align} $
Where $\delta_{ij}$ is the Kronecker delta: http://en.wikipedia.org/wiki/Kronecker_delta
To simplfy computing this using matrices, this can be written as $\Gamma=A^T A /S + \sigma^2 I$
import numpy
import nengo
from nengo.utils.ensemble import tuning_curves
from nengo.dists import Uniform
N = 100
model = nengo.Network(label='Neurons')
with model:
neurons = nengo.Ensemble(N, dimensions=1,
max_rates=Uniform(100,200)) #Defaults to LIF neurons,
#with random gains and biases for
#neurons between 100-200hz over -1,1
connection = nengo.Connection(neurons, neurons, #This is just to generate the decoders
solver=nengo.solvers.LstsqNoise(noise=0.2)) #Add noise ###NEW
sim = nengo.Simulator(model)
d = sim.data[connection].weights.T
x, A = tuning_curves(neurons, sim)
A_noisy = A + numpy.random.normal(scale=0.2*numpy.max(A), size=A.shape)
xhat = numpy.dot(A_noisy, d)
pyplot.plot(x, A_noisy)
xlabel('x')
ylabel('firing rate (Hz)')
figure()
plot(x, x)
plot(x, xhat)
xlabel('$x$')
ylabel('$\hat{x}$')
ylim(-1, 1)
xlim(-1, 1)
print('RMSE', np.sqrt(np.average((x-xhat)**2)))
RMSE 0.059589387758251525
$ \begin{align} E = {1 \over 2} \int_{-1}^1 \left(x-\sum_i a_i d_i \right)^2 \;dx + \sigma^2 \sum_i d_i^2 \end{align} $
$ \begin{align} E_{distortion} = {1 \over 2} \int_{-1}^1 \left(x-\sum_i a_i d_i \right)^2 dx \end{align} $
$ \begin{align} E_{noise} = \sigma^2 \sum_i d_i^2 \end{align} $
$ E = {1 \over 2} \int_{-1}^1 (x-\hat{x})^2 dx $
So that's actually a squared error term
Also, as number of neurons is greater than 100 or so, the error is dominated by the noise term ($1/N$).
From http://www.nature.com/nrn/journal/v3/n12/full/nrn986.html
There are also neurons whose response goes the other way. All of the neurons are directly connected to the muscle controlling the horizontal direction of the eye, and that's the only thing that muscle does, so we're pretty sure this is what's being repreesnted.
System Description
Design Specification
Implementation
#%pylab inline
import numpy
import nengo
from nengo.utils.ensemble import tuning_curves
from nengo.dists import Uniform
N = 10
tau_rc = 20
tau_ref = .001
lif_model = nengo.LIFRate(tau_rc=tau_rc, tau_ref=tau_ref)
model = nengo.Network(label='Neurons')
with model:
neurons = nengo.Ensemble(N, dimensions=1,
max_rates = Uniform(250,300),
neuron_type = lif_model)
sim = nengo.Simulator(model)
x, A = tuning_curves(neurons, sim)
plot(x, A)
xlabel('x')
ylabel('firing rate (Hz)');
#Have to run previous code cell first
noise = 0.2
with model:
connection = nengo.Connection(neurons, neurons, #This is just to generate the decoders
solver=nengo.solvers.LstsqNoise(noise=0.2)) #Add noise ###NEW
sim = nengo.Simulator(model)
d = sim.data[connection].weights.T
x, A = tuning_curves(neurons, sim)
A_noisy = A + numpy.random.normal(scale=noise*numpy.max(A), size=A.shape)
xhat = numpy.dot(A_noisy, d)
print('RMSE with %d neurons is %g'%(N, np.sqrt(np.average((x-xhat)**2))))
figure()
plot(x, x)
plot(x, xhat)
xlabel('$x$')
ylabel('$\hat{x}$')
ylim(-1, 1)
xlim(-1, 1);
RMSE with 10 neurons is 0.151211
System Description
Design Specification
Implementation
import numpy
import nengo
n = nengo.neurons.LIFRate()
theta = numpy.linspace(-numpy.pi, numpy.pi, 100)
x = numpy.array([numpy.sin(theta), numpy.cos(theta)])
e = numpy.array([-1.0, 0])
plot(theta*180/numpy.pi, n.rates(numpy.dot(x.T, e), bias=1., gain=0.2)) #bias 1->1.5
xlabel('angle')
ylabel('firing rate')
xlim(-180, 180)
show()
Does it match empirical data?
Interestingly, Georgopoulos suggests doing linear decoding: