lbfgs.c

来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,199 行 · 第 1/3 页

#include "f2c.h"
#include "netlib.h"
extern double sqrt(double); /* #include <math.h> */

static void mcsrch_(integer *n, doublereal *x, doublereal *f, doublereal *g, doublereal *s, doublereal *stp,
                    doublereal *ftol, doublereal *xtol, integer *maxfev, integer *info, integer *nfev, doublereal *wa);
static void mcstep_(doublereal *stx, doublereal *fx, doublereal *dx, doublereal *sty, doublereal *fy, doublereal *dy,
                    doublereal *stp, doublereal *fp, doublereal *dp, logical *brackt,
                    doublereal *stpmin, doublereal *stpmax, integer *info);
void lb1_(int *iprint, int *iter, int *nfun, double *gnorm, int *n, int *m,
          double *x, double *f, double *g, double *stp, int *finish);
void lbptf_(char* msg);
void lbp1d_(char* msg, int* i);

/* Common Block Declarations */

struct lb3_1_ {
    integer mp, lp;
    doublereal gtol, stpmin, stpmax, stpawf;
};

/*#define lb3_1 (*(struct lb3_1_ *) &lb3_)*/
#define lb3_1 lb3_

/* Initialized data */
struct lb3_1_ lb3_1 = { 6, 6, .9, 1e-20, 1e20, 1. };

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

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

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

