## Supplying 1st Derivatives

Once you start solving real (differentiable) problems you will want to supply *at least* exact first derivatives. These can make a large difference (orders of magnitude) in the convergence, accuracy and computation time of the specified solver, when compared to numerical (finite-difference) derivatives.

This section will detail how to obtain exact first derivatives using OPTI and MATLAB, and how to supply them to the solver.

### Mathematical Definition

#### Scalar Function

The first derivatives of a scalar, multivariable function `f`

(such as that used in a NLP objective) is defined as:

Noting the result is a *row vector*, with dimensions 1 x *n*.

#### Vector Function

The first derivatives of a vector, multivariable function `F`

(such as that used in a NLS objective, or vector of nonlinear constraints) is defined as:

Noting the result is a *matrix*, with dimensions *m* x *n*.

### Naming Conventions

OPTI uses a fairly standard naming convention when referring to first derivatives.

#### Gradient

The gradient is *always* the first derivative of the *objective* function. It is identified by the keyword `grad`

when supplying it as an argument to `opti`

. Note the gradient is always *dense*, as all solvers supplied with OPTI expect it this way.

#### Jacobian

The Jacobian is *always* the first derivative(s) of the *constraint* function(s). It is identified by the keyword `jac`

when supplying it as an argument to `opti`

. The Jacobian may be dense or sparse, depending on the solver being used.

Strictly speaking a gradient should be a vector, while the Jacobian should be a matrix. However using these rules, argument names could be flipped based on the problem being solved! Therefore to keep things standard (and hopefully simple), the above naming strategy has been adopted where names are based on whether they refer to an objective, or constraint.

### Example 1: Numerical Differentiation

As a baseline, let's consider the following scalar, multivariable nonlinear objective function from Hock & Schittkowski problem #71:

This can be created in MATLAB as follows:

fun = @(x) x(1)*x(4)*(x(1) + x(2) + x(3)) + x(3);

In order to numerically differentiate it, we will need to supply a point which we want to find the derivatives about. Then we can use the Intel `djacobi`

method to solve for the first derivatives:

x0 = [1;2;3;4];

% Numerically Differentiate

df = mklJac(fun,x0)

This is the method that OPTI uses automatically if no derivatives are supplied, and the solver requires them. It has the advantage that it works with any function, even if it uses MATLAB Toolbox functions, MEX functions, or other black-boxes. However it is typically slower, and much less accurate than achievable via methods described below.

### Example 2: Automatic Differentiation

One of the most widely used methods of supplying derivatives is via Automatic Differentiation (AD). It is used just like Numerical Differentiation (ND) above, but instead returns exact (to numerical precision) derivatives.

dfa = autoJac(fun,x0)

The OPTI `autoJac`

function uses the `ADIFF`

library to generate exact first derivatives. This is done using an object-orientated approach, overloading many common MATLAB operators and functions to calculate analytical first derivatives as it goes. `ADIFF`

implements a forward AD algorithm, successively applying the chain rule to each operator / function it encounters.

While the AD library supplied with OPTI is fairly comprehensive in the functions it covers, it does not support 2D indexing, and thus is limited to relatively simple functions. It also does not support second derivatives (without running it twice), or automatically generate sparse representations, which several commercial tools do.

### Example 3: Complex Step Differentiation

Another method for obtaining accurate derivatives is to use the complex-step method. This uses the built in complex number handling of MATLAB to accurately (to numerical precision) solve the first derivatives of a function.

dfc = cstepJac(fun,x0)

As with the automatic differentiation method, `cstepJac`

only works on pure MATLAB functions. Also be careful with functions like `abs`

, as well as relationship operators, as they are treated differently with complex numbers. In addition, watch out for using the complex conjugate operator (`'`

) rather than transpose (`.'`

), this is the most common cause of errors with this function.

When used correctly complex step differentiation is faster than all other methods and provides derivatives accurate to numerical precision. It is however a little tricky to implement, given that the underlying function must be valid for complex numbers as well.

### Example 4: Symbolic Differentiation

