加速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