Set and Cluster Similarity in R
The following set similarity functions are written in C++ for easy integration
with R using the sourceCpp command.
The first function computes the Jaccard index assuming inputs are unique, sorted integers where each integer denotes an element of a set. L1 and L2 are both sets.
/**
* The Jaccard Similarity Coefficient or Jaccard Index is used to compare the
* similarity/diversity of sample sets. It is defined as the size of the
* intersection of the sets divided by the size of the union of the sets.
*
* INPUTS:
* L1: A list of unique, sorted integers denoting set elements.
* L2: A list of unique, sorted integers denoting set elements.
*
* RETURNS:
* Size(Intersection(L1,L2))/Size(Union(L1,L2))
*
* SIDE-EFFECTS:
* None: this is the 21st-century.
*
* COMPLEXITY:
* Time: O(N1+N2)
* Space: O(1)
*/
// [[Rcpp::export]]
NumericVector JaccardSimilarity(const NumericVector L1, const NumericVector L2){
int n1 = L1.size(); //Number of elmeents in set 1
int n2 = L2.size(); //Number of elmeents in set 2
int i1 = 0; //Current element of interest in set 1
int i2 = 0; //Current element of interest in set 2
int size_intersection = 0; //Result: to be incremented in algorithm
int size_union = 0; //Result: to be incremented in algorithm
while(true){
if(i1==n1){
size_union += n2-i2;
break;
} else if (i2==n2){
size_union += n1-i1;
break;
}
if(L1[i1]==L2[i2]){ //Element is in both sets.
i1++; //Therefore, we discard both the current element
i2++; //from both sets and mark the current element as
size_union++; //being part of both the union and the
size_intersection++; //intersection
} else if(L1[i1]<L2[i2]){ //Element is only in the first set.
i1++; //Throw it out
size_union++; //and mark it as a member of only the union
} else if(L1[i1]>L2[i2]){ //Element is only in the second set.
i2++; //Throw it out
size_union++; //and mark it as a member of only the union
}
}
return NumericVector::create((double)size_intersection/(double)size_union);
}
The second function computes the Jaccard index assuming only that L1 and L2 are both sets whose integers indicate set members.
/**
* The Jaccard Similarity Coefficient or Jaccard Index is used to compare the
* similarity/diversity of sample sets. It is defined as the size of the
* intersection of the sets divided by the size of the union of the sets.
*
* INPUTS:
* L1: A list of integers denoting set elements.
* L2: A list of integers denoting set elements.
*
* RETURNS:
* Size(Intersection(L1,L2))/Size(Union(L1,L2))
*
* SIDE-EFFECTS:
* None: this is the 21st-century.
*
* COMPLEXITY:
* Time: O(N1+N2)
* Space: O(N1+N2) worst-case
*/
#include <unordered_map>
// [[Rcpp::export]]
NumericVector JaccardSimilarityUnsorted(const NumericVector L1, const NumericVector L2){
int n1 = L1.size(); //Number of elmeents in set 1
int n2 = L2.size(); //Number of elmeents in set 2
int size_intersection = 0; //Result: to be incremented in algorithm
int size_union = 0; //Result: to be incremented in algorithm
std::unordered_map<int, unsigned char> seen;
for(int i1=0;i1<n1;i1++){ //Loop through all elements in Set 1
if(!seen.count(L1[i1])){ //If we haven't seen the element before
seen[L1[i1]] = 1; //Mark it as having been seen in Set 1
size_union++; //And as being part of the union
}
}
for(int i2=0;i2<n2;i2++){ //Loop through all elements in Set 2
if(!seen.count(L2[i2])){ //If we have not seen the element before
size_union++; //then it is part of the union.
seen[L2[i2]] = 2; //Mark it as being done.
} else if(seen[L2[i2]]==1){ //Elsewise, we've seen the element in the
size_intersection++; //Set 1, it is part of the intersection
//But not the union, because that would
//be double counting.
seen[L2[i2]] = 2; //Mark it as being done.
}
}
return NumericVector::create((double)size_intersection/(double)size_union);
}
This function, due to Jamie Murdoch, quickly computes the Fowlkes-Mallows index of a cluster assignment.
/**
* Third iteration of the Fowlkes and Mallows similarity in C++.
* For each cluster i, track the number of elements in that cluster.
* For each pair of clusters (x, y), track the number of elements assigned to x
* in L1 and to y in L2. This is the "overlap" between x and y.
* <C1, C1> is the sum of squares of cluster sizes. To see this, note that
* for all pairs of elements i, j in cluster x, c_{ij} = 1. For all pairs
* of elements in different clusters, c_{ij} = 0. Then, for each cluster x,
* the number of elements with value 1 in C_{ij} is |x|^2, the square of the
* cluster size.
* Another way to see this is that C_{ij} is an adjacency matrix for a graph, in
* which the clusters represent separate cliques. We want to count the number
* of edges for each clique, counting self-edges once and undirected edges
* between elements twice.
* <C1, C2> is the sum of squares of overlap sizes. To see this, consider the
* overlap O_{xy} between clusters x and y. This is the number of elements that
* belong to x in L1 and to y in L2. For any pair (i,j) of these elements,
* C_{ij} = 1 in both C1 and C2. This ordered pair contributes 1 to <C1, C2>.
* We sum these contributions across all pairs in the overlap (of which there
* are |O_{xy}|^2) and all overlaps.
*
* INPUTS:
* L1: Each item of the list is a point. The number stored at the element
* is the cluster the point is assigned to.
* L2: The same as L1. Must be the same length as L1.
*
* RETURNS:
* The Fowlkes and Mallows index
*
* SIDE-EFFECTS:
* None
*
* COMPLEXITY:
* Time: O(K^2+n), where K = number of clusters
* Space: O(K^2)
*
* SOURCES:
* Asa Ben-Hur, Andre Elisseeff, and Isabelle Guyon (2001) A stability based
* method for discovering structure in clustered data. Biocomputing 2002: pp.
* 6-17.
*/
// [[Rcpp::export]]
NumericVector FowlkesMallows(const NumericVector L1, const NumericVector L2){
int n = L1.size();
int K = max(L1);
ind overlaps[K][K];
ind cluster_sizes1[K], cluster_sizes2[K];
for(int i = 0; i < K; i++){
cluster_sizes1[i] = 0;
cluster_sizes2[i] = 0;
for(int j = 0; j < K; j++)
overlaps[i][j] = 0;
}
//O(n) time. O(K^2) space. Determine the size of each cluster as well as the
//size of the overlaps between the clusters.
for(int i = 0; i < n; i++){
cluster_sizes1[(int)L1[i] - 1]++; // -1's account for zero-based indexing
cluster_sizes2[(int)L2[i] - 1]++;
overlaps[(int)L1[i] - 1][(int)L2[i] - 1]++;
}
// O(K^2) time. O(1) space. Square the overlap values.
int C1dotC2 = 0;
for(int j = 0; j < K; j++){
for(int k = 0; k < K; k++){
C1dotC2 += pow(overlaps[j][k], 2);
}
}
// O(K) time. O(1) space. Square the cluster sizes
int C1dotC1 = 0, C2dotC2 = 0;
for(int i = 0; i < K; i++){
C1dotC1 += pow(cluster_sizes1[i], 2);
C2dotC2 += pow(cluster_sizes2[i], 2);
}
return NumericVector::create((double)C1dotC2 / (sqrt(C1dotC1) * sqrt(C2dotC2)));
}
We can modify the above to quickly compute the cluster assignment similarity using the Jaccard index, as suggested by Ben-Hur, Elisseeff, and Guyon (2002).
/**
* The Jaccard Similarity Coefficient or Jaccard Index is used to compare the
* similarity/diversity of sample sets. It is defined as the size of the
* intersection of the sets divided by the size of the union of the sets. Here,
* it is used to determine how similar to clustering assignments are.
*
* INPUTS:
* L1: A list. Each element of the list is a number indicating the cluster
* assignment of that number.
* L2: The same as L1. Must be the same length as L1.
*
* RETURNS:
* The Jaccard Similarity Index
*
* SIDE-EFFECTS:
* None
*
* COMPLEXITY:
* Time: O(K^2+n), where K = number of clusters
* Space: O(K^2)
*
* SOURCES:
* Asa Ben-Hur, Andre Elisseeff, and Isabelle Guyon (2001) A stability based
* method for discovering structure in clustered data. Biocomputing 2002: pp.
* 6-17.
*/
// [[Rcpp::export]]
NumericVector JaccardIndex(const NumericVector L1, const NumericVector L2){
int n = L1.size();
int K = max(L1);
ind overlaps[K][K];
ind cluster_sizes1[K], cluster_sizes2[K];
for(int i = 0; i < K; i++){
cluster_sizes1[i] = 0;
cluster_sizes2[i] = 0;
for(int j = 0; j < K; j++)
overlaps[i][j] = 0;
}
//O(n) time. O(K^2) space. Determine the size of each cluster as well as the
//size of the overlaps between the clusters.
for(int i = 0; i < n; i++){
cluster_sizes1[(int)L1[i] - 1]++; // -1's account for zero-based indexing
cluster_sizes2[(int)L2[i] - 1]++;
overlaps[(int)L1[i] - 1][(int)L2[i] - 1]++;
}
// O(K^2) time. O(1) space. Square the overlap values.
int C1dotC2 = 0;
for(int j = 0; j < K; j++){
for(int k = 0; k < K; k++){
C1dotC2 += pow(overlaps[j][k], 2);
}
}
// O(K) time. O(1) space. Square the cluster sizes
int C1dotC1 = 0, C2dotC2 = 0;
for(int i = 0; i < K; i++){
C1dotC1 += pow(cluster_sizes1[i], 2);
C2dotC2 += pow(cluster_sizes2[i], 2);
}
return NumericVector::create((double)C1dotC2/(double)(C1dotC1+C2dotC2-C1dotC2));
}