[Rd] NaN and linear algebra

[Martin Maechler]

>>>>> "Bill" == Bill Northcott <w.northcott at unsw.edu.au>
>>>>>     on Wed, 23 Mar 2005 10:19:22 +1100 writes:

    Bill> 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.

    Bill> It may not violate the letter of IEEE-754 because matrix calculations 
    Bill> are not covered, but it certainly violates the spirit that arithmetic 
    Bill> 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.

    Bill> You may prefer the error, but it is not in the sprit of robust 
    Bill> arithmetic. ie
    >> d<-matrix(NaN,3,3)
    >> f<-solve(d)
    Bill> Error in solve.default(d) : Lapack routine dgesv: system is exactly 
    Bill> singular
    >> f
    Bill> 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?

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

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

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

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

    Bill> 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?

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

    Bill> Bill Northcott

Indeed!

I pretty much entirely agree with your points, Bill, and would
tend to declare that this Apple library is ``broken''
for building a correctly running R.

Let me ask one question I've been wondering about now for a
while:

  Did you run "make check" after building R,
  and "make check" ran to completion without an error?

If yes (which I doubt quite a bit), there *is* a bug in R's
quality control / quality assurance tools -- and I would want to
add a check for the misbehavior you've mentioned.

Martin Maechler, ETH Zurich