卷积神经网络加速计算

卷积加速方法

上一次实验中实现了一个简单的卷积层，但实际应用中当数据量大了之后不可能这么直接 for 循环执行代码，必须提高卷积的运行效率。

提高卷积运行效率有以下几种方法:
并行加速
选择更底层的语言，例如 C/C++/Fortran
算法加速
GPU 加速

关于第二点，使用 C/C++/Fortran 等更底层的语言不必多说，因为可以直接操作内存，相对于 Python 速度要快得多。但是这些语言开发效率要比 Python 等高级语言慢。不在本次课程内容范围，所以不会涉及。
下面主要讨论其他三点，首先介绍算法加速，其次是并行加速，最后则是 GPU 加速。

算法加速

上一次实验实现卷积层时，选择的方法简单粗暴，使用了很多 NumPy 方法来避免重复嵌套 for 循环，但是其本质仍然是一层一层 for 循环嵌套而来。

for w in 1..W
  for h in 1..H
    for x in 1..K
      for y in 1..K
        for m in 1..M
          for d in 1..D
            output(w, h, m) += input(w+x, h+y, d) * filter(m, x, y, d)
          end
        end
      end
    end
  end
end

这样做的效率当然很低，基本上没有深度学习框架采用这种方法。下面将会介绍 Caffe 中实现卷积的方法，并用 Python 代码实现。
Caffe 加速卷积的方法很简单，主要是利用 im2col 和 col2im 两种思想，以作者贾扬清在 在 Caffe 中如何计算卷积？ 的回答为例进行讲解。
im2col 方法

image.png

image.png

image.png

image.png

image.png

def im2col(data, kernel_size, stride):
    batch_size, input_channel, height, width = data.shape
    # 计算 Feature map 的大小
    feature_height = int((height-kernel_size)/stride)+1
    feature_width = int((width-kernel_size)/stride)+1

    # 初始化展平矩阵的大小, B*(H*W)*(C*K*K)
    col_data = np.zeros((batch_size, feature_height*feature_width,
                         kernel_size*kernel_size*input_channel), dtype=np.float32)

    # 卷积的滑窗
    for n in range(batch_size):
        for i in range(feature_height):
            for j in range(feature_width):
                # 将该窗口的数据展平保存
                col_data[n, i*feature_width+j, :] = np.ravel(
                    data[n, :, i*stride: i*stride+kernel_size, j*stride: j*stride+kernel_size])

    # 返回展平后的结果，和 Feature map 的高、宽
    return col_data, feature_height, feature_width

使用一个3×3 的矩阵，大小2×2 步长为 1 的卷积核进行测试。

import numpy as np

kernel = np.array([[1, 2], [2, 1]], dtype=np.float32).reshape((1, 1, 2, 2))
input = np.arange(1, 10).reshape((1, 1, 3, 3))

col_data, feature_height, feature_width = im2col(input, 2, 1)
col_data

结果是一致的，测试展开后的数据乘以展开后的卷积核的转置，是否是该结果。

np.dot(col_data[0], kernel.reshape(1, -1).T).reshape(2, 2)

答案是正确的。
col2im 方法
使用 im2col 之后乘积再 np.reshape 即可获得输出的 Feature map，难点在于如何将 Feature map 的梯度反向传播至该层，计算该层输入、权重、偏置的梯度。
以单个卷积为例，使用 im2col 之后再乘以卷积核，仍然是执行以下操作:

image.png

# 计算输入的梯度，将展开的梯度还原
def col2im(col_data, top_grad, weight, shape):
    # 参数
    batch_size, input_channel, width, height, feature_height, feature_width, kernel_size, stride= shape
    # 初始化原始梯度，和输入数据一样
    grad = np.zeros((batch_size, input_channel, width, height), dtype=np.float32)
    # 对每个样本的计算梯度
    grad_one_ = np.matmul(top_grad, weight)
    for n in prange(batch_size):
        for i in prange(feature_height):
            for j in prange(feature_width):
                # 每个样本的梯度累加
                # 被展开的梯度，还原成原始数据的梯度
                grad[n, :, i*stride:i*stride+kernel_size, j*stride:j*stride+kernel_size] += np.reshape(grad_one_[n, i*feature_width+j, :], (input_channel, kernel_size, kernel_size))
    return grad

并行加速

