class RBM(object):"""Restricted Boltzmann Machine (RBM)  """def __init__(self,input=None,n_visible=784,n_hidden=500,W=None,hbias=None,vbias=None,numpy_rng=None,theano_rng=None):"""RBM constructor. Defines the parameters of the model along withbasic operations for inferring hidden from visible (and vice-versa),as well as for performing CD updates.:param input: None for standalone RBMs or symbolic variable if RBM ispart of a larger graph.:param n_visible: number of visible units:param n_hidden: number of hidden units:param W: None for standalone RBMs or symbolic variable pointing to ashared weight matrix in case RBM is part of a DBN network; in a DBN,the weights are shared between RBMs and layers of a MLP:param hbias: None for standalone RBMs or symbolic variable pointingto a shared hidden units bias vector in case RBM is part of adifferent network:param vbias: None for standalone RBMs or a symbolic variablepointing to a shared visible units bias"""self.n_visible = n_visibleself.n_hidden = n_hiddenif numpy_rng is None:# create a number generatornumpy_rng = numpy.random.RandomState(1234)if theano_rng is None:theano_rng = RandomStreams(numpy_rng.randint(2 ** 30))if W is None:# W is initialized with `initial_W` which is uniformely# sampled from -4*sqrt(6./(n_visible+n_hidden)) and# 4*sqrt(6./(n_hidden+n_visible)) the output of uniform if# converted using asarray to dtype theano.config.floatX so# that the code is runable on GPUinitial_W = numpy.asarray(numpy_rng.uniform(low=-4 * numpy.sqrt(6. / (n_hidden + n_visible)),high=4 * numpy.sqrt(6. / (n_hidden + n_visible)),size=(n_visible, n_hidden)),dtype=theano.config.floatX)# theano shared variables for weights and biasesW = theano.shared(value=initial_W, name='W', borrow=True)if hbias is None:# create shared variable for hidden units biashbias = theano.shared(value=numpy.zeros(n_hidden,dtype=theano.config.floatX),name='hbias',borrow=True)if vbias is None:# create shared variable for visible units biasvbias = theano.shared(value=numpy.zeros(n_visible,dtype=theano.config.floatX),name='vbias',borrow=True)# initialize input layer for standalone RBM or layer0 of DBNself.input = inputif not input:self.input = T.matrix('input')self.W = Wself.hbias = hbiasself.vbias = vbiasself.theano_rng = theano_rng# **** WARNING: It is not a good idea to put things in this list# other than shared variables created in this function.self.params = [self.W, self.hbias, self.vbias]

   def propup(self, vis):'''This function propagates the visible units activation upwards tothe hidden unitsNote that we return also the pre-sigmoid activation of thelayer. As it will turn out later, due to how Theano deals withoptimizations, this symbolic variable will be needed to writedown a more stable computational graph (see details in thereconstruction cost function)'''pre_sigmoid_activation = T.dot(vis, self.W) + self.hbiasreturn [pre_sigmoid_activation, T.nnet.sigmoid(pre_sigmoid_activation)]def sample_h_given_v(self, v0_sample):''' This function infers state of hidden units given visible units '''# compute the activation of the hidden units given a sample of# the visiblespre_sigmoid_h1, h1_mean = self.propup(v0_sample)# get a sample of the hiddens given their activation# Note that theano_rng.binomial returns a symbolic sample of dtype# int64 by default. If we want to keep our computations in floatX# for the GPU we need to specify to return the dtype floatXh1_sample = self.theano_rng.binomial(size=h1_mean.shape,n=1, p=h1_mean,dtype=theano.config.floatX)return [pre_sigmoid_h1, h1_mean, h1_sample]def propdown(self, hid):'''This function propagates the hidden units activation downwards tothe visible unitsNote that we return also the pre_sigmoid_activation of thelayer. As it will turn out later, due to how Theano deals withoptimizations, this symbolic variable will be needed to writedown a more stable computational graph (see details in thereconstruction cost function)'''pre_sigmoid_activation = T.dot(hid, self.W.T) + self.vbiasreturn [pre_sigmoid_activation, T.nnet.sigmoid(pre_sigmoid_activation)] def sample_v_given_h(self, h0_sample):''' This function infers state of visible units given hidden units '''# compute the activation of the visible given the hidden samplepre_sigmoid_v1, v1_mean = self.propdown(h0_sample)# get a sample of the visible given their activation# Note that theano_rng.binomial returns a symbolic sample of dtype# int64 by default. If we want to keep our computations in floatX# for the GPU we need to specify to return the dtype floatXv1_sample = self.theano_rng.binomial(size=v1_mean.shape,n=1, p=v1_mean,dtype=theano.config.floatX)return [pre_sigmoid_v1, v1_mean, v1_sample]

  def gibbs_hvh(self, h0_sample):''' This function implements one step of Gibbs sampling,starting from the hidden state'''pre_sigmoid_v1, v1_mean, v1_sample = self.sample_v_given_h(h0_sample)pre_sigmoid_h1, h1_mean, h1_sample = self.sample_h_given_v(v1_sample)return [pre_sigmoid_v1, v1_mean, v1_sample,pre_sigmoid_h1, h1_mean, h1_sample]def gibbs_vhv(self, v0_sample):''' This function implements one step of Gibbs sampling,starting from the visible state'''pre_sigmoid_h1, h1_mean, h1_sample = self.sample_h_given_v(v0_sample)pre_sigmoid_v1, v1_mean, v1_sample = self.sample_v_given_h(h1_sample)return [pre_sigmoid_h1, h1_mean, h1_sample,pre_sigmoid_v1, v1_mean, v1_sample]

 def free_energy(self, v_sample):''' Function to compute the free energy '''wx_b = T.dot(v_sample, self.W) + self.hbiasvbias_term = T.dot(v_sample, self.vbias)hidden_term = T.sum(T.log(1 + T.exp(wx_b)), axis=1)return -hidden_term - vbias_term

 def get_cost_updates(self, lr=0.1, persistent=None, k=1):"""This functions implements one step of CD-k or PCD-k:param lr: learning rate used to train the RBM:param persistent: None for CD. For PCD, shared variablecontaining old state of Gibbs chain. This must be a sharedvariable of size (batch size, number of hidden units).:param k: number of Gibbs steps to do in CD-k/PCD-kReturns a proxy for the cost and the updates dictionary. Thedictionary contains the update rules for weights and biases butalso an update of the shared variable used to store the persistentchain, if one is used."""# compute positive phasepre_sigmoid_ph, ph_mean, ph_sample = self.sample_h_given_v(self.input)# decide how to initialize persistent chain:# for CD, we use the newly generate hidden sample# for PCD, we initialize from the old state of the chainif persistent is None:chain_start = ph_sampleelse:chain_start = persistent

   # perform actual negative phase# in order to implement CD-k/PCD-k we need to scan over the# function that implements one gibbs step k times.# Read Theano tutorial on scan for more information :# http://deeplearning.net/software/theano/library/scan.html# the scan will return the entire Gibbs chain([pre_sigmoid_nvs,nv_means,nv_samples,pre_sigmoid_nhs,nh_means,nh_samples],updates) = theano.scan(self.gibbs_hvh,# the None are place holders, saying that# chain_start is the initial state corresponding to the# 6th outputoutputs_info=[None, None, None, None, None, chain_start],n_steps=k,name="gibbs_hvh")

 # determine gradients on RBM parameters# note that we only need the sample at the end of the chainchain_end = nv_samples[-1]cost = T.mean(self.free_energy(self.input)) - T.mean(self.free_energy(chain_end))# We must not compute the gradient through the gibbs samplinggparams = T.grad(cost, self.params, consider_constant=[chain_end])

  # constructs the update dictionaryfor gparam, param in zip(gparams, self.params):# make sure that the learning rate is of the right dtypeupdates[param] = param - gparam * T.cast(lr,dtype=theano.config.floatX)if persistent:# Note that this works only if persistent is a shared variableupdates[persistent] = nh_samples[-1]# pseudo-likelihood is a better proxy for PCDmonitoring_cost = self.get_pseudo_likelihood_cost(updates)else:# reconstruction cross-entropy is a better proxy for CDmonitoring_cost = self.get_reconstruction_cost(updates,pre_sigmoid_nvs[-1])return monitoring_cost, updates

