综述：使用基于梯度的方法对深度学习模型进行解包（即模型结构的分析（教程）

《Analytica Chimica Acta》：Unboxing deep learning models using gradient-based methods – A tutorial

　　谢尔盖·库切里亚夫斯基
奥尔堡大学，尼尔斯·玻尔街8号，埃斯比约，6700，丹麦

摘要
深度学习模型，特别是卷积神经网络，在涉及光谱数据的回归和分类任务中表现出色。然而，它们的黑箱性质被认为是限制这些方法应用性的一个主要缺点。本文提供了一个全面的教程，介绍了基于梯度的方法来解决这个问题，这些方法提供了关于各个预测因子对预测值贡献的信息。这些方法基于对神经网络内部数据流动的理解，因此可以用来评估变量的重要性和整个模型的有效性。

本文介绍了神经网络中的梯度概念，详细推导了显著性图（普通梯度）和积分梯度的公式，并通过数值示例逐步说明了计算过程。理论部分之后是基于模拟数据和两个著名光谱数据集（Tecator和Beer）的例子。这些例子表明，当底层模型表现良好时，积分梯度提供的权重分数稳定且可解释，并且在不同的模型架构和训练迭代中保持一致。本文还提出了一种基于梯度响应乘积（GIR）的 slight 修改，以提高该方法对噪声数据集和具有系统性但不相关变异的数据集的信息量。

主要内容附带了 Jupyter 笔记本形式的 Python 代码，使这些方法能够用于实际应用。

1. 引言
深度学习（DL）是人工神经网络（ANN）的一个子领域，在2000年代和2010年代随着计算能力和数据可用性的提高而变得重要。与早期使用单个隐藏层和非线性激活函数的电位器模型不同，DL 架构采用多层不同类型的神经元来提取分层特征表示[1]。这些表示通过自动特征工程增强了回归和分类任务。其中最常用的 DL 架构之一是卷积神经网络（CNN），其特征提取主要通过卷积层来完成，后者捕获空间和时间模式[2]。尽管取得了成功，但 DL 模型常常因其缺乏透明度而受到批评。它们通常被视为“黑箱”，将输入映射到输出而不披露其内部决策过程。十多年前，人们提出了可视化方法来理解 CNN 中卷积层的作用和重要性[3]。此后，开发了额外的可解释性技术来突出输入数据中对预测起作用的关键区域，同时抑制不相关的背景信息[4][5][6][7]。这一机器学习（ML）领域变得非常重要，并且现在有一个专门的名称：可解释性 AI（XAI）[8]。

评估模型对预测响应值或类别分配（从而整体重要性）贡献的现代方法可以分为两类：模型特定的和模型不可知的。模型特定的方法直接依赖于模型架构、参数和相关结果。例如，回归模型的系数、PLS 中的变量重要性以及决策树中的基尼重要性[9][10]。相比之下，模型不可知的方法更为通用，可以用于任何分类或回归模型。其中最著名的方法是排列特征重要性（PFI）[11][12]，它主要基于使用原始数据和其中一个预测因子值被随机排列的数据时模型的性能指标（如均方根预测误差、决定系数或分类准确率）的比较。性能下降越大，该预测因子的就越重要。然而，PFI 在化学计量学中并不广泛使用，因为它不适用于包含共线预测因子的数据集。另一个值得注意的是，PFI 提供的是整个模型的全局重要性分数，而不评估局部重要性（即哪些变量对特定观测值的预测有贡献）。

最近一种迅速流行起来的模型不可知方法是 Shapley 加性解释（SHAP）[13]。简而言之，为了找到特定预测因子的贡献，SHAP 会查看将该预测因子添加到其他预测因子子集时模型输出的变化，并计算给定子集的边际贡献。通过对不同子集重复此过程并计算贡献值的加权平均值，得到最终的 SHAP 分数。所有 SHAP 分数之和（加上基线期望值）恰好等于该数据对象的模型输出。因此，SHAP 提供了直接可解释的贡献。

如上所述，SHAP 是一种计算要求较高的方法。虽然对于随机森林和梯度提升等模型存在高效的实现，但将其应用于深度学习模型仍然非常耗时。已经提出了更快的近似方法，如 DeepSHAP，但它们通常难以处理光谱数据中典型的高维共线性问题。对于深度学习模型，更有效的替代方法是基于梯度的归因方法[3][14][15][16][17]。这些方法是模型特定的，因为它们利用了神经网络的可微性，并使用反向传播来评估输入变化对模型输出的影响。本文旨在为基于梯度的预测因子重要性评估提供理论基础，并通过数值示例进行支持，并演示其在光谱数据中的应用。此外，本文还提出了一种简单的方法，以改善受噪声和系统性但不相关变异影响的光谱数据的基于梯度的重要性度量的信息量。

关于理论部分和基于模拟数据的示例的 Python 代码作为 Jupyter 笔记本提供在补充材料中。

2. 理论
理论背景通过数值示例进行说明，所有计算都可以手工完成。建议重现这些计算以更好地理解理论。

在下面的方程中，n 表示数据点（样本、观测值）的数量，而 m 表示预测因子（特征、描述符）的数量。因此，预测因子矩阵 X 有 n 行和 m 列。响应被组织为列向量 y，包含 n 个值。索引 i 用于指向 X 的特定观测值、行或 y 的值。如果只考虑单个观测值，则省略索引 i。索引 j 用于指向特定的预测因子（X 的列）。下面是一个简单的数据集示例，其中 n=4 和 m=3：
(1) X = [1.0, 2.0, 3.0, 4.0, 6.2, 2.0, 7.0, 1.0, 1.1, 5.3, 2.0, 4.0]
(2) y = [10.0, 9.0, 11.0, 14.0, 2.1]

