[R] speeding up aperm/ adding repmat

[Kevin Murphy]

Hi,

I have noticed that aperm is very slow, and I wondered if there was a
way of speeding it up.

Let me tell you a bit about the context of my problem, because perhaps
I shouldn't be using aperm at all.
The context is probabilistic inference in
graphical models. One of the most fundamental operations is two
compute an element-wise multiplication of two arrays of different
sizes, say A and B. This is illustrated below.


for (i in 1:Asz[1]) {
  for (j in 1:Asz[2]) {
    for (k in 1:Asz[3]) {
        C[i,j,k] <- A[i,j,k] * B[i,k]
    }
  }
}

Here, the "domain" of A is [X1 X2 X3] and the domain of B is [X1 X3],
which is why we write A[i,j,k,l] * B[i,k]. (Asz = dim(A).)

In general, we can implement this by making B the same size as A by
replicating its entries, then permuting the dimensions so they match
A's, and then doing element-wise multiplication.
Here is the R code:

B2 <- array(B, Asz) # duplicate elements of B to make it same size as A
diffdom <- Adom[-Bdom] # Adom \ Bdom
B2dom <- c(Bdom, diffdom)
perm <- match(Adom, B2dom)
B3 <- aperm(B2, perm) # make B2 have the same order as A (VERY SLOW)
C <- A * B3

In Matlab, there is a quicker way, using 'reshape' and 'repmat': one
can reshape B to have the same number of dimensions as A, where the new
dimensions have size 1; then one does a repmat to replicate B along
the new dimensions; finally one does the multiplication. This avoids the
cost of a permutation. This is much quicker than the R method,
especially when
the table gets large. Here is the Matlab code:


map = match(Bdom, Adom)
sz = ones(1, length(Adom));
sz(map) = Bsz;
B2 = reshape(B, sz); % add dimensions for the stuff not in A
sz = Asz;
sz(map) = 1; % don't replicate along A's dimensions
B3 = repmat(B2, sz(:)');
                           


So, ANOTHER QUESTION IS: is there something equivalent to Matlab's
multi-dimensional repmat? The recycling rule in R only replicates
along the final dimensions. Here is an example of the difference.
In Matlab, we have

>> B = reshape(1:4, [2 2]) % X1 X3
B =
     1     3
     2     4           

>> B2 = reshape(B, [2 1 2])  % X1 Y2 X3, where Y2 is a dummy dimension
B2(:,:,1) =
     1
     2
B2(:,:,2) =
     3
     4        

>> B3 = repmat(B2, [1 2 1]) % X1 X2 X2
B3(:,:,1) =
     1     1
     2     2
B3(:,:,2) =
     3     3
     4     4           



In R, we have


> B <- array(1:4, c(2,2))
> B
     [,1] [,2]
[1,]    1    3
[2,]    2    4



> B2 <- B; dim(B2) <- c(2,1,2)
> B2
, , 1

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

, , 2

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


> B3 <- array(B2, c(2,2,2))
> B3
, , 1

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

, , 2

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


So matlab "pumps up" the middle, whereas R just appends to the end.


Thanks!
Kevin
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._