[R] matrix row product and cumulative product

[Jeff Laake]

I spent a lot of time searching and came up empty handed on the 
following query. Is there an equivalent to rowSums that does product or 
cumulative product and avoids use of apply or looping? I found a rowProd 
in a package but it was a convenience function for apply. As part of a 
likelihood calculation called from optim, I’m computing products and 
cumulative products of rows of matrices with far more rows than columns. 
I started with apply and after some thought realized that a loop of 
columns might be faster and it was substantially faster (see below). 
Because the likelihood function is called many times I’d like to speed 
it up even more if possible.

Below is an example showing the cumprod.matrix and prod.matrix looping 
functions that I wrote and some timing comparisons to the use of apply 
for different column and row dimensions. At this point I’m better off 
with looping but I’d like to hear of any further suggestions.

Thanks –jeff

 > prod.matrix=function(x)
+ {
+ y=x[,1]
+ for(i in 2:dim(x)[2])
+ y=y*x[,i]
+ return(y)
+ }

 > cumprod.matrix=function(x)
+ {
+ y=matrix(1,nrow=dim(x)[1],ncol=dim(x)[2])
+ y[,1]=x[,1]
+ for (i in 2:dim(x)[2])
+ y[,i]=y[,i-1]*x[,i]
+ return(y)
+ }

 > N=10000000
 > xmat=matrix(runif(N),ncol=10)
 > system.time(cumprod.matrix(xmat))
user system elapsed
1.07 0.09 1.15
 > system.time(t(apply(xmat,1,cumprod)))
user system elapsed
29.27 0.21 29.50
 > system.time(prod.matrix(xmat))
user system elapsed
0.29 0.00 0.30
 > system.time(apply(xmat,1,prod))
user system elapsed
30.69 0.00 30.72
 > xmat=matrix(runif(N),ncol=100)
 > system.time(cumprod.matrix(xmat))
user system elapsed
1.05 0.13 1.18
 > system.time(t(apply(xmat,1,cumprod)))
user system elapsed
3.55 0.14 3.70
 > system.time(prod.matrix(xmat))
user system elapsed
0.38 0.01 0.39
 > system.time(apply(xmat,1,prod))
user system elapsed
2.87 0.00 2.89
 > xmat=matrix(runif(N),ncol=1000)
 > system.time(cumprod.matrix(xmat))
user system elapsed
1.30 0.18 1.46
 > system.time(t(apply(xmat,1,cumprod)))
user system elapsed
1.77 0.27 2.05
 > system.time(prod.matrix(xmat))
user system elapsed
0.46 0.00 0.47
 > system.time(apply(xmat,1,prod))
user system elapsed
0.7 0.0 0.7
 > xmat=matrix(runif(N),ncol=10000)
 > system.time(cumprod.matrix(xmat))
user system elapsed
1.28 0.00 1.29
 > system.time(t(apply(xmat,1,cumprod)))
user system elapsed
1.19 0.08 1.26
 > system.time(prod.matrix(xmat))
user system elapsed
0.40 0.00 0.41
 > system.time(apply(xmat,1,prod))
user system elapsed
0.57 0.00 0.56
 > xmat=matrix(runif(N),ncol=100000)
 > system.time(cumprod.matrix(xmat))
user system elapsed
3.18 0.00 3.19
 > system.time(t(apply(xmat,1,cumprod)))
user system elapsed
1.42 0.21 1.64
 > system.time(prod.matrix(xmat))
user system elapsed
1.05 0.00 1.05
 > system.time(apply(xmat,1,prod))
user system elapsed
0.82 0.00 0.81
 > R.Version()
$platform
[1] "i386-pc-mingw32"
.
.
.
$version.string
[1] "R version 2.7.1 (2008-06-23)"