平时我们基本用pytorch或者tensorflow框架时，基本对特别底层的函数实现关注不多，仅限于知道公式的原理。但是很多大佬往往自己会实现一些源码(比如ListNet复现)，在看这些源码时，经常出现各种有点难以理解的代码，本来很简单的东西，莫名其妙的各种转换，化简完之后可能感觉是一样的，这么费劲周折的折腾啥？殊不知还是对底层的实现原理了解少了，虽然有些源码不需要我们从底层造轮子(完全从底层造轮子也影响效率)，但是能理解其原理在我们debug以及看一些源码时不至于太多疑惑(毕竟国外很多大佬都喜欢实现一些底层utils)。

今天我们来重新认识一下我们经常用的Softmax、LogSumExp和Sigmoid

1. 背景概要

我们知道编程语言中的数值都有一个表示范围的，如果数值过大，超过最大的范围，就是上溢；如果过小，超过最小的范围，就是下溢。
今天要讨论的Softmax、LogSumExp和Sigmoid，就面临着上述溢出的问题，下面的一些梳理也主要用来解决计算Softmax或CrossEntropy时出现的上溢(overflow)或下溢(underflow)问题。

2. Softmax

在机器学习中，计算概率输出基本都需要经过Softmax函数，它的公式应该很熟悉了吧
$$\text{Softmax}(x_i)=\frac{\exp (x_i)}{{\sum }_{j=1}^n\exp (x_j)}\text{\hspace{1em}(1)}$$
但是Softmax存在上溢和下溢大问题。如果$x_i$太大，对应的指数函数也非常大，此时很容易就溢出，得到nan结果；如果$x_i$太小，或者说负的太多，就会导致出现下溢而变成0，如果分母变成0，就会出现除0的结果。
此时我们经常看到一个常见的做法是(其实用到的是指数归一化技巧, exp-normalize)，先计算x中的最大值$b={\max }_{i=1}^n{x}_i$，然后根据
$$\begin{array}{rl}\text{Softmax}(x_i)& =\frac{\exp (x_i)}{{\sum }_{j=1}^n\exp (x_j)}\\ & =\frac{\exp (x_i-b)\cdot \exp (b)}{{\sum }_{j=1}^n\left(\exp (x_j-b)\cdot \exp (b)\right)}\\ & =\frac{\exp (x_i-b)\cdot \exp (b)}{\exp (b)\cdot {\sum }_{j=1}^n\exp (x_j-b)}\\ & =\frac{\exp (x_i-b)}{{\sum }_{j=1}^n\exp (x_j-b)}\\ & =\text{Softmax}(x_i-b)\end{array}$$
这种转换是等价的，经过这一变换，就避免了上溢，最大值变成了$\exp (0)=1$；同时分母中也会有一个1，就避免了下溢。
我们通过实例来理解一下。

def bad_softmax(x):y = np.exp(x)return y / y.sum()x = np.array([1, -10, 1000])
print(bad_softmax(x)) #运行结果
#... RuntimeWarning: overflow encountered in exp
#... RuntimeWarning: invalid value encountered in true_divide
#array([ 0.,  0., nan])

接下来进行上面的优化，并进行测试：

def softmax(x):b = x.max()y = np.exp(x - b)return y / y.sum()print(softmax(x))
# array([0., 0., 1.])
x = np.array([-800, -1000, -1000])
print(bad_softmax(x))
# array([nan, nan, nan])
print(softmax(x))
# array([1.00000000e+00, 3.72007598e-44, 3.72007598e-44])

关于softmax的另外实现，参加下文。

3. LogSumExp

什么是LSE?
LSE被定义为参数指数之和的对数：
$$\text{LSE}(x_1,\cdots ,x_n)=\log \sum _{i=1}^n\exp (x_i)=\log \left(\exp (x_1)+\cdots +\exp (x_n)\right)$$
输入可以看成是一个n维的向量，输出是一个标量。

什么场景下要到LogSumExp呢？
交叉熵loss大家肯定不陌生，其中就有一项是log$p$，比如多分类场景里面就需要在对softmax的结果做对数处理：
$$\begin{array}{rl}\log \left(\text{Softmax}(x_i)\right)& =\log \frac{\exp (x_i)}{{\sum }_{j=1}^n\exp (x_j)}\\ & =x_i-\log \sum _{j=1}^n\exp (x_j)\end{array}$$

