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

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

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

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

где x, y - строки в матрице m.

Я экспериментировал с различными алгоритмами вычисления матрицы 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 и apply:

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
4 голоса | спросил Andres Kull 14 MarpmSat, 14 Mar 2015 18:03:12 +03002015-03-14T18:03:12+03:0006 2015, 18:03:12

2 ответа


0

Это моя реализация вашего 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,25 раза быстрее

> 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) кластеризации kmeans всех 50 000 компаний для идентификации, скажем, 1000 центроидов, и (3) с использованием интерполированных измерения и расстояние Дженсена-Шеннона между ними для кластеризации.

ответил Martin Morgan 15 MaramSun, 15 Mar 2015 04:28:19 +03002015-03-15T04:28:19+03:0004 2015, 04:28:19
0

Я уверен, что есть лучшие подходы, чем следующие, но сама ваша функция 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
ответил nrussell 14 MarpmSat, 14 Mar 2015 18:48:33 +03002015-03-14T18:48:33+03:0006 2015, 18:48:33

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132