> -----Original Message----- > From: r-help-boun...@r-project.org > [mailto:r-help-boun...@r-project.org] On Behalf Of Joshua Wiley > Sent: Friday, October 15, 2010 12:23 AM > To: Gregor > Cc: r-help@r-project.org > Subject: Re: [R] fast rowCumsums wanted for calculating the cdf > > Hi, > > You might look at Reduce(). It seems faster. I converted the matrix > to a list in an incredibly sloppy way (which you should not emulate) > because I cannot think of the simple way. > > > probs <- t(matrix(rep(1:10000000), nrow=10)) # matrix with > row-wise probabilites > > F <- matrix(0, nrow=nrow(probs), ncol=ncol(probs)); > > F[,1] <- probs[,1,drop=TRUE]; > > add <- function(x) {Reduce(`+`, x, accumulate = TRUE)} > > > > > > system.time(F.slow <- t(apply(probs, 1, cumsum))) > user system elapsed > 36.758 0.416 42.464 > > > > system.time(for (cc in 2:ncol(F)) { > + F[,cc] <- F[,cc-1,drop=TRUE] + probs[,cc,drop=TRUE]; > + }) > user system elapsed > 0.980 0.196 1.328 > > > > system.time(add(list(probs[,1], probs[,2], probs[,3], > probs[,4], probs[,5], probs[,6], probs[,7], probs[,8], > probs[,9], probs[,10]))) > user system elapsed > 0.420 0.072 0.539
One way to avoid the typing of list(probs[,1], probs[,2], ...) is to use split(probs, col(probs)). However, split() is surprisingly slow in this case. Another approach is to use a matrix multiply, by a ncol(probs) by ncol(probs) matrix with 0's in lower triangle and 1's elsewhere. Here are 4 approaches I've seen mentioned, made into functions that output the matrix requested: f1 <- function (x) t(apply(x, 1, cumsum)) f2 <- function (x){ if (ncol(x) > 1) for (i in 2:ncol(x)) x[, i] <- x[, i] + x[, i - 1] x } f3 <- function (x) matrix(unlist(Reduce(`+`, split(x, col(x)), accumulate = TRUE)), ncol = ncol(x)) f4 <- function (x) x %*% outer(seq_len(ncol(x)), seq_len(ncol(x)), "<=") I tested it as follows > probs <- matrix(runif(1e7), ncol=10, nrow=1e6) > system.time(v1 <- f1(probs)) user system elapsed 19.45 0.25 16.78 > system.time(v2 <- f2(probs)) user system elapsed 0.68 0.03 0.79 > system.time(v3 <- f3(probs)) user system elapsed 38.25 0.24 34.47 > system.time(v4 <- f4(probs)) user system elapsed 0.70 0.05 0.56 > identical(v1,v2) && identical(v1,v3) && identical(v1,v4) [1] TRUE If you have a fancy optimized/multithreaded BLAS linked to your version of R there you may find that %*% is much faster (I'm using the precompiled R for Windows). As ncol(x) gets bigger the %*% approach probably loses its advantage. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com > > > > .539 seconds for 10 vectors each with 1 million elements does not seem > bad. For 30000 each, on my system it took .014 seconds, which for > 200,000 iterations, I estimated as: > > > (200000*.014)/60/60 > [1] 0.7777778 > > (and this is on a stone age system with a single core processor and > only 756MB of RAM) > > It should not be difficult to get the list back to a matrix. > > Cheers, > > Josh > > > > On Thu, Oct 14, 2010 at 11:51 PM, Gregor <mailingl...@gmx.at> wrote: > > Dear all, > > > > Maybe the "easiest" solution: Is there anything that speaks > against generalizing > > cumsum from base to cope with matrices (as is done in matlab)? E.g.: > > > > "cumsum(Matrix, 1)" > > equivalent to > > "apply(Matrix, 1, cumsum)" > > > > The main advantage could be optimized code if the Matrix is > extreme nonsquare > > (e.g. 100,000x10), but the summation is done over the short > side (in this case 10). > > apply would practically yield a loop over 100,000 elements, > and vectorization w.r.t. > > the long side (loop over 10 elements) provides considerable > efficiency gains. > > > > Many regards, > > Gregor > > > > > > > > > > On Tue, 12 Oct 2010 10:24:53 +0200 > > Gregor <mailingl...@gmx.at> wrote: > > > >> Dear all, > >> > >> I am struggling with a (currently) cost-intensive problem: > calculating the > >> (non-normalized) cumulative distribution function, given > the (non-normalized) > >> probabilities. something like: > >> > >> probs <- t(matrix(rep(1:100),nrow=10)) # matrix with > row-wise probabilites > >> F <- t(apply(probs, 1, cumsum)) #SLOOOW! > >> > >> One (already faster, but for sure not ideal) solution - > thanks to Henrik Bengtsson: > >> > >> F <- matrix(0, nrow=nrow(probs), ncol=ncol(probs)); > >> F[,1] <- probs[,1,drop=TRUE]; > >> for (cc in 2:ncol(F)) { > >> F[,cc] <- F[,cc-1,drop=TRUE] + probs[,cc,drop=TRUE]; > >> } > >> > >> In my case, probs is a (30,000 x 10) matrix, and i need to > iterate this step around > >> 200,000 times, so speed is crucial. I currently can make > sure to have no NAs, but > >> in order to extend matrixStats, this could be a nontrivial issue. > >> > >> Any ideas for speeding up this - probably routine - task? > >> > >> Thanks in advance, > >> Gregor > >> > >> ______________________________________________ > >> R-help@r-project.org 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. > > > > ______________________________________________ > > R-help@r-project.org 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. > > > > > > -- > Joshua Wiley > Ph.D. Student, Health Psychology > University of California, Los Angeles > http://www.joshuawiley.com/ > > ______________________________________________ > R-help@r-project.org 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. > ______________________________________________ R-help@r-project.org 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.