我试图实现以下等式：
在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循环来实现。不幸的是，他们没有给出相同的答案，这使我怀疑以上内容是否正确。我将粘贴我想要/打算向量化的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个有效数字一致！甚至没有一个相同的数量级:(我猜这两个都是错误的，但是我似乎无法理解为什么或在哪里/如何。
在相关说明中，我还要求检查我顶部的导数是否实际上（数学上正确），因为在这一点上我不确定哪一部分是错误的，哪一部分是正确的。问题的链接在这里：
https://math.stackexchange.com/questions/1386958/partial-derivative-of-recursive-exponential-fx-sumk-2-k-2-1c-k-2-e
更新：
我用循环实现了新版本的导数，它几乎与我创建的一个小示例一致。
这是新的实现（在某处存在错误...）：
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);代替Dp=size(t1,1)在循环版本中，您忘了在总和后保留
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倍，这对于大型矩阵可能会有所不同。我的函数按原样通过将
k2分为两个和来计算(a-t2)和，请参见以(NB)结尾的注释。事实证明，对于大型矩阵（将2-2-3-5暗调的漂亮测试用例乘以100），这种分离可显着提高速度。当然，如果K2远大于t2的内部尺寸，那么您将失去这一技巧。为了完整性和测试，我在注释中添加了总和的“原始”版本。
最后，我们仅将导数的因素拼凑起来：总和，第二个指数和最后一项。请注意，如果最后一项包含
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代表矢量化代码的“原始”版本。事实证明，所得的导数与循环版本和朴素的矢量化版本完全吻合，而它们与优化的矢量版本之间最大的2e-14差异。这意味着我们很好。接近机器精度的差异仅仅是由于以不同顺序执行计算这一事实。为了评估性能，我将原始测试用例的尺寸乘以100：
%% dimensions
Dp = 300;
Np = 400;
Dd = 200;
K2 = 500;
K1 = Dd * Np;
而且我还设置了变量，以在每次函数调用之前和之后检查
cputime（因为tic/toc仅测量挂钟时间）。循环，优化和“原始”向量版本的测量时间分别为23秒，2秒和4秒。另一方面，后两个导数之间的最大差现在为1.8e-5。当然，我们的测试数据是随机的，至少可以说这不是最佳条件数据。在实际的应用程序中，这种差异可能不会成为问题，但是您应始终注意精度的损失（在优化版本中，我们具体减去了两个可能较大的数字）。当然，您可以尝试将公式划分为各种计算依据，这可能是一种更有效的方法。这也可能全部取决于阵列的大小。
半分析检查
您提到您尝试从定义中估计出导数，基本上是使用对称导数。您没有得到期望的结果，可能是由于原始功能的缺点。但是，我也想在这里注意几件事。您的
epsilon版本与您最初的尝试不一致的事实可能是由于您最初尝试中的实现错误
公式中的错误，即它实际上与
J的导数不对应（我知道您正在尝试在math.SE上调试这种情况）计算对称导数的神秘
J函数中的错误，仅在您的问题中提到如果一切顺利，您仍然可能存在纯粹的数学分歧：使用的
epsilon=1e-4因子完全是任意的。当您以这种方式检查导数时，基本上可以在给定点周围线性化函数。如果函数在半径epsilon附近变化太大（即非线性程度太大），则与精确值相比，对称导数将不准确。在进行这些检查时，应注意在导数中使用足够小的参数：足够小以期望函数具有线性行为，但又要足够大以避免由1/epsilon因子引起的数值噪声。最后说明：您可能应该避免在matlab中命名变量
eps，因为这是一个内置函数，告诉您“机器epsilon”（请查看help eps），与数字1的精度相对应默认值（即没有输入参数）。如果您有一个名为1i的变量，则可以调用复杂单元i，但如果可能的话，避免使用内置名称可能更安全。更新了最终的矢量化版本以对应于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 - 如何在MATLAB中很好地向量化以下关于向量的偏导数？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/31843929/