Raising a Power on Matrices

Question

I'm trying to create a function that will give me the value of a matrix once it has been raised to a power. This is what I've done so far:

A <- matrix(c(1,2,3,4,0,1,2,3,0,0,1,2,0,0,0,1),nrow=4,ncol=4)

power <- function(A,n){
+ if(n == 0){
+ return(diag(4))
+ }else{
+ return(A%*%A^(n-1))
+ }
+ }

OUTCOME:

> power(A,4)
[,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]   10    1    0    0
[3,]   46   10    1    0 
[4,]  146   46   10    1

This is giving a different value from what my calculator gets and I'm trying to figure what I'm doing wrong. Any help is appreciated!

@akrun Just use `(A%*%A)%*%(A%*%A)` in your R console...it really doesn't appear to be correct. — Tim Biegeleisen, Sep 16 '15 at 05:32
`^` works element-wise. See `A = matrix(1:4, nrow = 2); A^2; A^3`. — Gregor Thomas, Sep 16 '15 at 05:37

akrun · Answer 1 · 2015-09-16T06:07:45.233

We could use %^% from library(expm)

library(expm)
A%*%(A%^%3)

Using this in a function

power <- function(A,n){
   if(n == 0){
     return(diag(4))
   }else{
 return(A%*%(A%^%(n-1)))
 }
}

power(A,4)
#     [,1] [,2] [,3] [,4]
#[1,]    1    0    0    0
#[2,]    8    1    0    0
#[3,]   36    8    1    0
#[4,]  120   36    8    1

According to the description in ?matpow

Compute the k-th power of a matrix. Whereas ‘x^k’ computes element wise powers, ‘x %^% k’ corresponds to k - 1 matrix multiplications, ‘x %% x %% ... %*% x’.

Or a base R option is Reduce with %*% (but this would be slow compared to %^%.

Reduce(`%*%`,replicate(4, A, simplify=FALSE))

In a function,

power1 <- function(A,n){
   if(n == 0){
    return(diag(4))
   }else{
    Reduce(`%*%`,replicate(n, A, simplify=FALSE))
   }
}

power1(A,4)
#     [,1] [,2] [,3] [,4]
#[1,]    1    0    0    0
#[2,]    8    1    0    0
#[3,]   36    8    1    0
#[4,]  120   36    8    1

Tim Biegeleisen · Answer 2 · 2015-09-16T06:29:54.007

You have a problem with the way you are computing your matrix product. I use a while loop inside your power() function instead. It simply multiples the input matrix against itself n times and then returns the result. Here is a base R solution which is a continuation of the direction in which you were already going.

A <- matrix(c(1,2,3,4,0,1,2,3,0,0,1,2,0,0,0,1),nrow=4,ncol=4)

power <- function(A,n){
    B <- diag(nrow(A))
    if (n == 0) {
    return(diag(nrow(A)))
    } else {
        while (n > 0) {
            B <- A%*%B
            n <- n - 1
        }
        return(B)
    }
}

> power(A, 4)
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    8    1    0    0
[3,]   36    8    1    0
[4,]  120   36    8    1

sunilyadav0201 · Answer 3 · 2015-09-16T05:49:11.777

0

I assume you want to make the muliplication of the matrix.You have to first make the multiply the matrix and then try to multiply them using the same powers as you want ,so you can do two things

Write code to multiply the matrix .
Loop the code to multiply.

edited Sep 16 '15 at 05:49

answered Sep 16 '15 at 05:33

sunilyadav0201

45
5

Raising a Power on Matrices

3 Answers3