Более быстрое вычисление двойного цикла?

Question

У меня есть фрагмент рабочего кода, который занимает слишком много часов (дней?) Для вычисления. У меня есть разреженная матрица из 1 и 0, мне нужно вычесть каждую строку из любой другой строки во всех возможных комбинациях, умножить полученный вектор на другой вектор и, наконец, усреднить значения в нем, чтобы получить единственный скаляр, который мне нуженвставить в матрицу. То, что у меня есть:

m <- matrix( 
c(0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0), nrow=4,ncol=4,
byrow = TRUE)   

b <- c(1,2,3,4)

for (j in 1:dim(m)[1]){
 for (i in 1:dim(m)[1]){
    a <- m[j,] - m[i,]
    a[i] <- 0L
    a[a < 0] <- 0L
    c <- a*b
    d[i,j] <- mean(c[c > 0])
 }
}

Требуемый вывод - это матрица с такими же размерами m, где каждая запись является результатом этих операций. Этот цикл работает, но есть ли идеи, как сделать это более эффективным? Спасибо

Answer 1

1) создать тестовую разреженную матрицу:

nc <- nr <- 100
p <- 0.001
require(Matrix)
M <- Matrix(0L, nr, nc, sparse = T) # 0 matrix
n1 <- ceiling(p * (prod(dim(M)))) # 1 count
M[1:n1] <- 1L # fill only first column, to approximate max non 0 row count
# (each row has at maximum 1 positive element)
sum(M)/(prod(dim(M)))

b <- 1:ncol(M)

sum(rowSums(M))

Итак, если указанная пропорция правильная, то у нас будет не более 10 строк, которые содержат ненулевые элементы

На основании этого факта иВаши предоставленные вычисления:

# a <- m[j, ] - m[i, ]
# a[i] <- 0L
# a[a < 0] <- 0L
# c <- a*b
# mean(c[c > 0])

мы можем видеть, что результат будет значимым только для m[, j] строк, которые имеют по крайней мере 1 элемент, отличный от 0

==> мы можем пропустить вычисления длявсе m[, j], которые содержат только 0, поэтому:

minem <- function() { # write as function
  t1 <- proc.time() # timing
  require(data.table)
  i <- CJ(1:nr, 1:nr) # generate all combinations
  k <- rowSums(M) > 0L # get index where at least 1 element is greater that 0
  i <- i[data.table(V1 = 1:nr, k), on = 'V1'] # merge
  cat('at moust', i[, sum(k)/.N*100], '% of rows needs to be calculated \n')
  i[k == T, rowN := 1:.N] # add row nr for 0 subset
  i2 <- i[k == T] # subset only those indexes who need calculation
  a <- M[i2[[1]],] - M[i2[[2]],] # operate on all combinations at once
  a <- drop0(a) # clean up 0

  ids <- as.matrix(i2[, .(rowN, V2)]) # ids for 0 subset
  a[ids] <- 0L # your line: a[i] <- 0L
  a <- drop0(a) # clean up 0

  a[a < 0] <- 0L # the same as your line
  a <- drop0(a) # clean up 0

  c <- t(t(a)*b) # multiply each row with vector
  c <- drop0(c) # clean up 0

  c[c < 0L] <- 0L # for mean calculation
  c <- drop0(c) # clean up 0

  r <- rowSums(c)/rowSums(c > 0L) # row means
  i[k == T, result := r] # assign results to data.table
  i[is.na(result), result := NaN] # set rest to NaN
  d2 <- matrix(i$result, nr, nr, byrow = F) # create resulting matrix
  t2 <- proc.time() # timing
  cat(t2[3] - t1[3], 'sec \n')
  d2
}
d2 <- minem()
# at most 10 % of rows needs to be calculated 
# 0.05 sec 

Тест на меньшем примере, если результаты соответствуют

d <- matrix(NA, nrow(M), ncol(M))
for (j in 1:dim(M)[1]) {
  for (i in 1:dim(M)[1]) {
    a <- M[j, ] - M[i, ]
    a[i] <- 0L
    a[a < 0] <- 0L
    c <- a*b
    d[i, j] <- mean(c[c > 0])
  }
}
all.equal(d, d2)

Можем ли мы получить результаты для вашего реального размера данных?:

# generate data:
nc <- nr <- 6663L
b <- 1:nr
p <- 0.0001074096 # proportion of 1s
M <- Matrix(0L, nr, nc, sparse = T) # 0 matrix
n1 <- ceiling(p * (prod(dim(M)))) # 1 count
M[1:n1] <- 1L

object.size(as.matrix(M))/object.size(M)
# storing this data in usual matrix uses 4000+ times more memory

# calculation:
d2 <- minem()
# at most 71.57437 % of rows needs to be calculated 
# 28.33 sec 

Так что вам нужно преобразовать матрицу в разреженную с помощью

M <- Matrix(m, sparse = T)

Answer 2

Мое глупое решение - использовать функцию apply или sapply вместо цикла for для выполнения итераций:

sapply(1:dim(m)[1], function(k) {z <- t(apply(m, 1, function(x) m[k,]-x)); diag(z) <- 0; z[z<0] <- 0; apply(t(apply(z, 1, function(x) x*b)),1,function(x) mean(x[x>0]))})

Я пытался сравнить ваше решение и это с точки зрения выполнениявремя в моем компьютере, ваше занимает

t1 <- Sys.time()
d1 <- m
for (j in 1:dim(m)[1]){
  for (i in 1:dim(m)[1]){
    a <- m[j,] - m[i,]
    a[i] <- 0L
    a[a < 0] <- 0L
    c <- a*b
    d1[i,j] <- mean(c[c > 0])
  }
}
Sys.time()-t1

Ваши потребности Time difference of 0.02799988 secs. Для меня это немного, но не слишком сильно, например, Time difference of 0.01899815 secs, когда вы запускаете

t2 <- Sys.time()
d2 <- sapply(1:dim(m)[1], function(k) {z <- t(apply(m, 1, function(x) m[k,]-x)); diag(z) <- 0; z[z<0] <- 0; apply(t(apply(z, 1, function(x) x*b)),1,function(x) mean(x[x>0]))})
Sys.time()-t2

Вы можете попробовать его на своем компьютере с более крупной матрицей, удачи!