Thanks Thierry for the pointer, that's explains the problem. Is there anything I can do about the matrix instability or numerical inaccuracy? Mike On Wed, Feb 18, 2015 at 11:57 AM, Thierry Onkelinx <thierry.onkelinx at inbo.be> wrote:> 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: >> > >> 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 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. >> > >[[alternative HTML version deleted]]
Jeff Newmiller
2015-Feb-18 21:41 UTC
[R] Numerical stability of eigenvalue and hessian matrix in R
This question is getting pretty deep into numerical analysis theory. The usual
approach has already been mentioned... don't expect high accuracy in all
problems. Your specific problem could have a special technique somewhere, but
don't be surprised if we are not experts in your specific problem as such
tricks are not R-specific.
---------------------------------------------------------------------------
Jeff Newmiller The ..... ..... Go Live...
DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live
Go...
Live: OO#.. Dead: OO#.. Playing
Research Engineer (Solar/Batteries O.O#. #.O#. with
/Software/Embedded Controllers) .OO#. .OO#. rocks...1k
---------------------------------------------------------------------------
Sent from my phone. Please excuse my brevity.
On February 18, 2015 1:13:39 PM PST, C W <tmrsg11 at gmail.com>
wrote:>Thanks Thierry for the pointer, that's explains the problem.
>
>Is there anything I can do about the matrix instability or numerical
>inaccuracy?
>
>Mike
>
>On Wed, Feb 18, 2015 at 11:57 AM, Thierry Onkelinx
><thierry.onkelinx at inbo.be
>> wrote:
>
>> 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:
>>> >
>>>
>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 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.
>>>
>>
>>
>
> [[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.
David Winsemius
2015-Feb-18 21:54 UTC
[R] Numerical stability of eigenvalue and hessian matrix in R
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. -- david.> > Mike > > On Wed, Feb 18, 2015 at 11:57 AM, Thierry Onkelinx <thierry.onkelinx at 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 <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: >>>> >>> 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 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. >>> >> >> > > [[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.David Winsemius Alameda, CA, USA
Spencer Graves
2015-Feb-18 22:15 UTC
[R] Numerical stability of eigenvalue and hessian matrix in R
> On Feb 18, 2015, at 1:54 PM, David Winsemius <dwinsemius at 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.onkelinx at 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 <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: >>>>> >>>> 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 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. >>>> >>> >>> >> >> [[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. > > David Winsemius > Alameda, CA, USA > > ______________________________________________ > R-help at r-project.org <mailto:R-help at 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]]