> On Feb 18, 2015, at 1:54 PM, David Winsemius <dwinsem...@comcast.net> wrote: > > > On Feb 18, 2015, at 1:13 PM, C W wrote: > >> Thanks Thierry for the pointer, that's explains the problem. >> >> Is there anything I can do about the matrix instability or numerical >> inaccuracy? > > There are matrix methods in the Rmpfr package that support increased > precision, but it is implemented with S4 methods. You would probably need to > retool the numDeriv functions to use the mpfrMatrix-class.
How much numerical precision do you need? How important is the difference between 4121204 and (4121204-5.494803e-08)? Spencer > > -- > david. >> >> Mike >> >> On Wed, Feb 18, 2015 at 11:57 AM, Thierry Onkelinx <thierry.onkel...@inbo.be >>> wrote: >> >>> Have a look at FAQ 7.31 >>> >>> ir. Thierry Onkelinx >>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and >>> >>> 2015-02-18 17:44 GMT+01:00 C W <tmrs...@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.hu...@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: >>>>> >>>> http://r.789695.n4.nabble.com/Numerical-stability-of-eigenvalue-and-hessian-matrix-in-R-tp4703443p4703456.html >>>>> Sent from the R help mailing list archive at Nabble.com. >>>>> >>>>> ______________________________________________ >>>>> R-help@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@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@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. > > David Winsemius > Alameda, CA, USA > > ______________________________________________ > R-help@r-project.org <mailto:R-help@r-project.org> mailing list -- To > UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > <https://stat.ethz.ch/mailman/listinfo/r-help> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > <http://www.r-project.org/posting-guide.html> > and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] ______________________________________________ R-help@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.
