matlab - 使用softmax损失时，如何调试和矢量化径向基函数网络的偏导数？

Question

当人们使用标准交叉熵（softmax损失）时，我尝试运行（随机或批处理）梯度下降：



在将Radial Basis Function (RBF)网络用作模型时（如果需要，可以观看caltech here讲座），可以扩展为多类分类（只需将RBF网络的输出馈送到softmax层即可轻松扩展。请注意<通过将RBF网络的输出通过每个标签P(y=l|x)的softmax层来简单地计算cc>，如下所示：



其中，l为负责对标签\theta_l进行预测的参数建立索引。

在这方面，我想通过计算参数的导数来优化模型。回想一下，在径向基函数网络中要优化的参数是最后一层的权重l和第一层的中心c。我已经实现并调试了如何计算权重t的导数。该代码按预期工作，因为偏导数与数值导数匹配。您可以找到单元测试代码为here。

我还尝试编写关于中心实现导数的代码，但我似乎无法使导数的实现与数值导数相匹配。我正在尝试实现的相对于中心c的损耗J的导数公式如下：



其中，t_k对应于负责预测标签h_{\theta_l}的RBF的输出。实际上，l表示非常简单：



我的主要问题是计算相对于h_{\theta_l}的J导数（上面的等式）。为此，我实现了following function，可以在不向量化的情况下天真地计算它：

function [ dJ_dt ] = compute_dJ_dt(z,x,y,t,c)
%Computes dJ_dc
%   Input:
%       z = (K x 1)
%       x = data point (D, 1)
%       y = labels (1 x 1)
%       t = centers (D x K)
%       c = weights (K x L)
%   Output:
%       dJ_dc = (D x K)
[D,K] = size(t);
[~, L] = size(c);
dJ_dt = zeros(D, K);
for k=1:K
    dJ_dt_k = zeros(D, 1);
    for l=1:L
        c_l = c(:,l);
        dh_dt_l = compute_dh_dt(z,x,t,c_l); %(D x K)
        delta = (y==l);
        dJ_dt_k = dJ_dt_k + dh_dt_l(:,k) * delta;
    end
    dJ_dt(:,k) = -dJ_dt_k;
end
end

并且它与 the numerical derivatives code不匹配。

我尝试了其他方法来检查它是否有效，我将在此处进行解释。如果有人有其他想法，请随时与他们分享，我觉得我已经用尽了许多好的新想法来调试它。


首先，一个自然的好问题是，我试图实现的导数的数学推导正确吗？即使我没有明确地与某人核对过数学推导，我也非常确信它的正确性，因为模型中关于 t_k和 c的偏导数的推导是相同的，并且只更改符号 讨论有问题的参数。由于我已经实现了关于 t的导数，并且通过了我所有的派生测试，因此我认为关于 \theta的导数的推导或任何参数 c应该是正确的。可以在 math.stack exchange here中看到我对这个方程的推导。
一种选择可能是 t实际上没有实现我所期望的方程式。确实可能是这样，并检查我是否独立执行了更多的 vectorized version of that code来查看我是否确实在执行我在纸上记下的方程式。由于方程的两个版本输出相同的导数值，因此我非常确信它们正在计算，确实是我怀疑的方程（同样，如果有人能够进一步向量化该方程，那将非常棒！我添加了向量化如此琐碎，以至于看起来似乎没有那么有趣或没有太多性能提升，但是确实删除了一个for循环）。


由于我在纸上拥有的方程式（很有可能）是正确的，并且由于方程式的两个版本都输出相同的值，因此该方程式的实现似乎是正确的，因此使我得出结论，也许数值导数代码存在错误。


numerical derivative code非常简单，以至于很难检查到底是什么地方出了问题。我发生的唯一可能是错误的事情可能是我对 softmax cost J的实现是错误的，但是我对此表示高度怀疑，因为...我已经为此编写了单元测试！另外，我用它来检查关于 \theta的数值导数，并且总是通过 compute_dJ_dt的数值导数，因此我无法想象 c是错误的。
最后要检查的是， c的计算正确。我已经写了 units tests for dh_dt，并且由于它们在每次运行时都与它们对应的数值导数匹配，所以我怀疑代码是正确的。


在这一点上，我不是100％不确定我还能尝试些什么，希望某人有一个好主意，或者指出我正在做的愚蠢的事情？我不确定现在该怎么想。感谢您的帮助和时间社区！

Answer 1

这是一种反气候的解决方案，但是我想这是可以预期的，因为此代码似乎是由工作组件构建的，因此注定是一个愚蠢的小错误。错误是我上面粘贴的代码中的错误，我应该一直使用delta作为标签指示和该标签概率之间的差异，但是我忘了减去概率。所以上面的代码是：

    delta = (y==l);

应该是什么时候：

    prob_y_x_h_x = prob_y_x(h_x); % (L x 1)
    ind_y_l = (y==l);
    delta = ind_y_l - prob_y_x_h_x(l);

因此，固定代码现在通过了数值测试，其外观如下：

function [ dJ_dt ] = compute_dJ_dt(h_x,z,x,y,t,c)
%Computes dJ_dc
%   Input:
%       z = (K x 1)
%       x = data point (D, 1)
%       y = labels (1 x 1)
%       t = centers (D x K)
%       c = weights (K x L)
%   Output:
%       dJ_dc = (D x K)
[D,K] = size(t);
[~, L] = size(c);
dJ_dt = zeros(D, K);
for k=1:K
    dJ_dt_k = zeros(D, 1);
    for l=1:L
        c_l = c(:,l);
        dh_dt_l = compute_dh_dt(z,x,t,c_l); %(D x K)
        prob_y_x_h_x = prob_y_x(h_x); % (L x 1)
        ind_y_l = (y==l);
        delta = ind_y_l - prob_y_x_h_x(l);
        dJ_dt_k = dJ_dt_k + dh_dt_l(:,k) * delta;
    end
    dJ_dt(:,k) = -dJ_dt_k;
end
end

我仍然不知道如何进一步对上面的代码进行矢量化处理，因此，我仍然很高兴收到有关该部分问题的反馈！这是我到目前为止的向量化：

function [ dJ_dt ] = compute_dJ_dt_vec(h_x,z,x,y,t,c)
%Computes dJ_dc
%   Input:
%       z = (K x 1)
%       x = data point (D, 1)
%       y = labels (1 x 1)
%       t = centers (D x K)
%       c = weights (K x L)
%   Output:
%       dJ_dc = (D x K)
[D,K] = size(t);
[~, L] = size(c);
dJ_dt = zeros(D, K);
for l=1:L
    c_l = c(:,l);
    dh_dt = compute_dh_dt(z,x,t,c_l); %(D x K)
    ind_y_l = (y==l);
    prob_y_x_h_x = prob_y_x(h_x); % (L x 1)
    dJ_dh = repmat( ind_y_l - prob_y_x_h_x(l) , D, K); %(D x K)
    dJ_dt = dJ_dt + dJ_dh.*dh_dt;
end
dJ_dt = -dJ_dt;
end