# Understanding HyperLogLog

HyperLogLog is an algorithm for estimating cardinality of multisets. Notably, the algorithm only uses \(O(1)\) space. Join me on a journey to understand it by implementing it!

Estimating the number of unique elements in a set/stream has many applications,
especially in the field of databases. For example, if getting such estimates is
faster or cheaper than computing the true value, it can be used as part of
join ordering where the cost of tables is derived from the estimated row count.
This paper
describes the `HyperLogLog`

algorithm - an algorithm that only uses \(O(1)\)
space, and is only \(O(n)\) in runtime. It’s a very cool algorithm, and I
recommend reading the original paper, though there were definitely aspects of it
that I struggled to get an intuition for, so I decided to implement it myself as
a learning exercise.

## The HyperLogLog Algorithm

I highly recommend reading this
write-up
by Meta about how `HyperLogLog`

works. It helped me get a lot of intuition for
the process, but here’s my short summary:

`HyperLogLog`

works by mapping elements into a roughly uniformly distributed space
by using a hash function. If this space was between 0 and 1, we could use the
inverse of the lowest element as an approximation of the number of elements (if
you have \(n\) elements evenly distributed in \((0, 1)\), then the lowest
element will be \(\frac{1}{n}\)). We can get a value in \((0, 1)\) by
counting the number of leading 0’s (\(k\)) in every hashed value and then
using \(\frac{1}{2^{k}}\) as our value - we only care about the lowest value,
so we look for the element that maximizes \(k\). However, this would only give us
coarse, power-of-two sized estimates, so we could do better by using many
different hash functions and averaging the values estimated by each function.
This is computationally expensive, so instead, we take some number of bit
\(b\) from the prefix of the hashed value, and use that as an identifier to
split the input into \(m = 2^{b}\) streams. The remaining bits are used as
before to estimate the cardinality.

One interesting property to note is that `HyperLogLog`

“remembers” seen elements
by virtue of it’s usage of `max`

. The key idea is that `max`

is commutative, and
\(max(x, x) = x\), which means that adding the same element twice is the same
as adding it once. Since the commutativity means that order of operations
doesn’t matter, the order in which repeated elements are present doesn’t matter
either. This is mainly relevant to this discussion in terms of testing, since it
simplifies the kind of test data I’ll need to produce.

## Writing my implementation

Once I’d read the initial paper, I set about trying to implement it myself. The
algorithm itself is very straightforward, but was even cooler was that Github’s
Copilot seemed to understand what I was doing and suggested almost exactly what
I wanted! However, I was a little weary of relying on the auto-complete too
heavily since I wanted the experience of writing things out myself to learn
better. For my choice of hash function, I went with python’s built in `hash`

.
Since most of the paper relied on manipulating bits, the actual function I used
was: `bin(abs(hash(x)))[2:]`

which takes the absolute value of the hash,
converts it to a binary string, and removes the `0b`

prefix. However, running
my initial code produced wildly inaccurate results. Commence debugging!

## Debugging my implementation

My first step in the debugging process was to find a reference implementation to compare against. I quickly found this repo and began comparing. This repo seems to include a lot of optimizations that the original paper doesn’t mention, so I tried ignoring most of those, since this exercise is mostly about getting the core ideas down. Scrolling through the code seemed to indicate that I was generally on the right track, which was encouraging.

The most obvious difference to me was that the reference implementation tried to ignore registers that didn’t have any elements (empty streams ignored), however, that didn’t affect this implementation since every register was being filled. I added a debug assert to ensure that I wasn’t hitting that cast to be safe.

My first suspicion was the choice of hash function. At work we’d seen strange
behavior with the distribution of python’s hash function, albeit only on older
python versions. Looking at the reference implementation, they used `sha1`

instead. Could that really make a difference in the computation? To verify this,
I injected code into the reference and my own implementation to produce a
histogram of results from calls to the `rho`

function, which was responsible for
computing length of the 0 prefix of hashed data.

Here’s what my `rho`

function looked like with this additional collection:

```
class HyperLogLog:
...
observed_ps = {}
def rho(self, s: str) -> int:
v = s.index('1') + 1
self.observed_ps[v] = self.observed_ps.get(v, 0) + 1
return v
...
h = HyperLogLog()
# Simulate 100_000 unique values
for x in range(100000):
h.hyperloglog(f"foobar{x}")
print({k: h.observed_ps[k] for k in sorted(h.observed_ps.keys())})
```

And here’s some sample values produced:

```
{
1: 50102,
2: 25102,
3: 12383,
4: 6265,
5: 3093,
6: 1465,
7: 779,
8: 387,
9: 210,
10: 105,
11: 53,
12: 37,
13: 9,
14: 3,
15: 4,
16: 2,
18: 1
}
```

