softmax分类(多项逻辑回归)的Pyhon实现及其与SVM的比较

GREAT THANKS TO:http://cs231n.github.io/linear-classify/#softmax

• 
  
• 1..softmax分类
  
  • SVM和softmax是两大常用的分类方法，softmax分类器是二项逻辑回归在多分类问题上的应用.
  • 多项逻辑回归： $P(Y=k|x)=\frac{exp(\omega_kx)}{1+\sum\limits_{k=1}^{K-1}exp(\omega_kx)} \space,\hspace{0.2in}k=1,2,..,K$ $P(Y=K|x)=\frac{1}{1+\sum\limits_{k=1}^{K-1}exp(\omega_kx)}$
  • 在多分类支持向量机中，输出的评分是未经过修剪的，没有一个直观的解释。在softmax分类方法中,其输出是归一化的类的概率，有了直观的意义。
  • 其评分函数为:$\hspace{0.2in} f(\omega)=\omega x+b$
  • softmax分类方法的损失函数：$\hspace{0.2in} L_i = -\log\left(\frac{e^{f_{y_i}}}{ \sum_j e^{f_j} }\right) \hspace{0.2in} \text{也可写成} \hspace{0.2in} L_i = -f_{y_i} + \log\sum_j e^{f_j}$
  • 信息论观点的解释:实际的分布p和估计的分布q定义的交叉熵为:$\hspace{0.2in}H(p,q) = - \sum_x p(x)\log q(x)\hspace{0.2in}$,在softmax分类方法中，实际的分布p就是p=[0,1,...]只有在$y_i$处的概率为1其余为0,估计的分布就是$q = e^{f_{y_i}} / \sum_j e^{f_j}$,所以softmax分类方法就是要最小化上述分布组成的交叉熵。交叉熵可以写成熵和相对熵(KL散度)的和，$H(p,q) = H(p) + D_{KL}(p||q)$。
  • 概率观点的解释：$P(y_i \mid x_i; W) = \frac{e^{f_{y_i}}}{\sum_j e^{f_j} }\hspace{0.2in}$,从此式中可得,其输出是W与相应的X做运算得到的类别$y_i$的概率。从概率的角度来说，softmax损失函数是最小化正确类的负对数似然，也就是说，对正确的类进行最大似然估计，这样也可以把正则化项$R(W)$看成是关于W的先验概率分布，若使用的是$L_2$正则化，则W是服从高斯分布的。这样就是在通过最大化后验概率来分类的。
  • 实际中的问题：数值稳定性:因为这里面用到了指数运算$e^{f_{y_i}}$，有可能会产生很大的值，这会导致计算机在运算的时候不稳定。解决这个问题常用的方法是$\frac{e^{f_{y_i}}}{\sum_j e^{f_j}}= \frac{Ce^{f_{y_i}}}{C\sum_j e^{f_j}}= \frac{e^{f_{y_i} + \log C}}{\sum_j e^{f_j + \log C}}$,在这里，令$\log C = -\max_j f_j$。
  • 可能的困惑:在多分类支持向量机中使用的合页损失函数，也称为最大间隔损失(max-margin loss),而，softmax分类器使用的是交叉熵损失函数，将原始的类别评分压缩到了归一化的和为1的正值。不能称之为softmax损失，因为softmax函数在此只是相当于激励函数。