追踪进程

由于分割函数Z，RBMs训练起来很困难。我们无法在训练时估计指数相似log(P(x))，因此没有办法直接衡量以选择最优的超参数。

这里仅提供一些想法。

检查负样本

训练过程中得到的负样本可以图形化。随着训练进行，我们知道模型定义的RBM与实际的Ptrain(x)分布越来越接近。负样本因此应该与训练集中的样本相似。明显的不良超参数可以以这种方式舍弃。

图形化检视过滤器

模型学习到的过滤器可以图形化。相当于勾画每个单元权重的灰图（在重新调整为平方矩阵）。过滤器应该找出数据中的突出特征。虽然对于任意数据集很难说特征应该怎样，对于MNIST的训练过滤器一般能捕捉到笔画，而对于自然图像的训练如果结合稀疏标准则过滤器类似于Gabor。

相似性的代理

其他更能追踪的函数可以作为相似性的代理。当使用PCD训练RBM时，可以使用伪相似作为代理。伪相似（PL）计算代价小，因为它假设所有部分互相独立

这里x-i表示除了i的所有x。log-PL因此是每个Xi指数概率之和，给定其他所有部分状态。对于MNIST，涉及加总所有784输入维度，依然昂贵。因此我们使用log-PL的随机相似

期望由索引i的一致随机选择，N是可见单元的数量。为进行二进制操作，我们进一步引进指代x，i比特互反(1->0, 0->1)。RBM的log-PL写成

我们因此在RBM类的get_cost_updates函数返回这个成本以及RBM更新。 注意我们调整更新字典来递增i的索引。这会导致i遍历所有可能的值{0,1,...,N}。

注意对于输入和重构的（同降噪自动编码机）CD训练交叉熵成本比伪指数相似更可靠。以下是我们计算伪相似的代码：

