Gabor Grothendieck
2006-Jul-17 02:29 UTC
[Rd] multiplying multidimensional arrays (was: Re: [R] Manipulation involving arrays)
I am moving this to r-devel. The problem and solution below posted on r-help could have been a bit slicker if %*% worked with multidimensional arrays multiplying them so that if the first arg is a multidimensional array it is mulitplied along the last dimension (and first dimension for the second arg). Then one could have written: Tbar <- tarray %*% t(wt) / rep(wti, each = 9) which is a bit nicer than what had to be done, see below, given that %*% only works with matrices. I suggest that %*% be so extended to multidimensional arrays. Note that this is upwardly compatible and all existing cases would continue to work unchanged. On 7/16/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:> The double loop is the same as: > > Tbar[] <- matrix(tarray, 9) %*% t(wt) / rep(wti, each = 9) > > > On 7/16/06, RAVI VARADHAN <rvaradhan at jhmi.edu> wrote: > > Hi, > > > > I have the following piece of code that is part of a larger function. This piece is the most time consuming part of the function, and I would like to make this a bit more efficient. Could anyone suggest a way to do this faster? > > > > In particular, I would like to replace the nested "for" loop with a faster construct. I tried things like "kronecker" and "outer" combined with apply, but couldn't get it to work. > > > > > > Here is a sample code: > > > > ########################## > > n <- 120 > > sigerr <- 5 > > covmat <- diag(c(8,6,3.5)) > > mu <- c(105,12,10) > > mcsamp <- 10000 > > > > Tbar <- array(0, dim=c(3,3,n)) > > > > # theta is a mcsamp x 3 matrix > > theta <- mvrnorm(mcsamp, mu = mu, Sigma = covmat) > > > > wt <- matrix(runif(n*mcsamp),n,mcsamp) > > wti <- apply(wt,1,sum) > > > > tarray <- array(apply(theta,1,function(x)outer(x,x)),dim=c(3,3,mcsamp)) > > > > for (i in 1:n) { > > for (k in 1:mcsamp) { > > Tbar[,,i] <- Tbar[,,i] + wt[i,k] * tarray[,,k] > > } > > Tbar[,,i] <- Tbar[,,i] / wti[i] > > } > > >
RAVI VARADHAN
2006-Jul-17 13:18 UTC
[Rd] multiplying multidimensional arrays (was: Re: [R] Manipulation involving arrays)
Dear Gabor, Thank you very much for your solution. It speeded calculations by a factor of two. Now to your recommendation about making this solution a basic operation. I completely agree with your suggestion. That is exactly what I would have hoped for. In fact, my first try was to do: Tbar <- tarray %*% t(wt) Since tarray is (3,3,mcsamp) and wt is (n,mcsamp), I figured that the result would be a (3,3,n) array that would sum over mcsamp, as in the case of 2-dim array multiplication. But obviously that didn't work as "non-conformable". Your simple trick of forcing the 3-dim array to a 2-dim matrix using matrix(tarray, 3*3) is quite clever. Thanks again for your help. Best, Ravi. ----- Original Message ----- From: Gabor Grothendieck <ggrothendieck at gmail.com> Date: Sunday, July 16, 2006 10:29 pm Subject: multiplying multidimensional arrays (was: Re: [R] Manipulation involving arrays)> I am moving this to r-devel. > > The problem and solution below posted on r-help could have been > a bit slicker if %*% worked with multidimensional arrays multiplying > them so that if the first arg is a multidimensional array it is > mulitpliedalong the last dimension (and first dimension for the > second arg). > Then one could have written: > > Tbar <- tarray %*% t(wt) / rep(wti, each = 9) > > which is a bit nicer than what had to be done, see below, given > that %*% only > works with matrices. > > I suggest that %*% be so extended to multidimensional arrays. Note > that this is upwardly compatible and all existing cases would continue > to work unchanged. > > On 7/16/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote: > > The double loop is the same as: > > > > Tbar[] <- matrix(tarray, 9) %*% t(wt) / rep(wti, each = 9) > > > > > > On 7/16/06, RAVI VARADHAN <rvaradhan at jhmi.edu> wrote: > > > Hi, > > > > > > I have the following piece of code that is part of a larger > function. This piece is the most time consuming part of the > function, and I would like to make this a bit more efficient. > Could anyone suggest a way to do this faster? > > > > > > In particular, I would like to replace the nested "for" loop > with a faster construct. I tried things like "kronecker" and > "outer" combined with apply, but couldn't get it to work. > > > > > > > > > Here is a sample code: > > > > > > ########################## > > > n <- 120 > > > sigerr <- 5 > > > covmat <- diag(c(8,6,3.5)) > > > mu <- c(105,12,10) > > > mcsamp <- 10000 > > > > > > Tbar <- array(0, dim=c(3,3,n)) > > > > > > # theta is a mcsamp x 3 matrix > > > theta <- mvrnorm(mcsamp, mu = mu, Sigma = covmat) > > > > > > wt <- matrix(runif(n*mcsamp),n,mcsamp) > > > wti <- apply(wt,1,sum) > > > > > > tarray <- > array(apply(theta,1,function(x)outer(x,x)),dim=c(3,3,mcsamp))> > > > > for (i in 1:n) { > > > for (k in 1:mcsamp) { > > > Tbar[,,i] <- Tbar[,,i] + wt[i,k] * tarray[,,k] > > > } > > > Tbar[,,i] <- Tbar[,,i] / wti[i] > > > } > > > > > >