Small, unpredictable results in runs of a deterministic model

Question

I have a sizable model (~5000 lines) written in C. It is a serial program, with no random number generation anywhere. It makes use of the FFTW library for functions using FFT - I do not know the details of the FFTW implementation , but I assume the functions therein are deterministic also (correct me if I'm in err).

The issue I cannot understand is that I am getting small differences in the results for identical runs on the same machine (same compiler, same libraries).

I use double-precision variables, and to output the result in variable value for example, I issue: fprintf(outFID, "%.15e\n", value); or
fwrite(&value, 1, sizeof(double), outFID);

And I would constantly get differences such as:
2.078434696522064e-16 vs. 2.078434696522063e-16

I have spent much time trying to figure out why this is. I initially thought one of my memory chips have gone bad, and I've ordered and replaced them, to no avail. I subsequently also tried running my code on a colleague's Linux machine, and I get differences of the same nature.

What could be causing this? It is a small issue now, but I wonder if it is the "tip of the iceberg" (of a serious problem).

I thought I would post here instead of StackOverflow in case someone working with numerical models might have come across this issue. If anyone can shed light on this, I'd be much obliged.

Followup to comments:
Christian Clason and Vikram: first, thank you for your attention to my question. The articles you linked to suggest that: 1. rounding errors limits the accuracy, and 2. different code (such as introducing seemingly harmless print statements) can affect the results up to the machine epsilon. I should clarify that I am not comparing the effects fwrite and fprintf functions. I am using one OR the other. In particular, the same executable is used for both runs. I am simply stating the issue occurs whether I use fprintf OR fwrite.

So the code path (and executable) is the same, and the hardware is the same. With all these external factors held constant, where does the randomness come from, fundamentally? I suspected the bit flip happened due to faulty memory not retaining a bit correctly, which is why I replaced the memory chips, but that doesn't seem to be the issue here, I verified and you indicated. My program outputs thousands of these double precision numbers in a single run, and there are always a random handful that have random bit flips.

Followup to Christian Clason's first comment: Why is $2\cdot 10^{-16}$ the same as 0 within machine precision? The smallest positive number for a double is 2.22e-308, so shouldn't that be equal to 0? My program outputs thousands of values in the 10^-16 range (ranging from 1e-15 to 8e-17) and we've been seeing meaningful variations in our research project, so I'm hope we've not been looking at nonsensical numbers.

Followup #2:
This is a plot of the time series output by the model, to aid in the offshoot discussion in the comments.

Answer 1

There are aspects of modern computing systems that are inherently non-deterministic that can cause these kinds of differences. As long as the differences are very small in comparison with the required accuracy of your solutions, there probably isn't any reason to worry about this.

An example of what can go wrong based on my own experience. Consider the problem of computing the dot product of two vectors x and y.

$d=\sum_{i=1}^{n} x_{i}y_{i}$

Computing this dot product requires computing the products $x_{i}y_{i}$ and then adding up the results. The floating point multiplications should produce precisely the same results every time. If the floating point additions are computed in the same order each time then the sum should be identical. However, since floating point addition is not associative, you might get different results if the product of two vectors is computed in such a way that the additions are performed in different order.

For example, you might compute the product of two vectors first as

$d=((x_{1}y_{1})+(x_{2}y_{2}))+(x_{3}y_{3})$

and then as

$d=(x_{1}y_{1})+((x_{2}y_{2})+(x_{3}y_{3}))$.

How might this happen? Here are two possibilities.

1. Multithreaded computations on parallel cores. Modern computers typically have 2, 4, 8, or even more processor cores that can work in parallel. If your code is using parallel threads to compute a dot product on multiple processors, then any random perturbation of the system (e.g. the user moved his mouse and one of the processor cores has to process that mouse movement before returning to the dot product) could result in a change in the order of the additions.

2. ALignment of data and vector instructions. Modern Intel processors have a special set of instructions that can operate on (for example) for floating point numbers at a time. These vector instructions work best if the data is aligned on 16 byte boundaries. Typically, a dot product loop would break the data up into sections of 16 bytes (4 floats at a time.) If you rerun the code a second time the data might be aligned differently with the 16 byte blocks of memory so that the additions are performed in a different order, resulting in a different answer.

You can address point 1 by making your code run as a single thread and disabling all parallel processing. You can address point 2 by requiring memory allocation to align blocks of memory (typcially you'd do this by compiling the code with a switch such as -align.) If your code is still giving results that vary then there are other possibilities to look at.

This documentation from Intel discusses issues that can lead to non reproducibility of results with the Intel Math Kernel Library. Another document from Intel that discusses compiler switches to use with Intel's compilers.

Answer 2

The mentioned FFTW library might run in a non-deterministic mode.

If you're using FFTW_MEASURE or FFTW_PATIENT mode, the programs checks at runtime, which parameter values work the quickest and then will use those parameters throughout the whole program. Because the run time will obviously fluctuate a bit, the parameters will be different and the result of the Fourier transforms will be non-deterministic. If you want deterministic FFTW, use the mode FFTW_ESTIMATE.

Answer 3

While it's true that expression term evaluation order changes may very well occur due to multi-core/multi-thread processing scenarios, don't forget that there may be (even though it's a long shot) some kind of hardware design flaw at work. Remember the Pentium FDIV problem? (See https://en.wikipedia.org/wiki/Pentium_FDIV_bug ). Some time ago, I worked on pc-based analog circuit simulation sofware. Part of our methodology involved developing regression-test suites, which we would run against nightly builds of the software. With many of the models we developed, iterative methods (e.g. Newton-Raphson (https://en.wikipedia.org/wiki/Newton%27s_method) and Runge-Kutta) were used extensively in the simulation algorithms. With analog devices, it is often the case that internal artifacts, such as voltages, currents, etc., occur having extremely small numeric values. These values, as part of the simulation process, are incrementally varied over (simulated) time. The magnitude of these changes may be extremely small, and what we often observed was that subsequent FPU operations on such delta values bordered on the "noise" threshold of the FPU's precision (64-bit floating has a 53-bit mantissa, IIRC). That, coupled with the fact that we often had to introduce "PrintF" logging code into models to allow debugging (ah, the good ol' days!), practically guaranteed sporadic results, on a daily basis! So, what's all this mean? You have to expect to see differences under such circumstances, and the best thing to do is define and implement a way to decide (magnitude, frequency, trend etc) when/how to ignore them.

Answer 4

While floating point rounding from async operations may be the issue, I suspect that it is something more banal. The use of uninitialized variable that is adding randomness to your otherwise deterministic code. It is a common issue that is often overlooked by developers because when you run in debug mode all variables are initialized to 0 on declaration. When not running in debug mode, the memory assigned to a variable has whatever value the memory had before the assignment. The memory is not zeroed on assignment as an optimization. If this is happening in your code, it will be easy to fix, less so in the libraries code.