/* Subroutine */ void lbfgs_(n, m, x, f, g, diagco, diag, iprint, eps, xtol, w, iflag)
integer *n, *m;
doublereal *x, *f, *g;
logical *diagco;
doublereal *diag;
integer *iprint;
doublereal *eps, *xtol, *w;
integer *iflag;
{
    /* System generated locals */
    doublereal d__1;

    /* Local variables */
    static doublereal beta;
    static integer inmc;
    static integer info, iscn, nfev, iycn, iter;
    static doublereal ftol;
    static integer nfun, ispt, iypt;
    static integer 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;

/*        LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION */
/*                          JORGE NOCEDAL */
/*                        *** July 1990 *** */

/*     This subroutine solves the unconstrained minimization problem */

/*                      min F(x),    x= (x1,x2,...,xN), */

/*      using the limited memory BFGS method. The routine is especially */
/*      effective on problems involving a large number of variables. In */
/*      a typical iteration of this method an approximation Hk to the */
/*      inverse of the Hessian is obtained by applying M BFGS updates to */
/*      a diagonal matrix Hk0, using information from the previous M steps.  */
/*      The user specifies the number M, which determines the amount of */
/*      storage required by the routine. The user may also provide the */
/*      diagonal matrices Hk0 if not satisfied with the default choice. */
/*      The algorithm is described in "On the limited memory BFGS method */
/*      for large scale optimization", by D. Liu and J. Nocedal, */
/*      Mathematical Programming B 45 (1989) 503-528. */

/*      The user is required to calculate the function value F and its */
/*      gradient G. In order to allow the user complete control over */
/*      these computations, reverse  communication is used. The routine */
/*      must be called repeatedly under the control of the parameter */
/*      IFLAG. */

/*      The steplength is determined at each iteration by means of the */
/*      line search routine MCVSRCH, which is a slight modification of */
/*      the routine CSRCH written by More' and Thuente. */

/*      The calling statement is */

/*          CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG) */

/*      where */

/*     N       is an INTEGER variable that must be set by the user to the */
/*             number of variables. It is not altered by the routine. */
/*             Restriction: N>0. */

/*     M       is an INTEGER variable that must be set by the user to */
/*             the number of corrections used in the BFGS update. It */
/*             is not altered by the routine. Values of M less than 3 are */
/*             not recommended; large values of M will result in excessive */
/*             computing time. 3<= M <=7 is recommended. Restriction: M>0.  */

/*     X       is a DOUBLE PRECISION array of length N. On initial entry */
/*             it must be set by the user to the values of the initial */
/*             estimate of the solution vector. On exit with IFLAG=0, it */
/*             contains the values of the variables at the best point */
/*             found (usually a solution). */

/*     F       is a DOUBLE PRECISION variable. Before initial entry and on */
/*             a re-entry with IFLAG=1, it must be set by the user to */
/*             contain the value of the function F at the point X. */

/*     G       is a DOUBLE PRECISION array of length N. Before initial */
/*             entry and on a re-entry with IFLAG=1, it must be set by */
/*             the user to contain the components of the gradient G at */
/*             the point X. */

/*     DIAGCO  is a LOGICAL variable that must be set to .TRUE. if the */
/*             user  wishes to provide the diagonal matrix Hk0 at each */
/*             iteration. Otherwise it should be set to .FALSE., in which */
/*             case  LBFGS will use a default value described below. If */
/*             DIAGCO is set to .TRUE. the routine will return at each */
/*             iteration of the algorithm with IFLAG=2, and the diagonal */
/*              matrix Hk0  must be provided in the array DIAG. */

/*     DIAG    is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE., */
/*             then on initial entry or on re-entry with IFLAG=2, DIAG */
/*             it must be set by the user to contain the values of the */
/*             diagonal matrix Hk0.  Restriction: all elements of DIAG */
/*             must be positive. */

/*     IPRINT  is an INTEGER array of length two which must be set by the */
/*             user. */

/*             IPRINT(1) specifies the frequency of the output: */
/*                IPRINT(1) < 0 : no output is generated, */
/*                IPRINT(1) = 0 : output only at first and last iteration, */
/*                IPRINT(1) > 0 : output every IPRINT(1) iterations. */

/*             IPRINT(2) specifies the type of output generated: */
/*                IPRINT(2) = 0 : iteration count, number of function */
/*                                evaluations, function value, norm of the */
/*                                gradient, and steplength, */
/*                IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of */
/*                                variables and  gradient vector at the */
/*                                initial point, */
/*                IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of */
/*                                variables, */
/*                IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector.*/


/*    EPS     is a positive DOUBLE PRECISION variable that must be set by */
/*            the user, and determines the accuracy with which the solution*/
/*            is to be found. The subroutine terminates when */

/*                         ||G|| < EPS max(1,||X||), */

/*            where ||.|| denotes the Euclidean norm. */

/*    XTOL    is a  positive DOUBLE PRECISION variable that must be set by */
/*            the user to an estimate of the machine precision (e.g. */
/*            10**(-16) on a SUN station 3/60). The line search routine will*/
/*            terminate if the relative width of the interval of uncertainty*/
/*            is less than XTOL. */

/*     W       is a DOUBLE PRECISION array of length N(2M+1)+2M used as */
/*             workspace for LBFGS. This array must not be altered by the */
/*             user. */

/*    IFLAG   is an INTEGER variable that must be set to 0 on initial entry*/
/*            to the subroutine. A return with IFLAG<0 indicates an error, */
/*            and IFLAG=0 indicates that the routine has terminated without*/
/*            detecting errors. On a return with IFLAG=1, the user must */
/*            evaluate the function F and gradient G. On a return with */
/*            IFLAG=2, the user must provide the diagonal matrix Hk0. */

/*            The following negative values of IFLAG, detecting an error, */
/*            are possible: */

/*              IFLAG=-1  The line search routine MCSRCH failed. The */
/*                        parameter INFO provides more detailed information */
/*                        (see also the documentation of MCSRCH): */

/*                       INFO = 0  IMPROPER INPUT PARAMETERS. */

/*                       INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF */
/*                                 UNCERTAINTY IS AT MOST XTOL. */

/*                       INFO = 3  MORE THAN 20 FUNCTION EVALUATIONS WERE */
/*                                 REQUIRED AT THE PRESENT ITERATION. */

/*                       INFO = 4  THE STEP IS TOO SMALL. */

/*                       INFO = 5  THE STEP IS TOO LARGE. */

/*                       INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.*/
/*                                 THERE MAY NOT BE A STEP WHICH SATISFIES */
/*                                 THE SUFFICIENT DECREASE AND CURVATURE */
/*                                 CONDITIONS. TOLERANCES MAY BE TOO SMALL.  */

/*             IFLAG=-2  The i-th diagonal element of the diagonal inverse */
/*                       Hessian approximation, given in DIAG, is not */
/*                       positive. */

/*             IFLAG=-3  Improper input parameters for LBFGS (N or M are */
/*                       not positive). */


/*    ON THE DRIVER: */

/*    The program that calls LBFGS must contain the declaration: */

/*                       EXTERNAL LB2 */

/*    LB2 is a BLOCK DATA that defines the default values of several */
/*    parameters described in the COMMON section. */

/*    COMMON: */

/*     The subroutine contains one common area, which the user may wish to */
/*    reference: */

/* awf added stpawf */

/*    MP  is an INTEGER variable with default value 6. It is used as the */
/*        unit number for the printing of the monitoring information */
/*        controlled by IPRINT. */

/*    LP  is an INTEGER variable with default value 6. It is used as the */
/*        unit number for the printing of error messages. This printing */
/*        may be suppressed by setting LP to a non-positive value. */

/*    GTOL is a DOUBLE PRECISION variable with default value 0.9, which */
/*        controls the accuracy of the line search routine MCSRCH. If the */
/*        function and gradient evaluations are inexpensive with respect */
/*        to the cost of the iteration (which is sometimes the case when */
/*        solving very large problems) it may be advantageous to set GTOL */
/*        to a small value. A typical small value is 0.1.  Restriction: */
/*        GTOL should be greater than 1.D-04. */

/*    STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which */
/*        specify lower and uper bounds for the step in the line search.  */
/*        Their default values are 1.D-20 and 1.D+20, respectively. These */
/*        values need not be modified unless the exponents are too large */
/*        for the machine being used, or unless the problem is extremely */
/*        badly scaled (in which case the exponents should be increased).  */


/*  MACHINE DEPENDENCIES */

/*        The only variables that are machine-dependent are XTOL, */
/*        STPMIN and STPMAX. */


/*  GENERAL INFORMATION */

/*    Other routines called directly:  DAXPY, DDOT, LB1, MCSRCH */

/*    Input/Output  :  No input; diagnostic messages on unit MP and */
/*                     error messages on unit LP. */


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

/*   ---------------------------------------------------------- */
/*     DATA */
/*   ---------------------------------------------------------- */

/*      BLOCK DATA LB3 */

/*     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) {
            lbptf_("  GTOL IS LESS THAN OR EQUAL TO 1.D-04");
            lbptf_("  IT HAS BEEN RESET TO 9.D-01");
        }
        lb3_1.gtol = .9;
    }
    nfun = 1;
    point = 0;
    finish = FALSE_;
    if (*diagco) {
        for (i = 0; i < *n; ++i) {
            if (diag[i] <= 0.) {
                goto L195;
            }
        }
    } else {
        for (i = 0; i < *n; ++i) {
            diag[i] = 1.;
        }
    }

/*     THE WORK VECTOR W IS DIVIDED AS FOLLOWS: */
/*     --------------------------------------- */
/*     THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND */
/*         OTHER TEMPORARY INFORMATION. */
/*     LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO. */
/*     LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED */
/*         IN THE FORMULA THAT COMPUTES H*G. */
/*     LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH */
/*         STEPS. */
/*     LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M */
/*         GRADIENT DIFFERENCES. */

/*     THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A */
/*     CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT. */

    ispt = *n + (*m << 1);
    iypt = ispt + *n * *m;
    for (i = 0; i < *n; ++i) {
        w[ispt + i] = -g[i] * diag[i];
    }
    gnorm = sqrt(ddot_(n, g, &c__1, g, &c__1));
    stp1 = 1. / gnorm;

/*     PARAMETERS FOR LINE SEARCH ROUTINE */

    ftol = 1e-4;
    maxfev = 20;

    if (iprint[0] >= 0) {
        lb1_(iprint, &iter, &nfun, &gnorm, n, m, x, f, g, &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], &c__1, &w[ispt + npt], &c__1);
    if (! (*diagco)) {
        yy = ddot_(n, &w[iypt + npt], &c__1, &w[iypt + npt], &c__1);
        for (i = 0; i < *n; ++i) {
            diag[i] = ys / yy;
        }
    } else {
        *iflag = 2;
        return;
    }
L100:
    if (*diagco) {
        for (i = 0; i < *n; ++i) {
            if (diag[i] <= 0.) {
                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-1] = 1. / ys;
    for (i = 0; i < *n; ++i) {
        w[i] = -g[i];
    }
    cp = point;
    for (i = 0; i < bound; ++i) {
        --cp;