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

From: Paul Gilbert <pgilbert_at_bank-banque-canada.ca>
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] [,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.

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.