[深入推导]CS231N assignment 2#4 _ 卷积神经网络 学习笔记 & 解析

卷积神经网络

基本算法实现

卷积神经网络应该算是图像处理中绝对的主流了, 关于算法得基本思想我在之前也学的比较懂了, 这点如果不了解网上有很多教程.

不过我并没有用代码亲自实现它. 我们首先确定怎么编写. 前面搞全连接网络总是会想着怎么去简化运算, 现在我们接触了新的网络, 要实现基础版本反而又不大适应了. 但是既然不要求计算准确性,那不就是一个一个像素填充么. (后面也没有要求我们用没有for循环版本,估计应该是太难了吧)

我们已经明白, 卷积本质上就是矩阵得内积, 对图像三通道相加之后加上偏置项. 所以算法上就是矩阵切片和相乘. 下面我们确定输出得规模和细节. 在作业注释中,已经说明:

- x: Input data of shape (N, C, H, W)
    - w: Filter weights of shape (F, C, HH, WW)
    - b: Biases, of shape (F,)
    - conv_param: A dictionary with the following keys:
      - 'stride': The number of pixels between adjacent receptive fields in the
        horizontal and vertical directions.
      - 'pad': The number of pixels that will be used to zero-pad the input.
        Returns a tuple of:
    - out: Output data, of shape (N, F, H', W') where H' and W' are given by
      H' = 1 + (H + 2 * pad - HH) / stride  # 更准确应当是整除
      W' = 1 + (W + 2 * pad - WW) / stride
    - cache: (x, w, b, conv_param)

也就是说: 输入图像为: 图片数目,通道数, 高,宽; 卷积核为每个通道得卷积核数目, 通道数目, 卷积核得高和宽. stride和padding由参数指定. 随后, 我们从padding后得左上角开始移动, 直到移动stride长度后会超出边界为止.由此就可以得到规模得计算公式了.

下面就是算法本体. 

    N, C, H, W = x.shape # 获取shape信息
    F, C, HH, WW = w.shape
    stride, pad = conv_param['stride'], conv_param['pad'] # 读取参数
    # 整除
    H_out = 1 + (H + 2 * pad - HH) // stride
    W_out = 1 + (W + 2 * pad - WW) // stride
    # 预先分配
    out = np.zeros((N, F, H_out, W_out)) # 确定输出规模

    # 只对图像本体扩充,即第三和第四维度,利用0填充
    x_pad = np.pad(x, ((0, 0), (0, 0), (pad, pad), (pad, pad)), mode='constant', constant_values=0)

    for i in range(N):  # 对于第i幅图片
        for f in range(F):  # 对于第f个卷积核(per 通道)
            for j in range(H_out): # 高
                for k in range(W_out): # 宽
                    out[i, f, j, k] = np.sum(
                        x_pad[i, :, j * stride: HH + j * stride, k * stride: WW + k * stride] * w[f],axis=None) + b[f]
                    # 切片参数:第i张图片,全部通道,从j*stride开始到ans+HH, 与整个w逐个元素相乘并直接求和

卷积的图像意义

原始作业中用一种有趣的方式揭示了卷积运算的内涵. 在图像处理的学习中, 我们曾接触过sobel算子, 其可以用于边缘检测. 它就是通过特定的卷积核实现的. 

这个作业通过RGB通道权重1:6:3混合, 并设置w为Gy的sobel算子,. 强化了y向边缘. 结果如下:

梯度下降

别忘了我们是naive版本的, 所以我们可以逐个像素逐个像素求解的. 为此我们列出这个像素和w,x的那些位置有关, 结合表达式求导. 请注意,下面的分析我们都是先对X做padding处理, 在进行计算的.

所以代码也就出来了:

x, w, b, conv_param = cache # 取出数据
    N, C, H, W = x.shape
    F, C, HH, WW = w.shape
    stride, pad = conv_param['stride'], conv_param['pad']
    N, F, H_out, W_out = dout.shape

    # 先对X填充0
    x_pad = np.pad(x, ((0, 0), (0, 0), (pad, pad), (pad, pad)), mode='constant', constant_values=0)
    # 预分配
    dx = np.zeros_like(x)
    dx_pad = np.zeros_like(x_pad)
    dw = np.zeros_like(w)
    db = np.zeros_like(b)

    for i in range(N):  # 对输出的元素进行逐个遍历
        for f in range(F):  
            for j in range(H_out):
                for k in range(W_out):
                    window = x_pad[i, :, j * stride: HH + j * stride, k * stride: WW + k * stride] # 取出涉及的X区域
                    db[f] += dout[i, f, j, k] # db的结果就是dout
                    dw[f] += window * dout[i, f, j, k] # dw = x_i * dout
                    # dx_i = dw * dout
                    dx_pad[i, :, j * stride:HH + j * stride, k * stride: WW + k * stride] += w[f] * dout[i, f, j, k] 


    dx = dx_pad[:, :, pad: -pad, pad: -pad] # 截取

