R help - May 2007 - A function for raising a matrix to a power?

Hi,

Is there a function for raising a matrix to a power? For example if you like to
compute A%*%A%*%A, is there an abbreviation similar to A^3?

Atte Tenkanen
> A=rbind(c(1,1),c(-1,-2))
> A
     [,1] [,2]
[1,]    1    1
[2,]   -1   -2
> A^3
     [,1] [,2]
[1,]    1    1
[2,]   -1   -8

But:
> A%*%A%*%A
     [,1] [,2]
[1,]    1    2
[2,]   -2   -5

[Apologies: This will probably break the thread, but at the moment
I cannot receive mail since my remote mail-server is down, and so
am responding to the message as found on the R-help archives;
hence this is not a "Reply"]

From:
Atte Tenkanen attenka at utu.fi
Sun May 6 11:07:07 CEST 2007
> Is there a function for raising a matrix to a power? For example if you
> like to compute A%*%A%*%A, is there an abbreviation similar to A^3?
> 
> Atte Tenkanen
> 
> > A=rbind(c(1,1),c(-1,-2))
> > A
>      [,1] [,2]
> [1,]    1    1
> [2,]   -1   -2
> > A^3
>      [,1] [,2]
> [1,]    1    1
> [2,]   -1   -8
> 
> But:
> 
> > A%*%A%*%A
>      [,1] [,2]
> [1,]    1    2
> [2,]   -2   -5

Of course, "^" alone acts element-wise on the matrix A, so the result
of A^3 is the matrix B where B[i,j] = A[i,j]^3.

However, you can write your own, and call it for instance "%^%":

"%^%"<-function(A,n){
  if(n==1) A else {B<-A; for(i in (2:n)){A<-A%*%B}}; A
  }

Then:
> A
     [,1] [,2]
[1,]    1    1
[2,]   -1   -2
> A%^%1
     [,1] [,2]
[1,]    1    1
[2,]   -1   -2
> A%^%2
     [,1] [,2]
[1,]    0   -1
[2,]    1    3
> A%^%3
     [,1] [,2]
[1,]    1    2
[2,]   -2   -5
> A%^%100
              [,1]          [,2]
[1,] -1.353019e+20 -3.542248e+20
[2,]  3.542248e+20  9.273727e+20


Hoping this helps!
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 06-May-07                                       Time: 11:35:31
------------------------------ XFMail ------------------------------

>  Hi,
>
> Is there a function for raising a matrix to a power? For example if 
> you like to compute A%*%A%*%A, is there an abbreviation similar to A3?
>
> Atte Tenkanen 
Hi Atte,

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)
}

Then, for your example:
> > A=rbind(c(1,1),c(-1,-2))
> > mp(A,3)
>      [,1] [,2]
> [1,]    1    2
> [2,]   -2   -5
Cheers,
Ron

---
Ron E. VanNimwegen
Ph.D. Candidate, Division of Biology (EEB)
Kansas Cooperative Fish & Wildlife Research Unit
205 Leasure Hall
Kansas State University
Manhattan, KS 66506-3501
vanron at ksu.edu
---

Here is a recursive version of the same function:

"%^%" <- function(A, n) if(n == 1) A else A %*% (A %^% (n-1))
> a <- matrix(1:4, 2)
> a %^% 1
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> a %^% 2
     [,1] [,2]
[1,]    7   15
[2,]   10   22
> a %^% 3
     [,1] [,2]
[1,]   37   81
[2,]   54  118


-Christos
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Ron E. 
> VanNimwegen
> Sent: Sunday, May 06, 2007 10:48 AM
> To: r-help at stat.math.ethz.ch; attenka at utu.fi
> Subject: Re: [R] A function for raising a matrix to a power?
> 
> >  Hi,
> >
> > Is there a function for raising a matrix to a power? For example if 
> > you like to compute A%*%A%*%A, is there an abbreviation 
> similar to A3?
> >
> > Atte Tenkanen
> Hi Atte,
> 
> 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)
> }
> 
> Then, for your example:
> > > A=rbind(c(1,1),c(-1,-2))
> > > mp(A,3)
> >      [,1] [,2]
> > [1,]    1    2
> > [2,]   -2   -5
> Cheers,
> Ron
> 
> ---
> Ron E. VanNimwegen
> Ph.D. Candidate, Division of Biology (EEB)
> Kansas Cooperative Fish & Wildlife Research Unit
> 205 Leasure Hall
> Kansas State University
> Manhattan, KS 66506-3501
> vanron at ksu.edu
> ---
> 
> ______________________________________________
> 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.
> 
>

See mtx.exp  in the Malmig package  which is more
efficient than a simple recurvsive routine or alternatively,
%^%		in Lindsey's rmutil package

HTH

ken
> Hi,
>
> Is there a function for raising a matrix to a power? For example if  
> you like to compute A%*%A%*%A, is there an abbreviation similar to  
> A^3?
>
> Atte Tenkanen
>
> > A=rbind(c(1,1),c(-1,-2))
> > A
>      [,1] [,2]
> [1,]    1    1
> [2,]   -1   -2
> > A^3
>      [,1] [,2]
> [1,]    1    1
> [2,]   -1   -8
>
> But:
>
> > A%*%A%*%A
>      [,1] [,2]
> [1,]    1    2
> [2,]   -2   -5
--
Ken Knoblauch
Inserm U846
Institut Cellule Souche et Cerveau
D?partement Neurosciences Int?gratives
18 avenue du Doyen L?pine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.pizzerialesgemeaux.com/u846/

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

At 15:48 06/05/2007, Ron E. VanNimwegen wrote:
>>  Hi,
>>
>>Is there a function for raising a matrix to a power? For example if 
>>you like to compute A%*%A%*%A, is there an abbreviation similar to A3?
>>
>>Atte Tenkanen
I may be revealing my ignorance here, but is MatrixExp in the msm 
package (available from CRAN) not relevant here?

[snip other suggestion]



Michael Dewey
http://www.aghmed.fsnet.co.uk

You might try this, from 9/22/2006 with subject line Exponentiate a matrix:

I am getting a bit rusty on some of these things, but I seem to recall 
that there is a numerical advantage (speed and/or accuracy?) to 
diagonalizing:

 > expM <- function(X,e) { v <- La.svd(X); v$u %*% diag(v$d^e) %*% v$vt
}

 > P <- matrix(c(.3,.7, .7, .3), ncol=2)
 > P %*% P %*% P
      [,1]  [,2]
[1,] 0.468 0.532
[2,] 0.532 0.468
 > expM(P,3)
      [,1]  [,2]
[1,] 0.468 0.532
[2,] 0.532 0.468

I think this also works for non-integer, negative, large, and complex 
exponents:

 > expM(P, 1.5)
          [,1]      [,2]
[1,] 0.3735089 0.6264911
[2,] 0.6264911 0.3735089
 > expM(P, 1000)
     [,1] [,2]
[1,]  0.5  0.5
[2,]  0.5  0.5
 > expM(P, -3)
        [,1]    [,2]
[1,] -7.3125  8.3125
[2,]  8.3125 -7.3125
 > expM(P, 3+.5i)
                  [,1]              [,2]
[1,] 0.4713+0.0141531i 0.5287-0.0141531i
[2,] 0.5287-0.0141531i 0.4713+0.0141531i
 >

Paul Gilbert

Doran, Harold wrote:

 > Suppose I have a square matrix P
 >
 > P <- matrix(c(.3,.7, .7, .3), ncol=2)
 >
 > I know that
 >
 >
 >> P * P
 >
 > Returns the element by element product, whereas
 >
 >
 >
 >> P%*%P
 >>
 >
 > Returns the matrix product.
 >
 > Now, P2 also returns the element by element product. But, is there a
 > slick way to write
 >
 > P %*% P %*% P
 >
 > Obviously, P3 does not return the result I expect.
 >
 > Thanks,
 > Harold
 >


Atte Tenkanen wrote:
> Hi,
> 
> Is there a function for raising a matrix to a power? For example if you
like to compute A%*%A%*%A, is there an abbreviation similar to A^3?
> 
> Atte Tenkanen
> 
>> A=rbind(c(1,1),c(-1,-2))
>> A
>      [,1] [,2]
> [1,]    1    1
> [2,]   -1   -2
>> A^3
>      [,1] [,2]
> [1,]    1    1
> [2,]   -1   -8
> 
> But:
> 
>> A%*%A%*%A
>      [,1] [,2]
> [1,]    1    2
> [2,]   -2   -5
> 
> ______________________________________________
> 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.
===================================================================================
La version fran?aise suit le texte anglais.

------------------------------------------------------------------------------------

This email may contain privileged and/or confidential inform...{{dropped}}