Softmax函数及其导数

Softmax函数及其导数

本文翻译自The Softmax function and its derivative

基础概念

Softmax函数的输入是N维的随机真值向量，输出是另一个N维的真值向量，
且值的范围是$(0,1)$，和为1.0。即映射：$S(\textbf{a})=\mathbb{R}^N\rightarrow \mathbb{R}^N$：

$$\begin{equation*} S(\textbf{a}): \begin{bmatrix} a_1\\ a_2\\ ...\\ a_N \end{bmatrix} \rightarrow \begin{bmatrix} S_1\\ S_2\\ ...\\ S_N \end{bmatrix} \end{equation*}$$
其中每一个元素的公式为：
$$\begin{equation*} S_j=\frac{e^{a_j}}{\sum_{k=1}^{N}e^{a_k}} \qquad \forall j\in 1...N \end{equation*}$$

显然$S_j$总是正的~(因为指数)；因为所有的$S_j$的和为1，所以有$S_j<1$，因此它的范围是$(0,1)$。例如，一个含有三个元素的向量$\left [ 1.0,\,2.0,\,3.0 \right ]$被转化为$\left [ 0.09,\,0.24,\,0.67 \right ]$。
转化后的元素与原始的对应的元素位置上保持一致，且和为1。我们将原始的向量拉伸为$\left [ 1.0,\,2.0,\,5.0 \right ]$，得到变换后的$\left [ 0.02,\,0.05,\,0.93 \right ]$，同样具有前面的性质。注意此时因为最后一个元素(5.0)距离前面两个元素(1.0和2.0)较远，因此它的输出的softmax值占据了和1.0的大部分(0.93)。softmax并不是只选择一个最大元素，而是将向量分解为整个(1.0)的几部分，最大的输入元素得到一个比例较大的部分，但其他元素各自也获得对应的部分。

概率解释

softmax的性质(所有输出的值范围是$(0,1)$且和为1.0)使其在机器学习的概率解释中广泛使用。尤其是在多类别分类任务中，我们总是给输出结果对应的类别附上一个概率，即如果我们的输出类别有N种，我们就输出一个N维的概率向量且和为1.0。每一维的值对应一种类别的概率。我们可以将softmax解释如下：

$$\begin{equation*} S_j=P(y=j|a) \end{equation*}$$
其中， $y$是输出的N个类别中的某个(取值为$1...N$)。 $a$是任意一个N维向量。最常见的例子是多类别的逻辑斯谛回归，输入的向量$x$乘以一个权重矩阵W，且该结果输入softmax函数以产生概率。我们在后面会探讨这个结构。事实证明，从概率的角度来看，softmax对于模型参数的最大似然估计是最优的。 不过，这超出了本文的范围。有关更多详细信息，请参阅“深度学习”一书的第5章(链接：www.deeplearningbook.org)。

向量计算的准备

在深入理解计算softmax的导数之前，我们先了解向量计算的一些基础知识。
Softmax从根本上来说是一种向量函数。它将向量作为输入并输出另一个向量。换言之，它有多个输入和输出，因此我们不能直接就尝试求”softmax的导数”，我们首先要明确：

• 我们想要计算softmax的哪个组成成分(输出的某元素)的导数。
• 由于softmax具有多个输入，所以要计算关于哪个输入元素的偏导数。

听起来好像很复杂，但这正是为什么定义向量计算的原因。 我们正在寻找的偏导数是：

$$\begin{equation*} \frac{\partial S_i}{\partial a_j} \end{equation*}$$
这是第 $i$个输出关于第$j$个输入的偏导数。我们使用一个更简洁的式子来表示： $D_jS_i$。

因为softmax函数是一个$\mathbb{R}^N\rightarrow \mathbb{R}^N$的函数，所以我们计算得到的导数是一个雅可比矩阵：

$$\begin{equation*} DS= \begin{bmatrix} D_1S_1 & \cdots\ & D_NS_1\\ \vdots & \ddots & \vdots \\ D_1S_N & \cdots\ & D_NS_N\\ \end{bmatrix} \end{equation*}$$

在机器学习的文献中，常常用术语梯度来表示通常所说的导数。严格来说，梯度只是为标量函数来定义的，例如机器学习中的损失函数；对于像softmax这样的向量函数，说是“梯度”是不准确的；雅可比是一个向量函数的全部的导数，大多数情况下我们会说“导数”。

