Numerical instability in new R Windows development version

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Numerical instability in new R Windows development version

Hans W Borchers
I have a question concerning the new Windows toolchain for R >= 2.14.2.
When trying out my package 'pracma' on the win-builder development version
it will stop with the following error message:

  > f3 <- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2 <= 1))
  > dblquad(f3, -1, 1, -1, 1)     #   2.094395124 , i.e. 2/3*pi , err = 2e-8
  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2 <= 1)) : NaNs produced
  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2 <= 1)) : NaNs produced
  Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs,  :
    non-finite function value
  Calls: dblquad ...
         <Anonymous> -> f -> do.call -> mapply -> <Anonymous> -> integrate
  Execution halted
  ** running examples for arch 'x64' ... ERROR
  Running examples in 'pracma-Ex.R' failed

This probably means that the following expression got negative for some
values x, y:

  (1 - (x^2 + y^2)) * (x^2 + y^2 <= 1)

It appears to be an often used trick in numerical analysis. One advantage is
that a function using it is immediately vectorized while an expression such
as, e.g., "max(0, 1 - (x^2 + y^2))" is not.

The example runs fine on Debian Linux and Mac OS X 32-/64-bit architectures.
In my understanding the approach is correct and, as said above, often used in
numerical applications.

Can someone explain to me why this fails for the Windows 64-bit compiler and
what I should use instead. Thanks.

Hans Werner Borchers
ABB Corporate Research

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Numerical instability in new R Windows development version

Duncan Murdoch-2
On 12-01-27 7:23 AM, Hans W Borchers wrote:

> I have a question concerning the new Windows toolchain for R>= 2.14.2.
> When trying out my package 'pracma' on the win-builder development version
> it will stop with the following error message:
>
>    >  f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1))
>    >  dblquad(f3, -1, 1, -1, 1)     #   2.094395124 , i.e. 2/3*pi , err = 2e-8
>    Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>    Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>    Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs,  :
>      non-finite function value
>    Calls: dblquad ...
>           <Anonymous>  ->  f ->  do.call ->  mapply ->  <Anonymous>  ->  integrate
>    Execution halted
>    ** running examples for arch 'x64' ... ERROR
>    Running examples in 'pracma-Ex.R' failed
>
> This probably means that the following expression got negative for some
> values x, y:
>
>    (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)

I think you're right, it's a bug, hopefully easy to fix.  Here's a
simpler version:

x <- 0*(-1)
sqrt(x)

x is a "negative zero", and the sqrt() function incorrectly produces a
NaN in the new toolchain.

Duncan Murdoch

>
> It appears to be an often used trick in numerical analysis. One advantage is
> that a function using it is immediately vectorized while an expression such
> as, e.g., "max(0, 1 - (x^2 + y^2))" is not.
>
> The example runs fine on Debian Linux and Mac OS X 32-/64-bit architectures.
> In my understanding the approach is correct and, as said above, often used in
> numerical applications.
>
> Can someone explain to me why this fails for the Windows 64-bit compiler and
> what I should use instead. Thanks.
>
> Hans Werner Borchers
> ABB Corporate Research
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Numerical instability in new R Windows development version

Peter Dalgaard-2
In reply to this post by Hans W Borchers

On Jan 27, 2012, at 13:23 , Hans W Borchers wrote:

>  (1 - (x^2 + y^2)) * (x^2 + y^2 <= 1)
>
> It appears to be an often used trick in numerical analysis. One advantage is
> that a function using it is immediately vectorized while an expression such
> as, e.g., "max(0, 1 - (x^2 + y^2))" is not.

However, "pmax(0, 1 - (x^2 + y^2))" is (unless 0-length x,y is an issue).

But of course, Duncan is right: It is a bug if you can't take the square root of negative zero.

--
Peter Dalgaard, Professor
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Numerical instability in new R Windows development version

Prof Brian Ripley
In reply to this post by Duncan Murdoch-2
On 27/01/2012 13:26, Duncan Murdoch wrote:

