[Help-gsl] Re: Numerical Gradient?

[Richard Henwood]

David Doria wrote:

> I have a function f() that is not analytic and therefore I need to find a
> numerical gradient.  However to initialize a gsl_multimin_function_fdf, I
> still need a function my_df that simply returns the value of the gradient
> at the specified point.
> 
> For testing, lets try it with a function that we COULD have taken a normal
> gradient of (this example ships with the package)
> 
> double my_f (const gsl_vector *v, void *params)
> {
>     double x, y;
>     double *dp = (double *)params;
>     x = gsl_vector_get(v, 0);
>     y = gsl_vector_get(v, 1);
>     return 10.0 * (x - dp[0]) * (x - dp[0]) +
>     20.0 * (y - dp[1]) * (y - dp[1]) + 30.0;
> }
> 
> and its gradient:
> 
> void my_df (const gsl_vector *v, void *params, gsl_vector *df)
> {
>     double x, y;
>     double *dp = (double *)params;
>     x = gsl_vector_get(v, 0);
>     y = gsl_vector_get(v, 1);
>     gsl_vector_set(df, 0, 20.0 * (x - dp[0]));
>     gsl_vector_set(df, 1, 40.0 * (y - dp[1]));
> }
> 
> but instead i would like: (and I imagine I should keep the same function
> form (ie. same parameters))
> 
> void my_df (const gsl_vector *v, void *params, gsl_vector *df)
> {
>  NUMERICAL GRADIENT OF f at v;
> }
> 
> Is there a function I can call that will do that?
> 

My experience may be applicable to gsl_multimin_function_fdf.

I have implemented fitting (gsl_multifit_fdfsolver) using numerical
differentials. However, I did not find a function which would generically
differentiate a given function with respect to a given parameter. I wrote
the function I wanted to differentiate out multiple times, each time
exposing the parameter I wanted. see below for a terse example of fitting a
log Normal Probability Density Function:
[NOTE: log normal is not a good case as the analytical partial diffentials
exist!]

double myLogNormalPDF_x ( double x, void * p) {...};
double myLogNormalPDF_mu ( double , void * p) {...};
double myLogNormalPDF_sigma ( double , void * p) {...};

and then: 

int myLogNormalPDF_f (const gsl_vector * fitParams, void * extraParams,      
    gsl_vector * fn) {
...
    for (i = 0; i < n; i++) {
        current_x = xdata[i];
        if (current_x > 0.0) {
            current_eval = myLogNormalPDF_x (current_x, &params);
            gsl_vector_set (fn, i, current_eval - ydata[i]);
        }
    }
    return GSL_SUCCESS;
}

and:

int myLogNormalPDFnum_df ( const gsl_vector * fitParams, void * fitData,
    gsl_matrix * J) {
...
    for (i = 0; i < n; i++) {
        current_x = xdata[i];
        if (current_x > 0) {
            params.x = current_x;
            params.mu = fn_mu;
            params.sigma = fn_sigma;

            F.function = &myLogNormalPDF_mu;
            gsl_diff_central(&F, fn_mu, &df_dmu, &difAbserr);
            gsl_matrix_set(J, i, 0, df_dmu);

            F.function = &myLogNormalPDF_sigma;
            gsl_diff_central(&F, fn_sigma, &df_dsigma, &difAbserr);
            gsl_matrix_set(J, i, 1, df_dsigma);
        }
        else {
            gsl_matrix_set(J, i, 0, 0.0);
            gsl_matrix_set(J, i, 1, 0.0);
        }
    }
    return GSL_SUCCESS;
}