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]]
Thierry Onkelinx
2015-Feb-18 16:57 UTC
[R] 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: > > > 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]]
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]]