Числовой градиент в Matlab - вопросы округления

Question

Я пытаюсь вычислить числовой субградиент выпуклой функции. Моим подопытным является функция Вольфа. Он не должен быть очень точным, поэтому я попробовал нормальный конечный дифференциал в обоих направлениях: (f(xh)-f(x+h))/2h. В коде:

delta = 1e-10;

subgradient = zeros(length(xToEvaluate),1);

for i = 1 : length(xToEvaluate)
     deltaX = xToEvaluate;                     

     deltaX(i) = xToEvaluate(i) + delta;
     f1 = funct( deltaX );

     deltaX(i) = xToEvaluate(i) - delta;
     f2 = funct( deltaX );      

    subgradient(i,1) = (f1 - f2) / (2 * delta);  
end

При точном минимуме функции, при (-1,0), я получаю некоторые вещи на величину 1e-7так отлично. Когда я перехожу к чему-то вроде (-1, 0,1) или (-1, 1e-6), я получаю субградиент со вторым компонентом около 16,

Я знаю, что низкие дельты могут привести к ошибкам округления, но это не улучшается, когда я увеличиваю дельту.

Моя вторая попытка была одномерным трафаретом с пятью пунктами, но даже с дельтами приблизительно 1e-3 странный 16 продолжает появляться...

delta = 1e-3;

subgradient = zeros(length(xToEvaluate),1);

for i = 1 : length(xToEvaluate)

     xPlusTwo = xToEvaluate;
     xPlusOne = xToEvaluate;
     xMinusTwo = xToEvaluate;
     xMinusOne = xToEvaluate;

     xPlusTwo(i) = xToEvaluate(i) + 2*delta;
     xPlusOne(i) = xToEvaluate(i) + delta;
     xMinusTwo(i) = xToEvaluate(i) - 2*delta;
     xMinusOne(i) = xToEvaluate(i) - delta;

     subgradient(i,1) = (-funct(xPlusTwo) + 8*funct(xPlusOne) - 8*funct(xMinusOne) + funct(xMinusTwo))  / (12*delta);  
end

Кто-нибудь понял, что это такое?

Answer 1

Если вы определите градиент этой функции Вулфа, вы получите:

if x<=0;
    dfx = 9 - 81*x.^8;
    dfy = 16*sign(y);
elseif x>=abs(y);
    dfx = 5*0.5./sqrt(9*x.^2 + 16*y.^2)*9*2.*x;
    dfy  = 5*0.5./sqrt(9*x.^2 + 16*y.^2)*16*2.*y;
else
    dfx = 9;
    dfy  = 16*sign(y);
end

Итак, как вы можете видеть, второй компонент градиента для x<=0 является 16*sign(y)таким образом, это ноль, когда y==0, +-16 иначе.

Кстати, это не похоже на точный минимум лежит на [-1 0], а скорее на [-0.7598 0]
= [-(1/9)^(1/8) 0]