Vectorizing a complex nested for loop in R (running models on different subsets of a data set, subsetting the data differently for each loop) This is my first time posting so hopefully I manage to get my question across in a clear and concise manner. I'm not sure this is the right place but someone on StackOverflow suggested I ask here.

I am running a cross-validation experiment which involves training a model on a subset of a data set and then validating it against the rest of the data set. This is done by randomly subsetting the data into "folds" and then running the model for each fold. This is then replicated in order to produce standard errors. This is the code I have written in R:

reps = 10
folds = 10
rep.Res&lt;-list()

for(i in 1:reps){
  #i=1
  phen2&lt;-phen
  phen2$fold&lt;-sample(1:folds,nrow(phen),replace=T)

  fold.Res&lt;-list()

  for(j in 1:folds){
    #j=1

    print(paste(i," ",j))

    phen3&lt;-phen2
    phen3[which(phen3$fold==j),2]&lt;-NA

    ETA&lt;-list(list(X = geno_centred, model = 'BayesA'))

    model &lt;- BGLR(y = phen3[,2], ETA = ETA, nIter = 20000, 
                  burnIn = 1500,  thin=10) 

    model.output&lt;-data.frame(ID = phen3[which(phen2$fold == j), 1],
                             PBV = as.numeric(model$yHat[which(phen2$fold == j)]),
                             phen = phen2[which(phen2$fold == j), 2], rep = i, fold = j)


fold.Res[[j]]&lt;-model.output


  }

  rep.Res[[i]]&lt;-do.call(rbind,fold.Res)


}


"phen" is a data frame with a a column of IDs and a column of values. For each replicate, each sample is assigned to a fold randomly. Then, for each fold individuals in the fold are given an NA value, i.e. excluded from the model. 

The model is run using the Bayesian Generalised Linear Regression (BGLR) package. Unfortunately, I can't really give a reproducible example. X is a scaled centred genotype matrix of -1's, 0's and 1's (not critical for my problem here). The output of the model is then written out for each fold, and then finally for each replicate.

My question is, is it possible to vectorize these loops? Using 10 folds and 10 reps means that the model runs 100 time, at 20 000 MCMC iterations each time. This causes the system to slow down significantly before it's even half done. I know that there are many possible approaches using basic R or packages in the tidyverse suite, but I am really a bit lost with such a complex loop. I can imagine that I would try to write at least one of the loops into a function enclosing the other loop but I'm not even sure how that would work...

Again, first time posting so please let me know if I need to add anything to my question! Thanks!

Edit: Here are a small, random centred genotype matrix (X in the call to BGLR) and an accompanying phen file.

X (rows are sample IDs matching those in phen, columns are genomic locations (SNPs)): 

    A  B  C  D  E  F  G  H  I  J  K  L  M  N  O  P  Q  R  S  T
1   1  1  1 -1 -1  0 -1 -1  0  1 -1  0  0  0  0 -1  0  0  1  1
2  -1  1  0 -1  0  0  1 -1  1  0 -1  1 -1  0  0  0 -1 -1 -1 -1
3  -1 -1  1 -1  0  0  1  0  0 -1 -1  0 -1  0  0 -1  1 -1  0 -1
4   0  0 -1  1 -1  0 -1  1  1  1 -1  0  0  0  1  1 -1 -1  0  1
5   1  0  1  1  0  0  1  1  0  0 -1  1  1  1  1 -1 -1  0 -1  0
6   1 -1 -1  0  1 -1  1  1  1 -1  0  0  0 -1  0  1  0 -1 -1  0
7   1 -1 -1  1  1  0 -1  1 -1  1  0 -1  1  1 -1 -1 -1  1 -1  1
8   0  1  0  1 -1 -1  0 -1 -1 -1  1 -1  0  0  0  0  0 -1  0 -1
9  -1  1 -1 -1  1  0  0  0  0  0  1  0 -1 -1 -1  0  1 -1  0 -1
10  0  0 -1  1  0  0 -1  0 -1  1  1  1  0  0  1  1  0 -1  0  0


Phen:

ID  phenotype
1   1000
2   1500
3   1200
4   500
5   700
6   2000
7   1500
8   1000
9   1300
10  900
I'm not familiar with BGLR but I think I can give some pointers.

First to your direct question - yes you can vectorise.  A simple approach is to use one of the apply functions.  lapply usually is a good fit for such loops.  However this is not a magic bullet and I would strongly suggest to review the code.

I think the main thing that slows you down - besides the time BGLR needs - is how R deals with memory.  The relevant properties are:


most objects need to be placed continuous in memory
memory isn't immediately freed - but R waits till it is necessary and does garbage collect then
the time when R creates a copy of memory is a bit tricky.  It is not when one writes an assignment - at this point the variable just holds a reference to the original data.  But as soon as a change on the data is performed it creates a full copy.


All this effects don't appear in a "clean" R session, but as memory clogs things get significantly slower - agreeing with your observation. 

I see following issues in the code:

Growing lists
fold.Res and rep.Res are growing lists.  As grow them they need to be written anew.  This means first all the previous entries need to be copied (which can take some time) - but more important they can not be written in the same place as the new list takes more space.  This will create fragmented memory - meaning whenever R can't fit the new list where a whole from a previous was created it will first need to garbage collect and then maybe defragment memory.  I think this is the serious slow down you are experiencing.
Allocating lists unfortunately doesn't help, but with the apply functions like lapply this is also taken care of.

creating modified copies
You create phen2 and phen3 as copies - but you only modify and use column 2.  Meaning there is some unnecessary copying going on.  Also this will lead faster to some extra garbage collects as memory fills up.

possible side effects of BGLR
While the first runs are relatively free where to put their data, later may run into troubles finding suitable spots - causing again some superfluous memory cleanups.  This is not so easy to deal with.  I worked with stan which is quite the hog - and the best solution was to run it in separate R processes and save the results and the kill it - it was way faster than having these cleanups going on.  The main scripts would then just collect the files.  I described how I do it in my answer to: Using doParallel to cycle through *.rds files (specifying how to do it outside doParallel).  

to style 
It is a bit confusing that you sometimes refer to columns by name and sometimes by column number.  This makes things harder to read than necessary.  I would suggest always using the name approach.

Those apply functions are awkward to write.  Personally I would consider to rewrite it to use dplyr - in particular piping and do.  However that may be to much work.