comt#

Purpose#

Solve the nonlinear programming problem.

Format#

out = comt(&modelProc, par[, ..., c1])#

Parameters:

&modelProc (function pointer) – Pointer to a procedure that computes the function to be minimized.
par (struct) – An instance of a PV structure. The par instance is passed to the user-provided procedure pointed to by &modelProc. par is constructed using the “pack” functions.
... (various) – Optional input arguments. They can be any set of structures, matrices, arrays, strings required to compute the function. Can include GAUSS data types or a DS structure for dataset manipulation. Specific usage depends on the requirements of the modelProc.

c1 (struct) –

Optional input. Instance of a comtControl structure containing the following members:

c1.A	MxK matrix, linear equality constraint coefficients: `c1.A * p = c1.B` where `p` is a vector of the parameters.
c1.B	Mx1 vector, linear equality constraint constants: `c1.A * p = c1.B` where `p` is a vector of the parameters.
c1.C	MxK matrix, linear inequality constraint coefficients: `c1.C * p >= c1.D` where `p` is a vector of the parameters.
c1.D	Mx1 vector, linear inequality constraint constants: `c1.C * p >= c1.D` where `p` is a vector of the parameters.
c1.eqProc	Scalar, pointer to a procedure that computes the nonlinear equality constraints. It has two input arguments, an instance of a PV parameter structure, and an instance of a DS data structure, and one output argument, a vector of computed equality constraints. Default = {}, i.e., no equality procedure.
c1.ineqProc	Scalar, pointer to a procedure that computes the nonlinear inequality constraints. It has two input arguments, an instance of a PV parameter structure, and an instance of a DS data structure, and one output argument, a vector of computed inequality constraints. Default = {}, i.e., no inequality procedure.
c1.eqJacobian	Scalar, pointer to a procedure that computes the Jacobian of the equality constraints. It has two input arguments, an instance of a PV parameter structure, and an instance of a DS data structure, and one output argument, a matrix of derivatives of the equality constraints with respect to the parameters. Default = {}, i.e., no equality Jacobian procedure.
c1.ineqJacobian	Scalar, pointer to a procedure that computes the Jacobian of the inequality constraints. It has two input arguments, an instance of a PV parameter structure, and an instance of a DS data structure, and one output argument, a matrix of derivatives of the inequality constraints with respect to the parameters. Default = {}, i.e., no inequality Jacobian procedure.
c1.bounds	1x2 or Kx2 matrix, bounds on parameters. If 1x2 all parameters have same bounds. Default = {-1e256, 1e256}.
c1.algorithm	Scalar, descent algorithm, 0 - Modified BFGS (Default), 1 = BFGS, 2 = DFP, 3 = Newton.
c1.useThreads	Scalar, if nonzero threading is turned on, else off. Default = off.
c1.switch	4x1 or 4x2 vector, controls algorithm switching. If 4x1, specifics include algorithm number to switch to, COMT switches if function changes less than a certain amount, COMT switches if a certain number of iterations is exceeded, COMT switches if line search step changes less than a certain amount. Default = {1 3, .0001 .0001, 10 10, .0001 .0001}.
c1.lineSearch	Scalar, sets line search method, 0 augmented Lagrangian penalty method (requires constraints), 1 STEPBT (quadratic and cubic curve fit) (default), 2 Brent’s method, 3 half, 4 Strong Wolfe’s condition.
c1.trustRadius	Scalar, radius of the trust region. If scalar missing, trust region not applied. Sets a maximum amount of the direction at each iteration. Default = 0.1.
c1.penalty	Scalar, augmentation constant for augmented Lagrangian penalty line search method.
c1.active	Kx1 vector, set K-th element to zero to fix it to start value. Use the GAUSS function pvGetIndex to determine where parameters in the PV structure are in the vector of parameters. Default = {}, all parameters are active.
c1.maxIters	Scalar, maximum number of iterations. Default = 10000.
c1.dirTol	Scalar, convergence tolerance. Iterations cease when all elements of the direction vector are less than this value. Default = 1e-5.
c1.feasibleTest	Scalar, if nonzero, parameters are tested for feasibility before computing function in line search. If function is defined outside inequality boundaries, then this test can be turned off. Default = 1.
c1.maxTries	Scalar, maximum number of attempts in random search. Default = 100.
c1.randRadius	Scalar, if zero, no random search is attempted. If nonzero, it is the radius of the random search. Default = .001.
c1.gradMethod	Scalar, method for computing numerical gradient. 0 = central difference, 1 = forward difference (default), 2 = backward difference.
c1.hessMethod	Scalar, method for computing numerical Hessian. 0 = central difference, 1 = forward difference (default), 2 = backward difference.
c1.gradStep	Scalar or Kx1, increment size for computing numerical derivatives (gradient and/or Hessian). If scalar, stepsize will be gradStep value times parameter estimates for the numerical derivative. If Kx1, the step size for the numerical derivatives will be the literal values of the elements of the gradStep vector.
c1.gradCheck	Scalar, if nonzero and if analytical gradients and/or Hessian have been provided, numerical gradients and/or Hessian are computed and compared against the analytical versions.
c1.state	Scalar, seed for random number generator.
c1.title	String, title of run.
c1.printIters	Scalar, if nonzero, prints iteration information. Default = 0.
c1.disableKey	Scalar, if nonzero, keyboard input disabled.

Returns:

out (struct) –

An instance of a comtResults structure. Contains the results of the nonlinear programming problem solution, including parameter estimates, function evaluations, and detailed information about constraints handling and optimization process. The comtResults structure includes:

out.par	Instance of a PV structure containing the parameter estimates will be placed in the member matrix out.par.
out.fct	Scalar, function evaluated at parameters in par.
out.returnDescription	String, description of return values.
out.hessian	KxK matrix, Hessian evaluated at parameters in par.
out.xproduct	KxK matrix, cross-product of NxK matrix of first derivatives evaluated at parameters in par. Not available if log-likelihood function returns a scalar.
out.gradient	Kx1 vector, gradient evaluated at the parameters in par.
out.numIterations	Scalar, number of iterations.
out.elapsedTime	Scalar, elapsed time of iterations.
out.title	String, title of run.
out.lagr	An instance of a comtLagrange structure containing the Lagrangeans for the constraints. For an instance named lagr, the members are: out.lagr.lineq - Mx1 vector, Lagrangeans of linear equality constraints. out.lagr.nlineq - Nx1 vector, Lagrangeans of nonlinear equality constraints. out.lagr.linIneq - Px1 vector, Lagrangeans of linear inequality constraints. out.lagr.nlinineq - Qx1 vector, Lagrangeans of nonlinear inequality constraints. out.lagr.bounds - Kx2 matrix, Lagrangeans of bounds. out.Lagr.EqCov - (M+N)x(M+N) matrix, covariance matrix of equality constraints. out.Lagr.IneqCov - (P+Q)x(P+Q) matrix, covariance matrix of inequality constraints. Whenever a constraint is active, its associated Lagrangean will be nonzero. For any constraint that is inactive throughout the iterations as well as at convergence, the corresponding Lagrangean matrix will be set to a scalar missing value.
out.retcode	Return code: 0 - Normal convergence. 1 - Forced exit. 2 - Maximum number of iterations exceeded. 3 - Function calculation failed. 4 - Gradient calculation failed. 5 - Hessian calculation failed. 6 - Line search failed. 7 - Functional evaluation failed. 8 - Error with initial gradient. 9 - Error with constraints. 10 - Second update failed. 11 - Maximum time exceeded. 12 - Error with weights. 13 - Quadratic program failed. 14 - Equality constraint Jacobian failed. 15 - Inequality constraint Jacobian failed. 16 - Function evaluated as complex. 20 - Hessian failed to invert.

Examples#

Basic case with default optimization settings#

new;
cls;

// Make comt library available
library comt;

