Bloom Filters

Bloom filters are one of the really ingenious and simple building blocks for randomized data structures. A great summary is the paper by Broder and Mitzenmacher. In this post I will briefly review its key ideas since it forms the basis of the Count-Min sketch of Cormode and Muthukrishnan, it will also be necessary for an accelerated version of the graph kernel of Shervashidze and Borgwardt, and finally, a similar structure will be needed to compute data streams over time for a real-time sketching service.

At its heart a bloom filter uses a bit vector of length N and a set of k hash functions mapping arbitrary keys x into their hash values $$h_i(x) \in [1 .. N]$$ where $$i \in \{1 .. k\}$$ denotes the hash function. The Bloom filter allows us to perform approximate set membership tests where we have no false negatives but we may have a small number of false positives.

Initialize(b): Set all $$b[i] = 0$$

Insert(b,x): For all $$i \in \{1 .. k\}$$ set $$b[h_i(x)] = 1$$

Query(b, x): Return true if $$b[h_i(x)] = 1$$ for all $$i \in \{1 .. k\}$$, false otherwise

Furthermore, unions and intersections between sets are easily achieved by performing bit-wise OR and AND operations on the bloom hashes of the corresponding sets respectively.

It is clear that if we inserted x into the Bloom filter the query will return true, since all relevant bits in b are 1. To analyze the probability of a false positive take the probability of a bit being 1. After inserting m items using k hash functions on a range of N we have

$$\Pr(b[i] = \mathrm{TRUE}) = 1 - (1 - \frac{1}{N})^{k m} \approx 1 - e^{-\frac{km}{N}}$$

For a false positive to occur we need to have all k bits associated with the hash functions to be 1. Ignoring the fact that the hash functions might collide the probability of false positives is given by

$$p \approx (1 - e^{-\frac{km}{N}})^k$$

Taking derivatives with respect to $$\frac{km}{N}$$ shows that the minimum is obtained for $$\log 2$$, that is $$k = \frac{N}{m} \log 2$$.

One of the really nice properties of the Bloom filter is that all memory is used to store the information about the set rather than an index structure storing the keys of the items. The downside is that it is impossible to read out b without knowing the queries. Also note that it is impossible to remove items from the Bloom filter once they’ve been inserted. After all, we do not know whether some of the bits might have collided with another key, hence setting the corresponding bits to 0 would cause false negatives.

Real simple covariate shift correction

Imagine you want to design some algorithm to detect cancer. You get data of healthy and sick people; you train your algorithm; it works fine giving you high accuracy and you conclude that you’re ready for a successful career in medical diagnostics.

Not so fast …

Many things could go wrong. In particular, the distributions that you work with for training and those in the wild might differ considerably. This happened to an unfortunate startup I had the opportunity to consult for many years ago. They were developing a blood test for a disease that affects mainly older men and they’d managed to obtain a fair amount of blood samples from patients. It is considerably more difficult, though, to obtain blood samples from healthy men (mainly for ethical reasons). To compensate for that, they asked a large number of students on campus to donate blood and they performed their test. Then they asked me whether I could help them build a classifier to detect the disease. I told them that it would be very easy to distinguish between both datasets with probably near perfect accuracy. After all, the test subjects differed in age, hormone level, physical activity, diet, alcohol consumption, and many more factors unrelated to the disease. This was unlikely to be the case with real patients: Their sampling procedure had caused an extreme case of covariate shift that couldn’t be corrected by conventional means. In other words, training and test data were so different that nothing useful could be done and they had wasted significant amounts of money.

In general the situation is not quite so dire. Assume that we want to estimate some dependency $$p(y|x)$$ for which we have labeled data $$(x_i, y_i)$$. Alas, the observations $$x_i$$ are drawn from some distribution $$q(x)$$ rather than the ‘proper’ distribution $$p(x)$$. If we adopt a risk minimization approach, that is, if we want to solve