softmax的导数

对任意的$i$和$j$，让我们来计算$D_jS_i$：

$$\begin{equation*} D_jS_i=\frac{\partial S_i}{\partial a_j}=\frac{\partial \frac{e^{a_i}}{\sum_{k=1}^{N}e^{a_k}}}{\partial a_j} \end{equation*}$$

我们将使用链式法则来计算导数，即对于$f(x)=\frac{g(x)}{h(x)}$：

$$\begin{equation*} {f}'(x)=\frac{{g}'(x)h(x)-{h}'(x)g(x)}{[h(x)]^2} \end{equation*}$$

在我们的情况下，有：

$$\begin{equation*} \begin{split} g_i=e^{a_i} \\ h_i=\sum_{k=1}^{N}e^{a_k} \end{split} \end{equation*}$$

注意对于$h_i$，无论求其关于哪个$a_j$的导数，结果都是$e^{a_j}$，但是对于$g_i$就不同了。
$g_i$关于$a_j$的导数是$e^{a_j}$当且仅当$i=j$；否则结果为0。

让我们回到$D_jS_i$；我们先考虑$i=j$的情况。根据链式法则我们有：

$$\begin{equation*} \frac{\partial \frac{e^{a_i}}{\sum_{k=1}^{N}e^{a_k}}}{\partial a_j}= \frac{e^{a_i}\sum-e^{a_j}e^{a_i}}{\sum^2 } \end{equation*}$$
简单起见，我们使用 $\sum$表示 $\sum_{k=1}^{N}e^{a_k}$。继续化简下：
$$\begin{equation*} \begin{split} \frac{\partial \frac{e^{a_i}}{\sum_{k=1}^{N}e^{a_k}}}{\partial a_j}= \frac{e^{a_i}\sum-e^{a_j}e^{a_i}}{\sum^2 } \\ =\frac{e^{a_i}}{\sum }\frac{\sum-e^{a_j}}{\sum} \\ =S_i(1-S_j) \end{split} \end{equation*}$$
最后的公式使用其自身来表示 $(S_i$和 $S_j)$，这在包含指数函数时是一个常用的技巧。

类似的，考虑$i\neq j$的情况：

$$\begin{equation*} \begin{split} \frac{\partial \frac{e^{a_i}}{\sum_{k=1}^{N}e^{a_k}}}{\partial a_j}= \frac{0-e^{a_j}e^{a_i}}{\sum^2 } \\ =-\frac{e^{a_j}}{\sum }\frac{e^{a_i}}{\sum} \\ =-S_jS_i \end{split} \end{equation*}$$

总结如下：

