I am running into problems related, I think, to numerical precision of fast Fourier transforms. If I Fourier transform a gaussian distribution: absolute(fft(stats.norm.pdf(linspace(-10,10,512), loc=0, scale=1))) I find a floor of about 1e-16. Does anyone know of a way to improve the precision? Thanks, Darren
Darren Dale wrote:
I am running into problems related, I think, to numerical precision of fast Fourier transforms. If I Fourier transform a gaussian distribution:
absolute(fft(stats.norm.pdf(linspace(-10,10,512), loc=0, scale=1)))
I find a floor of about 1e-16. Does anyone know of a way to improve the precision?
1e-16 is the best you can do with double precision. -- Robert Kern robert.kern@gmail.com "In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
I have run into this as well. 1e-16 seems large to me for double precision. Robert can you explain this a bit more please. Is there an IEEE spec or something that specifies how much of the 64bits is for exponent and how much is for the mantissa (I think that is the right word)? I was playing with some FORTRAN code and it seems like there was a big difference with complex vs. double complex. It seems like 1e-16 was the magnitude floor of complex and double complex was better. Ryan On 12/21/05, Robert Kern <robert.kern@gmail.com> wrote:
Darren Dale wrote:
I am running into problems related, I think, to numerical precision of fast Fourier transforms. If I Fourier transform a gaussian distribution:
absolute(fft(stats.norm.pdf(linspace(-10,10,512), loc=0, scale=1)))
I find a floor of about 1e-16. Does anyone know of a way to improve the precision?
1e-16 is the best you can do with double precision.
-- Robert Kern robert.kern@gmail.com
"In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
Hi Ryan, Note that 1e-16 is the approximate _fractional_ precision of double precision floats. You can get absolute numbers much smaller than that. There are a bunch of good references to floating point accuracy etc. Some of them are linked here: http://babbage.cs.qc.edu/courses/cs341/IEEE-754references.html Note that the FFTW website lists FFT accuracies in both single and double precision here: http://www.fftw.org/accuracy/ Scott On Wednesday 21 December 2005 02:53 pm, Ryan Krauss wrote:
I have run into this as well. 1e-16 seems large to me for double precision. Robert can you explain this a bit more please. Is there an IEEE spec or something that specifies how much of the 64bits is for exponent and how much is for the mantissa (I think that is the right word)? I was playing with some FORTRAN code and it seems like there was a big difference with complex vs. double complex. It seems like 1e-16 was the magnitude floor of complex and double complex was better.
Ryan
On 12/21/05, Robert Kern <robert.kern@gmail.com> wrote:
Darren Dale wrote:
I am running into problems related, I think, to numerical precision of fast Fourier transforms. If I Fourier transform a gaussian distribution:
absolute(fft(stats.norm.pdf(linspace(-10,10,512), loc=0, scale=1)))
I find a floor of about 1e-16. Does anyone know of a way to improve the precision?
1e-16 is the best you can do with double precision.
-- Robert Kern robert.kern@gmail.com
"In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
-- Scott M. Ransom Address: NRAO Phone: (434) 296-0320 520 Edgemont Rd. email: sransom@nrao.edu Charlottesville, VA 22903 USA GPG Fingerprint: 06A9 9553 78BE 16DB 407B FFCA 9BFA B6FF FFD3 2989
Ryan Krauss wrote:
I have run into this as well. 1e-16 seems large to me for double precision. Robert can you explain this a bit more please. Is there an IEEE spec or something that specifies how much of the 64bits is for exponent and how much is for the mantissa (I think that is the right word)? I was playing with some FORTRAN code and it seems like there was a big difference with complex vs. double complex. It seems like 1e-16 was the magnitude floor of complex and double complex was better.
IEEE-754 double precision has 11 bits for exponent, 1 bit for sign, and 52 for the rest of the mantissa. 2.**-52 ~= 2.22e-16 Goldberg, David. What Every Computer Scientist Should Know About Floating Point Arithmetic. 1991. http://www.physics.ohio-state.edu/~dws/grouplinks/floating_point_math.pdf -- Robert Kern robert.kern@gmail.com "In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
Sorry to bring this back up again (this thread has been dormant since 12/21/05), but I think this is the key road block between me and really understanding this problem (and I think it went over my head when it was originally posted). If double precision numbers use 52 bits for the mantissa, what does it mean to say that 2**-52 is the approximate fractional precision of double precision floats? How does this affect optimization routines like fzero and fmin so that their accuracy seems limited to 2**-52 when the smallest possible floating point number should be more like 2**-1074. Does this somehow relate to how floating point subtraction is done (as relates to my post from earlier today). Thanks again, Ryan ========================= Hi Ryan, Note that 1e-16 is the approximate _fractional_ precision of double precision floats. You can get absolute numbers much smaller than that. There are a bunch of good references to floating point accuracy etc. Some of them are linked here: http://babbage.cs.qc.edu/courses/cs341/IEEE-754references.html Note that the FFTW website lists FFT accuracies in both single and double precision here: http://www.fftw.org/accuracy/ Scott - Show quoted text - On Wednesday 21 December 2005 02:53 pm, Ryan Krauss wrote:
I have run into this as well. 1e-16 seems large to me for double precision. Robert can you explain this a bit more please. Is there an IEEE spec or something that specifies how much of the 64bits is for exponent and how much is for the mantissa (I think that is the right word)? I was playing with some FORTRAN code and it seems like there was a big difference with complex vs. double complex. It seems like 1e-16 was the magnitude floor of complex and double complex was better.
Ryan
On 12/21/05, Robert Kern <robert.kern@gmail.com> wrote:
Darren Dale wrote:
I am running into problems related, I think, to numerical precision of fast Fourier transforms. If I Fourier transform a gaussian distribution:
absolute(fft(stats.norm.pdf(linspace(-10,10,512), loc=0, scale=1)))
I find a floor of about 1e-16. Does anyone know of a way to improve the precision?
1e-16 is the best you can do with double precision.
-- Robert Kern robert.kern@gmail.com
"In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
-- Scott M. Ransom Address: NRAO Phone: (434) 296-0320 520 Edgemont Rd. email: sransom@nrao.edu Charlottesville, VA 22903 USA GPG Fingerprint: 06A9 9553 78BE 16DB 407B FFCA 9BFA B6FF FFD3 2989 On 12/21/05, Robert Kern <robert.kern@gmail.com> wrote:
Ryan Krauss wrote:
I have run into this as well. 1e-16 seems large to me for double precision. Robert can you explain this a bit more please. Is there an IEEE spec or something that specifies how much of the 64bits is for exponent and how much is for the mantissa (I think that is the right word)? I was playing with some FORTRAN code and it seems like there was a big difference with complex vs. double complex. It seems like 1e-16 was the magnitude floor of complex and double complex was better.
IEEE-754 double precision has 11 bits for exponent, 1 bit for sign, and 52 for the rest of the mantissa. 2.**-52 ~= 2.22e-16
Goldberg, David. What Every Computer Scientist Should Know About Floating Point Arithmetic. 1991. http://www.physics.ohio-state.edu/~dws/grouplinks/floating_point_math.pdf
-- Robert Kern robert.kern@gmail.com
"In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
Ryan Krauss wrote:
Sorry to bring this back up again (this thread has been dormant since 12/21/05), but I think this is the key road block between me and really understanding this problem (and I think it went over my head when it was originally posted). If double precision numbers use 52 bits for the mantissa, what does it mean to say that 2**-52 is the approximate fractional precision of double precision floats? How does this affect optimization routines like fzero and fmin so that their accuracy seems limited to 2**-52 when the smallest possible floating point number should be more like 2**-1074. Does this somehow relate to how floating point subtraction is done (as relates to my post from earlier today).
I don't have time right now to go into a lecture on condition numbers and other such stuff, but let give you an example to chew on and a good reference: In [23]: a = rand(1024) In [24]: ifft(fft(a)) Out[24]: array([ 0.27261246 -4.33680869e-17j, 0.08520686 -1.04083409e-17j, 0.66654776 +6.60887690e-17j, ..., 0.05895311 +8.37004077e-17j, 0.40454606 -2.34796237e-17j, 0.38051768 +8.63024929e-17j]) In [25]: a Out[25]: array([ 0.27261246, 0.08520686, 0.66654776, ..., 0.05895311, 0.40454606, 0.38051768]) In [26]: a *= 1e16 In [27]: ifft(fft(a)) Out[27]: array([ 2.72612457e+15+0.21875j , 8.52068614e+14-0.359375j , 6.66547756e+15-0.11579423j, ..., 5.89531079e+14+1.16015625j, 4.04546057e+15+0.7886458j , 3.80517684e+15+0.71875j ]) The first volume of _Numerical Computation_ by Christoph W. Ueberhuber (_Computer-Numerik 1_ in the original German) covers these topics quite well and extensively. -- Robert Kern robert.kern@gmail.com "In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
Robert Kern wrote:
Ryan Krauss wrote:
Sorry to bring this back up again (this thread has been dormant since 12/21/05), but I think this is the key road block between me and really understanding this problem (and I think it went over my head when it was originally posted). If double precision numbers use 52 bits for the mantissa, what does it mean to say that 2**-52 is the approximate fractional precision of double precision floats?
I don't have time right now to go into a lecture on condition numbers and other such stuff, but let give you an example to chew on and a good reference:
If I can jump in and show something even simpler: In [1]: x,sum=0.1,0 In [2]: for m in range(10): ...: sum = sum + x ...: In [3]: sum Out[3]: 0.99999999999999989 In [4]: sum-1 Out[4]: -1.1102230246251565e-16 As for books, I've always liked Kahaner, Moler and Nash. Steve
Thanks, I am feeling much better about this. I actually found something like this on the internet that is really helpful to me: In [1]: m=1.0 In [2]: e=m In [3]: while m+e>m: ...: e=e/2.0 ...: In [4]: e Out[4]: 1.1102230246251565e-16 In [5]: m=1e-4 In [6]: e=m In [8]: while m+e>m: ...: e=e/2.0 ...: In [9]: e Out[9]: 5.551115123125783e-21 In [10]: m*(2.0**-52) Out[10]: 2.2204460492503132e-20 Ryan On 2/15/06, Stephen Walton <stephen.walton@csun.edu> wrote:
Robert Kern wrote:
Ryan Krauss wrote:
Sorry to bring this back up again (this thread has been dormant since 12/21/05), but I think this is the key road block between me and really understanding this problem (and I think it went over my head when it was originally posted). If double precision numbers use 52 bits for the mantissa, what does it mean to say that 2**-52 is the approximate fractional precision of double precision floats?
I don't have time right now to go into a lecture on condition numbers and other such stuff, but let give you an example to chew on and a good reference:
If I can jump in and show something even simpler:
In [1]: x,sum=0.1,0
In [2]: for m in range(10): ...: sum = sum + x ...:
In [3]: sum Out[3]: 0.99999999999999989
In [4]: sum-1 Out[4]: -1.1102230246251565e-16
As for books, I've always liked Kahaner, Moler and Nash.
Steve
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
Ryan Krauss wrote:
Thanks, I am feeling much better about this.
I actually found something like this on the internet that is really helpful to me: In [1]: m=1.0 In [2]: e=m In [3]: while m+e>m: ...: e=e/2.0 ...:
In [4]: e Out[4]: 1.1102230246251565e-16
In [5]: m=1e-4 In [6]: e=m
In [8]: while m+e>m: ...: e=e/2.0 ...:
In [9]: e Out[9]: 5.551115123125783e-21
In [10]: m*(2.0**-52) Out[10]: 2.2204460492503132e-20
Ryan
There is some code in numpy/lib/machar.py that you might find useful in this regard. -Travis
Robert Kern <robert.kern@gmail.com> writes:
Ryan Krauss wrote:
Sorry to bring this back up again (this thread has been dormant since 12/21/05), but I think this is the key road block between me and really understanding this problem (and I think it went over my head when it was originally posted). If double precision numbers use 52 bits for the mantissa, what does it mean to say that 2**-52 is the approximate fractional precision of double precision floats? How does this affect optimization routines like fzero and fmin so that their accuracy seems limited to 2**-52 when the smallest possible floating point number should be more like 2**-1074. Does this somehow relate to how floating point subtraction is done (as relates to my post from earlier today).
I don't have time right now to go into a lecture on condition numbers and other such stuff, but let give you an example to chew on and a good reference:
[snip example]
The first volume of _Numerical Computation_ by Christoph W. Ueberhuber (_Computer-Numerik 1_ in the original German) covers these topics quite well and extensively.
Another reference: I've been working through _Accuracy and Stability of Numerical Algorithms_ (2nd ed.), by Nicholas J. Higham, which I'm finding to be pretty good. Ryan mentioned that he was working through "What Every Computer Scientist Should Know About Floating Point" in another thread, which is good place to start. Looking at the stuff on William Kahan's site at http://www.cs.berkeley.edu/~wkahan/ can give you a good idea of what problems you can run into with floating point. Basically, when doing numerical computations, there are two levels of analysis: numerical analysis, where you're concerned about which algorithm to use, how stable it is, etc., and then the nitty-gritty things about floating point. For the last bit, you mainly worry about if rounding errors, etc. are big enough to worry about :) (In my experience, most people don't do that bit, but they're usually ok, except when they're not.) -- |>|\/|< /--------------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm@physics.mcmaster.ca
David M. Cooke wrote:
Basically, when doing numerical computations, there are two levels of
analysis: numerical analysis, where you're concerned about which algorithm to use, how stable it is, etc., and then the nitty-gritty things about floating point. For the last bit, you mainly worry about if rounding errors, etc. are big enough to worry about :) (In my experience, most people don't do that bit, but they're usually ok, except when they're not.)
Which, to bring the discussion back to SciPy, it is such an excellent idea that the basic algorithms in it are from the well tested NETLIB and related routine libraries. The people who wrote those did worry about rounding errors and the like, and the code is well crafted. Don't re-invent the wheel and all that.
participants (7)
-
cookedm＠physics.mcmaster.ca -
Darren Dale -
Robert Kern -
Ryan Krauss -
Scott Ransom -
Stephen Walton -
Travis Oliphant