Softmax Regression

作者：光火

邮箱：victor_b_zhang@163.com

• $SoftmaxRegression$，也称多项逻辑回归，是 $LogisticRegression$ 的一般化方法，可用于处理多分类问题，是机器学习的经典模型之一

Softmax

• 首要，咱们介绍 $Softmax$ 函数

  g(z)i=ezi∑j=1kezj=eziezi+∑j=iej∈(0,1)\Large{g(z)_i=\frac{e^{z_i}}{\sum_{j=1}^{k}e^{z_j}} = \frac{e^{z_i}}{e^{z_i} + \sum_{j\not = i} e^j} \in (0, 1)}g(z)i​=∑j=1k​ezj​ezi​​=ezi​+∑j=i​ejezi​​∈(0,1)

  多分类问题中，模型关于样本类别的打分值往往涣散在整个实数域。对此，咱们能够经过 $\max$ 函数获取得分最高的类别，并以此作为样本类别的猜测值

  $Softmax$ 函数也能够起到这样的效果，并且它还能额外给出样本属于各个类别的概率，并保证概率和为 $1$，因此刻常被用在网络的最终一层

  不过，$Softmax$ 函数也存在必定的缺陷，就比方它是 $over-parameterized$ 的（具有无穷最优解）。以及在编程实现时，由于 exe^xex 增长过快，很容易产生溢出。咱们会在下文结合实例，介绍相应的处理方案

• 机器学习中，咱们一般将 $Softmax$ 函数写成如下方法：

  P(y(i)=j∣x(i),)=exp(jTx(i))∑k=1Kexp(kTx(i))\Large P(y^{(i)}=j|x^{(i)},\theta) = \frac{exp(\theta_j^T x^{(i)})}{\sum_{k=1}^K exp(\theta_k^T x^{(i)})}P(y(i)=j∣x(i),)=∑k=1K​exp(kT​x(i))exp(jT​x(i))​

  该式就代表了将样本 x(i)x^{(i)}x(i) 分为第 $j$ 类的概率，其间 x(i)x^{(i)}x(i) 的上标 $(i)$ 表明的是 $x$ 在数据集抑或是当前 $mini-batch$ 中的位置。从几许视点分析，j\theta_jj​ 就相当于一个超平面，担任给出样本属于第 $j$ 类的分数

  接下来，咱们就结合 $Softmax$ 函数，经过最大似然估量 ($MLE$) ，推导出穿插熵丢失函数（$Cross-entropy$）

Cross-entropy

• 关于给定数据集 $X$，散布 $P(X)$ 和 $Q(X)$ 的穿插熵被界说为：

  H(P,Q)=Ep[−log⁡Q]=−∑x∈XP(x)log⁡Q(x)H(P,Q) = E_p[-\log Q] = – \sum_{x \in X} P(x) \log Q(x)H(P,Q)=Ep​[−logQ]=−x∈X∑​P(x)logQ(x)

  穿插熵的概念来自于信息论，它能够衡量散布 $P$ 与 $Q$ 之间的类似度，两者越附近，穿插熵就越小

• 咱们对上式进行改写，让 $Q(x)$ 除以一个 $P(x)$ 再乘以一个 $P(x)$

  −∑x∈XP(x)log⁡Q(x)=−∑x∈XP(x)log⁡Q(x)P(x)P(x)−∑x∈XP(x)log⁡Q(x)=−∑x∈XP(x)log⁡P(x)−∑x∈XP(x)log⁡Q(x)P(x)- \sum_{x \in X} P(x) \log Q(x) = – \sum_{x \in X}P(x) \log \frac{Q(x) }{P(x)} P(x) \\ – \sum_{x \in X} P(x) \log Q(x) = – \sum_{x \in X} P(x) \log P(x) – \sum_{x \in X} P(x)\log \frac{Q(x) }{P(x)}−x∈X∑​P(x)logQ(x)=−x∈X∑​P(x)logP(x)Q(x)​P(x)−x∈X∑​P(x)logQ(x)=−x∈X∑​P(x)logP(x)−x∈X∑​P(x)logP(x)Q(x)​

  其间，−∑x∈XP(x)log⁡P(x) – \sum_{x \in X} P(x) \log P(x)−∑x∈X​P(x)logP(x) 为概率散布 $P(x)$ 的熵，记作 $H(P)$，用于衡量该散布的混乱程度，并满足 $H(P)\ge 0$

  −∑x∈XP(x)log⁡Q(x)P(x) – \sum_{x \in X} P(x)\log \frac{Q(x) }{P(x)}−∑x∈X​P(x)logP(x)Q(x)​ 则被称为相对熵，抑或是 $KL$ 散度，记作 DKL(P∣∣Q)D_{KL}(P||Q)DKL​(P∣∣Q)

  所以，CE=H(P)+DKL(P∣∣Q)≥0CE = H(P) + D_{KL}(P||Q) \geq 0CE=H(P)+DKL​(P∣∣Q)≥0

