# Optimization with Analytic Gradient Computations#

This GAUSS optimization example demonstrates the use of COMT to estimate parameters using the sum of the squared residuals. This example moves beyond basic optimization and uses a user-defined function to compute the analytical first derivatives.

## Key example features#

• The user-defined objective function `ssq` for computing the sum of the squared residuals.

• The user-defined gradient function `grd` used for computing the first derivatives.

### Code for optimization#

```new;
cls;

library comt;

// Create dependent and independent variables
y = { 109,
149,
149,
191,
213,
224 };

x = { 1,
2,
3,
5,
7,
10 };

// Starting parameter values
b_strt =  { 100, 0.75 };

// Objective procedure, to compute sum of squared residuals, with 4 inputs:
//    i.      x       - The parameter vector.
//    ii-iii. y and x - Dependent and independent variable
//                      for nonlinear least squares estimation.
//    ii.     ind     - The indicator vector.

proc (1) = ssq(b, y, x, ind);
// Declare 'mm' to be a modelResults
// struct local to this procedure
struct modelResults mm;

local f,g;

// Calculate 'f' once, for use
// in computation of function value
f = b[1] * (1 - exp(-b[2]*x));

// If first element of 'ind' is non-zero,
// compute function evaluation
if ind[1];
mm.function = (y - f)'(y - f);
endif;

// If second element of 'ind' is non-zero,
if ind[2];
g = grd(b, y, x);
mm.gradient = -2 * g'(y - f);
endif;

// Return modelResults struct
retp(mm);
endp;

// Procedure used in gradient computation
// Note this procedure must take the same inputs as
// the objective function (excluding 'ind')
proc grd(b, y, x);
local f,g;

g = zeros(rows(x),rows(b));

f = exp(-b[2]*x);

g[.,1] = 1 - f;
g[.,2] = b[1]*x.*f;

retp(g);
endp;

// Declare 'ctl' to be a comtControl struct
// and fill with default settings
struct comtControl ctl;
ctl = comtControlCreate();

// Print numerical and analytical gradient
// on first iteration for comparison and debugging

// Declare 'out' to be a comtResults struct
// to hold the results from the optimization
struct comtResults out;

// Minimize objective function
// and print output to the screen
out = comtPrt(comt(&ssq, b_strt, y, x, ctl));
```

The code prints results to the Command Window.

### Results#

The first notable component of the results is the printed comparison of the analytical and numerical gradients. This is done because the ctl.gradCheck member of the `comtControl` structure was set to 1.

```analytical gradient        numerical gradient
-932.8196153813061073     -932.8196069540683766
-18609.6357937713837600     -18609.6345830073878460
```

## Convergence details#

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

```Return code    = 0
Function value = 1168.00888
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 `comt()` documentation.

## Parameter estimates#

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

```Parameters  Estimates    Gradient
---------------------------------------------------------------------
x[1,1]       213.8094      0.0000
x[2,1]         0.5472      0.0002
```

## Computation time#

The final section of the results reports the number of iterations and computation time.

```Number of iterations    2572
Minutes to convergence  0.00793
```

### 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 `grd` function and the results stored in the `comtResults` structure.

```// Extract estimated parameters
b_hat = pvGetParVector(out.par);

// Calculate gradient at the estimated parameters
g = grd(b_hat, y, x);

// Compute covariance matrix using
// the final function value and
// the gradient at estimated parameters
cov = (out.fct/rows(y))*invpd(g'*g);

// Find standard errors of estimates
sd = sqrt(diag(cov));
print "Standard errors of parameters:";
print sd';

// Find t-stats of estimates
t_stats = (b_hat./sd)';
print "Parameter t-statistics:";
print t_stats;
```

This prints the standard errors and t-statistics.

```Standard errors of parameters:
10.087419      0.085372828
Parameter t-statistics:
21.195650        6.4099726
```