[Rd] NaN and linear algebra

[stefano iacus]

On 23/mar/05, at 09:04, Martin Maechler wrote:

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

Martin,
the use vecLib was (since now?) the suggested way to build R on OS X 
with blas support.
This is the way OS X binaries are distributed and of course they pass 
make check (I've checked again just now and both 2.0.0 and r-devel pass 
make check) so probably the NaN is a special case with don't test for 
LAPACK?
stefano


>
> 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
>
> ______________________________________________
> R-devel at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>