Re: [R] A function for raising a matrix to a power?

From: Ravi Varadhan <rvaradhan_at_jhmi.edu>
Date: Mon, 07 May 2007 18:11:28 -0400

Paul,

Your solution based on SVD does not work even for the matrix in your example (the reason it worked for e=3, was that it is an odd power and since P is a permutation matrix. It will also work for all odd powers, but not for even powers).   

However, a spectral decomposition will work for all symmetric matrices (SVD based exponentiation doesn't work for any matrix).

Here is the function to do exponentiation based on spectral decomposition:

expM.sd <- function(X,e){Xsd <- eigen(X); Xsd$vec %*% diag(Xsd$val^e) %*% t(Xsd$vec)}

> P <- matrix(c(.3,.7, .7, .3), ncol=2)
> P

     [,1] [,2]
[1,] 0.3 0.7
[2,] 0.7 0.3
> P %*% P

     [,1] [,2]
[1,] 0.58 0.42
[2,] 0.42 0.58
> expM(P,2)

     [,1] [,2]

[1,] 0.42 0.58
[2,] 0.58 0.42

> expM.sd(P,2)
[,1] [,2]

[1,] 0.58 0.42
[2,] 0.42 0.58
>

Exponentiation based on spectral decomposition should be generally more accurate (not sure about the speed). A caveat is that it will only work for real, symmetric or complex, hermitian matrices. It will not work for asymmetric matrices. It would work for asymmetric matrix if the eigenvectors are made to be orthonormal.

Ravi.



Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan_at_jhmi.edu

Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html  



-----Original Message-----
From: r-help-bounces_at_stat.math.ethz.ch
[mailto:r-help-bounces_at_stat.math.ethz.ch] On Behalf Of Paul Gilbert Sent: Monday, May 07, 2007 5:16 PM
To: Atte Tenkanen
Cc: r-help_at_stat.math.ethz.ch
Subject: Re: [R] A function for raising a matrix to a power?

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



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.

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 Mon 07 May 2007 - 22:20:40 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 Tue 08 May 2007 - 13:31:30 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.