判别模型的玻尔兹曼机论文源码解读
前言
三號要去參加CAD/CG會議，投了一篇關于使用生成模型和判別模型的RBM做運動捕捉數據風格識別的論文。這段時間一直搞卷積RBM了，差點把原來的實驗內容都忘記了，這里復習一下判別式玻爾茲曼機的訓練流程。
國際慣例，貼幾個鏈接：
論文1——Energy Based Learning Classification task using Restricted Boltzmann Machines
鏈接：http://pan.baidu.com/s/1i5foeEx 密碼：flq7
論文2——Classification using Discriminative Restricted Boltzmann Machines
鏈接：http://pan.baidu.com/s/1qYGT9z2 密碼：hebj
代碼——RBM-on-Classification
github：https://github.com/atzenigor/RBM-on-Classification
鏈接：http://pan.baidu.com/s/1gf41MKz 密碼：z3jq
數據——mnist手寫數字數據庫
鏈接：http://pan.baidu.com/s/1boGxPxX 密碼：0l25
本文主要是結合論文1及其代碼對判別式限制玻爾茲曼機進行解讀，代碼主要使用在手寫數字識別mnist數據庫中。
先畫一下代碼對應的模型結構圖
第一步
從主函數入手：trainingClassRBM.m
首先是一系列的初始化，就不說了，主要包含學習率、動量項、權重衰減、分批大小、訓練次數，隨機初始化權重和偏置等。這里主要關注兩行程序：
discrim = 1; % 1 to activate the discriminative learning, 0 otherwisefreeenerg = 1; % if discrim is zero but one want to compute the free energy對于這兩句話，可以翻開論文1的第五章(第11頁)說過這樣一句話：
大概意思就是，為了解決特定的分類問題，我們可以簡單地定義一個目標函數去進行最小化，進而得到RBM的模型參數。而目標函數可以選擇三種：
1、生成模型的訓練目標函數
2、判別模型的訓練目標函數
3、混合模型的訓練目標函數
顯然，根據代碼，可以發現程序使用了變量discrim來控制是否使用第三種方法，如果discrim=1則使用第三種目標函數，如果discrim=0則使用第一種目標函數。
此后也進行了數據的分批處理和驗證集的選取
data_valid = data(end - valid_size + 1 : end, : ); label_valid = labels(end - valid_size + 1 : end, : );%convert y in the notetion one out of number of classes lab = boolean(set_y(labels(1:end - valid_size,:), num_class)); %convert the dataset in batches [batch_data, batch_labels , batch_size] = ...createBatches(data(1:end - valid_size,:), lab, num_batch);第二步
正式進入第一層for循環，用于控制訓練次數，每次訓練完畢都會輸出相應的幾個信息，包含重構誤差，驗證集誤差，自由能變化。
sum_rec_error = 0;sum_class_error = 0;sum_tot_free_energ = 0;下一句就是在進行第五次迭代時候，將動量項改變 if i_epoch == 5 momentum = final_momentum;end;并且使用模擬退火（simulated annealing）算法降低學習率 %simulated annealinggamma = (init_gamma) / (1 + (i_epoch - 1) / tau);第三步
正式進入內層循環對數據的分配處理，每次只取一批數據，大小為100（最后一批大小根據情況而定）
x = batch_data{iBatch}';y = batch_labels{iBatch}';batch_size=size(x,2);進入生成模型的計算，即計算第一步提到的第一個目標函數對應的模型參數更新。首先計算隱藏層的激活概率值P(h|x,y)，這就是在RBM中常說的positive phase根據文章可以得到
而相應的程序就是 hprob = prob_h_given_xy(x, y, w, u, b_h); 其中函數prob_h_given_xy為： function [ prob ] = prob_h_given_xy( x, y, w, u, b_h ) <span style="white-space:pre"> </span>prob = sigmf(bsxfun(@plus, w * x + u * y, b_h), [1 0]); end隨后進行隱單元的隨機激活，這是RBM中經常出現的一步驟，我們稱為激活狀態： h = sampling(hprob);激活函數實現如下： function [ samp ] = sampling( p ) <span style="white-space:pre"> </span>samp = p > rand(size(p)); end然后便是對應的negative phase，使用得到的隱單元狀態反向計算x和y單元數據根據論文介紹的激活函數
可以找到對應的激活函數程序寫法：
計算P(x|h)，很簡單，直接計算權重乘以隱單元的激活狀態(并非激活概率值)然后加上輸入層偏置即可。
function [ prob ] = prob_x_given_h( h, w, b_x ) <span style="white-space:pre"> </span>prob = sig(bsxfun(@plus, w' * h, b_x)); end計算P(y|h)，稍微復雜點，在分母的處理上，需要進行加和，u'*h計算得到的矩陣大小是類別數*批大小，實驗中大小是30*716，而sum是將所有行加到第一行，也就是說對每一個樣本計算得到的標簽層單元數據進行加和，整個過程類似于進行了歸一化。 function [ prob ] = prob_y_given_h( h, u , b_y ) <span style="white-space:pre"> </span>prob_not_norm = exp(bsxfun(@plus, u' * h, b_y)); <span style="white-space:pre"> </span>norm_factor = sum(prob_not_norm); <span style="white-space:pre"> </span>prob = bsxfun(@rdivide, prob_not_norm, norm_factor); end 計算得到x和y的激活概率值之后同樣需要激活狀態 xrprob = prob_x_given_h( h, w, b_x ); yrprob = prob_y_given_h( h, u, b_y ); xr = sampling(xrprob); yr = sampling_y(yrprob);進行negative phase的最后一步，利用新的x和y的狀態值(0或者1)去計算隱層的激活概率值 hrprob = prob_h_given_xy( xr, yr, w, u, b_h);根據以上計算的positive phase和negative phase相關值計算生成模型各參數的更新梯度，主要包含隱層至輸入層的連接權重w，隱層至標簽層的連接權重v，輸入層偏置b_x，隱藏層偏置b_h，標簽層偏置b_y。注意更新的時候用的是激活概率值去更新，而非激活狀態值去更新。 g_w_gen = (hprob * x' - hrprob * xrprob') / batch_size; g_u_gen = (hprob * y' - hrprob * yrprob') / batch_size; g_b_x_gen = (sum(x,2) - sum(xrprob,2)) / batch_size; g_b_y_gen = (sum(y,2) - sum(yrprob,2)) / batch_size; g_b_h_gen = (sum(hprob,2) - sum(hrprob,2)) / batch_size;第四步
進入判別模型的目標函數計算，上面第一節介紹的第二個目標函數。
當然首先就是判斷我們是否選擇使用判別模型，以及初始化判別模型計算中的相關參數，包含隱層偏置梯度，兩個權重w和u的梯度
if (discrim || freeenerg)%initialization of the gradientsg_b_h_disc_acc = 0;g_w_disc_acc = 0;g_u_disc_acc = zeros(h_size,num_class);特別注意，此處是沒有輸入層偏置的，詳細可看論文1中的這樣一句話：
意思就是b（文章的b是輸入層的偏置）的梯度為0，因為計算P(y|x)的時候并沒有用到輸入層的偏置。
讓我們看看P(y|x)的計算方法
其中F 就代表自由能量(free energy)。
還是寫一下比較好，后面的代碼會單獨計算這個自由能量，其中dy是標簽層偏置
后面的步驟便是逐步實現這個采樣函數
注意這里的softplus也是一個函數，而非簡單的加法，具體介紹點擊此處
①首先發現式子中始終有w*x這一項，所以提出來計算一下，避免后面的重復計算
%ausiliary variablesw_times_x = w * x;②接下來計算自由能量函數 for iClass= 1:num_classo_t_y(:,:,iClass) = bsxfun(@plus, b_h + u(:,iClass), w_times_x );neg_free_energ(iClass,:) = b_y(iClass) + sum(log(1 + exp(o_t_y(:,:,iClass)))); end;注意這里的兩個變量的維度：
o_t_y維度大小是800*716*30，分別代表隱單元神經元個數、一批樣本數、標簽總數，就像對隱層求了30個800*716的單元數值，然后后面計算的時候也確實如此，每次計算都用了一個for循環對第三個維度的800*716矩陣分別遍歷，應該就是為了計算此樣本屬于每一個標簽的概率。
neg_free_energ大小是30*716，分別代表標簽總數、樣本總數
計算的時候是依次計算標簽層每個神經元所代表的所有樣本之和。有點像先用x和y計算得到h，然后將h經過softplus映射，最后將每一個樣本的隱層800個單元求和，得到一個樣本在對應的標簽單元的值。比如o_t_y計算iClass=1時候得到的其實是一個800*716的矩陣，然后經過sum得到1*716的向量，直接作為標簽層的第一個單元的所有樣本的值即可。
然后記錄了一下當前批數據的總共的能量和，用于與下一批的自由能量和相加
sum_tot_free_energ = sum_tot_free_energ + sum(-neg_free_energ(y));至此“ if (discrim || freeenerg)”中的freeenerg計算完畢。③接下來就是計算discrim對應的參數更新梯度了
主要就是計算P(y|x)根據上面提到的式子。
程序中首先對求得的neg_free_energ進行了零均值化，并計算出了分子的值
med = mean2(neg_free_energ); p_y_given_x = exp(neg_free_energ - med);利用分子的值，計算分子/分母的值，注意這里的加是將每個樣本對應的30個標簽層值相加，即sum(p_y_given_x)的結果是將30*716加和成了1*716維度的向量 p_y_given_x = bsxfun(@rdivide, p_y_given_x, sum(p_y_given_x));④最后就是計算梯度了首先計算標簽層的梯度，按照原文的公式
說白了就是用原始標簽減去P(y|x)推理得到的標簽值，最后計算一下加和即可，便得到了標簽層偏置更新梯度
dif_y_p_y = y - p_y_given_x;g_b_y_disc_acc = sum(dif_y_p_y,2);隨后便是計算剩下的三個量權重w、v和隱層偏置b_h的梯度，查看論文1的計算方法
其中
相關的計算程序如下：
首先計算第一項和第二項的 O 函數
sig_o_t_y = sig(o_t_y(:,:,iClass)); sig_o_t_y_selected = sig_o_t_y(:,y(iClass,:)); 然后分別套入對b_h、w、u的計算中 g_b_h_disc_acc = g_b_h_disc_acc + sum(sig_o_t_y_selected, 2) - sum(bsxfun(@times, sig_o_t_y, p_y_given_x(iClass,:)),2); %update the gradient of w fot the class 'iClass' g_w_disc_acc = g_w_disc_acc + sig_o_t_y_selected * x(:,y(iClass,:))' - bsxfun(@times, sig_o_t_y, p_y_given_x(iClass,:)) * x'; %update the gradient of u fot the class 'iClass' g_u_disc_acc(:,iClass) = sum(bsxfun(@times, sig_o_t_y, dif_y_p_y(iClass,:)), 2);注意這兩部分代碼要用一個for循環包含起來 for iClass= 1:num_class主要也就是上面說過的，對o_t_y(維度為隱層單元數*一批樣本數*標簽類別總數)的第三維進行遍歷，每次遍歷得到的隱單元數*批樣本數作為隱層單元即sig_o_t_y(代碼維度是800*716)，然后從中抽取出標簽為iClass的樣本(比如716個樣本中，有29個樣本標簽為當前循環的iClass)得到sig_o_t_y_selected(iClass=1時文章維度為800*29)這樣就比較好理解了，就像生成模型的RBM一樣：計算隱層偏置b_h的時候，計算sig_o_t_y_selected與sig_o_t_y*P(y|x)的差值即可；
計算連接權重w的時候，計算sig_o_t_y_selected*x_selected(這個變量代表批樣本x中標簽為iClass的樣本矩陣)與sig_o_t_y*P(y|x)*x對應的差值即可；
計算連接權重u的時候，只需計算sig_o_t_y與dif_y_p_y(iClass,:)的乘積即可，因為在計算dif_y_p_y的時候，已經計算過差值了，所以沒有像上面計算b_h和w一樣去做差。
⑤全部計算完了，就更新一下判別模型中的幾個參數吧
g_b_y_disc = g_b_y_disc_acc / batch_size; g_b_h_disc = g_b_h_disc_acc / batch_size; g_w_disc = g_w_disc_acc / batch_size; g_u_disc = g_u_disc_acc / batch_size;第五步
利用第一步介紹的目標函數三進行全局梯度更新
delta_w = delta_w * momentum + ...gamma * (discrim * g_w_disc + alpha * g_w_gen - weight_cost * w); delta_u = delta_u * momentum + ...gamma * (discrim * g_u_disc + alpha * g_u_gen - weight_cost * u); delta_b_x = delta_b_x * momentum + ...gamma * alpha * g_b_x_gen; delta_b_y = delta_b_y * momentum + ...gamma * (discrim * g_b_y_disc + alpha * g_b_y_gen); delta_b_h = delta_b_h * momentum + ...gamma * (discrim * g_b_h_disc + alpha * g_b_h_gen); w = w + delta_w; u = u + delta_u; b_x = b_x + delta_b_x; b_y = b_y + delta_b_y; b_h = b_h + delta_b_h;后面做了一個平均，表示沒太懂是干什么用的，好像是削減梯度？而且在每次迭代輸出的驗證集誤差的時候使用的依舊是上面的b_y, b_h, w, u 而非做平均的b_y_avg, b_h_avg, w_avg, u_avg 。這個還有待探討。第六步
獲得所有參數以后需要進行識別，很簡單，有一個輸入testdata，對應的標簽testlabels，總類別數目num_class，對應的模型參數b_y, b_h, w, u，然后利用能量函數F的計算得到樣本屬于每個標簽的概率，前面也說了，neg_free_energ的大小是標簽數*樣本數，這樣就能得到每一個樣本對應的每一個標簽的概率是對少了，取最大的那個概率就是了。
function [ err ] = predict( testdata, testlabels, num_class, b_y, b_h, w, u )numcases= size(testdata, 1);w_times_x = w * testdata';neg_free_energ = zeros(num_class, numcases);for iClass= 1:num_classneg_free_energ(iClass,:) = b_y(iClass) + sum(log(1 + ...exp(bsxfun(@plus, b_h + u(:,iClass), w_times_x ))));end;[~, prediction] = max(neg_free_energ);err = sum(prediction ~= testlabels') / numcases * 100;end第七步
補充一下，如何看模型的效果好壞呢？
首先重構誤差不必多說，肯定是慢慢降低
其次就是自由能量了，根據RBM的描述，系統最穩定的時刻就是能量最少的時候。
當然咯，分類誤差肯定也是降低啦
還有，有時候能量會突然升高，這時候不要慌，先等等，看看是否一直升高，有可能在某次迭代完成就下降了呢，類似于這樣
總結
以上是生活随笔為你收集整理的判别模型的玻尔兹曼机论文源码解读的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 【caffe-Windows】微软官方c
- 下一篇： 游戏中实现鼠标拖尾效果