【机器学习与算法】python手写算法：softmax回归

算法原理

softmax回归用于解决多分类问题。它的基本思想是计算样本属于每一个类别的概率，属于哪个类别的概率最大，则预测输出为哪一类。softmax计算概率的方式为：
$$P(y^{(i)}=j|x^i;\theta )=\frac{e^{{\theta }_j^T\cdot x^{(i)}}}{{\sum }_{l=1}^K{e}^{{\theta }_l^T\cdot x^{(i)}}}$$
表示样本$x^{(i)}$属于第j个类别的概率。
softmax回归的损失函数如下：
$$J(\theta )=-\frac{1}{m}\sum _{i=1}^m\sum _{j=1}^{K}l(y^{(i)}=j)\ast ln\frac{e^{{\theta }_j^T\cdot x^{(i)}}}{{\sum }_{l=1}^K{e}^{{\theta }_l^T\cdot x^{(i)}}}$$
其中$l(\cdot )$是一个示性函数，括号里面的逻辑判断为真，则函数返回1，否则返回0。直观上理解就是让每个样本预测正确类别的联合概率最大化。这个损失函数叫交叉熵损失，为什么呢，我们看下交叉熵的公式：
$$H(p,q)=-\sum _{i=1}^{n}p(x_i)\ast ln(q(x_i))$$
其中p表示真实概率，q表示预测概率，和softmax损失函数是一致的。
当然，逻辑回归的损失函数$J(\theta )=-{\sum }_{i=i}^m(y\ast ln(\hat{y})+(1-y)\ast ln(1-\hat{y}))$也叫交叉熵损失，其实就是类别数K=2时候的一种交叉熵的特殊形式而已。
softmax损失函数的梯度如下：
$$\frac{\partial }{\partial {\theta }_j}J(\theta )=-\frac{1}{m}\sum _{i=1}^m[x^{(i)}(l(y^{(i)}=j)-P(y^{(i)}=j|x^i;\theta ))]+\lambda {\theta }_j$$
具体的推导不再赘述，如果对softmax原理的更详细解释和推导有兴趣，包括softmax损失函数梯度的求导，参数冗余的特性和softmax回归如何在K=2时推演为逻辑回归，请参阅这篇博客：softmax回归(Softmax Regression)

python实现算法

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib.patches as mpatches
from sklearn.preprocessing import PolynomialFeatures
from sklearn import datasets
import copy


class softmax:
    
    def __init__(self, K, alpha=0.0001, lamda=0.0001, tol = 1e-4, max_iter=1000):
        #类别数
        self.K = K
        #学习率
        self.alpha = alpha
        #正则项参数
        self.lamda = lamda
        #停止迭代的标准
        self.tol = tol
        #最大迭代轮数
        self.max_iter = max_iter

        
    def cal_p(self,W,X,k):
    	'''
    	计算所有样本在每个类别上的概率
    	'''

    	#分子：样本在该类别上线性方程的指数转换
        numerator = np.exp(np.dot(W[k],X.T))
        #分母：样本在所有类别上线性方程指数转换的加和
        denominator = np.sum(np.exp(np.dot(W,X.T)),axis=0)
        #返回概率
        return numerator/denominator
        
        
    def fit(self, X, Y):
    	'''
		训练模型
		'''

        n = len(X[0])
        #初始化参数W，是一个类别数 X 特征数的矩阵
        W = np.zeros((self.K,n))
        #初始化一个w_com用于和W进行比较，在差异小于tol的时候停止迭代
        w_com = copy.deepcopy(W)

        for times in range(self.max_iter):
            for k in range(self.K):
            	#计算示性函数项
                l = [1 if item==k else 0 for item in Y]
                #计算所有样本在该类别上的概率
                p = self.cal_p(W,X,k)
                #目标函数的梯度项
                g = np.sum((l-p)*X.T,axis=1)
                #沿着梯度进行下降
                W[k] = W[k] - (-self.alpha*g + self.lamda*W[k])
                
            #判断是否满足停止迭代的条件，满足即跳出循环
            if abs(np.sum(w_com)-np.sum(W))<self.tol:
                break
            else:
                w_com = copy.deepcopy(W)
        
        #把拟合好的参数W存到softmax类的属性中    
        self.__setattr__('W',W)
    
    
    def predict(self,X):
    	'''
    	预测
    	'''
        
        n = len(X)
        #计算样本属于每个类别k的概率
        prob = np.zeros((self.K, n))
        for k in range(self.K):
            prob[k] = self.cal_p(self.W, X, k)
        
        #返回类别中最大概率的下标，作为预测输出    
        res = np.argmax(prob,axis=0)
        
        return res

