qnewton.src

没有说明

/*
**  qnewton.src    Quasi-Newton minimization
**
** (C) Copyright 1996-1997  Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
*/


/* >proc QNewton
**
**  Format:  { x,f,g,retcode } = QNewton(&fct,x0)
**
**  Input:   &fct     pointer to a procedure that computes the function to
**                    be minimized.  This procedure must have one input
**                    argument, a vector of parameter values, and one
**                    output argument, the value of the function evaluated
**                    at the input vector of parameter values.
**
**             x0     Kx1 vector of start values
**
**  Output:     x     Kx1 vector of parameters at minimum
**
**              f     scalar function evaluated at x
**
**              g     Kx1 gradient evaluated at x
**
**        retcode     return code:
**
**                       0   normal convergence
**                       1   forced exit
**                       2   maximum number of iterations exceeded
**                       3   function calculation failed
**                       4   gradient calculation failed
**                       5   step length calculation failed
**                       6   function cannot be evaluated at initial
**                           parameter values
**
**  Globals:  _qn_RelGradTol   scalar, convergence tolerance for relative
**                             gradient of estimated coefficients.
**                             Default = 1e-5.
**
**            _qn_GradProc     scalar, pointer to a procedure that computes
**                             the gradient of the function with respect to
**                             the parameters. This procedure must have a
**                             single input argument, a Kx1 vector of parameter
**                             values, and a single output argument, a Kx1
**                             vector of gradients of the function with respect
**                             to the parameters evaluated at the vector of
**                             parameter values.
**
**            _qn_MaxIters     scalar, maximum number of iterations.
**                             Default = 1e+5. Termination can be forced
**                             by pressing C on the keyboard.
**
**            _qn_PrintIters   scalar, if 1, print iteration information.
**                             Default = 0.  Can be toggled during iterations
**                             by pressing P on the keyboard.
**
**            _qn_ParNames     Kx1 vector, labels for parameters
**
**          _qn_RandRadius     scalar, If zero, no random search is attempted.  If nonzero
**                             it is the radius of random search which is invoked whenever
**                             the usual line search fails.  Default = .01.
**
**            __output         scalar, if 1, prints results.
**
**
**  Remarks:  Pressing C on the keyboard will terminate iterations, and
**            pressing P will toggle iteration output.
**
**            QNewton is recursive, that is, it can call a version of itself
**            with another function and set of global variables,
*/

#include gauss.ext
#include qnewton.ext


proc (4) = QNewton(fnct,x0);
    local x,f,g,ret,iters;
    local lbl,mask,fmt,t0,t1;

    if __output;
        t0 = hsec;
    endif;

    { x,f,g,ret,iters } = _bfgs(fnct,x0,_qn_MaxIters,_qn_RelGradTol,
                             _qn_GradProc,_qn_PrintIters,_qn_RandRadius);


    if __output;
        t1 = hsec;
        print;
        call header("QNewton","",_rtl_ver);
        print;
        print "return code = " ftos(ret,"%*.*lf",4,0);
        if ret == 0;
           print "normal convergence";
        elseif ret == 1;
           print "forced termination";
        elseif ret == 2;
           print "maximum number of iterations exceeded";
        elseif ret == 3;
           print "function calculation failed";
        elseif ret == 4;
           print "gradient calculation failed";
        elseif ret == 5;
           print "step length calculation failed";
        elseif ret == 6;
           print "function cannot be evaluated at initial parameter values";
        elseif ret == 7;
           print "maximum time exceeded";
        endif;
        print;
        print "Value of objective function " ftos(f,"%*.*lf",15,6);
        print;
        print "Parameters    Estimates    Gradient";
        print "-----------------------------------------";

        if rows(g) /= rows(x);
             g = miss(zeros(rows(x),1),0);
        endif;

        if _qn_ParNames $== 0;
            lbl = 0 $+ "P" $+ ftocv(seqa(1,1,rows(x)),2,0);
        else;
            lbl = _qn_ParNames;
        endif;
        mask = 0~1~1;
        let fmt[3,3] = "-*.*s" 9 8 "*.*lf" 14 4 "*.*lf" 14 4;
        call printfm(lbl~x~g,mask,fmt);

        print;
        print "Number of iterations    " ftos(iters,"%*.*lf",5,0);
        print "Minutes to convergence  " ftos((t1-t0)/6000,"%*.*lf",10,5);
    endif;

    retp(x,f,g,ret);
