How to use the Finite Difference Method to get the gradient?

Hi there,

I need to calculate the gradient (partial derivative) of a function. I found that Matlab has got a '.p' file/function called 'finitedifferences' to do this.

However, I am not sure how to use it. For example, I found that Matlab do it like,

[gradFd,~,~,numEvals] = finitedifferences(x,funfcn{3},[],[],[],f,[],[], ...
        1:numberOfVariables,finDiffOpts,sizes,gradFd,[],[],finDiffFlags,[],varargin{:});

but what is funfcn{3} here?

Suppose I have a simple function like,

function F = myfun(x)

F = sin(x) + 3;

in which, x is a vector contains 6 elements. Then how to use the finitedifferences to get the gradient w.r.t each of this 6 elements? Thanks very much!

wbr, Aaron

Answer 1

[Andrew Newell]

Here is an example of a function that uses gradest correctly:

function [F, g] = myfun(x)

fun = @(y) 2*sum(sin(y)) + 3;
F = fun(x);

g = gradest(fun,x);

The key differences are:

1. the name of the function fun is different from myfun
2. fun and its gradient are both evaluated at x.

Also, I have made it more efficient in a couple of other ways. Now you can test this directly using:

x = rand(4,1);
ganalytic = 2*cos(x)'; % This is the analytic gradient
[F,g] = myfun(x);
max(abs(g-ganalytic))

The last line tells you the maximum difference between the analytic and numerical gradients.

Answer 2

[Arnaud Miege]

Have you looked at the gradient function?

Arnaud

Answer 3

[Andrew Newell]

That looks like an awkward way of doing it. I recommend downloading Adaptive Robust Numerical Differentiation from the FEX.

EDIT: Note that this package has functions for calculating gradient and Hessian.

Answer 4

[Aaronne]

Hi Andrew,

Thanks very much. I believe it is a very useful tool!

I wrote a test function as:

% myfun.m

function F = myfun(x)

F = 2*sin(x) + 3;

And the test:

% myfunTest.m
clear all;

xi = [4 2 3 4];

gradvec = derivest(@(x) myfun(x), xi(1), 'd', 1)

Then I can get the derivative at xi(1). That works perfectly!

However, another problem is I will use it in my optimization case. For example, my 'myfun.m' will return both function value and gradient at the same time, something like:

% myfun.m

function [g, F] = myfun(x)

F = 2*sin(x) + 3;

% Get the gradient here...(How to?)

g = [g_1, g_2, g_3, g_4]

and my test is:

% myfunTest.m
clear all;

xi = [4 2 3 4]; % Initialization

options = optimset('LargeScale', 'off', 'Display', 'iter-detailed', 'TolX', 0.0001, 'TolFun', 0.0000001, 'GradObj', 'on');

xMin = fminunc(@(x) myfun(x), xi, options);

Then may I ask how to get the g_1, g_2, g_3, g_4 using derivest please? Use the inline mode? Please give me some tips. Thanks a lot and have a nice day.

wbr, Aaron

Answer 5

[Aaronne]

Sorry Arnaud,

I didn't see you post before. Thanks for your reply. Please read my last post. The same question is if I use 'gradient' then how to get the g_1, g_2, g_3, g_4 in the function call?

Thanks very much.

wbr, Aaron

[Arnaud Miege]

I think you need to pass a vector of values to the gradient function, I'm not sure this would work within a function:

xi = [4 2 3 4]

F = 2*sin(xi) + 3;

gradient(F)

[Aaronne]

Hi Arnaud,

I've tested using 'gradient'; however, the result seems wrong:

First-order

Iteration Func-count f(x) Step-size optimality

0 1 2.07362 0

Optimization completed: The final point is the initial point.

The first-order optimality measure, 0.000000e+00, is less than

options.TolFun = 1.000000e-07.

Optimization Metric Options

relative first-order optimality = 0.00e+00 TolFun = 1e-07 (selected)

I have attached the .m files I run:

1.

% myfunTest.m

