From: Ravi Varadhan <rvaradhan_at_jhmi.edu>

Date: Mon, 07 May 2007 18:11:28 -0400

[1,] 0.58 0.42

[2,] 0.42 0.58

*>
*

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

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,] 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.
*