前言

卷積RBM相對RBM來說具有很多優(yōu)勢，詳細的我也不說了，看文章就行。主要還是為了加深自己對細節(jié)部分的理解吧。

國際慣例，貼幾個鏈接

卷積RBM的創(chuàng)始人Honglak Lee：http://web.eecs.umich.edu/~honglak/hl_publications.html#embc15

可以從他的主頁上找到兩篇與卷積RBM相關的文獻，我也順便傳到CSDN上了，賺兩個積分，嘿嘿：

第一篇：在2010年的那一欄的《Convolutional deep belief networks for scalable unsupervised learning of hierarchical representations》

文章及代碼下載：http://download.csdn.net/detail/zb1165048017/9683175

第二篇：在2010年的那一欄的《Unsupervised feature learning for audio classification using convolutional deep belief networks》

文章及代碼下載：http://download.csdn.net/detail/zb1165048017/9683181

本片博客所解讀的源碼地址：https://github.com/lonl/CDBN

【注】此代碼目前不支持步長大于1的卷積，因為matlab自帶的conv卷積函數(shù)未提供對不同步長的卷積操作，需要自己修改。

卷積RBM的模型與CNN的模型結構很類似：

第一步

分析主函數(shù)，主要是讀取數(shù)據(jù)集，總共9張圖片，這里訓練的時候采用第一張圖片，也可以同時讀入兩張圖片進行重構。無需隨機打亂圖片，看到后面會發(fā)現(xiàn)有隨機打亂zheyibuzhou

for i = 1:9str = ['./data/MITcoast/image_000',num2str(i)];str1 = strcat(str,'.jpg');image = imread(str1);image = double(image)/255;train_data(:,:,i) = image; endtrain_data = reshape(train_data,[256,256,1,9]); train_data = train_data(:,:,:,1:1);

這里注意一下train_data的每個維度分別代表什么，前面兩維表示輸入圖片的大小，第三個維度是通道數(shù)(比如彩色圖片的通道就是3)，第四個維度是輸入的圖片的個數(shù)(如果是輸入兩張，那么就可以寫成"1:2")。然后初始化兩層CRBM對應的參數(shù)(這里面重點注意n_map_v是指的通道數(shù)，而非圖片個數(shù))：

layer{1} = default_layer2D(); % DEFAULT PARAMETERS SETTING, % YOU CAN CHANGE THE PARAMETERS IN THE FOLLOWING LINESlayer{1}.inputdata = train_data;%輸入數(shù)據(jù) layer{1}.n_map_v = 1;%每張圖片的通道數(shù)，可以根據(jù)train_data的第三個維度調(diào)整 layer{1}.n_map_h = 9;%卷積核個數(shù) layer{1}.s_filter = [7 7];%卷積核大小 layer{1}.stride = [1 1]; %卷積步長，行方向和列方向 layer{1}.s_pool = [2 2];%池化 layer{1}.n_epoch = 10;%迭代次數(shù) layer{1}.learning_rate = 0.05;%學習率 layer{1}.sparsity = 0.03;%稀疏化 layer{1}.lambda1 = 5; layer{1}.lambda2 = 0.05; layer{1}.whiten = 1;%是否進行白化處理 layer{1}.type_input = 'Gaussian'; % OR 'Gaussian' 'Binary'輸入的是高斯還是二值數(shù)據(jù)% SECOND LAYER SETTING layer{2} = default_layer2D(); % DEFAULT PARAMETERS SETTING, % YOU CAN CHANGE THE PARAMETERS IN THE FOLLOWING LINESlayer{2}.n_map_v = 9; layer{2}.n_map_h = 16; layer{2}.s_filter = [5 5]; layer{2}.stride = [1 1]; layer{2}.s_pool = [2 2]; layer{2}.n_epoch = 10; layer{2}.learning_rate = 0.05; layer{2}.sparsity = 0.02; layer{2}.lambda1 = 5; layer{2}.lambda2 = 0.05; layer{2}.whiten = 1; layer{2}.type_input = 'Gaussian';關于白化處理的作用，請參考這篇博客： http://blog.csdn.net/whiteinblue/article/details/36171233
接下來便是根據(jù)設置的參數(shù)去訓練卷積DBN

tic;[model,layer] = cdbn2D(layer); save('./model/model_parameter','model','layer');toc;后面就是畫兩層C-DBN分別得到的重構圖

% THE FIRST LAYER POOLING MAPS figure(1); [r,c,n] = size(model{1}.output(:,:,:,1)); visWeights(reshape(model{1}.output(:,:,:,1),r*c,n)); colormap gray title(sprintf('The first Pooling output')) drawnow% THE SECOND LAYER POOLING MAPS figure(2); [r,c,n] = size(model{2}.output(:,:,:,1)); visWeights(reshape(model{2}.output(:,:,:,1),r*c,n)); colormap gray title(sprintf('The second Pooling output')) drawnow% ORIGINAL SAMPLE figure(3); imagesc(layer{1}.inputdata(:,:,:,1)); colormap gray; axis image; axis off title(sprintf('Original Sample'));重點就是關注訓練函數(shù)cdbn2D.m

