Ускорение функций с участием (ы) применяются

Question

Я профилировал свой код, используя lineprof пакет и определили узкие места, чтобы быть в трех функциях perm.stat.list, G.hat, а также emp.FDR, Общей темой, кажется, является использование (s)applyна основе выходных данных профилировщика.

Ниже приведена упрощенная версия моих функций, а также код для генерации воспроизводимого примера, включающего три функции. Я добавил комментарии, чтобы лучше объяснить, что делает каждая функция, и необходимые входные данные.

Я хотел бы значительно ускорить мой код, потому что даже с B=10процесс занимает почти полчаса вычислений. Ввод занимает большую матрицу (10000 x 10000), поэтому скорость важна. В идеале я бы хотел бежать B=5000 перестановки, что также увеличивает время вычислений.

Любые советы по улучшению моего кода с благодарностью.

### Functions ###
    perm.stat.list <- function(samp.dat,N1,N2,B){
      perm.list = NULL
      for (b in 1:B){
        #Permute the row "labels", preserving information across columns
        perm.dat.tmp = samp.dat[sample(nrow(samp.dat)),]

        #Compute the permutation-based test statistics
        #Need to save each (1 x M) permutation vector into a list
        perm.list[[b]] = apply(perm.dat.tmp,2,function(y) t.test(y[1:N1],y[(N1+1):(N1+N2)])$statistic)
      }
      return(perm.list)
    }

    G.hat = function(perm.mat,t){
      #Number of permutations
      B = nrow(perm.mat)
      #Compute an empirical distribution along each COLUMN of permutation matrix
      out = apply(perm.mat,2,function(x) sum(x>t,na.rm = TRUE))/B
      return(out)
    }

    emp.FDR <- function(t.vec,mat){
      #For each value in t.vec, apply G.hat function
      out = sapply(t.vec,function(i) sum(G.hat(mat,i),na.rm = TRUE)/max(sum(t.vec > i,na.rm = TRUE),1))
      return(out)
    }

,

### Generate reproducible example ###

### Global variables ###
#Sample sizes (rows)
N1=3000
N2=7000
#Number of columns
M = 10000
#Number of permutations
B = 10

### Data ###
set.seed(1)
X1 = matrix(rnorm(N1*M),ncol=M)
X2 = matrix(rnorm(N2*M),ncol=M)

### Combine data in one large matrix of size (N1+N2) rows x M columns ###
samp.dat = rbind(X1,X2)

### Compute statistic for each column of samp.dat ###
t.stats = apply(samp.dat,2,
               function(x) t.test(x[1:N1],x[(N1+1):(N1+N2)])$statistic)

### Sort t.stats in decreasing order (not necessarily needed for computation) ###
t.vec = sort(t.stats,decreasing=TRUE)

### Permutation matrix based on the data ###
perm.mat = perm.stat.list(samp.dat=samp.dat,N1=N1,N2=N2,B=B)

eFDR = emp.FDR(t.vec=t.vec,mat=perm.mat)