CODEWITHSUNDEEP
arraysrmatrixlinear-algebra
wolfsatthedoor
wolfsatthedoor

Reputation: 7313

Divide every slice of a matrix in an array by its own vector?

Suppose I have two arrays (or tensors if tensor package is needed)

dim(Xbeta)
  products      draws Households
        13      20        10

dim(denom)
        1       20        10


set.seed(1)
Xbeta=array(rnorm(13*20*10,0,1),dim=c(13,20,10))
denom=array(rnorm(1*20*10,0,1),dim=c(1,20,10))

Without looping, I want to do the following:

   for(i in 1:10){

    Xbeta[,,i]=t(t(Xbeta[,,i]) / denom[,,i])

       }

I want to to divide each column in Xbeta[,,i] slice by each corresponding number in denom[,,i].

For example...Xbeta[,1,i]/denom[,1,i]...etc

Upvotes: 2

Views: 417

Answers (3)

bgoldst
bgoldst

Reputation: 35324

You can avoid looping and replication by (1) 3-dimensionally transposing the numerator array and (2) flattening the denominator array to a vector, such that the division operation will naturally cycle the incomplete denominator vector across the entirety of the transposed numerator array in such a way that the data lines up the way you want. You then must 3-dimensionally "untranspose" the result to get back the original transposition.

aperm(aperm(Xbeta,c(2,3,1))/c(denom),c(3,1,2));

The first call to aperm() transposes columns to rows, z-slices to columns, and rows to z-slices. The c() call on denom flattens the denominator array to a vector, because when cycling, we don't care about dimensionality. The final call to aperm() reverses the original transposition.

To go into more detail about the logic of this solution, what you have with your inputs is basically a vector of divisors per z-slice of the numerator array, and you want to apply each divisor to every row of the corresponding z-slice and column. This means the vector of divisors must be applied across columns, first-and-foremost, and then, as each denominator z-slice is exhausted, applied across numerator z-slices. After a complete row (covering all z-slices in the row) of the numerator array has been exhausted, the entirety of the denominator vector has been exhausted, causing it to be cycled back to the beginning for the next row of the numerator array.

See https://stat.ethz.ch/R-manual/R-devel/library/base/html/aperm.html.

For a rough idea on performance:

r> set.seed(1);
r> Xbeta <- array(rnorm(13*20*10,0,1), dim=c(13,20,10) );
r> denom <- array(rnorm(1*20*10,0,1), dim=c(1,20,10) );
r> robert <- function() { result <- array(NA, dim=c(13,20,10) ); for (i in 1:10) { result[,,i] <- t(t(Xbeta[,,i]) / denom[,,i]); }; };
r> andre <- function() { denom_myVersion <- array(rep(c(denom), each=13 ), c(13,20,10) ); result <- Xbeta / denom_myVersion; };
r> bgoldst <- function() { result <- aperm(aperm(Xbeta,c(2,3,1))/c(denom),c(3,1,2)); };
r> N <- 99999;
r> system.time({ replicate(N, robert() ); });
   user  system elapsed
 25.421   0.000  25.440
r> system.time({ replicate(N, andre() ); });
   user  system elapsed
 12.578   0.594  13.283
r> system.time({ replicate(N, bgoldst() ); });
   user  system elapsed
  8.484   0.594   9.142

Also, as a general recommendation, it is helpful (for both questioners and answerers) to present these kinds of problems using minimal sample input, e.g.:

r> n <- array(1:12,dim=c(2,3,2)); n;
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

r> d <- array(1:6,dim=c(1,3,2)); d;
, , 1

     [,1] [,2] [,3]
[1,]    1    2    3

, , 2

     [,1] [,2] [,3]
[1,]    4    5    6

r> aperm(aperm(n,c(2,3,1))/c(d),c(3,1,2));
, , 1

     [,1] [,2]     [,3]
[1,]    1  1.5 1.666667
[2,]    2  2.0 2.000000

, , 2

     [,1] [,2]     [,3]
[1,] 1.75  1.8 1.833333
[2,] 2.00  2.0 2.000000

Upvotes: 1

Andre Michaud
Andre Michaud

Reputation: 128

I'm not sure why you choose a 3-dimensional array for the denon. Anyway, this can be done by paying close attention to how these numbers are stored in memory. In an array the first dimensions "moves the fastest". By replicating the denom values 13 times "each" then you create an array with the exact same dimensions as your numerator.

So, let's test it out: Let's save the ramdom values so we can use them for both methods:

set.seed(1)
Num_2600 <- rnorm(13*20*10,0,1)
Denom_200 <- rnorm(20*10,0,1)
Xbeta=array(Num_2600,dim=c(13,20,10))
denom=array(Denom_200,dim=c(1,20,10))
Your_result <- array(NA, dim=c(13,20,10))

Your code gives:

for(i in 1:10){
Your_result[,,i] <- t(t(Xbeta[,,i]) / denom[,,i])
   }

My code:

denom_myVersion <- array(rep( Denom_200 , each=13), c(13,20,10))

>  all(Your_result == Xbeta / denom_myVersion)
[1] TRUE
>

So we get the same results. The hard part is how to decide how to replicate so the numbers fall in the right spot. Notice:

denom_myVersion <- array(rep( Denom_200 , times=13), c(13,20,10))

>  all(Your_result == Xbeta / denom_myVersion)
[1] FALSE
>

With 'each' as a parameter in rep each element is repeated 13 times before going to the next element. With times, the whole vector is repeated 13 times. Compare:

> rep(1:3, each =3)
[1] 1 1 1 2 2 2 3 3 3
> rep(1:3, times=3)
[1] 1 2 3 1 2 3 1 2 3

Upvotes: 1

Rick
Rick

Reputation: 898

# Is this what you're looking for?
Xbeta <- array(rnorm(13*20*10,0,1),dim=c(13,20,10))
denom <- array(rnorm(1*20*10,0,1),dim=c(1,20,10))
div.list <- sapply(1:10, FUN = function(x) t(Xbeta[,,x]) / denom[,,x], simplify = FALSE) 
result <-   array(do.call(c, div.list), dim = dim(Xbeta)[c(2,1,3)])

Upvotes: 1

Related Questions