[Rd] NaN and linear algebra

Bill Northcott w.northcott at unsw.edu.au
Tue Mar 22 04:31:31 CET 2005

On 21/03/2005, at 10:09 PM, David Firth wrote:
> I am sorry that I wasn't clear.  All that I meant was that *this*
> problem can result in different behaviour in "ordinary" statistical
> applications.  For example, if the objective function in a call to
> optim() involves calling one of these linear algebra routines, the
> result may be NaN (on systems other than Mac OS X) --- which optim will
> typically handle sensibly --- or something else (an error, or perhaps
> some consequence of getting 0 for the determinant) under Mac OS X .
> Probably this was obvious to you.  Apologies if I misled you into
> thinking that there was some other problem I knew about.
> Best regards,
> David
>>> At 11:57 AM +0000 3/16/05, David Firth wrote:
>>>> I don't know whether this is a bug, or a problem with the way I
>>>> built R 2.0.1 (under Mac OS 10.3 on a G5), or something else.  Can
>>>> anyone else confirm (or otherwise) that this happens in their R
>>>> 2.0.1 on Mac OS X?
>>>>>  d<-matrix(NaN,3,3)
>>>>>   d
>>>>     [,1] [,2] [,3]
>>>> [1,]  NaN  NaN  NaN
>>>> [2,]  NaN  NaN  NaN
>>>> [3,]  NaN  NaN  NaN
>>>>>   solve(d)
>>>> Error in solve.default(d) : Lapack routine dgesv: system is exactly
>>>> singular
>>>>>   chol(d)
>>>> Error in chol(d) : the leading minor of order 1 is not positive
>>>> definite
>>>>>   det(d)
>>>> [1] 0
>>>> Doing the same thing on a Windows setup gave a different (and more
>>>> useful, I think) result

As I see it, the MacOS X behaviour is not IEEE-754 compliant.

I had a quick look at the IEEE web site and it seems quite clear that 
NaNs should not cause errors, but propagate through calculations to be 
tested at some appropriate (not too frequent) point.

It seems to me the most useful response is to file a bug on Apple's bug 
reporter complete with simple test case etc.  It might also help to 
bitch about it on Apple's scitech or hpc mail lists, which get 
monitored by Apple people.

It certainly is not up to R developers to fix it.

Bill Northcott

More information about the R-devel mailing list