lbfgs.c

Hieu Xuan Phan & Minh Le Nguyen 利用CRF统计模型写的可用于英文命名实体识别、英文分词的工具（开放源码）。CRF模型最早由Lafferty提出

/*
  MeCab -- Yet Another Part-of-Speech and Morphological Analyzer
 
  $Id: lbfgs.c,v 1.2 2005/05/29 07:27:10 taku-ku Exp $;

  Copyright (C) 2001-2004 Taku Kudo <[e-mail]>
  This is free software with ABSOLUTELY NO WARRANTY.
  
  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/  

#include <math.h>

typedef int integer;
typedef unsigned int uinteger;
typedef char *address;
typedef short int shortint;
typedef float real;
typedef double doublereal;
typedef int logical;
typedef short int shortlogical;
typedef char logical1;
typedef char integer1;

#define TRUE_ (1)
#define FALSE_ (0)
#define abs(x) ((x) >= 0 ? (x) : -(x))
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))

/* Common Block Declarations */
struct lb3_1_ {
  integer mp, lp;
  doublereal gtol, stpmin, stpmax;
};

/* Table of constant values */
static struct lb3_1_ lb3_1;
static integer c__1 = 1;

static doublereal ddot_  ();
static int daxpy_ ();
static int mcsrch_();
static int mcstep_();

/*     ---------------------------------------------------------------------- */
/*     This file contains the LBFGS algorithm and supporting routines */

/*     **************** */
/*     LBFGS SUBROUTINE */
/*     **************** */

/* Subroutine */ int lbfgs(n, m, x, f, g, diagco, diag, iprint, eps, xtol, w, iflag)
     integer *n, *m;
     doublereal *x, *f, *g;
     integer *diagco;
     doublereal *diag;
     integer *iprint;
     doublereal *eps, *xtol, *w;
     integer *iflag;
{
  /* Initialized data */
  static doublereal one = 1.0;
  static doublereal zero = 0.0;
   
  /* System generated locals */
  integer i__1;
  doublereal d__1;

  /* Builtin functions */
  double sqrt();

  /* Local variables */
  static doublereal beta;
  static integer inmc;
  static integer info, iscn, nfev, iycn, iter;
  static doublereal ftol;
  static integer nfun, ispt, iypt, i__, bound;
  static doublereal gnorm;

  static integer point;
  static doublereal xnorm;
  static integer cp;
  static doublereal sq, yr, ys;
  static logical finish;
  static doublereal yy;
  static integer maxfev;

  static integer npt;
  static doublereal stp, stp1;
   
  lb3_1.mp = 6;
  lb3_1.lp = 6;
  lb3_1.gtol = .9;
  lb3_1.stpmin = 1e-20;
  lb3_1.stpmax = 1e20;

  /* Parameter adjustments */
  --diag;
  --g;
  --x;
  --w;
  --iprint;

  /* Function Body */

  /*     INITIALIZE */
  /*     ---------- */

  if (*iflag == 0) {
    goto L10;
  }
  switch ((int)*iflag) {
  case 1:  goto L172;
  case 2:  goto L100;
  }
 L10:
  iter = 0;
  if (*n <= 0 || *m <= 0) {
    goto L196;
  }
  if (lb3_1.gtol <= 1e-4) {
    if (lb3_1.lp > 0) {}
    lb3_1.gtol = .9;
  }
  nfun = 1;
  point = 0;
  finish = FALSE_;
  if (*diagco != 0) {
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
      /* L30: */
      if (diag[i__] <= zero) {
	goto L195;
      }
    }
  } else {
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
      /* L40: */
      diag[i__] = 1.;
    }
  }

  ispt = *n + (*m << 1);
  iypt = ispt + *n * *m;
  i__1 = *n;
  for (i__ = 1; i__ <= i__1; ++i__) {
    /* L50: */
    w[ispt + i__] = -g[i__] * diag[i__];
  }
  gnorm = sqrt(ddot_(n, &g[1], &c__1, &g[1], &c__1));
  stp1 = one / gnorm;

  /*     PARAMETERS FOR LINE SEARCH ROUTINE */

  ftol = 1e-4;
  maxfev = 20;

  /*    if (iprint[1] >= 0) {
	lb1_(&iprint[1], &iter, &nfun, &gnorm, n, m, &x[1], f, &g[1], &stp, &
	finish);
	} */

  /*    -------------------- */
  /*     MAIN ITERATION LOOP */
  /*    -------------------- */

 L80:
  ++iter;
  info = 0;
  bound = iter - 1;
  if (iter == 1) {
    goto L165;
  }
  if (iter > *m) {
    bound = *m;
  }

  ys = ddot_(n, &w[iypt + npt + 1], &c__1, &w[ispt + npt + 1], &c__1);
  if (*diagco == 0) {
    yy = ddot_(n, &w[iypt + npt + 1], &c__1, &w[iypt + npt + 1], &c__1);
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
      /* L90: */
      diag[i__] = ys / yy;
    }
  } else {
    *iflag = 2;
    return 0;
  }
 L100:
  if (*diagco != 0) {
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
      /* L110: */
      if (diag[i__] <= zero) {
	goto L195;
      }
    }
  }

  /*     COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, */
  /*     "Updating quasi-Newton matrices with limited storage", */
  /*     Mathematics of Computation, Vol.24, No.151, pp. 773-782. */
  /*     --------------------------------------------------------- */

  cp = point;
  if (point == 0) {
    cp = *m;
  }
  w[*n + cp] = one / ys;
  i__1 = *n;
  for (i__ = 1; i__ <= i__1; ++i__) {
    /* L112: */
    w[i__] = -g[i__];
  }
  cp = point;
  i__1 = bound;
  for (i__ = 1; i__ <= i__1; ++i__) {
    --cp;
    if (cp == -1) {
      cp = *m - 1;
    }
    sq = ddot_(n, &w[ispt + cp * *n + 1], &c__1, &w[1], &c__1);
    inmc = *n + *m + cp + 1;
    iycn = iypt + cp * *n;
    w[inmc] = w[*n + cp + 1] * sq;
    d__1 = -w[inmc];
    daxpy_(n, &d__1, &w[iycn + 1], &c__1, &w[1], &c__1);
    /* L125: */
  }

  i__1 = *n;
  for (i__ = 1; i__ <= i__1; ++i__) {
    /* L130: */
    w[i__] = diag[i__] * w[i__];
  }

  i__1 = bound;
  for (i__ = 1; i__ <= i__1; ++i__) {
    yr = ddot_(n, &w[iypt + cp * *n + 1], &c__1, &w[1], &c__1);
    beta = w[*n + cp + 1] * yr;
    inmc = *n + *m + cp + 1;
    beta = w[inmc] - beta;
    iscn = ispt + cp * *n;
    daxpy_(n, &beta, &w[iscn + 1], &c__1, &w[1], &c__1);
    ++cp;
    if (cp == *m) {
      cp = 0;
    }
    /* L145: */
  }

  /*     STORE THE NEW SEARCH DIRECTION */
  /*     ------------------------------ */

  i__1 = *n;
  for (i__ = 1; i__ <= i__1; ++i__) {
    /* L160: */
    w[ispt + point * *n + i__] = w[i__];
  }

  /*     OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION */
  /*     BY USING THE LINE SEARCH ROUTINE MCSRCH */
  /*     ---------------------------------------------------- */
 L165:
  nfev = 0;
  stp = one;
  if (iter == 1) {
    stp = stp1;
  }
  i__1 = *n;
  for (i__ = 1; i__ <= i__1; ++i__) {
    /* L170: */
    w[i__] = g[i__];
  }
 L172:
  mcsrch_(n, &x[1], f, &g[1], &w[ispt + point * *n + 1], &stp, &ftol, xtol, 
	  &maxfev, &info, &nfev, &diag[1]);
  if (info == -1) {
    *iflag = 1;
    return 0;
  }
  if (info != 1) {
    goto L190;
  }
  nfun += nfev;

  /*     COMPUTE THE NEW STEP AND GRADIENT CHANGE */
  /*     ----------------------------------------- */

  npt = point * *n;
  i__1 = *n;
  for (i__ = 1; i__ <= i__1; ++i__) {
    w[ispt + npt + i__] = stp * w[ispt + npt + i__];
    /* L175: */
    w[iypt + npt + i__] = g[i__] - w[i__];
  }
  ++point;
  if (point == *m) {
    point = 0;
  }

  /*     TERMINATION TEST */
  /*     ---------------- */

  gnorm = sqrt(ddot_(n, &g[1], &c__1, &g[1], &c__1));
  xnorm = sqrt(ddot_(n, &x[1], &c__1, &x[1], &c__1));
  xnorm = max(1.,xnorm);
  if (gnorm / xnorm <= *eps) {
    finish = TRUE_;
  }

  /*    if (iprint[1] >= 0) {
	lb1_(&iprint[1], &iter, &nfun, &gnorm, n, m, &x[1], f, &g[1], &stp, &
	finish);
	}*/
  if (finish) {
    *iflag = 0;
    return 0;
  }
  goto L80;

  /*     ------------------------------------------------------------ */
  /*     END OF MAIN ITERATION LOOP. ERROR EXITS. */
  /*     ------------------------------------------------------------ */

 L190:
  *iflag = -1;
  return 0;
 L195:
  *iflag = -2;
  return 0;
 L196:
  *iflag = -3;

  return 0;
} /* lbfgs_ */

/*   ---------------------------------------------------------- */

/* Subroutine */ static int daxpy_(n, da, dx, incx, dy, incy)
     integer *n;
     doublereal *da, *dx;
     integer *incx;
     doublereal *dy;
     integer *incy;
{
  /* System generated locals */
  integer i__1;

  /* Local variables */
  static integer i__, m, ix, iy, mp1;


  /*     constant times a vector plus a vector. */
  /*     uses unrolled loops for increments equal to one. */
  /*     jack dongarra, linpack, 3/11/78. */


  /* Parameter adjustments */
  --dy;
  --dx;

  /* Function Body */
  if (*n <= 0) {
    return 0;
  }
  if (*da == 0.) {
    return 0;
  }
  if (*incx == 1 && *incy == 1) {
    goto L20;
  }

  /*        code for unequal increments or equal increments */
  /*          not equal to 1 */

  ix = 1;
  iy = 1;
  if (*incx < 0) {
    ix = (-(*n) + 1) * *incx + 1;
  }
  if (*incy < 0) {
    iy = (-(*n) + 1) * *incy + 1;
  }
  i__1 = *n;
  for (i__ = 1; i__ <= i__1; ++i__) {
    dy[iy] += *da * dx[ix];
    ix += *incx;
    iy += *incy;
    /* L10: */
  }
  return 0;

  /*        code for both increments equal to 1 */


  /*        clean-up loop */

 L20:
  m = *n % 4;
  if (m == 0) {
    goto L40;
  }
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__) {
    dy[i__] += *da * dx[i__];
    /* L30: */
  }
  if (*n < 4) {
    return 0;
  }
 L40:
  mp1 = m + 1;
  i__1 = *n;
  for (i__ = mp1; i__ <= i__1; i__ += 4) {
    dy[i__] += *da * dx[i__];
    dy[i__ + 1] += *da * dx[i__ + 1];
    dy[i__ + 2] += *da * dx[i__ + 2];
    dy[i__ + 3] += *da * dx[i__ + 3];
    /* L50: */
  }
  return 0;
} /* daxpy_ */



/*   ---------------------------------------------------------- */

static doublereal ddot_(n, dx, incx, dy, incy)
     integer *n;
     doublereal *dx;
     integer *incx;
     doublereal *dy;
     integer *incy;
{
  /* System generated locals */
  integer i__1;
  doublereal ret_val;

  /* Local variables */
  static integer i__, m;
  static doublereal dtemp;
  static integer ix, iy, mp1;


  /*     forms the dot product of two vectors. */
  /*     uses unrolled loops for increments equal to one. */
  /*     jack dongarra, linpack, 3/11/78. */


  /* Parameter adjustments */
  --dy;
  --dx;

  /* Function Body */
  ret_val = 0.;
  dtemp = 0.;
  if (*n <= 0) {
    return ret_val;
  }
  if (*incx == 1 && *incy == 1) {
    goto L20;
  }

  /*        code for unequal increments or equal increments */
  /*          not equal to 1 */

  ix = 1;
  iy = 1;
  if (*incx < 0) {
    ix = (-(*n) + 1) * *incx + 1;
  }
  if (*incy < 0) {
    iy = (-(*n) + 1) * *incy + 1;
  }
  i__1 = *n;
  for (i__ = 1; i__ <= i__1; ++i__) {
    dtemp += dx[ix] * dy[iy];
    ix += *incx;