Man Scilab

derivative

derivative - approximate derivatives of a function

Calling Sequence

derivative(F,x)
[J [,H]] = derivative(F,x [,h ,order ,H_form ,Q])

Parameters

Description

Numerical approximation of the first and second derivatives of a function F: R^n --> R^m at the point x. The Jacobian is computed by approximating the directional derivatives of the components of F in the direction of the columns of Q. (For m=1, v=Q(:,k) : grad(F(x))*v = Dv(F(x)).) The second derivatives are computed by composition of first order derivatives. If H is given in its default form the Taylor series of F(x) up to terms of second order is given by :


  F(x+dx) = F(x) + J(x)*dx + 1/2*H(x)*(dx .*. dx) + ...
   
    

(([J,H]=derivative(F,x,H_form='default'), J=J(x), H=H(x).)

Remarks

Numerical approximation of derivatives is generally an unstable process. The step size h must be small to get a low error but if it is too small floating point errors will dominate by cancellation. As a rule of thumb don't change the default step size. To work around numerical difficulties one may also change the order and/or choose different orthogonal matrices Q (the default is eye(n,n)), especially if the approximate derivatives are used in optimization routines. All the optional arguments may also be passed as named arguments, so that one can use calls in the form :


derivative(F, x, H_form = "hypermat")
derivative(F, x, order = 4) etc.

   
    

Examples


  function y=F(x)
   y=[sin(x(1)*x(2))+exp(x(2)*x(3)+x(1)) ; sum(x.^3)];
 endfunction
 function y=G(x,p) 
   y=[sin(x(1)*x(2)*p)+exp(x(2)*x(3)+x(1)) ; sum(x.^3)];
 endfunction

 x=[1;2;3];[J,H]=derivative(F,x,H_form='blockmat');
 disp(J)
 disp(H)

 n=3;
 // form an orthogonal matrix :   
 nu=0; while nu<n, [Q,nu]=colcomp(rand(n,n)); end  
 for i=[1,2,4]
    [J,H]=derivative(F,x,order=i,H_form='blockmat',Q=Q);
    mprintf("order= %d \n",i);
    disp(H);
 end

 p=1;h=1e-3;
 [J,H]=derivative(list(G,p),x,h,2,H_form='hypermat');
 disp(H);
 [J,H]=derivative(list(G,p),x,h,4,Q=Q);
 disp(H)

 // Taylor series example:
 dx=1e-3*[1;1;-1];
 [J,H]=derivative(F,x);
 F(x+dx)
 F(x+dx)-F(x)
 F(x+dx)-F(x)-J*dx
 F(x+dx)-F(x)-J*dx-1/2*H*(dx .*. dx)

 // A trivial example
 function y=f(x,A,p,w), y=x'*A*x+p'*x+w; endfunction
 // with Jacobian and Hessean given by J(x)=x'*(A+A')+p', and H(x)=A+A'.
 A = rand(3,3); p = rand(3,1); w = 1;
 x = rand(3,1);
 [J,H]=derivative(list(f,A,p,w),x,h=1,H_form='blockmat')
 // Since f(x) is quadratic in x, approximate derivatives of order=2 or 4 by finite
 // differences should be exact for all h~=0. The apparent errors are caused by
 // cancellation in the floating point operations, so a "big" h is choosen.
 // Comparison with the exact matrices:
 Je = x'*(A+A')+p'
 He = A+A'
 clean(Je - J)
 clean(He - H)
   
  

See Also

numdiff ,   derivat ,  

Author

Rainer von Seggern, Bruno Pincon

Back