clc; clf; clear all; close all;

xi = [4 2 3 4];

% gradvec = derivest(@(x) myfun(x), xi(1), 'd', 1)

options = optimset('LargeScale', 'off', 'Display', 'iter-detailed', ...

'TolX', 0.0001, 'TolFun', 0.0000001, 'GradObj', 'on');

xMin = fminunc(@(x) myfun(x), xi, options);

2.

% myfun.m

function [F, g] = myfun(x)

F = 2*sin(x(1)) + 2*sin(x(2)) + 2*sin(x(3)) + 2*sin(x(4)) + 3;

g = gradient(F);

[Arnaud Miege]

I think that's because your function myfun outputs a scalar and gradient needs a vector of values to work out what the gradient is at each point:

>> [F, g] = myfun([1 2 3 4])

F =

5.2702

g =

0

Answer 6

[Aaronne]

Sorry, the F should be

function F = myfun(x)

F = 2*sin(x(1)) + 2*sin(x(2)) + 2*sin(x(3)) + 2*sin(x(4)) + 3;

[Andrew Newell]

Generally your questions should be added to your original question or included as comments to one of the answers. This makes it more clear what is an answer and what is not.

[Aaronne]

Hi Andrew,

I am really a newbie to this forum (if it is). I am sorry that I am not familiar with the style. I always post 'answer' in normal forum.

I will do the 'comment' rather than 'answer' my own question :-) Sorry again.

[Andrew Newell]

Not a problem! You already format your code nicely, which really helps.

Answer 7

[Aaronne]

Now I have problems:

My function is:

% myfun.m

function [F, g] = myfun(x)

F = 2*sin(x(1)) + 2*sin(x(2)) + 2*sin(x(3)) + 2*sin(x(4)) + 3;

g1 = derivest(@(x) myfun(x), x(1), 'd', 1);
g2 = derivest(@(x) myfun(x), x(2), 'd', 1);
g3 = derivest(@(x) myfun(x), x(3), 'd', 1);
g4 = derivest(@(x) myfun(x), x(4), 'd', 1);

g = [g_1, g_2, g_3, g_4];

My test is:

% myfunTest.m

clc; clf; clear all; close all;

xi = [4 2 3 4];

% gradvec = derivest(@(x) myfun(x), xi(1), 'd', 1)

options = optimset('LargeScale', 'off', 'Display', 'iter-detailed', ...
    'TolX', 0.0001, 'TolFun', 0.0000001, 'GradObj', 'on');

xMin = fminunc(@(x) myfun(x), xi, options);

Then by running the optimization I got the errors:

??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N)
to change the limit.  Be aware that exceeding your available stack space can
crash MATLAB and/or your computer.

Error in ==> iscellstr

Caused by:
    Failure in initial user-supplied objective function evaluation. FMINUNC
    cannot continue.

I know it is a weird way to do this way. But thanks a lot for any help!

wbr, Aaron

[Andrew Newell]

You get a recursion error because MYFUN is calling itself instead of the function F. Instead of using derivest, you should use gradest from the same package (see the edit to my answer above).

[Aaronne]

Hi Andrew,

I still haven't found a clue how to do it. I have a simple test function:

1.

% myfunTest.m

clc; clf; clear all; close all;

xi = [4 2 3 4];

% gradvec = derivest(@(x) myfun(x), xi(1), 'd', 1)

options = optimset('LargeScale', 'off', 'Display', 'iter-detailed', ...

'TolX', 0.0001, 'TolFun', 0.0000001, 'GradObj', 'on');

xMin = fminunc(@(x) myfun(x), xi, options);

2.

% myfun.m

function [F, g] = myfun(x)

F = 2*sin(x(1)) + 2*sin(x(2)) + 2*sin(x(3)) + 2*sin(x(4)) + 3;

[grad, err] = gradest(@(x) sum(myfun(x)),[1 2 3 4]);

g = grad;

I still got the same error as before when I use derivest. Could you give more details about how to deal with this simple case please? Thanks a lot.