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

[Paul Gilbert]

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