$$\mathrm{minimize}_{f} \frac{1}{m} \sum_{i=1}^m l(x_i, y_i, f(x_i)) + \frac{\lambda}{2} \|f\|^2$$

we will need to re-weight each instance by the ratio of probabilities that it would have been drawn from the correct distribution, that is, we need to reweight things by $$\frac{p(x_i)}{q(x_i)}$$. This is the ratio of how frequently the instances would have occurred in the correct set vs. how frequently it occurred with the sampling distribution $$q$$. It is sometimes also referred to as the Radon-Nikodym derivative. Such a method is called importance sampling and the following derivation shows why it is valid:

$$\int f(x) dp(x) = \int f(x) \frac{dp(x)}{dq(x)} dq(x)$$

Alas, we do not know $$\frac{dp(x)}{dq(x)}$$, so before we can do anything useful we need to estimate the ratio. Many methods are available, e.g. some rather fancy operator theoretic ones which try to recalibrate the expectation operator directly using a minimum-norm or a maximum entropy principle. However, there exists a much more pedestrian, yet quite effective approach that will give almost as good results: logistic regression.

After all, we know how to estimate probability ratios. This is achieved by learning a classifier to distinguish between data drawn from $$p$$ and data drawn from $$q$$. If it is impossible to distinguish between the two distributions then it means that the associated instances are equaly likely to come from either oneof the two distributions. On the other hand, any instances that can be well discriminated should be significantly over/underweighted accordingly. For simplicity’s sake assume that we have an equal number of instances from both distributions, denoted by $$x_i \sim p(x)$$ and $$x_i’ \sim q(x)$$ respectively. Now denote by $$z_i$$ labels which are 1 for data drawn from $$p$$ and -1 for data drawn from $$q$$. Then the probability in a mixed dataset is given by

$$p(z=1|x) = \frac{p(x)}{p(x) + q(x)}$$

Hence, if we use a logistic regression approach which yields $$p(z=1|x) = \frac{1}{1 + e^{-f(x)}}$$, it follows (after some simple algebra) that

$$\frac{p(z=1|x)}{p(z=-1|x)} = e^{f(x)}.$$

Now we only need to solve the logistic regression problem

$$\mathrm{minimize}_f \frac{1}{2m} \sum_{(x,z)} \log [1 + \exp(-z f(x))] + \frac{\lambda}{2} \|f\|^2$$

to obtain $$f$$. Subsequently we can use $$e^{f(x_i)}$$ as covariate shift correction weights in training our actual classifier. The good news is that we can use an off-the-shelf tool such as logistic regression to deal with a decidedly nonstandard estimation problem.

Graphical Models for the Internet

Here are a few tutorial slides I prepared with Amr Ahmed for WWW 2011 in Hyderabad next week. They describe in fairly basic (and in the end rather advanced) terms how one might use graphical models for the amounts of data available on the internet. Comments and feedback are much appreciated.

Memory Latency, Hashing, Optimal Golomb Rulers and Feistel Networks

In many problems involving hashing we want to look up a range of elements from a vector, e.g. of the form $$v[h(i,j)]$$ for arbitrary $$i$$ and for a range of $$j \in \{1, \ldots, n\}$$ where $$h(i,j)$$ is a hash function. This happens e.g. for multiclass classification, collaborative filtering, and multitask learning.

While this works just fine in terms of estimation performance, traversing all values of j leads to an algorithm which is horrible in terms of memory access patterns. Modern RAM chips are much faster (over 10x) when it comes to reading values in sequence than when carrying out random reads. Furthermore, random access destroys the benefit of a cache. This leads to algorithms which are efficient in terms of their memory footprint, however, they can be relatively slow. One way to address this is to bound the range of $$h(i,j)$$ for different values of j. Here are some ways of how we could do this:

