lbfgs.c

NIST Handwriting OCR Testbed

  width = stpmax - stpmin;
  width1 = width / .5;
  for(j = 1; j <= n; j++)
    wa[j] = x[j];

  /* The variables stx, fx, dgx contain the values of the step,
  function, and directional derivative at the best step.  The
  variables sty, fy, dgy contain the value of the step, function, and
  derivative at the other endpoint of the interval of uncertainty.
  The variables stp, f, dg contain the values of the step, function,
  and derivative at the current step. */
  stx = 0.;
  fx = finit;
  dgx = dginit;
  sty = 0.;
  fy = finit;
  dgy = dginit;

  /* Start of iteration. */

  /* [Main loop:] */
  while(1) {

    /* Set the minimum and maximum steps to correspond to the present
    interval of uncertainty. */
    if(brackt) {
      stmin = min(stx, sty);
      stmax = max(stx, sty);
    }
    else {
      stmin = stx;
      stmax = *stp + XTRAPF * (*stp - stx);
    }

    /* Force the step to be within the bounds stpmax and stpmin. */
    *stp = max(*stp, stpmin);
    *stp = min(*stp, stpmax);

    /* If an unusual termination is to occur then let stp be the
    lowest point obtained so far. */
    if(brackt && (*stp <= stmin || *stp >= stmax) || *nfev >=
      maxfev - 1 || infoc == 0 || brackt && stmax - stmin <= xtol *
      stmax)
      *stp = stx;

    /* Evaluate the function and gradient at stp and compute the
    directional derivative.  We return to main program to obtain f and
    g. */
    for(j = 1; j <= n; j++)
      x[j] = wa[j] + *stp * s[j];
    *info = -1;
    return;

  L45:
    *info = 0;
    (*nfev)++;
    dg = 0.;
    for(j = 1; j <= n; j++)
      dg += g[j] * s[j];
    ftest1 = finit + *stp * dgtest;

    /* Test for convergence. */
    if(brackt && (*stp <= stmin || *stp >= stmax) || infoc == 0)
      *info = 6;
    if(*stp == stpmax && f <= ftest1 && dg <= dgtest)
      *info = 5;
    if(*stp == stpmin && (f > ftest1 || dg >= dgtest))
      *info = 4;
    if(*nfev >= maxfev)
      *info = 3;
    if(brackt && stmax - stmin <= xtol * stmax)
      *info = 2;
    if(f <= ftest1 && fabs((double)dg) <= gtol * (-dginit))
      *info = 1;

    /* Check for termination. */
    if(*info != 0)
      return;

    /* In the first stage we seek a step for which the modified
    function has a nonpositive value and nonnegative derivative. */
    if(stage1 && f <= ftest1 && dg >= min(ftol, gtol) * dginit)
      stage1 = FALSE;

    /* A modified function is used to predict the step only if we have
    not obtained a step for which the modified function has a
    nonpositive function value and nonnegative derivative, and if a
    lower function value has been obtained but the decrease is not
    sufficient. */
    if(stage1 && f <= fx && f > ftest1) {
      /* Define the modified function and derivative values. */
      fm = f - *stp * dgtest;
      fxm = fx - stx * dgtest;
      fym = fy - sty * dgtest;
      dgm = dg - dgtest;
      dgxm = dgx - dgtest;
      dgym = dgy - dgtest;

      /* Call cstep to update the interval of uncertainty and to
      compute the new step. */
      mcstep(&stx, &fxm, &dgxm, &sty, &fym, &dgym, stp, &fm, &dgm,
        &brackt, stmin, stmax, &infoc);

      /* Reset the function and gradient values for f. */
      fx = fxm + stx * dgtest;
      fy = fym + sty * dgtest;
      dgx = dgxm + dgtest;
      dgy = dgym + dgtest;
    }
    else {
      /* Call mcstep to update the interval of uncertainty and to
      compute the new step. */
      mcstep(&stx, &fx, &dgx, &sty, &fy, &dgy, stp, &f, &dg, &brackt,
        stmin, stmax, &infoc);
    }

    /* Force a sufficient decrease in the size of the interval of
    uncertainty. */
    if(brackt) {
/*      if((r1 = sty - stx, fabs((double)r1)) >= .66 * width1) */
      if(fabs((double)(sty - stx)) >= .66 * width1)
	*stp = stx + .5 * (sty - stx);
      width1 = width;
/*      width = (r1 = sty - stx, fabs((double)r1)); */
      width = fabs((double)(sty - stx));
    }

    /* End of iteration. */
  } /* while(1) [main loop] */

}

/* end of function mcsrch */

