Hello best helpers, I am a new user and I have been struggling for hours. So finally I decide to ask you: If I have a matrix P, and P.2 = P%*%P, and P.3=P.2%*%P is there a function to calculate the power of a matrix?? if not how can i do: for (i in 1:10) {P.i=P^i} after this I need to sum them up and my problem is to combine P and i to P.i can anyone help me please??? Thanks and have a nice day, Pascal. _________________________________________________________________ [[elided Hotmail spam]] [[alternative HTML version deleted]]
> I am a new user and I have been struggling for hours. So finally I > decide to ask you: > If I have a matrix P, and P.2 = P%*%P, and P.3=P.2%*%P > is there a function to calculate the power of a matrix?? if not how cani do:> for (i in 1:10) {P.i=P^i} > after this I need to sum them up and my problem is to combine P and > i to P.i can anyone help me please???Take a look at FAQ on R 7.21 ("How can I turn a string into a variable?"). This code will do what you ask (though the first matrix is named P.1, rather than P for convenience): paste0 <- function(x) paste("P.", x, sep="") P.1 <- matrix(rnorm(25), nrow=5) varnames <- paste0(1:10) for(i in 2:10) assign(paste0(i), get(paste0(i-1)) %*% P.1) Regards, Richie. Mathematical Sciences Unit HSL ------------------------------------------------------------------------ ATTENTION: This message contains privileged and confidential inform...{{dropped:20}}
pascal vrolijk <pascalvrolijk at hotmail.com> wrote in news:BAY120-W47837F2A708683682EB44AACD60 at phx.gbl:> Hello best helpers, > > I am a new user and I have been struggling for hours. So finally I > decide to ask you: If I have a matrix P, and P.2 = P%*%P, and > P.3=P.2%*%P is there a function to calculate the power of a matrix?? > if not how can i do: for (i in 1:10) {P.i=P^i}Three methods show up in a search "r-help power of a matrix": 1) Use a pre-built function: library(msm) ?MatrixExp 2) <https://stat.ethz.ch/pipermail/r-help/2007-May/131330.html> 3) Bill Venables offered this about a week ago in this list: -------------- This is probably as good a way as any way for this kind of problem. First define a binary operator:> "%^%" <- function(x, n)with(eigen(x), vectors %*% (values^n * t(vectors))) Your toy example then becomes> m <- matrix(c(1, 0.4, 0.4, 1), nrow = 2) > m[,1] [,2] [1,] 1.0 0.4 [2,] 0.4 1.0> m %^% (-0.5)[,1] [,2] [1,] 1.0680744 -0.2229201 [2,] -0.2229201 1.0680744 ------------- -- David Winsemius> after this I need to sum them up and my problem is to combine P and > i to P.i can anyone help me please??? > > Thanks and have a nice day, > > Pascal.
> 3) Bill Venables offered this about a week ago in this list: > -------------- > This is probably as good a way as any way for this kind of problem. > First define a binary operator: > > > "%^%" <- function(x, n) > with(eigen(x), vectors %*% (values^n * t(vectors))) >This example only works for _diagonalizable_ matrices. It crashes, for example, in cases like: m <- rbind(c(1,1,0), c(0,1,1), c(0,0,1)) m %^% 2 m %*% m Alberto Monteiro
Hi Pascal, I think the function could be better but try this: # Function: M is your matrix and n MUST be an integer>0 mat.pow<-function(M,n) { result<-M if(n>1){ for ( iter in 2:n) result<-M%*%result result } else {result} result } # The matrix m <- rbind(c(1,1,0), c(0,1,1), c(0,0,1)) # Goal m^2 = m x m goal=m%*%m # matpow res=mat.pow(m,2) # Check point all.equal(goal,res) See RSiteSearch('nth step transition matrices') HTH, Jorge On Tue, May 6, 2008 at 1:37 AM, pascal vrolijk <pascalvrolijk@hotmail.com> wrote:> Hello best helpers, > > I am a new user and I have been struggling for hours. So finally I decide > to ask you: > If I have a matrix P, and P.2 = P%*%P, and P.3=P.2%*%P > is there a function to calculate the power of a matrix?? if not how can i > do: > for (i in 1:10) {P.i=P^i} > after this I need to sum them up and my problem is to combine P and i to > P.i can anyone help me please??? > > Thanks and have a nice day, > > Pascal. > > > _________________________________________________________________ > [[elided Hotmail spam]] > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@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.[[alternative HTML version deleted]]
Alberto Monteiro
2008-May-07 11:49 UTC
[R] How do I write a power of matrices?? [was: How do I write a sum of matrixes??]
Jorge Ivan Velez wrote:> > I think the function could be better but try this: > > # Function: M is your matrix and n MUST be an integer>0 > mat.pow<-function(M,n) { > result<-M > if(n>1){ > for ( iter in 2:n) result<-M%*%result > result > } > else {result} > result > } >There are much more efficient ways to compute a power, just check: http://en.wikipedia.org/wiki/Exponentiation_by_squaring Or, translating the pseudo-code to R: mat.pow <- function(M, n) { result <- diag(1, ncol(M)) while (n > 0) { if (n %% 2 == 1) { result <- result %*% M n <- n - 1 } M <- M %*% M n <- n / 2 } return(result) } # The matrix m <- rbind(c(1,1,0), c(0,1,1), c(0,0,1)) # Goal m^6 = m x m x m x m x m x m goal= m %*% m %*% m %*% m %*% m %*% m # matpow res=mat.pow(m,6) # Check point all.equal(goal,res) This algorithm would be fast, unless n is a _very_ big number. Alberto Monteiro