1. Decompose $$h(i,j) = h(i) + j$$. This is computationally very cheap, it has good sequential access properties but it leads to horrible collisions should there ever be two $$i$$ and $$i’$$ for which $$|h(i) - h(i’)| \leq n$$.
2. Decompose $$h(i,j) = h(i) + h’(j)$$ where $$h’(j)$$ has a small range of values.
This is a really bad idea since now we have a nontrivial probability of collision as soon as the range of $$h’(j)$$ is less than $$n^2$$ due to the birthday paradox. Moreover, for adjacent values $$h(i)$$ and $$h(i’)$$ we will get many collisions.
3. Decompose $$h(i,j) = h(i) + g(j)$$ where $$g(j)$$ is an Optimal Golomb Ruler.
The latter is an increasing sequence of integers for which any pairwise distance occurs exactly once. In other words, the condition $$h(a) - h(b) = h(c) - h(d)$$ implies that $$a = c$$ and $$b = d$$. John Langford proposed this to address the problem. In fact, it solves our problem since there are a) no collisions for a fixed $$i$$ and b) for neighboring values $$h(i)$$ and $$h(i’)$$ we will get at most one collision (due to the Golomb ruler property). Alas, this only works up to $$n=26$$ since finding an Optimal Golomb Ruler is hard (it is currently unknown whether it is actually NP hard).
4. An alternative that works for larger n and that is sufficiently simple to compute is to use cryptography. After all, all we want is that the hash function $$h’(j)$$ has small range and that it doesn’t have any self collisions or any systematic collisions. We can achieve this by encrypting j using the key i to generate an encrypted message of N possible values. In other words we use
$$h(i,j) = h(i) + \mathrm{crypt}(j|i,N)$$
Since it is an encryption of j, the mapping is invertible and we won’t have collisions for a given value of j. Furthermore, for different i the encodings will be uncorrelated (after all, i is the key). Finally, we can control the range $$N>n$$ simply by choosing the encryption algorithm. In this case the random memory access is of bounded range, hence the CPU cache will not suffer from too many misses.

A particularly nice algorithm is the Feistel cipher. It works as follows: define the iterative map

$$f(x,y) = (y, x \mathrm{XOR} h(y))$$

Here $$h$$ is a hash function. After 4 iterations $$(x,y) \to f(x,y)$$ we obtain an encryption of $$(x,y)$$. Now use $$x=i$$ and $$y = j$$ to obtain the desired result. Basically we are trading off memory latency with computation (which is local).

Collaborative Filtering considered harmful

Much excellent work has been published on collaborative filtering, in particular in terms of recovering missing entries in a matrix. The Netflix contest has contributed a significant amount to the progress in the field.

Alas, reality is not quite as simple as that. Very rarely will we ever be able to query a user about arbitrary movies, books, or other objects. Instead, user ratings are typically expressed as preferences rather than absolute statements: a preference for Die Hard given a generic set of movies only tells us that the user appreciates action movies; however, a preference for Die Hard over Terminator or Rocky suggests that the user might favor Bruce Willis over other action heroes. In other words, the context of user choice is vital when estimating user preferences.

Hence if we attempt to estimate scores $$s_{ui}$$ of user $$u$$ regarding item $$i$$ it is important to use the context within which the ratings have been obtained. For instance, if we are given a session of items $$(i_1, \ldots, i_n)$$ out of which item $$i^*$$ was selected we might want to consider a logistic model of the form:

$$-\log p(i^*|i_1, \ldots, i_n) = \log \left[\sum_{i=1}^n e^{s_{ui}} \right] - s_{ui^*}$$

The option of no action is easy to add, simply by adding the null score $$s_{u0}$$ which captures the event of no action by a user.
Shuang Hong tried out this idea to get a significant performance improvement on a number of collaborative filtering datasets. Bottom line - make sure that the problem you’re solving is actually the one that a) generated the data and b) that will help you in practice. That is, in many cases matrix completion is not the problem you want to solve, even though it might win you benchmarks.

Why

