3

So, I have some R code that I am turning into C++. It reads a file, parses characters, and calculates tons and tons of means and standard deviations and returns them, along with counts of how many of each character occurred.

Now, there is a slight difference in the decimal values of the results R produces and those of C++. In the matrix of counts, as they are ints, the numbers are exactly the same. However, in the matrix of means, the values are the same up to the hundredths place, and they differ beyond that. With the standard deviations matrix, the values differ even more—to the tenths place.

What is causing this? I am assuming that there's some sort of precision difference in the ways that R and C++ handle numbers with decimals. I know that computers aren't exactly the best at representing floating point numbers to begin with, but how do I tell which output is the better?

...One thing I've tried is performing the calculation sqrt(41111.5/4522) in R, C++, and Calculator in Windows 7. They all produce the same result. Why then, when this exact same calculation is encountered during runtime, do they differ? In the runtime output, C++ agrees with Calculator, and R is the odd one out. I've also noticed that when performing these massive amounts of calculations, later output varies a slight bit more than earlier output. Does R just get tired when doing so many calculations and start to mess up? What's the deal?

Here are the outputs for the means:

C++:

38.6068 39.0122 38.633 38.5914 0
38.6159 38.7874 38.5053 38.7195 0
38.5205 38.7352 38.3694 38.5388 0
38.6331 38.7408 38.4588 38.5283 0
38.7503 38.6933 38.4173 38.6808 0
38.7637 38.7978 38.4967 38.603 0
38.7616 38.7384 38.4728 38.6946 0
38.6227 38.7689 38.4016 38.5352 0
38.5993 38.7334 38.3206 38.5514 0
38.6395 38.6598 38.43 38.4887 0
38.6414 38.746 38.4353 38.4908 0
38.4353 38.6767 38.3158 38.4694 0
38.35 38.5801 38.1486 38.3528 0
38.4122 38.6267 38.1731 38.3447 0
38.3751 38.5353 38.1782 38.2229 0
38.3373 38.6117 37.8952 38.2017 4.12443
38.332 38.4991 38.027 38.1984 0
38.2005 38.4417 38.0192 38.0446 4.12443
38.1719 38.4435 37.9727 38.0385 0
38.1346 38.3878 37.8634 37.9746 0
37.8505 38.2289 37.6202 37.6986 0
38.0932 38.142 37.7865 37.815 4.12443
37.9176 38.1381 37.5577 37.7273 0
37.7346 38.0934 37.4874 37.6546 0
37.6961 37.897 37.3342 37.4844 0
37.5534 37.9234 37.3341 37.3369 0
37.4914 37.7409 37.094 37.3211 0
37.2179 37.6653 36.9031 37.2592 0
37.0682 37.5625 36.6972 37.0218 4.12443
36.9713 37.4819 36.5387 36.8767 4.12443
36.8284 37.2411 36.223 36.6869 4.12443
36.7396 36.9682 36.0171 36.4556 4.12443
36.7874 36.9482 36.1641 36.5667 4.12443
36.695 36.9307 36.1856 36.3638 0
36.7224 36.9455 36.2212 36.695 4.12443
36.8983 37.1286 36.2652 36.8055 0
36.7835 36.8905 35.9562 36.4745 0
36.5364 36.9037 36.0927 36.4888 0
36.3959 36.6637 35.7378 36.323 0
35.9372 36.2034 35.452 35.6974 0

R:

            A        C        G        T N
[1,] 38.60573 39.01141 38.63195 38.59036 0
[2,] 38.61464 38.78523 38.50391 38.71826 0
[3,] 38.51908 38.73228 38.36774 38.53731 0
[4,] 38.63182 38.73834 38.45730 38.52657 0
[5,] 38.74903 38.69083 38.41585 38.67933 0
[6,] 38.76250 38.79534 38.49556 38.60156 0
[7,] 38.76039 38.73632 38.47145 38.69319 0
[8,] 38.62123 38.76703 38.40030 38.53354 0
[9,] 38.59810 38.73163 38.31917 38.55015 0
[10,] 38.63819 38.65792 38.42873 38.48740 0
[11,] 38.64002 38.74333 38.43387 38.48920 0
[12,] 38.43359 38.67401 38.31414 38.46783 0
[13,] 38.34827 38.57804 38.14686 38.35125 0
[14,] 38.41038 38.62463 38.17138 38.34302 0
[15,] 38.37329 38.53267 38.17653 38.22097 0
[16,] 38.33555 38.60949 37.89278 38.19956 4
[17,] 38.33024 38.49720 38.02496 38.19627 0
[18,] 38.19842 38.43880 38.01730 38.04205 4
[19,] 38.16998 38.44113 37.97058 38.03598 0
[20,] 38.13242 38.38488 37.86108 37.97245 0
[21,] 37.84771 38.22579 37.61745 37.69546 0
[22,] 38.09113 38.13806 37.78409 37.81250 4
[23,] 37.91487 38.13428 37.55473 37.72422 0
[24,] 37.73137 38.09007 37.48473 37.65181 0
[25,] 37.69295 37.89276 37.33098 37.48131 0
[26,] 37.54974 37.91984 37.33063 37.33263 0
[27,] 37.48773 37.73676 37.09027 37.31701 0
[28,] 37.21365 37.66051 36.89896 37.25519 0
[29,] 37.06418 37.55768 36.69254 37.01714 4
[30,] 36.96674 37.47745 36.53390 36.87150 4
[31,] 36.82324 37.23622 36.21721 36.68085 4
[32,] 36.73433 36.96207 36.01076 36.44930 4
[33,] 36.78201 36.94274 36.15842 36.56135 4
[34,] 36.68991 36.92524 36.17984 36.35769 0
[35,] 36.71720 36.94031 36.21548 36.68985 4
[36,] 36.89332 37.12322 36.25921 36.80057 0
[37,] 36.77870 36.88471 35.94958 36.46900 0
[38,] 36.53080 36.89801 36.08650 36.48348 0
[39,] 36.38996 36.65730 35.73058 36.31767 0
[40,] 35.93152 36.19707 35.44496 35.69141 0
11
  • 2
    it may help to post a little detail about the fp spec's/options of the C++ compiler and some minimal source code replicating the problem.
    – Necrolis
    Oct 28, 2011 at 19:06
  • By how much do they differ? Can you post the results of the calculations?
    – NickLH
    Oct 28, 2011 at 19:20
  • A simple calculation gives different results in R, C++ and Calculator, but only in the middle of a big run? Which one changes its behavior, and by how much?
    – Beta
    Oct 28, 2011 at 19:24
  • 3
    Typically, C implementations have two distinct floating point types: float and double, where double is more precise. I know only very little about R details, but it almost certainly uses one of those types itself (possibly both). If your C++ code uses float a lot, the R is probably more accurate, and you can change this by using double more. The same is true in reverse if the C++ code uses double. Oct 28, 2011 at 19:39
  • 2
    ...well, How do you do that in R then? R might use a different encoding, treat white space differently etc...
    – Tommy
    Oct 29, 2011 at 0:46

