[R] better way to iterate matrix multiplication?

[Martin Maechler]

>>>>> "EV" == Eik Vettorazzi <E.Vettorazzi at uke.uni-hamburg.de>
>>>>>     on Tue, 1 Feb 2011 18:02:42 +0100 writes:

    EV> yes there are. But Christofer doesn't need exp(A), but
    EV> A^n.

    EV> But there is a matpow-function %^% in this package, which is a little
    EV> bit slower, I think:

    EV> library(expm)
    EV> states<-1000
    EV> tm <- matrix(runif(states^2),nrow=states)  # random transition matrix
    EV> for illustration
    EV> tm <- t(apply(tm,1,function (x) x/sum(x)))  # make its rows sum to 1
    EV> p0<-pm <- c(0.5,0.5,rep(0,states-2))  # half of cases start in state 1,
    EV> half in state 2
    EV> n<-10000

    EV> system.time({dd<-eigen(tm,symmetric=F)
    EV> as.real(p0%*% dd$vectors%*% diag(dd$values^n)%*%solve(dd$vectors))})
    EV> User      System     elapsed
    EV> 15.20        0.09       15.57

    EV> system.time(p0%*%(tm%^%n))
    EV> User      System     elapsed
    EV> 38.61        0.00       39.62

Indeed, thank you for doing the experiment.

Note however that the eigen() method is only available in *some* cases for
the matrix power, e.g. only when the matrix is non-singular,
whereas  the `%^%` in package  expm  should work reliably in all
cases.

Also for (much) smaller powers than n=10'000
the cpu time needed is more comparable.

Martin Maechler, ETH Zurich


    EV> Am 01.02.2011 17:16, schrieb Ben Bolker:
    >> Eik Vettorazzi <E.Vettorazzi <at> uke.uni-hamburg.de> writes:
    >> 
    >>> 
    >>> if you have a homogeneous mc (= a constant transition matrix), your
    >>> state at time 10 is given by (chapman-kolmogorov)
    >>> p10=p0 %*% tm^(10)
    >>> so you need a matrix power function.
    >> 
    >> There are matrix exponential functions in the Matrix and expm
    >> packages ... don't know about their speed


    EV> -- 
    EV> Eik Vettorazzi
    EV> Institut für Medizinische Biometrie und Epidemiologie
    EV> Universitätsklinikum Hamburg-Eppendorf

    EV> Martinistr. 52
    EV> 20246 Hamburg

    EV> T ++49/40/7410-58243
    EV> F ++49/40/7410-57790