2.1 单个神经元的梯度和梯度下降
ANN 的核心元素是神经元，它通常有多个输入 x 和一个输出 y?。这里 x = [x1, …, xm] 是一个行向量，例如预测因子矩阵 X 的一行。输出分两步计算。首先，通过对输入进行加权求和并加上 bias：
(3) z = x1w1 + … + xmwm + b
然后应用激活函数 fa：
(4) y? = fa(z)

最常见的激活函数之一是修正线性单元（ReLU）：
(5) fa(z) = 0 如果 z ≤ 0；
z = 1 如果 z > 0

神经元的权重和 bias 是必须在训练过程中优化的参数。最初，它们被设置为随机值，然后使用梯度下降或更先进的基于 GD 的优化方法（例如 Adam[18]）逐渐调整。

训练的梯度计算需要一个损失函数，该损失函数又基于获得的输出值 y? 和预期的输出值 y 来计算。例如，在回归的情况下，最常见的损失函数是对一组观测值计算的均方误差（MSE）：
(6) L = 1/n ∑i=1n (yi ? y?i)2

如果只考虑一个数据点，可以省略索引 i 来简化表达式：
(7) L = (y ? y?)2

在这种情况下，特定权重 wj 的梯度可以通过对 L 关于 wj 的偏导数来计算：
(8) ?L/?wj = ?L/?y? ?y? ?z ?z ?wj

如果没有使用激活函数，则 y? = z，因此 ?y?/?z = 1，从而得到以下（简化的）公式：
(9) ?L/?wj = (y? ? y) × xj

2.2 反向传播
在前馈网络中，神经元组合成层，因此一层有多少个神经元就会有多少个输出。层之间是互联的——前一层的输出是下一层的输入，依此类推。例如，如果有 m=30 个预测因子，并且我们在第一层使用 10 个神经元，那么该层将有 30 个输入和 10 个输出，下一层也应该有 10 个输入，依此类推。这种架构也称为全连接层。

为了计算各层的梯度并优化每个神经元的权重，损失被反向传播[19]。考虑一个有 3 个预测因子和 1 个响应的问题。为了解决这个问题，将使用一个两层网络，第一层有 2 个神经元，第二层有 1 个神经元，如图 1 所示。为简单起见，这里不使用激活函数。假设所有权重和 bias 随机设置为以下值：

下载：下载高分辨率图像（194KB）
下载：下载全尺寸图像

图 1. 一个具有两层和三个神经元的人工神经网络。wlni 表示第 l 层的第 n 个神经元与第 i 个输入之间的权重；bln 表示第 l 层的第 n 个神经元的 bias。

第一层的两个神经元：
w11 = [1.0, 0.5, 2.0], b11 = 0.0
w12 = [2.0, 0.5, 1.0], b12 = 0.0

第二层的一个神经元：
w21 = [1.0, 0.5], b21 = 0.0

如果我们取一个数据点 {x, y}，其中 x = [1.0, 2.0, 3.0] 且 y = 10.0，那么前向传播过程如下：
第一层：
y?11 = x1w111 + x2w112 + x3w113 + b11 = 1.0 × 1.0 + 2.0 × 0.5 + 3.0 × 2.0 + 0.0 = 8.0
y?12 = x1w121 + x2w122 + x3w123 + b12 = 1.0 × 2.0 + 2.0 × 0.5 + 3.0 × 1.0 + 0.0 = 6.0

第二层：
y? = y?11w211 + y?12w212 + b21 = 8.0 × 1.0 + 6.0 × 0.5 + 0.0 = 11.0

误差为：
e = y ? y? = 10.0 ? 11.0 = ?1.0

输出层的梯度：
?L/?w211 = 2(y? ? y)
y?11 = 2 × 1.0 × 8.0 = 16
?L/?w212 = 2(y? ? y)
y?12 = 2 × 1.0 × 6.0 = 12
?L/?b21 = 2(y? ? y) = 2 × 1.0 = 2.0

现在来计算第一层神经元的梯度。首先反向传播损失：
?L/?y?11 = ?L/?y? ?y? ?y?11 = 2(y? ? y)
?L/?y?12 = ?L/?y? ?y? ?y?11
?L/?y?12 = ?L/?y? ?y?12
?L/?y?12 = 2(y? ? y)
?y?/?y?11 = w211
?y?/?y?12 = w212

然后可以将损失反向传播：
?L/?y?11 = 2(y? ? y) × w211 = 2 × 1 × 1.0 = 2.0
?L/?y?12 = 2(y? ? y) × w212 = 2 × 1 × 0.5 = 1.0

必须强调的是，这种误差表达式 2(y? ? y)w 是因为没有使用激活函数。否则，该表达式需要乘以该函数的导数。

现在可以计算第一层神经元的梯度。通用表达式如下：
?L/?w1ni = ?L/?y?1n
?L/?w1ni = ?L/?y?1n × xi
?L/?b1n = ?L/?y?1n
?L/?b1n = ?L/?y?1n

对于第一个神经元：
?L/?w111 = 2.0 × 1.0 = 2.0
?L/?w112 = 2.0 × 2.0 = 4.0
?L/?w113 = 2.0 × 3.0 = 6.0
?L/?b11 = 2.0 × 1.0 = 2.0

