[Rd] NaN and linear algebra

[apjaworski at mmm.com]

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
>

which, I guess, is what is expected.

Hope this helps,

Andy

__________________________________
Andy Jaworski
518-1-01
Process Laboratory
3M Corporate Research Laboratory
-----
E-mail: apjaworski at mmm.com
Tel:  (651) 733-6092
Fax:  (651) 736-3122


                                                                           
             Prof Brian Ripley                                             
             <ripley at stats.ox.                                             
             ac.uk>                                                     To 
             Sent by:                  Simon Urbanek                       
             r-devel-bounces at s         <simon.urbanek at r-project.org>       
             tat.math.ethz.ch                                           cc 
                                       Martin Maechler                     
                                       <maechler at stat.math.ethz.ch>,       
             03/23/2005 01:28          r-devel at stat.math.ethz.ch           
             PM                                                    Subject 
                                       Re: [Rd] NaN and linear algebra     
                                                                           
                                                                           
                                                                           
                                                                           
                                                                           
                                                                           




On Wed, 23 Mar 2005, Simon Urbanek wrote:

> On Mar 22, 2005, at 6:19 PM, Bill Northcott wrote:
>
>> On 23/03/2005, at 12:55 AM, Simon Urbanek wrote:
>>
>> 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
>
> For the record: this is the behavior on all platforms I tested (Mac,
Linux,
> IRIX) - the only platform with a different result is reportedly Windows.
So
> what I was saying is that this is not "IEEE-754 non-compliance of MacOS
X" as
> you put it.

Hmm: all the platforms I have tested using R's internal blas and some
with external BLASes pass this test.  Only those with ATLAS and Goto
BLAS' fail for me.

>> Clearly det(d) returning 0 is wrong.  As a result based on a computation

>> including a NaN, it should return NaN.
>
> That's exactly what I was pointing out, and, yes, this is a bug in
vecLib.
>
> @Martin: I guess the most simple test for this would be probably
> is.nan(det(matrix(NaN,2,2)))
>
>>> 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.
>
> Since you're referencing the "solve" problem here, again, it's not.
>
> As I was explaining the the previous e-mail, there are at least two
> completely separate issues - handling of NaNs in solve and determinant of
NaN
> matrices. The latter is a bug in vecLib, the first one is not an OS X
> specific problem.

No, but it does appears to be an external BLAS-specific problem.

> Many R functions will object to taking NaN inputs such as qr.solve which
> was the previous implementation of solve. At any rate this is error
> handling in R and not Lapack. The "example" on the top doesn't make any
> sense because when you type "f" you know that it doesn't exist and both
> are R errors. If you used that in a script you have both choices: abort
> and that point or continue with whatever fallback you choose.

--
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

______________________________________________
R-devel at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel