Numerical stability of eigenvalue and hessian matrix in R
Have a look at FAQ 7.31 ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance Kliniekstraat 25 1070 Anderlecht Belgium To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey 2015-02-18 17:44 GMT+01:00 C W <tmrsg11 at gmail.com>:
Hi Ben and JS, Thanks for the reply. I tried using: hessian(func = h_x, x, method = "complex"), it gives zero, that's good. # R code
hess.h <- hessian(func = h_x, x, method = "complex")
mat <- h_x(x)*hess.h - grad(h_x, x) %o% grad(h_x, x)
mat
[,1] [,2] [,3] [,4]
[1,] 2060602 0 0 0
[2,] 0 2060602 0 0
[3,] 0 0 -4039596 -816080
[4,] 0 0 -816080 4039596
But later I do,
eigen(mat)
$values
[1] -4121204 4121204 2060602 2060602
$vectors
[,1] [,2] [,3] [,4]
[1,] 0.00000000 0.00000000 1 0
[2,] 0.00000000 0.00000000 0 1
[3,] -0.99503719 0.09950372 0 0
[4,] -0.09950372 -0.99503719 0 0
And here is the problem,
eigen(mat)$values[2] == 4121204
[1] FALSE
eigen(mat)$values[2] - 4121204
[1] -5.494803e-08 Why doesn't the second value equal to 412104? How do I overcome that? Thanks, Mike On Wed, Feb 18, 2015 at 9:13 AM, JS Huang <js.huang at protective.com> wrote:
Hi, Since all entries in your hessian matrix and grad vector are integers,
I
suggest you execute the following for mat assignment.
mat <- round(h_x(x),digits=0)*round(hess.h,digits=0) - round(grad(h_x, x),digits=0) %o% round(grad(h_x, x),digits=0) mat
[,1] [,2] [,3] [,4] [1,] 0 0 0 -4080400 [2,] 0 7920000 0 -1600000 [3,] 0 0 12160400 0 [4,] -4080400 -1600000 0 -7920000 -- View this message in context:
Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
______________________________________________ R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.