Welford's online variance calculation algorithm for vectors

Question

I'm developing a face recognizing application using the face_recognition Python library.

The faces are encoded as 128-dimension floating-point vectors. In addition to this, each named known person has a variance value, which is refined iteratively for each new shot of the face along with the mean vector. I took the refining formula from Wikipedia.

I'm getting some false positives with the recognized faces, which I presume is because the library was developed primarily for Western faces whereas my intended audience are primarily Southern-Eastern Asian. So my primary concern with my code, is about whether or not I had gotten the mathematics correct.

Here's my refining algorithm in Python

import sys 
from functools import reduce
from math import hypot

# fake testing data. 

# new reference face. 
refenc = (0.2, 0.25, 0.4, 0.5) * 32

# previous face encoding and auxiliary info. 
baseenc = (0.2, 0.3, 0.4, 0.5) * 32
v = 0.01         # variance
n = 3            # current iteration

n = min(n, 28)   # heuristically limited to 28. 

vnorm = lambda v: reduce(hypot, v)
vdiff = lambda u, v: list(map(lambda s,t:s-t, u, v))

delta1 = vdiff(refenc, baseenc)
if( vnorm(delta1) > 0.4375 and n > 1 ):
    sys.exit() # possibly selected wrong face.
    pass
newenc = [ baseenc[i] + delta1[i] / n for i in range(128) ]

delta2 = vdiff(delta1, newenc)
v = v*(n-1)/n + vnorm(delta1)*vnorm(delta2)/n


print(repr((newenc, v, n)))

Irrelevant note : I used struct.(un)pack to serialize in binary to save space, because the repr of the data is too big.

Answer 1

I can't say I understand exactly what your algorithm does. Even when looking at Wikipedia it is hard to see. So, why don't you just use the reference implementation from there? It has functions with nice names and everything (well, not quite PEP8 compatible, but that could be changed).

Alternatively, once you have vectors of numbers and want to perform fast calculations on them, you probably want to use numpy, which allows you to easily perform some operation on the whole vector, like taking the difference between two vectors, scaling a vector by a scalar, taking the sum, etc.

import numpy as np

MAGIC_CUTOFF = 0.4375

class WrongFaceException(Exception):
    pass

def norm(x):
    return np.sqrt((x**2).sum())

def update(base_encoding, new_face_encoding, v, n):
    delta1 = new_face_encoding - base_encoding
    if norm(delta1) > MAGIC_CUTOFF and n > 1:
        raise WrongFaceException(f"norm(delta1) = {norm(delta1)} > {MAGIC_CUTOFF}")
    new_encoding = base_encoding + delta1 / n
    delta2 = delta1 - new_encoding
    v = v * (n - 1) / n + norm(delta1) * norm(delta2) / n
    return new_encoding, v, n + 1

if __name__ == "__main__":
    new_face_encoding = np.array((0.2, 0.25, 0.4, 0.5) * 32)
    base_encoding = np.array((0.2, 0.3, 0.4, 0.5) * 32)
    v = 0.01
    n = 3
    
    updated_encoding, v, n = update(base_encoding, new_face_encoding, v, n)
    print(updated_encoding, v, n)

In addition to using numpy, I expanded some of the names so it is clearer what they are, I made the cutoff value a global constant (you should give it a meaningful name), added a custom exception instead of just dying (the pass after it was unneeded and unreachable BTW), which makes it more clear what happened, wrapped the calling code under a if __name__ == "__main__": guard to allow importing from this script without running it and followed Python's official style-guide, PEP8, by having operators surrounded by one space.

Instead of this norm function, you can also use scipy.linalg.norm, which might be even faster.

---

To answer if you used the math correctly, I would say partially yes. It looks like you implemented the explicit formula given for the sample variance:

\$\sigma^2_n = \frac{(n - 1)\sigma^2_{n-1} + (x_n - \bar{x}_{n-1})(x_n - \bar{x}_n)}{n}\$

However, note that Wikipedia also mentions that this formula suffers from numerical instabilities, since you "repeatedly subtract a small number from a big number which scales with \$n\$".

Instead you want to keep the sum of the squared differences \$M_{2,n} = \sum_{i=1}^n (x_i - \bar{x}_n)^2\$ around and update that:

\$M_{2,n} = M_{2,n-1} + (x_n - \bar{x}_{n-1})(x_n - \bar{x}_n)\$

with

\$\bar{x}_n = \bar{x}_{n -1} + \frac{x_n - \bar{x}_{n-1}}{n}\$

Then the new population variance is given simply by:

\$\sigma^2_n = \frac{M_{2,n}}{n - 1}\$

This is a bit tricky to implement, since your \$x_i\$ are actually vectors of numbers, so you would need to find out how to properly reduce it in dimension (in your previous formula you used the norm for that).