第二步

主函數(shù)中調(diào)用的的cdbn2D函數(shù)解讀

首先就是去訓練第一層

layer{1} = preprocess_train_data2D(layer{1});fprintf('layer 1:\n');model{1} = crbm2D(layer{1});str1 = sprintf('./model/model_%s_%s_%d',layer{1}.type_input,layer{1}.cpu,1);save(str1,'model','layer');然后訓練第二層及以后的層數(shù)

for i = 2:length(layer)fprintf('layer %d:\n',i);layer{i}.inputdata = model{i-1}.output;layer{i} = preprocess_train_data2D(layer{i});model{i} = crbm2D(layer{i});str1 = sprintf('./model/model_%s_%s_%d',layer{i}.type_input,layer{i}.cpu,i);save(str1,'model','layer');end可以發(fā)現(xiàn)，最重要的函數(shù)包含兩部分：

第一部分是白化函數(shù)preprocess_train_data2D，第二部分是訓練函數(shù)crbm2D

第三步

白化函數(shù)的解讀：首先對輸入數(shù)據(jù)進行裁剪，使得運動在經(jīng)過卷積池化以后認為整數(shù)的大小，所以前面幾行代碼就是對圖片長寬裁剪

mod_1 = mod((size(layer.inputdata,1)-layer.s_filter(1))/layer.stride(1)+1,layer.s_pool(1)); if mod_1~=0layer.inputdata(1:floor(mod_1/2),:,:,:) =[];layer.inputdata(end-ceil(mod_1/2)+1:end,:,:,:) =[]; endmod_2 = mod((size(layer.inputdata,2)-layer.s_filter(2))/layer.stride(2)+1,layer.s_pool(2)); if mod_2~=0layer.inputdata(:,1:floor(mod_2/2),:,:) =[];layer.inputdata(:,end-ceil(mod_2/2)+1:end,:,:) =[]; end然后進入白化處理階段，兩個for循環(huán)，外層的控制選取圖片，內(nèi)層的控制通道。也就是說對每一張圖片的每一個通道單獨進行白化處理。

if layer.whiten if strcmp(layer.type_input, 'Gaussian')m = size(layer.inputdata,4);n = size(layer.inputdata,3);for i = 1 : mfor j =1 : nlayer.inputdata(:,:,j,i) = crbm_whiten(layer.inputdata(:,:,j,i));endendend end具體的白化方法如下：

function im_out = crbm_whiten(im)if size(im,3)>1im = rgb2gray(im); endim = double(im); im = im - mean(im(:)); im = im./std(im(:));N1 = size(im, 1); N2 = size(im, 2);[fx fy]=meshgrid(-N1/2:N1/2-1, -N2/2:N2/2-1); rho=sqrt(fx.*fx+fy.*fy)';f_0=0.4*mean([N1,N2]); filt=rho.*exp(-(rho/f_0).^4);If=fft2(im); imw=real(ifft2(If.*fftshift(filt)));im_out = 0.1*imw/std(imw(:)); % 0.1 is the same factor as in make-your-own-imagesend首先歸一化，然后傅里葉變換、反傅里葉變換等。具體過程以及函數(shù)可以去matlab論壇中找到相關使用方法。我們主要關注下一個函數(shù)crbm2D

第四步

核心部分crbm2D函數(shù)的解讀。跟RBM一樣，首先一堆初始化

