Specifying Derivatives

The function (RefLink[FindRoot, paclet:ref/FindRoot]) has a Jacobian option; the functions (RefLink[FindMinimum, paclet:ref/FindMinimum]), (RefLink[FindMaximum, paclet:ref/FindMaximum]), and (RefLink[FindFit, paclet:ref/FindFit]) have a (RefLink[Gradient, paclet:ref/Gradient]) option; and the Newton method has a method option Hessian. All these derivatives are specified with the same basic structure. Here is a summary of ways to specify derivative computation methods.

Methods for computing gradient, Jacobian, and Hessian derivatives.

The basic specification for a derivative is just the method for computing it. However, all of the derivatives take options as well. These can be specified by using a list {method,opts}. Here is a summary of the options for the derivatives.

Options for computing gradient, Jacobian, and Hessian derivatives.

A few examples will help illustrate how these fit together.

This loads a package that contains some utility functions:

<<Optimization`UnconstrainedProblems`

This defines a function that is only intended to evaluate for numerical values of the variables:

f[x_?NumberQ, y_?NumberQ] := Cos[x^2 - 3 y] + Sin[x^2 + y^2]

With just (RefLink[Method, paclet:ref/Method])->"Newton", (RefLink[FindMinimum, paclet:ref/FindMinimum]) issues an lstol message because it was not able to resolve the minimum well enough due to lack of good derivative information.

This shows the steps taken by (RefLink[FindMinimum, paclet:ref/FindMinimum]) when it has to use finite differences to compute the gradient and Hessian:

FindMinimumPlot[f[x,y],{{x,1},{y,1}}, Method->"Newton"]

The following describes how you can use the gradient option to specify the derivative.

This computes the minimum of f[x,y] using a symbolic expression for its gradient:

FindMinimumPlot[f[x,y],{{x,1},{y,1}}, Gradient->{2 x Cos[(x^2)+(y^2)]-2 x Sin[(x^2)-3 y],2 y Cos[(x^2)+(y^2)]+3 Sin[(x^2)-3 y]},Method->"Newton"]

Symbolic derivatives are not always available. If you need extra accuracy from finite differences, you can increase the difference order from the default of 1 at the cost of extra function evaluations.

This computes the minimum of f[x,y] using a second-order finite difference to compute the gradient:

FindMinimumPlot[f[x,y],{{x,1},{y,1}}, Gradient->{Automatic, "DifferenceOrder"->2},Method->"Newton"]

Note that the number of function evaluations is much higher because function evaluations are used to compute the gradient, which is used to approximate the Hessian in turn. (The Hessian is computed with finite differences since no symbolic expression for it can be computed from the information given.)

The information given from FindMinimumPlot about the number of function, gradient, and Hessian evaluations is quite useful. The (RefLink[EvaluationMonitor, paclet:ref/EvaluationMonitor]) options are what make this possible. Here is an example that simply counts the number of each type of evaluation. (The plot is made using (RefLink[Reap, paclet:ref/Reap]) and (RefLink[Sow, paclet:ref/Sow]) to collect the values at which the evaluations are done.)

This computes the minimum with counters to keep track of the number of steps and the number of function, gradient, and Hessian evaluations:

Block[{s = 0, e = 0, g = 0, h = 0},
{FindMinimum[Cos[x^2 - 3 y] + Sin[x^2 + y^2], {{ x,1},{y,1}},StepMonitor:>s++, EvaluationMonitor:>e++, Gradient-> {Automatic, EvaluationMonitor:>g++},  Method->{"Newton", "Hessian"->{Automatic, EvaluationMonitor:>h++}}], s, e, g, h}]

Using such diagnostics can be quite useful for determining what methods and/or method parameters may be most successful for a class of problems with similar characteristics.

When the Wolfram Language can access the symbolic structure of the function, it automatically does a structural analysis of the function and its derivatives and uses (RefLink[SparseArray, paclet:ref/SparseArray]) objects to represent the derivatives when appropriate. Since subsequent numerical linear algebra can then use the sparse structures, this can have a profound effect on the overall efficiency of the search. When the Wolfram Language cannot do a structural analysis, it has to assume, in general, that the structure is dense. However, if you know what the sparse structure of the derivative is, you can specify this with the "Sparse" method option and gain huge efficiency advantages, both in computing derivatives (with finite differences, the number of evaluations can be reduced significantly) and in subsequent linear algebra. This issue is particularly important when working with vector-valued variables. A good example for illustrating this aspect is the extended Rosenbrock problem, which has a very simple sparse structure.

