# PROBLEM 3: Hopfield Model¶

$$\newcommand{\T}{ {\raise{0.7ex}{\intercal}}} \newcommand{\sign}{\textrm{sign}}$$

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.show()

In [25]:
# Retrieving the patterns

image_names = ["battery.jpg", "computer.jpg", "coffee.jpg", "thumbsup.jpg", "toolbox.jpg", "copy.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.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.show()


## Storing 1 pattern:¶

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.show()


## Storing 2 patterns:¶

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$$

## 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.show()


## Storing 3 patterns:¶

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("__________________"))