Softmax函数及其导数

本文翻译自The Softmax function and its derivative

基础概念

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

S(a):⎡⎣⎢⎢⎢a1a2...aN⎤⎦⎥⎥⎥→⎡⎣⎢⎢⎢S1S2...SN⎤⎦⎥⎥⎥

\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*}
其中每一个元素的公式为：

Sj=eaj∑Nk=1eak∀j∈1...N

\begin{equation*} S_j=\frac{e^{a_j}}{\sum_{k=1}^{N}e^{a_k}} \qquad \forall j\in 1...N \end{equation*}

显然SjS_j总是正的~(因为指数)；因为所有的SjS_j的和为1，所以有Sj<1S_j，因此它的范围是(0,1)(0,1)。例如，一个含有三个元素的向量[1.0,2.0,3.0]\left [ 1.0,\,2.0,\,3.0 \right ]被转化为[0.09,0.24,0.67]\left [ 0.09,\,0.24,\,0.67 \right ]。
转化后的元素与原始的对应的元素位置上保持一致，且和为1。我们将原始的向量拉伸为[1.0,2.0,5.0]\left [ 1.0,\,2.0,\,5.0 \right ]，得到变换后的[0.02,0.05,0.93]\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)(0,1)且和为1.0)使其在机器学习的概率解释中广泛使用。尤其是在多类别分类任务中，我们总是给输出结果对应的类别附上一个概率，即如果我们的输出类别有N种，我们就输出一个N维的概率向量且和为1.0。每一维的值对应一种类别的概率。我们可以将softmax解释如下：

Sj=P(y=j|a)

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

向量计算的准备

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

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

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

∂Si∂aj

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

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

DS=⎡⎣⎢⎢D1S1⋮D1SN⋯ ⋱⋯ DNS1⋮DNSN⎤⎦⎥⎥

\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的导数

对任意的ii和jj，让我们来计算DjSiD_jS_i：

DjSi=∂Si∂aj=∂eai∑Nk=1eak∂aj

\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)=g(x)h(x)f(x)=\frac{g(x)}{h(x)}：

f′(x)=g′(x)h(x)−h′(x)g(x)[h(x)]2

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

在我们的情况下，有：

gi=eaihi=∑k=1Neak

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

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

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

∂eai∑Nk=1eak∂aj=eai∑−eajeai∑2

\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表示 ∑Nk=1eak\sum_{k=1}^{N}e^{a_k}。继续化简下：

∂eai∑Nk=1eak∂aj=eai∑−eajeai∑2=eai∑∑−eaj∑=Si(1−Sj)

\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*}
最后的公式使用其自身来表示 (Si(S_i和 Sj)S_j)，这在包含指数函数时是一个常用的技巧。

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

∂eai∑Nk=1eak∂aj=0−eajeai∑2=−eaj∑eai∑=−SjSi

\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*}

总结如下：

DjSi={Si(1−Sj)−SjSii=ji≠j

\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*}

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

δij={10i=ji≠j

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

于是我们有：

DjSi=Si(δij−Sj)

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

在雅可比矩阵中使用单位矩阵II来替换δ\delta，II使用元素的矩阵形式表示了δ\delta。
使用”1”作为函数名而不是克罗内克δ\delta，如下所示：DjSi=Si(1(i=j)−Sj)D_jS_i=S_i(1(i=j)-S_j)。这里1(i=j)意味着当i=ji=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，最大可表示数字的大小为1030810^{308}。
softmax函数中的求幂运算可以轻松超过这个数字，即使是相当适中的输入。避免这个问题的一个好方法是通过规范输入使其不要太大或者太小，通过观察我们可以使用任意的常量C，如下所示：

Sj=eaj∑Nk=1eak=Ceaj∑Nk=1Ceak

\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*}
然后将这个变量转换到指数上：

Sj=eaj+log(C)∑Nk=1eak+log(C)

\begin{equation*} S_j=\frac{e^{a_j+log(C)}}{\sum_{k=1}^{N}e^{a_k+log(C)}} \end{equation*}
因为C是一个随机的常量，所以我们可以写为：