$$\begin{equation*} D_jS_i=\left\{\begin{matrix} S_i(1-S_j) & i=j\\ -S_jS_i & i\neq j \end{matrix}\right. \end{equation*}$$

在文献中我们常常会见到各种各样的”浓缩的”公式，一个常见的例子是使用克罗内克函数：

$$\begin{equation*} \delta_{ij}=\left \{\begin{matrix} 1 & i=j\\ 0 & i\neq j \end{matrix}\right. \end{equation*}$$

于是我们有：

$$\begin{equation*} D_jS_i=S_i(\delta_{ij}-S_j) \end{equation*}$$
在文献中也有一些其它的表述：

• 在雅可比矩阵中使用单位矩阵$I$来替换$\delta$，$I$使用元素的矩阵形式表示了$\delta$。
• 使用”1”作为函数名而不是克罗内克$\delta$，如下所示：$D_jS_i=S_i(1(i=j)-S_j)$。这里1(i=j)意味着当$i=j$时值为1，否则为0。

当我们想要计算依赖于softmax导数的更复杂的导数时，“浓缩”符号会很有用; 否则我们必须在任何地方完整的写出公式。

计算softmax和数值稳定性

对于一个给定的向量，使用Python来计算softmax的简单方法是：

def softmax(x):
    """Compute the softmax of vector x."""
    exps = np.exp(x)
    return exps / np.sum(exps)

使用前面定义的softmax函数计算一个三维的向量：

In [146]: softmax([1, 2, 3])
Out[146]: array([ 0.09003057, 0.24472847,  0.66524096])

然而当我们使用该函数计算较大的值时(或者大的负数时)，会出现一个问题：

In [148]: softmax([1000, 2000, 3000])
Out[148]: array([ nan,  nan,  nan])

Numpy使用的浮点数的数值范围是有限的。对于float64，最大可表示数字的大小为$10^{308}$。
softmax函数中的求幂运算可以轻松超过这个数字，即使是相当适中的输入。避免这个问题的一个好方法是通过规范输入使其不要太大或者太小，通过观察我们可以使用任意的常量C，如下所示：

$$\begin{equation*} S_j=\frac{e^{a_j}}{\sum_{k=1}^{N}e^{a_k}}=\frac{Ce^{a_j}}{\sum_{k=1}^{N}Ce^{a_k}} \end{equation*}$$
然后将这个变量转换到指数上：
$$\begin{equation*} S_j=\frac{e^{a_j+log(C)}}{\sum_{k=1}^{N}e^{a_k+log(C)}} \end{equation*}$$
因为C是一个随机的常量，所以我们可以写为：
$$\begin{equation*} S_j=\frac{e^{a_j+D}}{\sum_{k=1}^{N}e^{a_k+D}} \end{equation*}$$
D也是一个任意常量。对任意D，这个公式等价于前面的式子，这让我们能够更好的进行计算。对于D，一个比较好的选择是所有输入的最大值的负数：
$$\begin{equation*} D=-max(a_1,\,a_2,\,...,a_n) \end{equation*}$$
假定输入本身彼此相差不大，这会使输入转换到接近于0的范围。最重要的是，它将所有的输入转换为负数(除最大值外，最大值变为0)。很大的负指数结果会趋于0而不是无穷，这就让我们很好的避免了出现NaN的结果。

def stablesoftmax(x):
    """Compute the softmax of vector x in a numerically
    stable way."""
    shiftx = x - np.max(x)
    exps = np.exp(shiftx)
    return exps / np.sum(exps)

现在我们有：

In [150]: stablesoftmax([1000, 2000, 3000])
Out[150]: array([ 0.,  0.,  1.])

请注意，这仍然是不完美的，因为数学上softmax永远不会真的产生零，但这比NaN好得多，且由于输入之间的距离非常大，所以无论如何都会得到非常接近于零的结果。

softmax层及其导数

softmax常用于机器学习中，特别是逻辑斯特回归：softmax层，其中我们将softmax应用于全连接层(矩阵乘法)的输出，如图所示。

在这个图中，我们有一个具有N个特征的输入x和T个可能的输出类别。权重矩阵W用于将x转换成具有T元素的向量(在机器学习的文献中称为“logits”)，并且softmax函数用于将logits转换成表示属于某一类别的概率。

我们如何计算这个“softmax层”的导数(先进行全连接矩阵乘法，然后是softmax)？当然是使用链式规则！

在我们开始之前的一个重要的观点：你可能会认为x是计算其导数的自然变量(natural variable)。但事实并非如此。实际上，在机器学习中，我们通常希望找到最佳的权重矩阵W，因此我们希望用梯度下降的每一步来更新权重。因此，我们将计算该层的关于W的导数。

我们首先将这个图改写为向量函数的组合。首先我们定义矩阵乘法$g(W)$，即映射：$\mathbb{R}^{NT}\rightarrow \mathbb{R}^T$。因为输入(矩阵W)$N\times T$个元素，输出有T个元素。

接下来我们来考虑softmax，如果我们定义logits的向量是$\lambda$，我们有：$\mathbb{R}^{T}\rightarrow \mathbb{R}^T$。总体来说，我们有：

$$\begin{equation*} P(W)=S(g(W))=(S \circ g)(W) \end{equation*}$$
使用多变量的链式法则，得到 $P(W)$的雅可比矩阵：
$$\begin{equation*} DP(W)=D(S \circ g)(W)=DS(g(W))\cdot Dg(W) \end{equation*}$$

我们之前已经计算过雅可比矩阵；只不过此时是对g(W)求解。因此g是一个非常简单的函数，因此计算雅可比矩阵很简单。唯一要注意的是
正确计算相应的索引。因为$g(W)$：$\mathbb{R}^{NT}\rightarrow \mathbb{R}^T$，所以它的雅可比矩阵是$T$行，$NT$列。

$$\begin{equation*} Dg= \begin{bmatrix} D_1g_1 & \cdots\ & D_{NT}g_1\\ \vdots & \ddots & \vdots \\ D_1g_T & \cdots\ & D_{NT}g_T\\ \end{bmatrix} \end{equation*}$$
在某种意义上，权重矩阵W被“线性化”为长度为NT的向量。 如果您熟悉多维数组的内存布局，应该很容易理解它是如何完成的。
在我们的例子中，我们可以做的一件事就是按照行主次序对其进行线性化处理，第一行是连续的，接着是第二行，等等。 $W_{ij}$
在雅可比矩阵中的列号是 $(i-1)N+j$。为了计算 $Dg$，让我们回顾 $g_1$：
$$\begin{equation*} g_1=W_{11}x_1+W_{12}x_2+\cdots+W_{1N}x_N \end{equation*}$$
因此：
$$\begin{equation*} \begin{split} D_1g_1=x_1 \\ D_2g_1=x_2 \\ \cdots \\ D_Ng_1=x_N \\ D_{N+1}g_1=0 \\ \cdots \\ D_{NT}g_1=0 \end{split} \end{equation*}$$
我们使用同样的策略来计算 $g_2\cdots g_T$，我们可以得到雅可比矩阵：
$$\begin{equation*} Dg= \begin{bmatrix} x_1 &x_1 &\cdots &x_N &\cdots &0 &0 &0 &0 \\ \vdots &\ddots &\ddots &\ddots &\ddots &\ddots &\ddots &\ddots &\vdots \\ 0 &0 &\cdots &0 &\cdots &x_1 &x_2 &\cdots &x_N \end{bmatrix} \end{equation*}$$
最后从另一个角度来这个问题，如果我们将W的索引分解为i和j，可以得到：

$$\begin{equation*} \begin{split} D_{ij}g_t=\frac{\partial (W_{t1}x_1+W_{t2}x_2+\cdots+W_{tN}x_N)}{\partial W_{ij}}\\ = \left\{\begin{matrix} x_j & i=t\\ 0 & i\neq t \end{matrix}\right. \end{split} \end{equation*}$$
在雅可比矩阵中表示第 $t$行，$(i-1)N+j$列。
最后，为了计算softmax层的完整的雅可比矩阵，我们只需要计算 $DS$和 $Dg$间的乘积。注意 $P(W)$： $\mathbb{R}^{NT}\rightarrow \mathbb{R}^T$，因此雅可比矩阵的维度可以确定。因此 $DS$是 $T\times T$， $Dg$是 $T\times NT$的，它们的乘积 $DP$是 $T\times NT$的。
在文献中，你会看到softmax层的导数大大减少了。因为涉及的两个函数很简单而且很常用。 如果我们仔细计算 $DS$的行和 $Dg$的列之间的乘积：
$$\begin{equation*} D_{ij}P_t=\sum_{k=1}^{T}D_kS_t\cdot D_{ij}g_k \end{equation*}$$
$Dg$大多数为0，所以最终的结果很简单，仅当 $i=k$时 $D_{ij}g_k$不为0；然后它等于 $x_j$。因此：
$$\begin{equation*} D_{ij}P_t=D_iS_tx_j=S_t(\delta_{ti}-S_i)x_j \end{equation*}$$
因此完全可以在没有实际雅可比矩阵乘法的情况下计算softmax层的导数; 这很好，因为矩阵乘法很耗时！由于全连接层的雅可比矩阵是稀疏的，我们可以避免大多数计算。

Softmax和交叉熵损失

我们刚刚看到softmax函数如何用作机器学习网络的一部分，以及如何使用多元链式规则计算它的导数。当我们处理这个问题的时候，经常看到损失函数和softmax一起使用来训练网络：交叉熵。

交叉熵有一个有趣的概率和信息理论解释，但在这里我只关注其使用机制。对于两个离散概率分布$p$和$q$，交叉熵函数定义为：

$$\begin{equation*} xent(p,q)=-\sum_{k}^{ }p(k)log(q(k)) \end{equation*}$$
其中 $k$遍历分布定义的随机变量的所有的可能的值。具体而言，在我们的例子中有$T$个输出类别，所以 $k$取值从1到$T$。

如果我们从softmax的输出$P$(一个概率分布)来考量。其它的概率分布是”正确的”类别输出，通常定义为$Y$，是一个大小为$T$的one-hot编码的向量，其中一个元素为1.0(该元素表示正确的类别)，其它为0。让我们重新定义该情况下的交叉熵公式：

$$\begin{equation*} xent(Y,P)=-\sum_{k=1}^{T}Y(k)log(P(k)) \end{equation*}$$
其中$k$遍历所有的输出类别，$P(k)$是模型预测的类别的概率。$Y(k)$是数据集提供的真正的类别概率。我们定义唯一的$Y(k)=1.0$的索引为$y$，因此对所有的$k\neq y$，都有$Y(k)=0$，于是交叉熵公式可以简化为：

$$\begin{equation*} xent(Y,P)=-log(P(y)) \end{equation*}$$
实际上，我们把 $y$当作一个常量，仅使用$P$来表示这个函数。进一步地，因为在我们的例子中 $P$是一个向量，我们可以将$P(y)$表示为 $P$的
第$y$个元素，即 $P_y$：
$$\begin{equation*} xent(P)=-log(P_y) \end{equation*}$$
$xent$的雅可比矩阵是 $1\times T$的矩阵(一个行向量)。因为输出是一个标量且我们有 $T$个输出(向量$P$有 $T$个元素)：
$$\begin{equation*} Dxent(Y,P)=\left [ D_1xent \quad D_2xent \quad \cdots \quad D_Txent\right ] \end{equation*}$$
现在回顾下 $P$可以表示为输入为权值的函数：$P(W)=S(g(W))$。所以我们有另一个函数表示：
$$\begin{equation*} xent(W)=(xent\circ P)(W)=xent(P(W)) \end{equation*}$$
我们可以再次使用多元链式法则来计算 $xent$关于 $W$的梯度：
$$\begin{equation*} \begin{split} Dxent(W)=D(xent\circ P)(W)\\ =Dxent(P(W))\cdot DP(W) \end{split} \end{equation*}$$
我们来检查一下雅可比行矩阵的维数。我们已经计算过了 $DP(W)$，它是 $T\times NT$的。 $Dxent(P(W))$是 $1\times T$的，所以得到的
雅可比矩阵 $Dxent(W)$是 $1\times NT$的。这是有意义的，因为整个网络有一个输出(交叉熵损失，是一个标量)和 $NT$个输入(权重)。
同样的，有一个简单的方式来找到 $Dxent(W)$的简单公式，因为矩阵乘法中的许多元素最终会被消除。注意到 $xent(P)$只依赖于 $P$的
第$y$个元素。因此，在雅可比矩阵中，只有 $D_yxent$是非0的：
$$\begin{equation*} Dxent=\left [ 0 \quad 0 \quad D_yxent \quad \cdots \quad 0 \right ] \end{equation*}$$
其中， $D_yxent=-\frac{1}{P_y}$。回到整个的雅可比矩阵 $Dxent(W)$，使 $Dxent(P)$乘以 $D(P(W))$的每一列，得到结果的行向量的每一个
元素。回顾用行向量表示的按行优先的“线性化”的整个权重矩阵W。清晰起见，我们将使用 $i$和$j$来索引它( $D_{ij}$)表示行向量的中的第
$iN+j$个元素)：
$$\begin{equation*} D_{ij}xent(W)=\sum_{k=1}^{T}D_kxent(P)\cdot D_{ij}P_k(W) \end{equation*}$$
因为在 $D_kxent(P)$中只有第 $y$个元素是非0的，所以我们可以得到下式：
$$\begin{equation*} \begin{split} D_{ij}xent(W)=D_yxent(P)\cdot D_{ij}P_y(W) \\ =-\frac{1}{P_y}\cdot S_y(\delta_{yi}-S_i)x_j \end{split} \end{equation*}$$
根据我们的定义， $P_y=S_y$，所以可得：
$$\begin{equation*} \begin{split} D_{ij}xent(W)=-\frac{1}{S_y}\cdot S_y(\delta_{yi}-S_i)x_j \\ =-(\delta_{yi}-S_i)x_j \\ =(S_i-\delta_{yi})x_j \end{split} \end{equation*}$$
即使最终的结果很简洁清楚，但是我们不一定非要这样做。公式 $D_{ij}xent(W)$可能最终成为一个和的形式(或者某些和的和)。关于雅可比矩阵的这些技巧可能并没有太大意义，因为计算机可以完成所有的工作。我们需要做的就是计算出单个的雅矩阵，这通常毕竟容易，因为它们是更简单的非复合函数。这技术体现了多元链式法则的美妙和实用性。