Locality Sensitive Hashing in R

Introduction

In the next series of posts I will try to explain base concepts Locality Sensitive Hashing technique.

Note, that I will try to follow general functional programming style. So I will use R’s Higher-Order Functions instead of traditional R’s *apply functions family (I suppose this post will be more readable for non R users). Also I will use brilliant pipe operator %>% from magrittr package. We will start with basic concepts, but end with very efficient implementation in R (it is about 100 times faster than python implementations I found).

The problem

Imagine the following interesting problem. We have two very large social netwotks (for example facebook and google+), which have hundreds of millions of profiles and we want to determine profiles owned by same person. One reasonable approach is to assume that these people have nearly same or at least highly overlapped sets of friends in both networks. One well known measure for determining degree of similarity of sets is Jaccard Index:

Set operations are computationally cheap and straightforward solution seems quite good. But let’s try to estimate computational time for duplicates detection for only people with name “John Smith”. Imagine that in average each person has 100 friends:

# for reproducible results
set.seed(seed = 17)
library('microbenchmark')
# we will use brilliant pipe operator %>%
library('magrittr')
jaccard <- function(x, y) {
  set_intersection <- length(intersect(x, y))
  set_union <- length(union(x, y))
  return(set_intersection / set_union)
}
# generate "lastnames"
lastnames <- Map(function(x) paste(sample(letters, 3), collapse = ''), 1:1e5) %>% unique
print(head(lastnames))

[[1]]
[1] "eyl"

[[2]]
[1] "ukm"

[[3]]
[1] "fes"

[[4]]
[1] "fka"

[[5]]
[1] "vuw"

[[6]]
[1] "ypg"

friends_set_1 <- sample(lastnames, 100, replace = F)
friends_set_2 <- sample(lastnames, 100, replace = F)
microbenchmark(jaccard(friends_set_1, friends_set_2))

Unit: microseconds
                                  expr    min      lq     mean median
 jaccard(friends_set_1, friends_set_2) 31.888 33.7735 37.29152 34.765
      uq     max neval
 36.1155 141.597   100

One operation takes 50 microseconds in average (on my machine). If we have 100000 of peoples with name John Smith and we have to compare all pairs, total computation will take more than 100 hours!

hours <- (50 * 1e-6) * 1e5 * 1e5 / 60 / 60
hours

[1] 138.8889

Of course this is unappropriate because of complexity of our brute-force algorithm.

Minhashing

To solve this kind problem we will use Locality-sensitive hashing – a method of performing probabilistic dimension reduction of high-dimensional data. It provides good tradeoff between accuracy and computational time and roughly speaking has complexity.
I will explain one scheme of LSH, called MinHash.
The intuition of the method is the following: we will try to hash the input items so that similar items are mapped to the same buckets with high probability (the number of buckets being much smaller than the universe of possible input items).
Let’s construct simple example:

set1 <- c('SMITH', 'JOHNSON', 'WILLIAMS', 'BROWN')
set2 <- c('SMITH', 'JOHNSON', 'BROWN')
set3 <- c('THOMAS', 'MARTINEZ', 'DAVIS')
set_list <- list(set1, set2, set3)

Now we have 3 sets to compare and identify profiles, related to same “John Smith”. From these sets we will construct matrix which encode relations between sets:

sets_dict <- unlist(set_list) %>% unique

m <- Map(f = function(set, dict) as.integer(dict %in% set), 
         set_list, 
         MoreArgs = list(dict = sets_dict)) %>% 
  do.call(what = cbind, .)

# This is equal to more traditional R's sapply call:
# m <- sapply(set_list, FUN = function(set, dict) as.integer(dict %in% set), dict = sets_dict, simplify = T)

dimnames(m) <- list(sets_dict, paste('set', 1:length(set_list), sep = '_'))
print(m)

         set_1 set_2 set_3
SMITH        1     1     0
JOHNSON      1     1     0
WILLIAMS     1     0     0
BROWN        1     1     0
THOMAS       0     0     1
MARTINEZ     0     0     1
DAVIS        0     0     1

Let’s call this matrix input-matrix. In our representation similarity of two sets from source array equal to the similarity of two corresponding columns with non-zero rows:

From table above we can conclude, that jaccard index between set_1 and set_2 is 0.75.
Let’s check:

print(jaccard(set1, set2))

[1] 0.75

column_jaccard <-  function(c1, c2) {
  non_zero <- which(c1 | c2)
  column_intersect <- sum(c1[non_zero] & c2[non_zero])
  column_union <- length(non_zero)
  return(column_intersect / column_union)
}
isTRUE(jaccard(set1, set2) == column_jaccard(m[, 1], m[, 2]))

[1] TRUE

All the magic starts here. Suppose random permutation of rows of the input-matrix m. And let’s define minhash function = # of first row in which column . If we will use independent permutations we will end with minhash functions. So we can construct signature-matrix from input-matrix using these minhash functions. Below we will do it not very efficiently with 2 nested for loops. But the logic should be very clear.

# for our toy example we will pick N = 4
N <- 4
sm <- matrix(data = NA_integer_, nrow = N, ncol = ncol(m))
perms <- matrix(data = NA_integer_, nrow = nrow(m), ncol = N)
# calculate indexes for non-zero entries for each column
non_zero_row_indexes <- apply(m, MARGIN = 2, FUN = function(x) which (x != 0) )
for (i in 1 : N) {
  # calculate permutations
  perm <- sample(nrow(m))
  perms[, i] <- perm
  # fill row of signature matrix
  for (j in 1:ncol(m))
    sm[i, j] <-  min(perm[non_zero_row_indexes[[j]]])
}
print(sm)

     [,1] [,2] [,3]
[1,]    3    3    1
[2,]    1    1    3
[3,]    1    1    2
[4,]    1    1    4

You can see how we obtain signature-matrix matrix after “minhash transformation”. Permutations and corresponding signatures marked with same colors:

You can notice that set_1 and set_2 signatures are very similar and signature of set_3 dissimilar with set_1 and set_2.

jaccard_signatures <-  function(c1, c2) {
  column_intersect <- sum(c1 == c2)
  column_union <- length(c1)
  return(column_intersect / column_union)
}
print(jaccard_signatures(sm[, 1], sm[, 2]))

[1] 1

print(jaccard_signatures(sm[, 1], sm[, 3]))

[1] 0

Intuition is very straighforward. Let’s look down the permuted columns and until we detect 1.

• If in both columns we find ones – (1, 1), then .
• In case (0, 1) or (1, 0) . So the probability over all permutations of rows that is the same as .

Moreover there exist theoretical guaranties for estimation of Jaccard similarity: for any constant there is a constant such that the expected error of the estimate is at most .

Implementation and bottlenecks

Suppose input-matrix is very big, say 1e9 rows. It is quite hard computationally to permute 1 billion rows. Plus you need to store these entries and access these values. It is common to use following scheme instead:

• Pick independent hash functions instead of premutations, .
• For each column and each hash function , keep a “slot” .
• will become the smallest value of for which column has 1 in row . I.e., gives order of rows for permutation.

So we end up with following ALGORITHM(1) from excellent Mining of Massive Datasets book:

for each row r do begin
  for each hash function hi do
    compute hi (r);
  for each column c
    if c has 1 in row r
      for each hash function hi do
        if hi(r) is smaller than M(i, c) then
          M(i, c) := hi(r);
end;

I highly recommend to watch video about minhashing from Stanford Mining Massive Datasets course.

Summary

Let’s summarize what we have learned from first part of tutorial:

• We can construct input-matrix from given list of sets. But actually we didn’t exploit the fact, that input-matrix is very sparse and construct it as R’s regular dense matrix. It is very computationally and RAM inefficient.
• We can construct dense signature-matrix from input-matrix. But we only implemented algorithm that is based on permutations and also not very efficient.
• We understand theorethical guaranties of our algorithm. They are proportional to number of independent hash functions we will pick. But how will we actually construct this family of functions? How can we efficiently increase number of functions in our family when needed?
• Our signature-matrix has small fixed number of rows. Each column represents input set and ~ . But we still have complexity, because we need to compair each pair to find duplicate candidates.

In the next posts I will describe how to efficently construct and store input-matrix in sparse format. Then we will discuss how to construct family of hash functions. After that we will implement fast vectorized version of ALGORITHM(1). And finally we will see how to use Locality Sensitive Hashing to determine candidate pairs for similar sets in time. Stay tuned!