我试图实现以下方程：
在MATLAB中。为了解释一些符号，df/dt^(1)_{i,j} 应该是一个向量，z^{(2)}_{k2} 是一个实数，a^{(2)}_{i,j} 是一个实数，[t^{(2)}_{k2}] 是一个向量，x_i 是一个向量，t^{(1)}_{i,j} 是一个向量。有关符号的更多澄清性评论，请查看相关的 math.stackexchange 问题。此外，我已经在代码中添加了大量注释，说明输入和输出应该是什么，以尽量减少对所讨论变量维度的混淆。
实际上，我确实有一个可能的实现（我认为是正确的），但有时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 DpK1 = Np * Dd;a_l2_col_vec = reshape(a_l2', [K1, 1]); %K1 x 1alpha = bsxfun(@minus, a_l2_col_vec, t2); %K1 x K2c_z_l2 = (c .* exp(-z_l2))'; % 1 x K2alpha = bsxfun(@times, c_z_l2, alpha); %K1 x K2alpha = bsxfun(@times, reshape(exp(-z_l1'),[K1, 1]) , alpha);alpha = sum(alpha, 2); %K1 x 1xi_t1 = bsxfun(@minus, x_parts', permute(t1, [1,3,2]));% alpha K1 x 1% xi_t1 Dp x Np x DddJ_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; %TODOend实际上，此时我决定再次使用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; endendend实际上，我甚至按照 Andrew Ng 建议 的方式来检查梯度下降方程的近似导数，如下所示：
为此，我甚至编写了相应的代码：
%% update t1 unit test%% dimensionsDp = 3;Np = 4;Dd = 2;K2 = 5;K1 = Dd * Np;%% fake data & paramsx = (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 gradientdJ_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);derivativedJ_dt1_ij_loops(1,1,1)dJ_dt1(1,1,1)但似乎两种导数都与“近似”导数不一致。一次运行的输出如下所示：
>> update_t1_gradient_unit_testderivative = 0.0027dJ_dt1_ij_loopsans = 0.0177dJ_dt1ans = -0.5182>> 这让我不清楚是否有错误…似乎它几乎与带循环的版本匹配，但这是否足够接近？
Andrew Ng 确实说：
然而，我没有看到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; endend这是计算数值导数的代码（这是正确的并且按预期工作）：
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 endendend我将提供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_l1for 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; endend%% Compute components of 2nd layer z_l2 and a_l2K1 = 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 layera_l3 = exp(-z_l2);% a_l3 = -z_l2;% a_l3 = sin(-z_l2);f = c' * a_l3;end这是我用于测试的数据：
%% Test 1: % dimensionsdisp('>>>>>>++++======--------> update t1 unit test');% fake data & paramsx = (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)% gradientdf_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_npend请注意，现在数值导数是正确的（我确定是因为我与Mathematica返回的值进行了比较，这些值是匹配的，而且f已经调试过，所以它按我希望的方式工作）。
这是一个输出的示例（其中数值导数的矩阵应该与使用我的方程的导数的矩阵匹配）：
---- Derivatives ----np = 1dJ_dt1_numerical_np = 7.4924 13.1801 14.9851 13.5230 22.4777 13.8660df_dt1_loops2_np = 7.4925 5.0190 14.9851 6.2737 22.4776 7.5285np = 2dJ_dt1_numerical_np = 11.4395 13.3836 6.9008 6.6363 2.3621 -0.1108df_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 K2x_parts = reshape(x, [Dp, Dd]); %Dp x DddJ_dt1 = zeros(Dp, Dd, Np); %Dp x Dd x Npfor 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; endendend需要注意的一些事情：
- 我将
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 K2x_parts = reshape(x, [Dp, Dd]); %Dp x Dd%reorder things to align for bsxfun latera_l2=a_l2'; %Dd x Np <-> i,jz_l1=z_l1'; %Dd x Np <-> i,jt2_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 efforttempvar_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 Nptempvar_lastterm = permute(tempvar_lastterm,[2 3 1]); %Dd x Np x Dp%put together what we havedJ_dt1 = bsxfun(@times,partialsum.*exp(-z_l1),tempvar_lastterm); %Dd x Np x DpdJ_dt1 = permute(dJ_dt1,[3 1 2]); %Dp x Dd x Np再次需要注意的一些事情：
- 我为纯粹依赖
k2的求和部分定义了一个临时变量，因为它在下一步中被使用了两次。 - 我还将净前置因子
-4附加到这个变量上，因为你只需要乘以K2次，而不是Dp*Dd*Np次，这对于大型矩阵可能有很大区别。 - 我的函数按原样计算
k2求和，通过将(a-t2)分成两个求和，请参见以(NB)结尾的注释。事实证明，对于大型矩阵（将你的2-3-4-5维度的测试案例乘以100），这种分离会显著加速。当然，如果K2远大于t2的内部维度，那么你就失去了这个技巧的好处。 - 我添加了求和的“原始”版本的注释以供完整性和测试之用。
- 最后，我们只是将导数的各个因子拼接在一起：求和、第二个指数和最后的项。请注意，如果你的最后一个项包含
x_j而不是x_i，那么维度需要相应地调整。
性能
我检查了循环版本和我的两个向量化版本的两个测试案例。首先，你的原始示例是
%% update t1 unit test%% dimensionsDp = 3;Np = 4;Dd = 2;K2 = 5;K1 = Dd * Np;%% fake data & paramsx = (1:Dp*Dd)';y = 3;c = (1:K2)';t2 = rand(K1, K2);t1 = rand(Dp, Dd, Np);%% update gradientdJ_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：
%% dimensionsDp = 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的精度。虽然你可以在有变量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 K2x_parts = reshape(x, [Dp, Np]); %Dp x Npt1 = permute(t1,[1 3 2]); %Dp x Np x Dda_l2=a_l2'; %Dd x Np <-> j,iz_l1=z_l1'; %Dd x Np <-> j,itempvar_k2 = -4*c.*exp(-z_l2); % K2 x 1partialsum = bsxfun(@minus,a_l2,t2_tensor); %Dd x Np x K2partialsum = permute(partialsum,[3 1 2]); %K2 x Dd x Nppartialsum = squeeze(sum(bsxfun(@times,tempvar_k2,partialsum),1)); %Dd x Nptempvar_lastterm = bsxfun(@minus,x_parts,t1); %Dp x Np x Ddtempvar_lastterm = permute(tempvar_lastterm,[3 2 1]); %Dd x Np x DpdJ_dt1 = bsxfun(@times,partialsum.*exp(-z_l1),tempvar_lastterm); %Dd x Np x Dptempout=tempvar_lastterm;dJ_dt1 = permute(dJ_dt1,[3 1 2]); %Dp x Dd x Np请注意，这几乎与原始向量化版本相同，只是x的维度发生了变化，并且一些索引已经进行了置换。
