Language Reference

NLPFDD Call

CALL NLPFDD (f, g, h, "fun", x0, <*>, par <*>, "grd" );

The NLPFDD subroutine uses the finite-differences method to approximate derivatives.

See the section Nonlinear Optimization and Related Subroutines for a listing of all NLP subroutines. See Chapter 14 for a description of the arguments of NLP subroutines.

The NLPFDD subroutine can be used for the following tasks:

  • If the module "fun" returns a scalar, the NLPFDD subroutine computes the function value f, the gradient vector g, and the Hessian matrix h, all evaluated at the point x0.

  • If the module "fun" returns a column vector of m function values, the subroutine assumes that a least squares function is specified, and it computes the function vector f, the Jacobian matrix $\mb{J}$, and the crossproduct of the Jacobian matrix $\mb{J'J}$ at the point x0. In this case, you must set the first element of the par argument to m.

If any of the results cannot be computed, the subroutine returns a missing value for that result.

You can specify the following input arguments with the NLPFDD subroutine:

  • The "fun" argument refers to a module that returns either a scalar value or a column vector of length m. This module returns the value of the objective function or, for least squares problems, the values of the m functions that the objective function comprises.

  • The x0 argument is a vector of length n that defines the point at which the functions and derivatives should be computed.

  • The par argument is a vector that defines options and control parameters. The par argument in the NLPFDD call is different from the one used in the optimization subroutines.

  • The "grd" argument is optional and refers to a module that returns a vector that defines the gradient of the function at x0. If the "fun" argument returns a vector of values instead of a scalar, the "grd" argument is ignored.

If the "fun" module returns a scalar, the subroutine returns the following values:

  • f is the value of the function at the point x0.

  • g is a vector that contains the value of the gradient at the point x0. If you specify the "grd" argument, the gradient is computed from that module. Otherwise, the approximate gradient is computed by a finite difference approximation by using calls of the function module in a neighborhood of x0.

  • h is a matrix that contains a finite difference approximation of the value of the Hessian at the point x0. If you specify the "grd" argument, the Hessian is computed by calls of that module in a neighborhood of x0. Otherwise, it is computed by calls of the function module in a neighborhood of x0.

If the "fun" module returns a vector, the subroutine returns the following values:

  • f is a vector that contains the values of the m functions that comprise objective function at the point x0.

  • g is the $m\times n$ Jacobian matrix $\mb{J}$, which contains the first-order derivatives of the functions with respect to the parameters, evaluated at x0. It is computed by finite difference approximations in a neighborhood of x0.

  • h is the $n\times n$ crossproduct of the Jacobian matrix, $\mb{J}^ T\mb{J}$. It is computed by finite difference approximations in a neighborhood of x0.

The par argument is a vector of length 3.

  • par[1] corresponds to the opt[1] argument in the optimization subroutines. This argument is relevant only to least squares optimization methods, in which case it specifies the number of functions returned by the module "fun". If par[1] is missing or is smaller than 1, it is set to 1.

  • par[2] corresponds to the opt[8] argument in the optimization subroutines. It determines what type of approximation is to be used and how the finite difference interval, h, is to be computed. See the section Finite-Difference Approximations of Derivatives for details.

  • par[3] corresponds to the par[8] argument in the optimization subroutines. It specifies the number of accurate digits in evaluating the objective function. The default is $-\log _{10}(\epsilon )$, where $\epsilon $ is the machine precision.

If you specify a missing value in the par argument, the default value is used.

The NLPFDD subroutine is particularly useful for checking your analytical derivative specifications of the "grd", "hes", and "jac" modules. You can compare the results of the modules with the finite difference approximations of the derivatives of f at the point x0 to verify your specifications.

In the unconstrained Rosenbrock problem (see the section Unconstrained Rosenbrock Function), the objective function is

\[ f(x) = 50(x_2-x_1^2)^2 + \frac{1}{2}(1-x_1)^2 \]

The gradient and the Hessian are

\begin{eqnarray*} g^{\prime }(x,y) & = & \left[ \begin{array}{c} \frac{\partial f}{\partial x_1} \\[0.1in] \frac{\partial f}{\partial x_2} \end{array} \right] = \left[ \begin{array}{c} 200x_1^3 - 200x_1x_2 + x_1 - 1 \\ -100x_1^2 + 100x_2 \end{array} \right]\\[0.10in] \mb{H}(x,y) & = & \left[ \begin{array}{cc} \frac{\partial ^2 f}{\partial x_1^2} & \frac{\partial ^2 f}{\partial x_1 \partial x_2} \\[0.10in] \frac{\partial ^2 f}{\partial x_2 \partial x_1} & \frac{\partial ^2 f}{\partial x_2^2} \end{array} \right] = \left[ \begin{array}{cc} 600x_1^2 - 200x_2 + 1 & -200x_1 \\ -200x_1 & 100 \end{array} \right] \end{eqnarray*}

At the point $x=(2,7)$, these matrices evaluate to

\begin{eqnarray*} g^{\prime }(2,7) & =& \left[ \begin{array}{r} -1199 \\ 300 \end{array} \right] \\[0.10in] \mb{H}(2,7) & = & \left[ \begin{array}{rr} 1001 & -400 \\ -400 & 100 \end{array} \right] \end{eqnarray*}

The following statements define the Rosenbrock function and use the NLPFDD call to compute the gradient and the Hessian: 

start F_ROSEN(x);
   y1 = 10 * (x[2] - x[1] * x[1]);
   y2 =  1 - x[1];
   f  = 0.5 * (y1 * y1 + y2 * y2);
   return(f);
finish F_ROSEN;
x = {2 7};
call nlpfdd(crit, grad, hess, "F_ROSEN", x);
print grad;
print hess;

Figure 24.234: Gradient and Hessian at a Point

grad
-1199 300.00001

hess
1000.9998 -400.0018
-400.0018 99.999993 




            

            


            
If the Rosenbrock problem is considered from a least squares perspective, the two functions are 

            

               
\begin{eqnarray*}  f_1(x) &  = &  10(x_2-x_1^2) \\ f_2(x) &  = &  1-x_1 \end{eqnarray*}

            

            
 The Jacobian and the crossproduct of the Jacobian are 

            

               
\begin{eqnarray*}  \mb{J} &  = &  \left[ \begin{array}{cc} \frac{\partial f_1}{\partial x_1} &  \frac{\partial f_1}{\partial x_2} \\[0.10in] \frac{\partial f_2}{\partial x_1} &  \frac{\partial f_2}{\partial x_2} \\ \end{array} \right] = \left[ \begin{array}{cc} -20x_1 &  10 \\ -1 &  0 \end{array} \right] \\[0.10in] \mb{J}^ T\mb{J} &  = &  \left[ \begin{array}{cc} 400x_1^2+1 &  -200x_1 \\ -200x_1 &  100 \end{array} \right] \end{eqnarray*}

            

            
 At the point $x=(2,7)$, these matrices evaluate to 
            

            

               
\begin{eqnarray*}  \mb{J}(2,7) &  = &  \left[ \begin{array}{rr} -40 &  10 \\ -1 &  0 \end{array} \right] \\[0.10in] \mb{J}^ T\mb{J}|_{(2,7)} &  = &  \left[ \begin{array}{rr} 1601 &  -400 \\ -400 &  100 \end{array} \right] \end{eqnarray*}

            

            
The following statements define the Rosenbrock problem in a least squares framework and use the NLPFDD call to compute the
               Jacobian and the crossproduct matrix. Since the value of the PARMS variable, which is used for the par argument, is 2, the NLPFDD subroutine allocates memory for a least squares problem with two functions, $f_1(x)$ and $f_2(x)$. 
            


start F_ROSEN(x);
   y = j(2, 1, 0);
   y[1] = 10 * (x[2] - x[1] * x[1]);
   y[2] =  1 - x[1];
   return(y);
finish F_ROSEN;
x     = {2 7};
parms = 2;
call nlpfdd(fun, jac, crpj, "F_ROSEN", x, parms);
print jac;
print crpj;



  

            


Figure 24.235: Jacobian and Crossproduct Matrix at a Point

               

                  

                     
                     

                        
                        
                        

                           
                           

                              
                              







jac
                                       
                                    




-40
                                       
                                       10
                                       
                                    



-1
                                       
                                       0
                                       
                                    






                           
                        

                        



                     
                     
                  

                  

                     
                     

                        
                        
                        

                           
                           

                              
                              







crpj
                                       
                                    




1601
                                       
                                       -400
                                       
                                    



-400
                                       
                                       100