> On 12-01-27 7:23 AM, Hans W Borchers wrote:
>> I have a question concerning the new Windows toolchain for R>= 2.14.2.
>> When trying out my package 'pracma' on the win-builder development
>> version
>> it will stop with the following error message:
>>
>> > f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1))
>> > dblquad(f3, -1, 1, -1, 1) # 2.094395124 , i.e. 2/3*pi , err = 2e-8
>> Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>> Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>> Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs, :
>> non-finite function value
>> Calls: dblquad ...
>> <Anonymous> -> f -> do.call -> mapply -> <Anonymous> -> integrate
>> Execution halted
>> ** running examples for arch 'x64' ... ERROR
>> Running examples in 'pracma-Ex.R' failed
>>
>> This probably means that the following expression got negative for some
>> values x, y:
>>
>> (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)
>
> I think you're right, it's a bug, hopefully easy to fix. Here's a
> simpler version:
>
> x <- 0*(-1)
> sqrt(x)
>
> x is a "negative zero", and the sqrt() function incorrectly produces a
> NaN in the new toolchain.

Well, for some definition of 'incorrectly'.  It is clearly what the
author of that piece of code intended.

It would be helpful if people would cite definitive references.  Someone
is going to have to report this on the bugtracker, and at present I
don't have enough evidence to do so: the C99/C11 standards do not seem
to mandate a particular value (they do say what happens for values less
than zero, but C compilers are allowed to have or not have signed
zeroes).  (Various Unix-alikes say what they do, usually -0, but that's
not evidence that other answers are 'incorrect'.)

> Duncan Murdoch
>
>>
>> It appears to be an often used trick in numerical analysis. One
>> advantage is
>> that a function using it is immediately vectorized while an expression
>> such
>> as, e.g., "max(0, 1 - (x^2 + y^2))" is not.
>>
>> The example runs fine on Debian Linux and Mac OS X 32-/64-bit
>> architectures.
>> In my understanding the approach is correct and, as said above, often
>> used in
>> numerical applications.
>>
>> Can someone explain to me why this fails for the Windows 64-bit
>> compiler and
>> what I should use instead. Thanks.
>>
>> Hans Werner Borchers
>> ABB Corporate Research
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel


--
Brian D. Ripley,                  [hidden email]
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Numerical instability in new R Windows development version

Duncan Murdoch-2
On 27/01/2012 12:32 PM, Prof Brian Ripley wrote:

> On 27/01/2012 13:26, Duncan Murdoch wrote:
> >  On 12-01-27 7:23 AM, Hans W Borchers wrote:
> >>  I have a question concerning the new Windows toolchain for R>= 2.14.2.
> >>  When trying out my package 'pracma' on the win-builder development
> >>  version
> >>  it will stop with the following error message:
> >>
> >>  >  f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1))
> >>  >  dblquad(f3, -1, 1, -1, 1) # 2.094395124 , i.e. 2/3*pi , err = 2e-8
> >>  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
> >>  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
> >>  Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs, :
> >>  non-finite function value
> >>  Calls: dblquad ...
> >>  <Anonymous>  ->  f ->  do.call ->  mapply ->  <Anonymous>  ->  integrate
> >>  Execution halted
> >>  ** running examples for arch 'x64' ... ERROR
> >>  Running examples in 'pracma-Ex.R' failed
> >>
> >>  This probably means that the following expression got negative for some
> >>  values x, y:
> >>
> >>  (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)
> >
> >  I think you're right, it's a bug, hopefully easy to fix. Here's a
> >  simpler version:
> >
> >  x<- 0*(-1)
> >  sqrt(x)
> >
> >  x is a "negative zero", and the sqrt() function incorrectly produces a
> >  NaN in the new toolchain.
>
> Well, for some definition of 'incorrectly'.  It is clearly what the
> author of that piece of code intended.
>
> It would be helpful if people would cite definitive references.  Someone
> is going to have to report this on the bugtracker, and at present I
> don't have enough evidence to do so: the C99/C11 standards do not seem
> to mandate a particular value (they do say what happens for values less
> than zero, but C compilers are allowed to have or not have signed
> zeroes).  (Various Unix-alikes say what they do, usually -0, but that's
> not evidence that other answers are 'incorrect'.)

Section 6.3 of IEEE 754-2008 says

Except that squareRoot(−0) shall be −0, every numeric squareRoot result
shall have a positive sign.

Duncan Murdoch

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Numerical instability in new R Windows development version

Duncan Murdoch-2
On 27/01/2012 1:26 PM, Duncan Murdoch wrote:

> On 27/01/2012 12:32 PM, Prof Brian Ripley wrote:
> >  On 27/01/2012 13:26, Duncan Murdoch wrote:
> >  >   On 12-01-27 7:23 AM, Hans W Borchers wrote:
> >  >>   I have a question concerning the new Windows toolchain for R>= 2.14.2.
> >  >>   When trying out my package 'pracma' on the win-builder development
> >  >>   version
> >  >>   it will stop with the following error message:
> >  >>
> >  >>   >   f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1))
> >  >>   >   dblquad(f3, -1, 1, -1, 1) # 2.094395124 , i.e. 2/3*pi , err = 2e-8
> >  >>   Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
> >  >>   Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
> >  >>   Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs, :
> >  >>   non-finite function value
> >  >>   Calls: dblquad ...
> >  >>   <Anonymous>   ->   f ->   do.call ->   mapply ->   <Anonymous>   ->   integrate
> >  >>   Execution halted
> >  >>   ** running examples for arch 'x64' ... ERROR
> >  >>   Running examples in 'pracma-Ex.R' failed
> >  >>
> >  >>   This probably means that the following expression got negative for some
> >  >>   values x, y:
> >  >>
> >  >>   (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)
> >  >
> >  >   I think you're right, it's a bug, hopefully easy to fix. Here's a
> >  >   simpler version:
> >  >
> >  >   x<- 0*(-1)
> >  >   sqrt(x)
> >  >
> >  >   x is a "negative zero", and the sqrt() function incorrectly produces a
> >  >   NaN in the new toolchain.
> >
> >  Well, for some definition of 'incorrectly'.  It is clearly what the
> >  author of that piece of code intended.
> >
> >  It would be helpful if people would cite definitive references.  Someone
> >  is going to have to report this on the bugtracker, and at present I
> >  don't have enough evidence to do so: the C99/C11 standards do not seem
> >  to mandate a particular value (they do say what happens for values less
> >  than zero, but C compilers are allowed to have or not have signed
> >  zeroes).  (Various Unix-alikes say what they do, usually -0, but that's
> >  not evidence that other answers are 'incorrect'.)
>
> Section 6.3 of IEEE 754-2008 says
>
> Except that squareRoot(−0) shall be −0, every numeric squareRoot result
> shall have a positive sign.

I believe the corresponding ISO standard is ISO/IEC/IEEE 60559:2011, but I don't have a copy, and I don't think my library does.

Duncan Murdoch

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Numerical instability in new R Windows development version

Duncan Murdoch-2
In reply to this post by Prof Brian Ripley
On 27/01/2012 12:32 PM, Prof Brian Ripley wrote:

> On 27/01/2012 13:26, Duncan Murdoch wrote:
> >  On 12-01-27 7:23 AM, Hans W Borchers wrote:
> >>  I have a question concerning the new Windows toolchain for R>= 2.14.2.
> >>  When trying out my package 'pracma' on the win-builder development
> >>  version
> >>  it will stop with the following error message:
> >>
> >>  >  f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1))
> >>  >  dblquad(f3, -1, 1, -1, 1) # 2.094395124 , i.e. 2/3*pi , err = 2e-8
> >>  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
> >>  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
> >>  Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs, :
> >>  non-finite function value
> >>  Calls: dblquad ...
> >>  <Anonymous>  ->  f ->  do.call ->  mapply ->  <Anonymous>  ->  integrate
> >>  Execution halted
> >>  ** running examples for arch 'x64' ... ERROR
> >>  Running examples in 'pracma-Ex.R' failed
> >>
> >>  This probably means that the following expression got negative for some
> >>  values x, y:
> >>
> >>  (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)
> >
> >  I think you're right, it's a bug, hopefully easy to fix. Here's a
> >  simpler version:
> >
> >  x<- 0*(-1)
> >  sqrt(x)
> >
> >  x is a "negative zero", and the sqrt() function incorrectly produces a
> >  NaN in the new toolchain.
>
> Well, for some definition of 'incorrectly'.  It is clearly what the
> author of that piece of code intended.
>
> It would be helpful if people would cite definitive references.  Someone
> is going to have to report this on the bugtracker, and at present I
> don't have enough evidence to do so: the C99/C11 standards do not seem
> to mandate a particular value (they do say what happens for values less
> than zero, but C compilers are allowed to have or not have signed
> zeroes).  (Various Unix-alikes say what they do, usually -0, but that's
> not evidence that other answers are 'incorrect'.)

This page:

http://pubs.opengroup.org/onlinepubs/9699919799/functions/sqrt.html

also says that the correct answer is -0.  I don't know if that has any
authority at all...

Duncan Murdoch

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Numerical instability in new R Windows development version

Stephen Weston
On Fri, Jan 27, 2012 at 1:40 PM, Duncan Murdoch
<[hidden email]> wrote:

> On 27/01/2012 12:32 PM, Prof Brian Ripley wrote:
>>
>> On 27/01/2012 13:26, Duncan Murdoch wrote:
>> >  On 12-01-27 7:23 AM, Hans W Borchers wrote:
>> >>  I have a question concerning the new Windows toolchain for R>= 2.14.2.
>> >>  When trying out my package 'pracma' on the win-builder development
>> >>  version
>> >>  it will stop with the following error message:
>> >>
>> >>  >  f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1))
>> >>  >  dblquad(f3, -1, 1, -1, 1) # 2.094395124 , i.e. 2/3*pi , err = 2e-8
>> >>  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>> >>  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>> >>  Error in integrate(function(y) f(x, y), ya, yb, subdivisions =
>> >> subdivs, :
>> >>  non-finite function value
>> >>  Calls: dblquad ...
>> >>  <Anonymous>  ->  f ->  do.call ->  mapply ->  <Anonymous>  ->
>> >>  integrate
>> >>  Execution halted
>> >>  ** running examples for arch 'x64' ... ERROR
>> >>  Running examples in 'pracma-Ex.R' failed
>> >>
>> >>  This probably means that the following expression got negative for
>> >> some
>> >>  values x, y:
>> >>
>> >>  (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)
>> >
>> >  I think you're right, it's a bug, hopefully easy to fix. Here's a
>> >  simpler version:
>> >
>> >  x<- 0*(-1)
>> >  sqrt(x)
>> >
>> >  x is a "negative zero", and the sqrt() function incorrectly produces a
>> >  NaN in the new toolchain.
>>
>> Well, for some definition of 'incorrectly'.  It is clearly what the
>> author of that piece of code intended.
>>
>> It would be helpful if people would cite definitive references.  Someone
>> is going to have to report this on the bugtracker, and at present I
>> don't have enough evidence to do so: the C99/C11 standards do not seem
>> to mandate a particular value (they do say what happens for values less
>> than zero, but C compilers are allowed to have or not have signed
>> zeroes).  (Various Unix-alikes say what they do, usually -0, but that's
>> not evidence that other answers are 'incorrect'.)
>
>
> This page:
>
> http://pubs.opengroup.org/onlinepubs/9699919799/functions/sqrt.html
>
> also says that the correct answer is -0.  I don't know if that has any
> authority at all...
>
> Duncan Murdoch

Section 5.4.1 "Arithmetic operations" of "IEEE Standard for Floating-
Point Arithmetic", IEEE Std 754-2008 says:

  "The operation squareRoot(x) computes √ x. It has a positive sign for
all operands ≥ 0, except that squareRoot(−0) shall be −0."

- Steve

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Numerical instability in new R Windows development version

Duncan Murdoch-2
In reply to this post by Hans W Borchers
This did turn out to be a bug in the new toolchain, and Brian Ripley has
devised a patch and put together a new one.  I've uploaded a new
Rtools215.exe, which should be available for download tomorrow, and
builds of R-patched and R-devel will soon use it.  Everything takes a
while to propagate to the volunteers and systems that build binaries and
the mirrors, but we should all be up to date by the end of the week or so.

Thanks for the report!

Duncan Murdoch


On 27/01/2012 7:23 AM, Hans W Borchers wrote:

> I have a question concerning the new Windows toolchain for R>= 2.14.2.
> When trying out my package 'pracma' on the win-builder development version
> it will stop with the following error message:
>
>    >  f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1))
>    >  dblquad(f3, -1, 1, -1, 1)     #   2.094395124 , i.e. 2/3*pi , err = 2e-8
>    Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>    Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>    Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs,  :
>      non-finite function value
>    Calls: dblquad ...
>           <Anonymous>  ->  f ->  do.call ->  mapply ->  <Anonymous>  ->  integrate
>    Execution halted
>    ** running examples for arch 'x64' ... ERROR
>    Running examples in 'pracma-Ex.R' failed
>
> This probably means that the following expression got negative for some
> values x, y:
>
>    (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)
>
> It appears to be an often used trick in numerical analysis. One advantage is
> that a function using it is immediately vectorized while an expression such
> as, e.g., "max(0, 1 - (x^2 + y^2))" is not.
>
> The example runs fine on Debian Linux and Mac OS X 32-/64-bit architectures.
> In my understanding the approach is correct and, as said above, often used in
> numerical applications.
>
> Can someone explain to me why this fails for the Windows 64-bit compiler and
> what I should use instead. Thanks.
>
> Hans Werner Borchers
> ABB Corporate Research
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel