asa_cg.cpp

数学计算程序

            printf ("Convergence tolerance for gradient satisfied\n") ;
        }
        else if ( status == 1 )
        {
            printf ("Terminating in cg since change in function value "
                    "<= feps*|f|\n") ;
        }
        else if ( status == 2 )
        {
            printf ("Number of iterations exceed specified limits "
                    "for cg routine\n") ;
            printf ("Iterations: %10.0f maxit: %10.0f totit: %10.0f\n",
                    (double) Com.cgiter, (double) Com.cgmaxit,
                    (double) cg_totit) ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, Com.tol) ;
        }
        else if ( status == 3 )
        {
            printf ("Slope always negative in cg line search\n") ;
            printf ("%s\n", mess1) ;
            printf ("   - your cost function has an error\n") ;
            printf ("%s\n", mess4) ;
        }
        else if ( status == 4 )
        {
            printf ("Line search fails in cg, too many secant steps\n") ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, Com.tol) ;
        }
        else if ( status == 5 )
        {
            printf ("Search direction not a descent direction in cg\n") ;
        }
        else if ( status == 6 ) /* line search fails */
        {
            printf ("Line search fails in cg iteration\n") ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, Com.tol) ;
            printf ("%s\n", mess4) ;
            printf ("%s\n", mess5) ;
        }
        else if ( status == 7 ) /* line search fails */
        {
            printf ("Line search fails in cg iteration\n") ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, Com.tol) ;
        }
        else if ( status == 8 ) /* line search fails */
        {
            printf ("Line search fails in cg iteration\n") ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, Com.tol) ;
            printf ("%s\n", mess4) ;
            printf ("%s\n", mess5) ;
        }
        else if ( status == 9 )
        {
            printf ("Debugger is on, function value does not improve in cg\n") ;
            printf ("new value: %25.16e old value: %25.16e\n",
                Com.f_debug, Com.f0) ;
        }
        else if ( status == 10 )
        {
            printf ("Insufficient memory\n") ;
        }
        else if ( status == 11 )
        {
            printf ("Number of iterations or function evaluation exceed "
                          "specified limits for cbb routine\n") ;
            printf ("Iterations: %i maxit: %i totit: %i\n",
                     Com.cbbiter, Com.pgmaxit, cbb_totit) ;
            printf ("Total function evaluations: %i maxfunc: %i\n",
                     Com.nf, Com.pgmaxfunc);
        }
        if ( status == 12 ) /* line search fails in cbb iteration */
        {
            printf ("Line search fails in cbb iteration\n") ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, Com.tol) ;
            printf ("%s\n", mess4) ;
        }

        if ( status == 13 )
        {
            printf ("Search direction not descent direction in "
                    "asa_grad_proj\n") ;
            printf ("directional derivative: %e\n", Com.gtd) ;
        }
        if ( status == 14 )
        {
             printf ("At cbb iteration %i function value became nan\n",
                      Com.cbbiter) ;
        }

        printf ("projected gradient max norm: %13.6e\n", Com.pgnorm) ;
        printf ("function value:              %13.6e\n", Com.f) ;
        printf ("\nTotal cg  iterations:           %10.0f\n",
                (double) Com.cgiter) ;
        printf ("Total cg  function evaluations: %10.0f\n",
                (double) Com.cgfunc) ;
        printf ("Total cg  gradient evaluations: %10.0f\n",
                (double) Com.cggrad) ;
        printf ("Total cbb iterations:           %10.0f\n",
                (double) Com.cbbiter) ;
        printf ("Total cbb function evaluations: %10.0f\n",
                (double) Com.cbbfunc) ;
        printf ("Total cbb gradient evaluations: %10.0f\n",
                    (double) Com.cbbgrad) ;
        printf ("------------------------------------------\n") ;
        printf ("Total function evaluations:     %10.0f\n",
                (double) Com.nf) ;
        printf ("Total gradient evaluations:     %10.0f\n",
                (double) Com.ng) ;
        printf ("==========================================\n\n") ;
    }
    free (ifree) ;
    if ( Work == NULL ) free (work) ;
    if ( Stat != NULL )
    {
        Stat->f = Com.f ;
        Stat->pgnorm = Com.pgnorm ;
        Stat->cgiter = Com.cgiter ;
        Stat->cgfunc = Com.cgfunc ;
        Stat->cggrad = Com.cggrad ;
        Stat->cbbiter = Com.cbbiter ;
        Stat->cbbfunc = Com.cbbfunc ;
        Stat->cbbgrad = Com.cbbgrad ;
    }
    return (status) ;
}

/* =========================================================================
   === asa_default ======================================================
   =========================================================================
   Set default parameter values for the ASA routine. The CG default
   parameter values are set by asa_cg_default.  If the parameter argument of
   asa_descent is NULL, this routine is called by asa_cg automatically.
   If the user wishes to set parameter values, then the asa_parameter structure
   should be allocated in the main program. The user could call asa_default
   to initialize the structure, and then individual elements in the structure
   could be changed, before passing the structure to asa_cg.
   =========================================================================*/