3 Answers 3

3

As you're no doubt aware, any number of things could be going on. As such, I can only posit one exotic possible cause.

One possibility is that R is performing computations in such a way as to minimize floating-point error; you wouldn't necessarily do this in C++ or when computing by hand, unless you knew better. In particular, you should sort your values in increasing order of exponent before you compute an aggregate sum (which should be the first step of any accurate averaging procedure). The reason for this is that floating-point arithmetic is not associative (unless you're using arbitary-precision libraries, which I assume is not the case). Due to rounding, (a + b) + c can equal c if a >> b, c whereas a + (b + c) would give a result larger than a (assuming a, b, c > 0). This is especially possible if R e.g. parallelizes its work, in which case you can reasonably expect to get a slightly different result every time!

Other less exotic possibilities include the following: the R and C++ code differ in some subtle though meaningful way (maybe there's an error in one where it misses the 72nd element, or you compute the STDEV in one using n-1 and using n in the other, etc.); there's a difference in the runtime between R and C++ which fundamentally results in this difference (different precision - doubles vs. floats vs. long doubles, etc., different library implementations, etc.).

I cannot tell whether this applies to your problem or not, but if nothing else somebody else might find this useful if they're having trouble understanding why floating point operations are not giving consistent results.

3
  • I doubt very much that R minimizes floating-point errors for the mean calculation by sorting. Calling sort on a large vector takes much longer than calling mean. And it currently (R 2.13.2) does not parallelize these operations.
    – Tommy
    Oct 28, 2011 at 20:29
  • 1
    @Tommy: While sorting could represent a significant overhead, I think that the underlying message - that R could be doing computations in a different order - remains a valid possibility. As I mention, parallelization can result in this. Also, depending on how R sorts the vector, the overhead could be pretty tolerable. I wouldn't put it past R to do some sophisticated processing logic. Then again, if R is doing anything wonky like this, it should be documented somewhere...?
    – Patrick87
    Oct 28, 2011 at 20:35
  • 2
    From R's source, summary.c:385, R does (paraphrasing) long double s = 0., t=0.; for (i = 0; i < n; i++) s += x[i]; s /= n; and then if (R_FINITE((double)s)) { for (i = 0; i < n; i++) t += x[i] - s; s += t/n; }. Oct 29, 2011 at 0:02
2

A few things I would check:

  • It may be a difference between using float and double or even long double. If you're using double it may be that R uses float. Try using float if its easy to switch between the two.
  • Check the floating point precision mode set in your compiler (for example, VC2010) and try different settings.
  • Ensure that all your calculations in C++ properly cast to double/float. For example, this code:

     double Test = 1.0 + 3/2;
    

    results in 2 not 2.5. R may cast such expressions differently leading to the differences seen in your results.

  • Double check that the functions in R and C++ are identical. For example, maybe cos() in R expects degrees while it is radians in C++. If in doubt make a quick test application in both to confirm.
  • If all else fails take one specific calculation and log/output detailed diagnostics for it in both the R and C++ applications. At some point you should begin to see the difference and trace back to where it originates. Try it with a smaller sample and see if you can replicate the behavior with just 6/60 samples instead of 6000.

One thing I just noticed is that in the last column of results C++ gives 4.12443 while R gives 4. Unless this is just a display issue see why this is the case. It may be that something in R is being rounded/cast to an integer but not in C++.

-1

All right, I'm just going to use the C++ values, mostly thanks to uesp's insight on that last column—R is probably casting to an integer during some intermediate step and losing some precision. I use doubles in every step in C++, so I trust it more (...not to mention that I just naturally mistrust interpreted languages a bit to begin with, hah)

1
  • 3
    R casting to an integer at an intermediate step seems really unlikely. I'm not convinced that you really know what's going on here ... If I were you, that would worry me. However, the problem seems sufficiently involved that debugging by committee, and by remote control, is just not that likely to succeed ...
    – Ben Bolker
    Oct 31, 2011 at 15:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.