Max-Pooling 池化层

池化层本质上就是对图像做降采样. 我们或许之前常常听到的是图像2*2降采样分辨率, 但实际上, 最通用的池化层是类似于前面的滑动窗口的方式选择最大值, 有窗口大小和步长. 因此, 我们也就很容易类似于前面的卷积得到公式: (作业有给出)

H' = 1 + (H - pool_height) / stride  # 准确来说还是整除
      W' = 1 + (W - pool_width) / stride

 所以结果也就得出(注意要善用切片,不要写四层循环了):

N, C, H, W = x.shape
    pool_height, pool_width, stride = pool_param['pool_height'], pool_param['pool_width'], pool_param['stride'] # 读取参数
    H_out = 1 + (H - pool_height) // stride
    W_out = 1 + (W - pool_width) // stride

    out = np.zeros((N, C, H_out, W_out))
    for j in range(H_out): # 对每个像素
        for k in range(W_out):
            out[:, :, j, k] = np.max(x[:, :, j * stride:pool_height + j * stride, k * stride: pool_width + k * stride],
                                     axis=(2, 3)) # 取axis2,3最大值

 反向推导, 因为池化层选取得是最大得那个元素, 所以梯度就和那个元素有关. 我们要做的, 就是得到池化层对应的那个最大是何元素, 然后对它的dx[...] += dout[...].

    x, pool_param = cache
    N, C, H, W = x.shape
    pool_height, pool_width, stride = pool_param['pool_height'], pool_param['pool_width'], pool_param['stride']
    H_out = 1 + (H - pool_height) // stride
    W_out = 1 + (W - pool_width) // stride
    dx = np.zeros_like(x)

    for i in range(H_out): # 逐个像素
        for j in range(W_out):
            x_masked = x[:, :, i * stride:i * stride + pool_width, j * stride:j * stride + pool_height] # 获取原始切片区域
            max_x_masked = np.max(x_masked, axis=(2, 3))[:, :, None, None] #得到的最大的结果,None相当于将维度扩大了以便于加回去
            # 反馈到原本的位置去. shape从(a,b) => (a,b,1,1) 可以利用广播机制
            outmask = (x_masked == max_x_masked) 
            dx[:, :, i * stride:i * stride + pool_height, j * stride:j * stride + pool_width] += outmask * dout[:,:,i,j][:,:,None,None] 
    pass

快速求解算法初探

虽然题目没有要求做到, 但还是看了一眼.其有一个大致思想, 即: 我们将卷积核变成一个向量, 输入数据变成矩阵(注意为了方便矩阵乘法产生了大量冗余), 随后矩阵乘法(还是那句话,内积可以通过恰当的方式变成矩阵乘法!)之后输出一个向量, reshape之后得到了标准输出.   我还是按照传统艺能给你解释一下, 你应该就明白了(这个解释被极度简化了)

那问题来了, W很好解释, X怎么办呢? 怎么得到这个大矩阵? 这就是np.lib.stride_tricks.as_strided这个函数所做的事了. 关于这个函数的技巧请参考numpy as_strided 入门教程（啰嗦版） - 知乎 (zhihu.com).

def conv_forward_strides(x, w, b, conv_param):
    N, C, H, W = x.shape
    F, _, HH, WW = w.shape
    stride, pad = conv_param['stride'], conv_param['pad']
    p = pad
    x_padded = np.pad(x, ((0, 0), (0, 0), (p, p), (p, p)), mode='constant') # 补0
    H += 2 * pad
    W += 2 * pad
    out_h = (H - HH) // stride + 1
    out_w = (W - WW) // stride + 1   # 到这里为止都还不难
    # ---- 反人类预警 ----
    shape = (C, HH, WW, N, out_h, out_w) # 首先我们初步确认规模
    strides = (H * W, W, 1, C * H * W, stride * W, stride) 
    strides = x.itemsize * np.array(strides)
    x_stride = np.lib.stride_tricks.as_strided(x_padded, # 这里的reshape有点特殊,它可以自由切片reshape同时冗余分配,具体请搜搜这个函数
                  shape=shape, strides=strides)
    x_cols = np.ascontiguousarray(x_stride) # 数据连续存放
    x_cols.shape = (C * HH * WW, N * out_h * out_w)  # reshape => X大矩阵的目标维度
    res = w.reshape(F, -1).dot(x_cols) + b.reshape(-1, 1)
    res.shape = (F, N, out_h, out_w)
    out = res.transpose(1, 0, 2, 3) # 结果reshape
    out = np.ascontiguousarray(out)
    cache = (x, w, b, conv_param, x_cols)
    return out, cache

