Just a correction to my previous posting - O(N^2.25) algorithm for multiplying two general NxN matrices was too optimistic! There exists an O(N^a) algorithm with a < 2.4, but the constant multiplying N^a is so big that for N around 1000 it seems that one will be unable to end up with significantly less than N^3 operations.
--- Moshe Olshansky <[EMAIL PROTECTED]> wrote: > 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 <[EMAIL PROTECTED]> > 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 ([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. > > > > ______________________________________________ > 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. > ______________________________________________ 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.