• 进一步，倘若 $P$ 是固定的，则最小化 $CE$ 就等价于最小化 DKL(P∣∣Q)D_{KL}(P||Q)DKL​(P∣∣Q)。换句话说，最小化穿插熵便是让 $P$ 与 $Q$ 尽可能地挨近。所以，咱们马上想到，假如让 $P$ 是 $label$ 的真实值，$Q$ 为咱们的猜测值，这不便是一个绝佳的丢失函数吗？

  ---

• 清晰了穿插熵的含义后，咱们经过最大似然估量，建立它与 $Softmax$ 的联系

• 多分类问题中，咱们以为类别的标签 $t$ 应当遵守 $Multinoulli$ 散布。这是一个合理的假设，就像在回归问题中，咱们一般假定输出 $y$ 遵守高斯散布，并运用均方丢失函数 ($MSE$) 一样

• 具体来说，关于一个 $k$ 分类问题，在 $Multinoulli$ 散布下，有如下联系式建立（该式的含义就相当于做了 $K$ 次随机试验）

  P(t∣x)=∏k=1KP(tk=1∣x)tk\Large P(t|x) = \prod_{k=1}^{K} P(t_k = 1|x)^{t_k}P(t∣x)=k=1∏K​P(tk​=1∣x)tk​

  其间，$t$ 代表类别的标签，是一个 $One-Hot$ 向量，即 t=(0,..,1,..,0)Tt = (0,..,1,..,0)^Tt=(0,..,1,..,0)T

  P(tk=1∣x)tkP(t_k = 1 |x)^{t_k}P(tk​=1∣x)tk​ 中呈现了两个 tkt_ktk​，坐落左边的 tkt_ktk​ 表明的是 $One-Hot$ 向量 $t$ 的第 $k$ 个分量，是一个随机变量。坐落指数位置的 tkt_ktk​ 则代表 tkt_ktk​ 的具体取值 （$0$ 或许 $1$）

  举个例子，假设有三个类别，且 $P(t=1)=0.2、P(t=2)=0.5,P(t=3)=0.3$

  则 P(t=(1,0,0)T)=P(t1=1)1P(t2=1)0P(t3=1)0=0.2P(t = (1, 0, 0)^T) =P(t_1 = 1)^1 P(t_2 = 1)^0 P(t_3 = 1)^0 = 0.2P(t=(1,0,0)T)=P(t1​=1)1P(t2​=1)0P(t3​=1)0=0.2

  可见，当且仅当 tk=1t_k = 1tk​=1 时，P(tk=1∣x)tkP(t_k = 1|x)^{t_k}P(tk​=1∣x)tk​ 才是有用概率，不然就会因指数部分为 $0$ 而不起效果

  这儿之所以选用独热编码，是由于咱们不希望人为地引入类别之间的数值与间隔联系（类与类之间应当是互相独立的）

• 关于 P(tk=1∣x)P(t_k = 1|x)P(tk​=1∣x)，咱们指定它为 $Softmax$ 函数（留意：这是人为设置的，假如换用其他的函数，就无法导出穿插熵）这也对应了推导的一般流程，即随机试验应当选用什么散布，以及这个散布的参数应当用什么函数迫临（这儿咱们选用了 $Multinoulli$ 散布和 $Softmax$ 函数）

  读者能够自行证明，假如是二分类问题，根据最大似然估量，使用 $Bernoulli$ 散布和 $Sigmoid$ 函数，就能够导出二分类下的 $Cross-entropy$

  LetP(tk=1∣x)=exp((k)Tx)∑j=1Kexp((j)Tx)=hk(x)∵P(t∣x)=∏k=1KP(tk=1∣x)tk∴关于样本间互相独立的数据集{(x(1),t(1)),..,(x(N),t(N))}有P(t(1),…,t(N)∣X)=∏n=1N∏k=1KP(tk(n)=1∣x(n))tk(n)Let \quad P(t_k = 1 | x) = \frac{exp(\theta^{(k)T}x)}{\sum_{j=1}^K exp(\theta^{(j)T}x)} = h_k(x) \\ \because P(t|x) = \prod_{k=1}^K P(t_k = 1|x)^{t_k} \\ \\ \therefore 关于样本间互相独立的数据集 \quad \{(x^{(1)}, t^{(1)}),..,(x^{(N)},t^{(N)})\} \quad 有 \\ P(t^{(1)},…,t^{(N)}|X) = \prod_{n = 1}^N \prod_{k = 1}^K P(t_k^{(n)}=1|x^{(n)})^{t_k^{(n)}}LetP(tk​=1∣x)=∑j=1K​exp((j)Tx)exp((k)Tx)​=hk​(x)∵P(t∣x)=k=1∏K​P(tk​=1∣x)tk​∴对于样本间彼此独立的数据集{(x(1),t(1)),..,(x(N),t(N))}有P(t(1),...,t(N)∣X)=n=1∏N​k=1∏K​P(tk(n)​=1∣x(n))tk(n)​

  最大化该似然函数，就等价于最小化其负对数似然函数

  E()=−1N∑n=1N∑k=1Ktk(n)ln⁡exp((k)Tx(n))∑j=1Kexp((j)Tx(n))E()=−1N∑n=1N∑k=1Ktk(n)ln⁡hk(n)(∗)E()=−1N∑n=1Nln⁡P(tk(n)=1∣x(n))E(\theta) = -\frac{1}{N} \sum_{n=1}^N \sum_{k=1}^K t_k^{(n)} \ln \frac{exp(\theta^{(k)T}x^{(n)})}{\sum_{j=1}^K exp(\theta^{(j)T}x^{(n)})} \\ E(\theta) = -\frac{1}{N} \sum_{n=1}^N \sum_{k=1}^K t_k^{(n)} \ln h_k^{(n)} \quad (*) \\ E(\theta) = -\frac{1}{N} \sum_{n = 1}^N \ln P(t_k^{(n)} = 1|x^{(n)})E()=−N1​n=1∑N​k=1∑K​tk(n)​ln∑j=1K​exp((j)Tx(n))exp((k)Tx(n))​E()=−N1​n=1∑N​k=1∑K​tk(n)​lnhk(n)​(∗)E()=−N1​n=1∑N​lnP(tk(n)​=1∣x(n))

  第二步中，咱们用 hk(n)h_k^{(n)}hk(n)​ 代表整个分式。第三步中，由于关于一个确定的样本 x(n)x^{(n)}x(n)，t(n)t^{(n)}t(n) 中有且仅有一个 tk(n)t_k^{(n)}tk(n)​ 为 $1$，所以内层的求和号就被消去了

  此刻，咱们回看 $(\ast )$ 式，它便是穿插熵丢失函数的方法，所有的 tk(n)t_k^{(n)}tk(n)​ 就构成了真实标签的 t(n)t^{(n)}t(n) 的 $One-Hot$ 向量，而 ln⁡hk\ln h_klnhk​ 正是咱们的猜测值

• 总结：H(P,Q)=−∑x∈XP(x)log⁡Q(x)H(P,Q) = – \sum_{x \in X} P(x) \log Q(x)H(P,Q)=−∑x∈X​P(x)logQ(x)，关于 $SoftmaxRegression$ 而言，$P(x)$ 即为真实标签（$label$）的 $One-Hot$ 编码，$Q(x)$ 则为经过 $Softmax$ 函数效果后的概率值

  ---

• 最终，咱们对丢失函数求梯度，获取权重 $W$ 的更新方法：

  W=W−[Softmax(WXT)−yT]XX=[x1Tx2T…xmT]Y=[y1Ty2T…ymT]其间,xi、yi均为列向量W = W – \alpha[Softmax(WX^T)-y^T]X \\ X = \begin{bmatrix} x_1^T \\ x_2^T \\ … \\ x_m^T \end{bmatrix} \quad Y = \begin{bmatrix} y_1^T \\ y_2^T \\ … \\ y_m^T \end{bmatrix} \\ 其间,x_i、y_i 均为列向量W=W−[Softmax(WXT)−yT]XX=⎣⎢⎢⎢⎡​x1T​x2T​...xmT​​⎦⎥⎥⎥⎤​Y=⎣⎢⎢⎢⎡​y1T​y2T​...ymT​​⎦⎥⎥⎥⎤​其中,xi​、yi​均为列向量

  至此，咱们已经根本了解 $SoftmaxRegression$，能够尝试用它处理一些简略的问题了

MNIST Classification

• 这儿，咱们结合 MNIST，做一个手写体数字辨认任务（下述代码只展现思路）

MNIST digits dataset is a widely used dataset for image classification in machine learning field. It contains 60,000 training examples and 10,000 testing examples. The digits have been size-normalized and centered in a fixed-size image. Each example is a 784 1 matrix, which is transformed from an original 28 28 gray-scale image. Digits in MNIST range from 0 to 9.

• 略去数据的读入部分

• 超参数设置，实际操作时应根据验证集进行调整

  max_epoch = 10
  batch_size = 100
  learning_reate = 0.01
  # L2 正则项的系数
  factor = 0.01
  