def get_pseudo_likelihood_cost(self, updates):"""Stochastic approximation to the pseudo-likelihood"""# index of bit i in expression p(x_i | x_{\i})bit_i_idx = theano.shared(value=0, name='bit_i_idx')# binarize the input image by rounding to nearest integerxi = T.round(self.input)# calculate free energy for the given bit configurationfe_xi = self.free_energy(xi)# flip bit x_i of matrix xi and preserve all other bits x_{\i}# Equivalent to xi[:,bit_i_idx] = 1-xi[:, bit_i_idx], but assigns# the result to xi_flip, instead of working in place on xi.xi_flip = T.set_subtensor(xi[:, bit_i_idx], 1 - xi[:, bit_i_idx])# calculate free energy with bit flippedfe_xi_flip = self.free_energy(xi_flip)# equivalent to e^(-FE(x_i)) / (e^(-FE(x_i)) + e^(-FE(x_{\i})))cost = T.mean(self.n_visible * T.log(T.nnet.sigmoid(fe_xi_flip -fe_xi)))# increment bit_i_idx % number as part of updatesupdates[bit_i_idx] = (bit_i_idx + 1) % self.n_visiblereturn cost

  # it is ok for a theano function to have no output# the purpose of train_rbm is solely to update the RBM parameterstrain_rbm = theano.function([index],cost,updates=updates,givens={x: train_set_x[index * batch_size: (index + 1) * batch_size]},name='train_rbm')plotting_time = 0.start_time = timeit.default_timer()# go through training epochsfor epoch in range(training_epochs):# go through the training setmean_cost = []for batch_index in range(n_train_batches):mean_cost += [train_rbm(batch_index)]print('Training epoch %d, cost is ' % epoch, numpy.mean(mean_cost))# Plot filters after each training epochplotting_start = timeit.default_timer()# Construct image from the weight matriximage = Image.fromarray(tile_raster_images(X=rbm.W.get_value(borrow=True).T,img_shape=(28, 28),tile_shape=(10, 10),tile_spacing=(1, 1)))image.save('filters_at_epoch_%i.png' % epoch)plotting_stop = timeit.default_timer()plotting_time += (plotting_stop - plotting_start)end_time = timeit.default_timer()pretraining_time = (end_time - start_time) - plotting_timeprint ('Training took %f minutes' % (pretraining_time / 60.))

 ##################################     Sampling from the RBM     ################################### find out the number of test samplesnumber_of_test_samples = test_set_x.get_value(borrow=True).shape[0]# pick random test examples, with which to initialize the persistent chaintest_idx = rng.randint(number_of_test_samples - n_chains)persistent_vis_chain = theano.shared(numpy.asarray(test_set_x.get_value(borrow=True)[test_idx:test_idx + n_chains],dtype=theano.config.floatX))

 plot_every = 1000# define one step of Gibbs sampling (mf = mean-field) define a# function that does `plot_every` steps before returning the# sample for plotting([presig_hids,hid_mfs,hid_samples,presig_vis,vis_mfs,vis_samples],updates) = theano.scan(rbm.gibbs_vhv,outputs_info=[None, None, None, None, None, persistent_vis_chain],n_steps=plot_every,name="gibbs_vhv")# add to updates the shared variable that takes care of our persistent# chain :.updates.update({persistent_vis_chain: vis_samples[-1]})# construct the function that implements our persistent chain.# we generate the "mean field" activations for plotting and the actual# samples for reinitializing the state of our persistent chainsample_fn = theano.function([],[vis_mfs[-1],vis_samples[-1]],updates=updates,name='sample_fn')# create a space to store the image for plotting ( we need to leave# room for the tile_spacing as well)image_data = numpy.zeros((29 * n_samples + 1, 29 * n_chains - 1),dtype='uint8')for idx in range(n_samples):# generate `plot_every` intermediate samples that we discard,# because successive samples in the chain are too correlatedvis_mf, vis_sample = sample_fn()print(' ... plotting sample %d' % idx)image_data[29 * idx:29 * idx + 28, :] = tile_raster_images(X=vis_mf,img_shape=(28, 28),tile_shape=(1, n_chains),tile_spacing=(1, 1))# construct imageimage = Image.fromarray(image_data)image.save('samples.png')

... loading data
Training epoch 0, cost is  -90.6507246003
Training epoch 1, cost is  -81.235857373
Training epoch 2, cost is  -74.9120966945
Training epoch 3, cost is  -73.0213216101
Training epoch 4, cost is  -68.4098570497
Training epoch 5, cost is  -63.2693021647
Training epoch 6, cost is  -65.99578971
Training epoch 7, cost is  -68.1236650015
Training epoch 8, cost is  -68.3207365087
Training epoch 9, cost is  -64.2949797113
Training epoch 10, cost is  -61.5194867893
Training epoch 11, cost is  -61.6539369402
Training epoch 12, cost is  -63.5465278086
Training epoch 13, cost is  -63.3787093527
Training epoch 14, cost is  -62.755739271
Training took 122.466000 minutes... plotting sample  0... plotting sample  1... plotting sample  2... plotting sample  3... plotting sample  4... plotting sample  5... plotting sample  6... plotting sample  7... plotting sample  8... plotting sample  9