[R] Speeding up resampling of rows from a large matrix

[Bill.Venables at csiro.au]

Here is a possibility.  The only catch is that if a pair of rows is
selected twice you will get the results in a block, not scattered at
random throughout the columns of G.  I can't see that as a problem.

### --- start code excerpt ---
nSNPs <- 1000
H <- matrix(sample(0:1, 120*nSNPs , replace=T), nrow=120)

# G <- matrix(0, nrow=3, ncol=nSNPs)

# Keep in mind that the real H is 120 x 65000

ij <- as.matrix(subset(expand.grid(i = 1:120, j = 1:120), i < j))

nResamples <- 3000
sel <- sample(1:nrow(ij), nResamples, rep = TRUE)
repf <- table(sel)                   # replication factors
ij <- ij[as.numeric(names(repf)), ]  # distinct choice made

G <- matrix(0, nrow = 3, ncol = nrow(ij))  # for now

for(j in 1:ncol(G))
  G[,j] <- rowSums(outer(0:2, colSums(H[ij[j, ], ]), "=="))

G <- G[, rep(1:ncol(G), repf)]         # bulk up the result

# _
# _
# _
# _pair <- replicate(nResamples, sample(1:120, 2))
# _
# _gen <- function(x){g <- sum(x); c(g==0, g==1, g==2)}
# _
# _for (i in 1:nResamples){
# _    G <- G + apply(H[pair[,i],], 2, gen)
# _}
### --- end of code excerpt ---

I did a timing on my machine which is a middle-of-the range windows
monstrosity...

> system.time({
+ 
+ nSNPs <- 1000
+ H <- matrix(sample(0:1, 120*nSNPs , replace=T), nrow=120)
+ 
+ # G <- matrix(0, nrow=3, ncol=nSNPs)
+ 
+ # Keep in mind that the real H is 120 x 65000
+ 
+ ij <- as.matrix(subset(expand.grid(i = 1:120, j = 1:120), i < j))
+ 
+ nResamples <- 3000
+ sel <- sample(1:nrow(ij), nResamples, rep = TRUE)
+ repf <- table(sel)                   # replication factors
+ ij <- ij[as.numeric(names(repf)), ]  # distinct choice made
+ 
+ G <- matrix(0, nrow = 3, ncol = nrow(ij))  # for now
+ 
+ for(j in 1:ncol(G))
+   G[,j] <- rowSums(outer(0:2, colSums(H[ij[j, ], ]), "=="))
+ 
+ G <- G[, rep(1:ncol(G), repf)]         # bulk up the result
+ 
+ # _
+ # _
+ # _
+ # _pair <- replicate(nResamples, sample(1:120, 2))
+ # _
+ # _gen <- function(x){g <- sum(x); c(g==0, g==1, g==2)}
+ # _
+ # _for (i in 1:nResamples){
+ # _    G <- G + apply(H[pair[,i],], 2, gen)
+ # _}
+ #
_#----------------------------------------------------------------------
---------
+ # _
+ })
   user  system elapsed 
   0.97    0.00    0.99 


Less than a second.  Somewhat of an improvement on the 80 minutes, I
reckon.  This will increase, of course when you step the size of the H
matrix up from 1000 to 65000 columns

Bill Venables
CSIRO Laboratories
PO Box 120, Cleveland, 4163
AUSTRALIA
Office Phone (email preferred): +61 7 3826 7251
Fax (if absolutely necessary):  +61 7 3826 7304
Mobile:                (I don't have one!)
Home Phone:                     +61 7 3286 7700
mailto:Bill.Venables at csiro.au
http://www.cmis.csiro.au/bill.venables/ 

-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Juan Pablo
Lewinger
Sent: Friday, 25 May 2007 4:04 PM
To: r-help at stat.math.ethz.ch
Subject: [R] Speeding up resampling of rows from a large matrix

I'm trying to:

Resample with replacement pairs of distinct rows from a 120 x 65,000 
matrix H of 0's and 1's. For each resampled pair sum the resulting 2 
x 65,000 matrix by column:

     0 1 0 1 ...
+
     0 0 1 1 ...
    _______
=  0 1 1 2 ...

For each column accumulate the number of 0's, 1's and 2's over the 
resamples to obtain a 3 x 65,000 matrix G.

For those interested in the background, H is a matrix of haplotypes, 
each pair of haplotypes forms a genotype, and each column corresponds 
to a SNP. I'm using resampling to compute the null distribution of 
the maximum over correlated SNPs of a simple statistic.


The code:
#-----------------------------------------------------------------------
--------
nSNPs <- 1000
H <- matrix(sample(0:1, 120*nSNPs , replace=T), nrow=120)
G <- matrix(0, nrow=3, ncol=nSNPs)
# Keep in mind that the real H is 120 x 65000

nResamples <- 3000
pair <- replicate(nResamples, sample(1:120, 2))

gen <- function(x){g <- sum(x); c(g==0, g==1, g==2)}

for (i in 1:nResamples){
    G <- G + apply(H[pair[,i],], 2, gen)
}
#-----------------------------------------------------------------------
--------
The problem is that the loop takes about 80 mins to complete and I 
need to repeat the whole thing 10,000 times, which would then take 
over a year and a half!

Is there a way to speed this up so that the full 10,000 iterations 
take a reasonable amount of time (say a week)?

My machine has an Intel Xeon 3.40GHz CPU with 1GB of RAM

 > sessionInfo()
R version 2.5.0 (2007-04-23)
i386-pc-mingw32

I would greatly appreciate any help.

Juan Pablo Lewinger
Department of Preventive Medicine
Keck School of Medicine
University of Southern California

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.