加速R算法来计算Hellinger距离的距离矩阵

Question

我正在寻找一种加速此算法的方法。

我的情况如下。我有一个包含6个习惯的25,000个用户的数据集。我的目标是为25,000个用户开发一个分层聚类。我在一个有16个内核，128GB RAM的服务器上运行它。 我花了3周时间才为在我的服务器上使用6个内核的10,000个用户计算这个距离矩阵。你可以想象这对我的研究来说太长了。

对于6种习惯中的每一种，我都创建了概率质量分布（PMF）。每个哈比特人的PMF可能大小（列）不同。一些习惯有10列大约256，全部取决于最不友好行为的用户。

我的算法的第一步是开发一个距离矩阵。我使用Hellinger距离来计算距离，这与使用的一些包相反。 cathersian /曼哈顿。我确实需要Hellinger距离，请参阅https://en.wikipedia.org/wiki/Hellinger_distance

我目前尝试的是通过应用多核处理器加速算法，每个核心都有6种习惯。两件事情，可能是加快

（1）C实现有益的 - 但我不知道如何做到这一点（我不是一个C程序员），你能帮助我在此C实现，如果这将是有益的？ （2）通过自己加入桌子制作一个carthesian产品，并让所有的行和所有的行进行一次行计算。 R点在例如默认情况下给出了一个错误。 data.table。对此有何建议？

还有其他想法吗？

此致Jurjen

# example for 1 habit with 100 users and a PMF of 5 columns 
Habit1<-data.frame(col1=abs(rnorm(100)), 
       col2=abs(c(rnorm(20),runif(50),rep(0.4,20),sample(seq(0.01,0.99,by=0.01),10))), 
       col3=abs(c(rnorm(30),runif(30),rep(0.4,10),sample(seq(0.01,0.99,by=0.01),30))), 
       col4=abs(c(rnorm(10),runif(10),rep(0.4,20),sample(seq(0.01,0.99,by=0.01),60))), 
       col5=abs(c(rnorm(50),runif(10),rep(0.4,10),sample(seq(0.01,0.99,by=0.01),30)))) 

    # give all users a username same as rowname 
    rownames(Habit1)<- c(1:100) 

    # actual calculation 
    Result<-calculatedistances(Habit1) 



     HellingerDistance <-function(x){ 
      #takes two equal sized vectors and calculates the hellinger distance between the vectors 

      # hellinger distance function 
      return(sqrt(sum(((sqrt(x[1,]) - sqrt(x[2,]))^2)))/sqrt(2)) 

     } 


     calculatedistances <- function(x){ 
     # takes a dataframe of user IID in the first column and a set of N values per user thereafter 

     # first set all NA to 0 
     x[is.na(x)] <- 0 



     #create matrix of 2 subsets based on rownumber 
     # 1 first the diagronal with 
     D<-cbind(matrix(rep(1:nrow(x),each=2),nrow=2),combn(1:nrow(x), 2)) 

     # create a dataframe with hellinger distances 
     B <<-data.frame(first=rownames(x)[D[1,]], 
         second=rownames(x)[D[2,]], 
         distance=apply(D, 2, function(y) HellingerDistance(x[ y,])) 
     ) 


     # reshape dataframe into a matrix with users on x and y axis 
     B<<-reshape(B, direction="wide", idvar="second", timevar="first") 

     # convert wide table to distance table object 
     d <<- as.dist(B[,-1], diag = FALSE) 
     attr(d, "Labels") <- B[, 1] 
     return(d) 

     } 

Answer 1

优化代码的第一件事情是仿形。通过分析您提供的代码，似乎主要瓶颈是HellingerDistance函数。

• 改进算法。在你的HellingerDistance函数中，可以看出在计算每对距离时，你每次重新计算平方根，这是一个总的浪费时间。所以这里是改进后的版本，calculatedistances1是新功能，它首先计算出x的平方根，并用新的HellingerDistanceSqrt来计算Hellinger距离，可以看出新版本加速了40％。

• 改善数据结构。我还注意到，您原来的calulatedistance函数中的x是一个data.frame，它的重载过多，所以我通过as.matrix将它转换为矩阵，这使得代码加快了一个数量级以上。

最后，新的calculatedistances1比我的机器上的原始版本快70多倍。

# example for 1 habit with 100 users and a PMF of 5 columns 
Habit1<-data.frame(col1=abs(rnorm(100)), 
        col2=abs(c(rnorm(20),runif(50),rep(0.4,20),sample(seq(0.01,0.99,by=0.01),10))), 
        col3=abs(c(rnorm(30),runif(30),rep(0.4,10),sample(seq(0.01,0.99,by=0.01),30))), 
        col4=abs(c(rnorm(10),runif(10),rep(0.4,20),sample(seq(0.01,0.99,by=0.01),60))), 
        col5=abs(c(rnorm(50),runif(10),rep(0.4,10),sample(seq(0.01,0.99,by=0.01),30)))) 

# give all users a username same as rowname 
rownames(Habit1)<- c(1:100) 

HellingerDistance <-function(x){ 
    #takes two equal sized vectors and calculates the hellinger distance between the vectors 

    # hellinger distance function 
    return(sqrt(sum(((sqrt(x[1,]) - sqrt(x[2,]))^2)))/sqrt(2)) 

} 

