CNLMath : Optimization : jacobian
jacobian
Approximates the Jacobian of m functions in n unknowns using divided differences.
Synopsis
#include <imsl.h>
void imsl_f_jacobian (void fcn(), int m, int n, float y[], float f[], float fjac[], , 0);
The type double function is imsl_d_jacobian.
Required Arguments
void fcn (int indx, float y[], float f[]) (Input/Output)
User-supplied function to compute the value of the function whose Jacobian is to be calculated using divided differences and/or the value of the analytic derivative of that function.
Required Arguments
int indx (Input)
Index of the variable whose derivative is to be computed. imsl_f_jacobian sets this argument to the index of the variable whose derivative is being computed. In those cases where finite differencing and direct analytic computation are combined to calculate a derivative (see optional argument IMSL_ACCUMULATE), imsl_f_jacobian makes an extra call to fcn (with argument indx negative) each time the derivative with respect to variable indx is calculated, in order to calculate the analytic component of that derivative. Note that indx runs from 1 to n, where n is the number of variables.
float y[] (Input)
Array of length n containing the point at which the function is to be computed.
float f[] (Output)
Array of length m, where m is the number of functions to be evaluated at point y, containing the function values at point y. Normally, the user returns the values of the functions evaluated at point y in f. However, when the function can be broken into two parts, a part which is known analytically and a part to be differenced, fcn is called by imsl_f_jacobian once with indx positive for the portion to be differenced and again with indx negative for the portion which is known analytically. In the case where optional argument IMSL_ACCUMULATE has been specified by the user, fcn must be written to handle the known analytic portion separately from the part to be differenced. (See Example 4 for an example where IMSL_ACCUMULATE is used.)
int m (Input)
The number of equations.
int n (Input)
The number of variables.
float y[] (Input)
Array of length n containing the point at which the Jacobian is to be evaluated.
float f[] (Output)
Array of length m containing the function values at point y.
float fjac[] (Input/Output)
m by n array which, on output, contains the estimated Jacobian. Note that if optional argument IMSL_METHOD, method, is used, then for each variable i for which method[i] is set to IMSL_DD_SKIP, array elements fjac[j=0,,m-1][i] are input arguments and must be set to the analytic derivatives with respect to variable i prior to calling imsl_f_jacobian. (See description of optional argument IMSL_METHOD and Example 3 below).
Synopsis with Optional Arguments
#include <imsl.h>
void imsl_f_jacobian (void fcn(), int m, int n, float y[], float f[], float fjac[],
IMSL_YSCALE, float scale[],
IMSL_METHOD, int method[],
IMSL_ACCUMULATE,
IMSL_FACTOR, float factor[],
IMSL_ISTATUS, int istatus[],
IMSL_FCN_W_DATA, void fcn_w_data(),
void *data,
0);
Optional Arguments
IMSL_YSCALE, float scale[] (Input)
An array of length n containing the diagonal scaling matrix for the variables. scale can also be used to provide appropriate signs for the increments. If the user sets scale[0] to 0.0, then differencing increment delj (for variable j = 1, , n) is set to σj * |y[j-1]| * factor[j-1]. Otherwise, delj is set to σj * |scale[j-1]| * factor[j-1]. (See the discussion of optional argument IMSL_FACTOR below for more information about the calculation of delj.)
Default: scale[i=0,,n-1] = 1.0.
IMSL_METHOD, int method[] (Input)
An array of length n containing the methods used to compute the derivatives. method[i] is the method to be used for variable i. method[i] can be one of the values in the following table:
Value
Description
IMSL_DD_ONE_SIDED
Indicates one-sided differences.
IMSL_DD_CENTRAL
Indicates central differences.
IMSL_DD_SKIP
Indicates that the user has set the input Jacobian fjac[j=0,,m-1][i] to the exact analytic derivative of the function with respect to variable i at point y[i], and that the calculation of the divided difference approximation is to be skipped.
See Example 2 and Example 3 below for demonstrations of how this optional argument is used.
Note that if (and only if) IMSL_DD_SKIP is specified for a variable i, the required array elements fjac[j=0,,m-1][i] must be set to the analytic derivatives with respect to variable i prior to calling imsl_f_jacobian. See Example 3 below.
Default: If optional argument IMSL_METHOD is not used, then one-sided differences are used for all variables.
IMSL_ACCUMULATE (Input)
Indicates that divided differences are to be accumulated with a Jacobian value previously initialized by the user with analytically calculated components of the derivatives via a call to fcn using negative values of fcn argument indx. See the description of indx above and Example 4 below.
IMSL_FACTOR, float factor[] (Input)
An array of length n containing the percentage factor for differencing.
For each divided difference for variable j = 1, n, the differencing increment used is delj. (See the Description below for a discussion of the differencing methods.) The value of delj is computed as follows: If scale[0] has been set to 0.0 (see the description of optional argument IMSL_YSCALE above), define ya,j = |y[j-1]| and σj = 1. Otherwise, if scale[j-1] {< , >} 0, define ya,j = |scale[j-1]| and σj = {-1 , 1}. Finally, compute delj = σj ya,j factor[j-1].
By changing the sign of scale[j-1], the difference delj can have any desired orientation, such as staying within bounds on variable j. For central differences, a reduced factor is used for delj that normally results in relative errors as small as ɛ2/3, where ɛ == machine precision = {imsl_f_machine(4), imsl_d_machine(4)} for {single, double} precision. The elements of factor must be such that: ɛ3/4  factor[j-1]  0.1.
Default: All elements of factor are set to ɛ1/2 .
IMSL_ISTATUS, int istatus[] (Output)
Array of length 10 containing status information that might prove useful to a user wanting to gain better control over the differencing parameters. This information can often be ignored. The following table describes the diagnostic information that is returned in each of the entries of array istatus[]:
Index
Description
0
The number of times a function evaluation was computed.
1
The number of columns in which three attempts were made to increase a percentage factor for differencing (i.e., a component in the factor array) but the computed delj (for j = 1, n) remained unacceptably small relative to y[j‑1] or scale[j‑1]. In such cases the percentage factor is set to the square root of machine precision.
2
The number of columns in which the computed delj was zero to machine precision because y[j‑1] or scale[j‑1] was zero. In such cases delj is set to the square root of machine precision.
3
The number of Jacobian columns that had to be recomputed because the largest difference formed in the column was close to zero relative to scale, where
and i denotes the row index of the largest difference in the column currently being processed. index = 9 gives the last column where this occurred.
4
The number of columns whose largest difference is close to zero relative to scale after the column has been recomputed.
5
The number of times scale information was not available for use in the round-off and truncation error tests. This occurs when
where i is the index of the largest difference for the column currently being processed.
6
The number of times the increment for differencing (del) was computed and had to be increased because
(scale[j‑1] + delj- scale[j‑1] was too small relative to y[j‑1]or scale[j‑1].
7
The number of times a component of the factor array was reduced because changes in function values were large and excess truncation error was suspected. index = 8 gives the last column in which this occurred.
8
The index of the last column where the corresponding component of the factor array had to be reduced because excessive truncation error was suspected.
9
The index of the last column where the difference was small and the column had to be recomputed with an adjusted increment (see index = 3). The largest derivative in this column may be inaccurate due to excessive round-off error.
IMSL_FCN_W_DATA, void fcn_w_data(), void *data[] (Input/Output)
User supplied function whose Jacobian is being calculated, and which can also accept a pointer to data that is supplied by the user. data is a pointer to the data to be passed to the user-supplied function. See Example 4 for a demonstration of how this optional argument is used and Passing Data to User-Supplied Functions in the “Introduction” chapter for more details on the use of IMSL_FCN_W_DATA.
Description
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 = { yjj=1,,n } and differenced points y + delj ej , where the ejj=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).
Examples
Example 1
In this example, the Jacobian matrix of
f1(x) = x1x-2
f2(x) = xx1x+ 1
is estimated by the finite-difference method at the point (1.0, 1.0).
 
