代码来源： UFLDL-Tutorial-notes-code，感谢其分享！
本人工作只是对代码进行讲解，如有改动会单独指出。
文章目录
-
-
- 任务
- 代码
- 思路
- 解析
-
- 计算目标函数 J ( θ ) J(\theta) J(θ)
- 计算梯度
-
任务
完成多分类softmax的目标函数 J ( θ ) J(\theta) J(θ)和梯度 ∇ θ J ( θ ) \nabla_{\theta}J(\theta) ∇θJ(θ)的计算。
代码
softmax_regression_vec.m 代码如下：
function [f,g] = softmax_regression(theta, X,y)
%
% Arguments:
% theta - A vector containing the parameter values to optimize.
% In minFunc, theta is reshaped to a long vector. So we need to
% resize it to an n-by-(num_classes-1) matrix.
% Recall that we assume theta(:,num_classes) = 0.
%
% X - The examples stored in a matrix.
% X(i,j) is the i'th coordinate of the j'th example.
% y - The label for each example. y(j) is the j'th example's label.
%
m=size(X,2);
n=size(X,1);
% theta is a vector; need to reshape to n x num_classes.
theta=reshape(theta, n, []);
num_classes=size(theta,2)+1;
% initialize objective value and gradient.
f = 0;
g = zeros(size(theta));
%
% TODO: Compute the softmax objective function and gradient using vectorized code.
% Store the objective function value in 'f', and the gradient in 'g'.
% Before returning g, make sure you form it back into a vector with g=g(:);
%
%%% YOUR CODE HERE %%%
A = exp([theta' * X; zeros(1,m)]); % 计算θ^T * X
B = bsxfun(@rdivide, A, sum(A)); % 对θ^T * X 归一化
C = log(B); % 取对数
I = sub2ind(size(C),y,1:size(C,2)); % 获取label对应的位置的值
f = (-1) * sum(C(I)); % 计算目标函数
% calculate g
Y = repmat(y',1,num_classes);
for i=1:num_classes
Y(Y(:,i)~=i,i) = 0; % 第i列中，仅标签为i的样本对应的位置不为0。
end
Y(Y~=0)=1; % 第i列非零元所在的位置，表示标签为i的元素的索引
disp(size(Y(:,1:(size(Y,2)-1))))
disp(size(B(1:(size(B,1)-1),:)'))
g = (-1) * X * (Y(:,1:(size(Y,2)-1))-B(1:(size(B,1)-1),:)');
%%% 别人的写法，两种写法效果一样,主要是稀疏矩阵生成不一样一点，他的速度略快%%
%%% 因为这里num_classes还很小，我耗时0.014272秒，他的耗时0.014249秒 %%%
% h = theta'*X;%h(k,i)第k个theta，第i个样本
% a = exp(h);
% a = [a;ones(1,size(a,2))];%加1行
% p = bsxfun(@rdivide,a,sum(a));
% c = log2(p);
% i = sub2ind(size(c), y,[1:size(c,2)]);
% values = c(i);
% f = -sum(values);
% d = full(sparse(1:m,y,1));
% d = d(:,1:(size(d,2)-1));%减1列
% p = p(1:(size(p,1)-1),:);%减1行
% g = X*(p'.-d);
g=g(:); % make gradient a vector for minFunc
思路
主要需要的就是这两个公式：
-
损失函数 J ( θ ) J(\theta) J(θ)
J ( θ ) = − [ ∑ i = 1 m ∑ k = 1 K 1 { y ( i ) = k } log exp ( θ ( k ) ⊤ x ( i ) ) ∑ j = 1 K exp ( θ ( j ) ⊤ x ( i ) ) ] = − [ ∑ i = 1 m ( 1 − y ( i ) ) log ( 1 − h θ ( x ( i ) ) ) + y ( i ) log h θ ( x ( i ) ) ] = − [ ∑ i = 1 m ∑ k = 0 1 1 { y ( i ) = k } log P ( y ( i ) = k ∣ x ( i ) ; θ ) ] \begin{aligned} J(\theta) &= - \left[ \sum_{i=1}^{m} \sum_{k=1}^{K} 1\left\{y^{(i)} = k\right\} \log \frac{\exp(\theta^{(k)\top} x^{(i)})}{\sum_{j=1}^K \exp(\theta^{(j)\top} x^{(i)})}\right] \\ &= - \left[ \sum_{i=1}^m (1-y^{(i)}) \log (1-h_\theta(x^{(i)})) + y^{(i)} \log h_\theta(x^{(i)}) \right] \\ &= - \left[ \sum_{i=1}^{m} \sum_{k=0}^{1} 1\left\{y^{(i)} = k\right\} \log P(y^{(i)} = k | x^{(i)} ; \theta) \right] \end{aligned} J(θ)=−[i=1∑mk=1∑K1{y(i)=k}log∑j=1Kexp(θ(j)⊤x(i))exp(θ(k)⊤x(i))]=−[i=1∑m(1−y(i))log(1−hθ(x(i)))+y(i)loghθ(x(i))]=−[i=1∑mk=0∑11{y(i)=k}logP(y(i)=k∣x(i);θ)]
-
梯度 ∇ θ J ( θ ) \nabla_{\theta}J(\theta) ∇θJ(θ)
∇ θ ( k ) J ( θ ) = − ∑ i = 1 m [ x ( i ) ( 1 { y ( i ) = k } − P ( y ( i ) = k ∣ x ( i ) ; θ ) ) ] \begin{aligned} \nabla_{\theta^{(k)}} J(\theta) = - \sum_{i=1}^{m}{ \left[ x^{(i)} \left( 1\{ y^{(i)} = k\} - P(y^{(i)} = k | x^{(i)}; \theta) \right) \right] } \end{aligned} ∇θ(k)J(θ)=−i=1∑m[x(i)(1{y(i)=k}−P(y(i)=k∣x(i);θ))]
但是有一点奇怪的地方：
在所给的框架下， θ \theta θ 的维度是 n × ( K − 1 ) n×(K-1) n×(K−1)，其中 K K K 为类别数，此处为10。
这样，我们的预测值 θ T X \theta^TX θTX 的维度就变成了 ( K − 1 ) × m (K-1)×m (K−1)×m，而不是 K × m K×m K×m，我个人认为还是挺奇怪的。
解决办法就是增加一个与 X X X 无关的偏置项（同样也不可优化），附加在 θ T X \theta^TX θTX 之后，从而将维度变为 K × m K×m K×m。
这样 n × ( K − 1 ) n×(K-1) n×(K−1) 维的参数确实也可以完成 K K K 类的分类，只不过有点难以理解。因为当 θ \theta θ 的每个元素都很小时，则偏置项相对较大，则被分类为第 K K K 类。
理解了这个之后，下面就很容易懂了。
解析
计算目标函数 J ( θ ) J(\theta) J(θ)
A = exp([theta' * X; zeros(1,m)]); % 计算θ^T * X
B = bsxfun(@rdivide, A, sum(A)); % 对θ^T * X 归一化
C = log(B); % 取对数
I = sub2ind(size(C),y,1:size(C,2)); % 获取label对应的位置的值
f = (-1) * sum(C(I)); % 计算目标函数
这部分对应目标函数的公式：
J ( θ ) = − [ ∑ i = 1 m ∑ k = 1 K 1 { y ( i ) = k } log exp ( θ ( k ) ⊤ x ( i ) ) ∑ j = 1 K exp ( θ ( j ) ⊤ x ( i ) ) ] = − [ ∑ i = 1 m ( 1 − y ( i ) ) log ( 1 − h θ ( x ( i ) ) ) + y ( i ) log h θ ( x ( i ) ) ] = − [ ∑ i = 1 m ∑ k = 0 1 1 { y ( i ) = k } log P ( y ( i ) = k ∣ x ( i ) ; θ ) ] \begin{aligned} J(\theta) &= - \left[ \sum_{i=1}^{m} \sum_{k=1}^{K} 1\left\{y^{(i)} = k\right\} \log \frac{\exp(\theta^{(k)\top} x^{(i)})}{\sum_{j=1}^K \exp(\theta^{(j)\top} x^{(i)})}\right] \\ &= - \left[ \sum_{i=1}^m (1-y^{(i)}) \log (1-h_\theta(x^{(i)})) + y^{(i)} \log h_\theta(x^{(i)}) \right] \\ &= - \left[ \sum_{i=1}^{m} \sum_{k=0}^{1} 1\left\{y^{(i)} = k\right\} \log P(y^{(i)} = k | x^{(i)} ; \theta) \right] \end{aligned} J(θ)=−[i=1∑mk=1∑K1{y(i)=k}log∑j=1Kexp(θ(j)⊤x(i))exp(θ(k)⊤x(i))]=−[i=1∑m(1−y(i))log(1−hθ(x(i)))+y(i)loghθ(x(i))]=−[i=1∑mk=0∑11{y(i)=k}logP(y(i)=k∣x(i);θ)]
说明比较重要的几点：
-
A = exp([theta' * X; zeros(1,m)]);
zeros(1,m)
A
-
I = sub2ind(size(C),y,1:size(C,2));
f = (-1) * sum(C(I));
完成的实际上是 1 { y ( i ) = k } 1\left\{y^{(i)} = k\right\} 1{y(i)=k} 这个函数。即筛选出预测类别中正确类别的概率，来计算目标函数。
其中前一句是得到预测类别中正确类别的概率所在的索引（
sub2ind
y
指明了获取的索引为预测类别中正确类别的概率。
后一句通过访问索引，得到 log P ( y ( i ) = k ∣ x ( i ) ; θ ) \log P(y^{(i)} = k | x^{(i)} ; \theta) logP(y(i)=k∣x(i);θ) ，计算目标函数。
变量的内容：
-
A
1.059044 1.016396 1.044481 1.094523 1.013291 …… 1.042245 1.009933 1.044206 1.097755 1.011094 …… 1.055637 1.012596 1.034192 1.084184 1.002852 …… 1.065303 1.014536 1.034789 1.09325 1.009064 …… 1.063292 1.014994 1.026204 1.085303 1.017293 …… 1.05112 1.016776 1.039049 1.088784 1.014322 …… 1.060112 1.019359 1.047857 1.10238 1.010763 …… 1.064268 1.006946 1.039117 1.097793 1.007203 …… 1.055498 1.020393 1.027037 1.076283 1.006618 …… 1 1 1 1 1 …… -
B
A
0.100703 0.100316 0.101044 0.101155 0.1004 …… 0.099106 0.099678 0.101017 0.101454 0.100183 …… 0.100379 0.099941 0.100048 0.100199 0.099366 …… 0.101298 0.100133 0.100106 0.101037 0.099982 …… 0.101107 0.100178 0.099276 0.100303 0.100797 …… 0.099949 0.100354 0.100518 0.100625 0.100503 …… 0.100804 0.100609 0.10137 0.101881 0.10015 …… 0.1012 0.099383 0.100525 0.101457 0.099797 …… 0.100366 0.100711 0.099356 0.099469 0.099739 …… 0.095088 0.098698 0.096741 0.092419 0.099083 …… -
C
B
-2.29558 -2.29943 -2.2922 -2.2911 -2.29859 …… -2.31157 -2.30581 -2.29247 -2.28815 -2.30076 …… -2.2988 -2.30317 -2.3021 -2.30059 -2.30894 …… -2.28969 -2.30126 -2.30153 -2.29227 -2.30277 …… -2.29158 -2.30081 -2.30986 -2.29956 -2.29465 …… -2.30309 -2.29905 -2.29742 -2.29636 -2.29757 …… -2.29457 -2.29652 -2.28898 -2.28395 -2.30109 …… -2.29066 -2.30877 -2.29735 -2.28812 -2.30462 …… -2.29893 -2.2955 -2.30905 -2.30791 -2.3052 …… -2.35295 -2.31569 -2.33572 -2.38142 -2.31179 …… -
I
6 14 21 37 44 ……
计算梯度
% calculate g
Y = repmat(y',1,num_classes);
for i=1:num_classes
Y(Y(:,i)~=i,i) = 0; % 第i列中，仅标签为i的样本对应的位置不为0。
end
Y(Y~=0)=1; % 第i列非零元所在的位置，表示标签为i的元素的索引
g = (-1) * X * (Y(:,1:(size(Y,2)-1))-B(1:(size(B,1)-1),:)');
对应公式：
∇ θ ( k ) J ( θ ) = − ∑ i = 1 m [ x ( i ) ( 1 { y ( i ) = k } − P ( y ( i ) = k ∣ x ( i ) ; θ ) ) ] \begin{aligned} \nabla_{\theta^{(k)}} J(\theta) = - \sum_{i=1}^{m}{ \left[ x^{(i)} \left( 1\{ y^{(i)} = k\} - P(y^{(i)} = k | x^{(i)}; \theta) \right) \right] } \end{aligned} ∇θ(k)J(θ)=−i=1∑m[x(i)(1{y(i)=k}−P(y(i)=k∣x(i);θ))]
写成向量格式：
∇ θ J ( θ ) = − [ X ( 1 { Y = K } − P ( Y = K ∣ X ; θ ) ) ] \begin{aligned} \nabla_{\theta} J(\theta) = - { \left[ X \left( 1\{ Y = K\} - P(Y = K | X; \theta) \right) \right] } \end{aligned} ∇θJ(θ)=−[X(1{Y=K}−P(Y=K∣X;θ))]
尤其注意一下几个变量的维度和内容：
-
X
-
Y
Y(:,1:(size(Y,2)-1)
…… …… 1 …… …… …… …… …… …… …… 1 …… …… …… 1 …… …… …… …… …… -
B
B(1:(size(B,1)-1),:)'
![吴恩达 coursera ML 第七课总结+作业答案前言目录正文模型表示作业答案[图]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)