Some readers might wonder why I’m writing this blog. Here’s an (incomplete) list:

• It’s fun.
• There are lots of fantastic blogs discussing the philosophy and big questions of machine learning (e.g. John Langford’s hunch.net) but I couldn’t find many covering simple tricks of the trade.
• Scientific papers sometimes obscure simple ideas. In the most extreme case, a paper will get rejected if the idea is presented in too simple terms (it happened to me more than once and the paper was praised once the simple parts had been obfuscated). Also, they need to come with ample evidence for why an idea works, strong theoretical guarantees and lots of experiments. This is all needed as a safeguard and it’s really really important. But it often hides the basic idea.
• Some ideas are really cute and useful but not big enough to write a paper about. It’s pointless to write 10 pages if the idea can be fully covered in 1 page. We’d need a journal of 1 page ideas to deal with this.
• Many practitioners are scared to pick up a paper with many equations but they might be willing to spend 10 minutes reading a blog post.
Hashing for Collaborative Filtering

This is a follow-up on the hashing for linear functions post. It’s based on the HashCoFi paper that Markus Weimer, Alexandros Karatzoglou and I wrote for AISTATS’10. It deals with the issue of running out of memory when you want to use collaborative filtering for very large problems. Here’s the setting:

Assume you want to do Netflix-style collaborative filtering, i.e. you want to estimate entries in a ratings matrix of (user, movie) pairs. A rather effective approach is to use matrix factorization, that is, to approximate $$M = U^\top V$$ where M is the ratings matrix, U is the (tall and skinny) matrix of features for each user, stacked up, and V is the counterpart for movies. This works well for the Netflix prize since the number of users and movies is comparatively small.

In reality we might have, say 100 million users for which we might want to recommend products. One option is to distribute all these users over several servers (similar to what a distributed hash table mapping does, e.g. for libmemcached). Alternatively, if we want to keep it all on one server, we’re facing the problem of having to store $$10^8 \cdot 100 \cdot 4 = 4 \cdot 10^10$$ bytes, i.e. 40 GB if we assume to allocate 400 Bytes per user (that’s a rather small footprint). That is 100 dimensions per user. Usually this is too big for all but the biggest servers. Even worse, suppose that we have user-churn. That is, new users might be arriving while old users disappear (obviously we don’t know whether they’ll ever come back again so we don’t really want to de-allocate the memory devoted to them). Obviously we cannot add more RAM. One possible solution is to store the data on disk and request it whenever a user arrives. This will cost us 5-10ms latency. An SSD will improve this but it still limits throughput. Moreover, it’ll require cache management algorithms to interact with the collaborative filtering code.

Here’s a simple alternative: apply the hashing trick that we used for vectors to matrices. Recall that in the exact case we compute matrix entries via

$$M[i,j] = \sum_{k=1}^{K} U[i,k] V[j,k]$$

Now denote by $$h_u$$ and $$h_v$$ hash functions mapping pairs of integers to a given hash range $$[1 \ldots N]$$. Moreover, let $$\sigma_u$$ and $$\sigma_v$$ be corresponding Rademacher hash functions which return a binary hash in $$\{\pm 1\}$$. Now replace the above sum via

$$M[i,j] = \sum_{k=1}^{K} u[h_u(i,k)] \sigma_u(i,k) v[h_v(j,k)] \sigma_v(j,k)$$

What happened is that now all access into U is replaced by access into a vector u of length N (and the same holds true for V). Why does this work: firstly, we can prove that if we construct u and v from U and V via

$$u[k] = \sum_{h_u(i,j) = k} \sigma_u(i,j) U[i,j] \text{ and } v[k] = \sum_{h_v(i,j) = k} \sigma_v(i,j) V[i,j]$$

then the approximate version of $$M[i,j]$$ converges to the correct $$M[i,j]$$ with variance $$O(1/N)$$ and moreover that the estimate is unbiased. Getting the exact expressions is a bit tedious and they’re described in the paper. In practice, things are even better than this rate: since we never use U and V but always u and v we simply optimize with respect to the compressed representation.

