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.