Python实现softmax回归
1.Softmax回归概念
Softmax回归可以用于多类分类问题，Softmax代价函数与logistic 代价函数在形式上非常类似，只是在Softmax损失函数中对类标记的 个可能值进行了累加。注意在Softmax回归中将 分类为类别 的概率为：
以下公式中， 是示性函数, 。举例来说，表达式 的值为1 ，的值为 0。我们的代价函数为：
对于 的最小化问题，目前还没有闭式解法。因此，我们使用迭代的优化算法（例如梯度下降法，或 L-BFGS）。经过求导，我们得到随机梯度下降公式如下(参考资料6有详解)：
以上求导公式有一个问题，公式中采用了示性函数，经我研究比较适合用在随机梯度下降中，并不适合用在批量梯度下降及mini-batch梯度下降中，即使实现了此示性函数也需要进行y从1到k类的循环判断，比较费时费力，在阅读了大量其他博主的文章后，我发现了另外一种解决方式：
先将数据集转化为逻辑分类可以处理的数据结构。即为对象添加值为1的属性x0（截距项b），将输出分类y转换为one-hot编码。分类1表示为[1,0,0]，分类2表示为[0,1,0]，分类3表示为[0,0,1]，将处理后的y带入，整理后的梯度下降公式如下（公式推导见参考资料9）：
数据处理方法，将分类y转换为one-hot编码形式：
# 随机产生多分类数据，n_samples=样本数，n_features=x数据维度，centers=y分类数
x,y = datasets.make_blobs(n_samples=100, n_features=2, centers=4, cluster_std=1.0, center_box=(-10.0, 10.0), shuffle=True, random_state=None)
plt.scatter(x[:,0],x[:,1],c=y,s=8) # 画图查看数据图像
# 转换数据为matrix类型
x=np.mat(x)
y=np.mat(y).T
# 数据处理，x增加默认值为1的b偏量，y处理为onehot编码类型
def data_convert(x,y):
b=np.ones(y.shape) # 添加全1列向量代表b偏量
x_b=np.column_stack((b,x)) # b与x矩阵拼接
K=len(np.unique(y.tolist())) # 判断y中有几个分类
eyes_mat=np.eye(K) # 按分类数生成对角线为1的单位阵
y_onehot=np.zeros((y.shape[0],K)) # 初始化y的onehot编码矩阵
for i in range(0,y.shape[0]):
y_onehot[i]=eyes_mat[y[i]] # 根据每行y值，更新onehot编码矩阵
return x_b,y,y_onehot
2.批量梯度下降算法
这里直接贴代码了，批量梯度下降算法使用所有数据进行梯度下降迭代。
def SoftmaxGD(x,y,alpha=0.05,max_loop=500):
# 梯度上升算法
# alpha=0.05 # 步长
# max_loop=500 # 循环次数
#x = StandardScaler().fit_transform(x) # 数据进行标准化
m=np.shape(x)[1] # x的特征数
n=np.shape(y)[1] # y的分类数
weights=np.ones((m,n)) # 权重矩阵
for k in range(max_loop):
# k=2
h=softmax(x*weights)
error=y-h
weights=weights+alpha*x.transpose()*error # 梯度下降算法公式
#print('k:',k,'weights:',weights.T)
return weights.getA()
3.随机梯度下降算法
这里直接贴代码了，随机梯度下降算法使用一条数据进行梯度下降迭代。
def SoftmaxSGD(x,y,alpha=0.05,max_loop=50):
# 随机梯度上升算法
# alpha=0.05
# max_loop=500
#x = StandardScaler().fit_transform(x) # 数据进行标准化
m=np.shape(x)[1]
n=np.shape(y)[1]
weights=np.ones((m,n))
for k in range(max_loop):
for i in range(0,len(x)):
# k=0;i=0
h=softmax(x[i]*weights)
error=y[i]-h[0]
weights=weights+alpha*x[i].T*error[0] # 随机梯度下降算法公式
#print('k:',k,'i:',i,'weights:',weights.T)
return weights.getA()
如果你看了，我前几篇线性回归、逻辑回归的代码，你会发现softmax回归只是在逻辑回归的基础上改成了softmax函数，除了需要将y进行改造外，其他代码基本都不用变，是不是很简单呢？
Softmax回归完整代码
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 1 15:46:06 2018
@author: admin
"""
import numpy as np
import pandas as pd
from sklearn.datasets import load_wine
import matplotlib.pyplot as plt
from sklearn import datasets
# 加载数据集，最后一列最为类别标签，前面的为特征属性的值
def loadDataSet(x,y):
# 生成X和y矩阵
#dataMat = np.mat(x)
y = np.mat(y)
b = np.ones(y.shape) # 添加全1列向量代表b偏量
X = np.column_stack((b, x)) # 特征属性集和b偏量组成x
X = np.mat(X)
labeltype = np.unique(y.tolist()) # 获取分类数目
eyes = np.eye(len(labeltype)) # 每一类用单位矩阵中对应的行代替，表示目标概率。如分类0的概率[1,0,0]，分类1的概率[0,1,0]，分类2的概率[0,0,1]
Y=np.zeros((X.shape[0],len(labeltype)))
for i in range(X.shape[0]):
Y[i,:] = eyes[int(y[i,0])] # 读取分类，替换成概率向量。这就要求分类为0,1,2,3,4,5这样的整数
return X, y,Y #X为特征数据集，y为分类数据集，Y为概率集
# 数据处理，x增加默认值为1的b偏量，y处理为onehot编码类型
def data_convert(x,y):
b=np.ones(y.shape) # 添加全1列向量代表b偏量
x_b=np.column_stack((b,x)) # b与x矩阵拼接
K=len(np.unique(y.tolist())) # 判断y中有几个分类
eyes_mat=np.eye(K) # 按分类数生成对角线为1的单位阵
y_onehot=np.zeros((y.shape[0],K)) # 初始化y的onehot编码矩阵
for i in range(0,y.shape[0]):
y_onehot[i]=eyes_mat[y[i]] # 根据每行y值，更新onehot编码矩阵
return x_b,y,y_onehot
# softmax函数，将线性回归值转化为概率的激活函数。输入s要是行向量
def softmax(s):
return np.exp(s) / np.sum(np.exp(s), axis=1)
# 逻辑回归中使用梯度下降法求回归系数。逻辑回归和线性回归中原理相同，只不过逻辑回归使用sigmoid作为迭代进化函数。
def gradAscent(x, y,alpha=0.05,max_loop=500):
# 梯度上升算法
#alpha=0.05 # 步长
#max_loop=500 # 循环次数
weights = np.ones((x.shape[1],y.shape[1])) #初始化权回归系数矩阵 系数矩阵的行数为特征矩阵的列数，系数矩阵的列数为分类数目
print('weights初始化值：',weights)
for k in range(max_loop):
# k=0
h = softmax(x*weights) #梯度上升矢量化公式，计算预测值（行向量）。每一个样本产生一个概率行向量
error = h-y #计算每一个样本预测值误差
weights = weights - alpha * x.T * error # 根据所有的样本产生的误差调整回归系数
#print('k:',k,'weights:',weights)
return weights # 将矩阵转换为数组，返回回归系数数组
def SoftmaxGD(x,y,alpha=0.05,max_loop=500):
# 梯度上升算法
# alpha=0.05 # 步长
# max_loop=500 # 循环次数
#x = StandardScaler().fit_transform(x) # 数据进行标准化
m=np.shape(x)[1] # x的特征数
n=np.shape(y)[1] # y的分类数
weights=np.ones((m,n)) # 权重矩阵
for k in range(max_loop):
# k=2
h=softmax(x*weights)
error=y-h
weights=weights+alpha*x.transpose()*error # 梯度下降算法公式
#print('k:',k,'weights:',weights.T)
return weights.getA()
def SoftmaxSGD(x,y,alpha=0.05,max_loop=50):
# 随机梯度上升算法
# alpha=0.05
# max_loop=500
#x = StandardScaler().fit_transform(x) # 数据进行标准化
m=np.shape(x)[1]
n=np.shape(y)[1]
weights=np.ones((m,n))
for k in range(max_loop):
for i in range(0,len(x)):
# k=0;i=0
h=softmax(x[i]*weights)
error=y[i]-h[0]
weights=weights+alpha*x[i].T*error[0] # 随机梯度下降算法公式
#print('k:',k,'i:',i,'weights:',weights.T)
return weights.getA()
# 多分类只能绘制分界区域。而不是通过分割线来可视化
def plotBestFit(dataMat,labelMat,weights):
# 获取数据边界值，也就属性的取值范围。
x1_min, x1_max = dataMat[:, 1].min() - .5, dataMat[:, 1].max() + .5
x2_min, x2_max = dataMat[:, 2].min() - .5, dataMat[:, 2].max() + .5
# 产生x1和x2取值范围上的网格点，并预测每个网格点上的值。
step = 0.02
xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, step), np.arange(x2_min, x2_max, step))
testMat = np.c_[xx1.ravel(), xx2.ravel()] #形成测试特征数据集
testMat = np.column_stack((np.ones(((testMat.shape[0]),1)),testMat)) #添加第一列为全1代表b偏量
testMat = np.mat(testMat)
# 预测网格点上的值
y = softmax(testMat*weights) #输出每个样本属于每个分类的概率
# 判断所属的分类
predicted = y.argmax(axis=1) #获取每行最大值的位置，位置索引就是分类
predicted = predicted.reshape(xx1.shape).getA()
# 绘制区域网格图
plt.pcolormesh(xx1, xx2, predicted, cmap=plt.cm.Paired)
# 再绘制一遍样本点，方便对比查看
plt.scatter(dataMat[:, 1].flatten().A[0], dataMat[:, 2].flatten().A[0],
c=labelMat.flatten().A[0],alpha=.5) # 第一个偏量为b，第2个偏量x1，第3个偏量x2
plt.show()
# 对新对象进行预测
def predict(weights,testdata):
y_hat=softmax(testdata*weights)
predicted=y_hat.argmax(axis=1).getA()
return predicted
if __name__ == "__main__":
# 随机产生多分类数据，n_samples=样本数，n_features=x数据维度，centers=y分类数
x,y = datasets.make_blobs(n_samples=100, n_features=2, centers=4, cluster_std=1.0, center_box=(-10.0, 10.0), shuffle=True, random_state=None)
plt.scatter(x[:,0],x[:,1],c=y,s=8) # 画图查看数据图像
# 转换数据为matrix类型
x=np.mat(x)
y=np.mat(y).T
# 调用数据预处理函数
X,Y,Y_onehot=data_convert(x,y)
# 批量梯度下降算法
weights1=SoftmaxGD(X, Y_onehot)
print('批量梯度下降算法')
print(weights1)
plotBestFit(X,Y,weights1)
y_hat1=predict(weights1,X)
print()
# 随机批量梯度下降算法
weights2=SoftmaxSGD(X, Y_onehot)
print('随机批量梯度下降算法')
print(weights2)
plotBestFit(X,Y,weights2)
y_hat2=predict(weights2,X)
参考资料：
1、Machine-Learning-With-Python
2、《机器学习实战》Peter Harrington著
3、《机器学习》西瓜书，周志华著
4、 斯坦福大学公开课 ：机器学习课程
5、机器学习视频，邹博
6、吴恩达_UFLDL中文教程
7、python 实现 softmax分类器（MNIST数据集）
8、python机器学习案例系列教程——逻辑分类/逻辑回归LR/一般线性回归（softmax回归）
9、机器学习 —— 基础整理（五）线性回归；二项Logistic回归；Softmax回归及其梯度推导；广义线性模型
