Какая реализация R дает самое быстрое вычисление матрицы JSD?

Question

JSD-матрица представляет собой матрицу подобия распределений на основе Jensen-Shannon divergence. В данной матрице m, в строках которой представлены распределения, мы хотели бы найти расстояние JSD между каждым распределением. Результирующая JSD-матрица представляет собой квадратную матрицу с размерами nrow (m) x nrow (m). Это треугольная матрица, в которой каждый элемент содержит значение JSD между двумя строками в m.

JSD можно рассчитать по следующей функции R:

JSD<- function(x,y) sqrt(0.5 * (sum(x*log(x/((x+y)/2))) + sum(y*log(y/((x+y)/2))))) 

, где х, у строк в матрице м.

Я экспериментировал с различными алгоритмами вычисления матрицы JSD в R, чтобы выяснить самый быстрый. Для моего удивления алгоритм с двумя вложенными циклами работает быстрее, чем разные векторизованные версии (распараллеливаются или нет). Я не доволен результатами. Не могли бы вы определить лучшие решения, чем те, в которых я играю?

library(parallel) 
library(plyr) 
library(doParallel) 
library(foreach) 

nodes <- detectCores() 
cl <- makeCluster(4) 
registerDoParallel(cl) 

m <- runif(24000, min = 0, max = 1) 
m <- matrix(m, 24, 1000) 
prob_dist <- function(x) t(apply(x, 1, prop.table)) 
JSD<- function(x,y) sqrt(0.5 * (sum(x*log(x/((x+y)/2))) + sum(y*log(y/((x+y)/2))))) 
m <- t(prob_dist(m)) 
m[m==0] <- 0.000001 

Алгоритм с двумя вложенными циклами:

dist.JSD_2 <- function(inMatrix) { 
    matrixColSize <- ncol(inMatrix) 
    resultsMatrix <- matrix(0, matrixColSize, matrixColSize) 
    for(i in 2:matrixColSize) { 
    for(j in 1:(i-1)) { 
     resultsMatrix[i,j]=JSD(inMatrix[,i], inMatrix[,j]) 
    } 
    } 
    return(resultsMatrix) 
} 

Алгоритм с наружным:

dist.JSD_3 <- function(inMatrix) { 
    matrixColSize <- ncol(inMatrix) 
    resultsMatrix <- outer(1:matrixColSize,1:matrixColSize, FUN = Vectorize(function(i,j) JSD(inMatrix[,i], inMatrix[,j]))) 
    return(resultsMatrix) 
} 

Алгоритм с combn и применения:

dist.JSD_4 <- function(inMatrix) { 
    matrixColSize <- ncol(inMatrix) 
    ind <- combn(matrixColSize, 2) 
    out <- apply(ind, 2, function(x) JSD(inMatrix[,x[1]], inMatrix[,x[2]])) 
    a <- rbind(ind, out) 
    resultsMatrix <- sparseMatrix(a[1,], a[2,], x=a[3,], dims=c(matrixColSize, matrixColSize)) 
    return(resultsMatrix) 
} 

Алгоритм с combn и aaply:

dist.JSD_5 <- function(inMatrix) { 
    matrixColSize <- ncol(inMatrix) 
    ind <- combn(matrixColSize, 2) 
    out <- aaply(ind, 2, function(x) JSD(inMatrix[,x[1]], inMatrix[,x[2]])) 
    a <- rbind(ind, out) 
    resultsMatrix <- sparseMatrix(a[1,], a[2,], x=a[3,], dims=c(matrixColSize, matrixColSize)) 
    return(resultsMatrix) 
} 

тест производительности:

mbm = microbenchmark(
    two_loops = dist.JSD_2(m), 
    outer = dist.JSD_3(m), 
    combn_apply = dist.JSD_4(m), 
    combn_aaply = dist.JSD_5(m), 
    times = 10 
) 
ggplot2::autoplot(mbm) 

> summary(mbm) 
     expr  min  lq  mean median 
1 two_loops 18.30857 18.68309 23.50231 18.77303 
2  outer 38.93112 40.98369 42.44783 42.16858 
3 combn_apply 20.45740 20.90747 21.49122 21.35042 
4 combn_aaply 55.61176 56.77545 59.37358 58.93953 
     uq  max neval cld 
1 18.87891 65.34197 10 a 
2 42.85978 48.82437 10 b 
3 22.06277 22.98803 10 a 
4 62.26417 64.77407 10 c 

Answer 1

Это моя реализация ваших dist.JSD_2

dist0 <- function(m) { 
    ncol <- ncol(m) 
    result <- matrix(0, ncol, ncol) 
    for (i in 2:ncol) { 
     for (j in 1:(i-1)) { 
      x <- m[,i]; y <- m[,j] 
      result[i, j] <- 
       sqrt(0.5 * (sum(x * log(x/((x + y)/2))) + 
          sum(y * log(y/((x + y)/2))))) 
     } 
    } 
    result 
} 

