In this post we’ll study the hyper log log algorithm and provide an implementation in Rust.

### Introduction

Philippe Flajolet was a French computer scientist at INRIA [5]. He introduced the field of Analytic Combinatorics and is known for the creation of a family of algorithms of probabilistic counting, including the HyperLogLog.

HyperLogLog is a probabilistic algorithm to determine the number of distinct elements (or cardinality) from a multi-set (a set that allows repeated values) with high accuracy using very low memory. It is then suitable for streaming applications or distributed databases in which we do not have the luxury of keeping all distinct values in memory at any time.

### The HyperLogLog

The intuition behind the HyperLogLog is the following. Imagine we have a basket full of balls, each containing a number, which can be repeated. We want to estimate how many distinct numbers we see in the basket with the limitation that we can only look at one ball at a time and we do not have paper and pen, so we have to rely on memory.

The number of balls in the basket can be very large, so it’s unrealistic to keep track of all possible values we see, and we need a proxy to estimate the number of distinct values. The trick is to find some property that only very rare numbers satisfy, and the intuition is that if such property was satisfied by any of the numbers we saw, then there’s a good chance we saw a lot of numbers (or we’re very lucky).

To start, say that the property is “multiple of 2”. That is not a rare property, so if we assume the numbers in the balls do not have any special pattern, on average, we just need to find 2 different numbers to see that property being satisfied. Now say the property is “multiple of 4”. Less numbers satisfy this property, so the average number of distinct values we need to look at would be 4. We can keep going on for higher and higher powers of 2, and the average number of distinct values we look at would need to be as big.

This means that if we keep track of the largest power of 2 that divides any number we drew, we could estimate that that was the number of distinct values! The power of 2 property is interesting because it doesn’t change if there are duplicated values (which is great because we only want distinct values) and it doesn’t rely on the magnitude of the values as long as the set of values are uniformly distributed.

This estimate is not great though because the error is proportional to the number of distinct values found. For example, if we found out that the largest power of 2 divider was 2048 (2^11), the actual number of distinct values could be between 2048 and 4095 (2^12 – 1). Also, we could get unlucky and have all the balls in basket be have a single number, say 1024, so the estimate would be off by a lot.

Averaging

To reduce the chances of drawing a ball with a value that turns out to be a large power of 2 and throwing off the estimates, we divide the balls into m groups. The important property of the groups is that balls with the same values should go to the same group and the distribution should be uniform.

We can then find the largest power of 2 divider for each group and average them out. This would help in the case where we saw a single number 1024, because while one group would estimate it to be 1024, all the other groups would be find that the largest divisor was 1 (empty group).

To assign each number to a group, we could use modular arithmetic. For example, if we had 12 groups, we use the remainder of a number by 12 as an index to the group. The problem is that the elements assigned to a given group have a property that will bias the result. For example, for the group corresponding to the reminder 0, all numbers in there are obviously a multiple of 12 and hence 2^2. To make the numbers in the group unbiased, we need to discard the information used to assign them to groups. An easy way to achieve that is to represent each number as binary, use the first bits to decide which group to assign it to, and then discard those bits when computing the power of 2.

We can see that we can increase the number of groups to reduce errors but the tradeoff is that it requires to keep more information in memory. This is the idea first proposed by [2] and was called stochastic averaging.

Harmonic Mean

The HyperLogLog algorithm uses a different combination of the powers of two to obtain a better estimate than using averages: the harmonic mean [1]. To recall, the harmonic mean consists in averaging the reverse of the elements and then reversing the result. For example, the harmonic mean of 1, 4 and 4 is

The bulk of the HyperLogLog paper [1] is actually proving that this metric yields an estimate with smaller errors than a simple average.

Hashing

Not all real world input will be numbers that are distributed uniformly. For example, we might be interested in counting the number of distinct values of a string column in a database. We then need to transform these values using a hash function which maps these inputs to numbers uniformly.

### Algorithm

We are now ready to outline the core of the algorithm to estimate the cardinality of a multi-set of values.

Given a stream of n elements:

• For each element
• Hash it
• Use the first bits of the hash to determine to which group j to assign it to
• Keep track of the largest power of 2 that divides some value in the group, and store the exponent in M[j]

Finally, estimate the number of distinct values based on a harmonic mean of the values in M.