• 初始化权重，并编写训练进程的代码

  # weight initialization
  w = np.random.randn(28 * 28, 10) * 0.001
  loss_list = []
  accuracy_list = []
  display_frequency = 100
  # training
  for epoch in range(0, max_epoch):
      counts = train_size // batch_size
      for count in range(0, counts):
          # 获取下一个 batch 
          x, label = train_set.next_batch(batch_size)
          # 调用 softmax 分类器
          loss, gradient, prediction = softmax_classifier(w, x, label, factor)
          # 核算准确率
          label = np.argmax(label, axis=1)
          accuracy = sum(prediction == label) / float(len(label))
          # 收集数据，方便可视化展现
          loss_list.append(loss)
          accuracy_list.append(accuracy)
          # 更新权重
          w -= learning_rate * gradient
          if count % display_frequency == 0:
              print("Epoch [{}][{}]\t Batch [{}][{}]\t Training Loss {:.4f}\t Accuracy {:.4f}".format(
                  epoch, max_epoch, count, counts, loss, accuracy))
  

• 实现 $Softmax$ 分类器

  def softmax_classifier(w, x, label, factor):
      """
        Softmax Classifier
        Inputs have dimension D, there are C classes, a mini-batch have N examples.
        Inputs:
        - w: A numpy array of shape (D, C) containing weights.
        - x: A numpy array of shape (N, D) containing a mini-batch of data.
        - label: A numpy array of shape (N, C) containing labels, label[i] is a
          one-hot vector, label[i][j]=1 means i-th example belong to j-th class.
        - factor: regularization strength, which is a hyper-parameter.
        Returns:
        - loss: a single float number represents the average loss over the mini-batch.
        - gradient: shape (D, C), represents the gradient with respect to weights W.
        - prediction: shape (N,), prediction[i]=c means i-th example belong to c-th class.
      """
      grade = softmax(np.matmul(x, w))
      prediction = np.argmax(grade, axis=1)
      gradient = np.matmul(x.T, (grade - label)) + 2 * factor * w
      loss = - np.sum(label * np.log(grade + 1e-5)) + factor * np.square(w).sum()
      return loss, gradient, prediction
  

  在书写代码时，咱们应当尽可能地运用矩阵乘法来替代 for 循环，由于前者在效率上要优胜得多

  具体到本例，$x$ 是由 $N$ 个样本组成的矩阵，它的每一行便是 $mini-batch$ 中的一个样本。设样本的维度为 $D$，则 $x$ 的形状为 $ND$

  $w$ 是权重矩阵，它的维度为 $DC$，其间 $C$ 为类别数。由于咱们是对手写数字 $0-9$ 进行分类，所以 $C=10$。$w$ 的第 $i$ 列就对应了第 $i$ 类的各个参数

  因此，咱们经过 np.matmul 将 $x$ 与 $w$ 相乘，得到一个 $NC$ 的矩阵。该矩阵的每一行就对应了该样本在各个类别上的打分值。随后，咱们将 $x\ast w$ 送入 $softmax$ 函数，将得分值转化为概率

  np.argmax(grade, axis=1) 用于取出矩阵 $grade$ 每行的最大值所对应的索引，它便是咱们对样本类别的猜测值

  $gradient$、$loss$ 分别对应在这个 $mini-batch$ 上核算得到的梯度和丢失函数值。咱们引入 $L2$ 正则项，防止 $Softmax$ 的多解性。这儿还有一个技巧，便是 np.log(grade + 1e-5)。为了防止 np.exp(x) 产生溢出，咱们会在必要时，让 $x\ast w$ 全体减去其最大值（见 $Softmax$ 函数代码），这就会导致矩阵中有 $0$ 值呈现，从而无法取对数，所以咱们就引入一个细小的正数

  def softmax(x):
      m = x.max()
      if m > 20:
          x -= x.max()
      return np.exp(x) / np.sum(np.exp(x), axis=1, keepdims=True)
  

  关于 np.sum 的参数：

  axis = 0 代表紧缩行，行将每一列的元素相加，将矩阵紧缩为一行

  axis = 1 代表紧缩列，行将每一行的元素相加，将矩阵紧缩为一列

  keepidims = True 结果保持原来的维数，默以为 False

• 进行测验

  correct = 0
  counts = test_size // batch_size
  # Test process
  for count in range(0, counts):
      data, label = test_set.next_batch(batch_size)
      # We only need prediction results in testing
      _, _, prediction = softmax_classifier(w, data , label, factor)
      label = np.argmax(label, axis=1)
      correct += sum(prediction == label)
  accuracy = correct * 1.0 / test_size
  print('Test Accuracy: ', accuracy)
  

• 绘制曲线

  # training loss curvep
  lt.figure()
  plt.plot(loss_list, 'b--')
  plt.xlabel('iteration')
  plt.ylabel('loss')
  # training accuracy curve
  plt.figure()
  plt.plot(accuracy_list, 'r--')
  plt.xlabel('iteration')
  plt.ylabel('accuracy')