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

C W tmrsg11 at gmail.com
Wed Feb 18 17:44:57 CET 2015


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



More information about the R-help mailing list