From: <Bill.Venables_at_csiro.au>

Date: Fri, 28 Mar 2008 08:58:30 +1000

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. Received on Thu 27 Mar 2008 - 23:02:11 GMT

Date: Fri, 28 Mar 2008 08:58:30 +1000

Ingmar Visser asks:

*>
*

> 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?
*

Perhaps, but that's not all that difficult for this kind of operation, surely? Even portability across platforms should be pretty manageable.

> 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.
*

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.

> Best, Ingmar

*>
**> On Mar 27, 2008, at 12:11 AM, <Bill.Venables_at_csiro.au>
**> <Bill.Venables_at_csiro.au> wrote:
**>
**> > 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.
**> >
**> >
**>
*

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. Received on Thu 27 Mar 2008 - 23:02:11 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Thu 27 Mar 2008 - 23:30:26 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*