Approximates the Jacobian of m functions in n unknowns using divided differences.

The type double function is imsl_d_jacobian.

The function imsl_f_jacobian computes the Jacobian matrix for a function f(y) with m components and n independent variables. imsl_f_jacobian uses divided finite differences to compute the Jacobian. This function is designed for use in numerical methods for solving nonlinear problems where a Jacobian is evaluated repeatedly at neighboring arguments. For example, this occurs in a Gauss-Newton method for solving non-linear least squares problems, or in a non-linear optimization method.

imsl_f_jacobian is suited for applications where the Jacobian is a dense matrix. All cases m < n, m = n, or m > n are allowed. Both one-sided and central divided differences can be used.

The design allows for computation of derivatives in a variety of contexts. Note that a gradient should be considered as the special case with m = 1, n ≥ 1. A derivative of a single function of one variable is the case m = 1, n = 1. Any non-linear solving routine that optionally requests a Jacobian or gradient can use imsl_f_jacobian. This should be considered if there are special properties or scaling issues associated with f(y). Use the optional argument IMSL_METHOD to specify different differencing options for numerical differentiation. These can be combined with some analytic subexpressions or other known relationships.

Two divided differences methods are implemented in imsl_f_jacobian for computing the Jacobian: one-sided and central differences.

One-sided differences are computed:

using values of the independent variables at the Jacobian evaluation point y = { yj, j=1,…,n } and differenced points y + delj ej , where the ej , j=1,…,n are the unit coordinate vectors.

Central differences are computed:

The value for each difference delj depends on the variable yj , the differencing method, and the scaling for that variable. This difference is computed internally. See IMSL_FACTOR for computational details. fi(y) is evaluated with user-supplied argument fcn, where index j, variable y, and output f == fi(y) are arguments to fcn.

There are five examples provided that illustrate various ways to use imsl_f_jacobian. For a discussion of the expected errors for the difference methods, see Ralston (1965).

Function imsl_f_jacobian is based upon the Fortran 77 program SJACG, which was designed and programmed by D. A. Salane, Sandia Labs (1986) and modified by R. J. Hanson, Rice University (June, 2002) with advice from F. T. Krogh. See Salane (1986).

In this example, the Jacobian matrix of

is estimated by the finite-difference method at the point (1.0, 1.0).

A simple use of imsl_f_jacobian is shown. The gradient of the function

is required for values

The analytic gradient of this function is:

Note that the comparison of the Jacobian estimates using one-sided and central differences with the exact analytic Jacobian results given in this example demonstrates the increased accuracy afforded by use of central differences. However, these estimates require up to twice as many function calculations as do the one-sided differences estimates for Jacobians with a large number of variables.

This example uses the same data as in Example 2. Here we assume that the second component of the gradient is analytically known. Therefore only the first gradient component needs numerical approximation. The input value IMSL_DD_SKIP of array element method[1] specifies that numerical differentiation with respect to y2 is skipped.

This example uses the same data as in Example 2, computing the Jacobian (gradient) of the function:

For this example, the analytic derivative of the second term with respect to y1, cy22, is provided by the user (in the example, see case -1: within the user-provided function fcn). This leaves only the first term, a exp(by1), to be evaluated in order to use direct differencing to calculate the first partial (see case 1: within the user-provided function fcn). Also, since the first term does not depend on the second variable, y2, it can be left out of the function evaluation when computing the partial with respect to y2 using differencing methods, potentially avoiding cancellation errors (see case 2: within the user-provided function fcn). Since the code does not specify the analytic derivative with respect to y2 for either of the two terms of f(y1, y2), fcn returns f[0] set to 0 for case -2:. The use of optional argument IMSL_ACCUMULATE thereby allows imsl_f_jacobian to use these facts to obtain greater accuracy using a minimum number of computations of the exponential function.

The Jacobian is:

1 2

This example uses CNL function imsl_f_bounded_least_squares to solve the nonlinear least-squares problem:

where: f0(x) = 10(x1-x02), f1(x) = 1-x0, an initial guess (-1.2,  1.0) is supplied, and the residual at the approximate solution is returned. This example is identical to Example 2 of imsl_f_bounded_least_squares, except that Example 2 uses an analytic Jacobian, and this example uses imsl_f_jacobian to approximate the Jacobian using the default one-sided differences.

Note that the function vector whose sum of squares is to be minimized, rosbck, is supplied directly (as a required argument) to imsl_f_bounded_least_squares and indirectly to imsl_f_jacobian, wrapped in function fcn. Function fcn is supplied to imsl_f_jacobian via optional argument IMSL_FCN_W_DATA, fcn, (void*) idata. imsl_f_jacobian is called from within function jacobian which is passed to imsl_f_bounded_least_squares via optional argument IMSL_JACOBIAN, jacobian.

Also note that the array size parameters m and n are passed to function rosbck (which is wrapped in function fcn for use by imsl_f_jacobian) via integer array idata, which is specified in the optional argument IMSL_FCN_W_DATA, fcn, (void*) idata. This is an example of how to pass necessary integer data to imsl_f_jacobian required argument function fcn using IMSL_FCN_W_DATA; an example of passing real data using IMSL_FCN_W_DATA is given in Example 4 above.

Example 5 demonstrates how imsl_f_jacobian can be used to supply estimates of the Jacobian matrix that are necessary for solving many optimization problems when the function to be minimized is complex and its Jacobian cannot be calculated analytically.