1

I am given the two-dimensional arrays A and B. They are identical, but obtained with two different methods. Consider the following lines:

In  [1]: (A==B).all()
Out [1]: True
In  [2]: A.shape
Out [2]: (500, 10805)
In  [3]: B.shape
Out [3]: (500, 10805)
In  [4]: numpy.mean(A,axis=1)[0]
Out [4]: -0.006108739586784807
In  [5]: numpy.mean(A[0,:])
Out [5]: -0.006108739586784786
In  [6]: numpy.mean(B,axis=1)[0]
Out [6]: -0.006108739586784786
In  [7]: numpy.mean(B[0,:])
Out [7]: -0.006108739586784786

As you can see, the result from line [4] differs from the results from lines [5], [6], and [7], but they should be identical. What is the reason for this?

The same problem occurs with numpy.sum() and numpy.std().

Háski
  • 41
  • 2
  • 1
    Probably floating arithmetic precision – Dani Mesejo Oct 18 '21 at 13:13
  • related https://stackoverflow.com/questions/39757559/working-with-floating-point-numpy-arrays-for-comparison-and-related-operations – Dani Mesejo Oct 18 '21 at 13:16
  • Most likely numpy.mean does multiple thread operation (OpenMP) and the reduction operation is not deterministic. In floating operation A+B+C is not necessarily strictly equal to A+C+B. You could try setting KMP_DETERMINISTIC_REDUCTION=TRUE as global environment variable before launching your python code if you judge it is an issue. – PilouPili Oct 18 '21 at 13:24
  • Is the output deterministic, i.e. does re-executing line 4 always get the *same* "wrong" result? – MisterMiyagi Oct 18 '21 at 13:26
  • There have been several similar scattered stackoverflow posts and github issues related to numpy float rounding with sum and mean over the years- see [this issue, which includes links to others](https://github.com/numpy/numpy/issues/9393). The wording used in the docs is that higher precision addition *may* be done if sums are over the 'fast' axis in memory. I also asked a [related question some time ago](https://stackoverflow.com/questions/69278758/are-np-sum-and-np-add-reduce-actually-equivalent). It does seem to be based on threading behavior – kcsquared Oct 18 '21 at 13:36
  • @MisterMiyagi Yes, the output appears deterministc, i.e., re-executing line [4] always produces the same "wrong" result. – Háski Oct 18 '21 at 13:48
  • 1
    @DaniMesejo You are probably right. But how can A and B be equal but not produce the same "wrong" result? I.e, numpy.mean(A,axis=1)[0] != numpy.mean(B,axis=1)[0] – Háski Oct 18 '21 at 13:49
  • It seems totally fine to me: there are 14 digits equal (once rounded). This is close to what you should get with a pair-wise algorithm using double-precision. You should not expect Numpy to compute the result with no ULP error. If the two algorithm are slighly different, they will result in different results. Using *harware SIMD units* do such a thing. If you want a very very accurate result, you need to use a specific algorithm like the Kahan-summation (or even better: the Kahan-Babushka-Klein one) but such algorithm are much slower. *Your algorithm must not amplify the sum error*. – Jérôme Richard Oct 18 '21 at 14:14

1 Answers1

0

You're doing numerical operations. They are expected to be not 100% accurate. E.g. numerical addition is strictly speaking not associative. I.e. it depends on order of the additions. Consider this example

(1e-20+1.0)-1.0 -> 0
1e-20+(1.0-1.0) -> 1e-20

Now if the data is ordered in memory such that it is faster to sum it up in a different order you can get slightly different results.

Here I transpose (mirror along the diagonal), copy the data and transpose it back so I know it must be equal but get different means.

np.random.seed(20881)
A = np.random.random((100,100))
B = A.T.copy().T
(A==B).all() --> True
B.mean()-A.mean() -> 3.3306690738754696e-16

Notice that np.random.random((100,100)) is usually random but after np.random.seed(20881) will always return the same array. The number/array is chosen such that even the 100 by 100 array gives you an error as large as you got with your much larger array.

Lukas S
  • 3,212
  • 2
  • 13
  • 25
  • Thank you for your answer. Your post could explain why there is a difference between the means of A and B; indeed, in my program, I have an original matrix A which I augment to matrix B with a series of transposes. However, I don't think your post answers the following question: why there is a difference in computing the mean of the first column of A via lines [4] and [5] when there is no difference in computing the mean of the first column of B via lines [6] and [7]? – Háski Oct 28 '21 at 09:01
  • @Háski It would be easier to tell if you provided the numbers in A and B and the result of `np.info(A),np.info(B)`. – Lukas S Oct 28 '21 at 14:10