Floating point errors

On the CCL mailing list Sina Türeli has posted this question:

“I am working on a project related to proteins and my precision for coordinates is there digits after the dot. When I do operations like rotation around a dihedral, the dihedrals which shouldn’t change change at about 0.01 angles and so. I am afraid though that is not much it might accumulate over time. So do you have any suggestions for reducing floating integer errors? Would’t be feasible if I turned lets say the coordinate 21.567 to 21567 and do my operations? Or maybe even 215670? “

I have spent a lot of time fighting similar problems so I know how annoying this can be. To understand what is hapenning, look at how floating point numbers are stored and operated upon according to the IEEE 754 standard. Because of the binary representation, our decimal fractional numbers do not match exactly to float/double numbers. Even though, the 23 fraction bits in a floating point binary number map roughly to 6 decimal digit precision, it still does not mean that all 3 digit decimal numbers can be represented precisely. On the other hand, integers can be stored exactly, so that gives basis to the idea to store the number 21.567 as 21567 or 215670. This would work for storing and also for applying some basic arithmetic (addition, subtraction and multiplication) to numbers accurately without any error. However, division starts a problem and any trigonometric functions or sqare root function blows up the problem to be much worse than what you have with floating point numbers. Those functions produce irrational numbers i.e. they cannot be represented by a division of two integers. Unfortunately, rotations are typically defined by angles and the coordinate transformation requires the sine and cosine of the angle — depending on the task, sometimes the transformation required can be expressed in other ways, e.g. by quaternions, and sometimes we can compute transformations by simpler arithmetic if the goal is to transform some atoms to specific positions, e.g. an overlay (rather than rotation by a given angle). So, in short using integer fixed point representation will not solve the problem of coordinate drifting errors of 3D transformations, especially if rotations by a given angle is required. On the other hand, it can solve simpler problems.
How bad is the floating point problem and what can be solved by fixed point integer arithmetic ? Let me give you an example: you learned in school that a+b+c = c+a+b. Well, this simple rule breaks even for single digit precision floating point numbers! Consider the following simple C code (you can download the source, or a linux binary):

#include  int main() {
float a, b, c, d, s;
for (a = 0.1; a < 9.9; a += 0.2) {
for (b=0.4; b < 9.9; b += 0.1) {
for (c=0.7; c < 9.9; c += 0.3) {
s = a + b;
s += c;
d = c + a;
d += b;
if ( s != d ) {
printf( "a=%f, b=%f, c=%f, a+b+c=%f, c+a+b=%f hex: a=0x%08x b=0x%08x c=0x%08x s=0x%08x d=0x%08xn",
a, b, c, s, d, *(int *) &a,*(int *) &b,*(int *) &c,*(int *) &s, *(int *) &d );
}
}
}
}
}

Sorry for the lack of indentation, I could not convince WordPress to keep it pre-formatted.

You can see that if the summation in two orders were the same for all tested cases it would never print anything. However, when you run it, you get plenty of examples (34585 cases on my computer, about one third of all cases tested) printed where the rule breaks. Here are the first few examples you get:

a=0.100000, b=0.400000, c=3.700000, a+b+c=4.200000, c+a+b=4.199999
hex: a=0×3dcccccd b=0×3ecccccd c=0×406ccccb s=0×40866666 d=0×40866665
a=0.100000, b=0.400000, c=7.600002, a+b+c=8.100002, c+a+b=8.100001
hex: a=0×3dcccccd b=0×3ecccccd c=0×40f33337 s=0×4101999c d=0×4101999b
a=0.100000, b=0.500000, c=2.500000, a+b+c=3.100000, c+a+b=3.100000
hex: a=0×3dcccccd b=0×3f000000 c=0×401fffff s=0×40466666 d=0×40466665

As you can see from the hexadecimal version, the difference is only in the last 1 bit. Nevertheless, it is enough to throw off the result. Imagine, if you sum up score components, sort them and select only the best N solutions. Suddenly, you may keep or lose a solution depending on the order of summation. Now, that is scary…This problem cannot be solved by using double precision, but it can be solved simply by using fixed point integer representation. However, fixed point does not help for the rotation problem I am afraid.

ZZ

5 Responses to “Floating point errors”

  1. ChemSpiderMan Says:

    Thanks for the education! ANd keep this type of post coming. It’s likely educational for the majority of the people in the domain who don’t understand the specifics.

  2. SimBioSys Blog » Blog Archive » Quality in chemical software - the debate continues Says:

    […] Are the docking and QSAR study results reproducible ? With eHiTS and LASSO, the answer is definitely YES! I understand that many tools on the docking/QSAR market use stochastic (read random) methods and therefore their results are inherently unreproducible. Again, I can only speak with authority about our own software, which uses strictly deterministic and reproducible techniques. So if a different researcher in a different location runs our software on the same input they will get the same result. However, I do not see how one could run the “same calculation” using a different software. By definition, if you are using a different software (which embodies the calculation) then you are not running the same calculation. I can assure you the same is true for QM software as well, for the simple floating point error reasons I have explained in a previous blog post. So any different QM implementation will necessarily involve computation steps in different orders (as simple as summation in different order will suffice) and therefore get slightly different results. PMR: […]

  3. Unilever Centre for Molecular Informatics, Cambridge - petermr’s blog » Blog Archive » Quality is emerging in chemical software Says:

    […] PMR: ZZ addresses this below in reporting a competition and I’ll continue there Are the docking and QSAR study results reproducible ? With eHiTS and LASSO, the answer is definitely YES! I understand that many tools on the docking/QSAR market use stochastic (read random) methods and therefore their results are inherently unreproducible. Again, I can only speak with authority about our own software, which uses strictly deterministic and reproducible techniques. So if a different researcher in a different location runs our software on the same input they will get the same result. However, I do not see how one could run the “same calculation” using a different software. By definition, if you are using a different software (which embodies the calculation) then you are not running the same calculation. I can assure you the same is true for QM software as well, for the simple floating point error reasons I have explained in a previous blog post. So any different QM implementation will necessarily involve computation steps in different orders (as simple as summation in different order will suffice) and therefore get slightly different results. […]

  4. Trevor Says:

    “This problem cannot be solved by using double precision”

    Actually, it can. If you change the floats to doubles, the sums will be identical for at least the first 6 or so decimal places.

  5. zsolt Says:

    Trevor,

    No, unfortunately, you are wrong for 2 reasons:

    1. First of all what you say is not true: differences like 4.200000 versus 4.199999 will happen in the same form for doubles, the only “advantage” is that you have a longer list of 9 digits for doubles, but the first difference still occurs at the first decimal place.

    2. OK, I know that is not what you meant. You mean, that the difference is smaller, so if I subtract the two results than I get 0 for at least 6 decimal places. However, that does not “solve” the problem either. The problem is not how accurate the value is, but whether or not a fundamental algebraic identity is broken. ANY difference, however small in a+b+c!=c+a+b can lead to horrible consequences. For example, if you are sorting partial solutions and want to keep the “best 100″ for further processing, then you selection set will be different depending on what order you added up numbers during calculations. Same problem if you use any kind of threshold decisions, e.g. “solutions under a given energy limit are accepted” — again you get different solution sets that may even be radically different in terms of conformation — once you drop a solution it is simply not present and there is no guarantee that you have another similar solution that made it through the threshold condition.

    ZZ

Leave a Reply