Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

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>