[R] Numerical stability of eigenvalue and hessian matrix in R

[C W]

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]]