> On Feb 18, 2015, at 2:15 PM, Spencer Graves <spencer.gra...@prodsyse.com> > wrote: > > >> On Feb 18, 2015, at 1:54 PM, David Winsemius <dwinsem...@comcast.net >> <mailto: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)?
p.s. If that difference is important, your best approach might be to use theory with e2 = 4121204 and d = (-5.494803e-08). That could give you answers — and insights — more accurate than any infinite precision arithmetic. > > > Spencer > >> >> -- >> david. >>> >>> Mike >>> >>> On Wed, Feb 18, 2015 at 11:57 AM, Thierry Onkelinx >>> <thierry.onkel...@inbo.be <mailto: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 >>>> <mailto: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 >>>>> <mailto: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 >>>>> >>>>> <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 <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 >>> 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.
