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/

Moshe Olshansky

2007-Nov-20 04:53 UTC

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

Hello St?phane, Since 20 = 4 + 16 you need 5 matrix multiplications to compute X^20 (2 for X^4, 2 more for X^16 and one more for X^20). If your matrix is NxN, one naive matrix multiplication requires about N^3 operations. In your case N is 900. If it were 1000, 1000^3 is one billion, so 5 matrix multiplications should take several seconds. To remedy the problem, one possibility is to use a faster matrix multiplication which (if I remember correctly - not sure) requires O(N^2.25) operations. Even better possibility is to exploit sparsity. I do not know what operations on sparse matrices the Matrix package supports. If it does not contain what you need you can do this yourself. It takes O(N^2) operations to find all non-zero elements in the matrix. Now if you want to multiply A by B and for each row of A and each column of B you have the indexes of non-zero elements, it will take only O(1) operations to decide whether element (i,j) of the product is non-zero (and to compute it), so matrix multiplication will take just O(N^2) operations (I believe that if N(A) is the number of non-zero elements in A and N(B) is that number for B, you will need O(N*(N(a)+N(B)) operations). So hoping that the powers of you matrix X do not fill in too fast you can compute the power much faster. I believe that there exist packages which can do this for you. Regards, Moshe. --- St?phane Dray <dray at biomserv.univ-lyon1.fr> wrote:> 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/ > > ______________________________________________ > R-help at 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. >