Caffe 代码阅读

Caffe代码阅读—卷积部分

Caffe卷积部分首先是将图像转换成一个列矩阵，然后把相应的卷积操作转换成了矩阵乘法。
CNN操作分为两步份，一部分是Forward操作，另一个是BP过程，要看懂BP首先得了解一下普通neural network的求导过程，还有在看一下《Notes on convolution networks 》看一下对于卷积网络的求导过程。首先看一下forward过程。
卷积网络是在文件”./src/caffe/layers/conv_layer.cpp”里面实现的。

template <typename Dtype>
void ConvolutionLayer<Dtype>::Forward_cpu(const vector<Blob<Dtype>*>& bottom, const vector<Blob<Dtype>*>& top) {
  const Dtype* weight = this->blobs_[0]->cpu_data();
  for (int i = 0; i < bottom.size(); ++i) {
    const Dtype* bottom_data = bottom[i]->cpu_data();
    Dtype* top_data = top[i]->mutable_cpu_data();
    for (int n = 0; n < this->num_; ++n) {
      this->forward_cpu_gemm(bottom_data + n * this->bottom_dim_, weight,
          top_data + n * this->top_dim_);
      if (this->bias_term_) {
        const Dtype* bias = this->blobs_[1]->cpu_data();
        this->forward_cpu_bias(top_data + n * this->top_dim_, bias);
      }
    }
  }
}

bottom是前一层数据，top是下一层数据。数据流bottom->top。cpu_data()返回blon的数据指针，和mutable_cpu_date()不一样的是，mutable_cpu_data()可以修改，cpu_data()返回是是const指针。然后具体运算在forward_cpu_data()里面，bottom和权值相乘，得到top的数据。
forward_cpu_gemm()在base_conv_layer.cpp文件里面。

template <typename Dtype>
void BaseConvolutionLayer<Dtype>::forward_cpu_gemm(const Dtype* input,
    const Dtype* weights, Dtype* output, bool skip_im2col) {
  const Dtype* col_buff = input;
  if (!is_1x1_) {
    if (!skip_im2col) {
      conv_im2col_cpu(input, col_buffer_.mutable_cpu_data());
    }
    col_buff = col_buffer_.cpu_data();
  }
  for (int g = 0; g < group_; ++g) {
    caffe_cpu_gemm<Dtype>(CblasNoTrans, CblasNoTrans, conv_out_channels_ /
        group_, conv_out_spatial_dim_, kernel_dim_,
        (Dtype)1., weights + weight_offset_ * g, col_buff + col_offset_ * g,
        (Dtype)0., output + output_offset_ * g);
  }
}

这里使用conv_im2col_cpu(input, col_buffer_.mutable_cpu_data());
对数据进行重排，简单说来，就是把本来的行列矩阵的图像数据转换成以适合做卷积运算的数据形式。
conv_im2col_cpu(input, col_buffer_.mutable_cpu_data());封装在vision_layer.hpp里面，是个内联函数。

inline void conv_im2col_cpu(const Dtype* data, Dtype* col_buff) {
    if (!force_nd_im2col_ && num_spatial_axes_ == 2) {
      im2col_cpu(data, conv_in_channels_,
          conv_input_shape_.cpu_data()[1], conv_input_shape_.cpu_data()[2],
          kernel_shape_.cpu_data()[0], kernel_shape_.cpu_data()[1],
          pad_.cpu_data()[0], pad_.cpu_data()[1],
          stride_.cpu_data()[0], stride_.cpu_data()[1], col_buff);
    } else {
      im2col_nd_cpu(data, num_spatial_axes_, conv_input_shape_.cpu_data(),
          col_buffer_shape_.data(), kernel_shape_.cpu_data(),
          pad_.cpu_data(), stride_.cpu_data(), col_buff);
    }
  }

可以看出，为了方便函数调用，这里省略了一些参数，将im2col_cpu()的参数封装了一下。具体实现在./src/caffe/util/im2col.cpp里面

template <typename Dtype>
void im2col_cpu(const Dtype* data_im, const int channels,
    const int height, const int width, const int kernel_h, const int kernel_w,
    const int pad_h, const int pad_w,
    const int stride_h, const int stride_w,
    Dtype* data_col) {
  int height_col = (height + 2 * pad_h - kernel_h) / stride_h + 1;
  int width_col = (width + 2 * pad_w - kernel_w) / stride_w + 1;
  int channels_col = channels * kernel_h * kernel_w;
  for (int c = 0; c < channels_col; ++c) {
    int w_offset = c % kernel_w;
    int h_offset = (c / kernel_w) % kernel_h;
    int c_im = c / kernel_h / kernel_w;
    for (int h = 0; h < height_col; ++h) {
      for (int w = 0; w < width_col; ++w) {
        int h_pad = h * stride_h - pad_h + h_offset;
        int w_pad = w * stride_w - pad_w + w_offset;
        if (h_pad >= 0 && h_pad < height && w_pad >= 0 && w_pad < width)
          data_col[(c * height_col + h) * width_col + w] =
            data_im[(c_im * height + h_pad) * width + w_pad];
        else
          data_col[(c * height_col + h) * width_col + w] = 0;
      }
    }
  }
}

const int channels, const int height, const int width是袁术数据const Dtype* data_im的信息，通道数还有高和宽。剩下的参数是kernel的大小，数据填充的大小，还有每次平移的长度。

int height_col = (height + 2 * pad_h - kernel_h) / stride_h + 1;
  int width_col = (width + 2 * pad_w - kernel_w) / stride_w + 1;

是算出每个channel的patch卷积操作后得到的数据大小。
int channels_col = channels * kernel_h * kernel_w;这是channel数乘以kernel的大小。首先想到的这是卷积层所有参数的个数。然后姑且认为这是得到数据的行数，每一行填充的是什么呢？
首先

int w_offset = c % kernel_w;
int h_offset = (c / kernel_w) % kernel_h;

这是每个行数索引相对于kernel的行和列位置
然后对所有的目标的数据（卷积后的）位置，得到

int h_pad = h * stride_h - pad_h + h_offset;
int w_pad = w * stride_w - pad_w + w_offset;

很显然，h_pad是卷积后h所对应的卷积前的c参数链接位置。
然后把数据通过最后两行复制过去。通过这很明显看到im2col_cpu()干了什么事情。
这个操作是为了方便计算卷积操作的，假设原数据是channels*height*width。卷积之后成了（channels*kernel_h*kernel_w）x（height_col*wight_col）。我们知道卷积操作有channels*kernel_h*kernel_w参数，代表kernel的weight，一共有channle个kernel。每个kernel的一个位置和（height_col*wight_col）有边相连。这样的话计算卷积的话，只需要算个矩阵乘法。
滴啊用cblas的矩阵乘法运算

void caffe_cpu_gemm<double>(const CBLAS_TRANSPOSE TransA,
    const CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
    const double alpha, const double* A, const double* B, const double beta,
    double* C) {
  int lda = (TransA == CblasNoTrans) ? K : M;
  int ldb = (TransB == CblasNoTrans) ? N : K;
  cblas_dgemm(CblasRowMajor, TransA, TransB, M, N, K, alpha, A, lda, B,
      ldb, beta, C, N);
}

函数定义为：

void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,

const enum CBLAS_TRANSPOSE TransB, const int M, const int N,

const int K, const float alpha, const float *A,

const int lda, const float *B, const int ldb,

const float beta, float *C, const int ldc)

得到的结果是:

C = alpha*op( A )*op( B ) + beta*C

const enum CBLAS_ORDER Order，这是指的数据的存储形式，在CBLAS的函数中无论一维还是二维数据都是用一维数组存储，这就要涉及是行主序还是列主序，在C语言中数组是用 行主序，fortran中是列主序。我还是习惯于是用行主序，所以这个参数是用CblasRowMajor，如果是列主序的话就是 CblasColMajor。

const int M，矩阵A的行，矩阵C的行

const int N，矩阵B的列，矩阵C的列

const int K，矩阵A的列，矩阵B的行

const float alpha， const float beta，计算公式中的两个参数值，如果只是计算C=A*B，则alpha=1,beta=0

const float *A， const float *B， const float *C，矩阵ABC的数据

const int lda， const int ldb， const int ldc，在BLAS的文档里，这三个参数分别为ABC的行数，但是实际使用发现，在CBLAS里应该是列数。

The following program computes the product of two matrices using the Level-3 BLAS function SGEMM,

[ 0.11 0.12 0.13 ] [ 1011 1012 ] [ 367.76 368.12 ]
[ 0.21 0.22 0.23 ] [ 1021 1022 ] = [ 674.06 674.72 ]
[ 1031 1032 ]
The matrices are stored in row major order but could be stored in column major order if the first argument of the call to cblas_sgemm was changed to CblasColMajor.

include

include