R help - Mar 2004 - recursive matrix multiplication question

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

Is the following what you want:

A1 <- array(1:20, dim=c(5, 4))
seq.mult <-
function(AO){
A<-AO
N <- dim(A)[1]
for(i in 1:(N-1)){
A[i+1, ] <- (matrix(AO[i,1:4],nrow=2,ncol=2,byrow=FALSE)
%*% matrix(AO[i+1,1:4],nrow=2,ncol=2,byrow=FALSE))
}
A
}
seq.mult(A1)
[,1] [,2] [,3] [,4]
[1,] 1 6 11 16
[2,] 79 124 199 344
[3,] 102 157 242 397
[4,] 129 194 289 454
[5,] 160 235 340 515

You got the error (or warning) by omitting the comma, writing A[t+1], 
where I wrote A[i+1,]. Any matrix is also a vector and can be addressed 
as such. If A is a 5 x 4 matrix, as in this example, then A[t+1] 
interprets it as a 1 x 20 vector.

I hope you won't object if I comment on a few other things I noticed: 
First, I assume your data reside in a 2541 x 4 matrix, not a 4 x 2541 
matrix. The modified version of the function should work for the former 
but not the latter. Second, I don't know where "h0m" is defined; I
assume you meant to write AO there.

Also, "t" is the function to transpose a matrix. While R is capable of
distinguishing matrices from functions in most contexts, there are some 
contexts where it is not obvious what is requested, and it is therefore 
usually not prudent to use the name of a function for some other 
purpose. This is done routinely for function arguments. For example, the 
density function for Student's t distribution, dt, has an argument df 
for degrees of freedom, but df is also the density function for the F 
distribution. However, I try to avoid that, to the point of trying names 
I plan to use to see if they are already defined.

Your subject line is misleading: I saw no recursion in your question. A 
"recursion", according to the third (and seemingly most relevant here)
definition in "www.m-w.com", is "a computer programming technique
involving the use of a procedure, subroutine, function, or algorithm 
that calls itself".

You may already know that the function to get eigenvalues and vectors is 
"eigen".

hope this helps.
spencer graves

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

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