[R] avoid loop with three-dimensional array

[Charles C. Berry]

On Mon, 21 Jul 2008, Michela Cameletti wrote:

> Dear R user,
> I'm trying to find a solution for optimizing my code. I have to run a 50.000
> iteration long simulation and it is absolutely necessary to have an
> optimized code.
>

What you have is just a quadratic form, so t( y ) %*% A %*% y will do it.

You need to rearrange your subscripts, this is one approach:

> X2 <- matrix( aperm(X, c(1,3,2) ),nr=d ) # X is as you defined it.
> nother.result <-  crossprod( matrix(A%*%X2,nc=k), matrix( X2, nc=k) )

Using dim( X2 ) <- ... rather than matrix() will probably speed this 
further.

HTH,

Chuck

p.s. This is the road to perdition. You can optimize the heck out of some 
piece of linear algebra in R only to find that it needs to be written in 
C, but the R code that you have written is so hard to read that you have 
to start from scratch. IIRC, writing block-trace algorithms for mixed 
models takes one down this path.


> I have to do this operation
> *sum_t ( t(X_t) %*% A %*% X_t )*
> where X_t is a (d*k) matrix which changes in time and A is a constant in
> time (d*d) matrix.
> I have put all my X_t in a three dimensional array X of dimension (d,k,T).
>
> At the moment for computing the sum over time I'm doing a for loop and
> saving the resulting (k*k) matrix in a list and at the end I sum the T
> matrices in this list. I'm wondering if there is a better way to do this.
>
> Here an example of what I have to do:
>
> *d=3
> k=2
> T=4
>
> X = array(rnorm(d*k*T),dim=c(d,k,T))
> A = matrix(rnorm(d*d),d,d)
>
> e1  = list()
> for (t in 1:T){ #I would like to avoid this
>      e1[[t]] = t(X[,,t])%*%A%*%X[,,t]
> }
>
> ##############################
> #Function for doing the sum of matrices in a list
> ##############################
>
> sumMatrices <- function(matrices){
>         if (length(matrices) > 2) matrices[[1]] + Recall(matrices[-1])
>         else matrices[[1]] + matrices[[2]]
> }
> ##############################
>
> result = sumMatrices(e1)
> *
>
>
> Thank you in advance for all your help,
> best
> Michela
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> 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.
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901