Обычных шагов заменить итеративные вычисления с векторизацией версий. Я переместил sqrt(0.5 * ...) изнутри петель, где он был применен к каждому элементу result, вне цикла, где он применяется к вектору result.

Я понял, что sum(x * log(x/(x + y)/2)) может быть написана как sum(x * log(2 * x)) - sum(x * log(x + y)). Первая сумма вычисляется один раз для каждой записи, но может быть рассчитана один раз для каждого столбца. Он также выходит из циклов с вектором значений (по одному элементу для каждого столбца), вычисленным как colSums(m * log(2 * m)).

Остающийся член внутри внутреннего контура sum((x + y) * log(x + y)). Для заданного значения i, мы можем разменять пространство для скорости по векторизации это во всех соответствующих y столбцов в качестве матричной операции

j <- seq_len(i - 1L) 
xy <- m[, i] + m[, j, drop=FALSE] 
xylogxy[i, j] <- colSums(xy * log(xy)) 

Конечный результат

dist4 <- function(m) { 
    ncol <- ncol(m) 

    xlogx <- matrix(colSums(m * log(2 * m)), ncol, ncol) 
    xlogx2 <- xlogx + t(xlogx) 
    xlogx2[upper.tri(xlogx2, diag=TRUE)] <- 0 

    xylogxy <- matrix(0, ncol, ncol) 
    for (i in seq_len(ncol)[-1]) { 
     j <- seq_len(i - 1L) 
     xy <- m[, i] + m[, j, drop=FALSE] 
     xylogxy[i, j] <- colSums(xy * log(xy)) 
    } 

    sqrt(0.5 * (xlogx2 - xylogxy)) 
} 

Что дает результаты, численно равны (хотя и не совсем идентичны) к исходному

> all.equal(dist0(m), dist4(m)) 
[1] TRUE 

и о 2.25x быстрее

> microbenchmark(dist0(m), dist4(m), dist.JSD_cpp2(m), times=10) 
Unit: milliseconds 
      expr  min  lq  mean median  uq  max neval 
     dist0(m) 48.41173 48.42569 49.26072 48.68485 49.48116 51.64566 10 
     dist4(m) 20.80612 20.90934 21.34555 21.09163 21.96782 22.32984 10 
dist.JSD_cpp2(m) 28.95351 29.11406 29.43474 29.23469 29.78149 30.37043 10 

Вы все равно будете ждать около 10 часов, хотя это, по-видимому, означает очень большую проблему. Алгоритм кажется квадратным по числу столбцов, но количество столбцов здесь мало (24) по сравнению с количеством строк, поэтому мне интересно, каков фактический размер обрабатываемых данных? Расчёты ncol * (ncol - 1)/2 для расчета.

Сырой подход к дальнейшему увеличению производительности параллельной оценки, которая следующие орудия, используя parallel::mclapply()

dist4p <- function(m, ..., mc.cores=detectCores()) { 
    ncol <- ncol(m) 

    xlogx <- matrix(colSums(m * log(2 * m)), ncol, ncol) 
    xlogx2 <- xlogx + t(xlogx) 
    xlogx2[upper.tri(xlogx2, diag=TRUE)] <- 0 

    xx <- mclapply(seq_len(ncol)[-1], function(i, m) { 
     j <- seq_len(i - 1L) 
     xy <- m[, i] + m[, j, drop=FALSE] 
     colSums(xy * log(xy)) 
    }, m, ..., mc.cores=mc.cores) 

    xylogxy <- matrix(0, ncol, ncol) 
    xylogxy[upper.tri(xylogxy, diag=FALSE)] <- unlist(xx) 

    sqrt(0.5 * (xlogx2 - t(xylogxy))) 
} 

Мой ноутбук имеет 8 номинальных ядер, и 1000 столбцов У меня

> system.time(xx <- dist4p(m1000)) 
    user system elapsed 
48.909 1.939 8.043 

предполагает, что Я получаю 48 секунд процессорного времени в 8 с тактового времени. Алгоритм по-прежнему квадратичен, поэтому это может сократить общее время вычисления до 1 часа для полной проблемы. Память может стать проблемой на многоядерной машине, где все процессы конкурируют за один и тот же пул памяти; может потребоваться выбрать mc.cores меньше, чем доступно.