Sj=eaj+D∑Nk=1eak+D

\begin{equation*} S_j=\frac{e^{a_j+D}}{\sum_{k=1}^{N}e^{a_k+D}} \end{equation*}
D也是一个任意常量。对任意D，这个公式等价于前面的式子，这让我们能够更好的进行计算。对于D，一个比较好的选择是所有输入的最大值的负数：

D=−max(a1,a2,...,an)

\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 numericallystable 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)g(W)，即映射：RNT→RT\mathbb{R}^{NT}\rightarrow \mathbb{R}^T。因为输入(矩阵W)N×TN\times T个元素，输出有T个元素。

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

P(W)=S(g(W))=(S∘g)(W)

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

DP(W)=D(S∘g)(W)=DS(g(W))⋅Dg(W)

\begin{equation*} DP(W)=D(S \circ g)(W)=DS(g(W))\cdot Dg(W) \end{equation*}

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

Dg=⎡⎣⎢⎢D1g1⋮D1gT⋯ ⋱⋯ DNTg1⋮DNTgT⎤⎦⎥⎥

\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的向量。如果您熟悉多维数组的内存布局，应该很容易理解它是如何完成的。
在我们的例子中，我们可以做的一件事就是按照行主次序对其进行线性化处理，第一行是连续的，接着是第二行，等等。 WijW_{ij}
在雅可比矩阵中的列号是 (i−1)N+j(i-1)N+j。为了计算 DgDg，让我们回顾 g1g_1：

g1=W11x1+W12x2+⋯+W1NxN

\begin{equation*} g_1=W_{11}x_1+W_{12}x_2+\cdots+W_{1N}x_N \end{equation*}
因此：

D1g1=x1D2g1=x2⋯DNg1=xNDN+1g1=0⋯DNTg1=0

\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*}
我们使用同样的策略来计算 g2⋯gTg_2\cdots g_T，我们可以得到雅可比矩阵：

Dg=⎡⎣⎢⎢x1⋮0x1⋱0⋯⋱⋯xN⋱0⋯⋱⋯0⋱x10⋱x20⋱⋯0⋮xN⎤⎦⎥⎥

\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，可以得到：

Dijgt=∂(Wt1x1+Wt2x2+⋯+WtNxN)∂Wij={xj0i=ti≠t

\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*}
在雅可比矩阵中表示第 tt行，(i−1)N+j(i-1)N+j列。
最后，为了计算softmax层的完整的雅可比矩阵，我们只需要计算 DSDS和 DgDg间的乘积。注意 P(W)P(W)： RNT→RT\mathbb{R}^{NT}\rightarrow \mathbb{R}^T，因此雅可比矩阵的维度可以确定。因此 DSDS是 T×TT\times T， DgDg是 T×NTT\times NT的，它们的乘积 DPDP是 T×NTT\times NT的。
在文献中，你会看到softmax层的导数大大减少了。因为涉及的两个函数很简单而且很常用。如果我们仔细计算 DSDS的行和 DgDg的列之间的乘积：

DijPt=∑k=1TDkSt⋅Dijgk

\begin{equation*} D_{ij}P_t=\sum_{k=1}^{T}D_kS_t\cdot D_{ij}g_k \end{equation*}
DgDg大多数为0，所以最终的结果很简单，仅当 i=ki=k时 DijgkD_{ij}g_k不为0；然后它等于 xjx_j。因此：

DijPt=DiStxj=St(δti−Si)xj

\begin{equation*} D_{ij}P_t=D_iS_tx_j=S_t(\delta_{ti}-S_i)x_j \end{equation*}
因此完全可以在没有实际雅可比矩阵乘法的情况下计算softmax层的导数; 这很好，因为矩阵乘法很耗时！由于全连接层的雅可比矩阵是稀疏的，我们可以避免大多数计算。

Softmax和交叉熵损失

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

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

xent(p,q)=−∑kp(k)log(q(k))

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

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

xent(Y,P)=−∑k=1TY(k)log(P(k))

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

xent(Y,P)=−log(P(y))