endp;






proc (5) = _bfgs(fnct,x0,Lmaxiters,Lrelgradtol,Lgradproc,Lprintiters,
                 Lrteps);

    local j,h,g0,g1,d,f0,f1,f2,oldt,np,relgrad,ky,gproc,vv,zz,iter,
          delta,ub,lb,cdelta,dg,s,tt,sprev,s2prev,rprev,r2prev,
          rteps,ab,a,b,sprev2,s2prev2,sp2,dsprev,qv,w1,v1,dx,maxtry,
          fnct:proc;

    vv = ones(2,2);
    vv[1,2] = -1;
    zz = zeros(2,1);
    sp2 = zeros(1,2);
    MaxTry = 100;

    if Lgradproc /= 0;
        gproc = Lgradproc;
        local gproc:proc;
    endif;

    np = rows(x0);
    f0 = fnct(x0);
    if scalmiss(f0) or (f0 $== __INFp) or (f0 $== __INFn) or (f0 $==
         __INDEFp) or (f0 $== __INDEFn);
            retp(error(0),error(0),error(0),6,0);
    endif;
    if Lgradproc /= 0;
        g0 = gproc(x0);
    else;
        g0 = gradp(&fnct,x0)';
    endif;

    if scalmiss(g0) or (g0 $== __INFp) or (g0 $== __INFn) or (g0 $==
         __INDEFp) or (g0 $== __INDEFn);
            retp(error(0),error(0),error(0),4,0);
    endif;

    h = eye(np)*maxc(sqrt(abs(f0))|1);

#ifUNIX
    for iter(1,Lmaxiters,1);
#else
    iter = 1; do until iter > Lmaxiters;
#endif

        /*  new direction */

        oldt = trapchk(1);
        trap 1,1;
        d = -cholsol(g0,h);
        trap oldt,1;
        if scalmiss(d);
            h = eye(np)*maxc(sqrt(abs(f0))|1);
            d = -g0 / maxc(sqrt(abs(f0))|1);
        endif;

        /* choose step length */


        delta = 1e-4;   /* must be in (0,1/2) interval */
        ub = 0.5;       /* Upper bound on acceptable reduction in s. */
        lb = 0.1;       /* Lower bound on acceptable reduction in s. */
        cdelta = 1 - delta;
        dg = d' * g0;
        s = 1;
        tt = s * dg;
        f1 = fnct(x0 + d);
        if scalmiss(f1) or (f1 $== __INFp) or (f1 $== __INFn) or (f1 $==
             __INDEFp) or (f1 $== __INDEFn);
              retp(error(0),error(0),error(0),3,iter);
        endif;

        if (f1 / tt - f0 / tt) < delta;

            s = maxc(-dg / (2 * (f1 - f0 - dg)) | lb);
            f2 = fnct(x0 + s*d);
            if scalmiss(f2) or (f2 $== __INFp) or (f2 $== __INFn) or (f2 $==
                 __INDEFp) or (f2 $== __INDEFn);
                    retp(error(0),error(0),error(0),5,iter);
            endif;


            tt = s * dg;
            w1 = f2 / tt - f0 / tt;

            if w1 < delta or w1 > cdelta;

                sprev = s;
                s2prev = 1;
                rprev = f2;
                r2prev = f1;
#ifUNIX
                for j(1,40,1);   /*  maxtries = 40 */
#else
                j = 1; do until j > MaxTry;
#endif
                    sprev2 = sprev * sprev;
                    s2prev2 = s2prev * s2prev;
                    sp2[1,1] = sprev2;
                    sp2[1,2] = s2prev2;
                    dsprev = sprev - s2prev;

                    vv[2,1] = -s2prev;
                    vv[2,2] = sprev;
                    vv = vv./sp2;

                    zz[1] = rprev - f0 - dg * sprev;
                    zz[2] = r2prev - f0 - dg * s2prev;
                    ab = (1 / dsprev) * vv * zz;
                    a = ab[1,1];
                    b = ab[2,1];

                    if a == 0;
                        s = -dg / (2 * b);
                    else;
                        qv = b * b - 3 * a * dg;
                        if qv < 0;
                            break;
                        endif;          /* terminate if not real root */
                        tt = 3 * a;
                        s = -b / tt + sqrt(qv) / tt;
                    endif;

                    if s > ub * sprev;
                        s = ub * sprev;
                    elseif s < lb * sprev;
                        s = lb * sprev;
                    endif;

                    f1 = fnct(x0 + s * d);
                    if scalmiss(f1) or (f1 $== __INFp) or (f1 $== __INFn) or
                       (f1 $==  __INDEFp) or (f1 $== __INDEFn);
                           s = error(0);
                           break;
                    endif;
                    tt = s * dg;
                    w1 = f1 / tt - f0 / tt;

                    if w1 >= delta and w1 <= cdelta;
                        break;
                    endif;

                    s2prev = sprev;
                    sprev = s;
                    r2prev = rprev;
                    rprev = f1;
