Small Python/NumPy snippet showcasing some of the inaccuracies of floating point calculation.

FloatPrecision.py

import numpy as np

# Creating some random data with 4 columns centered around 10
data = np.random.randn(4, 512) + 10
data = data.astype(np.float32) #np.float64
column_means = np.mean(data, axis=1, keepdims=True)
print(column_means.reshape(-1))
# [ 9.975992   9.9283495 10.003674  10.012635 ]

# Subtracting the column mean should yield zero-centered columns
data = data - column_means
column_means = np.mean(data, axis=1, keepdims=True)
print(column_means.reshape(-1))
# [-5.3644180e-07 -1.8626451e-08  2.7567148e-07  1.7695129e-07]
# In practice, the new columns' means are close to zero, but not truly zero.

# If we repeat the experiment with double precision (np.float64), we get:
# [ 9.96945238  9.93846489 10.04536377 10.05911843]
# [ 9.95731275e-16 -8.76035355e-16  6.73072709e-16 -6.60929644e-16]
# The results are still not zero but they are considerably closer.
# We went from e-7 to e-16. An improvement of 9 decimal digits.