\begin{equation*} xent(Y,P)=-log(P(y)) \end{equation*}
实际上，我们把 yy当作一个常量，仅使用PP来表示这个函数。进一步地，因为在我们的例子中 PP是一个向量，我们可以将P(y)P(y)表示为 PP的
第yy个元素，即 PyP_y：

xent(P)=−log(Py)

\begin{equation*} xent(P)=-log(P_y) \end{equation*}
xentxent的雅可比矩阵是 1×T1\times T的矩阵(一个行向量)。因为输出是一个标量且我们有 TT个输出(向量PP有 TT个元素)：

Dxent(Y,P)=[D1xentD2xent⋯DTxent]

\begin{equation*} Dxent(Y,P)=\left [ D_1xent \quad D_2xent \quad \cdots \quad D_Txent\right ] \end{equation*}
现在回顾下 PP可以表示为输入为权值的函数：P(W)=S(g(W))P(W)=S(g(W))。所以我们有另一个函数表示：

xent(W)=(xent∘P)(W)=xent(P(W))

\begin{equation*} xent(W)=(xent\circ P)(W)=xent(P(W)) \end{equation*}
我们可以再次使用多元链式法则来计算 xentxent关于 WW的梯度：

Dxent(W)=D(xent∘P)(W)=Dxent(P(W))⋅DP(W)

\begin{equation*} \begin{split} Dxent(W)=D(xent\circ P)(W)\\ =Dxent(P(W))\cdot DP(W) \end{split} \end{equation*}
我们来检查一下雅可比行矩阵的维数。我们已经计算过了 DP(W)DP(W)，它是 T×NTT\times NT的。 Dxent(P(W))Dxent(P(W))是 1×T1\times T的，所以得到的
雅可比矩阵 Dxent(W)Dxent(W)是 1×NT1\times NT的。这是有意义的，因为整个网络有一个输出(交叉熵损失，是一个标量)和 NTNT个输入(权重)。
同样的，有一个简单的方式来找到 Dxent(W)Dxent(W)的简单公式，因为矩阵乘法中的许多元素最终会被消除。注意到 xent(P)xent(P)只依赖于 PP的
第yy个元素。因此，在雅可比矩阵中，只有 DyxentD_yxent是非0的：

Dxent=[00Dyxent⋯0]

\begin{equation*} Dxent=\left [ 0 \quad 0 \quad D_yxent \quad \cdots \quad 0 \right ] \end{equation*}
其中， Dyxent=−1PyD_yxent=-\frac{1}{P_y}。回到整个的雅可比矩阵 Dxent(W)Dxent(W)，使 Dxent(P)Dxent(P)乘以 D(P(W))D(P(W))的每一列，得到结果的行向量的每一个
元素。回顾用行向量表示的按行优先的“线性化”的整个权重矩阵W。清晰起见，我们将使用 ii和jj来索引它( DijD_{ij})表示行向量的中的第
iN+jiN+j个元素)：

Dijxent(W)=∑k=1TDkxent(P)⋅DijPk(W)

\begin{equation*} D_{ij}xent(W)=\sum_{k=1}^{T}D_kxent(P)\cdot D_{ij}P_k(W) \end{equation*}
因为在 Dkxent(P)D_kxent(P)中只有第 yy个元素是非0的，所以我们可以得到下式：

Dijxent(W)=Dyxent(P)⋅DijPy(W)=−1Py⋅Sy(δyi−Si)xj

\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*}
根据我们的定义， Py=SyP_y=S_y，所以可得：

Dijxent(W)=−1Sy⋅Sy(δyi−Si)xj=−(δyi−Si)xj=(Si−δyi)xj

\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*}
即使最终的结果很简洁清楚，但是我们不一定非要这样做。公式 Dijxent(W)D_{ij}xent(W)可能最终成为一个和的形式(或者某些和的和)。关于雅可比矩阵的这些技巧可能并没有太大意义，因为计算机可以完成所有的工作。我们需要做的就是计算出单个的雅矩阵，这通常毕竟容易，因为它们是更简单的非复合函数。这技术体现了多元链式法则的美妙和实用性。