void asa_default
(
    asa_parm *Parm
)
{
    double eps, t ;

    /* T => print final statistics
       F => no printout of statistics */
    Parm->PrintFinal = TRUE ;

    /* Level 0  = no printing), ... , Level 4 = maximum printing */
    Parm->PrintLevel = 0 ;

    /* T => print parameters values
       F => do not display parmeter values */
    Parm->PrintParms = FALSE ;

    /* T => use approximate nonmonotone Armijo line search
       F => use ordinary nonmonotone Armijo line search, switch to
            approximate Armijo when |f_r-f| < AArmijoFac*|min (f_r, f_{max})| */
    Parm->AArmijo = FALSE ;
    Parm->AArmijoFac = 1.e-8 ;

    /* Stop Rules (these override the corresponding cg parameters):
       T => ||proj_grad||_infty <= max(grad_tol,initial ||grad||_infty*StopFact)
       F => ||proj_grad||_infty <= grad_tol*(1 + |f_k|) */
    Parm->StopRule = TRUE ;
    Parm->StopFac = 0.e-12 ;

    /* T => estimated error in function value = eps*|min (f_r, f_{max}) |
       F => estimated error in function value = eps */
    Parm->PertRule = TRUE ;
    Parm->eps = 1.e-6 ;

    /* T => only use gradient projection algorithm
       F => let algorithm decide between grad_proj and cg_descent */
    Parm->GradProjOnly = FALSE ;

    /* abort cbb after maxit_fac*n iterations in one pass through cbb */
    Parm->maxit_fac = INF ;

    /* abort cbb after totit_fac*n iterations in all passes through cbb */
    Parm->totit_fac = INF ;

    /* abort cbb iteration after maxfunc_fac*n function evaluations */
    Parm->maxfunc_fac = INF ;

    /* perturbation in bounds based on machine epsilon, which we now compute */
    eps = ONE ;
    t = ONE ;
    while ( t > 0 )
    {
        eps /= TWO ;
        t = ONE + eps ;
        t -= ONE ;
    }
    eps *= 2 ;                   /* machine epsilon */
    Parm->pert_lo = 1.e3*eps ;   /* perturbation of lower bounds */
    Parm->pert_hi = 1.e3*eps ;   /* perturbation of upper bounds */

    /* search for non nan function value by shrinking search interval
       at most nshrink times */
    Parm->nshrink = (int) 50 ;

    /* factor by which interval shrinks when searching for non nan value */
    Parm->nan_fac = 2.e-1 ;

    /* update fr if fmin was not improved after L iterations */
    Parm->L = 3 ;

    /* fmax = max (f_{k-i}, i = 0, 1, ..., min (k, m-1) ) */
    Parm->m = 8 ;

    /* update fr if initial stepsize was accepted in previous P iterations */
    Parm->P = 40 ;

    /* CBB cycle length */
    Parm->nm = 4 ;

    /* Reinitialize BB stepsize, if (s^t y)/(||s|| ||y||) >= gamma
       and ||s|| <= min (parm3*|f_k+1|/||g_k+1||_infty, 1) */
    Parm->gamma = 0.975e0 ;

    /* update reference value fr if (fr-fmin)/(fc-fmin) > gamma1 */
    Parm->gamma1 = (double) Parm->m / (double) Parm->L ;

    /* update fr if (fr-f)/(fmax-f) > gamma2, np > P, and fmax > f */
    Parm->gamma2 = (double) Parm->P / (double) Parm->m ;

    /* terminate Armijo line search when
       phi(alpha) <= phi_r + alpha * delta * phi'(0) where phi_r = fr or fcomp*/
    Parm->delta = 1.0e-4 ;   /* Armijo line search parameter */

    /* stepsize s in the line search must satisfy lmin <= s <= lmax */
    Parm->lmin = 1.0e-20 ;
    Parm->lmax = 1.0e+20 ;

    /* attempt a quadratic interpolation step in cg_descent if the
       provisional stepsize times parm1 <= stepsize to boundary */
    Parm->parm1 = 1.e-1 ;

    /* if quadratic interpolation step is attempted, the provisional step
       is at most parm2*stepsize to boundary */
    Parm->parm2 = 9.e-1 ;

    /* used in the the criterion of reinitializing the BB stepsize */
    Parm->parm3 = 1.e-1 ;

    /* maximum number of previous BB steps used when s^t y <= ZERO */
    Parm->parm4 = 6 ;

    /* if ginorm < tau1*pgnorm, continue gradient projection steps  */
    Parm->tau1 = 1.e-1 ;

    /* decay factor for tau1 */
    Parm->tau1_decay = 5.e-1 ;

    /* ginorm < tau2*pgnorm implies subproblem solved in cgdescent */
    Parm->tau2 = 1.e-1 ;

    /* decay factor for tau2 */
    Parm->tau2_decay = 5.e-1 ;

    /* if pgnorm < pgdecay*MAX (pgnorm0, ONE), check the undecided index set
                                pgnorm0 = pgnorm at starting point */
    Parm->pgdecay = 1.e-4 ;

    /* backtracking decay factor in the Armijo line search */
    Parm->armijo_decay = 5.e-1 ;

    /* use quadratic interpolation to compute Armijo step if it
       lies in the interval [.1 alpha, .9 alpha] */
    Parm->armijo0 = 1.e-1 ;
    Parm->armijo1 = 9.e-1 ;
}