对于第二个神经元：
?L/?w121 = 1.0 × 1.0 = 1.0
?L/?w122 = 1.0 × 2.0 = 2.0
?L/?w123 = 1.0 × 3.0 = 3.0
?L/?b12 = 1.0 × 1.0 = 1.0

请注意，在传播过程中首先计算原始梯度；缩放仅在更新步骤中进行，新权重根据 (12) 中所示进行计算。例如，如果 α = 0.001，则：
w111 = 1.0 ? 2.0 × 0.001 = 0.998
w112 = 0.5 ? 4.0 × 0.001 = 0.496
w113 = 2.0 ? 6.0 × 0.001 = 1.994
b11 = 0.0 ? 2.0 × 0.001 = ?0.002
w121 = 2.0 ? 1.0 × 0.001 = 1.999
w122 = 0.5 ? 2.0 × 0.001 = 0.498
w123 = 1.0 ? 3.0 × 0.001 = 0.997
b12 = 0.0 ? 1.0 × 0.001 = ?0.002

现在可以使用更新后的权重来计算新的损失值。第一层：
y?11=1.0×0.998+2.0×0.496+3.0×1.994?0.002=7.970
y?12=1.0×1.999+2.0×0.498+3.0×0.997?0.001=5.985

第二层：
y?=7.970×0.984+5.985×0.488?0.002=10.761

新的误差为：
e=y?y?=10.000?10.761=?0.761

为了计算损失值，需要对误差值进行平方，优化前的误差值为（?1.00）2=1.00，优化后的误差值为（?0.761）2=0.579——明显有所改善。

然而，训练集的损失值几乎总是会变小（前提是学习率不要太大）。因此，为了控制训练过程，通常会考虑使用验证集（一组未用于计算梯度的数据点）来计算损失值。权重优化过程会持续进行，直到观察到验证集的损失值几乎没有变化。

2.3. 突出性图
现在让我们看看基于梯度的方法如何帮助我们估计预测因子的贡献及其重要性，首先从突出性图（也称为普通梯度）[3]开始。首先，重点从更新权重转移到将最终预测值y?归因于原始输入x。因此，方程（8）可以重写为：
(16) ?L/?xj = ?L/?y? ?y?/?z ?z/?xj
或者，如果不使用激活函数，则为：
(17) ?L/?xj = ?L/?y? ?y?/?xj

如果使用MSE作为损失函数，则为：
(18) ?L/?xj = 2(y??y) / ?y?/?xj

由于模型已经训练完成，损失值已经不重要了，因为我们的目标是估计每个预测因子对预测值的贡献。此外，如果模型表现良好，那么y≈y?，因此?L/?y?≈0。这将导致即使预测因子很重要，梯度也很小（几乎为零）。为了解决这个问题，选择模型输出y?作为感兴趣的标量，我们只计算?y?/?xj。在计算上，这相当于从输出端反向传播一个单位信号（因此输出相对于输入的梯度被计算出来）：
(19) ?y?/?xj = ?y?/?z ?z/?xj

让我们以前面的例子为例，包括权重、预测因子向量或响应值，但现在假设模型已经训练完成，因此不需要优化权重。让我们计算每个xj的梯度?y?/?xj。这里计算简化了，因为没有使用激活函数（y?=z）。在最后一层的输出端，我们有初始信号1，该信号被反向传播到前一层：
?y?/?y11 = 1×w211 = 1.0
?y?/?y12 = 1×w212 = 0.5

然后可以对第一层的两个神经元做同样的事情，反向传播上面的计算结果（第一个神经元为1.0，第二个为0.5）：
?y?11/?x1 = 1.0×w111 = 1.00
?y?11/?x2 = 1.0×w112 = 0.50
?y?11/?x3 = 1.0×w113 = 2.00
?y?12/?x1 = 0.5×w121 = 1.00
?y?12/?x2 = 0.5×w122 = 0.25
?y?12/?x3 = 0.5×w123 = 0.50

每个xj的总梯度可以通过叠加各个神经元的梯度来计算：
g1 = ?y?/?x1 = 1.00 + 1.00 = 2.00
g2 = ?y?/?x2 = 0.50 + 0.25 = 0.75
g3 = ?y?/?x3 = 2.00 + 0.50 = 2.50

在上面的例子中，计算出的梯度值不依赖于预测因子值。这是因为没有使用激活函数来简化计算。实际上，必须先正向传递信号，保存激活函数的状态，然后再将1反向传播回去。这意味着计算出的?y?/?xj对于不同的输入会有所不同，因此梯度实际上是输入信号x的函数：gj(x) = ?y?/?xj(x)

现在可以计算突出性Sj，作为相应梯度的绝对值或直接使用梯度（即普通梯度突出性）。

突出性的一个缺点是没有考虑预测因子的尺度。例如，如果xj乘以100，相应的权重除以100，预测结果保持不变，但突出性会缩小100倍。为了克服这个问题，可以通过将梯度gj乘以输入值xj来计算重要性得分（这种方法称为梯度×输入）[20]。

最后，可以将整个数据集（例如训练集或测试集）的突出性值或梯度×输入值进行汇总，以获得每个预测因子的重要性得分：
(20) ISj = 1/n ∑i=1^n |xijgj(xi)|

普通梯度突出性和梯度×输入都是快速估计预测因子重要性的方法。然而，它们有两个关键局限性[20]：
1. 饱和问题。假设使用一个输出为Sigmoid函数的分类模型。如果预测值y?很大，接近最大值，那么Sigmoid函数的输出约为0.999，因此局部梯度几乎为零（当Sigmoid函数接近1时几乎是平的，其导数接近零）。因此，突出性值会很小，表明输入“不重要”，即使它是导致预测值很高的原因。
2. 零梯度问题：如果ReLU激活神经元在特定x处“关闭”，其导数为0。普通梯度方法会完全忽略该神经元之后的任何特征，即使这些特征只需微小的变化就能“激活”它。此外，对于较大的xj，当它乘以一个负权重并使总和低于零时，ReLU“门”会关闭，梯度变为零，从而隐藏特征的影响。

这两个问题可以通过更先进的基于梯度的方法——积分梯度来解决。

2.4. 积分梯度
积分梯度是一种解决普通梯度的饱和和零梯度问题的方法，通过使用积分[21]。不是计算一个梯度?y?/?xj，而是为一组位于感兴趣点x和基线x′之间的点x(k)计算一系列梯度，然后进行积分。基线是一个作为共同参考点的人造数据点。通常使用一个全为零的向量作为基线，即x′=[0,0,0]；然而，在光谱数据的情况下，平均光谱也是一个不错的选择（详见下面的某个部分）。

选择一个基线后，让我们在x′和x之间创建一组K个点（有时称为路径）。以前面的例子为例，在基线和原始数据点之间再创建一个点，得到以下路径（K=2）：
x(1) = [0.5,1.0,1.5]
x(2) = [1.0,2.0,3.0]
或者换句话说：
(21) x(k) = pkx + (1?pk)x′，其中pk=k/K。显然，当k=0时，x(k)=x′，因此x(0)也可以作为基线点的表示。

所有点x(k)的梯度按照前面描述的方法计算。我们的输出端初始信号为1。现在可以将其反向传播到前一层：
?y?/?y11? = 1×w211 = 1.0
?y?/?y12? = 1×w212 = 0.5

第一层的两个神经元：
?y?11/?x1 = 1.0×w111 = 1.00
?y?11/?x2 = 1.0×w112 = 0.50
?y?11/?x3 = 1.0×w113 = 2.00
?y?12/?x1 = 0.5×w121 = 1.00
?y?12/?x2 = 0.5×w122 = 0.25
?y?12/?x3 = 0.5×w123 = 0.50

每个xj的总梯度可以通过叠加各个神经元的梯度来计算：
g1 = ?y?/?x1 = 1.00 + 1.00 = 2.00
g2 = ?y?/?x2 = 0.50 + 0.25 = 0.75
g3 = ?y?/?x3 = 2.00 + 0.50 = 2.50

在上面的例子中，计算出的梯度值不依赖于预测因子值。这是因为为了简化计算，没有使用激活函数。实际上，必须先正向传递信号，保存激活函数的状态，然后再将1反向传播回去。这意味着对于不同的输入，计算出的?y?/?xj会有所不同，因此每个k的梯度是输入信号x的函数：gj(x) = ?y?/?xj(x)

现在可以计算突出性Sj，作为相应梯度的绝对值或直接使用梯度（即普通梯度突出性）。

突出性的一个缺点是没有考虑预测因子的尺度。例如，如果xj乘以100并且相应的权重除以100，预测结果保持不变，但突出性会缩小100倍。为了克服这个问题，可以通过将梯度gj乘以输入值xj来计算重要性得分（这种方法称为梯度×输入）[20]。

最后，可以汇总整个数据集（例如训练集或测试集）的突出性值或梯度×输入值，以获得每个预测因子的重要性得分：
(20) ISj = 1/n ∑i=1^n |xijgj(xi)|

普通梯度突出性和梯度×输入都是快速估计预测因子重要性的方法。然而，它们有两个关键局限性[20]：
1. 饱和问题。
2. 零梯度问题：如果ReLU激活神经元在特定x处“关闭”，其导数为0。普通梯度方法将完全忽略该神经元之后的任何特征，即使这些特征只需微小的变化就能“激活”它。此外，对于较大的xj，当它乘以一个负权重并将总和推到零以下时，ReLU“门”会关闭，梯度变为零，隐藏了特征的影响。

这两个问题可以通过使用更先进的基于梯度的方法——积分梯度来解决。

2.5. 积分梯度
积分梯度是一种通过使用积分来解决普通梯度的饱和和零梯度问题的方法[21]。不是计算一个梯度?y?/?xj，而是为位于感兴趣点x和基线x′之间的一组点x(k)计算一系列梯度，然后进行积分。基线是一个人工数据点，用作共同参考点。通常使用全为零的向量作为基线，即x′=[0,0,0]；然而，在光谱数据的情况下，平均光谱也是一个好的选择（详见下面的某个部分）。

选择一个基线后，让我们在x′和x之间创建一组K个点（有时称为路径）。以前面的例子为例，在基线和原始数据点之间再创建一个点，得到以下路径（K=2）：
x(1) = [0.5,1.0,1.5]
x(2) = [1.0,2.0,3.0]
或者换句话说：
(21) x(k) = pkx + (1?pk)x′，其中pk=k/K。显然，当k=0时，x(k)=x′，因此x(0)也可以用作基线点的表示。

所有点x(k)的梯度按照前面描述的方法计算。我们输出端的初始信号为1。现在可以将其反向传播到前一层：
?y?/?y11? = 1×w211 = 1.0
?y?/?y12? = 1×w212 = 0.5

