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 ([EMAIL PROTECTED] )
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/
______________________________________________
R-help@r-project.org 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.
