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

Question

我正在寻找一种方法来加速这个算法。

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

对于 6 种习惯中的每一种，我都创建了一个概率质量分布 (PMF)。每个习惯的 PMF 的大小(列)可能不同。有些习惯有 10 列，大约 256 列，所有这些都取决于最不习惯行为的用户。

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

我目前尝试的是通过应用多核进程来加速算法，每个进程在一个单独的核心上有 6 个习惯。可能对加速有好处的两件事

(1) C 实现 - 但我不知道该怎么做(我不是 C 程序员)如果这有帮助，你能帮我实现这个 C 实现吗？

(2) 通过单独加入表格并让所有行然后进行逐行计算来制作笛卡尔积。重点是 R 在默认情况下会给出错误，例如数据表。对此有什么建议吗？

还有其他想法吗？

最好的问候 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

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

下面是我将如何使用 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