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

Question

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

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

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

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

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

(1) C 实现 - 但我不知道如何执行此操作（我不是 C 程序员）如果有帮助，您能帮我解决这个 C 实现吗？

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

还有其他想法吗？

最好的问候尤尔詹

# 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

Answer 2

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

• 改进算法。在您的HellingerDistance函数中，可以看到在计算每对的距离时，您每次都重新计算平方根，这完全是浪费时间。所以这里是改进版，calculatedistances1是新的函数，它先计算x的平方根，然后用新的HellingerDistanceSqrt计算海灵格距离，可以看出新版本加速了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)