cpu = layer.cpu; batchsize = layer.batchsize; model.n_cd = layer.n_cd; model.momentum = layer.momentum; model.start_gau = layer.start_gau; model.stop_gau = layer.stop_gau; model.beginAnneal = Inf; layer.s_inputdata = [size(layer.inputdata,1),size(layer.inputdata,2)];% INITIALIZE THE WEIGHTS model.W = 0.01*randn(layer.s_filter(1), layer.s_filter(2), layer.n_map_v, layer.n_map_h);%卷積核大小*卷積核個數(shù)*輸入圖個數(shù) model.dW = zeros(size(model.W)); model.v_bias = zeros(layer.n_map_v, 1);%輸入層偏置，就一個1*1的值 model.dV_bias = zeros(layer.n_map_v, 1); model.h_bias = zeros(layer.n_map_h, 1);%隱藏層偏置是對每一個特征圖都是1*1的值，所以總共是卷積核個數(shù)*1個偏置 model.dH_bias = zeros(layer.n_map_h, 1); model.v_size = [layer.s_inputdata(1), layer.s_inputdata(2)]; model.v_input = zeros(layer.s_inputdata(1), layer.s_inputdata(2), layer.n_map_v,batchsize);%輸入層大小是圖片大小*輸入圖片個數(shù)*批大小 model.h_size = (layer.s_inputdata - layer.s_filter)./(layer.stride) + 1; model.h_input = zeros(model.h_size(1), model.h_size(2), layer.n_map_h, batchsize);%隱層特征圖大小*卷積核個數(shù)*每批數(shù)據(jù)個數(shù) model.error = zeros(layer.n_epoch,1);% ADD SOME OTHER PARAMETERS FOR TEST model.winc = 0; model.hinc = 0; model.vinc = 0;% NEED TO FIX THE STRIDE HERE H_out = ((size(layer.inputdata,1)-layer.s_filter(1))/layer.stride(1)+1)/layer.s_pool(1); W_out = ((size(layer.inputdata,2)-layer.s_filter(2))/layer.stride(2)+1)/layer.s_pool(2); model.output = zeros([H_out, W_out,layer.n_map_h,size(layer.inputdata,4)]);%模型的輸出就是特征圖大小*卷積核個數(shù)*(所有批合起來的)總圖片張數(shù)% CREATING BATCH DATA N = size(layer.inputdata,4);%總圖片張數(shù) numcases = size(layer.inputdata,4); numbatches = ceil(N/batchsize);%分批以后的批數(shù)目 groups = repmat(1:numbatches, 1, batchsize); groups = groups(1:N); perm = randperm(N); groups = groups(perm);%每次訓練都會隨機打亂輸入圖片 dW = zeros(size(model.dW));然后就可以獲得每次訓練所取的被隨機打亂的圖片了，上面的creating batch data的部分是對每一個數(shù)據(jù)都分組編號，假設有90張圖片，分為9批，每批10張圖片，那么每一張圖片的編號就是0~9，然后利用下面這個for循環(huán)分別取出隸屬于每一個編號的所有數(shù)據(jù)。我想說，好機智哇，我每次都是直接用一個隨機數(shù)，來打亂，都沒考慮這么多。o(╯□╰)o

for i = 1:numbatchesbatchdata{i} = layer.inputdata(:,:,:,groups == i); end

再對數(shù)據(jù)以及網(wǎng)絡參數(shù)都進行預處理以后，就正式進入訓練過程，正常訓練方法：兩個for循環(huán)，外層循環(huán)控制訓練次數(shù)，內(nèi)層循環(huán)控制每次訓練所使用的分批數(shù)據(jù)集。而且這里我們只對matlab版本進行解讀，有興趣的也可以看看mex和cuda版本

for epoch = 1:layer.n_epocherr = 0;original_err=0;sparsity = zeros(1,numbatches);tic;% FOR EACH EPOCH, ALL SAMPLES ARE COMPUTEDfor i = 1:numbatchesbatch_data = batchdata{i};%-----------------------------------------------------------------%switch cpucase 'mex'%----------------- HERE COMPLETE WITH MEX FUNCTION -----------%[model_new] = crbm2D_batch_mex(model,layer,batch_data);dW = model_new.dW;model.v_sample = model_new.v_sample;model.h_sample = model_new.h_sample;model.h_sample_init = model_new.h_sample_init;case 'cuda'%-------------- HERE COMPLETE WITH CUDA FUNCTION ---------%[model_new] = crbm2D_mex_cuda(model,layer,batch_data);dW = model_new.dW;model.v_sample = model_new.v_sample;model.h_sample = model_new.h_sample;model.h_sample_init = model_new.h_sample_init; case 'matlab'%-------------- HERE COMPLETE WITH MATLAB FUNCTION -------%[model, dW] = calc_gradient2D(model,layer,batch_data);end % switch好了，進入核心函數(shù)的核心部分calc_gradient2D，就是利用對比散度去計算參數(shù)更新的梯度了

第五步

前向推斷部分，在RBM中經(jīng)常稱為positive phase，直白點就是利用可見層計算隱層的激活狀態(tài)，函數(shù)名：crbm_inference2D

先回看一下model.h_input這個變量代表的是什么：

model.h_input = zeros(model.h_size(1), model.h_size(2), layer.n_map_h, batchsize);%隱層特征圖大小*卷積核個數(shù)*每批數(shù)據(jù)個數(shù)意思就是對于每批數(shù)據(jù)中的每一張圖片都有 (卷積核個數(shù)) 張?zhí)卣鲌D。
首先是將所有的特征圖初始化為0，并且記錄下當前處理的批的數(shù)據(jù)量，也就是圖片張數(shù)

model.h_input = zeros(size(model.h_input)); N = size(model.h_input,4);接下來進入三個嵌套的for循環(huán)，外層代表當前需要處理的批的數(shù)據(jù)個數(shù)，中層代表卷積核個數(shù)，內(nèi)層代表每張圖片的通道數(shù)，文章采用的是灰度圖，所以內(nèi)層循環(huán)可有可無，但是如果對于彩色圖的話，那就是三個通道了，彩色圖像的卷積我前面有篇博客提到過一個比較好的博客，o(╯□╰)o突然我也找不到了，回頭找到了再貼過來。反正就是對于三個通道分別進行卷積核的卷積，最后加和起來就行了。看看代碼也是這樣：

