Как я могу векторизовать цикл for, используемый для перестановки? - PullRequest
Новая
       10

Как я могу векторизовать цикл for, используемый для перестановки?

2 голосов
/ 26 октября 2019

Я использую R для анализа и хотел бы выполнить тест перестановки. Для этого я использую цикл for, который работает довольно медленно, и я хотел бы сделать код максимально быстрым. Я думаю, что векторизация является ключевой для этого. Однако после нескольких дней попыток я все еще не нашел подходящего решения, как перекодировать это. Я был бы очень признателен за вашу помощь!

У меня есть симметричная матрица с попарно экологическими расстояниями между популяциями ("dist.mat"). Я хочу случайным образом перемешать строки и столбцы этой матрицы расстояний, чтобы сгенерировать перестановочную матрицу расстояний ("dist.mat.mix"). Затем я хотел бы сохранить верхние треугольные значения в этой переставленной матрице расстояний (размером "nr.pairs"). Этот процесс следует повторить несколько раз ("nr.runs"). Результатом должна быть матрица ("result"), содержащая переставленные верхние треугольные значения нескольких прогонов, с размерами nrow=nr.runs и ncol=nr.pairs. Ниже приведен пример кода R, который делает то, что я хочу, используя цикл for: 

# example number of populations
nr.pops <- 20

# example distance matrix
dist.mat <- as.matrix(dist(matrix(rnorm(20), nr.pops, 5)))

# example number of runs
nr.runs <- 1000

# find number of unique pairwise distances in distance matrix
nr.pairs <- nr.pops*(nr.pops-1) / 2

# start loop
result <- matrix(NA, nr.runs, nr.pairs)
for (i in 1:nr.runs) {
  mix <- sample(nr.pops, replace=FALSE)
  dist.mat.mix <- dist.mat[mix, mix]
  result[i, ] <- dist.mat.mix[upper.tri(dist.mat.mix, diag=FALSE)]
}

# inspect result
result

Я уже сделал несколько неуклюжих попыток векторизации с помощью функции base::replicate, но это не ускоряет процесс. На самом деле это немного медленнее: 

# my for loop approach
my.for.loop <- function() {
  result <- matrix(NA, nr.runs, nr.pairs)
  for (i in 1:nr.runs){
    mix <- sample(nr.pops, replace=FALSE)
    dist.mat.mix <- dist.mat[mix ,mix]
    result[i, ] <- dist.mat.mix[upper.tri(dist.mat.mix, diag=FALSE)]
  }
}

# my replicate approach
my.replicate <- function() {
  results <- t(replicate(nr.runs, {
    mix <- sample(nr.pops, replace=FALSE)
    dist.mat.mix <- dist.mat[mix, mix]
    dist.mat.mix[upper.tri(dist.mat.mix, diag=FALSE)]
  }))
}

# compare speed
require(microbenchmark)
microbenchmark(my.for.loop(), my.replicate(), times=100L)

# Unit: milliseconds
# expr           min     lq      mean    median  uq      max       neval
# my.for.loop()  23.1792 24.4759 27.1274 25.5134 29.0666 61.5616   100
# my.replicate() 25.5293 27.4649 30.3495 30.2533 31.4267 68.6930   100

Я был бы очень признателен за вашу поддержку, если вы знаете, как ускорить мой цикл for, используя аккуратное векторизованное решение. Это вообще возможно?

1 Ответ

0 голосов
/ 28 октября 2019

Немного быстрее: 

minem <- function() {
  result <- matrix(NA, nr.runs, nr.pairs)
  ut <- upper.tri(matrix(NA, 4, 4)) # create upper triangular index matrix outside loop
  for (i in 1:nr.runs) {
    mix <- sample.int(nr.pops) # slightly faster sampling function
    result[i, ] <- dist.mat[mix, mix][ut]
  }
  result
}
microbenchmark(my.for.loop(), my.replicate(), minem(), times = 100L)
# Unit: microseconds
# expr               min      lq      mean   median       uq      max neval cld
# my.for.loop()   75.062  78.222  96.25288  80.1975 104.6915  249.284   100   a
# my.replicate() 118.519 122.667 152.25681 126.0250 165.1355  495.407   100   a
# minem()         45.432  48.000 104.23702  49.5800  52.9380 4848.986   100   a

Обновление: мы можем получить необходимые матричные индексы немного по-другому, поэтому мы можем поднастроить элементы сразу: 

minem4 <- function() {
  n <- dim(dist.mat)[1]
  ut <- upper.tri(matrix(NA, n, n))
  im <- matrix(1:n, n, n)
  p1 <- im[ut]
  p2 <- t(im)[ut]
  dm <- unlist(dist.mat)

  si <- replicate(nr.runs, sample.int(nr.pops))
  p <- (si[p1, ] - 1L) * n + si[p2, ]
  result2 <- matrix(dm[p], nr.runs, nr.pairs, byrow = T)
  result2
}

microbenchmark(my.for.loop(), minem(), minem4(), times = 100L)
# Unit: milliseconds
# expr                min        lq     mean    median        uq       max neval cld
# my.for.loop() 13.797526 14.977970 19.14794 17.071401 23.161867  29.98952   100   b
# minem()        8.366614  9.080490 11.82558  9.701725 15.748537  24.44325   100  a 
# minem4()       7.716343  8.169477 11.91422  8.723947  9.997626 208.90895   100  a

Обновление 2: некоторыедополнительное ускорение мы можем получить, используя dqrng функцию выборки: 

minem5 <- function() {
  n <- dim(dist.mat)[1]
  ut <- upper.tri(matrix(NA, n, n))
  im <- matrix(1:n, n, n)
  p1 <- im[ut]
  p2 <- t(im)[ut]
  dm <- unlist(dist.mat)

  require(dqrng)
  si <- replicate(nr.runs, dqsample.int(nr.pops))
  p <- (si[p1, ] - 1L) * n + si[p2, ]
  result2 <- matrix(dm[p], nr.runs, nr.pairs, byrow = T)
  result2
}

microbenchmark(my.for.loop(), minem(), minem4(), minem5(), times = 100L)
# Unit: milliseconds
# expr                min        lq      mean    median        uq      max neval  cld
# my.for.loop() 13.648983 14.672587 17.713467 15.265771 16.967894 36.18290   100    d
# minem()        8.282466  8.773725 10.679960  9.279602 10.335206 27.03683   100   c 
# minem4()       7.719503  8.208984  9.039870  8.493231  9.097873 25.32463   100  b  
# minem5()       6.134911  6.379850  7.226348  6.733035  7.195849 19.02458   100 a
...