Learning weights in a Boltzmann machine

I'm trying to understand how Boltzmann machines work, but I'm not quite sure how weights are learned, and haven't been able to find a clear description. Is the following correct? (Also, pointers to any good Boltzmann machine explanations would also be great.)

We have a set of visible units (e.g., corresponding to black/white pixels in an image) and a set of hidden units. Weights are initialized somehow (e.g., uniformly from [-0.5, 0.5]), and then we alternate between the following two phases until some stopping rule is reached:

1. Clamped phase - In this phase, all the values of the visible units are fixed, so we only update the states of the hidden units (according to the Boltzmann stochastic activation rule). We update until the network has reached equilibrium. Once we reach equilibrium, we continue updating $N$ more times (for some predefined $N$), keeping track of the average of $x_i x_j$ (where $x_i, x_j$ are the states of nodes $i$ and $j$). After those $N$ equilibrium updates, we update $w_ij = w_ij + \frac{1}{C} Average(x_i x_j)$, where $C$ is some learning rate. (Or, instead of doing a batch update at the end, do we update after we equilibrium step?)

2. Free phase - In this phase, the states of all units are updated. Once we reach equilibrium, we similarly continue updating N' more times, but instead of adding correlations at the end, we subtract: $w_{ij} = w_{ij} - \frac{1}{C} Average(x_i x_j)$.

So my main questions are:

1. Whenever we're in the clamped phase, do we reset the visible units to one of the patterns we want to learn (with some frequency that represents the importance of that pattern), or do we leave the visible units in the state they were in at the end of the free phase?

2. Do we do a batch update of the weights at the end of each phase, or update the weights at each equilibrium step in the phase? (Or, is either one fine?)

Intuitively, you can think of visible units as "what the model sees" and hidden units as "model's state of mind". When you set all visible units to some values, you "show the data to the model". Then, when you activate hidden units, the model adjusts it's state of mind to what it sees.

Next you let the model go free and fantasize. It will become shut-in and literally see some things it's mind generates, and generate new states of mind based on those images.

What we do by adjusting the weights (and biases) is making the model believe more in the data and less in it's own fantasies. This way after some training it will believe in some (hopefully) pretty good model of data, and we can for example ask "do you believe in this pair (X,Y)? How likely do you find it? What's your opinion mr. Boltzmann Machine?"

Finally here is a brief description of Energy Based Models, which should give you some intuition where do Clamped and Free phases come from and how we want to run them.

http://deeplearning.net/tutorial/rbm.html#energy-based-models-ebm

It is funny to see that the intuitively clear update rules come out form derivation of log-likelihood of generating data by the model.

With these intuitions in mind it's now easier to answer your questions:

1. We have to reset the visible units to some data we would like the model to believe in. If we use the values from the end of free phase it will just continue fantasizing, end enforce it's own misguided beliefs.

2. It's better to do updates after the end of the phase. Especially if it's the clamped phase, it's better to give the model some time to "focus" on the data. Earlier updates will slow down convergence, as they enforce the connections when the model hasn't adjusted it's state of mind to reality yet. Updating the weight's after each equilibrium step while fantasizing should be less harmful, although I have no experience with that.

If you want to improve your intuition on EBM, BM, and RBM I'd advise watching some of Geoffrey Hinton's lectures on the subject, he has some good analogies.

1. Yes, "we reset (clamp) the visible units to one of the patterns we want to learn (with some frequency that represents the importance of that pattern)."

2. Yes, "we do a batch update of the weights at the end of each phase." I don't think that updating "the weights at each equilibrium step in the phase" will lead to fast convergence because the network "gets distracted" by instantaneous errors — I have implemented Boltzmann machines that way, and I remember it not working very well until I changed it to a batch update.

Here is sample Python code for Boltzmann Machines based on Paul Ivanov's code from

http://redwood.berkeley.edu/wiki/VS265:_Homework_assignments

import numpy as np

def extract_patches(im,SZ,n):
imsize,imsize=im.shape;
X=np.zeros((n,SZ**2),dtype=np.int8);
startsx= np.random.randint(imsize-SZ,size=n)
startsy=np.random.randint(imsize-SZ,size=n)
for i,stx,sty in zip(xrange(n), startsx,startsy):
P=im[sty:sty+SZ, stx:stx+SZ];
X[i]=2*P.flat[:]-1;
return X.T

def sample(T,b,n,num_init_samples):
"""
sample.m - sample states from model distribution

function S = sample(T,b,n, num_init_samples)

T:                weight matrix
b:                bias
n:                number of samples
num_init_samples: number of initial Gibbs sweeps
"""
N=T.shape

# initialize state vector for sampling
s=2*(np.random.rand(N)<sigmoid(b))-1

for k in xrange(num_init_samples):
s=draw(s,T,b)

# sample states
S=np.zeros((N,n))
S[:,0]=s
for i in xrange(1,n):
S[:,i]=draw(S[:,i-1],T,b)

return S

def sigmoid(u):
"""
sigmoid.m - sigmoid function

function s = sigmoid(u)
"""
return 1./(1.+np.exp(-u));

def draw(Sin,T,b):
"""
draw.m - perform single Gibbs sweep to draw a sample from distribution

function S = draw(Sin,T,b)

Sin:      initial state
T:        weight matrix
b:        bias
"""
N=Sin.shape
S=Sin.copy()
rand = np.random.rand(N,1)
for i in xrange(N):
h=np.dot(T[i,:],S)+b[i];
S[i]=2*(rand[i]<sigmoid(h))-1;

return S

def run(im, T=None, b=None, display=True,N=4,num_trials=100,batch_size=100,num_init_samples=10,eta=0.1):
SZ=np.sqrt(N);
if T is None: T=np.zeros((N,N)); # weight matrix
if b is None: b=np.zeros(N); # bias

for t in xrange(num_trials):
print t, num_trials
# data statistics (clamped)
X=extract_patches(im,SZ,batch_size).astype(np.float);
R_data=np.dot(X,X.T)/batch_size;
mu_data=X.mean(1);

# prior statistics (unclamped)
S=sample(T,b,batch_size,num_init_samples);
R_prior=np.dot(S,S.T)/batch_size;
mu_prior=S.mean(1);

# update params
deltaT=eta*(R_data - R_prior);
T=T+deltaT;

deltab=eta*(mu_data - mu_prior);
b=b+deltab;

return T, b

if __name__ == "__main__":
A = np.array([\
[0.,1.,1.,0],
[1.,1.,0, 0],
[1.,1.,1.,0],
[0, 1.,1.,1.],
[0, 0, 1.,0]
])
T,b = run(A,display=False)
print T
print b

It works by creating patches of data, but this can be modified so the code works on all data all the time.