From simon.urbanek at r-project.org  Wed Mar 23 19:28:09 2005
From: simon.urbanek at r-project.org (Simon Urbanek)
Date: Wed Mar 23 19:28:28 2005
Subject: [Rd] NaN and linear algebra
In-Reply-To: <8d7b823ac46177052ba1cac59b6cba22@unsw.edu.au>
References: <200503211109.j2LB9Kbe018532@hypatia.math.ethz.ch>
	<098c2dce898207218d3e518868d7ffec@unsw.edu.au>
	<3E86E6D3-BB0E-42A1-A4BE-D00EF2F5C0EA@r-project.org>
	<8d7b823ac46177052ba1cac59b6cba22@unsw.edu.au>
Message-ID: <8F4C141E-635D-4447-ABBE-27175C0AC676@r-project.org>

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.

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

Cheers,
Simon


[Rd] NaN and linear algebra
Simon Urbanek