这一章我们来分析神经网络的稳健性优化，通过三种不同的策略来求解（或近似求解）内部最大化问题。内容较多，但读下来逻辑十分连贯，收获满满。
这次我们来考虑神经网络的情况。
再一次回到那个优化问题，即内部最大化问题
\[\DeclareMathOperator*{\maximize}{maximize} \maximize_{\|\delta\| \leq \epsilon} \ell(h_\theta(x + \delta), y)\]中的\(h_\theta(x)\)此时表示的是神经网络。为了进一步研究，我们需要给出神经网络\(h_\theta(x)\)更准确的定义。用如下等式来定义一个$d$层的网络\(h_\theta(x):\mathbb{R}^n \rightarrow \mathbb{R}^k\)：
\[% <![CDATA[ \begin{split} z_1 & = x \\ z_{i+1} & = f_i(W_i z_i + b_i), \;\; i,\ldots,d \\ h_\theta(x) & = z_{d+1} \end{split} %]]>\]其中$z_i$表示第$i$层的输出；$f_i$表示第$i$层的激活函数，通常前$d-1$层使用RELU 运算\(f_i(z) = \max\{0,z\}\)，第$d$层使用线性激活函数$f_i(z) = z$。网络的参数是\(\theta = \{W_1,b_1,\ldots,W_d,b_d\}\)。损失函数采用之前使用的多分类交叉熵损失：
\[\ell(h_\theta(x), y) = \log \left ( \sum_{j=1}^k \exp(h_\theta(x)_j) \right ) - h_\theta(x)_y.\]和线性模型相比，神经网络损失函数的解空间（按文章中直译的话应该叫做损失平面 loss surface）会更加“崎岖”。换句话说就是我们上一章提到的那条“分界线”会更“崎岖”。作者随便训练了一个神经网络来说明这个问题：
import torch
import torch.nn as nn
import torch.optim as optim
torch.manual_seed(0)
model = nn.Sequential(nn.Linear(1,100), nn.ReLU(),
nn.Linear(100,100), nn.ReLU(),
nn.Linear(100,100), nn.ReLU(),
nn.Linear(100,1))
opt = optim.SGD(model.parameters(),lr=1e-2)
for _ in range(100):
loss = nn.MSELoss()(model(torch.randn(100,1)), torch.randn(100,1))
opt.zero_grad()
loss.backward()
opt.step()
plt.plot(np.arange(-3,3,0.01), model(torch.arange(-3,3,0.01)[:,None]).detach().numpy())
plt.xlabel("Input")
plt.ylabel("Output")
这会带来两个主要的挑战。
第一，在我们通常使用深度神经网络这种高维情况中，很有可能在输入空间的任何一个点都存在在一些方向上损失平面会非常陡峭，也就是说会造成损失大幅增加或减少。这正是我们在Chapter 1中见到的那样，对输入做很小的扰动就会导致损失大幅增加。换句话说，神经网络由于其损失平面的性质，特别容易出现对抗性样本。
第二，不像线性模型，想要求解内部最大化问题并不容易。正如上图所示，神经网络的损失平面是非凸的、有很多极值点（这里只考虑针对输入，而不是针对参数）。当我们尝试最大化或最小化上述函数时，给定点处的初始梯度可能会或可能不会指向全局最值的方向。第二点在我们试图生成一些对抗样本的时候并不会构成问题，毕竟正如第一点所说，有很多方向都可以导致损失大幅增加，沿着这些方向就可以找到一些还不错的对抗样本，即使它不是最优的那个对抗样本。但是，当我们想要去进行稳健性优化时，就会发现出了大问题，因为Danskin’s theorem不再成立了，进而无法去求解真正的稳健性优化问题。
求解内部最大化问题的策略
当假设函数\(h_\theta\)是神经网络时，该如何（近似地）求解内部优化问题？
\[\maximize_{\|\delta\| \leq \epsilon} \ell(h_\theta(x), y)\]有三种主要的策略来解决，分别对应着下界，精确解和上界。
具体的，我们有如下三种选择：
- 我们可以找到优化目标的下界。因为任何可行的$\delta$都会给出一个下界，所以可以“经验地来求解优化问题”，即“找到某个对抗样本”。我对此的理解是，既然算不出就强行找出来，利用像SGD这种算法，也就是Chapter 1中所做的那样。这也是目前最常用的求解内部最大化问题的策略。它背后的理由是，神经网络的局部极值问题似乎并不如我们想象的那么严重。然而，为了找到足够好的对抗样本并用其来训练稳健的模型，需要我们足够好地来求解这个问题。（下界的这个说法是针对最大化问题的，而不是里面的那个函数。而且是下界，不是下确界。每有一个对抗样本，就能刷新优化目标的下界，即至少大于等于当前值）
- 我们可以精确地求解这个优化问题。这当然很有挑战性，但其实对于许多激活函数而言，我们都可以将最大化问题形式化为一个组合优化问题，并通过混合整数规划（mixed integer programming）这样的技巧来精确求解。虽然这些方法在扩展到大型模型时会面临巨大挑战，但是对于小的问题，它们强调了一个重要的观点，即在某些情况下可以构造出内部最大化问题的精确解。
- 最后可以找到优化目标的上界，和我们在多类别线性分类中所做的类似（这部分作者还没写）。其基本策略是考虑一个松弛的网络结构，使得这个松弛版本的网络既包含原本的网络，构建方式也更易于精确地优化。不同的是，这些方法通常不会给出真实网络的对抗样本（因为它们所操作的是松弛的模型，不同于原始模型），但用它们可以证明网络所具有抵抗攻击的能力。此外，当用于稳健性优化时，这些方法的效果是最好的。
一些神经网络示例
作者训练了三个简单的神经网络来作为后续实验的测试平台。这三个神经网络分别是：两层和四层的全连接网络，以及四层卷积层加一层全连接层的卷积神经网络。
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
class Flatten(nn.Module):
def forward(self, x):
return x.view(x.shape[0], -1)
model_dnn_2 = nn.Sequential(Flatten(), nn.Linear(784,200), nn.ReLU(),
nn.Linear(200,10)).to(device)
model_dnn_4 = nn.Sequential(Flatten(), nn.Linear(784,200), nn.ReLU(),
nn.Linear(200,100), nn.ReLU(),
nn.Linear(100,100), nn.ReLU(),
nn.Linear(100,10)).to(device)
model_cnn = nn.Sequential(nn.Conv2d(1, 32, 3, padding=1), nn.ReLU(),
nn.Conv2d(32, 32, 3, padding=1, stride=2), nn.ReLU(),
nn.Conv2d(32, 64, 3, padding=1), nn.ReLU(),
nn.Conv2d(64, 64, 3, padding=1, stride=2), nn.ReLU(),
Flatten(),
nn.Linear(7*7*64, 100), nn.ReLU(),
nn.Linear(100, 10)).to(device)
训练过程的代码就不放在这里，参考 Chapter 3。
这些模型将会在后续用于测试不同的求解内部问题的策略。
内部最大化的下界
这也许是求解\(\DeclareMathOperator*{\maximize}{maximize} \maximize_{\|\delta\| \leq \epsilon} \ell(h_\theta(x + \delta), y)\)最简单直接的方法了。做法也很基本，就是反向传播。通过计算损失函数关于$\delta$的梯度，不断让梯度下降从而最大化目标函数。当然，这里我们需要保证$\delta$满足范数的约束，因此可以在每一次更新后都需要将$\delta$再重新投射到约束范围内。
快速梯度下降（Fast Gradient Sign Method）
我们首先计算出调整$\delta$的方向：
\[g := \nabla_\delta \ell(h_\theta(x + \delta),y)\]如果计算的梯度是在$\delta=0$时（也就是第一次更新），那么它实际上等价于\(\nabla_x \ell(h_\theta(x),y)\)，只是我们为了表述的一致性依然称之为关于$\delta$的梯度。
为了最大化损失，需要让$\delta$在梯度下降的方向上以一定的步长$\alpha$迈出一步（take a step）：
\[\delta := \delta + \alpha g\]注意 ：因为我们的目标是最大化损失，所以更新取梯度方向，因此这里是“+”；而最小化损失是取负梯度方向。
当然还需要保证更新后的$\delta$在约束范围内\(\|\delta\| \leq \epsilon\)。这里我们取无穷范数\(\ell_\infty\).
现在的问题是，步长应该取多大？如果我们想要让损失尽可能大的话，就应该取尽可能大的$\alpha$。如果取了一个很大的$\alpha$，不仅会被重新投射回$[-\epsilon, \epsilon]$内，而且$g$对于更新的影响就只剩下决定符号正负了（$g$的大小相较于$\alpha$变得不再重要）。于是，当我们取一个很大的$\alpha$时，更新就为
\[\delta := \epsilon \cdot \mathrm{sign}(g).\]这就是所谓的快速梯度下降法（FGSM），它是深度学习领域中最早提出的构造对抗样本的方法之一。
其代码如下：
def fgsm(model, X, y, epsilon):
""" Construct FGSM adversarial examples on the examples X"""
delta = torch.zeros_like(X, requires_grad=True)
loss = nn.CrossEntropyLoss()(model(X + delta), y)
loss.backward()
return epsilon * delta.grad.detach().sign()
我们来看看构造的效果如何。
for X,y in test_loader:
X,y = X.to(device), y.to(device)
break
def plot_images(X,y,yp,M,N):
f,ax = plt.subplots(M,N, sharex=True, sharey=True, figsize=(N,M*1.3))
for i in range(M):
for j in range(N):
ax[i][j].imshow(1-X[i*N+j][0].cpu().numpy(), cmap="gray")
title = ax[i][j].set_title("Pred: {}".format(yp[i*N+j].max(dim=0)[1]))
plt.setp(title, color=('g' if yp[i*N+j].max(dim=0)[1] == y[i*N+j] else 'r'))
ax[i][j].set_axis_off()
plt.tight_layout()
### Illustrate original predictions
yp = model_dnn_2(X)
plot_images(X, y, yp, 3, 6)
### Illustrate attacked images
delta = fgsm(model_dnn_2, X, y, 0.1)
yp = model_dnn_2(X + delta)
plot_images(X+delta, y, yp, 3, 6)
仅仅是做了一些微小的变动，预测结果就从原来的只错了一个case变成只对了一个case，而这种变动对于人眼的判断并无影响。值得指出的是，全连接网络对这个问题是非常敏感的，相比之下，卷积网络的情况能好那么一些（当然还是很敏感）。我们来看看卷积神经网络上的表现如何。
### Illustrate attacked images
delta = fgsm(model_cnn, X, y, 0.1)
yp = model_cnn(X + delta)
plot_images(X+delta, y, yp, 3, 6)
直观上来看上面的对抗样本效果很不错，但是让我们来更加严格地评估攻击方法的性能——测试误差。
def epoch_adversarial(model, loader, attack, *args):
total_loss, total_err = 0.,0.
for X,y in loader:
X,y = X.to(device), y.to(device)
delta = attack(model, X, y, *args)
yp = model(X+delta)
loss = nn.CrossEntropyLoss()(yp,y)
total_err += (yp.max(dim=1)[1] != y).sum().item()
total_loss += loss.item() * X.shape[0]
return total_err / len(loader.dataset), total_loss / len(loader.dataset)
print("2-layer DNN:", epoch_adversarial(model_dnn_2, test_loader, fgsm, 0.1)[0])
print("4-layer DNN:", epoch_adversarial(model_dnn_4, test_loader, fgsm, 0.1)[0])
print(" CNN:", epoch_adversarial(model_cnn, test_loader, fgsm, 0.1)[0])
2-layer DNN: 0.9259
4-layer DNN: 0.8827
CNN: 0.4173
在继续之前，关于FGSM有几点需要说明。
第一，FGSM是针对无穷范数\(\ell_\infty\)下的攻击所设计的，因为FGSM只是在\(\ell_\infty\)约束下的一步更新。不过也可以容易地将其推广到其它范数情况下，后续会进一步说明。
第二，FGSM基于的假设是，在$x$点处计算的梯度方向上的线性近似是函数在整个区域\(\|\delta\|_\infty \leq \epsilon\)内的一个合理的近似，但对于神经网络而言，事实并不如此，因为即使在一块很小的区域内损失平面也不是线性的。换句话说就是，在$x$这一点上的确算出了梯度下降的方向，但朝着这个方向走一步后，这个方向很大可能已经不再是最优的方向了。如果想要得到更好的对抗样本，这种仅做一步更新并投射的方法还远远不够。
投影梯度下降（Projected gradient descent）
上述讨论引出了下一种求解内部问题的方法：投影梯度下降（PGD）。即和上述过程类似，只是从一步变成了迭代、从大步长变成了小步长。基本的PGD算法就是重复如下过程：
\[% <![CDATA[ \begin{split} & \mbox{Repeat:} \\ & \quad \delta := \mathcal{P}(\delta + \alpha \nabla_\delta \ell(h_\theta(x+\delta), y)) \end{split} %]]>\]其中\(\mathcal{P}\)表示约束至范围内的投影操作（例如针对无穷范数\(\ell_\infty\)的裁剪clip）。在PGD中可以指定诸如步长、迭代次数等。实际上PyTorch的优化器就可以完成这一工作，不过作者还是手工实现了如下代码，因为想要了解到底发生了什么。
def pgd(model, X, y, epsilon, alpha, num_iter):
""" Construct FGSM adversarial examples on the examples X"""
delta = torch.zeros_like(X, requires_grad=True)
for t in range(num_iter):
loss = nn.CrossEntropyLoss()(model(X + delta), y)
loss.backward()
delta.data = (delta + X.shape[0]*alpha*delta.grad.data).clamp(-epsilon,epsilon)
delta.grad.zero_()
return delta.detach()
然后看看在卷积神经网络上的表现：
### Illustrate attacked images
delta = pgd(model_cnn, X, y, 0.1, 1e4, 1000)
yp = model_cnn(X + delta)
plot_images(X+delta, y, yp, 3, 6)
对抗样本的效果似乎不如FGSM，甚至有个数字0的case看起来并没有加任何扰动。这里我们使用了一个较大的步长$\alpha=1e4$，主要是因为在$\delta=0$时的梯度实在太小了：
delta = torch.zeros_like(X, requires_grad=True)
loss = nn.CrossEntropyLoss()(model_cnn(X + delta), y)
loss.backward()
print(delta.grad.abs().mean().item())
1.8276920172866085e-06
可以看到，在任意像素上的梯度绝对值的平均值只有$10^{-6}$。这说明在$\delta=0$周围的小区域内较为“平坦”，所以我们需要用相对较大的$\alpha$来突破这块平坦区域。
附：最速下降（steepest descent）。为了解决上述问题，这里采用（规范化的）最速下降法。作者简单地介绍了最速下降法的概念。原本的梯度下降法是通过在梯度方向上反复迭代
\[z := z - \alpha \nabla_z f(z).\]但问题是对梯度值的尺度很敏感，需要根据梯度值的尺度来调整步长。相反，规范化最速下降法则是在满足约束的步径范围内（范数约束）找到的能够最大化与梯度方向內积的方向$v$，也即求梯度方向的对偶范数。（然后在方向$v$上进行线搜索，这部分作者没提）
\[\DeclareMathOperator*{\argmax}{argmax} z := z - \argmax_{\|v\| \leq \alpha} v^T \nabla_z f(z).\]如果使用$\ell_2$范数来约束$v$，该方向可以解析地求出（不过这里作者的推导中多一个比例因子$\alpha$，所以和《Convex Optimization》中的推导不太一样）
\[\argmax_{\|v\|_2 \leq \alpha} v^T \nabla_z f(z) = \alpha \frac{\nabla_z f(z)}{\|\nabla_z f(z)\|_2}\]可以看出该方向是梯度方向单位化乘以比例因子（如果没有比例因子的话，该方向就是梯度方向，与书中结论一致）。同理，如果我们使用$\ell_\infty$来约束$v$，那么该方向为
\[\argmax_{\|v\|_\infty \leq \alpha} v^T \nabla_z f(z) = \alpha \cdot \mathrm{sign}(\nabla_z f(z)).\]回到之前的神经网络对抗攻击上：迭代选择最速下降方向更新，相当于进行了多次“mini-FGSM”。这种方法被广泛称为“投影梯度下降”（投影体现在更新后保证扰动在约束范围内）
注意：上述公式中$\delta$的更新取减号，因为讨论的是“下降”，即最小化损失；而代码里$\delta$的更新依然是加号，因为我们的目标的最大化损失。
def pgd_linf(model, X, y, epsilon, alpha, num_iter):
""" Construct FGSM adversarial examples on the examples X"""
delta = torch.zeros_like(X, requires_grad=True)
for t in range(num_iter):
loss = nn.CrossEntropyLoss()(model(X + delta), y)
loss.backward()
delta.data = (delta + alpha*delta.grad.detach().sign()).clamp(-epsilon,epsilon)
delta.grad.zero_()
return delta.detach()
### Illustrate attacked images
delta = pgd_linf(model_cnn, X, y, epsilon=0.1, alpha=1e-2, num_iter=40)
yp = model_cnn(X + delta)
plot_images(X+delta, y, yp, 3, 6)
效果好多了。不仅效果比FGSM要好，而且关于步长的选择也比之前要容易。这是因为每次更新的步长$\alpha$（不再受梯度$g$模的影响）可以和扰动的界$\epsilon$放在同一尺度上比较，我们可以选择$\alpha$为$\epsilon$的一个合理的一小部分（即可以控制更新的幅度），同时也可以选择迭代次数为$\epsilon/\alpha$的某个较小倍数（即控制迭代轮数）。然后来看看这类攻击在整个测试集上的效果。
print("2-layer DNN:", epoch_adversarial(model_dnn_2, test_loader, pgd_linf, 0.1, 1e-2, 40)[0])
print("4-layer DNN:", epoch_adversarial(model_dnn_4, test_loader, pgd_linf, 0.1, 1e-2, 40)[0])
print("CNN:", epoch_adversarial(model_cnn, test_loader, pgd_linf, 0.1, 1e-2, 40)[0])
2-layer DNN: 0.9637
4-layer DNN: 0.9838
CNN: 0.7432
相比于FGSM，有了合理的效果提升。
还有最后一种策略可以进一步微小提升效果：随机化（randomization）。这类技术实际中并不常用，因为代价比较大，但它强调了一个重点。
PGD的效果其实依然受目标函数局部最优点的限制。而这个问题可以通过多次随机选择不同的起点来缓解。这种微小的效果提升说明了从0开始计算PGD的确会陷入到一些局部最优点中。当然，不利的一面就是会成倍地增加程序的运行时间，在许多情况中都不太可行。
def pgd_linf_rand(model, X, y, epsilon, alpha, num_iter, restarts):
""" Construct PGD adversarial examples on the samples X, with random restarts"""
max_loss = torch.zeros(y.shape[0]).to(y.device)
max_delta = torch.zeros_like(X)
for i in range(restarts):
delta = torch.rand_like(X, requires_grad=True)
delta.data = delta.data * 2 * epsilon - epsilon
for t in range(num_iter):
loss = nn.CrossEntropyLoss()(model(X + delta), y)
loss.backward()
delta.data = (delta + alpha*delta.grad.detach().sign()).clamp(-epsilon,epsilon)
delta.grad.zero_()
all_loss = nn.CrossEntropyLoss(reduction='none')(model(X+delta),y)
max_delta[all_loss >= max_loss] = delta.detach()[all_loss >= max_loss]
max_loss = torch.max(max_loss, all_loss)
return max_delta
print("CNN:", epoch_adversarial(model_cnn, test_loader, pgd_linf_rand, 0.1, 1e-2, 40, 10)[0])
CNN: 0.7648
针对性攻击
上面的演示都是非针对性攻击，即只要把让模型把标签预测错误即可。但正如在第一章介绍中所述，我们也可以让模型错误地预测为我们期望的类别。
\[\maximize_{\|\delta\| \leq \epsilon} \left ( \ell(h_\theta(x + \delta), y) - \ell(h_\theta(x + \delta), y_{\mathrm{targ}}) \right ) \equiv \maximize_{\|\delta\| \leq \epsilon} \left ( h_\theta(x + \delta)_{y_{\mathrm{targ}}} - h_\theta(x + \delta)_y \right )\]我们用PGD来进行攻击，这里为了使得目标类别可以是MNIST数据集中的大多数类别，使用了较大的可允许的扰动范围，$\epsilon=0.2$。
def pgd_linf_targ(model, X, y, epsilon, alpha, num_iter, y_targ):
""" Construct targeted adversarial examples on the examples X"""
delta = torch.zeros_like(X, requires_grad=True)
for t in range(num_iter):
yp = model(X + delta)
loss = (yp[:,y_targ] - yp.gather(1,y[:,None])[:,0]).sum()
loss.backward()
delta.data = (delta + alpha*delta.grad.detach().sign()).clamp(-epsilon,epsilon)
delta.grad.zero_()
return delta.detach()
我们试图让模型预测为数字2：
delta = pgd_linf_targ(model_cnn, X, y, epsilon=0.2, alpha=1e-2, num_iter=40, y_targ=2)
yp = model_cnn(X + delta)
plot_images(X+delta, y, yp, 3, 6)
结果相当好：只用了稍微大一些的$\epsilon$就能欺骗模型将所有被测样本预测为数字2（数字2本身并没有被预测错，是因为在这个case上的损失是0）。我们在将目标类别改为数字0试试：
delta = pgd_linf_targ(model_cnn, X, y, epsilon=0.2, alpha=1e-2, num_iter=40, y_targ=0)
yp = model_cnn(X + delta)
plot_images(X+delta, y, yp, 3, 6)
这一次尽管我们使模型在所有的非0数字上都预测错了，但却并没有让所有样本都被预测成目标类别。这是因为我们的最大化的目标是类别0对应的位（logit，理解成预测值即可）减去真实类别对应的位（这个量越大会使预测朝着数字0的预测值越大、真实类别的预测值越小的方向上优化）。但在这个过程中，我们并没有关注其他类别的位发生了什么变化。也许在某些样本上，其它类别的预测值比数字0的预测值还要大。我们可以通过修改目标函数，最大化目标类别预测值的同时，最小化其它类别的预测值：
\[\maximize_{\|\delta\| \leq \epsilon} \left ( h_\theta(x + \delta)_{y_{\mathrm{targ}}} - \sum_{y' \neq y_{\mathrm{targ}}} h_\theta(x + \delta)_{y'} \right )\]
def pgd_linf_targ2(model, X, y, epsilon, alpha, num_iter, y_targ):
""" Construct targeted adversarial examples on the examples X"""
delta = torch.zeros_like(X, requires_grad=True)
for t in range(num_iter):
yp = model(X + delta)
loss = 2*yp[:,y_targ].sum() - yp.sum()
loss.backward()
delta.data = (delta + alpha*delta.grad.detach().sign()).clamp(-epsilon,epsilon)
delta.grad.zero_()
return delta.detach()
delta = pgd_linf_targ(model_cnn, X, y, epsilon=0.2, alpha=1e-2, num_iter=40, y_targ=0)
yp = model_cnn(X + delta)
plot_images(X+delta, y, yp, 3, 6)
由于这个优化目标要困难的多，因此我们没有做到完全欺骗分类器。但在骗到分类器的样本上，大多数都被预测成为目标类别。
非$\ell_\infty$范数
截至目前，我们都在考虑对$\delta$的无穷范数$\ell_\infty$约束（在讨论$\ell_\infty$约束事时不用特地将$x+\delta$约束在$[0,1]$范围内，这是因为$\ell_\infty$将$\delta$约束在了$[-\epsilon,\epsilon]$内，而$\epsilon$取值较小，即使像PGD那样要迭代多轮也不容易令$x+\delta$超出$[0,1]$的范围）。在这一小节，我们终于可以考虑$\ell_2$扰动。但由于$\ell_2$扰动的约束范围$\epsilon$要比$\ell_\infty$的约束范围稍大一些（因为相同$\epsilon$下的$\ell_2$范数球的体积要比$\ell_\infty$范数球小，为了保证相似的扰动范围，需要给$\ell_2$稍大一些的$\epsilon$），因此需要额外地将$x+\delta$约束在$[0,1]$范围内。我们使用投影后的最速下降法，有如下形式：
\[\delta := \mathcal{P}_\epsilon \left(\delta - \alpha \frac{\nabla_\delta \ell(h_\theta(x+\delta),y)}{\|\nabla_\delta \ell(h_\theta(x+\delta),y)\|_2}\right )\]其中$\mathcal{P}_\epsilon$表示投影到半径为$\epsilon$的$\ell_2$球上的操作：
\[\mathcal{P}_\epsilon(z) = \epsilon\frac{z}{\max\{\epsilon, \|z\|_2\}}.\]于是，最终的攻击如下所示：
def norms(Z):
"""Compute norms over all but the first dimension"""
return Z.view(Z.shape[0], -1).norm(dim=1)[:,None,None,None]
def pgd_l2(model, X, y, epsilon, alpha, num_iter):
delta = torch.zeros_like(X, requires_grad=True)
for t in range(num_iter):
loss = nn.CrossEntropyLoss()(model(X + delta), y)
loss.backward()
delta.data += alpha*delta.grad.detach() / norms(delta.grad.detach())
delta.data = torch.min(torch.max(delta.detach(), -X), 1-X) # clip X+delta to [0,1]
delta.data *= epsilon / norms(delta.detach()).clamp(min=epsilon)
delta.grad.zero_()
return delta.detach()
delta = pgd_l2(model_cnn, X, y, epsilon=2, alpha=0.1, num_iter=40)
yp = model_cnn(X + delta)
plot_images(X+delta, y, yp, 3, 6)
print("2-layer DNN:", epoch_adversarial(model_dnn_2, test_loader, pgd_l2, 2, 0.1, 40)[0])
print("4-layer DNN:", epoch_adversarial(model_dnn_4, test_loader, pgd_l2, 2, 0.1, 40)[0])
print("CNN:", epoch_adversarial(model_cnn, test_loader, pgd_l2, 2, 0.1, 40)[0])
2-layer DNN: 0.9232
4-layer DNN: 0.9435
CNN: 0.7854
需要注意的是，这里$\epsilon$的取值比在$\ell_\infty$的情况下要大，因为$\ell_2$球的体积$\ell_\infty$的体积的$\sqrt{n}$倍成一定比例，其中$n$是输入的维度（确切地讲，其实还有一个额外的比例因子$\sqrt{2}/\sqrt{\pi e}$以使体积在高维中更精确地匹配）。所以对于$\ell_2$球若要有和$\epsilon=0.1$的$\ell_\infty$球相同的体积，近似等价于$\ell_2$球的半径为$\frac{\sqrt{2\cdot784}}{\sqrt{\pi e}} \cdot 0.1 \approx 1.35$，这里我们使用了稍大一点的半径$\epsilon=2$，因为从经验上讲，分类器能够处理比$\ell_\infty$扰动稍大的$\ell_2$扰动。
这个例子中的一个关键点是，$\ell_\infty$攻击导致的扰动几乎分布在整个图像上（这正是$\ell_\infty$球允许的扰动）；而$\ell_2$攻击造成的扰动发生在图像的局部，因为$\ell_2$可以“权衡”某些点处大的扰动和另一些点处小的扰动。
至于$\ell_1$攻击，可以使用类似步骤，即使用PGD算法，然后将更新步长投射到$\ell_1$范数球上（这一步要比映射到$\ell_2$和$\ell_\infty$范数球上更复杂）。$\ell_1$范数的特点是会更鼓励稀疏性，因此使用$\ell_1$范数球生成对抗样本的结果会是只有一些像素点做出了调整。
小结
作者再一次吐槽了该领域中论文关于攻击方法起的各式各样的名字。其实核心就两点（作者的观点）：1.用的哪种扰动范数球；2.用什么方法优化该范数扰动。
精确求解内部最大化问题（组合优化）
在上面的内容里我们考虑了在实际操作中生成对抗样本的“标准”方法，现在我们来考虑另一种思路，即使用组合优化来精确求解内部最大化问题。和线性模型的多分类情况类似，实际上我们并没有精确求解内部最大化问题，但通过考虑针对每一种可能类别的针对性攻击，我们可以准确确定在某个半径范围内是否存在一个对抗样本。我们的策略就是使用混合整数线性规划（mixed linear programming，一种组合优化的方法）。尽管我们不期望这些方法能够扩展到真正大尺度的问题上，但对于小的模型来说，这些方法的确提供了对神经网络的精确推理。
针对性攻击的约束形式
首先我们考虑对基于RELU的前馈神经网络的针对性攻击。
\[% <![CDATA[ \begin{split} z_1 & = x \\ z_{i+1} & = f_i(W_i z_i + b_i), \;\; i,\ldots,d \\ h_\theta(x) & = z_{d+1} \end{split} %]]>\]其中对于\(i=1,\cdots,d-1\)，\(f_{i}(z)=\operatorname{ReLU}(z)\)；\(f_{d}(x)=x\)，即最后一层是类别标签向量。
我们首先考虑如何将对抗攻击描述成优化问题。对于针对性攻击，试图去最小化真实类别所对应的标签位\(h_\theta(x+\delta)_y\)，同时最大化目标类别的标签位\(h_\theta(x+\delta)_{y_{\mathrm{targ}}}\)。我们\将把优化问题表述为关于输入$x$和所有中间激活变量$z_i$上的一个优化问题，但是有约束条件来约束这些激活值以遵循神经网络的行为。我们重新回到$\ell_\infty$攻击的情况，尽管其它的范数也有机会（但不是所有的都能写成混合整数规划，而且有一些更难去实现和求解）。具体来说，可以将问题写成：
\[% <![CDATA[ \DeclareMathOperator*{\minimize}{minimize} \DeclareMathOperator*{\subjectto}{subject to} \begin{split} \minimize_{z_{1,\ldots,d+1}} \;\; & (e_y - e_{y_{\mathrm{targ}}})^T z_{d+1} \\ \subjectto \;\; & \|z_1 - x\|_\infty \leq \epsilon \\ & z_{i+1} = \max\{0, W_i z_i + b_i\}, \;\; i=1,\ldots, d-1 \\ & z_{d+1} = W_d z_d + b_d \end{split} %]]>\]其中$e_i$表示单元向量，即一个向量只在第$i$个位置为1，其余位置都是0。并且我们移除了显式的$\delta$这一变量，相应地，只需要满足$z_1$和$x$的距离不超过$\epsilon$即可。
混合整数规划的形式
尽管这一形式化描述对于针对性攻击来说相对简单，但这不是求解器能够求解的形式。具体问题出在涉及max运算的等式约束上，其既不是一个凸的约束，也不是一个大多数求解器能够处理的约束。为了解决这一问题，我们需要将其转化成另一种形式，具体来说是（二元）混合整数线性规划（mixed integer linear program (MILP)）。二元MILP是指这样一个优化问题：包含线性目标（即$c^T z$，其中优化变量是$z$，系数向量是$c$）；线性等式和不等式约束（即$Az = b$和/或$Gz \leq h$，其中$A,G$是矩阵，$b,h$是向量）；以及在一些变量$z_i$上的二元约束$z_i \in {0,1}, \; \forall i \in \mathcal{I}$，其中$\mathcal{I}$是优化变量的某个子集。正是这些二元约束让问题变得困难，如果没有它们的话，问题就会变成能在多项式时间内求解的线性规划。但MILPs通常是NP难的，因此我们不指望能用其分析目前常见的神经网络模型。然而，对于小的问题，MILPs已经被广泛研究，而且与单纯地暴力尝试二元变量所有情况的方法相比，存在着更具扩展性的方法。
那如何使用线性约束和二元整数约束来表示\(z_{i+1} = \max\{0, W_i z_i + b_i\}\)呢？这里我们使用的方法被称为线性化（linearization），这是整数规划中的标准操作。假设我们已知$W_i z_i + b_i$取值的上界和下界，用$l_i$和$u_i$分别表示。此外，我们还引入一组和$z_{i+1}$尺寸相同的二元变量$v_i$。然后我们说下列不等式等价于约束\(z_{i+1} = \max\{0, W_i z_i + b_i\}\):
\[% <![CDATA[ \begin{split} z_{i+1} & \geq W_i z_i + b_i \\ z_{i+1} & \geq 0 \\ u_i \cdot v_i &\geq z_{i+1} \\ W_i z_i + b_i & \geq z_{i+1} + (1-v_i) l_i \\ v_i & \in \{0,1\}^{|v_i|} \end{split} %]]>\]想要理解上述约束，最好是通过一个简单的例子。
首先假设$W_i z_i + b_i > 0$。如果我们选择$v_i=0$，那么第三项约束则意味着$z_{i+1}\leq 0$，而第一项约束则意味着$z_{i+1} \geq W_i z_i + b_i > 0$，这会导致不可行解。因此，对于$W_i z_i + b_i > 0$，我们需要选择$v_i=1$。于是约束可以简化成
\[% <![CDATA[ \begin{split} z_{i+1} & \geq W_i z_i + b_i \\ z_{i+1} & \geq 0\\ u_i & \geq z_{i+1} \\ W_i z_i + b_i & \geq z_{i+1}, \end{split} %]]>\]其中第一项和第四项不等式共同表明$z_{i+1} = W_i z_i + b_i$，且第二项和第三项不等式永远成立。
若假设$W_i z_i + b_i < 0$，选择$v_i=1$会使得第四项约束表明$z_{i+1} \leq W_i z_i + b_i < 0$，这与第二项约束$z_{i+1}\geq 0$矛盾。因此，我们必须选择$v_i=0$，于是约束简化成
\[% <![CDATA[ \begin{split} z_{i+1} & \geq W_i z_i + b_i\\ z_{i+1} & \geq 0 \\ 0 & \geq z_{i+1} \\ W_i z_i + b_i & \geq z_{i+1} + l_i. \end{split} %]]>\]然后第二项和第三项约束表明$z_{i+1} = 0$，且第一项和第四项约束各自成立。
换言之，我们证明了在这些约束下，如果$W_i z_i + b_i > 0$则有$z_{i+1} = W_i z_i + b_i$；如果$W_i z_i + b_i < 0$则有$z_{i+1} = 0$。这正是\(z_{i+1} = \max\{0, W_i z_i + b_i\}\)的逻辑。
找到上下界
上述整数规划形式中还有一个问题待解决，即如何找到上下界$l_i$和$u_i$。从数学上讲，上述形式对于任意值的上下界都满足，例如你可以取一个特别小的下界和一个特别大的上界，然后丢给求解器去求解就完事了。但很不幸的是，整数规划求解器的实际求解时间高度依赖上下界的值。有趣的是，我们可以像计算针对性攻击那样如法炮制来计算出准确的上下界的值：正如针对性攻击会最小化$(e_y - e_{y_{\mathrm{targ}}})$乘以最有一层的输出，我们可以通过最小化中间层的激活值。例如$(l_k)_j$就是如下优化问题的解：
\[% <![CDATA[ \begin{split} \minimize_{z_{1,\ldots,k+1}} \;\; & e_j^T z_{k+1} \\ \mbox{subject to} \;\; & \|z_1 - x\|_\infty \leq \epsilon \\ & z_{i+1} = \max\{0, W_i z_i + b_i\}, \;\; i=1,\ldots, k-1 \\ & z_{k+1} = W_k z_k + b_k. \end{split} %]]>\]（这个很好理解，即使得第$k$层的输出向量在第$j$位上最小，求出的最小值即$l_k$在第$j$位下界）
尽管这样可以准确求出上下界，但是这是极不现实的，因为对每一个隐藏层单元都要求解两次整数规划（上下界各一次）。所以，更实际的做法是选择一组较为松弛的界使得能够快速完成计算。一种自然的选择是使用简单的区间界。如果我们有关于$z$的界，\(\hat{l} \leq z \leq \hat{u}\)（为了简洁起见，这里不标明层数下标），那么问题是：$W z + b$的取值能有多大或者多小？
让我们考虑其中的一位$(W z + b)_i$：
\[(W z + b)_i = \sum_{j} w_{ij} z_j + b_i\]如果想最小化该项，那么就是当$w_{ij}<0$时选择\(z_j = \hat{u}_j\)；当$w_{ij}>0$时选择\(z_j = \hat{l}_j\)。最大化的情况反之。于是，$W z + b$的边界为
\[\max\{W,0\} \hat{l} + \min\{W,0\} \hat{u} + b \leq (W z +b) \leq \max\{W,0\} \hat{u} + \min\{W,0\}\hat{l}.\]这个界其实相当松弛，因为对于不同的\(w_{i*}\)都可以选择不同的$z$（但实际上某一层的$z$取值是确定的）。但使用这个界却可以快速完成计算。根据这一策略，我们可以规范线性层的界，同理也可以给出RELU激活层的界（\(\max\{0,\hat{l}\}\leq z \leq \max\{0,\hat{u}\}\)）。此外，虽然上述推导是针对全连接层而言，但它同样适用于卷积层：卷积也是线性运算，可以用矩阵形式表示。
下面给出了计算上述边界的代码：
def bound_propagation(model, initial_bound):
l, u = initial_bound
bounds = []
for layer in model:
if isinstance(layer, Flatten):
l_ = Flatten()(l)
u_ = Flatten()(u)
elif isinstance(layer, nn.Linear):
l_ = (layer.weight.clamp(min=0) @ l.t() + layer.weight.clamp(max=0) @ u.t()
+ layer.bias[:,None]).t()
u_ = (layer.weight.clamp(min=0) @ u.t() + layer.weight.clamp(max=0) @ l.t()
+ layer.bias[:,None]).t()
elif isinstance(layer, nn.Conv2d):
l_ = (nn.functional.conv2d(l, layer.weight.clamp(min=0), bias=None,
stride=layer.stride, padding=layer.padding,
dilation=layer.dilation, groups=layer.groups) +
nn.functional.conv2d(u, layer.weight.clamp(max=0), bias=None,
stride=layer.stride, padding=layer.padding,
dilation=layer.dilation, groups=layer.groups) +
layer.bias[None,:,None,None])
u_ = (nn.functional.conv2d(u, layer.weight.clamp(min=0), bias=None,
stride=layer.stride, padding=layer.padding,
dilation=layer.dilation, groups=layer.groups) +
nn.functional.conv2d(l, layer.weight.clamp(max=0), bias=None,
stride=layer.stride, padding=layer.padding,
dilation=layer.dilation, groups=layer.groups) +
layer.bias[None,:,None,None])
elif isinstance(layer, nn.ReLU):
l_ = l.clamp(min=0)
u_ = u.clamp(min=0)
bounds.append((l_, u_))
l,u = l_, u_
return bounds
这里简单提一下（稍后会再此谈到），我们可以利用这一策略得到最后一层激活值（即最终的输出向量）的界，但我们很快就会看到，在标准模型下这些界太过松弛以至于不能提供太多信息。
以我们的卷积网络为例，以$\ell_\infty$做为初始的界（同时也需要将更新后的值约束在[0,1]之间），前向传播并观察每一层区间界：
epsilon = 0.1
bounds = bound_propagation(model_cnn, ((X - epsilon).clamp(min=0), (X + epsilon).clamp(max=1)))
print("lower bound: ", bounds[-1][0][0].detach().cpu().numpy())
print("upper bound: ", bounds[-1][1][0].detach().cpu().numpy())
lower bound: [-2998.6592 -2728.1714 -2460.193 -2534.023 -3073.9976 -2425.4656
-2874.6404 -2319.738 -2541.7148 -2524.9368]
upper bound: [2305.799 2738.133 2543.169 2636.1572 1995.432 2955.592 2503.6426
2955.3687 2904.1223 2960.8352]
这一结果告诉我们，对于扰动范围$\epsilon=0.1$，在MNIST第一个测试样本上，对应数字0类别的位，其值位于区间$[-2998, 2305]$。这个界对于确定预测值会如何变化并没有用处，但对于整数规划来说则是帮了大忙。从上述结果来看，对于深度网络，由于关于界的误差会逐层累积，因此这些界会变得越来越弱。例如，下面是针对两层全连接网络和四层全连接网络的结果：
epsilon = 0.1
bounds = bound_propagation(model_dnn_2, ((X - epsilon).clamp(min=0), (X + epsilon).clamp(max=1)))
print("lower bound: ", bounds[-1][0][0].detach().cpu().numpy())
print("upper bound: ", bounds[-1][1][0].detach().cpu().numpy())
lower bound: [-20.794655 -23.554483 -17.501461 -16.71551 -32.42357 -22.67279
-30.79689 -14.9981 -23.627762 -24.76149 ]
upper bound: [18.328127 14.5928 26.335987 28.887863 18.655684 28.70531 11.502836
32.54162 24.40407 25.302904]
epsilon = 0.1
bounds = bound_propagation(model_dnn_4, ((X - epsilon).clamp(min=0), (X + epsilon).clamp(max=1)))
print("lower bound: ", bounds[-1][0][0].detach().cpu().numpy())
print("upper bound: ", bounds[-1][1][0].detach().cpu().numpy())
lower bound: [-218.17267 -195.3326 -176.88179 -214.7243 -221.61171 -173.79388
-233.49634 -201.39543 -208.68443 -197.08893]
upper bound: [183.82648 207.23276 176.97774 205.42944 190.99425 209.2204 186.59421
231.4344 200.6849 213.99307]
比卷积网络的界要紧，但依旧没用。当然，我们稍后（在“基于边界传播的方法“中）会再次讨论这一点，因为如果训练过程得当的话，这些界实际上是有用的。目前，我们只需要关注将这些界用于整数规划即可。
最终的整数规划形式
将上述相关公式综合到一起来得到最终的整数规划形式：
\[% <![CDATA[ \begin{split} \minimize_{z_{1,\ldots,d+1}, v_{1,\ldots,d-1}} \;\; & (e_y - e_{y_{\mathrm{targ}}})^T z_{d+1} \\ \subjectto \;\; & z_{i+1} \geq W_i z_i + b_i, \;\; i=1\ldots,d-1 \\ & z_{i+1} \geq 0, \;\; i=1\ldots,d-1 \\ & u_i \cdot v_i \geq z_{i+1}, \;\; i=1\ldots,d-1 \\ & W_i z_i + b_i \geq z_{i+1} + (1-v_i) l_i, \;\; i=1\ldots,d-1 \\ & v_i \in \{0,1\}^{|v_i|}, \;\; i=1\ldots,d-1 \\ & z_1\leq x + \delta \\ & z_1 \geq x - \delta \\ & z_{d+1} = W_d z_d + b_d. \end{split} %]]>\]使用cvxpy库可以或多或少地求解上述问题。但需要注意的是，cvxpy需要与其兼容的后端求解器，这里我们使用的是Gurobi求解器，其对于学术用途是免费的。你也可以尝试兼容cvxpy的其它求解器，但在整数规划领域，免费的求解器在求解速度上远远落后于商业求解器。
这里出于简便的目的，我们针对全连接网络（仅由线性层和ReLU层组成，并假设每个线性层后面都跟随有一个ReLU层）进行求解。
import cvxpy as cp
def form_milp(model, c, initial_bounds, bounds):
linear_layers = [(layer, bound) for layer, bound in zip(model,bounds) if isinstance(layer, nn.Linear)]
d = len(linear_layers)-1
# create cvxpy variables
z = ([cp.Variable(layer.in_features) for layer, _ in linear_layers] +
[cp.Variable(linear_layers[-1][0].out_features)])
v = [cp.Variable(layer.out_features, boolean=True) for layer, _ in linear_layers[:-1]]
# extract relevant matrices
W = [layer.weight.detach().cpu().numpy() for layer,_ in linear_layers]
b = [layer.bias.detach().cpu().numpy() for layer,_ in linear_layers]
l = [l[0].detach().cpu().numpy() for _, (l,_) in linear_layers]
u = [u[0].detach().cpu().numpy() for _, (_,u) in linear_layers]
l0 = initial_bound[0][0].view(-1).detach().cpu().numpy()
u0 = initial_bound[1][0].view(-1).detach().cpu().numpy()
# add ReLU constraints
constraints = []
for i in range(len(linear_layers)-1):
constraints += [z[i+1] >= W[i] @ z[i] + b[i],
z[i+1] >= 0,
cp.multiply(v[i], u[i]) >= z[i+1],
W[i] @ z[i] + b[i] >= z[i+1] + cp.multiply((1-v[i]), l[i])]
# final linear constraint
constraints += [z[d+1] == W[d] @ z[d] + b[d]]
# initial bound constraints
constraints += [z[0] >= l0, z[0] <= u0]
return cp.Problem(cp.Minimize(c @ z[d+1]), constraints), (z, v)
为了能够快速求解MILP，我们需要考虑一个比之前使用的还要简单的网络。下列代码定义并训练了一个隐藏单元分别是50和20的三层网络（如果不想训练的话直接load即可）。
model_small = nn.Sequential(Flatten(), nn.Linear(784,50), nn.ReLU(),
nn.Linear(50,20), nn.ReLU(),
nn.Linear(20,10)).to(device)
# train model and save to disk
opt = optim.SGD(model_small.parameters(), lr=1e-1)
for _ in range(10):
train_err, train_loss = epoch(train_loader, model_small, opt)
test_err, test_loss = epoch(test_loader, model_small)
print(*("{:.6f}".format(i) for i in (train_err, train_loss, test_err, test_loss)), sep="\t")
torch.save(model_small.state_dict(), "model_small.pt")
0.203017 0.663095 0.086800 0.301833
0.078033 0.271283 0.063700 0.219485
0.058133 0.198609 0.050100 0.174730
0.046533 0.159953 0.045300 0.144941
0.039200 0.132439 0.039100 0.128218
0.034117 0.114285 0.037100 0.115507
0.029283 0.100041 0.035600 0.110929
0.026867 0.090279 0.031200 0.101059
0.024233 0.080470 0.030700 0.098896
0.022300 0.073255 0.032100 0.098225
# load model from disk
model_small.load_state_dict(torch.load("model_small.pt"))
现在我们来构建针对测试集中第一个测试样本的整数规划，通过针对性攻击将其标签判断为数字0（真实标签是7，从之前的章节的测试结果中可以找到）（代码中作者令c[2]=-1，意味着将试图将该样本预测为2，所以这里应该是作者笔误）
initial_bound = ((X[0:1] - epsilon).clamp(min=0), (X[0:1] + epsilon).clamp(max=1))
bounds = bound_propagation(model_small, initial_bound)
c = np.zeros(10)
c[y[0].item()] = 1
c[2] = -1
prob, (z, v) = form_milp(model_small, c, initial_bound, bounds)
最终使用Gurobi来求解该整数规划问题
prob.solve(solver=cp.GUROBI, verbose=True)
-4.966118334049208
算出的结果为负数，说明我们能够找到一个扰动，使得预测标签向量在目标类别上的值要比在真实类别上的值更大（这样相减之后才会是负值）。在本例中，求解MILP（含有70个二元变量）耗时不到1秒，但对于再大一点的模型来说，就会迅速变得不可解。可以观察输出向量的结果来验证：
print("Last layer values from MILP:", z[3].value)
Last layer values from MILP: [ 1.0796434 -2.62372054 8.66100952 7.01286408 -7.5661212 1.30622329
-9.93690445 3.69489119 7.47907522 -1.05974112]
的确在数字2的位上的值较大。观察扰动后的图像：
plt.imshow(1-z[0].value.reshape(28,28), cmap="gray")
验证稳健性
上述过程其实给我们提供了一种判断给定样本是否存在对抗性样本的方法，即，将每一种可能的类别作为针对性攻击类别，运行整数规划。如果在一些类别下的优化结果是负值，则说明存在对抗样本且优化过程可以将其提供给我们；相反，如果在所有类别下的优化结果都不是负值，那么说明该分类器对于该测试样本是稳健的。我们来看看如何验证在很小的扰动下不能改变类别标签：
epsilon = 0.05
initial_bound = ((X[0:1] - epsilon).clamp(min=0), (X[0:1] + epsilon).clamp(max=1))
bounds = bound_propagation(model_small, initial_bound)
for y_targ in range(10):
if y_targ != y[0].item():
c = np.eye(10)[y[0].item()] - np.eye(10)[y_targ]
prob, _ = form_milp(model_small, c, initial_bound, bounds)
print("Targeted attack {} objective: {}".format(y_targ, prob.solve(solver=cp.GUROBI)))
Targeted attack 0 objective: 5.810042236946861
Targeted attack 1 objective: 10.836173372728197
Targeted attack 2 objective: 2.039339255684599
Targeted attack 3 objective: 0.018872730723538567
Targeted attack 4 objective: 13.52836423718531
Targeted attack 5 objective: 5.146072840088203
Targeted attack 6 objective: 20.45439065167794
Targeted attack 8 objective: 3.5428745576535814
Targeted attack 9 objective: 4.274027214275465
上述优化结果都是正数，所以在$\epsilon=0.05$的$\ell_\infty$扰动下不能使得分类器误判。同理，对于那些存在对抗样本的情况，该过程可以告诉我们针对哪些类别的攻击是可能的：
epsilon = 0.1
initial_bound = ((X[0:1] - epsilon).clamp(min=0), (X[0:1] + epsilon).clamp(max=1))
bounds = bound_propagation(model_small, initial_bound)
for y_targ in range(10):
if y_targ != y[0].item():
c = np.eye(10)[y[0].item()] - np.eye(10)[y_targ]
prob, _ = form_milp(model_small, c, initial_bound, bounds)
print("Targeted attack {} objective: {}".format(y_targ, prob.solve(solver=cp.GUROBI)))
Targeted attack 0 objective: -2.545012509042315
Targeted attack 1 objective: 3.043376725812322
Targeted attack 2 objective: -4.966118334049208
Targeted attack 3 objective: -7.309837079755624
Targeted attack 4 objective: 4.741011344284811
Targeted attack 5 objective: -4.870590800081836
Targeted attack 6 objective: 7.980441149609722
Targeted attack 8 objective: -6.9536251954142365
Targeted attack 9 objective: -3.314370640372009
比如，$\epsilon=0.1$的$\ell_\infty$扰动无法使得分类器将图片预测为数字1。
内部最大化的上界（凸松弛，convex relaxations）
虽然上述使用整数规划求解内部最大化问题精确解很有价值（对于小型模型或是出于教学需要），但是该方法无法扩展到较大的网络上。这不仅仅是更多计算量的问题：这类组合优化问题的计算量在某点之后会呈指数增长。所以即使是中等大小的网络，也很容易出现优化过程永远不会结束，无论使用多么强大的计算资源。所以，精确求解是不可能的了，但如果我们能快速得到出内部最大化问题的上界，那还是有机会给出稳健性的理论保证的。例如，如果我们得到的上界表明无法实现针对性攻击，那这也是对稳健性的一种验证。
在本节，我们将会给出两种不同的方法来得到上界，其一是基于整数规划问题的凸松弛（convex relaxation）（可以给出更紧的上界，但对于较大的网络来说计算代价还是相对昂贵），另一种是基于边界传播（bound propagation）（更松弛的界，但可以更快地计算出来）。
验证问题的凸松弛
正是我们为了表示RELU运算而引入的二元整数约束才使问题变得难以求解 \(v_i \in \{0,1\}^{|v_i|}.\) 它不是一个凸集，因此在其上难以做优化，这是问题中唯一难以处理的约束。如果我们可以将其移除，那么问题就会变成一个可以快速求解的线性规划问题。于是，我们将该约束放松为位于0和1之间的小数值： \(0 \leq v_i \leq 1.\) 除此之外，问题的其他部分保持不变。这就是整数规划问题的凸松弛形式。需要强调的是，求解MILP的方法实际上都是基于此策略，也就是说求解器内部已经做了类似这样的操作。如果我们求解松弛的针对性攻击问题后得到的结果是正数，那我们就验证了该攻击是无效的；因为松弛集合严格大于原集合，所以如果在松弛集合中都不存在“对抗样本”，那么在原集合中就更不存在。
之所以这里给对抗样本加了引号，是因为松弛问题的解不再能给出真正的对抗样本（只能用来做验证）。允许$v_i$取小数意味着可以允许RELU “partially off and partially on”。通过下面的示例可以直观的理解这一点。原本的整数约束要求RELU运算前的\((W_i z_i + b_i)_j\)和运算后的\((z_{i+1})_j\)位于如下的“边界RELU”（bounded RELU）集中：
当我们放松了对$v_i$的限制使其可以取小数后，原先的边界RELU集被放松成凸集
换句话说，$z_{i+1}$不必是$(W_i z_i + b_i)$RELU运算的结果。例如，某个负值对应的RELU运算结果$z_{i+1}$可以是正数。这意味着如果我们将凸松弛问题下的解输入至神经网络，得到结果会不一样。所以，凸松弛的解不能用来构建对抗样本，相反它只是该松弛优化问题的解，能用来验证对抗样本的“不存在”。更进一步，如果凸松弛的优化结果是负值，也不能证明一定存在对抗样本。对于通过标准流程训练得到的网络，其边界相当松弛，但有趣的是，如果网络被专门训练去最小化这些界，那么我们通常可以给出关于分类器的稳健性的重要保证。
这是一个对上述内容的图解。当我们将$\ell_\infty$范数球输入到网络后，最后一层的结果（也被称为对抗多胞形，adversarial polytope）是一个非凸集。
在这个非凸集上找到最差的一点（即该点所在方向能最大化目标类别的位值和最小化真是类别的位值）就是MILP所做的事情。但当我们考虑该集合的凸松弛时，我们实际上是在考虑对抗多胞形的凸边界
在该集合上更容易优化（因为是凸的），并且由于它严格大于原集合，所以如果在这个大的集合上没有对抗样本的话，那么在对抗多胞形上也就不存在对抗样本。虽然这一凸边界可能在某些情况下非常松弛，但对于较大的网络来说，这是得到（神经网络）分类器稳健性保证的唯一可行的策略。
我们来看看凸松弛的实际效果。只需要稍微改动整数规划的代码，将$v_i$放松至0,1之间的小数值。
def form_lp(model, c, initial_bounds, bounds):
linear_layers = [(layer, bound) for layer, bound in zip(model,bounds) if isinstance(layer, nn.Linear)]
d = len(linear_layers)-1
# create cvxpy variables
z = ([cp.Variable(layer.in_features) for layer, _ in linear_layers] +
[cp.Variable(linear_layers[-1][0].out_features)])
v = [cp.Variable(layer.out_features) for layer, _ in linear_layers[:-1]]
# extract relevant matrices
W = [layer.weight.detach().cpu().numpy() for layer,_ in linear_layers]
b = [layer.bias.detach().cpu().numpy() for layer,_ in linear_layers]
l = [l[0].detach().cpu().numpy() for _, (l,_) in linear_layers]
u = [u[0].detach().cpu().numpy() for _, (_,u) in linear_layers]
l0 = initial_bound[0][0].view(-1).detach().cpu().numpy()
u0 = initial_bound[1][0].view(-1).detach().cpu().numpy()
# add ReLU constraints
constraints = []
for i in range(len(linear_layers)-1):
constraints += [z[i+1] >= W[i] @ z[i] + b[i],
z[i+1] >= 0,
cp.multiply(v[i], u[i]) >= z[i+1],
W[i] @ z[i] + b[i] >= z[i+1] + cp.multiply((1-v[i]), l[i]),
v[i] >= 0,
v[i] <= 1]
# final linear constraint
constraints += [z[d+1] == W[d] @ z[d] + b[d]]
# initial bound constraints
constraints += [z[0] >= l0, z[0] <= u0]
return cp.Problem(cp.Minimize(c @ z[d+1]), constraints), (z, v)
并使用和求解整数规划相同的样本来比较两者的不同
initial_bound = ((X[0:1] - epsilon).clamp(min=0), (X[0:1] + epsilon).clamp(max=1))
bounds = bound_propagation(model_small, initial_bound)
c = np.eye(10)[y[0].item()] - np.eye(10)[2]
prob, (z, v) = form_lp(model_small, c, initial_bound, bounds)
print("Objective value:", prob.solve(solver=cp.GUROBI))
print("Last layer from relaxation:", z[3].value)
Objective value: -26.119774634225454
Last layer from relaxation: [-0.54239629 -3.68197597 21.97748102 6.17979398 -6.16385997 -0.27719886
-2.25631931 -4.14229362 9.61796699 -5.31716916]
凸松弛问题的结果表明我们可以将数字2的预测值变得（通过扰动）比真实类别数字7的预测值大26.1198。求解松弛问题得到的初始扰动为：
plt.imshow(1-z[0].value.reshape(28,28), cmap="gray")
不同于整数规划的情况，如果我们将该扰动输入至神经网络，得到激活值于松弛问题下求解的激活值不同（虽然在这个例子下，它的确是一个对抗样本。一般情况下是无法保证这一点的）
print("Last layer from model:",
model_small(torch.tensor(z[0].value).float().view(1,1,28,28).to(device))[0].detach().cpu().numpy())
Last layer from model: [ 1.1390655 -2.760746 8.2507715 5.9955297 -5.9659457 0.492265
-8.71024 3.8318188 6.594626 -1.3704427]
作为最后一个例子，我们利用该方法来验证网络对于很小的扰动的稳健性。同样，和之前的代码相同，只是使用了凸松弛。
epsilon = 0.02
initial_bound = ((X[0:1] - epsilon).clamp(min=0), (X[0:1] + epsilon).clamp(max=1))
bounds = bound_propagation(model_small, initial_bound)
for y_targ in range(10):
if y_targ != y[0].item():
c = np.eye(10)[y[0].item()] - np.eye(10)[y_targ]
prob, _ = form_lp(model_small, c, initial_bound, bounds)
print("Targeted attack {} objective: {}".format(y_targ, prob.solve(solver=cp.GUROBI)))
Targeted attack 0 objective: 8.679119731697458
Targeted attack 1 objective: 12.233313266959202
Targeted attack 2 objective: 4.120827654631229
Targeted attack 3 objective: 1.857739424985871
Targeted attack 4 objective: 16.236514322580383
Targeted attack 5 objective: 7.964241745873003
Targeted attack 6 objective: 23.459304334257965
Targeted attack 8 objective: 6.10314963105955
Targeted attack 9 objective: 5.581247281946292
因此我们能验证在$\epsilon=0.02$的$\ell_\infty$扰动下无法欺骗到分类器。而且，不像整数规划，我们可以在较大的模型上使用这一方法。这里的验证对象是4层DNN，整数规划对此无能为力。
epsilon = 0.01
initial_bound = ((X[0:1] - epsilon).clamp(min=0), (X[0:1] + epsilon).clamp(max=1))
bounds = bound_propagation(model_dnn_4, initial_bound)
for y_targ in range(10):
if y_targ != y[0].item():
c = np.eye(10)[y[0].item()] - np.eye(10)[y_targ]
prob, _ = form_lp(model_dnn_4, c, initial_bound, bounds)
print("Targeted attack {} objective: {}".format(y_targ, prob.solve(solver=cp.GUROBI)))
Targeted attack 0 objective: 5.480101103982182
Targeted attack 1 objective: 2.196709825912893
Targeted attack 2 objective: 0.4646296745679468
Targeted attack 3 objective: -2.264958127925656
Targeted attack 4 objective: 9.8073503196978
Targeted attack 5 objective: 1.4951834763453595
Targeted attack 6 objective: 12.835283301376197
Targeted attack 8 objective: -0.39325769328255866
Targeted attack 9 objective: -1.9890586254866953
凸松弛的更快解法
虽然上述的凸松弛被认为是要比整数规划“快”，但要注意，上述代码都是在测试一个MNIST样本，在我们的机器上耗时数秒完成。如果想验证所有的MNIST样本，并且模型要是再大一点的话，那这个时间开销就无法接受了。
由于这一问题，大多数针对凸松弛方法的研究工作实际上都基于那些能快速求解凸松弛问题的方法。详细的细节超出了本教程的范围，但结果是，使用凸对偶（convex duality）以及对优化问题做一些处理，能够利用单次反向传播（single backward pass through the network）快速计算出松弛目标的可证下界（这反过来可以使原始问题更加“松弛”）。相比于我们之前提到的简单的边界传播，该过程非常快，我们甚至可以用其计算出更紧的区间界。然而，找到验证模型的复杂度和边界的松紧程度之间的权衡依然是一个开放问题。
基于间隔传播的边界(Interval-propagation-based)
在继续之前，我们想给出最后一种得到内部优化问题上界的方法，其更直接地基于我们之前讨论过的边界传播的技术。这些界要比那些基于凸松弛的更弱，并且这些界甚至不能像基于凸松弛方法那样给出“假的”对抗样本，但它们的优势是效率极高。实际上，最近的工作表明，在训练稳健的网络时（我们将在下一章中广泛讨论该主题），更值得使用这些边界传播的技术。不过关于这两类方法的对比依然是一个活跃的研究主题。
我们再次回到之前讨论过的在网络中传递区间界的方法。具体的，给定某个界$\hat{l} \leq z \leq \hat{u}$，则有
\[l \leq Wx + b \leq u\]其中
\[% <![CDATA[ \begin{split} l & = \max\{W,0\} \hat{l} + \min\{W,0\} \hat{u} + b \\ u & = \max\{W,0\} \hat{u} + \min\{W,0\}\hat{l} + b. \end{split} %]]>\]那么这些界本身就可以给出网络能取到的边界。然而，将这些界传播到最后一层其实是浪费的，因为我们最终希望最小化的是最后一层中的某个特定的线性函数（我的理解是，最后一层激活值的界是多少与优化目标并不直接相关）。因此，我们所做的就是将这些界传播到倒数第二层，然后解析地求解最后一层的最小化问题。
再具体一点，对于一个$d$层的网络，我们的目标是最小化最后一层的某个线性函数$c^T z_{d+1}$。如果我们有关于倒数第二层的区间界，那么我们可以简单地求解这个优化问题
\[% <![CDATA[ \begin{split} \minimize_{z_d,z_{d+1}} \;\; & c^T z_{d+1} \\ \subjectto \;\; & z_{d+1} = W_d z_d + b_d \\ & \hat{l} \leq z_d \leq \hat{u}. \end{split} %]]>\]具体地，如果我们只是将变量$z_{d+1}$用第一个约束替换掉，那么问题等价于
\[% <![CDATA[ \begin{split} \minimize_{z_d} \;\; & c^T (W_d z_d + b_d) \equiv (W_d^T c)^T z_d + c^T b_d \\ \subjectto \;\; & \hat{l} \leq z_d \leq \hat{u}. \end{split} %]]>\]这正是我们之前遇到过的有界约束下最小化线性函数的问题。其解析解为当$(W_d^T c)_j > 0$时选择$(z_d)_j = \hat{l}_j$，反之$(z_d)_j = \hat{u}_j$。于是，最优解为
\[\max\{c^T W_d c,0\} \hat{l} + \min\{c^T W_d,0\} \hat{u} + c^T b_d\]其非常好的一点是可以轻易实现同时计算多个$c$（即同时考虑所有可能的攻击）、多个样本的界。
def interval_based_bound(model, c, bounds):
# requires last layer to be linear
cW = c.t() @ model[-1].weight
cb = c.t() @ model[-1].bias
l,u = bounds[-2]
return (cW.clamp(min=0) @ l.t() + cW.clamp(max=0) @ u.t() + cb[:,None]).t()
用该界来验证$\epsilon=0.02$的情况
epsilon = 0.02
initial_bound = ((X[0:1] - epsilon).clamp(min=0), (X[0:1] + epsilon).clamp(max=1))
bounds = bound_propagation(model_small, initial_bound)
C = -torch.eye(10).to(device)
C[y[0].item(),:] += 1
print(interval_based_bound(model_small, C, bounds).detach().cpu().numpy())
[[ -7.426559 -2.0292485 -8.31501 -15.399408 2.3813596 -10.803006
5.575859 0. -13.850781 -11.277755 ]]
这些项是每一类的优化目标的下界（是原始最大化问题的上界）。正如我们所预计的那样，确实比解决凸松弛优化问题得到的界要松得多。和凸松弛不同的是，每个单独的激活值都假定前一层可以采用一组单独的值来最小化或最大化该激活值（这就是它更松的原因，同时也是它更快的原因）。这些界可以快速计算并且已经在PyTorch中实现，这意味着我们可以通过它们的反向传播，来训练以稳健性为标准的网络。这部分会在下一章中展开讲解。
总结
这一章中我们介绍了很多内容，但记住这些不同的方法都是为了优化（或近似优化）稳健性目标的内部优化问题（在一些例子中虽然只优化了某一位，但实质是一样的）。在一些情况中，这些方法给出了实际的攻击扰动应该是什么，但在另一些情况中，它们仅仅给出了优化目标的边界，可以用来验证分类器的稳健性。重要的是，我们希望透过这些方法来告诉大家内部优化问题的本质：其远不仅仅是找到对抗样本，而应该更多的强调对抗稳健性。
个人总结
针对神经网络的内部最大化问题提出了三种策略：
- 第一种方法的思路很简单，就是通过梯度下降来最大化目标函数（需要处理的是扰动约束的限制）。但这种方法依赖于对抗样本的结果，也就是神经网络优化的结果。为了找到更好的对抗样本，提出了不同的梯度下降方法。
- 对于小型的网络，可以将内部最大化问题转化为整数规划问题并精确求解。它可以告诉我们到底能不能针对某一类做针对性攻击（同理，分类器对于不同类别的攻击是否稳健）。最重要的是，这一结果是精确的。
- 遗憾的是，精确求解无法扩展到较大的模型上。因此考虑来其凸松弛问题（主要是针对整数规划中二元变量的放松）。虽然不再能给出精确解（因此也无法构建出对抗样本），但能够快速计算出优化目标的上界，从而用于验证对抗攻击的稳健性。不过求解凸松弛问题还是慢（虽然也有一些快速的解法），只能针对某一类别的攻击、某一个样本单独计算。因此还有一种基于边界传播的方法来得到上界：通过假设每一层的激活值都是根据前一层激活值的最值来解析地计算得出最大化问题的上界。因此其得到的界会更松，但速度会非常快（不需要再做整数规划了）。