This gets the extended Rosenbrock function with 1000 variables in symbolic form ready to be solved with (RefLink[FindRoot, paclet:ref/FindRoot]) using the UnconstrainedProblems` package:

n = 1000;Short[ pex =GetFindRootProblem[ExtendedRosenbrock, n], 20]

This solves the problem using the symbolic form of the function:

Timing[Norm[1 -( Array[(Subscript[X, #])&, n] /. ProblemSolve[pex])]]

For a function with simple form like this, it is easy to write a vector form of the function, which can be evaluated much more quickly than the symbolic form can, even with automatic compilation.

This defines a vector form of the extended Rosenbrock function, which evaluates very efficiently:

ExtendedRosenbrockResidual[X_List] := Module[{x1, x2},
x1 = Take[X, {1,-1,2}];
x2 = Take[X, {2,-1,2}];
Flatten[Transpose[{10 (x2 - x1^2), 1 - x1}]]]

This extracts the starting point as a vector from the problem structure:

Short[start = pex[[2, All, 2]]]

This solves the problem using a vector variable and the vector function for evaluation:

Timing[Norm[1 - (X /. FindRoot[ExtendedRosenbrockResidual[X], {X, start}])]]

The solution with the function, which is faster to evaluate, winds up being slower overall because the Jacobian has to be computed with finite differences since the x_List pattern makes it opaque to symbolic analysis. It is not so much the finite differences that are slow as the fact that it needs to do 100 function evaluations to get all the columns of the Jacobian. With knowledge of the structure, this can be reduced to two evaluations to get the Jacobian. For this function, the structure of the Jacobian is quite simple.

This defines a pattern (RefLink[SparseArray, paclet:ref/SparseArray]), which has the structure of nonzeros for the Jacobian of the extended Rosenbrock function. (By specifying _ for the values in the rules, the (RefLink[SparseArray, paclet:ref/SparseArray]) is taken to be a template of the (RefLink[Pattern, paclet:ref/Pattern]) type as indicated in the output form.)

sparsity = SparseArray[Flatten[Table[{{i, i}->_, {i, i+1}->_, {i +1, i}->_}, {i, 1, n-1, 2}]]]

This solves the problem with the knowledge of the actual Jacobian structure, showing a significant cost savings:

Timing[Norm[1 - (X /.FindRoot[ExtendedRosenbrockResidual[X], {X, start},Method->{"Newton"}, Jacobian->{Automatic, Sparse->sparsity}])]]

When a sparse structure is given, it is also possible to have the value computed by a symbolic expression that evaluates to the values corresponding to the positions given in the sparse structure template. Note that the values must correspond directly to the positions as ordered in the (RefLink[SparseArray, paclet:ref/SparseArray]) (the ordering can be seen using (RefLink[ArrayRules, paclet:ref/ArrayRules])). One way to get a consistent ordering of indices is to transpose the matrix twice, which results in a (RefLink[SparseArray, paclet:ref/SparseArray]) with indices in lexicographic order.

This transposes the nonzero structure matrix twice to get the indices sorted:

sparsity = Transpose[Transpose[sparsity]]

This defines a function that will return the nonzero values in the Jacobian corresponding to the index positions in the nonzero structure matrix:

ERJValues[X_List] := Module[{x1, zero},
x1 = Take[X, {1,-1,2}];
zero = 0. x1;
Flatten[Transpose[{-20 x1, 10. + zero, -1. + zero}]]]

This solves the problem with the resulting sparse symbolic Jacobian:

Timing[Norm[1 - (X /.FindRoot[ExtendedRosenbrockResidual[X], {X, start},Method->{"Newton"}, Jacobian->{ERJValues[X], Sparse->sparsity}])]]

In this case, using the sparse Jacobian is not significantly faster because the Jacobian is so sparse that a finite difference approximation can be found for it in only two function evaluations and because the problem is well enough defined near the minimum that the extra accuracy in the Jacobian does not make any significant difference.