Hi list,
I am running the following R code, the answer should be zero. But R gives
a very small negative number, what should I do?
##R code
library(numDeriv)
h_x <- function(x){
a = x[1]
b = x[2]
c = x[3]
d = x[4]
(a^2 + c^2 + d^2) * (b^2 + c^2 + d^2)
}
x1 = 10
x2 = 1
x3 = 0
x4 = 10
x = c(x1, x2, x3, x4)
hess.h <- hessian(func = h_x, x)
mat <- h_x(x)*hess.h - grad(h_x, x) %o% grad(h_x, x)
> mat
[,1] [,2] [,3] [,4]
[1,] 4.059961e-05 -3.038031e-04 -1.307195e-02 -4.080400e+06
[2,] -3.038031e-04 7.920000e+06 -2.663690e-02 -1.600000e+06
[3,] -1.307195e-02 -2.663690e-02 1.216065e+07 1.331894e-02
[4,] -4.080400e+06 -1.600000e+06 1.331894e-02 -7.920000e+06
# A lot of the elements should be zero, for example the first one. I will
change them manually.
m = ifelse(abs(mat)< 0.1, 0, mat)
eigen(m)$val
[1] 12160648 8103268 1665782 -9769050
> sum(eigen(m)$val[2:4])
[1] -0.0005430793
## End of R code
The last answer above should be zero, but it's appearing as a very small
negative number. How should I deal with it?
Thanks lots,
Mike
[[alternative HTML version deleted]]
JS Huang
2015-Feb-18 14:13 UTC
[R] Numerical stability of eigenvalue and hessian matrix in R
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.
Ben Bolker
2015-Feb-18 14:36 UTC
[R] Numerical stability of eigenvalue and hessian matrix in R
C W <tmrsg11 <at> gmail.com> writes:> > Hi list, > > I am running the following R code, the answer should be zero. But R gives > a very small negative number, what should I do? > > ##R code >library(numDeriv) h_x <- function(x){ a = x[1] b = x[2] c = x[3] d = x[4] (a^2 + c^2 + d^2) * (b^2 + c^2 + d^2) } x1 = 10 x2 = 1 x3 = 0 x4 = 10 x = c(x1, x2, x3, x4) hess.h <- hessian(func = h_x, x) mat <- h_x(x)*hess.h - grad(h_x, x) %o% grad(h_x, x) m <- zapsmall(mat) sum(eigen(m)$val[2:4]) I get a much smaller answer than you do (several orders of magnitude): [1] 2.421439e-08 This is with: R Under development (unstable) (2015-02-11 r67792) Platform: i686-pc-linux-gnu (32-bit) Running under: Ubuntu precise (12.04.5 LTS) numDeriv_2012.9-1 However, you mostly have to learn to deal with the inherent imprecision of floating-point calculations. It's hard to tell without any further context, but there may be a more numerically stable way to do what you want.> [1] 12160648 8103268 1665782 -9769050 > > > sum(eigen(m)$val[2:4]) > > [1] -0.0005430793 > > ## End of R code > > The last answer above should be zero, but it's appearing as a very small > negative number. How should I deal with it? > > Thanks lots, > > Mike > > [[alternative HTML version deleted]] > >
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]]