基于Softmax多分类+全批量梯度下降，实现三分类

本文章是笔者在大三时的《模式识别》课程作业，如今将此作业上传供网友学习。

1 $Softmax$多分类的基本原理

1.1 $Softmax$模型

 $softmax$多分类器，一种基于$softmax$函数的分类器，它可以预测一个样本属于每个样本的概率。$softmax$一般用于神经网络的输出层，叫做$softmax$层。
 $softmax$函数适用于处理多分类问题，应用广泛的逻辑函数就是$softmax$函数在二分类情形下的特例。$softmax$函数将一个$n$维的输入向量映射为n维的向量，使得输出向量的各元素取值在0到1之间，且所有元素之和为1，即所得到的向量可以作为事件发生的概率。
记函数的输入向量为:$Z=(z_1,z_2,\cdots ,z_n{)}^{\top }$，则函数值为：
$$softmax(X)=(\frac{e^{x_1}}{{\sum }_{i=1}^n{e}^{x_i}},\frac{e^{x_2}}{{\sum }_{i=1}^n{e}^{x_i}},\cdots ,\frac{e^{x_n}}{{\sum }_{i=1}^n{e}^{x_i}}{)}^{\top }$$

1.2 多分类模型

 对一个激活函数为$softmax$的单个神经元，记输入是一个有m个样本的数据，且每一个数据有n维：
$$X=\left[\begin{array}{ccccc}{x}_{10}& x_{11}& x_{12}& \cdots & x_{1n}\\ x_{20}& x_{21}& x_{22}& \cdots & x_{2n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_{m0}& x_{m1}& x_{m2}& \cdots & x_{mn}\end{array}\right]$$
 其中，$\mathbf{X}$第一列数据全部为1，是人为加入的bias项，表示模型带有偏置（即考虑截距项）。$X_{ij}(1\le i\le m,1\le j\le n)$表示第i个样本的第j个特征的值。
该样本对应的真实类别为：
$$Y=\left[\begin{array}{cccc}{y}_{11}& y_{12}& \cdots & y_{1k}\\ y_{21}& y_{22}& \cdots & y_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ y_{m1}& y_{m2}& \cdots & y_{mk}\end{array}\right]$$
 Y的每一行只有一个值为1，其它全为0。$y_{ij}=1(1\le i\le m,1\le j\le k)$表示第i个样本对应的类别为第j类。
记待估计的参数为
$$\Omega =\left[\begin{array}{cccc}{\omega }_{01}& {\omega }_{02}& \cdots & {\omega }_{0k}\\ {\omega }_{11}& {\omega }_{12}& \cdots & {\omega }_{1k}\\ {\omega }_{21}& {\omega }_{22}& \cdots & {\omega }_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ {\omega }_{n1}& {\omega }_{n2}& \cdots & {\omega }_{nk}\end{array}\right]$$
 记softmax函数的自变量为：
$$Z=\left[\begin{array}{cccc}{z}_{11}& z_{12}& \cdots & z_{1k}\\ z_{21}& z_{22}& \cdots & z_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ z_{m1}& z_{m2}& \cdots & z_{mk}\end{array}\right]$$
 Z是模型输入数据的加权求和，即$z_{ij}={\sum }_{k=0}^n{x}_{ik}{\omega }_{kj}$,$Z=X\Omega$。
记:
$$\hat{Y}=softmax(Z)=\left[\begin{array}{cccc}{\hat{y}}_{11}& {\hat{y}}_{12}& \cdots & {\hat{y}}_{1k}\\ {\hat{y}}_{21}& {\hat{y}}_{22}& \cdots & {\hat{y}}_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{y}}_{m1}& {\hat{y}}_{m2}& \cdots & {\hat{y}}_{mk}\end{array}\right]$$
 其中${\hat{y}}_{ij}=\frac{e^{z_{ij}}}{{\sum }_{p=1}^k{e}^{z_{ip}}}$，表示模型眼中第$i$个样本属于第$j$类的概率。

1.3 代价函数

 定义模型的总代价函数为
$$COST(X)=-\frac{1}{m}\sum _{i=0}^m(\sum _{j=1}^k{y}_{ij}\ln {\hat{y}}_{ij})$$
 将代价函数视为参数$\Omega$的函数$J_{\Omega }$,这就是我们要优化的目标。

 使用梯度下降求解目标函数极小值：
 由链式法则可知，
$$\frac{\partial J_{\Omega }}{\partial {\omega }_{pj}}=\sum _{i=1}^m\frac{\partial J_{\Omega }}{\partial z_{ij}}\frac{\partial z_{ij}}{\partial {\omega }_{pj}}$$
 而
$$\frac{\partial J_{\Omega }}{\partial z_{ip}}=\frac{\partial {\sum }_{q=0}^m(-{\sum }_{j=1}^k\ln {\hat{y}}_{qj})}{\partial z_{ip}}=-\sum _{j=1}^k\frac{y_{ij}}{{\hat{y}}_{ij}}\frac{\partial {\hat{y}}_{ij}}{\partial z_{ip}}$$
 由于
$$\begin{array}{c}\frac{\partial {\hat{y}}_{ij}}{\partial z_{ip}}=\{\begin{array}{rlr}& =\frac{e^{z_{ij}}}{{\sum }_{q=1}^k{e}^{z_{iq}}}-\frac{e^{z_{ij}}{e}^{z_{ip}}}{({\sum }_{q=1}^k{e}^{z_{iq}}{)}^2}={\hat{y}}_{ip}-{\hat{y}}_{ip}^2& j=p\\ & =-\frac{e^{z_{ij}}{e}^{z_{ip}}}{({\sum }_{q=1}^k{e}^{z_{iq}}{)}^2}=-{\hat{y}}_{ij}\widehat{y_{i}p}& j\ne p\end{array}\end{array}$$
将(3)式带入(2)式，由于${\sum }_{j=1}^k{y}_{ij}=1$,有
$$\frac{\partial J_{\Omega }}{\partial z_{ip}}={\hat{y}}_{ip}-y_{ip}$$
考虑$z_{ip}={\sum }_{q=0}^n{x}_{iq}{\omega }_{qp}$，并将(4)式带入(1)式，得到：
$$\frac{\partial J_{\Omega }}{\partial {\omega }_{pj}}=\sum _{i=1}^m\frac{\partial J_{\Omega }}{\partial z_{ij}}{x}_{ip}=\sum _{i=1}^m({\hat{y}}_{ij}-y_{ij})x_{ip}$$

写回矩阵的形式，有：
$$\nabla J_{\Omega }=\left[\begin{array}{cccc}\frac{J_{\Omega }}{{\omega }_{01}}& \frac{J_{\Omega }}{{\omega }_{02}}& \cdots & \frac{J_{\Omega }}{{\omega }_{0k}}\\ \frac{J_{\Omega }}{{\omega }_{11}}& \frac{J_{\Omega }}{{\omega }_{12}}& \cdots & \frac{J_{\Omega }}{{\omega }_{1k}}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{J_{\Omega }}{{\omega }_{n1}}& \frac{J_{\Omega }}{{\omega }_{n2}}& \cdots & \frac{J_{\Omega }}{{\omega }_{nk}}\end{array}\right]=\left[\begin{array}{cccc}{\sum }_{i=1}^m{x}_{i0}({\hat{y}}_{i1}-y_{i1})& {\sum }_{i=1}^m{x}_{i0}({\hat{y}}_{i2}-y_{i2})& \cdots & {\sum }_{i=1}^m{x}_{i0}({\hat{y}}_{ik}-y_{ik})\\ {\sum }_{i=1}^m{x}_{i1}({\hat{y}}_{i1}-y_{i1})& {\sum }_{i=1}^m{x}_{i1}({\hat{y}}_{i2}-y_{i2})& \cdots & {\sum }_{i=1}^m{x}_{i1}({\hat{y}}_{ik}-y_{ik})\\ ⋮& ⋮& \ddots & \cdots \\ {\sum }_{i=1}^m{x}_{in}({\hat{y}}_{i1}-y_{i1})& {\sum }_{i=1}^m{x}_{in}({\hat{y}}_{i2}-y_{i2})& \cdots & {\sum }_{i=1}^m{x}_{in}({\hat{y}}_{ik}-y_{ik})\end{array}\right]$$
梯度已经求出来了，指定步长使用梯度下降方法求解即可。

1.4算法流程

• 输入：训练数据集为 D = { ( x i , y i ) } i = 1 m D=\left\{(x_i,y_i) \right\}_{i=1}^{m} D={(xi​,yi​)}i=1m​,其中,$x_i\in \mathbf{X}\subseteq {\mathbf{R}}^n,y_i\in \mathbf{Y}$；学习率$\eta (0<\eta ⩽1)$
• 过程：
  梯度下降算法迭代更新模型参数直至收敛，每一轮具体流程如下：
  $(1)$选取初值$w$;
  $(2)$更新$w$: 首先通过正向传播求得本轮训练样本预测值，然后反向传播更新模型参数：
  $w_{pj}\leftarrow w_{pj}-\eta \frac{1}{m}\frac{\partial J_{\Omega }}{\partial {\omega }_{pj}}=w_{pj}-eta\frac{1}{m}{\sum }_{i=1}^m({\hat{y}}_{ij}-y_{ij})x_{ip}$
  $(3)$一直重复步骤2，直至模型参数收敛
• 输出：$w$

2 三分类流程

# 导入所需库
import csv
import numpy as np
from matplotlib import pyplot as plt
import pylab

• 读取Iris数据集，并查看部分数据
• 抽取样本：取前两类样本（共150条），将数据集的4个属性作为自变量X。将数据集的3个类别映射为{0,1,2}，作为因变量Y。

# 样本数据的抽取
with open('iris.data') as csv_file:
    data = list(csv.reader(csv_file, delimiter=','))

label_map = {
    'Iris-setosa': 0,
    'Iris-versicolor': 1,
    'Iris-virginica':2
}
# 使用Softmax解决多分类问题

# 抽取样本
X = np.array([[float(x) for x in s[:-1]] for s in data[:150]], np.float32) # X是一个四维数据
Y = np.array([[label_map[s[-1]]] for s in data[:150]], np.float32) 

使用One-Hot独热编码，并且抽取训练集和测试集

# 将 Y 转换为可以用0、1、2表示的向量
tmp = np.zeros((Y.shape[0],3)) # 寄存器
for i in range(Y.shape[0]):
    if Y[i] == 0:
        tmp[i] = [1,0,0]
    if Y[i] == 1:
        tmp[i] = [0,1,0]
    if Y[i] == 2:
        tmp[i] = [0,0,1]
Y = tmp

# 分割数据集

# 将数据集按照8：2划分为训练集和测试集
train_idx = np.random.choice(150, 120, replace=False)

test_idx = np.array(list(set(range(150)) - set(train_idx)))

# train-训练集 test-测试集
X_train, Y_train = X[train_idx], Y[train_idx]
b = np.ones((X_train.shape[0],1)) # 添加常数项
X_train = np.hstack((X_train, b))
X_test, Y_test = X[test_idx], Y[test_idx]
b = np.ones((X_test.shape[0],1)) # 添加常数项
X_test = np.hstack((X_test, b))

# print(X_train)
# print(Y_train)

定义$Softmax$多分类class

class Softmax(object):  # 定义Softmax多分类class
    def __init__(self, X, Y,TX,TY,epoch,eta):  # X为训练集样本,Y是训练集结果， TX为测试集样本,TY是测试集结果,epoch是训练次数,eta是学习率
        if X.shape[0] != Y.shape[0]:  # 要求X,Y中的数目一样
            raise ValueError('Error,X and Y must be same when axis=0 ')
        else:  # 在类中储存参数
            self.X = X
            self.Y = Y
            self.TX = TX
            self.TY = TY
            self.epoch = epoch
            self.eta = eta
            self.m=self.X.shape[0]#m个数据样本
            self.n=self.X.shape[1]#n个特征
            self.k=self.Y.shape[1]#k个类别
    def softmax(self,x,w):
        a=x@w #(m*n)*(n*k),得到的a是一个(m*k)维矩阵
        if x.ndim==1:
            sum=np.sum(np.exp(a))#如果输入的是向量则直接求和
        else:
            sum=np.sum(np.exp(a),axis=1).reshape(x.shape[0],1)
            #如果输入的是矩阵，sum对行求和后输出向量，reshape将向量转化为矩阵，然后除法需要用到矩阵的广播运算
        return np.exp(a)/(sum)#,按行求和，其实得到的a就是y的预测值
    

    def loss(self,w):#损失函数在迭代过程中根本就用不到
        sum=0
        for i in range(self.m):
            y_hat=Softmax.softmax(self,self.X[i],w)
            sum += np.dot(self.Y[i],np.log(y_hat))
        return -sum/self.m    
    
    def grad(self,w):#计算出梯度矩阵
        Y_hat=Softmax.softmax(self,self.X,w)#所有样本的预测值，y_hat是m*k大小的矩阵
        error=Y_hat-self.Y
        # print(self.X.shape)
        # print(error.shape)
        gradient=self.X.T@error#定义梯度矩阵
        # print(gradient)
        return gradient

    
    def test(self,x,y,w):
        num=0
        for i in range(x.shape[0]):
            max_hat_index=np.argmax(Softmax.softmax(self,x[i],w))
            max_y_index=np.argmax(y[i])
            if max_hat_index==max_y_index:
                num += 1
        return 1-num/x.shape[0]


    def train(self):
        train_err=[]
        test_err=[]
        w = np.ones((self.n,self.k))
        for i in range(self.epoch):
            grad_reg=Softmax.grad(self,w)/self.m
            w=w-self.eta*grad_reg# 损失值大于0，计算梯度，更新权值
            # print(w)
            train_err.append(Softmax.test(self,self.X,self.Y,w))
            test_err.append(Softmax.test(self,self.TX,self.TY,w))
        return w,train_err,test_err

训练函数

epoch=500#迭代次数
eta=0.05#学习率
SOFT=Softmax(X_train, Y_train,X_test,Y_test,epoch, eta)
w,train_err,test_err = SOFT.train() 

结果输出和可视化

print("Softmax多分类模型参数:")
print(w)
print('最终训练误差和测试误差:')
print('train error:',train_err[-1])
print('test error:',test_err[-1])
plt.figure(figsize=(12, 6))#设置图片大小
plt.style.use('seaborn-darkgrid') # 设置画图的风格 
# 绘制训练误差 train error
plt.subplot(1, 2, 1)
plt.plot(train_err)
plt.title('Softmax')
plt.xlabel('epoch')
plt.ylabel('train error')
# 绘制测试误差test error
plt.subplot(1, 2, 2)
plt.plot(test_err)
plt.title('Softmax')
plt.xlabel('epoch')
plt.ylabel('test error')

输出结果

Softmax多分类模型参数:
[[ 1.46855411  1.40304271  0.12840318]
 [ 2.25369077  0.69785336  0.04845587]
 [-0.71491953  1.07995137  2.63496816]
 [ 0.21334229  0.49745047  2.28920724]
 [ 1.23399573  1.22614849  0.53985578]]
最终训练误差和测试误差:
train error: 0.025000000000000022
test error: 0.0
Text(0, 0.5, 'test error')

4 实验总结

4.1实验中遇到的状况

 在实验中，借鉴他人的代码时，因为大多数的代码都是将用到的部分全都写为了函数，我主要花费的精力是理解算法过程并将这些函数整合为一个类。

• 而且遇到的代码主要是进行向量运算，我起初以为python和matlab一样，将向量也看作了矩阵，但是代码一直报错，最后通过输出样例，发现了python中向量和矩阵是不同的概念，基于此我将所有的函数都改为了既可以进行矩阵运算也可以进行向量运算，其实就是加了几句判断语句，但是发现这个问题耗费了大量的时间。
• 另外因为测试函数的return 语句的缩进不对，导致测试误差和训练误差的数据一直是凌乱的数据，最后想要逐个输出误差时，发现了该问题。

4.2 实验感悟

• 对$Softmax$完成多分类的问题有了一个初步的认识，对其算法编程过程有了系统性的理解。
• 了解道路三种梯度下降算法，明晰了三种算法的具体区别与优缺点(篇幅限制，未写在作业中)，并使用了全批量梯度下降算法实现了三分类问题。
• 对进行类的封装和增加函数的通用性方面的编程能力有了提升

5 借鉴博客

[1] 桃源梦 \quad 梯度下降处理softmax函数多分类问题 \quad https://www.cnblogs.com/Reader-Yu/p/10665239.html
[2] 明·煜 \quad 作业二·Softmax实现多分类 \quad https://blog.csdn.net/weixin_53195427/article/details/130358620?spm=1001.2014.3001.5502
[3] Twilight Sparkle. \quad 【机器学习笔记13】softmax多分类模型【上篇】完整流程与详细公式推导 \quad https://blog.csdn.net/qq_52466006/article/details/127525149
[4] Twilight Sparkle. \quad 机器学习笔记14】softmax多分类模型【下篇】从零开始自己实现softmax多分类器（含具体代码与示例数据集） \quad https://blog.csdn.net/qq_52466006/article/details/127555367