An alternative to computing derivatives at a point, during runtime, is symbolically calculate them once before solving, then supply them as complete functions to the solver. This requires a symbolic manipulator, and the MATLAB Symbolic Toolbox is used for this task.

grad = symJac(fun)

% Evaluate

dfs = grad(x0)

The above example shows how an anonymous function can be symbolically differentiated, then returned as a usable gradient function. It can then be evaluated at any point, and return an exact derivative.

The downside to this approach is that it can become very computationally expensive to evaluate the full derivative expression, especially second derivatives. On large, complex problems, ND or AD will be faster than evaluating the full, analytical expression.

Note the function `symJac`

is a simple demonstration function for toy problems only. It is severly restricted in the problems it can process. However I am in the development of SymBuilder, currently included (but not well documented) in the OPTI distribution. This allows complex expressions to be symbolically processed and differentiated.

### Example 5: Application to a NLP

Now we have some tools for generating first derivatives, lets apply it to the full HS#71 NLP:

As done in previous examples, let's code it up in MATLAB:

fun = @(x) x(1)*x(4)*(x(1) + x(2) + x(3)) + x(3);

% Nonlinear Constraints

nlcon = @(x) [ x(1)*x(2)*x(3)*x(4);

x(1)^2 + x(2)^2 + x(3)^2 + x(4)^2 ];

cl = [25;40];

cu = [Inf;40];

% Bounds (lb <= x <= ub)

lb = ones(4,1);

ub = 5*ones(4,1);

% Initial Guess

x0 = [1 5 5 1]';

Now let's solve it using each of the strategies described so far to generate the first derivatives.

#### Numerical Differentiation

While done automatically by OPTI, it is explicity written out below.

grad = @(x) mklJac(fun,x);

% Jacobian

jac = @(x) mklJac(nlcon,x);

% Build OPTI Problem

Opt = opti('fun',fun,'grad',grad,'nl',nlcon,cl,cu,'jac',jac,...

'bounds',lb,ub,'x0',x0)

% Solve NLP

[x,fval,exitflag,info] = solve(Opt)

#### Automatic Differentiation

Forward AD in MATLAB:

grad = @(x) autoJac(fun,x);

% Jacobian

jac = @(x) autoJac(nlcon,x);

% Build OPTI Problem

Opt = opti('fun',fun,'grad',grad,'nl',nlcon,cl,cu,'jac',jac,...

'bounds',lb,ub,'x0',x0)

% Solve NLP

[x,fval,exitflag,info] = solve(Opt)

#### Complex Step Differentiation

Complex Step in MATLAB:

grad = @(x) cstepJac(fun,x);

% Jacobian

jac = @(x) cstepJac(nlcon,x);

% Build OPTI Problem

Opt = opti('fun',fun,'grad',grad,'nl',nlcon,cl,cu,'jac',jac,...

'bounds',lb,ub,'x0',x0)

% Solve NLP

[x,fval,exitflag,info] = solve(Opt)

#### Symbolic Differentiation

Note function handles are returned by `symJac`

.

grad = symJac(fun);

% Jacobian

jac = symJac(nlcon);

% Build OPTI Problem

Opt = opti('fun',fun,'grad',grad,'nl',nlcon,cl,cu,'jac',jac,...

'bounds',lb,ub,'x0',x0)

% Solve NLP

[x,fval,exitflag,info] = solve(Opt)

#### Manual Differentiation

There is nothing stopping you entering the derivatives by hand!

grad = @(x) [x(1)*x(4) + x(4)*sum(x(1:3)), x(1)*x(4),...

x(1)*x(4) + 1, x(1)*sum(x(1:3))];

% Jacobian

