R devel - Mar 2005 - NaN and linear algebra

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

Bill,

On Mar 21, 2005, at 10:31 PM, Bill Northcott 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.

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.

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?

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?


Cheers,
Simon

On 24/03/2005, at 10:04 PM, Andy wrote:
> I just tested this problem on a Windows 2000 Pro P4-Xeon system.  I 
> tried
> this in patched 2.0.1 (2005-2-28) and 2.1.0-alpha (2005-3-23) with both
> generic BLAS and precompiled ATLAS BLAS dll (downloaded from CRAN 
> windows
> binaries).  In all 4 combinations I get
>
>> d<-matrix(NaN,3,3)
>> f<-solve(d)
>> f
>      [,1] [,2] [,3]
> [1,]  NaN  NaN  NaN
> [2,]  NaN  NaN  NaN
> [3,]  NaN  NaN  NaN
>> det(d)
> [1] NaN
>>
>>
I just want to clarify some of the discussion here for my own benefit.

As I understand it, there are two software libraries involved (three if 
you count R).  BLAS is the lowest layer providing basic matrix 
operations.  It comes in vanilla and platform optimised (ATLAS or Goto) 
flavours.  Presumably R's built in version is vanilla.

Relying on BLAS is LAPACK which provides the functions which R actually 
calls like dgesv.  LAPACK relies on BLAS to provide the platform 
optimisations.  As far as I can see the issue(s) here is(are) with 
LAPACK not BLAS.  It is LAPACK that issues the error and halts the 
program.  (The MacOS X Accelerate Framework contains BLAS and LAPACK as 
distinct libraries libBlas.dylib and libLapack.dylib)

If, as reported, the behaviour is fundamentally similar on MacOS X, 
Linux and IRIX then it would appear that this is a feature(bug?) of the 
LAPACK reference implementation and that whoever wrote the Windows 
version saw fit to improve on it.  I don't think LAPACK ever pretended 
to properly handle NaNs.

It seems to me that the Windows version is an improvement, because it 
complies with the robust arithmetic principles that programs should not 
halt but propagate NaN or other indicators that can tested as 
appropriate in the application.  So it might be good idea to press 
other implementers of LAPACK to do the same.

My view for what it is worth is that R should also be robust in this 
sense and not halt on NaNs.  For preference this should be handled in 
the LAPACK implementations, but if not R should should trap them.  As I 
see it, it is incorrect to say a matrix containing NaNs is singular, 
not positive definite or whatever, because since it contains NaNs the 
computations to make those determinations cannot be performed.

Is this the wrong way to look at it?

Bill Northcott

PS  This discussion has given me another thought.  On MacOS X, it would 
be much faster to use floats instead of doubles because they can be 
directly processed by the Altivec section of the CPU.  I can see no way 
to have R use floats.  Would that be possible as a future enhancement 
or is it just too hard?