for k = 1:Nfor i = 1 : layer.n_map_hfor j = 1 : layer.n_map_vmodel.h_input(:,:,i,k) = model.h_input(:,:,i,k) + conv2(data(:,:,j,k),model.W(end:-1:1,end:-1:1,j,i),'valid');endmodel.h_input(:,:,i,k) = model.h_input(:,:,i,k) + model.h_bias(i);end end這里注意一下，matlab的卷積會自動對卷積核進行180°的翻轉，有興趣可以去知乎上看看卷積和相干這兩個比較相似的概念。是不是因為這個原因所以在這里進行了類似的翻轉，其實matlab里面自帶了翻轉90°的方法，叫rot90(X,90)。驗證一下看看是不是一樣的。
>> A=[1 2 3 4;5 6 7 8]A =1 2 3 45 6 7 8>> A(end:-1:1,end:-1:1)ans =8 7 6 54 3 2 1>> rot90(A,2)ans =8 7 6 54 3 2 1在得到當前批的所有數(shù)據(jù)的卷積特征圖以后，進入池化相關的函數(shù)crbm_blocksum2D

首先判斷隱層是要求二值形式還是高斯形式，有不同的處理

if strcmp(layer.type_input, 'Binary')h_input = exp(model.h_input); elseh_input = exp(1.0/(model.start_gau^2).*model.h_input); end接下來就是對每一張?zhí)卣鲌D進行池化了，可以發(fā)現(xiàn)有三個for循環(huán)，外層是取出每個圖片對應的(卷積核個數(shù)張)特征圖，中層循環(huán)是提取出每次被池化的部分的行，比如特征圖大小是row*col，用[2 2]的卷積核，中層循環(huán)就是吧row*col切成row/2個2*col的特征圖，這樣內(nèi)層循環(huán)不言而喻，就是將2*col的特征圖切成col/2個2*2大小的特征圖，這樣就能提出來每次需要被池化的2*2塊

for k = 1:Nfor i = 1:floor(row/y_stride)offset_r = ((i-1)*y_stride+1):(i*y_stride);if length(offset_r)==1 % for extreme size like [1,1]offset_r = [offset_r, offset_r];endfor j = 1:floor(col/x_stride)offset_c = ((j-1)*x_stride+1):(j*x_stride);if length(offset_c)==1offset_c = [offset_c, offset_c];end接下來就是采用一種池化方法去進行池化操作，池化方法有很多種，比如均值池化，最大值池化，中值池化等，這里采用加和池化，將2*2塊中的所有元素加起來

block_val = squeeze(sum(sum(h_input(offset_r,offset_c,:,k))));注意一下這里控制住sum加和的個數(shù)，兩個sum函數(shù)的調(diào)用，剛好能保證對每張圖的每個特征圖分別加和，得到的block_val第一個維度是1(2*2的池化塊被加和的結果)，第二個維度是卷積核個數(shù)(池化后特征圖的數(shù)目是不變的)，第三個維度是1(當前選中的輸入圖片)，被squeeze以后去掉維度為1的維度，其實就是得到了9個值。然后再恢復2*2大小的特征，但是每個特征圖2*2塊里面的數(shù)據(jù)是一樣的，都是原來的2*2的加和。看看效果，我隨便輸出了其中一塊，改改代碼就能看到結果，調(diào)試的時候這樣改

block_val = squeeze(sum(sum(h_input(offset_r,offset_c,:,k))));block_valblock(offset_r,offset_c,:,k) = repmat(permute(block_val, [2,3,1]), numel(offset_r),numel(offset_c));block(offset_r,offset_c,:,k)然后便可以看到效果，窗口會彈出很多結果，這里就截取其中一部分：

block_val =3.99843.96214.01694.00073.98493.99793.96654.01753.9813ans(:,:,1) =3.9984 3.99843.9984 3.9984ans(:,:,2) =3.9621 3.96213.9621 3.9621ans(:,:,3) =4.0169 4.01694.0169 4.0169ans(:,:,4) =4.0007 4.00074.0007 4.0007ans(:,:,5) =3.9849 3.98493.9849 3.9849ans(:,:,6) =3.9979 3.99793.9979 3.9979ans(:,:,7) =3.9665 3.96653.9665 3.9665ans(:,:,8) =4.0175 4.01754.0175 4.0175ans(:,:,9) =3.9813 3.98133.9813 3.9813

至此，池化部分搞定。回到調(diào)用池化的函數(shù)crbm_inference2D中來，文章計算了一下一個分數(shù)，個人感覺這樣是計算了一種概率，與最后附錄介紹的最大概率池化方法很像，額。其實我想說好像就是概率最大池化方法，但是和最大概率池化方法最開始提出的作者的論文的代碼不一樣，所以說好像是。

