0

I'm trying to move some old code I have from Matlab to Python and I'm having an issue with the way data is calculated by each. As a simple example, I have a stack of 2000 2D images called X. I want to calculate the mean of the stack and then subtract that from every image in the stack. In Python I have the following:

X = np.double(X)
mean = np.mean(X)
delta = np.zeros(X.shape)
for i in range(2000):
    delta[:, :, i] = X[:, :, i] - mean
print(np.mean(delta), np.sum(delta), np.median(delta), np.std(delta, ddof=1))

The mean of delta is -5.900437149625761e-16, the sum is -1.2038071872666478e-08, the median is -0.18050000000005184, and the standard deviation is 31.64658195121922.

However, if I use the following approach in Matlab I get something different:

X = double(X);
delta = zeros(size(X));
m = mean(X, 3);
for i = 1:2000
    delta(:, :, i) = X(:, :, i) -m;
end

Here, the mean is -7.478523241320775e-16, the sum is -1.525768311694264e-08, the median is -0.180500000000052, and the standard deviation is 31.646581951219204.

I am aware that these differences are very small but such differences can still have a slight impact in my results. Considering the simplicity of the calculation I'm surprised that they are different at all. What can be causing this to happen?

Felipe Moser
  • 323
  • 4
  • 19
  • I don't know if that will fix your issue, but for performance concerns you should try `delta = X - mean` directly. Also, shouldn't we have `np.mean(delta) == 0` by definition? – Guimoute Jul 28 '20 at 12:01
  • 1
    Those are the same numbers. Floating point accuracy simply stops you to getting better values, simply adding things in an oposite order can cause this. If your applications needs this level of accuracy, then you can not use floating point arithmetic. Anything below 1e-15 *is* zero, for example. – Ander Biguri Jul 28 '20 at 12:11
  • 2
    What I mean about the ordering is that floating point math is not associative. This means that `(a+b)+c != a+(b+c)`, so simply by indexing things in a different order, you are bound to get the differences you see. And yes, both are "correct", or as correct as doing maths with floating points can get. Try `0.1+(0.2+0.3) == (0.1+0.2)+0.3` in MATLAB. – Ander Biguri Jul 28 '20 at 12:15
  • I am closing you question as duplicate, because while the question itself is not the same, the answer is exactly what you need: understanding of floating point maths. In essence, you are seing "how computers work", there is no solution. If your math depends on such tiny variability, your math is wrong, or you need to do this the very slow symbolic way (likely the first). – Ander Biguri Jul 28 '20 at 12:19
  • @Guimoute I have, but I wrote it like this for the example – Felipe Moser Jul 28 '20 at 12:42

0 Answers0