第一层的两个神经元：
?y?11/?x1 = 1.0×w111 = 1.00
?y?11/?x2 = 1.0×w112 = 0.50
?y?11/?x3 = 1.0×w113 = 2.00
?y?12/?x1 = 0.5×w121 = 1.00
?y?12/?x2 = 0.5×w122 = 0.25
?y?12/?x3 = 0.5×w123 = 0.50

每个xj的总梯度可以通过叠加各个神经元的梯度来计算：
g1 = ?y?/?x1 = 1.00 + 1.00 = 2.00
g2 = ?y?/?x2 = 0.50 + 0.25 = 0.75
g3 = ?y?/?x3 = 2.00 + 0.50 = 2.50

在上面的例子中，计算出的梯度值不依赖于预测因子值。这是因为为了简化计算，没有使用激活函数。实际上，如前所述，必须先正向传递信号，保存激活函数的状态，然后再将1反向传播回去。这意味着对于不同的输入，计算出的?y?/?xj会有所不同，因此对于每个k，梯度是输入x(k)的函数（类似于前面得到的普通梯度）：gj(k) = ?y?/?xj(x(k))

为了估计每个变量的重要性，需要沿着路径积分梯度（因此计算积分梯度IG）。在我们的例子中有两个步骤：
IGj = (xj(1)?xj(0))gj(1) + (xj(2)?xj(1))gj(2) = (x?x′)1/2
gj(1) + gj(2)

或者，考虑到我们的例子中没有使用激活函数且gj(1)=gj(2)=gj：
IGj = (xj?xj′)gj

在我们的例子中：
IG1 = (1.0?0)×2.00 = 2.0
IG2 = (2.0?0)×0.75 = 1.5
IG3 = (3.0?0)×2.50 = 7.5

有趣的是，积分梯度的总和给出了预测值y?：
IG1 + IG2 + IG3 = 2.0 + 1.5 + 7.5 = 11

换句话说，每个IGj给出了每个预测因子xj对预测响应值的直接贡献。然而，这是因为我们使用的基线为零。更准确地说，IG的总和等于我们的预测点与基线处的预测值之间的差（这被称为完整性公理）：
∑IG = f(x)?f(x′)

值得一提的是，由于没有使用激活函数，上面的例子已经简化了，即使我们在x和x′之间没有额外的点，得到的IG也是相同的。实际上，对于不同的x(k)，gj(k)的值可能会有很大差异，因此使用足够的点数K以确保积分在数值上准确是非常重要的，因为在这种情况下，不是计算真正的积分，而是基于一组离散值的数值近似（也称为黎曼和）。

这是最终的IG计算算法：
1. 定义基线x′和路径中的点数K（通常K≥50）。
2. 对于每个步骤k∈{1,…,K}：
- 计算x(k) = pkx + (1?pk)x′。
- 将x(k)通过网络传递并保存所有内部激活值(z)。
- 使用权重和保存的激活值反向传播信号1.0，以计算每个神经元在特定状态下的激活函数导数。
- 收集得到的输入梯度gj(k)。
3. 通过对所有步骤收集的输入梯度gj(k)进行平均，并将平均梯度按原始输入x与基线x′之间的差进行缩放，来计算每个预测因子的积分梯度。
数学上可以表示为：
(22) IGj = (xj?xj′)1/2 ∑k=1^Kgj(k)

现代的ANN库（如TensorFlow或PyTorch）将所有计算向量化，因此所有步骤可以用一行代码完成，无需使用循环。大多数库都实现了积分梯度方法，例如在PyTorch中，Captum库（https://captum.ai）提供了计算IG和其他类似指标所需的所有工具。

2.5. 全局重要性指标
积分梯度{IGj}是预测因子的赋值分数——它们告诉每个预测因子对单个数据点（样本/观测值）的预测响应值的贡献程度。对于不同的数据点，变量可能有不同的赋值分数。为了将此信息推广到整个数据集（训练集或测试集），并计算所谓的全局重要性指标GI，可以使用平均绝对赋值[22]对每个数据点获得的个别IG值进行汇总：
(23) GIj = 1/n ∑i=1^n |IGij|