/*******************************************************************/

  /* ---------------
     function mcstep
     ---------------

  The purpose of mcstep is to compute a safeguarded step for a
  linesearch and to update an interval of uncertainty for a minimizer
  of the function.

  The parameter stx contains the step with the least function value.
  The parameter stp contains the current step.  It is assumed that the
  derivative at stx is negative in the direction of the step.  If
  brackt is set to TRUE then a minimizer has been bracketed in an
  interval of uncertainty with endpoints stx and sty.

  The calling statement is:

    mcstep(stx, fx, dx, sty, fy, dy, stp, fp, dp, brackt, stpmin,
      stpmax, info);

  where:

    stx, fx, and dx are variables which specify the step, the
      function, and the derivative at the best step obtained so far.
      The derivative must be negative in the direction of the step,
      that is, dx and stp-stx must have opposite signs.  On output
      these parameters are updated appropriately.

    sty, fy, and dy are variables which specify the step, the
      function, and the derivative at the other endpoint of the
      interval of uncertainty.  On output these parameters are
      updated appropriately.

    stp, fp, and dp are variables which specify the step, the
      function, and the derivative at the current step.  If brackt is
      set to TRUE then on input stp must be between stx and sty.  On
      output stp is set to the new step.

    brackt is an int variable which specifies if a minimizer has been
      bracketed.  If the minimizer has not been bracketed then on
      input brackt must be set to FALSE.  If the minimizer is
      bracketed then on output brackt is set to TRUE.

    stpmin and stpmax are input variables which specify lower and
      upper bounds for the step.

    info is an int output variable set as follows: if info =
      1, 2, 3, 4, or 5, then the step has been computed according to
      one of the five cases below;  otherwise info = 0, and this
      indicates improper input parameters.


  Argonne National Laboratory.  MINPACK Project.  June 1983.
  Jorge J. More', David J. Thuente. */

void
mcstep(stx, fx, dx, sty, fy, dy, stp, fp, dp, brackt, stpmin, stpmax,
  info)
