16 ms·

The End of Numeric Error: An interview with John L. Gustafson

148 points by ozdave 10y ago

dang 10y ago
Url changed from https://motherboard.vice.com/read/a-new-number-format-for-computers-could-nuke-approximation-errors-for-good https://motherboard.vice.com/read/a-new-number-format-for-co..., which points to this.
robinhoodexe 10y ago
Not that relevant, but >Because the Julia language is particularly good at adding new data types without sacrificing speed, several people jumped on making a port to Julia, including the MIT founders of that language. I'm really impressed with the Julia language so far and considering implementing some quantum chemistry models for many-body systems (Hartree-Fock and probably Coupled Cluster if I ever get that far). Mostly because a) Julia looks fun, b) FORTRAN is still king of speed (which I'll need all I can get of) but is a pain to write and c) the alternative to FORTRAN is C(++), where being OO isn't that much of a bonus. And Julie looks fun to write.
- dnautics 10y ago
  Julia is really impressive. So, I uh, was working on a julia implementation of unums till I stopped being funded for it... I'm currently writing unit and continuous tests for a computer chip that is being simulated - basically julia does the computation in x86, which suggests what the IEEE result should be, and then generates the assembler and passes it on to the assembler. It probably could have been easily written in something like python, ruby, or perl, but the modularness of Julia, functional programming, and (optional) explicit typing result in easy-to-read code and a sense of security that there are no shenanigans going on under the hood with the type system. Another cool thing is for more complicated code (working on a demonstration fast fourier transform) is that I can write the core FFT algorithm as a function with macros, and by swapping a different include that defines the macros differently, it can either run it as a virtual processor using internal Julia arrays that represent the memory and register files, or build up a file that contains the assembler code. Again, this is obviously possible in other languages [0], but Julia makes it look nice and easy to comprehend. [0] I'm thinking that ruby (especially with it's DSL-eyness) would be really good for this, but I abandoned ruby when it became too difficult (I'm lazy and dense) to install the correct version on my own computer much less figure out how to install it on a server for other people to use.
  - trsohmers 10y ago
    Here's a link to the WIP Julia unum implementation... https://github.com/REX-Computing/unumjl https://github.com/REX-Computing/unumjl
- azag0 10y ago
  I started a DFT code in Julia two years ago but didn't get too far. I always felt that I have to think how to make the code fast in Julia (types/vectorization), whereas I just know it will be fast in Fortran. So I always end up writing modern Fortran C-binded to Python. Was it just me not paying enough time to get used to Julia or was/is it a real problem?
  - dnautics 10y ago
    probably the best way to do it would have been to bind the fortran functions to julia.
  - acidflask 10y ago
    Have you actually profiled the Julia code you wrote and ensured that you have followed the performance tips in the manual? If you did it correctly all the time should be spent in linear algebra operations, and so it shouldn't matter what language you implemented the code in since the performance is primarily determined by your underlying BLAS library.
    - azag0 10y ago
      It's not that I couldn't make it fast. I followed just the steps you mentioned and got to a really fast code. But the process didn't feel intuitive.
- santaclaus 10y ago
  > FORTRAN is still king of speed (which I'll need all I can get of) but is a pain to write I'm surprised that, with the restrict keyword, one would see better code gen from Fortran than C. Isn't the lack of aliasing Fortran's primary advantage?
  - neutronicus 10y ago
    Leveraging "restrict" is not a trivial thing to do. If you're a researcher first and a software developer second (which I suspect is true of gp), the last thing you want to do is get into the weeds of inspecting gcc assembly to figure out where you didn't put a "restrict" that you should've.
  - deleted 10y ago
    [deleted]
- eli_gottlieb 10y ago
  I dunno. I just started doing probabilistic programming in Haskell yesterday, and it's turned out to be faster at it than every other PPL I've tried. It didn't help that WebChurch, Anglican, and Venture are all interpreted languages. Treating a sampler as a lazy stream of samples is a nice abstraction, too.
  - neutronicus 10y ago
    As a counterpoint, my attempts to write fast Haskell PDE solvers have been underwhelming.
    - jfoutz 10y ago
      Once you wrap your head around the Haskell way, it seems easy to get something working. IMHO that something will be minimal and pretty as well. Once i try to make things fast, it gets weird. There are people who can write very fast Haskell code, but it seems like all of those people worked on ghc. I don't have a good mental model of when, say, a thunk will get created, so i liberally sprinkle strictness annotations, and fool around with list fusion operations. I really like Haskell. Even with all the effort i've put into it, it seems like i've barely scratched the surface if i want performance. sigh
  - marmaduke 10y ago
    Did you try Stan? It compiles your model to SSE2 enabled c++
    - eli_gottlieb 10y ago
      For my current problem, I don't want the available inference methods to impose any restrictions on what sort of model I can write, so I'm using general PPL only.
      - marmaduke 10y ago
        You can write you model in Stan and sample without inference too.
- acidflask 10y ago
  It would be really interesting if you'd like to implement some quantum chemistry! Rick Muller (of PyQuante fame) has implemented a port to Julia, but I think there's a lot of room to use more language features. https://github.com/rpmuller/pyquante2/tree/master/julia https://github.com/rpmuller/pyquante2/tree/master/julia Let me know if you do want to do some quantum chemistry coding. It would be good for my soul to revisit some of my PhD roots...
conistonwater 10y ago
Isn't this yet another implementation of interval arithmetic? What am I missing? Calling his book "The End of Error" seems a bit much.
- Someone 10y ago
  No, it is making a distinction between exact and inexact numbers, and uses arbitrary precision. I don't see it gain wide acceptance. If your problem is fixed-num, use a bignum library. If it is not, chances are you will meet a number that isn't representable in this format very, very soon (1/3, 1/5, sqrt(2), PI, etc.) The moment you do, you will have to decide to how many bits of precision you want to compute that inexact number. Problem is that you cannot make that choice, unless you know exactly what you intend to do with the number. I bet most users will choose a format that their hardware handles natively: IEEE. There may be edge cases in numerical computing that benefit from this, especially when problems are ill-conditioned, but even there, I expect people to convert to conventional IEEE on the outside of their API. Also, I'm not sure there is much room in this space below quad-width floats (but I am not an expert)
  - barsonme 10y ago
    wouldn't most of the issues be solved with use of bcd/dcd? you get the speed (compared to arbitrary-precision), fixed-precision, and exact answers up to the precision limit?
    - bsder 10y ago
      Unfortunately, no. BCD as "binary coded decimal" is, effectively, useless. Even arbitrary precision libraries tend to be faster than BCD. "decimal floating point", which is what I think you meant, just changes what isn't representable. Instead of 1/5 causing a problem, 1/3 causes a problem. The "classic" problem about this is matrix operations. Specifically, matrix inversion. You can use interval arithmetic on matrix inversion, and you find that the intervals blow up to almost infinity after just a couple of operations even though the center of the interval is actually pretty close to the solution. Of course, the reverse also holds: while most matricies invert relatively cleanly, there are always pathological cases that really do need humungous precision to invert. PS: Don't know why you were getting downvoted as you asked a legitimate question.
muizelaar 10y ago
Here's a python port of Gustafson's Mathematica unum implementation that I did: https://github.com/jrmuizel/pyunum https://github.com/jrmuizel/pyunum
- jepler 10y ago
  Hi, thanks for sharing your implementation. Am I totally misusing (py)unum? Using pyunum, repeatedly rotate the point (1,0) by 45 degrees. Print out the X coordinate every 8 rotations. Mathematically, the value printed each time should be exactly 1. Unfortunately, after just 79 iterations, the unum covers the range (-27352064, 27352064)! This includes the correct result but is useless. A similar program using Python floats (i.e., platform doubles) prints 0.9999999999999937 at step 79 and runs in 1/540 of the time on my oldish core2 notebook. (or 0.9999999999371962 at step 799999 in less than 1/10 of the time of the 79 unum steps) from unum import * a = divideu(x2u(1), sqrtu(x2u(2))) x, y = x2u(1), x2u(0) for i in range(80): x, y = timesu(minusu(x, y), a), timesu(plusu(x, y), a) if i % 8 == 7: print i, view(x)
  - AlexCoventry 10y ago
    Yes, you're misusing the concept, and pyunum makes it hard to use in this situation (probably because the author had a different purpose in mind.) Try re-running (in a fresh python invocation if you're in a shell) with this at the top of your script: import unum_config unum_config.e = 3 unum_config.f = 9 This gives 2^3=8 bits for the exponent, and 2^9=512 bits for the fraction. The default values for e and f are 3 (8 bits) and 4 (16 bits) respectively. The idea with unums is that your program sees the unacceptable bounds, and recomputes with higher precision. Comparing speeds is very unfair: one implementation is hardware which has received maybe $100B-worth of optimization over decades, the other is one guy's python script. In this implementation, 9 seems to be the maximum value for f: if I set it to 10 unum fails to import because it tries to divide a float by too large an integer. And 4 seems to be the maximum practical value for e. At e=5, initialization of unum takes over two minutes and 250M of memory, probably because it's computing the largest possible value in that floating point scheme.
ris 10y ago
Call me when there's a hardware implementation.
- trsohmers 10y ago
  <6 months for FPGA, ~12-18 months for silicon...
  - r0fls 10y ago
    Is that in the article? Or else, how do you know? I would like to write an implementation for fun in another language, but unum is getting mixed reviews here. Hence, it would be nice to see your source for it being incorporated into hardware.
    - trsohmers 10y ago
      At the bottom of the article, John mentions REX Computing, which I am the founder/CEO of. John has been a friend and advisor, and we have played around with minimal unum implementations, but haven't been able to commit the resources to a full one until after we tapeout our first chip in a little over a month (which will be using IEEE float). We also funded "dnautics" (a HN user who also commented here) to do a soft implementation of unum 1.0 in Julia.
      - oso2k 10y ago
        Has REX explored other alternative binary floating point formats? Say, dec64 [0]? [0] http://www.dec64.com http://www.dec64.com
        r0fls 10y ago
        How is dec64 different from Floating Point (IEEE 754)? It seems IEEE 754 allows for using base 10. What else do they do differently, do you know?
        10y ago
ori_b 10y ago
From what I'm understanding, this format uses a variable number of bits, automatically growing and shrinking as needed. How can this be implemented efficiently? If there's a cap on the maximum number of bits (eg, 64), why not just use that maximum number all the time, and not have to worry about shuffling things around (bitwise!) in memory?
- dnautics 10y ago
  if it's supported in the hardware, you could turn off those gates altogether and save energy.
- troymc 10y ago
  The article answers that question.
ambrop7 10y ago
Why bother? The precision of double is sufficient for the vast majority of engineering problems. Seldom do you have measurements that are more precice than the calculations you do with them.
- dnautics 10y ago
  A big problem is that chained floating point operations are a O(sqrt n) martingale. Also, some floating point ops that use pade approximants (most log algorithms, sin, cos) lose onoe or two precision points (versus half for arithmetic ops) of floating point accuracy.
  - ambrop7 10y ago
    Yeah, I know, there are cases where it matters. But I suppose many of these cases are also the ones where performance (and a hardware implementation) matters, and you're better off designing your algorithm to deal with roundoff error than use a magically precise FP implementation.
    - dnautics 10y ago
      sure. You could go the KSP route and write your own FP. There are still problems like the various krakens. KSP is obviously not so critical that krakens are life- or project- jeopardizing... But imagine being NASA for example, designing real world satellite guidance systems that must operate over many orders of magnitude. If you're using IEEE floating points, your systems will silently fail and (rarely) give an error like NaN. IEEE intervals will work, but often give a diagnostically useless answer. Unums have characteristic signatures that hint where to look to find the problematic calculation. You can then backtrack and redo the calculation at a higher precision and drop back to the standard precision as necessary. This is, of course, automatable.
  - acidflask 10y ago
    > A big problem is that chained floating point operations are a O(sqrt n) martingale. Mostly false; see https://www.cs.berkeley.edu/~wkahan/Mindless.pdf https://www.cs.berkeley.edu/~wkahan/Mindless.pdf It is true that an arbitrary, independent sequence of n floating point operations can potentially lose up to 2n bits of precision. However, it is false that most scientific computing code has such structure. On the contrary, the loss of precision in numerical algorithms can be estimated quite well using standard techniques of error analysis, especially when there is underlying matrix structure to exploit. The growth of error in properly written numerical code should NOT grow exponentially with the number of floating point operations. It should be bounded by intrinsic properties of the underlying matrix, such as its condition number or spectral norm.
aaibec 10y ago
Did anyone find peer-reviewed papers on the topic (e.g. in an IEEE journal given that it challenges the IEEE floating point standard)? I think the book is nice, but popularizing "unum"s with a book targeting non-experts feels like skipping an important step in due process, given the grand claims that are being made.
- jlgustafson 10y ago
  Plenty of experts refereed the book, believe me. William Kahan, David Bailey, Gordon Bell, Horst Simon, Ulrich Kulisch, John Gunnels, and anonymous reviewers as well. It was vetted for months before it was released. If you prefer a more formal treatment, Ulrich Kulisch has written up the unum concepts in a paper more palatable to mathematicians. Incidentally, I submitted a short version to the ARITH23 conference detailing how IEEE 754 could be repaired to give bit-identical results on different systems, and they had it reviewed by three reviewers, one who liked it, and the other two who served on the IEEE 754 committee and trashed it for some pretty silly reasons, like using the word "subnormal" instead of "denormal". There is no easy way to boil the ocean, but I'm trying.

148 points by ozdave 10y ago

dang 10y ago

Url changed from https://motherboard.vice.com/read/a-new-number-format-for-computers-could-nuke-approximation-errors-for-good https://motherboard.vice.com/read/a-new-number-format-for-co..., which points to this.

robinhoodexe 10y ago

Not that relevant, but >Because the Julia language is particularly good at adding new data types without sacrificing speed, several people jumped on making a port to Julia, including the MIT founders of that language. I'm really impressed with the Julia language so far and considering implementing some quantum chemistry models for many-body systems (Hartree-Fock and probably Coupled Cluster if I ever get that far). Mostly because a) Julia looks fun, b) FORTRAN is still king of speed (which I'll need all I can get of) but is a pain to write and c) the alternative to FORTRAN is C(++), where being OO isn't that much of a bonus. And Julie looks fun to write.

dnautics 10y ago
Julia is really impressive. So, I uh, was working on a julia implementation of unums till I stopped being funded for it... I'm currently writing unit and continuous tests for a computer chip that is being simulated - basically julia does the computation in x86, which suggests what the IEEE result should be, and then generates the assembler and passes it on to the assembler. It probably could have been easily written in something like python, ruby, or perl, but the modularness of Julia, functional programming, and (optional) explicit typing result in easy-to-read code and a sense of security that there are no shenanigans going on under the hood with the type system. Another cool thing is for more complicated code (working on a demonstration fast fourier transform) is that I can write the core FFT algorithm as a function with macros, and by swapping a different include that defines the macros differently, it can either run it as a virtual processor using internal Julia arrays that represent the memory and register files, or build up a file that contains the assembler code. Again, this is obviously possible in other languages [0], but Julia makes it look nice and easy to comprehend. [0] I'm thinking that ruby (especially with it's DSL-eyness) would be really good for this, but I abandoned ruby when it became too difficult (I'm lazy and dense) to install the correct version on my own computer much less figure out how to install it on a server for other people to use.
- trsohmers 10y ago
  Here's a link to the WIP Julia unum implementation... https://github.com/REX-Computing/unumjl https://github.com/REX-Computing/unumjl
azag0 10y ago
I started a DFT code in Julia two years ago but didn't get too far. I always felt that I have to think how to make the code fast in Julia (types/vectorization), whereas I just know it will be fast in Fortran. So I always end up writing modern Fortran C-binded to Python. Was it just me not paying enough time to get used to Julia or was/is it a real problem?
- dnautics 10y ago
  probably the best way to do it would have been to bind the fortran functions to julia.
- acidflask 10y ago
  Have you actually profiled the Julia code you wrote and ensured that you have followed the performance tips in the manual? If you did it correctly all the time should be spent in linear algebra operations, and so it shouldn't matter what language you implemented the code in since the performance is primarily determined by your underlying BLAS library.
  - azag0 10y ago
    It's not that I couldn't make it fast. I followed just the steps you mentioned and got to a really fast code. But the process didn't feel intuitive.
santaclaus 10y ago
> FORTRAN is still king of speed (which I'll need all I can get of) but is a pain to write I'm surprised that, with the restrict keyword, one would see better code gen from Fortran than C. Isn't the lack of aliasing Fortran's primary advantage?
- neutronicus 10y ago
  Leveraging "restrict" is not a trivial thing to do. If you're a researcher first and a software developer second (which I suspect is true of gp), the last thing you want to do is get into the weeds of inspecting gcc assembly to figure out where you didn't put a "restrict" that you should've.
- deleted 10y ago
  [deleted]
eli_gottlieb 10y ago
I dunno. I just started doing probabilistic programming in Haskell yesterday, and it's turned out to be faster at it than every other PPL I've tried. It didn't help that WebChurch, Anglican, and Venture are all interpreted languages. Treating a sampler as a lazy stream of samples is a nice abstraction, too.
- neutronicus 10y ago
  As a counterpoint, my attempts to write fast Haskell PDE solvers have been underwhelming.
  - jfoutz 10y ago
    Once you wrap your head around the Haskell way, it seems easy to get something working. IMHO that something will be minimal and pretty as well. Once i try to make things fast, it gets weird. There are people who can write very fast Haskell code, but it seems like all of those people worked on ghc. I don't have a good mental model of when, say, a thunk will get created, so i liberally sprinkle strictness annotations, and fool around with list fusion operations. I really like Haskell. Even with all the effort i've put into it, it seems like i've barely scratched the surface if i want performance. sigh
- marmaduke 10y ago
  Did you try Stan? It compiles your model to SSE2 enabled c++
  - eli_gottlieb 10y ago
    For my current problem, I don't want the available inference methods to impose any restrictions on what sort of model I can write, so I'm using general PPL only.
    - marmaduke 10y ago
      You can write you model in Stan and sample without inference too.
acidflask 10y ago
It would be really interesting if you'd like to implement some quantum chemistry! Rick Muller (of PyQuante fame) has implemented a port to Julia, but I think there's a lot of room to use more language features. https://github.com/rpmuller/pyquante2/tree/master/julia https://github.com/rpmuller/pyquante2/tree/master/julia Let me know if you do want to do some quantum chemistry coding. It would be good for my soul to revisit some of my PhD roots...

conistonwater 10y ago

Isn't this yet another implementation of interval arithmetic? What am I missing? Calling his book "The End of Error" seems a bit much.

Someone 10y ago
No, it is making a distinction between exact and inexact numbers, and uses arbitrary precision. I don't see it gain wide acceptance. If your problem is fixed-num, use a bignum library. If it is not, chances are you will meet a number that isn't representable in this format very, very soon (1/3, 1/5, sqrt(2), PI, etc.) The moment you do, you will have to decide to how many bits of precision you want to compute that inexact number. Problem is that you cannot make that choice, unless you know exactly what you intend to do with the number. I bet most users will choose a format that their hardware handles natively: IEEE. There may be edge cases in numerical computing that benefit from this, especially when problems are ill-conditioned, but even there, I expect people to convert to conventional IEEE on the outside of their API. Also, I'm not sure there is much room in this space below quad-width floats (but I am not an expert)
- barsonme 10y ago
  wouldn't most of the issues be solved with use of bcd/dcd? you get the speed (compared to arbitrary-precision), fixed-precision, and exact answers up to the precision limit?
  - bsder 10y ago
    Unfortunately, no. BCD as "binary coded decimal" is, effectively, useless. Even arbitrary precision libraries tend to be faster than BCD. "decimal floating point", which is what I think you meant, just changes what isn't representable. Instead of 1/5 causing a problem, 1/3 causes a problem. The "classic" problem about this is matrix operations. Specifically, matrix inversion. You can use interval arithmetic on matrix inversion, and you find that the intervals blow up to almost infinity after just a couple of operations even though the center of the interval is actually pretty close to the solution. Of course, the reverse also holds: while most matricies invert relatively cleanly, there are always pathological cases that really do need humungous precision to invert. PS: Don't know why you were getting downvoted as you asked a legitimate question.

muizelaar 10y ago

Here's a python port of Gustafson's Mathematica unum implementation that I did: https://github.com/jrmuizel/pyunum https://github.com/jrmuizel/pyunum

jepler 10y ago
Hi, thanks for sharing your implementation. Am I totally misusing (py)unum? Using pyunum, repeatedly rotate the point (1,0) by 45 degrees. Print out the X coordinate every 8 rotations. Mathematically, the value printed each time should be exactly 1. Unfortunately, after just 79 iterations, the unum covers the range (-27352064, 27352064)! This includes the correct result but is useless. A similar program using Python floats (i.e., platform doubles) prints 0.9999999999999937 at step 79 and runs in 1/540 of the time on my oldish core2 notebook. (or 0.9999999999371962 at step 799999 in less than 1/10 of the time of the 79 unum steps) from unum import * a = divideu(x2u(1), sqrtu(x2u(2))) x, y = x2u(1), x2u(0) for i in range(80): x, y = timesu(minusu(x, y), a), timesu(plusu(x, y), a) if i % 8 == 7: print i, view(x)
- AlexCoventry 10y ago
  Yes, you're misusing the concept, and pyunum makes it hard to use in this situation (probably because the author had a different purpose in mind.) Try re-running (in a fresh python invocation if you're in a shell) with this at the top of your script: import unum_config unum_config.e = 3 unum_config.f = 9 This gives 2^3=8 bits for the exponent, and 2^9=512 bits for the fraction. The default values for e and f are 3 (8 bits) and 4 (16 bits) respectively. The idea with unums is that your program sees the unacceptable bounds, and recomputes with higher precision. Comparing speeds is very unfair: one implementation is hardware which has received maybe $100B-worth of optimization over decades, the other is one guy's python script. In this implementation, 9 seems to be the maximum value for f: if I set it to 10 unum fails to import because it tries to divide a float by too large an integer. And 4 seems to be the maximum practical value for e. At e=5, initialization of unum takes over two minutes and 250M of memory, probably because it's computing the largest possible value in that floating point scheme.

ris 10y ago

Call me when there's a hardware implementation.

trsohmers 10y ago
<6 months for FPGA, ~12-18 months for silicon...
- r0fls 10y ago
  Is that in the article? Or else, how do you know? I would like to write an implementation for fun in another language, but unum is getting mixed reviews here. Hence, it would be nice to see your source for it being incorporated into hardware.
  - trsohmers 10y ago
    At the bottom of the article, John mentions REX Computing, which I am the founder/CEO of. John has been a friend and advisor, and we have played around with minimal unum implementations, but haven't been able to commit the resources to a full one until after we tapeout our first chip in a little over a month (which will be using IEEE float). We also funded "dnautics" (a HN user who also commented here) to do a soft implementation of unum 1.0 in Julia.
    - oso2k 10y ago
      Has REX explored other alternative binary floating point formats? Say, dec64 [0]? [0] http://www.dec64.com http://www.dec64.com
      - r0fls 10y ago
        How is dec64 different from Floating Point (IEEE 754)? It seems IEEE 754 allows for using base 10. What else do they do differently, do you know?
        10y ago

ori_b 10y ago

From what I'm understanding, this format uses a variable number of bits, automatically growing and shrinking as needed. How can this be implemented efficiently? If there's a cap on the maximum number of bits (eg, 64), why not just use that maximum number all the time, and not have to worry about shuffling things around (bitwise!) in memory?

dnautics 10y ago
if it's supported in the hardware, you could turn off those gates altogether and save energy.
troymc 10y ago
The article answers that question.

ambrop7 10y ago

Why bother? The precision of double is sufficient for the vast majority of engineering problems. Seldom do you have measurements that are more precice than the calculations you do with them.

dnautics 10y ago
A big problem is that chained floating point operations are a O(sqrt n) martingale. Also, some floating point ops that use pade approximants (most log algorithms, sin, cos) lose onoe or two precision points (versus half for arithmetic ops) of floating point accuracy.
- ambrop7 10y ago
  Yeah, I know, there are cases where it matters. But I suppose many of these cases are also the ones where performance (and a hardware implementation) matters, and you're better off designing your algorithm to deal with roundoff error than use a magically precise FP implementation.
  - dnautics 10y ago
    sure. You could go the KSP route and write your own FP. There are still problems like the various krakens. KSP is obviously not so critical that krakens are life- or project- jeopardizing... But imagine being NASA for example, designing real world satellite guidance systems that must operate over many orders of magnitude. If you're using IEEE floating points, your systems will silently fail and (rarely) give an error like NaN. IEEE intervals will work, but often give a diagnostically useless answer. Unums have characteristic signatures that hint where to look to find the problematic calculation. You can then backtrack and redo the calculation at a higher precision and drop back to the standard precision as necessary. This is, of course, automatable.
- acidflask 10y ago
  > A big problem is that chained floating point operations are a O(sqrt n) martingale. Mostly false; see https://www.cs.berkeley.edu/~wkahan/Mindless.pdf https://www.cs.berkeley.edu/~wkahan/Mindless.pdf It is true that an arbitrary, independent sequence of n floating point operations can potentially lose up to 2n bits of precision. However, it is false that most scientific computing code has such structure. On the contrary, the loss of precision in numerical algorithms can be estimated quite well using standard techniques of error analysis, especially when there is underlying matrix structure to exploit. The growth of error in properly written numerical code should NOT grow exponentially with the number of floating point operations. It should be bounded by intrinsic properties of the underlying matrix, such as its condition number or spectral norm.

aaibec 10y ago

Did anyone find peer-reviewed papers on the topic (e.g. in an IEEE journal given that it challenges the IEEE floating point standard)? I think the book is nice, but popularizing "unum"s with a book targeting non-experts feels like skipping an important step in due process, given the grand claims that are being made.

jlgustafson 10y ago
Plenty of experts refereed the book, believe me. William Kahan, David Bailey, Gordon Bell, Horst Simon, Ulrich Kulisch, John Gunnels, and anonymous reviewers as well. It was vetted for months before it was released. If you prefer a more formal treatment, Ulrich Kulisch has written up the unum concepts in a paper more palatable to mathematicians. Incidentally, I submitted a short version to the ARITH23 conference detailing how IEEE 754 could be repaired to give bit-identical results on different systems, and they had it reviewed by three reviewers, one who liked it, and the other two who served on the IEEE 754 committee and trashed it for some pretty silly reasons, like using the word "subnormal" instead of "denormal". There is no easy way to boil the ocean, but I'm trying.