Как ускорить цикл for в R для сопоставления вложенных матриц и colSums

Question

У меня, по-видимому, простая проблема, для которой мне требуется более быстрая реализация R, чем та, которую я разработал

Я инициализирую случайное начальное число и измерения для этого примера:

set.seed(1)
d1<-400
d2<-20000
d3<-50

У меня естьматрица X с размерами d1 x d2 :

X<-as.data.frame(matrix(rnorm(d1*d2),nrow=d1,ncol=d2))
rownames(X)<-paste0("row",1:nrow(X))
colnames(X)<-paste0("col",1:ncol(X))

и вектор u с d1 индексами строк:

u<-sample(rownames(X),nrow(X),replace=TRUE)

У меня также есть матрица C с именованными строками и размерами d3 x d2 :

C<-matrix(rnorm(d3*d2),nrow=d3,ncol=d2)
rownames(C)<-sample(rownames(X),nrow(C),replace=FALSE)

Теперь со следующим очень медленно цикл Я заполняю матрицу C суммами соответствующих строк X:

system.time(
    for(i in 1:nrow(C)){
        indexes<-which(u==rownames(C)[i])
        C[i,] <- colSums(X[indexes,])
    }
)

Эта операция занимает примерно 11,5 секунд на моем ПК, но я уверен, что ее можно ускорить, избегаяпетля.Есть идеи?Большое спасибо!

Answer 1

Вы можете попробовать использовать sapply для цикла.

system.time(
  C2 <- `dimnames<-`(t(sapply(match(rownames(C), u), function(x) 
    colSums(X[x, ]))), list(rownames(C), NULL))
)
#  user  system elapsed 
# 20.06    0.03   20.14 

stopifnot(all.equal(C, C2))

По сравнению с

system.time(
  for(i in 1:nrow(C)){
    indexes <- which(u == rownames(C)[i])
    C[i, ] <- colSums(X[indexes, ])
  }
)
#  user  system elapsed 
# 20.76    0.69   28.30  

На данный момент, однако, это только одно измерение.

Обновление

Кажется, чтобы работать немного быстрее ...

Unit: seconds
    expr      min       lq     mean   median       uq      max neval cld
 forloop 20.44852 20.57730 21.67771 20.74106 21.01723 29.63220    10   a
  sapply 19.86707 20.17126 21.34529 20.50283 20.81254 29.73764    10   a

Обновление 2

Но вы можете сделать это с parallel::parSapply.

system.time({
  library(parallel)
  cl <- makeCluster(detectCores() - 1)
  clusterExport(cl, c("C", "u", "X"))
  C3 <- parSapply(cl, match(rownames(C), u), function(x) colSums(X[x, ]))
  stopCluster(cl)
  C3 <- `dimnames<-`(t(C3), list(rownames(C), NULL))
})
# user  system elapsed 
# 0.81    3.16    9.82

stopifnot(all.equal(C, C3))

Теперь моя машина работает так же быстро, как ваша, с for -loop:)

Answer 2

Просто используйте matrixStats::colSums2 с возможностью передачи индексов строк и перемещения rownames() за пределы цикла (X необходимо преобразовать в матрицу):

Xm <- as.matrix(X)
names_of_rows <- rownames(C)
system.time(for (i in 1:nrow(C)) {
  indexes <- which(u == names_of_rows[i])
  C[i, ] <-  matrixStats::colSums2(Xm, rows = indexes)
})
# 0.03 sek

Answer 3

Вы можете найти решение data.table здесь. Если OP хочет только решение на основе R, я удалю этот пост:

library(data.table)
mtd_dt <- function() {
    setDT(dtX)[, u := as.integer(gsub("row","",u))]
    mX <- melt(dtX, id.var="u", variable.name="col")
    C2 <- data.table(rn=seq_len(nrow(C)), u=as.integer(gsub("row","",rownames(C))))
    dcast(mX[C2, on=.(u)][, sum(value), by=.(rn, col)], rn ~ col, value.var="V1")[,
        "NA" := NULL][,
            lapply(.SD, function(x) replace(x, is.na(x), 0))]
}

тайминги:

# A tibble: 2 x 14
  expression      min     mean   median      max `itr/sec` mem_alloc  n_gc n_itr total_time result                    memory                time    gc             
  <chr>      <bch:tm> <bch:tm> <bch:tm> <bch:tm>     <dbl> <bch:byt> <dbl> <int>   <bch:tm> <list>                    <list>                <list>  <list>         
1 mtd0()        59.1s    59.1s    59.1s    59.1s    0.0169     447MB    24     1      59.1s <dbl [50 x 20,000]>       <Rprofmem [44,515 x ~ <bch:t~ <tibble [1 x 3~
2 mtd_dt()       2.7s     2.7s     2.7s     2.7s    0.370      309MB     4     1       2.7s <data.table [50 x 20,001~ <Rprofmem [88,029 x ~ <bch:t~ <tibble [1 x 3~

временной код:

mtd0 <- function() {
    for (i in 1:nrow(C)) {
        indexes <- which(u==rownames(C)[i])
        C[i, ] <- colSums(X[indexes, ])
    }
    C
}

bench::mark(mtd0(), mtd_dt(), check=FALSE)

данные:

library(data.table)
set.seed(0)
#d1 <- 10
#d2 <- 10
#d3 <- 5
d1<-400
d2<-20000
d3<-50

X <- as.data.frame(matrix(rnorm(d1*d2),nrow=d1,ncol=d2))
rownames(X) <- paste0("row",1:nrow(X))
colnames(X) <- paste0("col",1:ncol(X))
dtX <- X

u <- sample(rownames(X),nrow(X),replace=TRUE)

C <- matrix(0,nrow=d3,ncol=d2)
rownames(C) <- sample(rownames(X),nrow(C),replace=FALSE)