Skip to content

avoiding loops

4 messages · Ingmar Visser, Bill Venables

#
Hi,
I need to compute an array from a matrix and an array:

A <- array(1:20,c(2,2,5))
B <- matrix(1:10,5)

And I would like the result to be an array consisting of the following:

rbind(A[1,,1]*B[1,],
A[2,,1]*B[1,])

rbind(A[1,,2]*B[2,],
A[2,,2]*B[2,])

rbind(A[1,,3]*B[2,],
A[2,,3]*B[2,])

etc.

Hence the result should have the same dimension as A, ie a series of  
2 by 2 matrices.

Short of a for loop over the the last index of A I have struggled  
with versions of apply but
to no avail ...

Any insights much appreciated,

Best, Ingmar
#
If you have lots of memory there is an obvious strategy:

d12 <- prod(dim(A)[1:2])
A <- A * array(rep(B, each = d12), dim = dim(A))

I don't really see much wrong with the obvious for() loop, though:

for(b in 1:length(B)) A[,,b] <- A[,,b] * B[b]


Bill Venables
CSIRO Laboratories
PO Box 120, Cleveland, 4163
AUSTRALIA
Office Phone (email preferred): +61 7 3826 7251
Fax (if absolutely necessary):  +61 7 3826 7304
Mobile:                         +61 4 8819 4402
Home Phone:                     +61 7 3286 7700
mailto:Bill.Venables at csiro.au
http://www.cmis.csiro.au/bill.venables/ 

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of Ingmar Visser
Sent: Thursday, 27 March 2008 7:58 AM
To: R-help at r-project.org
Subject: [R] avoiding loops

Hi,
I need to compute an array from a matrix and an array:

A <- array(1:20,c(2,2,5))
B <- matrix(1:10,5)

And I would like the result to be an array consisting of the following:

rbind(A[1,,1]*B[1,],
A[2,,1]*B[1,])

rbind(A[1,,2]*B[2,],
A[2,,2]*B[2,])

rbind(A[1,,3]*B[2,],
A[2,,3]*B[2,])

etc.

Hence the result should have the same dimension as A, ie a series of  
2 by 2 matrices.

Short of a for loop over the the last index of A I have struggled  
with versions of apply but
to no avail ...

Any insights much appreciated,

Best, Ingmar

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
#
Thanks for this.
I was afraid someone was going to say this ...
Does this mean the only way of getting this to run faster is by  
moving to C code?
The cases I'm thinking of applying this in have dimensions of A that  
are much larger than
the example, eg n by n by T where n has a max of 10 or so but T could  
be hundreds
or even thousands.
Best, Ingmar

On Mar 27, 2008, at 12:11 AM, <Bill.Venables at csiro.au>
<Bill.Venables at csiro.au> wrote:

            
#
Ingmar Visser asks:
Perhaps, but that's not all that difficult for this kind of operation,
surely?  Even portability across platforms should be pretty
manageable.
What I would try first, then is a loop which makes maximum use of
vectorisation. In this case it would be two loops, of course:

d <- dim(A)
for(i in 1:d[1])
  for(j in 1:d[2])
    A[i,j,] <- A[i,j,] * B

But the following non-loop solution might be your best option:

A <- aperm(aperm(A, c(3,1,2)) * as.vector(B), c(2,3,1))

If you can get your head around working with an A array where what is
now the third dimension becomes the first, then you can eliminate the
aperm()s and it becomes about as fast as C code (which it is, in
effect).  Permuting large arrays can be heavy on memory usage.

Bill Venables.