YOLOv1解析（二）

10-06

YOLOv1解析（二）

來自專欄目標檢測

轉載請註明出處！

YOLOv1解析（二）

————————

上一篇YOLOv1解析（一）討論了前向傳播，本篇繼續討論後向傳播。

後向傳播

從最後一層依次計算，直到第一層，見如下代碼，

network net = *netp; // netp: network parameterint i;network orig = net;for(i=net.n-1;i>=0;--i){ layer l = net.layers[i]; // 當前layer if(l.stopbackward) break; if(i==0){ net=orig; }else{ layer prev = net.layers[i-1]; // 上一層layer net.input=prev.output; // 上一層layer的output作為當前layer的input net.delta=prev.delta; // 用於存儲上一層layer的delta } net.index=i; l.backward(l,net);}

detection

從上一篇文章YOLOv1解析（一）最後的代碼片段中，我們知道目標函數L對network output的 $hat x_i, hat y_i, hat w_i, hat h_i$ （將這些變數看作本層輸出 $a_i$ ）的梯度存儲在l.delta中（即， $partial L / partial a_i$ )。detection層的backward傳播，僅僅將delta數據複製到net.delta中，也就是上一層的delta，畢竟，detection層就是直接計算損失，沒有對上一層的output作任何處理。

需要注意的是代碼中計算梯度均少了一個負號，所以在參數更新時本該使用梯度下降的全部改為梯度上升。比如損失對輸出值 $hat x_i$ 的梯度為，

$frac {partial L} {partial hat x_i} = -2(x_i - hat x_i)^2$

而代碼中的delta計算式為，

l.delta[box_index+0] = l.coord_scale*(net.truth[tbox_index + 0] - l.output[box_index + 0]);

發現少了一個負號。這裡我們保持與代碼一致，去掉這個負號，並切記參數更新時使用梯度上升，下文的梯度計算公式均不再保留負號。

2. connected

本層backward與forward處理順序相反，先處理激活函數，激活函數是線性激活函數，已知forward傳播公式為，

$z_i=sum_{j=1}^n W_{ij} cdot x_j+b_i quad (1) \ a_i=f(z_i) quad quad quad (2)$

其中下標i 表示本層第 i 個輸出單元，n表示本層input數量， $x_j$ 表示 j-th 輸入。

損失L對本層輸出值的梯度存儲在l.delta中。

根據導數的鏈式法則，

$frac {partial L} {partial z_i}=frac {partial L} {partial a_i} cdot frac {partial a_i} {partial z_i} quad (3)$

上式右邊第一項 $partial L / partial a_i$ 即本層的delta，乘以激活函數的梯度（導數） $partial a_i / partial z_i$ ，代碼如下

// 計算 dL/dz_i，並存儲在delta[i]處delta[i] *= gradient(x[i], a);// a指明激活函數類型

然後計算損失L對本層參數（權重和偏置）的梯度（導數）。根據全連接層的前向傳播公式，有

$frac {partial L} {partial b_i}=frac {partial L} {partial z_i}frac {partial z_i} {partial b_i}=frac {partial L} {partial z_i} quad (4)\ frac {partial L} {partial W_{ij}}=frac {partial L} {partial z_i}frac {partial z_i} {partial W_{ij}}=frac {partial L} {partial z_i}x_j quad (5) \ frac {partial L} {partial x_j}=sum_{i=1}^m frac {partial L} {partial z_i}frac {partial z_i} {partial x_j}=sum_{i=1}^mfrac {partial L} {partial z_i} W_{ij} quad (6)$

根據上面三個公式，分別計算各梯度。

(a)對bias的梯度

