您的位置:首页 > 其它

softmaxCost实现

2015-10-16 21:25 344 查看
本次练习是按照ufldl SoftMax进行实现的。

关于代价函数的由来可以看看《统计学习方法》里边的logistic回归那一章,因为logistic模型是概率模型,所以会使用似然函数,优化也就是最大化似然函数,代价便是-log似然函数,最后的目标也就是最小化代价函数了。

不过由代价函数到它的导数这一部分的推导我开始没想明白,所以下面对此做个记录,看看到底是如何求导的。

推导



J(θ)=−1m[∑i=1m∑j=1kI(y(i)=j)(θTjx(i)−log∑l=1keθTlx(i))]+λ2∑i=1k∑j=0nθ2ij=−1m{∑i=1m[−log∑l=1keθTlx(i)+∑j=1kI(y(i)=j)(θTjx(i))]}+λ2∑i=1kθ2i
J(\theta)\\
=-\frac{1}{m}[\sum_{i=1}^{m}\sum_{j=1}^{k}I({y}^{(i)}=j)(\theta_{j}^{T}x^{(i)}-log\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}})]+\frac{\lambda}{2}\sum_{i=1}^{k}\sum_{j=0}^{n}\theta_{ij}^{2} \\
=-\frac{1}{m}\{\sum_{i=1}^{m}[-log\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}}+\sum_{j=1}^{k}I({y}^{(i)}=j)(\theta_{j}^{T}x^{(i)})]\}+\frac{\lambda}{2}\sum_{i=1}^{k}\theta_{i}^{2}

这一步的转换说明如下,在∑kj=1I(y(i)=j)(log∑kl=1eθTlx(i))\sum_{j=1}^{k}I({y}^{(i)}=j)(log\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}})里边,log∑kl=1eθTlx(i)log\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}}是个常数项,如果没有I(y(i)=j)I({y}^{(i)}=j),那就相当于是k个log∑kl=1eθTlx(i)log\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}}相加了,可是呢,在k个j里边,只有1个j能使得I(y(i)=j)I({y}^{(i)}=j)为1,所以最后这个式子的值等于log∑kl=1eθTlx(i)log\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}}。



∂J(θ)∂θj=−1m[∑mi=1∂(−log∑kl=1eθTlx(i))∂θj+∑mi=1∑kj=1∂I(y(i)=j)(θTjx(i))∂θj]+λθi=−1m∑i=1m[−eθTjx(i)x(i)∑kl=1eθTlx(i)+I(y(i)=j)x(i)]+λθi=−1m∑i=1mx(i)[−eθTjx(i)∑kl=1eθTlx(i)+I(y(i)=j)]+λθi
\frac{\partial J(\theta)}{\partial \theta_{j}}\\
=-\frac{1}{m}[\frac{\sum_{i=1}^{m}\partial(-log\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}})
}{\partial \theta_{j}}+
\frac{\sum_{i=1}^{m}\sum_{j=1}^{k}\partial I({y}^{(i)}=j)(\theta_{j}^{T}x^{(i)})}
{\partial \theta_{j}}]+\lambda\theta_{i}\\
=-\frac{1}{m}\sum_{i=1}^{m}[-\frac{e^{\theta_{j}^{T}x^{(i)}}x^{(i)}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}}}
+
I({y}^{(i)}=j)x^{(i)}]+\lambda\theta_{i} \\
=-\frac{1}{m}\sum_{i=1}^{m}x^{(i)}[-\frac{e^{\theta_{j}^{T}x^{(i)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}}}
+
I({y}^{(i)}=j)]+\lambda\theta_{i}

上面的推导中要注意I(y(i)=j)I({y}^{(i)}=j)是个常数项,所以求导的时候不考虑它的。

练习中我们首先要实现的是softmaxCost函数。

内存不够的解决办法

要载入练习中的那个训练集,大概需要700MB大小的空间,所以要开启matlab的3GB开关。

在C:\Windows\System32下找到cmd.exe文件(或者用windows键+R,然后填cmd),点右键选择“以管理员身份打开”,

键入 BCDEdit /set increaseuserva 3072 回车后显示“操作成功完成”,这一步之后要重启计算机才能生效的!

softmaxCost实现

首先来看看函数说明

function [cost, grad] = softmaxCost(theta, numClasses, inputSize, lambda, data, labels)


theta也就是softmax的模型参数,根据教程,它是一个numClasses*inputSize大小的矩阵。

numClasses则是类别总数,inputSize即每个样本的特征大小。

lambda是权重衰减系数λ\lambda,在J(θ)J(\theta)里边负责惩罚过大的theta值,目的在于防止过拟合。

data是训练数据,N*M大小,N是inputSize,M是总样本数,所以data(:,i)是第i个样本。

labels是M*1的向量,与data一一对应。

练习里边已经提供的代码如下:

