1. 结合Logistic Regression 分析 Softmax
训练集由
m
m
m 个已标记的样本构成：
{
(
x
(
1
)
,
y
(
1
)
)
,
.
.
.
,
(
x
(
m
)
,
y
(
m
)
)
}
\{(x^{(1)}, y^{(1)}), ... , (x^{(m)}, y^{(m)})\}
{(x(1),y(1)),...,(x(m),y(m))}，其中输入的第
i
i
i 个样例
x
(
i
)
∈
ℜ
n
+
1
x^{(i)} \in \Re^{n+1}
x(i)∈ℜn+1，即特征向量
x
x
x 的维度为
n
+
1
n + 1
n+1 ，其中
x
0
=
1
x_0 = 1
x0=1 对应截距项。由于 logistic 回归是针对二分类问题的，因此类标记
y
(
i
)
∈
{
0
,
1
}
y^{(i)} \in \{0,1\}
y(i)∈{0,1}。假设函数(hypothesis function) 如下：
h
θ
(
x
)
=
1
1
+
e
(
−
θ
T
x
)
h_\theta(x) = \frac{1}{1+e^{(-\theta^Tx)}}
hθ(x)=1+e(−θTx)1
我们将训练模型参数
θ
\theta
θ，使其能够最小化代价函数 ：
J
(
θ
)
=
−
1
m
[
∑
i
=
1
m
y
(
i
)
log
h
θ
(
x
(
i
)
)
+
(
1
−
y
(
i
)
)
log
(
1
−
h
θ
(
x
(
i
)
)
)
]
J(\theta) = -\frac{1}{m} \left[ \sum_{i=1}^m y^{(i)} \log h_\theta(x^{(i)}) + (1-y^{(i)}) \log (1-h_\theta(x^{(i)})) \right]
J(θ)=−m1[i=1∑my(i)loghθ(x(i))+(1−y(i))log(1−hθ(x(i)))]
而在 Softmax回归中，我们解决的是多分类问题（相对于 logistic 回归解决的二分类问题），类标
y
\textstyle y
y 可以取
k
\textstyle k
k 个不同的值（而不是 2 个）。因此，对于训练集
{
(
x
(
1
)
,
y
(
1
)
)
,
…
,
(
x
(
m
)
,
y
(
m
)
)
}
\{ (x^{(1)}, y^{(1)}), \ldots, (x^{(m)}, y^{(m)}) \}
{(x(1),y(1)),…,(x(m),y(m))}，我们有
y
(
i
)
∈
{
1
,
2
,
…
,
k
}
y^{(i)} \in \{1, 2, \ldots, k\}
y(i)∈{1,2,…,k}。（注意此处的类别下标从 1 开始，而不是 0）。例如，在 MNIST 数字识别任务中，我们有
k
=
10
\textstyle k=10
k=10 个不同的类别。
对于输入的每一个样例
x
(
i
)
∈
ℜ
n
+
1
x^{(i)} \in \Re^{n+1}
x(i)∈ℜn+1 ，输出其属于每一类别
j
j
j 的概率值
p
(
y
=
j
∣
x
)
\textstyle p(y=j | x)
p(y=j∣x)，即输出一个
k
k
k 维向量（向量元素的和为1）来表示这
k
k
k 个估计的概率值。即 Softmax 的假设函数为
h
θ
(
x
(
i
)
)
h_{\theta}(x^{(i)})
hθ(x(i))
h
θ
(
x
(
i
)
)
=
[
p
(
y
(
i
)
=
1
∣
x
(
i
)
;
θ
)
p
(
y
(
i
)
=
2
∣
x
(
i
)
;
θ
)
⋮
p
(
y
(
i
)
=
k
∣
x
(
i
)
;
θ
)
]
=
1
∑
j
=
1
k
e
θ
j
T
x
(
i
)
[
e
θ
1
T
x
(
i
)
e
θ
2
T
x
(
i
)
⋮
e
θ
k
T
x
(
i
)
]
h_\theta(x^{(i)}) = \begin{bmatrix} p(y^{(i)} = 1 | x^{(i)}; \theta) \\ p(y^{(i)} = 2 | x^{(i)}; \theta) \\ \vdots \\ p(y^{(i)} = k | x^{(i)}; \theta) \end{bmatrix} =\frac{1}{ \sum_{j=1}^{k}{e^{ \theta_j^T x^{(i)} }} } \begin{bmatrix} e^{ \theta_1^T x^{(i)} } \\ e^{ \theta_2^T x^{(i)} } \\ \vdots \\ e^{ \theta_k^T x^{(i)} } \\ \end{bmatrix}
hθ(x(i))=⎣⎢⎢⎢⎡p(y(i)=1∣x(i);θ)p(y(i)=2∣x(i);θ)⋮p(y(i)=k∣x(i);θ)⎦⎥⎥⎥⎤=∑j=1keθjTx(i)1⎣⎢⎢⎢⎢⎡eθ1Tx(i)eθ2Tx(i)⋮eθkTx(i)⎦⎥⎥⎥⎥⎤
其中
θ
1
,
θ
2
,
…
,
θ
k
∈
ℜ
n
+
1
\theta_1, \theta_2, \ldots, \theta_k \in \Re^{n+1}
θ1,θ2,…,θk∈ℜn+1 是模型的参数。请注意
1
∑
j
=
1
k
e
θ
j
T
x
(
i
)
\frac{1}{ \sum_{j=1}^{k}{e^{ \theta_j^T x^{(i)} }} }
∑j=1keθjTx(i)1 这一项对概率分布进行归一化，使得所有概率之和为 1
其代价函数为：
J
(
θ
)
=
−
1
m
[
∑
i
=
1
m
∑
j
=
1
k
1
{
y
(
i
)
=
j
}
log
e
θ
j
T
x
(
i
)
∑
l
=
1
k
e
θ
l
T
x
(
i
)
]
J(\theta) = - \frac{1}{m} \left[ \sum_{i=1}^{m} \sum_{j=1}^{k} 1\left\{y^{(i)} = j\right\} \log \frac{e^{\theta_j^T x^{(i)}}}{\sum_{l=1}^k e^{ \theta_l^T x^{(i)} }}\right]
J(θ)=−m1[i=1∑mj=1∑k1{y(i)=j}log∑l=1keθlTx(i)eθjTx(i)]
其中
1
{
⋅
}
1\{\cdot\}
1{⋅} 为指示函数：
1
{
⋅
}
=
{
1
{
表
达
式
值
为
真
}
=
1
1
{
表
达
式
值
为
假
}
=
0
1\{\cdot\} = \begin{cases} 1\{表达式值为真\} = 1 \\ 1\{表达式值为假\} = 0 \\ \end{cases}
1{⋅}={1{表达式值为真}=11{表达式值为假}=0
梯度下降公式：
∇
θ
j
J
(
θ
)
=
−
1
m
∑
i
=
1
m
[
x
(
i
)
(
1
{
y
(
i
)
=
j
}
−
p
(
y
(
i
)
=
j
∣
x
(
i
)
;
θ
)
)
]
\nabla_{\theta_j} J(\theta) = - \frac{1}{m} \sum_{i=1}^{m}{ \left[ x^{(i)} \left( 1\{ y^{(i)} = j\} - p(y^{(i)} = j | x^{(i)}; \theta) \right) \right] }
∇θjJ(θ)=−m1i=1∑m[x(i)(1{y(i)=j}−p(y(i)=j∣x(i);θ))]
2. Softmax
结合上面理解此时的 Softmax函数：
h θ ( x ) = σ ( θ T x ) = σ ( z ) = ( σ 1 ( z ) , … , σ m ( z ) ) h_{\theta}(x) = \sigma(\theta^T x) = \sigma(z) = \left({\color{Red}{\sigma{1}(z)}}, \ldots, \sigma_{m}(z)\right) hθ(x)=σ(θTx)=σ(z)=(σ1(z),…,σm(z))
此时的 x x x 为一个样例输入 ∈ ℜ n + 1 \in \Re^{n+1} ∈ℜn+1，等价与上面的一个 x ( i ) x^{(i)} x(i)，其中共有 m m m 个类别， z i z_i zi 表示输入样例 x x x 是第 i i i 个类别的线性预测结果，等价于 z i = θ i T x z_i = \theta^T_i x zi=θiTx
Softmax 函数
σ
(
z
)
=
(
σ
1
(
z
)
,
…
,
σ
m
(
z
)
)
\sigma(z)=\left({\color{Red}{\sigma_{1}(z)}}, \ldots, \sigma_{m}(z)\right)
σ(z)=(σ1(z),…,σm(z)) 定义如下：
o
i
=
σ
i
(
z
)
=
exp
(
z
i
)
∑
j
=
1
m
exp
(
z
j
)
,
i
=
1
,
…
,
m
【
观
察
到
的
数
据
x
(
或
z
)
属
于
类
别
i
的
概
率
，
或
者
称
作
似
然
(
L
i
k
e
l
i
h
o
o
d
)
】
{\color{Green}{o_i}} = {\color{Red}{\sigma_{i}(z)}} =\frac{\exp \left(z_{i}\right)}{\sum_{j=1}^{m} \exp \left(z_{j}\right)}, \quad i=1, \ldots, m \\ {\color{Green}{【观察到的数据 \ x \ (或z)\ 属于类别\ i\ 的概率，或者称作似然 (Likelihood)】}}
oi=σi(z)=∑j=1mexp(zj)exp(zi),i=1,…,m【观察到的数据 x (或z) 属于类别 i 的概率，或者称作似然(Likelihood)】
它在 Logistic Regression 里其到的作用是讲线性预测值转化为类别概率：
m
m
m 代表类别数，假设
z
i
=
w
i
T
x
+
b
i
z_i = w_i^Tx + b_i
zi=wiTx+bi 是第
i
i
i 个类别的线性预测结果，带入
S
o
f
t
m
a
x
Softmax
Softmax 的结果其实就是先对每一个
z
i
z_i
zi 取 exponential 变成非负，然后除以所有项之和进行归一化，现在每个
o
i
=
σ
i
(
z
)
o_{i}=\sigma_{i}(z)
oi=σi(z) 就可以解释成：观察到的数据
x
x
x 属于类别
i
i
i 的概率，或者称作似然 (Likelihood)。
然后 Logistic Regression 的目标函数是根据最大似然原则来建立的，假设数据
x
x
x 所对应的类别为
y
y
y，则根据我们刚才的计算最大似然就是要最大化
o
y
o_y
oy 的值 (通常是使用 negative log-likelihood 而不是 likelihood，也就是说最小化
−
l
o
g
(
o
y
)
-log(o_y)
−log(oy) 的值，这两者结果在数学上是等价的)。后面这个操作就是 caffe 文档里说的 Multinomial Logistic Loss，具体写出来是这个样子：
ℓ
(
y
,
o
)
=
−
log
(
o
y
)
\ell(y, o)=-\log \left(o_{y}\right)
ℓ(y,o)=−log(oy)
3. Softmax-Loss = Softmax + Multinomial Logistic Loss
而 Softmax-Loss 其实就是把两者结合到一起，只要把
o
y
o_y
oy 的定义展开即可：
ℓ
~
(
y
,
z
)
=
−
log
(
e
z
y
∑
j
=
1
m
e
z
j
)
=
log
(
∑
j
=
1
m
e
z
j
)
−
z
y
【
将
观
察
到
的
数
据
x
(
或
z
)
属
于
类
别
y
的
概
率
做
最
大
似
然
估
计
，
或
最
小
n
e
g
a
t
i
v
e
l
o
g
似
然
】
\tilde{\ell}(y, z)=-\log \left(\frac{e^{z_{y}}}{\sum_{j=1}^{m} e^{z_{j}}}\right)=\log \left(\sum_{j=1}^{m} e^{z_{j}}\right)-z_{y} \\ {\color{Green}{\small{【将观察到的数据 \ x \ (或z)\ 属于类别 \ y \ 的概率做最大似然估计，或最小\ negative\ log\ 似然】}}} \\
ℓ~(y,z)=−log(∑j=1mezjezy)=log(j=1∑mezj)−zy【将观察到的数据 x (或z) 属于类别 y 的概率做最大似然估计，或最小 negative log 似然】
比如如果我们要写一个 Logistic Regression 的 solver，那么因为要处理的就是这个东西，比较自然地就可以将整个东西合在一起来考虑，或者甚至将 z i = w i T x + b i z_i = w_i^Tx + b_i zi=wiTx+bi 的定义直接一起带进去然后对 w w w 和 b b b 进行求导来得到 Gradient Descent 的 update rule.
反过来，如果是在设计 Deep Neural Networks 的库，则可能会倾向于将两者分开来看待：因为 Deep Learning 的模型都是一层一层叠起来的结构，一个计算库的主要工作是提供各种各样的 layer，然后让用户可以选择通过不同的方式来对各种 layer 组合得到一个网络层级结构就可以了。比如用户可能最终目的就是得到各个类别的概率似然值，这个时候就只需要一个 Softmax Layer，而不一定要进行 Multinomial Logistic Loss 操作；或者是用户有通过其他什么方式已经得到了某种概率似然值，然后要做最大似然估计，此时则只需要后面的 Multinomial Logistic Loss 而不需要前面的 Softmax 操作。因此提供两个不同的 Layer 结构比只提供一个合在一起的 Softmax-Loss Layer 要灵活许多。从代码的角度来说也显得更加模块化。但是这里自然地就出现了一个问题：numerical stability。
- Softmax-Loss 单层损失层的梯度计算
假设我们直接使用一层 Softmax-Loss 层，计算输入数据
z
k
z_k
zk 属于类别
y
y
y 的概率的极大似然估计。由于 Softmax-Loss 层是最顶层的输出层，则可以直接用最终输出 (loss):
ℓ
~
(
y
,
z
)
\tilde{\ell}(y, z)
ℓ~(y,z) 求对输入
z
k
z_k
zk 的偏导数：
∂
ℓ
~
(
y
,
z
)
∂
z
k
=
∂
∂
z
k
(
log
(
∑
j
=
1
m
e
z
j
)
−
z
y
)
=
exp
(
z
k
)
∑
j
=
1
m
exp
(
z
j
)
−
δ
k
y
=
σ
k
(
z
)
−
δ
k
y
\begin{aligned} \frac{\partial \tilde{\ell}(y, z)}{\partial z_{k}} & = \frac{\partial}{\partial z_{k}} \left(\log \left(\sum_{j=1}^{m} e^{z_{j}}\right)-z_{y} \right)\\ & = \frac{\exp \left(z_{k}\right)}{\sum_{j=1}^{m} \exp \left(z_{j}\right)}-\delta_{k y}=\sigma_{k}(z)-\delta_{k y} \end{aligned}
∂zk∂ℓ~(y,z)=∂zk∂(log(j=1∑mezj)−zy)=∑j=1mexp(zj)exp(zk)−δky=σk(z)−δky
其中
σ
k
(
z
)
\sigma_k(z)
σk(z) 是 Softmax-Loss 的中间步骤 Softmax 在 Forward Pass 的计算结果，而
δ
k
y
=
{
1
k
=
y
0
k
≠
y
\delta_{k y}=\left\{\begin{array}{ll}{1} & {k=y} \\ {0} & {k \neq y}\end{array}\right.
δky={10k=yk̸=y
即 Softmax-Loss 层的梯度为：
∂
ℓ
~
(
y
,
z
)
∂
z
k
=
{
σ
k
(
z
)
−
1
,
k
=
y
σ
k
(
z
)
,
k
≠
y
\begin{aligned} \frac{\partial \tilde{\ell}(y, z)}{\partial z_{k}} = \begin{cases} \sigma_{k}(z) - 1 , & k = y \\ \sigma_{k}(z) , &k \ne y \end{cases} \end{aligned}
∂zk∂ℓ~(y,z)={σk(z)−1,σk(z),k=yk̸=y
- Softmax + Multinomial Logistic Loss 两层分开叠加构造的损失层梯度的计算
接下来看，如果是 Softmax 层和 Multinomial Logistic Loss 层分成两层会是什么样的情况呢？继续回忆刚才的记号：我们把 Softmax 层的输出，也就是 Loss 层的输入记为
o
i
=
σ
i
(
z
)
o_{i}=\sigma_{i}(z)
oi=σi(z)，因此我们首先要计算顶层的 Multinomial Logistic Loss 层输出，对 Softmax 层输入的梯度：
∂
ℓ
(
y
,
o
)
∂
o
i
=
∂
∂
o
i
(
−
log
(
o
y
)
)
=
−
δ
i
y
o
y
\begin{aligned} \frac{\partial \ell(y, o)}{\partial o_{i}} & = \frac{\partial}{\partial o_{i}} \left(-\log \left(o_{y}\right)\right) \\ & = -\frac{\delta_{i y}}{o_{y}} \end{aligned}
∂oi∂ℓ(y,o)=∂oi∂(−log(oy))=−oyδiy
即 Multinomial Logistic Loss 层梯度为：
∂
ℓ
(
y
,
o
)
∂
o
i
{
−
1
o
y
,
i
=
y
0
,
i
≠
y
\begin{aligned} \frac{\partial \ell(y, o)}{\partial o_{i}} \begin{cases} -\frac{1}{o_{y}}, &i = y \\ 0 , &i \ne y \end{cases} \end{aligned}
∂oi∂ℓ(y,o){−oy1,0,i=yi̸=y
然后我们把这个导数向下传递，现在到达 Softmax 层，在 apply chain rule 之前，首先计算层内的导数
∂
o
i
∂
z
k
=
∂
∂
z
k
exp
(
z
i
)
∑
j
=
1
m
exp
(
z
j
)
=
δ
i
k
e
z
i
(
∑
j
=
1
m
e
z
j
)
−
e
z
i
e
z
k
(
∑
j
=
1
m
e
z
j
)
2
=
δ
i
k
o
i
−
o
i
o
k
\begin{aligned} \frac{\partial o_{i}}{\partial z_{k}} & = \frac{\partial}{\partial z_{k}} \frac{\exp \left(z_{i}\right)}{\sum_{j=1}^{m} \exp \left(z_{j}\right)} \\ & =\frac{\delta_{i k} e^{z_{i}}\left(\sum_{j=1}^{m} e^{z_{j}}\right)-e^{z_{i}} e^{z_{k}}}{\left(\sum_{j=1}^{m} e^{z_{j}}\right)^{2}} \\ & =\delta_{i k} o_{i}-o_{i} o_{k} \end{aligned}
∂zk∂oi=∂zk∂∑j=1mexp(zj)exp(zi)=(∑j=1mezj)2δikezi(∑j=1mezj)−eziezk=δikoi−oiok
即 Softmax 层的梯度为：
∂
o
i
∂
z
k
{
o
i
(
1
−
o
k
)
,
k
=
i
−
o
i
o
k
,
k
≠
i
\begin{aligned} \frac{\partial o_{i}}{\partial z_{k}} \begin{cases} o_i(1 - o_k), & k = i \\ -o_i o_k, & k \ne i \end{cases} \end{aligned}
∂zk∂oi{oi(1−ok),−oiok,k=ik̸=i
如果用 Chain Rule 带进去验算一下的话：
∑
i
=
1
m
∂
o
i
∂
z
k
⋅
∂
ℓ
(
y
,
o
)
∂
o
i
=
(
δ
i
k
o
i
−
o
i
o
k
)
⋅
(
−
δ
i
y
o
y
)
⇓
根
据
链
式
法
则
，
当
i
=
y
才
能
往
前
传
递
，
即
把
上
面
的
i
都
用
y
替
换
即
可
=
(
δ
y
k
o
y
−
o
y
o
k
)
⋅
(
−
1
o
y
)
=
o
k
−
δ
y
k
\begin{aligned} \sum_{i=1}^{m} \frac{\partial o_{i}}{\partial z_{k}} \cdot \frac{\partial \ell(y, o)}{\partial o_{i}} & = (\delta_{i k} o_{i}-o_{i} o_{k}) \cdot (-\frac{\delta_{i y}}{o_{y}}) \\ & {\color{Red}{\small{\Downarrow 根据链式法则，当 i = y \ 才能往前传递，即把上面的\ i\ 都用\ y \ 替换即可}}} \\ & = (\delta_{y k} o_{y}-o_{y} o_{k}) \cdot (-\frac{1}{o_{y}}) \\ & = o_{k}-\delta_{y k} \end{aligned}
i=1∑m∂zk∂oi⋅∂oi∂ℓ(y,o)=(δikoi−oiok)⋅(−oyδiy)⇓根据链式法则，当i=y 才能往前传递，即把上面的 i 都用 y 替换即可=(δykoy−oyok)⋅(−oy1)=ok−δyk
和刚才的结果一样的，看来我们求导没有求错。虽然最终结果是一样的，但是我们可以看出，如果分成两层计算的话，要多算好多步骤，除了计算量增大了一点，我们更关心的是数值上的稳定性。由于浮点数是有精度限制的，每多一次运算就会多累积一定的误差，注意到分成两步计算的时候我们需要计算 δ i y / o y \delta_{iy}/o_y δiy/oy 这个量，如果碰巧这次预测非常不准， o y o_y oy 的值，也就是正确的类别所得到的概率非常小（接近零）的话，这里会有 overflow 的危险。下面我们来实际试验一下，首先定义好两种不同的计算函数：
function softmax(z)
#z = z - maximum(z)
o = exp(z)
return o / sum(o)
end
function gradient_together(z, y)
o = softmax(z)
o[y] -= 1.0
return o
end
function gradient_separated(z, y)
o = softmax(z)
∂o_∂z = diagm(o) - o*o'
∂f_∂o = zeros(size(o))
∂f_∂o[y] = -1.0 / o[y]
return ∂o_∂z * ∂f_∂o
end
然后由于 float (Float32) 比 double (Float64) 的精度要小很多，我们就以 double 的计算结果为近似的“正确值”，然后来比较两种情况下通过 float 来计算得到的结果和正确值之差。绘图代码如下：
using DataFrames
using Gadfly
M = 100
y = 1
zy = vec(10f0 .^ (-38:5:38)) # float range ~ [1.2*10^-38, 3.4*10^38]
zy = [-reverse(zy);zy]
srand(12345)
n_rep = 50
discrepancy_together = zeros(length(zy), n_rep)
discrepancy_separated = zeros(length(zy), n_rep)
for i = 1:n_rep
z = rand(Float32, M) # use float instead of double
discrepancy_together[:,i] = [begin
z[y] = x
true_grad = gradient_together(convert(Array{Float64},z), y)
got_grad = gradient_together(z, y)
abs(true_grad[y] - got_grad[y])
end for x in zy]
discrepancy_separated[:,i] = [begin
z[y] = x
true_grad = gradient_together(convert(Array{Float64},z), y)
got_grad = gradient_separated(z, y)
abs(true_grad[y] - got_grad[y])
end for x in zy]
end
df1 = DataFrame(x=zy, y=vec(mean(discrepancy_together,2)),
label="together")
df2 = DataFrame(x=zy, y=vec(mean(discrepancy_separated,2)),
label="separated")
df = vcat(df1, df2)
format_func(x) = @sprintf("%s10<sup>%d</sup>", x<0?"-":"",int(log10(abs(x))))
the_plot = plot(df, x="x", y="y", color="label",
Geom.point, Geom.line, Geom.errorbar,
Guide.xticks(ticks=int(linspace(1, length(zy), 10))),
Scale.x_discrete(labels=format_func),
Guide.xlabel("z[y]"), Guide.ylabel("discrepancy"))
这里我们做的事情是保持
z
z
z 的其他坐标不变，而改变
z
y
z_y
zy 也就是对应于真是 label 的那个坐标的数值大小，我们刚才的推测是当
o
y
o_y
oy 很接近零的时候会有 overflow 的危险，而
o
y
=
σ
y
(
z
)
o_y = \sigma_y(z)
oy=σy(z)，忽略掉 normalization 的话，正比于
e
x
p
(
z
y
)
exp(z_y)
exp(zy)，所以我们需要把
z
y
z_y
zy 那个坐标设成绝对值很大的负数。在得到的图中我们可以看到以整个数值范围内的情况对比。 图中横坐标是
z
y
z_y
zy 的大小，纵坐标是分别用两种方法计算出来的结果和“真实值”之间的差距大小。
首先可以看到的是单层直接计算确实比分成两层算要好一点，不过从纵坐标上也可以看到两者差距其实非常小。往左边看的话，会发现黄色的点没有了，那是因为结果得到了 NaN了，比如
o
y
o_y
oy 由于求一个绝对值非常大的负数的 exponential，导致下溢超出 float 可以表示的小数点精度范围，直接变成 0 了，此时
1
/
o
y
1/o_y
1/oy 就是 Inf，当要乘以
o
y
o_y
oy 进行 cancel 的时候得到
0
×
∞
0 \times \infty
0×∞，对于浮点数这个操作会直接得到 NaN，也就是 Not a Number。反过来看蓝线的话，好像有点奇怪的是越往左边好像反而变得更加精确了，其实是因为我们的“真实值”也 underflow 了，因为 double 虽然比 float 精度高很多，但是也是有限制的。根据 Wikipedia，float 的精度范围大致是
1
0
−
38
∼
1
0
38
10^{-38} \sim 10^{38}
10−38∼1038，而 double 的精度范围大致是
1
0
−
308
∼
1
0
308
10^{-308} \sim 10^{308}
10−308∼10308，大了很多，但是我们不妨来看一下图中的
−
1
0
2
-10^2
−102 这个坐标点，注意到
e
x
=
1
0
x
/
log
10
e^{x}=10^{x / \log 10}
ex=10x/log10
所以
exp
(
−
1
0
2
)
≈
1
0
−
44
\exp \left(-10^{2}\right) \approx 10^{-44}
exp(−102)≈10−44，对于 float 来说已经下溢了，对于 double 来说还是可以表示的范围，但是和 0 的差别也已经如此小，在图上已经看不出区别来了。指数再移一格的话，
exp
(
−
1
0
3
)
≈
1
0
434
\exp \left(-10^{3}\right) \approx 10^{434}
exp(−103)≈10434，会直接导致 double 也 underflow，结果我们的“真实值”也会是零，所以“误差”直接变成零了。
比较有趣的是往右边的正数半轴看，发现到了
1
0
2
10^2
102 之后蓝线和黄线都没有了，说明他们都得到了 NaN，不过这里是另一个问题：对一个比较大的数求 exponential 非常容易发生 overflow。还是用刚才的式子可以看到 ，已经超过了 float 可以表达的最大上限，所以会变成 Inf，然后在 normalize 的一步会出现 Inf/Inf 这样的情况，于是就得到 NaN 了。
这个问题其实也是有解决办法的，我们刚才贴的代码里的 softmax 函数第一行有一行被注释掉的代码，就是 在求 exponential 之前将
z
z
z 的每一个元素减去
z
i
z_i
zi 的最大值。这样求 exponential 的时候会碰到的最大的数就是 0 了，不会发生 overflow 的问题，但是如果其他数原本是正常范围，现在全部被减去了一个非常大的数，于是都变成了绝对值非常大的负数，所以全部都会发生 underflow，但是 underflow 的时候得到的是 0，这其实是非常 meaningful 的近似值，而且后续的计算也不会出现奇怪的 NaN。
证明：将输入 z z z 的每一个元素都减去 z i z_i zi 元素中的最大值，然后求 Softmax 函数结果相等
σ i ( z ) = e z i ∑ j = 1 m e z j = e z i − z m a x ∑ j = 1 m e z j − z m a x = e z i / e z m a x ∑ j = 1 m e z j / e z m a x = e z i ∑ j = 1 m e z j \begin{aligned} \sigma_i(z) & = \frac{e^{z_i}}{\sum_{j = 1}^{m}e^{z_j}} \\ & = \frac{e^{z_i - z_{max}}}{\sum_{j = 1}^{m}e^{z_j - z_{max}}} \\ & = \frac{e^{z_i }/e^{z_{max}}}{\sum_{j = 1}^{m}e^{z_j}/e^{z_{max}}} \\ & = \frac{e^{z_i}}{\sum_{j = 1}^{m}e^{z_j}} \end{aligned} σi(z)=∑j=1mezjezi=∑j=1mezj−zmaxezi−zmax=∑j=1mezj/ezmaxezi/ezmax=∑j=1mezjezi
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
