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>
- Numerically computing partial derivatives