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

[David Winsemius]

On Feb 18, 2015, at 1:13 PM, C W wrote:

> Thanks Thierry for the pointer, that's explains the problem.
> 
> Is there anything I can do about the matrix instability or numerical
> inaccuracy?

There are matrix methods in the Rmpfr package that support increased precision, but it is implemented with S4 methods. You would probably need to retool the numDeriv functions to use the mpfrMatrix-class.

-- 
david.
> 
> 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
>> 
>> 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]]
> 
> ______________________________________________
> 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.

David Winsemius
Alameda, CA, USA