5

I've encountered a weird phenomenon when I compute the dot product of matrices using numpy. In short, for mathematically identical computations, numpy is giving me different results.

The following code snippet illustrates the issue:

import numpy as np
np.set_printoptions(precision=10)
np.random.seed(2)  # for reproducibility

# create matrices:
A = np.random.rand(300, 1)
A = A/np.sum(A)
A = np.repeat(A, 15, 1)
B = np.random.rand(300, 300)

# convert data to float32:
A_star = A.astype("float32")
B_star = B.astype("float32")

# do the matrix multiplication A'BA and take the diagonal:
res_star = np.diag(A_star.transpose().dot(B_star.dot(A_star)))

# print the results:
print(res_star)

Running this in python3.5 with numpy1.11.1 on a Windows machine prints the following array:

[0.5000683069  0.5000683069  0.5000683069  0.5000683069  0.5000683069
 0.5000683069  0.5000683069  0.5000683069  0.5000683069  0.5000683069
 0.5000683069  0.5000683069  0.5000681877  0.5000681877  0.5000683069]

Please note that the values in res_star[12:14] are different from the other elements of the array - although mathematically one would expect them to coincide.

My own inquiries into these differences weren't particularly sucessful, but I think I narrowed it down a little:

  • with dtype float64 the values in the resulting array are identical
  • in case of float32 values, differences in the values typically only start at the 7th decimal place (the region where float32 precision ends)

Yet, to me these findings can't really explain the differences within the array since the same mathematical operation is taking place (they do explain differences in the values if I were to compute the same dot products on the same matrices with dtype float64). Also, for different seeds and matrix sizes the resulting values might actually coincide - so this behavior is not consistent.

  • on a UNIX machine the code snippet given above returns identical figures (at least in this particular configuration)

Now,

1) what is causing this behavior and

2) how can I ensure identical results within the array without changing the dtype to float64?

apitsch
  • 1,532
  • 14
  • 31
  • I can't say with certainty yet I'd put general blame on the fact that floating point arithmetic is neither associative nor commutative. Expressed in binary, the float32-values of 0.5000683069 and 0.5000681877 differ by exactly one bit, and I would suspect that numpy executes the steps you describe in a slightly different order for some entries, causing rounding to happen in slightly different places (e.g. SIMD-path for most but not all elements). – user2722968 Apr 25 '19 at 12:06
  • @user2722968 Thanks for giving me this hint. Can I somehow direct numpy to change the order of execution steps? Or is there any other way to ensure identical results? – apitsch Apr 25 '19 at 14:46
  • I can't tell... – user2722968 Apr 25 '19 at 15:08
  • `1`: might be a good idea to review how [`computers do math`](https://www.youtube.com/watch?v=pQs_wx8eoQ8)... TLDW it's all in binary scientific notation; there will be rounding. `2`: I don't believe you can with Python alone, and there are so many processor types (all with _little qwerks_), that it may not be worth your time to sort out in a language that does provide that level of control... Eg. ARM has had soft and hard floating point CPUs, each generation has been slightly different, and that's just one architecture _family_; which is _messy_ without changing precision from 64 to 32 too. – S0AndS0 Apr 25 '19 at 15:46
  • 32-bit floats generally only provide 7 digits of precision. If you require more precision than that, you have to use 64-bit floats, or reconfigure your problem to use integers. See this for details about floats: https://stackoverflow.com/questions/13542944/how-many-significant-digits-have-floats-and-doubles-in-java (not specific to java btw) So short answer is, unless you don't mind rounding the numbers to 5 decimals, you can't get that kind of guarentee. – Nicolaj Rasmussen Apr 26 '19 at 13:33

0 Answers0