С большим ncol, чтобы получить лучшую производительность, следует избегать вычисления полного набора расстояний. В зависимости от характера данных может иметь смысл фильтровать дубликаты столбцов или фильтровать информационные столбцы (например, с наибольшей дисперсией) или ...Соответствующая стратегия требует дополнительной информации о том, что представляют собой столбцы и какова цель для матрицы расстояния. Вопрос: «Насколько похожа компания i на другие компании?» можно ответить, не вычисляя полную матрицу расстояний, всего одну строку, поэтому, если количество вопросов, заданных по отношению к общему числу компаний, невелико, то, возможно, нет необходимости вычислять полную матрицу расстояний? Другая стратегия может заключаться в сокращении числа компаний, подлежащих кластеризации, посредством (1) упрощения 1000 рядов измерений с использованием анализа основных компонентов, (2) кменизация кластеризации всех 50 тыс. Компаний для определения 1000 центроидов и (3) использование интерполированных измерений и расстояния Дженсена-Шэннона между ними для кластеризации.

Answer 2

Я уверен, что есть более эффективные подходы, чем ниже, но ваша сама JSD функция тривиальным может быть преобразована в функцию Rcpp просто заменив sum и log на их эквиваленты сахара Rcpp и используя std::sqrt в место R's base::sqrt.

#include <Rcpp.h> 

// [[Rcpp::export]] 
double cppJSD(const Rcpp::NumericVector& x, const Rcpp::NumericVector& y) { 
    return std::sqrt(0.5 * (Rcpp::sum(x * Rcpp::log(x/((x+y)/2))) + 
          Rcpp::sum(y * Rcpp::log(y/((x+y)/2))))); 
} 

---

Я тестировал только с dist.JST_2 подходом (так как это был самый быстрый вариант), но вы должны увидеть улучшение при использовании cppJSD вместо JSD независимо от реализации:

R> microbenchmark::microbenchmark(
    two_loops = dist.JSD_2(m), 
    cpp = dist.JSD_cpp(m), 
    times=100L) 
Unit: milliseconds 
     expr  min  lq  mean median  uq  max neval 
two_loops 41.25142 41.34755 42.75926 41.45956 43.67520 49.54250 100 
     cpp 36.41571 36.52887 37.49132 36.60846 36.98887 50.91866 100 

---

EDIT: На самом деле, ваш dist.JSD_2 функционировать сам по себе может быть легко преобразована в функцию Rcpp за дополнительное ускорение:

// [[Rcpp::export("dist.JSD_cpp2")]] 
Rcpp::NumericMatrix foo(const Rcpp::NumericMatrix& inMatrix) { 
    size_t cols = inMatrix.ncol(); 
    Rcpp::NumericMatrix result(cols, cols); 
    for (size_t i = 1; i < cols; i++) { 
    for (size_t j = 0; j < i; j++) { 
     result(i,j) = cppJSD(inMatrix(Rcpp::_, i), inMatrix(Rcpp::_, j)); 
    } 
    } 
    return result; 
} 

(где cppJSD был определен в том же файл .cpp, как указано выше).Вот тайминги:

R> microbenchmark::microbenchmark(
    two_loops = dist.JSD_2(m), 
    partial_cpp = dist.JSD_cpp(m), 
    full_cpp = dist.JSD_cpp2(m), 
    times=100L) 
Unit: milliseconds 
     expr  min  lq  mean median  uq  max neval 
    two_loops 41.25879 41.36729 42.95183 41.84999 44.08793 54.54610 100 
partial_cpp 36.45802 36.62463 37.69742 36.99679 37.96572 44.26446 100 
    full_cpp 32.00263 32.12584 32.82785 32.20261 32.63554 38.88611 100 

---

dist.JSD_2 <- function(inMatrix) { 
    matrixColSize <- ncol(inMatrix) 
    resultsMatrix <- matrix(0, matrixColSize, matrixColSize) 
    for(i in 2:matrixColSize) { 
    for(j in 1:(i-1)) { 
     resultsMatrix[i,j]=JSD(inMatrix[,i], inMatrix[,j]) 
    } 
    } 
    return(resultsMatrix) 
} 
## 
dist.JSD_cpp <- function(inMatrix) { 
    matrixColSize <- ncol(inMatrix) 
    resultsMatrix <- matrix(0, matrixColSize, matrixColSize) 
    for(i in 2:matrixColSize) { 
    for(j in 1:(i-1)) { 
     resultsMatrix[i,j]=cppJSD(inMatrix[,i], inMatrix[,j]) 
    } 
    } 
    return(resultsMatrix) 
} 

---

m <- runif(24000, min = 0, max = 1) 
m <- matrix(m, 24, 1000) 
prob_dist <- function(x) t(apply(x, 1, prop.table)) 
JSD <- function(x,y) sqrt(0.5 * (sum(x*log(x/((x+y)/2))) + sum(y*log(y/((x+y)/2))))) 
m <- t(prob_dist(m)) 
m[m==0] <- 0.000001