Re: Numerically computing partial derivatives

*To*: mathgroup at smc.vnet.net*Subject*: [mg47965] Re: [mg47940] Numerically computing partial derivatives*From*: Selwyn Hollis <sh2.7183 at misspelled.erthlink.net>*Date*: Tue, 4 May 2004 01:08:46 -0400 (EDT)*References*: <200405020850.EAA00333@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Hi Mark, Here's a numerical Hessian program that's based on a numerical gradient I wrote some time ago. I doubt that it's as efficient as it could be, but it might suit your needs. Options[NGrad] = (Options[NHessian] = {StepSize -> .005}); varspattern := {_Symbol..}|{Subscript[_Symbol,_]..}|{_Symbol[_Integer]..}; NGrad[f_, x:varspattern, opts___?OptionQ][p:{(_)..}]/; Length[x]==Length[p] := Module[{h=N[StepSize/.{opts}/.Options[NGrad]],n=Length[p],pp,hI}, pp = PadRight[{}, n, {p}]; hI = h*IdentityMatrix[n]; ((f/.Thread[x->#]&)/@(pp+hI)-(f/.Thread[x->#1]&)/@(pp-hI))/(2*h)] NHessian[f_, x:varspattern, opts___?OptionQ][p:{(_)..}]/; Length[x]==Length[p] := NGrad[NGrad[f, x, opts][x], x, opts][p] Example: vars = Array[x, 5]; fn = Plus @@ (vars^2)*(x[1]-x[2]) pt = Table[Random[], {5}]; NHessian[fn, vars][pt]// MatrixForm and for comparison, (D[D[fn,#]& /@ vars,#] &/@ vars /.Thread[vars->pt])//MatrixForm Note: Beware of using very small stepsizes! ----- Selwyn Hollis http://www.math.armstrong.edu/faculty/hollis (edit reply-to to reply) On May 2, 2004, at 4:50 AM, Mark Coleman wrote: > Greetings, > > I am working with a maximum likelihood problem in econometrics. With > some effort I can get Mathematica v5.0 to maximize the function. In > order to > derive standard errors of the estimates, however, I need to calculate > the Hessian of the function at the optimal solution. This requires, of > course, calculating the set of second derivatives of the function. Due > to the nature of the function, however, neither the built-in D or ND > operator seem to work (note: The function contains a term of the form > Log[Det[I-rho*W]], where I is the nxn Identity matrix, rho is a real, > and W is a non-symmetric (sparse) matrix of reals, or order n. In > practice, n > 1000 at times. As a result, symbolic differentiation is > not feasible. In addition, when I use ND [], I get nonsensical answers. > > As a result, I'd like to use a simple numerical differentiation scheme > based upon a Taylor-expansion, e.g., the simplest being df(x)/dx = > (f(x+h)-f(x)) / h. I thought that before I code or two or three-point > algorithm I would check on the list and see if anyone had already > written code for a numerical Hessian. > > Thanks, > > Mark > >

**References**:**Numerically computing partial derivatives***From:*Mark Coleman <mark@markscoleman.com>