jac = @(x) [prod(x')./x';

2.*x'];

% Build OPTI Problem

Opt = opti('fun',fun,'grad',grad,'nl',nlcon,cl,cu,'jac',jac,...

'bounds',lb,ub,'x0',x0)

% Solve NLP

[x,fval,exitflag,info] = solve(Opt)

### Example 6: Exploiting Sparsity

For large-scale NLP solving it is typical the Jacobian will be quite sparse. For this reason large-scale NLP solvers (such as IPOPT and it's MI cousin BONMIN) expect the Jacobian to be sparse, as well as the user to supply the sparsity structure (indicating nonzero entries). Exploiting sparsity in this way can substantially improve the efficiency of the solver for large problems.

To supply a sparse Jacobian you can simply wrap your function handle in the `sparse`

function:

jac = @(x) sparse([prod(x')./x';

2.*x']);

However this is not particularly efficient, as a sparse representation of the function will be generated at each function call. A more efficient method is to store it in the same way MATLAB does (sparse triplets), then convert it to the sparse datatype at the end. See the following MathWorks article on efficient handling of sparse matrices.

For the above problem you may note it is not rather sparse (in fact it is 100% dense), however it is just an example. Real examples of large scale nonlinear programming are too big to put here!

To generate the Jacobian structure callback, return a sparse matrix where a 1 indicates a nonzero entry, as follows:

jacstr = @() sparse(ones(2,4));

Note this function requires no input, as the structure is *static*. This means no matter what value of `x`

, the structure should indicate *all possible* nonzero entries. It may be that nonzero entries change depending on `x`

, however you must account for all possible entries in memory.

By default OPTI assumes the Jacobian Structure contains all ones, indicating a fully dense Jacobian. This is the most inefficient, but safest option.

Finally, to supply the entire problem to OPTI:

fun = @(x) x(1)*x(4)*(x(1) + x(2) + x(3)) + x(3);

% Gradient

grad = @(x) [x(1)*x(4) + x(4)*sum(x(1:3)), x(1)*x(4),...

x(1)*x(4) + 1, x(1)*sum(x(1:3))];

% Nonlinear Constraints

nlcon = @(x) [ x(1)*x(2)*x(3)*x(4);

x(1)^2 + x(2)^2 + x(3)^2 + x(4)^2 ];

cl = [25;40];

cu = [Inf;40];

% Jacobian

jac = @(x) sparse([prod(x')./x';

2.*x']);

% Jacobian Structure

jacstr = @() sparse(ones(2,4));

% Bounds (lb <= x <= ub)

lb = ones(4,1);

ub = 5*ones(4,1);

% Initial Guess

x0 = [1 5 5 1]';

% Options

opts = optiset('solver','ipopt','display','iter');

% Build OPTI Problem

Opt = opti('fun',fun,'grad',grad,'nl',nlcon,cl,cu,'jac',jac,'jacstr',jacstr,...

'bounds',lb,ub,'x0',x0,'options',opts)

% Solve NLP

[x,fval,exitflag,info] = solve(Opt)

You will note if you run the above example you will receive *no warnings*! You have now supplied the minimum amount of information required for OPTI to not make any assumptions to solve the problem using IPOPT.

### Example 7: Checking Derivatives

It is all too easy to miss an index, be out by one, or any of the other common programming caveats. This is made even worse with large, sparse, partial derivative matrices! In order to check for errors, OPTI now comes with a derivative checker.

To activate it, enable the `optiset`

option `derivCheck`

:

opts = optiset('solver','ipopt','display','iter','derivCheck','on');

% Build OPTI Problem

Opt = opti('fun',fun,'grad',grad,'nl',nlcon,cl,cu,'jac',jac,'jacstr',jacstr,...

'bounds',lb,ub,'x0',x0,'options',opts)

% Solve NLP

[x,fval,exitflag,info] = solve(Opt)

If you examine the command window print out, you will see it says (for the above problem):

`OPTI Derivative Checker detected no problems in 'Objective Gradient'`

`OPTI Derivative Checker detected no problems in 'Constraint Jacobian'`

Indicating no errors were detected.

### Summary

This section is a very brief introduction into generating first derivatives for optimization. However do not underestimate this topic, *every gradient based optimizer in OPTI requires first derivatives*. Just because OPTI will automatically approximate them using finite-difference does not mean you can ignore them!