采用以上算法对卷积的前向、后向传播进行加速之后，可以提高很多运行速度。但是虽然在关键步骤上加速了，上面介绍的算法还是有很多地方进行了多层 for 循环。这时候就需要优化 for 循环，可以采用并行的方法，使用多线程并行运算 for 循环，就不需要一个接一个地运行。
C++ 有类似于 OpenMP、OpenCL 等库进行 CPU 并行加速，Python 中也有类似的库，下面将以一个为例。
Numba 是 Python 的一个加速库，通过在编译时替换 NumPy 的函数成 Numba 实现的函数，并使用并行运算进行加速，支持多线程并行（parallel），无 Python （nopython），快速数学公式（fastmath）等。只需要在函数前添加 @jit 即可完成加速，下面以一个例子为例：

from numba import jit
import random

def monte_carlo_pi(nsamples):
    acc = 0
    for i in range(nsamples):
        x = random.random()
        y = random.random()
        if (x**2 + y**2) < 1.0:
            acc += 1
    return 4.0 * acc / nsamples

@jit(nopython=True, parallel=True)
def monte_carlo_pi_numba(nsamples):
    acc = 0
    for i in range(nsamples):
        x = random.random()
        y = random.random()
        if (x**2 + y**2) < 1.0:
            acc += 1
    return 4.0 * acc / nsamples

进行性能测试查看是否提高：

%timeit monte_carlo_pi(100)
%timeit monte_carlo_pi_numba(100)

代码运行速度提高了数十倍。

GPU 加速

相对于 CPU 来说，GPU 的计算能力要比 CPU 强得多。深度学习普遍计算量非常大，使用 CPU 训练基本不可能，只能使用一个、甚至多个 GPU 进行并行训练。
现有的深度学习框架都支持 GPU 加速，大部分都是在底层使用 C/C++ 编写 CUDA 核函数调用 GPU，也就是 NVIDIA 提供的一个工具包。极少数开源项目才会使用 AMD 显卡，通过 OpenCL 进行调用。OpenCL 虽然开源、支持的平台很多，但在开发上远没有 CUDA 方便，另外 NVIDIA 在深度学习基本上是垄断的，AMD 显卡一般被用于挖矿。更多原因可以参考 为什么在部分机器学习中训练模型时使用 GPU 的效果比 CPU 更好？。
下面将基于 PyTorch 的 GPU 版本，对比大数据量时 CPU 和 GPU 的运行速度，由于线上环境无 GPU，所以只提供图片显示:

image.png

重新实现卷积层

下面将会使用这三种加速方法重新实现卷积层:

import torch

# 是否有 GPU，需要配置 PyTorch GPU 环境
has_gpu = torch.cuda.is_available()
has_gpu

from numba import autojit, prange

@autojit
def im2col(data, kernel_size, stride):
    batch_size, input_channel, height, width = data.shape
    # 根据公式计算 feature map 的大小
    feature_height = int((height-kernel_size)/stride)+1
    feature_width = int((width-kernel_size)/stride)+1 
    
    # 初始化展平矩阵的大小, B*(H*W)*(C*K*K)
    col_data = np.zeros((batch_size, feature_height*feature_width, kernel_size*kernel_size*input_channel), dtype=np.float32)
    
    # 卷积的滑窗
    for n in prange(batch_size):
        for i in prange(feature_height):
            for j in prange(feature_width):
                # 将该窗口的数据展平保存
                col_data[n, i*feature_width+j, :] = np.ravel(data[n, :, i*stride: i*stride+kernel_size, j*stride: j*stride+kernel_size])
            
    # 返回展平后的结果，和 feature map 的高、宽
    return col_data, feature_height, feature_width

def matmul(input1, input2):
    assert input1.shape[0]==input2.shape[0], '必须相等'
    if has_gpu:
        grad = torch.sum(torch.einsum('ijk,ikl->ijl', (torch.from_numpy(input1).cuda(), torch.from_numpy(input2).cuda())), dim=0).cpu().numpy()
    else:
        grad = np.sum(np.einsum('ijk,ikl->ijl', input1, input2), axis=0)
    return grad

# 计算输入的梯度，将展开的梯度还原
@autojit
def col2im(col_data, top_grad, weight, shape):
    # 参数
    batch_size, input_channel, width, height, feature_height, feature_width, kernel_size, stride= shape
    # 初始化原始梯度，和输入数据一样
    grad = np.zeros((batch_size, input_channel, width, height), dtype=np.float32)
    # 对每个样本的计算梯度
    if has_gpu:
        grad_one_ = torch.matmul(torch.from_numpy(top_grad).cuda(), torch.from_numpy(weight).cuda()).cpu().numpy()
    else:
        grad_one_ = np.matmul(top_grad, weight)
    for n in prange(batch_size):
        for i in prange(feature_height):
            for j in prange(feature_width):
                # 每个样本的梯度累加
                # 被展开的梯度，还原成原始数据的梯度
                grad[n, :, i*stride:i*stride+kernel_size, j*stride:j*stride+kernel_size] += np.reshape(grad_one_[n, i*feature_width+j, :], (input_channel, kernel_size, kernel_size))
    return grad

class conv2d(object):
    '''
    output_channel: 该层卷积核的数量
    input_channel: 输入数据的通道数，例如灰度图为 1，彩色图为 3，又或者上一层的卷积数量为 12，那下一层卷积的输入通道为 12
    kernel_size: 卷积核大小
    stride: 卷积移动的步长
    padding: 补齐
    '''
    def __init__(self, input_channel, output_channel, kernel_size, stride=1, padding=0):
        # 初始化网络权重，卷积核数量*输入通道数*卷积核大小*卷积核大小，一个卷积核应为输入通道数*卷积核大小*卷积核大小，每一层有多个卷积核
        # self.weight = (np.random.randn(output_channel, input_channel, kernel_size, kernel_size)*0.01).astype(np.float32)
        self.weight = copy.deepcopy(torch.nn.Conv2d(input_channel, output_channel, kernel_size, stride).weight.data.numpy())
        # 偏置项，每一个卷积核都对应一个偏置项，y=xW+b
        self.bias = np.zeros((output_channel), dtype=np.float32)
        # 保存各项参数
        self.output_channel = output_channel
        self.input_channel = input_channel
        self.kernel_size = kernel_size
        self.stride = stride
        self.padding = padding
        
    def forward(self, input):
        # 输入数据各个维度代表的含义: batch_size, input_channel, height, width
        batch_size, input_channel, height, width = input.shape
        # 防止一些错误的调用
        assert len(input.shape)==4, '输入必须是四维的，batch_size*channel*width*height'
        # 上一层的卷积核数量必须等于这一层的输入通道数
        assert input_channel==self.input_channel, '网络输入通道数必须与网络定义的一致'
        
        # 对输入数据进行补齐
        input = np.pad(input, [(0, ), (0, ), (self.padding, ), (self.padding, )], 'constant', constant_values=0)
        
        # 保存输入数据，后续需要计算梯度
        self.data = input
        
        # 展开数据
        self.col_data, self.feature_height, self.feature_width = im2col(input, self.kernel_size, self.stride)
        # 计算输出，是否有 GPU，如果有则使用 GPU 计算
        if has_gpu:
            feature_maps = torch.matmul(torch.from_numpy(self.col_data).cuda(), torch.from_numpy(self.weight.reshape(self.output_channel, -1).T).cuda()) + torch.from_numpy(self.bias).cuda().unsqueeze(0)
            feature_maps = torch.reshape(feature_maps, (batch_size, self.feature_height, self.feature_width, self.output_channel)).permute(0, 3, 1, 2).cpu().numpy()
        else:
            # 计算输出之后同时 reshape, transpose 还原成本来应该的输出大小
            # feature_maps = np.dot(self.col_data, self.weight.reshape(self.output_channel, -1).T) + self.bias[np.newaxis, :]
            feature_maps = np.matmul(self.col_data, self.weight.reshape(self.output_channel, -1).T) + self.bias[np.newaxis, :]
            feature_maps = np.reshape(feature_maps, (batch_size, self.feature_height, self.feature_width, self.output_channel)).transpose(0, 3, 1, 2)
        return feature_maps
    
    '''
    top_grad: 上一层的梯度，维度为 batch_size, output_channel, feature_height, feature_width，与该层输出一致
    '''
    def backward(self, top_grad, lr):
        batch_size, output_channel, feature_height, feature_width = top_grad.shape
        # 将梯度展开
        top_grad = top_grad.transpose((0, 2, 3, 1)).reshape(batch_size, feature_height*feature_width, output_channel)
        # 计算权重、偏置的梯度
        self.grad_w = matmul(self.col_data.transpose(0, 2, 1), top_grad).T.reshape(output_channel, self.input_channel, self.kernel_size, self.kernel_size)
        self.grad_b = np.sum(top_grad, axis=(0, 1))
        # 各项参数
        shape = list(self.data.shape) + [self.feature_height, self.feature_width, self.kernel_size, self.stride]
        # 计算输入的梯度
        self.grad = col2im(self.col_data, top_grad, self.weight.reshape(self.output_channel, -1), shape)    
        # 更新参数
        self.weight -= lr*self.grad_w
        self.bias -= lr*self.grad_b

