(三)Softmax的原理及代码实现

Softmax

本系列重点在浅显易懂，快速上手。不进行过多的理论讲解：也就是不去深究what，而是关注how。全文围绕以下三个问题展开：

1）长什么样？

2）解决什么问题？

3）怎么实现？

​ 3.1）从数学讲，原理

​ 3.2）从代码上讲，如何掉包实现

长什么样

上节讲了逻辑回归实现2分类，那么多分类的问题，就要通过softmax来实现，在此之前回忆“线性回归”和“逻辑回归”

1）线性回归建模，用于预测y是连续值的问题，值域的范围是(-∞,+∞)最终的目的就是为了找到一组$${\theta }_1,{\theta }_2,...,{\theta }_{13},{\theta }_0$$（除了系数，还包含一个常数项$${\theta }_0$$,使得下式成立。

$${\theta }_1\ast x_1+{\theta }_2\ast x_2+......+{\theta }_k\ast x_k+{\theta }_0=y$$

2)逻辑回归建模，用于解决二分类问题，y值代表的是取1的概率。值域的范围是[0,1]最终的目的就是为了找到一组$${\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_0$$（除了4个系数，还包含一个常数项$${\theta }_0$$）,使得下式成立。
$$\frac{1}{1+e^{-({\theta }_1\ast x_1+{\theta }_2\ast x_2+{\theta }_3\ast x_3+{\theta }_k\ast x_k+{\theta }_0)}}=\frac{1}{1+e^{-{\theta }^{T}x}}=y$$

3)softmax模型，用于解决多分类问题，假如要解决3分类问题。模型的输出应该是什么形式的? 假如还是按逻辑回归的方式，输出一个[0,1]之间的值，那这个值应该代表取到类别1，类别2，类别的概率？所以，应该同时输出3个值，$$y_1,y_2,y_3$$，再配个分母，即：
$$\begin{gathered} \left[ \\ \frac{y_1}{y_1+y_2+y_3}, \\ \frac{y_2}{y_1+y_2+y_3}, \\ \frac{y_3}{y_1+y_2+y_3} \\ \right] \end{gathered}$$

分别代表是类别1，类别2，类别3的概率，加起来为1。具体的函数解析式，应该是如下形式：

$$\begin{gathered} {\theta }_{1\_1}\ast x_1+{\theta }_{1\_2}\ast x_2+......+{\theta }_{1\_n}\ast x_n+{\theta }_{1\_0}=y_1 \\ {\theta }_{2\_1}\ast x_1+{\theta }_{2\_2}\ast x_2+......+{\theta }_{2\_n}\ast x_n+{\theta }_{3\_0}=y_2 \\ {\theta }_{3\_1}\ast x_1+{\theta }_{3\_2}\ast x_2+......+{\theta }_{3\_n}\ast x_n+{\theta }_{3\_0}=y_3 \end{gathered}$$

相应的参数也是原来的3倍，即：

Expected node of symbol group type, but got node of type cr

但还有一个问题，$$y_1,y_2,y_3$$，的值有可能是负值，这是不合理的，所以取e指数，取了e指数后各个值之间的差异也会被放大。即：

$$\begin{gathered} \left[ \\ \frac{e^{y_1}}{e^{y_1+y_2+y_3}}, \\ \frac{e^{y_2}}{e^{y_1+y_2+y_3}}, \\ \frac{e^{y_3}}{e^{y_1+y_2+y_3}} \\ \right] \end{gathered}$$

写成一般形式如下:

$$p(y=k\mid x;\theta )=\frac{e^{{\theta }_k^{T}x}}{{\sum }_{l=1}^K{e}^{{\theta }_l^{T}x}},k=1,2\cdots ,K$$

即：

$$h_{\theta }(x)=\left[\begin{array}{c}p(y^{(i)}=1\mid x^{(i)};\theta )\\ p(y^{(i)}=2\mid x^{(i)};\theta )\\ ...\\ p(y^{(i)}=k\mid x^{(i)};\theta )\end{array}\right]=\frac{1}{{\sum }_{j=1}^k{e}^{{\theta }_j^T{x}^{(i)}}}\left[\begin{array}{c}{e}^{{\theta }_1^{T}x}\\ e^{{\theta }_2^{T}x}\\ ...\\ e^{{\theta }_k^{T}x}\end{array}\right]$$