这个确实很棒, 避免了对但是实际采用的却不是这种方法, 而是直接写循环, 啊这...

不过这个代码却不完全是python, 而是cython. 利用cython这个工具, 可以在python内调用C底层. 这个代码虽然是python风格的, 但是最终生成的却是c语言代码, 复杂度没有变化, 但是加速效果显著! 有关这个后续我可能会研究, 又要挖坑了...

sandwich层

这个也是题目已经实现好的, 就是直接凑出来,调用就行了. 我们实际上的操作就是  卷积 => 激活函数 => 池化.

a, conv_cache = conv_forward_fast(x, w, b, conv_param)
    s, relu_cache = relu_forward(a)
    out, pool_cache = max_pool_forward_fast(s, pool_param)
    cache = (conv_cache, relu_cache, pool_cache)
    return out, cache

三层卷积神经网络

我们想要复现的网络结构如下:

conv - relu - 2x2 max pool - affine - relu - affine - softmax

这个网络为一层卷积层和两层全连接层构成. 和前面构造两层神经网络一样, 写出代码:

# 初始化
        C,H,W = input_dim[0],input_dim[1],input_dim[2]
        self.params['W1'] = weight_scale * np.random.randn(num_filters,C,filter_size,filter_size) # 正方形卷积核
        self.params['b1'] = np.zeros(num_filters)
        # 全连接层,因为这里步长没有自己选择的余地, 根据默认参数,生成图像和输入图像维度一致,pool之后就是H/2*W/2
        self.params['W2'] = weight_scale * np.random.randn(num_filters*int(H/2)*int(W/2),hidden_dim) 
        self.params['b2'] = np.zeros(hidden_dim)
        self.params['W3'] = weight_scale * np.random.randn(hidden_dim,num_classes)
        self.params['b3'] = np.zeros(num_classes)        
        pass

        # 代码块1:前向推导
        out_conv, cache_conv = conv_relu_pool_forward(X, W1, b1, conv_param, pool_param)
        out_fc1, cache_fc1 = affine_relu_forward(out_conv, W2, b2)
        scores, cache_fc2 = affine_forward(out_fc1, W3, b3)
        #代码块2: loss & 反向推演
        loss, dout = softmax_loss(scores, y)
        loss += 0.5 * self.reg * (np.sum(W1 ** 2) + np.sum(W2 ** 2) + np.sum(W3 ** 2)) # 别忘记正则化
        dx3, dW3, db3 = affine_backward(dout, cache_fc2)
        dx2, dW2, db2 = affine_relu_backward(dx3, cache_fc1)
        dx1, dw1, db1 = conv_relu_pool_backward(dx2, cache_conv)
        grads['W3'] = dW3 + self.reg * W3
        grads['b3'] = db3
        grads['W2'] = dW2 + self.reg * W2
        grads['b2'] = db2
        grads['W1'] = dw1 + self.reg * W1
        grads['b1'] = db1

 随后为了验证正确性,仍然采取很小的数据集, 希望在训练集得到过拟合, 测试集的效果较差.

验证完毕之后就开始正式训练网络. 我们就此感受到了训练时间明显被拉长了, 这里仅仅训练了1个epoch, 在我的电脑上就花费了大约三分钟, 而之前的网络这个时间就可以跑至少50个epoch. 由此我们看出算力对于卷积网络(尤其是深度网络)的意义, 因为矩阵对GPU的天然优势, 后面我们就会转移阵地.但是效果确实不错, 仅仅1个epoch就达到了训练集50.4%,测试集49.9%的准确率了.

可视化的卷积核中, 看出其更加重视局部的特征匹配.

Spatial Batch Normalization

怎么将归一化用在卷积网络呢? 这里大概做法是: 对每个通道内部做正则化. 譬如我们的图片(或者上层输入)为N*C*H*W, 那我们对C个N*H*W内部去做正则化. 实际操作中, 我们希望直接用我们前面的成果而非再写一次, 就应当把图片化成矩阵的形式, 又因为batch norm是列向量的正则化, 所以生成矩阵应当是N*H*W,C的shape. 

N, C, H, W = x.shape
x_new = x.transpose(0, 2, 3, 1).reshape(N * H * W, C) # (N,C,H,W) => (N,H,W,C) => (N*H*W,C)
out, cache = batchnorm_forward(x_new, gamma, beta, bn_param) # 直接调用函数.这里cache用原来的目的: 反向也能用之前的backward函数
out = out.reshape(N, H, W, C).transpose(0, 3, 1, 2) # 变成原样