HellingerDistanceSqrt <-function(sqrtx){ 
    #takes two equal sized vectors and calculates the hellinger distance between the vectors 

    # hellinger distance function 
    return(sqrt(sum(((sqrtx[1,] - sqrtx[2,])^2)))/sqrt(2)) 

} 

calculatedistances <- function(x){ 
    # takes a dataframe of user IID in the first column and a set of N values per user thereafter 

    # first set all NA to 0 
    x[is.na(x)] <- 0 



    #create matrix of 2 subsets based on rownumber 
    # 1 first the diagronal with 
    D<-cbind(matrix(rep(1:nrow(x),each=2),nrow=2),combn(1:nrow(x), 2)) 

    # create a dataframe with hellinger distances 
    B <<-data.frame(first=rownames(x)[D[1,]], 
        second=rownames(x)[D[2,]], 
        distance=apply(D, 2, function(y) HellingerDistance(x[ y,])) 
    ) 


    # reshape dataframe into a matrix with users on x and y axis 
    B<<-reshape(B, direction="wide", idvar="second", timevar="first") 

    # convert wide table to distance table object 
    d <<- as.dist(B[,-1], diag = FALSE) 
    attr(d, "Labels") <- B[, 1] 
    return(d) 

} 


calculatedistances1 <- function(x){ 
    # takes a dataframe of user IID in the first column and a set of N values per user thereafter 

    # first set all NA to 0 
    x[is.na(x)] <- 0 

    x <- sqrt(as.matrix(x)) 



    #create matrix of 2 subsets based on rownumber 
    # 1 first the diagronal with 
    D<-cbind(matrix(rep(1:nrow(x),each=2),nrow=2),combn(1:nrow(x), 2)) 

    # create a dataframe with hellinger distances 
    B <<-data.frame(first=rownames(x)[D[1,]], 
        second=rownames(x)[D[2,]], 
        distance=apply(D, 2, function(y) HellingerDistanceSqrt(x[ y,])) 
    ) 


    # reshape dataframe into a matrix with users on x and y axis 
    B<<-reshape(B, direction="wide", idvar="second", timevar="first") 

    # convert wide table to distance table object 
    d <<- as.dist(B[,-1], diag = FALSE) 
    attr(d, "Labels") <- B[, 1] 
    return(d) 

} 

# actual calculation 
system.time(Result<-calculatedistances(Habit1)) 
system.time(Result1<-calculatedistances1(Habit1)) 
identical(Result, Result1) 

Answer 2

我知道这不是一个完整的答案，但是这个建议太长了评论。

以下是我如何使用data.table来加快此过程。它的方式，这个代码仍然没有达到你要求的，也许是因为我不完全确定你想要什么，但希望这将清楚地知道如何从这里开始。

此外，你可能想看看HellingerDist{distrEx}函数来计算Hellinger距离。现在

library(data.table) 

# convert Habit1 into a data.table 
    setDT(Habit1) 

# assign ids instead of working with rownames 
    Habit1[, id := 1:100] 

# replace NAs with 0 
    for (j in seq_len(ncol(Habit1))) 
    set(Habit1, which(is.na(Habit1[[j]])),j,0) 

# convert all values to numeric 
    for (k in seq_along(Habit1)) set(Habit1, j = k, value = as.numeric(Habit1[[k]])) 


# get all possible combinations of id pairs in long format 
    D <- cbind(matrix(rep(1:nrow(Habit1),each=2),nrow=2),combn(1:nrow(Habit1), 2)) 
    D <- as.data.table(D) 
    D <- transpose(D) 


# add to this dataset the probability mass distribution (PMF) of each id V1 and V2 
# this solution dynamically adapts to number of columns in each Habit dataset 
    colnumber <- ncol(Habit1) - 1 
    cols <- paste0('i.col',1:colnumber) 

    D[Habit1, c(paste0("id1_col",1:colnumber)) := mget(cols), on=.(V1 = id)] 
    D[Habit1, c(paste0("id2_col",1:colnumber)) := mget(cols), on=.(V2 = id)] 


# [STATIC] calculate hellinger distance 
D[, H := sqrt(sum(((sqrt(c(id1_col1, id1_col2, id1_col3, id1_col4, id1_col5)) - sqrt(c(id2_col1, id2_col2, id2_col3, id2_col4, id2_col5)))^2)))/sqrt(2) , by = .(V1, V2)] 

，如果你想使这个灵活的列在每个habit数据集数：

# get names of columns 
    part1 <- names(D)[names(D) %like% "id1"] 
    part2 <- names(D)[names(D) %like% "id2"] 

# calculate distance 
    D[, H2 := sqrt(sum(((sqrt(.SD[, ..part1]) - sqrt(.SD[, ..part2]))^2)))/sqrt(2) , by = .(V1,V2) ] 

现在，更快的距离计算

# change 1st colnames to avoid conflict 
    names(D)[1:2] <- c('x', 'y') 

# [dynamic] calculate hellinger distance 
    D[melt(D, measure = patterns("^id1", "^id2"), value.name = c("v", "f"))[ 
    , sqrt(sum(((sqrt(v) - sqrt(f))^2)))/sqrt(2), by=.(x,y)], H3 := V1, on = .(x,y)] 

# same results 
#> identical(D$H, D$H2, D$H3) 
#> [1] TRUE