chokkan / liblbfgs Public

The default line search for `liblbfgs` is specified by [1]. Its implementation had a bug
when searching for the minimum of a cubic interpolant used to identify the next trial
point during the line search. This patch fixes the bug.

In #28, user @jonasson2 identified a bug
in the `CUBIC_MINIMIZER2` macro, which finds the minimum of the cubic interpolant
used for trial step size selection. To be honest, @jonasson2 did most of the hard work
here, even providing a test case, I'm just writing up my understanding and submitting
the PR so the bug is fixed.

`CUBIC_MINIMIZER2` minimizes the cubic fit in case 3 of [1].

We can see similar logic in the `MCSTEP` routine of Nocedal's original [code](http://users.iems.northwestern.edu/~nocedal/lbfgs.html).
```
C
C        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
C        TO INFINITY IN THE DIRECTION OF THE STEP.
C
         GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S)))
         IF (STP .GT. STX) GAMMA = -GAMMA
         P = (GAMMA - DP) + THETA
         Q = (GAMMA + (DX - DP)) + GAMMA
         R = P/Q
         IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN
            STPC = STP + R*(STX - STP)
         ELSE IF (STP .GT. STX) THEN
            STPC = STPMAX
         ELSE
            STPC = STPMIN
            END IF
```

Above, the Fortran condition `STP .GT. STX` is equivalent to `d > 0` in C as `STP, STX` correspond to `v, u`.

[1] Jorge J. More and David J. Thuente. Line search algorithm with
guaranteed sufficient decrease. ACM Transactions on Mathematical
Software (TOMS), Vol 20, No 3, pp. 286-307, 1994.

At the previous commit 7fc7876 it is possible to trigger faulty behavior by simply
adding a `params.gtol = 0.1` to the sample Rosenbrock optimization function.

The patch at the end of this message below does so.

On my machine, this results in a failed optimization, ending as follows:

```
Iteration 7:
  fx = 95.209157, x[0] = -0.287044, x[1] = 0.032625
  xnorm = 2.042774, gnorm = 91.591482, step = 1.355833

L-BFGS optimization terminated with status code = -1001
  fx = 69.604508, x[0] = -0.287044, x[1] = 0.032625
```

After my change, optimization terminates successfully, showing

```
Iteration 23:
  fx = 0.000000, x[0] = 1.000000, x[1] = 1.000000
  xnorm = 10.000003, gnorm = 0.000078, step = 1.000000

L-BFGS optimization terminated with status code = 0
  fx = 0.000000, x[0] = 1.000000, x[1] = 1.000000
```

Testing out a couple of other settings, including default, results in normal optimization results.

```
diff --git a/sample/sample.c b/sample/sample.c
index 90e36cd..512bc7b 100644
--- a/sample/sample.c
+++ b/sample/sample.c
@@ -70,6 +70,7 @@ int main(int argc, char *argv[])
         Start the L-BFGS optimization; this will invoke the callback functions
         evaluate() and progress() when necessary.
      */
+    param.gtol = 0.1;
     ret = lbfgs(N, x, &fx, evaluate, progress, NULL, &param);

     /* Report the result. */
```