Matrix exponentiation in R
Why is there no matrix exponentiation in R?
I need it all the time, for dealing with continuous time Markov chains, but it just isn't there!
Well, it is probably found in a lot of packages (maybe even some I have installed), but I feel it should be there as a built in function.
I usually end up implementing it myself using eigen value decomposition
mexp <- function(Q,t) {
## calculates matrix exponential using a simple eigen-value decomposition
x <- eigen(Q)
U <- x$vectors
L <- diag(exp(t*x$values))
return (U %*% L %*% solve(U))
}
but that isn't that numerically stable, so it is a sub-optimal solution.
I guess it isn't trivial to get right, but damn I miss it...
--
29-48=-19