[R] Efficient way to compute power of a sparse matrix

[Stéphane Dray]

Dear all,
I would like to compute power of a square non symmetric matrix. This is 
a part of a simulation study. Matrices are quite large (e.g., 900 by 
900), and contains many 0 (more than 99 %). I have try the function 
mtx.exp of the Biodem package:

library(Biodem)
m <- matrix(0, 900, 900)
i <- sample(1:900, 3000, replace = T)
j <- sample(1:900, 3000, replace = T)

for(x in 1:3000)
    m[i[x],j[x]] <- runif(1)

system.time(mtx.exp(m, 20))

This returns (sorry, in french):

utilisateur     système      écoulé
      3.646       0.018       3.665

Then, I have try to use the Matrix package and construct my own Mtx.exp 
function:

library(Matrix)
Mtx.exp <- function (X, n)
{
    if (n != round(n)) {
        n <- round(n)
        warning("rounding exponent `n' to", n)
    }
    phi <- Matrix(diag(nrow(X)))
    pot <- X
    while (n > 0) {
        if (n%%2)
            phi <- phi %*% pot
        n <- n%/%2
        pot <- pot %*% pot
    }
    return(phi)
}
M <- Matrix(m)
system.time(Mtx.exp(M,20))

This approach is slower than the classical one:

utilisateur     système      écoulé
      9.265       0.053       9.348

I am quite sure that the problem comes the diagonal matrix used in 
Mtx.exp. it seems that different classes are used in the function which 
induces this lack of speed. I am not sure of which classes of Matrix 
must be used to improve the script. I wonder if someone could help me to 
obtain an efficient (i.e. fastest) procedure to compute the power of a 
matrix.

Please note that in this example, M^n could be a matrix of 0, but it is 
not the case with my real data.

Thanks in advances,
Cheers.

-- 
Stéphane DRAY (dray at biomserv.univ-lyon1.fr )
Laboratoire BBE-CNRS-UMR-5558, Univ. C. Bernard - Lyon I
43, Bd du 11 Novembre 1918, 69622 Villeurbanne Cedex, France
Tel: 33 4 72 43 27 57       Fax: 33 4 72 43 13 88
http://biomserv.univ-lyon1.fr/~dray/