hessp.src

来自「没有说明」· SRC 代码 · 共 118 行

/*
** hessp.src
** (C) Copyright 1988-1998 by 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.
**
**> hessp
**
**  Purpose:    Computes the matrix of second partial derivatives
**              (Hessian matrix) of a function defined by a procedure.
**
**  Format:     h = hessp( &f, x0 );
**
**  Inputs:     &f    pointer to a single-valued function f(x), defined
**                    as a procedure, taking a single Kx1 vector
**                    argument (f:Kx1 -> 1x1).  It is acceptable for
**                    f(x) to have been defined in terms of global
**                    arguments in addition to x:
**
**                       proc f(x);
**                           retp( exp(x'b) );
**                       endp;.
**
**            x0      Kx1 vector specifying the point at which the Hessian
**                    of f(x) is to be computed.
**
**  Output:   h       KxK matrix of second derivatives of f with respect
**                    to x at x0. This matrix will be symmetric.
**
**  Remarks:    This procedure requires K(K+1)/2 function evaluations. Thus
**              if K is large it may take a long time to compute the
**              Hessian matrix.
**
**              No more than 3 - 4 digit accuracy should be expected from
**              this function, though it is possible for greater accuracy
**              to be achieved with some functions.
**
**              It is important that the function be properly scaled, in
**              order to obtain greatest possible accuracy. Specifically,
**              scale it so that the first derivatives are approximately
**              the same size.  If these derivatives differ by more than a
**              factor of 100 or so, the results can be meaningless.
*/

proc hessp(&f,x0);
    local k, hessian, grdd, ax0, dax0, dh, xdh, ee, f0, i, j, eps;
    local f:proc;

    /* check for complex input */
    if iscplx(x0);
        if hasimag(x0);
            errorlog "ERROR: Not implemented for complex matrices.";
            end;
        else;
            x0 = real(x0);
        endif;
    endif;

/* initializations */
    k = rows(x0);
    hessian = zeros(k,k);
    grdd = zeros(k,1);
    eps = 6.0554544523933429e-6;

/* Computation of stepsize (dh) */
    ax0 = abs(x0);
    if x0 /= 0;
        dax0 = x0./ax0;
    else;
        dax0 = 1;
    endif;
    dh = eps*maxc((ax0~(1e-2)*ones(rows(x0),1))').*dax0;

    xdh = x0+dh;
    dh = xdh-x0;    /* This increases precision slightly */
    ee = eye(k).*dh;

/* Computation of f0=f(x0) */
    f0 = f(x0);

/* Compute forward step */
    i = 1;
    do until i > k;

        grdd[i,1] = f(x0+ee[.,i]);

        i = i+1;
    endo;

/* Compute "double" forward step */
    i = 1;
    do until i > k;
        j = i;
        do until j > k;

            hessian[i,j] = f(x0+(ee[.,i]+ee[.,j]));
            if i /= j;
                hessian[j,i] = hessian[i,j];
            endif;

            j = j+1;
        endo;
        i = i+1;
    endo;

    retp( ( ( (hessian - grdd) - grdd') + f0) ./ (dh.*dh') );
endp;