Similar modifications were done to the reference implementation. This turned out to be a great idea, and confirmed that my hash function was producing values that closely matched the distribution in the reference implementation. I was also concerned that the bit width of values I was working with not being large enough, but the histogram seemed to show that having at least 20 bits was more than sufficient for the error rate desired. This makes sense since the probability of encountering values with a specific 0-prefix length is exponentially inversely related to the length.

After some more poking around and comparing intermediate values from my solution vs the reference, I stumbled upon the bug. For some reason, copilot had implemented the sum term in the harmonic mean as \(2^{\sum{2^{-M_j}}}\) instead of \(\sum{2^{-M_j}}\). It goes to show that trusting AIs in domains where I’m not already well-informed is a bad idea. Lessons learned.

While this brought the estimates to a more reasonable value, it still seemed to
be underestimating by a factor of 2. This seemed suspiciously consistent and the
algorithm already includes a constant factor adjustment, so I concluded that the
general steps of the algorithm must be mostly correct. After reviewing the
reference solution again, I realized that my implementation of `rho`

was off by
1, which made sense, and it was easy to prove by looking at the original
equation that this would cause the estimate to be divided by 2 overall.

The working code is presented below:

## My code

Here’s the final code I ended up with, including some comments.

```
import math
from typing import List
class HyperLogLog:
error_rate = 0.1
# Determine the number of "buckets" based on the error_rate
b = int(math.ceil(math.log((1.04 / error_rate) ** 2, 2)))
m = 2**b
# Registers storing Max-values for each "bucket"
M: List[float]
def __init__(self) -> None:
# Initialize registers to 0 (-inf in the paper)
self.M = [0] * self.m
def rho(self, s: str) -> int:
"""Find the first 1 bit"""
v = s.index('1') + 1
return v
def h(self, s: str) -> str:
"""Compute the hash of s and return it as a binary string"""
x = abs(hash(s))
s = bin(x)[2:]
# pad to 32 bits
while len(s) < 32:
s = '0' + s
# truncate to 32 bits
s = s[-32:]
return s
def add(self, s: str):
"""Add an element to the estimator"""
x = self.h(s)
j = int(x[:self.b], 2)
w = x[self.b:]
self.M[j] = max(self.M[j], self.rho(w))
def z(self) -> float:
"""The sum portion (denominator) of the harmonic mean"""
return sum(2**-x for x in self.M)
def get_alpha(self):
"""Constant factor to weight results"""
# This is copied from https://github.com/svpcom/hyperloglog
if self.b == 4:
return 0.673
if self.b == 5:
return 0.697
if self.b == 6:
return 0.709
return 0.7213 / (1.0 + 1.079 / self.m)
def estimate(self):
"""Estimate the cardinality of all elements seen so far"""
# constant * m * harmonic mean of 2**-M
return self.get_alpha() * (self.m**2) / self.z()
h = HyperLogLog()
for x in range(100000):
h.add(f"foobar{x}")
# The implementation above doesn't handle unvisited registers well
# assert h.M.count(0) == 0, f"num registers/0 regs = {len(h.M)}/{h.M.count(0)}"
print("Estimate:", h.estimate())
```

## Playing with HyperLogLog

Once I had a working implementation, I was able to experiment with a few
things. I tried using `k`

hash functions instead of partition the input into `k`

streams. This was done by using the following logic:

```
for i in range(len(M)):
M[i] = max(M[i], rho(hash(f'hash_{i}: {s}')))
```

This performed worse (and was insanely slow), tending to overestime far more than the original implementation. I want to continue exploring this, and see if maybe it’s just a matter of adjusting the constants involved.

I also tried taking inspiration from older algorithms like `LogLog`

and I wanted
to see what would happen if I only took the 70% lowest values instead of all
values for the harmonic mean (bias away from outliers), this was done by
changing the `z`

function to:

```
def z(self) -> float:
def getLower70Percent():
return sorted(self.M)[:int(self.m * 0.7)]
return sum(2**-x for x in getLower70Percent())
```

and multiplying by \(0.7\) in the final computation. This performed way worse than the original implementation, but if I removed the \(0.7\) multiplied at the end, it only slightly underestimated compared to the original.

## Thoughts

This was a fun exercise, and implementing my own version of this algorithm greatly improved my understanding of the concepts involved. There’s definitely parts of the proof that invoke concepts from statistics that I’m not super familiar with anymore, but brushing up on some of those was also a good experience. I’m glad I got something working, and being able to tweak the algorithm and see how things change was also interesting. I would definitely recommend this process to anyone exploring a new area!