或者，可以计算均方根赋值：
(这是因为基线x′的选择构成了模型判断“重要”特征的参考点。虽然在许多数据集中（例如图像中的黑色像素）零值表示“无信息”，但在光谱学中并非如此。全为零的光谱在物理上往往是不可行的（根据测量方式的不同，它可能意味着所有波长的透射率为100%或反射率为0%）。由于卷积神经网络（CNN）模型从未在这样的状态下进行训练，因此在该基线附近计算的梯度往往带有噪声或在数学上没有意义。这导致某些波长被标记为重要，仅仅是因为它们的值“不为零”，而不是因为它们包含了化学上相关的特征。这个问题可以通过两种主要方法来解决：

1. 在训练之前对光谱进行均值中心化处理。均值中心化后，特定波长的值0代表训练集的平均光谱。使用零向量作为基线可以使积分重要性（IG）分析突出显示某个样本与平均光谱的差异。正的归因值表明偏离平均值使预测结果升高，而负的归因值则表明预测结果降低。
2. 通过从数据集中随机抽取几个真实光谱作为基线。然后对这些多个参考值进行IG计算的平均。尽管这种方法计算要求更高，但它可能提供更稳定的结果，并且对噪声和多重共线性的敏感性较低[24]。

在以下示例中，除非另有说明，否则都使用第一种方法。

3. 引导示例

3.1. 模拟数据
创建了一个简单的模拟数据集，其中包含3个预测变量和1个响应变量，具体步骤如下：
1. 生成训练集的预测变量，x1、x2和x3的范围是从-1到1，每个变量有5个等级。
2. 从U(-1,1)中生成测试集的预测变量，x1、x2和x3的100个随机值。
3. 计算两个集合的响应值，公式为y=b1x1^2+b2x2^2+b3x3^2+N(0,σ^2)。

使用这个训练集来训练一个具有两层的神经网络（结构为3-10-1）：
- 第一层包含10个神经元，激活函数为ReLU（因此该层有3个输入和10个输出）。
- 第二层包含1个神经元，没有激活函数（因此该层有10个输入和1个输出）。
使用随机梯度下降（Stochastic Gradient Descent）作为优化器，学习率为α=0.1，对网络进行1000个时代的训练。训练好的网络用于对训练集和测试集进行预测，并计算积分梯度。相应的Python代码可以在补充材料中找到。

下载：高分辨率图像（353KB）
下载：全尺寸图像

图2. 模拟数据集（蓝色圆圈 - 训练集，橙色十字 - 测试集）。图2显示了两个3D散点图，用于可视化b=[2.0,3.0,0.1]、σ=0.5时的模拟数据。这里提前知道x3对y值并不重要。

下载：高分辨率图像（499KB）
下载：全尺寸图像

图3. 每个观测值的积分梯度（上图）以及相应的平方预测值（下图）。左侧图表显示原始数据；右侧图表中的观测值按y值排序。上图的颜色代表IG值，下图的颜色代表平方预测值。图3包含了训练集中每个数据点的积分梯度值（上图）以及相应的平方预测值（下图）的热图。左侧图表针对原始训练集，右侧图表显示按y值从低到高排序的观测值。

下载：高分辨率图像（184KB）
下载：全尺寸图像

图4. 对四个数据集训练的模型计算的全局重要性指标。这些数据集使用不同的系数b生成，系数b显示在每个图的顶部。如图所示，x1和x2的IG值与其预测值之间存在明显的相关性，而对于x3，无论预测值如何，其IG值都保持较低——模型“知道”x3对预测结果不重要，因此相应的IG值很小。图4显示了使用四个数据集的个别IG值的平均绝对归因计算得出的全局重要性指标，这些数据集使用不同的系数b生成。系数显示在每个图的顶部。可以看出，GI值与系数之间存在直接关联。

下载：高分辨率图像（565KB）
下载：全尺寸图像

图5. 对包含不同噪声水平σ的数据集训练的模型计算的全局重要性指标。图标题显示了相应的σ值以及在所有迭代中计算出的R2值范围。与回归系数一样，只有当模型可靠时，IG值才值得信赖。图5显示了使用相同参数向量b=[2.0,0.1,3.0]但不同的σ值进行30个训练/测试周期后计算出的GI值。相应的σ值显示在每个图的顶部，同时还显示了使用模型对测试集预测计算出的决定系数R2的最小值和最大值。

下载：高分辨率图像（207KB）
下载：全尺寸图像

图6. 对基于两种不同架构的模型计算的全局重要性指标。图标题显示了架构标签以及在所有迭代中计算出的R2值范围。众所周知，由于权重初始化使用了随机值，因此ANN训练无法完全重现，因此在同一数据集上训练的多个模型的参数会有所不同。这种差异的幅度与预测变量和响应变量之间的关系强度直接相关。如图所示，当模型表现良好（σ=0.1）时，GI值稳定，一致地将x1和x3视为最重要的变量，同时保持x2的值较低。尽管它们对响应的贡献差异很小，但在所有情况下x3仍然比x1更受青睐。然而，当噪声增大时，不同训练/测试迭代中的GI值变得不稳定。尽管如此，即使对于σ=2.0（此时测试集的R2值非常低且IG值的变化较大），总体趋势仍然可以区分。

最后，必须指出，假设模型有效，IG值在不同架构上是稳定的。图6显示了使用原始模型（表示为“3-10-1”）和一个过于复杂的模型（结构为“3-10-100-10-1”）进行30个训练/测试周期后得到的GI值。在这两种情况下，数据集都是使用b=[2.0,0.1,3.0]和σ=0.1生成的。如图所示，复杂模型的R2值变化稍大，但GI值在两种情况下都相当一致。

3.2. Tecator数据
Tecator数据集包含了水分、蛋白质和脂肪含量不同的细切肉样本的近红外光谱。这些光谱在850–1050纳米范围内进行了100个波长的测量。测量使用的是基于近红外透射（NIT）原理的Tecator Infratec食品和饲料分析仪。该数据集从StatLib公共数据集档案（http://lib.stat.cmu.edu/datasets/）下载，包含总共215个样本，其中172个作为训练集，43个作为测试集。预测变量矩阵X包含光谱值，而脂肪含量作为响应变量。该数据集中的光谱形状（如图8左上图所示）受到散射、噪声等不同副作用的影响，这使得揭示相关信息变得困难。此外，预测变量与响应变量之间存在轻微的非线性关系，因此推荐进行适当的预处理或非线性建模[25]。

这里考虑了三组预测变量：原始光谱、通过标准正态变量（SNV）方法预处理的光谱，以及通过应用Savitzky-Golay滤波器（滤波器宽度为9，二次多项式和二次导数，然后进行SNV归一化）预处理的光谱。对每组数据分别训练了一个ANN模型。在所有情况下，预测变量都进行了均值中心化（以避免前面提到的零基线问题），同时响应值进行了标准化，尽管本文中的所有相关图表和统计都是基于原始未标准化的响应值显示的。

ANN模型的架构如下：
- 第一层为卷积层，具有一个输入通道、两个输出通道，滤波器宽度为7，填充量为3。该层对原始光谱应用两组不同的滤波器，从而生成两个特征向量（过滤后的光谱），每个向量包含100个值（与原始光谱相同）。
- 第二层为卷积层，具有两个输入通道、四个输出通道，滤波器宽度为3，填充量为1。该层从前一层获取两个信号，并对其各自应用两组不同的滤波器，从而生成四个特征向量，每个向量包含100个值。
- 三个全连接层，配置为400-100-10-1。除了输出层外，每层都使用了ReLU函数进行激活。

与模拟示例不同，近红外光谱包含局部光谱特征，如吸收峰和肩峰，其形状携带了关于样本化学成分的重要信息。卷积层非常适合提取这些特征，因为它们应用了可学习的滤波器，这些滤波器沿着光谱轴滑动并检测相邻波长通道中的局部模式。因此，卷积层作为一种可学习的光谱预处理步骤，类似于传统的Savitzky-Golay平滑或微分方法，但针对预测任务进行了优化。这三个ANN模型分别使用了0.0005的学习率和2000、1000、500个训练周期进行训练——光谱预处理得越充分，所需的训练周期就越少。优化器使用Adam，加上小的权重衰减来施加正则化，从而使模型在处理多重共线性数据集时更加稳定。训练过程在MacBook Air M2上耗时4到20秒。

部分最小二乘（Partial Least Squares）回归被用作基线，以便将积分梯度和相应的全局重要性指数与PLS回归系数进行比较。PLS模型使用SIMPLS算法开发，通过测试集性能确定了最佳组件数量，分别为6、4和3，对应于原始、SNV和Savgol + SNV三组预测变量。

下载：高分辨率图像（670KB）
下载：全尺寸图像

图7. 为每组预测变量（原始、SNV、Savgol + SNV）开发的CNN（上图）和PLS（下图）回归模型的预测值与测量值对比图。图例块内的数值表示均方根误差（RMSE）。

图7显示了预测值与测量值以及每个模型和每组预测变量的均方根误差（RMSE）值（在图例文本中括号中显示）。图8的上图显示了按脂肪含量着色的光谱（蓝色对应低脂肪含量，黄色对应高脂肪含量）。如该数据集所知，原始光谱存在各种失真，因此不能清楚地显示光谱形状与相应脂肪含量之间的关联。SNV预处理去除了最大的变异来源——散射效应——并揭示了这种关联，而使用二次导数后进行SNV归一化使这种关联更加明显，这一点通过图7中的模型性能得到了证实。

下载：高分辨率图像（882KB）
下载：全尺寸图像

图8. 光谱（上排）、积分梯度热图（中排）以及基于IG的全局重要性指数和PLS回归系数的绝对值（下排）针对每组预测变量。光谱图按脂肪含量着色（蓝色对应低脂肪含量，黄色对应高脂肪含量）。热图的颜色 represent 集积分梯度的值。图8第二排的图表显示了相应积分梯度的热图。这些行按脂肪含量排序（左排 - 低含量，右排 - 高含量）。可以看出，在第一种情况（原始数据）中，IG热图强调了三个波段：870至910纳米之间、910至940纳米之间以及950纳米左右。然而，在基于原始光谱训练的模型中，这种趋势并不一致，可以找到一些样本，例如在中间波段的积分梯度值较低，而在左右波段的值较高。在基于SNV预处理光谱的模型IG图中（中间列），中间波段的梯度变得更强，从而减少了来自左右波段的“补偿”效应。这些波段的宽度也有所减小。这些观察结果通过图底部显示的GI图得到了验证。在这两种情况下，大约930纳米处的峰值都变得最强（因此对预测值的贡献最大），但主要峰值两侧的两个峰值仍然明显存在。对于SNV处理过的模型，大约950纳米处的右侧峰值分裂成了两个更窄的峰值。最后一个案例（SavGol + SNV）使这些模式更加清晰。因此，波段变得更窄（这要归功于导数的作用），与响应值的关联也更强。大约950纳米处有两个窄波段而不是一个宽波段的现象也更加明显。还可以注意到，在这种情况下，IG热图的大部分区域呈现为“平坦”的浅绿色，这表明大多数波长都是无关紧要的，只有少数几个窄波段对预测有贡献。有趣的是，GI图与PLS回归系数的绝对值非常相似（在底部图中以灰色曲线显示）。最大的差异出现在SNV处理过的案例中，基于IG的指标显示波段更窄。这很可能是由于卷积层的作用，卷积层本质上实现了基于导数的过滤。

值得一提的是，IG值与预测值有直接的关系——IG矩阵的每一行之和（在图中以热图形式显示）等于预测的响应值，因此IG值直接反映了每个预测因子对模型结果的贡献。如前所述，由于权重初始化是随机的，ANN训练并不完全可复现，因此观察到如果重复训练过程（类似于前一节对模拟数据的研究），IG/GI的变化有多大是很有趣的。图9显示了从 scratch 开始训练的10个模型的GI图。各个GI值用灰色显示，而红色曲线显示了迭代平均值。

图9表明，不同模型的GI值确实有所变化，但它们一致地指出了光谱中相同的重要部分。基于SavGol + SNV处理光谱计算的模型的变化最小，这是预期之中的，因为这个模型具有最佳的测试性能，因此能够最有效地捕捉预测因子和响应之间的关系。

3.3. 啤酒数据
啤酒数据集被广泛认为是噪声和系统性但无关的变化如何妨碍PLS找到预测因子和响应之间正确关系的一个示例[9]。关于该数据集的详细描述可以在其他地方找到。简而言之，它包含了60个不同种类啤酒样品的可见光/近红外光谱。图10显示了包含在训练集中的40个样品的原始光谱和经过均值中心化的光谱（其余20个样品用于测试集），并按响应值进行了着色。

如图所示，样本在酒精含量和颜色上都有所不同，因此光谱的可见光部分（400–700纳米）变化很大。此外，1750纳米及以上的部分既与响应无关，又非常嘈杂。最重要的部分位于1100到1350纳米之间，这个区域在底部图中用浅灰色矩形标出。在这个例子中使用的ANN模型具有以下架构：
1. 一个卷积层，输入通道1个，输出通道2个，滤波器宽度11，填充5。
2. 一个卷积层，输入通道2个，输出通道4个，滤波器宽度7，填充3。
3. 一个卷积层，输入通道4个，输出通道8个，滤波器宽度3，填充1。
4. 三个全连接层，配置为7408-1024-64-1。与Tecator数据集的情况一样，除了最后一个层之外，所有层都使用了ReLU函数进行激活。

这种类型的数据集，其中训练样本的数量（40个）远小于预测因子的数量（926个），通常不被认为适合用于人工神经网络。然而，由于卷积层的作用，可以使用均值中心化的光谱获得一个相对良好的CNN模型，其RMSEP为0.32（相比之下，[9]中报告的相应PLS模型的RMSEP为0.65）。然而，如果观察积分梯度热图和对应的全局重要性指数图（图11中的两个左图），可以看出它们都表明了可见光部分对某些样本预测的重要性。虽然热图右侧部分看起来很嘈杂，但由于全局重要性指数是基于绝对值计算的，因此在GI图上也获得了相对较高的分数。

如果计算IG值与相应标准化响应值的乘积（IGR），然后通过对列取平均值得出全局重要性指数（图11中的两个右图），情况会有所不同。由于无关部分的IG值与响应值不相关，这减少了它们的影响，并在GIR图中清晰地标识出最重要的部分。值得注意的是，在这种情况下，GI图和GIR图类似于为PLS模型创建的VIP分数和选择性比率图[9]。

4. 结论
对于ANN模型（无论是深度模型还是浅层模型），基于梯度的归因方法提供了一种计算效率更高的替代方案，相比模型不可知的 方法（如SHAP）。这些方法利用了神经网络固有的可微性，在短时间内计算出归因分数。其中，积分梯度最为高效，因为它们克服了简单基于梯度的方法（如普通梯度显著性）的关键限制。通过沿着从基线到感兴趣输入（例如光谱）的路径积分梯度，该方法可以避免饱和和零梯度问题，这些问题可能导致误导性的重要性分数。使用模拟数据进行的实验表明，积分梯度能够正确识别预测因子的重要性，并且在模型表现良好时在不同模型架构中保持稳定性。然而，当模型性能因噪声或训练不足而下降时，重要性分数的稳定性和可靠性也会相应降低。这强调了一个基本原则，即像任何基于模型的指标一样，积分梯度依赖于底层模型的质量。

在使用IG处理光谱数据时，正确选择基线非常重要。在训练前对光谱进行均值中心化或使用真实光谱作为基线可以提供更有意义的参考点，并得出更易解释的结果。将积分梯度方法应用于真实的光谱数据集（Tecator和啤酒数据集）表明，它可以成功地识别化学上相关的光谱区域，同时抑制噪声和无关的系统变化。个别积分梯度值的热图揭示了不同样本如何利用光谱信息，而全局重要性指数将这些信息汇总成可解释的汇总统计量。最后，使用梯度-响应乘积为基于梯度的指标提供了额外的信息增强工具。通过结合归因分数与响应值之间的关系，这种方法可以更好地区分相关的光谱特征和噪声及无信息性的系统变化，正如啤酒数据集所示。

总体而言，基于梯度的归因方法为光谱学家和化学计量学家提供了一种强大的工具，用于解释ANN模型。它们可以帮助识别重要的波长区域，验证模型是否使用化学上有意义的特征，并提供关于模型如何处理单个样本的见解。这些方法的计算效率使它们特别适用于光谱分析的常规应用。积分梯度还可以作为一种间接的过拟合诊断工具。如果一个过参数化的模型拟合的是噪声而不是化学上有意义的模式，那么IG重要性剖面预计会显得随机或在光谱轴上均匀分布。相反，本研究中获得的IG剖面突出了与已知吸收带（例如与脂肪含量相关）一致的具体光谱区域，证明了模型学习了具有化学解释能力的特征，而不是偶然的相关性。然而，实践者应该意识到这些方法的局限性。基于梯度的方法是模型特定的，仅对特定的训练模型具有解释性。过拟合或训练不当的模型会产生不可靠的重要性分数。此外，尽管积分梯度满足完整性公理（归因之和等于预测值与基线之间的差异），但这并不意味着所标识的特征是响应变量的真正因果驱动因素——仅仅说明它们对模型预测很重要。