能看出来不，后面这一项就是我们上面说的LSE。

理清前因后果之后，我们就来看LSE的实现问题。

因为上面最后一项也有上溢的问题，所以应用跟softmax同样的技巧，得
$$\log \sum _{j=1}^n\exp (x_j)=\log \sum _{j=1}^n\exp (x_j-b)\exp (b)=b+\log \sum _{j=1}^n\exp (x_j-b)\text{\hspace{1em}(4)}$$

$b$同样是取$x$中的最大值。

这样，我们就得到了LSE的最终表示：
$$\text{LSE}(x)=b+\log \sum _{j=1}^n\exp (x_j-b)\text{\hspace{1em}(5)}$$
对LogSumExp求导就得到了exp-normalize(Softmax)的形式，
$$\frac{\partial \left(b+\log {\sum }_{j=1}^n\exp (x_j-b)\right)}{\partial x_j}=\frac{\exp (x_i-b)}{{\sum }_{j=1}^n\exp (x_j-b)}\text{\hspace{1em}(6)}$$

怎么实现LSE

def logsumexp(x):b = x.max()return b + np.log(np.sum(np.exp(x - b)))def softmax_lse(x):return np.exp(x - logsumexp(x))x1 = np.array([1, -10, 1000])
x2 = np.array([-900, -1000, -1000])
print(softmax_lse(x1))
# array([0., 0., 1.])
print(softmax(x1))
# array([0., 0., 1.])
print(softmax_lse(x2))
# array([1.00000000e+00, 3.72007598e-44, 3.72007598e-44])
print(softmax(x2))
# array([1.00000000e+00, 3.72007598e-44, 3.72007598e-44])

---

到此LSE的实现就清晰了，不过返回来再看Softmax，Softmax也可以这样表示:
$$\begin{array}{rl}\text{Softmax}(x_i)& =\text{exp}(\text{log}(\text{Softmax}(x_i)))\\ & =\text{exp}(x_i-\log \sum _{j=1}^n\exp (x_j))\\ & =\exp \left(x_i-b-\log \sum _{j=1}^n\exp (x_j-b)\right)\end{array}\text{\hspace{1em}(7)}$$

那我们是使用exp-normalize还是使用LogSumExp呢？

如果你需要保留Log空间，那么就计算$\log (\text{Softmax})$，此时使用LogSumExp技巧；如果你只需要计算Softmax，那么就使用exp-normalize技巧。

4. Sigmoid

我们知道Sigmoid函数是交叉熵之前常用的激活函数，公式为：
$$\sigma (x)=\frac{1}{1+\exp (-x)}\text{\hspace{1em}(8)}$$

对应的图像如下：

其中包含一个$\exp (-x)$，我们看一下$e^x$的图像：

从上图可以看出，如果x很大，$e^x$会非常大，而很小就没事，变成无限接近0。

当Sigmoid函数中的x负的特别多，那么$\exp (-x)$就会变成$\infty$​，就出现了上溢；

那么如何解决这个问题呢？$\sigma (x)$可以表示成两种形式：
$$\sigma (x)=\frac{1}{1+\exp (-x)}=\frac{\exp (x)}{1+\exp (x)}\text{\hspace{1em}(9)}$$

当$x\ge 0$时，我们根据$e^x$的图像，我们取$\frac{1}{1+\exp (-x)}$的形式；

当 x < 0 时，我们取$\frac{\exp (x)}{1+\exp (x)}$

# 原来的做法
def sigmoid_naive(x):return 1 / (1 + math.exp(-x))# 优化后的做法
def sigmoid(x):if x < 0:return math.exp(x) / (1 + math.exp(x))else:return 1 / (1 + math.exp(-x))print(sigmoid_naive(2000))
# 1.0
print(sigmoid(2000))
# 1.0
print(sigmoid_naive(-2000))
# OverflowError: math range error
print(sigmoid(-2000))
# 0.0

参考链接：

• The Log-Sum-Exp Trick
• 一文弄懂交叉熵损失
• Exp-normalize trick
• 一文弄懂LogSumExp技巧