[Bug-gsl] wrong formula for BFGS update in gsl_multimin?

[Martin Wiebusch]

I've had some trouble with non-converging minimizations using the
gsl_multimin_fdfminimizer_vector_bfgs minimizer, which finally led me to
look into the code. The comments in 'vector_bfgs.c' and 'vector_bfgs2.c'
read

/* This is the BFGS update: */
/* p' = g1 - A dx - B dg */
/* A = - (1+ dg.dg/dx.dg) B + dg.g/dx.dg */
/* B = dx.g/dx.dg */

where, as far as I can see, p' is the new search direction, g1 and g
both denote the new gradient, dx is the position change and dg the
change in the gradient after the line search.

All the literature I've found on the BFGS method defines the new search
direction as p' = B^(-1) g, where B^(-1) is the BFGS estimate of the
inverse Hessian, defined trough outer and inner products of dx and dg
and the inverse Hessian B_prev^(-1) of the _previous_ iteration. (See,
e.g. http://en.wikipedia.org/wiki/BFGS_method). Using these expressions,
I can only reproduce the formula from the comments for the special case
where B_prev^(-1) is the identity matrix. I've also compared, during an
actual minimization, the vectors B^(-1) g (with my own implementation of
the B^(-1) matrix update) and the vector p' from the formula in the
comments, and they are indeed different.

Is it possible that the BFGS update in GSL is incorrect?

Best regards,
Martin.