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

[Spencer Graves]

> On Feb 18, 2015, at 2:15 PM, Spencer Graves <spencer.graves at prodsyse.com> wrote:
> 
> 
>> On Feb 18, 2015, at 1:54 PM, David Winsemius <dwinsemius at comcast.net <mailto:dwinsemius at comcast.net>> wrote:
>> 
>> 
>> 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.
> 
> 
> How much numerical precision do you need?  
> 
> 
> How important is the difference between 4121204 and (4121204-5.494803e-08)?  


p.s.  If that difference is important, your best approach might be to use theory with e2 = 4121204 and d = (-5.494803e-08).  That could give you answers — and insights — more accurate than any infinite precision arithmetic.   
> 
> 
> Spencer 
> 
>> 
>> -- 
>> david.
>>> 
>>> Mike
>>> 
>>> On Wed, Feb 18, 2015 at 11:57 AM, Thierry Onkelinx <thierry.onkelinx at inbo.be <mailto: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 <mailto: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 <mailto: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 <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 <mailto:R-help at r-project.org> mailing list -- To UNSUBSCRIBE and more, see
>>> https://stat.ethz.ch/mailman/listinfo/r-help <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
>> 
>> ______________________________________________
>> R-help at r-project.org <mailto:R-help at r-project.org> mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help <https://stat.ethz.ch/mailman/listinfo/r-help>
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html <http://www.r-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
> 


	[[alternative HTML version deleted]]