[R] recursive matrix multiplication question

[Patrick Burns]

It would be much easier if you started out with a 3-dimensional
array. A large part of what makes R a good environment is that
it gives you the flexibility to have data structures that match how
you think about the data.

h0arr <- h0m
dim(h0arr) <- c(2, 2, ncol(h0m))

Then you can write a function to do the entire task that you want:

tseqeigen <- function(x, times=1:(nm-1)) {
dx <- dim(x)
if(length(dx) != 3) stop("x must be 3-dimensional array")
if(dx[1] != dx[2]) stop("first two dimensions of x must be equal")
nm <- dx[3]
ans <- array(NA, c(length(times), dx[1]))
for(i in seq(along=times)) {
this.t <- times[i]
this.mm <- x[,, this.t] %*% x[,, this.t + 1]
ans[i,] <- eigen(this.mm)$values
}
ans
}

This could then be used like:

tseqeigen(h0arr)

or

tseqeigen(h0arr, 1:4)

Enhancements to the function might include giving it a better
name, and dealing with dimnames.

There is some material in S Poetry on dealing with higher
dimensional arrays.


Patrick Burns

Burns Statistics
patrick at burns-stat.com
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of S Poetry and "A Guide for the Unwilling S User")

Ivan Kautter wrote:

> List serv subscribers,
>
> I am wondering if there is an efficient means of doing recursive 
> matrix multiplication in R. My data resides in a 4 X 2541 matrix from 
> which I need to extract 2541 2X2 matrices based on each row. If I 
> attempt something like this:
>
> function(AO)
> {A<-AO
> {for(t in 1:4)
> A[t+1]<-matrix(h0m[t,1:4],nrow=2,ncol=2,byrow=FALSE)%*%matrix(h0m[t+1,1:4],nrow=2,ncol=2,byrow=FALSE) 
>
> }
> }
>
> I get this type of error:
>
> 1: number of items to replace is not a multiple of replacement length
> 2: number of items to replace is not a multiple of replacement length
> 3: number of items to replace is not a multiple of replacement length
> 4: number of items to replace is not a multiple of replacement length
>
> I believe I understand the error here in that I am attempting to add 
> more elements to an object that is only a 2x2 matrix. So is there some 
> way to just update the value of the matrix recursively? If you are 
> asking yourself why I am bothering to do this if I mean to multiple 
> only 4 matrices, I actually intend to multiple up to all 2541 of them, 
> which seems no small task. Ultimately, I would like to define the time 
> period over which I am multiplyig these essentially time-dependent 
> matrices and obtain the eigenvalues of the final matrix product.
>
> If anyone can offer any assistance, I would greatly appreciate it. 
> Thanks.
>
> Ivan Kautter
>
> _________________________________________________________________
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
>
>