[matDL框架開發直播:2]全連接層(dense)的實現和優化

01-24

=====================================

往期:

用微信控制深度學習訓練：中國特色的keras插件 - 知乎專欄

[matDL框架開發直播:0]深度學習水平考試：手寫一個輕量級框架! - 知乎專欄

[matDL框架開發直播:1]matDL概述和基本使用 - 知乎專欄

=====================================

GitHub地址：QuantumLiu/matDL

在上一期我們了解了matDL的概況並搭建、訓練了第一個模型。

像我們第一個搭建的這種模型，每一層都是一個全連接層（full connected layer/FC/dense）.

全連接層是最基本、最基本的神經網路/深度神經網路層。

基礎知識

全連接層前向傳播的過程就是每一個樣本的輸入數據向量（矩陣）x，與每一個層的權值（weights）矩陣W做矩陣乘法（Hadamard product (matrices)）

對於一個a=p*m的矩陣和一個b=m*q的矩陣，他們的乘積c是一個p*q的矩陣。

$y=x*W+b$

特別的，代表一個樣本的單個向量可以視為矩陣列數為1時的特例。一般而言，現在深度學習普遍採用mini-batch方法，x往往是形如[input dim,minibatch]的矩陣張量而非單個向量。

在matDL中，輸入張量x形如[input dim,minibatch] 或 [input dim,timestep,minibatch]，權值矩陣W形如[output dim,input dim] 為了方便大家理解，我用win10畫圖畫了個靈魂圖解（敢說丑砍死你）

b代表bias是偏置項

反向傳播實際上是求x和W對於輸出y的偏導數

我們看，對於一個 $y_{i}=sum_{1}^{j}{w_{ij}*x_{j}}$ , $x_{j}$ 的導數就是 $W_{ij}$

我們對於y進行反傳求導得到的也是矩陣即dx和dw，形狀和x與w一致。

所以實際上，在得到y的梯度dy後,dW和dx可以通過

$dx=W^{T} *dy$

$dW=dy*x^{T}$ 輕鬆地求出來

如果我有沒有講清楚的地方歡迎指出或者提問~

matDL實現

既然全連接層本質就是一個矩陣乘法，那麼用MATLAB寫起了就非常簡單了。

創建層（init）

function layer=dense_init_gpu(prelayer,hiddensize ,flag,loss)%% Basic layer attributes%Input tensor sahpelayer.trainable=1;layer.flag=flag;if numel(prelayer.output_shape)>2 layer.timedistributed=1;else layer.timedistributed=0;endlayer.input_shape=prelayer.output_shape;dim=prelayer.output_shape(1);batchsize=prelayer.output_shape(end);layer.type="dense";layer.prelayer_type=prelayer.type;layer.output_shape=[hiddensize,layer.input_shape(2:end)];layer.hiddensize=hiddensize;layer.batchsize=batchsize;layer.batch=1;layer.epoch=1;%% Dense layer attributes%W contains weights biaslayer.weights_dim=dim+1;layer.W=(rand([hiddensize,layer.weights_dim],"single","gpuArray")-0.5)./100;layer.input=ones([layer.input_shape(1)+1,layer.input_shape(2:end)],"single","gpuArray");layer.output=zeros(layer.output_shape,"single","gpuArray");if ~strcmpi(layer.prelayer_type,"input")&&flag layer.dx=zeros(layer.input_shape,"single","gpuArray");endlayer.e=layer.output;if nargin>3&&flag [layer.loss_f,layer.loss_df]=loss_handle(loss); layer.loss=[];endlayer.ff=@(layer,prelayer)dense_ff_gpu(layer,prelayer);layer.bp=@(layer,next_layer)dense_bp_gpu(layer,next_layer);layer.configs.type=layer.type;layer.configs.input_shape=layer.input_shape;layer.configs.output_shape=layer.output_shape;layer.configs.hiddensize=layer.hiddensize;layer.configs.W=size(layer.W);end

第一個形參prelayer是所有matDL層共有的參數，代表前一個層，用來鏈式定義每一層的input。

我們重點來看這幾句

%W contains weights biaslayer.weights_dim=dim+1;layer.W=(rand([hiddensize,layer.weights_dim],"single","gpuArray")-0.5)./100;layer.input=ones([layer.input_shape(1)+1,layer.input_shape(2:end)],"single","gpuArray");

layer.weights_dim=dim+1;

正如之前所說，對於每一個輸出，還要加一個偏置項bias。如果單出來一個長度為hiddensize的向量來儲存bias會導致計算複雜。我們不妨把偏置項看做輸入數據多出來一個維度，這個維度所以數據都是常數1， $x_{bias}=1$ ，和其他數數據一樣在權值矩陣中有自己對應的 $W_{i,bias}$ 。於是我們的input多了一個全1維度，權值矩陣也多了一列bias參數。

由於MATLAB的面向對象不是很好，matDL中用struct類型來實現「偽」oop。

前向計算（feedfoward）

function layer=dense_ff_gpu(layer,prelayer)if isequal(class(prelayer),"struct") if ~isequal(size(prelayer.output),layer.input_shape) error("Shape unmatched!") end if layer.timedistributed layer.input(1:end-1,:,:)=prelayer.output; else layer.input(1:end-1,:)=prelayer.output; endelse if layer.timedistributed layer.input(1:end-1,:,:)=prelayer; else layer.input(1:end-1,:)=prelayer; endendif layer.timedistributed layer.output(:)=layer.W*sq(layer.input);else layer.output=layer.W*layer.input;endendfunction a=sq(a)a=reshape(a,size(a,1),[]);end

可以看到，我額外定義了一個函數sq(),這個函數將一個形如[input dim,timestep,batchsize]的張量展開為一個形如[input dim,timestep*batchsize]的向量。

之所以需要這個函數是為了實現TimeDistributed Dense功能，即對於每一個時間步的數據產生一個獨立的輸出，可參考包裝器Wrapper - Keras中文文檔。

我們注意到

layer.input(1:end-1,:)=prelayer;

這就對應著前面定義dense層時，layer.input矩陣多出來一維全1向量來作為偏置。

dense層的ff方法和其它層一致，除了層本身有一個形參prelayer（上一層），用來接收輸入張量。

dense層的輸出賦值給了layer.output

反向傳播（backpropagation）

function layer =dense_bp_gpu(layer,next_layer)if isequal(class(next_layer),"struct") if ~isequal(size(next_layer.dx),layer.output_shape) error("Shape unmatched!") end layer.e=next_layer.dx;endlayer.dW=sq(layer.e)*sq(layer.input)";if ~isequal(layer.prelayer_type,"input") if layer.timedistributed layer.dx(:)=layer.W(:,1:end-1)"*sq(layer.e); else layer.dx=layer.W(:,1:end-1)"*layer.e; endendendfunction a=sq(a)a=reshape(a,size(a,1),[]);end

dense層的bp方法和其它層一樣，有layer本身和next_layer兩個形參。next_layer即網路連接的下一層，本層的輸出張量是下一層的輸入張量。

我們使用

isequal(class(next_layer),"struct")

來判斷本層是否是最後一層（outputlayer）,對於最後一層，本層輸出對於loss的梯度layer.e將在eval_loss()時被賦值，否則獲得前一層的輸入（即本層輸出）的梯度dx的值。

如果本層不是首層

~isequal(layer.prelayer_type,"input")

那麼我們將計算本層的dx值傳給上一層。

於是，我們得到了權值矩陣的梯度dW，接下來優化器（optimizer）會根據dW和策略，對W進行參數更新，至此，我們完成了全連接層在一個batch數據上的迭代過程！

優化

全連接層的核心矩陣乘法使用的是MATLAB的gpuArray的built-in函數，效率極高，我寫了一個.cu的mex文件，希望通過調用NVIDIA爸爸的cublas庫來優化速度。

#include <cuda_runtime.h>#include <cublas_v2.h>#include <time.h>#include "mex.h"#include "gpu/mxGPUArray.h"typedef struct _matrixSize{ unsigned int uiWA, uiHA, uiWB, uiHB, uiWC, uiHC;} sMatrixSize;void matrixMultiply(float const* d_A, float const* d_B, float* d_C, sMatrixSize &matrix_size);void mexFunction(int nlhs, mxArray *plhs[], int nrhs, mxArray const *prhs[]){ mxGPUArray const *A; mxGPUArray const *B; mxGPUArray *C; mxInitGPU(); _matrixSize matrix_size; mwSize const *A_sz; mwSize const *B_sz; float const *d_A; float const *d_B; float *d_C; char const * const errId = "parallel:gpu:mexGPUExample:InvalidInput"; char const * const errMsg = "Invalid input to MEX file."; if (nrhs != 2) { mexErrMsgTxt("Need two inputs"); } A = mxGPUCreateFromMxArray(prhs[0]); B = mxGPUCreateFromMxArray(prhs[1]); A_sz=mxGPUGetDimensions(A); B_sz = mxGPUGetDimensions(B); matrix_size.uiWA = (unsigned int)A_sz[0]; matrix_size.uiHA = (unsigned int)A_sz[1]; matrix_size.uiWB = (unsigned int)B_sz[0]; matrix_size.uiHB = (unsigned int)B_sz[1]; mwSize C_sz[3] = { matrix_size.uiWA, matrix_size.uiHB, 1 }; d_A = (float const *)(mxGPUGetDataReadOnly(A)); d_B = (float const *)(mxGPUGetDataReadOnly(B)); C = mxGPUCreateGPUArray(mxGPUGetNumberOfDimensions(A), C_sz, mxGPUGetClassID(A), mxGPUGetComplexity(A), MX_GPU_DO_NOT_INITIALIZE); d_C = (float *)(mxGPUGetData(C)); matrixMultiply(d_A, d_B, d_C, matrix_size); plhs[0] = mxGPUCreateMxArrayOnGPU(C); mxFree((void*)A_sz); mxFree((void*)B_sz); mxGPUDestroyGPUArray(A); mxGPUDestroyGPUArray(B); mxGPUDestroyGPUArray(C);}void matrixMultiply(float const* d_A, float const* d_B, float* d_C, sMatrixSize &matrix_size){ cublasStatus_t status; cublasHandle_t handle; status=cublasCreate(&handle); if (status != CUBLAS_STATUS_SUCCESS) { if (status == CUBLAS_STATUS_NOT_INITIALIZED) { mexPrintf("CUBLAS initializing error"); } getchar(); return; } const float alpha = 1.0f; const float beta = 0.0f; cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, matrix_size.uiWA, matrix_size.uiHB, matrix_size.uiHA, &alpha, d_A, matrix_size.uiWA, d_B, matrix_size.uiWB, &beta, d_C, matrix_size.uiWA); cudaThreadSynchronize(); cublasDestroy(handle);}

然後用mexcuda編譯

mexcuda -lcublas cublasSegmm.cu -v COMPFLAGS="-arch=compute_50 -code=sm_50 -O2"

測試一下

cuseg=@cublasSegmm;a=ones(2048,2048,"single","gpuArray");b=ones(2048,2048,"single","gpuArray");tic;for i=1:100c_cu=cuseg(a,b);end;toctic;for i=1:100c_in=a*b;end;toc