In this exercise, we seek to construct a Hopfield network with $N = 64$ neurons, whose dynamics are given by the differential equation
$$\dot{\textbf{x}} = - \textbf{x} + f(W \textbf{x}) + σ η(t)$$
where
$\textbf{x}$ is the vector of firing rates
$W$ is the $N × N$ weight matrix
and $σ = 0.1$ is an external noise level.
For mathematical simplicity, we will set the activation function to
$$f(x) = \sign(x)$$
so that neural activity can run from $−1$ to $1$ (this step is simply for mathematical convenience, i.e. it makes the math simpler. If necessary, however, one can also simulate the Hopfield network with a more realistic input-output functions such as the one we used above.)
To store a pattern $\textbf{p}$ in the Hopfield network, we set the weight matrix to
$$W_{ij} = p_i p_j$$
i.e. the weight matrix is the outer product of two vectors:
$$W = \frac 1 N \textbf{p} \textbf{p}^\T$$
where by convention:
- $\textbf{p}$ is assumed to be a column vector
- the superscript $\T$ denotes the matrix transpose; here, it makes a column vector into a row vector.
In [22]:
# /!\ I'm using Python 3
import plotly
import numpy as np
import plotly.plotly as py
import plotly.graph_objs as go
import matplotlib.pyplot as plt
import scipy.io as sio
from scipy import stats
%matplotlib inline
plotly.offline.init_notebook_mode(connected=True)
from IPython.core.display import display, HTML, Markdown
# The polling here is to ensure that plotly.js has already been loaded before
# setting display alignment in order to avoid a race condition.
display(HTML(
'<script>'
'var waitForPlotly = setInterval( function() {'
'if( typeof(window.Plotly) !== "undefined" ){'
'MathJax.Hub.Config({ SVG: { font: "STIX-Web" }, displayAlign: "center" });'
'MathJax.Hub.Queue(["setRenderer", MathJax.Hub, "SVG"]);'
'clearInterval(waitForPlotly);'
'}}, 250 );'
'</script>'
))
# Colorscales
def colorscale_list(cmap, number_colors, return_rgb_only=False):
cm = plt.get_cmap(cmap)
colors = [np.array(cm(i/number_colors)) for i in range(1, number_colors+1)]
rgb_colors_plotly = []
rgb_colors_only = []
for i, c in enumerate(colors):
col = 'rgb{}'.format(tuple(255*c[:-1]))
rgb_colors_only.append(col)
rgb_colors_plotly.append([i/number_colors, col])
rgb_colors_plotly.append([(i+1)/number_colors, col])
return rgb_colors_only if return_rgb_only else rgb_colors_plotly
import urllib
(a) Invent a pattern $\textbf{p}$! Each component $p_i$ of this pattern should take either the value $-1$ or $1$. One way to visualize the “pattern” $\textbf{p}$ is to map this vector onto a matrix – for a network of $N = 64$ neurons, you can for instance map the $64$-dimensional vector $\textbf{p}$ onto a matrix of size $[8 × 8]$.¶
In [23]:
N = 64
sigma = 0.1
delta_t = 0.1
W = {}
ps_list = {}
In [24]:
# Display several images
def display_images(image_list, title='Step', cols=16, delta_step=None):
if delta_step is None:
delta_step = len(image_list)//100+1
fig = plt.figure()
rows = len(image_list[::delta_step])//cols+1
for i, img in enumerate(image_list[::delta_step]):
plt.subplot(rows+1, cols, i+1)
plt.imshow(img.reshape([8,-1]))
plt.axis('off')
if title is not None:
plt.title('{} {}'.format(title, i*delta_step+1))
#print((np.array(fig.get_size_inches())-np.array([0, 1]))*len(image_list[::delta_step]))
fig.set_size_inches((17,1.8*rows))
#plt.tight_layout(w_pad=0, h_pad=0)
plt.show()
In [25]:
# Retrieving the patterns
image_names = ["battery.jpg", "computer.jpg", "coffee.jpg", "thumbsup.jpg", "toolbox.jpg", "copy.jpg",\
"RSS.jpg", "phone.jpg", "arrow.jpg", "signal.jpg"]
thresholds = [100, 100, 180, 130, 100, 100, 180, 130, 100, 100]
images = []
patterns = []
for i, img in enumerate(image_names):
image = np.array(plt.imread(urllib.request.urlopen("http://younesse.net/images/Neuromodeling/"+img), 0)).dot(np.ones((3, 1))).reshape([8, -1])
image = 2*(np.array(image)>thresholds[i]).astype(float)-1
images.append(image)
patterns.append(image.flatten())
display(Markdown("### Available patterns"))
rows, cols = 4, 4
plt.figure(figsize=(15,15))
for k, img in enumerate(images):
plt.subplot(rows+1, cols, k+1)
plt.imshow(img)
plt.title('Pattern {}: {}'.format(k+1, image_names[k].replace(".jpg", "").capitalize()))
plt.axis('off')
plt.tight_layout(w_pad=0, h_pad=0)
plt.show()
(b) Create the weight matrix $W$ as described above and simulate the network with stepwidth $dt = 0.1$, using random initial conditions. To visualize what is going on, try to plot the network activity at each time step using the false-color-plot from above. Note that you will have to remap the $64$-dimensional network activity $\textbf{x}(t)$ at each time step onto the $8 × 8$ matrix! Repeat the simulation several times. What fixed point(s) does the network relax to?¶
In [5]:
nb_patterns = 1
W[nb_patterns] = sum(np.outer(patterns[i], patterns[i]) for i in range(nb_patterns))/N
display(Markdown("## Storing {} pattern:".format(nb_patterns)))
plt.figure(figsize=(3,3*nb_patterns))
for i in range(nb_patterns):
plt.subplot(1, nb_patterns, i+1)
plt.imshow(patterns[i].reshape([8, -1]))
plt.axis('off')
plt.tight_layout(w_pad=0, h_pad=0)
plt.show()
In [6]:
def network_simulation(p_0=None, delta_t=delta_t, W=W[nb_patterns]):
if p_0 is None:
p_0 = 2*np.random.randint(2, size=N).astype("float")-1
ps = [p_0]
sigma_sqrt_t = sigma*np.sqrt(delta_t)
converged = False
while not converged:
ps.append((1-delta_t)*ps[-1] + delta_t*np.sign(W.dot(ps[-1])) + sigma_sqrt_t*np.random.standard_normal(N))
if np.all(ps[-1] - ps[-2] < .05):
converged = True
return ps
In [7]:
ps_list[nb_patterns] = []
np.random.seed(nb_patterns)
for _ in range(5):
ps_list[nb_patterns].append(network_simulation(W=W[nb_patterns]))
In [8]:
display(Markdown("### Network simulations"))
for ps in (ps_list[nb_patterns]):
display_images(ps)
display(Markdown("__________________"))
(c) Try to create a second pattern $\textbf{q}$ and design the weight matrix according to:¶
$$W = \frac 1 N \left(\textbf{p}\textbf{p}^\T + \textbf{q}\textbf{q}^\T\right)$$
Repeat the experiment from (c). Verify that the network retrieves the right pattern by starting from initial conditions that are close to that pattern.¶
In [9]:
nb_patterns = 2
W[nb_patterns] = sum(np.outer(patterns[i], patterns[i]) for i in range(nb_patterns))/N
display(Markdown("## Storing {} patterns:".format(nb_patterns)))
plt.figure(figsize=(3*nb_patterns,3*nb_patterns))
for i in range(nb_patterns):
plt.subplot(1, nb_patterns, i+1)
plt.imshow(patterns[i].reshape([8, -1]))
plt.axis('off')
plt.tight_layout(w_pad=0, h_pad=0)
plt.show()
In [10]:
ps_list[nb_patterns] = []
np.random.seed(nb_patterns)
for _ in range(5):
ps_list[nb_patterns].append(network_simulation(W=W[nb_patterns]))
In [11]:
display(Markdown("### Network simulations"))
for ps in ps_list[nb_patterns]:
display_images(ps)
display(Markdown("__________________"))
(d) The general rule is that for a set of $M$ patterns $\textbf{p}_i$, the weight matrix is constructed as:¶
$$W = \frac 1 N \sum_{i=1}^M \textbf{p}_i \textbf{p}_i^\T$$
If you construct more and more patterns, how many patterns can you store in your network with $N = 64$ neurons? (More specifically, if you set the initial conditions close to one of the patterns, does your network relax back to that pattern or not?)¶
Investigate what happens if you replace the sign-function with a $\tanh$-function.¶
In [12]:
nb_patterns = 3
W[nb_patterns] = sum(np.outer(patterns[i], patterns[i]) for i in range(nb_patterns))/N
display(Markdown("## Storing {} patterns:".format(nb_patterns)))
plt.figure(figsize=(3*nb_patterns,3*nb_patterns))
for i in range(nb_patterns):
plt.subplot(1, nb_patterns, i+1)
plt.imshow(patterns[i].reshape([8, -1]))
plt.axis('off')
plt.tight_layout(w_pad=0, h_pad=0)
plt.show()
In [13]:
ps_list[nb_patterns] = []
np.random.seed(nb_patterns)
for _ in range(5):
ps_list[nb_patterns].append(network_simulation(W=W[nb_patterns]))
In [14]:
display(Markdown("### Network simulation"))
p_0 = np.array([-1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1., 1., -1.,
1., -1., -1., -1., -1., -1., -1., -1., 1., 1., -1., 1., -1.,
-1., 1., 1., -1., -1., 1., -1., 1., -1., 1., 1., -1., 1.,
-1., 1., -1., 1., 1., 1., 1., -1., -1., 1., 1., -1., -1.,
-1., -1., -1., -1., -1., 1., 1., -1., 1., -1., 1., 1.])
ps_list[nb_patterns].append(network_simulation(p_0=p_0, W=W[nb_patterns]))
display(Markdown("### Network simulations"))
for ps in ps_list[nb_patterns]:
display_images(ps)
display(Markdown("__________________"))
In [15]:
plt.imshow((patterns[0] - patterns[1] - patterns[2]).reshape([8,-1]))
plt.axis('off')
plt.show()
plt.imshow((patterns[0] - patterns[1] + patterns[2]).reshape([8,-1]))
plt.axis('off')
plt.show()
In [19]:
for nb_patterns in range(1, 11):
display(Markdown("## Storing {} patterns:".format(nb_patterns)))
np.random.seed(nb_patterns)
W[nb_patterns] = sum(np.outer(patterns[i], patterns[i]) for i in range(nb_patterns))/N
plt.figure(figsize=(3*nb_patterns,3*nb_patterns))
np.random.seed(nb_patterns)
for i in range(nb_patterns):
plt.subplot(1, nb_patterns, i+1)
plt.imshow(patterns[i].reshape([8, -1]))
plt.axis('off')
plt.tight_layout(w_pad=0, h_pad=0)
plt.show()
display(Markdown("### Network simulation starting from each pattern"))
for patter in patterns[:nb_patterns]:
ps = network_simulation(p_0=patter, W=W[nb_patterns])
delta_step2 = len(ps)//100+2
display_images(ps, delta_step=delta_step2)
display(Markdown("__________________"))