One of the advantages of the compressed representation is that we never really need to have any knowledge of all the rows of U. In particular, rather than mapping user IDs to rows in U we simply use the user ID as the hash key. If a new user appears, memory is effectively allocated to the new user by means of the hash function. If a user disappears, his parameters will simply get overwritten if we perform stochastic gradient descent with respect to the u and v vectors. The same obviously holds for movies or any other entity one would like to recommend.

Bottom line - we now can have fast (in memory) access to user parameters regardless of the number of users. The downside is that the latency is still quite high: remember that the hash function requires us to access $$u[h_u(i,k)]$$ for many different values of k. This means that each access in k is a cache miss, i.e. it’ll cost us 100-200ns RAM latency rather than the 10-20ns we’d pay for burst reads. How to break this latency barrier is the topic of one of the next posts.

Priority Sampling

Tamas Sarlos pointed out a much smarter strategy on how to obtain a sparse representation of a (possibly dense) vector: Priority Sampling by Duffield, Lund and Thorup (Journal of the ACM 2006).  The idea is quite ingenious and (surprisingly so) essentially optimal, as Mario Szegedy showed. Here’s the algorithm:

For each $$x_i$$ compute a priority $$p_i = \frac{x_i}{a_i}$$ where $$a_i \sim U(0, 1]$$ is drawn from a uniform distribution. Denote by $$\tau$$ the k+1 largest such priority. Then pick all k indices i which satisfy $$p_i > \tau$$ and assign them the value $$s_i = \mathrm{max}(x_i, \tau)$$. All other coordinates are set to $$s_i = 0$$.

This provides an estimator with the following properties:

1. The variance is no larger than that of the best k+1 sparse estimator.
2. The entries $$s_i$$ satisfy $$\mathbf{E}[s_i] = x_i$$
3. The covariance vanishes, i.e. $$\mathbf{E}[s_i s_j] = x_i x_j$$

Note that we assumed that all $$x_i \geq 0$$. Otherwise simply apply the same algorithm to $$|x_i|$$ and then return signed versions of the estimate.

Random elements from a stream

This is a classic trick when dealing with data streams. It shows how to draw a random element from a sequence of instances without knowing beforehand how long the sequence is and which symbols occur.

Let us first assume that we knew the identities of all symbols. In this case finding a random symbol would be easy. All we require is that for each symbol s we draw a random variable $$\xi_s \sim U[0,1]$$ from some distribution and subsequently we choose the symbol $$s^* = \mathrm{argmin}_s \xi_s.$$ Since each s has equal probability of being associated with the smallest value $$\xi_s$$ it follows that the draw is uniformly random.

Now assume that instead of requesting a random variable $$\xi_s$$ we simply compute the hash of s via $$h(s)$$ and we set $$s^* = \mathrm{argmin}_s h(s).$$ For a draw from the space of hash functions this again is uniform. The advantage is that we essentially determined all the random bits when selecting $$h$$ rather than at the time when we want to compute its value $$h(s)$$. The second advantage is that we can now simply keep track of what is the currently smallest value of $$h(s)$$ and update as we go along. We have the following algorithm:

INIT
hstar = MAXINT
n = 0
sstar = NONE
FOR ALL incoming s DO
IF h(s) = hstar:
n = n + 1
ELSE IF h(s) < hstar:
n = 1
hstar = h(s)
sstar = s
RETURN (sstar, n)


This algorithm will provide item counts for a random element of the sequence. If you want more than one sample, simply keep a list of the symbols with the k smallest hash values and their associated counts. Such algorithms can be used to compute the variance or other moments of a sequence.

Sparsifying a vector/matrix

