1. 结合Logistic Regression 分析 Softmax
训练集由 \(m\) 个已标记的样本构成：\(\{(x^{(1)}, y^{(1)}), ... , (x^{(m)}, y^{(m)})\}\)，其中输入的第 \(i\) 个样例 \(x^{(i)} \in \Re^{n+1}\)，即特征向量 \(x\) 的维度为 \(n + 1\) ，其中 \(x_0 = 1\) 对应截距项。由于 logistic 回归是针对二分类问题的，因此类标记 \(y^{(i)} \in \{0,1\}\)。假设函数(hypothesis function) 如下：
\[
h_\theta(x) = \frac{1}{1+e^{(-\theta^Tx)}}
\]
我们将训练模型参数 \(\theta\)，使其能够最小化代价函数 ：
\[
\begin{align}
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]
\end{align}
\]
而在 Softmax回归中，我们解决的是多分类问题（相对于 logistic 回归解决的二分类问题），类标 \(\textstyle y\) 可以取 \(\textstyle k\) 个不同的值（而不是 2 个）。因此，对于训练集 \(\{ (x^{(1)}, y^{(1)}), \ldots, (x^{(m)}, y^{(m)}) \}\)，我们有 \(y^{(i)} \in \{1, 2, \ldots, k\}\)。（注意此处的类别下标从 1 开始，而不是 0）。例如，在 MNIST 数字识别任务中，我们有 \(\textstyle k=10\) 个不同的类别。
对于输入的每一个样例 \(x^{(i)} \in \Re^{n+1}\) ，输出其属于每一类别 \(j\) 的概率值 \(\textstyle p(y=j | x)\)，即输出一个 \(k\) 维向量（向量元素的和为1）来表示这 \(k\) 个估计的概率值。即 Softmax 的假设函数为 \(h_{\theta}(x^{(i)})\)
\[
\begin{align}
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}
\end{align}
\]
其中 \(\theta_1, \theta_2, \ldots, \theta_k \in \Re^{n+1}\) 是模型的参数。请注意 \(\frac{1}{ \sum_{j=1}^{k}{e^{ \theta_j^T x^{(i)} }}}\) 这一项对概率分布进行归一化，使得所有概率之和为 1
其代价函数为：
\[
\begin{align}
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]
\end{align}
\]
其中 \(1\{\cdot\}\) 为指示函数：
\[
1\{\cdot\} =
\begin{cases}
1\{表达式值为真\} = 1 \\
1\{表达式值为假\} = 0 \\
\end{cases}
\]
梯度下降公式：
\[
\begin{align}
\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] }
\end{align}
\]
2. Softmax
结合上面理解此时的 Softmax函数：
\[ h_{\theta}(x) = \sigma(\theta^T x) = \sigma(z) = \left({\color{Red}{\sigma_{1}(z)}}, \ldots, \sigma_{m}(z)\right) \]
此时的 \(x\) 为一个样例输入 \(\in \Re^{n+1}\)，等价与上面的一个 \(x^{(i)}\)，其中共有 \(m\) 个类别，\(z_i\) 表示输入样例 \(x\) 是第 \(i\) 个类别的线性预测结果，等价于 \(z_i = \theta^T_i x\)
Softmax 函数 \(\sigma(z)=\left({\color{Red}{\sigma_{1}(z)}}, \ldots, \sigma_{m}(z)\right)\) 定义如下：
\[
{\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}{\small{【观察到的数据 \ x \ (或z)\ 属于类别\ i\ 的概率，或者称作似然 (Likelihood)】}}}
\]
它在 Logistic Regression 里起到的作用是讲线性预测值转化为类别概率：\(m\) 代表类别数，假设 \(z_i = w_i^Tx + b_i\) 是第 \(i\) 个类别的线性预测结果，带入 \(Softmax\) 的结果其实就是先对每一个 \(z_i\) 取 exponential 变成非负，然后除以所有项之和进行归一化，现在每个 \(o_{i}=\sigma_{i}(z)\) 就可以解释成：观察到的数据 \(x\) 属于类别 \(i\) 的概率，或者称作似然 (Likelihood)。
然后 Logistic Regression 的目标函数是根据最大似然原则来建立的，假设数据 \(x\) 所对应的类别为 \(y\)，则根据我们刚才的计算最大似然就是要最大化 \(o_y\) 的值 (通常是使用 negative log-likelihood 而不是 likelihood，也就是说最小化 \(-log(o_y)\) 的值，这两者结果在数学上是等价的)。后面这个操作就是 caffe 文档里说的 Multinomial Logistic Loss，具体写出来是这个样子：
\[
\ell(y, o)=-\log \left(o_{y}\right)
\]
3. Softmax-Loss = Softmax + Multinomial Logistic Loss
而 Softmax-Loss 其实就是把两者结合到一起，只要把 \(o_y\) 的定义展开即可：
\[
\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\ 似然】}}}
\]
比如如果我们要写一个 Logistic Regression 的 solver，那么因为要处理的就是这个东西，比较自然地就可以将整个东西合在一起来考虑，或者甚至将 \(z_i = w_i^Tx + b_i\) 的定义直接一起带进去然后对 \(w\) 和 \(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\) 属于类别 \(y\) 的概率的极大似然估计。由于 Softmax-Loss 层是最顶层的输出层，则可以直接用最终输出 (loss): \(\tilde{\ell}(y, z)\) 求对输入 \(z_k\) 的偏导数：
\[
\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}
\]
其中 \(\sigma_k(z)\) 是 Softmax-Loss 的中间步骤 Softmax 在 Forward Pass 的计算结果，而
\[
\delta_{k y}=\left\{\begin{array}{ll}{1} & {k=y} \\ {0} & {k \neq y}\end{array}\right.
\]
即 Softmax-Loss 层的梯度为：
\[
\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}
\]
- Softmax + Multinomial Logistic Loss 两层分开叠加构造的损失层梯度的计算
接下来看，如果是 Softmax 层和 Multinomial Logistic Loss 层分成两层会是什么样的情况呢？继续回忆刚才的记号：我们把 Softmax 层的输出，也就是 Loss 层的输入记为 \(o_{i}=\sigma_{i}(z)\)，因此我们首先要计算顶层的 Multinomial Logistic Loss 层输出，对 Softmax 层输入的梯度：
\[
\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}
\]
即 Multinomial Logistic Loss 层梯度为：
\[
\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}
\]
然后我们把这个导数向下传递，现在到达 Softmax 层，在 apply chain rule 之前，首先计算层内的导数
\[
\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}
\]
即 Softmax 层的梯度为：
\[
\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}
\]
如果用 Chain Rule 带进去验算一下的话：
\[
\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}
\]
和刚才的结果一样的，看来我们求导没有求错。虽然最终结果是一样的，但是我们可以看出，如果分成两层计算的话，要多算好多步骤，除了计算量增大了一点，我们更关心的是数值上的稳定性。由于浮点数是有精度限制的，每多一次运算就会多累积一定的误差，注意到分成两步计算的时候我们需要计算 \(\delta_{iy}/o_y\) 这个量，如果碰巧这次预测非常不准，\(o_y\) 的值，也就是正确的类别所得到的概率非常小（接近零）的话，这里会有 overflow 的危险。下面我们来实际试验一下，首先定义好两种不同的计算函数：
1 | function softmax(z) |
然后由于 float (Float32) 比 double (Float64) 的精度要小很多，我们就以 double 的计算结果为近似的“正确值”，然后来比较两种情况下通过 float 来计算得到的结果和正确值之差。绘图代码如下：
1 | using DataFrames |
这里我们做的事情是保持 \(z\) 的其他坐标不变，而改变 \(z_y\) 也就是对应于真是 label 的那个坐标的数值大小，我们刚才的推测是当 \(o_y\) 很接近零的时候会有 overflow 的危险，而 \(o_y = \sigma_y(z)\)，忽略掉 normalization 的话，正比于 \(exp(z_y)\)，所以我们需要把 \(z_y\) 那个坐标设成绝对值很大的负数。在得到的图中我们可以看到以整个数值范围内的情况对比。 图中横坐标是 \(z_y\) 的大小，纵坐标是分别用两种方法计算出来的结果和“真实值”之间的差距大小。
首先可以看到的是单层直接计算确实比分成两层算要好一点，不过从纵坐标上也可以看到两者差距其实非常小。往左边看的话，会发现黄色的点没有了，那是因为结果得到了 NaN了，比如 \(o_y\) 由于求一个绝对值非常大的负数的 exponential，导致下溢超出 float 可以表示的小数点精度范围，直接变成 0 了，此时 \(1/o_y\) 就是 Inf，当要乘以 \(o_y\) 进行 cancel 的时候得到 \(0 \times \infty\)，对于浮点数这个操作会直接得到 NaN，也就是 Not a Number。反过来看蓝线的话，好像有点奇怪的是越往左边好像反而变得更加精确了，其实是因为我们的“真实值”也 underflow 了，因为 double 虽然比 float 精度高很多，但是也是有限制的。根据 Wikipedia，float 的精度范围大致是 \(10^{-38} \sim 10^{38}\)，而 double 的精度范围大致是 \(10^{-308} \sim 10^{308}\)，大了很多，但是我们不妨来看一下图中的 \(-10^2\) 这个坐标点，注意到
\[
e^{x}=10^{x / \log 10}
\]
所以 \(\exp \left(-10^{2}\right) \approx 10^{-44}\)，对于 float 来说已经下溢了，对于 double 来说还是可以表示的范围，但是和 0 的差别也已经如此小，在图上已经看不出区别来了。指数再移一格的话，\(\exp \left(-10^{3}\right) \approx 10^{434}\)，会直接导致 double 也 underflow，结果我们的“真实值”也会是零，所以“误差”直接变成零了。
比较有趣的是往右边的正数半轴看，发现到了 \(10^2\) 之后蓝线和黄线都没有了，说明他们都得到了 NaN，不过这里是另一个问题：对一个比较大的数求 exponential 非常容易发生 overflow。还是用刚才的式子可以看到 ，已经超过了 float 可以表达的最大上限，所以会变成 Inf，然后在 normalize 的一步会出现 Inf/Inf 这样的情况，于是就得到 NaN 了。
这个问题其实也是有解决办法的，我们刚才贴的代码里的 softmax 函数第一行有一行被注释掉的代码，就是 在求 exponential 之前将 \(z\) 的每一个元素减去 \(z_i\) 的最大值。这样求 exponential 的时候会碰到的最大的数就是 0 了，不会发生 overflow 的问题，但是如果其他数原本是正常范围，现在全部被减去了一个非常大的数，于是都变成了绝对值非常大的负数，所以全部都会发生 underflow，但是 underflow 的时候得到的是 0，这其实是非常 meaningful 的近似值，而且后续的计算也不会出现奇怪的 NaN。
证明：将输入 \(z\) 的每一个元素都减去 \(z_i\) 元素中的最大值，然后求 Softmax 函数结果相等
\[ \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} \]