如何在 MATLAB 中很好地向量化以下关于向量的偏导数？
How does one vectorize nicely in MATLAB the following partial derivative with respect to a vector?
我试图实现以下等式：
在 MATLAB 中。解释一些记法df/dt^(1)_{i,j}应该是向量，z^{(2)}_{k2}是实数，a^{(2)}_{i,j}是实数，[t^{(2)}_{k2}]是向量，x_i是一个向量，t^{(1)}_{i,j} 是一个向量。有关符号的更多澄清注释，请参阅与此相关的 math.stackexchange question。此外，我试图通过对输入和输出应该是什么的注释来大量注释代码，以尽量减少对相关变量维度的混淆。
我确实有一个潜在的实现（我认为是正确的），但有时 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循环。不幸的是，他们没有给出相同的答案，这让我怀疑上面的内容是否正确。我将粘贴我 wanted/intended 进行矢量化的 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)
但似乎没有一个导数与 "approximated" 一致。一个 运行 的输出如下所示：
>> update_t1_gradient_unit_test
derivative =
0.0027
dJ_dt1_ij_loops
ans =
0.0177
dJ_dt1
ans =
-0.5182
>>
我不清楚是否有错误...似乎它与循环中的那个几乎匹配，但是足够接近吗？
Andrew Ng 说：
然而，我没有看到 4 位有效数字一致！甚至不在同一数量级:(我猜两者都是错误的，但我似乎无法理解为什么或 where/how.
在相关说明中，我还要求检查我在顶部的导数是否真的（数学上正确），因为在这一点上我不确定哪部分是错误的，哪部分是正确的.问题的 link 在这里：
更新:
我已经用循环实现了一个新版本的导数，它几乎与我创建的一个小例子一致。
这是新的实现（某处有错误...）：
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
更新：我对公式中某些数量的指数有一些误解，另请参阅更新的问题。我把原来的答案留在下面（因为矢量化应该以同样的方式进行），最后我添加了与 OP 的实际问题相对应的最终矢量化版本。
问题
您的代码和公式之间存在一些不一致。在您的公式中，您引用了 x_i，但是 x 数组的相应大小是索引 j 的大小。那么，这与您的 math.stackexchange 问题一致，其中 i 和 j 似乎与您在此处使用的符号互换...
无论如何，这是你函数的固定循环版本：
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
% Input:
% t1 = (Dp x Dd x Np)
% x = (D x 1)
% z_l1 = (Np x Dd)
% z_l2 = (K2 x 1)
% a_l2 = (Np x Dd)
% c = (K2 x 1)
% t2 = (K1 x K2)
%
% K1=Dd*Np
% D=Dp*Dd
% Dp,Np,Dd,K2 unique
%
% Output:
% dJ_dt1 = gradient (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); %Dd x Np x K2
x_parts = reshape(x, [Dp, Dd]); %Dp x Dd
dJ_dt1 = zeros(Dp, Dd, Np); %Dp x Dd x 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*exp(-z_l2_k2)*(a_l2_ij - t2_k2ij);
alpha_ij = alpha_ij + new_delta;
end
alpha_ij = -4*alpha_ij* exp(-z_l1_ij)*(xi - t_l1_ij);
dJ_dt1(:,i,j) = alpha_ij;
end
end
end
一些注意事项：
- 我将
x 的大小更改为 D=Dp*Dd 以保留公式的 i 索引。否则需要重新考虑更多的事情。
- 您可以使用
Dp=size(t1,1) 而不是 [Dp, ~, ~] = size(t1);
- 在您的循环版本中，您忘记在总和之后保留
alpha_ij，因为您用前置因子覆盖了旧值而不是将其相乘
如果我误解了你的意图，请告诉我，我会相应地更改循环版本。
矢量化版本
假设循环版本符合您的要求，这是一个矢量化版本，类似于您最初的尝试：
function [ dJ_dt1 ] = compute_t1_gradient_vect(t1,x,y,f,z_l1,z_l2,a_l2,c,t2)
%compute_t1_gradient_vect - computes the t1 parameter of a 2 layer HBF
% Input:
% t1 = (Dp x Dd x Np)
% x = (D x 1)
% y = (1 x 1)
% f = (1 x 1)
% z_l1 = (Np x Dd)
% z_l2 = (K2 x 1)
% a_l2 = (Np x Dd)
% c = (K2 x 1)
% t2 = (K1 x K2)
%
% K1=Dd*Np
% D=Dp*Dd
% Dp,Np,Dd,K2 unique
%
% Output:
% dJ_dt1 = gradient (Dp x Dd x Np)
Dp = size(t1,1);
[Np, Dd] = size(a_l2);
K2 = length(c);
t2_tensor = reshape(t2, Dd, Np, K2); %Dd x Np x K2
x_parts = reshape(x, [Dp, Dd]); %Dp x Dd
%reorder things to align for bsxfun later
a_l2=a_l2'; %Dd x Np <-> i,j
z_l1=z_l1'; %Dd x Np <-> i,j
t2_tensor = permute(t2_tensor,[3 1 2]); %K2 x Dd x Np
%the 1D part of the sum to be used in partialsum
%prefactors also put here to minimize computational effort
tempvar_k2 = -4*c.*exp(-z_l2); % K2 x 1
%compute sum(b(k)*(c-d(k)) as c*sum(b(k))-sum(b(k)*d(k)) (NB)
partialsum = a_l2*sum(tempvar_k2) ...
-squeeze(sum(bsxfun(@times,tempvar_k2,t2_tensor),1)); %Dd x Np
%alternative computation by definition:
%partialsum = bsxfun(@minus,a_l2,t2_tensor); %Dd x Np x K2
%partialsum = permute(partialsum,[3 1 2]); %K2 x Dd x Np
%partialsum = squeeze(sum(bsxfun(@times,tempvar_k2,partialsum),1)); %Dd x Np
%last part of the formula, (x-t1)
tempvar_lastterm = bsxfun(@minus,x_parts,t1); %Dp x Dd x Np
tempvar_lastterm = permute(tempvar_lastterm,[2 3 1]); %Dd x Np x Dp
%put together what we have
dJ_dt1 = bsxfun(@times,partialsum.*exp(-z_l1),tempvar_lastterm); %Dd x Np x Dp
dJ_dt1 = permute(dJ_dt1,[3 1 2]); %Dp x Dd x Np
再次提醒一下，一些注意事项：
- 我为总和的纯
k2 依赖部分定义了一个临时变量，因为它在下一步中使用了两次。
- 我还将净预因子
-4 附加到此变量，因为您只需要乘 K2 次而不是 Dp*Dd*Np 次，这对于大型矩阵可能会有所不同.
- 我的函数按原样通过将
(a-t2) 分成两个和来计算 k2 和，请参阅以 (NB) 结尾的注释。事实证明，对于大型矩阵（将 2-3-4-5 的漂亮测试用例乘以 100），这种分离会导致相当大的加速。当然，如果 K2 比 t2 的内部尺寸大得多，那么你就输了。
- 为了完整性和测试，我在评论中添加了总和的 "naive" 版本。
- 最后我们只是拼凑导数的因子：总和、第二个指数和最后一项。请注意，如果您的最终术语包含
x_j 而不是 x_i，则必须相应地调整维度。
性能
我检查了两个测试用例的循环版本和两个矢量化版本。首先，你原来的例子
%% update t1 unit test
%% dimensions
Dp = 3;
Np = 4;
Dd = 2;
K2 = 5;
K1 = Dd * Np;
%% fake data & params
x = (1:Dp*Dd)';
y = 3;
c = (1:K2)';
t2 = rand(K1, K2);
t1 = rand(Dp, Dd, Np);
%% update gradient
dJ_dt1_ij_loops = compute_t1_gradient_loops(t1,x,y,f,z_l1,z_l2,a_l2,c,t2);
dJ_dt1_vect = compute_t1_gradient_vect(t1,x,y,f,z_l1,z_l2,a_l2,c,t2);
dJ_dt1_vect2 = compute_t1_gradient_vect2(t1,x,y,f,z_l1,z_l2,a_l2,c,t2);
请注意，我再次更改了 x 的定义，..._vect2 代表向量化代码的 "naive" 版本。事实证明，对于循环版本和朴素矢量化，生成的导数与 完全一致 ，而它们与优化矢量版本之间存在最大 2e-14 差异。这意味着我们很好。接近机器精度的差异仅仅是因为计算是以不同的顺序执行的。
为了衡量性能，我将原始测试用例的维度乘以 100：
%% dimensions
Dp = 300;
Np = 400;
Dd = 200;
K2 = 500;
K1 = Dd * Np;
而且我还设置了变量以在每次函数调用前后检查 cputime（因为 tic/toc 仅测量挂钟时间）。循环、优化和 "naive" 矢量版本的测量时间分别为 23 秒、2 秒和 4 秒。另一方面，后两个导数之间的最大差异现在是 1.8e-5。当然，我们的测试数据是随机的，至少可以说不是最佳条件数据。可能在实际应用程序中，这种差异不会成为问题，但您应该始终小心精度损失（我们在优化版本中专门减去两个可能很大的数字）。
您当然可以尝试将您的公式划分为您计算它所依据的项，可能有一种更有效的方法。这也可能完全取决于数组的大小。
半解析检查
您提到您试图根据定义估计导数，主要是使用对称导数。你没有得到你所期望的，可能是因为你原来的功能有缺陷。不过，我也想在这里说明一些事情。您的 epsilon 版本与您的原始尝试不一致的事实可能是由于
- 您最初尝试中的实现错误
- 你的公式有误，即它实际上并不对应于
J 的导数（我知道你正试图在 math.SE 上调试这种情况）
- 你神秘的
J 函数计算你的对称导数时出错，这只在你的问题中提到
如果一切正常，您仍然可能有一个纯粹的数学分歧来源：您使用的 epsilon=1e-4 因子完全是任意的。当您以这种方式检查您的导数时，您基本上会围绕给定点线性化您的函数。如果您的函数在半径 epsilon 附近变化太大（即过于非线性），您的对称导数与精确值相比将不准确。在进行这些检查时，您应该小心地在导数中使用足够小的参数：小到足以预期函数的线性行为，但又大到足以避免 1/epsilon 因子引起的数值噪声。
最后说明：你应该避免在 matlab 中命名变量 eps，因为这是一个内置函数告诉你 "machine epsilon"（看看 help eps），默认情况下对应于数字 1 的精度（即没有输入参数）。如果您有一个名为 i 的变量，您可以调用复杂单元 1i，但尽可能避免使用内置名称可能更安全。
更新了最终矢量化版本以对应 OP 的更新问题：
function [ dJ_dt1 tempout] = compute_t1_gradient_vect(t1,x,z_l1,z_l2,a_l2,c,t2)
%compute_t1_gradient_vect - computes the t1 parameter of a 2 layer HBF
% Input:
% t1 = (Dp x Dd x Np)
% x = (D x 1)
% z_l1 = (Np x Dd)
% z_l2 = (K2 x 1)
% a_l2 = (Np x Dd)
% c = (K2 x 1)
% t2 = (K1 x K2)
%
% K1=Dd*Np
% D=Dp*Np
% Dp,Np,Dd,K2 unique
%
% Output:
% dJ_dt1 = gradient (Dp x Dd x Np)
Dp = size(t1,1);
[Np, Dd] = size(a_l2);
K2 = length(c);
t2_tensor = reshape(t2, Dd, Np, K2); %Dd x Np x K2
x_parts = reshape(x, [Dp, Np]); %Dp x Np
t1 = permute(t1,[1 3 2]); %Dp x Np x Dd
a_l2=a_l2'; %Dd x Np <-> j,i
z_l1=z_l1'; %Dd x Np <-> j,i
tempvar_k2 = -4*c.*exp(-z_l2); % K2 x 1
partialsum = bsxfun(@minus,a_l2,t2_tensor); %Dd x Np x K2
partialsum = permute(partialsum,[3 1 2]); %K2 x Dd x Np
partialsum = squeeze(sum(bsxfun(@times,tempvar_k2,partialsum),1)); %Dd x Np
tempvar_lastterm = bsxfun(@minus,x_parts,t1); %Dp x Np x Dd
tempvar_lastterm = permute(tempvar_lastterm,[3 2 1]); %Dd x Np x Dp
dJ_dt1 = bsxfun(@times,partialsum.*exp(-z_l1),tempvar_lastterm); %Dd x Np x Dp
tempout=tempvar_lastterm;
dJ_dt1 = permute(dJ_dt1,[3 1 2]); %Dp x Dd x Np
请注意，这与原始矢量化版本几乎相同，只是 x 的维度发生了变化，并且一些索引已被置换。
我试图实现以下等式：
在 MATLAB 中。解释一些记法df/dt^(1)_{i,j}应该是向量，z^{(2)}_{k2}是实数，a^{(2)}_{i,j}是实数，[t^{(2)}_{k2}]是向量，x_i是一个向量，t^{(1)}_{i,j} 是一个向量。有关符号的更多澄清注释，请参阅与此相关的 math.stackexchange question。此外，我试图通过对输入和输出应该是什么的注释来大量注释代码，以尽量减少对相关变量维度的混淆。
我确实有一个潜在的实现（我认为是正确的），但有时 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循环。不幸的是，他们没有给出相同的答案，这让我怀疑上面的内容是否正确。我将粘贴我 wanted/intended 进行矢量化的 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)
但似乎没有一个导数与 "approximated" 一致。一个 运行 的输出如下所示：
>> update_t1_gradient_unit_test
derivative =
0.0027
dJ_dt1_ij_loops
ans =
0.0177
dJ_dt1
ans =
-0.5182
>>
我不清楚是否有错误...似乎它与循环中的那个几乎匹配，但是足够接近吗？
Andrew Ng 说：
然而，我没有看到 4 位有效数字一致！甚至不在同一数量级:(我猜两者都是错误的，但我似乎无法理解为什么或 where/how.
在相关说明中，我还要求检查我在顶部的导数是否真的（数学上正确），因为在这一点上我不确定哪部分是错误的，哪部分是正确的.问题的 link 在这里：
更新:
我已经用循环实现了一个新版本的导数，它几乎与我创建的一个小例子一致。
这是新的实现（某处有错误...）：
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
更新：我对公式中某些数量的指数有一些误解，另请参阅更新的问题。我把原来的答案留在下面（因为矢量化应该以同样的方式进行），最后我添加了与 OP 的实际问题相对应的最终矢量化版本。
问题
您的代码和公式之间存在一些不一致。在您的公式中，您引用了 x_i，但是 x 数组的相应大小是索引 j 的大小。那么，这与您的 math.stackexchange 问题一致，其中 i 和 j 似乎与您在此处使用的符号互换...
无论如何，这是你函数的固定循环版本：
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
% Input:
% t1 = (Dp x Dd x Np)
% x = (D x 1)
% z_l1 = (Np x Dd)
% z_l2 = (K2 x 1)
% a_l2 = (Np x Dd)
% c = (K2 x 1)
% t2 = (K1 x K2)
%
% K1=Dd*Np
% D=Dp*Dd
% Dp,Np,Dd,K2 unique
%
% Output:
% dJ_dt1 = gradient (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); %Dd x Np x K2
x_parts = reshape(x, [Dp, Dd]); %Dp x Dd
dJ_dt1 = zeros(Dp, Dd, Np); %Dp x Dd x 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*exp(-z_l2_k2)*(a_l2_ij - t2_k2ij);
alpha_ij = alpha_ij + new_delta;
end
alpha_ij = -4*alpha_ij* exp(-z_l1_ij)*(xi - t_l1_ij);
dJ_dt1(:,i,j) = alpha_ij;
end
end
end
一些注意事项：
- 我将
x的大小更改为D=Dp*Dd以保留公式的i索引。否则需要重新考虑更多的事情。 - 您可以使用
Dp=size(t1,1) 而不是 - 在您的循环版本中，您忘记在总和之后保留
alpha_ij，因为您用前置因子覆盖了旧值而不是将其相乘
[Dp, ~, ~] = size(t1);
如果我误解了你的意图，请告诉我，我会相应地更改循环版本。
矢量化版本
假设循环版本符合您的要求，这是一个矢量化版本，类似于您最初的尝试：
function [ dJ_dt1 ] = compute_t1_gradient_vect(t1,x,y,f,z_l1,z_l2,a_l2,c,t2)
%compute_t1_gradient_vect - computes the t1 parameter of a 2 layer HBF
% Input:
% t1 = (Dp x Dd x Np)
% x = (D x 1)
% y = (1 x 1)
% f = (1 x 1)
% z_l1 = (Np x Dd)
% z_l2 = (K2 x 1)
% a_l2 = (Np x Dd)
% c = (K2 x 1)
% t2 = (K1 x K2)
%
% K1=Dd*Np
% D=Dp*Dd
% Dp,Np,Dd,K2 unique
%
% Output:
% dJ_dt1 = gradient (Dp x Dd x Np)
Dp = size(t1,1);
[Np, Dd] = size(a_l2);
K2 = length(c);
t2_tensor = reshape(t2, Dd, Np, K2); %Dd x Np x K2
x_parts = reshape(x, [Dp, Dd]); %Dp x Dd
%reorder things to align for bsxfun later
a_l2=a_l2'; %Dd x Np <-> i,j
z_l1=z_l1'; %Dd x Np <-> i,j
t2_tensor = permute(t2_tensor,[3 1 2]); %K2 x Dd x Np
%the 1D part of the sum to be used in partialsum
%prefactors also put here to minimize computational effort
tempvar_k2 = -4*c.*exp(-z_l2); % K2 x 1
%compute sum(b(k)*(c-d(k)) as c*sum(b(k))-sum(b(k)*d(k)) (NB)
partialsum = a_l2*sum(tempvar_k2) ...
-squeeze(sum(bsxfun(@times,tempvar_k2,t2_tensor),1)); %Dd x Np
%alternative computation by definition:
%partialsum = bsxfun(@minus,a_l2,t2_tensor); %Dd x Np x K2
%partialsum = permute(partialsum,[3 1 2]); %K2 x Dd x Np
%partialsum = squeeze(sum(bsxfun(@times,tempvar_k2,partialsum),1)); %Dd x Np
%last part of the formula, (x-t1)
tempvar_lastterm = bsxfun(@minus,x_parts,t1); %Dp x Dd x Np
tempvar_lastterm = permute(tempvar_lastterm,[2 3 1]); %Dd x Np x Dp
%put together what we have
dJ_dt1 = bsxfun(@times,partialsum.*exp(-z_l1),tempvar_lastterm); %Dd x Np x Dp
dJ_dt1 = permute(dJ_dt1,[3 1 2]); %Dp x Dd x Np
再次提醒一下，一些注意事项：
- 我为总和的纯
k2依赖部分定义了一个临时变量，因为它在下一步中使用了两次。 - 我还将净预因子
-4附加到此变量，因为您只需要乘K2次而不是Dp*Dd*Np次，这对于大型矩阵可能会有所不同. - 我的函数按原样通过将
(a-t2)分成两个和来计算k2和，请参阅以(NB)结尾的注释。事实证明，对于大型矩阵（将 2-3-4-5 的漂亮测试用例乘以 100），这种分离会导致相当大的加速。当然，如果K2比t2的内部尺寸大得多，那么你就输了。 - 为了完整性和测试，我在评论中添加了总和的 "naive" 版本。
- 最后我们只是拼凑导数的因子：总和、第二个指数和最后一项。请注意，如果您的最终术语包含
x_j而不是x_i，则必须相应地调整维度。
性能
我检查了两个测试用例的循环版本和两个矢量化版本。首先，你原来的例子
%% update t1 unit test
%% dimensions
Dp = 3;
Np = 4;
Dd = 2;
K2 = 5;
K1 = Dd * Np;
%% fake data & params
x = (1:Dp*Dd)';
y = 3;
c = (1:K2)';
t2 = rand(K1, K2);
t1 = rand(Dp, Dd, Np);
%% update gradient
dJ_dt1_ij_loops = compute_t1_gradient_loops(t1,x,y,f,z_l1,z_l2,a_l2,c,t2);
dJ_dt1_vect = compute_t1_gradient_vect(t1,x,y,f,z_l1,z_l2,a_l2,c,t2);
dJ_dt1_vect2 = compute_t1_gradient_vect2(t1,x,y,f,z_l1,z_l2,a_l2,c,t2);
请注意，我再次更改了 x 的定义，..._vect2 代表向量化代码的 "naive" 版本。事实证明，对于循环版本和朴素矢量化，生成的导数与 完全一致 ，而它们与优化矢量版本之间存在最大 2e-14 差异。这意味着我们很好。接近机器精度的差异仅仅是因为计算是以不同的顺序执行的。
为了衡量性能，我将原始测试用例的维度乘以 100：
%% dimensions
Dp = 300;
Np = 400;
Dd = 200;
K2 = 500;
K1 = Dd * Np;
而且我还设置了变量以在每次函数调用前后检查 cputime（因为 tic/toc 仅测量挂钟时间）。循环、优化和 "naive" 矢量版本的测量时间分别为 23 秒、2 秒和 4 秒。另一方面，后两个导数之间的最大差异现在是 1.8e-5。当然，我们的测试数据是随机的，至少可以说不是最佳条件数据。可能在实际应用程序中，这种差异不会成为问题，但您应该始终小心精度损失（我们在优化版本中专门减去两个可能很大的数字）。
您当然可以尝试将您的公式划分为您计算它所依据的项，可能有一种更有效的方法。这也可能完全取决于数组的大小。
半解析检查
您提到您试图根据定义估计导数，主要是使用对称导数。你没有得到你所期望的，可能是因为你原来的功能有缺陷。不过，我也想在这里说明一些事情。您的 epsilon 版本与您的原始尝试不一致的事实可能是由于
- 您最初尝试中的实现错误
- 你的公式有误，即它实际上并不对应于
J的导数（我知道你正试图在 math.SE 上调试这种情况） - 你神秘的
J函数计算你的对称导数时出错，这只在你的问题中提到
如果一切正常，您仍然可能有一个纯粹的数学分歧来源：您使用的 epsilon=1e-4 因子完全是任意的。当您以这种方式检查您的导数时，您基本上会围绕给定点线性化您的函数。如果您的函数在半径 epsilon 附近变化太大（即过于非线性），您的对称导数与精确值相比将不准确。在进行这些检查时，您应该小心地在导数中使用足够小的参数：小到足以预期函数的线性行为，但又大到足以避免 1/epsilon 因子引起的数值噪声。
最后说明：你应该避免在 matlab 中命名变量 eps，因为这是一个内置函数告诉你 "machine epsilon"（看看 help eps），默认情况下对应于数字 1 的精度（即没有输入参数）。如果您有一个名为 i 的变量，您可以调用复杂单元 1i，但尽可能避免使用内置名称可能更安全。
更新了最终矢量化版本以对应 OP 的更新问题：
function [ dJ_dt1 tempout] = compute_t1_gradient_vect(t1,x,z_l1,z_l2,a_l2,c,t2)
%compute_t1_gradient_vect - computes the t1 parameter of a 2 layer HBF
% Input:
% t1 = (Dp x Dd x Np)
% x = (D x 1)
% z_l1 = (Np x Dd)
% z_l2 = (K2 x 1)
% a_l2 = (Np x Dd)
% c = (K2 x 1)
% t2 = (K1 x K2)
%
% K1=Dd*Np
% D=Dp*Np
% Dp,Np,Dd,K2 unique
%
% Output:
% dJ_dt1 = gradient (Dp x Dd x Np)
Dp = size(t1,1);
[Np, Dd] = size(a_l2);
K2 = length(c);
t2_tensor = reshape(t2, Dd, Np, K2); %Dd x Np x K2
x_parts = reshape(x, [Dp, Np]); %Dp x Np
t1 = permute(t1,[1 3 2]); %Dp x Np x Dd
a_l2=a_l2'; %Dd x Np <-> j,i
z_l1=z_l1'; %Dd x Np <-> j,i
tempvar_k2 = -4*c.*exp(-z_l2); % K2 x 1
partialsum = bsxfun(@minus,a_l2,t2_tensor); %Dd x Np x K2
partialsum = permute(partialsum,[3 1 2]); %K2 x Dd x Np
partialsum = squeeze(sum(bsxfun(@times,tempvar_k2,partialsum),1)); %Dd x Np
tempvar_lastterm = bsxfun(@minus,x_parts,t1); %Dp x Np x Dd
tempvar_lastterm = permute(tempvar_lastterm,[3 2 1]); %Dd x Np x Dp
dJ_dt1 = bsxfun(@times,partialsum.*exp(-z_l1),tempvar_lastterm); %Dd x Np x Dp
tempout=tempvar_lastterm;
dJ_dt1 = permute(dJ_dt1,[3 1 2]); %Dp x Dd x Np
请注意，这与原始矢量化版本几乎相同，只是 x 的维度发生了变化，并且一些索引已被置换。