本文详细解析了多分类softmax回归的目标函数及梯度计算方法，通过向量化代码实现，介绍了如何利用矩阵运算高效计算目标函数值及其梯度，适用于大规模数据集的机器学习任务。
代码来源： UFLDL-Tutorial-notes-code，感谢其分享！
本人工作只是对代码进行讲解，如有改动会单独指出。
任务
完成多分类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)]);是计算 exp ( θ ⊤ X ) \exp(\theta^\top X) exp(θ⊤X)，其中zeros(1,m)就是上面提到的偏置项。A的维度是 [ ( K − 1 ) × n ] × [ n × m ] + 1 × m = [ ( K − 1 ) × m ] + [ 1 × m ] = K × m [(K-1)×n]×[n×m] + 1×m = [(K-1)×m] + [1×m] = K×m [(K−1)×n]×[n×m]+1×m=[(K−1)×m]+[1×m]=K×m，具体在此处为10×60000。 -
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： n × m n×m n×m = 758×60000。 -
Y：为 1 { Y = K } 1\{ Y = K\} 1{Y=K}，维度 n × ( K − 1 ) n×(K-1) n×(K−1)。对应代码中Y(:,1:(size(Y,2)-1)要去掉非参数的偏置项。0 0 0 …… 0 0 0 …… 1 0 0 …… 0 0 0 …… 0 0 0 …… 0 0 0 …… 0 0 0 …… 0 0 0 …… 0 0 0 …… 0 1 0 …… 0 0 0 …… 0 0 0 …… 0 0 1 …… …… …… …… …… -
B：样本被分类为各个类别的概率， K × n K×n K×n。对应代码中B(1:(size(B,1)-1),:)'要去掉非参数的偏置项，并转置。最后维度为 n × ( K − 1 ) n×(K-1) n×(K−1)。
扫一扫
2394
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
