Optimization with Analytic Gradient Computations

This GAUSS optimization example demonstrates nonlinear least squares estimation. The data for this example is taken from G.P.Y Cleark, “Approximate Confidence Limits for a Parameter Function in Nonlinear Regression”, JASA, 82:221-230, 1987.

The key feature of this example is the use of user-defined analytical first derivatives.

Code for optimization

The example:

• Creates a data matrix.

• Defines an objective function that takes 5 inputs.

• Specifies analytical computations for the gradients.

new;
cls;

// Make optmt library available
library optmt;

/*
** data  -  weight of cut grass from 10 randomly sited quadrants taken
**          each week for 13 weeks from a grazing pasture.
*/
y = { 3.183,
    3.059,
    2.871,
    2.622,
    2.541,
    2.184,
    2.110,
    2.075,
    2.018,
    1.903,
    1.770,
    1.762,
    1.550 };

x = seqa(1, 1, 13);

// number of observations
n = 10;

// The Mitcherlitz Equation
proc Mitcherlitz(b, x);

    // Calculate and return Mitcherlitz equation
    retp(b[3] + b[2]*exp(-b[1]*x));
endp;

// The Gradient
proc dMitcherlitz(b, x);
    local w;

    w = exp(-b[1]*x);
    retp(((-b[2]*x.*w)~w~ones(rows(w), 1)));
endp;


// The objective function
proc fct(b, x, y, nobs, ind);
    local dev;
    struct modelResults mm;

    dev = y - Mitcherlitz(b,x);

    // If first element of 'ind' is non-zero,
    // calculate and return objective function
    if ind[1];
        mm.function = dev'dev/nobs;
    endif;

    // If the second element of 'ind' is non-zero,
    // calculate and return gradient
    if ind[2];
        mm.gradient = -2*(dMitcherlitz(b,x)'dev)/nobs;
    endif;

    retp(mm);
endp;

// 3x1 vector of starting values
b0 = { 1, 1, 0 };

// Declare 'out' to be an optmtResults structure
// to hold the results of the optimization
struct optmtResults out;

// Optimize function
out = optmt(&fct, b0, x, y, n);

// Print summary report of optimization results
call optmtPrt(out);

The code prints results to the Command Window.

Results

Convergence details

The first portion of the results provide details about convergence and performance.

Return code    =    0
Function value =    0.00535
Convergence    :    normal convergence

These results indicate that the optimization converged normally, with a return code of 0. Any return Code other than 0 would indicate some issue with the convergence. The exact meaning of the return code can be found in the optmt() documentation.

Parameter estimates

The next section of the results reports the parameter estimates and the associated gradients.

Parameters  Estimates   Gradient
---------------------------------------------------------------------
x[1,1]      0.1031      0.0000
x[2,1]      2.5190      0.0000
x[3,1]      0.9631      0.0000

In this example, the gradients for parameters all parameters are zero, as is expected at or near an optimum.

Computation time

Number of iterations    28
Minutes to convergence     0.00013

Extension

Though the printed results do not include the standard errors of the parameters or the other key inference statistics, these can be easily computed the dMitcherlitz function and the results stored in the optmtResults structure.

// Extract estimated parameters
bhat = pvGetParvector(out.par);

// Compute gradients
grad = dMitcherlitz(bhat,x);

// Find covariance matrix
cov = out.fct*invpd(grad'*grad);

// Print results
print;
print "standard errors of parameters";
sd = sqrt(diag(cov));
print sd';
print;
print "Correlation matrix of parameters";
print cov./sd./sd';
print;
print "t-statistics";
print (bhat./sd)';

The results printed are:

standard errors of parameters
     0.025504306       0.26577209       0.32158853

Correlation matrix of parameters

       1.0000000     -0.92330571       0.98410596
     -0.92330571       1.0000000      -0.97238973
      0.98410596     -0.97238973        1.0000000

t-statistics
      4.0406689        9.4780502        2.9948710