% Unroll the parameters from theta
theta = reshape(theta, numClasses, inputSize);

numCases = size(data, 2);

groundTruth = full(sparse(labels, 1:numCases, 1));
cost = 0;

thetagrad = zeros(numClasses, inputSize);


首先把theta数据从向量转化成矩阵形式,我挺好奇的,为什么最后要把数据都转化为一个向量呢?计算出thetagrad之后也是这样做的,说是为了方便,可还是要记录下numClasses,inputSize这些值,在使用的时候还要进行转换,那方便之处在哪呢?

numCases就是样本总数啦。

groudTruth是后边做向量化要用到的一个矩阵,这个具体在向量化的时候进行介绍。

cost和thetagrad是函数的两个返回值,cost对应J(θ)J(\theta),

thetagrad(i,:)对应∂Jθ∂θ(i)\frac{\partial{J\theta}}{\partial{\theta(i)}}。

接着便是我们要写的代码了。

首先我写了个不怎么做向量化的版本,也就是for循环版本,不怎么用matlab的应该会更喜欢用for循环来解决问题吧!

无向量化版本



cost分析

来分析下这个cost该怎么求解。

对于第i个样本,按照公式,因为1{y(i)=j}1\{y^{(i)}=j\}是指示函数,当且仅当y(i)==jy^{(i)}==j的时候函数才为1,其余为0,所以我们并不需要做j=1:k的这个循环,只需要计算j = y(i)y^{(i)}时的值即可。

N is inputSize

θTlx(i)\theta_{l}^{T}x^{(i)} = theta(l,:)*data(:,i), theta(l,:) is a 1xN vector, data(:,i) is a

Nx1 vector

∑kl=1eθTlx(i)\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}} = sum(exp(theta*data(:,i))

theta is a numClasses x N matrix

求导分析



j是类别,所以对k个类别都要分别求下导数的。

对于第i个样本,第j个类的导数需要求得就是data(:,i)’*((l == j) - expJ/expAll)

下面代码中绿色部分不是注释啊,是因为它不识别matlab里边的这个矩阵转置”’”符号所以出了点问题。

expL也就是θTlx(i)\theta_{l}^{T}x^{(i)}

expAll也就是∑kl=1eθTlx(i)\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}}

cost1 = 0;
for i = 1:numCases
l = labels(i);
expL = exp(theta(l,:)*data(:,i));
expAll = sum(exp(theta*data(:,i)));
cost1 = cost1 - log(expL/expAll);
for j = 1:numClasses
expJ = exp(theta(j,:)*data(:,i));
thetagrad(j,:) = thetagrad(j,:) - data(:,i)'*((l == j) - expJ/expAll);
end
end
cost1 = cost1/numCases;
cost = cost1 + sum(sum(theta.*theta))*lambda/2;

for j = 1:numClasses
thetagrad(j,:) = thetagrad(j,:)/numCases + lambda*theta(j,:);
end


梯度检验ok!!!

可是在matlab里边用for很慢的所以还是要写出向量化形式加快速度。

怎么改呢?感觉挺有难度的

向量化版本

Implement softmaxCost

上面链接里的step2对向量化要用到的东西做了点介绍。

首先是groundTruth矩阵,它是个什么东西啊?

groundTruth = full(sparse(labels, 1:numCases, 1));

The command M = sparse(r, c, v) creates a sparse matrix such that M(r(i), c(i)) = v(i) for all i.

groundTruth是个numClasses*numCases大小的矩阵,第i列标记第i个样本的label情况,这个矩阵可以帮我们解决I(y(i)=j)I({y}^{(i)}=j)的表示问题。

full则是保证groundTruth是用矩阵表示,而不是因为稀疏性用其它方式保存。

%vectorization
Mat = theta*data;
Mat = bsxfun(@minus, Mat, max(Mat,[], 1));
Mat = exp(Mat);
Mat = bsxfun(@rdivide, Mat, sum(Mat));

%grad
thetagrad = -(groundTruth - Mat)*data'/numCases + lambda*theta;

Mat = log(Mat);
cost = -sum(sum(groundTruth.*Mat))/numCases + sum(sum(theta.*theta))*lambda/2;


下面对代码做解读:

假设样本数为numCases = 100;

类别数为k。

第一句Mat = theta*data;

已知θTix(j)\theta_{i}^{T}x^{(j)} = theta(i,:)*data(:,j);

所以theta(i,:)*data=

θTix(1)\theta_{i}^{T}x^{(1)}θTix(2)\theta_{i}^{T}x^{(2)}θTix(100)\theta_{i}^{T}x^{(100)}
所以theta*data=