if strcmp(layer.type_input, 'Binary')model.h_sample = exp(model.h_input)./(1+block); elsemodel.h_sample = exp(1.0/(model.start_gau^2).*model.h_input)./(1+block); end按照代碼可以發(fā)現(xiàn)對于每一個池化塊的隱層輸出應該是這樣計算，對于二值情況下，假設 代表其中一個特征圖中2*2的被池化塊block，那么新的 計算方法如下：

第六步

進入對比散度算法中的采樣階段。首先是利用隱層重構可見層crbm_reconstruct2D。

一開始就先計算一下隱層的激活狀態(tài)

h_state = double(rand(size(model.h_sample)) < model.h_sample);初始化可見單元v_input和記錄下當前批需要處理的圖片個數(shù)N。
model.v_input = zeros(size(model.v_input)); N = size(model.v_input,4);然后利用conv2進行反卷積，具體調(diào)用方法可以看我前面關于 conv2，filter2，imfilter區(qū)別的轉載，利用隱層狀態(tài)和權重來反推可見層數(shù)據(jù)

for k = 1:Nfor i = 1:layer.n_map_vfor j = 1:layer.n_map_hmodel.v_input(:,:,i,k) = model.v_input(:,:,i,k) + conv2(h_state(:,:,j,k),model.W(:,:,i,j),'full');endmodel.v_input(:,:,i,k) = model.v_input(:,:,i,k) + model.v_bias(i);end end三個循環(huán)，分別控制第幾張圖片的第幾個特征圖的重構是利用的第幾個卷積核。h_state的四個維度分別是特征圖大小*特征圖個數(shù)*輸入圖片個數(shù)，其中每一個通道的重構就是利用每個卷積核與當前特征圖連接權重的乘積的和來推理的，與彩色圖像卷積一樣，最后不要忘記了加上偏置項。這里再提示一下有心人有興趣可以看看偏置，誤差，方差的區(qū)別。

最后就是如果我們需要二值的可見單元，就用sigmoid函數(shù)激活一下，如果是高斯單元，就是直接將反卷積的輸出當做可見單元值即可

if strcmp(layer.type_input, 'Binary')model.v_sample = sigmoid(model.v_input);elsemodel.v_sample = model.v_input;end

第七步

就是反復進行吉布斯采樣，簡單點說就是不斷重復第五步、第六步，知道達到設置的k步對比散度算法的值。

至此我們能夠得到第一次進入模型時候的模型狀態(tài)和第k次采樣以后的模型狀態(tài)。利用這兩種狀態(tài)計算一下權重的更新梯度。需要注意的是計算的時候要先將隱單元的卷積核旋轉180°，原因是在positive phase階段，利用可見層推導隱藏層單元值的時候，對卷積核進行過旋轉。

dW = zeros(size(model.W)); for i = 1 : layer.n_map_hfor j = 1 : layer.n_map_vdW(:,:,j,i) = conv2(data(:,:,j,1),model.h_sample_init(end:-1:1,end:-1:1,i,1),'valid') ...- conv2(model.v_sample(:,:,j,1),model.h_sample(end:-1:1,end:-1:1,i,1),'valid');end end

第八步

更新模型參數(shù)，主要包含權重，可見層偏置，隱藏層偏置，對于高斯形式和二值形式都有不同的計算公式，先計算一下偏置更新梯度(權重更新梯度上一步計算了)

N = size(data,4); dV_bias = squeeze(sum(sum(sum(data - model.v_sample,1),2),4)); dH_bias = squeeze(sum(sum(sum(model.h_sample_init - model.h_sample,1),2),4));

可以發(fā)現(xiàn)，更新方法和RBM基本一模一樣。來一個小插曲，先看看matlab的sum函數(shù)是如何對一個高維矩陣進行工作的。

A(:,:,1)=[1 2 3;4 5 6] A(:,:,2)=[7 8 9;10 11 12] A(:,:,3)=[13 14 15;16 17 18]建立一個三維矩陣，先計算一次加和
sum(A)ans(:,:,1) =5 7 9ans(:,:,2) =17 19 21ans(:,:,3) =29 31 33

>> sum(A,2)ans(:,:,1) =615ans(:,:,2) =2433ans(:,:,3) =4251 >> sum(A,3)ans =21 24 2730 33 36

可以發(fā)現(xiàn)這個sum函數(shù)與第二個參數(shù)的選擇是有關的。如果是1，那么就行+行；如果3，列就列+列，如果是3，就第三維+第三維

回到程序中偏置的更新，都是對四維矩陣加和，對dV_bias分析一下(對于dH_bias的分析是一樣的)：

假設輸入的data是256*256*3*10，意思是輸入了10張圖片，每張圖片大小是256*256，且為三通道RGB彩色圖像

首先是計算一下原始圖與重構圖的差值，并進行加和，得到的矩陣應該是256*1*3*10：

sum(data - model.v_sample,1)然后進行第二個加和處理，應該是得到一個1*1*3*10大小的矩陣：

sum(sum(data - model.v_sample,1),2)繼續(xù)進行第三個加和處理：

sum(sum(sum(data - model.v_sample,1),2),4)這里的維度變成第四維了，為什么不是第三維呢？回想前面的理解中說過每一個特征圖都有一個偏置項，且每個偏置的值是1*1大小的共享偏置，也就是說有幾個特征圖就有幾個1*1的偏置值。 那么對應到輸入的彩色圖像中來，只有三個通道(RGB)對應有偏置值，而不是對每一個樣本的每個通道都有一個不同的偏置。因此進行加和的維度分別是第一維、第二維、第四維，分別代表輸入圖像大小，每次批處理的輸入圖像的個數(shù)。

接下來看看對于二值情況的更新方法：

if strcmp(layer.type_input, 'Binary')dW = dW/N;dH_bias = dH_bias/N;dV_bias = dV_bias/N;model.dW = model.momentum*model.dW + (1-model.momentum)*dW;model.W = model.W + layer.learning_rate*(model.dW - layer.lambda2*model.W);penalty = 0;model.dV_bias = model.momentum*model.dV_bias + (1-model.momentum)*dV_bias;model.v_bias = model.v_bias + layer.learning_rate*(model.dV_bias - penalty*model.v_bias);model.dH_bias = model.momentum*model.dH_bias + (1-model.momentum)*dH_bias;model.h_bias = model.h_bias + layer.learning_rate*(model.dH_bias - penalty*model.h_bias);model.h_bias = model.h_bias + layer.learning_rate*layer.lambda1*...(squeeze(layer.sparsity-mean(mean(mean(model.h_sample_init,1),2),4))); end處理步驟就是：

①計算完畢加和以后，除以每批樣本總數(shù)

②針對權重和偏置需要更新的梯度都是利用“動量項*①中的梯度+(1-動量項)*①中的梯度”進行計算的
③然后計算最終投入到下次迭代的模型參數(shù)：

首先是權重(有梯度懲罰項)，計算方法是：

下次迭代需要的權重=當前迭代使用的權重+學習率*(上一步計算的權重梯度 - 梯度懲罰率*當前迭代使用的權重)

然后是可見層偏置(簡單計算)，計算方法是：

下次迭代的可見層偏置=當前迭代使用的可見層偏置+學習率*(上一步計算的可見偏置梯度)

最后是隱藏層偏置(有稀疏項)，計算方法是：

中間變量=當前迭代使用的隱藏層偏置+學習率*(上一步計算的隱層偏置梯度)

下次迭代的隱藏層偏置=中間變量+學習率*稀疏懲罰*(稀疏項-原始輸入重構得到的隱層單元值的均值)

最后看看對于高斯輸入的梯度更新

if strcmp(layer.type_input, 'Gaussian')N = size(model.h_sample_init,1) * size(model.h_sample_init,2) * layer.batchsize;dW = dW/N - 0.01*model.W;%dh = (squeeze(sum(sum(model.h_sample_init,1),2)) - squeeze(sum(sum(model.h_sample,1),2)))/N - 0.05*(squeeze(mean(mean(model.h_sample_init,1),2)) - 0.002);dh = (squeeze(sum(sum(sum(model.h_sample_init,1),2),4)) - squeeze(sum(sum(sum(model.h_sample,1),2),4)))/N;dv = 0;model.winc = model.winc*model.momentum + layer.learning_rate*dW;model.W = model.winc + model.W;model.hinc = model.hinc*model.momentum + layer.learning_rate*dh;model.h_bias = model.hinc + model.h_bias;model.vinc = model.vinc*model.momentum + layer.learning_rate*dv;model.v_bias = model.vinc + model.v_bias; end分母除的不再是每批輸入樣本的個數(shù)了，而是采用另一種計算方法

N= size(model.h_sample_init,1) * size(model.h_sample_init,2) * layer.batchsize;代表的是隱單元每個特征圖單元數(shù)*每批的輸入數(shù)據(jù)數(shù)目

權重的更新也同樣是先計算被懲罰以后的梯度，然后直接計算被更新的最終權重

dW= dW/N - 0.01*model.W; model.winc= model.winc*model.momentum + layer.learning_rate*dW; model.W= model.winc + model.W;隱藏層偏置同樣也是計算對比散度算法前后的隱單元差值，感覺和上面一樣，沒什么變化，但是換了一種寫法，以及分母變化了而已

dh= (squeeze(sum(sum(sum(model.h_sample_init,1),2),4)) - squeeze(sum(sum(sum(model.h_sample,1),2),4)))/N; model.hinc= model.hinc*model.momentum + layer.learning_rate*dh; model.h_bias= model.hinc + model.h_bias;對于可見層偏置的梯度，直接置零了，用學習率多余計算了一下，可以刪掉學習率這一部分，直接把需要更新的梯度設置為動量項倍即可。
dv = 0; model.vinc = model.vinc*model.momentum + layer.learning_rate*dv; model.v_bias = model.vinc + model.v_bias;

第九步

輸出每次迭代的信息看看，包括稀疏程度、重構誤差

sparsity(i) = mean(model.h_sample_init(:)); err1= (batch_data - model.v_sample).^2; err= err + sum(sum(sum(sum(err1))));后面那個gau是控制如果隱層不是二值單元而是高斯單元的時候，進行隱層計算時候需要的高斯項，具體調(diào)用地方查看crbm_inference2D最好一行代碼。這里主要控制高斯項逐漸遞減，因為初始化時設置的(查看default_layer2D函數(shù))

layer.start_gau = 0.2; % GAUSSIAN START layer.stop_gau = 0.1; % GAUSSIAN END所以每次迭代完畢遞減一下：

if (model.start_gau > model.stop_gau)model.start_gau = model.start_gau*0.99;end

第十步

完成第一層卷積RBM的參數(shù)更新以后，計算一下進入下一層卷積RBM的輸入值。同樣只關注"matlab"部分，利用函數(shù)crbm_forward2D

很簡單就是，與crbm_inference2D處理方法類似，都是利用positive phase的計算流程，唯一不同的是，這里需要考慮到批數(shù)據(jù)的批數(shù)，而在crbm_inference2D中每次迭代已經(jīng)選定批數(shù)了。

先初始化相關參數(shù)，記錄一下所有訓練集的數(shù)量大小，當前層池化大小，以及隱層的特征大小。

n = size(data,4); x_stride = layer.s_pool(2); y_stride = layer.s_pool(1); row = size(model.h_sample,1); col = size(model.h_sample,2);然后由于上面介紹過隱層每個單元2*2，也就是被池化塊的部分元素都是一樣的，因為采用的是加和池化方法。所以輸入到下一層CRBM的數(shù)據(jù)大小是當前層隱層大小/池化大小，即：

output = zeros(floor(row/y_stride),floor(col/x_stride),layer.n_map_h,n);對應的分別是輸入大小*特征圖個數(shù)*原始圖片總數(shù)

初始化完畢以后就需要進行計算了，按照彩色圖像的計算方法，對每個通道卷積，然后求和

for i = 1:layer.n_map_h for j =1:layer.n_map_vmodel.h_input(:,:,i,1) = model.h_input(:,:,i,1) + conv2(batch_data(:,:,j,1),model.W(end:-1:1,end:-1:1,j,i),'valid'); endmodel.h_input(:,:,i,1) = model.h_input(:,:,i,1) + model.h_bias(i); end然后利用求和方法計算池化單元的值，利用公式計算單元值，但是在每個池化塊的重復元素中取一個，采用的方法是間隔采樣。

block = crbm_blocksum2D(model,layer); h_sample = 1-(1./(1+block)); output(:,:,:,k) = h_sample(1:y_stride:row, 1:x_stride:col,:);這里的output就是即將丟入到下一層訓練的輸入層的數(shù)值。

再后面就是DemoCDBN_Gaussian_2D.m主函數(shù)對每一層得到的特征圖可視化的方法啦，與卷積RBM關系不大，就不寫咯。

其實個人感覺這個第十步的核心部分(隱單元計算和池化層計算就是下面附錄介紹的最大池化方法)

附錄概率最大池化的理解

關于大牛Honglak Lee 提出的論文主要采用了一種概率最大池化的方法，表示看了代碼以后和原文可能有點對不上號，按照代碼說一下吧

主要函數(shù)就是multrand2.m，比較短，可以去開頭提供的地址下載，我這里也直接粘貼出來

function [S P] = multrand2(P) % P is 2-d matrix: 2nd dimension is # of choices概率最大池化的方法 %P是除以均值了 %S是進行激活以后的值，為0或者1 % sumP = row_sum(P); sumP = sum(P,2); P = P./repmat(sumP, [1,size(P,2)]);cumP = cumsum(P,2);%計算各行累加值 % rand(size(P)); unifrnd = rand(size(P,1),1); temp = cumP > repmat(unifrnd,[1,size(P,2)]);%輸出0,1 %發(fā)現(xiàn)很神奇的地方，Sindx每一行4個單元至多有一個為1 Sindx = diff(temp,1,2);%對列求一階偏導diff(A,m,n)對A按照n(1為列，2為行)求m階導數(shù) S = zeros(size(P)); S(:,1) = 1-sum(Sindx,2);%sum(A,2)是對行求和，得到一列數(shù)據(jù) S(:,2:end) = Sindx; %最終得到的S矩陣每一行有且僅有一個1 end 簡要摘取一下文章中關于概率最大池化(probabilistic max-pooling)方法的介紹吧。
一般來說高層特征的檢測需要逐步包含逐步增大的輸入?yún)^(qū)域的信息。現(xiàn)有的轉移不變性表示方法比如卷積神經(jīng)網(wǎng)絡經(jīng)常包含兩個步驟：

“檢測層“”對之前層進行卷積得到一個特征檢測器

“池化層”通過常數(shù)因子縮小檢測層的表示。

具體來說就是，每個池化層單元計算的是檢測層的一個小區(qū)域的最大激活值。利用最大池化方法可以使得高層的特征表示對于輸入的小變動具有不變性，并且降低了計算負擔。

最大池化一般用于前饋網(wǎng)絡。相反這里對一個圖片的生成模型比較感興趣，支持從上往下和從下往上的推斷。因此設計了一個生成模型，利用類似最大池化的推斷方法。

假設可見層是V，檢測層是H，池化層是P。那么H和P都有k個單元組，每個池化單元組大小是大小。

探測層的被池化塊和對應的池化單元有一個潛在約束：探測層至多有一個單元被打開(激活)，池化單元當且僅當一個探測單元激活的時候被激活。

另一種等價說法就是，我們可以考慮有個單元(被池化單元+對應的池化單元)作為一個獨立隨機變量：一個值代表所有探測單元被開啟，還有一個值代表所有的單元被關閉。

簡單地創(chuàng)建最大池化CRBM的能量函數(shù)就是：

約束條件里面的代表的是在探測層H中被池化的一個小塊，至多只有一個被激活為1。

探測層的第k組接受的輸入是從V由下往上的輸入信號，其實就是v乘以權重加上隱層偏置，沒啥好理解的：

然后就是對這個隱層單元分塊進行池化，假設是一個包含在塊α中的隱單元，條件概率就是這樣計算的

第一個公式就是利用可見層得到隱層的激活概率，第二個是計算對應的池化單元被抑制的概率。其實這個第二個公式很好理解，按照可見層被池化塊的幾個單元有任意一個單元沒激活的概率來計算。也就是說，第一個公式計算了其中一個單元被激活的概率，那么整個被池化塊的探測層單元至少有一個被激活的概率就是

這樣就計算出來了池化層被激活的概率，看看這個式子和上面第二個式子的和是等于1吧，代表被激活的概率加上被抑制的概率的和是1，因為只存在被激活和被抑制兩種狀態(tài)嘛。

然后看看代碼，發(fā)現(xiàn)，好像并沒有第一個公式提到的分母加1的情況。代碼先將隱單元均值化，然后依次計算累加值，這樣最后一行值就是1了，

sumP = sum(P,2); P = P./repmat(sumP, [1,size(P,2)]);cumP = cumsum(P,2);%計算各行累加值

這個累加函數(shù)cumsum類似于這樣

>> B=[1 2 3;4 5 6]B =1 2 34 5 6>> cumsum(B,2)ans =1 3 64 9 15>>

回到cump就會發(fā)現(xiàn)，每行的值都是依次增大的，然后在每一行利用一個概率去激活，可以得到在某個位置左邊的值全是0，右邊全是1，因為數(shù)值是依次增大。然后計算梯度，剛好在分界處得到一個1，其它地方都是0。

還是利用上面這個cumsum的例子，假設我們對結果的第一行用閾值2激活，第二行用閾值10激活，得到的結果就是

0 1 1 0 0 1接下來計算每行的梯度就是，而且總能確保每行至多只有一個值為1
1 0 0 1然后計算將每行加和起來，去判斷到底這一行是有一個1還是全0，對于全0的情況。

最后問題又來了，按照文章描述，隱層至少有一個單元為1的時候，池化單元才被激活，也就是說在這個包含C*C的被池化部分和對應的池化單元的個單元里面，要么有兩個1，要么全0，然而按照代碼

S(:,1) = 1-sum(Sindx,2);%sum(A,2)是對行求和，得到一列數(shù)據(jù) S(:,2:end) = Sindx;可以發(fā)現(xiàn)這兩句話控制了每行有且僅有一個1，與原文的限制條件有點區(qū)別。

用一個實例來體現(xiàn)這個最大池化函數(shù)multrand2的效用

>> A=rand(10,5)A =0.9604 0.3956 0.2853 0.3660 0.06910.7845 0.4872 0.2348 0.0293 0.32410.0900 0.3473 0.8310 0.2054 0.80440.5143 0.5906 0.9595 0.4494 0.64290.3923 0.7827 0.4778 0.3040 0.47320.9238 0.5114 0.5506 0.7581 0.00270.2364 0.5663 0.1936 0.4560 0.07490.2861 0.0362 0.5895 0.5049 0.29970.2368 0.1295 0.7861 0.6929 0.88480.3697 0.6487 0.8074 0.0170 0.7988>> S=multrand2(A)S =1 0 0 0 00 1 0 0 00 0 0 1 00 0 0 0 10 1 0 0 01 0 0 0 00 0 0 1 00 0 0 1 00 0 1 0 00 0 0 0 1>> 【結論】講道理的話，上面第十步的處理方法與最大概率池化的公式最近似，基本一樣。而原作者提供的代碼與論文描述稍微有點區(qū)別。讀者若有任何見解，請在評論區(qū)寫出，大家一起討論討論，謝謝。額，脖子、脖子，疼啊o(╯□╰)o

總結

以上是生活随笔為你收集整理的卷积RBM源码解读的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。