backward的代码如下:

N, C, H, W = dout.shape
dout_new = dout.transpose(0, 2, 3, 1).reshape(N * H * W, C) 
dx, dgamma, dbeta = batchnorm_backward_alt(dout_new, cache)
dx = dx.reshape(N, H, W, C).transpose(0, 3, 1, 2) # 直接对dout做reshape即可

Group Normalization

还记得我们曾在正则化的博客中提到, layer normalization需谨慎应用于CNN, 因为这么做会导致各个神经元之间的独特特征被抹平. 让我们设想一下, 原本的layer是行向量,也就是对C个通道内部正则化. 因为N*H*W之间没有关联, 那你想想这个图像不得乱套了. 这个方法就是从layer normalization所演化而来的. 这个部分在较早前版本的作业是没有的, 因为是2018年才发表的...

这个的做法是: 我们就缩小正则化的范围, 原本正则化的向量长度仅为C, 而现在被提高到了(C//G)*H*W, 向量数量为N*G, 也就是正则化的对象就是每张图片通道划分成数份去做,而不是武断地直接把所有图像全正则化了.

因为是比较近的论文, 所以直接就给出了gn的tensorflow实现方式... 

所以代码自然就清楚了

    N, C, H, W = x.shape
    # 将特征通道数分组，按照分组重新设置形状
    x_group = x.reshape((N, G, C // G, H, W))
    mean = np.mean(x_group, axis=(2, 3, 4), keepdims=True)
    var = np.var(x_group, axis=(2, 3, 4), keepdims=True)
    x_norm = (x_group - mean) / np.sqrt(var + eps)  # 归一化
    x_norm = x_norm.reshape((N, C, H, W))  # 还原维度
    out = x_norm * gamma + beta
    cache = (x, gamma, beta, G, eps, mean, var, x_norm)

可惜的是, 论文未给出backward地步骤. 我们回顾一下最开始实现的源码:

(x, mean, var, x_hat, eps, gamma, beta) = cache
    N = x.shape[0]
    dbeta = np.sum(dout, axis=0)
    dgamma = np.sum(dout * x_hat, axis=0)
    dx_hat = gamma * dout
    dvar = np.sum((x - mean) * dx_hat, axis=0) * (-0.5 / np.sqrt(var + eps) ** 3)
# 注意: sum的位置无关紧要,比如把-0.5 / np.sqrt(var + eps) ** 3)包含进sum不影响结果
    dmean = -1 / np.sqrt(var + eps) * np.sum(dx_hat, axis=0) + dvar * np.sum(-2 * (x - mean), axis=0) / N 
    # final gradient
    dx = 1 / np.sqrt(var + eps) * dx_hat + 1 / N * dmean + dvar * 2 * (x - mean) / N

当时, 从batch => layer只需要把axis变一下就行了. 现在, 我们从简单的gamma开始. 已经知道了规模:

Returns a tuple of:
- dx: Gradient with respect to inputs, of shape (N, C, H, W)
- dgamma: Gradient with respect to scale parameter, of shape (C,)
- dbeta: Gradient with respect to shift parameter, of shape (C,)

dbeta就是广播,反向对应求和,而dgamma是内积求和. 我们延续前面, 但考虑到规模问题, 所以axis=0应当变成axis=0,2,3.

dbeta = np.sum(dout, axis=(0, 2, 3), keepdims=True) 
    dgamma = np.sum(dout * x_norm, axis=(0, 2, 3), keepdims=True)
    dx_hat = dout * gamma

 后面怎么办? 我们还记得, mean和var是(N,G), 所以为了方便操作广播机制, 我们将得到的dhat转化成原本的格式, 随后一项一项求. 因为我们正则化是2,3,4,keepdim=True, 那现在就也是如此. 

    dx_hat_ = dx_hat.reshape((N, G, C // G, H, W))
    x_ = x.reshape((N, G, C // G, H, W)) # 将x和dxhat均reshape
    # 带入公式
    dvar = np.sum(dx_hat_ * (x_ - mean), axis=(2, 3, 4), keepdims=True) * -1.0 / 2 / np.sqrt(var + eps) ** (3)
    N_GROUP = C // G * H * W
    dmean = -1.0 / np.sqrt(var + eps) * np.sum(dx_hat_ , axis=(2, 3, 4), keepdims=True) \
            + dvar * -2.0 / N_GROUP * np.sum(x_ - mean, axis=(2, 3, 4), keepdims=True)

    dx_ = dx_hat_ * 1.0 / np.sqrt(var + eps)   \
                + dmean * 1.0 / N_GROUP  \
                + dvar * 2.0 / N_GROUP * (x_ - mean)

    dx = dx_.reshape((N, C, H, W))

-- おわり --