R help - Nov 2007 - Efficient way to compute power of a sparse matrix

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/

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.
>