$$⟹\theta =\left[\begin{array}{cccc}{\theta }_{11}& {\theta }_{12}& ...& {\theta }_{1n}\\ {\theta }_{21}& {\theta }_{22}& ...& {\theta }_{2n}\\ ...& ...& ...& ...\\ {\theta }_{k1}& {\theta }_{k2}& ...& {\theta }_{kn}\end{array}\right]$$

解决什么问题

先看一组数据（一组sklearn中公开的数据集，鸢尾花的种类预测，共有4个特征，最后一列y是花的种类，共有3种类型，这里只展示了6条数据，实际数据有150条）,读到这里，有读者会问，y不是应该是3个值吗，这里为什么是1个，这其实是一种简略写法，以第一个y为例子，0表示这个花是类别0，写成程序能理解的方式应该是[1,0,0]。同理，第二个y为1表示，这个花是类别是1，写成程序能理解的方式应该是[0,1,0]。以此类推。

x1	x2	x3	x4	y
5.1	3.7	1.5	0.4	0
5.8	2.7	3.9	1.2	1
5.1	3.8	1.9	0.4	0
6.7	3.3	5.7	2.5	2
5.5	2.4	3.7	1	1
5.9	3	5.1	1.8	2

采用softmax建模，最终的目的就是为了找到一组

$$\theta =\left[\begin{array}{ccccc}{\theta }_{11}& {\theta }_{12}& {\theta }_{13}& {\theta }_{14}& {\theta }_{10}\\ {\theta }_{21}& {\theta }_{22}& {\theta }_{23}& {\theta }_{24}& {\theta }_{20}\\ ...& ...& ...& ...& \\ {\theta }_{k1}& {\theta }_{k2}& {\theta }_{k3}& {\theta }_{k4}& {\theta }_{k0}\end{array}\right]$$

（除了4个系数，还包含一个常数项${\theta }_0$),使得下式成立。

$$\begin{gathered} {\theta }_{11}\ast x_1+{\theta }_{12}\ast x_2+{\theta }_{13}\ast x_3+{\theta }_{14}\ast x_4+{\theta }_{10}=y_1 \\ {\theta }_{21}\ast x_1+{\theta }_{22}\ast x_2+{\theta }_{23}\ast x_3+{\theta }_{24}\ast x_4+{\theta }_{20}=y_2 \\ {\theta }_{31}\ast x_1+{\theta }_{32}\ast x_2+{\theta }_{33}\ast x_3+{\theta }_{34}\ast x_4+{\theta }_{30}=y_3 \\ p(y=k\mid x;\theta )=\frac{e^{{\theta }_k^{T}x}}{{\sum }_{l=1}^K{e}^{{\theta }_l^{T}x}},k=1,2,3 \end{gathered}$$

怎么实现

数学原理

和上节讲到的逻辑回归一样的思路，通过最大似然估计法，得到损失函数。实际上分类问题的损失函数是一样的，称为交叉熵损失函数（也叫对数损失函数），这个不必再推到一边，直接可以写出来：

还是看上文的第一条样本，真实值是[1,0,0]，假设预测值是[0.7,0.2,0.1]。则，对应的交叉熵损失函数是：

$$J(\theta )=1\ast log(0.7)+0\ast log(0.2)+0\ast log(0.1)$$

是不是很简单，写成一般式，即:

$$J(\theta )=-\frac{1}{m}\sum _{i=1}^m\sum _{j=1}^{k}I(y^{(i)}=j)\ln \left(\frac{e^{{\theta }_j^T{x}^{(i)}}}{{\sum }_{i=1}^k{e}^{{\theta }_j^T{x}^{(i)}}}\right)I(y^{(i)}=j)=\{\begin{array}{ll}1,& y^{(i)}=j\\ 0,& y^{(i)}\ne j\end{array}$$

k是类别的个数，在上文中是3，m表示样本的个数，上文中是150。
同样的，有了损失函数(目标函数)，后求损失函数的最小值，先求导：

$$\begin{array}{r}\frac{\partial }{\partial {\theta }_j}J(\theta )=\frac{\partial }{\partial {\theta }_j}-I\left(y^{(i)}=j\right)\ln \left(\frac{e^{{\theta }_j^T{x}^{(i)}}}{{\sum }_{l=1}^k{e}^{{\theta }_l^T{x}^{(i)}}}\right)\\ =\frac{\partial }{\partial {\theta }_j}-I(y^{(i)}=j)({\theta }_j^T{x}^{(i)}-\ln (\sum _{l=1}^k{e}^{{\theta }_l^T{x}^{(i)}}\left)\right)\\ =-I(y^{(i)}=j)\left(1-\frac{e^{{\theta }_j^T{x}^{(i)}}}{{\sum }_{l=1}^k{e}^{{\theta }_l^T{x}^{(l)}}}\right)x^{(i)}\end{array}$$

写成矩阵的形式：
$$\frac{\partial }{\partial {\theta }_j}J(\theta )=-1\ast X^T\ast [(1-h_{\theta }(x))\cdot Y]$$

$\ast$表示矩阵的乘法，$\cdot$ 表示对应元素相乘。

求出倒数后，就利用梯度下降法求解$\theta$，不断的迭代，直到损失函数的值不再降低。

$$\theta =\theta -\alpha \cdot \frac{\partial }{\partial {\theta }_j}J(\theta )$$

手动代码实现

#_*_coding=utf-8_*_
import numpy as np
import random

class softmax():
    
    def __init__(self,X,Y,theta=None,epoch=10000,batch=100,eta=0.001,is_standard=True):
        '''
        self.X = m*n  m个样本，n个特征，
        self.Y = m*k  k个分类
        self.theta = n*k
        '''
        self.X = X
        self.Y = Y
        self.theta = theta
        self.epoch = epoch
        self.eta = eta
    
    def pre(self):
        predict = np.exp(np.dot(self.X, self.theta))
        predict = predict/np.sum(predict,axis=1).reshape(-1,1)
        return predict
    
    def gra(self):
        predict = self.pre()
        #predict = np.dot(self.X, self.theta)
        mark = np.zeros([self.X.shape[0],len(np.unique(self.Y))])
        mark[np.arange(self.X.shape[0]),self.Y.reshape(-1)] = 1
        result = -1*np.dot(self.X.T, np.multiply(1 - predict,mark))/self.X.shape[0]
        return result
   
    def cal_grad(self):
        '''
        梯度更新可以用矩阵相乘计算
        #grad = np.zeros(self.theta.shape)   #grad n*1      data  10*n   data的转置 n*10   predict 10*1
        #for i in range(len(grad)):
            #grad[i] = np.mean((predict - y)*self.data[:,i])
        #return grad'''
        predict = np.dot(self.X,self.theta)   
        y_hat = self.softmax(predict)
        
        y_onehot = np.array([[1 if int(m[0])==i else 0 for i in range(3)] for m in self.Y]) 
        return -np.dot(self.X.T,((1-y_hat)*y_onehot))/len(self.X)    
    
    
    
    
    def los(self):
        predict = self.pre()
        loss = -1/self.X.shape[0]*np.sum(np.log(predict)[np.arange(self.X.shape[0]),self.Y.reshape(-1)])
        return loss
    
    def  standard_scaler(self):
        mean = np.mean(self.X,axis=0)
        std = np.std(self.X,axis=0)
        self.X = (self.X-mean)/std
    
    @staticmethod    
    def cal(y,predict):
        result = []
        for i in range(3):
            TP = sum([1 for index,lable in enumerate(y) if lable==i and predict[index]==i])
            FN = sum([1 for index,lable in enumerate(y) if lable==i and predict[index]!=i])
            FP = sum([1 for index,lable in enumerate(y) if lable!=i and predict[index]==i])
            TN = sum([1 for index,lable in enumerate(y) if lable!=i and predict[index]!=i])
            a = round((TP+TN)/(TP+FN+FP+TN),4)
            p = round(TP/(TP+FP),4)
            r = round(TP/(TP+FN),4)
            f1 = round(2*p*r/(p+r),4)
            result.append([a,p,r,f1])
        return np.array(result)     
    
    
    def model(self):
        '''
        data m*n+1 包含标签
        theta  n*k
        return m*k
        '''
        predict = np.dot(self.X,self.theta)
        return self.softmax(predict)
    
    def softmax(self,values):
        '''
        values m*k  
        [[1,2,3,4],
         [2,3,4,5],
         [1,2,3,5]]
        
        retrun  y_hat m*k
        '''
        values = np.exp(values)
        sum_h = np.sum(values,axis=1)     # 按行求和 sum_h m*1
        sum_h = sum_h.reshape((-1,1))
        y_hat = values/sum_h
        return y_hat
    
    def cross(self):
        '''
        y_hat,predict 训练数据样本数*类别数  m*k
        theta   特征数*类别           n*k
        x 训练数据  训练样本数量*特征数 m*n
        self.data[:,-1:] m*1
        y_onehot m*k
        return 一个值
        '''
        predict = np.dot(self.X,self.theta)   
        y_hat = self.softmax(predict)
        # 转换成one_hot编码
        #y_onehot = np.array([[1 if int(m[0])==i else 0 for i in range(len(self.k))] for m in self.train_data[:,-1:]]) 
        #return -np.sum(y_onehot*np.log(y_hat))/len(self.train_data) 
        h = [m[int(self.Y[index])] for index,m in enumerate(y_hat)] 
        return -np.sum(np.log(h))/len(self.X)
        
    


    
if __name__ == "__main__":
    np.random.seed(2023)
    import sklearn.datasets as datasets
    data = datasets.load_iris()
    X = data.data
    
    Y = data.target.reshape(-1,1)
    #X = np.array(random.sample(X.tolist(),X.shape[0]))
    
    lr = 0.01
    sf = softmax(X, Y)
    sf.standard_scaler()
    ones = np.ones([sf.X.shape[0],1])
    sf.X = np.hstack([ones,sf.X]) 
    sf.theta = np.random.randn(sf.X.shape[1],len(np.unique(sf.Y)))  
    n = 1
    while n<100000:
        loss = sf.los()

        
        
        print(n,loss)

        gradient = sf.gra()
        #print(gradient)
        #gradient = sf.cal_grad()
        #print(gradient)        

        sf.theta = sf.theta - lr*gradient
        n+=1
        #sf.X = np.array(random.sample(sf.X.tolist(),len(sf.X)))
        #if n%10000==0:
            #predict = np.argmax(sf.pre(),axis=1)
            #print(sf.cal(sf.Y.reshape(-1), predict))
    #print(np.argmax(sf.pre(),axis=1))
    print(sf.theta)
    print(np.sum(np.argmax(sf.pre(),axis=1).reshape(-1,1)==sf.Y)/sf.X.shape[0])
    print('Done')

sklearn掉包实现

这里直接讲解，如何调用sklearn实现。实际上和二分类的代码完全一致。

from sklearn.linear_model import LogisticRegression #sklearn中，线性回归模型在linear_model模块中
# 调取sklearn中自带的数据集
from sklearn.datasets import load_iris #调用上文一开始提到大波士顿房价数据集

X, y = load_boston(return_X_y=True) #获取X，y数据
LR = LogisticRegression() #初始化一个l回归模型
LR.fit(X,y) #fit函数用于训练

LR.predict(X) #predict函数用于预测，直接返回类别值
LR.predict_proba(X) #返回每个类别的概率，
LR.coef_ #coef_用于返回θ的值
LR.intercept_ #intercept_用于返回常数项的值(截距)
LR.score(X, y) #给出预测的准确率。

更多内容可以参考sklearn的官方文档

https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression