2

I am multiplying two float64 matrices with the following values:

import numpy as np


# 4x5 matrix with identical columns.
x = np.zeros((4, 5,), dtype=np.float64)
x[1] = 1
x[3] = -3

w = np.array([1, 1, -1, 1 / 3], dtype=np.float64)

# The result should be an array of size 5 with equal values.
result = np.matmul(w, x)

print(x)

>>> 
[[ 0.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.]
 [ 0.  0.  0.  0.  0.]
 [-3. -3. -3. -3. -3.]]

print(w)

>>> [ 1.          1.         -1.          0.33333333]

print(result)

>>> [5.55111512e-17 5.55111512e-17 5.55111512e-17 5.55111512e-17 0.00000000e+00]

The result array should contain identical values, since each item is a dot product of the w array with an identical column. However, the last item is 0.0 unlike the other values which are very close to 0. This has a large effect over calculations downstream.

I am guessing this has something to do with the value 1/3, since replacing it with 1/2 gives a stable result. How can this instability be solved though?

Additional info since the problem doens't reproduce on all machines

I am using numpy 1.18.2 and Python 3.7.3, on MacOS. The problem reproduces on another machine which runs Ubuntu with the same Python and numpy versions.

Nur L
  • 791
  • 7
  • 15
  • I ran the same code on google colab and it gave expected value for `result`. – Balaji Sep 20 '21 at 15:18
  • 2
    While I think this is very interesting, though if this has a large effect downstream, then you actually have much bigger problems, and I would say your algorithm generally cannot be trusted. – lalala Sep 20 '21 at 15:21
  • This is by the nature of floating point numbers representation in computers (also why it is platform dependent). Related questions: [one](https://stackoverflow.com/questions/588004/is-floating-point-math-broken), [two](https://stackoverflow.com/questions/56820/round-doesnt-seem-to-be-rounding-properly) – Marat Sep 20 '21 at 15:54
  • Thank you @Marat. I don't fully understand though why this happens for only one item and not for all of the columns. – Nur L Sep 20 '21 at 16:17
  • 1
    @NurL bits in computer memory represent powers of 2. Some numbers, like 1/4, can be stored without loss of precision, some (1/3) have to be slightly rounded to fit this representation. Exact rounding error depends on the CPU and, sometimes, the OS (or, rather, the compiler it comes with). In some cases, these small rounding errors are enough to throw off the final results. Depending on the problem, there might be ways to counter that, e.g. by working with log-scaled values instead – Marat Sep 20 '21 at 16:33
  • 1
    @Marat: No, this is not by the nature of floating-point numbers. The nature of floating-point numbers does not cause identical calculations to produce different results. What must be happening is that `numpy` is not computing `matmul` in the simplistic by-definition way of doing a dot product of a row with a column. It must be doing something else that involves different calculations for the different positions, and that leads to different results. The same happens in integer arithmetic, as if somebody tried to calculate 7/3•3 with `7/3*3` versus `7*3/3`. – Eric Postpischil Sep 20 '21 at 19:04
  • On `numpy` 1.21.1 I don't get the near-zero values, even when I look at the `tolist()` form. `matmul` does pass the calculation on to BLAS kind of functions, so it's hard to identify exactly what algorithm is used. I suppose we could search the release-notes for any relevant change since version 1.18. – hpaulj Sep 20 '21 at 19:14
  • @EricPostpischil there's plenty of [anecdotic] evidence of [the opposite](https://www.google.com/search?q=floating+point+operations+determinism) – Marat Sep 20 '21 at 19:19
  • @Marat: Google results are not deterministic; they change for different users, locations, times, and so on, so you should not use them as citations. The first answer at one of the links in the results your link showed me states “Floating-point is deterministic.” As I wrote, it is not the **nature** of floating-point arithmetic that makes results non-deterministic. Performing operations in non-deterministic orders or groupings or with non-deterministic settings yields non-deterministic results… – Eric Postpischil Sep 20 '21 at 19:30
  • … That is truefor any kind of computation, such as integer arithmetic; it is not caused by the nature of floating-point arithmetic. – Eric Postpischil Sep 20 '21 at 19:31

1 Answers1

-1

Changed the array sizes and changed which rows are set to 1 and -3 & 1/3 on a i7 Mac running macOS 11.6 with python 3.9.7 and numpy 1.21.2 The 5.5e-17 values only occur when row 0 is 1 and row 2 is -3 and when row 1 is 1 and row 3 is -3. The number of 5.5e-17 and 0 values changes depending on the number of columns. Some number of columns such as 16 produce all 5.5e-17 values. Same behavior when -3 & 1/3 is replaced with values that are not a factor of 2 such as 12 & 1/12 though 5.5e-17 sometimes changes to a different e-17 value.

slasley
  • 1
  • 1
    As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers [in the help center](/help/how-to-ask). – Community Sep 21 '21 at 12:29