#include <imsl.h>
#include <stdio.h>
 
void fcn(int, float*, float*);
 
int main()
{
int n = 2, m = 2;
float fjac[4], y[2], f[2];
char *fmt="%14.5e";
 
y[0] = 1.0;
y[1] = 1.0;
 
/* Calculate and print
* Jacobian one-sided difference approximation: */
 
imsl_f_jacobian (fcn, m, n, y, f, fjac, 0);
 
imsl_f_write_matrix ("The Jacobian is:", m, n, fjac,
IMSL_WRITE_FORMAT, fmt, 0);
}
 
void fcn(int indx, float y[], float f[])
{
f[0] = y[0]*y[1] - 2.0;
f[1] = y[0] - y[0]*y[1] + 1.0;
}
Output
 
The Jacobian is:
1 2
1 1.00000e+000 1.00000e+000
2 0.00000e+000 -1.00000e+000
Example 2
A simple use of imsl_f_jacobian is shown. The gradient of the function
f(y1,y2) = a exp(by1) + cy1y22
is required for values
a = 2.5e6, b = 3.4, c = 4.5, y1 = 2.1, y2 = 3.2
The analytic gradient of this function is:
grad(f) = [ab exp(by1) + cy22, 2cy1y2]
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.
 

#include <imsl.h>
#include <stdio.h>
#include <math.h>
 
