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

[Ravi Varadhan]

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.