符号梯度与解析梯度大不相同

Question

我正在尝试模拟一个移动机器人网络，该网络使用人工势场来规划到共享目的地 xd 的移动。这是通过从符号表达式生成一系列 m 文件（每个机器人一个）来完成的，因为这似乎是计算时间和准确性方面的最佳方式。但是，我无法弄清楚我的梯度计算出了什么问题：正在计算的分析梯度似乎有问题，而数值梯度计算正确（参见下面发布的图片）。我写了一个下面列出的 MWE，它也出现了这个问题。我检查了代码的文件生成部分，它确实返回了具有正确梯度的正确函数文件。但我不明白为什么解析梯度和数值梯度如此不同。

（可以找到下图的放大版本here）

% create symbolic variables
xd = sym('xd',[1 2]);
x = sym('x',[2 2]);

% create a potential function and a gradient function for both (x,y) pairs
% in x
for i=1:size(x,1)

phi = norm(x(i,:)-xd)/norm(x(1,:)-x(2,:));          % potential field function

xvector = reshape(x.',1,size(x,1)*size(x,2));       % reshape x to allow for gradient computation
grad = gradient(phi,xvector(2*i-1:2*i));            % compute the gradient
gradx = grad(1);grady=grad(2);                      % split the gradient in two components

% create function file names
gradfun = strcat('GradTester',int2str(i),'.m');     
phifun = strcat('PotTester',int2str(i),'.m');       

% generate two output files
matlabFunction(gradx, grady,'file',gradfun,'outputs',{'gradx','grady'},'vars',{xvector, xd});
matlabFunction(phi,'file',phifun,'vars',{xvector, xd});

end

clear all               % make sure the workspace is empty: the functions are in the files

pause(0.1)              % ensure the function file has been generated before it is called

% these are later overwritten by a specific case, but they can be used for
% debugging
x = 0.5*rand(2);
xd = 0.5*rand(1,2);

% values for the Stackoverflow case
x = [0.0533    0.0023;
     0.4809    0.3875];
xd = [0.4087    0.4343];

xp = x;     % dummy variable to keep x intact

% compute potential field and gradient for both (x,y) pairs
for i=1:size(x,1)

    % create a grid centered on the selected (x,y) pair
    xGrid = (x(i,1)-0.1):0.005:(x(i,1)+0.1);
    yGrid = (x(i,2)-0.1):0.005:(x(i,2)+0.1);

    % preallocate the gradient and potential matrices
    gradx = zeros(length(xGrid),length(yGrid));
    grady = zeros(length(xGrid),length(yGrid));
    phi   = zeros(length(xGrid),length(yGrid));

    % generate appropriate function handles
    fun   = str2func(strcat('GradTester',int2str(i)));
    fun2  = str2func(strcat('PotTester',int2str(i)));

    % compute analytic gradient and potential for each position in the xGrid and
    % yGrid vectors
    for ii = 1:length(yGrid)
        for jj = 1:length(xGrid)

            xp(i,:) = [xGrid(ii) yGrid(jj)];                % select the position
            Xvec = reshape(xp.',1,size(x,1)*size(x,2));     % turn the input into a vector
            [gradx(ii,jj),grady(ii,jj)] = fun(Xvec,xd);     % compute gradients
            phi(jj,ii) = fun2(Xvec,xd);                     % compute potential value

        end
    end

    [FX,FY] = gradient(phi);                % compute the NUMERICAL gradient for comparison

    %scale the numerical gradient
    FX = FX/0.005;
    FY = FY/0.005;

    % plot analytic result
    subplot(2,2,2*i-1)
    hold all
    xlim([xGrid(1) xGrid(end)]);
    ylim([yGrid(1) yGrid(end)]);
    quiver(xGrid,yGrid,-gradx,-grady)
    contour(xGrid,yGrid,phi)
    title(strcat('Analytic result for position ',int2str(i)));
    xlabel('x');
    ylabel('y');

    subplot(2,2,2*i)
    hold all
    xlim([xGrid(1) xGrid(end)]);
    ylim([yGrid(1) yGrid(end)]);
    quiver(xGrid,yGrid,-FX,-FY)
    contour(xGrid,yGrid,phi)
    title(strcat('Numerical result for position ',int2str(i)));
    xlabel('x');
    ylabel('y');

end

在我的代码xd 中，我试图生成的势场由 (x,y) 位置定义。 x是维度N x 2的位置矩阵，其中第一列代表x1、x2等，第二列代表y1、y2等。 Xvec 只是将此向量重塑为 x1,y1,x2,y2,x3,y3 等等，因为我生成的 matlab 函数只接受向量输入。

机器人 i 的梯度是通过 w.r.t 的导数来计算的。 x_i 和 y_i，这两个分量一起产生一个单一的导数“向量”，如箭袋图中所示。导数应该看起来像this，我检查了 [gradx,grady] 的符号表达式确实看起来像在生成一个 m 文件之前。

Answer 1

为了解决问题中给出的特定问题，您实际上是在计算 phi，这意味着您执行 gradient(phi) 与符号梯度相比并没有给出正确的结果。我会试着解释一下。以下是您创建xGrid 和yGrid 的方法：

% create a grid centered on the selected (x,y) pair
xGrid = (x(i,1)-0.1):0.005:(x(i,1)+0.1);
yGrid = (x(i,2)-0.1):0.005:(x(i,2)+0.1);

但随后在for循环中，ii和jj的使用与phi(jj,ii)或gradx(ii,jj)类似，但对应的物理位置相同。这就是为什么你的结果不同的原因。您遇到的另一个问题是您错误地使用了gradient。 Matlab 假设[FX,FY]=gradient(phi) 意味着phi 是从phi=f(x,y) 计算出来的，其中x 和y 是使用meshgrid 创建的矩阵。您实际上将phi 的元素安排得与此不同，因此gradient(phi) 给出了错误的答案。在反转jj 和ii 以及不正确的渐变之间，错误被抵消了（我怀疑你在尝试phi(ii,jj) 之后尝试做phi(jj,ii) 并发现它不起作用）。

不管怎样，整理一下，在你创建xGrid和yGrid之后，把这个放进去：

[X,Y]=meshgrid(xGrid,yGrid);

然后在你加载fun和fun2之后把代码改成：

for ii = 1:length(xGrid) %// x loop
    for jj = 1:length(yGrid) %// y loop
        xp(i,:) = [X(ii,jj);Y(ii,jj)]; %// using X and Y not xGrid and yGrid
        Xvec = reshape(xp.',1,size(x,1)*size(x,2));
        [gradx(ii,jj),grady(ii,jj)] = fun(Xvec,xd);
        phi(ii,jj) = fun2(Xvec,xd);  
    end
end

[FX,FY] = gradient(phi,0.005); %// use the second argument of gradient to set spacing

subplot(2,2,2*i-1)
hold all
axis([min(X(:)) max(X(:)) min(Y(:)) max(Y(:))]) %// use axis rather than xlim/ylim
quiver(X,Y,gradx,grady)
contour(X,Y,phi)
title(strcat('Analytic result for position ',int2str(i)));
xlabel('x');
ylabel('y');

subplot(2,2,2*i)
hold all
axis([min(X(:)) max(X(:)) min(Y(:)) max(Y(:))])
quiver(X,Y,FX,FY)
contour(X,Y,phi)
title(strcat('Numerical result for position ',int2str(i)));
xlabel('x');
ylabel('y');

---

我还有一些关于您的代码的其他信息。我认为您的潜在功能定义不明确，这会导致各种问题。您在问题中说x 是一个 Nx2 矩阵，但您的潜在函数定义为

norm(x(i,:)-xd)/norm(x(1,:)-x(2,:));

这意味着如果N 是三个，您将拥有以下三个潜力：

norm(x(1,:)-xd)/norm(x(1,:)-x(2,:));
norm(x(2,:)-xd)/norm(x(1,:)-x(2,:));
norm(x(3,:)-xd)/norm(x(1,:)-x(2,:));

我认为第三个没有意义。我认为这可能会对渐变造成一些混淆。

另外，我不确定是否有理由在您的真实代码中创建 .m 文件函数，但它们对于您发布的代码不是必需的。