至此卷积加速完成，但即使是这样，相对于 PyTorch 实现的卷积层依然有数倍的差距，只能归咎于 Python 的性能问题和 PyTorch 的优化方法更有效了。
优化完成之后，使用一个示例进行介绍。

CIFAR-10 分类

在前面的实验中，选择的 MNIST 手写体数据，这个数据集训练集和测试集一共 <math xmlns="http://www.w3.org/1998/Math/MathML"><semantics><annotation encoding="application/x-tex">70000\times 1\times 28\times 28</annotation></semantics></math>70000×1×28×28 大小，在深度学习中是一个比较常用的数据集。但 MNIST 数据集是一个灰度图的数据集，所以还是比较简单的。
接下来介绍另一个比较经典的彩色图片分类问题。
CIFAR-10 数据集一共 10 类，60000 张32×32 的彩色图片，每个类别 6000 张，其中训练集 50000 张，每个类别 5000 张，其他的都是测试图片并且均匀分布。所以 CIFAR-10 数据集一共有60000×3×32×32，光是数据量就是 MNIST 的三倍左右，训练难度也更大。

image.png

训练之前下载 CIFAR10 数据集:

!wget -nc "http://labfile.oss.aliyuncs.com/courses/1213/cifar-10-batches-py.zip"
!unzip -o "cifar-10-batches-py.zip"

接下来，我们直接使用 PyTorch 来定义网络结构，并实现训练过程。

from torchvision import datasets, transforms
import torch.nn.functional as F
from tqdm.notebook import tqdm

class Net(torch.nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.conv = torch.nn.Sequential(
            torch.nn.Conv2d(3, 64, 5, 1),
            torch.nn.ReLU(inplace=True),
            torch.nn.MaxPool2d(2, 2),
            torch.nn.Conv2d(64, 128, 3, 1),
            torch.nn.ReLU(inplace=True),
            # BatchNorm 正则化方法
            torch.nn.BatchNorm2d(128),
            torch.nn.MaxPool2d(2, 2),
            torch.nn.Conv2d(128, 256, 2, 1),
            torch.nn.ReLU(inplace=True),
            torch.nn.Conv2d(256, 256, 2, 1),
            torch.nn.ReLU(inplace=True),
            # BatchNorm 正则化方法
            torch.nn.BatchNorm2d(256),
            torch.nn.MaxPool2d(2, 2),
        )
    
        self.classifier = torch.nn.Sequential(
            torch.nn.Linear(1024, 1024),
            torch.nn.ReLU(inplace=True),
            torch.nn.Linear(1024, 10),
        )

    def forward(self, x):
        x = self.conv(x)
        x = x.view(x.size(0), -1)
        x = self.classifier(x)
        return x

设置各项参数并加载数据集:

# 迭代次数、学习率等
batch_size = 120
base_lr = 0.1
EPOCHS = 20
step_size = 8
download = True
best_acc = -float('inf')
    
# 加载数据
train_loader = torch.utils.data.DataLoader(
    datasets.CIFAR10('.', train=True, download=download,
                   transform=transforms.Compose([
                       transforms.ToTensor(),
                   ])),
    batch_size=batch_size, shuffle=True, num_workers=1)
test_loader = torch.utils.data.DataLoader(
    datasets.CIFAR10('.', train=False, transform=transforms.Compose([
                       transforms.ToTensor(),
                   ])),
    batch_size=batch_size, shuffle=False, num_workers=1)

# 定义网络结构和优化器
model = Net()
if torch.cuda.is_available():
    model = model.cuda()

# weight_decay 表明使用权重衰减的系数，不宜过大，一般取小数点四位再慢慢调整
optimizer = torch.optim.SGD(model.parameters(), lr=base_lr, momentum=0.9, weight_decay=0.001)
exp_lr_scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=step_size, gamma=0.1)

开始训练:

# 开始迭代
for epoch in range(1, EPOCHS + 1):
    model.train()
    train_loss = 0
    correct = 0
    # 训练
    for data, target in tqdm(train_loader):
        if torch.cuda.is_available():
            data = data.cuda()
            target = target.cuda()
        optimizer.zero_grad()
        output = model(data)
        loss = F.cross_entropy(output, target)
        train_loss += loss.item()
        pred = output.max(1, keepdim=True)[1]
        correct += pred.eq(target.view_as(pred)).sum().item()
        loss.backward()
        optimizer.step()
    exp_lr_scheduler.step()
    train_loss /= len(train_loader.dataset)
    print('Epoch {}/{}:'.format(epoch, EPOCHS))
    print('Train set: Average loss: {:.4f}, Accuracy: {}/{} ({:.2f}%)'.format(
        train_loss, correct, len(train_loader.dataset),
        100. * correct / len(train_loader.dataset)))

    model.eval()
    test_loss = 0
    correct = 0
    # 测试
    with torch.no_grad():
        for data, target in test_loader:
            if torch.cuda.is_available():
                data = data.cuda()
                target = target.cuda()
            output = model(data)
            test_loss += F.cross_entropy(output, target, reduction='sum').item()
            pred = output.max(1, keepdim=True)[1]
            correct += pred.eq(target.view_as(pred)).sum().item()

    test_loss /= len(test_loader.dataset)
    acc = correct/len(test_loader.dataset)
    print('Test set: Average loss: {:.4f}, Accuracy: {}/{} ({:.2f}%)'.format(
        test_loss, correct, len(test_loader.dataset),
        100. * acc))
    if best_acc<acc: best_acc=acc
    print()

print('best accuracy: ', best_acc)

运行结束之后，迭代二十次获得的最高准确度大约 80%，对于一个如此简单的网络来说，这个表现已经非常好了。

实现一个简单卷积神经网络

实现一个简单卷积神经网络

前面的课程中，我们是这样实现一个简单的 LeNet-5 网络结构的：

class LeNet(object):
    def __init__(self):
        self.conv1 = conv2d(1, 20, 5, 1)
        self.relu1 = ReLU()
        self.pool1 = MaxPool2d(2, 2)
        self.conv2 = conv2d(20, 50, 5, 1)
        self.relu2 = ReLU()
        self.pool2 = MaxPool2d(2, 2)
        self.fc1 = Linear(800, 500)
        self.relu3 = ReLU()
        self.fc2 = Linear(500, 10)
        
    def forward(self, input):
        input = self.relu1.forward(self.conv1.forward(input))
        input = self.pool1.forward(input)
        input = self.relu2.forward(self.conv2.forward(input))
        input = self.pool2.forward(input)
        # 展开
        self.flatten_shape = input.shape
        input = np.reshape(input, (input.shape[0], -1))
        input = self.relu3.forward(self.fc1.forward(input))
        output = self.fc2.forward(input)
        return output
    
    def backward(self, top_grad, lr):
        self.fc2.backward(top_grad, lr)
        self.relu3.backward(self.fc2.grad)
        self.fc1.backward(self.relu3.grad, lr)
        unflattened_grad = np.reshape(self.fc1.grad, self.flatten_shape)
        self.pool2.backward(unflattened_grad)
        self.relu2.backward(self.pool2.grad)
        self.conv2.backward(self.relu2.grad, lr)
        self.pool1.backward(self.conv2.grad)
        self.relu1.backward(self.pool1.grad)
        self.conv1.backward(self.relu1.grad, lr)

一个网络结构，分别有 forward 和 backward 函数，对呀前向和后向传播过程。在某些框架中 backward 可能省略。其中一个网络结构又由许多基本层构成，每层也同样对应了 forward 和 backward 函数。
接下来请你模仿 LeNet-5 的实现，按照以下网络结构封装一个 SimpleNet 类:

image.png

各项参数如下:
输入数据为 $B110*10，，B$ 为批量大小，即输入通道数为 1。
第一层卷积层个数为 20，卷积核大小为 5，步长为 1。
ReLU 激活层。
步长为 2，核大小为 2 的最大池化层。
Flatten 之后接入一层线性层、ReLU 激活的模块，线性层输入为 180，输出为 100。
最后一层进行分类，无 ReLU，输出为 10.
在开始前，下载实验中封装好的网络层代码：

!wget -nc "http://labfile.oss.aliyuncs.com/courses/1213/nn.py"

请根据以上描述，完成设计网络结构。