[Help-gsl] Comments/questions about BFGS algorithm implemented for the GSL

[Alan W. Irwin]

Background: I have a GPLed fortran library (see freeeos.sf.net) that
requires a routine to minimize the thermodynamic free energy in order to
calculate the equation of state (pressure as a function of density and
temperature) for the physical conditions in the interiors of stars. I have
been experimenting with the BFGS technique to accomplish this task since
that technique comes highly recommended in Fletcher's book, "Practical
Methods of Optimization", 2nd edition.

For months I have been experimenting with the L-BFGS-B fortran
implementation (http://www.ece.northwestern.edu/~nocedal/lbfgsb.html) to
minimize the free-energy. The L-BFGS-B implementation is recommended
everywhere, and it usually works for me, but there are some problems.

(1) The free license is non-standard and has an obnoxious advertising clause
that means not only my scientific paper (which is a reasonable requirement
which I would do out of professional courtesy in any case), but also
unreasonably all papers that ever use my equation of state calculation must
include two references specified by the original authors of the L-BFGS-B
fortran implementation.

(2) the authors deliberately do not maintain the code now and instead prefer
users (according to recent e-mail from Nocedal to me) to switch to
proprietary (and quite costly) software for minimization. Presumably because
of the bad free license, nobody else has stepped forward to maintain the 
L-BFGS-B code.

(3) From my many experiments for a wide range of physical conditions,
L-BFGS-B usually does the job, but I have found two bugs:

  (a) the code sometimes asks for x values substantially outside the bounds
  specified by the user.  Fortunately, I was able to transform my problem to
  one which did not require bounds, but for general use this is a major
  deficiency of the L-BFGS-B implementation since its major claim to
  fame is the ability to specify simple bounds.

  (b) Even for the transformed unbounded case, I found one case in my
  experiments where the code ignored a good minimum solution found during
  the line minimization and demanded a bad solution should be used instead

Because of these licensing and non-robustness issues I have now given up
completely on the L-BFGS-B fortran implementation, and I have now
implemented a fortran alternative which follows exactly the BFGS technique
implemented (in C) in the GSL.  It is early days yet, and I have only tested
it on one standard problem (Rosenbrock's function), but for that my fortran
implementation gives essentially identical results to using the GSL.

Here are the results for the iteration number, the two components of x -
known minimum, f, and the two components of the gradient for the GSL case.
(My own fortran implementation provides essentially identical results).

    1    -2.03710e+00     8.38208e-02     4.15657e+00
   -6.49952e-01
    1.65088e+00
    2    -1.95205e+00    -1.28596e-01     3.93289e+00
   -1.72279e+01
   -6.99747e+00
    3    -1.95496e+00    -1.21321e-01     3.93253e+00
   -1.66174e+01
   -6.65343e+00
    4    -1.94768e+00    -1.18408e-01     3.82074e+00
   -1.01548e+01
   -3.30248e+00

...

  116    -4.77646e-05    -9.57178e-05     2.28510e-09
   -1.91717e-05
   -3.81806e-05
  117    -4.01129e-05    -8.03843e-05     1.61161e-09
   -1.61707e-05
   -3.20288e-05
  118    -2.48095e-05    -4.97174e-05     6.16493e-10
   -1.00287e-05
   -1.97957e-05
  119     5.79720e-06     1.16164e-05     3.36557e-11
    2.81497e-06
    4.38969e-06
  120     3.19460e-10    -1.72199e-10     6.58936e-17
    3.25087e-07
   -1.62224e-07
  121    -5.10669e-12    -1.02340e-11     2.61210e-23
   -1.95333e-12
   -4.13003e-12
  122    -1.11022e-16     0.00000e+00     4.94271e-30
   -8.90399e-14
    4.44089e-14
  123     0.00000e+00     0.00000e+00     0.00000e+00
    0.00000e+00
    0.00000e+00

Obviously it converges very nicely at the end with the maximum gradient
component rapidly going to exactly zero.

But why does it take so long (119 iterations) to get to the final rapid
convergence phase (both for the GSL C version and my fortran implementation
of the same algorithm)?  In Fletcher's book, "Practical Methods of
Optimization", 2nd edition, the problem only requires 21 iterations to
converge (see Table 3.5.1) with Fletcher's own BFGS code.  Although, my
experients showed the L-BFGS-B implementation is fundamentally non-robust,
for this test function it seems to be okay, and it converges in 39
iterations confirming there is something wrong with the efficiency of the
preliminary convergence steps in the GSL case.

I suspect there is a substantial inefficiency problem with the GSL line-
search implementation (since the BFGS update part is completely 
straightforward).  I am not any sort of expert in this field though so I
cannot pin down where the problem is occurring, but a superficial google
search reveals a number of complaints about the GSL minimization efficiency.

Thus, I urge someone with some experience of line-search algorithms should
have a look to make sure some easily rectifiable mistake is not being made.
A simple bug fix that results in a factor of 3 to 5 in efficiency
improvement in optimization would be well worth having.

Meanwhile, I am about to embark on a long series of tests of my fortran
implementation of the GSL BFGS technique for solving the equation of state
for a wide variety of physical conditions, and I will report back here how
robust it is.  Hopefully, it will be completely robust (always determine a
good solution for all physical conditions) for my equation of state problem.

Alan
__________________________
Alan W. Irwin

Astronomical research affiliation with Department of Physics and Astronomy,
University of Victoria (astrowww.phys.uvic.ca).

Programming affiliations with the FreeEOS equation-of-state implementation
for stellar interiors (freeeos.sf.net); PLplot scientific plotting software
package (plplot.org); the Yorick front-end to PLplot (yplot.sf.net); the
Loads of Linux Links project (loll.sf.net); and the Linux Brochure Project
(lbproject.sf.net).
__________________________

Linux-powered Science
__________________________