I'm passing various bits of data to a function that computes variance along the first dimension. Sometimes the variance is zero, which is fine, but then the following strange thing happens:
>> sigma = data.var(axis=0) + 1e-7 # data has zero variance so all entries should equal 1e-7
>> sigma
array([1.e-07, 1.e-07, 1.e-07, ..., 1.e-07, 1.e-07, 1.e-07], dtype=float16)
>> (sigma==1e-7).all()
True
>> sigma[0]==1e-7
False
On its own, the fourth line would be explained by the 16-bit precision, and indeed
>> np.float16(1e-7)==1e-7
False
But it seems to contradict the third line, which says they are equal. This was causing a bug in my code. I can redesign around it, but I want to understand why numpy is doing this so I'm not caught out again in the future.
This comes from the fact that numpy type promotion treats scalars and arrays differently. You can see this with np.result_type:
>>> np.result_type(sigma, 1E-7)
dtype('float16')
>>> np.result_type(sigma[0], 1E-7)
dtype('float64')
Essentially, when an array value is compared to a scalar value (the first case), the dtype of the array value takes precedence. When comparing two scalars or two arrays (the second case), the highest precision takes precedence.
What this means is that when you evaluate (sigma == 1E-7), both sides are first cast to float16 before comparison, whereas when you evaluate sigma[0] == 1E-7, both sides are first cast to float64 before comparison.
Because float16 cannot perfectly represent the value 1E-7, this causes a discrepency in the scalar-comparison case, where both values are cast to float64:
>>> np.float16(1E-7).astype(np.float64)
1.1920928955078125e-07
>>> np.float64(1E-7)
1e-07
Finally, please note that these scalar-specific type casting rules are being changed in NumPy 2.0 (see NEP 50: Promotion Rules for Scalars), so if you run your code with NumPy 2.0, both cases will promote to float16 and return True.