float *stx, *fx, *dx, *sty, *fy, *dy, *stp, *dp, stpmin, stpmax;
double *fp;
int *brackt, *info;
{
  static int bound;
  float r1, r2, r3;
  static float sgnd, stpc, stpf, stpq, p, q, gamma, r, s, theta;

  *info = 0;

  /* Check the input parameters for errors. */
  if(*brackt && (*stp <= min(*stx, *sty) || *stp >= max(*stx, *sty))
    || * dx * (*stp - *stx) >= 0. || stpmax < stpmin)
    return;

  /* Determine if the derivatives have opposite sign. */
  sgnd = *dp * (*dx / fabs((double)(*dx)));

  if(*fp > *fx) {

    /* First case.  A higher function value.  The minimum is
    bracketed.  If the cubic step is closer to stx than the quadratic
    step, the cubic step is taken, else the average of the cubic and
    quadratic steps is taken. */
    *info = 1;
    bound = TRUE;
    theta = (*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
    /* computing max */
    r1 = fabs((double)theta), r2 = fabs((double)(*dx)),
      r1 = max(r1, r2), r2 = fabs((double)(*dp));
    s = max(r1, r2);
    /* computing 2nd power */
    r1 = theta / s;
    gamma = s * sqrt((double)(r1 * r1 - *dx / s * (*dp / s)));
    if(*stp < *stx)
      gamma = -gamma;
    p = gamma - *dx + theta;
    q = gamma - *dx + gamma + *dp;
    r = p / q;
    stpc = *stx + r * (*stp - *stx);
    stpq = *stx + *dx / ((*fx - *fp) / (*stp - *stx) + *dx) / 2 *
      (*stp - *stx);
/*    if((r1 = stpc - *stx, fabs((double)r1)) < (r2 = stpq - *stx,
      fabs((double)r2))) */
    if(fabs((double)(stpc - *stx)) < fabs((double)(stpq - *stx)))
      stpf = stpc;
    else
      stpf = stpc + (stpq - stpc) / 2;
    *brackt = TRUE;
  }

  else if(sgnd < 0.) {

    /* Second case.  A lower function value and derivatives of
    opposite sign.  The minimum is bracketed.  If the cubic step is
    closer to stx than the quadratic (secant) step, the cubic step is
    taken, else the quadratic step is taken. */
    *info = 2;
    bound = FALSE;
    theta = (*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
    /* computing max */
    r1 = fabs((double)theta), r2 = fabs((double)(*dx)),
      r1 = max(r1, r2), r2 = fabs((double)(*dp));
    s = max(r1, r2);
    /* computing 2nd power */
    r1 = theta / s;
    gamma = s * sqrt((double)(r1 * r1 - *dx / s * (*dp / s)));
    if(*stp > *stx)
      gamma = -gamma;
    p = gamma - *dp + theta;
    q = gamma - *dp + gamma + *dx;
    r = p / q;
    stpc = *stp + r * (*stx - *stp);
    stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
/*    if((r1 = stpc - *stp, fabs((double)r1)) > (r2 = stpq - *stp,
      fabs((double)r2))) */
    if(fabs((double)(stpc - *stp)) > fabs((double)(stpq - *stp)))
      stpf = stpc;
    else
      stpf = stpq;
    *brackt = TRUE;
  }

  else if(fabs((double)(*dp)) < fabs((double)(*dx))) {

    /* Third case.  A lower function value, derivatives of the same
    sign, and the magnitude of the derivative decreases.  The cubic
    step is only used if the cubic tends to infinity in the direction
    of the step or if the minimum of the cubic is beyond stp.
    Otherwise the cubic step is defined to be either stpmin or stpmax.
    The quadratic (secant) step is also computed and if the minimum is
    bracketed then the the step closest to stx is taken, else the step
    farthest away is taken. */
    *info = 3;
    bound = TRUE;
    theta = (*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
    /* computing max */
    r1 = fabs((double)theta), r2 = fabs((double)(*dx)),
      r1 = max(r1, r2), r2 = fabs((double)(*dp));
    s = max(r1, r2);
    /* The case gamma = 0 only arises if the cubic does not tend to
    infinity in the direction of the step. */
    /* computing max */
    /* computing 2nd power */
    r3 = theta / s;
    r1 = 0., r2 = r3 * r3 - *dx / s * (*dp / s);
    gamma = s * sqrt((double)max(r1, r2));
    if(*stp > *stx)
      gamma = -gamma;
    p = gamma - *dp + theta;
    q = gamma + (*dx - *dp) + gamma;
    r = p / q;
    if(r < 0. && gamma != 0.)
      stpc = *stp + r * (*stx - *stp);
    else if(*stp > *stx)
      stpc = stpmax;
    else
      stpc = stpmin;
    stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
    if(*brackt) {
/*      if((r1 = *stp - stpc, fabs((double)r1)) < (r2 = *stp - stpq,
        fabs((double)r2))) */
      if(fabs((double)(*stp - stpc)) < fabs((double)(*stp - stpq)))
	stpf = stpc;
      else
	stpf = stpq;
    }
    else {
/*      if((r1 = *stp - stpc, fabs((double)r1)) > (r2 = *stp - stpq,
        fabs((double)r2))) */
      if(fabs((double)(*stp - stpc)) > fabs((double)(*stp - stpq)))
	stpf = stpc;
      else
	stpf = stpq;
    }
  }
  else {

    /* Fourth case.  A lower function value, derivatives of the same
    sign, and the magnitude of the derivative does not decrease.  If
    the minimum is not bracketed, the step is either stpmin or stpmax,
    else the cubic step is taken. */
    *info = 4;
    bound = FALSE;
    if(*brackt) {
      theta = (*fp - *fy) * 3 / (*sty - *stp) + *dy + *dp;
      /* computing max */
      r1 = fabs((double)theta), r2 = fabs((double)(*dy)),
        r1 = max(r1, r2), r2 = fabs((double)(*dp));
      s = max(r1, r2);
      /* computing 2nd power */
      r1 = theta / s;
      gamma = s * sqrt((double)(r1 * r1 - *dy / s * (*dp /s)));
      if(*stp > *sty)
	gamma = -gamma;
      p = gamma - *dp + theta;
      q = gamma - *dp + gamma + *dy;
      r = p / q;
      stpc = *stp + r * (*sty - *stp);
      stpf = stpc;
    }
    else if(*stp > *stx)
      stpf = stpmax;
    else
      stpf = stpmin;
  }

  /* Update the interval of uncertainty.  This update does not depend
  on the new step or the case analysis above. */
  if(*fp > *fx) {
    *sty = *stp;
    *fy = *fp;
    *dy = *dp;
  }
  else {
    if(sgnd < 0.) {
      *sty = *stx;
      *fy = *fx;
      *dy = *dx;
    }
    *stx = *stp;
    *fx = *fp;
    *dx = *dp;
  }

  /* Compute the new step and safeguard it. */
  stpf = min(stpmax, stpf);
  stpf = max(stpmin, stpf);
  *stp = stpf;
  if(*brackt && bound) {
    if(*sty > *stx) {
      /* computing min */
      r1 = *stx + (*sty - *stx) * .66;
      *stp = min(r1, *stp);
    }
    else {
      /* computing max */
      r1 = *stx + (*sty - *stx) * .66;
      *stp = max(r1, *stp);
    }
  }
  return;
}

/* end of function mcstep */

/*******************************************************************/