void fcn(int, float*, float*);
 
void fcn_drv(float*, float*);
 
int main()
{
int n = 2, m = 1, j;
int method[2];
float fjac[2], y[2], f[1], scale[2];
char *fmt="%14.5e";
 
y[0] = 2.1;
y[1] = 3.2;
scale[0] = 1.0;
scale[1] = 8000.0;
 
/* Use one-sided and central differences
* to approximate gradient and print results: */
 
for (j = 0; j <= 1; j++) {
if (j == 0) {
method[0] = IMSL_DD_ONE_SIDED;
method[1] = IMSL_DD_ONE_SIDED;
} else {
method[0] = IMSL_DD_CENTRAL;
method[1] = IMSL_DD_CENTRAL;
}
 
imsl_f_jacobian (fcn, m, n, y, f, fjac,
IMSL_YSCALE, scale,
IMSL_METHOD, method,
0);
 
if (j == 0) {
imsl_f_write_matrix ("One-Sided Jacobian:",
m, n, fjac, IMSL_WRITE_FORMAT, fmt, 0);
} else {
imsl_f_write_matrix ("Central Jacobian:",
m, n, fjac, IMSL_WRITE_FORMAT, fmt, 0);
}
}
 
/* Calculate analytic Jacobian: */
 
fcn_drv (y, fjac);
imsl_f_write_matrix ("Analytic Jacobian:",
m, n, fjac,IMSL_WRITE_FORMAT, fmt, 0);
}
 
void fcn(int indx, float y[], float f[])
{
float a, b, c;
 
a = 2500000.;
b = 3.4;
c = 4.5;
 
f[0] = a * exp (b * y[0]) + c * y[0] * y[1] * y[1];
}
void fcn_drv(float y[], float fjac[])
{
float a, b, c;
 
a = 2500000.;
b = 3.4;
c = 4.5;
 
fjac[0] = a * b * exp (b * y[0]) + c * y[1] * y[1];
fjac[1] = 2 * c * y[0] * y[1];
}
Output
 
One-Sided Jacobian:
1 2
1.07285e+010 9.26819e+001
 
Central Jacobian:
1 2
1.07225e+010 6.17690e+001
 
Analytic Jacobian:
1 2
1.07221e+010 6.04800e+001
Example 3
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.
 
#include <imsl.h>
#include <stdio.h>
#include <math.h>
 
void fcn(int, float*, float*);
 
int main()
{
int n = 2, m = 1;
int method[2];
float fjac[2], y[2], f[1], scale[2];
char *fmt="%14.5e";
 
y[0] = 2.1;
y[1] = 3.2;
scale[0] = 1.0;
scale[1] = 8000.0;
/* One-sided differences for fjac[0]: */
 
method[0] = IMSL_DD_ONE_SIDED;
 
/* Set method[1] to skip differencing for fjac[1]
* and initialize it to analytic derivative: */
 
method[1] = IMSL_DD_SKIP;
fjac[1] = 2.0 * 4.5 * y[0] * y[1];
 
/* Get Gradient approximation: */
 
imsl_f_jacobian (fcn, m, n, y, f, fjac,
IMSL_YSCALE, scale,
IMSL_METHOD, method,
0);
 
/* Print results: */
 
imsl_f_write_matrix ("The Jacobian is:", m, n, fjac,
IMSL_WRITE_FORMAT, fmt, 0);
}
 
