From: Alberto Vieira Ferreira Monteiro <albmont_at_centroin.com.br>

Date: Sun, 06 May 2007 16:19:49 +0000

result <- diag(1, ncol(mat))

while (n > 0) {

if (n %% 2 != 0) {

}

mat <- mat %*% mat

n <- n / 2

}

return(result)

}

R-help_at_stat.math.ethz.ch 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 Sun 06 May 2007 - 16:29:21 GMT

Date: Sun, 06 May 2007 16:19:49 +0000

Ron E. VanNimwegen wrote:

*>
*

> I was looking for a similar operator, because R uses scalar products

*> when raising a matrix to a power with "^". There might be something
**> more elegant, but this little loop function will do what you need for a
**> matrix "mat" raised to a power "pow":
**>
**> mp <- function(mat,pow){
**> ans <- mat
**> for ( i in 1:(pow-1)){
**> ans <- mat%*%ans
**> }
**> return(ans)
**> }
**>
*

This function is extremely inefficient for high powers of pow
[an unhappy choice of variables, because pow is the power
function in many languages]

A better method would keep track of the powers of two, and would optimize according to it.

For example, in order to compute mat^42, we just have to compute mat^2, mat^4, mat^8, mat^16, mat^32, and then mat^40 and mat^42, with "only" 7 multiplications, instead of 41.

See the technique in...

http://en.wikipedia.org/wiki/Exponentiation_by_squaring

matrix.power <- function(mat, n)

{

# test if mat is a square matrix # treat n < 0 and n = 0 -- this is left as an exercise # trap non-integer n and return an errorif (n == 1) return(mat)

result <- diag(1, ncol(mat))

while (n > 0) {

if (n %% 2 != 0) {

result <- result %*% mat n <- n - 1

}

mat <- mat %*% mat

n <- n / 2

}

return(result)

}

# Test case:

mat <- matrix(c(1,1,1,0), 2, 2)

mat42 <- matrix.power(mat, 42)

det(mat) # -1

det(mat42) # not 1, because of truncation errors, but close enough :-)
# I hate floating point arithmetics :-/

# Another test case:

mat <- matrix(c(1,1,0,1), 2, 2)

mat314159 <- matrix.power(mat, 314159)

# mat314159 is matrix(c(1,314159,0,1), 2, 2)

Alberto Monteiro

R-help_at_stat.math.ethz.ch 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 Sun 06 May 2007 - 16:29:21 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 Sun 06 May 2007 - 20:33:20 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.
*