derivtst.gms : How to test derivatives of functions

Description

This model demonstrates how to tests the first and second order derivatives
of a given function by comparing the (exact/analytical) derivatives with
derivatives calculated by finite differences.

GAMS provides through the function suffix .grad/.hess the (exact) derivatives
and the finite difference derivatives through function suffix gradn/.hessn.
These facilities work for intrinisc and extrinsic functions.

Reference

  • GAMS Development Corporation, Formulation and Language Example.

Small Model of Type : GAMS


Category : GAMS Model library


Main file : derivtst.gms

$Title How to test derivatives of functions (derivtst,SEQ=406)

$Ontext

This model demonstrates how to tests the first and second order derivatives
of a given function by comparing the (exact/analytical) derivatives with
derivatives calculated by finite differences.

GAMS provides through the function suffix .grad/.hess the (exact) derivatives
and the finite difference derivatives through function suffix gradn/.hessn.
These facilities work for intrinisc and extrinsic functions.

$Offtext

$if not set sampleSize $set sampleSize 100
$if not set tol        $set tol        1e-6
$if not set func       $set func       power

Set s     Sample             / s1*s%sampleSize% /
    ai    Argument Info      / lo, up, int, endo /

*** Function data
$iftheni %func%==pdfnormal
$  funclibin mylib stodclib
   function %func% /mylib.%func% /;
$  set args sp(s,'a1'),sp(s,'a2'),sp(s,'a3')
   Set a Function Arguments / a1*a3 /
   Table
       fi(a,ai) Function argument information
                   lo   up   int endo
             a1   -10   10          1
             a2    -5    5
             a3   0.5    1;
$elseifi %func%==power
$  set args sp(s,'a1'),sp(s,'a2')
   Set a Function Arguments / a1*a2 /
   Table
       fi(a,ai) Function argument information
                   lo   up   int endo
             a1  -100  100          1
             a2    -5    5     1;
$else
$  abort 'no function data for %func%
$endif
***

Parameter
    sp(s,a)      Sample points for argument
    grad(s,a)    Exact gradient
    gradn(s,a)   Numerical gradient
    hess(s,a,a)  Exact gradient
    hessn(s,a,a) Numerical gradient
Set
    ae(a)        Endogenos arguments
    badsp(s)     Sample points that trigger an execerror
;
Alias (ae,aep);

ae(a) = fi(a,'endo');
sp(s,a) =   uniformInt(fi(a,'lo'),fi(a,'up'))$(    fi(a,'int'))
          +    uniform(fi(a,'lo'),fi(a,'up'))$(not fi(a,'int'));

loop(s,
  grad (s,ae) = %func%.grad (ae.pos:%args%);
  gradn(s,ae) = %func%.gradn(ae.pos:%args%);
  hess (s,ae,aep) = %func%.hess (ae.pos:aep.pos:%args%);
  hessn(s,ae,aep) = %func%.hessn(ae.pos:aep.pos:%args%);
  if (execerror,
    badsp(s) = yes;
    execerror = 0;
  )
);

* Compare gradient and hessian
Set csp(s)           Sample points to compare
Parameter
    badargs(s,a)     Sample points that triggered an execution error
    gradError(s,a)   Error in gradient calculation
    HessError(s,a,a) Error in hessian  calculation;

csp(s) = not badsp(s);
badargs(badsp,a) = sp(badsp,a);
$macro Error(a,b) abs(a-b)/abs(a+1)
gradError(csp,ae)     = Error(grad(csp,ae),gradn(csp,ae));
hessError(csp,ae,aep) = Error(hess(csp,ae,aep),hessn(csp,ae,aep));
gradError(csp,ae    )$(gradError(csp,ae)    <%tol%) = 0;
hessError(csp,ae,aep)$(hessError(csp,ae,aep)<%tol%) = 0;

abort$(card(badsp)+card(gradError)+card(hessError)) badargs,gradError,hessError;