#ifUNIX
                endfor;
#else
                j = j + 1; endo;
#endif

            else;
                 f1 = f2;
            endif;

        endif;

        if scalmiss(s);
            j = 1;
            f1 = 1e250;
            do while f1 > f0;
                rteps = 10^trunc(log(meanc(abs(g0)))) * Lrteps;
                d = rteps*(2*rndu(np,1)-1).*x0;
                f1 = fnct(x0+d);
                if scalmiss(f1) or (f1 $== __INFp) or (f1 $== __INFn) or
                   (f1 $==  __INDEFp) or (f1 $== __INDEFn);
                   f1 = 1e250;
                   d = 1;
                endif;
                j = j + 1;
                if j > MaxTry;
                    break;
                endif;
            endo;
            x0 = x0 + d;
            f0 = f1;
        else;
            dx = s * d;
            x0 = x0 + dx;
            f0 = f1;
        endif;


        if Lprintiters;
            print;
            print " step length " s;
        endif;

        if Lgradproc /= 0;
            g1 = gproc(x0);
        else;
            g1 = gradp(&fnct,x0)';
        endif;
        if scalmiss(g1) or (g1 $== __INFp) or (g1 $== __INFn) or (g1 $==
             __INDEFp) or (g1 $== __INDEFn);
                retp(error(0),error(0),error(0),4,iter);
        endif;

        relgrad = (abs(g1).*maxc(abs(x0)'|ones(1,np)))/maxc(abs(f0)|1);

        if Lprintiters;
            print;
            print;
            print "---------------------------------------------------";
#ifUNIX
            print " iter "$+ftos(iter+0,"%*.*lf",1,0);;
#else
            print " iter "$+ftos(iter,"%*.*lf",1,0);;
#endif
            print "          function = "$+ftos(f0,"%*.*lf",15,8);
            print "---------------------------------------------------";
            print;
            print "       parameter      direction        gradient"\
                     "       relative gradient";
            call printfmt(x0~d~g1~relgrad,1);
#ifDLLCALL
            print /flush "";
#else
            print;
#endif
        endif;

        if abs(g1) < __macheps or relgrad < Lrelgradtol;
            retp(x0,f0,g1,0,iter);
        endif;

        /* BFGS update */

        v1 = g1'dx - g0'dx;
        if (v1 < 1e-22);
            h = eye(np)*maxc(sqrt(abs(f0))|1);
        else;
            oldt = trapchk(1);
            trap 1,1;
            h = cholup(h,(g1-g0)/sqrt(v1));
            trap oldt,1;
            if not scalmiss(h);
                oldt = trapchk(1);
                trap 1,1;
                h = choldn(h,g0/sqrt(-g0'd));
                trap oldt,1;
             endif;
            if scalmiss(h);
                h = eye(np)*maxc(sqrt(abs(f0))|1);
             endif;
        endif;

        g0 = g1;

        ky = key;
        do while ky;
            if upper(chrs(ky)) $== "C";
                retp(x0,f0,g0,1,iter);
            elseif upper(chrs(ky)) $== "P";
                Lprintiters = 1 - Lprintiters;
            endif;
            ky = key;
        endo;
#ifUNIX
    endfor;
#else
    iter = iter + 1; endo;
#endif

    retp(x0,f0,g0,2,Lmaxiters);

endp;


proc(0) = qnewtonset;
     gausset;
     _qn_MaxIters = 1e5;
     _qn_RelGradTol = 1e-5;
     _qn_GradProc = 0;
     _qn_PrintIters = 0;
     _qn_RandRadius = .01;
endp;