【发布时间】：2015-08-20 03:43:51
【问题描述】：
我试图实现以下等式：
在matlab中。解释一些符号df/dt^(1)_{i,j}应该是一个向量，z^{(2)}_{k2}是一个实数，a^{(2)}_{i,j}是一个实数，[t^{(2)}_{k2}]是一个向量，x_i是一个向量，t^{(1)}_{i,j}是一个向量.有关符号的更多澄清 cmets，请查看与此相关的 math.stackexchange question。此外，我尝试用 cmets 对代码进行大量注释，说明输入和输出应该是什么，以尽量减少对相关变量维度的混淆。
我确实有一个潜在的实现（我相信它是正确的），但有时 MATLAB 有一些很好的隐藏技巧，我想知道这是否是上述矢量化方程的一个很好的实现，或者是否有更好的实现。
目前这里是我的代码：
function [ dJ_dt1 ] = compute_t1_gradient(t1,x,y,f,z_l1,z_l2,a_l2,c,t2,lambda)
%compute_t1_gradient_loops - computes the t1 parameter of a 2 layer HBF
% Computes dJ_dt1 according to:
% dJ_dt1
% Input:
% t1 = centers (Dp x Dd x Np)
% x = data (D x 1)
% y = label (1 x 1)
% f = f(x) (1 x 1)
% z_l1 = inputs l2 (Np x Dd)
% z_l2 = inputs l1 (K2 x 1)
% a_l2 = activations l2 (Np x Dd)
% a_l3 = activations l3 (K2 x 1)
% c = weights (K2 x 1)
% t2 = centers (K1 x K2)
% lambda = reg param (1 x 1)
% mu_c = step size (1 x 1)
% Output:
% dJ_dt1 = gradient (Dp x Dd x Np)
[Dp, ~, ~] = size(t1);
[Np, Dd] = size(a_l2);
x_parts = reshape(x, [Dp, Np])'; % Np x Dp
K1 = Np * Dd;
a_l2_col_vec = reshape(a_l2', [K1, 1]); %K1 x 1
alpha = bsxfun(@minus, a_l2_col_vec, t2); %K1 x K2
c_z_l2 = (c .* exp(-z_l2))'; % 1 x K2
alpha = bsxfun(@times, c_z_l2, alpha); %K1 x K2
alpha = bsxfun(@times, reshape(exp(-z_l1'),[K1, 1]) , alpha);
alpha = sum(alpha, 2); %K1 x 1
xi_t1 = bsxfun(@minus, x_parts', permute(t1, [1,3,2]));
% alpha K1 x 1
% xi_t1 Dp x Np x Dd
dJ_dt1 = bsxfun(@minus, reshape(alpha,[Dd, Np]), permute(xi_t1, [3, 2, 1]));
dJ_dt1 = permute(dJ_dt1,[3,1,2]);
dJ_dt1 = -4*(y-f)*dJ_dt1;
dJ_dt1 = dJ_dt1 + lambda * 0; %TODO
end
实际上，此时我决定将上述函数再次实现为for循环。不幸的是，他们没有给出相同的答案，这让我怀疑上述是否正确。我将粘贴我想要/打算矢量化的 for 循环代码：
function [ dJ_dt1 ] = compute_t1_gradient_loops(t1,x,y,f,z_l1,z_l2,a_l2,c,t2)
%compute_t1_gradient_loops - computes the t1 parameter of a 2 layer HBF
% Computes t1 according to:
% t1 := t1 - mu_c * dJ/dt1
% Input:
% t1 = centers (Dp x Dd x Np)
% x = data (D x 1)
% y = label (1 x 1)
% f = f(x) (1 x 1)
% z_l1 = inputs l2 (Np x Dd)
% z_l2 = inputs l1 (K2 x 1)
% a_l2 = activations l2 (Np x Dd)
% a_l3 = activations l3 (K2 x 1)
% c = weights (K2 x 1)
% t2 = centers (K1 x K2)
% lambda = reg param (1 x 1)
% mu_c = step size (1 x 1)
% Output:
% dJ_dt1 = gradeint (Dp x Dd x Np)
[Dp, ~, ~] = size(t1); %(Dp x Dd x Np)
[Np, Dd] = size(a_l2);
K2 = length(c);
t2_tensor = reshape(t2, Dd, Np, K2);
x_parts = reshape(x, [Dp, Np]);
dJ_dt1 = zeros(Dp, Dd, Np);
for i=1:Dd
xi = x_parts(:,i);
for j=1:Np
t_l1_ij = t1(:,i,j);
a_l2_ij = a_l2(j, i);
z_l1_ij = z_l1(j,i);
alpha_ij = 0;
for k2=1:K2
t2_k2ij = t2_tensor(i,j,k2);
c_k2 = c(k2);
z_l2_k2 = z_l2(k2);
new_delta = c_k2*-1*exp(-z_l2_k2)*2*(a_l2_ij - t2_k2ij);
alpha_ij = alpha_ij + new_delta;
end
alpha_ij = 2*(y-f)*-1*exp(-z_l1_ij)*2*(xi - t_l1_ij);
dJ_dt1(:,i,j) = alpha_ij;
end
end
end
实际上，我什至用Andrew Ng suggests 的方式近似了导数，以检查梯度下降，例如：
为此，我什至为它编写了代码：
%% update t1 unit test
%% dimensions
Dp = 3;
Np = 4;
Dd = 2;
K2 = 5;
K1 = Dd * Np;
%% fake data & params
x = (1:Dp*Np)';
y = 3;
c = (1:K2)';
t2 = rand(K1, K2);
t1 = rand(Dp, Dd, Np);
lambda = 0;
mu_t1 = 1;
%% call f(x)
[f, z_l1, z_l2, a_l2, ~ ] = f_star(x,c,t1,t2,Np,Dp);
%% update gradient
dJ_dt1_ij_loops = compute_t1_gradient_loops(t1,x,y,f,z_l1,z_l2,a_l2,c,t2);
dJ_dt1 = compute_t1_gradient(t1,x,y,f,z_l1,z_l2,a_l2,c,t2,lambda);
eps = 1e-4;
e_111 = zeros( size(t1) );
e_111(1,1,1) = eps;
derivative = (J(y, x, c, t2, t1 + e_111, Np, Dp) - J(y, x, c, t2, t1 - e_111, Np, Dp) ) / (2*eps);
derivative
dJ_dt1_ij_loops(1,1,1)
dJ_dt1(1,1,1)
但似乎没有一个衍生物与“近似”的结果一致。一次运行的输出如下所示：
>> update_t1_gradient_unit_test
derivative =
0.0027
dJ_dt1_ij_loops
ans =
0.0177
dJ_dt1
ans =
-0.5182
>>
我不清楚是否有错误...似乎它几乎与循环匹配，但足够接近吗？
吴恩达确实说过：
但是，我没有看到 4 位有效数字表示同意！甚至不是相同的数量级:(我猜两者都是错误的，但我似乎无法理解为什么或在哪里/如何。
在相关说明中，我还要求检查我在顶部的导数是否实际上（数学上正确），因为此时我不确定哪一部分是错误的，哪一部分是正确的。问题的链接在这里：
更新：
我已经用循环实现了一个新版本的衍生产品，它几乎与我创建的一个小例子一致。
这是新的实现（在某处有错误...）：
function [ dJ_dt1 ] = compute_df_dt1_loops3(t1,x,z_l1,z_l2,a_l2,c,t2)
% Computes t1 according to:
% df/dt1
% Input:
% t1 = centers (Dp x Dd x Np)
% x = data (D x 1)
% z_l1 = inputs l2 (Np x Dd)
% z_l2 = inputs l1 (K2 x 1)
% a_l2 = activations l2 (Np x Dd)
% a_l3 = activations l3 (K2 x 1)
% c = weights (K2 x 1)
% t2 = centers (K1 x K2)
% Output:
% dJ_dt1 = gradeint (Dp x Dd x Np)
[Dp, Dd, Np] = size(t1); %(Dp x Dd x Np)
K2 = length(c);
x_parts = reshape(x, [Dp, Np]);
dJ_dt1 = zeros(Dp, Dd, Np);
for i=1:Np
xi_part = x_parts(:,i);
for j=1:Dd
z_l1_ij = z_l1(i,j);
a_l2_ij = a_l2(i,j);
t_l1_ij = t1(:,i,j);
alpha_ij = 0;
for k2=1:K2
ck2 = c(k2);
t2_k2 = t2(:, k2);
index = (i-1)*Dd + j;
t2_k2_ij = t2_k2(index);
z_l2_k2 = z_l2(k2);
new_delta = ck2*(exp(-z_l2_k2))*2*(a_l2_ij - t2_k2_ij);
alpha_ij = alpha_ij + new_delta;
end
alpha_ij = -1 * alpha_ij * exp(-z_l1_ij)*2*(xi_part - t_l1_ij);
dJ_dt1(:,i,j) = alpha_ij;
end
end
这是计算数值导数的代码（正确且按预期工作）：
function [ dJ_dt1_numerical ] = compute_numerical_derivatives( x, c, t1, t2, eps)
% Computes t1 according to:
% df/dt1 numerically
% Input:
% x = data (D x 1)
% c = weights (K2 x 1)
% t1 = centers (Dp x Dd x Np)
% t2 = centers (K1 x K2)
% Output:
% dJ_dt1 = gradeint (Dp x Dd x Np)
[Dp, Dd, Np] = size(t1);
dJ_dt1_numerical = zeros(Dp, Dd, Np);
for np=1:Np
for dd=1:Dd
for dp=1:Dp
e_dd_dp_np = zeros(Dp, Dd, Np);
e_dd_dp_np(dp,dd,np) = eps;
f_e1 = f_star_loops(x,c,t1+e_dd_dp_np,t2);
f_e2 = f_star_loops(x,c,t1-e_dd_dp_np,t2);
numerical_derivative = (f_e1 - f_e2)/(2*eps);
dJ_dt1_numerical(dp,dd,np) = numerical_derivative;
end
end
end
end
我将提供 f 的代码和我实际使用的数字，以防人们重现我的结果：
这是 f 所做的代码（这也是正确的并且按预期工作）：
function [ f, z_l1, z_l2, a_l2, a_l3 ] = f_star_loops( x, c, t1, t2)
%f_start - computes 2 layer HBF predictor
% Computes f^*(x) = sum_i c_i a^(3)_i
% Inputs:
% x = data point (D x 1)
% x = [x1, ..., x_np, ..., x_Np]
% c = weights (K2 x 1)
% t2 = centers (K1 x K2)
% t1 = centers (Dp x Dd x Np)
% Outputs:
% f = f^*(x) = sum_i c_i a^(3)_i
% a_l3 = activations l3 (K2 x 1)
% z_l2 = inputs l2 (K2 x 1)
% a_l2 = activations l2 (Np x Dd)
% z_l1 = inputs l1 (Np x Dd)
[Dp, Dd, Np] = size(t1);
z_l1 = zeros(Np, Dd);
a_l2 = zeros(Np, Dd);
x_parts = reshape(x, [Dp, Np]);
%% Compute components of 1st layer z_l1 and a_l1
for np=1:Np
x_np = x_parts(:,np);
t1_np = t1(:,:, np);
for dd=1:Dd
t1_np_dd = t1_np(:, dd);
z_l1_np_dd = norm(t1_np_dd - x_np, 2)^2;
a_l1_np_dd = exp(-z_l1_np_dd);
% a_l1_np_dd = -z_l1_np_dd;
% a_l1_np_dd = sin(-z_l1_np_dd);
% insert
a_l2(np, dd) = a_l1_np_dd;
z_l1(np, dd) = z_l1_np_dd;
end
end
%% Compute components of 2nd layer z_l2 and a_l2
K1 = Dd*Np;
K2 = length(c);
a_l2_vec = reshape(a_l2', [K1,1]);
z_l2 = zeros(K2, 1);
for k2=1:K2
t2_k2 = t2(:, k2); % K2 x 1
z_l2_k2 = norm(t2_k2 - a_l2_vec, 2)^2;
% insert
z_l2(k2) = z_l2_k2;
end
%% Output later 3rd layer
a_l3 = exp(-z_l2);
% a_l3 = -z_l2;
% a_l3 = sin(-z_l2);
f = c' * a_l3;
end
这是我用于测试的数据：
%% Test 1:
% dimensions
disp('>>>>>>++++======--------> update t1 unit test');
% fake data & params
x = (1:6)'/norm(1:6,2)
c = [29, 30, 31, 32]'
t2 = [(13:16)/norm((13:16),2); (17:20)/norm((17:20),2); (21:24)/norm((21:24),2); (25:28)/norm((25:28),2)]'
Dp = 3;
Dd = 2;
Np = 2;
t1 = zeros(Dp,Dd, Np); % (Dp, Dd, Np)
t1(:,:,1) = [(1:3)/norm((1:3),2); (4:6)/norm((4:6),2)]';
t1(:,:,2) = [(7:9)/norm((7:9),2); (10:12)/norm((10:12),2)]';
t1
% call f(x)
[f, z_l1, z_l2, a_l2, a_l3 ] = f_star_loops(x,c,t1,t2)
% gradient
df_dt1_loops = compute_df_dt1_loops3(t1,x,z_l1,z_l2,a_l2,c,t2);
df_dt1_loops2 = compute_df_dt1_loops3(t1,x,z_l1,z_l2,a_l2,c,t2);
eps = 1e-10;
dJ_dt1_numerical = compute_numerical_derivatives( x, c, t1, t2, eps);
disp('---- Derivatives ----');
for np=1:Np
np
dJ_dt1_numerical_np = dJ_dt1_numerical(:,:,np);
dJ_dt1_numerical_np
df_dt1_loops2_np = df_dt1_loops(:,:,np);
df_dt1_loops2_np
end
请注意，数值导数现在是正确的（我确信因为我比较了 mathematica 返回的匹配值，加上 f 已被调试，所以它可以按我的意愿工作）。
这是一个输出示例（数值导数的矩阵应该使用我的方程匹配导数的矩阵）：
---- Derivatives ----
np =
1
dJ_dt1_numerical_np =
7.4924 13.1801
14.9851 13.5230
22.4777 13.8660
df_dt1_loops2_np =
7.4925 5.0190
14.9851 6.2737
22.4776 7.5285
np =
2
dJ_dt1_numerical_np =
11.4395 13.3836
6.9008 6.6363
2.3621 -0.1108
df_dt1_loops2_np =
14.9346 13.3835
13.6943 6.6363
12.4540 -0.1108
【问题讨论】：
-
只是几个问题。为什么
reshape(x, [Dp, Np])'如果x是一个向量，为什么不直接用反向索引重塑它呢？为什么exp(-1 * z_l2)而不仅仅是exp(-z_l2)？ -
另外，您让我想知道，通过不逐步分配值，而是将
bsxfun调用相互嵌套，您对alpha的分配是否会明显更快。我知道这会变得不那么透明，而且你的问题是关于重组的，但它仍然让我感到疑惑。 -
感谢您的回答，您对嵌套有很好的看法。通过反转索引，我只是指使用
reshape(x,[Np Dp])而不是reshape(x, [Dp, Np])'。同样，不应该是性能问题，我只是好奇。 -
对不起，你是完全正确的:) 你的代码根本没有对我说“愚蠢”，所以并不奇怪。
-
现在更清楚了，谢谢。如果我是你，我会用粗体标注矢量（考虑到你有带有各种索引集的矢量和标量）。而且您的示例需要进行一些清理：一些术语已从您的公式中删除，而有些术语甚至从未存在过（例如对
a_l3的引用）。
标签： matlab machine-learning vectorization gradient-descent