θT1x(1)\theta_{1}^{T}x^{(1)}θT1x(2)\theta_{1}^{T}x^{(2)}θT1x(100)\theta_{1}^{T}x^{(100)}
θT2x(1)\theta_{2}^{T}x^{(1)}θT2x(2)\theta_{2}^{T}x^{(2)}θT2x(100)\theta_{2}^{T}x^{(100)}
θTkx(1)\theta_{k}^{T}x^{(1)}θTkx(2)\theta_{k}^{T}x^{(2)}θTkx(100)\theta_{k}^{T}x^{(100)}

第二句Mat = bsxfun(@minus, Mat, max(Mat,[], 1))

这一句举个例子会更清楚。

>> a = [1 5 3; 6 4 2]

a =

1     5     3
6     4     2
%max(a, [], 1)和max(a)的效果是一样的,都是返回每一列的最大值,所以是个行向量
>> max(a, [], 1)

ans =

6     5     3
%bsxfun既是对a的每一行都分别进行一个minus减操作,减的对象就是max(a,[],1)
>> bsxfun(@minus, a, max(a,[], 1))

ans =

-5     0     0
0    -1    -1
>> bsxfun(@minus, a, 1)

ans =

0     4     2
5     3     1

>> bsxfun(@minus, a, [1 2 3])

ans =

0     3     0
5     2    -1


所以这一句的操作是找出每一列最大值然后对这列的所有元素都减去这个最大值。

为什么要这么做呢?

When the products θTix(i)\theta_i^T x^{(i)} are large, the exponential function eθTix(i)e^{\theta_i^T x^{(i)}} will become very large and possibly overflow. When this happens, you will not be able to compute your hypothesis. However, there is an easy solution - observe that we can multiply the top and bottom of the hypothesis by some constant without changing the output

Hence, to prevent overflow, simply subtract some large constant value from each of the θTjx(i)\theta_j^T x^{(i)} terms before computing the exponential. In practice, for each example, you can use the maximum of the θTjx(i)\theta_j^T x^{(i)} terms as the constant. Assuming you have a matrix M containing these terms such that M(r, c) is θTrx(c)\theta_r^T x^{(c)}, then you can use the following code to accomplish this

如上所述就是为了预防θTix(i)\theta_i^T x^{(i)}太大的情况啦!

第三句Mat = exp(Mat)

这句没得说,求完了θTix(i)\theta_i^T x^{(i)},自然要求eθTix(i)e^{\theta_i^T x^{(i)}}

matlab倒是很方便哦,exp函数参数放个矩阵,便可以对矩阵中的每个元素都进行exp操作了。

第四句Mat = bsxfun(@rdivide, Mat, sum(Mat))

sum(Mat)是对Mat的每一列进行求和,最后得到是个行向量。

这一句求的便是eθTjx(i)∑kl=1eθTlx(i)\frac{e^{\theta_{j}^{T}x^{(i)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(i)}}}

所以得到的Mat如下:

eθT1x(1)∑kl=1eθTlx(1)\frac{e^{\theta_{1}^{T}x^{(1)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(1)}}}eθT1x(2)∑kl=1eθTlx(2)\frac{e^{\theta_{1}^{T}x^{(2)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(2)}}}eθT1x(100)∑kl=1eθTlx(100)\frac{e^{\theta_{1}^{T}x^{(100)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(100)}}}
eθT2x(1)∑kl=1eθTlx(1)\frac{e^{\theta_{2}^{T}x^{(1)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(1)}}}eθT2x(2)∑kl=1eθTlx(2)\frac{e^{\theta_{2}^{T}x^{(2)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(2)}}}eθT2x(100)∑kl=1eθTlx(100)\frac{e^{\theta_{2}^{T}x^{(100)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(100)}}}
eθTkx(1)∑kl=1eθTlx(1)\frac{e^{\theta_{k}^{T}x^{(1)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(1)}}}eθTkx(2)∑kl=1eθTlx(2)\frac{e^{\theta_{k}^{T}x^{(2)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(2)}}}eθTkx(100)∑kl=1eθTlx(100)\frac{e^{\theta_{k}^{T}x^{(100)}}}{\sum_{l=1}^{k}e^{\theta_{l}^{T}x^{(100)}}}

第五句求梯度

thetagrad = -(groundTruth - Mat)*data’/numCases + lambda*theta;



Mat便是p(y(i)=j|x(i);θ)∀iϵ[1,numCases],∀jϵ[1,k]p(y^{(i)}=j|x^{(i)};\theta) \forall i\epsilon [1,numCases],\forall j\epsilon [1,k]

groundTruth便是1{y(i)=j}∀iϵ[1,numCases],∀jϵ[1,k]1\{y^{(i)}=j\}\forall i\epsilon [1,numCases],\forall j\epsilon [1,k]

第六句Mat = log(Mat)

没什么好说的啦

第七句计算代价cost



cost = -sum(sum(groundTruth.*Mat))/numCases + sum(sum(theta.*theta))*lambda/2;

然后梯度检验ok!到此softmaxCost也就写完啦!

接下来完整的预测请看之后的啦!
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: