[R] A function for raising a matrix to a power?
Alberto Vieira Ferreira Monteiro
albmont at centroin.com.br
Sun May 6 18:19:49 CEST 2007
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 error
if (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
More information about the R-help
mailing list