梯度下降法实现softmax回归MATLAB程序

梯度下降法实现softmax回归MATLAB程序

版权声明：本文原创，转载须注明来源。

解决二分类问题时我们通常用Logistic回归，而解决多分类问题时若果用Logistic回归，则需要设计多个分类器，这是相当麻烦的事情。softmax回归可以看做是Logistic回归的普遍推广（Logistic回归可看成softmax回归在类别数为2时的特殊情况），在多分类问题上softmax回归是一个有效的工具。

关于softmax回归算法的理论知识可参考这两篇博文：http://deeplearning.stanford.edu/wiki/index.php/Softmax%E5%9B%9E%E5%BD%92 ；
http://blog.csdn.net/acdreamers/article/details/44663305 。

本文自编mysoftmax_gd函数用于实现梯度下降softmax回归，代码如下（链接：http://pan.baidu.com/s/1geF2WMJ 密码：9x3x）：

MATLAB程序代码：

function [theta,test_pre,rate] = mysoftmax_gd(X_test,X,label,lambda,alpha,MAX_ITR,varargin)
% 该函数用于实现梯度下降法softmax回归
% 调用方式：[theta,test_pre,rate] = mysoftmax_gd(X_test,X,label,lambda,alpha,MAX_ITR,varargin)
% X_test：测试输入数据
% X：训练输入数据，组织为m*p矩阵，m为案例个数，p为加上常数项之后的属性个数
% label：训练数据标签，组织为m*1向量（数值型）
% lambda：权重衰减参数weight decay parameter
% alpha：梯度下降学习速率
% MAX_ITR：最大迭代次数
% varargin：可选参数，输入初始迭代的theta系数，若不输入，则默认随机选取
% theta：梯度下降法的theta系数寻优结果
% test_pre：测试数据预测标签
% rata：训练数据回判正确率

% Genlovy Hoo，2016.06.29. genlovhyy@163.com
%% 梯度下降寻优
Nin=length(varargin);
if Nin>1
    error('输入太多参数') % 若可选输入参数超过1个，则报错
end
[m,p] = size(X);
numClasses = length(unique(label)); % 求取标签类别数
if Nin==0
    theta = 0.005*randn(p,numClasses); % 若没有输入可选参数，则随机初始化系数
else
    theta=varargin{1}; % 若有输入可选参数，则将其设定为初始theta系数
end
cost=zeros(MAX_ITR,1); % 用于追踪代价函数的值
for k=1:MAX_ITR
    [cost(k),grad] = softmax_cost_grad(X,label,lambda,theta); % 计算代价函数值和梯度
    theta=theta-alpha*grad; % 更新系数
end
%% 回判预测
[~,~,Probit] = softmax_cost_grad(X,label,lambda,theta);
[~,label_pre] = max(Probit,[],2);
index = find(label==label_pre); % 找出预测正确的样本的位置
rate = length(index)/m; % 计算预测精度
%% 绘制代价函数图
figure('Name','代价函数值变化图');
plot(0:MAX_ITR-1,cost)
xlabel('迭代次数'); ylabel('代价函数值')
title('代价函数值变化图');% 绘制代价函数值变化图
%% 测试数据预测
[mt,pt] = size(X_test);
Probit_t = zeros(mt,length(unique(label)));
for smpt = 1:mt
    Probit_t(smpt,:) = exp(X_test(smpt,:)*theta)/sum(exp(X_test(smpt,:)*theta));
end
[~,test_pre] = max(Probit_t,[],2);

function  [cost,thetagrad,P] = softmax_cost_grad(X,label,lambda,theta)
% 用于计算代价函数值及其梯度
% X：m*p输入矩阵，m为案例个数，p为加上常数项之后的属性个数
% label：m*1标签向量（数值型）
% lambda：权重衰减参数weight decay parameter
% theta：p*k系数矩阵，k为标签类别数
% cost：总代价函数值
% thetagrad：梯度矩阵
% P：m*k分类概率矩阵，P（i，j）表示第i个样本被判别为第j类的概率
m = size(X,1);
% 将每个标签扩展为一个k维横向量（k为标签类别数），若样本i属于第j类，则
% label_extend（i，j）= 1，否则label_extend（i，j）= 0。
label_extend = [full(sparse(label,1:length(label),1))]';
% 计算预测概率矩阵
P = zeros(m,size(label_extend,2));
for smp = 1:m
    P(smp,:) = exp(X(smp,:)*theta)/sum(exp(X(smp,:)*theta));
end
% 计算代价函数值
cost = -1/m*[label_extend(:)]'*log(P(:))+lambda/2*sum(theta(:).^2);
% 计算梯度
thetagrad = -1/m*X'*(label_extend-P)+lambda*theta;

clear
clc
close all
load fisheriris % MATLAB自带数据集
% 对标签重新编号并准备训练/测试数据集
index_train = [1:40,51:90,101:140];
index_test = [41:50,91:100,141:150];
species_train = species(index_train);
X=[ones(length(species_train),1),meas(index_train,:)];
label = zeros(size(species_train));
label(strcmp('setosa',species_train)) = 1;
label(strcmp('versicolor',species_train)) = 2;
label(strcmp('virginica',species_train)) = 3;
species_test = species(index_test);
X_test = [ones(length(species_test),1),meas(index_test,:)];
lambda = 0.004; % 权重衰减参数Weight decay parameter
alpha = 0.1; % 学习速率
MAX_ITR=500; % 最大迭代次数
[theta,test_pre,rate] = mysoftmax_gd(X_test,X,label,lambda,alpha,MAX_ITR)

clear
clc
close all
load MNISTdata % MNIST数据集
% 准备训练/测试数据集
label = labels(1:9000); % 训练集标签
X = [ones(length(label),1),[inputData(:,1:9000)]']; % 训练集输入数据
label_test = labels(9001:end); % 测试集标签
X_test = [ones(length(label_test),1),[inputData(:,9001:end)]']; % 测试输入数据
lambda = 0.004; % 权重衰减参数Weight decay parameter
alpha = 0.1; % 学习速率
MAX_ITR=100; % 最大迭代次数
[theta,test_pre,rate] = mysoftmax_gd(X_test,X,label,lambda,alpha,MAX_ITR)
index_t = find(label_test==test_pre); % 找出预测正确的样本的位置
rate_test = length(index_t)/length(label_test); % 计算预测精度

水平有限，敬请指正交流。genlovhyy@163.com 。

参考资料：
【1】：http://deeplearning.stanford.edu/wiki/index.php/Softmax%E5%9B%9E%E5%BD%92
【2】：http://blog.csdn.net/acdreamers/article/details/44663305