Re: numerical differentiation?

[Gerard Jungman <jungman at lanl dot gov>]

Dave Morrison wrote:
> 
> I had in mind the first item of your list - calculating derivatives of
> functions. 

Ok. That's good. We need that.

> And one thing that would
> make the `multimin' code easier to use would be the option of using an
> automatically calculated numerical gradient instead of requiring a
> user-supplied routine to return the gradient.

Yes, that's important.
 
> Anyway, back to derivatives.  What I'd really like to see are routines
> to calculate gradients for N-dimensional functions.  With those, and a
> few new multimin routines, I can see how to provide functionality that
> I think would benfit GSL.

I don't think we should introduce any new multimin functions. Brian
should look at this as well, but I think we can just use the current
interfaces. You can wrap the numerical differentiation inside
the "df" part of a gsl_function_fdf (or gsl_multimin_function_fdf,
whatever... Brian: why do we need gsl_multimin_function_??? as well);
so it's a client-side issue. The multimin interfaces seem complete.


> I don't think the point you raise about the domain for
> `gsl_diff_function' should be a problem.  The best step size to use
> when calculating the derivative using a divided difference depends on
> details of the function but it's typically O(GSL_SQRT_DBL_EPSILON), so
> there's no need to wander very far away from the point at which you
> want to evaluate the derivative.

Well, as Brian will affirm, I am big on correctness. I think
we do need it. It should not complicate the implementation
much. A related issue is what to do for something like
f(x) = x^(3/2), if somebody asks for f'(0). I would
like to see it return 0.0, without complaint. There
is some question about how this sort of thing should
be dealt with. I think we will want to have the notion
of one-sided derivatives for functions of one variable,
so there will have to be an interface of the form

typedef enum {
  GSL_DIFF_SYMM,
  GSL_DIFF_RIGHT,
  GSL_DIFF_LEFT
} gsl_diff_direction_t;

gsl_diff_function("function and domain", double x0, gsl_diff_direction_t
d);

Then gsl_diff_function(..., 0.0, GSL_DIFF_SYMM) for x^(3/2)
would be an error, but gsl_diff_function(..., 0.0, GSL_DIFF_RIGHT)
would produce 0.0.

Similarly, for functions of several variables, we ought
to support directional derivatives, as well as computation
of gradients, and there is similarly an issue of one-sided
directional computation there. That takes care of those
functions which have directional derivatives but are not C^1.
The implementation should be very easy; the one-dimensional
implementation could be used directly.

I know these things are more than what would be needed
just to get multimin with numerical gradients going.
If it gives you a headache, Brian and I can deal with
the extra stuff. I just want to make sure all of this
is thought of beforehand, so we don't have to undo
interfaces which are not sufficiently correct or
general.


-- 
G. Jungman