生成对抗网络(GAN)是一类非常有趣的神经网络。借助GAN，计算机能够生成逼真的图片。近年来有许多“AI绘画”的新闻，这些应用大多是通过GAN实现的。实际上，GAN不仅能做图像生成，还能辅助其他输入信息不足的视觉任务。比如SRGAN，就是把GAN应用在超分辨率(SR)任务上的代表之作。

在这篇文章中，我将主要面向深度学习的初学者，介绍SRGAN[1]这篇论文，同时分享以下知识：

• GAN的原理与训练过程
• 感知误差(Perceptual Loss)
• 基于的GAN的SR模型框架

讲完了知识后，我还会解读一下MMEditing的SRGAN的训练代码。看懂这份代码能够加深对SRGAN训练算法的理解。

SRGAN 核心思想

早期超分辨率方法的优化目标都是降低低清图像和高清图像之间的均方误差。降低均方误差，确实让增强图像和原高清图像的相似度更高。但是，图像的相似度指标高并不能代表图像的增强质量就很高。下图显示了插值、优化均方误差、SRGAN、原图这四个图像输出结果（括号里的相似度指标是PSNR和SSIM）。

从图中可以看出，优化均方误差虽然能让相似度指标升高，但图像的细节十分模糊，尤其是纹理比较密集的高频区域。相比之下，SRGAN增强出来的图像虽然相似度不高，但看起来更加清晰。

为什么SRGAN的增强结果那么清楚呢？这是因为SRGAN使用了一套新的优化目标。SRGAN使用的损失函数既包括了GAN误差，也包括了感知误差。这套新的优化目标能够让网络生成看起来更清楚的图片，而不仅仅是和原高清图像相似度更高的图片。

下面，我们来一步一步学习SRGAN的框架。

GAN 的原理

GAN[2]是一套搭建神经网络的框架。给定一个图片数据集$p_g$，GAN的目的是训练出一个生成网络$G$，使得G能够凭空生成出和$p_g$中大多数图片都类似的图片。比如说$p_g$是一个小猫图片数据集，那么$G$就应该能凭空生成出小猫图片。当然，$G$不是真的没有任何输入，真的能够凭空生成一幅图片。为了生成出不一样的图片，$G$要求输入一个随机量，这个随机量叫做噪声$z$。这样，只要输入的噪声$z$变了，$G$的输出$G(z)$就变了，就能画出长相不一样的小猫了。

为了指导图像生成，$G$应该有一个“老师”告诉它该怎么画出更像的图片。这个“老师”叫做判别网络$D$。$D$就是一个二分类网络，它能够严格地判定出一幅图片是否来自数据集$p_g$。如果$p_g$是一个小猫数据集，那么$D$就应该能判定一张图片是不是小猫。这样，如果$G$生成出来的图片$G(z)$已经非常逼真，连$D$都觉得$G(z)$来自数据集$p_g$，那么$G$就是一个很成功的网络了。

如果只是生成小猫，我们直接拿小猫图片和其他图片就能训练出一个$D$了。问题是，大多数情况下我们只有数据集$p_g$，而难以获得一个$p_g$的反例数据集。GAN的想法，则巧妙地解决了这个问题：刚开始，$G$生成出来的图片肯定是很差的，这些图片肯定不像$p_g$。所以，我们以$G(z)$为反例，和$p_g$一起训练出一个$D$来。等$D$的判定能力强了以后，又拿$D$回头训练$G$。这样，$D$的审美水平逐渐提高，$G$的绘画能力也逐渐提高。最终，$D$能成功分辨出一幅图片是否来自$p_g$，而$G$生成出来的图片和$p_g$中的看起来完全相同，连$D$也分辨不出来。就这样，我们得到了一个很棒的生成网络$G$。

规范地来说，给定一个数据集$p_g$，我们希望训练出两个网络$D, G$。$D$能够判断一幅输入图片是否来自$p_g$:

$G$则能够根据来自噪声分布$p_z$的$z$生成一个真假难辨的图片$G(z)$，使得$D(G(z))=1$。

为了达到这个目标，二分类器$D$应该最小化这样一个的交叉熵误差：

其中，$\hat{y}=D(x)$是预测结果为真的概率，$y$是0或1的标签。

对于来自数据集的图片$x \sim p_g$，$D$使用的标签$y$应该是1，误差公式化简为：

对于$G$生成的图片$G(z)$，$D$使用的标签$y$应该是0，误差公式化简为：

我们每步拿一张真图$x$和一张假图$G(z)$训练$D$。这样，每步的误差公式就是上面两个式子加起来：

反过来，$G$应该和$D$对抗，最大化上面那个误差，想办法骗过$D$。这个“对抗”就是GAN的名称“生成对抗网络”的由来。但是，$G$不能改变$D(x)$那一项。因此，$G$使用的误差函数是：

使用上面这两种误差，就可以训练神经网络了。训练GAN时，每轮一般会训练$k(k>=1)$次$D$，再训练1次$G$。这是为了先得到一个好的判别器，再用判别器去指导生成器。

GAN只是一套通用的框架，并没有指定神经网络$D, G$的具体结构。在不同任务中，$D, G$一般有不同的结构。

基于GAN的超分辨率网络

如前文所述，以优化均方误差为目标的超分辨率模型难以复原图像的细节。其实，超分辨率任务和图像生成任务类似，都需要一个“老师”来指导优化目标。SRGAN把GAN框架运用到了超分辨率任务上。原来的生成器$G$随机生成图像，现在用来输出高清图像；原来的判定器$D$用来判定图像是否属于某数据集，现在$D$用来判断一幅图像是否是高清图像。

具体来说，相比基础的GAN，在SRGAN中，$D$的输入是高清图像$I^{HR}$。而$G$的输入从随机噪声$z$变成了高清图像退化后的低清图像$I^{LR}$。这样，$G$就不是在随机生成图像，而是在根据一幅低清图像生成一幅高清图像了。它们的误差函数分别是：

借助GAN的架构，SRGAN能够利用$D$指导高清图像生成。但是，超分辨率任务毕竟和图像生成任务有一些区别，不能只用这种对抗误差来约束网络。因此，除了使用对抗误差外，SRGAN还使用了一种内容误差。这种内容误差用于让低清图片和高清图片的内容对齐，起到了和原均方误差一样的作用。

基于感知的内容误差

在介绍SRGAN的内容误差之前，需要对“内容误差”和“感知误差”这两个名词做一个澄清。在SRGAN的原文章中，作者把内容误差和对抗误差之和叫做感知误差。但是，后续的大部分文献只把这种内容误差叫做感知误差，不会把内容误差和对抗误差放在一起称呼。在后文中，我也会用“感知误差”来指代SRGAN中的“内容误差”。

在深度卷积神经网络（CNN）火起来后，人们开始研究为什么CNN能够和人类一样识别出图像。经实验，人们发现两幅图像经VGG（一个经典的CNN）的某些中间层的输出越相似，两幅图像从观感上也越相似。这种相似度并不是基于某种数学指标，而是和人的感知非常类似。

VGG的这种“感知性”被运用在了风格迁移等任务上。也有人考虑把这种感知上的误差运用到超分辨率任务上，并取得了不错的结果[3]。下图是真值、插值、基于逐像素误差、基于感知误差的四个超分辨率结果。

SRGAN也使用了这种感知误差，以取代之前常常使用的逐像素均方误差。这种感知误差的计算方法如下：VGG有很多中间层，用于计算感知误差的中间层$i$是可调的。假如我们用$\phi_{i}(I)$表示图像$I$经VGG的第$i$层的中间输出结果，$\phi_{i}(I)_{x, y}$表示中间输出结果在坐标$(x, y)$处的值，则感知误差的公式如下：

直观上解释这个公式，就是先把高清图像$I^{HR}$送入VGG，再把高清图像退化出来的低清图像$I^{LR}$送入生成器，并把生成器的输出$G(I^{LR})$也送入VGG。两幅图片经VGG第$i$层生成的中间结果的逐像素均方误差，就是感知误差。

算上之前的对抗误差，一个图像超分辨率网络的总误差如下：

这里的$w$用于调整两个误差的相对权重，原论文使用$w=10^{-3}$。

SRGAN的其他模块

定义好了误差函数，只要在决定好网络结构就可以开始训练网络了。SRGAN使用的生成网络和判别网络的结构如下：

判别网络就是一个平平无奇的二分类网络，架构上没有什么创新。而生成网络则先用几个残差块提取特征，最后用一种超分辨率任务中常用的上采样模块PixelShuffle对原图像的尺寸翻倍两次，最后输出一个边长放大4倍的高清图像。

SRGAN的这种网络结构在当时确实取得了不错的结果。但是，很快就有后续研究提出了更好的网络架构。比如ESRGAN[4]去掉了生成网络的BN层，提出了一种叫做RRDB的高级模块。基于RRDB的生成网络有着更好的生成效果。

不仅是网络架构，SRGAN的其他细节也得到了后续研究的改进。GAN误差的公式、总误差的公式、高清图像退化成低清图像的数据增强算法……这些子模块都被后续研究改进了。但是，SRGAN这种基于GAN的训练架构一直没有发生改变。有了SRGAN的代码，想复现一些更新的超分辨率网络时，往往只需要换一下生成器的结构，或者改一改误差的公式就行了。大部分的训练代码是不用改变的。

总结

SRGAN是把GAN运用在超分辨率任务上的开山之作。如正文所述，SRGAN中的部分设计虽然已经过时，但它的整体训练架构被一直沿用了下来。现在去回顾SRGAN这篇论文时，只需要关注以下几点即可:

• 如何把GAN套用在超分辨率任务上
• GAN误差
• 感知误差

通过阅读这篇论文，我们不仅应该学会GAN是怎样运用在SR上的，也应该能总结出如何把GAN应用在其他任务上。GAN的本质是去学习一个分布，令生成的$G(z)$看上去是来自分布$p_g$，而不是像图像分类等任务去学习一个$x \to y$的映射关系。因此，GAN会记忆一些和数据集相关的信息。在输入信息就已经比较完备的图像分类、目标检测等任务中，GAN可能没有什么用武之地。但是，在输入信息不足的超分辨率、图像补全等任务中，GAN记忆的数据集信息有很有用了。很多时候，GAN会“脑补”出输入图像中不够清楚的部分。

决定了要在某个任务中使用GAN时，我们可以在一个不使用GAN的架构上做以下改动：

• 定义一个分类网络$D$。
• 在原loss中加一项由$D$算出来的GAN loss。
• 在训练流程中，加入训练$D$的逻辑。

看完正文后，如果你对GAN在SR上的训练逻辑还是不太清楚，欢迎阅读附录中有关SRGAN训练代码的解读。

附录：MMEditing 中的 SRGAN

MMEditing中的SRGAN写在mmedit/models/restorers/srgan.py这个文件里。学习训练逻辑时，我们只需要关注SRGAN类的train_step方法即可。

以下是train_step的源代码（我的mmedit版本是v0.15.1）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

def train_step(self, data_batch, optimizer):
    """Train step.

    Args:
        data_batch (dict): A batch of data.
        optimizer (obj): Optimizer.

    Returns:
        dict: Returned output.
    """
    # data
    lq = data_batch['lq']
    gt = data_batch['gt']

    # generator
    fake_g_output = self.generator(lq)

    losses = dict()
    log_vars = dict()

    # no updates to discriminator parameters.
    set_requires_grad(self.discriminator, False)

    if (self.step_counter % self.disc_steps == 0
            and self.step_counter >= self.disc_init_steps):
        if self.pixel_loss:
            losses['loss_pix'] = self.pixel_loss(fake_g_output, gt)
        if self.perceptual_loss:
            loss_percep, loss_style = self.perceptual_loss(
                fake_g_output, gt)
            if loss_percep is not None:
                losses['loss_perceptual'] = loss_percep
            if loss_style is not None:
                losses['loss_style'] = loss_style
        # gan loss for generator
        fake_g_pred = self.discriminator(fake_g_output)
        losses['loss_gan'] = self.gan_loss(
            fake_g_pred, target_is_real=True, is_disc=False)

        # parse loss
        loss_g, log_vars_g = self.parse_losses(losses)
        log_vars.update(log_vars_g)

        # optimize
        optimizer['generator'].zero_grad()
        loss_g.backward()
        optimizer['generator'].step()

    # discriminator
    set_requires_grad(self.discriminator, True)
    # real
    real_d_pred = self.discriminator(gt)
    loss_d_real = self.gan_loss(
        real_d_pred, target_is_real=True, is_disc=True)
    loss_d, log_vars_d = self.parse_losses(dict(loss_d_real=loss_d_real))
    optimizer['discriminator'].zero_grad()
    loss_d.backward()
    log_vars.update(log_vars_d)
    # fake
    fake_d_pred = self.discriminator(fake_g_output.detach())
    loss_d_fake = self.gan_loss(
        fake_d_pred, target_is_real=False, is_disc=True)
    loss_d, log_vars_d = self.parse_losses(dict(loss_d_fake=loss_d_fake))
    loss_d.backward()
    log_vars.update(log_vars_d)

    optimizer['discriminator'].step()

    self.step_counter += 1

    log_vars.pop('loss')  # remove the unnecessary 'loss'
    outputs = dict(
        log_vars=log_vars,
        num_samples=len(gt.data),
        results=dict(lq=lq.cpu(), gt=gt.cpu(), output=fake_g_output.cpu()))

    return outputs

一开始，图像输出都在词典data_batch里。函数先把低清图lq和高清的真值gt从词典里取出。

1
2
3

# data
lq = data_batch['lq']
gt = data_batch['gt']

之后，函数计算了$G(I^{lq})$，为后续loss的计算做准备。

1
2

# generator
fake_g_output = self.generator(lq)

接下来，是优化生成器self.generator的逻辑。这里面有一些函数调用，我们可以不管它们的实现，大概理解整段代码的意思就行了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

losses = dict()
log_vars = dict()

# no updates to discriminator parameters.
set_requires_grad(self.discriminator, False)

if (self.step_counter % self.disc_steps == 0
        and self.step_counter >= self.disc_init_steps):
    if self.pixel_loss:
        losses['loss_pix'] = self.pixel_loss(fake_g_output, gt)
    if self.perceptual_loss:
        loss_percep, loss_style = self.perceptual_loss(
            fake_g_output, gt)
        if loss_percep is not None:
            losses['loss_perceptual'] = loss_percep
        if loss_style is not None:
            losses['loss_style'] = loss_style
    # gan loss for generator
    fake_g_pred = self.discriminator(fake_g_output)
    losses['loss_gan'] = self.gan_loss(
        fake_g_pred, target_is_real=True, is_disc=False)

    # parse loss
    loss_g, log_vars_g = self.parse_losses(losses)
    log_vars.update(log_vars_g)

    # optimize
    optimizer['generator'].zero_grad()
    loss_g.backward()
    optimizer['generator'].step()

为了只训练生成器，要用下面的代码关闭判别器的训练。

1
2

# no updates to discriminator parameters.
set_requires_grad(self.discriminator, False)

正文说过，训练GAN时一般要先训好判别器，且训练判别器多于训练生成器。因此，下面的if语句可以让判别器训练了self.disc_init_steps步后，每训练self.disc_steps步判别器再训练一步生成器。

1
2

if (self.step_counter % self.disc_steps == 0
    and self.step_counter >= self.disc_init_steps):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

if self.pixel_loss:
    losses['loss_pix'] = self.pixel_loss(fake_g_output, gt)
if self.perceptual_loss:
    loss_percep, loss_style = self.perceptual_loss(
        fake_g_output, gt)
    if loss_percep is not None:
        losses['loss_perceptual'] = loss_percep
    if loss_style is not None:
        losses['loss_style'] = loss_style
# gan loss for generator
fake_g_pred = self.discriminator(fake_g_output)
losses['loss_gan'] = self.gan_loss(
    fake_g_pred, target_is_real=True, is_disc=False)

# parse loss
loss_g, log_vars_g = self.parse_losses(losses)
log_vars.update(log_vars_g)

# optimize
optimizer['generator'].zero_grad()
loss_g.backward()
optimizer['generator'].step()

训练完生成器后，要训练判别器。和生成器的误差计算方法类似，判别器的训练代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

 # discriminator
set_requires_grad(self.discriminator, True)
# real
real_d_pred = self.discriminator(gt)
loss_d_real = self.gan_loss(
    real_d_pred, target_is_real=True, is_disc=True)
loss_d, log_vars_d = self.parse_losses(dict(loss_d_real=loss_d_real))
optimizer['discriminator'].zero_grad()
loss_d.backward()
log_vars.update(log_vars_d)
# fake
fake_d_pred = self.discriminator(fake_g_output.detach())
loss_d_fake = self.gan_loss(
    fake_d_pred, target_is_real=False, is_disc=True)
loss_d, log_vars_d = self.parse_losses(dict(loss_d_fake=loss_d_fake))
loss_d.backward()
log_vars.update(log_vars_d)

optimizer['discriminator'].step()

这段代码有两个重点：

1. 在训练判别器时，要用set_requires_grad(self.discriminator, True)开启判别器的梯度计算。
2. fake_d_pred = self.discriminator(fake_g_output.detach())这一行的detach()很关键。detach()可以中断某张量的梯度跟踪。fake_g_output是由生成器算出来的，如果不把这个张量的梯度跟踪切断掉，在优化判别器时生成器的参数也会跟着优化。

函数的最后部分是一些和MMEditing其他代码逻辑的交互，和SRGAN本身没什么关联。

1
2
3
4
5
6
7
8
9

self.step_counter += 1

log_vars.pop('loss')  # remove the unnecessary 'loss'
outputs = dict(
    log_vars=log_vars,
    num_samples=len(gt.data),
    results=dict(lq=lq.cpu(), gt=gt.cpu(), output=fake_g_output.cpu()))

return outputs

只要理解了本文的误差计算公式，再看懂了这段代码是如何训练判别器和生成器的，就算是完全理解了SRGAN的核心思想了。

参考资料

[1] (SRGAN): Photo-Realistic Single Image Super-Resolution Using a Generative Adversarial Network

[2] (GAN): Generative Adversarial Nets

[3] (Perceptual Loss)：Perceptual Losses for Real-Time Style Transfer and Super-Resolution

[4] (ESRGAN): ESRGAN: Enhanced Super-Resolution Generative Adversarial Networks

我之前的一篇文章介绍了如何给PyTorch添加CPU上的简单的加法算子。在这篇文章里，我将继续展示一个更具体的PyTorch自定义算子示例——自己动手复现二维卷积算子。这个示例是基于PyTorch Extension的，在迁移项目时，不需要自己生成动态库，只需要用setup.py重新编译一遍即可。我会同时介绍CPU版和CUDA版的实现。

许多前沿的神经网络都会对卷积进行一些修改。比如大名鼎鼎的可变形卷积(deformable convolution)。相信看完这篇文章后，大家能看懂PyTorch卷积的实现代码，并大概了解如何修改卷积的实现细节，并把新写好的卷积运用到自己的PyTorch项目中。

PyTorch Extension 实现二维卷积

搭建项目

在开始写代码前，要准备一个崭新的目录，在这个文件夹里搭建项目。

在根目录下，先创建一个setup.py，之后要填写这份安装文件。

之后，创建一个文件夹，其名字是项目名。在这个文件夹里合适的地方新建一个子文件夹，专门用来放和算子相关的文件。我的项目名叫做panoflow，算子相关文件放在了panoflow/core/op子文件夹下。

接下来，和算子实现相关的文件都应该放在算子文件夹里。使用和测试算子的文件可以放在项目文件夹的其他地方。

由于在实现中我借用了MMCV的代码，还要提前准备好一些头文件。首先新建一个文件pytorch_cpp_helper.hpp：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

#ifndef PYTORCH_CPP_HELPER
#define PYTORCH_CPP_HELPER
#include <torch/extension.h>

#include <vector>

using namespace at;

