[R] Matrix exponential

[Erin Hodgess]

How about this:
	
> expm
function(x,pow=2) 
{
xd <- diag((eigen(x)$values^pow))
xm <- eigen(x)$vector %*% xd %*% t(eigen(x)$vector)
return(xm)
}
> xa
     [,1] [,2]
[1,]    2    1
[2,]    1    3
> expm(xa,pow=3)
     [,1] [,2]
[1,]   15   20
[2,]   20   35
> xa %*% xa %*% xa
     [,1] [,2]
[1,]   15   20
[2,]   20   35
> 


From: Federico Calboli <f.calboli at ucl.ac.uk>
To: r-help <r-help at stat.math.ethz.ch>
Subject: [R] matrix exponential: M^0~~

Dear All, 

I would like to ask why the zeroeth power of a matrix gives me a matrix
of ones rather than the identity matrix:

> D<-rbind(c(0,0,0),c(0,0,0),c(0,0,0))
> D<-as.matrix(D)
> D
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

> D^0
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1

I would have expected the identity matrix here.

I find the same result with every other square matrix I used.
BTW, I am using R 1.8.1 on Linux Mandrake 9.1

Cheers,

Federico Calboli
-- 



=================================

Federico C. F. Calboli

PLEASE NOTE NEW ADDRESS

Dipartimento di Biologia
Via Selmi 3
40126 Bologna
Italy

tel (+39) 051 209 4187
fax (+39) 051 251 208

f.calboli at ucl.ac.uk