/* =========================================================================
   === asa_cg_default ======================================================
   =========================================================================
   Set default conjugate gradient parameter values. If the parameter argument
   of asa_cg is NULL, this routine is called by asa_cg automatically.
   If the user wishes to set parameter values, then the asa_parameter structure
   should be allocated in the main program. The user could call asa_cg_default
   to initialize the structure, and then individual elements in the structure
   could be changed, before passing the structure to asa_cg.
   =========================================================================*/
void asa_cg_default
(
    asacg_parm   *Parm
)
{
    /* Level 0 = no printing, ... , Level 4 = maximum printing */
    Parm->PrintLevel = 0 ;

    /* T => print parameters values
       F => do not display parmeter values */
    Parm->PrintParms = FALSE ;

    /* T => use approximate Wolfe line search
       F => use ordinary Wolfe line search, switch to approximate Wolfe when
                |f_k+1-f_k| < AWolfeFac*C_k, C_k = average size of cost */
    Parm->AWolfe = FALSE ;
    Parm->AWolfeFac = 1.e-3 ;

    /* T => estimated error in function value is eps*Ck,
       F => estimated error in function value is eps */
    Parm->PertRule = TRUE ;
    Parm->eps = 1.e-4 ;

    /* T => attempt quadratic interpolation in line search when
                |f_k+1 - f_k|/f_k <= QuadCutoff
       F => no quadratic interpolation step */
    Parm->QuadStep = TRUE ;
    Parm->QuadCutOff = 1.e-12 ;

    /* T => check that f_k+1 - f_k <= debugtol*C_k
       F => no checking of function values */
    Parm->debug = FALSE ;
    Parm->debugtol = 1.e-10 ;

    /* factor in [0, 1] used to compute average cost magnitude C_k as follows:
       Q_k = 1 + (Qdecay)Q_k-1, Q_0 = 0,  C_k = C_k-1 + (|f_k| - C_k-1)/Q_k */
    Parm->Qdecay = .7 ;

    /* if step is nonzero, it is the initial step of the initial line search */
    Parm->step = ZERO ;

    /* abort cg after maxit_fac*n iterations in one pass */
    Parm->maxit_fac = INF ;

    /* abort cg after totit_fac*n iterations in all passes */
    Parm->totit_fac = INF ;

    /* maximum number of times the bracketing interval grows or shrinks
       in the line search is nexpand */
    Parm->nexpand = (int) 50 ;

    /* maximum number of secant iterations in line search is nsecant */
    Parm->nsecant = (int) 50 ;

    /* conjugate gradient method restarts after (n*restart_fac) iterations */
    Parm->restart_fac = ONE ;

    /* stop when -alpha*dphi0 (estimated change in function value) <= feps*|f|*/
    Parm->feps = ZERO ;

    /* after encountering nan, growth factor when searching for
       a bracketing interval */
    Parm->nan_rho = 1.3 ;

    /* Wolfe line search parameter, range [0, .5]
       phi (a) - phi (0) <= delta phi'(0) */
    Parm->delta = .1 ;

    /* Wolfe line search parameter, range [delta, 1]
       phi' (a) >= sigma phi' (0) */
    Parm->sigma = .9 ;

    /* decay factor for bracket interval width in line search, range (0, 1) */
    Parm->gamma = .66 ;

    /* growth factor in search for initial bracket interval */
    Parm->rho = 5. ;

    /* conjugate gradient parameter beta_k must be >= eta*||d_k||_2 */
    Parm->eta = .01 ;

    /* starting guess for line search =
         psi0 ||x_0||_infty over ||g_0||_infty if x_0 != 0
         psi0 |f(x_0)|/||g_0||_2               otherwise */
    Parm->psi0 = .01 ;      /* factor used in starting guess for iteration 1 */

    /* for a QuadStep, function evalutated at psi1*previous step */
    Parm->psi1 = .1 ;

    /* when starting a new cg iteration, our initial guess for the line
       search stepsize is psi2*previous step */
    Parm->psi2 = 2. ;
}

/* =========================================================================
   === asa_descent =========================================================
   =========================================================================
   cg_descent conjugate gradient algorithm with modifications to handle the
   bound constraints.