#define CHECK_CUDA(x) \
  TORCH_CHECK(x.device().is_cuda(), #x " must be a CUDA tensor")
#define CHECK_CPU(x) \
  TORCH_CHECK(!x.device().is_cuda(), #x " must be a CPU tensor")
#define CHECK_CONTIGUOUS(x) \
  TORCH_CHECK(x.is_contiguous(), #x " must be contiguous")
#define CHECK_CUDA_INPUT(x) \
  CHECK_CUDA(x);            \
  CHECK_CONTIGUOUS(x)
#define CHECK_CPU_INPUT(x) \
  CHECK_CPU(x);            \
  CHECK_CONTIGUOUS(x)

#endif  // PYTORCH_CPP_HELPER

再创建一个文件pytorch_cuda_helper.hpp：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

#ifndef PYTORCH_CUDA_HELPER
#define PYTORCH_CUDA_HELPER

#include <ATen/ATen.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>

#include <ATen/cuda/CUDAApplyUtils.cuh>
#include <THC/THCAtomics.cuh>

#include "common_cuda_helper.hpp"

using at::Half;
using at::Tensor;
using phalf = at::Half;

#define __PHALF(x) (x)

#endif  // PYTORCH_CUDA_HELPER

还有一个common_cuda_helper.hpp：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120

#ifndef COMMON_CUDA_HELPER
#define COMMON_CUDA_HELPER

#include <cuda.h>

#define CUDA_1D_KERNEL_LOOP(i, n)                              \
  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < (n); \
       i += blockDim.x * gridDim.x)

#define CUDA_2D_KERNEL_LOOP(i, n, j, m)                             \
  for (size_t i = blockIdx.x * blockDim.x + threadIdx.x; i < (n);   \
       i += blockDim.x * gridDim.x)                                 \
    for (size_t j = blockIdx.y * blockDim.y + threadIdx.y; j < (m); \
         j += blockDim.y * gridDim.y)

#define CUDA_2D_KERNEL_BLOCK_LOOP(i, n, j, m)          \
  for (size_t i = blockIdx.x; i < (n); i += gridDim.x) \
    for (size_t j = blockIdx.y; j < (m); j += gridDim.y)

#define THREADS_PER_BLOCK 512

inline int GET_BLOCKS(const int N, const int num_threads = THREADS_PER_BLOCK) {
  int optimal_block_num = (N + num_threads - 1) / num_threads;
  int max_block_num = 4096;
  return min(optimal_block_num, max_block_num);
}

template <typename T>
__device__ T bilinear_interpolate(const T* input, const int height,
                                  const int width, T y, T x,
                                  const int index /* index for debug only*/) {
  // deal with cases that inverse elements are out of feature map boundary
  if (y < -1.0 || y > height || x < -1.0 || x > width) return 0;

  if (y <= 0) y = 0;
  if (x <= 0) x = 0;

  int y_low = (int)y;
  int x_low = (int)x;
  int y_high;
  int x_high;

  if (y_low >= height - 1) {
    y_high = y_low = height - 1;
    y = (T)y_low;
  } else {
    y_high = y_low + 1;
  }

  if (x_low >= width - 1) {
    x_high = x_low = width - 1;
    x = (T)x_low;
  } else {
    x_high = x_low + 1;
  }

  T ly = y - y_low;
  T lx = x - x_low;
  T hy = 1. - ly, hx = 1. - lx;
  // do bilinear interpolation
  T v1 = input[y_low * width + x_low];
  T v2 = input[y_low * width + x_high];
  T v3 = input[y_high * width + x_low];
  T v4 = input[y_high * width + x_high];
  T w1 = hy * hx, w2 = hy * lx, w3 = ly * hx, w4 = ly * lx;

  T val = (w1 * v1 + w2 * v2 + w3 * v3 + w4 * v4);

  return val;
}

template <typename T>
__device__ void bilinear_interpolate_gradient(
    const int height, const int width, T y, T x, T& w1, T& w2, T& w3, T& w4,
    int& x_low, int& x_high, int& y_low, int& y_high,
    const int index /* index for debug only*/) {
  // deal with cases that inverse elements are out of feature map boundary
  if (y < -1.0 || y > height || x < -1.0 || x > width) {
    // empty
    w1 = w2 = w3 = w4 = 0.;
    x_low = x_high = y_low = y_high = -1;
    return;
  }

  if (y <= 0) y = 0;
  if (x <= 0) x = 0;

  y_low = (int)y;
  x_low = (int)x;

  if (y_low >= height - 1) {
    y_high = y_low = height - 1;
    y = (T)y_low;
  } else {
    y_high = y_low + 1;
  }

  if (x_low >= width - 1) {
    x_high = x_low = width - 1;
    x = (T)x_low;
  } else {
    x_high = x_low + 1;
  }

  T ly = y - y_low;
  T lx = x - x_low;
  T hy = 1. - ly, hx = 1. - lx;

  // reference in forward
  // T v1 = input[y_low * width + x_low];
  // T v2 = input[y_low * width + x_high];
  // T v3 = input[y_high * width + x_low];
  // T v4 = input[y_high * width + x_high];
  // T val = (w1 * v1 + w2 * v2 + w3 * v3 + w4 * v4);

  w1 = hy * hx, w2 = hy * lx, w3 = ly * hx, w4 = ly * lx;

  return;
}
#endif  // COMMON_CUDA_HELPER

这些文件添加了CPU和CUDA实现时需要的头文件和定义，后面的C++源码会用到它们。

CPU

C++实现

在用C++实现一个算子时，我们要编写一个形如这样的文件：

1
2
3
4
5
6
7
8
9
10
11

#include <torch/torch.h>

torch::Tensor my_add(torch::Tensor t1, torch::Tensor t2)
{
	return t1 + t2;
}

TORCH_LIBRARY(my_ops, m)
{
	m.def("my_add", my_add);
}

这个C++文件主要包含两部分内容：算子的实现函数和C++接口绑定。在实现卷积时，也是要实现这两部分内容。

在修改一个现有的算子时，最好的方法不是从头写一个，而是去开源库里找一份实现，并在这个基础上进行修改。

我在MMCV的仓库里找到了可变形卷积的实现，并把它拆解回了普通的卷积。我参考了这篇教程：手把手教你如何高效地在 MMCV 中贡献算子。另外，这份笔记还参考了PyTorch官方Extension教程。

找到了卷积的实现后，在算子文件夹下新建一个cpp源文件。比如我的文件路径就是panoflow/core/op/my_conv.cpp。这样一个普通卷积的实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298

#include "pytorch_cpp_helper.hpp"

void my_conv_im2col_cpu(Tensor data_im,
                        const int channels, const int height,
                        const int width, const int ksize_h,
                        const int ksize_w, const int pad_h, const int pad_w,
                        const int stride_h, const int stride_w,
                        const int dilation_h, const int dilation_w,
                        const int parallel_imgs, Tensor data_col);

void my_conv_im2col_cuda(Tensor data_im,
                        const int channels, const int height,
                        const int width, const int ksize_h,
                        const int ksize_w, const int pad_h, const int pad_w,
                        const int stride_h, const int stride_w,
                        const int dilation_h, const int dilation_w,
                        const int parallel_imgs, Tensor data_col);

void my_conv_shape_check(at::Tensor input,
                         at::Tensor weight, int kH,
                         int kW, int dH, int dW, int padH, int padW,
                         int dilationH, int dilationW, int group)
{
    TORCH_CHECK(
        weight.ndimension() == 4,
        "4D weight tensor (nOutputPlane,nInputPlane,kH,kW) expected, but got: %s",
        weight.ndimension());

    TORCH_CHECK(weight.is_contiguous(), "weight tensor has to be contiguous");

    TORCH_CHECK(kW > 0 && kH > 0,
                "kernel size should be greater than zero, but got kH: %d kW: %d",
                kH, kW);

    TORCH_CHECK((weight.size(2) == kH && weight.size(3) == kW),
                "kernel size should be consistent with weight, ",
                "but got kH: %d kW: %d weight.size(2): %d, weight.size(3): %d",
                kH, kW, weight.size(2), weight.size(3));

    TORCH_CHECK(dW > 0 && dH > 0,
                "stride should be greater than zero, but got dH: %d dW: %d", dH,
                dW);

    TORCH_CHECK(
        dilationW > 0 && dilationH > 0,
        "dilation should be greater than 0, but got dilationH: %d dilationW: %d",
        dilationH, dilationW);

    int ndim = input.ndimension();
    int dimf = 0;
    int dimh = 1;
    int dimw = 2;

    if (ndim == 4)
    {
        dimf++;
        dimh++;
        dimw++;
    }

    TORCH_CHECK(ndim == 3 || ndim == 4,
                "3D or 4D input tensor expected but got: %s", ndim);

    long nInputPlane = weight.size(1) * group;
    long inputHeight = input.size(dimh);
    long inputWidth = input.size(dimw);
    long nOutputPlane = weight.size(0);
    long outputHeight =
        (inputHeight + 2 * padH - (dilationH * (kH - 1) + 1)) / dH + 1;
    long outputWidth =
        (inputWidth + 2 * padW - (dilationW * (kW - 1) + 1)) / dW + 1;

    if (outputWidth < 1 || outputHeight < 1)
        AT_ERROR(
            "Given input size: (%ld x %ld x %ld). "
            "Calculated output size: (%ld x %ld x %ld). Output size is too small",
            nInputPlane, inputHeight, inputWidth, nOutputPlane, outputHeight,
            outputWidth);

    TORCH_CHECK(input.size(1) == nInputPlane,
                "invalid number of input planes, expected: %d, but got: %d",
                nInputPlane, input.size(1));

    TORCH_CHECK((inputHeight >= kH && inputWidth >= kW),
                "input image is smaller than kernel");
}

void my_conv_forward(Tensor input, Tensor weight, Tensor bias,
                     Tensor output, Tensor columns, int kW,
                     int kH, int dW, int dH, int padW, int padH,
                     int dilationW, int dilationH, int group,
                     int im2col_step)
{
    bool isCuda = false;
    if (input.device().is_cuda())
    {
        CHECK_CUDA_INPUT(input);
        CHECK_CUDA_INPUT(weight);
        CHECK_CUDA_INPUT(bias);
        CHECK_CUDA_INPUT(output);
        CHECK_CUDA_INPUT(columns);
        isCuda = true;
    }
    else
    {
        CHECK_CPU_INPUT(input);
        CHECK_CPU_INPUT(weight);
        CHECK_CPU_INPUT(bias);
        CHECK_CPU_INPUT(output);
        CHECK_CPU_INPUT(columns);
    }

    my_conv_shape_check(input, weight, kH, kW, dH, dW, padH,
                        padW, dilationH, dilationW, group);
    at::DeviceGuard guard(input.device());

    int batch = 1;
    if (input.ndimension() == 3)
    {
        // Force batch
        batch = 0;
        input.unsqueeze_(0);
    }

    long batchSize = input.size(0);
    long nInputPlane = input.size(1);
    long inputHeight = input.size(2);
    long inputWidth = input.size(3);

    long nOutputPlane = weight.size(0);

    long outputWidth =
        (inputWidth + 2 * padW - (dilationW * (kW - 1) + 1)) / dW + 1;
    long outputHeight =
        (inputHeight + 2 * padH - (dilationH * (kH - 1) + 1)) / dH + 1;

    output = output.view({batchSize / im2col_step, im2col_step, nOutputPlane,
                          outputHeight, outputWidth});
    columns = at::zeros(
        {nInputPlane * kW * kH, im2col_step * outputHeight * outputWidth},
        input.options());

    input = input.view({batchSize / im2col_step, im2col_step, nInputPlane,
                        inputHeight, inputWidth});

    Tensor output_buffer = at::zeros({batchSize / im2col_step, nOutputPlane,
                                      im2col_step * outputHeight, outputWidth},
                                     output.options());

    output_buffer = output_buffer.view(
        {output_buffer.size(0), group, output_buffer.size(1) / group,
         output_buffer.size(2), output_buffer.size(3)});

    for (int elt = 0; elt < batchSize / im2col_step; elt++)
    {
        if (isCuda)
        {
            my_conv_im2col_cuda(input[elt], nInputPlane, inputHeight,
                            inputWidth, kH, kW, padH, padW, dH, dW, dilationH,
                            dilationW, im2col_step, columns);
        }
        else
        {
            my_conv_im2col_cpu(input[elt], nInputPlane, inputHeight,
                            inputWidth, kH, kW, padH, padW, dH, dW, dilationH,
                            dilationW, im2col_step, columns);
        }
        

        columns = columns.view({group, columns.size(0) / group, columns.size(1)});
        weight = weight.view({group, weight.size(0) / group, weight.size(1),
                              weight.size(2), weight.size(3)});

        for (int g = 0; g < group; g++)
        {
            output_buffer[elt][g] = output_buffer[elt][g]
                                        .flatten(1)
                                        .addmm_(weight[g].flatten(1), columns[g])
                                        .view_as(output_buffer[elt][g]);
        }
        columns =
            columns.view({columns.size(0) * columns.size(1), columns.size(2)});
        weight = weight.view({weight.size(0) * weight.size(1), weight.size(2),
                              weight.size(3), weight.size(4)});
    }

    output_buffer = output_buffer.view(
        {output_buffer.size(0), output_buffer.size(1) * output_buffer.size(2),
         output_buffer.size(3), output_buffer.size(4)});

    output_buffer = output_buffer.view({batchSize / im2col_step, nOutputPlane,
                                        im2col_step, outputHeight, outputWidth});
    output_buffer.transpose_(1, 2);

    output.copy_(output_buffer);
    output = output.view({batchSize, nOutputPlane, outputHeight, outputWidth});

    bias = bias.view({1, bias.size(0), 1, 1});
    output.add_(bias);

    input = input.view({batchSize, nInputPlane, inputHeight, inputWidth});

    if (batch == 0)
    {
        output = output.view({nOutputPlane, outputHeight, outputWidth});
        input = input.view({nInputPlane, inputHeight, inputWidth});
    }
}

template <typename T>
void my_conv_im2col_cpu_kernel(
    const int n, const T *data_im, 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,
    const int dilation_h, const int dilation_w,
    const int batch_size,
    const int num_channels, const int height_col,
    const int width_col, T *data_col)
{
    for (int index = 0; index < n; index++)
    {
        // index index of output matrix
        const int w_col = index % width_col;
        const int h_col = (index / width_col) % height_col;
        const int b_col = (index / width_col / height_col) % batch_size;
        const int c_im = (index / width_col / height_col) / batch_size;
        const int c_col = c_im * kernel_h * kernel_w;

        const int h_in = h_col * stride_h - pad_h;
        const int w_in = w_col * stride_w - pad_w;
        T *data_col_ptr =
            data_col +
            ((c_col * batch_size + b_col) * height_col + h_col) * width_col + w_col;
        const T *data_im_ptr =
            data_im + (b_col * num_channels + c_im) * height * width;

        for (int i = 0; i < kernel_h; ++i)
        {
            for (int j = 0; j < kernel_w; ++j)
            {
                T val = static_cast<T>(0);
                const int h_im = h_in + i * dilation_h;
                const int w_im = w_in + j * dilation_w;
                if (h_im > -1 && w_im > -1 && h_im < height && w_im < width)
                {
                    val = data_im_ptr[h_im * width + w_im];
                }
                *data_col_ptr = val;
                data_col_ptr += batch_size * height_col * width_col;
            }
        }
    }
}

void my_conv_im2col_cpu(Tensor data_im,
                        const int channels, const int height,
                        const int width, const int ksize_h,
                        const int ksize_w, const int pad_h, const int pad_w,
                        const int stride_h, const int stride_w,
                        const int dilation_h, const int dilation_w,
                        const int parallel_imgs, Tensor data_col)
{
    int height_col =
        (height + 2 * pad_h - (dilation_h * (ksize_h - 1) + 1)) / stride_h + 1;
    int width_col =
        (width + 2 * pad_w - (dilation_w * (ksize_w - 1) + 1)) / stride_w + 1;
    int num_kernels = channels * height_col * width_col * parallel_imgs;

    AT_DISPATCH_FLOATING_TYPES_AND_HALF(
        data_im.scalar_type(), "", [&]
        { my_conv_im2col_cpu_kernel<scalar_t>(
              num_kernels, data_im.data_ptr<scalar_t>(),
              height, width, ksize_h, ksize_w,
              pad_h, pad_w, stride_h, stride_w, dilation_h, dilation_w,
              parallel_imgs, channels,
              height_col, width_col, data_col.data_ptr<scalar_t>()); });
}

void my_conv_im2col_cuda(Tensor data_im,
                        const int channels, const int height,
                        const int width, const int ksize_h,
                        const int ksize_w, const int pad_h, const int pad_w,
                        const int stride_h, const int stride_w,
                        const int dilation_h, const int dilation_w,
                        const int parallel_imgs, Tensor data_col)
{
    
}

PYBIND11_MODULE(my_ops, m)
{
      m.def("my_conv_forward", my_conv_forward, "my_conv_forward",
            py::arg("input"), py::arg("weight"), py::arg("bias"),
            py::arg("output"), py::arg("columns"), py::arg("kW"),
            py::arg("kH"), py::arg("dW"), py::arg("dH"), py::arg("padW"),
            py::arg("padH"), py::arg("dilationW"), py::arg("dilationH"),
            py::arg("group"), py::arg("im2col_step"));
}

这份实现非常长，我挑一些重点的内容讲解。

从最下面的PYBIND11_MODULE(my_ops, m)看起。这里的my_ops是生成的库名，可以随便取名。待会要import这个库名。代码块里m.def用于定义C++函数的Python接口。"my_conv_forward"是Python调用时的函数名称，my_conv_forward是被Python代码调用的这份代码里的C++函数名称。也就是说，这份卷积实现的入口函数就是my_conv_forward。我们从这个函数看起。

1
2
3
4
5

void my_conv_forward(Tensor input, Tensor weight, Tensor bias,
                     Tensor output, Tensor columns, int kW,
                     int kH, int dW, int dH, int padW, int padH,
                     int dilationW, int dilationH, int group,
                     int im2col_step)

my_conv_forward就是卷积的主函数。它的参数除了PyTorch的Conv2d传入的参数外，还多了两个参数output, columus。这两个张量是保存中间结果的，在PyTorch侧是看不到的。output用于保存卷积输出，columns用于保存卷积时的列矩阵。底层实现卷积时，会先把图像转换成一个用列表示的矩阵，再把卷积操作当成一个矩阵乘法来完成。其中，第一步操作叫做”im2col”。对此原理不熟的话可以参考这篇文章：https://zhuanlan.zhihu.com/p/63974249。

my_conv_forward函数的大部分内容都是在做类型检查和张量形状转换。在修改卷积实现时，这些东西都可以不用改。整个卷积操作的核心都在这一部分：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

for (int elt = 0; elt < batchSize / im2col_step; elt++)
{
    if (isCuda)
    {
        my_conv_im2col_cuda(input[elt], nInputPlane, inputHeight,
                        inputWidth, kH, kW, padH, padW, dH, dW, dilationH,
                        dilationW, im2col_step, columns);
    }
    else
    {
        my_conv_im2col_cpu(input[elt], nInputPlane, inputHeight,
                        inputWidth, kH, kW, padH, padW, dH, dW, dilationH,
                        dilationW, im2col_step, columns);
    }
    

    columns = columns.view({group, columns.size(0) / group, columns.size(1)});
    weight = weight.view({group, weight.size(0) / group, weight.size(1),
                          weight.size(2), weight.size(3)});

    for (int g = 0; g < group; g++)
    {
        output_buffer[elt][g] = output_buffer[elt][g]
                                    .flatten(1)
                                    .addmm_(weight[g].flatten(1), columns[g])
                                    .view_as(output_buffer[elt][g]);
    }
    columns =
        columns.view({columns.size(0) * columns.size(1), columns.size(2)});
    weight = weight.view({weight.size(0) * weight.size(1), weight.size(2),
                          weight.size(3), weight.size(4)});
}

这段代码先做了im2col操作，再做了矩阵乘法。其实，包括可变形卷积在内，各种稀奇古怪的卷积操作通过靠修改im2col来完成的。CPU和CUDA版卷积的主要区别，也体现在im2col中（后面的矩阵乘法在CPU和CUDA上都能用）。

由于是讲CPU实现，这里的CUDA实现我暂时放了一个空函数。my_conv_im2col_cpu的内容如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

void my_conv_im2col_cpu(Tensor data_im,
                        const int channels, const int height,
                        const int width, const int ksize_h,
                        const int ksize_w, const int pad_h, const int pad_w,
                        const int stride_h, const int stride_w,
                        const int dilation_h, const int dilation_w,
                        const int parallel_imgs, Tensor data_col)
{
    int height_col =
        (height + 2 * pad_h - (dilation_h * (ksize_h - 1) + 1)) / stride_h + 1;
    int width_col =
        (width + 2 * pad_w - (dilation_w * (ksize_w - 1) + 1)) / stride_w + 1;
    int num_kernels = channels * height_col * width_col * parallel_imgs;

    AT_DISPATCH_FLOATING_TYPES_AND_HALF(
        data_im.scalar_type(), "", [&]
        { my_conv_im2col_cpu_kernel<scalar_t>(
              num_kernels, data_im.data_ptr<scalar_t>(),
              height, width, ksize_h, ksize_w,
              pad_h, pad_w, stride_h, stride_w, dilation_h, dilation_w,
              parallel_imgs, channels,
              height_col, width_col, data_col.data_ptr<scalar_t>()); });
}

这个函数其实只是处理了一下输入，真正的实现在my_conv_im2col_cpu_kernel里。AT_DISPATCH_FLOATING_TYPES_AND_HALF可以让实现兼容半精度和普通float，所以实现my_conv_im2col_cpu_kernel得写成一个模板函数。

my_conv_im2col_cpu_kernel的实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

template <typename T>
void my_conv_im2col_cpu_kernel(
    const int n, const T *data_im, 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,
    const int dilation_h, const int dilation_w,
    const int batch_size,
    const int num_channels, const int height_col,
    const int width_col, T *data_col)
{
    for (int index = 0; index < n; index++)
    {
        // index index of output matrix
        const int w_col = index % width_col;
        const int h_col = (index / width_col) % height_col;
        const int b_col = (index / width_col / height_col) % batch_size;
        const int c_im = (index / width_col / height_col) / batch_size;
        const int c_col = c_im * kernel_h * kernel_w;

        const int h_in = h_col * stride_h - pad_h;
        const int w_in = w_col * stride_w - pad_w;
        T *data_col_ptr =
            data_col +
            ((c_col * batch_size + b_col) * height_col + h_col) * width_col + w_col;
        const T *data_im_ptr =
            data_im + (b_col * num_channels + c_im) * height * width;

        for (int i = 0; i < kernel_h; ++i)
        {
            for (int j = 0; j < kernel_w; ++j)
            {
                T val = static_cast<T>(0);
                const int h_im = h_in + i * dilation_h;
                const int w_im = w_in + j * dilation_w;
                if (h_im > -1 && w_im > -1 && h_im < height && w_im < width)
                {
                    val = data_im_ptr[h_im * width + w_im];
                }
                *data_col_ptr = val;
                data_col_ptr += batch_size * height_col * width_col;
            }
        }
    }
}

Python封装

实现好了后，如果编译完了的话，刚刚的卷积接口可以通过以下方式在Python里调用：

1
2

import my_ops
my_ops.my_conv_forward(...)

这里的my_ops这个名称必须和开始PYBIND11_MODULE(my_ops, m)里面那个库名称对应。

基于这个接口，可以仿照PyTorch中Conv2d的接口，编写一个和Conv2d等价的torch.nn.Module出来。我的这个Python文件的路径是panoflow/core/op/my_conv.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

import torch
from torch.autograd import Function
import torch.nn as nn
from torch import Tensor
from torch.nn.modules.utils import _pair
from torch.nn.parameter import Parameter

import my_ops


class MyConvF(Function):

    @staticmethod
    def forward(ctx,
                input: torch.Tensor,
                weight,
                bias,
                stride=1,
                padding=0,
                dilation=1,
                groups=1,
                im2col_step=32):
        if input is not None and input.dim() != 4:
            raise ValueError(
                f'Expected 4D tensor as input, got {input.dim()}D tensor \
                  instead.')
        ctx.stride = _pair(stride)
        ctx.padding = _pair(padding)
        ctx.dilation = _pair(dilation)
        ctx.groups = groups
        ctx.im2col_step = im2col_step

        weight = weight.type_as(input)
        ctx.save_for_backward(input, weight)

        output = input.new_empty(MyConvF._output_size(ctx, input, weight))

        ctx.bufs_ = [input.new_empty(0), input.new_empty(0)]  # columns, ones

        cur_im2col_step = min(ctx.im2col_step, input.size(0))
        assert (input.size(0) % cur_im2col_step
                ) == 0, 'batch size must be divisible by im2col_step'

        my_ops.my_conv_forward(
            input,
            weight,
            bias,
            output,
            ctx.bufs_[0],
            kW=weight.size(3),
            kH=weight.size(2),
            dW=ctx.stride[1],
            dH=ctx.stride[0],
            padW=ctx.padding[1],
            padH=ctx.padding[0],
            dilationW=ctx.dilation[1],
            dilationH=ctx.dilation[0],
            group=ctx.groups,
            im2col_step=cur_im2col_step)
        return output

    @staticmethod
    def _output_size(ctx, input, weight):
        channels = weight.size(0)
        output_size = (input.size(0), channels)
        for d in range(input.dim() - 2):
            in_size = input.size(d + 2)
            pad = ctx.padding[d]
            kernel = ctx.dilation[d] * (weight.size(d + 2) - 1) + 1
            stride_ = ctx.stride[d]
            output_size += ((in_size + (2 * pad) - kernel) // stride_ + 1, )
        if not all(map(lambda s: s > 0, output_size)):
            raise ValueError(
                'convolution input is too small (output would be ' +
                'x'.join(map(str, output_size)) + ')')
        return output_size


my_conv = MyConvF.apply


class MyConv2d(nn.Module):

    def __init__(self,
                 in_channels: int,
                 out_channels: int,
                 kernel_size,
                 stride=1,
                 padding=0,
                 dilation=1,
                 groups: int = 1,
                 bias: bool = True):
        super().__init__()
        kernel_size_ = _pair(kernel_size)
        stride_ = _pair(stride)
        padding_ = _pair(padding)
        dilation_ = _pair(dilation)
        self.in_channels = in_channels
        self.out_channels = out_channels
        self.kernel_size = kernel_size_
        self.stride = stride_
        self.padding = padding_
        self.dilation = dilation_
        self.groups = groups
        self.weight = Parameter(
            torch.Tensor(out_channels, in_channels // groups, *kernel_size_))
        self.bias = Parameter(torch.Tensor(out_channels))

        # Useless attributes
        self.transposed = None
        self.output_padding = None
        self.padding_mode = None

    def forward(self, input: Tensor) -> Tensor:
        return my_conv(input, self.weight, self.bias, self.stride,
                       self.padding, self.dilation, self.groups)

以后，用自己的卷积MyConv2d就和用普通的Conv2d一样了。

编译

打开外面的setup.py，填写以下内容。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

from setuptools import setup
from torch.utils import cpp_extension
import os

src_root = 'panoflow/core/op'
cpp_src = ['my_conv.cpp']

if __name__ == '__main__':
    include_dirs = ['panoflow/core/op']
    cpp_path = [os.path.join(src_root, src) for src in cpp_src]

    setup(
        name='panoflow',
        ext_modules=[
            cpp_extension.CppExtension(
                'my_ops', cpp_path, include_dirs=include_dirs)
        ],
        cmdclass={'build_ext': cpp_extension.BuildExtension})

其中的路径要根据自己的实际情况修改。

和编译相关的内容都写在cpp_extension.CppExtension里。其中，源文件要写在第二个参数里，头文件目录要写在include_dirs。由于我的源文件放在panoflow/core/op里，我写了个源文件名数组cpp_src，在传参前把路径组合了一下。由于include_dirs和源文件在同一个目录下，我也填的是panoflow/core/op。

写完了setup.py后，运行python setup.py develop，就能一键编译和安装。如果运行后没有报编译错误，就可以把实现的卷积用起来了。

单元测试

用单元测试可以快速地验证卷积是否实现成功。我写了一个简单的单元测试文件，在任意一个文件夹下创建该文件即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

import torch
import torch.nn as nn
from panoflow.core.op.my_conv import MyConv2d

inc = 3
outc = 4
img_shaspe = (50, 50)

# device_name = 'cuda:0'
device_name = 'cpu'
open_bias = True


def test_one():
    ts = torch.ones([1, 1, 3, 3]).to(device_name)
    layer = nn.Conv2d(1, 1, 3, 1, 1, bias=open_bias).to(device_name)
    gt = layer(ts)
    my_layer = MyConv2d(1, 1, 3, 1, 1).to(device_name)
    my_layer.load_state_dict(layer.state_dict(), strict=False)
    res = my_layer(ts)
    res = res.to('cpu')
    gt = gt.to('cpu')
    assert torch.allclose(res, gt, 1e-3, 1e-5)


def test_two():
    ts = torch.rand([1, inc, *img_shaspe]).to(device_name)
    layer = nn.Conv2d(inc, outc, 3, 1, 1, bias=open_bias).to(device_name)
    gt = layer(ts)
    my_layer = MyConv2d(inc, outc, 3, 1, 1).to(device_name)
    my_layer.load_state_dict(layer.state_dict(), strict=False)
    res = my_layer(ts)
    res = res.to('cpu')
    gt = gt.to('cpu')
    assert torch.allclose(res, gt, 1e-3, 1e-5)


if __name__ == '__main__':
    test_one()
    test_two()

其中，panoflow.core.op.my_conv是我刚刚放MyConv2d的Python模块。

直接运行这个Python文件，如果没有任何输出（报错信息），就说明卷积实现成功了。

CUDA

C++实现

在刚刚的实现中，有一个my_conv_im2col_cuda的实现是空着的。在CUDA版本中，我们要实现这个函数。不过，这个函数要放在一个用nvcc编译的.cu文件里。注意！注意！注意！ 因此，my_conv.cpp里那个空的my_conv_im2col_cuda实现应该全部删掉。

新建一个文件my_conv_cuda.cu。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

// Modify from https://github.com/open-mmlab/mmcv/blob/my_conv/mmcv/ops/csrc/common/cuda/deform_conv_cuda_kernel.cuh
// Copyright (c) OpenMMLab. All rights reserved.
#include <torch/types.h>

#include "pytorch_cuda_helper.hpp"

template <typename T>
__global__ void my_conv_im2col_gpu_kernel(
    const int n, const T *data_im, 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,
    const int dilation_h, const int dilation_w,
    const int batch_size,
    const int num_channels, const int height_col,
    const int width_col, T *data_col)
{
    CUDA_1D_KERNEL_LOOP(index, n)
    {
        // index index of output matrix
        const int w_col = index % width_col;
        const int h_col = (index / width_col) % height_col;
        const int b_col = (index / width_col / height_col) % batch_size;
        const int c_im = (index / width_col / height_col) / batch_size;
        const int c_col = c_im * kernel_h * kernel_w;

        const int h_in = h_col * stride_h - pad_h;
        const int w_in = w_col * stride_w - pad_w;
        T *data_col_ptr =
            data_col +
            ((c_col * batch_size + b_col) * height_col + h_col) * width_col + w_col;
        const T *data_im_ptr =
            data_im + (b_col * num_channels + c_im) * height * width;

        for (int i = 0; i < kernel_h; ++i)
        {
            for (int j = 0; j < kernel_w; ++j)
            {
                T val = static_cast<T>(0);
                const int h_im = h_in + i * dilation_h;
                const int w_im = w_in + j * dilation_w;
                if (h_im > -1 && w_im > -1 && h_im < height && w_im < width)
                {
                    val = data_im_ptr[h_im * width + w_im];
                }
                *data_col_ptr = val;
                data_col_ptr += batch_size * height_col * width_col;
            }
        }
    }
}

void my_conv_im2col_cuda(Tensor data_im,
                         const int channels, const int height,
                         const int width, const int ksize_h,
                         const int ksize_w, const int pad_h, const int pad_w,
                         const int stride_h, const int stride_w,
                         const int dilation_h, const int dilation_w,
                         const int parallel_imgs, Tensor data_col)
{
    int height_col =
        (height + 2 * pad_h - (dilation_h * (ksize_h - 1) + 1)) / stride_h + 1;
    int width_col =
        (width + 2 * pad_w - (dilation_w * (ksize_w - 1) + 1)) / stride_w + 1;
    int num_kernels = channels * height_col * width_col * parallel_imgs;

    AT_DISPATCH_FLOATING_TYPES_AND_HALF(
        data_im.scalar_type(), "my_conv_im2col_gpu", [&]
        { my_conv_im2col_gpu_kernel<scalar_t><<<GET_BLOCKS(num_kernels),
                                                THREADS_PER_BLOCK, 0,
                                                at::cuda::getCurrentCUDAStream()>>>(
              num_kernels, data_im.data_ptr<scalar_t>(),
              height, width, ksize_h, ksize_w,
              pad_h, pad_w, stride_h, stride_w, dilation_h, dilation_w,
              parallel_imgs, channels,
              height_col, width_col, data_col.data_ptr<scalar_t>()); });

    AT_CUDA_CHECK(cudaGetLastError());
}

和CPU版的类似，my_conv_im2col_cuda也是预处理了输入，并调用核函数my_conv_im2col_gpu_kernel来实现im2col。

CUDA实现和CPU几乎一样，唯一的区别就是for循环变成了CUDA_1D_KERNEL_LOOP(index, n)。这个宏是头文件里帮我们定义的，它简化了CUDA的一维循环。

编译

修改setup.py：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

from setuptools import setup
from torch.utils import cpp_extension
import os

src_root = 'panoflow/core/op'
cpp_src = ['my_conv.cpp', 'my_conv_cuda.cu']

if __name__ == '__main__':
    include_dirs = ['panoflow/core/op']
    cpp_path = [os.path.join(src_root, src) for src in cpp_src]

    setup(
        name='panoflow',
        ext_modules=[
            cpp_extension.CUDAExtension(
                'my_ops', cpp_path, include_dirs=include_dirs)
        ],
        cmdclass={'build_ext': cpp_extension.BuildExtension})

首先，要把源文件加入cpp_src里。之后，把CppExtension改成CUDAExtension。这样，就能编译新写的CUDA文件了。

写完了之后，再次python setup.py develop编译即可。

编译小技巧：不拿IDE直接写C++和CUDA源代码是很容易出错误的。但如果你想只用setup.py来验证代码的正确性，可以python setup.py develop > tmp.txt把编译输出重定向到一个文件里来查看。由于编译时的信息过多，在命令行里很难从一堆编译warning里找到最重要的error。

测试

由于Python部分在之前都已经写好了，可以直接用刚刚的单元测试文件测试了。只要把刚刚那份文件的device_name改成cuda:0即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

import torch
import torch.nn as nn
from panoflow.core.op.my_conv import MyConv2d

inc = 3
outc = 4
img_shaspe = (50, 50)

device_name = 'cuda:0'
# device_name = 'cpu'
open_bias = True


def test_one():
    ts = torch.ones([1, 1, 3, 3]).to(device_name)
    layer = nn.Conv2d(1, 1, 3, 1, 1, bias=open_bias).to(device_name)
    gt = layer(ts)
    my_layer = MyConv2d(1, 1, 3, 1, 1).to(device_name)
    my_layer.load_state_dict(layer.state_dict(), strict=False)
    res = my_layer(ts)
    res = res.to('cpu')
    gt = gt.to('cpu')
    assert torch.allclose(res, gt, 1e-3, 1e-5)


def test_two():
    ts = torch.rand([1, inc, *img_shaspe]).to(device_name)
    layer = nn.Conv2d(inc, outc, 3, 1, 1, bias=open_bias).to(device_name)
    gt = layer(ts)
    my_layer = MyConv2d(inc, outc, 3, 1, 1).to(device_name)
    my_layer.load_state_dict(layer.state_dict(), strict=False)
    res = my_layer(ts)
    res = res.to('cpu')
    gt = gt.to('cpu')
    assert torch.allclose(res, gt, 1e-3, 1e-5)


if __name__ == '__main__':
    test_one()
    test_two()

同样，没报错就说明写对了。

ResNet是CV中的经典网络。在这篇文章中，我将按照阅读论文的通用方法由粗至精地解读这篇文章。如果你对ResNet不熟，最好对着原论文阅读本文。如果你已经很熟悉ResNet了，也可以通过这篇文章查缺补漏。

粗读

摘要

摘要的前三句话开门见山地介绍本文要解决的问题以及解决的方法。

• 问题：较深的神经网络很难训练。
• 本文的工作：提出了一种能够轻松训练比以往的网络要深得多的残差学习框架。
• 本文方法的进一步解释：神经网络的层将拟合一个基于输入的残差函数，而不是一个没有参考的函数。

随后，摘要展示了这种方法取得的成就：

• 经实验而非理论证明，增加深度后，模型的训练效果和测试效果都得到了提升。
• 精度打败了当时所有的分类模型。
• 深度高达152，大幅超过了当时较火的19层的VGG，同时并没有增加多少计算量。
• 精度打败了当时所有的检测、分割等任务的模型。这证明这种方法的泛化性强。

这篇摘要没有罗列贡献，而是在提出问题后轻轻一点，介绍了本文的方法及其作用。随后，所有的篇幅都在秀这种方法的成效。看完这段话，我们可能还不知道“残差”是怎么算的，但我们知道了这种方法很厉害，测试精度、训练难易度、泛化性都十分优秀，成功解决了较深的模型难以训练的问题。

在这轮粗读中，我们可以带着以下问题去概览论文：

• 文章要解决的具体是一个怎样的问题？为什么较深的神经网络很难训练？
• 文章提出的“残差”是什么？它为什么能解决深模型难训练的问题？
• 凭什么说深模型难训练的问题被解决了？实验是怎么做的？
• 这篇文章提出的模型比其他模型好在哪？

引言

引言用十分清晰的逻辑阐明了本文的核心思想——深度残差学习框架。引言先介绍“训练更深的神经网络”这一问题的由来，再抽丝剥茧地讨论该问题的解决方法，最后自然过渡到了本文的核心思想上。看完引言，我们就应该看懂这篇文章。

深度卷积神经网络的有效性，得益于它的“深度”。那么，是不是层数更深的网络就更好呢？过去的实验显示：不是的。

更深的网络面临的第一个问题是梯度弥散/爆炸，这些问题会令网络难以收敛。但是，通过参数归一初始化和归一化层，这一问题以及得到了有效缓解。证据就是，较深的网络能够收敛。

这时，较深网络又暴露出了第二个问题——退化：增加网络的深度，反而会降低网络的精度。这种退化和过拟合还不同，因为退化不仅导致精度降低，还提升了模型的训练误差。

这种退化，是不是说明较深的网络本质上就比不过较浅的网络呢？其实并不是，退化只是因为较深的网络不是那么容易优化。考虑一个较浅的网络，我们往它的后面增加几层全等映射(y=x)。这个变深的网络与原来的网络完全等价。从这个角度看，更深的网络最少不会比更浅的网络更差才对。因此，我们要想办法让学习算法能够更好地优化更深的网络，最起码优化出一个不比较浅网络更差的网络。也就是说，我们要想办法让学习算法能够轻松学会恒等映射。

文章提出了一种深度残差学习框架。假设一层网络表示的映射是$H(x)$，则该层的非线性层要拟合另一个表示残差的映射$F(x) := H(x) - x$。换个角度看，就是一层网络的输出由原来的$F(x)$变成了$F(x)+x$，多加上了一个$x$。这样，只要令$F(x)=0$，网络就能轻松描述恒等映射，较深的网络最起码不比较浅的网络更差了。

多个数据集上的实验表明，这种方法不仅容易训练，还能得到更高的精度。在多项任务的CV竞赛中，这种方法都独占鳌头。

看完了引言，我们基本能知道文章在解决怎样的问题，也大致明白了残差学习的原理。接下来，我们来过一过这篇文章的实验，看看这种方法的效果究竟如何。

实验

我们主要看几组ImageNet上的实验。为了方便称呼，作者把不使用残差连接的普通网络叫做“平坦(plain)”网络。第一组实验对比了不同深度的平坦网络的训练误差（细线）和验证误差（粗线）。

可以看出，更深的网络有更高的训练误差和验证误差。作者排除了梯度问题的影响。这些网络在训练时使用了Batch Normalization (BN)，经调试，它们的梯度值也一切正常。因此，这种问题不是梯度消失导致的。引言中对于退化问题的描述得到了验证。

第二组实验是对残差网络ResNet的实验。我们可以把它的结果和平坦网络的放在一起对比。

从图中可以看出，使用残差学习后，更深的网络果然有了更低的训练误差和验证误差。同时，从表中可以看出，残差网络的测试误差也降低了一大截。这说明残差学习很好地解决了前文提到的退化问题，且这种有效性在测试数据上依然能保持。

由于这轮粗读我们没读方法部分，和方法有关的实验结果得跳过。可以直接翻页到和其他模型的对比表格处：

ResNet 在各个评价指标上打败了当时的其他网络，确实很厉害。

后面的表格显示，不仅是图像分类，ResNet在检测和分割等任务中也取得了第一名。

一般看到这里，一轮粗读就差不多完成了。从这轮粗读中，我们看懂了残差学习的核心思想，基本上理解了本文的核心创新点。当然，粗读深度学习相关的论文时，还可以扫一眼网络的核心模块和模型结构。精读的时候，我们再仔细理解文章的方法。

精读

残差公式

设多个层构成的模块在拟合一个映射$H(x)$，其中第一层的输入是$x$，则我们希望网络的非线性层去拟合另一个残差函数$F(x) := H(x) - x$，或者说整个模块在拟合$F(x) + x$。由于神经网络的多个非线性层能拟合任意复杂的映射，拟合一个残差函数$H(x)-x$也是可以实现的。

残差块

有了整体的思路，我们来具体看看每一个带残差的模块是怎么搭建的。具体而言，一个模块可以由输入$x$和参数集合${W_i}$计算得到（为了简化表达，bias被省略了）:

我们主要关心$F(x, {W_i})$是怎么搭建的。文章中给出了一种简单的双层残差块示例，其中，$x$是以短路的形式连接到输出上的：

注意这里的$y = W_2\sigma(W_1x)$是一个未经激活的结果。整个模块送入下一层的输出是$\sigma(y)$，还要加上激活函数。

残差学习是一个框架，每个残差块可以有更多层。比如本文在实验部分还测试了另一种三层残差块（下图是一个示例，实际上通道数可以任意改变）。

多层的残差块都是可行的。但是，单层的残差块$y = W_1x + x$和线性层几乎一样，不能提升网络的有效性。

上述的输入$x$是直接加到激活前的输出上的。这个加法操作有一个前提：模块的输入输出通道数是一样的。在输入输出通道数不同时，可以在短路连接上加一个变换维度的线性运算$W_sx$。原来直接做加法的操作叫做全等映射，这里加上一个线性运算再做加法的操作叫做投影映射。

这里的符号标记都是基于全连接层的。这些计算对卷积层来说也一样，把矩阵乘法替换成卷积即可。

网络架构

本文先借鉴VGG的思路，先搭建了一个图像大小逐渐减半，深度逐渐翻倍的平坦网络。之后，在平坦网络连续的3x3卷积层上添加残差连接。下图中的实线代表普通的短路连接，虚线表示需要变换张量形状的短路连接。

在虚线表示的残差连接中，图像的大小和深度都发生了改变。对于深度的增加，即可以使用上一节提到的投影映射，也可以直接往多出来的通道数里填0（全等映射）。对于大小的减半，无论是用怎样的深度变化方法，都可以使用步幅为2的操作来实现大小减半。

在投影时，要进行卷积操作，卷积的步幅为2很好理解。但是，文章没有详细介绍步幅为2的全等映射是怎么做的。直观上理解，就是以步幅2跳着去输入张量里取值。

文章还提出了层数不同的ResNet架构。对于较深的网络，文章使用了3层残差块。

训练细节

训练的大部分配置都借鉴自AlexNet。如果你是刚入门图像识别且以后要从事这方面的研究，可以多关注这些细节。

数据处理：

• 随机缩放图像短边至 [256, 480] 中的均匀随机值
• 随机翻转
• 裁剪成 224x224
• 像素归一化（同AlexNet）
• 标准颜色增强（同AlexNet）

归一化

• 激活前使用BN
• 参数初始化使用 He Initialization

优化配置

• batch size 256 的 mini-batch 梯度下降
• 学习率初始化为0.1，发现误差不动了就除以10
• 迭代$60 \times 10^4$轮
• weight decay=0.0001, momentum=0.9
• 没有dropout

此外，为了提高比赛中的精度，在测试时使用了10-crop（对图像裁剪10次，把输入分别输入网络，取结果的平均值）。同时，图像以不同的尺寸输入进了网络，结果取所有运算的平均值。

实验

读懂了方法，我们再来详细读一遍实验部分。

再看一次平坦网络和残差网络的对比。

在这份对比实验中，残差网络相对平坦网络没有添加任何参数。当通道数改变时，残差网络用的是0填充的全等映射。这样，平坦网络和残差网络的唯一区别就是是否有短路连接，整个对比实验非常公平。实验的结果证明了残差连接的有效性。

这篇工作还做了另一个比较平坦网络的实验：在CIFAR-10数据集上统计了平坦网络和残差网络的3x3卷积输出的标准差。标准差大，能说明输出的数值较大。

从这份对比结果中，我们能看出残差网络的输出标准差小于平坦网络。这符合残差学习的设计初衷：残差块至少是一个不会使性能变差的全等映射，其效果能够在全等映射的基础上优化。也因此，残差网络的输出大小会比平坦网络更靠近0。

看完了和平坦网络的对比，再看一看不同配置的ResNet直接的对比。文章先比较了全等映射和投影映射的短路连接。本文探讨了短路连接的3种配置：A) 全部使用全等连接 B) 有通道数变动的地方才使用投影连接 C) 全部使用投影连接。它们的表现如下：

可以看出，投影用得越多，效果越好。作者认为，B相对A更好，是因为通道数变化时0填充的部分没有残差学习（也就是没有做$+x$的操作）；C相对B更好，是因为参数更多了。这些实验结果说明，投影连接并不能解决退化问题，出于性能的考虑，本文没有在其他实验中使用C。

后来，B是公认的ResNet标配。

文章还讨论了构筑更深网络的bottleneck结构。如前文所述，50层以上的ResNet在每个残差块里使用了3个1x1, 3x3, 1x3的卷积。这种设计主要是为了缩短训练时间。ResNet-50和ResNet-34有着同样的时间复杂度，却因为深度更大，精度更高。

这篇论文的主要实验就是这些。论文后面还展示了一个比较有趣的实验：在CIFAR-10数据集上训练超过1000层的ResNet。实验结果显示，1000多层的ResNet仍能成功被优化，训练误差也降到了很低，但是测试误差没有110层的网络好。作者认为，这是因为训练数据太少，网络过拟合了。

总结

ResNet是基于深度学习的计算机视觉中里程碑式的工作。通过阅读这篇论文，我们发现，ResNet使用的残差结构本身并不十分复杂。这篇工作真正出彩之处，是发现了深度神经网络的一个普遍存在的问题，并用精彩的推理设计出了一套能够解决此问题的方案。这种方案确实很有效，基于残差学习的ResNet击败了同时代的所有网络。残差连接被用在了后续几乎所有网络架构中，使用ResNet为backbone也成了其他CV任务较为常用的初始配置。

在这堂课里，我们要学习两个具体的应用：人脸识别、风格迁移。

相信大家已经在很多地方见识过人脸识别应用了：在火车站，要通过身份证和人脸核实身份；办理业务时，可以用手机完成人脸识别以核实身份；进办公楼时，员工只要刷脸通过就可以打开门禁。通过这节课的学习，我们能够学会如何用CNN完成人脸识别。

神经网络风格迁移是一项有趣的应用。它利用了CNN捕捉图像抽象信息的特点，能够把一幅图像的风格转移到另一幅图像上，从而生成一幅新的图像。这项技术不需要从头训练网络，学完这门课后，我们能快速地用代码实现神经网络风格迁移。

课堂笔记

人脸识别

准确来说，在人脸识别(Face Recognition)任务中，会给定一个有$K$个人的数据库。之后，每一次识别都会输入一张图片，输出这张图片是$K$个人中的哪一个，或者没有检测到相关人士。

有一个与这个相关的任务叫做人脸验证(Face Verification)。这个任务稍微简单一些，输入是一张图片和一个标记身份的数据（比如身份证号），要求输出图片中的人是否符合该身份。

单样本学习

按我们之前学的方法，假如我们在$K$个人的数据库上做识别（分类）任务，应该套用一个CNN模型，并在模型最后接一个$K+1$类的softmax，表示输入图片是K个人中的哪一个，或者都不是。

但是，这样的架构不适合人脸识别任务。以公司的门禁识别为例，这种方法有如下的缺点：

1. 每来一个新同事，模型就要重新训练一次。
2. 每个人都得上传大量的个人照片供网络训练。

理想情况下，我们希望模型能从一张人脸照片中学会分辨这张人脸，这样每个新同事只需要上传一张照片即可。这叫做单样本学习(One-shot Learning)。

为了完成单样本学习，我们可以从另一个角度来建模人脸识别问题：如果输入的人脸和数据库里某张人脸极为相似，我们就说识别出了这张人脸；否则，就说没有识别到有效的人脸。

这样，人脸识别问题就被转换为了一个求两张图片相似度的问题。我们可以让网络学习一个输入是两张图片，输出是二者相似度的一个映射。

孪生网络

在完成和相似度有关的问题时，一种常见的做法是使用孪生网络（Siamese Network)。

假设网络的倒数第二层有128个神经元。在普通分类网络中，这128个神经元输出的长度为128的向量会被输入进最后的softmax层。而在孪生网络中，我们要去掉softmax层，并用这个没有sofrmax的网络$f$分别输出两张图片$x^{(1)}, x^{(2)}$对应的128维向量$f(x^{(1)}), f(x^{(2)})$。这样，每张图片有了唯一对应的一个128维向量，这个向量可以看成是该图片的编码(encoding)。而我们又知道，对向量求相似度是很方便的。我们可以利用两张图片的编码求出相似度。

说起向量的相似度，最容易想到是向量间的距离$||v - u||^2$。因此，我们可以i设法让网络学会这样一种关系：

• 若$x^{(i)}, x^{(j)}$是同一人，则$||f(x^{(i)}) - f(x^{(j)})||^2$很小。
• 若$x^{(i)}, x^{(j)}$不是同一人，则$||f(x^{(i)}) - f(x^{(j)})||^2$很大。、

为了达成这个目标，我们来看看应该用怎样的误差来指导网络的优化方向。

三元组误差

在让网络学习基于距离的相似度时，一种常用的误差是三元组误差(Triplet Loss)。

在一轮训练中，我们要用到3张图片：一张基准图(anchor)$A$和与其相对的正负样本$P, N$各一张。设$d(A, B)$为图片$A, B$的编码的距离，则我们希望$d(A, P) \leq d(A, N)$。

左边的人是吴恩达的妻子，这几张图片已经出现了好几次。

移个项，即我们希望：

但这个条件太容易满足了，令$f(x)=0$恒成立的话无论怎样的输入都可以使上式左侧为0了。因此，我们要加一点小小的约束：

这个$\alpha$是一个较小的数，比如可以令$\alpha=0.2$。为了达成这个目标，我们可以设置以下的误差函数

这里取max的直观解释是：只要满足开始那个不等式，让正样本和负样本能够分清楚就行了。至于两类样本要不要分辨得那么分明，我们并不关心。

为了利用这个误差训练网络，我们要在训练集里为同一个人准备多张照片。比如1000个人，每人100张照片，共100000张照片构成训练集。

在提出此设计的FaceNet中，还有一些小细节：为了加大训练难度，让模型效果更好，每次训练时应该选择较难的三元组$A, P, N$。详情请阅读原论文。

人脸验证与二分类问题

其实，判断两张图片的编码是否相似的问题可以简单地转化成一个二分类问题：把两张图片的编码“拼起来”，输入进一个逻辑回归的单元，让网络输入两张图片是否相似。

这里把两张图片“拼起来”有很多的实现方式。一种简单的方式是求两个编码的绝对值（L1距离）。

另外，在部署时，由于比较的图像的编码是可以预处理的，只需要用神经网络跑一遍输入图片即可。

神经网络风格迁移

我之前写了一篇详细介绍神经网络风格迁移的文章，我认为那篇文章比这堂课更好理解，非常推荐大家阅读。因此，这部分笔记我会写得潦草一些。

什么是神经网络风格迁移？

如上图所示，在神经网络风格迁移中，输入一张表示内容的图(C)和一张表示画家风格的图(S)，我们可以借助神经网络生成一幅融合内容与风格的图片(G)。

接下来实现神经网络风格迁移时，我们会关注CNN浅层和深层提取出来的特征。因此，在具体介绍风格迁移的实现之前，我们先来看看CNN各层网络究竟学到了什么。

深度卷积网络学到了什么？

为了设法可视化神经网络每一层的输出，我们可以考虑把整个训练集喂入网络，找出使某一层某一神经元激活输出最大的几个像素块。

以AlexNet为例，令第一层某几个神经元激活输出最大的像素块是：

可以看出，浅层的神经网络更关注不同方向轮廓信息。

用同样的方式可视化更深的层，可以看出网络关注的信息越来越抽象。

损失函数

和优化一个网络的参数不同，在神经网络风格迁移中，我们要把一幅图像当成可以优化的参数，通过修改图像的像素值来最小化一个损失函数。让我们看看这个损失函数是怎么定义的。

损失函数由两部分构成：生成图像(G)和内容图像(C)的内容损失与和风格图像(S)的风格损失。二者之间的比例靠比例系数$\alpha, \beta$控制。

我们来看一个优化这个损失函数的步骤示例：

1. 生成一个随机图像$G: 100 \times 100 \times 3$
2. 使用梯度下降更新$G$，最小化$J(G)$

内容损失

内容损失的计算方式如下：

1. 使用一个预训练的网络（最常用的是VGG）。
2. 选择神经网络中的某一层$l$，这一层相关的数据会用来计算内容损失。$l$越浅，表示越具体的图像；$l$越深，表示越抽象的图像。一般选取适中的$l$。
3. 令$a^{l}, a^{[l][G]}$为图像$C, G$网络第$l$层的激活输出。我们认为，如果这两个值很相似，则两幅图像的内容就很相似。
4. 令$J_{content}(C, G) = ||a^{l} - a^{[l][G]}||^2$

风格损失

看完了内容损失，现在来看风格损失的计算方式。我们刚刚一直在用名词“风格”。这个词放在美术里，可以指画家的绘画风格。而对于神经网络来说，“风格”又是什么意思呢？

和刚刚的内容损失类似，计算风格损失时，我们也要选择CNN的某层$l$，计算这一层卷积激活输出的通道相关量。

这里“通道相关量”反映的是图像各个通道间两两的相关关系。这句话可能有点拗口，让我们看一个详细的通道相关量示例。

假设一个激活输出有5个通道，我们用颜色“红-黄-紫-绿-蓝”称呼它们。每个通道表示的是一类特征，比如红色的通道表示图像有没有竖着的条纹，黄色通道表示图像有没有橙色的色块。我们想计算红黄通道之间的相关量，其实就是计算图像每个像素处红色通道的值和绿色通道的值之间的相关关系，看看它们是会同时出现，还是一个出现了另一个就不出现。（两个数值的相关量就是它们的乘积，具体的数学表示会在后文中给出）

为什么这样的相关关系能够捕捉到风格信息呢？红色通道和黄色通道的相关关系，就是是不是有竖条纹的地方就有橙色色块。这样一种线条和颜色的关联，就可以代表图片的风格。（下图中红色的框和黄色的框分别圈出了两种特征）

从直觉上认识了风格，现在我们来看一看风格的具体计算方法。在计算两个向量的相关量时，可以直接求两个向量的积。因此，我们要用到一种叫做”风格矩阵“的中间量，它计算了所有通道向量的乘积：

其中，$a^{[l]}_{i, j, k}$是第$i$行$j$列第$k$个通道的激活输出，风格矩阵$G$的形状是$n_c^{[l]} \times n_c^{[l]}$，$G^{[l]}_{kk’}$描述的是第$l$层激活输出中，第$k, k’$个通道间的相关量。

在数学中这个矩阵叫做gram矩阵。

有了风格矩阵，就可以算风格损失了。为了简化表示，我们用$G^{[l]}(S), G^{[l]}(G)$分别表示风格图像S和生成图像G的风格矩阵。和刚刚一样，我们选择一层$l$，并计算风格误差：

这是一层上的风格误差。其实，我们还可以指定多个层上的风格误差，以使生成图像能够拟合风格图像在多个卷积层上（抽象程度不同）的风格。多个层的风格误差就是各层风格误差之和。

1D和3D上的卷积

在结束这门课之前，我们来学习最后一项内容——1D和3D数据上的卷积。之前，我们只关注2D的图像数据，当维度不是2时，卷积有怎样的变化呢？

1D数据可以是心电图。和可以用2D卷积捕捉2D图像的边缘一样，我们可以用下图中那个尖尖的1D卷积来捕获高频的数据。1D卷积前后的形状变化规律和2D一样，比如用16个长度为5的卷积卷一个$14 \times 1$的1D图像，可以得到$10 \times 16$的1D图像（第一维是图像长度，第二维是通道数）。

3D数据也是类似的道理。用16个$5 \times 5 \times 5$的卷积核卷一个单通道$14 \times 14 \times 14$的单通道图像，可以得到一个$5 \times 5 \times 5 \times 16$的图像。

总结

这节课主要介绍了人脸识别和神经网络风格迁移两项任务，最后顺便介绍了1D和3D上的卷积。

• 人脸识别
  • 人脸验证与人脸识别任务的定义
  • 如何对单样本学习建模
  • 孪生网络
  • 基于三元组误差和二分类误差的人脸识别网络
• 神经网络风格迁移
  • 神经网络风格迁移的输入输出
  • 利用小技巧可视化 CNN 学到的图像信息
  • 生成风格迁移图像时的内容误差与风格误差

人脸识别任务依赖于数据集，项目搭建起来比较麻烦。对于这一周的内容，我就不展示其他的代码实战笔记了，欢迎阅读我之前写的神经网络风格迁移的实现。

在这篇文章中，我将给出一份带运行示例的NMS Python脚本，并对算法和代码进行详细解说。相信大家看完这篇文章后，能够轻松地掌握NMS的底层原理。

如果你对目标检测的基础知识不太熟，欢迎先阅读我的上篇文章：目标检测简介

示例脚本（包括可视化的代码）链接：https://github.com/SingleZombie/DL-Demos/tree/master/dldemos/nms

算法介绍

在目标检测算法中，为了尽量不漏掉物体，会输出大量的检测结果（每一条结果由检测概率与检测框组成）。这些检测结果很可能有重复的，即会有多个框标出了同一个物体。下图就是一个例子，算法正确识别出了两辆车，但却有多个检测框标出了同一辆车。

我们需要一个算法来过滤多余的检测框。最常用的算法就是NMS(Non-Maximum Suppresion, 非极大值抑制)。该算法的思路很简单：只保留局部概率最大的检测框，与其重合的其他检测框都会被舍去。

算法的伪代码如下：

1
2
3
4
5
6
7

输入所有检测概率、检测框
当还有检测框没有处理：
    从剩下的框里挑出检测概率最大的框 bbox_a
    遍历所有没有处理的框bbox_i：
        如果 bbox_i != bbox_a 且 bbox_i 与 bbox_a 重合：
            舍去 bbox_i
    把 bbox_a 输出成一个检测结果

当然，这个算法的描述还不够准确：究竟该怎么定义两个检测框是“重合”呢？如果两个检测框有交集就说它们重合是行不通的，因为图片中很可能会有挨得很近的物体，它们的检测框就是相交的。因此，我们一般用IoU（交并比）来描述两个框的重合程度，如果IoU超出某个阈值，就说这两个框是“重合”的。

IoU的计算很简单，算出两个矩形的「交面积」，再用两个矩形面积之和减去「交面积」就可以得到「并面积」，「交面积」比「并面积」就是IoU。

在NMS之前，一般还会先做一步预处理：对于预测概率小于某个阈值的检测框，我们不信任它们，会在进行NMS之前就把它们舍去。在代码实现时这部分逻辑常常会放到NMS的函数里。这样，整个NMS的流程是这样的：

上述的NMS针对是识别一种物体的任务，在推广到多分类NMS时，只要把输入的检测概率换成有物体的概率乘上概率最大的类别的概率即可。

代码实现

理论上这段代码是兼容任何格式的张量的。经测试，NumPy, PyTorch 的张量都可以正确运行。

先看一下IoU的实现代码：

1
2
3
4
5

def iou(b1: Tuple[int, int, int, int], b2: Tuple[int, int, int, int]) -> float:
    intersection = box_intersection(b1, b2)
    inter_area = area(intersection)
    union_area = area(b1) + area(b2) - inter_area
    return inter_area / union_area

所有的检测框用一个长度为4的元组表示。box_intersection用于计算两个框的相交框，area用于计算框的面积。这段代码和之前描述得一样，先获取相交区域，计算交面积，再计算并面积，最后算除法。

box_intersection的实现如下：

1
2
3
4
5
6
7
8
9
10
11

def box_intersection(
        b1: Tuple[int, int, int, int],
        b2: Tuple[int, int, int, int]) -> Tuple[int, int, int, int]:
    x11, y11, x12, y12 = b1
    x21, y21, x22, y22 = b2

    xl = max(x11, x21)
    xr = min(x12, x22)
    yt = max(y11, y21)
    yb = min(y12, y22)
    return (xl, yt, xr, yb)

算相交区域，可以理解成分别求出x, y方向上的相交部分，再把两部分合成一个框。而求直线的相交，就是取最大的左端点和最小的右端点。

area的实现如下：

1
2
3
4
5

def area(box: Tuple[int, int, int, int]) -> float:
    x1, y1, x2, y2 = box
    width = max(x2 - x1, 0)
    height = max(y2 - y1, 0)
    return width * height

如果两个框不相交，x2 - x1和y2 - y1会出现小于0的情况。因此，要保证它们最小为0，再算面积就不会有错了。

有了iou，就可以实现了NMS了。我的NMS函数定义如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

def nms(predicts: np.ndarray,
        score_thresh: float = 0.6,
        iou_thresh: float = 0.3):
    """Non-Maximum Suppression

    Args:
        predicts (np.ndarray): Tensor of shape [n, 5]. The second demesion
            includes 1 probability and 4 numbers x, y, w, h denoting a bounding
            box.
        score_thresh (float): The boxes with probability lower than
            score_threash will be discarded.
        iou_thresh (float): The threshold determining whether two boxes are
            "overlapped".

    Returns:
        (np.ndarray, List[int]): The filtered predictions and the indices of
            remaining boxes.
    """

NMS算法需要输入检测结果predicts，输出过滤后的检测结果filtered_predicts。此外，NMS算法有两个输入属性score_thresh， iou_thresh，分别表示被选定的检测框最小需要的概率、判断两个框是否重合的IoU阈值。为了方便其他的计算（比如多分类NMS），我还输出了一个索引数组indices，表示被选中的框的索引。这方面的逻辑可以根据自己的项目要求进行优化，没有统一的规定。

NMS的实现和开始的伪代码几乎一模一样：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

def nms(predicts: np.ndarray,
        score_thresh: float = 0.6,
        iou_thresh: float = 0.3):
    n_remainder = len(predicts)
    vis = [False] * n_remainder

    # Filter predicts with low probability
    for i, predict in enumerate(predicts):
        if predict[0] < score_thresh:
            vis[i] = True
            n_remainder -= 1

    # NMS
    output_predicts = []
    output_indices = []
    while n_remainder > 0:
        max_pro = -1
        max_index = 0
        # Find argmax
        for i, p in enumerate(predicts):
            if not vis[i]:
                if max_pro < p[0]:
                    max_index = i
                    max_pro = p[0]

        # Append output
        max_p = predicts[max_index]
        output_predicts.append(max_p)
        output_indices.append(max_index)

        # Suppress
        for i, p in enumerate(predicts):
            if not vis[i] and i != max_index:
                if iou(p[1:5], max_p[1:5]) > iou_thresh:
                    vis[i] = True
                    n_remainder -= 1
        vis[max_index] = True
        n_remainder -= 1

    return output_predicts, output_indices

一开始，先声明一些辅助的变量。n_remainder表示还有多少个框没被访问，vis[i]表示第i个框是否被访问。

1
2

n_remainder = len(predicts)
vis = [False] * n_remainder

一开始，先过滤那些概率过小的框：

1
2
3
4

for i, predict in enumerate(predicts):
    if predict[0] < score_thresh:
        vis[i] = True
        n_remainder -= 1

之后，正式进入NMS，先准备好输出的列表：

1
2
3
4

# NMS
output_predicts = []
output_indices = []
    while n_remainder > 0:

在还有没访问的框时，先找出剩下的框中概率最大的那个框：

1
2
3
4
5
6
7
8
9

while n_remainder > 0:
    max_pro = -1
    max_index = 0
    # Find argmax
    for i, p in enumerate(predicts):
        if not vis[i]:
            if max_pro < p[0]:
                max_index = i
                max_pro = p[0]

之后，抑制掉和概率最大框“重合”的框。

1
2
3
4
5
6
7
8
9
10
11

while n_remainder > 0:
    # Find argmax
    ...

    max_p = predicts[max_index]
    # Suppress
    for i, p in enumerate(predicts):
        if not vis[i] and i != max_index:
            if iou(p[1:5], max_p[1:5]) > iou_thresh:
                vis[i] = True
                n_remainder -= 1

最后，把这个结果添加进输出，并维护好vis和n_remainder。

1
2
3
4
5
6
7
8
9
10
11
12

while n_remainder > 0:
    # Find argmax
    ...

    # Suppress
    ...

    # Append output
    output_predicts.append(max_p)
    output_indices.append(max_index)
    vis[max_index] = True
    n_remainder -= 1

这里可以把vis[max_index] = True那两行移到抑制操作的前面，这样判断里就不用加i != max_index了。我这样写是为了强调一下，判断重合的时候不需要判断自己。

实现完了，返回结果。

1

return output_predicts, output_indices

单元测试

做单元测试的时候，最好是找一份现有的实现用做对照。为了让整个测试过程更贴合实际一些，我用MMDetection的YOLOv3跑了一个NMS前和NMS后的检测框结果，并把NMS前的检测框输入进了自己实现的NMS，比较了两份NMS的输出。以下是具体的测试过程。

照着MMDetection的官方文档，下载好YOLOv3模型，用下面的代码即可对MMDetection的demo图片进行推理并保存结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14

from mmdet.apis import inference_detector, init_detector, show_result_pyplot
from mmdet.models.detectors import BaseDetector

# Choose to use a config and initialize the detector
config = 'configs/yolo/yolov3_mobilenetv2_320_300e_coco.py'
# Setup a checkpoint file to load
checkpoint = 'checkpoints/yolov3_mobilenetv2_320_300e_coco_20210719_215349-d18dff72.pth'
# initialize the detector
model: BaseDetector = init_detector(config, checkpoint, device='cuda:0')

img = 'demo/demo.jpg'
result = inference_detector(model, img)

model.show_result(img, result, score_thr=0.3, out_file='work_dirs/1.jpg')

为了得到NMS前的输入，我在mmdet/models/dense_head/YOLOV3Head.get_bboxes里面插入了一段输出结果的代码。

1
2
3
4
5
6
7
8
9
10
11

...

save_dict = {
    'bboxes': bboxes,
    'cls_scores': scores,
    'prob': objectness
}
torch.save(save_dict, 'work_dirs/bboxes.pt')

det_bboxes, det_labels = multiclass_nms(
...

得到NMS前的输入是

把这份数据输入进我们的NMS中，得到的可视化结果如下：

这跟开始那份输出一模一样。看来我们实现的这份NMS完全正确。

如果你对测试过程、可视化、多分类NMS感兴趣，欢迎直接阅读我的项目源码。

总结

在这篇文章中，我给出了一份NMS的简洁而靠近底层的Python实现，并对NMS算法进行了介绍。通过阅读这篇文章，相信大家已经完全理解了NMS的原理，并且能够用任何一种语言实现NMS。一般的NMS开源实现支持的参数更多，代码会更复杂一些，但它们的核心和我的这份实现是一样的。

这份NMS的实现还有很大的改进空间。比如每轮求概率最大的框时，可以先排好序，或者用优先队列，这样均摊下来每次获取概率最大的框的复杂度是O(logn)。但是后面判断重复的框一定有一个O(n)的计算，这部分的优化并不显著。大家有余力可以参考成熟的CV项目里NMS是怎么高效用C++实现的。

这节课中，我们要学习计算机视觉中最重要的任务之一——目标检测任务。我们会先认识目标定位和关键点检测这两个比较简单的任务，慢慢过度到目标检测任务。之后，我们会详细学习目标检测的经典算法YOLO。最后，我们会稍微认识一下语义分割任务及适用于此问题的U-Net架构。

课堂笔记

目标定位

在图像分类问题中，给定一幅图片，我们只要说出图片里的物体是什么就行了。在这堂课要讨论的任务中，我们还要多做一件事——定位。我们要先用边框圈出图中的物体，再说出框里的物体是什么。这叫做带定位(localization)的分类问题。更进一步，如果我们不再是只讨论一个物体，而是要把图片中所有物体都框出来，并标出每一个物体的类别，这就是目标检测问题，

我们对分类任务的神经网络结构已经很熟悉了。那么，带定位的分类该使用怎样的网络呢？实际上，一个边框可以用边框中心和边框宽高这四个量表示。除了softmax出来的分类结果外，我们只要让网络再多输出四个数就行了。如下图所示：

这里，要统一一下对于边框的定义。我们用$b_x, b_y$表示边框的中心坐标，$b_h, b_w$表示边框的高、宽。

来看一下标签$y$的具体写法。假设一共有四类物体：行人、汽车、摩托车、背景（没有物体）。那么，标签$y$应该用$y=[p_c, b_x, b_y, b_h, b_w, c_1, c_2, c_3]^T$表示。其中，$p_c$表示图中有没有物体。若$p_c=1$，则$c_1, c_2, c_3$分别表示物体属于除背景外的哪一类；若$p_c=0$，则其他值无意义。

这样，在算误差时，也需要分类讨论。若$p_c=1$，则算估计值与标签8个分量两两之间的均方误差；若$p_c=0$，只算$p_c$的均方误差，不用管其他量。

只要更换一下神经网络的输出格式，我们就能得到一个完成目标定位问题的网络。

关键点检测

我们刚刚学了用2个点表示一个边框。其实，拓展一下边框检测，就是一个关键点（英文有时叫做”landmark”，是“地标”的意思）检测问题。

比如，在人脸关键点检测中，我们可以定义一堆关键点，分别表示眼睛、鼻子、嘴巴……的位置。我们还是让网络先输出一个数，表示图中有没有人脸；再输出2n个数，表示n个人脸关键点。这样，网络就能学习怎么标出人脸关键点了。

很多应用都基于人脸关键点检测技术。比如我们检测到了眼睛周围的关键点后，就可以给人“戴上”墨镜。

总之，通过这一节的学习，我们要知道，目标定位中输出2个坐标只是关键点检测的一个特例。只要训练数据按照某种规律标出了关键点，不管这些关键点是表示一个框，还是人脸上各器官的位置，网络都能学习这种规律。

目标检测

有了之前的知识储备，现在我们来正式学习目标检测。目标检测可以用一种叫做“滑动窗口”的技术实现。

假设我们要构建一个汽车的目标检测系统。我们可以先构造一个汽车分类数据集——数据集的x是一些等大的图片，y表示图片里是不是有汽车。如果图片里有汽车，汽车应该占据图片的大部分位置。

通过学习，网络就能够判断一个框里的物体是不是汽车了。这样，我们可以用一个边框框出图片的一部分，裁剪下来，让网络看看图片的这一部分是不是汽车。只要我们尝试的次数够多，总能找出图中的汽车。

在遍历边框时，我们是通过“滑动”的方法：遍历边框的大小，选择好大小后把框放到左上角，先往右移，再往下移。所以这种方法叫做“滑动窗口”。

滑动窗口算法有一个缺点：如果我们移动窗口的步伐过小，则运行分类器的次数会很多；如果移动窗口的步伐过大，则算法的精度会受到影响。在深度学习时代之前，分类器都是一些简单的线性函数，能够快速算完，多遍历一些滑动窗口没有问题。而使用了深度CNN后，遍历滑动窗口的代价就很大了。

幸好，滑动窗口也可以通过卷积来生成，而不一定要遍历出来。让我们看下去。

基于卷积的滑动窗口

滑动窗口其实可以通过执行巧妙的卷积来生成。在那之前，我们先学一门前置技能：怎么把全连接层变成卷积层。

前两周学习CNN时，我们学过，卷积结束后，卷积的输出会被喂入全连接层中。实际上，我们可以用卷积来等价实现全连接层。比如下图中，一个$5 \times 5 \times 16$的体积块想变成一个长度为400的向量，可以通过执行400个$5 \times 5$的卷积来实现。

知道了这一点后，我们就可以利用卷积来快速实现滑动窗口了。

假设我们按照上一节的算法，先实现了对$14 \times 14$的小图片进行分类的分类器。之后，我们输入了一张$16 \times 16$的大图片。我们遍历滑动窗口，令步幅为2。这样，理论上，有4个合法的滑动窗口，应该执行4次分类器的运算，如下图所示：

可是，仔细一想，在执行4次分类器的过程中，有很多重复的运算。比如，对于4个滑动窗口中间那共有的$12 \times 12$个像素，它们的卷积结果被算了4次。理想情况下，只需要对它们做一次卷积就行了。这该怎么优化呢？

其实，很简单，我们可以利用卷积本身的特性来优化。卷积层只定义了卷积核，而没有规定输入图像的大小。我们可以拿出之前在$14 \times 14$的图像上训练好的卷积层，把它们用在$16 \times 16$的图片的卷积上。经过相同的网络，$16 \times 16$的图片会生成一个$2 \times 2$大小的分类结果，而不是$1 \times 1$的。这$2 \times 2$个分类结果，恰好就是那4个滑动窗口的分类结果。通过这样巧妙地利用卷积操作，我们规避了遍历滑动窗口带来的重复计算。

不过，这个方法还是有一些缺陷的。在刚才那个例子中，$16 \times 16$的图片其实可以放下9个$14 \times 14$大小的边框。但是，由于分类网络中max pool的存在，我们只能生成4个分类结果，也就是步幅为2的滑动窗口的分类结果。同时，最准确的检测框也不一定是正方形的，而可能是长方形的。为了让生成的滑动窗口更准确一些，我们要用到其他方法。

预测边框

在这一节，我们要使用YOLO(You Only Look Once)算法解决上一节中碰到的问题。还记得这周课开头学的目标定位问题吗？我们可以把滑动窗口和目标定位结合一下。

给定一幅图像，我们可以把图像分成$3 \times 3$个格子。训练模型前，我们要对训练数据做预处理。根据每个训练样本中物体的中心点所在的格子，我们把物体分配到每一个格子中。也就是说，不管一个物体的边框跨越了几个格子，它的中心点在哪，它就属于哪个格子。比如对于下图的训练样本，右边那辆车就属于橙色的格子。之后，我们给每个格子标上标签$y$。这个标签$y$就是目标定位中那个表示图片中是否有物体、物体的边框、物体的类别的标签向量。对于这个$3 \times 3$的格子，有9个标签向量，整个标签张量的形状是$3 \times 3 \times 8$。

这样，每一幅图像的输出和标签一样，也是一个$3 \times 3 \times 8$的张量了。输入一幅图片后，我们利用上一节学的卷积滑动窗口，同时预测出每个格子里的物体边框。

另外，这里要详细讨论一下$b_x, b_y, b_h, b_w$的表示方法。由于我们只关心框相对于格子的位置，因此我们可以把规定一个格子的边长为1。这样，就满足$0 \leq b_x, b_y \leq 1$了。不过，由于物体的边框可以超出小框，$b_h, b_w > 1$是很有可能的。

看到这，大家可能会有一些疑问：如果一个格子里有多个物体呢？的确，这个算法无法输出一个格子里的多个物体。一种解决方法是，我们可以把格子分得更细一点，比如$19 \times 19$个格子。这样，可以被检测到物体会多一些。但是，增加格子数又会引入一个新的问题——多个格子检测到了同一个物体。下面的两节里我们会尝试解决这个新的问题。

吴恩达老师说，YOLO这篇论文很难读懂，他和其他几个资深研究者都花了很大的功夫才读懂这篇论文。

IoU(交并比)

在目标检测中，有一个微妙的问题：框出一个物体的边框有无数个，想精确框出标签的边框是不可能的。怎么判定一个输出结果和标签里的边框“差不多”呢？这就要用到IoU(Intersection over Union，交并比) 这个概念。

IoU，顾名思义，二者的交集比上二者的并集，很好理解。比如下图中，网络的输出是紫框，真值是红框。二者的并集是绿色区域，交集是橙色区域。则IoU就是橙色比绿色。

依照惯例，如果IoU$\geq 0.5$，我们就认为网络的输出是正确的。当然，想更严格一点，0.6,0.7也是可以的。

IOU 也是 “I owe you(我欠了你的钱)”的缩写，哈哈哈。

NMS(非极大值抑制)

假设在YOLO中，我们用$19 \times 19$个小格来检测物体。可是，由于小格子太多了，算法得到了多个重复的检测框（以及每个框中有物体的概率）。这该怎么办呢？

NMS(Non-Maximum Suppresion，非极大值抑制)就是解决这个问题的算法。这个算法的名字听起来很奇怪，但大家理解了这个算法的实现后，就知道这个“抑制”是什么意思了。

讲起算法我就不困了。我会抛弃视频中的讲解思路，用我自己的逻辑讲一遍。讲算法，千万不能一上来就讲算法的步骤，一定要先讲清楚算法的思路。

在学NMS之前，我们先动动脑，看看在去掉重复的框时，我们期望得到怎样的去重输出结果。

首先，既然是去重，那么就不能出现两个框过度重合的情况。其次，我们希望留下来的框的预测概率尽可能大。

在这两个要求下，我们来看看上面那幅图的输出应该是怎样的。
我们一眼就能看出，对于左边那辆车，我们应该保留0.8的框；对于右边那辆车，我们应该保留0.9的框。

为什么我们能“一眼看出”呢？这是因为左边两个框、右边三个框恰好都分别表示了一辆车。我们能够快速地把这些框分成两类。但是，在情况比较复杂时，我们就难以快速找出最好的框了。比如下面这种情况中，两辆车很近，有些框甚至同时标出了两辆车：

为了处理这种复杂的情况，我们必须想出一种万全的算法，以筛选出那些概率比较大的框。

稍微思考一下，其实这样的算法非常简单：找出最大的框，去掉所有和它过度重合的框；在剩下的框中，找出最大的框，去掉所有和它过度重合的框；……。一直重复直到没有未处理的框为止。这就是NMS算法。

还是让我们来看看刚刚那个例子。使用NMS时，我们会先找到0.9这个框，“抑制”掉右边0.6和0.7的框。在剩下的框中，最大的是0.8这个框，它会“抑制”掉左边那个0.7的框。

接下来，让我们来严格描述一下这个算法。假设我们有$19 \times 19=361$个输出结果，每个输出结果是一个长度为5的向量$[p_c, b_x, b_y, b_h, b_w]$，分别表示有物体的概率、边框的中心和高宽（我们先不管检测多个物体的情况。事实上，当推广到多个物体时，只要往这个输出结果里多加几个概率就行了）。我们要用NMS输出应该保留的检测结果。“过度重合”，就是两个框的IoU大于0.5。

首先，先做一步初筛，扔掉概率$p_c$小于0.6的结果。

之后，对于没有遍历的框，重复执行：找出概率最大的框，把它加入输出结果；去掉所有和它IoU大于0.5的框。

这个过程用伪代码表示如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

# Input and preprocessing
input predicts of size [19, 19, 5]
resize predicts to [361, 5]

# Filter predicts with low probability
filtered_predicts = []
for predict in predicts:
    # drop p_c < 0.6
    if predict[0] >= 0.6:  
        filtered_predicts.append(predict)

# NMS
n_remainder  = len(filtered_predicts)
vis = [False] * n_remainder # False for unvisited item
output_predicts = []
while n_remainder > 0:
    max_pro = -1
    max_index = 0
    # Find argmax
    for i, p in enumerate(filtered_predicts):
        if not vis[i]:
            if max_pro < p[0]:
                max_index = i
                max_pro = p[0]

    # Append output
    max_p = filtered_predicts[max_index]
    output_predicts.append(max_p)

    # Suppress
    for i, p in enumerate(filtered_predicts):
        if not vis[i] and i != max_index:
            if get_IoU(p[1:5], max_p[1:5]) > 0.5:
                vis[i] = True
                n_remainder -= 1
    vis[max_index] = True
    n_remainder -= 1