[R] speeding up aperm/ adding repmat
Kevin Murphy
murphyk at cs.berkeley.edu
Wed Jul 11 01:26:59 CEST 2001
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list