[R] A function for raising a matrix to a power?
Ravi Varadhan
rvaradhan at jhmi.edu
Tue May 8 00:11:28 CEST 2007
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.
More information about the R-help
mailing list