The Expected Value and the Alpha Correction

The expectation of the number of distinct values is given by [1]:

E is equivalent to the harmonic mean times $\alpha_m m$, and the constant $\alpha_m$ is a constant to correct some bias. For implementation purposes, the authors provide approximations to this constant:

Small and Big Range Corrections

When the estimated number of distinct elements is relatively small compared to the number of groups m, the author propose a correction.

Let $V$ be the number of groups to which no element was assigned to. If $V > 0$, then experiments show that for $E \le \frac{5}{2} m$, we can change the estimate to

Conversely, if the number of distinct elements is very high, closer to 2^32, then the probability of hash collision are very high. The correction accounts for that:

Implementation in Rust

We’ll now present an implementation using Rust. The first part of the algorithm consists in computing M, which we name more clearly as first_non_zero_by_bucket. The following code implements the pseudo-code described above:

 let m: u32 = 1 << b; let mut first_non_zero_by_bucket: Vec = vec![0 as u32; m as usize]; for element in elements { let hash_value = hash(&element); // Extracts the first b bits from hash_value to determine the bucket let bucket_index: usize = (hash_value & first_b_bits_mask) as usize; // Finds the position of the first 1 bit in the remaining bits let mut first_non_zero: u32 = first_non_zero_bit_position(hash_value >> b); first_non_zero_by_bucket[bucket_index] = cmp::max( first_non_zero_by_bucket[bucket_index], first_non_zero ); }

view raw
collecting.rs
hosted with ❤ by GitHub

We don’t need to know much Rust to read the code above. It’s worth mentioning that Rust is very strict about types, so we need to perform explicit conversions. Also we use 2 bit operations: one is to obtain the least significant k bits of an integer by using a bit mask. It relies on the fact that (2^k)-1 is a number with k bits 1 and doing a bitwise AND with any number has the effect of only extracting the first k bits of that number. The other trick is to divide a number by 2^k, which can be done by shifting the bits to the right, via the >> operator.

The hash function we use here is from the package farmhash, which is a Rust implementation of Google’s Farmhash, which in turn is a variant of Murmurhash [6]. It basically takes a string and shuffles its bits in a hopefully uniform way, generating a 32-bit integer:

 fn hash(value: &String) -> u32 { farmhash::hash32(&value.as_bytes()) }

view raw
hash.rs
hosted with ❤ by GitHub

first_non_zero_bit_position() is given by:

 fn first_non_zero_bit_position(input: u32) -> u32 { let mut remaining: u32 = input; let mut first_non_zero: u32 = 1; while (remaining & 1) == 0 && remaining > 1 { remaining /= 2; first_non_zero += 1; } first_non_zero }

The formula to obtain the expected number of distinct values is given by

The code below implements this function:

 let alpha: f64 = match b { 4 => 0.673, 5 => 0.697, 6 => 0.709, // b >= 7 _ => 0.7213/(1.0 + 1.079/(m as f64)) }; let mut indicator: f64 = 0.0; let base: f64 = 2.0; for first_non_zero in &first_non_zero_by_experiment { indicator += base.powf(–(*first_non_zero as f64)); } let m_multiplier = m as f64; let mut estimate: f64 = (m_multiplier * m_multiplier * alpha) / indicator;

view raw
estimate.rs
hosted with ❤ by GitHub

The values of alpha were discussed in the The Expected Value and the Alpha Correction above.

The correction for small and large ranges can be implemented as:

 const LARGE_ESTIMATE_THRESHOLD: f64 = 143165576.53333333; const TWO_TO_32: f64 = 4294967296.0; // 2^32 if estimate <= 2.5 * m_multiplier { // Small range correction let mut buckets_with_zero = 0; for first_non_zero in first_non_zero_by_experiment { if first_non_zero == 0 { buckets_with_zero += 1; } } if buckets_with_zero > 0 { estimate = m_multiplier * (m_multiplier / (buckets_with_zero as f64)).ln(); } } else if estimate > LARGE_ESTIMATE_THRESHOLD { // Large range correction estimate = –TWO_TO_32 * (1.0 – estimate/TWO_TO_32).ln(); }

view raw
correction.rs
hosted with ❤ by GitHub

The complete code is available on Github.

### Experiments

For values of b ranging from 4 to 16, I ran the program 100 times for n=100k, with numbers randomly selected from 1 to 100k. Then I plotted the results using the box plot chart using a R script:

In the chart above, the x-axis represents the number of experiments we divided the input into, and the y-axis represents the relative error compared to the actual value. We can see that as we increase the number of experiments, the errors go down.

This chart was also helpful in finding a bug in my implementation: the initial plot had a sudden spike for values larger than 11, which was when the small range correction applied. After some debugging, I realized that algorithm should you the natural logarithm, not log2. It was great to spot the bug via data visualization!

### Conclusion

I’ve been wanting to study the HyperLogLog algorithm for a while, and also one of my resolutions for this year is to learn Rust. It was a productive exercise to implement it.

I’ll post some future impression on Rust from someone with experience in C++ in a future post.

### References

[1] HyperLogLog: the analysis of a near-optimal cardinality estimation algorithm
[2] Probabilistic Counting Algorithms for Database Applications
[3] Thu Trang Pham – Curiosity #2: How does Prestodb implement approx_distinct?
[4] Probabilistic Data Structures for Web Analytics and Data Mining
[5] Gödel’s Lost Letter and P=NP: Philippe Flajolet 1948–2011
[6] GitHub: seiflotfy/rust-farmhash

# Bloom Filters

Burton Howard Bloom is a MIT alumni, who, mwhile employed by the Computer Usage Company, developed and published a data structure that later became famous as “Bloom Filter” [1].

In this post we’ll talk about Bloom Filters. We’ll give a description of the data structure and its operations. We’ll then study the theory that explains the performance of this data structure. Next, we’ll describe an implementation in Python, run some experiments and finally discuss applications.

### Introduction

Bloom filter is a probabilist data structure that enables inserting elements in a set and test whether an element belongs this set, sacrificing accuracy for lower memory usage.

More precisely, it can say an element is in the set when in fact it’s not (false positive), but if it says it’s not in the set, then it’s true. Also, the original variant doesn’t allow removal of elements [2].

In its original form, it doesn’t allow removing elements.

### Algorithm

The bloom filter structure consists of a bit array of size $m$. We then define k independent hash functions, which can distribute a given value into any of the $m$ positions uniformly.

Inserting an element consists in applying the $k$ hash functions to the element and set the bit in the $k$ positions returned by these functions to 1.

Querying for an element consists in applying the same $k$ hash functions and testing whether all the bits in the k positions returned by these functions is set to 1.

If the test returns positive, there’s a chance the element is there, but it could be a false positive, since the bits could have been set by the insertion of other elements. On the other hand, if it returns negative, we’re sure the element is there, because we never unset a bit. Also, as the number of elements in the set grow, the probability of false positives increases.

Time complexity. Both insertion and querying for an element are constant-time operations, since they are proportional to $k$, assuming random access memory and $O(1)$ hash functions.

Probability of false positives. The probability $p$ that the algorithm will return true when the element is not there (false positive), can be described in terms of $m$, $k$, $n$, through the following equation:

$p = (1 -e^{(kn/m)})^k$

In the appendix A, we provide an approximated derivation of this result.

Optimal number of hash functions. For given values $n$ and $m$, we have a choice of how many hash functions to use, $k$. Intuitively, if we use too few hash functions, we may increase the chances of collision when querying for an element, whereas if we use too many, the array will be filled up earlier and thus the collisions will also increase.

It’s possible to prove (a proof is sketched in appendix B), that for fixed $m$, $n$, the value of $k$ that minimizes $p$ is given by:

$k = ln(2)m/n$

### Python Implementation

In this experiment, we use the bitarray python module. It’s a C implementation that represents a bit array efficiently. It has conveniently overloaded operators so most of the time it’s like working with arrays.

We can define a sufficiently large value of our bit array length, $m$, for example:

from bitarray import bitarray
m = 1 << 20
bit = bitarray(m)

# bitarray doesn't guarantee the bits are all set to 0
# upon initialization
bit.setall(False)


### The murmur algorithm

We need to pick a hash function to use with our Bloom filter. References such as [3] suggest using an algorithm called Murmur. This stack exchange thread has a nice survey about different hash functions.

Murmur is a hash algorithm developed by Austin Appleby. Even though it’s not suitable for cryptographic applications, in practice it is very fast and tends to distributed real world instances well.

According to the author, the name comes from an early implementation detail, which used multiply-rotate-multiply-rotate operations to mix the internal state of the hash (hence MuR-MuR).

From the source code, this algorithm seems to shuffle the bits from the input key by using multiplication, shifting and xor operations.

The author mentions using simulated annealing to find some of the constants of the final mix of algorithm. This final part is used to cause the avalanche effect, in which very similar inputs are hashed to very different values.

There’s is a Python library called pyhash that has an interface for several hash functions (implemented in C++).

To install it, we can use the easy_install command. Easy Install is a python module that is used to download packages. pyhash is particular is available in the default repository, PyPI (python package index), so all we need to do is:

> sudo easy_install pyhash


To use it within a Python script:

from pyhash import murmur3_x64_128

hasher = murmur3_x64_128()
hash_value = hasher(input)


### Families of hash functions

In Jonathan Ellis‘s post, he mentions a strategy for generating a family of hash functions. He cites a paper from Kirsch and Mitzenmacher [4] which shows we can generate a set of $k$ hash functions from 2 independent hash functions in the following way:

$h_k(x) = h_A(x) + i \cdot h_B(x), \, i = 0, \cdots, k-1$

In practice, since murmur3_x64_128() generates 128-bit hashes, we can use the first 64-bits as the first hash function, and the last 64-bits as the second.

from pyhash import murmur3_x64_128

hasher = murmur3_x64_128()
h128 = hasher(input)
h64l = h128 & ((1L << 64) - 1)
h64u = h128 >> 64

hashes = map(
lambda i: (h64l + i*h64u) % k,
range(k)
)


All the code is available on github.

### Experiments

Before running experiments with our Bloom filter implementation, let’s try to visualize the distribution of the Murmur hash functions. For this first experiment, we compute the hash function of keys ranging form 1 to 10k and module it 10k, so they’re all in the same range of values. Then we plot the keys vs. their hashes in a scatter plot, by exporting this data as a CSV file and rendering in R (using ggplot2):

ggplot(data, aes(x=i, y=h_i)) + geom_point(size=1.5)


From the chart, we can see it apparently distributes the keys uniformly, without any obvious pattern. This is not sufficient to prove it’s a good hash function, but it’s a necessary property.

To test the family of hashes, we can add a new column, which has i for the i-th hash function. Then, in ggplot, we render each class with a different color with the following command:

ggplot(data, aes(x=i, y=h_i, color=class))
+ geom_point(size=1.5) + facet_grid(class ~ .)


We can now compare the distributions between each of the hashes by plotting scatter plots for each pair i, j. Now we use another ggplot2 command to generate a 5×5 grid of charts:

ggplot(data, aes(x=i, y=h_i, color=class))
+ geom_point(size=1.5) + facet_wrap(~class, ncol=5)


We can see each chart has a uniform pattern, which is a good sign if we’re looking for independent hash functions. For comparison, let’s plot the same chart for a set of hash functions derived only from a single hash function, for example

$h_i(x) = i \cdot h_A(x), i = 1, \cdots, k$

Now, let’s analyze the false positive rate. We can do that by counting the number of false positives FP and the true negatives TN and evaluating:

$\dfrac{FP}{FP + TN}$

We fix $m$ and $k$ and shuffle an array from 0 to n-1 to simulate a sampling without replacement. After each inserted element of this array, we go over all non-inserted elements $(FP + TN)$ and count how many of those our Bloom filter thinks they are in the set $FP$. We then plot a line chart with the number of elements inserted vs. the false positive rate, using qplot() R function:

qplot(
data$n, data$prob,
ylab='Pr(false positive)',
xlab='N elements'
)


Let’s plot the effect of increasing the size of the bit array, using the optimal k and increasing the number of elements. In this particular experiment, we use n=250 and try different values of m: 100, 250, 500, 1000

qplot(
data$n, data$prob,
ylab='Pr(false positive)',
xlab='N elements',
size=I(1),
color=data\$m
) + scale_colour_manual(
values=c("red", "blue", "green", "black"),
name="m values"
)


One strange behaviour seems to happen towards the right end of the chart for some curves — the false probability ratio seems to drop. I’ve double checked the code and it looks sane. One possible explanation is that we are calculating the probability over the non-inserted elements and as it approaches the end, the sample size is smaller, so the noise is more significant. Other than that, the probability of false positives tend to increase as we add more and more elements.

Now, let’s analyze the effect of the number of hash functions used. We fix it in n=250, m=1000 and try out different values of k=1, 5, 10, 50.

We can see that using many hash functions can be bad because it fills up the array too fast (k=50). In general k=5 and k=10 perform best for most part of the spectrum.

### Applications

Bloom filters are suitable where false positives are tolerable. For example, it can be used as a first step to avoid lookup to a cache. Periodically a cache server could send a Bloom filter to a local machine, corresponding to items it has cached. If the Bloom filter says an element is on cache, it might not be true, but if it says the item is not there, the local client can avoid a request to the cache and talk directly to the underlying storage.

Bloom filters use little memory, so they are suitable for network transmission.

Broder and Mitzenmacher discuss a lot more examples [5], especially network applications.

### Conclusion

In this post we learned about Bloom filters. The idea itself is quite simple, by studying the theory we are able to review some probability and calculus concepts.

The implementation made us aware of details which are usually not discussed in theoretical texts, for example, which hash functions to use and how to generate a family of hash functions.

We used two Python libraries, pyhash and bitarray and learned a little bit about the Python packaging system, PyPI. We got some experience with the ggplot2 R library, which I plan to post about later.

What we didn’t cover was variants of Bloom filters, like the count version, which allows deletion. Chazelle et al., introduced the Bloomier Filters, which is a generalization of Bloom filters [6].

### References

[1] Quora – Where can one find a photo and biographical details for Burton Howard Bloom, inventor of the Bloom filter?
[2] Wikipedia – Bloom Filter
[3] Spyved – All you ever wanted to know about writing bloom filters
[4] Less hashing, same performance: Building a better Bloom filter – A Kirsch, M. Mitzenmacher.
[5] Network Applications of Bloom Filters – A. Broder, M. Mitzenmacher
[6] Bloomier Filters – B. Chazelle et al.
[7] Probability and Computing: Randomized Algorithms and Probabilistic Analysis – M. Mitzenmacher, E. Upfal

### Appendix A: Probability of false positives

Let’s consider an array of m bits and the insertion operation. The probability of a given bit to be set by one of the hash functions is $1/m$. Conversely, the probability of a bit not being set is $1 - 1/m$. By the hypothesis of independent hash functions, the probability of a bit not being set by any of the $k$ hash functions is

$\left(1 - \frac{1}{m}\right)^k$

Inserting $n$ elements is equivalent to running the $k$ functions $n$ times, so after inserting $n$ elements, the probability that a given bit is set to 0 is

$\left(1 - \frac{1}{m}\right)^{kn}$

or the probability it’s set to one as

$1 - \left(1 - \frac{1}{m}\right)^{kn}$

Now, the probability of a false positive $p$, is the probability $k$ positions being set for a element that was never inserted.

$p = \left(1 - \left(1 - \frac{1}{m}\right)^{kn}\right)^k$

According to [2] though, this derivation is not accurate, since assumes that the probability of each bit being set is independent from each other. The Wikipedia article sketches a more precise proof from Mitzenmacher and Upfal [7] which arrives at the same result.

To simplify the equation, we’ll use the following limit:

$\lim_{x\to\infty} \left( 1 - 1/x \right) = 1/e$

So, for large values of m, we can assume that

$(1 - 1/m)^{kn} = (1 - 1/m)^{(mkn)/m} \approx (1/e)^{(kn/m)}$

Finally we have:

$p = (1 -e^{(kn/m)})^k$

### Appendix B: Optimal number of hash functions

If we make $p$ a function of $k$, it’s possible to prove it has a global minimum for positive $k$‘s so we can use the derivative to find the value of $k$ that minimizes the false positive probability.

To find the derivative, we can use the generalized power rule to get

$\ln(1 - e^{-ck}) + \frac{k}{(1 - e^{-ck})} c e^{-ck}$

where $c=n/m$ is a constant. If we make $y = e^{-ck}$ (and use $k = \ln(y)/-c$)

$\ln(1 - y) - \frac{\ln(y)}{(1 - y)} y$

By making it equals to 0, we get the equation:

$\ln(1 - y) (1-y) = \ln(y) y$

To solve this for $y$, let $x = 1 - y$,

$\ln(x) x = \ln(y) y$

Since $\ln(n)$ and $n$ are both monotonic functions, so is $\ln(n)n$. Then if () holds, we can conclude that $x = y$, because otherwise, if $x > y$, $\ln(x) x > \ln(y) y$ and if $x < y$, $\ln(x) x < \ln(y) y$.

Replacing $x$ back, we find out that $y = 1/2$ and finally $k = ln(2)m/n$

# Skip Lists in Python

Skip list is a probabilistic data structure that allows efficient search, insertion and removal operations. It was invented by William Pugh [1] in 1989.

Other structures that have efficient operations are self-balancing binary trees, such as AVL, Red-black and splay tree. But they are often considered difficult to implement.

On the other hand, skip lists are much like multiple linked lists with some randomization.

Figure 1: Skip list

In the first level, we have a regular linked list with the elements sorted. Each element of this list has a probability $p$ to be also present in the level above. The second level will probably contain fewer elements and each of these elements will also have a chance $p$ to be on the third level, and so on. Figure 1 shows an example of a skip list.

We’ll implement a simple version of the skip list in python. To start, we define a skip list node. We’ll represent each level where the node appears by a list of pointers to the next nodes.

class SkipNode:
def __init__(self, height = 0, elem = None):
self.elem = elem
self.next = [None]*height


Our skip list is just a sentinel skip node with height initially set to 0 and that stores a null element.

class SkipList:
def __init__(self):


Now, let’s implement the search operation of this list.

### Search

To search for an element $q$ in a skip list we begin in the topmost level of the header. We go through the list in this level until we find node with the largest element that is smaller than $q$.

We then go to the level below and search again for node with the largest element that is smaller than $q$, but this time we began the search from the node we found in the level above.

When we find such node, we go down again and repeat this process until we reach the bottom level. The node $x$ found in the bottom level will be the largest element that is smaller than $q$ in the whole list and if $q$ is in this list, it will be to the right of $x$.

Also, we want to keep the nodes found in each level right before going down to the level below since it will make the insertion and deletion operations very simple.

This idea can be translated into the following code:

def updateList(self, elem):

while x.next[i] != None and \
x.next[i].elem < elem:
x = x.next[i]
update[i] = x

return update


It returns a list of nodes in each level that contains the greatest value that is smaller than elem.

The actual find function returns the node corresponding to the query element or None if it is not present in the skip list.

def find(self, elem, update = None):
if update == None:
update = self.updateList(elem)
if len(update) > 0:
candidate = update[0].next[0]
if candidate != None and candidate.elem == elem:
return candidate
return None


#### Complexity

The complexity of the search is given by the following Theorem [2]:

Theorem: The number of moves in a search is $O(\log n)$ with high probability.

By high probability we mean that we can set an arbitrarily high probability by increasing the constant hidden in the $O()$ notation. The proof of this Theorem is sketched in Appendix A.

### Insertion

The insertion consists in deciding the height of the new node, using randomHeight() and for each of the levels up to this height, insert this new node after the node specified in update.

def insert(self, elem):

node = SkipNode(self.randomHeight(), elem)

update = self.updateList(elem)
if self.find(elem, update) == None:
for i in range(len(node.next)):
node.next[i] = update[i].next[i]
update[i].next[i] = node


### Deletion

The deletion is pretty much like the insertion, but now we delete the node found using find() from all levels in which it appears.

def remove(self, elem):

update = self.updateList(elem)
x = self.find(elem, update)
if x != None:
for i in range(len(x.next)):
update[i].next[i] = x.next[i]


Note that for the sake of simplicity we do not resize head.next when the lists at top levels become empty. It does not change the theoretical complexity in the worst case, but in practice it may improve performance.

### Computational Experiments

The complete implementation of skip list is available at Github as well as some test cases and a simple implementation of a linked list.

Our computational experiments consist in comparing the execution time for linked lists ($O(n)$ per operation), our simple skip list ($O(\log n)$ with high probability) and an implementation of a red-black tree (worst-case $O(\log n)$ per operation).

We ran 10000 insertions in each of these structures with a random sequence, an increasing sequence and a decreasing sequence. We measure the CPU time in seconds:

Table 1: Times for insertion in different structures

As we can see, for random input Red-black tree and Skip list have similar performance. For increasing and decreasing sequences, Skip list performed better than Red-black tree because the former is unaffected by the ordering of insertions, while the latter has to make many balancing operations in such cases.

As for linked lists, we can verify it’s much slower than the other two structures since it has $O(n)$ worst case insertion time, but it outperforms when the elements are inserted in decreasing order since in this case insertion is $O(1)$ :)

### Conclusion

There is a combination of insertions and removals that may degenerate a skip list to a linked list, so the operations of searching, inserting and deleting become $O(n)$ in the worst case scenario, though it is very unlikely to happen.

Demaine’s analysis [2] is stronger than Pugh’s [1]. The latter proves that the cost of the search is $O(\log n)$ in average (i.e. expected cost) while the former proves it is $O(\log n)$ for most of the cases.

### References

[1] Skip Lists: A Probabilistic Alternative to Balanced Trees – W. Pugh.
[2] Introduction to Algorithms MIT – Lecture 12: Skip Lists – Erik Demaine

### Appendix A: Proofs

In this section we present the proof of the Theorem stated in the post. Let $L(n) = \log_{1/p} n$.

Before proving the Theorem, let’s prove the following Lemma:

Lemma: The number of levels in a skip list is $O(L(n))$ with high probability.

Proof:

Let’s consider the error-probability, that is, the probability that there are more than $c \cdot L(n)$ levels. We use the Boole’s inequality which says that for a set of events $\{E_1, E_2, \cdots, E_k\}$:

$Pr\{E_1 \cup E_2 \cup \cdots \cup E_k \} \le Pr\{E_1\} + Pr\{E_2\} + \cdots + Pr\{E_k \}$

Thus,

$Pr\{ \mbox{max level } \ge c \cdot L(n) \} \le$
$n \cdot Pr\{\mbox{node level } \ge c \cdot L(n) \}$

Since each node height is given by a geometric distribution, we have that for some given level $x$:

$Pr\{\mbox{node level } = x\} = p^{x - 1}(1 - p)$

$Pr\{\mbox{node level } < x\} = \sum_{i = 0}^{x - 1} Pr\{\mbox{node level } = x\} =$
$\, (1 - p)\sum_{i = 0}^{x - 1} p^i = 1 - p^{k}$

$Pr\{\mbox{node level } \ge x \} = p^{k}$

Also, $p^{c \cdot L(n)} = p^{c \cdot \log_{1/p} n} = \frac{1}{n^{c}}$

Finally,

$Pr\{ \mbox{max level } < c \cdot L(n) \} \ge 1 - \frac{1}{n^{c - 1}}$

Thus, for a sufficiently large constant $c$, we have a very high probability, which proves the Lemma.

Theorem: The number of moves in a search is $O(\log n)$ with high probability.

Proof:

Let’s prove that the number of moves in a search is $O(L(n))$ with high probability.

First, consider the reversed path to find the returned element. This reversed path consists of up and left movements along the skip list. Note first that the number of up movements is bounded by the levels of the skip list.

An up movement is done with probability $p$, which is the case that the current element has at least one more level. Otherwise a left movement is done with probability $1 - p$.

Then, the length of the path is given by the number of movements until we reach $c \cdot L(n)$ up movements. We claim that such number of movements is $O(L(n))$ with high probability.

To prove our claim, let the number of movements be $c'cL(n)$ for some other constant $c'$. In [2], some combinatoric relations are used to show that

$Pr\{\mbox{\# up moves} \le cL(n) \mbox{ among } c'cL(n) \mbox{ moves}\} \le \frac{1}{n^\alpha}$

where $\alpha = ((c'- 1) - L(c'e)) \cdot c$. We also have the converse:

$Pr\{\mbox{\# up moves} > cL(n) \mbox{ among } O(cL(n)) \mbox{ moves}\} = 1 - \frac{1}{n^\alpha}$

Since $(c'- 1) > L(c'e)$ for sufficiently large values of $c'$, we can choose $c'$ to make $\alpha$ arbitrarily large, which makes the probability above very high, proving the claim.

We conclude that the number of movements is $O(cL(n)) = O(L(n))$ with high probability.