void fcn(int indx, float y[], float f[])
{
float a, b, c;
 
a = 2500000.;
b = 3.4;
c = 4.5;
 
f[0] = a * exp (b * y[0]) + c * y[0] * y[1] * y[1];
}
Output
 
The Jacobian is:
1 2
1.07285e+010 6.04800e+001
 
Example 4
This example uses the same data as in Example 2, computing the Jacobian (gradient) of the function:
f(y1, y2) = a exp(by1) + cy1y22
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(y1y2), 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.
 
#include <imsl.h>
#include <stdio.h>
#include <math.h>
 
void fcn_w_data(int, float*, float*, void*);
 
int main()
{
int n = 2, m = 1;
float fjac[2], y[2], f[1], scale[2];
char *fmt="%14.5e";
 
float rdata[3];
 
y[0] = 2.1;
y[1] = 3.2;
scale[0] = 1.0;
scale[1] = 8000.0;
 
/* Set up to pass extra information to the function: */
 
rdata[0] = 2500000.0;
rdata[1] = 3.4;
rdata[2] = 4.5;
 
/* Use optional argument IMSL_ACCUMULATE so that the
* user can specify which function components are to be
* used for divided difference approximation of
* derivatives and which are to be replaced with exact
* analytically calculated derivatives. Both components
* are set to the default one-sided differences method.
*
* Calculate and print Jacobian approximation: */
 
imsl_f_jacobian (NULL, m, n, y, f, fjac,
IMSL_YSCALE, scale,
IMSL_ACCUMULATE,
IMSL_FCN_W_DATA, fcn_w_data, (void*) rdata,
0);
 
imsl_f_write_matrix ("The Jacobian is:", m, n, fjac,
IMSL_WRITE_FORMAT, fmt, 0);
}
 
void fcn_w_data(int indx, float y[], float f[],
void* data)
{
float a, b, c;
float *rdata = (float*) data;
 
a = rdata[0];
b = rdata[1];
c = rdata[2];
 
/* Handle both the differenced part and the part that is
* known analytically for each dependent variable: */
 
switch (indx) {
case 1:
f[0] = a * exp(b *y[0]);
break;
case -1:
f[0] = c * y[1] *y[1];
break;
case 2:
f[0] = c * y[0] * y[1] * y[1];
break;
case -2:
f[0] = 0.;
break;
}
}
Output
The Jacobian is:
1 2
1.07285e+010 6.04862e+001
Example 5
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.
 
#include <imsl.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
 
#define M 2
#define N 2
#define LDFJAC 2
 
void rosbck(int m, int n, float *x, float *f);
void jacobian(int m, int n, float *x, float *fjac, int fjac_col_dim);
void fcn(int indx, float *x, float *f, void* data);
 
void main()
{
int ibtype = 0;
float xlb[N] = {-2.0, -1.0};
float xub[N] = {0.5, 2.0};
float xguess[N] = {-1.2, 1.0};
float *fvec;
float *x;
 
x = imsl_f_bounded_least_squares (rosbck, M, N, ibtype, xlb, xub,
IMSL_JACOBIAN, jacobian,
IMSL_XGUESS, xguess,
IMSL_FVEC, &fvec,
0);
 
printf("x[0] = %f\n", x[0]);
printf("x[1] = %f\n\n", x[1]);
printf("fvec[0] = %f\n", fvec[0]);
printf("fvec[1] = %f\n\n", fvec[1]);
}
 
void rosbck (int m, int n, float *x, float *f)
{
f[0] = 10.0*(x[1] - x[0]*x[0]);
f[1] = 1.0 - x[0];
}
 
void jacobian (int m, int n, float *x, float *fjac, int fjac_col_dim)
{
int idata[2];
float *f = NULL;
f = (float*)malloc(m*sizeof(float));
idata[0] = m;
idata[1] = n;
imsl_f_jacobian(NULL, m, n, x, f, fjac,
IMSL_FCN_W_DATA, fcn, (void*) idata,
0);
if (f) free(f);
}
 
void fcn (int indx, float *x, float *f, void* data)
{
int *idata = (int*) data;
int m = idata[0];
int n = idata[1];
 
rosbck (m, n, x, f);
}
Output
 
x[0] = 0.500000
x[1] = 0.250000
 
fvec[0] = 0.000000
fvec[1] = 0.500000