We introduce simple, efficient algorithms for computing a MinHash of a probability distribution, suitable for both sparse and dense data, with equivalent running times to the state of the art for both cases. The collision probability of these algorithms is a new measure of the similarity of positive vectors which we investigate in detail. We describe the sense in which this collision probability is optimal for any Locality Sensitive Hash based on sampling. We argue that this similarity measure is more useful for probability distributions than the similarity pursued by other algorithms for weighted MinHash, and is the natural generalization of the Jaccard index. MinHashing BID0 is a popular Locality Sensitive Hashing algorithm for clustering and retrieval on large datasets. Its extreme simplicity and efficiency, as well as its natural pairing with MapReduce and key-value datastores, have made it a basic building block in many domains, particularly document clustering BID0 BID1 and graph clustering BID2 BID3.Given a finite set U, and a uniformly random permutation π, the map X → arg min i∈X π(i) provides a representation of any subset X of U that is stable under small changes to X. If X, Y are both subsets of U the well-known Jaccard index BID4 Practically, this random permutation is generated by applying some hash function to each i with a fixed random seed, hence "MinHashing."In order to hash objects other than sets, Chum et al. BID5 introduced two algorithms for incorporating weights in the computation of MinHashes. The first algorithm associates constant global weights with the set members, suitable for idf weighting. The collision probability that is ∑ i∈X∩Y wi ∑ i∈X∪Y wi. The second algorithm computesMinHashes of vectors of positive integers, yielding a collision probability of J W (x, y) = ∑ i min (x i, y i) ∑ i max (x i, y i) Subsequent work BID6 BID8 has improved the efficiency of the second algorithm and extended it to arbitrary positive weights, while still achieving J W as the collision probability. J W is one of several generalizations of the Jaccard index to non-negative vectors. It is useful because it is monotonic with respect to the L 1 distance between x and y when they are both L 1 normalized, but it is unnatural in many ways for probability distributions. If we convert sets to binary vectors x, y, with x i, y i ∈ {0, 1}, then J W (x, y) = J(X, Y). But if we convert these vectors to probability distributions by normalizing them so that x i ∈ {0, As a consequence, switching a system from an unweighted MinHash to a MinHash based on J W will generally decrease the collision probabilities. Furthermore, J W is insensitive to important differences between probability distributions. It counts all differences on an element in the same linear scale regardless of the mass the distributions share on that element. For instance, J W ((a, b, c, 0), (a, b, 0, c)) = J W ((a + c, b), (a, b + c)). This makes it a poor choice when the ultimate goal is to measure similarity using an expression based in information-theory or likelihood where having differing support typically in the worst possible score. For a drop-in replacement for the Jaccard index that treats its input as a probability distribution, we'd like it to have the following properties. 2) Not lower than the Jaccard Index when applied to discrete uniform distributions. 3) Sensitive to changes in support, in a similar way to information-based measures. 4) Easily achievable as a collision probability. J W fails all but the last.1) It isn't scale invariant, J W (αx, y) ̸ = J W (x, y).2) If the vectors are normalized to make it scale invariant, the values drop below the corresponding Jaccard index (equation 1.) 3) It is insensitive to changes in support. 4) Good algorithms exist, but they are non-trivial. Existing work has thoroughly explored improvements to Chum et al.'s second algorithm, while leaving their first untouched. In this work we instead take their first algorithm as a starting point. We extend it to arbitrary positive vectors (rather than sets with constant global weights) and analyze the . In doing so, we find that the collision probability is a new generalization of the Jaccard Index to positive vectors, which we here call J P. DISPLAYFORM0 The names used here, J W and J P, are chosen to reflect how each function interprets x and y, and the conditions under which they match the original Jaccard index. J W treats a difference in magnitude the same as any other difference, so treats vectors as "weighted sets." J P is scale invariant, so any input is treated the same as a corresponding probability distribution. The primary contribution of this work is to derive and analyze J P, and to show that in many situations where the objects being hashed are probability distributions, J P is a more useful collision probability than J W.We will describe the sense in which J P is an optimal collision probability for any LSH based on sampling. We will prove that if the collision probability of a sampling algorithm exceeds J P on one pair, it must sacrifice collisions on a pair that has higher J P.We will motivate J P's utility by showing experimentally that it has a tighter relationship to the JensenShannon divergence than J W, and is more closely centered around the Jaccard index than J W. We will even show empirically that in some circumstances, it is better for retrieving documents that are similar under J W than J W itself (and consequently, sometimes better for retrieving based on L 1 -distance.) Let h: [n] → be a pseudo-random hash mapping every element 1 ≤ i ≤ n to an independent uniform random value in. Over a non-negative vector x, define DISPLAYFORM0 For brevity, we will be using the extended real number system in which 1/∞:= 0. Each term is an exponentially distributed random variable with rate x i, DISPLAYFORM1 so it follows that DISPLAYFORM2 This well known and beautiful property of exponential random variables derives from the fact that DISPLAYFORM3 Proof. Any monotonic transform on DISPLAYFORM4 will not change the arg min, so multiplying each x i by a positive α won't either. H(αx) = H(x). Thus for x i, y i > 0 DISPLAYFORM5 ). ). Consequently, Repeating this process for all i in the intersection yields DISPLAYFORM0 DISPLAYFORM1 Continuing the notational convention, we will refer to hashing algorithms that achieve J P as their pair collision probability as P-MinHashes, and algorithms that achieve J W as W-MinHashes. While J P's expression is superficially awkward, we can aid intuition by representing it in other ways. The simplest interpretation is to view it as a variant of J W. We rewrite J W allowing one input to be rescaled before computing each term, DISPLAYFORM0 j ) and choose the vector α to maximize this generalized J W. If α i x i > y i, increasing α i raises only the denominator. If α i x i < y i, increasing α i raises the numerator more than the denominator. So the optimal α sets α i x i = y i, and in J P.We can derive a more powerful representation by viewing the P-MinHash algorithm itself geometrically. A vector of k + 1 exponential random variables, when normalized to sum to 1, is a uniformly distributed random point on the unit k-simplex. Every point in the unit k-simplex is also a probability distribution over k+1 elements. Using these two facts we can construct the PMinHash as a function of the simplex as illustrated in FIG2.For a probability distribution x, mark the point on the unit simplex corresponding to x, (x 1, . . ., x k+1), and connect it to each of the corners of the simplex. These edges divide it into k + 1 smaller simplices that fill the unit simplex. Each internal simplex has volume proportional to the coordinate of x opposite to its unique exterior face. As a , a uniformly chosen point on the unit simplex will land in one of the sub-simplices with probability given by x. P-MinHashing is equivalent to sampling in this fashion, but holding the chosen point constant when sampling from each new distribution. The match probability is then proportional to the sum of the intersections of simplices that share an external face. This representation makes several properties obvious on small examples that we prove generally in the next section. When MinHashing is used with a key-value store, high collision probabilities are generally more efficient than low collision probabilities, because as we discuss in section VI, it is much cheaper to lower them than to raise them. For this reason, we are interested in the question of the highest collision probability that can be achieved through sampling. The constraint that the samples follow each distribution forces the collision probability to remain discriminative, but given that constraint, we would like to make it as high as possible to maximize flexibility and efficiency. Suppose for two distributions, x and y, we want to choose a joint distribution that maximizes Pr[H(x) = H(y)]. If we were concerned with only these two particular distributions in isolation, the upper bound of Pr[H(x) = H(y)] is given by the Total Variation distance, or equivalently 1 − L 1 (x, y)/2. Meeting this bound requires the probability mass where x exceeds y to be perfectly coupled with the mass where y exceeds x. Both the mass they share and the mass they do not must be perfectly aligned. Rather than just two, we want to create a joint distribution (or coupling) of all possible distributions on a given space where the collision probability for any pair is as high as possible. It is always possible to increase the collision probability of one pair at the expense of another so long as the chosen pair has not hit the Total Variation limit, so the kind of optimality we are aiming for is Pareto optimality. This requires that no collision probability be able to exceed ours everywhere; any gain on one pair must have a corresponding loss on another. This by itself would not be a very consequential bound for its retrieval performance. We really only desire high collision probabilities for items that are similar, and we would happily lower the collision probability of a dissimilar pair to increase it for a similar pair. However we are able to prove something stronger by examining the pair whose collisions must be sacrificed. To increase the collision probability for one pair above its J P, you must always sacrifice collisions on a pair with even higher J P. To get better recall on one pair, you must always give up recall on an even better pair. The Jaccard index itself is optimal on uniform distributions, and the short proof is a model for the general case. Proof. Let Z be the symmetric difference of X and To prove the same claim on J P for all distributions, we need a few tools. We can rearrange J P to separate the two iteration indices within the max. To analyze J P we will need to refer to each term in its outer sum via subscript. DISPLAYFORM0 DISPLAYFORM1 DISPLAYFORM2 We will also use quantifiers in this subscript to indicate a partial sum, i.e. DISPLAYFORM3 Consider two linear combinations of these distributions with coefficients α and β. DISPLAYFORM0 Proof. For a given i, if every max chooses the same side, J P (x, y) i = min(x i, y i). x i /y i has both an arg max and arg min for which this is true. This gives us part 1. DISPLAYFORM1 Thus, we can form x ′, y ′ by merging the mass of x k, x l into one element and y k, y l into one element and have DISPLAYFORM2.. w n be distributions with disjoint support, and consider J P (α · w, β · w). Repeat the merging process until all elements of each w i are merged into one. This gives us.This ordering also lets us work more effectively with the z distributions we constructed in theorem II.1. This lemma contains all the algebra needed for the main proof. DISPLAYFORM3 Working first with the lower group of indices, FIG2, the green regions of x and y could be shifted to overlap more and improve the collision probability of this pair, but any modification that achieves that would worsen at least one of collision probabilities between x or y and z green (both of which have higher collision probability than the (x, y) pair.) DISPLAYFORM4 And since DISPLAYFORM5 we also know that DISPLAYFORM6 which gives us part 1. Now, continuing on to the upper group of indices, DISPLAYFORM7. By noting that the choices within the max are preserved, we conclude part 3. Finally, having bounded all indices, part 4: DISPLAYFORM8 We now have the tools to prove the optimality of J P. DISPLAYFORM9. This implies that no method of sampling from discrete distributions can have collision probabilities dominating J P. J P is Pareto optimal. Furthermore, J P (x, z) ≥ J P (x, y) and J P (y, z) ≥ J P (x, y). To exceed J P (x, y), G must sacrifice at least one pair that is closer under J P than (x, y).Proof. Let m be the number of elements i for which J P (x, y) i < min(x i, y i).In the base case where m = 0, J P (x, y) = ∑ i min(x i, y i) which cannot be improved. Assume the proposition to be proved is true ∀x, ∀y, and ∀p < m. By IV.2.1 we know that m ≤ n − 2, since at least the two endpoints of the sorted list have reached their upper bound. We proceed by induction on m. As in theorem II.1, for each a where DISPLAYFORM10. Reorder i according to the sorting of such that DISPLAYFORM11 ). Now consider a new sampling method, G. The following events are pairwise disjoint by inspection. DISPLAYFORM12 So their probabilities are constrained. DISPLAYFORM13 When these probabilities are given by J P they already sum to 1. DISPLAYFORM14 must be true. Since the two cases are symmetric, we will assume the first one: DISPLAYFORM15 and we are done. Otherwise, DISPLAYFORM16, and the terms i ≥ a must compensate for the loss on the terms i < a, so Pr[DISPLAYFORM17, so these terms have exhausted z a and cannot be increased. Using IV.3.1 and IV.3 we know that this adds at least one additional term that is fully consumed, i.e. the size of {i :, z). By IV.3 we know that J P (x, z a) ≥ J P (x, y) and by induction we know J P (x, z) ≥ J P (x, z a), so we conclude J P (x, z) ≥ J P (x, y) and symmetrically J P (y, z) ≥ J P (x, y). FIG5 shows the mechanism of the proof intuitively. On three element distributions, two of the elements are fully constrained, in this case the blue and red terms, so we construct our adversarial z around green. On three elements, no induction is necessary and the diagram itself proves the relationship. DISPLAYFORM18 Since J P is only Pareto optimal, we should be able to find a sampling method that exceeds it for some pairs but is below it on others. We can generalize our algorithm to construct such a method. Consider arranging the elements of the state space as the leaves of a tree. Internal nodes in the tree are given the weight of the sum of their children, and each is assigned its own exponential hash. Perform P-MinHash among the children of the root node. If the selected node is a leaf, emit it as the sample. If it is an internal node, recurse and repeat. In this generalization, our original algorithm is represented by making all elements direct children of the root. We can prioritize collisions on an index i by placing it closer to the root node than all others; in particular, the tree FIG2 ). Since i and its internal-node sibling form a two element distribution, by lemma 3 the probability of a collision on i will be min(x i, y i) for all x, y. What of J W then? Is it another Pareto optimum? It is not, it is dominated by J P. Proof. The lower bound becomes clear by rewriting J W in a similar form to J P. DISPLAYFORM19 To achieve this lower bound, we can transform the distributions by moving the "excess" mass to new elements. DISPLAYFORM20 Shifting the mass in this way has no effect on J W, but it decreases J P to equal J W. x ′ and y ′ can be expressed as linear combinations over the three sets of indices, so using lemma IV.2.2, DISPLAYFORM21 To achieve the upper bound, 1 − p, consider inverting this transformation, reallocating the p extra mass to maximize J P (x ′′, y ′′) while holding J W (x ′′, y ′′) constant. To avoid increasing J W, we must add the mass to disjoint elements, so divide the indices into two sets, X, Y. We find that if we distribute the mass proportional to the original value, J P reaches the total variation limit regardless of the choice of X, Y. Let |X|:= ∑ i∈X min(x i, y i). DISPLAYFORM22 We can express this as a linear combination of two distributions with disjoint support. DISPLAYFORM23 Since p is the total variation distance of x and y, 1 − p is the maximum collision probability that is possible between two distributions in any context, so it is the upper bound here as well. This gives us some insight into how J P and J W differ. J P ranks distributions as more similar than J W if their extra mass is on elements that both distributions share. Like 1 − J W, 1 − J P is a metric on probability distributions. Theorem IV.6. 1 − J P is a proper metric on P where P(Ω) is the space of probability distributions over a finite set Ω.Proof. Symmetry is obvious. Non-degeneracy over P follows from DISPLAYFORM24 The triangle inequality follows from being a collision prob- DISPLAYFORM25 for any distribution z. But by the union bound, DISPLAYFORM0 The algorithm we've presented so far is suitable for sparse data such as documents or graphs. It is linear in the number of non-zeros, equivalent to Ioffe 2010 BID6. On dense data (such a image feature histograms BID7) there's significant overlap in the supports of each distribution, so rehashing each element for every distribution wastes Algorithm 2: Dense and Continuous P-MinHash. "Global-Bound" A* Sampling BID9 with a fixed seed.input: sample space Ω, sigma-finite measure µ, proposal sigma-finite measure λ, finite upper bound, B:= max(µ(i)/λ(i)) shared random seed s output: Stable sample from (Ω, F, µ) DISPLAYFORM0 work. With a shared stream of sorted hashes, we expect the hash we select for each distribution to be biased towards the beginning of the stream, and closer to the beginning when the data is denser. Therefore one might expect that we could improve performance by searching only some prefix of the stream to find our sample. A* Sampling (Maddison, Tarlow, and Minka 2014 BID9) explores this idea thoroughly, and we lean on it heavily in this section. In particular, we use their "Global Bound" algorithm, and essentially just run it with a fixed random seed (algorithm 2.) We leave the proof of running time and of correctness as a sampling method to that work, and limit our discussion to the proof of the ing collision probability. (Their derivation uses Gumbel variables and maxima. We use exponential variables and minima to make the continuity with the rest of our work clear, which is achieved by a simple change of variables.)The key insight of A* Sampling is that when a (possibly infinite) stream of independent exponential random variables is ordered, the exact marginal distributions of each variable can be computed as a function of the rank and the previous variables. If the vector of sorted exponential variables is e with corresponding parameter vector x, then once e 1,..., e k−1 are all known, the distribution of e k is a truncated exponential with rate |x \ ∪ k i=1 x i | truncated from below at e k−1. The "statelessness" of exponential variables makes this truncation easy to accomplish. Simply generate the desired exponential, and add e k−1 to shift it. To change the parameters of those exponentials and find the new minimum element, only a small prefix of the list must be examined. The running time of finding the new minimum is a function of the difference between the two vectors of parameters, and has equivalent running time to rejection sampling. This gives algorithm 2 equivalent running time to the state of the art for computing a MinHash of dense data, Shrivastava 2016 BID7. Because algorithm 2 admits an unbounded list of random variables, it is also applicable to continuous distributions, as the paper BID9 describes in detail. Let's first show that algorithm 2 gives J P (µ, ν) when Ω is finite. Indeed, this construction is simply an alternative way of finding the minimum − log U i /µ i. A* sampling merely reads off the minimum of − log U i /λ i * λ i /µ i = − log U i /µ i, and similar for ν. To prove the general case for infinite Ω, we need to first define what we mean by J P (µ, ν) in that setting. One option is to replace all the summation by integrals in the formula BID1. This runs into two difficulties however: 1) A probability space Ω may not be a subset of R n. 2) Either µ or ν could be singular. Instead, we define it as a limit over increasingly finer finite partitions of Ω. More formally, Definition V.1. Assume J P (µ, ν) is defined as before when |Ω| < ∞, we define DISPLAYFORM1 where F ranges over finite partitions of the space Ω, and µ F denotes the push-forward of µ with respect to the map π: Ω → F, π(x) = Q ∈ F iff x ∈ Q. (µ F is simply a coarsified probability measure on the finite space F where it (tautologically) assigns probability µ(Q) to the element Q ∈ F.)First we verify that the definition above coincides with J P when Ω is finite. DISPLAYFORM2, with strict inequality if both µ, ν have nonzero masses on those two elements. Proof. By considering J P (µ, ν) as the probability that the argmin's of two lists of independent exponentials land on the same index 1 ≤ i * ≤ n, and using the fact that the minimum of two independent exponentials is an exponential with the sum of the parameters, we can couple the four argmin's arising from µ, ν, µ ′, ν ′ and conclude by inspection. The lemma shows that any partition of Ω will lead to a J P that's greater than or equal to the original J P. So the infimum is achieved with the most refined partition, namely Ω itself. J P Jensen−Shannon Divergence DISPLAYFORM3 count Fig. 3: The Jensen-Shannon divergence compared with J P and J W on pairs of normalized unigram term vectors of random web documents. JSD has a much tighter relationship with J P than J W. We show exact bounds for J W against JSD and approximate bounds for J P against JSD where DISPLAYFORM4 The curve that appears to lower bound JSD against J P is violated on 10 −7 of the pairs. The left graph shows the joint distribution of J P and J W and the bounds we prove. The right graph shows the conditional distributions of J P and J W against the (set) Jaccard index of the terms. We show the distribution of the log of their ratios against the Jaccard index and highlight the median. J P is generally centered around the Jaccard index, while J W is consistently centered below, as predicted by their behavior on uniform distributions. Finally we show that A* sampling applied to µ and ν simultaneously has a collision probability equal to J P (µ, ν) as defined above. Theorem V.3. Given two probability measures µ and ν on an arbitrary Polish space Ω, both absolutely continuous with respect to a common third measure λ, it is possible to apply A* sampling with base distribution λ, either in-order or with a hierarchical partition of Ω, to sample from µ and ν simultaneously. Further, the probability of the procedure terminating at the same point p ∈ Ω for both µ and ν is exactly J P (µ, ν).Proof. The first statement follow from the procedural definition of A* sampling described in BID9. For the second statement, since in-order A* is proven equivalent to hierarchical partition A* in BID9, we are free to choose any partition to our convenience. The natural choice is then the partition used in the definition of J P (µ, ν).More precisely, we know there is a finite partition: Precision/recall curves illustrating the typical case of retrieval using a key-value store. Each point represents outputting o independent sums of a hashes each, for a collision probability of 1 DISPLAYFORM5 DISPLAYFORM6 The cost in storage and CPU is dominated by o, so we connect these points to show the trade-offs possible at similar cost.representative of each part Q ∈ F are jointly distributed as exponentials with rate λ(Q) with common seed. Let U, V be the two coupled A* processes restricted to F. Either one of them does not terminate, or they both terminate and collide conditionally with probability J P (µ F, ν F). In other words, letting P(T) be the probability that both terminate at F level, and AC(U, V ; F) be the collision probability of U, V restricted to F, then AC(U, V ; DISPLAYFORM7 So AC(U, V ; F) is squeezed between J P (µ, ν)P(T) and J P (µ, ν) + ϵ. Since P(T) → 1 and ϵ → 0 under a refinement sequence F, we get in the limit DISPLAYFORM8 To determine whether the difference between J P and J W matters in practice and whether achieving J P as a collision probability is a useful goal, we computed both for a large sample of pairs of unigram term vectors of web documents. From an index of 6.6 billion documents, we selected pairs using a sum of 2 W-MinHashes to perform importance sampling. We computed several similarity scores for 100 million pairs of normalized unigram term vectors, and weighted them by the inverse of their sampling probability to simulate an unbiased sample of all non-zero pairs. The Jensen-Shannon divergence (JSD) defines the information loss that from representing two distributions using a model that is an equal mixture of them, and as such is the ideal criterion to form informationpreserving clusters of items of equal importance. Like both J W and J P it is bounded, symmetric, and monotonic in a metric distance. Due to these properties, as well as its popularity, we use it here as a basis for comparison. J P has a much tighter relationship with the JensenShannon divergence than J W as shown in figure 3. Tight bounds on JSD as a function of J W are given by J W's monotonic relationship with Total Variation, as described by BID10. Let p be the total variation distance, and d(p) = = p extends these to J W. These same bounds apply to J P as well, but J P has a much tighter relationship with what appears to be a much higher lower bound. We have approached finding this lower bound with large differential equations that we have only solved numerically, but small examples form good approximate bounds. On 2 element distributions, JSD has a direct relationship with J P, and only 1 × 10 −7 of the pairs fall below the ing curve, d(1 − J P). No pairs in our sample had JSD more than 0.0077 below it. In contrast, J W puts 7 × 10 −3 of the pairs below this curve, with the farthest point 0.16 below. We also compare both J P and J W to the Jaccard index of the set of terms, and compute a kernel density estimate of the log of their ratios. J P is generally centered around the Jaccard index, while J W is consistently centered below, as predicted by their behavior on uniform distributions. This makes P-MinHash less disruptive as a drop-in replacement for an unweighted MinHash. Parameters of the system such as the number of hashes or the length of concatenated hashes are likely to continue to function well. In the typical case of retrieval using a key-value store, performance is characterized by cheap ANDs and expensive ORs. BID11 To reduce the collision probability we can sum multiple hashes to form keys, but to raise the collision probability, we must output multiple independent keys. This lets us apply an asymmetric sigmoid to the collision probabilities, DISPLAYFORM0 with o independent outputs of a summed hashes each. Assuming that the cost of looking up hashes dominates the cost of generating them, ANDs are essentially free, while CPU and storage cost are both linear in the number of ORs. Furthermore, as a increases linearly, o must increase exponentially to keep the inflection point of the sigmoid in the same place. For instance, if the sigmoid passes through (0.5, 0.5), then o ≈ log2 a. This gives a significant performance advantage to algorithms with higher collision probabilities, and thus to J P over J W. Lowering the probability is much cheaper than raising it. The effect of this is demonstrated in FIG6. Unsurprisingly from the tightness of the joint distribution, PMinHash achieves better precision and recall retrieving low JSD documents for a given cost. More surprising is that it also achieves slightly better precision and recall on retrieving high J W documents when the cost is low, even though this is the task W-MinHashes are designed for. The reason for this can be seen from the upper bound, J P ≤ 2J W /(1 + J W). On items that achieve this bound, the collision probability when summing two hashes, (2x/(1+x)) 2, is similar to 4x 2 near 0 and similar to x near 1. This in effect gives it the recall of 1 hash with the precision of 2 hashes on this subset of items, and thus a better precision/recall trade-off overall. We've described a new generalization of the Jaccard index, and shown several qualities that motivate it as the natural extension to probability distributions. In particular, we proved that it is optimal on all distributions in the same sense that the Jaccard index is optimal on uniform distributions. We've demonstrated its utility by showing J P's similarity in practice to the Jensen-Shannon divergence, a popular clustering criterion. We've described two MinHashing algorithms that achieve this as their collision probability with equivalent running time to the state of the art on both sparse and dense data.
The minimum of a set of exponentially distributed hashes has a very useful collision probability that generalizes the Jaccard Index to probability distributions.