结果展示

为了方便用二维图片展示结果，只取‘花瓣长度’和’花瓣宽度’两个特征来训练。但是只选两个特征并不代表我们只有两个特征可用，我们可以使用sklearn里面的PolynomialFeatures来给这两个特征做升维。2个特征x1，x2升维后变6个：1，$x_1$，$x_2$，$x_1^2$，$x_2^2$，$x1\ast x2$。

if __name__ == '__main__':

	#用鸢尾花数据试验
    iris = datasets.load_iris()
    #为了方便用二维图片展示结果，只取‘花瓣长度'和'花瓣宽度'两个特征来训练
    X = iris.data[:,2:]
    Y = iris.target

    #用polynomial函数对特征进行二阶升维，升维后2个特征变6个：x1,x2升维到1，x1，x2，x1**2，x2**2, x1*x2
    poly = PolynomialFeatures(2)
    x_p = poly.fit_transform(X)

    #初始化softmax类并训练
    sf = softmax(K=3)
    sf.fit(x_p,Y)

    #这个是为了给最后的图片加个边
    def extend(a, b):
        return 1.05*a-0.05*b, 1.05*b-0.05*a

    N, M = 200, 200     # 横纵各采样多少个值
    x1_min, x1_max = extend(X[:, 0].min(), X[:, 0].max())   # 第0列的范围
    x2_min, x2_max = extend(X[:, 1].min(), X[:, 1].max())   # 第1列的范围
    t1 = np.linspace(x1_min, x1_max, N)
    t2 = np.linspace(x2_min, x2_max, M)
    x1, x2 = np.meshgrid(t1, t2)

    x_show = np.stack((x1.flat, x2.flat), axis=1)   # 测试点
    x_show_p = poly.fit_transform(x_show)

    y_hat = sf.predict(x_show_p)
    y_hat = y_hat.reshape(x1.shape)  # 使之与输入的形状相同


    cm_light = mpl.colors.ListedColormap(['#77E0A0', '#FF8080', '#A0A0FF'])
    cm_dark = mpl.colors.ListedColormap(['g', 'r', 'b'])
    mpl.rcParams['font.sans-serif'] = u'SimHei'
    mpl.rcParams['axes.unicode_minus'] = False
    plt.figure(facecolor='w')
    plt.pcolormesh(x1, x2, y_hat, cmap=cm_light)  # 预测值的显示
    plt.scatter(X[:, 0], X[:, 1], s=30, c=Y, edgecolors='k', cmap=cm_dark)  # 样本的显示
    x1_label, x2_label = '花瓣长度', '花瓣宽度'
    plt.xlabel(x1_label, fontsize=14)
    plt.ylabel(x2_label, fontsize=14)
    plt.xlim(x1_min, x1_max)
    plt.ylim(x2_min, x2_max)
    plt.grid(b=True, ls=':')

    patchs = [mpatches.Patch(color='#77E0A0', label='Iris-setosa'),
              mpatches.Patch(color='#FF8080', label='Iris-versicolor'),
              mpatches.Patch(color='#A0A0FF', label='Iris-virginica')]
    plt.legend(handles=patchs, fancybox=True, framealpha=0.8, loc='lower right')
    plt.title(u'鸢尾花softmax回归分类', fontsize=17)
    plt.show()

sklearn实现softmax回归

sklearn上没有softmax这样一个接口，那怎么使用softmax算法实现多分类呢？其实sklearn把softmax算法统一封装到LogisticRegression这个接口里了，我们只需要调用这个接口，然后改下multi_class参数为‘multinomial’就可以了，记得solver参数不要选用‘liblinear’：

from sklearn.linear_model import LogisticRegression

sf_sklearn = LogisticRegression(multi_class='multinomial', solver='saga',max_iter=1000)
sf_sklearn.fit(x_p,Y)
sf_sklearn.predict(x_p)

我们会发现multi_class这个参数不光‘multinomial’一个选择啊，还有一个‘ovr’，同样可以解决多分类问题。ovr的意思就是one versus all，就是预测每一类概率的时候，把这一类样本给个标签1，剩下的类别样本给个标签0，然后用二分类去预测，最后同样把概率值最高的类别作为输出。那这两种方法分别适用于什么情况呢，经验性来讲，softmax适用于各个类别相互独立的情况，而ovr适用于各个类别互有交集的情况。

---

最后，欢迎阅读其它算法的python实现：
【机器学习与算法】python手写算法：Cart树
【机器学习与算法】python手写算法：带正则化的逻辑回归
【机器学习与算法】python手写算法：xgboost算法
【机器学习与算法】python手写算法：Kmeans和Kmeans++算法