void backward_bias(float *bias_updates, float *delta, int batch, int n, int size){ int i,b; for(b = 0; b < batch; ++b){ // b-th image of this batch for(i = 0; i < n; ++i){ // i-th unit of this layers output // connected layer的size為1 // sum_array: 從第一個參數所指位置開始累加size個值，並返回這個累加值 // 疊加此batch中所有image在本layer第i個unit處的bias梯度（其實就是delta[i]） bias_updates[i] += sum_array(delta+size*(i+b*n), size); } }}

(b) 對weight的梯度

/** * M: output number * N: input number * K: batch * ALPHA: 1 * A: l.delta * lda: output number * B: net.input * ldb: input number * C: weight_updates * ldc: input number * */void gemm_tn(int M, int N, int K, float ALPHA, float *A, int lda, float *B, int ldb, float *C, int ldc){ int i,j,k; #pragma omp parallel for for(i = 0; i < M; ++i){ // i-th 輸出unit for(k = 0; k < K; ++k){ // batch中K個image data所產生的Weight梯度累加 // 獲取第k個image data所產生的損失L對output第i個unit的梯度，即dL/dz_i register float A_PART = ALPHA*A[k*lda+i]; for(j = 0; j < N; ++j){ // j-th 輸入unit // i*ldc+j：Wij的梯度存儲位置 // k*ldb+j：第k個image 對應的x_j C[i*ldc+j] += A_PART*B[k*ldb+j]; // 累加batch中K個Wij } } }}

(c)對本層input的梯度，即對上一層output的梯度，將存儲到上一層的delta中

/** * M: batch * N: input number * K: output number * ALPHA: 1 * A: l.delta * lda: output number * B: l.weights * ldb: input number * C: net.delta * ldc: input number * */void gemm_nn(int M, int N, int K, float ALPHA, float *A, int lda, float *B, int ldb, float *C, int ldc){ int i,j,k; #pragma omp parallel for for(i = 0; i < M; ++i){ // i-th image of this batch for(k = 0; k < K; ++k){ // k-th output unit // 獲取第i個image對應的損失L對第k個output unit的梯度 dL/dz_k register float A_PART = ALPHA*A[i*lda+k]; for(j = 0; j < N; ++j){ // j-th input unit // i*ldc+j：i-th image在本層的第j個輸入unit所對應的梯度存儲位置 // B[k*ldb+j]：Wkj，第k個output unit與第j個input unit之間的權重 // 累加所有K個乘積 C[i*ldc+j] += A_PART*B[k*ldb+j]; } } }}

3. drop

由於drop layer forward 傳播公式為

$a_i=egin{cases} 0 & rand() < prob \ x_i cdot scale & ext{otherwise} end{cases}$

故損失對本層input的梯度（記住，即損失L對上一層output的梯度，將存儲到上一層delta中），

$frac {partial L} {partial x_i}=frac {partial L} {partial a_i}frac {partial a_i} {partial x_i}=egin{cases} frac {partial L} {partial a_i} cdot 0 & rand() < prob \ frac {partial L} {partial a_i} cdot scale & ext{otherwise} end{cases}$

代碼也很簡單如下所示，

for(i = 0; i < l.batch * l.inputs; ++i){ float r = l.rand[i]; if(r < l.probability) net.delta[i] = 0; else net.delta[i] *= l.scale;}

4. local

根據前向傳播公式

$a_i=f(z_i)=f(x_i*W_i+b_i) quad (7)$

其中f是激活函數，*表示卷積， $W_i$ 是卷積核參數，因為local層不共享卷積核參數，所以帶下標i， $x_i$ 是 $a_i$ 輸出點對應於原輸入平面上視窗數據。

反向傳播先是經過激活函數，比較簡單，不再贅述，得到 ${partial L} / {partial z_i}$ 依然存儲在l.delta中。上一篇文章已經講過local層與convolution層類似，只是卷積核參數在spatial位置上不共享。

(a) 對bias的梯度

$frac {partial L} {partial b_i}=frac {partial L} {partial z_i}frac {partial z_i} {partial b_i}=frac {partial L} {partial z_i} quad (8)$

所以l.delta中就是損失L對bias的梯度，拷貝到l.bias_updates中，代碼如下，

for(i = 0; i < l.batch; ++i){ // 累加batch組bias梯度 axpy_cpu(l.outputs, 1, l.delta + i*l.outputs, 1, l.bias_updates, 1);}

(b) 對weight的梯度

前向傳播過程如下圖示意，由於卷積核參數不共享，所以在output spatial (ow*oh)中每一點的n個卷積核均不同，於是形成ow*oh個卷積權重矩陣B，輸入數據重組（參考卷積）得到矩陣A，由於卷積參數不共享，所以這裡將矩陣A按列分為ow*oh個列向量，每個列向量對應一個權重矩陣B，兩者相乘後得到一個列向量，共得到ow*oh個列向量，為矩陣C。

圖1 local層卷積示意圖。i表示輸出平面某點位置

由於卷積參數不共享，所以反向傳播公式為，

$frac {partial L} {partial B_{i,j,k}}=frac {partial L} {partial C_{j,i}} frac {partial C_{j,i}}{partial B_{i,j,k}}=frac {partial L} {partial C_{j,i}} A_{k,i}$

$frac {partial L} {partial A_{k,i}}=frac {partial L} {partial C_{j,i}} frac {partial C_{j,i}}{partial A_{k,i}}=frac {partial L} {partial C_{j,i}} B_{i,j,k}$

從上式中可以看出沒有求和項。代碼就不再貼出來分析了，與下文convolution層的代碼類似（可以先看卷積層的反向傳播部分），同樣的，由於input數據重組後有很多相同的數據點，所以需要對相同數據點的梯度求和得到相應的原數據點的梯度，使用函數col2im_cpu。

5. convolution

反向傳播先是經過激活函數，比較簡單，不再贅述，反向傳播經過激活函數後的值 ${partial L} / {partial z_i}$ 依然存儲在l.delta中，代碼如下，

gradient_array(l.output, l.outputs*l.batch, l.activation, l.delta);

由於卷積層的batch_normalize都設置為1，所以與前向傳播類似，執行後向傳播對應的批規範操作，

if(l.batch_normalize){ backward_batchnorm_layer(l, net);}

結合上一篇文章中對於Batch Normalization的前向傳播公式，這裡給出反向傳播公式，

$frac {partial L} {partial gamma}=sum_{i=1}^m frac {partial L} {partial z_i} frac {partial z_i} {partial gamma}=sum_{i=1}^m frac {partial L} {partial z_i} hat x_i quad (9)\ frac {partial L} {partial eta}=sum_{i=1}^m frac {partial L} {partial z_i} frac {partial z_i} {partial eta}=sum_{i=1}^m frac {partial L} {partial z_i} quad (10)$

其中m為batch中所有image的同一個channel平面上數據點之和，也就是說計算得到的 $gamma, eta$ （的梯度）數量均等於卷積層output channel數。因為這兩個參數也是變數，代碼中將Batch Normalization的處理放在專門的一個layer即batchnorm層，於是可以將這兩個參數類比為其他層的權重和偏置參數，所以需要計算其梯度。

有了(9)和(10)式，代碼就容易懂了，

(a) 對 $eta$ (bias)的梯度

void backward_bias(float *bias_updates, float *delta, int batch, int n, int size){ int i,b; for(b = 0; b < batch; ++b){ // b-th image in this batch for(i = 0; i < n; ++i){ // i-th channel // accumulate the i-th bias gradient, using the i-th channel data of all images in this batch bias_updates[i] += sum_array(delta+size*(i+b*n), size); } }}

(b) 對縮放值 $gamma$ 的梯度

// 求解 gamma 梯度backward_scale_cpu(l.x_norm, l.delta, l.batch, l.out_c, l.out_w*l.out_h, l.scale_updates);

其中l.x_norm就是 $hat x_i$ ，參考前向傳播的相關代碼不難印證這一點。

void backward_scale_cpu(float *x_norm, float *delta, int batch, int n, int size, float *scale_updates){ int i,b,f; for(f = 0; f < n; ++f){ // f-th channel float sum = 0; for(b = 0; b < batch; ++b){ // b-th image in this batch for(i = 0; i < size; ++i){ // i-th point of a spatial panel int index = i + size*(f + n*b); // index of a data point // accumulate batch*size items to get gamma gradient for channel f sum += delta[index] * x_norm[index]; // refer to (9) } } scale_updates[f] += sum; }}

也就是上一層的delta。相關反向傳播公式如下，

$frac {partial L} {partial hat x_i}=frac {partial L} {partial z_i}frac {partial z_i} {partial hat x_i}=frac {partial L} {partial z_i} gamma quad (11)$

$egin{align}frac {partial L} {partial sigma^2} & =sum_{i=1}^m frac {partial L} {partial hat x_i} frac {partial hat x_i} {partial sigma^2} \ &=sum_{i=1}^m frac {partial L} {partial hat x_i} frac { partial [(x_i-mu)(sigma^2+epsilon)^{-1/2}]}{partial sigma^2} \&=-frac 1 2sum_{i=1}^m frac {partial L} {partial hat x_i}(x_i-mu)(sigma^2+epsilon)^{-3/2} \&=-frac 1 2 (sigma^2+epsilon)^{-3/2}sum_{i=1}^m frac {partial L} {partial hat x_i}(x_i-mu) end{align} quad (12)$

$egin{align} frac {partial L}{partial mu}&=sum_{i=1}^m frac {partial L} {partial hat x_i} frac {partial hat x_i} {partial mu} +frac {partial L} {partial sigma^2} frac {partial sigma^2 }{partial mu}\ &=sum_{i=1}^m frac {partial L} {partial hat x_i} frac {-1} {sqrt{sigma^2+epsilon}}+frac {partial L} {partial sigma^2}frac {2 sum_{i=1}^m(mu-x_i)}m\&=sum_{i=1}^m frac {partial L} {partial hat x_i} frac {-1} {sqrt{sigma^2+epsilon}}+frac {partial L} {partial sigma^2}cdot 0 \&=frac {-1} {sqrt{sigma^2+epsilon}} sum_{i=1}^m frac {partial L} {partial hat x_i} end{align} quad (13)$

$egin{align} frac {partial L} {partial x_i}&=frac {partial L}{partial hat x_i}frac {partial hat x_i}{partial x_i}+frac {partial L}{partial sigma^2}frac {partial sigma^2}{partial x_i}+frac {partial L}{partial mu}frac {partial mu}{partial x_i} \ &=frac {partial L}{partial hat x_i} frac 1{sqrt{sigma^2+epsilon}}+frac {partial L}{partial sigma^2}frac {2(x_i-mu)} m + frac {partial L}{partial mu} frac 1 m end{align} quad (14)$

其中，l.x存儲了 $x_i$ ，m表示batch中所有image的同一channel上的數據點之和，大小為 $batch*ow*oh$ ，其中ow和oh表示卷積層的輸出寬和高。

(11-14)式對應的代碼

scale_bias(l.delta, l.scales, l.batch, l.out_c, l.out_h*l.out_w); // (11)mean_delta_cpu(l.delta, l.variance, l.batch, l.out_c, l.out_w*l.out_h, l.mean_delta); // (13)variance_delta_cpu(l.x, l.delta, l.mean, l.variance, l.batch, l.out_c, l.out_w*l.out_h, // (12) l.variance_delta);normalize_delta_cpu(l.x, l.mean, l.variance, l.mean_delta, l.variance_delta, l.batch, // (14) l.out_c, l.out_w*l.out_h, l.delta);

(d) 反向傳播經過卷積操作

上面(c)實際上算是batchnorm層了，只不過這一層沒有獨立出來，而是作為convolution層中的一步操作，所以(c)中的 $x_i$ 可以看作這個未獨立出來的batchnorm層的輸入，也就是convolution卷積後的輸出（相當於下圖的 $C_{j,i}$ ），根據下圖，

圖2 卷積示意圖

其中C表示輸出矩陣，B表示權重矩陣 (W)，A表示重組後的輸入矩陣，n為卷積核數量，v表示單個卷積核volume，那麼反向傳播公式為

$frac {partial L} {partial B_{j,k}}=sum_{i=1}^{ow*oh} frac {partial L} {partial C_{j,i}} frac {partial C_{j,i}} {partial B_{j,k}}=sum_{i=1}^{ow*oh} frac {partial L} {partial C_{j,i}}A_{k,i} quad (15)$

其中 $partial L / partial C_{j,i}$ 就是(c)中最後計算得到的 $partial L / partial x_l$ 。

所以先對輸入進行重組得到矩陣A，重組方法與前向傳播的情況一樣（參考《卷積》），代碼如下，

int i, j;int m = l.n/l.groups; // 單個int n = l.size*l.size*l.c/l.groups; // 分組卷積，卷積核沿channel方向切分成多個partint k = l.out_w*l.out_h; // output spatial 大小for(i = 0; i < l.batch; ++i){ for(j = 0; j < l.groups; ++j){ // 分組卷積 float *a = l.delta + (i*l.groups + j)*m*k; // dL/dC_ji float *b = net.workspace; // 存儲重組後的input data float *c = l.weight_updates + j*l.nweights/l.groups; // 存儲weight 梯度 float *im = net.input + (i*l.groups + j)*l.c/l.groups*l.h*l.w; float *imd = net.delta + (i*l.groups + j)*l.c/l.groups*l.h*l.w; if(l.size == 1){ b = im; } else { im2col_cpu(im, l.c/l.groups, l.h, l.w, l.size, l.stride, l.pad, b); // input data 重組 } gemm(0,1,m,n,k,1,a,k,b,k,1,c,n); // 調用gemm_nt

gemm_nt函數定義如下，

/** * M * N: 分組後的單個卷積核part的volume * A: dL/dC_ji * lda:output spatial size, 也就是圖2中的ow*oh * B: input data * ldb: 同lda * C: weight gradient **/void gemm_nt(int M, int N, int K, float ALPHA, float *A, int lda, float *B, int ldb, float *C, int ldc){ int i,j,k; #pragma omp parallel for for(i = 0; i < M; ++i){ for(j = 0; j < N; ++j){ register float sum = 0; for(k = 0; k < K; ++k){ // output spatial size. 參考圖2中的ow*oh sum += ALPHA*A[i*lda+k]*B[j*ldb + k]; // (15)，K個值求和 } // i*ldc+j：weight gradient 位置下標。發現此表達式不含輸出spatial信息(k,K) C[i*ldc+j] += sum; } }}

對輸入的梯度公式為，

$frac {partial L} {partial A_{k,i}}=sum_{j=1}^{n} frac {partial L} {partial C_{j,i}} frac {partial C_{j,i}} {partial A_{k,i}}=sum_{j=1}^{n} frac {partial L} {partial C_{j,i}}B_{j,k} quad (16)$

代碼部分為

if (net.delta) { a = l.weights + j*l.nweights/l.groups; // (16)中權重 B b = l.delta + (i*l.groups + j)*m*k; // (16)中dL/dC_ji c = net.workspace; // (16)中dL/dA_ki if (l.size == 1) { c = imd; } gemm(1,0,n,k,m,1,a,n,b,k,0,c,k); // 調用gemm_tn if (l.size != 1) { col2im_cpu(net.workspace, l.c/l.groups, l.h, l.w, l.size, l.stride, l.pad, imd); }}

其中gemm_tn函數定義如下，

/** * M:分組後的單個卷積核part的volume * N: output spatial size (ow*oh) * K: 分組後的每組卷積核數量void gemm_tn(int M, int N, int K, float ALPHA, float *A, int lda, float *B, int ldb, float *C, int ldc){ int i,j,k; #pragma omp parallel for for(i = 0; i < M; ++i){ for(k = 0; k < K; ++k){ register float A_PART = ALPHA*A[k*lda+i]; for(j = 0; j < N; ++j){ // 對K個值求和。參考(16)式，也是對n個卷積核求和 C[i*ldc+j] += A_PART*B[k*ldb+j]; } } }}

由於A是重組後的輸入矩陣，其中有重複的原始輸入點，所以還需要還原到原來的形狀，且對重複的數據點的梯度求和，就得到相應原始輸入點的梯度。代碼如下，

col2im_cpu(net.workspace, l.c/l.groups, l.h, l.w, l.size, l.stride, l.pad, imd);

6. maxpool

池化層沒有激活函數，其前向傳播，是取視窗中的最大值作為輸出值，由於前向傳播時，記錄了每個視窗的最大值位置（保存在l.indexes中），所以反向傳播，只有最大值位置處的梯度，其他位置梯度均為0。這一部分由於比較簡單，略。