• 2.SVM和softmax分类的比较
  
  
  • 上图中比较了使用同一个评分函数的，svm和softmax的输出值和损失值。值得注意的是这些损失值不具备可比性，只在当前的分类器中损失值的比较才有意义。
  • softmax函数提供了每个类的“概率”，SVM给出了未经裁剪不好解释的每个类的评分，但是softmax给出了每个类可能的“概率”。例如，SVM可能给出[12.5,0.6,-23.0]作为类“猫”，“狗”和“船”的评分。softmax给出的是这三个类对应的概率[0.9,0.09,0.01],可以将其当成对应的类的置信水平。但是，在这里之所以给概率加上了引号是因为，这里的概率还取决于正则参数$\lambda$的使用。例如，假如说某三个类的对数可能性对应输出为[1,-2,0],则其sotfmax计算的结果为: $[1, -2, 0] \rightarrow [e^1, e^{-2}, e^0] = [2.71, 0.14, 1] \rightarrow [0.7, 0.04, 0.26]$ 现在如果正则化参数$\lambda$变大，则权重参数W将被惩罚的更厉害，这样得到的权值也就更小，假如权重小了一半则f输出为[0.5,-1,0]，则softmax的计算结果为： $[0.5, -1, 0] \rightarrow [e^{0.5}, e^{-1}, e^0] = [1.65, 0.37, 1] \rightarrow [0.55, 0.12, 0.33]$ 可见正则化参数调大时，softmax输出的概率更趋于一致。所以啊，softmax输出的概率也同SVM的分数一样，只是表示一个置信水平，在原理上也没有可解释性。
  • 实际应用中SVM和sotfmax通常拿来比较，SVM和softmax的分类性能区别是很小的。SVM使用的是更局部的目标函数数， 可以把其当成缺点也可以当成特点。SVM不考虑单个类别的细节，像[10,-100,-100]和[10,9,9]则对于$\Delta = 1$的支持向量机是相同的，因为最大间隔已经满足了，其损失值就为0。但是，对softmax分类器来说这两种情况就不同了,很明显[10,9,9]的损失值更大。softmax分类器永远不会对其输出的分数满意，它会不停的要求提高正确类的概率降低错误类的概率，减小损失值。但是，SVM却不是这样，只要误分类的评分比正确分类的评分小最大间隔那么大就OK了，这也可以被当作一个特殊的功能，譬如为了训练一个货车和汽车的分类器，那么这个分类器应该致力于将汽车和货车分开而不受到已经评分很低的青蛙的影响(这很可能已经聚集为另一个不同的类了)。
• 3.程序实现：
  

数据获取地址同 上一篇多分类支持向量机

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Sun Jul 29 17:15:25 2018
@author: rd
"""
from __future__ import division
import numpy as np
"""
This dataset is part of MNIST dataset,but there is only 3 classes,
classes = {0:'0',1:'1',2:'2'},and images are compressed to 14*14 
pixels and stored in a matrix with the corresponding label, at the 
end the shape of the data matrix is 
num_of_images x 14*14(pixels)+1(lable)
"""
def load_data(split_ratio):
    tmp=np.load("data216x197.npy")
    data=tmp[:,:-1]
    label=tmp[:,-1]
    mean_data=np.mean(data,axis=0)
    train_data=data[int(split_ratio*data.shape[0]):]-mean_data
    train_label=label[int(split_ratio*data.shape[0]):]
    test_data=data[:int(split_ratio*data.shape[0])]-mean_data
    test_label=label[:int(split_ratio*data.shape[0])]
    return train_data,train_label,test_data,test_label

"""compute the cross entropy loss without using vector operation,
While dealing with a huge dataset,this will have low efficiency
X's shape [n,14*14+1],Y's shape [n,],W's shape [num_class,14*14+1]"""
def lossAndGradNaive(X,Y,W,reg):
    dW=np.zeros(W.shape)
    loss = 0.0
    num_class=W.shape[0]
    num_X=X.shape[0]
    for i in range(num_X):
        scores=np.dot(W,X[i])
        cur_scores=scores[int(Y[i])]
        loss+=-cur_scores+np.log(np.sum(np.exp(scores)))
        for j in range(num_class):
            if j==int(Y[i]):
                dW[int(Y[i]),:] += -X[i]+np.exp(cur_scores)/np.sum(np.exp(scores))*X[i]
            else:
                dW[j,:]+=np.exp(scores[j])/np.sum(np.exp(scores))*X[i]
    loss/=num_X
    dW/=num_X
    loss+=reg*np.sum(W*W)
    dW+=2*reg*W
    return loss,dW

def predict(X,W):
    X=np.hstack([X, np.ones((X.shape[0], 1))])
    Y_=np.dot(X,W.T)
    Y_pre=np.argmax(Y_,axis=1)
    return Y_pre

def accuracy(X,Y,W):
    Y_pre=predict(X,W)
    acc=(Y_pre==Y).mean()
    return acc

def model(X,Y,alpha,steps,reg):
    X=np.hstack([X, np.ones((X.shape[0], 1))])
    W = np.random.randn(3,X.shape[1]) * 0.0001
    for step in range(steps):
        loss,grad=lossAndGradNaive(X,Y,W,reg)
        W-=alpha*grad
        print"The {} step, loss={}, accuracy={}".format(step,loss,accuracy(X[:,:-1],Y,W))
    return W

train_data,train_label,test_data,test_label=load_data(0.2)
W=model(train_data,train_label,0.0001,25,0.5)
print"Test accuracy of the model is {}".format(accuracy(test_data,test_label,W))

---

refer
[1] https://zhuanlan.zhihu.com/p/31562236