Numerical Gradient in Matlab - Rounding Issues

Question

I'm trying to compute a numerical subgradient of a convex function. My test subject is the Wolfe function. It doesn't need to be super accurate, so I tried a normal finite differential in both directions: (f(x-h)-f(x+h))/2h. In code:

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

At the exact minimum of the function, at (-1 ,0), I get some things at the magnitude 1e-7, so perfectly fine. As I move to something like (-1, 0.1) or (-1, 1e-6), I get a subgradient with second component of about 16.

I'm aware that low deltas might introduce rounding errors, but it doesn't get better as I increase delta.

My second try was a one-dimensional five-point stencil, but even with deltas of around 1e-3 the weird 16 keeps popping up...

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

Anyone got an idea what this is all about?

It probably is because the gradient there _actually is `[-72 16]`_ — Gunther Struyf, Jul 23 '13 at 21:37

Gunther Struyf · Accepted Answer · 2013-08-07T14:19:09.030

0

If you work out the gradient of that Wolfe function, you come up with:

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

So as you can see, the second component of the gradient for x<=0 is 16*sign(y), thus it is zero when y==0, +-16 otherwise.

BTW, it doesn't look like the exact minimum lies at [-1 0], but rather at [-0.7598 0]
= [-(1/9)^(1/8) 0]

edited Aug 07 '13 at 14:19

answered Jul 23 '13 at 21:44

Gunther Struyf

11,158
2
34
58

Oh man, I totally derped on that... My mind is still set to continuus functions and differentials and said this can't be :O But I disagree on the minimum, [-1 0] has a function value of -8, [0.7598 0] has some positive value. – Florian Rohm Aug 07 '13 at 11:06
oops, forgot the negative sign. And I don't know how you calculate your functions, but `Wolfe(-1,0)=0` – Gunther Struyf Aug 07 '13 at 14:18
With x=-1, y=0 I'm in the first case, x \leq 0. Therefore Wolfe(-1,0)= 9(-1) + 16|0| - (-1)^9 = -9+1 = -8 – Florian Rohm Aug 07 '13 at 16:58
@Florian: that's not what [the link in your question](http://www.wolframalpha.com/input/?i=f%28x%2Cy%29+%3D+Piecewise%5B%7B%7B+9%2ax+%2B+16%2aabs%28y%29+-+9%2ax%5E9+%2C+x%3C%3D0%7D%2C%7B5%2asqrt%289%2ax%5E2+%2B+16+%2a+y%5E2%29%2C+x%3E%3Dabs%28y%29%7D%2C%7B9%2ax%2B16%2aabs%28y%29%2C+0%3Cx%3Cabs%28y%29%7D%7D%5D) says – Gunther Struyf Aug 07 '13 at 18:59
The code includes { 9*x + 16*abs(y) - 9*x^9 , x<=0} which is what I typed. Are you sure you didn't skid down to the alternate forms or the partial differential? – Florian Rohm Aug 08 '13 at 21:07
look at your formula in your previous comment, you're ommitting the `9*` before the `x^9` – Gunther Struyf Aug 08 '13 at 21:32

Numerical Gradient in Matlab - Rounding Issues

1 Answers1