Mean vs mean

I decided to start learning about machine learning a couple of weeks ago. Partly because it sounds interesting, partly because I think I have enough math to handle it, and partly because I felt it was important to have a deeper knowledge if I want to criticize AI.

I started with a library copy of Friedman and Mendhekar’s The Little Learner but decided to drop it after I was about halfway through. Although I was following along reasonably well, I kept struggling with the book’s notation—after 200 pages, I shouldn’t have to keep returning to the front of the book to see what things like t|i mean. If the book’s Scheme code were actually written in Scheme, I might have stuck with it, but its mixture of Scheme with math was too frustrating.

Yesterday I started Kneusel’s Practical Deep Learning. My library gives me free access to O’Reilly ebooks, so it’s easy to have the book open on my iPad while I work through example code on my MacBook. And yes, I certainly feel more comfortable working in Python. But I was surprised by a result pretty early on.

It had to do with calculating the mean values of the columns of a 2D array. The book is explaining a way to normalize (or standardize) a set of data whose features are on vastly different scales. Here’s the code (in an interactive Python session) that does it using the book’s example data:

python:
>>> import numpy as np
>>> x = np.array([
...  [6998, 0.1361, 0.3408, 0.00007350, 78596048],
...  [6580, 0.4908, 3.0150, 0.00004484, 38462706],
...  [7563, 0.9349, 4.3465, 0.00001003,  6700340],
...  [8355, 0.6529, 2.1271, 0.00002966, 51430391],
...  [2393, 0.4605, 2.7561, 0.00003395, 27284192],
...  [9498, 0.0244, 2.7887, 0.00008880, 78543394],
...  [4030, 0.6467, 4.8231, 0.00000403, 19101443],
...  [5275, 0.3560, 0.0705, 0.00000899, 96029352],
...  [8094, 0.7979, 3.9897, 0.00006691,  7307156],
...  [ 843, 0.7892, 0.9804, 0.00005798, 10179751],
...  [1221, 0.9564, 2.3944, 0.00007815, 14241835],
...  [5879, 0.0329, 2.0085, 0.00009564, 34243070],
...  [ 923, 0.4159, 1.7821, 0.00002467, 52404615],
...  [5882, 0.0002, 1.5362, 0.00005066, 18728752],
...  [1796, 0.7247, 2.3190, 0.00001332, 96703562],
... ])
...
>>> xs = (x - x.mean(axis=0)) / x.std(axis=0)
>>> xs
array([[ 0.69308862, -1.12598999, -1.53184967,  0.95258185,  1.18240196],
       [ 0.54647372, -0.01203876,  0.50510857, -0.01928358, -0.1141859 ],
       [ 0.89126426,  1.38267719,  1.51932212, -1.19969655, -1.14033263],
       [ 1.16906091,  0.49704356, -0.17121154, -0.53403994,  0.3047611 ],
       [-0.92213055, -0.10719726,  0.30790249, -0.38856532, -0.47533014],
       [ 1.56997199, -1.47678885,  0.33273415,  1.47140742,  1.18070086],
       [-0.34794732,  0.47757218,  1.88235192, -1.40315756, -0.73969021],
       [ 0.0887406 , -0.43538419, -1.73773921, -1.23496312,  1.7456197 ],
       [ 1.07751429,  0.95242255,  1.24754488,  0.72911384, -1.12072823],
       [-1.46579825,  0.92509981, -1.04466154,  0.42629603, -1.0279233 ],
       [-1.33321348,  1.4501989 ,  0.03239288,  1.11026413, -0.89668956],
       [ 0.30059562, -1.45009422, -0.26155005,  1.70335297, -0.25050967],
       [-1.43773798, -0.24726556, -0.43400063, -0.70325168,  0.33623535],
       [ 0.30164788, -1.55279003, -0.62130451,  0.1780736 , -0.75173073],
       [-1.1315303 ,  0.72253467, -0.02503987, -1.08813209,  1.7674014 ]])

The example data consists of 15 samples, each of which has 5 features. The third command rescales each feature by subtracting its mean and dividing by its standard deviation.

In the book, Kneusel checks the mean of the fourth column of rescaled data to confirm that it’s basically zero. He gets −1.33×10−16, which is effectively zero for 64-bit floating point values.

To make sure I was doing things right, I ran this:

python:
>>> xs.mean(axis=0)
array([ 1.48029737e-17, -8.88178420e-17,  6.77698638e-17, -1.62832710e-16,
        4.44089210e-17])

While the fourth mean, −1.63×10−16, is certainly as good a zero as the book’s answer, it isn’t the same. To double-check my work, I ran this:

python:
>>> xs[:, 3].mean()
-1.3322676295501878e-16

This is the book’s value. What’s going on?

You may think this is perfectly normal. Floating point numbers are inherently inexact, and as long as I get a result that’s effectively zero, there’s nothing to be bothered about. But I was bothered.

I’m using the same function, mean, from the same library, NumPy, on the same set of data, and I get different answers depending on how I invoke that function. Or maybe I do get the same answers. Here’s another command and its output:

python:
>>> for i in range(xs.shape[1]):
...     print(f'Column {i}: {xs.mean(axis=0)[i] == xs[:, i].mean()}')
...
Column 0: True
Column 1: True
Column 2: True
Column 3: False
Column 4: False

For the first three columns, it doesn’t matter which way I use mean; for the last two, it does.

I should mention here that while the different results are obvious when the mean is close to zero, this behavior is present when operating on any 2D array. I’ve tried it out on several 2D arrays made from randomly generated floating point numbers.

During these experiments with randomly generated arrays, I got the sense that the likelihood of getting the same mean both ways got smaller as the number of rows in the array got larger. I decided to check that out with this script:

python:
 1:  #!/usr/bin/env python3
 2:  
 3:  import numpy as np
 4:  
 5:  rng = np.random.default_rng()
 6:  
 7:  N = 5000
 8:  
 9:  for n in (10, 100, 1000, 10_000, 100_000):
10:    count = 0
11:    for i in range(N):
12:      x0 = rng.uniform(low=0, high=100, size=n)
13:      x1 = rng.uniform(low=100, high=200, size=n)
14:      x = np.column_stack((x0, x1))
15:      if x.mean(axis=0)[0] == x[:, 0].mean():
16:        count += 1
17:      if x.mean(axis=0)[1] == x[:, 1].mean():
18:        count += 1
19:    print(f'n = {n:,d}')
20:    print(f'  Same means: {count/(2*N):.2%}')

The results were

n = 10
  Same means: 73.97%
n = 100
  Same means: 28.47%
n = 1,000
  Same means:  7.47%
n = 10,000
  Same means:  2.50%
n = 100,000
  Same means:  1.03%

So arrays with more rows are less likely to yield the same mean values when the mean function is used in these two different ways.

My working theory is that when I use mean() on a 1D array extracted from a 2D array, the computations are done differently than when I use mean(axis=0) on the 2D array. Mathematically the same, but computationally different. Maybe it’s as simple as the summation being done in a different order so the partial sums round off in slightly different ways. The more partial sums, the greater the likelihood of a different answer.

Ultimately, −1.33×10−16 and −1.63×10−16 are both zero, so this was mainly an exercise in satisfying my curiosity. I worry when computer results that ought to be exactly the same are only nearly the same. It usually means I’ve done something wrong, but apparently not in this case.