From: Paul Gilbert <pgilbert_at_bank-banque-canada.ca>

Date: Tue, 08 May 2007 09:13:55 -0400

*> [2,] 0.7 0.3
*

> [2,] 0.42 0.58

> [2,] 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:
*

> like to compute A%*%A%*%A, is there an abbreviation similar to A^3?

> http://www.R-project.org/posting-guide.html

*> ========
*

*>
*

*> 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 Tue 08 May 2007 - 13:20:43 GMT

Date: Tue, 08 May 2007 09:13:55 -0400

Ravi Varadhan wrote:

> Paul,

*>
**> Your solution based on SVD does not work
*

Ooops. I really am getting rusty. The idea is based on eigenvalues which, of course, are not always the same as singular values.

Paul

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,] 0.3 0.7

> [,1] [,2]

>> P %*% P> [1,] 0.58 0.42

> [,1] [,2]

> [2,] 0.42 0.58

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

> [,1] [,2]

> [2,] 0.58 0.42

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

> [,1] [,2]

>> 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. Received on Tue 08 May 2007 - 13:20:43 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 Wed 09 May 2007 - 14:31:21 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.
*