CODEWITHSUNDEEP
rif-statementmatrixlinear-algebra
user5340702
user5340702

Reputation:

Raising a Power on Matrices

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!

Upvotes: 2

Views: 882

Answers (3)

Tim Biegeleisen
Tim Biegeleisen

Reputation: 521194

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

Upvotes: 1

akrun
akrun

Reputation: 887108

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

Upvotes: 3

sunilyadav0201
sunilyadav0201

Reputation: 45

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

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

Upvotes: 0

Related Questions