Sometimes we want to compress vectors to reduce memory footprint or to minimize computational cost. This occurs, e.g. when we want to replace a probability vector by a sample from it. Without loss of generality, assume that our vector v has only nonnegative entries and that its entries sum up to 1, i.e. that it is a probability vector. If not, we simply apply our sampling scheme to the following:

There are a few strategies one could use to obtain a sparse vector s:

Sampling uniformly at random

For an n-dimensional vector simply use the sparsification rule: pick a random number j in [1..n]. Set the j-th coordinate in the sparse vector to  and all other terms to 0. For k draws from this distribution, simply average the draws. This scheme is unbiased but it has very high variance since there’s a very high chance we’ll miss out on the relevant components of v. This is particularly bad if v only has a small number of nonzero terms. DO NOT USE.

Sampling according to v at random

A much better method is to draw based on the size of the entries in v. That is, treat v as a probability distribution and draw from it. If we draw coordinate j in [1..n] then set the j-th coordinate of the sparse vector to  and all other terms to 0. As before, for k draws from this distribution, simply average. The advantage of this method is that now the entries of s are within the range [0, 1] and that moreover we hone in on the nonzero terms with high weight. But we can do better …

Sampling according to v with replacement

The key weakness in the above methods can be seen in the following example: for a vector of the form v = (0.99, 0.005, 0.005, 0, 0, …., 0) we would need to perform many samples with replacement from v until we even draw a single instance from the second or third coordinate. This is very time consuming. However, once we’ve drawn coordinate 1 we can inspect the corresponding value in v at no extra cost and there’s no need to redraw it.

Enter sampling with replacement: draw from v, remove the coordinate, renormalize the remainder to 1 and repeat. This gives us a sample drawn without replacement (see the code below which builds a heap and then adds/removes things from it while keeping stuff normalized). However, we need to weigh things differently based on how they’re drawn. Here’s what you do when drawing k terms without replacement:

Here we initialize γ=1 and j is the index of the item drawn at the i-th step. Below is some (I hope relatively bug-free) sample code which implements such a sampler. As far as I recall, I found partial fragments of the heap class on stackoverflow but can’t quite recall where I found it and who to attribute this to. This could be made more efficient in C++ but I think it conveys the general idea.

class arrayHeap:
# List of terms on array
def __init__(self, probabilities):
self.m = len(probabilities)            # sample size
self.b = int(math.ceil(math.log(self.m,2))) # bits
self.limit = 1 << self.b
self.heap = numpy.zeros(self.limit << 1)    # allocate twice the size

# allocate the leaves
self.heap[self.limit:(self.limit + self.m)] = probabilities
# iterate up the tree (this is O(m))
for i in range(self.b-1,-1,-1):
for j in range((1 << i), (1 << (i + 1))):
self.heap[j] = self.heap[j << 1] + self.heap[(j << 1) + 1]

# remove term from index (this costs O(log m) steps)
def delete(self, index):
i = index + self.limit
w = self.heap[i]
for j in range(self.b, -1, -1):
self.heap[i] -= w
i = i >> 1

# add value w to index (this costs O(log m) steps)
i = index + self.limit
for j in range(self.b, -1, -1):
self.heap[i] += w
i = i >> 1

# sample from arrayHeap
def sample(self):
xi = self.heap[1] * numpy.random.rand()
i = 1
while (i < self.limit):
i <<= 1
if (xi >= self.heap[i]):
xi -= self.heap[i]
i += 1
return (i - self.limit)

#Input: normalized probabilities, sample size
def sampleWithoutReplacement(p, n):
heap = arrayHeap(p)
samples = numpy.zeros(n, dtype=int) # result vector
for j in range(n):
samples[j] = heap.sample()
heap.delete(samples[j])
return samples

#Input: dense matrix p, resampling rate n
def sparsifyWithoutReplacement(p, n=1):
(nr, nc) = p.shape
res = scipy.sparse.lil_matrix(p.shape) # allocate sparse container
for i in range(nr):
tmp = sampleWithoutReplacement(p[i,:], n)