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

[Moshe Olshansky]

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 <m_olshansky at yahoo.com> 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 <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.
> >
> 
> ______________________________________________
> 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.
>