// Create data needed by objective procedure
Q = { 0.78 -0.02 -0.12 -0.14,
    -0.02  0.86 -0.04  0.06,
    -0.12 -0.04  0.72 -0.08,
    -0.14  0.06 -0.08  0.74 };

b = { 0.76, 0.08, 1.12, 0.68 };

// Objective procedure with 4 inputs:
//    i.      x       - The parameter vector
//    ii-iii. Q and b - Extra data needed by the objective procedure
//    ii.     ind     - The indicator vector
proc  qfct(x, Q, b, ind);
    // Declare 'mm' to be a modelResults
    // struct local to this procedure, 'qfct'
    struct modelResults mm;

    // If the first element of the indicator
    // vector is non-zero, compute function value
    // and assign it to the 'function' member
    // of the modelResults struct
    if ind[1];
        mm.function = 0.5*x'*Q*x - x'b;
    endif;

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

// Starting parameter values
x_strt = { 1, 1, 1, 1 };

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

// Minimize objective function
out = comt(&qfct, x_strt, Q, b);

// Print output to the screen
call comtPrt(out);

Optimization with inequality constraints#

In this example two optional user-defined equations are specified, one to compute inequality constraints and one to compute the Jacobian of the inequality constraints.

new;
library comt;


// Load dataset
dat_file = __FILE_DIR $+ "nlin.dat";
nldat = loadd(dat_file);

// Assign dependent and independent variables
y = nldat[.,1];
x = nldat[.,2];

// 3x1 vector of starting values
b_strt = { 0.08, 1.1, 0.2 };

// Obective procedure to calculate
// the sum of the squared residuals
proc (1) = ssq(b, y, x, ind);

    struct modelResults mm;
    local dev;

    dev = y - b[1] + b[2] * exp(-b[3]*x);

    if ind[1];
        mm.function = dev'dev;
    endif;

    retp(mm);

endp;

// Procedure to compute inequality constraints.
// Constrains norm of coefficients to be less than 2
proc (1) = ineqp(b, y, x);
    retp(2 - b'b);
endp;

// Procedure to compute Jacobian of
// the inequality constraints
proc (1) = ineqj(b, y, x);
    retp(-2*b');
endp;

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

// Use 'HALF' line search method
ctl.lineSearch = 3;

// Assign pointers to procedures for
// computing inequality constraints
// and the Jacobian of the inequality
// constraints
ctl.ineqProc = &ineqp;
ctl.ineqJacobian = &ineqj;

// Declare 'out' to be a comtResults struct
// to hold the information returned by 'comt'
struct comtResults out;

// Compute parameter estimates
out = comt(&ssq, b_strt, y, x, ctl);

// Print output from optimization
call comtPrt(out);

Remarks#

The comt() function requires one user-provided procedure. This procedure:
- Must specify the computations of the objective function.
- Optionally can specify the computations for the first and/or second derivatives.
- Requires a vector of parameters or an instance of a PV structure containing the parameters as the first input.
- Any number of optional arguments including structures, matrices, arrays, strings, required to compute the objective function.
- A last input named ind.

There are four other optional procedures, one each for:

Equality constraints.

Inequality constraints.

The Jacobian of the equality constraints.

The Jacobian of the inequality constraints.

The optional procedures take just two arguments: the parameters and any optional arguments that were passed to comt().
The instance of the PV structure is set up using the PV pack procedures, pvPack(), pvPackm(), pvPacks(), and pvPacksm().
The optional arguments passed to the user-provided objective function procedure are untouched. This allows you to pass into your function any information it needs.
The procedures for nonlinear equality and inequality constraints take two input arguments, an instance of a PV parameters structure. For example, to constrain the sum of squares of the regression coefficients to be greater than one, provide the following procedure:
```
proc ineqConst(struct PV par1);
     local b;
     b = pvUnpack(p,"b");
     retp( sumc(b^2) - 1 );
endp;
```
If comt() has been called with optional arguments, then they must be included in the call to ineqConst() as well.