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

[Ben Bolker]

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