[Rd] NaN and linear algebra

[Bill Northcott]

On 23/03/2005, at 12:55 AM, Simon Urbanek wrote:
>> 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.
>
> This is not quite correct and in fact irrelevant to the problem you 
> describe. NaNs may or may not signal, depending on how they are used. 
> Certain operations on NaN must signal by the IEEE-754 standard. The 
> error you get is not a trap, and it's not a result of a signal check, 
> either. The whole problem is that depending on which algorithm is 
> used, the NaNs will be used different ways and thus may or may not use 
> signaling operations.

It may not violate the letter of IEEE-754 because matrix calculations 
are not covered, but it certainly violates the spirit that arithmetic 
should be robust and programs should not halt on these sorts of errors.
>
> I don't consider the `solve' error a bug - in fact I would rather get 
> an error telling me that something is wrong (although I agree that the 
> error is misleading - the error given in Linux is a bit more helpful) 
> than getting a wrong result.

You may prefer the error, but it is not in the sprit of robust 
arithmetic. ie
 > d<-matrix(NaN,3,3)
 > f<-solve(d)
Error in solve.default(d) : Lapack routine dgesv: system is exactly 
singular
 > f
Error: Object "f" not found

> If I would mark something in your example as a bug that would be 
> det(m)=0, because it should return NaN (remember, NaN==NaN is FALSE; 
> furthermore if det was calculated inefficiently using Laplace 
> expansion, the result would be NaN according to IEEE rules). det=0 is 
> consistent with the error given, though. Should we check this in R 
> before calling Lapack - if the vector contains NaNs, det/determinant 
> should return NaN right away?

Clearly det(d) returning 0 is wrong.  As a result based on a 
computation including a NaN, it should return NaN.  The spirit of 
IEEE-754 is that the programmer should choose the appropriate point at 
which to check for NaNs.  I would interpret this to mean the R 
programmer not the R library developer.  Surely that is why R provides 
the is.nan function.

 > d
      [,1] [,2] [,3]
[1,]  NaN  NaN  NaN
[2,]  NaN  NaN  NaN
[3,]  NaN  NaN  NaN
 > is.nan(solve(d))
Error in solve.default(d) : Lapack routine dgesv: system is exactly 
singular

This is against the spirit of IEEE-754 because it halts the program.

 > is.nan(det(d))
[1] FALSE

That is plain wrong.

>
> Many functions in R will actually bark at NaN inputs (e.g. qr, eigen, 
> ...) - maybe you're saying that we should check for NaNs in solve 
> before proceeding and raising an error?

However, this problem is in the Apple library not R.

Bill Northcott